適應橫向速度變化的波動方程偏移程序

　　摘 要： 本文實現了時間—空間域和頻率—波數域的地震波場外推，在時間—空間域實現了65度波動方程有限差分偏移，在頻率—波數域實現了只適合地震波速度，僅隨深度變化的相移法偏移和適應地震波速度橫向變化的分裂步有限差分偏移。
　　關鍵詞： 65度方程 有限差分 相移法偏移 分裂步傅立葉偏移（SSF）
　　
　　一、引言
　　1972年克萊鮑特（J.Clearbout）首先提出波動方程有限差分偏移法，給出了至今仍然廣泛應用的15度方程。雖然這種偏移技術只能對傾角很小地質構造進行成像，但是我們今天所用的各種偏移技術都借鑒了克萊鮑特對于波動方程的處理方法。1978年蓋茲達格（Gazdag）提出相移法偏移，同樣是在頻率—波數域偏移，相移法能夠適應地震波速度隨深度變化的情況，也是到目前為止速度僅隨深度變化地質模型的最佳偏移方法，其計算量遠小于波動方程有限差分偏移技術，沒有傾角限制，沒有假頻的危險，可是這種方法還不能適應地震波速度的橫向變化。1990年Stoffa提出分裂步傅立葉變換法（SSF）偏移，也叫相移加校正偏移，對相移法進行改進適應地震波速度橫向變化，在頻率域進行相位移動以后，在頻率—空間域再做一次相位校正。這種方法沒有傾角限制，計算量也較小，但是它要求速度函數是平滑的（一階可導），不能適應地震波速度橫向變化劇烈的地區。
　　我從地震波動理論出發，根據前人所做的工作，分別在時間—空間域和頻率—波數域，實現了單程波動方程地震波場延拓。運用爆炸反射模型的概念，在時間—空間域和頻率—空間域實現了零偏移距（疊后）65度波動方程有限差分偏移。在頻率—波數域實現了相移法偏移，進一步的，根據Stoffa的推導實現了分裂步傅立葉（SSF）偏移以適應地震波速度的橫向變化。在此基礎上，設計了向斜構造地層模型來驗證相移法偏移。
　　二、理論基礎與程序實現方法
　　1.爆炸反射面成像原理
　　圖1顯示了兩種波動傳播情形，第一種是假象的地下反射層突然爆炸的實驗，地表布置一排假想的檢波器，記錄從假想的爆炸反射面向上傳播到地表的波動。
　　圖1在地表所有的點布置炮—檢對觀測回聲信號（左邊）與“爆炸反射面”概念模型（右邊）
　　注意圖中野外記錄情況下的射線路徑與爆炸反射面情形的射線路徑看起來是一樣的。想象這兩種波場（觀測波場和假想波場）是一樣的，這將給我們帶來概念上的方便。如果它們是相同的，我們就可以忽略數以千計的業已做過的試驗，而把注意力集中到一個假想的試驗上。兩種情況的一個明顯的區別是野外觀測系統中的波必須先向下傳播然后又沿相同的路經向上傳播，而假想實驗只需要向上傳播。野外實驗的旅行時除以2，實際上，通常在分析野外實驗數據（雙程走時）時，都假設其聲波速度為實際速度的一半。類比零偏移數據和爆炸反射模型，我們用爆炸反射模型零時刻的波場來表示地下地層的構造圖像。
　　2.波動方程有限差分偏移
　　根據Clerbout的15度近似單程波方程
　　+=0 （1）
　　利用地面記錄的零偏移距地震波場外推得到地下任意點的地震波場，根據爆炸反射面概念，認為零時刻地下任意點的波場振幅即為該點的反射強度，即為地下反射界面圖像。
　　將方程（1）進行差分離散，并根據初始條件和邊界條件：
　　p=φ，z=0的邊界條件 （2）
　　并且，p=p=0側邊界條件（可以用吸收邊界等替代） （3）
　　還有初始條件p=0 （4）
　　我們從j=1，n=N開始，考慮所有的1≤i≤I，以及側邊界條件（3），就可以形成一個三對角方程組：
　　1+2a-a0… 0 0-a1+2a-a … 000 -a1+2a …0 0M MMO M M000… 1+2a -a000… -a 1+2a pppMpp
　　=cccMcc（5）
　　方程（5）可以通過追趕法求解，得到p，1，2，…，I。對所有j和n循環可求出所有p，i=1，2，…，I；j=1，2，…，J；n=1，2，…，N，應用成像條件t′=τ，p即為最終輸出圖像。
　　有限差分45度和65度方程偏移法：
　　考慮（10）式時我們不是直接把?鄣p/?鄣z′忽略，而是進一步對z′求導，利用（10）式消去?鄣p/?鄣z′，然后略去?鄣p/?鄣z′得到45度方程:
　　-+=0 （6）
　　用系數優化的方法把上式中的1/4和1/2分別由0.3767和0.4761代替可以得到65度方程，它們的差分求解方法類似于15度方程。
　　3.速度隨深度變化的相移法偏移
　　相移法用exp（ikz）直接進行向下外推，然后估算t=0時（反射面在t=0時激發）的波場。先對時間剖面進行二維傅立葉變換，然后把所有在（k，ω）平面內的變換后數據值乘以下式進行相位移：
　　е=exp{-i［1-（）］△z} （7）
　　總結起來，相移法的計算步驟為：
　　（1）對輸入疊加剖面進行二維傅立葉變化，把p（x，z，t）變換為p（k，0，ω）；
　　（2）對變換域波場進行相位移，即乘以因子C，即：
　　p（k，j△z，ω）=C×p（k，（j-1）△z，ω），j=1，2，…，J；
　　（3）對所有頻率ω求和，得到p（k，j△z，t=0）=p（k，j△z，ω），j=1，2，…，J；
　　（4）沿k進行反傅立葉變換，得到p（x，z）即為地層映像。
　　4.分裂步傅立葉偏移（SSF）
　　Stoffa把慢速度場（速度的倒數）分為兩部分，一部分稱為參考慢度u，另一部分稱為慢度擾動△u=u（，z）。與速度隨深度變化的相移法一樣，先用參考慢速做一次相移，再在頻率空間域在每一空間位置用慢速擾動做一次相位校正，從而適應地震波速的橫向變化。具體偏移步驟為：
　　（1）首先把地面波場變換到頻率-波數域p（，z=0，ω）；
　　（2）用參考慢度場把地面波場向下延拓，即p（，z+△z，ω）=p（，z，ω）е其中k=ω；
　　（3）把延拓后的波場變換到頻率空間域p（，z+△z，ω）；
　　（4）根據慢度擾動△u再次相移，即p（，z+△z，ω）=p（，z+△z，ω）е；
　　（5）對所有的頻率求和得到地下任意位置的地層圖像p（，z）。