王興+羅永月+苗春生+張皓

摘 要：提出一種多普勒天氣雷達(dá)回波圖像追蹤預(yù)測方法。該方法以拉格朗日力學(xué)作為動力學(xué)基礎(chǔ)，利用雷達(dá)探測到的基本反射率因子作為輸入場，結(jié)合分層抽稀技術(shù)捕捉回波在空間和時間上的位移矢量，并由此預(yù)測雷達(dá)回波的位置和強(qiáng)度，進(jìn)而為短臨天氣預(yù)報業(yè)務(wù)的實施提供必要支撐。實驗表明，該方法能夠較準(zhǔn)確地預(yù)測局地天氣演變過程，1小時預(yù)測結(jié)果與實況的相關(guān)性達(dá)到70%以上，相關(guān)性的個體差異在20%以下。同時，表明該方法對于輔助預(yù)報未來1小時內(nèi)的強(qiáng)對流事件是行之有效的。

關(guān)鍵詞：拉格朗日力學(xué)； 圖像預(yù)測； 圖像追蹤； 天氣雷達(dá)； 臨近預(yù)報

DOIDOI：10.11907/rjdk.162240

中圖分類號：TP317.4

文獻(xiàn)標(biāo)識碼：A文章編號：1672-7800（2016）012-0001-04

0 引言

多普勒天氣雷達(dá)（以下簡稱雷達(dá)）是現(xiàn)代氣象業(yè)務(wù)研究及應(yīng)用中不可或缺的重要工具，基于雷達(dá)回波及圖像的追蹤與外推是短時臨近預(yù)報的關(guān)鍵性技術(shù)，也是長期以來的難點問題。準(zhǔn)確而及時的雷達(dá)外推預(yù)報，可以為局地龍卷風(fēng)、雷暴、短時強(qiáng)降水等極端災(zāi)害性天氣提供預(yù)警，以便做好防御措施，最大限度保障人民及生命財產(chǎn)安全。

近半個世紀(jì)以來，很多學(xué)者在基于雷達(dá)資料的預(yù)報領(lǐng)域作出了大量卓有成效的貢獻(xiàn)，并逐步形成了一系列較為通行的方法。例如，TREC算法通過逐區(qū)域?qū)で笙噜彆r刻雷達(dá)反射率的最大相關(guān)，跟蹤整個回波區(qū)域的移動，并且假設(shè)回波具有一致的移動方向[1]。CTREC算法則利用交叉相關(guān)分析，跟蹤反射率因子大于一定閾值區(qū)域的移動，進(jìn)而推算回波的發(fā)展[2]。TITAN是由美國國家大氣研究中心（NCAR） 研發(fā)的一套風(fēng)暴識別、跟蹤、分析和預(yù)報系統(tǒng)，其利用雷達(dá)一次完整體掃所構(gòu)成的三維結(jié)構(gòu)數(shù)據(jù)對強(qiáng)回波中心進(jìn)行識別追蹤[3-4]。SCIT算法則更加側(cè)重對雷暴單體的有效追蹤和預(yù)測[5-7]。近年來，又有一些學(xué)者試從圖形圖像學(xué)中的光流技術(shù)入手，通過分析雷達(dá)回波時序圖像中的光流場特征進(jìn)行回波強(qiáng)度和位置的外推預(yù)測[8-11]。此外，還有基于神經(jīng)網(wǎng)絡(luò)、模式匹配等多種方法[12-14]。

總體來說，這些方法的共同之處是需要分析相鄰時刻雷達(dá)圖像（或基數(shù)據(jù)）中區(qū)域的相似性，而該相似性的度量，如最大相關(guān)法，雖然可以計算出最匹配的位置，但匹配和外推結(jié)果往往表現(xiàn)出發(fā)散性或多個最優(yōu)解。并且，基于窗口平移的模板匹配算法無法適應(yīng)區(qū)域圖像的旋轉(zhuǎn)和變形等情況。

為了克服相關(guān)性度量算法所遇到的問題，本文考慮對位移場的分析加以分層，也就是每個特征運(yùn)動被認(rèn)為是確定在相對粗糙的空間分辨率下的平穩(wěn)變化趨勢的總和，然后衍生為更高的空間分辨率下的小幅度局部修正，并對此過程進(jìn)行多次迭代。對于位移場的分析，將充分考慮雷達(dá)回波所指示風(fēng)暴的動力學(xué)特征，運(yùn)用拉格朗日力學(xué)相關(guān)理論構(gòu)建預(yù)測模型。該方法的研究意義在于進(jìn)一步提高中小尺度、強(qiáng)對流天氣事件的預(yù)測能力，且相對于光流等一些大運(yùn)算量算法，本方法運(yùn)算規(guī)模更小，進(jìn)而能夠更好地滿足短臨預(yù)報業(yè)務(wù)高時效性的要求。

