CN111679326B - 基于地质探测的瞬态宽带电磁探测信号处理方法 - Google Patents

Abstract

一种基于地质探测的宽带电磁信号（20Hz‑150kHz）处理技术，包括：在原始信号添加白噪声，归一化处理，得到新信号；根据带噪信号的包络特征，应用自动移位法、波形匹配法或比例延拓法做经验模式分解，得到IMF分量和残余分量，其中自动移位法首先移出变化最快的IMF分量，并且随着IMF级别的增加，可使宽带电磁信号的IMF分量快速平滑；消除IMF分量中的白噪声，获得真实的IMF分量；计算各IMF的时域一阶差分能量、相位一阶差分能量和归一化能量占比，选取三种特征均最大的IMF，即可作为辐射信号的主要特征；将选取出的IMF频率与特定矿物质所辐射的电磁波频率对比，若一致则说明矿石种含有该矿物质元素。

Description

基于地质探测的瞬态宽带电磁探测信号处理方法

技术领域

本发明涉及地质探测技术领域，具体涉及地壳中矿物质所产生的宽带磁场辐射信号提供数据处理方法。

背景技术

地下矿物质磁场信号包含矿物分布和地质结构等重要特征信息，可以为地质勘测提供重要参考保障，因此，其检测各种矿物质磁场辐射信号的技术备受重视。地壳运动中挤压矿物质，通常伴随着磁场辐射信号以波动方程规律向周围空间场辐射，该信号属于窄带宽频辐射信号，多表现为非线性、非平稳特征。宽带电磁辐射信号携带矿物种类、空间分布、地壳结构等重要信息，因此有效的检测、估计宽带电磁场辐射信号的波形、频率、相位等物理参数变化规律，是用于地质勘探矿物质的重要技术。实际矿物质中往往具有较复杂的成分，且辐射信号在传播过程中受到空间场背景辐射噪声影响，不可避免地致使接收到的目标信号出现能量衰减、形态扭曲、频率偏移、相位畸变等问题。可见，地壳中矿物质所产生的磁场辐射信号的信号形态及其参数估计，已成为矿物质检测领域面临的重要问题。

在地球深部探测领域常用的地球物理手段包括有重、磁、电、震、放射性等技术方法，在实际应用中又各有优势和不足之处。比如，电磁方法探测深度相对较大，野外仪器轻便，但分辨率和抗干扰能力较差；地震勘探分辨高、深度大，但其对地形条件要求较高，且费用昂贵。此外，传统方法以地下介质的物理差异为基础，属于间接找矿方法，大部分无法实现深部矿体的直接探测定位；发射性勘探虽属直接探测方法，但其勘探深度不足。

电磁声子探测以矿物晶体在地应力作用下的自发电磁辐射为物理基础，是一种被动源地球物理探测技术，具有分辨率高、探测深度大等特点，可实现深部地质结构的精细探测；仪器轻便，可搭载于汽车及无人机实现复杂地形的高效探测。

目前，窄带宽频辐射信号的常用检测包括以下几种方法：

(1)谱估计法

谱估计是以信号周期性和相关性为理论基础，根据谱峰位置确定信号频率。其基本思想是从观测数据估计自相关函数后，将自相关估计值乘以窗函数，最后对加窗延迟估计值做傅立叶变换得到功率谱估计。然而，谱估计法窗函数的截断与计算过程的离散化不可避免地出现能量泄露和栅栏效应，使得其频率、幅度和相位的估计都存在较大的误差，尤其在低信噪比情况下，估计精度会进一步恶化

(2)线性预测法

线性预测法是利用采样点相位的线性模型来提取信号频率的方法。其基本思想是对于含噪的辐射信号模型x(n)，可定义预测误差函数e(n)为

其中ω₀为估计频率。若代价函数

最小，即可得频率函数的估计值

但是，利用信号采样点线性预测模型估计频率的缺点是即使在信噪比较高时，性能也不能达到克拉美罗界，即检测准确性不高。

(3)最大熵谱估计法

最大熵谱估计法是对观测区间以外的点的自相关函数按最大熵进行外推。其基本思想是：已知p+1个相关函数r(m)(m＝0,1,…,p)，希望利用这p+1个值对m＞p时的自相关函数进行外推，其外推准则是外推后所对应的时间序列应具有最大的熵，因此要使熵率最大，即使得熵谱最大。该方法对于较短数据的信号序列，结果可以获取较高分辨率，但是较长数据会导致最大熵模型的计算量巨大，甚至无法得到结果。

(4)时频分析方法

