非水平海底情況下海底地震波時域有限差分數值模擬*

王穎 王學鋒 周士弘 趙晨 趙俊鵬 楊勇

1) (北京航天控制儀器研究所，研發中心，北京 100094)

2) (北京市光纖傳感系統工程技術研究中心，北京 100094)

3) (中國科學院聲學研究所，聲場聲信息國家重點實驗室，北京 100190)

研究復雜海洋環境中海底地震波激發及其傳播特性,對海底物理力學特性研究、資源勘探等具有重要的意義.目前針對時域海底地震波的研究大都局限于水平分層的情況,而實際的海底地質條件比較復雜,基于理想環境假設得出的數值解與實際差別較大.本文在考慮傾斜、隆起等非水平海底模型的情形下,采用時間2 階精度、空間10 階精度的交錯網格有限差分方法,同時結合多軸完全匹配層邊界條件,對復雜海洋環境下的海底地震波進行時域數值模擬與分析.利用計算得到的聲場時域波形,分析了復雜海洋環境下海底地震波的傳播特性.結果表明,采用空間高階精度的交錯網格有限差分方法,可改善數值計算中的頻散問題;同時采用多軸完全匹配層替代傳統的完全匹配層,解決了液-固介質中遠距離聲場數值模擬不穩定的問題.在含傾斜與隆起構造的復雜海底模型中,海底基巖隆起改變了Scholte 波的傳播方向,更有利于在較淺深度處接收到Scholte 波.

1 引言

海深、海底和海水聲速剖面等海洋聲學環境的復雜性,對聲傳播影響很大,導致海洋中的聲傳播過程非常復雜[1,2].特別是低頻/甚低頻聲波,由于其具有很強的穿透性,能量可透射到海底介質中,因此聲波的傳播會受到海底彈性基巖層特性的影響.建立符合實際海底環境的聲傳播模型,對于復雜的海洋環境中的聲傳播特性研究、海底地質反演、海底石油與天然氣的勘探等方面有重要的意義.

近幾十年來,壓縮波和剪切波在液-固界面處干涉形成的Scholte 界面波為檢測和反演介質屬性提供了新的思路[3].Scholte 波是一種沿液-固耦合界面傳播的面波[4],與體波相比,Scholte 界面波具有能量強、易識別、衰減慢、傳播距離遠、損耗小、抗干擾能力強等特點.液-固界面Scholte 波的頻散特征與海底介質特性之間存在密切聯系,利用界面波的頻散特性可有效反演海底沉積物參數以及海底淺層地質結構,因此,研究Scholte 界面波的激發、傳播、檢測與應用對海洋聲學與海洋地球物理等領域具有重要的理論意義和應用價值[5].盧再華等[6]將Scholte 界面波簡化為液-固多層半無限空間低頻點聲源引起的地震波動問題,基于波數積分方法,通過FFP 數值積分得到了海底表面聲壓、位移和加速度的頻率特性曲線,分析了不同淺海環境對點聲源海底地震波的波動特征的影響.盧再華等[7]采用大型有限元軟件ANSYS 中的液-固耦合計算模塊,結合黏彈性人工邊界條件,對水平分層海洋環境下低頻點聲源海底地震波進行了數值計算.韓慶邦等[8]分析了Scholte 波在兩相流體與多孔介質固體界面處的傳播特性,以及Scholte 波速與兩相流體積含量、粒徑等介質屬性的關系.左雷等[9]研究了不同海底介質下淺海海底地震波的傳播機理,分析了硬質/軟質海底條件下這種波的存在條件及其波動成分,給出了質點的位移分布和運動軌跡.

以上研究主要在頻域內進行,通過計算地震波場的傳播損失和頻率特性曲線獲取海底地震波的傳播特性,但頻域內的計算不能直觀地區分海底地震波中縱波、橫波和表面波等波動成分.時域的場量波場快照能比較直觀地顯示海底地震波的傳播及波動成分組成,而且時域波形也能為計算波場中各成分的傳播速度提供有效途徑[10].因此,有必要在時域對低頻點聲源激發的海底地震波進行波形的數值計算,以便更直觀地掌握在復雜海洋環境中海底地震波的波動特征和傳播規律.

