多孔介質(zhì)水平分層海底低頻地震波的數(shù)值模擬

盧再華,張志宏,顧建農(nóng)

(海軍工程大學理學院,湖北武漢 430033)

多孔介質(zhì)水平分層海底低頻地震波的數(shù)值模擬

盧再華,張志宏,顧建農(nóng)

(海軍工程大學理學院,湖北武漢 430033)

艦船低頻輻射噪聲在海底巖土層中引起的彈性波被稱為艦船地震波,可用于識別艦船目標。為獲取淺海厚沉積層環(huán)境下艦船地震波的傳播特性,基于多孔介質(zhì)波動理論建立海底地震波的三維交錯網(wǎng)格有限差分算法,對多孔介質(zhì)海底低頻聲源激勵的海底地震波進行數(shù)值計算。計算結(jié)果表明:近場范圍內(nèi),多孔介質(zhì)海底地震波的波動成分主要有透射快縱波、透射橫波、透射慢縱波和海底界面波；對于遠場地震波而言,海水層的波導效應(yīng)逐漸表現(xiàn)出來,海底地震波的波動成分主要是沿水平方向傳播的簡正波和海底界面波。

振動與波；多孔介質(zhì)；低頻；淺海；地震波；交錯網(wǎng)格

0 引言

航行艦船在海底巖土層中引起的彈性波被稱為艦船地震波,主要由艦船低頻輻射噪聲引起,可用于識別艦船目標[1]。為了獲取艦船地震波的傳播特性,研究人員對低頻聲源激勵的海底地震波進行了理論分析[2-4]。這些理論分析將淺海海底介質(zhì)視為單相的液體或固體彈性介質(zhì),難以較好地反映淺海沉積層的多孔介質(zhì)特性和波動傳播規(guī)律。實際淺海多數(shù)存在較厚的沉積層,中國南海北部大陸架淺海沉積層的試驗結(jié)果表明,沉積層厚約60 m,通常為淤泥、淤泥質(zhì)亞粘土、亞粘土、礫石亞粘土、粉砂等未固結(jié)或半固結(jié)巖土層,是典型的飽水兩相多孔介質(zhì)[5]。為了更好地掌握低頻聲信號在淺海厚沉積層環(huán)境中的傳播規(guī)律,基于多孔介質(zhì)波動理論進行研究是必要的。數(shù)值模擬是研究彈性波波場特征和傳播規(guī)律的重要手段,其中交錯網(wǎng)格有限差分法在石油地球物理勘探等研究領(lǐng)域已獲得較為廣泛的應(yīng)用[6-7]。由于低頻聲源激勵淺海地震波在應(yīng)用領(lǐng)域上存在一定的特殊性,相關(guān)的時域數(shù)值模擬比較少見。據(jù)此,本文采用交錯網(wǎng)格有限差分法,基于多孔介質(zhì)波動理論對低頻聲源激勵的淺海地震波的進行時域內(nèi)的數(shù)值模擬,以進一步明確艦船地震波的形成機理和傳播過程,為艦船地震波在水中目標探測領(lǐng)域的應(yīng)用和發(fā)展提供理論基礎(chǔ)。

1 交錯網(wǎng)格高階有限差分算法

將海底沉積層視為各向同性飽和的孔隙彈性介質(zhì),根據(jù)廣義胡克定律得到應(yīng)力應(yīng)變關(guān)系式為

式中:σij和eij=(ui,j+uj,i)/2分別為固相骨架的應(yīng)力張量和應(yīng)變張量;μ、A、Q、R為Biot多孔介質(zhì)彈性常數(shù);δij是Kronecker符號;固相骨架和孔隙液體的體積應(yīng)變分別為e=ui,i和ε=Ui,i,u和U分別為固相和液相的質(zhì)點位移;s=-φp為孔隙流體的有效壓力,與孔隙率φ和孔隙流體的壓強p有關(guān)。在不考慮外加力作用時,固相的應(yīng)力位移關(guān)系為

