黃建平， 李闖*， 李慶洋， 郭書娟， 段心彪， 李繼光， 趙勝天， 步長城

1 中國石油大學(xué)(華東)地球物理系， 青島 266580 2 中國石化南京物探研究院， 南京 210000 3 勝利油田物探研究院， 山東東營 257000

?

一種基于平面波靜態(tài)編碼的最小二乘逆時(shí)偏移方法

黃建平1， 李闖1*， 李慶洋1， 郭書娟2， 段心彪2， 李繼光3， 趙勝天3， 步長城3

1 中國石油大學(xué)(華東)地球物理系， 青島 266580 2 中國石化南京物探研究院， 南京 210000 3 勝利油田物探研究院， 山東東營 257000

平面波偏移是一種面炮偏移方法，相對(duì)于常規(guī)逐炮偏移，其具有較高的計(jì)算效率.然而常規(guī)平面波偏移方法成像精度低，且成像時(shí)會(huì)產(chǎn)生串?dāng)_噪音.為此，本文在實(shí)現(xiàn)常規(guī)平面波偏移算法基礎(chǔ)上，引入反演思想實(shí)現(xiàn)了基于靜態(tài)平面波編碼的最小二乘偏移理論方法及處理流程，在優(yōu)化算法基礎(chǔ)上對(duì)平層模型和復(fù)雜砂礫斷塊模型進(jìn)行了成像測試并與其他成像策略進(jìn)行對(duì)比.研究結(jié)果表明：基于時(shí)移編碼的平面波最小二乘偏移能有效抑制低頻成像噪音和串?dāng)_噪音，補(bǔ)償中深部成像能量，是一種較為有效的保幅成像策略.

平面波偏移； 巖性成像； 平面波編碼； 最小二乘偏移； 保幅成像

1 引言

近年來，油田勘探開發(fā)已進(jìn)入中、晚期階段，為了滿足目前油田開發(fā)的需要，小面元、高覆蓋次數(shù)的高密度地震采集方法開始得到關(guān)注，特別是在國外海上勘探區(qū)塊，已經(jīng)普遍推廣應(yīng)用高密度空間采樣技術(shù)(夏穎等，2008；趙會(huì)欣等，2008).高密度采集得到的海量數(shù)據(jù)符合平面波偏移對(duì)覆蓋次數(shù)的要求，同時(shí)其計(jì)算效率問題也能通過面炮合成壓縮炮記錄解決( Berkhout, 1992).另一方面，目前地震勘探成像研究逐漸步入精細(xì)成像階段，常規(guī)疊前深度偏移已不能滿足構(gòu)造成像向巖性成像過渡的需求，而基于反演思想的最小二乘偏移是高精度保幅成像的有效方法，但存在計(jì)算效率低的缺陷.因此，基于平面波編碼加速的最小二乘偏移方法具有一定的研究意義.

“面炮”偏移首先由Berkhout(1992)提出, 通過某種合成算子將所有的炮道集合成為一個(gè)面炮記錄,然后采用適于常規(guī)炮集的偏移算法進(jìn)行偏移.Monsher等(1997)提出應(yīng)用Randon變換的偏移距道集數(shù)據(jù)平面波偏移.隨后，張叔倫和孫沛勇(1999)將面炮技術(shù)應(yīng)用于傅里葉有限差分疊前深度偏移,提高了波動(dòng)方程疊前深度偏移的計(jì)算效率.陳生昌等(2002)、馮偉等(2004)通過平面波合成算子合成不同角度的平面波記錄，一系列方向的平面波偏移結(jié)果相疊加,進(jìn)而得到較高精度的地下復(fù)雜介質(zhì)的圖像.為了進(jìn)一步提高計(jì)算效率，孫沛勇和張叔倫(2000)提出了最大能量震源波場的平面波疊前深度偏移方法.此外，葉月明等(2007)、楊海生等(2011)將平面波偏移用于起伏地表的情況.

