波動方程疊前深度偏移直接產(chǎn)生角道集

劉禮農(nóng)，張 輝，張劍鋒

中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院油氣資源研究重點實驗室，北京 100029

1 引 言

隨著復(fù)雜構(gòu)造成像及巖性勘探的深入，地震資料疊前偏移成像處理中由不同炮或不同偏移距的偏移結(jié)果得到地下反射點對不同角度入射波場的反射強度，進而利用AVO或AVA響應(yīng)直接識別流體和含油氣情況成為一個日益受到關(guān)注的環(huán)節(jié).

波動方程疊前深度偏移是用較準確的波動理論來描述地震波的傳播過程，能正確模擬復(fù)雜上覆地層導(dǎo)致的多值走時等復(fù)雜的波傳播現(xiàn)象和地震波傳播的幾何擴散現(xiàn)象，因此既可對復(fù)雜構(gòu)造正確成像，又有望得到正確反映地下構(gòu)造反射強弱的成像幅值.但服務(wù)于復(fù)雜構(gòu)造成像的波動方程疊前深度偏移多是通過分炮偏移實現(xiàn)的，且每炮偏移實現(xiàn)時需要對原始炮集做鑲邊處理以滿足偏移孔徑的要求；這種數(shù)據(jù)集組織方式易于得到疊加剖面，但由這些炮集偏移結(jié)果抽取隨偏移距變化的共反射點（CRP）道集將產(chǎn)生巨大的硬盤存儲需求.此外，既然真正反映地下構(gòu)造物性參數(shù)對比的是隨入射角變化的角度域成像道集，人們就希望能由分炮偏移結(jié)果直接得到角道集.

已發(fā)展了一系列通過炮域波動方程疊前深度偏移產(chǎn)出角度域成像道集的方法.這些方法可歸納為如下兩類：一是在應(yīng)用成像條件前對入射和反傳波場進行局部平面波分解，對每對平面波分別成像［1－3］；二是在炮域偏移的成像條件中應(yīng)用空間移動或時間移動的成像條件，進而通過傾斜疊加等角度變換方法得到角度道集［4－6］.第一種方法需對逐個成像點進行局部平面波分解，將單點的成像轉(zhuǎn)換成一個成像矩陣，三維問題的計算量難以承受；第二種方法盡管計算量較第一種有所減少，但空間移動或時間移動的成像條件導(dǎo)致成像空間巨大，計算和內(nèi)存需求隨之增加，三維問題的成像空間將變?yōu)槲寰S.兩種方法均不能考慮炮點覆蓋不均勻?qū)е碌某上穹嫡`差，且波場延拓的空間采樣限制了角道集的有效頻率范圍.

本文在文獻［7］所發(fā)展的單程波算子地震波入射角計算方法的基礎(chǔ)上，通過計算震源入射波場與地層界面夾角，發(fā)展了由炮域偏移的結(jié)果直接通過部分疊加得到角道集的方法；該方法通過應(yīng)用“保幅”的反褶積成像條件和考慮角度區(qū)間內(nèi)累加的炮數(shù)，解決了炮點覆蓋不均勻?qū)е碌某上穹嫡`差問題.與前述兩類方法相比，這一方法所需的計算量及內(nèi)存幾乎可以忽略，角道集的有效頻率范圍也與偏移剖面相同，因而更適于現(xiàn)行的計算能力和工業(yè)資料的實際情況.

本文首先簡述基于單程波算子的地震波入射波之波前面法線、地層界面傾角估計以及以此為基礎(chǔ)的地震波入射角計算方法，然后建議了一種考慮炮點覆蓋不均勻的基于反褶積成像條件的直接累加生成角度域成像道集的方法，最后給出了實際應(yīng)用的計算流程和基于角道集改善地震成像效果的實用方法.文中算例表明，這一炮域波動方程深度偏移直接產(chǎn)生角道集方法既能剔除復(fù)雜上覆地層的影響，又可補償觀測系統(tǒng)非均勻覆蓋對成像幅值的影響；而與常規(guī)炮域波動方程疊前深度偏移相比，并不增加過多內(nèi)存需求與計算量.

