水平大定源回線瞬變電磁任意點一維正反演

石金晶，薛國強,4*，陳衛(wèi)營，武 欣

(1. 中國科學院地質(zhì)與地球物理研究所 中國科學院礦產(chǎn)資源研究重點實驗室，北京 100029; 2. 中國科學院地球科學研究院，北京 100029; 3. 中國科學院大學 地球與行星科學學院，北京 100049; 4. 西北有色地質(zhì)礦業(yè)集團有限公司，陜西 西安 710054)

0 引 言

回線源瞬變電磁法(TEM)回線裝置分為中心回線、重疊回線、分離回線、大定源回線等。中心回線裝置由于其理論計算的準確性和對復雜工作環(huán)境的適應(yīng)性被廣泛應(yīng)用于地下水和環(huán)境調(diào)查，并可作為煤炭和金屬礦產(chǎn)的勘探工具。在回線裝置中心某一區(qū)域接收到的信號沒有明顯的差異，因此，可以觀測大回線中心位置一定范圍內(nèi)的瞬變電磁響應(yīng)，這被稱為改進的中心回線裝置。對于改進的中心回線裝置，視電阻率的求取通常是利用中心回線裝置的公式進行近似計算。然而，這種處理方式往往會導致計算出現(xiàn)系統(tǒng)性的誤差而變得不適用。事實上，不僅線框內(nèi)部有感應(yīng)場出現(xiàn)，線框外部也有。為進一步提高探測效率，最好的解決方式是實現(xiàn)在一個大的固定回線框內(nèi)外的任意位置接收信號，即實現(xiàn)大定源回線裝置任意點測量。

大定源回線裝置一般采用較大邊長的矩形回線作為發(fā)射裝置。相較于中心回線裝置和重疊回線裝置只在中心點觀測的特點，大定源回線裝置一旦鋪設(shè)好回線裝置后，可以多臺接收機同時工作，不用頻繁挪動線框，工作效率較高。而且，大定源回線瞬變電磁探測技術(shù)成熟，已經(jīng)在瞬變電磁正反演研究方面取得了相當大的進展，尤其是一維正反演。在早期關(guān)于瞬變電磁的正演研究中多使用圓形發(fā)射裝置模擬野外多邊形(通常是矩形)發(fā)射線框。為方便計算，Singh等利用超幾何分布函數(shù)簡化了圓形發(fā)射裝置的電磁場表達式中的貝塞爾函數(shù)，從而得到了任意點圓形發(fā)射框的頻率域電磁場表達式。Xue等將積分范圍分段近似，利用雙貝塞爾函數(shù)得到了圓形發(fā)射回線內(nèi)外電磁場的表達式。

對于矩形發(fā)射線框，圓形線框的公式不再適用，尤其是計算早期響應(yīng)的時候，因此，必須考慮發(fā)射回線形狀給出新的計算方法。事實上，無論回線的形狀如何，回線內(nèi)任意一點的電磁場都可以看作回線各邊在該點的電磁場之和。為方便計算電磁場，需要對回線源進行剖分，以滿足水平電偶極子近似條件，即剖分后的線源到觀測點的距離遠大于剖分后線源本身的長度，再將水平電偶極子進行數(shù)值積分就可得到整個線源的電磁場。Poddar利用電磁互易原理和磁偶極子源得到了頻率域電磁場表達式，并對水平電偶極子的頻率域響應(yīng)進行反漢克爾變換，然后將結(jié)果進行傅里葉反變換轉(zhuǎn)化為時域。Raiche首先采用嵌套內(nèi)插值方法得到了一個頻域多邊形回線發(fā)射裝置的電磁場表達式，然后通過Gaver-Stehfest變換推導出相應(yīng)的時域電磁場表達式。Goldman等得到了水平電偶極子源雙層介質(zhì)表面的電磁場表達式，但利用水平電偶極子近似矩形回線時的場表達式比較復雜。此外，國內(nèi)學者也在電偶極子近似方面做了不少工作。例如，華軍等將貝塞爾函數(shù)的積分區(qū)間分成兩段，分別用不同的數(shù)值方法進行求解，明顯提高了對貝塞爾函數(shù)積分的計算效率；李建慧等沿用偶極子合成回線源的思想，獲得了大回線源的瞬變電磁響應(yīng)；李展輝等利用水平電偶極子合成原理，發(fā)展了一種適應(yīng)于任意形狀水平接地導線源瞬變電磁法的一維正演方法。趙越等利用貝塞爾函數(shù)展開法計算了大回線源非中心點垂直分量與水平分量的電磁響應(yīng)，對于框外存在過渡區(qū)“雙解”及“無解”的問題，提出了相應(yīng)的解決方案。

