高效的成像域最小二乘逆時(shí)偏移

陳生昌， 李代光， 金成玫

浙江大學(xué)地球科學(xué)學(xué)院, 杭州 310027

0 引言

地震數(shù)據(jù)偏移成像是得到地下構(gòu)造圖像的主要方法技術(shù)，特別是地下復(fù)雜構(gòu)造區(qū)的構(gòu)造成像(Gray et al., 2001；Yilmaz, 2001；Etgen et al.，2009；Robein, 2010; Leveille et al.，2011).隨著國內(nèi)“兩寬一高”地震數(shù)據(jù)采集技術(shù)的廣泛應(yīng)用和油氣勘探開發(fā)目標(biāo)的日趨復(fù)雜化和精細(xì)化，對(duì)于地震數(shù)據(jù)偏移成像方法技術(shù)要求能充分地利用地震數(shù)據(jù)的波形信息，得到的偏移成像結(jié)果不僅要滿足復(fù)雜區(qū)域構(gòu)造成像的要求，還要具有高分辨和保真性，使之能精細(xì)地反映地下反射面巖性變化的特征，以滿足地震數(shù)據(jù)巖性處理解釋，特別是復(fù)合型復(fù)雜油氣藏和非常規(guī)油氣藏地震數(shù)據(jù)巖性處理解釋的需要.逆時(shí)偏移是一種基于全波方程的偏移方法(Baysal et al., 1983; MeMECHAN, 1983; Whitmore, 1983)，被認(rèn)為是當(dāng)前在理論上最為先進(jìn)的偏移方法(Gray et al., 2001).相比于通常僅考慮單一射線路徑走時(shí)的Kirchhoff偏移(Schneider, 1978; Bleistein et al., 2005)，逆時(shí)偏移可有效地考慮地震波傳播過程中的多路徑與多走時(shí).相比于基于地震波方向分解的單程波偏移存在的大角度波(相對(duì)于深度方向的大角度波)限制(Ristow and Ruhl, 1994; Zhang et al., 2005)，逆時(shí)偏移無地震波傳播角度限制.

給定具有地震波運(yùn)動(dòng)學(xué)準(zhǔn)確的光滑偏移速度模型，常規(guī)逆時(shí)偏移可視為是一種利用一次波信息對(duì)有關(guān)地下反射率的線性反演問題的粗糙求解(Claerbout, 1985; Schuster, 2017).利用線性反演理論建立的有關(guān)地下反射率反演的最小二乘逆時(shí)偏移是對(duì)常規(guī)逆時(shí)偏移的改進(jìn)，它可有效地消除地震子波、地震數(shù)據(jù)觀測(cè)系統(tǒng)、成像目標(biāo)上覆介質(zhì)空間變化和地震波幾何擴(kuò)散對(duì)常規(guī)偏移成像的影響，改善地震波弱照明區(qū)和深部的偏移成像效果，得到具有高保真性和高分辨率的偏移成像結(jié)果(Tarantola, 1984；Yao and Jakubowicz, 2012; 黃建平等, 2013; Yang and Zhang, 2019).由于最小二乘逆時(shí)偏移涉及偏移結(jié)果的成像域到地震數(shù)據(jù)的時(shí)空域之間線性映射矩陣L的求逆，對(duì)于三維炮集地震數(shù)據(jù)的最小二乘逆時(shí)偏移，矩陣L的行數(shù)M和列數(shù)N，分別有M=Ny×Nx×Nt和列數(shù)N=Ny×Nx×Nz(其中，Ny、Nx和Nz分別為成像域的Y、X和Z方向的網(wǎng)格點(diǎn)數(shù)，Nt為地震道采樣點(diǎn)數(shù))，假設(shè)Ny=1000，Nx=1000，Nz=1000，Nt=10000，則L為1010×109的超大規(guī)模矩陣，致使利用矩陣L求逆的方法進(jìn)行最小二乘逆時(shí)偏移在現(xiàn)有計(jì)算機(jī)條件下難以實(shí)現(xiàn).因此，當(dāng)前對(duì)于最小二乘逆時(shí)偏移的求解有兩類方法，一是數(shù)據(jù)域方法(Schuster, 1993; Fletcher et al., 2012; Zhang et al., 2015; Zhao et al., 2015; Xue et al., 2016; Fletcher and Cavalca, 2018)；二是成像域方法(Yu et al., 2006; Aoki and Schuster, 2009; Tang, 2009; Fletcher et al., 2016).在數(shù)據(jù)域最小二乘逆時(shí)偏移方法中，以觀測(cè)的炮集地震記錄與數(shù)值計(jì)算的炮集地震記錄間的最佳擬合為準(zhǔn)則，采樣Landweber迭代方法(陳生昌和周華敏，2018)避免最小二乘逆時(shí)偏移中大規(guī)模矩陣求逆的計(jì)算復(fù)雜度.但在數(shù)據(jù)域最小二乘逆時(shí)偏移的迭代中需要進(jìn)行炮集地震數(shù)據(jù)的偏移和利用偏移成像結(jié)果的反偏移，所以數(shù)據(jù)域最小二乘逆時(shí)偏移方法的計(jì)算量巨大，約為常規(guī)逆時(shí)偏移方法計(jì)算量的2×Nite(Nite為迭代次數(shù))倍.在成像域最小二乘逆時(shí)偏移方法中，利用常規(guī)逆時(shí)偏移結(jié)果可構(gòu)建常規(guī)逆時(shí)偏移成像結(jié)果與最小二乘逆時(shí)偏移成像結(jié)果之間基于Hessian矩陣H(Hessian矩陣H為最小二乘逆時(shí)偏移的目標(biāo)函數(shù)關(guān)于地下反射率的二階導(dǎo)數(shù)，也是矩陣L的伴隨與L的乘積，即H=L*L，H中的列向量可視為成像點(diǎn)的點(diǎn)擴(kuò)展函數(shù))的空變線性積分關(guān)系，使成像域成為了最小二乘逆時(shí)偏移的計(jì)算空間.對(duì)于上述Ny=1000，Nx=1000，Nz=1000的三維最小二乘逆時(shí)偏移，矩陣H為109×109的超大規(guī)模矩陣，因此在現(xiàn)有計(jì)算機(jī)條件下也同樣難以利用矩陣H求逆的方法實(shí)現(xiàn)成像域的三維最小二乘逆時(shí)偏移.此外，由矩陣H的表達(dá)式可知，形成矩陣H的計(jì)算量對(duì)于三維最小二乘偏移也是十分巨大的，因此在成像域利用Landweber迭代方法求解也同樣存在巨大的計(jì)算量問題.鑒于成像域最小二乘逆時(shí)偏移的計(jì)算量與偏移成像區(qū)域大小有關(guān)，目前該方法多用于區(qū)域較小的目標(biāo)區(qū)最小二乘逆時(shí)偏移(Valenciano et al., 2006; Tang and Biondi, 2013; Zhao and Sen, 2018).

