相關測速聲吶模型參數矩陣化求解方法

崔曉康，薛鳳杰，趙二亮，王長紅

(1.中國科學院聲學研究所海洋聲學技術中心，北京 100190；2.中國科學院大學，北京 100190；3.北京市海洋聲學裝備工程技術研究中心，北京 100190)

0 引 言

隨著對海洋資源探索的不斷深入，潛航設備對測速導航功能的需求和設備所能提供的測速精度都在不斷提升。由于傳播問題和存在累積誤差，全球定位系統和慣性導航都不能滿足水下長時間的使用。目前主要的水下測速導航設備是測速聲吶。測速聲吶包含多普勒測速聲吶和相關測速聲吶兩種，按照不同的原理和功能進行劃分，可對應4種不同類型的設備：分別為聲相關計程儀(Acoustic Corre‐lation Log,ACL)、聲相關流速剖面儀(Acoustic Cor‐relation Current Profiler,ACCP)、聲多普勒計程儀(Acoustic Doppler Log,ADL)和聲多普勒流速剖面儀(Acoustic Doppler Current Profiler,ADCP)。多普勒測速聲吶主要依靠發射聲波和散射體之間的多普勒效應來解算載體速度信息[1]；相關測速聲吶依靠多個接收陣元接收信號的延時相關性來測量載體速度信息[2]。相對于多普勒測速聲吶相關測速聲吶具有以下兩點優勢：(1)相同尺寸的換能器，發射聲波頻率更低，作用距離更遠；(2)聲波垂直向下發射，受載體姿態影響較小。但是目前相關測速聲吶相對于多普勒測速聲吶測速精度較低。一直以來研究人員都在為提升測速精度不斷努力。馮雷博士提出了簡化的理論相關函數模型，使得參數解算更加簡便，便于實時運算[3]；蒯多杰根據修正后的理論相關函數值對最小二乘的目標函數計算公式進行了自適應加權[4]。易卉琴等[5]基于真實數據的相關模型將理論相關函數添加了方向幅值調節因子和旋轉因子，使得理論相關函數俯視圖形態由圓變為橢圓，更符合實際數據的相關形態。王映春等考慮了陣元坐標誤差給測速結果帶來的影響，利用相關延時為0 的數據和理論相關函數來估計陣元聲中心，進而提高測速精度[6]。文獻[7-8]中分析了載體縱橫搖對設備測速的影響，并提出了一種目標函數添加權重調節項來提升速度估計精度的方法。本文基于局部最小二乘算法，利用矩陣推導直接求解出目標的解析解，具有無須迭代和運算復雜度低的優點。

1 基本工作原理

波形不變性是聲相關技術的基本工作原理。波形不變性原理在常見的流體散射和海底散射等各種聲波散射過程中都適用[9-10]。

本文以最簡單的一發兩收的聲吶陣型來闡述相關測速聲吶的測速原理。相關測速聲吶采用的是寬波束發射換能器，聲波以一定的波束角垂直向下發射，以一定間隔連續發射兩次信號，對應的發射聲波海底散射區域圖如圖1 所示(圖1 中位置1 3 個陣元從右往左分別是發射陣元T、接收陣元R1和R2)。

圖1 聲波在海底散射區域的示意圖Fig.1 Schematic diagram of the sound scattering area on the sea floor

當以時間間隔τ發射兩次聲波，經過相似的傳播區域后以時間間隔τ被兩個接收陣元接收，此時兩個陣元接收的信號具有極大的相似性。根據幾何知識可以得出式(1)所示的關系式：

式(1)表明，載體速度V可以通過接收陣元R1和R2構成的矢量d12和兩次發射信號的時間間隔τ來確定。

波形不變性原理為載體測速提供了一種求解思路。該原理主要涉及載體相對于散射體的速度、接收陣元構成的空間相關矢量和發射信號相關延時3個參數。通過式(1)將求速度問題轉換成了求矢量d12和相關延時τ的問題，矢量和相關延時兩個參數如何獲得取決于相關測速聲吶的設計。有3種比較通用的方法[3]：

(1) 時間相關方法。固定矢量d12，通過相關函數搜索相關延時，然后得到速度信息。

(2) 空間相關方法。固定相關延時τ時，在較多的空間矢量中搜索出適合的矢量d12，進而得到速度信息。

(3) 時空相關方法。通過多種相關延時與陣型設計的多個矢量d12相結合，進而得到速度信息。

通常采用的速度解算方法是空間相關方法，即利用有限的陣元進行特定的布局，利用回波信號構建空間相關矩陣，借助理論相關函數來擬合數據相關矩陣的數據來獲取目標函數，最終求取極值來獲取水平速度信息。ACL數據處理流程如圖2所示。

