劉俊成， 陳祥忠 ， 郭井學(xué)， 王瑞星， 惠夢琳， 黃申碩

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026 2 北京桔燈地球物理勘探股份有限公司， 北京 102200 3 中國極地研究中心， 上海 200129

0 引言

地震偏移成像是地震數(shù)據(jù)處理過程中的重要環(huán)節(jié)(Schneider, 1978; Stolt et al., 1987; Baysal et al., 1983; Claerbout, 1992; Gray et al., 2000; 王華忠等，2009；孫小東等，2020)，其不但可以用于恢復(fù)地下的構(gòu)造信息，還可以估算地下的巖性和物性參數(shù)用于后續(xù)的地震反演和儲(chǔ)層預(yù)測.按照偏移速度的不同定義方式，地震偏移方法可以分為兩大類：時(shí)間偏移和深度偏移.深度偏移可以適應(yīng)地下介質(zhì)速度的劇烈橫向變化，是復(fù)雜介質(zhì)成像的首選技術(shù)手段，但其實(shí)際應(yīng)用效果嚴(yán)重依賴于速度建模的精度.時(shí)間偏移雖然無法對(duì)復(fù)雜構(gòu)造進(jìn)行準(zhǔn)確的成像，但當(dāng)速度橫向變化相對(duì)平緩時(shí)，其依然可以滿足大部分工區(qū)的成像處理需要.此外，相對(duì)于深度偏移，時(shí)間偏移的計(jì)算效率更高，所需偏移速度也更容易估算.因此，即便在計(jì)算機(jī)硬件高度發(fā)達(dá)的今天，時(shí)間偏移依然是業(yè)界廣泛使用的地震成像方法.

眾所周知，常規(guī)的偏移方法只是地震波場正演算子的共軛轉(zhuǎn)置(Tarantola, 1984; Claerbout, 1985)，對(duì)于有限觀測系統(tǒng)所采集的帶限地震數(shù)據(jù)，其只能生成模糊的地下構(gòu)造成像結(jié)果.此外，復(fù)雜的地下介質(zhì)構(gòu)造和非規(guī)則的地震數(shù)據(jù)空間采樣，還會(huì)導(dǎo)致偏移假象和非規(guī)則的成像照明，嚴(yán)重影響成像振幅信息的可靠性.基于線性化地震反演理論的最小二乘偏移(Least-Squares Migration, LSM)是解決上述問題的有效手段(Nemeth et al., 1999; Duquet et al., 2000; Dai and Schuster, 2013; Hu et al., 2016; Yue et al., 2019；周東紅等，2020；劉國章等，2020；劉玉敏等，2021).其將地震偏移視為線性反問題，并利用最優(yōu)化方法求取能夠準(zhǔn)確模擬采集地震數(shù)據(jù)的成像剖面作為地下真實(shí)反射率，在理論上能夠消除非規(guī)則采集、帶限子波等因素對(duì)成像結(jié)果的不利影響.但是，LSM每次迭代均需一次正演和偏移運(yùn)算，完整的LSM反演過程往往需要幾十倍于常規(guī)偏移的計(jì)算量，龐大的計(jì)算成本嚴(yán)重影響了該技術(shù)的應(yīng)用前景(張宇，2018).

本文將地震數(shù)據(jù)的局部平面波合成/分解策略(Hill, 2001; Gray, 2005; Liu and Palacharla, 2011; 岳玉波等，2012；李江等，2016；袁茂林等，2016；陳祥忠等，2020)應(yīng)用于Kirchhoff正演/偏移(KFM/KTM)中，推導(dǎo)了高效的Kirchhoff 射線束正演/偏移(KBFM/KBTM)算子，并以此為基礎(chǔ)發(fā)展了一種快速的最小二乘Kirchhoff射線束疊前時(shí)間偏移(LSKBTM)方法.同需要逐道映射運(yùn)算的最小二乘Kirchhoff疊前時(shí)間偏移(LSKTM)相比，LSKBTM中數(shù)據(jù)空間和成像空間之間的數(shù)據(jù)映射運(yùn)算僅需在稀疏的射線束中心位置進(jìn)行，因此其計(jì)算成本大幅度降低.文中利用二維模型和二維、三維實(shí)際數(shù)據(jù)對(duì)LSKBTM進(jìn)行了測試，其結(jié)果證明，LSKBTM可以有效的提高偏移結(jié)果的分辨率和照明度，其反演成像精度和殘差收斂速度同LSKTM相當(dāng)，但計(jì)算效率得到了大幅度提升.