成像域最小二乘逆時(shí)偏移中成像點(diǎn)的點(diǎn)擴(kuò)散函數(shù)的局部分布形態(tài)主要由地震數(shù)據(jù)觀測(cè)系統(tǒng)和用于偏移的速度模型決定，而地震數(shù)據(jù)觀測(cè)的多次覆蓋特點(diǎn)和偏移速度模型的空間光滑性，所以成像域的點(diǎn)擴(kuò)散函數(shù)具有一定的相似性.針對(duì)上述存在的問題和認(rèn)識(shí)，本文把常規(guī)逆時(shí)偏移成像結(jié)果與最小二乘逆時(shí)偏移成像結(jié)果之間在成像域的空變積分關(guān)系近似為局部區(qū)域的空不變褶積關(guān)系，提出基于局域空不變反褶積的成像域最小二乘逆時(shí)偏移方法.通過利用求解空不變反褶積的波數(shù)域偽廣義逆方法(陳生昌和王芳,1993；陳生昌和肖鵬飛，2007)，可得到一種穩(wěn)定高效的成像域最小二乘逆時(shí)偏移方法.把本文提出的成像域高效最小二乘逆時(shí)偏移方法應(yīng)用于模擬和實(shí)際地震數(shù)據(jù)均取得了理想的偏移成像結(jié)果.

1 最小二乘逆時(shí)偏移簡(jiǎn)介

給定具有地震波運(yùn)動(dòng)學(xué)準(zhǔn)確的光滑偏移速度模型，作為一種線性反演的地震數(shù)據(jù)偏移所對(duì)應(yīng)的地震數(shù)據(jù)線性正演方程可寫為

d(xr,xs,ω)=L(xr,x,xs,ω)m(x),

(1)

