羅騰騰 徐基祥 孫夕平

(中國石油勘探開發(fā)研究院，北京 100083)

0 引言

地層中的斷層、裂隙、尖滅和孔洞等小尺度地質(zhì)體對油氣勘探具有重要意義，這些小尺度地質(zhì)體在地震記錄上的響應(yīng)主要為能量較弱的繞射波。Klem等[1]指出地震數(shù)據(jù)中的繞射能量通常比反射能量弱1～2個數(shù)量級。在常規(guī)的地震數(shù)據(jù)處理流程中，來自反射界面以外的繞射波同相軸常被視作干擾噪聲而被濾除。即使同時將繞射波和反射波準(zhǔn)確偏移歸位，弱繞射能量仍會被大尺度地質(zhì)體反射能量淹沒。為了充分利用地震數(shù)據(jù)中的繞射信息，可從全波場中提取繞射波，實現(xiàn)繞射波單獨成像，進而提高繞射目標(biāo)的成像分辨率。

早在二十世紀(jì)五十年代，Krey[2]就利用繞射波信息研究斷層、斷點等，但直到二十世紀(jì)八十年代，學(xué)者們才真正意識到繞射波分離成像的重要性，并相繼提出不同的繞射波場分離方法。根據(jù)地震數(shù)據(jù)處理階段的不同，可以將繞射波分離與成像方法歸納為以下三類：

(1)基于疊前數(shù)據(jù)的繞射波分離與成像。該類方法是利用疊前道集中繞射波與反射波的運動學(xué)特征差異，從全波場中提取繞射波，進而實現(xiàn)繞射波單獨成像。Landa等[3]提出在共炮檢距道集上利用相位校正和定性分析檢測小尺度繞射體。Khaidu-kov等[4]在疊前共炮點道集中利用反射波與繞射波的運動學(xué)特征差異，即反射波聚焦于震源點的鏡像位置而繞射波仍保持發(fā)散的狀態(tài)，通過反射波聚焦、切除、再反聚焦方法成功分離出繞射波場。Nowak等[5]根據(jù)共炮點道集中反射波與繞射波同相軸橫向位置差異，借助于加權(quán)Radon變換分離反射繞射波場。Bansal等[6]討論了多種利用疊前數(shù)據(jù)壓制反射信號、增強繞射能量的方法，主要包括共炮檢距傾角濾波、正常時差校正傾角濾波、特征向量濾波、Radon濾波等，并且對比、分析了不同方法的優(yōu)缺點，認為特征向量法是較為有效的繞射波場分離方法。Taner等[7]利用平面波分解(Plan Wave Decomposition，PWD)濾波器生成平面波記錄，依據(jù)繞射波與反射波在平面波記錄中的同相軸差異，在疊前道集中分離繞射波場。黃建平等[8]對PWD解構(gòu)濾波器的使用方法和在疊前、疊后道集上分離成像繞射波的流程進行了系統(tǒng)總結(jié)?？籽┑萚9]和劉玉金等[10]證實了平面波記錄具有區(qū)分繞射波與反射波的特性，并利用平面波解構(gòu)濾波技術(shù)分離獲得繞射波場。朱生旺等[11]解決了PWD濾波器中繞射波場的低傾角信息估計不足問題，提高了繞射波成像結(jié)果的橫向分辨率。

(2)在疊加過程中進行繞射波分離與成像。Kanasewich等[12]直接對地震信號沿著繞射波走時曲線進行疊加，獲得了共炮檢距道集和共斷層點道集上的繞射波時間剖面。Berkovich等[13]、Dell等[14]分別利用繞射波多聚焦和共反射面的繞射走時計算方法，將選取的繞射波同相軸進行最優(yōu)疊加，得到包含繞射波的疊加剖面。Asgedom等[15]利用修正的共反射面(Common Reflection Surface，CRS)技術(shù)討論了相似度和多信號分類兩種方法在分離反射波與繞射波時的優(yōu)勢。