液相的應(yīng)力位移關(guān)系則滿足廣義Darcy定律,即

式中:ρ11=(1-φ)ρs+ρa;ρ12=-ρa;ρ22=φρf+ρa; ρa、ρs和ρf分別為固相與流相間耦合密度、固相密度和流相密度;耗散系數(shù)b=ηφ2/k,η為液體的粘性系數(shù),k為滲透系數(shù),變量上方的單圓點和雙圓點分別表示對時間的1階和2階導數(shù)。

由以上兩個應(yīng)力位移關(guān)系式可以得到質(zhì)點速度的1階波動方程,

(5)式～(8)式構(gòu)成速度-應(yīng)力表達的飽和孔隙介質(zhì)1階彈性波動方程。

根據(jù)交錯網(wǎng)格有限差分法原理,利用Taylor公式,可以將上述波動方程中速度或應(yīng)力對時間的1階導數(shù)展開成時間差分近似。以速度vx對時間的1階導數(shù)為例,具有2M階精度的時間差分表達式為

類似地,速度或應(yīng)力對空間的1階導數(shù)也可以用如下具有2N階精度的空間差分來表達,以x方向的導數(shù)為例:

一般地,速度或應(yīng)力對時間的導數(shù)采用2階精度的差分近似比較合適,對空間的導數(shù)則應(yīng)該采用較高的階數(shù),可以明顯地提高算法的精度和抑制數(shù)值計算過程中的頻散現(xiàn)象[8]。據(jù)此,本文采用時間2階、空間10階的交錯網(wǎng)格高階差分算法進行分析。在交錯網(wǎng)格技術(shù)中,速度和應(yīng)力等地震波波場分量在空間位置和時間節(jié)點上相互交錯分布。對于三維問題,波場分量和孔隙介質(zhì)彈性參數(shù)在空間位置上的交錯分布如圖1所示;在時間節(jié)點上應(yīng)力通常位于整數(shù)倍時間步長,而速度需和應(yīng)力錯開半個時間步長。

圖1 交錯網(wǎng)格中波場分量和彈性參數(shù)的分布示意圖Fig.1 Discretized variables and elastic constants on staggered grid

根據(jù)上述約定,以x方向速度為例,(5)式和(6)式的三維交錯網(wǎng)格差分格式如下:

類似地,以x方向正應(yīng)力和液相有效壓力為例, (7)式和(8)式的差分格式如下:

(11)式～(14)式中:速度分量vx()括號內(nèi)的上標n+代表離散的時間節(jié)點為(n+1/2)Δt,下標i+、j和k代表離散的空間位置為(i+1/2)Δx、jΔz和kΔy;Δx、Δy、Δz和Δt分別為空間、時間的離散步長;分別為σxx在x方向的空間差分和σxy在y方向的空間差分,差分算子的上標+和-號代表向前差分和向后差分。其他差分算子同上面約定,其余的速度分量、應(yīng)力分量的差分格式同理可得。

海底地震波在半無限區(qū)域中進行傳播,而在地震波的數(shù)值計算過程中,受計算機存儲容量和速度等運行環(huán)境限制,只能在有限區(qū)域內(nèi)進行計算,這種有限區(qū)域的邊界包括自由表面和人為截斷的邊界(人工邊界)。進行數(shù)值模擬時,必需對自由表面和人工邊界進行處理,否則影響數(shù)值模擬結(jié)果的可靠性。本文采用Collino等提出的完全匹配層(PML)方法[9]實現(xiàn)人工邊界,其基本思想是將地震波場分裂為不同坐標軸方向的幾個部分,在匹配層邊界內(nèi)按不同的坐標軸方向進行衰減。以vx的波動方程為例,根據(jù)PML的波場分裂思路,將vx分裂為、和3個部分,vx=++,上標代表該項只與相應(yīng)方向的空間導數(shù)有關(guān),得到不考慮耗散影響的PML吸收邊界方程如下:

式中:d(x)、d(y)和d(z)分別為匹配層內(nèi)在x、y和z方向的衰減因子;σxx,x表示σxx在x方向的空間導數(shù),s,x表示s在x方向的空間導數(shù),其他以此類推。然后,按照交錯網(wǎng)格有限差分法對上式進行離散,即可得到PML吸收邊界方程的差分格式,限于篇幅不再列出。

在自由表面上,應(yīng)力必須滿足自由表面邊界條件,即所有應(yīng)力分量為0.本文采用應(yīng)力鏡像法,通過在自由面上方設(shè)置虛擬層,將自由面以下的應(yīng)力鏡像到虛擬層中,實現(xiàn)自由面應(yīng)力為0的邊界條件。

海水層視為理想流體介質(zhì),其波動方程可以由多孔介質(zhì)彈性波動方程退化得到,即(1)式中,A取為流體的體積壓縮模量K=ρc2,ρ和c分別為流體密度和聲速;其他彈性常數(shù)μ、Q和R取為0.震源采用余弦包絡(luò)脈沖子波,波形函數(shù)為

式中:fc和tc分別為脈沖子波的中心頻率和帶寬參數(shù)。

2 數(shù)值算例及分析

為了考察本文有限差分算法的合理性,首先對多孔介質(zhì)海底中點源激勵的海底地震波進行了數(shù)值計算(算例1)。如圖2所示,整個計算域的長、寬均為1 200m,深度方向為1 000m,網(wǎng)格大小均為5m,時間步長為0.5ms.計算模型的前后左右和下部設(shè)置PML吸收邊界層,層厚均為10個網(wǎng)格。上部海水層厚度為100m,海水層密度1 000 kg/m3,聲速為1 500m/s.震源采用脹縮應(yīng)力源,波形函數(shù)為余弦包絡(luò)脈沖子波,中心頻率40 Hz,帶寬80 Hz,設(shè)置在模型中心,深度為400 m處。波場快照記錄時間為0.25 s.多孔介質(zhì)海底的孔隙率為0.39,物性參數(shù)如表1所示[10]。

圖2 淺海PML吸收層邊界示意圖Fig.2 Sketch of PML boundary in a shallow sea

表1 Biot雙相各向同性介質(zhì)物性參數(shù)Tab.1 Parameters of a homogeneous saturated porous sediment

圖3 固相豎直速度波場快照(t=0.25 s)Fig.3 Snapshot of solid particle vertical velocity wavefield(t=0.25 s)

圖3和圖4分別為固相和液相豎直速度在Oxz平面內(nèi)的波場快照。由于震源采用脹縮應(yīng)力源,直達波場中沒有產(chǎn)生橫波成分,只存在著兩類縱波,直達快縱波P2和慢縱波SP2(下標2代表波動所在的介質(zhì)層,下同)?？炜v波在固相和流相中的相位相同,而慢縱波在固相和流相中的相位相反。左右和下部的PML吸收邊界較好地實現(xiàn)了快縱波向計算域外的無反射傳播。直達快縱波P2向上傳播到海底后,產(chǎn)生向海水層傳播的透射縱波P2P1和向下傳播的3種轉(zhuǎn)換波,分別為反射快縱波P2P2、反射橫波P2S2和反射慢縱波P2SP2.反射慢縱波在流相中振幅大,而在固相中的振幅較小,即在相同條件下,從流相中更容易觀測到慢縱波。向海水層傳播的透射縱波P2P1經(jīng)由自由面反射后,形成向下傳播的縱波P2P1P1.上述傳播特征與Biot理論得出來的結(jié)論[11]基本一致,說明本文有限差分算法的有效性,對人工邊界、海底界面和自由面邊界的處理也是合理的。

圖4 液相豎直速度波場快照(t=0.25 s)Fig.4 Snapshot of liqiud particle vertical velocity wavefield(t=0.25 s)

