基于最小二乘-分步解耦的泄?jié)B漏帶與隧道圍巖滲流參數(shù)方法

張 鵬， 李彥霈， 田世寬， 萬 飛， 王占強(qiáng)， 胡 波

(1.南京工業(yè)大學(xué)交通運(yùn)輸工程學(xué)院， 南京 210009； 2.中交一公局第四工程有限公司， 南寧 530000； 3.交通運(yùn)輸部公路科學(xué)研究所， 北京 100088)

涌水突泥是巖溶區(qū)隧道工程的主要地質(zhì)災(zāi)害之一。巖溶巖體滲透能力主要取決于連接孔隙空間的通道大小[1-2]。由于巖溶發(fā)育具有隨機(jī)不均勻性，溶腔孔隙及其連接通道的分布、規(guī)模和組合導(dǎo)致其呈現(xiàn)顯著的隨機(jī)性、不均勻性[3]，從而致使工程中很難依據(jù)鉆孔資料對巖溶隧道涌水量做出可靠估計。盡管已有研究表明，巖溶巖體涌水量主要受巖腔及其周邊巖體的滲透系數(shù)控制[4]。然而，目前還沒有一種準(zhǔn)確、可靠地確定巖溶溶腔滲流性能的試驗(yàn)方法或解析方法[5]。

現(xiàn)場抽水試驗(yàn)和滲流參數(shù)反演是當(dāng)前確定巖體滲透參數(shù)的常用方法。盡管現(xiàn)場抽水試驗(yàn)是確定滲透參數(shù)最為經(jīng)典水文地質(zhì)方法之一[6]，但其以多孔介質(zhì)的達(dá)西定律為基礎(chǔ)，普遍存在數(shù)據(jù)離散和代表性不強(qiáng)以及費(fèi)用高等問題，且其對于高度非均勻非均質(zhì)的喀斯特巖溶的適用性尚未得到公認(rèn)[7-9]。同其他巖土工程問題一樣，采用基于實(shí)測數(shù)據(jù)進(jìn)行巖體滲流參數(shù)反演是確定滲透參數(shù)較為理想的方法。紀(jì)成亮等[10]利用壓水試驗(yàn)數(shù)據(jù)反演了裂隙巖體滲透參數(shù)；王媛等[11]通過動態(tài)觀測總結(jié)水頭、位移等多類型資料，建立滲流場與應(yīng)力場在裂隙巖體中動態(tài)全耦合的參數(shù)反演思路；張振昌等[12]對胡麻嶺隧道7號豎井進(jìn)行了野外現(xiàn)場定流量抽水試驗(yàn)，以實(shí)測數(shù)據(jù)為依據(jù),應(yīng)用布爾頓配線法、直線法以及水位恢復(fù)法分別得到了第三系粉細(xì)砂巖滲透系數(shù)；唐連松等[13]為得到模型中兩個關(guān)鍵層位的巖體滲透系數(shù)場，以PEST(parameter estimation)參數(shù)優(yōu)化算法對地下水封洞庫低滲巖體滲透系數(shù)進(jìn)行反演,；Gan等[14]基于非穩(wěn)定滲流試驗(yàn)成果提出了一種裂隙巖體滲流參數(shù)的綜合反分析技術(shù)；Zhao等[15]基于礦井涌水事件反演了巖體的滲透系數(shù)；Liu等[16]基于達(dá)西定律反演分析了煤巖的氣滲透系數(shù)；Zhang等[17]采用微震數(shù)據(jù)反演了頁巖的滲透系數(shù)；Liu等[18]反演預(yù)測了熱裂干巖的滲透率變化；Fang等[19]采用混合遺傳算法反演了巖體的孔隙率、飽和度及滲透系數(shù)；Wang等[20]基于流動電流檢測反演了巖體滲透系數(shù)。李永壽等[21]采用已有的損傷演化方程,通過滲透系數(shù)與損傷之間的定量關(guān)系,來刻畫巖體擾動過程中滲透系數(shù)的變化規(guī)律。盡管上述研究中有部分成果還處于理論研究階段，但亦有許多已應(yīng)用于實(shí)際工程中，這表明采用反演分析確定巖溶巖體滲透系數(shù)是可行的。

