起伏海面條件下的反射地震響應

孟祥羽 孫建國 魏脯力 徐則雙

(吉林大學地球探測科學與技術學院，吉林長春 130026)

0 引言

近年來，隨著海上時移地震的興起，隨機時變的起伏海面對反射地震響應的影響逐漸得到重視[1-2]。Laws等[3]分析時移地震成像觀測結果后認為：起伏海面造成的成像誤差不可忽略；隨后，系統地討論了北大西洋起伏海面對地震數據的影響，指出即使在適度平靜的海面之下，原始地震數據的振幅在某些頻段的誤差達1～3dB，相位誤差達10°～20°，疊后地震數據的均方根誤差達5%～10%。

為了研究起伏海面對地震波傳播的影響，人們普遍使海面“凍結”到某一個時刻上——“凍結”海面模型。有多種常用的求解不規則邊界的偏微分方程的方法，包括有限元法、譜元法、Kirchoff近似法、有限差分法等[4]；Robertsson等[5]同時采用譜元法、Kirchoff近似法、有限差分法進行起伏海面地震波數值模擬，并對比、分析了三種方法的優缺點；Ego-rov等[6]采用Kirchoff近似法模擬起伏海面數據并與實際地震數據對比。上述研究的區域和建模所用的海浪譜主要來自北海和北大西洋等海域。起伏海面條件下的反射地震響應主要與海面的起伏形態有關，而海面起伏形態一般受當地的重力加速度、水深、海風和潮汐等因素影響，故在不同的海域，海面起伏呈不同特征。

為了更真實地模擬起伏海面的影響，本文利用兼顧精度與計算效率的有限差分法，采用基于鬼點外推的不等距差分格式，利用中國的海浪譜建立起伏海面模型，以此研究中國海域內起伏海面條件下的反射地震響應。

1 基本原理

1.1 差分離散格式

在不規則邊界處求解聲波方程時，大部分邊界點與網格點并不重合，而是落在兩個網格點之間。這種沒有落在網格點上的邊界點稱之為非正則內點。不等距差分是有限差分的一種網格剖分格式，其算子包含了起伏邊界處的高程信息，算子形態靈活、網格起伏整潔而光滑。

圖1為以A0、B0、C0、D0為中心的不同形態的不等距差分算子(時間二階、空間二階)及其結構。由圖可見： 以D0為中心的算子只需在計算z方向波場二階導數時采用不等距差分格式；以A0、B0、C0為中心的算子需要在計算x和z方向波場二階導數時同時采用不等距差分格式(圖1a)。對于聲波方程

(1)

在x方向,不等距(式(2)上)和均勻網格(式(2)下)差分離散格式的波場二階導數分別為

(2)

可見，不等距差分離散格式二階導數中包含了真實界面的高程信息，均勻網格離散格式二階導數中不包含高程信息。則不等距差分的離散格式(時間二階、空間二階)為

u(i,j,2)=2u(i,j,1)-u(i,j,0)+

s(t)δ(i-i0)δ(j-j0)

(3)

式中u(i,j,0)、u(i,j,1)、u(i,j,2)分別為上一時刻、此時刻、下一時刻空間的波場值。對于不等距差分算子，當u1、u2、u3、u4為正則內點時，那么θ1、θ2、θ3、θ4等于1，此時不等距差分算子轉化為均勻網格差分算子(時間二階、空間二階)。

對于起伏海面這種小尺度不規則邊界的波動方程的求解，若采用均勻網格有限差分離散格式，則擬合出的實際界面為由點A3、B3、B0、C0、C3、D0連接的階梯狀起伏界面；若在起伏海面附近采用不等距差分離散格式，則擬合的實際反射界面為由點A4、A1、A2、B4、B1、C1、C2、D1連接的一條較光滑的折線(圖1a)。與均勻網格有限差分相比，不等距差分能更加真實地擬合復雜邊界形態。因此，在起伏界面附近，采用不等距差分格式，在遠離不規則邊界的計算區域，采用均勻網格差分格式。

