邢濤，袁偉，李建慧

(1.北京探創(chuàng)資源科技有限公司，北京 100071； 2.內(nèi)蒙古地質(zhì)工程有限責(zé)任公司，內(nèi)蒙古 呼和浩特 010010； 3.中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢 430074)

0 引言

瞬變電磁法一維反演始于20世紀(jì)80年代。如:Raiche等基于阻尼最小二乘法對瞬變電磁法重疊回線和電阻率法施倫貝謝裝置采集的數(shù)據(jù)開展了聯(lián)合反演[9]，黃皓平、王維中采用阻尼最小二乘法實現(xiàn)了時間域航空電磁數(shù)據(jù)的一維反演[10]；上述反演算法均采用了層狀介質(zhì)參數(shù)化反演，即同時反演層厚度和電阻率，反演結(jié)果依賴于模型層數(shù)和初始模型[11]。Farquharson和Oldenburg在目標(biāo)函數(shù)中加入待求模型與某個參考模型之間的偏離程度表征項，并設(shè)定反演模型層數(shù)多于觀測數(shù)據(jù)個數(shù)，且各層介質(zhì)層厚固定，只反演電阻率[11]。隨后，丹麥奧胡斯大學(xué)Auken課題組采用阻尼最小二乘法也開展了瞬變電磁法一維反演[12]，并引入橫向約束項用于解決反演斷面圖中電阻率和層厚橫向連續(xù)性差、薄層分辨率低等問題[13-14],引入空間約束項用于解決反演深度切片中層界面不光滑等問題[15]。德國科隆大學(xué)Tezkan課題組采用基于阻尼最小二乘法和Occam算法[16]的反演程序，開展回線源瞬變電磁法、長偏移距瞬變電磁法、射頻大地電磁法等方法的單一反演或聯(lián)合反演[17]。

國內(nèi)課題組也開展了類似研究。如：毛立峰等采用自適應(yīng)正則化方法開展了直升機(jī)航空瞬變電磁法數(shù)據(jù)一維反演[18]，齊彥福等采用Occam算法對m序列發(fā)射波形多道瞬變電磁法數(shù)據(jù)開展了反演研究[19]，類似工作也見于李海等[20]。除了上述線性化優(yōu)化方法外，模擬退火法[6]、粒子群優(yōu)化——阻尼最小二乘混合算法[21]等智能優(yōu)化算法也被成功用于瞬變電磁法一維反演。

綜上所述，瞬變電磁法一維反演發(fā)展比較完善成熟。本研究基于Constable等開發(fā)的Occam反演算法[16]和Key課題組開發(fā)的Dipole1D一維正演程序[22]，開展了回線源瞬變電磁法一維反演研究。

1 一維正演框架

Dipole1D適用于水平電偶源和垂直電偶源，并能通過方位角和傾角控制電偶源形態(tài)，且發(fā)射源和接收點可放置于層狀介質(zhì)任意層位。本節(jié)中只敘述電磁法勘探的一維正演框架，不再給出公式，具體公式詳見上述文獻(xiàn)。

圖1為回線源瞬變電磁法一維正演框架，該框架從水平電偶源和垂直電偶源激發(fā)的電磁場空間域、頻率域麥克斯韋方程組出發(fā)。需要說明的是，此處的電偶源并非點源，而是具有一定長度的接地線源，當(dāng)收發(fā)距遠(yuǎn)大于線源尺度時，其可視為電偶源。一個方形回線可簡單視為由4個首尾相連的接地線源組成，當(dāng)回線尺寸比較小時，每個接地線源還可進(jìn)一步細(xì)分。一個復(fù)雜形態(tài)發(fā)射回線可視為一系列首尾相連的接地線源，這些電偶源的方位角和傾角不盡相同。Dipole1D程序處理這種復(fù)雜形態(tài)場源時，如果該場源方位角與x和y軸斜交，需采用x和y兩個方向布設(shè)的電偶源組合求解；如果該場源傾角不為零，還需引入z方向布設(shè)的電偶源來共同求解。

