柯 璇,石 穎,2,張瑩瑩,張 振

(1.東北石油大學(xué)地球科學(xué)學(xué)院,黑龍江大慶163318;2.黑龍江省普通高?？萍紕?chuàng)新團(tuán)隊(duì)“斷層變形、封閉性及與流體運(yùn)移”,黑龍江大慶163318;3.中國(guó)石油塔里木油田分公司勘探開(kāi)發(fā)研究院,新疆庫(kù)爾勒841000)

地震疊前逆時(shí)偏移衰減隨機(jī)邊界條件研究

柯 璇1,石 穎1,2,張瑩瑩1,張 振3

(1.東北石油大學(xué)地球科學(xué)學(xué)院,黑龍江大慶163318;2.黑龍江省普通高?？萍紕?chuàng)新團(tuán)隊(duì)“斷層變形、封閉性及與流體運(yùn)移”,黑龍江大慶163318;3.中國(guó)石油塔里木油田分公司勘探開(kāi)發(fā)研究院,新疆庫(kù)爾勒841000)

逆時(shí)偏移的邊界條件及存儲(chǔ)需求往往相互制約,二者的選取直接影響逆時(shí)偏移的成像結(jié)果和適用條件,因此對(duì)逆時(shí)偏移邊界條件和對(duì)應(yīng)存儲(chǔ)策略的研究具有重要意義。提出了一種衰減隨機(jī)邊界條件,在真實(shí)波場(chǎng)外圍依次加入衰減層和隨機(jī)速度層,可有效降低隨機(jī)反射對(duì)真實(shí)波場(chǎng)的能量干擾。對(duì)波動(dòng)方程逆向運(yùn)算重建激發(fā)點(diǎn)波場(chǎng),結(jié)合互相關(guān)成像條件,完成逆時(shí)偏移計(jì)算。對(duì)比分析隨機(jī)邊界條件、完全匹配層邊界條件和衰減隨機(jī)邊界條件的測(cè)試結(jié)果,與隨機(jī)邊界條件相比,衰減隨機(jī)邊界條件在犧牲較少存儲(chǔ)效率的情況下,能夠獲得接近吸收邊界條件的偏移效果。鹽丘模型和大偏移距模型算例分析結(jié)果驗(yàn)證了以上結(jié)論的正確性。

逆時(shí)偏移;隨機(jī)邊界;存儲(chǔ);波場(chǎng)重建

逆時(shí)偏移[1]對(duì)復(fù)雜構(gòu)造成像具有良好的效果,因此成為地球物理研究的熱點(diǎn)之一[2-5]。而計(jì)算精度的提高和波場(chǎng)存儲(chǔ)需求的降低則一直是逆時(shí)偏移算法的重要研究方向。計(jì)算精度方面,除提升波場(chǎng)模擬精度外,還需要采用適合的邊界條件壓制邊界反射,削弱干擾噪聲。BéRENGER[6]提出了完全匹配層(Perfectly Matched Layer,PML)邊界條件,在波場(chǎng)外側(cè)設(shè)置衰減層并引入衰減因子,在衰減層內(nèi)采用帶衰減項(xiàng)的波動(dòng)方程進(jìn)行波場(chǎng)遞推計(jì)算,進(jìn)而達(dá)到消除邊界反射的目的,該方法易于實(shí)現(xiàn)且吸收效果較好[7-12]。張衡等[13]提出了VTI介質(zhì)非分裂式PML吸收邊界條件,人工反射吸收效率得到明顯提高。LIU等[14]提出一種混合邊界條件,衰減層內(nèi)的波場(chǎng)計(jì)算由雙程波方程逐漸轉(zhuǎn)化為單程波方程,再對(duì)波場(chǎng)進(jìn)行衰減吸收,該方法所需吸收層數(shù)較少,間接降低了波場(chǎng)的計(jì)算量和存儲(chǔ)量。在逆時(shí)偏移算法中,由于激發(fā)點(diǎn)和檢波點(diǎn)波場(chǎng)往往需進(jìn)行同時(shí)刻互相關(guān)計(jì)算,但全波場(chǎng)保存的存儲(chǔ)量占用在實(shí)際應(yīng)用中難以實(shí)現(xiàn),因此,逆時(shí)偏移的存儲(chǔ)需求成為限制逆時(shí)偏移發(fā)展的瓶頸之一。對(duì)此,SYMES[15]設(shè)置了檢查點(diǎn)存儲(chǔ)臨時(shí)波場(chǎng),節(jié)省了存儲(chǔ)空間,但需引入額外計(jì)算量用于波場(chǎng)重建;FENG等[16]提出一種波場(chǎng)存儲(chǔ)方法,存儲(chǔ)邊界處的單層波場(chǎng)及最后兩時(shí)刻的全波場(chǎng),采用變階數(shù)有限差分算法進(jìn)行波場(chǎng)重建,但由于缺少準(zhǔn)確的邊界波場(chǎng),這種方法在波場(chǎng)反傳時(shí)會(huì)引入誤差。CLAPP[17]提出了一種吸收邊界條件下的波場(chǎng)存儲(chǔ)策略,根據(jù)有限差分階數(shù)保存若干層波場(chǎng),可保證重建波場(chǎng)無(wú)精度損耗。但是該方法需進(jìn)行頻繁的數(shù)據(jù)交換,尤其是在對(duì)逆時(shí)偏移算法應(yīng)用GPU加速時(shí)[18-20],計(jì)算硬件之間頻繁的數(shù)據(jù)傳輸會(huì)導(dǎo)致較大的時(shí)間延遲,不利于提高計(jì)算效率;王保利等[21]針對(duì)該方法提出了一種有效邊界存儲(chǔ)策略,并將其應(yīng)用于GPU加速平臺(tái),引入檢查點(diǎn)思想,進(jìn)一步降低了數(shù)據(jù)存儲(chǔ)量,結(jié)合GPU高性能計(jì)算的優(yōu)勢(shì),在保證波場(chǎng)重建精度的同時(shí),計(jì)算耗時(shí)也得到良好的控制。

