李鵬舉， 王梓丞， 任 莉

(1.東北石油大學(xué) 地球科學(xué)學(xué)院， 黑龍江 大慶163318； 2.東北石油大學(xué) 非常規(guī)油氣成藏與開發(fā)省部共建國家重點實驗室培育基地， 黑龍江 大慶163318； 3.中國石油集團測井有限公司大慶分公司解釋評價中心， 黑龍江 大慶163000)

0 引 言

核磁共振測井技術(shù)能夠提供地層巖石物理性質(zhì)和流體性質(zhì)參數(shù)，如孔隙度、滲透率、可動水飽和度等[1,2]。這些參數(shù)能否計算準確與T2解譜技術(shù)息息相關(guān)。解譜技術(shù)是核磁共振測井技術(shù)的關(guān)鍵。目前，廣泛使用的解譜技術(shù)的實質(zhì)是多指數(shù)反演，反演方法有奇異值分解(SVD)[1,3-5]和最優(yōu)化及其改進方法[6]。例如：肖立志等[1,5]基于奇異值分解的反演方法；王為民等[7]提出基于解的自由度選取最優(yōu)反演模型的反演變換法；K.J.Kunn等[8]，王才志等[9]采用奇異值分解與正則化方法結(jié)合截去小的非零奇異值的方法。學(xué)者們普遍采用截斷法或正則化方法確保解的穩(wěn)定與平滑。截斷奇異值的選擇和正則化因子的選取與噪聲有關(guān)，一般采用信噪比作為噪聲的統(tǒng)計特征參數(shù)。而CPMG回波信號為多指數(shù)衰減信號，與穩(wěn)恒信號不同，信噪比并不是常數(shù)[3]。因此，以信噪比作為衰減信號的統(tǒng)計特征參數(shù)不適合。以噪聲的方差作為噪聲的統(tǒng)計特征參數(shù)，來確定截止奇異值，能夠最大限度提取自旋回波中的有用信息。根據(jù)解的物理意義，T2譜的幅值應(yīng)為非負實數(shù)。為了保證解的非負性，廣泛采用改進迭代法隨機搜索非負解，而這種解法可能會喪失解的連續(xù)性[4,10]。產(chǎn)生這個問題的原因是方程的求解與非負約束相分離。為此，筆者結(jié)合奇異值分解和最優(yōu)化方法，在保證解的穩(wěn)定性的同時進行非負約束，避開對解的隨機迭代搜索，以有限的計算次數(shù)快速求解連續(xù)的T2譜。

1 T2譜反演方法

1.1 儲層孔隙流體弛豫模型

巖石的孔隙度由巖石的礦物組分、孔隙中黏土分布與孔隙結(jié)構(gòu)等因素共同決定。核磁共振測井方法能夠提供巖石中各種流體的信息以及滲透率信息，并且不受巖性影響。測井儀器以氫核為目標采集包含孔隙流體信息的自旋回波串，通過自旋回波法得到橫向弛豫時間常數(shù)[1-2,11]，即

(1)

式中：bi——回波幅度序列；

fi——第i個分量初始幅度；

ti——采樣時間；

T2j——第j個橫向弛豫時間常數(shù)；

n——回波個數(shù)；

m——橫向弛豫分量個數(shù)；

εi——噪聲序列[4]。

將式(1)寫成矩陣形式，即

b=Af，

(2)

其中，b=g+ε；g∈Rn；b∈Rn；A∈Rn×m；f∈Rm；ε∈Rn。

由初始幅度f預(yù)設(shè)多個橫向弛豫時間常數(shù)擬合多指數(shù)衰減序列b的過程稱為T2譜反演。

1.2 固定截斷奇異值法

式(2)不存在唯一解，無法進行工程應(yīng)用，為得到最合適的最小二乘解，由式(2)得

f=A+b，

(3)

式中，A+——矩陣A的偽逆矩陣。

將矩陣A+奇異值分解代入式(3)，得

(4)

其中，r=rank(A)，Ur∈Rn×r，Vr∈Rm×r。

由于式(4)中奇異值由大到小排列，將小于預(yù)設(shè)截斷值δk的項全部取零值得

(5)

式中：δ1——最大奇異值；

δk——截斷奇異值；

g——無噪聲信號；

ε——噪聲信號。

在此步驟中難點在于截斷值的選取工作，若截斷奇異值δk選取偏大會導(dǎo)致結(jié)果誤差變大，降低準確性；相反則導(dǎo)致解的穩(wěn)定性變差。適當?shù)慕財嘀当仨毱胶饨獾姆€(wěn)定性與結(jié)果的平方誤差兩方面影響。

1.3 自適應(yīng)截斷奇異值法

