基于多道卷積信號盲分離的多次波自適應相減方法

李鐘曉，陸文凱，龐廷華，王 季

清華大學自動化系智能技術與系統國家重點實驗室清華信息科學與技術國家實驗室（籌），北京 100084

基于多道卷積信號盲分離的多次波自適應相減方法

李鐘曉，陸文凱＊，龐廷華，王 季

清華大學自動化系智能技術與系統國家重點實驗室清華信息科學與技術國家實驗室（籌），北京 100084

本文將多次波自適應相減問題表示為一個多道卷積信號的盲分離問題．利用2D卷積核來表示預測多次波和實際多次波之間的差異，并采用分離出的一次波信號的非高斯性最大化作為優化目標，我們提出一種基于多道卷積信號盲分離的多次波自適應相減算法．為了求解上述非線性優化問題，所提方法將其轉化為一個迭代線性優化問題，采用迭代最小二乘方法加以實現．由于采用了多道卷積信號盲分離模型，所提方法能夠適應預測和真實多次波之間在時間及空間上的變化．通過對簡單模型數據、Pluto數據和實際數據進行處理，驗證了所提算法的有效性．

多次波自適應相減，盲信號分離，多道卷積信號，迭代最小二乘

1 引 言

多次波是地震數據處理，特別是海洋地震數據處理中最突出的問題之一．多次波與有效波相互混疊會干擾對有效波的正確識別，降低地震資料的信噪比和分辨率，影響地震偏移的精度．多次波壓制方法主要有濾波方法和基于波動方程的預測相減法［1］．濾波方法主要利用多次波的周期性和一次波與多次波之間的可分離性［2－5］；基于波動方程的預測相減法首先通過模擬實際波場或者反演地震數據估計出多次波，然后從原始地震信號中將預測多次波自適應減去，從而得到一次波［6－9］．因此，基于波動方程的預測相減法分為多次波預測和多次波自適應相減兩個步驟．其中，多次波自適應相減的主要方法有基于匹配濾波器設計的方法［10］、基于模式的方法［11－16］和基于盲信號分離的獨立成分分析的方法［17－21］．

基于匹配濾波器設計的多次波自適應相減方法，通過優化與一次波有關的目標函數，得到有限長的匹配濾波器，從而利用預測多次波逼近真實多次波，達到壓制多次波的目的．常用的準則有基于一次波能量最小化的方法［22］，此方法隱含了一次波和多次波的正交性假設．因此，在一次波和多次波相互正交時，能夠很好地壓制地震數據中的多次波；而在非正交時，則常常會損傷有效信號．因此，為了克服這一問題，基于增廣匹配濾波器的一次波L2范數最小化方法［23－24］、基于一次波L1范數最小化的方法［25－26］和非平穩正則化L2范數方法［27］等被相繼提出．

基于模式的多次波自適應相減方法主要利用地震信號的一次波和多次波各自具有不同的空間模式，比如一次波和多次波的預測誤差濾波器系數不同或者一次波和多次波的曲線擬合系數不同等等，以此作為求解一次波和多次波的約束條件，并構造與預測誤差或擬合誤差相關的優化目標，從而得到多次波壓制結果［11－16］．

基于盲信號分離的ICA方法［17－21］，將多次波自適應相減問題表示為多個信號瞬時混合的盲分離問題，利用幾何ICA的方法估計得到混合矩陣，并根據地震信號的超高斯性構造優化目標，實現一次波和多次波的分離．基于盲信號分離的ICA方法利用了地震信號的超高斯性，避免了基于一次波能量最小化方法需要采用的多次波和一次波正交的假設．需要指出的是，上述方法［17－21］的實現過程包括子波差異消除和ICA分離兩步，且所采用的混合ICA模型只考慮預測和實際多次波在時間方向的差異．