(3)在偏移過程中進行繞射波分離與成像。Landa等[16]推導(dǎo)了深度偏移傾角域共成像點道集(Common Imaging Gathers，CIGs)中反射波和繞射波的形態(tài)，并分析了偏移速度對兩者的影響。Zhu等[17]在小波束偏移過程中，利用角度域局部成像矩陣中反射與繞射能量的分布差異，即反射能量沿著傾角方向呈線性分布而繞射能量彌散分布在整個矩陣中，分離出繞射能量并對其單獨成像。Klokov等[18]將相似度掃描反射頂點去除與混合Radon變換方法相結(jié)合，在傾角域CIGs上進行繞射波場分離。Decker 等[19]在偏移傾角道集中，對每個固定傾角的成像剖面應(yīng)用PWD解構(gòu)濾波技術(shù)，獲得單獨的繞射波成像剖面。Zhang等[20]提出在單炮的反射角域CIGs中通過切除菲涅耳帶提取繞射波場。黃建平等[21]和李振春等[22]將逆時偏移與最小二乘偏移的優(yōu)勢相結(jié)合，分別發(fā)展了適用于復(fù)雜構(gòu)造的保幅成像處理方法。劉夢麗等[23]提出一種滿足能量守恒的逆散射成像條件，并將其與最小二乘逆時偏移(Reverse Time Migration，RTM)相結(jié)合，實現(xiàn)了角度濾波成像。汪天池等[24]在傾角域中利用逆時偏移提取傾角道集，然后采用中值濾波方法實現(xiàn)反射波與繞射波的分離。

近些年來，大量學(xué)者嘗試在傾角域CIGs上進行繞射波與反射波的分離。根據(jù)傾角域CIGs中反射波響應(yīng)曲線呈雙曲規(guī)律而繞射波呈現(xiàn)擬線性規(guī)律的差異，發(fā)展了平面波解構(gòu)濾波、掃描相似度頂點去除與混合Radon變換相結(jié)合[18]、相似譜能量分析等繞射波分離方法。同時，由于在傾角域CIGs中繞射波同相軸對速度誤差的敏感性，Landa等[16]、Reshef等[25]、Li等[26]將其用于偏移速度分析以提高速度提取精度。

鑒于傳統(tǒng)混合域Radon變換采用迭代重加權(quán)算法進行求解，每一步迭代常常涉及到大規(guī)模矩陣的求逆問題，計算量巨大。為此，本文引入迭代收縮閾值算法(Iterative Shrinkage Thresholding Algorithm，ISTA)，發(fā)展了一種基于迭代收縮高分辨率Radon變換的繞射波分離與成像方法。首先采用迭代收縮閾值算法，通過時間域簡單的模型收縮步驟獲得Radon域模型的稀疏性，以解決常規(guī)Radon變換中反射波與繞射波在Radon域中能量團聚焦性差、繞射波場分離效果不理想的問題； 然后在編程實現(xiàn)算法的基礎(chǔ)上，通過模型數(shù)據(jù)驗證了該方法的有效性和適用性。

1 方法原理

1.1 傾角域共成像點道集波場特征

常速介質(zhì)疊后偏移(圖1)反射面的深度z(x)可以表示為[17]

z(x)=z1+xtanα0

(1)

式中：z1為傾斜反射界面與z軸的交點深度；α0為反射界面的傾角；x為橫坐標(biāo)。在地表y處的自激自收地震響應(yīng)可表示為

(2)

式中：v是介質(zhì)速度；w為旅行時間。

由模型坐標(biāo)(x，z)和數(shù)據(jù)坐標(biāo)(y，w)之間的映射關(guān)系可得偏移公式，即

(3)

τ=wcosα

(4)

圖1 零炮檢距反射示意圖

式中：vm表示偏移速度；α為偏移傾角；τ表示截距時間。將式(2)代入式(3)和式(4)中，消去y可得傾斜反射界面在傾角域的成像表達式，即

(5)

式中γ=vm/v。根據(jù)式(1)所定義的橫坐標(biāo)xm處的反射面深度為z0，可得

xmsinα0+z1cosα0=z0cosα0

