宋婉婷,陳衛(wèi)營*

(1.中國科學(xué)院地質(zhì)與地球物理研究所 中國科學(xué)院礦產(chǎn)資源研究重點實驗室,北京 100029;2.中國科學(xué)院地球科學(xué)研究院,北京 100029;3.中國科學(xué)院大學(xué) 地球與行星科學(xué)學(xué)院,北京 100049)

0 引 言

電性源短偏移瞬變電磁法(Short-offset Transient Electromagnetics Method,SOTEM)利用接地導(dǎo)線源向地下發(fā)射激勵電流，在0.3～2.0倍探測深度的偏移距范圍內(nèi)觀測純二次電磁場信號，并通過處理與解譯觀測信號達到探測地下目標(biāo)體的目的。大量研究和應(yīng)用表明，SOTEM具有探測深度大(一般可達2 km)、分辨率高、野外施工靈活等優(yōu)點，在深部金屬礦探測、煤田水文地質(zhì)調(diào)查等領(lǐng)域表現(xiàn)出巨大潛力。

受數(shù)值模擬手段和計算機計算效率的限制，SOTEM數(shù)據(jù)的多維反演尚未完全實用化，目前還以一維反演為主。傳統(tǒng)的一維反演方法是對一條測線上的測點逐一進行單點反演，反演結(jié)果只是利用層狀大地模型簡單的拼接來描述地下介質(zhì)的電性信息。在實際測量中，不同測點通常會具有不同的地形和噪聲環(huán)境，這種單點的反演方式通常會造成相鄰測點反演結(jié)果不連續(xù)，產(chǎn)生虛假的局部異常，造成反演解釋的困難。為了解決單點一維反演的電阻率橫向突變問題，同時保留一維反演速度快的優(yōu)勢，Auken等探索出對反演模型或數(shù)據(jù)施加約束條件的擬二維反演方法，并成功應(yīng)用于多種電磁數(shù)據(jù)的反演解釋中。Auken等最早將橫向約束用于直流電法數(shù)據(jù)的反演中，后又在最小二乘法二維反演中改進了約束矩陣，結(jié)果表明擬二維反演在提高分辨率和壓制噪聲方面都有所改善。Siemon等首次將橫向約束反演技術(shù)應(yīng)用于直升機載頻域電磁數(shù)據(jù)中。Santos等提出了一種基于橫向約束算法的直流電阻率和時域電磁數(shù)據(jù)的聯(lián)合反演方法，并驗證了該方法的實際應(yīng)用價值。殷長春等將橫向約束方法進一步改進，提出了加權(quán)橫向約束(WLCI)思想，并應(yīng)用于頻率域航空電磁數(shù)據(jù)的擬二維反演中，后又將該方法用于時間域航空電磁擬二維反演。近期，基于OCCAM反演的光滑擬二維反演方法和基于量子行為粒子群算法的擬二維反演技術(shù)也應(yīng)用于航空瞬變電磁數(shù)據(jù)的處理中。

在SOTEM勘探中，測線布置一般較長，測點多且密集，由同一發(fā)射源激發(fā)所得的數(shù)據(jù)穩(wěn)定性較好。其尤其是在煤礦領(lǐng)域應(yīng)用廣泛，沉積煤系地層橫向連續(xù)性好，同一測線上的電阻率和厚度在橫向上變化一般不大。將約束思想用于SOTEM數(shù)據(jù)處理中能大大提高反演結(jié)果的縱、橫向連續(xù)性，并且在一定程度上降低大深度探測中晚期信號易受干擾的影響。本文在前人研究基礎(chǔ)上，實現(xiàn)SOTEM數(shù)據(jù)的擬二維反演，并通過對理論模型數(shù)據(jù)和實測數(shù)據(jù)的反演評估擬二維反演在SOTEM數(shù)據(jù)處理中的有效性，旨在提供一種可靠、穩(wěn)定、實用的SOTEM數(shù)據(jù)反演方法。

