石 穎，李 瑩，王維紅，陸加敏

（1.東北石油大學(xué) 地球科學(xué)學(xué)院，黑龍江 大慶 163318； 2.大慶油田有限責(zé)任公司 勘探開發(fā)研究院，黑龍江 大慶163712）

線性Radon變換噪音壓制法及其在古龍斷陷中的應(yīng)用

石 穎1，李 瑩1，王維紅2，陸加敏2

（1.東北石油大學(xué) 地球科學(xué)學(xué)院，黑龍江 大慶 163318； 2.大慶油田有限責(zé)任公司 勘探開發(fā)研究院，黑龍江 大慶163712）

針對(duì)松遼盆地北部古龍斷陷地震資料信噪比低、線性干擾強(qiáng)特點(diǎn)，提出應(yīng)用線性Radon變換進(jìn)行疊前線性噪音壓制的預(yù)處理方法，Radon變換可在炮集和CMP道集上進(jìn)行運(yùn)算，算法簡(jiǎn)單，易于編程實(shí)現(xiàn)，其積分路徑的特點(diǎn)適合線性噪音壓制.模擬數(shù)據(jù)和實(shí)際地震資料應(yīng)用結(jié)果表明，線性Radon變換法能夠?qū)崿F(xiàn)保幅的線性噪音壓制，是疊前提高地震資料信噪比的實(shí)用方法，在地震資料預(yù)處理中具有應(yīng)用前景.

Radon變換；線性噪音；壓制；信噪比；古龍斷陷；地震資料

0 引言

當(dāng)前油氣勘探以復(fù)雜構(gòu)造和巖性油氣藏為主，要求高保真和高信噪比的地震成果數(shù)據(jù)體，以利于進(jìn)行綜合地質(zhì)研究和勘探部署.地震資料噪音主要包括隨機(jī)噪音和規(guī)則噪音兩類，其中規(guī)則噪音中的線性噪音在地震資料中普遍存在，往往具有較強(qiáng)的能量，所以線性噪音壓制效果在一定程度上決定著成像數(shù)據(jù)體質(zhì)量[].

根據(jù)有效信號(hào)與線性噪聲的特征差異，研究并提出了不同的線性噪音壓制方法，主要包括F-K濾波法、徑向?yàn)V波法、K-L變換法和Radon變換法等[2-6].其中F-K濾波法是頻率域方法，易于實(shí)現(xiàn)；但是在有效波和線性干擾之間的視速度相近時(shí)，難以有效實(shí)現(xiàn)有效波和噪音的分離.徑向?yàn)V波法根據(jù)選定的線性同相軸的視速度和濾波的徑向道數(shù)，實(shí)現(xiàn)線性干擾的壓制；但是該方法需要過多的人為參與，對(duì)于海量地震數(shù)據(jù)處理，計(jì)算效率低，計(jì)算結(jié)果在很大程度上取決于濾波道數(shù)的選擇，且地震數(shù)據(jù)體中往往有較多的殘余噪音.K-L變換法是基于目標(biāo)統(tǒng)計(jì)特征的正交變換，能夠?qū)崿F(xiàn)隨機(jī)噪音、線性噪音和多次波等規(guī)則噪音壓制；但是該方法噪音壓制效果取決于本征值選擇，本征值選擇需要反復(fù)的數(shù)值試驗(yàn)，若本征值選擇不當(dāng)將造成噪音的大量剩余，同時(shí)有效波的信號(hào)也受到一定的損失.

線性Radon變換法算法簡(jiǎn)單，在地震資料處理中被廣泛應(yīng)用，如資料預(yù)處理中的線性噪音壓制、地震數(shù)據(jù)的插值和重建、VSP資料處理中的波場(chǎng)分離和疊前偏移領(lǐng)域等[6-9].線性Radon變換法能夠較好地實(shí)現(xiàn)有效波和相干線性噪音的分離，計(jì)算效率高，線性同相軸在Radon正變換域可變換為點(diǎn)，濾波函數(shù)簡(jiǎn)單，同時(shí)經(jīng)噪音模型反變換后，從原始數(shù)據(jù)中將線性噪音減去，對(duì)有效波的振幅沒有任何傷害，是一種高效的相對(duì)保幅的線性噪音壓制技術(shù)[10-12].依據(jù)線性Radon變換原理，對(duì)大慶油田古龍斷陷區(qū)的地震資料進(jìn)行強(qiáng)線性干擾壓制，獲得良好成像效果，試算結(jié)果表明算法具有計(jì)算效率高、精度高和實(shí)用性強(qiáng)的特點(diǎn).

1 Radon變換原理與實(shí)現(xiàn)

對(duì)于二維地震資料，假設(shè)沿偏移距方向?yàn)榫鶆虿蓸?，則離散線性Radon正變換[9，13]可以寫為

