各向異性擴散濾波遠探測聲波測井圖像降噪方法*

傅艷莉 李 超 陳 浩 曹雪砷

(1 中國科學院大學 北京 100049)

(2 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)

(3 北京市海洋深部鉆探測量工程技術研究中心 北京 100190)

0 引言

近年來，隨著石油、天然氣及礦藏勘探和開發(fā)需求的增長及測井技術的發(fā)展，聲波遠探測測井成為井旁構造體識別的熱點技術[1]。常規(guī)測井項目探測井外信息較淺，難以反映距井較遠距離的地質(zhì)異常體的發(fā)育情況。聲波遠探測技術可以對井周圍數(shù)十米甚至上百米范圍進行探測，利用地層反射信息進行井周成像，探測到地層中裂縫、斷層等，大大豐富了研究領域，已經(jīng)發(fā)展成一種深部復雜油氣勘探領域中不可缺少的先進技術[2-3]。然而由于反射波幅度低、受幅度強的井筒直達波干擾等原因，反射波圖像具有低信噪比的特點，噪聲主要是殘留的井孔模式波和隨機噪聲(包括地層噪聲和濾波處理過程中引入的噪聲)，這些噪聲大都是非相干的。傳統(tǒng)的圖像去噪方法，如均值濾波、中值濾波等，在去除噪聲的同時會造成裂縫、斷層等地質(zhì)體結構的模糊。因此，為了解決去噪同時兼顧構造邊緣及細節(jié)的問題，前人引入了擴散濾波處理遠探測測井圖像[4]，取得了較好效果。但目前尚缺少對該方法適用性及其特點的系統(tǒng)性研究和分析，而采用各向異性擴散濾波對地層數(shù)據(jù)進行處理，在去除噪聲的同時可以更好地保留裂縫、斷層等地質(zhì)體結構[5]。

基于擴散方程的噪聲衰減方法是20 世紀90年代提出的一種數(shù)字圖像處理方法，是圖像降噪的研究熱點。與傳統(tǒng)的線性濾波方法相比，擴散濾波在去噪的同時保留圖像的邊緣細節(jié)，廣泛應用于圖像降噪領域的研究。其中最經(jīng)典的是Perona和Malik提出的各向異性擴散模型[6]。但是P-M模型的邊緣保持效果不明顯，后來又出現(xiàn)了許多改進模型，其中Weickert[7]模型用結構張量的形式表示擴散系數(shù)，使得各向異性擴散濾波廣泛應用于地震勘探的圖像處理中。Bakker[8]提出將基于圖像的各向異性擴散濾波應用到地震資料的預處理，并且提出幾何張量結構算法，用地層的幾何結構張量中包含的豐富地層信息，來描述地層的內(nèi)部結構，計算檢測斷裂系統(tǒng)。Fehmers等[9]把各向異性擴散濾波引入地震信息資料結構信息的邊緣保持。Lavialle等[10]提出地震斷層保持的擴散濾波方法，使得處理之后的地震數(shù)據(jù)中斷層信息更加突出。孫夕平等[11]、王緒松等[12]將二階導數(shù)引入擴散張量，探索了一致增強擴散濾波技術在二維地震剖面的保邊濾波應用。楊培杰等[13]提出一種方向性邊界增強技術，對于斷層信息起到了很好的保持作用。張爾華等[14]對三維地震資料非線性各向異性擴散濾波方法進行了初步探討，通過控制在斷層等不連續(xù)區(qū)域的濾波程度，使得各向異性擴散濾波器具有較好的保邊處理性能。各向異性擴散濾波在地震圖像的處理中已經(jīng)有了廣泛的應用，本文研究該降噪方法在聲波遠探測處理的適用性。

