齊彥福, 殷長(zhǎng)春*， 王若, 蔡晶

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 130021 2 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京 100029

?

多通道瞬變電磁m序列全時(shí)正演模擬與反演

齊彥福1, 殷長(zhǎng)春1*， 王若2, 蔡晶1

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 130021 2 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京 100029

傳統(tǒng)瞬變電磁方法主要用于金屬礦勘查，無(wú)法滿足油氣資源高阻目標(biāo)體的勘探需要.多通道瞬變電磁(MTEM)系統(tǒng)的出現(xiàn)解決了這一問(wèn)題.該方法采用偽隨機(jī)序列發(fā)射波形和擬地震觀測(cè)方式，測(cè)量同線電場(chǎng)分量，記錄全時(shí)發(fā)射電流及多道觀測(cè)數(shù)據(jù)，實(shí)現(xiàn)對(duì)高阻薄層的高精度探測(cè).鑒于國(guó)內(nèi)對(duì)此方法的研究還處于理論探索階段，尚未進(jìn)行相應(yīng)的仿真模擬和數(shù)據(jù)處理工作，本文針對(duì)m序列發(fā)射波形多通道瞬變電磁法的全時(shí)正演模擬和反演解釋進(jìn)行研究，為國(guó)內(nèi)正在進(jìn)行的MTEM儀器系統(tǒng)研發(fā)及數(shù)據(jù)解釋提供理論指導(dǎo).我們利用方波響應(yīng)移位疊加和電流導(dǎo)數(shù)與階躍響應(yīng)褶積兩種方法實(shí)現(xiàn)理論m序列和實(shí)際發(fā)射波形的全時(shí)正演模擬；再通過(guò)相關(guān)辨識(shí)技術(shù)，削弱噪聲影響，計(jì)算脈沖響應(yīng)；最后對(duì)積分得到的階躍響應(yīng)進(jìn)行共中心點(diǎn)道集數(shù)據(jù)聯(lián)合反演，獲取地下電性分布信息.

多通道瞬變電磁； m序列； 全時(shí)正演； 反演

The MTEM system adopts quasi-seismic observation, measuring in-line electric field component of pseudo-random binary sequence (PRBS) waveforms. Full-wave transmitting current and multi-channel observation data are recorded, so that the resistive thin layers can be detected. We perform full-wave forward modeling for a theoretical m-sequence and actual transmitter waveforms based on superposition of square wave responses and convolution of step responses with derivative of transmitting current, respectively. Taking advantages of the correlation identification technique, the impulse response is calculated while the noise is suppressed. Finally, we invert the CMP data by the Occam′s algorithm to obtain underground electrical information.

Firstly, Fourier transform is used to analyze the spectral components of square waveforms, 5-order 2n-sequence and 5-order m-sequence. It shows that the m-sequence has a larger frequency bandwidth with an equal interval between spectral lines than another two waveforms. Secondly, in order to verify the full-wave forward modelling results, a homogenous half-space of 30 Ωm is set up with a suit of 4-order m-sequence transmitting current. The results obtained by superposition of square wave responses, convolution of step responses with derivate of transmitting current and convolution of impulse response with transmitting current are compared. It shows that the first two results are consistent. We further add 30 dB Gaussian noise to the responses. The relative error after correlation between modeling response and transmitting waveforms is reduced to 1/4 of the original one. The step response obtained by integrating the impulse response is very close to the theoretical value with relative error of just 0.08%. Finally, we design two three-layer models to demonstrate the capability of the MTEM method to resolve shallow and deep hydrocarbon reservoirs. In the first model, a reservoir of 500 Ωm and 50 m thickness is set at the depth of 300 m in a homogenous half-space of 50 Ωm. The survey is done respectively at offsets of 900, 1000, 1100, 1200 and 1300 m. For the second model, we assume the hydrocarbon reservoir is moved to 2000 m depth, the survey is done at 5000, 5200, 5500, 5700 and 6000 m offsets, respectively. The inversion results for the two models by the Occam′s method show that the MTEM method can well resolve both shallow and deep hydrocarbon reservoirs. The CMP data have better resolution than one-offset data.

The MTEM method using m-sequence as transmitting waveforms, measuring multi-channel in-line electric fields, turns out to have a good resolution to hydrocarbon reservoirs. Spectral analysis shows that the m-sequence has a wide frequency band with an equal interval between spectrum lines. The techniques based on superposition of square-wave responses and convolution of step responses with derivate of transmitting current are proved to be very effective for modeling theoretical m-sequence and actual full-wave responses. Using correlation identification technique, the influence from transmitting waveforms is removed in the calculation of impulse response while the noise is suppressed. The Occam′s inversion of CMP data obtained by integration of impulse response can well resolve underground resistive targets like reservoirs.Keywords Multi-transient electromagnetic (MTEM); m-sequence; Full-wave forward modeling; inversion

1 引言

近年來(lái)，電磁方法在資源勘探中獲得越來(lái)越廣泛的應(yīng)用.瞬變電磁法作為一種重要的電磁勘探手段越來(lái)越受到重視.然而，由于電磁干擾的影響以及目標(biāo)體埋深的不斷加大，對(duì)瞬變電磁法勘探能力的要求不斷提高，新方法、新技術(shù)的研究成為解決這一問(wèn)題的有效途徑(薛國(guó)強(qiáng)等，2013).

