λ－f域拋物Radon變換多次波壓制方法

王維紅，張 振，石 穎，李 瑩

（東北石油大學(xué)地球科學(xué)學(xué)院，黑龍江大慶 163318）

λ－f域拋物Radon變換多次波壓制方法

王維紅，張 振，石 穎，李 瑩

（東北石油大學(xué)地球科學(xué)學(xué)院，黑龍江大慶 163318）

空間截?cái)嘈?yīng)降低傳統(tǒng)Radon變換壓制多次波的精度.提出一種λ－f域拋物Radon變換方法，與傳統(tǒng)最小二乘拋物Radon變換相比，地震數(shù)據(jù)變換到λ－f域后表現(xiàn)為直線的映射，更便于設(shè)計(jì)濾波函數(shù)，從而有效消除傳統(tǒng)多次波壓制過(guò)程中空間截?cái)嘈?yīng)的影響.另外，新變量λ的引入，使Radon正反算子成為與頻率無(wú)關(guān)的矩陣，可大幅度提高計(jì)算效率.理論模型和實(shí)際海洋地震數(shù)據(jù)多次波壓制試算結(jié)果表明，提出的λ－f域拋物Radon變換法可有效實(shí)現(xiàn)多次波壓制，提高地震數(shù)據(jù)的信噪比，具有較強(qiáng)的實(shí)用性.

Radon變換；多次波壓制；λ－f域；計(jì)算效率

0 引言

多次波壓制是海洋地震資料處理中的重要預(yù)處理環(huán)節(jié)［1－4］，因?yàn)槎啻尾ǖ拇嬖诮档偷卣饠?shù)據(jù)的分辨率，影響地震數(shù)據(jù)的成像質(zhì)量［5－7］和后續(xù)地震解釋的精度［8］，有效地壓制多次波在地震數(shù)據(jù)處理中顯得尤為重要.Radon變換是疊前地震數(shù)據(jù)處理中壓制多次波的有效方法，人們研究Radon變換壓制多次波的方法，Hampson D［9］提出頻率域最小二乘拋物Radon變換方法，但是在Radon域存在剪刀狀“拖尾”現(xiàn)象，致使濾波函數(shù)設(shè)計(jì)困難，從而導(dǎo)致多次波壓制效果不理想.Sacchi M等［10］改進(jìn)Hampson方法，提出頻率域高分辨率拋物Radon變換方法；該方法在處理簡(jiǎn)單數(shù)據(jù)時(shí)可以有效地解決變換域“拖尾”現(xiàn)象，但對(duì)復(fù)雜介質(zhì)的地震數(shù)據(jù)，無(wú)法得到高精度多次波壓制效果.也有學(xué)者提出時(shí)間域Radon變換方法［11－12］彌補(bǔ)頻率域Radon分辨率較低的缺陷，但時(shí)變矩陣的求解無(wú)法應(yīng)用Levinson遞推快速計(jì)算方法，導(dǎo)致計(jì)算效率大幅度降低.

筆者提出一種基于λ－f域拋物Radon變換壓制多次波的方法，通過(guò)變量λ的引入，該方法的Radon變換算子與頻率無(wú)關(guān)，從而比傳統(tǒng)的Radon變換具有更高的計(jì)算效率［13］.同時(shí)，經(jīng)λ－f域拋物Radon變換，拋物線型的地震同相軸在λ－f域變?yōu)榫哂幸欢ㄐ甭实闹本€，在一定程度上避免多次波壓制過(guò)程中“拖尾”現(xiàn)象的干擾［14］.經(jīng)過(guò)理論模型及實(shí)際地震資料的地震道集數(shù)據(jù)多次波壓制試算，該方法具有較好的多次波壓制效果.

1 最小二乘拋物Radon變換

拋物Radon正反變換時(shí)間域離散積分形式為

式（1－2）中：d（x，t）為時(shí)間—炮檢距域地震數(shù)據(jù)；u（q，τ）為對(duì)應(yīng)的Radon域的模型空間；q為曲率參數(shù)；τ為Radon域零偏移距截距時(shí)間；t為雙程旅行時(shí)；Nx為地震數(shù)據(jù)的道數(shù)；Nq為曲率參數(shù)q的個(gè)數(shù)；xn為偏移距.

由于時(shí)間t與截距時(shí)間τ存在線性關(guān)系，可經(jīng)由頻率域高效計(jì)算拋物Radon變換，所以方程式（1）、（2）沿時(shí)間方向做傅里葉變換后，頻率域形式表示為

式（3－4）中：u（q，f）、d（x，f）分別為u（q，τ）和d（x，t）對(duì)應(yīng)的傅里葉變換；f為頻率；j為虛數(shù)單位.對(duì)于單個(gè)頻率成分，方程式（4）的矩陣形式為

