一种基于水下机器人的构筑物场景复原优化方法

郭云翔，周 军

（河海大学 机电工程学院，江苏 常州 213022）

一种基于水下机器人的构筑物场景复原优化方法

郭云翔，周 军

（河海大学 机电工程学院，江苏 常州 213022）

为了克服水下构筑物场景复原过程中的累积误差，有效还原水下构筑物的表面结构信息，本文提出了一种基于水下机器人的构筑物场景复原优化方法。通过构建最大生成树，采用基本的图算法寻找基准坐标系以及完成坐标系的转换，并建立了误差方程，选取L-M迭代算法实现变换矩阵的最优化，采用对比实验实现了对该方法的可靠性检测。实验结果表明，采用该优化方法后可有效反映整个探测过程的图像信息，为后续的机器视觉处理建立了有效基础。

水下构筑物；累积误差；优化方法；最大生成树；变换矩阵

随着水下机器人的发展，传统的人工潜水作业已逐渐被水下机器人替代[1]。通过搭载在水下机器人上的各种传感器获得大量的数据，包括水深、机器人姿态等[2]，而视觉信息是机器人实现正确导航以及水下近距离作业的重要保障，是研究人员对水下环境以及水下构筑物研究的重要参考信息[3]。而视觉信息作为整个机器视觉系统的一个重要组成部分，有效的视觉信息可以提高后续判决执行部分的可靠性。

水下作业的重要一环为水下构筑物的探测任务，包括大坝裂缝、闸门漏洞、管道泄漏等[4]。通过对水下构筑物的场景重构，可以对水下构筑物表面的整体结构与细节部分有直接的对比观察，能更加清晰地判别出是否有安全隐患。场景重构技术的基本思想为：将零碎的细小的场景组合重构为一个完整的大的场景，包括二维场景重构和三围场景重构两部分[5-6]。二维场景重构中的图像拼接技术是将构筑物表面视频图像拼接复原为一副大的全景图像。

全景图拼接理念主要分为基于相位和特征相关两类图像拼接方法[7]。由于相位相关法具有一定的局限性，基于特征相关的图像拼接方法是当前研究的重点方向。图像拼接基本流程为：图像匹配、图像变换、图像融合。良好的视觉信息处理有利于更进一步的机器视觉方面的研究。

1 基于视频图像的水下构筑物场景重建优化方法

1.1 图像匹配

首先提取图像SURF特征点，构建k-d树建立特征向量索引树，采用BBF搜索算法实现特征点对的匹配，并通过RANSAC算法完成匹配图像之间的单应矩阵计算。

1.2 图像变换

1.2.1 单应矩阵

在计算机视觉领域，平面的单应性被定义为一个平面到另一个平面的投影映射。因此，常用单应矩阵将两幅图像变换到同一平面坐标系下。假设点二维平面上一点q（x,y）变换后为Q（x,y），利用齐次坐标系，由单应矩阵变换有：

式中：s——尺度因子；

H——单应矩阵；

m0，m1，m3，m4——旋转角度；

m2，m5——分别表示水平、垂直方向上的位移；

m6，m7——分别表示水平、垂直方向上的形变程度。

通常归一化尺度后，计算式（1）两边有：

1.2.2 最大生成树统一坐标系

在完成n副图像的拼接时，需要统一坐标变换矩阵，寻找的基准坐标系在具有稳定性的同时还应该减少变换过程中的运算量、累积误差等。可通过选择合适的准则构建最大生成树，采取广度优先搜索来寻找基准坐标系，采用深度优先搜索来完成各坐标系变换路径的搜索[8]。

最大生成树（Maximum Spanning Tree）：在无向连通图G=（V，E）中，选取一个边集E′使得集合V中所有点被连接，且E′中边权值和最小，其中E′∈E。本文构建最大生成树选择以下几条作为准则：

①以待拼接的n副图像作为顶点组成顶点集V；

②以两幅图像之间是否匹配来判别顶点之间是否存在边e；

③以两幅图像之间匹配的特征点对数目作为边的权重值w；

Kruskal算法常用来寻找最小生成树，是贪心算法的一种，相对应的也可以用来寻找最大生成树。此时，生成的最大生成树记为：G=（V，Emax）。整个算法步骤如下：

Step1：初始化图G=（V，E），输入镜头内所有的关键帧图像生成顶点集V，输入由所有关键帧之间匹配关系构成的边生成边集E；

