王鵬，譚茂金，2

（1.中國地質(zhì)大學(xué)地球物理與信息技術(shù)學(xué)院，北京100083；2.地下信息探測技術(shù)與儀器教育部重點實驗室，北京100083）

0 引 言

核磁共振（NMR）測井回波信號（CPMG Echo Trains）的采集不可避免地要受到儀器設(shè)備和井下環(huán)境噪聲的影響，記錄到的回波信號帶有很強的噪聲，給反演解釋工作帶來極大的不確定性。當(dāng)信噪比小于5時，直接反演便不能得出可靠的T2譜分布。Hansen等［1］針對噪聲影響的強度給出了離散的Picard條件，說明當(dāng)信噪比低于5時，回波信號與弛豫譜滿足的第1類Fredholm積分方程所反映的模型已經(jīng)完全不滿足離散Picard條件，也就是說因為噪聲的影響單從數(shù)學(xué)角度不能反演出真實可靠信息，必須對數(shù)據(jù)濾波或者重新采樣［1－2］。

實際生產(chǎn)中，數(shù)據(jù)采集采用不同相位信號疊加，壓制了很大一部分噪聲，對于解釋處理，疊加后的數(shù)據(jù)仍具有相當(dāng)高的噪聲。低信噪比核磁共振數(shù)據(jù)處理有翁愛華等［3］提出在反演中將基線偏移作為獨立的參數(shù)參與反演，并結(jié)合時間依賴濾波技術(shù)提高了長弛豫信息反演的準(zhǔn)確性。該方法將噪聲影響合理地抽象成反演參數(shù)并參與反演與校正，但沒能在小孔隙信息恢復(fù)上得到有效校正。Ma S等［4］提出了改進(jìn)的分窗小波降噪方法，并與中值濾波、普通小波濾波和FIR濾波的降噪效果進(jìn)行了對比，得出改進(jìn)的分窗小波方法能夠更好地提高帶噪信號的信噪比，但在反演信息準(zhǔn)確恢復(fù)上缺乏進(jìn)一步驗證。Ahmed O A W［5］提出了SLFT（Stable Linear Time－Frequence Transforms）方法，通過限制和平均變換系數(shù)實現(xiàn)了對分子藥物核磁共振數(shù)據(jù)成像上的降噪。對于低磁場高噪聲核磁共振測井?dāng)?shù)據(jù)，該方法目前僅具有借鑒意義。

本文針對低信噪比的核磁共振測井?dāng)?shù)據(jù)提出了基于奇異值分解法SVD（Singlar Value Decomposition）濾波的核磁共振數(shù)據(jù)變換反演方法BRD（Butler Reeds Dawson）。該方法利用Hankel矩陣與回波信號的相互構(gòu)建關(guān)系，把回波信號的Hankel矩陣作為SVD濾波處理的信號矩陣，實現(xiàn)了對回波信號的分解和濾波。結(jié)合BRD反演算法，利用SVD處理后的高信噪比全波列數(shù)據(jù)反演出流體的T2分布。在隨后的模擬數(shù)據(jù)試驗和實際數(shù)據(jù)處理中，可以看到這樣的分解和濾波在確保可動流體成分和提高束縛流體成分在反演信息準(zhǔn)確恢復(fù)上是有效的。

1 處理方法與原理

整個數(shù)據(jù)的處理流程包括回波信號分解為高信噪比部分和低信噪比部分；低信噪比數(shù)據(jù)濾波；低、高信噪比數(shù)據(jù)反演；反演譜疊加。

1.1 回波信號分解

設(shè)核磁共振采集的回波信號為X，對應(yīng)的T2譜為F。X可以分解為多項疊加的形式，如式（1）所示分解為X1和X2，相應(yīng)的T2譜為F1和F2，與F的關(guān)系如式（2）所示，有

式（1）稱為回波信號分解，這樣的分解處理意義在于，對于低信噪比回波信號，回波信號的有效分解能夠保證更多的有效信息與噪聲信號分離，提高主要信息在反演處理中恢復(fù)的準(zhǔn)確性。

圖1（a）給出了含有100個樣本的回波信號X，X分解為X1和X2的2個信號，它們之間滿足式（1）。圖1（b）給出了回波信號的T2譜，從左往右依次是X、X1和X2對應(yīng)T2分布，它們之間的關(guān)系滿足式（2）。

圖1 回波信號分解

1.2 SVD濾波處理

1.2.1 濾波原理

對于一個m×n維信號矩陣A，對其進(jìn)行SVD分解［6－7］，有