考慮到預測和實際多次波不只在時間上，而且在空間上也存在差異，本文提出了一種基于多道卷積信號盲分離的多次波自適應相減方法．在所提出的方法中，利用2D卷積模型，表示預測多次波和實際多次波之間的差異，并采用分離出的一次波信號的非高斯性最大化作為優化目標，估計2D卷積核系數．和傳統的基于L2范數的多道匹配濾波器相比，所提出的方法采用了非高斯性最大化準則，避免了多次波和一次波正交的假設．和基于瞬時混合ICA的多次波自適應相減方法不同，基于卷積信號盲分離的多次波自適應相減方法在實現一次波和多次波分離的過程中，同步消除預測和實際多次波在時間－空間上的差異．

下面首先介紹基于卷積信號盲分離的多次波自適應相減方法的基本原理，在給出優化目標的基礎上推導迭代求解算法，最后給出模型數據和實際數據的處理結果．

2 基本原理

基于濾波器設計的多次波自適應相減方法通常要對地震數據進行分塊處理，并假定同一分塊數據中的預測多次波與真實多次波之間，具有相同的尺度差異、子波差異、時延差異等，因此同一分塊中的數據只需要一個濾波器就能實現預測多次波和真實多次波之間的匹配．此時多次波自適應相減模型可以表示為［10］

其中，s為原始地震信號，p0為真實一次波，m是真實多次波，~m為預測多次波，＊表示卷積算子，h表示預測多次波與真實多次波之間差異的濾波器．這里，h可以是單道算子，也可以是多道算子．根據公式（1）表示的模型，多次波自適應相減問題就是如何設計匹配濾波器f來消除預測多次波與真實多次波之間的差異：

得到濾波器f后，一次波的估計p可以由原始地震信號、預測多次波和匹配濾波器獲得，即

其中，~M是由預測多次波構成的數據矩陣．

濾波器f設計的目標就是使預測多次波經過f濾波后與真實多次波之間的差異達到最小，從而使相減得到的一次波p的某個統計量最小．常用的基于能量最小化的多次波相減法就是假定一次波能量最小，即最小能量準則：

由式（4）得到f的最小二乘解

地震數據是超高斯分布的信號［28］，因此利用卷積信號盲分離中常用的非高斯性最大化準則來求取濾波器f和地震數據的統計特性更吻合．原始地震信號是由一次波和多次波相加得到的混合信號，而且在同一分塊中的一次波和多次波接近統計獨立且均符合超高斯分布［17］．根據概率論的中心極限定理，兩個滿足超高斯分布的相互獨立隨機變量的和，其分布比其中任何一個隨機變量的分布都要更接近高斯分布［29］．因此分塊數據中的原始地震信號減去多次波后，所得到的一次波具有最大的超高斯性．負熵是盲信號處理中一種常用的衡量非高斯性的準則，具有較高的魯棒性．根據負熵的定義，我們可建立如下的優化目標［29］：

3 求解算法

可利用多種方法對式（6）中的優化目標求解．我們借鑒文獻［4，30］中的求解思路，得到類似迭代加權最小二乘（IRLS）的求解算法．

附錄A中給出了單道濾波器求解的詳細推導過程，附錄B推導了多道濾波器的求解過程．針對多道匹配濾波器首先根據式（3）將濾波方程寫為［10］

其中，pi，j為濾波后得到的一次波估計，i＝1，2，…，T，T為時間方向分塊的大小，j＝1，2，…，S，S為分塊的最大道數，si，j為原始地震數據，m~i，j為相應的預測多次波數據，fk，l為多道濾波器系數，2K＋1為濾波器在時間方向的長度，2L＋1為濾波器在空間方向的道數．

將式（8）代入優化目標式（6）中并求導，最終可以得到如下方程（附錄B）：

其中，矩陣A中的第（k＋K）×（2L＋1）＋（l＋L＋1）行，第（k1＋K）×（2L＋1）＋（l1＋L＋1）列元素

列向量b中第（k＋K）×（2L＋1）＋（l＋L＋1）行元素bk，l為

列向量f中第（k＋K）×（2L＋1）＋（l＋L＋1）行元素為fk，l．式中σ為一次波估計p的方差．

