多波束系统安装偏差整体校准的地形特征匹配方法

李 铁，周丰年，赵建虎

1. 武汉大学测绘学院，湖北 武汉 430079； 2. 武汉大学海洋研究院，湖北 武汉 430079； 3. 长江水利委员会水文局长江口水文水资源勘测局，上海 200136

多波束回声测深系统(multibeam echo sounders，MBES)因其宽覆盖、高效率、高精度等优点，已成为目前水下地形测量的主要仪器[1-2]。换能器是MBES的核心声学单元，是测深结果的起算参考。受安装不精准、航行中支架松动等影响，换能器活性面不水平或其轴向与船体坐标系轴向不平行，导致测深断面与航迹方向不正交，严重影响了声线跟踪和测深点在地理坐标系下坐标的计算[3-5]。

多波束系统安装偏差校准方法通常利用特殊地形的海底数据以特定的顺序进行求取。Patch test方法主要基于去耦[6-7]原理校准偏差，其偏差校准顺序多样。文献[8—9]推荐的偏差校准顺序为时延—纵摇偏差—横摇偏差—艏摇偏差校准；文献[10]推荐的顺序为横摇偏差—时延—纵摇偏差—艏摇偏差校准；国际海道测量组织(International Hydrographic Organization，IHO)建议使用时延—纵摇偏差—艏摇偏差—横摇偏差的校准顺序[11]；交通运输行业标准JTT790—2010[12]中提出时延—横摇偏差—纵摇偏差—艏摇偏差的校准顺序；文献[13]提出了多波束系统安装校准的方法和横摇偏差—艏摇偏差—时延—纵摇偏差的校准顺序。不同的处理顺序主要是为了消除多波束不同安装偏差间的干扰或耦合效应。文献[14]系统地阐述了校准试验法(patch test)，该方法已成为多波束系统安装偏差校准通用方法。文献[15—16]研究了多波束系统安装偏差对海底地形测量的影响，提出了利用平面拟合校正多波束横向安装偏差的思想，但存在数据量大、运算速度慢和只能校准横向偏差等问题；文献[17]提出了一种多波束横摇角度偏差二次校准方法，有效地削弱了粗差和微小地形起伏对海底倾角的影响，但忽视了多波束系统安装偏差影响的综合性。文献[6—7，18]针对单项或顺序校准去耦合性原理不完善、校准工作量大、偏差校准结果不唯一问题开展了多波束系统安装偏差的整体校准研究。文献[18]提出了一种安装偏差整体校准校准方法。该方法通过波束点归位计算公式推导出高差与坡度、多波束系统安装偏差、平移量之间的函数关系，由于其假设的方位角很小，所以布设测线时要严格为南北向和东西向，在实际作业时存在很大的约束性。文献[6—7]提出使用MIBAC(multibeam-IMU boresight automatic calibration)进行安装偏差整体校准。此方法对光滑地形进行平面拟合，建立波束点坐标与安装偏差间的函数关系，以波束点到平面的距离和最小为约束条件，整体解算多波束系统安装偏差。由于多波束系统安装偏差与平面方程参数一次性解算，模型复杂，实施困难。为此，本文提出了一种安装偏差整体校准的新方法。该方法通过研究各偏差对测深影响，形成综合误差模型进而整体解算，实现多波束系统安装偏差的整体获取。下面详细介绍这种方法。

1 整体校准方法原理

测深点的地理坐标(XL,YL,ZL)可借助以下公式进行计算[13]

(1)

(2)

水平方向的分速度可由船速v和航向κ计算，垂直方向的分速度设为0

(3)

考虑ω是一个小量，可近似为ω≈0，对式(1)进行泰勒级数展开并保留一阶项可得

(4)

式中，(X0,Y0,Z0)表示测深点的概略地理坐标；(dXL,dYL,dZL)表示坐标偏差

(5)

式中，y、h分别表示声线跟踪横向距离和深度；dS/S是一个衡量斜距测量精度的尺度参数，与声速精度有关；dt、dκ、dω、dφ为时延、航向偏差、纵摇偏差和横摇偏差。

若条带1、2公共位置上测点的测量坐标分别为F1、F2，真实坐标为F0，两点点位偏差为dF1、dF2，则根据式(4)可得

(6)

则有

dF1-dF2-(F1-F2)=0

(7)

dF1、dF2按式(5)计算，F1、F2按照式(4)计算。若认为两条带测量中多波束的安装偏差相等，将式(7)展开，其矩阵形式为

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；由于时延dt不能为负，所以对dt进行约束，约束条件为GX≥W，即

(9)

利用平差准则[19-21]min{(BX-L)T(BX-L)|dt≥0}，对附有不等式约束的平差进行求解

(10)

得到求取的待估参数，利用式(5)计算各坐标的改正值，然后通过式(4)重新计算坐标的概略值。由于声线跟踪的水平距离y和深度h受波束俯角φ和尺度参数dS/S的影响，所以需要进行修正，计算公式如式(11)所示

