起伏海面條件下的反射地震響應(yīng)

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

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

0 引言

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

為了研究起伏海面對地震波傳播的影響，人們普遍使海面“凍結(jié)”到某一個(gè)時(shí)刻上——“凍結(jié)”海面模型。有多種常用的求解不規(guī)則邊界的偏微分方程的方法，包括有限元法、譜元法、Kirchoff近似法、有限差分法等[4]；Robertsson等[5]同時(shí)采用譜元法、Kirchoff近似法、有限差分法進(jìn)行起伏海面地震波數(shù)值模擬，并對比、分析了三種方法的優(yōu)缺點(diǎn)；Ego-rov等[6]采用Kirchoff近似法模擬起伏海面數(shù)據(jù)并與實(shí)際地震數(shù)據(jù)對比。上述研究的區(qū)域和建模所用的海浪譜主要來自北海和北大西洋等海域。起伏海面條件下的反射地震響應(yīng)主要與海面的起伏形態(tài)有關(guān)，而海面起伏形態(tài)一般受當(dāng)?shù)氐闹亓铀俣取⑺?、海風(fēng)和潮汐等因素影響，故在不同的海域，海面起伏呈不同特征。

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

1 基本原理

1.1 差分離散格式

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

圖1為以A0、B0、C0、D0為中心的不同形態(tài)的不等距差分算子(時(shí)間二階、空間二階)及其結(jié)構(gòu)。由圖可見： 以D0為中心的算子只需在計(jì)算z方向波場二階導(dǎo)數(shù)時(shí)采用不等距差分格式；以A0、B0、C0為中心的算子需要在計(jì)算x和z方向波場二階導(dǎo)數(shù)時(shí)同時(shí)采用不等距差分格式(圖1a)。對于聲波方程

(1)

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

(2)

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

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)分別為上一時(shí)刻、此時(shí)刻、下一時(shí)刻空間的波場值。對于不等距差分算子，當(dāng)u1、u2、u3、u4為正則內(nèi)點(diǎn)時(shí)，那么θ1、θ2、θ3、θ4等于1，此時(shí)不等距差分算子轉(zhuǎn)化為均勻網(wǎng)格差分算子(時(shí)間二階、空間二階)。

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

圖1 以A0、B0、C0、D0為中心的不同形態(tài)的不等距差分算子(時(shí)間二階、空間二階)(a)及其結(jié)構(gòu)(b)

1.2 非正則內(nèi)點(diǎn)處波場計(jì)算

不等距差分法的關(guān)鍵是計(jì)算非正則內(nèi)點(diǎn)處的波場，使其滿足自由邊界條件。在計(jì)算非正則內(nèi)點(diǎn)處的波場值時(shí)，鬼點(diǎn)外推法是一種有效方法。Zhang等[7]利用鬼點(diǎn)外推的不等距差分法實(shí)現(xiàn)了復(fù)雜起伏地表?xiàng)l件下的最小二乘逆時(shí)偏移。

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

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

構(gòu)建z方向的多項(xiàng)式

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

(4)

求取非正則內(nèi)點(diǎn)的波場值。式中：a、b、c、d為與正則內(nèi)點(diǎn)u1、u2、u3的波場值和一階導(dǎo)數(shù)有關(guān)的系數(shù)；u(z)為非正則內(nèi)點(diǎn)的波場值。

2 算例

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

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

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

圖3 層狀速度模型

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

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

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

圖5 圖4a的C區(qū)域(a)與圖4b的D區(qū)域(b)地震數(shù)據(jù)

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

2.2 單道直達(dá)波、反射波分析

抽取第80道直達(dá)波和反射波數(shù)據(jù)(圖7)直觀地分析起伏海面對觀測地震數(shù)據(jù)的影響。由圖7a、圖7c可見：起伏海面的直達(dá)波振幅與平靜海面差別較小，且這種差別與拖纜深度無關(guān)；起伏海面的直達(dá)波走時(shí)與平靜海面幾乎相同。原因?yàn)椋簩τ谥边_(dá)波的虛反射來說，聲波波長遠(yuǎn)大于海面的起伏高度，且入射方向平行于海面，由起伏海面產(chǎn)生的向其他方向的散射能量主要集中于鏡反射方向[9]，故虛反射影響較小，因此虛反射與直達(dá)波發(fā)生干涉后對直達(dá)波的影響較小。由圖7b、圖7d可見：起伏海面的振幅、走時(shí)與平靜海面差別較大，且這種差別與拖纜深度有關(guān)，這是造成起伏海面條件下鬼波壓制效果不理想的主要因素之一。原因?yàn)椋浩鸱Ｃ娓淖兞颂摲瓷涞穆窂?，因此虛反射和一次反射的互相干涉改變了一次反射的波形；隨著拖纜深度增加，虛反射和一次反射的走時(shí)時(shí)差增大，造成一次反射和虛反射的互相干涉減弱。綜上所述，隨著拖纜深度的不同，一次反射和不規(guī)則的虛反射產(chǎn)生了不同的干涉，導(dǎo)致起伏海面對地震數(shù)據(jù)的影響程度不同。

圖7 第80道地震數(shù)據(jù)

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

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

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

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

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

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

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

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

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

圖11 Marmousi速度模型

圖12 Marmousi模型地震記錄

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

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

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

3 結(jié)束語

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

(1)受起伏海面影響，鬼波與一次反射發(fā)生不規(guī)則的干涉而產(chǎn)生不同的反射地震響應(yīng)。

(2)由起伏海面產(chǎn)生的地震反射同相軸“抖動”現(xiàn)象在所有時(shí)刻都存在，且與海面的起伏形態(tài)相關(guān)。

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

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

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

(A-1)

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

(A-2)

表示。S(ω)可以看成角頻率為ω的海浪的能量，由S(ω)和ω構(gòu)成的描述海水能量分布的譜線即為海浪譜。某個(gè)區(qū)域的海浪譜通常是在海浪的不同形成時(shí)期和不同海洋環(huán)境下，經(jīng)過多次觀測、擬合得到的經(jīng)驗(yàn)性結(jié)果；在不同的海區(qū)，氣候、風(fēng)力、潮汐、水深、當(dāng)?shù)氐闹亓铀俣群秃Ｋ亩喾N物性參數(shù)等都會對海面的起伏形態(tài)產(chǎn)生一定的影響，因此，不同海域的海浪譜一般不同。

Neumann于20世紀(jì)50年代根據(jù)經(jīng)驗(yàn)提出了Neumann譜[17]

(A-3)

式中：c=3.05；g為當(dāng)?shù)氐闹亓铀俣?；w為高于海面7.5m處的風(fēng)速。

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

(A-4)

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

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

(A-5)

ω2(k)=gk

(A-6)

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

圖A1 海面起伏模型