近年來，最小二乘算法在地震資料處理解釋方面發(fā)揮了越來越重要的作用，其思想由LeBras和Beudoun(1988)等人提出，Lambare等(1992)在他們研究基礎(chǔ)上，使用最速下降法(王彥飛，2007)反演相對(duì)于背景速度的速度擾動(dòng).Nemeth等(1999)早期將最小二乘算法引入到地震偏移以求取地下反射系數(shù)，他通過Kirchhoff最小二乘偏移方法對(duì)不完整的反射地震數(shù)據(jù)進(jìn)行成像.此外，黃建平等(2013a)也應(yīng)用Kirchhoff最小二乘偏移對(duì)我國西部探區(qū)碳酸鹽巖裂縫型儲(chǔ)層進(jìn)行成像研究，劉玉金等(2013)應(yīng)用Kirchhoff最小二乘偏移解決了地震數(shù)據(jù)中采樣不規(guī)則、地震道缺失對(duì)成像造成的影響.隨后，Kuehl和Sachhi (2001, 2003)和黃建平等(2014b)將Stoffa等(1990)提出的裂步算子應(yīng)用到最小二乘偏移，賈曉峰和胡天躍(2005)也提出采用滑動(dòng)最小二乘窗求解波動(dòng)方程.黃建平等(2014c)將分頻編碼應(yīng)用于最小二乘裂步法偏移，在提高計(jì)算效率的同時(shí)壓制了串?dāng)_噪音.為了增加最小二乘偏移方法對(duì)劇烈橫向變速的適應(yīng)性，楊其強(qiáng)和張叔倫(2008)、黃建平等(2013b)實(shí)現(xiàn)了基于最小二乘的疊前深度偏移算法.近年來，國內(nèi)外研究學(xué)者也將逆時(shí)偏移算子應(yīng)用到最小二乘偏移中來進(jìn)行復(fù)雜構(gòu)造的保幅成像處理(Dai and Schuster, 2009, 2013; Dai et al., 2010, 2011, 2012; Li et al, 2014, 黃建平等, 2014a；郭振波和李振春, 2014).

本文在前人工作的基礎(chǔ)上實(shí)現(xiàn)了基于靜態(tài)平面波編碼的最小二乘逆時(shí)偏移方法.為了加快收斂速度，本文加入了基于照明補(bǔ)償預(yù)條件算子，并混合使用共軛梯度法與最速下降法進(jìn)行尋優(yōu).通過成像測試深入分析了該方法在保幅性、成像分辨率以及橫向能量均衡性等方面的優(yōu)勢，解決了常規(guī)成像時(shí)的復(fù)雜小斷塊邊界難以分辨的問題，有利于后期斷塊油氣藏的儲(chǔ)層劃分及油藏描述工作.

2 平面波最小二乘偏移方法原理

2.1 觀測記錄的平面波編碼

對(duì)于一個(gè)二維工區(qū)，對(duì)炮記錄編碼的過程可表示為

(1)

其中，d(xg,t;xs)為觀測炮記錄，δ(t-pxs)為編碼函數(shù)，編碼過程等價(jià)于通過tau-p變換合成平面波記錄，圖1為平面波編碼的原理圖(Zhangetal., 2005)，時(shí)移量pxs隨震源點(diǎn)位置xs線性變化，p為射線參數(shù)，

(2)

其中,θ是地表入射角，v是地表速度.

圖1 平面波編碼示意圖(Zhang et al.，2005)Fig.1 Diagram of plane-wave encoding

2.2 震源的平面波編碼

逆時(shí)偏移時(shí)需要進(jìn)行震源波場的正向延拓，同樣地，對(duì)第i個(gè)震源的編碼過程可表示為

(3)

其中，s(t)為震源函數(shù)，wi(p,t)為時(shí)移后的震源，該過程相當(dāng)于對(duì)各炮點(diǎn)進(jìn)行一次與其位置有關(guān)的時(shí)移，注意對(duì)炮記錄和震源的編碼是一一對(duì)應(yīng)的.編碼完成后，激發(fā)平面波震源，模擬得到正向延拓的震源波場.

平面波記錄的逆時(shí)偏移算子和炮域逆時(shí)偏移(Baysaletal., 1983)相同，由反傳的合成平面波記錄和平面波震源激發(fā)得到的震源波場相關(guān)成像.

2.3 平面波域反偏移算子

反偏移時(shí)仍然要對(duì)震源進(jìn)行平面波編碼，假設(shè)W(p,ω)為按式(3)編碼后的平面波震源，給定一個(gè)背景慢度場s0，亥姆霍茲方程的解可表示為

(4)

其中，G0(x;xp)為背景慢度s0對(duì)應(yīng)的格林函數(shù)，U0(p,x)為該平面波震源對(duì)應(yīng)的波場，xp對(duì)應(yīng)不同參數(shù)p的平面波震源.

假設(shè)慢度擾動(dòng)為δs(x)，真實(shí)速度場可表示為s(x)=s0+δs(x)，求解Helmholtz方程可求得全波場：

(5)

震源項(xiàng)為F(p,x,ω)=-δ(x-xp)W(p,ω)，將s(x)=s0+δs(x)代入(5)式，并舍去高階項(xiàng)O(δs2)得

(6)

將方程(6)的第三項(xiàng)移到右邊，兩邊同乘以格林函數(shù)G0(x;x′)，再在整個(gè)介質(zhì)中積分可得：

