基于模拟退火法的磁共振测深多源谐波噪声压制方法

陈亮，付立恒，蔡冻，李凡，李振宇，鲁恺

(1.中国电建集团 江西省电力设计院有限公司，江西 南昌 330096； 2.中国地质大学(武汉) 地球物理与空间信息学院，湖北 武汉 430074)

0 引言

1 MRS方法的基本原理及多源谐波电磁噪声特征

1.1 MRS方法的基本原理

MRS方法以NMR理论为基础，即地下水中的氢核磁矩受到一定频率的交变电流脉冲激发而发生偏转，磁矩弛豫过程中产生核磁共振信号，位于地面的接收线圈接收到一个自由感应衰减信号，即NMR信号，其过程如图1所示。图中：0～t1为仪器激发前的纯噪声采集时间；t1～t2为仪器发射拉莫尔频率的激发脉冲持续时间，脉冲包络线为矩形；t2～t3为发射转接收时间，称为死区时间；t3～t4为接收NMR信号时间，NMR信号的包络线按指数规律衰减。NMR信号ENMR为：

(1)

图1 MRS方法激发/接收NMR信号时序图Fig.1 Transmitting and receiving NMR signal in time domain of MRS

1.2 干扰噪声的种类及其特征

1.2.1 噪声种类

MRS方法接收到的NMR信号通常非常微弱，一般在10～1 000 nV之间，由于仪器接收灵敏度高易受到外部环境电磁噪声的干扰，故接收到的信号信噪比通常较低，提取NMR信号的难度较大[17]。按照噪声成分及来源，电磁噪声分为尖峰噪声、谐波噪声和随机噪声等，可以将接收到的NMR信号和电磁噪声表示为

Eobs=ENMR+Espike+Eharmonic+Erandom,

(2)

式中:Eobs为仪器接收到的信号数据,ENMR为磁共振信号,Espike为尖峰噪声,Eharmonic为谐波噪声,Erandom为随机噪声。

1.2.2 多源谐波电磁噪声特征

在经典谐波噪声模型条件下，2013年Larsen等[14]提出假设谐波噪声的基频在信号采集持续期间内是单一且稳定的，由此写出如下谐波模型式：

(3)

式中：Eharmonic(f0,n)表示采样时间n处的谐波模型的振幅，Ak、φk分别为第k次谐波分量的振幅和相位，fs为采样频率，f0为谐波噪声的基频。

随着MRS勘查工作逐步向城市靠近，面对的电磁噪声情况也愈发复杂，在测量区域可能存在2个或2个以上的谐波干扰源。本文分析存在2个谐波源情形下的谐波噪声特征，多源谐波与之相通[16]。

假设有2个基频不同的谐波干扰源，且在测量时间内不变，可以写出2个谐波源的谐波模型式：

(4)

式中：f1和f2分别表示要搜索的2个不同的基频，Ak和Bk分别表示每个基频的第k次谐波的振幅，φk和φk分别表示每个基频的第k次谐波的相位。

含有双基频的谐波噪声频谱相对于单一且稳定基频的谐波噪声频谱会出现2种特殊的峰值频谱特征：“弯曲”和“双峰”(图2)，“双峰”峰值频谱特征出现的条件较“弯曲”峰值频谱要更为严苛，要求频率间隔大且幅值相近，因此“双峰”峰值频谱的出现概率要低于“弯曲”峰值频谱。

2 同步删除法和模拟退火法原理

2.1 同步删除法原理

由双基频谐波噪声的特征可知[17]，主要问题是谐波的频率彼此接近但不完全相同，使得其频率内容重叠但不完全一致。为解决双基频谐波噪声的相互干扰问题，可以拓展模型去噪的搜索方式，在一个二维频率空间内同步搜索2个基频，求取观测信号与双基频谐波模型的最小值：

‖Eobs-Eharmonic(f1,f2,n)‖→min 。

(5)

由于式(4)中的f1、f2、Ak、Bk、φk和φk的确定是一个非线性优化问题，可以通过标准化改写式(4)中的余弦项：

图2 双基频的谐波噪声频谱特征示意Fig.2 Spectrum characteristics of harmonic noise with double fundamental frequencies

(6)

式中变量αk、βk、γk和δk可通过以下关系与对应谐波源的振幅和相位相关联：

(7)

(8)

将式(6)中的余弦项展开变换，假设f1和f2已知，观测信号与谐波模型的差值最小值问题就变成线性，可以用矩阵表示法写成：

(9)

式中：[k0,k1,…,kN]表示构建谐波模型的阶数，[n0,n1,…,nf]为采样时间样本。同时，可将式(9)整理成矩阵形式：

Ax=b，

(10)

式中：x=[αk0,βk0,γk0,δk0,…,αkN,βkN,γkN,δkN]T,b=[Eobs(n0),Eobs(n1),…,Eobs(nf-1),Eobs(nf)]T。利用最小二乘法求解该标准线性方程可得：

