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

(1.河南理工大學(xué) 資源環(huán)境學(xué)院，河南 焦作 454000；2.河南省地質(zhì)局生態(tài)環(huán)境地質(zhì)服務(wù)中心，河南鄭州 450011；3.煤炭安全生產(chǎn)與清潔高效利用省部共建協(xié)同創(chuàng)新中心，河南 焦作 454000)

地?zé)崮茏鳛橐环N蘊藏在地球內(nèi)部儲量巨大的綠色清潔能源，已成為21世紀(jì)能源發(fā)展中不可忽視的可再生能源之一。中國賦存有豐富的地?zé)豳Y源，估算2 000 m深度范圍內(nèi)的地?zé)崮軆α肯喈?dāng)于2 500億t標(biāo)準(zhǔn)煤(2022年中國煤炭產(chǎn)量約40億t)[1]，其開發(fā)利用對提高清潔能源比重、實現(xiàn)“雙碳”目標(biāo)具有重要作用。通常情況下，地?zé)崤c水伴生賦存于巖石介質(zhì)中且以水為載體而開采出來，受埋藏深、影響因素多、時空多變等影響，如開采不合理則會出現(xiàn)儲層結(jié)構(gòu)損害、溫度下降、水質(zhì)惡化、水量減少、水位下降等一系列問題。因此，開展深部地?zé)崴嚅_采井聯(lián)合優(yōu)化研究，對于維持地?zé)豳Y源可持續(xù)開發(fā)、充分發(fā)揮地?zé)崴б?、降低工程運行成本、避免環(huán)境水文地質(zhì)問題發(fā)生具有重要意義[2]。

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

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

上述關(guān)于工程實踐優(yōu)化模型建立及灰狼算法求解，為提高工程效率、降低運行成本做出了貢獻(xiàn)。但在建立地?zé)崴_發(fā)管理模型時，一般以需水條件下的區(qū)域水位降深最小為目標(biāo)，而較少涉及井間干擾最小的問題，導(dǎo)致地?zé)峋_發(fā)效益低下。另外，灰狼算法也存在對初始種群過度依賴、容易造成局部最優(yōu)、過早收斂等問題[13]。筆者以河南開封市城區(qū)1 200～1 400 m熱儲層為研究對象，以地?zé)崴唤瞪钭钚?、各地?zé)峋g水位降深影響值最小、運行費用最低為目標(biāo)建立優(yōu)化管理模型，應(yīng)用Pareto-灰狼算法進(jìn)行求解，以期為深部地?zé)崴_采井聯(lián)合優(yōu)化提供一種高效可行的方法。

1 多目標(biāo)管理模型建立

1.1 地?zé)崴鲾?shù)學(xué)模型

1.1.1 水文地質(zhì)模型

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

圖1 開封市地質(zhì)構(gòu)造平面圖Fig.1 Planar view of geological structures in Kaifeng City

開封凹陷區(qū)地?zé)嶂饕植荚?00～2 600 m深度范圍，尤以600～1 600 m范圍內(nèi)開采最多，屬新近系明化鎮(zhèn)組以及館陶組熱儲層，研究區(qū)目前開采該熱儲層地?zé)峋?8眼，如圖2所示。按埋深分屬于600～800、800～1 000、1 000～1 200、1 200～1 400和 1 400～1 600 m熱儲層。據(jù)調(diào)查，目前地?zé)峋畬嵭袉为氝\行、按需供水模式，長期開采會導(dǎo)致熱儲層水位顯著下降、水井吊泵、降落漏斗擴大及地面沉降等問題，因此實施多開采井聯(lián)合運行模式優(yōu)化研究勢在必行。

圖2 研究區(qū)地?zé)峋植技? 200～1 400 m熱儲層地?zé)峋恢肍ig.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

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

圖3 研究區(qū)柱狀圖Fig.3 Stratigraphic column of Kaifeng City

