回線源瞬變電磁法的一維Occam反演

邢濤，袁偉，李建慧

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

0 引言

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

國(guó)內(nèi)課題組也開展了類似研究。如：毛立峰等采用自適應(yīng)正則化方法開展了直升機(jī)航空瞬變電磁法數(shù)據(jù)一維反演[18]，齊彥福等采用Occam算法對(duì)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適用于水平電偶源和垂直電偶源，并能通過(guò)方位角和傾角控制電偶源形態(tài)，且發(fā)射源和接收點(diǎn)可放置于層狀介質(zhì)任意層位。本節(jié)中只敘述電磁法勘探的一維正演框架，不再給出公式，具體公式詳見上述文獻(xiàn)。

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

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

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

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

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

2 Occam反演方法

采用經(jīng)典的Occam算法[16]開展一維反演計(jì)算。該算法的基本思想是：通過(guò)在目標(biāo)函數(shù)中引入模型粗糙度的正則化項(xiàng)，來(lái)獲得能夠擬合數(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ù)，對(duì)其在某個(gè)模型mk處進(jìn)行泰勒展開，有：

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

(2)

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

(3)

矩陣元素為

(4)

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

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

(5)

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

(6)

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

3 理論模型

3.1 層狀介質(zhì)模型

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

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

圖2 對(duì)于四層模型，回線中心點(diǎn)反演地電模型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，沿走向方向長(zhǎng)度為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)。采用中心回線裝置，觀測(cè)量為傾斜界面法向方向的感應(yīng)電動(dòng)勢(shì)。采用課題組開發(fā)的時(shí)域矢量有限單元法[25]計(jì)算磁感應(yīng)強(qiáng)度脈沖響應(yīng)三分量，并將其合成觀測(cè)量。觀測(cè)時(shí)間序列范圍為0.025～1.995 ms，共20個(gè)時(shí)刻。

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

將該傾斜界面模型逆時(shí)針旋轉(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，觀測(cè)量為感應(yīng)電動(dòng)勢(shì)垂直分量。接下來(lái)對(duì)實(shí)際模型和等效模型開展一維反演。Occam反演中，地電模型(含空氣層)共含26層，其中地下25層介質(zhì)電阻率待反演，層厚度均為5 m。反演時(shí)電阻率范圍為100～104Ω·m，初始模型為均勻半空間(101.5Ω·m)，無(wú)參考模型，這種設(shè)置我們稱之為反演方案1。反演中，最大迭代次數(shù)為40，目標(biāo)擬合差為3?？紤]到此處反演主要用于驗(yàn)證上述兩種模型反演結(jié)果的等價(jià)性，我們對(duì)待反演的感應(yīng)電動(dòng)勢(shì)只添加了1%的隨機(jī)噪音。

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

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

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

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

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

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

4 實(shí)例

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

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

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

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

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

a、c—不同測(cè)點(diǎn)的感應(yīng)電動(dòng)勢(shì)；b、d—擬合差與實(shí)測(cè)數(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時(shí)的一維反演結(jié)果Fig.9 The data computed from the inversion model for x=1 000 m and 2 000 m

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

5 結(jié)論

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