王曉毅，陳璽,楊振,張江杰

1 中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京 100029 2 中國科學(xué)院地球科學(xué)研究院，北京 100029 3 南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室(廣州)，廣州 511458 4 廣州海洋地質(zhì)調(diào)查局，廣州 510075 5 中國科學(xué)院大學(xué)，北京 100049

0 引言

地下介質(zhì)中廣泛發(fā)育的斷層、尖滅、溶洞等中小尺度不連續(xù)體，對油氣生成、運(yùn)移和存儲具有重要意義.由于繞射波的形成與裂縫等構(gòu)造相關(guān)聯(lián)(黃建平等，2015；欒錫武和李繼光，2018)，其可以作為偏移處理對象來對目標(biāo)地質(zhì)體進(jìn)行精細(xì)成像.

目前，相關(guān)研究主要集中在如何利用繞射波與反射波在運(yùn)動學(xué)上的區(qū)別進(jìn)行波場分離與成像.常規(guī)的繞射波成像技術(shù)以射線類偏移方法為主.其中，針對繞射波較強(qiáng)的情形，對地震數(shù)據(jù)采用繞射增強(qiáng)的疊加算法是一種有效的策略(Faccipieri et al.,2016)，但其難以處理能量較弱的繞射波.利用反射波和繞射波不同的運(yùn)動學(xué)特征，在地震資料中對二者進(jìn)行分離，同樣能夠達(dá)到成像的目的(Merzlikin and Fomel,2017;Rad et al.,2018).然而，當(dāng)?shù)叵陆橘|(zhì)構(gòu)造復(fù)雜時(shí)，反射波在變換空間的形態(tài)不再簡單，極大地提高了分離的難度.此外，還可以在偏移道集中區(qū)分繞射波和反射波(Zhang and Zhang,2014)，該類方法對復(fù)雜反射具有很好的適應(yīng)性，但對于深度域三維成像而言，其運(yùn)算成本較高.Li和Zhang(2019)提出基于反穩(wěn)相偏移產(chǎn)生的時(shí)差道集來甄別和增強(qiáng)繞射波信號，該方法要求地層傾角信息精確，不需要過大內(nèi)存就可實(shí)現(xiàn)繞射波成像.

射線類方法在提取地震波的傳播方向和走時(shí)上具有靈活簡便的優(yōu)點(diǎn)，但在復(fù)雜構(gòu)造區(qū)域的成像上往往效果有限.基于此類方法的繞射波成像雖然提高了繞射體的成像分辨率，但其成像位置的準(zhǔn)確性難以滿足要求.利用波動類方法進(jìn)行繞射波成像精度較高，需要的計(jì)算成本也較大.例如，Liu等(2016)借助相移插值(PSPI)法提取傾角道集，并通過中值濾波實(shí)現(xiàn)了繞射點(diǎn)成像.逆時(shí)偏移作為一種經(jīng)典的波動方程類成像方法，具有諸多優(yōu)點(diǎn)，如不受地層傾角的限制，能夠靈活處理回轉(zhuǎn)波和棱柱波等.根據(jù)波傳播射線的幾何特征，選擇不同的波場組分進(jìn)行互相關(guān)運(yùn)算，可以實(shí)現(xiàn)對某一類構(gòu)造的單獨(dú)成像(王一博等，2016;Zhang et al.,2019a,b).當(dāng)入射波沿右下方傳播而出射波沿左上方傳播時(shí)，反射界面是正傾角的地層；入射波沿左下方傳播并且出射波沿右上方傳播時(shí)，反射界面是負(fù)傾角地層.由于在不連續(xù)點(diǎn)處產(chǎn)生的繞射是沿各個(gè)方向的，因此在多種成像條件下繞射波的能量都可以得到聚焦.Zhang等(2019a,b)通過一種乘法成像條件(the multiplication imaging condition)，利用正傾角地層剖面和負(fù)傾角地層剖面相乘得到單獨(dú)的繞射點(diǎn)位置剖面.這一方法較為簡單，但偏移剖面只包括位置信息，忽略了同相軸的幅值和相位信息.

