基于奇異值分解的核磁共振T2譜反演研究

李 輝， 肖忠祥， 劉曉娟

(西安石油大學(xué)，陜西 西安 710065)

1946年哈佛大學(xué)的Purcell 和斯坦福大學(xué)的Bloch 兩位教授分別用吸收法和感應(yīng)法發(fā)現(xiàn)核磁共振(NMR)現(xiàn)象以來(lái)，經(jīng)過(guò)近60年的發(fā)展，核磁共振T2 譜已經(jīng)能夠提供許多的巖石儲(chǔ)集物性參數(shù)及流體物性參數(shù)(肖立志，1998)，如巖石孔隙度、有效孔隙度、孔隙尺寸分布、滲透率、可動(dòng)流體與束縛流體飽和度、粘度、流體擴(kuò)散系數(shù)等。而這些參數(shù)都是由測(cè)量的巖石樣品及流體樣品的核磁共振弛豫信號(hào)經(jīng)過(guò)多指數(shù)反演后的T1、T2譜計(jì)算得到的。這些巖石物理信息為儲(chǔ)層評(píng)價(jià)、產(chǎn)量預(yù)測(cè)等提供了重要的理論依據(jù)。

目前，國(guó)內(nèi)外先后提出的核磁共振T2譜反演方法主要有三種:罰函數(shù)法(BRD)(Butler et al.，1981)、奇異值分解法(SVD)(王忠東等，2003)以及聯(lián)合迭代重建算法(SIRT)。

早期解譜方法多為單指數(shù)或多指數(shù)擬合，Whitall 等(1989)提出NNLS 算法，美國(guó)Numar 公司提出的商業(yè)化的基于正則化方法的UPEN 方案，Gy.Steinbrecher 提出的基于Laplas 變化的反卷積法(姚緒剛等，2003)，這些方法都存在著各種缺點(diǎn)，反演結(jié)果不理想，需要做進(jìn)一步改進(jìn)。

本文通過(guò)對(duì)傳統(tǒng)奇異值算法進(jìn)行改良，實(shí)現(xiàn)低信噪比的可操作性。

1 核磁共振T2 譜反演原理

1.1 基本原理

核磁共振測(cè)井的原始數(shù)據(jù)是幅度隨時(shí)間衰減的回波信號(hào)，這些回波信號(hào)隱含重要的巖石物理和油氣信息，并且能夠反映橫向弛豫過(guò)程，必須經(jīng)過(guò)一個(gè)基本的多指數(shù)變化，得到反映孔徑特征的T2譜分布，這一過(guò)程稱為核磁共振測(cè)井的解譜。如果核磁共振檢測(cè)到的信號(hào)很微弱，便會(huì)出現(xiàn)信噪比低的情況，在解譜過(guò)程中必須考慮這一問(wèn)題，提高信噪比，突出有用信號(hào)。

理論和實(shí)驗(yàn)研究都表明，單個(gè)孔隙內(nèi)磁化強(qiáng)度信號(hào)的橫向弛豫過(guò)程是按照單指數(shù)規(guī)律遞減的(Coates et al.，2007)。測(cè)井過(guò)程中，由于探測(cè)區(qū)域存在一系列大小不同的孔隙群體，孔隙中存在多種弛豫組分，測(cè)量得到的信號(hào)就是這些不同弛豫組分信號(hào)的迭加，即總磁化強(qiáng)度(Total Magnetization)是一系列指數(shù)衰減信號(hào)的迭加(肖立志等，2010)(圖1，2)。

1.2 計(jì)算過(guò)程

圖1 指數(shù)衰減信號(hào)曲線Fig.1 Exponential decay signal curve

圖2 NMR 原始信號(hào)—回波串Fig.2 The original NMR signals-Echo string

T2譜橫向弛豫信號(hào)分布滿足Fredholm 第一類積分方程，但實(shí)際過(guò)程中，由于噪聲的存在，弛豫過(guò)程滿足下面的方程:

其離散化模型為:

