基于非二次冪Curvelet變換的最小二乘匹配算法及其應(yīng)用

袁艷華，王一博，劉伊克，常 旭

中國(guó)科學(xué)院地質(zhì)與地球物理研究所，北京 100029

1 引 言

Curvelet變換于1999年由Candès［1］等人首次提出，此后被廣泛應(yīng)用于圖像處理領(lǐng)域［2－5］．在2004年的SEG 年會(huì)上，Herrmann等做了關(guān)于Curvelet域多次波衰減、基于Curvelet變換和小波變換的AVO 反 演等報(bào)告［6－7］，自此Curvelet變換開(kāi)始被 廣泛應(yīng)用于地震數(shù)據(jù)處理中［8－18］，其中信噪分離是Curvelet變換一個(gè)重要的應(yīng)用方向［13－18］．一般而言，信號(hào)和噪聲可以在諸如波長(zhǎng)、能量、位置、方向因素的一個(gè)或幾個(gè)方面加以區(qū)分，這些因素在Curvelet分解中都是非常重要的參數(shù)，經(jīng)Curvelet變換后，各個(gè)因素均會(huì)在Curvelet域得到明顯體現(xiàn)．Curvelet變換在信噪分離上的獨(dú)特優(yōu)勢(shì)，為預(yù)測(cè)匹配濾波方法提供了強(qiáng)有力的保證．

Curvelet變換根據(jù)最大分解尺度及低頻成分的方向分辨率等參數(shù)進(jìn)行自動(dòng)二進(jìn)分解，具體來(lái)講，每增加兩個(gè)分解尺度，角度方向網(wǎng)格加密一倍．該分解模式的優(yōu)點(diǎn)在于簡(jiǎn)單、易于理解和實(shí)現(xiàn)．對(duì)于能量在各個(gè)頻率區(qū)間和方向上分布比較均勻的圖像或數(shù)據(jù)而言，Curvelet變換是一種較為理想的分析和處理工具．然而，在分析那些對(duì)于特定頻率范圍或方向上有明顯偏向的數(shù)據(jù)時(shí)，經(jīng)典二進(jìn)Curvelet變換遇到了瓶頸，地震數(shù)據(jù)就是個(gè)比較典型的例子．震源性質(zhì)、地應(yīng)力方向、地層對(duì)不同頻率和方向地震波的衰減作用、檢波器在接收頻段和方向上的設(shè)計(jì)等因素，導(dǎo)致我們難以得到在各個(gè)頻率和方向上能量分布均衡的地震數(shù)據(jù)．經(jīng)典Curvelet變換可以高效率地對(duì)具有一般特征的數(shù)據(jù)或圖像進(jìn)行自動(dòng)二進(jìn)分解，但難以做到針對(duì)地震數(shù)據(jù)進(jìn)行自適應(yīng)分析．針對(duì)這個(gè)問(wèn)題，本文采用了一種改進(jìn)的Curvelet變換，即自適應(yīng)的非二次冪Curvelet變換［14］．該變換的主要思想在于根據(jù)數(shù)據(jù)本身的特征和先驗(yàn)信息，設(shè)計(jì)靈活的、有針對(duì)性的濾波窗函數(shù)，并根據(jù)各區(qū)間方向分辨率進(jìn)行自動(dòng)剖分．該方法繼承了經(jīng)典Curvelet變換在高維數(shù)據(jù)分析上的優(yōu)勢(shì)，并彌補(bǔ)了其在處理地震資料上的局限性．就效果而言，非二次冪變換方法靈活性高、針對(duì)性強(qiáng)、分解結(jié)果分辨率高；就效率而言，相對(duì)于全局的二次冪分解，局部的非二次冪算法可以節(jié)省對(duì)信息量稀疏區(qū)域的繁瑣處理時(shí)間，當(dāng)然，效率提高的程度取決于數(shù)據(jù)本身的特征．

基于非二次冪Curvelet變換，本文提出了一種改進(jìn)的最小二乘匹配濾波算法（Non－Dyadic Curvelet based Least－Squares Method，NDCLSM）．該方法首先將輸入數(shù)據(jù)和預(yù)測(cè)的噪聲模型變換到非二次冪Curvelet域，然后在噪聲能量集中的非二次冪Curvelet子記錄上進(jìn)行最小二乘匹配濾波．該方法具有非二次冪Curvelet變換多尺度、多方向、根據(jù)數(shù)據(jù)特征進(jìn)行自適應(yīng)參數(shù)選擇等特點(diǎn)，適用性強(qiáng)，可廣泛應(yīng)用于諸如面波［18］、層間多次波和自由表面多次波［19－25］等的自適應(yīng)匹配相減．本文對(duì)非二次冪Curvelet變換和NDCLSM 算法的原理進(jìn)行了介紹，并使用三維實(shí)際地震資料驗(yàn)證了本文方法的有效性．

