利用隨機多項式展開的海底聲學參數反演方法

李風華 王翰卓

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

2) (中國科學院大學電子電氣與通信工程學院, 北京 100049)

為了提高地聲反演算法的計算效率, 探索克服地聲反演結果多值性問題, 本文利用寬帶、多收發位置的傳播損失數據結合傳播損失在地聲參數先驗搜索區間內的隨機多項式展開系數矩陣, 反演得到海底縱波聲速、吸收率和密度比重.使用隨機多項式展開近似傳播損失時, 展開系數的自變量為聲波頻率、收發位置等參數, 隨機多項式的自變量為表示聲速、吸收率、比重在各自搜索區間內均勻分布的隨機變量.傳播損失的展開系數通過嵌入隨機多項式的聲學寬角拋物方程結合蓋遼金投影、最小角度回歸算法計算求得.在低頻、一定聲傳播水平距離以內和地聲參數搜索區間長度適中時, 使用隨機多項式展開近似傳播損失的相對誤差在1%以下.仿真發現, 在淺海環境中使用低頻、一定聲傳播水平距離以內的傳播損失數據, 在接收信號信噪比較高、聲源和水聽器相對位置誤差較小時, 選擇合適的隨機多項式展開截斷冪次可較準確地反演海底聲速、吸收率和密度比重, 且計算效率比網格遍歷搜索方法提高一個數量級以上.

1 引 言

海洋聲學中的逆問題是利用海水中聲波所攜帶的聲源位置及海洋信道環境特性信息, 對聲源位置和介質參數進行求解.地聲反演問題是快速獲取局部海域等效海底聲學參數的方法, 是水聲逆問題的重要研究方向[1].地聲反演過程需要進行以下幾個方面的內容: 海底環境模型選取, 待反演參數搜索范圍的確定, 拷貝聲場計算, 匹配物理量及代價函數的構造, 最優化算法的選取, 反演結果的評價.為求解待反演參數的全局最優值, 需要采用遍歷或最優化算法在待反演參數空間內反復選取待反演參數向量、計算其對應的拷貝聲場以及代價函數.傳統基于遍歷搜索方法所需計算量大, 計算時間長;代價函數對部分待反演參數的不敏感以及待反演參數在代價函數中的耦合效果使得代價函數存在多解性問題[1].提高反演速率以及克服多解性是地聲反演中需要解決的問題.

為了減少求解反演參數空間拷貝聲場的計算規模, Gerstoft[2,3]和Dosso[4,5]分別提出了基于遺傳算法和模擬退火的貝葉斯推斷法, 該方法可快速求解作為代價函數的后驗概率密度; Sambridge[6]提出了鄰域插值近似方法, 該方法利用集合信息指導參數空間的重采樣, 從而減少了計算量.

為了解決地聲反演中的多值性問題, Holland和Osler[7]提出了組合反演方法, 將空間-時間域和空間-頻率域的數據結合起來, 采用多種不同的地聲模型擬合同一數據集.李整林等[8?13]在綜合分析簡正波的頻散特性、海底反射系數、傳播損失等聲場特性的基礎上, 應用匹配場處理器、自適應時頻分析算法和并行遺傳算法提出多物理量聯合地聲反演方法.

隨機多項式展開法[14]是一種表示方程中的多維隨機參數在方程解中非線性傳遞函數關系的方法.近些年引入到海水聲速起伏[15?23]、海底參數不確知[22?25]、聲源和接收點深度不確知[25,26]的聲場建模中.隨機多項式展開方法將頻域復聲壓信號表示成以隨機多項式為基底的級數展開形式, 隨機多項式的自變量即為描述不確知(或隨機)環境的隨機變量.在地聲反演問題中, 假設待反演參數為在搜索區間內均勻分布的隨機變量[24], 通過隨機多項式展開可以得到聲場關于待反演海底參數的解析表達式.該方法的計算效率較遍歷方法高, 可以作為計算拷貝聲場的“代理”模型.

