基于最小二乘-分步解耦的泄滲漏帶與隧道圍巖滲流參數方法

張 鵬， 李彥霈， 田世寬， 萬 飛， 王占強， 胡 波

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

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

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

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

1 工程概況與模型構建

1.1 雞冠山隧道概況

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

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

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

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

1.2 巖溶隧道模型構建

為進一步明確雞冠山隧道上覆巖溶發育情況，采用高密度電法和地質雷達探查了天坑落水洞至隧道圍巖之間的下泄滲漏帶分布狀況，然后再綜合區域工程地質條件實現了巖溶隧道地質模型的構建。

1.2.1 下泄滲漏帶探查

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

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

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

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

1.2.2 數值分析模型及邊界條件

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

圖4 計算模型Fig.4 Computational model

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

2 最小二乘-分步耦合優化反演

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

最小二乘法目標函數為

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

(1)

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

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方向矢量。

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

(AK+_I)DK=GK

(4)

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

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

因此，采用基于最小二乘-分步耦合優化方法的巖溶隧道圍巖滲流參數反演法。以3個涌水事件的涌水量、涌水速率、涌水與降雨關聯時間、天坑積水水頭為反演邊界條件，降雨發生至涌水出現的間隔時間與隧道涌水速率作為反演搜索的關鍵指標，關聯性解耦下泄滲漏帶與隧道圍巖滲透參數。下泄滲漏帶滲透參數主要控制降雨發生至隧道涌水出現的間隔時間；下泄滲漏帶滲透參數與隧道圍巖區滲透參數的綜合滲透效應控制隧道涌水速率。故反演的流程可按照兩階段展開。

3 巖體滲透參數反演分析

3.1 參數反演依據

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

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

圖6 左線水頭壓力監測點示意圖Fig.6 Schematic diagram of monitoring points about water pressure on the left tunnel

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

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

3.2 涌水規律分析

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

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

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

圖9所示為3個涌水事件滲流5 h后的涌水壓力分布，顯然，處于下泄滲漏帶側的隧道涌水壓力要大于另一側。圖10所示為3個涌水事件隧道涌水量速率曲線。經多次試算，計算峰值涌水速率與實際涌水峰值吻合時的下泄滲漏帶滲透系數在1×10-3～3×10-3cm/s，隧道圍巖的滲透系數在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 滲透系數反演

為進一步保證反演參數的準確性，采用逐步掃描法來優化滲透參數取值。簡單起見，以計算下滲時間與實際涌水事件中隧道收集到水的時間，數值模擬的隧道左線涌水速率與實際涌水事件中涌水速率的差值構建目標函數，即

(5)

(6)

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

3.3.1 下泄滲漏帶滲透系數反演

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

圖11 目標函數與反演滲透系數關系Fig.11 Relationship between objective function and permeability coefficient of water infiltration zone

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

3.3.2 反演隧道圍巖滲透系數

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

圖13 目標函數與隧道圍巖滲透系數關系Fig.13 Relationship between objective function and permeability coefficient of surrounding rock mass

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

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

4 結論

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