朱生旺，李 佩，寧俊瑞

中國石油化工股份有限公司石油勘探開發(fā)研究院，北京 100083

1 引 言

在反射地震數(shù)據(jù)中，繞射波是反映地下介質(zhì)非連續(xù)性的重要信息，對于斷裂和地質(zhì)異常體的解釋，特別是針對碳酸鹽巖縫洞儲層這類復(fù)雜非均質(zhì)儲層的認(rèn)識，要求地震資料處理中保護繞射波信息并使其精確成像是十分重要的．

如何有效利用繞射波信息來反映地下介質(zhì)的橫向變化，前人已做過很多研究和嘗試，如：Landa等［1］根據(jù)繞射波運動學(xué)特性，提出通過相關(guān)來增強繞射點位置的地震信號振幅，以探測局部非均質(zhì)異常體的方法；Pozdnyakov［2］提出的估計散射強度屬性的聚焦變換法；Khaidukov［3］提出的聚焦－去焦（focusing－defocusing）繞射成像方法，等等．在實際中反映局部地質(zhì)異常或儲層的繞射波信息通常淹沒在強反射信息中，因此人們又提出直接從疊前數(shù)據(jù)中分離出繞射波并單獨對其進行偏移成像，目的是使成像結(jié)果能夠凸顯非均質(zhì)地質(zhì)體的邊界和形體，從而幫助更準(zhǔn)確地預(yù)測與認(rèn)識儲層．如Nowak等［4］嘗試了從地震記錄中分離繞射波的方法，即利用炮集記錄中反射波與繞射波同相軸軌跡的幾何差別，通過雙曲線Radon變換，將繞射波從地震記錄中分離出來．

實際上，繞射與反射沒有明確的界限，因此很難從反射記錄中嚴(yán)格地將繞射波分離出來．進行繞射波分離只能是加強源于橫向非連續(xù)性的反射信息，而壓制連續(xù)反射界面的反射信號．提取繞射波信息的方法基本上是傳統(tǒng)的信號分離方法，假設(shè)反射界面傾角不大，這樣在一定的道集（如共偏移距道集）或某些變換結(jié)果（如平面波分解結(jié)果）上繞射波與反射波同相軸存在較大的傾角差異，依此實現(xiàn)反射與繞射主要信息的分離．Reeshidev［5］總結(jié)了增強繞射波的幾種方法，包括分別在炮點和檢波點域壓制正向散射能量的雙濾波、共偏移距道集上的傾角濾波、特征向量濾波、Harlan信噪分離［6］等方法，這些方法有各自的優(yōu)勢和局限性；Sergey［7］給出了平面波去反射濾波，也就是利用傾角差異壓制反射波的局部傾角濾波方法，該方法的最大優(yōu)點是局部適應(yīng)性好，但要求輸入數(shù)據(jù)要有較高信噪比，因此適合在平面波分解結(jié)果上應(yīng)用；M．Turhan Taner［8］利用Sergey的方法給出了繞射波分離的理論算例．

分離出的繞射波信息的完整性直接影響繞射波成像結(jié)果的質(zhì)量．即使假設(shè)反射波是來自接近水平的低傾角界面，我們?nèi)圆荒芡ㄟ^單一的傾角濾波方法分離出相對完整的繞射波信息，這是由于在低傾角附近，繞射波與反射波不能由同相軸傾角區(qū)分．因此，基于傾角差異的分解方法都普遍存在低傾角繞射波信息失真或缺失問題．低傾角信息缺損使空間有效帶寬變窄，這將導(dǎo)致繞射波成像結(jié)果橫向分辨率降低．為了彌補繞射波分離結(jié)果的帶寬損失，我們提出局部傾角濾波和頻率－空間域預(yù)測反演聯(lián)合的繞射波分離方法，其實現(xiàn)步驟是：（1）對單炮記錄進行平面波分解，并經(jīng)過數(shù)據(jù)重排得到對應(yīng)不同斜率參數(shù)的剖面；（2）對每一個共斜率剖面用平面波濾波方法實現(xiàn)局部傾角濾波，提取繞射波高傾角信息成分，并對剔除高傾角成分后的數(shù)據(jù)通過預(yù)測反演得到繞射波的低傾角信息成分；（3）將提取的繞射波高、低傾角信息相加得到最終繞射波估計；（4）進行平面波合成和數(shù)據(jù)重排得到繞射波炮集數(shù)據(jù)．理論試算結(jié)果表明，這種繞射波分離方法能比較有效地克服靠單一的傾角差異進行繞射波分離時損失低傾角信息這一問題，所得到的繞射信息相對完整，從而使繞射波成像結(jié)果有較高的橫向分辨率．

