多道瞬變電磁法共中心點(diǎn)道集數(shù)據(jù)聯(lián)合反演

李海, 薛國(guó)強(qiáng), 鐘華森， 底青云

1 中國(guó)科學(xué)院礦產(chǎn)資源研究重點(diǎn)實(shí)驗(yàn)室，中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 2 中國(guó)科學(xué)院頁(yè)巖氣與地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室，中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 中國(guó)科學(xué)院大學(xué)， 北京 100049

?

多道瞬變電磁法共中心點(diǎn)道集數(shù)據(jù)聯(lián)合反演

李海1,3, 薛國(guó)強(qiáng)1, 鐘華森1,3， 底青云2

1 中國(guó)科學(xué)院礦產(chǎn)資源研究重點(diǎn)實(shí)驗(yàn)室，中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 2 中國(guó)科學(xué)院頁(yè)巖氣與地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室，中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 中國(guó)科學(xué)院大學(xué)， 北京 100049

多道瞬變電磁法是目前地面資源勘查領(lǐng)域的研究熱點(diǎn)，開(kāi)展該方法觀測(cè)數(shù)據(jù)的一維反演對(duì)完成裝備研發(fā)及技術(shù)推廣具有積極的推動(dòng)作用.針對(duì)多道瞬變電磁法特殊的源波形和數(shù)據(jù)采集方式，本文對(duì)多道瞬變電磁數(shù)據(jù)的響應(yīng)提取和一維反演進(jìn)行了研究.首先，采用反卷積方法從偽隨機(jī)二進(jìn)制序列全波形響應(yīng)數(shù)據(jù)中提取大地脈沖響應(yīng).為了避免脈沖響應(yīng)中空氣波的干擾，提出采用反向積分算法，把脈沖響應(yīng)轉(zhuǎn)換成下降沿階躍響應(yīng)，并將其整理成共中點(diǎn)道集.最后，為充分利用各個(gè)偏移距下的響應(yīng)數(shù)據(jù)對(duì)不同深度范圍目標(biāo)體的分辨能力，采用基于OCCAM算法的聯(lián)合反演方法對(duì)共中心點(diǎn)道集數(shù)據(jù)進(jìn)行聯(lián)合反演.結(jié)果表明在不同目標(biāo)體深度以及有噪聲干擾的情況下，共中心點(diǎn)道集數(shù)據(jù)聯(lián)合反演計(jì)算均可獲得較好的結(jié)果.

多道瞬變電磁法； OCCAM反演； PRBS； 聯(lián)合反演

1 引言

瞬變電磁法是一種時(shí)間域人工源電磁法，廣泛應(yīng)用于金屬礦、石油、煤炭等化石能源探測(cè)中(薛國(guó)強(qiáng)等, 2007;Xue et al., 2012a, 2012b).目前淺部資源開(kāi)發(fā)殆盡，探測(cè)深部資源成為趨勢(shì).傳統(tǒng)的瞬變電磁法的探測(cè)深度和探測(cè)精度難以滿(mǎn)足深部資源探測(cè)的要求，發(fā)展大深度探測(cè)的瞬變電磁法探測(cè)新技術(shù)和新裝備迫在眉睫(滕吉文, 2010;薛國(guó)強(qiáng)等, 2013).英國(guó)愛(ài)丁堡大學(xué)的Wright等(2001)提出了多道瞬變電磁法(Multi-channel Transient Electromagnetic method, MTEM)新技術(shù).與傳統(tǒng)瞬變電磁法相比，該方法具有采用大功率接地源發(fā)射、陣列式接收、偽隨機(jī)二進(jìn)制序列(Pseudo-random binary sequence, PRBS)作為激勵(lì)等特點(diǎn)(Hobbs et al., 2006;Ziolkowski et al., 2007a, 2007b;薛國(guó)強(qiáng)等, 2015;Li et al., 2016)，已成功應(yīng)用于油藏監(jiān)控和油氣資源探測(cè)領(lǐng)域(Ziolkowski et al., 2010)，而且探測(cè)深度和探測(cè)精度可望達(dá)到對(duì)深部礦產(chǎn)資源進(jìn)行探測(cè)的要求.

然而，由于與傳統(tǒng)瞬變電磁法存在的眾多不同，多道瞬變電磁法的數(shù)據(jù)處理及解釋與傳統(tǒng)瞬變電磁法存在明顯區(qū)別.一方面，傳統(tǒng)瞬變電磁法在激勵(lì)源的關(guān)斷期間接收純二次場(chǎng)信號(hào)，而多道瞬變電磁法采用偽隨機(jī)二進(jìn)制序列作為激勵(lì)源，采集的是源發(fā)射期間的全波形響應(yīng).為此，傳統(tǒng)的基于純二次場(chǎng)衰減曲線(xiàn)反演方法不能直接應(yīng)用于多道瞬變電磁法數(shù)據(jù)，需要研究適用于PRBS信號(hào)激勵(lì)源數(shù)據(jù)的處理及解釋方法.另一方面，多道瞬變電磁法采用接收陣列進(jìn)行數(shù)據(jù)采集，數(shù)據(jù)量大且對(duì)于測(cè)線(xiàn)上的每個(gè)觀測(cè)點(diǎn)均可獲得多個(gè)偏移距下的多次覆蓋數(shù)據(jù)(Ziolkowski et al., 2007b).由于探測(cè)深度與偏移距、電阻率、時(shí)間等多個(gè)因素相關(guān)，為充分利用不同偏移距數(shù)據(jù)對(duì)地下不同深度目標(biāo)體的覆蓋，需要發(fā)展適用于多道瞬變電磁法多次覆蓋數(shù)據(jù)的反演方法.

