什么是波場延拓 波動方程偏移
反射地震資料的偏移處理方法，波動方程偏移，基于二維地質建模的兩種地震數(shù)值模擬方法的應用及分析。
本文導航
反射地震資料的偏移處理方法
反射地震資料的偏移校正、射線偏移和波動方程偏移等方法統(tǒng)稱偏移處理。偏移處理可使傾斜界面的反射波，斷層面上的斷面波，彎曲界面上的回轉波以及斷點、尖滅點上的繞射波收斂和歸位，得到地下反射界面的真實位置和構造形態(tài)，得到清晰可辨的斷點和尖滅點。因此，偏移處理對提高地震勘探的橫向分辨率具有重要的作用。
偏移處理通常又可稱為偏移歸位、偏移成像、波場延拓成像等。從偏移原理考慮，可將偏移處理分為射線偏移和波動方程偏移兩大類。若以偏移處理的流程，也可分為疊前偏移和疊后偏移。一般對三維資料，應用三維偏移處理，對二維資料則用二維偏移方法處理。以下討論主要以二維偏移方法為主，對三維偏移處理方法，可在二維方法的基礎上擴展即可。
10.4.1 偏移的概念
反射波水平疊加剖面相當于自激自收記錄剖面，在疊加剖面上的反射波同相軸與地下的反射界面有關。當反射界面水平時，反射波同相軸與地下界面形態(tài)一致，如圖10-11（a）所示。當反射界面傾斜時，反射波同相軸則與反射界面形態(tài)不一致，若直接將反射時間作時深轉換，所得視界面為，如圖10-11（b）。與地下真實反射界面 R1 R2 比較，不論是界面長度、界面位置及界面傾角兩者均不一致。視界面相對界面 R1 R2，向界面下傾方向偏移，而且傾角變小。我們稱這種現(xiàn)象為偏移現(xiàn)象，R1 至的水平距離稱為偏移距。偏移現(xiàn)象隨反射界面的埋深和陡度增加而越嚴重。由于偏移現(xiàn)象的存在，當?shù)叵聵嬙鞆碗s時，自激自收剖面上反映的視界面因位置不正確可能產(chǎn)生界面交叉重疊或出現(xiàn)空白帶。如圖10-12，背斜界面出現(xiàn)空白帶，而向斜界面出現(xiàn)界面交叉重疊。
圖10-11 偏移的反射分析
另外，根據(jù)繞射理論，在斷點，尖滅點等巖性突變點還會產(chǎn)生繞射波，這些繞射波再與偏移后的反射波疊加，就使得水平疊加地震剖面上的反射波變得很復雜。若直接用水平疊加剖面解釋地下界面，很難得出正確的結論?？梢娖片F(xiàn)象使地震剖面的橫向分辨率降低。若能使偏移后的波場歸位，繞射波收斂到繞射點，就可恢復反射界面的真實形態(tài)，因此偏移處理就是針對偏移現(xiàn)象的反偏移方法。偏移處理經(jīng)常也簡稱為偏移，但含義和前面提到偏移是不一樣的。
圖10-12 反射界面位置不正確造成空白或干涉
對非零炮檢距地震記錄的偏移處理（疊前偏移），可理解為實現(xiàn)地震波傳播的逆過程，讓波場反向傳播（稱為延拓）。相當于將激發(fā)點和接收點平面逐漸向地下移動，隨之炮檢距也變小。當激發(fā)點和接收點移至反射點時，炮檢距為零，這時的激發(fā)點和接收點位置為偏移處理后的反射點真實位置。
10.4.2 射線偏移方法
射線偏移是建立在幾何地震學基礎上的一類偏移方法?？蓪崿F(xiàn)疊后偏移，也可實現(xiàn)疊前偏移。偏移的基本原理可用地震脈沖的偏移響應來說明，偏移方法可分為圓法偏移、繞射掃描疊加偏移以及橢圓法偏移。
10.4.2.1 偏移脈沖響應
偏移脈沖響應可分為輸入剖面脈沖響應和輸出剖面的脈沖響應。輸入脈沖響應指在輸入的時間剖面（水平疊加剖面）或共炮集記錄上的一個脈沖，在深度剖面上可能存在的界面位置軌跡。例如在均勻介質中，水平疊加時間剖面上的一個脈沖在深度剖面上的偏移脈沖響應是一個圓；而炮集記錄中某一道中的一個脈沖對應的偏移響應是一個橢圓。輸出脈沖響應指在輸出的深度剖面上的一個脈沖，對應時間剖面上的時間軌跡。例如在均勻介質，目標空間有一個脈沖（或繞射點），在自激自收時間剖面上脈沖響應為一繞射雙曲線；在非零炮檢距的炮集記錄上時距曲線也為雙曲線，但兩個雙曲線的計算公式不同。根據(jù)偏移脈沖響應可形成一系列射線偏移方法。
10.4.2.2 繞射掃描疊加偏移
在多種射線偏移方法中，以繞射掃描疊加偏移為例，介紹如下。
利用這一方法作偏移處理時，只考慮幾何關系，將繞射雙曲線上的能量會聚于其頂點。首先，將地下空間劃分為網(wǎng)格，認為每個網(wǎng)格點都是繞射點。根據(jù)網(wǎng)格點坐標計算出它的繞射波時距曲線
勘查技術工程學
式中：（x，z）是繞射點R（地下網(wǎng)格點）的坐標；而（xi，0）是接收點坐標；t0為繞射點至地表的雙程垂直傳播時間。然后按此繞射雙曲線的時距關系ti-xi在實際記錄道上取對應的振幅值，將它們相加后放置在繞射點R處，作為偏移后該點的輸出振幅（圖10-13）。依次對每個網(wǎng)格點都作如上處理就完成了繞射掃描疊加偏移的工作。
如果R點是真正的繞射點（界面點），則按繞射雙曲線取出的各道記錄振幅應當是同相的，它們相加是同相疊加，能量增強，偏移后R點處振幅突出。若R點不是真正的繞射點（非界面點），則參與疊加的幅值是隨機的，疊加結果必然會相互完全抵消或部分抵消，從而使R點處振幅相對較小。因此，偏移后的剖面上，繞射波自動收斂到其繞射點處，在有反射界面處振幅變大，無界面處振幅自然相對減小，顯示出了真實反射界面的位置（繞射雙曲線頂點連線）。
圖10-13 繞射掃描疊加的雙曲線示意圖
10.4.3 波動方程偏移方法
射線偏移是一種近似的幾何偏移，雖然地震波的運動學特點得以恢復，但波的動力學特點（如振幅、波形、相位等）卻受到畸變。因此，射線偏移已逐漸被高精度的波動方程偏移所代替。波動方程偏移是以波動理論為基礎的偏移處理方法，其基本思路是，當?shù)乇懋a(chǎn)生彈性波向下傳播（稱為下行波），遇到反射界面時將產(chǎn)生反射，這時可將反射界面看作新的波源，又有新的波以波動理論向上傳播（稱為上行波），在地表接收到的地震記錄就可看作反射界面產(chǎn)生的波場效應。偏移就是將地表接收到的波場按波動方程的傳播規(guī)律反向向下傳播，通常稱為波場反向延拓，當波場反向延拓到反射界面時成像（成像剖面為偏移剖面），從而找到了真實反射界面，達到了偏移處理的目的?？梢姴▌臃匠唐浦饕刹▓鲅油睾统上駜刹糠纸M成。波場延拓可用多種不同的方法實現(xiàn)，隨之形成了多種不同的波動方程偏移方法。成像也有成像的原理，疊前和疊后偏移各有不同的成像條件。
10.4.3.1 波動方程偏移的成像原理
（1）爆炸反射界面成像原理
該原理屬疊后偏移成像原理。疊加剖面相當自激自收剖面，若將剖面中時間除2，或將傳播速度減一半，就可將自激自收剖面看作在反射界上同時激發(fā)的地震波沿界面法線傳播到地表所接收的記錄，即可將界面看作爆炸源，稱為爆炸反射界面。若用波動方程將地表接收的波場（疊加剖面）作反時間方向傳播（向下延拓），當波場延拓到時間t為零（t=0）時，該波場的所在位置就是反射界面位置。因此，t=0成為疊后波動方程偏移的成像條件。從延拓的結果（地下各點的波場）中取出地下各點處零時刻的波場值組成的剖面就為成像剖面，該剖面為疊后波動方程偏移結果。
（2）波場延拓的時間一致性成像原理
時間一致性成像原理適用于疊前偏移。此成像原理可描述為：在地下某一深度存在一反射界面R（如圖10-14（a）），在地面S點激發(fā)的下行波D到達界面R時產(chǎn)生反射上行波U，到達G點被接收。下行波D到達界R面的時間（或空間位置）與上行波U產(chǎn)生的時間（或空間位置）是一致的，即稱為時間（或空間位置）一致性。設波從S點到R的傳播時間為ts，從R至G的傳播時間為tg，從S到G的總時間為tsg=ts+tg。在疊前偏移中，若模擬一震源函數(shù)D自S點正向（向下）延拓，而將G點接收到的上行波U反向延拓，當D和U延拓深度為Z1時，D的正向傳播時間和U的反向傳播時間分別為ts1和tg1。因Z1ts1說明上行波和下行波所在的時間（或空間位置）不一致（如圖10-14（b）），當D和U延拓深度為zz=ZR時，下行波正向傳播時間為ts1=ts，上行波反向傳播時間為tg2=tg，即有tsg-tg2=ts2，或tsg-tg=ts。這時上、下行波所在的時間（或空間位置）是一致的。再將D、U延拓到Z3，Z3>ZR，即當延拓深度Z>ZR以后，不會再出現(xiàn)時間（或深度位置）一致的現(xiàn)象。在上、下行波延拓過程中，若求下行波場D和上行波場U的零移位互相關，在滿足時間（或空間位置）一致性條件時，相關值最大，而在其他情況相關值很小或為零，延拓過程中的相關結果就為疊前偏移成像剖面。
圖10-14 時間一致性成像原理示意圖
10.4.3.2 有限差分法波動方程偏移
有限差分法波動方程偏移是以地面上獲得的水平疊加時間剖面作為邊界條件，用差分代替微分，對只包含上行波的近似波動方程求解以得到地下界面的真實圖像。這也是一個延拓和成像的過程。
（1）延拓方程
由二維波動方程出發(fā)，
勘查技術工程學
經(jīng)數(shù)學推導，即可得到下面方程：
勘查技術工程學
式中：uxx，uττ，uτt分別表示波場u（x，z，t）的二次導數(shù)。注意，此方程仍然包含了上行波和下行波，仍不能用來進行延拓。
當上行波的傳播方向與垂直方向之間的夾角較小時（小于15°），uττ可以忽略，而對下行波來說，uττ不能忽略。忽略掉uττ項，就得到只包含上行波的近似方程
勘查技術工程學
此即15°近似方程（因為它只適用于夾角小于15°的上行波，或者只有傾角小于15°的界面形成的上行波才能滿足它），為常用的延拓方程。
為了求解此方程還必須給出定解條件。由于震源強度有限，可給出如下定解條件：
1）測線兩端外側的波場為零，即
勘查技術工程學
2）記錄最大時間以外的波場為零，即
勘查技術工程學
3）自激自收記錄（水平疊加剖面）為給定的邊界條件，即時間深度τ=0處的波場值u（x，0，t）已知。
有了這些定解條件就可對方程（10.4-3）求解得到地下任意深度處的波場值u（x，τ，t），這是延拓過程。再根據(jù)前述成像原理，取時間t=0時刻時的波場值，即為時間t=τ時刻的波場值u（x，τ，t）就組成了偏移后的輸出剖面。
（2）差分方程
為了求解微分方程（10.4-3），用差分近似微分，采用如圖10-15所示的12點差分格式，將uxx、uτt表示為差分表達式，可得差分方程
勘查技術工程學
圖10-15 12點差分格式
式中：I和T為向量。
勘查技術工程學
α、β為標量
勘查技術工程學
（3）計算關系和偏移結果
圖10-16畫出了偏移時的計算關系及結果取值位置。A 表示地面觀測到的疊加剖面。由 A 計算下一個深度Δτ處的波場值B，計算 B 時先算第1′排的數(shù)值（只用到 A 中第1排數(shù)值），再算第2′排數(shù)值（要用 A 中第1、2排和 B 中第1′排數(shù)值），依此類推，直到 t=τ為止。再由 B 算下一個深度2Δτ處波場值C……在二維空間（x，t=τ）上呈現(xiàn)出需要的結果剖面信息。
圖10-16 偏移結果取值位置圖
當延拓計算步長Δτ與地震記錄的采樣間隔Δt一樣時，由圖10-16的幾何關系可以看到，偏移剖面是該圖中45°對角線上的值。實際工作中Δτ不一定要與Δt相等，可根據(jù)界面傾角大小確定Δτ，傾角較大時應取較小的Δτ，傾角較小時Δτ可取的較大些，以減少計算工作量。中間值可用插值求得。
與其他波動方程偏移方法相比，有限差分法有能適應橫向速度變化，偏移噪聲小，在剖面信噪比低的情況下也能很好地工作等優(yōu)點。但15°有限差分法對傾角太大的情況不能得到好的偏移效果。因此，相繼又研究發(fā)展了45°、60°有限差分偏移方法和適應更大傾角的高階有限差分分裂算法。
10.4.3.3 克?；舴蚍e分偏移
克?；舴蚍e分偏移是一種基于波動方程克?；舴蚍e分解的偏移方法。
三維縱波波動方程的克?；舴蚍e分解（可見原理部分）為
勘查技術工程學
式中：Q為包圍點（x，y，z）的閉曲面；n為Q的外法線；r為由（x，y，z）點至Q面上各點的距離；〔 〕表示延遲位；〔u〕=u（t-r/v）。
此解的實質是由已知的閉曲面Q上各點波場值計算面內(nèi)任一點處的波場值。它正是惠更斯原理的嚴格數(shù)學形式。
選擇閉曲面Q由一個無限大的平面Q0和一個無限大的半球面Q1所組成。Q1面上各點波場值的面積分對面內(nèi)一點波場函數(shù)的貢獻為零。因此，僅由平地面Q0上各點的波場值計算地下各點的波場值
勘查技術工程學
此時，原公式中的項消失，積分號前的負號也因 z 軸正向與n 相反而變?yōu)檎?/p>
以上是正問題的克?；舴蚍e分計算公式。偏移處理的是反問題，是將反射界面的各點看作為同時激發(fā)上行波的源點，將地面接收點看作為二次震源，將時間“倒退”到t=0時刻，尋找反射界面的源波場函數(shù)，從而確定反射界面。反問題也能用上式求解，差別僅在于〔 〕
不再是延遲位而是超前位，〔u〕=u（t+）。根據(jù)這種理解，克?；舴蚍e分延拓公式應為
勘查技術工程學
按照成像原理，此時t=0時刻的波場值即為偏移結果。只考慮二維偏移，忽略掉y坐標，將空間深度z轉換為時間深度t0=2z/v，得到克?；舴蚍e分偏移公式
勘查技術工程學
式中：τ=[]1/2，xl 為地面記錄道橫坐標，x 為偏移后剖面道橫坐標，r=〔z2+（x-xl）2〕1/2（見圖10-17）。
由=-cosθ，得
勘查技術工程學
由此可見，克?；舴蚍e分偏移與繞射掃描疊加十分相似，都是按雙曲線取值疊加后放在雙曲線頂點處。不同之處在于：①不僅要取各道的幅值，還要取各道的幅值對時間的導數(shù)值參加疊加。②各道相應幅值疊加時不是簡單相加，而是按（10.4-10）式的加權疊加。
圖10-17 克?；舴蚱乒街懈髁渴疽鈭D
正因如此，雖然形式上克?；舴蚍e分法與繞射掃描疊加類似，但二者有著本質區(qū)別。前者的基礎是波動方程，可保留波的動力學特性，后者屬幾何地震學范疇，只保留波的運動學特征。
與其他波動方程偏移法相比，克?；舴蚍e分法具有容易理解，能適應大傾角地層等優(yōu)點。它在速度橫向變化較大的地區(qū)難以使用，且偏移噪聲較大。
波動方程偏移
波動方程偏移與繞射掃描疊加偏移相比有本質上的重大改進，是目前實際生產(chǎn)中使用的主要偏移方法，其中又以15°有限差分偏移最為典型。
1.15°度有限差分法波動方程偏移
15°度有限差分法波動方程偏移是以地面上獲得的水平疊加時間剖面作為邊界條件，用差分代替微分，對只包含上行波的近似波動方程求解得到地下各點波場值，并進而獲得地下界面真實圖像的一種偏移方法。其偏移的過程也是一個延拓和成像的過程。
1）延拓方程的推導
由下述二維波動方程出發(fā)
地震波場與地震勘探
根據(jù)爆炸反射面模型，將速度縮小一半，即用v/2代替v，可得：
地震波場與地震勘探
此方程有二個解，分別對應于上行波和下行波。地震記錄是單純的上行波記錄，故不能用此方程進行延拓，必須將它化為單純的上行波方程才能利用。通常采用的方法是進行坐標變換后取近似。第一步是坐標變換，令
地震波場與地震勘探
上式中第一個變換無任何改變；第二個變換只是將空間深度z換成時間深度τ，也無實質性變化。關鍵是第三個變換，它表示不再用傳統(tǒng)的舊時鐘計時，而是用一個運行速度與舊鐘一樣，但起始時刻各深度不同的新時鐘計時。采用新時鐘計時時，上、下行波就表現(xiàn)出差異。
因為坐標變換并不會改變實際波場，故原坐標系中的波場u（x，z，t）與新坐標系中的波場
是完全一樣的，即
地震波場與地震勘探
由復合函數(shù)微分法，得：
地震波場與地震勘探
將上述二階偏微分結果代入方程（4-4-3），整理后得：
地震波場與地震勘探
為書寫方便，以u、x、t分別代替
，則（4-4-5）式可寫為
地震波場與地震勘探
式中uxx、uττ、uτt分別表示u的二次導數(shù)。注意。此方程仍然包含了上行波和下行波，仍不能用來進行延拓，故還有第二步。
經(jīng)過了坐標變換，雖然波場不變，但在新的坐標系下上、下行波表現(xiàn)出差異，此差異主要表現(xiàn)為uττ的大小不同：當上行波的傳播方向與垂直方向之間的夾角較小時（小于15°），uττ可以忽略；而對下行波來說，uττ不能忽略。忽略掉uττ項，就得到只包含上行波的近似方程
地震波場與地震勘探
此即15°上行波近似方程（因為它只適用于運行方向與垂直方向間的夾角小于15°的上行波，或曰只有傾角小于15°的界面形成的上行反射波才能滿足它），為常用的延拓方程。
為了求解此方程還必須給出定解條件。由于震源強度有限，可以給出如下定解條件：
a.測線兩端外側的波場為零，即
u（x，τ，t）≡0 當 x＞xmax或 x ＜ xmin時
b.記錄最大時間以外的波場為零，即
u（x，τ，t）≡0 當 t＞tmax時
c.自激自收記錄（水平疊加剖面）為給定的邊界條件，即時間深度τ＝0處的波場值u（x，0，t）已知。
圖4-4-5 12點差分格式
有了這些定解條件就可以對方程（4-4-7）式求解得到地下任意深度處的波場值u（x，τ，t），這是延拓過程。再根據(jù)前述成像原則，取傳統(tǒng)舊時鐘零時刻時的波場值，即新時鐘時間t＝τ時刻的波場值u（x，τ，τ）就組成了偏移后的輸出剖面。
2）差分方程的建立
為了求解微分方程（4-4-7）式，用差分近似微分，采用如圖4-4-5所示的12點差分格式，可得：
地震波場與地震勘探
將（4-4-8）式和（4-4-9）式代入（4-4-7）式中得：
地震波場與地震勘探
地震波場與地震勘探
定義向量I、T：
I＝［0，1，0］ T＝［-1，2，-1］
令向量u （x，j，l）為
u（x，j，l）＝［u（i-1，j，l），u（i，j，l），u（i＋1，j，l）］
則（4-4-10）式可簡寫為
地震波場與地震勘探
又令
地震波場與地震勘探
則（4-4-11）式可寫成如下形式：
［I-（α＋β）T］u（x，j＋1，l＋1）-［I＋（α-β）T］u（x，j，l＋1）＋［I-（α＋β）T］u（x，j，l）＝［I＋（α-β）T］u（x，j＋1，l）
因此有：
地震波場與地震勘探
此即適合計算機計算的差分方程。
3）計算步驟和偏移結果
差分方程（4-4-12）形式上是一個隱式方程，即時間深度τ＝（j＋1）Δτ處的波場值不能單獨地用時間深度τ＝jΔτ處的波場值組合得到，方程右邊仍然有τ＝（j＋1）Δτ的項。如圖4-4-6所示，為了求得一排數(shù)據(jù)u（x，j＋1，l），必須用到三排數(shù)據(jù)u（x，j＋1，l＋1），u（x，j，l＋1）和u（x，j，l）。一般來說，隱式方程的求解必須用求解聯(lián)立方程的方法進行，比較麻煩，但這里可以利用有利的定解條件，無須復雜的聯(lián)立運算。
利用定解條件b，在計算新的深度τ＝（j＋1）Δτ處的波場值時，由最大時間開始，首先計算t＝tmax的那一排值。因u（x，j＋1，tmax＋Δt）≡0和u（x，j，tmax＋Δt）≡0，有：
地震波場與地震勘探
計算u（x，j＋1，tmax）只用到已知的u（x，j，tmax）值，十分容易。然后再利用（4-4-12）式遞推地求τ＝（j＋1）Δτ深度處任何時刻的波場值就沒有任何困難了。
具體計算時由地面向下延拓，計算深度Δτ處的波場值：首先計算此深度處在t＝tmax時的波場，然后向t減小的方向進行計算直至本深度處的全部波場值計算完。一個深度的波場值計算結束后，再向下延拓一個步長Δτ繼續(xù)計算。依此類推，可以得到地下所有點在不同時刻的波場值。
如前所述，在新時鐘t＝τ時刻的波場值正是所欲求的“像”。因此，每次遞推計算某一深度τ處的波場值時，由t＝tmax向t減小的方向計算至t＝τ時就可以結束了，u （x，τ，τ）為該深度處的“像”。不同深度處的“像”組成偏移后的輸出剖面。
圖4-4-6 有限差分法偏移求解中的一步
①u（x，j，l＋1），②u（x，j，l），③u（x，j＋1，l＋1），④u（x，j＋1，l）
圖4-4-7 偏移結果取值位置圖
圖4-4-7畫出了偏移時的計算關系及結果取值位置。A表示地面觀測到的疊加剖面，由A計算下一個深度Δτ處的波場值B，計算B時先算第1′排的數(shù)值（只用到A中第1排數(shù)值），再算第2′排數(shù)值（要用到A中第1、2排和B中第1′排的數(shù)值），依此類推進行計算，直到算出t＝Δτ的值為止，再由B計算下一個深度2Δτ處的波場值C，到算出t＝2Δτ的值為止……在二維空間（x，t＝τ）上呈現(xiàn)出需要的結果剖面信息。
當延拓計算步長Δτ與地震記錄的采樣間隔Δt一樣時，由圖4-4-7的幾何關系可以看到，偏移剖面是該圖中45°對角線上的值。實際工作中Δτ不一定要與Δt相等，應當根據(jù)界面傾角大小確定Δτ，傾角較大時應取較小的Δτ，傾角較小時Δτ可取得大一些，以減少計算工作量，中間值用插值方法求得。
與其他波動方程偏移方法相比，有限差分法有能適應橫向速度變化、偏移噪聲小、在剖面信噪比低的情況下也能很好地工作等優(yōu)點，但15°有限差分法在界面傾角太大時不能得到好的偏移效果。因此，又發(fā)展了45°、60°甚至90°的有限差分偏移方法，有興趣的讀者可參閱有關文獻。
2.頻率波數(shù)域波動方程偏移
有限差分偏移方法是在時間空間域中進行計算的。利用傅里葉變換也可以使偏移在頻率波數(shù)域中實現(xiàn)。
與有限差分法偏移的思想完全一樣，認為水平疊加剖面是由界面上無數(shù)震源同時向上發(fā)出的上行波在地面處的波場值u（x，0，t），用它反求地下任一點的波場值u（x，z，t）是延拓過程。再根據(jù)成像原理，取其在t＝0時刻的值u（x，z，0），就組成了偏移后的輸出剖面。
仍由速度減半后的波動方程（4-4-3）出發(fā)，對方程兩邊作關于x和t的二維傅里葉變換，得到一個常微分方程：
地震波場與地震勘探
式中：U＝U（kx，z，ω）是波場函數(shù)u（x，z，t）的二維傅里葉變換，ω＝2πf為角頻率，kx為x方向上的空間波數(shù)。
（4-4-13）式是常微分方程，很容易求解。其解有二個，分別對應于上行波和下行波。偏移研究的是上行波的向下延拓問題，故只考慮上行波解
地震波場與地震勘探
其中U（kx，0，ω）是解的初值，即上行波在地面（z＝0）處記錄的傅里葉變換。因此，式（4-4-14）表示由z＝0 處波場的傅里葉變換求出地下任何深度處波場傅里葉變換的過程，是頻率波數(shù)域中的波場延拓。
通過傅里葉反變換可以由U（kx，z，ω）求出地下任何深度處的波場值：
地震波場與地震勘探
根據(jù)成像原理，偏移結果應是該深度處t＝0時刻的波場值：
地震波場與地震勘探
這就是頻率波數(shù)域偏移的數(shù)學模型。其具體實現(xiàn)步驟就不贅述了。
如果求解常微分方程（4-4-13）時初值不取z＝0處波場值的傅里葉變換，而取任一較淺處的波場傅里葉變換值，則可得到：
地震波場與地震勘探
從而得到相移法偏移的數(shù)學模型：
地震波場與地震勘探
利用此式可逐步向下延拓成像，每延拓一次所用的速度均可改變，所以相移法能夠適應速度的縱向變化。
由于快速傅里葉變換的應用，頻率波數(shù)域偏移法效率十分高，運行時間少，是波動方程偏移算法中最經(jīng)濟的方法，且適用于大傾角地區(qū)。因為計算在頻率波數(shù)域中進行，需要注意假頻問題，且此法對橫向速度變化的地區(qū)不太適應。
3.克?；舴蚍e分偏移
克?；舴蚍e分偏移是一種基于波動方程克希霍夫積分解的偏移方法。
三維縱波波動方程的克希霍夫積分解（見第一章）為
地震波場與地震勘探
式中Q為包圍點（x，y，z）的閉曲面，n為Q的外法線，r為由（x，y，z）點至Q面上各點的距離，［ ］表示延遲位，
。
此解的實質是由已知的閉曲面Q上各點波場值計算面內(nèi)任一點處的波場值，它正是惠更斯原理的嚴格數(shù)學形式。
選擇閉曲面Q由一個無限大的平地面Q0 和一個無限大的半球面Q1 所組成。Q1 面上各點波場值的面積分對面內(nèi)一點波場函數(shù)的貢獻為零，因此僅由平地面Q0 上各點的波場值計算地下各點的波場值。在此條件下，地下任一點的波場值為
地震波場與地震勘探
此時，原公式中的
項消失，積分號前的負號也因z軸正向與n相反而變?yōu)檎?/p>
已知源函數(shù)，求取波傳到某點的波場值是正問題。以上是正問題的克?；舴蚍e分計算公式。偏移處理的是反問題，是將地面接收到的波場值看作為二次震源，將時間“倒退”尋找地下波場值，取t＝0 時刻的波場值確定反射界面的問題。反問題也能用上式求解，差別僅在于［ ］不再是延遲位而是超前位，
。根據(jù)這種理解，克?；舴蚍e分延拓公式應為
地震波場與地震勘探
按照成像原理，t＝0時刻的波場值即為偏移結果。只考慮二維偏移，忽略掉y坐標，將空間深度z轉換為時間深度t0＝2z /v，得到克希霍夫積分偏移公式
地震波場與地震勘探
式中：
；xl為地面記錄道橫坐標；x為偏移后剖面道橫坐標（見圖4-4-8）。
圖4-4-8 克?；舴蚱乒街懈髁渴疽鈭D
由
，得：
地震波場與地震勘探
由此可見，克希霍夫積分偏移與繞射掃描疊加十分相似，都是按照繞射雙曲線取值疊加后放在雙曲線頂點處。不同之處在于：
a.不僅要取各道的幅值，還要取各道幅值對時間的導數(shù)值
參加疊加；
b.各道相應幅值疊加時不是簡單的相加，而是按（4-4-22）式的加權疊加。
正因如此，所以雖然形式上克希霍夫積分偏移法與繞射掃描疊加偏移類似，但二者有著本質的區(qū)別。前者的基礎是波動方程，可保留波的動力學特性；后者屬幾何地震學范疇，只保留波的運動學特征。
與其他波動方程偏移法相比，克希霍夫積分法具有容易理解，能適應大傾角地層等優(yōu)點。它在速度橫向變化較大的地區(qū)難以使用，且偏移噪聲較大。
地震波實際上是在三維空間中傳播的，故要實現(xiàn)完全的偏移必須是三維偏移。目前三維偏移方法已經(jīng)得到了極大的發(fā)展，從二步法到分裂法，到目前已經(jīng)實現(xiàn)了沒有近似的一步法完全三維偏移。
上面介紹的疊后偏移有一個基本假設，即水平疊加剖面是自激自收剖面，實際在地下界面較復雜時這一假設是不成立的。為了實現(xiàn)真正的疊加和偏移，發(fā)展了疊前偏移方法，它將偏移與疊加同時進行，保證了能達到真正的共反射點疊加。但是，一次就完成偏移和疊加的任務對于求取速度參數(shù)是不利的。為此又發(fā)展了疊前部分偏移方法，它只進行部分偏移，使共中心點道集成為真正的共反射點道集，以保證實現(xiàn)共反射點疊加。疊前部分偏移方法加疊后偏移就等于疊前偏移。有了共反射點道集，就好進行速度分析。
目前大部分偏移還是時間偏移，即偏移后得到的是時間剖面。時間偏移沒有考慮地震波穿過界面時的彎折現(xiàn)象。如果考慮這一現(xiàn)象，偏移后得到的就是深度剖面，這種偏移稱為深度偏移。目前，三維疊前深度偏移是最先進的偏移方法。
地震波實際上是彈性波。目前的偏移方法均使用聲波方程，它只是一種聲波偏移，要實現(xiàn)真正完全的偏移應當是彈性波偏移。因為需要多波地震資料，目前彈性波偏移還很少使用。
基于二維地質建模的兩種地震數(shù)值模擬方法的應用及分析
趙忠泉
（廣州海洋地質調查局 廣州 510760）
作者簡介：趙忠泉，男，（1983—），碩士，主要從事海洋油氣資源調查研究工作，E-mail：zzqhello＠163.com。
摘要 利用地震數(shù)值模擬技術結合實際資料，可以建立各種地質體的地震識別模型，有效地避免地震現(xiàn)象的多解性，從而可以提高解釋的精度。本文介紹了二維地質建模的方法流程及兩種模擬方法-褶積法和PSPI波動方程法，前者無邊界條件約束和頻率域中的信號損失，簡潔易行，計算穩(wěn)定，應用廣泛，是最早的地震波場模擬方法；后者通過求解波動方程，包含豐富的波場信息，能夠充分反映地震波的動力學和運動學特征。實際應用中利用褶積法對三維潮道模型及簡化的碳酸鹽巖多旋回傾斜薄互層沉積模型進行了模擬；利用零炮檢距的頻率波數(shù)域的波動方程法模擬了生物礁的地震響應，結果對于碳酸鹽巖生物礁識別有一定指導意義。
關鍵詞 地質建模 數(shù)值模擬 褶積法PSPI法
不同地質體由于其巖性、物性、含油氣性、內(nèi)部結構和巖石組合等的差異，在地震上具有不同的反射特征，包括內(nèi)部結構、外部形態(tài)、振幅、頻率等參數(shù)。由于地震波在地下地質體中傳播的復雜性，加上各種干擾，造成了地震剖面中的各種反射現(xiàn)象存在多解性，大大增加了地震解釋的難度。利用地震數(shù)值模擬技術結合實際資料，在建立不同地質體的地震識別模型的同時也有效地避免地震現(xiàn)象的多解性，從而可以提高解釋的精度。
1 地質建模
地震數(shù)值模擬技術的基礎是地質地球物理模型的建立，可歸結為對地質及地球物理模型結構的數(shù)學描述。
二維封閉結構模型用于建立復雜地質模型。二維封閉結構模型就是定義相同地質屬性為一獨立封閉的地質單元，按照地質屬性將地質模型劃分成多個獨立封閉的地質單元，把所有獨立封閉地質單元按照空間分別有序地排列起來，這樣組成的集合體就構建了一個二維地質模型。封閉結構模型是以積木方式定義地下地質結構，可以描述非常復雜的地質體。二維封閉結構模型被描述為具有相同地質屬性（速度、密度等），并被地層界面、斷層界面或模型邊界所圍成的地質單元的有機組合。對封閉結構模型的描述，實際上就是描述封閉地質單元和封閉地質單元之間的關系，前者包括對封閉地質單元屬性和封閉地質單元邊界的描述；后者是對地質單元空間關系的描述，也就是描述封閉地質單元邊界相接關系及地層屬性［1］。
在進行數(shù)值模擬過程中，為了驗證某些復雜地質體的波場特征，需要繪制多種不同的地質模型，通?？山柚Ｒ?guī)繪圖軟件（繪圖板、Photoshop、CorelDraw，AutoCAD等）繪制好二維封閉結構面，再根據(jù)圖像處理中的區(qū)域填充算法（種子填充和掃描轉換填充），對不同二維封閉結構面進行不同顏色的填充。其中不同顏色代表不同的二維封閉結構面屬性（速度、密度等）；合并相同屬性的封閉面，形成最終的二維封閉結構模型［1］。為了得到二維封閉結構模型的屬性（速度、密度等）模型，需要對二維封閉結構模型的彩色圖進行速度像素空間和屬性空間轉換，根據(jù)顏色空間和屬性空間的相互映射，就可以得到復雜地質體的屬性（速度、密度等）模型，如圖1為模型創(chuàng)建流程圖。
圖1 二維封閉結構模型建立流程圖
2 兩種數(shù)值模擬方法
2.1 褶積模型
在褶積模型中，我們把地震反射信號s（t）看作是地震子波w（t）與地下反射率r（t）的褶積。地震子波w（t），使用實際地震系統(tǒng)記錄到的地下一個單獨的反射界面反射的波形（如圖2，理想的無噪聲褶積過程）。反射率r（t）則代表理想的無噪聲地震記錄。記錄到的地震道s（t）可看作是地震信號w（t）* r（t）與可加噪聲n（t）之和，因此可以把地震道看作是一種有噪聲干擾的，經(jīng)過了濾波的地下反射率的變形。
在無噪聲褶積模型中，我們把地震信號S（t）看作是地震子波w（t）和地下反射系數(shù)r（t）的褶積：
南海地質研究.2012
式中：s（t）——合成地震記錄；
r（t）——反射系數(shù)；
w（t）——地震子波。
圖2 褶積過程
2.2 PSPI波動方程法
通過求解波動方程的數(shù)值模擬方法，能夠充分反映地震波的動力學和運動學特征，波場信息豐富，模擬結果較為準確。這里僅介紹適合橫向速度劇烈變化的頻率-波數(shù)域相移加插值的波場延拓方法［2］。
相位移加插值的波場延拓方法，簡稱PSPI法，基本思想是在波場向下延拓的每個深度步長Δz之內(nèi)，將波場的延拓分成兩部進行，首先用L個參考速度V1，V2，…VL，將位于深度zi處的波場p（x，zi，ω）延拓到zi＋1＝zi＋Δz處，得到L個參考波場p1（x，zi＋1，ω），p2（x，zi＋1，ω），…，PL（x，zi＋1，ω）。第二步，按實際的偏移速度V（x，z）同參考速度V1，V2…，VL的關系，用波場插值的方法求出zi＋1處的波場p（x，zi＋1，ω），按同樣的步驟，可將zi＋1處的波場值p（x，zi＋1，ω）延拓到深度zi＋2，得p（x，zi＋2，ω），直到延拓到最大的深度zmax為止。
對于各向同性介質，取二維標量聲波方程作為延拓的基本方程：
南海地質研究.2012
式中，p＝p（x，z，t）為二維地震波場值；x，z分別為水平方向和垂直方向坐標軸；t為時間軸；v（x，z）為縱、橫向都可變的地震波傳播速度。將式（2）分別對x、t作傅氏變換，考慮到并考慮到?2/?x2和與（-ikx）2和（iw）2的對應關系，可得：
南海地質研究.2012
式中， 是p（x，z，t）的二維傅氏變換；v為地震波速度；w為圓頻率；kx為水平波數(shù)；kz為垂直波數(shù)。零炮檢距情況下的地震記錄模擬只考慮單程波，因此可得到相位移波場延拓公式如下：
南海地質研究.2012
式中， （kx，zi，w）為頻率波數(shù)域波場值；Δz為深度延拓步長；kx為測線方向波數(shù)；kz為深度方向波數(shù)。式（4）為二維波場正演公式，其延拓方向為由地下向地面延拓；式（5）為二維波場偏移公式，其延拓方向為由地面向地下延拓。
為了適應地下地震波場速度在縱橫向均可變的要求，在同一延拓深度內(nèi)用幾個不同地震波速度分別作相移，再用拉格朗日插值公式進行插值，就可求出所有的以不同速度傳播的延拓波場值P（x，zi＋1，t），從而近似地解決了橫向變速時的波場延拓問題［3］。
3 模擬實例
3.1 三維潮道數(shù)值模擬
運用褶積原理建立了一個簡單三維潮道模型，此三維潮道事實上為多個（128）二維剖面排列而成，三維模型的采樣點為128×128×128，利用MATLAB實現(xiàn)。選用子波為雷克（Ricker）子波，其公式為：
南海地質研究.2012
其中fp為主頻。在處理過程中選用主頻為fp＝40 Hz、采樣間隔2 ms，對稱采樣點數(shù)為24，子波波形如圖3。
圖3 雷克子波
圖4 潮道平面圖
圖4為潮道平面圖，該圖僅反映了潮道的平面形態(tài)，作為計算機實現(xiàn)三維建模的邊界控制，橫坐標代表inline線，縱坐標代表xline（crossline）線，圖5為三維地質模型示意圖，模型較簡單，整體由三個水平層疊置而成，在第二層和第三層之間鑲嵌了形如圖4的潮道，此潮道沒有考慮進水方向，根據(jù)此地質模型進行計算機地震正演模擬，可得到相應三維地震數(shù)據(jù)體，從圖中可以看到，黃色虛線（上）和藍色虛線（下）位置上，分別橫跨了三個潮道分支和兩個潮道分支，就是說在相應兩條虛線位置上的兩條測線應該分別有三個和兩個潮道顯示，提取相應的兩條剖面如下圖6和圖7：
圖5 三維地質模型
圖6 xline＝100（黃色虛線）剖面
圖7 xline＝100（藍色虛線）剖面
再在三維數(shù)據(jù)體中沿水平方向做切片，即提取時間切片。圖8為時間切片在地震剖面上的位置示意圖，圖中五條標示線從上到下依次為白色實線、黃色虛線、白色實線、紅色虛線和白色實線，與之對應的時間分別為70 ms、85 ms、95 ms、99 ms和110 ms（時間范圍是0～128 ms），圖9～圖13為相應切片，從圖中可以看出，隨著所做切片時間的增大（深度的增加），潮道的展布范圍逐漸減小，由于地層是水平層狀的，使得時間切片等同于地層切片和沿層切片，其切片效果非常明顯，切片中潮道形態(tài)得到了很好的展示，但是在多個切片中發(fā)現(xiàn)，從可以見到潮道形態(tài)一直到潮道消失的時間范圍是在70～110 ms之間，而潮道的真實范圍是在80～100 ms之間，顯然依據(jù)切片所圈定的潮道的范圍相比真實的范圍擴大了，究其原因是由于不管選取哪一種子波，子波都有一定的延續(xù)長度和有限頻寬，這就限制了合成地震記錄本身的分辨率并不能達到等時厚度反射系數(shù)序列的分辨率。因此在對實際地震資料進行解釋的時候，對地質異常體邊界的識別應該考慮地震子波并非脈沖波所帶來的影響。
圖8 剖面示意圖
圖9 切片t＝70 ms
圖10 切片t＝85 ms
圖11 切片t＝95 ms
圖12 切片t＝99 ms
圖13 切片t＝110 ms
3.2 薄互層沉積模型
圖14為簡化的碳酸鹽巖多旋回傾斜薄互層沉積模型（Zeng，2003），模型簡化是為了更好地突出由巖相控制的波阻抗結構和地震信號之間的相互關系。該模型所有傾斜的傾角都相同，每層都有相同的垂直時間厚度（5 ms或15 m，速度為6000 m/s），泥巖與低孔隙度顆粒灰?guī)r的波阻抗差，以及低孔隙度顆粒灰?guī)r與高孔隙度顆?；?guī)r的波阻抗差都相同，所有高孔隙度顆?；?guī)r具有相同的深度范圍，綜合起來形成了一個水平的巖性地層單元。
其時間域地震響應（圖15）中，高頻情況下（60 Hz雷克子波），地震反射被建設性地調諧到時間地層單元，因此地震同相軸沿著時間地層單元分布（圖15a）。當子波頻率減到40 Hz時，地震反射對時間地層單元和巖性地層單元都有響應（圖15b）。當用30 Hz雷克子波時（圖15c），地震同相軸破壞性地調諧到時間地層單元和建設性地調諧到巖性地層單元，因而時間地層單元的反射進一步變?nèi)?，地震同相軸被巖相反射所控制［4］。
這個模擬過程強調了了解地質格架和時間地層單元以及巖性地層相帶厚度尺度的重要性。時間地層（圖15a）和巖性地層（圖15c）成像都是有用的，前者用于對比，后者用于粗略的儲層評價。然而，這兩種響應不能混淆在一起。圖15b中的兩組相互矛盾的地震同相軸會造成地震假象［4］。
圖14 簡化的碳酸鹽巖多旋回傾斜薄互層沉積模型
3.3 生物礁數(shù)值模擬［5～7］
頻率—波數(shù)域的相移加插值偏移（PSPI）在每一個深度間隔內(nèi)使用多個參考速度進行偏移，由多個偏移結果插值生成最終的偏移剖面，所用插值的速度越多，越能反映實際介質的速度變化情況，此方法在成像精度及橫向變速適應性上具有很大的優(yōu)越性，但處理所需的時間稍長，鑒于本文的二維疊后建模對處理時間沒有過高要求，因此應用PSPI方法做正演、偏移。
圖16為某區(qū)塊過生物礁的原始地震剖面，圖17為根據(jù)此剖面建立的生物礁速度模型：模型速度變化范圍是5600 m/s到5980 m/s，從圖16中可以看出生物礁的底界面清晰可辨，圍巖有披覆現(xiàn)象，內(nèi)部呈雜亂反射。為了檢驗該地質建模的正確性，先采用PSPI方法對該模型進行了波場正演模擬計算，其模擬剖面如圖18所示。由于生物礁埋藏深，生物礁頂?shù)追瓷涞幕《容^大，不規(guī)則點的繞射波雜亂，因此用圖15的速度模型對其進行疊后時間偏移，得到了偏移剖面（圖19），橫向表示256個地震道，縱向表示零偏移距反射時間，礁體最大時間厚度約40 ms。從圖19可以看出，模擬記錄中的礁體頂界與原始剖面有一定差距，但是生物礁底界反射和內(nèi)幕反射以及側翼反射與原始剖面基本一致，其他的地層界面形態(tài)與原始剖面也吻合較好，在一定程度上驗證了地質模型的正確性，說明當生物礁與圍巖之間存在一定波阻抗差異時，在地震剖面上必然出現(xiàn)異常反映，經(jīng)過有效的構造和參數(shù)反演，能夠將其分辨出來。相信通過模型改進以及算法中參數(shù)的調整，能夠與原始剖面更好地吻合，從而為生物礁的地震解釋提供一種有力的驗證工具。
圖15 圖14模型時間域地震響應
4 結論
地震數(shù)值模擬（正演）技術基于地球物理模型的建立，運用概念二維封閉結構地質模型的建立方法，得到復雜地質體的數(shù)學模型，結合各種算法對其進行模擬從而可以驗證相應地質體的地震波場特征；結合實際資料建立不同地質體的地震識別模型，可以有效地減少地震現(xiàn)象的多解性，從而提高解釋的精度；褶積法無邊界條件約束和頻率域中的信號損失，簡潔易行，計算穩(wěn)定，應用廣泛，本文用此方法模擬的偽三維潮道模型及傾斜薄互層模型取得了較好的效果；通過求解波動方程的數(shù)值模擬方法，包含豐富的波場信息，能夠充分反映地震波的動力學和運動學特征，PSPI波場沿拓方法為其中之一，利用正演與偏移相結合的流程模擬了生物礁的地震響應特征，檢驗解釋成果的正確性，為生物礁的地震解釋提供了一種有力的檢驗工具。
圖16 原始剖面
圖17 生物礁地質速度模型（256×256）
圖18 正演記錄（子波主頻30Hz）
圖19 偏移剖面（子波主頻30Hz）
參考文獻
［1］劉遠志.碳酸鹽巖地震相分析與數(shù)值模擬［D］.成都：成都理工大學，2009.
［2］韓建彥.復雜地質體地震正演與偏移［D］.成都：成都理工大學，2008.
［3］賀振華，王才經(jīng)等.反射地震資料偏移處理與反演方法［M］.重慶：重慶大學出版社，1989.
［4］Zeng Hongliu ＆Kerans，C.Seismic frequency control on carbonate seismic stratigraphy；a case study of the Kingdom Abosequence，West Texas，American Association of Petroleum Geologists Bulletin，2003.87，273～293.
［5］賀振華，黃德濟，文曉濤，等.碳酸鹽巖礁灘儲層多尺度高精度地震識別技術［R］.成都：成都理工大學地球探測與信息技術教育部重點實驗室，2009.
［6］熊曉軍，賀振華，黃德濟.生物礁地震響應特征的數(shù)值模擬［J］.石油學報，2009，30（1）：7～65.
［7］熊忠，賀振華，黃德濟.生物礁儲層的地震數(shù)值模擬與響應特征分析［J］.石油天然氣學報，2008，30（1）：75～78
The application and analysis of two kinds of seismic numerical simulation method based on the2D-geological modeling
Zhao Zhongquan
（Guangzhou Marine Geological Survey，Guangzhou，510760）
Abstract：Pick to using seismic numerical simulation technology combined with the actual seismicdata，we can build all kinds of seismic recognition model of geologic body and effectively avoidthe multiple solutions of seismic phenomenon，which can improve the precision of the explana-tion.This paper describes the method of the process of 2D geological modeling and two simulationmethods，seismic convolution method and PSPI wave equation method，the former has no bounda-ry condition and the signal loss in frequency domain，is concise and easy，it can be calculatedsteadily and be applied widely，is the earliest simulation method in seismic wave field，the latterbased on the wave equation，it contains the rich information in wave field，can fully reflect thedynamics and kinematics characteristics of seismic wave.In the practical application，we use theconvolution model in 3D-tidal channel model and the multi-cyclic simplified deposition model oftilt thin interbed layer of carbonate；We simulate the seismic response of reefs using the method ofzero-offset wave equation in frequency and wave number domain，it is confirmed that the resulthas definite significance for the identify of the reef.
Key words：Geological modeling Numerical simulation Convolution PSPI method
