應用迭代收縮高分辨率Radon變換的繞射波分離與成像方法

羅騰騰 徐基祥 孫夕平

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

0 引言

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

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

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

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

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

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

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

1 方法原理

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

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

z(x)=z1+xtanα0

(1)

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

(2)

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

由模型坐標(x，z)和數據坐標(y，w)之間的映射關系可得偏移公式，即

(3)

τ=wcosα

(4)

圖1 零炮檢距反射示意圖

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

(5)

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

xmsinα0+z1cosα0=z0cosα0

(6)

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

(7)

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

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

(8)

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

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

(9)

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

τ(α)=τ0

(10)

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

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

(11)

地下繞射點在傾角域的響應表達式為

(12)

式中ρ=x-x0。

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

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

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

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

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

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

1.2 常規Radon變換

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

(13)

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

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

(14)

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

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

d=Lm

(15)

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

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

對于式(15)，可構造如下目標函數，即

(16)

(17)

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

m=(LTL+βI)-1LTd

(18)

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

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

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

(19)

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

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

LF(mk-1)]}

(20)

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

(21)

(22)

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

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

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

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

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

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

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

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

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

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

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

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

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

2 模型試算

2.1 有效性

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

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

圖5 模擬數據的CIG記錄

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

2.2 抗噪能力

為了測試本文方法的抗噪能力，在圖5的模擬數據中加入了高斯白噪，得到如圖7所示數據。圖8左為圖7的常規最小二乘Radon變換結果，包含大量隨機噪聲，圖8右為分離后得到的時空域繞射波波場，由于在Radon域對噪聲的壓制不足，嚴重影響分離后的繞射波場質量。圖9為迭代收縮高分辨率Radon變換結果。可以看出：迭代收縮高分辨率Radon變換方法能夠壓制原始數據中的隨機噪聲，有效提高原始數據的信噪比(圖9左)；本文方法從低信噪比的數據中提取繞射波的能力較強，可很好地壓制模擬數據中的隨機噪聲(圖9右)。

圖6 不同迭代次數下迭代收縮高分辨率Radon變換結果

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

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

2.3 計算效率

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

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

3 實際資料應用

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

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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