周輝，陳漢明，李緒宣，范廷恩

1中國石油大學油氣資源與探測國家重點實驗室，CNPC物探重點實驗室，北京 102249

2中國海洋石油集團公司研究總院，北京 100027

1 引言

偏移，特別是疊前偏移是目前地震資料處理的一個重要環(huán)節(jié).當前實現(xiàn)偏移成像的途徑大致分為兩類，即基于射線理論的Kirchhoff積分法偏移和基于波動方程的偏移.Kirchhoff積分法偏移是目前工業(yè)界廣泛使用的方法，它的基本理論幾何意義明確，這類偏移方法具有較靈活的目標處理能力和較高的計算效率，還可以適應不規(guī)則觀測系統(tǒng).但是，Kirchhoff積分類偏移方法的基礎是射線追蹤，在復雜介質(zhì)條件下，射線追蹤存在多路徑、焦散等問題，因為這些問題，積分類偏移方法有時得不到令人滿意的成像結(jié)果（Han，1998；金勝汶等，2002；劉喜武和劉洪，2002）.波動方程相對于射線追蹤技術，能更準確地描述波在復雜介質(zhì)中的傳播以及波場能量的變化情況，這使得基于波動方程的偏移方法更受青睞.波動方程偏移方法又可以進一步分為基于單程波方程的偏移方法和基于雙程波方程的偏移方法.精確的單程波方程只能在頻率-波數(shù)域表示，如相移法（Gazdag，1978）.但是相移法無法適應速度的橫向變化，因此一些學者在相移法的基礎上相繼提出了一些改進措施以適應速度的橫向變化（Gazdag and Sguazzero，1984；Stoffa etal.，1990；Ristow and Rühl，1994；Biondi，2002）.另外一種途徑是直接獲得在空間域表述的單程波方程，采用有限差分法進行數(shù)值計算，從而隱式地適應速度的任意變化.最早的時空域單程波方程由Claerbout（1976）提出，它是基于對雙程波方程頻散關系近似得到的.然而，在空間域表述的單程波方程有明顯的傾角限制，無法對高陡構(gòu)造準確成像.為了進一步擴大單程波的傳播角度以滿足大傾角地層成像的要求，一些學者采用了不同的高階近似方法，獲得了一系列高階單程波方程（馬在田，1983；Bamberger etal.，1988；Collins and Westwood，1991）.但這些方程往往包含高階偏導數(shù)項（二階以上），使得數(shù)值計算成為令人頭疼的問題.

基于雙程波的逆時偏移能有效地克服單程波成像傾角的限制，并且在保幅方面具有優(yōu)勢.逆時偏移（Baysal etal.，1983）包含波場的正向外推和逆時外推，其成像過程需要存儲各個時刻的波場，非常耗時.即便采取一些措施避免波場存儲的問題（Clapp，2009），但也是以計算換存儲的策略，會額外增加計算量.此外，適用于波動方程偏移的互相關成像條件最初是針對單程波方程提出的，將其直接用于逆時偏移，會產(chǎn)生低頻噪音干擾（Fletcher etal.，2006）.不可否認，逆時偏移是波動方程偏移方法的發(fā)展趨勢，但三維逆時偏移的計算量和對速度模型的高要求使得它在目前工業(yè)界的應用效果并不是十分理想.

本文的研究對象是Guddati（2006）提出的任意廣角波動方程（AWWE），它是一種在空間域表述的高精度單程波方程.AWWE的優(yōu)點是精度高，要提高AWWE的成像傾角，只需增加參考速度的個數(shù)即可.AWWE的形式簡單，只包含二階偏導數(shù)項，數(shù)值計算非常方便.因為這些優(yōu)點，AWWE已成功用于偏移成像（Guddati and Heidari，2005；何兵壽等，2008a；何兵壽等，2008b；孫歧峰和杜啟振，2011；Zhou etal.，2012）.Chen等（2013）將完美匹配層（PML）用于AWWE，有效地解決了AWWE數(shù)值計算中的吸收邊界問題.盡管AWWE具有高精度的特點，但是其計算過程包含重復的矩陣求逆和矩陣相乘，使得其計算代價比其他單程波方程高出許多，特別是增加參考速度個數(shù)時，其計算量顯著增加.