本文首先分析了遠探測聲波測井的殘余噪聲問題，通過梳理近年來各向異性擴散濾波的發(fā)展脈絡，從基本原理入手并針對其在遠探測圖像處理中的應用研究展開分析，對比具有不同擴散張量的濾波方法及濾波參數(shù)的處理效果，建立了適合于遠探測聲波測井信號的處理方法流程，最后通過實際數(shù)據(jù)的處理分析驗證了方法的效果，并對該技術今后改進和發(fā)展方向進行了展望與討論。

1 遠探測聲波測井殘余噪聲分析

由于反射波幅度低、受幅度強的井筒直達波干擾等原因，反射波圖像具有低信噪比的特點，原始波形經(jīng)過反射波分離處理后，廣義上的噪聲主要是殘留的井孔模式波和隨機噪聲(包括地層噪聲和濾波處理過程中引入的噪聲)，除了殘余模式波外，噪聲大都是非相干的，對后續(xù)的處理和地質(zhì)體識別造成干擾。以深部地層采集的偶極遠探測實際數(shù)據(jù)為例，如圖1(a)所示，原始的反射波分離后的剖面信噪比較低，部分反射體特征淹沒在較強的噪聲中，除了10 ms以內(nèi)的在軸向上表現(xiàn)為非相干的殘余彎曲波外，圖中可見兩種斜率的反射波，如黑色圓圈和紅色圓圈所示，其特征表現(xiàn)為在時間上延伸較長，頂部和井孔模式波重合。為了說明噪聲特征并區(qū)別其來源是否為井外反射，取其中3 個深度點陣列波形進行相關速度分析。圖1(b)為對圖像6350 m深度的8道波形(見圖1(a)區(qū)域A)做的相干分析，由于未見明顯反射，分析表明其是非相干的隨機噪聲。圖1(c)為對深度6245 m 的陣列波形(見圖1(a)區(qū)域B)進行相關分析，正向陣列(從1 到8)未見相干性，但反向陣列(從8 到1)得到其視速度在1400 m/s 左右，由于反射波在接收陣列上的入射角較小，一般具有較高的視速度，因此確定其為沿井軸反向傳播的反射斯通利波(下行)。同理，圖1(d)是對另一斜率的6310 m 深度波形(見圖1(a)中區(qū)域C)相干分析的結果，得到其下行視速度為3100 m/s 左右，驗證其為沿井軸傳播的反射橫波。

圖1 反射波分離后剖面及3 個深度點的相干分析Fig.1 Reflection wave separation profile and coherence analysis of three depth points

以上分析表明，隨機噪聲、殘余模式波和井孔模式波反射存在于遠探測聲波測井的波形中，尤其是沿井軸傳播的反射干擾常被誤以為是井外反射。這些噪聲的存在，使獲取的圖像難以表示真正的井外反射體特征。而一般而言，經(jīng)過疊加后，反射信號得到一定程度的增強，這些噪聲特征會減弱，總體表現(xiàn)為非相干性，但其仍然會掩蓋部分真實反射波。因此采取適當?shù)臑V波方法進一步壓制噪聲并增強反射信號是實際中需要解決的關鍵難題。針對該問題，本文引入各向異性擴散濾波方法來進行處理分析。

2 各向異性擴散濾波

2.1 各向異性擴散濾波原理

各向異性擴散濾波是一種基于偏微分方程的圖像處理方法。在物理學中，擴散定義是指物質(zhì)分子從高熱量區(qū)域向低熱量區(qū)域轉移，最終達到平衡狀態(tài)。當把擴散過程應用到圖像處理中，研究對象就變成了圖像的像素值。最簡單的擴散為線性各向同性擴散，Koenderink[15]證明圖像與高斯核函數(shù)的卷積與熱傳導方程是一致的，將物理學的熱擴散偏微分方程應用到圖像處理中可以表示為

式(1)中，U是待處理的數(shù)據(jù)；a為擴散系數(shù)，在各向同性擴散中a是一個常數(shù)；t為擴散時間。由于上述濾波器的擴散系數(shù)為常量，即對圖像各個方向的擴散強度是一樣的，不能起到保護圖像邊緣和細節(jié)信息的作用。為了克服這個缺點，Perona等[6]提出了非線性的各向異性擴散濾波模型：

