一種通用的巖心核磁共振實驗數據分析軟件設計與實現

覃瑩瑤 張 宮 何宗斌 張家成

(長江大學油氣資源與勘探技術教育部重點實驗室 湖北 武漢 430100)

0 引 言

核磁共振測井是一種有效的測井方法，巖心核磁共振實驗分析是NMR測井資料采集、處理、解釋及應用的基礎[1]。測井前做少量的實驗室測量能節省實際測井時間，提高測井曲線質量。而且實驗室的巖心核磁共振測量實驗能夠實現巖心與測井之間的標定,建立核磁共振測井解釋模型[2]。利用巖心核磁共振實驗可以直接觀測到巖石孔隙的流體信號，得到與巖性無關的孔隙度、孔徑分布、滲透率、自由流體指數及流體飽和度等參數[3-4]。巖心核磁共振實驗中，原始采集數據的處理與分析非常關鍵，但各廠家核磁共振儀器提供的數據處理配套軟件的反演方法、參數均不相同。例如美國NUMAR公司生產的CoreSpec-1000型巖心核磁共振分析儀由Echofit軟件利用非負最小二乘法(NNLS)和基于奇異值分解的MAP-Ⅱ算法進行反演[5]；英國Resonance Instruments公司的MARAN系列核磁共振巖心分析儀沒有提供配套的數據處理分析軟件[6]；NIUMAG公司生產的巖心核磁共振分析儀配套軟件使用罰函數法(BRD)和聯合重建迭代法(SIRT)進行反演[7]。儀器性能和反演算法都不相同導致同樣的巖心測量得到的回波信息無法對比分析，無法確定反演的T2譜有差異的原因，很大程度上限制了核磁巖心分析儀的作用。王忠東等研發的巖心核磁實驗解釋分析軟件CoreMR可以對MARAN系列譜儀所采集巖心和流體核磁信號進行解釋分析[6]，但是這款軟件只能處理一種譜儀采集的數據，沒有實現對不同廠家儀器測量數據的統一處理分析。

研發一種通用的巖心核磁共振實驗數據分析軟件可以解決這一問題，使用C#語言進行軟件開發，支持硬件加速的WPF通用框架進行繪圖，采用格式統一、易于交互的XML序列化存儲反演參數，易于傳輸、解析的JSON序列化存儲采樣數據和處理結果。通過采用統一的模式進行數據存儲、統一的快速高分辨反演算法，設置不同參數進行核磁參數計算，本軟件可以使所有儀器測量的回波數據得到統一處理和分析，在消除軟件處理方法差異的影響上使同樣巖心的回波信息進行對比分析。

1 原 理

巖心核磁共振實驗數據分析軟件一般包括加載數據、反演、計算參數和導出報告四個步驟。本文設計的軟件力圖將不同儀器所采集的數據進行統一處理，業務流程如圖1所示，具體步驟為：① 數據加載要實現把所有儀器所測的不同格式數據解析成統一的模式；② 進行T2譜反演；③ 計算參數要根據反演后的T2譜計算出孔隙度、束縛水飽和度等核磁參數；④ 數據的處理和分析完成后需要導出指定格式的實驗報告。

圖1 軟件業務流程設計

1.1 數據解析

軟件需要加載各類巖心分析儀測量得到不同格式的原始回波串數據，首先要對數據進行解析。以CoreSpec-1000型儀器測量得到的RES數據和NIUMAG公司儀器測量得到的PEA數據為例。首先進行RES文件或PEA、PAR文件的格式分析，如圖2所示RES文件中有記錄參數信息的頭文件和原始回波數據；PEA文件中是原始回波數據，參數信息在PAR文件中。然后對文件進行解析，讀取RES文件或PEA、PAR文件的數據后，利用C#面向對象程序設計“CPMGCoreData”類提供的一系列方法將數據都解析成統一的模式。

圖2 RES文件數據

1.2 反演方法

