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

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

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

0 引言

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

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

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

1 理論推導(dǎo)

1.1 波數(shù)積分理論

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

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

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

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

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

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

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

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

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

1.3.1 聲源處于第一層時的矩陣方程

為了避免系數(shù)矩陣的系數(shù)項(xiàng)是奇異矩陣，需要對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)不存在趨于無窮大的項(xiàng)Bi(ζ)以及Bi'(ζ)，是無條件穩(wěn)定的，可以求得未知變量，然后通過式(6)～(9)求得聲源在第一層時，各個層的深度格林函數(shù)，最后利用式(3)可以求得聲場的位移勢能，從而求得各個層的聲場。

1.3.2 聲源處于第三層時的矩陣方程

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

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

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

1.3.3 聲源處于第二層時的矩陣方程

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

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

此時：

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

1.4 聲場計算

通過前面的推導(dǎo)，可以得到聲源位于三個層中的無條件穩(wěn)定的深度格林函數(shù)，再利用式(3)可以得到聲場的位移勢。在實(shí)際的仿真中，對格林函數(shù)尖峰的欠采樣，容易引入較大的誤差。為了解決這一問題，采用復(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ù)躍層淺海聲場計算實(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é)果差異：對于數(shù)值不穩(wěn)定算法來說，水平波數(shù)最大取值只能取到第三層海水中波數(shù)k3附近，其中k3=ω/cd，當(dāng)水平波數(shù)的取值超過k3值時，系數(shù)矩陣存在數(shù)值不穩(wěn)定的現(xiàn)象，線性方程求解器在求解時會發(fā)生預(yù)警；而對于優(yōu)化后的算法來說，水平波數(shù)可以取到一個極大值。在實(shí)際的仿真計算中，水平波數(shù)的取值應(yīng)該根據(jù)算例的具體情況來確定：理論上，為保證格林函數(shù)的收斂性，水平波數(shù)的取值應(yīng)該越大越好，但是為保證計算效率，水平波數(shù)取值又應(yīng)該越小越好。對于本算例來說，1.2k3是一個可以保證格林函數(shù)收斂到0值附近的極小取值，因此本文中優(yōu)化后算法的水平波數(shù)一律取到1.2k3。圖5是聲源位于20 m(聲源在第一層)時，各層的深度格林函數(shù)。在本算例中，數(shù)值不穩(wěn)定算法的深度格林函數(shù)的水平波數(shù)最大值只能取到k3(即0.215)，優(yōu)化后的算法的水平波數(shù)最大取值為1.2k3(即0.258)時，就可以保證該算例下深度格林函數(shù)的收斂性。圖6給出相應(yīng)的傳播損失結(jié)果。

圖6 聲源深度為20 m，不同接收深度處的傳播損失計算結(jié)果（穩(wěn)定與不穩(wěn)定波數(shù)積分算法對比）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ù)，可以看到，當(dāng)接收深度為40 m時，數(shù)值不穩(wěn)定算法的深度格林函數(shù)不能收斂到0附近，而優(yōu)化后的算法可以保證深度格林函數(shù)的足夠收斂性。圖8給出了相應(yīng)的傳播損失結(jié)果。

圖7 聲源深度為40 m，不同接收深度處的深度格林函數(shù)計算結(jié)果（穩(wěn)定與不穩(wěn)定波數(shù)積分算法對比）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，不同接收深度處的傳播損失的計算結(jié)果(穩(wěn)定與不穩(wěn)定波數(shù)積分算法對比)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)確性，比較本文所提方法的計算結(jié)果與簡正波模型KRAKENC的計算結(jié)果[14]，圖9給出聲源深度在20 m(聲源位于第一層)時各個接收深度的傳播損失結(jié)果，圖10給出聲源深度在40 m(聲源位于第二層)時的各個接收深度的傳播損失結(jié)果。結(jié)果顯示，兩種方法的計算結(jié)果吻合得較好。由于簡正波方法在某些低頻淺海環(huán)境中的近場計算上存在誤差[9]，因此本文所提方法可以作為負(fù)躍層淺海環(huán)境中聲場計算的有效補(bǔ)充。

圖9 聲源深度為20 m，不同接收深度處的傳播損失的計算結(jié)果(本文算法與KRAKENC算法對比)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，不同接收深度處的傳播損失的計算結(jié)果(本文算法與KRAKENC算法對比)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é)論

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

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