一种联合EMD与TFPF算法的全波磁共振信号随机噪声消减方法

摘要

本发明涉及一种联合EMD与TFPF算法的全波磁共振信号随机噪声消减方法，是一种不需要设计滤波区间的“盲”滤波方法。首先利用EMD算法的分解特性，将全波磁共振信号分解成不同本征模态分量，再使用TFPF算法将信号主导模态分量编码为单位幅度解析信号的瞬时频率，利用解析信号的时频分布沿着瞬时频率集中的特性来抑制随机噪声。该方法需要较少的滤波约束条件，操作简单，不需要在时频域设计滤波区间，对于低信噪比的全波磁共振信号具有较强的适应性。显著提高探测效率，一次测量即可获得较好的消噪效果，有效消减随机噪声的同时不损失信号成分，可显著增强信噪比，提高了后期反演解释的准确性。

说明书

技术领域：

本发明涉及一种核磁共振测深(Magnetic Resonance Sounding,MRS)信号噪声消减方法，尤其是利用基于EMD与TFPF联合算法原理消减全波磁共振信号中所含随机噪声的数据处理方法。

背景技术：

磁共振地下水探测技术是一项新兴的地球物理探测技术，由于其具有直接探测和定量评估的优点，近年来受到了广泛的关注，并取得了快速的发展。目前，已经由理论研究逐步发展到研制比较成熟的仪器探测系统，并在一些地区进行了实地测量，成功探测到浅层地下水资源，对改善目前国家贫水现状以及预测由地下水引起的地质灾害具有重要的意义。但由于地下水中受激氢质子弛豫产生的磁共振信号幅度很小(纳伏级1nV＝10^-9V)，极易受到各种环境噪声干扰。对于随机噪声干扰，由于其频带分布范围广、与有效信号交织在一起，难以去除，严重影响了信号的信噪比，导致在后期进行反演解释时出现较大误差。因此，有效的消除随机噪声获取可靠的磁共振信号将为提高反演解释结果的准确性奠定重要的基础。

专利CN104459809A公开了一种“基于独立成分分析的全波核磁共振信号噪声滤除方法”，主要针对全波核磁共振信号中的工频谐波干扰或某一单频干扰。首先利用核磁共振测深探水仪采集MRS信号，通过频谱分析获得采集信号中含有的工频谐波干扰或某一单频干扰的频率，采用数字正交法构造输入通道信号解决欠定盲源分离问题；然后，将构造的输入通道信号与采集的MRS信号一并作为输入信号进行独立成分分析，得到分离MRS信号；最后采用频谱校正法解决分离后MRS信号幅值不定问题，进而提取去噪后MRS信号。专利CN104898172A公开了一种“基于互相关的核磁共振全波信号噪声滤除方法”，是利用噪声与拉莫尔频率的正弦信号不相关，而FID幅度衰减正弦信号与拉莫尔频率的正弦信号具有相关性的特点，通过互相关运算滤除噪声，然后拟合互相关波形的包络，并重构不含噪声的互相关波形，最后利用解卷积算法提取核磁共振全波数据中的FID信号。专利CN104614778A公开了一种“基于ICA的核磁共振地下水探测信号噪声消除方法”，该方法首先录入三组核磁共振响应数据，分别对这三组数据进行傅里叶变换，确定每组数据核磁共振中心频率附近所含工频谐波，然后构造与工频谐波同频率，与核磁共振相应数据同长度的正弦函数、余弦函数，并与核磁共振响应数据组成观测信号，采用独立分量分析算法对每组观测信号进行分离得到解混信号，进行数据重构以消除工频谐波的干扰，将三组去除工频谐波的核磁共振数据作为观测信号，再利用ICA算法处理，削弱剩余随机噪声的干扰。蒋川东在Near Surface Geophysics[2011,9(5),459-468]上发表的论文“Statistical stacking and adaptive notch filter to remove high-level electromagnetic noise from MRS measurements”中采用叠加方法抑制随机噪声，利用磁共振地下水探测系统多次(N_s次)采集信号，将所有单次采集的信号进行叠加处理，信噪比可提高倍。

