CN105652325A - 基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法 - Google Patents

Abstract

本发明涉及一种基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法。本发明针对地空测量的单点电磁数据，根据地空时域电磁数据近似e指数规律衰减的特征，按照对数等间隔确定时窗将数据进行分段，提取每时段内数据的时间常数值作为拟合参数，采用e指数拟合方法对每时段内的数据进行处理，将拟合输出结果作为预测值输入到滤波器中，再应用自适应标量卡尔曼滤波方法对数据进行电磁噪声滤除，将滤波后数据进行电阻率-深度成像。本发明的滤波方法与现有地空时域电磁数据滤波算法相比，不仅能够有效地抑制地空时域电磁数据中的电磁噪声，而且充分保留和增强了实测数据中的有效信息，提高了地空时域电磁数据的信噪比和质量。

Description

基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法

技术领域

本发明涉及一种时域电磁数据滤波方法，适用于符合e指数衰减规律的电磁数据滤波，尤其适用于地空时域电磁数据滤波。

背景技术

地空时域电磁法(Ground-airborneTime-domainelectromagneticsystem，简称GATEM)是利用地表铺设的长导线源发射双极性脉冲电流激励大地，在空间形成大面积的电磁波，在空中采用无人机或飞艇携带接收线圈对磁场响应进行采集。地空时域电磁法结合了地面瞬变电磁和航空时域电磁的优点。此方法不仅具有比地面瞬变电磁系统高效、应对复杂地形能力强等优势，而且比航空时域电磁系统空间分辨率更好、探测深度更大、信噪比更高、安全性更高。但是地空时域电磁系统在飞行测量过程中，受风向气流、飞行不平稳、传感器姿态等产生运动噪声，以及工频噪声、天电干扰等电磁噪声，导致地空电磁测量信号中含有较为复杂的噪声干扰，因此严重影响后期反演结果的准确性，所以对地空电磁探测数据的去噪方法研究具有重要意义。

目前已经应用到地空电磁数据噪声消除的方法有很多，其中Wang,Y.(2013)设计了一种基于综合小波变换的基线漂移校正方法，可以去除地空电磁测量数据的基线漂移。Bouchedda,A.(2010)提出天电噪声提取和小波系数叠加方法去除天电噪声，也就是应用小波变换能够有效地抑制天电噪声，但是当局部出现较大的天电干扰时，将改变衰减曲线特性。Chen,B.(2014)应用独立主成分分析法对于特征明显的工频强干扰去噪效果较好以及滤除一部分人文噪声。

CN201110197887.4公开了长导线源瞬变电磁地空探测方法，采用长导线源瞬变电磁地空探测装置对深部地质目标体进行探测，获得瞬变电磁的观测数据，将其转换成瞬变电磁虚拟波数据，采用多点数据合成获得瞬变电磁合成孔径数据体，并对合成孔径数据体进行处理解释，获得深部地质目标体信息。

CN201510193706.9公开了基于数学形态学的航空电磁数据去噪方法，提出自适应多尺度复合形态滤波方法，根据信号局部特征和噪声特点自适应的选择结构元素类型和尺寸大小，滤除信号中的正、负脉冲噪声和随机噪声。

CN201110197887.4公开了长导线源瞬变电磁地空探测方法以及处理解释方法，但是并没有给出对地空探测数据的滤波方法；CN201510193706.9公开了一种航空去噪方法，针对航空电磁信号特点进行滤波，虽然可以采用数学形态学这种方法对地空电磁探测数据进行去噪处理，去除信号中的正、负脉冲噪声和随机噪声，但是实际地空电磁探测数据中同时含有白噪声、工频干扰、发动机、运动等电磁噪声，混合存在很难进行分离。目前针对地空电磁探测数据的去噪方法大多针对某一类电磁噪声能够进行有效地去除，但是未考虑地空电磁信号特点，导致有较大噪声干扰时衰减曲线特性被改变。因此需要设计一种考虑地空电磁信号特征并且能够有效滤除多种噪声的滤波方法。

发明内容

本发明的目的在于针对现有去噪方法的不足，根据地空电磁数据衰减特征，提供一种能够同时有效滤除多种电磁噪声，基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法。

本发明的主要思想是：在提取地空电磁数据衰减特征的基础上，建立指数拟合预报方程，将拟合输出结果作为预测值，将拟合输出结果作为预测值输入到滤波器中，滤除实测数据的电磁噪声，提取有效信号，将滤波后地空电磁数据进行电阻率-深度参数计算、成像。

本发明是这样实现的，一种基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法包括：

1)首先从实测连续数据中进行分割，获得单点的地空电磁测量数据；

2)利用加权均值滤波方法对电磁数据进行预处理；