在公式（9）中，矩陣A和列向量b都需要知道一次波估計p．而一次波估計p需要通過濾波器f得到．這里我們采用迭代最小二乘方法得到濾波器f的估計：

（1）設置最大迭代次數，并給定濾波器f的初始值；

（2）由式（8）得到當前濾波器下的一次波估計p，若一次波估計p的能量為零，則停止迭代；

（3）由p，f和s構造式（9）；

（4）求解式（9）得到更新后的濾波器f，然后轉入步驟（2）進行下一輪的迭代，直到濾波器f的系數基本不變或者迭代次數達到最大迭代次數．

在第（1）步中，通常濾波器的初始值可以設為隨機數，或者采用基于能量最小化的方法得到初始濾波器系數．第（2）步中為了減少迭代求解的次數，可以給定一個閾值，使得當濾波后的一次波能量的變化小于給定的閾值時，就停止迭代．方程（9）可以通過直接矩陣求逆求解f＝A－1b，也可采用Leivensen遞歸求解方程（ATA）f＝ATb，或者采用共軛梯度法進行求解，以減少由于矩陣求逆運算而帶來的計算量問題．在下面的實驗中，卷積信號盲分離方法的最大迭代次數設為5次，就可以得到滿意的結果，所以，和采用同樣濾波器參數的能量最小方法相比，所提出方法的計算量大概是前者的5倍．

4 應用實例

為了驗證所提方法的有效性，我們處理了人工合成數據和實際數據，并和基于能量最小化的多道匹配濾波器方法進行了比較．

首先處理了一個簡單的人工合成模型數據．模型數據是由一個一次波和一個多次波組成的共中心點道集，用來合成地震數據的子波在空間域上是變化的．時間采樣率是4ms．預測多次波由下式得到：

圖1a顯示了原始地震數據，圖1b為預測多次波，圖1c和圖1d分別是單道和多道濾波器的能量最小方法得到的一次波估計．圖1e和圖1f分別是單道和多道濾波器的卷積信號盲分離方法得到的一次波估計．在這個實驗中，單道濾波器參數為K＝2，多道濾波器參數為K＝2，L＝1．數據分塊參數為100采樣點和40道．可以看出，當預測和實際多次波之間存在空間上的差異時，單道的方法得不到滿意的結果．需要指出的是，單道濾波器可以利用更長的濾波器（即參數K取更大的值）來減少多次波殘余，但會造成一次波的畸變．多道濾波器的采用，可以在時間方向使用較少的濾波樣點，從而能更好地保護一次波．

為了定量評價基于卷積信號盲分離的多次波自適應相減方法的性能，采用如下定義的信噪比：

其中，p0是真實一次波，p是估計出的一次波．表1給出了簡單合成模型數據的信噪比結果，從圖1和表1可以看出，采用單道濾波器和多道濾波器時，本文所提出的方法的性能都優于基于能量最小的方法．

表1 簡單合成模型的信噪比Table 1 SNR of the simple model

圖1 一個簡單模型數據實驗Fig．1 A simple model data test

將卷積信號盲分離方法應用于Pluto數據，同時考慮到預測多次波和真實多次波之間可能具有空間上的差異，因此采用多道濾波器的設計方法．對整個數據采用分塊處理，以保證其滿足匹配濾波器設計的一致性假設條件．作為對比，多道能量最小化方法也用來處理同樣的數據．在本實驗中，兩種方法采用同樣的處理參數．每塊數據中包含30道，每道包含40個采樣點．多道濾波器參數為K＝7，L＝1．圖2是對Pluto數據分別用基于能量最小化的方法和本文方法得到的結果．圖2a是原始地震數據，圖2b是預測出的多次波．圖2c是基于能量最小化方法得到的一次波估計，圖2d是本文方法得到的一次波估計．圖3給出了圖2所給出結果的局部放大圖，圖3a和圖3b分別是圖2c和圖2d左側方框內數據的放大圖，圖3c和圖3d是圖2c和圖2d右側方框內數據的放大圖．從圖3中箭頭所指區域可以看出，本文方法的處理結果較基于能量最小化方法，能夠更好地壓制多次波，且估計出的一次波具有更高的連續性和信噪比．

