基于改进光流法的雷达图像运动估计

王志斌，肖艳姣，吴 涛

(1.中国气象局 武汉暴雨研究所，湖北 武汉 430205；2.武汉中心气象台，湖北 武汉 430074)

基于改进光流法的雷达图像运动估计

王志斌1，肖艳姣1，吴 涛2

(1.中国气象局 武汉暴雨研究所，湖北 武汉 430205；2.武汉中心气象台，湖北 武汉 430074)

介绍了一种改进的光流方法对天气雷达图像运动进行估计，计算天气雷达图像运动的光流变化，从而得到矢量场。用传统的Horn & Schunck方法中的全局平滑化算子估计时，其对数据的噪声敏感，且难以达到全场最优条件。通过加入改进后的高价平滑算子，可以减少靠向二级最小值的可能。另外，通过加入Lucas & Kanade方法计算局部光流得到小范围的光流场，在小的区域内光流场容易满足最优条件。适当运用平滑处理得到整场的矢量场，并以它作为改进Horn & Schunck方法的初猜场进行计算，且应用二维连续无辐散质量方程约束，同时对得到的风场进行风速以及风向的质量控制，得到比较完整的雷达运动矢量场。通过实验分析发现，改进的光流法优于气象上传统的交叉相关法，且优于单一使用某种光流方法的结果。

雷达图像;运动估计;光流法;约束;改进

1 概 述

多普勒天气雷达图像中包含了丰富的天气现象，通过图像可以对大风、冰雹、雷电、暴雨等强天气进行监测，也可通过多幅图像的变化对它们进行预报[1]。目前气象预报业务上使用自动外推预报技术，包括两类：交叉相关法[2-3]和单体质心法[4-7]。单体质心法是将雷达的三维数据进行分析；交叉相关法利用不同时次天气雷达图像求最优空间的相关，选取不同时次相关最好的空间匹配块，来确定图像的移动矢量特征，并进行图像外推预报。交叉相关算法计算简单，可以求不同时刻的相关系数，也可以求欧氏距离。但对形状变化较快的天气雷达图像，通过交叉相关法计算出的运动矢量场质量会降低[8]。

为了克服这些缺点，引入了在计算机领域运用较多的光流方法。采用改进的光流法计算得到天气雷达图像的光流场来代替交叉相关法得到的风场(后面提到的光流场和风场及矢量场是等价的)。光流法最早由Gibson[9]提出，Horn & Schunck[10]提出了一种基于全局约束的求解方法(简称HS方法)。Lucas-Kanade[11]提出了一种局部的约束求解方法(简称LK方法)。基于光流的基本约束方程，还产生了其他光流方法，但大多数可以用变分模型来描述，即最小化某个能量泛函的过程。这个泛函可以由两项组成，即平滑项和数据项。数据项和图像本身的数据有关，如雷达回波数据、图像的灰度等，数据项决定了运动的真实程度。平滑项是泛函的附加项，进行条件约束，主要解决方程的不适定性。确定了泛函方程后，可以运用欧拉-拉格朗日方程进行求解，其数值解一般采用松弛迭代或最优下降算法求出。光流法在计算机图像处理、模式识别、交通、医学等领域[12-14]应用广泛。韩雷[8]运用HS方法进行了外推预报，其不足之处是全场不容易满足最优条件，容易陷入局部最小化，特别是在天气雷达图像的场能量不集中时。曹春燕[15]使用LK方法，运用局部最优进行求解计算，对局部容易满足最优条件[16-18]，但对整场不容易满足最优条件。

文中将两者结合起来，先运用LK求局部最优；再以局部最优结果作为初始风场，利用改进HS方法计算全场光流值，得到相对合理的运动矢量场；最后通过二维连续质量方程进行约束，得到更加完善的运动矢量场。与传统的交叉相关法和单一的光流法相比较，结果表明该方法更合理有效，在气象的短时预报业务中得到了较好的应用。

2 光流方法改进及天气雷达数据的质量控制

