基于射線理論點(diǎn)擴(kuò)散函數(shù)的成像域最小二乘偏移

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

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

0 引言

作為地震數(shù)據(jù)處理中的關(guān)鍵技術(shù)環(huán)節(jié)之一，地震成像技術(shù)在日趨復(fù)雜化和精細(xì)化的油氣勘探過程中發(fā)揮著重要作用.理論上來說，地震成像方法不僅需要具備精確構(gòu)造成像的能力，還需要具備巖性成像能力以保證地震反演和儲(chǔ)層預(yù)測的精度(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).此外，復(fù)雜的地下介質(zhì)構(gòu)造和非規(guī)則的地震數(shù)據(jù)采集，還會(huì)導(dǎo)致偏移假象和非均勻的成像照明，嚴(yán)重影響成像振幅的可靠性.

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

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

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

1 方法原理

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

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

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

×G(x0,xs,ω)dx0，

(1)

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

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

(2)

其中，m′(x)代表地下成像點(diǎn)x處的成像結(jié)果，符號(hào)*代表復(fù)共軛運(yùn)算.

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

(3)

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

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

(4)

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

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

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

(5)

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

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矩陣的能量主要集中在主對角線附近，因此可以只保留成像點(diǎn)x(PSF中心點(diǎn))附近滿足x0=x+Δx，|Δx|

(8)

其中，H(x,x+Δx)是中心點(diǎn)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)即為射線理論P(yáng)SF直接算法公式，其計(jì)算過程同Kirchhoff偏移非常接近，區(qū)別主要在于走時(shí)項(xiàng)的計(jì)算方式.由于計(jì)算過程相互獨(dú)立，上述算法能夠適應(yīng)任意的PSF場空間采樣，其計(jì)算成本正比于PSF場的總采樣點(diǎn)數(shù)，耗時(shí)至少在一次常規(guī)Kirchhoff偏移以上.

1.2 PSF快速算法

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

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

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

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

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

(10)

其中，ps和pr分別為中心點(diǎn)x處的震源和接收點(diǎn)射線參數(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處的中心點(diǎn)射線參數(shù).

(13)

(14)

(15)

其中，x′r和x′s分別為對應(yīng)的震源和接收點(diǎn).

1.3 算法效率對比

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

(16)

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

(17)

1.4 成像域最小二乘偏移

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

(18)

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

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

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

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è)了三個(gè)反射系數(shù)為1的水平反射界面.通過計(jì)算反射時(shí)距曲線并同主頻為25 Hz的雷克子波進(jìn)行褶積合成了120炮記錄，每炮121道，炮間隔和道間隔為均為30 m.圖2展示了該數(shù)據(jù)的Kirchhoff深度偏移結(jié)果，其中沿其反射界面提取的振幅曲線如圖3所示，可以看到雖然常規(guī)偏移能夠?qū)Ψ瓷浣缑孢M(jìn)行正確歸位，但由于受到幾何擴(kuò)散、非均勻數(shù)據(jù)覆蓋等因素的影響，成像振幅嚴(yán)重失真.

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

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

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

圖4 不同方法求取的PSF場

圖5 中心點(diǎn)在不同空間位置的局部PSF對比

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

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

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

2.2 Marmousi模型試算

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

圖9 Marmousi模型速度場

圖10 Marmousi模型PSF場

圖11 中心點(diǎn)在不同空間位置的局部PSF對比

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

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

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

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

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

圖14 均方根速度場

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

圖15 PSF場

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

圖17 成像剖面局部對比

圖18 成像頻譜對比

3 結(jié)論

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