深部地熱水多開采井聯合運行模式優化

平 宇，李 放，季文清，袁筱瑩，齊玉峰，鄧曉穎，王心義,3,*

(1.河南理工大學 資源環境學院，河南 焦作 454000；2.河南省地質局生態環境地質服務中心，河南鄭州 450011；3.煤炭安全生產與清潔高效利用省部共建協同創新中心，河南 焦作 454000)

地熱能作為一種蘊藏在地球內部儲量巨大的綠色清潔能源，已成為21世紀能源發展中不可忽視的可再生能源之一。中國賦存有豐富的地熱資源，估算2 000 m深度范圍內的地熱能儲量相當于2 500億t標準煤(2022年中國煤炭產量約40億t)[1]，其開發利用對提高清潔能源比重、實現“雙碳”目標具有重要作用。通常情況下，地熱與水伴生賦存于巖石介質中且以水為載體而開采出來，受埋藏深、影響因素多、時空多變等影響，如開采不合理則會出現儲層結構損害、溫度下降、水質惡化、水量減少、水位下降等一系列問題。因此，開展深部地熱水多開采井聯合優化研究，對于維持地熱資源可持續開發、充分發揮地熱水效益、降低工程運行成本、避免環境水文地質問題發生具有重要意義[2]。

近年來，地熱水資源開發受到社會重視，開展了大量研究，目前國內外學者都是以地熱資源可持續開發利用為研究目的，考慮經濟效益、社會效益以及環境效益等方面為研究目標，開展優化管理工作。N.John等[3]基于地理信息系統(Geographic Information System，GIS)下層次分析法(Analytic Hierarchy Process，AHP)和加權乘積模型(Weighted Aggregated Sum Product Assessment，WASPA)的多準則決策方法，確定了可開采地熱資源量條件下的最佳利用方案，為科學管理地熱資源奠定了基礎。Hou Shengjun等[4]協同地熱能、生物質能和空氣能，依據能量平衡原理，以年運行成本最小為目標，構建了TRNSYS軟件支持下的動態仿真模型，為多能源供熱比的動態優化提供了支撐。國內也有相關研究，付銀環等[5]基于多目標規劃、模糊規劃和區間規劃原理，以經濟效益、社會效益和環境效益最大為目標，構建一種適用于多水源、多子區、多用戶的多目標模糊規劃模型，實現系統運行綜合效益最大。王一杰等[6]以安徽省泗縣為研究區，以經濟效益最大、社會缺水量最小和污染物排放量最小為目標函數，建立了水資源多目標優化配置模型，為水資源合理配置提供了依據。

Mirijalili受狼群捕獵啟發，沿用仿生學原理，于2014年提出了灰狼算法。該優化算法具有群體智能、調節參數少、結構簡單、易實現等特點，學者對其原理進行深入研究，對算法進行優化改進，克服此算法在處理復雜問題上存在的種種問題，并被廣泛應用于工程實踐中[7]。Li Su等[8]基于多目標灰狼優化算法，以經濟效益最大和社會缺水量最小為目標函數，優化了領導狼的選擇機制，為水資源多目標優化配置問題提供了一種新方法。鄧飛等[9]針對目前灰狼優化算法存在的缺陷，引入遺傳算法中的交叉操作和大規模領域搜索的破壞修復，進一步加強了局部搜索能力。吳新忠等[10]基于礦井智能控制系統實驗平臺驗證應急調風方案，應用差分灰狼算法得到了最優的可調分支集和風阻調節范圍。胡軍等[11]建立了支持向量機(SVM)預測模型，引入多策略改進灰狼優化算法確定了模型參數，對尾礦壩地下水位預測具有較好的適用性和可靠性。門飛等[12]針對露天礦煤炭運輸調度問題，應用改進灰狼算法求解了優化調度模型，實現了全局和局部尋優能力的均衡，有效降低了煤炭運輸費用。

上述關于工程實踐優化模型建立及灰狼算法求解，為提高工程效率、降低運行成本做出了貢獻。但在建立地熱水開發管理模型時，一般以需水條件下的區域水位降深最小為目標，而較少涉及井間干擾最小的問題，導致地熱井開發效益低下。另外，灰狼算法也存在對初始種群過度依賴、容易造成局部最優、過早收斂等問題[13]。筆者以河南開封市城區1 200～1 400 m熱儲層為研究對象，以地熱水位降深最小、各地熱井之間水位降深影響值最小、運行費用最低為目標建立優化管理模型，應用Pareto-灰狼算法進行求解，以期為深部地熱水開采井聯合優化提供一種高效可行的方法。