考慮到發(fā)射波形對(duì)瞬變電磁法勘探能力的影響，前人對(duì)傳統(tǒng)的方波、梯形波、半正弦波和三角波等發(fā)射波形進(jìn)行了深入研究(Liu，1998；陳曙東等，2012；關(guān)珊珊等，2012)，以期在原有裝置和觀測(cè)系統(tǒng)基礎(chǔ)上，通過(guò)調(diào)整發(fā)射波形及相關(guān)參數(shù)達(dá)到改善勘探效果的目的.為獲得更加豐富的頻率信息，何繼善(1998)提出了偽隨機(jī)脈沖電磁觀測(cè)方法和技術(shù).在儀器設(shè)備、裝置型式和數(shù)據(jù)處理手段上前人也做了許多改進(jìn).陳曉東等(2012)研制出適合瞬變電磁法的高溫超導(dǎo)磁強(qiáng)計(jì)，利用其直接測(cè)量磁場(chǎng)分量，與傳統(tǒng)感應(yīng)線圈傳感器相比有效地提高了探測(cè)深度.李貅等(2012)將合成孔徑應(yīng)用到瞬變電磁數(shù)據(jù)處理當(dāng)中，利用相關(guān)疊加技術(shù)，實(shí)現(xiàn)多孔徑數(shù)據(jù)合成，將單點(diǎn)處理方式發(fā)展成為逐點(diǎn)推移多次覆蓋的處理方法，大大提高了瞬變電磁法的分辨率.薛國(guó)強(qiáng)等(2013)提出短偏移距瞬變電磁技術(shù)(SOTEM).隨著偏移距的縮短，接地導(dǎo)線源的場(chǎng)對(duì)地層的反映更加靈敏，通過(guò)增加觀測(cè)時(shí)長(zhǎng)，加大探測(cè)深度，使該方法成為大深度、高分辨探測(cè)地下礦產(chǎn)資源的一種新技術(shù)手段.

Wright 等(2002)設(shè)計(jì)出一種全新的電磁觀測(cè)系統(tǒng)——多通道瞬變電磁系統(tǒng)(MTEM)，此方法利用偽隨機(jī)序列作為發(fā)射波形，測(cè)量同線電場(chǎng)分量，同時(shí)觀測(cè)多道數(shù)據(jù).Wright將其成功應(yīng)用于油氣資源勘查和油田開(kāi)采狀態(tài)監(jiān)測(cè)，大大拓寬了瞬變電磁法的應(yīng)用領(lǐng)域.Ziolkowski 等(2007)利用MTEM裝置對(duì)法國(guó)西南部天然氣儲(chǔ)層進(jìn)行測(cè)量，獲得了理想的勘探效果.Wright 和Ziolkowski(2007)分析了MTEM數(shù)據(jù)噪聲的類型及相應(yīng)壓制方法.Ziolkowski 等(2010)再次對(duì)北海Harding油氣田進(jìn)行了成功勘探.Ziolkowski等(2011)在挪威Peon 油氣田利用瞬變電磁方法進(jìn)行淺?？碧?，并將方波和偽隨機(jī)兩種信號(hào)源進(jìn)行了對(duì)比，結(jié)果顯示偽隨機(jī)信號(hào)擁有更強(qiáng)的分辨能力.

m序列偽隨機(jī)碼作為目前應(yīng)用最廣泛的一種偽隨機(jī)系列已被應(yīng)用于通信領(lǐng)域，如擴(kuò)頻通信、衛(wèi)星通信的碼分多址、數(shù)字?jǐn)?shù)據(jù)中的加密等領(lǐng)域(林智慧等，2009).鑒于m序列具有諸多優(yōu)良的性質(zhì)，逐漸受到地球物理工作者的關(guān)注.湯井田和羅維斌(2008)對(duì)基于相關(guān)辨識(shí)的逆重復(fù)m序列偽隨機(jī)電磁法的原理進(jìn)行分析，討論參數(shù)選取對(duì)壓制噪聲、提高信噪比的影響.逆重復(fù)m序列是由m序列隔位取反得到的，它們具有十分相近的性質(zhì)，對(duì)m序列的研究具有重要的參考價(jià)值.Li 等(2012)對(duì)相關(guān)識(shí)別技術(shù)在時(shí)間域譜激電中的應(yīng)用做了更加深入的研究.淳少恒等(2014)對(duì)偽隨機(jī)m序列及其在電法勘探中的應(yīng)用現(xiàn)狀做了比較全面的介紹.

與國(guó)外日趨成熟的技術(shù)相比，國(guó)內(nèi)對(duì)m序列瞬變電磁法的研究還停留在理論探索階段，并沒(méi)有進(jìn)行仿真模擬和數(shù)據(jù)處理研究(淳少恒等，2014).本文在前人已有工作基礎(chǔ)上，將現(xiàn)有技術(shù)集成應(yīng)用到一個(gè)新的領(lǐng)域——多通道瞬變電磁法，并給出了針對(duì)MTEM裝置m序列發(fā)射波形正演模擬及反演解釋方法和完整技術(shù)流程，為國(guó)內(nèi)正在進(jìn)行的MTEM系統(tǒng)儀器研究及其后續(xù)數(shù)據(jù)解釋工作提供理論指導(dǎo).我們首先簡(jiǎn)要介紹m序列的產(chǎn)生及其重要性質(zhì)，然后著重闡述m序列響應(yīng)的全時(shí)正演模擬以及脈沖辨識(shí)方法，最后對(duì)積分獲得的階躍響應(yīng)進(jìn)行Occam反演，獲得地下電性分布信息.本文的研究成果必將對(duì)國(guó)內(nèi)m序列多通道瞬變電磁勘查系統(tǒng)研發(fā)起到積極的促進(jìn)作用.

2 MTEM裝置及測(cè)量方式

根據(jù)Ziolkowski等(2007)，MTEM裝置采用電偶極發(fā)射源，通過(guò)電極AB向地下發(fā)射偽隨機(jī)序列(PRBS)，在同線方向觀測(cè)電場(chǎng)ex分量(如圖1).在整個(gè)工作過(guò)程中，觀測(cè)并記錄全部時(shí)刻的發(fā)射電流和電磁響應(yīng).

圖1 MTEM裝置示意圖Fig.1 Schematic diagram of MTEM

在實(shí)際工作中，一般連續(xù)發(fā)射多個(gè)周期的偽隨機(jī)序列，以增大信號(hào)強(qiáng)度.對(duì)于同一套偽隨機(jī)編碼序列進(jìn)行重復(fù)發(fā)射，可以實(shí)現(xiàn)觀測(cè)數(shù)據(jù)的多次疊加，從而壓制噪聲.通過(guò)移動(dòng)發(fā)射源和接收機(jī)位置，可以獲得共中心點(diǎn)道集多道觀測(cè)數(shù)據(jù).對(duì)其同時(shí)反演處理，可以有效地提高對(duì)高阻儲(chǔ)集層的識(shí)別能力.

3 m序列

3.1 產(chǎn)生機(jī)理

