二維核磁共振測井正反演模擬軟件設計與應用

張家成 張 宮* 黃若坤 覃瑩瑤

1(長江大學油氣資源與勘探技術教育部重點實驗室 湖北 武漢 430100)2(塔里木油田分公司勘探開發研究院 新疆 庫爾勒 841000)

0 引 言

核磁共振(NMR)測井通過觀測儲層中含氫核流體的NMR信號來獲得用于評價儲層參數的信息[1-2]。自20世紀50年代被引入到石油測井領域以來，核磁共振技術在飽和度、滲透率計算、流體識別、孔隙度、孔隙結構評價等方面均得到了成功應用[3-5]。2002年Sun等[6]和Hurlimann等[7]在石油測井領域引入了核磁共振波譜學中的二維核磁共振概念，進一步拓展了核磁共振測井的使用范圍，提高了儲層測井解釋評價的準確率[8]。核磁反演方法也從對單個回波串數據進行反演發展到同時對多個回波串數據進行反演[9]，但大量的數據加重了反演過程中的運算負擔。一直以來，國內外學者為滿足反演精度和反演速度上的要求，對二維核磁反演方法進行了不斷的探索。顧兆斌等[10]利用改進的奇異值分解法對二維核磁共振進行反演；謝然紅等[9]針對二維的多回波串數據處理問題提出了多回波串聯合反演方法；李鵬舉等[11]針對BRD和TSVD兩種反演方法的不足提出了變參量迭代的快速核磁反演方法。在二維核磁反演方法研究過程中，繁瑣的數據構造和結果評測過程增加了研究的時間成本。

本文從精度和效率出發，研發了一套二維核磁共振正反演模擬軟件，實現了核磁共振二維弛豫及反演的全過程，降低了二維核磁反演方法的研究難度。

1 原 理

1.1 核磁共振弛豫測量

核磁共振測井通過觀測流體中氫原子核的NMR信號，處理得到地層孔隙中流體的含量及其性質。巖石物理信息，如滲透率、孔徑分布、孔隙度和束縛水，也可以通過核磁共振弛豫測量得到。核磁共振弛豫過程有兩種，其中縱向弛豫計算如下：

(1)

式中：Pi是第i個弛豫信號占總的縱向弛豫信號的百分比；t表示時間；T1表示縱向弛豫時間，反映了自旋系統從外界吸收或者釋放能量的效率，主要通過反轉恢復脈沖序列或飽和恢復法脈沖序列進行測量[13]。

橫向弛豫計算如下：

(2)

式中：M(t)表示核磁共振測量得到的橫向弛豫信號；Pi是第i個弛豫信號占總的橫向弛豫信號的百分比；t表示時間；T2表示橫向弛豫時間，可以通過CPMG(Carr-Purcell-Meiboom-Gill)脈沖序列進行測量[13]。在一系列CPMG序列之間，等待一段時間確保所有的自旋完全恢復，再準備下一個激勵，這段時間稱為 “等待時間”。實際測量中，不同地層除了橫向弛豫時間T2不同外，縱向弛豫時間T1、擴散系數D、內部梯度G等均存在差異，會直接影響回波信號的衰減，用一組回波串進行反演往往不能充分地挖掘流體之間的差異。二維核磁共振采用多組回波串組合測量，通過聯合反演，得到二維分布譜，使不同流體在處理成果中有更好的區分度。

1.2 二維核磁弛豫構造

構造二維弛豫譜是采用對數坐標下的二維高斯分布公式進行計算[14],具體實現如下：

(3)

式中：j=1,2,…,m，m是構造譜的布點數目；Aj是第j組分的幅度；F是孔隙度刻度因子；T2j是第j組分的T2弛豫時間；T2mid是T2弛豫時間峰值的位置；Dj是第j組分的擴散系數；Dmid是擴散系數峰值的位置；σ1、σ2可以控制分布的展布寬度；ρ是T2j和Dj的相關系數。

為了更加真實地模擬出實際二維譜，二維構造譜必須能夠反映多種不同性質流體的疊加效果，所以式(3)需要進一步改進為：

(4)

式中：j=1,2,…,m，m是構造譜的布點數目；設具有4類不同性質的流體，Fg是第g種流體的孔隙度分量；HIg是第g種流體的含氫指數；T2g,mid是第g種流體的T2弛豫中心值；Dg,mid是第g種流體的擴散系數峰值；σ1g、σ2g分別是第g種流體T2峰值的中心展布寬度和擴散系數峰值的展布寬度；ρ是T2j和Dj的相關系數。

1.3 常見的二維核磁共振脈沖序列