1 SOTEM一維響應(yīng)計算

SOTEM一維正演從水平電偶極子源出發(fā)，獲得具有積分形式的接地長導(dǎo)線源的頻率域響應(yīng)，再通過傅氏變換轉(zhuǎn)化為時間域的響應(yīng)。置于地表的電偶極源在直角坐標(biāo)系下可以產(chǎn)生三磁(、、)、三電(、、)共6個電磁場分量。已有研究表明水平電場分量和垂直磁場強度分量的場值分布、衰減更利于觀測和處理。本文研究選用野外實際觀測中最為常用的垂直磁場強度分量(或其時間導(dǎo)數(shù)d/d)。直角坐標(biāo)系下，置于地表的水平電偶極值在層狀大地表面產(chǎn)生的垂直磁場強度計算公式為

(1)

當(dāng)大地磁導(dǎo)率等于自由空間中的磁導(dǎo)率時，反射系數(shù)表達式為

(2)

對于層大地，有遞推公式為

(3)

(4)

通過偶極子疊加得到有限長導(dǎo)線源的電磁場，其表達式為

(5)

式中：為長導(dǎo)線源分割成的偶極子個數(shù)；Δ為每個偶極子的長度；為接收點到每個偶極子中心的距離。

2 擬二維反演方法

擬二維反演的主要思想是通過在相鄰測點間施加約束條件(電阻率、厚度或深度)，以保證電阻率斷面圖的縱、橫向連續(xù)性，其出發(fā)點仍為一維反演。基于正則化思想，一維反演的目標(biāo)函數(shù)可歸結(jié)為

()=()+()→min

(6)

式中：()為數(shù)據(jù)目標(biāo)函數(shù)；()為模型先驗信息目標(biāo)函數(shù)；為正則化調(diào)節(jié)因子。

在一維正則化反演理論的基礎(chǔ)上，引入模型橫向粗糙度目標(biāo)函數(shù)(())和縱向粗糙度目標(biāo)函數(shù)(())，確定擬二維整體反演的目標(biāo)函數(shù)。其表達式為

()=()+[()+()+()]

(7)

式中：、、分別為模型先驗信息、模型橫向粗糙度、模型縱向粗糙度的正則化系數(shù)。

整條測線各測點的大地模型參數(shù)按一列排布，=[,…,1,,,…,1,-1,…,,1,…,,,,1,…,,-1]，其中,為第個測點第層的電導(dǎo)率，,-1為第個測點第-1層的厚度。

觀測數(shù)據(jù)的目標(biāo)函數(shù)為

(8)

(9)

(10)

式中：Δ為相鄰測點間距離；,為第個測點第層厚度。

在反演迭代過程中，通過線性搜索自適應(yīng)迭代的方式自動確定，即在每次迭代過程中，通過線性搜索的方法，自動選取最優(yōu)的正則化因子。、、在初始反演時根據(jù)實際情況選擇合適取值，在反演過程中保持不變。

為求得目標(biāo)函數(shù)的最優(yōu)解，以第+1次迭代為例，令目標(biāo)函數(shù)對待求模型參數(shù)+1的偏導(dǎo)等于0。其表達式為

(11)

+1=

(12)

求解上述線性方程組，即可得到第+1次迭代所求得的新模型向量+1。反演的終止條件則為()≤1，其中為測線電磁數(shù)據(jù)個數(shù)。

3 理論模型擬二維反演

首先通過一個如圖1所示的二維地電理論模型來驗證擬二維反演處理SOTEM數(shù)據(jù)的效果。模型共分為3層：第一層厚度為200 m，電阻率()為200 Ω·m；第二層底界面起伏變換，厚度為200～400 m，電阻率()為50 Ω·m；第三層為基底，電阻率()為500 Ω·m。利用單點一維正演獲得整條剖面的響應(yīng)數(shù)據(jù)，正演中發(fā)射源長度為1 000 m，測線長度為1 000 m，點距為20 m，偏移距為500 m，發(fā)射電流為1 A，計算時間為0.1～100.0 ms，共分41個時間道，計算場量為垂直磁場強度分量。正演計算中收發(fā)布置如圖2所示。

