时间域编码电磁勘探方法研究

张建国 武 欣 齐有政② 黄 玲 方广有

①(中国科学院电磁辐射与探测技术重点实验室 北京 100190)

②(中国科学院大学 北京 100049)

1 引言

随着我国经济的快速发展，能源需求量急剧加大，油气资源勘探的攻关研究已经迫在眉睫。油田钻井代价高昂且只能确定钻井附近的地质参数，不适合大规模普查油气资源[1,2]；对于地表及地质条件复杂的地区，电磁勘探方法能从不同于地震方法的角度反映油气层的物理特性，可实现对油气藏的有效勘探。然而，由于空气波等外界干扰，如何提高信噪比(Signal to Noise Ratio,SNR)已然成为电磁勘探方法的首要问题[3,4]。

长偏移距瞬变电磁测深法[5](Long Offset Transient ElectroMagnetics,LOTEM)是一种典型的油气资源电磁勘探方法，它以电偶源瞬变测深理论为基础，对地下深部油气层有较强的勘探能力。为了实现大勘探深度，LOTEM采用长达几百甚至几千米的发射线将大功率激励信号注入大地，并在远区(早期)观测含有地电信息的大地电磁响应，通过加大发射能量来提高 SNR。由于仪器设备以及地质条件限制，加大发射能量来提高SNR有较大的局限性。

文献[6]基于LOTEM方法，提出了MTEM(Multi-Transient ElectroMagnetics)方法，在实验后期为了提高SNR，将激励信号由双极性方波换成伪随机编码(Pseudo-Random Binary Sequences,PRBS)信号。由于 PRBS具有良好的相关性，能够提升SNR。MTEM方法能分辨出大深度下的油气层，实验结果达到了预期目标。然而，MTEM方法只是在实验后期对编码信号进行了一次尝试，对于编码信号的诸多参数选择和实验分析等关键性问题，都未进行理论阐述和实验验证。

为了抑制空气波干扰，本文提出时间域编码电磁勘探方法，采用伪随机编码作为激励信号，基于同步记录的编码信号和接收电压信号，通过时间域反卷积信号复原方法(Time-domain Deconvolution Signal Recovery,TDSR)，获得含有地下异常体信息的大地冲激响应，同时分析了时钟频率与码元长度的选择依据，并进行了外场实验验证。PRBS编码信号扩展了发射信号的频带宽度，提高了接收数据的SNR。TDSR方法能抑制空气波干扰，进而实现对地下异常体的有效勘探。

2 编码电磁勘探原理

2.1 层状大地模型的冲激响应

对于电磁勘探方法，如果地下异常体为导电性介质，如金属矿产，采用水平磁性源发射(水平闭合线圈作为发射源)最有利于电磁勘探；如果期望探测的地下异常体呈现高阻特性，如油气层，则最佳的激励源为水平电性源(两端接地的馈电导线)。针对地下高阻油气层勘探，本文采用水平电性源作为发射源。电磁勘探方法常用的水平层状大地模型如图1所示，地下各层介质的电阻率和厚度分别记为ρn和hn。由于在一般情况下发射源长度小于偏移距(发射机与接收机之间的距离)的1/4，因此可把发射源当作水平电偶极子处理，如图1中A,B所示。在距发射源一定距离(图1中r)处，采用水平长导线(图1中M,N)接收来自地下的水平感应电场电压。

图1 层状大地模型Fig.1 Layered earth model

把图1所示的大地模型看作一个大地系统，任意激励信号的系统响应均可分解为冲激响应之和。在忽略位移电流前提下，采用谐变电流作为激励信号，则大地冲激响应(对应于水平同线电场分量)的频率域表达式为[7]：

J0,J1分别为 0阶和 1阶 Bessel函数，P=I⋅ L为偶极矩，I=I0⋅ex p(−i ω⋅t)为谐变电流，L为发射源长度，μ为真空磁导率，λ为空间频率，r为偏移距，ω为角频率，为第n层的波数，角度Φ为发射点和接收点的连线与x轴的夹角，N为大地模型的总层数。