傳統(tǒng)截斷奇異值方法衡量噪聲統(tǒng)計特征時，以信噪比作為參數(shù)[1-2]，并以信噪比作為截斷奇異值選取依據(jù)。然而，回波序列是一種不同于穩(wěn)恒信號的多指數(shù)衰減信號，在噪聲環(huán)境不變的情況下，信噪比是隨時間變化的變量。在實際處理過程中影響信噪比的因素很多，單純以信噪比作為截斷奇異值選取依據(jù)，無法保證解的穩(wěn)定性和一致性。本文提出的方法在衡量噪聲統(tǒng)計特征時，以噪聲的方差作為參數(shù)衡量噪聲幅度。由式(5)得：

(6)

式中：fs——有用信號項；

fn——噪聲項。

引入?yún)?shù)T2snr，有

(7)

利用T2snr評價解的準確性，T2snr值與解的準確性呈正相關(guān)關(guān)系，與解的穩(wěn)定性呈負相關(guān)。利用自適應(yīng)搜索法自動選取奇異值截斷值，平衡解的穩(wěn)定性與準確性相矛盾現(xiàn)象。

首先對經(jīng)過奇異值分解得到的對角陣中奇異值最大值進行搜索，判斷由經(jīng)驗選取的目標函數(shù)M是否小于T2snr。若不滿足，計算下一個奇異值；若滿足，則停止搜索，并以上一個奇異值為奇異值截斷值。自適應(yīng)搜索法具體流程圖如圖1所示，步驟如下：

圖1 自適應(yīng)搜索流程Fig. 1 Adaptive search flow

步驟一以對角陣中奇異值最大值δ1作為搜索初始值，即令k=1。

步驟二以第k個奇異值作為奇異值截斷值進行反演計算，得到f。

步驟三計算fs和|fn|。根據(jù)研究可知，反演結(jié)果中噪聲項均值為零，因此，fs=f。隨機信號ε的方差由正交方程組性質(zhì)得：

(8)

式中,Vi、Ui——標準正交向量組。

噪聲項的方差為

(9)

步驟五令k=k+1，執(zhí)行步驟二。

步驟六以δk-1作為奇異值截斷值輸出。

1.4 有效集非負約束法

在井下核磁儀器采集信號過程中，難免會記錄無用信號，例如噪聲。在解譜過程中由于無用信號的存在，反演結(jié)果無法準確反映井下真實流體信息，尤其是解出現(xiàn)負值的情況，不符合實際物理意義。傳統(tǒng)的非負約束方法，例如迭代法非負約束等均采用迭代的方法去除負解，其缺點有：迭代次數(shù)較多時，計算速度較慢；迭代次數(shù)較少時約束效果較差。采用有效集最小二乘非負約束方法可保證解的非負性。

根據(jù)最優(yōu)化理論的思想結(jié)合T2譜反演構(gòu)造目標函數(shù)，以解的非負性作為約束條件，得

(10)

(11)

由二次規(guī)劃一般形式得

此外，我國養(yǎng)老服務(wù)還存在發(fā)展不平衡問題，居家社區(qū)與機構(gòu)養(yǎng)老服務(wù)發(fā)展不均衡，城鄉(xiāng)養(yǎng)老服務(wù)業(yè)發(fā)展不平衡，區(qū)域間養(yǎng)老服務(wù)業(yè)發(fā)展不平衡。公辦機構(gòu)建設(shè)水平總體好于民辦機構(gòu)，城市養(yǎng)老機構(gòu)明顯高于農(nóng)村敬老院的保障水平。

(12)

式中,H——n階對角陣。

有效集最小二乘非負約束方法獲取最優(yōu)解α具體流程如圖2所示，步驟如下：

圖2 有效集非負約束流程Fig. 2 Active set non-negativity constraint flow

步驟一確定初始值。以可行點αk(k=0)為起點。

步驟二求解子問題。確定有效約束的等式條件，得

(13)

步驟三終止準則的檢驗。引入新的參數(shù)λk，即拉格朗日乘子，表達式為

λk=Bkrk，

(14)

令(λk)t=min{(λk)i}，i∈I(αk)。若(λk)t≥0，終止運算，并以αk為結(jié)果輸出；若(λk)t<0，則將有效集中第t個約束條件去除，令Sk=Sk/{t}，執(zhí)行步驟二。

步驟五令k=k+1，執(zhí)行步驟二。

2 仿真實驗

為了檢驗方法的有效性，參考實際地質(zhì)狀況，運用數(shù)值模擬技術(shù)構(gòu)造T2譜得到構(gòu)造譜，參照P型核磁系統(tǒng)，按照式(1)正演出不同信噪比的回波數(shù)據(jù)。對這些數(shù)據(jù)做T2譜反演，與構(gòu)造譜對比，研究方法的有效性。

2.1 不同截斷奇異值方法的對比

正演過程中，分別給回波信號加入信噪比αSNR分別為10、50、200、105的高斯噪聲，分別用自適應(yīng)截斷奇異值和固定截斷奇異值方法進行反演，固定截斷奇異值δ分別為1.389 2、0.313 1、0.064 7，如圖3所示。