圖1 回線源瞬變電磁法一維正演框架Fig.1 The framework of 1D forward modeling for loop-source transient electromagnetic methods

借助勢表示電場E和磁場H是求解麥克斯韋方程組的有效途徑之一。引入矢量勢和標(biāo)量勢后，由麥克斯韋方程組推導(dǎo)出關(guān)于矢量勢的雙旋度方程，再結(jié)合一些限定條件，比如洛倫茲規(guī)范，可得到關(guān)于矢量勢在空間域的亥姆霍茲方程。對于一維介質(zhì)，介質(zhì)分界面在水平坐標(biāo)x和y無限延伸，并與垂直坐標(biāo)面(等z平面)重合，可利用二維傅里葉變換將亥姆霍茲方程從空間域轉(zhuǎn)換至波數(shù)域；伴隨著這種轉(zhuǎn)換，將關(guān)于x、y和z的偏微分方程轉(zhuǎn)變?yōu)橐子谇蠼?、關(guān)于z的常微分方程。結(jié)合電磁場衰減的邊界條件和在界面處電磁場切向分量的連續(xù)性，可以唯一求解上述關(guān)于矢量勢的微分方程，獲得波數(shù)域矢量勢；再利用二維傅里葉逆變換得到空間域矢量勢，并根據(jù)勢與場的轉(zhuǎn)換關(guān)系，獲得頻率域電磁場。這一轉(zhuǎn)換過程中，需數(shù)值求解多個關(guān)于發(fā)射源長度和波數(shù)的二重積分。對關(guān)于發(fā)射源長度的積分，通常采用高斯—勒讓德積分求解，該方法可通過增大積分點數(shù)來提高求解精度。關(guān)于波數(shù)的積分含有零階或一階形式的第一類貝塞爾函數(shù)，該類積分通常采用漢克爾變換求解。關(guān)于漢克爾變換及常用的變換系數(shù)，可參考Guptasarma和Singh[23]等已發(fā)表的文獻(xiàn)。

逐個計算并累加每個接地線源激發(fā)的頻率域磁場，即可獲得回線源激發(fā)的磁場，再采用正弦變換或余弦變換求得磁場階躍響應(yīng)或脈沖響應(yīng)，進(jìn)而得到感應(yīng)電動勢。本研究采用Anderson研發(fā)的787個系數(shù)的正弦變換和余弦變換[24]。

課堂教學(xué)主要分為兩個層次，第一個是理解性層次，即對新知識的基礎(chǔ)性掌握，這是每個學(xué)生都需要達(dá)到的層次。不得不引起重視的是，學(xué)生對新知識的掌握程度是不一致的，這將反應(yīng)在他們的應(yīng)用能力上。第二個是操練鞏固層次，基于第一個層次帶來的差距，操練鞏固階段的分層應(yīng)兼顧橫向分層和縱向分層，橫向分層包括聽、說、讀、寫等不同題型的訓(xùn)練，縱向分層則是在每類題型中設(shè)置不同難度層次的題目以適應(yīng)A、B、C類學(xué)生的不同需要。習(xí)題的設(shè)置是對教師教學(xué)水平的一大考驗，要避免兩個極端。一是無用題，即對所有學(xué)生來說答案一目了然的題目；二是超綱題，一旦題目難度超出了A類生的發(fā)展區(qū)水平，在課堂上缺乏受眾，也是毫無意義的。

2 Occam反演方法

采用經(jīng)典的Occam算法[16]開展一維反演計算。該算法的基本思想是：通過在目標(biāo)函數(shù)中引入模型粗糙度的正則化項，來獲得能夠擬合數(shù)據(jù)的最光滑模型。設(shè)定反演的目標(biāo)函數(shù)為

U=‖?m‖2+‖P(m-m*)‖2+

(1)

為了使得目標(biāo)函數(shù)最小化，通常求取它關(guān)于模型參數(shù)m的一階導(dǎo)數(shù)，并令其為零。由于正演函數(shù)F(m)是模型參數(shù)m的非線性函數(shù)，對其在某個模型mk處進(jìn)行泰勒展開，有：