对式(1)进行傅里叶逆变换，得到大地冲激响应的时间域形式g(t)，其波形随时间变化曲线如图2所示。图中t0是空气波到达接收点的时刻，也是数据的时间零点。空气波是全空间条件下发射源激励的电磁场在接收机产生的响应，它不带有任何地下信息。在发射瞬间出现，并在偏移距较大时掩盖来自地下的有用电磁响应。大地冲激响应是发射源激励的电磁场[8]，经大地系统作用后产生的电磁响应带有地下异常体信息，且响应持续时间较长，其峰值时刻为。为了识别出地下异常体，必须将空气波从大地冲激响应数据中分离，即tp和t0之间的时间间隔必须足够大以满足数据解释的要求，下文将给出定量分析。

2.2 PRBS编码信号特性

图2 时间域大地冲激响应波形Fig.2 The impulse response waveform of the earth in time-domain

传统电磁勘探方法一般采用双极性方波作为激励信号，比如CSAMT (Controlled Source Audio Magneto Telluric)和CSEM (Controlled Source Electro Magnetic)，双极性方波的特点为：频谱只含有基频和奇次谐波，且N次谐波的振幅是基频振幅的 1/N。此特点表明双极性方波的功率主要分布于3次谐波内的3个频点，其信号频带较窄，因此，双极性方波不是宽频带信号。

为了提高电磁勘探的分辨率，采用 m 序列PRBS编码作为发射信号来扩展信号频带宽度，进而提高信号分辨率和数据质量[9,10]。PRBS编码信号的特点：(1)宽频带：信号能量在fsNP(fs为发射信号的时钟频率、NP为PRBS编码的码长)至频带内呈现出近似的均匀分布(如图3(a))，采用PRBS编码信号作为激励源相当于频率域电磁法采用多个频点同时发射，可将更多信号能量注入地下；(2)良好的相关性：PRBS编码信号的自相关函数(如图3(b))呈近似的冲激函数特性，这表明PRBS编码信号具有良好的自相关特性。根据离散型的维纳-霍夫(Wiener-Hopf)方程：

我们知道，如果输入信号x(n)的自相关函数(式(2)中φ)是冲激函数形式，那么输入信号x(n)和输出信号y(n)的互相关函数是系统冲激响应g∗(p)的最小均方估计，即在本文方法中，如采用PRBS编码信号作为激励源，则大地脉冲响应g(t)可以获得最佳近似值；因此，基于PRBS编码信号的良好相关性，采用PRBS编码信号作为电磁勘探的激励信号，同步记录PRBS编码信号和接收处的大地电磁响应，可获得精确的大地脉冲响应。

2.3 TDSR方法

把时变的PRBS编码电流信号注入大地，激励大地产生电磁感应，同步记录发射电流信号i(t)和一定偏移距处的观测电压信号V(t)，通过i(t)和V(t)的反卷积处理，可以获得带有地下异常体信息的大地冲激响应g(t)。图 4为大地系统的脉冲响应示意图，接收处的观测电压V(t)可表示为：

式中：i(t)为PRBS编码发射电流，g(t)为大地冲激响应，n(t)代表没有相关性的噪声，*表示卷积。

对于式(3)所示的反问题，通过反卷积方法进行求解，其主要处理步骤为：

(a) 利用傅里叶变换，将式(3)转换到频率域，略去频率域标记，得到：

其中： I ,G ,V 和n分别是i(t),g(t),V(t)和n(t)对应的频率域信号。

(b) 由于在接收大地响应时，也同步记录了发射的PRBS编码信号；因此，可对式(4)两边同时乘以PRBS编码信号的复共轭∗I：

图3 PRBS编码信号特点Fig.3 The characteristics of the PRBS signal

图4 大地系统的脉冲响应示意图Fig.4 Schematic diagram of the impulse response of the earth system

(c) 通过傅里叶变换将I∗⋅I转换到时间域，在理论上可得到“类冲激信号”δ(t)；但由于接收数据的长度有限及噪声的干扰，类冲激信号的长度不是无限窄，其长度由发射信号决定。因此，将I∗⋅I傅里叶逆变换后的表达式记为=F−1(I∗⋅I)，则通过傅里叶逆变换得到式(5)的时间域表达式为：

