Radon變換及其在地震矢量數(shù)據(jù)處理中的應(yīng)用研究現(xiàn)狀

邱新明,汪 超,苑益軍,王 赟

(1.中國(guó)地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院,北京100083;2.中國(guó)科學(xué)院地球化學(xué)研究所,礦床地球化學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,貴州貴陽(yáng)550081)

邱新明,汪超,苑益軍,等.Radon變換及其在地震矢量數(shù)據(jù)處理中的應(yīng)用研究現(xiàn)狀[J].石油物探,2017,56(6):-914

QIU Xinming,WANG Chao,YUAN Yijun,et al.Research status of the application of Radon transform in seismic vector field processing[J].Geophysical Prospecting for Petroleum,2017,56(6):-914

Radon變換及其在地震矢量數(shù)據(jù)處理中的應(yīng)用研究現(xiàn)狀

邱新明1,2,汪 超2,苑益軍1,王 赟1

(1.中國(guó)地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院,北京100083;2.中國(guó)科學(xué)院地球化學(xué)研究所,礦床地球化學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,貴州貴陽(yáng)550081)

在系統(tǒng)調(diào)研Radon變換及其各種改進(jìn)方法、應(yīng)用效果的基礎(chǔ)上,對(duì)目前Radon變換的研究成果進(jìn)行了總結(jié)。首先,詳細(xì)介紹了Radon變換基本原理及適用性,重點(diǎn)介紹了在Radon變換中采用稀疏正則化約束提高分辨率、在受限模型空間中求解提高計(jì)算效率的一些具體應(yīng)用與改進(jìn);其次對(duì)Radon變換在地震矢量數(shù)據(jù)處理中的應(yīng)用進(jìn)行了概要介紹。調(diào)研發(fā)現(xiàn)Radon變換在地震矢量場(chǎng)處理中的應(yīng)用仍處于探索階段,如何獲得最佳的處理效果并保持其矢量場(chǎng)特性是Radon變換應(yīng)用的發(fā)展方向,針對(duì)多分量地震數(shù)據(jù),研發(fā)或改進(jìn)現(xiàn)有的Radon變換方法使其適用于多分量矢量地震數(shù)據(jù)的處理將是以后該領(lǐng)域攻關(guān)的重點(diǎn)方向之一。

Radon變換;分辨率;矢量場(chǎng);多分量;稀疏正則化

Radon變換在地震信號(hào)處理中備受關(guān)注,被廣泛應(yīng)用于提取波場(chǎng)信息和波場(chǎng)分離等。Radon變換是二維地震數(shù)據(jù)在特定路徑t=φ(p,τ,x)上的曲線積分[1-2],其目的是將數(shù)據(jù)中規(guī)律排列信號(hào)分解為Radon域內(nèi)稀疏的系數(shù),以便識(shí)別和分離信號(hào)。依據(jù)積分路徑,Radon變換可分為線性Radon變換、拋物Radon變換、雙曲Radon變換等[3];依據(jù)算子類型,Radon變換又可分為時(shí)變Radon變換和時(shí)不變Radon變換[3]。線性Radon變換、拋物Radon變換屬于時(shí)不變Radon變換,而雙曲Radon變換屬于時(shí)變Radon變換[4]。

實(shí)際處理過(guò)程中往往先定義Radon反變換,再由反變換推出正變換[5]。Radon反變換公式[6-7]為d=Lm。式中:d為時(shí)空域地震數(shù)據(jù);m為Radon變換系數(shù);L為變換算子。Radon反變換可以理解為模型m在算子L的作用下生成數(shù)據(jù)響應(yīng)d,那么,相應(yīng)的Radon變換是反問(wèn)題:由數(shù)據(jù)求解Radon模型[8]。起初,采用最小平方反演Radon模型,此方法稱為最小二乘Radon變換[9-10]。后來(lái),為提高Radon模型的分辨率,引入稀疏約束,提出了高分辨率Radon變換[6,11]。高分辨率Radon變換是非線性反演問(wèn)題,通常采用迭代重加權(quán)共軛梯度算法求解[12],并在受限模型空間中求解以減少計(jì)算量。近年來(lái),結(jié)合正交多項(xiàng)式變換對(duì)信號(hào)振幅變化信息的描述,提出了高階Radon變換,重建信號(hào)的同時(shí)保持了信號(hào)的AVO特性[13]。

Radon變換在傳統(tǒng)縱波勘探信號(hào)處理中已有很多應(yīng)用,如地震信號(hào)去噪[14-15]、多次波消除[16-17]、地震道重建[18-19]等。線性Radon變換提供了不同速度下縱波、轉(zhuǎn)換橫波、面波的到時(shí)信息,這些信息可用于轉(zhuǎn)換波靜校正[20]、面波頻散曲線提取[21]、鬼波壓制[22]等;拋物Radon變換常用于去除經(jīng)NMO處理后的CMP道集數(shù)據(jù)中的多次波;而雙曲Radon變換常用于地震反射波場(chǎng)分離[17]。隨著多波地震勘探和微震監(jiān)測(cè)技術(shù)在非常規(guī)油氣藏勘探開發(fā)中的迅猛發(fā)展,地震數(shù)據(jù)矢量場(chǎng)處理技術(shù)也得到快速發(fā)展。但是Radon變換在矢量場(chǎng)處理中的應(yīng)用仍處于探索階段,如何通過(guò)Radon變換獲得良好的處理效果并保持信號(hào)的矢量場(chǎng)特性是研究的重點(diǎn)之一。近年來(lái),國(guó)內(nèi)外的很多學(xué)者從不同的角度對(duì)Radon變換進(jìn)行改進(jìn)以適應(yīng)矢量場(chǎng)處理,取得了很多成果[23-24]。本文梳理了Radon變換研究成果,主要討論了其在稀疏正則化、高效計(jì)算等方面的研究歷史與現(xiàn)狀,最后介紹了Radon變換在地震矢量數(shù)據(jù)處理中的應(yīng)用情況,希望研發(fā)出更適于矢量場(chǎng)數(shù)據(jù)處理的新方法。