式(2)中，將擴散系數(shù)設定為關于圖像梯度的非負單減函數(shù)g(‖?U ‖)，滿足g(0)=1，g(∞)=0。但是P-M 模型由于選取的擴散系數(shù)是一個標量，沒有考慮到圖像的紋理特征，使得擴散結果存在階梯效應。

Weickert[7]將結構分析引入擴散濾波中，將擴散系數(shù)變?yōu)榻Y構張量，使擴散隨方向變化而變化，在去噪的同時保護了圖像的線性紋理特征。構建了如下擴散模型：

式(4)中，*為卷積運算符，G為高斯核函數(shù)，σ為噪聲尺度，用以控制濾波過程中需要保護對象的最小尺寸，ρ為整合尺度，一般整合尺度應大于噪聲尺度。結構張量Sρ是2×2 的對稱半正定矩陣，對其進行特征分解：

式(5)中，v1和v2為梯度結構張量的兩個特征向量,v1是梯度變化最大的方向，v2是垂直于梯度變化的方向；λ1和λ2為梯度結構張量相對應的非負特征值，且λ1≥λ2，特征值反映了信號在局部特征方向的變化強度。

擴散的方向由結構張量的特征向量決定，要保證擴散張量的特征向量與結構張量的特征向量相一致：

式(6)中，μ1和μ2是兩個擴散方向的擴散強度，主要有3 種設置方法：相干增強擴散模型[16]、邊緣增強擴散模型[17]和混合擴散模型[18]。

相干增強擴散模型在執(zhí)行過程中，會使得待處圖像的每個像素都沿著相干方向進行擴散，它的特征值取值與結構張量Sρ的兩個特征值λ1和λ2的關系有關，設定為

式(7)中，參數(shù)α為一很小的正實數(shù)，控制擴散濾波過程中對不連續(xù)結構特征的保持程度，一般取0.001，使圖像梯度變化較大方向的擴散強度很小，保護圖像的結構信息；k表征圖像不同方向的相干性，表示為k=(λ1-λ2)2；c1為閾值，通常取為1。

在圖像相干性較高的區(qū)域，μ2約等于1，也就是說沿著特征向量v1方向的擴散強度很小，近似于不進行擴散，主要沿著特征向量v2的方向進行擴散；而在各向同性區(qū)域，v1和v2兩個方向的擴散強度都很小，因此可以有效地保護圖像的結構信息。

邊緣增強擴散體現(xiàn)在沿邊緣方向擴散大而垂直邊緣方向擴散小，以保護圖像邊緣結構。其擴散張量的特征值設定為

式(8)中，c2和e為正實數(shù)，分別取3 和0.02；gradA為圖像U兩個方向梯度的平方和。μ2=1即在垂直于梯度方向擴散強度始終為1；而在梯度方向的擴散強度取決于圖像兩個方向梯度的平方和gradA,這個值越大，μ1越小，說明梯度方向的擴散越小。有效地識別圖像邊緣，對邊緣起到保護與增強作用。

混合各向異性擴散是將相干增強擴散與邊緣增強擴散結合起來。其擴散張量的特征值設定為

式(11)中，h為一正實數(shù)，取0.5；其他變量與式(7)、式(8)中取相同值。

式(3)可以用有限差分進行求解，主要有兩種離散方法：半隱式離散形式[16]和顯式離散形式[16]。利用差分代替微分，實現(xiàn)對圖像U對擴散時間的求導，得到非線性擴散方程的半隱式離散形式：

式(13)中，Un+1和Un分別是圖像在第n+1 和第n次迭代得到的擴散結果；Δt為迭代步長，取值越小表示擴散越慢，為了保證圖像計算精度和穩(wěn)定性，一般取一個很小的正實數(shù)。

