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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式(3)

(1)

(2)

|TVi|=1

(3)

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

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

(4)

φi=90-arcsin(ZTVi)

(5)

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

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

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

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

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

具體實施如下：

步驟1 在扇區內，等間隔選取5個波束作為插值結點，其他波束作為待插值點。

步驟2 插值結點位置計算。

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

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

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

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

(5) 若MR1-R2

(6) 根據tD1、tD2、t0賦權，由D1、D2等深面內四組截點坐標插值得到波束腳印位置。

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

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

同理完成全部待插值點歸算，即實現全扇區波束腳印位置歸算。

2.1 插值結點波束腳印位置歸算

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

不考慮截線的幾何特性，當迭代終止條件設為MR1-R2<ε0、MT1-T2<ε0時，可得到嚴密結果。但為了減少迭代次數，上述算法結合截線的最小曲率半徑，計算滿足ε0時的截點間距上限MR、MT，將迭代終止條件設為MR1-R2

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

圖5 D0等深面內收/發截線Fig.5 Receiving and sending intercept lines in the D0 isobaric plane

(6)

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

(7)

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

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

2.1.2 迭代區間初值和矯正

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

(8)

設置αR1、αR2、αT1、αT2、D1、D2的初值為-rαR、rαR、-rαT、rαT、D0-rd、D0+rd，對αT1、αT2和αR1、αR2對應的4個矢量聲線跟蹤，保留其在D1和D2等深面內截點，其可能的分布如圖6所示。圖6(a)和圖6(b)中發射截線和接收截線不相交，需調整初始角度α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 待插值波束腳印位置歸算

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

2.2.1 聲線跟蹤參數插值函數建立

2.2.1.1 φT、tR/tT插值函數建立

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

圖7 扇區內參數關系Fig.7 Parameter relationship in the sector

(9)

(10)

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

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

(11)

(12)

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

2.2.2 參數插值及位置歸算

3 試驗分析及討論

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

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

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

3.1 交叉線公共覆蓋區格網點差異對比

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

表1 交叉線公共覆蓋區交叉點深度差異

由表1可知，交叉點深度絕對差異均值、絕對差異標準差、相對差異均值3種統計結果的表現為：在淺水區，本文算法與Caris相近，差異小于12%，且本文算法絕對差異均值更小，但標準差和平均相對差異略大；而傳統算法差異均值大于本文算法60.0%，差異標準差大于本文算法10.0%。在中深水區，本文算法與Caris結果相近，差異小于2.5%；傳統算法差異均值大于本文算法60.1%，差異標準差高于本文算法124.1%。

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

3.2 同號波束歸算差異對比

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

由表2可知，在淺水區，偏差均值為5～7 cm，偏差標準差為3～5 cm，相對偏差小于0.35‰Z；在中深水區，偏差均值為7～13 cm，偏差標準差為5～14 cm，相對偏差小于0.11‰Z，符合相關規范[25]。由圖10可知，淺水區，90%以上的偏差集中在0～10 cm；中深水區，90%以上的偏差集中在0～30 cm。這表明本文方法計算結果與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 位置歸算效率分析

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

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

4 結 論

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