提升小波濾波在隨鉆核磁共振測井中的應用

徐志星，程晶晶，王光偉，張嘉偉

(1.華中科技大學 控制系，湖北 武漢430074;2.中海油田服務股份有限公司，河北廊坊065201)

隨鉆核磁共振測井以磁共振技術為基礎，在泥漿濾液侵入之前或者侵入很淺時對地層屬性進行探測，與傳統(tǒng)的電纜核磁共振測井相比，能更真實地反映儲層特征，評價準確性更高［1］。同時，在對低孔低滲、復雜巖性、特殊巖性等世界復雜油氣藏的勘探和開發(fā)方面，具有其他測井方式無法比擬的優(yōu)勢。隨鉆核磁共振測井被國際上三大測井公司列為最具發(fā)展?jié)摿Φ臏y井方法之一［2］。

受水平井或者大斜度井的限制，隨鉆核磁共振測井系統(tǒng)無法采用通信電纜傳輸數(shù)據(jù)，而僅能選取最大速率只有50 b/s的泥漿脈沖傳輸方式。受測井時序和傳輸速率的雙重限制，隨鉆核磁測井無法同電纜測井一樣傳輸大量的原始數(shù)據(jù)，而是經(jīng)過濾波、T2譜反演后，直接上傳數(shù)據(jù)量較少的井下測量結果。而通過有效的濾波算法提高性噪比是獲取準確T2譜的關鍵前提。因此，提高隨鉆核磁共振測井效果，需要設計一種計算量少、信噪比高的井下濾波算法。

基于此，筆者的實驗根據(jù)核磁共振測井回波串信號的頻譜特征，提出具有多分辨率分析的提升小波算法［3－4］，避免傳統(tǒng)小波算法的卷積運算，降低運算復雜度和內存的需求量，并可實現(xiàn)整數(shù)小波變換，便于電路實現(xiàn)。研究結果表明，提升小波濾波后，回波串的信噪比提高了65%，反演結果符合實際地質特征，孔隙度信息受噪聲干擾波動減小、測量一致性明顯提高。

1 核磁共振回波串信號的時頻特征

核磁共振測井的原始信號由一串幅度隨時間衰減的回波構成，如圖1所示，其信號模型為［5］:

式中:T2i為弛豫時間分量;φi為對應的孔隙度分量;t(n)=iTE，i=1，2，…，n，為采集第 i個回波時的時間;noise為噪聲。

圖1 回波串與回波信號示意圖

回波串信號幅度的指數(shù)衰減規(guī)律由T2描述，而單個回波信號的指數(shù)衰減由T*2來表征。T2為本征弛豫時間，T*2為實際橫向弛豫時間，它們之間的關系如式(2)所示，其中γ為旋磁比，ΔB0為靜磁場梯度，它反映了靜磁場的不均勻性。由式(2)可知T*2＜T2，因此它們的頻率范圍不同。

假設自旋回波信號的幅值為A(t)，它滿足式(3)，其中Mp為極化磁場的磁化矢量強度。從式(3)可知，回波信號幅值是指數(shù)衰減的曲線，發(fā)射180°激勵脈沖后，自旋回波在中心點達到最大值，每組回波信號的峰值點幅值的集合就是原始回波串信息。

假定原始回波信號經(jīng)濾波和峰值檢測后，獲得的回波串信息的集合為表示低通濾波器的截止頻率為wc，Ni為噪聲，Peak{}為取回波峰值，可得回波信號如式(4)所示:

由式(4)和式(5)可知，回波串信號與回波信號的噪聲同源。

回波串信號反映了本征弛豫過程，回波信號既包含本征弛豫過程，又包含散相弛豫過程，而散相弛豫過程比本征弛豫過程快，因此，散相弛豫過程的頻譜信號占據(jù)信號的高頻部分。另外，在對回波信號進行濾波時，為了避免信號的過度損失，保留了較多的散相弛豫占據(jù)的高頻噪聲。因此從頻譜上看，噪聲頻譜和散相弛豫頻譜是重疊的，而與本征弛豫頻譜有明顯區(qū)分，回波串頻譜并不包含散相弛豫頻譜和高頻噪聲頻譜，因此有必要對回波串進行濾波。回波串與回波信號頻譜如圖2所示。