盧再華等[11]基于有限元數值計算模型,對二維半無限空間層狀海洋環境極低頻點聲源作用下的海底地震波進行了數值模擬,得到了遠場海底表面處加速度和水壓的時域信號.盧再華等[12]基于多孔介質波動理論建立了海底地震波的三維交錯網格有限差分算法,計算了淺海厚沉積層環境下多孔介質低頻聲源激勵的海底地震波.任波等[13]采用有限元軟件ANSYS/LS-DYNA,對簡單水平分層海洋環境下的海底地震波進行了數值模擬,并分析了基巖介質下海底地震波在液-固界面產生的Scholte 波的特性和傳播規律.但是這些時域海底地震波的研究都僅限于水平分層的情況,而實際的海底地質條件比較復雜[14],多數為不均勻厚度分層的楔形或傾斜海底環境,有時海底還存在不規則洼地或臺地[15,16],上述研究尚未考慮這些比較復雜的實際海底環境.

本文在考慮傾斜、隆起等非水平海底模型的情形下,對復雜海洋環境下的低頻地震波進行時域數值模擬與分析,重點研究Scholte 波在該海洋環境下的傳播.計算過程中采用時間2 階精度、空間10 階精度的交錯網格有限差分方法[17,18],同時結合在液-固介質長時間模擬中具有獨特穩定性優勢的多軸完全匹配層(MPML)邊界條件[19,20],并對計算得到的聲場時域合成波形進行分析,以進一步明確復雜海洋環境下海底地震波的傳播特性,評價海洋環境對Scholte 界面波信號的影響.

2 理論和方法

以波動理論為基礎,建立平面直角坐標系,采用一階速度-應力形式的彈性波方程描述二維液-固多層半無限空間海洋環境模型.在交錯網格中,利用空間10 階、時間2 階精度的有限差分算子對該彈性波方程進行求解.同時,由于計算區域是有限的,在波場模擬中采用MPML 作為吸收邊界條件[21],以壓制截斷邊界造成的人工反射,同時能夠保證長時間數值模擬的穩定性[22].

2.1 波動方程及其有限差分

在直角坐標系中,設海水層表面位于x 軸,H為海水深度,0 ≤z ≤H 的區域為海水介質,z >H的半無限區域為彈性海底介質.對于二維液-固多層半無限空間海洋環境模型,其場量滿足如下所示的一階速度-應力形式的二維彈性波方程(假定體力為零):

其中,x 和z 分別為水平方向和垂直深度方向,vx和vz為x 和z 方向的速度,τxx和 τzz為正應力,τxz為剪應力,λ和μ為拉梅系數,bx和bz表示密度的倒數.

研究采用董良國提出的交錯網格高階差分格式進行網格的離散剖分,采用空間2L 階、時間二階精度差分格式時,方程(1)對應的離散格式為

其中,am為 2L 階空間差分精度的差分系數,且

當使用時間二階、空間十階差分格式時,在最短波長達到6 個網格以上時可以保證在地震波采集時間內不會發生頻散[18].

2.2 自由邊界條件

在自由表面上,應力必須滿足自由表面邊界條件,即所有應力分量為0.本文采用Zeng 等[23]改進的真空法,即對 μ,bx,bz等參數作如下設置:

通過(5)式—(7)式保證了模型上邊界處自由邊界條件的實現.而且通過該方法,在海水與海底交界處,對應的液-固界面處剪應力為零的邊界條件也能自動實現.

2.3 吸收邊界條件

