岳玉波，秦寧，楊哲，陳祥忠，徐云貴，曹衛(wèi)平，唐靜

1 頁巖油氣富集機理與有效開發(fā)國家重點實驗室，北京 102206 2 西南石油大學地球科學與技術(shù)學院，成都 610500 3 中國石化勝利油田物探研究院，東營 257022 4 中國石油勘探院西北分院，蘭州 730020 5 山東理工大學資源與環(huán)境工程學院，淄博 255000

0 引言

作為地震數(shù)據(jù)處理中的關(guān)鍵技術(shù)環(huán)節(jié)之一，地震成像技術(shù)在日趨復雜化和精細化的油氣勘探過程中發(fā)揮著重要作用.理論上來說，地震成像方法不僅需要具備精確構(gòu)造成像的能力，還需要具備巖性成像能力以保證地震反演和儲層預測的精度(Claerbout, 1971; Bleistein et al., 2001).然而，傳統(tǒng)的地震成像方法只是線性正演算子的共軛轉(zhuǎn)置，對于有限觀測系統(tǒng)采集的帶限地震數(shù)據(jù)，只能得到模糊的構(gòu)造成像結(jié)果(Tarantola, 1984; Yu et al., 2006; 岳玉波等, 2012; 劉玉敏等, 2021).此外，復雜的地下介質(zhì)構(gòu)造和非規(guī)則的地震數(shù)據(jù)采集，還會導致偏移假象和非均勻的成像照明，嚴重影響成像振幅的可靠性.

基于線性地震反演理論的最小二乘偏移技術(shù)(Least-Squares Migration, LSM)是解決上述問題的有效途徑(Tarantola, 1984; Nemeth et al., 1999).該技術(shù)將地震偏移成像視為線性反問題并利用最優(yōu)化方法進行迭代求解，理論上能夠消除非規(guī)則采集、帶限子波等因素對成像結(jié)果造成的不良影響，提高成像分辨率和振幅保真度.然而，該技術(shù)的理論優(yōu)勢并沒有充分轉(zhuǎn)化成為實際的應用效果，其主要原因之一在于計算成本過于高昂.經(jīng)典的數(shù)據(jù)域LSM(Nemeth et al., 1999; Duquet et al., 2000; Dai et al., 2011; Yue et al., 2019, 2021a,b; 楊宏偉等, 2022)通過反演迭代使模擬數(shù)據(jù)逐步逼近實際地震數(shù)據(jù)，進而實現(xiàn)地下反射率的準確估計.該方法每次迭代都需要一次正演和偏移運算，整體反演過程所需的計算成本往往在常規(guī)偏移的幾十倍以上，嚴重制約了該方法的推廣應用前景.成像域迭代算法也是目前常用的LSM實現(xiàn)方法，該方法假定Hessian矩陣可以描述地震成像系統(tǒng)對地下介質(zhì)的模糊效應(Jansson, 1997; Lecomte, 2008; Jensen et al., 2021)，在通過反演迭代校正Hessian矩陣的模糊效應后便可獲得地下真實的反射率.成像域LSM的核心是Hessian矩陣的求取，然而由于Hessian矩陣的規(guī)模十分龐大，其計算和存儲過程依然是一項巨大的挑戰(zhàn)(Valenciano et al., 2009; Tang, 2009; Jiang and Zhang, 2019).

Hessian矩陣是一個主對角占優(yōu)的帶狀矩陣(Ren et al., 2011; 任浩然等，2013)，因此可以在不失精度的前提下，通過保留主對角線附近的部分非對角元素來降低Hessian矩陣的計算和存儲成本，這種局部化的Hessian矩陣也稱作點擴散函數(shù)(Point Spread Function, PSF).Fletcher等(2016)、Osorio等(2021)、段偉國等(2022)、Mao等(2022)、陳生昌等(2022)提出通過串聯(lián)的正演+偏移運算來獲取地下稀疏空間位置的PSF場，然后通過插值構(gòu)建任意空間位置的PSF，該過程的計算成本約為兩次常規(guī)偏移，但在實際應用中卻存在如下問題.PSF是同地震速度、觀測系統(tǒng)等因素相關(guān)的非平穩(wěn)空間信號，當?shù)叵陆橘|(zhì)速度變化平緩、地震觀測系統(tǒng)較為規(guī)則時，使用較大的PSF采樣間隔即可滿足插值精度；當?shù)叵陆橘|(zhì)速度變化劇烈、觀測系統(tǒng)非規(guī)則性很強時，則需要使用較為密集的PSF空間采樣，但密集的空間采樣往往會導致相鄰PSF間的交叉串擾，因此在實際應用中往往需要多次計算PSF場，使得計算成本大幅提高.Jiang和Zhang(2019)、Xu等(2022)提出利用地震射線理論求解格林函數(shù)，并以此為基礎(chǔ)進行PSF的直接解析計算.該算法的計算實現(xiàn)過程類似于Kirchhoff偏移，能夠適應任意的PSF場空間采樣，但計算成本依然至少在一次常規(guī)Kirchhoff偏移以上.