本文基于逆時(shí)偏移對繞射波能量進(jìn)行聚焦，運(yùn)用多種方法得到單獨(dú)的繞射點(diǎn)剖面.首先，本文詳細(xì)分析了繞射波和反射波的運(yùn)動學(xué)特征以及各種成像方式的特殊性.借助Hilbert變換波場和Fourier變換，實(shí)現(xiàn)波場分離.基于此，我們獲取了單一傾向構(gòu)造的成像剖面以及更高質(zhì)量的傾角域共成像點(diǎn)道集.隨后，我們通過乘法成像條件獲取了單獨(dú)的繞射點(diǎn)剖面并在傾角域拾取了繞射波能量.此外，我們還對無需波場方向分解的正負(fù)傾角構(gòu)造成像公式進(jìn)行了推導(dǎo)，并提出了一種在成像剖面上快速獲取繞射點(diǎn)信息的方案.

1 基本理論

目前，基于逆時(shí)偏移對繞射點(diǎn)和間斷點(diǎn)實(shí)現(xiàn)定位的策略分為以下兩種：選擇不同的震源波場和檢波器波場組分互相關(guān)以得到正傾角和負(fù)傾角地層構(gòu)造，相乘獲取單獨(dú)的繞射點(diǎn)剖面(Zhang et al,2019a,b,2020)；逆時(shí)偏移產(chǎn)生共成像點(diǎn)傾角道集，然后在角道集中甄別、拾取繞射波能量(汪天池等,2020).以上兩種方案中，選擇合適組分和提高角道集質(zhì)量都可以利用Hilbert變換實(shí)現(xiàn)波場分離.

1.1 Hilbert變換與波場分解

在頻率-波數(shù)域中，根據(jù)圓頻率ω和波矢量k的符號，可以將原始波場按照傳播方向分解為若干個(gè)分量.在二維情形下，將時(shí)間-空間域的波場轉(zhuǎn)化到頻率-波數(shù)域需要三維傅里葉變換.其中，沿時(shí)間軸方向的傅里葉變換需要用到所有時(shí)刻的波場快照，這極大地增加了存儲成本并制約了運(yùn)行效率.

通過Hilbert變換對原始信號進(jìn)行改造，可以只保留頻譜和振幅譜中的正頻率部分.圖1a是均值為0.25，方差為0.02的高斯子波信號，其振幅譜和相位譜分別如圖1b藍(lán)色線和橙色線所示.將原始信號進(jìn)行Hilbert變換之后，得到以原始信號為實(shí)部，圖1c中橙色線段為虛部的復(fù)數(shù)信號.圖1d是復(fù)數(shù)信號所對應(yīng)的振幅譜(藍(lán)色線)和相位譜(橙色線)，可見復(fù)數(shù)信號在負(fù)頻率部分振幅為0，正頻率部分振幅擴(kuò)大為原來2倍.換言之，Hilbert變換使信號只保留正頻率部分.將這一性質(zhì)應(yīng)用于波場分解，可以避免圓頻率ω符號的判斷，從而減少沿時(shí)間軸方向波場的存儲.

圖1 原始信號與其經(jīng)Hilbert變換之后得到的復(fù)數(shù)信號(a)原始信號：均值為0.25，方差為0.02的高斯子波；(b)與信號(a)對應(yīng)的振幅(藍(lán)色線)和相位譜(橙色線)；(c)復(fù)數(shù)信號：以原始信號為實(shí)數(shù)部，以Hilbert變換后信號為虛數(shù)部(如橙色線段所示)；(d)與復(fù)數(shù)信號(c)對應(yīng)的振幅(藍(lán)色線)和相位譜(橙色線).注意振幅在負(fù)頻率部分為0.Fig.1 The original signal and the complex signal obtained after the Hilbert transform to it(a)The original signal:a Gaussian wavelet with mean value of 0.25 and variance of 0.02;(b)The amplitude and phase spectrum (blue and orange lines)corresponding to the signal (a);(c)The complex signal:taking the original signal as the real part and the Hilbert transformed signal (the orange line)as the imaginary part;(d)The amplitude and phase spectrum (blue and orange lines)corresponding to the signal (c).Note that the amplitude is zero in the negative frequency part.

