基于局部均值分解的地面核磁共振信号包络提取方法

田宝凤，孙士聪，刘隆昌，蒋川东

1.地球信息探测仪器教育部重点实验室(吉林大学)，长春 130026 2.吉林大学仪器科学与电气工程学院，长春 130026

0 引言

地面核磁共振(surface nuclear magnetic resonance, SNMR)又称为磁共振测深(magnetic resonance sounding, MRS)，是一种先进的能够定量进行地下水测量的地球物理勘探方法。与其他的地球物理方法相比，SNMR是一种定量、直接的找水方法，通过对地下水中氢质子的响应来实现地下水探测，无需钻孔就能够反演解释地下水的含量、位置以及储水介质的孔隙大小等水文地质信息[1]。因此，近年来，该方法在水资源探测、评价以及地质灾害隐患监测预警等方面均获得了广泛应用[2-3]。但在实际的信号采集与包络提取过程中，噪声环境较为复杂，往往会淹没纳伏量级的SNMR信号，导致目标信号特征无法辨识，严重制约了SNMR信号包络的提取[4]。

针对SNMR有效信号检测的问题，通常的技术手段是首先进行噪声压制，然后再通过数据拟合方法提取包络。普遍采用尖峰建模方法去除尖峰噪声[5]，采用基于最小二乘法理论的工频谐波建模等方法消除工频谐波噪声[6-7]，采用双奇异值分解等方法同时抑制工频谐波和随机噪声干扰[8]，实现参数提取。此外，田宝凤[9]等提出了以独立成分分析为基础的快速固定点算法实现信噪分离，王琦等[10]对随机噪声背景下多弛豫SNMR信号提出了基于稀疏表示的参数提取方法。尽管如此，由于噪声环境的高复杂性、不确定性与随机性，SNMR信号的包络提取问题依然是地面磁共振探测领域的前沿核心技术难题[11]。

因此，根据信号特征直接获得SNMR信号的包络提取方法吸引了诸多学者的目光。Ghanati等[12]基于经验模态分解(empirical mode decomposition，EMD)以及完全集合经验模态分解方法对有效SNMR信号衰减趋势进行提取，但只有信噪比较高时才能获得较好的包络提取效果。Liu Lichao等[13]通过对SNMR信号使用谱分析和滑动窗口的方式提取信号包络。上述方法虽然在一定前提条件下获得了一定效果，但依然存在一定的局限性。

为实现SNMR信号与噪声的有效分离，进行多种工况下信号包络的高精度提取，本文提出了基于局部均值分解(local mean decomposition, LMD)的SNMR信号包络提取方法[14]。LMD算法能够将信号分解成一系列乘积(product function，PF)分量，PF分量由包络信号与纯调频信号相乘获得，因而瞬时频率是正的、连续的、具有物理意义的[15]。相较于常用的EMD算法，LMD算法采用滑动平均法来处理目标信号，避免了EMD算法采用三次样条函数形成上下包络时产生的过包络和欠包络现象[16]。此外，LMD算法比EMD算法的端点效应小且扩散较慢[17]。本文首先采用LMD算法将原始含噪SNMR信号逐级分离，得到一系列的PF分量；其次通过PF分量逐一判定去除其他高频成分；然后通过余量提取SNMR信号的包络信息，进一步采用非线性拟合方法提取关键参数，进行误差分析[18]；再通过仿真实验，论证方法的提取效果与可行性；最后，通过野外实测数据进一步证明所提LMD算法的有效性和实用性。

1 方法原理

1.1 局部均值分解

LMD算法能够自适应地将非平稳、非线性的信号分解为一系列包含信号特征的PF分量，凸显信号不同时间尺度的波动特征，保留原始信号的幅频特性。

首先，对目标信号f(t)进行处理，得到第一个PF分量Z1(t)；接着，从原始信号f(t)中分离出Z1(t)，得到余量q1(t)，并将其作为新的待分解信号，循环分解；然后直至经过k次循环后，分离出全部的PF分量，得到残余分量qk(t)为一个单调函数，LMD算法分解结束。

(1)

最终，原始信号f(t)分解得到k个PF分量和残余分量：

(2)

式中，Zi(t)为信号的第i个PF分量。

1.2 PF分量计算

PF分量是包络信号与纯调频函数的乘积，这些乘积函数构成的有限集合中元素的总和就是复杂信号的时频分布，每一个元素称为一个PF分量。PF分量的本质是一个幅度、频率调制信号，它的瞬时频率可以用纯频率调制信号直接计算得出，反映了每个波动特征对应的物理意义。而包络信号就是PF分量的幅度调制函数。

对于频率调制函数和幅度调制函数的计算，则通过以下过程实现。

首先，需要依次比较信号中某一采样点与跟它相邻的两个采样点的大小关系，得到信号包络局部极值点ni，计算相邻两个极值点(包含极大值点和极小值点)的算术平均值mi为

