唐連松，陳 剛，胡 成

(中國地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北 武漢 430074)

在地下水封洞庫的建設(shè)和運(yùn)營中，地下水壓力的分布直接影響洞庫的水封效果，因此準(zhǔn)確地刻畫地下水滲流場對于保證地下水封洞庫的安全和穩(wěn)定運(yùn)行至關(guān)重要[1]。由于地下水封洞庫工程的復(fù)雜性，在其設(shè)計(jì)、施工和運(yùn)營階段，一般都需要借助數(shù)值模擬方法對地下洞庫所在水文地質(zhì)單元的地下水滲流場特征、水位動態(tài)變化及涌水量進(jìn)行模擬和預(yù)測。在數(shù)值模擬的過程中，模型參數(shù)的選取是決定模型計(jì)算準(zhǔn)確度的重要影響因素之一。以往的數(shù)值模擬工作中，往往是依據(jù)當(dāng)?shù)氐刭|(zhì)條件將研究區(qū)簡單地劃分為若干均質(zhì)區(qū)，分區(qū)內(nèi)部為均質(zhì)，分區(qū)間為非均質(zhì)，以整個(gè)分區(qū)為單個(gè)參數(shù)進(jìn)行參數(shù)估計(jì)，這種方法被稱為傳統(tǒng)參數(shù)分區(qū)法或常規(guī)分區(qū)法[2]，該方法雖然可以極大地減少模型參數(shù)的數(shù)量，但無法精確地刻畫含水層的滲透性分布，尤其是對于強(qiáng)烈非均質(zhì)各向異性的裂隙巖體含水層，由于裂隙巖體含水層的滲透系數(shù)表現(xiàn)出強(qiáng)烈的空間變異性，僅利用幾個(gè)鉆孔試驗(yàn)數(shù)據(jù)進(jìn)行常規(guī)統(tǒng)計(jì)，根本無法全面地反映出整個(gè)研究區(qū)含水層滲透性的分布特征。而采用參數(shù)反演分析法對影響地下水滲流場模擬的重要參數(shù)——滲透系數(shù)進(jìn)行優(yōu)化可以提高模擬的準(zhǔn)確度。反演分析法實(shí)際上是一種參數(shù)優(yōu)化方法，通過擬合地下水水位與實(shí)際地下水水位的接近程度來確定裂隙巖體含水層的滲透系數(shù)[3]，接近程度越高，說明擬合得到的滲透系數(shù)越精確。參數(shù)反演分析法分為直接反演法和間接反演法兩種[4]。直接反演法就是將實(shí)際的地下水水位分布作為已知數(shù)(已有觀測值)，將含水層滲透系數(shù)作為未知數(shù)進(jìn)行求解，當(dāng)模擬得出的地下水水位分布與實(shí)際地下水水位分布很接近時(shí)，求解完成。間接反演法是在前期水文地質(zhì)勘探及調(diào)查的基礎(chǔ)上，綜合研究區(qū)巖體結(jié)構(gòu)和水文地質(zhì)特性，從總體上確定裂隙巖體含水層的相對滲透性，并圈定滲透系數(shù)的大致范圍，再以此作為已知條件去擬合地下水水位分布，并通過比較地下水水位的計(jì)算值和實(shí)際值，逐步修正含水層滲透系數(shù)以達(dá)到較好的模擬效果[5]。

本文利用FEFLOW軟件建立了煙臺某地下水封洞庫運(yùn)營階段的研究區(qū)地下水滲流場數(shù)值模型，由于花崗巖體含水層的強(qiáng)烈非均質(zhì)性，其模擬結(jié)果與實(shí)測數(shù)據(jù)存在一定的誤差，尤其是洞庫地下水涌水量模擬值與實(shí)測數(shù)據(jù)有較大偏差，因此，進(jìn)一步以現(xiàn)場水力試驗(yàn)數(shù)據(jù)作為先驗(yàn)信息，利用地下壓力計(jì)的壓力監(jiān)測數(shù)據(jù)及工程結(jié)構(gòu)內(nèi)部的流量監(jiān)測數(shù)據(jù)，對建立的數(shù)值模型參數(shù)進(jìn)行了反演，在提高模型模擬精度的同時(shí)刻畫了強(qiáng)烈非均質(zhì)含水介質(zhì)的滲透系數(shù)場。

