卫星云图为天气预报提供云参数、大气流场和各种大气物理过程等重要的气象信息，能够直接明了地看出云系的发展和形态变化。弥补了海洋、高原和沙漠地区气象观测点少的不足，还能直观的监测到各种天气系统的变化。很多气象领域研究者利用卫星产品数据研发定点、定量、定时的预报方法，尤其是使用卫星资料和雷达产品相结合进行研究来提升预报水平，获得了大量研究成果。然而当前气象台站对卫星资料的应用水平还不高，卫星云图仅作为日常天气分析和数值预报辅助工具，甚至被忽视。针对卫星云图图像外推研究却很少，云图的图像外推工作可提供气象台和气象领域研究者所需要的下个时段云图发展的资料，具有很高的利用价值。以往的数值天气预报最新进展使我们能够以很高的分辨率预报大气动力学[1]，但对于更新周期频繁的预报应用要求来说，计算成本通常很高，效率很低。因此我们需要进一步开发计算时间短，精度高的算法。光流法作为一种新的预报算法，外推时间较短、成本较低。光流法是假设云图亮度强度保持恒定，从一系列卫星云图中跟踪云图特征点运动，基于该运动将特征点移动到未来即将到来的位置从而达到外推的目的。
根据对云图特征的假设和文献[2]提出的问题——外推方法扩展到气象卫星图像的亮度温度。将卫星云图外推大致分为两大类，第一类是欧拉持续性，只需将当前观测值等同于下时刻卫星云图数据，不用推断系统在空间和时间上的演变。第二类是拉格朗日持续性，在假定云图特征亮度强度和运动场是持续性的前提下对最近的云图特征进行外推[2-3]。此外，我们还可以根据预测不确定性对现有的预测方法进行分类：与确定性方法相反，整体临近预报试图通过运动场的不同实现和云图特征亮度强度本身的演变来说明预测的不确定性[4]。在本研究中，我们将模型开发的重点放在提供确定性的拉格朗日持续性模型组上。
拉格朗日方法包括2个步骤：跟踪和预测(外推)[5]。在跟踪步骤中，从一系列连续的卫星云图中以每个像素为基础计算速度场[2-6]。在外推步骤中，使用该速度场平移最近时刻的云图特征，即根据云图的运动场将其移动到即将到来的位置，该步骤基于半拉格朗日方案[2]、插值算法[6]、基于网格的模型实现[4]。可以在每个步骤(跟踪和预测)中使用不同的算法，从而实现云图外推过程。
“光流”是很杰出的跟踪算法，光流通常被理解为从连续图像帧中推断运动模式或速度场的一组方法[6-7]。对于速度场估计，假定亮度是恒定的，并且附加一组光流约束(OFC)条件。OFC的空间属性标志着光流模型分为两个主要类别：局部(微分)和全局(变分)[6-8]。局部模型尝试仅在某些邻域中设置OFC，而全局模型将OFC应用于整幅图像。其中还有一组独特的光谱方法，是将傅立叶变换应用于输入，并在光谱(傅立叶)域中解析OFC[9]。Bowler等首次将局部光流算法应用到预报中，为模型的发展提供了新的方向[7]，Liu等提出对卫星图像的每个像素独立使用局部Lucas-Kanade光流算法[6]，发现在红外图像外推中优于全局的HS光流算法[10]，Yeung等使用了不同的全局光流算法在香港建立了业务预报的SWIRLS产品[8]。
在过去几十年里，光流算法在临近预报中应用广泛，很多研究者把研究重点放在降水、雷暴等研究上，而忽视了大范围的气象卫星云图外推研究。本研究中我们基于用于跟踪步骤的光流公式Sparse(稀疏光流模型)和Dense(稠密光流模型)，以及图像仿射变换和空间插值的简化外推技术，通过调整其参数开发了卫星云图的外推程序，比较这2组(每组2种)光流模型在FY- 4A红外图像的外推精度，采用表现较好的光流组进行后续的预报工作。
1 研究数据
FY- 4A是我国新一代静止气象卫星，装载多种观测仪器，包括多通道扫描成像辐射计(AGRI)、干涉式大气垂直探测仪、闪电成像仪和空间环境监测仪器等。其中AGRI包含14个通道。本文研究数据是2018年FY- 4A多通道扫描成像辐射计所探测的4 km分辨率全圆盘数据，时间间隔为1 h，读取13通道的红外图像数据(地理范围为：2°30′N～62°30′N，67°30′E～150°E)进行图像外推工作。
卫星观测数据和地理信息数据是分别存放的，本文采用公式法将FY- 4A标称上的经纬度与行列号进行转换。将定标级数与定标表中每一级观测数据一一对应能读取数据[11]。从指定地理范围读取相应通道数据存储于1 201×1 651矩阵中，FY- 4A红外卫星云图示例如图1所示。
图1 FY- 4A红外卫星云图
在实验之前，对读取的每组数据依次按照[T-23，T-22，…，T-1，T，T+1，…，T+12]重新命名编排，每组数据中选用前24个时间间隔为1 h的全圆盘数据作为红外图像外推的实验数据，每组后12个数据做验证数据(实测数据)，外推出的图像是12幅时间间隔为1 h的图像。选取6组不同日期的实验数据用于本次研究，每组数据是36个时间间隔为1 h FY- 4A数据，实验数据具体如表1所示。
表1 FY- 4A基于光流模型图像外推实验组数据
起始时间结束时间持续时间时间间隔12018年4月23日01∶002018年4月24日12∶0036 h1 h22018年5月09日02∶002018年5月10日13∶0036 h1 h32018年6月10日03∶002018年6月11日14∶0036 h1 h42018年7月19日11∶002018年7月20日22∶0036 h1 h52018年9月12日21∶002018年9月14日08∶0036 h1 h62018年12月06日12∶002018年12月07日23∶0036 h1 h
2 光流模型
2.1 光流法原理介绍
2.1.1 LK光流法
光流法本质上是将真实世界的物体运动转换成相机平面的二维图像，并求解目标在图像上的运动参数u、v的方法[12]。Lucas-Kanada稀疏光流法是一种2帧差分的光流估计算法，计算得到的一种稀疏光流场[13]。其假设相邻帧之间图像亮度不变；图像中的运动随时间变化较小；像素点与领域内的点具有相似的运动，求解得到光流方程为：
Ixu+Iyv+It=0
(1)
式(1)中：I是整幅图像；Ix、Iy、It分别是对x、y、t方向求导；u、v分别为x、y方向的速度。
2.1.2 DIS光流法
Dense Inverse Search(DIS)光流法关注时间复杂度，提出了一个非常低的时间复杂度和具有竞争力的精度解决方案用于计算稠密光流[17]，它包括2个部分：
1) 快速逆向搜索。在此认为W(X;u)是在像素X上由X=(x，y)T参数化决定的图像变换，使得W(X;u)=(x+u，y+u)，逆搜索的目标函数则为
∑x[It+1(W(X;u))-T(X)]2
(2)
It+1是参考图像，T(X)是给定的图像块。通过迭代最小化找到图像变换参数，将图像变换参数更新为u←u+Δu得到:
∑x[It+1(W(X;u+Δu))-T(X)]2
(3)
使用高斯-牛顿梯度下降法，将式(3)一阶泰勒展开得到：
(4)
其中▽It+1是W(X;u)处的图像梯度。
表示图像变换的雅可比(Jacobian)矩阵。
对于参数更新为Δu，则存在一最小二乘法的闭合形式解。将式(3)的偏导数设为0得到：
(5)
其中
求解得到：
Δu=H-1∑xST·[T(X)-It+1(W(X;u))]
(6)
其中H=∑xSTS是一个n×n大小近似于海森(Hessian)矩阵。
由于S依赖于图像梯度，因此在每次迭代时必须重新计算位移u、S、海森矩阵H，所以逆向搜索算法是将参考图像与给定图像块的角色颠倒一下，从而避免每次迭代中重新计算Jacobiao和Hessian，目标函数如下：
∑x[T(W(X;Δu))-It+1(W(X;u))]2
(7)
式(7)一阶泰勒展开得到：
(8)
式(8)中：W(X;0)是初始的图像变换，▽T是W(X;0)在的图像梯度。使用最小二乘法得到：
(9)
其中
雅可比是在(X;0)处的求值，S′和H′不依赖于u，可以预先计算，没必要每次迭代中重新计算。
2) 通过多尺度图像块聚集生成稠密位移场。为了加强前后一致性，我们从It→It+1和It+1→It这两个方向并行运行算法，算法的所有步骤独立进行向前和向后光流的计算，最后合并2个方向，对应的稀疏光流的加权求和得到稠密光流场并在每个像素X上创建一个密集的向前流场
(10)
归一化后得到：
(11)
其中，当是向前位移时图像块i与参考图像位置X重叠时
当是向后位移时图像块i与参考图像位置X重叠时
和
分别表示给定图像块和变换图像之间向前和向后图像变换操作的强度差异。
和
表示图像块的向前或者向后位移。
2.2 稀疏光流模型组
这组方法的核心思想是识别FY- 4A红外云图中用于跟踪的典型特征，这些特征是卫星图像中云特征强度梯度变化显著的点，俗称”角点”。该方法不依赖于尺度大小，比将风暴看作连续对象进行跟踪经典算法更具有普遍性[14]，因为无需指定任意和尺度相关的云图特征，而角点的识别仅取决于单元领域内梯度的变化。在这一组中，使用了2个在跟踪和外推方面略有不同的模型。
2.2.1 SparseSD光流模型
SparseSD(Sparse Single Delta)读取两幅连续时刻的FY- 4A红外云图进行识别、跟踪和外推图像特征，假设T时刻是FY- 4A图像起始时间，模型参数如表2所示，具体步骤如下：
1) 使用Shi-Tomasi角点检测器[15]识别T-1时刻卫星云图中的特征。该检测器基于角质量度量公式min(λ1，λ2)来确定图像中强角点。原理是Shi-Tomasi在Harris检测算法上进行了改进，其自相关矩阵的2个特征值中较小的一个大于设定的阈值，λ1、λ2是自相关矩阵的特征值，则可以得到强角点。
2) 使用局部Lucas-Kanade光流算法[16]跟踪T时刻卫星云图上特征点(T-1时刻识别出的)并计算出大气运动矢量(AMV)。该算法使用最小二乘法求解局部特征邻域中的一组光流方程，从而在T时刻识别我们先前在T-1时刻FY- 4A红外图像中识别出的特征的位置。
3) 使用线性外推法外推特征点运动，以便预测出在每个步长时间为n的特征点位置。
4) 使用最小二乘法，根据T到T+n时刻红外图像所有已识别特征点的位置，计算每个时间步长n的仿射变换矩阵[17]。
5) 通过仿射变换外推T时刻处的FY- 4A红外图像，具体过程是将T时刻的FY- 4A红外云图每个像素点的位置通过仿射变换矩阵转换为将来T+1～T+n时刻的红外云图像素点位置。对预测图像中的剩余像素值进行线性插值，获得与T时刻相对应的FY- 4A红外图像云图特征强度图[18]。
2.2.2 Sparse光流模型
Sparse模型读取24幅时间间隔一小时的FY- 4A红外图像，仅考虑在整个时间段(24个时间步长)中保持不变的特征。模型参数如表2所示，具体步骤如下：
1) 使用Shi-Tomasi角点检测器识别T-23时刻红外图像上的特征点[15]。
2) 使用局部Lucas-Kanade光流算法跟踪T-22到T时刻红外图像上这些特征点(T-23时刻识别出的)[16]。
3) 依据时间(T-23到T)为每个成功跟踪的特征独立地建立线性回归模型。
4) 同SparseSD外推的3-5步骤。
2.3 稠密光流模型组
稠密模型组使用Dense Inverse Search(DIS)全局光流算法[17]。该算法通过对两幅连续的FY- 4A红外图像的分析估计每个图像像素点的移动速度。DIS算法与其他全局光流算法(如DeepFlow、PCAFlow)相比有更高的准确性和更高的计算效率[18-19]。
该组中的两个模型仅在外推(或平流)步骤方有所不同。第1个模型(Dense)使用恒定矢量平流方案[7]，而第2个模型(DenseRotation)使用半拉格朗日平流方案[2]。2种方法之间的主要区别在于，假定运动场本身是持续的(如图2所示)，恒定矢量方案不允许旋转运动。半拉格朗日方案允许表示大尺度旋转运动。
图2 4种平流方案的位移向量图，P是一个地块的网格点
2种平流方案的实施方式有2种可能的选择：时间上向前(在空间的下游)或时间上向后(在空间的上游)(如图2所示)。在向前方案中，从网格点P开始，将地块向下游平流直至点Q。因此，“向前”表示在时间上向前，在空间上在下游。而在向后方案中，向上游移动并确定地块的原点O，该原点将在网格点P处结束。通常，O和Q都不完全与网格点重合。在向前方案中，平流值被重新分配到相邻的网格点，而在向后方案中，需要插值方法来确定点O处的值。目前尚不清楚哪种方案可以被视为基于FY- 4A数据图像外推的最合适和通用的解决方案，一方面是关于质量守恒，另一方面是小尺度的能量损失。因此，进行了实验，将向前与向后，恒定矢量与半拉格朗日平流进行了任何可能的组合。将向后方案用作Dense和DenseRotation模型的参数。
Dense和DenseRotation模型都使用线性插值方法，以便将预测位置的卫星云图像素值插值到对应的FY- 4A卫星云图位置。稠密组模型的模型参数如表2所示。
表2 FY- 4A基于光流模型图像外推参数
模型名称输入FY-4A图像跟踪算法外推法SparseSD2幅Shi-Tomasi corner detec-tor,Lucas-Kanade optical flowConstant delta-change,affine warpingSparse24幅Shi-Tomasi corner detec-tor,Lucas-Kanade optical flowLinear regression,affine warpingDense2幅DIS optical flowBackward constant-vector advection schemeDenseRotation2幅DIS optical flowBackward semi-Lagrangian ad-vection scheme
具体实现步骤如下：
1) 基于T-1和T时刻的FY- 4A红外图像，使用全局DIS光流算法[17]计算速度场。
2) 对于T+n时刻来说根据速度场使用向后恒定向量[7]或向后半拉格朗日方案[2]逐步外推每个像素值，对于半拉格朗日方案，我们通过将速度场线性插值到该时间步的像素位置来更新每个预测时间步n处位移像素的速度。
3) 平流步骤的结果，获得了不规则速度场，由从原始位置偏移的原图像像素点组成的。使用在T+n时刻预测位置处每个位移像素的强度，对原始FY- 4A红外图像每个网格点处的像素强度进行反距离加权插值[6]。在插值之前对时间步长内每个像素进行平流处理从而将数值扩散降到最低(如Germann和Zawadzki在2002年提出的”一次插值”方法)[2]。这样，以避免插值效应使空间中的图像特征变得平滑。
2.4 持续性(Persistence)基准模型
欧拉持续性基准模型是一种启发式方法，假设云图运动场在时间和空间上是不变的，即对于任何时间步长n，云图运动场与T时刻运动场相同。这种预测技术比较粗糙，但很简单，对于很短的步长时间却是一个强大的模型，这使得欧拉持续性模型在现实世界的应用程序中是一个易于使用的选项。在目前的分析背景下，欧拉持续性型是一个基线，用来评估增加的复杂性水平是否可以提高预测技能。同时，它可以很好地衡量不同事件在时间上的去相关性。
3 光流模型评估
为了进行验证模型效果，使用4种类型的评价方式：平均绝对误差(MAE)用以评价外推图像与实际图像之间的差异，值越小差异越小。均方根误差(RMSE)用以评价预测值与真实值之间的差异，值越小模型精度越高。平均绝对百分比误差(MAPE)用以评价预测值与真实值逼近程度，值越小预测值与真实值越接近，模型精度越高。与欧拉基准模型(Persistence)比较，光流模型外推效果整体上得优于基准模型。
3.1 评价指标
使用Mean Absolute Error(MAE)用以评价预测值与真实值之间的差异。公式如下：
(14)
使用Root Mean Square Error(RMSE)用来衡量预测值与真实值之间的偏差。公式如下：
(15)
使用Mean Absolute Percentage Error(MAPE)用以计算预测值与真实值的逼近程度。公式如下：
(16)
其中simi和obsi分别指外推的卫星云图第i个像素值和实际云图的第i个像素值，n是像素数。
3.2 结果分析
研究时间内，所有模型(Sparse，SparseSD，Dense，DenseRotation，Persistence)均用程序计算1～12 h(时间间隔为1 h)的12幅外推图像。选择时间步长为3、6、9和12 h的外推图像示例如图3所示。
从图3来看，2组模型(每组2种)在外推结果上总体较好，对锋区及台风重叠区域的位置均能较为准确的推演，随着时间步长增加，外推出的图像与实际相比整体图像亮度有损失，外推云图与实际相比部分轮廓不清晰，整体上不影响实际工作中对锋区和台风影响范围的预报。由于卫星云图边缘特征点不明显，Sparse模型组在使用Shi-Tomasi角点检测器难检测出云图左上角边缘特征，并且局部检测出的特征无法代表整体云图的运动模式，随着外推步长时间增加，Sparse模型组外推的图像边缘有像素点缺失并随时间步长增加缺失数目增多。而Dense模型组外推图像整体完整，外推云图里视觉上与实际差距不大，内部轮廓不清晰。
图3 2018年7月20日09∶00 Sparse和Dense模型组图像外推图像
图4表示外推图像像素值与实际图像的差，就空间上来看外推效果，Sparse模型组和Dense模型组在图像y轴600～800区间和x轴0～500区间区域及图像右下角外推较不好，随着时间步长增加外推效果下降。由于实际图像中云团云系类型较多较杂导致两组模型都外推效果均不佳。Sparse和Dense模型组在步长时间为3时外推效果整体较好。
图4 2018年7月20日09∶00模型组图像外推图与实际图差
3.3 模型比较
本文使用评价指标为平均绝对误差(MAE)、均方根误差(RMSE)、平均绝对百分比误差(MAPE)对两组(每组两种)模型进行精度评估，选定6组连续时间的FY- 4A数据对两组(每组两种)模型图像外推效果进行评估，MAE如图5所示，RMSE如图6所示，MAPE如图7所示。
图5 Sparse模型组和Dense模型组平均绝对误差(MAE)曲线
图6 Sparse模型组和Dense模型组均方根误差(RMSE)曲线
图7 Sparse模型组和Dense模型组平均绝对百分比误差(MAPE)曲线
图5和图6分别表示平均绝对误差(MAE)，均方根误差(RMSE)。其表示外推卫星云图图像像素值与实际的差异，MAE和RMSE越小说明与实际图像差异越小，从而外推精度越高。从图6中看出：Dense模型组与欧拉持续基准对比明显优于Sparse模型组。
图7是MAPE，从图7中可以看出2组光流模型总体精度在11%之内，整体精度较高，Dense组2种模型精度优于Sparse组2种模型。
从误差图7中可以看出Dense模型组在欧拉持续性方面明显优于Sparse模型组，Sparse模型组所使用的Shi-Tomasi角点角检测器识别出的局部特征(角点)可能无法代表云图的整体运动：该检测器集中于云图特征有明显的强度和梯度变化的特征点，其运动可能不代表占主导的中尺度运动模式。Dense模型组将图像中每个像素与位移关联能够更好的外推出云图图像，Dense模型和DenseRotation之间的差异(换句话说是恒定向量和半拉格朗日方案之间的差异)可以忽略不计。半拉格朗日方案的理论优越性可能会通过持续旋转运动在其他事件(台风、气旋等预测)中体现出来。因此，应在以后的研究中进行更全面的分析。Dense和DenseRotation具有相同的插值过程，插值包括图像变形的后处理(稀疏模型)和网格化临近预报的计算中(Dense模型的一部分)。通常，这样的插值步骤会导致数值扩散，从而导致小尺度特征强度的退化或损失[2]。然而，通过在特定时间步长对预测都执行插值，我们就能够对稀疏模型和稠密模型组都包含此类影响。正如文献[2]所说，对于较长的步长时间，这些影响可能很明显，取决于所使用的外推技术。
对于图像外推而言，模型计算性能可能是满足频繁更新的用户需求的一个重要标准。在具有Intel(R)Core(TM)i7-9750H CPU(16 G内核，2.6 GHz)的笔记本上使用Jupyter notebook运行的图像外推模型程序，得到稀疏组生成12幅外推图像(1 h的时间步长)的平均时间是3.93～7.28 s，稠密组为13.6～19.1 s。因此Dense模型组的计算成本高。
4 结论
对Sparse和Dense光流模型组每组两种模型在FY- 4A红外图像外推中应用与欧拉持续性基准模型对比进行了精度评价。结果显示光流模型比以往的数值预报成本低，计算效率高，Dense模型组采用的DIS光流算法和向后平流方案外推效果优于Sparse模型组采用的Shi-Tomasi角点检测和LK光流算法方案。Dense模型组通过将图像中的每个像素点与前后两帧的图像的位移关联起来，从而外推效果较好，Dense模型组比欧拉基准模型效果更好。更能够为缺少气象站点的区域提供下个时段云图发展资料。
从平均绝对误差、均方根误差、平均绝对百分比误差评价指标上看，Dense模型组明显优于Sparse模型组，两组模型各有其优点，Sparse模型组在局部外推较好，而边缘有信息缺失，计算时间短。Dense模型组外推较为完整，计算时间长。因此综合考虑在对FY- 4A红外图像进行图像外推工作时，优先考虑Dense模型组。
[1] BAUERP，THORPE A，BRUNET G.The Quiet Revolution of Numerical Weather Prediction[J].Nature，2015，525(7567)：47-55.
[2] GERMANN U，ZAWADZKI I.Scale-Dependence ofthe Predictability of Precipitation from Continental Radar Images.Part I：Description of the Methodology[J].Monthly Weather Review，2002，130(12)：2859-2873.
[3] WOO W，WONG W.Operational Application of Optical Flow Techniques to Radar-Based Rainfall Nowcasting[J].Atmosphere，2017，8(3)：48.
[4] BERENGUER M，SEMPERE-TORRES D，PEGRAMG G S.SBMcast-An Ensemble Nowcasting Technique to Assess the Uncertainty in Rainfall Forecasts by Lagrangian Extrapolation[J].Journal of Hydrology，2011，404(3/4)：226-240.
[5] AUSTIN G L，BELLON A.The Use of Digital Weather Radar Records for Short-Term Precipitation Forecasting[J].Quarterly Journal of the Royal Meteorological Society，1974，100(426)：658-664.
[6] LIU Y，XI D G，LI Z L，et al.A New Methodology for Pixel-Quantitative Precipitation Nowcasting Using a Pyramid Lucas Kanade Optical Flow Approach[J].Journal of Hydrology，2015，529(11)：354-364.
[7] BOWLER N E H，PIERCE C E，SEED A.Development of a Precipitation Nowcasting Algorithm Based upon Optical Flow Techniques[J].Journal of Hydrology，2004，288(1/2)：74-91.
[8] CHEUNG P，YEUNG H Y.Application of Optical-Flow Techniqueto Significant Convection Nowcast for Terminal Areas in Hong Kong[C]//The 3rd WMO International Symposium on Nowcasting and VeryShort-Range Forecasting (WSN12).2012：6-10.
[9] RUZANSKI E，CHANDRASEKAR V，WANG Y.The CasaNowcasting System[J].Journal of Atmospheric and Oceanic Technology，2011，28(5)：640-655.
[10] HORN B K P，SCHUNCK B G.Determining Optical Flow[C]//Techniques and Applications of Image Understanding.International Society for Optics and Photonics，1981，281：319-331.
[11] 王清平，吴晓京，陈阳权，等.FY- 4A卫星数据可视化及应用[J].气象科技，2019 (3)：18.
Wang Q P，Wu X J，Chen Y Q，et al.FY- 4A Satellite Data Visualization and Application[J].Meteorological Science and Technology，2019 (3)：18.
[12] 周伟.基于光流法的目标提取与跟踪算法研究[D].赣州：江西理工大学，2017.
Zhou W.Research on Target Extraction and Tracking Algorithm Based on Optical Flow Method[D].Ganzhou：Jiangxi University of Science and Technology，2017.
[13] 李成美，白宏阳，郭宏伟，等.一种改进光流法的运动目标检测及跟踪算法[J].仪器仪表学报，2018，39(05)：249-256.
Li C M，Bai H Y，Guo H W，et al.A Moving Target Detection and Tracking Algorithm Based on Improved Optical Flow Method[J].Chinese Journal of Scientific Instrument，2018，39(05)：249-256.
[14] WILSON J W，EBERT E E，SAXEN T R，et al.Sydney 2000 Forecast Demonstration Project：Convective Storm Nowcasting[J].Weather and Forecasting，2004，19(1)：131-150.
[15] SHI J.Good Features to Track[C]//Proceedings of IEEE Conference on Computer Vision and Pattern Recognition.IEEE，1994：593-600.
[16] LUCASBD.An Iterative Image Registration Technique with an Application toStereoVision[C]//Proceedings DARPA Image Understanding Workshop.Vancouver，BC，Canada，24-28 August 1981，674-679，1981.
[17] SCHNEIDER P，EBERLY D H.Geometric Tools for Computer Graphics[M].Elsevier，2002.
[18] WOLBERGG.Digital Image Warping[M].IEEE Computer Society Press，Los Alamitos，CA，USA，1990.
[19] KROEGER T，TIMOFTE R，DAI D，et al.Fast Optical Flow Using Dense Inverse Search[C]//EuropeanConference on Computer Vision.Springer，Cham，2016：471-488.
[20] WEINZAEPFELP，REVAUD J，HARCHAOUI Z，et al.DeepFlow：Large Displacement Optical Flow with Deep Matching[C]//IEEE International Conference on Computer Vision.IEEE，2014.
[21] WULFF J，BLACKM J.Efficient Sparse-To-Dense Optical Flow Estimation Using a Learned Basis and Layers[J].2015(4)：120-130.
[22] BARRON J L，FLEET D J，BEAUCHEMIN S S.Performanceof Optical Flow Techniques[J].InternationalJournal of Computer Vision，1994，12(1)：43-77.
[23] OTTE M，NAGEL H H.Estimation of Optical Flow Based on Higher-Order Spatiotemporal Derivatives in Interlaced and Non-Interlaced Image Sequences[J].Artificial Intelligence，1995，78(1/2)：5-43.
[24] 孔丁丁.图像序列光流计算评估理论及方法研究[D].南昌：南昌航空大学，2015.
Kong D D.Research on the Theory and Method ofImage Sequence Optical Flow Calculation and Evaluation[D].Nanchang：Nanchang Hangkong University，2015.