光流方法的原理是，光流是空间运动物体在平面上投影时运动的速度[19-20]，也是研究图像灰度在时间上的变化，因此可以将光流矢量定义为投影平面特定坐标点上的灰度瞬时变化率，代表二维矢量场。二维光流场实际上是实际物体的三维运动场的投影[21-22]。光流法的实质是由二维光流场重构三维运动场。下面介绍两种经典的光流方法：HS和LK。

设在平面坐标上有一点M(x,y)，它代表了三维空间中一点V(x,y,z)在平面坐标中的投影。该点在t时刻的灰度为I(x,y,t)。另假设该点在t+Δt时刻运动到(x+Δx,y+Δy),在一定时间间隔Δt内灰度值保存不变。即

I(x,y,t)=I(x+Δx,y+Δy,t+Δt)

(1)

将式(1)按照泰勒公式展开：

(2)

其中，Rn包含了Δx,Δy,Δt的二次以上的无穷阶小项。

变换式(2)即为：

(3)

式(3)除以Δt，并且取Δt→0的极限，忽略掉高阶项Rn可得：

(4)

可以简写为：

Ixu+Iyv+It=0

(5)

式(5)就是光流方程。Ix,Iy为图像灰度的空间梯度；It为图像灰度的时间变化率；u,v为光流场。在方程中，Ix,Iy,It可以通过相邻的两幅图像的灰度值计算出来，但式中有两个变量，而只有一个方程，如果不做任何假设，方程式无法求解。因此，需要引入更多的约束条件，才能求解u,v。

2.1 基于全局的HS算法及其改进

如果相隔时间很短的两幅图像且灰度变化也很小，在光照不变的条件下，Horn与Schunk等导出了光流的基本方程。巧妙地将二维运动和图像灰度相关联。这个方程有两个假设，灰度守恒和平滑约束，并且提出了全局约束条件方程。通过变分方法可以求解该方程。全局约束方程如下：

J=J0+α2JHS

(6)

其中，α为平滑系数；J0为灰度守恒项；JHS为平滑约束项。

J0=∬(Ixu+Iyv+II)2dxdy

(7)

(8)

根据HS约束方程，用欧拉-拉格朗日方程求解,可得：

Ix(Ixu+Iyv+It)-α2Δu=0

(9)

Iy(Ixu+Iyv+It)-α2Δv=0

(10)

其中，Δ是二阶的拉普拉斯算子。

(11)

用雅克比迭代法可以求解式(9)和式(10)。

(12)

(13)

用上述方法求解光流场，由于采用全局优化方案，很容易陷入局部最小值，而且对图像噪声敏感。因此改平滑项为其他方法，即改用Wahba and Wendelberger[23]方法(简称Jww)。这可以减少靠向二级最小值的可能。即JHS改为Jww：

(14)

同样可以用高阶欧拉-拉格朗日方程求解。

Ix(Ixu+Iyv+It)-α2Δ2u=0

(15)

Iy(Ixu+Iyv+It)-α2Δ2v=0

(16)

其中，Δ2u,Δ2v为：

(17)

(18)

同样可以用雅克比迭代法求解。方程如下：

(19)

(20)

图1中描述了各格点。其中f0表示为ui,j,i,j为对应的格点坐标。

i-2i-1ii+1i+2j+2f12j+1f8f4f5jf11f3f0f1f9j-1f7f2f6j-2f10

图1 格点表示

2.2 LK方法

在HS方法中，采用Jww来改进数据平滑项，但仍然有可能不是全局最优，特别是数据噪声大时。因此LK提出了局部最优化的方法。基本思路是把图像分为多个小区域，因为在小范围内容易满足光流条件，抗噪声也强。

假设一个小区域范围内的光流恒定，根据光流方程有：

(21)

式中有两个未知变量u,v，但方程多余两个，是一个超定方程，需要用最小二乘法求解。

将式(15)记为如下行列式：

(22)

(23)

(24)

得到u,v的解为：

(25)

如果在选择的小区域范围内使用一个窗口权重函数，使得靠近中心区域的影响比其他范围大，得到的速度场为：