對于聲波而言，標(biāo)量波動方程可以表述為

(1)

(2)

(3)

在這里，i=lu,ru,ld或者rd分別表示左-上行波，右-上行波，左-下行波和右-下行波.四者對應(yīng)的積分區(qū)域(kmin,kmax)分別為((0,-∞),(+∞,0))，((-∞,0),(0,-∞))，((+∞,0),(0,+∞))和((0,+∞),(-∞,0)).Re{·}表征取實(shí)數(shù)部算子.顯然，公式(3)也可以通過空間域Hilbert變換實(shí)現(xiàn).需要特別指出的是，盡管我們在論述的過程中選用的是標(biāo)量聲波方程，使用一階速度-應(yīng)力方程同樣可以先對震源進(jìn)行Hilbert變換后延拓得到伴隨波場，再實(shí)現(xiàn)波場分解.

1.2 針對不同地層結(jié)構(gòu)的成像條件

逆時(shí)偏移作為一種基于雙程波動方程延拓的方法，成像結(jié)果往往受低頻噪音以及偏移假象的影響.解決這一問題的途徑之一是在成像條件中只保留下行震源波場和上行檢波器波場的互相關(guān)(Claerbout,1971;Fei et al.,2015)，即

(4)

其中，s和r分別表示震源波場和檢波器波場，下標(biāo)u和d則分別代表上行和下行方向，tmax是最大記錄時(shí)間.按照波的傳播方向，震源波場和檢波器波場除了可以分為上下組分之外，還可以進(jìn)一步分為左右組分,即sd=sld+srd,ru=rlu+rru.因此，(4)式又可以進(jìn)一步寫為

(5)

利用震源波場和檢波器波場不同組分之間的互相關(guān)，可以實(shí)現(xiàn)對單一傾向的地質(zhì)構(gòu)造的成像.圖2a是一正傾角地層，地震波從震源激發(fā)向右下方向傳播至界面，沿左上方向發(fā)生反射并被地表檢波器接收.與之對應(yīng)，對于負(fù)傾角地層，地震波是向左下方傳播而反射波則沿右上方向傳播至檢波器，如圖2b所示.不同于傾斜的連續(xù)反射界面，繞射點(diǎn)可以同時(shí)滿足兩種情況下的入射波和反射波的方向特征.換言之，如果我們利用(5)式的第一項(xiàng)I1(x)，則只有正傾角地層和繞射點(diǎn)被成像；如果利用(5)式的第二項(xiàng)I2(x)，則只有負(fù)傾角地層和繞射點(diǎn)被成像.進(jìn)一步地，針對繞射點(diǎn)定位的乘法成像條件被提出(Zhang et al.,2019a,b)，即

圖2 針對正傾角(a)和負(fù)傾角(b)的地層構(gòu)造，入射波與出射波之間的關(guān)系Fig.2 The relationship between the incident wave and the outgoing wave for the formation structure with positive dip angle (a)and negative dip angle (b)

(6)

在這里，ns表示總炮數(shù).

不同于前兩者，采用I3(x)和I4(x)得到的構(gòu)造剖面沒有特定性，無論正、負(fù)傾角地層、平界面或者繞射點(diǎn)，都可能被成像.并且由互易性定理可知，當(dāng)觀測系統(tǒng)對地下介質(zhì)能夠充分覆蓋時(shí)，I3(x)和I4(x)的成像結(jié)果具有很好的相似性.

為了更好地利用多次波、回轉(zhuǎn)波、棱柱波等波型，我們還可以采用Causal成像條件(Liu et al.,2011;Li et al.,2019)，即

(7)

與之對應(yīng)，有

(8)

不同于(5)式，上式中的成像條件可以在不顯式對波場方向分解的情況下實(shí)現(xiàn).經(jīng)過推導(dǎo)(詳見附錄)，(8)式可以寫為

(9)

其中

(10)

(11)