本文在上述射線理論PSF直接算法的基礎(chǔ)上，提出了一種更為高效的PSF快速算法.該算法利用地震波走時一階近似對PSF的計算過程進行優(yōu)化，大幅提升了計算效率.以此為基礎(chǔ)，本文進一步發(fā)展了基于PSF的成像域LSM，能夠高效、靈活地實現(xiàn)地下反射率的迭代更新，獲得分辨率更高、照明更加均衡的反演成像結(jié)果.文中通過模型和實際數(shù)據(jù)的試算對PSF快速算法的計算精度和效率以及成像域LSM的應用效果進行了驗證.

1 方法原理

本節(jié)首先對PSF直接算法的基本原理進行簡要介紹，接下來給出本文快速算法的推導實現(xiàn)過程并對兩種算法的計算效率進行對比分析，最后基于PSF構(gòu)建成像域LSM迭代反演流程.

1.1 基于射線理論的PSF直接算法

基于地震散射理論(Tarantola, 1984)，線性Born正演算子可以表示為：

×G(x0,xs,ω)dx0，

(1)

其中，xs和xr分別代表震源和接收點的空間位置，F(xiàn)(ω)是震源子波函數(shù)，m(x0)是散射點x0處的散射強度(Hu et al., 2016)，V代表地下散射點的集合，G(x,x′,ω)是震源位于x′、觀測點位于x的格林函數(shù).通過求取式(1)的共軛轉(zhuǎn)置，可以得到相應的共軛偏移算子：

×G*(x,xs,ω)d(xr,xs,ω)dxrdxs,

(2)

其中，m′(x)代表地下成像點x處的成像結(jié)果，符號*代表復共軛運算.

將式(1)代入式(2)，可以得到成像結(jié)果m′(x)同地下真實散射強度m(x0)之間的關(guān)系：

(3)

其中，H(x,x0)為Hessian矩陣，具有如下的形式：

×G(x0,xr,ω)G(x0,xs,ω)dxrdxs，

(4)

Hessian矩陣描述了復雜速度、非規(guī)則地震采集、帶限子波等因素造成的地震成像系統(tǒng)對地下介質(zhì)的模糊效應，因此常規(guī)偏移只能得到模糊化的構(gòu)造成像結(jié)果.

直接采用波動方程的數(shù)值解計算式(4)中的格林函數(shù)需要耗費高昂的計算和存儲成本.為提高計算效率，可以采用高效、靈活的地震射線理論進行格林函數(shù)的計算：

G(x,x′,ω)=a(x,x′)exp[iωt(x,x′)]，

(5)

其中，a(x,x′)和t(x,x′)分別為地震波的振幅和走時，深度域的走時和振幅信息可以分別通過運動學和動力學射線追蹤求取，在時間域計算則可以使用Zhang等(2000)提出的簡化算法.將式(5)代入式(2)即可得到常規(guī)Kirchhoff積分偏移公式.將式(5)代入式(4)，并應用傅里葉逆變換可得：

H(x,x0)=?a*(x,xr)a*(x,xs)a(x0,xr)a(x0,xs)

-t(x0,xr))dxrdxs，

(6)

其中：

(7)

式(7)為地震子波的自相關(guān)函數(shù).考慮到Hessian矩陣的能量主要集中在主對角線附近，因此可以只保留成像點x(PSF中心點)附近滿足x0=x+Δx，|Δx|

(8)

其中，H(x,x+Δx)是中心點x處的局部化Hessian矩陣，也稱作PSF(Jansson, 1997; Lecomte, 2008)，具有如下的形式：

H(x,x+Δx)=?a*(x,xr)a*(x,xs)a(x+Δx,xr)a(x+Δx,xs)

(9)

式(9)即為射線理論PSF直接算法公式，其計算過程同Kirchhoff偏移非常接近，區(qū)別主要在于走時項的計算方式.由于計算過程相互獨立，上述算法能夠適應任意的PSF場空間采樣，其計算成本正比于PSF場的總采樣點數(shù)，耗時至少在一次常規(guī)Kirchhoff偏移以上.

1.2 PSF快速算法

為進一步提高PSF的計算效率，本文在上述直接算法的基礎(chǔ)上，利用地震波走時一階近似對PSF的計算過程進行優(yōu)化，發(fā)展了如下的快速算法.

由于PSF的有效樣點x+Δx分布在中心點x附近，那么可以將x+Δx處的地震波走時用x處走時的一階泰勒展開來近似表示(如圖1所示)：

