傾斜彈性海底條件下淺海聲場的簡正波相干耦合特性分析*

張士釗 樸勝春1)2)?

1) (哈爾濱工程大學,水聲技術(shù)重點實驗室,哈爾濱 150001)

2) (哈爾濱工程大學,海洋信息獲取與安全工信部重點實驗室,哈爾濱 150001)

3) (哈爾濱工程大學水聲工程學院,哈爾濱 150001)

建立了一種適用于復本征值情況下的彈性耦合簡正波聲場模型,給出了包含泄漏模態(tài)時彈性簡正波的歸一化和耦合系數(shù)表達式,耦合系數(shù)滿足聲能流守恒.利用該模型分析了傾斜彈性海底條件下聲場的簡正波相干耦合特性,發(fā)現(xiàn)考慮泄漏模態(tài)時,簡正波耦合不僅會導致簡正波幅度變化,而且也會帶來附加相移.仿真計算表明:其一,在傾斜彈性海底條件下,考慮泄漏模態(tài)耦合附加相移影響后得到的聲傳播損失更接近有限元商用軟件的計算結(jié)果;其二,耦合大大提升了界面波的幅度.此外,本文還分析了海洋環(huán)境參數(shù)變化對聲傳播損失的影響.

1 引言

在彈性海底情況下,聲波在海洋中傳播時,海底會激發(fā)起縱波和橫波.而在海洋環(huán)境隨著距離發(fā)生變化的情況下,海洋中傳播的各階簡正波之間會發(fā)生簡正波耦合效應,帶來簡正波之間的能量轉(zhuǎn)移.

Pekeris[1]首次提出簡正波模型,解決了液態(tài)海底條件下兩層水平分層介質(zhì)中的水中點源爆炸聲傳播問題.Press 和Ewing[2]將Pekeris 提出的簡正波模型擴展到包含彈性海底的情況,解決了彈性海底條件下的水平分層介質(zhì)中的爆炸聲傳播問題.Porter 和Reiss[3]采用有限差分方法,數(shù)值計算彈性海底水平分層介質(zhì)中的簡正波聲傳播.Ainslie[4]計算了多層彈性海底條件下水平分層介質(zhì)中的縱波和橫波反射和傳播系數(shù).Schneiderwind 等[5]研究了水平分層彈性Pekeris 波導條件下的聲場簡正波解,并通過水池實驗加以驗證.

實際海洋環(huán)境隨距離變化,導致各階簡正波之間的耦合效應.Pierce[6]首次將簡正波聲場分析方法擴展到距離依賴海洋環(huán)境中,提出了簡正波耦合模型.Chwieroth 等[7]將耦合簡正波方法用于隨距離變化的拋物線聲速模型,分析了該環(huán)境下的簡正波能量轉(zhuǎn)移.Williams[8]分析了隨距離線性變化海深條件下,上坡和下坡時的前兩階簡正波之間的能量耦合.Rutherford 和Hawker[9]分析了傾斜海底環(huán)境的邊界條件,用修正邊界條件對簡正波耦合理論做了一階能量守恒修正.McDaniel[10]分析了常量坡度傾斜液態(tài)海底條件下的簡正波能量交換,發(fā)現(xiàn)了靠近截止頻率的簡正波對連續(xù)譜簡正波的強烈耦合效應.Evans[11]采用距離步進簡正波耦合模型,考慮了水深隨距離變化帶來的反向散射,以及有限階傳播簡正波之間的耦合.Porter 等[12]采用改進邊界條件,解決了傾斜海底階梯近似中的單向拋物方程模型能量守恒問題.Abawi 和Kuperman[13]采用耦合簡正波拋物方程模型,計算距離依賴海洋環(huán)境下的聲場,處理簡正波之間以及不同方位角之間的能量耦合.Knobles 等[14]采用雙向積分耦合簡正波方法,分析大陸架海域的聲傳播,既包括離散譜,也包括連續(xù)譜,并且用泄漏模態(tài)近似表示連續(xù)譜.莫亞梟等[15]提出了一種滿足能量守恒的雙向簡正波耦合方法,實現(xiàn)了非水平分層波導中的雙向場的準確計算.

彈性海底情況下,簡正波耦合問題變得更加復雜.Kennett[16]分析了水平變化海洋環(huán)境下的彈性簡正波耦合,包括了前向和后向傳播的簡正波.Odom[17]研究了隨距離變化海深彈性海底條件下,簡正波之間的耦合,包括與連續(xù)譜的耦合.Hall[18]通過測量垂直于海底的深度,解決了傾斜海底時簡正波之間耦合系數(shù)的計算錯誤;發(fā)現(xiàn),淺水傾斜海底環(huán)境下,負聲速梯度導致了簡正波耦合系數(shù)增大.Maupin[19]采用局地簡正波耦合方法處理距離依賴的彈性海底海洋環(huán)境,分析了二維聲場中的Love 波和Rayleigh 波之間的簡正波能量轉(zhuǎn)移.Collins[20]采用能量守恒拋物方程方法,解決涉及傾斜海底界面波傳播的問題,以及涉及簡正波截止和耦合至海底剪切波束的問題.Collins 和Siegmann[21]對彈性拋物方程進行聲能流密度守恒修正,處理距離依賴海深時的聲場.Jeroen[22]分析了海底地形變化條件下,界面波傳播中的簡正波耦合,而Abawi 和Porter[23]將該成果應用于傾斜彈性海底條件下的耦合簡正波傳播,得到了該條件下的縱波和橫波聲壓.Odom 等[24]將彈性海底中的傾斜邊界條件轉(zhuǎn)換為與距離導數(shù)無關的項,從而簡化了彈性海底時的簡正波耦合.Abawi 和Porter[25]采用虛源法計算楔形彈性海底條件下的聲場,得到了和彈性拋物方程方法一致的結(jié)果.Xie 等[26]采用拉格朗日譜元法,分析半無限空間中的流體-固體耦合問題,處理距離有關海洋環(huán)境中的聲傳播問題.