AWWE與其他單程波方程都存在的一個問題是倏逝波干擾，這是由單程波方程近似的本質(zhì)造成的.推導單程波方程的過程實際上是用一個有理函數(shù)近似平方根算子，當波數(shù)較大，平方根算子出現(xiàn)復數(shù)值時，近似的有理函數(shù)卻仍然是一個實數(shù)，這樣的誤差造成了倏逝波干擾.抑制倏逝波干擾有兩種途徑，一種是采用復有理函數(shù)近似平方根算子（Milinazzo etal.，1997），這種方法已在傅里葉有限差分偏移方法中得到應用，但該方法會因參數(shù)選擇不當產(chǎn)生不穩(wěn)定的問題（Amazonas etal.，2009）.另一種途徑是在倏逝波產(chǎn)生后進行壓制，如頻率-波數(shù)（f-k）濾波.Kosloff和Baysal（1983）在進行波場深度延拓時，逐層采用f-k濾波抑制倏逝波干擾.用f-k濾波方法抑制倏逝波的關鍵在于合理選取濾波所需要的視速度.Kosloff和Baysal（1983）采用延拓所在層的最大速度作為濾波器的門檻，會損失一部分高角度波場信息，從而降低對高陡構(gòu)造的成像能力，這已被Sandberg和Beylkin（2009）證實.此外，逐層進行濾波也會增加額外的計算量.

本文的研究內(nèi)容針對上述兩個問題，即AWWE計算量大和存在明顯的倏逝波干擾.首先，本文改進現(xiàn)有的有限差分計算方案，包括內(nèi)部區(qū)域和PML吸收邊界區(qū)域，使其計算效率大幅度提高.其次，采用一種有效的方法，能在f-k域比較準確地區(qū)分倏逝波和非倏逝波，從而獲得一個合理的視速度作為f-k濾波器的門檻，利用f-k濾波方法壓制倏逝波干擾，提高成像質(zhì)量.

2 任意廣角波動方程深度延拓成像的數(shù)值實現(xiàn)

2.1 AWWE疊前偏移成像

采用AWWE進行疊前深度延拓成像時，分別采用下行波方程和上行波方程將震源波場和記錄波場沿深度方向延拓，在同一個深度上將兩個波場互相關得到該深度上的像.因此，震源波場和記錄波場分別滿足定解問題（1）和（2），

其中uS（x，t）和uR（x，t）分別為震源波場和記錄波場.上述方程中包含的其他變量定義如下：

c為地震波在介質(zhì)中的傳播速度，ci（i＝1，2，…，n）為參考速度，d＝ （1，0，…，0）1×n，u＝ （u，u1，u2，…，un-1）T，標量u為位移變量，ui（i＝1，2，…，n-1）為輔助變量.參考速度與傳播速度的比例影響AWWE的精度，事實上cosθi＝c／ci，θi為AWWE傳播算子與精確的單程波算子相匹配的角度，本文采用ci＝λic（其中λi為常數(shù)）這種方法設置參考速度，而不是像Guddati和Heidari（2005）部分偏移算例中，取參考速度為常數(shù)，而參考速度與傳播速度的比值隨空間變化.后續(xù)數(shù)值實驗將證實ci＝λic（其中λi為常數(shù)）這種選取參考速度的方法的合理性.

此外，參考速度的個數(shù)也直接影響AWWE的精度.增加參考速度的個數(shù)，同時適當選取參考速度與傳播速度的比值，可以顯著地提高AWWE的精度，從而增大成像傾角.為方便敘述，文中簡記帶n個參考速度的AWWE為AWWEn.

圖1描述了精確的單程波算子，AWWE2、AWWE3以及AWWE5歸一化的慢度曲線，圖1a為實部，圖1b為虛部，考慮到近似單程波算子的虛部都為零，圖中用AWWEn概括之.從圖1a可以看出，參考速度為（c，4c）的AWWE2已經(jīng)具有很高的精度，其成像傾角大約為80°.圖1b表明，當水平慢度大于1時，精確的單程波算子的垂直慢度為純虛數(shù)，而近似的單程波算子的垂直慢度卻仍為實數(shù)，這種近似誤差造成了單程波數(shù)值計算中的倏逝波干擾.