Step2：新建图Gmax=（V，Emax），其中Emax=NULL，此时Gmax中所有顶点都未连接；

Step3：对E中所有边按权重值进行降序排列，排列后边集记为Esort；

Step4：loop：遍历Esort中所有边

广度优先搜索（Breadth First Search）为基本图搜索算法的一种，算法步骤如下：

Step1：初始化输入max_span_tree，叶节点集合span_tree_leafs；初始化集合max_dist（n，0），n为顶点数目，0为初始化值；

Step2：loop1：遍历span_tree_leafs中所有叶节点由叶节点开始广度优先搜索，搜索结果记录在cur_dist中；

Loop2：遍历cur_dist和max_dist中元素，当前元素记为I，以cur_dist[i]和max_dist[i]中较大的元素更新max_dist[i]；

Step3：遍历max_dist中所有元素，选择其中最小的元素记为min_max_dist；

Step4：记max_dist中所有等于min_max_dist的元素为候选基准坐标系集合，若集合中候选基准坐标系大于1，则选择其中邻接点最多的元素作为基准坐标系。

如图1所示为采用Kruskal算法寻找的最大生成树，左图为输入镜头内所有关键帧图像之间的匹配关系，右图为寻找的最大生成树。

图1 Kruskal算法构建最大生成树过程

以图1为例，采用广度优先搜索Step2中描述方法，由每个叶节点进行广度优先搜索，各叶节点cur_dist记录值数据如下：

表1 cur_dist中记录数据

比较后可得到到max_dist数据由a至f分别为：2、2、3、3、3、3，显而易见max_dist中最小的值为2，此时候选的基准坐标系集合为：a、b。由Step4中描述方法，最后选择b为基准坐标系。简单的验证可知：

选a时，变换总次数：1+2+2+2=7次；

选b时，变换总次数：1+2+1+1=5次；

由于深度优先搜索（Breadth First Search）的目的为记录从其他坐标系变换到基准坐标系的路径，因此采用改进的深度优先搜索方法：从基准坐标系以外的点出发，依次递归深度优先访问它的邻接点，记录搜索过程中的每一步路径，已访问的不再记录，直到搜索到基准坐标系为止。

1.2.3 基于L-M算法最优化变换矩阵

若直接将输入的n幅图像变换到同一坐标系下会发现n幅图像之间的配准误差较大，在图像拼接缝隙处有明显的错位现象。这是由于，单个变换矩阵只是图像两两匹配之间的结果，尽管其局部误差较小，对于n幅图像的统一变换而言，其累计误差较大，从而影响变换后结果。

以1.2.2节中统一坐标变换系的过程中产生的累积误差建立数学模型，即误差方程。如图2所示，将B变换到C时，先将B变换到A，再变换到C。假设：（x1，y1）、（x2，y2）、（x3，y3）分别为A、B、C中提取的特征点，其中（x1，y1）->（x2，y2）匹配，点（x1，y1）->（x3，y3）匹配。由1.2.1节中式（1）计算变换后点对之间的误差可得：

图2 变换误差示意图

式中：r为变换误差，Hba为图像B变换到A的单应矩阵，Hca为图像C变换到A的单应矩阵。

因此，对于n幅图像变换，计算其所有匹配关系以及所有匹配点对之间的累积误差模型建立如下：

式中：E——累积误差；

i——第i副图像；

n——输入带拼接的图像总数；

g（i）——所有与图像i匹配的集合；

j——集合g（i）中第j副与i匹配的图像；

m（i，j）——图像i、j匹配的特征点对集合；

k——第k对匹配的特征点；

Hi、Hj——图像i、j变换到基准坐标系的单应矩阵。

最优化问题为：求解变量{H1，H2，……，Hn}，使得累积误差E最小。由1.2.1节中式（1）知单应矩阵为8参数变换矩阵，考虑到各图像变换尺度不同，对于n副图像未知参数最多为9×n。又由式知（2）计算一对匹配点误差，本质为计算其x、y方向上的误差，即可由两个未知方程组成的方程组f表示。因此，对于所有匹配点对误差方程组合为向量函数f：R2×K→R9×n。最优化问题转为求解：

式中：f：R2×K→R9×n，2×K＞9×n；m为n个单应矩阵未知参数组成的向量，最多为9×n；K为所有匹配点对总数。

迭代法为解非线性最小二乘问题的基本方法，迭代条件为：f（xk+1）＜f（xk），也称为下降法。常用的有：最速下降法、牛顿法、高斯-牛顿法、Levenberg-Marquardt算法。