=-∫Uδ(x-x′)dx′

=U(x)

右邊 = ∫G0(x;x′)Fdx′-2ω2∫s0δs(x′)U(x′;xp)

×G0(x;x′)dx′

=U0(x)+ω2∫m(x′)U(x′;xp)G0(x;x′)dx′,

(7)

其中，m(x′)=-2s0δs(x′)代表反射率.而全波場U(x′;xp)=W(p,ω)G(x′;xp)，假設(shè)慢度擾動(dòng)足夠小，即G(x′;xp)≈G0(x′;xp)，則Born近似下的散射波場為

U1(p,x) =U(p,x)-U0(p,x)

≈ω2∫m(x′)W(p,ω)G0(x′;xp)G0(x;x′)dx′.

(8)

為了簡化公式，用矢量矩陣符號(hào)來表示Born正傳算子：

(9)

其中，m為偏移剖面或反射系數(shù)的矩陣形式；d是反偏移得到的平面波記錄的矩陣形式；L為Born近似下的平面波域反偏移算子矩陣形式.

2.4 預(yù)條件算子

為了更好地補(bǔ)償深部照明的不足，本文使用一種照明算子進(jìn)行補(bǔ)償(Beydoun and Mendes, 1989; Luo and Schuster, 1991)，該算子可表達(dá)為

(10)

其中，U(p;x,z;t)為射線參數(shù)為p時(shí)、t時(shí)刻的震源波場，I(p,x,z)為射線參數(shù)為p的平面波記錄的照明算子.

2.5 平面波最小二乘逆時(shí)偏移

假設(shè)有Np個(gè)平面波道集，平面波域的誤差函數(shù)可表示為

(11)

其中，di代表第i個(gè)平面波道集，而Li是與該平面波道集有關(guān)的正演算子，mi是與第i個(gè)平面波道集有關(guān)的偏移剖面.本文使用共軛梯度法求解平面波域的誤差函數(shù)：

(12)

其中，I即為2.4節(jié)討論的照明補(bǔ)償預(yù)條件算子.

當(dāng)目標(biāo)函數(shù)為n維的二次可微函數(shù)時(shí),利用共軛梯度法理論上最多只要n次迭代即可達(dá)到極小值點(diǎn)，但在實(shí)際計(jì)算中存在舍入誤差、計(jì)算誤差等因素，n次迭代后往往不能收斂，而n維函數(shù)問題的共軛方向最多只有n個(gè)，所以n次以后的迭代將失去意義，同時(shí)誤差的積累會(huì)對(duì)收斂不利(陳開周，1985).因而，本文使用的共軛梯度法會(huì)在n次迭代之后將方向重設(shè)為最速下降方向重新啟動(dòng)算法.本文所使用的共軛梯度法的流程為：

(1)令迭代步數(shù)k=0，k1=0,設(shè)定允許誤差ε和一輪搜索的最大次數(shù)n，初始搜索方向?yàn)樽钏傧陆捣较颍?/p>

(2)計(jì)算最速下降方向g(k)，若‖g(k)‖≤ε，則停止迭代，輸出當(dāng)前模型m(k)；

(3)判斷是否滿足條件k1

(5)更新反射率模型m(k+1)=m(k)-α(k+1)z(k+1)，k=k+1，k1=k1+1.

平面波最小二乘逆時(shí)偏移的實(shí)現(xiàn)流程如圖2所示，其中初始模型賦0值，即第一次迭代的結(jié)果為逆時(shí)偏移的偏移剖面，在隨后的迭代中通過公式(11)

圖2 平面波最小二乘逆時(shí)偏移實(shí)現(xiàn)流程Fig.2 Flow chart of plane-wave LSRTM

使其向真實(shí)的反射率模型收斂.通過平面波編碼壓縮數(shù)據(jù)將Ns個(gè)炮記錄壓縮為Np個(gè)平面波記錄，最小二乘逆時(shí)偏移的計(jì)算效率可提高Ns/(Np×2)倍.需要強(qiáng)調(diào)的是，在反偏移時(shí)仍然要對(duì)震源進(jìn)行平面波編碼，且與炮記錄的編碼是一一對(duì)應(yīng)的.

圖3 平層模型速度場Fig.3 Horizontal layered velocity model

3 模型試處理

3.1 平層模型

