武欣， 薛國強(qiáng)， 底青云， 張一鳴， 方廣有

1 中國科學(xué)院電子學(xué)研究所， 北京 100190 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 北京工業(yè)大學(xué)， 北京 100124

?

偽隨機(jī)編碼源電磁響應(yīng)的精細(xì)辨識

武欣1， 薛國強(qiáng)2， 底青云2， 張一鳴3， 方廣有1

1 中國科學(xué)院電子學(xué)研究所， 北京 100190 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 北京工業(yè)大學(xué)， 北京 100124

與傳統(tǒng)階躍源激勵方式相比，采用m序列偽隨機(jī)編碼對發(fā)射源波形進(jìn)行編碼，提高了電磁探測的深度和分辨能力.然而受這種編碼源發(fā)射波形自相關(guān)旁瓣效應(yīng)的影響，使得對大地沖激響應(yīng)的精細(xì)辨識效果受到一定限制.為了解決這一問題，在以往相關(guān)辨識方法研究的基礎(chǔ)上，進(jìn)一步考慮發(fā)射自相關(guān)旁瓣的影響，首先提出一種由收發(fā)互相關(guān)中高精度提取大地沖激響應(yīng)的數(shù)學(xué)方法；然后通過數(shù)值模擬給出了m序列編碼源大地沖激響應(yīng)的精細(xì)辨識結(jié)果；同時對以m序列為發(fā)射波形的勘探系統(tǒng)相關(guān)參數(shù)選擇進(jìn)行了分析；最后利用本文提出的方法對野外實測數(shù)據(jù)進(jìn)行辨識處理，通過與其他EM方法結(jié)果進(jìn)行對比，證明了本文提出方法的可靠性.

m序列； 偽隨機(jī)二進(jìn)制序列； 系統(tǒng)辨識； 大地沖激響應(yīng)； 自相關(guān)旁瓣效應(yīng)

The fundamental of system identification using the m sequence is the Wiener-Hopf function, which reveals the relation between auto-correlation of input signal AR, identified system g and cross-correlation between input and output signal CR that is the convolution of AR andg. Because of the complicated effect of the sidelobes of AR, even if the main lobe of AR could be regarded as a delta function, g cannot yet be read out directly from CR with high precision. From the Wiener-Hopf function, it can be recognized that the effect of the sidelobes comes into the CR through the convolution process. This means that the AR could be regarded as a filter, and therefore, a deconvolution could be performed to extract g from CR by inverting the effect of the convolution with AR. For this purpose, CR, AR and g are primarily discretized and reformed individually into vectors and matrixes, and the convolution process is consequently reformed into vector=matrix×vector, where the left vector is CR, matrix is AR and right vector is g. The AR matrix could be regarded as a linear operator. Therefore the least square procedure could be introduced for the g vector estimation. Through several iterations and smoothing, g could be extracted with high precision.

To check the effectiveness of this method, a numerical simulation is performed. A grounded line source on a homogeneous half-space is computed as the model, and the transmitting waveforms are coded by the m sequence. The result shows that the relative identification error is smaller than 2% in the time interval from 0 to 10tpeak, wheretpeakis the arrival time of the peak of the earth impulse response. On contrary the relative identification error with the traditional method in the same interval is much worse, and the biggest relative error is more than 538%. In another numerical simulation, it is proved that this method can provide identification result with high quality even if the choice of m sequence parameters is considered not fully enough from the perspective of some existing theories. Furthermore, we organized a field test in 2014 and compared the identified earth impulse response with that of other EM methods. The comparison shows that the identified result is matched well with other EM methods.

The sidelobes can be regarded as a disturbing factor in identification, and its effect comes into cross-correlation through the convolution process. The new algorithm proposed in this paper realized the high precision identification of earth impulse response by deconvolution. Numerical simulation shows that the relative identification error using this method is smaller than 2% in the main section of the earth impulse response. This method is also applied to field data processing already, and the identified earth impulse response is in good agreement with that of other EM methods. All these tests prove that this method is an applicable method in the future EM surveys.Keywords m sequence；PRBS；System Identification；Electromagnetic impulse response of the earth；Sidelobe effect of auto-correlation

1 引言

勘探地球物理電磁方法往往通過對大地電磁響應(yīng)的觀測、數(shù)據(jù)處理以及反演等流程，建立地下電性結(jié)構(gòu)模型，以實現(xiàn)對地下目標(biāo)體的識別.對于所研究對象，其電磁感應(yīng)過程從整體上可以看作為一個線性、時不變過程，將系統(tǒng)辨識的相關(guān)理論與方法引入到電磁法精細(xì)勘探中，成為目前研究的熱點.

盡管國內(nèi)外對基于m序列發(fā)射波形的EM方法研究方興未艾，然而，如何能夠在時域高精度地辨識出大地系統(tǒng)的電磁沖激響應(yīng)依然是一個研究難點.本文在基本相關(guān)辨識方法的基礎(chǔ)上，進(jìn)一步考慮發(fā)射自相關(guān)旁瓣的影響，提出一種數(shù)學(xué)方法實現(xiàn)大地沖激響應(yīng)的高精度辨識.在此基礎(chǔ)上，本文對實際觀測中的m序列碼型參數(shù)選擇展開研究，提出選擇m序列主要參數(shù)的方法.最后，利用本文提出的方法對野外實測數(shù)據(jù)進(jìn)行辨識處理，通過與其他方法結(jié)果進(jìn)行對比，驗證本文提出方法的可靠性.