在波場模擬中由于計算區域是有限的,為了壓制截斷邊界造成的人工邊界反射,本文采用MPML 作為吸收邊界條件.該邊界條件作為傳統PML 邊界條件的擴展,在多個正交方向同時引入衰減因子,通過調整不同方向衰減因子的比例系數達到改善MPML 穩定性的目的.以x 方向為例,傳統的PML 方程中,只有一個衰減因子是空間變量x 的函數,即

而在MPML 中,3 個正交方向上的衰減因子同為坐標x 的函數,即

時間域的分裂式MPML 方程表達式如下:

當比例系數均為零時,MPML 方程(10)—(14)退化為傳統的PML 形式.在高泊松比介質中,MPML能保證數值模擬的穩定性.

同樣地,采用空間2L 階精度、時間2 階精度的有限差分方法對上述二維分裂式MPML 方程求解,表達式如下:

3 計算方法驗證

設計了橫向距離較大的兩層介質模型,分別采用不同的空間差分精度(4 階和10 階)以及PML和MPML 邊界條件進行該模型中聲場傳播計算,對比驗證空間高階精度方法抑制頻散的能力以及MPML 邊界條件在穩定性方面的優勢.

3.1 模型介紹

計算采用的模型為橫向狹長結構,橫向距離為10 km,縱向深度為150 m.為了更加清楚地展示該模型結構,僅顯示了橫向距離為500 m 的情形,如圖1 所示,模型上層為海水介質,縱橫波速度分別為vp1=1500 m/s,vs1=0,海水深度為50 m;下層為半無限空間,縱橫波速度分別為vp2=3200 m/s,vs2=1800 m/s,厚度100 m.計算時所用的網格大小為dx=2 m,dz=1 m,時間步長為0.2 ms,時間長度為10 s.震源子波主頻為50.0 Hz,其位置坐標為(10 m,45 m) (圖1 中紅色星形),分別在10 m,30 m,45 m,50 m(海水與海底分界面),55 m,70 m等深度處布置接收點(圖1 中紅色倒三角).

圖1 半無限空間海洋環境模型Fig.1.Ocean enviroment model with infinite half-space.

3.2 空間4 階和10 階差分精度數值計算對比

首先,采用空間4 階精度的有限差分方法對該半無限空間模型進行計算,橫向5000 m、深度50 m 處接收到的波形記錄及時頻分析結果如圖2所示.從時域波形記錄可以看出,在記錄的結尾出現了波形抖動,而且時頻分析結果中,Scholte 波所在的能量集中區域也有頻散現象,這與Scholte 波在高頻段非頻散的特性是相悖的.

圖2 空間4 階精度的波形記錄和時頻分析結果 (a)時域波形;(b)時頻分析Fig.2.Seismogram and frequency domain result by method with 4 th order spatial accuracy:(a) Time domain waveform;(b) Timefrequency analysis.

保持空間網格大小不變,將時間步長改變為dt=0.1 ms 或0.05 ms,對模型重新計算,圖3 給出了計算的結果,從圖中可以看出,頻散的現象仍然存在.

圖3 源距5 km 深度 為50 m 的波形記錄和頻域結果 (a)時域波形,時間步長為0.1 ms;(b)時頻 分析,時間步長為0.1 ms;(c)時域波形,時間步長為0.05 ms;(d)時頻分析,時間步長為0.05 msFig.3.Seismogram and frequency domain result of offset 5 km and depth 50 m:(a) Time domain waveform,time step is 0.1 ms;(b) Time-frequency analysis,time step is 0.1 ms;(c) Time Domain waveform,time step is 0.05 ms;(d) Time-frequency analysis,time step is 0.05 ms.

嘗試改變空間網格大小之后,設置dx=1.0 m,dz=1.0 m,圖4 展示了重新計算后的模擬結果,從圖中可以看出,時域波形末端干凈清晰,之前波形中的“尾巴”消失,而頻域分析中高頻段的頻散現象也得到了很好地控制,模擬結果更加符合Scholte 波的頻散特性.總結以上不同參數的計算結果可以發現,空間網格的減小而非時間步長的減小對高頻頻散的改善起到了關鍵的作用,因此,該頻散現象主要歸根于空間網格離散,可定性為空間離散造成的數值頻散,而該頻散問題可通過提高空間差分精度來解決.

