基于总体最小二乘-旋转不变算法的地表核磁共振信号参数估计

于晓辉 冯 海 田宝凤 孙海欣 孙晓东*

①(吉林大学通信工程学院 长春 130022)

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

③(长春大学电子信息学院 长春 130022)

1 引言

地表核磁共振(Surface Nuclear Magnetic Resonance, SNMR)技术是目前唯一可以直接估算地表100 m下含水量和孔隙结构的方法[1]，在过去的几十年里，随着测量技术和数据处理技术的发展和成熟，SNMR技术广泛应用于地下水资源调查与水文地质评估等领域[2–4]。应用地表磁共振探测技术探测地下储水层时，其探测对象是目标水层中的氢质子。地下水中的氢质子受地球的磁场影响时会以拉莫尔频率绕静态磁场进动，当受到拉莫尔频率的激励脉冲时，氢质子会转为激发态，由激发态再次回到平衡状态时会发出一个衰减信号，即SNMR信号，通过采集SNMR信号并求得其参数即可获取该地区地下水的有关信息[5,6]。

表1 SNMR信号参数实际值以及使用TLS-ESPRIT算法求得的参数估计值

核磁找水仪器在现场采集到的SNMR信号主要包含3类噪声：随机噪声是一种在时间上随机产生的噪声，任意时刻其幅值、波形、相位都无法预测，但总体服从一定的统计分布规律。尖峰噪声由持续时间短、幅值大、能量高的不规则脉冲组成。工频谐波噪声是一种频率为基频的整数倍且持续时间长的噪声，主要来源为附近的电力系统、发电机组等工业及生活用电设[7]。

基于这3种噪声的特点与性质，目前针对SNMR信号抑制噪声方法的研究有很多成果。随机噪声常用的抑制方法是采用多次测量并进行叠加消除的方法，叠加次数越多效果越好。尖峰噪声抑制方法主要有：Legchenko等人[8]采用的加权叠加算法；Costabel等人[9]采用的小波变换；万玲等人[10]采用的基于能量运算的磁共振信号尖峰噪声抑制方法。工频谐波噪声是对SNMR信号影响最大的噪声，工频谐波噪声的信号特征与SNMR信号相近，因此抑制工频谐波噪声难度较大，最为常见的去除方法是谐波建模法，通过建立工频谐波噪声模型以实现抑制工频谐波干扰的目的。针对谐波建模法有很多改进，如Legchenko等人[11]采用的正弦减法；Larsen等人[12]采用基于模型的噪声抵消；田宝凤等人[13]采用基于谐波建模和自相关的去噪方法。

上述传统的方法以抑制或者消除信号中的噪声为目的，对得到的高信噪比的信号进行拟合以求得参数，参数的提取精度易受去噪效果的影响。3种噪声中，工频谐波噪声和随机噪声对信号的影响更大，是目前研究的重点。本文采用一种基于总体最小二乘-旋转不变法(Total Least Squares-Estimation of Signal Parameters via Rotational Invariance Technique, TLS-ESPRIT)的地表核磁共振信号参数提取方法[7]，可以从混有随机噪声和工频谐波噪声的信号中直接获取SNMR信号的参数。该方法利用了SNMR信号与工频谐波噪声相似的信号特征，以此构造出混合信号模型，原始信号变为随机噪声和混合信号之和，使用TLS-ESPRIT算法将混合信号参数提取问题转换为旋转不变矩阵的广义特征值求解[14–16]，以此计算出SNMR信号的弛豫时间和拉莫尔频率，并结合最小二乘法求得SNMR信号的初始振幅和相位。在使用TLS-ESPRIT算法前，引入Hankel矩阵构建子阵，子阵间的关系即对应着旋转不变矩阵[17,18]。通过仿真和实测数据实验，验证了本文基于TLS-ESPRIT算法的地表核磁共振信号参数估计方法在谐波噪声和随机噪声背景下的有效性，与谐波建模法进行的对比实验证明了其在参数提取精度上具有优势。

2 混合噪声模型

核磁共振找水仪工作时，采集到的SNMR信号常混有工频谐波噪声和随机噪声，采集信号可用式(1)所示的数学模型表示。

3 基于Hankel矩阵的TLS-ESPRIT算法

3.1 旋转不变矩阵

ESPRIT算法利用接收数据协方差矩阵的信号子空间的旋转不变性获取信号的参数，其核心是旋转不变矩阵的求取。

构造N维向量xn和ωn