1 Radon變換及其特點(diǎn)

Radon變換離散形式定義為[25]:

(1)

式中:d(x,φ)為時(shí)空域地震數(shù)據(jù);x是偏移距,單位m;φ是地震波雙程旅行時(shí),單位s;τ為時(shí)間截距,單位s;m(p,τ)為Radon變換系數(shù);p是斜率或者曲率。此定義為伴隨Radon變換,記為:m=LHd,LH稱為伴隨算子。

1.1 不同域Radon變換

Radon變換可以在時(shí)間域、頻率域、時(shí)頻混合域中進(jìn)行[8],三者在分辨率和計(jì)算效率上差異明顯,見(jiàn)表1。

表1 不同域Radon變換的對(duì)比

時(shí)間域與頻率域Radon變換的計(jì)算效率不同:時(shí)間域Radon變換所有時(shí)間分量相互耦合,疊加算子矩陣過(guò)大,致使計(jì)算成本較高;頻率域Radon變換不同頻率分量相互解耦,降低了變換矩陣的維數(shù),且可以只在有效信號(hào)的頻帶內(nèi)計(jì)算,提高了計(jì)算效率[26]。再者Radon模型分辨率也有差異:SACCHI等[11]提出的頻率域高分辨率Radon變換,只在斜率方向(或曲率方向)有約束[27-28],CARY[28]稱其為高分辨率的必要非充分條件;而THORSON[29]提出的時(shí)間域高分辨率Radon變換,在斜率方向和時(shí)間方向都有約束,Radon模型更稀疏,分辨率更高[8,28]。

TRAD等[8]結(jié)合時(shí)間域方法與頻率域方法,提出時(shí)頻混合域Radon變換:在時(shí)間域內(nèi)進(jìn)行稀疏約束,利用迭代重加權(quán)法求高分辨率解;在頻率域內(nèi)進(jìn)行算子矩陣與單頻向量乘積,避免大矩陣運(yùn)算。時(shí)頻混合域方法具有較高的分辨率和計(jì)算效率,受到國(guó)內(nèi)外學(xué)者廣泛關(guān)注[8,26-27]。但是需要指出,時(shí)頻混合域求解方法并不適用于所有積分路徑,例如時(shí)變Radon變換就只能在時(shí)間域求解[26]。

1.2 不同積分路徑Radon變換

另外,Dix動(dòng)校正公式在水平多層介質(zhì)近偏移距假設(shè)下近似成立[35],在遠(yuǎn)偏移距時(shí)不能正確描述反射波的到時(shí),為解決這一問(wèn)題,CASTLE[36]提出了相移雙曲NMO公式:

(2)

圖1 Dix常規(guī)動(dòng)校正雙曲線(a)與相移雙曲線(b)分別與反射同相軸的對(duì)比

BICKEL[37]對(duì)拋物Radon變換變形,提出偽雙曲Radon變換。取路徑:

(3)

式中:z=vresT0,稱為參考深度,與波速和層深度相關(guān),vres為剩余時(shí)差速度;p為慢度。參數(shù)z的取值控制反射同相軸在Radon域的稀疏性,z過(guò)大或過(guò)小都不能將反射同相軸稀疏表示[37]。建立如圖2所示的反射波模型,其偽雙曲Radon變換結(jié)果見(jiàn)圖3,可見(jiàn)此例中只有當(dāng)z取500m左右時(shí),Radon模型是稀疏的;而z取其它值時(shí),Radon模型不稀疏。通過(guò)速度掃描,求取目標(biāo)層反射波的對(duì)應(yīng)z值,可將目標(biāo)層反射波在Radon域稀疏表示。由于CMP道集中不同反射界面的多種反射波進(jìn)行Radon變換所需最優(yōu)的z值都不相同,賈春梅等[3]提出在變換中加入動(dòng)態(tài)時(shí)窗,在不同截距處選取不同的參考深度的方法加以解決。偽雙曲Radon變換積分路徑具有時(shí)不變性,可以在頻率域或時(shí)頻混合域中求解,相對(duì)于前述時(shí)變雙曲Radon變換計(jì)算效率更高。

1.3 高階Radon變換

伴隨Radon變換無(wú)法實(shí)現(xiàn)信號(hào)保真[38]。在標(biāo)準(zhǔn)同相軸假設(shè)下,高分辨率Radon變換有很好的保幅性。但是,Radon變換的分辨率和AVO特性相互矛盾[39-40]。為解決此問(wèn)題,李晶晶等[41]基于拋物Radon變換,根據(jù)AVO特征對(duì)不同偏移距的地震道采用不同的均衡系數(shù),實(shí)現(xiàn)相對(duì)保幅地震道重建;薛亞茹等[13]利用正交多項(xiàng)式譜表征AVO與Radon變換結(jié)合得到高階Radon變換。高階Radon變換包含了不同方向同相軸的振幅信息與振幅變化信息,常用于地震數(shù)據(jù)保幅重建[13,42]。高階Radon變換定義式:

圖2 反射波模型

圖3 對(duì)反射波模型進(jìn)行偽雙曲Radon變換 a z=100m; b z=500m; c z=800m; d z=1300m

(4)

式中:Pj(x)表示j階單位正交多項(xiàng)式;m(j)表示j階的Radon變換系數(shù),j=0,1,2。高階Radon反變換的時(shí)間域形式:

(5)

簡(jiǎn)記為:

(6)

式中:L0,L1,L2分別為疊加算子、梯度算子、曲率算子[13]。與常規(guī)Radon變換求解思路一致,高階Radon變換也采用稀疏反演的方法提高分辨率。不同的是,高階Radon變換采用Radon域的能量來(lái)確定模型加權(quán)矩陣:

(7)

2 Radon變換的正則化問(wèn)題

對(duì)于反問(wèn)題d=Lm,其解不全滿足存在性、唯一性和穩(wěn)定性,故反問(wèn)題是不適定的[11]。借助于某些先驗(yàn)信息而構(gòu)建一個(gè)唯一且穩(wěn)定的解,這種方法就是正則化[11]。最常見(jiàn)的是二次正則化,先驗(yàn)信息以固定形式給出,作用是平滑模型,避免觀測(cè)隨機(jī)誤差被放大。利用二次正則化得到的反問(wèn)題d=Lm的解,就是最小二乘Radon變換。最小二乘Radon變換分辨率不高,并不適用于波場(chǎng)分離。為此,引入模型相關(guān)信息來(lái)得到稀疏正則解,這就是高分辨率Radon變換。由于數(shù)據(jù)的高分辨率Radon變換系數(shù)更稀疏,所以波場(chǎng)分離效果比較理想[43]。本文總結(jié)了高分辨率Radon變換的兩種推導(dǎo)方法。

2.1 基于貝葉斯原理的高分辨率Radon變換

此方法基于貝葉斯原理以概率形式給出先驗(yàn)信息,常用到的先驗(yàn)概率分布為廣義高斯分布和柯西分布[6,44]。

柯西分布:

(8)

式中:σc為柯西分布參數(shù)。

廣義高斯分布:

(9)

式中:p是形狀參數(shù);σp是尺度參數(shù)。當(dāng)p=1時(shí),為L(zhǎng)aplace雙指數(shù)分布,表示稀疏模型;當(dāng)p=2時(shí),為正態(tài)分布,表示平滑模型;當(dāng)p=時(shí),為均勻分布。

SACCHI等[11]推導(dǎo)了基于柯西正則化的高分辨率Radon變換;1996年,SACCHI[44]提出了基于廣義高斯分布和柯西分布的高分辨率Radon變換;TRAD等[6]利用雙指數(shù)分布和柯西分布估計(jì)稀疏模型。

利用貝葉斯原理[11,29],可得模型后續(xù)分布:

(10)

式中:S(m)是模型的全局約束[11];λ0,λ1分別為拉格朗日乘數(shù)。假設(shè)噪聲滿足正態(tài)分布,則(10)式中K為歸一化常數(shù),Cn為噪聲的協(xié)方差矩陣[44],使后續(xù)分布取得最大值的解m稱為最大先驗(yàn)(maximum a posterior,MAP)解[11],也是稀疏解。MAP解等價(jià)于求如下目標(biāo)函數(shù)的最小化問(wèn)題:

(11)

式中:λ=2λ1。對(duì)目標(biāo)函數(shù)(11)最小化,得到方程:

(12)

式中:對(duì)角陣Qm是與先驗(yàn)信息相關(guān)的模型加權(quán)矩陣。當(dāng)先驗(yàn)信息滿足柯西分布時(shí),有:

(13)

當(dāng)先驗(yàn)信息滿足Laplace雙指數(shù)分布時(shí),有:

(14)

由于高分辨率解對(duì)參數(shù)σc,σp的取值不靈敏,可采用估計(jì)方法取值:先將模型向量m歸一化,然后選擇其為0到1之間的值。σc,σp越小,模型越稀疏[6]。

2.2 基于L1范數(shù)稀疏約束的高分辨率Radon變換

在反演理論中,對(duì)于不適定問(wèn)題,為獲得唯一且穩(wěn)定的解,更直接的方法是加入正則項(xiàng)[44]。對(duì)于d=Lm,可構(gòu)造目標(biāo)函數(shù):

(15)

為求稀疏解,正則項(xiàng)和誤差項(xiàng)都使用L1范數(shù)。取p=1,q=1,則:

其中,對(duì)角矩陣Wd稱為數(shù)據(jù)加權(quán)矩陣,反映數(shù)據(jù)的標(biāo)準(zhǔn)差[8];對(duì)角矩陣Wm稱為模型加權(quán)矩陣,決定模型稀疏性[8,46]。Wm和Wd元素分別為:

則d=Lm可表述為新的目標(biāo)函數(shù):

(20)

(20)式優(yōu)化問(wèn)題的解即為Radon變換高分辨率解:

(21)

(22)

式中:參數(shù)λ用于平衡數(shù)據(jù)擬合和模型稀疏的比重[8]。λ越大,模型越稀疏;反之,越接近最小二乘解。