(6)

將式(6)代入式(5)，可得傾角域CIGs中固定橫向位置的反射波響應(yīng)表達式，即

(7)

式中τ0=2z0/v為反射點處的雙程垂直旅行時。

將式(7)對α求偏導(dǎo)可得到

(8)

從式(8)可以看出，當(dāng)偏移速度正確，即γ=1時，反射響應(yīng)為具有穩(wěn)相頂點的“笑臉”狀曲線，且穩(wěn)相頂點位于α=α0處。

類似地，可以得到地下繞射點的傾角域響應(yīng)表達式

(9)

式中ξ=(x-x0)/z0，ξ是繞射點x0和成像點x之間的橫向距離與繞射點深度z0的比值。當(dāng)偏移速度準(zhǔn)確(vm=v)且成像點在繞射點處(ξ=0)，式(9)簡化為

τ(α)=τ0

(10)

由式(10)可以看出，此時繞射響應(yīng)為一水平直線。

對于式(7)和式(9)，利用τ(α)=2z(α)/vm和τ0=2z0/v，即可從時間域轉(zhuǎn)換到偽深度域。深度域反射界面在傾角域的響應(yīng)表達式為

(11)

地下繞射點在傾角域的響應(yīng)表達式為

(12)

式中ρ=x-x0。

通過一個簡單的常速介質(zhì)模型(圖2)可直觀地闡明傾角域CIGs響應(yīng)曲線的特點，并且可以進一步分析兩者之間的差異。該模型由一個垂直斷層和一個繞射體組成，記錄共150炮，最大炮檢距為500m，時間采樣間隔為1ms，檢波點間距為10m。

圖2 常速介質(zhì)模型垂直斷層和繞射體分別位于0.8km和1.2km

圖3和圖4分別顯示了0.8km處斷層和1.2km處繞射體在不同偏移速度下產(chǎn)生的傾角道集。可以看出，在傾角域CIGs上，位于繞射點正上方的繞射響應(yīng)可分為三種情況：當(dāng)偏移速度小于準(zhǔn)確的介質(zhì)速度時，繞射波表現(xiàn)為向上彎曲的“笑臉”狀(凹形)(圖3a、圖4a)；當(dāng)偏移速度準(zhǔn)確時，繞射波同相軸表現(xiàn)為線性形態(tài)(圖3b、圖4b)；當(dāng)偏移速度大于準(zhǔn)確的介質(zhì)速度時，繞射波同相軸則為向下彎曲的“哭臉”狀(凸形)(圖3c、圖4c)。而反射界面響應(yīng)曲線無論偏移速度準(zhǔn)確與否，始終保持為具有穩(wěn)相頂點的“笑臉”狀，且穩(wěn)相頂點橫坐標(biāo)指示界面的傾角。

綜上所述，在傾角域CIGs中，反射響應(yīng)與繞射響應(yīng)曲線之間存在著明顯的幾何形態(tài)差異，因而可以利用兩者之間形態(tài)差異作為波場分離的依據(jù)。

圖3 不同偏移速度誤差時斷層處的傾角域CIGs(a)偏移速度為介質(zhì)速度的90%； (b)偏移速度等于介質(zhì)速度； (c)偏移速度為介質(zhì)速度的110%藍色和紅色箭頭分別指示反射波和繞射波同相軸，圖4同

圖4 不同偏移速度誤差時繞射體處的傾角域CIGs(a)偏移速度為介質(zhì)速度的90%； (b)偏移速度等于介質(zhì)速度； (c)偏移速度為介質(zhì)速度的110%

1.2 常規(guī)Radon變換

二維拋物線Radon變換的定義式可以表示為

(13)

式中：d(t，u)為時間—空間域地震數(shù)據(jù)，其中t為時間，u為炮檢距；m(τ，q)為Radon域數(shù)據(jù)，其中q表示曲率參數(shù)。

因此，拋物線Radon正變換的離散形式可以表示為

(14)