圖4 源距5 km 深度為50 m 的波形記錄和頻域結果,空間網格大小為1 m (a)時域波形;(b)時頻分析Fig.4.Seismogram and frequency domain result of offset 5 km and depth 50 m,with spatial interval 1 m:(a) Time domain waveform;(b) Time-frequency analysis.

采用空間高階(取2L=10)精度差分算法對半無限空間模型重新計算,其余參數保持不變,圖5給出了空間10 階精度計算結果,從圖中可以看出,空間10 階精度計算的時域波形更加干凈清晰,無數值頻散造成的“尾巴”,時頻分析的結果也表明采用空間高階差分,Scholte 波的高頻頻散現象也得到了較好的控制.

圖5 空間10 階精度的波形記錄和頻域結果 (a)時域波形;(b)時頻分析Fig.5.Seismogram and frequency domain result by method with 10 th order spatial accuracy:(a) Time domain waveform;(b) Time-frequency analysis.

3.3 PML 與MPML 邊界條件

當交錯網格差分方法的精度和參數設置不變,只改變邊界條件,圖6 給出了分別采用PML 和MPML 邊界條件時,在50 m 深度處接收到的不同橫向位置的時域波形.從圖6 可以看出,在模擬時間較長時(時間長度40 s),采用PML 作為邊界條件,數值模擬出現了不穩定現象,表現為波形幅值特別大,掩蓋了真實的波形信號,導致數值計算失敗.采用MPML 替代傳統的PML,不穩定現象得到控制,因此MPML 有助于改善長時間數值模擬的穩定性.

圖6 深度50 m 處的接收記錄 (a) PML;(b) MPMLFig.6.Record at depth 50 m:(a) PML;(b) MPML.

4 含傾斜構造與隆起的復雜海底模型

4.1 模型介紹

設計了含傾斜構造與隆起的復雜海底模型,以考察傾斜構造及隆起對Scholte 波產生及接收的影響,尤其是海底隆起是否更有利于Scholte 波的接收.模型最大橫向距離為50 km,整體速度模型如圖7(a)所示,其中,在橫向5—35 km 處有傾角3°左右的傾斜構造,在41—44 km 處有海底隆起穿透軟泥層一直向上延伸到海水層底界面(隆起構造放大顯示見圖7(a)).模型自上而下有四層介質組成,各層介質的彈性參數見表1,其中由于非水平界面的存在,各層介質的厚度僅給出了對應的范圍.海水層聲速隨深度呈分段線性變化,在水面處聲速為1540 m/s 然后線性減小,到1000 m 深度處達到最低值1480 m/s,然后海水聲速線性增大,到2000 m 深度處達到1490 m/s,對應的聲速曲線見圖7(b).

表1 含傾斜與隆起的模型參數Table 1.Parameters of model including slope and uplift.

圖7 復雜海底模型 (a) 模型及隆起構造放大顯示;(b)海水層聲速曲線Fig.7.Model with complex seafloor topography:(a) Whole model and zoomed display of uplift structure;(b) Sound speed profile of sea water.

仿真計算時,接收點沿3 個非水平界面(如圖7(a)中白色虛線所示)布置,分別是海水與軟泥層分界面(簡記為界面1)、軟泥層與海底基巖1 分界面(隆起處稍有不同,簡記為界面2)、海底兩層基巖之間分界面(簡記為界面3),模型整個橫向范圍每隔1 km 間距設置接收點,如圖7(a)中白色倒三角所示.由于該模型橫向距離遠大于深度,為了提高計算效率,在橫向和縱向設置不同的網格間距,計算時采用的網格大小為dx=5.0 m,dz=1.0 m,所用Ricker 子波的主頻為15.0 Hz.

4.2 數值結果與分析

