姚舟磊，于曉林，許偉杰

(1.中國(guó)科學(xué)院聲學(xué)研究所東海研究站，上海 201815；2.中國(guó)科學(xué)院大學(xué)，北京 100049)

0 引言

Pekeris 波導(dǎo)是Pekeris 提出的由一個(gè)均勻海水層和一個(gè)液體半無(wú)限空間組成的波導(dǎo)環(huán)境[1]。Pe‐keris波導(dǎo)可以看作是淺海環(huán)境的一個(gè)合理的簡(jiǎn)化模型，研究該模型中的聲傳播問(wèn)題具有重要的理論和實(shí)際意義[2-3]。

波數(shù)積分方法是聲場(chǎng)建模的主要方法之一，適用于水平分層的波導(dǎo)環(huán)境。Pekeris[1]在1948年首次將波數(shù)積分方法引入到水聲學(xué)中，利用波數(shù)積分方法來(lái)處理Pekeris波導(dǎo)中的聲傳播問(wèn)題。駱文于等[4]通過(guò)合理的歸一化，提出了一種數(shù)值穩(wěn)定的波數(shù)積分方法。

國(guó)內(nèi)學(xué)者對(duì)水下靜止聲源的聲場(chǎng)已有較多的研究[5-7]，但對(duì)運(yùn)動(dòng)聲源的聲場(chǎng)研究較少。如今魚雷等水下武器裝備的運(yùn)動(dòng)速度能達(dá)到80 m?s-1[8]，這對(duì)水下目標(biāo)跟蹤提出了更高的要求。對(duì)水下運(yùn)動(dòng)聲源的聲場(chǎng)研究能為聲吶測(cè)速和測(cè)距提供理論基礎(chǔ)，具有現(xiàn)實(shí)意義。國(guó)外對(duì)水下運(yùn)動(dòng)聲源的聲場(chǎng)特性研究已有不少成果，但多數(shù)是用簡(jiǎn)正波方法[9-10]和射線方法[11]進(jìn)行分析。Schmidt 等[12]給出了水平不變環(huán)境中聲源和接收器聯(lián)合運(yùn)動(dòng)時(shí)的波數(shù)積分和簡(jiǎn)正波表達(dá)式，并用修改的SAFARI模型進(jìn)行仿真。

本文對(duì)Pekeris 波導(dǎo)中的運(yùn)動(dòng)聲源聲場(chǎng)進(jìn)行分析，給出運(yùn)動(dòng)聲源二維波數(shù)積分的一維近似表示，并利用文獻(xiàn)[4]提出的數(shù)值穩(wěn)定波數(shù)積分方法對(duì)聲場(chǎng)進(jìn)行求解。首先在自由空間驗(yàn)證了近似的一維波數(shù)積分表示的準(zhǔn)確性，再分析Pekeris 波導(dǎo)中運(yùn)動(dòng)聲源的聲場(chǎng)特性。

1 理論推導(dǎo)

1.1 Pekeris波導(dǎo)中運(yùn)動(dòng)聲源的波數(shù)積分解

Pekeris波導(dǎo)環(huán)境示意圖如圖1所示，海水的密度和聲速為ρ1,c1,海水深度為D，海底的密度和聲速為ρ2,c2。聲源為單頻點(diǎn)源，以速度vs勻速運(yùn)動(dòng)。由于聲源的運(yùn)動(dòng)會(huì)破壞柱面的對(duì)稱性，因此引入笛卡爾坐標(biāo)系。

圖1 Pekeris波導(dǎo)環(huán)境示意圖Fig.1 Schematic diagram of Pekeris waveguide environment

為了計(jì)算方便，令x-y平面為海平面且x軸正方向?yàn)槁曉催\(yùn)動(dòng)方向，在t=0時(shí)刻聲源正處于z軸，位置記為(0,0,zs)。該聲源滿足如下的波動(dòng)方程[9]：

式中：ψ表示t時(shí)刻的聲壓，Ω為聲源的圓頻率，δ函數(shù)為狄拉克函數(shù)。利用傅里葉變換對(duì)：

