-
受侧扫声呐扫描成图机理和拖曳作业模式[1]制约，其瀑布图像存在倾斜几何畸变、辐射畸变、压缩畸变、坐标不准、时变局部畸变[2]等不足，并给其精处理及镶嵌工作带来了困难。长期以来，侧扫声呐图像仅用于海床目标的定性分析。然而，侧扫声呐图像作为建设“数字海洋”“透明海洋”的重要数据源之一，其处理质量和成果形式不仅影响海床目标的正确判读与识别，更决定了其能否与其他地理空间数据准确融合和再利用。如何形成高质量、大区域侧扫声呐镶嵌图像已成为海洋数字化工程中一项亟待解决的难题。
在拖鱼位置、姿态信息难以准确获得的情况下，为实现相邻条带图像公共区域的目标级拼接，国内外学者基于尺度不变特征转换(scale-invariant feature transform, SIFT)、快速稳健特征(speeded-up robust features, SURF) [3]、目标阴影分布网形[4-6]等特征和统计参数开展了侧扫声呐条带图像的配准工作。然而，由于侧扫声呐图像的噪声、灰差变化以及大小和类型变化的局部畸变，上述方法特征点对提取较少、相似度测量困难、相邻条带配准精度低，难以实现相邻条带间共视目标的保形性镶嵌。
本文研究了侧扫声呐条带数据处理理论[7]，精处理后的条带图像有效地保证了SIFT/SURF特征的提取数量和质量[8]，但SURF算法用于侧扫声呐条带图像配准时，仍然存在特征检测遍历空间大、效率低、局部变化的畸变引起的图像各处变换参数不一致、整体配准精度低等问题。因此，本文基于SURF，将薄板样条函数引入分块弹性配准机制，以实现侧扫声呐相邻条带之间的特征级镶嵌。
-
侧扫声呐瀑布图像的预处理工作主要包括海底追踪、灰度增益、斜距改正、拖鱼归位估算、地理编码、缝隙填补、消噪增强等环节。原始侧扫声呐数据经过上述各环节处理后，可最大程度消除各种畸变，并去除水柱区域，最终可形成一幅具有直观坐标信息的高质量无缝地貌条带图像。但受拖曳模式和拖鱼位置估算模型制约，精处理后的地貌图像仍存在系统位置偏差和局部畸变，从而造成相邻条带镶嵌后共视目标的错位和畸变(图 1)。
图 1 条带图像共视目标的畸变及错位
Figure 1. Distortion and Dislocation Between Corresponding Features in Side-Scan Sonar Mosaic
预处理后的侧扫声呐图像可有效地保证SURF特征点对提取的数量和质量，为开展侧扫声呐条带图像特征级镶嵌工作提供了条件。
-
侧扫声呐图像为固定宽度、长条状扫描图像，图像畸变大小沿航迹方向不断变化，该畸变机理决定相邻条带间公共区域的配准参数也随着航迹方向变化；另一方面，侧扫声呐图像地理编码后图幅较大，而相邻条带间仅有部分区域重叠，基于整幅图像进行特征点检测，计算量非常大。因此，相邻条带图像配准时，应沿航迹方向对重叠区域进行分区处理。
本文根据有误差的地理坐标信息，采用分块方法对重叠区域进行特征点检测。分块原则如图 2所示，相邻条带重叠区域可以由条带1的顶点d1、d2及条带2的顶点p1、p2共4个点确定，上述4点在纵、横向上长度为Lx、Ly。在图幅较长的方向，即图 1中的纵向，按等距离分为n段，各区块的横向宽度Δw可以通过航向a、扫幅S及设计重叠度b进行设置：
图 2 相邻条带图像重叠区域的分块
Figure 2. Block Scheme in Overlap Region of Adjacent Sonar Images
$$ \left\{ \begin{array}{l} \Delta h = \Delta w = S \times b\% /\left| {\cos a} \right|\\ n = {\mathop{\rm int}} \left( {{L_x}/\Delta h} \right) + 1\\ \Delta l = \left( {{L_y} - \Delta w} \right)/n \end{array} \right. $$ (1) 式中，纵向高度Δh取值可与Δw相等；Δl为相邻区块在横向上的平移间隔。
-
对所有区块提取的SURF特征点，按最近邻方法进行特征点匹配；然后，采用随机抽样一致性算法(random sample consensus, RANSAC)[9]剔除误匹配特征点对，设置合适的RANSAC模型阈值d，当一组特征点对的配准模型偏差小于d时，即认为其为正确的特征点对；最后，根据所有正确的特征点对，按最小二乘方法求解相邻条带图像间的刚性变换参数。为了计算方便，按式(2)确定相邻条带间的变换关系。
$$ \left( {\begin{array}{*{20}{c}} {{x_2}}\\ {{y_2}}\\ 1 \end{array}} \right) = \left[ {\begin{array}{*{20}{c}} a&b&c\\ { - b}&a&d\\ 0&0&1 \end{array}} \right] \times \left( {\begin{array}{*{20}{c}} {{x_1}}\\ {{y_1}}\\ 1 \end{array}} \right) $$ (2) 上述变换关系中共有4个未知数，其中a和b为旋转参数，c和d为平移参数。仅需2对特征点对(x1, y1)、(x2, y2)就可以解算出上述参数的初值t0=[a0 b0 c0 d0]T。当有n对特征点对时，可以得到关于各参数改正量的方程：
$$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{x_1}}&{{y_1}}&1&0\\ {{y_1}}&{ - {x_1}}&0&1\\ \vdots &{}&{}& \vdots \\ {{x_n}}&{{y_n}}&1&0\\ {{y_n}}&{ - {x_n}}&0&1 \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {\hat a}\\ {\hat b\hat c}\\ {\hat d} \end{array}} \right] = }\\ {\left[ {\begin{array}{*{20}{c}} {{{x'}_1}}\\ {{{y'}_1}}\\ \vdots \\ {{{x'}_n}}\\ {{{y'}_n}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{x_1}}&{{y_1}}&1&0\\ {{y_1}}&{ - {x_1}}&0&1\\ \vdots &{}&{}& \vdots \\ {{x_n}}&{{y_n}}&1&0\\ {{y_n}}&{ - {x_n}}&0&1 \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{a_0}}\\ {{b_0}}\\ {{c_0}}\\ {{d_0}} \end{array}} \right]} \end{array} $$ (3) 式(3)的一般形式为$ \mathit{\boldsymbol{B\hat t = l}}$，按最小二乘间接平差方法可以求得各参数的改正量$\mathit{\boldsymbol{\hat t}} $为：
$$ \mathit{\boldsymbol{\hat t}} = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}} \times \mathit{\boldsymbol{B}}} \right)^{ - 1}} \times {\mathit{\boldsymbol{B}}^{\rm{T}}} \times \mathit{\boldsymbol{l}} $$ (4) 对t0添加改正量$\mathit{\boldsymbol{\hat t}} $便可得到各变换参数。
-
受局部畸变的影响，声呐图像上局部目标呈现非刚性变化。若能根据有限控制点集P(x′i, y′i)和Q(xi, yi)(i=1, 2...n)构建一个随空间变化的光滑连续变换函数f(Q)，使得式(5)成立，就可以实现由Q到P的弹性配准[10]。
$$ \begin{array}{*{20}{c}} {E\left( f \right) = \sum\limits_{i = 1}^n {{{\left\| {{P_i} - f\left( {{Q_i}} \right)} \right\|}^2}} + } \\ {\lambda \iint {\left[ {{{\left( {\frac{{{\partial ^2}f}}{{\partial {x^2}}}} \right)}^2} + {{\left( {\frac{{{\partial ^2}f}}{{\partial x\partial y}}} \right)}^2} + {{\left( {\frac{{{\partial ^2}f}}{{\partial {y^2}}}} \right)}^2}} \right]{\text{d}}x{\text{d}}y} = \min } \end{array} $$ (5) 其中，λ为光滑参数，当λ=0时，式(5)的一般解形式为[11]：
$$ f\left( Q \right) = \sum\limits_{v = 1}^3 {{a_v}{\varphi _v}\left( Q \right)} + \sum\limits_{i = 1}^n {{w_i}{U_i}\left( {P,Q} \right)} $$ (6) 式中，第一项为关于Q的仿射变换，φ1(Q)=1，φ2(Q)=x，φ3(Q)=y；av为仿射变化系数；第二项基函数U(P, Q)为根据控制点个数张成的一个n维空间；wi为权值系数。为了使弹性配准后的光滑曲面弯曲能量最小，可选择具有该性质的薄板样条函数作为基函数[12]，其径向基函数为U(r)=r2lgr2，r为任意Pi和Qi之间的欧氏距离，当r为零时，U(r)=0。式(6)中有n+3个未知数，n对特征点对仅能列出n个方程，要求解各未知数还需要添加3个边界条件：
$$ \sum\limits_{i = 1}^n {{w_i} = 0} ,\sum\limits_{i = 1}^n {{x_i}{w_i} = 0} ,\sum\limits_{i = 1}^n {{y_i}{w_i} = 0} $$ (7) 上述边界条件表明，当x、y→∞时，式(6)趋向为一般的仿射变换。若分别用Fi和Gi代替x和y方向的基函数权值系数wi，根据P和Q的n对特征点，结合式(6)和(7)，可以列出方程组：
$$ \left[ {\begin{array}{*{20}{c}} 1&{{x_1}}&{{y_1}}&0&{U\left( {{r_{12}}} \right)}& \cdots &{U\left( {{r_{1n}}} \right)}\\ 1&{{x_2}}&{{y_2}}&{U\left( {{r_{21}}} \right)}&0& \cdots &{U\left( {{r_{2n}}} \right)}\\ \vdots &{}&{}&{}& \vdots&\ddots&\vdots \\ 1&{{x_n}}&{{y_n}}&{U\left( {{r_{n1}}} \right)}&{U\left( {{r_{n2}}} \right)}& \cdots &0\\ 0&0&0&1&1& \cdots &1\\ 0&0&0&{{x_1}}&{{x_2}}&{}&{{x_n}}\\ 0&0&0&{{y_1}}&{{y_2}}& \cdots &{{y_n}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{a_0}}\\ {{a_1}}\\ {{a_2}}\\ {{F_1}}\\ {{F_2}}\\ \vdots \\ {{F_n}} \end{array}} \right]\left( x \right) = \left[ {\begin{array}{*{20}{c}} {{{x'}_1}}\\ {{{x'}_2}}\\ \vdots \\ {{{x'}_n}}\\ 0\\ 0\\ 0 \end{array}} \right] $$ (8) $$ \left[ {\begin{array}{*{20}{c}} 1&{{x_1}}&{{y_1}}&0&{U\left( {{r_{12}}} \right)}& \cdots &{U\left( {{r_{1n}}} \right)}\\ 1&{{x_2}}&{{y_2}}&{U\left( {{r_{21}}} \right)}&0& \cdots &{U\left( {{r_{2n}}} \right)}\\ \vdots &{}&{}&{}& \vdots&\ddots&\vdots \\ 1&{{x_n}}&{{y_n}}&{U\left( {{r_{n1}}} \right)}&{U\left( {{r_{n2}}} \right)}& \cdots &0\\ 0&0&0&1&1& \cdots &1\\ 0&0&0&{{x_1}}&{{x_2}}&{}&{{x_n}}\\ 0&0&0&{{y_1}}&{{y_2}}& \cdots &{{y_n}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{b_0}}\\ {{b_1}}\\ {{b_2}}\\ {{G_1}}\\ {{G_2}}\\ \vdots \\ {{G_n}} \end{array}} \right]\left( y \right) = \left[ {\begin{array}{*{20}{c}} {{{y'}_1}}\\ {{{y'}_2}}\\ \vdots \\ {{{y'}_n}}\\ 0\\ 0\\ 0 \end{array}} \right] $$ (9) 式(8)和(9)可以分别解算出特征点对在x和y方向的弹性变换参数[a0 a1 a2 F1 F2...Fn]T和[b0 b1 b2 G1 G2...Gn]T。根据上述参数，Q空间上的任意点(x, y)可以按式(10)弹性变换至P空间上。
$$ \left\{ \begin{array}{l} x' = {a_0} + {a_1}x + {a_2}y + \sum\limits_{i = 1}^n {{F_i}r_i^2\lg r_i^2} \\ y' = {b_0} + {b_1}x + {b_2}y + \sum\limits_{i = 1}^n {{G_i}r_i^2\lg r_i^2} \end{array} \right. $$ (10) 式(8)和(9)的一般形式可写为：
$$ \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{P}}}\limits_{n \times 3} }&{\mathop {\mathit{\boldsymbol{U}}}\limits_{n \times n} }\\ {\mathop {\mathit{\boldsymbol{0}}}\limits_{3 \times 3} }&{{{\mathop {\mathit{\boldsymbol{P}}^{\rm{T}}}\limits_{3 \times n} }}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{A}}}\limits_{3 \times 1} }\\ {\mathop {\mathit{\boldsymbol{w}}}\limits_{n \times 1} } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{Q}}}\limits_{n \times 1} }\\ {\mathop {\mathit{\boldsymbol{\mathit{\pmb{0}}}}}\limits_{3 \times 1} } \end{array}} \right] $$ (11) 由于式(5)中的λ=0，式(11)矩阵块U的对角元素均为零，薄板样条弯曲能量最小，薄板函数严格经过各参考点。当不严格要求薄板函数过形值点时，矩阵U为：
$$ \mathit{\boldsymbol{U}} = \mathit{\boldsymbol{U}} + \lambda \mathit{\boldsymbol{I}} $$ (12) 式中，I为单位向量；λ的值决定了薄板样条函数的局部弹性变形大小，当λ→∞时，式(11)就变为一般的仿射变换。
薄板样条的上述网格弹性变化可以用于解决侧扫声呐图像上目标的变形及扭曲问题。本文采用RANSAC设定阈值d来代替式(12)中的λ，以此来减弱特征点对匹配误差对局部弹性配准的影响；另一方面，添加一系列边界点条件，并且令矩阵U相应对角线元素为零，从而保证变换前后这些点的位置不发生变化，以此来解决区域的接边问题。上述处理后，式(12)将变为：
$$ \mathit{\boldsymbol{U}} = \left[ {\begin{array}{*{20}{c}} {\mathop {{\mathit{\boldsymbol{U}}_b}}\limits_{m \times m} }&{\mathop {{\mathit{\boldsymbol{U}}_{b0}}}\limits_{m \times n} }\\ {\mathop {{\mathit{\boldsymbol{U}}_{0b}}}\limits_{n \times m} }&{\mathop {{\mathit{\boldsymbol{U}}_0}}\limits_{n \times n} } \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\pmb{0}}}\limits_{m \times m} }&{\mathop {\mathit{\pmb{0}}}\limits_{m \times n} }\\ {\mathop {\mathit{\pmb{0}}}\limits_{n \times m} }&{{\rm{d}}\mathit{\boldsymbol{I}}} \end{array}} \right] $$ (13) 式(13)中，方阵U对角线元素的大小决定了相关特征点对网格变形的贡献大小。当U对角线上第i个元素为零时，薄板样条函数将严格通过第i个特征点；当第i个元素逐渐增大时，第i个特征点对网格变形贡献量逐渐减小。图 3为采用薄板样条这种特性来处理边界情况前后的效果对比。
图 3 薄板样条局部弹性变换的边界处理
Figure 3. Boundary Treatment in Local Elastic Transformation
-
为了检验本文方法的性能，本文选取了两幅具有重叠区域的相邻侧扫声呐条带图像进行镶嵌实验。所选取的两幅条带图已经过预处理，消除了大多数畸变，并具有坐标信息，每个条带的扫幅为400 m，像素尺寸为0.6 m，条带间设计重叠度为50%，其他参数如表 1所示。
表 1 声呐条带图像的相关参数
Table 1. Relative Parameters of the Sonar Images
图像 参数 数值 扫幅/m 400 像素尺寸/m 0.6 条带1 测线长度/m 2 382.862 Ping数 3 889 图幅(宽×高)/(m×m) 1 480.8×2 320.2 分辨率(宽×高)/像素 2 468×3 864 条带2 测线长度/m 2 468.524 Ping数 3 643 图幅(宽×高)/(m×m) 1 519.2×2 433.6 分辨率(宽×高)/像素 2 532×4 056 相邻条带重叠度(%) 50 根据地理坐标信息，对上述两个相邻的条带图像的重叠区域进行分块处理。两个条带重叠区域宽度约为200 m，分块边长为300 m，区块大小为500×500像素，采用SURF算法对每个区块图像进行特征点检测，图 4(a)、4(b)分别为两个条带分块及特征点检测情况。
图 4 相邻条带重叠区域的分块及SURF特征点检测
Figure 4. Results of SURF Feature Extraction by Block in Overlap Region
按最近邻法提取条带1与条带2的特征点对，并采用RANSAC算法剔除误匹配特征点对。本文设定变换模型的阈值d=15像素，特征检测及配对统计情况如表 2所示。可以看出，采用重叠区域分块特征检测，忽略了重叠区域以外的非相关区域，从而节约了特征点检测时间，减少了错误特征点对的数量，提高了RANSAC提取正确特征点对的效率。
表 2 SURF特征点检测及匹配相关参数统计表
Table 2. Statistical Parameters of SURF Feature Extraction and Registration
统计项目 全局 分块 CPU时间/s 1 304.261 258.14 条带1特征点数 2 778 1 274 条带2特征点数 3 305 1 444 RANSAC配对个数 178 174 根据上述所有特征点对，按最小二乘方法求解条带1和2之间的刚性变换参数，并采用加权融合法得到两个条带的镶嵌图像，如图 5(b)所示。与直接根据坐标信息得到的镶嵌图像(图 5(a))进行比较，基于SURF特征点对的刚性变换方法消除了相邻条带图像之间的整体位置及角度偏差，镶嵌图像上基本上看不出共视目标的错位。
图 5 不同配准方法的镶嵌效果
Figure 5. Mosaic Results by Different Registration Methods
然而受侧扫声呐图像局部畸变的影响，经刚性变换并融合后，仍有部分共视目标的轮廓模糊，本文根据特征点的分布情况，在各区块内采用薄板样条对共视目标进行局部弹性配准。如果各区块特征点对中没有误匹配情况，薄板样条弹性配准可以很好地消除局部目标的扭曲变形。为了清楚地观察不同配准方法的效果，从图 5中不同方法形成的镶嵌图像上，分别选取了3块对应区域的区块图像放大显示，如图 6所示。
图 6 不同配准方法的局部区域效果比较
Figure 6. Results of Local Region by Different Registration Methods
从图 6(b)中可以看出，刚性变换难以实现特征点对之间的精确对准，局部畸变引起共视目标融合后轮廓模糊；从图 6(c)中可以看出，弹性变换使各特征点对精确对准，局部的变形得到较好的修正，融合后的共视目标轮廓清晰。
根据特征点对集合配准前后的坐标信息，可对各镶嵌方法进行定量分析，如表 3所示，主要统计参数包括坐标偏差的最大绝对值、坐标偏差的平均值、坐标偏差的标准差、特征点对标准点位差。综合图 7及表 3可以看出：
图 7 采用不同配准方法时特征点对的坐标偏差分布
Figure 7. Coordinate Error Distribution of Feature Points by Different Registration Methods
表 3 不同配准方法相关参数统计
Table 3. Statistical Parameters of Different Registration Methods
配准方法 偏差类型 最大绝对值/m 平均值/m 标准差/m 标准点位差/m 计算公式 横向偏差dE $\max \left\{ {|{x_i}-{x_i}'|, i = 1, 2 \cdots n} \right\} $ $ {\rm{d}}\bar x = \frac{1}{n}\sum\limits_{i = 1}^n {\left( {{x_i}-{x_i}'} \right)} $ $ {\sigma _x} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {\left( {{x_i}-{x_i}'-{\rm{d}}\bar x} \right)} } $ $\sigma = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {\left\{ {{{\left( {{x_i}-{x_i}'} \right)}^2} + {{\left( {{y_i}-{y_i}'} \right)}^2}} \right\}} } $ 纵向偏差dN $ \max \left\{ {|{y_i}-{y_i}'|, i = 1, 2 \cdots n} \right\}$ $ {\rm{d}}\bar y = \frac{1}{n}\sum\limits_{i = 1}^n {\left( {{y_i}-{y_i}'} \right)} $ ${\sigma _y} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {\left( {{y_i}-{y_i}'-{\rm{d}}\bar y} \right)} } $ 坐标镶嵌 横向偏差d 10.79 -2.19 2.69 8.91 纵向偏差dN 14.83 -7.76 2.61 刚性变换 横向偏差dE 6.27 0 2.69 3.74 纵向偏差dN 7.21 0 2.61 弹性变换 横向偏差dE 6.59 0 1.11 1.21 纵向偏差dN 4.92 0 1.17 1) 实验所选相邻条带之间存在较大的位置偏差，这种偏差很大一部分是由系统误差引起的。由于测量要求不高，实测时GPS至拖曳点之间的位置关系及拖缆长度估算并不准确。当船往返测量时，这种系统误差引起的偏差将在航向方向上表现更为明显。表 3中，特征点对的原始地理坐标在北方向(north, N)上的平均偏差达到了-7.76 m，东方向(east, E)上为-2.19 m。
2) 由于船速变化、拖鱼姿态不稳定，地理编码后的声呐图像各局部区域仍有残留畸变，各特征点对东北方向上原始坐标偏差的标准差分别为2.69 m和2.61 m。根据SURF特征点对进行条带间的整体刚性变换，可以很好地消除条带间的系统误差，特征点对之间的标准点位差由原来的8.91 m降为3.74 m，但刚性变换难以消除。由于拖鱼瞬时运动状态变化引起的局部畸变，各特征点对东北方向上的坐标偏差标准差仍为2.69 m和2.61 m。
3) 经过各区块的薄板样条弹性变换，各特征点对的对准精度相较于刚性变换有了很大的提高，特征点的坐标偏差标准差提高至1.11 m和1.17 m，特征点对之间的标准点位差也达到了1.21 m(2像素)，极大地消弱了局部残留畸变对侧扫声呐图像镶嵌的影响，得到了较好的配准结果。
-
本文针对侧扫声呐拖曳作业模式引起的相邻条带图像间的系统误差和局部畸变，基于SURF及薄板样条函数，提出了一种适用于侧扫声呐条带图像的分块弹性镶嵌方法。该方法首先基于坐标信息，在相邻条带重叠区域进行分块SURF提取；然后，基于所有特征点对对相邻条带图像进行整体刚性配准；最后，采用薄板样条在各区块内实现特征点对弹性配准。实验结果表明该方法减少了侧扫条带图像特征提取时间，削弱了系统误差和局部畸变，实现了公共区域目标的保形性，显著提高了侧扫声呐镶嵌图像的质量。
Elastic Mosaic Method in Block for Side-Scan Sonar Image Based on Speeded-Up Robust Features
-
摘要: 拖曳式侧扫声呐图像沿航迹方向存在类型和大小不断变化的局部畸变，造成条带图像镶嵌时共视目标扭曲和错位。为了实现公共区域融合后目标的保形性，提出了一种基于加速稳健特征（speeded-up robust features，SURF）的分块弹性镶嵌方法。该方法通过对相邻条带重叠区域自动分块提取特征、整体刚性变换、局部薄板样条弹性变换等处理，既可消除相邻条带图像间的系统性位置和航向偏差，又能消弱各处随机性局部畸变，从而实现共视目标的特征级保形镶嵌。实验结果表明，特征点对整体配准精度达到2个像素，验证了该方法的有效性。Abstract: The side-scan sonar towed working mode leads to changing local distortion in types and magnitude along the track, and further results in features' disposition and distortion in mosaic image. In order to guarantee clarity edge of the local features in overlap region of adjacent images, this paper proposes an elastic blocking mosaic method based on speeded-up robust features. Firstly, combining with the information of track line and swath, each stripe image can be sliced to several blocks. Then the SURF image registration is done in each block including features extraction, matching and refining. With the registration point-pairs between blocks from adjacent stripe image, a rigid transformation in global image and a series of elastic transformation in each block image can be done. The method not only eliminates systematic errors of position and heading between two adjacent side-scan sonar images, but also weakens random local distortion and achieves high-accurate registration in local area by corresponding features. At last, the experiment demonstrates the validity of this method, the overall registration accuracy reaches two pixels.
-
Key words:
- side-scan sonar image /
- SURF match /
- block mosaic /
- elastic mosaic
-
图 1 条带图像共视目标的畸变及错位
Figure 1. Distortion and Dislocation Between Corresponding Features in Side-Scan Sonar Mosaic
图 4 相邻条带重叠区域的分块及SURF特征点检测
Figure 4. Results of SURF Feature Extraction by Block in Overlap Region
图 6 不同配准方法的局部区域效果比较
Figure 6. Results of Local Region by Different Registration Methods
图 7 采用不同配准方法时特征点对的坐标偏差分布
Figure 7. Coordinate Error Distribution of Feature Points by Different Registration Methods
表 1 声呐条带图像的相关参数
Table 1. Relative Parameters of the Sonar Images
图像 参数 数值 扫幅/m 400 像素尺寸/m 0.6 条带1 测线长度/m 2 382.862 Ping数 3 889 图幅(宽×高)/(m×m) 1 480.8×2 320.2 分辨率(宽×高)/像素 2 468×3 864 条带2 测线长度/m 2 468.524 Ping数 3 643 图幅(宽×高)/(m×m) 1 519.2×2 433.6 分辨率(宽×高)/像素 2 532×4 056 相邻条带重叠度(%) 50 
下载: 导出CSV
表 2 SURF特征点检测及匹配相关参数统计表
Table 2. Statistical Parameters of SURF Feature Extraction and Registration
统计项目 全局 分块 CPU时间/s 1 304.261 258.14 条带1特征点数 2 778 1 274 条带2特征点数 3 305 1 444 RANSAC配对个数 178 174
下载: 导出CSV
表 3 不同配准方法相关参数统计
Table 3. Statistical Parameters of Different Registration Methods
配准方法 偏差类型 最大绝对值/m 平均值/m 标准差/m 标准点位差/m 计算公式 横向偏差dE $\max \left\{ {|{x_i}-{x_i}'|, i = 1, 2 \cdots n} \right\} $ $ {\rm{d}}\bar x = \frac{1}{n}\sum\limits_{i = 1}^n {\left( {{x_i}-{x_i}'} \right)} $ $ {\sigma _x} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {\left( {{x_i}-{x_i}'-{\rm{d}}\bar x} \right)} } $ $\sigma = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {\left\{ {{{\left( {{x_i}-{x_i}'} \right)}^2} + {{\left( {{y_i}-{y_i}'} \right)}^2}} \right\}} } $ 纵向偏差dN $ \max \left\{ {|{y_i}-{y_i}'|, i = 1, 2 \cdots n} \right\}$ $ {\rm{d}}\bar y = \frac{1}{n}\sum\limits_{i = 1}^n {\left( {{y_i}-{y_i}'} \right)} $ ${\sigma _y} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {\left( {{y_i}-{y_i}'-{\rm{d}}\bar y} \right)} } $ 坐标镶嵌 横向偏差d 10.79 -2.19 2.69 8.91 纵向偏差dN 14.83 -7.76 2.61 刚性变换 横向偏差dE 6.27 0 2.69 3.74 纵向偏差dN 7.21 0 2.61 弹性变换 横向偏差dE 6.59 0 1.11 1.21 纵向偏差dN 4.92 0 1.17
下载: 导出CSV
-
[1] 赵建虎, 刘经南.多波束测深及图像数据处理[M].武汉:武汉大学出版社, 2008 Zhao Jianhu, Liu Jingnan. Principles, Methods and Applications of Sounding and Sonar Image Data Processing of Multibeam Bathmetric System[M]. Wuhan:Wuhan University Press, 2008 [2] Zhao Jian, Wang Aixue, Zhang Hongmei, et al. Mosaic Method of Side-scan Sonar Strip Images Using Corresponding Features[J]. IET Image Processing, 2013, 7(6):616-623 doi: 10.1049/iet-ipr.2012.0468 [3] Vandrish P, Vardy A, Walker D, et al. Side-Scan Sonar Image Registration for AUV Navigation[C]. IEEE Symposium Underwater Technology (UT) and 2011 Workshop on Scientific Use of Submarine Cables and Related Technologies (SSC), Tokyo, Japan, 2011 http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=5774096 [4] Daniel S, Leannec F L, Roux C, et al. Side-scan Sonar Image Matching[J]. IEEE Journal of Oceanic Engineering, 1998, 23(3):245-259 doi: 10.1109/48.701197 [5] Chailloux C, Caillec J M L, Gueriot D, et al. Intensity-Based Block Matching Algorithm for Mosaicing Sonar Images[J]. IEEE Journal of Oceanic Engineering, 2011, 36(4):627-645 doi: 10.1109/JOE.2011.2141850 [6] Chailloux C, Zerr B. Non Symbolic Methods to Register SONAR Images[C]. IEEE OCEAN'05 Europe, Brest, France, 2005 http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=1511725 [7] 王爱学. 基于侧扫声纳图像的三维海床地形恢复[D]. 武汉: 武汉大学, 2014 Wang Aixue. Research on 3D Seafloor Terrian Recovery from the Side Scan Sonar Image[D]. Wuhan: Wuhan University, 2014 [8] 赵建虎, 王爱学, 王晓, 等.侧扫声纳条带图像分段拼接方法研究[J].武汉大学学报·信息科学版, 2013, 38(9):1034-1038 http://ch.whu.edu.cn/CN/abstract/abstract2746.shtml Zhao Jianhu, Wang Aixue, Wang Xiao, et al. A Segmented Mosaic Method for Side Scan Sonar Strip Images Using Corresponding Features[J]. Geomatics and Information Science of Wuhan University, 2013, 38(9):1034-1038 http://ch.whu.edu.cn/CN/abstract/abstract2746.shtml [9] Chum O. Two-View Geometry Estimation by Random Sample and Consensus[D]. Prague: Czech Technical University, 2005 https://www.researchgate.net/publication/240366623_Two-view_geometry_estimation_by_random_sample_and_consensus [10] Galway L. Spline Models for Observational Data[M]. Siam:Society for Industrial and Applied Mathematics, 1990 [11] Sprengel R, Rohr K, Stiehl H S. Thin-Plate Spline Approximation for Image Registration[C]. The 18th Annual International Conference on Bridging Disciplines for Biomedicine, Amsterdam, Netherlands, 1996 http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=652767 [12] Turk G, O'brien J F. Shape Transformation Using Variational Implicit Functions[C]. Annual Conference of Computer Graphics, New York, USA, 1999 http://dl.acm.org/citation.cfm?id=311580 -
点击查看大图
计量
- 文章访问数: 1768
- HTML全文浏览量: 102
- PDF下载量: 281
- 被引次数: 0
