一種可穩(wěn)定計(jì)算負(fù)躍層淺海環(huán)境下聲場(chǎng)的波數(shù)積分方法

于曉林，許偉杰，楊春梅，駱文于

(1.中國(guó)科學(xué)院聲學(xué)研究所東海研究站，上海 201800；2.自然資源部第一海洋研究所，山東青島 266100；3.中國(guó)科學(xué)院聲學(xué)研究所聲場(chǎng)聲信息國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100190)

0 引言

在夏、秋季節(jié)，淺海中常存在較強(qiáng)的負(fù)躍層。在負(fù)躍層海域，海水可以分成三層，最上層以及最下層的聲速隨深度變化緩慢，兩層之間的中間層的聲速變化劇烈。淺海負(fù)躍層對(duì)水下聲監(jiān)測(cè)、漁業(yè)生產(chǎn)、艦艇航行等有極大的影響[1]，因此對(duì)淺海負(fù)躍層聲傳播規(guī)律的研究極為重要[2]。

波數(shù)積分方法首先是由Pekeris引入到水聲學(xué)中的[3]，簡(jiǎn)正波[4-5]方法是通過(guò)復(fù)圍線積分法計(jì)算該積分，得到的是不同階簡(jiǎn)正波對(duì)應(yīng)的留數(shù)之和，而波數(shù)積分方法是通過(guò)直接對(duì)波數(shù)進(jìn)行數(shù)值積分求解的，國(guó)內(nèi)外許多學(xué)者對(duì)波數(shù)積分方法進(jìn)行了相應(yīng)的研究[6-8]。在某些低頻淺海環(huán)境下，波數(shù)積分算法的計(jì)算結(jié)果比簡(jiǎn)正波方法更準(zhǔn)確[9]。

文獻(xiàn)[10]提出了一個(gè)負(fù)躍層淺海環(huán)境中的波數(shù)積分解法，該解法在理論上是正確的，但是在實(shí)際的仿真中發(fā)現(xiàn)，當(dāng)水平波數(shù)的取值過(guò)大時(shí)，會(huì)造成艾里(Airy)函數(shù)的數(shù)值溢出，系數(shù)矩陣存在數(shù)值不穩(wěn)定的現(xiàn)象。為了解決文獻(xiàn)[10]中給出的解法中存在的數(shù)值不穩(wěn)定的問(wèn)題，本文通過(guò)合理地歸一化，提出了一種負(fù)躍層淺海環(huán)境下無(wú)條件穩(wěn)定的波數(shù)積分方法。

1 理論推導(dǎo)

1.1 波數(shù)積分理論

1.2 負(fù)躍層淺海環(huán)境下不穩(wěn)定的波數(shù)積分算法

文獻(xiàn)[10]給出了計(jì)算負(fù)躍層淺海環(huán)境下的深度格林函數(shù)的線性方程組，該線性方程組在理論推導(dǎo)上是沒(méi)問(wèn)題的，但是由于未對(duì)Airy函數(shù)做任何的變換，因此在實(shí)際的仿真中，系數(shù)矩陣存在數(shù)值不穩(wěn)定現(xiàn)象。下面簡(jiǎn)單介紹下文獻(xiàn)[10]中給出的線性方程組并分析其不穩(wěn)定的原因，然后通過(guò)對(duì)Airy函數(shù)合理的變換，從而求得穩(wěn)定的深度格林函數(shù)。

負(fù)躍層淺海環(huán)境中的聲速如式(5)所示。

其中：cu為第一層中的聲速，cd為第三層中的聲速，cb為海底聲速，ρw為海水密度，ρb為海底密度，H為海水深度，H1、H2分別為負(fù)躍層的起、止深度。

以聲源位于第一層為例，各層中的線性方程可以表示如下[11]：

該矩陣方程的系數(shù)矩陣?yán)锎嬖谮呌跓o(wú)窮大的項(xiàng)Bi(ζ)以及Bi'(ζ)，因此可能存在數(shù)值溢出的現(xiàn)象。在仿真中發(fā)現(xiàn)，當(dāng)水平波數(shù)取值大于第三層海水波數(shù)k3時(shí)，仿真計(jì)算會(huì)提示矩陣接近奇異值、結(jié)果可能不準(zhǔn)確的警告?？梢?jiàn)，由于Airy函數(shù)的漸近特性，該深度格林函數(shù)的解在數(shù)值上是不穩(wěn)定的，為了解決該問(wèn)題，需要在系數(shù)項(xiàng)里避免趨于無(wú)窮大的項(xiàng)Bi(ζ)以及Bi'(ζ)，因此本文提出了一種無(wú)條件穩(wěn)定的波數(shù)積分算法。

1.3 負(fù)躍層淺海環(huán)境下穩(wěn)定的波數(shù)積分算法

