基于改進光流法的雷達圖像運動估計

王志斌，肖艷姣，吳 濤

(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方法)?；诠饬鞯幕炯s束方程，還產生了其他光流方法，但大多數可以用變分模型來描述，即最小化某個能量泛函的過程。這個泛函可以由兩項組成，即平滑項和數據項。數據項和圖像本身的數據有關，如雷達回波數據、圖像的灰度等，數據項決定了運動的真實程度。平滑項是泛函的附加項，進行條件約束，主要解決方程的不適定性。確定了泛函方程后，可以運用歐拉-拉格朗日方程進行求解，其數值解一般采用松弛迭代或最優下降算法求出。光流法在計算機圖像處理、模式識別、交通、醫學等領域[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]提出的算法進行。首先針對不同雷達型號對數據進行規整和缺測填補，對原始數據進行有效性檢查，對孤立回波進行過濾。然后選取了五個回波特征量作為候選因子，對超折射回波進行抑制，計算出超折射回波和正?；夭ǖ母怕史植紙D，找出兩者的差別，建立各因子的隸屬函數。最后通過模糊邏輯算法，把這些因子進行綜合運算，得到控制值，大于某值則去掉該回波，并用周圍正?；夭ㄟM行替代，否則該值正常不進行處理。

進行單站質量控制后，把極坐標的單站雷達資料轉換成三維的直角坐標格式，采用文獻[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