随着经济快速发展和现代化交通工具的不断普及, 机场、公路交通以及航道对能见度的依赖日趋突出^[1]. 尤其近年, 我国各地区雾霾天气不断加剧, 持续低能见度天气导致的恶性交通事故时有发生. 可见, 能见度的实时观测与预报日趋成为不可或缺的重要内容. 但能见度的观测现状却与实际需求相差甚远.

目前, 国内外对能见度的探测虽有成熟产品, 但价格和性能之间很难调和^[2-5]: 观测效果较佳的仪器价格过高. 反之, 价格较低的仪器观测效果又不理想.

而基于CCD数字摄像技术, 在我国始于上世纪末研制的数字摄像能见度仪(Digital Photography Visio-meter System, DPVS)^[6-11], 不但观测效果较为理想, 而且物美价廉. DPVS仿照人眼观测能见度原理, 完全遵照能见度定义研制, 在国内外尚属首创, 是能见度观测的理想仪器. 但也因为是创新产品, 因此存在一些研发障碍.

本文正是基于图像配准技术, 解决DPVS中观测目标物自动搜索和准确定标问题. 这是视频能见度观测中急待解决的首要科学难题, 也是多年来视频能见度观测研究中急待解决的主要科学难题. 因为定标略有偏差, 观测准确率就会受到严重影响. 这是因为能见度的观测值大小取决于目标物和天空背景之间的亮度差异(或亮度对比度)^[12-15]. 当目标物的亮度出现偏差的时候, 将会直接影响观测的准确性. 因此, 只有实现目标物自动和准确定标, 才能实时定位目标物, 始终保持定标准确. 这是确保数字能见度仪观测准确率的关键.

模板匹配是匹配算法当中的经典算法^[16,17]. 但是这种匹配技术存在匹配耗时长, 且无法抵抗旋转的问题. 但DPVS在观测时, 其摄像头在地面震动或者大风吹动的情况下, 会存在一定小角度的旋转以及抖动等情况. 因此, 针对实际情况, 本文将模板匹配算法改动, 提出针对小角度旋转匹配不变的一种快速算法, 该方法使匹配时间呈指数级数减少, 极大地提高了计算速度, 并能抵抗一定角度的旋转.

1 卷积与模板相关运算

卷积运算简单快速, 是图像处理中的常用算法. 本文利用卷积运算与模板运算之间的相似性, 将模板相关算法和卷积进行有效融合, 得到快速的目标自动搜索与定位算法.

1.1 卷积运算

针对二维图像, 设 $x$ 与 $y$ 为两独立变量, 其卷积公式如下所示^[18]:

$ \begin{split} h(x,y)& = f*g \\ &= \int_{ - \infty }^\infty {\int_{ - \infty }^\infty {f(u,v)g(x - u,y - v){\rm d}u{\rm d}v} } \end{split}$

(1)

而对于图像卷积公式, 因为其自变量取整数, 故双重积分, 可改为双重求和. 则有:

$\left\{\begin{split} &H = F*G\\ &H(i,j) = \sum\limits_m {\sum\limits_n {F(m,n)G(i - m,j - n)} } \end{split}\right.$

(2)

式中,F是图像, G是卷积核, H是图像进行卷积后的结果.

1.2 互相关算法

进行图像匹配定位时, 设大小为 ${M_r} \times {N_r}$ 的基准图为 ${{{G}}_r}$ , 大小为 ${M_s} \times {N_s}$ 的模板为 ${{{G}}_s}$ , 并且基准图大小要大于模板, 即 ${M_s} < {M_r}$ , ${N_s} < {N_r}$ . 则模板与其覆盖下的对应子图 ${{{G}}_r}\left( {u,v} \right)$ 之间的去均值归一化互相关度量 $\rho \left( {u,v} \right)$ ( $\left( {u,v} \right)$ 为基准图左上角), 可用下式计算^[19]:

$ \rho \left( {u,v} \right) = \frac{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {[{{{G}}_r}\left( {i + u,j + v} \right) - \overline {{{{G}}_r}\left( {u,v} \right)} ]} } \times [{{{G}}_s}\left( {i,j} \right) - \overline {{{{G}}_s}} ]}}{{\sqrt {{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {[{{{G}}_r}\left( {i + u,j + v} \right) - \overline {{{{G}}_r}\left( {u,v} \right)} ]} } }^2}} \sqrt {{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {[{{{G}}_s}\left( {i,j} \right) - \overline {{{{G}}_s}} ]} } }^2}} }}$

(3)

其中, $ - 1 \le \rho \left( {u,v} \right) \le 1$ , $\overline {{{{G}}_s}} $ 和 $\overline {{{{G}}_r}\left( {u,v} \right)} $ 分别为 ${{{G}}_s}$ 和 ${{{G}}_r} $ $\left( {u,v} \right)$ 的灰度均值. $\rho $ 最大之时, 正是模板与基准图中所对应子图最相关之处, 即为所求.

2 互相关算法的改进

光受卷积启发, 针对互相关算法公式进行简化, 使之具备与卷积相同的运算模式, 极大降低算法的计算量.

