一种基于信噪分类的音频大地电磁信号去噪方法

摘要

本发明公开了一种基于信噪分类的音频大地电磁信号去噪方法，包括：提取音频大地电磁信号样本；分别计算每个音频大地电磁信号样本多重分形谱的不均匀度和不规则度；利用每个音频大地电磁信号样本多重分形谱的不均匀度、不规则度以及每个音频大地电磁信号样本的类别值来训练预设分类模型得到信噪分类数学模型；依据信噪分类数学模型对待处理的实测音频大地电磁信号进行分类得到非强干扰的信号段和受强干扰的信号段；对受强干扰的电磁信号段进行匹配追踪去噪处理；将去噪处理后的信号段与非强干扰的信号段进行合并得到重构的音频大地电磁有用信号。通过上述方法获取质量更高的音频大地电磁信号，能有效避免非强干扰信号被滤除的情况，提高去噪精度。

说明书

技术领域

本发明属于电磁勘探数据处理领域，具体涉及一种基于信噪分类的音频大地电磁信号去噪方法。

背景技术

目前，音频大地电磁测深(Audio magnetotelluric，AMT)是基于大地电磁测深原理，是一种平面波卡尼亚电阻率频率域电磁测深法。由于人文电磁干扰源复杂，导致有些干扰在音频大地电磁时间序列上表现出明显特征，如脉冲干扰、工频噪声等；有些则不明显，但在频谱上有较强的特征，如地磁噪声等；有些则在时间序列和频谱上均显示不出任何特征，如地质噪声等。正是因为这些噪声与电磁信号叠加在一起，严重影响了音频大地电磁数据质量，制约了地电结构的可解释性。

随着国民经济和重工业生产的飞速发展，音频大地电磁测深适用于深度的工程勘探，其中，音频大地电磁信号采集时间短，受强干扰影响情况急剧加深。因此，如何挖掘音频大地电磁信号中更鲁棒性的特征参数，将音频大地电磁有用信号和强干扰进行定量区分，精细地保留更丰富的低频段音频大地电磁有用数据将是其研究热点。目前，地球物理学家们相继提出了一系列技术，主要集中在四个方面：远参考法、稳健阻抗估计、时频转换技术和时间域去噪。这些方法在一定程度上均可改善音频大地电磁数据质量。然而，由于这些技术仅是对采集的数据进行整体处理，缺少对信号和干扰的分类，且无法划分非强干扰的信号和强干扰，导致非强干扰的信号也被随之滤除，在矿集区要获取高置信度的音频大地电磁有用信号面临巨大困难。

发明内容

本发明的目的是提供一种基于信噪分类的音频大地电磁去噪方法，通过对音频大地电磁信号进行先信噪分类后去噪的方式来获取质量更高的音频大地电磁信号，解决现有的将非强干扰的信号也进行滤除的问题。

一种基于信噪分类的音频大地电磁去噪方法，包括如下步骤：

步骤1：从采集的音频大地电磁信号中提取电磁信号样本；

其中，提取的所述音频大地电磁信号样本分为非强干扰的电磁信号和受强干扰的电磁信号；

步骤2：分别计算每个电磁信号样本的多重分形谱，并计算出各自多重分形谱的不均匀度和不规则度；

步骤3：利用每个电磁信号样本多重分形谱的不均匀度、不规则度以及每个电磁信号样本的类别值来训练预设分类模型得到信噪分类数学模型；

其中，每个电磁信号样本的类别值表示电磁信号是否受强干扰来编码，所述信噪分类数学模型的输入数据为电磁信号多重分形谱的不均匀度和不规则度，输出数据为对应电磁信号的类别值；

步骤4：依据步骤3中的信噪分类数学模型对待处理的实测音频大地电磁信号进行分类得到非强干扰的电磁信号段和受强干扰的电磁信号段；

其中，将实测的音频大地电磁信号划分为J段，并分别计算每段电磁信号多重分形谱的不均匀度和不规则度，再将每段电磁信号的不均匀度和不规则度输入信噪分类数学模型得到每段电磁信号的类别值，J为正整数；

步骤5：对步骤4中划分为受强干扰的音频大地电磁信号段进行匹配追踪去噪处理；