(11)

通过不断迭代计算，对应坐标差值小于给定限差即终止迭代。迭代终止条件设置为dφ、dω和dκ均小于0.01×π/180。由于每次得到的dφ、dω和dκ是相对上次的偏差量，因此对各次的dφ、dω和dκ分别叠加，最终得到多波束的安装偏差。

2 地形匹配点对的提取

若两测深条带存在公共部分，则公共位置存在两次测量点对，即匹配点对[22]。借助这些匹配点对可利用以上模型，获得多波束系统安装偏差。匹配点对可借助人工获得，但人为选择地形影响匹配精度和最终的安装偏差精度。水下地形最直观的表现形式为三维模型，由于测深点数量大，点云匹配耗时长。为提高匹配效率和精度，将三维数据利用二维灰度图像表述，借助图像匹配技术实现匹配点对的快速获取。地形灰度图像是将离散的测深数据网格化，将高程转化为0～255灰度级，形成地形图像。当测区地形的高程(深度)变化较小时，这种转换会进一步凸显地形变化，当地形变化较大时，255个色阶不能保证地形的精度，可以转化为16级灰度地形图。为确保地形的精细度，网格化时取多波束纵向和横向测深数据间隔平均值作为网格大小，如式(12)所示

(12)

(13)

获得了地形图像后，利用图像中反映地物起伏的特征点对形成由公共覆盖的多波束条带匹配点对。匹配点对借助SIFT(scale-invariant feature transform)来寻找。SIFT算法所查找的关键点是不因光照、仿射变换、噪声等因素而变化的突出点，具有尺度不变、多量性、独特性好、信息量丰富和消除边缘响应等优点[23-24]。SIFT将图像利用不同标准差进行高斯模糊构建尺度空间。通过降采样构建高斯金字塔，金字塔每一组为不同大小图像，同一组包含多幅不同尺度空间的图像。同一组相邻两层不同尺度空间的图像生成高斯差分图像，通过邻域比较获得局部极值点，并对极值点进行曲线拟合获得关键点，然后进行关键点描述。SIFT图像匹配步骤如下：

(1) 极值检测。搜索所有尺度上的图像位置，识别对于尺度和旋转不变的兴趣点。

(2) 关键点定位。在每个候选位置上，通过一个拟合精细的模型来确定位置和尺度。

(3) 方向确定。基于图像局部梯度方向，分配给每个关键点一个或多个方向。所有对图像操作均是相对关键点方向、尺度和位置的变换。

(4) 关键点的描述。在每个关键点周围邻域内，在选定的尺度上测量图像局部的梯度。

SIFT匹配获得的点对中会存在误匹配点对，进而会影响后续的安装偏差计算精度。RANSAC算法可确保匹配点对的可靠性[25-26]。本文借助RANSAC算法对SIFT提取出的匹配点对进行检核，剔除异常点对。通过SIFT匹配，可以获得两幅图像的匹配点像素坐标，对像素坐标变换，得到各匹配点的地理坐标(Xm,Ym,Zm)

很多疾病都是因胸痛就诊而被检出发现[6]，心血管系统造成的急性胸痛是临床常见病因[7-8]，如急性肺动脉栓塞、急性主动脉夹层等疾病。急性胸痛具有较高的发病率和死亡率，临床采用X线检查、心电图诊断急性胸痛[9-10]，前者对肺实变性的显示存在不足之处，而心电图虽然空白期较为短暂，但仍会导致误诊或漏诊情况发生，因此应选择一种更加准确有效的方法明确急性胸痛的病因，便于尽早实施对症治疗。

(14)

(15)

由此，匹配点对在各自图像中的三维坐标求得，进而求得匹配点对之间的坐标差、平均值和单位权中误差，利用2σ原则对坐标差进行误差剔除。基于这些点对和对应的坐标差，构建式(8)方程组，借助式(8)—(11)整体解算多波束的安装偏差。

3 试验及分析

为了验证本文提出算法的正确性，借助Sonic 2024多波束在水深11～18 m的水域开展了4条测线测量。试验区地形较平坦；测量时，Sonic 2024的ping更新频率设置为25 Hz、扇开角为140°、波束个数为256。测量时，未进行多波束系统安装偏差校正，完成了4条长度均约为450 m的测线测量，其中条带1、3为同向重复测线测量结果、条带2、4为反向同测线测量，条带1(或3)与条带2(或4)间公共覆盖度约为40%。测量期间条带3的船速约为3.5 m/s，其余约为2 m/s。4条测线的分布及形成的地形图如图1所示。

通过解算4条测线数据，可得测点在换能器坐标系和地理坐标系下的平面坐标和水深。借助这些数据根据如下过程实现多波束偏差的整体解算：

(1) 地形图像生成。由船速和ping更新频率可得沿航迹方向的相邻ping间隔为0.14 m，垂直于航迹方向距中央波束2/3开角处的间隔为0.29 m，格网大小最终设置为二者的均值0.22 m。对测深数据按照0.22 m的格网插值，并将测点高程变换为灰度，形成4个条带的地形灰度图像如图2和图3所示。

