隨機多項式展開多特征向量約束-匹配場聲源定位算法*

王翰卓 李風華

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

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

0 引言

自適應匹配場聲源定位算法的加權向量由拷貝場向量和實驗數據互譜矩陣聯合得到。與傳統線性匹配器相比，自適應匹配場定位方法具有對旁瓣更好的抑制能力。但自適應聲源定位算法在實際應用中難以達到理想性能，主要原因有：(1)算法對水文環境參數失配和水聽器位置誤差比較敏感；(2)目標運動和水文環境的非平穩性導致的相干信號拍數的不足，難以獲得真實的數據互譜矩陣；(3)水面強干擾以及定位二維網格空間劃分過于稀疏將降低定位的準確率[1]。環境的失配是主導因素之一，因而，尋找耐受環境失配的聲源定位算法是匹配場定位研究的重點。對此，研究者們對應開發了環境寬容匹配場處理器[2-5]、最優不確定場處理器[6-7]、降階自適應匹配場處理器[8]、基于扇區特征向量約束的匹配場處理器[9-10]等。以上耐受環境失配的聲源定位算法通常需要使用蒙特卡洛統計方法對隨機或不確知的聲傳播環境參數進行反復抽樣后計算拷貝聲場。如上述環境擾動約束的匹配場處理器需要對隨機或不確知環境下拷貝場復聲壓的互譜矩陣進行估計，最優不確定場處理器在求取關于聲源位置的貝葉斯邊緣后驗概率密度時需要進行蒙特卡洛積分。針對二維空間聲源定位，在空間網格劃分完成后，依次假設聲源位于其中的某個網格點上求取匹配器的輸出。當定位網格點劃分較稠密時，算法是比較耗時的。

為了克服蒙特卡洛統計方法計算量大、耗時高的問題，學者們提出了系列蒙特卡洛統計的代理方法，其中有代表性的為聲場位移法[11-13]和隨機多項式展開法[14]。其中隨機多項式展開法由于其通用性較廣泛受到部分研究者的青睞。隨機多項式展開法將參數為隨機變量的微分方程的解表示成以相互正交的隨機多項式為基底的級數形式，隨機多項式基底是方程中隨機參數的函數。隨機多項式展開方法提供了一種由輸入端(方程參數)到輸出端(方程解)的解析表達式，方便了有關方程輸入參數的隨機性在輸出量中傳遞性的研究。對于方程解中隨機多項式展開系數的求解，有基于蓋遼金投影的嵌入型方法[14-15]和基于聯合最小角度回歸[16]、概率配點[17]、隨機響應面[18-20]等一系列非嵌入型方法。近些年來，隨機多項式展開方法被應用到環境參數隨機或不確知時不確定聲場的計算中：文獻[21-24]使用Karhunen Loève展開[25]描述了水體聲速水平方向的擾動量，使用嵌入隨機多項式的拋物方程計算了聲強均值、方差和復聲壓的水平互相關等統計量。通過與蒙特卡洛統計結果的比較，證實了該方法的準確性和高效性。相同近似誤差下，聲速起伏強度越強、復聲壓統計矩階次越高時，復聲壓需要的展開截斷冪次越高；文獻[26]結合微分形式的單向耦合簡正波聲場計算模型，將各號模態的相位表示為隨機多項式展開的形式。相比于直接將模態復幅度表示為隨機多項式展開的方法，其使用較低的展開截斷冪次即可獲得準確的聲場統計結果。文獻[27-31]在淺海聲源深度、潮汐、初始溫躍層結構、沉積層聲速等存在多種隨機或不確定因素時，使用拉丁超立方抽樣采集訓練樣本，以多元線性回歸求解了復聲壓的隨機多項式展開系數。在此基礎上，計算了聲強均值、方差和概率密度，并進行了敏感性分析。結果證實了環境存在多類型隨機和不確知參數時隨機多項式展開法的高效性和準確性；文獻[32-34]使用海洋動力學模型模擬了海水聲速的時空變化，引入深度上的經驗正交函數降低海水聲速不確定參數維度，使用隨機多項式展開方法預報了聲傳播損失的概率密度并進行了實驗驗證。