吸收邊界條件雖然可以壓制邊界反射干擾,但往往需要進(jìn)行大量的波場(chǎng)存儲(chǔ),CLAPP[22]提出的隨機(jī)邊界條件,在速度場(chǎng)外圍引入隨機(jī)速度層,降低了邊界反射的相關(guān)性,從而以隨機(jī)噪聲的形式保存了波場(chǎng)能量,對(duì)原方程進(jìn)行逆向運(yùn)算,即可獲得歷史波場(chǎng),僅需保存最后兩時(shí)刻波場(chǎng)切片,所需存儲(chǔ)量較少。然而,隨機(jī)噪聲是初至波場(chǎng)直接經(jīng)由隨機(jī)速度層反射而來(lái),其能量往往與有效反射波場(chǎng)的能量處于同一量級(jí),會(huì)在偏移剖面上產(chǎn)生較大的影響。SHEN等[23]提出一種改進(jìn)的隨機(jī)邊界條件,對(duì)低頻波場(chǎng)的隨機(jī)散射效率更高,更有利于進(jìn)行波形反演計(jì)算,但仍無(wú)法解決隨機(jī)反射的干擾問(wèn)題。綜上所述,目前關(guān)于逆時(shí)偏移的邊界條件和存儲(chǔ)策略的研究均為針對(duì)解決邊界噪聲和波場(chǎng)存儲(chǔ)量二者當(dāng)中的一個(gè)問(wèn)題,并沒(méi)有一個(gè)更靈活的折中方案。

本文方法綜合運(yùn)用PML邊界條件和隨機(jī)邊界條件,旨在保證逆時(shí)偏移成像精度的基礎(chǔ)上,降低波場(chǎng)存儲(chǔ)空間。邊界條件方面,在真實(shí)波場(chǎng)外圍依次加入衰減層和隨機(jī)速度層,降低隨機(jī)反射波振幅,從而減少隨機(jī)噪聲對(duì)逆時(shí)偏移結(jié)果的影響。存儲(chǔ)策略方面,與隨機(jī)邊界條件的存儲(chǔ)策略類似,保存最后時(shí)刻波場(chǎng),采用波動(dòng)方程的逆向運(yùn)算重建激發(fā)點(diǎn)波場(chǎng)信息,進(jìn)而與接收點(diǎn)波場(chǎng)完成互相關(guān)計(jì)算,實(shí)現(xiàn)逆時(shí)偏移成像。在此基礎(chǔ)上,本文方法需設(shè)置檢查點(diǎn),保存臨時(shí)波場(chǎng),以降低波場(chǎng)重建過(guò)程中的累積誤差。最后,對(duì)鹽丘模型和大偏移距模型進(jìn)行逆時(shí)偏移試算,驗(yàn)證了本文方法的有效性。

1 衰減隨機(jī)邊界條件

根據(jù)逆時(shí)偏移的算法需求,在進(jìn)行地震數(shù)據(jù)偏移成像過(guò)程中,需獲得歷史時(shí)刻激發(fā)點(diǎn)波場(chǎng)信息,但保存全部激發(fā)點(diǎn)波場(chǎng)信息會(huì)消耗巨大的存儲(chǔ)空間[18],為節(jié)約數(shù)據(jù)存儲(chǔ)空間,研究人員常使用CLAPP提出的隨機(jī)邊界條件[22],該方法在有效速度場(chǎng)外加入隨機(jī)速度場(chǎng),不進(jìn)行邊界處波場(chǎng)能量的吸收,然后根據(jù)波動(dòng)方程進(jìn)行反向運(yùn)算即可構(gòu)建歷史波場(chǎng)。但是由于初至波場(chǎng)能量較強(qiáng),入射至隨機(jī)邊界層后,產(chǎn)生的隨機(jī)反射會(huì)形成強(qiáng)能量干擾,往往會(huì)在偏移剖面上形成較強(qiáng)的虛假同相軸。我們結(jié)合隨機(jī)邊界條件思想,提出一種速度場(chǎng)的構(gòu)建方法:首先在有效速度場(chǎng)(圖1 中區(qū)域a)外進(jìn)行速度等效拉伸(圖1中區(qū)域b),然后在外層加入隨機(jī)速度場(chǎng)(圖1中區(qū)域c)。逆時(shí)偏移算法中,后續(xù)的波場(chǎng)模擬計(jì)算,即采用上述方法構(gòu)建的速度場(chǎng)。

圖1 隨機(jī)邊界速度場(chǎng)

逆時(shí)偏移算法可采用聲波波動(dòng)方程進(jìn)行波場(chǎng)模擬,在二維空間內(nèi),假定密度為常數(shù),通常采用(1)式所示的聲波方程進(jìn)行數(shù)值計(jì)算:

(1)

式中:P(x,z,t)為地震波場(chǎng);V(x,z)為地下介質(zhì)速度;x和z分別為空間水平方向和垂直方向的坐標(biāo);t為時(shí)間。

本文采用PML吸收方程進(jìn)行波場(chǎng)計(jì)算[12]:

(2)

式中:A(x,z)為衰減系數(shù)矩陣,衰減系數(shù)值隨邊界層外法線方向呈遞增趨勢(shì)分布,可通過(guò)線性函數(shù)、指數(shù)函數(shù)和三角函數(shù)等方法進(jìn)行構(gòu)建[10-12]。當(dāng)A(x,z)為0時(shí),方程(2)可化簡(jiǎn)為(1)式所示的聲波波動(dòng)方程。

