孫維薔 王華忠 符瓊琲

（同濟大學海洋與地球科學學院波現(xiàn)象與反演成像研究組，上海200092）

空間方向譜估計CLEAN算法及其在地震數據處理中的應用

孫維薔*王華忠 符瓊琲

（同濟大學海洋與地球科學學院波現(xiàn)象與反演成像研究組，上海200092）

傳統(tǒng)的射線束形成方法存在泄漏噪聲和分辨率低的問題。為此，引入CLEAN算法，以空間方向譜作為反演解的評價準則，實現(xiàn)高分辨率的射線束形成。利用聚類和定位方法，進一步提高了常規(guī)CLEAN算法的空間分辨率，使該方法能更好地應用于實際地震資料去噪。理論數據和實際資料的應用結果表明，該方法可以有效壓制泄漏噪聲，同時提高射線束的分辨率。與傳統(tǒng)的局部傾斜疊加和最小二乘意義下的射線束形成方法相比，文中方法具有明顯的理論優(yōu)勢和較好的應用價值。

空間方向譜 協(xié)方差矩陣 CLEAN算法 聚類和定位 面波壓制

1 引言

射線束在很多學科領域中都有重要應用，譬如雷達、射電天文、地震勘探等。應用射線束方法能夠壓制與有效反射波具有不同射線參數的噪聲，如隨機噪聲及面波和多次波等線性噪聲［1-3］。除了地震資料去噪，射線束的另一個重要應用是射線束偏移成像，通過射線束篩選出有效的局部方向反射波場，實現(xiàn)高效和高信噪比成像［4，5］。

射線束形成方法本質上是一種線性信號空間方向譜估計方法［6］。局部傾斜疊加是最常見的射線束形成方法［7］，它是線性 Radon變換（τ-P變換）的一種具體實現(xiàn)方式。由于τ-P變換的基函數非正交、非完備，導致其結果中包含大量的泄漏噪聲，嚴重降低了τ-P變換的分辨率。提高射線束分辨率的基本思想是建立參數化模型預測局部線性信號，在不同類型的范數下反演估計射線束，這就將射線束形成問題轉化為參數估計問題。該類方法通常用實際數據與預測數據的差作為誤差范函，在最小二乘意義下估計射線束［8-10］。最小二乘意義下的估計結果在信噪比方面有比較顯著的改善，但其分辨率依然受測不準原理限制。為了提高分辨率，Sacchi等［8］提出用柯西范數代替L2范數，利用有效信號的稀疏性作為約束條件，實現(xiàn)了頻率域高分辨率的線性Radon變換；Cary［9］指出時間域的線性Radon變換結果具有更強的稀疏性，更加適合描述地震數據中的時變信號，可以更準確地估計線性信號。對于局部線性信號而言，在高斯白噪聲的假設下，其對應的自相關矩陣的特征向量就是使預測誤差最小的線性信號。在非高斯白噪聲情況下，構造類似L1／L0范數意義下的估計可以得到更稀疏的線性信號預測結果。Hu等［10］用L1／L0范數下的稀疏求解得到了更高分辨率的射線束?？傊?，高分辨率的射線束形成方法可以有效地估計地震數據中的線性信號，更好地壓制數據中的隨機噪聲，或實現(xiàn)更好的插值等，是一項實用價值較高的處理技術。

基于自相關矩陣分解及合成，建立誤差泛函，再求解最佳逼近自相關矩陣的方向譜，是一種很有啟發(fā)性的方向譜估計方法。CLEAN算法就是基于這一思想發(fā)展起來的半參數化高分辨率空間方向譜估計方法［11，12］。近年來，一類低秩近似的插值方法也是基于自相關矩陣的特征而設計的［13］。本文將CLEAN算法用于地震數據處理中，提出了一種數據驅動類射線束形成方法。在傳統(tǒng)CLEAN算法估計空間方向譜的基礎上，從協(xié)方差矩陣模型出發(fā)，結合射線束形成的正問題、協(xié)方差矩陣的正交分解和協(xié)方差矩陣擬合的反問題，并利用空間方向譜作為反演解的評價準則，估計真實的射線束。類似匹配追蹤的求解方法，它能在壓制泄漏噪聲的同時，提高射線束估計結果的空間分辨率。引入聚類和定位算法優(yōu)化了常規(guī)CLEAN算法所形成的射線束結果，可進一步提高射線束的空間分辨率。

2 理論與方法

射線束形成可以抽象地理解為：局部線性信號預測器的建立和基于該預測器的不同誤差范數意義下的最優(yōu)估計兩個關鍵步驟?；谠撍枷?，可以發(fā)展出各種不同的射線束形成方法。

2.1 L2誤差范數下射線束形成方法