式中：D(ω，u)和M(ω，q)分別為d(t，u)和m(τ，q)的傅里葉變換結(jié)果，其中ω為頻率；n為Radon域的曲率個數(shù)。

由式(14)可知，Radon變換可以表示為一個線性方程組的求解問題，即

d=Lm

(15)

式中：d為地震數(shù)據(jù)；m為模型數(shù)據(jù)；L是Radon變換算子，由地震數(shù)據(jù)的采集參數(shù)和Radon變換參數(shù)共同確定。

求式(15)的反問題，不能完全滿足存在性、唯一性和穩(wěn)定性，故該反問題是不適定的。對于不適定問題，為獲得唯一且穩(wěn)定的解，最直接方法就是在反問題中引入與模型信息相關(guān)的正則化項以求解稀疏正則解，即高分辨率Radon模型解。

對于式(15)，可構(gòu)造如下目標(biāo)函數(shù)，即

(16)

(17)

對式(17)求極小值，可以得到m的常規(guī)最小二乘解，即

m=(LTL+βI)-1LTd

(18)

式中I為單位矩陣。由式(18)可以看出，在傳統(tǒng)最小二乘意義下求得的Radon模型解分辨率較低。由于對每一個頻率所采用的阻尼因子相同，且阻尼因子的實際作用是對Radon模型數(shù)據(jù)進行平滑，導(dǎo)致模型空間的能量團無法較好地聚焦，降低了Radon變換的分辨率。因此，常規(guī)最小二乘Radon變換不能滿足繞射波場分離的高分辨率需求。

1.3 迭代收縮高分辨率Radon變換

為了獲得稀疏的高分辨率Radon模型解，一般為正則項選用L1范數(shù)、數(shù)據(jù)誤差項選用L2范數(shù)，即p=1、s=2，則式(18)可轉(zhuǎn)換為

(19)

在求解稀疏反演問題(式(19))時，本文利用Lu[27]提出的混合頻率—時間域迭代收縮閾值算法，即在時間域通過模型收縮步驟獲得Radon模型的稀疏性，矩陣與向量的乘積操作在頻率域?qū)崿F(xiàn)。具體來說，在第k步迭代過程中，更新的模型解為

mk=Tσ{mk-1+2t0F-1{(LTL)-1LT[F(d)-

LF(mk-1)]}

(20)

式中：k為當(dāng)前迭代次數(shù)；t0是迭代步長；F、F-1分別表示正、反傅里葉變換；Tσ是一個收縮算子，定義為

(21)

(22)

本文基于迭代收縮算法的高分率Radon變換處理步驟如下。

(1)給定閾值σ、最大迭代次數(shù)K和二維中值濾波器的大小以及初始迭代步長t1，將當(dāng)前迭代次數(shù)k設(shè)置為0。

(2)選取Radon模型的最小二乘解作為模型迭代的初始值，m0=F-1[(LTL)-1LTd]，同時設(shè)置k=1。

(3)在迭代變量k=1，2，...，K時，執(zhí)行下述迭代過程：

①將上次的模型值mk-1傅里葉變換到頻率域得到F(mk-1)；

②計算Radon算子L與頻率域模型值的乘積即LF(mk-1)；

③計算頻率域中數(shù)據(jù)擬合誤差F(d)-LF(mk-1)；

④計算算子的偽逆矩陣(LTL)-1LT與數(shù)據(jù)擬合誤差的乘積，即(LTL)-1LT[F(d)-LF(mk-1)]；

⑤利用傅里葉反變換將頻率域數(shù)據(jù)變換到時間域，即F-1{(LTL)-1LT[F(d)-LF(mk-1)]}；

⑥利用帶有收縮運算的方程Tσ計算模型更新值mk。

(4)循環(huán)增加迭代次數(shù)k，檢查是否達到了最大迭代次數(shù)，如果“否”，則轉(zhuǎn)到第(3)步。

(5)輸出Radon模型的高分辨率解。

利用上述迭代收縮高分辨Radon變換算法，在求解過程中僅需計算一次偽逆矩陣，相比于傳統(tǒng)的迭代重加權(quán)最小二乘(Iterative Reweighted Least Squares，IRLS)求解算法，極大地降低了計算量，更適用于實際地震資料處理。

2 模型試算

2.1 有效性

下面用模型試算驗證本文算法的有效性。模擬CIGs中包含三條反射波和兩條繞射波同相軸(圖5)，反射波頂點分別位于深度600、1000、1500m處，兩條繞射波同相軸則分別處于深度1300、2000m的位置。

圖6為模型數(shù)據(jù)在不同迭代次數(shù)下的試算結(jié)果。測試過程中，設(shè)置最大迭代次數(shù)K=20、步長t0=0.5、閾值算子σ=0.05，二維均值濾波器的大小與模型相同。

圖5 模擬數(shù)據(jù)的CIG記錄

從圖6可以看出，隨著迭代次數(shù)的增加，Radon域模型的分辨率逐漸提高，即反射波和繞射波能量團聚焦性逐漸變好，剪刀狀拖尾能量逐漸變?nèi)?，分離后的繞射波場也越來越干凈。在此模擬數(shù)據(jù)的測試過程中，當(dāng)?shù)螖?shù)K=20時，Radon域中的反射波和繞射波能量聚焦性很好，分辨率高，切除反射后的繞射波場中幾乎不存在殘余反射能量。