本文方法在速度等效拉伸區(qū)域(圖1中區(qū)域b)構(gòu)建衰減系數(shù),當(dāng)波場(chǎng)入射至此區(qū)域時(shí),產(chǎn)生能量衰減,其余區(qū)域衰減系數(shù)設(shè)置為0,無(wú)振幅衰減,即(2)式中A(x,z)項(xiàng)為0,波場(chǎng)的傳播遵循(1)式。

本文采用有限差分?jǐn)?shù)值方法對(duì)方程(2)進(jìn)行離散,整理可得如下形式:

(3)

應(yīng)用方程(3)可進(jìn)行激發(fā)點(diǎn)波場(chǎng)的正向傳播模擬,且當(dāng)?shù)卣鸩ㄟM(jìn)入邊界吸收層后,其振幅會(huì)逐漸衰減。

2 波場(chǎng)重建

(4)

應(yīng)用方程(4)可實(shí)現(xiàn)激發(fā)點(diǎn)波場(chǎng)的逆時(shí)傳播,當(dāng)衰減系數(shù)di,j為0時(shí),激發(fā)點(diǎn)波場(chǎng)無(wú)振幅衰減,當(dāng)di,j為非零值時(shí),即在正傳波場(chǎng)的衰減區(qū)域,地震波振幅會(huì)逐漸增強(qiáng)。

對(duì)利用方程(3)計(jì)算得到的正傳波場(chǎng),僅保存最后兩時(shí)刻激發(fā)點(diǎn)波場(chǎng),由方程(4)進(jìn)行激發(fā)點(diǎn)波場(chǎng)的逆時(shí)重建,再對(duì)其同接收點(diǎn)波場(chǎng)應(yīng)用互相關(guān)成像條件,即可完成逆時(shí)偏移成像。

為驗(yàn)證本文方法對(duì)隨機(jī)噪聲的衰減作用,分別利用隨機(jī)邊界條件和衰減隨機(jī)邊界條件進(jìn)行地震波場(chǎng)的正傳模擬和重建實(shí)驗(yàn)。模型大小為301×301個(gè)網(wǎng)格點(diǎn),網(wǎng)格間距為8m,設(shè)置50層邊界,其中內(nèi)側(cè)40層做速度等效拉伸,最外側(cè)的10層為隨機(jī)速度層。采用衰減隨機(jī)邊界條件時(shí),邊界層內(nèi)側(cè)40層為吸收層,設(shè)置衰減因子。假定模型為3500m/s的勻速介質(zhì),震源位于模型的中央位置,震源子波為主頻30Hz的雷克子波,時(shí)間采樣率為0.5ms,波場(chǎng)傳播最大時(shí)刻為0.55s,采用兩種邊界條件計(jì)算的地震波場(chǎng)在0.55s的波場(chǎng)快照如圖2所示。

抽取1.2km深度處(圖2中虛線位置)波場(chǎng)數(shù)據(jù)可見(jiàn),假定地震子波最大振幅為10000時(shí),在有效波場(chǎng)范圍內(nèi)(0.4～2.0km),采用隨機(jī)邊界條件在虛線處產(chǎn)生的隨機(jī)噪聲最大振幅值為2588,如圖3中虛線所示;采用衰減隨機(jī)邊界條件,由于隨機(jī)噪聲需經(jīng)過(guò)衰減層傳播,進(jìn)入有效波場(chǎng)區(qū)域的隨機(jī)噪聲得以衰減,最大振幅值為669,如圖3中實(shí)線所示。由此可知,采用衰減隨機(jī)邊界條件,經(jīng)過(guò)衰減層后,隨機(jī)反射噪聲被有效壓制,振幅強(qiáng)度下降明顯。

為驗(yàn)證本文方法波場(chǎng)重建的準(zhǔn)確性,進(jìn)行了地震波場(chǎng)的重建實(shí)驗(yàn),并在同參數(shù)條件下,與隨機(jī)邊界條件的波場(chǎng)重建結(jié)果進(jìn)行了對(duì)比。圖4a和圖4b分別為采用衰減隨機(jī)邊界條件在0.25s時(shí)刻波場(chǎng)快照和波場(chǎng)傳至0.55s后,反傳至0.25s時(shí)的波場(chǎng)快照,可以看出,正、反傳波場(chǎng)無(wú)明顯差異,抽取1.2km深度處的波場(chǎng)數(shù)據(jù)進(jìn)行對(duì)比(圖5)可以看出,二者振幅無(wú)明顯差異。

圖2 不同邊界條件波場(chǎng)快照對(duì)比(t=0.55s)a 隨機(jī)邊界條件; b 衰減隨機(jī)邊界條件

圖3 隨機(jī)邊界條件與衰減隨機(jī)邊界條件隨機(jī)噪聲振幅對(duì)比

圖4 衰減隨機(jī)邊界條件正傳與反傳波場(chǎng)快照對(duì)比(t=0.25s)a 正傳波場(chǎng)快照; b 反傳波場(chǎng)快照

圖6a和圖6b分別為隨機(jī)邊界條件下0.25s時(shí)刻波場(chǎng)快照和波場(chǎng)傳至0.55s后,再反傳至0.25s時(shí)的波場(chǎng)快照。由圖6可見(jiàn),二者并無(wú)明顯差別。進(jìn)一步抽取1.2km處的波場(chǎng)數(shù)據(jù)進(jìn)行對(duì)比,如圖7所示,二者振幅大小相當(dāng)。

圖5 圖4a和圖4b所示波場(chǎng)在1.2km處振幅對(duì)比

圖6 隨機(jī)邊界條件正傳與反傳波場(chǎng)快照對(duì)比(0.25s)a 正傳波場(chǎng)快照; b 反傳波場(chǎng)快照

圖7 圖6a和圖6b所示波場(chǎng)在1.2km處振幅對(duì)比