1 理論與方法基礎(chǔ)

1.1 雷達(dá)回波外推預(yù)測

大量研究表明，合理的外推預(yù)測算法可以為降水、雷暴、冰雹等對流天氣的預(yù)報提供重要支撐[15-17]?；谕馔祁A(yù)報的一般性描述為：

其中，pt（x，y）表示任一位置的回波強(qiáng)度，U和V分別表示回波在水平和垂直方向上的偏移量，由U和V共同組成回波移動的速度矢量。g表示一個函數(shù)，用來計算單位時間間隔后回波強(qiáng)度的變化。根據(jù)式（1），Δpt（x，y）反映了任一點（x， y） 在t 時刻回波強(qiáng)度的變化情況，U和V反映了回波移動的方向和速度。此外，函數(shù)g代表一個拉格朗日動力學(xué)過程，在此過程中雷達(dá)回波的強(qiáng)度是由其沿回波路徑移動時在拉格朗日坐標(biāo)系統(tǒng)中的歷史變化推導(dǎo)出的，也就是用當(dāng)前回波演變的趨勢預(yù)測回波未來的位置和強(qiáng)度。究其趨勢預(yù)測的方法，多年來諸多學(xué)者作出很多研究，本文主要從拉格朗日力學(xué)角度進(jìn)行分析，提出一種雷達(dá)回波圖像追蹤預(yù)測的方法。

1.2 拉格朗日力學(xué)

拉格朗日力學(xué)是由Joseph Lagrange[18-19]最早提出的一種力學(xué)分析方法。由于該方法引用了廣義坐標(biāo)的概念，使得對力學(xué)相關(guān)問題的研究更具普適性。

在雷達(dá)回波圖像預(yù)測研究中，如何準(zhǔn)確得到回波運(yùn)動矢量是預(yù)測需要解決的關(guān)鍵問題。在不考慮天氣系統(tǒng)的非線性變化時，拉格朗日力學(xué)模型能夠滿足構(gòu)建回波發(fā)展演變過程的算法要求，式（1） 可以改寫為：

有研究表明，對整個回波圖像采用統(tǒng)一的U和V所構(gòu)成的位移矢量，可適用于對大尺度天氣系統(tǒng)的預(yù)測分析，如對衛(wèi)星圖像的外推預(yù)測和云導(dǎo)風(fēng)的分析[20]。但對于局地強(qiáng)對流天氣系統(tǒng)，預(yù)測結(jié)果往往與實際偏差較大。因此，這也是本文將重點闡述解決的問題。

2 雷達(dá)圖像追蹤預(yù)測

2.1 基于拉格朗日力學(xué)的追蹤算法

在上述理論基礎(chǔ)上，根據(jù)大氣運(yùn)動演變發(fā)展的規(guī)律及其在雷達(dá)回波圖像上的表征特點，構(gòu)建基于拉格朗日力學(xué)的追蹤算法模型，如式（3） 所示。

該模型假定所預(yù)測的回波圖像是當(dāng)前和過去若干個回波圖像以固定時間間隔而變化的函數(shù)。f2表征一個用于估測單位時間間隔前后回波圖像各相應(yīng)網(wǎng)格點回波強(qiáng)度值變化率的函數(shù)，即回波的演變趨勢。在不斷生消、發(fā)展的對流系統(tǒng)中，ΔP可以為正，也可以為負(fù)數(shù)。式（4）～式（6） 進(jìn)一步表明了函數(shù)f1和f2，即回波移動矢量及回波強(qiáng)度的計算方法。

2.2 中心極值濾波

為了降低雷達(dá)雜波對位移矢量計算的不良影響，本節(jié)提出采用一種濾波器對雷達(dá)基數(shù)據(jù)進(jìn)行濾波處理。其基本思想是：逐網(wǎng)格分析回波強(qiáng)度特征，如果某格點的值大于周邊最相鄰一圈（共8個網(wǎng)格）的最大值，或者該值小于周邊最相鄰一圈的最小值，則將當(dāng)前網(wǎng)格點回波強(qiáng)度值用上述8個網(wǎng)格的最大值或最小值替代。