從圖2可以看出，要想有效濾除回波串中的噪聲，必須選擇合適的截止頻率wc，使它滿足式(6)。

圖2 回波串與回波信號頻譜

由于回波頻率是變化的，導致wc難以確定。而小波變換濾波與其他濾波方式相比的優(yōu)勢在于它可以根據(jù)信號的頻譜特征，在時域和頻域內進行多分辨率分析，在保證時頻窗面積不變的條件下，改變窗口形狀以確定不同時間域下合理的截止頻率，從而將噪聲和信號區(qū)分開。因此，小波濾波的特點決定了它在核磁共振測井回波串濾波方面具有獨特的優(yōu)勢。

母乳是早產(chǎn)兒的最佳食糧。目前主張早產(chǎn)兒盡早母乳喂養(yǎng)，但是暖箱內寶寶不適宜母乳喂養(yǎng)，一般是給予特制的早產(chǎn)兒配方奶粉。早產(chǎn)兒配方奶粉不僅保留了母乳的許多優(yōu)點，使蛋白質、糖、脂肪等營養(yǎng)素易于消化和吸收，和其他配方奶粉相比，還適當?shù)靥岣吡藷崃?，強化了多種維生素和礦物質，以補充母乳對早產(chǎn)兒營養(yǎng)需要的不足。但早產(chǎn)兒配方奶粉缺乏母乳中的許多生長因子、酶和IgA等。

2 提升小波算法

2.1 提升小波算法的原理

提升小波方案將小波分解分成分裂、預測和更新 3 個階段［6－7］。

(1)分裂。將輸入信號si按照奇偶性分解成兩個子集 si－1和 di－1，其中 si－1為小波的低頻近似分量，di－1為小波的高頻細節(jié)分量。

(2)預測。根據(jù)原始信號數(shù)據(jù)之間的相關性，利用與數(shù)據(jù)無關的濾波器P對偶數(shù)序列作用后作為奇數(shù)序列的預測值，奇數(shù)序列的預測值與實際值之間的差為誤差信號，它反映了兩者之間的逼近程度。一般用該差值作為di－1，預測表達式為 di－1=di－1－ P(si－1)，如果預測合理，差值數(shù)集將比原始的di－1包含更少的信息，這樣可以用更小的子集si－1和小波子集di－1來代替原始信號si，di－1在小波多分辨率分析中被稱為小波系數(shù)。

(3)更新。為使預測得到的信號的一些全局特性(如均值、消失矩等)繼續(xù)保持，可以使用上述的預測值對si－1進行更新，因此需要通過算子U 更新 si－1，si－1被稱為尺度系數(shù)，更新表達式如式(7)所示:

由小波分解得到小波系數(shù)和尺度系數(shù)后，對小波系數(shù)進行閾值處理。閾值的設定有硬閾值和軟閾值兩種方式［8］，硬閾值將大于閾值的小波系數(shù)保留，小于閾值的小波系數(shù)置零，雖然會帶來間斷點，但是能較好地保留信號的細節(jié)部分，表達式如式(8)所示:

軟閾值將大于閾值的小波系數(shù)與閾值做差，小于閾值的小波系數(shù)置零，避免了硬閾值方式造成的間斷點，但是細節(jié)系數(shù)會丟失，表達式如式(9)所示。

基于以上分析，采用硬閾值方式，根據(jù)Birge－Massart策略［9］，假設分解層數(shù)為 m，對于第 k層系數(shù)，按照式nk=M/(m+2－k)α保留絕對值最大的nk個系數(shù)，其中，M為第m層近似系數(shù)的長度，α取值為3?；夭ù饕堑皖l信號，經(jīng)小波變換后成為高頻系數(shù)，而Birge－Massart閾值設置策略能有效保留高頻系數(shù)，濾除絕對值較小的低頻系數(shù)，保留有效信號。

