基于能量比的Curvelet閾值迭代面波壓制

李繼偉 劉曉兵 周俊驊 簡世凱 曾 珍

(①東方地球物理公司西南物探研究院，四川成都 610000； ②中國石油西南油氣田勘探開發(fā)研究院，四川成都 610000； ③中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580)

0 引言

在實(shí)際地震勘探中，面波是影響原始地震數(shù)據(jù)信噪比的重要因素。在道集記錄上主要表現(xiàn)為頻率低、速度低、能量強(qiáng)，呈掃帚狀發(fā)散分布。針對(duì)面波的壓制，現(xiàn)今已研發(fā)出多種成熟的配套技術(shù)，包括頻率域?yàn)V波、F-K濾波、F-X濾波、小波變換、Curvelet變換等[1-5]?？偟陌l(fā)展趨勢(shì)是由單一特征差異向多特征差異發(fā)展，由頻域向時(shí)頻域發(fā)展。

在頻率域，主要依據(jù)面波與有效波的頻率差異壓制面波，這是單一特征差異的面波壓制方法。F-K濾波、F-X濾波相對(duì)于頻率域?yàn)V波而言，增加了一項(xiàng)差異特征，如速度差異、空間位置差異等。頻率域?yàn)V波、F-K濾波和F-X濾波都具有“頻率時(shí)不變”的共同特性。小波變換是一種多尺度分析方法，它在將數(shù)據(jù)變換到頻率域的同時(shí)還保留了時(shí)間特性，可呈現(xiàn)信號(hào)在各尺度下的頻率特征。由于小波變換角度分辨率不高[6]，它不能很好地表達(dá)圖像邊緣的方向特征。

Curvelet變換是一種多尺度、多方向的分析方法。 1999年，Candès等[7]基于Rideglet變換提出了第一代Curvelet變換。第一代Curvelet變換的數(shù)字實(shí)現(xiàn)比較繁瑣，且具有很大的數(shù)據(jù)冗余。隨后，Candès等[8-10]進(jìn)一步提出了一種改進(jìn)的快速離散傅里葉變換方法，即第二代Curvelet變換。Curvelet變換具有較強(qiáng)的方向特征，針對(duì)高維信號(hào)的線奇異特征具有最優(yōu)的非線性表達(dá)，能克服小波變換對(duì)于圖像邊緣方向特征表達(dá)不足的欠缺。

Curvelet變換的多尺度、多方向特性被廣泛應(yīng)用于地震資料去噪，主要包括隨機(jī)噪聲壓制、面波分離、多次波預(yù)測(cè)等[11-18]。張恒磊等[13]通過Curvelet多尺度分解，在曲線變化特征良好逼近的基礎(chǔ)上采用非線性閾值壓制隨機(jī)噪聲，取得了良好效果。董烈乾等[15]利用面波與有效波在Curvelet域不重疊的特性，對(duì)Curvelet系數(shù)進(jìn)行處理，在Curvelet域成功地壓制了面波； 他們還提出利用曲波多方向特性，自動(dòng)選取遵循斯奈爾定律的多次波貢獻(xiàn)道集，對(duì)多次波模型道進(jìn)行構(gòu)建，彌補(bǔ)了自由表面多次波壓制方法(SRME)預(yù)測(cè)多次波時(shí)產(chǎn)生的振幅誤差。

針對(duì)基于Curvelet變換壓制面波的方法，已有許多學(xué)者提出了多種不同的閾值計(jì)算公式[19-21]，且具有各自的優(yōu)缺點(diǎn)。但其壓制效果都受到面波與有效波在Curvelet域的重疊程度影響。鑒于此，本文提出一種基于能量比值的Curvelet閾值迭代壓制面波方法，通過將面波與有效波進(jìn)行多次閾值迭代，可以更徹底地分離面波與有效波。

1 基于二代Curvelet變換的面波壓制

1.1 第二代Curvelet變換原理

