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

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

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

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

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

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

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

(1) 坐標(biāo)參考系。為構(gòu)建傳播曲面模型，需定義相關(guān)坐標(biāo)參考系。本文傳感器陣列坐標(biāo)系：以陣列中心為坐標(biāo)原點(diǎn)，x軸指向船左舷方向，y軸指向航向，z軸與x、y軸構(gòu)成右手正交坐標(biāo)系；當(dāng)?shù)厮阶鴺?biāo)系：以傳感器中心為坐標(biāo)原點(diǎn)，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)換至船中心為原點(diǎn)的水平坐標(biāo)系或地理坐標(biāo)系下。

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

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

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

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

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

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

式(3)

(1)

(2)

|TVi|=1

(3)

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

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

(4)

φi=90-arcsin(ZTVi)

(5)

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

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

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

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

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

具體實(shí)施如下：

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

步驟2 插值結(jié)點(diǎn)位置計(jì)算。

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

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

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

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

(5) 若MR1-R2

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

步驟3 待插值點(diǎn)位置計(jì)算。

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

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

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

2.1.1 截點(diǎn)間距上限MR、MT確定

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

若不考慮折射，收/發(fā)傳播曲面是由收/發(fā)矢量錐面延伸形成的圓錐面，任意等深面內(nèi)的收/發(fā)截線均為圓錐曲線[22]。如圖5所示，A點(diǎn)為實(shí)際收/發(fā)截線的交點(diǎn)，B點(diǎn)為收/發(fā)截點(diǎn)連線的交點(diǎn)；εR、εT表示收/發(fā)截線以直代曲的偏差，當(dāng)εR、εT的最大值均小于0.5ε0時(shí)，以B點(diǎn)代替A點(diǎn)的偏差小于ε0。εR、εT最大值為：半徑為截線最小曲率半徑，弦長(zhǎng)等于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)

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

(7)

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

式(6)—式(7)的推導(dǎo)基于傳播曲面為圓錐面，因此MT和MR并非精度ε0下的嚴(yán)密推導(dǎo)結(jié)果。但式(6)中的不等號(hào)可在一定程度上保證估算的有效性。大量試驗(yàn)表明，一般情況下，當(dāng)截點(diǎn)間距小于MT和MR時(shí)，以上計(jì)算相對(duì)嚴(yán)密迭代計(jì)算偏差小于ε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，對(duì)αT1、αT2和αR1、αR2對(duì)應(yīng)的4個(gè)矢量聲線跟蹤，保留其在D1和D2等深面內(nèi)截點(diǎn)，其可能的分布如圖6所示。圖6(a)和圖6(b)中發(fā)射截線和接收截線不相交，需調(diào)整初始角度αR1、αR2、αT1、αT2。對(duì)于圖6(a)情況，按比例向αT2方向平移(αT1，αT2)；對(duì)于圖6(b)，向αR1方向平移(αR1，αR2)；最終收斂情況下截點(diǎn)分布情況如圖6(c)所示。

圖6 截點(diǎn)分布矯正角度范圍Fig.6 Correction angle range of intercept

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

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

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

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

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

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

(9)

(10)

2.2.1.2 不同地形條件下多項(xiàng)式擬合殘差

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

(11)

(12)

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

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

3 試驗(yàn)分析及討論

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

圖9 測(cè)線地形圖Fig.9 Survey line topographic map

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

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

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

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

由表1可知，交叉點(diǎn)深度絕對(duì)差異均值、絕對(duì)差異標(biāo)準(zhǔn)差、相對(duì)差異均值3種統(tǒng)計(jì)結(jié)果的表現(xiàn)為：在淺水區(qū)，本文算法與Caris相近，差異小于12%，且本文算法絕對(duì)差異均值更小，但標(biāo)準(zhǔn)差和平均相對(duì)差異略大；而傳統(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計(jì)算精度接近，傳統(tǒng)算法的計(jì)算精度略低；在中深水區(qū)，本文算法和Caris計(jì)算精度相近，傳統(tǒng)算法精度顯著偏低。這可能是由于隨著測(cè)量深度增加，收發(fā)傳感器間距離增大，傳統(tǒng)模型下的波束入射角和發(fā)射中心誤差增大，導(dǎo)致計(jì)算精度偏低。統(tǒng)計(jì)結(jié)果中，中深水區(qū)差異并非0均值分布，這主要是由于此地區(qū)的聲速剖面數(shù)據(jù)質(zhì)量不高，導(dǎo)致邊緣波束向兩側(cè)翹起，以及較大的地形起伏導(dǎo)致的精度不均一。

3.2 同號(hào)波束歸算差異對(duì)比

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

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

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

表2 同號(hào)波束腳印位置較差

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

3.3 位置歸算效率分析

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

表3 單個(gè)波束平均效率對(duì)比

4 結(jié) 論

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