式中，d(xr,xs,ω)為炮道集的頻率域一次波地震數(shù)據(jù)；L(xr,x,xs,ω)為偏移成像結(jié)果到地震數(shù)據(jù)的線性映射算子；m(x)為偏移成像結(jié)果；xr為接收點(diǎn)坐標(biāo)；xs為震源點(diǎn)坐標(biāo)；ω為頻率；x為地下成像點(diǎn)坐標(biāo).關(guān)于方程(1)中線性算子L(xr,x,xs,ω)的具體表達(dá)式，當(dāng)前有三種不同的版本：一是平面波無窮大平反射面的經(jīng)典反射版本(Claerbout, 1971; Bleistein et al., 2001)；二是基于Born近似的散射版本(Schuster, 2017；Yang and Zhang, 2019)；三是局部平面波局部平反射面的局部反射版本(陳生昌和周華敏，2018).對(duì)于經(jīng)典反射版本的線性算子L(xr,x,xs,ω),其具體表達(dá)式有

L(xr,x,xs,ω)=Gr(xr,x,ω)Gi(x,xs,ω)s(ω)，

(2)

對(duì)于散射版本的線性算子L(xr,x,xs,ω),其具體表達(dá)式有

(3)

對(duì)于局部反射版本的線性算子L(xr,x,xs,ω),其具體表達(dá)式有

(4)

上述式中，vm(x)表示用于偏移的光滑速度模型；Gr(xr,x,ω)表示散射(或反射)傳播的Green函數(shù)；Gi(x,xs,ω)表示入射傳播的Green函數(shù)；s(ω)表示震源子波時(shí)間函數(shù)的Fourier變換.

把方程(1)簡(jiǎn)寫為矩陣方程形式，有

d=Lm，

(5)

式中，d為地震數(shù)據(jù)d(xr,xs,ω)所對(duì)應(yīng)的矩陣；L為線性算子L(xr,x,xs,ω)所對(duì)應(yīng)的矩陣；m為偏移成像結(jié)果m(x)所對(duì)應(yīng)的矩陣.利用線性最小二乘反演方法，由方程(5)可得到有關(guān)m的最小二乘解估計(jì)，即

mls=(L*L)-1L*d,

(6)

式中L*為L(zhǎng)矩陣的伴隨.對(duì)于方程(6)的解表達(dá)，由于矩陣L的超大規(guī)模，致使難以通過直接求矩陣逆的方式得到解mls.針對(duì)這樣的大規(guī)模矩陣方程，在計(jì)算數(shù)學(xué)中常采用近似迭代的方式進(jìn)行求解，如Landweber迭代(Kirsch, 2011).對(duì)于方程(6)，有如下的Landweber迭代格式，

(7)

式中，k為迭代次數(shù)；τ為解的修正步長(zhǎng).由迭代格式(7)可知，其反演運(yùn)算主要是在數(shù)據(jù)域進(jìn)行，因此應(yīng)用迭代格式(7)的最小二乘偏移被稱為數(shù)據(jù)域最小二乘偏移方法.

對(duì)于方程(6)可寫為

mls=H-1mc,

(8)

其中，H=L*L,根據(jù)線性最小二乘反演理論可知，H矩陣為線性反演問題(5)的Hessian矩陣；對(duì)于mc=L*d，L*d表示對(duì)地震數(shù)據(jù)d進(jìn)行常規(guī)逆時(shí)偏移運(yùn)算，mc表示常規(guī)逆時(shí)偏移結(jié)果.

由公式(8)可知，把Hessian矩陣H的逆作用于常規(guī)逆時(shí)偏移結(jié)果就可得到最小二乘逆時(shí)偏移結(jié)果.由于公式(8)的運(yùn)算僅在成像域進(jìn)行，通常把下述三步運(yùn)算組成的最小二乘逆時(shí)偏移稱為成像域最小二乘逆時(shí)偏移方法，即

(9)

如果偏移成像的范圍比較大，直接按公式(9)進(jìn)行最小二乘逆時(shí)偏移存在形成Hessian矩陣H和求取逆矩陣H-1的兩大計(jì)算困難，它們涉及的計(jì)算量巨大.因此，由公式(9)組成的成像域最小二乘逆時(shí)偏移比較適合偏移成像范圍比較小的區(qū)域，如面向目標(biāo)區(qū)的最小二乘逆時(shí)偏移(Tang and Biondi, 2013; Zhao and Sen, 2018).為了進(jìn)行大區(qū)域三維地震數(shù)據(jù)的成像域最小二乘逆時(shí)偏移，我們需要研究可實(shí)現(xiàn)公式(9)所表示的成像域最小二乘逆時(shí)偏移的近似算法.

2 高效的成像域最小二乘逆時(shí)偏移方法

公式(8)所表示的反演求解形式有與其對(duì)應(yīng)的正演表示式，即

mc=Hmls.

(10)