m序列是由線性反饋移位寄存器(圖2)生成的最長(zhǎng)線性反饋移位寄存器序列，其周期為2n-1，n為移位寄存器的階數(shù).通過(guò)調(diào)整移位寄存器起始狀態(tài)和反饋線鏈接狀態(tài)，可以實(shí)現(xiàn)對(duì)m序列進(jìn)行編碼(吳明捷和杜天蒼，2002).

圖2 線性反饋移位寄存器示意圖Fig.2 Schematic diagram of linear feedback shift register

圖2為線性反饋移位寄存器示意圖，其中a為移位寄存器的狀態(tài)，C為反饋系數(shù)，C=0表示斷開(kāi)，而C=1表示鏈接.碼元的新?tīng)顟B(tài)由n階串接寄存器線性抽頭經(jīng)模2加法反饋獲得(淳少恒等，2014).

3.2 m序列特性

m序列具有很多優(yōu)良性質(zhì)：(1)均衡性：在一個(gè)周期內(nèi)“1”元素出現(xiàn)的次數(shù)比“0”元素出現(xiàn)的次數(shù)多1；(2)游程特性：一個(gè)周期內(nèi)共有2n-1個(gè)游程，且“0”元素的游程和“1”元素的游程各占一半，長(zhǎng)度為k(1≤k≤n-2)的游程數(shù)占游程總數(shù)的1/2k(林智慧等， 2009)；(3)優(yōu)良的自相關(guān)特性：m序列具有尖銳的自相關(guān)函數(shù)，其表達(dá)形式為(謝雪康和楊曉蓉，2010)

(1)式中N是m序列碼長(zhǎng)，t0是碼元寬度，k為周期數(shù)(k=0,1,2,…).由(1)式可以看出碼長(zhǎng)越長(zhǎng)且碼元寬度越小，m序列的自相關(guān)函數(shù)越接近脈沖函數(shù).

3.3 頻譜成分分析

在傳統(tǒng)電磁法中，一般發(fā)射方波作為激勵(lì)源.這種波形在電子元件上容易實(shí)現(xiàn)，然而其頻率成分比較單一，能量主要集中在基頻及其奇次諧波上(Ziolkowski et al.，2011)，且隨頻率增加能量迅速衰減，基波及3次諧波的功率占方波總功率的90%以上(湯井田和羅維斌，2008).為了豐富發(fā)射波形的頻譜成分，何繼善(1998)提出了偽隨機(jī)觀測(cè)方案，設(shè)計(jì)并發(fā)射2n系列偽隨機(jī)信號(hào)(何繼善等，2009).2n系列的能量主要集中在主頻2k(k=0，1，2，…，n-1)上，且主頻上的幅值基本相等，隨著n的增大頻帶變寬(湯井田和羅維斌，2008).與前兩種發(fā)射波形不同，m序列偽隨機(jī)碼具有較寬的頻帶，譜線之間是等間隔的，包含更多的頻率信息，而且其包含直流成分，這是前兩種波形所沒(méi)有的(湯井田等，2007).

將m序列作為發(fā)射電流波形，影響其勘探能力的主要參數(shù)包括m序列的階數(shù)、碼元寬度、重復(fù)發(fā)射周期數(shù)以及電流強(qiáng)度(湯井田和羅維斌，2008；淳少恒等，2014).碼元寬度和階數(shù)影響頻帶寬度及能量分布.碼元寬度越大，階數(shù)越低，頻帶就越窄，能量主要集中在低頻部分，勘探深度越大，信噪比也越高.然而由于高頻成分的缺失必將造成分辨率降低.碼元寬度越小，階數(shù)越高，頻帶就越寬，高頻成分增加，分辨率將得到提高.然而由于能量被均勻分配到整個(gè)頻帶范圍內(nèi)，勢(shì)必造成單個(gè)頻率成分能量降低，降低信噪比.電流強(qiáng)度越大，重復(fù)發(fā)射次數(shù)越多，壓制噪聲的能力就越強(qiáng)，信噪比越高.隨之對(duì)發(fā)射機(jī)的硬件要求就越高，工作效率降低.在實(shí)際工作當(dāng)中，應(yīng)綜合考慮信噪比、分辨率和勘探深度，根據(jù)勘探目標(biāo)合理調(diào)整m序列的相關(guān)參數(shù)，才能達(dá)到理想的勘探效果.

圖3 發(fā)射波形及頻譜特征 (a1) 方波波形；(b1) 5階2n序列；(c1) 5階m序列； (a2)方波振幅譜；(b2) 5階2n序列振幅譜；(c2) 5階m序列振幅譜.Fig.3 Transmitted waveform and spectrum (a1) Square wave; (b1) 5 series 2n-sequence; (c1) 5 series m-sequence; (a2) Amplitude spectrum for square wave; (b2) Amplitude spectrum for 5 series 2n-sequence; (c2) Amplitude spectrum for 5 series m-sequence.

4 m序列全時(shí)正演理論

4.1 理論m序列電磁響應(yīng)模擬

在數(shù)字信號(hào)處理中，方波可以利用兩個(gè)階躍波反向移位疊加合成(圖4).參照此原理，在忽略二次電流耦合的條件下，可采用階躍響應(yīng)反向移位疊加計(jì)算方波響應(yīng)(圖5).

圖4 階躍波合成方波示意圖Fig.4 Schematic diagram of square-wave obtained from superposed step-waves

圖5 階躍響應(yīng)合成方波響應(yīng)Fig.5 Square-wave response obtained from step-waves

圖5中帶有圓圈符號(hào)的灰色實(shí)線為第一個(gè)正向階躍波激發(fā)的電磁響應(yīng)，帶有三角符號(hào)的灰色實(shí)線是第二個(gè)負(fù)向階躍波激發(fā)的電磁響應(yīng)，黑色實(shí)線是兩個(gè)階躍響應(yīng)疊加合成的方波響應(yīng).

m序列實(shí)際上是由一系列方波組成的，所以m序列可以利用方波移位疊加進(jìn)行合成(圖6).與階躍響應(yīng)合成方波響應(yīng)的方法類似，可采用方波響應(yīng)合成m序列響應(yīng).

圖6 方波合成m序列示意圖Fig.6 Schematic diagram of m-sequence obtained from square-waves

4.2 實(shí)際發(fā)射波形電磁響應(yīng)模擬

