CN109885906B

CN109885906B - 一种基于粒子群优化的磁共振测深信号稀疏消噪方法

Info

Publication number: CN109885906B
Application number: CN201910088618.0A
Authority: CN
Inventors: 田宝凤; 王亮; 王子强; 刘健楠
Original assignee: Jilin University
Current assignee: Jilin University
Priority date: 2019-01-30
Filing date: 2019-01-30
Publication date: 2022-05-13
Anticipated expiration: 2039-01-30
Also published as: CN109885906A

Abstract

本发明为一种基于粒子群优化的磁共振测深信号稀疏消噪方法，主要面向处理磁共振信号中的工频谐波干扰和随机白噪声。首先，对磁共振测深找水仪采集到的MRS信号利用带通滤波的方式对其进行预处理，通过频谱分析获得采集信号中含有的工频谐波干扰及MRS信号的频率，分别提出了构建针对MRS信号和工频谐波噪声特征的振荡原子库；然后，采用粒子群算法记录个体极值和群体极值来更新粒子群中各粒子的速度与位置，在工频谐波振荡原子库中挑选最佳原子重构工频以去除谐波干扰；最后，利用粒子群算法在MRS信号振荡原子库挑选出最佳原子重构MRS信号，若不满足实验要求，则计算残差信号，反复迭代直到符合条件。本发明不仅构建了针对MRS信号的新型原子库，且实现了含噪MRS信号中工频谐波干扰和随机白噪声的有效滤除，与传统MRS信号去噪方法相比，具有运算速度快、精度高，实用性强等优点。

Description

一种基于粒子群优化的磁共振测深信号稀疏消噪方法

技术领域

本发明属于磁共振测深(Magnetic Resonance Sounding,MRS)信号噪声滤除领域，具体是一种基于粒子群优化的磁共振测深信号稀疏消噪方法。

背景技术

与探地雷达、电磁法和电阻率测深等只能提供含水层间接岩性信息的地球物理方法相比，磁共振测深(magnetic resonance sounding，MRS)作为目前国际上唯一一种直接探测地下水的地球物理勘探手段，能够对地下水含量、储水层位置和厚度、地下介质孔隙度等信息给出定量解释，因而被广泛应用在水资源勘查、评价以及堤坝渗漏、矿井/隧道突涌水等领域。

由于MRS信号十分微弱，为纳伏级。在实际的工程应用中，采集到的MRS信号，经常和环境中的随机噪声、尖峰噪声和工频谐波噪声混合在一起，影响了我们需要获取的MRS信号的质量。针对MRS信号中的噪声滤除，国内外专家学者都进行了一些研究。从最初的“8字”型线圈和叠加的硬件方法提升信噪比，但提升幅度有限。到后来采用算法对采集到的含噪MRS信号进行噪声滤除。2003年Legchenko和Valla在论文《Removal of power-lineharmonics from proton magnetic resonance measurements》(《Journal of AppliedGeophysics》)中提出了采用区块对消法消除工频谐波噪声。2014年，Larsen等人在论文《Noise cancelling of MRS signals combining model-based removal of powerlineharmonics and multichannel Wienerfiltering》(《Geophysical JournalInternational》)中提出基于工频谐波建模和多通道维纳滤波相结合的方法来进行噪声抑制。但是上述两种方法主要是对工频谐波噪声进行滤除，不能消除随机噪声。2014年，Ghanati等人在论文《Filtering and parameter estimation of surface-NMR datausing singularspectrum analysis》(《Journal of Applied Geophysics》)中提出采用奇异谱分析的方法进行MRS噪声抑制及MRS信号提取，在去除工频谐波噪声和随机噪声方面取得了不错的效果。但是在MRS信号附近含有能量较大的随机噪声时，奇异谱分析算法效果不太理想。虽然以上算法在MRS信号噪声抑制方面取得了不错的效果，但是由于环境中含有噪声的以及机器采集时自身的不确定性干扰，使得以上的方法都具有一定的局限性。

