胡江濤，王華忠，王雄文

（同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院，上海200092）

射線束在地震道集或地震剖面上對(duì)應(yīng)一段線性同相軸，它在地震資料處理中被廣泛應(yīng)用于壓制采集過(guò)程中引入的隨機(jī)噪聲［1］，同時(shí)它還能用于壓制與有效反射波射線參數(shù)不同的線性噪聲，比如面波［2］、多次波［3］等。另外，射線束偏移［4－6］可以通過(guò)射線束形成篩選出有效的反射波，實(shí)現(xiàn)高效率和高信噪比的成像。

傳統(tǒng)的射線束形成方法［7］是局部τ－p變換，存在泄漏噪聲和低分辨率的問(wèn)題。在局部τ－p變換中引入相似系數(shù)濾波能夠壓制泄漏噪聲［8］，但不能提高射線束的分辨率；采用一般的最小二乘反演方法形成射線束能夠有效提高射線束的分辨率，但是不能壓制泄漏噪聲。這主要是由于上述方法中的約束條件沒有考慮波的方向性?？紤]方向信息的先驗(yàn)信息的加入能夠有效改善反演解的傾向性，是一種很有效的正則化方法。

我們將拉東譜能量作為先驗(yàn)信息加入到最小二乘反演過(guò)程中，改善射線束形成的分辨率低和泄漏噪聲的問(wèn)題。通過(guò)理論模型和實(shí)際資料測(cè)試證明了該方法的有效性，并將其應(yīng)用于面波的壓制。

1 最小二乘反演射線束形成方法

1.1 頻率空間域中射線束傳播正過(guò)程描述

為了得到高精度的射線束數(shù)據(jù)，一般采用反演方法形成射線束。而要進(jìn)行反演，首先需要研究相應(yīng)的射線傳播的正過(guò)程。射線束形成的正過(guò)程可以在時(shí)間空間域表示［9－10］，也可以在頻率空間域表示［11－12］。時(shí)間空間域表示正過(guò)程的優(yōu)點(diǎn)是反演時(shí)便于在時(shí)間和射線參數(shù)兩個(gè)方向上加入稀疏約束，缺點(diǎn)是效率較低。頻率空間域表示正過(guò)程的優(yōu)點(diǎn)是并行計(jì)算效率高，便于在射線參數(shù)方向上加入稀疏約束等。后者已得到廣泛應(yīng)用。

在頻率空間域，單個(gè)頻率的空間局部數(shù)據(jù)可以用數(shù)據(jù)矩陣表示：

射線束源（局部平面波源）用源矩陣表示：

射線束的時(shí)移矩陣為

式中：ω為圓頻率；p為射線參數(shù)；x為各道的空間位置；x0為參考道的位置，通常取中心道的位置。矩陣中每一列代表一個(gè)方向局部平面波傳播的相位移。

空間局部數(shù)據(jù)可以理解為射線束源經(jīng)過(guò)時(shí)移矩陣作用，在各局部空間點(diǎn)疊加產(chǎn)生的結(jié)果。因此，射線束形成的正過(guò)程可以表示為

類似地，我們可以寫出三維情形下的射線束形成的正過(guò)程，形成三維射線束。

1.2 頻率空間域中射線束形成的反問(wèn)題描述

從局部數(shù)據(jù)中估計(jì)射線束源可以表達(dá)為如下反演問(wèn)題：給定空間局部數(shù)據(jù)X和時(shí)移矩陣A，估計(jì)射線束源S。如果時(shí)移矩陣的逆矩陣存在，那么可以直接估計(jì)出射線束數(shù)據(jù)。由于時(shí)移矩陣通常不是滿秩的，上述反問(wèn)題通常是不定解。因此，一般通過(guò)最小化如下目標(biāo)函數(shù)來(lái)估計(jì)射線束源：

其中，λ為正則化因子，通常為較小的常數(shù)，它的主要作用是使反演數(shù)值計(jì)算過(guò)程穩(wěn)定。上述目標(biāo)泛函的最小二乘解為

將 Hessian矩陣AHA的逆矩陣項(xiàng)（AHA＋λI）－1忽略，得到近似解：

