伪随机编码源电磁响应的精细辨识

武欣， 薛国强， 底青云， 张一鸣， 方广有

1 中国科学院电子学研究所， 北京 100190 2 中国科学院地质与地球物理研究所， 北京 100029 3 北京工业大学， 北京 100124



伪随机编码源电磁响应的精细辨识

武欣1， 薛国强2， 底青云2， 张一鸣3， 方广有1

1 中国科学院电子学研究所， 北京 100190 2 中国科学院地质与地球物理研究所， 北京 100029 3 北京工业大学， 北京 100124

与传统阶跃源激励方式相比，采用m序列伪随机编码对发射源波形进行编码，提高了电磁探测的深度和分辨能力.然而受这种编码源发射波形自相关旁瓣效应的影响，使得对大地冲激响应的精细辨识效果受到一定限制.为了解决这一问题，在以往相关辨识方法研究的基础上，进一步考虑发射自相关旁瓣的影响，首先提出一种由收发互相关中高精度提取大地冲激响应的数学方法；然后通过数值模拟给出了m序列编码源大地冲激响应的精细辨识结果；同时对以m序列为发射波形的勘探系统相关参数选择进行了分析；最后利用本文提出的方法对野外实测数据进行辨识处理，通过与其他EM方法结果进行对比，证明了本文提出方法的可靠性.

m序列； 伪随机二进制序列； 系统辨识； 大地冲激响应； 自相关旁瓣效应

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 引言

勘探地球物理电磁方法往往通过对大地电磁响应的观测、数据处理以及反演等流程，建立地下电性结构模型，以实现对地下目标体的识别.对于所研究对象，其电磁感应过程从整体上可以看作为一个线性、时不变过程，将系统辨识的相关理论与方法引入到电磁法精细勘探中，成为目前研究的热点.

尽管国内外对基于m序列发射波形的EM方法研究方兴未艾，然而，如何能够在时域高精度地辨识出大地系统的电磁冲激响应依然是一个研究难点.本文在基本相关辨识方法的基础上，进一步考虑发射自相关旁瓣的影响，提出一种数学方法实现大地冲激响应的高精度辨识.在此基础上，本文对实际观测中的m序列码型参数选择展开研究，提出选择m序列主要参数的方法.最后，利用本文提出的方法对野外实测数据进行辨识处理，通过与其他方法结果进行对比，验证本文提出方法的可靠性.

2 基本相关辨识方法与旁瓣效应问题

利用m序列进行系统辨识的基本原理基于Wiener-Hopf方程，对于如图1所示的一个线性时不变系统，表达式可为

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

(1)

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

(2)

其中u(t)为辨识输入信号，g(t)为系统的冲激响应，z(t)为u(t)经过系统的输出信号，v(t)为噪声，y(t)是包含噪声的输出信号.基于图1的系统描述，首先仅考虑u(t)与z(t)，Wiener-Hopf方程可写为

图1 典型线性系统数据流程图Fig.1 Data flow chart of a typical linear system

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

(3)

其中，CR(l,m)表示信号l与m的互相关，AR(l)表示信号l的自相关.公式(3)表示输入信号的自相关与系统冲激响应的卷积为输入信号与输出信号的互相关.考虑噪声n，公式可为

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

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

(4)

公式(4)即为基本相关辨识方法的数学描述，以公式(4)为基础选择输入信号，其要求包括：辨识输入信号应具有随机性，使式中的CR(n,u)可忽略；辨识输入信号自相关近似δ(t)函数，则使g*AR(u)近似等于g(t)成立.

m序列能够较好地满足上述要求，因此成为一种常见的系统辨识输入信号.利用m序列对电法发射电流波形进行编码，并按照系统辨识方法提取大地冲激响应，其信号流程如图2.假设观测到的实际发射源电流波形为Tx(t)(后称“实际发射电流观测波形”)，大地冲激响应为g(t)，在观测点观测到的大地响应为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)为发射机实际输出的发射电流波形(即实际对大地产生激励的波形)，htr(t)为用于记录Iw(t)的接收机的系统响应，hr(t)为用于观测大地响应的接收机的系统响应，v为噪声.

在实际工作中，一次独立的发射过程主要包含三个参数：阶数n、码元宽度ts以及循环发射次数Ncyc，即循环发射Ncyc次参数为n和ts的m序列.在发射过程中，一方面由于发射系统的外部阻抗环境(尤其对于电性源)对于发射波形中不同频率成分的阻抗特性不一致，导致不同频率成分与地下的耦合特性不同；另一方面，不同的发射时间和地点，系统外部阻抗环境也不同.因此，对每一次独立发射过程，均需要对实际发射电流波形Iw(t)与大地响应Iw(t)*g(t)进行完整观测，并将两者观测数据Tx(t)与Rx(t)整体应用于数据处理中(Wright et al., 2002; Ziolkowski et al., 2006).

