李鵬舉， 王梓丞， 任 莉

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

0 引 言

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

1 T2譜反演方法

1.1 儲(chǔ)層孔隙流體弛豫模型

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

(1)

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

fi——第i個(gè)分量初始幅度；

ti——采樣時(shí)間；

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

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

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

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

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

b=Af，

(2)

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

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

1.2 固定截?cái)嗥娈愔捣?/h3>

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

f=A+b，

(3)

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

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

(4)

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

由于式(4)中奇異值由大到小排列，將小于預(yù)設(shè)截?cái)嘀郸膋的項(xiàng)全部取零值得

(5)

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

δk——截?cái)嗥娈愔担?/p>

g——無噪聲信號(hào)；

ε——噪聲信號(hào)。

在此步驟中難點(diǎn)在于截?cái)嘀档倪x取工作，若截?cái)嗥娈愔郸膋選取偏大會(huì)導(dǎo)致結(jié)果誤差變大，降低準(zhǔn)確性；相反則導(dǎo)致解的穩(wěn)定性變差。適當(dāng)?shù)慕財(cái)嘀当仨毱胶饨獾姆€(wěn)定性與結(jié)果的平方誤差兩方面影響。

1.3 自適應(yīng)截?cái)嗥娈愔捣?/h3>

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

(6)

式中：fs——有用信號(hào)項(xiàng)；

fn——噪聲項(xiàng)。

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

(7)

利用T2snr評(píng)價(jià)解的準(zhǔn)確性，T2snr值與解的準(zhǔn)確性呈正相關(guān)關(guān)系，與解的穩(wěn)定性呈負(fù)相關(guān)。利用自適應(yīng)搜索法自動(dòng)選取奇異值截?cái)嘀担胶饨獾姆€(wěn)定性與準(zhǔn)確性相矛盾現(xiàn)象。

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

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

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

步驟二以第k個(gè)奇異值作為奇異值截?cái)嘀颠M(jìn)行反演計(jì)算，得到f。

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

(8)

式中,Vi、Ui——標(biāo)準(zhǔn)正交向量組。

噪聲項(xiàng)的方差為

(9)

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

步驟六以δk-1作為奇異值截?cái)嘀递敵觥?/p>

1.4 有效集非負(fù)約束法

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

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

(10)

(11)

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

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

(12)

式中,H——n階對(duì)角陣。

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

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

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

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

(13)

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

λk=Bkrk，

(14)

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

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

2 仿真實(shí)驗(yàn)

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

2.1 不同截?cái)嗥娈愔捣椒ǖ膶?duì)比

正演過程中，分別給回波信號(hào)加入信噪比αSNR分別為10、50、200、105的高斯噪聲，分別用自適應(yīng)截?cái)嗥娈愔岛凸潭ń財(cái)嗥娈愔捣椒ㄟM(jìn)行反演，固定截?cái)嗥娈愔郸姆謩e為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)造譜。以某一固定截?cái)嗥娈愔档玫降腡2譜曲線在不同信噪比條件下曲線差異巨大，并且不能保證準(zhǔn)確性，這說明傳統(tǒng)固定截?cái)嗥娈愔捣椒ㄔ谠肼暩蓴_的影響下穩(wěn)定性較差。自適應(yīng)截?cái)嗥娈愔捣椒軌蚋鶕?jù)噪聲干擾的影響自動(dòng)選取截?cái)嘀担サ粜〉姆橇闫娈愔狄栽鰪?qiáng)穩(wěn)定性，不會(huì)因噪聲干擾程度不同對(duì)結(jié)果產(chǎn)生較大影響，抗噪能力強(qiáng)，并且峰值大小與峰值所對(duì)應(yīng)的弛豫時(shí)間較為準(zhǔn)確。

2.2 非負(fù)約束方法

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

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

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

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

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

3 實(shí) 例

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

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

圖5中第一道為深度道；第二道為自然伽馬曲線；第三道SPECTA為利用自適應(yīng)截?cái)嗥娈愔到Y(jié)合有效集非負(fù)約束法反演得到的T2譜；第四道DISTA為CIFLog測(cè)井處理解釋平臺(tái)反演得到的T2譜。

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

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

4 結(jié) 論

(1)自適應(yīng)截?cái)嗥娈愔捣椒軌蛟诒ＷCT2譜穩(wěn)定的前提下，最大限度提取回波串中的有用信息。在噪聲水平變化時(shí)，能夠根據(jù)噪聲自動(dòng)調(diào)節(jié)奇異值個(gè)數(shù)，使反演T2譜能夠在平均值意義下達(dá)到較高的精度。

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

(3)自適應(yīng)截?cái)嗥娈愔到Y(jié)合有效集非負(fù)約束方法能夠均衡反演結(jié)果的準(zhǔn)確性與穩(wěn)定性，較為精準(zhǔn)地反映各弛豫時(shí)間信號(hào)大小。該方法具有準(zhǔn)確度高、穩(wěn)定性好、抗噪能力較強(qiáng)和反演速度快的特點(diǎn)，能夠滿足野外核磁共振測(cè)井資料處理的要求。