Φ为旋转不变矩阵，其中包含SNMR信号的拉莫尔频率fL和弛豫时间T2∗，求取Φ的特征值后，即可估计出二者的值，初始幅值E0和初始相位φ0在获得前两个参数的估计值后可采用最小二乘法求得。

3.2 改进TLS-ESPRIT算法

在使用TLS-ESPRIT算法对采集信号进行数据处理时引入Hankel矩阵，将采集到的数据序列x(0),x(1),...,x(L) 构造两个M×N维Hankel矩阵如式(11)和式(12)所示

式(11)和式(12)中M,N与L的关系为

将式(11)所示的X1和式(12)所示的X2合并成一个新矩阵如式(14)所示

式中US为信号子空间，对应着混合信号，即SNMR信号和谐波噪声，UN为噪声子空间，对应着随机噪声，二者根据特征值的大小划分。协方差矩阵的特征值代表了功率谱密度，高斯白噪声的平均能量接近零，而SNMR信号和谐波噪声的能量远大于零，因此信号子空间对应的特征值大于噪声子空间对应的特征值，信号子空间和噪声子空间存在一个边界[20]。

通过式(20)、式(21)、式(24)、式(25)求得了混合信号的每个信号分量的频率、弛豫时间、振幅与相位，SNMR信号的4个参数也包含在其中。在混合信号模型中，工频谐波噪声事实上不存在弛豫时间参数，但可认为其弛豫时间为无穷大，因此Tl(l=1,2,...,K)中绝对值最小的值即为SNMR信号的弛豫时间T2∗，fl(l=1,2,...,K)中与之配对的为SNMR信号的拉莫尔频率fL，则使用最小二乘法求得的E和θ中与T2∗和fL对应的参数为SNMR信号初始振幅E0和初始相位φ0。

本文通过TLS-ESPRIT算法估计SNMR信号的拉莫尔频率和弛豫时间，然后使用最小二乘法求得SNMR信号的初始振幅和初始相位，参数估计的总流程如图1所示。

图1 基于TLS-ESPRIT算法的SNMR信号参数估计的流程图

4 仿真实验

为了分析TLS-ESPRIT算法对SNMR信号的参数估计效果，进行了如下仿真实验。(1) 分别改变随机噪声和工频谐波噪声的强度，在不同的噪声背景下使用TLS-ESPRIT算法估计SNMR信号的参数；(2) 改变SNMR信号的弛豫时间，对比TLS-ESPRIT算法对不同弛豫时间的SNMR信号的参数估计效果；(3) 将TLS-ESPRIT算法与谐波建模法进行对比实验，验证TLS-ESPRIT算法具有的优势。

4.1 不同噪声水平下的SNMR信号参数估计

本小节将分析不同噪声强度的随机噪声和工频谐波噪声对使用TLS-ESPRIT算法估计SNMR信号的影响，仿真SNMR信号的参数设为：E0=100 nV，T2∗=200 ms，fL=2 360 Hz，φ0= π/3，仿真信号添加的随机噪声为零均值，标准差σ为50或100的高斯白噪声，添加的工频谐波噪声幅值在(100 nV,200 nV, 300 nV, 400 nV)随机分布，基频为50 Hz，频率范围在2000～3 000 Hz，相位在[–π，π]rad范围内均匀分布。当SNMR信号处在标准差为100的高斯白噪声和幅值在400 nV分布的谐波噪声环境中时，SNMR信号与两种噪声之和的信噪比约为–18 dB。在8组不同噪声背景下分别进行100次蒙特卡罗实验，实验结果的分布如图2—图5所示。

图2 不同噪声水平下初始振幅统计

图3 不同噪声水平下弛豫时间统计

图4 不同噪声水平下拉莫尔频率统计

图5 不同噪声水平下初始相位统计

对图2—图5的SNMR参数提取统计结果进行分析，虽然每组SNMR信号处在不同的噪声背景下，但使用本算法估计得到的参数结果中位数与预设值基本一致。随着高斯白噪声标准差的增大，25%～75%范围框和最大最小范围变大，但工频谐波幅值增大时，25%～75%范围框和最大最小范围变化不明显。使用TLS-ESPRIT算法估计得到的SNMR信号两个参数中，拉莫尔频率的估计效果更好，弛豫时间的离散程度较大。使用最小二乘法求得的初始振幅和初始相位的精度会受到前两个参数的影响。

4.2 不同弛豫时间下的SNMR信号参数估计

