劉杰，宋美鈺，于彥，邊宗斌，康楠，秦莉紅

（1.天津地?zé)峥辈殚_發(fā)設(shè)計院，天津300250；2.天津地質(zhì)資料館，天津30000）

天津市地?zé)豳Y源廣泛應(yīng)用于供暖、洗浴、生活和工農(nóng)業(yè)用水等諸多領(lǐng)域。自上世紀(jì)八十年代以來，地?zé)豳Y源的開發(fā)利用進(jìn)入了快速發(fā)展時期，開發(fā)利用形式也從粗放型逐漸發(fā)展到集約化階段。面對地?zé)豳Y源需求的快速增長和開采規(guī)模的不斷擴(kuò)大，如何維持可持續(xù)開發(fā)利用，為地?zé)豳Y源最優(yōu)化開發(fā)利用的采灌布局提供依據(jù)，是目前我市地?zé)嵝袠I(yè)急需解決的問題。本文以天津潘莊凸起區(qū)為例，說明應(yīng)用優(yōu)化管理模型進(jìn)行地?zé)豳Y源優(yōu)化管理。

1 模型建立及識別

1.1 水文地質(zhì)模型

選?、艏墭?gòu)造單元潘莊凸起區(qū)為本次的研究區(qū)，具體邊界：北部至漢沽斷裂，西部以天津斷裂為界，東部邊界為滄東斷裂[1]，南部為自然流量邊界，總面積374 km2（圖1）。研究區(qū)內(nèi)主要以基巖裂隙型薊縣系霧迷山組熱儲層為主要開采層位，共有地?zé)峋?8眼。

研究區(qū)內(nèi)霧迷山組熱儲層普遍分布，巖性以灰褐色含灰泥晶白云巖、泥晶白云巖為主，含石英粉砂泥晶砂屑。鉆井揭露頂板埋深為1 752～2 016 m，揭露厚度為480～1 032 m，揭露的最大深度3 634 m，尚無鉆孔鉆穿，故為了模型計算，熱儲層底板統(tǒng)一以4 000 m計算。

根據(jù)本區(qū)的地?zé)岬刭|(zhì)、水文地質(zhì)條件和地?zé)崃黧w開發(fā)利用特點、同位素分析結(jié)果：計算區(qū)西側(cè)天津斷裂為隔水邊界、但不隔熱；北部漢沽斷裂為定流量流入邊界；東部滄東斷裂為熱流體導(dǎo)水導(dǎo)熱通道；南部為自然流量邊界[2]。

1.2 數(shù)學(xué)模型

對所研究的流體系統(tǒng)封閉邊界Гn包圍的任意子域Vn積分。TOUGH2求解的質(zhì)量和能量平衡方程的一般形式可以寫為：

式中：MK代表單位體積的質(zhì)量或能量（κ=1、2、…… NK，標(biāo)記物質(zhì)組分水、空氣、CO2、溶質(zhì)、……）；φ代表孔隙度；sβ代表狀態(tài)β的飽和度（即β態(tài)占據(jù)的毛孔體積部分）；ρβ代表狀態(tài)β的密度（kg/代表組分κ出現(xiàn)在狀態(tài)β中質(zhì)量比；F代表表示質(zhì)量或能量通量（流體流量）；q代表源匯項；n代表表面單元dГn上的指向dVn的法向矢量。同時根據(jù)Darcy定律多相形式給出各狀態(tài)的通量：

式中：ρβ代表狀態(tài)β的密度(kg/m3)；υβ代表狀態(tài)β的Darcy速率（體積通量）（m/s）；k代表絕對滲透率（mD）；krβ代表狀態(tài)β的相對滲透率（mD）；μβ代表粘度（kg/m/s）；Pβ代表狀態(tài)β中的流體壓力，它是相對其他狀態(tài)（通常參照氣態(tài)）壓力P和毛細(xì)壓力Pcβ（≤0）的總和；g代表重力加速度向量。

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

1.3.1 研究區(qū)剖分、初始條件及時間離散

根據(jù)霧迷山組熱儲層的埋藏分布條件、地質(zhì)構(gòu)造分布特征、富水性特征及地?zé)峋牟季趾退粍討B(tài)監(jiān)測情況[3]，利用TOUGH2.0模擬軟件將研究區(qū)剖分為4 580個單元（圖2）。

研究區(qū)有較系統(tǒng)的水位常觀資料始于2002年，所以以2002年1月為模型模擬的初始時刻。

擬合時段分供暖期（120天）與非供暖期（245天），時間步長為30天，模型識別時間到2014年10月。

1.3.2 邊界補(bǔ)給量

邊界的具體補(bǔ)給條件根據(jù)水位等值線圖及歷史水位擬合確定，相應(yīng)的單寬補(bǔ)給量可根據(jù)水位等值線圖計算，并在擬合過程中調(diào)整（表1）。