2 非二次冪Curvelet變換

Curvelet變換是一種多尺度、多角度的時(shí)頻分析方法，通常在頻率域內(nèi)實(shí)現(xiàn)，Curvelet正變換可表示為：

其中Uj，l（w）為尺度j角度l定義下的Curvelet基函數(shù)對(duì)應(yīng)的頻率域窗函數(shù)，通常Uj，l（w）的支撐區(qū)域是由Wj（w）（徑向窗）和Vl（w）（角度窗）的支撐區(qū)間聯(lián)合限定的楔形區(qū)域．為Curvelet域坐標(biāo)為k＝（k1，k2）的Curvelet系數(shù)在時(shí)間空間域內(nèi)對(duì)應(yīng)的位置，其中Rθl＝

Curvelet變換基于二進(jìn)的尺度角度分解，尺度窗函數(shù)定義為：Wj（w）＝W（2－jw）；角度窗函數(shù)定義為：Vj（w）＝V（2｜＿j/2＿｜w2/w1），其中為求整符號(hào)．

圖1顯示了Curvelet變換的多分辨率、多角度和空間定位等特性，與小波變換等變換方法相比，Curvelet變換更適于高維信號(hào)處理，特別是當(dāng)信號(hào)具有交叉和重疊等現(xiàn)象的時(shí)候．

圖1 幾個(gè)Curvelet函數(shù)在時(shí)間空間域（左）和頻率波數(shù)域（右）的顯示，表明Curvelet具有多尺度、多方向、空間頻率定位特性，以及滿足拋物尺度關(guān)系．Fig．1 Several Curvelets in time－spatial domain（left）and frequency－wavenumber domain（right）indicative of Curvelet′s multi－scale，multi－directional and localized in T－X and F－K domains，satisfying parabolic scaling relationship．

非二次冪Curvelet變換基于目標(biāo)特征進(jìn)行自適應(yīng)的尺度角度分解．假設(shè)目標(biāo)區(qū)域在頻率域可定義為：（w1，w2）∈Ej，則可引入

引入

對(duì)應(yīng)于目標(biāo)區(qū)域Ej的方向范圍，假設(shè)目標(biāo)Ej被分解為L(zhǎng)j個(gè)方向，則該區(qū)域的方向分辨率為．定義：

則頻率域窗函數(shù)為：Uj，l∶＝Wj×Vjl．對(duì)非二次冪Curvelet變換的細(xì)節(jié)描述可參考文獻(xiàn)［14］．

經(jīng)典Curvelet變換的分解方式未能充分考慮數(shù)據(jù)本身特點(diǎn)，從而影響了變換的效果和效率，非二次冪Curvelet變換繼承了Curvelet變換的優(yōu)點(diǎn)，并克服了其局限性和非自適應(yīng)性．這里我們分別使用兩種方法對(duì)一個(gè)“米”字形模型進(jìn)行分解，如圖2a所示，該模型包含四個(gè)方向的分量：水平、垂直和兩個(gè)對(duì)角．圖2b為二次冪Curvelet變換結(jié)果，從分解結(jié)果可以看出，第三、第四尺度均得到合適的分解，但第二尺度方向分解不足，未將交叉分量區(qū)分開(kāi)．而增加第二尺度角度分解總數(shù)則會(huì)導(dǎo)致高頻分量分解過(guò)量而模糊化．圖3為非二次冪Curvelet變換的分解結(jié)果，各個(gè)尺度、各個(gè)方向均得到恰當(dāng)剖分．

3 NDCLSM 算法

自適應(yīng)匹配濾波的目的是使去噪后的殘差最小化，而使殘差的L2范數(shù)最小的自適應(yīng)濾波稱(chēng)為最小二乘匹配濾波（LSM）．最小二乘自適應(yīng)匹配濾波過(guò)程可以表示為：

