劉 巍，肖汶斌，程興華，王勇獻(xiàn)，張理論

(國防科技大學(xué) 氣象海洋學(xué)院， 湖南 長沙 410073)

聲波能夠在海水介質(zhì)中遠(yuǎn)距離傳播，是目前獲取水下信息的主要載體，因此水下聲傳播特性的精確、高效計(jì)算始終是各軍事強(qiáng)國的重要研究方向[1]。水聲傳播的物理過程受波動(dòng)方程控制，由于聲壓是空間與時(shí)間的四維變量，直接求解計(jì)算量巨大。潛艇等水下兵器的聲頻譜一般具有窄帶特性，因而通常將時(shí)域內(nèi)的波動(dòng)方程采用傅立葉變換轉(zhuǎn)為頻域內(nèi)空間三維Helmholtz方程，每次只針對特定頻率的聲場進(jìn)行計(jì)算，使聲場求解難度有所降低[2]。由于Helmholtz方程屬于橢圓型方程，直接采用有限差分、有限元法離散需要建立方程組迭代求解，這對于追求時(shí)效性的水聲軍事應(yīng)用來講計(jì)算量仍顯過大，因此實(shí)際中經(jīng)常通過各類假設(shè)對Helmholtz方程進(jìn)行簡化，衍生出波數(shù)積分法、簡正波法、拋物方程法、射線法等水聲模型[3]。

對于水平分層的理想聲學(xué)介質(zhì)，波數(shù)積分法沒有模型誤差，被稱為“精確解”，廣泛應(yīng)用于海洋聲場仿真、海洋聲學(xué)參數(shù)反演與目標(biāo)定位匹配場反演[4-6]。然而，波數(shù)積分采用傳遞函數(shù)矩陣法計(jì)算積分核函數(shù)時(shí)會(huì)遇到計(jì)算不穩(wěn)定問題(出現(xiàn)NaN(not a number)而崩潰),雖然減小無限波數(shù)譜的截止波數(shù)可以提高穩(wěn)定性，但波數(shù)積分過程的數(shù)值精度也將隨之降低。數(shù)值計(jì)算穩(wěn)定性問題較難采用理論分析方法事先準(zhǔn)確判斷，故目前尚無自動(dòng)算法能夠獲得穩(wěn)定的、盡可能大的截止波數(shù)。

目前關(guān)于波數(shù)積分法計(jì)算穩(wěn)定性的研究較少，駱文于等[7]針對Pekeris波導(dǎo)問題(單層水體)，對深度采用局部坐標(biāo)，提出了一種穩(wěn)定的波數(shù)積分算法。本文對波數(shù)積分傳遞矩陣法自下而上(單向)、上下對進(jìn)(雙向)兩種計(jì)算順序進(jìn)行了分析，提出了基于預(yù)估-校正思想的最大截止波數(shù)自動(dòng)選取算法。

1 理論模型

在介質(zhì)水平分層假設(shè)下(介質(zhì)密度、聲速、吸收系數(shù)在各層內(nèi)部均勻，隱含聲壓二維軸對稱性質(zhì))，可采用Hankel變換將聲場Helmholtz方程轉(zhuǎn)換成如式(1)所示波數(shù)積分形式。

(1)

其中：P為頻域相對聲壓(參考點(diǎn)位于距聲源點(diǎn)1 m處)；J0(krr)為Bessel函數(shù)；φ(kr,z)為波數(shù)積分核函數(shù)[3]，且φ(kr,z)滿足深度方程

(2)

式中，k為聲場介質(zhì)波數(shù)，z為豎直方向坐標(biāo)，zs為聲源深度，δ為狄拉克函數(shù)。深度方程的求解是波數(shù)積分法的關(guān)鍵內(nèi)容，求解方法包括傳遞函數(shù)矩陣法、直接全局矩陣法、不變嵌入法等[3]。其中傳遞函數(shù)矩陣法具有解析形式，計(jì)算精度高、算法形式簡單，本文即針對該方法的不穩(wěn)定問題提出最大截止波數(shù)自動(dòng)選取算法。

1.1 傳遞函數(shù)矩陣

首先將聲場按照需要?jiǎng)澐殖扇舾伤綄?，并保證在聲源點(diǎn)深度zs存在一個(gè)分界面，各層與其分界面編號(hào)如圖1所示。

圖1 聲場介質(zhì)分層示意圖Fig.1 Sketch of acoustic multilayers

由于在各層聲介質(zhì)內(nèi)部聲學(xué)參數(shù)均勻且無聲源，因此在第m層的深度方程形式為

(3)

(4)

(5)

(6)

式中，ρm為第m層的密度，ω=2πf為角速度。若令Γm=kz,m/(ρmω)，則

(7)

進(jìn)而寫出第m層聲矢量

(8)

其中，

(9)

(10)

根據(jù)第m層聲矢量表達(dá)式可得第m層的上界面zm-1處聲矢量為

vm(kr,zm-1)=cm(kr,zm-1)am(kr)

(11)

第m層的下界面zm處聲矢量為

vm(kr,zm)=cm(kr,zm)am(kr)

(12)

