地空时域电磁系统浅部低阻异常体高精度数据解释方法与流程

本发明属于地球物理勘探技术领域，具体涉及一种时域电磁系统浅部低阻异常体高精度数据解释方法，尤其适用于地空时间域电磁系统探测方法。

背景技术：

长导线源地空时域电磁系统在地面铺设长导线作为发射源，将便携接收设备搭载在旋翼无人机、飞艇等平台上作为空中接收装置，可实现滩涂、山区、海路交互带等复杂区域的高效、快速地质勘察任务。它相对于航空时域电磁探测系统具有勘探深度大、飞行成本低、安全性高等优势，近年来已逐渐成为时域电磁探测系统的研究热点，被广泛的应用在生产实践中。

为了实现对低阻异常体的识别，长导线源时域电磁系统在实际应用中需要对数据进行解释处理。长导线源时域电磁解释方法通过计算视电阻率与视深度，可实现浅层以及深层异常体的数据解释与成像，它具有简单、高效的特点，因此被广泛使用。但是该方法局限于对地面数据的高精度解释，在处理过程中忽略了飞行高度对解释带来的影响，这一影响对于浅部异常体尤其明显。

中国专利201510039201.7公开了一种频率域地空电磁勘探方法，该方法采用地面发射，空中接收电磁波信号的工作模式，接收系统搭载在飞行器上，通过全区视电阻率法反演解释地下电性结构，是一种新型的电磁勘探方法。在测区上空测量磁场，能够适应地表结构复杂环境的同时减弱了近场影响引起的静态效应，拓展了电磁勘探的探测范围。系统可在测量多个磁场分量的情况下对被测磁场分量进行校正和补偿，提高了磁场测量的分辨能力。

中国专利201510916688.2公开了一种地面参考区磁场延拓的地空协同电磁数据校正方法。采用低温超导磁传感器在地面参考区进行磁场测量后，采用fft方法对磁场进行向上延拓、滤波取样，以获得参考区在空中飞行高度下的磁场测量基准值，再将参考区的实际飞行测量数据进行基线校正、滤波、叠加取样、电阻率-深度成像。采用svd奇异值分解方法实现磁场实测数据与基准值的拟合，确定测量系统的固有误差、基线漂移量、电阻率-深度参数误差等，实现地空电磁数据的高精度测量。

中国专利201410081433.4公开了一种长导线源瞬变电磁地空探测方法，将观测数据转换成瞬变电磁虚拟波数据，采用多点数据合成的方法，获得瞬变电磁合成孔径数据体，完成对深部地质目标体的精细解释的数据储备；采用克希霍夫偏移成像方法及速度建模方法，获得瞬变电磁合成孔径成像剖面，并对瞬变电磁合成孔径数据体进行处理及解释，完成对深部地质目标体的精细探测，获得深部地质目标体的信息。

上述三个专利，其中专利201510039201.7公布了频率域地空电磁勘探方法，它采用全区视电阻率法反演解释地下电性结构。由于频率域电磁法的主要关注点集中在深部异常体的探测，因此并没有提及空中接收装置距离地面的高度对磁场测量分辨能力的影响。

专利201510916688.2公布了地面参考区磁场延拓的地空协同电磁数据校正方法，为了避免fft向下延拓方法的不稳定性，该方法的本质是采用地面参考区测得的磁场进行向上延拓，以实现空中磁场实测数据的数据校正，实现地空电磁数据的高精度测量。

专利201410081433.4发布了一种长导线源瞬变电磁地空探测方法，它采用瞬变电磁合成孔径方法对数据体进行处理和解释，避免了飞行高度的影响，完成了深部地质目标体的精细探测。

技术实现要素：

为了提高地空时域电磁探测系统对浅部低阻异常体的数据解释精度，本专利提出一种时域电磁系统浅部低阻异常体高精度数据解释方法，通过浅部拟二维异常体长导线源的时域电磁响应数值模拟、随机噪声加载、基于正则化方法的含噪电磁响应的向下延拓、噪声抑制，视电阻率-视深度的计算，提高数据的解释及成像精度。

本发明是这样实现的，提供一种地空时域电磁系统浅部低阻异常体高精度数据解释方法，所述时域电磁系统为长导线源时域电磁系统，包括如下步骤：

1)基于长导线源时域电磁系统的浅部拟二维低阻异常体电磁响应数值模拟；

