基于Backus-Gilbert理論的孔隙介質核磁共振弛豫反演

肖立志，張恒榮，廖廣志，傅少慶，李 奎

油氣資源與探測國家重點實驗室中國石油大學，北京 102249

1 引 言

巖石的NMR弛豫測量能夠提供巖石物理特性和流體特性的許多信息，如巖石孔徑分布、孔隙度、束縛水、可動流體以及滲透率等.NMR測井從井下采集到的是一系列反映地層弛豫特性的回波串，需要經(jīng)過回波串多指數(shù)擬合才能得到T2譜.回波串擬合是一個數(shù)學反演過程[1-3］.國內(nèi)外學者提出了多種反演方法，如Butler等[4］基于非負約束的模平滑方法（BRD）；Prammer[5］基于奇異值分解的反演方法（SVD）；Dunn等[6］對比了模平滑方法和曲率平滑方法的優(yōu)缺點，發(fā)展了奇異值分解反演方法；Liaw等[7］和 Miller等[8］分別提出基于B樣條函數(shù)和Gamma函數(shù)的基函數(shù)方法；Borgia和Brown[9］設計了一個“曲率平滑”的復雜方法；Rafael和Boqin[10］等提出了蒙特卡羅最優(yōu)反演方法.國內(nèi)王才志[11］，姜瑞忠等[12］嘗試了奇異值分解方法；王忠東等[13］嘗試了基于整體迭代修正的反演方法；翁愛華[14］嘗試了高分辨率反演方法，譚茂金等[15］嘗試了基于遺傳算法的反演方法，馬建海等[16］嘗試了迭代Tikhonov正則化反演方法.這些算法嘗試從不同角度解決核磁共振弛豫測量的多指數(shù)反演問題，但通常會受低信噪比、正則化參數(shù)、迭代收斂等的影響.

本文基于Backus-Gilbert理論，嘗試一種新的孔隙介質NMR弛豫反演方法，簡稱BGW，其核心是從解的非唯一性出發(fā)，構造出一系列解估計，然后通過引入解估計分辨率和解估計方差對解進行評價，最終在其中找出一個最佳折衷的解估計.實踐表明，該方法不僅穩(wěn)定，并且適應于低信噪比數(shù)據(jù).

2 孔隙介質NMR弛豫反演問題描述

NMR測井采集包含數(shù)百個乃至數(shù)千個回波的CPMG自旋回波串信號.經(jīng)過一些預處理，如相位旋轉和累加等，然后進行多指數(shù)反演.回波測量值M（ti）可以寫成：

式中M（ti）表示第i個回波幅度（i=1，…，n），ti是第n個相等回波間隔的時間，εi代表第n個回波的噪聲，并有m個已知的預先選擇的弛豫時間T2在對數(shù)刻度上等間隔分布，a（T2j）是由方程（1）解出的弛豫時間T2j對應的幅度，可以刻度成區(qū)間孔隙度.

從數(shù)學角度，求解a（T2j）分布是一個反拉普拉斯變換[17-18］.噪聲的存在會導致反演結果的非唯一性.一個可行的方法是，利用最小二乘擬合使下面的求和最小化：

3 基于BG線性評價理論的反演方法

BG反演理論是地球物理反演的重要成果，其核心思想是從反問題解的非唯一性出發(fā)，對有限而精確數(shù)據(jù)的反演，反問題的解雖不唯一，但解的平均可以是唯一的；對有限帶誤差數(shù)據(jù)的反演，反問題的解可能不存在，即使存在也不唯一，這時要用分辨率和方差來評價解估計.由于解估計的方差和分辨率不可能同時達到極小，只能通過修改折衷參數(shù)的辦法求在分辨率和方差之間取最佳折衷的解估計[19-21］.

BG反演理論討論的是連續(xù)模型，與BG理論對應的離散模型反演由Wiggins和Jackson先后提出[22］.

NMR數(shù)據(jù)反演問題是離散非線性反問題，采集到的回波串是有限帶誤差的數(shù)據(jù).方程（1）的數(shù)值解將該方程變換為線性方程組的形式如下：

其中Aij，fj和gi是它們各自相應的連續(xù)函數(shù)的離散值，gi為帶誤差的回波幅度值，Aij為預先設置基函數(shù)組成的系數(shù)矩陣，fj為待定值，即待求的T2分布.

根據(jù)奇異值分解定理，n×m實矩陣A可分解為n×n的正交陣U和m×m的正交陣V和一個對角陣S：

由于系數(shù)矩陣A的病態(tài)性，其中sq+1，…，sm很小，根據(jù)離散BG理論，模型的協(xié)方差矩陣為

此時數(shù)據(jù)分辨率矩陣為

對應分辨率矩陣的解估計展布為

從理論上講，希望方差最小，分辨率最高.根據(jù)BG理論，由于解估計的方差和分辨率不可能同時達到極小和極大，只能通過修改折衷參數(shù)的辦法求在分辨率和方差之間取最佳折衷的解估計，而保留奇異值個數(shù)q則為折衷參數(shù).隨著折衷參數(shù)q的變化，解估計以及解估計的方差和展布也隨之變化，q從小到大逐漸變化，解估計的方差越來越大，展布越來越小，然后用（6）和（8）式計算方差和展布，以展布為橫坐標，方差為縱坐標繪出最佳折衷曲線圖，最佳折衷的q值應位于最接近原點的左下方，因為原點處展布和方差都為零.然后根據(jù)求出解估計中第k個分量，解估計中每一個分量都對應一個最佳折衷q值，而q值可能相同，也可能不相同.

傳統(tǒng)正則化方法是對整個fj向量直接求取，在選擇正則化參數(shù)時無論是人工還是自動都是對整個向量進行處理.例如，SVD方法是取一個截斷奇異值應用到每個分量；BRD方法是取一個平滑值加到每個分量；SIRT是基于整個向量的迭代修改.而BGW方法是分步求取每個分量，對每個分量都構造出一系列解估計，然后引入解估計分辨率和解估計方差來評價所有解估計，找出最佳折衷解估計，因此本方法在精度上更勝一籌，只是由于每個分量都要依次解出，反演時間上會有所增加.

4 基于BG理論的NMR弛豫反演方法

首先給定一個具有雙峰特征的T2分布模型，如圖1中的實線，根據(jù)這個模型生成回波串，設回波間隔TE=1.2ms，回波個數(shù)NE=1024，對數(shù)均勻布點32個，T2=0.3～3000ms，加入信噪比為100的高斯分布隨機噪聲（信噪比定義為第一個回波幅度與噪聲標準差的比值），用此模型考察BGW反演方法.

根據(jù)BG理論的反演方法，求出每一個弛豫分量的方差和展布，本文以第5個弛豫分量為例，將其方差和展布繪成折衷曲線，如圖2所示.

圖2中橫坐標（線性）為解估計的展布，縱坐標（對數(shù)）為解估計的方差；最佳折衷參數(shù)為6，計算此點與原點的距離最小.32個弛豫分量可以依次計算出與最佳折衷參數(shù)對應的最佳折衷解，即可得到T2分布.T2反演結果如圖1的圈實線.可見反演結果與模型吻合很好，從而證實了該方法的有效性.

下面通過構造不同模型，如大孔單峰、小孔單峰、大孔占優(yōu)雙峰和小孔占優(yōu)雙峰，并加不同噪聲，來考察BGW方法的適應性.

從圖3中看出，從單一孔隙結構到復雜孔隙結構，從信噪比高達100到信噪比低達10，BGW均能給出較好的結果.

5 方法對比

將BGW方法分別與SVD、BRD以及SIRT方法反演結果進行對比，仍然采用大孔占優(yōu)T2分布模型，信噪比SNR分別為100，30，20，10（見圖4）.

各個反演方法計算結果的相對誤差如表1表示，BGW方法在各信噪比下的相對誤差都是最小的.

圖4 大孔占優(yōu)模型不同信噪比不同方法反演結果（a）BGW；（b）SVD；（c）BRD；（d）SIRT.Fig.4 Inversion results of dominant large pore model with different SNR（a）BGW；（b）SVD；（c）BRD；（d）SIRT.

圖5是模擬某地層NMR T2分布.第一道是模型，在此模型基礎上產(chǎn)生回波串，并加入噪聲，使其信噪比為10，然后分別用四種方法進行反演.第二道是BGW反演結果；第三道是SIRT反演結果；第四道是SVD反演結果；第五道是BRD反演結果.可見，在信噪比為10的情況下，BGW反演結果與模型最為接近，說明在低信噪比情況下，BGW方法能得到更為有效的反演結果.

表1 4種方法反演結果的相對誤差比較（%）Table 1 Comparison of the relative error of inversion results（%）

圖5 四種方法反演結果與地層模型的對比Fig.5 Comparison of formation model computed by four methods

6 BGW方法影響因素研究

由式（1）知，不同噪聲、回波間隔、回波個數(shù)及T2布點都可能影響反演結果.

6.1 噪聲

信噪比低是NMR測井的重要特征，所以，反演方法必須考慮它對噪聲的依賴性.

噪聲影響可以用 Hansen等[23］給出的離散Picard條件進行分析.離散Picard準則認為，若方程組Ax=y的傅里葉系數(shù)，〈ui，y〉趨于零的速度在平均意義下快于矩陣A的奇異值σi趨于零的速度，則稱該方程組滿足離散Picard準則（條件）.

Picard準則是檢驗模型的不適定性以及模型受噪聲污染情況的重要準則.對于一個給定的模型，若其滿足或部分滿足Picard準則，則可以通過適當?shù)奶幚砬蠼庵?；當模型的不適定性較強，或者數(shù)據(jù)受噪聲污染較嚴重，以至于完全違背這一準則，那么再好的方法對該模型也是徒勞的.

圖6顯示了信噪比分別為106、30、10、5、2.5時傅里葉系數(shù)的變化.可以看出，傅里葉系數(shù)隨信噪比的降低而逐步偏離奇異值曲線，即傅里葉系數(shù)趨于零的速度逐步慢于奇異值趨于零的速度.從本文大孔占優(yōu)雙峰模型中，信噪比低于5后，已經(jīng)不滿足離散Picard準則，也就是說在這個模型下，由于噪聲的影響，從數(shù)學角度已經(jīng)不能反演出真實可靠信息，必須要對數(shù)據(jù)濾波或者重新采集.

圖7是BGW方法反演結果，證明信噪比小于5以后反演結果不再具有實際意義.

6.2 回波間隔影響

NMR回波串數(shù)據(jù)在時間域采集，按回波間隔TE記錄N 個離散回波幅度，總時間為T=N×TE.分辨T2譜的能力依賴于信噪比（SNR）和采集模式.定義一個有效采集時間ΔT（T2）為T2組分從信號最強衰減到噪聲水平的時間間隔.則相應的有效回波個數(shù)NE（T2）=ΔT（T2）／TE，T2譜的分辨能力可以用下面公式表示[24］：

從式中可以看出：（1）信噪比越低，T2譜分辨能力越差；（2）相對于長弛豫組分，短弛豫組分更難分辨；（3）在信噪比和T2組分確定時，短的TE能提高T2譜分辨能力.

圖8是不同信噪比下BGW方法反演結果，顯示回波間隔對T2譜反演的影響，同時證明了短回波間隔能提高T2譜分辨能力，并且信噪比越低時效果越明顯.

6.3 回波個數(shù)與布點影響

回波個數(shù)會影響數(shù)據(jù)采集的質量，反演的T2譜可能出現(xiàn)不收斂.引起這種不收斂的主要原因是由于所設置的回波間隔較小而采集的回波個數(shù)又不夠，導致長弛豫組分的信號采集不充分.因此，信號采集時要設置合理的等待時間和足夠的回波個數(shù)，使長弛豫組分的極化信號充分采集.

T2譜的布點個數(shù)，即反演時預先設置的基函數(shù)個數(shù)，對多指數(shù)反演也會有影響.當布點個數(shù)較少時，弛豫時間譜的分辨率較低，不能精細的反映儲層孔隙結構；當布點個數(shù)較多時，可以提高弛豫時間譜對孔隙結構的分辨率，但是布點數(shù)越多計算速度就越慢.實驗表明，巖心分析時，布點數(shù)在32～256個之間較合適[25］.

7 結 論

（1）根據(jù)BG反演理論，引入解估計分辨率和解估計方差，每一個解分量會對應一系列折衷參數(shù)，其中必定有一個最佳折衷參數(shù)使得分辨率和方差達到最佳平衡，此時該點在折衷曲線圖中離原點最近.

（2）與以往正則化方法每個解分量都對應相同的正則化參數(shù)不同，BGW方法在求解每個解分量時，最佳折衷參數(shù)不一定都相同，使求解結果更精確.

（3）通過構造符合不同地層T2譜的模型，并與其它算法相比較，驗證了BGW方法在孔隙介質NMR弛豫反演中的有效性.

（4）BGW方法在求解T2譜時，雖然不存在迭代收斂問題，但每個弛豫分量都是依次計算，因此整體計算速度有待優(yōu)化.

（5）利用離散Picard準則分析了噪聲對反演的影響，發(fā)現(xiàn)BGW方法在本文大孔占優(yōu)雙峰模型中信噪比低于5后，難以從回波串中提取有效信息.短回波間隔有利于提高T2譜分辨能力，采集回波個數(shù)盡量使長弛豫組分采集完全，T2布點需要包含長、短弛豫組分，盡量橫跨三個數(shù)量級，使得反演結果更能反映粘土束縛水、毛管束縛水和可動流體信息.

（References）

[1]肖立志.核磁共振成像測井與巖石核磁共振及其應用.北京：科學出版社，1998.Xiao L Z.NMR Image Logging and NMR in Rock Experiments（in Chinese）.Beijing：Science Press，1998.

[2]Coates G R，Xiao L Z，Prammer M G.NMR Logging Principles and Applications.Texas：Gulf Publishing Company，1999.

[3]Dunn K J，Bergman D J，Latorraca G A.Nuclear Magnetic Resonance：Petrophysical and Logging Applications.Pergamon：Elsevier Science，2002.

[4]Butler J P，Reeds J A，Dawson S V.Estimating solutions of first kind integral equations with nonnegative constraints and optimal smoothing.SIAMJ Numer.Anal.，1981，18（3）：381-397.

[5]Prammer M G.NMR Pore Size Distributions and Permeability at the Well Site.SPE 28368SPE Annual Technical Conference and Exhibition，25-28September 1994，New Orleans，Louisiana.

[6]Dunn K J，Latorraca G A，Warner J L.On the calculation and interpretation of NMR relaxation time distribution.SPE28367，69th Annual SPE Technical Conference and Exhibition，New Orleans，1994：45-54.

[7]Liaw H K，Kulkarni R，Chen S.Characterization of fluid distributions in porous media by NMR techniques.AICHE J.，1994，42（2）：538-548.

[8]Miller A，Chen S H，Georgi D T，Vozoff K.A new method for estimating T2distributions from NMR measurements.Mag.Res.Imaging，1998，16（5-6）：617-619.

[9]Borgia G C，Brown R J，F(xiàn)antazzini P.Uniform-penalty inversion of multiexponential decay data.Journal of Magnetic Resonance，1998，132（1）：65-77.

[10]Salazar-Tio R，Sun B Q.Monte Carlo Optimization-inversion Methods for NMR，SPWLA 50th Annual Logging Symposium，June 21-24，2009.

[11]王才志，尚衛(wèi)忠.應用奇異值分解算法的核磁共振測井解譜方法.石油地球物理勘探，2003，38（1）：91-94.Wang C Z，Shang W Z.Application of singular value decomposition （SVD）in solution of T2relaxation spectra from nuclear magnetic resonance （NMR）log data.Oil Geophysical Prospecting （in Chinese），2003，38（1）：91-94.

[12]姜瑞忠，姚彥平，苗盛等.核磁共振T2譜奇異值分解反演改進算法.石油學報，2005，26（6）：57-59.Jiang R Z，Yao Y P，Miao S，et al.Improved algorithm for singular value decomposition inversion of T2spectrum in nuclear magnetic resonance.Acta Petrolei Sinica（in Chinese），2005，26（6）：57-59.

[13]王忠東，肖立志，劉堂宴.核磁共振弛豫信號多指數(shù)反演新方法及其應用.中國科學（G），2003，33（4）：323-332.Wang Z D，Xiao L Z，Liu T Y.New inversion method with multi-exponent and its application.Science in China （G）（in Chinese），2003，33（4）：323-332.

[14]翁愛華，李舟波，莫修文等.核磁共振測井數(shù)據(jù)高分辨反演方法研究.測井技術，2002，26（6）：455-459.Weng A H，Li Z B，Mo X W，et al.On high resolution inversion of NMR logging data.WLT （in Chinese），2002，26（6）：455-459.

[15]譚茂金，石耀霖，謝關寶.基于遺傳算法的核磁共振T2分布反演.測井技術，2007，31（5）：413-416.Tan M J，Shi Y L，Xie G B.NMR T2distribution inversion based on genetic algorithm.Well Logging Technology（in Chinese），2007，31（5）：413-416.

[16]馬建海，孫建孟，孫萌等.基于迭代Tikhonov正則化的核磁測井解譜方法研究.西南石油大學報（自然科學版），2009，31（1）：37-40.Ma J H，Sun J M，Sun M，et al.Tikhonov iteration method of echo-stripping T2relaxation spectra from nuclear magnetic resonance log data.Journal of Southwest Petroleum Institute （Science &Technology Edition） （in Chinese），2009，31（1）：37-40.

[17]Song Y Q，Venkataramanan L，Burcaw L.Determining the resolution of Laplace inversion spectrum.The Journal of Chemical Physics，2005，122：104104.

[18]Song Y Q.Resolution and uncertainty of Laplace inversion spectrum.Magnetic Resonance Imaging，2007，25（4）：445-448.

[19]Backus G E，Gilbert J F.The resolving power of gross earth data.Geophys.J.R.Astron.Soc.，1968，16（2）：169-205.

[20]Backus G E，Gilbert J F.Uniqueness in the inversion of inaccurate gross earth data.Philosophical Transactions，1970，266（1173）：123-192.

[21]楊文采.地球物理反演的理論和方法.北京：地質出版社，1996.Yang W C.Theory and Methods of Geophysical Inversion（in Chinese）.Beijing：Geological Publishing House，1996.

[22]Jackson D D.Interpretation of inaccurate，insufficient and inconsistent data.Geophys.，1972，28（2）：97-109.

[23]Hansen P C.Rank-Deficient and Discrete Ill-Posed Problems：Numerical Aspects of Linear Inversion.Philadelphia：SIAM，1998.

[24]Chen S.Improving Estimation of NMR Log T2cutoff Value.SCA，2000.

[25]廖廣志，肖立志，謝然紅等.孔隙介質核磁共振弛豫測量多指數(shù)反演影響因素研究.地球物理學報，2007，50（3）：932-938.Liao G Z，Xiao L Z，Xie R H，et al.Influence factors of multi-exponential inversion of NMR relaxation measurement in porous media.Chinese J.Geophys.（in Chinese），2007，50（3）：932-938.