可以得到運(yùn)動(dòng)點(diǎn)源波動(dòng)方程對(duì)應(yīng)的亥姆霍茲(Helmholtz)方程：

式中：k=ω/c，是頻率ω對(duì)應(yīng)的波數(shù)。

使用二維空間傅里葉變換可將Helmholtz 方程在空間-波數(shù)域上進(jìn)行轉(zhuǎn)換，其中x域的傅里葉變換對(duì)為

對(duì)式(4)進(jìn)行空間傅里葉變換可將x域轉(zhuǎn)換到kr波數(shù)域：

利用δ函數(shù)的性質(zhì)可對(duì)式(7)進(jìn)行化簡(jiǎn)，可以得到：

再將式(8)從y域變換到ky域可得：

式(9)是三維笛卡爾坐標(biāo)系下運(yùn)動(dòng)點(diǎn)源深度分離的波動(dòng)方程，其解的表達(dá)式為

其中：G(kx,ky,z,ω)是頻率為ω時(shí)的深度格林函數(shù)。對(duì)式(10)進(jìn)行傅里葉逆變換就能得到聲場(chǎng)的時(shí)域解：

將式(11)代入式(10)并化簡(jiǎn)，可得：

式(12)為運(yùn)動(dòng)點(diǎn)源聲場(chǎng)的二維波數(shù)積分解。式(12)中的雙重積分一般只能通過(guò)數(shù)值方法求解，而在遠(yuǎn)場(chǎng)的情況下，為了提高解的精度需要增大水平波數(shù)kx、ky的采樣數(shù)量，但會(huì)大大增加運(yùn)算時(shí)間。

因此在遠(yuǎn)場(chǎng)條件下，假設(shè)聲源與接收器都位于x-z平面，且聲源與接收器的水平距離遠(yuǎn)大于垂直距離，可認(rèn)為聲源與接收器距離是勻速變化的。利用穩(wěn)相法[16]可將式(12)的二維波數(shù)積分解近似為遠(yuǎn)場(chǎng)條件下的一維波數(shù)積分：

式中：Ω±kxvs表示聲源運(yùn)動(dòng)產(chǎn)生的多普勒頻移對(duì)深度格林函數(shù)造成的影響，當(dāng)聲源靠近接收器時(shí)多普勒頻移大于0，符號(hào)為正；相反當(dāng)聲源遠(yuǎn)離接收器時(shí)取負(fù)號(hào)。式中R(t)=R0?vst為t時(shí)刻時(shí)聲源與接收器的水平距離，R0為t=0時(shí)聲源與接收器的初始水平距離，當(dāng)聲源接近接收器時(shí)取負(fù)號(hào)，聲源遠(yuǎn)離接收器時(shí)取正號(hào)。

1.2 數(shù)值穩(wěn)定格林函數(shù)的表示

在之前的研究中，文獻(xiàn)[13]和文獻(xiàn)[14]給出了兩種Pekeris 波導(dǎo)中點(diǎn)源聲場(chǎng)的深度分離格林函數(shù)表達(dá)式：

文獻(xiàn)[4]指出，式(17)和式(18)理論上是正確的，但是在求解待定系數(shù)矩陣方程時(shí)，由于上下行波選取的參考點(diǎn)不合理，導(dǎo)致在實(shí)際仿真中系數(shù)矩陣容易產(chǎn)生數(shù)值溢出的情況，會(huì)使求解得到的深度格林函數(shù)不穩(wěn)定。文獻(xiàn)[4]通過(guò)合理的歸一化，將深度格林函數(shù)表示為

代入邊界條件(14)～(16)，可得線性方程組：

通過(guò)求解式(20)中的方程組，可得到Pekeris波導(dǎo)中數(shù)值穩(wěn)定的深度格林函數(shù)，在計(jì)算過(guò)程中再將頻率替換成由聲源運(yùn)動(dòng)所產(chǎn)生的多普勒頻移Ω±kxvs并代入式(13)，便可得到運(yùn)動(dòng)聲源數(shù)值穩(wěn)定的一維時(shí)域波數(shù)積分近似解。