去除不相关噪声 n1(t)(I∗⋅n 对应的时间域噪声)，∗ g(t)≈g(t)即是含有地下异常体信息的大地冲激响应，即TDSR方法的目标。由式(6)可知，TDSR方法的核心思想是采用“类冲激信号”的近似值对大地冲激响应进行卷积；由于=F−1(I∗⋅I)，故的时间域长度由PRBS编码信号的时钟频率fs决定；fs越大，编码信号I中的码元间隔Δts越小，则δ～(t)的时间域波形越窄，TDSR方法获得的大地冲激响应越精确。

因此，对上述分析进行总结：TDSR方法相当于采用 PRBS编码的码元间隔Δts对大地冲激响应进行时间域采样，以恢复大地冲激响应的真实值；Δts越小(fs越大)，图 2中(tp−t0)时间段内的采样点越多，则TDSR方法对空气波的抑制能力越强，得到的结果越精确。

根据式(6)得到大地冲激响应的方式即为本文的时间域编码电磁勘探方法。

3 基于阵列式接收的编码电磁勘探系统

图5为基于阵列式接收的PRBS编码电磁勘探系统布置图。通过接地良好的发射电极，将PRBS编码电流信号注入大地系统，在其激发下，地下介质中激励起的感应涡流将产生随时间变化的感应电磁场(即2次场)。在平行于发射源方向的测线上，通过多通道接收机记录接收处的大地响应电压值，同时由于TDSR方法的需要，同步记录PRBS编码电流信号。根据获得的大地响应数据和PRBS编码信号，通过2.3节所示的反卷积复原方法，即可获得大地冲激响应。基于大地冲激响应，通过后续的数据解释和反演计算，可获得地下异常体的地电特性。

图5 基于阵列式接收的PRBS编码电磁勘探系统布置图Fig.5 The layout of the PRBS coded electromagnetic exploration system based on the array receivers

编码电磁勘探系统的发射部分采用大功率激电仪，最大发射电压达1000 V，最大工作电流为40 A。PRBS编码的时钟频率在0.022～8192 Hz范围内可选，编码长度可根据编码的阶数n(最大阶数为 15)进行调节，码长计算公式为2n−1，时间基准由GPS模块提供。

接收部分采用阵列式接收方式，如图5所示，接收阵列由多条测线组成(图中为 9条测线)，每条测线均平行于发射源方向，并等间隔分布；每条测线上有一定数量的测点，每个测点布置一个多通道接收站采集大地响应产生的感应电压。接收站为一个24位的多通道接收机，其动态范围为212 dB，最小接收电平为 ±0 .03 μ V，最大输入电压为 ±32 V。

控制部分主要负责发射部分的参数选择，如PRBS编码的时钟频率、编码阶数以及编码的周期重复次数。控制部分实现了发射信号和接收数据的时间同步，保证大地脉冲响应数据解释的正确性。

4 系统参数优化方法

4.1 PRBS编码时钟频率选取

如2.3节的分析，TDSR方法相当于采用PRBS编码的码元间隔Δts对大地冲激响应进行时间域采样，以恢复大地冲激响应的真实值；Δts越小(fs越大)，图 2中(tp−t0)时间段内的采样点越多，越有助于该方法去除空气波，进而恢复大地冲激响应的真实值；经过多次实验验证，只要码元间隔Δts达到大地冲激响应峰值时刻与空气波到达时刻两者差值的1/10，

则得到的大地冲激响应即可满足数据解释的需要，式(7)确定了PRBS编码的时钟频率fs的下限值。从去除空气波的角度考虑，fs越大越有利于结果的准确性；但是，根据大地滤波效应，随着fs增大，接收处的接收电压(式(1))会相应降低，为了获得接收处的大地冲激响应，接收电压必须大于多通道接收机的最小分辨电平，这成为PRBS编码时钟频率fs的上限值。fs的下限值和上限值是时钟频率选择的依据。

由式(1)所示的大地冲激响应表达式可知，大地感应的2次场信号与偏移距的5次方成反比，如果发射源只采用一个固定的时钟频率fs进行PRBS编码发射，那么各个测点的接收电压值会有较大数量级的差别。以图5中与发射线共线的测线5为例(中间位置与发射源同线的测线)，假定勘探深度为d，要分辨出深度d处目标的顶部和底部，最大偏移距rmax=4d，最小偏移距 rmin=0.1⋅rmax=0.4d；那么最大偏移距和最小偏移距处的接收电压比值Vmax/Vmin=10−5，二者相差5个数量级，不利于数据接收和实时处理。

