船载折反射全景视觉系统电子稳像算法

朱齐丹，徐从营，蔡成涛

（哈尔滨工程大学 自动化学院，哈尔滨150001）

0 引 言

折反射全景视觉系统以其360°视场观察能力，被广泛应用于车载／舰载视频监控、无人飞行器的姿态计算、机器人导航、图像测量等领域［1－2］。但全景视觉系统在船载、机载、车载等运动载体上进行拍摄时，受载体运动的影响，采集到的图像序列会产生晃动，严重影响人们观测与识别目标。

电子 稳 像（Electronic image stabilization，EIS）是指使用图像处理技术从视频图像序列中去除因摄像机随机运动而引入的图像抖动，使图像序列平滑、稳定的技术。在机载、车载、舰载电视摄像等系统得到了广泛应用。比较成熟的电子稳像算法有块匹 配算法［3－4］、相位匹配算法［5］、特征 点 匹 配 算 法［6－7］、位 平 面 匹 配 算 法［8］、投 影 算法［9］以及专门针对舰载摄像系统的水天线稳像算法［10－11］等。文献［4］提出的稳像算法采用宏块预判算法和改进的SSDA 宏块匹配算法快速计算运动矢量，结合最小二乘法与LMedS 算法求解最优的稳像模型参数，具有亚像素级的稳像精度。文献［7］基于自适应滤波和特征点集运动矢量提取技术，提出了带特征匹配验证的电子稳像算法（FMVA）。这两种方法均可以有效处理摄像机的平移和旋转抖动，实现实时稳定图像序列。文献［12］提出了一种稳像图像修补技术，判断相邻图像帧的光流变化，并利用光流法填补当前稳像图像的丢失像素，有效实现了稳像无效区域的修补。然而，以上文献研究对象多是针对普通摄像机拍摄的视频，而针对折反射全景视觉系统进行电子稳像的研究较少。这些文献虽然对于全景视频稳像有一定的参考作用，但折反射式全景摄像机与普通摄像机的成像原理存在很大差别，所以对于后者来说，仍需要深入研究适应其特点的电子稳像方法。

本文针对舰载折反射全景视觉系统拍摄的图像序列进行了电子稳像算法研究，提出了一种使用海天线进行电子稳像的方法。根据全景图像中海天线成像特点，使用基于最优边缘估计的海天线提取算法计算海天线成像方程，提出了基于海天线的简化全景图像电子稳像模型及图像无效区域重建方法。使用海上拍摄的全景序列图像进行了稳像实验，结果表明该算法可实时地进行电子稳像，有效地减小全景图像的晃动，极大地改善了视觉效果。

1 全景视觉系统成像原理

折反射全景摄像机由反射镜和普通摄像机组成，空间中的光线经过反射镜反射后进入摄像机成像，生成环形的全景图像，如图1（a）所示，两个白色环形区域是全景视觉系统的有效成像区域。当系统安装满足单视点要求时［13］，可方便地建立系统的成像模型，计算所有空间成像点对应的入射光线向量。

图1 双曲面折反射摄像机与全景图Fig.1 Hyperbolical catadioptic camera and panoramic image

图1 （b）为理想情况下拍摄的全景海洋图片，不考虑摄像机镜头畸变时可将海天线近似为圆形。当载体发生晃动时，空间景物成像将发生形变，海天线成像为椭圆，如图1（c）所示。与普通稳像图像不同，在全景图像中除载体外，背景区域由海面和天空组成。海面会不断发生变化，而天空则比较平滑，具有极少的不变特征。因此，传统的稳 像 算 法，例 如 基 于 块 匹 配［3－4］和 特 征 点 匹配［6－7］的稳像算法，在这种海面不断变化的背景模式下，会造成严重的误匹配，无法实现全局运动矢量的准确估计。同时现有的稳像算法均无法解决全景图像稳像过程中初始参考帧的选取。

与普通图像不同，全景图像发生形变导致各像素点的移动量均不一样，因此现有电子稳像算法中全局运动矢量补偿方法无法矫正全景图像中景物的形变。若建立精确的矫正模型，则必然要计算每个图像点的移动量，算法复杂，不具有实用性。因此现有的全景图像电子稳像算法均针对展开图像，使用传统稳像算法进行电子稳像，以降低算法难度［14－15］。