2 仿真算例

第1 節(jié)給出了運(yùn)動(dòng)點(diǎn)源聲場(chǎng)的二維波數(shù)積分解。但由于二維解的計(jì)算量較大，因此在遠(yuǎn)場(chǎng)情況下引入一定的假設(shè)，將解近似為一維波數(shù)積分解。本節(jié)先模擬自由空間下運(yùn)動(dòng)點(diǎn)源的聲場(chǎng)，通過(guò)對(duì)比遠(yuǎn)場(chǎng)條件下的一維波數(shù)積分解與解析解，來(lái)驗(yàn)證式(13)的準(zhǔn)確性；再根據(jù)式(13)和式(20)的方程組對(duì)Pekeris 波導(dǎo)環(huán)境下的水下運(yùn)動(dòng)聲源的聲場(chǎng)進(jìn)行研究。

2.1 運(yùn)動(dòng)聲源聲場(chǎng)一維波數(shù)積分?jǐn)?shù)值解的驗(yàn)證

自由空間中單頻運(yùn)動(dòng)點(diǎn)源所產(chǎn)生的聲場(chǎng)存在解析解。假設(shè)點(diǎn)源頻率為f0，以速度v沿x軸勻速運(yùn)動(dòng)，t=0時(shí)刻位于原點(diǎn)，則聲場(chǎng)解析解為[15]

其中：

式中：M為聲源運(yùn)動(dòng)速度v與介質(zhì)聲速c的比v/c。

模擬自由空間參數(shù)如下：空間的介質(zhì)聲速為1 500 m?s-1，介質(zhì)密度為1 000 kg?m-3,運(yùn)動(dòng)點(diǎn)源的頻率為100 Hz，從原點(diǎn)(0,0,0)開(kāi)始沿x軸正方向、以100 m?s-1速度勻速運(yùn)動(dòng)，接收器坐標(biāo)為(10 000,0,50)。

圖2(a)、圖2(b)、圖2(c)分別是在0、20、40 s時(shí)計(jì)算解析解和數(shù)值解得到的聲壓幅值，對(duì)應(yīng)的是聲源與接收器距離10、8和6 km的時(shí)刻。圖2(d)是0～80 s接收器接收聲壓級(jí)的變化曲線。圖2(d)中解析解與數(shù)值解得到的聲壓級(jí)變化曲線吻合，誤差為0.3 dB左右。圖2中藍(lán)色實(shí)線是解析解，黃點(diǎn)是式(13)計(jì)算的數(shù)值解。可以看出在自由空間中，聲源與接收器相距較遠(yuǎn)的遠(yuǎn)場(chǎng)條件下，一維波數(shù)積分?jǐn)?shù)值解和解析解的計(jì)算結(jié)果符合的比較好。

圖2 不同時(shí)刻的接收波形和接收聲壓級(jí)的解析解與數(shù)值解對(duì)比Fig.2 Comparison between analytical and numerical solutions of the received waveforms and sound pressure levels at different times

因此，在遠(yuǎn)場(chǎng)條件下，一維波數(shù)積分的數(shù)值解的精度較高，能滿足運(yùn)動(dòng)聲源聲場(chǎng)的計(jì)算。在自由空間下運(yùn)動(dòng)聲源多普勒頻移的理論值為

式中：f0為聲源的頻率，θ(t)為t時(shí)刻聲源、接收器連線方向與聲源運(yùn)動(dòng)方向的夾角。當(dāng)水平距離遠(yuǎn)大于垂直距離時(shí)可忽略角度θ的影響，計(jì)算得到距離較遠(yuǎn)時(shí)的多普勒頻移理論值為107.14 Hz。圖3對(duì)比了0～99 s接收信號(hào)解析解、近似解的瞬時(shí)頻率和上式計(jì)算得到的多普勒頻移理論值，可以看到接收信號(hào)的瞬時(shí)頻率與理論值吻合較好。