為了進(jìn)行誤差的量化對(duì)比,分別對(duì)兩種邊界條件的正、反傳波場(chǎng)進(jìn)行了差剖面測(cè)試。圖8a為采用隨機(jī)邊界條件在0.25s時(shí)刻正、反傳波場(chǎng)的差剖面,最大誤差值為0.023,圖8b為采用衰減隨機(jī)邊界條件在0.25s時(shí)刻正、反傳波場(chǎng)的差剖面,最大誤差值為0.310?？梢钥闯?應(yīng)用衰減隨機(jī)邊界條件進(jìn)行波場(chǎng)恢復(fù)所產(chǎn)生的誤差值雖然較隨機(jī)邊界條件更大,但誤差主要分布在衰減層區(qū)域，對(duì)有效成像區(qū)域的影響可忽略。主要原因在于:當(dāng)波場(chǎng)進(jìn)入衰減區(qū)后,隨著波場(chǎng)的傳播,振幅衰減程度逐漸增大,采用本文方法進(jìn)行波場(chǎng)重建時(shí),會(huì)有部分波場(chǎng)由于振幅過(guò)小而無(wú)法完全得到保真恢復(fù),進(jìn)而產(chǎn)生波場(chǎng)畸變[17]。本文方法雖然在衰減區(qū)引入了部分誤差，但與對(duì)隨機(jī)噪聲的衰減效果相比，波場(chǎng)的整體噪聲水平有所下降。

圖8 不同邊界條件正傳與反傳波場(chǎng)誤差對(duì)比(t=0.25s)a 隨機(jī)邊界條件; b 衰減隨機(jī)邊界條件

保存波場(chǎng)的時(shí)間周期的選取,受衰減系數(shù)影響,可根據(jù)正傳波場(chǎng)與反傳波場(chǎng)之間的一致性選擇波場(chǎng)保存周期。

本文采用公式(5)進(jìn)行波場(chǎng)一致性估算:

(5)

式中:pf(i,j)為正傳波場(chǎng);pb(i,j)為反傳重建波場(chǎng);I,J為空間水平、垂直方向網(wǎng)格點(diǎn)數(shù);i,j為空間坐標(biāo);Eτ為兩種波場(chǎng)的一致性。公式(5)實(shí)質(zhì)上是估算每個(gè)網(wǎng)格點(diǎn)上正傳波場(chǎng)、反傳波場(chǎng)的能量差異與正傳波場(chǎng)能量的平均比值,可間接反映誤差對(duì)真實(shí)波場(chǎng)的影響程度。

圖9顯示了當(dāng)Eτ不超過(guò)0.2時(shí),波場(chǎng)迭代時(shí)間隨衰減系數(shù)的變化趨勢(shì)。當(dāng)衰減矩陣中最大衰減系數(shù)小于30時(shí),由于衰減程度較弱,振幅恢復(fù)較準(zhǔn)確,更接近于隨機(jī)邊界條件,波場(chǎng)持續(xù)時(shí)間均接近于本次實(shí)驗(yàn)最大采樣時(shí)間。隨著衰減系數(shù)的增大,重建波場(chǎng)在滿足精度要求的基礎(chǔ)上，迭代時(shí)間下降?？紤]到實(shí)際應(yīng)用時(shí)對(duì)噪聲影響和波場(chǎng)存儲(chǔ)量的需求,建議最大衰減系數(shù)選取20～30,既可保證對(duì)隨機(jī)噪聲的衰減,又能增加迭代時(shí)間，減少檢查點(diǎn)個(gè)數(shù)，控制波場(chǎng)存儲(chǔ)量。

3 存儲(chǔ)量分析

波場(chǎng)存儲(chǔ)量和計(jì)算量是評(píng)判逆時(shí)偏移邊界條件和存儲(chǔ)策略的重要指標(biāo),本文對(duì)比分析了隨機(jī)邊界條件、衰減隨機(jī)邊界條件和PML吸收邊界條件的計(jì)算量和存儲(chǔ)量。

如前所述，索緒爾從未提及語(yǔ)言本體論，自然不可能指出本體及其屬性是其理論推理的基點(diǎn)，但他卻將“任意性作為構(gòu)造理論的第一原理或公理，其他論斷則由此演繹出來(lái)” [35]21。那么，任意性和本體是什么關(guān)系呢？“任意性原則之所以是‘頭等重要的’，在于它集中地反映了索緒爾語(yǔ)言學(xué)的本體論特征，它充分地說(shuō)明了語(yǔ)言實(shí)體的本原性” [36]54。這也就是說(shuō)，任意性是本體的根本屬性。索緒爾雖然沒(méi)有明確提及本體論，但他特別強(qiáng)調(diào)任意性，并把它高調(diào)定位為語(yǔ)言學(xué)的第一原則，實(shí)質(zhì)上也就確立了本體論的基礎(chǔ)地位。

設(shè)定模型水平方向網(wǎng)格點(diǎn)數(shù)為Nx,垂直方向的網(wǎng)格點(diǎn)數(shù)為Nz,邊界層數(shù)為L(zhǎng),空間差分階數(shù)為2K,波場(chǎng)傳播在時(shí)間方向的采樣點(diǎn)數(shù)為T,每隔C個(gè)時(shí)間點(diǎn)進(jìn)行一次臨時(shí)波場(chǎng)存儲(chǔ),為便于統(tǒng)計(jì),假設(shè)Nx=Nz=N，則二維空間內(nèi),隨機(jī)邊界條件的存儲(chǔ)需求為2(N+2L)2·4字節(jié),衰減隨機(jī)邊界條件的存儲(chǔ)需求為2T·(N+2L)2·4/C字節(jié),采用PML吸收邊界條件及有效邊界存儲(chǔ)策略的存儲(chǔ)需求為4KNT·4字節(jié),易知通常情況下PML吸收邊界條件存儲(chǔ)量占用最大,衰減隨機(jī)邊界條件次之,隨機(jī)邊界條件的存儲(chǔ)需求最小。