在海洋環(huán)境隨著距離發(fā)生變化的情況下,簡正波耦合的計算過程中會受到一些因素的影響.1 個典型的影響是,特征方程中,垂直波數(shù)的計算往往會選取Pekeris 割線法,而特定環(huán)境下,算出的本征值過于接近割線位置,導致本征值出現(xiàn)較大的誤差,甚至搜根失敗,嚴重影響最終的簡正波歸一化和耦合結(jié)果.為了消除Pekeris 割線對聲場計算的影響,一些學者采用了多種方法來處理.駱文于等[27]采用穩(wěn)定的深度格林函數(shù),結(jié)合波數(shù)積分方法,得到了簡正波水平波數(shù)接近海底波數(shù)時的準確聲壓.而莫亞梟等[28]采用復等效深度法,提高了波導簡正波本征值位于割線枝點附近情況下聲傳播損失的計算精度.Chapman 等[29]采用等效深度法處理彈性Pekeris 波導中的聲傳播問題,得到了精確的簡正波相速度.McCollom 和Collis[30]采用格林函數(shù)處理彈性海底Pekeris 波導中的簡正波求解問題.Zhang 和Tindle[31]將采用復等效深度法,處理淺海彈性海底時的簡正波求解問題,包含了“泄漏”模態(tài).Stickler[32]采用割線積分方法處理“泄漏”簡正波,發(fā)現(xiàn)了割線積分對淺水聲場的顯著貢獻.Westwood 和Koch[33]采用弱正梯度聲速半無限空間近似方法,在液態(tài)海底環(huán)境下,消除了Pekeris割線附近的簡正波計算誤差;Westwood 等[34]還將此方法應用于彈性海底情況,結(jié)合等效反射系數(shù)和等效深度方法修正本征值.

前人的工作,多從位移和應力張量的角度,將耦合簡正波分為局地簡正波,和與距離有關的簡正波幅度,并采用耦合差分方程來計算簡正波幅度,從而得到彈性海底時,海洋環(huán)境變化情況下的聲場.然而,局地簡正波階梯近似的方法,多見于液態(tài)海底的情況;彈性海底時,這方面的分析較少,而且多局限于“受限”的簡正波,很少有對“泄漏”簡正波的分析.而本文將簡正波耦合擴展到包含“泄漏”簡正波的情況,建立了復本征值彈性耦合簡正波模型.在距離依賴彈性海底的海洋環(huán)境中,該模型提高了對聲傳播損失的預報精度.此外,還將能量守恒近似擴展以包含復數(shù)本征值(泄漏簡正波),來保證結(jié)果的準確性.本文還分析了傾斜彈性海底條件下的簡正波耦合效應,發(fā)現(xiàn)了泄漏模態(tài)對簡正波耦合系數(shù)的影響,以及由此導致的簡正波聲壓相位和幅度的復雜變化.

2 復本征值彈性耦合簡正波模型

在二維直角坐標系下,考慮時間因子為exp(iωt)的點聲源,不考慮外力作用,也不考慮y方向的聲場,在彈性海底水平變化波導中激發(fā)起聲場,采用階梯近似來處理水平變化彈性海底邊界,如圖1所示.

圖1 水平變化彈性海底的分段階梯近似Fig.1.Stepwise approximation of a horizontally varying elastic seabed.

2.1 彈性波導中聲場的運動方程

在彈性海底水平變化波導中激發(fā)起的聲場,其應力t和質(zhì)點位移w滿足的運動方程如下:

這里ρ為介質(zhì)密度,應力矢量t定義如下:

其中Cij為彈性模量有關系數(shù),詳見附錄A.假設聲源為線源,只考慮xOz平面,y方向保持不變,

將位移w做時域傅里葉變換,得到

其中,x表示接收器到聲源的距離,z表示接收深度.運動方程:

其中頻域位移w=(wx,wz),wx和wz分別為水平位移和垂直位移分量.t1=(t11,t13),t3=(t31,t33),ti=Ci1?w/?x+Ci3?w/?z.

下文中,定義t=t1=(τxx,τxz),位移應力矢量u=(w,t)T.從而得到

其中差分算子:

矩陣Qgh定義如下:

邊界條件要求界面處的位移和應力連續(xù),對于形如z=h(x)的第m層界面,法向為n,界面斜率dzm/dx=?hm/?x=h˙m,界面傾角θm=arctan(h˙m),法向矢量可以寫成n=sinθmi–cosθmk,第m層界面上的應力T=tjnj連續(xù)可以寫成:

方括號表示第m層界面處變量的值的跳變,從底部(z＞hm)到頂部(z＜hm),下同.根據(jù)(8)式,得到邊界條件

由此可將運動方程(5)和邊界條件(9)合寫成:

相應的邊界條件為

需要說明的是,(10)式中的δ(z-hm(x)) 項作用于u矢量這個整體,而當考慮下文中的局地簡正波un(x0,z)exp{-ikn(x0)x}時,δ(z-hm(x)) 項作為邊界條件存在,用于推導正交歸一化,不直接出現(xiàn)在局地簡正波un(x0,z) 表達式中.

2.2 復本征值簡正波歸一化

二維距離依賴介質(zhì)中的聲傳播問題可以通過距離無關介質(zhì)中的本征函數(shù)來解決.運動方程(10)的各向同性局部簡正波解可以寫成:

其中un滿足:

在(13)式中,

矩陣算子A中各個參數(shù)含義見附錄A.

我們發(fā)現(xiàn),在本征值為復數(shù)的情況下,根據(jù)特征方程[5]