2 基本相關(guān)辨識方法與旁瓣效應(yīng)問題

利用m序列進(jìn)行系統(tǒng)辨識的基本原理基于Wiener-Hopf方程，對于如圖1所示的一個線性時不變系統(tǒng)，表達(dá)式可為

y(t)=z(t)+v(t)，

(1)

z(t)=g(t)*u(t)，

(2)

其中u(t)為辨識輸入信號，g(t)為系統(tǒng)的沖激響應(yīng)，z(t)為u(t)經(jīng)過系統(tǒng)的輸出信號，v(t)為噪聲，y(t)是包含噪聲的輸出信號.基于圖1的系統(tǒng)描述，首先僅考慮u(t)與z(t)，Wiener-Hopf方程可寫為

圖1 典型線性系統(tǒng)數(shù)據(jù)流程圖Fig.1 Data flow chart of a typical linear system

CR(z,u)=g*AR(u)，

(3)

其中，CR(l,m)表示信號l與m的互相關(guān)，AR(l)表示信號l的自相關(guān).公式(3)表示輸入信號的自相關(guān)與系統(tǒng)沖激響應(yīng)的卷積為輸入信號與輸出信號的互相關(guān).考慮噪聲n，公式可為

CR(y,u) =CR(z,u)+CR(n,u)

=g*AR(u)+CR(n,u),

(4)

公式(4)即為基本相關(guān)辨識方法的數(shù)學(xué)描述，以公式(4)為基礎(chǔ)選擇輸入信號，其要求包括：辨識輸入信號應(yīng)具有隨機(jī)性，使式中的CR(n,u)可忽略；辨識輸入信號自相關(guān)近似δ(t)函數(shù)，則使g*AR(u)近似等于g(t)成立.

m序列能夠較好地滿足上述要求，因此成為一種常見的系統(tǒng)辨識輸入信號.利用m序列對電法發(fā)射電流波形進(jìn)行編碼，并按照系統(tǒng)辨識方法提取大地沖激響應(yīng)，其信號流程如圖2.假設(shè)觀測到的實際發(fā)射源電流波形為Tx(t)(后稱“實際發(fā)射電流觀測波形”)，大地沖激響應(yīng)為g(t)，在觀測點觀測到的大地響應(yīng)為Rx(t).Tx(t)和Rx(t)可被表述為

Rx(t)=g(t)*Iw(t)*hr(t)+v，

(5)

Tx(t)=Iw(t)*htr(t)，

(6)

其中Iw(t)為發(fā)射機(jī)實際輸出的發(fā)射電流波形(即實際對大地產(chǎn)生激勵的波形)，htr(t)為用于記錄Iw(t)的接收機(jī)的系統(tǒng)響應(yīng)，hr(t)為用于觀測大地響應(yīng)的接收機(jī)的系統(tǒng)響應(yīng)，v為噪聲.

在實際工作中，一次獨立的發(fā)射過程主要包含三個參數(shù)：階數(shù)n、碼元寬度ts以及循環(huán)發(fā)射次數(shù)Ncyc，即循環(huán)發(fā)射Ncyc次參數(shù)為n和ts的m序列.在發(fā)射過程中，一方面由于發(fā)射系統(tǒng)的外部阻抗環(huán)境(尤其對于電性源)對于發(fā)射波形中不同頻率成分的阻抗特性不一致，導(dǎo)致不同頻率成分與地下的耦合特性不同；另一方面，不同的發(fā)射時間和地點，系統(tǒng)外部阻抗環(huán)境也不同.因此，對每一次獨立發(fā)射過程，均需要對實際發(fā)射電流波形Iw(t)與大地響應(yīng)Iw(t)*g(t)進(jìn)行完整觀測，并將兩者觀測數(shù)據(jù)Tx(t)與Rx(t)整體應(yīng)用于數(shù)據(jù)處理中(Wright et al., 2002; Ziolkowski et al., 2006).

圖2 辨識系統(tǒng)信號流程圖Fig.2 Flow chart of signal identification system

然而，與常見m序列自相關(guān)函數(shù)(隱含周期性條件)不同，有限次循環(huán)m序列的整體自相關(guān)為非周期函數(shù)(如圖3)，其包含若干幅度不同的類δ(t)函數(shù)尖峰，同時在各尖峰兩側(cè)存在具有震蕩特性的旁瓣.在實際工作中，由于大地對發(fā)射電流的畸變作用以及受接收機(jī)htr(t)特性的影響，發(fā)射自相關(guān)旁瓣的震蕩將更為復(fù)雜.由公式(4)可見，收發(fā)互相關(guān)是大地沖激響應(yīng)g(t)與發(fā)射自相關(guān)的卷積，故自相關(guān)旁瓣的影響也會隨著卷積過程進(jìn)入收發(fā)互相關(guān)，從而導(dǎo)致無法按照基本相關(guān)辨識方法直接從收發(fā)互相關(guān)中高精度地提取出大地沖激響應(yīng).