(26)

其中，W为权重函数，用最小二乘法求解可得：

(27)

2.3 HS-LW混合方法求解

LW方法对一个小区域内进行光流计算，HS是对整场进行最优计算，将两者结合起来，先用LW计算出全场的光流场，再用这个光流场作为HS方法的初始值并用HS方法计算出全场的光流。计算出的光流场相对合理。

利用HS和LW方法进行灰度导数计算时，用Sobel算子进行计算。Sobel算子有对噪声处理的鲁棒性。LK方法中的窗口函数为杨辉三角系数，一般取n=9。

2.4 无辐散质量方程约束

最后，用二维连续无辐散质量方程对HS-LW方法求出的矢量场进行约束。

二维连续无辐散质量方程[24]为：

(28)

同样用变分的方法，对HS-LW求出的矢量场u0,v0和实际的矢量场u,v用一个代价函数表示，并最终求解出实际矢量场。

代价函数为：

(29)

根据欧拉-拉格朗日方程求解得:

(30)

(31)

边界条件为λ(x,y)=0。

将式(30)和式(31)代入式(28)得到：

(32)

式(32)可以通过有限差分方法计算出λ。同理也可以通过式(30)和式(31)计算出u,v。

2.5 数据处理质量控制及程序计算流程

文中使用的资料为湖北省内的6部多普勒天气雷达资料。在进行拼图前，对单部天气雷达图像数据进行了质量控制，方法采用文献[25]提出的算法进行。首先针对不同雷达型号对数据进行规整和缺测填补，对原始数据进行有效性检查，对孤立回波进行过滤。然后选取了五个回波特征量作为候选因子，对超折射回波进行抑制，计算出超折射回波和正常回波的概率分布图，找出两者的差别，建立各因子的隶属函数。最后通过模糊逻辑算法，把这些因子进行综合运算，得到控制值，大于某值则去掉该回波，并用周围正常回波进行替代，否则该值正常不进行处理。

进行单站质量控制后，把极坐标的单站雷达资料转换成三维的直角坐标格式，采用文献[26]提出的方法进行。取3 km高度的CAPPI拼图作为光流场的输入灰度值。在进行了单站质量控制的基础上，为了进一步减少数据噪声，仍用9点平滑方法对3 km的CAPPI场进行质量控制。

通过上述各种光流方法计算出的光流场，有些格点的值仍有不合理的地方，需要进一步加上附加条件进行约束和控制。首先对计算出的光流场某些格点值大于50 m/s的值进行剔除。将某一格点和其周围格点(一般是9*9格点范围)的平均值进行比较。当超过平均速度6 m/s，以及夹角超过30°，这个值也剔除，并用9*9格点范围内的平均值进行替换，同时这个值和整场的平均值偏差大于15 m/s时给予剔除。

3 实 验

为了评价各种光流方法的效果，把各种光流算法和TREC计算出的风场对天气雷达图像运用欧拉后向外推算法进行外推预报，并用气象上三个通用指标进行判断。三个指标分别为击中率(POD)、虚警率(FAR)和临界成功指数(CSI)。

(33)

(34)

(35)

将预报结果和实况位置的雷达图像数据进行比较，预报出现而实况也出现为成功；实况出现但预报没有出现为失败；预报出现但实况没有出现为虚假。外推时间为30 min和60 min。把雷达图像数据数值分为5个等级(5～15 dbz，15～25 dbz，25～35 dbz，35～45 dbz，45 dbz)进行评分。选取3个不同的影响湖北省的降水天气过程，分别是2016年4月15日到16日的层状云降水，2016年6月30日到7月1日以及2016年7月3日到7月4日对流性降水过程。

从表1和表2可以看出，HS+LW+二维约束的效果略好，无论是击中率还是成功指数都比其他方法好，虚警率也较低，其外推30 min的成功指数比TREC高5个百分点。其他光流方法和TREC的效果相当，指标表现上各有千秋，外推60 min的效果明显不如外推30 min的，说明预报时间越长，效果越差。