1 研究區(qū)概況

研究區(qū)位于山東半島北部、煙臺市西側(cè)，為某公司地下水封洞庫所在場地。該洞庫工程包括1個(gè)丙烷洞庫、1個(gè)丁烷洞庫、1個(gè)LPG洞庫，3個(gè)洞庫在平面上呈“品”字形分布。其中，位于北部的丙烷洞庫總庫容為50萬m3，主洞室所在高程為-120～-146 m(黃海高程)，水幕巷道所在高程為-94～-100 m(黃海高程)；位于西側(cè)的丁烷洞庫和位于東側(cè)的LPG洞庫總庫容均為25萬m3，主洞室所在高程為-90～-116 m(黃海高程)，水幕巷道所在高程為-64～-70 m(黃海高程)。

研究區(qū)含水介質(zhì)主要為第四系覆蓋層和中生代燕山早期中粗粒黑云母二長花崗巖，地下水的主要類型為松散巖類孔隙水、淺層基巖網(wǎng)狀裂隙水、深層基巖脈狀裂隙水。研究區(qū)地下水主要接受大氣降水的入滲補(bǔ)給，沿地下裂隙網(wǎng)絡(luò)向中西部低洼處徑流，排泄于第四系殘坡積或沖洪積層中，或被人工開采。研究區(qū)地下水動態(tài)監(jiān)測數(shù)據(jù)表明，淺層地下水隨季節(jié)動態(tài)變化明顯，深層裂隙水長期保持相對穩(wěn)定。研究區(qū)的基本地質(zhì)概況見圖1。

圖1 研究區(qū)地質(zhì)概況Fig.1 Geology introduction of the study area

2 數(shù)值模型的建立

2. 1 模型的邊界和數(shù)學(xué)模型

根據(jù)研究區(qū)的地質(zhì)條件，將研究區(qū)概化為非均質(zhì)各向異性三維潛水流，模型上邊界為入滲補(bǔ)給、蒸發(fā)排泄邊界；下邊界(標(biāo)高為-176 m)為隔水邊界，東部和西部邊界是沿流域分水嶺劃分的邊界，在模型中概化為隔水邊界；南部邊界部分為隔水邊界，部分為第二類邊界(定流量邊界)，流量根據(jù)區(qū)域地下水徑流模數(shù)計(jì)算給出；北部為黃海，是研究區(qū)的最低排泄基準(zhǔn)面，定為第一類邊界。洞室和水幕系統(tǒng)等工程結(jié)構(gòu)按第一類邊界處理。依照上述條件建立數(shù)學(xué)模型如下：

式中：Ω為研究區(qū)域；Kxx、Kyy、Kzz分別為x、y、z主方向的滲透系數(shù)(m/d)；H0為初始地下水水位(m)；H1為洞庫或水幕水位(m)；S0為地下水自由表面；S1為第一類邊界，表示黃海、河流、洞庫和水幕的位置；μs為單位儲水系數(shù)(1/m)；μd為給水度；ε為源匯項(xiàng);W為潛水面上的蒸發(fā)排泄，降雨入滲補(bǔ)給，井的抽水量和泉的排泄量等；S2為第二類邊界；q為第二類邊界的水流通量(m/d)。

2. 2 模型的網(wǎng)格劃分和參數(shù)賦值

2.2.1 模型的網(wǎng)格劃分

采用三角網(wǎng)格對模型進(jìn)行剖分，對洞室及水幕系統(tǒng)進(jìn)行局部加密，在平面上將4 993.7 m×5 165 m的研究區(qū)剖分為49 308個(gè)有限單元格，在垂向上按現(xiàn)場水力試驗(yàn)分段(現(xiàn)場壓水試驗(yàn)按10 m一段進(jìn)行)及洞庫工程結(jié)構(gòu)將研究區(qū)共劃分為18層，得到研究區(qū)地下水滲流場數(shù)值模型的三維結(jié)構(gòu)見圖2。

圖2 研究區(qū)地下水滲流場數(shù)值模型三維結(jié)構(gòu)Fig.2 Three-dimensional structure of the model in the study area

2.2.2 模型的參數(shù)賦值

