深部地?zé)崴嚅_(kāi)采井聯(lián)合運(yùn)行模式優(yōu)化

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

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

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

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

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

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

1.1.1 水文地質(zhì)模型

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

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

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

圖2 研究區(qū)地?zé)峋植技? 200～1 400 m熱儲(chǔ)層地?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ǔ)層產(chǎn)狀平緩、厚度穩(wěn)定、巖性均勻，頂板、底板泥巖較厚且穩(wěn)定，具有良好的隔水作用，使得上、下熱儲(chǔ)層間均無(wú)明顯的水力聯(lián)系[16]，如圖3所示。各熱儲(chǔ)層熱水的補(bǔ)給來(lái)源在70 km以上，補(bǔ)給微弱、徑流緩慢。地?zé)崴浑S開(kāi)采量的變化而急劇變化，呈明顯的非穩(wěn)定流特征。因此，各熱儲(chǔ)層單獨(dú)概化為無(wú)垂向補(bǔ)給的承壓二維非穩(wěn)定流運(yùn)動(dòng)水文地質(zhì)模型。

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

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

表1 1 200～1 400 m熱儲(chǔ)層地?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é)模型

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

初始條件為：

邊界條件為：

1.2 管理模型

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

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

1) 水位降深最小

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

2) 井間干擾最小

針對(duì)現(xiàn)有地?zé)峋季郑钚【g距函數(shù)為：

3) 運(yùn)行費(fèi)用最低

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

在同一揚(yáng)程下，當(dāng)電價(jià)及水泵效率不同時(shí)，提取1 m3水體所消耗的電費(fèi)Cd為：

1.2.2 約束條件

1) 水位降深約束

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

2) 供需平衡約束

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

2 求解方法

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

2.1.1 Galerkin有限單元法

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

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

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

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

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

a、、c采用以下公式計(jì)算：

2.1.2 模型識(shí)別與驗(yàn)證

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

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

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

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

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

數(shù)學(xué)模型識(shí)別與驗(yàn)證的關(guān)鍵是水文地質(zhì)參數(shù)的確定。不同地?zé)峋幍臒醿?chǔ)層水文地質(zhì)參數(shù)初值，可根據(jù)穩(wěn)定流和非穩(wěn)定流抽水試驗(yàn)求得[24]。依據(jù)參數(shù)初值并結(jié)合巖性空間變化將研究區(qū)劃分為3個(gè)參數(shù)分區(qū)，不同分區(qū)的參數(shù)變化范圍見(jiàn)表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ù)，正演計(jì)算2021年4月15日至2023年4月15日地?zé)峋?。將?jì)算水位與實(shí)測(cè)水位進(jìn)行對(duì)比分析，以檢驗(yàn)擬合參數(shù)是否符合實(shí)際。

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

圖5 2022年4月15日計(jì)算與實(shí)測(cè)水位等值線Fig.5 Calculated and measured water level contours of April 15 in 2022

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

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

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

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

和計(jì)算式為：

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

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

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

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

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

2.2.3 罰函數(shù)引入

在各種罰函數(shù)法中，死亡罰函數(shù)法是處理約束最簡(jiǎn)單、最高效的一種方法，因?yàn)樗恍枰浪愕玫降慕獾姆强尚谐潭?，只要是非可行解就直接丟棄[30]。在灰狼算法中，每個(gè)個(gè)體計(jì)算出的罰函數(shù)即為個(gè)體的適應(yīng)度值，公式如下：

2.2.4 Pareto解集引入

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

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

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

2) 決策灰狼選擇

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

2.3 模型求解流程及驗(yàn)證

2.3.1 求解流程

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

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

2.3.2 有效性驗(yàn)證分析

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

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

3 結(jié)果與分析

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

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

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

3.2 求解結(jié)果

以2022年3月15日水位為初始值，在側(cè)向補(bǔ)給量及水文地質(zhì)參數(shù)相同條件下，優(yōu)化2022年3月16日至2022年11月14日非取暖期、2022年11月15日至2023年3月15日取暖期的地?zé)崴_(kāi)采量，所求8眼地?zé)峋?2組Pareto最優(yōu)解(序號(hào)P1?P12)見(jiàn)表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。針對(duì)在 Pareto前沿上選取的12個(gè)最優(yōu)解，采用TOPSIS法[36]進(jìn)行排序，其結(jié)果見(jiàn)表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é)果，對(duì)應(yīng)的各節(jié)點(diǎn)水位變化如圖10所示。

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

3.3 結(jié)果分析

3.3.1 地?zé)崴_(kāi)采量

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

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

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

3.3.2 水位降深

由圖10結(jié)果可以看出，優(yōu)化管理前后相比，水位降深顯著減少，優(yōu)化管理對(duì)抑制水位降低起到較好作用。針對(duì)3個(gè)分區(qū)，分別計(jì)算節(jié)點(diǎn)水位下降速率的均值和均方差，箱線圖如圖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%。表明本文形成的多開(kāi)采井聯(lián)合運(yùn)行優(yōu)化方法在降低地?zé)峋幌陆捣矫鎯?yōu)勢(shì)明顯。

圖11 分區(qū)Ⅰ優(yōu)化前后水位對(duì)比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é)果有效減緩了降落漏斗的擴(kuò)散。這對(duì)于防止熱儲(chǔ)層疏干、水源枯竭、地面沉降等環(huán)境水文地質(zhì)問(wèn)題的發(fā)生具有重要作用。

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

4 結(jié)論

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

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

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

符號(hào)注釋?zhuān)?/p>

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