在編程實(shí)現(xiàn)方法的基礎(chǔ)上，本文建立了如圖3所示的平層模型進(jìn)行成像方法正確性測試.其中，模型參數(shù)為：深度方向80個(gè)采樣點(diǎn)，水平方向100個(gè)采樣點(diǎn)，網(wǎng)格間距10 m.其次，計(jì)算參數(shù)如下：震源為主頻30 Hz 的雷克子波，時(shí)間采樣間隔0.5 ms，時(shí)間采樣點(diǎn)2101個(gè)，使用時(shí)間2階、空間8階有限差分正演模擬.觀測系統(tǒng)的設(shè)計(jì)：總炮數(shù)101炮，101道檢波器接收，炮間距和道間距都為10 m.

使用如上參數(shù)計(jì)算得到炮記錄后，根據(jù)公式(3)合成不同參數(shù)p的平面波記錄，p的取值范圍為-0.25～0.25 ms·m-1(對(duì)應(yīng)的傾角范圍-30°～ 30°)，呈線性變化，共合成24個(gè)平面波記錄.圖4展示了3個(gè)不同角度的平面波道集(已切除直達(dá)波).

圖5為對(duì)如圖4b所示的平面波記錄(p=0.0 ms·m-1)RTM成像結(jié)果，剖面中含有較強(qiáng)的低頻噪音、串?dāng)_噪音以及震源效應(yīng)，剖面信噪比及分辨率不高.

該平面波記錄(p=0 ms·m-1)的P-LSRTM成像結(jié)果如圖6所示，從圖6可知，隨著迭代次數(shù)增加震源效應(yīng)減弱；深部反射同相軸能量逐漸增強(qiáng)；成像剖面信噪比和分辨率得到改善；但對(duì)于由稀疏采集產(chǎn)生的串?dāng)_噪音的壓制效果較差.

為了進(jìn)一步壓制串?dāng)_噪音，本文將24個(gè)平面波記錄偏移結(jié)果疊加起來，結(jié)果如圖7所示.由圖7可以看到，串?dāng)_噪音基本被壓制，剖面分辨率和信噪比也得到了改善.對(duì)于不同傾角的平面波記錄，其偏移產(chǎn)生的串?dāng)_噪音相干性差，因此通過疊加不同傾角的平面波偏移結(jié)果，可以使反射同相軸能量增強(qiáng)，壓制串?dāng)_噪音.此外，平面波源具有方向性，其照射方向與目標(biāo)結(jié)構(gòu)方向垂直時(shí)成像分辨率最高(馮偉等，2004).對(duì)于平層模型，p=0 ms·m-1時(shí)成像分辨率最高，來自其他p參數(shù)的成像結(jié)果的疊加可以壓制串?dāng)_噪音，但在一定程度上降低了疊加結(jié)果的分辨x率.當(dāng)模型復(fù)雜時(shí)，這一分辨率差異的現(xiàn)象將不再明顯.

圖4 p取值為-0.25 ms·m-1 (a)， 0 ms·m-1 (b)， 0.25 ms·m-1 (c) 時(shí)的平面波記錄Fig.4 Plane-wave records with p= -0.25 ms·m-1 (a), p=0 ms·m-1 (b), p=0.25 ms·m-1 (c)

圖5 P-RTM偏移結(jié)果Fig.5 Image of P-RTM

為了更清楚地觀察平面波LSRTM對(duì)反射同相軸的能量補(bǔ)償效果，本文提取了偏移距為500 m時(shí)(圖7a實(shí)線所示)P-RTM和P-LSRTM的單道記錄對(duì)比，對(duì)比結(jié)果如圖8a所示，隨迭代次數(shù)增加，P-LSRTM結(jié)果振幅能量增強(qiáng)，更接近雷克子波的形態(tài)，其保幅性優(yōu)于常規(guī)P-RTM.

圖8b展示了P-RTM和P-LSRTM深層反射同相軸(圖7a中虛線位置所示)的振幅對(duì)比，主要探究P-LSRTM對(duì)橫向振幅均衡性的改善效果.由圖8b可知，由于深部的照明效果較差，P-RTM結(jié)果反射振幅弱，且橫向振幅能量不均衡；而P-LSRTM的橫向振幅能量及均衡性均隨迭代次數(shù)增加而變強(qiáng).

P-LSRTM的殘差收斂曲線如圖9所示，由圖9可知，成像反演結(jié)果較為穩(wěn)定；隨著迭代次數(shù)的增加，殘差逐漸減小，成像結(jié)果逐漸趨近于真實(shí)模型；在第10次迭代時(shí)結(jié)果已基本收斂.

通過平層模型正演及偏移成像結(jié)果，驗(yàn)證了本文方法的正確性.

3.2 砂礫斷塊模型

