吳奎鋒，王福亮，丁建芳，趙香玲,寇春陽

(1.中鐵西南科學(xué)研究院有限公司，四川 成都 611731；2.四川綿九高速公路有限責(zé)任公司，四川 綿陽 621700)

1 引 言

隨著我國交通設(shè)施建設(shè)的加快，將會(huì)修建更多的隧道工程，含水構(gòu)造超前探測是保障隧道施工安全的重要環(huán)節(jié)，在隧道含水體的超前預(yù)報(bào)物探方法中，應(yīng)用最為廣泛的是地震波反射法(Tunnel Seismic Prediction,TSP)[1]、電磁波反射法、直流電法[2-4]、瞬變電磁法和復(fù)頻電導(dǎo)法(Complex Frequency Conductivity,CFC)[5,6]。激發(fā)極化法因受隧道環(huán)境限制，相對于其他方法，在隧道中的應(yīng)用研究較晚；德國Geohydrauuc Data公司推出了一種全新聚焦電流頻率域激電隧道超前預(yù)報(bào)方法(Bore-tunning Electrical Ahead Monitoring,BEAM)[7]，該技術(shù)應(yīng)用在TBM環(huán)境下，通過探測儀器、傳感器與TBM裝備集成，實(shí)現(xiàn)了自動(dòng)測量，它是一種充分利用電流同性相斥的原理，通過外圍環(huán)狀電極的頻域誘導(dǎo)極化特征來探測掌子面前方的地質(zhì)異常體，這種方法在歐洲許多國家得到應(yīng)用，但在國內(nèi)報(bào)道或應(yīng)用相對較少；山東大學(xué)李術(shù)才院士團(tuán)隊(duì)[8,9]研發(fā)的以激發(fā)極化法為理論基礎(chǔ)的隧道超前探測系統(tǒng)，其從視電阻率值和半衰時(shí)差著手，通過三維數(shù)據(jù)反演預(yù)測前方含水異常體，該方法在國內(nèi)部分隧道中得到應(yīng)用，并取得了一定的成果。

前人的研究工作取得了階段性的成果，但也存在如下問題：一是激發(fā)極化法三維反演時(shí)，大多采用的反演思路均先進(jìn)行電阻率參數(shù)反演，再依據(jù)電阻率反演結(jié)果進(jìn)行極化率反演，致使反演結(jié)果過多依賴于電阻率反演[10]；二是在進(jìn)行電阻率和極化率參數(shù)反演時(shí)，未考慮兩種物性參數(shù)之間的耦合關(guān)系。針對上述問題，本文對反演目標(biāo)函數(shù)進(jìn)行優(yōu)化，采用有限內(nèi)存擬牛頓法設(shè)計(jì)快速有效的自動(dòng)反演方法，反演過程中同時(shí)迭代更新電阻率和極化率兩種物性參數(shù)；并依據(jù)電阻率和極化率之間的結(jié)構(gòu)耦合關(guān)系，引入交叉梯度函數(shù)約束兩種參數(shù)反演[11]，提高反演解譯的精度與準(zhǔn)確性，最后通過理論模型合成數(shù)據(jù)驗(yàn)證正反演算法的可行性和有效性。

2 激發(fā)極化法探測理論及隧道環(huán)境探測裝置選擇

激發(fā)極化法超前探測是以巖礦石的激發(fā)極化效應(yīng)為基礎(chǔ)，通過觀測隧道位置的視電阻率和視極化率變化規(guī)律來探測掌子面前方地質(zhì)情況的一種地質(zhì)預(yù)報(bào)方法[12-14]，由于節(jié)理裂隙、斷層破碎帶、溶洞等易富水地質(zhì)構(gòu)造的富水量、分布規(guī)模等對圍巖的激發(fā)極化效應(yīng)影響較大，依據(jù)激發(fā)極化法探測富水地質(zhì)構(gòu)造具有天然的優(yōu)勢，地層含水量越大，地層激電效應(yīng)越強(qiáng)，其極化電場越明顯、越容易探測，進(jìn)而通過分析獲取掌子面前方含水地層的位置和規(guī)模等情況[15]。