圖3 連續(xù)循環(huán)5個周期m序列的整體自相關(guān)Fig.3 Auto-correlation of an m sequence including 5 continuous cycles

3 大地沖激響應(yīng)的精確辨識方法

3.1 精確辨識算法

作為一種基于系統(tǒng)辨識的方法，本文并未使用術(shù)語“系統(tǒng)響應(yīng)”描述使用接收機(jī)觀測實際發(fā)射電流波形的過程.在早期MTEM研究(項目THERMIE project OG/0305/92/NL-UK)中使用術(shù)語“系統(tǒng)響應(yīng)”源自Strack的著作(Strack, 1992)，認(rèn)為大地沖激響應(yīng)是輸入信號，而由實際發(fā)射波形與數(shù)據(jù)采集設(shè)備整體引起的對大地沖激響應(yīng)觀測的失真總效果是系統(tǒng).此系統(tǒng)可被看作一個濾波器，其響應(yīng)即“系統(tǒng)響應(yīng)”.由于要求上述濾波器中僅包含收發(fā)系統(tǒng)的特性，在早期Duncan的研究中將之稱作“零大地異常濾波器”.可見，使用術(shù)語“系統(tǒng)響應(yīng)”，即認(rèn)為對大地沖激響應(yīng)的提取過程是對輸入信號的重構(gòu).而本文所討論的方法認(rèn)為大地為待辨識系統(tǒng)，實際發(fā)射電流觀測波形Tx(t)是對辨識輸入信號Iw(t)的觀測.因此，盡管實現(xiàn)“系統(tǒng)響應(yīng)”觀測和“實際發(fā)射電流波形”觀測的方法類似，但其體現(xiàn)的大地沖激響應(yīng)提取邏輯并不相同.

對于公式(5)，基于Wiener-Hopf方程計算Rx(t)與Tx(t)的互相關(guān)，得到：CR(Rx(t),Tx(t))=CR(g(t)*Iw(t)*hr(t),Tx(t)) +CR(v,Tx(t)).

(7)

根據(jù)互相關(guān)的性質(zhì)，因為Tx(t)為實信號，有：CR(Rx(t),Tx(t))=Rx(t)*Tx(-t)=g(t)*Iw(t)

*hr(t)*Tx(-t)+v*Tx(-t)，

(8)其中，Tx(-t)為Tx(t)的反折.利用公式(8)，將公式(7)中不含噪聲的部分改寫為

g(t)*Iw(t)*hr(t)*Iw(-t)*htr(-t)

=g(t)*AR(Iw(t))*CR(htr(t),hr(t)).

(9)

由公式(9)可見，收發(fā)互相關(guān)中除了大地沖激響應(yīng)外，實際還包含了：(1)實際發(fā)射波形Iw(t)的自相關(guān)；(2)用于記錄發(fā)射波形的接收機(jī)與用于觀測大地響應(yīng)的接收機(jī)系統(tǒng)響應(yīng)之間的互相關(guān).其中，htr(t)和hr(t)可以在實驗室中觀測得到，但I(xiàn)w(t)無法獨立獲得.在實際操作中，一種簡便的方法是使用實際觀測到的發(fā)射波形自相關(guān)AR(Tx(t))來替代AR(Iw(t))*CR(htr(t),hr(t))，即：

AR(Tx(t)) =Iw(t)*Iw(-t)*htr(t)*htr(-t)

=Iw(t)*Iw(-t)*hr(t)*htr(-t).

(10) 將公式(10)帶入公式(9)，則有：

g(t)*Iw(t)*hr(t)*Iw(-t)*htr(-t)

=g(t)*AR(Tx(t))，

(11)

公式(10)能夠成立的條件是htr(t)=hr(t)——要求用于發(fā)射波形觀測的系統(tǒng)與用于大地響應(yīng)觀測系統(tǒng)具有相同的系統(tǒng)響應(yīng).如果在實際中不能滿足htr(t)=hr(t)假設(shè)，則需要提前對兩者的系統(tǒng)響應(yīng)進(jìn)行觀測，獲得兩者的關(guān)系函數(shù)f(t)，從而構(gòu)成公式(12)為

htr(t)=f(t)*hr(t).

(12)

根據(jù)公式(12)，可在相關(guān)計算前使用關(guān)系函數(shù)f(t)對響應(yīng)Rx(t)進(jìn)行處理，得到公式為Rxf(t)=Rx(t)*f(t)=g(t)*Iw(t)*htr(t)

+v*f(t)=g(t)*Iw(t)*htr(t)+v，

(13)其中噪聲v與f(t)的卷積依然是噪聲，故繼續(xù)使用v來表示.計算Rxf(t)與Tx(t)的互相關(guān)，其中不含噪聲的部分為

g(t)*Iw(t)*htr(t)*Iw(-t)*htr(-t)