x=(ATA)-1ATb。

(11)

每次代入1对假设基频值f1和f2，可以得到1组谐波参数。利用式(5)寻找1组能使得该差值最小的假设的基频值作为构建谐波模型的基频值，再利用测量数据Eobs减去构建的Eharmonic，即可达到同步压制2个干扰源谐波噪声的效果。

2.2 模拟退火法原理

模拟退火法的发展受到金属物理退火的启发[18]：当金属从红热状态缓慢冷却时，金属内部形成有序的最小能量晶体结构。最初，当金属处于红热状态时，原子的运动完全由随机的热涨落控制，但是，随着温度的缓慢降低，原子间的作用力变得越来越重要，最后，原子变成了一个代表最小能量结构的晶格[19]。在模拟退火算法中，参数T代表模拟温度，误差E代表模拟能量，T值很大时算法的行为类似于Monte Carlo搜索，T值很小时算法的行为更有方向性。在物理退火中，一开始是一个大值的T，然后随着越来越多的试验被检验，T值慢慢降低；最初，大量的模型空间是随机采样的，但是随着搜索的进行，搜索变得越来越有方向性[20]。

模拟退火法从一个初始模型m(p)和其对应的误差E(m(p))开始。然后生成一个在m(p)附近的测试模型m(*)，计算其对应的误差E(m(*))，m(*)可以通过向m(p)增加一个从高斯分布中提取的增量Δm。当E(m(*))≤E(m(p))时，m(*)会替代m(p)成为新的比较值。但是,有时E(m(*))>E(m(p))，m(*)也可以接受。为了决定后面的情况，这个测试参数t可通过

(12)

计算。在区间[0,1]生成一个均匀分布的随机数r，如果t>r，则接受m(*)。当T值很大时，测试参数t趋近于1，故不管误差值是多少，m(*)几乎总是被接受的，这对应着“热运动”情况下，模型参数在空间上以无向的方式探索；当T值很小时，测试参数t趋近于0，m(*)几乎不被接受，这与定向搜索情况相对应，只有使误差E减小的测试模型才被接受。

3 仿真数据验证

3.1 网格搜索同步删除法

前文论述过，若要式(5)成立需要多次求解线性问题，同时优化二维网格搜索，每个维度对应不同的干扰源谐波基频。为了加快搜索速度，可以将搜索分2步进行：首先，覆盖相对宽且粗糙的网格，每次计算模型减去后的剩余信号；然后，从在第一步中获得的最小值开始探索更细的网格，并且通过第二步获得的最小剩余信号产生最佳频率估计值。

图3a是在一个相对粗糙的网格内搜索2个基频值，即图中2个相互对称的深蓝色色块;图3b是选取图3a中的一个色块范围进一步精细化搜索2个基频,可以明显看出找到的基频值与仿真设定的一致。搜索基频频率的精度要达到0.001 Hz需要进行600多次的搜索，处理时间为60.36±1.27 s。图4显示，同步删除法能够压制双基频谐波噪声，经过处理的信号与真实NMR信号相似程度较高，呈指数衰减趋势。图中纵坐标变量A表示振幅。

图3 频率空间谐波基频网格搜索Fig.3 Harmonic fundamental frequency grid search in frequency space

图4 网格搜索同步删除法处理双基频谐波时域Fig.4 Time domain diagram of double fundamental frequency harmonics processed by grid search simultaneous removal method

3.2 模拟退火同步删除法

为了加快基频的搜索和避免陷入局部极值点，利用模拟退火法来搜索干扰源谐波的基频值。仿真条件同上，迭代次数200次，搜索路径如图5所示。图5中黑色圆圈代表初始模型(50.000；50.000)，白色圆圈代表每次迭代被接受的模型，绿色圆圈代表迭代结束最终模型(50.031 7；49.984 6)，红色线段代表基频搜索路径。可以看出，每次被接受的模型都在逐步逼近真实模型，虽然模拟退火法最终能找到最小值，但在有限的迭代次数时搜索路径会存在一定的偏差。为了更好地评价去噪效果，特意引入均方根误差(root mean square error，RMSE)这个指标参数[21]：

(13)

式中:N为仪器采集到的信号长度，Edenoised(i)是去噪后的第i个数据，ENMR(i)是真实NMR信号的第i个数据。

图5 模拟退火同步删除法的基频搜索路径Fig.5 Fundamental frequency search path of simulated annealing simulated annealing

