黃建平，薛志廣，步長城，李振春，王常波，高國超，曹曉莉，李國磊

1.中國石油大學（華東）地球科學與技術(shù)學院，山東 青島 266580

2.勝利油田物探研究院，山東 東營 257022

0 引言

偏移是利用反射地震數(shù)據(jù)對地下構(gòu)造進行成像的重要手段。對炮域偏移方法而言，經(jīng)典的成像條件認為：正向傳播的震源波場與反向傳播的接收波場互相關(guān)后，可確定反射層的位置[1]。但是事實上，這種成像準則僅僅是正演算子的共軛轉(zhuǎn)置[2]，而不是它的逆。此外，由于采集孔徑的限制，速度模型復(fù)雜以及波場頻寬有限，常規(guī)的偏移方法通常對地下結(jié)構(gòu)模糊成像，僅能夠提供較準確的構(gòu)造信息，而對中深部界面反射系數(shù)的刻畫不夠準確，這明顯無法滿足巖性油氣藏勘探開發(fā)的要求。

為了解決這些問題，將成像視為最小二乘意義下的反演問題，通過共軛梯度法迭代使誤差函數(shù)達到最小，得到橫向分辨率更高、振幅保真性更好的反演成像結(jié)果。最小二乘偏移早期發(fā)展由LeBras等[3]和 G.Lambare等[4]完成。Nemeth等[5]提出了基于Kirchhoff的最小二乘偏移方法，并證明了最小二乘偏移（LSM）方法可以減少偏移假頻。Duquet等[6]證實了最小二乘偏移在處理起伏地表照明和由不規(guī)則粗采樣的地震波場引起的成像誤差時比Kirchhoff偏移具有更大優(yōu)勢。通過引入數(shù)據(jù)協(xié)方差矩陣，Kuehl等[7]在理論上證明了最小二乘反演的思想可以應(yīng)用在相移偏移上。Clapp[8]驗證了經(jīng)過幾次迭代后，LSM比常規(guī)偏移具有更好的成像精度。

國內(nèi)也有許多學者將最小二乘思想用于偏移、去噪以及數(shù)據(jù)規(guī)則化的研究工作中：賈曉峰等[9]實現(xiàn)波動方程算子的最小二乘偏移；楊其強等[10]進一步發(fā)展了基于傅里葉有限差分的最小二乘偏移方法，并通過模型驗證了成像效果；沈雄君等[11]詳細介紹了基于分步傅里葉變換的最小二乘偏移方法，并給出了Marmousi模型的計算結(jié)果；黃建平等[12]推導并實現(xiàn)了基于Kirchhoff成像算子的疊前最小二乘偏移算法。

在前人研究的基礎(chǔ)上，筆者推導并實現(xiàn)了基于裂步雙均方根（DSR）的最小二乘偏移算法，并將算法用于平層和Marmousi模型進行偏移試算，以測試LSM對復(fù)雜地質(zhì)構(gòu)造成像的能力，并同常規(guī)的裂步DSR偏移算法進行成像效果對比，證實LSM算法的優(yōu)越性。

1 方法原理

1.1 LSM基本原理

已知線性正演算子L，由地下反射率模型m，可得到合成地震數(shù)據(jù)d。正演過程可描述為

未加約束條件的最小二乘偏移問題等價于優(yōu)化特定的目標函數(shù)S（m）：

其中：dobs為觀測到的記錄數(shù)據(jù)；‖＊‖為L2范數(shù)。

能夠使目標函數(shù)S（m）最小化的最小二乘偏移解可表示為

其中：LT為偏移算子（L的共軛轉(zhuǎn)置）；mmig為偏移得到的像，mmig＝LTdobs；H＝LTL為目標函數(shù)S（m）的Hessian矩陣。

由（3）式可知求解最小二乘偏移解～m的2種方法：一是直接求解H，二是間接計算，通過共軛梯度法（或最速下降法）不斷地優(yōu)化模型m。Hessian矩陣方法計算速度快，但占用大量內(nèi)存，且直接求取Hessian矩陣的逆較困難；間接計算方法在每次迭代中均要進行一次正演和偏移計算，計算時間比一般偏移的計算時間多出一個數(shù)量級，為正演過程和偏移過程時間和的n倍，其中n為迭代次數(shù)，但成像算子計算過程中不額外增加存儲空間。筆者采用間接計算方法，用共軛梯度方法來擬合目標函數(shù)，通過迭代的次數(shù)來控制擬合程度，最終得到符合精度要求的成像結(jié)果。

1.2 裂步DSR正傳播和反傳播算子

