楊城增 張宣堂 盛同杰 黎小偉 馬運(yùn)利 劉 韜

(①中國(guó)石化華北油氣分公司勘探開(kāi)發(fā)研究院，河南鄭州 466400；②中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083；③中國(guó)石化石油勘探開(kāi)發(fā)研究院，北京 100083)

0 引言

繞射波作為地下局部不連續(xù)體的地震響應(yīng)，攜帶了斷裂、溶洞、尖滅、陷落柱等小尺度非均質(zhì)體的重要信息，是提高地震勘探分辨率的重要工具[1]。然而，在地震記錄中，繞射波能量比反射波能量低1～2個(gè)數(shù)量級(jí)，在偏移成像時(shí)，與繞射相關(guān)的小尺度異常體常被強(qiáng)能量的反射體掩蓋，難以識(shí)別[2-3]；同時(shí)，在以反射波為有效信號(hào)的傳統(tǒng)地震數(shù)據(jù)處理流程中，繞射波也通常作為噪聲被壓制[4]。因此，對(duì)繞射波與反射波進(jìn)行分離，并對(duì)繞射波單獨(dú)成像，成為地下小尺度非均質(zhì)體高分辨率成像的關(guān)鍵。

20世紀(jì)50年代，地球物理學(xué)家Krey[5]在研究含有斷層或尖滅的反射地震記錄時(shí)發(fā)現(xiàn)了繞射波。近年來(lái)，為提高小尺度地質(zhì)體的成像分辨率，各種繞射波分離方法相繼被提出。根據(jù)分離繞射波時(shí)所處的地震數(shù)據(jù)處理階段，可以將繞射波分離方法分為以下兩類(lèi)。

(1)在偏移前進(jìn)行繞射波分離的方法。Nowak等[6]根據(jù)共炮點(diǎn)道集中繞射波與反射波同相軸橫向位置差異，使用加權(quán)拉東變換，在變換域中提取繞射波。平面波破壞濾波器(PWD)在繞射波分離中應(yīng)用廣泛，F(xiàn)omel[7]改進(jìn)了該方法，只需局部平面波場(chǎng)斜率一個(gè)參數(shù)?？籽┑萚8]根據(jù)繞射波與反射波在平面波記錄上的幾何形態(tài)差異，應(yīng)用PWD壓制反射波。朱生旺等[9]對(duì)共炮點(diǎn)道集記錄進(jìn)行平面波分解，首先利用局部?jī)A角濾波方法分離繞射波大傾角信息成分，再利用預(yù)測(cè)反演方法分離繞射波低傾角信息成分，將兩次分離結(jié)果相加得到最終分離的繞射波。魏巍等[10]在共炮點(diǎn)道集上通過(guò)動(dòng)校正拉平反射波，增加其相干性，而繞射波無(wú)法被拉平，相干性弱；然后根據(jù)繞射波與反射波相干性差異，利用奇異值分解技術(shù)提取繞射波。Rad等[11]將共反射波面元方法擴(kuò)展至疊前域，提出了一種基于波前屬性的疊前繞射分離流程。Gong等[12]結(jié)合Radon變換、Stolt偏移和反偏移，在炮域分離繞射波。此外，多道奇異譜分析[13]、相干波場(chǎng)疊加方法[14]、在線(xiàn)字典學(xué)習(xí)方法[15]等也被用于繞射波分離。