步骤6：将去噪处理后的音频大地电磁信号段与步骤4中划分为非强干扰的音频大地电磁信号段进行合并得到重构的音频大地电磁有用信号。

本发明针对音频大地电磁信号选用的不均匀度和不规则度两个特征进行分类处理，具体是将非强干扰和受强干扰的音频大地电磁信号多重分形谱的不均匀度和不规则度作为特征参数进行训练，得到信噪分类数学模型，再利用信噪分类数学模型对实测音频大地电磁信号进行分类，划分受强干扰段和非强干扰段的电磁信号段，随后将划分为受强干扰的电磁信号段进行匹配追踪去噪处理，并与划分为非强干扰的音频大地电磁信号段进行合并得到重构的音频大地电磁信号，进而有效避免了现有技术对音频大地电磁数据进行整体处理的同时，一并滤除非强干扰信号的弊端。

且不均匀度和不规则度不仅能从整体上对音频大地电磁信号的本质特征进行全面描述，而且能对强干扰信号段的局部特征进行更为细致的刻画。综合对比现有衡量复杂度特征的样本熵和模糊熵等进行信号分类的方案，本发明选用不均匀度和不规则度是从音频大地电磁信噪的角度上提升了区分精度，同样达到了精确分类的效果。

进一步优选，多重分形谱为分形维数f(α)与标度指数α的关系函数，每个多重分形谱的不均匀度和不规则度的计算公式如下：

Δα＝α_max-α_min

Δf(α)＝f(α)_max-f(α)_min

式中，Δα为多重分形谱的不均匀度，α_max、α_min分别为多重分形谱中标度指数α的最大值、最小值，Δf(α)为多重分形谱的不规则度，f(α)_max、f(α)_min分别为多重分形谱中分形维数f(α)的最大值、最小值。

进一步优选，电磁信号多重分形谱的计算过程如下：

首先，基于音频大地电磁信号的时间序列沿着时间轴划分为N个尺度为r的盒子，并获取每个盒子的概率测度；

其次，基于每个盒子概率测度采用如下公式计算出配分函数；

式中，χ_q(r)为配分函数，P_i(r)为第i个盒子的概率测度，r为盒子的尺度；

然后，基于配分函数按照如下公式计算质量指数；

式中，τ(q)为质量指数；

最后，根据勒让德变换依次计算出标度指数α和分形维数f(α)得到电磁信号的多重分形谱；

α＝dτ(q)/dq

f(α)＝qα-τ(q)。

进一步优选，步骤5中采用匹配追踪算法对受强干扰的电磁信号段进行去噪处理，具体过程如下：

步骤5.1：利用sin/dct原子和sym/db小波原子构造冗余字典；

式中，D为冗余字典，分别为冗余字典D中的原子，n₁为冗余字典D中原子的数量，所述冗余字典将自适应地选取原子进行匹配；其中，冗余字典D中的原子类型包含了sin/dct原子和sym/db小波原子；

步骤5.2：从所述冗余字典中选取与受强干扰的电磁信号段最优匹配的原子；

其中，原子最优匹配规则如下：

式中，f为受强干扰的电磁信号段，g_r为最优匹配的原子，g_k为冗余字典中第k个原子，<>

步骤5.3：基于最优匹配的原子对受强干扰的电磁信号段按照如下公式进行分解；

其中，电磁信号段分解公式如下：

f＝<f,g_r>g_r+R¹f；

式中，<f,g_r>g_r为受强干扰的电磁信号段f在最优匹配的原子gr上的分量，R1f为最优匹配的原子g_r对受强干扰的电磁信号段f进行匹配后的残差；

步骤5.4：基于最优匹配的原子对步骤5.3得到的残差进行分解，并对分解后得到的残差重复迭代分解直至残差小于预设精度后得到去噪后的电磁信号段；

其中，去噪后的电磁信号段如下所示：

式中，R^kf表示得到的第k个残差，n₂为得到残差的个数。