(3)

计算相邻两个极值点的局部包络函数值ai为

(4)

再连接mi序列得到一条折线，对ai采用相同方式得到另一条折线。

随后，采用滑动平均法分别对2条折线进行平滑处理，计算出信号Z1(t)第一次循环时的局部均值函数和包络估计函数，记为m11(t)与a11(t)。从信号f(t)中分离出m11(t)得到余量h11(t)为

h11(t)=f(t)-m11(t)。

(5)

利用a11(t)对h11(t)进行解调，计算得频率调制信号为

(6)

式中：s11(t)为信号Z1(t)第一次循环下的频率调制信号。

接着，以s11(t)作为新的待分解信号，进行多次循环分解。当分解了n次后，得到n个频率调制函数与n个包络估计函数。Z1(t)的分解计算过程如下：

(7)

(8)

当Z1(t)的第n个频率调制信号s1n(t)的包络估计函数a1n+1(t)满足式(9)条件时停止迭代，此时，s1n(t)为Z1(t)的纯调频信号。

1-Δ≤a1n+1(t)≤1+Δ。

(9)

(10)

Z1(t)=s1n(t)a1(t)。

(11)

最后，从信号中分离出得到的PF分量，对剩余信号重复上述计算过程，循环迭代，最终得到k个PF分量和单调的残余分量qk(t)。

1.3 SNMR信号包络提取流程

在分解SNMR信号的过程中，信号频率与环境中掺杂噪声的频率往往十分接近。为此，首先对采集到SNMR信号采用Hilbert变换，转换成复包络信号，提取信号的实部与虚部；再对信号的实部与虚部分别进行局部均值分解处理，得到各自的PF分量与余量。分解后信号f(t)可以表示为

(12)

式中：Re[·]表示实部数据；Im[·]表示虚部数据。

根据SNMR信号频率的特征，分离出原始信号中的高频成分，可以得到SNMR信号的包络曲线：

x(t)=Re[qk(t)]+Im[qk(t)]·i。

(13)

对得到的包络曲线进行非线性拟合提取特征参数，以此还原出理想的SNMR信号。

我国国库集中支付制度起步较晚，从2001年开始确立，紧接着在几所本科高校实行改革试点，并在几年内在全国高校进行推广。发展到今天，国库集中支付制度已经得到了很多高校的认可，各项工作的开展正有条不紊的进行，对高校财务管理系统产生了较大影响。笔者根据实际工作和学习情况，对国库集中支付对高校财务管理体系的构建进行研究，希望促进我国高校财务管理体系的健康发展。

2 仿真实验与结果分析

针对采集SNMR信号时噪声环境较为复杂、随机性高的问题，本文基于LMD算法对仿真SNMR信号进行包络提取，并量化提取结果的信噪比、初始振幅拟合误差和平均横向弛豫时间拟合误差，来评估方法的有效性与优越性。信噪比RSN计算方式为

(14)

式中：fE为理想SNMR信号；fx为含噪信号；PE为理想SNMR信号的功率；PN为噪声的功率。

对提取到的理想SNMR信号趋势项进行非线性拟合，得特征参数的初始振幅和平均横向弛豫时间的拟合误差分别为：

(15)

(16)

式中：ΔE0为初始振幅拟合误差；E0′为所提取包络结果的初始振幅；E0为理想SNMR信号的初始振幅；ΔT2*为平均横向弛豫时间拟合误差；T2*′为所提取包络结果的平均横向弛豫时间；T2*为平均横向弛豫时间。

2.1 迭代精度设置

在对信号进行局部均值分解的过程中，内层循环通过解调源数据与局部均值函数之差得到调频信号，分析调频信号的包络是否为1来判断该信号是否为纯调频信号，进而决定是否跳出此次循环进行下一个PF分量的求取。但在实际计算过程中，为了避免迭代次数过多，造成运算量增大、程序运行时间过长的问题，采取式(9)进行判断，即纯调频信号包络信号的值域区间为[1-Δ，1+Δ]。Δ取值为10-7。

同时，在外层循环中，还需要判断原始含噪信号分离出所有PF分量后的余量qk(t)是否为一单调函数，决定是否继续迭代分解下一个PF分量。当qk(t)的极值点数目小于2时，可认为分解过程结束，原始含噪信号表示为所有PF分量与余量之和。

2.2 基于LMD的SNMR信号包络提取结果

理想的全波核磁共振信号表达式为

(17)

式中：fL为采集地点的拉莫尔频率；φ0为初始相位。

利用本文方法对目标信号进行处理。首先，对仿真含噪信号进行Hilbert变换，将得到的复包络信号分解为实部与虚部信号(图2)。然后，采用LMD算法，分别对实部信号和虚部信号进行分解，分别得到其各阶PF分量(图3)。以实部和虚部信号的最后2个PF分量作为特征分量合成目标信号包络(图4)。在前0.05 s内由于端点效应导致包络信号失真，但是在后续时间里包络信号与理想SNMR信号保持一致的趋势，包络信号在理想包络线处上下波动，经过算法处理后提取出的信号损失与失真程度很小。

