隧道瞬變電磁克希霍夫偏移成像與地質(zhì)災(zāi)害探測(cè)

戚志鵬，郭建磊,2，孫乃泉，李 貅，樊亞楠

(1.長(zhǎng)安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安 710054；2.中煤科工集團(tuán)西安研究院有限公司，陜西 西安 710077)

隧道與巷道超前預(yù)報(bào)不僅要預(yù)報(bào)災(zāi)害體的存在，更應(yīng)該確定災(zāi)害體的幾何結(jié)構(gòu)，判斷災(zāi)害的成災(zāi)性質(zhì)，為后續(xù)的災(zāi)害治理提供有效的信息[1]。瞬變電磁法作為一種電磁感應(yīng)方法，在低阻災(zāi)害體隧道超前預(yù)報(bào)中起到了預(yù)警的作用，避免了重大災(zāi)害事故的發(fā)生。但是隨著綠色工程施工的提出，在地下工程建設(shè)過程中不僅需要探測(cè)引起突水、突泥等低阻災(zāi)害體的物性信息，更需要獲得地質(zhì)災(zāi)害源的詳細(xì)結(jié)構(gòu)，為災(zāi)害治理提供勘探信息[2]，降低施工與災(zāi)害治理對(duì)環(huán)境的破壞。傳統(tǒng)的隧道瞬變電磁解釋方法僅給出了地下介質(zhì)的電性信息，對(duì)災(zāi)害體的幾何結(jié)構(gòu)特征描述并不明顯，因此，需要引入成像技術(shù)，提高瞬變電磁隧道超前預(yù)報(bào)方法的解釋精度[3]。

由于瞬變電磁具有對(duì)低阻敏感的特征，李貅等[4]將該方法引入到隧道超前地質(zhì)預(yù)報(bào)與檢測(cè)領(lǐng)域；薛國(guó)強(qiáng)等[5-6]利用適合全空間模型的等效導(dǎo)電平面法進(jìn)行瞬變電磁隧道超前預(yù)報(bào)解釋，使該方法進(jìn)一步推廣；Sun Huaifeng 等[7]利用視電阻率定義方法實(shí)現(xiàn)陣列式瞬變電磁隧道超前預(yù)報(bào)解釋；戚志鵬等[8]基于反函數(shù)原理實(shí)現(xiàn)全空間隧道瞬變電磁視電阻率解釋。傳統(tǒng)的瞬變電磁隧道超前預(yù)報(bào)方法均從擴(kuò)散方程出發(fā)，根據(jù)場(chǎng)與電阻率的轉(zhuǎn)換關(guān)系，實(shí)現(xiàn)地下介質(zhì)的電阻率成像。擴(kuò)散方程表征的是場(chǎng)的感應(yīng)與擴(kuò)散特征，主要描述物質(zhì)性質(zhì)的變化，對(duì)介質(zhì)結(jié)構(gòu)反應(yīng)不敏感[9]。而地震勘探或探地雷達(dá)等滿足的波動(dòng)方程主要描述波在介質(zhì)中的傳播特征，其中波的反射與折射特征對(duì)介質(zhì)幾何結(jié)構(gòu)反應(yīng)十分敏感。因此，在波動(dòng)方程條件下實(shí)現(xiàn)瞬變電磁解釋是提高隧道超前預(yù)報(bào)解釋精度的有效手段。

波場(chǎng)變換為瞬變電磁擬波動(dòng)場(chǎng)成像提供了前提基礎(chǔ)。瞬變電磁擴(kuò)散場(chǎng)方程與波動(dòng)方程存在數(shù)學(xué)變換關(guān)系[10-12]，將滿足擴(kuò)散方程的瞬變電磁場(chǎng)轉(zhuǎn)換為滿足波動(dòng)方程的虛擬波場(chǎng)后，可以借助地震勘探中發(fā)展起來的一些成像方法，解決復(fù)雜情況下瞬變電磁場(chǎng)的成像問題[13-14]；K.H.Lee 等[10]建立了滿足時(shí)域擴(kuò)散方程的電磁場(chǎng)與虛擬波場(chǎng)的關(guān)系式，通過這個(gè)變換式，可以將對(duì)應(yīng)地電模型的波動(dòng)方程模擬結(jié)果變換成時(shí)域電磁響應(yīng)，但該變換的意義并不在于對(duì)電磁場(chǎng)的正演模擬，而在于其逆變換過程；李貅等[15]提出了從瞬變電磁場(chǎng)到波場(chǎng)的優(yōu)化算法。在波場(chǎng)的正變換中，在保證瞬變電磁場(chǎng)的計(jì)算精度的同時(shí)，應(yīng)用兩步最優(yōu)化和正則化算法，成功地解決了波場(chǎng)反變換問題。戚志鵬等[16]等利用預(yù)條件正則化實(shí)現(xiàn)了波場(chǎng)的全時(shí)域變換算法。上述研究表明，變換得到的虛擬波場(chǎng)滿足波動(dòng)方程，具有傳播、反射、透射特征，這為瞬變電磁場(chǎng)的偏移成像創(chuàng)造了條件。