本文根据海天线成像特点，使用海天线成像方程进行全景图像形变矫正，提出了简化的稳像模型。该模型虽然不能完全矫正图像形变，但算法速度快，矫正效果好，可极大地改善视觉效果，并且不需要参考图像，实用性较强。

2 海天线提取

2.1 分区自适应阈值Canny边缘检测

全景图像容易受光照影响造成整幅图像对比度不均匀，传统Canny算子在进行边缘检测时使用同一对高低阈值，会导致边缘信息丢失或出现伪边缘。本文将图像均匀分为大小相同的图像块，对每个图像块使用文献［16］提出的最大类间方差法（Otsu算法）确定Canny算子的高低阈值。每个子图像阈值更能反映其所在子区域的边缘特性，检测结果既保留了局部的细节边缘，也不会产生虚假边缘。修改后Canny 边缘检测算法具体流程为：

（1）利用二维高斯函数对图像进行低通平滑滤波，去除图像中的噪声。

（2）使用2×2邻域一阶偏导的有限差分计算梯度幅值和梯度方向，记每个像素点的梯度幅值和梯度方向分别为M［i，j］和θ［i，j］。

（3）对梯度幅值进行非极大值抑制，在梯度的方向上互相比较邻接像素的梯度幅值，保留幅值局部变化最大的点。

（4）将图像均匀分为大小相同的图像块（64×64像素，也可根据图像尺寸设定），每个方块使用Otsu算法确定Canny算子的高低阈值并进行双阈值化检测和边缘连接。

（5）经过非极大值抑制获取的边缘图像还达不到单像素级边缘，为方便进行边缘长度统计，使用文献［17］提出的形态学算子对边缘检测后的图像进行边缘细化。

2.2 双阈值梯度方向过滤

与普通图像不同，全景图像中海天线边缘上像素点梯度方向指向图像中心。因此可通过判断梯度方向去除载体天线和全景装置支撑杆等造成的径向边缘并保留海天线边缘。但由于噪声、图像采样量化和载体剧烈地晃动等原因，造成部分海天线边缘点梯度方向偏离圆心，采用单一阈值会造成海天线边缘的断裂，因此本文借鉴Canny算法中双阈值化检测的思想，提出了双阈值梯度方向过滤的方法，即设置高低阈值分别为θth1、θth2，其中θth1＞θth2。对所有边缘点计算其与圆心的夹角θ′［i，j］：

其中［u0，v0］为全景图像中心。如果θ′［i，j］＜0，则θ′［i，j］＝θ′［i，j］＋360°，将角度转到0～360°范围内。若｜θ［i，j］－θ′［i，j］｜≤θth2，则保留此边缘点。这样得到的海天线边缘会因为间断而变的不连续，为了最大程度地保留海天线边缘点，采用递归边界跟踪方法，将八邻域区域内｜θ［i，j］－θ′［i，j］｜≤θth1的点判定为边缘点并保留。

图2（a）为边缘检测效果，海天线边缘全部被检测出来。图2（b）为有效区域内边缘双阈值过滤结果，取θth1＝5，θth2＝15。从图中可以看出，径向边缘被全部消除，而海天线边缘被完整地保留下来。

图2 边缘检测与过滤Fig.2 Edge detection and filtering

2.3 最优边缘估计算法

由图2（b）中可以发现边缘图像有如下几个特点：①环形海天线被分割为断续的曲线段；②图像中最长的曲线段为海天线边缘；③相邻的海天线曲线段首尾端点在切向方向上距离最近，即到图像中心点距离相近的同时保证端点距离最近。根据以上特点，本文提出了一种最优边缘估计算法提取海天线边缘，建立了最优边缘估计方程，并使用一种简单的方法判断海天线提取是否成功。

首先统计所有边缘曲线的起始点、起始角、曲线长度等信息。假设有效区域内共有m 条边缘，将所有边缘信息表示为：