本文根據(jù)區(qū)域地質(zhì)調(diào)查和前期勘探結(jié)果，選取明化鎮(zhèn)組下段到館陶組1 200～1 400 m熱儲層為對象，區(qū)域為開封市城區(qū)東至東平路、西至四大街，南到隴海鐵路線、北到連霍高速公路，地?zé)峋植荚诔菂^(qū)112 km2范圍內(nèi)，目前有G1?G8共8眼地?zé)峋?，參?shù)見表1[17]。鑒于開采條件下，地?zé)崴蜷_封市城區(qū)的漏斗中心徑流，故四周邊界定性為第二類(定流量)邊界。

表1 1 200～1 400 m熱儲層地?zé)峋畢?shù)[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 數(shù)學(xué)模型

對應(yīng)于承壓二維非穩(wěn)定流運動的數(shù)學(xué)模型由基本微分方程、初始條件和邊界條件構(gòu)成[18]。其中區(qū)域Ω中的微分方程為：

初始條件為：

邊界條件為：

1.2 管理模型

1.2.1 目標(biāo)函數(shù)

開封市地?zé)崴_發(fā)管理模型的目標(biāo)函數(shù)包括以下3個。

1) 水位降深最小

在滿足地?zé)崴杷織l件下，實現(xiàn)全區(qū)地?zé)崴唤瞪钭钚。浜瘮?shù)為：

2) 井間干擾最小

針對現(xiàn)有地?zé)峋季?，最小井間距函數(shù)為：

3) 運行費用最低

為降低地?zé)峋_采運行成本，應(yīng)保證最低的運行費用，其函數(shù)為：

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

1.2.2 約束條件

1) 水位降深約束

管理期末各剖分節(jié)點的水位降深不得超過地?zé)崴O(shè)計允許降深：

2) 供需平衡約束

整個管理期內(nèi)，地?zé)峋┧繎?yīng)滿足用戶需水量，同時應(yīng)小于地?zé)崴畲笤试S開采量：

2 求解方法

2.1 數(shù)學(xué)模型識別

2.1.1 Galerkin有限單元法

有限單元法(Finite Element Method，F(xiàn)EM)是指在研究區(qū)劃分有限個互不重疊的網(wǎng)格單元，用插值函數(shù)表示節(jié)點水頭，然后運用加權(quán)剩余法來離散求解微分方程的算法[19]。依據(jù)離散方法不同，可分為Galerkin和Rayleigh-Ritz有限單元法[20]。

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

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

如果節(jié)點上存在流量為Qi的地?zé)峋?，則式(11)右側(cè)對應(yīng)的節(jié)點應(yīng)加入：

如果地?zé)岢樗?xw,yw)位于三角形單元e內(nèi)，則式(12)右側(cè)對應(yīng)的三角形3個節(jié)點應(yīng)該分別加入：

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

2.1.2 模型識別與驗證

根據(jù)開封市1 200～1 400 m熱儲層的埋藏分布條件[21]、地質(zhì)構(gòu)造展布特征[22]、富水特性[23]、地?zé)峋季趾退膭討B(tài)監(jiān)測情況，將研究區(qū)剖分為90個三角形單元，共162個節(jié)點，其中內(nèi)節(jié)點124個、二類邊界節(jié)點38個，如圖4所示。

圖4 研究區(qū)剖分及節(jié)點分布Fig.4 Partition and node distribution of Kaifeng City

根據(jù)地?zé)峋_采水量，將時間離散為非取暖期和取暖期。非取暖期自3月16日至11月14日，共245 d；取暖期自11月15日至次年3月15日，共計120 d。

根據(jù)地?zé)崴畯搅魈卣?，開封城區(qū)1 200～1 400 m熱儲層地?zé)崴边?、西邊以及南邊CD段概化為定流量補給邊界，東邊和南邊BC段概化為定流量排泄邊界，如圖4所示。根據(jù)水位等值線圖可以計算相應(yīng)的單寬流量，見表2。

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