大定源回線裝置雖然提出已久，但是目前國內(nèi)實際生產(chǎn)中主要還是利用中心回線裝置，或者改進的中心回線裝置，觀測區(qū)域局限于發(fā)射線框內(nèi)中間區(qū)域的一定范圍內(nèi)，少有對大定源回線框外響應(yīng)特性及其反演效果的研究。基于此，本文利用畢奧-薩法爾(Biot-Savart)定理計算了大定源回線的一次場，利用偶極子疊加原理計算了均勻半空間大定源回線框內(nèi)外的二次場響應(yīng)，通過二次場響應(yīng)發(fā)現(xiàn)晚期段二次場的分布趨于均勻，從而得出在實際觀測中可以選擇在框內(nèi)外一定范圍觀測二次場；在考慮發(fā)射框的形狀和觀測的實際位置情況下，使用該方法對一個三層地質(zhì)模型進行正反演，獲得較好的效果；最后，將該方法應(yīng)用到湖北省荊州市某機場的大定源回線勘探數(shù)據(jù)中，取得了較好的效果，反演結(jié)果與已有資料具有較好的一致性。

1 理論背景

1.1 一維正演

對于具有一定尺寸的矩形回線發(fā)射源，可以采用偶極子疊加的方式來計算其產(chǎn)生的電磁場。疊加方式一般有兩種：一種是采用磁偶極子源進行面積分；另一種是采用電偶極子進行線積分。由于采用電偶極子源沿發(fā)射線框進行線積分更易于模擬矩形線框，故本文采用這種疊加方式。

圖1給出了大定源回線裝置電磁響應(yīng)計算的示意圖?；鼐€框的中心在(0，0)處，在方向和方向上的邊長分別為2和2，接收裝置可以位于回線框的內(nèi)部或者外部的任意位置?；鼐€框的4個頂點分別用、、和來表示。

圖1 大定源回線裝置電磁響應(yīng)計算示意圖Fig.1 Schematic Diagram of Electromagnetic Response Calculation of Large-loop Fixed Source Device

納比吉安給出了水平電偶極子在均勻?qū)訝罱橘|(zhì)的大地表面上的頻率域垂直磁場分量(())表達式。其表達式為

(1)

將式(1)進行積分，可以求得長度為2的線源垂直磁場分量為

(2)

式(2)為關(guān)于貝塞爾函數(shù)的無窮積分，無解析解，必須使用數(shù)值方法求解。選用濾波系數(shù)法求其數(shù)值解，將積分形式轉(zhuǎn)化為累加形式。其表達式

(3)

式中：為線源剖分的偶極子總數(shù)；為第個偶極子的收發(fā)距；偶極子剖分長度按照d=′50的規(guī)則進行確定，其中′為觀測點到發(fā)射源某個邊長的垂直距離。

由式(3)可獲得回線框各邊的垂直磁場。將各邊的垂直磁場相加即可獲得整個回線框產(chǎn)生的總垂直磁場為

(4)

圖1所示回線框的總垂直磁場可表示為

(5)

求得頻率域場值后，可通過余弦變換將其轉(zhuǎn)換到時間域，得到時刻的階躍電磁響應(yīng)。其表達式為

(6)

式中：()為時間域垂直磁場響應(yīng)。

針對式(4)中涉及的貝塞爾函數(shù)積分，本文采用Key等給出的201點濾波系數(shù)進行計算。針對式(6)中涉及的余弦變換，本文采用Anderson等給出的787點濾波系數(shù)進行計算。

1.2 一維反演

OCCAM反演結(jié)果受初始模型影響小，運算穩(wěn)定收斂。本文采用一階粗糙度定義模型光滑程度的粗糙度。模型向量數(shù)為時的一階粗糙度()為

(7)

式中：=(,,…,),代表模型參數(shù)向量；是向量的第個分量。

引入粗糙度矩陣，定義為

(8)

模型粗糙度可表示為

=‖‖

(9)

定義模型數(shù)據(jù)與實測數(shù)據(jù)的擬合殘差()，使其遵從加權(quán)最小二乘法的原理。其表達式為

=‖-F[]‖

(10)

式中：=diag{1,1,,1}，是誤差加權(quán)矩陣，1為第個數(shù)據(jù)點的方差；=(,,…,)代表實測數(shù)據(jù)向量；[]是模型參數(shù)向量對應(yīng)的正演數(shù)據(jù)。

