李 鐵，周豐年，趙建虎

1. 武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079； 2. 武漢大學(xué)海洋研究院，湖北 武漢 430079； 3. 長(zhǎng)江水利委員會(huì)水文局長(zhǎng)江口水文水資源勘測(cè)局，上海 200136

多波束回聲測(cè)深系統(tǒng)(multibeam echo sounders，MBES)因其寬覆蓋、高效率、高精度等優(yōu)點(diǎn)，已成為目前水下地形測(cè)量的主要儀器[1-2]。換能器是MBES的核心聲學(xué)單元，是測(cè)深結(jié)果的起算參考。受安裝不精準(zhǔn)、航行中支架松動(dòng)等影響，換能器活性面不水平或其軸向與船體坐標(biāo)系軸向不平行，導(dǎo)致測(cè)深斷面與航跡方向不正交，嚴(yán)重影響了聲線(xiàn)跟蹤和測(cè)深點(diǎn)在地理坐標(biāo)系下坐標(biāo)的計(jì)算[3-5]。

多波束系統(tǒng)安裝偏差校準(zhǔn)方法通常利用特殊地形的海底數(shù)據(jù)以特定的順序進(jìn)行求取。Patch test方法主要基于去耦[6-7]原理校準(zhǔn)偏差，其偏差校準(zhǔn)順序多樣。文獻(xiàn)[8—9]推薦的偏差校準(zhǔn)順序?yàn)闀r(shí)延—縱搖偏差—橫搖偏差—艏搖偏差校準(zhǔn)；文獻(xiàn)[10]推薦的順序?yàn)闄M搖偏差—時(shí)延—縱搖偏差—艏搖偏差校準(zhǔn)；國(guó)際海道測(cè)量組織(International Hydrographic Organization，IHO)建議使用時(shí)延—縱搖偏差—艏搖偏差—橫搖偏差的校準(zhǔn)順序[11]；交通運(yùn)輸行業(yè)標(biāo)準(zhǔn)JTT790—2010[12]中提出時(shí)延—橫搖偏差—縱搖偏差—艏搖偏差的校準(zhǔn)順序；文獻(xiàn)[13]提出了多波束系統(tǒng)安裝校準(zhǔn)的方法和橫搖偏差—艏搖偏差—時(shí)延—縱搖偏差的校準(zhǔn)順序。不同的處理順序主要是為了消除多波束不同安裝偏差間的干擾或耦合效應(yīng)。文獻(xiàn)[14]系統(tǒng)地闡述了校準(zhǔn)試驗(yàn)法(patch test)，該方法已成為多波束系統(tǒng)安裝偏差校準(zhǔn)通用方法。文獻(xiàn)[15—16]研究了多波束系統(tǒng)安裝偏差對(duì)海底地形測(cè)量的影響，提出了利用平面擬合校正多波束橫向安裝偏差的思想，但存在數(shù)據(jù)量大、運(yùn)算速度慢和只能校準(zhǔn)橫向偏差等問(wèn)題；文獻(xiàn)[17]提出了一種多波束橫搖角度偏差二次校準(zhǔn)方法，有效地削弱了粗差和微小地形起伏對(duì)海底傾角的影響，但忽視了多波束系統(tǒng)安裝偏差影響的綜合性。文獻(xiàn)[6—7，18]針對(duì)單項(xiàng)或順序校準(zhǔn)去耦合性原理不完善、校準(zhǔn)工作量大、偏差校準(zhǔn)結(jié)果不唯一問(wèn)題開(kāi)展了多波束系統(tǒng)安裝偏差的整體校準(zhǔn)研究。文獻(xiàn)[18]提出了一種安裝偏差整體校準(zhǔn)校準(zhǔn)方法。該方法通過(guò)波束點(diǎn)歸位計(jì)算公式推導(dǎo)出高差與坡度、多波束系統(tǒng)安裝偏差、平移量之間的函數(shù)關(guān)系，由于其假設(shè)的方位角很小，所以布設(shè)測(cè)線(xiàn)時(shí)要嚴(yán)格為南北向和東西向，在實(shí)際作業(yè)時(shí)存在很大的約束性。文獻(xiàn)[6—7]提出使用MIBAC(multibeam-IMU boresight automatic calibration)進(jìn)行安裝偏差整體校準(zhǔn)。此方法對(duì)光滑地形進(jìn)行平面擬合，建立波束點(diǎn)坐標(biāo)與安裝偏差間的函數(shù)關(guān)系，以波束點(diǎn)到平面的距離和最小為約束條件，整體解算多波束系統(tǒng)安裝偏差。由于多波束系統(tǒng)安裝偏差與平面方程參數(shù)一次性解算，模型復(fù)雜，實(shí)施困難。為此，本文提出了一種安裝偏差整體校準(zhǔn)的新方法。該方法通過(guò)研究各偏差對(duì)測(cè)深影響，形成綜合誤差模型進(jìn)而整體解算，實(shí)現(xiàn)多波束系統(tǒng)安裝偏差的整體獲取。下面詳細(xì)介紹這種方法。