圍繞雞冠山隧道涌水?dāng)?shù)據(jù)資料，首先釆用電法和地質(zhì)雷達(dá)探測方法查明了隧道上覆巖溶溶腔的分布狀態(tài)，從而結(jié)合工程地質(zhì)勘察資料構(gòu)建了上覆天坑降雨入滲的滲流地質(zhì)模型，然后采用基于最小二乘-分步耦合優(yōu)化方法的巖溶隧道圍巖滲流參數(shù)反演法，以3個涌水事件的涌水量、涌水速率、涌水與降雨關(guān)聯(lián)時間、天坑積水水頭為反演邊界條件，關(guān)聯(lián)性解耦下泄?jié)B漏帶與隧道圍巖滲流規(guī)律，最終通過分步優(yōu)化探討了下泄?jié)B漏帶和隧道圍巖的滲透系數(shù)，以期為工程安全性評價提供參數(shù)取值依據(jù)。

1 工程概況與模型構(gòu)建

1.1 雞冠山隧道概況

雞冠山隧道位于貴州省威寧縣爐山鎮(zhèn)境內(nèi)，為典型巖溶地層條件下的特長分離式越嶺隧道，左線隧道起訖里程樁號為ZK70+238～ZK74+010，長3 772 m，右線隧道起訖里程樁號為YK70+235～YK74+035，長3 800 m，隧道總體呈S形展布，洞身主體呈北東-南西向展布，洞身走向方位角為248°,進(jìn)出口段均略呈弧形展布，其中進(jìn)口段路線走向約274°，出口段走向約265°。隧道建筑限界為10.25 m×5 m(寬×高)，設(shè)計縱坡1.4%，時速度80 km/h，采用燈光照明，機(jī)械通風(fēng)。

雞冠山隧道區(qū)域構(gòu)造屬溶蝕峰叢洼地地貌，峰叢基座相連，溶丘與洼地相間展布，漏斗、落水洞等地貌發(fā)育,如圖1所示。穿越區(qū)溶蝕丘峰峰頂多呈長橢圓-渾圓狀，山脊狹長，脊頂?shù)匦屋^緩，受本區(qū)構(gòu)造控制，多呈近東西向展布。雞冠山隧道穿越段自然地面標(biāo)高2 040～2 244 m。巖溶天坑中心樁號為K72+280,沿隧道方向范圍為K72+220～K72+340,距離隧道洞頂76 m,地層走向和隧道軸線之間大角度斜交。隧道區(qū)出露基巖為二疊系下統(tǒng)茅口組(P1m)灰?guī)r、二疊系下統(tǒng)棲霞組(P1q)灰?guī)r、二疊系下統(tǒng)梁山組(P1l)泥巖、石炭系上統(tǒng)馬平組(C3mp)灰?guī)r，石炭系中統(tǒng)黃龍組(C2hl)灰?guī)r。其中，進(jìn)洞口地層為茅口組(P1m)灰?guī)r，巖層產(chǎn)狀為15°～35°∠21°～27°，出洞口地層為石炭系中統(tǒng)黃龍群(C2hl),巖層產(chǎn)狀為73°～80°∠32°～49°。區(qū)內(nèi)陡傾構(gòu)造及溶蝕節(jié)理裂隙發(fā)育，節(jié)理主要以走向NE35°～45°、傾向北西和南東、傾角65°～85°和走向NW30°～65°、傾向北東和南西、傾角60°～80°的張裂隙為主，裂隙線密度可達(dá)2～4條/m，延長2 m、切深3～5 m，裂面起伏，透水性好，是地表水下滲的重要通道。

圖1 天坑、落水洞Fig.1 Karst pit and catavothre

雞冠山隧道位于亞熱帶季風(fēng)濕潤氣候區(qū)，雨量充沛，氣候宜人。據(jù)威寧氣象站資料，多年平均降水量為960.6 mm，蒸發(fā)量1 407 mm，日最大降水量為166 mm，極端最高氣溫35.7 ℃，最低-13.8 ℃。全年平均日照百分率為33%。5—10月份為降雨高發(fā)季，11月—次年4月份為枯季，西部降水量一般比東部多150 mm，氣溫比東部低3～4 ℃。兩者高差差距較大，地域差異性在氣溫、降水量等氣象要素上較為明顯。