1.3 高質(zhì)量傾角道集的提取

在地震資料處理中，角度域共成像點(diǎn)道集具有十分重要的作用.在張角域，反射波和繞射波信號之間沒有明顯的區(qū)別，因此區(qū)分二者能量比較困難.但是在傾角域，反射波能量只集中在地層的真實(shí)傾角附近，而繞射波的能量則連續(xù)分布在一定的范圍內(nèi).所以，通過先提取傾角道集再拾取繞射波能量的方式，可以凸顯繞射點(diǎn)的正確位置.Jin等(2014)將基于逆時(shí)偏移計(jì)算角道集的方法總結(jié)為三類：方向矢量法(DVB)，局部平面波分解法(LPWD)和局部相移成像條件法(LSIC).其中，方向矢量法運(yùn)算量少，方便快捷，并且結(jié)果分辨率高.為此，我們采用Poynting矢量法計(jì)算共成像點(diǎn)傾角道集.

本文選擇聲波方程的一階速度-應(yīng)力形式對波場進(jìn)行延拓，因此Poynting矢量可以直接由速度和應(yīng)力之積得到，即

p=-τv,

(12)

在這里，p表示Poynting矢量，τ表示應(yīng)力張量.v表示質(zhì)點(diǎn)的運(yùn)動速度矢量.Poynting矢量的幅值代表能流密度的大小，其方向指向波的傳播方向.由(12)式可知，對于任意一個(gè)空間位置，Poynting矢量只能給出一個(gè)傳播方向.所以在波發(fā)生重疊的區(qū)域，例如強(qiáng)反射界面附近，計(jì)算得到的角度準(zhǔn)確性不足.

利用Hilbert變換伴隨波場，首先將原始聲壓場和質(zhì)點(diǎn)速度分量分離為沿不同方向(左上，左下，右上和右下)的組分，再使用(12)式分別獲取不同波精確的傳播方向.雖然波場分解增加了方案的復(fù)雜度和運(yùn)算量，但能夠避免由波重疊引起的Poynting矢量不準(zhǔn)確的問題，從而提升傾角道集的質(zhì)量.

圖3描述了波由震源出發(fā)，在界面處發(fā)生反射最終被檢波器接收的過程.假設(shè)在成像點(diǎn)處發(fā)生反射，利用震源波場Poynting矢量ps和檢波器波場Poynting矢量pg可以計(jì)算出地層傾角α，即

圖3 傾角計(jì)算的幾何示意Fig.3 Geometric diagram of the dip-angle calculation

(13)

在偏移過程中，將每個(gè)時(shí)刻震源波場和檢波器波場的互相關(guān)值疊加到成像點(diǎn)對應(yīng)的α處即可得到共成像點(diǎn)傾角道集.

在傾角道集中，甄別、拾取繞射點(diǎn)能量既可以人工操作，也可以利用中值濾波等方法來自動實(shí)現(xiàn).值得注意的是，繞射能量也會受到偏移孔徑的影響，對于有限的空間采樣，繞射點(diǎn)埋深越大，其在傾角域分布的范圍越窄.而窄化的繞射點(diǎn)能量會弱化和反射點(diǎn)之間的區(qū)別，增大繞射點(diǎn)成像的難度.

1.4 成像后繞射點(diǎn)定位

前兩種方案都是以逆時(shí)偏移為基礎(chǔ)，根據(jù)地震波遇到不同地質(zhì)構(gòu)造時(shí)射線路徑的特征，實(shí)現(xiàn)繞射點(diǎn)定位和能量聚集的.因此，它們的過程相對復(fù)雜，存儲成本較高，運(yùn)算量成倍增加.如果能夠直接在成像剖面上拾取繞射點(diǎn)和間斷點(diǎn)，方法將不局限于逆時(shí)偏移技術(shù)，且效率也會有效提升.目前，已有相關(guān)研究展開.例如，奚先和黃江清(2020)通過卷積神經(jīng)網(wǎng)絡(luò)(CNN)實(shí)現(xiàn)了在地震剖面上散射體的定位和成像.本文提出另一種簡單直接的繞射點(diǎn)定位方案.

