陳生昌， 李代光， 金成玫

浙江大學地球科學學院, 杭州 310027

0 引言

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

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

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

1 最小二乘逆時偏移簡介

給定具有地震波運動學準確的光滑偏移速度模型，作為一種線性反演的地震數據偏移所對應的地震數據線性正演方程可寫為

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

(1)

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

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

(2)

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

(3)

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

(4)

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

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

d=Lm，

(5)

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

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

(6)

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

(7)

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

對于方程(6)可寫為

mls=H-1mc,

(8)

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

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

(9)

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

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

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

mc=Hmls.

(10)

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

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

(11)

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

(12)

式中，x為PSF在成像域中的計算點位置；x0代表成像域中的成像點位置；PSF(x,x0)為成像點x0的點擴散函數；Re表示取實部運算.

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

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

(13)

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

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

(14)

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

(15)

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

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

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

圖2 9個不同網格點位置成像點的PSF分布圖Fig.2 The PSF distribution of 9 imaging points at different grid-point positions

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

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

(16)

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

(17)

(18)

(19)

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

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

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

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

3 數值試驗

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

3.1 Marmousi模型合成地震數據

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

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

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

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

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

3.2 二維實際地震數據

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

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

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

圖10 二維實際數據的成像域最小二乘逆時偏移成像結果Fig.10 The IDLSRTM result of 2D field data

3.3 三維實際地震數據

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

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

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

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

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

圖13 縱測線131的成像域最小二乘逆時偏移成像結果Fig.13 The IDLSRTM result of In-Line 131

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

圖15 縱測線231的成像域最小二乘逆時偏移成像結果Fig.15 The IDLSRTM result of In-Line 231

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

圖17 橫測線100的成像域最小二乘逆時偏移成像結果Fig.17 The IDLSRTM result of Cross-Line 100

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

4 結論

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

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