由于Airy函數(shù)的漸近特性，上述得到的矩陣方程存在數(shù)值不穩(wěn)定的問(wèn)題，下面給出一種數(shù)值穩(wěn)定的計(jì)算方法。負(fù)躍層淺海環(huán)境中，聲源位于第一層、第三層、第二層時(shí)的示意圖如圖1～3所示。

圖1 聲源處于第一層Fig.1 The sound source in the first layer

1.3.1 聲源處于第一層時(shí)的矩陣方程

為了避免系數(shù)矩陣的系數(shù)項(xiàng)是奇異矩陣，需要對(duì)Airy函數(shù)進(jìn)行變換，將Airy函數(shù)表示成式(22)中的形式[12]：

圖2 聲源處于第三層Fig.2 The sound source in the third layer

圖3 聲源處于第二層Fig.3 The sound source in the second layer

該矩陣方程的系數(shù)項(xiàng)不存在趨于無(wú)窮大的項(xiàng)Bi(ζ)以及Bi'(ζ)，是無(wú)條件穩(wěn)定的，可以求得未知變量，然后通過(guò)式(6)～(9)求得聲源在第一層時(shí)，各個(gè)層的深度格林函數(shù)，最后利用式(3)可以求得聲場(chǎng)的位移勢(shì)能，從而求得各個(gè)層的聲場(chǎng)。

1.3.2 聲源處于第三層時(shí)的矩陣方程

當(dāng)聲源處于第三層時(shí)，同樣可以得到聲源位于第三層時(shí)的無(wú)條件穩(wěn)定的矩陣方程：

此時(shí)：A、x與式(30)中取值相同，b的取值如下：

通過(guò)該矩陣方程，可以求得7個(gè)未知變量，進(jìn)而求得總聲場(chǎng)。

1.3.3 聲源處于第二層時(shí)的矩陣方程

當(dāng)聲源處于第二層時(shí)，需要在聲源深度處構(gòu)建一個(gè)虛擬界面，同時(shí)在該虛擬邊界z=zs處添加兩個(gè)邊界條件：

從而得到聲源位于第二層的矩陣方程：

此時(shí)：

利用該矩陣方程，可以求出9個(gè)未知變量，進(jìn)而求得總聲場(chǎng)。

1.4 聲場(chǎng)計(jì)算

通過(guò)前面的推導(dǎo)，可以得到聲源位于三個(gè)層中的無(wú)條件穩(wěn)定的深度格林函數(shù)，再利用式(3)可以得到聲場(chǎng)的位移勢(shì)。在實(shí)際的仿真中，對(duì)格林函數(shù)尖峰的欠采樣，容易引入較大的誤差。為了解決這一問(wèn)題，采用復(fù)圍線積分的方法[11]：

式中：M是采樣點(diǎn)數(shù)，kmax為數(shù)值積分的上限，kmin為數(shù)值積分的下限。

2 數(shù)值算例

考慮圖4的海洋環(huán)境，海洋環(huán)境為負(fù)躍層淺海環(huán)境，海水深度為100 m，海水密度為1.0 g·cm-3，0～25 m的聲速為1 500 m·s-1，25～45 m的聲速由1 500 m·s-1變?yōu)? 460 m·s-1，45～100 m的聲速為1 460 m·s-1，海底聲速為1 800 m·s-1，海底密度為1.8 g·cm-3，海底吸收系數(shù)為0.15 dB·λ-1。

圖4 負(fù)躍層淺海聲場(chǎng)計(jì)算實(shí)例的海洋聲速剖面示意圖Fig.4 A typical sound speed profile for sound field calculation in the shallow-water with a negative thermoclines

分別將聲源置于第一層等聲速層以及第二層負(fù)躍層中，比較數(shù)值不穩(wěn)定的算法與本文所提出的數(shù)值穩(wěn)定的解法的結(jié)果差異：對(duì)于數(shù)值不穩(wěn)定算法來(lái)說(shuō)，水平波數(shù)最大取值只能取到第三層海水中波數(shù)k3附近，其中k3=ω/cd，當(dāng)水平波數(shù)的取值超過(guò)k3值時(shí)，系數(shù)矩陣存在數(shù)值不穩(wěn)定的現(xiàn)象，線性方程求解器在求解時(shí)會(huì)發(fā)生預(yù)警；而對(duì)于優(yōu)化后的算法來(lái)說(shuō)，水平波數(shù)可以取到一個(gè)極大值。在實(shí)際的仿真計(jì)算中，水平波數(shù)的取值應(yīng)該根據(jù)算例的具體情況來(lái)確定：理論上，為保證格林函數(shù)的收斂性，水平波數(shù)的取值應(yīng)該越大越好，但是為保證計(jì)算效率，水平波數(shù)取值又應(yīng)該越小越好。對(duì)于本算例來(lái)說(shuō)，1.2k3是一個(gè)可以保證格林函數(shù)收斂到0值附近的極小取值，因此本文中優(yōu)化后算法的水平波數(shù)一律取到1.2k3。圖5是聲源位于20 m(聲源在第一層)時(shí)，各層的深度格林函數(shù)。在本算例中，數(shù)值不穩(wěn)定算法的深度格林函數(shù)的水平波數(shù)最大值只能取到k3(即0.215)，優(yōu)化后的算法的水平波數(shù)最大取值為1.2k3(即0.258)時(shí)，就可以保證該算例下深度格林函數(shù)的收斂性。圖6給出相應(yīng)的傳播損失結(jié)果。