(14)

與式(6)相同，我們采用乘法成像條件獲取單獨(dú)的繞射點(diǎn)成像剖面，即

Idif(x)=I+(x)·I-(x).

(15)

顯然，本方案與方案一相比所需的運(yùn)算量較少、可以對工區(qū)剖面做局部化運(yùn)算、也不局限于逆時(shí)偏移所產(chǎn)生的地震剖面.

這里我們給出一個(gè)簡單的實(shí)例.圖4a是輸入圖像，包括左傾、右傾、水平和豎直四條線段和一個(gè)孤點(diǎn).圖4b是與4a相對應(yīng)的振幅譜，按照kx·kz≥0和kx·kz<0，我們分別得到了圖4c和4d.進(jìn)一步，經(jīng)過反傅里葉變換，我們得到了如圖4e和4f的構(gòu)造.需要注意的是，在圖4e中，除了傾角為正值的線段和孤點(diǎn)，還包括連續(xù)的水平線段和垂直線段，這是因?yàn)槲覀冊谄洳〝?shù)域中保留了kx·kz=0的區(qū)域.對于實(shí)際的地震剖面，為了避免在較平緩界面或者接近垂直的構(gòu)造處成像，我們會在波數(shù)域kx和kz接近0的位置將值設(shè)置為零.最后采用(15)式乘法運(yùn)算，我們得到了孤立點(diǎn)以及唯一的交叉點(diǎn)的位置，如圖4g中紅圈所示.

圖4 原始輸入圖像，不同構(gòu)造分離和繞射點(diǎn)成像圖(a)是原始圖像，包括不同傾向的線段和繞射點(diǎn).圖(b)是與圖(a)相對應(yīng)的幅度譜，按照(14)式可分為圖(c)和(d).經(jīng)過反傅里葉變換，可得圖(e)和(f).進(jìn)一步，按照乘法成像條件可得同時(shí)包含交叉點(diǎn)和孤點(diǎn)的最終剖面(g).Fig.4 Original input image,separation of different structures,and imaging of diffraction pointsPanel (a)is the original image,including line segments with different inclinations and an isolated point.Panel (b)is the amplitude spectrum corresponding to figure (a),which can be divided into panels (c)and (d)according to equation (14).After the inverse Fourier transform,panels (e)and (f)can be obtained.Furthermore,by using the multiplication imaging condition,we have a final profile (g)containing both an intersection point and an isolated point.

2 數(shù)值算例

為了驗(yàn)證上述理論的正確性和方案的可行性，我們在簡單模型和Marmousi2模型上分別進(jìn)行了實(shí)驗(yàn)并對結(jié)果進(jìn)行了展示.

2.1 簡單模型

在此模型中，我們設(shè)置了正、負(fù)傾角反射面、平界面以及三個(gè)繞射點(diǎn)，如圖5所示.40個(gè)炮點(diǎn)在地表均勻分布，間隔設(shè)為75 m，檢波器則分布在整個(gè)地表，間隔設(shè)為5 m.震源選擇為主頻30 Hz，時(shí)間采樣間隔0.5 ms的雷克子波.地震數(shù)據(jù)由時(shí)間二階、空間八階的有限差分算法模擬得到.該實(shí)驗(yàn)選擇3000 m·s-1的均勻模型作為偏移背景速度.

圖5 簡單模型，其中包括正負(fù)傾角界面，平界面以及三個(gè)繞射點(diǎn)Fig.5 A simple model.It includes positive-and negative-dip interfaces,a flat interface,and three diffraction points