將地震記錄向下半空間延拓，可求出地下任何一點的波場，進而實現(xiàn)地震波偏移。本文中，假設(shè)向下延拓記為反傳播，向上延拓記為正傳播。其中，正傳播算子為反傳播算子的共軛算子。為此，在推導二者表示形式的過程中可先推導反傳播算子，在得到反傳播算子的基礎(chǔ)上進行共軛轉(zhuǎn)置求取正傳播算子，并用點積標定[13]驗證二者的共軛性。

為了滿足DSR疊前偏移在橫向變速介質(zhì)中精確成像的要求，Popovici[14]對雙均方根方程進行了改進，提出了裂步DSR偏移的算法。從深度z到z＋Δz的延拓過程分為2步：在頻率-波數(shù)域運用相移項P；在頻率-空間域采用時移校正項T。二者可分別表示為

其中：v（z）為參照速度，在層內(nèi)Δz上是一個常數(shù)；ky為中心點波數(shù)；kh為偏移距波數(shù)；ω為頻率。

其中：v（y-h，z）為深度z上炮點位置的速度；v（y＋h，z）為深度z上接收點位置的速度。將此兩項結(jié)合起來，可得到反傳播算子表達式：

其中：F、F-1分別為中心點、偏移距的二維傅里葉正逆變換，且二者互為共軛；U（y，h，ω，z）為正傳播算子。

正傳播算子為反傳播算子的共軛算子，則正傳播算子可表示為

其中：PT、TT分別為P和T的共軛轉(zhuǎn)置，在計算過程中共軛項之間的冪指數(shù)符號相反。

為了確保正傳播算子和反傳播算子互為共軛，需要用點積實驗進行檢驗，算子應(yīng)滿足

其中：輸入向量x，y可取任意值。通過大量模型計算證實：式（6）、（7）所示算子滿足式（8）的共軛關(guān)系。為此，可根據(jù)式（6）、（7）進行裂步DSR的偏移和正演計算過程。最終綜合式（6）、（7）與（2），即可實現(xiàn)基于共軛梯度法的裂步DSR最小二乘偏移方法。

2 模型試算

2.1 簡單平層模型試算

圖1a為該模型速度場，水平采樣間隔10m，深度采樣間隔10m，最大深度為4km。通過反偏移方法得到平層偏移所需要的疊前炮記錄數(shù)據(jù)。

圖1 簡單平層模型試算Fig.1 Simple flat layer model trial

圖1b為常規(guī)裂步DSR偏移的結(jié)果，圖1c為LSM迭代15次的結(jié)果，模型初始值mini＝0，二者的試算結(jié)果采用相同的增益顯示。對比深度3km處的分界面成像效果可以發(fā)現(xiàn)：經(jīng)共軛梯度法15次迭代的最小二乘偏移得到的同向軸更細，成像分辨率較高，能量收斂性較好，較常規(guī)裂步DSR偏移具有更好的保幅性。

圖1d給出了15次迭代過程中對應(yīng)的殘差變化關(guān)系?？梢钥吹剑S著迭代次數(shù)的增加，成像結(jié)果與真實模型的殘差在逐步縮減，這意味著在迭代過程中模型和記錄數(shù)據(jù)匹配的精度越來越高，成像效果越來越好。在迭代次數(shù)較小時，殘差下降迅速，成像結(jié)果收斂較快，成像效果改善明顯；當?shù)螖?shù)較大時，殘差變化趨緩，對成像效果的改善效果減弱。

2.2 Marmousi模型試算

在通過簡單的平層模型驗證了本文方法正確性的基礎(chǔ)上，進一步通過國際標準的SEG／EAGE Marmousi模型數(shù)據(jù)，來檢驗本文最小二乘偏移方法對復(fù)雜模型的適應(yīng)性。Marmousi模型的速度場如圖2a所示，具體參數(shù)為：橫向497個采樣點，縱向750個采樣點，速度場水平采樣間隔12.5m，深度采樣間隔為4m，最大深度3km?；跇藴实?D聲波有限差分法算法模擬得到單炮記錄，并抽取出炮記錄對應(yīng)的共中心點道集如圖2b所示。

為了減少計算成本，偏移算法和LSM算法均采用25個頻率進行偏移成像，頻率為10～60Hz，成像結(jié)果如圖2c、d。對比圖2c和d中圓圈標示的位置可知：裂步DSR偏移算法可以對復(fù)雜地質(zhì)剖面精確成像，而LSM經(jīng)過5次迭代后的效果要優(yōu)于常規(guī)裂步DSR偏移算法?？偟恼f來，LSM成像剖面的地層信息更豐富，分辨率更高。

