相關(guān)測速聲吶模型參數(shù)矩陣化求解方法

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

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

0 引 言

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

1 基本工作原理

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

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

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

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

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

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

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

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

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

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

圖2 ACL數(shù)據(jù)處理流程圖Fig.2 Flow chart of ACL data processing

理論相關(guān)函數(shù)的準(zhǔn)確性對速度估計結(jié)果的精度起著決定性作用。在確定了理論相關(guān)函數(shù)后，模型擬合的精確程度決定了速度估計結(jié)果的精度。

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

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

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

針對以局域最小二乘法為核心的參數(shù)估計方法，當(dāng)局域范圍確定后，以最小二乘確定目標(biāo)函數(shù)，搜索使目標(biāo)函數(shù)最小的參數(shù)變量。最優(yōu)參數(shù)求解流程如圖3所示。

圖3 最優(yōu)參數(shù)求解流程圖Fig.3 Flow chart of solving optimal parameters

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

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

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

2 基于矩陣推導(dǎo)的模型參數(shù)求解方法

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

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

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

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

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

為了將運(yùn)算轉(zhuǎn)變?yōu)榫仃囘\(yùn)算，對式(8)進(jìn)行變量分離，結(jié)果為

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

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

C為6×1型的列向量：

式中：

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

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

誤差項為

令：

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

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

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

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

3 數(shù)據(jù)驗證

3.1 仿真數(shù)據(jù)驗證

通過一個二維旋轉(zhuǎn)高斯模型產(chǎn)生樣本數(shù)據(jù)，將樣本數(shù)據(jù)輸入新的模型參數(shù)求解方法，根據(jù)結(jié)果進(jìn)行驗證。

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

表1 仿真數(shù)據(jù)驗證結(jié)果Table 1 Verification results of simulation data

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

從表1可以看出，在不同采樣點(diǎn)下，仿真測試結(jié)果都能正確地解算出目標(biāo)參數(shù)，在沒有噪聲的情況下偏差基本為0。測試結(jié)果驗證了基于矩陣推導(dǎo)求解最優(yōu)參數(shù)的可行性。

3.2 實(shí)驗數(shù)據(jù)驗證

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

表2 不同航速下湖試數(shù)據(jù)的矩陣運(yùn)算結(jié)果Table 2 Matrix calculation results of lake-trial data at different speeds

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

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

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

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

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

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

4 結(jié) 論

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

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