激發(fā)極化法在隧道地質(zhì)預(yù)報(bào)探測中，針對隧道全空間環(huán)境，應(yīng)對探測布極、聚焦、供電、數(shù)據(jù)處理、解釋等技術(shù)進(jìn)行改進(jìn)，方能實(shí)現(xiàn)隧道超前預(yù)報(bào)的目的。傳統(tǒng)激發(fā)極化法探測裝置有：對稱四級測深裝置、中間梯度裝置、固定電源排列裝置、聯(lián)合剖面裝置等；由于隧道環(huán)境下可布極空間有限，測線只能沿隧道洞身方向布設(shè)；中間梯度裝置、對稱四級裝置在隧道環(huán)境下無法布設(shè)；聯(lián)合剖面裝置沿隧道洞身方向移動(dòng)電極，主要觀測隧道洞身段旁側(cè)地質(zhì)信息，而隧道超前地質(zhì)預(yù)報(bào)主要是探測掌子面前方地質(zhì)信息，因此適合隧道環(huán)境探測裝置只能選擇定源排列裝置。

定源排列裝置可分為定源二級裝置(AM)和定源三級裝置(AMN)[16,17]，定源三級裝置相較于定源二級裝置，具有縱向分辨率高，受旁側(cè)地質(zhì)信息影響小的特點(diǎn)，本文選用定源三級裝置作為觀測裝置[18]，測線沿隧道走向布置在隧道底板或邊墻，供電點(diǎn)源A極固定在掌子面，B極放置在無窮遠(yuǎn)處，測量電極M、N沿測線移動(dòng)以采集數(shù)據(jù)，在隧道底板和左右邊墻上布置多條測線，一方面多條測線數(shù)據(jù)之間可相互校核，提高數(shù)據(jù)的可信度，另一方面增加了觀測數(shù)據(jù)量，有利于抑制反演的多解性，每條測線可反映不同位置的異常信息，有利于實(shí)現(xiàn)掌子面前方異常體的三維成像[19]。

3 三維正演數(shù)值模擬

3.1 正演理論

隧道施工激發(fā)極化法三維正演數(shù)值模擬就是在已知圍巖電阻率和極化率分布的情況下，求解隧道腔體接收器位置的視電阻率和視極化率，三維介質(zhì)正演問題無解析解，只能通過數(shù)值模擬方法來進(jìn)行求解，正演從點(diǎn)電源微分方程出發(fā)，采用矩形網(wǎng)格離散地下介質(zhì)，并在邊界處設(shè)置邊界條件，最后求解大型稀疏矩陣線性方程組。含隧道腔體地質(zhì)模型如圖1所示。點(diǎn)電源A位于隧道掌子面處，電性參數(shù)Ω1介質(zhì)中存在電性參數(shù)為Ω2的不均勻體，?！薇硎緹o窮遠(yuǎn)處邊界，Γ表示兩種介質(zhì)分界面。

對于一次電位求解，在均勻全空間模型中，若忽略地表反射界面對地下介質(zhì)電場的影響或點(diǎn)電源距離地表足夠遠(yuǎn)，則一次電位可表示為：

(1)

式中，u0為一次電位，單位為V;I為電流強(qiáng)度，單位為A;r為場點(diǎn)到源點(diǎn)的距離，單位為m;σ0為均勻全空間電導(dǎo)率，單位為S/m。對于二次電位的求解，采用異常電位法求解含隧道腔體的全空間三維點(diǎn)源場[20]，二次電位所滿足的微分方程如下：

?·(σ?u)=-?·[(σ-σ0)?u0]

(2)

式中，u為二次電位，單位為V;σ為異常體電導(dǎo)率，單位為S/m。離散區(qū)域邊界可分為內(nèi)邊界和外邊界，在計(jì)算過程中內(nèi)邊界條件會(huì)自動(dòng)滿足，外邊界是指網(wǎng)格剖分區(qū)域模擬無窮遠(yuǎn)處的邊界及地表與空氣界面，在這些邊界處采用第三類邊界條件。

(3)

式中，n為邊界外法線向量;n1和n2分別為Ω1和Ω2的外法線方向;u1和u2分別為Ω1和Ω2的電位，單位為V。

考慮到計(jì)算精度和計(jì)算效率，選取有限差分法進(jìn)行數(shù)值模擬，有限差分法解直流電法正演問題時(shí)，首先需要離散偏微分方程式(2)，再通過數(shù)值方法求解離散后的方程，從而得到地下的電位分布。直流電法正演方程可表示為如下線性方程組：

Ku=-(K-K0)u0

(4)

其中，K是一個(gè)大型稀疏對稱正定的系數(shù)矩陣;K0為一次電位全空間系數(shù)矩陣;u是剖分區(qū)域各節(jié)點(diǎn)的電位向量;u0為均勻全空間一次電位。