對原始回波信號進行反演處理得到T2分布譜，這個過程叫做反演[8]。T2譜的反演方法較多。其中，最常用的算法有：罰函數法(BRD)、非負最小二乘法(NNLS)、聯合迭代重建算法(SIRT)、奇異值分解法(SVD)以及針對SVD法的各種改進型反演算法。

軟件的反演模塊主要實現了罰函數法(BRD)和聯合迭代重建算法(SIRT)改進的快速高分辨反演算法，可以在不抽樣壓縮回波串的情況下，實時得到T2譜。

在巖心核磁共振實驗中，核磁共振橫向弛豫信號y(t)描述為如下的形式：

其矩陣方程的形式為:

Y=A·P

(1)

Y為測量的橫向弛豫信號；A為n×m系數矩陣；P為所要反演計算的T2譜中各個弛豫時間所對應的譜幅度值。在此利用聯合迭代重建算法(SIRT)實現反演[9]，首先給定譜的初始模型P′，求出預測弛豫信號Y′與實測弛豫信號Y的誤差ΔY，即:

ΔY=Y-Y′

(2)

(3)

利用Δyi求出弛豫譜幅度的修正量Δpj，這是SIRT算法實現的主要思想。pj的修正量為:

(4)

通過式(4)計算的修正公式為:

P(q+1)=P(q)+ΔP

(5)

SIRT算法簡單容易實現，在利用算法處理時不需要用戶干預，也不需要設置一些復雜的反演控制參數，減少了人為因素方面的反演結果誤差。

1.3 參數計算

(1) 核磁孔隙度的確定 用標準樣品對飽和巖樣測得的T2譜進行刻度，將核磁信號強度轉換成孔隙度[10]，轉換公式如下：

(6)

式中：φnmr為巖樣核磁孔隙度值(%)；M為標準樣品T2譜的總幅度；V為標準樣品總含水量(cm3)；S為標準樣品在核磁共振數據采集時的累積次數；G為標準樣品在核磁共振數據采集時的接受增益；mi為第i個巖樣T2分量的核磁共振T2譜幅度；v為巖樣的體積(cm3)；s為巖樣在核磁共振數據采集時的累積次數；g為樣品在核磁共振數據采集時的接受增益。

(2)T2截止值的確定 采用孔隙度累加法確定T2截止值[11]。先將巖樣用水100%飽和得到飽水樣，進行核磁共振測量得到T2分布、累加孔隙度曲線和有效孔隙度值(用MPHI表示)。然后對巖樣做脫水處理，在一定壓力的條件下，使自由水脫出巖樣，得到孔隙空間中只剩下束縛水的離心樣，再做巖心核磁共振測量，測得T2分布、累加孔隙度曲線和束縛水孔隙體積值(用MBVI表示)。根據測量結果，以MBVI作一條與橫軸平行的平行線，此線與100%飽和水時的累加孔隙度曲線有一交點，以該交點作一條與縱軸平行的平行線，此線與橫軸交點對應的T2值即T2cutoff值。

(3)T2平均值的計算 在分析過程中，常用T2的平均值來表征T2分布，算術平均值T2s和幾何平均值T2g按下式計算[12]：

(7)

(8)

式中：φi為對應分量T2i的孔隙度分量。

(4) 核磁束縛水飽和度的確定 確定核磁共振束縛水飽和度的方法是以小孔隙束縛水體積模型作為理論基礎的T2截止值法[13]。在巖心核磁共振實驗測得的T2譜曲線中，核磁束縛水飽和度為T2譜中小于T2截止值T2cutoff的不可動峰下包面積和整個T2譜曲線下包面積的比值。束縛水體積等于核磁束縛水飽和度和孔隙度的乘積，可動水體積等于巖樣的孔隙體積與束縛水體積的差值。

1.4 數據報告導出