數(shù)學(xué)模型識別與驗證的關(guān)鍵是水文地質(zhì)參數(shù)的確定。不同地?zé)峋幍臒醿铀牡刭|(zhì)參數(shù)初值，可根據(jù)穩(wěn)定流和非穩(wěn)定流抽水試驗求得[24]。依據(jù)參數(shù)初值并結(jié)合巖性空間變化將研究區(qū)劃分為3個參數(shù)分區(qū)，不同分區(qū)的參數(shù)變化范圍見表3。

表3 不同分區(qū)水文地質(zhì)參數(shù)變化范圍Table 3 Variation ranges of hydrogeological parameters in different zones

利用2018年4月15日至2021年4月15日地?zé)峋?，反演擬合得到不同分區(qū)的水文地質(zhì)參數(shù)；再利用擬合參數(shù)，正演計算2021年4月15日至2023年4月15日地?zé)峋?。將計算水位與實測水位進(jìn)行對比分析，以檢驗擬合參數(shù)是否符合實際。

2022年4月15日全區(qū)計算水位與實測水位等值線如圖5所示；驗證時段內(nèi)，G1、G5、G6、G7地?zé)峋坏挠嬎阒蹬c實測值對比如圖6所示?？梢钥闯?，計算水位和實測水位總體趨勢保持一致且差值較小，說明擬合的水文地質(zhì)參數(shù)符合實際。

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

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

擬合和檢驗后的各分區(qū)水文地質(zhì)參數(shù)值列于表4，這些參數(shù)可以用來表征開封市城區(qū)1 200～1 400 m熱儲層地?zé)崴倪\移特征。

表4 水文地質(zhì)參數(shù)取值Table 4 Values of hydrogeological parameters

2.2 改進(jìn)的灰狼算法

2.2.1 傳統(tǒng)灰狼算法

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

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

和計算式為：

在灰狼種群中，α、β、δ狼是最接近獵物、最能感知獵物信息的群體，其余的灰狼個體的位置都是根據(jù)這3種狼來確定的[27]?；依歉聜€體位置的數(shù)學(xué)描述如下：

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

2.2.2 收斂因子改進(jìn)

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

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

2.2.3 罰函數(shù)引入

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

2.2.4 Pareto解集引入

采用 Pareto 解集對其全區(qū)地?zé)崴唤瞪睢⒏鞯責(zé)峋g水位降深影響值和運行費用進(jìn)行優(yōu)化，通過外部存儲Archive保存每次迭代后的解集，計算非支配解集以獲得最優(yōu)解[31]。

1) 外部種群準(zhǔn)入

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

2) 決策灰狼選擇

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

2.3 模型求解流程及驗證

2.3.1 求解流程

利用改進(jìn)的灰狼算法求解優(yōu)化管理模型的主要流程如圖8所示。

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

2.3.2 有效性驗證分析

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

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

3 結(jié)果與分析

3.1 參數(shù)設(shè)置

本文應(yīng)用Matlab語言編寫多目標(biāo)灰狼優(yōu)化算法程序，基本參數(shù)見表5。

表5 算法基本參數(shù)設(shè)置Table 5 Settings of the general parameters in the algorithm

3.2 求解結(jié)果

以2022年3月15日水位為初始值，在側(cè)向補給量及水文地質(zhì)參數(shù)相同條件下，優(yōu)化2022年3月16日至2022年11月14日非取暖期、2022年11月15日至2023年3月15日取暖期的地?zé)崴_采量，所求8眼地?zé)峋?2組Pareto最優(yōu)解(序號P1?P12)見表6。

表6 多目標(biāo)優(yōu)化Pareto最優(yōu)解Table 6 Pareto optimal solutions for multi-objective optimization

為直觀地給出最終優(yōu)化方案，參考前人成果[35]，將目標(biāo)函數(shù)h1、h2、C3權(quán)重分別設(shè)置為0.4、0.3、0.3。針對在 Pareto前沿上選取的12個最優(yōu)解，采用TOPSIS法[36]進(jìn)行排序，其結(jié)果見表7。