2 反射角計算

2.1 基于單程波算子的入射波之波前面法線計算

單程波算子是通過對標準的聲波方程進行分解而得到的，一般情況人們總是沿垂直方向分解單程波算子.單程波算子可正確模擬復(fù)雜介質(zhì)中地震波的單向傳播，與標準的聲波（雙程波）方程相比，單程波算子沿其分解方向不發(fā)生地震波反射.由于這一特征，單程波算子所得到的簡諧波場就對應(yīng)著地震波傳播過程中的波前面.地震波場的波數(shù)是與波前面的法線方向相聯(lián)系的，而波場的空間導(dǎo)數(shù)又與波數(shù)存在著對應(yīng)關(guān)系，因此，我們可利用求取空間導(dǎo)數(shù)前后簡諧波場的變化，簡單求得全場的地震波入射角方向.

由單程波方程可知，頻率－波數(shù)域下行的波場P＋（x，y，z；ω）滿足

其中j是單位虛數(shù)，kz是垂直波數(shù)，ω是角頻率，x和y是水平位置坐標，z是深度坐標，由頻散關(guān)系得

式中k＝ω／c（x，y，z）為波數(shù)，c（x，y，z）是空間各位置處的介質(zhì)波速，角度θx，θy，θz分別是空間位置（x，y，z）處波前面法線與水平坐標軸x，y以及深度軸z的夾角.以z方向為例，將（2a）式代入（1）式并做離散的時域傅氏反變換，且令時間t＝0，可得

式中符號Re代表求復(fù)波場值的實部.若令～P＋（x，y，z＝0，ω）＝δ（x－xs，y－ys）f（ω），其中f（ω）是地震子波的譜，則（3）式中Σω＞0Re（～P＋（x，y，z；ω））對應(yīng)于震源點置于 （xs，ys，0）處時的脈沖響應(yīng)，而波前面法線就是 （xs，ys，0）處炮點的出射波在 （x，y，z）點的入射波方向.

脈沖響應(yīng)Σω＞0Re（～P＋（x，y，z；ω））在空間中僅存在一個波前面，對沒有波前面的空間位置，不能用波前面來確定炮點入射波的方向.若僅對主頻分量計算脈沖響應(yīng)，即僅計算Re（～P＋（x，y，z；ωd）），則波前面將分布于整個區(qū)域，易于求得所有成像點的入射波方向，而此時計算脈沖響應(yīng)的計算量也大幅減少.僅考慮主頻分量，按照上述思路，對于x，y，z三個方向分別有：

由此可得到 （xs，ys，0）處炮點的出射波在 （x，y，z）點的入射波之波前面法線［7］：

2.2 基于局部結(jié)構(gòu)張量分析的地層傾角估算

基于局部結(jié)構(gòu)的張量分析（local structure tensors）為簡單信號提供了一個較簡潔的描述其方向的方法.這里所謂的簡單信號，是指該信號在局部范圍內(nèi)只能沿一個方向變化.對于這類信號，可以把它的局部方向定義成一個張量Ts，該張量的大小是由信號的維數(shù)來確定的.比如，對于3D信號，其局部方向可以表示為3×3的張量.Van Vliet和Verbeek［8］利用一階梯度張量計算了2D情況下的局部地層傾角.在這里，我們沿用一階梯度張量的方法給出3D情況下一階梯度張量表達式：