1 方法原理

1.1 LSM基本原理

LSM可被視為一個(gè)最優(yōu)化反演問題，其通過迭代的方式求取能夠準(zhǔn)確模擬接收記錄的成像剖面，并將其作為地下真實(shí)的反射率(Tarantola, 1984).LSM的目標(biāo)函數(shù)可以定義為

J=‖LR-d‖2+μ‖R‖2，

(1)

式中右側(cè)第一項(xiàng)為數(shù)據(jù)匹配項(xiàng)，用于對(duì)地下反射率R的反演更新，第二項(xiàng)為保證反演穩(wěn)定引入的正則化項(xiàng)，μ為控制正則化程度的參數(shù)，d為接收地震記錄，L為線性化波場正演算子.式(1)所示的LSM目標(biāo)函數(shù)一般通過梯度類算法(如共軛梯度法)進(jìn)行迭代求解，其核心在于正演算子L以及共軛算子LT的構(gòu)建.接下來，本文在時(shí)間域KFM/KTM算子引擎的基礎(chǔ)上，通過引入地震數(shù)據(jù)的局部平面波分解/合成策略，推導(dǎo)高效的KBFM/KBTM算子，并以此為基礎(chǔ)實(shí)現(xiàn)LSKBTM.

1.2 KFM/KTM算子

基于體積分表征的Kirchhoff正演公式(Bleistein, 2001)可以表示為

D(rd,rs,ω)=

其中，Ω0代表地下散射點(diǎn)的集合，D(rd,rs,ω)為所模擬的單炮地震記錄，F(xiàn)(ω)是震源頻譜，R(r0)為地下散射點(diǎn)r0=(x0,y0,t0)處的反射率，G(r0,r′,ω)為震源在r′=(x′,y′,0)、觀測點(diǎn)在r0處的格林函數(shù)，其在高頻射線理論下可以表示為

G(r0,r′,ω)=A(r′)exp[iωT(r′)]，

(3)

其中，T(r′)為地震波由震源到觀測點(diǎn)的單程走時(shí)，在時(shí)間偏移中一般應(yīng)用單平方根公式進(jìn)行計(jì)算：

(4)

A(r′)為傳播振幅，其在時(shí)間域的簡化表達(dá)形式可以參照Zhang等(2000)中的相關(guān)內(nèi)容.

×exp{iω[T(rs)+T(rd)]}，

(5)

同時(shí)，求取式(5)的共軛轉(zhuǎn)置可以得到相應(yīng)的KTM算子：

(6)

由式(5)和式(6)不難看出，KFM/KTM需要逐道進(jìn)行地震數(shù)據(jù)和成像剖面間的數(shù)據(jù)映射運(yùn)算，其計(jì)算量與地震數(shù)據(jù)的總道數(shù)呈正比關(guān)系.由于采集地震數(shù)據(jù)往往存在密集的空間采樣，因此根據(jù)局部平面波合成/分解策略可將其轉(zhuǎn)化為稀疏射線束位置的局部平面波，以此為基礎(chǔ)實(shí)現(xiàn)KBFM/KBTM可以顯著降低地震數(shù)據(jù)和成像結(jié)果間的映射運(yùn)算次數(shù)，從而降低LSM的計(jì)算量.

1.3 KBFM算子

首先，應(yīng)用基于高斯窗的單位分解算法(Hill, 2001; Hu et al., 2016; Yue et al., 2019)對(duì)式(5)兩側(cè)進(jìn)行加窗處理：

(7)

其中，L代表高斯窗中心位置，也稱為射線束中心，其采樣間隔ΔL=(ΔLx,ΔLy)，ωr為參考頻率.此時(shí)，式(5)可以轉(zhuǎn)化為