传统的时频分析方法虽然同时兼顾了时域和频域两方面，但也普遍存在一些致命缺陷，影响实际使用。如以基于核函数分解的线性方法，由于固定的窗口或基函数自适应性不佳，不适合分析多尺度和突变信号，而且有限长度的基函数会不可避免地产生能量泄漏，影响分辨力；以基于Wigner分布的二次时频分布方法存在严重交叉干扰项，为抑制二次时频分布所固有的交叉项干扰，通常需要在分辨率和交叉项干扰之间取得某种折衷，因此不可避免地破坏了信号自项成分的集中，导致其分辨率降低。

地球内部岩石层能对地球电磁场产生扰动，这种扰动能够通过电磁场的变化到达地球表面，进而在地球上形成了复杂多变的电磁场，这是我们研究地球内部性质和结构的地球物理依据。以矿物晶体在地应力作用下的自发电磁辐射为物理基础，是一种被动源地球物理探测技术，具有分辨率高、探测深度大等特点，可实现深部地质结构的精细探测；仪器轻便，可搭载于汽车及无人机实现复杂地形的高效探测。

发明内容

磁法勘探是地球物理勘探的一个重要分支，本专利提出的矿物质晶体的辐射检测方法，能够利用地下岩矿石的磁性差异，达到探测地下地质结构，寻找矿产资源和解决水文、工程地质问题和环境监测等重要问题。本专利考虑了已有信号处理方法在上述几方面的不足，提出了一种基于自动移位法的经验模式分解和多维特征提取IMF方法实现瞬态宽带电磁波的探测，以期更好地揭示矿物质辐射信号本质，提高勘测矿物质的水平。

本专利具有获取宽带信号稳定特征的检测能力：根据宽带电磁信号的不同形态，自动选取恰当的经验模式分解方法，抑制了信号两端偏离信号真实的包络线而产生发散现象；由于宽带电磁信号具有多种特征，利用时域、相位和能量三种特征作为选取可以代表原信号特征的IMF的依据，进而实现了宽带电磁辐射信号随机成分和确定性成分的分离，该随机信号即为电磁辐射信号的特征，可在矿物质勘探等领域具有广泛的用途。

本发明主要为地壳中受到地应力作用下矿物的宽带电磁信号检测的方法设计：

一种基于地质探测的宽带电磁信号(20Hz-150kHz)处理技术，包括：在原始信号添加白噪声，归一化处理，得到新信号；根据带噪信号的包络特征，应用自动移位法、波形匹配法或比例延拓法做经验模式分解，得到IMF分量和残余分量，其中自动移位法首先移出变化最快的IMF分量，并且随着IMF级别的增加，可使宽带电磁信号的IMF分量快速平滑；消除IMF分量中的白噪声，获得真实的IMF分量；计算各IMF的时域一阶差分能量、相位一阶差分能量和归一化能量占比，选取三种特征均最大的IMF，即可作为辐射信号的主要特征；将选取出的IMF频率与特定矿物质所辐射的电磁波频率对比，若一致则说明矿石种含有该矿物质元素。

快速迭代经验模式分解和多维特征提取IMF方法首先利用了白噪声频谱均匀分布的特性，对信号添加白噪声，进而可以将不同时间尺度的信号分布到合适的参考尺度上；其次，依次应用自动移位法、波形匹配法或比例延拓法快速对宽带电磁信号进行分解；然后利用白噪声的零均值统计特性，采用足够的测试量并取全体均值，消除噪声，从而避免了模态混叠；最后利用时域、相位和能量三种特征可以获取有效的IMF分量。具体步骤如下：

第一步：在原始辐射信号s(t)多次添加等长度等幅度的高斯分布白噪声n(t)，并进行归一化处理，得到新的信号x(t)。通过高斯白噪声的引入，使极值点分布更趋均匀，避免了传统经验模式分解方法由于间歇性髙频分量造成的模式混叠。

第二步：情况1：对带噪信号x(t)应用自动移位法实现信号的经验模式分解，步骤如下：

1)设定h(t)＝x(t)，h_old(t)＝h(t)；

2)计算h_old(t)的各极大值与极小值；

3)使用三次样条函数，在各极大值/极小值之间做插值，进而得到上包络e_max(t)/下包络e_min(t)；

4)计算上包络和下包络的均值，如

5)如果h_new(t)＝h_old(t)-m(t)满足IMF的两个条件(其一为局部极大值以及局部极小值的数目之和必须与零交越点的数目相等或是最多只能差1，其二为在任何时间点，局部最大值所定义的上包络线与局部极小值所定义的下包络线，取平均要接近为零)，则h_new(t)为IMF，即c_n＝h_new(t)，否则设定h_old(t)＝h_new(t)，重复2)到5)过程，最后会得到c₁＝h_new(t)；