专利CN105277973A公开了“一种基于匹配追踪的子波分解优化方法”，属于灾害预防领域；专利CN105513056A公开了“车载单目红外相机外参自动标定方法”，属于图像处理领域；专利CN106291677A公开了一种“一种基于匹配追踪方法的叠后声波阻抗反演方法”，属于油气开发领域；专利CN108507789A公开了“基于平均随机弱正交匹配追踪的滚动轴承故障稀疏诊断方法”，属于故障检测领域。可见匹配追踪算法已被成功应用到了信号处理的各个领域，但尚未见其应用于MRS信号的噪声滤除中。

目前，在利用匹配追踪处理信号的过程中，由于算法是建立在过完备原子库上，所以在寻找信号时需要和原子库内所有原子进行适应度计算，由于原子库的庞大导致每一轮计算都需要很长的时间，这导致该算法在很多场合的不适用。此类问题无疑将限制匹配追踪算法在磁共振测深信号噪声滤除中的应用。此外，如若含噪MRS信号经匹配追踪算法后，重构MRS信号得不到准确估计，则会导致后续反演解释中地层含水量估计错误率增高。

发明内容

针对直接应用传统匹配追踪算法所导致的运行效率低下、计算精度粗糙、MRS信号的信息不完整等不足，本发明提供了一种基于粒子群优化的磁共振测深信号稀疏消噪方法，该方法不仅能够实现多个工频谐波及随机白噪声的有效滤除，而且有效加快了匹配追踪的寻优过程，继而通过最终得到的最佳原子重构出MRS信号。

本发明是这样实现的，

一种基于粒子群优化的磁共振测深信号稀疏消噪方法，包括以下步骤：

步骤1：对磁共振测深探水仪采集到的一组观测MRS信号X(t)，利用带通滤波的方式对其进行预处理，得到目标频带范围内含噪MRS信号x(t)；

步骤2：建立工频谐波振荡原子库，采用粒子群算法记录个体极值和群体极值来更新粒子群中各粒子的速度与位置，优化匹配追踪挑选工频谐波特征的最佳原子，重构工频谐波干扰，并将其从含噪MRS信号x(t)中去除；

步骤3：构建MRS信号振荡原子库，采用粒子群算法记录个体极值和群体极值来更新粒子群中各粒子的速度与位置，优化匹配追踪挑选MRS信号特征的最佳原子，重构MRS信号，直到满足误差精度要求停止迭代，实现将随机噪声从含噪MRS信号x(t)的去除。

进一步地，采用粒子群算法优化匹配追踪挑选工频谐波特征的最佳原子以及采用粒子群算法优化匹配追踪挑选MRS信号特征的最佳原子包括：

根据MRS信号特点以及工频谐波噪声特点建立振荡原子库表达式：

其中，f为振荡频率，

为振荡相位，范围为[0,2π/rad]，ρ为振荡衰减系数，ρ≥0。当ρ＞0时为衰减振荡，对应MRS信号；当ρ＝0时为谐波干扰，对应工频谐波干扰；将g_r归一化后得到系数K_r；

将参数组

中的各参数进行离散化处理，其中w＝2πf，f为振荡频率，

为振荡相位，ρ为振荡衰减系数；

得到原子各参数和离散化后参数的对应关系；

得到第一离散参数p与振荡频率f之间关系；

得到第二离散参数r与平均弛豫时间

之间的关系。

进一步地，参数组

中的各参数进行离散化处理，各参数按照以下式(2)进行离散化:

进一步地，原子各参数和离散化后参数的对应关系为式(3)：

进一步地，第一离散参数p与信号频率f之间关系式为：

进一步地，第二离散参数r与

之间的关系式为：

进一步地，所述粒子群算法，采用由Rberhart和Kennedy从生物学角度模拟鸟群捕食的过程所提出的寻优算法，具体步骤为：

初始化粒子群中各粒子的速度与位置；

计算粒子群中最适应个体；

记录个体极值与群体极值；

更新各粒子的速度与位置；

判断是否满足粒子群算法终止条件；

若是，则结束；

若否则重新计算粒子群中最适应个体。