本文利用寬帶多收發位置的聲傳播損失數據,結合地聲參數搜索區間內傳播損失函數的隨機多項式展開系數矩陣, 求解隨機多項式基底的數值,利用基底中的一次冪項, 唯一地反演海底聲學參數.該方法比傳統遍歷搜索法計算速度快, 且對信號噪聲和收發距離隨機誤差有一定的穩健性.

2 海底參數不確知下聲場的隨機多項式展開

在頻率為f, 二維柱坐標水平和深度位置(r,z)處, 海底參數在各自區間內取值對應的向量為時, 頻域復聲壓可以用隨機多項式展開方法表示為

其中 k0是參考波數;{γq(f,r,z);q=0,1,···,Q?1}為隨機多項式展開的系數, 是頻率f和空間位置(r,z) 的函數;是隨機多項式展開基, 是隨機變量的函數; N是復聲壓的隨機多項式展開中展開基的截斷冪次; Q代表展開截斷冪次為N時, 復聲壓的隨機多項式展開項數.在均方意義下, 在截斷冪次 N →∞ 時, 頻率復聲壓的隨機多項式展開近似結果可以收斂到真實值, 誤差隨著截斷冪次N的增加指數減小[14].隨機多項式基底有勒讓德多項式的函數形式, 滿足加權正交歸一化條件:

其中 μ 為常數, 約為0.0366; κe(f) 為復波數平方的隨機多項式展開系數; 隨機多項式展開基與(2)式中復聲壓的展開基為同一組，M是復波數的隨機多項式展開中, 展開基的截斷冪次, 要求M < N,在計算中通常取M = 2; E代表展開截斷冪次為M時復波數的隨機多項式展開項數.

為了求解(2)式中的復聲壓隨機多項式展開系數 { γq(f,r,z);q=0,1,···,Q?1} , 將(2)式代入到聲學拋物方程(4)式后采用蓋遼金投影法[14],在等式兩端依次正交上隨機多項式展開基的每一項, 獲得展開系數 { γq(f,r,z);q=0,1,···,Q?1} 滿足的方程組,可求得不同位置 ( r,z) 處, 復聲壓隨機多項式展開系數 γq(f,r,z) 的數值[14,28].利用(1)式中海底聲速、吸收率、比重與隨機變量的關系, 將(2)式中的替換成得到復聲壓關于海底參數的解析表達式

淺海聲傳播中傳播損失是對地聲參數較敏感且易于觀測和計算的物理量, 傳播損失定義為接收位置 ( r,z) 處聲壓與距離聲源位置1 m處的聲壓 P0幅度之比的分貝形式,表達式為

如(8)式所示, 聲傳播損失同樣可表示為隨機多項式展開形式, 其中展開系數為 { ιq(f,r,z) ; q=0,1,···,Q?1}.關于傳播損失展開系數的計算,需對(2)式中的海底參數(對應到各隨機變量)進行拉丁超立方采樣[29], 利用(2)式中復聲壓的隨機多項式展開近似結果以及(7)式中聲傳播損失的表達式, 獲取大量地聲參數-傳播損失的訓練樣本, 采用稀疏自適應最小角度回歸[30]計算得到傳播損失的隨機多項式展開系數 ιq(f,r,z).

3 基于聲傳播損失隨機多項式展開的地聲反演方法

依照以上的理論建模, 可以在一定條件下獲取聲傳播損失關于聲速、吸收率、比重在一定不確知區間內的解析表達式.本節討論利用多頻點、多收發水平距離的聲傳播損失隨機多項式展開系數矩陣, 獲取地聲參數的反演算法.