(11)

在反演迭代時，采用局部線性化思想把非線性問題轉(zhuǎn)化為線性問題，對模型參數(shù)進行泰勒展開，得到局部近似式子。其表達式為

[+Δ]≈[]+()Δ

(12)

(13)

式中：Δ為向量的修正量；()是×階雅克比矩陣，和分別為模型向量和數(shù)據(jù)向量的個數(shù)；[]為正演數(shù)據(jù)向量的第個分量，=1,2，…，。

利用差分思想對式(13)進行計算，將最小化目標函數(shù)轉(zhuǎn)化為

(14)

當系數(shù)矩陣滿秩時，可求得方程組的解為

(15)

2 響應(yīng)特性

2.1 初始場分布特征

瞬變電磁二次場以及一次場的分布與初始場的分布密切相關(guān)。在研究二次場特性之前，首先計算矩形回線源在框內(nèi)外產(chǎn)生的初始場分布特征。初始場(磁感應(yīng)強度)與模型的電阻率無關(guān)，其場值可根據(jù)畢奧-薩法爾定理計算，設(shè)回線框中電流為1 A，計算網(wǎng)格為5×5 m，結(jié)果如圖2所示。圖2中，“+”表示區(qū)域值為正，“-”表示區(qū)域值為負。從圖2可以看出，框內(nèi)外初始場的方向相反，初始場幅值最高的區(qū)域集中在發(fā)射線框兩側(cè)很近的距離范圍內(nèi)，框內(nèi)中心區(qū)域的場值分布較為均勻，框外區(qū)域隨著距離的增大場值迅速衰減，框外150 m處初始場的幅值相比框內(nèi)中心點的幅值已下降近2個數(shù)量級。

圖2 一次場平面分布Fig.2 Plane Distribution of Primary Field

2.2 二次場響應(yīng)特征

以電阻率為100 Ω·m的均勻半空間模型為例分析二次場的分布與擴散特性，發(fā)射源尺寸和計算區(qū)域與圖3一致。發(fā)射電流為1 A，計算時窗為0.01～100.00 ms，共51個時間道，計算場量為垂直感應(yīng)電動勢d/d。其中，為垂直磁感應(yīng)強度。

圖3 不同時刻二次場平面分布Fig.3 Plane Distributions of Secondary Field at Different Times

首先分析二次場的分布特性，選擇所有計算點在0.01、0.1、1和10 ms等4個時刻的響應(yīng)值，繪制成平面等值圖(圖3)。圖3中，“+”表示區(qū)域值為正，“-”表示區(qū)域值為負。早期時刻(=0.01 ms)，二次場高值區(qū)域主要集中在框內(nèi)，框外距線框較近處出現(xiàn)符號反轉(zhuǎn)現(xiàn)象；隨著時間推移，高值區(qū)域逐漸變大，信號反轉(zhuǎn)位置也向外移動[圖3(b)、(c)]；當觀測時刻為10 ms時，計算區(qū)域內(nèi)的二次場分布趨于均勻，不同位置處的響應(yīng)差別非常小(均方差僅為0.007)。

圖4給出了框內(nèi)外6個測點的衰減曲線，可以更好地觀察二次場隨時間的變化特征。當測點在框內(nèi)時，二次場不發(fā)生符號反轉(zhuǎn)現(xiàn)象，但測點位置不同會導致信號在早期幅值不同，靠近中心點處的幅值更高。當測點在框外但距離線框很近時(如=105點)，在計算時窗范圍內(nèi)信號也未發(fā)生符號反轉(zhuǎn)，其衰減趨勢與框內(nèi)基本一致；但當測點離發(fā)射線框較遠時(如=150點和=200點)，早期二次場響應(yīng)為負，并在某個時刻信號符號轉(zhuǎn)變?yōu)檎?。距離線框位置越遠，負值響應(yīng)的時間范圍越長，信號符號反轉(zhuǎn)的時刻也越久。這種框外測點響應(yīng)符號反轉(zhuǎn)的現(xiàn)象可用煙圈的擴散來解釋。無論在框內(nèi)還是框外，隨著觀測時間推移，各測點的衰減曲線趨于重合，表明晚期二次場在空間上的分布趨于均勻，發(fā)射源和接收點的幾何關(guān)系帶來的影響已非常小。

x=0,50,95點為框內(nèi)點；x=105,150,200點為框外點圖4 框內(nèi)外不同測點的二次場衰減曲線Fig.4 Secondary Field Attenuation Curves of Different Measuring Points Inside and Outside the Frame