=g(t)*AR(Tx(t)).

(14)

基于公式(12)得到公式(14)的形式與公式(11)相同，以此為基礎(chǔ)，公式(8)可改寫為

CR(Rx(t),Tx(t))=g(t)*AR(Tx(t))

+CR(v,Tx(t))，

(15)

簡化代換可得：

(16)

則公式(15)可被改寫為

a(t)=g(t)*b(t)+cv，

(17)

從本質(zhì)上講cv依然是噪聲，但其體現(xiàn)了使用m序列的噪聲抑制能力.從辨識的角度，也可以認(rèn)為cv體現(xiàn)了由于噪聲存在而產(chǎn)生的辨識偏差.a(t)、g(t)及b(t)的關(guān)系在離散系統(tǒng)可寫為

(18)

其中，Ng為辨識大地沖激響應(yīng)的采樣點數(shù).

將公式(18)矩陣化，可得：

A=BG+cv，

(19)

其中：

式中n1和n2分別為a(n)和b(n)序列中最大值采樣點的序列號.矩陣B主對角線上的元素為自相關(guān)最大值，其他元素為旁瓣值.公式(19)表明矩陣B描述了大地沖激響應(yīng)與收發(fā)互相關(guān)之間的映射關(guān)系，其實質(zhì)是自相關(guān)主瓣與旁瓣對大地沖激響應(yīng)的加權(quán)和作用.矩陣B滿足線性算子條件，以公式(19)為模型，建立上述映射的逆過程，利用最小二乘辨識方法并經(jīng)過若干次迭代，即可實現(xiàn)對大地沖激響應(yīng)的高精度辨識，處理流程如圖4.

3.2 數(shù)值模擬與辯識效果

假設(shè)使用m序列對電性源EM方法發(fā)射波形進(jìn)行編碼，發(fā)射電極間距1000 m，發(fā)射電流10 A.在發(fā)射電極軸向延長線上布設(shè)觀測點(X軸向)，觀測點距收發(fā)電極中心1000 m.觀測電場X分量，觀測電極間距為單位長度，其他計算參數(shù)如表1.

表1 數(shù)值模擬計算參數(shù)Table 1 Parameters of simulation

為了更好地模擬真實觀測，使用理想大地階躍響應(yīng)ex(t)[公式(20)]與實際發(fā)射波形時間導(dǎo)數(shù)?Iw(t)/?t卷積生成觀測數(shù)據(jù)，得到的觀測數(shù)據(jù)時間序列(經(jīng)電流和發(fā)射電極間距歸一化)如圖5a.發(fā)射電流波形Iw(t)如圖5b，Iw(t)自相關(guān)(經(jīng)電流歸一化)如圖5c，而電流波形與觀測波形互相關(guān)(經(jīng)電流和發(fā)射電極間距歸一化)如圖5d.利用公式(19)描述的過程對大地沖激響應(yīng)進(jìn)行辨識，辨識過程使用最小二乘算法，并經(jīng)過了4次迭代.辨識結(jié)果與由公式(21)計算的理論沖激響應(yīng)比較如圖5e.為了體現(xiàn)辨識算法高精度的辨識能力，圖5e中的辨識曲線并沒有經(jīng)過平滑處理.圖5f給出利用公式(22)計算的辨識誤差，其中g(shù)r是辨識大地沖激響應(yīng).公式(20)、(21)和(22)分別為

(20)

(21)

(22)

由上述數(shù)值模擬計算可以看出，辨識算法自辨識起點至約10倍tpeak(8.38 ms)時段內(nèi)的辨識誤差小于2%，這對于多數(shù)情況是可以接受的.如果進(jìn)一步對辨識算法進(jìn)行優(yōu)化，并對辨識后數(shù)據(jù)進(jìn)行適當(dāng)?shù)钠交?，可得到精度更高的辨識結(jié)果.

圖4 基于自相關(guān)旁瓣去除的大地沖激響應(yīng)辨識方法流程Fig.4 Block diagram of the identification method for earth impulse response based on cancellation of auto-correlation sidelobes

圖5d表示收發(fā)互相關(guān)，在此算例中，也表示理論大地沖激響應(yīng)與電流波形自相關(guān)的卷積結(jié)果.由圖5c我們可以看到在電流波形自相關(guān)各極大值外側(cè)均存在劇烈震蕩的旁瓣，這種現(xiàn)象也會隨卷積過程對收發(fā)互相關(guān)產(chǎn)生影響.如果不對這種影響進(jìn)行消除，而直接將數(shù)據(jù)帶入公式(4)中進(jìn)行辨識，則這種影響將最終作用于辨識結(jié)果中：圖6a為利用圖5d的數(shù)據(jù)按照公式(4)進(jìn)行辨識得到的結(jié)果，圖6b為圖6a的辨識誤差.由圖6b可見，圖6a的辨識誤差在約10倍tpeak(8.38 ms)時已達(dá)到538.6%.可見，在辨識過程中需要考慮自相關(guān)的旁瓣效應(yīng)，否則受其影響的辨識結(jié)果將給后續(xù)處理造成很大的困難.3.3 以m序列為發(fā)射波形的勘探系統(tǒng)相關(guān)參數(shù)選擇