在震源波場和檢波器波場中分別選取不同的組分采用互相關(guān)成像條件，可以得到不同的成像結(jié)果.圖6a和6b分別對應(yīng)式(5)中的I2(x)和I1(x)，即分別對負(fù)和正傾角構(gòu)造進(jìn)行了成像.除了連續(xù)界面，繞射點(diǎn)或者不連續(xù)點(diǎn)在剖面上也有所顯示.通過I1(x)和I2(x)剖面相乘，可以獲取只包含繞射點(diǎn)位置的成像結(jié)果，即圖6f.另外，我們在圖6c和6d中展示了I3(x)和I4(x)所對應(yīng)的剖面，由于震源點(diǎn)和檢波器點(diǎn)之間的互易性，兩者有一定的相似性.特別地，我們計(jì)算了以(4)式作為成像條件的剖面Isum(x)，如圖6e所示.按照第三種方案，我們采用獲取到的完整剖面作為輸入，分別提取了負(fù)傾角構(gòu)造和正傾角構(gòu)造，如圖7a和7b所示.應(yīng)當(dāng)注意的是，我們在波數(shù)域的分離避開了波數(shù)分量為0的區(qū)域，從而避免了單一傾角構(gòu)造剖面上的平界面的出現(xiàn).最后，通過乘法成像條件，獲取了圖7c中的繞射點(diǎn)剖面.從圖中可見，方案一和方案三在簡單的模型上的繞射點(diǎn)定位結(jié)果基本一致.圖8a和8b是按照成像條件(11)得到的負(fù)傾角和正傾角剖面，在求取過程中需對震源波場和檢波器波場進(jìn)行Hilbert變換，而無需伴隨波場的延拓和原始波場的方向分解.圖8c是由圖8a和8b得到的繞射點(diǎn)剖面.

圖6 簡單模型不同成像條件下的成像結(jié)果(a)是負(fù)傾角構(gòu)造，與成像條件I2(x)對應(yīng)；(b)是正傾角構(gòu)造，與成像條件I1(x)對應(yīng)；(c)和(d)分別是在成像條件I3(x)和I4(x)下的剖面；(e)完整的成像剖面；(f)繞射點(diǎn)剖面,為(a)和(b)的乘積.Fig.6 Imaging results of the simple model under different imaging conditions(a)is the negative-dip structure,which corresponds to the imaging condition I2(x);(b)is the positive-dip structure,which corresponds to the imaging condition I1(x);(c)and (d)are the profiles under imaging conditions I3(x)and I4(x),respectively;(e)is the complete profile;(f)is the profile of diffraction points,which is the product of (a)and (b).

圖7 簡單模型按照方案三得到的成像剖面(a)負(fù)傾角構(gòu)造；(b)正傾角構(gòu)造；(c)繞射點(diǎn)剖面,為(a)和(b)的乘積.Fig.7 Imaging profiles of the simple model obtained from scheme 3(a)The negative-dip structure;(b)The positive-dip structure;(c)The profile of diffraction points,which is the product of (a)and (b).

圖8 簡單模型按照式(11)中的成像條件得到的成像剖面(a)負(fù)傾角構(gòu)造；(b)正傾角構(gòu)造；(c)繞射點(diǎn)剖面,為(a)和(b)的乘積.Fig.8 Imaging profiles of the simple model according to image conditions in equation (11)(a)The negative-dip structure;(b)The positive-dip structure;(c)The profile of diffraction points,which is the product of (a)and (b).

在計(jì)算疊加剖面的同時(shí)，我們提取了共成像點(diǎn)傾角道集，并在圖9中展示了分別與I1、I2、I3、I4和Isum所對應(yīng)的角道集在橫向1500 m處的片段.從圖中可以看出，平界面在以I1和I2為成像條件提取的道集上沒有能量聚集，但是繞射點(diǎn)在四種成像條件下均可被呈現(xiàn).在疊加后的角道集9e上，繞射波能量連續(xù)分布在一段區(qū)域內(nèi)，而反射波的能量則聚焦在0°附近.以這些特征為依據(jù)，我們在傾角域手動拾取了繞射波能量.圖10展示了拾取的結(jié)果，繞射點(diǎn)和界面不連續(xù)點(diǎn)都得到了很好的呈現(xiàn)，且保留了同相軸的振幅和相位特征.至此，三種繞射點(diǎn)成像方案在簡單模型上都得到了良好的效果.