瞬變電磁虛擬波場(chǎng)成像，可以借鑒地震勘探處理中比較成熟的偏移成像技術(shù)。1978 年W.A.Schneider提出Kirchhoff 積分法偏移，并得到了廣泛的應(yīng)用[17]。Kirchhoff 積分法偏移建立在波動(dòng)方程積分解的基礎(chǔ)上，它利用Kirchhoff 繞射積分公式把分散在地表的各個(gè)測(cè)點(diǎn)上來自同一繞射點(diǎn)的能量匯集在一起，置于地下相應(yīng)的物理繞射點(diǎn)上，從而實(shí)現(xiàn)繞射點(diǎn)偏移成像。李貅等[18]將Kirchhoff 積分成像引入到地面瞬變電磁三維數(shù)據(jù)虛擬波場(chǎng)成像之中，得到了較好的應(yīng)用效果。

筆者在實(shí)現(xiàn)瞬變電磁波場(chǎng)變換基礎(chǔ)上，利用克?；舴蚍e分法對(duì)瞬變電磁虛擬波場(chǎng)進(jìn)行延拓，實(shí)現(xiàn)瞬變電磁場(chǎng)的波場(chǎng)成像。不同于地面瞬變電磁解釋方法，筆者將該方法應(yīng)用于隧道超前預(yù)報(bào)二維數(shù)據(jù)偏移之中。在算法中，首先實(shí)現(xiàn)擴(kuò)散場(chǎng)到虛擬波場(chǎng)的變換；然后忽略走向影響，基于二維克?；舴蚍e分實(shí)現(xiàn)隧道偏移成像，分別對(duì)理論模型與實(shí)測(cè)資料進(jìn)行處理驗(yàn)證文中算法。

1 基本原理

擴(kuò)散方程主要反映的是低頻電磁場(chǎng)傳播過程中的擴(kuò)散和感應(yīng)特征，直接利用瞬變電磁擴(kuò)散場(chǎng)數(shù)據(jù)進(jìn)行延拓成像，結(jié)果對(duì)電性界面分辨能力差，體積效應(yīng)強(qiáng)。而波動(dòng)方程反映的是波的傳播特征，基于波動(dòng)方程的處理方法能準(zhǔn)確地確定介質(zhì)的物性界面，給出異常體的幾何信息。由于波動(dòng)方程和擴(kuò)散方程之間的差異，將滿足擴(kuò)散方程的瞬變電磁數(shù)據(jù)轉(zhuǎn)換為滿足波動(dòng)方程的波場(chǎng)數(shù)據(jù)，借鑒地震勘探領(lǐng)域成熟技術(shù)進(jìn)行電磁數(shù)據(jù)解釋，是提高隧道超前預(yù)報(bào)質(zhì)量的有效途徑。

1.1 瞬變電磁波場(chǎng)變換原理

將瞬變電磁擴(kuò)散場(chǎng)轉(zhuǎn)化成虛擬波場(chǎng)的前提條件是瞬變電磁場(chǎng)與波場(chǎng)的各個(gè)分量具有數(shù)學(xué)上的對(duì)應(yīng)關(guān)系?，F(xiàn)以磁場(chǎng)的垂直分量為例，它們滿足如下關(guān)系：

引入平滑約束，構(gòu)建目標(biāo)函數(shù)φ如下所示：

式中：‖‖2為向量的二范數(shù)；α為正則化因子；D為2 階差分矩陣。

通過目標(biāo)函數(shù)的最小化，即可獲得U的近似估計(jì)，其等價(jià)方程如下式所示。