在塔河碳酸鹽巖縫洞發(fā)育區(qū)，孔洞結(jié)構(gòu)與分布復(fù)雜，孔洞尺度不一、填充物多樣性、幾何形狀不規(guī)則、空間變化劇烈，其大小、充填、形態(tài)同繞射波特征有著密切的聯(lián)系［9］，將地震波場中的繞射波分離出來單獨成像，有利于提高縫洞的識別精度，將本文提出的繞射波分離方法應(yīng)用在塔河地區(qū)，壓制了風(fēng)化殼的強反射能量，突出了縫洞的串珠信息，其成像結(jié)果能夠更好地識別縫洞．

2 方法原理

繞射波與反射波的一個顯著差異是時距曲線的不同，地面地震記錄中的繞射時距曲線由雙平方根方程描述，而反射時距曲線為雙曲線．這種時距關(guān)系上的差異反映在平面波域是，在平面波分解（τ－p）域的共斜率（p）剖面上，繞射波同相軸能量大部分包含在高傾角信息成分中［8］．因此，假設(shè)反射界面傾角較小，則在τ－p域的共p剖面上，反射波與繞射波同相軸存在較大的傾角差異，二者具有較好的可分離性．

將地震炮記錄變換到τ－p域，用sp（x，τ）表示斜率為p的共p剖面，這里x為炮點地面坐標(biāo)．設(shè)sp（x，τ）中反射波同相軸的時間傾角（相鄰道時差）在（－σm，σm）范圍內(nèi)，則通過局部傾角濾波可將sp（x，τ）分解為兩部分，即

最后，由

給出共p剖面上的繞射波分離結(jié)果．

現(xiàn)用理論模型數(shù)據(jù)來說明上述方法的實現(xiàn)過程．在圖1中，（a）為生成基于聲波方程的數(shù)值模擬數(shù)據(jù)的速度模型；（b）為平面波分解的共傾角剖面；（c）為從（b）中分離出的繞射波高傾角成分，即（x，τ）；（d）為（b）與（c）相減的結(jié)果，即（x，τ）；（e）為從（d）中分離出的繞射波低傾角成分，即（x，τ）；（f）為（d）與（e）相減的結(jié)果，即反射波估計；（g）為（c）與（e）相加得到的最終繞射波分離結(jié)果dp（x，τ）；（h）為原始數(shù)據(jù)的深度偏移剖面；（i）為分離出的繞射波的深度偏移剖面．從該模型數(shù)據(jù)處理結(jié)果可以看到，分離出的繞射波信息相對完整，與原始數(shù)據(jù)成像剖面比較，在分離出的繞射波成像剖面上，連續(xù)界面的反射波得到壓制后，洞的反射變得清晰．

圖1 繞射波分離過程與結(jié)果（a）速度模型；（b）平面波分解的共傾角剖面；（c）用傾角濾波方法從（b）中分離出的繞射波高傾角成分；（d）（b）與（c）相減的結(jié)果；（e）用預(yù)測反演方法從（d）中分離出的繞射波低傾角成分；（f）（d）與（e）相減的結(jié)果；（g）（c）與（e）相加得到的最終繞射波分離結(jié)果；（h）原始數(shù)據(jù)的深度偏移剖面；（i）分離出的繞射波的深度偏移剖面．Fig．1 A theoretical model example of diffraction separation（a）Velocity model；（b）Planar decomposition section；（c）High－dip diffraction data separated from（b）by local dip filtering；（d）Result of subtracting（c）from（b）；（e）Low－dip diffraction data using prediction inversing；（f）Result of subtracting（e）from（d）；（g）Result of adding（e）to（c）；（h）PSDM result of original data；（i）PSDM result of diffraction data．

下面給出繞射波高、低傾角成分提取的具體方法．

2．1 局部傾角濾波及繞射波高傾角信息提取