其中d為含噪地震記錄，m為預(yù)測(cè)的噪聲干擾模型，G為自適應(yīng)濾波因子．LSM 在時(shí)間空間域得到了廣泛的應(yīng)用，但是該算法要求有效信號(hào)與噪聲互不相關(guān)，而實(shí)際地震數(shù)據(jù)很難符合這個(gè)要求．以面波為例，一般情況下面波干擾分布范圍廣，能量強(qiáng)，與有效信號(hào)重疊嚴(yán)重，因此如果直接在時(shí)間空間域進(jìn)行最小二乘匹配，會(huì)損傷有效信號(hào)并導(dǎo)致面波能量殘留．而非二次冪Curvelet變換的高度稀疏性和區(qū)分信噪的能力使得基于這種變換的最小二乘匹配算法能在充分壓制噪聲的同時(shí)，有效保留反射信號(hào)，并大幅提高算法的穩(wěn)定性和準(zhǔn)確性．

NDCLSM 的算法流程可以歸納為以下幾點(diǎn)：通過(guò)分析輸入數(shù)據(jù)和噪聲干擾的特征，對(duì)記錄及預(yù)測(cè)的噪聲干擾進(jìn)行非二次冪Curvelet變換；然后分析噪聲干擾能量主要集中的子域，并實(shí)施最小二乘意義下的匹配算法，最小化噪聲干擾和實(shí)際資料中噪聲干擾所對(duì)應(yīng)的非二次冪Curvelet系數(shù)的殘差能量．

在非二次冪Curvelet域的最小二乘匹配濾波可以表示為：

其中dc為實(shí)際地震記錄對(duì)應(yīng)的非二次冪Curvelet系數(shù)，mc為預(yù)測(cè)的噪聲干擾記錄對(duì)應(yīng)的非二次冪Curvelet系數(shù)，Gc為將預(yù)測(cè)噪聲干擾對(duì)應(yīng)的非二次冪Curvelet系數(shù)進(jìn)行旋轉(zhuǎn)、尺度化和平移化的變換或?yàn)V波矩陣．求解（7）式殘差的L2范數(shù)最小可以得到最小二乘意義下的最佳濾波因子Gc．為了保證算法的穩(wěn)定性，我們引入阻尼因子ε2，并同時(shí)使能量殘差和非二次冪Curvelet系數(shù)達(dá)到最小，即

上述目標(biāo)函數(shù)可以使用共軛梯度等最優(yōu)化方法求解．在噪聲干擾能量集中的非二次冪Curvelet子記錄上使用最小二乘算法得到對(duì)應(yīng)的濾波因子及面波的非二次冪Curvelet系數(shù)cmc，然后通過(guò)非二次冪Curvelet變換的伴隨算子CT，得到匹配后的面波記錄m＝CTcmc．

4 三維實(shí)際地震資料的面波提取

在實(shí)際地震資料噪聲處理中，NDCLSM 算法可廣泛用于諸如面波、表面多次波和層間多次波等的自適應(yīng)匹配相減，本文將該方法應(yīng)用于面波提?。娌ㄊ顷懮系卣鹂碧街凶畛Ｒ?jiàn)的噪聲干擾之一，特別是在中國(guó)西部山區(qū)、高原地區(qū)、沙漠地區(qū)等廣泛存在．面波能量強(qiáng)，大范圍地覆蓋在有效信號(hào)之上，通常存在頻散現(xiàn)象，嚴(yán)重影響著后續(xù)的處理和解釋工作．傳統(tǒng)的去噪方法主要利用面波頻率低、視速度小的特點(diǎn)，很難做到在消除面波的同時(shí)保護(hù)有效信號(hào)，并且難以保存對(duì)后續(xù)地質(zhì)構(gòu)造研究和解釋極為重要的低頻數(shù)據(jù)．此外，在天然地震學(xué)研究中，有效提取面波成分對(duì)反演S波速度、解釋地殼和上地幔結(jié)構(gòu)具有重要意義．

為了測(cè)試本文方法在面波提取中的效果，我們對(duì)實(shí)際三維地震資料進(jìn)行了處理．圖4a為某地區(qū)三維地震數(shù)據(jù)中的一個(gè)單炮記錄，可以看出該數(shù)據(jù)中面波能量非常強(qiáng)，覆蓋在有效反射信號(hào)同相軸之上，使得反射信號(hào)中斷、模糊甚至畸變，這種數(shù)據(jù)無(wú)法直接用于后續(xù)的處理和解釋?zhuān)覀兲崛?shù)據(jù)中低頻低視速度的成分作為初始面波模型，如圖4b所示．圖4c為本文方法壓制面波干擾后的結(jié)果．可以看出面波成分壓制得比較徹底，保證了反射信號(hào)的完整性，使隱藏在面波能量下的反射信號(hào)清晰可見(jiàn)，并與周?chē)耐噍S較為自然地連為一體．