2.2 有限差分離散方案

考慮到方程（1）和（2）的離散方案類似，這里只以上行波方程為例進行說明.定義差分算子如下：

圖1 帶不同參考速度的AWWE精度分析，（a）為實部，（b）為虛部Fig.1 Accuracy of AWWE with various numbers of reference velocity including real part（a）and imaginary part（b）

其中 （i，j，k）分別為x，z和t方向離散網(wǎng)格的節(jié)點號，（Δx，Δz，Δt）分別為空間和時間的網(wǎng)格間距.采用上述差分算子離散方程（2）得到

利用Crank-Nicolson方法沿深度外推波場，即滿足

將方程（7）代入方程（6）進一步整理得到上行波方程沿深度外推的逆時計算格式為

其中，

該差分格式由Guddati和Heidari（2005）給出，精度為O（Δt2+Δz2+Δx2），其穩(wěn)定性條件為cmaxΔt≤Δx，cmax為波在介質(zhì)中傳播的最大速度.下行波方程（1）對應的深度延拓時間正向外推差分格式也可以類似得到.利用方程（9）進行計算時，需要計算逆矩陣，考慮到αz是隨空間變化的，因此計算過程中需要反復求逆，并需要反復計算矩陣H2和H3，這極大地降低了AWWE的計算效率，并且隨著參考速度個數(shù)的增加，其計算量顯著地增加.

本文從方程（8）出發(fā)，修改計算格式，將求逆的次數(shù)減少為一次，從而能顯著地提高計算效率.注意到方程（3）和（4），在網(wǎng)格剖分確定的情況下，如果參考速度ci＝λic，其中λi為恒定的常數(shù)，則矩陣Λ1和Λ2為恒定的矩陣.事實上λi為恒定的常數(shù)這一假設是合理的，后述數(shù)值實驗將予以證實.在上述假設的前提下，將方程（8）修改為

該方程可進一步寫為

利用矩陣D（只有一個非零元素）和向量d的特點，方程（12）可以簡化為

其中，

方程（13）和（14）中的矩陣A和H2不隨空間變化，因此只需要計算一次.表1列出了方程（13）和方程（9）的計算量，單位是15°單程波方程計算量的倍數(shù).從表中可以看出，采用改進后的計算方案AWWE的計算效率顯著提高，AWWE2的計算量僅僅為15°單程波方程的1.67倍，AWWE5的計算量甚至比采用原方案的AWWE3的計算量還要小.

2.3 AWWE的PML匹配方程及其差分格式

Chen等（2013）推導了AWWE的PML匹配方程并給出了相應的有限差分計算格式，解決了AWWE正演模擬（偏移）中截斷邊界的問題.本文不重述PML的推導過程，而是直接給出分裂式的PML匹配方程.與內(nèi)部區(qū)域類似，PML區(qū)域的差分計算公式也可作相應改進，以提高PML區(qū)域的計算效率.下行波方程（1）對應的PML匹配方程為

其中，

dx（x）為PML區(qū)域的衰減函數(shù)，.直接給出方程（16）改進之后的差分計算公式為

其中，

方程（17）和（18）為下行波方程深度延拓時間正向外推時PML區(qū)域的更新公式，該公式中也避免了矩陣的反復求逆，比Chen等（2013）之前給出的差分公式的計算效率高.上行波方程對應的PML匹配方程及其差分計算公式也可作類似推導.大量的數(shù)值實驗證實內(nèi)部區(qū)域和PML區(qū)域差分格式的穩(wěn)定性條件為，二維情形cmaxΔt≤Δx，三維情形cmaxΔt≤min（Δx，Δy）.

3 傾角濾波抑制倏逝波干擾