式中：RAi、RBi为第i条边缘起始点和终止点到图像中心点的距离；θAi、θBi为边缘起始点和终止点的方向角；ni为边缘长度。

所有边缘信息按逆时针方向统计排列，即边缘终止点在起始点的逆时针方向。

然后使用最优边缘估计算法提取海天线边缘。根据边缘图像特点②将最长的边缘设为第1条海天线边缘，并记为Lj，将已经搜索到的海天线边缘起始／终止角，起始／终止半径分别记为：

最优边缘估计算法搜索海天线边缘具体流程为：

（1）逆时针搜索。根据边缘图像特点③，假设最优边缘估计方程为：

式中：α和β 为比例系数。由于相邻端点距离非常近，即使海天线发生严重变形，半径变化差值也会非常小，因此选取α＝0.3，β＝0.7；θth为角度偏差阈值；Rth为半径偏差阈值，端点相邻越远则半径相差阈值应该越大，因此取Rth＝1.5θth，或根据图像分辨率设置为0.03R，其中R 为理想情况下海天线成像圆半径。

选取计算值最小，即与现有边缘曲线半径相似并且距离最近的边缘作为海天线。同时将θ2和R2替换为选取的海天线边缘的θBi和RBi，并删除起始角和终止角位于θ1至θ2范围内的所有边缘曲线信息以提高下一次搜索效率。重复步骤（1）直至没有边缘满足条件公式（3）。

（2）顺时针搜索。与逆时针搜索原理相同，条件公式及偏差计算公式替换为公式（4）：

重复步骤（2）直至没有边缘满足条件公式（4）。

（3）终止条件判断。为提高椭圆拟合精度，搜索完成后海天线起始与终止角度若满足如下条件：

且已搜索到的边缘角度范围大于240°，则转到步骤（4），否则重新计算最长边缘，转到步骤（1）。

（4）椭圆拟合。根据搜索出的海天线边缘像素点位置直接使用椭圆拟合方法求解海天线椭圆成像方程。本文直接使用OPenCV 提供的椭圆拟合函数求解椭圆方程。

图3显示了海天线搜索的结果及海天线椭圆方程拟合结果，其中图像分辨率为1024×1024像素，角度偏差阈值θth＝10，半径偏差阈值Rth＝15。从图中可以看出，白线代表的椭圆方程较好地拟合出了海天线的真实位置，具有较高的精度。

图3 海天线提取Fig.3 Sea－sky line extraction

大部分情况下由于载体的晃动会使海天线成像为椭圆。多次手动提取海天线实验结果表明，假设理想情况下海天线成像圆半径为rc，海天线变形时成像的椭圆长短轴分别为a，b，则它们之间的关系满足rc≈（a＋b）／2，因此可将其作为海天线提取是否成功的判断依据，其中rc可手动确定。若提取的海天线满足如下条件公式，则认为海天线提取成功：

式中：dth为判断阈值。

3 基于海天线的全景图像电子稳像

3.1 电子稳像模型

本文所提出的简化稳像模型基于如下两个假设：①载体不发生晃动时，全景图像可由一系列相邻的密集同心圆表示，当图像发生晃动时，同心圆变为同心椭圆并发生中心点偏移；②同一圆周上所有的景物点在变形后仍会位于同一椭圆上。如图4所示，其中黑色粗实线表示海天线。实际情况下全景图像中海天线上的点对应于无穷远处景物，只有海天线上的景物点满足条件②。当载体发生晃动时，除与全景相机固定连接的载体外，所有空间景物点的投影角度都将发生变化，成像将发生形变，假设条件②并不能精确满足。

基于两个假设条件下的全景图像稳像算法即将同心椭圆矫正为同心圆，如图5（a）所示，然后将圆心Oc平移至理想情况下海天线中心点O′c，如图5（b）所示。

图4 全景图像稳像模型Fig.4 Omnidirectional image stabilization model

图5 图像坐标系设置Fig.5 Image coordinate system

