航空图像光流场的逆向金字塔计算方法

李佳田，李显凯，李应芸，钱堂慧，李果家，林　艳

1. 昆明理工大学国土资源工程学院，云南 昆明 650093； 2. 中国人民公安大学警务信息工程学院，北京 100038



航空图像光流场的逆向金字塔计算方法

李佳田1，李显凯1，李应芸1，钱堂慧1，李果家1，林艳2

1. 昆明理工大学国土资源工程学院，云南 昆明 650093； 2. 中国人民公安大学警务信息工程学院，北京 100038

航空图像光流场是低空运动目标检测与变化信息获取的基础，通常将图像金字塔结构引入数值过程以增强全局收敛性。然而，金字塔结构往往是由底层至顶层的递进方式构建，其忽略图像的几何成像过程，造成微小光流或不能得到光流的问题，导致难以支撑后续建模与分析。本文提出了一种以顶层图像为基准的逆向金字塔结构，首先依据中心投影定量地计算出顶层图像的降采样因子，使得顶层图像光流能够反映所设定的地面目标位移阈值；其次，结合顶层与原始图像，以等比方式确定中间层降采样因子；最后，利用高斯平滑与图像插值得到中间层图像，并形成金字塔。对比试验与分析表明，逆向金字塔可准确地计算航空图像光流场，在抑制地面微小位移方面具有优势。

航空图像；光流场；逆向金字塔；中心投影；降采样

光流场(optical flow field)是空间运动目标速度场在二维图像平面的投影[1]。由于光流场能够反映出图像的变化，为观察者提供目标的运动信息和几何结构信息，所以其在摄影测量与遥感领域得到越来越多的关注。文献[2—3]利用光流场进行图像变化检测分析；文献[4—5]提出用光流向量与图像特征相结合的低空视频运动车辆检测方法；文献[6]利用光流的运动信息描述地物变化，设计高分辨率遥感影像土地利用与土地覆盖变化的信息自动获取方法；文献[7]则根据光流运动向量能够反映出影像间不同物体的变化，实现车载影像序列的全景匹配。

为寻求全局最优解及提高算法效率，目前多是将金字塔结构嵌入光流的数值计算过程，其核心在于通过降采样构造出图像的多个比例尺空间，然后由小比例尺到大比例尺以传递的方式计算光流场[8-11]。当光流场算法用于近景场景时，可以回避图像的几何成像过程，通常是以原始图像为基础的等比降采样方式来构建图像金字塔。然而，航空图像是在距离地物较远处获得，直接采用等比图像金字塔缺少降采样的理论依据，即地面目标位移量与光流向量之间的定量关系，进而会导致地面目标移动不能有效地在光流场中反映的问题。以中心投影成像模型为基础，本文通过建立成像大小与航高之间的关系来描述降采样比例因子，提出适用于航空图像的逆向金字塔结构，使得光流场能够定量地反映地面目标运动变化。

1　逆向金字塔

1.1降采样因子的中心投影描述

航空相片是地面经过中心投影后在承影面的构象，即通过曝光瞬间将收集到的地物反射光线直接在感光胶片上感光，得到与地面地物亮度一致的正像相片的过程。如图1(a)所示，f为焦距，H为航高，S是投影中心，地面上的线段aobo经过中心投影后的构象是线段aibi。根据中心投影模型，则物方、像方具有如下的成像关系

aibi/aobo=f/H

(1)

aobo>H/f

(2)

式(2)说明在焦距为f、航高为H的情况之下，理论上像方所能分辨的物方最小位移大小为H/f。更进一步，由式(2)可知，地面目标的位移量在航向或侧向方向的投影值大于H/f时，其在图像上才能有大于1个像素的反映，即可得到光流向量。

通常，将光流应用于某类地物目标，因此，物方位移植在一定范围之内，需要确定光流向量所能区分的位移下限。与此同时，要求所对应的像方位移不能过小。当像方位移过小时，一是光流计算模型会快速收敛，并认为没有变化发生，二是环境光变化会与像方微小位移相互混淆，得到的光流向量不便于观察与分析。为此，在原始图像中由物方单位位移do与式(1)可得所对应的像方单位位移值di如下

di=fdo/H

(3)

设物方最小可分辨位移为λ，则有物方位移比为

λ/do=λf/diH=c0

(4)