倏逝波是單程波方程數(shù)值計算中都存在的問題，在引言中已經(jīng)說明了其產(chǎn)生的原因，并用圖1b進行了解釋.抑制倏逝波的常用方法是傾角濾波，即f-k濾波，在設計f-k濾波器時需要一個視速度作為門檻.基于觀察和實驗發(fā)現(xiàn)，當AWWE的參考速度有微小擾動時，非倏逝波場幾乎不發(fā)生變化.將微小擾動參考速度前的波場與微小擾動之后的波場相減，殘留波場主要反映的是倏逝波的能量.將殘差波場（后稱殘差記錄）變換到f-k域得到其振幅譜，該振幅譜能清晰地反映倏逝波的能量分布，因此能提供一個合理的視速度作為f-k濾波器的門檻.

表1 差分格式改進前后的計算量比較（二維情形）Table 1 Computational efficiency comparison between the improved scheme and the existing scheme in 2Dcase

基于以上考慮，在AWWE疊前偏移的過程中采用如下方法壓制倏逝波的干擾：先將波場延拓一次得到記錄u1（x，t），然后微小擾動參考速度（選取的擾動量為10-4），重新做一次延拓，得到記錄u2（x，t），這兩個記錄的差記為 Δu（x，t），稱為殘差記錄.根據(jù)殘差記錄Δu（x，t）的振幅譜獲得一個合理的視速度，然后設計光滑的扇形濾波器對記錄u1（x，t）濾波以壓制倏逝波成分，再用濾波后的記錄進行后續(xù)的深度延拓，其流程如圖2所示.

嚴格意義上，在深度延拓的過程中需要逐層進行濾波，但逐層濾波會額外增加計算量.出于以上考慮，本文只在某些深度上加入上述濾波處理，盡管減少了濾波次數(shù)，但抑制倏逝波干擾的效果仍然很明顯，后續(xù)數(shù)值實驗將予以證實.

圖2 加入f-k濾波處理后的偏移流程Fig.2 The migration flow with an f-k filtering procedure

4 數(shù)值實驗

4.1 AWWE正演模擬

圖3為二維SEG／EAGE鹽丘速度模型，網(wǎng)格大小 （Δx，Δz）＝ （5m，4m），由于高速鹽丘體的存在，橫向上速度存在劇烈變化.圖4a和4b為AWWE2模擬鹽丘模型中的下行波場的波場切片，參考速度（c1，c2）＝ （c，4c），震源靠近左邊界，左邊界采用PML吸收邊界（Chen etal.，2013）.圖4a和4b的時刻相同，分別采用改進前的差分計算公式（9）和改進后的計算公式（13）進行計算，兩圖在相同的色標范圍內(nèi)顯示.圖4說明改進之后的計算方案仍能穩(wěn)定、精確地模擬下行波場，與改進前的計算結(jié)果完全一致.此外，在該算例中，采用改進后的計算方案，使用相同層數(shù)的PML，PML區(qū)域的計算效率也提高了7.6%.

圖5為AWWE3模擬均勻介質(zhì)中的下行波場切片，以x＝1500m為界限（圖中黑線），左側(cè)區(qū)域選取的三個參考速度依次為傳播速度的1倍，2倍和4倍，而右側(cè)區(qū)域三個參考速度依次為傳播速度的0.8倍，1.5倍和2.4倍，因此參考速度與傳播速度的比例隨空間是變化的.很明顯，在分界處產(chǎn)生了虛假反射，該反射并非由波在介質(zhì)中的傳播速度差異造成，而是由于左右區(qū)域參考速度與傳播速度的比值不一致造成的.同理，當波的傳播速度劇烈變化，而參考速度取固定值時，也會產(chǎn)生虛假反射.因此，前文假設參考速度是傳播速度的常數(shù)倍是合理的.

4.2 AWWE偏移試算

4.2.1 SEG／EAGE鹽丘速度模型