在實際工作中，選擇合適的m序列參數(shù)對于高效完成對大地沖激響應(yīng)辨識具有重要意義.為了能夠有效地辨識大地沖激響應(yīng)，Ziolkowski提出一種簡單確定m序列碼元寬度ts的方法(Ziolkowski, 2007).假設(shè)沖激響應(yīng)的峰值出現(xiàn)時間為tpeak，為了保證m序列能夠?qū)peak與之前的空氣波脈沖清晰地分辨開，Ziolkowski提出碼元寬度ts應(yīng)至少等于tpeak/10.但如此估算出的ts可能導(dǎo)致實際中的發(fā)射機(jī)無法輸出設(shè)計的波形.以前述的算例為例，如圖4e所示，tpeak出現(xiàn)在約0.838 ms，則按要求ts應(yīng)為83.8 μs.目前市場上常見的電性源發(fā)射機(jī)，最大升壓水平約為1000 V，對于一般接地條件，當(dāng)ts小于100 μs時則很難保證發(fā)射波形的完整.然而，使用本文所介紹的辨識方法，利用實際發(fā)射波形的自相關(guān)對旁瓣效應(yīng)進(jìn)行去除，則ts的選擇將可突破tpeak/10的限制.

為此，可以在前述模擬計算的基礎(chǔ)上添加一個空氣波，其脈沖寬度為3個采樣點(基于24000 Hz采樣率)，脈沖幅度為沖激響應(yīng)峰值的15倍.我們依然使用碼元寬度ts為1/6000 s≈167 μs，采用同樣的辨識過程，辨識結(jié)果如圖7a.由圖7a可見，盡管此時的tpeak與ts并不滿足10倍關(guān)系，但依然實現(xiàn)了對大地沖激響應(yīng)主體部分的辨識.當(dāng)然，選擇的ts越短，能夠辨識出的大地沖激響應(yīng)的范圍也越大.使用相同的模型，改變采樣率為48000 Hz，設(shè)計ts=1/12000 s≈83 μs≈tpeak/10，保持空氣波脈沖時長不變，則辨識出的大地沖激響應(yīng)如圖7b.對比圖7a和7b可見，采用更短的ts可以辨識出大地沖激響應(yīng)曲線上更早的部分.事實上，如有必要，取更窄的ts甚至可以在一定程度上實現(xiàn)對空氣波的辨識.另一方面，ts的取值也不能過寬.使用相同的模型，采樣率為24000 Hz，設(shè)計ts=1/1200 s≈833 μs ≈tpeak，此時的辨識結(jié)果如圖7c.顯然，取此碼元寬度，系統(tǒng)已喪失對大地沖激響應(yīng)的辨識能力.

通過上述討論可知，ts的選取需要綜合考慮多種因素.在能夠滿足發(fā)射機(jī)發(fā)射能力的條件下，同時考慮觀測系統(tǒng)的存儲能力，可以盡可能選擇短的ts，以獲得更完整的大地沖激響應(yīng)辨識結(jié)果.

圖6 簡單引用公式(4)的辨識結(jié)果(a) 辨識結(jié)果；(b)辨識誤差.Fig.6 Identification of impulse response of the earth with Eq(4)(a) Identification of impulse response of the earth; (b) Error of Identification.

對編碼長度N(或者編碼階數(shù)n)的選擇也同樣需要考慮多種因素.僅從理論上講，N越長越好，因為利用m序列對系統(tǒng)信噪比的提升大約在到N的范圍.然而N越大，觀測時間越長，觀測過程并不經(jīng)濟(jì)；另外，當(dāng)N過大時，m序列對信噪比的提升能力會下降.對于N的選擇主要考慮噪聲水平ηr，其不僅包括硬件系統(tǒng)的噪聲，同時也考慮環(huán)境噪聲的影響.假設(shè)在發(fā)射機(jī)停止發(fā)射后tg時段內(nèi)，觀測信號幅度衰減至ηr以下(也可以將m序列對信噪比的大致提升效果考慮在內(nèi))，則發(fā)射結(jié)束后采樣數(shù)Ng=tg×Fsr，其中Fsr為觀測系統(tǒng)的采樣率.可以根據(jù)公式(23)來估計編碼長度，公式(23)為

(23)

上述參數(shù)確定過程中，不能僅依賴?yán)碚撚嬎?，還需要結(jié)合測試數(shù)據(jù)以實現(xiàn)最優(yōu)參數(shù)選擇.

4 野外實測數(shù)據(jù)處理

為驗證本文提出的大地沖激響應(yīng)精確辨識方法的可靠性，我們于2014年初在河北省固安縣進(jìn)行了方法驗證試驗.試驗主要參考MTEM系統(tǒng)，也依據(jù)我們已有的一些研究成果(Wu et al., 2012; 張建國等，2014)完成發(fā)射波形相關(guān)參數(shù)確定、觀測裝置選擇等.

