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

陳詩揚

（中國地質(zhì)大學（武漢）地球物理與空間信息學院，武漢 430074）

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

陳詩揚

（中國地質(zhì)大學（武漢）地球物理與空間信息學院，武漢 430074）

當測區(qū)磁場分布不均勻時，地面核磁共振技術(shù)所反映出的地下孔隙度會相對真實情況偏小。而將自旋回波信號應(yīng)用于地面核磁共振方法中，則會有效克服磁場的不均勻分布對該方法的影響。因此，本文嘗試提出一種利用自旋回波信號求取地面核磁共振方法中橫向弛豫時間譜的反演算法，為今后將自旋回波技術(shù)應(yīng)用于地面核磁共振方法的研究提供參考資料。本算法主要數(shù)學模型為多指數(shù)反演模型，反演算法為奇異值分解法。

地面核磁共振；自旋回波；多指數(shù)反演；奇異值分解

1 地面核磁共振基本原理

在靜磁場B0的作用下，具有磁矩μ的原子核將受到力矩μ×B0的作用，沿磁場B0方向為軸（設(shè)為Z軸方向），以ω0=γ|B0|做拉莫爾（Larmer）進動，γ為磁旋比，其中μ與B0方向之間的夾角保持固定不變[1]。

現(xiàn)在垂直于B0方向上施加一個射頻場B1，其振動角頻率為ω0。在B1作用下，自旋磁矩μ繞射頻磁場B1方向的進動，且與Z軸方向之間的夾角逐漸增大，這個夾角又叫扳轉(zhuǎn)角，其大小與γ、B1和它持續(xù)作用的時間τ（也叫脈沖寬度）的乘積呈正相關(guān)。射頻場B1的施加，使得原自旋系統(tǒng)磁化強度的縱向分量減小，同時磁化強度的水平分量增加。由于吸收了射頻場B1的能量，原子由穩(wěn)定的低態(tài)變成了不穩(wěn)定的激發(fā)態(tài)。因此當射頻場B1撤去后，原子會發(fā)生躍遷，向外輻射電磁波，而這個輻射出來的電磁波信號，即是我們通常所要測量的自由感應(yīng)衰減信號E（t，q）[2]。

當交變電流脈沖結(jié)束作用后，由同一線圈接收天線接收到自由感應(yīng)衰減信號按指數(shù)規(guī)律衰減。以下是地面核磁共振方法的基本方程：

其中，M0是氫質(zhì)子的磁化強度；是激發(fā)場垂直于地磁場的分量；r為激發(fā)半徑；I0為電流脈沖的振幅；α為地磁場傾角。根據(jù)式（1），可由以上參數(shù)與觀測值E（t，q），求出表示地下含水量的與表示地下孔隙度的

2 自旋回波[3]信號介紹

當穩(wěn)定磁場B0為不均勻磁場時，不同位置上自旋的同種原子核受到的磁場強度不同，因此它們以不同的進動頻率繞著B0作進動，具有一個總磁矩M。通過用90°—τ—180°—2τ—180°—2τ—……的脈沖序列，使總磁矩M沿XZ平面多次翻轉(zhuǎn)，其中τ為相鄰脈沖之間的間隔。設(shè)進動頻率較快的原子核的相位矢量為A，較慢的相位矢量為B，當向自選系統(tǒng)施以90°脈沖之后，A和B因為進動頻率不同而散開。τ時間之后施以180°脈沖，自旋系統(tǒng)繞y軸翻轉(zhuǎn)了180°，因此原來進動的快的A現(xiàn)在落到了B后面追趕B，又經(jīng)過τ的間隔A與B重聚，重聚時便產(chǎn)生自旋回波信號。隨著脈沖序列的不斷進行，該自選系統(tǒng)還會以同樣的形式產(chǎn)生一系列自旋回波信號，這一系列自旋回波信號稱之為CPMG回波串。

在不均勻磁場中，由于不同原子核進動速率不同，使得自旋系統(tǒng)發(fā)生散相，在通過180°翻轉(zhuǎn)脈沖后，分離的相位仍可在重聚后得到恢復。因此，我們可以通過一連串脈沖序列得到一組CPMG回波串，其T2*衰減曲線上每一個相位重聚時刻所對應(yīng)的點就是沒有因為磁場不均勻而產(chǎn)生散相的點，即真實T2衰減曲線上的點。因此，我們可以通過將這些零散的有用點擬合成曲線，得到真實的T2衰減曲線。

3 磁場梯度對地面核磁共振信號的影響

當存在磁場梯度時，地面核磁共振方法中的 并不等于真實值，因此無法反映出地下孔隙度的實際情況。與真實橫向弛豫時間 的關(guān)系如式（2）所示[4]：

由于實際測區(qū)磁場分布通常具有不均勻性，一般很難通過（2）式對所測T2*值進行高精度校正。而自旋回波信號中的橫向弛豫時間不會受到磁場不均勻性的影響，因此將自旋回波技術(shù)應(yīng)用到地面核磁共振方法中，能夠有效拓寬該方法的應(yīng)用范圍，例如可在火成巖發(fā)育地區(qū)進行找水工作。

