逆時(shí)偏移計(jì)算中的邊界處理分析及應(yīng)用

胡 昊，劉伊克，常 旭，王一博，杜向東，楊仁虎

1中國(guó)科學(xué)院地質(zhì)與地球物理研究所，北京 100029

2中海石油（中國(guó)）有限公司研究總院，北京 100027

3中國(guó)石油集團(tuán)東方地球物理勘探有限責(zé)任公司，河北涿州 072751

1 引 言

逆時(shí)偏移歷史悠久，早在20世紀(jì)80年代就已經(jīng)提出［1－2］，其成像基本原理是基于反射地震成像演變而來(lái)［3］.逆時(shí)偏移是基于雙程波算子，精度高，不受傾角和速度劇烈變化限制［4］，并且可以利用回轉(zhuǎn)波成像［5］，比常規(guī)偏移更具有成像優(yōu)勢(shì).當(dāng)然，逆時(shí)偏移也有著計(jì)算量巨大，低頻成像噪音，和波場(chǎng)延拓方向不一致帶來(lái)的存儲(chǔ)等問(wèn)題需要解決［6］.最近幾年計(jì)算機(jī)技術(shù)的飛速發(fā)展，讓逆時(shí)偏移重新獲得生命力，它所存在的問(wèn)題也得到了長(zhǎng)足的改善.在處理計(jì)算量的問(wèn)題上，通過(guò)計(jì)算機(jī)的并行計(jì)算［7］，應(yīng)用GPU計(jì)算技術(shù)可以大幅提高計(jì)算效率［8－9］.在處理低頻成像噪音的問(wèn)題上，有在波傳播過(guò)程中壓制產(chǎn)生噪音的方法［10－12］，有對(duì)成像條件做出改進(jìn)的方法［13－14］，有成像后濾波的方法［15－17］等.在解決由波場(chǎng)延拓方向不一致帶來(lái)的波場(chǎng)傳播歷史的存儲(chǔ)問(wèn)題上，有設(shè)置檢查點(diǎn)（Checkpoint）方法［18］，記錄波場(chǎng)邊界信息的吸收邊界方法［19］，以及隨機(jī)邊界方法［20］等.本文著重討論對(duì)比記錄波場(chǎng)邊界的吸收邊界和隨機(jī)邊界這兩種邊界處理方式的計(jì)算效率和成像效果，提出吸收邊界中只記錄邊界范圍內(nèi)未衰減單層波場(chǎng)的方法并驗(yàn)證了其可行性，采用不同的邊界策略應(yīng)用于模型數(shù)據(jù)和實(shí)際資料的處理中.

在疊前逆時(shí)偏移的過(guò)程中，需要將源波場(chǎng)從零時(shí)刻正傳，檢波點(diǎn)波場(chǎng)從最大時(shí)刻逆?zhèn)鳎⒃谕粫r(shí)刻成像，波場(chǎng)延拓方向的不一致導(dǎo)致必須將某一波場(chǎng)的傳播歷史對(duì)另一波場(chǎng)的傳播過(guò)程可見(jiàn)，最直觀的處理方法就是將其中一波場(chǎng)傳播的歷史全部保存，在另一波場(chǎng)傳播過(guò)程中載入.但這種方法對(duì)存儲(chǔ)要求太高，目前的計(jì)算機(jī)硬件無(wú)法承受大型區(qū)域的計(jì)算.針對(duì)波場(chǎng)延拓方向不一致造成的存儲(chǔ)問(wèn)題，目前工業(yè)上主流的處理方法是設(shè)置檢查點(diǎn)方法.設(shè)置檢查點(diǎn)方法，是在波場(chǎng)正傳的中間過(guò)程中設(shè)置若干個(gè)檢查點(diǎn)，將檢查點(diǎn)處的波場(chǎng)記錄下來(lái)，需要某一時(shí)刻點(diǎn)的波場(chǎng)的時(shí)候再通過(guò)最近的檢查點(diǎn)正推即可得到.通過(guò)設(shè)置最優(yōu)檢查點(diǎn)［18］，可以使得獲得源波場(chǎng)傳播歷史的計(jì)算量控制在O（NlogM）的量級(jí)上，其中N是正演波場(chǎng)n步的計(jì)算量，M是設(shè)置的檢查點(diǎn)數(shù).一般情況下，當(dāng)n＝10000，M取36時(shí)，計(jì)算量和內(nèi)存的控制達(dá)到平衡，但計(jì)算量是O（3.1N）量級(jí).在二維問(wèn)題的時(shí)候，檢查點(diǎn)時(shí)刻的波場(chǎng)可以存放在內(nèi)存中而不需要進(jìn)行磁盤I／O操作，但在處理三維問(wèn)題的時(shí)候，就必須將檢查點(diǎn)時(shí)刻的波場(chǎng)存入磁盤或者采取有損壓縮的方法存儲(chǔ).同時(shí)，3倍正演波場(chǎng)的計(jì)算量也明顯增加了計(jì)算成本.