3)提取地空电磁数据的衰减特征，建立指数拟合预报方程，对电磁数据进行e指数拟合；

4)建立自适应标量卡尔曼滤波器，将拟合输出结果作为预测值输入到滤波器中，滤除实测数据的电磁噪声，提取有效信号；

5)将滤波后地空电磁数据进行电阻率-深度参数计算、成像。

进一步地，步骤3中，按照对数等间隔方式将数据进行分段，提取每时段数据的时间常数τ值作为拟合参数，采用e(自然对数)指数拟合方法对每时段数据进行拟合。拟合方程为：

$v_k=a_1+a_2\×e^{-a_3{t}_k}---\left(1\right)$

其中，v_k为第k道的电磁响应预测值，a₁为调整系数，a₂为系统增益，a₃为时间常数τ，t_k为第k道的延迟时间。

进一步地，步骤4中包含以下步骤：

Ⅰ、从实测连续数据中，通过计算噪声均值，提取测量系统的电磁噪声统计特性；

Ⅱ、输入参数p_k，b_k和R_k初始值，其中p_k为误差协方差，b_k遗忘因子,R_k为测量噪声的方差；

Ⅲ、计算加权系数估计值其中v_k为步骤3中第k取样道电磁数据的拟合结果；

Ⅳ、计算自适应标量卡尔曼滤波器的基本递归表达式参数；

Ⅴ、计算时变噪声统计估计量；

Ⅵ、对单点地空电磁测量数据的所有时间道重复Ⅲ-Ⅴ步骤，并判断是否完成全部取样道，如果未完成，则计算遗忘因子b_k，返回步骤Ⅲ。

步骤4所述的的自适应标量卡尔曼滤波器，其中自适应标量卡尔曼滤波器的基本递归表达式：

k_k+1＝α_kp_k(cα_kp_k+R_k+1)^-1(2)

p_k+1＝(1-k_k+1c)α_kp_k(3)

${\hat{x}}_{k+1}={\mathrm{\α}}_k{\hat{x}}_k+k_{k+1}(y_{k+1}-{\mathrm{c\α}}_k{\hat{x}}_k)---\left(4\right)$

其中k_k+1为卡尔曼增益，c为信号估计系数，y_k+1为第k取样道的地空测量值，为滤波结果。

进一步地，步骤5中包含以下步骤：

a.输入测量点坐标；

b.计算收发距r，确定计算长导线源时所需剖分数N；

c.计算每个电偶极子响应的核函数；

d.计算N个电偶极子的总核函数；

e.对总核函数f(θ)进行最大值归一化处理，其中μ₀为真空中磁场导率，σ为地下电阻率，t为时间自变量；

f.实测数据计算的核函数f(θ)是否为最大，确定电阻率值；

g.计算深度参数，进行电阻率-深度成像。

有益效果：本发明与现有技术相比，不仅能够有效地抑制实际地空飞行探测数据中同时含有的白噪声、工频干扰、天电噪声等电磁感应混合噪声，而且充分保留和增强了实测数据中的有效信息，提高了地空时域电磁数据的信噪比和质量，完全满足地空时域电磁要求。本方法为我国特殊地区开展地空飞行电磁探测寻找矿产资源提供新的技术保障，有利于地空时域电磁探测技术的实用化。

附图说明

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

图2是基于指数拟合-自适应卡尔曼滤波方法流程图；

图3是指数拟合方法流程图；

图4是自适应标量卡尔曼滤波方法流程图；

图5是电磁数据计算电阻率-深度参数方法示意图；

图6是本发明一个实施例的指数拟合效果图；

图7是本发明一个实施例的地空时域电磁数据滤波效果图；

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合吉林省公主岭市怀德镇地空时域电磁探测为例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例

图1是地空时域电磁探测系统示意图，在吉林省公主岭市怀德镇地表铺设的长导线源发射双极性脉冲电流，在空中采用滑翔机携带接收线圈对磁场响应进行采集。基于指数拟合-自适应卡尔曼对地空电磁测量数据进行滤波处理。

参见图2结合图1所示，一种基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法，包括：

1)首先从实测连续数据中进行分割，获得单点的地空电磁测量数据；

从滑翔机携带的接收系统接收的连续地空电磁测量数据中，提取出接收信号中的衰减曲线，获得单点的地空电磁测量数据及其对应坐标，并对单点的地空电磁测量数据进行以下处理。

2)利用加权均值滤波方法对电磁数据进行预处理；

为了提高拟合精度，对单点的地空电磁测量数据先进行加权均值滤波，公式为：

$y_k=\frac{y_{k-1}+2y_k+y_{k+1}}{4}---\left(1\right)$