图5 为稳像坐标系设置，图像坐标系原点位于图像左上角。其中XeOeYe为失稳图像同心椭圆坐标系，XcOcYc为矫正后同心椭圆坐标系，X′cO′cY′c为稳像图像同心圆坐标系，θ为海天线椭圆方程倾斜角。基于假设条件的图像稳像虽然不能完全矫正景物形变，但能极大地简化稳像算法，提高计算速度。

记理想情况下海天线成像圆心为［uc，vc］，半径为R；a、b为失稳图像中海天线成像椭圆长短轴；［ue，ve］为椭圆中心点；θ为椭圆旋转角度。假设点p ＝［i，j］为稳像图像上任意一点的图像坐标，稳像过程即计算出点p 对应失稳图像上的点坐标p′＝［i′，j′］。记r为该点到图像中心点［uc，vc］的距离，矩阵E和C 分别表示该点对应稳像前后的椭圆和圆，a′，b′为对应的椭圆长短轴，在图5（a）坐标系下椭圆和圆二次曲线方程可表示为：

式中：矩阵E 和C 表示为：

已知椭圆旋转角度为θ，则点m 在XcOcYc坐标系下可表示为：

式中：

由公式（6）知椭圆长短轴a′、b′与r 存在如下关系：

由比例关系可得：

假设矩阵K 为3×3矩阵，令：

将公式（12）带入公式（7）得：

对比公式（7）中两方程有：

对公式（14）求解可得：

计算出m′后，用下式算出最终p′点坐标：

3.2 背景模板

全景图像中部分背景不会随着载体摇晃而发生变化，这部分图像为固定背景，可进行手动提取，并设置为非稳像区域，如图6所示。

图6 背景区域Fig.6 Background region

4 实验结果

4.1 全景图像采集设备

实验时使用自主设计制造的全景图像采集处理设备，如图7所示。设备采用模块化分体设计，可方便地根据任务需要拆分使用或者换装不同分辨率的相机和反射镜。图像采集后即可直接进行本地处理，也可使用光纤模块传输至处理终端进行处理。相机具有曝光时间自动调整功能，整套装置采用三防设计，可满足长时间海洋环境使用要求，前端处理器为UNO－3083 嵌入式计算机，主要配置为intel酷睿i7 处理器，主频为2.2 GHz，内存大小为4G，操作系统为Windows XP，实验以VC6.0软件为平台结合OpenCV 库编写。

图7 实验装置Fig.7 Experimental device

4.2 算法复杂度分析和海天线提取实验

在最优边缘估计算法中，已知有效区域内共有m 条边缘，则每次搜索中计算最长边缘选取所需的最大计算次数为m，顺／逆时针搜索最优边缘估计方程最大计算次数为m－1。假设最大循环搜索次数为t，每次逆／顺时针搜索的边缘个数最大为s，则算法的最大运算时间复杂度可表示为：

在实际测试中，t和s 均为个位数，算法只需数次计算即可完成海天线边缘搜索。

为测试海天线提取算法性能，本文使用分别在港口、航道、远海3种不同实验环境下拍摄的大量图片进行了实际测试。图像分辨率为1024×1024像素，取θth1＝5，θth2＝15，θth＝10，Rth＝15。部分测试图片及海天线计算结果如图8所示。

图8中（a1）拍摄环境为远海，（a2）和（a3）拍摄环境分为航道和港口，使用单测垂直视角为120°的双曲面反射镜，较大的垂直视野使全景摄像机对载体的晃动具有更强的适应性。图中除设置有效检测区域外（两白色圆圈之间部分），根据遮挡情况设置了遮挡区域（两条白线之间部分）。

分区自适应阈值边缘检测算法有效地检测出了不同环境下的海天线区域（图8（b1）～（b3）），双阈值梯度方向过滤法则有效地过滤掉了径向边缘，并将干扰物分割为较短的边缘（图8（c1）～（c3））。图8（d1）～（d3）和图8（e1）～（e3）则分别显示了海天线提取结果及最终的椭圆拟合效果。即使在有效区域小于180°的环境下仍然精确地拟合出了海天线（图8（e3））。

图8 海天线提取实验Fig.8 Sea－sky line extraction experiment

为测试算法的成功率与实时性，本文使用连续拍摄的1700张图片进行了测试，设置不同的图像分辨率与检测参数，测试结果如表1所示。