进一步地，采用粒子群算法优化匹配追踪挑选MRS信号特征的最佳原子包括：将数据化后的鸟群作为MRS信号的原子数据群，设在4维度空间飞行，设MRS信号的原子数据群的四个参数分别为S₁、S₂、S₃和S₄，鸟群整体的飞行空间的大小为：

S＝S₁·S₂·S₃·S₄ (6)

初始化时，原子在数据范围内随机取值。

进一步地，将鸟群作为原子，鸟群速度为算法运行的参数值，鸟群位置为原子的评价值，食物为信号；越靠近食物即原子和信号越相似，评价值越高，反之则越小，其中，评价值与信号的关系如式(7)所示：

Value＝|<Signal,Atom>| (7)

式(7)中，Value表示评价值，Signal表示信号，Atom表示原子；

在粒子群算法中，鸟群的的飞行速度由评价值，经验值，随机参数共同决定的，通过式(8)、式(9)得到最适应个体：

v_id＝ω·v_id+c₁·e₁·(p_id-x_id)+c₂·e₂·(p_gd-x_id) (8)

x_id＝x_id+v_id (9)

式(8)与(9)中，v_id表示飞行速度，在飞行空间内速度越大，位置的改变就越大，飞行的单位时间为程序运行的次数。c₁和c₂表示惯性权值，代表了上一次的飞行速度对于下一次飞行速度的影响，权值越大在大范围的参数中寻找的最优个体，权值越小，则搜寻的结果越精准；e₁和e₂表示[0-1]的随机数，使所有时避免陷入局部最优；x_id表示原子位置，为此时所选择取的参数，；p_id表示个体极值，为在一个原子的k次飞行中的最佳位置，p_gd表示群体极值，代表整个飞行粒子群在k次飞行中的最佳位置。

本发明与现有技术相比，有益效果在于：本发明提出了一种基于粒子群优化的磁共振测深信号稀疏消噪方法，不仅实现了MRS信号中多个工频谐波干扰以及随机白噪声的有效滤除，而且该方法在建立振荡原子库的基础上，利用粒子群算法有效加快了匹配追踪的寻优过程，无需按照传统MP算法的贪婪思想，在寻找信号时需要和原子库内所有原子进行适应度计算，使得算法的运行效率和结果的计算精度均得到了加强。与传统MRS信号去噪方法相比，本发明方法能够减少MRS信号信息缺失，运算速度快，信噪比高。此外，本发明方法还开辟了粒子群算法优化匹配追踪在磁共振信号消噪领域的新天地，为其后续的的应用拓宽了思路。

附图说明

图1本发明实施例提供的基于粒子群优化的磁共振测深信号稀疏消噪方法的流程框图；

图2本发明实施例提供的粒子群算法流程框图；

图3理想MRS信号(图3A)及其频谱(图3B)；

图4含噪MRS信号(图4A)及其频谱(图4B)；

图5本发明实施例提供的重构工频(图5A)及其频谱(图5B)；

图6本发明实施例提供的重构MRS信号(图6A)及其频谱(图6B)；

图7本发明实施例提供的实测MRS信号(图7A)及其频谱(图7B)；

图8本发明实施例提供的带通滤波后MRS信号(图8A)及其频谱(图8B)；

图9本发明实施例提供的实测数据重构工频(图9A)及其频谱(图9B)；

图10本发明实施例提供的实测数据重构MRS信号(图10A)及其频谱(图10B)。

具体实施方式

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

如图1所示，一种基于粒子群优化的磁共振测深信号稀疏消噪算法，包括以下步骤：

步骤1：对磁共振测深(MRS)探水仪采集到的一组观测MRS信号X(t)，利用带通滤波的方式对其进行预处理，得到目标频带范围内含噪MRS信号x(t)；

步骤2：建立工频谐波振荡原子库，采用粒子群算法优化匹配追踪挑选工频谐波特征的最佳原子，重构工频谐波干扰，并将其从x(t)中去除。

步骤3：构建MRS信号振荡原子库，采用粒子群算法优化匹配追踪挑选MRS信号特征的最佳原子，重构MRS信号，直到满足误差精度要求停止迭代，实现将随机噪声从x(t)的去除；