其中，y_k为第k道的地空测量值。

3)提取地空电磁数据的衰减特征，建立指数拟合预报方程，对电磁数据进行e指数拟合；

如图3所示的步骤3中，由于电磁信号在早期衰减快，晚期衰减慢，因此，在进行数据拟合时，按照对数等间隔方式将数据进行分段，提取每时段数据的时间常数τ值作为拟合参数。采用e(自然对数)指数拟合方法对每时段数据进行拟合，拟合方程为：

$v_k=a_1+a_2\×e^{-a_3{t}_k}---\left(2\right)$

其中，v_k为第k道的电磁响应预测值，a₁为调整系数，a₂为系统增益，a₃为时间常数τ，t_k为第k道的延迟时间。

4)建立自适应标量卡尔曼滤波器，将拟合输出结果作为预测值输入到滤波器中，滤除实测数据的电磁噪声，提取有效信号；

如图4所示的步骤4中包含以下步骤：

Ⅰ、从实测连续数据中，通过计算噪声均值，提取测量系统的电磁噪声统计特性；

地空时域电磁测量过程中，其数学模型可以由如下测量方程给出：

x_k＝αx_k-1+w_k-1(3)

y_k＝cx_k+η_k(4)

其中x_k为k时刻的估计信号；y_k为k时刻的测量信号；w_k-1为测量噪声；η_k为系统噪声，α、c为信号估计系数，与信号本身有关；系统噪声和测量噪声的统计特性：

E[η_kη_l]＝R_kδ_kl(5)

其中R_k为测量噪声的方差，

Ⅱ、根据电性源地空数据特点，以及系统噪声和测量噪声的统计特性，给出参数p_k，b_k和R_k初始值，其中p_k为误差协方差，b_k遗忘因子。

Ⅲ、计算加权系数估计值其中v_k为步骤3中第k取样道电磁数据的拟合结果；

Ⅳ、计算自适应标量卡尔曼滤波器的基本递归表达式参数；

自适应标量卡尔曼滤波器的基本递归表达式：

k_k+1＝α_kp_k(cα_kp_k+R_k+1)^-1(6)

p_k+1＝(1-k_k+1c)α_kp_k(7)

${\hat{x}}_{k+1}={\mathrm{\α}}_k{\hat{x}}_k+k_{k+1}(y_{k+1}-{\mathrm{c\α}}_k{\hat{x}}_k)---\left(8\right)$

其中k_k+1为卡尔曼增益，y_k+1为第k道的地空测量值，为滤波结果；

Ⅴ、计算时变噪声统计估计量；

时变噪声统计估计量通过下式计算：

$R_{k+1}=(1-d_k)R_k+d_k(y_{k+1}-{\mathrm{c\α}}_k{\hat{x}}_k-c_k{P}_{k+1,k})---\left(9\right)$

其中b_k为遗忘因子

Ⅵ、对单点地空电磁测量数据的所有时间道重复Ⅲ-Ⅴ步骤，并判断是否完成全部取样道，如果未完成，则计算遗忘因子b_k，返回步骤Ⅲ。

5)将滤波后地空电磁数据进行电阻率-深度参数计算、成像。

参见图5所示并结合图2，步骤5中感应电动势计算电阻率方法。长导线源激励响应可以采用多个偶极源响应的叠加等效，根据电偶极子的均匀半空间的解析解，设计了求解全区视电阻率算法。包含以下步骤：

a.输入测量点坐标；

b.计算收发距r，确定计算长导线源时所需剖分数N；

c.计算每个电偶极子响应的核函数；

在地面接收电磁响应时，电偶极源时域感应电动势写为：

$V_z\left(t\right)=\frac{3{\mathrm{IS\μ}}_{0}dl}{8\mathrm{\π}t}\frac{y}{{\mathrm{\θ}}^2{r}^5}\[erf\left(\mathrm{\θ}r\right)-\frac{2}{\sqrt{\mathrm{\π}}}\mathrm{\θ}r(1+\frac{2}{3}{\mathrm{\θ}}^2{r}^2)e^{-{\mathrm{\θ}}^2{r}^2}\]---\left(10\right)$

式中S为接收线圈有效面积，I为发射电流，t为时间自变量，dl为电偶极子长度，r为收发距，x、y、z为测量点坐标，其中σ为地下电阻率，μ₀为真空中磁场导率，erf为高斯误差函数。

d.计算N个电偶极子的总核函数；

长导线源感应电动势响应写为：

$H_c=\underset{i=1}{\overset{N}{\Σ}}{H}_{di}---\left(11\right)$

式中H_c为长导线源电磁响应，H_di为第i个电偶极源的电磁响应。把(10)式代入到(11)式中，并化简可得核函数：