圖2 研究區(qū)剖分圖Fig.2 Subdivision graph maps in the study area

1.3.3 參數(shù)介質(zhì)分區(qū)

計算水文地質(zhì)參數(shù)選用是依據(jù)三維模型的計算需要，主要采取分區(qū)賦值的方法。根據(jù)霧迷山組熱儲層的埋藏分布和控?zé)針?gòu)造分布特征，利用所收集的降壓試驗資料進(jìn)行綜合統(tǒng)計分析，按所求水文地質(zhì)參數(shù)進(jìn)行初始分區(qū)。此次共收集9眼霧迷山組地?zé)峋祲涸囼灁?shù)據(jù)（表2），根據(jù)資料[4]將該組熱儲層水文地質(zhì)參數(shù)劃分為3個初始分區(qū)（圖2）。

1.3.4 模型識別與驗證

模型擬合的開始時間選擇在2009年10月，模擬終止于2014年10月，劃分為61個時間步長（每個月為一個時間步長）。

對霧迷山熱儲層的觀測水位按液面20℃溫度進(jìn)行統(tǒng)一換算，利用模型計算出的水位值與觀測值進(jìn)行擬合對比（圖3）?？梢钥闯觯河嬎阒蹬c觀測值在變化趨勢和數(shù)值都很接近，大部分實際觀測值點都在模擬曲線上[5]，計算值與觀測值平均相對誤差DL-19為1.38%；DL-19B為2.39%；DL-40為2.12%；DL-40B為1.98%，計算值和觀測值較為一致。

表1 邊界單寬補(bǔ)給量(m3/d·m)Tab.1 Single wide recharge boundary(m3/d·m)

表2 熱儲模型主要參數(shù)表Tab.2 Main parameter of the thermal reservoir model

從上述綜合分析可以看出，研究區(qū)數(shù)值模擬模型所對應(yīng)的地?zé)岬刭|(zhì)模型與實際地?zé)岬刭|(zhì)模型是基本吻合的，同時也說明該模型能基本上合理地反映示范區(qū)的地?zé)岬刭|(zhì)條件，為下一步建立地?zé)豳Y源管理模型奠定基礎(chǔ)。

2 管理模型的建立

研究區(qū)地?zé)豳Y源優(yōu)化管理的目標(biāo)[6]：一是在盡可能滿足地?zé)嵘a(chǎn)需求的條件下，實現(xiàn)全區(qū)水位降幅最??；二是針對目前的地?zé)衢_發(fā)利用現(xiàn)狀，盡量減少地?zé)峋g相互干擾，充分發(fā)揮模型調(diào)控功能。

2.1 目標(biāo)函數(shù)的建立

研究區(qū)地?zé)豳Y源管理面對的是整個地?zé)豳Y源系統(tǒng)，目的是規(guī)劃、調(diào)整地?zé)豳Y源系統(tǒng)的時空人為輸出，使地?zé)崃黧w輸出趨近于既定目標(biāo)[7]，同時使地?zé)豳Y源系統(tǒng)的天然補(bǔ)給朝著有利于實現(xiàn)系統(tǒng)功能的方向發(fā)展以使有限的地?zé)豳Y源發(fā)揮最大的經(jīng)濟(jì)、環(huán)境和社會效益。

故取示范區(qū)中地?zé)峋拈_采/回灌量為模型的決策變量，即：該模型的目標(biāo)函數(shù)為使n年的規(guī)劃階段、m個地?zé)峋牡責(zé)豳Y源總開采量最大：

式中：maxZ代表最大地?zé)豳Y源總開采量；m代表地?zé)峋當(dāng)?shù)量；n代表規(guī)劃階段（年）；Q代表地?zé)峋_采量。

本次規(guī)劃階段分為30年，其中計算A1階段（1年）、A5階段（5年）和A30階段（30年），共有12眼地?zé)峋畢⑴c計算，即可得總開采量為：

式中：maxZ代表最大地?zé)豳Y源總開采量；Q代表地?zé)峋_采量。

2.2 約束條件的建立

2.2.1 約束條件的選擇

約束條件首要的是水位降深約束[8]，要求各井的水位降深不能超過某個給定的值。如果是給定的值，該值要考慮三個因素：一是不應(yīng)超過地?zé)嵋?guī)劃中規(guī)定的熱儲層中水位最大埋深；二是經(jīng)濟(jì)效益，開采井降深過大會導(dǎo)致能耗過大；三是不應(yīng)超過抽水試驗達(dá)到的最大降深。

圖3 研究區(qū)地?zé)峋宦裆睿?0℃）計算值與實測值擬合圖Fig.3 Water depth(20℃)calculated values and measured values fitting figure of the geothermal wells in the study area