圖3 信號(hào)瞬時(shí)頻率解析解、近似解和理論值的對(duì)比Fig.3 Comparison of the analytical and approximate solutions of signal instantaneous frequency and the theoretical values

2.2 Pekeris波導(dǎo)運(yùn)動(dòng)聲源數(shù)值穩(wěn)定解的仿真

將Pekeris 波導(dǎo)中數(shù)值穩(wěn)定深度格林函數(shù)代入式(13)，并進(jìn)行數(shù)值仿真。Pekeris波導(dǎo)仿真環(huán)境如圖4所示，其中海水聲速為1 500 m?s-1，密度為1 000 kg?m-3，海深為100 m，海底聲速為1 800 m?s-1，密度為1 800 kg?m-3，海底的吸收系數(shù)為0.2 dB?λ-1。在初始時(shí)刻聲源位于z軸，深度為36 m，坐標(biāo)為(0,0,36)。聲源頻率為50 Hz，以100 m?s-1速度沿x軸正方向運(yùn)動(dòng)。接收器深度為46 m，與聲源水平距離10 000 m，坐標(biāo)為(10 000,0,46)。

圖4 Pekeris波導(dǎo)環(huán)境的仿真Fig.4 Simulation of Pekeris waveguide environment

當(dāng)聲源靜止時(shí)，直接求解式(20)中的方程組得到該波導(dǎo)環(huán)境下靜止聲源的深度格林函數(shù)的幅度，結(jié)果如圖5所示。用KRAKENC聲場(chǎng)仿真軟件仿真該波導(dǎo)環(huán)境，計(jì)算得到該波導(dǎo)環(huán)境共會(huì)產(chǎn)生4階簡(jiǎn)正波。表1 是圖5 的4 個(gè)峰值幅度對(duì)應(yīng)的波數(shù)kx與KRAKENC求解得到的各階簡(jiǎn)正波水平波數(shù)值。表1顯示該波導(dǎo)環(huán)境下靜止聲源的深度格林函數(shù)峰值對(duì)應(yīng)的波數(shù)與各階簡(jiǎn)正波的水平波數(shù)相符合。

表1 各階簡(jiǎn)正波水平波數(shù)與深度格林函數(shù)峰值對(duì)應(yīng)波數(shù)對(duì)比Table 1 Comparison of the horizontal wavenumber of each order of normal waves and the corresponding peak wavenumber of depth Green's function

圖5 靜止聲源的深度格林函數(shù)幅度Fig.5 The amplitudes of the depth Green's function for the stationary sound source

當(dāng)聲源運(yùn)動(dòng)時(shí)，深度格林函數(shù)需要加上多普勒頻移，即式(13)中的Ω±kxvs，當(dāng)聲源靠近接收器時(shí)取負(fù)號(hào)，遠(yuǎn)離接收器時(shí)取正號(hào)。圖6是該波導(dǎo)環(huán)境下聲源不同運(yùn)動(dòng)狀態(tài)時(shí)深度格林函數(shù)的幅值，其中實(shí)線是聲源靜止時(shí)的深度格林函數(shù)的幅值，虛線是聲源以100 m?s-1速度遠(yuǎn)離接收器時(shí)深度格林函數(shù)的幅值，點(diǎn)劃線是以100 m?s-1速度接近接收器時(shí)深度格林函數(shù)的幅值。由圖6可以看出，當(dāng)聲源接近接收器時(shí)多普勒頻移使深度格林函數(shù)峰值對(duì)應(yīng)的波數(shù)增大；聲源遠(yuǎn)離接收器時(shí)深度格林函數(shù)的峰值波數(shù)減小。

圖6 聲源不同運(yùn)動(dòng)狀態(tài)時(shí)深度格林函數(shù)的幅度Fig.6 Amplitudes of the depth Green's function under different motion states of a sound source