(22)式與貝葉斯原理導(dǎo)出的Laplace雙指數(shù)正則解形式一致。SACCHI[44]認(rèn)為:兩種推導(dǎo)方法實(shí)質(zhì)聯(lián)系密切,且有近似的目標(biāo)函數(shù)?；谪惾~斯原理導(dǎo)出的高分辨率解更能說(shuō)明模型分布特性,但是沒(méi)有解釋清楚“為何先驗(yàn)概率使用廣義高斯分布或者柯西分布而不是其它分布”[44];推導(dǎo)過(guò)程中假定噪聲的協(xié)方差矩陣Cn已知,但實(shí)際數(shù)據(jù)中Cn不容易獲取。基于L1范數(shù)稀疏約束的高分辨率解更實(shí)用,每次迭代容易計(jì)算加權(quán)矩陣Wm和Wd。

3 Radon變換的計(jì)算效率問(wèn)題

計(jì)算高分辨率Radon變換((12)式、(21)和(22)式)可歸納為求解非線性方程:

(23)

式中:系數(shù)矩陣A是稀疏度高的厄米特(Hermite)矩陣,且系數(shù)矩陣A與模型和數(shù)據(jù)相關(guān)。此類方程需要迭代求解,且隨Radon系數(shù)的元素個(gè)數(shù)增多,計(jì)算成本急劇升高。目前此方程的高效解法是在受限模型空間[47]中運(yùn)用迭代重加權(quán)共軛梯度(iteratively reweighted conjugate gradient,IRCG)算法[6,8]。

3.1 迭代重加權(quán)共軛梯度算法

在每次迭代中,更新加權(quán)矩陣,將方程(23)轉(zhuǎn)化為線性方程,并利用共軛梯度法求解。內(nèi)層共軛梯度迭代次數(shù)可由廣義交叉校驗(yàn)(generalized cross validation,GCV)函數(shù)設(shè)定[8,26],外層迭代一般為3～7次[6]。以方程(21)的時(shí)頻混合域解法為例,總結(jié)迭代IRCG算法。

1) 給定初始平衡參數(shù)λ,最大迭代次數(shù)n,初始迭代變量k=1。

2) 二維地震剖面做關(guān)于時(shí)間的傅里葉變換,提取單頻向量,并計(jì)算Radon算子矩陣L及伴隨算子矩陣LH。先求頻率域最小二乘解,并對(duì)其做傅里葉反變換到時(shí)間域作為初始解m0。

3) 在迭代變量k=1,2,…,K時(shí),執(zhí)行下述迭代過(guò)程:

a) 更新模型加權(quán)矩陣Wm和數(shù)據(jù)加權(quán)矩陣Wd;

c) 共軛梯度法解方程Amk=y,得到第k次迭代解mk;

d) 修改平衡參數(shù)λ。

4) 得到Radon系數(shù)向量m=mk。

5) 對(duì)所有頻率的數(shù)據(jù)做前述相同的運(yùn)算,得到數(shù)據(jù)的高分辨率Radon變換系數(shù)。

3.2 受限模型空間Radon變換

前述經(jīng)典高分辨率Radon變換使用的是Radon變換的全部系數(shù),計(jì)算成本高[47]。LIU等[48]提出壓縮計(jì)算域,提高了計(jì)算效率,但每次迭代時(shí)更新算子,影響算法的收斂性。LU[27]采用迭代收縮閾值算法提高反演的收斂速度,但只有閾值可以控制系數(shù)的稀疏性,需要屢次嘗試才能選到合適的閾值。WANG等[49]運(yùn)用貪婪算法,依據(jù)每次迭代模型的殘差幅值選取子空間。SABBIONE等[47]直接在受限模型空間中運(yùn)用共軛梯度算法,保證了目標(biāo)函數(shù)不會(huì)隨迭代而改變,并通過(guò)試驗(yàn)證明了受限模型空間中的Radon變換系數(shù)足以表示地震數(shù)據(jù)。

受限模型空間在迭代前就被選取,它是指伴隨Radon變換系數(shù)高于閾值的部分構(gòu)成的空間,不隨迭代進(jìn)行而改變,這使得受限模型空間方法計(jì)算效率優(yōu)于貪婪最小二乘方法。而在結(jié)果分辨率方面,受限模型空間方法優(yōu)于迭代收縮方法,稀疏性不僅由閾值控制,而且受概率先驗(yàn)信息的約束;從另一個(gè)角度看,此閾值選取較為隨意[23]。受限模型空間定義[47]為:

(24)

式中:T是閾值,取值范圍0

(25)

方程(25)同樣采用IRCG算法求解。

4 Radon變換在地震矢量數(shù)據(jù)處理中的應(yīng)用

Radon變換已廣泛應(yīng)用于傳統(tǒng)縱波勘探地震數(shù)據(jù)處理的多次波壓制、地震數(shù)據(jù)保幅重建、去噪、鄰炮干擾壓制以及波場(chǎng)走時(shí)信息提取等。隨著多波地震勘探和微震監(jiān)測(cè)技術(shù)在非常規(guī)油氣藏勘探開發(fā)中的大量應(yīng)用,急需研發(fā)針對(duì)多分量地震數(shù)據(jù)的Radon變換方法,如何通過(guò)Radon變換獲得良好的處理效果并保持信號(hào)的矢量場(chǎng)特性是研究的重點(diǎn)。

4.1 縱橫波分離