在實(shí)際工作中，發(fā)射電流在供電瞬間并不能迅速達(dá)到穩(wěn)定，斷電時(shí)出現(xiàn)延遲，發(fā)射波形實(shí)際上為一系列梯形波.為模擬實(shí)際發(fā)射波形的電磁響應(yīng)，可以采用褶積方法進(jìn)行計(jì)算.Smith 和Lee(2004)指出：在航空系統(tǒng)中，當(dāng)t→0時(shí)半空間的脈沖響應(yīng)趨近于t-1/2，出現(xiàn)奇異性，而階躍響應(yīng)趨近一個(gè)常數(shù).針對(duì)這一問(wèn)題，殷長(zhǎng)春等(2013)利用階躍響應(yīng)與電流的時(shí)間導(dǎo)數(shù)褶積代替脈沖響應(yīng)與電流褶積的方法成功實(shí)現(xiàn)了任意波形全時(shí)響應(yīng)的計(jì)算.研究發(fā)現(xiàn)，在地面電磁法中，t→0時(shí)電場(chǎng)脈沖響應(yīng)同樣出現(xiàn)奇異性，所以本文采用階躍響應(yīng)與電流導(dǎo)數(shù)褶積的方法計(jì)算實(shí)際發(fā)射波形電磁響應(yīng)，即

(2)

其中I(t)是發(fā)射電流，而Bs(t)是下階躍響應(yīng).(2)式中只考慮了由于電流變化感應(yīng)出的電磁場(chǎng)，并沒(méi)有考慮電流自身產(chǎn)生的直流場(chǎng)，為計(jì)算任意發(fā)射波形全時(shí)響應(yīng)的總場(chǎng)，需在(2)式中加入直流項(xiàng)，

(3)

式中fDC是單位電流強(qiáng)度的直流場(chǎng).

4.3 結(jié)果對(duì)比

理論上，當(dāng)發(fā)射波形供電穩(wěn)定和斷電延遲時(shí)間足夠短時(shí)，梯形波可以近似方波，利用移位疊加和褶積兩種方法計(jì)算的m序列電磁響應(yīng)具有可比性.本文設(shè)計(jì)一個(gè)30 Ωm的均勻半空間模型，電偶源采用單位長(zhǎng)度，收發(fā)距為1000 m，電流強(qiáng)度為30 A，發(fā)射波形為4階m序列(圖7)，正負(fù)電平轉(zhuǎn)換時(shí)間為40.96 μs.

圖7 4階m序列電流波形Fig.7 Transmitting waveform of 4-order m-sequence

圖8 移位疊加與褶積結(jié)果對(duì)比圖Fig.8 Comparison of shift-superposition with convolution results

圖8為電場(chǎng)ex響應(yīng)，其中黑色實(shí)線為方波響應(yīng)移位疊加計(jì)算結(jié)果，十字星為階躍響應(yīng)與電流導(dǎo)數(shù)褶積計(jì)算結(jié)果，而灰色實(shí)線為脈沖響應(yīng)與電流波形褶積結(jié)果.可以看出，三種方法在off-time會(huì)得到一致的結(jié)果，然而在on-time，由于脈沖響應(yīng)在t→0時(shí)出現(xiàn)奇異性，使得結(jié)果出現(xiàn)明顯偏差.而前兩組結(jié)果吻合較好，均方相對(duì)誤差小于3%.

5 m序列電磁響應(yīng)反演

由于在電磁反演當(dāng)中需要多次重復(fù)正演運(yùn)算，

如果直接對(duì)含有發(fā)射波形的電磁響應(yīng)進(jìn)行反演，正演計(jì)算量過(guò)大，嚴(yán)重影響反演效率.因此本文利用反卷積方法，消除發(fā)射波形的影響，提高反演效率.

5.1 脈沖辨識(shí)

根據(jù)Ziolkowski等(2007)，m序列的電磁響應(yīng)可以寫為

(4)

式中I(t)是發(fā)射電流，BI(t)是脈沖響應(yīng)，rCD(t)是接收機(jī)響應(yīng)，n(t)是噪聲.通常情況下，將發(fā)射電流和接收機(jī)響應(yīng)的褶積稱為儀器系統(tǒng)響應(yīng)，可將(4)式化簡(jiǎn)為

(5)

式中s(t)是儀器系統(tǒng)響應(yīng)，具有與發(fā)射電流相同的性質(zhì).根據(jù)Wiener-Hopf方法(Wiener, 1949)，將(5)式兩端分別與儀器系統(tǒng)響應(yīng)做互相關(guān)，可以得到

(6)

其中Rss(t)是s(t)的自相關(guān)矩陣.由于m序列具有優(yōu)良相關(guān)特性，噪聲是隨機(jī)的，所以n(t)?s(t)≈0，(6)式可化簡(jiǎn)為

(7)

式中B′(t)=B(t)?s(t).當(dāng)Rss(t)為脈沖函數(shù)時(shí)可以得到

(8)

但實(shí)際工作中Rss(t)并不是脈沖函數(shù)，所以將(8)式寫成離散矩陣形式

(9)

(9)式可采用最小二乘方法進(jìn)行求解.

圖9展示了利用m序列的相關(guān)特性壓制噪聲的效果.圖9a為電場(chǎng)響應(yīng)ex，其中黑色實(shí)線是理論電場(chǎng)響應(yīng)，灰色實(shí)線是加入高斯噪聲后電場(chǎng)響應(yīng)，其信噪比為30 dB.圖9b為電場(chǎng)ex響應(yīng)與發(fā)射波形(圖7)互相關(guān)后的結(jié)果，其中黑色實(shí)線和灰色實(shí)線分別為圖9a中理論和加入噪聲后響應(yīng)與發(fā)射波形的互相關(guān)結(jié)果，兩組結(jié)果幾乎完全重合，而均方相對(duì)誤差減少1/4多.可以看出利用m序列的相關(guān)特性可以有效地壓制噪聲干擾.

本文針對(duì)圖8中的m序列響應(yīng)進(jìn)行相關(guān)辨識(shí)，獲得脈沖響應(yīng)辨識(shí)結(jié)果(圖10)，其中黑色實(shí)線為理論脈沖響應(yīng)，十字星為相關(guān)辨識(shí)結(jié)果，可以看出在早期兩組結(jié)果吻合較好；在晚期(20 ms之后)，辨識(shí)結(jié)果出現(xiàn)微弱震蕩.