在驗(yàn)證了偏移方法正確性的基礎(chǔ)上，本文引入了勝利某探區(qū)典型的砂礫斷層模型，對(duì)方法的適應(yīng)性進(jìn)行測試.計(jì)算用到的模型如圖10所示，模型水平方向大小為5.7 km，深度方向?yàn)?.0 km，網(wǎng)格間距10 m，上覆巖層厚度約700 m，模型右側(cè)為一個(gè)高速的高陡傾構(gòu)造.模型自上而下依次分布5個(gè)高速小斷塊體，且埋藏較深，傾角較大，斷塊間間距小，分界不明顯.由于高速斷塊體造成的深部照明不足、強(qiáng)烈的橫向變速以及模糊的斷塊邊界都將給斷塊體的高精度成像帶來巨大的挑戰(zhàn).

成像測試的計(jì)算參數(shù)為：震源為主頻30 Hz的雷克子波，網(wǎng)格間距10 m，時(shí)間采樣間隔0.5 ms，采集時(shí)長2 s，使用時(shí)間2階、空間8階有限差分正演模擬.觀測系統(tǒng)的設(shè)計(jì)：總炮數(shù)571炮，571道檢波器接收，炮間距和道間距均為10 m.

使用如上參數(shù)計(jì)算得到炮記錄后，根據(jù)公式(3)合成不同參數(shù)p的平面波記錄，p的取值范圍為-0.1742～0.1742 ms·m-1，呈線性變化，共合成24個(gè)平面波記錄.圖11展示了3個(gè)不同角度的平面波道集(已切除直達(dá)波)，可以看到反射同相軸十分復(fù)雜、錯(cuò)亂.

圖12為24個(gè)平面波記錄P-LSRTM成像后的疊加結(jié)果，由圖12可知，在殘差收斂過程中，由成像產(chǎn)生的低頻噪音和串?dāng)_噪音得到壓制；反射同相軸能量增強(qiáng)；深部反射同相軸能量均衡性變好；但對(duì)成像結(jié)果右側(cè)平層上的強(qiáng)低頻噪音不能完全壓制，這可能是模型右側(cè)的平層反射系數(shù)過大造成的.根據(jù)Zhang等(2005)給出的計(jì)算公式，需要90個(gè)平面波記錄偏移結(jié)果的疊加才能完全壓制串?dāng)_噪音，但是LSRTM本身具有壓制噪音的作用，因此本文只疊加了24個(gè)平面波記錄，串?dāng)_噪音已基本壓制.

圖6 P-LSRTM迭代5次(a)、10次(b)、15次(c)、20次(d)、25次(e)、30次(f)偏移結(jié)果(p=0 ms·m-1)Fig.6 Image of P-LSRTM after 5(a),10(b),15(c),20(d),25(e),30(f) iterations (p=0 ms·m-1)

圖7 P-LSRTM迭代5次(a)、10次(b)、15次(c)、20次(d)、25次(e)、30次(f)的疊加結(jié)果Fig.7 The stacked image of P-LSRTM after 5 (a),10 (b),15 (c),20 (d),25 (e),30 (f) iterations

圖8 保幅性及能量均衡性對(duì)比(a) 單道振幅對(duì)比(Offset=500 m)； (b) 橫向振幅對(duì)比(圖2紅線所示).Fig.8 Comparison of amplitude preservation and energy equilibrium(a) Comparison of vertical amplitude at CDP=50； (b) Comparison of horizontal amplitude (the red line showed in Fig.2).

圖9 殘差收斂曲線Fig.9 Convergence curve of residual

圖10 砂礫斷層模型速度場Fig.10 Gravel fault model velocity model

P-LSRTM的殘差收斂曲線如圖13所示，由圖13可知，成像反演結(jié)果穩(wěn)定；對(duì)于復(fù)雜模型，10次迭代后數(shù)據(jù)殘差基本收斂.

圖14為幾種常用疊前深度偏移算法的成像對(duì)比，由圖14(a—d)可知，RTM成像剖面中含有較強(qiáng)的低頻偏移噪音，在一定程度上掩蓋了地下的真實(shí)構(gòu)造，尤其是深層構(gòu)造難以辨認(rèn)；SSF對(duì)高陡傾構(gòu)造成像效果差；LSRTM成像質(zhì)量最高，低頻噪音得到壓制，深層反射軸能量得到補(bǔ)償，高陡構(gòu)造和小斷塊都準(zhǔn)確地成像；平面波LSRTM成像質(zhì)量接近LSRTM.

圖11 p取值為 -0.1742 ms·m-1 (a)、0 ms·m-1 (b)、0.1742 ms·m-1 (c) 時(shí)的平面波記錄Fig.11 Plane-wave records with p=-0.1742 ms·m-1 (a), p=0 ms·m-1 (b), p=0.1742 ms·m-1 (c)