圖4 對(duì)一個(gè)共炮點(diǎn)道集記錄的處理結(jié)果（a）一個(gè)帶有面波干擾的炮集；（b）預(yù)測(cè)的面波模型；（c）應(yīng)用本文方法得到的處理結(jié)果．Fig．4 Processing result of a common shot record（a）Shoot gather with surface waves interferences；（b）Predicted ground－roll model；（c）Result using the proposed filtering method．

由于NDCLSM 算法效率較高，我們可以實(shí)現(xiàn)對(duì)大量數(shù)據(jù)的快速處理．圖5a為一個(gè)三維數(shù)據(jù)體的部分顯示，我們對(duì)該數(shù)據(jù)體的單炮記錄重復(fù)圖4中的去噪過(guò)程．可以看出，處理前面波幾乎覆蓋了整個(gè)數(shù)據(jù)體，有效信號(hào)難以識(shí)別，處理后的結(jié)果如圖5b所示，面波能量得到了有效去除，有效反射信號(hào)變得清晰、自然而且連續(xù)，表明了本文方法在處理大數(shù)據(jù)量實(shí)際資料中的有效性．由于面波能量隨著距離的增大而迅速衰減，因此其能量主要集中在近偏移距道集．為了更清楚的顯示去噪效果，我們將數(shù)據(jù)按共偏移距道集排列，并顯示近偏移距部分，如圖6a所示，圖6b為圖6a對(duì)應(yīng)的處理結(jié)果．圖6c—6d為圖6a中矩形區(qū)域處理前和處理后的顯示結(jié)果．

圖5 三維實(shí)際資料的處理結(jié)果（a）處理前的記錄；（b）應(yīng)用本文方法處理后的結(jié)果．Fig．5 Processing results of a 3－D field data（a）Records before filtering；（b）Result after filtering by the proposed method．

圖6 對(duì)三維數(shù)據(jù)體共偏移距道集近偏移距部分的處理結(jié)果（a）處理前記錄；（b）應(yīng)用本文方法處理后的結(jié)果；（c）（a）中矩形區(qū)域的放大顯示；（d）對(duì)（c）進(jìn)行去噪處理的結(jié)果．Fig．6 Processing results of near offset parts of 3Dcommon offset gather（a）Data before filtering；（b）Result filtering by the proposed method；（c）Zoom view of the rectangular part in（a）；（d）Result of de－noise processing to（c）．

5 結(jié) 論

通過(guò)分析地震數(shù)據(jù)的頻率和方向特征，本文對(duì)地震數(shù)據(jù)進(jìn)行自適應(yīng)非二次冪Curvelet分解，最大程度地將噪聲干擾和有效反射信號(hào)分離開(kāi)來(lái)，在此基礎(chǔ)上，結(jié)合最小二乘匹配濾波算法，實(shí)現(xiàn)了自適應(yīng)信噪分離，提高了常規(guī)最小二乘匹配算法的穩(wěn)定性和匹配結(jié)果的準(zhǔn)確度．三維實(shí)際數(shù)據(jù)的處理結(jié)果驗(yàn)證了本文方法的可行性和有效性．

致 謝 感謝兩位匿名評(píng)審人提出的寶貴意見(jiàn)．

（References）

［1］ Candès E J，Donoho D L．Curvelets：A Surprisingly Effective Nonadaptive Representation for Objects with Edges，Curves and Surfaces．TN：Vanderbilt University Press，1999．

［2］ Candès E J，Guo F．New multiscale transforms，minimum total variation synthesis：Applications to edge－preserving image reconstruction．SignalProcessing，2002，82（11）：1519－1543．

［3］ Starck J，Murtagh F，Candes E，Donoho D．Gray and color image contrast enhancement by the Curvelet transform．IEEETransactionsonImageProcessing，12：706－717．

［4］ Candès E J，Donoho D L．Curvelets， multiresolution representation，and scaling laws．Wavelet Applications in Signal and Image Processing VIII，A．Aldroubi，A．F．Laine，M．A．Unser eds．，Proc．SPIE 4119．

［5］ Starck J L，Candès E J，Donoho D L．The Curvelet transform for image denoising．IEEETrans．Image Processing，2002，11（6）：670－684．

