趙廣學(xué)，阮帥，吳肅元

(1.恒達(dá)新創(chuàng)(北京)地球物理技術(shù)有限責(zé)任公司，北京 100020； 2.中國(guó)地質(zhì)科學(xué)院，北京 100037)

0 引言

在鐵路隧道施工前期，利用音頻大地電磁測(cè)深(以下簡(jiǎn)稱(chēng)AMT)對(duì)施工線進(jìn)行斷層破碎帶、溶洞等地質(zhì)災(zāi)害隱患評(píng)估時(shí)，勘探數(shù)據(jù)的確定性和分辨率都非常重要，但二者往往不可兼得，提升反演分辨率意味著增加假異常的出現(xiàn)幾率，造成過(guò)度風(fēng)險(xiǎn)評(píng)估，反之則會(huì)增加施工風(fēng)險(xiǎn)，造成安全事故。因此，對(duì)鐵路隧道AMT勘探數(shù)據(jù)的反演進(jìn)行更深入的研究極其重要。

1 復(fù)雜地電模型及二維非線性共軛梯度反演

在石油勘探領(lǐng)域，研究多在地質(zhì)條件接近一維的情況下進(jìn)行，AMT一維反演縱向分辨率較高，往往通過(guò)淺層電阻率調(diào)查對(duì)觀測(cè)數(shù)據(jù)進(jìn)行靜態(tài)位移校正，然后再進(jìn)行一維反演直接解釋?zhuān)艿玫讲诲e(cuò)的地層分辨率。但這種反演方案只適應(yīng)于沉積巖地層，在復(fù)雜地形地質(zhì)條件下，靜位移校正并無(wú)簡(jiǎn)單規(guī)律可循，盲目地使用一維反演方案會(huì)出現(xiàn)電阻率失真和中、淺層假異常的問(wèn)題。隧道勘探的數(shù)據(jù)反演應(yīng)該使用二維算法，避免對(duì)輸入數(shù)據(jù)的靜態(tài)效應(yīng)校正，除非因觀測(cè)電極距測(cè)量誤差出現(xiàn)明顯的曲線平移現(xiàn)象。

目前，主流的AMT二維反演算法是Rodi和Mackie(2001)改進(jìn)的非線性共軛梯度法(以下簡(jiǎn)稱(chēng)NLCG)，它避免了直接求解雅克比矩陣，一次反演迭代只需3次正演，從而大大節(jié)省了計(jì)算量，具有迭代穩(wěn)定，內(nèi)存需求小的優(yōu)點(diǎn)，非常適應(yīng)二、三維MT/AMT數(shù)據(jù)反演。其目標(biāo)函數(shù)定義為[1]：

Ψ(m)=[d-F(m)]TV-1[d-F(m)]+

λmTLTLm,

(1)

式中：m為反演模型;F(m)為模型的二維正演響應(yīng);d為觀測(cè)數(shù)據(jù);λ為正則化因子，控制模型修正沿著光滑方向還是數(shù)據(jù)最佳擬合方向；V是正定矩陣，為誤差向量的方差；L為拉爾朗日算子矩陣。

若使用線性共軛梯度法求解式(1)的極小值，對(duì)于每次迭代，需進(jìn)行以下模型修正：

k=0,1,…,k-1。

式中：ml為第l次迭代的模型，pl,k控制了第l次迭代中的模型修正方向，αl,k為修正步長(zhǎng)因子。

自2001年以來(lái)，這一反演方案成為大地電磁(MT)與AMT數(shù)據(jù)反演研究和應(yīng)用最多的方法[3]。不過(guò)，在實(shí)際應(yīng)用中，若改變不同的反演參數(shù)(如初始模型、正則化因子)，同一組數(shù)據(jù)的反演結(jié)果存在較大的差別，不同的參數(shù)控制了模型的修正方向，側(cè)重點(diǎn)可以是放在反演結(jié)果的平滑度上，也可以是觀測(cè)數(shù)據(jù)的擬合度。因此對(duì)于具體的地質(zhì)任務(wù)，往往需要通過(guò)多種組合實(shí)驗(yàn)的方法獲取最佳選擇方案。