求解式(13)得到時(shí)間區(qū)間為0～80 s時(shí)接收聲壓信號(hào)的波形，結(jié)果如圖7所示。其中藍(lán)線是聲源靜止時(shí)的接收信號(hào)波形，紅線是聲源以100 m?s-1速度勻速靠近接收器時(shí)接收的信號(hào)。

圖7 在0～80 s時(shí)間內(nèi)的接收波形Fig.7 The received waveform from 0 to 80 s

圖8(a)～8(c)分別對(duì)應(yīng)0、20、40 s時(shí)刻接收器的聲壓時(shí)域波形。由圖8可以看出，由于聲源與接收器水平距離隨時(shí)間發(fā)生變化，接收聲壓的幅度會(huì)發(fā)生起伏，并隨著聲源的接近幅度增大。由于聲源向接收器運(yùn)動(dòng)產(chǎn)生多普勒頻移，接收信號(hào)的頻率大于聲源靜止時(shí)接收信號(hào)的頻率。

圖8 不同時(shí)刻淺海Pekeris波導(dǎo)接收信號(hào)的時(shí)域波形Fig.8 Received signal waveforms in the shallow-sea Pekeris waveguide environment at different times

在淺海Pekeris 波導(dǎo)中，聲源產(chǎn)生的聲波會(huì)沿不同的路徑到達(dá)接收器，除了直達(dá)聲波外，還存在經(jīng)過(guò)海面和海底多次反射后的聲波。由于聲源的運(yùn)動(dòng)，不同傳播路徑的聲波會(huì)產(chǎn)生不同的多普勒頻移，2.1 節(jié)的多普勒頻移求解公式僅適用于對(duì)直達(dá)聲波的求解。

圖9 為計(jì)算得到聲源通過(guò)接收器正上方的50～150 s 時(shí)間區(qū)間接收信號(hào)的瞬時(shí)頻率與直達(dá)聲波的多普勒頻移理論值的比較，可以看到當(dāng)聲源靠近接收器時(shí)接收信號(hào)頻率大于原始頻率，遠(yuǎn)離時(shí)小于原始頻率。由于直達(dá)聲波具有最大的多普勒頻移，沿其他路徑傳播的聲波多普勒頻移小于直達(dá)聲波，因此當(dāng)聲源接近接收器時(shí)仿真得到接收信號(hào)的瞬時(shí)頻率略小于直達(dá)聲波的多普勒頻移。由于在Pekeris波導(dǎo)中運(yùn)動(dòng)聲源與接收器的水平距離隨時(shí)間發(fā)生變化，接收信號(hào)的幅度會(huì)隨時(shí)間發(fā)生起伏。在幅度較小處的瞬時(shí)頻率會(huì)發(fā)生較大的變化，因此在圖8中可以看到計(jì)算得到的瞬時(shí)頻率存在突變。

圖9 一維波數(shù)積分解得到的信號(hào)瞬時(shí)頻率與理論值對(duì)比Fig.9 Comparison between the signal instantaneous frequency obtained by the one-dimensional wavenumber integration solution and the theoretical value

截取并放大70～80 s 時(shí)刻的時(shí)域波形和瞬時(shí)頻率，如圖10所示。由圖10中的仿真結(jié)果可以看到瞬時(shí)頻率突變的時(shí)刻與波形幅度極小值相吻合。

對(duì)0～80 s的接收信號(hào)進(jìn)行傅里葉變換得到接收信號(hào)的頻譜，如圖11所示?？梢钥吹浇邮招盘?hào)的頻譜有多個(gè)譜峰，這是由于Pekeris 波導(dǎo)中存在多途效應(yīng)，聲波沿多條路徑到達(dá)接收器。通過(guò)KRAKENC模擬得到該波導(dǎo)環(huán)境下50 Hz的聲源會(huì)產(chǎn)生4 階簡(jiǎn)正波，不同的簡(jiǎn)正波有不同的相速度，因此也就有不同的多普勒頻移[9]，當(dāng)聲源的運(yùn)動(dòng)軌跡通過(guò)接收器的正上方，各階簡(jiǎn)正波的頻移不會(huì)發(fā)生改變，因此接收的信號(hào)為多個(gè)單頻信號(hào)，在圖11中顯示了多個(gè)譜峰，譜峰的位置是由波導(dǎo)條件決定的。