使用正向差分近似對式(3)離散化，這樣可以直接從前一個時間級別計算出新的時間級別上的值，而不用求解線性或非線性方程，得到顯式離散形式：

式(14)中，An*Un是div(D?u)的離散形式，即將圖像與一個時變與空變的模板An進行卷積。An如表1[17]所示，其中，a、b和c是擴散張量中的3 個時、空變的值。

表1 卷積模板Table 1 Convolutional template

關于離散化方式，在同樣的精度要求下，顯式有限差分方法比半隱式格式更加有效；而半隱式擴散可以采用較大的步長，穩(wěn)定性要求比顯式低。

2.2 三種擴散模型的對比分析

相干增強擴散模型將μ1設置為一個極小的常數(shù)α，即在特征向量v1方向的擴散很小，而在v2方向的擴散取決于圖像的相干性，圖像相干性強的區(qū)域(圖像邊緣)μ2接近于1，該模型適用于圖像噪聲比較強的情況，用于增強線狀紋理；邊緣增強擴散模型正好相反，將μ2設置為常數(shù)1，即在垂直梯度方向的擴散始終保持最強，而在v1方向的擴散取決于圖像兩個方向梯度的平方和，該模型適用于圖像噪聲較弱的情況，用于平滑噪聲，同時增強邊緣；混合擴散模型是將以上兩種模型結合起來，根據(jù)不同的噪聲強度和圖像特征選擇不同的擴散模型。為了更形象地說明3 種擴散張量模型的差異，以指紋圖像為例對其進行進一步對比分析。如圖2(c)、圖2(f)中紅色箭頭所示位置，相比于其他兩種擴散模型，相干增強擴散濾波后間斷圖像的連續(xù)性增強，圖像的紋理特征得到加強。

圖2 指紋圖像的濾波效果Fig.2 The filtering effect of fingerprint image

2.3 算法實現(xiàn)及參數(shù)選擇

根據(jù)以上分析，設計的圖像擴散迭代算法為

(1)輸入原始地層圖像U0，作為迭代算法的初始值，選擇迭代次數(shù)K，設定當前迭代次數(shù)k=0；

(2)對圖像做高斯平滑處理，得到Uσ；

(3)求Uσ的梯度結構張量，并對其做高斯平滑處理，得到Sρ；

(4)對Sρ做特征分解，得到兩個特征值λ1和λ2及其對應的特征向量v1和v2；

(5)選擇一種擴散張量模型計算擴散張量D；

(6)選擇一種迭代方式得到迭代k+1次結果；

(7)如果k ＜K，則令k=k+1，返回步驟(2)；否則，結束迭代；

(8)輸出濾波結果。

在步驟(2)中對高斯圖像做尺度為σ的高斯卷積，來控制濾波過程中要保護對象的最小尺寸；在步驟(3)中對梯度結構張量做尺度為ρ的高斯卷積，可以消除噪聲對計算結果的干擾。兩個尺度參數(shù)的選擇要根據(jù)待處理的數(shù)據(jù)特征確定，但要保證ρ ＞σ，當?shù)貙訑?shù)據(jù)的裂縫結構比較豐富時，可以將兩個尺度參數(shù)取小一點；當數(shù)據(jù)受到噪聲干擾比較嚴重時，可以將尺度參數(shù)取大一點。在步驟(5)中利用梯度結構張量的特征值和特征向量構建擴散張量D，根據(jù)圖像特性選擇合適的擴散張量；在步驟(6)步中根據(jù)要處理的數(shù)據(jù)選擇合適的迭代方法。迭代次數(shù)的選取要根據(jù)實際數(shù)據(jù)和預期的處理效果選定，保證處理后的數(shù)據(jù)信噪比有很大提高，裂縫得到很好的保持。同時，為了保證圖像的收斂性，迭代步長Δt應該盡量小。

3 模擬數(shù)據(jù)處理分析