為了書寫方便, 假設接收深度z是統一的,(8)式中, 隨機多項式展開方法把頻率為 fi, 收發水平距離為 rj, 在地聲參數搜索區間內的聲傳播損失 T L(fi,rj;表示成了Q項隨機多項式展開項的疊加.對于有F個頻點、R個收發距離的聲傳播損失, 可以將(8)式寫成矩陣形式, 如(9)式.其中的展開系數 { ιq(fi,rj);q=0,1,···,Q?1 ;i=1,2,···,F;j=1,2,···,R}組成系數矩陣隨機多項式 { pq;q=0,1,···,Q?1} 記成向量; 對應的傳播損失 { TL(fi,rj) ; i=1,2,···,F;j=1,2,···,R}記為向量

其中 p0是常數項, p1p2p3分別是聲速、吸收率、比重的一次冪函數.q >3 時,是各地聲參數的高次冪項及交叉項.

若實驗區域海底真實的聲速、吸收率、比重參數為c,α,ρ, 實驗獲取的傳播損失向量為在先驗搜索區間內, 計算得到的隨機多項式展開系數矩陣為可使用線性方程組求解方法計算隨機多項展開基的數值:

4 反演仿真計算與評估

4.1 聲場隨機多項式展開的計算精度和效率

上述算法可以成功反演海底參數的前提是(2)式和(7)式中復聲壓和傳播損失的隨機多項式展開級數可以精確地近似其真值.為了方便討論(2)式中不同截斷冪次下隨機多項式展開表達聲信號的精度[31], 首先假設聲壓為平面波的情況,(14)式是聲速和吸收率隨機均勻分布的介質中, 平面波復聲壓的隨機多項式展開.

研究得到, (14)式中當聲速、吸收率的分布區間,隨機多項式展開的截斷冪次一定時, 平面波復聲壓的相對誤差正比于頻率、水平傳播距離; 誤差上限一定時, 隨機多項式展開截斷冪次越大, 可預報的頻率f越高、距離r越遠.圖1給出隨機多項式展開在不同截斷冪次下, 頻率-距離平面上復聲壓相對誤差等于1%的等值線圖.可以看出, 誤差小于1%的最遠距離r (單位為km)與最高頻率f近似滿足 r ∝f?1, 如頻率 f =50Hz , 誤差小于1%的

圖1 不同截斷冪次下, 在“頻率f-距離r”平面上隨機多項式展開平面波聲壓相對誤差1%的等值線.其中聲速的范圍是1645?1655 m/s, 吸收率的范圍是0.55?0.65 dB/λFig.1.Isolines of 1% relative error about sound pressure for plane wave expanded by the polynomial chaos in frequencyrange space.The intervals of sound speed and attenuation are 1645?1655 m/s and 0.55?0.65 dB/λ, respectively.

最遠距離r為2.9 km.

在實際波導下的聲傳播環境中, 聲壓由多個模態干涉疊加產生.以淺海Pekeris波導為例, 對隨機多項式展開的計算精度和計算效率進行分析.其中海底聲速、吸收率的范圍同上, 比重的范圍在1.4—2.0.波導的水文環境和聲源參數如表1所列.

表1 Pekeris波導的水文環境和聲源參數Table 1.Hydrological conditions of Pekeris waveguide and acoustic source parameters.

在展開截斷冪次N = 4時, Pekeris波導中地聲參數取區間中值時的距離-傳播損失曲線、地聲參數區間內隨機多項式展開近似傳播損失的距離-驗證集平均誤差曲線分別[30]如圖2和圖3所示.在傳播水平距離10 km以內, 隨機多項式展開近似傳播損失的誤差小于1%; 在趨勢上誤差仍隨著水平傳播距離r的增加而增加, 但不再單調, 如圖2和圖3中圈點部分, 傳播損失的部分極大值點位置上近似誤差出現極大值.

圖2 Pekeris波導中100 m深度不同水平距離上的聲傳播損失, 圈點為傳播損失的部分極大值點Fig.2.Acoustic transmission loss at different horizontal ranges with the depth of 100 m in Pekeris waveguide.Circled points are partial local maximum points.

圖3 Pekeris波導中100 m深度不同水平距離上隨機多項式展開近似聲傳播損失的驗證集誤差, 圈點為部分誤差極大值點Fig.3.Validation set errors of acoustic transmission loss expanded by polynomial chaos at different horizontal ranges with the depth of 100 m in Pekeris waveguide.Circled points are partial local maximum points.

在計算效率方面, 嵌入隨機多項式的聲學拋物方程的計算復雜度與算法在深度、水平方向的差分網格數和隨機多項式展開截斷冪次有關.若聲場計算中, 在深度方向上差分離散的網格數為正整數Z, 水平方向上差分離散的網格數為正整數X, 隨機多項式展開基個數為Q (對應隨機多項式展開截斷冪次為 N, Q與N的關系見(2)式), 算法的時間復雜度為

其中參數a正比于截斷冪次N, 當 0 ≤N≤10 時,a的變化范圍是 1.00

使用基于聲學拋物方程的遍歷方法, 取海底聲速、吸收率、比重取各自搜索區間不同網格點上的數值并計算對應的拷貝聲場.方法的計算復雜度與地聲參數搜索區間上的網格劃分數目有關.假設海底聲速、吸收率、比重區間上劃分的網格數目相同，記為正整數Y, 基于網格遍歷搜索的地聲反演算法的計算復雜度為

圖4為遍歷方案中地聲參數的搜索網格數Y取值為50, 隨機多項式展開截斷冪次N取不同值時遍歷搜索法與隨機多項式展開法運算復雜度之比隨拋物方程算法在深度上的離散網格數Z的變化曲線.可見, 若遍歷搜索法的參數搜索網格數Y一定, 隨機多項式展開方法的截斷冪次N和拋物方程深度方向上的離散網格數 Z 越小時, 隨機多項式展開法的計算效率優勢越明顯.在淺海、低頻環境下, 寬角拋物方程聲場計算在深度方向上的離散網格數 Z 較小, 隨機多項式展開法地聲反演的計算效率比遍歷搜索法會有極大提升.如在截斷冪次N≤6,Z≤6000時, 隨機多項式方法相比遍歷法的計算效率提高至少一個數量級.

圖4 不同隨機多項式截斷冪次N、不同深度方向差分網格數Z下, 遍歷法與隨機多項式展開法計算復雜度之比Fig.4.Calculation complexity ratio of the ergodic method and the polynomial chaos expansion method under different polynomial truncated power N and difference grid number in depth direction Z.

4.2 地聲反演仿真計算與評估

仿真采用淺海負梯度聲速剖面的聲傳播環境,海水聲速剖面如圖5所示.聲源深度50 m, 采用20—200 Hz頻段, 1/3倍頻程間隔共計12個中心頻點的數據.接收水聽器位于海底100 m深度, 收發水平距離為2000—5000 m, 均勻間隔共計61個接收水聽器.海底縱波聲速、吸收率、比重的搜索區間為

圖5 仿真計算使用的海水聲速剖面Fig.5.Sound speed profile used in the simulation case.

以嵌入隨機多項式的拋物方程模型計算展開系數矩陣, 用簡正波模型[32]計算的傳播損失作為仿真的實驗傳播損失, 以此反演海底的聲學參數.

考慮實際的實驗條件, 聲信號的信噪比條件、水文的不確知和水聽器位置的誤差都會影響反演結果的準確度.仿真驗證了在接收信號不同信噪比, 收發水平距離存在不同大小的隨機誤差時的反演結果.表2和表3為在海底聲學參數搜索區間隨機取樣5000組聲速、吸收率、比重作為真實值, 使用隨機多項式展開反演結果之誤差均值.海底聲速、吸收率、比重的反演平均誤差隨著信噪比的增加近似指數衰減, 隨著收發水平距離誤差的增加近似指數增加, 在實驗中應避免選擇信噪比低或收發水平位置存在較大誤差(如垂直陣中靠近海表的水聽器)的信號作為數據, 以提高反演精度.

表2 接收信號不同信噪比下反演誤差均值Table 2.Average inversion errors under different signal to noise ratios.

表3 收發水平距離不同誤差下反演誤差均值Table 3.Average inversion errors under different horizontal distances.

在信號存在噪聲、收發水平距離存在隨機誤差時, 存在最優的展開冪次 Nb使得反演的誤差最小.圖6給出了無噪聲、收發距離準確, 以及信噪比15 dB、收發距離存在 ± 20 m以內的隨機誤差時, 反演誤差均值與隨機多項式截斷冪次N的變化.

圖6 不同隨機多項式展開截斷冪次下, (a)海底聲速、(b)吸收率、(c)比重反演結果的誤差Fig.6.Errors of geo-acoustic inversion results for (a) sound speed, (b) attenuation and (c) ratio of density under different truncated powers of polynomial chaos.

在不考慮信號噪聲和收發位置隨機因素時, 反演結果隨著隨機多項式展開截斷冪次的增加而愈加準確; 但是噪聲和收發水平距離誤差使得反演存在最優截斷冪次 Nb.

1) 當截斷冪次 N