1 整體校準(zhǔn)方法原理

測(cè)深點(diǎn)的地理坐標(biāo)(XL,YL,ZL)可借助以下公式進(jìn)行計(jì)算[13]

(1)

(2)

水平方向的分速度可由船速v和航向κ計(jì)算，垂直方向的分速度設(shè)為0

(3)

考慮ω是一個(gè)小量，可近似為ω≈0，對(duì)式(1)進(jìn)行泰勒級(jí)數(shù)展開(kāi)并保留一階項(xiàng)可得

(4)

式中，(X0,Y0,Z0)表示測(cè)深點(diǎn)的概略地理坐標(biāo)；(dXL,dYL,dZL)表示坐標(biāo)偏差

(5)

式中，y、h分別表示聲線(xiàn)跟蹤橫向距離和深度；dS/S是一個(gè)衡量斜距測(cè)量精度的尺度參數(shù)，與聲速精度有關(guān)；dt、dκ、dω、dφ為時(shí)延、航向偏差、縱搖偏差和橫搖偏差。

若條帶1、2公共位置上測(cè)點(diǎn)的測(cè)量坐標(biāo)分別為F1、F2，真實(shí)坐標(biāo)為F0，兩點(diǎn)點(diǎn)位偏差為dF1、dF2，則根據(jù)式(4)可得

(6)

則有

dF1-dF2-(F1-F2)=0

(7)

dF1、dF2按式(5)計(jì)算，F(xiàn)1、F2按照式(4)計(jì)算。若認(rèn)為兩條帶測(cè)量中多波束的安裝偏差相等，將式(7)展開(kāi)，其矩陣形式為

V=BX-L

(8)

式中，a11=h2sinκ2-h1sinκ1；a21=h1cosκ1-h2cosκ2；a31=y2-y1；a12=h1cosκ1-h2cosκ2；a22=h1sinκ1-h2sinκ2；a32=0；a13=y1cosκ1-y2cosκ2；a23=y2sinκ2-y1sinκ1；a33=0；a14=vx1-vx2；a24=vy1-vy2；a34=0；a15=y2sinκ2-y1sinκ1；a25=y1cosκ1-y2cosκ2；a35=h1-h2；由于時(shí)延dt不能為負(fù)，所以對(duì)dt進(jìn)行約束，約束條件為GX≥W，即

(9)

利用平差準(zhǔn)則[19-21]min{(BX-L)T(BX-L)|dt≥0}，對(duì)附有不等式約束的平差進(jìn)行求解

(10)

得到求取的待估參數(shù)，利用式(5)計(jì)算各坐標(biāo)的改正值，然后通過(guò)式(4)重新計(jì)算坐標(biāo)的概略值。由于聲線(xiàn)跟蹤的水平距離y和深度h受波束俯角φ和尺度參數(shù)dS/S的影響，所以需要進(jìn)行修正，計(jì)算公式如式(11)所示

(11)

通過(guò)不斷迭代計(jì)算，對(duì)應(yīng)坐標(biāo)差值小于給定限差即終止迭代。迭代終止條件設(shè)置為dφ、dω和dκ均小于0.01×π/180。由于每次得到的dφ、dω和dκ是相對(duì)上次的偏差量，因此對(duì)各次的dφ、dω和dκ分別疊加，最終得到多波束的安裝偏差。

2 地形匹配點(diǎn)對(duì)的提取

