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

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

1 逆向金字塔

图像金字塔是一系列以金字塔形状排列的分辨率逐步降低的图像集合，金字塔的底部是原始图像的高分辨率表示，而顶部是低分辨率的近似^[12-13]。当向金字塔的上层移动时，尺寸和分辨率降低。在不考虑图像的几何成像过程时，图像金字塔可依据分辨率因子γ以γ_f⁰→γ_f¹→γ_f²→…→γ_fⁿ形式对原始图像降采样并形成n+1层金字塔，此处，建立过程是由下至上的，并且分辨率是逐步减少的，称这种方式为正向金字塔(forward pyramid)，如图 1(c)所示。根据中心投影模型描述地物目标的物方(object side)、像方(image side)成像关系，以金字塔结构能够有效地反映物方设定的最小位移为条件，首先需要定量地得出金字塔中最小尺寸图像的降采样因子，进而根据等比关系确定中间层图像分辨率。由于在这种金字塔的构建过程中需要首先确定顶层图像的分辨率，即γ_b⁰→γ_bⁿ→γ_b¹→…→γ_b^n-1，所以被称为逆向金字塔(backward pyramid)，如图 1(b)所示。

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

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

设地面某点在前后两时刻的位置分别为a_o、a′_o，依据式(1)，则该点在像方的位置分别为a_i、a′_i。相片成像是以像素(pixel)为单位并且光流向量由像素点移动构成，因此，光流向量在像方的最小长度应大于1。式(1)中，置a_ib_i的长度值大于1，则有式(2)

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

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

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

金字塔中图像是原始图像的降采样结果，因此，图像所对应的物方位移比c′₀=αc₀，0 < α < 1，则式(4)变为

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

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

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

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

1.2 逆向金字塔构建

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

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

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

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

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

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

双线性插值计算公式如下

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

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

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

2 试验与分析 2.1 光流场计算过程

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

式中，I_i为逆向金字塔中的第i层，I₀为底层，即原始图像，而I_n为顶层。

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

输入：图像I¹、图像I²，物方最小可分辨位移值λ。

输出：光流向量。

(1) 根据λ分别构建I¹与I²的逆向金字塔BP¹与BP²。

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

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

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

(a) 用有限差分法计算I_i²在x和y方向的梯度图像I_x与I_y，I_i²∈BP²。

(b) 计算I_i¹和I_i²在时域上的差值I_t，I_i¹∈BP¹，I_i²∈BP²。

(c) 用3×3模板，计算每个像素点光流向量的局部平均值u、v。

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

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

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

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

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

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

(6) 算法结束。

2.2 逆向金字塔作用分析

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

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

设λ为11.1 m，即光流向量可以反映出地面目标11.1 m及以上位移大小的变化。因为曝光时间间隔为2 s，所以对于车辆目标来说，相当于如果其行驶速度大于20 km/h，即可以计算出光流。设ε等于0.001，根据式(4)与式(7)，得γ_min为0.01。由于航空图像分辨率较高，设置金字塔为5层，由式(9)得γ_b⁴=γ_min=0.01、γ_b¹=0.316、γ_b²=0.1、γ_b³=0.032，用γ¹-γ⁴降采样因子构建逆向金字塔。通常正向金字塔是按0.5倍的方式进行图像降采样^{[7-8, 19-20]}，得出各图层的采样比率为γ_f¹=0.5、γ_f²=0.25、γ_f³=0.125、γ_f⁴=0.062 5。原始图像如图 2(a)、(b)所示，逆向金字塔与正向金字塔的光流计算结果如图 2(c)、(d)所示。

原始图像中包含对向行驶的两辆汽车，速度均大于20 km/h，由于γ_f⁴>γ_b⁴，所以在几何成像过程上，正向金字塔比逆向金字塔具有更小的地面目标位移分辨率。因此，逆、正向金字塔均能够计算得到汽车位移的光流向量，如图 2(c)、(d)中的光流线段所示。在γ_f⁴>γ_b⁴情况之下，逆、正向金字塔光流计算的差异主要表现在抑制地面目标微小位移。产生微小位移的原因主要有两方面构成，一是风会造成树木等地面目标摆动并在图像上产生位移，二是环境光照的变化会在以像素点为单位的光流计算过程中产生位移。本文算例中，环境光照无变化，主要是以风造成的树木枝叶摆动位移为主，如图 2(d)所示。经统计，非车辆位移在原始图像中的光流表现为，光流向量最小长度3像素，最大长度5像素，平均长度3.5像素，这说明当地面微小位移变化能够被正向金字塔顶层图像的像素值所反映。与此相反，在逆向金字塔中，γ_b⁴由几何成像过程严格计算得出，其基础是以反映大位移λ为条件，由于γ_min仅为0.01，经过高斯平滑与插值过程使得原始图像中的微小位移被充分过滤掉，在顶层图像中没有反映，如图 2(c)所示。逆向金字塔在地面目标大位移情况之下的表现是优于正向金字塔的，其能够有效地规避地面微小位移的影响，更便于后续分析与处理。

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

设λ为0.5 m，即光流向量可以反映出地面目标0.5 m及以上位移大小的变化。因为曝光时间间隔为2 s，所以对于车辆目标来说，相当于如果其行驶速度大于0.9 km/h，即可以计算出光流。根据式(4)与式(7)，得γ_min为0.2。同样设置金字塔为5层，得γ_b⁴=γ_min=0.2、γ_b¹=0.669、γ_b²=0.447、γ_b³=0.299，用γ¹-γ⁴降采样因子构建逆向金字塔。用降采样比率0.5建立正向金字塔，即γ_f¹=0.5、γ_f²=0.25、γ_f³=0.125、γ_f⁴=0.062 5。分别用逆向金字塔与正向金字塔进行光流计算，结果如图 3所示。

原始图像是位于清水江边的某处煤炭转运站，货车正处于装煤的小位移位置调整中，原始图像如图 3(a)、(b)所示。逆向金字塔计算结果如图 3(c)所示。经统计，光流向量最小长度为7像素，最大长度为11像素，平均长度为8.7像素，光流向量可以正确地反映发生小位移变化的3辆货车。与此相反，正向金字塔并未计算出光流。由于γ_f⁴为0.062 5，由中心投影模型可知，能够被正向金字塔顶层图像所反映的地面目标最小位移为1.6 m，因为货车的最大位移仅为1.1 m，因此，理论分析与试验结果得到相互印证。逆向金字塔在地面目标小位移情况之下的表现同样是优于正向金字塔的，其能够准确地生成地面目标小位移的光流向量，更有利于变化检测分析。

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

3 结论

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