祝捍皓，肖 瑞，朱 軍，唐 駿

(1.浙江海洋大學(xué) 海洋科學(xué)與技術(shù)學(xué)院，浙江 舟山 316022；2.浙江大學(xué) 浙江省海洋觀測(cè)-成像試驗(yàn)區(qū)重點(diǎn)實(shí)驗(yàn)室，浙江 舟山 316021；3.中國(guó)科學(xué)院聲學(xué)研究所 聲場(chǎng)聲信息國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100190；4.浙江海洋大學(xué) 船舶與機(jī)電工程學(xué)院，浙江 舟山 316022；5.挪威科技大學(xué) 聲學(xué)組，挪威 特隆赫姆 7491)

淺海中的聲傳播問(wèn)題是水聲領(lǐng)域研究的熱點(diǎn)之一，是理解、預(yù)測(cè)及應(yīng)用各種淺海聲學(xué)現(xiàn)象的基礎(chǔ)[1-3].隨著潛艇隱身技術(shù)的發(fā)展和對(duì)各類物理海洋現(xiàn)象的監(jiān)測(cè)需求，各類淺海水下聲學(xué)探測(cè)手段逐漸轉(zhuǎn)向低頻(小于1 kHz)，這顯然與低頻信號(hào)在淺海中的傳播問(wèn)題緊密聯(lián)系，因而對(duì)淺海低頻聲傳播問(wèn)題的研究獲得了越來(lái)越廣泛的關(guān)注[4-6].但目前在研究淺海環(huán)境下的聲傳播問(wèn)題時(shí)，多數(shù)只關(guān)注海水層中聲場(chǎng)的計(jì)算建模和傳播規(guī)律，且通常假設(shè)海底為水平分層的液態(tài)介質(zhì)，而有關(guān)考慮海底為彈性介質(zhì)下地形變化復(fù)雜的淺海低頻聲場(chǎng)建模的研究則相對(duì)較少[7].

淺海海底結(jié)構(gòu)一般可視為由近孔隙介質(zhì)的沉積層和近彈性介質(zhì)的基巖層構(gòu)成，沉積層厚度一般在10 m左右[8-10].由于低頻水聲信號(hào)波長(zhǎng)長(zhǎng)且有較強(qiáng)的穿透力，其在淺海環(huán)境中傳播時(shí)能穿透沉積層到達(dá)基巖層，所以，在研究淺海低頻聲傳播問(wèn)題時(shí)將海底視為彈性介質(zhì)更為合理.此外，實(shí)際淺海海底受內(nèi)陸河流沖刷及大陸板塊擠壓等因素影響，其海底地形、海底底質(zhì)均存在較大差異.以海底地形為例，有的海域較平整，有的則成斜坡?tīng)?，有的甚至為陡峭的山狀?qiáng)切割地形.若仍采用傳統(tǒng)水平分層的淺海波導(dǎo)模型將無(wú)法準(zhǔn)確描述低頻聲信號(hào)在淺海中的傳播特性，更無(wú)法為淺海水聲工程應(yīng)用提供合理的理論參考.因此，在對(duì)淺海低頻聲傳播問(wèn)題的研究中，亟待開(kāi)展在淺海水平變換波導(dǎo)下的討論.

任何聲傳播問(wèn)題的研究基礎(chǔ)是合理的聲場(chǎng)計(jì)算方法，經(jīng)過(guò)近五、六十年的發(fā)展，目前在水聲傳播研究領(lǐng)域已建立起多種聲場(chǎng)計(jì)算方法，如簡(jiǎn)正波方法、射線方法、拋物方程方法、快速場(chǎng)方法以及一系列衍生計(jì)算方法[11].但上述方法都是在對(duì)波動(dòng)方程和海洋環(huán)境作出一些假設(shè)和近似的基礎(chǔ)上推導(dǎo)得到，其通用性均受到了一定的限制，尤其是對(duì)淺海水平變化波導(dǎo)下低頻聲傳播問(wèn)題的計(jì)算[12].隨著對(duì)淺海水平變換波導(dǎo)下低頻聲傳播問(wèn)題的日益關(guān)注，積極建立通用性更強(qiáng)的淺海聲場(chǎng)計(jì)算方法日益重要.有限元法(FEM)[13]把環(huán)境劃分為離散的單元進(jìn)行求解，因而更適用于描述聲場(chǎng)中的各類變化，但由于其計(jì)算量大，以往除用于提供基準(zhǔn)解外，一般在海洋聲傳播問(wèn)題中很少使用.近年來(lái)，計(jì)算機(jī)硬件水平的快速提高為FEM在海洋聲傳播問(wèn)題中的應(yīng)用提供了有效的支持.