如圖1所示，位于當(dāng)前中心點的數(shù)值39大于其最鄰近一圈8個網(wǎng)格的最大值。因此，使用數(shù)值23替換當(dāng)前網(wǎng)格的39。

從圖像上看，該濾波方法可顯著降低圖像中的椒鹽噪聲，從實際效果上看，該方法可以有效過濾單點的雜波奇異值，較傳統(tǒng)的均值濾波和中值濾波更好地保留了回波細(xì)節(jié)[21-22]，特別是回波中梯度變化較大的邊緣區(qū)域。

2.3 分層外推預(yù)測算法

由于天氣系統(tǒng)復(fù)雜多變，特別是尺度較小的局地強(qiáng)對流，其生命周期短的只有幾分鐘到幾十分鐘，由于其空間尺度小，生消速度快，因此，包括基于拉格朗日力學(xué)在內(nèi)的各種線性關(guān)系外推算法，其預(yù)測準(zhǔn)確性都存在一定的局限性。為了改善這一問題，本節(jié)提出采用分層的外推預(yù)測方法。該思想最早由Bellerby等 [20]提出，并研究應(yīng)用于衛(wèi)星圖像的云頂平流場分析中。

該算法的關(guān)鍵流程為：先將當(dāng)前雷達(dá)回波圖像逐級抽稀，降低圖像的分辨率，以模糊回波細(xì)節(jié)，由此粗略估算出回波主體的移動趨勢；然后再反向逐級提高圖像分辨率，在較粗的移動趨勢基礎(chǔ)上，細(xì)化和訂正位移矢量的細(xì)節(jié)。從而產(chǎn)生一個在空間上連續(xù)和平滑的且不受模板邊界不連續(xù)性影響的矢量場。計算方法如式（7）所示：

對于每一級抽稀計算，都是將當(dāng)前一級各網(wǎng)格點的回波強(qiáng)度值經(jīng)由公式（7）計算，并往復(fù)迭代。其中P表示某一點的回波強(qiáng)度，L和L-1代表抽稀的層級，在本文下述實驗中，采用的最高層級為4。

在計算兩個相鄰時刻圖像中回波的位移時，可以借鑒交叉相關(guān)法，計算方法如式（8），在每個選定的匹配窗口遍歷出最大相關(guān)矩陣的位置，從而輸出位移矢量。

式（8） 中，P和P 分別表示相鄰兩個時刻（如t-Δt與t）的回波，（x， y） 表示圖像中的某一點，X和Y表示匹配窗口的大小。再將兩幅回波圖像之間的網(wǎng)格還原或內(nèi)插到其先前空間分辨率的兩倍，重復(fù)上述匹配。該迭代過程還考慮到了由非矩形網(wǎng)格代表的局部扭曲，結(jié)合這些局部扭曲，使外推預(yù)測算法能夠適應(yīng)旋轉(zhuǎn)、擴(kuò)展、縮小等回波圖像形態(tài)上的變化。如此插值和匹配計算，迭代直到網(wǎng)格分辨率達(dá)到原始雷達(dá)圖像分辨率。

3 實驗與結(jié)果分析

3.1 實驗數(shù)據(jù)說明

為檢驗所述方法的預(yù)測效果，本實驗數(shù)據(jù)使用2016年6月南京地區(qū)多普勒天氣雷達(dá)的基數(shù)據(jù)文件。該雷達(dá)使用VCP-21體掃模式，探測周期為6分鐘。實驗選用1.5°和2.4°仰角的基本折射率數(shù)值。為方便計算，實驗前將原始數(shù)據(jù)由極坐標(biāo)系統(tǒng)轉(zhuǎn)換為平面直接坐標(biāo)系統(tǒng)，數(shù)據(jù)的圖像分辨率為920×920。為減少樣本數(shù)量，從全部7199個基數(shù)據(jù)文件中篩選出以230庫長為半徑，其覆蓋區(qū)域內(nèi)具有大面積強(qiáng)回波的數(shù)據(jù)文件，共計880個。

3.2 實驗結(jié)果分析

為檢驗雷達(dá)回波圖像預(yù)測的準(zhǔn)確性，使用與預(yù)測同一時刻的雷達(dá)實際探測數(shù)據(jù)作比對分析，計算過程采用交叉相關(guān)檢驗法。

