CN113640891A - 一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法 - Google Patents

Abstract

本发明属于地球物理勘探中瞬变电磁信号处理技术领域，公开了一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法，该方法包括：利用瞬变电磁探测系统获取一组实测含噪的瞬变电磁数据，对原始数据进行叠加；将含噪的瞬变电磁数据重新嵌入成轨迹矩阵Y，并利用奇异值分解法进行分解得到特征值向量λ和奇异值向量σ；提取奇异值向量σ中明显大于零的奇异值作为有效奇异值，根据有效奇异值重构信号并估计期望信号，将不同的估计信号对应到相应的奇异值中。解决地质勘探过程中采集的瞬变电磁数据被噪声干扰的问题，该方法实用性强、消噪后信噪比提升高，能够清晰地反映低阻地层分布；同时能够提取噪声，为后续的噪声分析做准备。

Description

一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法

技术领域

本发明属于地球物理勘探中瞬变电磁信号处理技术领域，公开了一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法。

背景技术

瞬变电磁探测系统是一种广泛应用于地球物理勘探领域的方法。利用其进行探测时，接收的电磁信号容易受到多种噪声影响，导致电磁数据中含有较为复杂的噪声干扰，尤其在高压线或者城区附近时，由于工频干扰很大，加之观测的是能量较弱的二次场，数据信噪比降低，影响数据解释结果的准确性。

CN110850482A公开了一种基于变分模态分解原理的瞬变电磁信噪分离方法，用于消除采集到的瞬变电磁信号中的噪声。该方法通过VMD参数选取，分解后得到一系列本征模态分量，实现瞬变电磁信号的消噪。但方法需预先定义模态数K不断寻找最优解，容易存在模态混叠的问题。

CN109885906A公开了一种基于粒子群优化的磁共振测深信号稀疏消噪方法，主要面向处理磁共振信号中的工频谐波干扰和随机白噪声。该方法重构信号库，并实现随机噪声的去除。但消噪过程中参数的选择对消噪结果影响较大，很难调整到最适宜的参数，并且容易陷入局部最优。

CN107045149A公开了一种基于双奇异值分解的全波核磁共振信号噪声滤除方法，针对单通道采集的全波MRS数据,可以同时去除工频谐波干扰和随机噪声的影响。但消噪过程中双奇异值分解相对繁琐，同时第二次奇异值分解时容易受到第一次分解结果的影响。

Li.Y("The baseline wander correction based on the improved ensembleempirical mode decomposition(EEMD)algorithm for grounded electrical sourceairborne transient electromagnetic signals,"Geoscientific Instrumentation,Methods and Data Systems,vol.9,(2),443-450,2020.)使用的基于EEMD的自适应对消滤波法可以抑制瞬变电磁信号中的基线漂移，但需要采集数据附近的噪声数据作自适应对消，噪声数据与原始数据相关性过大时消噪效果不佳。

发明内容

本发明所要解决的技术问题在于提供一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法，解决地质勘探过程中采集的瞬变电磁数据被噪声干扰的问题，该方法实用性强、消噪后信噪比提升高，能够清晰地反映低阻地层分布；同时能够提取噪声，为后续的噪声分析做准备。

本发明是这样实现的，

一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法，该方法包括：

利用瞬变电磁探测系统获取一组实测含噪的瞬变电磁数据，对原始数据进行叠加；

将含噪的瞬变电磁数据重新嵌入成轨迹矩阵Y，并利用奇异值分解法进行分解得到特征值向量λ和奇异值向量σ；

提取奇异值向量σ中明显大于零的奇异值作为有效奇异值，根据有效奇异值重构信号并估计期望信号，将不同的估计信号对应到相应的奇异值中；

其中根据噪声滤除的要求，采用粒子群优化算法求取奇异谱分解法中的两个重要参数：窗口长度M与重构信号时奇异值的阶数。

进一步，将含噪瞬变电磁探测数据重新嵌入，并利用奇异值分解法进行分解，具体包括：

将一维瞬变电磁数据y＝[y(1)，y(2)，...，y(N)]映射形成M×L轨迹矩阵Y，并重新排列；其中M为窗口长度，L＝N-M+1；

计算轨迹矩阵Y的自协方差矩阵C＝YY^T，并进行奇异值分解C＝U∑V^T，U为左特征向量矩阵、V为右特征向量矩阵，并得到其特征值向量λ，特征值向量表示如下：

λ＝[λ₁，λ₂，…，λ_M]，|λ₁|≥|λ₂|≥…≥|λ_M| (1)

同时，自协方差矩阵C的奇异值向量σ表示成如下形式：

进一步，利用粒子群优化算法获取窗口长度M，具体包括：

设有m个粒子在一个S维的搜索空间中组成群体，第i个粒子表示为一个S维的向量X_i，将X_i代入目标函数计算出对应的适应值；

令第i个粒子当前迭代的最优个体极值为P₁＝(p_i1，p_i2，...，p_iS)，速度为V_i＝(v_i1，v_i2，...，v_iS)，粒子群当前迭代的最优全局极值为P_g＝(p_g1，p_g2，...，p_gS)，利用公式(3)、公式(4)更新S维上第i个粒子的速度和位置：