根據(jù)研究區(qū)地質(zhì)圖、水文地質(zhì)圖，并結(jié)合現(xiàn)場水文地質(zhì)試驗(yàn)結(jié)果，對模型參數(shù)進(jìn)行賦值，其中第四系覆蓋層和巖體強(qiáng)風(fēng)化層即模型的第1層和第2層按參數(shù)分區(qū)法賦值，其參數(shù)分區(qū)情況見圖3。模型第1層為第四系覆蓋層，在第1層分區(qū)中,區(qū)域1、2為低山區(qū)域且第四系覆蓋層較薄，區(qū)域3、4為殘坡積物堆積形成的第四系覆蓋層，區(qū)域5、6為九曲河沖積物形成，區(qū)域7為九曲河河道，區(qū)域8、9為水田；模型第2層為巖體強(qiáng)風(fēng)化層，其中區(qū)域3的風(fēng)化程度稍強(qiáng)于區(qū)域1、2。模型中其他各層的參數(shù)利用現(xiàn)場鉆孔壓水試驗(yàn)結(jié)果的平均值進(jìn)行賦值。由此得到研究區(qū)地下水滲流場數(shù)值模型各層各分區(qū)具體參數(shù)的初始值見表1。

圖3 研究區(qū)地下水滲流場數(shù)值模型第1、2層參數(shù)分區(qū)圖Fig.3 Distribution of hydrogeological parameters in layer 1 and layer 2 of the numerical model of the groundwater seapage field in the study area

分層分區(qū)號水平滲透系數(shù)/(m·d-1)垂直滲透系數(shù)/(m·d-1)給水度單位儲水系數(shù)/m-111.555 200.155 5200.101×10-421.296 000.129 6000.101×10-4312.960 001.296 0000.201×10-4第四系覆48.640 000.864 0000.201×10-4蓋層(模50.728 000.072 8000.201×10-4型第1層)62.592 000.259 2000.251×10-4786.400 008.640 0000.301×10-4817.280 001.728 0000.251×10-499.750 000.975 0000.251×10-4強(qiáng)風(fēng)化層10.355 600.035 5600.025×10-5(模型第20.112 320.011 2320.025×10-52層)30.864 000.086 4000.055×10-5丁烷、LPG洞庫水幕層(模型第6層)0.003 000.000 3007×10-7丙烷洞庫水幕層(模型第10層)0.002 050.000 2057×10-7

2.3 地下洞庫及其周圍區(qū)域地下水壓力及水位的監(jiān)測

在地下洞庫運(yùn)營過程中，為了保證洞庫的水封條件長期有效，需要對地下洞庫及其周圍區(qū)域的地下水壓力進(jìn)行全方位監(jiān)測，因此在部分水幕孔中安裝有地下壓力計(jì)。地下壓力計(jì)安裝位置為3個(gè)洞庫的水幕層，其中丁烷、LPG洞庫水幕層(標(biāo)高為-64～-70 m)位于模型第6層，丙烷洞庫水幕層(標(biāo)高為-94～-100 m)位于模型第10層，模型在這兩個(gè)層位的初始水文地質(zhì)參數(shù)見表1。同時(shí)，將部分勘察鉆孔改造為地下水水位監(jiān)測井，對地下外圍的地下水水位進(jìn)行長期監(jiān)測。

2. 4 模型的模擬精度分析

利用讀數(shù)較為穩(wěn)定、數(shù)據(jù)較為可靠的地下壓力計(jì)及地下水水位監(jiān)測井(共31個(gè)監(jiān)測點(diǎn)位)獲得的數(shù)據(jù)對地下水滲流場數(shù)值模型計(jì)算結(jié)果進(jìn)行校驗(yàn)，其結(jié)果見圖4。

圖4 地下水水位模擬值與實(shí)測值的擬合情況Fig.4 Fitting of the simulated and measured groundwater level

由圖4可見，受研究區(qū)花崗巖含水層強(qiáng)烈非均質(zhì)性的影響，數(shù)值模型的模擬精度一般，地下水水位31個(gè)監(jiān)測點(diǎn)位模擬值與實(shí)測值的殘差平方和為264.86，地下水水位模擬值與實(shí)測值的最大誤差為7.4 m。

地下洞庫地下水涌水量模擬值與實(shí)測值的對比，見表2。