圖1 以A0、B0、C0、D0為中心的不同形態的不等距差分算子(時間二階、空間二階)(a)及其結構(b)

1.2 非正則內點處波場計算

不等距差分法的關鍵是計算非正則內點處的波場，使其滿足自由邊界條件。在計算非正則內點處的波場值時，鬼點外推法是一種有效方法。Zhang等[7]利用鬼點外推的不等距差分法實現了復雜起伏地表條件下的最小二乘逆時偏移。

根據流體力學理論，海面和海面以上的波場值為0。為了使地震波場滿足自由邊界條件，在不等距差分格式中，往往假設非正則內點uz或海面以上的點u0的波場值不為0，稱uz或u0這類點為鬼點。以z方向為例，鬼點外推法的主要思想是假設鬼點的波場值不為0，以此求解在海面附近的波場二階導數(圖2a)。根據不規則邊界運算區域內的u1、u2、u3等點的波場值及其一階導數，通過多項式外推法計算鬼點的波場值。在數值模擬過程中，鬼點的波場值僅僅用于計算波場二階導數，而與波場的空間分布無關。由鬼點外推法的波場快照(圖2b)可見，基于鬼點外推法的不等距差分格式的模擬結果符合物理規律。

圖2 鬼點外推示意圖(a)與波場快照(t=0.3s)(b)

構建z方向的多項式

u(z)=az3+bz2+cz+d

(4)

求取非正則內點的波場值。式中：a、b、c、d為與正則內點u1、u2、u3的波場值和一階導數有關的系數；u(z)為非正則內點的波場值。

2 算例

利用不等距差分的離散格式和波場外推方法進行數值模擬，通過分析模擬結果闡述中國海域內起伏海面對地震數據的影響。

2.1 基于起伏海面的層狀速度模型

建立層狀速度模型(圖3)，采用文圣常等[8]提出的海浪譜計算模型建立了有效波高為2m的適度平靜海面起伏模型(見附錄A的圖A1)。由不同海面模型地震記錄發現：①與平靜海面模型地震記錄(圖4b)相比，起伏海面模型地震記錄的同相軸存在“抖動”的現象(圖4a)，二次反射的“抖動”更明顯，這是由于一次反射經過起伏海面的不規則反射造成的。②起伏海面模型與平靜海面模型地震記錄存在差別(圖4c)，這種差別是由海面起伏造成的，具體表現為：對于直達波(圖4c的A區域)，誤差主要分布在炮檢距較小的位置，隨著炮檢距增大，誤差逐漸減小；對于反射波(圖4c的B區域)，誤差均勻分布，且誤差大小與海面起伏形態有關。造成上述誤差分布的原因為：在炮檢距較大處，地震波的入射和反射方向幾乎平行于海面，直達波和鬼波發生相對規則的干涉，起伏海面影響較小；對于反射波，地震波的入射方向幾乎垂直于海面，一次反射和鬼波發生不規則干涉，起伏海面影響較大。

圖3 層狀速度模型

模型橫向有401個采樣點，縱向有301個采樣點，空間中正則內點采樣間距為5m，非正則內點的采樣間距主要與海面高程有關，介于0～5m。拖纜位于水下10m處，震源位于橫向第150個采樣點處，主頻為30Hz，采用PML吸收邊界條件

圖4 不同海面模型地震記錄

通過對比發現： 在時間和空間方向上，平靜海面模型地震反射同相軸光滑且能量均勻(圖5b)，起伏海面模型地震反射同相軸粗糙且能量不均勻(圖5a)，這是由于規則的一次反射與經過起伏海面產生的不規則二次反射相互干涉所致。在起伏海面環境下，由于存在地震數據的“抖動”和能量不均勻分布現象，往往造成壓制鬼波效果不理想、成像信噪比和分辨率低、地下照明不均勻等問題。頻譜分析結果表明：與平靜海面模型的地震數據頻譜(圖6b)相比，起伏海面模型地震數據頻譜的頻帶寬窄和能量分布發生變化(圖6a)，這種變化與海面起伏形態有關。