圖6 聲源深度為20 m，不同接收深度處的傳播損失計(jì)算結(jié)果（穩(wěn)定與不穩(wěn)定波數(shù)積分算法對(duì)比）Fig.6 The transmission losses at different depths for the sound source at the depth of 20 m(comparison between stable and unstable wave-number integral algorithms)

圖7是聲源位于40 m(聲源在第二層)時(shí)，不同深度的深度格林函數(shù)，可以看到，當(dāng)接收深度為40 m時(shí)，數(shù)值不穩(wěn)定算法的深度格林函數(shù)不能收斂到0附近，而優(yōu)化后的算法可以保證深度格林函數(shù)的足夠收斂性。圖8給出了相應(yīng)的傳播損失結(jié)果。

圖7 聲源深度為40 m，不同接收深度處的深度格林函數(shù)計(jì)算結(jié)果（穩(wěn)定與不穩(wěn)定波數(shù)積分算法對(duì)比）Fig.7 The magnitudes of the Green’s function at different depths for the sound source at the depth of 40 m(comparison between stable and unstable wave-number integral algorithms)

圖8 聲源深度為40 m，不同接收深度處的傳播損失的計(jì)算結(jié)果(穩(wěn)定與不穩(wěn)定波數(shù)積分算法對(duì)比)Fig.8 The transmission losses at different depths for the sound source at the depth of 40 m(comparison between stable and unstable wave-number integral algorithms)

為了進(jìn)一步驗(yàn)證本文所提方法的準(zhǔn)確性，比較本文所提方法的計(jì)算結(jié)果與簡(jiǎn)正波模型KRAKENC的計(jì)算結(jié)果[14]，圖9給出聲源深度在20 m(聲源位于第一層)時(shí)各個(gè)接收深度的傳播損失結(jié)果，圖10給出聲源深度在40 m(聲源位于第二層)時(shí)的各個(gè)接收深度的傳播損失結(jié)果。結(jié)果顯示，兩種方法的計(jì)算結(jié)果吻合得較好。由于簡(jiǎn)正波方法在某些低頻淺海環(huán)境中的近場(chǎng)計(jì)算上存在誤差[9]，因此本文所提方法可以作為負(fù)躍層淺海環(huán)境中聲場(chǎng)計(jì)算的有效補(bǔ)充。

圖9 聲源深度為20 m，不同接收深度處的傳播損失的計(jì)算結(jié)果(本文算法與KRAKENC算法對(duì)比)Fig.9 The transmission losses at different depths for the sound source at the depth of 20 m(comparison between the proposed algorithm and KRAKENC algorithm)

圖10 聲源深度為40 m，不同接收深度處的傳播損失的計(jì)算結(jié)果(本文算法與KRAKENC算法對(duì)比)Fig.10 The transmission losses at different depths for the sound source at the depth of 40 m(comparison between the proposed algorithm and KRAKENC algorithm)

3 結(jié)論

本文通過(guò)對(duì)Airy函數(shù)的處理，得到了負(fù)躍層淺海環(huán)境下無(wú)條件穩(wěn)定的波數(shù)積分方法。該方法保持了系數(shù)矩陣方程中的數(shù)值穩(wěn)定性，避免了由于數(shù)值不穩(wěn)定造成的深度格林函數(shù)的收斂性不足的問(wèn)題，能夠提高負(fù)躍層淺海環(huán)境下聲場(chǎng)計(jì)算的穩(wěn)定性。

將本文所提方法與簡(jiǎn)正波模型KRAKENC的計(jì)算結(jié)果進(jìn)行比較，結(jié)果顯示，兩種方法有很好的一致性。本文所提的方法在數(shù)值上是穩(wěn)定的，并且可以直接通過(guò)穩(wěn)定的系數(shù)矩陣方程求得相應(yīng)的解析解。因此，本文所提的方法可以作為負(fù)躍層淺海環(huán)境中聲場(chǎng)計(jì)算的有效補(bǔ)充。