上述发明的基于独立成分分析的全波核磁共振信号噪声滤除方法，可解决独立成分分析中惯有的欠定盲源问题和幅值不定问题，与传统MRS信号去噪方法相比，具有运算速度快、信噪比高、实用性强等优点，但该方法只能实现工频干扰和某一频率干扰的滤除，不能抑制随机噪声这类频带分布范围大的噪声干扰；而基于互相关的核磁共振全波信号噪声滤除方法可同时压制工频及其谐波噪声、随机噪声和尖峰噪声，但主要利用噪声频率和信号频率不相关的特性，对于随机噪声中频率成分和信号频率重叠的部分无法去除；基于ICA的核磁共振地下水探测信号噪声消除方法，对源信号和传输通道的先验知识没有要求，在试验过程中无需铺设参考线圈，操作简单、方便、效率高，但该方法要求数据至少为三组以上，且计算过程复杂，计算量大，对于非技术人员的操控难度大；叠加法是目前最简单且普遍采用的磁共振信号随机噪声消减方法，在相同强度的电磁噪声环境下，利用较少的叠加次数可获取较高信噪比的地下水磁共振信号，但当噪声强度大时需要进行多次叠加，将增加测量的时间，导致仪器的工作效率降低，且消噪效果有限。

发明内容：

本发明的目的就在于针对上述现有技术的不足，提供一种联合经验模态分解(Empirical Mode Decomposition,EMD)与时频峰值滤波(Time-frequency peak filtering,TFPF)算法的全波磁共振信号随机噪声消减方法，有效消减随机噪声的同时不损失信号成分，该方法对受随机噪声干扰的低信噪比全波磁共振信号具有较强的适应性，一次探测即可获得较好的消噪效果，是能够全面提高探测效率的一种联合EMD与TFPF算法的全波磁共振信号随机噪声消减方法。

本发明目的是通过以下技术方案实现的：

一种联合EMD与TFPF算法的全波磁共振信号随机噪声消减方法，包括以下步骤：

a、对核磁共振地下水探测仪采集的全波磁共振观测信号X(n)进行EMD分解(n为离散样本点)，从高频到低频分解出i个不同的本征模态分量C₁(n)，…，C_i(n)以及一个趋势项R_i(n)，X(n)＝C₁(n)+L+C_i(n)+R_i(n)；

b、对本征模态分量C₁(n)，…，C_i(n)以及趋势项R_i(n)进行傅里叶变换和自相关分析，综合傅里叶变换结果和自相关分析结果得到噪声主导的模态分量和信号主导的模态分量，提取出信号主导的模态分量c₁(n)，…，c_j(n)(j为信号主导的模态分量个数)；

c、使用TFPF方法对信号主导的模态分量c₁(n)，…，c_j(n)分别进行处理，消除各模态分量中所含的随机噪声，得到信号分量s₁(n)，…，s_j(n)；

d、对处理后得到的信号分量s₁(n)，…，s_j(n)进行相加，得到最终不含随机噪声干扰的全波磁共振信号s(n)，s(n)＝s₁(n)+L+s_j(n)。

所述步骤a中的EMD算法的具体步骤为：

首先，识别观测信号X(n)的所有极大值点和所有极小值点，分别对所有极大值点和所有极小值点进行三次样条插值得到数据的上包络曲线E_max(n)和下包络曲线E_min(n)；

其次，计算上、下包络曲线对应的平均值，得到均值曲线F₁(n)；

第三，将观测信号X(n)与平均值F₁(n)作差得到细节分量H₁(n)，判断该细节分量是否为本征模态函数，判断条件为：

①、函数的均值为零、局部对称，并且过零点个数和极值点个数相同或至多相差一个；

②、任意时刻函数上、下包络值之和恒为零；

如果不满足条件，继续对H₁(n)重复前两个步骤，直到满足条件，得到第一个本征模态分量C₁(n)；