求得本征值的方法會有較大的誤差,其中kp和ks為海底縱波和橫波波數(shù).為了解決這個問題,對特征方程進行黎曼面上的分支處理:

可以得到滿足(13)式的復數(shù)kn(x).而如果不進行黎曼面上的分支處理,將無法準確求得方程的復數(shù)根.

2.3 復本征值簡正波正交性

在復本征值條件下,以往的簡正波正交性條件[19]不再成立,需要給出新的適用于復數(shù)本征值kn(x)的正交性推導.定義局地本征函數(shù)um和un的標量積如下:

經(jīng)過一些推導(見附錄B),得到適用于復本征值kn(x)的正交性表達式:

對un實現(xiàn)歸一化,適用于所有的復數(shù)kn(x).而參考文獻[19]的(A2)式(被參考文獻[24]用于推導簡正波歸一化)中給出的簡正波正交性表達式i(km-kn)〈um,un〉=0,僅適用于實數(shù)kn(x),和本文推導出的適用于復數(shù)kn(x)的簡正波正交歸一化表達式(15)并不相同.

根據(jù)(14)—(16)式進行簡正波歸一化,采用了簡正波幅度Am的共軛相乘,歸一化除法因子含有,歸一化前的簡正波聲壓幅度在歸一化后變?yōu)?相位得到保留.在經(jīng)典的聲壓計算方法中,通常采用朗斯基歸一化,簡正波聲壓幅度在歸一化后將變?yōu)?.為了和經(jīng)典的聲場計算方法所的結(jié)果進行對比,我們將經(jīng)過簡正波正交歸一化的聲壓除以,以便使簡正波聲壓幅度變?yōu)?,消去幅度的相位,得到修正后的歸一化所得聲壓,且和未修正的聲壓進行了對比(圖3).

2.4 包含泄漏模態(tài)的耦合簡正波模型

基于上述歸一化,只考慮沿著x軸正方向傳播的簡正波,有

其中幅度am滿足:

實際數(shù)值計算中采用離散的差分形式,將(18)式寫成:

其中Bmn(xj+1,xj)為考慮泄漏模態(tài)的耦合系數(shù),表達式見附錄C 中的(C1)式.

由此得到海洋環(huán)境隨距離變化時,彈性海底條件下的簡正波耦合系數(shù)表達式,用以計算該條件下的水中聲壓.需要特別說明的是,(C1)式的結(jié)果,和參考文獻[24]中的耦合系數(shù)表達式(A1)有所不同.這個不同的結(jié)果,使得kn(x)為復數(shù)時,仍滿足參考文獻[24]中給出的能量守恒表達式(27),而文獻[24]中的耦合系數(shù)表達式(A1),僅適用于實數(shù)kn(x)的情況,而對復數(shù)kn(x)失效.文獻[24]中的結(jié)果不能用于復數(shù)本征值,也就是相應的泄漏模態(tài),在利用耦合簡正波方法計算聲場的時候,僅考慮實數(shù)本征值,會帶來較大的聲壓計算誤差.而本文提出的改進方法,則能夠?qū)⒃擇詈虾喺ǚ椒ㄍ茝V到包含復數(shù)本征值的情況,得到的聲壓更加接近較為成熟的彈性聲場計算方法,如RAMS 和有限元商用軟件COMSOL.這些結(jié)果將會在仿真分析部分予以展示.

以上推導了包含泄漏模態(tài)時的彈性海底條件下的簡正波歸一化,修正了能量守恒,推導了耦合簡正波表達式,得到了距離依賴彈性海底海洋環(huán)境下的聲場模型.接下來,將對這些模型進行仿真,與經(jīng)典聲場計算方法進行對比檢驗.

3 仿真分析

3.1 模型對標準檢驗問題的計算

3.1.1 彈性海底和液態(tài)海底的水中聲壓對比

分別用本文提出的簡正波耦合方法(圖2 中標注為“耦合矩陣”,下同.下文簡稱為“本文方法”)、RAMS、有限元商用軟件COMSOL 計算彈性海底時的水中聲場(本文方法和RAMS 計算聲場的水平步長均為1 m,下同),用Couple07,RAM 計算液態(tài)海底時的水中聲場.zs和zr分別表示聲源和接收深度,zs=zr=30 m,聲波頻率f=50 Hz,海洋環(huán)境和所得聲壓如圖2 所示.

從圖2(b)可見,Couple07 和RAM 程序計算液態(tài)海底時的水中聲壓,結(jié)果基本一致,說明這兩種方法所得液態(tài)海底時的水中聲壓正確.從圖2(d)可見,本文方法和RAMS 以及COMSOL,計算彈性海底時的水中聲壓基本一致,證明了本文方法計算彈性海底時的水中聲壓的正確性.此外,對比圖2(b)和圖2(d)可見,彈性和液態(tài)海底時的水中聲壓有較大差異,體現(xiàn)了彈性海底對水中聲壓的顯著影響.后續(xù)計算中,如無特別說明,f均為50 Hz,zs和zr均為30 m,聲壓計算距離均為3400 m,海洋環(huán)境均如圖2(c)所示.

3.1.2 彈性海底標準檢驗問題的計算結(jié)果

考慮如圖2(c)所示的彈性標準楔形海底條件下的海洋環(huán)境.分別采用本文方法、RAMS 以及COMSOL,計算水中聲壓,所得結(jié)果如圖3 所示.

由圖3(a)可知,未經(jīng)幅度相位修正的〈um,un〉歸一化所得聲壓,和COMSOL 所得聲壓基本一致,但是局部差異較大,特別是x=3100 m 附近出現(xiàn)了較大的不一致.而修正后的聲壓,如圖3 所示,和RAMS 以及COMSOL 所得聲壓基本一致,說明了幅度相位修正的必要性.

