雷達圖像運動矢量場計算

伊凡

（中國海洋大學 信息科學與工程學院，山東 青島 266100）

多普勒天氣雷達是強對流天氣監(jiān)測和預警的主要工具。因此對雷達圖像的運動矢量場的準確計算在天氣預報中起了很重要的作用，也對減少雷雨暴風等自然災害給國家和人民帶來的損失有重要貢獻。對于雷達圖像運動矢量場的計算有多種方法，主要方法有：基于梯度的方法、基于匹配的方法、基于能量的方法、基于相位的方法[1]。它是一種基于梯度的全局平滑約束算法。而在1981年Lucas和Kanade又提出了一種的基于局部平滑約束的算法，后來的研究中稱其為L-K算法。在前人工作的基礎(chǔ)上，本文將詳細分析如何利用L-K算法具體計算雷達圖像的運動矢量場的問題。

1 算 法

算法部分將主要簡要介紹一下計算雷達圖像運動矢量場過程中所涉及的主要算法并對如何計算雷達圖像矢量場提出了一種新的詳細算法。

1.1 全局約束光流算法

1981年，Horn等人在相鄰圖像的間隔時間很小，其圖像中灰度變化也很小的前提下，導出了灰度圖像光流場計算的基本等式。設(shè)圖像平面上點（x，y）在 t時刻的圖像灰度為 I（x，y，t）。當模式運動的時候一個特殊點的灰度將是一個常量。即

利用泰勒展開式展開，最終會得到：

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

1.2 L-K算法

L-K算法的目的就是使誤差和的平方最小化[3]。Lucas和Kanade假設(shè)在一個小的空間鄰域Ω上運動矢量保持恒定[4]，然后使用加權(quán)最小二乘法（weighed least-squares）估計光流。在一個小的空間鄰域Ω上，光流估計誤差定義為：

上式中Ω代表以P點為中心的一個小區(qū)域，W2（x）為窗函數(shù)，代表區(qū)域中各點的權(quán)重，離P點越近，權(quán)重越高[5]，V=（u，v）T是中心點P的光流。上式的解可以由下面的方程得到：

對于在鄰域Ω內(nèi)的n個點xi，其中，

最后，方程的解為：

上式中所有的求和都是在Ω的所有點上進行的。對于Ω的選取需要經(jīng)過多次實驗找出最佳值。

1.3 雷達圖像運動矢量場的計算

如圖1所示將雷達圖像劃分ng×ng個子塊，假設(shè)圖像的大小為M×N，則每個子塊的大小為，每個子塊的頂點如圖中黑點所示。取Ω區(qū)域的大小為相鄰上下左右4個子塊的大小，如圖1中4種顏色所圈的4個不同區(qū)域，每個顏色代表一個Ω區(qū)域。

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

首先，需要計算一下前一時刻和當前時刻這兩幅相鄰時刻輸入圖像之間的平均運動矢量。然后需要計算分塊后的圖像每個網(wǎng)格點上的運動矢量。

第一，將輸入的兩幅相鄰時刻雷達圖像所有像素點都進行拉普拉斯平滑處理，平滑次數(shù)設(shè)為10次。然后對這兩幅輸入圖像進行一個簡單的閾值處理，即將灰度小于某一閾值的像素點的灰度設(shè)為零。此處閾值可以設(shè)為15 dbz。這樣就可以將一些不必要的點進行濾除，從而減少運算次數(shù)，節(jié)省運算時間。此時就完成了對輸入圖像的預處理。

第二，計算區(qū)域內(nèi)非零像素點的個數(shù)，若非零像素點的個數(shù)小于100，將區(qū)域內(nèi)的光流全部設(shè)為零。若非零點的個數(shù)大于100，則使用L-K局部光流算法計算每個區(qū)域中心點處的光流。 圖1中 4個不同區(qū)域的中心點分別為 （x1，y1），（x2，y2），（x3，y3），（x4，y4）。 其中利用使小區(qū)域上下左右各移動 5 個像素點的方法來計算x和y方向上的亮度梯度。由此可獲得除邊界子塊外其他子塊的頂點上的光流。

第三，邊界上點的值等于邊界上右下角頂點的值，比如邊界點K的值等于對角點L的值，邊界點H的值等于對角點I的值。對于缺失點的光流速度用平均光流速度來代替。如果計算出的光流的值為無效值，也將其用前面計算的平均運動矢量來代替。

最后對整幅圖像的光流速度進行1次平滑處理。然后用雙線性插值法對整幅圖像的運動矢量進行填充，假設(shè)某一子塊的大小為 m×n，4 個頂點的光流值分別為：V（x，y），V（x，y+m），V（x+n，y），V（x，n，y+m），則子塊中任一點（x′，y′）的光流為：

其中，

利用此方法既保持了圖像的規(guī)律又有利于保護邊緣，這樣即可求出圖像中每一點的光流，然后各點的光流乘以兩幅圖像的相隔時間進一步計算出圖像運動的矢量場[3]。

2 實驗結(jié)果分析

經(jīng)過多次實驗分析，最終得出以下結(jié)果：關(guān)于將圖像分塊的ng的大小，實驗過程中采用了多個不同的ng進行實驗，最終確定當ng=35時，實驗結(jié)果最為理想。本次試驗是采用的天津地區(qū)雷達站2005年6月14的數(shù)據(jù)。圖2.a和圖2.b分別為00:15和02:21時刻的雷達矢量運動場圖像。

圖2 天津雷達2005年6月14日00:15和02:21的雷達圖像矢量場的計算結(jié)果Fig.2 Tianjin radar on June 14， 2005，00:15 and 02:21calculation results of radar image vector fields

由圖2可以看出，計算出的雷達圖像矢量場非常平滑，同時光流場的方向也具有較高的準確性。

3 結(jié)束語

文中總結(jié)了在計算雷達圖像矢量場的過程中所用到的兩種主要光流算法即全局光流算法和L-K算法的主要內(nèi)容。在此基礎(chǔ)上提出了一種新的計算雷達運動矢量場的詳細算法。因為在計算過程中采用將圖像分塊的方式進行計算光流[7]，不需要對圖像上每一個像素點計算，因此，本文的算法可以有效節(jié)省運算時間。并經(jīng)過試驗表明，這種算法能很好的計算出雷達圖像[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，F(xiàn)leet 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]林存花，陳海峰.光流場模型用于非剛性醫(yī)學圖像配準[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]劉強，劉忠義，楊澤剛，等.基于可編程渲染管線的雷達圖像分層模型設(shè)計與實現(xiàn)[J].現(xiàn)代電子技術(shù)，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.