PP波和PS波的分離是多波多分量地震勘探技術(shù)中一個(gè)重要的處理環(huán)節(jié)?；赑P波和PS波的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)特征,目前分離方法主要有:Radon變換法[50-52]、極化濾波法[24]、波動(dòng)方程法[53]、散度和旋度法[54]等。Radon變換法分離PP波和PS波直接利用了二者的視速度差異,基本思想為:地震剖面上相互交叉的PP波和PS波雙曲同相軸,在雙曲Radon域中可以稀疏表示,易于分離。但是,Radon變換只考慮了PP波和PS波的速度和走時(shí)信息,對(duì)于復(fù)雜介質(zhì),地震反射波的時(shí)距曲線不再是標(biāo)準(zhǔn)的雙曲線,即使采用高分辨率高階雙曲Radon變換也未必能將波場(chǎng)稀疏表示,此時(shí)無(wú)法有效分離波場(chǎng)。對(duì)此,WANG等[24]提出基于τ-p域下的空間矢量旋轉(zhuǎn)波場(chǎng)分離方法,利用了PP波和PS波的速度差異和偏振特性,波場(chǎng)分離精度更高。

4.2 微震信號(hào)提取與噪聲壓制

Radon變換在反射同相軸的識(shí)別、壓制隨機(jī)噪聲方面有良好的效果[14]。薛昭等[38]證明了高分辨率Radon變換去噪方法的相對(duì)保幅性。目前,針對(duì)多分量地震數(shù)據(jù)去噪問(wèn)題,普遍采用分量單獨(dú)處理,而聯(lián)合處理的技術(shù)應(yīng)用較少[55]。實(shí)際上,多分量地震數(shù)據(jù)是一個(gè)矢量波場(chǎng),接收器接收到的各個(gè)分量數(shù)據(jù)是位移矢量在不同坐標(biāo)方向的投影。若不考慮多分量地震數(shù)據(jù)的矢量特性,各個(gè)分量單獨(dú)做Radon變換去噪處理會(huì)導(dǎo)致弱分量信號(hào)被壓制。FORNASIER等[56]認(rèn)為各個(gè)分量有效信號(hào)的稀疏表示系數(shù)具有相同的稀疏分布模式,并提出組稀疏閾值準(zhǔn)則;VERA RODRIGUEZ等[57]證明了基于組稀疏表示的去噪方法可以有效地保持多分量地震數(shù)據(jù)分量之間的振幅相對(duì)關(guān)系。對(duì)于三分量微震信號(hào),由變換域同一位置的各分量稀疏系數(shù)組成向量g,其L2范數(shù)‖g‖2的大小決定系數(shù)是否為組稀疏。

與組稀疏的思想相近,三分量微震信號(hào)的均方根包絡(luò)[58]為:

(26)

式中:g={ex(x,t),ey(x,t),ez(x,t)},ex,ey,ez為3個(gè)分量的包絡(luò)。均方根包絡(luò)常被用于識(shí)別微震事件。為避免震源輻射花樣的影響[23],SABBIONE等[59]對(duì)微震信號(hào)的均方根包絡(luò)進(jìn)行頂偏Radon變換,識(shí)別微震事件。各個(gè)分量分別進(jìn)行高分辨率頂偏Radon變換,根據(jù)微震事件截取信號(hào)、壓制噪聲。為高效計(jì)算,可采用兩種措施:在時(shí)空窗中掃描事件以限制變換的道數(shù);在受限空間中計(jì)算Radon變換[59]。即便用于處理低信噪比的微震數(shù)據(jù),此方法依然能成功檢測(cè)微震事件并壓制噪聲。

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

Radon變換在地震學(xué)領(lǐng)域獲得了廣泛的應(yīng)用,并已在多次波壓制、縱橫波分離、地震數(shù)據(jù)保幅重建、去噪以及微震檢測(cè)等方面取得良好的應(yīng)用效果。本文總結(jié)了Radon變換目前的研究成果,著重討論了稀疏正則化、高效計(jì)算等方面的研究歷史與現(xiàn)狀;并簡(jiǎn)單介紹了Radon變換在地震矢量數(shù)據(jù)處理中的應(yīng)用。高分辨率Radon變換常用于地震波場(chǎng)分離;受限模型空間Radon變換在海量地震數(shù)據(jù)處理提高計(jì)算效率方面獲得了應(yīng)用嘗試。但這些方法改進(jìn)和應(yīng)用是以地震標(biāo)量波場(chǎng)單分量數(shù)據(jù)處理為主,未考慮到多分量地震信號(hào)的矢量特征,也無(wú)法保持多分量地震信號(hào)處理中的矢量特性不畸變,真正針對(duì)Radon變換在地震矢量場(chǎng)處理中的應(yīng)用研究很少。

本文認(rèn)為可將組稀疏約束的概念與高分辨率Radon變換相結(jié)合,發(fā)展適合于多分量矢量場(chǎng)的Radon變換處理方法。但基于組稀疏約束的Radon變換也只是利用了地震矢量場(chǎng)的振幅信息,未兼顧矢量場(chǎng)的方向性。因此,探索基于矢量場(chǎng)數(shù)據(jù)處理的Radon變換方法,利用矢量場(chǎng)的方向特性、偏振特性等改進(jìn)現(xiàn)有的Radon變換方法將是以后重點(diǎn)研究方向之一。此外,隨著Radon變換方法的發(fā)展,變換系數(shù)的維度也越來(lái)越高,運(yùn)算成本過(guò)高限制了Radon變換的實(shí)際應(yīng)用,因此,基于CPU或GPU并行加速的Radon變換方法也是未來(lái)的發(fā)展方向之一。