圖2(a) 液態(tài)海底海洋環(huán)境模型;(b) 液態(tài)海底時的水中聲壓;(c) 彈性海底海洋環(huán)境模型;(d) 彈性海底時的水中聲壓Fig.2.(a) Marine environment model with liquid seabed;(b) sound pressure in water with liquid seabed;(c) marine environment model with elastic seabed;(d) sound pressure in water with elastic seabed.

圖3 是否考慮幅度相位處理的 〈um,un〉 歸一化時本文方法所得聲壓,及其和RAMS 以及COMSOL 所得聲壓對比 (a) f=50 Hz,不考慮;(b) f=50 Hz,考慮;(c) f=25 Hz 時,不考慮;(d) f=25 Hz,考慮;(e) f=25 Hz,cs=2000 m/s,不考慮;(f) f=25 Hz,cs=2000 m/s,考慮Fig.3.Sound pressure obtained by the method in this paper when (a) f=50 Hz,without,(b) f=50 Hz,with,(c) f=25 Hz,without,(d) f=25 Hz,with,(e) f=25 Hz,cs=2000 m/s,without,(f) f=25 Hz,cs=2000 m/s,with considering the normalization of amplitude and phase processing,compared with the sound pressure obtained by RAMS and COMSOL.

本文選擇f=50 Hz,是因為此時海深H=200 m 處算出多達14 階簡正波,單階簡正波的局部計算誤差對聲場的影響較小;而常規(guī)的標準楔形海底算例采用f=25 Hz.為進一步驗證COMSOL計算結(jié)果的準確性,計算25 Hz 時的水中聲壓,如圖3(c)和圖3(d)所示,可看出該條件下,RAMS和COMSOL 所得聲壓基本一致,從而驗證了COMSOL 所得聲壓的準確性,可以用COMSOL來檢驗本文方法聲壓計算結(jié)果的準確性.從圖3(c)可見,不考慮簡正波幅度的相位處理時,本文方法和COMSOL 所得聲壓在x＞ 2.2 km 時相差很大;而從圖3(d)可見,考慮簡正波幅度的相位處理后,本文方法所得聲壓與COMSOL 所得聲壓的接近程度在x＞ 2.2 km 時大幅提升,表明這樣的處理能夠得到更準確的聲壓.此外,從圖3(b)和圖3(d)可見,f=25 Hz 時,本文方法所得聲壓在一些位置處和COMSOL 計算結(jié)果相差較大,而50 Hz 時則與COMSOL 計算結(jié)果基本一致.原因是,本征值接近割線位置ks時計算不準確,導致局部聲壓計算錯誤.這樣的錯誤出現(xiàn)在越低階次的簡正波,對聲壓計算的影響越大,而在本文的計算中,25 Hz時,第2 階簡正波本征值在x=3271—3350 m 處靠近ks而計算出錯,而50 Hz 時,從第3 階簡正波才開始出現(xiàn)x＜ 3400 m 時的本征值計算錯誤,造成的聲壓計算誤差小得多.而當cs從1800 m/s 增加到2000 m/s 時,圖3(f)中的本文聲壓計算結(jié)果和COMSOL 更符合,其原因在于,此時ks變小,低階簡正波本征值靠近ks割線計算出錯的距離增大,對前3400 m 距離的聲壓計算帶來的誤差變小,由此導致了更為準確的聲壓計算結(jié)果.因此,本文方法在聲波頻率更高以及海底橫波聲速更大的情況下,水中聲壓計算結(jié)果更準確.下文所有計算中,均考慮簡正波幅度相位修正后的簡正波正交歸一化所得聲壓.

除了對比本文方法和RAMS 以及COMSOL的聲壓計算精度,還對比了三者的計算速度,結(jié)果如表1 所列.

從表1 可見,耦合矩陣法計算聲場耗時和RAMS 較為接近,均比COMSOL 耗時短得多.

表1 圖2(c)環(huán)境下聲場計算的平均耗時Table 1.Average time consumption for acoustic field computation under the environment shown in Fig.2(c).

3.2 泄漏模態(tài)對彈性耦合簡正波的影響分析

3.2.1 泄漏模態(tài)對聲壓傳播損失的影響分析

本文討論的泄漏模態(tài),指的是本征值km實部小于ks的模態(tài).考慮和不考慮泄漏模態(tài)時本文方法所得聲壓,及其和RAMS,COMSOL 的結(jié)果對比,如圖4 所示.

從圖4(a)可見,不考慮泄漏模態(tài)時,本文方法所得水中聲壓和RAMS,COMSOL 結(jié)果相差較大.而從圖4(b)可見,考慮泄漏模態(tài)時,本文方法所得聲壓和RAMS,COMSOL 結(jié)果較為接近.這說明,在計算標準楔形彈性海底時的水中聲壓時,必須考慮泄漏模態(tài).這個結(jié)論,也可以通過波數(shù)域聲壓來佐證.

圖4(a) 不考慮和(b)考慮泄漏模態(tài)時的耦合矩陣法所得聲壓,及其和RAMS,COMSOL 程序所得聲壓的對比;(c) RAMS 所得聲壓的FK 變換結(jié)果Fig.4.Sound pressure obtained by the coupling matrix method when the leaky mode is (a) not considered or (b) considered,compared with the sound pressure obtained by the RAMS and COMSOL program;(c) the FK transformation result of the sound pressure obtained by RAMS.

對RAMS 所得聲壓p(x,zr) (Δx=1 m)進行頻率-波數(shù)域濾波(frequency-wavenumber filtering,FK)變換,得到波數(shù)域聲場

如圖4(c)所示,在kx＜ks時聲波幅度依然較強,不可忽略,因而泄漏模態(tài)對彈性海底條件下的水中聲壓的貢獻不可忽略.

3.2.2 聲源深度對聲場的影響分析