首先以自然日為單位，統(tǒng)計逐日樣本數(shù)據(jù)中每份預(yù)測結(jié)果與實況交叉檢驗的相關(guān)系數(shù)的平均值，如圖2所示。

圖2中3種圖案標(biāo)記分別表示預(yù)測6分鐘、30分鐘和60分鐘的檢驗結(jié)果，每個值代表當(dāng)日所有樣本檢驗結(jié)果的平均值。橫坐標(biāo)為2016年6月的逐個日期，縱坐標(biāo)為相關(guān)系數(shù)，其中橫坐標(biāo)4、5、9、10等日期沒有標(biāo)記圖案，原因是這些日期的當(dāng)日為晴天或少云，體現(xiàn)在雷達(dá)上沒有強(qiáng)的大面積回波，因此沒有列入樣本數(shù)據(jù)進(jìn)行分析。從圖2中還可以看出，本方法預(yù)測未來6分鐘的結(jié)果與實況相比，相關(guān)系數(shù)超過87%，平均達(dá)到93%以上，而隨著預(yù)測時效的延長，預(yù)測準(zhǔn)確率逐步下降，在未來60分鐘的預(yù)測中，全月平均相關(guān)系數(shù)為70%左右。

進(jìn)一步統(tǒng)計分析每次預(yù)測準(zhǔn)確率的穩(wěn)定性。以6月19日全天樣本數(shù)據(jù)為例，統(tǒng)計每批樣本所預(yù)測6、12、18至60分鐘結(jié)果分別與實況交叉相關(guān)檢驗的情況，如圖3所示。

圖3中，每個柱狀條的頂端和底端分別表示檢驗的相關(guān)系數(shù)的最大值和最小值，柱狀條中間的黑色方形表示相關(guān)系數(shù)的均值。可以看出，隨著預(yù)測時效的增長，其預(yù)測準(zhǔn)確率的個體差異也隨之增大。在前6分鐘的預(yù)測中，該差異約為3%，30分鐘時約為8%，而到預(yù)測60分鐘時，差異進(jìn)一步增大到20%。結(jié)果與強(qiáng)對流天氣系統(tǒng)具有生命史短、突發(fā)性強(qiáng)，水氣生消發(fā)展變化快的特點是相一致的。

4 結(jié)語

由于天氣系統(tǒng)復(fù)雜多變，特別是對于中小尺度的對流系統(tǒng)，其生消、發(fā)展時間短、變化快，如何進(jìn)行準(zhǔn)確、有效的預(yù)報是提升當(dāng)今短時臨近預(yù)報的關(guān)鍵環(huán)節(jié)之一?？紤]到大氣中水氣等物質(zhì)的移動變化應(yīng)遵循一般力學(xué)規(guī)律，而拉格朗日力學(xué)正是表征和計算動力學(xué)問題的普適性方法，因此，本文的預(yù)測動力模型建立在拉格朗日力學(xué)關(guān)系基礎(chǔ)之上。又由于天氣系統(tǒng)的變化表現(xiàn)在雷達(dá)圖像上，其回波圖形具有相當(dāng)?shù)牟淮_定性，因此，本文提出采用分層的位移場分析方法，先假定位移矢量是在相對粗糙的空間分辨率下的平穩(wěn)變化趨勢的總和，然后在更高空間分辨率下作小幅度局部修正，并如此進(jìn)行多次迭代。為了減少低仰角雜波對實驗結(jié)果的影響，提出采用中心極值濾波對實驗數(shù)據(jù)進(jìn)行處理。結(jié)合上述理論構(gòu)建起基于拉格朗日力學(xué)的追蹤預(yù)測模型及算法流程，以雷達(dá)基本反射率因子作為輸入場，追蹤和預(yù)測回波在空間和時間上的位移矢量，并由此預(yù)測未來一段時間雷達(dá)回波的位置和強(qiáng)度。

實驗部分采用1個月樣本數(shù)據(jù)對本算法模型進(jìn)行檢驗，通過預(yù)測結(jié)果與同時刻實況的比對分析，得出兩者的相關(guān)性和個例穩(wěn)定性等評價指標(biāo)。結(jié)果表明，該方法能夠較好地預(yù)測局地天氣系統(tǒng)的演變過程，在未來30分鐘的預(yù)測中準(zhǔn)確率平均超過80%，且對于輔助預(yù)報未來60分鐘內(nèi)的局地龍卷風(fēng)、強(qiáng)降水、雷暴等災(zāi)害性事件具有實踐應(yīng)用的價值。