式中：d（xi，t）為x-t域地震數(shù)據(jù)體，共有N 道數(shù)據(jù)，xi為偏移距；m（τ，p）為Radon正變換域數(shù)據(jù)體，τ為地震數(shù)據(jù)中零偏移距道對(duì)應(yīng)的截距時(shí)間；p為射線參數(shù).

對(duì)Radon變換的離散公式（1）兩端進(jìn)行Fourier變換，對(duì)任一頻率分量進(jìn)行計(jì)算[9，11]：

式中：M（p，ω），D（xi，ω）分別為m（τ，p）和d（xi，t）相對(duì)應(yīng)的Fourier變換域的表示形式.將方程（2）寫成矩陣形式：

式中：d為時(shí)間空間域地震數(shù)據(jù)矩陣；m為Radon域矩陣；L為數(shù)據(jù)變換矩陣，L中的元素可表示為

一般來說，方程（3）是欠定或超定的，采用最小二乘方法求其最優(yōu)解，其一般形式為

式中：LH為矩陣L的共軛轉(zhuǎn)置；m，d分別為模型空間（Radon正變換域）和數(shù)據(jù)空間（時(shí)間空間域）的任一頻率成分的向量.

式（5）和式（3）共同構(gòu)成線性Radon正反變換對(duì).在離散計(jì)算時(shí)，Radon變換進(jìn)行參數(shù)的合理采樣.線性Radon變換的參數(shù)τ和時(shí)間域的參數(shù)t的采樣相等，其射線參數(shù)p的采樣間隔Δp為

最大射線參數(shù)pmax的表達(dá)式為

式（6-7）中：fmax為地震資料的最大有效頻率；Xmax為原始數(shù)據(jù)單炮或CMP道集的最大偏移距；Δx為道集中的道間距[14].

線性Radon變換求解方法是矩陣求逆，式（5）可能是超定方程（未知數(shù)個(gè)數(shù)少于方程個(gè)數(shù)），也可能是欠定方程（未知數(shù)個(gè)數(shù)多于方程個(gè)數(shù)），對(duì)于文中線性噪音壓制問題，單個(gè)頻率成分所形成的矩陣為Toeplitz矩陣，可以應(yīng)用共軛梯度法、Cholesky分解法或Levinson遞推法進(jìn)行求解，其中Levinson遞推算法的計(jì)算效率最高[13].

2 模擬數(shù)據(jù)試驗(yàn)

為驗(yàn)證Radon變換線性噪音壓制方法的有效性，應(yīng)用理論模擬的單炮數(shù)據(jù)進(jìn)行噪音壓制試算.數(shù)值模擬含一個(gè)有效波同相軸和一個(gè)線性同相軸的單炮數(shù)據(jù)（見圖1（a））.由圖1（a）可以看出，模擬線性噪音的能量很強(qiáng).模擬采用主頻為20 Hz的Ricker子波，共84道地震記錄，地震道采樣間隔為40 m，時(shí)間方向采樣間隔為4 ms，最大偏移距為3 320 m.采用式（5）計(jì)算最小二乘Radon域的正變換域見圖1（b）.由圖1（b）可以看出，模擬炮記錄中線性同相軸在線性Radon正變換域表現(xiàn)為相對(duì)聚焦的剪刀狀特征，而有效波表現(xiàn)為橢圓形式.在Radon正變換的模型空間域設(shè)計(jì)濾波函數(shù)，可以容易地濾除點(diǎn)狀形式的線性同相軸，將濾波的結(jié)果進(jìn)行反變換，即可得到線性同相軸（見圖1（c）），將圖1（c）從原始地震數(shù)據(jù)中減去，得到線性噪音壓制后的結(jié)果剖面（見圖1（d））.因此，線性Radon變換法可以有效壓制地震數(shù)據(jù)中的線性同相軸.

圖1 Radon變換線性噪音壓制模擬數(shù)據(jù)試驗(yàn)

線性Radon變換法的基本思想是對(duì)某個(gè)函數(shù)在給定的線性路徑上進(jìn)行積分運(yùn)算，基于平面波分解原理，也考慮地震波場(chǎng)的性質(zhì).根據(jù)線性Radon域的特點(diǎn)，可以設(shè)計(jì)濾波器，在Radon正變換域?yàn)V波，以切除不需要的同相軸或者相干噪音，進(jìn)而大幅提高地震資料的信噪比，表明Radon變換法在線性噪音壓制計(jì)算中的有效性和實(shí)用性.

3 現(xiàn)場(chǎng)資料應(yīng)用

