SOTEM數(shù)據(jù)擬二維反演研究與應(yīng)用

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

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

0 引 言

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

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

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

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

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

(1)

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

(2)

對(duì)于層大地，有遞推公式為

(3)

(4)

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

(5)

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

2 擬二維反演方法

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

()=()+()→min

(6)

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

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

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

(7)

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

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

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

(8)

(9)

(10)

式中：Δ為相鄰測(cè)點(diǎn)間距離；,為第個(gè)測(cè)點(diǎn)第層厚度。

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

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

(11)

+1=

(12)

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

3 理論模型擬二維反演

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

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

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

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

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

為研究噪聲水平對(duì)反演結(jié)果的影響，本文在所有測(cè)點(diǎn)的正演數(shù)據(jù)中按照時(shí)間道早晚依次添加不同程度的白噪聲。第一組噪聲水平為1%～50%，第二組為1%～100%。分別對(duì)兩組加噪后的數(shù)據(jù)進(jìn)行單點(diǎn)一維反演和擬二維反演，反演參數(shù)與前述無噪情況一致，獲得整條剖面的電阻率-深度斷面圖(圖4、5)。從圖4、5可以看出，當(dāng)反演數(shù)據(jù)含噪時(shí)，單點(diǎn)一維反演結(jié)果雖仍能反映出模型的大致電性分布，但受噪聲影響，電性層分界面變得毛糙，橫向連續(xù)性變差。特別是當(dāng)噪聲水平較高時(shí)，出現(xiàn)了相鄰測(cè)點(diǎn)電阻率和層厚度的突變，空間連續(xù)性變差，產(chǎn)生了虛假的局部異常(圖5)。而擬二維反演結(jié)果受噪聲的影響相對(duì)較弱，仍能很好地反映模型的真實(shí)電性結(jié)構(gòu)，對(duì)模型電性界面的刻畫清晰且光滑，很大程度上降低了噪聲引起的這種電阻率縱、橫向突變，使得僅由個(gè)別測(cè)點(diǎn)引起的假異?，F(xiàn)象得到很好的抑制。通過觀察擬合殘差可以發(fā)現(xiàn)，在相同的迭代次數(shù)情況下，擬二維反演并不一定能減小最終的擬合殘差，且擬二維反演的擬合殘差普遍大于單點(diǎn)一維反演結(jié)果。這是因?yàn)樵跀M二維反演中，反演模型受橫向電阻率約束，并不能完全擬合曲線的實(shí)際形態(tài)。當(dāng)個(gè)別測(cè)點(diǎn)的曲線出現(xiàn)嚴(yán)重畸變時(shí)，擬二維反演的擬合殘差可能會(huì)很大，但會(huì)避免追求低擬合殘差帶來的錯(cuò)誤結(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í)，由于擬二維反演同時(shí)處理的是整條剖面的數(shù)據(jù)，當(dāng)總迭代次數(shù)相同(7次)時(shí)，相較于單點(diǎn)一維反演，擬二維反演需要進(jìn)行的正演計(jì)算和雅克比矩陣求解次數(shù)更少，從而反演所需時(shí)間更少，效率更高。表1給出了上述3種噪聲水平下，采用單點(diǎn)一維反演和擬二維反演處理整條剖面數(shù)據(jù)的耗時(shí)。計(jì)算環(huán)境為：Intel(R)Xeon(R)CPU E5-1603 v4@ 2.80 GHz處理器，16 GB內(nèi)存。由表1可知，擬二維反演可顯著提升數(shù)據(jù)反演的效率。

表1 單點(diǎn)一維反演與擬二維反演耗時(shí)Table 1 Time Consuming of Single-site 1D Inversion and Quasi-2D Inversion

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

圖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í)測(cè)數(shù)據(jù)擬二維反演

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

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

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

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

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

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

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

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

5 結(jié) 語

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