圖3 不同信噪比下T2譜反演結(jié)果Fig. 3 Comparison of T2 spectrum inversion results with different αSNR

從圖3中可以看出，自適應(yīng)譜隨著信噪比的變化分別與構(gòu)造譜、δ1、δ2、δ3譜重合，并且自適應(yīng)譜總能最接近構(gòu)造譜。以某一固定截斷奇異值得到的T2譜曲線在不同信噪比條件下曲線差異巨大，并且不能保證準確性，這說明傳統(tǒng)固定截斷奇異值方法在噪聲干擾的影響下穩(wěn)定性較差。自適應(yīng)截斷奇異值方法能夠根據(jù)噪聲干擾的影響自動選取截斷值，去掉小的非零奇異值以增強穩(wěn)定性，不會因噪聲干擾程度不同對結(jié)果產(chǎn)生較大影響，抗噪能力強，并且峰值大小與峰值所對應(yīng)的弛豫時間較為準確。

2.2 非負約束方法

對比在信噪比αSNR分別為10、100時，非負最小二乘約束、簡單迭代法非負約束(迭代次數(shù)n=1000)和正則化因子非負約束下的反演效果，如圖4所示。

圖4 不同信噪比下非負約束對比 Fig. 4 Non-negative constraint comparison of different αSNR

與構(gòu)造譜對比，在迭代次數(shù)較少時，簡單的迭代法非負約束和正則化因子非負約束在T2譜尾端均出現(xiàn)幅度值小于零值的情況，并且隨信噪比的降低，迭代法非負約束出現(xiàn)負值幅度變大的情況。有效集非負最小二乘法非負約束T2譜尾端未出現(xiàn)負幅度值，且曲線平滑，與構(gòu)造譜曲線重合度較高。

在迭代次數(shù)較少時，簡單迭代方法無法完全實現(xiàn)非負約束，正則化因子法為了實現(xiàn)非負約束，犧牲了解的準確性。非負最小二乘法在進一步保證解的準確性的同時實現(xiàn)了完全的非負約束。

簡單迭代法和正則化法要實現(xiàn)完全非負約束需要犧牲大量的計算時間。而由有效集非負最小二乘法步驟可知，該方法無需迭代，有計算速度快的優(yōu)點。

3 實 例

某井在2 170～2 230 m井段核磁共振測井，利用新舊兩種T2譜反演方法得出的結(jié)果對比如圖5所示。

圖5 不同方法下T2譜反演結(jié)果對比Fig. 5 Comparison of T2 spectrum inversion results of different methods

圖5中第一道為深度道；第二道為自然伽馬曲線；第三道SPECTA為利用自適應(yīng)截斷奇異值結(jié)合有效集非負約束法反演得到的T2譜；第四道DISTA為CIFLog測井處理解釋平臺反演得到的T2譜。

由圖5可以看出，SPECTA和DISTA反演結(jié)果整體基本一致。在深度2 170～2 182 m井段，SPECTA在弛豫時間30～3 000 ms處相比DISTA能夠清晰地反映出完整的峰值形態(tài)。在深度2 200～2 210 m井段，SPECTA和DISTA兩者形態(tài)基本相同，但SPECTA較為平緩，未出現(xiàn)較大幅度波動。深度2 217～2 226 m井段為泥巖段，自然伽馬曲線幅度較高，SPECTA和DISTA兩者在低弛豫時間處均出現(xiàn)明顯幅度。

利用自適應(yīng)截斷奇異值結(jié)合有效集非負約束法反演出的T2譜具有精細的結(jié)構(gòu)，曲線平滑且連續(xù)。該方法在處理速度上也有較好的表現(xiàn)，能夠在幾秒鐘內(nèi)完成對整口井的T2譜反演處理工作。

4 結(jié) 論

(1)自適應(yīng)截斷奇異值方法能夠在保證T2譜穩(wěn)定的前提下，最大限度提取回波串中的有用信息。在噪聲水平變化時，能夠根據(jù)噪聲自動調(diào)節(jié)奇異值個數(shù)，使反演T2譜能夠在平均值意義下達到較高的精度。

(2)有效集非負約束方法能以有限的運算次數(shù)，快速實現(xiàn)解的非負約束，得到連續(xù)的非負解，該方法相比迭代非負約束方法和正則化因子非負約束方法準確性高，且迭代次數(shù)少，運算速度快。

(3)自適應(yīng)截斷奇異值結(jié)合有效集非負約束方法能夠均衡反演結(jié)果的準確性與穩(wěn)定性，較為精準地反映各弛豫時間信號大小。該方法具有準確度高、穩(wěn)定性好、抗噪能力較強和反演速度快的特點，能夠滿足野外核磁共振測井資料處理的要求。