高精度的射線束，可以根據射線束形成的正問題，求解反問題的稀疏解得到。對于不同的頻率，射線束形成正問題可以描述為

式中：ω為圓頻率；向量X（ω）和S（ω）分別為單頻空間局部數據向量和射線束向量；時移矩陣A（ω）中每一列向量表示射線束形成的基函數，可以表示為

式中：m為空間局部數據的道數；n為射線參數的個數；Pi為射線參數；向量a（Pi）為射線束形成的第i個基函數；x i為局部空間數據的位置；x0為參考道的位置，通常取中心道的位置。

以上僅為二維頻率—空間域正問題表達形式，更高維空間數據的排列參考螺旋坐標排列方式［14］。射線束形成方法可以從如下反演框架來解釋。反問題可以描述為已知局部空間數據和時移矩陣求解射線束的過程。通過最小化目標泛函估計射線束

式中λ為權系數。該泛函的解可表示為

式中：“H”表示共軛轉置；I表示單位矩陣。將Hessian逆矩陣項（AHA+λI）-1忽略，可得到近似解

該近似解就是局部傾斜疊加形成的射線束結果。由于時移矩陣A不是正交矩陣，所以AHA不等于單位矩陣，說明局部傾斜疊加形成的射線束結果中存在泄漏噪聲和低分辨率問題。為了得到高分辨率的結果，Hessian逆矩陣項不能被忽略。但在數據孔徑小的情況下最小二乘估計的射線束同樣分辨率不高，需要加入先驗信息約束。分辨率更高的解需要在射線束稀疏性的假設下用L1／L0范數求解。

2.2 CLEAN算法高分辨率射線束形成方法

CLEAN算法是一種半參數化的空間方向譜估計方法，它基于協(xié)方差矩陣模型擬合觀測樣本的協(xié)方差矩陣，實現(xiàn)對真實射線束的估計。求解過程類似于匹配追蹤，因此它能提高射線束的空間分辨率。根據射線束形成的正問題，局部空間數據的協(xié)方差矩陣可以表示為

式中：R表示協(xié)方差矩陣；E｛·｝表示期望；“T”表示轉置。CLEAN算法假設具有不同射線參數的射線束彼此不相關，則射線束的協(xié)方差矩陣退化為對角陣，即，其中是射線參數為P i的射線束的功率。因此，式（6）所示的協(xié)方差矩陣R可以分解為一系列與方向角度相關的子協(xié)方差矩陣

由于不能從有限的觀測數據樣本中精確地確定理論協(xié)方差矩陣R，因此需要由協(xié)方差矩陣的估值代替R，即式（7）為CLEAN算法所建立的協(xié)方差矩陣半參數化模型，是算法的基礎，其中的每一項對應于射線參數為Pi的射線束的協(xié)方差矩陣?？梢钥闯?，所謂半參數化，主要是引入了基于局部數據估計得到的協(xié)方差矩陣。CLEAN算法根據上述協(xié)方差矩陣模型，對觀測記錄的協(xié)方差矩陣進行擬合的同時，估計參數，從而實現(xiàn)真實射線束的估計。在最小二乘意義下建立目標泛函擬合協(xié)方差矩陣

CLEAN算法通過迭代反演來逐步實現(xiàn)參數估計。第一步迭代假設簡化式（8），建立估計參數的目標泛函

上式等價于

其解可以表示為

與第一次迭代類似，式（12）的解可以表示為

其中第二次迭代的空間方向譜為

由式（13）和式（14）可知，第二次迭代是利用上次迭代所估計的參數更新協(xié)方差矩陣，由新的協(xié)方差矩陣計算空間方向譜的更新量從上一步的空間方向譜中減去該更新量的一部分，實現(xiàn)對空間方向譜的更新。ρ為步長，0＜ρ＜1，通常取值為0.1～0.25。ρ＜1的原因是：在空間方向譜的峰值處附近可能還存在若干個源信號，若只減去對應于參數為的空間方向譜的一部分，則在1處及其附近將剩余一些功率。在后續(xù)的迭代過程中，CLEAN算法還可以反演位于峰值1處附近的其他射線束。若將空間方向譜全部減去，可能一些射線束的譜會出現(xiàn)負值，這不符合物理現(xiàn)實。綜上所述，重復以上迭代過程，直至迭代次數達到預設值或某個射線束的空間方向譜出現(xiàn)負值為止。