式中：β（i，j，t）代表單位脈沖響應(yīng)函數(shù)，表示由在第j點脈沖的作用，在第t時段末第i點的水位降深響應(yīng)；Q（i，j）代表第t時刻作用于j點施加的脈沖量(開采/回灌量)；Si代表第i點的水位埋深約束。

②水位降深極限約束：

式中：S(i，t)代表第t時刻i點的水位降深；S0(i)代表第i點的水位降深約束。

③決策變量上下限約束：

④非負(fù)性約束，表示各決策變量為非負(fù)：區(qū)內(nèi)的地?zé)豳Y源開發(fā)利用更趨于采灌平衡狀態(tài)（控制每年降幅不大于1 m）。

式中：Q(j，t)代表第t時刻作用于j點施加的脈沖量（開采/回灌量）；S(i，t)代表第t時刻i點的水位降深。

2.2.2 約束條件的具體數(shù)值

（1）水位埋深約束

根據(jù)天津市地?zé)豳Y源管理辦法，霧迷山組熱儲水位埋深總體不超過200 m。

（2）決策變量上下限約束

決策變量上下限約束是指在滿足生產(chǎn)需求的前提下，對開采量進(jìn)行合理的上下限控制。針對示范區(qū)內(nèi)12眼參與計算的地?zé)峋M(jìn)行調(diào)查，根據(jù)近三年的開發(fā)利用情況制定約束數(shù)值（表3）。

（3）水位降深約束

水位降深約束條件指在分布于研究區(qū)內(nèi)的18個水位降深約束點處，把每個約束點的水位降深約束在一定水平上的一組數(shù)據(jù)。以基本維持目前水位為基準(zhǔn)，不超過控制開采區(qū)標(biāo)準(zhǔn)為要求，并且分別計算水位降深為0 m/a，1 m/a，3 m/a的優(yōu)化方案。

為了下文方便描述，將水位降深約束為0 m/a定為優(yōu)化方案一、水位降深約束為1 m/a定為優(yōu)化方案二、水位降深約束為3 m/a定為優(yōu)化方案三。為了體現(xiàn)優(yōu)化管理模型在地?zé)豳Y源管理上的管理能力，另外選取一種方案（優(yōu)化方案四），即保持示范區(qū)內(nèi)開采總量和回灌總量不變的情況下，通過管理模型調(diào)整各個地?zé)峋拈_采量或回灌量，進(jìn)而讓整個示范

表3 決策變量約束統(tǒng)計表Tab.3 Statistical data of the decision variables constraints

3 計算方法與程序設(shè)計

3.1 計算方法

本項目采用Matlab編程進(jìn)行優(yōu)化計算[9]，對管理模型進(jìn)行最優(yōu)求解。利用Matlab編程來進(jìn)行線性優(yōu)化，采用Matlab中單純形法（Linprog函數(shù)）進(jìn)行求解[10]。

3.2 程序設(shè)計

現(xiàn)以第一時段（即2015年管理模型）為例，以下為Matlab程序設(shè)計：

Linprog線性規(guī)劃函數(shù)

x=Linprog（f，A，b）解如下形式的線性規(guī)劃問題

minfx

x

Subject to：Ax≤b

x=Linprog（f，A，b，vlb，vub）參數(shù)vlb、vub給出設(shè)計變量（開采/回灌量）的上下邊界約束，即vlb≤x≤vub。

x=Linprog（f，A，b，vlb，vub，x0）設(shè)置初始值

[x,fval,exitflag,output,lambda]=linprog（c,a,b,[],[],vlb,vub）

返回目標(biāo)函數(shù)最優(yōu)解x，和在x處的值：fval=fT*x，是否存在exitflag標(biāo)志，優(yōu)化解結(jié)構(gòu)體output，拉格朗日乘子結(jié)構(gòu)體lambda。

4 優(yōu)化結(jié)果及模型檢驗

4.1 優(yōu)化結(jié)果

對12眼地?zé)峋煌瑮l件下各時段分別建立MATLAB公式進(jìn)行求解，結(jié)果見表4。

表4 不同優(yōu)化方式的總量（×104m3）Tab.4 The quantity of different ways of the optimization（×104m3）

從表4可以看出：1）2014年薊縣系霧迷山組熱儲層年開采量為213.22×104m3，各時段的開采量基本能滿足示范區(qū)地?zé)豳Y源的開發(fā)利用及管理模型建立的限制條件。2）在該區(qū)域每年降深為0 m或1 m時，優(yōu)化管理后回灌量都比開采量大，只有當(dāng)降深為3 m的時候，優(yōu)化管理后的開采量比回灌量大，證明該區(qū)域地?zé)豳Y源開發(fā)利用相對密集，地?zé)峋嗷ブg影響較大，必須增加回灌來保證水位降深不斷變小。3）從表中可以看出水位降深越大，優(yōu)化開采量越大；在同一開采方式下，年份越長，優(yōu)化開采量也越大。所以可以利用管理模型對區(qū)域內(nèi)地?zé)豳Y源進(jìn)行資源調(diào)控，對降深大的區(qū)域按照優(yōu)化管理后的開采量進(jìn)行控制，把優(yōu)化開采量定為年開采指標(biāo)。