(1)傾斜構造下兩種激發-接收模式對比

在傾斜構造存在的情況下,首先測試震源位置不同對Scholte 波接收的影響,以確定有利的震源激發位置.設計了兩種激發模式:1)在模型最左側界面1 上方50 m 深度處激發,即激發-接收沿傾斜界面下傾方向;2)在模型最右側界面1 上方50 m深度處激發,即激發-接收沿傾斜界面上傾方向.接收記錄總長度為40 s.

圖8 展示了兩種激發模式下不同界面處接收到的波形記錄,其中圖8 左欄(圖8(a),圖8(c),圖8(e))、右欄(圖8(b),圖8(d),圖8(f))分別對應震源在模型左、右側的情形.對比相同界面處,左欄中的波場幅值均小于右欄中波場幅值,說明震源設置在模型右側更有利于接收到幅值較強的聲場信號,有助于更好地觀察與分析Scholte 波的傳播.

圖8 不同界面處的接收記錄 (a) 界面1,震源在左;(b) 界面1,震源在右;(c) 界面2,震源在左;(d) 界面2,震源在右;(e) 界面3,震源在左;(f) 界面3,震源在右Fig.8.Records at different interfaces:(a) Interface 1,the source is on the left;(b) Interface 1,the source is on the right;(c) Interface 2,the source is on the left;(d) Interface 2,the source is on the right;(e) Interface 3,the source is on the left;(f) Interface 3,the source is on the right.

另外,觀察每種激發方式中不同界面處接收到的波形記錄,可以看出,每個深度處都可以接收到橫波、縱波等到時較早的波形,而Scholte 波在界面2 處接收到的波形幅值最大,界面1、界面3 兩個深度的信號幅值略小.可見,Scholte 波在界面2(液-固界面)產生,即海底基巖上表面處,向上和向下傳播10 m 離開界面后Scholte 波的幅值迅速減小;而且,觀察同一深度不同源距的Scholte 波,可以看到沿橫向距離傳播,Scholte 波幅值沒有明顯變化,而體波傳播較近的距離之后幅值就衰減較嚴重,這與Scholte 界面波具有衰減慢、傳播距離遠等特點是吻合的.在近源距處,由于水聲反射強度較大,與此處的Scholte 波反射強度差異相對較小,導致Scholte 波在時間或幅值上都不易識別.而傳播較遠距離后,體波衰減大,此時Scholte 波能量相對更高,更容易與地層中的縱橫波反射區分開來.因此,遠距離接收更有利于Scholte 波的檢測與識別.

(2)隆起構造對Scholte 波傳播影響

通過上述仿真與分析得知,相比聲源放置在模型最左側的情形,當聲源放置在模型右側界面1 上方時,更有利于聲信號接收.在此基礎上,采用同樣的震源位置設置方式,仿真分析隆起對Scholte波傳播的影響.隆起構造放大顯示見圖7(a),其橫向范圍為41—44 km,在42—43 km 處隆起頂部與軟泥沉積層上表面平齊.

仿真時,震源位置設置在模型最右側界面1 上方20 m 深度處,沿3 個非水平界面分別布置接收點,接收水平間距為1 km,如圖7 中白色倒三角所示,接收排列總水平長度為50 km,接收記錄時間長度為45 s.由于界面3 處接收到的波形幅值較弱,且該界面是兩層海底基巖間的分界面,對于觀察Scholte 波傳播規律影響較小,因此,僅展示了界面1 和界面2 處各接收點采集的波形記錄,如圖9(a)和(b)所示.由于隆起構造所處橫向距離對應偏移距6—9 km,單獨展示了偏移距為0—10 km時界面1 和界面2 接收到的波形記錄,見圖9(c)和9(d),以便更清晰地觀察隆起構造對Scholte 波傳播的影響.