[1] 王有新.應(yīng)用地震數(shù)據(jù)處理方法[M].北京:石油工業(yè)出版社,2009:81-85

WANG Y X.Applied seismic data processing method[M].Beijing:Publishing House of Oil Industry,2009:81-85

[2] 劉法啟,張關(guān)泉.借助廣義Radon變換進(jìn)行方位角校正[J].石油物探,1996,35(4):43-51

LIU F Q,ZHANG G Q.Azimuth correction with the aid of the generalized Radon transform[J].Geophysical Prospecting for Petroleum,1996,35(4):43-51

[3] 賈春梅,姜國(guó)慶,劉志成,等.頻域稀疏雙曲Radon變換去噪方法[J].物探與化探,2016,40(3):527-533

JIA C M,JIANG G Q,LIU Z C,et al.Denoising method based on sparse hyperbolic Radon transform in the frequency domain[J].Geophysical and Geochemical Exploration,2016,40(3):527-533

[4] 鞏向博,韓立國(guó),王升超.混合域高分辨率雙曲Radon變換及其在多次波壓制中的應(yīng)用[J].石油物探,2016,55(5):711-718

GONG X B,HAN L G,WANG S C.High-resolution hyperbolic Radon transform and its application in multiple suppression[J].Geophysical Prospecting for Petroleum,2016,55(5):711-718

[5] HAMPSON D.Inverse velocity stacking for multiple elimination[J].Expanded Abstractes of 56thAnnual Internat SEG Mtg,1986:422-424

[6] TRAD D,ULRYCH T,SACCHI M.Accurate interpolation with high-resolution time-variant Radon transforms[J].Geophysics,2002,67(2):644-656

[7] 謝俊法,孫成禹,韓文功.迭代拋物Radon變換法分離一次波與多次波[J].石油地球物理勘探,2014,49(1):76-81

XIE J F,SUN C Y,HAN W G.Primary and multiple separation based on iterative parabolic Radon transform[J].Oil Geophysical Prospecting,2014,49(1):76-81

[8] TRAD D,ULRYCH T,SACCHI M.Latest views of the sparse Radon transform[J].Geophysics,2003,68(1):386-399

[9] KOSTOV C.Toeplitz structure in slant-stack inversionJ].Expanded Abstractes of 60thAnnual Internat SEG Mtg,1990:1618-1621

[10] ANDERSON J E.Parabolic and linear 2D tau-p transforms using the generalized radon transform:center for wave phenomena[R].CWP Report,1993:109-120

[11] SACCHI M D,ULRYCH T J.High-resolution velocity gathers and offset space reconstruction[J].Geophysics,1995,60(4):1169-1177

[12] SACCHI M D,PORSANI M.Fast high resolution parabolic Radon transform[J].Expanded Abstracts of 69thAnnual Internat SEG Mtg,1999:1477-1480

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

XUE Y R,TANG H H,CHEN X H.Seismic data reconstruction based on high order high resolution Radon transform[J].Oil Geophysical Prospecting,2014,49(1):95-100

[14] 徐彥凱,曹思遠(yuǎn),何元.雙曲Radon-ASVD方法壓制疊前地震數(shù)據(jù)隨機(jī)噪聲[J].石油地球物理勘探,2017,52(3):451-457

XU Y K,CAO S Y,HE Y.Prestack seismic random noise attenuation with a hyperbolic Radon-ASVD[J].Oil Geophysical Prospecting,2017,52(3):451-457

[15] 王立歆,李強(qiáng),姬小兵,等.用Radon變換法消除沙丘鳴震的應(yīng)用及效果分析[J].石油物探,2002,41(1):88-91

WANG L X,LI Q,JI X B,et al.Application of Radon transformation to elimination of dune ringing and its effect analysis[J].Geophysical Prospecting for Petroleum,2002,41(1):88-91

[16] 范景文,李振春,宋翔宇,等.各向異性高分辨率Radon變換壓制多次波[J].石油地球物理勘探,2016,51(4):665-669

FAN J W,LI Z C,SONG X Y,et al.Multiple attenuation with anisotropic high-resolution Radon transform[J].Oil Geophysical Prospecting,2016,51(4):665-669

[17] 謝玉洪,李列,劉玉金,等.基于共散射點(diǎn)道集的多次波壓制方法研究[J].石油地球物理勘探,2015,50(2):232-237

XIE Y H,LI L,LIU Y J,et al.Multiple elimination in common scattering gathers[J].Oil Geophysical Prospecting,2015,50(2):232-237

[18] 張紅梅,劉洪.基于稀疏離散τ-p變換的非均勻地震道重建[J].石油物探,2006,45(2):141-145

ZHANG H M,LIU H.Sparseness discreteτ-ptransform in irregular seismic trace reconstruction[J].Geophysical Prospecting for Petroleum,2006,45(2):141-145

[19] 王升超,韓立國(guó),鞏向博.基于各向異性Radon變換的疊前地震數(shù)據(jù)重建[J].石油物探,2016,55(6):808-815

WANG S C,HAN L G,GONG X B.Prestack seismic data reconstraction by anisotropic Radon transform[J].Geophysical Prospecting for Petroleum,2016,55(6):808-815

[20] COVA R,HENLEY D,INNANEN K.Addressing shear wave static corrections in the ray parameter domain:a non-stationary interferometric approach[J].Expanded Abstracts of 85thAnnual Internat SEG Mtg,2015:2129-2133