圖2 Marmousi模型試算Fig.2 Marmousi model trial

圖3 裂步DSR＿LSM和常規(guī)偏移相同位置處成像道振幅對比圖Fig.3 The comparison of image trace amplitude between split-step DSR＿LSM and common migration

從圖2c和d上分別抽取150、250和350網(wǎng)格點位置處的成像道，如圖2d中豎線標示位置，提取常規(guī)法偏移和裂步DSR＿LSM成像剖面的振幅結(jié)果進行對比，成像結(jié)果如圖3所示。對于淺部構(gòu)造，裂步DSR偏移和最小二乘偏移得到的振幅幅值和變化趨勢基本吻合，首先證明了本文方法的正確性。同時，針對中深層，最小二乘偏移具有較好的保幅性，具體體現(xiàn)在LSM單道偏移結(jié)果振幅幅值大于常規(guī)裂步法偏移得到的振幅幅值，如圖3中虛線框所示。通過成像剖面及單道結(jié)果對比可知：LSM可以補償由于斷層和背斜引起的深層能量損失，突顯深層局部結(jié)構(gòu)，對中深部構(gòu)造成像有一定保幅性。裂步法最小二乘偏移之所以能夠?qū)崿F(xiàn)對中深部儲層的較好保幅性，可能是因為常規(guī)偏移方法偏移算子為正演算子的伴隨矩陣而非真正的逆矩陣，隨著傳播路徑的增加，伴隨矩陣與逆矩陣的差異較大，所以常規(guī)偏移中深部不能很好地實現(xiàn)保幅；而裂步法最小二乘偏移在迭代過程中，不斷優(yōu)化伴隨矩陣使之與逆矩陣算子較為接近，從而可從理論上使得中深部儲層具有更好的保幅性[15-17]①黃建平，曹曉莉，李振春，等，最小二乘逆時偏移在近地表高精度成像中的應(yīng)用.石油地球物理勘探，2013，待刊。。

3 結(jié)論與討論

本研究實現(xiàn)了基于雙平方根波場延拓算子的最小二乘偏移方法，通過2個理論模型成像試算，驗證了LSM經(jīng)過幾次迭代后的成像效果優(yōu)于常規(guī)偏移。LSM成像分辨率優(yōu)于常規(guī)偏移方法主要體現(xiàn)在中深層能量得到了補償以及中深層構(gòu)造信息得到了更為準確的刻畫。這一優(yōu)勢，對我國西部儲層深度在3km以下的碳酸鹽巖儲層[18-20]具有非常高的實用性。因此，本方法有望在未來西部勘探開發(fā)中發(fā)揮重要作用。

裂步DSR最小二乘偏移方法相對于Kirchhoff等射線類的最小二乘偏移方法，成像精度更高，不需要計算反射率模型，受地下介質(zhì)密度誤差影響較小。而相對于雙程波的最小二乘逆時偏移（LSRTM）方法，成像效率更高，迭代收斂速度更快，同時對于初始速度模型的依賴性更弱，不易陷入局部極小值。再次，裂步DSR最小二乘偏移其反偏移算子也較雙程波LSRTM易于求取，但其成像精度低于LSRTM方法?？傊?，相對于射線類和雙程波類的最小二乘偏移方法，裂步DSR兼顧計算效率和成像精度。

當然也應(yīng)該看到裂步法最小二乘偏移仍然存在一些不足：1）應(yīng)用LSM對地下構(gòu)造進行成像時，要考慮殘差的變化趨勢。在較小迭代次數(shù)時，殘差變化明顯；而迭代次數(shù)較大時，殘差變化較小，對成像效果的改善能力降低。2）由于每次迭代都要進行一次正演和偏移計算，LSM算法的計算量較大，在應(yīng)用LSM時，要綜合考慮成像精度和計算效率，將二者達到最佳的平衡點。3）在下一步的研究中，也將進一步引入相位編碼技術(shù)[21-22]，以提高最小二乘偏移的成像效率，同時壓制串擾噪音。4）在今后的工作中，逐步將算法用于實際資料的試處理，為復(fù)雜構(gòu)造中深部儲層保幅成像研究服務(wù)。

（References）：

[1]Claerbout J F.Towards a Unified Theory of Reflector Mapping[J].Geophysics，1971，36：467-481.

[2]Lailly P.The Seismic Inverse Problem as a Sequence of Before Stack Migration[C]／／Conference on Inverse Scattering，Theory and Applications.Philadelphia：Soc Industr Appl Math，1983：206-220.