在波場(chǎng)反傳過(guò)程中,需要對(duì)整體波場(chǎng)進(jìn)行計(jì)算,但由于衰減隨機(jī)邊界條件采用帶有衰減項(xiàng)的波動(dòng)方程,因此在網(wǎng)格點(diǎn)數(shù)相同的情況下,衰減隨機(jī)邊界條件的計(jì)算量略大于隨機(jī)邊界條件,由于PML吸收邊界條件在波場(chǎng)反傳過(guò)程中并非全波場(chǎng)計(jì)算,因此無(wú)法進(jìn)行計(jì)算量的量化對(duì)比。

4 逆時(shí)偏移算例

為驗(yàn)證本文方法在逆時(shí)偏移應(yīng)用中的有效性,分別采用隨機(jī)邊界條件、PML邊界條件(有效邊界存儲(chǔ)策略)和衰減隨機(jī)邊界條件進(jìn)行逆時(shí)偏移計(jì)算。

測(cè)試的模型水平方向網(wǎng)格點(diǎn)數(shù)為800、垂直方向網(wǎng)格點(diǎn)數(shù)為600,空間網(wǎng)格間距為10m,激發(fā)100炮,炮點(diǎn)位于0～8000m,炮間距80m,接收點(diǎn)位于地表0～8000m,共800個(gè)檢波器,道間距10m,時(shí)間采樣間隔為0.8ms,采樣點(diǎn)數(shù)為6000,共記錄4.8s,成像條件為震源歸一化互相關(guān)成像條件。

4.1 隨機(jī)邊界條件

考慮到常規(guī)聲波波動(dòng)方程無(wú)振幅能量衰減的特性,隨機(jī)邊界條件[22]在速度場(chǎng)外圍加入隨機(jī)速度層,當(dāng)波場(chǎng)傳播至隨機(jī)速度層后,并不會(huì)產(chǎn)生衰減,即波場(chǎng)能量并無(wú)損失,而是以隨機(jī)散射噪聲的形式得以保存。因此,理論上對(duì)隨機(jī)邊界條件下的任意時(shí)刻波場(chǎng),采用波動(dòng)方程的逆向運(yùn)算,即可獲得歷史波場(chǎng)。在逆時(shí)偏移應(yīng)用中,保存最大時(shí)刻的激發(fā)點(diǎn)波場(chǎng)信息,根據(jù)原波動(dòng)方程即可實(shí)現(xiàn)激發(fā)點(diǎn)波場(chǎng)的逆時(shí)反傳,從而與檢波點(diǎn)波場(chǎng)進(jìn)行相關(guān)成像。

本文測(cè)試設(shè)定邊界層數(shù)為100,外側(cè)30層為隨機(jī)速度層,內(nèi)側(cè)70層根據(jù)模型速度場(chǎng)做等效拉伸。存儲(chǔ)策略為保存最后兩時(shí)刻激發(fā)點(diǎn)波場(chǎng),同時(shí)進(jìn)行激發(fā)點(diǎn)波場(chǎng)反傳和接收點(diǎn)波場(chǎng)正傳,完成逆時(shí)偏移成像。

4.2 PML吸收邊界條件(有效邊界存儲(chǔ)策略)

PML吸收邊界條件[6]通過(guò)修改波動(dòng)方程,引入衰減項(xiàng),當(dāng)波場(chǎng)進(jìn)入邊界層后,會(huì)由內(nèi)向外逐漸衰減,且不會(huì)產(chǎn)生劇烈反射。該邊界條件可有效解決波場(chǎng)的邊界反射問(wèn)題,但由于在邊界層內(nèi),波場(chǎng)振幅被衰減,所以常規(guī)的波動(dòng)方程逆向運(yùn)算無(wú)法精確重建激發(fā)點(diǎn)波場(chǎng)。有效邊界存儲(chǔ)策略是根據(jù)有限差分法中空間差分的計(jì)算特性,在激發(fā)點(diǎn)波場(chǎng)正傳時(shí),保存一定層數(shù)的內(nèi)部波場(chǎng),從而進(jìn)行波場(chǎng)逆推,理論上,層數(shù)為空間差分階數(shù)的一半時(shí),即可實(shí)現(xiàn)波場(chǎng)的精確重建[21]。

本文測(cè)試設(shè)定邊界層數(shù)為100,邊界層內(nèi)為PML吸收邊界條件,在激發(fā)點(diǎn)波場(chǎng)正傳過(guò)程中,根據(jù)本文采用的10階空間差分,保存5層邊界波場(chǎng)及最后兩時(shí)刻波場(chǎng);在激發(fā)點(diǎn)波場(chǎng)反傳過(guò)程中,每時(shí)刻均在有效邊界層中進(jìn)行波場(chǎng)加載,保證內(nèi)部有效區(qū)域內(nèi)的波場(chǎng)信息完整。理論上,該方法對(duì)有效區(qū)域內(nèi)激發(fā)點(diǎn)波場(chǎng)的恢復(fù)精度最高,等同于保存全部波場(chǎng)所獲得的計(jì)算精度,因此在本文中,將采用該方法的逆時(shí)偏移剖面作為評(píng)判偏移效果的參考。

4.3 衰減隨機(jī)邊界條件(檢查點(diǎn)存儲(chǔ)策略)

