自旋回波在地面核磁共振找水技术中的反演算法初探

陈诗扬

（中国地质大学（武汉）地球物理与空间信息学院，武汉 430074）

自旋回波在地面核磁共振找水技术中的反演算法初探

陈诗扬

（中国地质大学（武汉）地球物理与空间信息学院，武汉 430074）

当测区磁场分布不均匀时，地面核磁共振技术所反映出的地下孔隙度会相对真实情况偏小。而将自旋回波信号应用于地面核磁共振方法中，则会有效克服磁场的不均匀分布对该方法的影响。因此，本文尝试提出一种利用自旋回波信号求取地面核磁共振方法中横向弛豫时间谱的反演算法，为今后将自旋回波技术应用于地面核磁共振方法的研究提供参考资料。本算法主要数学模型为多指数反演模型，反演算法为奇异值分解法。

地面核磁共振；自旋回波；多指数反演；奇异值分解

1 地面核磁共振基本原理

在静磁场B0的作用下，具有磁矩μ的原子核将受到力矩μ×B0的作用，沿磁场B0方向为轴（设为Z轴方向），以ω0=γ|B0|做拉莫尔（Larmer）进动，γ为磁旋比，其中μ与B0方向之间的夹角保持固定不变[1]。

现在垂直于B0方向上施加一个射频场B1，其振动角频率为ω0。在B1作用下，自旋磁矩μ绕射频磁场B1方向的进动，且与Z轴方向之间的夹角逐渐增大，这个夹角又叫扳转角，其大小与γ、B1和它持续作用的时间τ（也叫脉冲宽度）的乘积呈正相关。射频场B1的施加，使得原自旋系统磁化强度的纵向分量减小，同时磁化强度的水平分量增加。由于吸收了射频场B1的能量，原子由稳定的低态变成了不稳定的激发态。因此当射频场B1撤去后，原子会发生跃迁，向外辐射电磁波，而这个辐射出来的电磁波信号，即是我们通常所要测量的自由感应衰减信号E（t，q）[2]。

当交变电流脉冲结束作用后，由同一线圈接收天线接收到自由感应衰减信号按指数规律衰减。以下是地面核磁共振方法的基本方程：

其中，M0是氢质子的磁化强度；是激发场垂直于地磁场的分量；r为激发半径；I0为电流脉冲的振幅；α为地磁场倾角。根据式（1），可由以上参数与观测值E（t，q），求出表示地下含水量的与表示地下孔隙度的

2 自旋回波[3]信号介绍

当稳定磁场B0为不均匀磁场时，不同位置上自旋的同种原子核受到的磁场强度不同，因此它们以不同的进动频率绕着B0作进动，具有一个总磁矩M。通过用90°—τ—180°—2τ—180°—2τ—……的脉冲序列，使总磁矩M沿XZ平面多次翻转，其中τ为相邻脉冲之间的间隔。设进动频率较快的原子核的相位矢量为A，较慢的相位矢量为B，当向自选系统施以90°脉冲之后，A和B因为进动频率不同而散开。τ时间之后施以180°脉冲，自旋系统绕y轴翻转了180°，因此原来进动的快的A现在落到了B后面追赶B，又经过τ的间隔A与B重聚，重聚时便产生自旋回波信号。随着脉冲序列的不断进行，该自选系统还会以同样的形式产生一系列自旋回波信号，这一系列自旋回波信号称之为CPMG回波串。

在不均匀磁场中，由于不同原子核进动速率不同，使得自旋系统发生散相，在通过180°翻转脉冲后，分离的相位仍可在重聚后得到恢复。因此，我们可以通过一连串脉冲序列得到一组CPMG回波串，其T2*衰减曲线上每一个相位重聚时刻所对应的点就是没有因为磁场不均匀而产生散相的点，即真实T2衰减曲线上的点。因此，我们可以通过将这些零散的有用点拟合成曲线，得到真实的T2衰减曲线。

3 磁场梯度对地面核磁共振信号的影响

当存在磁场梯度时，地面核磁共振方法中的 并不等于真实值，因此无法反映出地下孔隙度的实际情况。与真实横向弛豫时间 的关系如式（2）所示[4]：

由于实际测区磁场分布通常具有不均匀性，一般很难通过（2）式对所测T2*值进行高精度校正。而自旋回波信号中的横向弛豫时间不会受到磁场不均匀性的影响，因此将自旋回波技术应用到地面核磁共振方法中，能够有效拓宽该方法的应用范围，例如可在火成岩发育地区进行找水工作。

4 自旋回波信号的T2解谱算法