基于以上分析，本文對(duì)多道瞬變電磁法共中心點(diǎn)道集數(shù)據(jù)一維反演進(jìn)行了研究.首先，通過(guò)反卷積處理從PRBS激勵(lì)源響應(yīng)中提取直接包含地下介質(zhì)信息的大地脈沖響應(yīng).針對(duì)現(xiàn)有的對(duì)大地脈沖響應(yīng)積分獲得上升沿階躍響應(yīng)存在的技術(shù)難點(diǎn)，本文提出采用反向積分算法從大地脈沖響應(yīng)中進(jìn)一步恢復(fù)出下降沿階躍響應(yīng).最后，對(duì)于所提取的共中心點(diǎn)道集數(shù)據(jù)，采用多偏移距OCCAM聯(lián)合反演，得到地下地質(zhì)的地電信息.反演結(jié)果表明共中心點(diǎn)道集數(shù)據(jù)的反演在不同目標(biāo)體深度以及有噪聲干擾的情況下，均可獲得較好的反演結(jié)果.

2 多道瞬變電磁法原理

MTEM方法的關(guān)鍵技術(shù)在于采用電偶極源發(fā)射源信號(hào)，采用電偶極陣列在源的軸向采集響應(yīng)數(shù)據(jù)(Wright et al., 2002;Ziolkowski et al., 2007b;Zhdanov, 2010).MTEM數(shù)據(jù)采集系統(tǒng)的平面簡(jiǎn)圖如圖1所示.

假設(shè)大地為線(xiàn)性時(shí)不變系統(tǒng)，接收電極對(duì)采集的響應(yīng)信號(hào)可表示為源激勵(lì)信號(hào)與大地脈沖響應(yīng)的卷積的形式：

v(t)=ΔxsΔxriAB(t)*g(t)*r(t)+n(t),

(1)

其中Δxs是源電極距，Δxr是接收電極距，iAB(t)是源電流強(qiáng)度，*表示卷積，g(t)是大地脈沖響應(yīng)，r(t)為接收儀器有關(guān)的系統(tǒng)響應(yīng)，n(t)為不相關(guān)的噪聲信號(hào).與常規(guī)瞬變電磁法在源信號(hào)關(guān)斷期間采集純二次場(chǎng)信號(hào)不同，在多道瞬變電磁法中，采集系統(tǒng)記錄的為源信號(hào)發(fā)射和關(guān)斷期間的全波形信號(hào).通過(guò)在源附近采用與接收系統(tǒng)相同的記錄設(shè)備采集源電流信號(hào)，并根據(jù)式(1)，采用反卷積方法提取與地下介質(zhì)電性結(jié)構(gòu)相關(guān)的大地脈沖響應(yīng).

3 純二次場(chǎng)響應(yīng)提取

由于多道瞬變電磁法采集的為全波形信號(hào)，傳統(tǒng)基于純二次場(chǎng)響應(yīng)的反演方法不能直接應(yīng)用于多道瞬變電磁法數(shù)據(jù)的反演.為此，本文首先采用反卷積的方法從PRBS全波形數(shù)據(jù)中提取大地脈沖響應(yīng)，并提出采用反向積分算法從大地脈沖響應(yīng)恢復(fù)得到下降沿階躍響應(yīng).

圖1 MTEM 數(shù)據(jù)采集裝置示意圖Fig.1 Schematic diagram of source-receiver configuration of MTEM acquisition

3.1 大地脈沖響應(yīng)

由于采用反卷積方法從接收響應(yīng)中提取大地脈沖響應(yīng)，多道瞬變電磁法可采用PRBS信號(hào)作為發(fā)射源信號(hào).PRBS信號(hào)能夠壓制隨機(jī)噪聲，增加有用信號(hào)的頻帶寬度，獲得更好信噪比的響應(yīng).當(dāng)采用PRBS信號(hào)作為激勵(lì)源時(shí)，得到的全波形響應(yīng)與常規(guī)瞬變電磁法所得的衰減曲線(xiàn)不同.圖2b為正演模擬得到的PRBS響應(yīng)信號(hào)，對(duì)應(yīng)的源波形信號(hào)如圖2a所示，正演模型為電阻率為100 Ωm的均勻半空間，偏移距為1000 m.由圖可知，不同于常規(guī)瞬變電磁純二次場(chǎng)響應(yīng)，偽隨機(jī)碼激勵(lì)源響應(yīng)包含源的一次場(chǎng)響應(yīng)和感應(yīng)產(chǎn)生的二次場(chǎng)響應(yīng).因此，須采用反卷積方法提取大地脈沖響應(yīng)，以獲得地下介質(zhì)信息.

PRBS信號(hào)具有二值自相關(guān)特性(Mutagi, 1996; 陳海龍和李宏，2005)，考慮長(zhǎng)度為N的PRBS信號(hào)，其自相關(guān)函數(shù)Rss為

(2)