第四，将第一个本征模态分量C₁(n)从观测信号X(n)中减去，得到残差量R₁(n)，即：

X(n)-C₁(n)＝R₁(n)

把该残差量R₁(n)作为新的待处理信号重复前三个步骤继续筛选，直到获得第二个本征模态分量C₂(n)，依此类推，得到i个不同的本征模态分量C₁(n)，…，C_i(n)以及一个趋势项R_i(n)：

R₁(n)-C₂(n)＝R₂(n)，…，R_i-1(n)-C_i(n)＝R_i(n)

X(n)＝C₁(n)+L+C_i(n)+R_i(n)；

所述步骤c中TFPF算法是对信号主导的模态分量进行处理，得到消除随机噪声后的信号分量，具体步骤包括：

1)、对信号主导的模态分量c₁(n)进行边缘数据延拓，得到延拓后的信号c₁'(m)，即

其中，N是原始观测信号的长度，p是两端延拓的数据点数；

2)、将边缘数据延拓后的信号c₁'(m)进行尺度变换：

其中，d₁'(m)是尺度变换后的信号，系数a和b分别是变换后信号的最大值和最小值，满足0.5≥a＝max[d₁'(m)]，b＝min[d₁'(m)]≥0；

3)、将尺度变换后信号d₁'(m)进行频率调制编码，得到单位幅度的解析信号：

其中，z₁(m)是频率调制编码后的解析信号；

4)、对频率调制编码后的解析信号z₁(m)进行离散伪Wigner-Ville变换，计算z₁(m)的时频分布：

其中，w(l)是窗函数，其宽度为2L+1；

5)、取W_z1(m,k)的峰值作为信号的有效估计，即

其中，arg_kmax[·]表示沿着频率取最大值算子；

6)、对有效信号的估计值进行反尺度变换，恢复信号幅度：

7)、舍去滤波后结果的拓展边缘得到消除随机噪声后的信号分量s₁(n)：

8)、接下来，按照上述步骤1)到步骤7)依次对信号主导的模态分量c₂(n)，…，c_j(n)进行处理，得到消除随机噪声后的信号分量s₂(n)，…，s_j(n)。

有益效果：本发明是一种不需要设计滤波区间的“盲”滤波方法，首先利用EMD算法的分解特性，将全波磁共振信号分解成不同的本征模态分量，再使用TFPF算法将信号主导的模态分量编码为单位幅度解析信号的瞬时频率，利用解析信号的时频分布沿着瞬时频率集中的特性来抑制随机噪声。有效消减随机噪声的同时不损失信号成分，对准确反演解释水文地质参数具有重要意义。该方法滤波约束条件较少，操作简单，不需要在时频域设计滤波区间，对于低信噪比的全波磁共振信号具有较强的适应性。显著提高了探测效率，一次测量即可获得较好的消噪效果，显著增强信噪比，提高了后期反演解释的准确性，具有较大的实际应用价值。

附图说明：

图1为联合EMD与TFPF算法的全波磁共振信号随机噪声消减方法流程图；

图2为含随机噪声的全波磁共振观测信号；

图3为EMD算法分解观测信号后得到的本征模态分量；

图4为本征模态分量的傅里叶变换；

图5为本征模态分量的自相关分析；

图6为信号主导的本征模态分量；

图7为TFPF算法消减随机噪声后的信号分量；

图8为消减随机噪声后的全波磁共振信号。

具体实施方式：

下面结合附图和实施例对本发明作进一步的详细说明：

一种联合EMD与TFPF算法的全波磁共振信号随机噪声消减方法，包括以下步骤：

a、对核磁共振地下水探测仪采集的全波磁共振观测信号X(n)进行EMD分解(n为离散样本点)，从高频到低频分解出i个不同的本征模态分量C₁(n)，…，C_i(n)以及一个趋势项R_i(n)，X(n)＝C₁(n)+L+C_i(n)+R_i(n)；