其中，Z(tj)為tj觀測(cè)到的回波幅度;T2i為第i個(gè)弛豫分量的橫向弛豫時(shí)間;T2min和T2max分別為測(cè)量的自由感應(yīng)衰減信號(hào)所能辨別的最短和最長(zhǎng)弛豫時(shí)間;xi為第i 種弛豫分量零時(shí)刻的信號(hào)大小;n 為弛豫分量個(gè)數(shù);tj為時(shí)間(通常為回波間隔的整數(shù)倍);εi為噪聲信號(hào)。

解譜過(guò)程就是如何從(2)中分解出T2i和xi。

弛豫分量的橫向弛豫時(shí)間T2i一般取2i+1(i =1，2，…)ms。假設(shè)弛豫分量的個(gè)數(shù)為n，觀測(cè)的回波個(gè)數(shù)為m，則由式(2)寫(xiě)成矩陣形式如下:

其中

在式(3)左右兩邊同時(shí)左乘矩陣的廣義逆矩陣A+，可得

其中，A+為n × m 矩陣，且由最小二乘法知A+=(ATA)-1AT。

根據(jù)矩陣?yán)碚摰钠娈愔刀ɡ?肖立志，1998)可知，對(duì)于A ∈，r(A)= r，存在m 階酉矩陣U和n 階酉矩陣V，使

其中，∑r= diag(δ1，δ2，…，δr)，δ1，δ2，…，δr是A 的全部非零奇異值，且呈遞減排列。為m × n矩陣，于是有，

不論A 是否奇異，總可以分解成(5)形式，并且∑r是唯一的。但通常情況下，由式(7)解出的X 值不穩(wěn)定，且常有負(fù)數(shù)出現(xiàn)。根據(jù)核磁共振的解譜結(jié)果的物理意義可知，xi為第i 種弛豫分量在零時(shí)刻的初始信號(hào)，不能為負(fù)值，而一般純碎的數(shù)學(xué)過(guò)程不能保證這一點(diǎn)。

1.3 數(shù)值解穩(wěn)定性控制

本文采用對(duì)系數(shù)矩陣A 加阻尼項(xiàng)的方法使數(shù)值解穩(wěn)定，并結(jié)合迭代算法進(jìn)行非負(fù)性約束。

1.3.1 阻尼方法

根據(jù)式(3)構(gòu)造函數(shù)

其中，λ2稱為阻尼因子。為求解F 最小值，令F 梯度為零，有

所以

又因

所以

對(duì)式(12)求逆矩陣得

綜上可得

另外，給定x 一個(gè)初始解X = X0，得到Z0=AX0，與Z = AX 相減，可得Z - Z0= A(X - X0)，記為

與上述推理類似，能夠得出

同時(shí)，需對(duì)ΔX 進(jìn)行截?cái)嗵幚?，?/p>

信噪比SNR 是由原始資料通過(guò)下述公式確定的

而由X = X0+ ΔX 即可得到X。

最終，當(dāng)SNR 降低時(shí)，對(duì)角線上保留較少項(xiàng)數(shù)，反之，SNR 提高時(shí)，對(duì)角線上保留較多項(xiàng)數(shù)。從而減少了噪聲對(duì)結(jié)果的影響，增加了解的穩(wěn)定性，在不同信噪比情況下都可以得到相對(duì)穩(wěn)定的T2譜。

1.3.2 非負(fù)性約束實(shí)現(xiàn)

(1)給定初始值X(k)，令k = 0;

(4)計(jì)算X(k+1)= X(k)+ ΔX(k);

(5)若X(k+1)中分量全不為負(fù)，則終止，輸出X(k+1)，否則轉(zhuǎn)(6);

(6)令X(k+1)中負(fù)值為0，k = k +1，轉(zhuǎn)(2)再次迭代，直到滿足條件(5)。

2 反演效果檢驗(yàn)

為驗(yàn)證該方法效果，用微機(jī)編制相關(guān)程序，并對(duì)其進(jìn)行數(shù)值模擬。圖2 中實(shí)驗(yàn)T2譜的布點(diǎn)數(shù)為32，其余實(shí)驗(yàn)中T2譜的布點(diǎn)數(shù)為64，布點(diǎn)區(qū)間為0.1 ～10 000 ms。