基于此二值自相關(guān)特性，在最小二乘意義上采用式(1)進(jìn)行反卷積計(jì)算提取大地脈沖g(t)響應(yīng)問(wèn)題可轉(zhuǎn)換為求解如下維納-何甫方程(Pipkin, 1991; Noor, 1993)：

(3)

其中Rvs為源電流與接收響應(yīng)的互相關(guān)函數(shù)，Ts為大地脈沖響應(yīng)的最晚時(shí)間道.通過(guò)源電流和接收電壓的互相關(guān)計(jì)算，接收響應(yīng)中的不相關(guān)噪聲得到了壓制.式(3)可以寫(xiě)成矩陣形式：

(4)

矩陣方程(4)左端的系數(shù)矩陣為對(duì)稱(chēng)矩陣，且沿主對(duì)角線(xiàn)的平行線(xiàn)排列的元素全部相等.因此，這個(gè)N+1階方陣實(shí)際上由N+1個(gè)元素確定，這種方陣被稱(chēng)為托布列茲(Toepliz)陣(Golub and Van Loan, 2012).相應(yīng)的方程組被稱(chēng)為托布列茲方程組，它的解可以采用萊文森遞歸算法進(jìn)行快速求解，得到大地脈沖響應(yīng)(鄒謀炎, 2001).圖2a為采用反饋移位寄存器所得到的7階PRBS信號(hào)，圖2b為正演所得到的電阻率為100 Ωm，偏移距為1000 m時(shí)的偽隨機(jī)碼響應(yīng)，圖2c為采用上述算法從圖2a和圖2b中恢復(fù)得到的大地脈沖響應(yīng).當(dāng)采用不含噪聲數(shù)據(jù)時(shí)，所恢復(fù)的大地脈沖響應(yīng)與真實(shí)大地脈沖響應(yīng)之間的相對(duì)誤差小于1%.

3.2 下降沿階躍響應(yīng)

圖2 理想情況下反卷積提取大地脈沖響應(yīng)(a) 7階PRBS信號(hào)； (b) 7階PRBS信號(hào)激勵(lì)源響應(yīng)； (c) 對(duì)(a)和(b)反卷積得到的大地脈沖響應(yīng).Fig.2 Extracting earth impulse responses by deconvolution(a) PRBS sequence of 7 order; (b) Response exited by 7_order PRBS sequence; (c) Impulse responses from deconvolution of (a) and (b).

(6)

在任意時(shí)刻，上階躍響應(yīng)與下階躍響應(yīng)之和為直流響應(yīng)值F(0)，圖3給出了電阻率為20 Ωm的均勻半空間下，上升沿階躍響應(yīng)與下降沿階躍響應(yīng)之間的關(guān)系.由圖可知，在早期上升沿階躍響應(yīng)與下降沿階躍響應(yīng)具有相等的幅值，均為直流響應(yīng)F(0)的一半.在晚期，上升沿階躍響應(yīng)逐漸趨近于直流響應(yīng)值，而下降沿階躍響應(yīng)逐漸衰減至零.

圖3 上階躍響應(yīng)與下階躍響應(yīng)之間的關(guān)系Fig.3 Relationship between step-up and step-down responses

(7)

其中tmax為積分起點(diǎn)，即最晚時(shí)間道，M為對(duì)應(yīng)的時(shí)間道序號(hào).采用該算法計(jì)算下降沿階躍響應(yīng)無(wú)需直流響應(yīng)值，同時(shí)能避開(kāi)早期空氣波的干擾.為了驗(yàn)證該算法，對(duì)電阻率分別為20 Ωm和100 Ωm的均勻半空間，偏移距為1000 m時(shí)的大地脈沖響應(yīng)進(jìn)行積分，并將積分結(jié)果與直接采用解析解計(jì)算的下降沿階躍響應(yīng)進(jìn)行對(duì)比.如圖4所示，反向積分獲得的衰減曲線(xiàn)與解析解之間的相對(duì)誤差很小，幾乎可忽略不計(jì).

圖4 反向積分獲得的階躍響應(yīng)與解析解的對(duì)比Fig.4 Comparison of analytical solution and step-off response from backward integration

4 共中心點(diǎn)道集

由于多道瞬變電磁法采用接收陣列采集響應(yīng)數(shù)據(jù)，對(duì)于每個(gè)記錄點(diǎn)均可獲得不同偏移距的多次覆蓋數(shù)據(jù).如圖5a所示，當(dāng)在發(fā)射點(diǎn)A處向大地注入源信號(hào)，在接收點(diǎn)A采集響應(yīng)信號(hào)時(shí)，記錄點(diǎn)A為發(fā)射點(diǎn)與接收點(diǎn)的中心位置(Ziolkowski et al., 2007b).記錄點(diǎn)位置與偏移距分別采用公式(8)和公式(9)進(jìn)行計(jì)算.

(8)

圖5 (a)陣列式數(shù)據(jù)采集系統(tǒng)及(b)共中心點(diǎn)道集示意圖Fig.5 Schematic diagrams of data collection array (a) and resulting common middle-point gather (b)

(9)