為了檢驗AWWE的成像能力以及f-k濾波方法壓制倏逝波的效果，采用二維SEG／EAGE鹽丘模型（圖3）進行驗證.離散網(wǎng)格大小為 （Δx，Δz）＝（5m，4m），利用該模型合成12炮地震數(shù)據(jù)，炮間距250m，每炮501道接收，道間距5m，記錄時間2.5s，震源采用主頻為25Hz的雷克子波.利用參考速度為（c，4c）的AWWE2的偏移算子進行深度延拓偏移成像，其成像結(jié)果為圖6a，該成像結(jié)果有很強的倏逝波干擾.為了抑制倏逝波干擾，在深度延拓開始之前分別對地表的震源波場和記錄波場進行一次f-k濾波，用地表所在網(wǎng)格層的最大速度作為濾波器的視速度，其成像結(jié)果為圖6b，對比圖6a，倏逝波得以抑制，成像質(zhì)量明顯改善.為了進一步抑制倏逝波的干擾，在深度分別為4m、8m和12m的網(wǎng)格層上對震源和記錄波場各進行了一次f-k濾波，濾波器的視速度利用殘差記錄的振幅譜獲得，其流程如圖2所示，最終的成像結(jié)果如圖6c所示.對比圖6b，在200～400m的深度范圍內(nèi)倏逝波的干擾進一步被壓制，如圖6b中白色箭頭所指.在圖6c的基礎上，在440m深度上對震源和記錄波場再次進行f-k濾波，濾波器的視速度利用殘差記錄的振幅譜獲得，其成像結(jié)果為圖6d.與圖6c相比，鹽丘正上方的小斷塊內(nèi)（圖6d中白色箭頭所指）的倏逝波干擾更弱.圖6e與圖6d不同，在440m深度上進行f-k濾波時采用橫向上的最大速度作為濾波器的視速度，由于視速度選取不合理，導致440m以下深度不能準確成像.

圖3 二維SEG／EAGE鹽丘速度模型，網(wǎng)格大小 （Δx，Δz）＝ （5m，4m）Fig.3 2DSEG／EAGE salt velocity model with a grid size of（5m，4m）

圖4 AWWE2模擬二維SEG／EAGE鹽丘模型中的下行波場切片（a）采用改進前的計算方案，（b）采用改進后的計算方案.Fig.4 Snapshots propagated by downward propagating AWWE2in 2DSEG／EAGE salt model（a）and（b）are respectively calculated by using the original calculation scheme and the improved scheme.

圖5 不合理的參考速度設置引起的虛假反射Fig.5 Spurious artifacts caused by unreasonable setting of reference velocities

圖6 SEG／EAGE鹽丘速度模型的偏移結(jié)果（a）無濾波；（b）在地表進行f-k濾波；（c）在4m、8m和12m的深度上進行f-k濾波；（d）在4m、8m、12m和440m的深度上進行f-k濾波；（e）在4m、8m、12m和440m的深度上進行f-k濾波，第四次濾波的視速度取該深度上速度的最大值.Fig.6 Migration results of the SEG／EAGE salt model（a）Without f-kfiltering；（b）f-kfiltering at surface；（c）f-kfiltering at depth of 4m，8m，and 12m；（d）f-kfiltering at depth of 4m，8m，12m，and 440m；（e）f-kfiltering at depth of 4m，8m，12m，and 440m，with the lateral maximum velocity being the apparent velocity in the last f-kfiltering.

圖7 是利用殘差記錄的振幅譜獲取440m深度上f-k濾波器視速度門檻值的示意圖.圖7a是將地表的雷克子波震源延拓至深度440m得到的記錄，利用參考速度為（c，4c）的下行AWWE2延拓算子；圖7b是將地表的雷克子波震源延拓至深度440m得到的記錄，利用參考速度為 （1.0001c，4.0001c）的下行AWWE2延拓算子；圖7c是圖7a與圖7b的差，即殘差記錄；圖7d是殘差記錄的振幅譜，其中虛線為拾取的視速度1300m·s-1，實線為所在深度橫向上的最大速度4185m·s-1，用最大速度作為濾波器的門檻會濾掉有效的高傾角成分，從而損害濾波層以下地層的成像，如圖6e所示.

4.2.2 Marmousi模型