圖1 理論模型示意圖Fig.1 View of Theoretical Model

圖2 一維正演收發(fā)布置圖Fig.2 Layout of Transmitter and Receivers for 1D Forward

對正演獲得的響應(yīng)數(shù)據(jù)分別進行單點一維反演和擬二維反演。初始模型為100 Ω·m的均勻半空間，反演最大深度取1 100 m，首層厚度10 m，往下各層厚度以上一層1.1倍數(shù)遞增，共包含30層，迭代次數(shù)設(shè)置為7次。擬二維反演中，正則化系數(shù)、、取值都設(shè)為1，正則化因子初始值設(shè)為1 000，在迭代過程中利用線性搜索算法自適應(yīng)選取最優(yōu)的正則化因子。首先對正演得到的無噪聲數(shù)據(jù)進行反演，從而得到單點一維反演和擬二維反演電阻率-深度斷面圖(圖3)。由圖3(a)、(b)可以看出，當(dāng)數(shù)據(jù)不含噪聲時，一維反演與擬二維反演結(jié)果基本吻合，都很好地恢復(fù)了真實模型的電阻率分布，對層界面刻畫也非常清晰。圖3(c)給出了兩種反演模式下整條測線的擬合殘差曲線，可以看出兩種反演的擬合殘差都收斂到較低水平(最大值不超過3.6%)，且擬二維反演情況下擬合殘差更小。

圖3 無噪聲反演結(jié)果Fig.3 Inversion Results Without Noise

為研究噪聲水平對反演結(jié)果的影響，本文在所有測點的正演數(shù)據(jù)中按照時間道早晚依次添加不同程度的白噪聲。第一組噪聲水平為1%～50%，第二組為1%～100%。分別對兩組加噪后的數(shù)據(jù)進行單點一維反演和擬二維反演，反演參數(shù)與前述無噪情況一致，獲得整條剖面的電阻率-深度斷面圖(圖4、5)。從圖4、5可以看出，當(dāng)反演數(shù)據(jù)含噪時，單點一維反演結(jié)果雖仍能反映出模型的大致電性分布，但受噪聲影響，電性層分界面變得毛糙，橫向連續(xù)性變差。特別是當(dāng)噪聲水平較高時，出現(xiàn)了相鄰測點電阻率和層厚度的突變，空間連續(xù)性變差，產(chǎn)生了虛假的局部異常(圖5)。而擬二維反演結(jié)果受噪聲的影響相對較弱，仍能很好地反映模型的真實電性結(jié)構(gòu)，對模型電性界面的刻畫清晰且光滑，很大程度上降低了噪聲引起的這種電阻率縱、橫向突變，使得僅由個別測點引起的假異?，F(xiàn)象得到很好的抑制。通過觀察擬合殘差可以發(fā)現(xiàn)，在相同的迭代次數(shù)情況下，擬二維反演并不一定能減小最終的擬合殘差，且擬二維反演的擬合殘差普遍大于單點一維反演結(jié)果。這是因為在擬二維反演中，反演模型受橫向電阻率約束，并不能完全擬合曲線的實際形態(tài)。當(dāng)個別測點的曲線出現(xiàn)嚴(yán)重畸變時，擬二維反演的擬合殘差可能會很大，但會避免追求低擬合殘差帶來的錯誤結(jié)果。

圖4 加入1%～50%白噪聲的反演結(jié)果Fig.4 Inversion Results with 1%-50% White Noise

圖5 加入1%～100%白噪聲的反演結(jié)果Fig.5 Inversion Results with 1%-100% White Noise