因此，利用隨機多項式展開表示環境隨機或不確知時的不確定聲場，可以快速計算拷貝場復聲壓的互譜矩陣，從而提高環境寬容匹配場聲源定位算法的計算效率。如文獻[35]利用隨機多項式展開進行了不確知海洋中環境寬容的貝葉斯波束形成，在實驗中實現了淺海不確定場中聲源位置的定位。本文將以環境寬容的匹配場聲源定位算法為例，討論受內波、湍流等因素影響，淺海水體聲速存在起伏的情況下，使用隨機多項式展開表示拷貝場隨機復聲壓、并估計其互譜矩陣時定位算法的性能表現。仿真討論了輸入為單拍、不同信噪比下算法的定位準確率和輸出峰均比，比較了蒙特卡洛統計和隨機多項式展開估計拷貝場聲壓互譜矩陣的計算耗時。結果表明，在淺海海水深度-聲速剖面為負梯度且聲傳播路徑上的二維聲速存在起伏的環境中，在低頻、水平聲傳播一定距離內，隨機多項式展開法可以準確估計垂直陣接收聲壓的互譜矩陣，且計算效率較蒙特卡洛統計方法提高十倍以上；多特征向量約束的最小方差匹配場處理器可以實現比線性匹配器和最小方差無失真匹配器更高的定位準確率和輸出模糊度峰均比。

1 隨機多項式展開多特征向量約束匹配器

假設二維柱坐標系中，r是聲傳播水平方向的坐標，z是海水深度方向的坐標。c(r,z,ξ)是線性內波擾動下聲傳播路徑上的海水聲速，其中矢量ξ= [ξ1,ξ2,··· ,ξL]，ξ中的元素為相互獨立且服從標準正態分布的隨機變量。二維聲傳播路徑上的海水聲速表示為平均聲速ˉc(z)與擾動量δc(r,z,ξ)之和：

在深度上只保留第一階經驗正交函數[36]時，采用Karhunen Loève展開[25]表示聲速擾動量為

其中， eof(z)為海水聲速的主經驗正交函數；λ1,λ2,··· ,λL為經驗正交函數展開系數的Karhunen Loève 展開特征值，L為展開項數；φ1(r),φ2(r),··· ,φL(r)為展開基。

本文采用了結合寬角拋物方程的嵌入型方法[21,24]計算頻域復聲壓的隨機多項式展開系數。假設某頻點空間位置(r,z)處復聲壓P(r,z,ξ)的隨機多項式展開為

其中，k0為參考波數；U(r,z,ξ)為復聲壓的包絡；{pq(ξ)q=0,1,···,Q-1}為展開基底函數集合，函數自變量為描述海水聲速起伏的隨機變量；γ0(r,z),γ1(r,z),··· ,γQ-1(r,z)為展開系數，其自變量為位置坐標；Q為隨機多項式的展開項數，L為獨立隨機變量的個數，N為展開的截斷冪次，三者關系為Q=(L+N)!/L!/N!。

由于隨機變量為高斯型隨機變量，展開基對應為埃爾米特多項式，當隨機變量的個數L=2、展開截斷冪次N= 2 時，埃爾米特隨機多項式的函數形式如表1所示。

表1 埃爾米特隨機多項式的函數形式Table 1 Hermite polynomials

隨機多項式展開基在概率空間上相互正交，

其中，ρ(ξ)為高斯型獨立隨機變量的聯合概率密度函數，δqq′為Kronecker函數。

對于隨機多項式展開基的具體構造方法及其性質，讀者可參考文獻[14]。

在寬角聲學拋物方程中， 復聲壓的包絡U(r,z,ξ)滿足方程

其中，算子X在聲速存在隨機擾動時近似為

將式(6)代入到聲學寬角拋物方程(5)中，并使用Padé近似[37]表示式(5)中的寬角拋物方程算子。

之后，將復聲壓包絡的隨機多項式展開級數式(3)代入式(5)中，使用蓋遼金投影，對式(5)式兩端依次乘上p0(ξ),p1(ξ),··· ,pQ-1(ξ)并作期望。投影保證了復聲壓包絡的展開誤差正交于隨機多項式基底張成的函數空間。計算結果得到位置(r,z)處， 復聲壓包絡展開系數γ0(r,z),γ1(r,z),··· ,γQ-1(r,z)的Q個耦合微分方程組，結合RAM模型的數值方法可對之求解[37-38]。

假設聲源的真實位置為(rs,zs)，固定位置的水聽器陣列接收聲信號某頻點的復聲壓列向量為P(rs,zs,ξ)，復聲壓互譜矩陣Re(rs,zs)定義為

等式最右端為使用M拍實驗數據對聲壓互譜矩陣的估計，Pm(rs,zs,ξ)代表第m拍的拷貝場復聲壓向量，其中上標*代表共軛轉置。

其中，q1,q2,··· ,qD代表各約束點的期望響應。問題轉化為優化式(10)中的函數，其中c1,c2,··· ,cD為拉格朗日乘子。