鐘華森[19]、魯凱亮[20]等給出了式(4)的波場(chǎng)變換算法，即設(shè)定一個(gè)時(shí)間窗口，使之在整個(gè)時(shí)間序列上滑動(dòng)，每個(gè)窗口內(nèi)矩陣階數(shù)較小，可利用奇異值分解法或精細(xì)積分法進(jìn)行波場(chǎng)變換獲得波場(chǎng)值；在全時(shí)間段采用正則化迭代法進(jìn)行波場(chǎng)變換；判斷各窗口變換所得波場(chǎng)值與全時(shí)段波場(chǎng)值之間的相關(guān)性，若相關(guān)系數(shù)大于一定閾值，則將小時(shí)窗內(nèi)波場(chǎng)值與全時(shí)段的變換結(jié)果進(jìn)行疊加，最終得到高分辨的波場(chǎng)[20]。

1.2 Kirchhoff 偏移成像原理

將擴(kuò)散場(chǎng)變換到波場(chǎng)后，可以利用地震波場(chǎng)的資料處理方法進(jìn)行數(shù)據(jù)處理。通過對(duì)波場(chǎng)數(shù)據(jù)的處理分析，即可獲得地下介質(zhì)分界面的特征，進(jìn)而借以研究地下地質(zhì)巖層結(jié)構(gòu)及其物性特征。地下介質(zhì)分界面的求取問題，可以看作一種反散射問題，而地震偏移技術(shù)就是求解反散射問題、獲取反射界面分布的一種技術(shù)?？讼；舴蚍e分法是在地震資料偏移處理中應(yīng)用較為廣泛的一種方法，本文利用此方法進(jìn)行波場(chǎng)數(shù)據(jù)的偏移成像。

在直角坐標(biāo)系中，地下無源區(qū)波場(chǎng)所滿足的標(biāo)量波動(dòng)方程為：

當(dāng)對(duì)邊界(圖1)進(jìn)行偏移時(shí)，可忽略外邊界的積分，且地表與xoy平面重合時(shí)，克希霍夫解為：

圖1 積分區(qū)域與邊界Fig.1 Schematic diagram of regions and boundaries for integration

式中：[U]為波場(chǎng)的延遲值，A/m；r為計(jì)算點(diǎn)到積分掃描點(diǎn)的距離，m；v為虛擬波速度，s為面積，m2；n為邊界的法向方向。

若已知xoy平面上波場(chǎng)值和波場(chǎng)導(dǎo)數(shù)值，即可利用此式計(jì)算空間中任一點(diǎn)處任意時(shí)刻的波場(chǎng)值。

在進(jìn)行偏移時(shí)，測(cè)量數(shù)據(jù)定義在z=0 平面上。對(duì)于自激自收情況，二次波場(chǎng)產(chǎn)生時(shí)刻為地面接收到反射波時(shí)刻的一半，若取偏移速度為真實(shí)速度的一半，則二次場(chǎng)產(chǎn)生的時(shí)刻對(duì)應(yīng)為t=0。將式(6)中的延遲式替代為超前值，即令：

式中：下標(biāo)p表示偏移點(diǎn)p的坐標(biāo)。

1.3 隧道瞬變電磁偏移成像技術(shù)

與地面瞬變電磁三維偏移不同，在進(jìn)行隧道或巷道的超前預(yù)報(bào)時(shí)，經(jīng)常處理來自一條測(cè)線上的數(shù)據(jù)。傳統(tǒng)擴(kuò)散場(chǎng)瞬變電磁不能直接對(duì)y坐標(biāo)進(jìn)行忽略，而是將y坐標(biāo)轉(zhuǎn)化成波數(shù)域進(jìn)行求解；但是在波場(chǎng)中可以直接將y坐標(biāo)忽略，從而簡(jiǎn)化數(shù)據(jù)解釋。此外，克?；舴蚍e分的邊界與區(qū)域如圖2 所示，將隧道邊界視為無窮遠(yuǎn)邊界，并忽略無窮遠(yuǎn)邊界的影響，則有二維偏移場(chǎng)如下式所示。

圖2 隧道偏移區(qū)域與邊界Fig.2 Schematic diagram of the tunnel area and boundary

在計(jì)算隧道工作面前方虛擬波場(chǎng)時(shí)不僅要已知地表的波場(chǎng)，還需已知在地下介質(zhì)中虛擬波場(chǎng)的傳播速度。根據(jù)波場(chǎng)變換關(guān)系可知波場(chǎng)的傳播速度與地下介質(zhì)的電阻率分布有關(guān)，如下式所示。