(2)在偏移后進(jìn)行繞射波分離的方法。Khai-dukov等[16]提出了聚焦—反聚焦方法分離繞射波。Moser等[17]利用反穩(wěn)相濾波器修正Kirchhoff偏移核函數(shù)，實(shí)現(xiàn)繞射波分離與成像。Silvestrov 等[18]利用逆時(shí)偏移生成共成像點(diǎn)道集，根據(jù)反射波與繞射波在共成像點(diǎn)道集上的傾角差異分離反射波。Zhao等[19]根據(jù)馬氏距離構(gòu)建了振幅衰減函數(shù)，據(jù)此對(duì)Kirchhoff核函數(shù)進(jìn)行修正，有效分離出繞射波。Yu等[20]基于局部平面波方程提出了數(shù)據(jù)規(guī)則化方法，實(shí)現(xiàn)了繞射波精確、穩(wěn)定提取。在方位傾角道集中，反射波具有稀疏特性，而繞射波則表現(xiàn)為低秩分量，根據(jù)這種差異，Zhao等[21]提出一種稀疏低秩矩陣分解方法提取繞射波。Li等[22]構(gòu)建了垂直旅行時(shí)差道集，在該道集中，繞射波是平坦的，而反射波表現(xiàn)為上凸形態(tài)，根據(jù)形態(tài)差異分離繞射波。汪天池等[23]證明了逆時(shí)偏移產(chǎn)生的傾角道集比Kirchhoff積分偏移產(chǎn)生的傾角道集具有更強(qiáng)的反射聚焦能力，更適合繞射波分離。羅騰騰等[24]在傾角域共成像點(diǎn)道集中，引入迭代收縮高分辨率Radon變換分離繞射波，具有良好的抗噪能力。劉培君等[25]通過(guò)構(gòu)建反穩(wěn)相濾波算子，在深度域進(jìn)行繞射波分離與成像。

傳統(tǒng)疊前繞射波分離方法大多在共炮檢距道集開(kāi)展，但由于深層繞射波能量分布范圍大，因此造成繞射波能量損失。針對(duì)此問(wèn)題，本文提出一種疊前繞射波分離方法，通過(guò)共虛震源變換將共炮點(diǎn)道集轉(zhuǎn)換成共虛震源道集，并在共虛震源道集中進(jìn)行偏移—反射波去除—反偏移處理，提取繞射波信號(hào)，解決在地震波場(chǎng)交叉和相切情況下繞射波弱信號(hào)難以分離的問(wèn)題，最后通過(guò)共虛震源逆變換得到共炮點(diǎn)道集繞射波場(chǎng)。

1 方法原理

1.1 共虛震源變換

以均勻介質(zhì)模型描述反射波和繞射波的旅行時(shí)特征(圖1)。根據(jù)反射定律，反射波旅行時(shí)可表示為

圖1 反射波和繞射波射線(xiàn)路徑及聚焦特性示意圖

(1)

式中：v為地震波傳播速度；h為炮檢距；d為炮點(diǎn)到反射界面法向距離；θ為界面傾角。在共炮點(diǎn)道集中，炮點(diǎn)坐標(biāo)固定，h變化，上式可寫(xiě)為

(2)

該旅行時(shí)曲線(xiàn)的頂點(diǎn)坐標(biāo)為(-2dsinθ，2dcosθ/v)，漸近線(xiàn)為±1/v。

繞射波旅行時(shí)可以表示為

(3)

式中：d1表示炮點(diǎn)至繞射體D的距離；φ為炮點(diǎn)至繞射體的射線(xiàn)角度。繞射波旅行時(shí)曲線(xiàn)的頂點(diǎn)坐標(biāo)為(-d1sinφ，(d1+d1cosφ)/v)，漸近線(xiàn)為±1/v。

根據(jù)式(2)和式(3)，繞射波(圖2的藍(lán)色箭頭所示)和反射波(圖2的紅色箭頭所示)都表現(xiàn)為雙曲線(xiàn)形態(tài)，在共炮點(diǎn)道集中分離繞射波具有一定難度。

圖2 共炮點(diǎn)道集中繞射波和反射波形態(tài)示意圖[26]