基于上述研究現(xiàn)狀，針對(duì)淺海水平波導(dǎo)下低頻聲信號(hào)的傳播特性的研究需求，本文采用COMSOL Multiphysics(以下簡(jiǎn)稱COMSOL)[14]有限元軟件平臺(tái)建立了淺海低頻聲傳播的有限元計(jì)算模型.針對(duì)當(dāng)前淺海聲場(chǎng)建模研究中多數(shù)只關(guān)注海水層中聲場(chǎng)的計(jì)算建模和傳播規(guī)律且通常假設(shè)海底為液態(tài)介質(zhì)的研究現(xiàn)狀，以聲能流為研究對(duì)象[15]，同步討論低頻聲波在淺海水體/彈性海底中的傳播變化規(guī)律，以期為低頻聲信號(hào)在淺海中的探測(cè)應(yīng)用提供理論基礎(chǔ).

1 理論概述

1.1 聲場(chǎng)的有限元形式

針對(duì)本文所關(guān)注的淺海環(huán)境下的低頻聲傳播問(wèn)題，在三維柱坐標(biāo)系下建立了符合淺海環(huán)境特點(diǎn)的波導(dǎo)模型.如圖1(a)所示，模型中深度為zs的簡(jiǎn)諧點(diǎn)聲源位于柱坐標(biāo)對(duì)稱軸上，海水層與海底層分別近似為均勻各向同性的流體介質(zhì)和彈性體介質(zhì)；ρw、cw分別為流體層中密度和聲速；ρb、cp及cs分別為海底層中的密度、縱波聲速和橫波聲速；αp、αs分別為海底層中的縱波聲速衰減和橫波聲速衰減；海水層深度H(r)隨水平距離r變化；Sw為海水層，在模型中表示為均勻流體；Sb為海底層，在模型中表示為均勻彈性介質(zhì).由于柱坐標(biāo)系的軸對(duì)稱性，可以將三維問(wèn)題轉(zhuǎn)化為二維(r,z)平面上求解[14].

圖1 淺海聲場(chǎng)模型Fig.1 Acoustic field model in shallow water

有限元法求解過(guò)程中的定義域如圖1(b)所示.圖中：L1邊界為柱坐標(biāo)中心對(duì)稱軸邊界；L2為海水層與淺海海表交互作用邊界；L3為海水層與彈性海底層的交互作用邊界；模型四周為完美匹配層(Perfectly Matched Layer，PML)[13].

對(duì)Helmholtz方程進(jìn)行權(quán)重積分并結(jié)合高斯理論，Sw中的聲學(xué)有限元方程以及Sb中的有限元振動(dòng)方程可表示為[13-14]

(Kw+iωCw-ω2Mw)pi=[Fwi]

(1)

(Kb+iωCb-ω2Mb)si=[Fbi]

(2)

式中：Kw、Cw及Mw分別為剛度矩陣、阻尼矩陣及質(zhì)量矩陣；ω為模型中簡(jiǎn)諧點(diǎn)聲源發(fā)出聲波的角頻率；pi為海水層節(jié)點(diǎn)i處的聲壓；Fwi為聲學(xué)激勵(lì)；Kb、Cb及Mb分別為彈性海底層中上沒(méi)有受到邊界約束部分的剛度矩陣、阻尼矩陣及質(zhì)量矩陣；si為海底層節(jié)點(diǎn)i處的位移；Fbi為結(jié)構(gòu)體上的激勵(lì)載荷.