6)若得到第n个IMFc_n，设定h_old(t)＝h_old(t)-c_n，重复2)到5)过程，最后会得到所有IMFc_n+1。

自动移位法首先移出变化最快的IMF分量，并且随着IMF级别的增加，可使宽带电磁信号的IMF分量快速平滑，这种快速平滑特点可以适用于绝大部分的宽带电磁波信号的IMF分解。

情况2：对于受到噪声干扰的少类电磁波信号情况，当自动移位法分解过程可能会导致端点效应和/或模式混叠现象，此时采用波形匹配法使用信号内部波形与端点处变化趋势最为相似的子波来对端点处的信号进行延拓。波形匹配法使用信号内部波形和端点处变化趋势最为相似的子波来对端点处的信号进行延拓。该方法充分考虑了信号边缘处的变化趋势，使端点处的延拓更加合理。以原始信号x(t)向右延拓为例，波形匹配的过程：

(1)求信号x(t)的极值点，极大值集合记为P，极小值集合记为Q；

(2)原始信号的终点为m₀，向左取一段曲线段使之至少包含一个局部极大值点、一个局部极小值点和一个过零点，求s(t)最左端的两个极值点极大值P₁、极小值Q₁，从起始点到极小值Q₁的波形记为x'₁(t)，且其长度记为l；

(3)以极大值集合{P_i}中的每一个极大值P_i作为参考点，计算该段相同长度的波形x'_i(t)和x'₁(t)的匹配度，记为m[x'₁(t),x_i'(t),P_i]，匹配度公式为：

N为序列长度；

(4)使用

确定延拓波形

若/>

(其中α为一常数)，则取波形/>

左侧包含了一个极大值和极小值的子波，作为原始x(t)左端的延拓，延拓完毕；否则转(4)；

(5)若没有找到匹配的波形则直接指定端点处的极大和极小值：取原始信号最右端的两个相邻极大值点的均值作为右端点的极大值，取信号最右端的两个相邻极小值点的均值作为右端点的极小值。

情况3：如果信号内部没有和端点处波形相匹配的波形，那么波形匹配方法将不再有效，不能很好的体现波形匹配的优点。当信号内部没有和端点处波形相匹配的波形时，本专利结合使用比例延拓法对端点进行处理，从而更有效的解决端点效应问题。

本专利针对不同情况，分别使用自动移位法、波形匹配法或者比例延拓法寻找IMF，从而更有效的解决端点效应和模式混叠问题。这样，我们可以得到一组IMF分量{c_i(t)|1≤i≤N}和一个残余分量r(t)。

第三步：为了获得真实的来自目标电磁波的IMF分量，需要消除前面所加入的白噪声，假设加入了N次，则真实分解结果可以表示为:

所加入的白噪声遵循如下规律：

式(3)中：e为标准差，即输入信号与相应IMF分量重构结果的偏离；k为白噪声的幅值。可以看出N值取值与幅值k有关，但白噪声幅值对于分解结果影响不大，因此总体来说集合经验模式分解分析与人的主观介入不存在必然联系，仍然具有自适应性。

第四步：计算各IMF的时域一阶差分、相位一阶差分和归一化能量分布三种特征，选取各IMF的时域变化的强度、相位变化的强度以及各IMF能量的权重均较大的一个IMF，根据此IMF的频率值确定信号中是否含有目标电磁信号。所述IMF的时域一阶差分、相位一阶差分和归一化能量分布三种特征分别如下：

1)IMF的时域变化的强度D_t,i：

2)IMF的相位变化的强度D_p,i：

3)IMF能量的权重E_norm,i：

考量此三种特征的均为最大值的IMF，即作为辐射信号的主要特征。不同的矿物质在挤压过程中会释放不同频率的电磁波，通过对比选定的IMF与某矿物质电磁波频率是否一致，可判定该矿石中是否含有此矿物质。

本发明的有益效果为：

在地壳中挤压矿物质晶体过程会产生宽带电磁辐射信号，(1)本专利对该信号在经验模式分解过程中的端点效应和特征提取问题，提出了分解IMF的新依据。本专利的分解IMF方法可以有效抑制信号两端偏离信号真实的包络线而产生发散现象。(2)本专利根据时域差分、相位差分以及能量权重三种特征，提供了选取目标IMF的新依据，这实现了电磁辐射信号随机成分和确定性成分的分离，该随机信号即为电磁辐射信号的特征，进而可在了解矿物质晶体的辐射特点领域应用前景广泛。

附图说明