表1 2016年4月15日-16日30 min外推结果平均

表2 2016年4月15日-16日60 min外推结果平均

从表3和表4可以看出，HS+LW+二维约束在击中率、成功指数以及虚警率方面比其他方法好，其30 min击中率比TREC高9个百分点，成功指数高6个百分点。除HS方法和TREC相当外，其他光流方法比TREC效果略好。

表3 2016年6月30日-7月1日30 min外推结果平均

表4 2016年6月30日-7月1日60 min外推结果平均

从表5和表6可以看出，HS+LW+二维约束在击中率、成功指数以及虚警率方面都比其他方法略好，其30 min击中率比HS高8个百分点，成功指数高7个百分点。除HS方法和TREC相当外，其他光流方法比TREC效果略好 。

表5 2016年7月3日-7月4日30 min外推结果平均

表6 2016年7月3日-7月4日60 min外推结果平均

天气雷达图像光流估计处理流程见图2。

图2 雷达图像光流计算处理流程

4 结束语

介绍了光流方法的基本原理，提出了一种改进的HS方案。结合HS及LW各自优点联合计算出光流场，通过二维连续方程进行约束。通过天气雷达图像的变化计算出实用的风场，并利用外推方法对各种结果进行了比较。

改进方法在对流性降水时优于TREC，而在层状云降水时差别不明显，只是HS+LW+二维约束方法略好于其他方法。通过改进的HS方法，进一步减少了计算结果趋向局部最小的可能性，结果优于HS。通过使用LW方法产生的初始光流的结果，作为改进的HS方法的初始场，使得原来使用单一的LW方法的结果得到了进一步优化。通过结合改进的HS、LW和二维连续方程约束等方法，得到了效果较好的光流结果。

与传统的TREC方法相比，光流法侧重于变化，因为在计算光流场时不仅考虑了天气雷达图像在空间上的变化，同时也考虑了时间上的变化。该中和方法计算耗时小，在一般微机上几秒即可完成，因此有较大的应用推广价值。

[1] Wong W K.Merging radar nowcast with convection-permitting rapidly-updated NWP model for significant convection forecast[C]//4th international symposium on nowcasting and very-short-range forecast.[s.l.]:[s.n.],2016.

[2] Rinehart R E,Garvey E T.Three-dimensional storm motion detection by conventional weather radar[J].Nature,1978,273(5660):287-289.

[3] Wilson J W,Crook N A,Mueller C K,et al.Nowcasting thunderstomes:a status report[J].Bulletin of the American Meteorological Society,1998,79(10):2079-2100.

[4] Crane R K.Automatic cell detection and tracking[J].IEEE Transactions on Geoscience Electronics,1979,17(4):250-262.

[5] Dixon M, Wiener G. TITAN: Thunderstorm identification, tracking,analysis,and nowcasting:a radar-based methodology[J].Journal of Atmospheric & Oceanic Technology,1993,10(6):785-797.

[6] Rasmussen R,Dixon M,Hage F,et al.Weather support to deicing decision making (WSDDM):a winter weather nowcasting system[J].Bulletin of the American Meteorological Society,2001,82(4):579-595.

[7] Mueller C,Saxen T,RobertsR,et al.NCAR auto-nowcast system[J].Weather & Forecasting,2003,18(4):545-561.

[8] 韩 雷,王洪庆,林隐静.光流法在强对流天气临近预报中的应用[J].北京大学学报:自然科学版,2008,44(5):751-755.

[9] Gibson J J.The perception of the visual world[M].[s.l.]:[s.n.],1950.

[10] Horn B K P,Schunck B G.Determining optical flow[J].Artificial Intelligence,1981,17(1-3):185-204．

[11] Lucas B D,Kanade T.An iterative image registration technique with an application to stereo vision[C]//Proceedings of seventh international joint conference on artificial intelligence.Vancouver:IEEE,1981.

[12] 史金龙,白素琴,施金宛,等.基于CLG光流算法的云的运动分析与研究[J].计算机技术与发展,2011,21(12)：135-138.

