顧及傳播曲面的多波束波束腳印高精度快速歸位算法

畢自軍，趙建虎，鄭 根，劉美琴

1. 武漢大學(xué)測繪學(xué)院，湖北 武漢 430079； 2. 武漢大學(xué)海洋研究院，湖北 武漢 430079

海底地形是海洋基礎(chǔ)地理信息的重要組成部分，多波束測深系統(tǒng)是高效獲取海底地形的典型設(shè)備之一[1-3]，位置歸算是獲取高精度多波束測深點的重要環(huán)節(jié)[4-5]，為此，國內(nèi)外學(xué)者進(jìn)行了大量研究。文獻(xiàn)[6]提出了多波束精密聲線跟蹤算法和測深點歸位計算模型，并指出姿態(tài)角對歸位計算精度影響顯著；文獻(xiàn)[7]分析了姿態(tài)誤差對歸位計算的影響，并利用傅里葉變換對其進(jìn)行去除；文獻(xiàn)[7—8]推導(dǎo)了顧及姿態(tài)的聲線跟蹤模型，并在淺水區(qū)驗證了其有效性；文獻(xiàn)[8—9]通過坐標(biāo)旋轉(zhuǎn)建立了顧及船姿的波束入射向量計算模型，提高了聲線跟蹤和歸位計算精度。以上研究都假設(shè)波束發(fā)射和接收過程共路線或路線對稱。文獻(xiàn)[10—12]指出，受測量船運動和姿態(tài)瞬時變化的影響，上述假設(shè)與實際存在差異，這種影響在淺水區(qū)并不明顯，但隨著船速增大和水深增加，歸位計算誤差越來越大，嚴(yán)重影響多波束測深精度。

文獻(xiàn)[11]提出虛擬同心陣(virtual concentric array，VCCA)模型，根據(jù)收/發(fā)傳感器姿態(tài)和波束指向角建立收/發(fā)矢量錐面，且假設(shè)波束的傳播路徑為：從收/發(fā)時刻傳感器位置的中點處發(fā)出并在海底散射后沿原路徑返回。VCCA顧及收/發(fā)位置差異，一定程度上改善中深水的波束腳印歸位計算精度，但其所假設(shè)的波束傳播路徑降低了歸算精度。文獻(xiàn)[12]在VCCA模型基礎(chǔ)上建立了非同心陣(non concentric array，NCCA)模型，即在平行于收/發(fā)陣的平面族內(nèi)以非同心雙曲線交點估計波束腳印，根據(jù)聲線跟蹤的雙程時間與觀測時間，迭代調(diào)整目標(biāo)平面與換能器間距離，完成波束腳印位置歸算。NCCA顧及了波束收/發(fā)陣列位置的不同，改善VCCA的理論精度，但仍存在多次迭代導(dǎo)致的效率較低等問題。文獻(xiàn)[13]對VCCA模型進(jìn)行了推導(dǎo)和簡化，提高了計算效率，但計算精度與VCCA模型近似。

為此，本文基于多波束測量原理和波束傳播理論，建立波束傳播曲面模型，進(jìn)而提出一種高精度高效率的波束腳印位置歸算方法。

1 位置歸算基本原理及傳播曲面模型

文獻(xiàn)[12]指出波束腳印為收/發(fā)波束在海底形成各自波束腳印的交集，歸算步驟如下。

(1) 坐標(biāo)參考系。為構(gòu)建傳播曲面模型，需定義相關(guān)坐標(biāo)參考系。本文傳感器陣列坐標(biāo)系：以陣列中心為坐標(biāo)原點，x軸指向船左舷方向，y軸指向航向，z軸與x、y軸構(gòu)成右手正交坐標(biāo)系；當(dāng)?shù)厮阶鴺?biāo)系：以傳感器中心為坐標(biāo)原點，X軸指向地北子午線方向，Y軸指向東，Z軸與X軸、Y軸構(gòu)成右手正交坐標(biāo)系[14-15]。以下步驟(2)—步驟(3)基于當(dāng)?shù)厮阶鴺?biāo)系，各矢量均需經(jīng)傳感器坐標(biāo)系旋轉(zhuǎn)至當(dāng)?shù)厮阶鴺?biāo)系[11]；步驟(4)中收/發(fā)傳播曲面需轉(zhuǎn)換至船中心為原點的水平坐標(biāo)系或地理坐標(biāo)系下。