αc0=αλf/(diH)=λf/(diHα-1)

(5)

如图1(a)所示，式(5)的物理意义是，相对于原始图像，当图像的物方位移比为αc0时，相当于焦距减小为αf，或者航高H增大为H/α，即图像成像的分辨率大小为原始图像的α2倍，如果置原始图像为金字塔底层，则顶层图像的边长应降采样为原始图像边长的α倍。

需要金字塔中各层图像均能够识别出物方最小可分辨位移λ，即可用式(6)来确定α

αλ/do>1或αλf/(diH)>1

(6)

即α>1/c0。α越大，降采样图像越接近于原始图像；α越小，则降采样图像的分辨率越小即尺寸越小。取α的下界1/c0以建立金字塔顶层图像，则其最小降采样分辨率因子γmin为

γmin=1/c0+ε

(7)

式中，ε为小的正数，即可用γmin得到图像金字塔中具有最小分辨率的顶层图像，并且地面上大于λ的位移所构成的光流可由顶层图像计算得出。

图1　结合中心投影描述的逆向金字塔构建Fig.1　Construction of backward pyramid using central projection

1.2逆向金字塔构建

根据γmin得到金字塔顶层图像大小后，即可依据等比方式计算金字塔中间层图像的降采样因子γ，以w表示原始图像的宽度，则有如下关系

wγn=wγmin

(8)

式中，n为中间层图像数量，如果置原始相片为金字塔的底层，则n+1表示金字塔的层数。以3层金字塔为例，则n为2，即有

(9)

确定采样因子后即可建立金字塔，金字塔构建主要包括高斯滤波和降采样。首先采用高斯滤波作用于原始图像以消除噪声，降低异质点影响，然后通过降采样为金字塔的每一层获取图像数据。

高斯滤波属于线性平滑滤波，常被用于图像去噪，其实质是对图像进行加权平均的过程[14-15]。定义滤波器gfilter，根据正态分布与3σ原则可知，位于当前像素点3σ个像素外的像素点不会对当前像素点产生影响[16]，则设置滤波器宽度k=3σ×2+1

(10)

式中，i表示滤波器中像素的位置；di表示其对应像素点与当前像素点的距离；为避免图像过度平滑，σ通常设置为1。高斯滤波有两种作用方式，一是采用滤波模板与图像作卷积，另一种是通过傅里叶变换实现。σ为1时，模板大小为7×7，根据高斯函数的可分离性，可采用两个方向的分离滤波器进行叠加，即把二维窗口卷积分离为两次一维卷积运算。

在降采样过程中，亚像素(subpixel)值的求取可通过图像插值进行，常用的图像插值方式有双线性插值[17]和双三次插值[18]。以p表示当前像素点的像素值，p1、p2、p3、p4表示当前像素点的四邻域像素点(方位顺序为左上、左下、右上、右下)的像素值，x、y表示当前像素点的坐标，[·]表示取整数，则有

(11)

双线性插值计算公式如下p(dx，dy)=(1-dx)(1-dy)p1+(1-dx)dyp2+

(1-dy)dxp3+dxdyp4

(12)

双三次插值采样需要提供其相邻的16个像素点的像素值，其获得的图像具有相对更好的平滑性，其插值公式为

(13)

式中，aij是插值系数，可以通过对式(13)求偏导数之后代入相邻像素点的像素值解出。

逆向高斯金字塔构建过程为：首先依据式(4)与式(7)确定金字塔顶层图像的降采样因子γmin，并根据γmin完成金字塔顶层图像的采样；然后用式(8)由γmin计算中间层等比降采样因子γ；当γ的值确定后以原始图像为基础，分别以γ1、γ2、…、γn-1为降采样因子完成中间层图像的采样。在采样之前需要先对原始图像进行高斯滤波，然后通过图像插值即可获取降采样图像中每一像素点的像素值。

2　试验与分析

2.1光流场计算过程

航空图像光流场的实质是计算图像间同名像素点的偏移向量。首先建立两帧图像的逆向金字塔，形式化为式(14)

BP={Ii，i=0，1，…，n}

(14)

式中，Ii为逆向金字塔中的第i层，I0为底层，即原始图像，而In为顶层。

将逆向金字塔嵌入HS(Horn-Schunck)光流场算法[8，18]，描述如下：