若兩測(cè)深條帶存在公共部分，則公共位置存在兩次測(cè)量點(diǎn)對(duì)，即匹配點(diǎn)對(duì)[22]。借助這些匹配點(diǎn)對(duì)可利用以上模型，獲得多波束系統(tǒng)安裝偏差。匹配點(diǎn)對(duì)可借助人工獲得，但人為選擇地形影響匹配精度和最終的安裝偏差精度。水下地形最直觀的表現(xiàn)形式為三維模型，由于測(cè)深點(diǎn)數(shù)量大，點(diǎn)云匹配耗時(shí)長(zhǎng)。為提高匹配效率和精度，將三維數(shù)據(jù)利用二維灰度圖像表述，借助圖像匹配技術(shù)實(shí)現(xiàn)匹配點(diǎn)對(duì)的快速獲取。地形灰度圖像是將離散的測(cè)深數(shù)據(jù)網(wǎng)格化，將高程轉(zhuǎn)化為0～255灰度級(jí)，形成地形圖像。當(dāng)測(cè)區(qū)地形的高程(深度)變化較小時(shí)，這種轉(zhuǎn)換會(huì)進(jìn)一步凸顯地形變化，當(dāng)?shù)匦巫兓^大時(shí)，255個(gè)色階不能保證地形的精度，可以轉(zhuǎn)化為16級(jí)灰度地形圖。為確保地形的精細(xì)度，網(wǎng)格化時(shí)取多波束縱向和橫向測(cè)深數(shù)據(jù)間隔平均值作為網(wǎng)格大小，如式(12)所示

(12)

(13)

獲得了地形圖像后，利用圖像中反映地物起伏的特征點(diǎn)對(duì)形成由公共覆蓋的多波束條帶匹配點(diǎn)對(duì)。匹配點(diǎn)對(duì)借助SIFT(scale-invariant feature transform)來(lái)尋找。SIFT算法所查找的關(guān)鍵點(diǎn)是不因光照、仿射變換、噪聲等因素而變化的突出點(diǎn)，具有尺度不變、多量性、獨(dú)特性好、信息量豐富和消除邊緣響應(yīng)等優(yōu)點(diǎn)[23-24]。SIFT將圖像利用不同標(biāo)準(zhǔn)差進(jìn)行高斯模糊構(gòu)建尺度空間。通過(guò)降采樣構(gòu)建高斯金字塔，金字塔每一組為不同大小圖像，同一組包含多幅不同尺度空間的圖像。同一組相鄰兩層不同尺度空間的圖像生成高斯差分圖像，通過(guò)鄰域比較獲得局部極值點(diǎn)，并對(duì)極值點(diǎn)進(jìn)行曲線(xiàn)擬合獲得關(guān)鍵點(diǎn)，然后進(jìn)行關(guān)鍵點(diǎn)描述。SIFT圖像匹配步驟如下：

(1) 極值檢測(cè)。搜索所有尺度上的圖像位置，識(shí)別對(duì)于尺度和旋轉(zhuǎn)不變的興趣點(diǎn)。

(2) 關(guān)鍵點(diǎn)定位。在每個(gè)候選位置上，通過(guò)一個(gè)擬合精細(xì)的模型來(lái)確定位置和尺度。

(3) 方向確定?；趫D像局部梯度方向，分配給每個(gè)關(guān)鍵點(diǎn)一個(gè)或多個(gè)方向。所有對(duì)圖像操作均是相對(duì)關(guān)鍵點(diǎn)方向、尺度和位置的變換。

(4) 關(guān)鍵點(diǎn)的描述。在每個(gè)關(guān)鍵點(diǎn)周?chē)徲騼?nèi)，在選定的尺度上測(cè)量圖像局部的梯度。

SIFT匹配獲得的點(diǎn)對(duì)中會(huì)存在誤匹配點(diǎn)對(duì)，進(jìn)而會(huì)影響后續(xù)的安裝偏差計(jì)算精度。RANSAC算法可確保匹配點(diǎn)對(duì)的可靠性[25-26]。本文借助RANSAC算法對(duì)SIFT提取出的匹配點(diǎn)對(duì)進(jìn)行檢核，剔除異常點(diǎn)對(duì)。通過(guò)SIFT匹配，可以獲得兩幅圖像的匹配點(diǎn)像素坐標(biāo)，對(duì)像素坐標(biāo)變換，得到各匹配點(diǎn)的地理坐標(biāo)(Xm,Ym,Zm)