將式(11)、式(12)結(jié)合可得聲矢量傳遞迭代式

vm(kr,zm-1)=Mm(kr)vm(kr,zm)

(13)

式中，Mm(kr)即為第m層的傳遞矩陣，

(14)

式中，hm=zm-zm-1代表第m層厚度，雙曲余弦、正弦函數(shù)表達(dá)式為

cosh(ikzhm)=(eikzhm+e-ikzhm)/2

(15)

sinh(ikzhm)=(eikzhm-e-ikzhm)/2

(16)

1.2 聲源界面條件與上、下邊界條件

1.2.1 聲源界面條件

聲源界面處的聲壓連續(xù)，但垂直振速不連續(xù)，存在點(diǎn)源條件

(17)

1.2.2 上邊界條件

在上邊界外側(cè)(第0層)

(18)

(19)

第0層下邊界的垂直振速為

(20)

令w0(kr,z0)與φ0(kr,z0)之比為

(21)

則從聲矢量中提出振速后形式為

(22)

需要指出的是，計(jì)算域上邊界以上部分(第0層)取不同密度值時(shí)，可表示不同的上邊界類型：①當(dāng)ρ0=0時(shí)(真空)，表示絕對軟邊界；②當(dāng)上邊界以上密度取第一層水體密度(即ρ0=ρ1)時(shí)，表示自由空間的開邊界；③當(dāng)上邊界以上密度取無窮大值(即ρ0→∞，實(shí)際取很大的正數(shù))時(shí)，表示絕對硬邊界，聲波不能進(jìn)入該介質(zhì)，邊界上質(zhì)點(diǎn)法向振速為0(當(dāng)ρ0→∞時(shí)，-Γ0→0、w0(kr,z0)→0)。

1.2.3 下邊界條件

在計(jì)算域下邊界外側(cè)(第N+1層)

(23)

(24)

在第N+1層上邊界zN處的垂直振速為

(25)

令φN+1(kr,zN)與wN+1(kr,zN)二者之比為

(26)

則從聲矢量中提出振速后形式為

(27)

與上邊界同理，下邊界以下部分(第N+1層)取不同密度值時(shí)，可表示不同的下邊界類型。

1.3 自下而上(單向)傳遞求解方法

傳遞函數(shù)矩陣法的聲矢量傳遞順序主要有兩種：一種是自下而上(單向)傳遞計(jì)算各界面聲矢量；第二種是分別從下邊界與上邊界(雙向)傳遞求解至聲源界面。本小節(jié)討論單向傳遞。

1)通過聲矢量迭代式從下至上計(jì)算聲源點(diǎn)界面緊下方聲矢量。在聲源下方區(qū)域，各分層界面處的聲壓與垂直振速應(yīng)連續(xù)，因此下層上界面的聲矢量與相鄰上層下界面的聲矢量相等，即

vm(kr,zm)=vm+1(kr,zm)

(28)

這樣即可通過聲矢量迭代式(13)一層一層地將最下層信息傳遞到聲源點(diǎn)界面緊下側(cè)，獲得vs+1(kr,zs)。

2)根據(jù)聲源界面條件計(jì)算出聲源點(diǎn)界面緊上側(cè)聲矢量。在聲源點(diǎn)界面，聲壓保持連續(xù)，垂直振速存在間斷，如式(17)，令

(29)

則聲源點(diǎn)界面緊上側(cè)聲矢量為

vs(kr,zs)=vs+1(kr,zs)-Δv(kr,zs)

(30)

3)從聲源點(diǎn)界面緊上側(cè)開始，繼續(xù)向上傳遞計(jì)算聲矢量，直至到達(dá)第一層上界面，獲得v1(kr,z0)。

4)根據(jù)最上層下界面聲向量v0(kr,z0)=v1(kr,z0)求出未定量wN+1(kr,zN)與w0(kr,z0)，進(jìn)而可計(jì)算出各位置的φ(kr,z)。

1.4 上下對進(jìn)(雙向)傳遞求解方法

單向與雙向傳遞求解方法實(shí)質(zhì)上都是以上邊界振速w0(kr,z0)與下邊界振速wN+1(kr,zN)為未知數(shù)，以分層界面聲壓連續(xù)、振速連續(xù)(除聲源界面外)建立兩個(gè)方程，進(jìn)行封閉求解。雙向傳遞求解方法為：

1)通過聲矢量傳遞迭代式，從下至上計(jì)算聲源界面緊下側(cè)聲矢量vs+1(kr,zs)。

2)通過聲矢量逆迭代式從上至下計(jì)算聲源界面緊上側(cè)聲矢量。逆?zhèn)鬟f迭代式為

(31)

通過自上而下逐層迭代可得vs(kr,zs)。

3)根據(jù)聲源界面條件(點(diǎn)源條件)計(jì)算出聲矢量未知數(shù)。在聲源界面上

vs+1(kr,zs)-vs(kr,zs)=Δv(kr,zs)

(32)

通過該式可計(jì)算出未知數(shù)wN+1(kr,zN)與w0(kr,z0)，進(jìn)而可計(jì)算出各位置的φ(kr,z)。