圖12 P-LSRTM迭代5次(a)、10次(b)、15次(c)、20次(d)、25次(e)、30次(f)的疊加結(jié)果Fig.12 The stacked image of P-LSRTM after 5(a),10(b),15(c),20(d),25(e),30(f) iterations

圖14 各偏移方法成像效果對(duì)比(a) 逆時(shí)偏移； (b) 分步傅里葉法偏移； (c) 最小二乘逆時(shí)偏移(迭代30次)； (d) 平面波最小二乘逆時(shí)偏移(迭代30次).Fig.14 Comparison of images of different pre-stack migration methods(a) RTM； (b) SSF； (c) LSRTM(after 30 iterations)； (d) P-LSRTM(after 30 iterations).

LSRTM迭代一次所用的時(shí)間是RTM的2倍，計(jì)算效率低.通過平面波編碼壓縮數(shù)據(jù)，P-LSRTM的計(jì)算時(shí)間是炮域RTM的24×30×2/571≈2.52倍.RTM偏移時(shí)需要先讀取炮記錄并存入內(nèi)存，P-LSRTM存入內(nèi)存的數(shù)據(jù)量為RTM的24/571≈0.042倍.

4 結(jié)論與討論

本文發(fā)展了基于靜態(tài)平面波編碼的最小二乘逆時(shí)偏移方法，在實(shí)現(xiàn)算法的基礎(chǔ)上，通過對(duì)簡單平層模型和復(fù)雜砂礫斷塊模型的成像試處理以及與其他幾種疊前偏移方法的對(duì)比，得到了如下認(rèn)識(shí)：

(1)基于平面波編碼的LSRTM成像算法能夠較好地抑制RTM算法中的低頻噪音，減弱稀疏采集造成的震源效應(yīng)，補(bǔ)償深層反射同相軸能量，對(duì)復(fù)雜小構(gòu)造的高精度成像具有一定的優(yōu)勢；

(2)串?dāng)_噪音可以通過疊加不同傾角的偏移剖面壓制，此外迭代次數(shù)的增加也對(duì)串?dāng)_噪音有一定的壓制作用；

(3)通過平面波編碼壓縮數(shù)據(jù)體，提高了計(jì)算效率，減小了I/O消耗，提高了最小二乘高精度成像理論方法在生產(chǎn)中應(yīng)用的可能性.

基于靜態(tài)平面波編碼的最小二乘逆時(shí)偏移對(duì)不同角度的平面波記錄分別進(jìn)行偏移，計(jì)算效率高，有利于平面波域共成像點(diǎn)道集和角道集的提取，可用于成像質(zhì)量的控制和偏移速度分析.但是，本文所實(shí)現(xiàn)的平面波編碼只適用于固定的觀測系統(tǒng)，無法處理炮點(diǎn)和接收點(diǎn)分布不固定的情況.因此，角道集的提取以及適用于非固定觀測系統(tǒng)的P-LSRTM方法都將是作者下一步研究的重點(diǎn).

Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration.Geophysics, 48(11): 1514-1524.

Berkhout A J. 1992. Areal shot-record technology.JournalofSeismicExploration, 1(3): 251-264.

Beydoun W B, Mendes M. 1989. Elastic ray-born L2-migration/inversion.GeophysicalJournalInternational, 97(1): 151-160.

Dai W, Schuster J. 2009. Least-squares migration of simultaneous sources data with a deblurring filter. ∥ SEG Houston International Exposition and Annual Meeting. 2990-2994. Dai W, Boonyasiriwat C, Schuster G T. 2010. 3D multi-source least-squares reverse time migration. ∥ 2010 SEG Technical Program Expanded Abstracts 2010: 3120-3124.

Dai W, Wang X, Schuster G T. 2011. Least-squares migration of multisource data with a deblurring filter.Geophysics, 76(5): R135-R146.

Dai W, Fowler P, Schuster G T. 2012. Multi-source least-squares reverse time migration.GeophysicalProspecting, 60(4): 681-695.

Dai W, Schuster G T. 2013. Plane-wave least-squares reverse-time migration.Geophysics, 78(4): S165-S177. Kuehl H, Sachhi M D. 2001. Split-step WKBJ least-square migration / inversion of incomplete data. ∥5th SEGJ International Symposium Imaging Technology. Kuehl H, Sacchi M D. 2003. Least-squares wave-equation migration for AVP/AVA inversion.Geophysics, 68(1): 262-273. Lambare G, Virieux J, Madariaga R, et al. 1992. Iterative asymptotic inversion in the acoustic approximation.Geophysics, 57(9): 1138-1154. LeBras R, Clayton R W. 1988. An iterative inversion of back-scattered acoustic waves.Geophysics, 53(4): 501-508.