1 多目標管理模型建立

1.1 地熱水流數學模型

1.1.1 水文地質模型

河南開封市城區面積約225 km2。在地貌上位于濟源?開封凹陷區東南隅，先后經歷了印支期、早燕山期、晚燕山期和喜馬拉雅期等構造運動，多期運動的疊加使區內形成了一系列次級斷(凹)陷盆地和斷塊凸起，如濟源?開封凹陷和通許凸起等[14]，如圖1所示。

圖1 開封市地質構造平面圖Fig.1 Planar view of geological structures in Kaifeng City

開封凹陷區地熱主要分布在300～2 600 m深度范圍，尤以600～1 600 m范圍內開采最多，屬新近系明化鎮組以及館陶組熱儲層，研究區目前開采該熱儲層地熱井88眼，如圖2所示。按埋深分屬于600～800、800～1 000、1 000～1 200、1 200～1 400和 1 400～1 600 m熱儲層。據調查，目前地熱井實行單獨運行、按需供水模式，長期開采會導致熱儲層水位顯著下降、水井吊泵、降落漏斗擴大及地面沉降等問題，因此實施多開采井聯合運行模式優化研究勢在必行。

圖2 研究區地熱井分布及1 200～1 400 m熱儲層地熱井位置Fig.2 Distribution of geothermal wells and the locations of geothermal wells for geothermal reservoirs at burial depths ranging from 1 200 to 1 400 m in Kaifeng City

根據以往獲得的鉆探資料[15]表明，各熱儲層產狀平緩、厚度穩定、巖性均勻，頂板、底板泥巖較厚且穩定，具有良好的隔水作用，使得上、下熱儲層間均無明顯的水力聯系[16]，如圖3所示。各熱儲層熱水的補給來源在70 km以上，補給微弱、徑流緩慢。地熱水水位隨開采量的變化而急劇變化，呈明顯的非穩定流特征。因此，各熱儲層單獨概化為無垂向補給的承壓二維非穩定流運動水文地質模型。

圖3 研究區柱狀圖Fig.3 Stratigraphic column of Kaifeng City

本文根據區域地質調查和前期勘探結果，選取明化鎮組下段到館陶組1 200～1 400 m熱儲層為對象，區域為開封市城區東至東平路、西至四大街，南到隴海鐵路線、北到連霍高速公路，地熱井分布在城區112 km2范圍內，目前有G1?G8共8眼地熱井，參數見表1[17]。鑒于開采條件下，地熱水向開封市城區的漏斗中心徑流，故四周邊界定性為第二類(定流量)邊界。

表1 1 200～1 400 m熱儲層地熱井參數[17]Table 1 Parameters of geothermal wells for geothermal reservoirs at burial depths ranging from 1 200 to 1 400 m[17]

1.1.2 數學模型

對應于承壓二維非穩定流運動的數學模型由基本微分方程、初始條件和邊界條件構成[18]。其中區域Ω中的微分方程為：

初始條件為：

邊界條件為：

1.2 管理模型

1.2.1 目標函數

開封市地熱水開發管理模型的目標函數包括以下3個。

1) 水位降深最小

在滿足地熱水需水量條件下，實現全區地熱水位降深最小，其函數為：

2) 井間干擾最小

針對現有地熱井布局，最小井間距函數為：

3) 運行費用最低

為降低地熱井開采運行成本，應保證最低的運行費用，其函數為：

在同一揚程下，當電價及水泵效率不同時，提取1 m3水體所消耗的電費Cd為：

1.2.2 約束條件

1) 水位降深約束

管理期末各剖分節點的水位降深不得超過地熱水設計允許降深：

2) 供需平衡約束

整個管理期內，地熱井供水量應滿足用戶需水量，同時應小于地熱水最大允許開采量：

2 求解方法

2.1 數學模型識別

2.1.1 Galerkin有限單元法

有限單元法(Finite Element Method，FEM)是指在研究區劃分有限個互不重疊的網格單元，用插值函數表示節點水頭，然后運用加權剩余法來離散求解微分方程的算法[19]。依據離散方法不同，可分為Galerkin和Rayleigh-Ritz有限單元法[20]。

對式(1)?式(3)，運用Galerkin法得到的有限元方程為：

由三角元離散式(10)并以矩陣的形式表現為：

如果節點上存在流量為Qi的地熱井，則式(11)右側對應的節點應加入：

如果地熱抽水井(xw,yw)位于三角形單元e內，則式(12)右側對應的三角形3個節點應該分別加入：

a、、c采用以下公式計算：