另一個(gè)直觀的方法是將某一波場(chǎng)正傳到時(shí)間最大值后，與另一波場(chǎng)同時(shí)逆?zhèn)?，即可得到兩個(gè)波場(chǎng)在同一時(shí)刻的值.這種方法的關(guān)鍵是如何才能保證波場(chǎng)正傳到最大時(shí)刻后能正確逆?zhèn)骰貋?lái).在正傳的過(guò)程中，為了消除計(jì)算區(qū)域邊界帶來(lái)的反射在成像過(guò)程中形成反射假象，必須要在邊界處做處理，減少邊界反射，而這樣會(huì)帶來(lái)波場(chǎng)能量的損失，波場(chǎng)不能正確逆?zhèn)骰貋?lái).記錄波場(chǎng)邊界信息的吸收邊界方法，就是在波場(chǎng)正傳的同時(shí)，記錄各個(gè)時(shí)刻邊界處未衰減的波場(chǎng)值，到達(dá)最大時(shí)刻后，在逆?zhèn)鞯倪^(guò)程中逐步載入.依據(jù)計(jì)算的可逆性，能保證回傳波場(chǎng)的準(zhǔn)確性.這種方法獲得波場(chǎng)傳播歷史所需的計(jì)算量是正演波場(chǎng)的計(jì)算量的兩倍，邊界處波場(chǎng)的存儲(chǔ)也是一個(gè)值得討論的問(wèn)題.一般來(lái)講，有限差分求解波場(chǎng)方程的精度為n階的時(shí)候，需要記錄n／2層的邊界層的波場(chǎng)值［19］，除去計(jì)算機(jī)數(shù)值計(jì)算的誤差外，可以使得波場(chǎng)完全正確地逆?zhèn)骰貋?lái).這種方法獲得的波場(chǎng)傳播歷史是準(zhǔn)確的，同保存波場(chǎng)傳播歷史的傳統(tǒng)方法獲得的成像效果一致.在處理二維問(wèn)題時(shí)，只需要記錄四個(gè)邊上邊界層的波場(chǎng)信息，所需內(nèi)存不大，但涉及到大區(qū)域三維問(wèn)題的時(shí)候，需存儲(chǔ)計(jì)算區(qū)域的六個(gè)面上的波場(chǎng)信息，會(huì)帶來(lái)一定的存儲(chǔ)問(wèn)題，影響計(jì)算效率.

2009年Clapp提出了隨機(jī)邊界方法［20］.在正常的介質(zhì)周圍加上一層隨機(jī)生成的速度介質(zhì)，使邊界處產(chǎn)生明顯的漫反射，波場(chǎng)傳播到最大時(shí)刻后能正確地逆?zhèn)骰貋?lái)，再通過(guò)多次疊加減輕隨機(jī)邊界帶來(lái)的邊界漫反射影響.此方法使得在獲得波場(chǎng)傳播歷史的計(jì)算過(guò)程中不需要額外存儲(chǔ)任何波場(chǎng)信息，計(jì)算量是兩倍于正演波場(chǎng)的計(jì)算量.但這種方法也有明顯的缺點(diǎn)，即隨機(jī)介質(zhì)造成的漫反射效果需要多次疊加才能減輕，對(duì)觀測(cè)系統(tǒng)提出了更高的要求.邊界層的厚度會(huì)明顯影響到總的計(jì)算量.