图1为本专利所提供的基于地质探测的瞬态宽带电磁探测信号处理方法流程示意图。

图2为轴信号及其自动移位法实现信号的经验模式分解结果。图中：(图2a)IMF1；(图2b)IMF2；(图2c)IMF3；(图2d)IMF4；(图2e)IMF5；(图2f)IMF6；(图2g)IMF7；(图2h)IMF8；(图2i)IMF9；(图2j)IMF10；(图2k)IMF11；(图2l)IMF12；(图2m)IMF13；(图2n)原信号。

图3为各层IMF的频谱。图中：(图3a)IMF1；(图3b)IMF2；(图3c)IMF3；(图3d)IMF4；(图3e)IMF5；(图3f)IMF6；(图3g)IMF7；(图3h)IMF8；(图3ia)IMF9全局；(图3ib)IMF9局部；(图3ja)IMF10全局；(图3jb)IMF10局部；(图3ka)IMF11全局；(图3kb)IMF11局部；(图3la)IMF12全局；(图3lb)IMF12局部；(图3ma)IMF13全局；(图3mb)IMF13局部。

图4为IMF8-IMF13频谱叠加图；图4a频谱叠加全局图；图4b频谱叠加局部图。

具体实施方式

挤压矿石的辐射信号由X、Y、Z三个维度方向的传感器拾取，采样频率500K Hz。其中Y轴方向的信号能量较为明显，故对其做分解，过程如下

第一步：在原始辐射信号s(t)添加了3次高斯分布白噪声n(t)，并进行归一化处理，得到新的信号x(t)；

第二步：对带噪信号x(t)应用自动移位法实现信号的经验模式分解，步骤如下：

1)设定h(t)＝x(t)，h_old(t)＝h(t)；

2)计算h_old(t)的各极大值与极小值；

3)使用三次样条函数，在各极大值/极小值之间做插值，进而得到上包络e_max(t)/下包络e_min(t)；

4)计算上包络和下包络的均值，如

5)如果h_new(t)＝h_old(t)-m(t)满足IMF的两个条件，则h_new(t)为IMF，即c₁＝h_new(t)，否则设定h_old(t)＝h_new(t)，重复2)到5)过程，最后会得到c₁＝h_new(t)；

6)若得到第n个IMFc_n，设定h_old(t)＝h_old(t)-c_n，重复2)到5)过程，最后得到13个IMF分量{c_i(t)|1≤i≤N}和1个残余分量r(t)。所得到的13个IMF分量避免了模态函数的畸变，尽可能保证原信号的固有特征，因而可有效抑制端点效应和模式混叠问题，

第三步：为了获得真实的来自目标电磁波的IMF分量，消除前面所加入的白噪声；

第四步：计算13个IMF的时域一阶差分能量、相位一阶差分能量和归一化能量占比三种特征，第12个IMF三种特征的数值均是最大的，即可作为辐射信号的主要特征，且与挤压石英所辐射的电磁频率400Hz一致。这种选取函数的方法可以将最能代表原信号特征的本征模态提取出来，作为区分其他信号特征的重要依据。

Y轴原信号波形及经验模式分解的分解结果如图2所示，为了进一步了解各IMF信息，前13个IMF的频谱如图3所示。

由图2、图3可知，IMF12与IMF13为近似正弦波，代表挤压矿物质的所辐射信号的基频成分，约400Hz，且在图4所示的IMF8-IMF13频谱叠加图可以看出，且IMF12的基频能量最为显著，为10的5次幂能级。IMF9-IMF11分别以2、4倍频成分为主。而IMF1-IMF8中杂波成分较多，且能量集中于5KHz以上的高频部分，属于低频部分的基频信号及其谐波被淹没。

Claims (5)

1.一种基于地质探测的瞬态宽带电磁探测信号处理方法，其特征在于，包括以下步骤：

S1)在原始宽带电磁辐射信号s(t)添加高斯分布白噪声n(t)，进行归一化处理，得到带噪信号x(t)；

S2)按照步骤S2.1)～S2.3)，对带噪信号x(t)应用自动移位法、波形匹配法或比例延拓法，得到一组IMF分量{c_i(t)|1≤i≤N}和一个残余分量r(t)：

S2.1)对带噪信号x(t)应用自动移位法实现信号的经验模式分解，得到IMF分量{c_i(t)|1≤i≤N}和残余分量r(t)；

S2.2)判断S2.1)分解过程是否导致端点效应和/或模式混叠现象：

当分解过程不导致端点效应和/或模式混叠现象时，S2)结束；

当分解过程导致端点效应和/或模式混叠现象时，应用波形匹配法对带噪信号x(t)的端点处的信号进行延拓；