該近似解即為常規(guī)局部τ－p變換（線性拉東變換）的結(jié)果。由于方向矩陣A不是正交矩陣，所以AHA不等于單位陣，這就是局部τ－p變換結(jié)果存在泄漏噪聲和低分辨率問(wèn)題的原因。為了得到高分辨率的結(jié)果，Hessian矩陣的逆矩陣項(xiàng)不能忽略。但最小二乘解同樣存在泄漏噪聲的問(wèn)題，其結(jié)果可以通過(guò)加入先驗(yàn)信息來(lái)進(jìn)行改善。

2 拉東譜約束的射線束形成方法

先驗(yàn)信息可以通過(guò)正則化方式加入到反演過(guò)程中。Liu等［13］和 Wang等［14］通過(guò)將先驗(yàn)信息加入到正則化中來(lái)改善地震數(shù)據(jù)插值效果。Wang等［14］還分析了根據(jù)先驗(yàn)信息選擇正則化的原因。但是，不同情況下的先驗(yàn)信息是不同的。在射線束形成過(guò)程中，線性拉東變換結(jié)果Sa作為最小二乘反演的初始值，雖然受到泄漏噪聲和低分辨率的影響，但它的能量大小仍然能夠決定其作為有效射線束數(shù)據(jù)的可能性。能量越大，作為有效射線束的可能性就越大；能量越小，作為泄漏噪聲的可能性就越大。

用于約束的能量計(jì)算公式為

式中：ω1，ω2為有效頻帶的起始和終止頻率；i＝1，2，…，np。在實(shí)際計(jì)算時(shí)，需要將該能量歸一化到0～1。該能量約束因子對(duì)隨機(jī)噪聲不敏感，隨機(jī)噪聲對(duì)反演結(jié)果的影響小。通常，局部地震數(shù)據(jù)僅包含少數(shù)幾個(gè)射線參數(shù)，因此該能量也是一個(gè)比較稀疏的曲線。

拉東譜能量約束的目標(biāo)函數(shù)為

其中，R為np×np對(duì)角陣，其對(duì)角線上的元素為

其中，ε取一個(gè)很小的數(shù)，以使數(shù)值計(jì)算穩(wěn)定。拉東譜約束的目標(biāo)函數(shù)的最小二乘解為

通過(guò)共軛梯度方法求解（11）式，可以得到高質(zhì)量的射線束數(shù)據(jù)。該射線束可以作為射線束偏移的基礎(chǔ)數(shù)據(jù)，也可以用于壓制與有效信號(hào)具有不同射線參數(shù)的噪聲或隨機(jī)噪聲。在陸上地震資料中，比較常見的這類線性噪聲是面波。

3 數(shù)值試驗(yàn)

3.1 高分辨率的局部射線束產(chǎn)生

為了測(cè)試?yán)瓥|譜約束的射線束形成方法的有效性，合成一個(gè)具有兩個(gè)射線參數(shù)的局部數(shù)據(jù)，如圖1所示。圖2為局部τ－p變換得到的射線束數(shù)據(jù)。理論上，局部數(shù)據(jù)只包含兩個(gè)射線參數(shù)，理想的射線束合成結(jié)果應(yīng)該只包含兩個(gè)波形，但由圖2可見，局部τ－p變換結(jié)果除了兩個(gè)波形外還包含泄漏的噪聲，同時(shí)分辨率較低。圖3為常規(guī)最小二乘反演結(jié)果，正則化因子為0.01。對(duì)比圖2和圖3可以看出，采用反演方法得到的射線束分辨率有一定程度的提高，但仍然存在泄漏噪聲的影響。圖4為利用公式（8）計(jì)算的拉東譜能量曲線，可以看出，能量越大作為有效射線束的可能性越大；能量越小作為泄漏噪聲的可能性越大。因此該曲線可以用于約束最小二乘反演。圖5為理論局部數(shù)據(jù)拉東譜能量約束的最小二乘反演結(jié)果，λ＝0.01。對(duì)比圖2，圖3和圖5可以看出，拉東譜能量約束結(jié)果具有更高的分辨率，同時(shí)，泄漏的噪聲也得到了很好的壓制。因此，拉東譜約束的射線束形成方法能夠獲得更高質(zhì)量的射線束數(shù)據(jù)，為基于射線束數(shù)據(jù)的處理方法奠定了良好的基礎(chǔ)。

圖1 局部合成記錄

3.2 在線性噪聲壓制中的應(yīng)用

圖7 實(shí)際數(shù)據(jù)局部τ－p變換結(jié)果（a）與拉東譜能量約束的射線束形成結(jié)果（b）對(duì)比