式中，U為m×m維正交矩陣；S為m×n維似對角矩陣；V為n×n維正交矩陣。當(dāng)m＝n，S＝Σ；m＜n，S＝［Σm×mOm×（n－m）］；m＞n，S＝［Σn×n O（m－n）×n］T。Om×（n－m）是m×（n－m）維 零 矩 陣，O（m－n）×n是（m－n）×n維零矩陣。Σ＝diag（σ1，σ2，…，σp，…，σr）為非負(fù)對角矩陣，它的對角元素是由從大到小排列的奇異值構(gòu)成。

通常，沒有噪聲的信號矩陣的秩很小，而含噪聲的信號矩陣是一個滿秩矩陣。受噪聲的影響，原來排列靠后很小的奇異值會變得很大，SVD濾波原理就是將排列在某一奇異值之后所有奇異值置0，通過重構(gòu)信號矩陣達(dá)到降噪的目的［8－9］。

為了描述置0奇異值分界點的選取，定義重構(gòu)信號矩陣的能量比為

式中，p作為置0奇異值的分界點，奇異值σp以后的奇異值將被置0，輸出信號占總能量比例為η。在實際處理中，為獲取最優(yōu)濾波信號，往往根據(jù)信號矩陣奇異值分布曲線的拐點選取p的位置，此時得到的能量比稱為最優(yōu)能量比。

1.2.2 信號矩陣

對于核磁共振回波數(shù)據(jù)X，它由N個等間隔采樣點構(gòu)成，表示為

考慮噪聲影響，令S表示期望信號，N0表示噪聲，S和N0與X有相同的采樣形式，式（5）可寫為

根據(jù)回波信號與Hankel矩陣的構(gòu)建關(guān)系［10］，令回波信號的Hankel矩陣作為反映該回波信息的矩陣，稱為信號矩陣?；夭ㄐ盘朮的信號矩陣H如式（7）所示，信號S和噪聲N0也可以寫出自己的信號矩陣，它們之間的關(guān)系H＝HS＋HN；HS和HN是S和N0對應(yīng)的信號矩陣，有

信號矩陣大小的選取應(yīng)滿足min（R，L）＞rank（H），對于數(shù)據(jù)大樣本信號，構(gòu)建信號矩陣應(yīng)盡可能減少行或列的維數(shù)，以減少計算量。

信號矩陣通過式（8）計算可恢復(fù)相應(yīng)回波信號

式中，g＝max（1，i－R＋1）；k＝min（i，N－R＋1）。

1.2.3 算法實現(xiàn)

回波信號的SVD濾波處理包含3個步驟。①構(gòu)建信號矩陣，利用式（7）構(gòu)建待處理回波信號的信號矩陣；②重構(gòu)矩陣，根據(jù)濾波目的選擇合適的能量比，得到降噪后的信號矩陣；③重構(gòu)信號，由回波信號濾波后的矩陣經(jīng)過式（8）得到重構(gòu)信號，此時的重構(gòu)信號即為經(jīng)過SVD濾波處理后的回波信號。

1.3 反演弛豫信息

核磁共振測井依據(jù)地層流體在外加磁場作用下的弛豫特性反映地層流體孔隙度和區(qū)分地層中可動流體與束縛流體。根據(jù)核磁共振測井原理，從核磁共振測井采集的CPMG回波信號反演T2分布是典型的病態(tài)問題［2，11］，該問題被描述為式（9）及考慮噪聲影響的離散形式（10），有

式中，xi為記錄的第i個回波幅度的回波信號；ti（其中i＝1，2，…，N）為第i個相等回波間隔的衰減時間；Tj對應(yīng)第j種流體的弛豫特性（其中j＝1，2，…，M），即弛豫時間T2；fj為第j種流體的T2幅度。TW為等待時間，R為T1／T2，方程（1）為第1類Fredholm方程。

測量信號xi不可避免地要受到噪聲的影響。小的噪聲干擾可能造成fj很大程度的不穩(wěn)定，使得T2的分布強度fj不能有效、唯一地確定。目前提出的很多反演算法和正則化方法旨在抑制解的不穩(wěn)定性和加強解的可靠性［12－15］，其中較為有效的算法是Venkataramanan等［13］提出BRD反演算法，該算法在目標(biāo)方程中增加一項范數(shù)比例懲罰項抑制噪聲影響。如目標(biāo)方程（11）所示，方程中增加了帶有平滑因子α的懲罰項