(8)

接下來，用射線束中心振幅A(L)來近似式(8)中的振幅項(xiàng)A(rd)，并走時(shí)T(L)的一階泰勒展開來近似(8)中的走時(shí)項(xiàng)的T(rd)：

(9)

此時(shí)，式(8)可以轉(zhuǎn)化為

×exp{iω[T(rs)+T(L)]}

(10)

如果將對(duì)應(yīng)接收點(diǎn)射線參數(shù)為pd的散射點(diǎn)的集合定義為Ω0(pd)，那么式(10)可以表示為

(11)

其中：

×exp{iω[T(rs)+T(L)]}.

(12)

最后，分別對(duì)式(11)和式(12)應(yīng)用傅氏反變換，得到最終的時(shí)間域KBFM算子：

?s(L,pd,t-pd·(rd-L))，

(13)

其中，d(rd,rs,t)為模擬炮集地震記錄，符號(hào)?代表褶積運(yùn)算，γ(rd,t)為高斯窗的時(shí)間域響應(yīng)函數(shù)，s(L,pd,t)為在射線束中心處合成的局部平面波：

×A(L)δ(t-T(rs)-T(L))，

(14)

其中，δ為delta函數(shù)，f(t)為式(15)所定義的時(shí)間域子波信號(hào)：

(15)

利用KBFM進(jìn)行單炮模擬的計(jì)算過程可以大致概括為：

(1)確定射線束中心寬度ΔL，初始束寬度w0，射線參數(shù)采樣間隔Δpd等計(jì)算參數(shù).

(2)對(duì)于每一個(gè)射線束中心L，進(jìn)行地下散射點(diǎn)的循環(huán)，計(jì)算式(14)中的走時(shí)、振幅等參數(shù)，并以此為基礎(chǔ)將地下反射率映射轉(zhuǎn)換到τ-p域.

(3)重復(fù)步驟(2)直到所有射線束中心處理完成，然后計(jì)算τ-p域地震道同子波f(t)的褶積得到局部平面波分量s(L,pm,t).

(4)利用式(13)所示的逆傾斜疊加將所合成的平面波分量累加到地震記錄d(rd,rs,t)中，待所有射線束中心處理完成后，即可得到最終的模擬炮記錄.

1.4 KBTM算子

對(duì)式(6)所示的KTM算子應(yīng)用類似的走時(shí)近似和加窗單位分解算法，可以得到如下的KBTM算子：

+T(L)),

(16)

(17)

(18)

KBTM的單炮偏移計(jì)算過程可以簡要概括為：

(1)沿射線束中心循環(huán)，并在不同的射線束中心處利用式(17)和式(18)將單炮地震記錄分解為局部平面波分量.

(2)對(duì)于每一個(gè)成像點(diǎn)，求取由震源經(jīng)地下成像點(diǎn)到射線束中心的雙程地震波走時(shí)和接收點(diǎn)射線參數(shù)等信息，根據(jù)式(16)拾取平面波振幅并累加到成像點(diǎn)上.

(3)重復(fù)步驟(1)和(2)，待所有射線束中心處理完成后，即可得到單炮成像結(jié)果.

2 模型及實(shí)際資料試算

本節(jié)首先應(yīng)用Marmousi模型進(jìn)行LSKTM和LSKBTM的測試，證明兩者具有非常接近的成像精度，并驗(yàn)證LSKBTM在計(jì)算效率上的優(yōu)勢.接下來，應(yīng)用二維和三維探區(qū)地震數(shù)據(jù)測試LSKBTM的實(shí)際應(yīng)用效果.

2.1 Marmousi模型

首先，將Marmousi模型的層速度場通過Dix公式轉(zhuǎn)換為圖1a所示的均方根速度場，為便于后續(xù)的成像對(duì)比，在此直接將層速度場計(jì)算得到反射率剖面作為時(shí)間域反射率(圖1b)；然后，將時(shí)間域均方根速度和反射率作為輸入數(shù)據(jù)，利用KFM計(jì)算了200炮地震記錄(部分炮記錄如圖2所示)，炮間隔為40 m，每炮接收地震道數(shù)為151，道間距為20 m，震源信號(hào)為主頻20Hz的雷克子波.