圖2 ACL數據處理流程圖Fig.2 Flow chart of ACL data processing

理論相關函數的準確性對速度估計結果的精度起著決定性作用。在確定了理論相關函數后，模型擬合的精確程度決定了速度估計結果的精度。

理論相關函數也被稱作混響時空相關模型，最早由國外的Dickey 推導給出[2]，國內最早由Zhu等[11]根據Faure，Ol’shevskii 和Middleton 提出的FOM 模型推導建立，隨后經過馮雷博士、蒯多杰博士以及易卉琴等研究人員的不斷改進后[3-5]，目前采用的理論相關函數為

其中：τ是相關時延，V是三維速度，Φ(τ)是發射信號相關函數，k=ω0/c是波數。

在確定了底回波的理論相關函數后，速度估計方法就是決定性能和精度的關鍵。目前采用的速度解算方法是以局域最小二乘法為核心的參數分布估計法。該方法采用了水平速度和垂向速度分步估計的思路，利用局域最小二乘法估計水平速度，利用解析方法求解垂向速度，是目前最有效的速度解算方法。

針對以局域最小二乘法為核心的參數估計方法，當局域范圍確定后，以最小二乘確定目標函數，搜索使目標函數最小的參數變量。最優參數求解流程如圖3所示。

圖3 最優參數求解流程圖Fig.3 Flow chart of solving optimal parameters

搜索方向和步長采用擬牛頓法，用近似黑塞(Hessian)矩陣來避免Hessian 矩陣求取時的復雜運算。近似Hessian 矩陣采用的是1970 年前后由Broyden、Fletcher、Goldfarb 以及Shanno 提出BFGS算法，具有不用求逆和不用計算二階導數等優點，在此用Bk代替，Bk的迭代公式為

式中：sk=θk-θk-1，yk=gk-gk-1。sk代表目標參數θk的差值，yk代表一階梯度gk的差值。根據簡化后的近似Hessian矩陣即根據式(3)，不斷迭代來獲得一定相對偏差Δgk下的目標參數，進而解算出速度信息。

此算法在參數的求解過程中需要多次迭代，參數多，過程復雜，而且需要考慮眾多的邊界條件限制。為了簡化求解流程，降低算法復雜度，本文經過系統的理論推導，提出了基于矩陣推導的模型參數求解方法。

2 基于矩陣推導的模型參數求解方法

基本思路：分離變量使得模型函數可以通過矩陣表示，然后基于最小二乘進行理論推導，最終根據參與運算的數據點直接推導出可以使誤差最小的參數，整個過程不需要迭代運算。

基于矩陣的推導方法的目的同迭代方法一樣，也是獲得最優的模型參數。因此分析還是要基于底回波混響相關模型。目前所采用的底回波模型是經過添加橢圓系數以及旋轉角后的模型，模型如式(5)所示。首要目標是把底回波模型中的參數分離，轉化為矩陣形式：

對式(5)進行歸一化處理后，可得：

為了運算清晰，令x0=2Vxτ，y0=2Vyτ，替換后如式(7)所示：

下面基于式(7)從理論上分析求解最優參數的可行性。首先對式(7)取對數，可得

為了將運算轉變為矩陣運算，對式(8)進行變量分離，結果為

式(9)已經完成了變量的分離。下一步將基于式(9)推導最優參數。假如參與擬合的矢量點有N個，則將這N個矢量數據點寫成矩陣的形式為A=BC，其中：A為N×1的列向量，其元素為

B為N×6型的矩陣，其行向量為

C為6×1型的列向量：

式中：

N點數據的誤差列向量為：E=A-BC，采用最小二乘擬合，使N個數據點的均方誤差(Mean Square Erro,MSE)最小，即：

為了簡化計算，將矩陣B進行QR分解，即B=QR，分解后Q是一個N×N型的正交矩陣，R是一個N×6型的上三角矩陣，對E=A-BC兩邊同時乘以QT，推導可得：

誤差項為

令：

其中，S為6 維列向量，T為N×6 維列向量，R1為6×6型的上三角矩陣，則式(14)可以進一步轉化為

從式(18)很容易可以看出，當S=R1C時，式(14)可以取得最小值，因此只需要求得最小值時的C即可，由S=R1C可得：

根據式(13)以及式(19)得到的結果解線性方程組，即可求得參數x0和y0：

已知x0=2Vxτ、y0=2Vyτ，再根據式(20)即可解算出速度信息。

3 數據驗證

3.1 仿真數據驗證

通過一個二維旋轉高斯模型產生樣本數據，將樣本數據輸入新的模型參數求解方法，根據結果進行驗證。

采用的二維高斯旋轉模型如式(21)所示。分析結果如表1所示。