采用穩(wěn)定雙共軛梯度法(BICGSTAB法)求解式(4)線性方程組，即可得到所有網(wǎng)格節(jié)點(diǎn)處的二次電位u。通過整理電位或者電位差，即可得到視電阻率參數(shù)：

(5)

式中，ρa(bǔ)為視電阻率，單位為Ω·m;K為裝置系數(shù);Um為接收點(diǎn)M的電位(單接收極)，單位為V;ΔUmn為接收電極M和N之間的電位差，單位為V。

激發(fā)極化法正演包含兩部分內(nèi)容，即無激電效應(yīng)所產(chǎn)生的場和含激電效應(yīng)產(chǎn)生的場的疊加，在進(jìn)行激發(fā)極化法正演模擬時(shí)，需進(jìn)行兩次直流電法正演，將式(2)中電導(dǎo)率σ替換成帶激電效應(yīng)的電導(dǎo)率：

ση=σ(1-η)

(6)

其中，η是地下地質(zhì)體的極化率，單位為C·m2/V;ση表示帶激電效應(yīng)的電導(dǎo)率，單位為S/m。

利用直流電法的正演方法，先正演計(jì)算不含激電效應(yīng)的視電阻率，再正演計(jì)算含激電效應(yīng)的視電阻率，最后得到視極化率：

(7)

其中，ηa是地下地質(zhì)體的視極化率，單位為C·m2/V;ρη和ρ分別為含激電效應(yīng)和不含激電效應(yīng)接收點(diǎn)處的視電阻率，單位為Ω·m。

3.2 隧道腔體影響分析

在全空間模型中，植入一高阻體來模擬隧道腔體，模型參數(shù)如下：全空間默認(rèn)電阻率為1 000 Ω·m，隧道腔體默認(rèn)電阻率為1.0×106Ω·m；掌子面前方為均勻介質(zhì)；采用定源三極裝置，點(diǎn)源A位于掌子面，接收極M、N沿隧底均勻布極，極距MN=3 m，收發(fā)距AO從3 m至117 m,以3 m為點(diǎn)距均勻增大。通過正演計(jì)算即可得到圖2中純隧道腔體的視電阻率響應(yīng)曲線和純隧道腔體的視極化率響應(yīng)曲線。

圖2 純隧道腔體正演響應(yīng)對比曲線Fig.2 Contrastive curve of forward response of pure tunnel cavity

對比視電阻率響應(yīng)曲線和視極化率響應(yīng)曲線可得出如下結(jié)論：

1)均勻全空間視電阻率響應(yīng)基本與圍巖電阻率一致，在收發(fā)距較小時(shí)略有所偏差，為場源效應(yīng)所致。

2)存在隧道腔體時(shí)，視電阻率曲線受隧道腔體影響嚴(yán)重，視電阻率數(shù)據(jù)無法反映圍巖的真實(shí)電阻率；視極化率曲線不受隧道腔體影響，視極化率數(shù)據(jù)能反映圍巖的真實(shí)極化率。

3.3 正演模型算例

隧道施工中，富水溶洞是常見的致災(zāi)地質(zhì)構(gòu)造，設(shè)定隧道前方為電性介質(zhì)均勻的全空間，默認(rèn)為電阻率為1 000 Ω·m，在均勻空間中植入低電阻率高極化率異常體，如圖3所示，異常體模型參數(shù)如下：電阻率為10 Ω·m，極化率為20 %，來模擬隧道中富水溶洞構(gòu)造。

圖3 富水溶洞模型Fig.3 Water-rich karst cave model

在點(diǎn)電源A的激發(fā)下，通過正演計(jì)算出所有節(jié)點(diǎn)的視電阻率和視極化率，選取其中一剖面(Z=0 m)，其視電阻率和視極化率響應(yīng)如圖4所示。

圖4 富水溶洞模型正演響應(yīng)Fig.4 Forward response of water-rich karst cave model

從視電阻率響應(yīng)曲線和視極化率響應(yīng)可以看出，視電阻率及視極化率對三個(gè)異常體都有明顯的響應(yīng)，且對隧道開挖方向(X方向)分辨率較高；視電阻率響應(yīng)存在局部的旁側(cè)效應(yīng)，出現(xiàn)高阻假異常，視極化率響應(yīng)不存在旁側(cè)效應(yīng)。

4 三維反演數(shù)值模擬

4.1 反演理論