同時，由于擬二維反演同時處理的是整條剖面的數(shù)據(jù)，當(dāng)總迭代次數(shù)相同(7次)時，相較于單點一維反演，擬二維反演需要進行的正演計算和雅克比矩陣求解次數(shù)更少，從而反演所需時間更少，效率更高。表1給出了上述3種噪聲水平下，采用單點一維反演和擬二維反演處理整條剖面數(shù)據(jù)的耗時。計算環(huán)境為：Intel(R)Xeon(R)CPU E5-1603 v4@ 2.80 GHz處理器，16 GB內(nèi)存。由表1可知，擬二維反演可顯著提升數(shù)據(jù)反演的效率。

表1 單點一維反演與擬二維反演耗時Table 1 Time Consuming of Single-site 1D Inversion and Quasi-2D Inversion

最后，本文測試了擬二維反演對復(fù)雜模型的反演效果。設(shè)計如圖6所示的一個雙異常體模型。對該模型進行一維正演獲得垂直感應(yīng)電動勢d/d響應(yīng)，并對數(shù)據(jù)添加和圖5一樣的高水平噪聲。對含噪數(shù)據(jù)分別進行單點一維反演和擬二維反演，反演參數(shù)與前述一致，結(jié)果如圖7所示。從圖7可以看出，對于含兩個異常體的復(fù)雜模型，當(dāng)噪聲水平較高時，采用單點一維反演很難清晰地反演出兩個異常體的形態(tài)，反演結(jié)果出現(xiàn)了大量的假異常[圖7(a)]；而當(dāng)采用擬二維反演時，雖然也存在個別的假異常，但目標(biāo)異常體的橫向范圍得到了很好的恢復(fù)[圖7(b)]；從擬合殘差曲線來看，擬二維反演的擬合殘差在200～400號點之間存在一個高值區(qū)域[圖7(c)]，這是因為該范圍內(nèi)的響應(yīng)畸變較嚴(yán)重，在橫向約束下反演迭代中程序未完全按照原始數(shù)據(jù)的形態(tài)進行擬合。

圖6 雙異常體復(fù)雜模型示意圖Fig.6 View of Complex Model with Two Anomaly Bodies

圖7 含噪復(fù)雜模型反演結(jié)果Fig.7 Inversion Results of Complex Model with Noise

4 實測數(shù)據(jù)擬二維反演

為進一步驗證擬二維約束反演在處理SOTEM數(shù)據(jù)中的有效性和優(yōu)越性，本文對山西晉城地區(qū)某煤礦的實測SOTEM數(shù)據(jù)進行反演處理。測區(qū)位于沁水煤田高平礦區(qū)，探測目的是調(diào)查3#煤富水采空區(qū)分布和區(qū)內(nèi)富水導(dǎo)水通道，并評價奧陶系灰?guī)r地層含水性。測區(qū)主要地層及電性特征如表2所示。

表2 測區(qū)地層綜合電性特征Table 2 Comprehensive Electrical Properties of Strata in the Survey Area

本次SOTEM探測工作共布置測線42條，編號為L1～L42，測線方向為EW向，測線間距為50 m，測點間距為20 m。本文以其中兩條測線L4和L18為例，分析擬二維反演效果。發(fā)射源與測線布置及相關(guān)幾何參數(shù)如圖8所示，發(fā)射電流為12 A，基頻為12.5 Hz，觀測垂直感應(yīng)電壓()分量，接收線圈有效面積為100 m。觀測儀器為加拿大鳳凰公司的V8多功能電法工作站，測點干擾較小且數(shù)據(jù)較為光滑時，單點觀測時間為2 min，測點電磁干擾較大時，觀測時間為5 min。實測數(shù)據(jù)質(zhì)量整體尚好，但由于L4測線靠近村莊，受附近民用電線及居民電器影響，數(shù)據(jù)質(zhì)量相對較差，大部分測點的晚期信號信噪比較低，個別測點甚至在全時間段內(nèi)出現(xiàn)嚴(yán)重畸變。兩條測線的實測多測道感應(yīng)電壓曲線如圖9所示。本文在進行反演處理時，對晚期畸變嚴(yán)重的時間道數(shù)據(jù)進行了剔除。