表1 海天线提取算法性能Table1 Algorithm performance of sea－sky line extraction

部分图片由于晃动较大，海天线被舰船遮挡或超出视野，造成海天线检测失败。同时图像分辨率的下降也会影响最终椭圆拟合结果的准确性，造成海天线提取成功率下降，但由表1中可以看出影响较小。为保证海天线提取精度，进行海天线提取时图像尺寸不宜太小。

4.3 全景图像稳像实验

本文使用低分辨率256×256像素的图像进行快速海天线成像方程计算，然后对分辨率为1024×1024像素的图像进行图像稳像实验，稳像效果如图9所示。

图9 全景图像稳像Fig.9 Omnidirectional image stabilization

图9 （a）（b）为稳像矫正前后的图像对比，图9（b）中白色部分为稳像无效区域。如果不填充合适的像素值将会严重影响图像的视觉效果，因此本文提出了关键帧补偿法对稳像无效区域进行重建，改善视觉效果。

4.4 无效区域重建像实验

由于稳像对象为连续的图像序列，载体的晃动会导致周期性地出现晃动较小的图像，因此可通过提取的椭圆参数进行判断。如果椭圆长短轴非常接近，并且椭圆中心点与理想情况下海天线圆中心点非常接近，则可以认为该图像晃动较小，如图10所示。

图10 理想海天图像及稳像结果Fig.10 Ideal omnidirectional image and stabilization result

图10 为较理想的海天图像及其稳像图像。从图10（a）中可以看出图像晃动极小，因此图10（b）中的无效区域较小，使用与原图相同位置处的像素值进行填充，重建效果如图11所示。

图11 无效区域重建Fig.11 Invalid region reconstruction

由图11中可以看出，填充后的图像较为完整，视觉效果大为改善。此时可设置此稳像图像为关键帧图像，其后稳像图像中的无效区域全部填充为关键帧图像中相同位置的像素值，直至出现新的关键帧图像，关键帧图像被更新。图12示出了对1700张连续全景视频图像序列进行稳像实验的部分图像稳像效果。

图12 稳像实验结果Fig.12 Results of image stabilization

从图12中可以看出，晃动不但被抑制，而且无效区域被填充为相邻关键帧中的像素，即使图像晃动使海天线超出视野，重建后的稳像图像中仍然能保证海天区域的完整性，极大地改善了视觉效果。基于简化稳像模型的稳像算法计算简单，对分辨率为1024×1024像素的全景图像矫正和重建平均计算时间为10 ms，基于海天线的全景图像电子稳像算法可控制在50ms以内，具有较高的实时性。由于未见国内外其他全景图像整体稳像算法的相关文献报道，同时全景图像中海浪时刻变化，普通图像电子稳像算法中的帧差法和峰值信噪比等评价方法均不适用，因此未对稳像效果作进一步评价。

5 结束语

本文针对船载折反射全景视觉系统采集的图像序列进行了研究，并提出了基于海天线的电子稳像算法。首先使用最优边缘估计算法提取海天线边缘，进行椭圆拟合得到海天线成像椭圆方程；然后建立稳像模型并进行图像稳像；最后根据全景海天图像特点对稳像无效区域进行重建。实验结果证明，基于海天线的电子稳像算法可有效抑制海天线晃动，并具有计算速度快、成功率高的优点，在实际应用中取得了良好的使用效果。

［1］Mondragod I F，Pascual C，Caro M，et al.Omnidirectional vision applied to unmanned aerial vehicles（UAVs）attitude and heading estimation［J］.Robotics and Autonomous Systems，2010，58（6）：809－819.

［2］Mariottini G L，Scheggi S，Morbidi F，et al.An accurate and robust visual－compass algorithm for robot－mounted omnidirectional cameras［J］.Robotics and Autonomous Systems，2012，60（9）：1179－1190.

［3］蔡成涛，梁小龙，谭吉来，等.自适应最优块匹配视频稳像算法［J］.系统工程与电子技术，2013，35（6）：1324－1329.Cai Cheng－tao，Liang Xiao－long，Tan Ji－lai，et al.Adaptive optimal block matching video stabilization algorithm［J］.Systems Engineering and Electronics，2013，35（6）：1324－1329.