式中：μ為真空磁導(dǎo)率；H/m，σ為空間電導(dǎo)率，S/m。

當(dāng)μ=μ0時(shí)，虛擬波場(chǎng)的速度只與大地介質(zhì)的電導(dǎo)率有關(guān)，能否準(zhǔn)確地獲得地下介質(zhì)的電導(dǎo)率，成為速度分析準(zhǔn)確與否的關(guān)鍵。電阻率(電導(dǎo)率的倒數(shù))是瞬變電磁解釋方法的重要參數(shù)，可以根據(jù)隧道瞬變電磁電阻率成像進(jìn)行虛擬波場(chǎng)速度提取。當(dāng)獲得介質(zhì)電性分布后可根據(jù)具體情況采用不同的速度分析方法。原則上只要能求得準(zhǔn)確的地層電導(dǎo)率，便能實(shí)現(xiàn)精確的速度分析，當(dāng)前隧道瞬變電磁電阻率計(jì)算可以根據(jù)反函數(shù)原理進(jìn)行計(jì)算[8]。

2 理論模型分析

在隧道建設(shè)過程中突水突泥災(zāi)害主要由斷層、溶洞等水源災(zāi)害導(dǎo)致，因此，在超前預(yù)報(bào)中實(shí)現(xiàn)斷層與溶洞的預(yù)報(bào)十分重要。為了驗(yàn)證文中方法，分別設(shè)計(jì)隧道低阻夾層與含水溶洞模型，采用瞬變電磁虛擬波場(chǎng)克希霍夫偏移算法對(duì)災(zāi)害體進(jìn)行預(yù)報(bào)。

2.1 隧道充水?dāng)鄬幽Ｐ?/h3>

設(shè)計(jì)充水?dāng)鄬幽Ｐ?，如圖3 所示。其中隧道寬度為12 m，背景電阻率ρ1=100 Ω·m；傾斜斷層異常電阻率ρ2=1 Ω·m，寬5 m；中心點(diǎn)距工作面埋深為18 m，與水平線的傾斜角為75°。采用中心回線裝置進(jìn)行探測(cè)，發(fā)射線框3 m×3 m，發(fā)射電流1 A，點(diǎn)距0.5 m。

圖3 充水?dāng)鄬幽Ｐ虵ig.3 Synthetic model of the water-filling fault

利用三維正演進(jìn)行數(shù)值模擬，采用波場(chǎng)變換將衰減電壓轉(zhuǎn)換為虛擬波場(chǎng)，采用全時(shí)域視電阻率計(jì)算方法計(jì)算工作面前介質(zhì)電性分布，如圖4a 所示。采用本文提出的2D 克希霍夫變換方法進(jìn)行偏移成像，成像結(jié)果如圖4b 所示，圖中同相軸與模型邊界吻合且能夠顯示異常體的傾斜方向；視電阻率結(jié)果能夠顯示出低阻區(qū)域與異常傾向。將電阻率圖與偏移成像結(jié)果組合，如圖4c 所示，能夠有效顯示低阻異常邊界與區(qū)域電性特征，為災(zāi)害預(yù)報(bào)與治理提供有效的幾何邊界信息。

圖4 隧道充水?dāng)鄬游闹蟹椒ǖ木C合解釋Fig.4 Interpretation of tunnel water-filled fault

2.2 充水溶洞模型

隧道施工中常會(huì)遇到充水溶洞等情況，文中設(shè)計(jì)了如圖5 所示的充水溶洞隧道模型，隧道寬12 m，背景電阻率為1 000 Ω·m；含水溶洞尺寸6 m ×6 m×5 m、中心點(diǎn)距工作面35 m、電阻率50 Ω·m。采用中心回線裝置進(jìn)行探測(cè)，裝置參數(shù)與充水?dāng)鄬幽Ｐ鸵恢隆?/p>

圖5 充水溶洞模型Fig.5 Synthetic model of water-filling karst cave

模型中溶洞截面為矩形而電阻率斷面圖中則顯示為類似的圓形低阻區(qū)域(圖6a)，偏移成像結(jié)果(圖6b)準(zhǔn)確顯示模型的邊界位置。將電阻率與偏移成像二者結(jié)合(圖6c)，既能解釋異常的邊界結(jié)構(gòu)又能分辨電性分布特征。

圖6 充水溶洞模型本文方法的綜合解釋Fig.6 Interpretation of water-filled model