在式(3)中, 将其分母的两项分别设为:

$\left\{ \begin{array}{l} {{{G}}_{{S_{i,j}}}} = \sqrt {{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {[{{{G}}_s}\left( {i,j} \right) - \overline {{{{G}}_s}} ]} } }^2}} \\ {{{G}}_{{R_{i,j}}}} = \sqrt {{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {[{{{G}}_r}\left( {i + u,j + v} \right) - \overline {{{{G}}_r}\left( {u,v} \right)} ]} } }^2}} \\ \end{array} \right.$

代入(3)式中, 则有:

$\rho \left( {u,v} \right) = \frac{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {[{{{G}}_r}\left( {i + u,j + v} \right) - \overline {{{{G}}_r}\left( {u,v} \right)} ]} } \times [{{{G}}_s}\left( {i,j} \right) - \overline {{G_s}} ]}}{{{{{G}}_{{R_{i,j}}}}{{{G}}_{{S_{i,j}}}}}}$

(4)

为了降低计算复杂度, 考虑以模板为基准, 即其他图像与模板的相似程度, 而非与模板的互相似程度. 因此式(4)可用模板的各像素平方之和来进行归一化. 同时, 假设在此计算的是二值图像, 二值化图像通常仅存0与255或0与1, 则可以消去灰度均值的计算. 因此, 式(4)可简化如下:

${\rho '}\left( {u,v} \right) = \frac{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {{{{G}}_r}\left( {i + u,j + v} \right)} } \times {{{G}}_s}\left( {i,j} \right)}}{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {{{{G}}_s}{{\left( {i,j} \right)}^2}} } }}$

(5)

显然, 式(5)中的分母是固定常数, 是模板各像素平方的和, 设为K, 则公式可进一步简化为下式:

${\rho '}\left( {u,v} \right) = \frac{1}{K}\left( {\sum\limits_{i = 1}^{{M_s}} {\sum\limits_{j = 1}^{{N_s}} {{{{G}}_r}\left( {i + u,j + v} \right)} } \times {{{G}}_s}\left( {i,j} \right)} \right)$

(6)

其中, $0 \le {\rho '} \le 1$ , 此时, 上式已经和卷积计算式(2)高度相似.

因为是二值图像, 若图像中像素是0与1, 则上式计算本质上, 和两图像对应像素“与运算”之和的结果等效. 否则就与两图像对应像素“与运算”结果平方之和的结果等效. 用数学表达式, 可描述如下:

$\left\{ \begin{array}{l} {\rho '}\left( {u,v} \right) = \dfrac{1}{K}\left( {\displaystyle\sum\limits_{i = 1}^{{M_s}} {\displaystyle\sum\limits_{j = 1}^{{N_s}} {{{{G}}_r}\left( {i + u,j + v} \right) \cdot {{{G}}_s}\left( {i,j} \right)} } } \right)\\ {{{G}}_r}\left( {i + u,j + v} \right) \cdot {{{G}}_s}\left( {i,j} \right) = \left\{ \begin{array}{l} 0{\rm{ }}\\ 1{\rm{ }} \end{array} \right.{\rm{ }}{\text{当像素为}}0{\text{和}}1{\text{时}}\\ {{{G}}_r}\left( {i + u,j + v} \right) \cdot {{{G}}_s}\left( {i,j} \right) = \left\{ \begin{array}{l} 0{\rm{ }}\\ {n^2} \end{array} \right.{\rm{ }}{\text{当像素为}}0{\text{和}}n{\text{时}} \end{array} \right.$

(7)

式中, “ $ \cdot $ ”代表“与运算”. 因为 ${G_r}$ 与 ${G_s}$ 仅在有限范围内非零, 所以计算仅需在非零部分重叠的区域进行即可. 这大大降低了计算复杂度, 无论对算法硬件的实现还是对计算的快捷性均具有极大的价值.

需要指出的是, ${\rho '}$ 的计算仅仅基于模板为基准. 因此当基准图中模板覆盖对应下的子图 ${{{G}}_r}\left( {u,v} \right)$ 像素值全部都是 $n$ 时, 将会导致 ${\rho '} = 1$ 的错误匹配. 针对此种情况, 设模板对应的子图像素值之和为 $\Sigma $ , 令 $\Delta < \Sigma \le { T}$ , T是模板像素的和, $\Delta $ 可以根据需要调节, 本文取 $\Delta = { T}/2$ . 这样, 实际匹配的过程中, 仅计算模板对应的子图像素值之和在此范围之内的即可, 从而进一步加速算法的计算速度.

3 目标自动捕获和识别流程

以目标黑体为例(目标光源自动定位与黑体类似, 过程不再赘述), 其流程如图1所示.

4 实验分析

以为了验证算法的有效性, 将算法整合在DPVS系统中, 透过在DPVS中使用和去除本文算法来进行能见度观测结果的对比, 来进行分析本文算法的有效性.

实验观测场所选在北京市观象台. 观象台是国家基准气候站, 功能完备, 拥有专业的人工观测员与各类观测仪器和高分辨率的能见度人工观测数据, 以及各类丰富完整的气象要素观测数据, 是开展能见度仪观测分析的理想站点.

现将利用本文算法定标前后的效果图, 显示如图2~图5.

相机抖动会造成背景图发生相对定标框的偏移或旋转. 在实际中, 定位的目标是保证定标框能自动定位到正确位置. 因相机抖动造成的旋转和移动不宜量化. 为验证算法的有效性, 一方面可透过视频窗直接观测定位效果进行判断, 如图5所示. 另一方面, 可以利用单帧背景图, 以1度为单位, 进行旋转仿真实验, 定标框只要落在光源区即认为定位正确, 即为1. 否则, 定位失败, 即为0, 则有如图6关系图.

可以看到, 当角度旋转在一个较小的角度内时(大约[–10, 10]之间), 本文算法能够保持比较稳定的定位准确性. 而上图是显得较小的第二光源所得角度与定位间的关系图, 第一光源因为面积较大, 故抗旋转性略强, 角度范围大致在[–15, 15]之间. 黑体与光源类似, 不再赘述.

系统于2017年7月12~13日和22~23日, 将本文快速定位算法在能见度观测系统上进行了功能去除. 获取的实验观测系统显示图, 如图7所示.

现将加入本文算法, 现阶段系统的观测效果图显示如图8~图10所示.

图8~图10中, DPVS_77为DPVS系统观测, DNQ2-1为国产的(无锡)前向散射仪观测, PWD22为进口前向散射仪观测, LT31为透射仪观测.

为了进一步检验本算法的有效性, 考虑到手工目标定位的难度, 考虑从2018年8月的数据中, 按观测距离进行随机抽样. 为了便于分析, 将能见度(简称V, 以下同)观测数据从低到高分为4组, 即: 小于2 km、2–5 km、5–10 km、10–15 km. 首先以DPVS、LT31、PWD22及DNQ2-1等4种观测数据的平均值作为标准值, 分别计算每个数据段DPVS和LT31与标准值之间的相对误差和相对标准差. 相对标准差计算公式如式(8)所示:

${\sigma _{rd}} = \sqrt {\dfrac{{\displaystyle\sum\limits_{i = 1}^N {{{[({X_i} - {Y_i})/{Y_i}]}^2}} }}{{N{\rm{ - }}1}}} \times 100\% $

(8)

其中, X为DPVS或LT31的测量值, Y为DPVS、LT31、PWD22以及DNQ2-1这4种观测数据的平均值, N为样本个数. 具体分析结果详见表1. 基准图大小皆为640×480. 普通算法, 即互相关算法(以下同). 手动定位精度和手动定位的频率密切相关, 考虑到实际可操作性, 这里的频次为1小时重新定位一次. 因为手动定位没有办法统计定位时间, 也无实际意义. 因此, 这里仅将本文算法单次定位耗时与普通算法耗时对比如表2所示.

本文算法的优越性显而易见, 观测准确度比普通算法提升7倍多, 比手动定位观测提升将近5倍, 而匹配速度比普通算法平均提升200多倍. 而本文算法与LT31和标准值的相对标准差都在20%之内, 均符合世界气象组织对能见度仪研制标准的要求.

普通算法目标定位之所以引起的观测误差最大, 是因为算法虽然能够进行自动定位, 但匹配速度不够, 无法在有效观测时间内(单次观测还需要亮度调节、视频图像抽取和处理等其他耗时过程), 及时自动定位, 导致观测失败. 同时, 原有算法的抗旋转能力较弱, 会陷入局部最优, 导致错误定位, 进而加剧观测的误差.

手动定位, 在镜头相对比较稳定的情况下(无大风或者车辆等外在原因引发的镜头震动), 人工定位的代价还可以接受, 且观测准确度还是比较高的, 但是这种无震动的条件不符合现实, 在实际中往往需要不断的进行定位以对抗震动引发的定位偏离. 这种人工操作, 代价太大.

综合对比分析数据, 可知本文算法效果最佳. 在进行自动定位时, 本文算法单次定位耗时平均在24.2 ms左右, 远远小于最小的单次观测时间间隔: 1 min. 而且本文算法定位精度高, 一方面是因为增加了定位策略, 使搜索集中在最大可能区域, 另一方面因为修改了互相关的算法, 大大消弱了定位时的相关性, 从而提升了对镜头晃动所引发的图片微小旋转的抵抗力.

5 结论与展望

通本文针对数字摄像能见度仪因相机震动造成的目标物定位偏离问题, 提出了一种快速自动定位解决技术, 确保准确获取数字摄像能见度仪实际观测时的目标物, 即黑体和光源的亮度, 为准确计算大气能见度提供了有效支撑. 大量实验表明该方法的有效性. 作为一种新型的、唯一符合能见度定义的大气能见度观测仪, DPVS具有多种优点: 仿照人眼观测能见度原理, 完全遵照能见度定义研制设计, 其安装便捷、成本低廉, 是最适合取代人工观测能见度的仪器. 而随着CCD-CAMERA技术的迅猛发展, 其成本不断下降, 而分辨率却不断提升, 早已达到与人眼相当的千万像素级. 因此, DPVS具有广阔良好的应用前景.