其中xs和xr分別為發(fā)射點(diǎn)和接收點(diǎn)的坐標(biāo)，xc為觀測(cè)點(diǎn)坐標(biāo)，r為偏移距.由于采用接收陣列采集數(shù)據(jù)，在某一發(fā)射點(diǎn)處發(fā)射時(shí)，可以得到不同記錄點(diǎn)、不同偏移距的數(shù)據(jù).對(duì)于圖5a所示數(shù)據(jù)采集系統(tǒng)，在發(fā)射點(diǎn)A所采集到的數(shù)據(jù)為圖中實(shí)心圓點(diǎn)所示.當(dāng)完成發(fā)射點(diǎn)A處的數(shù)據(jù)采集后，整個(gè)發(fā)射-接收陣列將沿測(cè)線(xiàn)向前平移以完成此條測(cè)線(xiàn)的數(shù)據(jù)采集.因此，對(duì)于任一記錄點(diǎn)，如記錄點(diǎn)B，將獲得不同偏移距下的響應(yīng)數(shù)據(jù).將該記錄點(diǎn)各個(gè)偏移距下的數(shù)據(jù)整理為一個(gè)數(shù)據(jù)體，稱(chēng)之為共中心點(diǎn)道集.共中心點(diǎn)道集內(nèi)的響應(yīng)數(shù)據(jù)的個(gè)數(shù)取決于觀測(cè)系統(tǒng)的布設(shè)方式，對(duì)于圖5a中的數(shù)據(jù)采集系統(tǒng)，在記錄點(diǎn)B可以獲得r1，r2，r3，r4四個(gè)偏移距下的數(shù)據(jù)，如圖中虛線(xiàn)所示.圖5b給出了100 Ωm均勻半空間下，1000 m, 1200 m, 1400 m, 1600 m四個(gè)偏移距下的下降沿階躍響應(yīng)所組成的共中心點(diǎn)道集示意圖.

5 聯(lián)合反演算法

由于對(duì)共中心點(diǎn)道集進(jìn)行反演，在每個(gè)記錄點(diǎn)均有多個(gè)偏移距下的數(shù)據(jù)，需要進(jìn)行聯(lián)合反演以充分利用各個(gè)偏移距對(duì)不同深度的覆蓋.反演的基本算法采用OCCAM算法，通過(guò)在目標(biāo)函數(shù)中引入模型粗糙度的正則化項(xiàng)，最終得到滿(mǎn)足目標(biāo)擬合殘差的最光滑模型.OCCAM算法的目標(biāo)函數(shù)U為

(10)

其中m=(m1,m2,…,mN)為模型向量，d=(d1,d2,…,dM)為數(shù)據(jù)向量，F(xiàn)為正演算子，*為目標(biāo)擬合殘差，?為粗糙度矩陣，用于計(jì)算模型的粗糙度，為數(shù)據(jù)誤差加權(quán)矩陣.μ為拉格朗日乘子，用于平衡模型粗糙度和擬合殘差對(duì)目標(biāo)函數(shù)的影響.

OCCAM方法是一種將非線(xiàn)性問(wèn)題線(xiàn)性化的反演算法.根據(jù)泰勒定理，采用局部線(xiàn)性化的思想，有

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

(11)

其中，mk+Δm=mk+1，J(mk)為雅可比矩陣，

(12)

通過(guò)如式(11)所示近似，正則化最小二乘問(wèn)題(10)可以通過(guò)(13)式進(jìn)行迭代求解，

(13)

對(duì)于多道瞬變電磁法共中心點(diǎn)道集數(shù)據(jù)，本文采用多個(gè)偏移距下的響應(yīng)進(jìn)行聯(lián)合反演，共同產(chǎn)生一個(gè)地電模型.假設(shè)共中心點(diǎn)道集內(nèi)，不同偏移距響應(yīng)數(shù)量為k，則在反演過(guò)程中，數(shù)據(jù)向量、正演算子、雅可比矩陣均應(yīng)進(jìn)行拓展，

Wd=diag(1/δ1,1/δ2,…,1/δM×k),

(14)

(15)

其中數(shù)據(jù)向量d和正演函數(shù)F的維度為M×k，加權(quán)矩陣W的維度為(M×k)×(M×k)，雅克比矩陣J的維度為(M×k)×N.

在多道瞬變電磁法的觀測(cè)方式中，探測(cè)深度與偏移距呈正相關(guān)(Ziolkowski et al., 2007b;Li et al., 2015)，同時(shí)與電阻率相關(guān).另外，由于多道瞬變電磁法最終得到的是多個(gè)時(shí)間道的電磁響應(yīng)，隨著電磁場(chǎng)的擴(kuò)散，時(shí)間與深度之間本身存在對(duì)應(yīng)關(guān)系.因此，為了充分利用共中心點(diǎn)道集內(nèi)不同偏移距對(duì)不同深度范圍的覆蓋，提高多道瞬變電磁法在各個(gè)深度范圍的分辨率，我們根據(jù)共中心點(diǎn)道集內(nèi)偏移距不同，對(duì)反演中的靈敏度矩陣進(jìn)行加權(quán)處理.

構(gòu)造一個(gè)靈敏度加權(quán)矩陣Wm，對(duì)不同偏移距范圍的靈敏度函數(shù)進(jìn)行加權(quán).在多道瞬變電磁法中，對(duì)不同深度范圍內(nèi)的目標(biāo)體存在一個(gè)最佳探測(cè)偏移距.對(duì)于深度為d處的目標(biāo)體，其最佳探測(cè)偏移距r約為2dr/2，對(duì)應(yīng)權(quán)值為50%.最后，共中心點(diǎn)道集內(nèi)偏移距最大的響應(yīng)所對(duì)應(yīng)的靈敏度函數(shù)不加權(quán).最終，得到的靈敏度矩陣為：