将向量函数f：R2×K→R9×n，在x处泰勒展开[9]有：

式中：迭代步长h=[h1，h2，……，h9×n]，J（x）为f（x）的雅克比矩阵

则式（7）中F（m）在x处的梯度为：

F（m）在x处的Hessian矩阵为：

将式（7）在x处泰勒展开有：

分析可知，已知x，可求得矩阵f（x）、J（x）分别记为f、J。此时式（12）可看作以h为变量的向量函数L（h）：

L（h）的梯度与Hessian矩阵分别为：

高斯-牛顿法以L（h）的极小值点hgn为步长进行迭代，提取SRUF特征点时采用Hessian矩阵[10]来检测极值点，其较小值判定条件为：△f（m）=0，且H（M）是正定矩阵。该方法同样适用于检测L（h）的极值点，由式（14）知L″（h）只与J相关，若J为满秩矩阵，则L″（h）是正定矩阵，因此对高斯-牛顿法步长hgn求解得：

高斯-牛顿迭代法公式有：

Levenberg-Marquardt算法是高斯-牛顿法的改进算法[11]，Levenberg-Marquardt算法中加入阻尼因子[12]来约束迭代步长，其迭代公式为：

式中：I为单位矩阵，μ＞0。

μ初始化与A0=J（x0）TJ（x0）中各元素值的范围有关，常取：

式中：τ由用户选择，与初始化点相关，在本文中初始化点x0为统一变换矩阵后的各单应矩阵参数组成的向量。在1.1节中采用RANSAC算法计算得到的各单应矩阵具有一定的精度，因此本文取τ=1。

μ随着迭代次数而更新，由式（12）与式（13）的变化程度的比值θ确定：

②当μ很小时，hlm≈hgn，适用于当前跌打参数离最优值较近的情况；

迭代终止条件：g=F′（x*）=J（x*）Tf（x*）=0，然而实际应用中往往不符合这一要求，因此迭代终止条件常用：||g||∞≤□1，其中□1是足够小的正数，有用户决定的。另外，为了避免无限循环，还应设置最大迭代次数，当迭代次数超过最大迭代次数时终止迭代。

Levenberg-Marquardt算法伪码如下：

2 实验结果与分析

2.1 图像匹配

采用1.1节中基于SURF特征点的图像匹配算法，快速准确的实现了镜头中提取的9幅关键帧图像之间的匹配[13]。实验场景图如图3所示。

图3 实验场景图

2.2 图像变换

以两幅图像匹配的特征点对数目为权重值构建的最大生成树[14]。如图4所示，有：

图4 镜头中关键帧图像构成的最大生成树

编号2、3最大变换次数都为3，但此时编号2总变换次数：1×4+2×3+3=13，编号3总变换次数：1× 2+2×4+3×2=16。

由基准坐标系寻找规则，编号2、3最大变换次数相同，但是编号2总变换次数较少且邻接点较多，因此选择编号2为基准坐标系。最大生成树树中边的各参数值如表2所示，weight为两幅图像匹配的特征点对数目。图5为深度搜索统一变换坐标系，以及统一后的单应矩阵，push后编号表示当前图像到基准坐标系的变换路径。

表2 最大生成树中各边对应相关值

图5中homography is后面矩阵为未优化之前的各图像到编号2图像的单应矩阵值，直接变换后结果如图6a所示，有着明显的误差存在。

以各图像变换模型为例，迭代次数与迭代时间实验结果如表3所示。

图5 深度优先搜索统一坐标系变换路径

图6 拼接效果图

表3 各算法迭代次数与时间实验结果

由表3可以看出高斯-牛顿法与L-M算法在400次范围内迭代收敛，最速下降法在400次迭代范围内不能达到收敛条件，且L-M算法较高斯-牛顿法有较少的迭代次数和耗时。L-M算法综合了高斯-牛顿法迭代精度高与最速下降法迭代速度快的优点[15]，且L-M算法较高斯-牛顿法有较少的迭代次数和耗时[16]。因此，选取L-M算法可以优化统一变换矩阵过程中产生的累积误差。

由1.2.3节中L-M算法迭代优化后的参数如表4所示，total num为所有匹配点对数目，且对应式中K的值。

表4 L-M迭代优化后的参数

2.3 图像融合