式中fx，fy，fz分別表示偏移數(shù)據(jù)體沿x，y，z方向上的一階導(dǎo)數(shù)，〈〉表示高斯平滑.我們依據(jù)這一張量可以獲取地下介質(zhì)某點所在地層的切平面及其法線方向：矩陣TS的最大特征值對應(yīng)的特征向量描述了我們所需要的地層法線方向；而其它兩個特征向量組成的平面即為地層的切平面.在炮域波動方程疊前深度偏移提取角道集的研究中，首先將所有單炮的疊前深度偏移結(jié)果疊加得到整個工區(qū)的偏移數(shù)據(jù)體，然后應(yīng)用該偏移數(shù)據(jù)體采用這一方法計算得到地下介質(zhì) （x，y，z）點處地層界面法向量為n（x，y，z）＝（cosαx，cosαy，cosαz），則結(jié)合2.1節(jié)得到的炮點 （xs，ys，0）處出射波在 （x，y，z）點的入射波之波前面法線，可計算入射地震波與地層界面法線的夾角：

3 “保幅”角道集的形成

就反射波偏移成像而言，震源入射波場與反射界面法線的夾角，即入射角，就是反映偏移道集AVA特征的角度.基于此角度的角道集完全反映了由界面波阻抗和泊松比對比決定的隨角度變化的界面反射特征［9－10］.獲得此類角道集是正確應(yīng)用反射幅值特性的基礎(chǔ)，也一直是各種偏移方法的努力方向.

基于第2節(jié)獲得的反射角，可以按反射角對炮域偏移結(jié)果直接累加來得到角道集.有學(xué)者從理論上證明［11］，若按反射角對炮域偏移結(jié)果直接累加，基于相關(guān)成像條件的炮域偏移結(jié)果恰可得到“保幅”的角道集.這里的保幅指的是正確補償了地震波的幾何擴散效應(yīng)，并沒有考慮補償黏性吸收和薄層散射［12］.但實際應(yīng)用中由于炮、檢點的離散分布和非均勻覆蓋，因而難以滿足文獻［11］的應(yīng)用條件.實際上，在地震資料的整個處理流程中，如何更好地解決觀測系統(tǒng)非均勻覆蓋對幅值影響，一直是一個重要的研究課題.

本文因此建議了一種可同時處理幾何擴散效應(yīng)和觀測系統(tǒng)非均勻覆蓋的“保幅”的角道集生成方法.角道集是按反射角對炮域偏移結(jié)果直接累加形成的；在偏移成像中用反褶積成像條件來補償幾何擴散，同時應(yīng)用統(tǒng)計按角度累加過程中落入不同角度區(qū)間內(nèi)的炮數(shù)的方法來剔除觀測系統(tǒng)非均勻覆蓋的影響.

炮域偏移的反褶積成像條件可表達為

式中PU（x，y，z，mΔω）和PD（x，y，z，mΔω）分別是深度延拓后頻率域炮點和檢波點波場，l1Δω和l2Δω分別對應(yīng)地震數(shù)據(jù)有效頻帶的下、上限；令T0是地震記錄的時長，則頻率采樣Δω＝2π／T0.引入?yún)?shù)β是為了保持地震子波不變，三維偏移時β＝j(luò) mΔω，二維偏移時β＝.式（8）的反褶積成像條件可正確處理地震波傳播的幾何擴散效應(yīng).類似式（8）的反褶積成像條件在波動方程偏移發(fā)展的初期就被提出了，由于實際應(yīng)用中除法經(jīng)常導(dǎo)致不穩(wěn)定，因此更多使用的還是相關(guān)成像條件.文獻［13］建議了一個震源波場平滑的反褶積成像條件應(yīng)用方法.本文則通過引入穩(wěn)定的除法實現(xiàn)算法，保證了反褶積成像條件的穩(wěn)定性.

式（8）中的波場相除可進一步表達為

式中上標＊代表共軛.為避免式（9）中除法計算的不穩(wěn)定，我們在除法中應(yīng)用了如下的級數(shù)展開近似式：

當(dāng)x趨于零，式（10）仍可得到穩(wěn)定的結(jié)果，而當(dāng)0.6＜x≤1.2時，式（10）的近似有很好的精度.利用式（10），我們發(fā)展了穩(wěn)定的相除算法，令式（9）中實數(shù)的分母為A（x，y，z，mΔω），由（10）式可得