输入：图像I1、图像I2，物方最小可分辨位移值λ。

输出：光流向量。

(1) 根据λ分别构建I1与I2的逆向金字塔BP1与BP2。

(2) 假设各像素点均未发生位移，设顶层图像中各像素点的光流向量初始值为0，即水平、垂直分量为0，u=0，v=0。

(3) 计算金字塔第i层的光流向量，每个像素点光流分量的求解通过雅克比迭代法得到，公式如下

(15)

式中，k表示迭代次数；β为平滑系数。

(d) 由式(15)计算每个像素点的光流向量。

(e) 判断是否满足迭代收敛条件，如满足转至步骤(4)，否则，重复(c)—(e)。

(4) 光流的由顶层至底层传递计算。

(a) 将步骤(3)中解出的ui、vi按金字塔采样因子进行层间扩大，计算如下式

(16)

(b)i=i-1，判断i层图像和原始图像的关系，当i>0时，以u、v作为第i层光流初值，重复步骤(3)—(4)；当i=0时，以u、v作为第i层光流初值，重复步骤(3)。

(5) 输出像素点光流向量。

(6) 算法结束。

2.2逆向金字塔作用分析

试验地点为云南省曲靖市罗平县，试验区域面积为2 km2。航拍相机为PHASE-ONE IXU-R 1000，镜头焦距为50 mm，相对航高1087 m，图像大小为11 608×8708像素，地面分辨率为0.1 m/像素。采用连续曝光方式获取多张图像，曝光时间间隔为2 s。经过内业处理后得到正摄图像，选取前后时刻相同的图像区域进行光流场计算。分别采用逆向金字塔与正向金字塔，对地面目标大位移、地面目标小位移两种情况的光流场计算进行对比分析如下。

2.2.1地面目标大位移对比试验

图2　地面目标大位移光流计算对比Fig.2　Comparison of optical flow computation for ground target with large displacement

2.2.2地面目标小位移对比试验

图3　地面目标小位移光流计算对比Fig.3　Comparison of optical flow computation for ground target with small displacement

需要说明的是，分别用双线性插值与双三次插值进行图像降采样，在本文试验中并未表现出明显差异，从计算效率上考虑，采用双线性插值为宜。

3　结　论

用航空图像计算光流场时，首先需要解决的问题是降采样因子的确定，然而，目前仍然缺少定量依据。本文的创新之处是，以中心投影几何成像过程为基础定量地得到降采样因子，提出了一种逆向金字塔结构，对比试验结果表明，其在精确性与鲁棒性方面具有优势。此外，逆向金字塔建立过程简单，已知条件为航空图像地面分辨率。

[1]宋爽, 杨健, 王涌天. 全局光流场估计技术及展望[J]. 计算机辅助设计与图形学学报, 2014, 26(5): 841-850. SONG Shuang, YANG Jian, WANG Yongtian. Technology and Prospect of Global Optical Flow[J]. Journal of Computer-aided Design & Computer Graphics, 2014, 26(5): 841-850.

[2]李秀智, 贾松敏, 尹晓琳, 等. 视觉光流矢量场估计算法综述[J]. 北京工业大学学报, 2013, 39(11): 1638-1643.

LI Xiuzhi, JIA Songmin, YIN Xiaolin, et al. Review on Optical Flow Vector Fields Estimation Algorithm[J]. Journal of Beijing University of Technology, 2013, 39(11): 1638-1643.

[3]涂志刚, 谢伟, 熊淑芬, 等. 一种高精度的TV-L1光流算法[J]. 武汉大学学报(信息科学版), 2012, 37(4): 496-499.

TU Zhigang, XIE Wei, XIONG Shufen, et al. An Efficient TV-L1Optical Flow Method[J]. Geomatics and Information Science of Wuhan University, 2012, 37(4): 496-499.

[4]刘慧, 李清泉, 高春仙, 等. 利用C_SURF配准的空基视频运动目标检测[J]. 武汉大学学报(信息科学版), 2014, 39(8): 951-955.

LIU Hui, LI Qingquan, GAO Chunxian, et al. Moving Target Detection Using C_SURF Registration[J]. Geomatics and Information Science of Wuhan University, 2014, 39(8): 951-955.