3 實(shí)測(cè)數(shù)據(jù)分析

以西南地區(qū)實(shí)際隧道工程為例，利用瞬變電磁進(jìn)行隧道超前預(yù)報(bào)，采用虛擬波場(chǎng)克希霍夫偏移成像方法進(jìn)行數(shù)據(jù)處理。測(cè)線位于隧道橫洞內(nèi)靠近出水孔一側(cè)，采用中心回線裝置進(jìn)行數(shù)據(jù)采集，發(fā)射線圈邊長(zhǎng)3 m，發(fā)射頻率8.3 Hz，發(fā)射電流5 A。橫洞前方10 m左右有一已知出水孔且地質(zhì)資料已知。為了驗(yàn)證方法有效性，在橫洞側(cè)壁采用中心回線裝置進(jìn)行測(cè)量如圖7所示，引水隧洞無襯砌，且遠(yuǎn)離工作面，電磁干擾影響較小，實(shí)測(cè)數(shù)據(jù)如圖8 所示。

圖7 橫洞側(cè)壁施工布置Fig.7 Operation layout in auxiliary tunnel sidewall

圖8 測(cè)點(diǎn)瞬變電磁原如曲線Fig.8 Original transient electromagnetic curve of measuring point

采用文中所述的虛擬波場(chǎng)克?；舴蚱品椒ㄟM(jìn)行數(shù)據(jù)處理，處理結(jié)果如圖9 所示。由圖9 可知，在深度為10 m 位置有一低電阻率薄層，根據(jù)同相軸圖可知層厚較小推斷為涌水裂隙。這一結(jié)果與已知出水點(diǎn)位置相符。在深度為10～20 m 處，電阻率介于巖石與水體之間，推斷為裂隙周邊受水體浸泡的巖層；在測(cè)點(diǎn)4?測(cè)點(diǎn)10、深度大于20 m 后電阻率整體較大，推斷為完整的巖體。成像解釋結(jié)果與實(shí)際情況一致，說明方法是穩(wěn)定的，能夠應(yīng)用于野外數(shù)據(jù)處理。

圖9 實(shí)測(cè)數(shù)據(jù)及其克希霍夫偏移成像解釋結(jié)果Fig.9 Kirchhoff migration imaging result of survey data

4 結(jié) 論

a.瞬變電磁場(chǎng)經(jīng)積分逆變換后可轉(zhuǎn)換為滿足波動(dòng)方程的虛擬波場(chǎng)。利用克?；舴蚍e分方法可將虛擬波場(chǎng)從地面向地下進(jìn)行逆時(shí)偏移，得到清晰的異常邊界。將得到電阻率的斷面分布與偏移成像結(jié)果進(jìn)行綜合實(shí)現(xiàn)聯(lián)合反演解釋。

b.在擴(kuò)散場(chǎng)瞬變電磁解釋中，對(duì)于二維地電斷面情況需要將走向方向轉(zhuǎn)化為波數(shù)域，增加了數(shù)據(jù)解釋的難度，而轉(zhuǎn)化成虛擬波場(chǎng)后二維解釋可將走線方向坐標(biāo)直接忽略，不僅降低了二維數(shù)據(jù)解釋難度，而且獲得了準(zhǔn)確的剖面構(gòu)造特征。

c.將虛擬波場(chǎng)克?；舴蚨S偏移成像方法應(yīng)用于理論模型數(shù)據(jù)處理，綜合解釋結(jié)果與理論模型相符，不僅能夠判斷異常電阻率分布區(qū)域，同時(shí)能夠揭示災(zāi)害體的傾向、電性邊界等構(gòu)造特征，充分驗(yàn)證了方法的正確性。將方法應(yīng)用于實(shí)際資料處理，綜合解釋結(jié)果圈定了低阻異常位置，波場(chǎng)成像有效顯示導(dǎo)水通道的邊界與走向特征，解釋結(jié)構(gòu)與隧道開挖揭露結(jié)果一致。理論結(jié)果與實(shí)測(cè)數(shù)據(jù)處理結(jié)果表明瞬變電磁虛擬波場(chǎng)克?；舴蚱瞥上窨梢詫?shí)現(xiàn)隧道超前預(yù)報(bào)，能夠?yàn)榈刭|(zhì)災(zāi)害預(yù)警與防治提供豐富的災(zāi)害體地質(zhì)信息。