2.2 抗噪能力

為了測試本文方法的抗噪能力，在圖5的模擬數(shù)據(jù)中加入了高斯白噪，得到如圖7所示數(shù)據(jù)。圖8左為圖7的常規(guī)最小二乘Radon變換結(jié)果，包含大量隨機噪聲，圖8右為分離后得到的時空域繞射波波場，由于在Radon域?qū)υ肼暤膲褐撇蛔?，?yán)重影響分離后的繞射波場質(zhì)量。圖9為迭代收縮高分辨率Radon變換結(jié)果?？梢钥闯觯旱湛s高分辨率Radon變換方法能夠壓制原始數(shù)據(jù)中的隨機噪聲，有效提高原始數(shù)據(jù)的信噪比(圖9左)；本文方法從低信噪比的數(shù)據(jù)中提取繞射波的能力較強，可很好地壓制模擬數(shù)據(jù)中的隨機噪聲(圖9右)。

圖6 不同迭代次數(shù)下迭代收縮高分辨率Radon變換結(jié)果

圖8 常規(guī)最小二乘Radon變換結(jié)果左為Radon域全波場，右為時空域繞射波波場

圖9 迭代收縮高分辨率Radon變換結(jié)果左為Radon域全波場，右為時空域繞射波波場

2.3 計算效率

為了說明本文方法的計算效率相較于IRLS算法的優(yōu)勢，不同迭代次數(shù)下本文方法和IRLS算法在模擬數(shù)據(jù)上的運算時間如表1所示。從表中可以看出，IRLS算法的運算時間與迭代次數(shù)具有緊密的聯(lián)系。當(dāng)?shù)?次時，本文方法運算時間比IRLS算法略大，這是由于基于迭代收縮算法的高分辨率Radon變換需要在迭代前對變換算子L進行一次廣義逆(LTL)-1LT的計算，并將其儲存于內(nèi)存中，使后續(xù)迭代只需要計算矩陣向量相乘以及向量之間的運算即可，這一特性也是本文方法能夠?qū)崿F(xiàn)快速運算的核心所在。隨著迭代次數(shù)的增加，本文方法運算所需時間變化平緩，而IRLS算法的運算時間隨迭代次數(shù)的增加而陡然變化。通過對比可以看出，針對此模擬數(shù)據(jù)，當(dāng)?shù)螖?shù)超過20次后，本文方法的計算效率是常規(guī)IRLS算法的兩倍左右。

表1 兩種Radon變換方法運算時間 s

3 實際資料應(yīng)用

選取中國西部M探區(qū)實際地震資料測試本文方法。M探區(qū)溶蝕孔洞和裂縫儲層十分發(fā)育，利用傳統(tǒng)的成像處理流程難以對儲層內(nèi)部的多種繞射目標(biāo)體精確成像。