從視電阻率和視極化率曲線很難直觀得到異常體的信息，需要借助反演技術(shù)對異常體進(jìn)行解釋，地球物理反演是依據(jù)觀測數(shù)據(jù)去反推地球物理模型的過程，在三維反演中，比較常見的反演方法有最小二乘法、非線性共軛梯度法(Non-linear Conjugate Gradient, NLCG)、OCCAM法、有限內(nèi)存擬牛頓法(L-BFGS)等。

表1 不同反演方法優(yōu)缺點(diǎn)對比

本文選用有限內(nèi)存擬牛頓法(L—BFGS)實(shí)現(xiàn)電阻率和極化率雙參數(shù)反演[11]。三維激發(fā)極化法反演目標(biāo)函數(shù)可以表示為：

(8)

Wd=diag{1/σ1,1/σ2,...1/σj,...1/σm}

(9)

其中，σj為第j個(gè)觀測數(shù)據(jù)的誤差。引入數(shù)據(jù)協(xié)方差矩陣主要有兩個(gè)目的：

1)不同偏移距時(shí)，對數(shù)據(jù)有一定的歸一化作用，可以有效地平衡不同偏移距的觀測數(shù)據(jù)，以達(dá)到穩(wěn)定反演的效果。

2)數(shù)據(jù)質(zhì)量參差不齊時(shí)，可以降低質(zhì)量較差觀測數(shù)據(jù)的權(quán)重，避免反演過度的去擬合質(zhì)量較差的數(shù)據(jù)。

在多參數(shù)反演中，介質(zhì)的物性變化規(guī)律可以用一個(gè)梯度矢量進(jìn)行描述，既包含變化數(shù)值的信息，又含有方向信息。對于不同物性參數(shù)的反演，提出了一種基于交叉梯度約束的多參數(shù)反演算法，在反演的目標(biāo)函數(shù)中加入交叉梯度項(xiàng)來進(jìn)行模型參數(shù)的耦合約束。在模型更新過程中，兼顧數(shù)據(jù)擬合、模型正則化、不同的物性參數(shù)等信息同步迭代下去，直到達(dá)到反演終止條件。

三維模型交叉梯度函數(shù)的定義如下：

(10)

式中，mσ代表電導(dǎo)率參數(shù)；mη代表極化率參數(shù)。交叉梯度函數(shù)有以下幾點(diǎn)性質(zhì)：

1)兩種電性參數(shù)進(jìn)行交叉梯度計(jì)算時(shí)，當(dāng)任意一種電性參數(shù)不發(fā)生變化時(shí)，交叉梯度值為零。

2)交叉梯度值不受電性參數(shù)數(shù)值改變幅度的影響。

3)交叉梯度值與兩種電性參數(shù)結(jié)構(gòu)耦合性相關(guān)，物性參數(shù)變化方向不同時(shí)，交叉梯度值不為零，結(jié)構(gòu)越相似，交叉梯度數(shù)值越小。

在式(8)引入交叉梯度項(xiàng)后的目標(biāo)函數(shù)：

(11)

對應(yīng)的電導(dǎo)率參數(shù)目標(biāo)函數(shù)梯度可表示為：

(12)

對應(yīng)的極化率參數(shù)目標(biāo)函數(shù)梯度可表示為：

(13)

交叉梯度項(xiàng)關(guān)于電導(dǎo)率的偏導(dǎo)數(shù)只與極化率參數(shù)有關(guān)，關(guān)于極化率參數(shù)的偏導(dǎo)數(shù)只與電導(dǎo)率有關(guān)，目標(biāo)函數(shù)梯度直接影響海森矩陣的計(jì)算及模型更新的步長及方向等因素，引入交叉梯度項(xiàng)，可以讓兩種參數(shù)在反演過程中結(jié)構(gòu)有意識地去趨于一致。

有限內(nèi)存擬牛頓法(L—BFGS)主要思路是采用當(dāng)前迭代之前幾次的目標(biāo)函數(shù)梯度去逼近當(dāng)前迭代過程中的海森矩陣逆，L—BFGS逼近海森矩陣大致的流程如下所示：

1)選擇好雙參數(shù)初始模型m0，給出反演過程中最大的迭代次數(shù)IterMax，并給出反演迭代最小擬合差Rms；

3)設(shè)定好初始步長，在初始步長的基礎(chǔ)上進(jìn)行一維線搜索，尋求相對最佳的步長α：

mk+1=mk+αpk

(14)