其中，振荡原子库建立方法如下：

当处理对象为MRS信号，其表达式为：

其中，E₀表示初始振幅，代表地下水含量的大小；T₂ ^*表示平均弛豫时间，代表地下介质孔隙度尺寸的大小；f_l表示拉莫尔频率，与地磁场有关，它在全球的变化范围大致为1.3kHz～3.7kHz；θ表示初始相位，代表地下介质的导电性。由公式(10)可见，MRS信号是呈e指数衰减特性的余弦信号。

工频谐波噪声是由电力线产生的一组正弦波，表达式为：

其中，A_n和

分别是第n个谐波的幅度和相位；N是谐波个数，f₀为工频基频，理论上f₀＝50Hz，但存在一定程度的波动。结合公式(10)和(11)，提出采用振荡原子库：

式(1)中，f为振荡频率，

为振荡相位，范围为[0,2π/rad]，ρ为振荡衰减系数，ρ≥0。当ρ＞0时为衰减振荡，此时对应MRS信号；当ρ＝0时为谐波干扰，此时对应工频谐波干扰；将g_r归一化后得到系数K_r。

在进行原子库构建时，需要将参数组

中的各参数进行离散化处理，其中w＝2πf，各参数按照以下方式进行离散化:

原子各参数和离散化后参数的对应关系为：

化简式(2)和(3)得到离散参数p与信号频率f之间关系式为：

离散参数r与

之间的关系式为：

如图2，本发明采用由Rberhart和Kennedy从生物学角度模拟鸟群捕食的过程所提出的粒子群算法来寻找最佳原子，具体步骤为：具体步骤为：

初始化粒子群中各粒子的速度与位置；

计算粒子群中最适应个体；

记录个体极值与群体极值；

更新各粒子的速度与位置；

判断是否满足粒子群算法终止条件；

若是，则结束；

若否则重新计算粒子群中最适应个体。

对MRS信号包括：

1)、在鸟群捕食最初时，鸟群应该是随机分布在一个三维空间内，且都具有一个随机速度。在将鸟群捕食数据化的工程中，首先鸟群即为初代原子数据群，鸟群所处的三维空间就相当于是三个参数变量，而对MRS信号来讲，由于MRS信号为4参数变量，所以数据化后的鸟群即原子数据群，理论上是在一个4维度空间飞行。假设四个参数分别为S₁、S₂、S₃和S₄，那么鸟群整体的飞行空间的大小为：

S＝S₁·S₂·S₃·S₄ (6)

在初始化时，鸟群随机在整个空间内分布，相当于数据化之后，原子在这个数据范围内随机取值。

2)、在鸟群捕食的过程中，鸟群的飞行策略是向距离食物最近的鸟飞行。主要具有两个属性：速度与位置。速度是具有惯性的，且收到位置的影响，即如果位置好，那么根据惯性继续飞行，如果位置不好，那么加速或减速向更好的位置飞行。距离食物的远近即为优劣的判断条件。将上面的属性与飞行策略数据化，即鸟群为原子，鸟群速度为算法运行的参数值，鸟群位置为原子的评价值，食物即为信号。越靠近食物即原子和信号越相似，那么评价值越高，反之则越小。其中，评价值与信号的关系如式(7)所示。

Value＝|<Signal,Atom>| (7)

式(7)中，Value表示评价值，Signal表示信号，Atom表示原子。

3)、在粒子群算法中，鸟群的的飞行速度是由评价值，经验值，随机参数等共同决定的。具体公式如式(8)、式(9)所示。

v_id＝ω·v_id+c₁·e₁·(p_id-x_id)+c₂·e₂·(p_gd-x_id) (8)

x_id＝x_id+v_id (9)