图2 辨识系统信号流程图Fig.2 Flow chart of signal identification system

然而，与常见m序列自相关函数(隐含周期性条件)不同，有限次循环m序列的整体自相关为非周期函数(如图3)，其包含若干幅度不同的类δ(t)函数尖峰，同时在各尖峰两侧存在具有震荡特性的旁瓣.在实际工作中，由于大地对发射电流的畸变作用以及受接收机htr(t)特性的影响，发射自相关旁瓣的震荡将更为复杂.由公式(4)可见，收发互相关是大地冲激响应g(t)与发射自相关的卷积，故自相关旁瓣的影响也会随着卷积过程进入收发互相关，从而导致无法按照基本相关辨识方法直接从收发互相关中高精度地提取出大地冲激响应.

图3 连续循环5个周期m序列的整体自相关Fig.3 Auto-correlation of an m sequence including 5 continuous cycles

3 大地冲激响应的精确辨识方法

3.1 精确辨识算法

作为一种基于系统辨识的方法，本文并未使用术语“系统响应”描述使用接收机观测实际发射电流波形的过程.在早期MTEM研究(项目THERMIE project OG/0305/92/NL-UK)中使用术语“系统响应”源自Strack的著作(Strack, 1992)，认为大地冲激响应是输入信号，而由实际发射波形与数据采集设备整体引起的对大地冲激响应观测的失真总效果是系统.此系统可被看作一个滤波器，其响应即“系统响应”.由于要求上述滤波器中仅包含收发系统的特性，在早期Duncan的研究中将之称作“零大地异常滤波器”.可见，使用术语“系统响应”，即认为对大地冲激响应的提取过程是对输入信号的重构.而本文所讨论的方法认为大地为待辨识系统，实际发射电流观测波形Tx(t)是对辨识输入信号Iw(t)的观测.因此，尽管实现“系统响应”观测和“实际发射电流波形”观测的方法类似，但其体现的大地冲激响应提取逻辑并不相同.

对于公式(5)，基于Wiener-Hopf方程计算Rx(t)与Tx(t)的互相关，得到：CR(Rx(t),Tx(t))=CR(g(t)*Iw(t)*hr(t),Tx(t)) +CR(v,Tx(t)).

(7)

根据互相关的性质，因为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)可见，收发互相关中除了大地冲激响应外，实际还包含了：(1)实际发射波形Iw(t)的自相关；(2)用于记录发射波形的接收机与用于观测大地响应的接收机系统响应之间的互相关.其中，htr(t)和hr(t)可以在实验室中观测得到，但Iw(t)无法独立获得.在实际操作中，一种简便的方法是使用实际观测到的发射波形自相关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)——要求用于发射波形观测的系统与用于大地响应观测系统具有相同的系统响应.如果在实际中不能满足htr(t)=hr(t)假设，则需要提前对两者的系统响应进行观测，获得两者的关系函数f(t)，从而构成公式(12)为

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

(12)

根据公式(12)，可在相关计算前使用关系函数f(t)对响应Rx(t)进行处理，得到公式为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)的卷积依然是噪声，故继续使用v来表示.计算Rxf(t)与Tx(t)的互相关，其中不含噪声的部分为

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

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

(14)

基于公式(12)得到公式(14)的形式与公式(11)相同，以此为基础，公式(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)

从本质上讲cv依然是噪声，但其体现了使用m序列的噪声抑制能力.从辨识的角度，也可以认为cv体现了由于噪声存在而产生的辨识偏差.a(t)、g(t)及b(t)的关系在离散系统可写为

(18)

其中，Ng为辨识大地冲激响应的采样点数.

将公式(18)矩阵化，可得：

A=BG+cv，

(19)

其中：

式中n1和n2分别为a(n)和b(n)序列中最大值采样点的序列号.矩阵B主对角线上的元素为自相关最大值，其他元素为旁瓣值.公式(19)表明矩阵B描述了大地冲激响应与收发互相关之间的映射关系，其实质是自相关主瓣与旁瓣对大地冲激响应的加权和作用.矩阵B满足线性算子条件，以公式(19)为模型，建立上述映射的逆过程，利用最小二乘辨识方法并经过若干次迭代，即可实现对大地冲激响应的高精度辨识，处理流程如图4.

3.2 数值模拟与辩识效果