5.2 積分計(jì)算階躍響應(yīng)

根據(jù)裴易峰等(2014)，采用簡(jiǎn)單的矩形數(shù)值積分算法可以實(shí)現(xiàn)由dB/dt計(jì)算B，通過(guò)對(duì)晚期道數(shù)據(jù)進(jìn)行指數(shù)擬合可以有效壓制噪聲.

(10)

B(0)是上階躍響應(yīng)的起始值，為0.針對(duì)圖10中利用相關(guān)辨識(shí)方法獲得的脈沖響應(yīng)，采用此方法計(jì)算階躍響應(yīng).

圖9 相關(guān)壓制噪聲效果Fig.9 Noise suppression via correlation

圖10 脈沖辨識(shí)結(jié)果Fig.10 Impulse response identified by deconvolution

圖11是理論階躍響應(yīng)與積分結(jié)果對(duì)比圖，其中黑色實(shí)線是理論階躍響應(yīng)，十字星是積分結(jié)果.可以發(fā)現(xiàn)兩組結(jié)果十分吻合，平均相對(duì)誤差為0.08%.積分運(yùn)算有效地壓制了辨識(shí)脈沖響應(yīng)晚期震蕩產(chǎn)生的影響.

5.3 數(shù)據(jù)反演

本文對(duì)積分獲得的階躍響應(yīng)利用Occam方法(Constable et al., 1987)進(jìn)行反演，此方法在目標(biāo)函數(shù)中加入光滑度約束，使得在滿足數(shù)據(jù)擬合的同時(shí)，獲得最光滑的反演模型.其目標(biāo)函數(shù)定義為

(11)

式中?是粗糙度矩陣，m是模型參數(shù)，μ是拉格朗日乘子，W是均方差矩陣，F(xiàn)(m)是正演響應(yīng).通過(guò)求目標(biāo)函數(shù)的極值，可以獲得新模型

(12)

圖11 理論階躍響應(yīng)與積分結(jié)果對(duì)比Fig.11 Comparison of theoretical step response with integrated result

其中J是雅克比矩陣.采用迭代技術(shù)，對(duì)反演模型進(jìn)行更新.在每次迭代過(guò)程中搜索合適的拉格朗日乘子，調(diào)整光滑約束項(xiàng)和數(shù)據(jù)擬合項(xiàng)的相對(duì)權(quán)重，使反演模型逐漸逼近真實(shí)解，最終獲得最光滑的反演模型.此方法反演穩(wěn)定，受初始模型影響小.

6 算例

為了檢驗(yàn)MTEM系統(tǒng)對(duì)高阻儲(chǔ)集層的勘探能力，本文設(shè)計(jì)了儲(chǔ)層埋深不同的兩個(gè)模型.針對(duì)淺部高阻模型，發(fā)射高階、碼元寬度較小的m序列，增加高頻發(fā)射成分，提高分辨率；對(duì)于深部高阻模型，發(fā)射低階、碼元寬度較大的m序列，增強(qiáng)低頻發(fā)射成分，加大勘探深度.

6.1 模型一

為檢驗(yàn)MTEM裝置對(duì)淺部高阻儲(chǔ)層的反映能力，本文首先設(shè)計(jì)了一個(gè)三層介質(zhì)模型(圖12)，在50 Ωm的均勻半空間中存在500 Ωm的高阻儲(chǔ)集層，埋深為300 m.發(fā)射源長(zhǎng)度為50 m，發(fā)射波形為8階m序列偽隨機(jī)碼，碼元寬度為1 ms(圖13).根據(jù)Ziolkowski等(2007)，當(dāng)收發(fā)距大于2倍儲(chǔ)集層埋深時(shí)，脈沖響應(yīng)的峰值強(qiáng)度和出現(xiàn)時(shí)刻會(huì)出現(xiàn)明顯的異常.基于這一結(jié)論，本文選擇了900，1000 m，1100 m，1200 m和1300 m五組收發(fā)距進(jìn)行觀測(cè).

圖14是收發(fā)距為1000 m的m序列電場(chǎng)響應(yīng).其中黑色實(shí)線為不含有高阻層的均勻半空間響應(yīng)，灰色實(shí)線為含有高阻儲(chǔ)層的電磁響應(yīng).由圖可以看

圖12 淺部高阻儲(chǔ)集層模型Fig.12 Model for a shallow hydrocarbon reservoir

圖13 8階m序列發(fā)射波形Fig.13 8-order m-sequence transmitting waveform

出高阻儲(chǔ)層對(duì)電場(chǎng)的影響非常大.在on-time和off-time與均勻半空間模型均有較大差別，在off-time出現(xiàn)一個(gè)較高的負(fù)向峰值.在正演數(shù)據(jù)中加入高斯噪聲，信噪比為65 dB，然后利用相關(guān)辨識(shí)技術(shù)獲得脈沖響應(yīng)，并積分計(jì)算階躍響應(yīng).圖15為淺層高阻儲(chǔ)層存在與否對(duì)脈沖和階躍響應(yīng)影響對(duì)比圖.圖15a是利用相關(guān)辨識(shí)技術(shù)獲得的脈沖響應(yīng)，圖15b是對(duì)圖15a中脈沖響應(yīng)進(jìn)行積分計(jì)算的階躍響應(yīng).其中黑色實(shí)線為不含高阻儲(chǔ)層的均勻半空間響應(yīng)，而灰色實(shí)線為含有高阻儲(chǔ)層的響應(yīng)曲線.由圖可以看出，由于噪聲的存在使得脈沖響應(yīng)出現(xiàn)了震蕩，而積分方法有效地壓制了噪聲干擾，獲得光滑的階躍響應(yīng)曲線.與均勻半空間情況相比，高阻儲(chǔ)層的出現(xiàn)使得脈沖響應(yīng)出現(xiàn)較高的峰值，而且峰值時(shí)間向早期道偏移.除此之外，階躍響應(yīng)的振幅也明顯增強(qiáng).

