雷达图像运动矢量场计算

伊凡

（中国海洋大学 信息科学与工程学院，山东 青岛 266100）

多普勒天气雷达是强对流天气监测和预警的主要工具。因此对雷达图像的运动矢量场的准确计算在天气预报中起了很重要的作用，也对减少雷雨暴风等自然灾害给国家和人民带来的损失有重要贡献。对于雷达图像运动矢量场的计算有多种方法，主要方法有：基于梯度的方法、基于匹配的方法、基于能量的方法、基于相位的方法[1]。它是一种基于梯度的全局平滑约束算法。而在1981年Lucas和Kanade又提出了一种的基于局部平滑约束的算法，后来的研究中称其为L-K算法。在前人工作的基础上，本文将详细分析如何利用L-K算法具体计算雷达图像的运动矢量场的问题。

1 算 法

算法部分将主要简要介绍一下计算雷达图像运动矢量场过程中所涉及的主要算法并对如何计算雷达图像矢量场提出了一种新的详细算法。

1.1 全局约束光流算法

1981年，Horn等人在相邻图像的间隔时间很小，其图像中灰度变化也很小的前提下，导出了灰度图像光流场计算的基本等式。设图像平面上点（x，y）在 t时刻的图像灰度为 I（x，y，t）。当模式运动的时候一个特殊点的灰度将是一个常量。即

利用泰勒展开式展开，最终会得到：

也就是光流基本等式[2]。其中，

1.2 L-K算法

L-K算法的目的就是使误差和的平方最小化[3]。Lucas和Kanade假设在一个小的空间邻域Ω上运动矢量保持恒定[4]，然后使用加权最小二乘法（weighed least-squares）估计光流。在一个小的空间邻域Ω上，光流估计误差定义为：

上式中Ω代表以P点为中心的一个小区域，W2（x）为窗函数，代表区域中各点的权重，离P点越近，权重越高[5]，V=（u，v）T是中心点P的光流。上式的解可以由下面的方程得到：

对于在邻域Ω内的n个点xi，其中，

最后，方程的解为：

上式中所有的求和都是在Ω的所有点上进行的。对于Ω的选取需要经过多次实验找出最佳值。

1.3 雷达图像运动矢量场的计算

如图1所示将雷达图像划分ng×ng个子块，假设图像的大小为M×N，则每个子块的大小为，每个子块的顶点如图中黑点所示。取Ω区域的大小为相邻上下左右4个子块的大小，如图1中4种颜色所圈的4个不同区域，每个颜色代表一个Ω区域。

图1 L-K算法原理图Fig.1 L-k algorithm principle diagram

首先，需要计算一下前一时刻和当前时刻这两幅相邻时刻输入图像之间的平均运动矢量。然后需要计算分块后的图像每个网格点上的运动矢量。

第一，将输入的两幅相邻时刻雷达图像所有像素点都进行拉普拉斯平滑处理，平滑次数设为10次。然后对这两幅输入图像进行一个简单的阈值处理，即将灰度小于某一阈值的像素点的灰度设为零。此处阈值可以设为15 dbz。这样就可以将一些不必要的点进行滤除，从而减少运算次数，节省运算时间。此时就完成了对输入图像的预处理。

第二，计算区域内非零像素点的个数，若非零像素点的个数小于100，将区域内的光流全部设为零。若非零点的个数大于100，则使用L-K局部光流算法计算每个区域中心点处的光流。 图1中 4个不同区域的中心点分别为 （x1，y1），（x2，y2），（x3，y3），（x4，y4）。 其中利用使小区域上下左右各移动 5 个像素点的方法来计算x和y方向上的亮度梯度。由此可获得除边界子块外其他子块的顶点上的光流。

第三，边界上点的值等于边界上右下角顶点的值，比如边界点K的值等于对角点L的值，边界点H的值等于对角点I的值。对于缺失点的光流速度用平均光流速度来代替。如果计算出的光流的值为无效值，也将其用前面计算的平均运动矢量来代替。

最后对整幅图像的光流速度进行1次平滑处理。然后用双线性插值法对整幅图像的运动矢量进行填充，假设某一子块的大小为 m×n，4 个顶点的光流值分别为：V（x，y），V（x，y+m），V（x+n，y），V（x，n，y+m），则子块中任一点（x′，y′）的光流为：

其中，

利用此方法既保持了图像的规律又有利于保护边缘，这样即可求出图像中每一点的光流，然后各点的光流乘以两幅图像的相隔时间进一步计算出图像运动的矢量场[3]。

2 实验结果分析

经过多次实验分析，最终得出以下结果：关于将图像分块的ng的大小，实验过程中采用了多个不同的ng进行实验，最终确定当ng=35时，实验结果最为理想。本次试验是采用的天津地区雷达站2005年6月14的数据。图2.a和图2.b分别为00:15和02:21时刻的雷达矢量运动场图像。

图2 天津雷达2005年6月14日00:15和02:21的雷达图像矢量场的计算结果Fig.2 Tianjin radar on June 14， 2005，00:15 and 02:21calculation results of radar image vector fields

由图2可以看出，计算出的雷达图像矢量场非常平滑，同时光流场的方向也具有较高的准确性。

3 结束语

文中总结了在计算雷达图像矢量场的过程中所用到的两种主要光流算法即全局光流算法和L-K算法的主要内容。在此基础上提出了一种新的计算雷达运动矢量场的详细算法。因为在计算过程中采用将图像分块的方式进行计算光流[7]，不需要对图像上每一个像素点计算，因此，本文的算法可以有效节省运算时间。并经过试验表明，这种算法能很好的计算出雷达图像[8]的运动矢量场。

[1]杨国亮，王志良，牟世堂，等.一种改进的光流算法[J].计算机工程，2006，32（15）：.YANG Guo-liang， WANG Zhi-liang， MU Shi-tang， et al.An improved optical flow algorithm[J].Computer Engineering，2006，32（15）：.

[2]Horn B K P，Schunck B G.Determining Optical Flow[J].Artificial Intelligence，1981（17）:185-204.

[3]Bakeb S，Matthews I.Lucasdetermining optical flow.artificial intelligencekanade 20 years On:a unifying frameworke[J].Internationl Journal of Computer Vision，2004，56 （3）：221-255.

[4]Lucas B，Kanade T.An interative image registration technique with an application to stereo vision[J].Proceedings of the 7th International Jion conference on Articial Intelligencee，1981（2）:121-130.

[5]Barron J L，Fleet D J，Beauchemin S S.Systems and experiment performance of optical flow technique[J].International Journal of Computer Visione，1994，12（1）:43-77.

[6]王敏.基于多普勒雷达图像的风暴运动分析 [D].青岛:中国海洋大学，2009.

[7]林存花，陈海峰.光流场模型用于非刚性医学图像配准[J].电子科技，2011（4）：53-56.LIN Cun-hua，CHEN Hai-feng.Optical flow field model for non-rigid medical image registration[J].Science and Technology，2011(4)：53-56.

[8]刘强，刘忠义，杨泽刚，等.基于可编程渲染管线的雷达图像分层模型设计与实现[J].现代电子技术，2009（5）：50-52.LIU Qiang，LIU Zhong-yi，YANG Ze-gang，LIU Ping.Design and realization of radar image delaminaton based on programmable render pipline[J].Modern Electronics Technique，2009（5）:50-52.