上述分析表明，觀測點的空間位置不同會導致響應(yīng)的幅值和形態(tài)發(fā)生變化，尤其是框外的測點在早期會發(fā)生信號符號反轉(zhuǎn)現(xiàn)象，但隨著時間推移，觀測點空間位置帶來的影響逐漸減弱，晚期二次場的分布趨于均勻。因此，在實際觀測中，可以選擇在框外一定范圍觀測二次場。

3 一維反演

針對大定源回線裝置的數(shù)據(jù)處理，不能套用中心回線裝置的處理方式，必須考慮發(fā)射源的實際形狀和觀測的具體位置。在考慮這兩個因素前提下，基于一維正演方法，利用OCCAM算法對數(shù)據(jù)進行一維反演，并評估反演效果。OCCAM反演是一種帶約束性條件的反演方法，它在追求模擬數(shù)據(jù)與實測數(shù)據(jù)最大擬合的同時，要求模型數(shù)據(jù)最平滑，跳出了傳統(tǒng)梯度方法簡單求取模型增量的模式，而是在對拉格朗日乘子值的搜索過程中尋找靠近極值點的目標模型，保證了算法始終具有很好的穩(wěn)定性。反演過程中，本文用一維線性搜索法求取拉格朗日算子。

建立一個三層地電模型，各層厚度分別為300和100 m以及無窮大，各層電阻率分別為100、30和80 Ω·m。發(fā)射源尺寸為200 m×100 m，測線長度為500 m，點距為5 m，發(fā)射源和接收點的幾何關(guān)系如圖5所示。計算電磁場分量為垂直感應(yīng)電動勢d/d，計算時窗范圍為0.01～100 ms，共51個時間道。

圖5 發(fā)射源-接收點裝置示意圖Fig.5 Schematic Diagram of Transmitting-receiving Device

一維正演完成后，對整條線的正演數(shù)據(jù)進行OCCAM一維反演。設(shè)置反演模型最大深度為800 m，共分50層，層厚由淺至深按指數(shù)等間隔增加，迭代次數(shù)設(shè)為10次。結(jié)果表明：反演過程中，經(jīng)過5次迭代后各點的擬合殘差都已經(jīng)下降到2%以內(nèi)。圖6給出了框內(nèi)外6個測點的數(shù)據(jù)擬合情況。從圖6可以看出，經(jīng)過10次迭代，無論是框內(nèi)還是框外的測點都擬合得非常好。圖7給出了上述6個測點的反演結(jié)果，其表現(xiàn)出了很好的一致性，說明考慮了發(fā)射源和接收點實際幾何位置的反演，可以很好地恢復地層的真實電性結(jié)構(gòu)。

x=0,50,95點為框內(nèi)點；x=105,150,200點為框外點圖6 不同測點的數(shù)據(jù)擬合情況Fig.6 Data Fitting Curves of Different Measuring Points

圖7 框內(nèi)外典型測點一維OCCAM反演結(jié)果Fig.7 One-dimensional Occam Inversion Results of Typical Measuring Points Inside and Outside the Frame

4 應(yīng)用實例

在湖北省荊州市某機場進行了大定源回線裝置的探測試驗，并利用本文提出的正反演方法對實測數(shù)據(jù)進行了處理。根據(jù)已有地質(zhì)地球物理資料，試驗區(qū)上覆較厚的第四系和第三系沉積層(厚度大于500 m)，巖性主要為黏土、泥質(zhì)粉砂、粉砂、礫石層等，地層電阻率整體偏低。本次測量的目的是探測區(qū)內(nèi)淺表含水層及下部隔水層的深度，并評估其橫向連續(xù)性。試驗工作中，測線總長度為1 200 m，點距為20 m，共布置4個矩形發(fā)射線框，尺寸為200 m×100 m，每個發(fā)射線框控制測線長度為300 m，發(fā)射源-接收點布置如圖8所示。受工作區(qū)茂密的植被影響，發(fā)射線框未能布置成標準的矩形，根據(jù)GPS記錄的航跡繪制了實際的布置形狀。工作儀器為TerraTEM，發(fā)射電流為5～9 A，基頻為5 Hz，觀測垂直感應(yīng)電壓分量，接收探頭有效面積為10 000 m，疊加次數(shù)為256次。

圖8 發(fā)射源-接收點布置圖Fig.8 Transmitter-receiver Layout