圖9 不同界面處的接收記錄 (a) 界面1,偏移距 ≤ 50 km;(b) 界面2,偏移距 ≤ 50 km;(c) 界面1,偏移距 ≤ 10 km;(d) 界面2,偏移距 ≤ 10 kmFig.9.Records at different interfaces:(a) Interface 1,offset ≤ 50 km;(b) Interface 2,offset ≤ 50 km;(c) Interface 1,offset ≤10 km;(d) Interface 2,offset ≤ 10 km.

從圖9 可以看出,在界面1 的接收記錄中,相比其他源距,源距7 km(對應模型橫向距離43 km)和8 km (對應模型橫向距離42 km)處的波形幅值顯著增大,而界面2 的接收記錄中,以上兩個源距處的波形幅值則顯著減小.

對模型橫向41 km,42 km,43 km,44 km 四個位置處(分別對應源距為9 km,8 km,7 km,6 km)接收到的波形局部放大進行分析,波形圖如圖10 所示.

圖10 隆起處的波形記錄 (a) 偏移距9 km;(b) 偏移距8 km;(c) 偏移距7 km;(d) 偏移距6 kmFig.10.Seismograms at uplit structure:(a) Offset 9 km;(b) Offset 8 km;(c) Offset 7 km;(d) Offset 6 km.

對于不同界面處接收到的記錄,沿著橫向傳播,大體上都有兩組波形,第一組到時較早,該組波形對應海底基巖1 中的縱橫波、軟泥沉積層中的縱波以及海水中的縱波,隨著傳播距離增大,該組波的能量衰減比較迅速;第二組波形到時較晚,但能量較為集中.從圖10 中可以看出,橫向41 km,42 km,43 km,44 km 四個位置處(對應偏移距分別為9 km,8 km,7 km,6 km)能量最強的波形到時分別為6.9 s、6.2 s、5.4 s、4.6 s,計算出該波形對應的速度大致為1300 m/s,與Scholte 波速度基本一致,從而可以確定圖10 中接收到的波形大部分能量為Scholte 波.

而且,41 km 和44 km 處界面2 接收到的信號幅度較大,42 km 和43 km 處由于海底沉積層的隆起,界面1 在該橫向范圍上成為海水與海底基巖的分界面,取代了界面2 成為Scholte 波產生的有利位置,因此在界面1 處接收到的Scholte 波能量比界面2 處更強;而界面2 在此處相當于已離開Scholte 波產生的界面一段距離,因而Scholte波能量迅速衰減,檢測到的信號幅值減小.在縱向上,Scholte 波在海底基巖上表面產生后,離開界面向上傳播后幅值迅速減小;而沿橫向距離傳播,不同源距處的Scholte 波幅值變化較小.這與Scholte波水平和垂直方向衰減特性是相符的.

形象地說,Scholte 波到了橫向41 km 處開始沿著隆起構造傳播,改變了之前的水平傳播方向,到了44 km 處,隆起消失,波的傳播又回到了隆起右側界面2 的水平方向上.總結來說,海底基巖隆起有利于在較淺深度接收到Scholte 波.

5 結論

本文借助有限差分方法,對傾斜與隆起等復雜構造存在時的低頻海底地震波進行時域數值模擬與分析,主要實現了以下的方法:

1)在對遠距離海洋模型計算的過程中,發現當計算時間較長時,數值模擬容易陷入不穩定的問題.采用多軸完全匹配層替代傳統的PML,將其應用于低頻點源激發的海洋聲場模型的數值計算,解決了遠距離低頻聲場數值模擬不穩定的問題.

2)在低頻聲場計算過程中采用時間2 階精度、空間高階精度的交錯網格有限差分方法,改善數值計算中的頻散問題.

針對含傾斜與隆起構造的復雜海底模型,根據仿真計算結果與分析,可以得出結論:

1)相對于沿傾斜界面下傾方向的激發-接收模式,沿傾斜界面上傾方向激發-接收更有利于聲波的檢測;

2)海底基巖的隆起可改變Scholte 波的傳播方向,更有利于在較淺深度處接收到Scholte 波.