采用去平面波濾波方法實現(xiàn)局部傾角濾波，進行繞射波高傾角信息提?。陬l率－空間域，相鄰道間時差為σ的二維平面波信號（n，ω）滿足關(guān)系

式中，n為道序下標(biāo)，ω為角頻率．Sergey（2002）給出了時移算子eiωσ的時間域近似TSσ（Zt）如下：

這里b－1＝（1－σ）（2－σ）/12，b0＝（2＋σ）（2－σ）/6，b1＝（1＋σ）（2＋σ）/12．顯然，利用（7）式進行時移運算可采用計算效率很高的追趕法．由式（6）、（7）可構(gòu)建時空域局部傾角濾波的二維濾波器

濾波器Fσ（Zx，Zt）可用于消除輸入數(shù)據(jù)中時間傾角為σ的信息成分．

進行局部傾角濾波，首先要給出傾角值，信號傾角σ由傾角掃描得到．記sp（x，τ）的離散形式為sp（j，k），這里k、j分別為時間和空間離散采樣下標(biāo)．設(shè)信號傾角σmin≤σ≤σmax，給定傾角掃描增量Δσ和時、空方向上的掃描時窗長度Lt、Lx．取σn＝σmin＋（n－1）Δσ∈［σmin，σmax］，en（j，k）為用Fσn（Zx，Zt）對sp（j，k）作濾波的結(jié)果，則由使

成立的n確定sp（j，k）的局部傾角σn（j，k）．

式（9）是單傾角掃描方法，可擴展到雙傾角掃描，以估計可能存在的兩個不同的信號傾角．設(shè)enm（j，k）為用Fσn（Zx，Zt）Fσm（Zx，Zt）對sp（j，k）作濾波的結(jié)果，則由使

成立的n、m確定sp（j，k）的兩個局部傾角σn（j，k）和σm（j，k）．雙傾角掃描計算量較大，為提高計算效率，可由單傾角掃描得到雙傾角中主信號傾角的估計值（j，k），然后以（j，k）為中心，縮小主傾角掃描范圍．

有了sp（j，k）的局部傾角參數(shù)σ（j，k），則在反射波同相軸的時間傾角限于（－σm，σm）范圍內(nèi)的假設(shè)下，由

即得到繞射波的高傾角信息成分的估計，這里Fσ是由（8）式給出的局部傾角濾波器．

濾波器（8）用于分離不同時間傾角的信號具有良好的效果，現(xiàn)給出一個用上述方法進行傾角分離的理論數(shù)據(jù)算例．圖2是一理論數(shù)據(jù)，其中（a）為振幅在空間上漸變的一組水平同相軸，（b）為振幅不變的三組正弦波形狀的彎曲同相軸，（c）是（a）與（b）疊加的結(jié)果．用上述的局部傾角濾波方法對（c）進行水平同相軸和彎曲同相軸分離，其結(jié)果如圖3所示．圖3a是分離出的水平同相軸信息，圖3b是分離出的彎曲同相軸信息，圖3c是圖3b與圖2b之差，即分離誤差．從結(jié)果可看到，水平與彎曲同相軸的大部分能量都得到了較好的分離，特別是在兩組同相軸傾角差較大時，分離精度較高．

2．2 預(yù)測反演及繞射波低傾角信息提取

即使是假設(shè)反射界面傾角較小，依靠局部傾角濾波也僅能獲得繞射波同相軸高傾角部分信息，而繞射波同相軸頂點附近的低傾角能量不能由反射波與繞射波同相軸的傾角差異進行分離．在共p剖面上，剔出高傾角信息成分后的繞射波的殘留低傾角信息（x，τ）以孤立的振幅異常與反射波疊加在一起，在頻率－空間域，利用反射波的線性可預(yù)測性可以把這種殘留的繞射波低傾角信息提取出來，將其加入到通過局部傾角濾波得到的繞射波高傾角信息成分中，就能夠得到相對完整的繞射波信息估計．