定義期望響應向量q= [q1,q2,··· ,qD]T，求解式(10)得到匹配器的加權向量為

對應匹配器輸出歸一化后的平面模糊度函數為

期望響應向量q的選擇決定了匹配器的定位準確率和輸出增益。通常期望響應向量的選擇有最大化平均的空間白噪聲增益和最大化最小的空間白噪聲增益兩類[5]。假設靜態加權向量為wq，前者使得匹配器輸出增益

最大，對應wq=h1，對應q=[1,0,··· ,0]T；后者選擇靜態加權向量以獲得最差最優增益

且wq*wq= 1。式(14)意味著在聲速隨機變量ξ的所有的抽樣中，選擇wq使其最小的輸出增益最大化。該方法沒有顯式形式的靜態加權向量解。

由于對匹配器的定位準確率、輸出峰均比的評價定義在聲速隨機變量集平均意義上，且為了最大化發揮隨機多項式展開方法對匹配場算法效率的提升，本文研究選擇最大化平均的空間白噪聲增益下的期望響應，即q=[1,0,··· ,0]T。

上述算法需要對二維空間劃分的每個網格點上對拷貝場復聲壓互譜矩陣進行估計。在一定條件下，隨機多項式展開方法可以替代蒙特卡洛統計方法，更高效地估計互譜矩陣。蒙特卡洛統計方法如式(7)，需要多次抽樣生成聲傳播路徑上隨機的海水聲速，并計算陣列位置處復聲壓的互譜矩陣。式(7)中，當抽樣次數M足夠多時，聲壓互譜矩陣的統計估計值將接近真實值，估計誤差的收斂速度正比于M-1/2，即抽樣次數的-1/2 次冪；采用隨機多項式展開方法，若第i個水聽器位置處拷貝場隨機復聲壓為Pi，第j個水聽器位置處拷貝場隨機復聲壓為Pj，利用式(4)中隨機多項式的正交性，兩處復聲壓的互譜即拷貝場互譜矩陣第i行第j列元素數值為隨機多項式展開法對復聲壓互譜矩陣的估計誤差隨展開截斷冪次N的增加指數衰減。若式(5)聲學拋物方程的計算在深度上的網格點數為Z，蒙特卡洛統計方法和隨機多項式展開方法的時間復雜度分別為O(Z×M)和O(Z3+Q3)。在淺海、低頻和一定聲傳播距離內，隨機多項式展開方法較蒙特卡洛統計方法存在一個數量級的效率優勢[24]。

2 數值模擬及結果

為了驗證多特征向量約束的匹配場聲源定位算法中，隨機多項式展開方法對蒙特卡洛統計方法在計算拷貝場互譜矩陣上的效率優勢以及算法性能，采用南中國海北部聲傳播起伏實驗SWSF2015[39-40]中的觀測水文數據作為仿真起伏環境。實驗期間內潮波、線性隨機內波和孤立子內波導致了聲傳播路徑上海水聲速的劇烈起伏[39-40]，海水平均聲速剖面和聲速均方差剖面見圖1、圖2。聲源距離海底2 m，距離海表約107 m，坐底發射10 s脈寬、中心頻率200 Hz、帶寬50 Hz的線性調頻脈沖，收發水平距離為14.73 km。信號由距離海表23～77 m 近似等間隔垂直陣列的16 個水聽器接收。匹配場聲源定位中，深度上的搜索范圍為0～109 m，水平距離上的搜索范圍為0～30 km，深度、距離上網格跨度分別為2 m、10 m。

圖2 統計均方差聲速剖面Fig.2 The depth-dependent standard deviations of the sound speeds

如式(1)～(2)，本文采用深度上的第一階經驗正交函數表示海水聲速起伏在深度方向上的分布，并使用Karhunen Loève 展開[25]近似不同水平位置處隨機系數。其中，隨機變量ξ的維度L=10，以保證各深度位置處模擬近似聲速的均方差近似等于實際聲速均方差。以式(5)中嵌入隨機多項式展開的聲學拋物方程計算聲源位于不同搜索網格點時接收拷貝場復聲壓的隨機多項式展開。取展開截斷冪次N= 2，此時展開項數目Q= 66。以樣本數目M充分大時的蒙特卡洛統計方法估計的復聲壓互譜矩陣為參考，驗證上述展開參數設置下，式(15)中隨機多項式展開方法對互譜矩陣估計的準確性。圖3為聲源深度=107 m，聲傳播水平距離= 15 km 時，靠近海面最近的23 m 第一個水聽器深度的復聲壓與0～109 m 深度復聲壓的互譜函數E[P(z0=23 m)P*(z)] 0 m ≤z≤109 m(E代表數學期望)：實線代表隨機多項式展開方法的計算結果，圓點代表接收垂直陣陣元深度上蒙特卡洛統計結果。可知，隨機多項式展開法可較為準確地估計復聲壓在深度方向上的互譜函數。