式中A0（z）是A（x，y，z，mω）在深度z上的平均值.該式既保證了穩(wěn)定性，又可以得到光滑的相除結(jié)果.將式（11）代入（9）再代入（8），可保證反褶積成像條件的穩(wěn)定.

疊加是壓制隨機噪音的有效手段，對地震資料整個處理流程，疊加起到了關(guān)鍵的作用；但疊加也改變了地震資料的幅值特征，特別是在觀測系統(tǒng)非均勻覆蓋的情況下.解決觀測系統(tǒng)非均勻覆蓋是地震資料處理中一直受到關(guān)注的問題，目前尚沒有非常有效的方法.在生成角度道集的過程中，為獲得較高信噪比的角道集，也需要采用疊加處理，即保證有足夠多的反射角（每個角度對應(yīng)一個炮點）落在一個角度區(qū)間內(nèi).對等間距的角度區(qū)間而言，即使是規(guī)則的三維觀測系統(tǒng)，由于（各向同性情況下）不考慮方位角，落在各個角度區(qū)間內(nèi)的炮點數(shù)也是不相同的.

由于采用反褶積成像條件，可認為每個單炮偏移結(jié)果對應(yīng)某一個角度的（有限帶寬）的反射系數(shù)；因此，在累加過程中統(tǒng)計落入不同角度區(qū)間內(nèi)的炮數(shù)，可籍此剔除非均勻覆蓋對角道集幅值的影響.具體實現(xiàn)如下：定義角度采樣間距Δγ，對每一成像點定義一個一維數(shù)組，數(shù)組長度小于90／Δγ；對全部炮循環(huán)，對每炮的偏移結(jié)果，讀取由上節(jié)得到的該炮點震源對各成像點的地層界面法向量夾角γ，計算l＝1＋int｛γ／Δγ｝，將各成像點的幅值累加到對應(yīng)數(shù)組的l元素上，同時統(tǒng)計累加到元素上的炮數(shù)；完成對全部炮循環(huán)，然后用統(tǒng)計得到的各元素上的累加炮數(shù)除各成像點對應(yīng)數(shù)組中的各元素值.這樣就得到各成像點隨角度變化的成像幅值.需要指出的是，這一操作并不能簡單推廣應(yīng)用到常規(guī)的疊加處理；以炮域偏移的疊加過程為例，由于不能保證任一炮記錄的偏移結(jié)果對某一同相軸有實質(zhì)的貢獻，簡單的統(tǒng)計、剔除并不能消除非均勻覆蓋的影響.

4 計算流程和改善波動方程成像效果

結(jié)合炮域波動方程偏移流程以及集群式并行機的軟硬件分布特點，提出了如圖1所示的一種直接由炮域偏移得到角道集的實現(xiàn)流程.這一流程具有如下幾個特點：首先，就各成像點震源入射波場與地層界面的夾角的計算環(huán)節(jié)而言，因其使用了與波動方程炮域疊前深度偏移處理完全相同的偏移速度模型和單程波算子，因此得到的角度更真實地反映了入射波與反射波之間的夾角；而其所需計算量僅為單個頻率的震源波場延拓加上求導(dǎo)及平滑、插值等處理，與疊前偏移的全頻帶波場延拓相比，幾乎可以忽略；其次，這一環(huán)節(jié)主要的計算（圖1中兩虛線箭頭中間部分所示）均是在集群并行機計算節(jié)點上并行進行的，因而可獲得較高的計算效率.就這一實現(xiàn)流程的內(nèi)存需求而言，各成像點震源入射波場與地層界面的夾角的計算環(huán)節(jié)所需內(nèi)存與單炮偏移基本相同，且這一操作與單炮偏移是先后進行的，不需要額外的內(nèi)存，這克服了空間移動或時間移動的成像條件產(chǎn)生角道集方法巨大內(nèi)存需求導(dǎo)致不能應(yīng)用于大規(guī)模工業(yè)數(shù)據(jù)的障礙.本文流程是先對單炮偏移結(jié)果直接在各個計算節(jié)點部分疊加，然后在主節(jié)點對部分疊加結(jié)果進行疊加產(chǎn)生各成像點的角道集，這一實現(xiàn)方法規(guī)避了首先由炮集偏移結(jié)果抽取共成像點道集，然后再結(jié)合入射角信息將共成像點道集轉(zhuǎn)化為角道集時產(chǎn)生的巨大文件存儲及數(shù)據(jù)交換需求.