本文將簡(jiǎn)單介紹記錄波場(chǎng)邊界信息的吸收邊界方法和隨機(jī)邊界方法，提出吸收邊界只需記錄邊界范圍內(nèi)未衰減單層波場(chǎng)的方法并論證其可行性，并對(duì)比分析選取不同的記錄波場(chǎng)邊界信息的吸收邊界和隨機(jī)邊界的邊界處理策略的成像效果和計(jì)算效率.

2 記錄波場(chǎng)邊界信息的吸收邊界方法和隨機(jī)邊界方法

2.1 吸收邊界方法

本文采用的是交錯(cuò)網(wǎng)格求解聲波方程［21－23］，三維情況下，公式如下所示：

其中，K＝ρV2，采用時(shí)間一階差分，空間n階中心差分格式.正傳的過(guò)程如下：

因此，需要記錄t時(shí)刻的邊界n／2＋1層的Vx，Vy，Vz和P，才能完全逆?zhèn)?

表1 交錯(cuò)網(wǎng)格不同差分階數(shù)的差分系數(shù)［22］Table 1 The difference coefficients of different differential orders for staggered grid［22］

分析差分系數(shù)，以交錯(cuò)網(wǎng)格一階導(dǎo)數(shù)為例.如表1所示，a1值明顯較大，即某點(diǎn)處的導(dǎo)數(shù)是由該點(diǎn)緊鄰的左右兩點(diǎn)主要貢獻(xiàn)的，這種現(xiàn)象在求取規(guī)則網(wǎng)格二階導(dǎo)數(shù)時(shí)也存在.在逆?zhèn)鬟^(guò)程中，波場(chǎng)邊界值是作為邊界條件每個(gè)時(shí)刻均加載入的，因此，若只記錄邊界處未衰減的一層波場(chǎng)，也應(yīng)當(dāng)能保證回傳的波形信息不受破壞，僅能量有所丟失.

如圖1a所示，震源在中心處，波向兩邊傳播，采用的是空間四階交錯(cuò)網(wǎng)格和時(shí)間一階差分形式，邊界處采用海綿吸收處理，記錄各個(gè)時(shí)刻邊界處未衰減單層波場(chǎng)信息，并在逆?zhèn)鬟^(guò)程中載入，同時(shí)記錄正傳經(jīng)過(guò)0.5s和正傳到最大時(shí)刻2s后回傳到0.5s的波形，其中正傳的波形峰值為1，網(wǎng)格數(shù)為400，網(wǎng)格間距dx＝10m，步長(zhǎng)dt＝1ms.正傳的波形用圓圈表示，逆?zhèn)鞯牟ㄐ斡镁€條表示，逆?zhèn)鞯牟ㄐ瓮齻鞯牟ㄐ巫鰵w一化后用圓點(diǎn)表示.圖1b表示為正傳波形，逆?zhèn)鞑ㄐ?，以及逆?zhèn)鞑ㄐ瓮齻鞑ㄐ螝w一化后的頻譜.可以看出正傳到最大時(shí)刻后逆?zhèn)鞯牟ㄐ危ň€條）同正傳的波形（圓圈）在同一時(shí)刻時(shí)波形信息保持非常好，僅僅振幅有所衰弱.同時(shí)，將逆?zhèn)鞯牟ㄐ瓮齻鞯牟ㄐ巫鰵w一化后，可看到其形態(tài)同正傳的波形（圓圈）幾乎重合，表明了逆?zhèn)鞯牟ㄐ瓮齻鞯牟ㄐ蝺H僅是振幅上的衰弱，波形形態(tài)幾乎沒(méi)有發(fā)生改變.

