施 恩,李 騫,顧大權,趙章明
(國防科學技術大學 氣象海洋學院,南京 211101)
臨近預報主要指0~3h的高時空分辨率的天氣預報,主要預報對象包括強降水、大風、冰雹等災害性天氣[1]。目前,臨近預報的主要手段是基于天氣雷達資料的雷達回波外推技術,即根據當前時刻雷達觀測結果,推測雷達回波未來的位置和強度,以實現(xiàn)對強對流系統(tǒng)的跟蹤預報[2-3]。目前常用的雷達回波外推方法是質心跟蹤法和交叉相關法(Tracking Radar Echoes by Correlation, TREC)[4]。質心跟蹤法通過預測雷達回波的質心來推測下一時刻的回波位置,該方法對采集到的信息的利用更加充分,但它依賴閾值來識別風暴單體,僅適用于對風暴的追蹤,難以應用于預測大范圍降水的回波變化[5]; TREC通過對雷達回波劃分區(qū)域,并求得各個區(qū)域當前時刻與前一時刻的相關系數,相關系數最大的區(qū)域即是回波移動矢量的終點[6]。研究人員在TREC算法的基礎上,進一步發(fā)展了COTREC(Continuity Of TREC vectors)[7]和DITREC(Difference Image-based TREC)[8]等方法,但是此類方法僅根據幾個時刻的回波特征推測下一時刻的回波分布,數據利用率較低,并且假設回波是線性演變的,而實際情況下回波的變化更為復雜。
針對上述傳統(tǒng)雷達回波外推方法存在的問題,本文在傳統(tǒng)卷積神經網絡(Convolutional Neural Network, CNN)的基礎上,提出了一種基于輸入的動態(tài)卷積神經網絡結構(Dynamic Convolutional Neural Networks based on Input,DCNN-I),該網絡從歷史的時序雷達回波圖像集中學習回波變化規(guī)律,從而實現(xiàn)對雷達回波圖像的預測與外推。DCNN-I通過改進傳統(tǒng)CNN的結構和網絡訓練過程中的參數,使網絡中直接影響圖像變換的卷積核在網絡測試階段仍然能夠隨著輸入的不同而變化,更加適用于雷達回波外推這一類輸入圖像與輸出圖像之間存在較強相關性的問題,能夠取得較好的雷達回波圖像預測結果。
本文首先介紹DCNN-I的結構和特點,進而介紹利用該網絡進行雷達回波外推的全過程。將大量有序的雷達回波強度等高平面顯示 (Constant Altitude Plan Position Indicator, CAPPI) 圖像集構造成訓練樣本集和測試樣本集,利用訓練樣本集訓練網絡,使網絡收斂,并用測試樣本集測試網絡,得到雷達回波外推圖像。在對比實驗中采用降水預報業(yè)務中常用的COTREC算法和DITREC算法作為對比方法,并從預測圖像的準確率和外推時效兩方面進行對比。
卷積神經網絡是深度學習的重要分支,由于CNN具有特殊的網絡結構,在圖像處理領域具有深遠的研究前景,近年來得到高度重視并引起廣泛研究[9-12]。CNN是一種具有多層結構的網絡,其低層由卷積層和下采樣層交替組成,是網絡提取圖像特征的重要環(huán)節(jié),高層為分類器,一般為全連接層[12]。
考慮到雷達回波外推問題的輸入與輸出之間存在較強相關性,本文在傳統(tǒng)卷積神經網絡基礎上提出一種基于輸入的動態(tài)卷積神經網絡結構(DCNN-I),增加了動態(tài)子網絡(Dynamic Sub-Network, DSN)和概率預測層(Probability Prediction Layer, PPL),PPL中的卷積核由DSN計算得出,訓練階段結束之后卷積核可根據輸入圖像序列的不同而動態(tài)變化。DCNN-I的整體結構如圖1所示。
圖1 DCNN-I的結構 Fig. 1 Architecture of DCNN-I
其中DSN與傳統(tǒng)的CNN具有相同的結構,用于計算概率向量,PPL中將輸入圖像序列中的最后一幅雷達回波圖像與卷積核進行兩次卷積操作得到最終的外推圖像。DSN中包含DCNN-I的全部權值參數,其輸出作為PPL的卷積核。
在傳統(tǒng)的卷積神經網絡中,卷積核在網絡訓練階段通過反向傳播算法不斷得到更新,而在測試階段保持不變。與傳統(tǒng)的CNN不同,DCNN-I中PPL的卷積核是DSN對輸入圖像序列進行處理的結果,因此在網絡測試階段仍會隨著輸入樣本的不同而變化,使網絡具有動態(tài)特性,對于雷達回波外推這一類輸出圖像與輸入圖像強相關的問題具有更好的適用性。
在DCNN-I中,輸入的雷達回波圖像首先輸入動態(tài)子網絡DSN,經過卷積、采樣以及激活函數處理后,得到一維列向量VPV和一維行向量HPV。DSN可以視為一個獨立的卷積神經網絡,其低層包含5個卷積層和4個下采樣層,高層為1個全連接層,DSN的結構如圖2所示。
DSN中包含C1、C2、C3、C4和C5共5個卷積層,圖2中給出了各個卷積層包含的卷積核個數(在數值上與輸出特征圖數量相同)以及輸出特征圖的分辨率,其中C5層的卷積核大小為7×7,其余各個卷積層的卷積核大小均為9×9。在卷積層中,輸入特征圖與卷積核相卷積,卷積結果加上一個偏置項后作為激活函數的輸入,經過激活函數處理后得到該層的輸出特征圖。DSN中的卷積層的計算公式為:
(1)
圖2 DSN的結構 Fig. 2 Architecture of DSN
DSN中包含S1、S2、S3和S4共4個下采樣層,圖2中給出了各個下采樣層輸出特征圖的數量和分辨率,各層中采樣核的大小為2×2,采樣步長為2。在下采樣層中,對輸入特征圖降分辨率,在提取圖像特征的同時降低計算復雜度,采樣時以步長為2進行連續(xù)采樣,采樣過后輸出特征圖的分辨率降為原輸入特征圖的1/4。DSN中的下采樣層的計算公式為:
(2)
在DSN的全連接層F1中,首先將C5層的32個分辨率為16×16的輸出特征圖按列順序展開,得到大小為521×1的列向量。分別計算該向量與垂直參數矩陣WV和水平參數矩陣WH的外積,計算結果分別加上一個偏置向量,再經過Softmax函數處理之后,得到一維列向量VPV和一維行向量HPV。DSN中的F1層的計算公式為:
(3)
其中:aC5表示C5層輸出的大小為512×1的特征向量,WV和WH均為41×512的參數矩陣,BV和BH均為41×1的偏置向量,計算VPV時需要將計算結果轉置,最終得到大小為1×41的VPV和41×1的HPV。由于Softmax函數的特性,得到的一維向量VPV和HPV中所有元素均為正值且元素之和為1,因此可以將VPV視為垂直概率向量,將HPV視為水平概率向量。
PPL包含兩個輸入,輸入圖像序列中的最后一幅圖像作為PPL的輸入特征圖,DSN的輸出概率向量VPV和HPV作為PPL的卷積核,在PPL中,將對其輸入特征圖進行兩次卷積操作,得到外推圖像,PPL的結構如圖3所示。
PPL包含DC1和DC2兩層結構,垂直概率向量VPV為DC1層的卷積核,水平概率向量HPV為DC2層的卷積核。在PPL中進行兩次卷積操作,在DC1層中,將輸入圖像序列中的最后一幅雷達回波圖像與VPV相卷積,得到分辨率為240×280的輸出特征圖;在DC2層中,將DC1層的輸出特征圖與HPV相卷積,得到分辨率為240×240的預測圖像。由于概率向量的特殊性(元素均為正,元素和為1),可以將DC1層和DC2層中的卷積操作視為垂直方向的預測和水平方向的預測。
圖3 PPL的結構 Fig. 3 Architecture of PPL
實驗中采用南京、杭州、廈門三市的2016年7至9月份CINRAD-SA型多普勒天氣雷達資料訓練DCNN-I,數據源中大約包含60 000個樣本,可滿足神經網絡訓練需求。由于雷達系統(tǒng)中為用戶提供可視化界面的PUP軟件沒有批量生成雷達回波圖像的功能,因此本文需要研究雷達回波強度CAPPI圖像生成方法,并構造數據集。
CINRAD-SA型多普勒天氣雷達基數據為三維極坐標下的散亂數據,需要依次經過坐標轉換、數據插值、水平采樣和繪制才能夠得到一幅回波強度CAPPI圖像, 圖4為多普勒天氣雷達基數據預處理過程。
圖4 多普勒天氣雷達基數據預處理過程 Fig. 4 Preprocessing process of Doppler weather radar data
首先根據文獻[13]中提供的方法,將三維極坐標下的數據轉換到三維笛卡爾直角坐標系中,采用反距離加權法[14]進行數據插值,得到三維笛卡爾直角坐標系下的規(guī)整網格數據。然后對數據進行水平采樣,提取某一高度下的二維平面數據,將數據映射至0~255,便能夠得到回波強度CAPPI灰度圖像,每幅圖像的分辨率為2 000×2 000。為了降低CNN的訓練復雜度,不對圖像進行顏色映射。
雷達基數據經過數據預處理之后得到回波強度CAPPI灰度圖像,此時圖像的邊緣部分包含大量的空白信息,需要進行裁剪和壓縮處理,得到分辨率為280×280的灰度圖。然后再根據圖像之間的時間順序進行采樣,構造網絡訓練和測試需要的樣本,得到的每一組樣本包含5幅圖像{x1,x2,x3,x4,y},每幅圖像之間的時間間隔恒定(6 min)。其中{x1,x2,x3,x4}是分辨率為280×280的輸入樣本圖像序列,{y}是再次經過裁剪、分辨率為240×240的對照標簽(網絡的期望輸出),此時對第5幅圖像進行裁剪的目的是保證對照標簽與DCNN-I的輸出圖像的分辨率相一致。最終得到包含48 000幅樣本圖像的訓練樣本集和包含4 800幅樣本圖像的測試樣本集,圖5為數據集中的一組樣本圖像。
圖5 數據集中的一組樣本圖像 Fig. 5 A group of sample images in the dataset
為了驗證本文方法的有效性,在2.40 GHz CPU、內存大小為4 GB的PC上實現(xiàn)了算法。其中卷積核和參數矩陣(VPV和HPV)的初始化依據Xavier初始化方法[15],將偏置初始化為0向量。訓練過程采用反向傳播算法計算誤差,并采取帶動量的梯度下降法更新網絡參數,學習率為0.000 1,動量系數為0.5。網絡訓練階段采取批訓練的方式,每次輸入10組樣本,網絡迭代次數為40。
對比實驗中,對2016年7月7日南京地區(qū)、2016年7月16日杭州地區(qū)以及2016年8月11日廈門地區(qū)三個大范圍降水過程進行回波強度的外推,并分別從回波圖像和預測指標兩方面分析不同雷達回波外推方法的優(yōu)劣。
對比實驗中分別對3個降水過程進行雷達回波強度外推,結果如圖6所示,圖中的第1行為輸入的雷達回波圖像,第2行為實際觀測到的雷達回波圖像,第3行和第4行分別為采用DCNN-I和COTREC的外推結果,每一行中相鄰的圖像之間的時間間隔為6 min。
由圖6可知,與COTREC算法相比,基于DCNN-I的方法外推圖像與觀測到的實況圖像相似程度更高;尤其是在圖像邊緣部分,采用COTREC算法很難準確計算出觀測區(qū)域邊緣部分的矢量場,而DCNN-I可以在訓練過程中學習觀測區(qū)域邊緣部分的回波變化特征,從而給出較為可靠的預測。此外,圖6中采用DCNN-I方法得到的外推圖像都較為模糊,這是因為雷達回波外推本身存在的不確定性,其輸入和輸出并不是簡單的對應關系,DCNN-I在訓練過程中學習了各種天氣過程回波變化規(guī)律,網絡在預測時幾乎不可能對整個觀測區(qū)域內雷達回波作出精準的預測,導致外推的回波圖像較為模糊。
在3.1節(jié)中,本文從圖像層面分析外推的雷達回波圖與觀測到的雷達回波實況圖的相似性,對比分析了外推結果。除了根據圖像的相似性來判斷外推準確性以外,還可以通過預測指標來評價外推結果。本文中采用的預測指標包括: 臨界成功指數(Critical Success Index,CSI)、誤報率(False Alarm Rate,F(xiàn)AR)以及探測概率(Probability Of Detection,POD)[16]。
圖6 三個降水過程的預測結果 Fig. 6 Prediction results of the three example
首先設定一個回波強度閾值(對于回波圖像,即代表某一灰度閾值),若雷達回波圖像中某個像素的灰度值超過該閾值,則判定該像素對應區(qū)域的預報值是活躍的,反之,則判定該像素的預報值為不活躍。因此對于單個像素的預測分為成功(預報值與實際值均活躍,標記為S)、漏報(預報值不活躍,實際值活躍,標記為M)和空報(預報值活躍,實際值不活躍,標記為F)。本文以nS、nM和nF分別表示預測圖像中成功、漏報和空報的格點數,則定義CSI、FAR和POD的計算公式如下:
表1 不同雷達回波外推方法的結果對比Tab. 1 Comparison of mean and variance for different method over 15 prediction steps
(4)
CSI、FAR和POD能夠較好地評估預測的雷達回波圖像的預報效果,但是仍然不能直觀地反映預測的雷達回波的準確率,因此本文進一步采用均方誤差(Mean Square Error, MSE)來衡量雷達回波外推的準確率,其計算形式為:
(5)
進行計算之前,需要先根據Z-R關系[17],將圖像中的像素每個像素的值換算成每個格點的降水強度。上式中,N為觀測區(qū)域內Ω的格點總數,實際上就是指圖像包含的像素點的總數,F(xiàn)為實際觀測的雷達回波圖像,F(xiàn)′為預測的雷達回波圖像,t0為預報時刻,τ為預測的時效,F(xiàn)′(t0+τ,x)表示在t0+τ時刻,預測的雷達回波在x格點處的降水強度。
實驗中,采用COTREC算法和DITREC算法作為對比方法,此外還將輸入圖像中的最后一幅( Last Image)作為對比的依據。對2016年8月11日廈門地區(qū)大范圍降水過程進行連續(xù)15次外推預測,并分別計算4種方法預測結果的MSE、CSI、FAR和POD預測指標,對比結果如圖7所示。
實驗中回波強度閾值設置為10 dBZ,對應圖像的灰度值為25。由圖7可知,隨著預測時間的增長,實況回波圖像與Last Image之間的差別越來越大,COTREC算法與DITREC算法的預測結果較為接近,基于DCNN-I的雷達回波外推方法預測結果最好,尤其是在準確率(以MSE表示)和誤報率方面,采用DCNN-I優(yōu)化效果明顯。該降水過程的15次外推過程中,各項預測指標的平均值和方差如表1所示。
從表1可得,采用DCNN-I的雷達回波外推方法能夠得到更好的預測結果,預測結果也較為穩(wěn)定,預測準確率約有10%的提升。
本文通過去相關時間L來描述不同雷達回波外推方法的最大預測時效[18]。去相關時間L與預測的雷達回波圖像與觀測到的實況回波圖像之間的相關系數c的值有關,相關系數c的計算公式為:
(6)
其中,G與G′分別表示實況雷達回波圖像和預測的雷達回波圖像中某個像素點的灰度值。相關系數的取值為c∈[0,1],c的值越大,表明兩幅圖像越相似,雷達回波外推的效果越好。由于相關系數c符合指數規(guī)律,因此將相關系數下降到1/e時對應的時間定義為去相關時間L。當相關系數c低于1/e時,可以判定預測的雷達回波不再具有有效性,因此去相關時間L能夠作為衡量雷達回波外推時效的依據。圖8為分別采取COTREC、DITREC和DCNN-I方法對三個降水過程進行外推時,相關系數關于預測時間的變化。
圖7 不同方法外推結果的預測指標對比 Fig. 7 Comparison of four precipitation nowcasting metrics on different method
圖8 三次降水過程中相關系數的變化 Fig. 8 Correlation coefficient changes over forecast time during three precipitation processes
由圖8可知:在三次降水過程中,相關系數隨著預測時間呈指數遞減;與COTREC算法和DITREC算法相比,采用DCNN-I的外推方法對應的相關系數隨時間下降更為緩慢。從去相關時間L來看,采用傳統(tǒng)方法進行雷達回波外推得到的L值約為2.4~3 h,采用DCNN-I的雷達外推方法得到的L則大于4 h,證明基于DCNN-I的雷達回波外推方法能夠有效延長外推時效。
卷積神經網絡在二維圖像的識別和預測上具有深遠的應用前景,目前卷積神經網絡在天氣預報領域的應用較少,本文嘗試利用卷積神經網絡進行雷達回波外推,在傳統(tǒng)卷積神經網絡的基礎之上,通過增加DSN和PPL,提出了基于輸入的動態(tài)卷積神經網絡DCNN-I,該網絡中PPL的卷積核由DSN計算得出,卷積核中包含當前輸入圖像的特征,在網絡測試階段仍然能夠根據輸入圖像的不同而變化,從而保證網絡輸入圖像和輸出圖像之間的強相關性。與傳統(tǒng)的雷達回波外推方法相比,基于DCNN-I的雷達回波外推方法通過網絡訓練從大量的雷達數據中學習雷達回波變化特征,有效提高了數據的利用率,對于“陌生”的雷達回波具有較強的適應性。實驗中利用南京、杭州、廈門三個地區(qū)的雷達回波圖訓練DCNN-I,并將外推結果與降水預報業(yè)務中常用的COTREC算法和DITREC算法相比,結果表明,基于DCNN-I的雷達回波外推方法能夠提高外推回波圖像的準確率,并有效提高外推時效。
在下一步的研究工作中將進一步優(yōu)化DCNN-I,降低預測圖像的 “模糊” 程度,并將其應用于其他的圖像識別與預測領域。
參考文獻(References)
[1] DOVIAK R J, ZRNIC D S. Doppler Radar & Weather Observations [M]. Waltham, MA: Academic Press, 2014: 1-9.
[2] 張沛源, 楊洪平, 胡紹萍. 新一代天氣雷達在臨近預報和災害性天氣警報中的應用[J].氣象,2008,34(1):3-11.(ZHANG P Y, YANG H P, HU S P. Applications of new generation weather radar to nowcasting and warning of severe weather [J]. Meteorological Monthly, 2008, 34(1): 3-11.)
[3] OTSUKA S, TUERHONG G, KIKUCHI R, et al. Precipitation nowcasting with three-dimensional space-time extrapolation of dense and frequent phased-array weather radar observations [J]. Weather and Forecasting, 2016, 31(1): 329-340.
[4] RINEHART R E, GARVEY E T. Three-dimensional storm motion detection by conventional weather radar [J]. Nature, 1978, 273(5660): 287-289.
[5] 李英俊,韓雷.基于三維雷達圖像數據的風暴體追蹤算法研究[J].計算機應用,2008,28(4):1078-1080.(LI Y J, HAN L. Storm tracking algorithm developmentbased on the three-dimensional radar image data [J]. Journal of Computer Applications, 2008, 28(4): 1078-1080.)
[6] LIANG Q Q, FENG Y R, DENG W J, et al. A composite approach of radar echo extrapolation based on TREC vectors in combination with model-predicted winds [J]. Advances in Atmospheric Sciences, 2010, 27(5): 1119-1130.
[7] FLETCHER T D, ANDRIEU H, HAMEL P. Understanding, management and modelling of urban hydrology and its consequences for receiving waters: a state of the art [J]. Advances in Water Resources, 2013, 51(1): 261-279.
[8] 張亞萍,程明虎,夏文梅,等.天氣雷達回波運動場估測及在降水臨近預報中的應用[J].氣象學報,2006,64(5):631-646.(ZHANG Y P, CHENG M H, XIA W M, et al. Estimation of weather radar echo motion field and its application to precipitation nowcasting [J]. Acta Meteor Sinica, 2006, 64(5): 631-646.)
[9] KRIZHEVSKY A, SUTSKEVER I, HINTON G E. Imagenet classification with deep convolutional neural networks [C]// Proceedings of the 25th International Conference on Neural Information Processing Systems. [S.l.]: Curran Associates Inc., 2012,1: 1097-1105.
[10] 常亮,鄧小明,周明全,等.圖像理解中的卷積神經網絡[J].自動化學報,2016,42(9):1300-1312.(CHANG L, DENG X M, ZHOU M Q, et al. Convolutional neural networks in image understanding [J]. Acta Automatic Sinica, 2016, 42(9): 1300-1312.)
[11] 李倩玉,蔣建國,齊美彬.基于改進深層網絡的人臉識別算法[J].電子學報,2017,45(3):619-625.(LI Q Y, JIANG J G, QI M B. Face recognition algorithm based on improved deep networks [J]. Acta Electronica Sinica, 2017, 45(3): 619-625.)
[12] LeCUN Y, HUANG F J. Large-scale learning with SVM and convolutional for generic object categorization [C]// CVPR ’06: Proceedings of the 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Washington, DC: IEEE Computer Society, 2006, 1: 284-291.
[13] 張志強,劉黎平,王紅艷.三維可視化技術在雷達三維組網產品顯示中的運用[J].氣象科技,2010,38(5):605-608.(ZHANG Z Q, LIU L P, WANG H Y. Application of 3D visualization technology to display of Doppler radar networking products [J]. Metrological Science and Technology, 2010, 38(5): 605-608.)
[14] LU G Y, WONG D W. An adaptive inverse-distance weighting spatial interpolation technique [J]. Computers & Geosciences, 2008, 34(9): 1044-1055.
[15] GLOROT X, BENGIO Y. Understanding the difficulty of training deep feedforward neural networks [EB/OL]. [2017- 03- 01]. http://www.weblio.jp/redirect?url=http%3A%2F%2Fjmlr.org%2Fproceedings%2Fpapers%2Fv9%2Fglorot10a%2Fglorot10a.pdf&etd=3e03a9174d723be0.
[16] 王改利, 趙翠光, 劉黎平,等. 雷達回波外推預報的誤差分析[J]. 高原氣象, 2013, 32(3):874-883. (WANG G L, ZHAO C G, LIU L P, et al. Error analysis of radar echo extrapolation [J]. Plateau Meteorology, 2013, 32(3):874-883.)
[17] UIJLENHOET R, POMEROY J H. Raindrop size distribution and radar reflectivity-rain rate relationships for radar hydrology [J]. Hydrology & Earth System Sciences & Discussions, 2001, 5(4): 3012-3018.
[18] SHI X, CHEN Z, WANG H, et al. Convolutional LSTM network: a machine learning approach for precipitation nowcasting [C]// Proceedings of the 28th International Conference on Neural Information Processing Systems. Cambridge, MA: MIT Press, 2015,1: 802-810.
This work is partially supported by the National Natural Science Foundation of China (41305138, 61473310).
SHIEn, born in 1992, M. S. candidate. His research interests include neural network, pattern recognition.
LIQian, born in 1980, Ph. D., associate professor. His research interests include pattern recognition.
GUDaquan, born in 1959, professor. His research interests include intelligent information system, neural network.
ZHAOZhangming, born in 1993, M. S. candidate. His research interests include pattern recognition, distributed network.