0 引言
随着地下水探寻、地下洞体探测、考古和未爆弹排查等应用需求的增加，浅层地质构造的高分辨探测越来越引起人们的关注。时间域航空电磁法是一种以飞机为运载工具进行地下目标体探测的方法，具有速度快、成本低、通行性好等优势[1]。目前国际上航空电磁浅层探测系统，如SkyTEM等，一般不改变发射基频与发射波形，而采用减小发射磁矩、缩短发射过冲时间、提早采样时间的方法获取浅层大地电磁响应。本文采用高频伪随机序列作为激励源，提高发射频率，通过系统辨识方法获取浅层目标体的电磁响应，实现航空电磁浅层探测。
目前，伪随机序列(pseudo-random binary sequence，PRBS)在电磁勘探中已经开始应用。Wright等[2]将PRBS应用到多通道瞬变电磁法(multi-channel transient electromagnetic method，MTEM)中进行电磁勘探。Ziolkowski等应用MTEM对法国西南部天然气储层测量数据进行了分析处理和反演[3]，并在挪威气田的勘探应用中发现PRBS较方波具有更强的分辨能力[4]。汤井田等[5]分析了PRBS的波形特征、相关函数及频谱特性。罗维斌等[6]在地电勘探系统中应用PRBS, 发现使用PRBS对大埋深薄层目标体具有精细的分辨能力。齐彦福等[7]实现了理论m序列和实际发射波形的全时正演模拟。武欣等[8]研究了伪随机编码源电磁响应的精细辨识问题，提出相关旁瓣效应会影响辨识结果。王若等[9]研究了伪随机编码源激发下的时间域电磁信号合成。王显祥等[10]研究了应用m序列的电磁法抗噪能力与m序列参数之间关系，并提出时间域辨识方法的抗噪能力要优于频率域辨识方法。欧阳涛等[11]通过采用伪随机编码发射技术的多道瞬变电磁法对内蒙古兴安盟某铅锌银矿进行了实例勘查，实现了深部低阻体信息的有效提取。陈超健等[12]研究了基于灰色判别准则和有理函数滤波的伪随机多频电磁信号中强干扰的压制方法。
在国内外学者对PRBS应用研究的基础上，本文针对时间域航空电磁法中应用PRBS高精度辨识大地脉冲响应的方法进行研究，并分析这种方法对噪声的压制能力。
1 伪随机源航空电磁正演在航空电磁系统中，垂直磁偶极子(水平线圈)是常用的发射源装置。设直角坐标系的原点位于垂直磁偶极正下方的地表，z坐标向下为正，磁偶极子发射源的高度为h。收发线圈为圆形中心回线式结构，接收观测装置位于发射线圈的中心，发射线圈半径为r，如图 1所示。
|
| ρl、hl分别为第l层的电阻率、厚度。 图 1 航空电磁探测层状大地模型 Fig. 1 Layered earth model for airborne electromagnetic exploration |
|
|
沿z方向的磁场垂直分量频率域表达式为[13]
(1) 其中：
式中：ω为角频率；μ0为自由空间磁导率；I(ω)为发射电流的频谱表达式；Z(j)、Zj分别为第j层的输入阻抗、固有阻抗；λ为积分变量；J0(x)、J1(x)分别为零阶、一阶贝塞尔函数；r0为接收线圈半径；l为地层数。
令发射电流为单位脉冲电流，即I(ω)=1，得到频率域电磁响应后，利用频-时转换公式即可得到时间域航空电磁系统大地单位脉冲响应：
(2) 式中：B(ω)为当发射电流为单位脉冲电流时的时间域航空电磁系统频率域响应；t为时间。
PRBS是具有一些随机序列特征的非随机序列[14]。m序列是一种常见的PRBS，由线性反馈移位寄存器生成，以m序列为源的电流由正负电平构成(图 2a)，主要参数包括阶数、码元频率、重复周期数，它的自相关函数是接近δ函数的三角脉冲[15](图 2b)。
|
| 图 2 4阶m序列的波形(a)及归一化自相关函数(b) Fig. 2 Waveform (a) and normalized autocorrelation function (b) of m-sequence generated by 4-rank shift register |
|
|
伪随机源电流产生的航空电磁响应根据全时正演模拟方法[1, 16]进行计算：
(3) 式中：*表示卷积；I(t)为伪随机源电流；BS(t)为时间域航空电磁系统大地单位阶跃响应。
2 伪随机源时间域航空电磁响应辨识 2.1 基本相关辨识方法相关辨识算法是系统辨识的经典算法，根据系统的输入信号和输出信号辨识得到系统的脉冲响应，要求辨识输入的信号具有随机性，因此具有一定的抗噪声干扰能力[17]。
典型的线性时不变系统模型如图 3所示。输入信号、输出信号和噪声信号间的数学关系为
(4)
|
| h(t). 单位脉冲响应；x(t). 输入信号；y(t). 输出信号；w(t).噪声信号。 图 3 线性时不变系统模型 Fig. 3 Linear time invariant system model |
|
|
对式(4)等号左右两边分别与x(t)进行相关运算得到
(5) 式中：Rxy(t)为输入信号和输出信号的互相关函数；Rxx(t)为输入信号的自相关函数；Rxw(t)为输入信号与噪声信号的互相关函数。
当输入信号为m序列时，可以近似认为m序列与噪声不相关，忽略Rxw(t)后得到
(6) 式(6)就是时间域相关辨识的核心公式——维纳-霍夫(Wiener-Hopf)方程。对式(6)进行反卷积运算就可以得到待辨识系统的单位脉冲响应。
2.2 时间域航空电磁响应辨识根据基本相关辨识方法，对时间域航空电磁响应的辨识可以将大地看作线性时不变系统。因此，忽略观测系统的影响后，观测到电压U(t)的表达式为
(7) 式中：S为线圈等效面积；B(t)为线圈接收到的大地电磁响应。
由于航空电磁系统只能观测到随时间变化的感应电磁场，无法观测到直流分量产生的响应，因此在相关辨识中需要去除m序列中的直流分量。去除直流分量后的自相关函数为
(8) 式中：RII(τ)为电流I(t)的自相关函数；τ为相关计算时移变量；a为电流幅值；N为m序列周期长度，N=2n-1(n为阶数)。
对式(7)两端分别与I(t)做相关运算并去除直流分量得到
(9) 式中：RIU(τ)为I(t)与U(t)的互相关函数；RIw(τ)为I(t)与w(t)的互相关函数。
计算相关函数时，若采用线性相关法，在移位过程中需要在序列的头部和尾部补零，计算结果在相关峰两侧将产生具有震荡特性的旁瓣[18]，影响最终辨识结果。而采用循环相关计算方法将序列首尾环状相连，则在计算中序列循环移位无需补零[19]，可以消除相关峰两侧的旁瓣。计算公式为
(10) 式中：T为U(t)的观测时间；Uc(t)为将U(t)首尾环状相连后的信号。T通常为m序列周期的倍数：
(11) 式中：fc为码元频率；q为重复周期数。
公式(9)反卷积的运算可以在时间域和频率域下进行。根据文献[10]，在频率域，当输出信号存在噪声干扰时，其转换的频率域信号会受到干扰，无法获得十分精确的值，且对于超出m序列频带宽度的噪声信号，噪声反而会被放大，因此无法抗噪。在时间域下，可以充分利用m序列与噪声不相关这一特点，实现对噪声的压制。
对于离散时间序列，将式(9)中电流自相关函数转换为自相关矩阵A，并设B为互相关序列、H为待求解的序列，根据维纳-霍夫方程，近似忽略噪声影响，则式(9)转换为矩阵形式：
(12) m序列的自相关矩阵为对称正定矩阵，因此式(12)可采用线性方程组的求解方法。本文采用共轭梯度法进行迭代求解。共轭梯度法是一种经典的优化算法，算法求解速度快且较为简单，步骤如下。
步骤1，给定H的迭代初始值H0和迭代精度ε，计算残量r0=B-A×H0，初始化p0=r0；
步骤2，迭代(k=0, 1, …)：
当Hk+1-Hk的标准差小于ε时迭代结束，记H=Hk+1。
获得反卷积运算结果后还需一步积分运算，可进一步压制随机噪声。积分运算采用矩形数值积分方法[20]。最终得到的大地脉冲响应辨识结果为
(13) 建立三层层状大地模型，根据伪随机源航空电磁正演进行数值模拟仿真，参数如表 1所示。
| ρ1/(Ω·m) | ρ2/(Ω·m) | ρ3/(Ω·m) | h1/m | h2/m | h3/m | r/m | h/m | S/m2 | f0/MHz | n | fc/kHz | a/A | q |
| 200 | 5 | 200 | 50 | 20 | ∞ | 3 | 30 | 1 | 1 | 6 | 10 | 1 | 10 |
| 注：f0为采样频率。 | |||||||||||||
不加入噪声，伪随机源电流信号如图 4a所示，通过式(7)合成电压信号如图 4b所示。采用循环相关计算的电流归一化自相关函数如图 4c所示，自相关峰两侧无旁瓣影响。作为对比，采用线性相关方法计算相关函数，其优势是通过快速傅里叶变换和逆变换可大大加快相关运算的速度，采用这种方法计算的电流归一化自相关函数如图 4d所示，尖峰两侧含有较大旁瓣。分别采用循环相关函数和线性相关函数进行辨识，辨识结果和理论正演响应如图 4e所示，计算相对误差如图 4f所示。相对误差的计算公式如下：
(14)
|
| a. 伪随机源电流信号；b. 观测到的电压信号；c. 循环相关计算的电流归一化自相关函数；d. 线性相关计算的电流归一化自相关函数；e. 两种辨识结果对比；f.两种辨识结果的相对误差。 图 4 数值模拟仿真数据及辨识结果 Fig. 4 Numerical simulation data and identification results |
|
|
式中，δ(n)为辨识结果相对误差的时间序列，其均方根为
(15) 式中，M为辨识结果相对误差的时间序列长度。
对比两种辨识的相对误差和相关计算耗费时间：使用循环相关函数进行辨识，辨识结果相对误差均方根为2.2×10-7，相关计算耗时为0.55 s；使用线性相关函数进行辨识，辨识结果相对误差均方根为2.4×105，相关计算耗时为0.01 s。因此，采用循环相关函数方法进行辨识的精度很高。本文提出的辨识方法是一种高精度辨识方法，尤其对于早期的辨识结果，其相对误差小于10-10，有利于准确提取浅层大地的响应，其缺点是运算速度较慢，对于庞大的航空电磁探测数据，更适用于部分含异常数据的高精度辨识。
3 航空电磁响应辨识压噪能力分析伪随机源时间域航空电磁响应高精度辨识方法的噪声压制能力来自于两部分：第一部分基于m序列与噪声不相关的这一特点实现噪声压制，主要体现在高精度辨识方法中的相关运算和反卷积运算；第二部分是积分运算中的多次累加起到了低通滤波的功能，进一步对噪声进行压制。
对噪声压制的能力以信噪比的提升作为评价依据，提升的信噪比计算公式如下：
(16) 式中：w为原始噪声；w′为经过辨识后的噪声；Var表示方差。SNRi是辨识算法第一部分提升的信噪比SNRi1和第二部分提升的信噪比SNRi2之和。
本文主要研究对随机噪声的压制能力，并采用正态分布随机数进行模拟[9]：
(17) 式中：wr为符合正态分布的随机数；μ为均值；σ为正态分布的标准差(本文表示噪声水平与观测电压值标准差之积)；Ni为0到1之间均匀分布的随机数；m为随机数的数目。
按照表 1中的参数进行数值模拟仿真，根据公式(7)，在观测电压信号中加入均值为0、噪声水平5%的随机噪声干扰，即令w=wr，辨识结果如图 5a所示，辨识结果的相对误差如图 5b所示。在时间域上对比原始噪声(图 5c)和经过辨识后的噪声(图 5d)，可以看出辨识后的噪声幅值明显减小。通过计算，提高的总信噪比为8.7 dB，说明这种高精度辨识算法对噪声具有较强的压制能力。
|
| a. 含噪辨识结果；b. 含噪辨识误差；c. 加入的原始噪声；d. 辨识后的噪声。 图 5 含噪时的辨识结果及噪声变化 Fig. 5 Identification results with noise and noise variation |
|
|
辨识对噪声的压制能力与伪随机序列的码元频率、阶数和重复周期数等参数有关[10, 21]，选择不同参数的m序列进行辨识，提升的信噪比如表 2所示。
| fc/kHz | n | q | T/ms | SNRi1/dB | SNRi2/dB | SNRi/dB |
| 10 | 6 | 1 | 6.300 | 6.2 | 2.5 | 8.7 |
| 10 | 6 | 2 | 12.600 | 8.1 | 3.3 | 11.4 |
| 10 | 6 | 4 | 25.200 | 10.4 | 2.8 | 13.6 |
| 10 | 6 | 6 | 37.800 | 12.2 | 2.9 | 15.1 |
| 10 | 6 | 8 | 50.400 | 13.6 | 2.7 | 16.3 |
| 10 | 6 | 10 | 63.000 | 14.2 | 2.7 | 16.9 |
| 10 | 7 | 1 | 12.700 | 7.9 | 3.3 | 11.2 |
| 10 | 8 | 1 | 25.500 | 11.0 | 2.8 | 13.8 |
| 5 | 6 | 1 | 12.600 | 7.6 | 3.2 | 10.8 |
| 20 | 7 | 1 | 6.350 | 6.0 | 3.2 | 9.2 |
| 40 | 8 | 1 | 6.375 | 5.1 | 3.8 | 8.9 |
根据表 2，高精度辨识算法第一部分提升的信噪比与相关运算的观测时间直接相关，当观测时间相同或相近时，尽管阶数、码元频率和重复周期数不相同，提升的信噪比也是相近的。
对SNRi1与T进行一元回归分析，建立多项式回归方程，采用最小二乘法计算最佳拟合模型，利用拟合优度R2评估模型性能，得到最佳的拟合回归线结果如图 6所示。
|
| 图 6 信噪比提升与观测时间的回归分析 Fig. 6 Regression analysis of SNRi and observed time |
|
|
由图 6可知，拟合优度接近于1，说明模型性能很好，可以根据观测时间对第一部分提升的信噪比进行较为准确的预测。
辨识算法第二部分提升的信噪比不受m序列参数影响，计算其平均值为3.0 dB。
综上，高精度辨识提升总信噪比的公式为
(18) 可据此对辨识算法能够提高的信噪比进行估计。
随着观测时间的延长，信噪比提升的增长幅度不断减小，并且存在极值，表明延长观测时间并不能一直提升信噪比，应合理选择伪随机序列的参数。
作为对比，采用双极性方波作为激励源的时间域航空电磁法常采用多周期线性叠加方法进行降噪，加入与上述仿真相同的噪声，并采用公式(16)计算提升的信噪比，结果如表 3所示。
| q | SNRi/dB |
| 1 | 6.1 |
| 2 | 12.1 |
| 4 | 18.1 |
| 6 | 21.6 |
| 8 | 24.1 |
| 10 | 26.0 |
对比表 2和表 3，由于单周期双极性方波的电磁响应可进行两次叠加，在多重复周期条件下，辨识提升的信噪比普遍低于线性叠加法。而在浅层高分辨的航空电磁探测中，多周期线性叠加降噪是不适用的，尤其对于早期的大地响应，往往需要采用单周期的数据进行处理。此时，根据表 2，单周期辨识提升的信噪比最小为8.7 dB，最大为13.8 dB，均大于双极性方波单周期叠加的6.1 dB，说明这种辨识方法对早期大地响应的压噪能力更好。
以上分析表明，本文提出的辨识方法具有较好的压噪能力，通过设置合适的伪随机序列参数，能够对浅层航空电磁探测中早期大地响应的噪声进行有效压制。
4 结论本文主要研究了伪随机源时间域航空电磁探测大地响应的辨识方法，数值模拟仿真和对比分析表明本文提出的辨识方法精度很高，是一种高精度辨识方法。进一步分析高精度辨识对噪声的压制能力，得出以下结论：
1) 本文采用循环相关计算相关函数方法消除了相关峰两侧的旁瓣干扰，并应用共轭梯度迭代法在时间域下进行反卷积运算，大大提高了大地脉冲响应的辨识精度，尤其早期大地响应辨识结果的相对误差小于10-10，非常有利于浅层地质体高分辨的探测。
2) 本文提出的辨识算法具有较好的噪声压制能力，且与相关运算的信号观测时间长度有关，对于实现浅层高分辨探测目的的时间域航空电磁法，在观测时间较短时，辨识算法较线性叠加法能提升更多的信噪比。
| [1] |
殷长春, 黄威, 贲放. 时间域航空电磁系统瞬变全时响应正演模拟[J]. 地球物理学报, 2013, 56(9): 3153-3162. Yin Changchun, Huang Wei, Ben Fang. The Full-Time Electromagnetic Modeling for Time-Domain Airborne Electromagnetic Systems[J]. Chinese Journal of Geophysics, 2013, 56(9): 3153-3162. |
| [2] |
Wright D, Ziolkowski A, Hall G. Improving Signal-to-Noise Ratio Using Pseudo Random Binary Sequences in Multi Transient Electromagnetic (MTEM) Data[C]//68th EAGE Conference and Exhibition. Vienna: EAGE, 2006: P065.
|
| [3] |
Ziolkowski A, Wright D. Comparison of PRBS and Square-Wave Transient CSEM Data over Peon Gas Discovery, Norway[J]. SEG Technical Program Expanded Abstracts, 2011, 30(1): 583-588. |
| [4] |
Ziolkowski A, Hobbs B, Wright D. Multitransient Electromagnetic Demonstration Survey in France[J]. Geophysics, 2007, 72(4): F197-F209. DOI:10.1190/1.2735802 |
| [5] |
汤井田, 罗维斌. 基于相关辨识的逆重复m序列伪随机电磁法[J]. 地球物理学报, 2008, 51(4): 1226-1233. Tang Jingtian, Luo Weibin. Pseudo-Random Electromagnetic Exploration Based on Invert-Repeated m-Sequence Correlation Identification[J]. Chinese Journal of Geophysics, 2008, 51(4): 1226-1233. DOI:10.3321/j.issn:0001-5733.2008.04.034 |
| [6] |
罗维斌, 李庆春, 汤井田. 编码电磁测深[J]. 地球物理学报, 2012, 55(1): 341-349. Luo Weibin, Li Qingchun, Tang Jingtian. Coded Source Electromagnetic Sounding Method[J]. Chinese Journal of Geophysics, 2012, 55(1): 341-349. DOI:10.6038/j.issn.0001-5733.2012.01.035 |
| [7] |
齐彦福, 殷长春, 王若, 等. 多通道瞬变电磁m序列全时正演模拟与反演[J]. 地球物理学报, 2015, 58(7): 2566-2577. Qi Yanfu, Yin Changchun, Wang Ruo, et al. Multi-Transient EM Full-Time Forward Modeling and Inversion of m-Sequences[J]. Chinese Journal of Geophysics, 2015, 58(7): 2566-2577. |
| [8] |
武欣, 薛国强, 底青云, 等. 伪随机编码源电磁响应的精细辨识[J]. 地球物理学报, 2015, 58(8): 2792-2802. Wu Xin, Xue Guoqiang, Di Qingyun, et al. Accurate Identification for the Electromagnetic Impulse Response of the Earth with Pseudo Random Coded Waveforms[J]. Chinese Journal of Geophysics, 2015, 58(8): 2792-2802. |
| [9] |
王若, 王妙月, 底青云, 等. 伪随机编码源激发下的时间域电磁信号合成[J]. 地球物理学报, 2016, 59(12): 4414-4423. Wang Ruo, Wang Miaoyue, Di Qingyun, et al. Signal Synthesis of Time-Domain Electromagnetic Response Excited by a Pseudo Random Binary Sequence Source[J]. Chinese Journal of Geophysics, 2016, 59(12): 4414-4423. DOI:10.6038/cjg20161203 |
| [10] |
王显祥, 底青云, 王妙月, 等. 基于m伪随机序列的电磁法抗噪能力分析[J]. 地球物理学报, 2016, 59(5): 1861-1874. Wang Xianxiang, Di Qingyun, Wang Miaoyue, et al. A Study on the Noise Immunity of Electromagnetic Methods Based on m Pseudo-Random Sequence[J]. Chinese Journal of Geophysics, 2016, 59(5): 1861-1874. |
| [11] |
欧阳涛, 底青云, 薛国强, 等. 利用多通道瞬变电磁法识别深部矿体: 以内蒙兴安盟铅锌银矿为例[J]. 地球物理学报, 2019, 62(5): 1981-1990. Ouyang Tao, Di Qingyun, Xue Guoqiang, et al. Identifying Deep Ore Bodies Using the Multi-Channel Transient Electromagnetic Method (MTEM): An Example of a Lead-Zinc-Silver Mine in Inner Mongolia[J]. Chinese Journal of Geophysics, 2019, 62(5): 1981-1990. |
| [12] |
陈超健, 蒋奇云, 莫丹, 等. 基于灰色判别准则和有理函数滤波的伪随机电磁数据去噪[J]. 地球物理学报, 2019, 62(10): 3854-3865. Chen Chaojian, Jiang Qiyun, Mo Dan, et al. De-Noising Pseudo-Random Electromagnetic Data Using Gray Judgment Criterion and Rational Function Filtering[J]. Chinese Journal of Geophysics, 2019, 62(10): 3854-3865. DOI:10.6038/cjg2019M0415 |
| [13] |
牛之琏. 时间域电磁法原理[M]. 长沙: 中南大学出版社, 2007. Niu Zhilian. Principle of Time-Domain Electromagnetic Method[M]. Changsha: Central South University Publishing House, 2007. |
| [14] |
淳少恒, 陈儒军, 耿明会. 伪随机m序列及其在电法勘探中的应用进展[J]. 地球物理学进展, 2014, 29(1): 439-446. Chun Shaoheng, Chen Rujun, Geng Minghui. Review of the Pseudo-Random m Sequence and Its Electrical Prospecting of Exploration Geophysics[J]. Progress in Geophysics, 2014, 29(1): 439-446. |
| [15] |
李白男. 伪随机信号及相关辨识[M]. 北京: 科学出版社, 1987. Li Bainan. Pseudo-Random Signal and Correlation Identification[M]. Beijing: Science Press, 1987. |
| [16] |
韩雪, 殷长春, 任秀艳, 等. 时间域航空电磁系统探测深度研究[J]. 吉林大学学报(地球科学版), 2019, 49(5): 1448-1456. Han Xue, Yin Changchun, Ren Xiuyan, et al. Research on Exploration Depth of Time-Domain Airborne EM System[J]. Journal of Jilin University (Earth Science Edition), 2019, 49(5): 1448-1456. |
| [17] |
甄晓丹. 相关辨识的算法优化[D]. 北京: 中国地质大学(北京), 2018. Zhen Xiaodan. The Algorithm Optimization of Correlation Identification[D]. Beijing: China University of Geosciences (Beijing), 2018. |
| [18] |
罗维斌, 李梅, 陈儒军. 相关辨识谱激电法及野外实验[J]. 地球物理学进展, 2017, 32(3): 1393-1398. Luo Weibin, Li Mei, Chen Rujun. Induced Polarization Method Based on Circular Cross-Correlation Method and Field Experiment[J]. Progress in Geophysics, 2017, 32(3): 1393-1398. |
| [19] |
Roinila T, Helin T, Vilkko M, et al. Circular Correlation Based Identification of Switching Power Converter with Uncertainty Analysis Using Fuzzy Density Approach[J]. Simulation Modelling Practice and Theory, 2009, 17(6): 1043-1058. DOI:10.1016/j.simpat.2009.03.006 |
| [20] |
裴易峰, 殷长春, 刘云鹤, 等. 时间域航空电磁磁场计算与应用[J]. 地球物理学进展, 2014, 29(5): 2191-2196. Pei Yifeng, Yin Changchun, Liu Yunhe, et al. Calculation and Application of B-Field for Time-Domain Airborne EM[J]. Progress in Geophysics, 2014, 29(5): 2191-2196. |
| [21] |
张建国, 武欣, 齐有政, 等. 时间域编码电磁勘探方法研究[J]. 雷达学报, 2014, 3(2): 158-165. Zhang Jianguo, Wu Xin, Qi Youzheng, et al. Investigation on Time Domain Coded Electromagnetic Exploration Method[J]. Journal of Radars, 2014, 3(2): 158-165. |