為了說明各向異性擴散濾波去除噪聲同時保持圖像細節(jié)的作用，利用模擬裂縫數(shù)據(jù)進行驗證。原始模型數(shù)據(jù)為井外兩列平行的裂縫的聲波遠探測偏移成像結果(圖3(a))，在原始數(shù)據(jù)模型加入強噪聲來模擬實際情況，噪聲是使信噪比為0.1 的高斯噪聲(圖3(b))。可以看到，由于噪聲的影響，裂縫特征變得模糊，給裂縫識別造成干擾。

圖3(c)是用均值濾波方法處理的，在濾除噪聲的同時，裂縫邊緣變得不清晰，部分裂縫被消除，不利于裂縫結構的識別與解釋。圖4 是利用相干增強擴散張量模型、邊緣增強擴散張量模型和混合擴散張量模型處理加噪數(shù)據(jù)得到的結果，采用的都是顯式擴散。對比3 種模型的處理效果可以看出，3 種模型都能濾除圖像中的大部分隨機噪聲，同時裂縫結構得到很好的顯示。在深度15 m 附近的右邊裂縫處，相干增強擴散張量模型比其余兩種擴散張量模型更加清晰，處理效果更好，信噪比更高。為了更加清晰地看出3 幅處理后圖像的區(qū)別，將其分別與原始圖像做差得到圖5(a)、圖5(b)、圖5(c)，可以看到5(a)右側裂縫顏色更淺，說明裂縫結構得到了更好的保留，誤差更小。而右側裂縫相對左側幅度較弱，說明相干增強擴散模型更適合于弱反射結構在強噪聲環(huán)境下的處理，這是由于相對于梯度來說，相干性受噪聲影響更小。圖6(a)、圖6(b)分別是采用相干增強擴散張量的半隱式擴散模型和顯式擴散模型對加噪數(shù)據(jù)的處理結果，顯式擴散模型的濾噪效果明顯優(yōu)于半隱式擴散模型。不僅有效濾除了圖像中的噪聲，保持了圖像的分辨率，提高了信噪比，而且圖像中的裂縫區(qū)域更加明顯清晰，增強了裂縫特征的細節(jié)，具有很好的保持邊緣的特性，增加了對裂縫的識別能力。一般而言，隱式差分比顯示的穩(wěn)定性要求要低，但是為了保持圖像的細節(jié)特征，一般采用較小的時間步長，處理結果表明此時顯式模型的效果要優(yōu)于隱式的。

圖3 模擬數(shù)據(jù)及均值濾波效果Fig.3 Simulation data and mean filtering result

圖4 不同擴散張量模型的濾波效果Fig.4 Filtering results of different diffusion tensor models

圖5 不同擴散張量模型的濾波圖像與原圖做差Fig.5 Filtering results of different diffusion tensor models are subtracted from the original image

圖6 不同迭代模型的濾波效果Fig.6 Filtering results of different iteration models

為了定量地比較處理前后圖像質(zhì)量的改善效果，計算加噪圖像與處理后圖像的信噪比：

式(15)中，u0表示處理前的圖像，u表示處理后的圖像。選取不同的參數(shù)進行實驗，并計算處理前后的信噪比和能量保持的關系。實際計算表明，處理后信噪比都增加了5～10 dB，能量保持60%以上,從圖7 可以看出，在一定范圍內(nèi)信噪比隨著迭代次數(shù)的增加而增加。相干增強顯式擴散模型雖然信噪比增加較慢，但是隨著迭代次數(shù)增加最終達到的信噪比最高。當?shù)螖?shù)取35 時，圖3 三種模型處理后的信噪比分別為9.03985 dB、8.6109 dB 和8.8968 dB。說明選取合適的參數(shù)，不僅可以抑制噪聲還可以保留邊界信息。

圖7 不同擴散模型信噪比隨迭代次數(shù)的變化Fig.7 The signal-to-noise ratio of different iteration models varies with the number of iterations