事實(shí)上，繞射波和反射波具有不同的聚焦性。如圖1a所示，反射波在虛震源SV處聚焦，其旅行時(shí)等于從虛震源至接收點(diǎn)的時(shí)間。繞射波則在真正的繞射體位置聚焦，其旅行時(shí)等于從炮點(diǎn)到繞射體、再到接收點(diǎn)的時(shí)間(圖1b)。在一個(gè)共炮點(diǎn)道集中，同一個(gè)反射界面產(chǎn)生的反射波可以視為來(lái)自于一個(gè)公共虛震源，而繞射波則不能，因此可以利用反射波旅行時(shí)，將反射波映射到對(duì)應(yīng)的虛震源位置上。通過(guò)這種方式，共炮點(diǎn)道集被轉(zhuǎn)換成一個(gè)新的道集——共虛震源道集(CVSG)。

假設(shè)地下某一位置點(diǎn)是虛震源，將共炮點(diǎn)道集轉(zhuǎn)換為CVSG可以表示為

δ{t-tR[S,SV(x,z),R]}dt

(4)

式中：x、z分別代表虛震源的水平位置和深度；R表示接收點(diǎn)；U表示記錄時(shí)間為t的共炮點(diǎn)道集數(shù)據(jù)；V是對(duì)應(yīng)的CVSG數(shù)據(jù)；在各向同性或等效介質(zhì)中，由震源S激發(fā)、R點(diǎn)接收、對(duì)應(yīng)(x,z)處虛震源的反射波走時(shí)可以簡(jiǎn)單地表示為

tR[S,SV(x,z),R]=

(5)

其中：t0是零炮檢距雙程旅行時(shí)，等于從虛震源SV到震源S的時(shí)間；W0為波前半徑；α為射線(xiàn)的出射角；vrms為均方根速度。在直射線(xiàn)假設(shè)條件下，有

(6)

式中xS表示震源的水平坐標(biāo)。根據(jù)式(6)，可求得式(4)中的反射波旅行時(shí)。值得注意的是，SV是一個(gè)假想的震源，并沒(méi)有真實(shí)的速度信息，使用均方根速度vrms計(jì)算從虛震源SV到接收點(diǎn)的旅行時(shí)。在簡(jiǎn)單情況下，反射波的出射角等于界面的真傾角，因此，可以用已知的實(shí)際傾角范圍限制SV的位置。

在共炮點(diǎn)道集向CVSG轉(zhuǎn)換的過(guò)程中存在一個(gè)問(wèn)題：如圖1a所示，當(dāng)接收點(diǎn)位于界面延伸至地表位置Rm左側(cè)時(shí)，可以計(jì)算從虛震源SV至接收點(diǎn)的反射旅行時(shí)；但是當(dāng)接收點(diǎn)位于Rm右側(cè)時(shí)，無(wú)法觀(guān)察到反射波。一般情況下，可以用S和SV計(jì)算Rm的近似位置處理這一問(wèn)題。轉(zhuǎn)換過(guò)程中存在的另一個(gè)問(wèn)題是：當(dāng)虛震源和接收點(diǎn)相距較遠(yuǎn)時(shí)，會(huì)產(chǎn)生拉伸畸變現(xiàn)象。如圖3和式(4)所示，CSVG中的數(shù)據(jù)點(diǎn)和共炮點(diǎn)道集中的數(shù)據(jù)點(diǎn)具有一一對(duì)應(yīng)的映射關(guān)系。當(dāng)時(shí)差Δt小于時(shí)間采樣間隔dt時(shí)，將會(huì)出現(xiàn)拉伸變形的現(xiàn)象。在均勻介質(zhì)中，時(shí)差Δt可以表示為

圖3 炮集U和共虛震源道集V的映射關(guān)系示意圖[23]j是時(shí)間t的離散序號(hào)

Δt=tR[S,SV(x,zi+1),R]-tR[S,SV(x,zi),R]

(7)

式中：i為深度z離散序號(hào)；LRSV(x,zi)表示接收點(diǎn)R到虛震源SV(x,zi)的距離。

當(dāng)虛震源和接收點(diǎn)相距越遠(yuǎn)，Δt的值越小，拉伸變形就越嚴(yán)重，在界面陡峭和淺層情況下尤為突出。參考動(dòng)校正的處理方法，當(dāng)拉伸畸變無(wú)法接受時(shí)，切除變形嚴(yán)重的部分。

