多項式混沌展開的淺海水聲環(huán)境參數(shù)敏感度分析*

張 鵬 吳立新 胡治國 張雪冬

(1 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)

(2 中國科學院大學 北京 100049)

0 引言

水聲信道是一個時空頻變化的復雜多途信道。海洋環(huán)境特性隨著時間、季節(jié)、地域和海洋動力過程等因素而產(chǎn)生不確定的變化，這些不確定的變化必然會導致由海洋環(huán)境參數(shù)決定的水下聲場的不確定分布，從而引起水聲通信和探測系統(tǒng)性能的降低[1-2]。因此需要深入了解海洋環(huán)境參數(shù)不確定情況下聲場分布統(tǒng)計特性以及各個環(huán)境參數(shù)對聲場分布的影響程度，來實現(xiàn)聲吶系統(tǒng)的性能評估預報或者對聲吶系統(tǒng)的設計進行指導。

近些年來，海洋聲學中的不確定性研究受到了廣泛的關注，研究的重點在于環(huán)境信息不完全的情況下基于模擬的聲場不確定性預測，同時定量地評估不確定的環(huán)境參數(shù)對聲場的影響[3-16]。傳統(tǒng)的蒙特卡洛方法在存在多個不確定環(huán)境參數(shù)時所需計算量很大，且收斂較慢，而且該方法在進行環(huán)境參數(shù)敏感度分析時僅考慮單個變量的影響，不能計算不確定環(huán)境參數(shù)之間的交互效應[4-7]。還有一些學者應用聲場位移法來分析環(huán)境參數(shù)的敏感度，但該方法有嚴格的適用條件，適用范圍較窄[8-10]。代理模型方法(多項式混沌展開(Polynomial chaos expansion，PCE)法、Kriging法等)是一種常用的處理多維和非線性問題的有效方法，近年來被廣泛用于水聲不確定傳播問題的研究。已有的研究大多應用PCE或者Kriging法來分析預測不確定環(huán)境下的聲場分布，但是研究選取的不確定環(huán)境參數(shù)數(shù)量較少，隨著環(huán)境參數(shù)數(shù)量的增加，需要更多的多項式展開項數(shù)和訓練點數(shù)使得計算結果收斂，這會導致計算效率大幅度降低[11-16]。

本文考慮一個3 層模型的淺海海洋環(huán)境，根據(jù)設定的“淺海負梯度溫躍層”信道，定量分析了代表聲速剖面、水深和沉積層的8 個不確定環(huán)境參數(shù)對水下聲傳播的影響，同時根據(jù)射線理論從海底反射的角度解釋了各個環(huán)境因素影響程度差異的具體原因。

1 聲場分布的不確定性分析原理

水下聲場的不確定性分析包含隨機海洋環(huán)境參數(shù)概率分布的不確定性量化、聲場模型輸出響應集合的不確定性傳播以及聲場分布對于隨機海洋環(huán)境參數(shù)的敏感性分析。圖1顯示了基于PCE法的聲場不確定性分析過程。

圖1 基于PCE 方法的聲場不確定性量化分析Fig.1 Quantitative analysis of sound field uncertainty based on PCE model

1.1 多項式混沌展開

由于海洋環(huán)境參數(shù)空間分布的不均勻性、海洋動力過程導致的參數(shù)不確定性以及觀測數(shù)據(jù)的測量誤差，環(huán)境參數(shù)表現(xiàn)出服從一定概率密度函數(shù)的分布。將N個不確定的環(huán)境參數(shù)用集合來表示，對于給定的概率空間，聲場模型的輸出值可以展開為

其中，第一項為多項式混沌展開的均值，Γβ(ξ)為各階混沌展開項，αβ為對應各階混沌展開項的系數(shù)，ξi1，ξi2，ξi3，···，ξiN表示一組互相獨立的環(huán)境參數(shù)隨機變量。

如果用標準正態(tài)隨機變量來描述不確定參數(shù)變量ξ，則PCE 方法可以將模型響應描述為不確定參數(shù)的Hermite 多項式函數(shù)的展開式。對于其他類型的隨機變量，可以使用不同的多項式基或進行適當?shù)淖儞Q[17]。表1總結了多項式的選擇和隨機變量的分布類型之間的對應關系。

表1 隨機變量分布與正交多項式之間的關系Table 1 Classical univariate polynomial families used in PCE

多項式混沌展開的項數(shù)可由展開式的最高階數(shù)M以及隨機變量的個數(shù)N確定為