表2 洞庫地下水涌水量模擬值與實(shí)測值的對比(單位：m3/d)

由表2可知，地下洞庫地下水涌水量的模擬值與實(shí)測值也存在一定的誤差，尤其是丁烷洞庫地下水涌水量的模擬值與實(shí)測值的誤差較大，說明丁烷洞庫處可能存在滲透系數(shù)較大的強(qiáng)滲透帶，這種情況下需要對模型參數(shù)進(jìn)行進(jìn)一步的優(yōu)化以提高模型模擬的精度。

3 模型參數(shù)反演

3.1 PEST參數(shù)反演方法

由于該地下洞庫修建于較為完整的花崗巖體中，工程場地的地層巖性較為單一，而且工程建設(shè)會避開導(dǎo)水性較強(qiáng)或規(guī)模較大的地質(zhì)構(gòu)造，因此傳統(tǒng)的參數(shù)分區(qū)法很難應(yīng)用于這種情況。為此，本文采用PEST程序?qū)τ绊憯?shù)值模型計(jì)算結(jié)果的主要參數(shù)——滲透系數(shù)進(jìn)行參數(shù)反演，以提高模型模擬的精度。

PEST參數(shù)優(yōu)化算法是一種較為成熟、應(yīng)用較為廣泛的參數(shù)優(yōu)化方法。最初PEST是由Doherty[6]開發(fā)出的一種通用參數(shù)優(yōu)化算法，之后Tokin等[7]提出了一種結(jié)合Tikhonov正則化及T-SVD正則化方法優(yōu)點(diǎn)的混合正則化方法(SVD-Assist),進(jìn)一步優(yōu)化了PEST算法,使其更加完善。PEST算法在國內(nèi)外地下水?dāng)?shù)值模型參數(shù)優(yōu)化工作中取得了較好的效果，董艷輝等[8]利用并行化的PEST算法對大尺度的地下水?dāng)?shù)值模型進(jìn)行了參數(shù)優(yōu)化，得到了較為理想的結(jié)果;王禮恒等[9]利用假想模型進(jìn)一步驗(yàn)證了PEST算法在優(yōu)化地下水流動模型參數(shù)上的可靠性;姜蓓蕾等[10]對PEST算法所采用的向?qū)c(diǎn)-正則化方法進(jìn)行了深入的研究，證明了PEST是一種有效的參數(shù)反演方法，同時(shí)證明該方法中并不是向?qū)c(diǎn)越多反演精度越高，向?qū)c(diǎn)過多反而會影響計(jì)算精度且增加計(jì)算量。

PEST程序采用的是Gauss-Marquardt-Levenberg非線性參數(shù)估計(jì)算法，該算法是一種改進(jìn)的Gauss-Newton法[11]。該算法在每一次迭代過程中，計(jì)算模型結(jié)果對各待估參數(shù)的導(dǎo)數(shù)，利用這些導(dǎo)數(shù)將非線性模型線性化，并更新一次參數(shù)向量，計(jì)算新的目標(biāo)函數(shù)后，再進(jìn)入下一次迭代過程，通過反復(fù)迭代直至目標(biāo)函數(shù)達(dá)到最小值。

實(shí)際工作中，大部分模型都是非線性的，假設(shè)模型的n個(gè)初始參數(shù)值儲存在向量x0中，模型的m個(gè)實(shí)際觀測值儲存在向量y0中，則x0與y0的關(guān)系為

y0=M(x0)

(1)

根據(jù)泰勒公式，任意一組參數(shù)向量x與其對應(yīng)的觀測向量y之間的關(guān)系可以表示為

(2)

當(dāng)x→x0時(shí)，忽略高階無窮小項(xiàng)，則式(2)可簡化為

y=y0+J(x-x0)

(3)

式中:J為m行n列的偏導(dǎo)數(shù)雅克比矩陣:

(4)

由此可得出非線性模型參數(shù)反演過程中的目標(biāo)函數(shù)為

Φ=(y-y0-J(x-x0))TQ(y-y0-J(x-x0))

(5)

式中:Q為m行n列的實(shí)際觀測值的權(quán)重矩陣，參數(shù)優(yōu)化過程就是不斷迭代更新參數(shù)向量x使得目標(biāo)函數(shù)Φ最小化的過程[12]。

