一種通用的巖心核磁共振實(shí)驗(yàn)數(shù)據(jù)分析軟件設(shè)計(jì)與實(shí)現(xiàn)

覃瑩瑤 張 宮 何宗斌 張家成

(長江大學(xué)油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室 湖北 武漢 430100)

0 引 言

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

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

1 原 理

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

圖1 軟件業(yè)務(wù)流程設(shè)計(jì)

1.1 數(shù)據(jù)解析

軟件需要加載各類巖心分析儀測量得到不同格式的原始回波串?dāng)?shù)據(jù)，首先要對數(shù)據(jù)進(jìn)行解析。以CoreSpec-1000型儀器測量得到的RES數(shù)據(jù)和NIUMAG公司儀器測量得到的PEA數(shù)據(jù)為例。首先進(jìn)行RES文件或PEA、PAR文件的格式分析，如圖2所示RES文件中有記錄參數(shù)信息的頭文件和原始回波數(shù)據(jù)；PEA文件中是原始回波數(shù)據(jù)，參數(shù)信息在PAR文件中。然后對文件進(jìn)行解析，讀取RES文件或PEA、PAR文件的數(shù)據(jù)后，利用C#面向?qū)ο蟪绦蛟O(shè)計(jì)“CPMGCoreData”類提供的一系列方法將數(shù)據(jù)都解析成統(tǒng)一的模式。

圖2 RES文件數(shù)據(jù)

1.2 反演方法

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

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

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

其矩陣方程的形式為:

Y=A·P

(1)

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

ΔY=Y-Y′

(2)

(3)

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

(4)

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

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

(5)

SIRT算法簡單容易實(shí)現(xiàn)，在利用算法處理時不需要用戶干預(yù)，也不需要設(shè)置一些復(fù)雜的反演控制參數(shù)，減少了人為因素方面的反演結(jié)果誤差。

1.3 參數(shù)計(jì)算

(1) 核磁孔隙度的確定 用標(biāo)準(zhǔn)樣品對飽和巖樣測得的T2譜進(jìn)行刻度，將核磁信號強(qiáng)度轉(zhuǎn)換成孔隙度[10]，轉(zhuǎn)換公式如下：

(6)

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

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

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

(7)

(8)

式中：φi為對應(yīng)分量T2i的孔隙度分量。

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

1.4 數(shù)據(jù)報(bào)告導(dǎo)出

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

圖3 巖心數(shù)據(jù)報(bào)告

2 軟件設(shè)計(jì)及實(shí)現(xiàn)

根據(jù)原理和軟件業(yè)務(wù)流程設(shè)計(jì)，本軟件基于.NET Framework 4.6開發(fā)框架在Visual Studio 2017上使用C#語言進(jìn)行開發(fā)。

2.1 軟件架構(gòu)

軟件架構(gòu)設(shè)計(jì)結(jié)構(gòu)如圖4所示，整體分為三大部分：算法、應(yīng)用模塊、數(shù)據(jù)及圖形顯示。

圖4 軟件設(shè)計(jì)結(jié)構(gòu)

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

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

圖5 軟件主界面

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

2.2 功能實(shí)現(xiàn)

2.2.1巖心分析數(shù)據(jù)加載模塊

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

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

圖6 數(shù)據(jù)加載模塊界面

圖7 原始回波串顯示

2.2.2反演模塊

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

圖8 反演模塊及反演結(jié)果

2.2.3參數(shù)計(jì)算模塊

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

圖9 參數(shù)計(jì)算模塊及計(jì)算結(jié)果

3 可靠性驗(yàn)證及應(yīng)用實(shí)例

3.1 可靠性驗(yàn)證

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

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

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

本軟件和NIUMAG公司生產(chǎn)的儀器配套的軟件進(jìn)行比對。用NIUMAG公司生產(chǎn)的儀器測量一塊飽和水砂巖巖心，圖11是NIUMAG儀器配套軟件反演的T2譜和用本軟件進(jìn)行反演處理得到的譜?？梢杂^察到兩個T2譜的形態(tài)相似，峰值位置也非常接近。

這兩個驗(yàn)證說明本文研發(fā)的軟件和其他軟件處理結(jié)果一致，驗(yàn)證了本軟件的可靠性。

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

3.2 應(yīng)用實(shí)例

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

第一步：用CoreSpec-1000型儀器和NIUMAG儀器對一塊飽和水的砂巖巖心進(jìn)行核磁共振實(shí)驗(yàn)，得到兩組巖心數(shù)據(jù)。

第二步：將兩儀器的兩組巖心數(shù)據(jù)導(dǎo)出，用研發(fā)的巖心核磁共振實(shí)驗(yàn)數(shù)據(jù)分析軟件處理這兩組數(shù)據(jù)，進(jìn)行反演得到T2譜。

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

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

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

圖13 多回波間隔數(shù)據(jù)處理結(jié)果

4 結(jié) 語

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