b、对本征模态分量C₁(n)，…，C_i(n)以及趋势项R_i(n)进行傅里叶变换和自相关分析，综合傅里叶变换结果和自相关分析结果得到噪声主导的模态分量和信号主导的模态分量，提取出信号主导的模态分量c₁(n)，…，c_j(n)(j为信号主导的模态分量个数)；

c、使用TFPF方法对信号主导的模态分量c₁(n)，…，c_j(n)分别进行处理，消除各模态分量中所含的随机噪声，得到信号分量s₁(n)，…，s_j(n)；

d、对处理后得到的信号分量s₁(n)，…，s_j(n)进行相加，得到最终不含随机噪声干扰的全波磁共振信号s(n)，s(n)＝s₁(n)+L+s_j(n)。

所述步骤a中的EMD算法的具体步骤为：

首先，识别观测信号X(n)的所有极大值点和所有极小值点，分别对所有极大值点和所有极小值点进行三次样条插值得到数据的上包络曲线E_max(n)和下包络曲线E_min(n)；

其次，计算上、下包络曲线对应的平均值，得到均值曲线F₁(n)；

第三，将观测信号X(n)与平均值F₁(n)作差得到细节分量H₁(n)，判断该细节分量是否为本征模态函数，判断条件为：

①、函数的均值为零、局部对称，并且过零点个数和极值点个数相同或至多相差一个；

②、任意时刻函数上、下包络值之和恒为零；

如果不满足条件，继续对H₁(n)重复前两个步骤，直到满足条件，得到第一个本征模态分量C₁(n)；

第四，将第一个本征模态分量C₁(n)从观测信号X(n)中减去，得到残差量R₁(n)，即：

X(n)-C₁(n)＝R₁(n)

把该残差量R₁(n)作为新的待处理信号重复前三个步骤继续筛选，直到获得第二个本征模态分量C₂(n)，依此类推，得到i个不同的本征模态分量C₁(n)，…，C_i(n)以及一个趋势项R_i(n)：

R₁(n)-C₂(n)＝R₂(n)，…，R_i-1(n)-C_i(n)＝R_i(n)

X(n)＝C₁(n)+L+C_i(n)+R_i(n)；

所述步骤c中TFPF算法是对信号主导的模态分量进行处理，得到消除随机噪声后的信号分量，具体步骤包括：

1)、对信号主导的模态分量c₁(n)进行边缘数据延拓，得到延拓后的信号c₁'(m)，即

其中，N是原始观测信号的长度，p是两端延拓的数据点数；

2)、将边缘数据延拓后的信号c₁'(m)进行尺度变换：

其中，d₁'(m)是尺度变换后的信号，系数a和b分别是变换后信号的最大值和最小值，满足0.5≥a＝max[d₁'(m)]，b＝min[d₁'(m)]≥0；

3)、将尺度变换后信号d₁'(m)进行频率调制编码，得到单位幅度的解析信号：

其中，z₁(m)是频率调制编码后的解析信号；

4)、对频率调制编码后的解析信号z₁(m)进行离散伪Wigner-Ville变换，计算z₁(m)的时频分布：

其中，w(l)是窗函数，其宽度为2L+1；

5)、取W_z1(m,k)的峰值作为信号的有效估计，即

其中，arg_kmax[·]表示沿着频率取最大值算子；

6)、对有效信号的估计值进行反尺度变换，恢复信号幅度：

7)、舍去滤波后结果的拓展边缘得到消除随机噪声后的信号分量s₁(n)：

8)、接下来，按照上述步骤1)到步骤7)依次对信号主导的模态分量c₂(n)，…，c_j(n)进行处理，得到消除随机噪声后的信号分量s₂(n)，…，s_j(n)。

联合EMD与TFPF算法的全波磁共振信号随机噪声消减方法的详细说明：

1)、如图1所示，读取核磁共振地下水探测仪采集的受随机噪声干扰的全波磁共振观测信号X(n)，观测信号见图2。