改變聲源深度,得到的聲壓見圖5.從圖5 可見,考慮泄漏模態(tài),在zs=50,70 和100 m 時,本文方法和COMSOL 所得聲壓基本一致.說明在改變zs時,本文算法所得聲壓能夠和COMSOL 相符合.而不考慮泄漏模態(tài)時,3 種不同zs時,本文方法與COMSOL 所得聲壓均有較大差異,說明在這些zs條件下,泄漏模態(tài)能量較高,對聲壓的貢獻不可忽略.

圖5 zr=30 m,zs=(a) 50,(b) 70 和(c) 100 m 時的水中聲壓;m=(d) 1,(e) 5 和(f) 9 時,x=0 m 處的耦合系數(shù)BmnFig.5.Sound pressure in water when zr=30 m,zs=(a) 50,(b) 70,and (c) 100 m;the coupling coefficient Bmn at x=0 m when m=(d) 1,(e) 5 和(f) 9.

此外,觀察到,不考慮泄漏模態(tài)時,聲壓起伏較小,聲場干涉不太明顯.而當考慮泄漏模態(tài)時,聲場干涉起伏明顯增強.為了解釋這個現(xiàn)象,以x=0 m 處的耦合系數(shù)Bmn為例來說明:

圖6(a) zs=50,(b) 70 和(c) 100 m 時,不含波數(shù)差的耦合系數(shù)取實部和保持為復數(shù)時的聲壓對比Fig.6.Acoustic pressure when the real part of the coupling coefficient without wavenumber difference is taken when zs=(a) 50,(b) 70,(c) 100 m,compared with the sound pressure when the coupling coefficient is kept as a complex number.

為了分析局地簡正波帶來的耦合系數(shù)Bmn的附加相位,去掉簡正波水平波數(shù)差帶來的相位因子,只考慮剩下的部分.以x=0 m 處為例,共計算得到14 階簡正波,其中前8 階為受限模態(tài)(第一階為界面波),后6 階為泄漏模態(tài).分別以界面波,某一受限模態(tài)(第5 階)和某一泄漏模態(tài)(第9 階)為例,畫出這些階和其他階簡正波之間的耦合系數(shù),如圖5(d)—(f)所示.從圖5(d)—(f)可見,界面波(第1 階)和其他受限模態(tài)(第2—8 階)之間的耦合系數(shù)Bmn是純虛數(shù),而受限模態(tài)(第5 階)和其他受限模態(tài)(第2—4,6—8 階)之間的耦合系數(shù)Bmn是實數(shù),這樣的耦合只帶來各階簡正波間的能量轉(zhuǎn)移,不會引入附加的相位變化.而考慮泄漏模態(tài)時,泄漏模態(tài)之間(第9階和第10—14 階之間),以及泄漏模態(tài)(第9 階)和受限模態(tài)(第1—8 階)之間的耦合系數(shù)Bmn是復數(shù)(實部和虛部均不為零),不僅導致各階簡正波間的能量轉(zhuǎn)移,也帶來附加相位變化,導致泄漏模態(tài)之間,以及泄漏模態(tài)和受限模態(tài)之間的干涉,從而帶來了比只考慮受限模態(tài)時更加復雜的聲場干涉,聲場干涉起伏明顯增強.

為了驗證上述分析的正確性,將不含波數(shù)差的耦合系數(shù)Bmn取實部,繪制聲壓傳播損失,且和Bmn保持為復數(shù)時的聲傳播損失進行對比,如圖6所示.

從圖6 可見,3 個不同zs情況下,無論不含波數(shù)差的耦合系數(shù)Bmn取實部,還是保持復數(shù),本文方法所得聲壓結(jié)果都基本一致,僅在某些局部略有差別.這說明,復數(shù)耦合系數(shù)導致的相干耦合以及附加相移,對聲傳播損失的影響較小,并不是考慮泄漏模態(tài)后聲場干涉結(jié)構(gòu)變得復雜的主要原因.

為了分析泄漏模態(tài)導致聲場干涉結(jié)構(gòu)變得復雜的原因,繪制各階簡正波本征值實部隨著距離的變化,如圖7(c)所示.從圖7(c)可見,第9—14 階簡正波在x=0 m 處就是泄漏模態(tài),它們只對近場聲壓有貢獻,而對遠場聲壓的貢獻可以忽略.而第2—8 階簡正波在x=0 m 處是傳播模態(tài),傳播過程中變?yōu)樾孤┠B(tài),它們變?yōu)樾孤┠B(tài)的距離較遠,從該位置處才開始指數(shù)衰減,指數(shù)衰減距離并不長,因而它們對聲壓的貢獻不可忽略,與傳播模態(tài)的干涉作用不可忽略,帶來了比只考慮受限模態(tài)時更復雜的聲場干涉.

3.2.3 耦合對簡正波聲壓幅度am的幅度和相位的影響

根據(jù)(18)式,簡正波耦合只影響各階簡正波的幅度am,而am是1 個復數(shù),因而分別分析耦合對復數(shù)am的幅度和相位的影響.絕熱近似和考慮耦合時,各階簡正波am的模|am|以及相位隨距離的變化如圖7(a)和圖7(b)所示.