很多疾病都是因胸痛就診而被檢出發(fā)現(xiàn)[6]，心血管系統(tǒng)造成的急性胸痛是臨床常見(jiàn)病因[7-8]，如急性肺動(dòng)脈栓塞、急性主動(dòng)脈夾層等疾病。急性胸痛具有較高的發(fā)病率和死亡率，臨床采用X線(xiàn)檢查、心電圖診斷急性胸痛[9-10]，前者對(duì)肺實(shí)變性的顯示存在不足之處，而心電圖雖然空白期較為短暫，但仍會(huì)導(dǎo)致誤診或漏診情況發(fā)生，因此應(yīng)選擇一種更加準(zhǔn)確有效的方法明確急性胸痛的病因，便于盡早實(shí)施對(duì)癥治療。

(14)

(15)

由此，匹配點(diǎn)對(duì)在各自圖像中的三維坐標(biāo)求得，進(jìn)而求得匹配點(diǎn)對(duì)之間的坐標(biāo)差、平均值和單位權(quán)中誤差，利用2σ原則對(duì)坐標(biāo)差進(jìn)行誤差剔除。基于這些點(diǎn)對(duì)和對(duì)應(yīng)的坐標(biāo)差，構(gòu)建式(8)方程組，借助式(8)—(11)整體解算多波束的安裝偏差。

3 試驗(yàn)及分析

為了驗(yàn)證本文提出算法的正確性，借助Sonic 2024多波束在水深11～18 m的水域開(kāi)展了4條測(cè)線(xiàn)測(cè)量。試驗(yàn)區(qū)地形較平坦；測(cè)量時(shí)，Sonic 2024的ping更新頻率設(shè)置為25 Hz、扇開(kāi)角為140°、波束個(gè)數(shù)為256。測(cè)量時(shí)，未進(jìn)行多波束系統(tǒng)安裝偏差校正，完成了4條長(zhǎng)度均約為450 m的測(cè)線(xiàn)測(cè)量，其中條帶1、3為同向重復(fù)測(cè)線(xiàn)測(cè)量結(jié)果、條帶2、4為反向同測(cè)線(xiàn)測(cè)量，條帶1(或3)與條帶2(或4)間公共覆蓋度約為40%。測(cè)量期間條帶3的船速約為3.5 m/s，其余約為2 m/s。4條測(cè)線(xiàn)的分布及形成的地形圖如圖1所示。

通過(guò)解算4條測(cè)線(xiàn)數(shù)據(jù)，可得測(cè)點(diǎn)在換能器坐標(biāo)系和地理坐標(biāo)系下的平面坐標(biāo)和水深。借助這些數(shù)據(jù)根據(jù)如下過(guò)程實(shí)現(xiàn)多波束偏差的整體解算：

(1) 地形圖像生成。由船速和ping更新頻率可得沿航跡方向的相鄰ping間隔為0.14 m，垂直于航跡方向距中央波束2/3開(kāi)角處的間隔為0.29 m，格網(wǎng)大小最終設(shè)置為二者的均值0.22 m。對(duì)測(cè)深數(shù)據(jù)按照0.22 m的格網(wǎng)插值，并將測(cè)點(diǎn)高程變換為灰度，形成4個(gè)條帶的地形灰度圖像如圖2和圖3所示。

(2) 借助SIFT算法對(duì)條帶1—2和2—4進(jìn)行圖像匹配，并借助RANSAC算法對(duì)異常的匹配點(diǎn)對(duì)剔除，并將正確的匹配點(diǎn)對(duì)連線(xiàn)如圖2和圖3所示。1—2條帶總匹配點(diǎn)對(duì)數(shù)為181，有效匹配點(diǎn)對(duì)數(shù)為136。2—4條帶總匹配點(diǎn)數(shù)為951，有效匹配點(diǎn)對(duì)數(shù)為844。