(3) 根據(jù)聲速剖面，建立收/發(fā)波束傳播曲面。傳播聲線為各RV/TV對應(yīng)的聲波傳播路徑，可根據(jù)聲速、聲線方位角θ、俯角φ，借助聲線跟蹤獲得(圖1(b))。收/發(fā)矢量錐面上所有矢量對應(yīng)的傳播聲線集合構(gòu)成收/發(fā)傳播曲面(圖1(c))。下稱傳播曲面與等深面F的交線為截線，聲線與等深面F的交點為截點。

(4) 確定實際波束腳印。各個等深面內(nèi)可確定發(fā)射截線與接收截線交點(圖1(d))，所有等深面內(nèi)該交點的連線即為收/發(fā)傳播曲面的交線，由雙程傳播時間可在該交線上確定波束腳印位置。

圖1 波束腳印位置歸算原理Fig.1 The principle of beam footprint position reduction

以上步驟中，傳播曲面的建立是確定截線并最終獲得波束腳印的關(guān)鍵。以發(fā)射為例，該過程在當(dāng)?shù)厮阶鴺?biāo)系下可分為4步：①確定發(fā)射錐面主軸方向單位矢量(TX)和錐面半頂角；②計算遍歷初始矢量(Tfirst)；③遍歷全部TV；④計算各TV對應(yīng)θ、φ并根據(jù)聲線跟蹤形成聲線，聲線集合即為發(fā)射傳播曲面。

圖2 遍歷波束發(fā)射/接收矢量Fig.2 Traverse transmit/receive vector

則對于錐面上任意矢量TVi滿足式(1)—

式(3)

(1)

(2)

|TVi|=1

(3)

根據(jù)式(1)—式(3)，對每一個αi均可計算其對應(yīng)矢量TVi為(XTVi，YTVi，ZTVi)。因此，通過在[-180°,180°]內(nèi)遍歷α即可計算錐面上每一個TV的矢量坐標(biāo)，實現(xiàn)遍歷。

矢量TVi對應(yīng)的方位角θi和聲線俯角φi由式(4)—式(5)計算

(4)

φi=90-arcsin(ZTVi)

(5)

基于θi和φi，根據(jù)常梯度聲線跟蹤[16-20]，保留所有收/發(fā)聲線路徑，即可形成收/發(fā)傳播曲面。在收/發(fā)傳播曲面交線上找到收/發(fā)傳播時間等于觀測雙程時間的點，即為波束腳印，圖3描述了該過程。

圖3 傳播曲面模型Fig.3 The propagation surface model

2 顧及傳播曲面模型的高效歸算算法

傳播曲面模型理論嚴(yán)密，但建立收/發(fā)傳播曲面耗時較長，為此下面給出一種高效、高精度歸算算法。

為提高計算效率，本文采用如下策略：將扇區(qū)內(nèi)波束分為插值結(jié)點和待插值點；前者用迭代搜索，替代完整傳播曲面的建立；后者用插值結(jié)點歸算結(jié)果開展參數(shù)插值，避免多次迭代(圖4)。

具體實施如下：

步驟1 在扇區(qū)內(nèi)，等間隔選取5個波束作為插值結(jié)點，其他波束作為待插值點。

步驟2 插值結(jié)點位置計算。

(1) 利用VCCA模型計算起始矢量Tfirst、概略深度D0和概略聲程L0。

(2) 由D0和平面精度ε0，計算等深面內(nèi)，兩個發(fā)射或接收截點間的直線距離(下稱發(fā)射或接收截點間距)的上限MT、MR。

(3) 結(jié)合ε0、D0、L0設(shè)置迭代收/發(fā)遍歷角度α的初始區(qū)間(αR1，αR2)、(αT1，αT2)和深度區(qū)間(D1，D2)。