松遼盆地北部深層天然氣勘探在徐家圍子斷陷取得成功，獲得儲(chǔ)量的巨大發(fā)現(xiàn)[15]，當(dāng)前大慶油田探區(qū)的古龍斷陷是重要的儲(chǔ)量接替區(qū).該區(qū)針對(duì)深層火山巖勘探的研究程度相對(duì)較低，其中影響因素主要包括火山巖目的層埋藏深、地震資料品質(zhì)差和信噪比低等.為提高古龍斷陷深層復(fù)雜構(gòu)造和復(fù)雜地質(zhì)體的成像精度，采用Kirchhoff積分的疊前深度偏移技術(shù)[16].該技術(shù)要求輸入信噪比高的疊前預(yù)處理數(shù)據(jù)體，所以疊前精細(xì)噪音壓制是提高成像精度的關(guān)鍵環(huán)節(jié).

以古龍斷陷的葡南工區(qū)為例，原始地震資料分析表明，該地區(qū)地表?xiàng)l件復(fù)雜，折射干擾嚴(yán)重，發(fā)育兩組折射波，速度分別為1 750 m/s和2 200 m/s.線性噪音具有高能量、強(qiáng)振幅和分布范圍廣的特點(diǎn)，且有較強(qiáng)的規(guī)律性.全區(qū)測(cè)線或多或少含有這種干擾，幾乎分布于整個(gè)單炮記錄，從而掩蓋有效波組.一般而言，根據(jù)線性干擾與有效波在速度、頻率、時(shí)空上的差別，進(jìn)行線性干擾的識(shí)別和壓制.在古龍斷陷區(qū)，線性干擾波具有相對(duì)穩(wěn)定的視速度，所以可以應(yīng)用線性Radon變換法實(shí)現(xiàn)線性干擾的有效壓制.另外，在實(shí)現(xiàn)過程中同時(shí)采用噪音壓制的減去法，使得有效波的保幅性得到很大程度的增強(qiáng).

Radon變換法線性干擾壓制前后的單炮效果見圖2.由圖2（a）可以看出，線性噪音呈排狀自上而下分布，幾乎掩蓋所有的有效波；由圖2（b）可以看出，壓制干擾后的炮記錄較為清晰地顯示拋物線型的有效波同相軸，單炮記錄的信噪比得到提高；由圖2（c）可以看出，該方法能夠有效壓制古龍斷陷的強(qiáng)能量線性噪音，同時(shí)有效波的振幅不受損失，即文中方法具有很好的保幅效果.

圖2 線性干擾壓制前后炮記錄對(duì)比

線性干擾壓制前后疊加剖面的效果見圖3.由圖3可以看出，經(jīng)過線性噪音壓制后，地震剖面信噪比明顯提高.通過圖2和圖3的對(duì)比，無論是從單炮還是從疊加剖面都可以看出，線性干擾得到有效壓制，突出有效波能量，地震資料的信噪比大幅提高，為疊前深度偏移成像提供資料保障.古龍斷陷的實(shí)際地震資料線性噪音壓制表明，文中Radon變換法具有很強(qiáng)的實(shí)用性，能夠在信噪比低的同類地震資料處理中推廣應(yīng)用.

圖3 線性干擾壓制前后疊加剖面對(duì)比

4 結(jié)論

（1）線性Radon變換法是對(duì)給定函數(shù)沿一定路徑的積分計(jì)算，實(shí)際上也考慮地震波場(chǎng)的性質(zhì)，其變換形式可從平面波分解的原理進(jìn)行求解和推導(dǎo).線性Radon變換原理簡(jiǎn)單，在地震資料預(yù)處理中應(yīng)用廣泛.

（2）線性Radon變換法能夠較好地實(shí)現(xiàn)有效波和相干線性噪音的分離.該方法計(jì)算效率高，可以用于炮集或CMP道集數(shù)據(jù)，線性同相軸在Radon正變換域可變換為點(diǎn)，濾波函數(shù)簡(jiǎn)單，同時(shí)經(jīng)噪音模型反變換后，從原始數(shù)據(jù)中減去，能夠?qū)崿F(xiàn)線性噪音的保幅壓制.

（3）松遼盆地古龍斷陷地震資料的線性噪音較為發(fā)育，在部分炮記錄上具有很強(qiáng)的能量，應(yīng)用線性Radon變換法進(jìn)行線性噪音壓制，實(shí)際地震數(shù)據(jù)線性噪音壓制結(jié)果表明，文中給出的線性噪音壓制方法具有有效性和實(shí)用性，對(duì)后續(xù)地震波精確成像具有指導(dǎo)意義.

[1]胡天躍.地震資料疊前去噪技術(shù)的現(xiàn)狀與未來[J].地球物理學(xué)進(jìn)展，2002，17（2）：35-40.

[2]閆立志，景新義，李剛.F-K濾波在噪音減去法中的應(yīng)用[J].海洋地質(zhì)動(dòng)態(tài)，2006，22（10）：28-32.