2)对步骤1)的数值模拟结果进行基于不同信噪比的随机噪声加载；

3)采用正则化方法实现步骤2)中含噪信号的长导线源时间域电磁响应的向下延拓；

4)通过调节正则化参数确定不同信噪比条件下正则化延拓方法的最优精度；

5)采用自适应卡尔曼滤波算法对步骤4)中不符合容差条件的电磁响应进行消噪处理；

6)计算视电阻率，视深度，并完成视电阻率-视深度成像。

进一步地，上述方法具体步骤如下：

步骤1)中，采用层状介质模型电磁响应表达式对浅部拟二维低阻异常体电磁响应进行数值模拟：通过测点下方的大地模型确定数据模拟时的大地层数、各层电导率及厚度值，接地长导线源n层大地模型的z方向磁场频率域响应表达式为：

其中，l为接地导线的半长度，i为发射电流，x为观测点的x坐标，y为观测点的y坐标，z为观测点的z坐标，r＝[(x-x′)²+y²]^1/2，λ、x′均为被积变量，j1为贝塞尔函数，为反射系数，i²＝-1，ω为角频率，n为大地模型的层数，hn为每一层的厚度，σn为电导率，μ0为真空介质的磁导率，μn为大地各层模型磁导率，当层数不为1时，可通过逐层迭代的方式最终推出地表的反射系数；

步骤2)中，对测线上不同测点的长导线源时域电磁响应数值模拟结果分别进行随机噪声的加载，所加载的噪声为仪器本底噪声和高斯白噪声，它们的平均值分别为40db、50db、60db、70db、80db和90db，具体计算式如下：

其中，k为时间道，ntime为总的时间道数量，vz为长导线源时域电磁响应，vn为噪声信号；

步骤3)中，采用正则化方法实现步骤2)中含噪信号的长导线源时域电磁响应数值的向下延拓，正则化算子在低频时近似原始向下延拓算子，能够对高频噪声起到抑制作用，基于正则化方法的向下延拓计算式如下：

vz＝0(kx,ky,z＝0)＝htik(kx,ky)·vz＝h(kx,ky,z＝h)(3)

其中vz＝h(kx,ky,z＝h)为空中长导线源时域电磁响应数值模拟结果经二维傅里叶变换后的值，vz＝0(kx,ky,z＝0)为向下延拓后的地面电磁响应波数域值，htik(kx,ky)为正则化向下延拓算子，其表达式为kx、ky为x、y方向的波数，α为正则化参数；

步骤4)中，采用步骤3)中的方法对信噪比分别为40db、50db、60db、70db、80db和90db的含噪信号进行向下延拓，通过调节正则化参数确定出不同信噪比条件下的最优延拓精度，如表1所示；

表1、不同信噪比下的最优延拓精度与正则化参数

步骤5)中，采用自适应卡尔曼滤波算法对步骤4)中不符合容差条件的电磁响应数值模拟结果进行消噪处理，自适应卡尔曼滤波算法的具体步骤如下：

a、计算接收机本底噪声和高斯白噪声的均值qk、rk与方差qk、rk，k为时间道；

b、根据本底噪声和高斯白噪声的均值与方差，结合接地长导线源地空时域电磁数据的特点确定参数pk、bk、qk、rk、qk和rk的初始值，其中pk为误差协方差，bk为遗忘因子；

c、构造自适应卡尔曼滤波器的基本递归表达式：

αk＝α(1-lk+1)(4)

pk+1＝(1-lk+1c)pk+1,k(7)

其中αk为加权系数估计值，pk+1,k为预测协方差，lk+1为卡尔曼增益，c为信号估计参数，ek+1为预测误差，为滤波结果，yk+1为k+1时刻的实测信号；

d、噪声的自适应统计估计量-均值与方差为：

其中

e、判断是否完成全部时间道的计算，

步骤6)中，均匀半空间模型的长导线源时域电磁响应表达式如下所示：

其中收发距i为发射电流，dl为电偶极子长度，σ为电导率，μ为磁导率，t为采样时间，s为线圈面积，ncf为长导线源所拆分的电偶极子的数量；

