最小二乘Kirchhoff射線束疊前時間偏移

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

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

0 引言

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

眾所周知，常規的偏移方法只是地震波場正演算子的共軛轉置(Tarantola, 1984; Claerbout, 1985)，對于有限觀測系統所采集的帶限地震數據，其只能生成模糊的地下構造成像結果.此外，復雜的地下介質構造和非規則的地震數據空間采樣，還會導致偏移假象和非規則的成像照明，嚴重影響成像振幅信息的可靠性.基于線性化地震反演理論的最小二乘偏移(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).其將地震偏移視為線性反問題，并利用最優化方法求取能夠準確模擬采集地震數據的成像剖面作為地下真實反射率，在理論上能夠消除非規則采集、帶限子波等因素對成像結果的不利影響.但是，LSM每次迭代均需一次正演和偏移運算，完整的LSM反演過程往往需要幾十倍于常規偏移的計算量，龐大的計算成本嚴重影響了該技術的應用前景(張宇，2018).

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

1 方法原理

1.1 LSM基本原理

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

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

(1)

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

1.2 KFM/KTM算子

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

D(rd,rs,ω)=

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

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

(3)

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

(4)

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

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

(5)

同時，求取式(5)的共軛轉置可以得到相應的KTM算子：

(6)

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

1.3 KBFM算子

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

(7)

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

(8)

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

(9)

此時，式(8)可以轉化為

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

(10)

如果將對應接收點射線參數為pd的散射點的集合定義為Ω0(pd)，那么式(10)可以表示為

(11)

其中：

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

(12)

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

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

(13)

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

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

(14)

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

(15)

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

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

(2)對于每一個射線束中心L，進行地下散射點的循環，計算式(14)中的走時、振幅等參數，并以此為基礎將地下反射率映射轉換到τ-p域.

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

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

1.4 KBTM算子

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

+T(L)),

(16)

(17)

(18)

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

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

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

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

2 模型及實際資料試算

本節首先應用Marmousi模型進行LSKTM和LSKBTM的測試，證明兩者具有非常接近的成像精度，并驗證LSKBTM在計算效率上的優勢.接下來，應用二維和三維探區地震數據測試LSKBTM的實際應用效果.

2.1 Marmousi模型

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

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

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

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

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

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

圖7 歸一化目標函數收斂曲線其中藍色曲線代表LSKTM，綠色曲線代表LSKBTM.Fig.7 Normalized misfits versus number of iterationsBlue is LSKTM, and green is LSKBTM.

2.2 二維實際數據

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

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

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

圖9 KBTM成像結果Fig.9 Imaging result using KBTM

圖10 LSKBTM成像結果Fig.10 Imaging result using LSKBTM

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

2.3 三維實際數據

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

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

圖13 KBTM成像結果Fig.13 Image result of KBTM

圖14 LSKBTM成像結果Fig.14 Image result of LSKBTM

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

3 結論

本文基于地震數據的局部平面波分解/合成策略，在時間域KFM/KTM的基礎上，實現了高效的KBFM/KBTM算子，并以此為基礎發展了快速的LSKBTM方法.同常規LSKTM相比，LSKBTM不但具備相當的反演成像精度和迭代收斂速度，還大幅度降低了計算成本.不同于深度域的束偏移方法，本文方法中平面波的合成/分解僅與成像點的速度有關，因此其對偏移速度的敏感性同常規時間偏移完全相同.雖然本文是通過常規直射線走時算法進行推導，但很容易引入彎曲射線、各向異性等因素來提高走時計算精度.此外，本文方法還可以發展為共炮檢距算法，從而利用地震數據在炮點方向的冗余度進一步提高計算效率.