本小节将分析SNMR信号的弛豫时间不同时，对使用TLS-ESPRIT算法估计其参数的影响，仿真建立4组弛豫时间为(0.2 s, 0.4 s, 0.6 s, 0.8 s)的SNMR信号，其他参数设为：E0=100 nV，fL=2 360 Hz，φ0= π/3，仿真信号添加的随机噪声为零均值、标准差50的高斯白噪声，添加的工频谐波噪声幅值在200 nV随机分布，基频为50 Hz，频率范围在2000～3 000 Hz，相位在[–π，π] rad范围内均匀分布。对4组不同弛豫时间的数据分别进行100次蒙特卡罗实验，提取结果的分布如图6—图9所示。

图6 不同 T2∗下初始振幅统计

图7 不同 T2∗下弛豫时间统计

图8 不同 T2∗下拉莫尔频率统计

图9 不同 T2∗下初始相位统计

图6—图9展示了弛豫时间参数改变对使用TLSESPRIT算法估计SNMR信号参数结果的统计。随着弛豫时间增大，初始振幅、拉莫尔频率和初始相位的25%～75%范围框和最大最小范围基本没有变化，而弛豫时间的提取结果的25%～75%范围框和最大最小范围明显增大，异常值出现。这是由于SNMR信号弛豫时间越大，衰减越慢，在相同时间内信号变化趋势不明显，因此弛豫时间估计效果随着弛豫时间的增大而降低。

4.3 与谐波建模方法的对比实验

谐波建模法是最常用的去除SNMR信号中工频谐波噪声的方法，它通过建立谐波噪声模型，将谐波噪声模型从含噪信号中减去而实现去除SNMR信号受到的工频谐波干扰，之后可通过非线性拟合估计SNMR信号参数。本小节将基于TLS-ESPRIT算法的SNMR参数估计方法与谐波建模法进行对比实验。对比实验的工频谐波噪声条件同4.1节实验一致，随机噪声取均值0、标准差为50的高斯白噪声，每组噪声背景下分别使用两种方法对SNMR信号进行参数估计，实验结果通过计算均方根误差(RMSE)衡量，如图10和图11所示。

图10 初始振幅RMSE

图11 弛豫时间RMSE

从图10和图11可以看出当谐波噪声的幅值增大，使用谐波建模法结合非线性拟合得到SNMR信号初始振幅E0和弛豫时间T2∗的精确度下降。但本文采用的TLS-ESPRIT算法受谐波噪声幅值变化的影响非常小，初始振幅E0和弛豫时间T2∗的RMSE在4种不同幅值的谐波噪声环境下均低于谐波建模法。在谐波噪声幅值较大时，使用TLS-ESPRIT算法对SNMR信号进行参数估计具有明显的优势。

5 实测数据分析

为了进一步验证使用T LS-ESPRIT算法提取SNMR信号参数的效果，使用信号发生器模拟地下水发出的SNMR信号，产生一组弛豫时间不同，振幅E0=400 nV，拉莫尔频率为fL=2 355 Hz的信号，并由吉林大学自行研制的核磁共振找水仪JLMRS-1进行采集。表1展示了SNMR信号参数实际值以及使用TLS-ESPRIT算法求得的参数估计值，图12展示了其中弛豫时间为0.2 s的采集信号和信号发生器生成的SNMR信号的时频域图，从图中可以看出信号完全淹没在噪声里，此时实测信号的信噪比约为–15 dB。

图12 实测信号和SNMR信号时频域图

实测数据实验结果表明，本文采用的TLS-ESPRIT算法能在实际噪声环境下有效估计核磁共振找水仪采集的信号参数，对于SNMR信号弛豫时间和拉莫尔频率的估计精度较高，初始幅值的误差相对较大，实测数据实验与仿真实验结果基本一致。

6 结论

针对SNMR信号在采集时易受到工频谐波噪声和随机噪声影响的特点，本文基于TLS-ESPRIT算法对SNMR信号进行参数提取，通过仿真对比实验和实测数据实验得出如下结论：

(1) TLS-ESPRIT算法能够从混有工频谐波噪声和随机噪声干扰的信号中估计出SNMR信号的关键参数，具有较好的效果。实测数据实验表明在实际噪声背景下此算法仍有效。

(2) 弛豫时间的改变会对使用TLS-ESPRIT算法估计SNMR信号参数造成影响，初始振幅、拉莫尔频率和初始相位受弛豫时间改变的影响较小，但随着弛豫时间的增大，弛豫时间自身的估计效果变差。

(3) 与谐波建模法相比，本文算法对SNMR信号的参数提取精度具有明显的优势，在工频谐波噪声强度较高时此优势尤为明显。