一般情況下,聲源遠離海底時,界面波(圖中標注為第一階,下同)的激發(fā)強度很弱,幾乎無法被探測到.而從圖7(a)可見,zs=50 m 時,簡正波耦合對界面波的幅度影響最為明顯,將界面波的聲壓幅度|a1|提高了110 dB,使得其幅度達到和第5階簡正波近似的水平,使界面波有可能被探測到.而第5 階簡正波的幅度|a5|受耦合的影響,僅次于界面波,最大起伏達到32.4 dB.其原因在于,在聲源位置處,zs=50 m 位于第5 階簡正波的波節(jié)處附近,因而第5 階簡正波激發(fā)強度較小,容易受到其他階簡正波的耦合影響而帶來幅度的劇烈變化.而其他階簡正波的幅度起伏相對較小,均不超過5 dB,由此可見,耦合雖然改變了各階簡正波的幅度的模|am|,但是大部分簡正波幅度變化并不大.而從圖7(b)可見,耦合改變了所有階簡正波am的相位,這種改變,不僅對泄漏模態(tài)有影響,對受限模態(tài)也不例外.其中,對界面波幅度a1的相位變化最明顯,呈現(xiàn)出相位的鋸齒狀類周期性變化,說明界面波的相位受到以相鄰階簡正波為主的各階簡正波的強烈耦合影響.而第2—4 階簡正波am的相位也呈現(xiàn)類周期性變化.此外,當簡正波從受限模態(tài)逐步轉(zhuǎn)變?yōu)樾孤┠B(tài)后,其幅度am的相位隨距離的變化逐漸變小,最終相位趨于穩(wěn)定.其原因在于,轉(zhuǎn)變?yōu)樾孤┠B(tài)后,由于衰減系數(shù)的存在,耦合系數(shù)Bmn中的隨著距離而指數(shù)衰減,Bmn隨著距離而指數(shù)衰減,導致該階簡正波和其他簡正波之間的耦合作用隨著距離而指數(shù)衰減,因而其相位隨距離的變化逐漸減慢,最終趨于停滯,|am|的變化也存在類似規(guī)律.

圖7 zs=50 m,各階簡正波的(a) |am|和(b) am 的相位隨距離的變化;(c) f=50 Hz,海洋環(huán)境如圖2(c)時,各階簡正波本征值的實部隨距離的變化Fig.7.When zs=50 m,(a) |am| and (b) the phases of am of each order of normal mode varying with distance.(c) When f=50 Hz and the marine environment is shown in Fig.2(c),the real part of the normal mode eigenvalues of each order varying with distance.

為了進一步分析簡正波耦合對各階簡正波的影響,對各階簡正波的聲壓做FK 變換,得到圖8.

從圖8 可見,zs=50 m 時,相比絕熱近似,耦合顯著改變了界面波和第5 階簡正波聲壓的幅度和相位分布,而對其他階簡正波影響不大.而zs=100 m 時,耦合顯著改變了界面波和第3,5,7 階簡正波聲壓的幅度和相位分布.其原因在于,zs=100 m 時,聲源深度位于局地全海深的一半,第3,5,7 階簡正波的激發(fā)強度較小,容易受到其他階簡正波的耦合影響而改變幅度.

圖8 zs=(a) 50 和(b) 100 m 時,各階簡正波聲壓的FK 變換結(jié)果圖Fig.8.FK transformation results of normal mode sound pressure of each order when zs=(a) 50 and (b) 100 m.

3.2.4 海底坡度對聲場的影響

保持聲源處海深為200 m 不變,僅改變上坡海底的坡度,得到坡度為0.0125,0.025,0.05,0.075時的上坡海底時的聲壓,如圖9 所示,計算結(jié)束的距離分別為3400,3400,3400 和2267 m.

從圖9(a)—(d)可見,考慮泄漏模態(tài)時,各個海底坡度下,本文方法和COMSOL 所得聲壓基本一致;而不考慮泄漏模態(tài)時,海底坡度為0.0125時,在R＜ 1.2 km 時本文方法和COMSOL 所得聲壓有較大差異,而R≥ 1.2 km 時,本文方法和COMSOL 所得聲壓逐漸趨于一致,說明在該坡度下,泄漏模態(tài)對于近場有顯著影響,而在遠場逐步減弱,可以忽略.

而當坡度增加到0.025 時,不考慮泄漏模態(tài)時,除了R＜ 1 km 時本文方法和COMSOL 所得聲壓相差較大以外,x=2.9—3.3 km 處,本文方法和COMSOL 所得聲壓相差也較大,僅在x=1—2.9 km 段與COMSOL 所得聲壓基本一致.坡度進一步增加到0.05 和0.075 時,不考慮泄漏模態(tài)條件下,本文方法和COMSOL 所得聲壓在各個距離處都相差較大.其原因在于,以坡度為0.05 為例,結(jié)合前文的分析,第9—14 階簡正波在x=0 m 處就是泄漏模態(tài),只對近場聲壓有貢獻.而第2—8 階簡正波變?yōu)樾孤┠B(tài)的距離較遠,遠場處指數(shù)衰減距離并不長,因而它們對遠場聲壓的貢獻不可忽略.

海底坡度對簡正波耦合強度也有影響.對比這4 種不同海底坡度下的聲壓發(fā)現(xiàn),考慮泄漏模態(tài)時,R＜ 500 m 時,4 種坡度時所得聲場基本一致;而當R＞ 500 m 時,隨著海底坡度的增加,聲壓隨著距離的起伏變得越來越劇烈.其原因在于,隨著海底坡度的增加,簡正波耦合強度增強,簡正波耦合帶來的能量轉(zhuǎn)移的速度加快,從而帶來更為復雜的聲起伏,導致聲起伏的速度加快,表現(xiàn)為隨著海底坡度的增加,聲壓隨著距離的起伏變得越來越劇烈.從圖9(e)和圖9(f)可見,海底坡度從0.0125 增至0.025,x=1 km 處的簡正波耦合強度峰值增加了一倍,證實了海底坡度對簡正波耦合強度的影響.

圖9 上坡彈性海底的坡度為(a) 0.0125,(b) 0.025,(c) 0.05 和(d) 0.075 時的水中聲壓;考慮泄漏模態(tài),坡度為(e) 0.0125 和(f) 0.025,x=1 km 時的|Bmn(x)|Fig.9.Sound pressure in water when the seabed slope is (a) 0.0125,(b) 0.025,(c) 0.05 and (d) 0.075 on the upslope elastic seabed.The |Bmn(x)| considering the leaky mode when x=1 km when the slope is (e) 0.0125 and (f) 0.025.