Curvelet變換屬于多尺度分析理論范疇，其實(shí)質(zhì)是設(shè)定一個(gè)基函數(shù)與信號(hào)進(jìn)行內(nèi)積，實(shí)現(xiàn)信號(hào)的多尺度表達(dá)。假設(shè)一個(gè)二維空間，其空間域參數(shù)為x，頻率域參數(shù)為w，r、θ為極坐標(biāo)參數(shù)，l為方向參數(shù)。

(1)

(2)

對(duì)于每一個(gè)尺度參數(shù)j≥j0(j0為起始值)，都可利用傅里葉變換定義一個(gè)頻率窗

(3)

定義Curvelet“母”函數(shù)為Φj(x),通過傅里葉變換，可得

(4)

則其他2-j尺度的Curvelet函數(shù)都可由該母函數(shù)Φj(x)經(jīng)過旋轉(zhuǎn)、平移得到。

(5)

(6)

由此可知，給定一個(gè)Curvelet函數(shù)，經(jīng)過伸縮、平移和旋轉(zhuǎn)可生成L2(R2)(平方可積函數(shù)空間)的緊標(biāo)架，這意味著它具有重構(gòu)公式

(7)

式中：Φj,l,k表示由指標(biāo)(j,l,k)確定的Curvelet簇；M表示指標(biāo)。

圖1是Curvelet變換頻率域尺度分割示意圖，圖中陰影區(qū)為楔形，每個(gè)楔形對(duì)應(yīng)一個(gè)Curvelet基函數(shù)。由Curvelet定義可知，每個(gè)Curvelet基函數(shù)有一個(gè)特定的尺度、方向和位置，這些特定因素可由式(6)中的下標(biāo)(j,l,k)表征。

圖1 Curvelet變換頻率域分割示意圖

整個(gè)Curvelet變換按照頻率由低到高劃分尺度，由最內(nèi)層粗尺度(Coarse)、中間層精細(xì)尺度(Detail)和最外層最佳尺度(Fine)組成。每個(gè)尺度沿著順時(shí)針劃分方向，每個(gè)方向由相應(yīng)頻率的Curvelet系數(shù)組成。其中： 粗尺度包含了低頻系數(shù)信息，不具有方向性； 最佳尺度包含了高頻系數(shù)信息，具有方向特性； 精細(xì)尺度主要用于方向分割，具有方向特性。

在利用Curvelet變換進(jìn)行信號(hào)處理過程中，通常會(huì)將信號(hào)進(jìn)行尺度和方向劃分。其中粗尺度與最佳尺度下的數(shù)據(jù)包含了信號(hào)的原本概貌，不作為處理對(duì)象，而精細(xì)尺度下的數(shù)據(jù)包含了待處理的信息，是主要研究對(duì)象。

1.2 Curvelet域面波壓制原理

面波具有頻率低、速度低、能量強(qiáng)，在道集記錄中呈掃帚狀發(fā)散分布等特點(diǎn)，根據(jù)面波與有效波在頻率、速度、方向上的特征差異，利用Curvelet變換的多尺度、多方向等特性，將面波與有效波變換到Curvelet域，可更細(xì)致地描述面波與有效波的差異。與其他傳統(tǒng)面波壓制方法相比，Curvelet變換的多尺度、多方向特性具有明顯優(yōu)勢(shì)。

假設(shè)一個(gè)地震數(shù)據(jù)體由下式表示

S=x1+x2+n

(8)

式中：x1為有效信號(hào)；x2為面波；n為高斯白噪聲。

對(duì)地震數(shù)據(jù)體做Curvelet變換，可將x1、x2和n用Curvelet系數(shù)表示

(9)

式中：CT表示Curvelet變換；C1、C2、C3對(duì)應(yīng)有效波、面波及高斯白噪聲的Curvelet系數(shù)。

由式(9)可知，通過合適的尺度、方向等變換參數(shù)，可將有效波與面波變換到Curvelet域中，將面波系數(shù)置零，實(shí)現(xiàn)面波與有效波的分離。