圖1 直接由炮集成像結(jié)果得到角道集的實現(xiàn)流程圖Fig.1 Flow chart of angle gather obtained directly from shots gather prestack migration result

與現(xiàn)行各類波動方程偏移生成角道集的方法相比，本文提出的方法幾乎不增加額外計算量和內(nèi)存需求，可容易與現(xiàn)行炮域偏移流程結(jié)合，獲得有較好幅值特征的角道集.本文方法獲得的角道集，除服務(wù)于疊前巖石物理參數(shù)反演，還可用于進一步改善偏移成像效果.

就改進偏移成像效果而言，現(xiàn)行積分法偏移流程都包括一個后續(xù)的剩余動校及拉伸切除步驟；這一步驟將存在微小彎曲的同相軸拉平并切除無法拉平和存在拉伸畸變的部分.盡管剩余動校量較小，但它使得同相軸的波峰和波谷對得更齊，因而在疊加過程中更好地保持了高頻成分；而切除也提高了成像的信噪比.但對常規(guī)的炮域偏移而言，由于產(chǎn)出的是隨炮點變化的CIG道集，同相軸的彎曲沒有一定的規(guī)律，有效成像段并不一定在中心部分，將工業(yè)界經(jīng)常使用的剩余動校及拉伸切除應(yīng)用于炮域波動方程偏移就產(chǎn)生了困難.而提取角道集后使得在炮域波動方程疊前深度偏移處理中應(yīng)用剩余動校和拉伸切除變得與Kirchhoff積分偏移一樣容易.就本文提取的角道集而言，其剩余動校量可定量描述為

式中z是零角度處同相軸的深度，γ是角度.在1.0附近對ρ進行掃描并將角道集沿角度疊加，可得到不同ρ值對應(yīng)的能量譜，由能量譜中拾取最大能量對應(yīng)的ρ，由式（12）計算不同角度處同相軸的移動量Δz，同時將大角度處的無效信息切除，這樣就簡單地實現(xiàn)了剩余動校及拉伸切除.

5 數(shù)值算例

5.1 三維層狀介質(zhì)模型

為驗證本文發(fā)展的方法能夠得到相對保幅的角道集，文中設(shè)計了一個四層的水平層狀均勻速度模型，將層位按照從上到下排列，分別編號為1，2，3，4.各層位深度－厚度及速度密度信息見表1，其中第3層為一個低速薄層.

先應(yīng)用聲波方程正演方法生成理論炮集記錄［14］，然后按照圖1所示流程，可由各單炮偏移結(jié)果直接得到角道集.圖2a給出了典型空間位置上本文方法求得的入射角的等值線圖.圖2b展示了某CDP點未進行覆蓋不均勻補償?shù)慕堑兰粚@一角道集進行覆蓋不均勻補償后則示于圖2c.圖2b的角道集沒有進行覆蓋不均勻補償，因此這一角道集的幅值特征完全受到非均勻覆蓋的控制，小角度照明弱，中等角度照明強；可容易地在三維積分法偏移的CRP道集中觀察到這一特征.而覆蓋不均勻補償后的圖2c則正確反映了界面的AVA反射特征.

表1 水平層狀模型層位深度－厚度、速度及密度信息Table 1 The depth－thickness，velocity and density of 3－D layered model

圖3進一步給出在三個界面處，圖2c所示角道集反映的AVA特性與理論值的對比，可見兩者擬合得很好；大角度出現(xiàn)偏差是因為本文正演結(jié)果的長偏移距數(shù)據(jù)有限導(dǎo)致的.

5.2 三維野外數(shù)據(jù)實例