首先，利用本文方法對某一成像點下Kirchhoff疊前深度偏移(PSDM)產(chǎn)生的傾角域CIGs進行反射波與繞射波的分離測試，如圖10所示。

圖10a和10b分別為分離前實際地震資料的傾角域CIGs和經(jīng)過本文方法分離后的繞射波CIGs。從圖中可以看出，全波場CIGs中的強反射能量被明顯壓制，盡管在分離后的繞射波CIGs中(圖10b)無法清晰地觀察到分離出的弱繞射能量，但是強能量的反射同相軸被明顯削弱，說明了本文方法分離反射波與繞射波是有效的。

圖10 基于迭代高收縮分辨率Radon變換的CIGs波場分離(a)實際地震資料的全波場CIGs； (b)波場分離后的繞射波CIGs

進一步將本文方法應(yīng)用于該地區(qū)Inline1320測線的傾角域CIGs，圖11顯示了測線Inline1320的0～60°傾角域CIGs。由圖可見，強反射波能量呈雙曲線形態(tài)，繞射能量則呈水平直線或者曲率較小的擬線性同相軸。

圖12為經(jīng)過Radon變換壓制反射波后分離得到的繞射波傾角域CIGs。由圖可見，經(jīng)過本文方法處理后，大部分強反射能量同相軸被去除，而同一深度點的繞射波被較好地保留。

該測線波場分離前的PSDM剖面如圖13a所示，利用本文方法分離的繞射波成像剖面如圖13b所示。在原始疊加剖面(圖13a)中，由于強反射能量的干擾，繞射體的位置難以識別。從圖13b中可以看出，上覆地層反射能量殘余多，中深層反射殘余較弱(幾乎全部去除)，主要原因在于上覆地層采集腳印明顯，前期CDP道集預(yù)測去除反射能量不多，而下伏地層反射波能量得到準(zhǔn)確壓制。相較于原始偏移剖面，經(jīng)過繞射分離后得到的繞射波成像剖面(圖13b)中，強反射同相軸被濾除，淹沒的繞射體位置凸顯出來(圖13b中黃色圓圈所示)。因此，經(jīng)過迭代收縮高分辨率Radon變換方法的處理，繞射體的成像結(jié)果被凸顯出來，說明本文方法適用于地下小尺度非均質(zhì)體的識別和精確定位。

圖11 Inline1320測線傾角域CIGs(0～60°)

圖12 Inline1320測線Radon變換后的傾角域CIGs(0～60°)

圖13 中國西部M探區(qū)實際資料(Inline1320測線)處理結(jié)果(a)波場分離前PSDM剖面； (b)波場分離后繞射波成像剖面

4 結(jié)論

本文基于迭代收縮高分辨率Radon變換的繞射波分離成像方法，利用反射波與繞射波在傾角域CIGs中的的運動學(xué)特征差異，并通過構(gòu)建一個常速介質(zhì)模型以說明兩者之間的差異，進一步利用反射波與繞射波的形態(tài)差異在傾角域CIGs中分離繞射波場。根據(jù)模擬數(shù)據(jù)和實際地震資料處理結(jié)果得出以下結(jié)論：

(1)在傾角域CIGs中，反射波同相軸總是呈現(xiàn)開口向上的“笑臉”狀曲線，繞射同相軸則表現(xiàn)為擬線性，根據(jù)兩者曲線形態(tài)的差異可以實現(xiàn)波場分離；

(2)迭代收縮高分辨率Radon變換對反射和繞射波同相軸在Radon域的聚焦能力強，分辨率高，且在計算過程中采用收縮閾值算法，不僅提高了反演的收斂速度，而且較常規(guī)方法具有更好的抗噪性；

(3)本文方法雖然能夠有效分離反射波與繞射波，并利用繞射信息揭示了地下小規(guī)模繞射體，可提高地震資料解釋的精度，但由于算法中需要較多的人工干預(yù)(如迭代步長、閾值算子等的選取)，存在一定的局限性。