$f\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{N}{\Σ}}\frac{y}{{\mathrm{\θ}}^2{r_i}^5}\[erf\left({\mathrm{\θr}}_i\right)-\frac{2}{\sqrt{\mathrm{\π}}}{\mathrm{\θr}}_i(1+\frac{2}{3}{\mathrm{\θ}}^2{r_i}^2)e^{-{\mathrm{\θ}}^2{r_i}^2}\]=\frac{8\mathrm{\π}V\left(t\right)}{3{\mathrm{IS\μ}}_0{\mathrm{dl}}_i}---\left(12\right)$

式中N为长导线源剖分个数，dl_i为第i个电偶极子长度，L为长导线源半长度，式中V(t)为测量感应电动势，

e.对总核函数f(θ)进行最大值归一化处理；

f.实测数据计算的核函数f(θ)是否为最大，确定电阻率值；

通过对核函数进行求解θ值来反求视电阻率为了保证解存在，在核函数中引入校正系数，采用核函数的最大值进行归一处理，根据核函数曲线的单调性进行划分为早期和晚期曲线，在每一段曲线上采用二分法进行搜索求解θ值计算电阻率，当求出两个电阻率值，需要和早期、晚期的近似解进行判断，最后确定视电阻率值。

g.计算深度参数，进行电阻率-深度成像。

根据步骤f计算的视电阻率值，计算第i时间道所在的深度

图6为采用图3所示的本发明一个实施例地空理论模型衰减曲线的指数拟合效果图；

图7为本发明一个实施例的地空时域电磁数据滤波后的结果，滤波后的衰减曲线信噪比明显提高，充分验证了指数拟合-自适应卡尔曼滤波方法的有效性，为地空时域电磁测量数据高精度解释提供了新的思路和方法。

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

Claims (5)

1.一种基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法，其特征在于，包括如下步骤：

1)首先从实测连续数据中进行分割，获得单点的地空电磁测量数据；

2)利用加权均值滤波方法对电磁数据进行预处理；

3)提取地空电磁数据的衰减特征，建立指数拟合预报方程，对电磁数据进行e指数拟合；

4)建立自适应标量卡尔曼滤波器，将拟合输出结果作为预测值输入到滤波器中，滤除实测数据的电磁噪声，提取有效信号；

5)将滤波后地空电磁数据进行电阻率-深度参数计算、成像。

2.按照权利要求1所述的一种基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法，其特征在于：

步骤3中，按照对数等间隔方式将数据进行分段，提取每时段数据的时间常数τ值作为拟合参数，采用e(自然对数)指数拟合方法对每时段数据进行拟合。

3.按照权利要求1所述的一种基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法，其特征在于，步骤4中包含以下步骤：

Ⅰ、从实测连续数据中，通过计算噪声均值，提取测量系统的电磁噪声统计特性；

Ⅱ、输入参数p_k，b_k和R_k初始值，其中p_k为误差协方差，b_k遗忘因子,R_k为测量噪声的方差；

Ⅲ、计算加权系数估计值其中v_k为步骤3中第k取样道电磁数据的拟合结果；

Ⅳ、计算自适应标量卡尔曼滤波器的基本递归表达式参数；

Ⅴ、计算时变噪声统计估计量；

Ⅵ、对单点地空电磁测量数据的所有时间道重复Ⅲ-Ⅴ步骤，并判断是否完成全部取样道，如果未完成，则计算遗忘因子b_k，返回步骤Ⅲ。

4.按照权利要求3所述的自适应标量卡尔曼滤波器，其特征在于，自适应标量卡尔曼滤波器的基本递归表达式：

k_k+1＝α_kp_k(cα_kp_k+R_k+1)^-1(1)

p_k+1＝(1-k_k+1c)α_kp_k(2)

${\hat{x}}_{k+1}={\mathrm{\α}}_k{\hat{x}}_k+k_{k+1}(y_{k+1}-{\mathrm{c\α}}_k{\hat{x}}_k)---\left(3\right)$

其中k_k+1为卡尔曼增益，c为信号估计系数，y_k+1为第k取样道的地空测量值，为滤波结果。

5.按照权利要求1所述的一种基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法，其特征在于，步骤5中包含以下步骤：

a.输入测量点坐标；

b.计算收发距r，确定计算长导线源时所需剖分数N；

c.计算每个电偶极子响应的核函数；

d.计算N个电偶极子的总核函数；

e.对总核函数f(θ)进行最大值归一化处理，其中μ₀为真空中磁场导率，σ为地下电阻率，t为时间自变量；

f.实测数据计算的核函f(θ)是否为最大，确定电阻率值；

g.计算深度参数，进行电阻率-深度成像。