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

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

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

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

0 引言

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

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

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

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

1 方法原理

1.1 LSM基本原理

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

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

其中：dobs為觀測(cè)到的記錄數(shù)據(jù)；‖＊‖為L(zhǎng)2范數(shù)。

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

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

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

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

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

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

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

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

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

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

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

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

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

2 模型試算

2.1 簡(jiǎn)單平層模型試算

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

圖1 簡(jiǎn)單平層模型試算Fig.1 Simple flat layer model trial

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

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

2.2 Marmousi模型試算

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

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

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

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

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

3 結(jié)論與討論

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

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

當(dāng)然也應(yīng)該看到裂步法最小二乘偏移仍然存在一些不足：1）應(yīng)用LSM對(duì)地下構(gòu)造進(jìn)行成像時(shí)，要考慮殘差的變化趨勢(shì)。在較小迭代次數(shù)時(shí)，殘差變化明顯；而迭代次數(shù)較大時(shí)，殘差變化較小，對(duì)成像效果的改善能力降低。2）由于每次迭代都要進(jìn)行一次正演和偏移計(jì)算，LSM算法的計(jì)算量較大，在應(yīng)用LSM時(shí)，要綜合考慮成像精度和計(jì)算效率，將二者達(dá)到最佳的平衡點(diǎn)。3）在下一步的研究中，也將進(jìn)一步引入相位編碼技術(shù)[21-22]，以提高最小二乘偏移的成像效率，同時(shí)壓制串?dāng)_噪音。4）在今后的工作中，逐步將算法用于實(shí)際資料的試處理，為復(fù)雜構(gòu)造中深部?jī)?chǔ)層保幅成像研究服務(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]賈曉峰，胡天躍.滑動(dòng)最小二乘法求解地震波波動(dòng)方程[J].地球物理學(xué)進(jìn)展，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]楊其強(qiáng)，張叔倫.最小二乘傅立葉有限差分偏移[J].地球物理學(xué)進(jìn)展，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].地球物理學(xué)進(jìn)展，2012，27（2）：761-770.Shen Xiongjun，Liu Nengchao.Split-Step Least-Squares Migration[J].Progress in Geophysics，2012，27（2）：761-770.

[12]黃建平，李振春，劉玉金，等.碳酸鹽巖裂縫型儲(chǔ)層最小二乘偏移成像方法研究[J].地球物理學(xué)報(bào)，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ī)妥?，?縫洞型儲(chǔ)層地震識(shí)別理論與方法[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]姚姚，唐文榜.深層碳酸鹽巖巖溶風(fēng)化殼洞縫型油氣藏可檢測(cè)性的理論研究[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.