(4) 分別對αR1、αR2對應(yīng)的兩個接收矢量和αT1、αT2對應(yīng)的兩個發(fā)射矢量聲線跟蹤，統(tǒng)計這4根聲線在D1、D2兩個等深面內(nèi)的收發(fā)傳播時間和tD1、tD2及D2等深面內(nèi)收/發(fā)截點間距MR1-R2、MT1-T2。

(5) 若MR1-R2

(6) 根據(jù)tD1、tD2、t0賦權(quán)，由D1、D2等深面內(nèi)四組截點坐標(biāo)插值得到波束腳印位置。

步驟3 待插值點位置計算。

(4) 由φTi、θTi聲線跟蹤至tTi耗盡，即為第i個待求波束腳印位置。

同理完成全部待插值點歸算，即實現(xiàn)全扇區(qū)波束腳印位置歸算。

2.1 插值結(jié)點波束腳印位置歸算

2.1.1 截點間距上限MR、MT確定

不考慮截線的幾何特性，當(dāng)?shù)K止條件設(shè)為MR1-R2<ε0、MT1-T2<ε0時，可得到嚴(yán)密結(jié)果。但為了減少迭代次數(shù)，上述算法結(jié)合截線的最小曲率半徑，計算滿足ε0時的截點間距上限MR、MT，將迭代終止條件設(shè)為MR1-R2

若不考慮折射，收/發(fā)傳播曲面是由收/發(fā)矢量錐面延伸形成的圓錐面，任意等深面內(nèi)的收/發(fā)截線均為圓錐曲線[22]。如圖5所示，A點為實際收/發(fā)截線的交點，B點為收/發(fā)截點連線的交點；εR、εT表示收/發(fā)截線以直代曲的偏差，當(dāng)εR、εT的最大值均小于0.5ε0時，以B點代替A點的偏差小于ε0。εR、εT最大值為：半徑為截線最小曲率半徑，弦長等于MR、MT的圓弧的拱高。推導(dǎo)可得MR/MT滿足式(6)，式中ρR和ρT表示收/發(fā)截線最小曲率半徑

圖5 D0等深面內(nèi)收/發(fā)截線Fig.5 Receiving and sending intercept lines in the D0 isobaric plane

(6)

圓錐曲線最小曲率半徑為圓錐頂點到截平面的距離與圓錐半頂角正切值的乘積[22]，因此，ρR、ρT為收/發(fā)傳感器到波束腳印深度與收/發(fā)錐面半頂角的乘積

(7)

式中，D表示傳感器吃水深度。將式(7)代入式(6)即可得到MR、MT。

式(6)—式(7)的推導(dǎo)基于傳播曲面為圓錐面，因此MT和MR并非精度ε0下的嚴(yán)密推導(dǎo)結(jié)果。但式(6)中的不等號可在一定程度上保證估算的有效性。大量試驗表明，一般情況下，當(dāng)截點間距小于MT和MR時，以上計算相對嚴(yán)密迭代計算偏差小于ε0。

2.1.2 迭代區(qū)間初值和矯正

合理設(shè)置迭代角度初始區(qū)間(αR1，αR2)、(αT1，αT2)和深度初始區(qū)間(D1，D2)，可減少迭代次數(shù)。因此，本文根據(jù)MT和MR，估算角度區(qū)間半徑rαR、rαT和深度區(qū)間半徑rd

(8)

設(shè)置αR1、αR2、αT1、αT2、D1、D2的初值為-rαR、rαR、-rαT、rαT、D0-rd、D0+rd，對αT1、αT2和αR1、αR2對應(yīng)的4個矢量聲線跟蹤，保留其在D1和D2等深面內(nèi)截點，其可能的分布如圖6所示。圖6(a)和圖6(b)中發(fā)射截線和接收截線不相交，需調(diào)整初始角度αR1、αR2、αT1、αT2。對于圖6(a)情況，按比例向αT2方向平移(αT1，αT2)；對于圖6(b)，向αR1方向平移(αR1，αR2)；最終收斂情況下截點分布情況如圖6(c)所示。

圖6 截點分布矯正角度范圍Fig.6 Correction angle range of intercept

2.2 待插值波束腳印位置歸算