圖5 圖4a的C區域(a)與圖4b的D區域(b)地震數據

圖6 圖4a(a)、圖4b(b)的頻譜

2.2 單道直達波、反射波分析

抽取第80道直達波和反射波數據(圖7)直觀地分析起伏海面對觀測地震數據的影響。由圖7a、圖7c可見：起伏海面的直達波振幅與平靜海面差別較小，且這種差別與拖纜深度無關；起伏海面的直達波走時與平靜海面幾乎相同。原因為：對于直達波的虛反射來說，聲波波長遠大于海面的起伏高度，且入射方向平行于海面，由起伏海面產生的向其他方向的散射能量主要集中于鏡反射方向[9]，故虛反射影響較小，因此虛反射與直達波發生干涉后對直達波的影響較小。由圖7b、圖7d可見：起伏海面的振幅、走時與平靜海面差別較大，且這種差別與拖纜深度有關，這是造成起伏海面條件下鬼波壓制效果不理想的主要因素之一。原因為：起伏海面改變了虛反射的路徑，因此虛反射和一次反射的互相干涉改變了一次反射的波形；隨著拖纜深度增加，虛反射和一次反射的走時時差增大，造成一次反射和虛反射的互相干涉減弱。綜上所述，隨著拖纜深度的不同，一次反射和不規則的虛反射產生了不同的干涉，導致起伏海面對地震數據的影響程度不同。

圖7 第80道地震數據

2.3 單道直達波、反射波振幅譜、相位譜分析

為了驗證起伏海面對地震波頻譜的影響，利用正弦函數建立海面速度模型(圖8)，對單道地震記錄的直達波和反射波分別進行振幅譜(圖9)、相位譜(圖10)分析。由圖9可見，振幅譜的影響主要體現在不同頻率地震波的振幅差異，表現為：①波峰地震記錄的直達波能量較大，頻帶變寬；波谷地震記錄的直達波能量較小，頻帶變窄(圖9a)。②當目標檢波器上方為波峰時，反射波頻譜面積增大，總能量增大，即相對于平靜海面，頻譜低頻能量增大，高頻能量減小，向低頻移動；當目標檢波器上方為波谷時，反射波頻譜面積減小，總能量減小，即相對于平靜海面，頻譜低頻能量減小，高頻能量增大，向高頻移動(圖9b)。

相位譜的影響主要體現在不同頻率地震波的走時差異。由圖10可見：起伏海面對直達波的相位譜影響較小，只在高頻部分有一定影響(圖10a)；無論低、高頻成分，起伏海面對反射波的相位譜影響較大(圖10b)。

圖8 海面起伏呈正弦函數的速度模型

圖9 直達波(a)、反射波(b)振幅譜分析

圖10 直達波(a)、反射波(b)相位譜分析

2.4 基于起伏海面的Marmousi速度模型

建立Marmousi速度模型(圖11)，圖12為Marmousi模型地震記錄。由圖可見：由起伏海面產生的同相軸抖動現象在所有時刻都存在，且與海面的起伏形態相關(圖12a)。圖13為Marmousi模型地震記錄振幅譜分析。由圖可見：與平靜海面地震數據相比，起伏海面地震數據在某些地震道的低頻段(圖13c的橢圓區域)能量較大，這種分布與海面起伏形態具有一定相關性，在其他區域能量較小；與平靜海面地震數據相比，起伏海面地震數據的高頻能量損失大于低頻能量，致使頻帶變窄(圖13d)，意味著地震子波分辨率降低，因此對拖纜資料的寬頻處理和鬼波壓制提出了新的挑戰[10-15]。