表7 地?zé)崴嗄繕?biāo)優(yōu)化TOPSIS排序結(jié)果Table 7 TOPSIS ranking results for multi-objective optimization of geothermal water

依據(jù)表7排序結(jié)果，選擇方案8作為地?zé)崴嗄繕?biāo)優(yōu)化配置結(jié)果，對應(yīng)的各節(jié)點水位變化如圖10所示。

圖10 3個分區(qū)優(yōu)化前后水位對比Fig.10 Comparison of water levels in three zones before and after optimization

3.3 結(jié)果分析

3.3.1 地?zé)崴_采量

優(yōu)化管理前后的地?zé)崴_采量對比見表8，目標(biāo)函數(shù)值對比見表9。研究區(qū)內(nèi)8眼地?zé)峋_采時，運行費用由188 055.01元減少至134 083.50元，下降28.7%，成本顯著降低。

表8 優(yōu)化前后地?zé)峋_采量對比Table 8 Comparison of geothermal well production before and after optimization

表9 優(yōu)化前后目標(biāo)函數(shù)值對比Table 9 Comparison of objective function values before and after optimization

3.3.2 水位降深

由圖10結(jié)果可以看出，優(yōu)化管理前后相比，水位降深顯著減少，優(yōu)化管理對抑制水位降低起到較好作用。針對3個分區(qū)，分別計算節(jié)點水位下降速率的均值和均方差，箱線圖如圖11所示。顯然，優(yōu)化后各分區(qū)地?zé)崴陆邓俾识嫉陀趦?yōu)化前，分區(qū)Ⅰ由優(yōu)化前的7.86 m變至6.85 m，變化率為?12.85%；分區(qū)Ⅱ由優(yōu)化前的9.31 m變至6.69 m，變化率為?28.14%；其中分區(qū)Ⅲ變化最大，由優(yōu)化前的11.79 m變至6.63 m，變化率為?43.77%。表明本文形成的多開采井聯(lián)合運行優(yōu)化方法在降低地?zé)峋幌陆捣矫鎯?yōu)勢明顯。

圖11 分區(qū)Ⅰ優(yōu)化前后水位對比Fig.11 Comparison of water levels in zone I before and after optimization

3.3.3 降落漏斗面積

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

圖12 分區(qū)Ⅱ優(yōu)化前后水位對比Fig.12 Comparison of water levels in zone II before and after optimization

4 結(jié)論

a.構(gòu)建了河南開封市城區(qū)埋深1 200～1 400 m熱儲層水文地質(zhì)模型及數(shù)學(xué)模型，基于Galerkin有限單元法、應(yīng)用近5 a地?zé)崴粍討B(tài)資料進(jìn)行模型識別與驗證，為表征深埋地?zé)崴\移特征奠定了基礎(chǔ)。

b.選取多目標(biāo)函數(shù)，多約束條件，建立了地?zé)崴嚅_采井運行優(yōu)化模型。以收斂因子、死亡罰函數(shù)、Pareto解集為切入點，改進(jìn)了基本灰狼算法，標(biāo)準(zhǔn)測試函數(shù)的測試效果良好，形成一種地?zé)崴_采多目標(biāo)函數(shù)管理模型優(yōu)化配置新途徑。

c.優(yōu)化前后各節(jié)點水位降深之和由640.7 m降至190.3 m，井間水位降深影響值之和由1 465.25 m降至1 306.59 m，運行費用由188 055.01元降至134 083.50元。優(yōu)化管理期末，分區(qū)Ⅰ、Ⅱ、Ⅲ水位降深較現(xiàn)狀開采分別下降12.85%、28.14%、43.77%。?80和?90 m水位等值線包圍的降落漏斗面積分別下降15.5%和28.7%。研究成果對于保證深部地?zé)崴沙掷m(xù)開發(fā)利用具有重要意義。

符號注釋：

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