通過(guò)移動(dòng)發(fā)射源位置，獲得多組共中心點(diǎn)道集(CMP)數(shù)據(jù)，利用Occam方法進(jìn)行反演.圖16為不同收發(fā)距的反演結(jié)果，其中黑色粗實(shí)線、黑色細(xì)實(shí)線、灰色粗實(shí)線、黑色虛線和灰色細(xì)實(shí)線分別是收發(fā)距為900、1000、1100、1200 m和1300 m的反演結(jié)果，灰色虛線是對(duì)以上5個(gè)收發(fā)距數(shù)據(jù)同時(shí)反演結(jié)果.可以看出單個(gè)收發(fā)距數(shù)據(jù)的反演結(jié)果均可以反映出地下電性分布的總體趨勢(shì)，然而對(duì)高阻儲(chǔ)層位置和電阻率的反映不是很準(zhǔn)確.與其相比，利用5個(gè)收發(fā)距數(shù)據(jù)同時(shí)反演獲得的結(jié)果更加接近真實(shí)模型，界面位置和高阻儲(chǔ)層電阻率都比較準(zhǔn)確.因此，利用MTEM裝置同時(shí)測(cè)量多道數(shù)據(jù)探測(cè)淺部高阻儲(chǔ)層十分有效.

6.2 模型二

為檢驗(yàn)MTEM裝置對(duì)深部高阻儲(chǔ)層的反映能力，本文再次設(shè)計(jì)了一個(gè)三層介質(zhì)模型(圖17)，在50 Ωm的均勻半空間中存在500 Ωm的高阻儲(chǔ)集層，埋深為2000 m.為增加信號(hào)強(qiáng)度發(fā)射源長(zhǎng)度設(shè)為100 m，發(fā)射波形為4階m序列偽隨機(jī)碼，碼元寬度為30 ms(圖18).我們選擇收發(fā)距分別為5000、5200、5500、5700 m和6000 m.

圖15 淺部高阻儲(chǔ)層對(duì)脈沖及階躍響應(yīng)影響 (a) 相關(guān)辨識(shí)脈沖響應(yīng)dex/dt；(b)積分階躍響應(yīng)ex.Fig.15 Effects of shallow hydrocarbon reservoir on impulse and step response (a) Impulse response by correlation identification; (b) Step response obtained by integration.

圖16 不同收發(fā)距反演結(jié)果Fig.16 Inversion results for different offsets

圖17 深部高阻儲(chǔ)集層模型Fig.17 Model for a deep hydrocarbon reservoir

圖18 4階m序列發(fā)射波形Fig.18 4-order m-sequence transmitting waveform

圖19是收發(fā)距為5000 m的m序列電場(chǎng)ex響應(yīng).其中黑色實(shí)線為不含有高阻層的均勻半空間響應(yīng)，灰色實(shí)線為含有高阻儲(chǔ)層的電磁響應(yīng).與圖14相比，由于高阻層埋深加大，收發(fā)距增加，高阻儲(chǔ)層的出現(xiàn)對(duì)電場(chǎng)ex的影響減弱，但在on-time和off-time與均勻半空間模型均有一定差異.在正演數(shù)據(jù)中加入高斯噪聲，信噪比為65 dB，然后利用與模型一相同的數(shù)據(jù)處理方法進(jìn)行脈沖響應(yīng)辨識(shí)，并積分計(jì)算階躍響應(yīng).

圖20為深部高阻儲(chǔ)層存在與否對(duì)脈沖和階躍響應(yīng)影響對(duì)比圖.圖20a和圖20b分別是利用相關(guān)辨識(shí)技術(shù)獲得的脈沖響應(yīng)以及積分計(jì)算的階躍響應(yīng).其中黑色和灰色實(shí)線分別為不含高阻儲(chǔ)層的均勻半空間響應(yīng)以及含有高阻儲(chǔ)層的響應(yīng)曲線.與淺部高阻層的影響(圖15)相比，深部高阻層的出現(xiàn)同樣使得脈沖響應(yīng)產(chǎn)生較高的峰值，峰值出現(xiàn)時(shí)刻發(fā)生偏移，但幅度明顯減弱，且寬度變大.對(duì)階躍響應(yīng)影響的幅度也明顯減弱.

利用Occam方法對(duì)不同收發(fā)距的共中心點(diǎn)道集(CMP)數(shù)據(jù)進(jìn)行反演，可以獲得多組反演結(jié)果(圖21).圖中黑色細(xì)實(shí)線、灰色細(xì)實(shí)線、黑色粗實(shí)線、黑色虛線和灰色粗實(shí)線分別是收發(fā)距為5000、5200、5500、5700 m和6000 m的反演結(jié)果，灰色虛線是對(duì)以上5個(gè)收發(fā)距數(shù)據(jù)同時(shí)反演結(jié)果.由圖可以看出，高阻層上下層的電阻率得到很好的反映，高阻層的深度和厚度也得到很好的反映.與淺部高阻層的反演結(jié)果相似(參見(jiàn)圖16)，利用5個(gè)收發(fā)距數(shù)據(jù)同時(shí)處理可以獲得更加理想的反演結(jié)果.因此，利用MTEM裝置同時(shí)測(cè)量多道數(shù)據(jù)探測(cè)深部高阻儲(chǔ)層同樣十分有效.

圖19 m序列電磁響應(yīng)Fig.19 EM responses for m-sequence

圖20 深部高阻儲(chǔ)層對(duì)脈沖及階躍響應(yīng)影響 (a)相關(guān)辨識(shí)脈沖響應(yīng)dex/dt；(b)積分階躍響應(yīng)ex.Fig.20 Effects of deep hydrocarbon reservoir on impulse and step responses (a) Impulse response by correlation identification; (b) Step response obtained by integration.

圖21 不同收發(fā)距反演結(jié)果Fig.21 Inversion results for different offsets

7 結(jié)論