具體計算中，解的可行域范圍內(nèi)預(yù)設(shè)弛豫基確保解的合理性以及計算效率。例如，均勻分布在0.1～104ms范圍內(nèi)的30個離散值構(gòu)成弛豫基。為方便起見，方程（12）給出了矩陣形式，待求參數(shù)為弛豫時間T2的分布強度向量F

1.4 處理流程

核磁共振測井采集的回波信號受噪聲的影響主要表現(xiàn)，①在大、中孔隙可動流體成分造成基線偏移，在反演的弛豫譜上基線偏離表現(xiàn)為異常的長弛豫組分，這種影響在核磁數(shù)據(jù)信噪比較低時尤其明顯；②在微、小孔隙束縛流體成分產(chǎn)生畸變，在反演的弛豫譜上表現(xiàn)為跳變的尖銳假峰。提高低信噪比回波信號反演弛豫信息的準(zhǔn)確性，就是確保可動流體成分和提高束縛流體成分在反演信息準(zhǔn)確恢復(fù)。

首先，應(yīng)用SVD濾波高效的信噪分離能力進(jìn)行回波信號分解，得到高信噪比的回波信號X1和低信噪比的回波信號X2，高信噪比的回波信號X1含有大量的T2弛豫信息可以直接反演得到相應(yīng)的譜分布。其次，低信噪比的回波信號X2嚴(yán)重受到噪聲影響，需要再次經(jīng)過SVD濾波提取有用信息，然后反演得到相應(yīng)譜分布。最后，根據(jù)回波信號分解關(guān)系，疊加2部分反演的T2譜作為最終要求取的T2分布。

反演數(shù)據(jù)的處理過程中2次用到SVD濾波。

（1）反演的回波信號經(jīng)過SVD濾波實現(xiàn)回波信號的分解。SVD濾波算法采用低于最優(yōu)能量比的設(shè)置，保證分解得到高信噪比數(shù)據(jù)X1。

（2）由回波信號分解理論，由X2＝X－X1計算得到相應(yīng)的低信噪比回波信號X2。

（3）低信噪比數(shù)據(jù)X2的SVD濾波，提取有用信息?；夭ㄐ盘朮2含有較高噪聲，SVD濾波采用高于最優(yōu)能量比的設(shè)置，保證有用信息不被丟失。

（4）回波信號反演。2部分經(jīng)SVD濾波后的數(shù)據(jù)，反演得到相應(yīng)T2分布、F2和F1。

（5）反演譜的疊加。由F＝F1＋F2計算回波信號X的T2分布、F。

2 方法驗證與分析

根據(jù)鄧克俊等［11］提出的 FCD（Fluid Component Decomposition）回波串構(gòu)建弛豫信號的方法，取Gaussian函數(shù)作為基函數(shù)，高斯寬度（Full Width at Half Maximum）為0.9，考慮不同強度的噪聲影響，正演得到不同信噪比的回波信號。對于含有束縛水和中等黏度油的雙峰特征T2譜模型，束縛水的T2峰值在20ms處，峰值幅度為0.4，含水飽和度為40%；油的T2峰值在300ms處，峰值幅度0.6，含油飽和度60%；模擬采集間隔TE＝0.9ms，回波數(shù)N＝500；模擬得到的回波信號見圖2。

圖2 模擬信號處理結(jié)果對比

圖2中藍(lán)色實線為沒有加入噪聲的回波信號，即期望信號，黑色點線為帶噪聲的回波信號，紅色虛線為經(jīng)過SVD濾波處理反演前2部分疊加的回波信號。圖2（a）到圖2（e）分別為信噪比50、30、20、10和5的帶噪回波信號，相應(yīng)的經(jīng)過SVD濾波處理后的重構(gòu)信號與期望信號。從中可以看出濾波后的信號有效地恢復(fù)到了期望信號，但隨著信噪比下降過低時，SVD濾波處理仍舊保證了大弛豫成分信息的準(zhǔn)確性，小弛豫信息會造成少量的丟失。

同時，反演了相應(yīng)的弛豫譜。反演計算中，弛豫基分布區(qū)間0.3～3 000ms，布點數(shù)30（見圖3），藍(lán)色實線表示弛豫模型，黑色點實線表示帶噪回波信號直接反演得到的弛豫譜（Fn），紅色星實線表示經(jīng)本文濾波處理得到的弛豫譜（F）。

如圖3所示，應(yīng)用本文方法處理得到的T2譜更接近模型，反演譜峰的位置與模型近乎一致。隨著信噪比的降低，小弛豫成分峰值會有所丟失，但還是保證了弛豫位置的準(zhǔn)確性；而直接反演得到T2譜只在油峰位置與模型相近，隨著信噪比的降低，小弛豫位置會出現(xiàn)假峰。從對比的結(jié)果看，本文方法在低信噪比數(shù)據(jù)的處理上更接近真實情況。