a. 幅值-时间图；b. 幅值-频率图。

图2 含噪SNMR包络信号实部与虚部

对LMD算法处理结果采用最小二乘法进行非线性拟合，目标函数公式为

(18)

式中，fobs为 LMD算法处理后提取的包络信号。图1中也给出了LMD算法提取信号的拟合结果。由图1可见，拟合信号SNMR呈现明显的e指数衰减趋势，与理想SNMR信号重合程度较高，能量损失较小。局部均值分解后提取的SNMR信号初始振幅拟合误差为1.52%，平均横向弛豫时间拟合误差为1.93%，信噪比提升了29.0 dB。

为了进一步讨论LMD算法的适用范围，分析算法在不同环境条件下的处理效果，开展了不同工频谐波幅值、不同噪声幅度、不同平均横向弛豫时间的算法实验，结果如表1、表2、表3所示。用LMD算法处理SNMR信号，当信号初始信噪比大于-18.0 dB时，初始振幅和平均横向弛豫时间的拟合误差分别控制在±4.17%和±5.63%范围内，信噪比提升30.3～37.2 dB;当初始信噪比低于-18.0 dB时，参数拟合误差偏大，如当初始信噪比为-21.4 dB时，初始振幅的拟合误差可达15.02%，平均横向弛豫时间的拟合误差达到17.41%。

2.3 方法对比

为了更加全面、具象地说明算法的处理效果，通过仿真实验将LMD算法与经典的EMD算法和工频谐波建模方法进行对比，结果如图5所示。从图5中可以看出，对于3种方法的处理结果，LMD算法提取的SNMR信号包络可以更好地逼近理想SNMR信号，而且受到端点效应和模态混叠的影响更小。

经过3种方法处理前所得后的参数提取误差如

a. 实部；b. 虚部。

图4 LMD算法提取SNMR信号包络

表4所示，相比之下，LMD算法的初始振幅拟合误差、平均横向弛豫时间拟合误差均为最小，信噪比提升程度最大，证明了LMD算法的处理效果更加准确，提取的包络与理想SNMR信号包络趋势最为接近。

3 实测数据处理结果

为验证LMD算法对于实际复杂噪声环境中SNMR信号包络的提取效果，于长春市文化广场实地采集信号源产生的SNMR信号和环境噪声的合成数据开展算法实测数据实验验证。采用吉林大学自主研发全波磁共振探测仪器系统采集了4组不同参数的实测数据，采集线圈尺寸为25 m×25 m，采样率为25 kHz，采样时长为0.5 s，数据1、数据2与数据4的拉莫尔频率为2 355 Hz，数据3的拉莫尔频率为2 360 Hz。从图6中可以看出，含噪信号中都存在较大工频和随机噪声，且在拉莫尔频率点附近有噪声的密集分布。经过LMD算法处理后，由拟合后SNMR信号与含噪SNMR信号的对比可以发现，本文方法可以有效去除噪声，提取目标信号包

表1 不同工频谐波幅值下LMD算法处理结果的信噪比与参数拟合误差

表2 不同随机白噪声幅值下LMD算法处理结果的信噪比与参数拟合误差

表3 不同平均横向弛豫时间下LMD算法处理结果的信噪比与参数拟合误差

图5 LMD算法、EMD算法和工频谐波建模方法处理结果对比

表4 LMD、EMD与工频谐波建模参数拟合误差与信噪比提升程度

络。4组数据处理后平均横向弛豫时间拟合误差分别为1.53%、4.80%、3.65%、2.30%，可见，基于LMD的SNMR信号包络提取方法有着很好应用性与精确性，可以有效地从复杂噪声中提取目标的SNMR信号，反映水体特征。

a. 数据1；b. 数据2；c. 数据3；d. 数据4。

4 结论

1) 基于LMD算法实现了对不同信噪比下含噪信号复包络的实部与虚部的分解，基于余量组合实现了对SNMR信号包络的有效提取。算法仿真结果表明，该方法能够有效去除随机噪声、工频谐波等环境噪声，初始振幅的拟合误差在±4.17%之内，平均横向弛豫时间的拟合误差在±5.63%之内，信噪比提高了30.3～37.2 dB。

2) 对比分析了LMD算法与经典EMD包络提取算法和谐波建模方法的性能，表明LMD算法有着较好的优越性与准确程度，提取的包络与理想信号包络更加吻合。实测数据处理结果进一步验证了方法的有效性。该方法可以显著提升参数提取精度，降低运算时间，减少环境噪声对勘查结果的影响。 对于拓展磁共振探测仪器的应用范围具有借鉴意义。