參考文獻(xiàn)：

[1] 劉紅艷，魏鳴.多普勒雷達(dá)風(fēng)場資料在臨近預(yù)報中的應(yīng)用[J].大氣科學(xué)學(xué)報，2015（4）：483-491.

[2] 鄭永光，林隱靜，朱文劍，等.強(qiáng)對流天氣綜合監(jiān)測業(yè)務(wù)系統(tǒng)建設(shè)[J].氣象，2013（2）：234-240.

[3] 周康輝，鄭永光，藍(lán)渝.基于閃電數(shù)據(jù)的雷暴識別、追蹤與外推方法[J].應(yīng)用氣象學(xué)報，2016，（2）：173-181.

[4] 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.

[5] SHAH S， NOTARPIETRO R， BRANCA M.Storm identification，tracking and forecasting using high-resolution images of short-range X-band radar[J].Atmosphere， 2015， 6（5）：579-606.

[6] 莊旭東，胡勝，陳榮，等.“雨燕”中風(fēng)暴算法與新一代雷達(dá)SCIT產(chǎn)品的對比分析[J].熱帶氣象學(xué)報，2011（3）：299-306.

[7] DUAN Y， XU Y， ZHI S.Application analysis of the hail suppression operation based on the improved SCIT Algorithm[J].Meteorology & Disaster Reduction Research， 2014（15）：23-29.

[8] 曹春燕，陳元昭，劉東華，等.光流法及其在臨近預(yù)報中的應(yīng)用[J].氣象學(xué)報，2015（3）：471-480.

[9] GARCIA F， CERRI P， BROGGI A， et al.Data fusion for overtaking vehicle detection based on radar and optical flow[J].2012， 7（2272）：494-499.

[10] 王興，王新，苗春生，等.基于GPU加速的雷暴追蹤外推方法研究[J].南京師范大學(xué)學(xué)報：工程技術(shù)版，2015（1）：35-42.

[11] STAINVAS OLSHANSKY I， BILIK I， BIALER O.Doppler-Based Segmentation and Optical Flow in Radar Images： US20160084953[P].2016.

[12] 盛仲飆.BP神經(jīng)網(wǎng)絡(luò)在數(shù)據(jù)預(yù)測中的應(yīng)用[J].軟件導(dǎo)刊，2016（1）：147-148.

[13] 王利卿，黃松杰.基于多尺度卷積神經(jīng)網(wǎng)絡(luò)的圖像檢索算法[J].軟件導(dǎo)刊，2016（2）：38-40.

[14] WANG X， GU Y H， MIAO C S， et al.Parallelization and performance optimization of radar extrapolation algorithm with OpenCL[J].Journal of Internet Technology， 2016（17）：323-330.

[15] 王丹.雷達(dá)外推預(yù)報與暴雨數(shù)值模式融合預(yù)報降水方法研究[D].北京：中國氣象科學(xué)研究院，2013.

[16] 張蕾.多普勒雷達(dá)回波演變的動力學(xué)分析及臨近預(yù)報算法改進(jìn)[D].南京：南京信息工程大學(xué)，2015.

[17] FOX N I， WEBB R， BALLY J， et al.The impact of advanced nowcasting systems on severe weather warning during the sydney 2000 forecast demonstration project：3 November 2000[J].Weather & Forecasting， 2004， 19（1）：97-114.

[18] 李艷艷.相似空間中不變的歐拉-拉格朗日方程[J].河南大學(xué)學(xué)報：自然科學(xué)版，2014，03：273-276.

[19] ERICKSEN R E， GUITERAS J J， LARRIVEE J A， et al.A parachute recovery system dynamic analysis [J].Journal of Spacecraft & Rockets， 1967， 4（3）：321-326.

[20] BELLERBY T J.High-resolution 2-D cloud-top advection from geostationary satellite imagery[J].IEEE Transactions on Geoscience & Remote Sensing， 2006， 44（12）：3639-3648.

[21] 莫曉齊，何愛.基于自適應(yīng)開關(guān)中值濾波算法的工程圖像處理[J].軟件導(dǎo)刊，2014（3）：55-59.

[22] 李佐勇，湯可宗，胡錦美，等.椒鹽圖像的方向加權(quán)均值濾波算法[J].中國圖象圖形學(xué)報，2013（11）：1407-1415.

（責(zé)任編輯：陳福時）