由于公式(14)非常复杂，难以从中求解出对应的θ值，通过给定电导率σ的取值范围为0.001s/m-10s/m以及时间t的取值范围10^-5s-10^-2s，计算出对应的一组值，再将其带回到公式(14)中求出对应于不同θ的vz均匀，最后采用步骤1)中计算出的z＝0平面的理论感应电动势vz值与这一组vz均匀值逐个进行比较，将拟合程度最好的vz均匀所对应的θ作为所求得的值，再将它带入到视电阻率ρ和θ的关系式中，即可计算出视电阻率值：

将视电阻率带入到视深度计算式中：

最终通过计算出的视电阻率和视深度即可实现地空时域电磁数据的解释与成像。

进一步地，表达式(1)中的积分采用汉克尔积分变换算法进行求解，之后再采用guptasarm数字滤波方法将其从频率域转换到时间域，即可得到z＝0平面理论上的时域电磁响应vz。

进一步地，步骤3)中，正则化参数α的取值范围为0.1-10，取值越小，转折频率越低，对高频噪声的抑制能力越强，但同时也消除了更多的有效信号；取值越大，转折频率越高，对有效信号的保留越多，但同时降低了高频噪声的抑制能力；实际中根据信号所含噪声具体情况选择正则化参数：噪声大时减小正则化参数，噪声小时增大正则化参数，参数的具体选取方法采用对向下延拓结果再进行向上延拓，求出误差，得到最优的正则化参数值。

与现有技术相比，本发明的优点在于：

本发明针对采用视电阻率-视深度数据解释方法对地空时域电磁探测系统的浅部异常体进行成像时精度低的问题，并且当信号中存在高频噪声时，采用迭代法进行向下延拓，随着延拓次数的增加，延拓精度不仅不会降低还会逐渐增加，直至不收敛，本专利提出采用正则化方法实现地空时域电磁信号的向下延拓，正则化参数的具体选取方法可以通过对向下延拓后的结果再进行向上延拓，求出误差，不断尝试以得到最优的参数值，而非l曲线选取法，本方法可避免飞行高度的影响，实现浅部异常体的含噪信号时域电磁响应的高精度数据解释与成像。

附图说明

图1是地空长导线源时域电磁探测系统示意图；

图2是基于地空时域电磁探测系统的浅部异常体高精度数据解释方法流程图；

图3是浅部拟二维低阻异常体模型示意图；

图4是正则化向下延拓方法及参数确定方法流程图；

图5是视电阻率、视深度计算过程流程图；

图6是低阻浅部异常体视电阻率-视深度成像图，其中图6(a)为理论模型示意图、图6(b)为向下延拓后数据解释结果、图6(c)为空中电磁响应数据解释结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

实施例、

参考图2结合图1，本发明提供一种地空时域电磁系统浅部低阻异常体高精度数据解释方法，地空时域电磁系统为长导线源时域电磁系统，包括如下步骤：

1)基于长导线源时域电磁系统的浅部拟二维低阻异常体电磁响应数值模拟；

采用层状介质模型电磁响应表达式对浅部拟二维低阻异常体电磁响应进行数值模拟：参考图3(σ′2为低阻异常体电导率)，通过测点下方的大地模型确定数据模拟时的大地层数、各层电导率及厚度值，接地长导线源n层大地模型的z方向磁场频率域响应表达式为：

其中，l为接地导线的半长度，i为发射电流，x为观测点的x坐标，y为观测点的y坐标，z为观测点的z坐标，r＝[(x-x′)²+y²]^1/2，λ、x′均为被积变量，j1为贝塞尔函数，为反射系数，i²＝-1，ω为角频率，n为大地模型的层数，hn为每一层的厚度，σn为电导率，μ0为真空介质的磁导率，μn为大地各层模型磁导率，当层数不为1时，可通过逐层迭代的方式最终推出地表的反射系数；

表达式(1)中的积分采用汉克尔积分变换算法进行求解，之后再采用guptasarm数字滤波方法将其从频率域转换到时间域，即可得到z＝0平面理论上的时域电磁响应vz。

2)对步骤1)的数值模拟结果进行基于不同信噪比的随机噪声加载；

对测线上不同测点的长导线源时域电磁响应数值模拟结果分别进行随机噪声的加载，所加载的噪声为仪器本底噪声和高斯白噪声，它们的平均值分别为40db、50db、60db、70db、80db和90db，具体计算式如下：

其中，k为时间道，ntime为总的时间道数量，vz为长导线源时域电磁响应，vn为噪声信号；