圖14為Marmousi模型地震記錄相位譜分析。由圖14c可見，相位譜差異主要分布于較小炮檢距道。原因為：當炮檢距較小時，地震波近似垂直入射于海面，對地震波走時的影響較大；當炮檢距較大時，射向海面的地震波入射角較大，對地震波走時的影響較小。由圖14d可見，起伏海面對地震數據的平均相位譜也有一定影響。圖15為起伏海面與平靜海面地震數據的振幅差和相位差。由圖可見：在有效波高為2m的條件下，起伏海面在某些頻段引橫向有401個采樣點，縱向有301個采樣點，空間采樣間距為5m，拖纜位于水下10m處，震源位于橫向第5個采樣點處，主頻為30Hz，采用PML吸收邊界條件起的振幅差可達1.0～2.5dB(圖15a)，在頻率為50Hz時相位差可達20°(圖15b)。

圖11 Marmousi速度模型

圖12 Marmousi模型地震記錄

圖13 Marmousi模型地震記錄振幅譜分析

圖14 Marmousi模型地震記錄相位譜分析

圖15 起伏海面與平靜海面地震數據的振幅差(a)和相位差(b)

3 結束語

采用文圣常等[8]提出的海浪譜計算模型建立了有效波高為2m的適度平靜海面起伏模型，通過對簡單的層狀模型和Marmousi模型的數值模擬，得到以下認識：

(1)受起伏海面影響，鬼波與一次反射發生不規則的干涉而產生不同的反射地震響應。

(2)由起伏海面產生的地震反射同相軸“抖動”現象在所有時刻都存在，且與海面的起伏形態相關。

(3)由于二次反射的不均勻照明及其與一次反射的不規則干涉，造成地震反射能量的空間分布不均勻現象，影響直達波和反射波的走時、波形、振幅、頻譜、相位和帶寬等，造成鬼波壓制效果不佳、成像信噪比低、分辨率低、局部反演結果不收斂等一系列實際問題。

附錄A 基于海浪譜的起伏海面建模

針對起伏海面的建模與仿真問題，人們常用海浪譜模擬海面的起伏形態。早在20世紀40年代，人們就把海面的起伏看成由不同頻率、不同振幅、不同相位的波動的疊加結果。海面起伏是一個時變過程，具有一定的隨機性。Pierson等[16]用隨機過程分析海浪(譜分析法)。在某一點

(A-1)

式中：η(t)為質點在t時刻的振動幅度；an、ωn分別為振幅和頻率；εn為在0～2π內隨機分布的初始相位。海浪的能量由

(A-2)

表示。S(ω)可以看成角頻率為ω的海浪的能量，由S(ω)和ω構成的描述海水能量分布的譜線即為海浪譜。某個區域的海浪譜通常是在海浪的不同形成時期和不同海洋環境下，經過多次觀測、擬合得到的經驗性結果；在不同的海區，氣候、風力、潮汐、水深、當地的重力加速度和海水的多種物性參數等都會對海面的起伏形態產生一定的影響，因此，不同海域的海浪譜一般不同。

Neumann于20世紀50年代根據經驗提出了Neumann譜[17]

(A-3)

式中：c=3.05；g為當地的重力加速度；w為高于海面7.5m處的風速。

Pierson等[18]對北大西洋充分成長型海浪進行觀測，并提出PM譜

(A-4)

式中：a=0.0081；β=0.74；g為重力加速度；w為高于海面19.5m處的風速。

中國研究者通過對渤海資料的觀測，提出了一種符合中國海域的海浪譜，并與JONSWAP譜進行對比[19-20]。文圣常等[8]在廣泛研究中國海域資料的基礎上，對前人提出的譜進行改進和簡化，提出了一種新海浪譜

(A-5)

ω2(k)=gk

(A-6)

把海浪譜寫成波數域的形式，將波數域的海浪譜乘以一組高斯分布的隨機復數，然后對其進行傅里葉逆變換，便得到不同海況下的隨機起伏海面模型(圖A1)。

圖A1 海面起伏模型