(3) 利用這些匹配點(diǎn)對(duì)，借助式(14)和式(15)計(jì)算匹配點(diǎn)對(duì)地理坐標(biāo)。求取匹配點(diǎn)對(duì)坐標(biāo)差，統(tǒng)計(jì)各坐標(biāo)分量的平均值、標(biāo)準(zhǔn)差，并繪制偏差直方圖，如圖4和圖5所示。其中，1—2條帶X、Y、Z方向的均值分別為0.190、-0.306和0.007 m，中誤差分別為0.303、0.380和0.07 m；2—4條帶X、Y、Z方向的均值分別為-0.089、-0.302和-0.007 m，中誤差分別為0.230、0.384和0.031 m。由圖5可以看出，匹配點(diǎn)對(duì)的數(shù)量足夠大時(shí)，其坐標(biāo)偏差服從正態(tài)分布。

(4) 利用計(jì)算得到的匹配點(diǎn)對(duì)、聲線(xiàn)跟蹤橫坐標(biāo)與水深、速度和航向數(shù)據(jù)，建立式(8)誤差模型，通過(guò)對(duì)時(shí)延約束和加權(quán)迭代，對(duì)未知參數(shù)進(jìn)行求解，得到多波束系統(tǒng)安裝偏差及時(shí)延。

圖1 試驗(yàn)區(qū)域水下地形 圖2 1—2條帶地形SIFT匹配 圖3 2—4條帶地形SIFT匹配Fig.1 Terrain map of test area Fig.2 Terrain matching of line 1 & 2 Fig.3 Terrain matching of line 2 & 4

圖4 匹配點(diǎn)對(duì)坐標(biāo)差折線(xiàn)圖Fig.4 Coordinate differences curves of matching point pairs

圖5 匹配點(diǎn)對(duì)坐標(biāo)差直方圖Fig.5 Histograms of matching point pairs coordinate differences

多波束邊緣波束點(diǎn)精度較低，邊緣波束匹配點(diǎn)對(duì)的對(duì)應(yīng)坐標(biāo)偏差較大，對(duì)多波束系統(tǒng)安裝偏差解算結(jié)果會(huì)產(chǎn)生影響。為此裁減掉條帶左右邊緣波束各40個(gè)波束點(diǎn)，利用中央波束點(diǎn)通過(guò)上述方法再次進(jìn)行多波束系統(tǒng)安裝偏差解算，并與Patch test結(jié)果比較，結(jié)果如表2所示?？梢钥闯?，相對(duì)Patch test結(jié)果，本文方法與Caris結(jié)果相近。分析認(rèn)為，裁剪掉邊緣波束，雖提高了參與計(jì)算的點(diǎn)對(duì)精度，但因參與計(jì)算的數(shù)據(jù)量減少，解算的冗余度也會(huì)隨之降低。

表1 本文方法計(jì)算的多波束系統(tǒng)安裝偏差與Caris計(jì)算結(jié)果對(duì)比

方法rolld?/(°)pitchdω/(°)yawdκ/(°)latencydt/sratedS/S本文方法 0.162-0.1550.5100-0.0012Patch test0.200-0.1800.4800—差值-0.0380.0250.0310—

表2 裁剪后對(duì)應(yīng)條帶平差結(jié)果與Caris計(jì)算結(jié)果對(duì)比

由圖2可以發(fā)現(xiàn)，原始條帶1、2的覆蓋度為40%，但匹配點(diǎn)對(duì)數(shù)量不多；由圖3可知，原始條帶2、4為往返測(cè)量線(xiàn)，匹配點(diǎn)對(duì)多，且多分布于靠近中央測(cè)線(xiàn)附近。所以，裁減邊緣波束對(duì)于往返或同向測(cè)線(xiàn)匹配點(diǎn)對(duì)數(shù)量的減少影響較小，對(duì)計(jì)算的精度也影響較小。因此，在多波束系統(tǒng)安裝偏差校準(zhǔn)時(shí)，可利用同一測(cè)線(xiàn)同向或相向測(cè)量數(shù)據(jù)開(kāi)展多波束系統(tǒng)安裝偏差校準(zhǔn)。