V_i(t+1)＝ωV_i(t)+c₁v₁(P_i(t)-X_i(t))+c₂v₂(P_g(t)-X_i(t)) (3)

X_i(t+1)＝X_i(t)+V_i(t+1) (4)

其中t为迭代次数；ω为惯性权重，取ω∈[0.4，1.0]；c₁和c₂为学习因子，取c₁，c₂∈[0，2]；r₁和r₂为相互独立且服从[0,1]上均匀分布的伪随机数。

进一步，获取重构阶数的过程采用二分法及均值截断法重构信号，在对瞬变电磁数据重构时重构阶数选择为2。

进一步，重构信号包含以下步骤：

利用分解后的特征值对应的特征向量表示出轨迹矩阵Y的第j个分量轨迹矩阵：

将轨迹矩阵分量Y_j分成p个不相交的子集合J₁，J₂，…J_p，则合成矩阵Y_j表示如下：

计算各分量对应的合成矩阵，表示出原始矩阵Y：

重新分组后进行对角平均化，轨迹矩阵被转换成时间道为N的一维数据，令

为L×K的轨迹矩阵，根据对角平均化公式将

转换为一维时间序列；

将不同的估计信号s_i对应到相应的奇异值σ_i，其中的奇异值σ_i越大，期望信号和需要得到的消噪后信号近似程度越高。

本发明与现有技术相比，有益效果在于：

本发明提出基于奇异谱分析(Singular Spectrum Analysis，SSA)的瞬变电磁探测数据噪声滤除方法，解决了地质勘探中的实测瞬变电磁数据被噪声干扰的问题，采用粒子群优化算法求取奇异谱分析中的参数，弥补了参数选取依靠经验的不足之处，噪声滤除后的数据能够清晰地反映地层分布，同时能够提取噪声，为后续的噪声分析做准备。

附图说明

图1是本发明实施例提供的基于奇异谱分析的瞬变电磁探测数据噪声滤除方法流程框图；

图2是不同窗口选择对信号分离程度的影响；

图3为粒子群优化算法适应度曲线窗口长度(a)和均方根误差(b)；

图4不同含噪情况的TEM信号及其奇异谱；图4(a)为含白噪声为主的信号，图4(b)为含工频噪声为主的信号，图4(c)实测信号，图4(d)、(e)、(f)为对应的奇异谱图；

图5是加入工频噪声后的瞬变电磁信号；

图6是奇异谱分析后分离出的信号(a)和噪声(b)；

图7是对工频噪声消噪前后对比图；

图8为不同消噪算法对实测噪声消噪对比图。

具体实施方式

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

一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法，参见图1所示，包括：

1)本实施例中将理想的瞬变电磁信号加入工频噪声得到如图5所示的含噪瞬变电磁数据，验证本发明的有益效果，对含噪瞬变电磁数据重新排列数据。然后重新嵌入成轨迹矩阵Y；

2)计算轨迹矩阵Y的自协方差矩阵C，并进行奇异值分解(SVD)：

C＝U∑V^T (1)

式中U、V分别为左右特征向量矩阵，并得到得其特征值向量λ，特征值向量表示如下：

λ＝[λ₁，λ₂，…，λ_M]，|λ₁|≥|λ₂|≥…≥|λ_M| (2)

同时，矩阵C的奇异值向量σ可以表示成如下形式(其中M为窗口长度)：

3)选取奇异谱分析法中的窗口长度M：

3)a：设有m个粒子在一个S维的搜索空间中组成群体，第i个粒子表示为一个S维的向量X_i，将X_i代入目标函数计算出对应的适应值。

3)b：令第i个粒子当前迭代的最优个体极值为P_i＝(p_i1，p_i2，...，p_iS)，速度为V_i＝(v_i1，v_i2，...，v_iS)，粒子群当前迭代的最优全局极值为P_g＝(p_g1，p_g2，…，p_gS)，利用公式(4)、(5)更新S维上第i个粒子的速度和位置。

V_i(t+1)＝ωV_i(t)+c₁v₁(P_i(t)-X_i(t))+c₂v₂(P_g(t)-X_i(t)) (4)

X_i(t+1)＝X_i(t)+V_i(t+1) (5)

其中t为迭代次数；ω为惯性权重，一般取ω∈[0.4，1.0]；c₁和c₂为学习因子，一般取c₁，c₂∈[0，2]；r₁和r₂为相互独立且服从[0,1]上均匀分布的伪随机数。

3)c：设置参数为ω＝0.7，进化因子c₁＝c₂＝2，种群规模为20，同时限定初始种群位置范围[10,1000]，重复步骤II，经粒子群优化算法迭代计算50次，得到如图3所示的窗口长度(a)及均方根误差曲线(b)。由图3可以看出，算法收敛在窗口长度为150，与图2结果一致。因此选择窗口长度M＝150。