F(mk+Δm)≈F(mk)+JkΔm,

(2)

式中Jk為正演函數(shù)對模型參數(shù)的一階導(dǎo)數(shù)(稱為雅可比矩陣)：

(3)

矩陣元素為

(4)

Dipole 1D程序在正演時，可采用解析法一并計算靈敏度矩陣。

Occam反演的模型更新公式與傳統(tǒng)的高斯牛頓法(Gauss-Newton，GN)相同，即：

(5)

在本研究中，使用均方根擬合差(RMS misfit)來衡量數(shù)據(jù)的擬合程度，其表達(dá)式為：

(6)

式中：M為數(shù)據(jù)數(shù)目，sm為第m個數(shù)據(jù)對應(yīng)的標(biāo)準(zhǔn)差，[dm-Fm(mk+1(μ))]/sm為單個數(shù)據(jù)的擬合差。

3 理論模型

3.1 層狀介質(zhì)模型

采用中心回線裝置，發(fā)射源為邊長50 m的方形回線。該層狀介質(zhì)共含4層，第一至第四層介質(zhì)厚度依次為50、50、100 m和∞，電阻率ρ依次為200、20、100 Ω·m和50 Ω·m。觀測時間序列共含31個時刻，其范圍在10-5～10-2s。待反演感應(yīng)電動勢由Dipole 1D計算，并添加了5%的隨機(jī)噪音。Occam反演中，地電模型(含空氣層)共含52層，其中地下51層介質(zhì)電阻率待反演，層厚度均為10 m。反演時電阻率范圍為100～103Ω·m，初始模型為均勻半空間(101.5Ω·m)，無參考模型。反演中，最大迭代次數(shù)為40，目標(biāo)擬合差為3。對于該模型，迭代7次后，均方根擬合差小于3，迭代終止。

如圖2所示，對于該模型，反演結(jié)果能夠較好地反映真實地電結(jié)構(gòu)。圖3中，由反演結(jié)果計算的感應(yīng)電動勢與待反演數(shù)據(jù)吻合良好，二者擬合差絕對值多數(shù)在3以內(nèi)，只有個別大于3。通過這一模型，驗證了程序的正確性。

圖2 對于四層模型，回線中心點反演地電模型Fig.2 The inversion geo-electric model for the loop-center point and for 4-layer model

圖3 四層模型的反演結(jié)果Fig.3 The data computed from the inversion model

3.2 傾斜界面模型

如圖4a所示，發(fā)射回線鋪設(shè)在傾角為20°的斜坡之上。該模型表層介質(zhì)厚度為30 m，電阻率為50 Ω·m，沿走向方向長度為200 m；表層下伏介質(zhì)為100 Ω·m的半空間。發(fā)射回線為10 m×10 m的正方形。圖4a中，發(fā)射回線Tx1、Tx2和Tx3的中心位置分別為(-93.969, 0.0,-34.302)、(0.0, 0.0,0.0)和(93.969, 0.0, 34.302)。采用中心回線裝置，觀測量為傾斜界面法向方向的感應(yīng)電動勢。采用課題組開發(fā)的時域矢量有限單元法[25]計算磁感應(yīng)強(qiáng)度脈沖響應(yīng)三分量，并將其合成觀測量。觀測時間序列范圍為0.025～1.995 ms，共20個時刻。

圖4 傾斜界面模型示意Fig.4 The sketch map for the tilted-surface model

將該傾斜界面模型逆時針旋轉(zhuǎn)20°可得到其等效模型(圖4b)。該模型中，發(fā)射回線Tx1、Tx2和Tx3均鋪設(shè)于水平界面之上，其中心位置分別為(-100,0.0,0.0) m、(0.0,0.0,0.0) m和(100,0.0,0.0) m，觀測量為感應(yīng)電動勢垂直分量。接下來對實際模型和等效模型開展一維反演。Occam反演中，地電模型(含空氣層)共含26層，其中地下25層介質(zhì)電阻率待反演，層厚度均為5 m。反演時電阻率范圍為100～104Ω·m，初始模型為均勻半空間(101.5Ω·m)，無參考模型，這種設(shè)置我們稱之為反演方案1。反演中，最大迭代次數(shù)為40，目標(biāo)擬合差為3?？紤]到此處反演主要用于驗證上述兩種模型反演結(jié)果的等價性，我們對待反演的感應(yīng)電動勢只添加了1%的隨機(jī)噪音。