邊界層數(shù)為100,外側(cè)30層為隨機(jī)速度層,對(duì)波場(chǎng)無(wú)衰減;內(nèi)側(cè)70層根據(jù)模型速度場(chǎng)做等效拉伸且設(shè)置PML衰減層。存儲(chǔ)策略方面,根據(jù)檢查點(diǎn)間隔,每隔一定時(shí)間進(jìn)行一次波場(chǎng)保存,本次測(cè)試檢查點(diǎn)設(shè)置為600,即每隔600個(gè)時(shí)間點(diǎn)保存一次波場(chǎng)。在運(yùn)用方程(4)進(jìn)行激發(fā)點(diǎn)波場(chǎng)反傳過(guò)程中,當(dāng)傳播至檢查點(diǎn)時(shí)刻時(shí),用保存的波場(chǎng)進(jìn)行波場(chǎng)替換,然后繼續(xù)反傳。

圖10給出了速度模型及采用不同邊界條件得到的偏移剖面。由圖10可見(jiàn),采用隨機(jī)邊界條件的偏移剖面,受隨機(jī)噪聲影響較大,在中淺層由于入射至邊界層的波場(chǎng)能量較強(qiáng),因此隨機(jī)噪聲影響較為明顯,出現(xiàn)了雜亂的同相軸,影響對(duì)真實(shí)構(gòu)造形態(tài)的判斷;采用衰減隨機(jī)邊界條件與PML邊界條件的偏移剖面無(wú)明顯差異。在本文逆時(shí)偏移測(cè)試中,PML吸收邊界條件波場(chǎng)存儲(chǔ)量約為320.4MB,衰減隨機(jī)邊界條件(檢查點(diǎn)間隔為600)波場(chǎng)存儲(chǔ)量?jī)H約為54.9MB,較前者節(jié)約了5倍存儲(chǔ)量,隨機(jī)邊界條件波場(chǎng)存儲(chǔ)量占用最小,約為6.1MB,實(shí)驗(yàn)結(jié)果與本文第3部分結(jié)論一致,3種邊界條件的波場(chǎng)存儲(chǔ)量及噪聲影響情況如表1所示。

圖10 速度模型及采用不同邊界條件得到的逆時(shí)偏移剖面a 速度模型; b 隨機(jī)邊界條件; c PML吸收邊界條件; d 衰減隨機(jī)邊界條件

邊界條件及存儲(chǔ)策略波場(chǎng)存儲(chǔ)量噪聲影響隨機(jī)邊界條件小大衰減隨機(jī)邊界條件中小PML吸收邊界條件(有效邊界存儲(chǔ)策略)大小

4.4 大偏移距模型測(cè)試

在橫縱比較大的模型中,為完成遠(yuǎn)偏移距炮集數(shù)據(jù)的成像計(jì)算,通常采用較大的采樣時(shí)長(zhǎng),即逆時(shí)偏移的時(shí)間采樣時(shí)長(zhǎng)遠(yuǎn)超過(guò)近偏移距數(shù)據(jù)的成像時(shí)長(zhǎng),從而能更全面地測(cè)試近偏移距區(qū)域邊界條件的應(yīng)用效果。為進(jìn)一步檢測(cè)衰減隨機(jī)邊界條件的適用性,本文進(jìn)行了大偏移距模型的逆時(shí)偏移測(cè)試。測(cè)試模型水平方向網(wǎng)格點(diǎn)數(shù)為1501,垂直方向網(wǎng)格點(diǎn)數(shù)為151,橫縱比約為10∶1,邊界層數(shù)為100,3種邊界條件的設(shè)置與4.1到4.3節(jié)一致。水平和垂直方向的空間網(wǎng)格間距為10m,精確速度模型如圖11a所示,逆時(shí)偏移所用的平滑速度模型如圖11b所示。激發(fā)11炮,炮點(diǎn)位于4000～10000m,炮間距600m,接收點(diǎn)位于地表0～15000m,共1501個(gè)檢波器,道間距10m,時(shí)間采樣間隔為1.0ms,采樣點(diǎn)數(shù)為5000,共記錄5.0s,成像條件為互相關(guān)成像條件。

圖11c到圖11e分別為隨機(jī)邊界條件、PML邊界條件(有效邊界存儲(chǔ)策略)和衰減隨機(jī)邊界條件(檢查點(diǎn)間隔為600)的逆時(shí)偏移測(cè)試結(jié)果。對(duì)比可知,在該模型測(cè)試中,采用隨機(jī)邊界條件時(shí),逆時(shí)偏移剖面仍存在較大的隨機(jī)噪聲干擾;采用PML邊界條件和本文提出的衰減隨機(jī)邊界條件的結(jié)果相近,無(wú)隨機(jī)噪聲干擾。

由于遠(yuǎn)偏移距數(shù)據(jù)受覆蓋次數(shù)限制和折射波影響,逆時(shí)偏移效果不盡理想,因此本文在逆時(shí)偏移疊加成像過(guò)程中,對(duì)遠(yuǎn)偏移距的成像結(jié)果進(jìn)行了衰減,分別采用上述3種邊界條件得到的測(cè)試結(jié)果如圖11f到圖11h所示。對(duì)比可知,經(jīng)遠(yuǎn)偏移距成像數(shù)據(jù)衰減后,采用3種邊界條件的逆時(shí)偏移結(jié)果得到明顯改善,同時(shí)也可發(fā)現(xiàn),采用隨機(jī)邊界條件仍存在隨機(jī)噪聲干擾;采用PML邊界條件和衰減隨機(jī)邊界條件得到的結(jié)果則不存在隨機(jī)噪聲干擾。

圖11 速度模型及不同邊界條件逆時(shí)偏移和遠(yuǎn)偏移距衰減逆時(shí)偏移剖面a 精確速度模型; b 平滑速度模型; c 隨機(jī)邊界條件逆時(shí)偏移剖面; d PML邊界條件逆時(shí)偏移剖面; e 衰減隨機(jī)邊界條件逆時(shí)偏移剖面; f 隨機(jī)邊界條件遠(yuǎn)偏移距衰減逆時(shí)偏移剖面; g PML吸收邊界條件遠(yuǎn)偏移距衰減逆時(shí)偏移剖面; h 衰減隨機(jī)邊界條件遠(yuǎn)偏移距衰減逆時(shí)偏移剖面