圖1a數據由兩個常射線參數的同相軸組成，紅框圈定區(qū)域是作為輸入的空間局部數據。根據采樣定理，波數方向的分辨率為，其中Δx為空間采樣間隔。由于波數和射線參數之間存在定量關系：k=ωP，射線束分辨率可以表示為ΔP=因此，空間局部數據窗口越小，分辨率越低。由于地震數據中的有效信號僅在局部范圍內才滿足線性假設，這就決定了空間局部數據的窗口不宜太大，因此局部傾斜疊加只能得到較低分辨率的射線束。圖1b和圖1c為局部傾斜疊加和CLEAN算法的射線束形成結果。常規(guī)局部傾斜疊加所形成的射線束中含有很強的泄漏噪聲，而且分辨率較低。CLEAN算法所形成的射線束壓制了大部分的泄漏噪聲，同時也提高了空間分辨率。圖2展示了迭代過程中的空間方向譜的變化情況。原始空間方向譜的兩個峰值都集中在真實射線束成分的射線參數附近。說明空間方向譜可以作為評價真實射線束成分的有效標準。每次迭代會反演出一個真實的射線束成分，同時用來更新上次迭代的空間方向譜，直至達到終止條件。

圖1 包含兩個常射線參數同相軸的理論數據（a）及局部傾斜疊加（b）和CLEAN算法（c）的射線束形成結果對比

圖2 迭代過程中空間方向譜的變化

為了進一步優(yōu)化CLEAN算法得到的射線束的結果，提高空間分辨率，引入聚類和定位算法。在聚類算法中，根據射線束分辨率ΔP，定義了最大可分離距離εΔP，其中ε為經驗參數，通常ε=0.5。聚類算法的作用是對常規(guī)CLEAN算法所確定的射線束進行分類，所有距離小于最大可分離距離的射線束歸為同一個射線束。當空間局部窗內波現(xiàn)象較復雜時，常規(guī)CLEAN算法所確定的射線束存在假象，進而降低了空間分辨率。據此，聚類算法可以用于優(yōu)化常規(guī)CLEAN算法，進一步提高射線束的空間分辨率。

定位算法的作用是確定每一類射線束中的真實射線束的位置。利用加權均方差建立泛函來估計每一類射線束中真實射線束成分的位置，即

式中：l ij表示第i類中的第j個射線束的位置，hij是其振幅值；ˉh i是第i類中所有射線束的平均振幅值；Ni為第i類射線束的射線個數。求解泛函式（15），得到第i類射線束中的真實射線束位置的估計值

3 數值試驗

3.1 理論數據測試

圖3a為有限差分正演所得的單炮記錄，其中紅框區(qū)域為用于形成射線束的空間局部數據。圖3b～圖3d比較了局部傾斜疊加、最小二乘反演和高分辨率CLEAN算法所形成的射線束。由于局部傾斜疊加的基函數非完備正交，導致結果中包含大量的泄漏噪聲，這些噪聲會嚴重降低射線束的空間分辨率；最小二乘意義下的射線束估計結果在信噪比方面有比較顯著的改善，但其分辨率依然受測不準原理限制，空間積分孔徑越小，分辨率越低；高分辨率的CLEAN反演方法不僅壓制了泄漏噪聲，同時也提高了空間分辨率。圖4描述了迭代過程中空間方向譜的變化情況。原始空間方向譜的峰值都集中在對應于真實射線束成分的射線參數附近，說明空間方向譜可以作為評價真實射線束的有效標準。每次迭代會反演出一個真實的射線束，同時更新了上次迭代的空間方向譜，直至達到迭代終止條件為止?？臻g方向譜在迭代過程中的變化證明了新方法可以通過逐步迭代，有效地反演出真實射線束。

圖3 理論記錄不同方法的結果對比

圖4 迭代過程中空間方向譜的變化

圖5 實際單炮記錄

3.2 實際數據測試

將基于CLEAN算法的高分辨率射線束形成方法用于實際數據壓制面波噪聲。面波壓制過程中，首先利用基于CLEAN算法的高分辨率射線束形成方法將時空域中的觀測記錄變換到局部平面波域；在平面波域中根據有效信號與面波的可分離性，對面波進行壓制；最后將平面波域壓制面波后的記錄反變換回時空域。圖5為野外觀測的單炮記錄，其中面波能量很強。圖6對比了圖5紅框區(qū)域中數據的局部傾斜疊加和基于CLEAN算法所形成的射線束。由于基函數非完備正交且觀測孔徑有限，局部傾斜疊加形成的射線束中存在大量泄漏噪聲且空間分辨率低，模糊了面波在平面波域的特征，難以實現(xiàn)面波與有效波的分離（圖6a）。在CLEAN算法產生的射線束中，由于泄漏的噪聲得到了壓制且空間分辨率較高，使得面波和有效波可以明顯區(qū)分（圖6b，紅色箭頭所指即為面波區(qū)域）。在基于CLEAN算法的高分辨率射線束形成方法的基礎上，可以通過切除壓制面波。圖7為壓制面波后的炮記錄及壓制掉的面波，可以發(fā)現(xiàn)，大部分面波已經被壓制，證明了本文方法的有效性。

4 結論與討論