在計(jì)算時(shí)，設(shè)定L1為Dirichlet邊界條件，即邊界L1上聲壓p滿足恒等于0的物理假設(shè)，其表達(dá)式為

p=0,在L1上

(3)

在海水層與彈性海底層的相互耦合邊界L4上建立聲-彈耦合方程，即滿足法向位移、法向應(yīng)力連續(xù)，切向應(yīng)力為0，其表達(dá)式為

(4)

(5)

(6)

為模擬聲信號(hào)在無(wú)限大的海洋環(huán)境中的傳播，在模型四周增加完美匹配層(Perfectly Matched Layer，PML)[13]，通過(guò)在有限元方程中增加吸收系數(shù)轉(zhuǎn)化為PML方程.

通過(guò)耦合上述有限元法建立的物理場(chǎng)方程、流-彈耦合方程和PML方程，便可實(shí)現(xiàn)對(duì)淺海環(huán)境下海水層中聲壓p和海底層中位移場(chǎng)s的計(jì)算，通過(guò)p、s與各層其他物理場(chǎng)間的對(duì)應(yīng)關(guān)系，可進(jìn)一步得到模型中如振速場(chǎng)v、應(yīng)力場(chǎng)T等其它物理場(chǎng)數(shù)值[15].在本文研究中，上述計(jì)算過(guò)程將在有限元軟件COMSOL平臺(tái)上實(shí)現(xiàn).

1.2 波導(dǎo)中聲能流的計(jì)算

在傳統(tǒng)標(biāo)量聲場(chǎng)的建模中，對(duì)聲場(chǎng)能量的研究多是基于傳播過(guò)程中聲壓的變化特性展開(kāi).相比聲壓，聲能流更適合揭示聲波能量傳播中的一般性“流動(dòng)”特征且更適合用于對(duì)水體/海底中聲能量特性的統(tǒng)一討論.因此本文研究中，也將借助聲能流密度這一矢量對(duì)低頻聲信號(hào)在淺海環(huán)境下的傳播特性展開(kāi)討論.在本文所定義的(r,z)平面內(nèi)，頻域上海水層與彈性海底層中的聲能流Iw、Ib定義為

(7)

(8)

(9)

(10)

下文以In為研究對(duì)象，對(duì)淺海水平變化波導(dǎo)下的低頻聲信號(hào)傳播特性進(jìn)行研究.

2 有限元法聲場(chǎng)計(jì)算結(jié)果驗(yàn)證

為驗(yàn)證FEM對(duì)海洋聲傳播問(wèn)題計(jì)算結(jié)果的正確性，本節(jié)首先使用FEM仿真計(jì)算了水平[15]及楔形上坡[11]兩類海底地形下的聲壓傳播特性，并與已有聲場(chǎng)計(jì)算程序所得結(jié)果進(jìn)行對(duì)比.所設(shè)定兩類類海洋環(huán)境模型與參數(shù)示意圖如圖2所示.

兩類模型中，均設(shè)定水體中cw=1 500 m/s，ρw=1 000 kg/m3.圖2(a)為一類水平分層海底模型，水深在r=0至r0=4 km (r0為模型最遠(yuǎn)距離)處恒為100 m.圖2(b)上坡海底海洋模型設(shè)定為美國(guó)聲學(xué)學(xué)會(huì)(Acousticcal Society of America，ASA)提出的楔形海域檢驗(yàn)?zāi)Ｐ?，即在水深在r=0至r0=4 km處由200 m線性減小為0，上坡傾斜角度約為α1=2.68°.兩類模型中均考慮海底為彈性介質(zhì)，圖2(a)中ρb=1 500 kg/m3，cp=2 000 m/s，cs=1 000 m/s，αpλ=αsλ=0.1 dB (λ為波長(zhǎng))；圖2(b)中ρb=1 500 kg/m3，cp=1 700 m/s，cs=800 m/s，αpλ=λαs=0.5 dB.圖2(a)中zs=20 m，頻率100 Hz.圖2(b)中zs=100 m，頻率25 Hz.