Jw=WdJWm.

(16)

將得到的靈敏度矩陣代入迭代方程(13)，通過(guò)一維線(xiàn)性搜索算法選擇μ.在擬合殘差達(dá)到目標(biāo)擬合殘差后，引入模型粗糙度計(jì)算，最終得到滿(mǎn)足目標(biāo)擬合殘差的最光滑模型.

6 模擬數(shù)據(jù)反演算例

根據(jù)以上算法，采用Fortran語(yǔ)言編制了反演程序.反演的基本模型為一維層狀模型，采用如前所述算法提取下降沿階躍響應(yīng).瞬變電磁響應(yīng)的時(shí)間道數(shù)為51道，采用對(duì)數(shù)等間距的時(shí)間序列.在OCCAM反演中，由于反演模型層數(shù)較多，采用厚度相等的層狀模型進(jìn)行反演是一種常見(jiàn)策略(Key, 2009; Sudha et al., 2014).另外，在多道瞬變電磁法中，模型分辨率不僅與深度有關(guān)，且與數(shù)據(jù)點(diǎn)的偏移距也存在較大的關(guān)系.為此，在本文的研究中，采用厚度相等的層狀模型進(jìn)行反演，反演模型的層數(shù)為40層，每層的厚度為50 m，反演參數(shù)為電阻率.反演所采用的數(shù)據(jù)為整理的共中心點(diǎn)道集，同時(shí)對(duì)下降沿階躍響應(yīng)單支曲線(xiàn)進(jìn)行反演以作為對(duì)比.反演過(guò)程中的擬合殘差是反演結(jié)果模型的響應(yīng)與觀測(cè)數(shù)據(jù)(此處為真實(shí)模型的響應(yīng))之間的相對(duì)誤差，計(jì)算公式為：

(17)

其中di為真實(shí)模型的響應(yīng)，fi(ms)為反演結(jié)果模型的響應(yīng)數(shù)據(jù)，N為數(shù)據(jù)點(diǎn)的個(gè)數(shù)，σi為每個(gè)數(shù)據(jù)點(diǎn)的標(biāo)準(zhǔn)偏差.在本文的數(shù)值試驗(yàn)及多道瞬變電磁法實(shí)測(cè)數(shù)據(jù)采集中，均沒(méi)有對(duì)每個(gè)數(shù)據(jù)點(diǎn)的誤差進(jìn)行估算.因此，本文參考了Ziolkowski等(2007b)的處理方式，將每個(gè)時(shí)間道的標(biāo)準(zhǔn)偏差σi設(shè)置為其幅值的1%.

6.1 低阻薄層模型

為了對(duì)比單支衰減曲線(xiàn)與共中心點(diǎn)道集聯(lián)合反演的反演效果，首先采用低阻薄層模型進(jìn)行反演計(jì)算.反演的初始模型為電阻率為1 Ωm的均勻半空間，H型薄層模型的低阻層位于500～600 m深度，電阻率為5 Ωm，背景電阻率為100 Ωm，目標(biāo)擬合殘差為1.反演的目標(biāo)為恢復(fù)深度為500 m處的低阻層的存在，為此采用偏移距為目標(biāo)體深度2倍左右的數(shù)據(jù)，分別采用了800 m、1000 m、1200 m偏移距下的單支曲線(xiàn)，同時(shí)采用這三個(gè)偏移距的數(shù)據(jù)組成的共中心點(diǎn)道集進(jìn)行了反演.圖6a表示共中心點(diǎn)道集數(shù)據(jù)，圖6b為共中心點(diǎn)道集內(nèi)偏移距1000 m響應(yīng)的擬合情況.反演結(jié)果如圖6c所示，對(duì)于H型模型，偏移距為800 m, 1000 m和1200 m的單支曲線(xiàn)的反演結(jié)果均可以反映低阻目標(biāo)層的存在，但是與共中心點(diǎn)道集的反演結(jié)果對(duì)比，單支曲線(xiàn)的反演結(jié)果的目標(biāo)層的電阻率吻合更好.另外，由于聯(lián)合反演算法中減弱了偏移距較小的響應(yīng)對(duì)反演深部目標(biāo)體的影響，單支曲線(xiàn)對(duì)目標(biāo)層以下的背景層的電阻率反演效果較差，而共中心點(diǎn)道集的反演效果較好.6.2 五層典型模型