2.1.2 模型識別與驗證

根據開封市1 200～1 400 m熱儲層的埋藏分布條件[21]、地質構造展布特征[22]、富水特性[23]、地熱井布局和水文動態監測情況，將研究區剖分為90個三角形單元，共162個節點，其中內節點124個、二類邊界節點38個，如圖4所示。

圖4 研究區剖分及節點分布Fig.4 Partition and node distribution of Kaifeng City

根據地熱井開采水量，將時間離散為非取暖期和取暖期。非取暖期自3月16日至11月14日，共245 d；取暖期自11月15日至次年3月15日，共計120 d。

根據地熱水徑流特征，開封城區1 200～1 400 m熱儲層地熱水北邊、西邊以及南邊CD段概化為定流量補給邊界，東邊和南邊BC段概化為定流量排泄邊界，如圖4所示。根據水位等值線圖可以計算相應的單寬流量，見表2。

表2 邊界單寬補給量Table 2 Recharge per unit width at boundaries 單位：m3/(d·m)

數學模型識別與驗證的關鍵是水文地質參數的確定。不同地熱井處的熱儲層水文地質參數初值，可根據穩定流和非穩定流抽水試驗求得[24]。依據參數初值并結合巖性空間變化將研究區劃分為3個參數分區，不同分區的參數變化范圍見表3。

表3 不同分區水文地質參數變化范圍Table 3 Variation ranges of hydrogeological parameters in different zones

利用2018年4月15日至2021年4月15日地熱井水位，反演擬合得到不同分區的水文地質參數；再利用擬合參數，正演計算2021年4月15日至2023年4月15日地熱井水位。將計算水位與實測水位進行對比分析，以檢驗擬合參數是否符合實際。

2022年4月15日全區計算水位與實測水位等值線如圖5所示；驗證時段內，G1、G5、G6、G7地熱井水位的計算值與實測值對比如圖6所示。可以看出，計算水位和實測水位總體趨勢保持一致且差值較小，說明擬合的水文地質參數符合實際。

圖5 2022年4月15日計算與實測水位等值線Fig.5 Calculated and measured water level contours of April 15 in 2022

圖6 G1、G5、G6、G7井的計算與實測水位動態Fig.6 Dynamic curves of calculated and measured water levels in wells G1,G5,G6 and G7

擬合和檢驗后的各分區水文地質參數值列于表4，這些參數可以用來表征開封市城區1 200～1 400 m熱儲層地熱水的運移特征。

表4 水文地質參數取值Table 4 Values of hydrogeological parameters

2.2 改進的灰狼算法

2.2.1 傳統灰狼算法

灰狼算法(Grey Wolf Optimizer，GWO)是Mirjalili模仿灰狼的領導階層和狩獵機制提出的一種新的啟發式算法。灰狼嚴格遵守等級關系。α狼位于最高等級，負責決策等總體事務；β狼位于第二級，協助α狼決策及指揮其他低等級狼；δ狼位于第三級，服從α、β狼的命令，并負責偵查、放哨等事務；ω狼等級最低，負責平衡種族內部關系。在 GWO 算法中，每只灰狼代表種群中的一個候選解。根據灰狼的等級制度，α為最優解，β為次優解，δ為第三優解，其余解均稱為ω[25]。

灰狼群體狩獵時，狼群在α、β和δ的帶領下對獵物進行包圍、獵捕和攻擊。狼群對獵物的包圍可用以下公式描述[26]：

和計算式為：

在灰狼種群中，α、β、δ狼是最接近獵物、最能感知獵物信息的群體，其余的灰狼個體的位置都是根據這3種狼來確定的[27]。灰狼更新個體位置的數學描述如下：

說明：式(21)定義了α、β、δ灰狼與當前灰狼之間的距離(、、)，式(22)定義了ω向α、β、δ移動的步長和距離；式(23)則是ω的最終位置。

2.2.2 收斂因子改進

GWO 算法中，當協同系數向量A>1，種群進行全局搜索；當A<1，種群進行局部搜索。由此可見，參數的設置直接影響了算法搜索的性能。傳統 GWO 算法中參數隨迭代過程由2線性減小至0，不能貼合實際的尋優過程。為此本文引入了一種非線性參數策略，考慮到余弦函數具有周期性，其波形可以在全局范圍內進行搜索，有助于算法跳出局部最優解，并在整個搜索空間中進行探索，提出了一種基于余弦規律變化的收斂因子，可以更好地對全局搜索和局部搜索進行調整，有效保持了種群的多樣性[28]，其修正表達式為：