3.2.5 海底地形對聲場的影響

以往的標準楔形海底考慮上坡海底的聲場,從上坡海底改為坡度為0.05 的下坡海底,海洋環(huán)境如圖10(a)所示,畫聲壓圖時以聲源位置為x=0 m.該條件下,本文方法和RAMS 以及COMSOL所得聲場如圖10(b)所示.

圖10(a) 下坡海底海洋環(huán)境;(b)下坡和(c)上坡的坡度為0.05 時的彈性海底聲場Fig.10.(a) Downslope seabed marine environment;sound field in water with (b) downslope and (c) upslope elastic seabed with a slope of 0.05.

從圖10(b)可見,下坡海底條件下,考慮泄漏模態(tài)時,本文方法和COMSOL 所得聲壓基本一致;而不考慮泄漏模態(tài)時,本文方法和COMSOL 所得聲壓相差較大.從圖10(c)也可以看到上坡海底時的相似結(jié)果,顯示了泄漏模態(tài)對水中聲壓不可忽略的貢獻.而比較圖10(b)和圖10(c)可見,上坡海底時泄漏模態(tài)對聲場的貢獻更為顯著,其原因在于,下坡海底情況下,聲源位置處海深只有50 m,只激發(fā)起了4 階簡正波,而上坡海底時聲源位置處海深200 m,激發(fā)起14 階簡正波.雖然耦合對聲場起到一定的影響,但是下坡時的4 階簡正波主導了聲場,而這4 階簡正波中并不包含泄漏模態(tài),只是在距離增加的過程中,能量通過簡正波耦合進入泄漏模態(tài).而上坡海底時,泄漏模態(tài)一開始就存在于14 階簡正波中,因而比下坡海底時對聲場有更大的貢獻.

此外,對比圖10(b)和圖10(c)可見,坡度絕對值都是0.05 的情況下,上坡情況下的聲壓傳播損失隨距離的變化起伏比下坡時更為劇烈.原因是,上坡海底時激發(fā)起的簡正波階數(shù)遠少于下坡海底時,因而上坡海底時簡正波之間的干涉比下坡時復雜得多,聲壓傳播損失隨距離的起伏更為劇烈.同時也說明,下坡海底時,即便因為簡正波耦合效應的存在,而在海深隨著距離變深的情況下激發(fā)起了更多階局地簡正波,這些受激發(fā)而產(chǎn)生的新簡正波的強度,也比原始的4 階簡正波弱得多,因而無法對聲場的干涉起到主要貢獻.

另外圖10(b)特別畫出了下坡海底時的RAMS聲壓計算結(jié)果.其目的在于展示出,下坡海底情況下,RAMS 的聲壓結(jié)果,和COMSOL 以及本文方法有較大差別,說明此時RAMS 的聲壓計算結(jié)果錯誤.再結(jié)合表1 可知,本文提出的耦合矩陣法,在3 種方法中,計算聲壓速度相對較快,而且能在下坡海底情況下,RAMS 計算聲壓出錯時,依然得到較為準確的水中聲壓,表明了提出本文方法的必要性.

3.2.6 海底橫波聲速對聲場的影響

改變cs,計算水中聲壓,所得結(jié)果如圖11 所示.從圖11(a)—(c)可見,考慮泄漏模態(tài)時,本文方法所得聲壓,在3 種不同海底橫波聲速的情況下,均與COMSOL 計算結(jié)果基本一致,表明了結(jié)果的準確性.而不考慮泄漏模態(tài)時,本文方法所得聲壓在3 種海底橫波聲速時均與COMSOL 計算結(jié)果產(chǎn)生了較大差異,從而顯示了泄漏模態(tài)對聲壓的顯著貢獻.此外,對比可見,cs=2000 m/s 時,泄漏模態(tài)對聲場的影響最小.其原因在于,在海底橫波聲速大于水中聲速的前提下,海底橫波聲速cs越大,一部分高階模態(tài)因為cs變大而從泄漏模態(tài)變?yōu)槭芟弈B(tài),泄漏模態(tài)的階數(shù)減少,因而泄漏模態(tài)對聲場的影響隨著海底橫波聲速的增大而減小.

圖11 cs=(a) 1800,(b) 1900,(c) 2000 m/s 時的水中聲壓;(d) 考慮泄漏模態(tài)條件下,cs=1800,1900,2000 m/s 時采用耦合簡正波法得到的聲壓對比Fig.11.Sound pressure in water when cs=(a) 1800,(b) 1900,(c) 2000 m/s;(d) sound pressure comparison obtained by mode coupling method while condidering leaky modes when cs=1800,1900,2000 m/s.

而從圖11(d)可見,在最初的500 m 距離,3 種橫波聲速情況下的水中聲壓基本一致,而后續(xù)的聲壓呈現(xiàn)出相似的起伏特性.然而,對比可見,對于這些相似的起伏,cs=2000 m/s 時聲壓隨距離的起伏最快(峰值最先出現(xiàn)),而cs=1800 m/s 時則是起伏最慢,聲壓隨距離的速度隨著橫波聲速的變小而變慢.此外,聲壓傳播損失隨著海底橫波聲速的增加而減小.其原因在于,上文提到,一部分高階模態(tài)因為cs變大而從泄漏模態(tài)變?yōu)槭芟弈B(tài),能量增強,表現(xiàn)為聲壓傳播損失減小;而這部分變?yōu)槭芟弈B(tài)的高階簡正波,對聲場干涉的影響增強,而高階模態(tài)之間的干涉的空間周期比低階模態(tài)之間的干涉空間周期短,呈現(xiàn)出聲壓隨著距離的快速起伏.