為了更進(jìn)一步對(duì)比共中心點(diǎn)道集聯(lián)合反演與單支曲線(xiàn)的反演對(duì)由淺及深整個(gè)深度范圍內(nèi)的目標(biāo)層的反演效果.設(shè)計(jì)如下五層模型.五層模型的電阻率分別為100 Ωm，70 Ωm，30 Ωm，60 Ωm和20 Ωm，各層的深度范圍分別為0～200 m, 200～500 m, 500～700 m，700～1000 m，1000～底層.模型的電阻率滿(mǎn)足ρ1>ρ2>ρ3<ρ4>ρ5，在這個(gè)模型中Q型模型(ρ1>ρ2>ρ3)，K型模型(ρ3<ρ4>ρ5)，H型模型ρ2>ρ3<ρ4等三層模型均得到了體現(xiàn).該五層模型包含了三種典型的三層模型，能在一定程度上反映反演方法對(duì)各種模型和各個(gè)深度范圍的反演效果.與前述反演中目標(biāo)層深度為500 m情形不同，此處反演的目標(biāo)為識(shí)別深度為200 m，500 m，700 m，1000 m收發(fā)距的電性分界面.為了提高反演對(duì)700 m和1000 m深度的電性層的分辨能力，采用了1600 m偏移距處的數(shù)據(jù).反演共采用了800 m、1000 m、1600 m偏移距下的單支曲線(xiàn)，同時(shí)采用這三個(gè)偏移距組成的共中心點(diǎn)道集進(jìn)行了反演，反演的目標(biāo)擬合殘差為1.

圖6 低阻薄層模型反演結(jié)果(a) 共中心點(diǎn)道集數(shù)據(jù)； (b) 共中心點(diǎn)道集內(nèi)1000 m偏移距數(shù)據(jù)擬合情況； (c) 反演結(jié)果.Fig.6 Inversion results of model H with low-resistivity thin layer(a) Common middle point gather data; (b) Fitted offset data at 1000 m; (c) Inversion results.

圖7 典型五層模型反演結(jié)果Fig.7 Inversion result for typical 5-layer model

圖7給出了反演結(jié)果，從圖中可以看出，800 m偏移距單支曲線(xiàn)的反演結(jié)果較好反映出了淺部地層，但是深部電阻率偏離真實(shí)模型.1000 m和1600 m偏移距下的反演結(jié)果無(wú)法反映出淺部高阻層，對(duì)第三層和第四層的電阻率的反映效果好，對(duì)于較深部區(qū)域的低阻基底，1600 m偏移距下的反演結(jié)果要好于800 m和1000 m偏移距下的反演結(jié)果.由此可見(jiàn)，單支曲線(xiàn)的反演結(jié)果無(wú)法單獨(dú)反演出所有五層模型，然而共中心點(diǎn)道集的聯(lián)合反演結(jié)果基本給出了各層的電阻率和深度信息，在不同目標(biāo)體深度情況下，均可獲得較好的探測(cè)效果.因此，采用多道瞬變電磁法觀測(cè)方式獲得單個(gè)測(cè)點(diǎn)的共中心點(diǎn)道集進(jìn)行聯(lián)合反演解釋?zhuān)梢垣@得較好的探測(cè)效果.

6.3 反演算法穩(wěn)定性檢驗(yàn)

為了研究反演算法的穩(wěn)定性，在模擬數(shù)據(jù)中添加噪聲后再進(jìn)行反演.首先，在模擬數(shù)據(jù)中添加均值為零、標(biāo)準(zhǔn)偏差為3%的高斯噪聲.采用上述反卷積算法和反向積分算法獲得下降沿階躍響應(yīng).采用Munkholm和Auken (1996)提出的噪聲模型，對(duì)得到的下降沿階躍響應(yīng)采用對(duì)數(shù)等間距積分采樣，并對(duì)采樣得到的數(shù)據(jù)按時(shí)間道進(jìn)行多次疊加.因此，生成噪聲數(shù)據(jù)隨時(shí)間衰減，符合瞬變電磁響應(yīng)曲線(xiàn)的衰減特性.由于在模擬數(shù)據(jù)中添加了噪聲，當(dāng)擬合殘差過(guò)小時(shí)，迭代過(guò)程開(kāi)始擬合所添加的噪聲.為此，將反演的目標(biāo)擬合殘差設(shè)為2，其他反演參數(shù)與圖7中的反演參數(shù)一致.

圖8 對(duì)數(shù)據(jù)添加噪聲后典型模型反演結(jié)果Fig.8 Inversion result of data added with noisy on typical 5-layer model

圖8給出了對(duì)添加噪聲后的共中心點(diǎn)道集數(shù)據(jù)的反演結(jié)果.圖中紅色曲線(xiàn)表示擬合殘差為2時(shí)的反演結(jié)果，灰色曲線(xiàn)簇表示當(dāng)擬合殘差低于2時(shí)得到的模型.當(dāng)擬合噪聲的反演結(jié)果擬合殘差小于2時(shí)，有兩支偏離深部地層的電阻率，其余的反演結(jié)果均能較好地反映深部的電阻率.這反映了采用不同偏移距組合的共中心點(diǎn)道集數(shù)據(jù)進(jìn)行反演具備一定的抗干擾能力.從圖中可以看到，共中心點(diǎn)道集反演結(jié)果淺部地層受到噪聲的影響較大，對(duì)深部地層的電阻率的反演效果較好.這是由于所采用的三個(gè)偏移距范圍內(nèi)，在擬合殘差小于2的模型空間中的模型較多且淺層電阻率差異較大.因此，當(dāng)目標(biāo)擬合殘差為2時(shí)(紅色曲線(xiàn)所示)，所得到的淺部反演結(jié)果較差.當(dāng)目標(biāo)擬合殘差小于2時(shí)(灰色曲線(xiàn)所示)，淺部的反演結(jié)果得到改善.因此對(duì)于多道瞬變電磁法觀測(cè)方式而言.為提高整個(gè)深度范圍內(nèi)對(duì)目標(biāo)體的分辨能力，需要獲得多個(gè)偏移距下信噪比較高的數(shù)據(jù).

