基于高階Radon域的同時震源波場分離

朱翔宇

1.中國科學(xué)院測量與地球物理研究所計算與勘探地球物理研究中心；大地測量與地球動力學(xué)國家重點(diǎn)實(shí)驗室，湖北 武漢 430077；2.中國科學(xué)院大學(xué)，北京 100049

0 引言

常規(guī)地震數(shù)據(jù)采集，不同震源之間需要足夠大的時間間隔以避免相互之間的干擾[1]，而同時震源技術(shù)突破了常規(guī)采集技術(shù)的瓶頸。在震源數(shù)固定時，同時震源技術(shù)可以降低采樣時間，提高采集效率；而在相同的采樣時間內(nèi)，通過增加震源密度，可以提升成像質(zhì)量[2]。并且該技術(shù)是寬方位高密度地震數(shù)據(jù)采集的有效方法，為深層地球物理勘探提供基礎(chǔ)。20世紀(jì)末，Beasley[3]首次提出同時震源（Simultaneous Sources）這一概念，2008年Berkhout[4]詳細(xì)介紹了同時震源在采集效率及成像質(zhì)量上的優(yōu)勢，使得該方法迅速的發(fā)展起來。

同時震源數(shù)據(jù)處理一般可以分為直接偏移[5]和先進(jìn)行分離然后偏移兩種方法。分離類方法將混疊在一起的炮記錄進(jìn)行分離，得到單炮記錄，隨后可以利用常規(guī)的方法進(jìn)行處理。近幾年研究表明分離的方法可以分為兩大類即主動分離[6]和被動分離[7]。被動分離[8]一般是從偽分離開始，利用相干濾波進(jìn)行處理。主動分離[9]是將炮分離看成一個反演問題，利用反演的方法進(jìn)行處理，因此主動分離方法也可以稱之為反演類方法，是將同時震源數(shù)據(jù)轉(zhuǎn)換到稀疏域，在稀疏域內(nèi)利用反演進(jìn)行分離重構(gòu)。該方法的關(guān)鍵就在于稀疏域的選取。

處理同時震源數(shù)據(jù)常用的稀疏域有f-k域[10]，Seislet域[11]，Curvelet域[12]，以及 Radon 域[9][13]。由于共檢波點(diǎn)域中主炮的數(shù)據(jù)是相干的，其他炮的數(shù)據(jù)相對于主炮是不相干的，在Radon域中這兩者分別變現(xiàn)為集中的能量點(diǎn)與分散的噪聲，借助這個性質(zhì)可以達(dá)到分離的效果，并且Radon變換本身就有去噪能力[14]。

本文利用一種迭代估計的方法進(jìn)行炮分離，我們首先將共檢波點(diǎn)道集數(shù)據(jù)變換到Radon域，然后利用閾值算子[15]約束模型，通過迭代不斷地修正模型，提升分離的信噪比。我們使用的稀疏域為Radon域，但傳統(tǒng)的Radon變換正反變換通常會對原始地震數(shù)據(jù)產(chǎn)生損傷，當(dāng)?shù)卣饠?shù)據(jù)含有較強(qiáng)AVO特征時，地震數(shù)據(jù)經(jīng)Radon變換后在振幅上會出現(xiàn)較大誤差；高階Radon變換引入高階多項式，增強(qiáng)對同相軸保護(hù)能力[16]，改善分離結(jié)果。整形正則化可以加快模型的收斂速度，提高反演效率[11]。通過簡單模型和復(fù)雜模型數(shù)據(jù)的處理可以驗證本文方法的正確性以及實(shí)用性。

1 基于整形正則化的反演分離

海上多源地震可以分為兩種激發(fā)方式：高密度震源覆蓋以及寬范圍震源覆蓋[17]。其中高密度覆蓋時，主震源船正常激發(fā)，輔震源船延遲激發(fā)，然后主震源船與輔震源船同向行駛，在下一個位置，仍然保持主震源船正常激發(fā)，輔震源船延遲激發(fā)，依此類推得到混合數(shù)據(jù)。在這種激發(fā)模式下，混合數(shù)據(jù)在共檢波點(diǎn)道集表現(xiàn)出以下特征：主震源反射信號是連續(xù)、相干的，而輔震源的反射信號由于時間延遲的存在而不能形成連續(xù)、相干的同相軸[18]。本文的分離方法利用上述特征，將共檢波點(diǎn)域混合數(shù)據(jù)進(jìn)行分離，雖然該方法可以分離多炮混合的數(shù)據(jù)，但為方便描述以雙源混合為例，介紹數(shù)據(jù)混合的過程，混合數(shù)據(jù)的表達(dá)為：