本发明构建的全新冗余字典能自适应匹配受强干扰段的音频大地电磁信号。其中，通过多次实验论证，本发明构建的冗余字典可以更有针对性地对受强干扰的音频大地电磁信号段进行去噪处理，尤其是匹配原子后进行线性组合时，冗余字典中的sin/dct原子适用于匹配受强干扰信号段中的平稳成分，sym/db小波原子适用于匹配受强干扰信号段中的突变成分。即针对受干扰信号段本发明构建的冗余字典将自适应地选取原子进行匹配，若受干扰信号段为突变成分，则自适应匹配到sym/db小波原子，若为平稳成分，则自适应匹配到sin/dct原子。

其中，sin/dct原子和sym/db小波原子构造冗余字典来匹配音频大地电磁信号中的强干扰，较好地描述了受强干扰信号段中的平稳及突变成分；对于划分为受强干扰的信号段，能有的放矢地剔除干扰，更为丰富地保留原始微弱的音频大地电磁有用信号。

其中所述预设精度为经验值，选取为趋于0的范围。

进一步优选，步骤1中提取的音频大地电磁信号样本中受强干扰的电磁信号为含类三角波干扰和含类脉冲干扰的电磁信号。

本发明将选择含类三角波干扰和含类脉冲干扰的电磁信号作为受强干扰的电磁信号，无干扰的电磁信号作为非强干扰的电磁信号。

进一步优选，步骤3中所述预设分类模型为支持向量机或BP神经网络模型。

有益效果

1、本发明基于非强干扰的电磁信号以及受强干扰的电磁信号样本的特征来训练得到信噪分类数学模型，进而利用信噪分类数学模型将实测音频大地电磁信号进行分类得到非强干扰和受强干扰的电磁信号段，再对受强干扰的电磁信号段进行匹配追踪去噪处理，最后将去噪处理后的电磁信号段与非强干扰的电磁信号段组合重构得到音频大地电磁信号。其中，通过对电磁信号先分类再去噪的方式可以有效避免现有技术中基于整体处理而将非强干扰的信号进行滤除的情况。

2、本发明将音频大地电磁信号多重分形谱的不均匀度和不规则度作为分类特征，其中，不均匀度特征反映了音频大地电磁非强干扰信号段与受强干扰信号段的区别，为了利用多重分形谱准确描述音频大地电磁非强干扰信号段和受强干扰信号段，结合不规则度特征，从多重分形谱图和特征参数值上进一步分析了音频大地电磁信噪区别，综合对比现有衡量复杂度特征的样本熵和模糊熵，从音频大地电磁信噪的角度上提升了区分精度，有利于后续支持向量机的信噪分类；多重分形谱的不均匀度和不规则度不仅从整体上对音频大地电磁信号的本质特征进行了全面描述，而且对强干扰信号段的局部特征进行了更为细致的刻画，进而最终提高了信噪分类数学模型的可靠性，提升了分类的准确度。

3、本发明提供了匹配追踪算法，利用sin/dct原子和sym/db小波原子构造冗余字典来匹配受强干扰的音频大地电磁信号，较好地描述了受强干扰的电磁信号段中的平稳及突变成分；对于划分为受强干扰的信号段，能有针对性地剔除干扰，最大限度地保留原始微弱的音频大地电磁有用信号。

4、本发明提供的预设分类模型可选为支持向量机，基于统计学习理论的支持向量机分类方法通过将能表征音频大地电磁信噪特征的多重分形谱作为支持向量机的输入来分类音频大地电磁信号中的非强干扰的电磁信号和受强干扰的电磁信号。

5、本发明将多重分形谱、基于统计学习的分类模型和匹配追踪相结合，充分利用它们在特征提取、模式识别、信噪分离等方面的优势。将实测音频大地电磁信号划分为非强干扰的有用信号段和强干扰信号段，从而有的放矢地从划分为强干扰的信号段中提取出微弱的音频大地电磁有用信号，相对于整体去噪方法保留了更多的有用信息不被过处理，为后续电磁法反演、解释提供了更为真实、可靠的音频大地电磁测深数据。

附图说明

图1为本发明实施例提供的流程图。

图2为本发明提供的样本库中电磁信号样本的示意图，其中(a)、(b)、(c)图分别为未受干扰的音频大地电磁信号、含类三角波干扰和含类脉冲干扰的音频大地电磁信号的示意图。

图3为本发明提供的样本库中三类信号的多重分形谱，其中(a)、(b)、(c)图分别为未受干扰的音频大地电磁信号、含类三角波干扰的音频大地电磁信号、含类脉冲干扰的音频大地电磁信号的多重分形谱的示意图。