7 結(jié)論

為了加速多道瞬變電磁法裝備的研發(fā)進(jìn)程，本文對(duì)多道瞬變電磁法數(shù)據(jù)進(jìn)行了一維反演研究.首先，針對(duì)多道瞬變電磁法所采集的偽隨機(jī)碼激勵(lì)源全波形數(shù)據(jù)，采用反卷積方法提取大地脈沖響應(yīng).由于上升沿階躍響應(yīng)早期幅值難以獲取以及大地脈沖響應(yīng)早期的空氣波的存在，對(duì)脈沖響應(yīng)的直接積分難以獲得準(zhǔn)確幅值的上升沿階躍響應(yīng).為此，通過(guò)分析上升沿階躍響應(yīng)與下降沿階躍響應(yīng)的對(duì)偶關(guān)系，提出采用反向積分算法恢復(fù)下降沿階躍響應(yīng).通過(guò)計(jì)算表明，反向積分恢復(fù)的下降沿階躍響應(yīng)能夠很好地與解析解對(duì)應(yīng)，該方法為一種簡(jiǎn)單有效的提取下降沿階躍響應(yīng)的積分算法.

其次，針對(duì)多道瞬變電磁法陣列式采集方式所得到的多次覆蓋數(shù)據(jù)，本文的聯(lián)合反演方法能夠更充分地利用不同偏移距對(duì)不同深度的覆蓋.單支曲線(xiàn)反演結(jié)果和共中心點(diǎn)道集聯(lián)合反演結(jié)果的對(duì)比表明，所采用的聯(lián)合反演策略在不同目標(biāo)體深度情況下，均可獲得更好的結(jié)果.最后，對(duì)多道瞬變電磁法數(shù)據(jù)進(jìn)行加噪后進(jìn)行反演，反演結(jié)果能夠恢復(fù)出模型特征.本文的研究表明，在不同目標(biāo)體深度及有噪聲干擾的情況下，共中心點(diǎn)道集數(shù)據(jù)聯(lián)合反演計(jì)算能獲得較好的結(jié)果.

Chen H L, Li H. 2005. Generation and analysis of PN sequence based on MATLAB.ComputerSimulation(in Chinese), 22(5): 98-100.Connell D. 2011. A comparison of marine time-domain and frequency-domain controlled source electromagnetic methods [Master′s thesis]. San Diego: University of California.Golub G H, Van Loan C F. 2012. Matrix Computations. Vol.3. Baltimore: JHU Press. Hobbs B, Ziolkowski A, Wright D. 2006. Multi-Transient Electromagnetics (MTEM)-controlled source equipment for subsurface resistivity investigation.∥ 18th EM Induction Workshop. El Vendrell, Spain. Key K. 2009. 1D inversion of multicomponent, multifrequency marine CSEM data: Methodology and synthetic studies for resolving thin resistive layers.Geophysics, 74(2): F9-F20.

Li H, Xue G Q, Zhou N N, et al. 2015. Appraisal of an array tem method in detecting a mined-out area beneath a conductive layer.PureandAppliedGeophysics, 172(10): 2917-2929.

Li H, Xue G Q, Zhao P, et al. 2016. The Hilbert-Huang transform-based Denoising Method for the TEM response of a PRBS source signal.PureandAppliedGeophysics, 173(8): 2777-2789. Munkholm M S, Auken E. 1996. Electromagnetic noise contamination on transient electromagnetic soundings in culturally disturbed environments.JournalofEnvironmentalandEngineeringGeophysics, 1(2): 119-127. Mutagi R N. 1996. Pseudo noise sequences for engineers.Electronics&CommunicationEngineeringJournal, 8(2): 79-87. Noor M A. 1993. Wiener-Hopf equations and variational inequalities.JournalofOptimizationTheoryandApplications, 79(1): 197-206.Pipkin A C. 1991. Wiener-Hopf equations.∥ A Course on Integral Equations. New York: Springer, 160-172.

Sudha, Tezkan B, Siemon B. 2014. Appraisal of a new 1D weighted joint inversion of ground based and helicopter-borne electromagnetic data.GeophysicalProspecting, 62(3): 597-614.Teng J W. 2010. Strengthening exploration of metallic minerals in the second depth space of the crust, accelerating development and industralization of new geophysical technology and instrumental equipment.ProgressinGeophys. (in Chinese), 25(3): 729-748, doi: 10.3969/j.issn.1004-2903.2010.03.001.

Wright D, Ziolkowski A, Hobbs B. 2002. Hydrocarbon detection and monitoring with a multicomponent transient electromagnetic

(MTEM) survey.TheLeadingEdge, 21(9): 852-864.Wright D A, Ziolkowski A, Hobbs B A. 2001. Hydrocarbon detection with a multi-channel transient electromagnetic survey.∥ SEG Technical Program Expanded Abstracts. SEG, 1435-1438.

Xue G Q, Li X, Di Q Y. 2007. The progress of TEM in theory and application.ProgressinGeophys. (in Chinese), 22(4): 1195-1200.

Xue G Q, Bai C Y, Yan S, et al. 2012a. Deep sounding TEM investigation method based on a modified fixed central-loop system.JournalofAppliedGeophysics, 76: 23-32.