本節(jié)測(cè)試中,隨機(jī)邊界條件存儲(chǔ)量占用約為2.27MB,PML吸收邊界條件波場(chǎng)存儲(chǔ)量約為315.09MB,衰減隨機(jī)邊界條件(檢查點(diǎn)點(diǎn)數(shù)為600)波場(chǎng)存儲(chǔ)量約為38.0MB,進(jìn)一步驗(yàn)證了本文方法在保證成像精度的基礎(chǔ)上,降低了波場(chǎng)存儲(chǔ)量。

5 結(jié)論與展望

本文方法綜合利用了隨機(jī)邊界條件對(duì)反射波場(chǎng)隨機(jī)化的特性和吸收邊界條件對(duì)波場(chǎng)振幅的衰減特性,在隨機(jī)速度層內(nèi)側(cè)加入波場(chǎng)衰減層,對(duì)隨機(jī)反射波場(chǎng)進(jìn)行衰減,從而降低隨機(jī)噪聲干擾,并且通過(guò)對(duì)波動(dòng)方程的逆向運(yùn)算重建激發(fā)點(diǎn)波場(chǎng),降低逆時(shí)偏移過(guò)程中波場(chǎng)的存儲(chǔ)量。對(duì)于波場(chǎng)重建過(guò)程中產(chǎn)生的波形畸變,結(jié)合檢查點(diǎn)波場(chǎng)存儲(chǔ)策略可將其有效壓制。模型數(shù)據(jù)逆時(shí)偏移試算表明,本文方法能夠在犧牲較小存儲(chǔ)效率的基礎(chǔ)上,有效壓制隨機(jī)噪聲的影響。

相比于隨機(jī)邊界條件方法,本文方法存在一定不足,例如計(jì)算量略大,而且因需要在檢查點(diǎn)時(shí)間處進(jìn)行波場(chǎng)保存,會(huì)帶來(lái)額外的波場(chǎng)存儲(chǔ)量和數(shù)據(jù)傳輸延遲。

相比于PML吸收邊界條件及其存儲(chǔ)策略,本文方法波場(chǎng)存儲(chǔ)需求較小。而且由于具有較高的成像精度,噪聲引入少,因此本文方法同樣適用于最小二乘逆時(shí)偏移及全波形反演等算法中的梯度計(jì)算,這也是作者下一步的研究工作。

[1] BAYSAL E,KOSLOFF D D,SHERWOOD J W C.Reverse time migration[J].Geophysics,1983,48(11):1514-1524

[2] 楊勤勇,段心標(biāo).逆時(shí)偏移技術(shù)發(fā)展現(xiàn)狀與趨勢(shì)[J].石油物探,2010,49(1):92-98 YANG Q Y,DUAN X B.Development status and trend of reverse time migration[J].Geophysical Prospecting for Petroleum,2010,49(1):92-98

[3] 畢麗飛,秦寧,楊曉東,等.彈性多波高斯束逆時(shí)偏移方法[J].石油物探,2015,54(1):64-70 BI L F,QIN N,YANG X D,et al.Gauss beam reverse time migration method for elastic multiple wave[J].Geophysical Prospecting for Petroleum,2015,54(1):64-70

[4] 郭書娟,馬方正,段心標(biāo),等.最小二乘逆時(shí)偏移成像方法的實(shí)現(xiàn)與應(yīng)用研究[J].石油物探,2015,54(3):301-308 GUO S J,MA F Z,DUAN X B,et al.Research of least-squares reverse-time migration imaging method and its application[J].Geophysical Prospecting for Petroleum,2015,54(3):301-308

[5] 段心標(biāo),王華忠,白英哲,等.基于GPU的三維起伏地表單程波疊前深度偏移[J].石油物探,2016,55(2):223-230 DUAN X B,WANG H Z,BAI Y Z,et al.3D one-way wave equation prestack depth migration from topography based on the acceleration of GPU[J].Geophysical Prospecting for Petroleum,2016,55(2):223-230

[6] BéRENGER J P.A perfectly matched layer for the absorption of electromagnetic waves[J].Journal of Computational Physics,1994,114(2):185-200

[7] WENG C C,WEEDON W H.A 3D perfectly matched medium from modified maxwell’s equations with stretched coordinates[J].Microwave & Optical Technology Letters,1994,7(13):599-604

[8] 吳國(guó)忱,梁鍇.VTI介質(zhì)準(zhǔn)P波頻率空間域組合邊界條件研究[J].石油物探,2005,44(4):301-307 WU G C,LIANG K.Combined boundary conditions of quasi-P wave within frequency-space domain in VTI media[J].Geophysical Prospecting for Petroleum,2005,44(4):301-307

[9] 單啟銅,樂(lè)友喜.PML邊界條件下二維粘彈性介質(zhì)波場(chǎng)模擬[J].石油物探,2007,46(2):126-130 SHAN Q T,YUE Y X.Wavefield simulation of 2-D viscoelastic medium in perfectly matched layer boundary[J].Geophysical Prospecting for Petroleum,2007,46(2):126-130

[10] 陳可洋.完全匹配層吸收邊界條件研究[J].石油物探,2010,49(5):472-477 CHEN K Y.Study on perfectly matched layer absorbing boundary condition[J].Geophysical Prospecting for Petroleum,2010,49(5):472-477

[11] 高剛,賀振華,黃德濟(jì),等.完全匹配層人工邊界條件中的衰減因子分析[J].石油物探,2011,50(5):430-433 GAO G,HE Z H,HUANG D J,et al.Analysis on attenuation factor in the processing of artificial boundary conditions of PML[J].Geophysical Prospecting for Petroleum,2011,50(5):430-433