根據(jù)稀疏效應原則，只有輸入變量之間的低階相互作用才是最重要的，絕大部分高階相互作用的影響可以忽略。因此可以利用Q 范數(shù)來定義一個雙曲截斷方案：

式(3)中，

PCE方法的核心是展開系數(shù){αβ}β∈A的求解。根據(jù)公式(1)可以得到

式(5)中，P(ξ)是實際輸出，PPCE(ξ)為PCE 方法預測輸出，ε為殘差。

PCE 方法的系數(shù)求解問題可以構造為最小二乘優(yōu)化問題，使得殘差最小化：

1.2 Sobol敏感度分析

本文使用Sobol 敏感度指數(shù)估計全局非線性敏感度[18]。Sobol 指數(shù)和PCE 具有相似的多項式展開，并且都使用正交項。因此，一旦計算出PCE 的多項式系數(shù)，就可以直接得到各階Sobol 指數(shù)和總階Sobol 指數(shù)的解析解，而不需要任何額外的模型計算。根據(jù)混沌多項式的正交性可以得到

式(7)中，

則Sobol敏感度指數(shù)為

2 數(shù)值計算與分析

2.1 環(huán)境參數(shù)設置

本文選取如圖2所示的3層水聲信道參數(shù)模型，圖2中確定性參數(shù)為聲源頻率500 Hz，溫躍層下限深度D2= 70 m，海底附近聲速為1500 m/s，其他不確定參數(shù)及其分布如表2所示，其中沉積層的吸收系數(shù)根據(jù)Hamilton經(jīng)驗公式[19]來確定

表2 數(shù)值計算中使用的不確定性環(huán)境參數(shù)分布Table 2 The uncertain environmental parameters in numerical caculations

圖2 水聲信道環(huán)境參數(shù)模型Fig.2 Environmental parameter model of shallow water acoustic channel

2.2 聲傳播的不確定性分析

拋物方程近似聲場模型(Range-dependent acoustic model-parabolic equation，RAM-PE)是一種常用的聲場計算模型，本文采用該模型與多項式混沌展開法相結合來研究聲場的不確定性分布。

圖3給出了上發(fā)下收(聲源深度：溫躍層內(nèi)；接收深度：80 m)和下發(fā)下收(聲源深度：溫躍層以下；接收深度：80 m)這兩種收發(fā)情況在100 km 接收距離處PCE 方法與10000 次蒙特卡羅方法計算得到的傳播損失概率密度對比，兩種方法計算的結果較吻合，但是PCE 方法僅需調(diào)用聲場模型200 次即可滿足給定的精度閾值1%的均方根誤差，相較于蒙特卡羅方法的計算效率大幅提升。可以看到聲源位于溫躍層以下時，傳播損失均值和標準差較小，分布相對集中；而當聲源位于溫躍層內(nèi)時，由于溫躍層聲速的分布特性，聲線與海底的相互作用明顯增加，因此海洋環(huán)境的不確定性對傳播損失的影響更加劇烈。接下來將利用Sobol 敏感度指數(shù)來定量分析表2中各個序號對應的環(huán)境參數(shù)對聲場分布的影響程度。

圖3 不同聲源深度的傳播損失分布Fig.3 Transmission loss distribution of different sound source depth

2.3 環(huán)境參數(shù)敏感度分析

應用1.1 節(jié)中Sobol 敏感度分析方法，圖4給出了聲源分別位于溫躍層內(nèi)和溫躍層以下時在100 km 接收距離80 m 接收深度處的環(huán)境參數(shù)敏感度指數(shù)，兩種情況下沉積層聲速的不確定性對于該接收位置處的傳播損失的影響都是占絕對主導地位，但是聲源位于溫躍層內(nèi)時，聲源深度對聲場分布的不確定性也有一定的貢獻，其他因素的影響基本可以忽略。同時觀察到僅有環(huán)境參數(shù)的低階相互作用對聲場不確定性分布影響顯著，高階相互作用的影響可以忽略，這也符合稀疏效應原則，因此后續(xù)主要從總敏感度指數(shù)的角度來分析環(huán)境參數(shù)的影響程度。

圖4 不同聲源位置在100 km 距離上接收點處一階、二階及總敏感度指數(shù)Fig.4 The sensitivity indexes of different sound source depth at 100 km receiving distance

圖5為上發(fā)下收和下發(fā)下收時所有隨機參數(shù)的總敏感度指數(shù)隨傳播距離的變化。觀察發(fā)現(xiàn)兩種收發(fā)情況下在20 km 傳播距離內(nèi)，各個環(huán)境參數(shù)都對聲場的不確定性分布有較重要的影響，且隨著傳播距離的增加，除了沉積層聲速外的各個參數(shù)影響程度逐漸減小。