在公式(10)中，mc和mls分別表示由常規(guī)逆時(shí)偏移結(jié)果和最小二乘逆時(shí)偏移結(jié)果組成的列矩陣，即mc=[mc1,mc2,…,mcN]T，mls=[mls1,mls2,…,mlsN]T,它們分別是正演表達(dá)式中的數(shù)據(jù)和模型，N代表成像域的網(wǎng)格點(diǎn)數(shù)；H為N×N的方陣.

對(duì)于矩陣方程(9)，如果矩陣mls中僅第i個(gè)元素值等于1，其他的元素值都等于0，則矩陣mc的值等于矩陣H中第i個(gè)列向量的值，即

(11)

因此，H中的列向量值也稱為成像域中成像點(diǎn)的點(diǎn)擴(kuò)散函數(shù)(PSF)值.根據(jù)矩陣H的表達(dá)式，H=L*L，利用公式(2)或(3)或(4)，可得到下述多炮多道的PSF計(jì)算式(利用公式(3))，

(12)

式中，x為PSF在成像域中的計(jì)算點(diǎn)位置；x0代表成像域中的成像點(diǎn)位置；PSF(x,x0)為成像點(diǎn)x0的點(diǎn)擴(kuò)散函數(shù)；Re表示取實(shí)部運(yùn)算.

式(12)所表示的計(jì)算實(shí)質(zhì)是在偏移速度模型vm(x)下對(duì)成像點(diǎn)x0處單位脈沖式散射體的地震波場(chǎng)多炮模擬和多炮逆時(shí)偏移成像，即從炮點(diǎn)xs激發(fā)的震源波場(chǎng)傳播到成像點(diǎn)x0，與該點(diǎn)處的單位脈沖式散射體發(fā)生作用形成散射虛源激發(fā)散射波，散射波場(chǎng)傳播到接收點(diǎn)xr，得到記錄波場(chǎng)，然后多道記錄波場(chǎng)經(jīng)伴隨傳播到PSF計(jì)算點(diǎn)x，同時(shí)從炮點(diǎn)xs激發(fā)的震源波場(chǎng)也經(jīng)伴隨傳播到計(jì)算點(diǎn)x，最后再利用成像條件得到偏移成像結(jié)果.因此，式(12)所得到的PSF(x,x0)即為成像點(diǎn)x0處單位脈沖式散射體的逆時(shí)偏移成像結(jié)果.由地震波的傳播理論和偏移成像理論可知，PSF(x,x0)具有以成像點(diǎn)x0為中心的局部分布特征.該局部分布特征為多成像點(diǎn)的PSF同時(shí)計(jì)算提供了可能.

根據(jù)地震波傳播與成像理論、矩陣H的表達(dá)式和PSF的計(jì)算式(12)，公式(10)對(duì)應(yīng)一個(gè)線性空變積分，因此可把式(10)寫為下述的空變積分形式，即

(13)

式中k(x,y)為與式(10)中矩陣H對(duì)應(yīng)的空變積分核.與式(10)對(duì)應(yīng)的反演轉(zhuǎn)化為式(13)的空變線性積分方程求解.由于積分核k(x,y)的空變，對(duì)于積分方程(12)的求解沒有快速算法.

k(x,y)的空變特征與偏移速度模型vm(x)的空間變化和地震數(shù)據(jù)觀測(cè)系統(tǒng)范圍的有限性、覆蓋的不均勻性密切相關(guān).如果vm(x)為均勻模型、觀測(cè)系統(tǒng)具有無限大的觀測(cè)范圍且均勻覆蓋，則k(x,y)退化為空不變，相應(yīng)的式(13)退化為線性空不變褶積型積分方程，即

(14)

對(duì)于式(14)的線性空不變褶積型積分方程可利用波數(shù)域偽廣義逆方法(陳生昌和王芳，1993；陳生昌和肖鵬飛，2007)進(jìn)行穩(wěn)定高效地求解，即

(15)

由于偏移速度模型vm(x)為空間宏觀變化的光滑速度模型和地震數(shù)據(jù)采集的規(guī)則多次覆蓋特性，我們認(rèn)為光滑模型vm(x)和多次覆蓋下k(x,y)的空間變化是很緩慢的，即k(x,y)在一定的空間范圍內(nèi)可近似為空不變的.對(duì)應(yīng)的在一定空間范圍內(nèi)的各成像點(diǎn)的PSF(x,x0)的變化形態(tài)具有很大的相似性.