1.2 巖溶隧道模型構(gòu)建

為進(jìn)一步明確雞冠山隧道上覆巖溶發(fā)育情況，采用高密度電法和地質(zhì)雷達(dá)探查了天坑落水洞至隧道圍巖之間的下泄?jié)B漏帶分布狀況，然后再綜合區(qū)域工程地質(zhì)條件實(shí)現(xiàn)了巖溶隧道地質(zhì)模型的構(gòu)建。

1.2.1 下泄?jié)B漏帶探查

下泄?jié)B漏帶探查采用高密度電法和地質(zhì)雷達(dá)進(jìn)行。其中，高密度電法采用四極法裝置(電極排列)2維測量，反演深度約70 m。地質(zhì)雷達(dá)采用的是意大利IDS公司的探地雷達(dá)，在隧道襯砌內(nèi)壁掃描，每個測區(qū)測線長度約100 m，每個測區(qū)測線6條，共3個測區(qū)。測線布置如圖2所示。

圖2 勘探區(qū)域及測線布置Fig.2 Exploration area and survey line layout

根據(jù)高密度電法及探地雷達(dá)探測數(shù)據(jù)對雞冠山隧道下泄?jié)B漏帶分布情況進(jìn)行可視化處理，如圖3所示。可知，測區(qū)內(nèi)存在兩條較大斷層SE22°傾角52°，SE23°傾角59°，存在多處溶腔向下延伸至隧道涌水點(diǎn)，形成多處貫通性涌水帶，結(jié)合地表洞室分布形成了空間多處通道，隧道建設(shè)及運(yùn)營存在較大安全隱患。

圖3 測區(qū)空穴、破碎帶、斷層隨隧道里程空間分布Fig.3 Distributions of karst cave, fracture zone, fault in the space

1.2.2 數(shù)值分析模型及邊界條件

研究內(nèi)容主要基于涌水資料數(shù)據(jù)以及積水水頭邊界條件，反演隧道上覆圍巖的綜合滲透系數(shù)。因此，模擬采用二維平面，隧道縱向取1 m。根據(jù)以往的經(jīng)驗(yàn)和精度要求，隧道計算模型根據(jù)上一節(jié)地質(zhì)物探成果及雞冠山隧道設(shè)計的實(shí)際開挖尺寸確定建立。以雞冠山隧道ZK72+220～ZK72+340下穿巖溶天坑段落水洞位置截取二維剖面，構(gòu)建數(shù)值模型，隧道圍巖賦予綠色屬性，簡化后的下泄?jié)B漏帶賦予藍(lán)色屬性，計算模型如圖4所示。該模型寬320 m，高185 m，采用四邊形網(wǎng)格剖分，共有網(wǎng)格單元7 021個，節(jié)點(diǎn)14 008個。

圖4 計算模型Fig.4 Computational model

考慮到雞冠山隧道巖溶天坑段上覆巖體存在落水洞，且溶腔周邊存在斷層及破碎帶，降雨積水主要沿巖體溶蝕通道及裂隙擴(kuò)散，因此，采用達(dá)西飽和滲流模型開展數(shù)值模擬。為方便試算，參考地質(zhì)手冊中建議的關(guān)于巖溶地區(qū)巖體滲透系數(shù)的取值范圍1×10-4～1×10-2cm/s，巖溶溶腔下滲通道初始滲透系數(shù)取5×10-3cm/s，隧道圍巖區(qū)初始滲透系數(shù)取5×10-4cm/s。襯砌的滲透系數(shù)參數(shù)依據(jù)工程設(shè)計取1.0×10-10cm/s。模型邊界條件：x方向的邊界x=0和x=320固定，y方向的邊界y=0和y=1固定，z方向的邊界為z=0，z方向上邊界延伸至地表。模型左右及底邊界默認(rèn)為不透水固定位移邊界，地表水頭邊界以天坑匯水計算確定的天坑地表降雨積水深度為準(zhǔn)。并參考涌水報告考慮地下水位情況，模型初始設(shè)置以模型底邊界向上40 m作為地下水位線。

2 最小二乘-分步耦合優(yōu)化反演