4.2 模型檢驗

根據(jù)優(yōu)化方案四（以2014年開采/回灌量為依據(jù)，1 m/a水位降深為約束）得出的優(yōu)化方案以及現(xiàn)有地?zé)豳Y源開發(fā)利用情況，基于建立好的TOUGH2.0數(shù)值模擬平臺對霧迷山組熱儲層壓力場動態(tài)變化情況進(jìn)行預(yù)測，具體預(yù)測數(shù)值如表5。

分別從優(yōu)化方案四中取出相同的四眼地?zé)峋―L-19、DL-19B、DL-40、DL-40B），與按2014年開采方案開采預(yù)測值進(jìn)行對比（圖4）。

從上述數(shù)據(jù)可以看出：按優(yōu)化方案四（水位降幅約束為1 m，按2014年開采/回灌量）計算的30年優(yōu)化方案代入數(shù)值模型中，出來的結(jié)果每眼地?zé)峋骄瞪钭畲鬄?.8 m，最小為0.6 m。

圖4 優(yōu)化方案四與現(xiàn)用方案對比曲線圖Fig.4 Comparison chart of the optimized scheme four and current solutions

綜合分析，將三個優(yōu)化方案結(jié)果分別代入數(shù)值模型進(jìn)行模擬，得出每眼地?zé)峋畬嶋H地?zé)崃黧w水位降深與優(yōu)化方案前提降深值基本吻合，說明該管理模型能合理計算出優(yōu)化開采量，比維持現(xiàn)狀開采更為合理。

表5 按優(yōu)化方案四計算的各年份數(shù)值模擬結(jié)果Tab.5 The numerical simulation from the optimized scheme four calculating results in each year

5 結(jié)論

天津地?zé)豳Y源優(yōu)化管理模型的研究表明：面對天津地?zé)豳Y源需求的快速增長和開采規(guī)模的不斷擴(kuò)大，對其進(jìn)行優(yōu)化管理是使其能夠可持續(xù)利用的有效方法。

科學(xué)調(diào)度合理開采地?zé)豳Y源，天津地區(qū)地?zé)豳Y源豐富，但如何可持續(xù)開采面臨很大困難。在開采地?zé)崃黧w保證環(huán)境安全的同時，嚴(yán)格控制地區(qū)水位埋深，使地?zé)豳Y源開發(fā)利用達(dá)到“采補(bǔ)平衡”。今后地?zé)豳Y源開發(fā)利用都按本研究對各區(qū)的約束條件進(jìn)行控制，在對約束條件論證后管理部門應(yīng)考慮制定有關(guān)條例規(guī)定，保證天津市地?zé)豳Y源更精細(xì)化管理。

[1]天津市地質(zhì)礦產(chǎn)局.天津市區(qū)域地質(zhì)志[M].北京：地質(zhì)出版社，1992.

[2]GB/T 11615-2010地?zé)豳Y源地質(zhì)勘查規(guī)范[S].

[3]胡燕，林黎，等.天津市地?zé)豳Y源可持續(xù)開發(fā)潛力評價報告[R].天津:天津地?zé)峥辈殚_發(fā)設(shè)計院，2015.

[4]宗振海，劉杰，等.2015年天津地?zé)豳Y源開發(fā)利用動態(tài)監(jiān)測項目[R].天津:天津地?zé)峥辈殚_發(fā)設(shè)計院，2015.

[5]寧立波，董少剛，馬傳明，等.地下水?dāng)?shù)值模擬的理論與實踐[M].武漢：中國地質(zhì)大學(xué)出版社，2010.

[6]國家技術(shù)監(jiān)督局.地下水資源管理模型工作要求(GB/T 14497-93)[Z].1993.

[7]舒艷，王紅旗，等.地下水管理模型研究進(jìn)展[J].水文地質(zhì)工程地質(zhì)，2005，32（6）:85-90.

[8]盧文喜.地下水系統(tǒng)的模擬預(yù)測和優(yōu)化管理[M].北京：科學(xué)出版社，1999.

[9]王沫然.MATLAB與科學(xué)計算[M].北京：電子工業(yè)出版社，2005，第2版.

[10]尹尚先.地下水管理模型中響應(yīng)函數(shù)解法探討[J].水文地質(zhì)工程地質(zhì)，2001，29（2）：33-35.