利用隨機多項式展開的海底聲學(xué)參數(shù)反演方法

李風(fēng)華 王翰卓

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

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

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

1 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 反演仿真計算與評估

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

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

研究得到, (14)式中當(dāng)聲速、吸收率的分布區(qū)間,隨機多項式展開的截斷冪次一定時, 平面波復(fù)聲壓的相對誤差正比于頻率、水平傳播距離; 誤差上限一定時, 隨機多項式展開截斷冪次越大, 可預(yù)報的頻率f越高、距離r越遠(yuǎn).圖1給出隨機多項式展開在不同截斷冪次下, 頻率-距離平面上復(fù)聲壓相對誤差等于1%的等值線圖.可以看出, 誤差小于1%的最遠(yuǎn)距離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.

最遠(yuǎn)距離r為2.9 km.

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

表1 Pekeris波導(dǎo)的水文環(huán)境和聲源參數(shù)Table 1.Hydrological conditions of Pekeris waveguide and acoustic source parameters.

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

圖2 Pekeris波導(dǎo)中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波導(dǎo)中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.

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

其中參數(shù)a正比于截斷冪次N, 當(dāng) 0 ≤N≤10 時,a的變化范圍是 1.00

使用基于聲學(xué)拋物方程的遍歷方法, 取海底聲速、吸收率、比重取各自搜索區(qū)間不同網(wǎng)格點上的數(shù)值并計算對應(yīng)的拷貝聲場.方法的計算復(fù)雜度與地聲參數(shù)搜索區(qū)間上的網(wǎng)格劃分?jǐn)?shù)目有關(guān).假設(shè)海底聲速、吸收率、比重區(qū)間上劃分的網(wǎng)格數(shù)目相同，記為正整數(shù)Y, 基于網(wǎng)格遍歷搜索的地聲反演算法的計算復(fù)雜度為

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

圖4 不同隨機多項式截斷冪次N、不同深度方向差分網(wǎng)格數(shù)Z下, 遍歷法與隨機多項式展開法計算復(fù)雜度之比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 地聲反演仿真計算與評估

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

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

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

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

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

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

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

圖6 不同隨機多項式展開截斷冪次下, (a)海底聲速、(b)吸收率、(c)比重反演結(jié)果的誤差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.

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

1) 當(dāng)截斷冪次 N

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

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

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

圖8 信噪比為15 dB, 收發(fā)水平距離誤差為 ± 20 m時,(a)模擬聲速、(b)吸收率、(c)比重的反演結(jié)果及其置信區(qū)間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.

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

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

5 結(jié) 論

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

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

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