提升小波重構是小波分解的逆過程，分為反更新、反預測和合并3個步驟，合并得到的信號就是原始信號si。

2.2 基于FPGA的提升小波算法的硬件實現(xiàn)

首先對9/7小波多相位矩陣因式分解，得到分解因式，如式(10)所示:因此，可得小波變換提升分解算法為:

為方便在FPGA中對小波系數(shù)移位，實現(xiàn)乘法運算，提高運算速度，降低資源量，對小波系數(shù)近似處理。理想小波系數(shù)、近似小波系數(shù)及對應的二進制表示如表1所示。

可由提升小波分解可逆向推導出提升小波重構表達式:

綜上所述，小波提升算法FPGA實現(xiàn)框圖如圖3所示。

圖4為二維核磁共振測井儀事件控制采集電路板實物圖，采用FPGA+DSP架構，以該電路板為硬件平臺，實現(xiàn)了激勵信號的發(fā)射以及回波信號的采集和處理，回波串信號采用提升小波算法濾波，該算法工作頻率為100 M，采樣頻率為8倍回波信號頻率，最大值約為 40 M，共占用A3P1000約40%的資源量。

表1 提升小波參數(shù)表

圖3 提升小波算法的FPGA實現(xiàn)框圖

圖4 二維核磁共振測井儀事件控制采集電路板實物圖

3 實驗

該實驗采用自行研發(fā)的隨鉆核磁共振測井系統(tǒng)樣機分別在實驗室刻度筒和標準油井中進行測量，并采用筆者設計的井下提升小波變換算法對其做濾波處理，并同時對比分析了濾波前后的回波串、反演結果和孔隙度。

圖5為實驗室原始回波串與小波濾波后的回波串曲線。圖6為標準油井原始回波串與小波濾波后的回波串曲線。從圖5和圖6中可以看出，小波濾波前，原始的回波串曲線的噪聲大，經(jīng)計算，圖5中原始回波串SNR=16，而圖6中SNR=20，小波濾波以后，圖5中回波串SNR=30，而圖6中SNR=33，回波信號被完整提取出來，噪聲得到抑制，利用小波去噪能較好地去除回波信號中夾雜的隨機噪聲，提高信噪比。

圖5 實驗室原始回波串與小波濾波后的回波串曲線

圖6 標準油井原始回波串與小波濾波后的回波串

目前隨鉆核磁共振測井儀中通常采用奇異值分解法做譜反演。奇異值分解法在信噪比較低時，會大量刪減原始數(shù)據(jù)，漏掉很多重要信息［10］，但在信噪比較高時，能得到較好的反演結果。

圖7、圖8分別為實驗室刻度筒和標準油井數(shù)據(jù)濾波前后的反演結果對比圖。濾波前，T2譜都出現(xiàn)了畸變。由于信噪比較低，部分有用數(shù)據(jù)被刪除，同時某些奇異點反演后不穩(wěn)定，造成部分反演數(shù)據(jù)丟失，沒有顯示完整的T2譜信息。由于實驗室與標準油井數(shù)據(jù)相比衰減明顯，并且信噪比比標準井高，因此畸變比標準井小。濾波后，回波串信噪比均大于30?？梢钥吹?，在高信噪比的情況下，T2譜反演的結果與低信噪比相比，沒有出現(xiàn)畸變，反演曲線光滑，對稱性好。提升小波濾波提高信噪比后，反演效果良好。該反演結果與地層解釋軟件中的解釋結果一致，圖7和圖8均正確反映了實際地層的油、水特征。

核磁共振自旋回波串的初始幅度或者T2譜分布曲線圍成的面積與探測范圍內孔隙流體中的氫原子核數(shù)量成正比［11］，該值經(jīng)過刻度以后就是孔隙度值。該實驗在T2譜反演后，可計算該曲線圍成的面積從而得出孔隙度信息。

圖7 實驗室原始回波串T2譜和小波濾波后T2譜

圖8 標準油井原始回波串T2譜和小波濾波后T2譜