圖1為用于成像點(diǎn)PSF計(jì)算的Marmousi模型的光滑偏移速度模型.圖2為9個(gè)不同網(wǎng)格點(diǎn)((340,100)、(370,100)、(400,100)，(340,130)、(370,130)、(400,130)，(340,160)、(370,160)、(400,160))位置成像點(diǎn)的PSF分布圖.計(jì)算PSF的觀測(cè)系統(tǒng)為地面均勻全覆蓋.由圖2所示的PSF分布圖可看出，成像點(diǎn)的PSF分布圖不僅具有很好的以成像點(diǎn)為中心的局部分布特征，而且x和z方向距離30個(gè)網(wǎng)格點(diǎn)(網(wǎng)格點(diǎn)間距dx=dz=10 m)的成像點(diǎn)的PSF還具有很強(qiáng)的相似性，因此可認(rèn)為在30×30網(wǎng)格點(diǎn)的區(qū)域內(nèi)各個(gè)成像點(diǎn)的PSF形態(tài)變化是很小的，也即在該范圍內(nèi)的k(x,y)可近似為空不變的.圖3為圖1 Marmousi模型上間距30×30網(wǎng)格的所有成像點(diǎn)的PSF分布圖(圖3的顯示加了perc=99的增益).由圖3同樣可以看出圖中相鄰成像點(diǎn)PSF分布圖的相似性.

圖1 用于成像點(diǎn)PSF計(jì)算的Marmousi模型的光滑速度模型Fig.1 The smooth velocity model of Marmousi model for the PSF calculation of imaging points

圖2 9個(gè)不同網(wǎng)格點(diǎn)位置成像點(diǎn)的PSF分布圖Fig.2 The PSF distribution of 9 imaging points at different grid-point positions

圖3 圖1速度模型上網(wǎng)格點(diǎn)間距30×30的所有成像點(diǎn)的PSF分布圖Fig.3 The PSF distribution of all imaging points with 30×30 grid-points interval in velocity model of Fig.1

根據(jù)上述有關(guān)k(x,y)的局部分布特性和光滑速度模型、規(guī)則多次覆蓋下k(x,y)在一定空間范圍內(nèi)可近似為空不變的認(rèn)識(shí)，式(13)的空變積分方程可近似為局域空不變褶積的疊加，即

(16)

式中，Lb表示可把k(x,y)近似為空不變的局部區(qū)塊；kLb(x)表示以局部區(qū)塊的中心點(diǎn)xb為成像點(diǎn)的PSF近似的空不變核函數(shù).把式(15)表示的線性空不變褶積型積分方程求解方法，應(yīng)用于式(16)，有

(17)

(18)

(19)

式中，變量β、γ和φ的意義見示意圖4.

圖4 區(qū)塊邊界拼接示意圖Fig.4 The sketch of block boundary tapering

把上述的常規(guī)逆時(shí)偏移、多炮多道的PSF計(jì)算和式(17)的分塊局域空不變褶積型積分方程求解方法結(jié)合起來，就可得到一種穩(wěn)定高效的成像域最小二乘逆時(shí)偏移方法，我們也把這種最小二乘逆時(shí)偏移方法稱之為成像域的局域空不變反褶積方法.常規(guī)數(shù)據(jù)域最小二乘逆時(shí)偏移的計(jì)算量約為常規(guī)逆時(shí)偏移的2×Nite(一般情況下，Nite>10)倍，而本文的最小二乘逆時(shí)偏移的計(jì)算量約為常規(guī)逆時(shí)偏移的2倍，所以本文方法相對(duì)于常規(guī)數(shù)據(jù)域最小二乘逆時(shí)偏移方法約提高計(jì)算效率Nite倍.

在本文提出的成像域最小二乘逆時(shí)偏移方法中，局域區(qū)塊Lb尺寸的大小不僅關(guān)系到方法的計(jì)算量，也關(guān)系到方法的最終偏移成像效果.局域區(qū)塊Lb尺寸大，雖然方法的計(jì)算量小，但對(duì)k(x,y)的空變性考慮不足，會(huì)影響方法的最終偏移成像效果.如果局域區(qū)塊Lb取為整個(gè)偏移成像區(qū)域，雖然方法的計(jì)算量最小，但相當(dāng)于把k(x,y)近似為全成像域空不變(如方程(13)).局域區(qū)塊Lb尺寸小，雖然方法的計(jì)算量大，但能較好地考慮k(x,y)的空變性，有利于最終的偏移成像結(jié)果.如果局域區(qū)塊Lb取為一個(gè)點(diǎn)，則方法就等同于第1節(jié)公式(8)的成像域最小二乘逆時(shí)偏移成像方法.因此，區(qū)塊Lb尺寸大小的選擇原則為，在保證區(qū)塊中心點(diǎn)xb的PSF(x,xb)完整的條件下，使區(qū)塊Lb的尺寸盡量小，即在兼顧計(jì)算量的同時(shí)盡可能考慮k(x,y)的空變性以保證偏移成像效果.PSF(x,xb)的完整也就是褶積算子的完整，這是得到高分辨率反褶積結(jié)果的必要條件之一.