基于CLEAN算法的高分辨率射線束形成方法從協(xié)方差矩陣的半參數化模型出發(fā)，以射線束形成的方向譜的分辨率作為反演解的評價標準，實現(xiàn)真實射線束的估計，提高了射線束的空間分辨率。數值試驗結果證明了本方法在提高射線束空間分辨率方面的有效性和可行性，未考慮提高射線束時間分辨率是本文方法的不足之處。將高分辨率射線束形成方法用于實際資料的面波壓制中，取得了較好的效果。除了面波壓制，高分辨率的射線束形成方法還可以應用于線性與隨機噪聲壓制、地震數據規(guī)則化、地震數據壓縮和射線束偏移等地震數據處理中，具有重要的應用價值。

圖6 實際數據不同方法形成射線束對比

圖7 實際數據面波壓制后的單炮記錄（a）及壓制的面波（b）

［1］ 張軍華，呂寧，田連玉等.地震資料去噪方法綜合評述.石油地球物理勘探，2005，40（S1）：121-127.Zhang Junhua，LüNing，Tian Lianyu et al.Overview of seismic data de-noising methods.OGP，2005，40（S1）：121-127.

［2］ 朱生旺，魏修成，李鋒等.用拋物線Radon變換稀疏解分離和壓制多次波.石油地球物理勘探，2002，37（2）：110-116.Zhu Shengwang，Wei Xiucheng，Li Feng et al.Separating and suppressing multipies by sparse solution of parabolic Radon transform.OGP，2002，37（2）：110-116.

［3］ 吳艷輝，唐博文，柯本喜等.GeoEast多次波壓制技術在陸上地震資料處理中的應用.石油地球物理勘探，2014，49（S1）：93-98.Wu Yanhui，Tang Bowen，Ke Benxi et al.Multiple suppression technique based on GeoEast and its application on onshore seismic data processing.OGP，2014，49（S1）：93-98.

［4］ Gray S H.Gaussian beam migration of common-shot records.Geophysics，2005，70（4）：S71-S77.

［5］ 王華忠，馮波，劉少勇等.疊前地震數據特征波場分解、偏移成像與層析反演.地球物理學報，2015，58（6）：2024-2034.Wang Huazhong，F(xiàn)eng Bo，Liu Shaoyong et al.Characteristic wavefield decomposition，imaging and inversion with prestack seismic data.Chinese Journal of Geophysics，2015，58（6）：2024-2034.

［6］ Stoica P，Mosese R L.Spectral Analysis of Signals.Prentice-Hall，Upper Saddle River，New Jersey，2005.

［7］ Neidell N S，Taner M T.Semblance and other coherency measures for multichannel data.Geophysics，1971，36（3）：483-497.

［8］ Sacchi M，Ulrych T.High-resolution velocity gathers and offset space reconstruction.Geophysics，1995，60（4）：1169-1177.

［9］ Cary P.The simplest discrete Radon transform.SEG Technical Program Expanded Abstracts，1998，17：1999-2002.

［10］ Hu J，Wang X，Wang H.Beam forming by least square inversion with radon spectrum constraints.SEG Technical Program Expanded Abstracts，2012，31：PAPR1140.

［11］ H?gbom J A.Aperture synthesis with a non-regular distribution of interferometer baselines.Astronomy＆Astrophysics Supplement，1974，15（15）：417-426.

［12］ Cornwell T and Bridle A.Deconvolution Tutorial（Technical report）.National Radio Astronomy Observatory，1996.

［13］ Yang Y，Ma J，Osher S.Seismic data reconstruction via matrix completion.Inverse Problems and Imaging，2013，4（4）：1379-1392.

［14］ Rickett J E and Claerbout J F.Calculation of the Sun＇s acoustic impulse response by multi-dimensional spectral factorization.Solar Physics，2000，192（1）：203-210.

P631

A

10.13810／j.cnki.issn.1000-7210.2017.01.004

孫維薔，王華忠，符瓊琲.空間方向譜估計CLEAN算法及其在地震數據處理中的應用.石油地球物理勘探，2017，52（1）：20-26.

1000-7210（2017）01-0020-07

*上海市同濟大學海洋與地球科學學院波現(xiàn)象與反演成像研究組，200092。Email：tswq_123456＠m(xù)sn.com

本文于2016年1月28日收到，最終修改稿于同年12月3日收到。

本項研究受國家重點基礎研究“973”發(fā)展計劃項目（2011CB201002）和國家科技重大專項（2011ZX05005-005-008HZ）聯(lián)合資助。

（本文編輯：宜明理）

孫維薔 博士，1988年生；2010年獲同濟大學海洋與地球科學學院地球物理學專業(yè)學士學位；2016年獲該校海洋與地球科學學院地球物理學專業(yè)博士學位，主要從事多次波壓制與成像方法研究。