圖7 收斂因子變化規律Fig.7 Variation laws of convergence factors

2.2.3 罰函數引入

在各種罰函數法中，死亡罰函數法是處理約束最簡單、最高效的一種方法，因為它不需要估算得到的解的非可行程度，只要是非可行解就直接丟棄[30]。在灰狼算法中，每個個體計算出的罰函數即為個體的適應度值，公式如下：

2.2.4 Pareto解集引入

采用 Pareto 解集對其全區地熱水位降深、各地熱井之間水位降深影響值和運行費用進行優化，通過外部存儲Archive保存每次迭代后的解集，計算非支配解集以獲得最優解[31]。

1) 外部種群準入

算法迭代時，判斷新灰狼是否可以加入Archive的條件如下：如果新灰狼被所有原始灰狼支配，則其不能加入；若新灰狼被其中一只及以上而非所有灰狼支配，則新灰狼加入，被支配灰狼踢出；若新灰狼與原始灰狼互相不支配，在未達到上限情況下新灰狼加入[32]。

2) 決策灰狼選擇

采用輪盤賭的方式從 Archive 中選擇決策狼。在Archive 中尋找擁擠度最小的網格，隨機選擇α、β和δ對應的解決方案，沒有優劣之分。若該段數量不夠，則順延至擁擠度次之的網格進行選擇。

2.3 模型求解流程及驗證

2.3.1 求解流程

利用改進的灰狼算法求解優化管理模型的主要流程如圖8所示。

圖8 算法流程圖Fig.8 Flow chart of the improved Grey Wolf Optimizer

2.3.2 有效性驗證分析

為驗證算法的有效性，需用標準測試函數進行驗證。本文選擇Zhang Qingfu等[33]在CEC2009上提出的常用多目標測試函數UF9和CF9[34]。測試函數的決策變量設為30，算法迭代300次。改進灰狼算法得到的Pareto前沿(簡稱計算Pareto前沿)和真實Pareto前沿對比如圖9所示。顯然，測試函數的計算Pareto前沿較好地逼近真實Pareto前沿，算法有效性得到證實。

圖9 測試函數UF9和CF9的Pareto前沿Fig.9 Pareto frontiers of test function UF9 and CF9

3 結果與分析

3.1 參數設置

本文應用Matlab語言編寫多目標灰狼優化算法程序，基本參數見表5。

表5 算法基本參數設置Table 5 Settings of the general parameters in the algorithm

3.2 求解結果

以2022年3月15日水位為初始值，在側向補給量及水文地質參數相同條件下，優化2022年3月16日至2022年11月14日非取暖期、2022年11月15日至2023年3月15日取暖期的地熱水開采量，所求8眼地熱井的12組Pareto最優解(序號P1?P12)見表6。

表6 多目標優化Pareto最優解Table 6 Pareto optimal solutions for multi-objective optimization

為直觀地給出最終優化方案，參考前人成果[35]，將目標函數h1、h2、C3權重分別設置為0.4、0.3、0.3。針對在 Pareto前沿上選取的12個最優解，采用TOPSIS法[36]進行排序，其結果見表7。

表7 地熱水多目標優化TOPSIS排序結果Table 7 TOPSIS ranking results for multi-objective optimization of geothermal water

依據表7排序結果，選擇方案8作為地熱水多目標優化配置結果，對應的各節點水位變化如圖10所示。

圖10 3個分區優化前后水位對比Fig.10 Comparison of water levels in three zones before and after optimization

3.3 結果分析

3.3.1 地熱水開采量

優化管理前后的地熱水開采量對比見表8，目標函數值對比見表9。研究區內8眼地熱井開采時，運行費用由188 055.01元減少至134 083.50元，下降28.7%，成本顯著降低。

表8 優化前后地熱井開采量對比Table 8 Comparison of geothermal well production before and after optimization

表9 優化前后目標函數值對比Table 9 Comparison of objective function values before and after optimization

3.3.2 水位降深

由圖10結果可以看出，優化管理前后相比，水位降深顯著減少，優化管理對抑制水位降低起到較好作用。針對3個分區，分別計算節點水位下降速率的均值和均方差，箱線圖如圖11所示。顯然，優化后各分區地熱水下降速率都低于優化前，分區Ⅰ由優化前的7.86 m變至6.85 m，變化率為?12.85%；分區Ⅱ由優化前的9.31 m變至6.69 m，變化率為?28.14%；其中分區Ⅲ變化最大，由優化前的11.79 m變至6.63 m，變化率為?43.77%。表明本文形成的多開采井聯合運行優化方法在降低地熱井水位下降方面優勢明顯。