基于以上事實(shí)，為適應(yīng)鐵路隧道AMT勘探的解釋需求，按照西南某鐵路隧道線的地形、地質(zhì)情況建立地電斷面模型(圖1)。該地電模型斷面背景為類(lèi)K型地電斷面，背景地層電阻率自淺至深分別為20、100、500、800 Ω·m和100 Ω·m；在距離600～800 m、1 200～1 600 m處各存在一個(gè)低阻斷層破碎帶，電阻率10 Ω·m；距離2 000～2 200 m處存在一塊電阻率為1 000 Ω·m的高阻空腔；分別在200～400 m、1 000～2 000 m、2 000～2 200 m的地形起伏處設(shè)置電阻率為500、500、100 Ω·m的淺層不均勻體以模擬靜態(tài)效應(yīng)。

圖1 典型地電斷面模型的正演網(wǎng)格Fig.1 Typical forward grid of geological model

對(duì)此模型的合成AMT數(shù)據(jù)(頻帶10 000～10 Hz)進(jìn)行不同參數(shù)組合的二維NLCG反演[4]，以分析參數(shù)選取對(duì)反演結(jié)果的影響，從而建立鐵路隧道勘探二維NLCG反演的參數(shù)選取準(zhǔn)則，適應(yīng)地質(zhì)解釋需要。

2 反演網(wǎng)格的影響

圖1的模型高程起伏劇烈，海拔最高的Q測(cè)點(diǎn)和最低的E測(cè)點(diǎn)之間高程差接近180 m，實(shí)際勘探中還有更劇烈的高程變化。只有精細(xì)的網(wǎng)格剖分才能準(zhǔn)確模擬靜態(tài)效應(yīng)以擬合實(shí)測(cè)數(shù)據(jù)，但過(guò)密的反演網(wǎng)格會(huì)增加反演問(wèn)題的病態(tài)，成倍增加計(jì)算時(shí)間。因此討論反演網(wǎng)格的疏密對(duì)結(jié)果的影響非常必要。

選擇2種網(wǎng)格進(jìn)行對(duì)比研究。第一種是讓測(cè)點(diǎn)位于網(wǎng)格中心、橫向不做加密、縱向以1.2倍因子增長(zhǎng)的粗糙網(wǎng)格(簡(jiǎn)稱(chēng)粗糙網(wǎng)格)，本文設(shè)置的網(wǎng)格為64×16；另一種反演網(wǎng)格與合成數(shù)據(jù)的正演網(wǎng)格一致(簡(jiǎn)稱(chēng)精細(xì)網(wǎng)格)，文中設(shè)置的網(wǎng)格為64×64。為不失一般性，本文使用統(tǒng)一的迭代次數(shù)，在相同的迭代次數(shù)下進(jìn)行對(duì)比。初始模型為100 Ω·m均勻空間，最大迭代次數(shù)100次，正則化因子取3，合成數(shù)據(jù)的各種反演組合的結(jié)果如圖2～圖7所示，反演收斂RMS方差如表1所示。

比較圖2和圖3可知，只選用TE模式數(shù)據(jù)進(jìn)行NLCG反演時(shí)，網(wǎng)格的粗糙度對(duì)反演結(jié)果的影響不是太大，兩個(gè)反演結(jié)果的斷層破碎帶和高阻空腔的位置基本一致，地層電阻率也基本正確[4]。

和TE模式數(shù)據(jù)反演不同，參照?qǐng)D4和圖5可見(jiàn)，若只選取TM模式數(shù)據(jù)進(jìn)行反演，粗糙網(wǎng)格的反演結(jié)果對(duì)破碎帶的分辨能力下降，二者的電阻率值有較大差別。

圖2 粗糙網(wǎng)格TE模式反演結(jié)果(圖例同圖1)Fig.2 Rough mesh TE mode inversion results (The legend is the same as Figure 1)

圖3 精細(xì)網(wǎng)格TE模式反演結(jié)果(圖例同圖1)Fig.3 Fine mesh TE mode inversion results (The legend is the same as Figure 1)

圖4 粗糙網(wǎng)格TM模式反演結(jié)果(圖例同圖1)Fig.4 Rough mesh TM mode inversion results (The legend is the same as Figure 1)

表1 不同網(wǎng)格、不同模式的反演RMS誤差

若同時(shí)對(duì)TE和TM模式的數(shù)據(jù)進(jìn)行反演(參考圖6和圖7)，二者的橫向分辨率差不多，不同之處在于粗糙網(wǎng)格的剖面整體電阻率相對(duì)較高，對(duì)于鐵路勘探而言這不是重要問(wèn)題，對(duì)高阻空腔的反映能力還是精細(xì)網(wǎng)格更好。