多通道瞬變電磁方法采用m序列發(fā)射波形，同時(shí)觀測(cè)同線電場(chǎng)分量多道數(shù)據(jù)，在油氣儲(chǔ)層勘探中可取得較好的應(yīng)用效果.頻譜分析結(jié)果顯示m序列具有較寬的頻帶，譜線之間等間距，而且包含直流成分.通過(guò)改變碼元寬度和階數(shù)可以調(diào)整頻帶寬度以及能量分布，從而達(dá)到提高信噪比和分辨率的目的.根據(jù)實(shí)際勘探目標(biāo)合理調(diào)整m序列的相關(guān)參數(shù)，可以獲得理想的勘探效果.本文利用方波響應(yīng)移位疊加以及電流對(duì)時(shí)間導(dǎo)數(shù)與階躍響應(yīng)褶積的方法實(shí)現(xiàn)了對(duì)理論m序列和實(shí)際發(fā)射波形電磁響應(yīng)的全時(shí)模擬.采用相關(guān)辨識(shí)技術(shù)去除波形，計(jì)算脈沖響應(yīng).對(duì)其進(jìn)行積分可以獲得階躍響應(yīng).m序列尖銳的相關(guān)特性，可以有效壓制噪聲干擾.直接對(duì)處理得到的階躍響應(yīng)數(shù)據(jù)進(jìn)行反演，大大提高了工作效率.本文給出的詳細(xì)m序列發(fā)射波形正演模擬以及后續(xù)數(shù)據(jù)處理流程同樣適用于2D和3D模型，為MTEM電磁勘查系統(tǒng)研發(fā)提供理論保障，必將推動(dòng)該項(xiàng)技術(shù)的發(fā)展和應(yīng)用.

致謝 作者向中國(guó)科學(xué)院電子研究所武欣及本電磁研究團(tuán)隊(duì)的其他相關(guān)人員在文章準(zhǔn)備過(guò)程中提供的幫助致謝.特別向?qū)徃迦撕途庉媽?duì)本文提出的建設(shè)性意見(jiàn)表示感謝.

Chen S D, Lin J, Zhang S. 2012. Effect of transmitter current waveform on TEM response.ChineseJ.Geophys. (in Chinese), 55(2): 709-716, doi: 10.6038/j.issn.0001-5733.2012.02.035.Chen X D, Zhang Y, Zhang J, et al. 2012. The applications of HTc SQUID magnetometer to TEM.ChineseJ.Geophys. (in Chinese), 55(2): 702-708, doi: 10.6038/j. issn.0001-5733.2012.02.034.Chun S H, Chen R J, Geng M H. 2014. Review of the pseudo-random m sequence and its application in electrical prospecting of exploration geophysics.ProgressinGeophysics(in Chinese), 29(1): 439-446, doi: 10.6038/pg20140164.Constable S C, Parker R L, Constable C G. 1987. Occam′s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data.Geophysics, 52(3): 289-300.

Guan S S, Lin J, Ji Y J, et al. 2012. Analysis of excitation signal airborne transient effect on grounded electrical-source electromagnetic response.ChineseJournalofRadioScience(in Chinese), 27(4): 766-772.He J S. 1998. 2nPseudo-random signal and its application. ∥ Annual of the Chinese Geophysical Society (in Chinese), Xi′an: Xi′an Map Press.

He J S, Tong T G, Liu J X. 2009. Mathematical analysis and realization ofansequence pseudo-random multi-frequencies signal.JournalofCentralSouthUniversity(ScienceandTechnology) (in Chinese), 40(6): 1666-1671.

Li M, Wei W B, Luo W B, et al. 2012. Time-domain spectral induced polarization based on pseudo-random sequence.PureandAppliedGeophysics, 170(12): 2257-2262, doi: 10.1007/ s00024-012-0624-z.

Li X, Xue G Q, Liu Y A, et al. 2012. A research on TEM imaging method based on synthetic-aperture technology.ChineseJ.Geophys. (in Chinese), 55(1): 333-340, doi: 10.6038/j.issn.0001-5733.2012.01.034.

Lin Z H, Chen S Y, Wang Y Y. 2009. m sequence and its application in communication.ModernElectronicsTechnique(in Chinese), 32(9): 49-51, doi: 10.3969/j.issn.1004-373X.2009.09.014.

Liu G M. 1998. Effect of transmitter current waveform on airborne TEM response.ExplorationGeophysics, 29(2): 35-41.

Pei Y F, Yin C C, Liu Y H, et al. 2014. Calculation and application of B-field for time-domain airborne EM.ProgressinGeophysics(in Chinese), 29(5): 2191-2196, doi: 10.6038/pg20140530.Smith R S, Lee T J. 2004. Asymptotic expansions for the calculation of the transient electromagnetic fields induced by a vertical magnetic dipole source above a conductive halfspace.PureandAppliedGeophysics, 161(2): 385-397.

Tang J T, Li F, Luo W B. 2007. Electrical fine-exploration transmitter design based on invert-repeated m-sequence.ProgressinGeophysics(in Chinese), 22(3): 994-1000.Tang J T, Luo W B. 2008. Pseudo-random electromagnetic exploration based on Invert-Repeated m-Sequence correlation identification.ChineseJ.Geophys. (in Chinese), 51(4): 1226-1233, doi: 10.3321/j.issn:0001-5733.2008.04.034.Wright D, Ziolkowski A, Hobbs B. 2002. Hydrocarbon detection and monitoring with a multicomponent transient electromagnetic (MTEM) survey.TheLeadingEdge, 21(9): 852-864, doi: 10.1190/1.2144389.Wright D, Ziolkowski A. 2007. Suppression of noise in MTEM data. ∥ SEG Technical Program Expanded Abstracts, 549-553.

Wiener N. 1949. Extrapolation, Interpolation, and Smoothing of Stationary Time Series. New York: Wiley.

Wu M J, Du T C. 2002. The pseudonoise sequence code and the birth for computer.JournalofLiaoningTechnicalUniversity(NaturalScience) (in Chinese), 21(2): 203-206.

Xie X K, Yang X R. 2010. Analysis of anti-jamming to Pseudo-Random code fuse signal.ElectronicInformationWarfareTechnology(in Chinese), 25(3): 48-53.

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.

Yin C C, Huang W, Ben F. 2013. The full-time electromagnetic modeling for time-domain airborne electromagnetic system.ChineseJ.Geophys. (in Chinese), 56(9): 3153-3162, doi: 10.6038/cig20130928.

Ziolkowski A, Hobbs B A, Wright D. 2007. Multitransient electromagnetic demonstration survey in France.Geophysics, 72(4): F197-F209, doi: 10.1190/1.2735802.

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, doi: 10.1111/j.1365-2478.2010.00882.x.