其中，Radon算子為

由于L與LH不是互逆的，一般不直接采用算子LH計(jì)算Radon變換，而是采用最小二乘方法定義拋物Radon正變換，即

式中：I為單位矩陣；μ通常取LHL矩陣主對(duì)角線數(shù)據(jù)的0.01倍［15］，增強(qiáng)矩陣求解的穩(wěn)定性.

式（5）、（7）即為正反Radon變換的關(guān)系式，由于式（5）、（7）中的Radon算子與頻率有關(guān)，所以在式（5）、（7）求解過(guò)程中，每計(jì)算一個(gè)頻率分量，需要計(jì)算一次Radon算子及算子的逆，降低計(jì)算效率.

2 λ－f域拋物Radon變換

定義新變量λ＝fq，式（4）在λ－f域的Radon變換定義為

式中：M為變量λ的個(gè)數(shù).

與式（7）對(duì)應(yīng)的λ－f域拋物Radon正變換為

其中，Radon算子L（λ）表達(dá)式為

λ－f域拋物Radon變換的算子L（λ）是關(guān)于偏移距參數(shù)xn及λ的復(fù)數(shù)矩陣，與頻率無(wú)關(guān)，在整個(gè)λf域拋物Radon正反變換過(guò)程中僅需計(jì)算一次，相比于傳統(tǒng)τ－q域拋物Radon變換，計(jì)算效率有很大程度的提高.并且LHL為主對(duì)角線元素相等的Toeplitz矩陣，所以Levinson遞推法快速求解依然適用.

理論模擬炮記錄的Radon變換域見(jiàn)圖1.圖1（a）為含有3個(gè)同相軸的單炮地震記錄，共100道，道間距為20m，時(shí)間采樣間隔為4ms，每道600個(gè)采樣點(diǎn).對(duì)圖1（a）的地震數(shù)據(jù)采用最小二乘τ－q域拋物Radon正變換，理論上在Radon域得到3個(gè)點(diǎn)，然而實(shí)際上是3個(gè)含有很大假象的剪刀狀的結(jié)果（“拖尾”現(xiàn)象）（見(jiàn)圖1（b））.這是因?yàn)樵赗adon變換離散運(yùn)算中，地震數(shù)據(jù)的有限孔徑引起截?cái)嘈?yīng).在Radon變換壓制多次波過(guò)程中，由于存在截?cái)嘈?yīng)的干擾，有效波與多次波在Radon域的假象可能交織在一起，使得有效波與多次波的假象得不到完全分離，進(jìn)而導(dǎo)致波場(chǎng)分離不徹底.

為消除空間截?cái)嘈?yīng)的影響，對(duì)圖1（a）的地震數(shù)據(jù)進(jìn)行λ－f域拋物Radon變換，得到λ－f域剖面（見(jiàn)圖1（c）），地震數(shù)據(jù)在經(jīng)過(guò)λ－f域拋物Radon變換后，能量分布在以q為斜率的直線上，并未因空間截?cái)嘈?yīng)在λ－f域產(chǎn)生假象，在一定程度上避免“拖尾”現(xiàn)象的干擾.

3 離散采樣及濾波函數(shù)

3.1 離散采樣

為避免在λ－f域拋物Radon變換過(guò)程中出現(xiàn)假頻，對(duì)采樣間隔參數(shù)Δλ和λ的范圍進(jìn)行合理離散采樣，其中：

式（11－12）中：λmax、λmin分別為λ的最大和最小值；xmax、xmin分別為數(shù)據(jù)的最大和最小偏移距；ΔTmax、ΔTmin分別為地震數(shù)據(jù)在時(shí)空域的最大與最小動(dòng)校正時(shí)差；fmax為數(shù)據(jù)最大有效頻率.

3.2 濾波函數(shù)

根據(jù)λ定義，地震數(shù)據(jù)d在經(jīng)過(guò)λ－f域拋物Radon變換后，一次波與多次波呈現(xiàn)過(guò)原點(diǎn)以不同曲率參數(shù)q為斜率的直線，引入濾波矩陣G（λ，f），對(duì)有效波和多次波進(jìn)行分離.

式中：ΔT0為有效波與多次波分離所對(duì)應(yīng)的動(dòng)校正時(shí)差.有效波在λ－f域可表示為

變換到x－f域的有效波可表示為

4 數(shù)據(jù)試算

4.1 理論模型