假设使用m序列对电性源EM方法发射波形进行编码，发射电极间距1000 m，发射电流10 A.在发射电极轴向延长线上布设观测点(X轴向)，观测点距收发电极中心1000 m.观测电场X分量，观测电极间距为单位长度，其他计算参数如表1.

表1 数值模拟计算参数Table 1 Parameters of simulation

为了更好地模拟真实观测，使用理想大地阶跃响应ex(t)[公式(20)]与实际发射波形时间导数∂Iw(t)/∂t卷积生成观测数据，得到的观测数据时间序列(经电流和发射电极间距归一化)如图5a.发射电流波形Iw(t)如图5b，Iw(t)自相关(经电流归一化)如图5c，而电流波形与观测波形互相关(经电流和发射电极间距归一化)如图5d.利用公式(19)描述的过程对大地冲激响应进行辨识，辨识过程使用最小二乘算法，并经过了4次迭代.辨识结果与由公式(21)计算的理论冲激响应比较如图5e.为了体现辨识算法高精度的辨识能力，图5e中的辨识曲线并没有经过平滑处理.图5f给出利用公式(22)计算的辨识误差，其中gr是辨识大地冲激响应.公式(20)、(21)和(22)分别为

(20)

(21)

(22)

由上述数值模拟计算可以看出，辨识算法自辨识起点至约10倍tpeak(8.38 ms)时段内的辨识误差小于2%，这对于多数情况是可以接受的.如果进一步对辨识算法进行优化，并对辨识后数据进行适当的平滑，可得到精度更高的辨识结果.

图4 基于自相关旁瓣去除的大地冲激响应辨识方法流程Fig.4 Block diagram of the identification method for earth impulse response based on cancellation of auto-correlation sidelobes

图5d表示收发互相关，在此算例中，也表示理论大地冲激响应与电流波形自相关的卷积结果.由图5c我们可以看到在电流波形自相关各极大值外侧均存在剧烈震荡的旁瓣，这种现象也会随卷积过程对收发互相关产生影响.如果不对这种影响进行消除，而直接将数据带入公式(4)中进行辨识，则这种影响将最终作用于辨识结果中：图6a为利用图5d的数据按照公式(4)进行辨识得到的结果，图6b为图6a的辨识误差.由图6b可见，图6a的辨识误差在约10倍tpeak(8.38 ms)时已达到538.6%.可见，在辨识过程中需要考虑自相关的旁瓣效应，否则受其影响的辨识结果将给后续处理造成很大的困难.3.3 以m序列为发射波形的勘探系统相关参数选择

在实际工作中，选择合适的m序列参数对于高效完成对大地冲激响应辨识具有重要意义.为了能够有效地辨识大地冲激响应，Ziolkowski提出一种简单确定m序列码元宽度ts的方法(Ziolkowski, 2007).假设冲激响应的峰值出现时间为tpeak，为了保证m序列能够将tpeak与之前的空气波脉冲清晰地分辨开，Ziolkowski提出码元宽度ts应至少等于tpeak/10.但如此估算出的ts可能导致实际中的发射机无法输出设计的波形.以前述的算例为例，如图4e所示，tpeak出现在约0.838 ms，则按要求ts应为83.8 μs.目前市场上常见的电性源发射机，最大升压水平约为1000 V，对于一般接地条件，当ts小于100 μs时则很难保证发射波形的完整.然而，使用本文所介绍的辨识方法，利用实际发射波形的自相关对旁瓣效应进行去除，则ts的选择将可突破tpeak/10的限制.

为此，可以在前述模拟计算的基础上添加一个空气波，其脉冲宽度为3个采样点(基于24000 Hz采样率)，脉冲幅度为冲激响应峰值的15倍.我们依然使用码元宽度ts为1/6000 s≈167 μs，采用同样的辨识过程，辨识结果如图7a.由图7a可见，尽管此时的tpeak与ts并不满足10倍关系，但依然实现了对大地冲激响应主体部分的辨识.当然，选择的ts越短，能够辨识出的大地冲激响应的范围也越大.使用相同的模型，改变采样率为48000 Hz，设计ts=1/12000 s≈83 μs≈tpeak/10，保持空气波脉冲时长不变，则辨识出的大地冲激响应如图7b.对比图7a和7b可见，采用更短的ts可以辨识出大地冲激响应曲线上更早的部分.事实上，如有必要，取更窄的ts甚至可以在一定程度上实现对空气波的辨识.另一方面，ts的取值也不能过宽.使用相同的模型，采样率为24000 Hz，设计ts=1/1200 s≈833 μs ≈tpeak，此时的辨识结果如图7c.显然，取此码元宽度，系统已丧失对大地冲激响应的辨识能力.