2 基于能量比值的閾值迭代方法

Curvelet變換壓制面波的效果受面波與有效波在Curvelet域中的重疊程度影響，兩者重疊程度越低，壓制效果越好。在實(shí)際地震采集資料中，因不同的地質(zhì)環(huán)境、采集工藝，造成面波分布特征不盡相同。針對(duì)面波與有效波重疊程度較高的資料，Curvelet變換也難以較徹底地分離面波與有效波。本文嘗試提出一種基于能量比的閾值迭代壓制面波方法：將面波與有效波在Curvelet域進(jìn)行不同頻帶的多次分解；在不同尺度和方向上，根據(jù)面波與有效波的分布特征，以不同的閾值進(jìn)行多次迭代，力求更徹底地分離面波與有效波。

2.1 能量比值

在信號(hào)處理過程中，分析兩信號(hào)的不同點(diǎn)及相同點(diǎn)時(shí)，通常采用對(duì)比法，如能量對(duì)比法等。本文基于能量對(duì)比法，通過面波模型與地震數(shù)據(jù)體中的面波在Curvelet域中進(jìn)行能量對(duì)比，根據(jù)相似程度，將地震數(shù)據(jù)體中的面波分離出來。

假設(shè)面波x2已知，將地震數(shù)據(jù)體S與面波x2進(jìn)行相同參數(shù)的Curvelet變換，得到CT(S)和CT(x2)。

將相同尺度和方向下的CT(S)與CT(x2)進(jìn)行能量對(duì)比

(10)

式中A(j,l)1和A(j,l)2分別為地震數(shù)據(jù)體S和面波x2在j尺度、l方向分量的振幅。

由式(10)可知，通過分析能量比值ζ的大小，便可知地震數(shù)據(jù)體S在j尺度和l方向分量下的面波與面波模型x2在相應(yīng)尺度和方向分量下的面波的相似程度。其中，ζ越大，表明地震數(shù)據(jù)體在該尺度和方向分量下的面波與面波模型在相應(yīng)尺度和方向分量下的面波的相似性越小，表征面波含量越小； 反之越大。

2.2 閾值選取

Curvelet變換壓制面波，很重要的一個(gè)環(huán)節(jié)閾值處理。在閾值去噪中，通常分為軟閾值處理和硬閾值處理兩種。由于硬閾值處理存在間斷點(diǎn)，容易造成重構(gòu)失真現(xiàn)象，本文采用軟閾值處理方法。

由式(10)可知，根據(jù)不同尺度j和方向l下ζ的大小，采用基于能量比值的軟閾值面波壓制方法。重構(gòu)系數(shù)設(shè)定為

(11)

2.3 迭代流程

上述處理方法，只要給定一個(gè)與地震數(shù)據(jù)體中的面波相似的含面波數(shù)據(jù)體，便可由式(11)分離出等同或者更多、更少的面波。

圖2為處理流程圖，可細(xì)分為以下具體步驟：

(1)將原始地震數(shù)據(jù)體進(jìn)行Curvelet變換；

(2)選取含有面波的尺度和方向分量做Curvelet反變換，得到預(yù)估面波模型；

(3)將預(yù)估面波模型與原始地震數(shù)據(jù)做相同參數(shù)的Curvelet變換，并進(jìn)行閾值處理；

(4)通過多次設(shè)置不同的Curvelet變換參數(shù)(尺度、方向等)、尺度頻率范圍，根據(jù)面波與有效波的分布特征，針對(duì)不同區(qū)域設(shè)置不同的閾值進(jìn)行多次迭代，直至達(dá)到滿意效果。

圖2 面波壓制迭代處理流程

3 實(shí)際應(yīng)用

實(shí)際資料來源于四川盆地M工區(qū)(圖3)。主要采集參數(shù)分別為：道距30m，接收道數(shù)180，采樣間隔2ms，中點(diǎn)炸藥激發(fā)，記錄長度5s，本文記錄顯示長度為3.5s。