圖3所示為分別采用FEM以及現(xiàn)有聲場(chǎng)計(jì)算方法模擬得到的圖2中兩類海洋環(huán)境下的單接收點(diǎn)(z=30 m)聲壓傳播損失(Transmission Loss，TL)對(duì)比圖[13]，TL計(jì)算公式為

(11)

式中：pref│r=1 m為距離聲源1 m處的參考聲壓[1].

圖3中虛線為FEM計(jì)算結(jié)果，實(shí)線對(duì)應(yīng)現(xiàn)有聲場(chǎng)計(jì)算方法計(jì)算結(jié)果.計(jì)算中，根據(jù)環(huán)境模型與現(xiàn)有聲場(chǎng)計(jì)算方法適用性，分別選擇快速場(chǎng)程序(Scotter)和拋物方程程序(RAMs)進(jìn)行求解[16].

可以看出，F(xiàn)EM與兩種聲場(chǎng)計(jì)算程序仿真計(jì)算得到的傳播損失曲線完全一致，驗(yàn)證了本文所研究FEM方法在計(jì)算水聲場(chǎng)結(jié)果的準(zhǔn)確性.從仿真結(jié)果也可以看到，在聲源頻率及聲源深度一致的前提下，海底地形、海底聲學(xué)參數(shù)的變化將對(duì)波導(dǎo)中聲能量的傳播產(chǎn)生重要影響，而海底地形在水平方向上起伏、變化又廣泛存在于實(shí)際淺海環(huán)境中.因此，深入研究海底地形變化對(duì)低頻聲信號(hào)的影響規(guī)律，才能準(zhǔn)確將低頻聲傳播特性用于水聲工程實(shí)踐應(yīng)用中.

針對(duì)上述目標(biāo)，在本節(jié)研究基礎(chǔ)上，下一節(jié)將通過(guò)數(shù)值模擬進(jìn)一步討論海底地形變化對(duì)低頻聲傳播的影響規(guī)律，以期為低頻聲信號(hào)的應(yīng)用提供理論基礎(chǔ).考慮FEM在使用中對(duì)波動(dòng)方程不做任何假設(shè)和近似，可以處理各種幾何形狀下的物理場(chǎng)問(wèn)題，因此更適用于對(duì)淺海實(shí)際水平復(fù)雜變化波導(dǎo)環(huán)境下的計(jì)算，下文研究中在COMSOL平臺(tái)上利用FEM方法開(kāi)展仿真研究.

圖2 兩類海底地形下海洋環(huán)境模型示意圖Fig.2 Schematic diagram for marine environment models of two types of seabeds

圖3 接收深度z=30 m時(shí)兩類海洋環(huán)境下TL曲線對(duì)比圖Fig.3 TL versus r at z=30 m in two marine environments

3 海底地形變化對(duì)聲信號(hào)傳播特性的影響

3.1 楔形海底

圖4～6分別為當(dāng)海底聲學(xué)參數(shù)不變，海底地形分別為水平、上坡及下坡時(shí)水下20 m處100 Hz點(diǎn)聲源激發(fā)起聲能流在波導(dǎo)中的傳播特征，聲能流強(qiáng)度用聲強(qiáng)級(jí)SIL表示，其計(jì)算式為

Iref=6.76×10-19(W/m2)

(12)

式中：Iref為參考聲強(qiáng)[1].

以下圖中各點(diǎn)聲能流的大小和傳播方向分別用箭頭的長(zhǎng)度和方向來(lái)表示.除地形變化外，圖4～6其余仿真條件均與圖2(a)中一致，ΔH1為水平方向4 km處海面與海底的垂直距離.

圖4 水平海底仿真環(huán)境下聲能流傳播圖(ΔH1=100 m)Fig.4 Distribution of acoustic energy flux in simulation of horizontal seabed (ΔH1=100 m)