將驗證算例1中的脹縮應(yīng)力源改為海水層中的脈沖聲源,聲源深度取為88 m,計算區(qū)域大小再相應(yīng)調(diào)整,即可對低頻聲源在多孔介質(zhì)海底中激勵的海底地震波進行數(shù)值模擬。首先,對聲源附近的近場波動進行了算例2的數(shù)值計算。圖5和圖6分別為固相和液相豎直速度在Oxz平面內(nèi)的波場快照?？梢?脈沖聲源在海水層中激勵出直達縱波P1后, P1經(jīng)海底表面反射產(chǎn)生上行縱波,P1經(jīng)自由面反射形成下行縱波P1P1.在多孔介質(zhì)海底中,直達縱波P1在海底表面發(fā)生透射,產(chǎn)生的轉(zhuǎn)換波主要有透射快縱波P1P2、透射橫波P1S2和透射慢縱波P1SP2.另外,在海底表面附近還形成了海底界面波,其波動能量局限于海底表面附近,隨著遠離海底表面波幅迅速衰減?？梢?近場范圍內(nèi)傳播距離較短,各種波動成分的衰減不大,多孔介質(zhì)海底中地震波的波動成分主要有透射快縱波、透射橫波、透射慢縱波和海底界面波,海水層內(nèi)則主要是經(jīng)海底表面、自由面的反射后形成的上行縱波、下行縱波。

圖5 固相豎直速度波場快照(t=0.13 s)Fig.5 Snapshot of solid particle vertical velocity wavefield(t=0.13 s)

圖6 液相豎直速度波場快照(t=0.13 s)Fig.6 Snapshot of liqiud particle vertical velocity wavefield(t=0.13 s)

算例3對低頻聲源激勵的遠場地震波進行了數(shù)值模擬,聲源的中心頻率為20 Hz,帶寬40 Hz,聲源深度為80m,計算域的水平橫距取為1 800m.圖7和圖8分別為固相的水平速度和豎直速度在Oxz平面內(nèi)的波場快照?？梢?當聲源激勵的波動傳播到較遠距離后,波場快照中主要存在兩個波群。對于橫距1 200m附近的波群,幅值隨著遠離海底表面迅速衰減,應(yīng)為海底界面波。對于橫距1 500 m附近的波群,幅值空間分布規(guī)律和左邊的界面波明顯不同,水平速度的波形基本限制在海水層內(nèi),在海水層中間深度達到最大幅值;豎直速度在海水層中間深度發(fā)生反相,這些幅值的空間分布規(guī)律符合海水層內(nèi)1階簡正波的傳播特征。實際上,當海水層內(nèi)的上行波、下行波傳播到遠場后,波陣面曲率變得很小,沿豎直方向基本失去傳播特性,形成主要沿水平方向進行傳播的波群,通常稱為簡正波或模式波[12-13]。

圖7 固相水平速度波場快照(t=1.1 s)Fig.7 Snapshot of solid particle horizontal velocity wavefield(t=1.1 s)

圖8 固相豎直速度波場快照(t=1.1 s)Fig.8 Snapshot of solid particle vertical velocity wavefield(t=1.1 s)

根據(jù)Biot多孔介質(zhì)理論[11],慢縱波衰減大,只能在近距離傳播。多孔介質(zhì)海底中的透射快縱波、透射橫波為球面波,空間衰減較大,也不利于遠距離傳播。海底界面波局限于海底表面附近傳播,簡正波主要在海水層波導內(nèi)傳播,相對而言能傳播到較遠距離。上述模擬結(jié)果表明,對于遠場地震波而言,海水層的波導效應(yīng)逐漸表現(xiàn)出來,海底地震波的波動成分主要是沿水平方向傳播的簡正波和海底界面波。圖9和圖10分別為固相的水平速度和豎直速度在海底表面Oxy平面內(nèi)的波場快照,圖中兩種波動成分表現(xiàn)得更加明顯。根據(jù)上述3個算例,和單相彈性海底中的地震波傳播相比,多孔介質(zhì)海底的直達波和轉(zhuǎn)換波場中都出現(xiàn)了慢縱波成分,但慢縱波衰減大,只能在近距離傳播。在遠場范圍,兩種海底介質(zhì)條件下地震波的波動成分差別不大,均為沿水平方向傳播的簡正波和海底界面波。