3. 2 參數(shù)反演范圍

圖5 模型參數(shù)反演范圍Fig.5 Range of parameter inversion

由于前期勘察工作及洞庫運(yùn)營階段的監(jiān)測工作都是圍繞著洞庫所在區(qū)域開展的，缺少洞庫工程外圍的實(shí)測資料，因此參數(shù)反演的范圍為以洞庫區(qū)域?yàn)橹行牡? 000 m×1 000 m正方形區(qū)域(見圖5)。地下水水位監(jiān)測點(diǎn)和巖體滲透性原位試驗(yàn)點(diǎn)主要集中在丁烷、LPG洞庫水幕層(-64～-70 m)和丙烷洞庫水幕層(-94～-100 m)這兩個(gè)關(guān)鍵層位，這是由于在洞庫運(yùn)營過程中，水幕巷道及水幕孔的持續(xù)供水對保證洞庫水封效果起著至關(guān)重要的作用，而且在深部修建這種特殊的工程結(jié)構(gòu)也會使這兩個(gè)層位的地下水動態(tài)變得更加復(fù)雜，對這兩個(gè)層位的模型參數(shù)進(jìn)行優(yōu)化是客觀評價(jià)洞庫水封效果、提高水幕系統(tǒng)工作效率的必要工作，因此本次以巖體滲透系數(shù)原位試驗(yàn)點(diǎn)和地下水水位監(jiān)測點(diǎn)的監(jiān)測數(shù)據(jù)對丁烷、LPG洞庫水幕層和丙烷洞庫水幕層這兩個(gè)關(guān)鍵層位的巖體滲透系數(shù)(K)進(jìn)行參數(shù)反演，丁烷、LPG洞庫水幕層和丙烷洞庫水幕層的原位試驗(yàn)點(diǎn)、監(jiān)測點(diǎn)及水位分布情況見圖6和圖7。

圖6 丁烷、LPG洞庫水幕層的原位試驗(yàn)點(diǎn)、監(jiān)測點(diǎn) 及水位分布Fig.6 In-situ test points,monitoring points and ground- water level distribution in layer of water curtains of Butane cavern and LPG cavern

圖7 丙烷洞庫水幕層的原位試驗(yàn)點(diǎn)、監(jiān)測點(diǎn)及水位分布Fig.7 In-situ test points,monitoring points and ground- water level distribution in layer of water curtain of Propane cavern

3. 3 參數(shù)反演結(jié)果和模型優(yōu)化結(jié)果

運(yùn)行PEST程序?qū)Χ⊥椤PG洞庫和丙烷洞庫水幕層的巖體滲透系數(shù)進(jìn)行參數(shù)反演后，得到洞庫巖體滲透系數(shù)的反演結(jié)果(以lgK表示)以及地下水水位模擬值與實(shí)測值的擬合結(jié)果見圖8和圖9。丁烷、LPG洞庫地下水水位模型模擬值與實(shí)測值的殘差平方和為39.82，丙烷洞庫地下水水位模擬值與實(shí)測值的殘差平方和為19.15。

由圖8和圖9可見，經(jīng)參數(shù)反演后的丁烷、LPG洞庫水幕層的巖體滲透系數(shù)最大值為0.006 1 m/d，最小值為0.000 03 m/d，丙烷洞庫水幕層的巖體滲透系數(shù)最大值0.059 m/d，最小值為0.000 105 m/d，其滲透系數(shù)在空間上表現(xiàn)出強(qiáng)烈的非均質(zhì)性。

圖8 丁烷、LPG洞庫水幕層巖體滲透系數(shù)反演結(jié)果(以lgK表示)以及地下水水位模擬值與實(shí)測值的擬合結(jié)果Fig.8 Parameter inversion result of layer of water curtains of Butane cavern and LPG cavern and the fitting result of the groundwater level simulated values and the meosured values

圖9 丙烷洞庫水幕層巖體滲透系數(shù)反演結(jié)果(以lgK表示)以及地下水水位模擬值與實(shí)測值的擬合結(jié)果Fig.9 Parameter inversion result of layer of the water curtain of Propane cavern and the fitting result of the groundwater level simulated values and the measured values