图4为支持向量机对样本库的信噪分类效果图。

图5为匹配追踪算法中不同原子的去噪效果，其中(a)、(b)分别为sin/dct原子构造冗余字典的去噪效果、sym/db小波原子构造冗余字典的去噪效果。

图6为模拟信号信噪分类和去噪效果，其中(a)、(b)分别为模拟三角波干扰的信噪分类和去噪效果、模拟脉冲干扰的信噪分类和去噪效果。

图7为实测音频大地电磁信噪分类和去噪效果以及与匹配追踪整体处理对比效果示意图，其中(a)图为含类三角波干扰的实测音频大地电磁信噪分类和去噪效果以及与匹配追踪整体处理的对比效果图，图7(b)为含类脉冲干扰的实测音频大地电磁信噪分类和去噪效果以及与匹配追踪整体处理的对比效果图。

图8为音频大地电磁实测点1365C21B的卡尼亚电阻率-相位(R-P)曲线对比图，其中左图表示R_xy-P_xy分量(XY方向)，右图表示R_yx-P_yx分量(YX方向)；曲线1为实测点原始数据的卡尼亚电阻率-相位曲线，曲线2为匹配追踪整体处理的卡尼亚电阻率-相位曲线，曲线3为本发明处理的卡尼亚电阻率-相位曲线。

图9为音频大地电磁实测点1365C21B的极化方向对比图，其中(a)、(b)分别为电道4Hz数据的极化方向对比、磁道4Hz数据的极化方向对比。

具体实施方式

下面将结合实施例对本发明做进一步的说明。

如图1所示，本发明公开了一种基于信噪分类的音频大地电磁信号去噪方法，包括如下步骤：

步骤1：从采集的音频大地电磁信号中提取电磁信号样本并构建成样本库。

本实施例中，从采集的音频大地电磁时间序列中提取50个未受干扰的AMT信号、50个含类三角波干扰的AMT信号和50个含类脉冲干扰的AMT信号构造样本库，每个电磁信号样本的采样长度为240。再将此三类信号中150个样本进行整合生成一个样本库，即150个采样长度为240的信号段。本实施例中，将未受干扰的AMT信号作为非强干扰信号的样本。如图2所示，(a)、(b)、(c)分别为样本库中未受干扰的音频大地电磁信号、含类三角波干扰和含类脉冲干扰的音频大地电磁信号样本示意图。

步骤2：分别计算每个电磁信号样本的多重分形谱，并计算出每个电磁信号样本多重分形谱的不均匀度和不规则度。

其中，如图3所示为样本库中上述三类信号的多重分形谱，如(a)图为未受干扰的音频大地电磁信号的多重分形谱，(b)图为含类三角波干扰的音频大地电磁信号的多重分形谱，(c)图为含类脉冲干扰的音频大地电磁信号的多重分形谱。

其中，每个电磁信号的多重分形谱的计算过程如下：

2.1)将电磁信号样本中的音频大地电磁时间序列沿时间轴划分为N个尺度为r＝2的盒子，在空间上将P_i(r)定义为归一化盒子数的分布概率测度；当分布概率测度P_i(r)在盒子上均匀分布，则定义不同的标度指数α来表示分布概率测度P_i(r)：P_i(r)～rα，标度指数α为描述分布概率测度P_i(r)局部奇异性强度的参数；此外，定义具有不同标度指数α的分形子集叫做分形维数f(α)。

2.2)对分布概率测度P_i(r)采用q次方进行加权求和得到配分函数χ_q(r)，配分函数χ_q(r)与尺度r可表示为：

其中，q∈[-10,10]

式中，χ_q(r)为配分函数，P_i(r)为第i个盒子的概率测度，r为盒子的尺度；

2.3)计算曲线lnχ_q(r)和曲线lnr的斜率，即质量指数τ(q)：

2.4)根据勒让德变换，得到标度指数α和分形维数f(α)：

α＝dτ(q)/dq；

f(α)＝qα-τ(q)。

通过上述过程，将每个电磁信号样本依次进行处理得到多重分形谱，如图3所示，多重分形谱为分形维数f(α)与标度指数α的关系函数。