4 自旋回波信號的T2解譜算法

對于復雜地質(zhì)條件中的自旋回波串，其曲線不再是單一的一條簡單的衰減曲線，而是多個指數(shù)衰減函數(shù)疊加而成的復雜曲線。因此，我們要從這個復雜的衰減曲線里提取出來的不是單個T2值，而是組成該曲線的每一個指數(shù)衰減函數(shù)所對應(yīng)的T2值的集合：T2譜。

（1）反演模型

由于原始回波串是由多個指數(shù)衰減函數(shù)所疊加而成的，因此本文中選取的數(shù)學模型為多指數(shù)反演模型。又因為往往原始回波串中是存在噪聲的，所以擬合的式子可以如下表達：

其中Mj（t）是tj時刻回波信號觀測幅度，m為每個回波弛豫分量數(shù)，Xi為第i種弛豫分量零時刻的信號大小，noise為噪聲[5]。

（2）T2時間解譜算法

假設(shè)觀測到的回波數(shù)為n個，弛豫分量數(shù)為m，那么原擬合式子就可以寫成M=AX的矩陣乘積形式[6-8]，如上所示。其中，M為所有回波觀測數(shù)據(jù)的矩陣，X為所有弛豫分量的矩陣，A為系數(shù)矩陣。首先考慮AX=M的攝動問題。當M因為存在噪聲而攝動時，必然要引起X的攝動，設(shè)M經(jīng)過攝動后變?yōu)镸+ΔM，解X經(jīng)過攝動后變?yōu)閄+ΔX，可以推知：

取T2i=2i+1，考慮到自擴散的影響，回波間隔不能太大，因此取Δt=0.5ms，最后求得cond（A）≈105，由（4）式可知當M因為存在噪聲而存在一絲毫攝動時，都會引起X產(chǎn)生巨大的變化。因此，這是一個病態(tài)問題。為了解決這個病態(tài)問題[9]，現(xiàn)利用奇異值分解法求解ΔX的最小二乘解[7]，這里的ΔX指的是真實解X與假定的已知解之差，并在求解時對S矩陣求逆的過程中引入阻尼系數(shù)f來使解穩(wěn)定。引入阻尼系數(shù)f將原本對S矩陣求逆的式子Si（k，k）=1/S（k，k）變?yōu)镾i（k，k）=S（k，k）/（S（k，k）^2+f^2），Si是所求的近似逆矩陣。

以下是解譜算法的主要步驟：

（1）對矩陣A進行SVD分解，獲得U、V和S矩陣；

（2）引入阻尼系數(shù)f，并根據(jù)Sk≦S1/SNR來確定矩陣S的大小；（SNR為信噪比）

（3）對S求逆時引入f，得到近似逆矩陣Si；

（4）設(shè)已知假定解初始值為X0=0；

（5）計算ΔMk=M-A*Xk，k是第k次迭代標號；

（6）根據(jù)M、U、V、Si直接計算||A*ΔXk-ΔMk||2的最小二乘解ΔXsk；

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

（8）將Xk+1中小于零的項都改為零，并重復第4項步驟，直到假定解Xk所有元素都滿足非負約束條件時，就將假定解Xk定為所要求取的反演解X；

（9）根據(jù)反演解X的每一個弛豫分量對應(yīng)的T2值，得到T2時間譜。

5 結(jié)束語

目前國內(nèi)地面核磁共振方法中自旋回波信號的提取過程尚存在問題，而國際上也僅有美國和法國少數(shù)學者對此方向進行了理論探索。正因為自旋回波信號在地面核磁共振方法中的應(yīng)用尚屬空白，所以一旦該技術(shù)成功投入到實際應(yīng)用中，其意義必將重大。而本文初步提出的反演算法依然存在不足之處，就是對去噪后的原始數(shù)據(jù)信噪比的要求非常高，在低信噪比的條件下可用奇異值數(shù)量十分有限，因此對于算法的改進還有待深入研究。

[1]王鵬，李振宇.地面核磁共振-垂向電測深組合找水模式[J].地質(zhì)科技情報，2006（03）：105-108.

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

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

[4]吳云.磁性不均勻?qū)NMR信號的影響及克服方法研究[D].中國地質(zhì)大學（武漢），2013.

[5]于蕾.利用SE求取SNMR T2時間方法研究[D].中國地質(zhì)大學（武漢），2012.

[6]王才志，尚衛(wèi)忠.應(yīng)用奇異值分解算法的核磁共振測井解譜方法[J].石油地球物理勘探，2003，38（01）：91-94.

[7]姜瑞忠，姚彥平，苗盛等.核磁共振T2譜奇異值分解反演改進算法[J].石油學報，2005，26（06）：57～59.

[8]李慶謀，成秋明.分形奇異（特征）值分解方法與地球物理和地球化學異常重建[J].地球科學-中國地質(zhì)大學學報，2004，29（1）：109～118.

[9]鄧克俊，謝然紅等.核磁共振測井理論及應(yīng)用[M].山東省東營：中國石油大學出版社，2010.