圖9 簡單模型在不同成像條件下的傾角道集(a),(b),(c)和(d)分別與成像條件 I1,I2,I3和I4對應(yīng)；(e)是(a),(b),(c)和(d)之和.Fig.9 Dip gathers for the simple model under different imaging conditions(a),(b),(c),and (d)correspond to the imaging conditions I1,I2,I3,and I4,respectively.(e)is the sum of (a),(b),(c),and (d).

圖10 簡單模型在傾角域手動拾取能量得到的繞射點(diǎn)剖面Fig.10 Diffraction point profile of the simple model obtained by manually picking up energy in the dip-angle domain

2.2 Marmousi2模型

相對復(fù)雜的Marmousi2模型的速度變化如圖11a所示.該模型在水平方向上有1325個(gè)采樣點(diǎn)，豎直方向上有467個(gè)采樣點(diǎn).震源埋深為5 m,總炮數(shù)為40，炮間距設(shè)置為165 m.而檢波器間隔設(shè)為5 m,在整個(gè)地表均有分布.地震數(shù)據(jù)同樣由時(shí)間二階、空間八階的有限差分算法模擬得到.偏移速度由真實(shí)模型經(jīng)過高斯光滑之后得到，如圖11b所示.

圖11 Marmousi2模型真實(shí)速度(a)與偏移所用的背景速度(b)Fig.11 True velocity of the Marmousi2 model (a)and background velocity used for migration (b)

利用震源波場的左下組分和檢波器波場的右上組分之間的互相關(guān)運(yùn)算I2，我們得到了如圖12a所示的負(fù)傾角構(gòu)造；與之相反，利用震源波場的右下組分和檢波器波場的左上組分之間的互相關(guān)運(yùn)算I1，我們得到了如圖12b所示的正傾角構(gòu)造.然后，借助正負(fù)傾角剖面的乘積，實(shí)現(xiàn)了繞射點(diǎn)和間斷點(diǎn)的定位，如圖12c所示.同樣，我們計(jì)算了采用I3和I4成像條件時(shí)的結(jié)果，即13a和13b，二者具有很好的相似性.

圖12 Marmousi2模型的負(fù)傾角剖面(a)和正傾角剖面(b)，以及最終的繞射點(diǎn)剖面(c)Fig.12 The negative-and positive-dip profiles of the Marmousi2 model,and the final diffraction point profile

圖13 Marmousi2模型在成像條件I3(a)和I4(b)下的結(jié)果，以及完整的剖面(c)Fig.13 (a)and (b)are the imaging profiles of the Marmousi2 model under imaging conditions I3 and I4,respectively.(c)is the complete profile

圖14 Marmousi2模型按照方案三得到的成像剖面(a)負(fù)傾角構(gòu)造；(b)正傾角構(gòu)造；(c)繞射點(diǎn)剖面,為(a)和(b)的乘積.Fig.14 Imaging profiles of the Marmousi2 model obtained from scheme 3(a)The negative-dip structure;(b)The positive-dip structure;(c)The profile of diffraction points,which is the product of (a)and (b).

圖15 Marmousi2模型按照式(11)中的成像條件得到的成像剖面(a)負(fù)傾角構(gòu)造；(b)正傾角構(gòu)造；(c)繞射點(diǎn)剖面,為(a)和(b)的乘積.Fig.15 Imaging profiles of the Marmousi2 model according to imaging conditions in equation (11)(a)The negative-dip structure;(b)The positive-dip structure;(c)The profile of diffraction points,which is the product of (a) and (b).

在圖16中，我們分別展示了與I1，I2，I3和I4相對應(yīng)的傾角道集的片段，其位于地表2500 m處.相較于簡單模型，Marmmousi2模型提取到的角道集噪音更加發(fā)育，繞射點(diǎn)的識別特征也不夠強(qiáng).在此，我們采用中值濾波對繞射波信號進(jìn)行了甄別和判定，得到了如圖17的結(jié)果.當(dāng)繞射點(diǎn)埋深較大且觀測系統(tǒng)孔徑有限時(shí)，繞射波信號與反射波信號并無顯著區(qū)別，為能量拾取帶來了一定的困難.因此，與圖12c相比，成像結(jié)果除在正確位置附近(如斷層處)比較清晰和準(zhǔn)確外，在其他區(qū)域也有部分能量殘余.