因此，本文对偏移距进行分段，每个偏移距段选择一个时钟频率用于发射PRBS编码信号。采用相同时钟频率的接收点属于同一偏移距段，整个工区根据时钟频率划分为多个子分区，兼顾了各接收点的信号质量和分辨率。随着子分区与发射源之间距离的增大，时钟频率会相应降低，但每个子分区内部的时钟频率保持不变。

4.2 PRBS编码码长选取

由PRBS编码信号的特性可知：在时钟频率fs固定时，其码长越大，发射信号向地下传送的能量越多，反卷积数据处理后的SNR增益越高。确定了每个偏移距段的时钟频率后，选择码长较大的PRBS编码，有利于提高电磁勘探系统的数据质量。当时钟频率fs为定值时，接收电压数据和记录的发射电流信号进行反卷积处理，在理论上SNR增益G1等于 PRBS编码码长NP；但是，由于噪声等因素干扰，实际操作中反卷积处理获得的SNR增益G1：。

另外循环发射PRBS编码，再对得到的目标数据进行叠加也可以提高SNR。循环发射PRBS编码NC次，再对反卷积后的目标数据进行叠加，带来的SNR增益G2=。

通过增大编码长度和循环发射编码两个操作，可以显著提高数据的SNR，总的SNR增益为两个操作分别进行时SNR提升倍数的乘积，即G= G1⋅G2。

4.3 PRBS编码游程选取

基于PRBS编码信号的接收信号，可以看作是编码信号的各个上升沿和下降沿激励的大地阶跃响应的组合，如图6(a)所示。阶数为n的PRBS编码最大游程是n，最小游程是1,PRBS编码以不同的游程结尾，其对应的2次场数据SNR差别明显，如图6(b)所示，红色实线和蓝色虚线分别表示PRBS编码以最大和最小游程结尾。当以最大游程结尾时，发射信号关断后的 2次场数据的 SNR最优；当PRBS编码信号以最小游程结尾时，2次场数据的SNR最差。因此，在电磁勘探系统实际工作中，可视工区所需的 SNR条件，选择以合适游程结尾的PRBS编码信号。

5 实验结果与分析

为了验证本文编码电磁勘探方法的有效性，在河北省涿鹿县开展了外场实验，如图7和图8所示。发射线长度取500 m，发射电流10 A，接收线长度取200 m，接收间距为200 m。PRBS编码时钟频率根据4.1节所述的方法在32～8192 Hz之间选择，码元长度选取4095(12阶)和16383(14阶)两个较大的码长，每次编码发射均让PRBS编码信号循环发射5个周期，并同步记录编码发射电流信号和接收处电压信号，接收机采样率为24 kHz。

在不同时钟频率的对比实验中，采用偏移距1 km，码元长度16383(编码结尾游程为8)，选择1024 Hz和4096 Hz两个频点作为编码信号的时钟频率，在其它收发参数相同条件下，进行了对比实验。根据 2.1节中大地冲激响应的计算公式，可求解出偏移距1 km的大地冲激响应的峰值时刻为11 ms，则根据式(7)，编码时钟频率应该大于910 Hz，对比实验选择的1024 Hz和4096 Hz符合要求。编码信号发射信号和相应的接收波形分别如图 9和图 10所示，编码发射的数据信号质量较好，能有效地抑制大地电磁噪声等干扰，适合用于数据解释。基于图9和图10所示的数据，TDSR方法得到的大地冲激响应波形如图 11所示，在两个不同时钟频率条件下，得到的大地冲激响应的波形曲线基本一致，均有效地抑制了空气波干扰，且4096 Hz的大地冲激响应的峰值接近于1024 Hz峰值的1/4，与理论计算相符。从波形曲线中，可清晰地分辨出大地冲激响应的峰值时刻便于后续的数据解释和反演计算。

图7 野外实验的工区位置Fig.7 The surveying area location of the field experiments

图8 野外一个接收站的现场图 Fig.8 The figure of a receiving station in the field

图9 PRBS编码发射信号波形(时钟频率为1024 Hz)Fig.9 The transmitting waveform of the PRBS coded signal (clock frequency is 1024 Hz)