圖1說(shuō)明，采用吸收邊界并只記錄單層波場(chǎng)邊界的方法會(huì)使正傳波場(chǎng)傳播到最大時(shí)刻后回傳過(guò)程中損失部分能量，但每一時(shí)刻都更新邊界處的波場(chǎng)，保證了逆?zhèn)鞑▓?chǎng)只是在振幅上有一些衰減，波形信息依舊保持很好，不影響最終的成像結(jié)果.因此，所需記錄的波場(chǎng)邊界大大減少，存儲(chǔ)問(wèn)題也得到有效改善.此方法能獲得同存儲(chǔ)傳統(tǒng)波場(chǎng)傳播歷史的方法和設(shè)置檢查點(diǎn)方法一樣的成像效果，但只增加一倍波場(chǎng)正演的計(jì)算量.

圖1 正傳波形與記錄單層邊界逆?zhèn)鞑ㄐ蔚膶?duì)比（a）正傳至0.5s波形和回傳至0.5s波形比較.圓圈表示正傳波形，線條表示逆?zhèn)鞑ㄐ危瑘A點(diǎn)表示同正傳波形做歸一化后的逆?zhèn)鞑ㄐ?；（b）正傳至0.5s波形和回傳至0.5s波形的頻譜比較.圓圈表示正傳波形頻譜，線條表示逆?zhèn)鞑ㄐ晤l譜，圓點(diǎn)表示同正傳波形做歸一化后的逆?zhèn)鞑ㄐ晤l譜.Fig.1 The contrast of waveform between forward and backward propagation by recording one layer of boundary（a）The contrast of waveform between forward and backward propagation at the same time（0.5s）.The circle indicates forward propagated waveform，the line indicates backward propagated waveform and the spot indicates normalized backward propagated waveform with forward propagated waveform；（b）The contrast of waveform spectrum between forward and backward propagation at the same time（0.5s）.The circle indicates forward propagated waveform spectrum，the line indicates backward propagated waveform spectrum and thespot indicates normalized backward propagated waveform with forward propagated waveform spectrum.

2.2 隨機(jī)邊界方法

隨機(jī)邊界的目的在于保證波場(chǎng)能夠傳播到最大時(shí)刻后正確回傳并且邊界反射不會(huì)影響正常區(qū)域的成像.因此，控制好邊界漫反射是關(guān)鍵.隨機(jī)層越厚，隨機(jī)效果越好，反射強(qiáng)度越弱，反射的紊亂性越好，對(duì)靠近隨機(jī)邊界處的成像影響越小，但是層數(shù)的增加也將明顯增加總的計(jì)算量.如圖2所示，計(jì)算有效區(qū)域的網(wǎng)格數(shù)為400×400，在邊界外設(shè)定的隨機(jī)邊界的厚度是50層，整個(gè)計(jì)算區(qū)域的網(wǎng)格數(shù)為500×500，網(wǎng)格間距dx＝dz＝10m，采用的是空間四階交錯(cuò)網(wǎng)格和時(shí)間一階差分形式，時(shí)間步長(zhǎng)dt＝1ms，震源在區(qū)域中心，是標(biāo)準(zhǔn)的Ricker子波，峰值為1×10－3，正傳到最大時(shí)刻3s后逆?zhèn)鞯?.5s的波場(chǎng)和正傳波場(chǎng)在同一時(shí)刻處的絕對(duì)誤差的精度為10－11的量級(jí)，計(jì)算采用的是單精度浮點(diǎn)型，而單精度浮點(diǎn)型的有效位數(shù)是7位，因此，該量級(jí)的誤差僅僅是計(jì)算機(jī)精度的誤差.

隨機(jī)邊界和吸收邊界的設(shè)定也是值得討論的問(wèn)題.隨機(jī)邊界會(huì)帶來(lái)邊界漫反射影響，這在上表面的運(yùn)用中很明顯.當(dāng)上表面采用自由反射邊界條件的時(shí)候，又會(huì)帶來(lái)一些邊界反射假象，當(dāng)采用吸收邊界時(shí)，必須記錄波場(chǎng)上表面信息才能使得波場(chǎng)正確逆?zhèn)?因此，本文將在后面的模型和實(shí)際資料的運(yùn)用中對(duì)比分析設(shè)置四種不同邊界處理策略：第一種是四周均為隨機(jī)邊界；第二種是上表面為自由反射邊界，其它為隨機(jī)邊界；第三種是上表面是吸收邊界，其它為隨機(jī)邊界，記錄每時(shí)刻上表面處單層波場(chǎng)并在逆?zhèn)鬟^(guò)程中載入；第四種是四周均為吸收邊界，記錄每時(shí)刻四個(gè)邊界處單層波場(chǎng)并在逆?zhèn)鬟^(guò)程中載入.