圖1 地震射線參數(shù)與走時近似

t(xs,x+Δx)≈t(xs,x)+ps·Δx,t(xr,x+Δx)

≈t(xr,x)+pr·Δx，

(10)

其中，ps和pr分別為中心點x處的震源和接收點射線參數(shù).將式(10)代入式(9)，并假定：

a(x+Δx,xs)≈a(x,xs),a(x+Δx,xr)≈a(x,xr)，

(11)

那么式(9)可以近似表示為：

×dxrdxs，

(12)

其中，pm=ps+pr是x處的中心點射線參數(shù).

(13)

(14)

(15)

其中，x′r和x′s分別為對應的震源和接收點.

1.3 算法效率對比

接下來對比分析兩種PSF算法的計算效率.假設(shè)中心點x處PSF的采樣點數(shù)為Npsf，計算孔徑內(nèi)的地震記錄道數(shù)為Ntr，那么PSF直接算法的計算成本可以估算為：

(16)

Cfast=Ntr×Oamp+Npsf×Np×Oker，

(17)

1.4 成像域最小二乘偏移

在求取PSF場后，可以通過線性插值高效計算地下任意空間位置的PSF.以此為基礎(chǔ)，我們構(gòu)建了如下的成像域LSM目標函數(shù)(Osorio et al., 2021)：

(18)

其中，m為待求的地下反射率，H為求取的PSF場，m′為Kirchhoff偏移結(jié)果，J(m)是以L2模定義的數(shù)據(jù)匹配項，R(m)是用來保證反演過程穩(wěn)定的正則化項，本文在此使用的是Tikhonov正則化(Nemeth et al., 1999)，μ是正則化權(quán)系數(shù)，一般需要多次試驗來確定.式(18)的求解過程一般通過共軛梯度法等梯度類迭代算法進行實現(xiàn)，其計算成本相比于m′的計算過程可以忽略不計.因此，以PSF快速算法為基礎(chǔ)發(fā)展的成像域LSM具備極高的計算效率，所需的計算成本(包括PSF計算和反演迭代過程)遠低于一次常規(guī)Kirchhoff偏移(通常在20%以內(nèi)).

2 模型及實際資料試算

本節(jié)首先采用水平層狀模型和Marmousi模型進行測試，對比兩種PSF算法的計算精度和效率，并驗證成像域LSM的正確性和有效性；接下來采用某探區(qū)二維數(shù)據(jù)進行測試，檢驗成像域LSM對于實際地震數(shù)據(jù)的應用效果.

2.1 層狀模型試算

層狀模型網(wǎng)格大小為480×375，縱橫向采樣間隔分別為10 m和8 m.假設(shè)模型速度為常速2000 m·s-1，在深度為0.8 km、1.6 km、2.4 km處分別布設(shè)了三個反射系數(shù)為1的水平反射界面.通過計算反射時距曲線并同主頻為25 Hz的雷克子波進行褶積合成了120炮記錄，每炮121道，炮間隔和道間隔為均為30 m.圖2展示了該數(shù)據(jù)的Kirchhoff深度偏移結(jié)果，其中沿其反射界面提取的振幅曲線如圖3所示，可以看到雖然常規(guī)偏移能夠?qū)Ψ瓷浣缑孢M行正確歸位，但由于受到幾何擴散、非均勻數(shù)據(jù)覆蓋等因素的影響，成像振幅嚴重失真.

圖2 Kirchhoff深度偏移結(jié)果

圖3 沿圖2反射界面提取的振幅曲線

分別應用直接算法和快速算法進行PSF場的計算，其結(jié)果如圖4a和圖4b所示.圖5a—d進一步放大對比了中心點在x=2.1 km,z=1.3 km處和x=3.3 km,z=1.9 km處的局部PSF計算結(jié)果，可以看到不論對于整體PSF場還是對于單點PSF，兩種算法的計算結(jié)果都非常接近，兩者的相對誤差僅為4%左右，且主要集中在PSF的旁瓣，并不會對LSM反演產(chǎn)生影響.統(tǒng)計對比兩種PSF算法的計算時間，其中直接算法的耗時約為251.3 s，高效算法的耗時為24.7 s(不足直接算法的1/10)，由此可見本文提出PSF快速算法可以在保證計算精度的同時大幅提升計算效率.利用常規(guī)偏移剖面和PSF場作為輸入進行成像域LSM處理，20次迭代后的歸一化數(shù)據(jù)殘差約為3%，最終的反演成像結(jié)果如圖6所示(圖中 “振鈴”狀的成像旁瓣是拓展成像頻帶導致的吉布斯現(xiàn)象).為更好地評估其效果，圖7展示了沿反射界面提取的振幅曲線，圖8對比了反演前后的成像波數(shù)譜，可以看到成像域LSM取得了預期的應用效果，其不但準確恢復了界面的真實反射系數(shù)，還明顯拓寬了成像波數(shù)，提高了成像分辨率.