式(8)与(9)中，v_id表示飞行速度，在飞行空间内速度越大，位置的改变就越大，飞行的单位时间即为程序运行的次数。c₁和c₂表示惯性权值，代表了上一次的飞行速度对于下一次飞行速度的影响。理论上，权值越大可以在大范围的参数中寻找的最优个体，权值越小，则搜寻的结果越精准。e₁和e₂表示[0-1]的随机数，可以使所有时避免陷入局部最优。x_id表示原子位置，即为此时所选择取的参数。p_id表示个体极值，它的含义是在一个原子的k次飞行中的最佳位置，p_gd表示群体极值，代表整个飞行粒子群在k次飞行中的最佳位置。

实施例1

本实施例是在MATLAB 7.0编程环境下开展的本发明方法的仿真实验。

基于粒子群优化的磁共振测深信号稀疏消噪方法的仿真算法，参照图1，包括以下步骤：

步骤(1)：利用式

构造拉莫尔频率为2345Hz，幅值e₀为150nV，弛豫时间

为0.18s的理想MRS信号，如图3所示理想MRS信号(图3A)及其频谱(图3B)。在该信号拉莫尔频率附近添加2200Hz、2250Hz、2300Hz、2350Hz、2400Hz、2450Hz和2500Hz的工频干扰和幅值为100nV的随机噪声，经过某种线性组合形成信噪比为-6.7471dB的观测MRS信号x(t)(为行向量)，如图4所示含噪MRS信号(图4A)及其频谱(图4B)；

步骤(2)：建立工频谐波振荡原子库，对工频干扰重构后进行去除。根据图2所描述的粒子群算法分别挑选工频谐波2200Hz、2250Hz、2300Hz、2350Hz、2400Hz、2450Hz和2500Hz相应的原子，并对其进行重构后从原始含噪信号中去除，如图5所示重构工频(图5A)及其频谱(图5B)；

步骤(3)：针对经步骤(2)去除工频后的残余信号中含有MRS信号和随机噪声，对MRS信号表达式中的四个重要特征参数E₀、

f_L、

构建MRS信号振荡原子库，根据图2所描述的粒子群算法挑选MRS信号相应的原子并重构MRS信号，完成对随机噪声的去除，如图6所示重构MRS信号(图6A)及其频谱(图6B)；

为了验证本发明方法的实用性，将去噪后MRS信号s(t)进行了信噪比(SNR)估计。经计算，其SNR＝26.04dB，较分离前的SNR提高了44.21dB；接着对s(t)进行了包络提取和数据拟合，以获得分离信号的关键参数初始振幅E₀和弛豫时间T₂ ^*，计算可得，E₀＝149.34nV，T₂ ^*＝179.0s，相对误差分别为-0.44％、-0.53％，均控制在±5％以内，满足应用要求。

实施例2

本实施例是以长春市文化广场实地采集的实测MRS信号作为本发明方法的处理对象。

基于粒子群优化的磁共振测深信号稀疏消噪方法，如图1所示，包括以下步骤：

步骤(1)：对核磁共振测深(MRS)探水仪采集到的一组观测MRS信号X(t)，如图7所示实测MRS信号(图7A)及其频谱(图7B)；利用带通滤波的方式对其进行预处理，得到目标频带范围内含噪MRS信号x(t)(为行向量)，如图8A所示，计算其信噪比为SNR＝-7.11dB；对其进行傅里叶变换得到其频谱，如图8B所示，可以看到该信号在f₁＝2300Hz、f₂＝2400Hz处均有较强的工频谐波干扰，此外，在f₃＝2311Hz、f₄＝2347Hz、f₅＝2368Hz处受到了强烈的单频干扰；

步骤(2)：建立工频谐波振荡原子库，根据图2所描述的粒子群算法挑选2300Hz、2400Hz的工频谐波以及2311Hz、2347Hz和2368Hz三个单频干扰的原子，对工频谐波和单频干扰进行重构后从原始含噪信号中去除，如图9所示给出了重构的工频谐波和单频干扰(图9A)及其频谱(图9B)；

f_L、

构建MRS信号振荡原子库，根据图2所描述的粒子群算法挑选MRS信号相应的原子并重构MRS信号，实现随机噪声的去除，得到最终去噪后信噪比为SNR＝14.82dB的目标MRS信号s(t)，如图10所示实测数据重构MRS信号(图10A)及其频谱(图10B)。

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