图10 编码发射的接收波形(时钟频率为1024 Hz)Fig.10 The receiving waveform corresponding to the PRBS coded signal (clock frequency is 1024 Hz)

在不同码元长度的对比实验中，仍采用上述参数，固定时钟频率为1024 Hz，分别采用4095和16383编码长度进行编码信号发射，同步记录编码发射电流信号和接收处电压信号，通过TDSR方法获得大地冲激响应波形，如图12所示。采用两个不同码元长度，得到的大地冲激响应的波形曲线基本一致，均有效地抑制了空气波干扰，从波形中可清晰地分辨出峰值时刻，且16383码长的大地冲激响应的峰值与4095码长的峰值相比，有接近于两倍的提升，与理论计算基本一致。

对于图11和图12中大地冲激响应曲线，如果成倍增加时钟频率，波形将会变得光滑，但是，随着时钟频率的增加，接收机采样率也要成倍增加，这会给数据存储和实时处理带来诸多不便，如时钟频率为1024 Hz、编码长度为4095，编码循环发射5个周期，则单点的接收数据量达到2.3 GByte。因此，采用目前编码时钟频率和接收机采样率，虽然恢复的大地冲激响应曲线不是太光滑，但仍满足数据解释的需要，为综合考虑各方面因素后的折衷选择。

图11 基于时钟频率的对比实验(码元长度为16383)Fig.11 The comparative experiments based on the clock frequency (the code length is 16383)

图12 基于码元长度的对比实验(时钟频率为1024 Hz)Fig.12 The comparative experiments based on the code length (the clock frequency is 1024 Hz)

6 结论

本文针对如何抑制电磁勘探中的空气波干扰，提出了一种时间域编码电磁勘探方法，并进行了外场验证性实验。其具体方法是采用伪随机编码序列作为激励信号，基于发射信号与接收信号之间的相关性，通过TDSR方法，获得含有地下异常体信息的大地冲激响应。实验结果表明本文方法可抑制空气波干扰，能有效识别出地下异常体信息。对比实验证实，降低时钟频率与增大编码长度可显著提高电磁数据的信号质量。在未来的工作中将继续研究编码游程对勘探结果的影响。

[1]何继善.海洋电磁法原理[M].北京: 高等教育出版社,2012:3-4.He J S.Marine Electromagnetic Principle[M].Beijing: Higher Education Press,2012: 3-4.

[2]Kaufman A A and Keller G V.频率域和时间域电磁测深[M].北京: 地质出版社,1987: 144-152.Kaufman A A and Keller G V.Frequency and Transient Soundings[M].Beijing: Geological Publishing House,1987:144-152.

[3]Ziolkowski A and Carson R.Maximising signal-to-noise ratio in transient EM data[C].EAGE 69th Conference &Exhibition,London,UK,2007: E048-E053.

[4]Ziolkowski A.Optimisation of MTEM Parameters[P]:International Patent WO 2007/104949 A1,2007.

[5]Strack K M.Exploration with Deep Transient Electromagnetics[M].USA: Elsevier Science Publishers B.V.,1992:21-24.

[6]Ziolkowski A,Hobbs B A,and Wright D.Multi-transientelectromagnetic demonstration survey in France[J].Geophysics,2007,72(4): F197-F209.

[7]汤井田,罗维斌,刘长生.海底油气藏地质模型的冲激响应[J].地球物理学报,2008,51(6): 1929-1935.Tang J T,Luo W B,and Liu C S.Impulse response of seafloor hydrocarbon reservoir model[J].Chinese Journal of Geophysics,2008,51(6): 1929-1935.

[8]牛之琏.时间域电磁法原理[M].长沙: 中南大学出版社,2007:1-10.Niu Z L.Theory of Transient Electromagnetic Methods[M].Changsha: Central South University,2007: 1-10.

[9]汤井田,罗维斌.基于相关辨识的逆重复m序列伪随机电磁法[J].地球物理学报,2008,51(4): 1226-1233.Tang J T and Luo W B.Pseudo-random electromagnetic exploration based on Invert-Repeated m-Sequence correlation identification[J].Chinese Journal of Geophysics,2008,51(4): 1226-1233.

[10]Streich R,Becken M,and Ritter O.Robust processing of noisy land-based controlled-source electromagnetic data[J].Geophysics,2013,78(5): E237-E247.