首先對(duì)采集數(shù)據(jù)中的面波做速度分析和頻率分析。通過速度掃描得知，面波視速度最高約達(dá)2100m/s。圖4為采集數(shù)據(jù)的頻率掃描結(jié)果，采用高通濾波器掃描，掃描間隔為8Hz。圖中頻率掃描到24Hz時(shí)仍能看到部分面波(圖4c紅色直線下方區(qū)域)，該部分面波能量強(qiáng)，線性相關(guān)性低，與有效波重疊程度高。但在實(shí)際應(yīng)用中，利用傳統(tǒng)方法壓制這部分面波時(shí)，由于其線性相關(guān)性低，且頻率高，往往難以達(dá)到理想壓制效果。

針對(duì)該區(qū)面波分布特征，基于第二代Curvelet變換，采用能量比值的閾值迭代方法進(jìn)行壓制。即首先對(duì)原始數(shù)據(jù)做Curvelet變換，得到各尺度下T-X域的地震數(shù)據(jù)體(圖5)。尺度分析所采用的頻帶范圍是0～40Hz，尺度依次為1～7。

從圖5可見：尺度1～3主要是一些超低頻信息，這部分信息屬于低頻系數(shù)部分，涵蓋了原始數(shù)據(jù)的概貌(圖5a～圖5c)；尺度4～7中都含有面波(圖5d～圖5g中紅色直線下方區(qū)域)。含有面波的尺度是重點(diǎn)研究對(duì)象，需進(jìn)行方向細(xì)分。

圖3 四川盆地某工區(qū)采集數(shù)據(jù)

圖4 采集數(shù)據(jù)不同低截頻率時(shí)的頻率掃描結(jié)果

圖5 原始數(shù)據(jù)在各尺度下的T-X域數(shù)據(jù)體

將尺度4～7下的數(shù)據(jù)體再做方向細(xì)分，根據(jù)面波分布的視速度范圍，確定面波在Curvelet域所處角度方向，選取含有面波的角度范圍做Curvelet反變換，所得數(shù)據(jù)體作為預(yù)估面波模型。

根據(jù)工區(qū)內(nèi)面波分布特點(diǎn)，遵循先強(qiáng)后弱、先低頻后高頻、先低速后高速的原則，對(duì)原始地震數(shù)據(jù)體進(jìn)行了3次閾值迭代。此3次Curvelet變換參數(shù)如表1所示。

圖6是多次迭代后的道集記錄對(duì)比?？梢娫疾杉瘮?shù)據(jù)經(jīng)3次閾值迭代后面波含量逐漸減少，同相軸連續(xù)性增強(qiáng)； 被面波掩蓋的有效波在3次迭代后，逐漸突顯出來，信噪比得到顯著提高。

表1 多次迭代的Curvelet變換參數(shù)表

圖6 多次迭代后的道集記錄對(duì)比圖

圖7是多次迭代前、后的頻譜曲線對(duì)比?？梢姷皖l面波在每一次設(shè)定的尺度頻率范圍內(nèi)被多次衰減，中高頻有效波頻譜特征與原先基本一致，同時(shí)還保留了部分低頻信息。

圖7 多次迭代前、后的頻譜曲線對(duì)比圖

為了對(duì)比驗(yàn)證本文方法的可靠性和優(yōu)越性，分別對(duì)原始數(shù)據(jù)再做傳統(tǒng)Curvelet變換閾值面波壓制和F-K濾波，三種方法參數(shù)如表2所示。

表2 三種方法參數(shù)對(duì)照表

圖8是三種方法壓制面波后的道集記錄及噪聲記錄對(duì)比圖。從總體上看，三種方法都具有一定的面波壓制作用，且基于Curvelet變換的面波壓制方法(圖8b和圖8c)效果明顯優(yōu)于F-K濾波(圖8a)。