复杂的电磁噪声可能存在3个或更多基频的多源谐波噪声。为了验证模拟退火同步删除法能够适用于多源谐波电磁噪声的情况，选取了3个干扰源进行仿真，谐波基频分别为49.950、50.010、50.050，迭代次数300次，初始基频模型(50.000；50.000；50.000)，其他仿真条件同上，图6为仿真结果。从图6a可以看出模拟退火同步删除法在3个基频谐波的情况下能够较好地压制多源谐波噪声，去噪后的蓝色曲线与真实NMR信号的红色曲线基本吻合。图6b中黑色圆圈显示了迭代初始模型的频率空间位置，蓝色圆圈是300次迭代结束最终模型的频率空间位置(50.009 2；49.951 2；50.049 9)，红色圆圈是3个真实谐波基频不同排列顺序的频率空间位置，灰色曲线代表基频搜索路径；经过一定迭代次数的搜索，最终结果会归于6个红色圆圈中任意一个的附近。

表1 不同迭代次数的模拟退火同步删除法去噪效果

图6 模拟退火同步删除法压制多源谐波噪声Fig.6 Simulated annealing simultaneous removal method suppressing Multi- source harmonic noise

4 实测数据检验

为了验证本文提出的噪声压制方法在实测数据中的应用效果，选择天门市多宝镇汉江边的一处场地进行实验。该场地地下含水层已由场地内的钻孔所揭露，场地附近有多个人居村庄，符合多源谐波的噪声特点，实测数据可以很好地应用于去噪处理的实验工作中。数据采集使用法国IRIS公司的NUMISpoly仪器。本次测量采集技术参数：方形线圈边长100 m；拉莫尔频率fL=2 135.2 Hz；激发脉冲矩范围为0.12～10.866 A·s，分为16个；采样率19 200 Hz；采样时间1 s；死区时间0.04 s。

通过分析野外实际工作环境和噪声，发现2个较强的谐波干扰源。根据数据处理流程对噪声数据进行了尖峰噪声的处理后，采用100次迭代的模拟退火同步删除法来压制双谐波噪声(图7)，再将去噪后的信号进行叠加、希尔伯特变换和非线性拟合，结果见图8。

图7是利用模拟退火同步删除法来压制实测双谐波噪声。图7a中经过模拟退火同步删除法压制后的去噪曲线(红色曲线)有明显的指数衰减趋势，仅剩余随机噪声的干扰；图7b中相对于实测含噪信号(灰色曲线)，去噪曲线在50 Hz整数附近的峰值被较好压制，拉莫尔频率的信号保存较为完整。

图8中灰色曲线是去噪后的信号经过叠加处理和希尔伯特变换得到的包络曲线，红色曲线是灰色曲线经过非线性拟合得到的拟合曲线。包络曲线和拟合曲线重合性较好，表明模拟退火同步删除法能够较好地压制多谐波噪声，也能看出残余谐波噪声对包络曲线的影响在可接受的范围，且对拟合曲线的参数提取影响较小。

图7 模拟退火同步删除法在实测数据处理中的效果Fig.7 Effect of simulated annealing simultaneous removal method in measured data processing

图8 去噪后各脉冲矩的包络线及其拟合曲线Fig.8 Envelope and fitting curve of each pulse moment after denoising

图9 天门河滩MRS方法某测深点实测数据及其反演结果与钻探结果对比Fig.9 Measured data, inversion results and drilling results of a sounding point by MRS method in River Beach, Tianmen

5 结论

研究与实践表明，对于单一且规律性的谐波电磁噪声，利用模型去噪的方法能够有效地压制，而变化、多源的谐波噪声的压制是磁共振测深方法信号处理中的难点和研究热点，压制效果的好坏直接影响后续反演解释工作的准确性。根据多源谐波电磁噪声的特点出发，验证了网格搜索同步删除法的可行性，为了节约计算时间成本，提出了模拟退火同步删除法这一新的磁共振测深方法信号处理手段，通过仿真数据验证和实测数据检验，本文论述的信号处理手段取得了较好的噪声压制效果。

1)针对磁共振测深信号中模型去噪方法无法有效压制多源谐波噪声的问题，利用网格搜索同步删除法和模拟退火同步删除法均能有效找到较为准确的基频值，达到压制多源谐波的目的。

2)通过模拟仿真实验，证明网格搜索同步删除法和模拟退火同步删除法均能压制多谐波噪声。在双基频谐波情况下比较不同迭代次数下的压制效果，模拟退火同步删除法较网格搜索同步删除法的计算效率提高了2.35倍，大大降低了时间成本，且在多源谐波噪声压制方面也能有较好的效果，为开展更加复杂的谐波情况的压制提供了可能。

3)通过野外MRS方法的实测资料证明了模拟退火同步删除法压制多源谐波噪声的有效性，反演得到的含水量与钻探结果基本相符，极大地拓展了磁共振测深方法在人居区等电磁情况复杂地区作业的有效性，也为多通道磁共振多维探测提供新的数据处理思路。

仿真和实测结果表明，模拟退火同步删除法能够有效压制磁共振测深方法中的多源谐波噪声，为后期反演资料的解释提供了技术保证和选择方法。本文主要讨论单通道情况下多源谐波噪声的压制，未来将继续研究二维三维多通道变化谐波噪声的压制方法。