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

陳詩(shī)揚(yáng)

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

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

陳詩(shī)揚(yáng)

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

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

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

1 地面核磁共振基本原理

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

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

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

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

2 自旋回波[3]信號(hào)介紹

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

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

3 磁場(chǎng)梯度對(duì)地面核磁共振信號(hào)的影響

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

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

4 自旋回波信號(hào)的T2解譜算法

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

（1）反演模型

由于原始回波串是由多個(gè)指數(shù)衰減函數(shù)所疊加而成的，因此本文中選取的數(shù)學(xué)模型為多指數(shù)反演模型。又因?yàn)橥蓟夭ù惺谴嬖谠肼暤模詳M合的式子可以如下表達(dá)：

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

（2）T2時(shí)間解譜算法

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

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

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

（1）對(duì)矩陣A進(jìn)行SVD分解，獲得U、V和S矩陣；

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

（3）對(duì)S求逆時(shí)引入f，得到近似逆矩陣Si；

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

（5）計(jì)算ΔMk=M-A*Xk，k是第k次迭代標(biāo)號(hào)；

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

（7）計(jì)算Xk+1=Xk+ΔXsk；

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

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

5 結(jié)束語(yǔ)

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

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

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

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

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

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

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

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

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

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