反演結(jié)果和構(gòu)造T2譜對(duì)比結(jié)果(圖3，4)。其中圖3 為32個(gè)布點(diǎn)的反演無(wú)噪回波數(shù)據(jù)T2譜與構(gòu)造譜對(duì)比圖，圖3 為64個(gè)布點(diǎn)反演無(wú)噪回波數(shù)據(jù)T2譜與構(gòu)造譜對(duì)比圖?？梢钥闯觯囱莸那€和數(shù)據(jù)點(diǎn)與構(gòu)造譜的曲線和數(shù)據(jù)點(diǎn)基本重合。這表明，本文的奇異值分解法可以較準(zhǔn)確反映真實(shí)的T2譜。

分別加入信噪比SNR =30，SNR =60，SNR =100 的高斯白噪聲后的數(shù)據(jù)反演結(jié)果對(duì)比(圖5)。從圖中可看出，三種情況下反演結(jié)果與構(gòu)造譜分布趨勢(shì)一致，譜線光滑并連續(xù)，具有較高吻合度。同時(shí)還可以看出，信噪比越高，反演T2譜與構(gòu)造T2譜的符合率越高。該算法具有較高的穩(wěn)定性和收斂性。

用該方法處理某巖心實(shí)驗(yàn)室NMR 數(shù)據(jù)，解譜結(jié)果與實(shí)驗(yàn)室解譜結(jié)果對(duì)比(圖6)。可以看出，巖心奇異值分解法反演結(jié)果與實(shí)驗(yàn)室反演結(jié)果吻合性較好，能夠得到連續(xù)的T2譜。

3 結(jié)論

此種SVD 法能夠獲得連續(xù)性較好的反演T2譜，算法較穩(wěn)定，收斂性好，非負(fù)約束條件容易實(shí)現(xiàn)。反演T2譜與構(gòu)造T2譜吻合性較高。由NMR實(shí)驗(yàn)室數(shù)據(jù)反演結(jié)果可知，該方法可以應(yīng)用于實(shí)驗(yàn)室?guī)r心的T2譜反演。

圖3 無(wú)噪聲反演結(jié)果對(duì)比(32個(gè)布點(diǎn))Fig.3 Inversion results of NMR data without noise and origin T2 spectrum(32 stationing)

圖4 無(wú)噪聲反演結(jié)果對(duì)比(64個(gè)布點(diǎn))Fig.4 Inversion results of NMR data without noise and origin T2 spectrum(64 stationing)

圖5 有噪聲反演結(jié)果對(duì)比Fig.5 Inversion results of three groups of NMR data with noise

圖6 巖心反演結(jié)果對(duì)比Fig.6 Inversion results of NMR data from core using Improved SVD and Laboratory solution spactrum

王忠東，肖立志，劉堂宴.2003.核磁共振馳豫信號(hào)多指數(shù)反演新方法及其應(yīng)用［J］.中國(guó)科學(xué)(G 輯)，33(4):323-324.

肖立志，柴細(xì)元，孫寶喜，等.2010.核磁共振測(cè)井資料解釋與應(yīng)用導(dǎo)論［M］.北京:石油工業(yè)出版社:15-38.

肖立志.1998.核磁共振成像測(cè)井與巖石核磁共振及其應(yīng)用［M］.北京:科學(xué)出版社:11-18.

姚緒剛，王忠東.2003.一種新的核磁共振弛豫譜反演算法［J］.測(cè)井技術(shù)，27(5):373-376.

Butler J P，Reeds J A，Dawson S V.1981.Estimating Solutions of First Kind Integral Equations with Nonnegative Constraints and Optimal Smoothing［J］.SIAMJ Numer Anal.，18(3):381-397.

Coates G，肖立志，Prammer M. 2007. 核磁共振測(cè)井原理與應(yīng)用［M］.北京:石油工業(yè)出版社:39-41.

Whitall K P，Mackay A I.1989.Quantitative Interpretation of NMR Relaxation Data［J］. Journal of Megnetic Resonance，4(5):134-152.