3 數(shù)值試驗(yàn)

為了驗(yàn)證本文提出的把空變的PSF近似為局域空不變PSF的可行性、正確性和局域空不變反褶積成像域最小二乘逆時(shí)偏移方法對(duì)實(shí)際數(shù)據(jù)的有效性，我們進(jìn)行了基于二維Marmousi模型的合成地震數(shù)據(jù)、二維和三維實(shí)際地震數(shù)據(jù)的數(shù)值試驗(yàn).在試驗(yàn)中，要求已得到地震數(shù)據(jù)的常規(guī)逆時(shí)偏移結(jié)果和用于逆時(shí)偏移的光滑偏移速度模型.根據(jù)PSF的局部分布特性，在給定的光滑偏移速度模型上，按照確定的間距規(guī)則地放置單位脈沖式點(diǎn)散射體，點(diǎn)散射體間的間距所確定區(qū)域即認(rèn)為是可把k(x,y)近似為空不變的局部區(qū)塊，點(diǎn)散射體(即成像點(diǎn))位于局部區(qū)塊中心；然后根據(jù)地震數(shù)據(jù)采集的觀測(cè)系統(tǒng)利用公式(12)計(jì)算點(diǎn)散射體(成像點(diǎn))在對(duì)應(yīng)的局部區(qū)塊的PSF分布；最后把常規(guī)逆時(shí)偏移結(jié)果和各個(gè)局部區(qū)塊的PSF分布(對(duì)應(yīng)的Fourier變換譜)應(yīng)用于本文提出的穩(wěn)定高效的成像域高效最小二乘逆時(shí)偏移計(jì)算式(17)，得到最小二乘逆時(shí)偏移成像結(jié)果.由于各個(gè)局部區(qū)塊的PSF僅分布于其對(duì)應(yīng)的區(qū)塊內(nèi)，因此利用公式(12)的一次計(jì)算就可以計(jì)算出成像域內(nèi)各個(gè)局部區(qū)塊的PSF，計(jì)算代價(jià)約為常規(guī)逆時(shí)偏移的2倍.

3.1 Marmousi模型合成地震數(shù)據(jù)

Marmousi模型為二維模型，試驗(yàn)中，用于成像點(diǎn)PSF計(jì)算的局部區(qū)塊尺寸為30×30網(wǎng)格點(diǎn)，即成像點(diǎn)在x方向和z方向的間距均為30個(gè)網(wǎng)格點(diǎn).圖3為利用圖1所示的光滑偏移速度模型和公式(12)計(jì)算得到的全成像域間距30×30的所有成像點(diǎn)的PSF分布圖.由于k(x,y)的本質(zhì)空變特征，所以各個(gè)局部區(qū)塊的成像點(diǎn)PSF還是存在差異的，尤其是深部和邊部的成像點(diǎn)PSF.圖5為Marmousi模型合成數(shù)據(jù)利用圖1所示的光滑偏移速度模型得到的常規(guī)逆時(shí)偏移成像結(jié)果.

圖5 Marmousi模型數(shù)據(jù)的常規(guī)逆時(shí)偏移成像結(jié)果Fig.5 The conventional RTM result of Marmousi model data

圖7的試驗(yàn)結(jié)果表明，即使k(x,y)在理論上是空變的，由于用于偏移和PSF計(jì)算的速度模型的光滑性和多次覆蓋的數(shù)據(jù)采集，在實(shí)踐中把k(x,y)近似為局域空不變，得到的成像域最小二乘逆時(shí)偏移效果是可行的和正確的.圖6的試驗(yàn)結(jié)果表明，用某個(gè)局部區(qū)塊的PSF分布作為全成像域的PSF分布，雖然也能得到不錯(cuò)的成像效果，計(jì)算量也很小，但由于不能考慮k(x,y)的空變特性，我們不建議使用全域空不變反褶積進(jìn)行成像域最小二乘逆時(shí)偏移.

圖6 全域空不變反褶積得到的 Marmousi模型數(shù)據(jù)的成像域最小二乘逆時(shí)偏移成像結(jié)果Fig.6 The IDLSRTM result of Marmousi model data by global spatial-invariant de-convolution

圖7 局域空不變反褶積得到的 Marmousi模型數(shù)據(jù)的成像域最小二乘逆時(shí)偏移成像結(jié)果Fig.7 The IDLSRTM result of Marmousi model data by local spatial-invariant de-convolution

3.2 二維實(shí)際地震數(shù)據(jù)