[12] 王維紅,柯璇,裴江云.完全匹配層吸收邊界條件應(yīng)用研究[J].地球物理學(xué)進(jìn)展,2013,28(5):2508-2514 WANG W H,KE X,PEI J Y.Application investigation of perfectly matched layer absorbing boundary condition[J].Progress in Geophysics,2013,28(5):2508-2514

[13] 張衡,劉洪,李博,等.VTI介質(zhì)聲波方程非分裂式PML吸收邊界條件研究[J].石油物探,2016,55(6):781-792 ZHANG H,LIU H,LI B,et al.The research on unsplit PML absorbing boundary conditions of acoustic equation for VTI media[J].Geophysical Prospecting for Petroleum,2016,55(6):781-792

[14] LIU Y,SEN M K.A hybrid scheme for absorbing edge reflections in numerical modeling of wave propagation[J].Geophysics,2010,75(2):A1-A6

[15] SYMES W W.Reverse time migration with optimal checkpointing[J].Geophysics,2007,72(5):213-221

[16] FENG B,WANG H Z.Reverse time migration with source wavefield reconstruction strategy[J].Journal of Geophysics & Engineering,2012,9(1):69-74

[17] CLAPP R G.Reverse time migration:saving the boundaries[R].Stanford Exploration Project,2008:136-144

[18] 李博,劉紅偉,劉國(guó)峰,等.地震疊前逆時(shí)偏移算法的CPU/GPU實(shí)施對(duì)策[J].地球物理學(xué)報(bào),2010,53(12):2938-2943 LI B,LIU H W,LIU G F,et al.Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J].Chinese Journal of Geophysics,2010,53(12):2938-2943

[19] 劉紅偉,李博,劉洪,等.地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn)[J].地球物理學(xué)報(bào),2010,53(7):1725-1733 LIU H W,LI B,LIU H,et al.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J].Chinese Journal of Geophysics,2010,53(7):1725-1733

[20] 劉紅偉,劉洪,鄒振,等.地震疊前逆時(shí)偏移中的去噪與存儲(chǔ)[J].地球物理學(xué)報(bào),2010,53(9):2171-2180 LIU H W,LIU H,ZOU Z,et al.The problems of denoise and storage in seismic reverse time migration[J].Chinese Journal of Geophysics,2010,53(9):2171-2180

[21] 王保利,高靜懷,陳文超,等.地震疊前逆時(shí)偏移的有效邊界存儲(chǔ)策略[J].地球物理學(xué)報(bào),2012,55(7):2412-2421 WANG B L,GAO J H,CHEN W C,et al.Efficient boundary storage strategies for seismic reverse time migration[J].Chinese Journal of Geophysics,2012,55(7):2412-2421

[22] CLAPP R G.Reverse time migration with random boundaries[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:2809-2813

[23] SHEN X K,CLAPP R G.Random boundary condition for memory-efficient waveform inversion gradient computation[J].Geophysics,2015,80(6):R351-R359

[24] CLAERBOUT J F.Toward a unified theory of reflector mapping[J].Geophysics,1971,36(3):467-481

(編輯：陳 杰)

A damped random boundary condition for prestack reverse time migration

KE Xuan1,SHI Ying1,2,ZHANG Yingying1,ZHANG Zhen3

(1.EarthScienceCollegeofNortheastPetroleumUniversity,Daqing163318,China;2.HeilongjiangProvincecollegescienceandtechnologyinnovationteam‘faultsdeformation,sealingabilityandfluidmigration’,Daqing163318,China;3.ResearchInstituteofExploration&Development,PetroChinaTarimOilfieldCompany,KorlaXinjiang841000,China)

The boundary conditions and storage requirements of reverse time migration (RTM) are often restricted to each other.The selection of these two parameters directly affect the imaging results and application conditions of RTM.Therefore,it is very important to study the boundary conditions and corresponding storage strategies for RTM.We proposed a damped random boundary condition,in which the energy attenuation layer and the random velocity layer were added outside of the real wavefield in turn.This method could effectively reduce the energy interference from random reflection.The inverse operation of the wave equation could reconstruct the source wavefield.Combined with cross-correlation imaging conditions,we conducted the RTM calculation.The results of comparison and analysis show that,compared with the random boundary condition,at the expense of less storage efficiency,the proposed method could obtain a better migration result,and the accuracy is very close to the results obtained by the perfectly matched layer boundary condition.The correctness of the above conclusion is verified by the example of salt model and big offset model.

reverse time migration (RTM),random boundary,storage,wavefield reconstruction

2016-09-27;改回日期：2017-04-15。

柯璇(1989—),男,博士在讀,主要從事地震波場(chǎng)正演模擬與偏移成像等方面的研究。

石穎(1976—),女,教授,博士生導(dǎo)師,主要從事地震資料處理方面的研究。

國(guó)家自然科學(xué)基金項(xiàng)目(41574117,41474118)、黑龍江省杰出青年科學(xué)基金項(xiàng)目(JC2016006)、大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金項(xiàng)目(LP1509)、黑龍江省自然科學(xué)基金項(xiàng)目(D2015011)和東北石油大學(xué)研究生創(chuàng)新科研項(xiàng)目(YJSCX2016-002NEPU)聯(lián)合資助。

P631

A

1000-1441(2017)04-0523-11

10.3969/j.issn.1000-1441.2017.04.008

This research is financially supported by the National Natural Science Foundation of China (Grant Nos.41574117,41474118),Heilongjiang Province Natural Science Fund for Distinguished Young Scholar (Grant No.JC2016006),the Open Project of State Key Laboratory of Coastal and Offshore Engineering of Dalian University of Technology (Grant No.LP1509),the Natural Science Foundation of Heilongjiang Province (Grant No.D2015011) and the Northeast Petroleum University Innovation Foundation for Postgraduate (Grant No.YJSCX2016-002NEPU).