表1 仿真數據驗證結果Table 1 Verification results of simulation data

驗證方案為固定模型參數，通過改變樣本點數，針對不同樣本點數得到的結果驗證求解方案的可行性。仿真參數設定x0=0.5，y0=1，θ=30°，a=0.8，從中心點向外隨機取點，分別在選取樣本點為10×10、8×8、5×5、3×3 的情況下進行了測試，運算結果如表1所示。

從表1可以看出，在不同采樣點下，仿真測試結果都能正確地解算出目標參數，在沒有噪聲的情況下偏差基本為0。測試結果驗證了基于矩陣推導求解最優參數的可行性。

3.2 實驗數據驗證

采用湖上試驗數據進行結果驗證，試驗設備為75 kHz聲相關計程儀(ACL)，比測設備為GPS(型號為ProPak-G2)，湖況兩級左右，測試水域水深65 m左右，各航速情況下的試驗數據都在4 200幀左右。矩陣化處理程序和原程序處理結果如表2、3所示。

表2 不同航速下湖試數據的矩陣運算結果Table 2 Matrix calculation results of lake-trial data at different speeds

由于矩陣化運算中舍棄了幅值調節因子，因此所使用的測速范圍會有一定程度的下降。速度在7 kn(1 kn=1.852 km·h-1)以上測速偏差明顯增大，表2中只列出7 kn及以下航速下的運算結果。

對比表2、表3 中不同航速下湖試數據的矩陣運算和原程序運算的結果可以看出，雖然矩陣運算結果在航速5 kn和7 kn時相對偏差比原程序大，但是偏差方向一致，而且整體測速結果穩定，測速偏差范圍在0.2%以內。原程序運算結果在不同航速下偏差有正有負。在測速結果整體一致性上矩陣運算優于原程序求解。

表3 不同航速下湖試數據的原程序處理結果Table 3 Original program processing results of lake-trial data at different speeds

另外，以上運算所使用的樣本點數(接收陣元兩兩之間構成相關矢量，對應的相關值被稱為樣本點)是基于大量數據驗證，是原程序求解方法最優的選擇。考慮到矩陣求解法舍棄了幅值調節因子，造成一部分樣本點的模型失配，對此采取了減少樣本點和增加樣本點的嘗試。減少樣本點是盡量利用滿足模型的樣本點進行運算。增加樣本點是在樣本點不足的情況下，樣本點新增和帶來的模型失配之間會有一個最佳樣本點數，使得測速取得最優結果。通過大量的數據驗證，本文獲得了以下經驗結果：

(1) 基于式(1)可以看出，測速范圍與相關延時和相關矢量相關，陣元布局決定了相關矢量，當相關延時加倍時，測速范圍就會減半，也就是說相同的相關矢量數量，所能表示的測速范圍降低為原來的一半，這樣會使得同樣速度范圍內可以利用的樣本點數增多，因此通過將相關延時更改為2τ來使相關運算的衰減速度降低，增加可以利用的相關點數；

(2) 原程序是在原模型下進行運算，在12個樣本點時結果最優。而本文基于矩陣運算的方法，對使用的模型進行了一定調整，因此需要通過調整參與運算的樣本點數來搜索更適合矩陣運算的樣本點數。通過對比不同點數下的測速結果，發現2τ延時下采用15點做速度解算得到的測速結果最優(速度小于7 kn)，而且測速均方差也降低到原來的70%左右。另一方面，對航速7～10 kn情況也針對不同點數做了進一步分析。發現當參與相關運算的點數為13 時測速結果穩定。不同速度下的矩陣運算測速結果如表4所示。整體偏差方向一致，不同速度下的相對偏差在0.3%以內。以上結果進一步驗證了矩陣化求解算法的可行性。

4 結 論

本文通過理論分析得到了基于矩陣推導的最優參數的解析解，并通過仿真驗證和湖上試驗數據，驗證了該方案作為一種速度解算方法的可行性。另外通過對時間復雜度分析，原程序由于需要求逆運算和多次迭代，時間復雜度為M×N3量級，其中M代表迭代運算次數。基于矩陣推導的運算方法雖然也涉及矩陣求逆，但是為上三角矩陣求逆，整體時間復雜度為N2量級，其中乘法運算和加法運算次數相當。因此，基于矩陣推導的最優參數求解方法相較于基于擬牛頓法的迭代參數求解，在時間復雜度上具有一定的優勢。

本文由于無法應用幅值調節因子，因此帶來了一定的模型失配問題，導致文中所提測速方案沒有達到最優的測速結果。后續研究將考慮從提升可以利用的樣本點數和提高矩陣推導算法的容錯性兩方面進一步深入研究，以獲得更加精確的估計參數。