基于奇異值分解的核磁共振測(cè)井T2 譜反演方法的改進(jìn)

李鵬舉，葛 成，孫國(guó)平，陳 新，王彥凱

（1.大慶石油學(xué)院地球科學(xué)學(xué)院，黑龍江大慶163318；2.大慶鉆探工程公司測(cè)井二分公司，吉林松原138000；3.大慶油田勘探開發(fā)研究院，黑龍江大慶163712）

0 引 言

核磁共振（NMR）測(cè)井弛豫信號(hào)的多指數(shù)在NMR巖心分析與測(cè)井解釋反演結(jié)果的好壞直接影響到后續(xù)計(jì)算儲(chǔ)層物性參數(shù)（如孔隙度、滲透率、孔徑分布等）和T2截止值的準(zhǔn)確性［1］。Butler等［2］、Dunn等［3］、王忠東和肖立志［4］先后討論了罰函數(shù)法（BRD）、奇異值分解算法（SVD）以及聯(lián)合迭代重建算法（SIRT）。Borgia等在罰函數(shù)法的基礎(chǔ)上討論了均勻懲罰的反演方法。Jonathan等用線性正則化的方法分析了由短弛豫組分構(gòu)成的多指數(shù)弛豫數(shù)據(jù)的擬合問題，這些算法從不同的正則化角度給出了多指數(shù)反演的設(shè)計(jì)思路。此外，Mohnke等［5］利用模擬退火法優(yōu)化反演中的隨機(jī)搜索問題；Dunn等［6］討論了對(duì)不同測(cè)量誤差的原始回波數(shù)據(jù)進(jìn)行組合反演的方法。本文基于矩陣分析理論，進(jìn)一步改進(jìn)奇異值分解反演方法。

1 核磁共振測(cè)井T2譜解譜方法

1.1 建立儲(chǔ)層流體核磁共振模型

求解弛豫過程的Bloch方程可知，氫原子核系統(tǒng)磁化強(qiáng)度矢量的橫向分量是按指數(shù)規(guī)律衰減的，由CPMG脈沖序列測(cè)得的回波串也按指數(shù)規(guī)律衰減。儲(chǔ)層巖石通常存在1個(gè)孔隙尺寸分布，并且常常含有多種流體成分，此時(shí)孔隙中存在多種弛豫組分，即橫向弛豫時(shí)間常數(shù)（T2）不是單值，而是1個(gè)T2譜分布。因此，由CPMG脈沖序列測(cè)量記錄的自旋回波串按多指數(shù)規(guī)律衰減，即各單指數(shù)衰減的疊加

式中，m、n分別為測(cè)量到的回波個(gè)數(shù)和弛豫分量個(gè)數(shù)，i＝1，2，…，m；j＝1，2，…，n；ti為采集時(shí)間（通常是回波間隔（Te）的整數(shù)倍］；gi為測(cè)量得到的回波幅度；T2，j為預(yù)先選擇的弛豫時(shí)間分布；μi為測(cè)量的誤差；f（T2，i）為第i個(gè)弛豫分量對(duì)總的橫向弛豫信號(hào)的貢獻(xiàn)。展開式（1）可以得到方程組

把式（2）寫成向量的形式為

令

1.2 矩陣奇異值分解與方程組的解

在式（4）左右兩邊同乘A＋得

根據(jù)矩陣的奇異值分解定理［7］，對(duì)任意m×n階矩陣A，存在m階列正交矩陣U和n階正交矩陣V，滿足

把式（7）代入式（5）可得

1.3 奇異值分解反演算法的改進(jìn)

通常情況下，由式（8）解出的fi數(shù)值不穩(wěn)定且會(huì)有負(fù)值出現(xiàn)，而fi的物理意義為第i種弛豫分量在0時(shí)刻的初始幅度值，是非負(fù)值。本文采取對(duì)系數(shù)矩陣加阻尼項(xiàng)的處理方法使其數(shù)值解穩(wěn)定，應(yīng)用迭代法對(duì)解進(jìn)行非負(fù)性約束。

1.3.1 阻尼方法

根據(jù)式（4）建立式（10）目標(biāo)函數(shù)求其極小值

可得

所以

對(duì)式（15）求逆

把式（6）、式（16）代入式（13）可得

1.3.2 解的非負(fù)性約束

（1）置初始值fk，k＝0；（k代表進(jìn)行第k次迭代計(jì)算）

（2）dk?Afk；Δdk?d－dk；

（4）計(jì)算fk＋1?fk＋Δfk；

（5）如果（fk＋1中的分量全部≥0）停機(jī)，輸出fk＋1，否則進(jìn)行下一步；

（6）置fk＋1中的負(fù)值為0，k＝k＋1轉(zhuǎn)入第2步再次迭代，直到解滿足非負(fù)約束條件為止。

2 反演實(shí)例及效果分析

2.1 無(wú)噪聲理論數(shù)據(jù)反演及效果檢驗(yàn)

為了檢驗(yàn)該方法的解譜效果，構(gòu)造T2譜，按方程（1）正演模型生成了無(wú)噪聲理論回波串衰減曲線（見圖1）。反演結(jié)果和構(gòu)造的T2譜對(duì)比（見圖2、圖3）?？梢钥闯觯娈惤刂怪禐?0－6和奇異截止值為10－7反演的T2譜幾乎完全一致，與構(gòu)造的T2譜幾乎完全吻合。這表明，SVD法可以準(zhǔn)確地反演出真實(shí)的T2譜，奇異截止值對(duì)解譜效果影響小。