利用軟件對巖心原始回波串數據設置合適參數進行處理和分析后，可以將文件保存成易于傳輸和解析的JSON格式文件，還能夠將分析結果導出成所選擇的固定格式的EXCLE格式報告，如圖3所示。巖心測量報告中有具體的實驗結果和回波數據，實驗結果中包括了巖心的井號、樣號、體積、回波間隔等參數信息，孔隙度、飽和度等核磁參數信息，回波衰減曲線和反演的T2分布譜；回波數據中包括了采集時間、原始回波和T2分布的時間及幅度。

圖3 巖心數據報告

2 軟件設計及實現

根據原理和軟件業務流程設計，本軟件基于.NET Framework 4.6開發框架在Visual Studio 2017上使用C#語言進行開發。

2.1 軟件架構

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

圖4 軟件設計結構

軟件使用的所有算法統一放在算法庫中，獨立于交互界面。本文軟件實現了罰函數法和聯合迭代重建算法改進的反演算法，利用外部算法接口，可以對算法進行持續改進和擴展；應用模塊包括數據加載模塊、反演模塊和參數計算模塊；繪圖模塊和數據管理模塊同樣相對獨立，并專門設置有外部平臺接口，可以很容易地將研究成果遷移到其他平臺進行使用。其中最核心的三個模塊設計了四個類：CPMGCoreData類實現了數據的輸入與輸出；T2Inversion類實現反演的功能；Parameter類實現參數計算，這三個類都引用了通用算法庫的CommonAnalysis類。

軟件主界面如圖5所示，包括菜單欄、繪圖區、參數查閱區和狀態欄四部分。菜單欄有文件(File)、工具(Tool)、數據分析(Data Analysis)和幫助(Help)。

圖5 軟件主界面

File菜單主要是加載實驗數據、輸出報告、打開文件和保存文件；Tool菜單中有不同儀器的工具箱，功能主要有批量實驗數據處理、導出測井曲線文件；Data Analysis菜單主要實現反演和參數計算功能；幫助菜單中是本軟件的一些開發過程及人員信息。

2.2 功能實現

2.2.1巖心分析數據加載模塊

本軟件將各類巖心分析儀測量得到的原始回波串數據經過CPMGCoreData類提供的一系列方法將數據都解析后，用統一的模式進行存儲。軟件目前能實現的主要包括：(1) CoreSpec-1000型巖心核磁共振分析儀測量得到的RES數據；(2) NIUMAG公司生產的各型巖心核磁共振分析儀測量得到的PEA數據；(3) OXFORD公司生產的各型巖心核磁共振分析儀測量得到的RIDAT數據；(4) Magritek公司生產的儀器測量得到的數據；(5) EXCEL通用數據格式。

軟件自動解析并加載原始回波數據和對應的參數信息，如圖6所示。通過參數表格的形式，顯示出了巖心分析儀采集數據的參數及采集得到的具體回波串數據，并在繪圖區域繪制原始回波串(圖7)。參數表列中主要參數有：等待時間、回波間隔、采集數目等。

圖6 數據加載模塊界面

圖7 原始回波串顯示

2.2.2反演模塊

反演模塊利用罰函數法(BRD)和聯合迭代重建算法(SIRT)改進的快速高分辨反演算法對加載的回波串數據進行反演。核心算法是原理中式(1)-式(5)實現的聯合迭代重建算法。反演參數設置的參數表區域中，可以直接在里面進行反演參數的修改。主要參數包括：是否需要進行基線漂移校正、反演平滑級別、起始反演回波、終止反演回波、是否進行邊界約束、反演布點數目、布點最小值、布點最大值。在反演參數設置區域填寫合適的反演參數后，單擊右下角反演按鈕，即可進行快速高分辨反演，反演結果(圖8)會顯示在左下角的繪圖區域。反演的同時，會對原始回波串進行擬合觀察反演結果的好壞。

圖8 反演模塊及反演結果

2.2.3參數計算模塊