我們將本文算法應用于一個海上數據集．在本實驗中，每塊數據中包含30道，每道包含40個采樣點．多道濾波器參數為K＝3，L＝1．圖4a是原始的共偏移距道集（偏移距為208m），圖4b是預測出的多次波，圖4c和圖4d分別是用某商業軟件的多道基于能量最小化方法和本文方法進行自適應相減所得的結果，黑色箭頭所指區域表明本文方法相對于能量最小化方法可更好地壓制多次波，白色箭頭所指區域表明本文方法可更好地保護一次波，壓制結果中一次波具有更好的連續性．

圖2 Pluto模型數據實驗Fig．2 Pluto data test

圖3 圖2的局部放大Fig．3 Areas zoomed

5 結 論

在基于匹配濾波器設計的多次波自適應相減問題中，原始地震信號可以表示為一次波和預測多次波的卷積混合，因此多次波自適應相減問題可以看作是從兩個觀測信號中分離出一個源信號的卷積盲分離問題．非高斯性最大化是盲信號分離中一個重要準則，認為獨立非高斯性信號在混合后的非高斯性要低于混合前信號的非高斯性．本文中，我們將盲信號分離理論中的非高斯性最大化技術和基于濾波器設計的多次波自適應相減技術相結合，提出一種基于多道卷積信號盲分離的多次波自適應相減算法．由于避免了采用一次波與多次波的正交性假設，因此與基于能量最小化的方法相比，所提出的方法能夠更有效地壓制多次波和保留一次波的能量，使單道的一次波估計結果具有更高的信噪比和連續性．對人工合成數據和實際數據的處理結果表明了算法的有效性．

一般預測多次波和真實多次波之間可能會存在空間上的差異，僅利用單道濾波器的設計方式很難校正空間差異，因此我們推薦使用多道濾波器進行多次波的自適應相減．多道濾波器參數K，L的取值會影響處理效果，對于某一待處理數據，我們首先利用幾個道集進行不同參數的實驗，然后選取效果好的參數來處理整個數據體．

非線性函數G的選取對衡量非高斯有重要的作用，本文中選取的是負指數函數．需要指出的是，非線性函數G也可以選取雙曲函數、對數雙曲函數等多種類型的函數．

附錄A

對于單道濾波器，濾波方程為

其中，pi為一次波估計，si為原始地震信號，為預測多次波．fk是單道濾波器，2 K＋1為濾波器在時間方向的長度．

由式（6）中的優化目標

對f求導，得

將上式展開，并求導：

由非線性函數式（7）得

另一方面，由濾波方程式（A1）得

由于地震數據的均值近似為零，因此我們通常只對一次波估計p進行方差歸一化后得到pˉ，σ表示一次波估計p的方差：

于是有

將式（A7）和式（A10）代入式（A6），于是鏈式求導為

零的結果變為

再將式（A1）代入上式，得到

將上式改寫為

整理上式就能得到如下方程組

其中

附錄B

附錄A中給出了單道濾波器的求解算法，多道濾波器的求解算法與此類似，只是迭代求解的方程組中的元素有所不同而已，附錄A中1D信號變為2D信號，如pi，分別變為pi，j和．

對于多道濾波器，濾波方程為

其中，pi，j為濾波后得到的一次波估計，si，j為原始地震數據，為預測多次波，i＝1，2，…，T，T為時間方向分塊的大小，j＝1，2，…，S，S為分塊的最大道數．fk，l為多道濾波器系數，2 K＋1為濾波器在時間方向的長度，2L＋1為濾波器在空間方向的道數．

由式（6）中的優化目標

對f求導，得

將上式展開，并求導：

由非線性函數式（7）得

另一方面，由濾波方程式（B1）得

由于地震數據的均值近似為零，因此我們通常只對一次波估計p進行方差歸一化后得到pˉ，即

于是有

將式（B7）和式（B10）代入式（B6），于是鏈式求導為零的結果變為