正演和反演的參數(shù)T2范圍為（1×10－3）～3s，按對(duì)數(shù)分別布11個(gè)和30個(gè)T2點(diǎn)；1 000個(gè)回波串；回波間隔為0.6×10－3s。

2.2 不同信噪比數(shù)據(jù)反演及效果檢驗(yàn)

給圖1回波串信號(hào)分別加入信噪比SNR＝5、25、50、100的高斯白噪聲，其中加入SNR＝5高斯白噪聲后的回波曲線見圖4。圖5給出不同SNR數(shù)據(jù)反演結(jié)果與構(gòu)造T2譜的對(duì)比。從圖5可以看出，信噪比SNR＝25、50、100的情況下，解譜結(jié)果與構(gòu)造譜吻合得很好；信噪比SNR＝5的情況下反演的T2譜和構(gòu)造譜分布趨勢(shì)一致，譜線光滑連續(xù)，大組分反演效果好，但小組分反演結(jié)果與構(gòu)造譜相差較大，分辨率低；信噪比越高，反演所得的T2譜與構(gòu)造譜的符合率越高。這表明，該方法反演高信噪比數(shù)據(jù)效果顯著，而對(duì)低信噪比數(shù)據(jù)的反演可以保持波形特點(diǎn)大體一致。

圖1 由構(gòu)造譜生成的無(wú)噪聲理論回波串衰減曲線

圖2 11個(gè)布點(diǎn)反演無(wú)噪回波數(shù)據(jù)T2譜與構(gòu)造譜對(duì)比圖

圖3 30個(gè)布點(diǎn)反演無(wú)噪回波數(shù)據(jù)T2譜與構(gòu)造譜對(duì)比圖

圖4 由構(gòu)造譜生成的信噪比SNR＝5的理論回波曲線

表1 古89-47井巖心NMR實(shí)驗(yàn)參數(shù)表

表2 古89-47井巖心SVD解譜反演參數(shù)表

2.3 實(shí)驗(yàn)室NMR數(shù)據(jù)反演及效果檢驗(yàn)

用該方法處理了大慶油田古89-47井1塊巖心實(shí)驗(yàn)室NMR數(shù)據(jù)，解譜結(jié)果與實(shí)驗(yàn)室解譜結(jié)果對(duì)比見圖6。巖心NMR實(shí)驗(yàn)參數(shù)和SVD解譜反演參數(shù)如表1和表2。

圖5 不同信噪比理論回波數(shù)據(jù)反演T2譜與構(gòu)造T2譜對(duì)比圖

圖6 古89-47井巖心SVD解譜與實(shí)驗(yàn)室解譜結(jié)果對(duì)比圖

圖7 古89-47井巖心測(cè)量的回波曲線及回波串的擬合曲線

從圖6可以看出，巖心SVD法反演結(jié)果與實(shí)驗(yàn)室反演結(jié)果吻合得很好，譜的形狀和變化趨勢(shì)一致。從圖7可以看出，SVD法多指數(shù)擬合曲線較好地?cái)M合了巖心實(shí)驗(yàn)室采集的回波串曲線和實(shí)驗(yàn)室方法擬合的回波串曲線。上述結(jié)果表明改進(jìn)的SVD法可以快速、有效地反演核磁共振T2譜，對(duì)提高核磁共振測(cè)井精細(xì)解釋符合率有重要意義。

3 結(jié) 論

經(jīng)過改進(jìn)，SVD法可以有效反演出真實(shí)的T2譜，加阻尼項(xiàng)可以使數(shù)值解穩(wěn)定。在布點(diǎn)數(shù)為11～128、信噪比大于25的情況下擬合精度較高，適用于高信噪比回波數(shù)據(jù)的反演，并且信噪比越高，反演的T2譜與構(gòu)造譜的符合率越高。大組分比小組分反演分辨率高。同時(shí)新方法采用迭代的思想進(jìn)行譜的非負(fù)性限制，而未刪減原始測(cè)量數(shù)據(jù)，保證了T2譜的真實(shí)性和連續(xù)性。核磁共振實(shí)驗(yàn)室數(shù)據(jù)反演結(jié)果表明，該方法可以應(yīng)用于實(shí)驗(yàn)室?guī)r心T2譜反演。

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

［2］ Butler J P，Reeds J A，Dawson S V.Estimating Solutions of First Kind Integral Equations with Nonnegative Constraints and Optimal Smoothing［J］.SIAMJ Numer Anal，1981，18（3）：381－397.

［3］ Dunn K J，Latorraca G A，Warner J L.On the Calculation and Interpretation of NMR Relaxation Time Distribution［C］∥SPE28367，69th Annual SPE Technical Conference and Exhibition，New Orleans，1994：45－54.

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

［5］ Mohnke O，Yaramanci U.Smooth and Block Inversion of Surface NMR Amplitudes and Decay Times Usingsimulated Annealing［J］.Journal of Applied Geophysics，2002，50：163－177.

［6］ Dunn K J，Latorraca G A.The Inversion of NMR Log Data Sets with Different Measurement Errors［J］.Journal of Magnetic Resonance，1999，140：153－161.

［7］ 張賢達(dá).矩陣分析與應(yīng)用［M］.北京：清華大學(xué)出版社，2004：341－398.