由于上述試驗(yàn)所求取的多波束系統(tǒng)安裝偏差較小，但實(shí)際安裝偏差一般位于[-5°，5°]，所以將多波束系統(tǒng)安裝偏差人為增大，roll角偏差增大2.8°，pitch角偏差增加-3°，yaw角不變，再利用本文方法對(duì)多波束系統(tǒng)安裝偏差進(jìn)行校準(zhǔn)。其中1—2條帶的有效匹配點(diǎn)對(duì)數(shù)量為35；2—4條帶的有效匹配點(diǎn)對(duì)數(shù)量為182。得到的結(jié)果如表3所示。

表3 增大換能器安裝偏差后對(duì)應(yīng)條帶平差結(jié)果與Caris計(jì)算結(jié)果對(duì)比

由表3可以看出，對(duì)于換能器安裝偏差角較大的情況，此方法也可以適用，但是本文提出的方法與Patch test方法差值變大。出現(xiàn)此種情況的主要原因是換能器安裝偏差大會(huì)導(dǎo)致匹配點(diǎn)對(duì)的位置變化較大，通過(guò)RANSAC算法和2σ原則對(duì)匹配點(diǎn)對(duì)篩選時(shí)刪除的匹配點(diǎn)對(duì)多，導(dǎo)致正確的匹配點(diǎn)對(duì)數(shù)量下降。因此，需要通過(guò)增加測(cè)線(xiàn)的長(zhǎng)度來(lái)保證足夠的匹配點(diǎn)對(duì)數(shù)量。

由表1—3可知，校正的時(shí)延量很小，主要由于MBES測(cè)量中采用了GNSS的1PPS同步技術(shù)。若測(cè)量中均采用該技術(shù)，則可對(duì)式(1)右側(cè)關(guān)于速度與時(shí)間乘積的第二項(xiàng)刪除，從而實(shí)現(xiàn)模型的簡(jiǎn)化。

對(duì)于兩種多波束系統(tǒng)安裝偏差校準(zhǔn)方式，從多波束系統(tǒng)安裝偏差校準(zhǔn)原理上說(shuō)，本文方法充分考慮了多波束系統(tǒng)安裝偏差角度之間的耦合性，相對(duì)于Patch test方法更為嚴(yán)謹(jǐn)。多波束系統(tǒng)安裝偏差校正后相應(yīng)格網(wǎng)點(diǎn)的高差如圖6所示。

由圖6中的數(shù)據(jù)標(biāo)簽可以看出，兩種方法在邊緣波束相差不超過(guò)4 cm，中央波束的差值極小，而本文方法相對(duì)于Patch test方法相應(yīng)點(diǎn)位的高差更小。圖6(a)橫向上高程差表現(xiàn)為左負(fù)右正的趨勢(shì)，說(shuō)明Patch test方法求取的多波束系統(tǒng)安裝偏差有一定的殘差。圖6(b)整體高程差較均勻，沒(méi)有明顯的趨勢(shì)，說(shuō)明本文方法對(duì)多波束系統(tǒng)安裝偏差的校正更為徹底。

4 結(jié)論與建議

本文提出的多波束系統(tǒng)安裝偏差整體校準(zhǔn)方法，考慮了多波束系統(tǒng)安裝偏差對(duì)測(cè)深結(jié)果影響的耦合性，理論上更加完備，獲得的偏差精度更高。本文方法適用于較平坦但有微地貌的地形(沙波、亂石區(qū)等除外)，無(wú)需刻意尋找特定地形開(kāi)展測(cè)量，可利用多波束作業(yè)中的同測(cè)線(xiàn)往返或同向測(cè)量結(jié)果即可實(shí)現(xiàn)多波束系統(tǒng)安裝偏差的計(jì)算，因此簡(jiǎn)化了校準(zhǔn)作業(yè)流程。

為進(jìn)一步提高校準(zhǔn)精度，建議測(cè)量時(shí)布設(shè)至少2條長(zhǎng)度大于500 m的測(cè)線(xiàn)開(kāi)展往返或同向測(cè)量。此外，計(jì)算時(shí)建議對(duì)波束入射角大于60°的邊緣波束實(shí)施裁剪，利用中央波束測(cè)深點(diǎn)對(duì)開(kāi)展多波束系統(tǒng)安裝偏差計(jì)算。

圖6 2—4測(cè)線(xiàn)格網(wǎng)點(diǎn)高差Fig.6 Elevation difference of 2 & 4 survey route line grid nodes