共炮點(diǎn)道集和CVSG具有精確的映射關(guān)系，可以通過(guò)逆變換由CVSG恢復(fù)炮集數(shù)據(jù)

U(S,R,t)=

(8)

假設(shè)地下每一個(gè)點(diǎn)都是虛震源點(diǎn)，則可逐點(diǎn)給定坐標(biāo)對(duì)共炮點(diǎn)道集數(shù)據(jù)進(jìn)行共虛震源變換，來(lái)自同一反射面的反射事件將被校正到同一個(gè)虛震源。本文并不在虛震源位置進(jìn)行反射波的聚焦疊加，而是依次排列成為線(xiàn)性同相軸(圖4紅色箭頭所示)。繞射波無(wú)虛震源，所以不會(huì)在共虛震源道集中表現(xiàn)為水平線(xiàn)性，而是表現(xiàn)為彎曲狀(圖4藍(lán)色箭頭所示)。

圖4 CVSG中繞射波和反射波形態(tài)示意圖[26]

1.2 繞射波分離方法

利用Stolt偏移、正則化PWD方法和反偏移在CVSG中提取繞射波弱信號(hào)。

Stolt偏移方法是目前最快速的偏移算法，廣泛應(yīng)用于地震數(shù)據(jù)處理。為了便于Stolt偏移，對(duì)共虛震源道集進(jìn)行深時(shí)轉(zhuǎn)換，轉(zhuǎn)換后的CVSG在三維情況下可以寫(xiě)為V(x,y,z=0,t)，其傅里葉變換為

exp(ikxx+ikyy-iωt)dxdydt

(9)

則成像數(shù)據(jù)M(x,y,z)的傅里葉變換為

(10)

式中vm為偏移速度。

式(10)即為Stolt偏移公式。圖5為圖4的偏移結(jié)果，反射波表現(xiàn)為平面波，繞射波收斂為繞射點(diǎn)。

圖5 共虛震源道集偏移結(jié)果

在CVSG經(jīng)Stolt偏移后，反射波表現(xiàn)出線(xiàn)性連續(xù)特征具有平面波的性質(zhì)，符合平面波假設(shè)，可以通過(guò)PWD濾波器預(yù)測(cè)[7]。然而，在偏移數(shù)據(jù)中，繞射波收斂為繞射點(diǎn)，繞射波不符合平面波假設(shè)，難以通過(guò)PWD濾波器預(yù)測(cè)，在預(yù)測(cè)過(guò)程中，繞射波會(huì)作為預(yù)測(cè)誤差得到。因此，可以使用PWD方法得到偏移剖面中的繞射波

C(σ)M=Q(σ)

(11)

式中：M=(m1,m2,…,mN)，其中mi為第i道偏移數(shù)據(jù)，N為總道數(shù)；Q(σ)為預(yù)測(cè)誤差，σ為斜率矢量；C(σ)是PWD算子，定義為

(12)

其中：I為單位矢量；Bl,n(σl)表示由第l道偏移數(shù)據(jù)預(yù)測(cè)第n道偏移數(shù)據(jù)的預(yù)測(cè)誤差算子[7]。

從式(12)中可以看出，PWD算子所需要的唯一參數(shù)就是局部斜率σ，只需要求解局部斜率，再將其代入PWD算子，就可以利用式(11)得到繞射波。因此，獲得繞射波的關(guān)鍵就是估計(jì)局部斜率。

求解局部斜率是一個(gè)迭代過(guò)程，可以使用高斯—牛頓法[7]，在Q(σ)近似為零條件下，有

C′[σ(K)]Δσ(K)M+C[σ(K)]M=0

(13)

式中：K為迭代次數(shù)；C′[σ(K)]為C[σ(K)]關(guān)于σ(K)的微分；Δσ(K)為每次迭代中傾角斜率矢量的更新步長(zhǎng)，局部斜率矢量的更新公式為σ(K+1)=σ(K)+Δσ(K)。

