基于裂步DSR的最小二乘偏移方法

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

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

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

0 引言

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

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

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

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

1 方法原理

1.1 LSM基本原理

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

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

其中：dobs為觀測到的記錄數據；‖＊‖為L2范數。

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

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

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

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

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

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

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

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

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

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

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

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

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

2 模型試算

2.1 簡單平層模型試算

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

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

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

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

2.2 Marmousi模型試算

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

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

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

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

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

3 結論與討論

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

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

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

（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，Fowler 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]撒利明，姚逢昌，狄幫讓，等.縫洞型儲層地震識別理論與方法[M].北京：石油工業出版社，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]吳俊峰，姚姚，撒利明.碳酸鹽巖特殊孔洞型構造地震響應特征分析[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.