區(qū)內(nèi)整體干擾水平較低，實測衰減曲線在大部分計算時窗范圍內(nèi)都較為光滑，僅在晚期出現(xiàn)一定程度的震蕩。圖9給出了框4控制的16個測點的實測歸一化(電流和接收有效面積歸一)信號曲線，其中儀器記錄的關(guān)斷時間約為220 μs。在關(guān)斷期間，同時存在一次場和二次場，但一次場的強度要遠高于二次場，因此，在220 μs之前的時刻，響應(yīng)曲線主要表現(xiàn)的是一次場的特性。從圖9可以看出：當測點位于框內(nèi)時(960～1140點)，未關(guān)斷的一次場為正值；而當測點位于框外時(900～940點、1160～1200點)，未關(guān)斷的一次場為負值；這與圖2、3給出的結(jié)論相吻合。圖10給出了所有測點的多測道曲線。為避免信號符號反轉(zhuǎn)給數(shù)據(jù)分析和處理帶來的不便，將起始時間道(288.5 μs)選在框外測點出現(xiàn)負響應(yīng)的時刻之后。從圖10可以看出：不考慮晚期受干擾信號情況下，4個發(fā)射框控制的測點響應(yīng)整體上較為均勻；但當測點在框外時，響應(yīng)值會出現(xiàn)明顯的降低，如圖中黑色虛線框圈出的位置。

對所有測點的數(shù)據(jù)進行掐頭去尾處理并進行五點圓滑濾波處理，得到光滑衰減的信號曲線，然后利用OCCAM算法對整條測線的數(shù)據(jù)進行一維反演。反演最大深度取600 m，共分40層，迭代次數(shù)為10次。正演部分中，發(fā)射源按照實際布設(shè)形狀進行路徑積分。圖11給出了所有測點的最終擬合殘差。從圖11可以看出：數(shù)據(jù)擬合程度較好，除個別測點因干擾程度較大導致擬合殘差較大外，其他測點的擬合殘差都在5%以內(nèi)。

圖9 框4控制測點的實測歸一化響應(yīng)曲線Fig.9 Measured Normalized Response Curves of Control Measuring Points in Loop 4

圖10 所有測點的多測道響應(yīng)曲線Fig.10 Multi-channel Response Curves of All Measuring Points

圖11 所有測點迭代10次的擬合殘差Fig.11 Fitting Residuals of All Measuring Points with 10 Iterations

圖12 電阻率-深度反演斷面Fig.12 Inversion Section of the Resistivity-depth

圖12給出了整條測線的反演電阻率-深度斷面圖。從圖12可以看出，測區(qū)整體電阻率較低，符合該區(qū)厚新生代沉積層的實際地質(zhì)情況。根據(jù)反演結(jié)果，可以確定區(qū)內(nèi)含水層深度范圍為20～60 m，下伏表現(xiàn)為相對高阻的隔水層(厚度約為60 m)。含水層在橫向上厚度較為均勻，表現(xiàn)出較好的連續(xù)貫通性。隔水層電阻率在橫向上出現(xiàn)一定程度的突變，推測可能與地層的完整性、空隙度、含水性等因素有關(guān)。在測點300處有一個鉆進深度約150 m的鉆孔，鉆孔揭示含水層的頂界面深度為22.5 m，隔水層頂界面深度為56.8 m，隔水層底界面深度為121.5 m，該結(jié)果與反演結(jié)果基本一致，驗證了方法的有效性和準確性。

5 結(jié) 語

本文針對回線瞬變電磁法中的大定源回線裝置開展了理論分析和應(yīng)用研究。通過對矩形回線框內(nèi)外一次場、二次場的計算以及對三層地電模型和野外實測數(shù)據(jù)進行反演，得到以下結(jié)論。

(1)利用電偶極子沿發(fā)射線框進行線性積分，可以準確獲得大定源回線裝置發(fā)射線框內(nèi)外任意點的瞬變電磁場。

(2)不同觀測點空間位置會導致響應(yīng)的幅值和形態(tài)發(fā)生變化，尤其是框外的測點在早期會發(fā)生信號符號反轉(zhuǎn)現(xiàn)象，但隨著時間推移，觀測點空間位置帶來的影響逐漸減弱，晚期二次場的分布趨于均勻。因此，在實際觀測中，可以選擇在框外一定范圍觀測二次場。

(3)在進行一維反演時，考慮發(fā)射線框的形狀和實際的觀測位置等因素，可以很好地反映地層的電性結(jié)構(gòu)，使反演結(jié)果更加可靠。