若考慮反演的RMS誤差，參考表1，所有粗糙網(wǎng)格的最終RMS誤差都在2以上，而精細(xì)網(wǎng)格的RMS誤差都小于1。精細(xì)網(wǎng)格不但有利于模擬靜態(tài)位移、實(shí)現(xiàn)更準(zhǔn)確的異常分辨，更有利于達(dá)到更小的RMS誤差。

圖5 精細(xì)網(wǎng)格TM模式反演結(jié)果(圖例同圖1)Fig.5 Fine mesh TM mode inversion results (The legend is the same as Figure 1)

圖6 粗糙網(wǎng)格TE+TM模式聯(lián)合反演結(jié)果(圖例同圖1)Fig.6 Rough mesh TE+TM mode inversion results (The legend is the same as Figure 1)

圖7 精細(xì)網(wǎng)格TE+TM模式聯(lián)合反演結(jié)果(圖例同圖1)Fig.7 Fine mesh TM+TE mode inversion results (The legend is the same as Figure 1)

綜上所述，實(shí)際資料反演中，若只使用TE模式反演，網(wǎng)格可相對(duì)粗糙；若TM模式數(shù)據(jù)參與反演，則要求地表網(wǎng)格相對(duì)精細(xì)，以盡量精確地模擬地形變化。圖6和圖7是最佳的鐵路隧道勘探反演方案，即使用相對(duì)精細(xì)網(wǎng)格的TE+TM模式聯(lián)合反演。

3 不同頻帶的反演實(shí)驗(yàn)

由于AMT的高頻部分主要場(chǎng)源為高空雷電，一般在晴朗的上午，即使工區(qū)干擾較小，也會(huì)在5 000～1 000 Hz頻帶出現(xiàn)較差的信噪比，造成觀測(cè)曲線脫節(jié)。這個(gè)頻帶也被稱(chēng)為“死帶(death band)”。由于“死帶”影響，野外數(shù)據(jù)高于1 000 Hz的功率譜取舍很難，有些學(xué)者使用插值計(jì)算出“死帶”數(shù)據(jù)，不過(guò)這并非真實(shí)的二維響應(yīng)。

為評(píng)估“死帶”對(duì)深部異常體的二維NLCG反演影響，對(duì)圖1的合成數(shù)據(jù)使用1 000～10 Hz的中、低頻數(shù)據(jù)進(jìn)行精細(xì)網(wǎng)格反演，其結(jié)果如圖8～10所示，其他控制參數(shù)和本文第2節(jié)一致，反演的最終RMS誤差如表2所示。

比較圖3和圖8可見(jiàn)，對(duì)于TE模式反演，舍棄高頻數(shù)據(jù)(10 000～1 000 Hz)之后，反演結(jié)果基本相同，不同的是圖3由于有更多的高頻數(shù)據(jù)參與反演，對(duì)1 200～1 400 m處的低阻破碎帶分辨更好。

比較圖5和圖9可見(jiàn)，TM模式數(shù)據(jù)的二維NLCG反演對(duì)高頻數(shù)據(jù)更依賴(lài)，兩個(gè)斷層破碎帶的形狀和深度都發(fā)生了較大的變化，這是因?yàn)門(mén)M數(shù)據(jù)的高頻部分更易受到地形和地表不均勻影響，若不擬合高頻部分則模型等效性更強(qiáng)。

對(duì)比圖7和圖10，若使用TE和TM數(shù)據(jù)聯(lián)合反演，無(wú)論是否保留高頻數(shù)據(jù)，深部模型的結(jié)果差別不大，對(duì)3個(gè)橫向異常都有明顯反映，且地層電阻率無(wú)明顯差別。

從反演的RMS誤差看，舍棄高頻部分更有利于數(shù)據(jù)擬合，采取1 000～10 Hz的數(shù)據(jù)進(jìn)行反演的最終RMS誤差均小于0.5，因此若實(shí)測(cè)資料有比較嚴(yán)重的“死帶”現(xiàn)象，只要使用TE+TM聯(lián)合模式的精細(xì)網(wǎng)格反演，不妨將10 000～1 000 Hz的高頻數(shù)據(jù)舍去。該對(duì)比由于頻點(diǎn)數(shù)不一樣，頻點(diǎn)數(shù)少的擬合更容易，誤差會(huì)更小，因此這個(gè)對(duì)比主要是對(duì)比反演結(jié)果的精細(xì)度。綜上所述，實(shí)際資料的二維NLCG反演若單獨(dú)使用TM數(shù)據(jù)，最好盡可能多地使用高頻數(shù)據(jù)，以對(duì)靜態(tài)位移現(xiàn)象進(jìn)行最佳擬合，若反演輸入數(shù)據(jù)是TE模式或TE+TM模式聯(lián)合，可以考慮只使用1 000 Hz以下的數(shù)據(jù)進(jìn)行反演。