当波形匹配法在信号内部没有找到与端点处波形相匹配的波形时，使用比例延拓法对带噪信号x(t)的端点处的信号进行延拓；

S2.3)应用经验模式分解，重新得到IMF分量{c_i(t)|1≤i≤N}和残余分量r(t)；

S3)应用公式(Ⅰ)消除IMF分量{c_i(t)|1≤i≤N}中的白噪声，获得真实的IMF分量：

式(Ⅰ)中，N为所加入的白噪声序列的个数，c_i(t)为最终得到的第i个分量IMF，c_ij(t)为添加第j个白噪声处理后的第i个分量IMF，1≤j≤N，i＝1,2,3,……正整数；

所加入的白噪声遵循以下公式：

式(Ⅱ)中，e为标准差，即输入信号与相应IMF分量重构结果的偏离；k为白噪声的幅值；

S4)按照步骤S4.1)～S4.3)，分别计算真实的IMF分量的三种特征，包括时域一阶差分、相位一阶差分及归一化能量分布，得到各IMF分量的时域变化的强度、相位变化的强度以及能量的权重，判断信号中是否含有目标电磁信号；第i个IMF分量c_i的时域变化的强度、相位变化的强度以及能量的权重的计算方式为：

S4.1)第i个IMF分量c_i的时域变化的强度D_t,i：

S4.2)第i个IMF分量c_i的相位变化的强度D_p,i：

式(Ⅳ)中，

是第i个IMF分量c_i的相位；

S4.3)第i个IMF分量c_i的能量的权重E_norm,i：

式(Ⅲ)-(Ⅴ)中，

是IMF的相位；

S5)根据S4)判断结果，选取真实的IMF分量中的特定分量重构信号，得到目标电磁信号。

2.如权利要求1所述的方法，其特征在于，应用波形匹配法处理信号，包括如下步骤：

S2.2.1)求带噪信号x(t)的极值点，极大值集合记为P，极小值集合记为Q；

S2.2.2)以带噪信号x(t)进行延拓的一端的终点为m₀，以m₀为端点在带噪信号x(t)上取一段曲线段x₁(t)，使之至少包含一个局部极大值点、一个局部极小值点和一个过零点，求x₁(t)上的两个极值点极大值P₁、极小值Q₁，m₀到Q₁的距离记为l₁；

S2.2.3)以极大值集合P中的每一个极大值P_i作为参考点将带噪信号x(t)分为若干个波形曲线片段x_i'(t)，选取每一个波形曲线计算该线段的波形x_i'(t)与x₁(t)的匹配度，记为m[x'₁(t),x_i'(t),P_i]，匹配度公式为：

式(Ⅵ)中，N为序列长度；

S2.2.4)令

其中/>

是式(Ⅵ)的波形曲线片段，/>

是式(Ⅵ)的极大值的解：

若

则紧邻/>

在/>

与m₀间选取包含了一个极大值和一个极小值的子波作为原始宽带电磁辐射信号s(t)在m₀处的延拓，延拓完毕，其中α为常数；

若没有找到匹配的波形，则采用比例延拓法指定端点处的极大和极小值。

3.如权利要求2所述的方法，其特征在于，所述采用比例延拓法指定端点处的极大和极小值，是指以m₀临近的两个相邻极大值点的均值作为m₀端的极大值；以m₀临近的两个相邻极小值点的均值作为m₀端的极小值。

4.如权利要求1所述的方法，其特征在于，时域变化的强度D_t,i、相位变化的强度D_p,i、能量的权重E_norm,i均最大的IMF作为辐射信号的主要特征；通过对比IMF与某矿物质电磁波频率是否一致，判定矿石中是否含有对应矿物质。

5.如权利要求4所述的方法，其特征在于，包括如下步骤：

S4.2.1)对每个IMF进行傅里叶变换，得到：

C_i(t)＝F{c_i(t)}

(Ⅶ)；

式(Ⅶ)中，C_i(t)＝{C_i,1,C_i,2,...C_i,k...,C_i,M}是第i个IMF在频域空间的一组复系数；第i个IMF的第k个系数C_i,k满足式(Ⅷ)：

式(Ⅷ)中，k代表频率

处的复系数角标，M是IMF频谱的复系数个数；

S4.2.2)计算相位谱，将系数的虚部和实部的比值运用于反正切函数，如下式所示：

式(Ⅸ)中，I[C_i(t)]与R[C_i(t)]分别表示i个IMF的实部与虚部；

S4.2.3)计算一阶差分的相位变化的强度D_p,i：

Images