L-M优化后，采用4层拉普拉斯金字塔融合方法结果如图6b所示。对比图6a、6b可以看出，在图6b中n副图像之间的配准误差已经较小，在图像拼接缝隙处的错位现象也得到了一定的改善，且过渡平滑。

除了对于图像的视觉效果的评价外，本文分别对于融合图像的信息熵、平均梯度和均方根误差作为评价标准对融合后的图像进行评估如表5所示。

表5 优化前后的评价指标

由表5可知，优化融合后的图像的均方根误差小于未优化的图像，均方根误差越小，匹配误差越小，融合效果越好。且优化融合后的图像信息熵和平均梯度都增大，融合后的图像拥有更多的信息量。对比优化前后，对于n副图像的统一变换而言，其累计误差得以减小，从而优化了水下构筑物场景复原的结果。

3 结束语

文章介绍了一种基于水下机器人的构筑物场景复原优化方法，并进行了测试。实验结果表明，该变换方法可以很好地减少n幅图像拼接过程中的累计误差。最后的水下构筑物的复原结果能够较好的反映整个探测过程的图像信息，为后续的机器视觉处理奠定了可靠的基础。

[1]徐玉如，李彭超.水下机器人发展趋势[J].自然杂志，2011，33（3）：125-132.

[2]唐旭东，庞永杰.基于单目视觉的水下机器人管道检测[J].机器人，2010，32（5）：592-600.

[3]吴国帅.观测型ROV视觉系统设计［D］.青岛：中国海洋大学，2014.

[4]张洪欣，马 龙.水下机器人在海洋观测领域的应用进展[J].遥测遥控，2015，36（5）：23-27.

[5]赵次郎.基于激光视觉数据融合的三维场景重构与监控[D].大连：大连理工大学，2014.

[6]赵魏钰.基于压缩感知的图像超分辨率重构算法研究[D].吉林：吉林大学，2014.

[7]王 娟，师 军，吴宪祥.图像拼接技术综述[J].计算机应用研究，2008，25（7）：1940-1943.

[8]Thomas H.cormen.算法导论[M].北京：机械工业出版社，2013-01.

[9]Richard Szeliski.Computer Vision：Algorithms and Applications[M]. Berlin：Springer，2010.

[10]张 桦.场景图像拼接关键技术研究[D].天津：天津大学，2008.

[11]高 隽，谢 昭.图像理解理论与方法[M].北京：科学出版社，2009.

[12]李珊珊.计算机视觉中特征与相似性度量研究[D].合肥：中国科学技术大学，2010.

[13]李 珅，廖华丽.基于虚幻引擎的ROV实时运动仿真[J].计算机与现代化，2014，（7）：150-153.

[14]Lowe，D.G.Distinctive Image Featuresfrom Scale-Invariant Keypoints[J].International Journal of Computer Vison，2004.60（2）：91-110.

[15]David G.Lowe.Distinctive Image Features from Scale-Invariant Keypoints[J].International Journal of Computer Vision，2004.60（2）：91-110.

[16]王 娟，师 军，吴宪祥.图像拼接技术综述[J].计算机应用研究，2008，25（7）：1940-1943.

A kind of scene reconstruction optimizationmethod for building on the basis of remotely operated vehicles

GUO Yunxiang,ZHOU Jun
（School of Mechanical and Electrical Engineering,Hehai University,Changzhou 213022,Jiangsu China）

In order to overcome the cumulative error in the reconstruction process of underwater building and effectively restore the surface structure information of underwater building,a kind of scene reconstruction optimizationmethod for building on the basis of remotely operated vehicles has been put forward in the text.The maximum spanning tree has been established.The basic graph algorithm has been used to search reference coordinate system and complete the transformation of coordinate system.The error equation has been established.The optimum transformation matrix has been realized by use of L-M iterative algorithm. The reliable detection has been conducted to this method through contrast experiment.The experimental results show that the image information of the whole detection process has been effectively reflected by adopting this optimized method.It provides effective basis for follow-up machine vision processing.

Underwater building;Cumulative error;Optimized method;Max.spanning tree;Transformation matrix

TP391.7

A

10.16316/j.issn.1672-0121.2016.06.030

1672-0121（2016）06-0118-06

2016-07-16；

2016-09-20

江苏省科技支撑项目（BE2012096）

周 军（1961-），男，教授，从事微机测控自动化研究；郭云翔（1992-），男，硕士在读，主攻计算机视觉研究。E-mail：1028596390@qq.com