圖3 模擬信號處理結(jié)果對比

為定量說明反演效果，令反演得到的T2譜用列向量f表示，T2譜模型用列向量m表示，定義相對誤差Er作為反演精度的度量，有

根據(jù)式（12），相應(yīng)反演譜的相對誤差見表1。由表1可知，隨著回波信號信噪比的降低，直接反演方法得到譜分布的相對誤差增加的很快，說明了當(dāng)信噪比低或過低時，直接反演結(jié)果的不可靠性；與直接反演方法相比，在較低信噪比情況下，采用本文濾波反演方法得到譜的相對誤差更小，誤差變化幅度也更小。

當(dāng)回波信號信噪比較高（SNR≥20）時，直接反演譜和經(jīng)濾波反演譜的流體弛豫位置都與模型相一致，但濾波反演譜的整體相對誤差要略大于直接反演結(jié)果，說明較高信噪比回波信號濾波處理時，由于信噪的分離難度加大，在設(shè)置濾除低強度噪聲的同時，部分有用信息也被濾除了。因此，較高信噪比回波信號處理時應(yīng)采用設(shè)置為全通能量的濾波反演模式或跳過濾波處理直接反演。

表1 反演譜相對誤差

3 油田實例分析

某區(qū)×井核磁共振測井?dāng)?shù)據(jù)，測量儀器為MRIL－Prime，采用雙TW測量模式識別可動流體，采集參數(shù)為，A組采用長磁化參數(shù)TW＝12.7s，B組采用短磁化參數(shù)TW＝2.0s，2組采用相同的采集間隔TE＝0.9ms和回波數(shù)NE＝500。2組回波信號的統(tǒng)計信噪比為，A組大約在2～19，B組大約在1～19，屬于較低的信噪比數(shù)據(jù)。

圖4為本文方法處理結(jié)果與商業(yè)軟件DPP處理結(jié)果的對比圖。圖4中，第1道GR、SP；第2道為地層電阻率Rt、沖洗帶電阻率Rxo；第3道為深度道；第4、5道為長TW，L和短TW，S回波信號的 DPP反演結(jié)果；第6道為DPP軟件處理的長TW，L與短TW，S弛豫差譜；第7、8道為本文濾波反演方法處理的長TW，L和短TW，S結(jié)果；第9道為本文濾波反演方法處理的長TW，L與短TW，S弛豫差譜?？梢钥闯?，在××55～××73m深度段與DPP處理結(jié)果相比，在小弛豫位置弛豫譜沒有出現(xiàn)噪聲引起的跳變尖銳峰；在大弛豫位置壓制了噪聲引起的基線偏移造成的大弛豫成分；在A組和B組的弛豫譜上能夠看到該層段存在大量的束縛流體和可動流體；從第6道和第9道的差譜信息上看，兩者處理的差譜信號都很弱，但文中方法處理結(jié)果仍能看出可動流體中含有一定量的烴。

圖4 實測數(shù)據(jù)處理結(jié)果對比

圖5 巖性柱狀圖

圖5為該層段錄井巖性及含油級別測試結(jié)果。錄井巖性顯示，該層段巖性為深灰、灰黑色泥巖、灰色粉砂質(zhì)泥巖與淺灰、灰褐色細(xì)砂巖呈略等厚互層，含油測試結(jié)果為油跡、油斑級別。該層段的測井響應(yīng)資料顯示，自然伽馬曲線值為64.7～89.7API，均值為78.0API。三孔隙度曲線變化平緩，聲波時差、密度均值分別為244.8μs／m、2.41g／cm3，中子16.8%～27.9%，說明地層孔隙發(fā)育，物性好。該層電阻率值略有變化，深感應(yīng)電阻率值為5.0～15.0 Ω·m，均值為6.8Ω·m，陣列感應(yīng)電阻率曲線明顯分異，120in＊非法定計量單位，1ft＝12in＝0.304 8m，下同陣列感應(yīng)電阻率值為4.9～14.9 Ω·m，孔隙度為13.8%，滲透率為0.62×10－3μm2，表明地層含油性較好，該層解釋為油水同層。

4 結(jié)論與認(rèn)識

（1）受噪聲影響，核磁共振測井采集的回波信號表現(xiàn)為在大、中孔隙可動流體成分造成基線偏移，在微、小孔隙束縛流體成分產(chǎn)生畸變；在反演的弛豫譜上表現(xiàn)為異常的長弛豫和小弛豫位置跳變的尖銳假峰。尤其當(dāng)信噪比低于5時，單從數(shù)學(xué)角度不能反演出真實可靠信息，必須對數(shù)據(jù)濾波或重新采樣。