[3]余波，黃中玉，談大龍，等.徑向道濾波法去線性干擾[J].石油物探，2005，44（2）：109-112.

[4]黃雪繼，劉來祥，王永勝.分頻徑向道中值濾波在地震資料處理中的應(yīng)用[J].物探與化探，2012，36（2）：317-320.

[5]劉洪林，張春堂，朱秋影，等.K-L變換在地震資料去噪中的應(yīng)用[J].大慶石油學(xué)院學(xué)報(bào)，2007，31（4）：19-21.

[6]吳律，武克奮，孫立.變換方法及其在地震資料處理中的應(yīng)用[J].石油物探，1996，1（1）：37-52.

[7]Kabir M N，Verschuur D J.Restoration of missing offsets by parabolic Radon transform[J].Geophys.Prosp.，1995，43：347-368.

[8]王維紅，劉洪.拋物 Radon變換法近偏移距波場(chǎng)外推[J].地球物理學(xué)進(jìn)展，2005，20（2）：289-293.

[9]Schultz P S，Claerbout J F.Velocity estimation and downward continuation by wavefront synthesis[J].Geophysics，1978，43（3）：691-714.

[10]林文，魏大力，王建民，等.Radon變換多次波壓制方法及應(yīng)用研究[J].物探化探計(jì)算技術(shù)，2009，31（4）：344-348.

[11]王維紅，首皓，劉洪，等.線性同相軸波場(chǎng)分離的高分辨率τ-p變換法[J].地球物理學(xué)進(jìn)展，2006，21（1）：74-78.

[12]Cambois G.Preserved amplitude processing in the presence of noise[C].66th Annual Internat.Mtg.，Soc.Expl.Geophys.，SEG Expanded Abstracts，1996：1595-1598.

[13]Kostov C.Toeplitz structure in slant-stack inversion[C].60th Annual Internat.Mtg.，Soc.Expl.Geophys.，SEG Expanded Abstracts，1990：1618-1621.

[14]張旭東.幾種多次波壓制方法討論[J].內(nèi)蒙古石油化工，2011（1）：84-85.

[15]胡明，遲建功，呂延防，等.徐家圍子斷陷徐深21區(qū)塊蓋層封閉能力演化與天然氣成藏的關(guān)系[J].大慶石油學(xué)院學(xué)報(bào)，2010，34（1）：24-28.

[16]王維紅，林春華，陳志德，等.古龍斷陷深層火山巖地震資料成像方法及應(yīng)用研究[J].地球物理學(xué)報(bào)，2011，54（2）：310-319.

Prestack reverse time migration based on GPU parallel accelerating algorithm/2012，36（4）：111-115

SHI Ying1，LU Jia-min2，KE Xuan1，TIAN Dong-sheng1，WANG Fei3
（1.School of Geosciences，Northeast Petroleum University，Daqing，Heilongjiang 163318，China；2.Exploration and Development Research Institute，Daqing Oilfield Co.Ltd.，Daqing，Heilongjiang 163712，China；3.School of Earth Science，China Petroleum University （Beijing），Beijing 102200，China）

In order to improve the complex subsurface imaging accuracy and computational efficiency of the algorithm，this paper presents an algorithm of prestack reverse time migration based on GPU（Graphic Processing Unit）accelerating which can image the underground complex structure effectively and accurately.By two-way wave equation to calculate wave field extrapolation，prestack reverse-time migration can overcome the dip limit，and the imaging algorithm is performed by high order finite difference in the paper.Wave field extrapolation and imaging condition are calculated by GPU parallel accelerating technology，comparing to conventional algorithm，its computation efficiency has been greatly improved，and it meets large amount of computation requirement in prestack reverse-time migration.The random boundary condition approach is adopted to obtain wavefield information，which reduces the memory demand but sacrifices the computation cost，and it solves the massy memory problem in reverse time migration.The tests on model illustrate that this approach can imaging complicated geological body efficiently and precisely.

reverse time migration；GPU；acceleration；high order finite difference；random boundary condition；complicated structure

TE132.1

A

2095-4107（2012）04-0116-05

DOI 10.3969/j.issn.2095-4107.2012.04.021

2012-05-07；編輯：任志平

國家自然科學(xué)基金青年基金項(xiàng)目（41104088，41004057）；國家“863”高技術(shù)研究發(fā)展計(jì)劃項(xiàng)目（2012AA061202）；中國博士后科學(xué)基金項(xiàng)目（2011M501009）；中國石油科技創(chuàng)新基金項(xiàng)目（2011D-5006-0304）；黑龍江省教育廳科學(xué)技術(shù)研究項(xiàng)目（12511025）；黑龍江省博士后科學(xué)基金項(xiàng)目（LBH-Z11272）

石 穎（1976-），女，博士，副教授，主要從事地震資料處理方面的研究.