由于聲線跟蹤至少需要發(fā)射矢量方位角θT和俯角φT、發(fā)射聲線單程傳播時間tT[12,23]。本文利用插值計算以上參數(shù)，避免多次迭代和聲線跟蹤，提高算法效率。

2.2.1 聲線跟蹤參數(shù)插值函數(shù)建立

2.2.1.1 φT、tR/tT插值函數(shù)建立

為了分析影響發(fā)射矢量方位角φT和單程傳播時間tT的主要參數(shù)，進(jìn)行如下假設(shè)：發(fā)射扇面垂直向下；接收傳感器位置不變；接收傳感器無姿態(tài)；水體無折射。

圖7 扇區(qū)內(nèi)參數(shù)關(guān)系Fig.7 Parameter relationship in the sector

(9)

(10)

2.2.1.2 不同地形條件下多項式擬合殘差

本文等間隔選取n+1個點，建立扇區(qū)內(nèi)tR/tT、φT的擬合函數(shù)如式(11)、式(12)所示

(11)

(12)

圖8 平坦、傾斜、曲線地形下多項式擬合殘差Fig.8 Polynomial fitting residuals under flat,sloping,and curved terrain

2.2.2 參數(shù)插值及位置歸算

3 試驗分析及討論

本文在南海某區(qū)中深水和淺水各選擇兩條交叉線，測量儀器為Kongsberg EM302，地形如圖9所示。其中圖9(a)為淺水區(qū)，測線1水深180～200 m，東南-西北方向，面積約4.5×106m2，共2209 ping；測線2(圖9為部分)水深160～240 m，西南-東北方向，共16 064 ping；兩條測線每ping均432個波束，4扇區(qū)，交叉區(qū)域約1×106m2。圖9(b)為中深水區(qū)，測線3水深900～1400 m，西南-東北方向，面積約1.8×108m2，共2630 ping，測線4(圖9為部分)水深900～1600 m，西北-東南方向，共2 362 ping；兩條測線每ping均432個波束，8扇區(qū)，交叉區(qū)域約為5×107m2。

圖9 測線地形圖Fig.9 Survey line topographic map

數(shù)據(jù)預(yù)處理步驟包括：聲速剖面等數(shù)據(jù)質(zhì)量控制和船文件編輯等。歸算使用相同的原始數(shù)據(jù)、聲速剖面、船文件，并分為3種方法：①Caris使用HIPS and SIPS11.1中Georeference Bathymetry模塊的有聲速無潮位模式；②本文算法；③顧及姿態(tài)及聲線彎曲的歸算模型[8-9](下稱傳統(tǒng)算法)。后兩種方法均經(jīng)過波束腳印相對傳感器偏移量計算，坐標(biāo)轉(zhuǎn)換等步驟。下文對比3種算法在交叉測線的公共覆蓋區(qū)內(nèi)交叉點精度差異；本文算法與Caris同號波束點的互差；本文算法與Caris平均單個波束的運算時間，以驗證本文算法的有效性和運行效率。

3.1 交叉線公共覆蓋區(qū)格網(wǎng)點差異對比

在測線1和測線2、測線3和測線4的公共覆蓋區(qū)中各選擇13 200個交叉點(一條測線的中央波束與同組另一條測線的交點)?；?種方法的歸算結(jié)果，在淺水區(qū)和中深水區(qū)分別計算各測線在交叉點處深度，統(tǒng)計測線1-2，測線3-4交叉點深度差異(表1)。

表1 交叉線公共覆蓋區(qū)交叉點深度差異

由表1可知，交叉點深度絕對差異均值、絕對差異標(biāo)準(zhǔn)差、相對差異均值3種統(tǒng)計結(jié)果的表現(xiàn)為：在淺水區(qū)，本文算法與Caris相近，差異小于12%，且本文算法絕對差異均值更小，但標(biāo)準(zhǔn)差和平均相對差異略大；而傳統(tǒng)算法差異均值大于本文算法60.0%，差異標(biāo)準(zhǔn)差大于本文算法10.0%。在中深水區(qū)，本文算法與Caris結(jié)果相近，差異小于2.5%；傳統(tǒng)算法差異均值大于本文算法60.1%，差異標(biāo)準(zhǔn)差高于本文算法124.1%。

