李風(fēng)華 王翰卓
1) (中國(guó)科學(xué)院聲學(xué)研究所, 聲場(chǎng)聲信息國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 100190)
2) (中國(guó)科學(xué)院大學(xué)電子電氣與通信工程學(xué)院, 北京 100049)
為了提高地聲反演算法的計(jì)算效率, 探索克服地聲反演結(jié)果多值性問(wèn)題, 本文利用寬帶、多收發(fā)位置的傳播損失數(shù)據(jù)結(jié)合傳播損失在地聲參數(shù)先驗(yàn)搜索區(qū)間內(nèi)的隨機(jī)多項(xiàng)式展開(kāi)系數(shù)矩陣, 反演得到海底縱波聲速、吸收率和密度比重.使用隨機(jī)多項(xiàng)式展開(kāi)近似傳播損失時(shí), 展開(kāi)系數(shù)的自變量為聲波頻率、收發(fā)位置等參數(shù), 隨機(jī)多項(xiàng)式的自變量為表示聲速、吸收率、比重在各自搜索區(qū)間內(nèi)均勻分布的隨機(jī)變量.傳播損失的展開(kāi)系數(shù)通過(guò)嵌入隨機(jī)多項(xiàng)式的聲學(xué)寬角拋物方程結(jié)合蓋遼金投影、最小角度回歸算法計(jì)算求得.在低頻、一定聲傳播水平距離以內(nèi)和地聲參數(shù)搜索區(qū)間長(zhǎng)度適中時(shí), 使用隨機(jī)多項(xiàng)式展開(kāi)近似傳播損失的相對(duì)誤差在1%以下.仿真發(fā)現(xiàn), 在淺海環(huán)境中使用低頻、一定聲傳播水平距離以內(nèi)的傳播損失數(shù)據(jù), 在接收信號(hào)信噪比較高、聲源和水聽(tīng)器相對(duì)位置誤差較小時(shí), 選擇合適的隨機(jī)多項(xiàng)式展開(kāi)截?cái)鄡绱慰奢^準(zhǔn)確地反演海底聲速、吸收率和密度比重, 且計(jì)算效率比網(wǎng)格遍歷搜索方法提高一個(gè)數(shù)量級(jí)以上.
海洋聲學(xué)中的逆問(wèn)題是利用海水中聲波所攜帶的聲源位置及海洋信道環(huán)境特性信息, 對(duì)聲源位置和介質(zhì)參數(shù)進(jìn)行求解.地聲反演問(wèn)題是快速獲取局部海域等效海底聲學(xué)參數(shù)的方法, 是水聲逆問(wèn)題的重要研究方向[1].地聲反演過(guò)程需要進(jìn)行以下幾個(gè)方面的內(nèi)容: 海底環(huán)境模型選取, 待反演參數(shù)搜索范圍的確定, 拷貝聲場(chǎng)計(jì)算, 匹配物理量及代價(jià)函數(shù)的構(gòu)造, 最優(yōu)化算法的選取, 反演結(jié)果的評(píng)價(jià).為求解待反演參數(shù)的全局最優(yōu)值, 需要采用遍歷或最優(yōu)化算法在待反演參數(shù)空間內(nèi)反復(fù)選取待反演參數(shù)向量、計(jì)算其對(duì)應(yīng)的拷貝聲場(chǎng)以及代價(jià)函數(shù).傳統(tǒng)基于遍歷搜索方法所需計(jì)算量大, 計(jì)算時(shí)間長(zhǎng);代價(jià)函數(shù)對(duì)部分待反演參數(shù)的不敏感以及待反演參數(shù)在代價(jià)函數(shù)中的耦合效果使得代價(jià)函數(shù)存在多解性問(wèn)題[1].提高反演速率以及克服多解性是地聲反演中需要解決的問(wèn)題.
為了減少求解反演參數(shù)空間拷貝聲場(chǎng)的計(jì)算規(guī)模, Gerstoft[2,3]和Dosso[4,5]分別提出了基于遺傳算法和模擬退火的貝葉斯推斷法, 該方法可快速求解作為代價(jià)函數(shù)的后驗(yàn)概率密度; Sambridge[6]提出了鄰域插值近似方法, 該方法利用集合信息指導(dǎo)參數(shù)空間的重采樣, 從而減少了計(jì)算量.
為了解決地聲反演中的多值性問(wèn)題, Holland和Osler[7]提出了組合反演方法, 將空間-時(shí)間域和空間-頻率域的數(shù)據(jù)結(jié)合起來(lái), 采用多種不同的地聲模型擬合同一數(shù)據(jù)集.李整林等[8?13]在綜合分析簡(jiǎn)正波的頻散特性、海底反射系數(shù)、傳播損失等聲場(chǎng)特性的基礎(chǔ)上, 應(yīng)用匹配場(chǎng)處理器、自適應(yīng)時(shí)頻分析算法和并行遺傳算法提出多物理量聯(lián)合地聲反演方法.
隨機(jī)多項(xiàng)式展開(kāi)法[14]是一種表示方程中的多維隨機(jī)參數(shù)在方程解中非線性傳遞函數(shù)關(guān)系的方法.近些年引入到海水聲速起伏[15?23]、海底參數(shù)不確知[22?25]、聲源和接收點(diǎn)深度不確知[25,26]的聲場(chǎng)建模中.隨機(jī)多項(xiàng)式展開(kāi)方法將頻域復(fù)聲壓信號(hào)表示成以隨機(jī)多項(xiàng)式為基底的級(jí)數(shù)展開(kāi)形式, 隨機(jī)多項(xiàng)式的自變量即為描述不確知(或隨機(jī))環(huán)境的隨機(jī)變量.在地聲反演問(wèn)題中, 假設(shè)待反演參數(shù)為在搜索區(qū)間內(nèi)均勻分布的隨機(jī)變量[24], 通過(guò)隨機(jī)多項(xiàng)式展開(kāi)可以得到聲場(chǎng)關(guān)于待反演海底參數(shù)的解析表達(dá)式.該方法的計(jì)算效率較遍歷方法高, 可以作為計(jì)算拷貝聲場(chǎng)的“代理”模型.
本文利用寬帶多收發(fā)位置的聲傳播損失數(shù)據(jù),結(jié)合地聲參數(shù)搜索區(qū)間內(nèi)傳播損失函數(shù)的隨機(jī)多項(xiàng)式展開(kāi)系數(shù)矩陣, 求解隨機(jī)多項(xiàng)式基底的數(shù)值,利用基底中的一次冪項(xiàng), 唯一地反演海底聲學(xué)參數(shù).該方法比傳統(tǒng)遍歷搜索法計(jì)算速度快, 且對(duì)信號(hào)噪聲和收發(fā)距離隨機(jī)誤差有一定的穩(wěn)健性.
在頻率為f, 二維柱坐標(biāo)水平和深度位置(r,z)處, 海底參數(shù)在各自區(qū)間內(nèi)取值對(duì)應(yīng)的向量為時(shí), 頻域復(fù)聲壓可以用隨機(jī)多項(xiàng)式展開(kāi)方法表示為
其中 k0是參考波數(shù);{γq(f,r,z);q=0,1,···,Q?1}為隨機(jī)多項(xiàng)式展開(kāi)的系數(shù), 是頻率f和空間位置(r,z) 的函數(shù);是隨機(jī)多項(xiàng)式展開(kāi)基, 是隨機(jī)變量的函數(shù); N是復(fù)聲壓的隨機(jī)多項(xiàng)式展開(kāi)中展開(kāi)基的截?cái)鄡绱? Q代表展開(kāi)截?cái)鄡绱螢镹時(shí), 復(fù)聲壓的隨機(jī)多項(xiàng)式展開(kāi)項(xiàng)數(shù).在均方意義下, 在截?cái)鄡绱?N →∞ 時(shí), 頻率復(fù)聲壓的隨機(jī)多項(xiàng)式展開(kāi)近似結(jié)果可以收斂到真實(shí)值, 誤差隨著截?cái)鄡绱蜰的增加指數(shù)減小[14].隨機(jī)多項(xiàng)式基底有勒讓德多項(xiàng)式的函數(shù)形式, 滿足加權(quán)正交歸一化條件:
其中 μ 為常數(shù), 約為0.0366; κe(f) 為復(fù)波數(shù)平方的隨機(jī)多項(xiàng)式展開(kāi)系數(shù); 隨機(jī)多項(xiàng)式展開(kāi)基與(2)式中復(fù)聲壓的展開(kāi)基為同一組,M是復(fù)波數(shù)的隨機(jī)多項(xiàng)式展開(kāi)中, 展開(kāi)基的截?cái)鄡绱? 要求M < N,在計(jì)算中通常取M = 2; E代表展開(kāi)截?cái)鄡绱螢镸時(shí)復(fù)波數(shù)的隨機(jī)多項(xiàng)式展開(kāi)項(xiàng)數(shù).
為了求解(2)式中的復(fù)聲壓隨機(jī)多項(xiàng)式展開(kāi)系數(shù) { γq(f,r,z);q=0,1,···,Q?1} , 將(2)式代入到聲學(xué)拋物方程(4)式后采用蓋遼金投影法[14],在等式兩端依次正交上隨機(jī)多項(xiàng)式展開(kāi)基的每一項(xiàng), 獲得展開(kāi)系數(shù) { γq(f,r,z);q=0,1,···,Q?1} 滿足的方程組,可求得不同位置 ( r,z) 處, 復(fù)聲壓隨機(jī)多項(xiàng)式展開(kāi)系數(shù) γq(f,r,z) 的數(shù)值[14,28].利用(1)式中海底聲速、吸收率、比重與隨機(jī)變量的關(guān)系, 將(2)式中的替換成得到復(fù)聲壓關(guān)于海底參數(shù)的解析表達(dá)式
淺海聲傳播中傳播損失是對(duì)地聲參數(shù)較敏感且易于觀測(cè)和計(jì)算的物理量, 傳播損失定義為接收位置 ( r,z) 處聲壓與距離聲源位置1 m處的聲壓 P0幅度之比的分貝形式,表達(dá)式為
如(8)式所示, 聲傳播損失同樣可表示為隨機(jī)多項(xiàng)式展開(kāi)形式, 其中展開(kāi)系數(shù)為 { ιq(f,r,z) ; q=0,1,···,Q?1}.關(guān)于傳播損失展開(kāi)系數(shù)的計(jì)算,需對(duì)(2)式中的海底參數(shù)(對(duì)應(yīng)到各隨機(jī)變量)進(jìn)行拉丁超立方采樣[29], 利用(2)式中復(fù)聲壓的隨機(jī)多項(xiàng)式展開(kāi)近似結(jié)果以及(7)式中聲傳播損失的表達(dá)式, 獲取大量地聲參數(shù)-傳播損失的訓(xùn)練樣本, 采用稀疏自適應(yīng)最小角度回歸[30]計(jì)算得到傳播損失的隨機(jī)多項(xiàng)式展開(kāi)系數(shù) ιq(f,r,z).
依照以上的理論建模, 可以在一定條件下獲取聲傳播損失關(guān)于聲速、吸收率、比重在一定不確知區(qū)間內(nèi)的解析表達(dá)式.本節(jié)討論利用多頻點(diǎn)、多收發(fā)水平距離的聲傳播損失隨機(jī)多項(xiàng)式展開(kāi)系數(shù)矩陣, 獲取地聲參數(shù)的反演算法.
為了書寫方便, 假設(shè)接收深度z是統(tǒng)一的,(8)式中, 隨機(jī)多項(xiàng)式展開(kāi)方法把頻率為 fi, 收發(fā)水平距離為 rj, 在地聲參數(shù)搜索區(qū)間內(nèi)的聲傳播損失 T L(fi,rj;表示成了Q項(xiàng)隨機(jī)多項(xiàng)式展開(kāi)項(xiàng)的疊加.對(duì)于有F個(gè)頻點(diǎn)、R個(gè)收發(fā)距離的聲傳播損失, 可以將(8)式寫成矩陣形式, 如(9)式.其中的展開(kāi)系數(shù) { ιq(fi,rj);q=0,1,···,Q?1 ;i=1,2,···,F;j=1,2,···,R}組成系數(shù)矩陣隨機(jī)多項(xiàng)式 { pq;q=0,1,···,Q?1} 記成向量; 對(duì)應(yīng)的傳播損失 { TL(fi,rj) ; i=1,2,···,F;j=1,2,···,R}記為向量
其中 p0是常數(shù)項(xiàng), p1p2p3分別是聲速、吸收率、比重的一次冪函數(shù).q >3 時(shí),是各地聲參數(shù)的高次冪項(xiàng)及交叉項(xiàng).
若實(shí)驗(yàn)區(qū)域海底真實(shí)的聲速、吸收率、比重參數(shù)為c,α,ρ, 實(shí)驗(yàn)獲取的傳播損失向量為在先驗(yàn)搜索區(qū)間內(nèi), 計(jì)算得到的隨機(jī)多項(xiàng)式展開(kāi)系數(shù)矩陣為可使用線性方程組求解方法計(jì)算隨機(jī)多項(xiàng)展開(kāi)基的數(shù)值:
上述算法可以成功反演海底參數(shù)的前提是(2)式和(7)式中復(fù)聲壓和傳播損失的隨機(jī)多項(xiàng)式展開(kāi)級(jí)數(shù)可以精確地近似其真值.為了方便討論(2)式中不同截?cái)鄡绱蜗码S機(jī)多項(xiàng)式展開(kāi)表達(dá)聲信號(hào)的精度[31], 首先假設(shè)聲壓為平面波的情況,(14)式是聲速和吸收率隨機(jī)均勻分布的介質(zhì)中, 平面波復(fù)聲壓的隨機(jī)多項(xiàng)式展開(kāi).
研究得到, (14)式中當(dāng)聲速、吸收率的分布區(qū)間,隨機(jī)多項(xiàng)式展開(kāi)的截?cái)鄡绱我欢〞r(shí), 平面波復(fù)聲壓的相對(duì)誤差正比于頻率、水平傳播距離; 誤差上限一定時(shí), 隨機(jī)多項(xiàng)式展開(kāi)截?cái)鄡绱卧酱? 可預(yù)報(bào)的頻率f越高、距離r越遠(yuǎn).圖1給出隨機(jī)多項(xiàng)式展開(kāi)在不同截?cái)鄡绱蜗? 頻率-距離平面上復(fù)聲壓相對(duì)誤差等于1%的等值線圖.可以看出, 誤差小于1%的最遠(yuǎn)距離r (單位為km)與最高頻率f近似滿足 r ∝f?1, 如頻率 f =50Hz , 誤差小于1%的
圖1 不同截?cái)鄡绱蜗? 在“頻率f-距離r”平面上隨機(jī)多項(xiàng)式展開(kāi)平面波聲壓相對(duì)誤差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.
在實(shí)際波導(dǎo)下的聲傳播環(huán)境中, 聲壓由多個(gè)模態(tài)干涉疊加產(chǎn)生.以淺海Pekeris波導(dǎo)為例, 對(duì)隨機(jī)多項(xiàng)式展開(kāi)的計(jì)算精度和計(jì)算效率進(jìn)行分析.其中海底聲速、吸收率的范圍同上, 比重的范圍在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.
在展開(kāi)截?cái)鄡绱蜰 = 4時(shí), Pekeris波導(dǎo)中地聲參數(shù)取區(qū)間中值時(shí)的距離-傳播損失曲線、地聲參數(shù)區(qū)間內(nèi)隨機(jī)多項(xiàng)式展開(kāi)近似傳播損失的距離-驗(yàn)證集平均誤差曲線分別[30]如圖2和圖3所示.在傳播水平距離10 km以內(nèi), 隨機(jī)多項(xiàng)式展開(kāi)近似傳播損失的誤差小于1%; 在趨勢(shì)上誤差仍隨著水平傳播距離r的增加而增加, 但不再單調(diào), 如圖2和圖3中圈點(diǎn)部分, 傳播損失的部分極大值點(diǎn)位置上近似誤差出現(xiàn)極大值.
圖2 Pekeris波導(dǎo)中100 m深度不同水平距離上的聲傳播損失, 圈點(diǎn)為傳播損失的部分極大值點(diǎn)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深度不同水平距離上隨機(jī)多項(xiàng)式展開(kāi)近似聲傳播損失的驗(yàn)證集誤差, 圈點(diǎn)為部分誤差極大值點(diǎn)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.
在計(jì)算效率方面, 嵌入隨機(jī)多項(xiàng)式的聲學(xué)拋物方程的計(jì)算復(fù)雜度與算法在深度、水平方向的差分網(wǎng)格數(shù)和隨機(jī)多項(xiàng)式展開(kāi)截?cái)鄡绱斡嘘P(guān).若聲場(chǎng)計(jì)算中, 在深度方向上差分離散的網(wǎng)格數(shù)為正整數(shù)Z, 水平方向上差分離散的網(wǎng)格數(shù)為正整數(shù)X, 隨機(jī)多項(xiàng)式展開(kāi)基個(gè)數(shù)為Q (對(duì)應(yīng)隨機(jī)多項(xiàng)式展開(kāi)截?cái)鄡绱螢?N, Q與N的關(guān)系見(jiàn)(2)式), 算法的時(shí)間復(fù)雜度為