為驗(yàn)證λ－f域拋物Radon變換在壓制多次波中的有效性，合成含表面多次波的模型數(shù)據(jù)，見(jiàn)圖2（a）.理論模型數(shù)據(jù)共100道，道間距為20m，時(shí)間采樣間隔為4ms，每道601個(gè)采樣點(diǎn).分段部分動(dòng)校正后的剖面見(jiàn)圖2（b），一次波同相軸向上彎曲，曲率參數(shù)q為負(fù)值，為校正過(guò)量；多次波同相軸向下彎曲，曲率參數(shù)q為正值，為校正不足.地震數(shù)據(jù)在經(jīng)過(guò)部分動(dòng)校正后，有效波與多次波在變換域內(nèi)得到更加有效的區(qū)分.對(duì)圖2（b）的地震數(shù)據(jù)應(yīng)用文中方法得到λ－f域剖面（見(jiàn)圖2（c））.根據(jù)λ定義，一次波與多次波在λ－f域剖面上為經(jīng)過(guò)原點(diǎn)且分布在原點(diǎn)兩側(cè)不同斜率的直線上.根據(jù)式（13）選取合適的濾波算子對(duì)圖2（c）進(jìn)行多次波切除，得到濾波后λ－f域剖面（見(jiàn)圖2（d））.λ－f域剖面變換到時(shí)空域得到的單炮記錄見(jiàn)圖2（e），壓制多次波后炮記錄只有5個(gè)有效波同相軸，在圖2（a）的其他多次波得到有效壓制，即地震資料的信噪比得到有效提高.對(duì)理論模型數(shù)據(jù)應(yīng)用最小二乘拋物Radon變換方法壓制多次波后的剖面見(jiàn)圖2（f），模擬的地震數(shù)據(jù)在1.8s處，有效波同相軸與多次波同相軸在時(shí)間域難以分開(kāi)，在Radon變換域由于受空間截?cái)嘈?yīng)的影響，兩者的像在變換域也難以有效分開(kāi)，所以影響濾波函數(shù)的設(shè)計(jì)，使最小二乘拋物Radon變換多次波壓制的效果不理想，即多次波有較大殘余能量，同時(shí)有效波的振幅也受到一定程度的損失.

圖2 合成炮記錄的多次波壓制效果Fig.2 Multiple suppression result on synthetic shot record

4.2 實(shí)際地震數(shù)據(jù)多次波壓制

為進(jìn)一步測(cè)試文中方法在壓制多次波中的穩(wěn)定性和效果，對(duì)海洋地震數(shù)據(jù)進(jìn)行多次波壓制試算（見(jiàn)圖3）.圖3（a）為某工區(qū)海洋地震資料的一個(gè)CMP道集，數(shù)據(jù)共60道，道間距為50m，時(shí)間采樣間隔為4 ms，時(shí)間方向8×103個(gè)采樣點(diǎn).圖3（a）箭頭在0.8s處的一階表面多次波較為發(fā)育，應(yīng)用λ－f域拋物Radon變換方法對(duì)該地震數(shù)據(jù)進(jìn)行濾波處理，得到多次波壓制后的道集記錄見(jiàn)圖3（b），能量很強(qiáng)的一階表面多次波得到有效的壓制，一次波的同相軸清晰可見(jiàn)，地震數(shù)據(jù)的信噪比得到有效提高，可為后續(xù)地震資料成像處理提供數(shù)據(jù)基礎(chǔ)，實(shí)際數(shù)據(jù)的計(jì)算表明文中方法在壓制海洋地震資料的多次波方面具有有效性和實(shí)用性.

圖3 實(shí)際數(shù)據(jù)多次波壓制結(jié)果Fig.3 Multiple suppression result on field data

5 結(jié)論

（1）λ－f域拋物Radon變換通過(guò)引入新變量λ，Radon算子的特點(diǎn)和傳統(tǒng)最小二乘方法不同，其正反算子僅需計(jì)算一次，大幅提高該方法的計(jì)算效率.

（2）拋物線型同相軸在λ－f域拋物Radon變換域?yàn)檫^(guò)原點(diǎn)的具有一定斜率的直線，在一定程度上避免傳統(tǒng)拋物Radon變換壓制多次波過(guò)程中空間截?cái)嘈?yīng)的影響，并且其多次波壓制濾波函數(shù)的設(shè)計(jì)也更為簡(jiǎn)單.

（3）應(yīng)用λ－f域拋物Radon變換可有效提高地震數(shù)據(jù)的信噪比，具有穩(wěn)定、有效和實(shí)用的特點(diǎn).

（References）：

［1］ 薛東川.幾種疊前逆時(shí)偏移成像條件的比較［J］.石油地球物理勘探，2013，48（2）：222－227.

Xue Dongchuan.Comparison of several imaging condition for prestack reverse－time migration［J］.Oil Geophysical Prospecting，2013，48（2）：222－227.

［2］ 馬繼濤，Sen M K，陳小宏.平面波域反數(shù)據(jù)處理壓制多次波方法研究［J］.地球物理學(xué)報(bào)，2009，52（3）：808－816.