4)获取重构阶数的过程可以采用二分法、均值截断法等重构信号，采用含噪仿真和实测数据进行了研究。图4给出了三组不同TEM信号及其奇异谱，其中图4(a)是含白噪声为主的仿真TEM信号；图4(b)是含工频噪声为主的仿真TEM信号；图4(c)是实测信号，图4(d)(e)(f)为对应的奇异谱图。可以看出瞬变电磁信号的奇异值主要集中在前两个，因此在采用SSA算法对地空瞬变电磁数据的消噪中重构阶数选择为2。

5)利用步骤2)分解后的特征值对应的特征向量表示出轨迹矩阵Y的第j个分量轨迹矩阵：

同时将轨迹矩阵分量Y_j分成p个不相交的子集合J₁，J₂，…J_p，则合成矩阵Y_j表示如下：

此时计算各分量对应的合成矩阵，表示出原始矩阵Y：

6)将数据重新分组后进行对角平均化，轨迹矩阵被转换成时间道为N的一维数据。令

为L×K的轨迹矩阵，根据对角平均化公式将

转换为一维时间序列；

将步骤6)中不同的估计信号s_i对应到相应的奇异值σ_i，其中的奇异值σ_i越大，期望信号和需要得到的消噪后信号近似程度越高。

利用奇异谱分析加入工频噪声的数据得到分离出的信号图6(a)和噪声图6(b)；消噪后如图7所示，实现瞬变电磁信号与噪声的分离。

另一实施例中与上述的实施例不同之处在于，利用不同消噪算法对实测瞬变电磁信号进行去噪，效果对比如图8所示，由此可见奇异谱分析法去噪性能更高。

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

Claims (5)

1.一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法，其特征在于，该方法包括：

利用瞬变电磁探测系统获取一组实测含噪的瞬变电磁数据，对原始数据进行叠加；

将含噪的瞬变电磁数据重新嵌入成轨迹矩阵Y，并利用奇异值分解法进行分解得到特征值向量λ和奇异值向量σ；

提取奇异值向量σ中明显大于零的奇异值作为有效奇异值，根据有效奇异值重构信号并估计期望信号，将不同的估计信号对应到相应的奇异值中；

其中根据噪声滤除的要求，采用粒子群优化算法求取奇异谱分解法中的两个重要参数：窗口长度M与重构信号时奇异值的阶数。

2.按照权利要求1所述的方法，其特征在于，将含噪瞬变电磁探测数据重新嵌入，并利用奇异值分解法进行分解，具体包括：

将一维瞬变电磁数据y＝[y(1)，y(2)，...，y(N)]映射形成M×L轨迹矩阵Y，并重新排列；其中M为窗口长度，L＝N-M+1；

计算轨迹矩阵Y的自协方差矩阵C＝YY^T，并进行奇异值分解C＝U∑V^T，U为左特征向量矩阵、V为右特征向量矩阵，并得到其特征值向量λ，特征值向量表示如下：

λ＝[λ₁，λ₂，…，λ_M]，|λ₁|≥|λ₂|≥…≥|λ_M| (1)

同时，自协方差矩阵C的奇异值向量σ表示成如下形式：

3.按照权利要求2所述的方法，其特征在于，利用粒子群优化算法获取窗口长度M，具体包括：

设有m个粒子在一个S维的搜索空间中组成群体，第i个粒子表示为一个S维的向量X_i，将X_i代入目标函数计算出对应的适应值；

令第i个粒子当前迭代的最优个体极值为P_i＝(p_i1，p_i2，...，p_iS)，速度为V_i＝(v_i1，v_i2，...，v_iS)，粒子群当前迭代的最优全局极值为P_g＝(p_g1，p_g2，...，p_gS)，利用公式(3)、公式(4)更新S维上第i个粒子的速度和位置：

V_i(t+1)＝ωV_i(t)+c₁v₁(P_i(t)-X_i(t))+c₂v₂(P_g(t)-X_i(t)) (3)

X_i(t+1)＝X_i(t)+V_i(t+1) (4)

其中t为迭代次数；ω为惯性权重，取ω∈[0.4，1.0]；c₁和c₂为学习因子，取c₁，c₂∈[0，2]；r₁和r₂为相互独立且服从[0，1]上均匀分布的伪随机数。

4.按照权利要求1所述的方法，其特征在于，获取重构阶数的过程采用二分法及均值截断法重构信号，在对瞬变电磁数据重构时重构阶数选择为2。

5.按照权利要求1所述的方法，其特征在于，重构信号包含以下步骤：

利用分解后的特征值对应的特征向量表示出轨迹矩阵Y的第j个分量轨迹矩阵：

将轨迹矩阵分量Y_j分成p个不相交的子集合J₁，J₂，…J_p，则合成矩阵Y_j表示如下：

计算各分量对应的合成矩阵，表示出原始矩阵Y：

重新分组后进行对角平均化，轨迹矩阵被转换成时间道为N的一维数据，令

为L×K的轨迹矩阵，根据对角平均化公式将

转换为一维时间序列；

将不同的估计信号s_i对应到相应的奇异值σ_i，其中的奇异值σ_i越大，期望信号和需要得到的消噪后信号近似程度越高。

Images