再將式（B1）代入上式，得到

將上式改寫為

整理上式可得如下方程組：

其中，矩陣A中的第（k＋K）×（2L＋1）＋（l＋L＋1）行，第（k1＋K）×（2L＋1）＋（l1＋L＋1）列元素

列向量b中第（k＋K）×（2L＋1）＋（l＋L＋1）行元素bk，l分別為

列向量f中第（k＋K）×（2L＋1）＋（l＋L＋1）行元素為fk，l．

（References）

［1］ Weglein A B．Multiple attenuation：an overview of recent advance and the road ahead．The Leading Edge，1999，18（1）：40－44．

［2］ Foster D J，Mosher C C．Suppression of multiple reflections using the Radon transform．Geophysics，1992，57（3）：386－395．

［3］ Lu W K，Zhang X G，Li Y D．Multiple removal based on detection and estimation of localized coherent signal．Geophysics，2003，68（2）：745－750．

［4］ Liu J，Lu W K．An improved predictive deconvolution based on maximization of non－Gaussianity．Applied Geophysics，2008，5（3）：189－196．

［5］ 薛亞茹，陳小宏，陸文凱．壓制多次波的正交多項式譜減法．地球物理學報，2009，52（3）：817－823．Xue Y R，Chen X H，Lu W K．Orthogonal polynomial spectrum subtraction for multiple attenuation．Chinese J．Geophys．（in Chinese），2009，52（3）：817－823．

［6］ Wiggins J W．Attenuation of complex water－bottom multiples by wave－equation－based prediction and subtraction．Geophysics，1988，53（12）：1527－1539．

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

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

［9］ Ma J T，Sen M K，Chen X H．Free－surface multiple attenuation using inverse data processing in the coupled planewave domain．Geophysics，2009，74（4）：V75－V81．

［10］ 龐廷華．基于匹配濾波器的多次波自適應相減［碩士論文］．北京：清華大學自動化系，2009．Pang T H．Adaptive multiple subtraction based on matching filters［Master′s thesis］（in Chinese）．Beijing：Automation Department of Tsinghua University，2009．

［11］ Manin M，Spitz S．3Dattenuation of targeted multiples with apattern recognition technique．∥57th Annual International Meeting，European Association of Geoscientists and Engineers，Extended Abstracts．1995，B046．

［12］ Guitton A，Cambois G．Prestack elimination of complex multiples：A Gulf of Mexico subsalt example．68th Annual International Meeting，SEG，Expanded Abstracts，1998，17：1329－1332．

［13］ Brown M，Clapp R G．T－x domain，pattern－based ground roll removal．∥70th Annual International Meeting，SEG，Expanded Abstracts．2000：2103－2106．

［14］ Guitton A．Coherent noise attenuation using inverse problems and prediction－error filters．First Break，2002，20（3）：161－167．

［15］ Fomel S．Applications of plane－wave destruction filters．Geophysics，2002，67（6）：1946－1960．

［16］ Guitton A．Multiple attenuation in complex geology with a pattern－based approach．Geophysics，2005，70（4）：V97－V107．

［17］ 陸文凱，駱毅，趙波等．基于獨立分量分析的多次波自適應相減技術．地球物理學報，2004，47（5）：886－891．Lu W K，Luo Y，Zhao B，et al．Adaptive multiple wave subtraction using independent component analysis．Chinese J．Geophys．（in Chinese），2004，47（5）：886－891．

［18］ Lu W K，Mao F．Adaptive multiple subtraction using independent component analysis．The Leading Edge，2005，24（3）：282－284．

［19］ Lu W K．Adaptive multiple subtraction using independent component analysis．Geophysics，2006，71（5）：S179－S184．

［20］ 李艷東．約束盲信號分離算法及應用研究［博士論文］．北京：清華大學自動化系，2006．Li Y D．Studies on algorithm and application of constrained blind source separation［Ph．D．thesis］（in Chinese）．Beijing：Automation Department of Tsinghua University，2006．

［21］ Lu W K，Liu L．Adaptive multiple subtraction based on constrained independent component analysis．Geophysics，2009，74（1）：V1－V7．