上述是傳統(tǒng)的求解局部斜率的過(guò)程，在求解過(guò)程中，并未考慮與斷層等不連續(xù)地質(zhì)體有關(guān)的地震數(shù)據(jù)對(duì)反演過(guò)程干擾作用和地震反射信號(hào)平滑特征，反演求解不穩(wěn)定，局部斜率估計(jì)不準(zhǔn)確，應(yīng)用不準(zhǔn)確的局部斜率分離繞射波，可能會(huì)導(dǎo)致反射波殘留、分離出的繞射波不完整等問(wèn)題。因此，為了避免上述問(wèn)題，本文采用Yu等[20]提出的正則化局部斜率估計(jì)模型

ψ〈L1σ,σ〉+β〈L2σ,σ〉→min

(14)

式中：ψ、β均為正則化因子，用于平衡預(yù)測(cè)地震道和觀(guān)測(cè)地震道的擬合情況；〈·,·〉表示兩個(gè)向量的內(nèi)積；L1為負(fù)拉普拉斯算子，即

(15)

L2為差分算子，即

(16)

上述正則化傾角估計(jì)模型保證了反演求解穩(wěn)定性，也確保三維地震強(qiáng)反射信號(hào)估計(jì)準(zhǔn)確性[20]。

為快速求解上述局部斜率估計(jì)模型，目標(biāo)函數(shù)JΨ,β(σ)的梯度可表示為[20]

(17)

模型更新方案為

σ(K+1)=A{σ(K)-ξKg[σ(K)]}

(18)

其中：ξK為更新步長(zhǎng)；對(duì)于任意ξK>0，A(·)表示斜率矢量的可行集投影。

在得到局部斜率后，可由式(11)得到預(yù)測(cè)誤差Q(σ)，即偏移數(shù)據(jù)中的繞射波。

在利用上述方法去除偏移數(shù)據(jù)中的反射波后，再進(jìn)行反偏移

(19)

通過(guò)逆傅里葉變換可得CVSG繞射波

exp(-ikxx-ikyy+iωt)dkxdkydω

(20)

最后，可以通過(guò)式(8)將繞射波CVSG變換為繞射波共炮點(diǎn)道集。

本文方法的處理流程如圖6所示。

圖6 疊前繞射波分離流程

2 模型和實(shí)際數(shù)據(jù)應(yīng)用

構(gòu)建巖溶—沖溝模型驗(yàn)證本文繞射波分離方法的效果，并與常規(guī)的PWD方法對(duì)比，闡明本文方法的優(yōu)勢(shì)；通過(guò)實(shí)際資料處理驗(yàn)證本文方法的適用性。

2.1 巖溶—沖溝模型

巖溶-沖溝模型(圖7)尺寸為9000m×4000m，1000～2300m深度范圍有兩個(gè)薄層，2000～3500m深度范圍內(nèi)隨機(jī)分布著一些形狀、尺寸不同的巖溶體，最小的直徑為4m，最大的直徑為100m，填充速度為4200m/s，圍巖的速度為4700m/s。圖8為模型的反射系數(shù)剖面，其中小尺度巖溶構(gòu)造清晰，可用于對(duì)比分析不同分離方法的效果。對(duì)該模型使用Kirchoff正演算法生成合成數(shù)據(jù)，震源使用主頻為25Hz的Ricker子波；共901炮，炮點(diǎn)位置從0到9000m變化，炮間距為10m；每炮901道，道間距為10m，采樣間隔為4ms，其共炮檢距剖面如圖9a所示。

圖7 巖溶—沖溝速度模型

圖8 巖溶—沖溝模型的反射系數(shù)剖面

應(yīng)用PWD法對(duì)圖9a所示的共炮檢距剖面進(jìn)行繞射波分離，結(jié)果如圖9b所示。PWD方法有效去除了平滑、連續(xù)的反射波，但繞射波也遭到嚴(yán)重的損傷，尤其在繞射波頂點(diǎn)，即繞射波與反射波相干處。