（2）提出基于SVD濾波的BRD反演處理方法，在低信噪比（SNR≤20）保證了可動流體和束縛流體信息在反演處理中的有效恢復(fù)，提高低信噪比核磁共振數(shù)據(jù)反演弛豫信息的準(zhǔn)確性，在數(shù)值模擬和油田實例中都取得了令人滿意的效果。從油田實例的應(yīng)用效果看，與DPP處理結(jié)果相比，在小弛豫位置弛豫譜沒有出現(xiàn)噪聲引起的跳變尖銳峰，在大弛豫位置壓制了噪聲引起的基線偏移造成的大弛豫成分，在弛豫譜形態(tài)上能更精細(xì)地識別束縛流體和可動流體。

（3）SVD濾波處理利用Hankel矩陣與回波信號的相互構(gòu)建關(guān)系，根據(jù)噪聲強度選取適當(dāng)?shù)木仃嚹芰勘葔褐坪徒財嘁蛟肼曇鸬钠娈愔档淖兓哂泻軓姷男旁敕蛛x能力。

（4）實現(xiàn)了對回波信號的分解和濾波。在處理過程中，回波信號分解采用低于最優(yōu)能量比的SVD濾波設(shè)置，保證信號分解的有效性和回波信號的反演精度，回波信號濾波采用高于最優(yōu)能量比的SVD濾波設(shè)置，盡可能保留有效信息使之能夠在反演中恢復(fù)。

［1］ Hansen P C.The Discrete Picard Condition for Discrete Ill－Posed Problems［J］.BIT Numerical Mathematics，1990，30（4）：658－672.

［2］ 肖立志，謝然紅，廖廣志.中國復(fù)雜油氣藏核磁共振測井理論與方法 ［M］.北京：科學(xué)出版社，2012.

［3］ 翁愛華，李舟波，莫修文，等.低信噪比核磁共振測井資料的處理技術(shù) ［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，2003，33（2）：232－235.

［4］ MA S，KONG L，CHEN J.An Improved NMR Signal De－noising Algorithm Based on Wavelet Transform［J］.Computational Information Systems，2011，7（13）：4651－4659.

［5］ Ahmed O A W.Method for Removing Noise from Nuclear Magnetic Resonance Signals and Images：U.S.Patent 7，253，627［P］.2007－08－07.

［6］ Kalman D.A Singularly Valuable Decomposition：The SVD of a Matrix［C］∥College Math Journal，1996.

［7］ Aharon M，Elad M，Bruckstein A.SVD：An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation［C］∥Signal Processing，IEEE Transactions，2006，54（11）：4311－4322.

［8］ Di F，Changzhi L，Qinguang C.SVD Filter Based on Noise Singular Values Clustering ［C］∥Intelligent Computation Technology and Automation （ICICTA），2010International Conference，IEEE，2010，3：665－668.

［9］ Sanliturk K Y，Cakar O.Noise Elimination from Measured Frequency Response Functions ［J］.Mechanical Systems and Signal Processing，2005，19（3）：615－631.

［10］ 劉志鵬，趙偉，陳小宏，等.局部頻率域SVD壓制隨機(jī)噪聲方法 ［J］.石油地球物理勘探，2012，47（2）：202－206.

［11］ 鄧克俊，謝然紅.核磁共振測井理論及應(yīng)用 ［M］.東營：中國石油大學(xué)出版社，2010.

［12］ Butler J P，Reeds J A，Dawson S V.Estimating Solutions of the First Kind Integral Equations with Nonnegative Constraints and Optimal Smoothing［J］.SIAM Journal on Numerical Analysis，1981，18（3）：381－397.

［13］ Venkataramanan L，Song Y Q，Hurlimann M D.Solving Fredholm Integrals of the First Kind with Tensor Product Structure in 2and 2.5Dimensions［C］∥Signal Processing，IEEE Transactions，2002，50（5）：1017－1026.

［14］ LIAO G Z，XIAO L Z，XIE R H，et al.Influence Factors of Multi－exponential Inversion of NMR Relaxation Measurement in Porous Media ［J］.Chinese Journal of Geophysics，2007，50（3）：796－802.

［15］ Wang W，Li P，Ye C.Multi－exponential Inversions of Nuclear Magnetic Resonance Relaxation Signal［J］.Science in China Series A：Mathematics，2001，44（11）：1477－1484.