Li C, Huang J P, Li Z C, et al. 2014. Application of plane-wave least square migration in fault block reservoirs-a case study. ∥76th EAGE Conference and Exhibition.

Luo Y, Schuster G T. 1991. Wave-equation traveltime inversion.Geophysics, 56(5): 645-653.

Mosher C C, Foster D J, Hassanzadeh S. 1997. Common angle imaging with offset plane waves. ∥Expanded Abstract of 67th Annual Internat SEG Mtg., 1379-1382.

Nemeth T, Wu C J, Schuster G T. 1999. Least-squares migration of incomplete reflection data.Geophysics, 64(1): 208-221.

Stoffa P L, Fokkema J T, de Luna Freire R M, et al. 1990. Split-step Fourier migration.Geophysics, 55(4): 410-421.

Zhang Y, Sun J, Notfors C, et al. 2005. Delayed-shot 3D depth migration.Geophysics, 70(5): E21-E28.

附中文參考文獻(xiàn)

陳開周. 1985. 最優(yōu)化計(jì)算方法. 西安: 西北電訊工程學(xué)院出版社.

陳生昌, 曹景忠, 馬在田. 2002. 平面波偏移. 勘探地球物理進(jìn)展, 25(3): 37-41.

馮偉, 王華忠, 吳如山等. 2004. 面向目標(biāo)控制照明的合成波源偏移. 石油物探, 43(3): 223-228.

郭振波, 李振春. 2014. 最小平方逆時(shí)偏移真振幅成像. 石油地球物理勘探, 49(1): 113-120.

黃建平, 李振春, 孔雪等. 2013a. 碳酸鹽巖裂縫型儲(chǔ)層最小二乘偏移成像方法研究. 地球物理學(xué)報(bào), 56(5): 1716-1725.

黃建平, 李振春, 劉玉金等. 2013b. 復(fù)雜介質(zhì)最小二乘疊前深度偏移成像方法. 地球物理學(xué)進(jìn)展, 28(6): 2977-2983.

黃建平, 曹曉莉, 李振春等. 2014a. 最小二乘逆時(shí)偏移在近地表高精度成像中的應(yīng)用. 石油地球物理勘探, 49(1): 107-112.

黃建平, 薛志廣, 步長城等. 2014b. 基于裂步DSR的最小二乘偏移方法. 吉林大學(xué)學(xué)報(bào): 地球科學(xué)版, 44(1): 369-374.

黃建平, 孫鄖松, 李振春等. 2014c. 一種基于分頻編碼的最小二乘裂步偏移方法. 石油地球物理勘探, 49(4): 702-707.

賈曉峰, 胡天躍. 2005. 滑動(dòng)最小二乘法求解地震波波動(dòng)方程. 地球物理學(xué)進(jìn)展, 20(4): 920-924.

劉玉金, 李振春, 吳丹等. 2013. 局部傾角約束最小二乘偏移方法研究. 地球物理學(xué)報(bào), 56(3): 1003-1011.

孫沛勇, 張叔倫. 2000. 平面波最大能量疊前深度偏移. 石油地球物理勘探, 35(3): 283-289.

王彥飛. 2007. 反演問題的計(jì)算方法及其應(yīng)用. 北京: 高等教育出版社.

夏穎, 祝彩霞, 孫靈群. 2008. 地震勘探儀器在高密度采集中的應(yīng)用. 物探裝備, 18(1): 7-10, 21.

楊海生, 王必金, 陳玉明, 等. 2011. 基于起伏地表的合成平面波疊前深度偏移. 石油物探, 50(2): 129-138.

楊其強(qiáng), 張叔倫. 2008. 最小二乘傅立葉有限差分法偏移. 地球物理學(xué)進(jìn)展, 23(2): 433-437.

葉月明, 李振春, 仝兆岐等. 2007. 起伏地表?xiàng)l件下的合成平面波偏移及其并行實(shí)現(xiàn). 石油地球物理勘探, 42(6): 622-628.

張叔倫, 孫沛勇. 1999. 基于平面波合成的傅里葉有限差分疊前深度偏移. 石油地球物理勘探, 34(1): 1-7.

趙會(huì)欣, 晉志剛, 張宇生等. 2008. 高密度空間采樣地震采集覆蓋次數(shù)的選擇. 天然氣工業(yè), 27(增刊A): 68-69.

(本文編輯 胡素芳)

Least-squares reverse time migration with static plane-wave encoding

HUANG Jian-Ping1, LI Chuang1*, LI Qing-Yang1, GUO Shu-Juan2, DUAN Xin-Biao2,LI Ji-Guang3, ZHAO Sheng-Tian3, BU Chang-Cheng3