根據(jù)傅里葉理論我們知道，時空域的任何信號都能由該域上的有限個數(shù)（N個）的平面波信號疊加來逼近，這里需要的N越大，則反映信號在空間上的變化越劇烈．N個平面波的疊加結(jié)果在頻率空間域上具有線性預(yù)測關(guān)系，且預(yù)測算子長為N．頻率空間域的預(yù)測算子長度和預(yù)測濾波誤差反映信號在空間上的復(fù)雜程度，因此，信號的連續(xù)性可用線性可預(yù)測性來反映．在頻率空間域，反射波一般具有較強的線性可預(yù)測性，不可預(yù)測的成分可認(rèn)為是非反射的局部異常信息．

將（3）式變換到頻率空間域，得

對給定的頻率ω，記sk＝，S＝（s1，s2，…，sNx）T，rk＝Rp（xk，ω），R＝（r1，r2，…，rNx）T，dk＝（xk，ω），D＝（d1，d2，…，dNx）T，這里Nx為一個計算窗在空間方向上的道數(shù)．則有

設(shè)pl（l＝1，…，L）為S的線性預(yù)測算子，其由極小化目標(biāo)函數(shù)

得到，這里＊表示復(fù)共軛．用算子pl（l＝1，…，L）對S作濾波可得到R的估計，再由D＝S－R即得到局部異常信息D的估計，這就是通常的異常信息估計的預(yù)測濾波方法．預(yù)測濾波方法最大的缺陷是局部異常能量向鄰近道泄漏，使得異常信息估計精度不高，鑒于此，我們提出一種更有效的反演方法．

記Nx×Nx矩陣

對于待求的反射波R，一方面要求其要盡可能滿足線性預(yù)測關(guān)系，另一方面又要與輸入S接近，因此，取目標(biāo)函數(shù)

極小化上述目標(biāo)函數(shù)，可得到方程

這里PH為P的共軛轉(zhuǎn)置，I為單位矩陣，λ是權(quán)衡R的可預(yù)測性和R與S的背離程度的參數(shù)，λ越小，則越強調(diào)待求信號R的線性預(yù)測關(guān)系的滿足程度，相反，λ越大則越強調(diào)縮小R與輸入S之間的差異．求解（16）式便可得到反射信號R的一個估計．有了R，則由D＝S－R即得到D，最終將得到繞射波的低傾角信息成分（x，ω）的估計．

圖4是通過理論數(shù)據(jù)測試預(yù)測反演方法有效性的結(jié)果．從圖上可以看到，較傳統(tǒng)的預(yù)測濾波比，預(yù)測反演方法提取的局部異常具有較高橫向分辨率，更接近真實結(jié)果．

圖4 反演與濾波提取局部異常信息的效果對比（a）連續(xù)同相軸；（b）局部異常；（c）（a）與（b）相加的結(jié)果；（d）用預(yù)測反演方法從（c）中提取的局部異常；（e）用預(yù)測濾波方法從（c）中提取的局部異常．Fig．4 Comparison of local anomalies from inversion and filtering（a）Continuous events；（b）Abnormal signal；（c）Result of adding（b）to（a）；（d）Separated signal using inversion method；（e）Separated signal using filtering method．

在圖2給出的水平同相軸與彎曲同相軸分離測試中，如果對其中分離出的水平同相軸再用預(yù)測反演方法做局部異常提取，并將提取的局部異常加到圖2b的結(jié)果上，則得到用局部傾角濾波與預(yù)測反演聯(lián)合方法分離水平和彎曲同相軸的結(jié)果，如圖5所示．顯然，聯(lián)合方法得到的分解結(jié)果的精度有了顯著提高，預(yù)測反演方法的應(yīng)用使彎曲同相軸的零傾角位置的能量損失得到了有效補償．

圖5 采用局部傾角濾波與預(yù)測反演聯(lián)合方法得到的水平和彎曲同相軸分離結(jié)果（a）分離出的水平同相軸成分；（b）分離出的彎曲同相軸成分；（c）分離出的彎曲同相軸成分與真實結(jié)果（圖2b）之差．Fig．5 Results of separation combining local dip filtering and prediction inversion（a）Aclinal events；（b）Flexural events；（c）Difference between Fig．5band Fig．2b．

3 方法試算