4)計(jì)算更新后模型的擬合差是否達(dá)到給定的精度，若達(dá)到精度則終止反演，若迭代次數(shù)達(dá)到上限也終止反演，否則繼續(xù)回到第二步進(jìn)行循環(huán)。

4.2 理論合成數(shù)據(jù)反演算例

在進(jìn)行激發(fā)極化法超前探測前需明確該方法的有效探測范圍，實(shí)際工作中，觀測數(shù)據(jù)信噪比大于3時(shí)方可視為有效數(shù)據(jù)，激發(fā)極化法探測中最大誤差水平為5 %，即超過15 %的觀測異常數(shù)據(jù)方可視為有效數(shù)據(jù)，以視電阻率數(shù)據(jù)為例，視電阻率的異常幅度是指含水體存在時(shí)的視電阻率數(shù)據(jù)極小值相對于均勻圍巖情況下的視電阻率數(shù)據(jù)的變化程度，表示式如下：

(15)

均勻全空間圍巖電阻率默認(rèn)為1 000 Ω·m，隧道腔體電阻率默認(rèn)為1.0×106Ω·m，假定隧道寬12m，高12m；在均勻空間中植入低阻高極化異常體，參數(shù)如下：電阻率為10 Ω·m，極化率為20 %，長8 m，寬8 m，高6 m，來模擬隧道中富水溶洞構(gòu)造，如圖5所示。

數(shù)據(jù)采集采用定源三極裝置，點(diǎn)電源A位于掌子面處，沿隧道左、右邊墻、隧底分別布置測線，接收極距MN=3.0 m，收發(fā)距AO=3～117 m，通過正演模擬計(jì)算觀測數(shù)據(jù)，其視電阻率及視極化率響應(yīng)如圖6所示。

圖6 純隧道腔體正演響應(yīng)對比曲線Fig.6 Contrast curve of forward response of pure tunnel cavity

對視電阻率和視極化率數(shù)據(jù)分別附加5 %誤差，初始模型設(shè)定為電阻率1.0×106Ω·m均勻全空間，分別進(jìn)行雙參數(shù)同步反演和交叉梯度約束反演計(jì)算，反演結(jié)果取YZ剖面(X=12 m)及XZ剖面(Y=6 m)，如圖所示。

圖7 XZ剖面(Y=6 m)模型成果Fig.7 XZ profile (Y=6 m) model result drawing

圖8 YZ剖面(X=12 m)模型成果Fig.8 YZ profile (X=12 m) model result drawing

通過電阻率和極化率反演成果圖可得出：

1)同步反演電阻率和極化率時(shí)，電阻率反演結(jié)果受測線布置影響較大，深度信息缺失嚴(yán)重；極化率反演受測線布置影響較小，能大致圈定異常體形態(tài)。

2)引入交叉梯度約束反演電阻率和極化率時(shí)，可有效消除測線對電阻率反演的影響，電阻率和極化率兩種物性參數(shù)相互約束，均能大致圈定異常體形態(tài)。

3)隧道環(huán)境下觀測數(shù)據(jù)信息量存在局限性，測線主要包含X方向(隧道開挖方向)深度信息，對Y方向和Y方向深度信息缺失較大，導(dǎo)致X方向(隧道開挖方向)異常體分辨率較高，Y方向和Z方向分辨率較低。

5 結(jié) 論

本文以實(shí)現(xiàn)隧道環(huán)境下激發(fā)極化法三維正演和反演算法為核心任務(wù)，取得如下主要結(jié)論：

1)通過正演對隧道腔體影響進(jìn)行分析，得出視電阻率曲線受隧道腔體影響嚴(yán)重，視極化率曲線不受隧道腔體影響。

2)推導(dǎo)電阻率和極化率對應(yīng)的目標(biāo)函數(shù)梯度公式，依據(jù)有限內(nèi)存擬牛頓法實(shí)現(xiàn)電阻率和極化率同步反演，相對于傳統(tǒng)電阻率和極化率分布反演，提高了反演計(jì)算效率。

3)基于不同物性參數(shù)結(jié)構(gòu)耦合性，反演過程中引入交叉梯度函數(shù)約束電阻率和極化率，可有效提高反演精度。

4)實(shí)際探測中地質(zhì)條件更為復(fù)雜，數(shù)據(jù)干擾更為嚴(yán)重，沿隧道腔體布設(shè)測線得到的數(shù)據(jù)局限性較大，可優(yōu)化布極方式，也可通過綜合超前地質(zhì)預(yù)報(bào)獲得更多的先驗(yàn)信息，提高超前探測精度。