(2) 借助SIFT算法对条带1—2和2—4进行图像匹配，并借助RANSAC算法对异常的匹配点对剔除，并将正确的匹配点对连线如图2和图3所示。1—2条带总匹配点对数为181，有效匹配点对数为136。2—4条带总匹配点数为951，有效匹配点对数为844。

(3) 利用这些匹配点对，借助式(14)和式(15)计算匹配点对地理坐标。求取匹配点对坐标差，统计各坐标分量的平均值、标准差，并绘制偏差直方图，如图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可以看出，匹配点对的数量足够大时，其坐标偏差服从正态分布。

(4) 利用计算得到的匹配点对、声线跟踪横坐标与水深、速度和航向数据，建立式(8)误差模型，通过对时延约束和加权迭代，对未知参数进行求解，得到多波束系统安装偏差及时延。

图1 试验区域水下地形 图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 匹配点对坐标差折线图Fig.4 Coordinate differences curves of matching point pairs

图5 匹配点对坐标差直方图Fig.5 Histograms of matching point pairs coordinate differences

多波束边缘波束点精度较低，边缘波束匹配点对的对应坐标偏差较大，对多波束系统安装偏差解算结果会产生影响。为此裁减掉条带左右边缘波束各40个波束点，利用中央波束点通过上述方法再次进行多波束系统安装偏差解算，并与Patch test结果比较，结果如表2所示。可以看出，相对Patch test结果，本文方法与Caris结果相近。分析认为，裁剪掉边缘波束，虽提高了参与计算的点对精度，但因参与计算的数据量减少，解算的冗余度也会随之降低。

表1 本文方法计算的多波束系统安装偏差与Caris计算结果对比

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

表2 裁剪后对应条带平差结果与Caris计算结果对比

由图2可以发现，原始条带1、2的覆盖度为40%，但匹配点对数量不多；由图3可知，原始条带2、4为往返测量线，匹配点对多，且多分布于靠近中央测线附近。所以，裁减边缘波束对于往返或同向测线匹配点对数量的减少影响较小，对计算的精度也影响较小。因此，在多波束系统安装偏差校准时，可利用同一测线同向或相向测量数据开展多波束系统安装偏差校准。

由于上述试验所求取的多波束系统安装偏差较小，但实际安装偏差一般位于[-5°，5°]，所以将多波束系统安装偏差人为增大，roll角偏差增大2.8°，pitch角偏差增加-3°，yaw角不变，再利用本文方法对多波束系统安装偏差进行校准。其中1—2条带的有效匹配点对数量为35；2—4条带的有效匹配点对数量为182。得到的结果如表3所示。

表3 增大换能器安装偏差后对应条带平差结果与Caris计算结果对比

由表3可以看出，对于换能器安装偏差角较大的情况，此方法也可以适用，但是本文提出的方法与Patch test方法差值变大。出现此种情况的主要原因是换能器安装偏差大会导致匹配点对的位置变化较大，通过RANSAC算法和2σ原则对匹配点对筛选时删除的匹配点对多，导致正确的匹配点对数量下降。因此，需要通过增加测线的长度来保证足够的匹配点对数量。

由表1—3可知，校正的时延量很小，主要由于MBES测量中采用了GNSS的1PPS同步技术。若测量中均采用该技术，则可对式(1)右侧关于速度与时间乘积的第二项删除，从而实现模型的简化。

对于两种多波束系统安装偏差校准方式，从多波束系统安装偏差校准原理上说，本文方法充分考虑了多波束系统安装偏差角度之间的耦合性，相对于Patch test方法更为严谨。多波束系统安装偏差校正后相应格网点的高差如图6所示。

由图6中的数据标签可以看出，两种方法在边缘波束相差不超过4 cm，中央波束的差值极小，而本文方法相对于Patch test方法相应点位的高差更小。图6(a)横向上高程差表现为左负右正的趋势，说明Patch test方法求取的多波束系统安装偏差有一定的残差。图6(b)整体高程差较均匀，没有明显的趋势，说明本文方法对多波束系统安装偏差的校正更为彻底。

4 结论与建议

本文提出的多波束系统安装偏差整体校准方法，考虑了多波束系统安装偏差对测深结果影响的耦合性，理论上更加完备，获得的偏差精度更高。本文方法适用于较平坦但有微地貌的地形(沙波、乱石区等除外)，无需刻意寻找特定地形开展测量，可利用多波束作业中的同测线往返或同向测量结果即可实现多波束系统安装偏差的计算，因此简化了校准作业流程。

为进一步提高校准精度，建议测量时布设至少2条长度大于500 m的测线开展往返或同向测量。此外，计算时建议对波束入射角大于60°的边缘波束实施裁剪，利用中央波束测深点对开展多波束系统安装偏差计算。

图6 2—4测线格网点高差Fig.6 Elevation difference of 2 & 4 survey route line grid nodes