本節(jié)將文中提出的波動方程疊前深度偏移直接生成角道集流程應(yīng)用于我國東部油田某地區(qū)滿覆蓋11.2km2的三維開發(fā)地震數(shù)據(jù).這一數(shù)據(jù)集觀測系統(tǒng)如下：地表接收數(shù)據(jù)約為3600個共炮點道集，每炮16×168＝2688道，接收點間距inline方向為40m，crossline方向為160m；炮點間距inline方向為120m，crossline方向為40m，記錄長度為6s.

利用圖1所示流程，由該數(shù)據(jù)集3600個炮集偏移結(jié)果直接生成了滿覆蓋區(qū)各成像點的角道集.出于計算效率以及計算精度的考慮，流程中的炮集深度偏移及入射角計算中的波場延拓均采用了最優(yōu)分裂Fourier方法［15－17］.對這一實際資料，我們繼續(xù)應(yīng)用了文中建議的改進偏移成像效果方法，即進行了剩余動校及拉伸切除.

圖4a展示了某成像點初始的角道集，在1.0附近對ρ進行掃描并將角道集沿角度疊加，可得如圖4b所示的能量譜，由能量譜中拾取最大能量對應(yīng)的ρ，由式（12）計算不同角度處同相軸的移動量Δz，可將圖4a所示的角道集轉(zhuǎn)換到圖4c，比較兩個圖可見同相軸變得更加平直.由圖4c可見，大角度處已無有效信息，因此可將它們?nèi)鐖D4d所示那樣切除.圖5展示了滿覆蓋區(qū)由剩余動校前后角道集疊加得到的數(shù)據(jù)體抽取的深度為750m的水平切片對比.比較可見在圖5a顯示的經(jīng)過剩余動校后的水平切片上，信噪比得到較大提高，斷層及斷點顯示更為清晰.

我們將所生成的角道集與這一地區(qū)的17口測井?dāng)?shù)據(jù)進行了AVA特性對比.圖6給出了某一井點處所生成的角度域成像道集與典型層位的AVA特性曲線對比，其中理論曲線的計算參數(shù)中，速度與密度信息是由測井?dāng)?shù)據(jù)得到的，泊松比參數(shù)是由有關(guān)地層的經(jīng)驗參數(shù)給出的.圖6a為該點所有深度（自地面至深度為5000m）的角度域成像道集；圖6b與圖6c為圖中所示的深度處典型層位的AVA特性曲線與依據(jù)測井信息計算的理論曲線對比.考察可見，所生成的角度域道集與理論計算曲線變化趨勢吻合較好，從而表明所生成的角度域成像道集展現(xiàn)了良好的AVA特性.

6 結(jié) 論