[3]LeBras R，Clayton R W.An Iterative Inversion of Back-Scattered Acoustic Waves[J].Geophysics，1988，53：501-508.

[4]Lambare G，Virieus J，Madariaga R，et al.Iterative Asymptotic Inversion in the Acoustic Approximation[J].Geophysics，1992，57：1138-1154.

[5]Nemeth T，Wu C，Schuster G T.Least-Squares Migration of Incomplete Reflection Data[J].Geophysics，1999，64：208-221.

[6]Duquet B，Marfurt J K，Dellinger J A.Kirchhoff Modeling，Inversion for Reflectivity，and Subsurface Illumination[J].Geophysics，2000，65：1195-1209.

[7]Kuhl H，Sacchi W D.Least-Squares Wave-Equation Migration for AVP／AVA Inversion[J].Geophysics，2003，68（1）：262-273.

[8]Clapp M L.Imaging Under Salt：Illumination Compensation by Regularized Inversion[D].Stanford：Stanford University，2005.

[9]賈曉峰，胡天躍.滑動最小二乘法求解地震波波動方程[J].地球物理學進展，2005，20（4）：920-924.Jia Xiaofeng，Hu Tianyue.Solving Seismic Wave Equation by Moving Least Squares（MLS）Method[J].Progress in Geophysics，2005，20（4）：920-924.

[10]楊其強，張叔倫.最小二乘傅立葉有限差分偏移[J].地球物理學進展，2008，23（2）：433-437.Yang Qiqiang，Zhang Shulun.Least-Squares Fourier Finite-Difference Migration[J].Progress in Geophysics，2008，23（2）：433-437.

[11]沈雄君，劉能超.裂步法最小二乘偏移[J].地球物理學進展，2012，27（2）：761-770.Shen Xiongjun，Liu Nengchao.Split-Step Least-Squares Migration[J].Progress in Geophysics，2012，27（2）：761-770.

[12]黃建平，李振春，劉玉金，等.碳酸鹽巖裂縫型儲層最小二乘偏移成像方法研究[J].地球物理學報，2013，56（5）：1716-1725.Huang Jianping，Li Zhenchun，Liu Yujin，et al.A Study of Least-Squares Migration Imaging Method for Fractured-Type Carbonate Reservoir[J].Chinese Journal of Geophysics，2013，56（5）：1716-1725.

[13]Claerbout J F.Earth Soundings Analysis：Processing Versus Inversion[M].Oxford：Blackwell Scientific Publications，2004.

[14]Popovici A M.Prestack Migration by Split-Step DSR[J].Geophysics，1995，59：1412-1416.

[15]Yao G，Jakubowicz H.Least-Squares Reverse Time Migration[C]／／74th EAGE Conference ＆ Exhibition Incorporating SPE.Copenhagen：[s.n.]，2012.

[16]Dai W，F(xiàn)owler P，Schuster G T.Multi-Source Least-Squares Reverse Time Migration[J].Geophysical Prospecting，2012，60：681-695.

[17]Schuster G T，Wang X，Huang Y，et al.Theory of Multisource Crosstalk Reduction by Phase-Encoded Statics[J].Geophysical Journal International，2011，184：1289-1303.

[18]撒利明，姚逢昌，狄?guī)妥?，?縫洞型儲層地震識別理論與方法[M].北京：石油工業(yè)出版社，2009.Sa Liming，Yao Fengchang，Di Bangrang，et al.Vuggy Reservoir Seismic Recognition Theory and Methods[M].Beijing：Petroleum Industry Press，2009.

[19]姚姚，唐文榜.深層碳酸鹽巖巖溶風化殼洞縫型油氣藏可檢測性的理論研究[J].石油地球物理勘探，2003，38（6）：623-629.Yao Yao，Tang Wenbang.Carbonate Karst Weathering Crust Deep Hole Seam Reservoirs Can Be Detected Theoretical Study[J].Oil Geophysical Prospecting，2003，38（6）：623-629.

[20]吳俊峰，姚姚，撒利明.碳酸鹽巖特殊孔洞型構(gòu)造地震響應(yīng)特征分析[J].石油地球物理勘探，2007，42（2）：180-185.Wu Junfeng，Yao Yao，Sa Liming.The Feature Analysis About Seismic Response to Carbonate Special Hole Type Structure[J].Oil Geophysical Prospecting，2007，42（2）：180-185.

[21]Romero L A，Ghiglia D C.Phase Encoding of Shot Records in Prestack Migarion[J].Geophysics，2000，426-436.

[22]Dai W，Wang X.Least-Squares Migration of Multisource Data with a Deblurring Filter[J].Geophysics，2011，76（5）：135-146.