伊凡

（中國(guó)海洋大學(xué) 信息科學(xué)與工程學(xué)院，山東 青島 266100）

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

1 算 法

算法部分將主要簡(jiǎn)要介紹一下計(jì)算雷達(dá)圖像運(yùn)動(dòng)矢量場(chǎng)過(guò)程中所涉及的主要算法并對(duì)如何計(jì)算雷達(dá)圖像矢量場(chǎng)提出了一種新的詳細(xì)算法。

1.1 全局約束光流算法

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

利用泰勒展開(kāi)式展開(kāi)，最終會(huì)得到：

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

1.2 L-K算法

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

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

對(duì)于在鄰域Ω內(nèi)的n個(gè)點(diǎn)xi，其中，

最后，方程的解為：

上式中所有的求和都是在Ω的所有點(diǎn)上進(jìn)行的。對(duì)于Ω的選取需要經(jīng)過(guò)多次實(shí)驗(yàn)找出最佳值。

1.3 雷達(dá)圖像運(yùn)動(dòng)矢量場(chǎng)的計(jì)算

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

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

首先，需要計(jì)算一下前一時(shí)刻和當(dāng)前時(shí)刻這兩幅相鄰時(shí)刻輸入圖像之間的平均運(yùn)動(dòng)矢量。然后需要計(jì)算分塊后的圖像每個(gè)網(wǎng)格點(diǎn)上的運(yùn)動(dòng)矢量。

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

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

第三，邊界上點(diǎn)的值等于邊界上右下角頂點(diǎn)的值，比如邊界點(diǎn)K的值等于對(duì)角點(diǎn)L的值，邊界點(diǎn)H的值等于對(duì)角點(diǎn)I的值。對(duì)于缺失點(diǎn)的光流速度用平均光流速度來(lái)代替。如果計(jì)算出的光流的值為無(wú)效值，也將其用前面計(jì)算的平均運(yùn)動(dòng)矢量來(lái)代替。

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

其中，

利用此方法既保持了圖像的規(guī)律又有利于保護(hù)邊緣，這樣即可求出圖像中每一點(diǎn)的光流，然后各點(diǎn)的光流乘以?xún)煞鶊D像的相隔時(shí)間進(jìn)一步計(jì)算出圖像運(yùn)動(dòng)的矢量場(chǎng)[3]。

2 實(shí)驗(yàn)結(jié)果分析

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

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

由圖2可以看出，計(jì)算出的雷達(dá)圖像矢量場(chǎng)非常平滑，同時(shí)光流場(chǎng)的方向也具有較高的準(zhǔn)確性。

3 結(jié)束語(yǔ)

文中總結(jié)了在計(jì)算雷達(dá)圖像矢量場(chǎng)的過(guò)程中所用到的兩種主要光流算法即全局光流算法和L-K算法的主要內(nèi)容。在此基礎(chǔ)上提出了一種新的計(jì)算雷達(dá)運(yùn)動(dòng)矢量場(chǎng)的詳細(xì)算法。因?yàn)樵谟?jì)算過(guò)程中采用將圖像分塊的方式進(jìn)行計(jì)算光流[7]，不需要對(duì)圖像上每一個(gè)像素點(diǎn)計(jì)算，因此，本文的算法可以有效節(jié)省運(yùn)算時(shí)間。并經(jīng)過(guò)試驗(yàn)表明，這種算法能很好的計(jì)算出雷達(dá)圖像[8]的運(yùn)動(dòng)矢量場(chǎng)。

[1]楊國(guó)亮，王志良，牟世堂，等.一種改進(jìn)的光流算法[J].計(jì)算機(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á)圖像的風(fēng)暴運(yùn)動(dòng)分析 [D].青島:中國(guó)海洋大學(xué)，2009.

[7]林存花，陳海峰.光流場(chǎng)模型用于非剛性醫(yī)學(xué)圖像配準(zhǔn)[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]劉強(qiáng)，劉忠義，楊澤剛，等.基于可編程渲染管線(xiàn)的雷達(dá)圖像分層模型設(shè)計(jì)與實(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.