［4］仲训昱，朱齐丹，张 智.数字稳像中的快速和鲁棒运动估计研究［J］.电子学报，2010，38（1）：251－256.Zhong Xun－yu，Zhu Qi－dan，Zhang Zhi.Study of fast and robust motion estimation in the digital image stabilization［J］.Acta Electronica Sinica，2010，38（1）：251－256.

［5］邵向鑫，郭树旭，王朗.基于边缘扩展相位相关的图像拼接算法［J］.吉林大学学报：信息科学版，2010，28（1）：95－99.Shao Xiang－xin，Guo Shu－xu，Wang Lang.Image mosaic algorithm based on extended phase correlation of edge［J］.Journal of Jilin University：Information Science Edition，2010，28（1）：95－99.

［6］王洪，戴明.区域快速鲁棒特征跟踪电子稳像［J］.吉林大学学报：工学版，2012，42（2）：451－458.Wang Hong，Dai Ming.Area speed up robust feature tracking for electric image stablization［J］.Journal of Jilin University：Information Science Edition，2012，42（2）：451－458.

［7］朱娟娟，郭宝龙.电子稳像的特征点跟踪算法［J］.光学学报，2006，26（4）：516－521.Zhu Juan－juan，Guo Bao－long.Features tracking algorithm for video stabilization［J］.Acta Optica Sinica，2006，26（4）：516－521.

［8］Zhao H，Wang T.Image stabilization algorithm based on multi－bitplane［C］∥Proc of the International Symposium on Photo－electronic Detection and Imaging，Beijing，2009：73842Z.

［9］Tian S，Wang N，Zhao P，et al.Aims at moving objects improvement based on gray projection of algorithm of the electronic image stabilization［C］∥Proc of the 3th International Congress on Image and Signal Processing，Yantai，2010：2483－2487.

［10］赵红颖，晏磊，熊经武.船载摄像系统的一种电子稳像算法［J］.光学技术，2003，29（5）：628－633.Zhao Hong－ying，Yan Lei，Xiong Jing－wu.A new electronic image stabilizing algorithm for the image sequence of camera system on board［J］.Optical Technique，2003，29（5）：628－633.

［11］Xue Y T，Wu T.A new horizon electronic image stabilization algorithm based on SVM［C］∥2010International Conference on Optoelectronics and Image Processing，Haiko 2010：59－63.

［12］Matsushita Y，Ofek E，Ge W N，et al.Full－frame video stabilization with motion inpainting［J］.IEEE Trans Pattern Analysis and Machine Intelligence，2006，7（28）：1150－1163.

［13］马子岭，王建中.折反射全景侦察系统电子稳像技术研究［J］.北京理工大学学报，2010，30（12）：1475－1479.Ma Zi－ling，Wang Jian－zhong.Research on electronic image stabilizing algorithm of catadioptric panoramic reconnaissance system［J］.Transactions of Beijing Institute of Technology，2010，30（12）：1475－1479.

［14］朱齐丹，马宏业，左小祥.全景摄像系统中展开图像的电子稳像方法研究［J］.计算机应用研究，2009，26（3）：1192－1194.Zhu Qi－dan，Ma Hong－ye，Zuo Xiao－xiang.Electronic image stabilizing method research for developing image in omnidirectional camera system［J］.Application Research of Computers，2009，26（3）：1192－1194.

［15］Baker S，Nayar S.A theory of catadioptric image formation［C］∥Sixth International Conference on Computer Vision，Bombay，1998：35－42.

［16］Mei F，Yue G X，Yu Q C.The study on an application of otsu method in canny operator［C］∥Proceedings of the 2009International Symposium on Information Processing，Huangshan，2009：109－112.

［17］谭郁松，周兴铭.保持连通的边缘细化算法［J］.国防科技大学学报，2004，26（4）：51－56.Tan Yu－song，Zhou Xing－ming.Connectivity preserved edge thinning algorithm［J］.Journal of National University of Defense Technology，2004，26（4）：51－56.