基于多重分形谱计算出的不均匀度与不规则度，计算公式如下：

Δα＝α_max-α_min

Δf(α)＝f(α)_max-f(α)_min

式中，Δα为多重分形谱的不均匀度，α_max、α_min分别为多重分形谱中标度指数α的最大值、最小值，Δf(α)为多重分形谱的不规则度，f(α)_max、f(α)_min分别为多重分形谱中分形维数f(α)的最大值、最小值。

其中，如下表1所示为三类样本信号中三个电磁信号样本多重分形谱的不均匀度和不规则度参数值，并对比已有特征参数如样本熵和模糊熵，进一步从多重分形谱的不均匀度和不规则度特征出发，提升了后续支持向量机的信噪分类精度。

表1

步骤3：利用每个电磁信号样本多重分形谱的不均匀度、不规则度以及每个电磁信号样本的类别值来训练支持向量机，得到信噪分类数学模型。

本实施例中选用的预设分类模型为支持向量机模型，并以此为例说明；其他可行的实施例中，预设的分类模型还可以是其他分类模型，如BP神经网络模型。

其中，利用步骤2计算出每个电磁信号样本多重分形谱的不均匀度和不规则度来构建训练样本集。本实施例中，训练样本集表示为：{(Δα_j,Δf(α)_j,yj),j＝1,2,...Q}，Δα_j,Δf(α)_j,yj分别为输入的第j个训练样本点的不均匀度、不规则度以及类别值。Δα＝[Δα₁,Δα₂,...,Δα_Q]和Δf(α)＝[Δf(α)₁,Δf(α)₂,...,Δf(α)_Q]分别是训练样本点的不均匀度、不规则度集合，Q是三类样本总数，本实施例中Q等于150，y_j∈{-1,1},j＝1,2,...,Q，y_j＝-1表示受强干扰信号段，y_j＝1表示非强干扰的信号段。

根据现有支持向量机分类，将训练样本集作为训练阶段的输入数据，计算样本库三类样本中每个样本数据离超平面最近的距离，以及与超平面之间的最大距离，得到最优分类面，并通过最优分类面划分样本库中的音频大地电磁信号，生成一个信噪分类数学模型。应当理解，支持向量机的训练过程为现有过程，因此不对其进行具体描述；其中，信噪分类数学模型包含了现有支持向量机中拉格朗日乘子和核函数，训练后的信噪分类数学模型的输入为实测音频大地电磁信号段的多重分形谱的不均匀度和不规则度，输出为信噪分类数学模型的类别值，其输出类别值可用于实测音频大地电磁数据进行信噪分类。

步骤4：依据步骤3中的信噪分类数学模型对待处理的实测音频大地电磁信号进行分类得到非强干扰的电磁信号段和受强干扰的电磁信号段。

将实测音频大地电磁信号等间距划分为J段，并按照2.1)～2.4)提取每段信号多重分形谱的不均匀度和不规则度，利用生成的信噪分类数学模型对实测音频大地电磁信号中的J段信号均进行分类得到为非强干扰的信号段或受强干扰信号段的分类结果。当信噪分类数学模型输出类别值为1时，表示此时通过信噪分类数学模型划分的信号段为非强干扰的信号段；反之，当信噪分类数学模型输出类别值为-1时，表示此时划分的信号段为受强干扰的信号段，即强干扰信号段。如图4所示为利用支持向量机的分类效果。

步骤5：对步骤4中受强干扰的电磁信号段进行去噪处理。

具体的，采用匹配追踪算法进行去噪处理。其过程如下：

步骤5.1：利用sin/dct原子和sym/db小波原子构造冗余字典；

式中，D为冗余字典，分别为冗余字典D中的原子，n₁为冗余字典D中原子的数量，所述冗余字典将自适应地选取原子进行匹配；其中，冗余字典D中的原子类型包含了sin/dct原子和sym/db小波原子。

其中，对于强干扰信号段f，都能自适应地从冗余字典中匹配最优原子进行线性组合。

步骤5.2：从所述冗余字典中选取与受强干扰的电磁信号段最优匹配的原子；

其中，原子最优匹配规则如下：