為了進(jìn)一步驗(yàn)證本文提出的局域空不變反褶積最小二乘逆時(shí)偏移方法在實(shí)際地震數(shù)據(jù)上的效果，我們首先用某地的一條二維線進(jìn)行試驗(yàn).圖8為該線的常規(guī)逆時(shí)偏移成像結(jié)果.試驗(yàn)中計(jì)算各個(gè)局部區(qū)塊中心點(diǎn)PSF分布的速度模型是用于該線逆時(shí)偏移成像的偏移速度模型，觀測(cè)系統(tǒng)采用實(shí)際地震數(shù)據(jù)采集的觀測(cè)系統(tǒng).用于成像點(diǎn)PSF計(jì)算的局域區(qū)塊尺寸為20×40網(wǎng)格點(diǎn)(網(wǎng)格點(diǎn)間距dx=20 m，dz=10 m)，即成像點(diǎn)在x方向和z方向的間距分別為20和40個(gè)網(wǎng)格點(diǎn).圖9為計(jì)算得到的成像域網(wǎng)格點(diǎn)間距20×40的所有成像點(diǎn)的PSF分布圖.

圖8 二維實(shí)際數(shù)據(jù)的常規(guī)逆時(shí)偏移成像結(jié)果Fig.8 The conventional RTM result of 2D field data

圖9 成像域網(wǎng)格點(diǎn)間距20×40的所有成像點(diǎn)的PSF分布圖Fig.9 The PSF distribution of all imaging points with 20×40 grid-points interval in imaging domain

圖10 二維實(shí)際數(shù)據(jù)的成像域最小二乘逆時(shí)偏移成像結(jié)果Fig.10 The IDLSRTM result of 2D field data

3.3 三維實(shí)際地震數(shù)據(jù)

為了驗(yàn)證本文方法對(duì)三維實(shí)際地震數(shù)據(jù)最小二乘逆時(shí)偏移成像的有效性，我們對(duì)某工區(qū)的三維實(shí)際地震數(shù)據(jù)進(jìn)行了試驗(yàn).與上面的二維實(shí)際數(shù)據(jù)試驗(yàn)一樣，我們首先有該實(shí)際數(shù)據(jù)的三維逆時(shí)偏移成像結(jié)果，并利用實(shí)際地震數(shù)據(jù)采集的觀測(cè)系統(tǒng)和逆時(shí)偏移的偏移速度模型計(jì)算三維全成像域各個(gè)局部區(qū)塊中心點(diǎn)的PSF分布.局域區(qū)塊尺寸為20×20×50網(wǎng)格點(diǎn)(網(wǎng)格點(diǎn)間距dx=30 m，dy=30 m，dz=10 m)，即成像點(diǎn)在x方向和y方向的間距均為20個(gè)網(wǎng)格點(diǎn)，在z方向的間距為50個(gè)網(wǎng)格點(diǎn).

圖11為其中一個(gè)局部區(qū)塊中心點(diǎn)的三維PSF分布的xz平面切片圖.由于該P(yáng)SF的成像點(diǎn)位于成像域的中淺部，PSF的完整性比較好.

圖11 三維PSF分布的xz平面切片圖Fig.11 The xz plane slice of 3D PSF distribution

圖12、13分別為縱測(cè)線131的常規(guī)逆時(shí)偏移成像結(jié)果和本文提出的成像域最小二乘逆時(shí)偏移成像結(jié)果，圖14、15分別為縱測(cè)線231的常規(guī)逆時(shí)偏移成像結(jié)果和本文提出的成像域最小二乘逆時(shí)偏移成像結(jié)果，圖16、17分別為橫測(cè)線100的常規(guī)逆時(shí)偏移成像結(jié)果和本文提出的成像域最小二乘逆時(shí)偏移成像結(jié)果.通過比較圖12與圖13、圖14與圖15和圖16與圖17可看出，成像域最小二乘逆時(shí)偏移成像結(jié)果相對(duì)于常規(guī)逆時(shí)偏移成像結(jié)果在空間分辨率方面有很大的提高.

圖12 縱測(cè)線131的常規(guī)逆時(shí)偏移成像結(jié)果Fig.12 The conventional RTM result of In-Line 131

圖13 縱測(cè)線131的成像域最小二乘逆時(shí)偏移成像結(jié)果Fig.13 The IDLSRTM result of In-Line 131

圖14 縱測(cè)線231的常規(guī)逆時(shí)偏移成像結(jié)果Fig.14 The conventional RTM result of In-Line 231

圖15 縱測(cè)線231的成像域最小二乘逆時(shí)偏移成像結(jié)果Fig.15 The IDLSRTM result of In-Line 231