圖11 在0～80 s時(shí)間內(nèi)接收信號(hào)頻譜Fig.11 The spectrum of the received signal from 0 to 80 s

改變接收器的深度，并計(jì)算接收信號(hào)的頻譜，得到的結(jié)果如圖12所示，圖中橙色虛線和紫色點(diǎn)劃線分別為接收器在66 m、86 m 深度時(shí)接收信號(hào)的頻譜?？梢钥吹阶V峰的高低發(fā)生了變化，但是對(duì)應(yīng)的頻率并沒(méi)有改變。將聲源頻率改為35 Hz，接收信號(hào)的頻譜如圖13所示。利用Krakenc計(jì)算該波導(dǎo)環(huán)境的簡(jiǎn)正波函數(shù)，可以得到3階簡(jiǎn)正波，與圖13中的3個(gè)譜峰相對(duì)應(yīng)。

圖12 不同深度接收信號(hào)的頻譜Fig.12 Spectrums of the received signals at different depths

圖13 聲源頻率為35 Hz時(shí)接收信號(hào)的頻譜Fig.13 The spectrum of the received signal at the sound source frequency of 35 Hz

目前水下航行器和水下武器的運(yùn)行速度還難以達(dá)到100 m?s-1，因此在相同波導(dǎo)環(huán)境下對(duì)相對(duì)低速運(yùn)動(dòng)的聲源進(jìn)行仿真。聲源聲速為50 m?s-1和20 m?s-1時(shí)仿真得到0～80 s 接收器接收信號(hào)的時(shí)域波形如圖14所示。與圖7 相比，聲源運(yùn)動(dòng)速度減小，接收信號(hào)幅度的變化速度也變緩。

圖14 聲源不同速度運(yùn)動(dòng)時(shí)接收器的時(shí)域波形Fig.14 The time domain waveforms of the received signals when the sound source moves at different speeds

圖15 為不同運(yùn)動(dòng)速度時(shí)接收信號(hào)的時(shí)頻圖，可以看到當(dāng)聲源運(yùn)動(dòng)速度減小時(shí)，接收信號(hào)的多普勒頻移也會(huì)減小。從圖15 中的顏色深淺變化也能看出速度越快接收信號(hào)變化得越快。

圖15 聲源不同速度運(yùn)動(dòng)時(shí)接收信號(hào)的時(shí)頻圖Fig.15 Time-frequency diagrams of the received signals when the sound source moves at different speeds

3 結(jié)論

本文為求解Pekeris 波導(dǎo)中運(yùn)動(dòng)聲源的聲場(chǎng)提供一種新思路，將Pekeris 波導(dǎo)中靜止聲源的數(shù)值穩(wěn)定波數(shù)積分解應(yīng)用到運(yùn)動(dòng)聲源的聲場(chǎng)分析中。在仿真實(shí)驗(yàn)中，驗(yàn)證了二維波數(shù)積分的一維近似數(shù)值解的正確性，并仿真分析了Pekeris 波導(dǎo)下運(yùn)動(dòng)的單頻點(diǎn)源產(chǎn)生的聲場(chǎng)。波數(shù)積分方法是直接數(shù)值求積分的方法，本文只在遠(yuǎn)場(chǎng)條件下引入較小的誤差。利用數(shù)值穩(wěn)定的計(jì)算方法來(lái)求解運(yùn)動(dòng)聲源的深度格林函數(shù)，不存在數(shù)值溢出的問(wèn)題。因此本文的方法可以作為研究Pekeris 波導(dǎo)中運(yùn)動(dòng)單頻點(diǎn)源的遠(yuǎn)場(chǎng)特性的參考模型來(lái)應(yīng)用。