假設(shè)巖體滲透系數(shù)服從對數(shù)正態(tài)分布，將其反演結(jié)果繪制成巖體滲透系數(shù)概率密度曲線，同時(shí)將現(xiàn)場16個(gè)鉆孔中對應(yīng)層位的壓水試驗(yàn)數(shù)據(jù)與巖體滲透系數(shù)概率密度曲線進(jìn)行對比，其結(jié)果見圖10。

由圖10可見，丁烷、LPG洞庫和丙烷洞庫水幕層兩個(gè)層位巖體滲透系數(shù)的反演結(jié)果與現(xiàn)場實(shí)測數(shù)據(jù)具有相似的統(tǒng)計(jì)特征，說明其反演結(jié)果能夠較為客觀地反映研究區(qū)的實(shí)際情況。根據(jù)概率密度函數(shù)可知，丁烷、LPG洞庫水幕層巖體的滲透系數(shù)以0.000 5 m/d為主，丙烷洞庫水幕層的滲透系數(shù)以0.003 9 m/d為主。

圖10 丁烷、LGP洞庫和丙烷洞庫水幕層巖體滲透系數(shù) 反演結(jié)果與壓水試驗(yàn)數(shù)據(jù)的對比Fig.10 Comparison of inversion results of rock mass permeability coefficients of Butane,LGP and Propane cavern with pressure water test data

利用參數(shù)優(yōu)化后的數(shù)值模型計(jì)算得到的各洞庫地下水涌水量見表3。

表3 參數(shù)優(yōu)化后洞庫地下水涌水量模擬值與實(shí)測值的對比(單位：m3/d)

由表3可知，與初始數(shù)值模型的計(jì)算結(jié)果(見表2)相比較，參數(shù)優(yōu)化后的數(shù)值模型計(jì)算結(jié)果與實(shí)測數(shù)據(jù)更加相符，說明經(jīng)過參數(shù)優(yōu)化后的數(shù)值模型具有更高的擬合精度。

4 結(jié)論與展望

結(jié)合研究區(qū)水文地質(zhì)條件分析和現(xiàn)場水文地質(zhì)試驗(yàn)數(shù)據(jù)，建立了煙臺某地下水封洞庫運(yùn)營階段研究區(qū)地下水滲流場數(shù)值模型，通過對比監(jiān)測數(shù)據(jù)發(fā)現(xiàn)該均質(zhì)數(shù)值模型的計(jì)算結(jié)果與實(shí)測數(shù)據(jù)有較大偏差，本研究進(jìn)一步采用參數(shù)反演與現(xiàn)場試驗(yàn)數(shù)據(jù)相結(jié)合的方法來刻畫地下水封洞庫特定層位的巖體滲透系數(shù)，經(jīng)過參數(shù)優(yōu)化后，數(shù)值模型地下洞庫擬合精度得到提高，模擬結(jié)果表明：丁烷、LPG洞庫水幕層的巖體滲透系數(shù)以0.000 5 m/d為主，滲透性最強(qiáng)處巖體滲透系數(shù)為0.006 1 m/d，但這一層位上沒有形成明顯的強(qiáng)滲透帶；丙烷洞庫水幕層的巖體滲透系數(shù)以0.003 9 m/d為主，滲透性最強(qiáng)滲透系數(shù)為0.059 m/d，這一層位上出現(xiàn)了明顯的強(qiáng)滲透帶，這也解釋了為什么丁烷洞庫從開始運(yùn)行起就出現(xiàn)了地下水涌水量偏大的現(xiàn)象。由此說明，經(jīng)參數(shù)反演后的地下水滲流場數(shù)值模型能夠更加客觀、準(zhǔn)確地評價(jià)水封洞庫的運(yùn)行狀況。由于工程實(shí)際情況較為復(fù)雜，本文只對巖體滲透系數(shù)這一個(gè)參數(shù)進(jìn)行了反演分析，實(shí)際上數(shù)值模型的計(jì)算結(jié)果往往是由模型計(jì)算所用到的所有參數(shù)共同決定的，因此數(shù)值模型的參數(shù)反演工作還有待進(jìn)一步完善，對數(shù)值模型所用到的所有參數(shù)進(jìn)行綜合分析與協(xié)同反演將是未來研究的發(fā)展趨勢。