參數計算模塊采用了參數計算原理中描述的確定孔隙度、T2截止值等核磁參數的公式和方法進行計算。圖9所示的參數計算設置區域中，主要需要設置的參數包括：樣品體積，T2截止值，滲透率計算模型參數等。反演完成后，在參數設置區域填寫正確的巖心分析參數，然后單擊分析按鈕即可完成分析，其中最重要的參數是巖心的體積。參數計算能夠得到T2幾何均值、有效孔隙度、可動孔隙度、滲透率、束縛水飽和度等參數。

圖9 參數計算模塊及計算結果

3 可靠性驗證及應用實例

3.1 可靠性驗證

為了驗證軟件的可靠性，將本文研發的軟件和其他軟件對回波信息處理的結果進行比對，以CoreSpec-1000型儀器和NIUMAG公司儀器的配套軟件為例。

本軟件和CoreSpec-1000型儀器配套的軟件進行比對。用CoreSpec-1000型儀器測量硫酸銅標樣，圖10為CoreSpec-1000型儀器的配套軟件處理的T2譜和本軟件對原始回波數據進行反演處理得到的T2譜，可以看出峰值位置非常接近，都在200 ms附近。

圖10 CoreSpec-1000配套軟件和本軟件反演的T2譜比對

本軟件和NIUMAG公司生產的儀器配套的軟件進行比對。用NIUMAG公司生產的儀器測量一塊飽和水砂巖巖心，圖11是NIUMAG儀器配套軟件反演的T2譜和用本軟件進行反演處理得到的譜。可以觀察到兩個T2譜的形態相似，峰值位置也非常接近。

這兩個驗證說明本文研發的軟件和其他軟件處理結果一致，驗證了本軟件的可靠性。

圖11 NIUMAG儀器配套軟件和本軟件反演的T2譜比對

3.2 應用實例

研發的巖心核磁共振實驗數據分析軟件可以使不同儀器測量同樣巖心得到的回波信息進行對比分析。以對一塊飽和水的砂巖巖心進行不同儀器的核磁共振比對實驗為例，具體步驟為：

第一步：用CoreSpec-1000型儀器和NIUMAG儀器對一塊飽和水的砂巖巖心進行核磁共振實驗，得到兩組巖心數據。

第二步：將兩儀器的兩組巖心數據導出，用研發的巖心核磁共振實驗數據分析軟件處理這兩組數據，進行反演得到T2譜。

圖12是本文研發的軟件對NIUMAG和CoreSpec-1000兩種儀器測得同樣巖心兩組回波數據進行反演得到的T2譜。在利用本軟件對兩組數據處理的反演方法和參數相同的情況下，觀察到譜的形態有差異，峰值位置不同。圖12(a)中主峰峰值位置在100 ms附近，圖12(b)主峰峰值位置在200 ms附近，是由于圖12(a)的儀器主頻較高，圖12(b)的儀器主頻較低，反映出同樣巖心T2譜的不同是儀器性能不同所導致的。說明了本軟件用統一的反演算法處理巖心數據，能對比同樣巖心的回波信息反映出儀器性能的不同。

圖12 軟件對NIUMAG儀器和CoreSpec-1000型儀器

軟件可以一次處理多回波間隔的數據得到長短TE的核磁T2譜。以CoreSpec-1000型儀器所測量的兩個回波間隔分別為0.1 ms和0.2 ms的數據為例，用本軟件一次進行處理，得到了如圖13所示的回波串和長短TE的T2譜。

圖13 多回波間隔數據處理結果

4 結 語

本文開發的軟件處理巖心核磁共振實驗數據結果精確，操作簡單方便，且具有通用性，實現了對所有儀器測量的回波數據進行統一處理和分析，可以對比分析不同儀器測量同樣巖心的回波信息，反映出儀器性能的差異。另外軟件預留有外部接口，方便擴展，其中其他系統接口可以對接外部平臺，外部算法接口可以迅速實現新的算法。