基于以上數(shù)值模型的建立，擬采用反演分析方法確定巖溶巖體滲透系數(shù)。當(dāng)前，巖土體滲流參數(shù)識別主要可分為直接法和間接法兩類[22]。其中，直接法主要是借助邊界條件求解基于控制方程(Darcy定律)的含待識別參數(shù)的線性一階偏微分方程來獲得滲透參數(shù)，但其前提是整個滲流研究區(qū)域的水頭導(dǎo)數(shù)必須已知；間接法則是通過假定滲透參數(shù)初值，然后通過正演分析獲得滲流特征量，再將其與實(shí)測值進(jìn)行對比，如不符合實(shí)際情況，則調(diào)整參數(shù)，重新計算，直至符合要求。顯然，間接分析法可適用于監(jiān)測數(shù)據(jù)較少的情況，且不需要計算監(jiān)測點(diǎn)數(shù)據(jù)的微分。當(dāng)然，間接分析法也存在求極小化殘差過程的非線性問題，但該問題可通過許多優(yōu)化方法得以解決，如經(jīng)典高斯-牛頓法、正則化最小二乘法、最小二乘支持向量機(jī)、人工神經(jīng)網(wǎng)絡(luò)等。其中，最小二乘是通過最小化誤差平方和尋找數(shù)據(jù)的最佳函數(shù)匹配，較為方法簡便。

最小二乘法目標(biāo)函數(shù)為

minE=[hc-hm]T[hc-hm]

(1)

式(1)中：hc為計算矢量，hc={hc1,hc2,…,hcm}；hm為觀測矢量，hm={hm1,hm2,…,hmM},M為觀測點(diǎn)數(shù)目與觀測次數(shù)之積。若假設(shè)T為所求參數(shù)矢量，T={T1,T2,…,TL},L為要識別參數(shù)維數(shù)，對于非約束極小化問題，其高斯-牛頓迭代公式為

TK+1=TK-DK,K=1,2,3,…

(2)

AKDK=GK,K=1,2,3,…

(3)

式中：AK=JTJ(L×L),GK=JK[hc-hm](L×1)；J為M×L階Jacobian矩陣；D為Gauss-Newton方向矢量。

需要注意的是，當(dāng)系數(shù)矩陣A奇異時，D矩陣構(gòu)造將破壞，當(dāng)系數(shù)矩陣A病態(tài)時，式(3)的解將嚴(yán)重失真。采用正則化處理可有效解決這一問題。正則化最小二乘方程為

(AK+_I)DK=GK

(4)

式(4)中：_為正則化因子；I為單位矩陣。

圖5 正則化最小二乘反分析流程Fig.5 Flow chart of regularized least square inverse analysis

因此，采用基于最小二乘-分步耦合優(yōu)化方法的巖溶隧道圍巖滲流參數(shù)反演法。以3個涌水事件的涌水量、涌水速率、涌水與降雨關(guān)聯(lián)時間、天坑積水水頭為反演邊界條件，降雨發(fā)生至涌水出現(xiàn)的間隔時間與隧道涌水速率作為反演搜索的關(guān)鍵指標(biāo)，關(guān)聯(lián)性解耦下泄?jié)B漏帶與隧道圍巖滲透參數(shù)。下泄?jié)B漏帶滲透參數(shù)主要控制降雨發(fā)生至隧道涌水出現(xiàn)的間隔時間；下泄?jié)B漏帶滲透參數(shù)與隧道圍巖區(qū)滲透參數(shù)的綜合滲透效應(yīng)控制隧道涌水速率。故反演的流程可按照兩階段展開。

3 巖體滲透參數(shù)反演分析

3.1 參數(shù)反演依據(jù)

依據(jù)隧道涌水與降雨關(guān)系分析，下穿巖溶天坑雞冠山隧道段的降雨下滲路徑與隧道涌水時間主要由下泄?jié)B漏帶滲透系數(shù)來控制。反演邊界條件如表1所示，并以降雨發(fā)生至涌水出現(xiàn)的間隔時間與隧道涌水量作為關(guān)鍵反演指標(biāo)。圖6為涌水流量與壓力監(jiān)測點(diǎn)位置。

表1 隧道下穿天坑洞段的涌水事件反演指標(biāo)表Table 1 Inversion index of water inrush event of tunnel underlying karst sinkhole