圖6是根據(jù)塔里木盆地某探區(qū)實際解釋的地質(zhì)模型而建立的一個二維速度模型．模型中3700m 深度附近的界面代表奧陶系風(fēng)化面，其下部隨機分布了一些尺度、形狀不同的溶洞．溶洞的橫向尺度在5～100m，高度在10～200m 范圍內(nèi)變化，填充速度為3200～3900m/s，圍巖的速度為4400 m/s．對該模型進行數(shù)值模擬計算得到模擬疊前地震數(shù)據(jù)，再對模擬數(shù)據(jù)進行繞射波分離和疊前深度偏移處理．

圖6 速度模型Fig．6 Velocity model

圖7給出的是一個炮集記錄結(jié)果．從分離結(jié)果看，繞射波與反射波的絕大部分能量得到了正確的分離．

圖8是對該模型數(shù)據(jù)及分離出的繞射波進行疊前深度偏移的結(jié)果．從原始數(shù)據(jù)的成像結(jié)果可以看到，部分靠近強反射面，且繞射能量較弱的洞體成像后淹沒在強反射同相軸中，難以識別．在分離出的繞射波成像結(jié)果上，連續(xù)反射波能量得到了壓制，而斷點和孤立的洞體反射能量保留下來，顯然，該結(jié)果有助于對縫洞等地質(zhì)異常體的識別．

對溶洞分布區(qū)域，圖9給出了本文方法與Harlan信噪分離方法的結(jié)果對比．從結(jié)果可看到，在用Harlan信噪分離方法得到的繞射波成像結(jié)果上，與溶洞對應(yīng)的“串珠”狀反射在橫向上出現(xiàn)較強的旁瓣，這說明所提取的繞射波的空間帶寬（即低傾角信息）損失較大，而采用本文方法得到的繞射波成像結(jié)果與精確的繞射波成像結(jié)果更為接近．

圖10給出的是三維應(yīng)用實例，圖中（a）、（b）兩剖面分別取自原始數(shù)據(jù)和提取的繞射波信息的三維疊前時間偏移數(shù)據(jù)體，這里僅顯示了反映奧陶系風(fēng)化面反射（3．5s附近）的一段數(shù)據(jù)．顯然，在繞射波成像剖面上，來自連續(xù)反射面的反射信號得到了壓制，保留的反射信息與斷層、反射界面形態(tài)突變、溶洞等地質(zhì)現(xiàn)象相關(guān)，尤其是使疑似來自溶洞的孤立“串珠”狀反射更被突顯出來．顯然，這種繞射波成像結(jié)果可用于更好地進行縫洞儲集體的識別．

4 結(jié) 語

將繞射波從原始地震數(shù)據(jù)中分離出來單獨成像，其結(jié)果有助于提高碳酸鹽巖縫洞型儲層這類復(fù)雜非均質(zhì)儲層的識別能力．實現(xiàn)繞射波的有效分離是得到可靠的反映非均質(zhì)體特征的成像結(jié)果的關(guān)鍵，在地下反射界面傾角較小的前提下，局部傾角濾波和預(yù)測反演相結(jié)合是在平面波域提取繞射波信息的有效方法．僅依賴傾角差異分離繞射波會嚴(yán)重?fù)p失低傾角信息，配合以預(yù)測反演方法則可得到相對完整的繞射波信息，減少繞射波空間帶寬的損失，從而更好地保證繞射波成像結(jié)果的橫向分辨率．

（References）

［1］ Landa E，Keydar S．Seismic monitoring of diffraction images for detection of local heterogeneities．Geophysics，1998，63（3）：1093－1100．

［2］ Tcheverda V，Pozdnyakov V A．3Dfocusing transformation：reliable tool for imaging of scattering objects．∥Abstracts of SEG/New Orleans 2006Annual Meeting，2006：2564－2568．

［3］ Khaidukov V．Diffraction imaging by a focusing－defocusing approach．∥Expanded Abstracts of 73rdAnnual Internat SEG Mtg，2003：26－31．

［4］ Nowak E J，Imhof M G，Tech V．Diffractor localization via weighted Radon translations．∥Expanded Abstracts of 74thAnnual Internat SEG Mtg，2004：10－15．

［5］ Bansal R，Imhof M G．Diffraction enhancement in prestack seismic data．Geophysics，2005，70（3）：V73－V75．

［6］ Harlan W S，Claerbout J F，Rocca F．Signal/noise separation and velocity estimation．Geophysics，1984，49（11）：1869－1880．