圖16 Marmousi2模型按照不同成像條件得到的傾角道集(a),(b),(c)和(d)分別與成像條件 I1,I2,I3和I4對應(yīng)；(e)是(a),(b),(c)和(d)之和.Fig.16 Dip-angle gathers for the Marmousi2 model under different imaging conditions(a),(b),(c),and (d)correspond to the imaging conditions I1,I2,I3,and I4,respectively.(e)is the sum of (a),(b),(c),and (d).

圖17 Marmousi2模型在傾角道集通過中值濾波得到的繞射點(diǎn)剖面Fig.17 Diffraction point profile of the Marmousi2 model is obtained by median filtering the dip-angle gathers

3 結(jié)論

本文以繞射點(diǎn)為成像目標(biāo)，分析了基于逆時(shí)偏移的各種繞射波成像方案.既利用了逆時(shí)偏移方法對地下構(gòu)造歸位準(zhǔn)確的優(yōu)點(diǎn)，又盡量減少計(jì)算成本.通過理論分析和數(shù)值算例驗(yàn)證，我們將三種方案總結(jié)如下：

(1)借助震源波場和檢波器波場的不同組分之間的互相關(guān)，可以實(shí)現(xiàn)不同構(gòu)造的成像.利用繞射點(diǎn)在正傾角剖面和負(fù)傾角剖面均能呈現(xiàn)的特點(diǎn)，通過乘法條件可以實(shí)現(xiàn)繞射點(diǎn)的定位.這一方法經(jīng)過與Hilbert變換和伴隨波場延拓相結(jié)合，存儲成本大大降低，計(jì)算效率得以提高；

(2)從Causal成像條件出發(fā)，推導(dǎo)正負(fù)傾角構(gòu)造成像公式，從而可以利用空間域的Hilbert變換避免震源波場和檢波器波場的方向分解；

(3)Poynting矢量法提取傾角道集進(jìn)而拾取繞射能量.為了提升角道集的質(zhì)量，我們采用先波場分離再計(jì)算傾角的策略，避免了Poynting矢量在波場重疊區(qū)域求取不準(zhǔn)的問題.這一方法的運(yùn)算量和存儲量都較大，并且當(dāng)偏移孔徑有限時(shí)，繞射波信號和反射波信號之間的區(qū)別會被弱化，從而在拾取上造成困難；

(4)通過地震剖面在波數(shù)域的分解，實(shí)現(xiàn)正負(fù)構(gòu)造剖面的提取.然后采取乘法運(yùn)算，確定繞射點(diǎn)和間斷點(diǎn)的位置.其中，為了避免在比較平緩或接近垂直界面上的錯誤成像，我們犧牲少許繞射點(diǎn)信息，將波數(shù)分量接近于0的區(qū)域值設(shè)置為0.相較于前兩種方案，這種圖像處理的方法具有以下優(yōu)點(diǎn)：運(yùn)算和存儲成本極低，不受偏移方法的限制，針對大工區(qū)可以只對目標(biāo)區(qū)域做局部化成像.當(dāng)然，本方法也需要信噪比較高的剖面作為輸入.

附錄A

由正文可知

(A1)

在頻率-波數(shù)域，上式可表示為

(A2)

因?yàn)?/p>

(A3)

在此，我們定義中間變量

(A4)

顯然，有

(A5)

同理，有

(A6)

(A7)

(A8)

將上式變換回空間-時(shí)間域，可得

(A9)

其中，

(A10)

上標(biāo)“±”表示取“+”或者取“-”.顯然，我們可以通過對空間進(jìn)行Hilbert變換實(shí)現(xiàn)式(A9),即

(A11)

其中，H算子定義為Hx(s)=s+iHx(s)，Hz(s)=s+iHz(s),Re{·}表示取實(shí)數(shù)部運(yùn)算.

與上述推導(dǎo)類似，我們得出

(A12)

這與Zhang等(2020)的結(jié)論一致.