圖6 左線水頭壓力監(jiān)測點(diǎn)示意圖Fig.6 Schematic diagram of monitoring points about water pressure on the left tunnel

為綜合考慮降雨涌水時間間隔和隧道涌水量影響，關(guān)聯(lián)性解耦下泄?jié)B漏帶與隧道圍巖滲透參數(shù)，從而采用兩步反演法開展巖溶巖體滲透參數(shù)分析，即首先以降雨發(fā)生至涌水出現(xiàn)的間隔時間為判定標(biāo)尺，設(shè)計大量計算樣本開展正演分析，初步確定各地層滲透系數(shù)反演參數(shù)取值區(qū)間；然后以隧道涌水速率為判定標(biāo)尺，分步優(yōu)化所得下泄?jié)B漏帶和隧道圍巖區(qū)的滲透參數(shù)。具體流程如圖7所示。

圖7 滲透系數(shù)反演分析流程圖Fig.7 Flow chart of permeability coefficient inversion analysis

3.2 涌水規(guī)律分析

隧道涌水規(guī)律與上部巖溶地層性質(zhì)密切相關(guān)，尋找合理滲透參數(shù)取值區(qū)間則是獲得理想反演結(jié)果的前提。本文隧道涌水滲流數(shù)值模擬采用有限差分FLAC3D軟件進(jìn)行，地表水頭邊界以天坑匯水計算確定的天坑地表積水深度為準(zhǔn)。考慮到實(shí)際降雨半小時隧道即出現(xiàn)涌水，因此，以0.5 h作為判定標(biāo)尺。

分別對照涌水事件1、2、3開展隧道涌水滲流數(shù)值模擬。隧道左洞監(jiān)測點(diǎn)水頭壓力曲線如圖8所示，處于下泄?jié)B漏帶的左線右拱腰位置的水頭壓力在0.5 h左右率先開始增加，其他監(jiān)測點(diǎn)位置的水頭壓力在隨后才開始增加，因此，以左線右拱腰位置的水頭壓力開始增加的時間作為降雨發(fā)生至涌水出現(xiàn)的間隔時間。

圖8 涌水事件隧道左線監(jiān)測點(diǎn)水頭壓力曲線圖Fig.8 Changes of pore water pressure at monitoring points in left tunnel during water gushing

圖9所示為3個涌水事件滲流5 h后的涌水壓力分布，顯然，處于下泄?jié)B漏帶側(cè)的隧道涌水壓力要大于另一側(cè)。圖10所示為3個涌水事件隧道涌水量速率曲線。經(jīng)多次試算，計算峰值涌水速率與實(shí)際涌水峰值吻合時的下泄?jié)B漏帶滲透系數(shù)在1×10-3～3×10-3cm/s，隧道圍巖的滲透系數(shù)在1×10-4～3×10-4cm/s。

圖9 涌水事件水頭壓力分布Fig.9 Distribution of pore water pressure during different water gushing events

圖10 涌水事件涌水速率曲線Fig.10 Curves of water gushing rate

3.3 滲透系數(shù)反演

為進(jìn)一步保證反演參數(shù)的準(zhǔn)確性，采用逐步掃描法來優(yōu)化滲透參數(shù)取值。簡單起見，以計算下滲時間與實(shí)際涌水事件中隧道收集到水的時間，數(shù)值模擬的隧道左線涌水速率與實(shí)際涌水事件中涌水速率的差值構(gòu)建目標(biāo)函數(shù)，即

(5)

(6)

式中：k1為下泄?jié)B漏帶滲透系數(shù)；k2為隧道圍巖滲透系數(shù)；N為樣本總數(shù)；xt為第t次涌水事件的計算下滲時間；Xt為第t次實(shí)際涌水事件中隧道收集到水的時間；yt為第t次涌水事件涌水速率；Yt為第t次實(shí)際涌水事件中單位長度左洞涌水速率。

3.3.1 下泄?jié)B漏帶滲透系數(shù)反演