圖7 隨機多項式展開截斷冪次 N =1 時的海底聲速的反演結果Fig.7.Geo-acoustic inversion results of sound speed of seabed when the truncated power of polynomial chaos N =1.

2) 當截斷冪次 N >Nb時, 隨著截斷冪次N的增加, 隨機多項式展開基的個數快速增加, 系數矩陣的條件數快速增長, 方程(13)的解對傳播損失的誤差愈加敏感, 反演方法的魯棒性喪失,使得反演結果不準確.

Nb的選擇受到待反演物理量、信噪比、收發水平距離誤差大小的影響.仿真算例中, 信噪比15 dB、收發水平距離誤差為 ± 20 m時, 海底聲速、吸收率、比重的最優截斷冪次 Nb分別為 Nb=3,2,3.在反演的區間內隨機抽取聲速、吸收率、比重的三者測試樣本200組, 反演結果與95%置信區間如圖8所示.在此條件下, 海底聲速、吸收率、比重的反演平均誤差為1.34 m/s, 0.02 dB/λ, 0.03.在實際的反演過程中, 為了保證反演結果的可靠性, 可以采用后驗方法評估反演的結果, 從采用不同截斷冪次的反演結果中挑選出泛化效果最優的結果.

圖8 信噪比為15 dB, 收發水平距離誤差為 ± 20 m時,(a)模擬聲速、(b)吸收率、(c)比重的反演結果及其置信區間Fig.8.Geo-acoustic inversion results and its confidential intervals of (a) sound speed, (b) sound attenuation and (c) ratios of density between seabed and sea water when the signal to noise ratio is 15 dB and error of horizontal distance is ± 20 m.