圖5 楔形上坡海底仿真環(huán)境下聲能流傳播圖Fig.5 Distribution of acoustic energy flux in simulation of up-sloping wedge-shaped seabed

圖6 楔形下坡海底仿真環(huán)境下聲能流傳播圖Fig.6 Distribution of acoustic energy flux in simulation of down-sloping wedge-shaped seabed

在圖2(a)仿真條件下，100 Hz聲波在其水體將激發(fā)起共9階波導(dǎo)簡(jiǎn)正波，其對(duì)應(yīng)在水平海底界面上的掠射角θ∈[4°,40°][1].對(duì)比分析圖4和5的仿真結(jié)果可以看到，相比水平海底，由于受α1的影響，水體層中i′次反射后聲波與海底接觸時(shí)的掠射角將增加(2i′-1)α1，原本小掠射角對(duì)應(yīng)的低階簡(jiǎn)正波傳播特性不斷向大掠射角、易衰減的高階簡(jiǎn)正波耦合，聲能量向海底的“泄漏”效應(yīng)顯著增強(qiáng)且起伏跨度縮短.α1越大，低階簡(jiǎn)正波向高階簡(jiǎn)正波的耦合越多，高階聲能量向海底“泄漏”效應(yīng)越強(qiáng)，對(duì)應(yīng)水體層中總聲能量相應(yīng)的衰減也越快，傳輸至遠(yuǎn)場(chǎng)后所保留的簡(jiǎn)正波階數(shù)也越少，干涉特性越簡(jiǎn)單.在圖5仿真條件下，隨著ΔH1從100 m變化至0，水層中100 Hz聲能量傳輸至r=4 km處衰減為20 dB左右.

對(duì)比分析圖4和6可以看到，由于有α2的存在，聲源位于淺水處，隨著水平傳播距離的增大水深不斷增加，水體層中i′次反射后聲波與海底接觸時(shí)的掠射角將減小(2i′-1)α2，與楔形上坡海底的影響規(guī)律相反，原本大掠射角對(duì)應(yīng)的高階簡(jiǎn)正波傳播特性不斷向小掠射角的低階簡(jiǎn)正波耦合，聲能量向海底的“泄漏”效應(yīng)明顯減弱且起伏跨度增大.α2越大，高階簡(jiǎn)正波向低階簡(jiǎn)正波的耦合越多，聲能量向遠(yuǎn)處傳播時(shí)的衰減越慢，聲場(chǎng)的干涉特性也越簡(jiǎn)單，即楔形下坡海底的存在對(duì)聲場(chǎng)形成了增強(qiáng)效應(yīng).

3.2 小型海底山，小型海溝

圖7 海底山、海溝地形下海洋環(huán)境模型示意圖Fig.7 Schematic diagram of marine environment model of seamount and trench

圖8和9分別為當(dāng)海底聲學(xué)參數(shù)不變，海底存在小型海底山與小型海溝時(shí)100 Hz點(diǎn)聲源激發(fā)起聲能流在波導(dǎo)中的傳播特征，仿真中設(shè)定海底山與海溝中心位置位于水平方向2 km處.除地形變化外，圖8和9中其余仿真條件均與圖2(a)中一致.

圖8中設(shè)定在海底水平距離1.5～2.5 km處存在一小型海底山.

圖8 海底山仿真環(huán)境下聲能流傳播圖Fig.8 Distribution of acoustic energy flux in simulation of seamount