這表明，在淺水區(qū)，本文算法、Caris計算精度接近，傳統(tǒng)算法的計算精度略低；在中深水區(qū)，本文算法和Caris計算精度相近，傳統(tǒng)算法精度顯著偏低。這可能是由于隨著測量深度增加，收發(fā)傳感器間距離增大，傳統(tǒng)模型下的波束入射角和發(fā)射中心誤差增大，導(dǎo)致計算精度偏低。統(tǒng)計結(jié)果中，中深水區(qū)差異并非0均值分布，這主要是由于此地區(qū)的聲速剖面數(shù)據(jù)質(zhì)量不高，導(dǎo)致邊緣波束向兩側(cè)翹起，以及較大的地形起伏導(dǎo)致的精度不均一。

3.2 同號波束歸算差異對比

淺水區(qū)使用測線1，中深水區(qū)使用測線3，統(tǒng)計本文算法與Caris的同號波束計算結(jié)果互差(表2)，并形成偏差頻率分布條形圖(圖10)。

由表2可知，在淺水區(qū)，偏差均值為5～7 cm，偏差標(biāo)準(zhǔn)差為3～5 cm，相對偏差小于0.35‰Z；在中深水區(qū)，偏差均值為7～13 cm，偏差標(biāo)準(zhǔn)差為5～14 cm，相對偏差小于0.11‰Z，符合相關(guān)規(guī)范[25]。由圖10可知，淺水區(qū)，90%以上的偏差集中在0～10 cm；中深水區(qū)，90%以上的偏差集中在0～30 cm。這表明本文方法計算結(jié)果與Caris具有較好的一致性，驗證了本文算法的有效性。

圖10 本文算法與Caris同號波束腳印偏差Fig.10 The beam footprint calculate by algorithm in this paper and Caris deviation

表2 同號波束腳印位置較差

圖11 歸算差異與接收角度Fig.11 Reduction difference and receiving angle

3.3 位置歸算效率分析

為驗證本文算法有較高的運算效率，分別統(tǒng)計淺水區(qū)測線1和中深水區(qū)測線3的全部波束，計算本文算法與Caris單個波束平均運算時間，以及本文算法的平均聲線跟蹤次數(shù)(表3)。由表3可知，聲線跟蹤作為算法中耗時占比較高的部分，本文算法的平均次數(shù)略高于傳統(tǒng)算法的1次，低于NCCA模型的6～8次。在淺水區(qū)，本文算法效率相對Caris提高8.22%；中深水區(qū)，相對Caris提高35.21%。兩區(qū)域效率均有一定程度提高，且中深水區(qū)相對淺水區(qū)效率提高更明顯。

表3 單個波束平均效率對比

4 結(jié) 論

相比于傳統(tǒng)模型、VCCA模型、NCCA模型，本文算法提出的傳播曲面模型準(zhǔn)確還原了波束腳印位置歸算過程， 顧及了收發(fā)傳感器的位置差異及收發(fā)傳播時間不相等的問題， 理論基礎(chǔ)更為嚴(yán)謹(jǐn)。在此基礎(chǔ)上，通過迭代搜索和參數(shù)插值的方式顯著減少聲線跟蹤次數(shù)，提高計算效率，實現(xiàn)基于傳播曲面模型的高效位置歸算。經(jīng)實測數(shù)據(jù)驗證，本文算法交叉線公共覆蓋區(qū)交叉點誤差與Caris相近；與Caris同號波束腳印歸算結(jié)果較差，淺水區(qū)差異小于0.3‰Z，中深水區(qū)差異小于0.1‰Z；且兩測區(qū)內(nèi)運算效率相較Caris均有所提高。隨著深度的增加，本文結(jié)果與Caris的相對偏差值，有一定程度的降低，且運行效率相對更高，深水區(qū)有更好的適用性。在測量原理相同，且已記錄發(fā)射、接收指向角、雙程時間等參數(shù)的多波束測深數(shù)據(jù)中，本文算法適用。