圖11 分區Ⅰ優化前后水位對比Fig.11 Comparison of water levels in zone I before and after optimization

3.3.3 降落漏斗面積

圖12為優化前后管理期末(2023年4月15日)地熱水水位等值線。優化前后相比，?80和?90 m等值線包圍的降落漏斗面積分別下降15.5%和28.7%，優化結果有效減緩了降落漏斗的擴散。這對于防止熱儲層疏干、水源枯竭、地面沉降等環境水文地質問題的發生具有重要作用。

圖12 分區Ⅱ優化前后水位對比Fig.12 Comparison of water levels in zone II before and after optimization

4 結論

a.構建了河南開封市城區埋深1 200～1 400 m熱儲層水文地質模型及數學模型，基于Galerkin有限單元法、應用近5 a地熱水水位動態資料進行模型識別與驗證，為表征深埋地熱水運移特征奠定了基礎。

b.選取多目標函數，多約束條件，建立了地熱水多開采井運行優化模型。以收斂因子、死亡罰函數、Pareto解集為切入點，改進了基本灰狼算法，標準測試函數的測試效果良好，形成一種地熱水開采多目標函數管理模型優化配置新途徑。

c.優化前后各節點水位降深之和由640.7 m降至190.3 m，井間水位降深影響值之和由1 465.25 m降至1 306.59 m，運行費用由188 055.01元降至134 083.50元。優化管理期末，分區Ⅰ、Ⅱ、Ⅲ水位降深較現狀開采分別下降12.85%、28.14%、43.77%。?80和?90 m水位等值線包圍的降落漏斗面積分別下降15.5%和28.7%。研究成果對于保證深部地熱水可持續開發利用具有重要意義。

符號注釋：

a、b、c為擬合參數；B為優化管理期年限，a；C3為管理年限內地熱井運行費用，元；C1(i,b)、C2(i,b)分別為第b年非取暖期和取暖期內，開采單位第i眼井地熱水量的運行費用，元/m3；、，、，、分別為α、β、δ的導數向量；和為系數向量；D為導水矩陣；g為迭代次數；h1為管理期末各節點水位降深之和，m；h2為管理年限末各地熱井之間水位降深影響值之和，m；h1(j,b)、h2(j,b)分別為第b年取暖期末和非取暖期末，第j個節點的水位降深值，m(i,j,b)、(i,j,b)分別為第b年取暖期末和非取暖期末，第i眼井開采時對第j眼井的水位降深影響值，i≠j，m；f為現行電價，元/(kw·h)；F為常數項矢量；Fτ為無約束條件下的目標函數值；h(j)、hp(j)分別為管理年限末第j個節點的水位降深和設計允許降深，m；H為熱儲層水位高程，m；Hj為熱儲層第j個節點水位高程，m；H0為水頭初始值，m；He為水位埋深，m；為水頭近似解，m；K為滲透系數，m/d；M為熱儲層厚度，m；n為邊界Γ2的外法線方向；N為管理區內剖分節點總數，個；Nj為地熱井數；?為懲罰項，一般取很大的正數；P為貯水矩陣；q為第二類邊界(Γ2)上的單位寬度側向補給量，補給取“+”，排泄取“?”，m3/(d?m)；Q為注水井流量，m3；Q1(i,b)、Q2(i,b)分別為第b年非取暖期和取暖期內，第i眼井開采的地熱水量，m3；Qk(i)為第i眼井的最大允許開采量，m3；Qp(i)為第i眼井用戶需水量，m3；r為貼現率(折扣率)；s為斷面；S為熱儲層彈性釋水系數；t為時間，s；T=KM為熱儲層導水系數，m2/d；和為位于[0,1]的隨機向量；W為源匯項，補給取正，排泄(含開采量)取負，m/d；x、y為坐標，m；、、分別為當前灰狼與α、β、δ灰狼進行擾動后所對應的位置；(g)為獵物的位置；(g)為第g代灰狼個體的位置；、、分別為灰狼α、β、δ當前的位置；(g+1)為當前目標位置向量；v為當前迭代次數；V為總迭代次數；為狼群與獵物之間的距離，m；γ為遞減參數，取值0～1，本文取0.7；η為水泵效率，%；ψ為基函數；μ為個體不滿足約束條件的個數；Πτ(x)為個體適應度值，τ=1,2,3；Δe為單元e的面積，m2；?w為井管面積，m2；Γw為井周界；σw為井徑，m；為收斂因子，隨著迭代次數從 2 線性減少到 0。