針對二維核磁共振測井數據采集的特點，國內外學者均提出了不同的CMPG脈沖序列。Hurlimann等[7]提出了擴散編輯脈沖序列，該序列中時序采集分為兩個窗口，第一個窗口內僅采集兩個回波，通過調整這兩個回波的回波間隔來采集孔隙中流體的擴散弛豫信息，而第二個窗口則與改良式CPMG脈沖序列一樣，都是采用最小的回波間隔來采集孔隙流體的橫向弛豫信息，如圖1所示。

圖1 擴散編輯T2-D測量序列

在回波采集過程中，回波的幅度衰減計算如下：

(5)

式中：bik表示回波間隔為TElk時第一個窗口內的第i個回波的幅度；γ表示氫核的旋磁比；G表示磁場梯度；TE表示第二個窗口內的回波間隔；當擴散系數為Dp，橫向弛豫時間為T2j時的孔隙度分量用f(Dp,T2j)表示。該脈沖序列通過調整第一個窗口內的回波間隔，分別采集多道CPMG回波信號，從而實現T2-D模式下的二維數據采集，利用式(5)將采集獲得的回波數據進行聯合反演，便可得到地層孔隙中流體的T2-D分布。

謝然紅等[14]則利用回波間隔不同的多條CPMG回波串來完成T2-D二維數據的采集，相對于前一種方法更加簡便，利用現有的一維核磁共振測井儀器便可以完成該類數據的采集：

(6)

式中：bik是回波間隔為TEk時第i個回波的幅度[9]。

除此之外，還有通過調整第一個窗口內的延時來采集多條CPMG回波信號的T1編輯序列[4]以及測量多組不同等待時間的多TW序列等其他常見的脈沖序列[15]，這里不再贅述。

1.4 二維譜反演方法

目前常見的二維譜反演方法主要有奇異值分解法(SVD)、非負約束正則化法(BRD)和截斷奇異值分解法(TSVD)等。

SVD法是通過反演測量誤差的核函數矩陣中的最小奇異值，以實現正則化，求取合理解。顧兆斌等[10]針對SVD法容易導致譜圖不連續的問題，提出一種適用于對信噪比較高的數據進行反演的奇異值分解法。譚茂金等[15]提出了同時考慮TSVD的分辨率和精度的混合反演算法。

BRD法是先將求解方程Tikhonov正則化，再加入罰函數方程的解進行限制，使最終的解在約束條件范圍內。但是BRD法反演結果的優劣取決于選用的正則化因子值是否合適[11]。一般利用L曲線等方法來確定正則化因子，但這些方法在處理低信噪比的數據時常常出現結果沒有意義或在反演時不收斂的問題。所以國內外學者提出了很多解決辦法，如：文獻[16]利用L曲線的斜率來選擇正則化參數；文獻[17]提出一種相對誤差較小的基于SVD和BRD的正則化反演算法等。

2 軟件設計與實現

本文模擬軟件基于.NET4. 6框架，使用C#語言在Visual Studio 2017環境下進行開發。

2.1 結構設計

軟件設計結構如圖2所示，整體分為四大部分：算法、應用模塊、數據及圖形顯示。

圖2 軟件設計結構

軟件設計時采用模塊化編程方式，將所有數學算法函數集中打包放在統一的算法庫中，從而獨立于交互界面，便于算法的持續改進和擴展；應用模塊包括二維譜構造模塊(用于根據用戶設置的流體類型及參數，構造出滿足二維高斯分布的構造譜)，模擬回波采集模塊(用于根據用戶設置的采集模式信息，模擬儀器進行回波采集的過程，同時保存采集的回波串數據)，二維譜反演模塊(實現了SVD和BRD二維譜反演算法，用于將采集得到的回波串數據反演成二維譜)；數據管理模塊和繪圖模塊同樣相對獨立，并設置有外部平臺接口，用于導入實驗室巖心分析數據或實際測井數據。

主界面分為菜單欄、繪圖區、參數區三部分，如圖3所示。菜單欄有文件(NMR File)、參數設置(Parameters Set)和關于(About)。其中：NMR File菜單主要功能包括新建模擬文件、保存模擬文件、打開模擬文件和退出；Parameters Set菜單主要功能包括加載各個應用模塊的參數設置界面、設置參數和導入其他平臺的數據；About菜單用于顯示本文軟件的開發過程及人員信息。

圖3 軟件主界面

2.2 功能實現

(1) 二維譜構造模塊。根據回波串組合序列完成“二維譜構造”模塊，該模塊是基于二維高斯分布進行設計的。通過選擇相應流體類型、填寫流體T2值及其展布、分量孔隙度(各類流體含量)、含氫指數和另一維度參數(縱向弛豫時間T1、擴散系數D、或內部梯度G)，以及最后的布點方案：T2及另一維度參數的布點的起始值、終止值和布點數目，來進行二維譜的構造。

