基于Bayes Bootstrap统计降噪方法的磁共振测深信号检测
张海如1，欧阳缮2，王国富2，张法全2
(1. 西安电子科技大学 电子工程学院，陕西 西安，710071；
2. 桂林电子科技大学 信息与通信学院，广西 桂林，541004)
摘要：为从受强噪声干扰的实测磁共振测深(MRS)信号中提取有效特征参数，针对导电层状模型MRS信号，利用各导电层MRS信号分量取对数后与采样时间具有严格线性关系的特性，提出从各导电层MRS信号分量中提取2个最优估计点来重建MRS信号的思想，即用2个最优估计点代替整体统计误差最小时的估计值，以提高重建精度。该方法先将多周期各导电层MRS信号分量序列转化为二维矩阵，再用Bayes Bootstrap模型对各列信息进行误差估计，最后，根据误差估计结果，统计出2个最优有效数据点。研究结果表明：本文方法提取MRS信号特征参数的平均误差率远比传统方法的小，能从信噪比为1 dB的强噪声环境中，以5.5%的平均误差率提取出有效特征参数；能从同时受到工频干扰、奇异噪声干扰和平稳随机噪声干扰的含噪MRS信号中提取有效特征参数。此外，本文方法能从较少周期的实测数据中提取出MRS信号的有效特征参数，可提高MRS仪器地质探测效率。
关键词：磁共振测深；统计降噪；信号检测；Bayes Bootstrap模型
中图分类号：P631 文献标志码：A 文章编号：1672-7207(2014)09-3144-06
Magnetic resonance sounding signal detection based on statistical noise reduction method of Bayes Bootstrap
ZHANG Hairu1, OUYANG Shan2, WANG Guofu2, ZHANG Faquan2
(1. School of Electronic Engineering, Xidian University, Xi’an 710071, China;
2. School of Information and Communication Engineering, Guilin University of Electronic Technology,
Guilin 541004, China)
Abstract: In order to extract a set of feature parameters of magnetic resonance sounding (MRS) signal, an algorithm of statistical noise reduction was presented over layered conductive earth, which took full use of the linear relationship aong the logarithmic MRS signal of each layered electrical conductive media and its sampling time. It estimated two optimal points from the real inspected signal sequence to reconstruct the MRS signal. Since it used only two optimal statistical points instead of the set of parameters with the minimum global statistical error, the algorithm could enhance the reconstruction accuracy of the MRS signal. The algorithm translated the multi-period MRS signal of each layered electrical conductive media sequence into two-dimensional matrix at first. Then, it structured the Bayes Bootstrap model to estimate the error of each column vector. Finally, according to these results of error analysis, the two optimal points were acquired. The results show the algorithm has a much lower average error rate than the traditional ones. It can extract effective feature parameters from the acquired MSR signal under poor SNR (1 dB) at the average error rate 5.5%. Furthermore, it also works for complicated original MSR signal, which contains the stable stochastic noise, the power frequency disturbing signal and the spike noise. With it the time for the sample period can be reduced greatly, inversion precision of the MRS systems can be improved.
Key words: magnetic resonance sounding; statistical noise reduction; signal detection; Bayes Bootstrap model
磁共振测深(MRS)技术具有操作简单、信息量丰富和唯一解等特性，而且是一种直接找水的物探方 法，随着对该项技术研究的逐渐深入和成熟，已在探测地下水、考古、地下水污染检测等领域得到了一定应用[1-4]。MRS技术通过向地下发射不同大小的脉冲矩，以激发不同深度范围内水中氢核，使其发生能量跃迁；当发射机停止发射信号时，已发生能量跃迁的氢核将恢复自由状态，接收天线能感应到该过程氢核释放的能量，接收机采集接收天线感应的信号，通过对接收信号反演解释，以判定探测范围内的水文地质构造[5-7]。但是，MRS信号为纳伏级信号，而实际探测环境复杂，MRS信号很容易被自然噪声和人为噪声淹没，这严重影响了MRS技术的应用范围。为提高MRS技术的实用性，近年来该技术得到了国内外专家学者的广泛关注和研究[8-10]，并且取得了一定成果。其中，MRS接收机采用了正交矢量型锁定放大技术，其探测精度得到了很大提高；设计特殊形状的天线能对噪声起到一定的抑制作用，但限制了设备的探测深度。此外，随着正、反演技术的逐渐发展使得MRS信号一维正、反演理论趋于成熟，研究者正在深入研究二维、三维正反演模型[11-14]。然而，目前关于MRS信号预处理方面的研究还较少，对强噪声背景下采集的MRS信号，通常采用多周期叠加的方式来提高信噪比，这大大增加了信号采集时间，而且只能在一定程度上削弱平稳噪声。为消除工频干扰，田宝凤等[15]提出基于参考线圈和变步长自适应的MRS信号噪声压制方法。该方法在理论上能较好地消除工频干扰，但是实用中对参考线圈的设计要求较高，影响了该技术的实用性。针对奇异噪声干扰，王中兴等[16]提出差阈值替代叠加方法。该方法与多周期叠加法类似，与其相比只是减弱了奇异噪声干扰，并没有消除其干扰。磁共振测深精确反演取决于是否能从接收到的MRS信号中提取有效特征参数[17]，而各导电层MRS信号分量随着采样时间的增加成指数规律衰减，对各导电层MRS信号分量取对数后，与采样时间具有线性关系，这样将MRS信号特征参数提取问题转化为求直线的斜率和截距的问题，因此，只需要2个有效点便能求出目标信息，而对MRS信号的多个周期采样包含了较多的信息冗余。为从采样得到的多周期各导电层MRS信号分量中估计出2个有效点信息，本文提出一种适合于小样本的统计降噪算法来提取MRS信号的特征参数信息，将多周期各导电层MRS信号分量序列转化为二维矩阵，矩阵各列为同一目标点的多次观测值，剔除各列向量中的野点，之后用Bayes Bootstrap方法估计各列信息的观测误差，最后，从估计结果中选出2个估计误差最小的列向量来估计有效数据点，以求取MRS信号特征参数。由于Bayes Bootstrap模型是小样本分析方法[18-19]，而且只需选出2个有效点信息，所以该方法大大减少了MRS信号的采样周期数，提高了MRS仪器地质探测效率；同时，该方法用最优估计值代替整体平均误差值，提高了MRS信号特征参数的提取精度。
1 实测各导电层MRS信号分量模型
MRS接收机采用正交矢量型锁定放大技术，其接收到的MRS信号模型为
(1)
其中：S(t)为实测信号；E(t)为实际目标信号；Npow(t)为工频噪声；Npul(t)为奇异噪声；N(t)为其他环境随机噪声。
接收到的MRS信号各个周期夹杂的工频噪声干扰信号的相位是随机的，即
(2)
其中：Apow为工频谐波信号的幅值；fpow为工频谐波信号的频率；
为工频谐波信号的随机相位。
由于地磁场强度从赤道附近的30 μT到南北极的60 μT，对应的拉莫尔频率为1.278 kHz到2.556 kHz，接收到的MRS信号要通过带通滤波器进行滤波处理，以去除目标频带外的干扰信号。经过滤波后的奇异噪声转化为半正弦干扰噪声，即
(3)
其中：Ap为幅值；Tw为干扰脉冲宽度。通常奇异噪声宽度Tw≤25 ms。
对于导电层状模型，将有效探测深度划分为K层，则相应需要发射K个不同的激发脉冲矩。随着激发脉冲矩的增大能依次激发到不同深度层的水中氢离子。当接收机采集M个周期的MRS信号，每个周期采样N个数据点时，得到接收数据向量
(k=1, 2, …,K)，将S(qk)转化为二维数据矩阵
(4)
令发射机关断时刻为0，每个采样点相对于发射机关断时间是恒定的，即
的所有行第j列数据
采样时刻相同，N个采样点的采样时间为
，令各导电层MRS信号分量为
，与对应的实际目标信号数据矩阵为
(5)
各列数据相等，则：
(6)
T与
(i=1, 2, …, M)满足的关系式为
(7)
计算式(7)中的2个未知参数
和
，只需要已知2个不同点坐标
和
。为此，只需要从中的N列数据中选出2列最有效的数据列，再分别对其所在列中的数据选择相对最有效的部分数据，估计出对应的和即可。
2 Bayes Bootstrap统计降噪算法用于MRS信号检测
2.1 各导电层MRS信号分量的Bayes Bootstrap估计
在矩阵中，
可等价为对
的M个观测记录，由于受到随机噪声的干扰，可以等价为离散型随机变量的1个样本空间，该样本空间的理论均值为，方差为0。令二值数据矩阵
表示中数据状态，
的取值为0或者1，1表示对应位置数据有效，0表示对应位置数据无效，且已被删除，将中全部元素初始化为1。
步骤1：剔除数据序列中的不稳定点。
1) 求中M个数的均值
和方差
。
2) 去除中与距离最远的数据点
，剩下的数据继续求均值
和方差
。若
＞
，则
，
，
，执行2)；否则，执行3)。通常，
。
3) 对中N列数据分别执行①和②得到中数据状态矩阵。
步骤2：用Bayes Bootstrap方法估计各列数据的误差。
1) 将中有效数据存入X中，即
。X的样本均值为
，用
估计，估计误差为
。
2) 对X中数据进行排序得
，对排序后的数据构造经验分布函数：
(8)
3) 从
中随机选取m个子样本
，其平均值为
。
4)设计统计量
，其中，
是满足Dirichlet分布的随机变量，即~
。
5) 用Dm分布来估计Tm分布，重复执行n次步骤3)和和4)，得到n个Dm，对其取平均值即为该组数据的估计误差，估计误差用
表示，i表示所在列序列。
6) 对中N列数据分别执行1)~5)。
2.2 各导电层MRS信号分量的特征参数估计
令
为
的估计值；
为各列数据之和。
(9)
找到
中数值最小的2个数据所对应的列记为l1和l2，则
和
即为接收信号最有效的2个估计值。将2点坐标
和
代入式(7)求取E0k和
2个对应的特征参数。
3 仿真验证和实验结果
为验证统计降噪算法的性能，进行统计降噪算法与自适应降噪[15]、差阈值替代叠加降噪[16]的对比实验，经过多次仿真实验，以特征参数E0和T2*的平均误差率来评价算法性能；同时，还以实测数据对3种方法进行实用性验证。信噪比fSNR公式为
(10)
其中：xi为信号序列；n i为噪声序列。E0和T2*的平均误差率fMER公式为
(11)
其中：Nnum为实验次数；Ei和
为从每次实验数据中提取到的特征参数值；E0和T2*为实际特征参数值。
3.1 仿真实验
在信噪比fSNR=1.5 dB，T2*=400 ms时，E0取不同值，每种取值实验1 000次，每次实验生成20个周期的单一导电层的含噪MRS信号，3种特征参数提取方法性能对比结果如表1所示。实验中，随机噪声强度为E0的4倍；工频干扰：幅值为1 nV，频率为2.2 kHz，相位随机产生；奇异噪声：脉冲强度为20~600 nV之间的随机数，脉宽为5~25 ms之间的随机数，起始时间为0~400 ms之间的随机数，单周期信号中脉冲个数为0~15之间的随机数。从表1可知：统计降噪算法的性能远优于自适应降噪和差阈值替代叠加降噪方法，证明在强噪声背景下，现有的特征参数提取方法无效时，统计降噪算法仍能较精确的提取特征参数。
在E0=350 nV，T2*=400 ms时，fSNR取不同值，每种取值实验1 000次，每次实验生成20个周期的单一导电层的含噪MRS信号，3种特征参数提取方法性能对比如表2所示。实验中，根据信噪比不同，调整随机噪声、工频干扰噪声和奇异噪声的参数。从表2可知：在信噪比较低时，差阈值替代叠加降噪方法略优于自适应降噪方法；在信噪比较高时，自适应降噪方法优于差阈值替代叠加降噪方法；统计降噪算法在不同信噪比下性能都优于自适应降噪和差阈值替代叠加降噪方法，而且性能较稳定，整体精度较高；当信噪比较高时，自适应降噪方法性能趋近于统计降噪算法，但仍比统计降噪算法精度低。
3.2 实测数据验证
为验证统计降噪算法的实用效果，在广西百色市进行野外实验，实验场地附近有高压线通过，且多云有雷电干扰。该探测区域附近的水文钻孔显示垂深25~40 m为第四系含水层，其含水质量分数为38%，孔隙度为42%。实验中，铺设100 m×100 m的线圈，发射频率为1 918 Hz，发射10个脉冲矩，分别为87，169，300，455，656，965，1 440，2 054，2 830和3 806 A·ms。当发射3 806 A·ms脉冲矩时，接收到的单周期MRS信号如图1所示。从图1可见：实验中接收到的MRS信号受到严重干扰。
表1 E0取不同值时3种特征参数提取方法性能对比表
Table 1 Performance comparisons about feature parameter extraction when E0 takes different values
表2 fSNR取不同值时3种特征参数提取方法性能对比表
Table 2 Performance comparisons about feature parameter extraction when fSNR takes different values
对不同发射脉冲矩下接收到的MRS信号，分别用自适应降噪、差阈值替代叠加降噪和统计降噪算法提取特征参数，之后进行一维反演，反演采用文献[20] 中已经成熟的反演算法，其中E0用于反演各地质层的含水质量分数，T2*用于反演各含水层的孔隙度。图2所示为反演得到的该测区含水量垂直分布图对比结果，其中各矩形条的长度对应含水量，颜色对应孔隙度。从图2可见：由于实测数据受噪声干扰严重，用自适应降噪和差阈值替代叠加降噪方法提取的特征参数进行反演，反演结果不能反映实际水文地质构造；用统计降噪算法提取的特征参数进行反演，反演结果与实际水文地质构造吻合，进一步验证了统计降噪算法的实用性。
图1 实测单周期MRS信号
Fig. 1 Monocyclic original MRS signal
图2 含水质量分数垂直分布图
Fig. 2 Vertical distributions of water mass fraction
4 结论
1) 充分挖掘导电层状模型下MRS信号的特性，构建了同时具有工频干扰、奇异噪声干扰和平稳随机噪声干扰的含噪MRS信号模型。
2) 利用取对数后的各导电层MRS信号分量与采样时间具有严格线性关系的特性，建立了从含噪各导电层MRS信号分量中提取2个最优估计点来重建MRS信号的模型。将一维多周期各导电层MRS信号分量序列转化为二维矩阵，并用统计方法从二维矩阵中挖掘有效信息，以提高MRS信号的重建精度。该方法的实质是用2个最优估计值代替整体平均误差最小时的估计值，提高了MRS信号特征参数的提取精度。
3) 与现有方法相比，本文方法从采样信号中提取的信息量较少，而且用小样本误差估计模型Bayes Bootstrap来处理信号，只需较少的采样周期数便能从实测数据中提取出MRS信号的有效特征参数，提高了MRS仪器地质探测效率。
参考文献：
[1] Greben J M, Meyer R, Kimmie Z. The underground application of Magnetic Resonance Soundings[J]. Journal of Applied Geophysics, 2011, 75(2): 220-226.
[2] Legchenko A, Vouillamoz J M, Roy J. Application of the magnetic resonance sounding method to the investigation of aquifers in the presence of magnetic materials[J]. Geophysics, 2010, 75(6): 91.
[3] Nils P, Kamhaeng W, Khamphouth P, et al. Characterization of aquifers in the Vientiane basin, Laos, using magnetic resonance sounding and vertical electrical sounding[J]. Journal of Applied Geophysics, 2011, 73: 207-220.
[4] Chalikakis K, Nielsen M R, Legchenko A, et al. Investigation of sedimentary aquifers in Denmark using the magnetic resonance sounding method (MRS)[J]. Comptes Rendus Geoscience, 2009, 341(10): 918-927.
[5] Hertrich M, Yaramanci U. Joint inversion of surface nuclear magnetic resonance and vertical electrical sounding[J]. Journal of Applied Geophysics, 2002, 50(1): 179-191.
[6] Yaramanci U, Legchenko A, Roy J. Magnetic resonance sounding[J]. Journal of Applied Geophysics, 2008, 66(3): 71-72.
[7] Legchenko A, Valla P. A review of the basic principles for proton magnetic resonance sounding measurements[J]. Journal of Applied Geophysics, 2002, 50(1/2): 3-17.
[8] 翁爱华, 李舟波, 王雪秋. 层状导电介质中地面核磁共振响应特征理论研究[J]. 地球物理学报, 2004, 47(1): 156-163.
WENG Aihua, LI Zhoubo, WANG Xueqiu. A study on surface nuclear magnetic resonance over layered conductive earth[J]. Chinese Journal of Geophysics, 2004, 47(1): 156-163.
[9] Lubczynski M, Roy J. Hydrogeological interpretation and potential of the new magnetic resonance sounding (MRS) method[J]. Journal of Hydrology-Amsterdam, 2003, 28(1): 19-40.
[10] Vouillamoz J M, Favreau G, Massuel S, et al. Contribution of magnetic resonance sounding to aquifer characterization and recharge estimate in semiarid Niger[J]. Journal of Applied Geophysics, 2008, 64(3): 99-108.
[11] Weng A H. Occam’s inversion of magnetic resonance sounding on a layered electrically conductive earth[J]. Journal of Applied Geophysics, 2010, 70(1): 84-92.
[12] Hertrich M, Braun M, G
nther T, et al. Surface nuclear magnetic resonance tomography[J]. IEEE Transactions on Geoscience and Remote sensing, 2007, 45(11): 3752-3759.
[13] Braun M, Yaramanci U. Evaluation of the Influence of 2-D Electrical Resistivity on Magnetic Resonance Sounding[J]. Journal of Environmental & Engineering Geophysics, 2011, 16(3): 95-103.
[14] Legchenko A, Descloitres M, Vincent C, et al. Three- dimensional magnetic resonance imaging for groundwater[J]. New Journal of Physics, 2011, 13, 025022 (1-17).
[15] 田宝凤, 林君, 段清明, 等. 基于参考线圈和变步长自适应的磁共振信号噪声压制方法[J]. 地球物理学报, 2012, 55(7): 2462-2472.
TIAN Baofeng, LIN Jun, DUAN Qingming, et al. Variable step adaptive noise cancellation algorithm for magnetic resonance sounding signal with a reference coil[J]. Chinese Journal of Geophysics, 2012, 55(7): 2462-2472.
[16] 王中兴, 荣亮亮, 林君. 地面核磁共振找水信号中的奇异干扰抑制[J]. 吉林大学学报(工学版), 2009, 39(5): 1282-1287.
WANG Zhongxing, RONG Liangliang, LIN Jun. Spike noise elimination from surface nuclear magnetic resonance signal for underground water[J]. Journal of Jilin University (Engineering and Technology Edition), 2009, 39(5): 1282-1287.
[17] Legchenko A, Baltassat J M, Bobachev A, et al. Magnetic resonance sounding applied to aquifer characterization[J]. Ground Water, 2004, 42(3): 363-373.
[18] Abraham D A, Lyons A P. Reliable methods for estimating the K-distribution shape parameter[J]. IEEE Journal of Oceanic Engineering, 2010, 35(2): 288-302.
[19] Candy J V. Bootstrap particle filtering[J]. IEEE Signal Processing Magazine, 2007, 24(4): 73-85.
[20] Mohnke O, Yaramanci U. Smooth and block inversion of surface NMR amplitudes and decay times using simulated annealing[J]. Journal of Applied Geophysics. 2002, 50(1/2): 163-177.
(编辑 邓履翔)
收稿日期：2013-09-06；修回日期：2013-11-28
基金项目：国家自然科学基金资助项目(61102115，61371186，61362020)；广西自然科学基金资助项目(2012GXNSFBA053177, 2013GXNSFFA019004)；桂林电子科技大学广西信息科学实验中心经费资助(2013112)
通信作者：张海如(1987-)，男，天津宝坻人，博士研究生，从事磁共振测深信号建模、自适应信号处理、统计信号分析研究；电话：13557133936；E-mail: zhanghairu66@126.com