圖5 所有環(huán)境參數(shù)的總敏感度指數(shù)隨傳播距離的變化Fig.5 The Sobol sensitivity index varies with propagation distance

當傳播距離大于20 km 時，兩種收發(fā)情況下聲場分布對環(huán)境參數(shù)的敏感度呈現(xiàn)一定的差異：對于溫躍層內(nèi)的聲源，聲源深度的敏感度指數(shù)隨著傳播距離增加逐漸增大，并且在60 km 以后穩(wěn)定在0.25左右，沉積層聲速的敏感度指數(shù)總體保持在0.7 以上。這是由于躍層中聲速的負梯度使得聲線向海底傳播，進而增加了與海底的相互作用次數(shù)，本文并未考慮海面粗糙度以及波浪等海表面不確定性參數(shù)對聲傳播的影響，因此海底環(huán)境參數(shù)對聲場影響最大，其中海底沉積層聲速的影響最大。其他環(huán)境參數(shù)的影響非常小，基本可以忽略。當聲源位于溫躍層以下時，各個環(huán)境參數(shù)中，沉積層聲速的敏感度指數(shù)基本保持在0.85 以上，是聲場不確定性分布的最主要貢獻參數(shù)，其他環(huán)境參數(shù)的敏感度指數(shù)都低于0.1。因此對于20 km 以后的聲場分布，給定環(huán)境下僅需要重點考慮沉積層聲速的影響。

根據(jù)上述計算分析結果，選擇下發(fā)下收的情況研究不確定環(huán)境參數(shù)對不同頻率聲源聲場分布的影響。從圖6可以看到，當聲源位于溫躍層以下時，聲源深度對聲場分布的影響隨著頻率的升高不斷增加，與之相反的是沉積層聲速敏感度指數(shù)逐漸降低。在較低頻段內(nèi)(0～1700 Hz)，聲源深度的敏感度指數(shù)在0.8 以上，其他參數(shù)敏感度指數(shù)在絕大多數(shù)頻率下均保持在0.2以下。

圖6 所有環(huán)境參數(shù)的Sobol 敏感性指數(shù)隨頻率的變化Fig.6 The Sobol sensitivity index varies with source frequency

2.4 環(huán)境參數(shù)敏感度差異解釋

海水中聲信號經(jīng)過遠距離傳輸后的能量損失主要由聲傳播過程中波陣面的擴展帶來的擴展衰減、不均勻海水介質造成的吸收衰減以及海底海面邊界反射的損失組成。由于給定的信道環(huán)境未考慮海面的影響，而且通過敏感度指數(shù)的分析發(fā)現(xiàn)海底底質聲學參數(shù)對聲場分布影響較大，因此這里重點從射線角度分析海底沉積層參數(shù)的影響。

考慮如圖7所示的3 層模型[20]，圖中Rij和Tij分別表示反射系數(shù)和透射系數(shù)，下標對應傳播的方向。除了邊界處幅度的變化，還需要考慮相位的變化：

圖7 聲波在兩層海底模型中的反射Fig.7 Reflection and transmission by a layered halfspace

則總的反射系數(shù)

式(12)中，

其中，Zi為有效阻抗，

則公式(12)轉化為用有效阻抗來表示：

海底反射損失可以表示為

在考慮沉積層吸收時，公式(13)中有效阻抗中沉積層聲速和基底聲速可以分別表示為[20]

將表2中的不確定環(huán)境參數(shù)帶入式(15)求解海底反射損失，統(tǒng)計其在不同頻率和角度分布的標準差和均值結果，如圖8所示。

圖8 海底反射損失標準差和均值分布Fig.8 The standard deviation and mean distribution of bottom reflection loss

為了方便解釋，通過Bellhop 射線模型計算上發(fā)下收和下發(fā)下收兩種情況下典型聲源深度(50 m和80 m)對應的本征聲線和時間到達結構，如圖9所示，計算結果也可以反映兩種情況下不同聲源深度對應聲線和到達結構的總體趨勢。圖9中：藍綠色線表示掠射角小于7.5°的本征聲線和到達結構，紅色表示掠射角大于7.5°小于15°的本征聲線和到達結構，其他掠射角對應的聲線和到達結構用黑色線表示，結果表明兩種收發(fā)情況下掠射角均在15°以內(nèi)。本文中的掠射角特指聲線與海底相互作用時對應的掠射角，7.5°以內(nèi)稱為小掠射角，7.5°～15°之間稱為大掠射角。