圖9 固相水平速度波場快照(t=1.1 s)Fig.9 Snapshot of solid particle horizontal velocity wavefield(t=1.1 s)

圖10 固相豎直速度波場快照(t=1.1 s)Fig.10 Snapshot of solid particle vertical velocity wavefield(t=1.1 s)

3 結(jié)論

為獲取淺海厚沉積層下艦船地震波的傳播特性,本文基于多孔介質(zhì)波動理論,建立了海底地震波的三維交錯網(wǎng)格有限差分算法。對多孔介質(zhì)海底低頻聲源激勵的海底地震波進行了時域內(nèi)的正演數(shù)值模擬,根據(jù)數(shù)值模擬的結(jié)果,得出以下結(jié)論:

1)近場范圍內(nèi),多孔介質(zhì)海底中地震波的波動成分主要有透射快縱波、透射橫波、透射慢縱波和海底界面波,海水層內(nèi)則主要是經(jīng)海底表面、自由面反射后形成的上行縱波、下行縱波。

2)當海水層內(nèi)的上行波、下行波傳播到遠場后,波陣面曲率變得很小,沿豎直方向基本失去傳播特性,形成主要沿水平方向在海水層中進行傳播的簡正波。

3)慢縱波衰減大,只能在近距離傳播;多孔介質(zhì)海底中的透射快縱波、透射橫波為球面波,空間衰減較大,也不利于遠距離傳播;海底界面波局限于海底表面附近傳播,簡正波主要在海水層波導內(nèi)傳播,相對而言能傳播到較遠距離。

4)對于遠場地震波而言,海水層的波導效應(yīng)逐漸表現(xiàn)出來,海底地震波的波動成分主要是沿水平方向傳播的簡正波和海底界面波。

低頻聲源激勵淺海地震波的波動成分和傳播規(guī)律比較復(fù)雜,本文主要是建立了三維問題的交錯網(wǎng)格有限差分算法,初步分析了近場和遠場海底地震波的波動成分。聲源頻率、聲源深度、海水層厚度、海底軟硬程度和衰減大小等對多孔介質(zhì)海底地震波特性的影響規(guī)律,還需進一步進行數(shù)值計算和分析。典型的三維海底為楔形海底,是海底低頻地震波在工程應(yīng)用領(lǐng)域的一種主要海洋環(huán)境,將本文計算模型推廣應(yīng)用到楔形或傾斜海底的地震波傳播問題,對艦船地震波的理論研究具有重要意義。

References)

[1] 盧再華,張志宏,顧建農(nóng).艦船海底地震波及其在水雷引信中的應(yīng)用[J].水雷戰(zhàn)與艦船防護,2004(4):25-28.

LU Zai-hua,ZHANG Zhi-hong,GU Jian-nong.The ship-induced seismic wave on the seafloor and it's application in water mine fuze[J].Mine Warfare and Ship Self-defence,2004(4):25-28.(in Chinese)

[2] 顏冰,周偉,龔沈光.淺海地震波傳播的簡正波模型[J].武漢理工大學學報,2006,30(5):805-807.

YAN Bing,ZHOU Wei,GONG Shen-guang.A normal mode model for ocean seismo-acoustics in shallow water[J].Journal of Wuhan University of Technology,2006,30(5):805-807.(in Chinese)

[3] 盧再華,張志宏,顧建農(nóng).基于波數(shù)積分的點聲源海底地震波計算方法研究[J].武漢理工大學學報,2011,35(1): 118-121.

LU Zai-hua,ZHANG Zhi-hong,GU Jian-nong.An algorithm for seismic wave on seafloor caused by low frequency point sound source based on wave-number integration technique[J].Journal of Wuhan University of Technology,2011,35(1):118-121.(in Chinese)

[4] 張海剛,樸勝春,楊士莪.水中甚低頻聲源激發(fā)海底地震波的傳播[J].哈爾濱工程大學學報,2010,31(7):879-887.