在發(fā)射回線Tx1和Tx3激發(fā)情況下，反演結(jié)果如圖5所示。圖中可見實際模型和等效模型的反演結(jié)果具有良好的一致性，這說明了對于鋪設(shè)于傾斜界面的發(fā)射回線模型，如果觀測量是界面法向方向的感應(yīng)電動勢，那么可按水平界面情形開展一維反演。值得注意的是，盡管反演結(jié)果中深部結(jié)果與真實模型吻合較好，但淺層結(jié)果與真實模型存在較大差異，如反演結(jié)果中低阻層主要集中于地下15～25 m，其電阻率遠(yuǎn)小于50 Ω·m，低阻層的上下介質(zhì)均為高阻層，地表電阻率甚至超過1 000 Ω·m。

圖5 回線中心測點數(shù)據(jù)反演獲取的地電模型(反演方案1)Fig.5 Geoelectric model obtained by inversion of data at the center of the loop( inversion schemes 1)

為了使反演結(jié)果與真實模型一致，反演算法中考慮了參考模型以及結(jié)構(gòu)約束，即反演方案2和反演方案3。方案2中，將模型中地下第一、第二層介質(zhì)電阻率設(shè)置為50 Ω·m；方案3在方案2的基礎(chǔ)上，允許電阻率結(jié)構(gòu)在地下30 m處不連續(xù)。如圖6a所示，以發(fā)射回線Tx2激發(fā)時情形為例，采用3種反演方案獲取的地電模型結(jié)構(gòu)基本一致，即使考慮參考模型和結(jié)構(gòu)約束也未能改變低阻層上下介質(zhì)均為高阻層的這一現(xiàn)象。再次對相應(yīng)的一維模型(即表層低阻層在水平方向足夠延伸)開展反演計算(圖6b)，結(jié)果顯示，真實模型為一維情形時反演獲取的地電模型結(jié)構(gòu)與三維情形的結(jié)果基本一致。

圖6 采用不同反演方案時，實際模型(Tx2發(fā)射回線)和一維情形的反演地電模型Fig.6 The inversion geo-electric models obtained from different inversion schemes

圖7是在反演數(shù)據(jù)中添加5%的隨機(jī)誤差并經(jīng)過10次反演計算的結(jié)果。由于每一次反演添加的誤差并不一樣，反演結(jié)果也略有差異，但地電模型結(jié)構(gòu)基本一致，這也驗證了反演算法的穩(wěn)定性。

黑色曲線為一維模型，彩色曲線為反演結(jié)果The black line denotes the realistic model, the color lines denote the inversion models圖7 隨機(jī)噪聲影響下，實際模型(Tx2發(fā)射回線)的反演地電模型Fig.7 Under the influence of random noise, the inversion geoelectric model of the actual model (Tx2 loop)

4 實例

內(nèi)蒙古自治區(qū)那仁寶力格煤田分布區(qū)內(nèi)地勢平坦，新近系噴出巖發(fā)育。從火山口和火山錐的存在和分布形態(tài)看，以中心式噴發(fā)為主，裂隙式噴發(fā)次之；巖性組合主要有橄欖拉斑玄武巖和玄武巖?？碧絽^(qū)內(nèi)，玄武巖體電阻率可達(dá)上千歐姆米，而砂巖和泥巖地層電阻率僅為幾十歐姆米。這種明顯的電性差異也是采用瞬變電磁法圈定玄武巖在沉積巖中分布范圍的物性基礎(chǔ)。本次野外工作布設(shè)了11條測線，線距100 m，點距40 m。儀器采用加拿大Phoenix公司的V8多功能電法儀，觀測時間序列為5 Hz的30個時間窗口。發(fā)射回線邊長為600 m×300 m，其中600 m邊長沿測線方向布置。