收發(fā)裝置布設(shè)如圖8，其中AB點為發(fā)射電極埋設(shè)位置，P點是測點位置.發(fā)射電極布設(shè)沿用國家公益性行業(yè)專項“深部探測技術(shù)與實驗研究”課題SinoProbe-09-02野外試驗所采用的發(fā)射電極布設(shè)位置，東西向，極間距1 km.試驗過程中發(fā)射電極間接地電阻約為40 Ω.發(fā)射機(jī)使用Zonge GGT-30，使用改進(jìn)型的XMT-32驅(qū)動GGT-30實現(xiàn)m序列編碼的輸出.在試驗中測試了多種m序列參數(shù)組合，其碼元寬度ts包括1/2048 s、1/4096 s、1/5765 s以及1/8192 s；階數(shù)n包括12、13、15以及18.此外，還使用中國科學(xué)院電子所自主研發(fā)的碼型驅(qū)動器實現(xiàn)對GGT-30的驅(qū)動輸出，除了測試與上述改進(jìn)型XMT-32相同的m序列參數(shù)，還測試了其他碼元寬度.對實際輸出電流波形的觀測采用寬帶寬、大動態(tài)范圍的霍爾傳感器，并使用與用于大地響應(yīng)觀測相同的接收機(jī)記錄電流波形.

觀測點P位于發(fā)射電極東側(cè)的軸向延長線上，距發(fā)射電極中心的距離為1500 m，進(jìn)行電場觀測，接收電極與發(fā)射電極線方向平行，極距160 m，此收發(fā)方式與MTEM方法相同.根據(jù)Wright的研究(Wright, 2003)，除上述與發(fā)射電極共線的電場觀測，其他場量的觀測并不能提供更多的信息，因此MTEM在后期的野外工作中也僅觀測與發(fā)射共線的電場分量(Ziolkowski et al., 2007).觀測用的接收機(jī)為SinoProbe-09-02項目研發(fā)的SEP接收機(jī)，采樣率24000 Hz.

以一組觀測結(jié)果為例說明利用本文方法對大地沖激響應(yīng)的辨識過程，其使用的m序列參數(shù)為：階數(shù)n=18，碼元寬度ts=1/8192 s，發(fā)射電流幅度16.7A.根據(jù)先驗知識，碼元寬度以小于1/2000 s為宜.如前所述，試驗中測試了多組碼元寬度，包括1/2048 s、1/4096 s、1/5765 s以及1/8192 s.不同的碼元寬度意味著m序列帶寬不同，為能夠更完整地實現(xiàn)大地沖激響應(yīng)辨識，在條件允許的情況下應(yīng)盡量選擇更窄的碼元寬度.經(jīng)測試，在當(dāng)?shù)亟拥貤l件下發(fā)射機(jī)對碼元寬度為1/8192 s的編碼波形輸出質(zhì)量滿足處理要求，且實測數(shù)據(jù)未出現(xiàn)飽和現(xiàn)象，因此可以用于大地沖激響應(yīng)辨識.

圖7 碼元寬度選擇對辨識結(jié)果的影響(a) ts=1/6000 s時的辨識結(jié)果；(b) ts=1/12000 s時的辨識結(jié)果；(c) ts=1/1200 s時的辨識結(jié)果.Fig.7 The Identification results with different ts(a) Identification result with ts=1/6000 s；(b) Identification result with ts=1/12000 s；(c) Identification result with ts=1/1200 s.

圖8 野外試驗收發(fā)裝置布設(shè)Fig.8 Survey layout in the field experiment

圖9 野外測試數(shù)據(jù)與辨識結(jié)果(a) 野外觀測數(shù)據(jù)實測波形(片段)；(b) 野外實際發(fā)射波形(片段)；(c) 野外數(shù)據(jù)辨識結(jié)果(含空氣波)；(d) 野外數(shù)據(jù)辨識結(jié)果(沖激響應(yīng)主體部分).Fig.9 Field test data and Identification result(a) Real time series of field recorded signal (segment)；(b) Real time series of transmitting current (segment)；(c) Identification result (airwave included)；(d) Identification result (main part of earth impulse response).

圖10 相同測區(qū)MT視電阻率圖Fig.10 Apparent resistivity (MT) of the same survey area

觀測數(shù)據(jù)飽和是應(yīng)該盡量避免的.由于高頻分量隨偏移距增大而衰減較塊，因此當(dāng)偏移距較小時容易首先在低頻分量上出現(xiàn)飽和現(xiàn)象.對m序列與基頻等于1/(Nts)的方波進(jìn)行比較，在發(fā)射能量相同條件下，由于m序列將發(fā)射能量均勻地分布于其帶寬內(nèi)，而方波將發(fā)射能量集中于基頻與諧波頻點上，故m序列低頻頻點所分得的能量一般要比方波各頻點(尤其是基頻)所分得的能量少很多.因此，使用m序列出現(xiàn)低頻飽和現(xiàn)象的最小收發(fā)偏移距也將比使用方波出現(xiàn)相同現(xiàn)象的最小偏移距小很多.本次試驗中，采用階數(shù)參數(shù)n=18來配合偏移距、碼元寬度及循環(huán)次數(shù)等參數(shù)，使整個m序列帶寬內(nèi)能量分布適中，既保證高頻分量能夠被有效觀測到，也確保低頻分量不會因為偏移距較小而飽和.可見，在野外工作中，m序列的參數(shù)需要綜合考慮.若需要進(jìn)一步減小或增大偏移距，為保證實現(xiàn)有效探測，則m序列及發(fā)射各參數(shù)均需進(jìn)行適當(dāng)調(diào)整.