1.5 最大截止波數(shù)自動(dòng)選取算法

在數(shù)值計(jì)算過程中，需要對無限波數(shù)譜進(jìn)行適當(dāng)截止，并將積分離散為有限項(xiàng)求和

(33)

式中，kmax稱為截止波數(shù)。若kmax過大，可能導(dǎo)致計(jì)算發(fā)散[3]；若kmax過小，波數(shù)積分誤差就會(huì)增大(在聲源深度附近界面上誤差更大)。因此為了準(zhǔn)確計(jì)算各個(gè)深度的聲壓，就需要在保持計(jì)算穩(wěn)定的情況下，選取盡可能大的截止波數(shù)(下文稱其為最大截止波數(shù))。

2 算例測試

為便于直觀判斷計(jì)算結(jié)果的正確性與精確度，采用具有解析解的自由空間球面波算例進(jìn)行算法測試，聲壓解(R為到聲源點(diǎn)的距離)為

P=eikR/R

(34)

2.1 單向求解

采用單向求解方法，最大截止波數(shù)自動(dòng)選取，得kmax=2.9k=1.84。由于自由空間球面波聲源上、下兩側(cè)具有對稱性，因此聲源下方50 m與上方50 m深度上的積分核函數(shù)幅值隨水平波數(shù)變化曲線應(yīng)一致，如圖2所示，其中圖2(a)所示的聲源下方50 m處核函數(shù)幅值A(chǔ)bs(φ)是正確值，圖2(b)所示的聲源上方50 m處核函數(shù)“發(fā)散”。事實(shí)上，單向傳遞算法的核函數(shù)在通過聲源深度后很快就出現(xiàn)“發(fā)散”，使該位置以上的所有深度核函數(shù)的“發(fā)散”程度越來越大。

(a) 聲源下方50 m(z-zs=50 m)(a) 50 m under source depth(z-zs=50 m)

(b) 聲源上方50 m(z-zs = -50 m)(b) 50 m above source depth(z-zs= -50 m)圖2 球面波積分核函數(shù)幅值，單向求解(f=150 Hz)Fig.2 Wavenumber kernel function amplitude of free space spherical wave, bottom-to-top propagator(f=150 Hz)

如果各深度采用相同截止波數(shù)的傳統(tǒng)方法進(jìn)行積分，則單向求解的傳播損失云圖如圖3(a)所示，顯然在核函數(shù)發(fā)散的區(qū)域，聲場計(jì)算不正確。若在積分過程中，對不同深度采用不同的截止波數(shù)，可解決單向求解由核函數(shù)“發(fā)散”帶來的問題。在Re(kr)>kmm區(qū)間內(nèi)，正常情況下積分核函數(shù)絕對值應(yīng)隨Re(kr)的增加而減小，若在某一深度上積分核函數(shù)不斷增加，則表明出現(xiàn)了“發(fā)散”，可在此“發(fā)散”的水平波數(shù)位置停止積分。采用該技術(shù)后，單向求解所得傳播損失云圖如圖3(b)所示。

由于本算例有解析解，可以分析波數(shù)積分的數(shù)值誤差。在不同深度上采用不同的截止波數(shù)后(圖3(b)的狀態(tài))，單向求解的聲壓均方根誤差為RMSP=2.67E-3、傳播損失(TL=-20lg|P|, 單位dB)均方根誤差為RMSTL=3.59E-2。

(a) 所有深度均采用相同的截止波數(shù)(a) Same truncate wavenumbers at different depths

(b) 不同深度采用不同的截止波數(shù)(b) Different truncate wavenumbers at different depths圖3 球面波傳播損失云圖，單向求解(f=150 Hz)Fig.3 Acoustic transmission loss of free space spherical wave, bottom-to-top propagator(f=150 Hz)

2.2 雙向求解

采用雙向求解方法計(jì)算相同的算例，最大截止波數(shù)自動(dòng)選取算法得到kmax=5.0k=3.14(取到了初始預(yù)估值)，積分核函數(shù)在各深度上均未出現(xiàn)“發(fā)散”現(xiàn)象，聲場計(jì)算正確。聲壓均方根誤差為RMSP=2.62E-3、傳播損失均方根誤差為RMSTL=3.36E-2。與單向求解法相比，雙向求解法穩(wěn)定性更好、數(shù)值誤差更小。

3 結(jié)論

1)針對波數(shù)積分傳遞函數(shù)矩陣法計(jì)算不穩(wěn)定的問題(出現(xiàn)NaN而崩潰)，提出了最大截止波數(shù)自動(dòng)選取算法，可處理任意多層聲介質(zhì)，且具有簡單可靠、附帶計(jì)算量小的優(yōu)點(diǎn)；

2)深度方程單向求解法在越過聲源點(diǎn)深度繼續(xù)向上傳遞的過程中，易出現(xiàn)積分核函數(shù)非物理發(fā)散(未出現(xiàn)NaN，但數(shù)值異常)，需要在不同深度上采用不同的截止波數(shù)才能積分出正確結(jié)果；

3)與單向求解法相比，雙向求解法穩(wěn)定性更好，數(shù)值誤差更小。