对于复杂地质条件中的自旋回波串，其曲线不再是单一的一条简单的衰减曲线，而是多个指数衰减函数叠加而成的复杂曲线。因此，我们要从这个复杂的衰减曲线里提取出来的不是单个T2值，而是组成该曲线的每一个指数衰减函数所对应的T2值的集合：T2谱。

（1）反演模型

由于原始回波串是由多个指数衰减函数所叠加而成的，因此本文中选取的数学模型为多指数反演模型。又因为往往原始回波串中是存在噪声的，所以拟合的式子可以如下表达：

其中Mj（t）是tj时刻回波信号观测幅度，m为每个回波弛豫分量数，Xi为第i种弛豫分量零时刻的信号大小，noise为噪声[5]。

（2）T2时间解谱算法

假设观测到的回波数为n个，弛豫分量数为m，那么原拟合式子就可以写成M=AX的矩阵乘积形式[6-8]，如上所示。其中，M为所有回波观测数据的矩阵，X为所有弛豫分量的矩阵，A为系数矩阵。首先考虑AX=M的摄动问题。当M因为存在噪声而摄动时，必然要引起X的摄动，设M经过摄动后变为M+ΔM，解X经过摄动后变为X+ΔX，可以推知：

取T2i=2i+1，考虑到自扩散的影响，回波间隔不能太大，因此取Δt=0.5ms，最后求得cond（A）≈105，由（4）式可知当M因为存在噪声而存在一丝毫摄动时，都会引起X产生巨大的变化。因此，这是一个病态问题。为了解决这个病态问题[9]，现利用奇异值分解法求解ΔX的最小二乘解[7]，这里的ΔX指的是真实解X与假定的已知解之差，并在求解时对S矩阵求逆的过程中引入阻尼系数f来使解稳定。引入阻尼系数f将原本对S矩阵求逆的式子Si（k，k）=1/S（k，k）变为Si（k，k）=S（k，k）/（S（k，k）^2+f^2），Si是所求的近似逆矩阵。

以下是解谱算法的主要步骤：

（1）对矩阵A进行SVD分解，获得U、V和S矩阵；

（2）引入阻尼系数f，并根据Sk≦S1/SNR来确定矩阵S的大小；（SNR为信噪比）

（3）对S求逆时引入f，得到近似逆矩阵Si；

（4）设已知假定解初始值为X0=0；

（5）计算ΔMk=M-A*Xk，k是第k次迭代标号；

（6）根据M、U、V、Si直接计算||A*ΔXk-ΔMk||2的最小二乘解ΔXsk；

（7）计算Xk+1=Xk+ΔXsk；

（8）将Xk+1中小于零的项都改为零，并重复第4项步骤，直到假定解Xk所有元素都满足非负约束条件时，就将假定解Xk定为所要求取的反演解X；

（9）根据反演解X的每一个弛豫分量对应的T2值，得到T2时间谱。

5 结束语

目前国内地面核磁共振方法中自旋回波信号的提取过程尚存在问题，而国际上也仅有美国和法国少数学者对此方向进行了理论探索。正因为自旋回波信号在地面核磁共振方法中的应用尚属空白，所以一旦该技术成功投入到实际应用中，其意义必将重大。而本文初步提出的反演算法依然存在不足之处，就是对去噪后的原始数据信噪比的要求非常高，在低信噪比的条件下可用奇异值数量十分有限，因此对于算法的改进还有待深入研究。

[1]王鹏，李振宇.地面核磁共振-垂向电测深组合找水模式[J].地质科技情报，2006（03）：105-108.

[2]熊国欣，李立本.核磁共振成像原理[M].科学出版社，2007.

[3]E.L.Hahn.Spin Echoes.PHYSICAL REVIEW，1950（04）：580-594.

[4]吴云.磁性不均匀对SNMR信号的影响及克服方法研究[D].中国地质大学（武汉），2013.

[5]于蕾.利用SE求取SNMR T2时间方法研究[D].中国地质大学（武汉），2012.

[6]王才志，尚卫忠.应用奇异值分解算法的核磁共振测井解谱方法[J].石油地球物理勘探，2003，38（01）：91-94.

[7]姜瑞忠，姚彦平，苗盛等.核磁共振T2谱奇异值分解反演改进算法[J].石油学报，2005，26（06）：57～59.

[8]李庆谋，成秋明.分形奇异（特征）值分解方法与地球物理和地球化学异常重建[J].地球科学-中国地质大学学报，2004，29（1）：109～118.

[9]邓克俊，谢然红等.核磁共振测井理论及应用[M].山东省东营：中国石油大学出版社，2010.