3)采用正则化方法实现步骤2)中含噪信号的长导线源时间域电磁响应的向下延拓；

如图4所示，采用正则化方法实现步骤2)中含噪信号的长导线源时域电磁响应数值的向下延拓，正则化算子在低频时近似原始向下延拓算子，能够对高频噪声起到抑制作用，基于正则化方法的向下延拓计算式如下：

vz＝0(kx,ky,z＝0)＝htik(kx,ky)·vz＝h(kx,ky,z＝h)(3)

其中vz＝h(kx,ky,z＝h)为空中长导线源时域电磁响应数值模拟结果经二维傅里叶变换后的值，vz＝0(kx,ky,z＝0)为向下延拓后的地面电磁响应波数域值，htik(kx,ky)为正则化向下延拓算子，其表达式为kx、ky为x、y方向的波数，α为正则化参数；

正则化参数α的取值范围为0.1-10，取值越小，转折频率越低，对高频噪声的抑制能力越强，但同时也消除了更多的有效信号；取值越大，转折频率越高，对有效信号的保留越多，但同时降低了高频噪声的抑制能力；实际中根据信号所含噪声具体情况选择正则化参数：噪声大时减小正则化参数，噪声小时增大正则化参数，参数的具体选取方法采用对向下延拓结果再进行向上延拓，求出误差，得到最优的正则化参数值。

4)通过调节正则化参数确定不同信噪比条件下正则化延拓方法的最优精度；

采用步骤3)中的方法对信噪比分别为40db、50db、60db、70db、80db和90db的含噪信号进行向下延拓，通过调节正则化参数确定出不同信噪比条件下的最优延拓精度，如表1所示；

表1、不同信噪比下的最优延拓精度与正则化参数

5)采用自适应卡尔曼滤波算法对步骤4)中不符合容差条件的电磁响应进行消噪处理；

采用自适应卡尔曼滤波算法对步骤4)中不符合容差条件的电磁响应数值模拟结果进行消噪处理，自适应卡尔曼滤波算法的具体步骤如下：

a、计算接收机本底噪声和高斯白噪声的均值qk、rk与方差qk、rk，k为时间道；

b、根据本底噪声和高斯白噪声的均值与方差，结合接地长导线源地空时域电磁数据的特点确定参数pk、bk、qk、rk、qk和rk的初始值，其中pk为误差协方差，bk为遗忘因子；

c、构造自适应卡尔曼滤波器的基本递归表达式：

αk＝α(1-lk+1)(4)

pk+1＝(1-lk+1c)pk+1,k(7)

其中αk为加权系数估计值，pk+1,k为预测协方差，lk+1为卡尔曼增益，c为信号估计参数，ek+1为预测误差，为滤波结果，yk+1为k+1时刻的实测信号；

d、噪声的自适应统计估计量-均值与方差为：

其中

e、判断是否完成全部时间道的计算，

6)计算视电阻率，视深度，并完成视电阻率-视深度成像；

参见图5，均匀半空间模型的长导线源时域电磁响应表达式如下所示：

其中收发距i为发射电流，dl为电偶极子长度，σ为电导率，μ为磁导率，t为采样时间，s为线圈面积，ncf为长导线源所拆分的电偶极子的数量；

由于公式(14)非常复杂，难以从中求解出对应的θ值，通过给定电导率σ的取值范围为0.001s/m-10s/m以及时间t的取值范围10^-5s-10^-2s，计算出对应的一组值，再将其带回到公式(14)中求出对应于不同θ的vz均匀，最后采用步骤1)中计算出的z＝0平面的理论感应电动势vz值与这一组vz均匀值逐个进行比较，将拟合程度最好的vz均匀所对应的θ作为所求得的值，再将它带入到视电阻率ρ和θ的关系式中，即可计算出视电阻率值：

将视电阻率带入到视深度计算式中：

最终通过计算出的视电阻率和视深度即可实现地空时域电磁数据的解释与成像。

如图6所示，图6(a)为理论模型示意图、图6(b)为向下延拓后数据解释结果、图6(c)为空中电磁响应数据解释结果。从图中可以看出，解释结果均识别出了低阻层的变化趋势。然而，空中电磁响应数据的解释结果与理论模型存在较大偏差，向下延拓后数据的解释结果和理论模型更加接近，提高了浅部低阻异常体的识别分辨率。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。