1DepartmentofGeophysics,ChinaUniversityofPetroleum(EastChina),Qingdao266580,China2SINOPECGeophysicalResearchInstitute,Nanjing210000,China3ShengliGeophysicalResearchInstituteofSINOPEC,DongyingShandong257000,China

To solve the exploration problem of widely distributed fault block reservoirs in China, high-density seismic acquisition methods with small surface element and high coverage times become more popular recently which results in the huge dataset for seismic processing. Plane-wave migration is a kind of areal shot migration method which can process enormous data of high-density acquisition with much higher computational efficiency compared with shot domain migration. But conventional plane-wave migration produces low quality images with crosstalk. To solve these problems, the paper introduces the plane-wave migration to the inversion framework and presents the theory and work flow of least-squares reverse time migration with static plane-wave encoding.The key points of plane-wave least squares reverse time migration can be divided into two parts. Firstly, the plane-wave decomposition is applied to the shot data which can be performed by using some delayed-shot variant of slant-stack processing. This involves applying a linear time delay to the shot records. And the time delay is also applied to the corresponding sources. Note that there is one-to-one correspondence between source encoding and encoding to common shot gathers. The input data of migration is compressed a lot after plane-wave encoding which improves the computational efficiency. Secondly, the misfit function of plane-wave least squares reverse time migration is defined as the data difference between predicted plane-waves by Born modelling and observed plane-waves in the pre-stack migration. The conjugate gradient method is implemented to find the solution that minimizes the misfit function. However, the searching of CG method becomes very slow after several iterations in practical application, so we set the steepest descent direction as the gradient after several iterations to restart the CG method. One other thing to note is that the migration results are updated separately with different plane-waves and stacked together after the iterations are finished.The plane-wave least squares reverse time migration is firstly tested via the synthetic data set with the three layer medium model to verify the validity of the proposed method. The plane-wave reverse time migration is also applied to the model for better comparison. The comparison of imaging quality, vertical amplitude and transverse amplitude balance are presented which illustrates the advantages of the proposed method. The stable convergence also proves the robustness of the method. Then, the proposed method is applied to a complex fault block model in Shengli Oil-field to test its applicability. Several fault blocks are distributed in the mid-deep part and a high steep structure locates on the right of the model. The strong velocity variation, the burial depth and fuzzy boundary of fault blocks make it a good velocity model to test the resolution of the imaging method. The imaging results of the proposed method is compared with the reverse time migration, the split step migration and least-squares revere time migration results, which show that the proposed method can produce high quality images similar to least-squares reverse time migration but onlynp/nstimes computation is needed (nsis the number of the shots andnpis the number of the encoded plane-waves).The paper presents the theory of prestack least-squares reverse time migration with a static plane-wave encoding technique. From the imaging test we can get the conclusions including: (1) the proposed method has the advantages of resolution enhancement and amplitude compensation in mid-deep part with the increase of iterations; (2) the crosstalk introduced by plane-waves can be suppressed by the stacking of images from different plane-waves and the optimization; (3) the computation and input/output cost is reduced a lot with the plane-waves encoding. Besides, the common image gathers can be extracted efficiently from the migration results of the proposed method which is also available for the migration velocity analysis and our next work will be focused on this part. However, the proposed method is only implemented with the fixed geometry and another point of our next work is to modify the method to fit irregular geometry.

Plane-wave migration； Lithology imaging； Plane-wave encoding； Least-square migration； Amplitude persevered migration

10.6038/cjg20150619.

國家973課題(2014CB239006，2011CB202402)，國家自然科學(xué)基金(41104069，41274124)，山東省自然科學(xué)基金(ZR2011DQ016)，中央高?？蒲袠I(yè)務(wù)費(fèi)專項(xiàng)基金(R1401005A)聯(lián)合資助.

黃建平，男，1982年生，副教授，主要從事地震波正演及偏移成像工作.E-mail：jphuang@mail.ustc.edu.cn

*通訊作者 李闖,男, 1992年生, 博士生，主要研究方向?yàn)榈卣鸩ㄗ钚《四鏁r(shí)偏移. E-mail:chli0409@126.com

10.6038/cjg20150619

P631

2014-05-28，2015-03-19收修定稿

黃建平， 李闖， 李慶洋等. 2015. 一種基于平面波靜態(tài)編碼的最小二乘逆時(shí)偏移方法.地球物理學(xué)報(bào)，58(6):2046-2056,

Huang J P, Li C, Li Q Y, et al. Least-squares reverse time migration with static plane-wave encoding.ChineseJ.Geophys. (in Chinese),58(6):2046-2056,doi:10.6038/cjg20150619.