當海底橫波聲速小于水中聲速時,考慮典型的彈性海底環(huán)境參數(shù),取cp=1700 m/s,cs=600,100,800 m/s,其余環(huán)境參數(shù)如圖2(c)所示,zs=zr=30 m,計算f=50 Hz 時的水中聲場,如圖12所示.此時,由于cs較小,ks較大,界面波之外都是泄漏模態(tài),只考慮所有模態(tài)而不再去掉泄漏模態(tài),只對比不同算法得到的聲壓,以及不同cs對聲壓的影響.

圖12 cp=1700 m/s,cs=(a) 600,(b) 700 和(c) 800 m/s 時的水中聲壓;(d) cs=600,700,800 m/s 時,考慮泄漏模態(tài)情況下,采用耦合簡正波方法所得聲壓對比;(e) cs=600,700,800 m/s 時第2 階簡正波衰減系數(shù)隨距離的變化Fig.12.Sound pressure in water when cs=(a) 600,(b) 700,(c) 800 m/s;(d) sound pressure obtained by mode coupling method considering leaky modes and when cs=600,700,800 m/s;(e) the attenuation coefficient of the 2nd normal mode varying with range when cs=600,700,800 m/s.

從圖12(a)—(c)可見,當海底橫波聲速小于水中聲速時,耦合矩陣法所得聲壓與COMSOL 計算結(jié)果基本一致,在局部出現(xiàn)了一定的誤差.而從圖12(d)可見,隨著cs的增大,聲起伏隨距離的變化加快,其原因是和cs大于水中聲速時的原因一致.有所不同的是,聲傳播損失卻隨著cs的增大而變大,與cs大于水中聲速時的情況相反.其原因在于,在cs小于水中聲速時,隨著cs的增大,計算得到的各階簡正波本征值的衰減系數(shù)變大:以第2 階簡正波為例,3 種cs情況下,第2 階簡正波的衰減系數(shù)如圖12(e)所示,圖中可見,隨著cs的增大,第2 階簡正波的衰減系數(shù)變大,其余各階簡正波也存在此規(guī)律,此處不再畫圖顯示.

此外,我們注意到,圖12(a)中,大約在r=3050 m 處,模型結(jié)果出現(xiàn)奇異值,其原因在于,cp=1700 m/s,f=50 Hz 時,海底縱波波數(shù)kp≈0.1847996 rad/m,而本文模型計算得到的第3 階簡正波局地本征值在該點的值k3(r=3050 m) ≈(0.1847996—0.013628 i) rad/m,它的實部十分接近kp,這樣的情況下,模型的本征值計算結(jié)果出現(xiàn)了微小的誤差,但是這個微小的誤差,導致特征方程在本征值靠近kp時產(chǎn)生巨大誤差,從而導致了聲場計算結(jié)果奇異值的產(chǎn)生.圖12(b)中存在類似的問題,聲場奇異值處對應的也是第3 階簡正波本征值實部接近kp的情況,目前該問題尚未得到解決.

4 結(jié)論

采用RAMS,COMSOL 和簡正波耦合模型,本文分析了等聲速彈性傾斜海底條件下,聲場中的簡正波耦合以及由此帶來的能量轉(zhuǎn)移,發(fā)現(xiàn),本文模型所得聲壓計算結(jié)果和拋物方程所得結(jié)果基本一致;標準楔形海底算聲壓,必須考慮泄漏模態(tài),否則會帶來較大的誤差;海底橫波聲速會通過影響泄漏模態(tài)的階數(shù),影響了聲壓的起伏,以及傳播損失的大小;泄漏模態(tài)導致耦合系數(shù)變?yōu)閺蛿?shù),使得簡正波耦合不僅包含能量轉(zhuǎn)移,還包含相位變化,影響了簡正波聲壓的相位和幅度變化規(guī)律,從而導致了復雜的聲場干涉;同等坡度情況下,泄漏模態(tài)在上坡海底比下坡海底時對水中聲壓的貢獻更顯著.

此外,Pekeris 割線附近的本征值求解問題,本文的解仍存在一定的誤差.本文算法所得聲壓傳播損失和COMSOL 方法相比仍有一定的誤差,需要解決.這些問題,將會在今后的研究中進行解決.

感謝COMSOL 公司對本文中的聲場計算提供的技術(shù)支持.

附錄A彈性系數(shù)矩陣

矩陣算子A中各個參數(shù)含義:

附錄B復數(shù)本征值正交性推導

根據(jù)(13)式和(14)式,得到標量積〈um,Aun〉表達式:

采用分部積分處理(B1)式,得到

采用和(B1)式與(B2)式相同的處理方法,得到

將(B2)式和(B3)式疊加,得到

根據(jù)邊界條件(11)式可知,[wn]=[t3,n]=0,可知(B4)式右側(cè)為零,得到

再結(jié)合(13)式,由(B5)式得到適用于復本征值kn(x)的正交性表達式:

附錄C耦合系數(shù)推導與計算

1)耦合系數(shù)的推導

其中位移和應力局地簡正波本征函數(shù)表達式通過經(jīng)典的簡正波方法[35]求得.

2)能量守恒條件的推導

對于能量守恒條件[21]

得到,在垂直邊界x=xj處,

由此得到

3)耦合系數(shù)的簡化計算

根據(jù)分部積分得到

同理可得

這里 [·] 表示括號內(nèi)值在界面處的躍變,從底部到頂部.由此得到

附錄D復數(shù)本征值對應的各階簡正波的位移和應力

其中下角標1 和2 分別表示水中和海底的參量;u和w表示水平和垂直位移;τzz表示正應力,τzx表示剪應力;系數(shù)An,Dn和En根據(jù)經(jīng)典的簡正波理論[35]求得.