圖16 橫測(cè)線100的常規(guī)逆時(shí)偏移成像結(jié)果Fig.16 The conventional RTM result of Cross-Line 100

圖17 橫測(cè)線100的成像域最小二乘逆時(shí)偏移成像結(jié)果Fig.17 The IDLSRTM result of Cross-Line 100

Marmousi模型合成數(shù)據(jù)的成像域最小二乘逆時(shí)偏移結(jié)果，驗(yàn)證了把空變PSF近似為局域空不變PSF的可行性與正確性；而在后續(xù)的兩個(gè)實(shí)際數(shù)據(jù)的成像域最小二乘逆時(shí)偏移中, 對(duì)成像域進(jìn)行了分塊，使用各個(gè)局部區(qū)塊中心點(diǎn)的PSF(即多點(diǎn)的PSF)，驗(yàn)證了本文提出的局域空不變反褶積的最小二乘逆時(shí)偏移的有效性.局域空不變反褶積對(duì)于不同位置的局域有不同的PSF，計(jì)算效率不如全域空不變反褶積，但它可以有效地考慮PSF在成像域的空變特性，有助于提高成像域最小二乘逆時(shí)偏移的成像效果,這是一種既兼顧PSF的(緩慢)空變特性，又具有計(jì)算高效性的方法.

4 結(jié)論

給定具有地震波運(yùn)動(dòng)學(xué)準(zhǔn)確的光滑偏移速度模型，地震數(shù)據(jù)的一次地震波偏移成像可視為一個(gè)有關(guān)地下反射率分布的線性反演問題.利用最小二乘反演方法，在成像域常規(guī)逆時(shí)偏移結(jié)果是最小二乘反演的Hessian矩陣與最小二乘逆時(shí)偏移結(jié)果的空變積分.Hessian矩陣中的列向量也稱為成像點(diǎn)的點(diǎn)擴(kuò)散函數(shù)，是成像點(diǎn)單位脈沖體所產(chǎn)生的地震響應(yīng)的逆時(shí)偏移結(jié)果，因此點(diǎn)擴(kuò)散函數(shù)具有空間局部性.由于偏移速度模型的光滑性和地震數(shù)據(jù)采集多次覆蓋，使得一定空間范圍內(nèi)成像點(diǎn)的點(diǎn)擴(kuò)散函數(shù)具有很大的相似性.利用三維常規(guī)逆時(shí)偏移結(jié)果和以所有成像點(diǎn)的點(diǎn)擴(kuò)散函數(shù)為列向量的Hessian矩陣求取最小二乘逆時(shí)偏移結(jié)果，不僅涉及形成Hessian矩陣的巨大計(jì)算量，而且還是一個(gè)超大規(guī)模的線性反演問題.利用一定空間范圍內(nèi)成像點(diǎn)的點(diǎn)擴(kuò)散函數(shù)所具有的相似性，可把該空間范圍內(nèi)的常規(guī)逆時(shí)偏移結(jié)果近似為Hessian矩陣與最小二乘逆時(shí)偏移結(jié)果的空不變褶積，再利用求解空不變反褶積的波數(shù)域偽廣義逆方法可穩(wěn)定高效地得到最小二乘逆時(shí)偏移結(jié)果.利用點(diǎn)擴(kuò)散函數(shù)的空間局部性，可同時(shí)計(jì)算成像域內(nèi)以一定空間間隔均勻分布的多個(gè)成像點(diǎn)的點(diǎn)擴(kuò)散函數(shù)，這也是成像域高效最小二乘逆時(shí)偏移方法的組成部分.在本文所提出的最小二乘逆時(shí)偏移方法中，局域區(qū)塊Lb尺寸的大小是計(jì)算效率與偏移成像效果之間的折中.由于深部成像點(diǎn)的點(diǎn)擴(kuò)散函數(shù)的分布范圍大于淺部的點(diǎn)擴(kuò)散函數(shù)的分布范圍，因此局域區(qū)塊Lb尺寸大小的選擇應(yīng)以深部點(diǎn)擴(kuò)散函數(shù)的完整性為原則.Marmousi模型數(shù)據(jù)和實(shí)際數(shù)據(jù)的數(shù)值試驗(yàn)結(jié)果驗(yàn)證了本文所提出的成像域最小二乘逆時(shí)偏移方法的有效性.

致謝本文的研究工作得到了中石化科技部的資助，也得到了中石化石油物探技術(shù)研究院段心標(biāo)、白英哲、陶永慧、孫敏傲等專家的大力幫助，在此表示感謝.