2)、识别观测信号X(n)的所有极大值点和所有极小值点，分别对所有极大值点和所有极小值点进行三次样条插值得到数据的上包络曲线E_max(n)和下包络曲线E_min(n)。

3)、计算上、下包络曲线对应的平均值，得到均值曲线F₁(n)，

4)、将观测信号X(n)与平均值F₁(n)作差得到细节分量H₁(n)，H₁(n)＝X(n)-F₁(n)，判断该细节分量是否为本征模态函数，判断条件为：

第一、函数的均值为零、局部对称，并且过零点个数和极值点个数相同或至多相差一个；

第二、任意时刻函数上、下包络值之和恒为零；

如果不满足条件，继续对H₁(n)重复步骤2)和步骤3)，直到满足条件，得到第一个本征模态分量C₁(n)。

5)、将第一个本征模态分量C₁(n)从观测信号X(n)中减去，得到残差量R₁(n)，即X(n)-C₁(n)＝R₁(n)。把该残差量R₁(n)作为新的待处理信号重复步骤2)～步骤4)继续筛选，直到获得第二个本征模态分量C₂(n)。依此类推，参见图3，得到十四个不同的本征模态分量C₁(n)，…，C₁₄(n)。

6)、对图3所示的本征模态分量C₁(n)，…，C₁₄(n)进行傅里叶变换和自相关分析，在Matlab中使用fft命令可得到本征模态分量的频谱，如图4所示；使用xcorr命令计算本征模态分量的归一化自相关函数，如图5所示。综合傅里叶变换结果和自相关分析结果得到噪声主导的模态分量和信号主导的模态分量，提取出三个信号主导的模态分量c₁(n)，c₂(n)和c₃(n)，如图6所示。

7)、使用TFPF方法对第一个信号主导的模态分量c₁(n)进行处理，该模态分量接下来会被编码为解析信号的瞬时频率，由于解析信号的时频分布边缘能量不集中，从时频分布的峰值恢复信号时会产生端点失真现象。因此，需要先进行边缘数据延拓，得到延拓后的信号c₁'(m)，即

其中，N是原始观测信号的长度，p是两端延拓的数据点数。

8)、为避免频率溢出，将信号调制为频率信息时，必须确保它位于限定的频率范围内，因此，需要将信号c₁'(m)进行尺度变换：

其中，d₁'(m)是尺度变换后的信号，系数a和b分别是变换后信号的最大值和最小值，满足0.5≥a＝max[d₁'(m)]，b＝min[d₁'(m)]≥0。

9)、将尺度变换后信号d₁'(m)进行频率调制编码，得到单位幅度的解析信号：

其中，z₁(m)是频率调制编码后的解析信号。

10)、对频率调制编码后的解析信号z₁(m)进行离散伪Wigner-Ville变换。计算z₁(m)的时频分布：

其中，w(l)是窗函数，其宽度为2L+1。信号经过Wigner-Ville变换后可看作是其能量在联合的时间域和频率域中的分布，是分析非平稳时变信号的重要工具。当信号淹没于噪声中，并且是时间的线性函数时，可得到信号的无偏估计，由于全波磁共振信号不是时间的线性函数，根据瞬时频率估计的方法，采用了加窗的Wigner-Ville变换，即伪Wigner-Ville变换计算时频分布，进行局部线性化处理，使信号在一个窗长内近似满足线性瞬时频率条件，减小信号非线性引起的偏差。

11)、取W_z1(m,k)的峰值作为信号的有效估计，即

其中，arg_kmax[·]表示沿着频率取最大值算子。

12)、对有效信号的估计值进行反尺度变换，恢复信号幅度。

13)、舍去滤波后结果的拓展边缘得到消除随机噪声后的信号分量s₁(n)。

14)、接下来，按照上述步骤7)到步骤13)依次对信号主导的模态分量c₂(n)和c₃(n)进行处理，得到消除随机噪声后的信号分量，如图7所示。

15)、对处理后得到的信号分量进行相加，得到最终消除随机噪声干扰的全波磁共振信号，如图8所示。