ZHANG Hai-gang,PIAO Sheng-chun,YANG Shi-e.Propagation of seismic waves caused by underwater infrasound[J].Journal of Harbin Engineer University,2010,31(7):879-887.(in Chinese)

[5] 盧博,李趕先.用聲學三參數(shù)判別海底沉積物性質(zhì)[J].聲學技術(shù),2007,26(1):6-14.

LU Bo,LIGan-xian.Discrimination of seafloor sediment properties from sound velocitywaveform and amplitude[J].Technical A-coustics,2007,26(1):6-14.(in Chinese)

[6] Dai N,Vafidis A,Kanasewich E R.Wave propagation in heterogeneous,porousmedia:a velocity-stress,finite-differencemethod [J].Geophysics,1995,60(2):327-340.

[7] Zeng Y Q,He JQ,Liu Q H.The application of the perfectly matched layer in numericalmodeling of wave propagation in poroelastic media[J].Geophysics,2001,66(4):1258-1266.

[8] 董良國,馬在田,曹景忠.一階彈性波方程交錯網(wǎng)格高階差分解法[J].地球物理學報,2000,43(3):411-419.

DONG Liang-guo,MA Zai-tian,CAO Jing-zhong.A staggeredgrid high-order difference method of first-order elastic equation [J].Chinese Journal of Geophysics,2000,43(3):411-419. (in Chinese)

[9] Collino F,Sogka T.Application of the perfectly matched absorbing layermodel to the linear elasto-dynamic problem in anisotropic hetero-geneousmedia[J].Geophysics,2001,66(1):294-307.

[10] Bouzidi Y.The acoustic reflectivity and transmissivity of liquid saturated Porousmedia:experimental tests of the theoretical concepts[D].Edmonton:University of Alberta,2003.

[11] Biot M A.Theory of propagation of elastic waves in a fluid-saturated porous solid(Ⅰ):low-frequency range[J].JAcoust Soc Am,1956,28:168-178.

[12] Jensen F B,Kuperman W A,Porter M B.Computational ocean acoustics[M].New York:AIP Press,1993.

[13] 李整林,王耀俊,馬力.海底沉積物參數(shù)對淺海中低頻聲傳播的影響[J].聲學學報,2000,25(3):242-247.

LIZheng-lin,WANG Yao-jun,MA Li.Effects of sediment parameters on the low frequency acoustic wave propagation in shallow water[J].Acta Acustica,2000,25(3):242-247.(in Chinese)

A Numerical Simulation of Seism ic W ave Caused by Low Frequency Sound Source in Shallow Sea w ith Thick Porous Sediment by Staggered-grid Finite Difference M ethod

LU Zai-hua,ZHANG Zhi-hong,GU Jian-nong
(College of Science,Naval University of Engineering,Wuhan 430033,Hubei,China)

Elastic wave in the seabed caused by low frequency noise radiated from ship is the so-called ship seismic wave which can be used to identify ship target.In order to obtain the propagation characteristics of ship seismic wave in shallow seawith thick porous sediment,a 3D staggered grid finite difference algorithm for calculating seismic wave at seafloor is presented based on Biotwave theory for porous elastic sediment.Numerical calculation of seismic wave caused by low frequency sound source in a typical shallow sea environment is carried out.The calculated results show that the components of seismic wave in the near field include not only fast P-wave,slow P-wave and S-wave transmitted into the sediment,but also an interface wave which propagates nearby the seafloor.However in the far field,the wave components of seismic wave aremainly normalmodes and interface wave,the former propagatesmainly in the waveguide of sea water and the latter propagates along the seafloor.

vibration and wave;porous sediment;low frequency;shallow sea;seismic wave;staggered grid

P733.23

A

1000-1093(2014)12-2065-07

10.3969/j.issn.1000-1093.2014.12.019

2014-02-24

國家自然科學基金項目(51179195);國防基金項目(513030203-02)

盧再華(1972—),男,副教授,碩士生導師。E-mail:luzaihua01@163.com