圖1 時(shí)間域Marmousi模型 (a) 均方根速度場; (b) 時(shí)間域反射率剖面.Fig.1 Marmousi model in time domain(a) RMS velocity field； (b) Reflectivity profile.

圖2 不同位置處的KFM單炮記錄Fig.2 Single-shot gathers using KFM at different source locations

分別使用LSKTM和LSKBTM(初始束寬度設(shè)為240 m)對(duì)模擬記錄進(jìn)行偏移處理，經(jīng)過1次迭代得到的KTM和KBTM結(jié)果分別如圖3a和圖3b所示，圖4同時(shí)展示了兩者在CDP=420處的成像波形同真實(shí)反射率的疊合對(duì)比(KTM藍(lán)色，KBTM綠色，真實(shí)反射率黑色).可以看到KTM和KBTM均較好的恢復(fù)了地下的構(gòu)造信息，成像效果非常接近.但由于KTM和KBTM只是正演過程的共軛轉(zhuǎn)置，同圖1b所示的真實(shí)反射率相比，其成像分辨率較低，成像振幅也不準(zhǔn)確.經(jīng)過12次迭代后的LSKTM和LSKBTM結(jié)果分別如圖5a和圖5b所示，相應(yīng)的成像波形對(duì)比結(jié)果如圖6所示，可以看到兩種LSM方法的成像結(jié)果幾乎完全相同，均較好的恢復(fù)了地下真實(shí)的反射率信息，成像分辨率和照明度較常規(guī)偏移結(jié)果有了大幅度改善.由圖7所示的反演收斂曲線對(duì)比可以看出，兩種方法也具備相似的收斂速度.使用相同的計(jì)算節(jié)點(diǎn)對(duì)比兩種方法的計(jì)算時(shí)間，其中，LSKTM的總計(jì)算時(shí)間為361.9 s，LSKBTM的計(jì)算時(shí)間為64.1 s，因此LSKBTM相對(duì)于常規(guī)LSKTM的加速比約為5.7，也就是說，LSKBTM僅需要約4次常規(guī)KTM的計(jì)算量即可完成12次LSKBTM的迭代運(yùn)算.

圖3 常規(guī)偏移成像剖面對(duì)比(a) KTM成像結(jié)果; (b) KBTM成像結(jié)果.Fig.3 Comparison of migrated images(a) KTM image； (b) KBTM image.

圖4 CDP=420處成像波形對(duì)比其中紅色曲線代表真實(shí)反射率，藍(lán)色曲線代表KTM，綠色曲線代表KBTM.Fig.4 Comparison of waveforms at CDP=420Red curve is true reflectivity. Blue curve is KTM. Green curve is KBTM image.

圖5 最小二乘偏移成像剖面對(duì)比(a) LSKTM成像剖面； (b) LSKBTM成像剖面.Fig.5 Comparison of least-squares migration images(a) LSKTM image； (b) LSKBTM image.

圖6 CDP=420處成像波形對(duì)比其中紅色曲線代表真實(shí)反射率，藍(lán)色曲線代表LSKTM，綠色曲線代表LSKBTM.Fig.6 Comparison of waveforms at CDP=420Red is true reflectivity, blue is LSKTM image, and green is LSKBTM image.

圖7 歸一化目標(biāo)函數(shù)收斂曲線其中藍(lán)色曲線代表LSKTM，綠色曲線代表LSKBTM.Fig.7 Normalized misfits versus number of iterationsBlue is LSKTM, and green is LSKBTM.

2.2 二維實(shí)際數(shù)據(jù)

使用某二維探區(qū)實(shí)際地震數(shù)據(jù)進(jìn)行LSKBTM的測試，圖8展示了時(shí)間域偏移速度分析后得到的均方根速度場.該數(shù)據(jù)共有220炮，每炮271道，炮間隔為40 m，道間距為20 m，其中第100炮接收記錄如圖9a所示.測試使用的震源信號(hào)為主頻為20 Hz的雷克子波，選擇的束初始寬度為240 m.