Ziolkowski A, Wright D, Mattsson J. 2011. Comparison of PRBS and square-wave transient CSEM data over Peon gas discovery, Norway. ∥ SEG Technical Program Expanded Abstracts, 583-588, doi: 10.1190/1.3628149.

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

陳曙東, 林君, 張爽. 2012. 發(fā)射電流波形對(duì)瞬變電磁響應(yīng)的影響. 地球物理學(xué)報(bào), 55(2): 709-716, doi: 10.6038/j.issn.0001-5733.2012 02.035.

陳曉東, 趙毅, 張杰等. 2012. 高溫超導(dǎo)磁強(qiáng)計(jì)在瞬變電磁法中的應(yīng)用. 地球物理學(xué)報(bào), 55(2): 702-708, doi: 10.6038/j.issn.0001-5733.2012.02.034.

淳少恒, 陳儒軍, 耿明會(huì). 2014. 偽隨機(jī)m序列及其在電法勘探中的應(yīng)用進(jìn)展. 地球物理學(xué)進(jìn)展, 29(1): 439-446, doi: 10.6038/pg20140164.

關(guān)珊珊, 林君, 稽艷鞠等. 2012. 激勵(lì)信號(hào)對(duì)地-空瞬變電磁響應(yīng)的影響分析. 電波科學(xué)學(xué)報(bào), 27(4): 766-772.

何繼善. 1998. 2n系列偽隨機(jī)信號(hào)及應(yīng)用. ∥ 中國(guó)地球物理學(xué)會(huì)年刊. 西安: 西安地圖出版社.

何繼善, 佟鐵鋼, 柳建新. 2009.an序列偽隨機(jī)多頻信號(hào)數(shù)學(xué)分析及實(shí)現(xiàn). 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 40(6): 1666-1671.

李貅, 薛國(guó)強(qiáng), 劉銀愛(ài)等. 2012. 瞬變電磁合成孔徑成像方法研究. 地球物理學(xué)報(bào), 55(1): 333-340, doi: 10.6038/j.issn.0001-5733.2012.01.034.

林智慧, 陳綏陽(yáng), 王元一. 2009. m序列及其在通信中的應(yīng)用. 現(xiàn)代電子技術(shù), 32(9): 49-51, doi: 10.3969/j.issn.1004-373X.2009.09.014.

裴易峰, 殷長(zhǎng)春, 劉云鶴等. 2014. 時(shí)間域航空電磁磁場(chǎng)計(jì)算與應(yīng)用. 地球物理學(xué)進(jìn)展, 29(5): 2191-2196, doi: 10.6038/pg20140530.

湯井田, 李飛, 羅維斌. 2007. 基于逆重復(fù)m序列的精細(xì)探測(cè)電法發(fā)送機(jī)設(shè)計(jì). 地球物理學(xué)進(jìn)展, 22(3): 994-1000.

湯井田, 羅維斌. 2008. 基于相關(guān)辨識(shí)的逆重復(fù)m序列偽隨機(jī)電磁法. 地球物理學(xué)報(bào), 51(4): 1226-1233, doi: 10.3321/j.issn:0001-5733.2008.04.034.

吳明捷, 杜天蒼. 2002. 偽隨機(jī)碼及計(jì)算機(jī)的產(chǎn)生. 遼寧工程技術(shù)大學(xué)學(xué)報(bào)(自然科學(xué)版), 21(2): 203-206.

謝雪康, 楊曉蓉. 2010. 復(fù)雜電磁環(huán)境下偽隨機(jī)信碼引信信號(hào)抗干擾性分析. 電子信息對(duì)抗技術(shù), 25(3): 48-53.

薛國(guó)強(qiáng), 陳衛(wèi)營(yíng), 周楠楠等. 2013. 接地源瞬變電磁短偏移深部探測(cè)技術(shù). 地球物理學(xué)報(bào), 56(1): 255-261, doi: 10.6038/cjg20130126.

殷長(zhǎng)春, 黃威, 賁放. 2013. 時(shí)間域航空電磁系統(tǒng)瞬變?nèi)珪r(shí)響應(yīng)正演模擬. 地球物理學(xué)報(bào), 56(9): 3153-3162, doi: 10.6038/cig20130928.

(本文編輯 何燕)

Multi-transient EM full-time forward modeling and inversion of m-sequences

QI Yan-Fu1, YIN Chang-Chun1*, WANG Ruo2, CAI Jing1

1CollegeofGeo-explorationSciencesandTechnology,JilinUniversity,Changchun130021,China2InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China

The traditional transient electromagnetic method is mainly used in mineral exploration. It is difficult to meet the demand for oil and gas exploration. The multi-transient electromagnetic (MTEM) system can solve this problem due to its grounded wire as transmitter and in-line electric field component, resulting in the underground resistive targets can be more easily detected. However, the domestic research for this method is still at the stage of theoretical derivation, lacking forward modeling and data processing. We work on full-time forward modelling and inversion for m-sequence MTEM method to provide theoretical basis for the on-going development of MTEM system as well as its data processing and interpretation.

10.6038/cjg20150731.

國(guó)家自然科學(xué)基金項(xiàng)目(41274121，41174111)、國(guó)家重大科研裝備研制項(xiàng)目(ZDYZ2012-1-03，ZDYZ2012-1-05-04)聯(lián)合資助.

齊彥福，男，1989年生，博士，主要從事地球物理電磁正反演理論和方法技術(shù)研究. E-mail：jdqiyanfu@126.com

*通訊作者 殷長(zhǎng)春，男，1965年生，教授，國(guó)家“千人計(jì)劃”特聘專家，主要從事電磁勘探理論，特別是航空和海洋電磁方面的研究. E-mail: yinchangchun@jlu.edu.cn

10.6038/cjg20150731

P631

2014-11-24，2015-03-19收修定稿

齊彥福, 殷長(zhǎng)春，王若等. 2015. 多通道瞬變電磁m序列全時(shí)正演模擬與反演.地球物理學(xué)報(bào)，58(7):2566-2577,

Qi Y F, Yin C C, Wang R, et al. 2015. Multi-transient EM full-time forward modeling and inversion of m-sequences.ChineseJ.Geophys. (in Chinese),58(7):2566-2577,doi:10.6038/cjg20150731.