本文給出了基于炮域波動方程疊前深度偏移結(jié)果直接生成角道集的方法.該方法可將炮域偏移結(jié)果直接轉(zhuǎn)換成反映界面AVA特征的角道集.文中也分別建議了校正因炮點覆蓋不均勻?qū)е碌某上穹嫡`差的方法和基于角道集改善波動方程深度偏移成像效果的方法.與現(xiàn)行各類提取角道集方法相比，本文方法計算效率高，內(nèi)存需求及存儲需求甚少，可靈活應(yīng)用于二維和三維問題，更適用于工業(yè)規(guī)模資料.文中給出了該方法的實際應(yīng)用流程和三維資料的應(yīng)用實例.理論及實際數(shù)值算例表明，本文方法提取的角道集真實反映了界面的AVA特征，而改善偏移成像效果的方法提高了偏移剖面的信噪比.較大規(guī)模的實際三維資料應(yīng)用表明，本文方法已具備了較好的工業(yè)資料應(yīng)用能力.本文方法可為直接識別油氣的AVA反演提供高質(zhì)量的數(shù)據(jù).

圖6 某井點處角道集與典型層位的AVA特性曲線（a）某井點處本文方法得到的角道集；（b）深度為0.95km處由該角道集提取的AVA特性曲線（實線）與測井?dāng)?shù)據(jù)計算得到的AVA特性曲線（虛線）對比；（c）深度為1.07km處由該角道集提取的AVA特性曲線（實線）與測井?dāng)?shù)據(jù)計算得到的AVA特性曲線（虛線）對比.Fig.6 Angle gather at some well and its AVA characteristic curves（a）Angle gather at some well computed with the proposed scheme；（b）Comparisons of AVA characteristic curves of a typical layer at the depth of 0.95km between the angle gather（solid line）and the well logging data（dashed line）；（c）Comparisons of AVA characteristic curves of a typical layer at the depth of 1.07km between the angle gather（solid line）and the well logging data（dashed line）.

（References）

［1］Xie X，Wu R.Extracting angle domain information from migrated wavefield.72nd Annual International Meeting，SEG，Expanded Abstracts，2002：1360－1363.

［2］Wu R，Chen L.Prestack depth migration in angle－domain using beamlet decomposition：Local image matrix and local AVA.73rd Annual International Meeting，SEG，Expanded Abstracts，2003：973－976

［3］Soubaras R.Angle gathers for shot－record migration by local harmonic decomposition.73rd Annual International Meeting，SEG，Expanded Abstracts，2003：889－892.

［4］Rickett J E，Sava P C.Offset and angle－domain common image point gathers for shot－profile migration.Geophysics，2002，67（3）：883－889.

［5］Sava P C，F(xiàn)omel S.Angle－domain common－image gathers by wavefield continuation methods.Geophysics，2003，68（3）：1065－1074.

［6］Biondo B.Angle－domain common－image gathers from anisotropic migration.Geophysics，2007，72（2）：S81－S91.

［7］張劍鋒，張輝，劉禮農(nóng).單程波算子地震波入射角計算.地球物理學(xué)報，2013，56（3）：953－960.Zhang J F，Zhang H，Liu L N.Incident angle field estimation：a one－way propagator approach.Chinese J.Geophys.（in Chinese），2013，56（3）：953－960.

［8］van Vliet L J，Verbeek P W.Estimators for orientation and anisotropy in digitized images.Proceedings of the First Annual Conference of the Advanced School for Computing and Imaging ASCI′95，Heijen（The Netherlands），1995：442－450.

［9］Aki K，Richards P G.Quantitative Seismology.San Francisco：W.H.Freeman ＆Co.，1980.

［10］de Bruin C G M，Wapenaar C P A，Berkhout A J.Angle dependent reflectivity by means of prestack migration.Geophysics，1990，55：1223－1234.

［11］Zhang Y，Xu S，Bleistein N，et al.True－amplitude，angle－domain，common－image gathers from one－way wave－equation migrations.Geophysics，2007，72（1）：S49－S58.

［12］Zhang J F，Wapenaar C P A.Wavefield extrapolation and amplitude－variation－with－angle migration in highly discontinuous media.Geophysical Journal International，2003，155：327－339.

［13］Antoine G，Alejandro V，Dimitri B，et al.Smoothing imaging condition for shot－profile migration.Geophysics，2007，72（3）：S149－S154.

［14］Gao H，Zhang J.Parallel 3－D simulation of seismic wave propagation in heterogeneous anisotropic media：agrid method approach.Geophysical Journal International，2006，165：875－888.

［15］Liu L N，Zhang J F.3Dwavefield extrapolation with optimum split－step Fourier method.Geophysics，2006，71（3）：T95－T108.

［16］Zhang J F，Liu L N.Optimum split－step Fourier 3Ddepth migration：Developments and practical aspects.Geophysics，2007，72（3）：S167－S175.

［17］張劍鋒，盧寶坤，劉禮農(nóng).波動方程深度偏移的頻率相關(guān)變步長延拓方法.地球物理學(xué)報，2008，51（1）：221－227.Zhang J F，Lu B K，Liu L N.Frequency－dependent varyingstep depth extrapolation scheme for wave equation based migration.Chinese J.Geophys.（in Chinese），2008，51（1）：221－227.