[13] 马 龙,王鲁平,李 飚,等.自然场景图像的光流场估计[J].系统工程与电子技术,2012,34(6):1278-1282.

[14] 柳士俊,张 蕾.光流法及其在气象领域里的应用[J].气象科技进展,2015,5(4):16-21.

[15] 曹春燕,陈元昭,刘东华,等.光流法及其在临近预报中的应用[J].气象学报,2015，73(3)：471-480.

[16] 杨国亮,王志良,牟世堂,等.一种改进的光流算法[J].计算机工程,2006,32(15)：187-188.

[17] 陈晓静，敬忠良，张 军.基于卫星图像数据的大范围风场重建[J].计算机工程,2013，39(12)：233-236.

[18] 魏国剑，侯志强，李 武,等.融合光流检测与模板匹配的目标跟踪算法[J].计算机应用研究,2014，31(11)：3498-3501.

[19] 郑来芳,孙 炜,欧阳明华,等.结合光流和人工势场的风管机器人避障方法[J].计算机工程与应用,2016，52(9)：243-247.

[20] 吴振杰，毛晓波.一种改进的帧间差光流场算法[J].计算机工程与应用,2013,49(18):200-203.

[21] 陈 震,高满屯,曲仕茹,等.基于直线光流场的三维运动和结构重建[J].计算机学报,2004,27(8):1046-1055.

[22] 宋 爽,杨 健,王涌天.全局光流场估计技术及展望[J].计算机辅助设计与图形学学报,2014,26(5):841-850.

[23] Wahba G,Wendelberger J.Some new mathematical methods for variational objective analysis using splines and cross validation[J].Monthly Weather Review,1980,108:1122-1143.

[24] Li L,Schmid W,Joss J.Nowcasting of motion and growth of precipitation with radar over a complex orography[J].Journal of Applied Meteorology,1995,34:1286-1300.

[25] 吴 涛,万玉发,沃伟锋,等.SWAN系统中雷达反射率因子质量控制算法及其应用[J].气象科技,2013,41(5):809-817.

[26] 肖艳姣,刘黎平.新一代天气雷达网资料的三维格点化及拼图方法研究[J].气象学报,2006,64(5):647-657.

MotionEstimationforRadarImageBasedonImprovedOpticalFlowMethod

WANG Zhi-bin1，XIAO Yan-jiao1，WU Tao2

(1.Institute of Wuhan Heavy Rain,China Meteorological Administration,Wuhan 430205,China；2.Wuhan Central Meteorological Observatory,Wuhan 430074,China)

An improved optical flow method is introduced to estimate the radar images motion and calculate the optical flow change of motion for radar image,obtaining the vector field.It is sensitive to data noise and difficult to meet the optimal condition when global smoothing operator from traditional Horn & Schunck method is used to estimate.After adding the improved smoothing operator with high level,it can reduce the possibility for the secondary minimum.In addition,through Lucas & Kanade method for calculation of the local optical flow,a small range of optical flow field is obtained,which can meet the optimal requirements.Smooth processing is applied to get the entire vector field appropriately,with it as first guess filed of improved Horn & Schunck method for calculation.At the same time,through the two-dimensional continuous quality equation constraints without radiation,controlling the wind speed and direction of wind field acquired,the more complete motion vector field for radar is obtained.Experimental analysis shows that the improved optical flow method is better than traditional cross-correlation method and a single optical flow method.

radar images;motion estimation;optical flow method;constraints;improvement

TP39

A

1673-629X(2017)12-0170-06

10.3969/j.issn.1673-629X.2017.12.037

2016-10-20

2017-02-23 < class="emphasis_bold">网络出版时间

时间：2017-08-01

国家自然科学基金资助项目(41275106)；科技部行业专项公益性行业(气象)科研专项(GYHY201306008)；湖北省雷电专项( FL-Z-201401)

王志斌(1965-)，男，正研高工，研究方向为计算机图像处理、数据库开发。

http://kns.cnki.net/kcms/detail/61.1450.TP.20170801.1550.020.html