圖3 深度方向上復聲壓的互譜函數Fig.3 The correlation function of the complex pressures

基于式(1)～(2)，多次抽樣海水聲速并計算對應的垂直陣列接收的復聲壓信號，以此為模擬的測量場，分別采用線性匹配器(Bartlett)[41]，對角加載的最小方差無失真匹配器(MV)[42-43](約束白噪聲增益小于6 dB)和基于蒙特卡洛統計(MV-EPC-MC)和隨機多項式展開(MV-EPC-PC)的多特征向量約束最小方差匹配器模擬定位聲源位置。比較4 種匹配場聲源定位算法在不同輸入信噪比下的定位準確率和輸出模糊度的峰均比，同時比較基于蒙特卡洛統計方法和隨機多項式展開方法的多特征向量約束最小方差匹配器的計算耗時。

線性匹配器和最小方差無失真匹配器輸出模糊度平面分別為式(16)和式(17)，

定位準確率中判定為定位正確的“窗口”定義為與真實聲源位置誤差在深度方向上為±5 m，水平方向上為±500 m，在本算例中為深度落在102～109 m 且水平距離落在14.23～15.23 km。匹配器輸出的峰均比(Peak to background ratio,PBR)[44]定義為

其中，Zmax和Zmean代表各類算法輸出模糊度的峰值和均值。

四種匹配器在單水聽器的輸入信噪比為-10～60 dB 時，定位準確率和峰均比的結果如圖4、圖5 所示。圖6 顯示了在定位性能上，由于海水聲速的起伏，線性匹配器與白噪聲增益約束的對角加載最小方差無失真匹配器在定位準確率上最高只能達到29%，與之相比，高信噪比下基于蒙特卡洛統計的多特征向量約束最小方差匹配場處理器在定位準確率和輸出峰均比上明顯優于線性匹配器和最小方差無失真匹配器。在輸入信噪比在25 dB 以上時，定位準確率在84%，輸出峰均比在42 dB。基于隨機多項式展開方法的多特征向量約束最小方差匹配場處理器在性能上略遜于前者，定位準確率和輸出峰均比達到78% 和38 dB。式(2)中近似海水聲速起伏的Karhunen Loève 展開截斷和式(3)中近似隨機復聲壓的隨機多項式展開冪次截斷引起的拷貝場復聲壓互譜矩陣的估計誤差是導致定位算法性能下降的原因；在計算時間上，如表2 所示，使用多特征向量約束的匹配場算法中，在估計拷貝場互譜矩陣時，隨機多項式展開方法的時間開銷不到蒙特卡洛統計方法的7%。

圖4 不同信噪比下Bartlett、MV、MV-EPC-MC 以及MV-EPC-PC 四種匹配器的聲源定位準確率Fig.4 Probability of correct localization vs SNR per snapshot as a function of four matched filed processors: Bartlett,MV,MV-EPC-MC,and MV-EPC-PC

圖5 不同信噪比下Bartlett、MV、MV-EPC-MC 以及MV-EPC-PC 四種匹配器輸出模糊度的峰均比Fig.5 PBRs vs SNR per snapshot as a function of four matched filed processors: Bartlett,MV,MVEPC-MC,and MV-EPC-PC

圖6 基于隨機多項式展開的多特征向量約束最小方差匹配器(MV-EPC-PC)的輸出模糊度Fig.6 MV-EPC-PC normalized ambiguity surface

表2 多特征向量約束匹配場聲源定位算法的時間開銷Table 2 The time cost of the MV-EPC source localization algorithms

3 結論

本文提出了一種基于復聲壓隨機多項式展開的多特征向量約束不確定場聲源定位方法。通過將隨機聲壓表示成隨機多項式展開的形式，提高了估計拷貝場復聲壓互譜矩陣的計算效率。在淺海、低頻、海水聲速存在起伏時，多特征向量約束的匹配場聲源定位算法在定位性能(定位準確率和輸出峰均比)上優于線性匹配器和最小方差無失真匹配器，算法在估計拷貝場復聲壓互譜矩陣時，與蒙特卡洛統計方法相比，使用隨機多項式展開方法會使得算法在定位性能略遜于前者，但計算效率比之提升一個數量級。