通过上述讨论可知，ts的选取需要综合考虑多种因素.在能够满足发射机发射能力的条件下，同时考虑观测系统的存储能力，可以尽可能选择短的ts，以获得更完整的大地冲激响应辨识结果.

图6 简单引用公式(4)的辨识结果(a) 辨识结果；(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(或者编码阶数n)的选择也同样需要考虑多种因素.仅从理论上讲，N越长越好，因为利用m序列对系统信噪比的提升大约在到N的范围.然而N越大，观测时间越长，观测过程并不经济；另外，当N过大时，m序列对信噪比的提升能力会下降.对于N的选择主要考虑噪声水平ηr，其不仅包括硬件系统的噪声，同时也考虑环境噪声的影响.假设在发射机停止发射后tg时段内，观测信号幅度衰减至ηr以下(也可以将m序列对信噪比的大致提升效果考虑在内)，则发射结束后采样数Ng=tg×Fsr，其中Fsr为观测系统的采样率.可以根据公式(23)来估计编码长度，公式(23)为

(23)

上述参数确定过程中，不能仅依赖理论计算，还需要结合测试数据以实现最优参数选择.

4 野外实测数据处理

为验证本文提出的大地冲激响应精确辨识方法的可靠性，我们于2014年初在河北省固安县进行了方法验证试验.试验主要参考MTEM系统，也依据我们已有的一些研究成果(Wu et al., 2012; 张建国等，2014)完成发射波形相关参数确定、观测装置选择等.

收发装置布设如图8，其中AB点为发射电极埋设位置，P点是测点位置.发射电极布设沿用国家公益性行业专项“深部探测技术与实验研究”课题SinoProbe-09-02野外试验所采用的发射电极布设位置，东西向，极间距1 km.试验过程中发射电极间接地电阻约为40 Ω.发射机使用Zonge GGT-30，使用改进型的XMT-32驱动GGT-30实现m序列编码的输出.在试验中测试了多种m序列参数组合，其码元宽度ts包括1/2048 s、1/4096 s、1/5765 s以及1/8192 s；阶数n包括12、13、15以及18.此外，还使用中国科学院电子所自主研发的码型驱动器实现对GGT-30的驱动输出，除了测试与上述改进型XMT-32相同的m序列参数，还测试了其他码元宽度.对实际输出电流波形的观测采用宽带宽、大动态范围的霍尔传感器，并使用与用于大地响应观测相同的接收机记录电流波形.

观测点P位于发射电极东侧的轴向延长线上，距发射电极中心的距离为1500 m，进行电场观测，接收电极与发射电极线方向平行，极距160 m，此收发方式与MTEM方法相同.根据Wright的研究(Wright, 2003)，除上述与发射电极共线的电场观测，其他场量的观测并不能提供更多的信息，因此MTEM在后期的野外工作中也仅观测与发射共线的电场分量(Ziolkowski et al., 2007).观测用的接收机为SinoProbe-09-02项目研发的SEP接收机，采样率24000 Hz.

以一组观测结果为例说明利用本文方法对大地冲激响应的辨识过程，其使用的m序列参数为：阶数n=18，码元宽度ts=1/8192 s，发射电流幅度16.7A.根据先验知识，码元宽度以小于1/2000 s为宜.如前所述，试验中测试了多组码元宽度，包括1/2048 s、1/4096 s、1/5765 s以及1/8192 s.不同的码元宽度意味着m序列带宽不同，为能够更完整地实现大地冲激响应辨识，在条件允许的情况下应尽量选择更窄的码元宽度.经测试，在当地接地条件下发射机对码元宽度为1/8192 s的编码波形输出质量满足处理要求，且实测数据未出现饱和现象，因此可以用于大地冲激响应辨识.

图7 码元宽度选择对辨识结果的影响(a) ts=1/6000 s时的辨识结果；(b) ts=1/12000 s时的辨识结果；(c) ts=1/1200 s时的辨识结果.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 野外试验收发装置布设Fig.8 Survey layout in the field experiment

图9 野外测试数据与辨识结果(a) 野外观测数据实测波形(片段)；(b) 野外实际发射波形(片段)；(c) 野外数据辨识结果(含空气波)；(d) 野外数据辨识结果(冲激响应主体部分).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 相同测区MT视电阻率图Fig.10 Apparent resistivity (MT) of the same survey area

观测数据饱和是应该尽量避免的.由于高频分量随偏移距增大而衰减较块，因此当偏移距较小时容易首先在低频分量上出现饱和现象.对m序列与基频等于1/(Nts)的方波进行比较，在发射能量相同条件下，由于m序列将发射能量均匀地分布于其带宽内，而方波将发射能量集中于基频与谐波频点上，故m序列低频频点所分得的能量一般要比方波各频点(尤其是基频)所分得的能量少很多.因此，使用m序列出现低频饱和现象的最小收发偏移距也将比使用方波出现相同现象的最小偏移距小很多.本次试验中，采用阶数参数n=18来配合偏移距、码元宽度及循环次数等参数，使整个m序列带宽内能量分布适中，既保证高频分量能够被有效观测到，也确保低频分量不会因为偏移距较小而饱和.可见，在野外工作中，m序列的参数需要综合考虑.若需要进一步减小或增大偏移距，为保证实现有效探测，则m序列及发射各参数均需进行适当调整.

实测观测数据(片段，经发射电流与发射电极距归一化)如图9a，实测发射波形(片段)如图9b，辨识出的大地冲激响应(包含空气波与大地冲激响应主体部分)如图9c，而辨识出的大地冲激响应主体部分的细节如图9d. SinoProbe-09-02项目的研发人员在相同区域进行过多次MT、CSAMT试验，获得该区域的视电阻率图(MT方法)如图10.参考图9d，根据Ziolkowski提出的快速估计均匀半空间大地冲激响应峰值时间的公式(24)反推大地电阻率(Ziolkowski, 2007)，可快速判断地下中浅部视电阻率约为25 Ωm左右，这与图10的结果基本吻合.公式(24)为

(24)

5 结语

传统EM由于受激励源带宽的限制无法对大地系统的电磁感应过程进行持续激励，这导致了传统EM方法在建立大地电性结构分布模型过程中存在局限性.按照系统辨识的思想，使用m序列对EM方法的激励波形进行编码，可以在一定程度上打破上述局限性，提高EM勘探的深度与分辨能力.本文基于基本相关辨识提出的新方法，充分考虑了发射波形自相关，尤其是其旁瓣在收发互相关中产生的复杂作用，并使用反卷积的方法对其影响予以消除，从而有效提高辨识精度.通过对野外实测数据的处理对本方法做出验证，其辨识结果与其他地球物理EM方法得到的结果相吻合，证明了本文提出方法的可靠性.在野外工作开始前，需要结合工作区域的实际情况对m序列的相关参数进行优化选择.

(1)对于码元宽度ts的选择主要依据对大地冲激响应峰值时刻tpeak的估计以及发射机的实际发射能力.对tpeak的估计需要以收发偏移距以及大地浅层电阻率做参数，因此对tpeak的估计在一定程度上依赖先验知识.

(2)在ts选定后可对接收机的采样率Fsr进行选择，保证接收机的采样率约为4到5倍的1/ts即可.这也即意味着，当收发偏移距增大导致ts延长时，可选择较小的采样率.

(3)对编码长度N或阶数n的选择主要依赖噪声水平，其包括观测系统的硬件噪声水平也包括环境噪声水平.根据噪声水平大体确定发射结束后采样数Ng，再通过Ng、Fsr以及ts来估算m序列的最小编码长度.

致谢 感谢“千人计划”吉林大学殷长春对本文撰写所提供的帮助，感谢吴凯、黄江杰、张乐、张建国参与了方法试验的野外工作，感谢黄玲博士提供了MT对比数据，感谢中国科学院地质与地球物理研究所的安志国博士对野外试验进行了实地指导.

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.

附中文参考文献

何继善. 2010. 广域电磁法和伪随机信号电法. 北京: 高等教育出版社.

罗维斌, 李庆春, 汤井田. 2012. 编码电磁测深. 地球物理学报, 55(1): 341-349, doi: 10.6038/j.issn.0001-5733.2012.01.035.

罗延钟. 2015. 伪随机信号电法仪器的“第二方案”算法. 地球物理技术委员会2015年学术年会.

汤井田, 罗维斌. 2008. 基于相关辨识的逆重复m序列伪随机电磁法. 地球物理学报, 51(4): 1226-1233.

泽珀尼克 H J, 费尔阁 A. 2007. 伪随机信号处理理论与应用. 甘良才等译. 北京: 电子工业出版社.

张建国, 武欣, 齐有政等. 2014. 时间域编码电磁勘探方法研究. 雷达学报, 3(2): 158-165.

赵碧如, 赵健, 张洪魁等. 2006. PS100型IP到端可控源高精度大地电测仪系统—CDMA技术首次在地电阻测量中的应用. 地球物理学进展, 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年生，助理研究员，主要从事时间域电磁法理论与系统研发.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收修定稿

武欣，薛国强，底青云等. 2015. 伪随机编码源电磁响应的精细辨识.地球物理学报，58(8):2792-2802,