Xue G Q, Qin K Z, Li X, et al. 2012b. Discovery of a large-scale porphyry molybdenum deposit in Tibet through a modified TEM exploration method.JournalofEnvironmental&EngineeringGeophysics, 17: 19-25.

Xue G Q, Chen W Y, Zhou N N, et al. 2013. Short-offset TEM technique with a grounded wire source for deep sounding.ChineseJ.Geophys. (in Chinese), 56(1): 255-261, doi: 10.6038/cjg20130126.

Xue G Q, Yan S, Di Q Y, et al. 2015. Technical analysis of multitransient electromagnetic method.JournalofEarthSciencesandEnvironment(in Chinese), 37(1): 94-100.

Zhdanov M S. 2010. Electromagnetic geophysics: Notes from the past and the road ahead.Geophysics, 75(5): 75A49-75A66.

Ziolkowski A, Carson R, Wright D. 2007a. New technology to acquire, process, and interpret transient EM data.∥ EGM 2007 International Workshop. Italy.Ziolkowski A, Hobbs B A, Wright D. 2007b. Multitransient electromagnetic demonstration survey in France.Geophysics, 72(4): F197-F209.

Ziolkowski A, Parr R, Wright D, et al. 2010. Multi-transient electromagnetic repeatability experiment over the North Sea Harding field.GeophysicalProspecting, 58(6): 1159-1176.

Zou M Y.2001. Deconvolution and Signal Recovery (in Chinese). Beijing: National Defence Industry Press, 143-144.

附中文參考文獻(xiàn)

陳海龍, 李宏. 2005. 基于MATLAB的偽隨機(jī)序列的產(chǎn)生和分析. 計(jì)算機(jī)仿真, 22(5): 98-100.

滕吉文. 2010. 強(qiáng)化第二深度空間金屬礦產(chǎn)資源探查, 加速發(fā)展地球物理勘探新技術(shù)與儀器設(shè)備的研制及產(chǎn)業(yè)化. 地球物理學(xué)進(jìn)展, 25(3): 729-748, doi: 10.3969/j.issn.1004-2903.2010.03.001.

薛國(guó)強(qiáng), 李貅, 底青云. 2007. 瞬變電磁法理論與應(yīng)用研究進(jìn)展. 地球物理學(xué)進(jìn)展, 22(4): 1195-1200.薛國(guó)強(qiáng), 陳衛(wèi)營(yíng), 周楠楠等. 2013. 接地源瞬變電磁短偏移深部探測(cè)技術(shù). 地球物理學(xué)報(bào), 56(1): 255-261, doi: 10.6038/cjg20130126.

薛國(guó)強(qiáng), 閆述, 底青云等. 2015. 多道瞬變電磁法(MTEM)技術(shù)分析. 地球科學(xué)與環(huán)境學(xué)報(bào), 37(1): 94-100.

鄒謀炎. 2001. 反卷積和信號(hào)復(fù)原. 北京: 國(guó)防工業(yè)出版社, 143-144.

(本文編輯 何燕)

Joint inversion of CMP gathers of multi-channel transient electromagnetic data

LI Hai1,3, XUE Guo-Qiang1, ZHONG Hua-Sen1,3, DI Qing-Yun2

1KeyLaboratoryofMineralResources,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2KeyLaboratoryofShaleGasandGeoengineering,InstituteofGeologyandGeophysics,ChineseAcademyofSciences，Beijing100029,China3UniversityofChineseAcademyofSciences,Beijing100049,China

The multi-channel transient electromagnetic (MTEM) method is one of the research focuses of prospecting for underground mineral resources. The study of the 1D inversion of MTEM data will promote the development of the instrument and theory of the methods. Based on the pseudo-random binary sequence (PRBS) source waveforms and configurations of the MTEM method, this paper presents the processing and inversion methods of MTEM data. Firstly, the earth impulse response of PRBS is obtained by deconvolution. Then step-off responses recovered from the impulse responses by reversed integration are sorted into common mid-point (CMP) gathers. Finally, the joint inversion scheme based on OCCAM inversion is tested on the CMP gather data. The results show that the inversion of CMP gathers can obtain good results at different depths with presence of noise.

Multi-channel transient electromagnetic method; OCCAM inversion; PRBS; Joint inversion

10.6038/cjg20161206.

國(guó)家重大科研裝備研制項(xiàng)目“深部資源探測(cè)核心裝備研發(fā)”(ZDYZ2012-1)-05子項(xiàng)目“多通道大功率電法勘探儀”-04課題“M-TEM資料處理及偏移成像軟件研制”資助.

李海，男，1988年生，博士，主要從事電磁法探測(cè)理論與應(yīng)用研究. E-mail: thinda@163.com

10.6038/cjg20161206

P631

2015-10-30，2016-07-22收修定稿

李海, 薛國(guó)強(qiáng), 鐘華森等. 2016. 多道瞬變電磁法共中心點(diǎn)道集數(shù)據(jù)聯(lián)合反演. 地球物理學(xué)報(bào)，59(12):4439-4447,

Li H, Xue G Q, Zhong H S, et al. 2016. Joint inversion of CMP gathers of multi-channel transient electromagnetic data.ChineseJ.Geophys. (in Chinese),59(12):4439-4447,doi:10.6038/cjg20161206.