［6］ Herrmann F J，Verschuur E．Curvelet－domain multiple elimination with sparseness constraints．SocietyofExploration Geophysicists，ExpandedAbstracts，2004，23（1）：1333－1336．

［7］ Hennenfent G，Herrmann F J．Three－term amplitude－versusoffset（AVO）inverse on revisited by Curvelet and wavelet transforms．SocietyofExplorationGeophysicists，ExpandedAbstracts，2004，23：211－214．

［8］ Chauris H，Nguyen T．Seismic demigration/migration in the Curvelet domain．Geophysics，2008，73（2）：S35－S46．

［9］ Douma H，de Hoop M．Leading－order seismic imaging using Curvelets．Geophysics，2007，72（6）：S231－S248．

［10］ Hermann F J，Wang D，Hennenfent G，et al．Curvelet－based seismic data processing：A multiscale and nonlinear approach．Geophysics，2008，73（1）：706－716．

［11］ Herrmann F G，Hennenfent G．Non－parametric seismic data recovery with Curvelet frames．Geophys．J．Int．，2008，173（1）：233－248．

［12］ Naghizadeh M，Sacchi M D．Beyond alias hierarchical scale Curvelet interpolation of regularly and irregularly sampled seismic data．Geophysics，2010，75（6）：WB189－WB202．

［13］ Neelamani R，Baumstein A I，Gillard D G，et al．Coherent and random noise attenuation using the Curvelet transform．TheLeadingEdge，2008，27（2）：240－248．

［14］ 袁艷華，王一博，劉伊克，常旭．非二次冪Curvelet變換及其在地震噪聲壓制中的應(yīng)用．地球物理學(xué)報(bào)，2013，56（3）：1023－1032．

Yuan Y H，Wang Y B，Liu Y K，Chang X．Non－dyadic Curvelet transform and its application in seismic noise elimination．ChineseJ．Geophys．（in Chinese），2013，56（3）：1023－1032．

［15］ Wang Y B，Yuan Y H，Liu Y K，Chang X．Geophysical data oriented Curvelet－like transform．SEGAbstract，2011：3663－3667．

［16］ Yuan Y H，Wang Y B，Liu Y K，Chang X．Curvelet domain adaptive least－squares subtraction of internal multiple．JournalofSeismicExploration，2011，20（3）：273－288．

［17］ 包乾宗，高靜懷，陳文超．Curvelet域垂直地震剖面波場(chǎng)分離．西安交通大學(xué)學(xué)報(bào)，2007，41（6）：650－654．

Bao Q Z，Gao J H，Chen W C．Wave field separation of vertical seismic profiling in Curvelet domain．JournalofXi′anJiaotongUniversity（in Chinese），2007，41（6）：650－654．

［18］ Wang Y B，Xue Y W，Dong S Q．Surface waves suppression using interferometric prediction and Curvelet domain hybrid L1/L2norm subtraction．SEGExpandedAbstracts，2009，28：3292．

［19］ Verschuur D J，Berkhout A J，Wapenaar C P A．Adaptive surface－related multiple elimination．Geophysics，1992，57（9）：1166－1177．

［20］ Spitz S． Pattern recognition，spatial predictability and subtraction of multiple events．TheLeadingEdge，1999，18（1）：55－58．

［21］ Guitton A． Multiple attenuation with multidimensional prediction－error filter．73th Ann．Internat．Mtg．，Soc．Expl．Geophys．，Expanded Abstracts，2003：57－74．

［22］ 金德剛，常旭，劉伊克．逆散射級(jí)數(shù)法預(yù)測(cè)層間多次波的算法改進(jìn)及其策略．地球物理學(xué)報(bào)，2008，51（4）：1209－1217．

Jin D G，Chang X，Liu Y K．Algorithm improvement and strategy of internal multiples prediction based on inverse scattering series method．ChineseJ．Geophys．（in Chinese），2008，51（4）：1209－1217．

［23］ Lin D，Young D，Huang Y，et al．3－D SRME application in the Gulf of Mexico．74th Ann．Internat．Mtg．，Soc．Expl．Geophys．，Expanded abstracts，2004，23：1257－1260．

［24］ Berkhout A J，Verschuur D J．Estimation of multiple scattering by iterative inversion，part I：theoretical considerations．Geophysics，1997，62（5）：1586－1595．

［25］ Weglein A B，Carvalho F A，Stolt P M．An inverse scattering series method for attenuating multiples in seismic reflection data．Geophysics，1997，62（6）：1975－1989．