以涌水事件1為始，下泄?jié)B漏帶滲透系數(shù)基準(zhǔn)值取1×10-3cm/s，然后同步增加或減少k1，每次計算步長為1×10-4cm/s，最終獲得下泄?jié)B漏帶滲透系數(shù)k1與f(k1)的關(guān)系曲線如圖11(a)所示，對應(yīng)左線右拱腰監(jiān)測點(diǎn)涌水壓力如圖12(a)所示。顯然，目標(biāo)函數(shù)值隨下泄?jié)B漏帶滲透系數(shù)增加先降后升。當(dāng)目標(biāo)函數(shù)最小時，所取的下泄?jié)B漏帶滲透系數(shù)與實(shí)際滲透系數(shù)越接近，因此，下泄?jié)B漏帶滲透系數(shù)最終可確定為1.3×10-3cm/s。同理，可以涌水事件2和3進(jìn)一步修正下泄?jié)B漏帶滲透系數(shù)。以涌水事件1所得下泄?jié)B漏帶滲透系數(shù)取值為基準(zhǔn)，以1×10-4cm/s為掃描步長進(jìn)行數(shù)值模擬，所得涌水事件2和3下泄?jié)B漏帶滲透系數(shù)k1與f(k1)的關(guān)系曲線如圖11(b)和圖11(c)所示。對應(yīng)的左線右拱腰監(jiān)測點(diǎn)涌水壓力曲線如圖12(b)和圖12(c)所示。可知，下泄?jié)B漏帶滲透系數(shù)分別為1.4×10-3cm/s和4.5×10-3cm/s。為減小誤差，對3個涌水事件反演出的下泄?jié)B漏帶滲透系數(shù)求平均值，計算所得滲透系數(shù)2.4×10-3cm/s，即下泄?jié)B漏帶滲透系數(shù)。

圖11 目標(biāo)函數(shù)與反演滲透系數(shù)關(guān)系Fig.11 Relationship between objective function and permeability coefficient of water infiltration zone

圖12 左線右拱腰監(jiān)測點(diǎn)孔壓變化曲線Fig.12 Water gushing pressure of monitoring points at the right arch lumbar of the left tunnel

3.3.2 反演隧道圍巖滲透系數(shù)

采用相同分析步驟，開展隧道圍巖滲透系數(shù)反演。具體實(shí)施時，依次涌水事件1、2、3為依據(jù)，隧道圍巖滲透系數(shù)基準(zhǔn)值取1.3×10-4cm/s，然后同步增加或減少k2，每次計算步長為0.2×10-4cm/s，獲得隧道圍巖滲透系數(shù)k2與f(k2)的關(guān)系曲線和對應(yīng)單位長度隧道峰值監(jiān)測涌水速率曲線，如圖13、圖14所示，確定三次涌水事件所對應(yīng)的隧道圍巖滲透系數(shù)如表2所示。為減小誤差，對3個涌水事件反演出隧道圍巖滲透系數(shù)求平均值，所得滲透系數(shù)3.2×10-4cm/s，即為隧道圍巖滲透系數(shù)。

圖13 目標(biāo)函數(shù)與隧道圍巖滲透系數(shù)關(guān)系Fig.13 Relationship between objective function and permeability coefficient of surrounding rock mass

圖14 單位長度隧道峰值監(jiān)測涌水速率Fig.14 Peak water gushing rate in the left tunnel within unit length

表2 隧道圍巖滲透系數(shù)表Table 2 Table of permeability coefficient of surrounding rock of tunnel

4 結(jié)論

雞冠山巖溶隧道涌水量主要受下泄?jié)B漏帶的滲透系數(shù)控制，合理確定隧道圍巖滲透參數(shù)是估計隧道涌水量與防排水設(shè)計的關(guān)鍵。釆用綜合勘察法得到雞冠山隧道區(qū)工程地質(zhì)條件和水文地質(zhì)條件、巖溶天坑的下泄?jié)B漏帶位置，構(gòu)建飽和達(dá)西滲流計算模型，然后，以3個涌水事件對應(yīng)的積水水頭邊界、降雨涌水時間差及涌水速率為依據(jù)，采用最小二乘-分步耦合優(yōu)化方法反演巖體滲流參數(shù)，最終獲得下泄?jié)B漏帶滲透系數(shù)建議取值為2.4×10-3cm/s、隧道圍巖滲透系數(shù)取3.2×10-4cm/s，可為后續(xù)降雨條件下巖溶隧道襯砌結(jié)構(gòu)安全性分析研究提供參數(shù)取值依據(jù)。