表2 不同頻帶的反演RMS誤差

圖8 1 000～10 Hz 精細(xì)網(wǎng)格TE模式反演結(jié)果(圖例同圖1)Fig.8 1 000～10 Hz Fine mesh TE mode inversion results(The legend is the same as Figure 1)

圖9 1 000～10 Hz 精細(xì)網(wǎng)格TM模式反演結(jié)果(圖例同圖1)Fig.9 1 000～10 Hz Fine mesh TM mode inversion results(The legend is the same as Figure 1)

圖10 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式聯(lián)合反演結(jié)果(圖例同圖1)Fig.10 1 000～10 Hz Fine mesh TE+TM mode inversion results (The legend is the same as Figure 1)

4 不同正則化因子的反演實(shí)驗(yàn)

正則化因子實(shí)際上可視為一個(gè)加權(quán)系數(shù)，將數(shù)據(jù)擬合的目標(biāo)函數(shù)和模型約束(粗糙度)目標(biāo)函數(shù)聯(lián)合起來(lái)，構(gòu)成總目標(biāo)函數(shù)[5]。在NLCG二維反演程序中，該值是一個(gè)經(jīng)驗(yàn)常數(shù)，它既能讓模型的修正在一定程度上擬合數(shù)據(jù)，挖掘數(shù)據(jù)的分辨率，又要求模型具有較好光滑度，以減少冗余構(gòu)造。一般來(lái)說(shuō)，正則化因子越大，模型越光滑，實(shí)際資料的最佳正則化因子往往采用實(shí)驗(yàn)的方式獲得，文中分別選用0.3、3、30、300的正則化因子對(duì)1 000～10 Hz的合成數(shù)據(jù)進(jìn)行TE+TM模式聯(lián)合反演，其他反演參數(shù)同前，反演結(jié)果如圖11～14所示，最終RMS誤差見(jiàn)表3。

表3 不同正則化因子的反演RMS誤差

由反演結(jié)果(圖11～13)可見(jiàn)，在其他參數(shù)不變的情況下，隨著正則化因子的增大，模型的分辨率越來(lái)越差，接近30～300的正則化因子的二維NLCG反演結(jié)果基本對(duì)深部斷層破碎帶和高阻空腔完全沒(méi)有反映。

從表3可知，隨著正則化因子的增大，為了讓模型修正向更光滑的方向，犧牲了數(shù)據(jù)的擬合能力，當(dāng)正則化因子達(dá)到300時(shí)，最終RMS誤差達(dá)到 3.46。因此，對(duì)于本文討論的鐵路隧道的AMT數(shù)據(jù)二維NLCG反演應(yīng)使用較小的正則化因子，一般可以在0.1～10之間進(jìn)行反演實(shí)驗(yàn)，選取一個(gè)地層清晰、RMS誤差又不太大的即可。

圖11 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式反演結(jié)果(正則化因子0.3)Fig.11 1 000～10 Hz Fine mesh TE+TM mode inversion results (smooth factor 0.3)

圖12 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式反演結(jié)果(正則化因子3)Fig.12 1 000～10 Hz Fine mesh TE+TM mode inversion results (smooth factor 3)

圖13 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式反演結(jié)果(正則化因子30)Fig.13 1 000～10 Hz Fine mesh TE+TM mode inversion results (smooth factor 30)

圖14 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式反演結(jié)果(正則化因子300)Fig.14 1 000～10 Hz Fine mesh TE+TM mode inversion results (smooth factor 300)

5 不同初始模型的反演實(shí)驗(yàn)

二維反演的初始模型會(huì)對(duì)最終結(jié)果產(chǎn)生影響，很多學(xué)者認(rèn)為應(yīng)該使用一維反演結(jié)果作為二維反演的初始模型[6]。但這一觀點(diǎn)必須考慮二維NLCG算法的特點(diǎn)，若初始模型太粗糙(比如靜態(tài)效應(yīng)嚴(yán)重、或?qū)崪y(cè)曲線不夠圓滑的一維反演結(jié)果)，則很有可能導(dǎo)致誤差不收斂。