其中d為混合數(shù)據(jù)，d1和d2是需要分離的單炮數(shù)據(jù)，而T代表時間延遲算子。在這種混合方式下，震源1的數(shù)據(jù)d1是連續(xù)相干的，而震源2的數(shù)據(jù)d2作用了一個時間延遲算子T，變?yōu)椴贿B續(xù)、不相干的信號。其中時間延遲算子T的表達(dá)式為：

方程（2）中F與F-1分別代表傅里葉正反變換，P代表時間延遲的編碼方式，其表達(dá)式為：

方程（3）中δtn代表時間延遲。對方程（1）左右兩邊同時乘以時間算子的逆T-1，可以得到：

結(jié)合方程（1）、（4）可以得到一個增廣方程：

其中

上式中I為單位陣。方程（5）即我們需要求解的方程。

方程（5）可以利用基于整形正則化的迭代方法求解[11]（Chenet al，2014），迭代公式可以表示為：

方程（7）中，S是整形算子，與選取的稀疏域有關(guān)，B是逆映射算子，B的近似值為二分之一單位陣，即：

方程（5）的解可以表示成L1約束下目標(biāo)函數(shù)的解：

方程（9）中A-1代表變換算子，μ為正則化參數(shù)。方程（7）中整形則正則化參數(shù)s的表達(dá)式為：

結(jié)合我們使用的方法，方程（10）中的A、A-1分別代表Radon變換及其逆變換。共檢波點(diǎn)道集中，主炮的信號在Radon域是聚焦的能量點(diǎn)，而輔炮則是分散的噪聲，Radon變換可以為數(shù)據(jù)提供良好稀疏性，提高反演的精度。離散Radon變換算子L的表達(dá)式為：

方程（11）中ω是頻率，qi為慢度，χm為道間距，當(dāng)n=1時代表線性Radon變換，n=2為拋物Radon變換。頻率域離散Radon可以表示為：

D是頻率域地震數(shù)據(jù)，U是Radon域地震數(shù)據(jù)，L是變換算子。由于L一般不是方陣，傳統(tǒng)的方式是通過最小二乘求其廣義逆即：

同樣給出高階拋物Radon變換的反變換在頻率域的表達(dá)式：

其中pj(x)表示第j階正交多項式基函數(shù)，上述方程改寫成矩陣形式：

式中，D(f)為頻率域地震數(shù)據(jù)；LHPRT(f)為高階拋物Radon變換算子矩陣；變換算子矩陣LHPRT(f)的構(gòu)造形式如（3-20）式所示：

L0,L1,L2分別為疊加算子，梯度算子與曲率算子，與常規(guī)Radon相比，高階變換的變換算子LHPRT(f)將Radon變化進(jìn)行延拓，補(bǔ)充振幅的高階多項式特征，有利于保護(hù)同相軸的真振幅[23]。

利用最小二乘方法得到高階拋物Radon變換的正變換計算公式為：

方程（10）中另一個參數(shù)Tγ指的是帶有輸入?yún)?shù)γ的閾值算子，閾值算子是將稀疏域中較大的值保留，較小的值舍去。常用到的閾值算子包含兩類：硬閾值和軟閾值。硬閾值是將大于閾值的系數(shù)保留，小于閾值的系數(shù)舍去；軟閾值是將大于閾值的系數(shù)做收縮處理后保留，小于閾值的參數(shù)直接舍去[21]。硬閾值的表達(dá)式為：

軟閾值的表達(dá)式為：

結(jié)合后續(xù)的算例本文采用硬閾值進(jìn)行處理。常用的閾值參數(shù)γ的選擇方式有常數(shù)、線性衰減、指數(shù)衰減、百分位等方式，不同方式的閾值參數(shù)選取會影響迭代的收斂速度[22]。本文采用百分位法。

2 算例分析

2.1 拋物算例