從圖9可以看到，上發(fā)下收時，大掠射角本征聲線數(shù)量遠大于小掠射角聲線，而下發(fā)下收時的小掠射角本征聲線數(shù)量更多。小掠射角聲線具有更高的幅度，因此圖3(b)中下發(fā)下收時的傳播損失均值遠小于圖3(a)中上發(fā)下收的情況。由圖8(a)可見，聲源頻率為500 Hz 時，大掠射角聲線對應更大的海底反射損失標準差，因此圖3(a)中下發(fā)下收時傳播損失的分布標準差更大。

圖9 兩個不同深度聲源的本征聲線和時間到達結構Fig.9 Comparison of the eigenrays and arrivals at two different source depth

對比圖5中上發(fā)下收和下發(fā)下收兩種情況，通過Bellhop 射線模型計算不同接收距離處掠射角度。當收發(fā)距離在20 km 以內(nèi)時，各個環(huán)境因素不確定性對公式(14)中的掠射角度和阻抗等都具有一定程度的影響，且影響程度隨距離增加逐漸減小。

在20 km 以后較遠的傳播距離上，聲線到達結構與圖9中的計算結果較類似。上發(fā)下收時，溫躍層中聲速的負梯度使得聲線向海底傳播，聲源深度的變化會顯著地改變掠射角的分布位置，從而造成海底反射損失的較大波動，所以在圖5(a)中20 km以后聲源深度的敏感度指數(shù)較大。

對于下發(fā)下收的情況，遠距離處聲源深度的變化會一定程度影響小掠射角聲線的分布，從圖8中提取5°和11°這兩個典型掠射角對應的標準差、均值與頻率的關系，如圖10 所示。觀察圖8和圖10 可以看到，不同頻率聲源在小掠射角下的海底反射損失波動比大掠射角聲線的反射損失波動小很多，因此聲場分布的不確定性程度主要取決于大掠射角聲線海底反射損失波動。低頻段內(nèi)大掠射角對應的海底反射損失要比高頻段小很多，聲能量可以傳播到較遠距離處，所以低頻段內(nèi)海底沉積層對聲場分布的不確定性影響程度最大。隨著聲源頻率的增大，海底反射損失的均值也逐漸增加，大掠射角的聲線很難到達遠距離處，因此小掠射角的聲線貢獻了主要的聲場能量，此時聲源深度的影響程度也逐漸增強，對應圖6中隨頻率變化的環(huán)境參數(shù)敏感度。

圖10 掠射角為5°和11°時海底反射損失標準差、均值與頻率的關系Fig.10 The relationship between the standard deviation，mean value and frequency of bottom reflection loss

3 結論

本文將Sobol 敏感度指標的計算與混沌多項式展開方法相結合，根據(jù)設定的“水下負梯度溫躍層”信道模型，利用Q 范數(shù)約束的雙曲截斷方案來減少多項式展開項數(shù)，有效地分析了8 個海洋環(huán)境參數(shù)及其交互效應對不同聲源位置、傳播距離和聲源頻率的聲場不確定性分布的影響及Sobol 敏感度指數(shù)代表的影響程度，計算結果表明：(1)對于20 km 以內(nèi)較近的傳播距離內(nèi)，上發(fā)下收和下發(fā)下收兩種情況下各個環(huán)境參數(shù)都對聲場不確定性分布有較大的影響，在20 km后的遠場，沉積層聲速是聲場不確定性分布的最主要貢獻，同時上發(fā)下收時，聲源深度對聲場分布也有較大的影響。其他環(huán)境參數(shù)的影響基本可以忽略。(2)對于下發(fā)下收的情況，聲源深度對聲場分布的影響隨著頻率的升高不斷增加，與之相反的是沉積層聲速敏感度指數(shù)逐漸降低。

最后利用射線理論從海底反射的角度解釋了各個環(huán)境因素影響程度差異的具體原因。理論分析與計算結果較吻合，結果表明該方法可以有效地處理大量隨機環(huán)境參數(shù)情況下的聲場不確定性量化問題，并識別出對傳播損失最有影響力的參數(shù)。在聲吶應用中，根據(jù)實際使用條件，綜合考慮環(huán)境參數(shù)對聲場分布的敏感性，忽略敏感度較低的環(huán)境參數(shù)的變化性、關注對聲場不確定性影響較大的環(huán)境參數(shù)，可更精確預報和提升聲吶系統(tǒng)的性能。下一步將在更加復雜的情況，如三維海洋環(huán)境和內(nèi)波等水平變化環(huán)境，對文本所用分析方法的有效性進行驗證。