[21] 邵廣周,李慶春.聯(lián)合應(yīng)用τ-p變換法和相移法提取面波頻散曲線[J].石油地球物理勘探,2010,45(6):836-840

SHAO G Z,LI Q C.S Joint application ofτ-pand phase-shift stacking method to extract ground wave dispersion curve[J].Oil Geophysical Prospecting,2010,45(6):836-840

[22] 張威,韓立國(guó),李洪建.基于起伏海水表面的拖纜鬼波壓制方法[J].石油物探,2017,56(4):500-506

ZHANG W,HAN L G,LI H J.Deghosting method based on a variable sea surface for conventional streamer seismic data[J].Geophysical Prospecting for Petroleum,2017,56(4):500-506

[23] SABBIONE J I,VELIS D R,SACCHI M D.Microseismic data denoising via an apex-shifted hyperbolic Radon transform[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:2155-2161

[24] WANG Y,SINGH S C,BARTON P J.Separation of P-and SV-wavefields from multi-component seismic data in theτ-pdomain[J].Geophysical Journal International,2002,151(2):663-672

[25] IBRAHIM A,SACCHI M D.Simultaneous source separation using a robust Radon transform[J].Geophysics,2014,79(1):V1-V11

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

XIONG D,ZHAO W,ZHANG J F.Hybrid-domain parabolic Radon transform and its application to demultiple[J].Chinese Journal of Geophysics,2009,52(4):1068-1077

[27] LU W.An accelerated sparse time-invariant Radon transform in the mixed frequency-time domain based on iterative 2D model shrinkage[J].Geophysics,2013,78(4):V147-V155

[28] CARY P W.The simplest discrete Radon transform[J].Expanded Abstracts of 68thAnnual Internat SEG Mtg,1998:1999-2002

[29] THORSON J R.Velocity-stack and slant-stack stochastic inversion[J].Geophysics,1985,50(12):2727-2741

[30] HARGREAVES N,VERWEST B,WOMBELL R,et al.Multiple attenuation using an apexshifted radon transform[J].Expanded Abstracts of 73rdAnnual Internat SEG Mtg,2003:1929-1932

[31] TRAD D.Interpolation and multiple attenuation with migration operators[J].Geophysics,2003,68(6):2043-2054

[32] TRAD D,SILIQI R,POOLE G,et al.Fast and robust deblending using Apex Shifted Radon transform[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-5

[33] IBRAHIM A,SACCHI M D.Fast simultaneous seismic source separation using Stolt migration and demigration operators[J].Geophysics,2015,80(6):WD27-WD36

[34] GONG X B,YU C X,WANG Z H.Separation of prestack seismic diffractions using an improved sparse apex-shifted hyperbolic Radon transform[J].Exploration Geophysics,2016:https:∥doi.org/10.1071/EG16031

[35] MOLDOVEANU-CONSTANTINESCU C,SACCHI M D.Enhanced resolution in Radon domain using the shifted hyperbola equation[J].Expanded Abstracts of 75thAnnual Internat SEG Mtg,2005:2277-2280

[36] CASTLE R J.A theory of normal moveout[J].Geophysics,1994,59(6):983-999

[37] BICKEL S H.Focusing aspects of the hyperbolic Radon transform[J].Geophysics,2000,65(2):652-655

[38] 薛昭,董良國(guó),單聯(lián)瑜.Radon變換去噪方法的保幅性理論分析[J].石油地球物理勘探,2012,47(6):858-867

XIE Z,DONG L G,SHAN L Y.Amplitude preservation theoretical analysis of Radon transform de-noise method[J].Oil Geophysical Prospecting,2012,47(6):858-867

[39] 薛亞茹,陳小宏,馬繼濤.多方向正交多項(xiàng)式變換壓制多次波[J].地球物理學(xué)報(bào),2012,55(10):3450-3458

XUE Y R,CHEN X H,MA J T.Multiples attenuation based on directional orthogonal polynomial transform[J].Chinese Journal of Geophysics,2012,55(10):3450-3458

[40] SCHONEWILLE M,ZWARTJES P.High-resolution transforms and amplitude preservation[J].Expanded Abstracts of 72ndAnnual Internat SEG Mtg,2002:2066-2069

[41] 李晶晶,孫成禹,謝俊法,等.相對(duì)保幅的拋物線Radon變換法地震道重建[J].石油物探,2014,53(2):181-187

LI J J,SUN C Y,XIE J F,et al.Seismic trace reconstruction by relative amplitude-preserved parabolic Radon transform[J].Geophysical Prospecting for Petroleum,2014,53(2):181-187

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

TANG H H,MAO W J.Amplitude preserved seismic data reconstruction by 3D high-order parabolic Radon transform[J].Chinese Journal of Geophysics,2014,57(9):2918-2927

[43] 曾有良,樂(lè)友喜,單啟銅.基于高分辨率Radon變換的VSP波場(chǎng)分離方法[J].石油物探,2007,46(2):115-119

ZENG Y L,LE Y X,SHAN Q T.VSP wavefield separation based on high-resolution Radon transform[J].Geophysical Prospecting for Petroleum,2007,46(2):115-119

[44] SACCHI M D.Aperture compensated Radon and Fourier transforms[D].Vancouver:The University of British Columbia,1996

[45] SCALES J A,GERZTENKORN A,TREITEL S.FastIpsolution of large,sparse,linear systems:application to seismic travel time tomography[J].Journal of Computational Physics,1988,75(2):314-333

[46] GORODNITSKY I F,RAO B D.Sparse signal reconstruction from limited data using FOCUSS:a re-weighted minimum norm algorithm[J].IEEE Transactions on Signal Processing,1997,45(3):600-616

[47] SABBIONE J I,SACCHI M D.Restricted model domain time Radon transforms[J].Geophysics,2016,81(6):A17-A21

[48] LIU Y,SACCHI M D.De-multiple via a fast least squares hyperbolic radon transform[J].Expanded Abstracts of 72ndAnnual Internat SEG Mtg,2002:2182-2185

[49] WANG J F,NG M,PERZ M.Fast high-resolution Radon transforms by greedy least-squares method[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:3128-3132

[50] 劉財(cái),楊慶節(jié),耿美霞.多波多分量地震資料速度掃描法轉(zhuǎn)換波識(shí)別技術(shù)[J].石油物探,2013,52(4):383-387

LIU C,YANG Q J,GENG M X.P-SV wave identification technique by velocity scanning for multi-wave multi-component seismic data[J].Geophysical Prospecting for Petroleum,2013,52(4):383-387

[51] 馮晅,張先武,劉財(cái),等.帶有多道相關(guān)的拋物線Radon變換法分離P-P、P-SV波[J].地球物理學(xué)報(bào),2011,54(2):304-309

FENG X,ZHANG X W,LIU C,et al.Saparating P-P and P-SV wave by parabolic Radon transform with multiple coherence[J].Chinese Journal of Geophysics,2011,54(2):304-309

[52] 陶春輝,何樵登.地面多分量地震資料縱橫波分離方法[J].石油物探,1993,32(2):47-55

TAO C J,HE Q D.Separation of P- and S-waves in multicomponent surface seismic data[J].Geophysical Prospecting for Petroleum,1993,32(2):47-55

[53] 馬德堂,朱光明.彈性波波場(chǎng)P波和S波分解的數(shù)值模擬[J].石油地球物理勘探,2003,38(5):482-486

MA D T,ZHU G M.Numerical modeling of P-wave and S-wave separation in elastic wavefield[J].Oil Geophysical Prospecting,2003,38(5):482-486

[54] DELLINGER J,ETGEN J.Wave-field separation in two-dimensional anisotropic media[J].Geophysics,1990,55(7):914-919

[55] 尋超,汪超,王赟.多方向矢量中值濾波在多分量地震數(shù)據(jù)中的應(yīng)用[J].石油物探,2016,55(5):703-710

XUN C,WANG C,WANG Y.The application of multi-directional vector median filtering in multi-component seismic data[J].Geophysical Prospecting for Petroleum,2016,55(5):703-710

[56] FORNASIER M,RAUHUT H.Recovery algorithms for vector valued data with joint sparsity constraints[J].SIAM Journal on Numerical Analysis,2008,46(2):577-613

[57] VERA RODRIGUEZ I,BONAR D,SACCHI M.Microseismic data denoising using a 3C group sparsity constrained time-frequency transform[J].Geophysics,2012,77(2):21-29

[58] MICHAUD G,LEANEY S.Continuous microseismic mapping for real-time event detection and location[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008:1357-1361

[59] SABBIONE J I,SACCHI M D,VELIS D R.Radon transform-based microseismic event detection and signal-to-noise ratio enhancement[J].Journal of Applied Geophysics,2015,113(1):51-63

(編輯：朱文杰)

ResearchstatusoftheapplicationofRadontransforminseismicvectorfieldprocessing

QIU Xinming1,2,WANG Chao2,YUAN Yijun1,WANG Yun1

(1.SchoolofGeophysicsandInformationTechnology,ChinaUniversityofGeosciences,Beijing100083,China;2.StateKeyLaboratoryofOreDepositGeochemistry,InstituteofGeochemistry,ChineseAcademyofSciences,Guiyang550081,China)

Through the extensive investigation of Radon transform with its various improvement and application effects,current research results on Radon transform are summarized.On the basis of introducing the basic principles and applicability of various improved Radon transform methods in detail,this paper focuses on the introduction of using sparse regularization constraints in Radon transform to improve the resolution,and of deriving solution in the restricted model space to improve the iteration efficiency for massive seismic data.Moreover,this paper demonstrates that,the application of Radon transform in seismic vector field processing is still in the experimental stage,and that its development trend is to obtain satisfactory processing effects as well as to keep the vector characteristics.Moreover,there is an urgent need to design and improve Radon transform methods for multi-component seismic data processing in the immediate future.

Radon transform,resolution,vector field,multi-component,sparse regularization

2017-04-01;改回日期2017-06-05。

邱新明(1994—),男,博士在讀,主要從事多分量地震數(shù)據(jù)處理方法研究。

汪超(1982—),男,博士,副研究員,主要從事地球物理矢量場(chǎng)數(shù)字信號(hào)處理研究。

國(guó)家自然科學(xué)基金項(xiàng)目(41425017,41504107,41574126和41604119)和國(guó)家科技重大專項(xiàng)“薄層地震波場(chǎng)特征與反演研究”(2016ZX05002-005-008)聯(lián)合資助。

This research is financially supported by the National Natural Science Foundation of China (Grant No.41425017,41504107,41574126,41604119) and the National Major Scientific and Technological Special Project (Grant No.2016ZX05002005-008).

P631

A

1000-1441(2017)06-0905-10

10.3969/j.issn.1000-1441.2017.06.017