用簡單的拋物算例來測試分離方法的有效性，拋物算例包含四條拋物型同相軸，四條同相軸的能量隨深度依此減弱，其中圖1中（a）為未混合數(shù)據(jù)，對應(yīng)公式（1）中的d1，圖1中（b）為混合數(shù)據(jù)，公式（1）對應(yīng)于d。拋物算例中未混合信號由于沒有時間延遲因此是是連續(xù)相干的，而混合信號中輔炮的信號由于存在時間延遲變的不連續(xù)、不相干。圖2、圖3均是是利用方程（7）的迭代方法的分離結(jié)果，其中圖2基于傳統(tǒng)的Radon變換，圖3基于高階Radon變換。兩種分離方法應(yīng)用的迭代方程相同，只是稀疏域選取的不同，從分離結(jié)果看，兩種方法均將主炮信號與輔炮信號進(jìn)行了分離，但是從殘差可以看出，高階Radon變換對有效信號具有更好的保幅效果，為了評價分離效果，本文利用信噪比進(jìn)行評估，信噪比的定義為[23]∶

方程（17）中d代表未混合數(shù)據(jù)，dn代表第n次迭代的結(jié)果。

圖1 拋物算例Fig.1 Parabolic example

圖2 拋物算例傳統(tǒng)Radon變換分離結(jié)果Fig.2 Separation result of traditional Radon transform for parabolic example

圖3 拋物算例高階Radon變換分離結(jié)果Fig.3 Separation result of high order Radon transform for parabolic example

利用方程（17）的信噪比公式分別評估兩種方法的分離效果，結(jié)果見圖4，其中藍(lán)色實(shí)線為高階Radon變換的分離效果，橙色虛線為傳統(tǒng)Radon變換的分離效果，兩種方法在迭代十五次時基本收斂，且高階Radon變換的分離結(jié)果收斂于25.92 dB，而傳統(tǒng)Radon變換的分離結(jié)果收斂于20.77 dB。從信噪比來看，高階Radon變換的分離效果要優(yōu)于傳統(tǒng)Radon變換的分離效果。

2.2 SEAM_2B模型算例

第二個算例，我們用SEAM_2B模型的雙源混合數(shù)據(jù)來模擬實(shí)際地震記錄。圖5是SEAM_2B模型共檢波點(diǎn)道集地震記錄，其中圖5中（a）為未混合數(shù)據(jù)，對應(yīng)公式（1）中的d1，圖5中（b）為混合數(shù)據(jù)，對公式（1）應(yīng)于d。

圖4 拋物算例兩種分離方法信噪比曲線Fig.4 S/N ratio for parabolic example

圖5、圖6分別是兩種方法利用方程（7）的迭代方法的分離結(jié)果，其中圖5基于傳統(tǒng)的Radon變換，圖6基于高階Radon變換。利用方程（17）的信噪比公式評估分別評價分離效果，結(jié)果見圖7，SEAM_2B模型算例的信噪比也在十五次迭代后基本趨于穩(wěn)定，高階Radon變換的分離結(jié)果收斂于24dB，而傳統(tǒng)Radon變換的分離結(jié)果收斂于19.33dB。從信噪比來看，高階Radon變換的分離效果要優(yōu)于傳統(tǒng)Radon變換的分離效果，并且，高階Radon變換的分離結(jié)果對有效信號有很強(qiáng)的保護(hù)效果。

圖5 SEAM_2B模型算例Fig.5 SEAM_2B example

圖6 SEAM_2B模型傳統(tǒng)Radon變換分離結(jié)果Fig.6 Separation result of traditional Radon transform for SEAM_2B example

圖7 SEAM_2B模型高階Radon變換分離結(jié)果Fig.7 Separation result of high order Radon transform for SEAM_2B example

3 結(jié)論

本文通過高階Radon變換將共檢波點(diǎn)域數(shù)據(jù)轉(zhuǎn)換到稀疏域，利用閾值算子約束模型，通過迭代不斷地修正模型，提升分離的信噪比。通過兩個模型的試算，可以證明高階Radon選取的合理性以及迭代分離方法的準(zhǔn)確性。并且我們的方法只需要較小的迭代次數(shù)就能得到高信噪比的分離結(jié)果，說明該算法有很好的實(shí)用性。同時對比分析可以發(fā)現(xiàn)，兩個算例的結(jié)果顯示高階Radon變換對有信號具有更好的保幅效果，分離效果也要優(yōu)于傳統(tǒng)Radon變換。