二維NLCG的初始模型最好比較光滑，可以考慮使用較大正則化因子的二維反演結(jié)果，但使用均勻半空間更有利于模型研究。選取10、100、500、1 000 Ω·m 4種均勻半空間作為初始模型，采用1 000～10 Hz的TE+TM聯(lián)合數(shù)據(jù)，正則化因子取0.3，其他控制參數(shù)同前，反演結(jié)果如圖15～18所示，最終RMS結(jié)果如表4所示。

由反演結(jié)果可知，所有模型的反演結(jié)果對(duì)于1 200～1 400 m處的斷層破碎帶基本正確，這證明二維NLCG反演對(duì)初始模型的依賴(lài)性相對(duì)較弱。不過(guò)，當(dāng)初始模型增加到500 Ω·m后，600～800 m處的斷層破碎帶已經(jīng)出現(xiàn)了錯(cuò)誤的傾向，當(dāng)初始模型變?yōu)? 000 Ω·m時(shí)，傾向錯(cuò)誤更明顯，而且反演結(jié)果的整體電阻率變得過(guò)大[7]。

圖15 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式反演結(jié)果(初始模型10 Ω·m)Fig.15 1 000～10 Hz Fine mesh TE+TM mode inversion results (initial model 10 Ω·m)

圖16 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式反演結(jié)果(初始模型100 Ω·m)Fig.16 1 000～10 Hz Fine mesh TE+TM mode inversion results (initial model 100 Ω·m)

圖17 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式反演結(jié)果(初始模型500 Ω·m)Fig.17 1 000～10 Hz Fine mesh TE+TM mode inversion results (initial model 500 Ω·m)

圖18 1 000～10 Hz 精細(xì)網(wǎng)格TE+TM模式反演結(jié)果(初始模型1 000 Ω·m)Fig.18 1 000～10 Hz Fine mesh TE+TM mode inversion results (initial model 1 000 Ω·m)

表4 不同初始模型的反演RMS誤差

參考表4給出的反演RMS誤差可知，當(dāng)初始模型的電阻率取得和中淺部地層電阻率接近時(shí)，RMS誤差最小。因此，實(shí)際資料反演最好首先嘗試選取和中、淺層電阻率較接近的初始模型，這一電阻率可用較大的正則化因子進(jìn)行二維反演得到。

6 結(jié)論

進(jìn)行了不同參數(shù)組合的典型鐵路隧道AMT勘探的模型非線性共軛梯度反演試驗(yàn)，以確定實(shí)際資料反演的最佳參數(shù)選取方式。結(jié)論如下：

1) TM模式數(shù)據(jù)的反演對(duì)網(wǎng)格剖分要求更高。實(shí)際資料若只使用TE模式反演，可使用粗糙網(wǎng)格，但模型分析的結(jié)果建議鐵路隧道勘探的AMT數(shù)據(jù)最好使用TE+TM模式聯(lián)合反演，反演網(wǎng)格的淺層應(yīng)盡量精細(xì)以精確模擬地形。

2) 類(lèi)似的，TM模式對(duì)于參與反演的模型頻帶更敏感，實(shí)際數(shù)據(jù)若高頻帶較差需要舍去1 000 Hz以上的數(shù)據(jù)時(shí)，最好使用TE模式或TE+TM模式的聯(lián)合反演。

3) 鐵路隧道勘探結(jié)果需要相對(duì)精細(xì)的橫向變化，所以二維NLCG反演不能片面追求模型光滑，在誤差收斂、結(jié)果不影響背景地層的情況下，應(yīng)盡量使用較小的正則化因子。

4) 二維NLCG反演的初始模型一般是均勻半空間，其電阻率值最好和中、淺層電阻率接近，若初始模型選擇不合適，在TE+TM模式聯(lián)合反演時(shí)可能會(huì)導(dǎo)致在深部斷層傾向不明確的結(jié)果。有條件的情況下，可以參考較大正則化因子的二維NLCG反演結(jié)果來(lái)確定中、淺部的電阻率以正確選擇初始模型。

必須提出，文中所有結(jié)論是基于特定模型得到的，對(duì)于更復(fù)雜的模型或者類(lèi)A形、H型或Q型的地電斷面，可能局部結(jié)論和本文不同，這些斷面有待更多學(xué)者展開(kāi)更多的細(xì)致研究，以提供更多鐵路勘探AMT數(shù)據(jù)二維NLCG反演參數(shù)選取參考。