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

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

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

DOIDOI：10.11907/rjdk.162240

中圖分類號(hào)：TP317.4

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

0 引言

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

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

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

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

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

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

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

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

1.2 拉格朗日力學(xué)

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

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

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

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

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

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

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

2.2 中心極值濾波

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

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

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

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

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

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

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

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

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

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

3.1 實(shí)驗(yàn)數(shù)據(jù)說明

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

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

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

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

圖2中3種圖案標(biāo)記分別表示預(yù)測(cè)6分鐘、30分鐘和60分鐘的檢驗(yàn)結(jié)果，每個(gè)值代表當(dāng)日所有樣本檢驗(yàn)結(jié)果的平均值。橫坐標(biāo)為2016年6月的逐個(gè)日期，縱坐標(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ù)測(cè)未來6分鐘的結(jié)果與實(shí)況相比，相關(guān)系數(shù)超過87%，平均達(dá)到93%以上，而隨著預(yù)測(cè)時(shí)效的延長(zhǎng)，預(yù)測(cè)準(zhǔn)確率逐步下降，在未來60分鐘的預(yù)測(cè)中，全月平均相關(guān)系數(shù)為70%左右。

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

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

4 結(jié)語

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

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

參考文獻(xiàn)：

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

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

[3] 周康輝，鄭永光，藍(lán)渝.基于閃電數(shù)據(jù)的雷暴識(shí)別、追蹤與外推方法[J].應(yīng)用氣象學(xué)報(bào)，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)品的對(duì)比分析[J].熱帶氣象學(xué)報(bào)，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ù)報(bào)中的應(yīng)用[J].氣象學(xué)報(bào)，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é)報(bào)：工程技術(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ù)測(cè)中的應(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ù)報(bào)與暴雨數(shù)值模式融合預(yù)報(bào)降水方法研究[D].北京：中國氣象科學(xué)研究院，2013.

[16] 張蕾.多普勒雷達(dá)回波演變的動(dòng)力學(xué)分析及臨近預(yù)報(bào)算法改進(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é)報(bào)：自然科學(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é)報(bào)，2013（11）：1407-1415.

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