圖2 隨機(jī)邊界速度模型和正傳逆?zhèn)髡`差（a）生成的隨機(jī)邊界速度模型，背景速度為2000m／s；（b）隨機(jī)邊界速度模型的橫向剖面；（c）使用隨機(jī)邊界正傳波場(chǎng)和逆?zhèn)鞑▓?chǎng)在0.5s的誤差；（d）利用記錄n／2＋1層邊界正傳波場(chǎng)和逆?zhèn)鞑▓?chǎng)在0.5s的誤差.Fig.2 The random velocity boundary model and the error of waveform between forward and backward propagation at the same time（a）The random velocity boundary model with background as 2000m／s；（b）One section of the random velocity boundary；（c）The error of wavefield between forward and backward propagation with random boundary at 0.5s；（d）The error of wavefield between forward and backward propagation with absorbing boundary by recording n／2＋1layers of boundary at 0.5s.

3 實(shí)例分析

本文將對(duì)二維崎嶇海底模型（圖3）和某海域?qū)嶋H資料（圖4）分別利用四種不同的邊界處理策略實(shí)現(xiàn)疊前逆時(shí)偏移成像處理，并分析對(duì)比其成像效果和計(jì)算效率，其中，吸收邊界采用的是20層的完全匹配層（PML）吸收邊界［24］，采用吸收邊界時(shí)記錄邊界處未衰減單層波場(chǎng)信息在逆?zhèn)鲿r(shí)載入.隨機(jī)邊界的方法是在邊界均采用50層的隨機(jī)速度層.其差分格式均采用的是空間四階交錯(cuò)網(wǎng)格和時(shí)間一階差分形式，成像條件均為零延遲互相關(guān)成像條件，在成像后均采用Laplace算子去除低頻成像噪音，震源均為標(biāo)準(zhǔn)的Ricker子波.本文測(cè)試的集群節(jié)點(diǎn)采用的是2個(gè)Quad－Core AMD Opteron（tm）Processor 2378 CPU，每個(gè)節(jié)點(diǎn)16G內(nèi)存，采用的是MPI并行環(huán)境.

3.1 模型資料對(duì)比

本文采用二維崎嶇海底模型的網(wǎng)格數(shù)為900×400，如圖3a所示，同一成像點(diǎn)的最大疊加次數(shù)為120，網(wǎng)格間距為dx＝dz＝12m，時(shí)間步長(zhǎng)dt＝1ms，時(shí)間上的總步數(shù)為6000.

如圖3b1和圖3b2中箭頭所示，當(dāng)四周均為隨機(jī)邊界時(shí)，上邊界的隨機(jī)速度介質(zhì)同人工震源距離太近，導(dǎo)致漫反射帶來(lái)的反射假象嚴(yán)重影響成像質(zhì)量.當(dāng)上表面采用自由反射邊界，其余采用隨機(jī)邊界時(shí)，成像如圖3c1所示，不需要額外存儲(chǔ)波場(chǎng)信息，計(jì)算時(shí)也不需要額外計(jì)算上邊界內(nèi)區(qū)域的波場(chǎng)傳播，內(nèi)存最少，計(jì)算效率最高，但會(huì)形成一些自由表面反射假象，在圖3c2中箭頭所示位置較為明顯.當(dāng)上表面采用吸收邊界其余三邊界采用隨機(jī)邊界時(shí)，需要記錄每時(shí)刻波場(chǎng)上邊界信息，成像結(jié)果如圖3d1所示，其細(xì)節(jié)圖3d2同圖3c2相比，沒(méi)有自由表面產(chǎn)生的邊界反射假象.如圖3e1所示，四邊界均采用吸收邊界，記錄波場(chǎng)邊界逆?zhèn)鲿r(shí)載入的成像結(jié)果是最好的.圖3e2的成像效果比圖3d2更加清晰干凈，如圖中箭頭所示，說(shuō)明隨機(jī)邊界帶來(lái)的漫反射影響不能通過(guò)疊加得到完全的抵消；總體能量相比圖3b1和圖3c1弱，這也證明了記錄一層邊界會(huì)帶來(lái)一定的能量損失，但不會(huì)影響成像質(zhì)量.當(dāng)模擬資料采用帶有多次波的記錄時(shí)，不管逆時(shí)偏移正傳逆?zhèn)鬟^(guò)程中上表面采用何種邊界時(shí)，成像均會(huì)產(chǎn)生多次波假象［25］，但若上表面采用自由反射邊界，會(huì)比采用吸收邊界多出一部分上表面反射假象，成像效果也是采用吸收邊界時(shí)更好.

崎嶇海底模型運(yùn)用四種不同邊界處理策略所需內(nèi)存，計(jì)算時(shí)間和成像效果定性分析如表2所示.四邊界均采用PML吸收邊界的策略用時(shí)最少，成像效果最好，但其內(nèi)存需求也最高.不過(guò)以現(xiàn)在計(jì)算機(jī)硬件條件，二維情況下的單層波場(chǎng)邊界的信息可以直接放入內(nèi)存中，不會(huì)帶來(lái)存儲(chǔ)問(wèn)題.

表2 崎嶇海底模型運(yùn)用不同邊界處理策略所需內(nèi)存，時(shí)間對(duì)比和成像效果定性分析Table 2 The contrast of memory and costing time for rugged water－bottom model by different boundary treatment and qualitative analysis of imaging results

3.2 實(shí)際資料對(duì)比

利用同崎嶇海底模型所采用的四種不同的邊界處理策略對(duì)某海域?qū)嶋H資料做疊前逆時(shí)偏移，其成像結(jié)果如圖4所示.

如圖4a1和圖4a2所示，當(dāng)四周均采用隨機(jī)邊界時(shí)，上表面隨機(jī)邊界的漫反射嚴(yán)重影響成像質(zhì)量.圖4b1中存在一些上邊界反射假象，背景速度簡(jiǎn)單的淺層比較明顯，圖4b2中虛線圓框內(nèi)所示，但由于已對(duì)該數(shù)據(jù)處理過(guò)表面相關(guān)多次波，所以上表面反射假象得到了一定的消除［25］，在深層時(shí)，通過(guò)疊加覆蓋假象掩藏在有效信號(hào)中，總體成像質(zhì)量影響不明顯.圖4c1和圖4c2中的隨機(jī)邊界漫反射的影響通過(guò)多次疊加得到有效壓制，幾乎不影響成像質(zhì)量，同時(shí)也證明了隨機(jī)邊界的可行性；同圖4b2相比，解決了上表面反射的影響.如圖4d1和圖4d2所示，當(dāng)四周均采用吸收邊界時(shí)成像效果最好.當(dāng)對(duì)實(shí)際資料處理采用上表面為自由表面的策略時(shí)，必須對(duì)表面相關(guān)多次波進(jìn)行處理，以減少假象的影響.

3.3 邊界處理策略分析

邊界層的厚度直接影響到總的計(jì)算量，經(jīng)數(shù)值試驗(yàn)分析，本文認(rèn)為給定隨機(jī)邊界的厚度為50層是比較合理的，吸收邊界采用的是20層的完全匹配層（PML）吸收邊界，與模型大小無(wú)關(guān).當(dāng)然，也可以選取更少的隨機(jī)邊界的厚度和吸收層的厚度，或者采用旁軸近似吸收等其它只需要更少層的吸收邊界條件，但其代價(jià)是邊界反射的增強(qiáng).

在計(jì)算二維區(qū)域時(shí)，假設(shè)計(jì)算區(qū)域網(wǎng)格數(shù)為1000×1000，時(shí)間上的步數(shù)為10000，采用單精度計(jì)算，四邊界均為吸收邊界，記錄單層波場(chǎng)邊界，則每炮計(jì)算所需記錄的邊界為：（4×1000）×10000×（4byte）≈152M，以現(xiàn)在的計(jì)算機(jī)硬件設(shè)施，完全可以全部放在內(nèi)存中存儲(chǔ)，不需要進(jìn)行額外的磁盤I／O操作.同時(shí)，三邊采用50層隨機(jī)邊界的計(jì)算區(qū)域也明顯比四邊采用20層吸收邊界的大，而計(jì)算區(qū)域的大小同計(jì)算量成正比.因此，本文認(rèn)為二維情況下疊前逆時(shí)偏移采用四邊界均為吸收邊界并記錄單層波場(chǎng)邊界逆?zhèn)鲿r(shí)載入的策略更具有成像效果和計(jì)算效率優(yōu)勢(shì).

當(dāng)計(jì)算三維區(qū)域時(shí)，若計(jì)算區(qū)域比較大時(shí)，假設(shè)計(jì)算區(qū)域網(wǎng)格數(shù)為1000×1000×500，時(shí)間步數(shù)為10000步.當(dāng)六個(gè)邊界均采用吸收邊界時(shí)，每炮需要記錄的邊界為：（2×1000×1000＋4×1000×500）×10000×（4byte）≈149.01G；當(dāng)上表面采用吸收邊界，其它邊界均為隨機(jī)邊界時(shí)，每炮需要記錄的邊界為：（1000×500）×10000×（4byte）≈18.6G.波場(chǎng)邊界的存儲(chǔ)將明顯影響計(jì)算效率，這時(shí)運(yùn)用上表面為自由反射邊界，其它五面為隨機(jī)邊界就更具有效率優(yōu)勢(shì).若計(jì)算區(qū)域較小時(shí)，也可適當(dāng)選取上表面為吸收邊界，或者六邊均為吸收邊界的方式.

4 結(jié) 論

疊前逆時(shí)偏移中由波場(chǎng)延拓方向不一致帶來(lái)的波場(chǎng)傳播歷史的存儲(chǔ)問(wèn)題，可以通過(guò)記錄波場(chǎng)邊界信息的吸收邊界和隨機(jī)邊界的方法得到有效解決.本文論證了吸收邊界只需記錄邊界處未衰減的單層波場(chǎng)信息即能保證逆?zhèn)鞑▓?chǎng)的波形完整性，只是能量上有少許缺失，不影響成像效果，極大地減少了邊界的存儲(chǔ)量.分析了四種不同的邊界處理策略，結(jié)果表明：四周均采用記錄單層波場(chǎng)邊界的吸收邊界的策略成像效果最好；上表面采用自由反射邊界，其它三邊采用隨機(jī)邊界的策略不需要額外存儲(chǔ)波場(chǎng)信息，經(jīng)過(guò)多次疊加之后隨機(jī)邊界漫反射對(duì)成像質(zhì)量影響很小，但會(huì)產(chǎn)生一定的上表面反射假象；當(dāng)上表面采用吸收邊界后，需記錄波場(chǎng)上表面信息，可消除上表面反射假象.在二維情況下，四邊界均采用記錄單層波場(chǎng)邊界的吸收邊界的策略更具有成像效果和計(jì)算效率優(yōu)勢(shì)；在三維情況下，特別是大區(qū)域逆時(shí)偏移的計(jì)算中，隨機(jī)邊界會(huì)更具有計(jì)算效率優(yōu)勢(shì).針對(duì)不同的計(jì)算問(wèn)題和計(jì)算機(jī)硬件條件，可以靈活使用吸收邊界和隨機(jī)邊界，兼顧計(jì)算效率和成像效果.

（

）

［1］ Baysal E，Kosloff D D，Sherwood J W C.Reverse time migration.Geophysics，1983，48（11）：1514－1524.

［2］ Whitmore N D.Iterative depth migration by backward time propagation.53rd Annual International Meeting.SEG，Expanded Abstracts，1983：382－385.

［3］ Claerbou J F.Toward a unified theory of reflector mapping.Geophysics，1971，36（3）：467－481.

［4］ Zhu J M，Lines L R.Comparison of Kirchhoff and reversetime migration methods with applications to prestack depth imaging of complex structures.Geophysics，1998，63（4）：1166－1176.

［5］ Biondi B，Shan G.Prestack imaging of overturned reflections by reverse time migration.72nd Annual International Meeting.SEG，Expanded Abstracts，2002.

［6］ Yoon K，Marfurt K J，Starr W.Challenges in reverse－time migration.74th Annual International Meeting.SEG，Expanded Abstracts，2004：1057－1060.

［7］ Fricke J R.Reverse－time migration in parallel：a tutorial.Geophysics，1988，53（9）：1143－1150.

［8］ Zhang J H，Wang S Q，Yao Z X.Accelerating 3DFourier migration with Graphics Processing Units.Geophysics，2009，74（6）：129－139.

［9］ 劉紅偉，李博，劉洪等.地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn).地球物理學(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.ChineseJ.Geophys.（in Chinese），2010，53（7）：1725－1733.

［10］ Baysal E，Kosloff D D，Sherwood J W C.A 2－way nonreflecting wave equation.Geophysics，1984，49（2）：132－141.

［11］ Loewenthal D，Mufti I R.Reversed time migration in spatial frequency domain.Geophysics，1983，48（5）：627－635.

［12］ Loewenthal D，Stoffa P L，F(xiàn)aria E L.Suppressing the unwanted reflections of the full wave equation.Geophysics，1987，52（7）：1007－1012.

［13］ Liu F Q，Zhang G Q，Morton S A，et al.Reverse－time migration using one－way wavefield imaging condition.77th Annual International Meeting.SEG，Expanded Abstracts，2007：2170－2174.

［14］ Liu F Q，Zhang G Q，Morton S A，et al.An effective imaging condition for reverse－time migration using wavefield decomposition.Geophysics，2011，76（1）：S29－S39.

［15］ Zhang Y，Sun J.Practical issues of reverse time migration：true amplitude gathers，noise removal and harmonic－source encoding.70th EAGE Conference ＆Exhibition，Rome，2008.

［16］ 楊仁虎，常旭，劉伊克.疊前逆時(shí)偏移影響因素分析.地球物理學(xué)報(bào)，2010，53（8）：1902－1913.

Yang R H，Chang X，Liu Y K.The influence factors analyses of imaging precision in pre－stack reverse time migration.ChinsesJ.Geophys.（in Chinese），2010，53（8）：1902－1913.

［17］ 楊仁虎.復(fù)雜介質(zhì)地震波傳播與逆時(shí)偏移成像方法研究［博士學(xué)位論文］.北京：中國(guó)科學(xué)院研究生院，2010.

Yang R H.The study on seismic wave propagation and reverse time migration in complex media［Ph.D.thesis］（in Chinese）.Beijing：Graduate School of Chinese Academy of Sciences，2010.

［18］ Symes W W.Reverse time migration with optimal checkpointing.Geophysics，2007，72（5）：SM213－SM221.

［19］ Dussaud E，Symes W W，Williamson P，et al.Computational strategies for reverse－time migration.78th Annual International Meeting.SEG，Expanded Abstracts，2008，27（1）：2267－2271.

［20］ Clapp R G.Reverse time migration with random boundaries.79th Annual International Meeting.SEG，Expanded Abstracts，2009：2809－2813.

［21］ Graves R W.Simulating seismic wave propagation in 3D elastic media using staggered－grid finite differences.Bulletin oftheSeismologicalSocietyofAmerica，1996，86（4）：1091－1106.

［22］ Liu Y，Sen M K.An implicit staggered－grid finite－difference method for seismic modelling.GeophysicalJournalInternational，2009，179（1）：459－474.

［23］ Kindelan M，Aamel A，Sguazzero P.On the construction and efficiency of staggered numerical differentiators for the wave equation.Geophysics，1990，55（1）：107－110.

［24］ Collino F，Tsogka C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media.Geophysics，2001，66（1）：294－307.

［25］ Liu Y K，Chang X，Jin D G，et al.Reverse time migration of multiples for subsalt imaging.Geophysics，2011，76（5）：WB209－WB216.