在反演的計算時間方面, 采用遍歷法在海底參數搜索區間網格數Y = 40與采用隨機多項式展開法在不同截斷冪次N下反演所需要的計算時間如表4所列.在本算例采用的最優截斷冪次下, 隨機多項式展開法的計算用時在遍歷法的1%以內。此外, 仿真發現: 在截斷冪次 N ≤6 時, 使用隨機多項式展開方法反演地聲參數比基于代價函數搜索的遺傳算法具有效率優勢.

表4 使用遍歷法與隨機多項式展開法反演海底參數所需計算時間Table 4.Calculation time of geo-acoustic inversion using method of traversing and polynomial chaos expansion.

5 結 論

本文提出一種利用寬帶、多收發位置的傳播損失和地聲參數搜索區間內傳播損失的隨機多項式展開系數矩陣快速反演海底聲速、吸收率、比重的方法。展開系數矩陣由嵌入隨機多項式的寬角拋物方程聲場計算模型求解獲得.根據仿真實驗的結果，得到以下主要結論.

1)隨機多項式展開法在中低聲波頻率、一定的聲傳播水平距離以內、海底參數搜索區間長度適中時, 相較于遍歷搜索的地聲反演方法, 提高了反演的計算效率: 在深度上的差分網格數小于6000,隨機多項式展開截斷冪次小于6時, 計算效率至少可以提高一個數量級.在仿真算例中, 使用最優的隨機多項式展開截斷冪次時, 隨機多項式展開法反演地聲參數的時間開銷不到遍歷法的1%.

2)在高信噪比、收發相對位置誤差較小時，使用寬帶多收發位置傳播損失隨機多項式展開系數矩陣, 可準確且唯一地反演海底聲速、吸收率、比重.仿真算例中: 在信噪比大于15 dB, 收發距離誤差小于20 m時, 聲速、吸收率、比重的反演誤差小于1.34 m/s, 0.02 dB/λ, 0.03.