對(duì)模擬數(shù)據(jù)應(yīng)用本文方法分離繞射波，并將分離結(jié)果(共炮點(diǎn)道集)抽成共炮檢剖面，如圖9c所示。本文方法有效去除了剖面中平滑、連續(xù)的反射波，比PWD方法更完整地保留了繞射波信息，為繞射波的清晰成像奠定了基礎(chǔ)。

圖9 巖溶—沖溝模型兩種繞射波分離方法的共炮檢距剖面對(duì)比

圖10a是PWD法分離的繞射波偏移成像結(jié)果，雖然繞射波弱信號(hào)在一定程度上突顯，但由于分離的繞射波不完整，繞射點(diǎn)成像不清晰，信噪比較低。圖10b是使用本文方法分離的繞射波偏移成像結(jié)果，得益于反射波去除干凈、繞射波保留完整，繞射點(diǎn)成像清晰完整，小尺度巖溶能夠突顯。

圖10 巖溶—沖溝模型兩種方法分離的繞射波偏移剖面對(duì)比

2.2 實(shí)際資料應(yīng)用

實(shí)際資料采集于內(nèi)蒙古杭錦旗。該區(qū)奧陶系縫洞型儲(chǔ)集體發(fā)育，具有埋藏深、尺度小、非均質(zhì)性強(qiáng)等特點(diǎn)，繞射波響應(yīng)較弱。在常規(guī)反射成像剖面(圖11a)上縫洞型儲(chǔ)集體反映不清晰，精細(xì)識(shí)別與刻畫(huà)難度大，目的層的沖溝和巖溶體有顯示(紅色圓圈所示)。由于小尺度巖溶體的存在，平滑線(xiàn)性的反射層呈現(xiàn)向下彎曲的形態(tài)，但由于背景強(qiáng)反射的影響，分辨率較低，很難識(shí)別其中的小尺度巖溶體。對(duì)實(shí)際資料使用本文方法分離繞射波，并對(duì)分離結(jié)果進(jìn)行成像(圖11b)。由于強(qiáng)反射已被移除，原來(lái)被強(qiáng)反射掩蓋的小尺度巖溶突顯(紅色圓圈內(nèi))，分辨率和聚焦性顯著提高。相干體屬性可用于斷裂識(shí)別，圖12a是實(shí)際資料的反射波相干屬性切片，雖然可以看到巖溶體和斷裂，但小尺度巖溶體和斷裂難以識(shí)別。圖12b為繞射波能量屬性切片，可識(shí)別的小尺度巖溶體和斷裂明顯增多，小尺度巖溶體和斷裂的邊界、形態(tài)也更清晰(紅色圓圈內(nèi))。

圖11 實(shí)際數(shù)據(jù)成像剖面對(duì)比

圖12 反射波相干屬性(a)與繞射波能量屬性(b)的切片對(duì)比

3 結(jié)論

本文通過(guò)共虛震源變換利用共炮點(diǎn)道集生成CVSG，對(duì)CVSG應(yīng)用Stolt偏移使雙曲形態(tài)的繞射波能量收斂，使用正則化PWD去除反射波，再利用反偏移重新得到CVSG，最后利用共虛震源反變換得到共炮點(diǎn)道集繞射波場(chǎng)。本文方法解決了反射與繞射波相切或相交時(shí)繞射波難以分離的問(wèn)題，避免了傳統(tǒng)方法在共炮檢距道集上進(jìn)行分離造成的深層繞射波能量損失。數(shù)值模型測(cè)試驗(yàn)證了本文提出的疊前繞射波分離方法可以提升小尺度地質(zhì)體成像分辨率，實(shí)際資料的應(yīng)用進(jìn)一步證明了本文方法對(duì)探明縫洞型儲(chǔ)集體具有積極意義。本文方法中共虛震源正、反變換雖然數(shù)學(xué)上是可逆的，但實(shí)現(xiàn)過(guò)程中有能量損失，如何減小這種損失是下一步的研究方向。