圖4 不同方法求取的PSF場

圖5 中心點在不同空間位置的局部PSF對比

圖6 成像域LSM結(jié)果

圖7 沿圖6反射界面提取的振幅曲線

圖8 成像波數(shù)譜對比

2.2 Marmousi模型試算

Marmousi模型速度場如圖9所示，網(wǎng)格大小為737×750，縱橫向采樣間隔分別為10 m和4 m.利用有限差分法正演了240炮地震記錄，炮點位于地表，初始炮點位于x=2.5 km處.炮間隔為20 m，每炮120道單邊接收，道間隔也為20 m，震源函數(shù)為主頻20 Hz的雷克子波.圖10a和圖10b分別展示了利用直接算法和快速算法求取的PSF場，其中中心點位于x=3.0 km,z=0.8 km處和x=5.5 km,z=1.3 km處的局部PSF對比結(jié)果如圖11所示.可以看到即便對于復雜的Marmousi模型，快速算法也可以得到媲美直接算法的計算結(jié)果，但其計算效率(耗時81.6 s)相比直接算法(耗時493.2 s)具有明顯的優(yōu)勢.

圖9 Marmousi模型速度場

圖10 Marmousi模型PSF場

圖11 中心點在不同空間位置的局部PSF對比

接下來對正演記錄進行Kirchhoff深度偏移處理，所得到的成像結(jié)果如圖12a所示.可以看到雖然其較為準確地恢復了模型的主要構(gòu)造信息，但成像分辨率較差、成像照明不均.利用PSF場對偏移剖面進行成像域LSM處理，25次迭代后的歸一化數(shù)據(jù)殘差約為10%，最終的反演成像結(jié)果如圖12b所示.可以看到其不但明顯提升了成像分辨率，還有效改善了地下成像照明.圖13提取了x=3.1 km處的成像道波形(紅色曲線)并同真實反射率(藍色曲線)進行對比，可以看到常規(guī)偏移結(jié)果同真實反射率存在很大的偏差(圖13a)，成像波形和振幅嚴重失真，而成像域LSM(圖13b)則大幅提升了成像分辨率和振幅保真性.

圖12 Marmousi模型成像結(jié)果對比

圖13 x=3.1 km處成像道波形(紅色曲線)同理論反射率(藍色曲線)對比

2.3 實際數(shù)據(jù)試算

實際數(shù)據(jù)的偏移測試在時間域進行，其均方根速度場如圖14所示，網(wǎng)格大小為800×750，縱橫向采樣間隔分別為4 ms和12.5 m.將主頻為25 Hz的雷克子波作為震源信號進行PSF場的計算，所得到的結(jié)果如圖15所示.可以看到PSF場的能量中間強、兩側(cè)弱，同地震數(shù)據(jù)的空間覆蓋分布具有很好的對應性.

圖14 均方根速度場

圖16a展示了該數(shù)據(jù)的Kirchhoff時間偏移結(jié)果，可以看到其能量分布特點同圖15所示的PSF場基本一致.利用PSF場對偏移剖面進行成像域LSM處理，20次迭代后的反演成像結(jié)果如圖16b所示.可以看到成像域LSM有效補償了剖面兩側(cè)的成像照明，成像分辨率也得到了明顯的提升.在圖17所示的局部成像對比以及圖18所示的成像頻譜對比中可以看到，成像域LSM明顯拓展了成像頻帶，不但具有更高的薄層分辨能力(圖17紅色箭頭)，其補償?shù)牡皖l信息還有效增強了成像連續(xù)性(圖17紅色方框).

圖15 PSF場

圖16 成像結(jié)果對比

圖17 成像剖面局部對比

圖18 成像頻譜對比

3 結(jié)論

LSM是當前地震成像技術(shù)重要的研究和發(fā)展方向，但龐大的計算成本嚴重制約了其的應用前景.本文基于射線理論格林函數(shù)和地震波走時一階近似，提出了一種快速的PSF計算方法，不但可以獲得媲美原有算法的PSF計算結(jié)果，還具有明顯的計算效率優(yōu)勢.以此為基礎(chǔ)發(fā)展的成像域LSM能以明顯低于常規(guī)偏移的計算成本實現(xiàn)地震成像剖面的迭代更新，獲得高分辨率、高保真的反演成像結(jié)果，有望成為一種切實可行的LSM實施方案.接下來，我們會將本文工作進一步拓展到三維，并研究使用L1、TV等正則化算法以及合理的預條件算子，進一步優(yōu)化本文LSM的成像效果和計算效率.