采用某陸上實(shí)際地震數(shù)據(jù)對(duì)拉東譜約束的射線束形成方法壓制線性噪聲的效果進(jìn)行了測(cè)試（圖6）。該數(shù)據(jù)具有較強(qiáng)的隨機(jī)噪聲和面波，壓制面波的方法是切除面波所在的射線束數(shù)據(jù)。圖7為圖6數(shù)據(jù)中紅色框選定區(qū)域的射線束變換結(jié)果，其中圖7a為局部τ－p變換結(jié)果；圖7b為拉東譜能量約束的射線束變換結(jié)果。對(duì)比圖6和圖7可以看出，拉東譜約束的射線束形成方法在分辨率和泄漏噪聲的壓制方面明顯好于局部τ－p變換結(jié)果。圖7b中紅色箭頭所指區(qū)域?yàn)槊娌ㄕ共紖^(qū)域，將其切除后進(jìn)行反變換，即可得到壓制面波后的記錄。對(duì)圖6逐框進(jìn)行處理，得到壓制面波后的結(jié)果如圖8所示，可見實(shí)際數(shù)據(jù)中的面波得到了很好的壓制，被面波壓制的有效信號(hào)得到了很好的恢復(fù)。此外，隨機(jī)噪聲也得到了壓制，資料信噪比顯著提升。

圖8 面波壓制后的實(shí)際地震數(shù)據(jù)

4 結(jié)束語(yǔ)

高質(zhì)量射線束數(shù)據(jù)對(duì)于地震資料射線束偏移、線性同相軸提取、線性噪聲壓制、隨機(jī)噪聲壓制等基于射線束的處理方法至關(guān)重要。研究表明：與常規(guī)方法相比，拉東譜約束的最小二乘反演射線束形成方法能提高射線束分辨率和壓制泄漏噪聲；同時(shí)，拉東譜能量對(duì)隨機(jī)噪聲不敏感，隨機(jī)噪聲對(duì)射線束形成的影響較小。因此，該方法能提取出高質(zhì)量的射線束數(shù)據(jù)，有效壓制面波類線性噪聲。

［1］Ozbek A.Adaptive beamforming with generalized linear constraints［J］.Expanded Abstracts of 70thAnnual Internat SEG Mtg，2000，2081－2084

［2］Panea I，Drijkoningen G.The spatial data－adaptive minimum－variance distortionless－response beamformer on seismic single－sensor data［J］.Geophysics，2008，73（5）：Q29－Q42

［3］胡天躍，王潤(rùn)秋，White R E.地震資料處理中的聚束濾波方法［J］.地球物理學(xué)報(bào)，2000，43（1）：105－115 Hu T Y，Wang R Q，White R E.Beamforming in seismic data processing［J］.Chinese Journal of Geophysics，2000，43（1）：105－115

［4］Hill R N.Prestack Gaussian－beam depth migration［J］.Geophysics，2001，66（4）：1240－1250

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

［6］Gao F，Zhang P，Wang B，et al.Fast beam migrationa step toward interactive imaging［J］.Expanded Abstracts of 76thAnnual Internat SEG Mtg，2006，2470－2473

［7］Capon J.High－resolution frequency－wavenumber spectrum analysis［J］.Proceeding of the IEEE，1969，57（8）：1408－1418

［8］Stoffa P L，Buhl P，Diebold J B，et al.Direct mapping of seismic data to the domain of intercept time and ray parameter－aplane wave decomposition［J］.Geophysics，1981，46（3）：255－267

［9］Thorson J R，Claerbout J F.Velocity－stack and slant－stack stochastic inversion［J］.Geophysics，1985，50（12）：2727－2741

［10］Cary P W.The simplest discrete Radon transform［J］.Expanded Abstracts of 68thAnnual Internat SEG Mtg，1998，1999－2002

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

［12］Hu J T，Wang X W，Wang H Z.Beam forming by least square inversion with Radon spectrum constraints［J］.Expanded Abstracts of 82ndAnnual Internat SEG Mtg，2012，1－5

［13］Liu B，Sacchi M D.Minimum weighted norm interpolation of seismic records［J］.Geophysics，2004，69（6）：1560－1568

［14］Wang X W，Wang H Z.Minimizing ellipsoidal norm seismic data interpolation with Radon spectrum constraints［J］.Expanded Abstracts of 73rdEAGE Annual Conference，2011，A046