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

陳詩揚

（中國地質大學（武漢）地球物理與空間信息學院，武漢 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.