圖8 均方根速度場Fig.8 RMS velocity field

應(yīng)用LSKBTM對(duì)數(shù)據(jù)進(jìn)行了10次迭代成像運(yùn)算，對(duì)應(yīng)的數(shù)據(jù)殘差得到了較好的收斂(降低約83%).經(jīng)過1次(KBTM)和10次LSKBTM反演迭代后的成像結(jié)果分別如圖9和圖10所示，可以看到，同KBTM相比，LSKBTM大幅提升了成像層位的分辨率和連續(xù)性.圖11進(jìn)一步對(duì)比了兩者的成像頻譜，可以看到同KBTM相比，LSKBTM有效拓寬了成像剖面的頻帶范圍.KBTM的成像效果同常規(guī)KTM近乎效同(未展示)，但計(jì)算時(shí)間僅為KTM的1/6，因此在本例中，LSKBTM相對(duì)于LSKTM的加速比約為6.

圖9 KBTM成像結(jié)果Fig.9 Imaging result using KBTM

圖10 LSKBTM成像結(jié)果Fig.10 Imaging result using LSKBTM

圖11 成像頻譜對(duì)比Fig.11 Spectrum comparison between KBTM and LSKBTM

2.3 三維實(shí)際數(shù)據(jù)

使用某三維探區(qū)實(shí)際地震數(shù)據(jù)進(jìn)行LSKBTM的測試，圖12展示了時(shí)間域偏移速度分析后得到的均方根速度場.應(yīng)用LSKBTM對(duì)數(shù)據(jù)進(jìn)行了10次迭代成像運(yùn)算，其中經(jīng)過1次迭代的KBTM結(jié)果和10次迭代的LSKBTM結(jié)果分別如圖13和圖14所示，可以看到同KBTM偏移剖面相比，LSKBTM大幅提升了成像的分辨率，時(shí)間切片中展示的地質(zhì)構(gòu)造更加的清晰聚焦.圖15進(jìn)一步對(duì)比了兩個(gè)成像剖面的頻譜(單位為dB)，可以看到同KBTM(紅色曲線)相比，LSKBTM(藍(lán)色曲線)不但大幅度提高了成像結(jié)果的主頻，有效頻帶寬度也得到大幅度拓寬.由于三維炮記錄在crossline方向的采樣稀疏(200 m)，本例僅在inline方向(采樣間隔20 m)進(jìn)行了平面波的合成和分解(束中心間隔為200 m)，同常規(guī)LSKTM相比，LSKBTM的加速比為5.3.

圖12 均方根速度場Fig.12 3D view of RMS velocity field

圖13 KBTM成像結(jié)果Fig.13 Image result of KBTM

圖14 LSKBTM成像結(jié)果Fig.14 Image result of LSKBTM

圖15 KBTM(藍(lán)色曲線)和LSKBTM(紅色曲線)的成像頻譜對(duì)比Fig.15 Spectrum comparison between KBTM (blue) and LSKBTM (red)

3 結(jié)論

本文基于地震數(shù)據(jù)的局部平面波分解/合成策略，在時(shí)間域KFM/KTM的基礎(chǔ)上，實(shí)現(xiàn)了高效的KBFM/KBTM算子，并以此為基礎(chǔ)發(fā)展了快速的LSKBTM方法.同常規(guī)LSKTM相比，LSKBTM不但具備相當(dāng)?shù)姆囱莩上窬群偷諗克俣?，還大幅度降低了計(jì)算成本.不同于深度域的束偏移方法，本文方法中平面波的合成/分解僅與成像點(diǎn)的速度有關(guān)，因此其對(duì)偏移速度的敏感性同常規(guī)時(shí)間偏移完全相同.雖然本文是通過常規(guī)直射線走時(shí)算法進(jìn)行推導(dǎo)，但很容易引入彎曲射線、各向異性等因素來提高走時(shí)計(jì)算精度.此外，本文方法還可以發(fā)展為共炮檢距算法，從而利用地震數(shù)據(jù)在炮點(diǎn)方向的冗余度進(jìn)一步提高計(jì)算效率.