式中，f为受强干扰的电磁信号段，g_r为最优匹配的原子，g_k为冗余字典中第k个原子，<>

当强干扰信号段匹配原子进行线性组合时，冗余字典中的原子将自适应地选取；其中，sin/dct原子自适应匹配平稳成分，sym/db小波原子自适应匹配突变成分。

步骤5.3：基于最优匹配的原子对受强干扰的电磁信号段按照如下公式进行分解；

其中，电磁信号段分解公式如下：

f＝<f,g_r>g_r+R¹f；

式中，<f,g_r>g_r为受强干扰的电磁信号段f在最优匹配原子gr上的分量，R1f为最优匹配原子g_r对受强干扰的电磁信号段f进行匹配后的残差；

步骤5.4：基于最优匹配的原子对步骤5.3得到的残差进行分解，并对分解后得到的残差重复迭代分解直至残差趋于0，即得到去噪后的电磁信号。其中，强干扰信号段f可由最优原子的线性组合再加上残差表示：

其中，由于残差呈指数衰减，去噪后的电磁信号段可近似表示如下：

其中，以图5中的(a)图和(b)图分别表示sin/dct原子适用于匹配信号中的平稳成分，sym/db小波原子适用于匹配信号中的突变成分。

步骤6：将去噪处理后的音频大地电磁信号段与步骤4中非强干扰的音频大地电磁信号段进行合并得到重构的音频大地电磁有用信号。

其中，如图6所示，(a)、(b)图分别表示模拟三角波干扰、模拟脉冲干扰的电磁信号分类处理过程。如图7所示，(a)、(b)图分别表示实测类三角波干扰、实测类脉冲干扰的电磁信号分类处理后的效果与现有的匹配追踪整体处理效果的对比。

本发明针对音频大地电磁数据易受强干扰影响，将多重分形谱的不均匀度和不规则度特征作为支持向量机的输入，通过支持向量机训练生成一个信噪分类数学模型来定量区分非强干扰的信号和强干扰，在此基础上利用匹配追踪算法有针对性地从辨识为强干扰的数据段中提取出微弱的音频大地电磁有用信号。本发明通过先对音频大地电磁时间序列进行信噪分类，后仅对划分为强干扰的信号段有的放矢地进行去噪，避免了整体方法处理时损失的有用信号。因此，本发明获取了更为真实、可靠的音频大地电磁数据。

通过对比原始数据、整体去噪方法和本发明方法处理后的卡尼亚电阻率-相位曲线和电磁场极化方向来对本发明的效果进行评价。如图8和图9所示，图8为音频大地电磁实测点1365C21B的卡尼亚电阻率-相位(R-P)曲线对比图，其中左图表示R_xy-P_xy分量(XY方向)，右图表示R_yx-P_yx分量(YX方向)；曲线1为实测点原始数据的卡尼亚电阻率-相位曲线，曲线2为匹配追踪整体处理的卡尼亚电阻率-相位曲线，曲线3为本发明处理的卡尼亚电阻率-相位曲线。图9为音频大地电磁实测点1365C21B的极化方向对比图，其中(a)图为电道4Hz数据的极化方向对比，(b)图为磁道4Hz数据极化方向对比。由于原始音频大地电磁实测点在时间域中受强干扰的影响，导致了R_xy分量呈45°上升的近源效应，相应的相位也趋近于0°，其极化方向在40°～80°之间表现出明显的聚集性和规律性；经整体去噪方法处理，强干扰虽能压制但同时损失了有用信号，导致卡尼亚电阻率-相位曲线在低频段跳变混乱；不难发现，经本发明方法处理后，实测点的卡尼亚电阻率-相位曲线更为光滑、连续，尤其在低频段，电阻率值处于合理的范围，近源效应得到了有效压制；同时对比观测电磁场极化方向可知，本发明方法处理后的极化点随机趋势增大，极化方向的分布逐渐呈离散、无序状态。为此，本发明方法表明：处理后实测点的特征更接近于天然电磁场的特征规律。

以上结合附图和具体实施方式，对本发明的技术领域、背景、目的、方案和有益效果做了进一步的详细说明，所应理解的是，本实施方式仅为本发明的优选方式而已，并不用于限制本发明，凡是本领域普通技术人员所具备的知识范围内，做出的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。