Ma Jitao，Sen M K，Chen Xiaohong.Multiple attenuation using inverse data processing in the plane－wave domain［J］.Chinese Journal of Geophysics，2009，52（3）：808－816.

［3］ 劉伊克，Sun Hongchuan，常旭.基于波射線路徑偏移壓制多次波［J］.地球物理學(xué)報(bào)，2004，47（4）：697－701.

Liu Yike，Sun Hongchuan，Chang Xu.Multiple removal by wavepath migration［J］.Chinese Journal of Geophysics，2004，47（4）：697－701.

［4］ 熊登，趙偉，張劍鋒.混合域高分辨率拋物Radon變換及在衰減多次波中的應(yīng)用［J］.地球物理學(xué)報(bào)，2009，52（4）：1068－1077.

Xiong Deng，Zhao Wei，Zhang Jianfeng.Hybrid－domain high－resolution parabolic Radon transform and its application to demultiple［J］.Chinese Journal of Geophysics，2009，52（4）：1068－1077.

［5］ 郭雪豹，王建民，王維紅，等.基于GPU并行加速的VSP數(shù)據(jù)逆時(shí)偏移［J］.東北石油大學(xué)學(xué)報(bào)，2014，38（2）：58－62.

Guo Xuebao，Wang Jianmin，Wang Weihong，et al.VSP data reverse－time migration based on GPU parallel acceleration［J］.Journal of Northeast Petroleum University，2014，38（2）：58－62.

［6］ 石穎，陸加敏，柯璇，等.基于GPU并行加速的疊前逆時(shí)偏移方法研究［J］.東北石油大學(xué)學(xué)報(bào)，2012，36（4）：111－115.

Shi Ying，Lu Jiamin，Ke Xuan，et al.Prestack reverse time migration based on GPU parallel accelerating algorithm［J］.Journal of Northeast Petroleum University，2012，36（4）：111－115.

［7］ 田東升，王云專，李義鵬，等.單程和雙程波動(dòng)方程疊前深度偏移方法［J］.東北石油大學(xué)學(xué)報(bào)，2014，38（4）：39－44.

Tian Dongsheng，Wang Yunzhuan，Li Yipeng，et al.One－way and two－way wave equation pre－stack depth migration approaches［J］.Journal of Northeast Petroleum University，2014，38（4）：39－44.

［8］ 張進(jìn)鐸.地震解釋技術(shù)現(xiàn)狀及發(fā)展趨勢(shì)［J］.地球物理學(xué)進(jìn)展，2006，49（2）：578－587.

Zhang Jinduo.Present status and future trend of seismic data interpretation techniques［J］.Progress in Geophysics，2006，49（2）：578－587.

［9］ Hampson D.Inverse velocity stacking for multiple elimination［J］.Journal of the Canadian Society of Exploration Geophysicists，1986，22（1）：44－55.

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

［11］ Trad D，Ulrych J，Sacchi M.Accurate interpolation with high－resolution time－variant Radon transform［J］.Geophysics，2002，67（2）：644－656.

［12］ Ng M，Perz M.High resolution Radon transform in the t－x domain using"intelligent"prioritization of the Gauss－Seidel estimation sequence［J］.G xandd Abra，2004，18（3）：260－263.

［13］ 石穎，張振，李婷婷，等.λ－f域加權(quán)拋物Radon變換地震數(shù)據(jù)重建方法研究［J］.地球物理學(xué)進(jìn)展，2014，29（4）：1752－1757.

Shi Ying，Zhang Zhen，Li Tingting，et al.Investigation on seismic data reconstruction by domain weighted parabolic Radon transform method［J］.Progress in Geophysics，2014，29（4）：1752－1757.

［14］ 謝俊法，孫成禹，韓文功.迭代拋物Radon變換法分離一次波與多次波［J］.石油地球物理勘探，2014，49（1）：76－81.

Xie Junfa，Sun Chengyu，Han Wengong.Iterative parabolic Radon transform for primary and multiple separations［J］.Oil Geophysical Prospecting，2014，49（1）：76－81.

［15］ Schonewille M A，Duijndam A J W.Parabolic Radon transform，sampling and efficiency［J］.Geophysics，2001，66（2）：667－678.

DOI 10.3969／j.issn.2095－4107.2015.01.003

TE132.1，P631

A

2095－4107（2015）01－0017－06

2014－12－01；編輯：陸雅玲

國(guó)家自然科學(xué)基金項(xiàng)目（41474118）；國(guó)家高技術(shù)研究發(fā)展計(jì)劃（863計(jì)劃）項(xiàng)目（2012AA061202）

王維紅（1975－），男，博士，副教授，主要從事地震資料數(shù)字處理方面的研究.