圖8 SEAM_2B模型兩種分離方法信噪比曲線Fig.8 S/N ratio for SEAM_2B example

[1]Zhou Y, Gao J, Chen W, et al.Seismic Simultaneous Source Separation via Patchwise Sparse Representation[J].IEEE Transactions on Geoscience and Remote Sensing,2016, 54(9)∶5271-5284.

[2]Li C, Mosher C C, Morley L C, et al.Joint source deblending and reconstruction for seismic data[M]//SEG Technical Program Expanded Abstracts 2013.Society of Exploration Geophysicists,2013∶ 82-87.

[3]Beasley C J, Chambers R E, Jiang Z.A new look at simultaneous sources[M]//SEG Technical Program Expanded Abstracts 1998.Society of Exploration Geophysicists, 1998∶133-135.

[4]Berkhout A J.Changing the mindset in seismic data acquisition[J].The Leading Edge,2008, 27(7)∶ 924-938.

[5]Godwin J, Sava P.2010.Simultaneous source imaging by amplitude encoding[J].CWP report,2010, 645.

[6]Moore I.Simultaneous Sources-Continuous Acquisition and Processing for Marine Data[C]//76th EAGE Conference and Exhibition 2014.2014.

[7]Liu Z, Wang B, Specht J, et al.Enhanced adaptive subtraction method for simultaneous source separation[M]//SEG Technical Program Expanded Abstracts 2014.Society of Exploration Geophysicists, 2014∶115-119.

[8]韓立國,譚塵青,呂慶田等.基于迭代去噪的多源地震混合采集數(shù)據(jù)分離[J].地球物理學(xué)報, 2013,56(7)∶2402-2412.

[9]Moore I, Dragoset B, Ommundsen T, et al.Simultaneous source separation using dithered sources[M]//SEG Technical Program Expanded Abstracts 2008.Society of Exploration Geophysicists,2008∶2806-2810.

[10]Mahdad A.Deblending of seismic data[D].TU Delft, Delft University of Technology,2012.

[11]Chen Y, Fomel S, Hu J.Iterative deblending of simultaneoussource seismic data using seislet-domain shaping regularization[J].Geophysics, 2014,79(5)∶ V179-V189.

[12]Zu S, Zhou H, Chen Y, et al.A periodically varying code for improving deblending of simultaneous sources in marine acquisition[J].Geophysics, 2016,81(3)∶ V213-V225.

[13]Sacchi M D.Sparse inversion of the Radon coefficients in the presence of erratic noise with application to simultaneous seismic source processing[C]//Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on.2014∶2386-2389.

[14]唐歡歡,毛偉建.3D高階拋物Radon變換地震數(shù)據(jù)保幅重建[J].地球物理學(xué)報,2014,57(9)∶2918-2927.

[15]Fomel S.Shaping regularization in geophysical-estimation problems[J].Geophysics, 2007,72(2)∶ R29-R36.

[16]薛亞茹, 唐歡歡, 陳小宏.高階高分辨率Radon變換地震數(shù)據(jù)重建方法[J].石油地球物理勘探, 2014, 49(1)∶95-100.

[17]Berkhout A J, Blacquière G, Verschuur D J.The concept of double blending∶ Combining incoherent shooting with incoherent sensing[J].Geophysics, 2009,74(4)∶ A59-A62.

[18]祖紹環(huán),周輝,陳陽康，等.2016.不規(guī)則采樣的多震源數(shù)據(jù)整形正則化分離方法[J].石油地球物理勘探,2016,51(2)∶247-253.

[19]鞏向博.高精度.Radon變換及其應(yīng)用研究[D].吉林∶吉林大學(xué),2008.

[20]Aster R C, Borchers B, Thurber C H.Parameter Estimation and Inverse Problems[M].New York∶ Academic Press,2011.

[21]王典,劉財,劉洋，等.基于提升算法和百分位數(shù)軟閾值的小波去噪技術(shù)[J].地球物理學(xué)進(jìn)展,2008,23(4)∶1124-1130.

[22]Yang P, Gao J, Chen W.On analysis-based two-step interpolation methods for randomly sampled seismic data[J].Computers & Geosciences, 2013,51∶ 449-461.

[23]Hennenfent G, Herrmann F J.Seismic denoising with nonuniformly sampled curvelets[J].Computing in Science &Engineering, 2006,8(3)∶ 16-25.