為了進一步驗證AWWE對陡傾角構(gòu)造的成像能力以及f-k濾波方法的有效性，用Marmousi模型（圖8）進行偏移試驗.速度模型的網(wǎng)格大小為（Δx，Δz）＝ （10m，5m），使用參考速度為 （c，4c）的AWWE2算子進行偏移成像，數(shù)值計算采用文中改進之后的差分計算方案.51炮合成地震數(shù)據(jù)的偏移結(jié)果（沒有濾波處理）為圖9a，僅在前四個網(wǎng)格深度上加入f-k濾波處理之后的成像結(jié)果為圖9b，兩圖在相同的色標范圍內(nèi)顯示，兩圖的成像結(jié)果都證實了AWWE偏移算子對高陡構(gòu)造的成像能力.在深度延拓過程中只進行了數(shù)次f-k濾波，其增加的計算量基本可以忽略，但其壓制倏逝波干擾進而提高成像質(zhì)量的效果十分明顯.

圖7 利用殘差記錄的振幅譜獲取440m深度上f-k濾波器的視速度（a）參考速度為（c，4c）的AWWE2算子延拓后的記錄；（b）參考速度為（1.0001c，4.0001c）的AWWE2算子延拓后的記錄；（c）殘差記錄；（d）殘差記錄的振幅譜，虛線代表視速度為1300m·s-1，實線代表視速度為4185m·s-1.Fig.7 Obtaining apparent velocity threshold at depth of 440musing the amplitude spectrum of residual wavefield（a）The wavefield downward-propagated by AWWE2with reference velocities（c，4c）；（b）The wavefield downward-propagated by AWWE2with reference velocities（1.0001c，4.0001c）；（c）Residual source wavefield；（d）f-k amplitude spectrum of（c），the dashed line represents an apparent velocity 1300m·s-1，and the solid line 4185m·s-1.

圖8 Marmousi速度模型Fig.8 Marmousi velocity model

圖9 Marmousi模型偏移結(jié)果，（a）無濾波處理，（b）進行4次f-k濾波處理Fig.9 The migration results of Marmousi model，（a）without and（b）with four times f-k filtering

5 結(jié)論

AWWE是一種高精度的單程波方程，其成像傾角接近90°，可以對高陡構(gòu)造準確成像，但是其計算量較大，特別是增加參考速度的個數(shù)時，計算量會顯著增加.本文通過分析和數(shù)值實驗證實，參考速度必須為背景速度的常數(shù)倍以避免虛假反射，基于此分析，本文改進了現(xiàn)有的有限差分計算方案，能避免重復的矩陣求逆并減少矩陣相乘運算的次數(shù)，從而顯著地提高了AWWE的計算效率.在分析倏逝波產(chǎn)生原因的基礎上，采用了f-k濾波方法來抑制倏逝波對成像結(jié)果影響，在設計扇形濾波器時，根據(jù)殘差記錄的振幅譜來獲取視速度的門檻值.數(shù)值實驗說明，即使速度橫向變化，殘差記錄的振幅譜也能很好地描述倏逝波的能量分布，從而提供合理的視速度門檻信息.合成數(shù)據(jù)的偏移試算證實，進行少數(shù)幾次f-k濾波處理就能很好地抑制倏逝波干擾，提高成像質(zhì)量.

Amaz onas D，Costa J C，Schleicher J，etal.2009.Wide-angle FD and FFD migration using complex padéapproximations.Geophysics，72（6）：S215-S220.

Bamberger A，Engquist B，Halpern L，etal.1988.High-order paraxial wave-equation approximations in heterogeneous media.SIAM J.Appl.Math.，48（1）：129-154.

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

Biondi B.2002.Stable wide-angle Fourier finite-difference downward extrapolation of 3-D wavefields.Geophysics，67（3）：872-882.

Chen H M，Zhou H，Lin H，etal.2013.Application of perfectly matched layer for scalar arbitrarily wide-angle wave equations.Geophysics，78（1）：T29-T39.

Claerbout J F.1976.Fundamentals of Geophysical Data Processing.Oxford：Scientific Publications.