圖8a是利用傳統(tǒng)F-K濾波壓制面波后的道集及其噪聲記錄。可見面波壓制欠徹底，尚存較多低頻的殘留面波(紅線下方區(qū)域)。這部分強(qiáng)能量面波掩蓋了有效波，使波組特征表現(xiàn)不突出，嚴(yán)重影響了整個(gè)資料的信噪比及后期資料處理和解釋。

由于傳統(tǒng)F-K濾波需滿足線性相關(guān)的前提條件，該局限性使其對(duì)在道集記錄中表現(xiàn)為低頻、低線性相關(guān)、能量強(qiáng)，且與有效波重疊程度較高的這部分面波，難以達(dá)到理想效果。

圖8b是利用傳統(tǒng)Curvelet變換壓制面波后的道集及其噪聲記錄。可見在F-K濾波后(圖8a)仍殘留的面波也得到較好壓制(紅線下方區(qū)域)，其有效波較突顯，同相軸連續(xù)性變強(qiáng)，波組特征明顯，資料的整體信噪比顯著提高。充分利用Curvelet變換具有的多尺度、多方向特性，通過合適的尺度、方向等Curvelet變換參數(shù)，在Curvelet域?qū)⒌卣饠?shù)據(jù)體中有效波與低線性相關(guān)的殘留面波進(jìn)行有效分隔，進(jìn)而消除。

圖8c是基于本文方法壓制面波后的道集及其噪聲記錄。可見在圖8b(傳統(tǒng)Curvelet變換閾值法)道集記錄上仍然殘留的部分能量強(qiáng)、與有效波重疊程度高的面波被較徹底地消除(圖8c紅線下方區(qū)域)，該去噪后道集記錄信噪比更高，同相軸連續(xù)性更強(qiáng)，波組特征更明顯。

圖8 三種方法去噪后的道集(左)及其噪聲(右)記錄對(duì)比

圖9為圖8b和圖8c去噪后道集記錄的局部放大圖，對(duì)比可見圖9a的背景記錄上存在很多反向線性噪聲，即Curvelet系數(shù)腳印。因此，本文方法不僅可很好地分離面波與有效波，還能消除Curvelet變換遺留下來的Curvelet系數(shù)腳印。

圖10是三種方法壓制面波前、后的頻譜曲線對(duì)比，可見其對(duì)面波頻帶的壓制效果有所差異。基于本文的壓制面波方法，面波頻帶的壓制效果優(yōu)于傳統(tǒng)的面波壓制方法。

圖9 對(duì)應(yīng)圖8b(a)和圖8c(b)去噪記錄的局部放大

圖10 三種方法去噪前、后的頻譜曲線對(duì)比

為了驗(yàn)證本文方法應(yīng)用效果，將三種壓制面波方法處理后的剖面(圖11)進(jìn)行對(duì)比。從圖11中可以看出：利用F-K濾波方法壓制面波后，存在面波壓制不干凈，有效波被掩蓋的不足； 利用傳統(tǒng)Curvelet變換閾值法壓制面波后，剖面整體信噪比明顯提高，有效波出露；基于本文方法壓制面波后，剖面有效波連續(xù)性更強(qiáng)，信噪比更高。

4 結(jié)論

Curvelet變換是一種多尺度、多方向、多分辨率的分析方法，它將信號(hào)進(jìn)行多尺度分解，得到多方向上的Curvelet域信息； 通過分析面波與有效波的分布特征，可在Curvelet域有效分離面波與有效波。但是，Curvelet變換壓制面波受有效波與面波在Curvelet域重疊程度影響。本文提出的基于能量比值的閾值迭代壓制面波方法，將面波與有效波在Curvelet域進(jìn)行不同頻帶的多次分隔，在不同的尺度和方向上，根據(jù)面波與有效波的分布特征，采用不同的閾值進(jìn)行多次迭代，對(duì)面波的壓制達(dá)到了很好的效果。理論分析與實(shí)際應(yīng)用效果對(duì)比均表明，相對(duì)于傳統(tǒng)面波壓制方法，本文方法面波壓制更徹底，所獲資料信噪比更高，同相軸更連續(xù)。