選取測線1作為本次研究的實例。該測線橫穿噴出于地表的玄武巖露頭，其長度為2 600 m，共66個測點，測點x坐標(biāo)范圍為680～3 280 m。Occam反演中，地電模型(含空氣層)共含52層，其中地下51層介質(zhì)電阻率待反演，層厚度均為20 m。反演時電阻率范圍為100.5～104Ω·m，反演數(shù)據(jù)為感應(yīng)電動勢，初始模型為均勻半空間(101.5Ω·m)，無參考模型。圖8為該工區(qū)5個測點的一維反演電阻率曲線，其曲線形態(tài)大致為電阻率隨深度的增加而先減小再增大。測點x坐標(biāo)由1 000 m增大至3 000 m時，反演結(jié)果中電阻率變化趨勢是先增大后減小，其中淺層幅值變化最為明顯。

圖8 5個測點數(shù)據(jù)分別反演獲取的地電模型Fig.8 The inversion geo-electric models for five survey points

圖9a、c中，由反演結(jié)果計算的感應(yīng)電動勢與實測數(shù)據(jù)曲線形態(tài)基本一致。當(dāng)x=1 000 m時，由反演結(jié)果計算的感應(yīng)電動勢與實測數(shù)據(jù)吻合程度良好(圖9a)；x=2 000 m時，兩類數(shù)據(jù)在8 ms之后吻合程度較低(圖9c)。這一現(xiàn)象與其晚期時刻實測數(shù)據(jù)的標(biāo)準(zhǔn)差(standard deviation，SD)偏大有關(guān)。如圖9d所示，該測點標(biāo)準(zhǔn)差與實測數(shù)據(jù)比值在8 ms之后快速上升，最大值接近28%。另外，如圖9b、d所示，由反演結(jié)果計算的感應(yīng)電動勢與實測數(shù)據(jù)之間的擬合差絕對值在大多數(shù)時刻都小于3，僅在少數(shù)時刻偏大。

圖10是由單點反演結(jié)果拼接的測線1地電模型斷面。該測線共經(jīng)過5個鉆孔，圖中數(shù)字為鉆孔揭露的玄武巖厚度。斷面圖中，上部高阻區(qū)域(藍(lán)色和綠色)為玄武巖分布范圍，下部低阻區(qū)域(紅色)為沉積巖分布范圍，二者分界面連續(xù)、清晰，并呈現(xiàn)出噴出巖典型的“鍋底狀”分布形態(tài)，但未見火山通道。對比反演結(jié)果與鉆孔資料，發(fā)現(xiàn)兩種方法獲取的玄武巖厚度在鉆孔ZK1、ZK2、ZK4和ZK5處基本一致，但在鉆孔ZK3處差異明顯。由此推斷，鉆孔ZK3可能位于火山通道分布范圍之內(nèi)。進(jìn)一步推斷，由于火山通道隨著深度增加，分布范圍逐漸縮小，導(dǎo)致一維反演結(jié)果與真實情況偏差較大。

a、c—不同測點的感應(yīng)電動勢；b、d—擬合差與實測數(shù)據(jù)標(biāo)準(zhǔn)差a、c—EMF curves at different measuring points； b、d—Misfit curver and the standard deviation curves of measured data圖9 x=1 000 m和2 000 m時的一維反演結(jié)果Fig.9 The data computed from the inversion model for x=1 000 m and 2 000 m

圖10 由單點反演結(jié)果拼接的測線1地電模型斷面Fig.10 The section view stitched from the single-point 1D inversion models for survey line 1

5 結(jié)論

本研究實現(xiàn)了基于Occam算法的回線源瞬變電磁法一維反演，并采用理論模型和野外實例驗證了程序的正確性，為瞬變電磁法資料處理解釋提供了有力工具，也為后續(xù)反演研究考慮更多影響因素奠定了基礎(chǔ)。