［22］ Luo Y，Kelamis P G，Wang Y．Simultaneous inversion of multiples and primaries：inversion versus subtraction．The LeadingEdge，2003，22（9）：814－818．

［23］ Wang Y H．Multiple subtraction using an expanded multichannel matching filter．Geophysics，2003，68（1）：346－354．

［24］ 李鵬，劉伊克，常旭等．均衡擬多道匹配濾波法在波動方程法壓制多次波中的應用．地球物理學報，2007，50（6）：1844－1853．Li P，Liu Y K，Chang X，et al．Application of the equipoise pseudomultichannel matching filter in multiple elimination using wave－equation method．Chinese J．Geophys．（in Chinese），2007，50（6）：1844－1853．

［25］ Guitton A，Verschuur D J．Adaptive subtraction of multiples using the L1－norm．Geophysical Prospecting，2004，52（1）：27－38．

［26］ Pang T H，Lu W K，Ma Y J．Adaptive multiple subtraction using a constrained L1－norm method with lateral continuity．Applied Geophysics，2009，6（3）：241－247．

［27］ Fomel S．Adaptive multiple subtraction using regularized nonstationary regression．Geophysics，2009，74（1）：25－33．

［28］ Walden A T．Non－Gaussian reflectivity，entropy，and deconvolution．Geophysics，1985，50（12）：2862－2888．

［29］ Hyvarinen A．Fast and robust fixed－point algorithms for independent component analysis．IEEE Transactions on Neural Networks，1999，10（3）：626－634．

［30］ Liu J．Echo cancellation using predictive deconvolution based on high order statistics．∥International Congress on Image and Signal Processing．2008：324－327．

Adaptive multiple subtraction based on multi－traces convolutional signal blind separation

LI Zhong－Xiao，LU Wen－Kai＊，PANG Ting－Hua，WANG Ji
State Key Laboratory of Intelligent Technology and Systems，Tsinghua National Laboratory for Information Science and Technology，Department of Automation，Tsinghua University，Beijing100084，China

This paper represents the adaptive multiple subtraction problem as a blind signal separation problem using multi－traces convolutional signal blind separation model．By expressing the difference between the predicted and true multiples using a 2Dconvolutional kernel，we propose an adaptive multiple subtraction method based on the multi－traces convolutional signal blind separation technique，which adopts maximization of the non－Gaussianity of the recovered primaries as the objective function．To solve the above non－linear optimization problem，we transfer it to an iterative linear one，which is realized by the iterative least squares algorithm．Taking advantage of the multi－traces convolutional signal blind separation model，the proposed method is applicable to the situation that there are differences in the time－space domain between the predicted and true multiples．Through the processing of the simple model data，the Pluto data and the real seismic data，the validity of the proposed method is demonstrated．

Adaptive multiple subtraction，Blind signal separation，Multi－traces convolutional signal，Iterative least squares

P631收修定稿2011－03－01，2011－12－26收修定稿

國家自然科學

（40874056）、國家科技重大專項（2011ZX05023－005－007）共同資助．

李鐘曉，男，1987年生，清華大學自動化系博士生，主要從事地震信號處理研究．E－mail：lizhongxiao09＠mails．tsinghua．edu．cn；thulzx＠gmail．com

＊通訊作者陸文凱，男，1969年生，清華大學自動化系研究員，主要從事地震信號處理研究．E－mail：lwkmf＠mail．tsinghua．edu．cn

李鐘曉，陸文凱，龐廷華等．基于多道卷積信號盲分離的多次波自適應相減方法．地球物理學報，2012，55（4）：1325－1334，

10．6038／j．issn．0001－5733．2012．04．028．

Li Z X，Lu W K，Pang T H，et al．Adaptive multiple subtraction based on multi－traces convolutional signal blind separation．Chinese J．Geophys．（in Chinese），2012，55（4）：1325－1334，doi：10．6038／j．issn．0001－5733．2012．04．028．

（本文編輯 胡素芳）