[5]刘慧, 李清泉, 曾喆, 等. 利用低空视频检测道路车辆[J]. 武汉大学学报(信息科学版), 2011, 36(3): 316-320. LIU Hui, LI Qingquan, ZENG Zhe, et al. Vehicle Detection in Low-altitude Aircraft Video[J]. Geomatics and Information Science of Wuhan University, 2011, 36(3): 316-320.

[6]闫利, 巩翼龙, 张毅, 等. 光流动态纹理在土地利用/覆盖变化检测研究中的应用[J]. 光谱学与光谱分析, 2014, 34(11): 3056-3061.

YAN Li, GONG Yilong, ZHANG Yi, et al. Application of Optical Flow Dynamic Texture in Land Use/Cover Change Detection[J]. Spectroscopy and Spectral Analysis, 2014, 34(11): 3056-3061.

[7]张正鹏, 江万寿, 张靖. 光流特征聚类的车载全景序列影像匹配方法[J]. 测绘学报, 2014, 43(12): 1266-1273. DOI: 10.13485/j.cnki.11-2089.2014.0172.

ZHANG Zhengpeng, JIANG Wanshou, ZHANG Jing. An Image Match Method Based on Optical Flow Feature Clustering for Vehicle-borne Panoramic Image Sequence[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(12): 1266-1273. DOI: 10.13485/j.cnki.11-2089.2014.0172.

[8]ALVAREZ L, WEICKERT J, SNCHEZ J. Reliable Estimation of Dense Optical Flow Fields with Large Displacements[J]. International Journal of Computer Vision, 2000, 39(1): 41-56.

[9]BRUHN A, WEICKERT J, SCHNÖRR C. Lucas/Kanade Meets Horn/Schunck: Combining Local and Global Optic Flow Methods[J]. International Journal of Computer Vision, 2005, 61(3): 211-231.

[10]李秀智, 尹晓琳, 贾松敏, 等. 改进的TV-L1平滑光流估计[J]. 光学学报, 2013, 33(10): 1015002.

LI Xiuzhi, YIN Xiaolin, JIA Songmin, et al. Improved TV-L1Algorithm for Smooth Optical Flow[J]. Acta Optica Sinica, 2013, 33(10): 1015002.

[11]陈震, 张聪炫, 晏文敬, 等. 基于图像局部结构的区域匹配变分光流算法[J]. 电子学报, 2015, 43(11): 2200-2209.

CHEN Zhen, ZHANG Congxuan, YAN Wenjing, et al. Region Matching Variational Optical Flow Algorithm Based on Image Local Structure[J]. Acta Electronica Sinica, 2015, 43(11): 2200-2209.

[12]张聪炫, 陈震, 黎明. 金字塔光流三维运动估计与深度重建直接方法[J]. 仪器仪表学报, 2015, 36(5): 1093-1105. ZHANG Congxuan, CHEN Zhen, LI Ming. Direct Method for 3D Motion Estimation and Depth Reconstruction Based on Pyramid Optical Flow[J]. Chinese Journal of Scientific Instrument, 2015, 36(5): 1093-1105.

[13]胡正华, 孟令奎, 张文. 面向关系数据库扩展的自适应影像金字塔模型[J]. 测绘学报, 2015, 44(6): 678-685. DOI: 10.11947/j.AGCS.2015.20140279.

HU Zhenghua, MENG Lingkui, ZHANG Wen. Relational Database Extension Oriented, Self-adaptive Imagery Pyramid Model[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(6): 678-685. DOI: 10.11947/j.AGCS.2015.20140279.

[14]谢剑斌, 王晖, 程江华, 等. 基于HS约束与轮廓条件的光流场计算[J]. 系统工程与电子技术, 2009, 31(4): 761-763. XIE Jianbin, WANG Hui, CHENG Jianghua, et al. Optical Flow Computation Method Based on HS Constraint and Outline Condition[J]. Systems Engineering and Electronics, 2009, 31(4): 761-763.

[15]高银, 云利军, 石俊生, 等. 基于各向异性高斯滤波的暗原色理论雾天彩色图像增强算法[J]. 计算机辅助设计与图形学学报, 2015, 27(9): 1701-1706.GAO Yin, YUN Lijun, SHI Junsheng, et al. Enhancement Dark Channel Algorithm of Color Fog Image Based on the Anisotropic Gaussian Filtering[J]. Journal of Computer-aided Design & Computer Graphics, 2015, 27(9): 1701-1706.

[16]肖进胜, 杜康华, 涂超平, 等. 基于多聚焦图像深度信息提取的背景虚化显示[J]. 自动化学报, 2015, 41(2): 304-311.

XIAO Jinsheng, DU Kanghua, TU Chaoping, et al. Bokeh Display Based on Depth Information Extraction of Multi-focus Images[J]. Acta Automatica Sinica, 2015, 41(2): 304-311.

[17]王昊京, 王建立, 王鸣浩, 等. 采用双线性插值收缩的图像修复方法[J]. 光学 精密工程, 2010, 18(5): 1234-1241.

WANG Haojing, WANG Jianli, WANG Minghao, et al. Efficient Image Inpainting Based on Bilinear Interpolation Downscaling[J]. Optics and Precision Engineering, 2010, 18(5): 1234-1241.

[18]庞志勇, 谭洪舟, 陈弟虎. 一种改进的低成本自适应双三次插值算法及VLSI实现[J]. 自动化学报, 2013, 39(4): 407-417. PANG Zhiyong, TAN Hongzhou, CHEN Dihu. An Improved Low-cost Adaptive Bicubic Interpolation Arithmetic and VLSI Implementation[J]. Acta Automatica Sinica, 2013, 39(4): 407-417.

[19]BAKER S, SCHARSTEIN D, LEWIS J P, et al. A Database and Evaluation Methodology for Optical Flow[J]. International Journal of Computer Vision, 2011, 92(1): 1-31.

[20]SEVILLA-LARA L, SUN Deqing, LEARNED-MILLER E G, et al. Optical Flow Estimation with Channel Constancy[C]∥Proceedings of the 13th European Conference on Computer Vision. Switzerland: Springer, 2014: 423-438.

(责任编辑：丛树平)

A Backward Pyramid Oriented Optical Flow Field Computing Method for Aerial Image

LI Jiatian1，LI Xiankai1，LI Yingyun1，QIAN Tanghui1，LI Guojia1，LIN Yan2

1. Faculty of Land Resource Engineering，Kunming University of Science and Technology，Kunming 650093, China； 2. Policing Information Technology College，People’s Public Security University of China，Beijing 100038, China

Aerial image optical flow field is the foundation for detecting moving objects at low altitude and obtaining change information. In general，the image pyramid structure is embedded in numerical procedure in order to enhance the convergence globally. However，more often than not，the pyramid structure is constructed using a bottom-up approach progressively，ignoring the geometry imaging process.In particular，when the ground objects moving it will lead to miss optical flow or the optical flow too small that could hardly sustain the subsequent modeling and analyzing issues. So a backward pyramid structure is proposed on the foundation of top-level standard image. Firstly，down sampled factors of top-level image are calculated quantitatively through central projection，which making the optical flow in top-level image represent the shifting threshold of the set ground target. Secondly，combining top-level image with its original，the down sampled factors in middle layer are confirmed in a constant proportion way. Finally，the image of middle layer is achieved by Gaussian smoothing and image interpolation，and meanwhile the pyramid is formed. The comparative experiments and analysis illustrate that the backward pyramid can calculate the optic flow field in aerial image accurately，and it has advantages in restraining small ground displacement.

aerial image；optical flow field；backward pyramid；central projection；down sample

The National Natural Science Foundation of China (Nos. 41561082；41161061)

LI Jiatian(1975—)，male，PhD，associate professor，majors in numerical optimization and scene understanding for robot.

LI Xiankai

李佳田，李显凯，李应芸，等.航空图像光流场的逆向金字塔计算方法[J].测绘学报，2016,45(9)：1059-1064.

10.11947/j.AGCS.2016.20150367.

LI Jiatian，LI Xiankai，LI Yingyun，et al.A Backward Pyramid Oriented Optical Flow Field Computing Method for Aerial Image[J]. Acta Geodaetica et Cartographica Sinica,2016,45(9):1059-1064. DOI:10.11947/j.AGCS.2016.20150367.

P237

A

1001-1595(2016)09-1059-06

国家自然科学基金(41561082；41161061)

2015-07-13

李佳田(1975—)，男，博士，副教授，研究方向为数值最优化方法与机器场景理解。

E-mail： ljtwcx@163.com

李显凯

E-mail： 1132468376@qq.com

修回日期： 2016-06-14