對(duì)比分析圖4和8的仿真結(jié)果可以看到，相比圖4水平海底環(huán)境下的聲能量傳播特性，由于小型海底山的“阻擋”，整個(gè)仿真?zhèn)鞑ミ^(guò)程中聲能量的衰減增加，且海底山越高，聲能量的衰減越大.從圖8的對(duì)比中可以看出，當(dāng)海底存在海底山影響時(shí)，海底山前、后聲能量的傳播規(guī)律不相同.由于仿真條件下的海底山地形具有底邊短、變化快的起伏特點(diǎn)，低頻聲波受其影響主要集中于海山前后，在點(diǎn)聲源激發(fā)起的聲能量到達(dá)海底山前受海底山地形的影響并不明顯.隨著聲波不斷臨近海底山，受海底山前沿的反射作用，使部分原沿水平方向傳播的聲能量向海底山上方水層反射，聲能量在海底山上方升高，即對(duì)應(yīng)z<ΔH2處的聲能量受海底山的“阻擋”作用并不顯著.而同時(shí)在此過(guò)程中，部分聲能量由海底山上坡處透射進(jìn)入海底，該部分聲能量主要對(duì)應(yīng)掠射角較大，易衰減的中、高階簡(jiǎn)正波，因此當(dāng)聲能量在“翻越”海底山的過(guò)程中其包含的中、高階簡(jiǎn)正波能量更易快速衰減，在“翻越”海底山后只余少量低階簡(jiǎn)正波存在，因此海底山后聲能量起伏穩(wěn)定且衰減較慢；而當(dāng)對(duì)接收點(diǎn)深度z>ΔH2時(shí)，海底山對(duì)聲能量傳播產(chǎn)生的明顯的“阻擋”作用，聲能量衰減明顯.

與圖8仿真結(jié)果相對(duì)應(yīng)，圖9所示為當(dāng)海底存在海溝時(shí)的聲波傳播特性的仿真圖.仿真中設(shè)定在海底r=1.5～2.5 km處存在一海溝.

圖9 海溝仿真環(huán)境下聲能流傳播圖Fig.9 Distribution of acoustic energy flux in simulation of trenches on seabed

4 結(jié)論

本文基于有限元法，結(jié)合具體仿真算例討論了淺海水平變化波導(dǎo)下聲能量的傳播變化特性及機(jī)理，以期為淺海水聲工程設(shè)備的研發(fā)，測(cè)試及使用等提供一些理論參考.對(duì)比淺海水平分層波導(dǎo)下的聲能量傳播規(guī)律，結(jié)果表明：

(1)利用有限元法可以準(zhǔn)確實(shí)現(xiàn)對(duì)各類復(fù)雜海洋環(huán)境下聲場(chǎng)分布的求解.

(2)無(wú)論是楔形上坡或是下坡海底，二者均會(huì)對(duì)淺海中聲能量的傳播產(chǎn)生重要影響.相比水平海底，上坡海底傾角α1越大，高階聲能量向海底“泄漏”效應(yīng)越強(qiáng)，聲能量向遠(yuǎn)處傳播時(shí)的衰減越快；而楔形下坡海底地形的影響則正好相反，傾斜角度α2越大，聲能量向遠(yuǎn)處傳播時(shí)的衰減越慢.正波的轉(zhuǎn)移越多，聲能量向遠(yuǎn)處傳播時(shí)的衰減越慢，即楔形下坡海底的存在對(duì)聲場(chǎng)形成了增強(qiáng)效應(yīng).

(3)當(dāng)海底存在小型海底山影響時(shí)，由于海底山上坡面的反射作用，點(diǎn)聲源激發(fā)起的聲能量在“翻越”海底山的過(guò)程中，將在海底山上方產(chǎn)生提升，海底山后聲能量起伏穩(wěn)定且衰減較慢.

(4)與海底山的影響類似，在點(diǎn)聲源激發(fā)起的聲能量到達(dá)海溝前其受地形的影響并不明顯，當(dāng)海溝深度較淺，聲能量在“跨越”海溝過(guò)程中將產(chǎn)生部分聲能量的額外衰減.隨著海溝深度的增加，當(dāng)存在夾角大于波導(dǎo)中部分低階簡(jiǎn)正波的掠射角時(shí)，聲能量在海溝上方水層產(chǎn)生提升，隨著海溝深度的繼續(xù)增加，“跨越”后聲能量的起伏特性與水平海底基本一致.

致謝在此感謝哈爾濱工程大學(xué)水聲工程學(xué)院樸勝春教授課題組為本論文的撰寫(xiě)提供有限元軟件支持.