提高任意廣角波動(dòng)方程計(jì)算效率及壓制倏逝波干擾方法

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

1中國石油大學(xué)油氣資源與探測(cè)國家重點(diǎn)實(shí)驗(yàn)室，CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249

2中國海洋石油集團(tuán)公司研究總院，北京 100027

1 引言

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

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

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

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

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

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

2.1 AWWE疊前偏移成像

采用AWWE進(jìn)行疊前深度延拓成像時(shí)，分別采用下行波方程和上行波方程將震源波場和記錄波場沿深度方向延拓，在同一個(gè)深度上將兩個(gè)波場互相關(guān)得到該深度上的像.因此，震源波場和記錄波場分別滿足定解問題（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，標(biāo)量u為位移變量，ui（i＝1，2，…，n-1）為輔助變量.參考速度與傳播速度的比例影響AWWE的精度，事實(shí)上cosθi＝c／ci，θi為AWWE傳播算子與精確的單程波算子相匹配的角度，本文采用ci＝λic（其中λi為常數(shù)）這種方法設(shè)置參考速度，而不是像Guddati和Heidari（2005）部分偏移算例中，取參考速度為常數(shù)，而參考速度與傳播速度的比值隨空間變化.后續(xù)數(shù)值實(shí)驗(yàn)將證實(shí)ci＝λic（其中λi為常數(shù)）這種選取參考速度的方法的合理性.

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

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

2.2 有限差分離散方案

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

圖1 帶不同參考速度的AWWE精度分析，（a）為實(shí)部，（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é)點(diǎn)號(hào)，（Δx，Δz，Δt）分別為空間和時(shí)間的網(wǎng)格間距.采用上述差分算子離散方程（2）得到

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

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

其中，

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

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

該方程可進(jìn)一步寫為

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

其中，

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

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

Chen等（2013）推導(dǎo)了AWWE的PML匹配方程并給出了相應(yīng)的有限差分計(jì)算格式，解決了AWWE正演模擬（偏移）中截?cái)噙吔绲膯栴}.本文不重述PML的推導(dǎo)過程，而是直接給出分裂式的PML匹配方程.與內(nèi)部區(qū)域類似，PML區(qū)域的差分計(jì)算公式也可作相應(yīng)改進(jìn)，以提高PML區(qū)域的計(jì)算效率.下行波方程（1）對(duì)應(yīng)的PML匹配方程為

其中，

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

其中，

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

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

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

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

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

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

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

4 數(shù)值實(shí)驗(yàn)

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的時(shí)刻相同，分別采用改進(jìn)前的差分計(jì)算公式（9）和改進(jìn)后的計(jì)算公式（13）進(jìn)行計(jì)算，兩圖在相同的色標(biāo)范圍內(nèi)顯示.圖4說明改進(jìn)之后的計(jì)算方案仍能穩(wěn)定、精確地模擬下行波場，與改進(jìn)前的計(jì)算結(jié)果完全一致.此外，在該算例中，采用改進(jìn)后的計(jì)算方案，使用相同層數(shù)的PML，PML區(qū)域的計(jì)算效率也提高了7.6%.

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

4.2 AWWE偏移試算

4.2.1 SEG／EAGE鹽丘速度模型

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

圖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）采用改進(jìn)前的計(jì)算方案，（b）采用改進(jìn)后的計(jì)算方案.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 不合理的參考速度設(shè)置引起的虛假反射Fig.5 Spurious artifacts caused by unreasonable setting of reference velocities

圖6 SEG／EAGE鹽丘速度模型的偏移結(jié)果（a）無濾波；（b）在地表進(jìn)行f-k濾波；（c）在4m、8m和12m的深度上進(jìn)行f-k濾波；（d）在4m、8m、12m和440m的深度上進(jìn)行f-k濾波；（e）在4m、8m、12m和440m的深度上進(jìn)行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，實(shí)線為所在深度橫向上的最大速度4185m·s-1，用最大速度作為濾波器的門檻會(huì)濾掉有效的高傾角成分，從而損害濾波層以下地層的成像，如圖6e所示.

4.2.2 Marmousi模型

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

圖7 利用殘差記錄的振幅譜獲取440m深度上f-k濾波器的視速度（a）參考速度為（c，4c）的AWWE2算子延拓后的記錄；（b）參考速度為（1.0001c，4.0001c）的AWWE2算子延拓后的記錄；（c）殘差記錄；（d）殘差記錄的振幅譜，虛線代表視速度為1300m·s-1，實(shí)線代表視速度為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）進(jìn)行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°，可以對(duì)高陡構(gòu)造準(zhǔn)確成像，但是其計(jì)算量較大，特別是增加參考速度的個(gè)數(shù)時(shí)，計(jì)算量會(huì)顯著增加.本文通過分析和數(shù)值實(shí)驗(yàn)證實(shí)，參考速度必須為背景速度的常數(shù)倍以避免虛假反射，基于此分析，本文改進(jìn)了現(xiàn)有的有限差分計(jì)算方案，能避免重復(fù)的矩陣求逆并減少矩陣相乘運(yùn)算的次數(shù)，從而顯著地提高了AWWE的計(jì)算效率.在分析倏逝波產(chǎn)生原因的基礎(chǔ)上，采用了f-k濾波方法來抑制倏逝波對(duì)成像結(jié)果影響，在設(shè)計(jì)扇形濾波器時(shí)，根據(jù)殘差記錄的振幅譜來獲取視速度的門檻值.數(shù)值實(shí)驗(yàn)說明，即使速度橫向變化，殘差記錄的振幅譜也能很好地描述倏逝波的能量分布，從而提供合理的視速度門檻信息.合成數(shù)據(jù)的偏移試算證實(shí)，進(jìn)行少數(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.

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

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

何兵壽，張會(huì)星，張晶.2008b.任意廣角波動(dòng)方程疊前逆時(shí)深度偏移.地震學(xué)報(bào)，30（5）：491-499.

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

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

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

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