圖9 為實驗室孔隙度信息，圖10為標準油井孔隙度信息。圖9為刻度筒中重復測量兩次的孔隙度信息，圖10為井下457～467 m測量兩次的孔隙度信息。濾波前可以看出:①兩次測量的孔隙度受到噪聲干擾波動較大，并且出現(xiàn)了一些異常點;②在同一深度處孔隙度值相差較大，重復性差。這些問題會導致地質解釋出現(xiàn)偏差。濾波后孔隙度波動明顯減小，并且沒有出現(xiàn)奇異數(shù)據(jù)，曲線比原始數(shù)據(jù)計算的孔隙度光滑，兩次測量的結果一致性明顯提高，可以看出，實驗室中孔隙度曲線經(jīng)過濾波后基本重合，標準井孔隙度曲線濾波后從463～467 m基本是重合的，而從457～463 m兩次測量的孔隙度值相差也比濾波前小，說明小波濾波以后，噪聲對孔隙度的影響降低。兩次測量的一致性，使得實際地層的真實特征更有說服力。對于標準油井，從圖10中可以看出，在井下463 m和458 m附近，含有較大孔隙，應該含有油氣水等物質，實際的T2譜顯示在463 m和458 m附近，確實含有束縛水和油。

4 結論

隨鉆核磁共振測井中信號受噪聲干擾較嚴重，導致地層分析結果受到很大影響，因此需要對原始的回波串進行濾波。筆者根據(jù)核磁共振回波串信號的頻譜特征，提出利用小波變換的多分辨

圖9 實驗室原始回波串孔隙度和小波濾波后孔隙度

圖10 標準井原始回波串孔隙度和小波濾波后孔隙度

率的特點對信號進行濾波，為方便硬件實現(xiàn)，將提升小波算法引入到隨鉆核磁測井中，并將實驗室刻度筒數(shù)據(jù)、標準油井測井數(shù)據(jù)濾波前后的回波串、T2譜和孔隙度計算的結果進行了對比。結果表明，小波濾波既保留了原始回波串特征，又將信噪比提高了65%，有效抑制了噪聲對測井結果的影響。

［1］邵才瑞，唐海全，張福明.隨鉆測井解釋評價技術現(xiàn)狀和研究方向［J］.大慶石油地質與開發(fā)，2010，29(5):159－162.

［2］肖立志，謝然紅，于娛嬌，等.核磁共振測井儀器的最新進展與未來發(fā)展方向［J］.測井技術，2003，10(2):130－136.

［3］盧春暉，王英健，陳麗.提升小波算法的FPGA硬件實現(xiàn)［J］.微機處理，2010，38(2):45－47.

［4］李美珊，張宏.9/7提升小波算法結合SPIHT算法的遙感圖像壓縮［J］.科技傳播，2012，50(23):240 －242.

［5］GEORGE R，XIAO L Z，MANFRED G.NMR logging principles and applications［M］.［s.n］:Halliburton Energy Services Pulication，1999:53 －54.

［6］李玲玲，丁明躍，周成平，等.一種基于提升小波變換的快速圖像融合方法［J］.小型微型計算機系統(tǒng)，2005，30(6):140 －141.

［7］陳升來，黃廉卿，郭靜寰.基于整型小波變換的圖像處理及 DSP實現(xiàn)［J］.光學精密工程，2006，14(3):498－500.

［8］JIANGM Y.Wavelet threshold optimization with artifi－catial fish swarm algorithm［J］.Neural Networks and Brain，2005(1):569 －572.

［9］唐炬，李偉，歐陽有鵬.采用小波變換奇異值分解方法的局部放電模式識別［J］.高電壓技術，2010，36(7):1687－1688.

［10］劉宏強，王祝文.基于奇異值分解法的核磁測井解譜算法［J］.吉林大學學報，2007，37(3):135 －136.

［11］HENK V A，TOM S ，F(xiàn)RANK J.MRIof intact plants［J］.Photos Ynthesis Research，2009(102):213 －222.