實測觀測數(shù)據(jù)(片段，經(jīng)發(fā)射電流與發(fā)射電極距歸一化)如圖9a，實測發(fā)射波形(片段)如圖9b，辨識出的大地沖激響應(yīng)(包含空氣波與大地沖激響應(yīng)主體部分)如圖9c，而辨識出的大地沖激響應(yīng)主體部分的細(xì)節(jié)如圖9d. SinoProbe-09-02項目的研發(fā)人員在相同區(qū)域進(jìn)行過多次MT、CSAMT試驗，獲得該區(qū)域的視電阻率圖(MT方法)如圖10.參考圖9d，根據(jù)Ziolkowski提出的快速估計均勻半空間大地沖激響應(yīng)峰值時間的公式(24)反推大地電阻率(Ziolkowski, 2007)，可快速判斷地下中淺部視電阻率約為25 Ωm左右，這與圖10的結(jié)果基本吻合.公式(24)為

(24)

5 結(jié)語

傳統(tǒng)EM由于受激勵源帶寬的限制無法對大地系統(tǒng)的電磁感應(yīng)過程進(jìn)行持續(xù)激勵，這導(dǎo)致了傳統(tǒng)EM方法在建立大地電性結(jié)構(gòu)分布模型過程中存在局限性.按照系統(tǒng)辨識的思想，使用m序列對EM方法的激勵波形進(jìn)行編碼，可以在一定程度上打破上述局限性，提高EM勘探的深度與分辨能力.本文基于基本相關(guān)辨識提出的新方法，充分考慮了發(fā)射波形自相關(guān)，尤其是其旁瓣在收發(fā)互相關(guān)中產(chǎn)生的復(fù)雜作用，并使用反卷積的方法對其影響予以消除，從而有效提高辨識精度.通過對野外實測數(shù)據(jù)的處理對本方法做出驗證，其辨識結(jié)果與其他地球物理EM方法得到的結(jié)果相吻合，證明了本文提出方法的可靠性.在野外工作開始前，需要結(jié)合工作區(qū)域的實際情況對m序列的相關(guān)參數(shù)進(jìn)行優(yōu)化選擇.

(1)對于碼元寬度ts的選擇主要依據(jù)對大地沖激響應(yīng)峰值時刻tpeak的估計以及發(fā)射機(jī)的實際發(fā)射能力.對tpeak的估計需要以收發(fā)偏移距以及大地淺層電阻率做參數(shù)，因此對tpeak的估計在一定程度上依賴先驗知識.

(2)在ts選定后可對接收機(jī)的采樣率Fsr進(jìn)行選擇，保證接收機(jī)的采樣率約為4到5倍的1/ts即可.這也即意味著，當(dāng)收發(fā)偏移距增大導(dǎo)致ts延長時，可選擇較小的采樣率.

(3)對編碼長度N或階數(shù)n的選擇主要依賴噪聲水平，其包括觀測系統(tǒng)的硬件噪聲水平也包括環(huán)境噪聲水平.根據(jù)噪聲水平大體確定發(fā)射結(jié)束后采樣數(shù)Ng，再通過Ng、Fsr以及ts來估算m序列的最小編碼長度.

致謝 感謝“千人計劃”吉林大學(xué)殷長春對本文撰寫所提供的幫助，感謝吳凱、黃江杰、張樂、張建國參與了方法試驗的野外工作，感謝黃玲博士提供了MT對比數(shù)據(jù)，感謝中國科學(xué)院地質(zhì)與地球物理研究所的安志國博士對野外試驗進(jìn)行了實地指導(dǎo).

Duncan P M. 1978. Electromagnetic deep crustal sounding with a controlled Pseudo-noise source [Ph. D. thesis]. Canada: University of Toronto.

Duncan P M, Hwang A, Edwards R N, et al. 1980. The development and application of a wide-band electromagnetic sounding system using a pseudo noise source.Geophysics, 45(8): 1276-1296.

Edwards R N, Holladay J S. 1991. Pseudo-Random Binary Sequence (PRBS) Techniques. ∥ Nabighian M N, et al. Electromagnetic Methods in Applied Geophysics—Applications Part A and Part B. USA: Society of Exploration Geophysicists. Luo W B, Li Q C, Tang J T. 2012. Coded source electromagnetic sounding method.ChineseJ.Geophys. (in Chinese), 55(1): 341-349, doi: 10.6038/j.issn.0001-5733.2012.01.035.

He J S. 2010. Wide field electromagnetic methods and Pseudo-random signal electro-prospecting. Beijing: higher education press.