對話框中的各個參數用戶可根據實際模擬的流體進行設置，如圖4所示，完成后軟件會根據用戶填寫的結果自動計算出總孔隙度和核磁孔隙度。參數設置分為流體參數(各類流體的核磁屬性信息及孔隙度分量)和全局參數(模擬模式構造譜分布參數)，所有參數設置完成后，可以實時地顯示出二維構造譜，如圖5所示。

圖4 二維譜構造參數設置對話框

圖5 構造的T2-T1二維分布譜

(2) 模擬回波串采集模塊。該模塊根據上一步設置構造的二維譜數據，結合給定的模擬采集指令模擬采集回波串，并存儲回波數據。參數包括：需要采集的回波組數，各組回波串對應的回波間隔(TE)值、等待時間(TW)值、回波采集個數(NE)、磁場梯度，以及設置的信噪比(S/N)。

參數設置完成后，軟件會自動顯示出模擬采集得到的回波數據及對應的噪音數據，用戶可以根據圖形實時對參數進行修訂，如圖6所示，模擬多TE二維采集得到的回波串數據如圖7所示。另外，該模塊允許用戶導入實驗室巖心分析得到的二維回波數據，以及實際測井得到的二維回波數據，并對數據進行后續的反演處理和分析。

圖6 回波采集參數設置對話框

圖7 回波采集結果

(3) 二維譜反演模塊。該模塊的主要功能是將上一步模擬采集得到的回波串數據進行二維反演，得到最終的二維反演譜。需要設置的參數包括：T2維度布點參數(布點的數量，起始值和終止值)、第二維參數譜布點的方式(布點的數量，起始值和終止值)、反演使用的回波串選擇、平滑因子以及反演方法的選擇，參數設置對話框如圖8所示。圖9為反演得到的T2-T1二維譜。

圖8 二維譜反演參數設置對話框

除了軟件內置的二維反演方法(SVD及BRD)外，還可以利用軟件中提供的C#接口(該接口命名為：INMR_InversionMethod)來編寫自己的反演算法，從而研究不同二維譜反演算法的適用性及影響因素。反演完成后，軟件會再次將反演得到的二維譜進行正演(擬合)，觀測擬合結果與采集得到的回波串是否吻合，從而判斷反演效果好壞。二維回波串擬合效果如圖10所示。

圖10 二維譜回波串擬合效果

3 可靠性驗證及應用實例

3.1 可靠性驗證

軟件編寫完成后，以硫酸銅溶液的核磁物理實驗為例，對正反演模擬軟件的可靠性進行驗證，具體實施步驟如下：

第一步配置三種不同濃度的硫酸銅溶液。

第二步分別取三種不同濃度的硫酸銅溶液3 ml放入玻璃質的集氣瓶中，然后將三只集氣瓶封口后放入試管內進行二維核磁共振(T2-T1)測量。實驗結果顯示：三種濃度的硫酸銅溶液的弛豫時間分別約為2 500、80、2 ms，如圖11所示。

圖11 二維反演譜(物理實驗)

第三步利用本文研發的正反演模擬軟件采用與物理實驗相同的參數進行二維譜構造、回波采集以及二維譜反演，得到構造譜和反演譜的數值模擬結果分別如圖12、圖13所示。

圖12 二維構造譜(軟件模擬)

圖13 二維反演譜(軟件模擬)

將軟件模擬硫酸銅實驗得到的二維構造譜和二維反演譜與實際測量得到的二維反演譜進行對比?？梢钥吹綌抵的M結果與物理實驗結果無論是分布大小還是分布位置都十分接近，驗證了軟件的可靠性。

3.2 應用實例

研發的二維核磁共振正反演模擬軟件可以進行多種二維模式的正反演數值模擬?，F以T2-T1二維模式下核磁共振方法在含有束縛水的油層中識別流體的效果為例進行模擬，驗證軟件的可靠性。

油水的特征參數如表1所示，根據該表參數構造出油水模型的理想分布譜即構造譜，如圖14所示，然后設置不同的等待時間模擬出15組回波串數據：TE均設為1 ms，TW分別為10、30、50、70、100、300、500、700、1 000、2 000、4 000、8 000、10 000、12 000、15 000 ms。

表1 油水模型參數表[14]

圖14 二維構造譜

同時利用軟件的反演模塊將模擬數據進行反演處理，結果如圖15所示。通過對比二維反演譜與構造譜，可見油水的分布位置基本重合,表明了軟件對于T2-T1模式下的油水模型的正反演數值模擬具有較高的可信度。

圖15 二維反演譜

4 結 語

本文研發的二維核磁共振測井正反演模擬軟件可以進行多種二維模式的正反演數值模擬，其模擬結果與實際巖心核磁共振儀器實驗結果非常接近。該軟件可以同時模擬4種不同的核磁共振弛豫過程，并且能夠對各類核磁共振屬性進行單參量因素分析，為二維核磁共振相關研究工作提供了高效的技術手段。