該實例說明，基于相干增強擴散張量的顯式各向異性擴散濾波可以明顯改善資料品質(zhì)，有效識別裂縫結構。在處理中，選擇σ= 2、ρ= 5、Δt= 0.1、K= 35，明顯看出選擇合適的參數(shù)與模型，各向異性擴散濾波處理之后的圖像，信噪比得到提高，圖像更加清晰，有利于裂縫結構的識別。因此，本文建立了相干增強擴散張量結合顯式擴散模型的處理步驟，并依據(jù)該組合進一步對實際數(shù)據(jù)進行處理。

4 實際數(shù)據(jù)結果分析

下面是對遠探測實測數(shù)據(jù)進行各向異性擴散濾波的一個實例。由于地下條件復雜，數(shù)據(jù)采集深度較大，在深部高溫環(huán)境下?lián)Q能器性能受到影響，信號質(zhì)量下降，噪聲增加。原始的反射波分離后的疊加剖面圖像信噪比非常低,如圖8(b)所示，雖然有部分反射體特征存在，但是均淹沒在較強的噪聲里面，與具有連續(xù)性的裂縫帶相比，噪聲在圖像中的分布和大小不規(guī)則，更加雜亂且具備隨機性和不連續(xù)性。圖8(c)和圖8(d)分別是采用相干增強擴散張量經(jīng)過半隱式和顯式各向異性擴散濾波之后的結果，3 幅圖的色標范圍一致。對比處理前后的剖面圖像可以看出，處理之后的剖面相對于原始剖面，信噪比得到明顯提高，裂縫帶連續(xù)性也得到增強，特別是深度6310～6340 m 范圍效果更明顯。沿井壁附近傳播并在時間上衰減較小的殘留的反射斯通利波和反射橫波(如圖8(b)中黑色圓圈和紅色圓圈所示)，由于經(jīng)過反射波分離并疊加后其相干性較弱，也均得到了有效壓制。圖8(d)相比于圖8(c)，更好地去除了反射斯通利波，如圖8(b)第二個黑色圓圈所示。

圖8 反射波分離結果剖面的各向異性擴散濾波效果對比Fig.8 Comparison of anisotropic diffusion filtering results of reflected wave separation result profile

圖9 是另一深度段的反射波疊加剖面的處理結果，進一步驗證了基于相干增強擴散張量的顯式擴散模型在處理遠探測實測數(shù)據(jù)上的優(yōu)越性，信噪比更高，裂縫更加清晰，尤其是井壁附近的殘留彎曲波噪聲得到了很好的壓制。本文中的數(shù)據(jù)降噪方法應用于實際數(shù)據(jù)，取得較好效果，改善了成像圖的視覺質(zhì)量，提高了成像中反射體的清晰度，增強了對井外構造的解釋度。

圖9 反射波分離結果剖面的各向異性擴散濾波效果對比Fig.9 Comparison of anisotropic diffusion filtering results of reflected wave separation result profile

5 結論

聲波遠探測作為一種有效的測井手段，在井外地質(zhì)體探測和隱蔽油藏識別等方面有重要的作用。本文將圖像處理中的各向異性擴散濾波技術引入聲波遠探測的圖像處理中，經(jīng)過對比將相干增強張量模型和顯示擴散模型相結合，可以有效抑制深部地層中反射波圖像中的非相干噪聲，同時保持圖像的地質(zhì)結構信息，突出裂縫的邊界特征。處理后的資料品質(zhì)有較為明顯的提高，對細微裂縫的分辨率更高，更有利于對小裂縫和斷裂的識別，證實了各向異性擴散濾波在反射波圖像處理中的有效性和實用性。各向異性擴散在強噪聲干擾下的圖像降噪方面可能會出現(xiàn)將有效信息誤認為噪聲的情況，這方面還有待進一步的完善和提高。為了進一步突顯有效地質(zhì)信息，可以將相干體技術和各向異性擴散濾波結合起來，進行屬性提取分析，有望更好地表示井外反射體信息。