Luo Y Z. 2015. The second scheme algorithm of Instruments for Pseudo-random signal electro-prospecting. ∥ 2015 annual symposium of academic conference of geophysical technology.Strack K M. 1992. Exploration with Deep Transient Electromagnetics. Amsterdam: Elsevier.Tang J T, Lou W B. 2008. Pseudo-random electromagnetic exploration based on Invert-Repeated m Sequence correlation identification.ChineseJ.Geophys. (in Chinese), 51(4): 1226-1233.Wright D, Ziolkowski A, Hall G. 2006. Improving Signal-to-noise ratio using pseudo random binary sequences in Multi Transient Electromagnetic (MTEM) Data. ∥ 68th EAGE Conference and Exhibition, Vienna 2006. Expanded Abstracts, P065.

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. 2003. Detection of hydrocarbons and their movement in a reservoir using time-lapse multi-transient electromagnetic (MTEM) data [Ph. D. thesis]. UK: University of Edinburgh. Wu X, Zhu W H, Fang G Y, et al. 2012. Optimization of PRBS Coded Transmitting current Parameters for Electric Source and Array receiving TEM. ∥ 18th European Meeting of Environmental and Engineering Geophysics Paris, France.Zepernick H J, Filger A. 2007. Pseudo Random Signal Processing Theory and Application (in Chinese). Gan L C, et al. Trans. Beijing: Publishing House of Electronics Industry.

Zhang J G, Wu X, Qi Y Z, et al. 2014. Investigation on time domain coded electromagnetic exploration method.JournalofRadars, 3(2): 158-165.

Zhao B R, Zhao J, Zhang H K, et al. 2006. The PS100 high precision earth-electricity instrument system (IP to IP) with controllable source-application of CDMA technology to the measurement of earth-resistivity for the first time.ProgressinGeophysics(in Chinese), 21(2): 675-682, doi: 10.3969/j.issn.1004-2903.2006.02.053.

Ziolkowski A. 2007. Development in the transient electromagnetic method.FirstBreak, 25: 99-106.

Ziolkowski A, Hobbs B, Dawes G, et al. 2006. True Amplitude Transient electromagnetic system response Measurement. International Patent: WO 2006/114561 A2.Ziolkowski A, Hobbs B, Wright D. 2007. Multitransient electromagnetic demonstration survey in France.Geophysics, 72(4): F197-F209.

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

何繼善. 2010. 廣域電磁法和偽隨機(jī)信號電法. 北京: 高等教育出版社.

羅維斌, 李慶春, 湯井田. 2012. 編碼電磁測深. 地球物理學(xué)報, 55(1): 341-349, doi: 10.6038/j.issn.0001-5733.2012.01.035.

羅延鐘. 2015. 偽隨機(jī)信號電法儀器的“第二方案”算法. 地球物理技術(shù)委員會2015年學(xué)術(shù)年會.

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

澤珀尼克 H J, 費爾閣 A. 2007. 偽隨機(jī)信號處理理論與應(yīng)用. 甘良才等譯. 北京: 電子工業(yè)出版社.

張建國, 武欣, 齊有政等. 2014. 時間域編碼電磁勘探方法研究. 雷達(dá)學(xué)報, 3(2): 158-165.

趙碧如, 趙健, 張洪魁等. 2006. PS100型IP到端可控源高精度大地電測儀系統(tǒng)—CDMA技術(shù)首次在地電阻測量中的應(yīng)用. 地球物理學(xué)進(jìn)展, 21(2): 675-682, doi: 10.3969/j.issn.1004-2903.2006.02.053.

(本文編輯 張正峰)

Accurate identification for the electromagnetic impulse response of the earth with pseudo random coded waveforms

WU Xin1, XUE Guo-Qiang2, DI Qing-Yun2, ZHANG Yi-Ming3, FANG Guang-You1

1InstituteofElectronics,ChineseAcademyofSciences,Beijing100190,China2InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China3BeijingUniversityofTechnology,Beijing100124,China

The electromagnetic method with m sequence, one kind of pseudo random binary sequences, that codes transmitting waveforms has several advantages over the traditional electromagnetic method with square-form transmitting waveforms. However, because of the existing sidelobes of auto-correlation of transmitting waveforms and its influence in a complicated way, the precision of identification would be limited if the identification process is only based on a basic cross-correlation algorithm. A new algorithm is proposed to overcome this difficulty.

國家重大科研裝備研制項目(ZDYZ2012-1-05-01)資助.

武欣，男，1982年生，助理研究員，主要從事時間域電磁法理論與系統(tǒng)研發(fā).E-mail: wu_xin18@mail.ie.ac.cn

10.6038/cjg20150814.Wu X, Xue G Q, Di Q Y, et al. 2015. Accurate identification for the electromagnetic impulse response of the earth with pseudo random coded waveforms.ChineseJ.Geophys. (in Chinese),58(8):2792-2802,doi:10.6038/cjg20150814.

10.6038/cjg20150814

P631

2014-12-02，2015-06-24收修定稿

武欣，薛國強(qiáng)，底青云等. 2015. 偽隨機(jī)編碼源電磁響應(yīng)的精細(xì)辨識.地球物理學(xué)報，58(8):2792-2802,