圖8 測線及發(fā)射源布置圖Fig.8 Layout of Survey Line and Emission Source

圖9 實測感應(yīng)電壓多測道曲線Fig.9 Measured Multichannel Curves of Induced Voltage

利用SOTEM單點一維反演和擬二維反演技術(shù)對經(jīng)預(yù)處理后的每條測線數(shù)據(jù)進行反演。初始模型為100 Ω·m的均勻半空間，反演最大深度為600 m，首層厚度為10 m，以下各層厚度按上層厚度的1.05倍增加，共分為29個電性層，擬合次數(shù)設(shè)置為10次。擬二維反演中，先驗信息采用一維等效源反演結(jié)果，正則化系數(shù)、、取值為1、10、1，正則化因子初始值設(shè)為1 000。圖10和圖11為得到的高程-電阻率斷面圖。由圖10、11可以看出，對于每條測線，單點一維反演和擬二維反演結(jié)果對地層電性結(jié)構(gòu)的反映在整體上基本一致。尤其是對于淺部地層，兩種反演方法對電性結(jié)構(gòu)變化細節(jié)的反映也能很好的吻合。但是隨著深度的增大，電磁信號分辨率降低，單點一維反演的不確定性導(dǎo)致其個別測點的結(jié)果在橫向上出現(xiàn)較明顯的不連續(xù)性。特別是當(dāng)原始數(shù)據(jù)受到干擾程度較高時，深部電阻率橫向突變現(xiàn)象尤為嚴(yán)重。L4測線60、160和260 m處出現(xiàn)縱向延伸的局部低阻帶[圖10(a)]，L18測線0～350 m范圍內(nèi)也出現(xiàn)多處這種“掛面條”式低阻帶[圖11(a)]，這些異?，F(xiàn)象都有悖于煤礦區(qū)域地下結(jié)構(gòu)橫向連續(xù)性高的特點，給解釋帶來困難。而經(jīng)過橫向約束的擬二維反演結(jié)果展現(xiàn)了更接近真實情況的地下電性結(jié)構(gòu)，在很大程度上降低了這種電阻率橫向上的突變，使得僅由個別測點引起的假異常現(xiàn)象得到很好的抑制[圖10(b)、11(b)]。為了驗證反演結(jié)果的可靠性，本文將礦方在L18測線320號點附近實施的一個鉆孔(ZK-206)結(jié)果繪制于圖11中。根據(jù)鉆孔揭示，該位置在鉆進深度245 m范圍內(nèi)未發(fā)現(xiàn)采空區(qū)。煤系地層深度范圍內(nèi)地層應(yīng)呈高阻反映，這與擬二維反演結(jié)果更為接近。另外，根據(jù)L18測線反演結(jié)果推斷的兩條斷層(圖11中F1和F2斷層)與真實斷層的位置基本一致，進一步驗證了反演結(jié)果的可靠性。

圖10 L4測線反演結(jié)果Fig.10 Inversion Results of Line L4

圖11 L18測線反演結(jié)果Fig.11 Inversion Results of Line L18

5 結(jié) 語

本文針對SOTEM的特點和需求，研究了SOTEM數(shù)據(jù)擬二維約束反演，采用自適應(yīng)正則化算法，在單點一維反演的基礎(chǔ)上引入模型粗糙度目標(biāo)函數(shù)，對相鄰測點之間模型參數(shù)進行橫向及縱向約束。通過對理論模型數(shù)據(jù)和野外實測數(shù)據(jù)的處理，驗證了該方法可克服單點一維反演結(jié)果中電阻率橫向連續(xù)性差的問題，有效壓制噪聲對反演結(jié)果的影響，大大提高了反演精度，且計算效率更高。