Claims

1.一种基于粒子群优化的磁共振测深信号稀疏消噪方法，其特征在于，包括以下步骤：

步骤3：构建MRS信号振荡原子库，采用粒子群算法记录个体极值和群体极值来更新粒子群中各粒子的速度与位置，优化匹配追踪挑选MRS信号特征的最佳原子，重构MRS信号，直到满足误差精度要求停止迭代，实现将随机噪声从含噪MRS信号x(t)的去除；

采用粒子群算法优化匹配追踪挑选工频谐波特征的最佳原子以及采用粒子群算法优化匹配追踪挑选MRS信号特征的最佳原子包括：

其中，f为振荡频率，

为振荡相位，范围为[0,2π/rad]，ρ为振荡衰减系数，ρ≥0，当ρ＞0时为衰减振荡，对应MRS信号；当ρ＝0时为谐波干扰，对应工频谐波干扰；将g_r归一化后得到系数K_r；

将参数组

中的各参数进行离散化处理，其中w＝2πf，f为振荡频率，

为振荡相位，ρ为振荡衰减系数；

得到原子各参数和离散化后参数的对应关系；

得到第一离散参数p与振荡频率f之间关系；

得到第二离散参数r与平均弛豫时间

之间的关系。

2.按照权利要求1所述的方法，其特征在于，参数组

3.按照权利要求2所述的方法，其特征在于，原子各参数和离散化后参数的对应关系为式(3)：

4.按照权利要求1所述的方法，其特征在于，第一离散参数p与信号频率f之间关系式为：

5.按照权利要求3所述的方法，其特征在于，第二离散参数r与

之间的关系式为：

6.按照权利要求1所述的方法，其特征在于，所述粒子群算法，采用由Rberhart和Kennedy从生物学角度模拟鸟群捕食的过程所提出的寻优算法，具体步骤为：

初始化粒子群中各粒子的速度与位置；

计算粒子群中最适应个体；

记录个体极值与群体极值；

更新各粒子的速度与位置；

判断是否满足粒子群算法终止条件；

若是，则结束；

若否则重新计算粒子群中最适应个体。

7.按照权利要求6所述的方法，其特征在于，采用粒子群算法优化匹配追踪挑选MRS信号特征的最佳原子包括：将数据化后的鸟群作为MRS信号的原子数据群，设在4维度空间飞行，设MRS信号的原子数据群的四个参数分别为S₁、S₂、S₃和S₄，鸟群整体的飞行空间的大小为：

S＝S₁·S₂·S₃·S₄ (6)

初始化时，原子在数据范围内随机取值。

8.按照权利要求6或7所述的方法，其特征在于，将鸟群作为原子，鸟群速度为算法运行的参数值，鸟群位置为原子的评价值，食物为信号；越靠近食物即原子和信号越相似，评价值越高，反之则越小，其中，评价值与信号的关系如式(7)所示：

Value＝|<Signal,Atom>| (7)

式(7)中，Value表示评价值，Signal表示信号，Atom表示原子；

在粒子群算法中，鸟群的飞行速度由评价值，经验值，随机参数共同决定的，通过式(8)、式(9)得到最适应个体：

v_id＝ω·v_id+c₁·e₁·(p_id-x_id)+c₂·e₂·(p_gd-x_id) (8)

x_id＝x_id+v_id (9)

式(8)与(9)中，v_id表示飞行速度，在飞行空间内速度越大，位置的改变就越大，飞行的单位时间为程序运行的次数；c₁和c₂表示惯性权值，代表了上一次的飞行速度对于下一次飞行速度的影响，权值越大在大范围的参数中寻找的最优个体，权值越小，则搜寻的结果越精准；e₁和e₂表示[0-1]的随机数，使所有时避免陷入局部最优；x_id表示原子位置，为此时所选择取的参数；p_id表示个体极值，为在一个原子的k次飞行中的最佳位置，p_gd表示群体极值，代表整个飞行粒子群在k次飞行中的最佳位置。