Clapp R G.2009.Reverse time migration with random boundaries.79th Ann.Inernat Mtg.，Soc.Expi.Geophys.Expanded Abstracts，2809-2813.

Collins M D，Westwood E K.1991.A higher-order energyconserving parabolic equation for range-dependent ocean depth，sound speed，and density.J.Acoust.Soc.Am.，89（3）：1068-1175.

Fletcher R P，F(xiàn)owler P J，Kitchenside P，etal.2006.Suppressing unwanted internal reflections in prestack reverse-time migration.Geophysics，71（6）：E79-E82.

Gazdag J，Sguazzero P.1984.Migration of seismic data by phaseshift plus interpolation.Geophysics，49（2）：124-131.

Gazdag J.1978.Wave equation migration with the phase-shift method.Geophysics，43（7）：1342-1351.

Guddati M N，Heidari A H.2005.Migration with arbitrarily wideangle wave equations.Geophysics，70（3）：S61-S70.

Guddati M N.2006.Arbitrarily wide-angle wave equations for complex media.Comput.Methods Appl.Mech.Eng.，195（1-3）：65-93.

Han B A.1998.Comparison of four depth migration methods.66th Ann.Internat Mtg.，Soc.Expi.Geophys.Expanded Abstracts，1104-1107.

He B S，Zhang H X，Han L H.2008a.Prestack reverse-time depth migration of two acoustic wave equation and comparison between both.Oil Geophys.Prosp.（in Chinese），43（6）：661-684.

He B S，Zhang H X，Zhang J.2008b.Prestack reverse-time depth migration of arbitrarily wide-angle wave equations.Acta Seismologica Sinica（in Chinese），30（5）：491-499.

Jin S W，Xu S R，Wu R S.2002.Wave equation based prestack depth migration using generalized screen propagator.Chinese J.Geophys.（in Chinese），45（5）：684-690.

Kosloff D D，Baysal E.1983.Migration with the full acoustic wave equation.Geophysics，48（6）：677-687.

Liu X W，Liu H.2002.Status and progress on wave equation migration methods.Progress in Geophysics（in Chinese），17（4）：582-591.

Ma Z T.1983.A splitting-up method for solution of higher-order migration equation by finite-difference scheme.Chinese J.Geophys.（Acta Geophysica Sinica）（in Chinese），26（4）：377-389.

Milinazzo F A，Zala C A，Brooke G H.1997.Rational square-root approximations for parabolic equation algorithms.J.Acoust.Soc.Am.，101（2）：760-766.

Ristow D，Rühl T.1994.Fourier finite-difference migration.Geophysics，59（12）：1882-1893.

Sandberg K，Beylkin G.2009.Full-wave equation depth extrapolation for migration.Geophysics，74（6）：WCA121-WCA128.

Stoffa P L，F(xiàn)okkema J T，de Luna Freire R M，etal.1990.Splitstep Fourier migration.Geophysics，55（4）：410-421.

Sun Q F，Du Q Z.2011.The frequency-space domain prestack depth migration with arbitrarily wide-angle wave equation.Oil Geophys.Prosp.（in Chinese），46（6）：890-896.

Zhou H，Lin H，Sheng S B，etal.2012.High angle prestack depth migration with absorption compensation.Applied Geophysics，9（3）：293-300.

附中文參考文獻

何兵壽，張會星，韓令賀等.2008a.兩種聲波方程的疊前逆時深度偏移及比較.石油地球物理勘探，43（6）：661-684.

何兵壽，張會星，張晶.2008b.任意廣角波動方程疊前逆時深度偏移.地震學報，30（5）：491-499.

金勝汶，許士勇，吳如山.2002.基于波動方程的廣義屏疊前深度偏移.地球物理學報，45（5）：684-690.

劉喜武，劉洪.2002.波動方程地震偏移成像方法的現(xiàn)狀與進展.地球物理學進展，17（4）：582-591.

馬在田.1983.高階方程偏移的分裂算法.地球物理學報，26（4）：377-389.

孫歧峰，杜啟振.2011.任意廣角波動方程頻率-空間域疊前深度偏移成像.石油地球物理勘探，46（6）：890-896.