基于ConvLSTM的氣象雷達(dá)回波外推

李英睿，李 瓊，魏加華，3*，喬 禛，申惟文

(1.青海大學(xué)土木工程學(xué)院，青海 西寧 810016；2.省部共建三江源生態(tài)與高原農(nóng)牧業(yè)國家重點(diǎn)實(shí)驗(yàn)室，青海大學(xué)，青海 西寧 810016；3.清華大學(xué)水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

雷達(dá)回波外推是臨近降水預(yù)報(bào)的主要手段，原理是利用雷達(dá)探測到的回波數(shù)據(jù)，確定回波的強(qiáng)度分布及回波體的移動(dòng)速度和方向，通過對回波體進(jìn)行線性或非線性的外推，預(yù)報(bào)一定時(shí)間段后的雷達(dá)回波狀態(tài)。傳統(tǒng)的雷達(dá)回波外推中交叉相關(guān)法[1-2]對數(shù)據(jù)利用率低，對于速度梯度大、演變快的回波，反演誤差大，難以推演實(shí)際天氣狀況。質(zhì)心跟蹤法[3-4]在穩(wěn)定性降水預(yù)報(bào)中可以取得較好效果[5]，但在局部地區(qū)對流天氣中，回波發(fā)展演變較快，無法滿足守恒條件，預(yù)報(bào)效果會(huì)隨時(shí)間的變化快速下降[6-7]。光流法[8]是從連續(xù)的圖像序列中計(jì)算光流場[9]，利用圖像序列中像素在時(shí)間域上的變化及相鄰幀之間的相關(guān)性，建立上一幀跟當(dāng)前幀之間存在的對應(yīng)關(guān)系，再計(jì)算出相鄰幀之間物體的運(yùn)動(dòng)信息。為提高預(yù)測效果，基于神經(jīng)網(wǎng)絡(luò)方法進(jìn)行臨近預(yù)報(bào)是目前研究的熱點(diǎn)[10]。Ayzel等[11]提出各層卷積核數(shù)量不同的6層全卷積Dozhdya.Net模型，采用零填充保持空間分辨率，與光流法相比在預(yù)報(bào)上具有更高的精度。陳家慧等[12]將BP神經(jīng)網(wǎng)絡(luò)模型應(yīng)用于混合型降水回波的外推，從歷史回波數(shù)據(jù)中學(xué)習(xí)得到回波運(yùn)動(dòng)的特征。施恩等[13]提出基于輸入的動(dòng)態(tài)卷積神經(jīng)網(wǎng)絡(luò)模型(DCNN-1)，與傳統(tǒng)回波外推方法相比，該方法有更高的預(yù)測精度且能有效延長外推時(shí)效。Shi等[14]提出基于LSTM的ConvLSTM(Convolutional Long Short-Term Memory，ConvLSTM)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，預(yù)測精度高于光流法。相比于傳統(tǒng)方法[1-4]，使用神經(jīng)網(wǎng)絡(luò)進(jìn)行回波預(yù)測，具有數(shù)據(jù)利用效率高、預(yù)報(bào)精度準(zhǔn)確、優(yōu)化潛力大等優(yōu)勢。

本文利用2019—2020年雷達(dá)體掃數(shù)據(jù)，通過數(shù)據(jù)預(yù)處理，建立樣本庫。利用ConvLSTM回波外推模型，實(shí)現(xiàn)基于雷達(dá)回波數(shù)據(jù)的預(yù)測，并與光流法對比，評價(jià)回波在不同時(shí)間上的預(yù)測精度，分析誤差隨外推時(shí)間及海拔高度的變化趨勢，探討深度學(xué)習(xí)模型在雷達(dá)回波預(yù)報(bào)上的適用性。

1 實(shí)驗(yàn)數(shù)據(jù)

雷達(dá)位于青海大學(xué)校內(nèi)(101.74°E，36.72°N)，海拔2 301 m。雷達(dá)有效探測半徑50 km，24 h以VOL方式運(yùn)行，5 min內(nèi)可完成從0.5°～19.4°的14個(gè)仰角360°掃描。雷達(dá)采用抬升仰角和波束掃描的觀測方式，回波由波束抬升和距離變化發(fā)生衰減。將三維極坐標(biāo)系的數(shù)據(jù)轉(zhuǎn)至笛卡爾坐標(biāo)系，并進(jìn)行波束遮擋以及雜波抑制等預(yù)處理工作，以提高雷達(dá)數(shù)據(jù)的利用率，最大可能地保留雷達(dá)原始數(shù)據(jù)的結(jié)構(gòu)特征。

1.1坐標(biāo)轉(zhuǎn)換假設(shè)目標(biāo)在雷達(dá)站A的極坐標(biāo)系中坐標(biāo)為(ρ，α，β)，其中ρ為斜距，α為方位角，θ為仰角。目標(biāo)在雷達(dá)站A直角坐標(biāo)系中的坐標(biāo)為(x，y，z)，則變換公式：

(1)

由式(1)投影后得到雷達(dá)的三維空間格點(diǎn)數(shù)據(jù)[15]。經(jīng)過坐標(biāo)轉(zhuǎn)換后的雷達(dá)網(wǎng)格數(shù)據(jù)分布散亂，為便于后續(xù)計(jì)算，采用最鄰近插值算法生成標(biāo)準(zhǔn)的網(wǎng)格數(shù)據(jù)，水平和垂直方向網(wǎng)格大小均為500 m×500 m。

1.2波束遮擋雷達(dá)波束在傳播過程中遇到障礙物(如高山、建筑物等)或地物干擾，導(dǎo)致雷達(dá)發(fā)射的電磁波束受到局部或全部遮擋。波束被遮擋會(huì)導(dǎo)致觀測信號失真，探測數(shù)據(jù)質(zhì)量受到影響。圖1a為雷達(dá)半徑50 km的掃描范圍中的DEM數(shù)據(jù)，低仰角時(shí)波束會(huì)被附近山體遮擋，同時(shí)地物雜波對雷達(dá)回波產(chǎn)生較大干擾。為減少波束遮擋影響，實(shí)測時(shí)調(diào)整雷達(dá)俯仰縮減探測范圍，以提高數(shù)據(jù)質(zhì)量，并以 20 km范圍內(nèi)的高仰角掃描數(shù)據(jù)(圖1b)為研究對象。此外，雷達(dá)回波預(yù)報(bào)時(shí)效以1 h內(nèi)最為理想，且回波與700 hPa的氣流有強(qiáng)相關(guān)性[16]。經(jīng)過綜合分析，選取海拔高3 000 m高度以上的數(shù)據(jù)進(jìn)行雷達(dá)回波外推。

圖1 雷達(dá)掃描區(qū)域及不同海拔波束遮擋情況Fig.1 Radar scanning area and beam occlusion at different altitudes

1.3地物雜波抑制地物雜波是雷達(dá)觀測中由地表物體(如植被、巖石、山體等)引發(fā)的異常噪聲信號，是復(fù)雜地形條件下經(jīng)常常遇到的雷達(dá)觀測誤差。固定地物雜波引起的噪聲會(huì)長期存在，對降雨的估計(jì)帶來一定誤差。地物雜波通常發(fā)生在距離雷達(dá)較近的區(qū)域內(nèi),且固定地物引起的雜波信號具有徑向速度接近為零等特征[17]。經(jīng)過綜合分析，在以雷達(dá)為中心的20 km的掃描范圍內(nèi)，將徑向速度等于零的值標(biāo)記為雜波干擾值，反射率因子的值視為無效值。圖2為2020年5月9日的晴空雷達(dá)回波數(shù)據(jù)，通常晴天雷達(dá)掃描過程中不會(huì)出現(xiàn)回波，故將晴空回波認(rèn)為是雜波干擾，圖2a和2b分別為雜波去除前和除去后的對比圖。

圖2 2020年5月9日晴空雷達(dá)回波去除雜波前后對比Fig.2 Comparison of radar echo of clear sky before and after removing the clutter on May 9，2020

2 實(shí)驗(yàn)方法

圖3 ConvLSTM結(jié)構(gòu)示意圖Fig.3 Schematic diagram of ConvLSTM structure

2.1ConvLSTM長短期記憶網(wǎng)絡(luò)(Long Short-Term Memory，LSTM)是一種常用來處理序列數(shù)據(jù)的時(shí)間循環(huán)神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Network，RNN)。本文選用Shi等[14]提出的ConvLSTM網(wǎng)絡(luò)結(jié)構(gòu),與用于一維時(shí)間序列數(shù)據(jù)的LSTM模型結(jié)構(gòu)不同，該模型將卷積與LSTM相結(jié)合，將信息傳遞時(shí)的全連接計(jì)算改為卷積計(jì)算，使其在時(shí)間和空間維度上均具有捕捉特征的能力。圖3為ConvLSTM的結(jié)構(gòu)，采用與LSTM相似的門控單元控制當(dāng)前輸入特征是否傳遞下去。ConvLSTM的3個(gè)門結(jié)構(gòu)分別為遺忘門(ft)、輸入門(it)和輸出門(Ot)，Ct是記憶細(xì)胞。記憶細(xì)胞中不僅保留了當(dāng)前的輸入特征，并通過控制前一時(shí)刻的信息是否繼續(xù)傳遞，信息傳遞關(guān)系用式(2)表示：

it=σ(Wxi*Xt+Whi*Ht-1+Wci°Ct-1+bi)ft=σ(Wxf*X+Whf*Ht-1+Wcf°Ct-1+bf)Ct=ft°Ct-1+it°tanh(Wxc*Xt+Whc*Ht-1+bc)Ot=σ(Wxo*Xt+Who*Ht-1+Wco°Ct+bo)Ht=Ot°tanh(Ct)

(2)

其中：σ表示sigmoid激活函數(shù)，(°)表示矩陣對應(yīng)元素相乘，*表示卷積操作。本文采用卷積操作對特征進(jìn)行提取，并通過sigmoid層和點(diǎn)乘運(yùn)算操作構(gòu)成的門結(jié)構(gòu)對信息進(jìn)行選擇。信息經(jīng)過遺忘門sigmoid單元，決定Xt和Ht-1中哪些信息需要丟棄，哪些信息需要繼續(xù)向后傳遞；繼續(xù)傳遞的信息流入輸入門通過sigmoid層決定需要更新的信息，通過tanh層得到新的細(xì)胞信息，進(jìn)行細(xì)胞更新；將通過輸出門sigmoid層的信息與通過tanh層的記憶細(xì)胞中的信息相乘得到模型的最終輸出Ht。

2.2光流法光流法將雷達(dá)回波映射到0～255的數(shù)值區(qū)間內(nèi)，基于光流法的短時(shí)回波預(yù)測，假設(shè)雷達(dá)在短時(shí)間內(nèi)的移動(dòng)趨勢不變，計(jì)算有位移區(qū)域的平均位移矢量，并將空白區(qū)域的位移設(shè)為平均值，得到各點(diǎn)位移矢量后，將矩陣中各點(diǎn)都按照位移矢量平移，得到新的雷達(dá)回波矩陣。因此，要預(yù)測 30 min后目標(biāo)點(diǎn)的像素值，需要迭代重復(fù)上述步驟6次，得到的坐標(biāo)點(diǎn)像素值即為預(yù)測目標(biāo)。

2.3對比方法及指標(biāo)選取2019年8月—2020年6月中22場降水的雷達(dá)回波數(shù)據(jù)對模型進(jìn)行訓(xùn)練，6場時(shí)間較長的降水過程作為模型的驗(yàn)證，并選取典型降水過程，分析誤差在時(shí)空上的變化。設(shè)計(jì)不同模型輸入時(shí)長，分別采用ConvLSTM與光流法進(jìn)行不同時(shí)間的回波預(yù)測，并對外推效果進(jìn)行對比分析。

由于難以直接量化回波圖像預(yù)測的準(zhǔn)確性，本文采用均方誤差MSE評價(jià)預(yù)測值和實(shí)際值的偏差，MSE數(shù)值越低，預(yù)測效果越好。計(jì)算公式見式(3)：

(3)

3 結(jié)果與分析

本實(shí)驗(yàn)通過改變模型的輸入時(shí)間(30、60、90、120 min)，采用ConvLSTM和光流法兩種方法分別進(jìn)行2 h(時(shí)間間隔30 min)的回波外推；固定ConvLSTM模型的輸入時(shí)間為30 min時(shí)，進(jìn)行2 h內(nèi)(每30 min)的回波預(yù)測；對選取典型過程回波外推效果和誤差進(jìn)行分析。

3.1回波外推時(shí)間過程分析選取2020年6月30日降水中不同時(shí)刻的雷達(dá)回波數(shù)據(jù)及評價(jià)指標(biāo)進(jìn)行對比，雷達(dá)回波可以直觀地反應(yīng)預(yù)測和實(shí)際觀測數(shù)據(jù)的整體變化趨勢；評價(jià)指標(biāo)反應(yīng)整場降雨中誤差隨時(shí)間的變化情況。

模型對輸入30、60、90和120 min的雷達(dá)回波數(shù)據(jù)分別進(jìn)行不同時(shí)間的外推。對2020年6月30日降水過程中19∶10、19∶40、20∶10、20∶40和21∶10時(shí)的回波進(jìn)行可視化(圖4)。第1行為觀測圖像，第2行為光流法外推后續(xù)第30 min在5個(gè)相應(yīng)時(shí)刻的外推回波圖像，第3行至第6行分別為對30(ConvLSTM-6)、60(ConvLSTM-12)、90(ConvLSTM-18)及120 min(ConvLSTM-24)的輸入時(shí)長，進(jìn)行第30分鐘、第60分鐘、第90分鐘和第120分鐘外推在相應(yīng)5個(gè)時(shí)刻的回波圖像。由圖4可知，基于ConvLSTM模型預(yù)測結(jié)果與實(shí)際觀測圖像相似程度更高；在回波數(shù)值較低的區(qū)域中，光流法很難準(zhǔn)確計(jì)算出回波強(qiáng)度較弱區(qū)域的矢量場，使得部分結(jié)果被低估。但該模型在一定程度上仍存在誤差，這是因?yàn)槔走_(dá)回波外推本身有不確定性，雖然ConvLSTM在訓(xùn)練過程中可以捕捉回波在時(shí)空上的變化規(guī)律，但無法對區(qū)域內(nèi)所有格點(diǎn)做出完全精準(zhǔn)的預(yù)測，某些位置只能反映出變化趨勢。此外，隨著外推時(shí)間的增加，模型的預(yù)測誤差逐步累計(jì)，從21∶10回波圖像中可以直觀地看出ConvLSTM-18和ConvLSTM-24對回波預(yù)測準(zhǔn)確度在不同程度上有所下降，雖然總體趨勢與觀測圖像類似，但仍有不一樣的細(xì)節(jié)特征。

圖4 各時(shí)刻回波圖像對比Fig.4 Comparison of images of echo at different times

選取2020年6月的6場降水過程作為驗(yàn)證集，均方誤差隨時(shí)間的變化如圖5a所示。從圖中可以看出，隨著外推時(shí)間的增加，預(yù)測誤差也增加，ConvLSTM預(yù)測結(jié)果相對穩(wěn)定，該模型對特征的學(xué)習(xí)及預(yù)測的誤差逐步累積，而光流法的整體預(yù)測誤差均相對較大。

圖5 光流法和ConvLSTM的預(yù)測誤差對比Fig.5 Comparison of prediction errors between optical flow method and ConvLSTM

綜上所述，ConvLSTM預(yù)測結(jié)果優(yōu)于光流法，且可以有效延長預(yù)測時(shí)間。ConvLSTM模型不僅具有傳統(tǒng)全連接LSTM網(wǎng)絡(luò)學(xué)習(xí)時(shí)間特征上的能力，而且還具備捕捉空間特征的能力。光流法在捕捉時(shí)間和空間變化特征的能力上稍顯不足，原因在于雷達(dá)回波數(shù)據(jù)并非一一對應(yīng)，存在較強(qiáng)的非線性、時(shí)滯性。光流法要求遵循圖像灰度不變的假設(shè)，雷達(dá)回波各時(shí)刻都存在升消變化，在不結(jié)合云層運(yùn)動(dòng)變化的情況下，會(huì)出現(xiàn)反射率因子不守恒帶來的誤差，導(dǎo)致模型難以對測試集進(jìn)行準(zhǔn)確預(yù)測。

3.2回波外推空間預(yù)測精度分析由于雷達(dá)波束抬升，探測高度增大，探測距離和不同海拔回波的預(yù)測精度存在差異，本實(shí)驗(yàn)選取2020年6月25日的降水過程進(jìn)行分析回波外推空間預(yù)測的精度。ConvLSTM模型與光流法預(yù)測結(jié)果在不同海拔的誤差分布見圖6和圖7，從左到右、從上到下分別為海拔3 000～5 500 m外推誤差空間分布圖，層高間隔為500 m。圖中ConvLSTM模型在各高度層的誤差值均小于光流法，呈現(xiàn)較穩(wěn)定的預(yù)測精度，隨著探測高度增加，誤差逐漸增大。假定實(shí)際值和預(yù)測值≤5 dBz時(shí)滿足精度要求，ConvLSTM模型在1～3層中約有96%的預(yù)測誤差處于0～5 dBz，4～6層中預(yù)測誤差處于0～5 dBz的概率下降，第6層中誤差≤5 dBz的概率僅有49%。盡管光流法的預(yù)測精度較低，68%的格點(diǎn)預(yù)測誤差在該范圍內(nèi)，但隨著高度的增加預(yù)測誤差仍舊呈現(xiàn)上升趨勢?？梢?，隨著探測高度的增加，預(yù)報(bào)的誤差逐步累積。

圖6 ConvLSTM模型第30分鐘預(yù)測在不同海拔高度預(yù)測誤差累積Fig.6 Cumulation of prediction errors at different altitudes in the 30min prediction of ConvLSTM

圖7 光流法第30分鐘預(yù)測在不同海拔高度預(yù)測誤差累積Fig.7 Cumulation of prediction errors at different altitudes in the 30min prediction of optical flow method

圖8 平均預(yù)測誤差對比Fig.8 Comparison of average prediction error

兩種方法在不同海拔高度中不同預(yù)測時(shí)長的平均誤差如圖8所示。海拔3 000～4 000 m高度中的誤差均相對較低；海拔3 000 m高度時(shí)的預(yù)測效果最好，海拔高于4 000 m誤差呈現(xiàn)明顯的上升趨勢。光流法各層的預(yù)測誤差均值相對較高。在ConvLSTM模型中，利用海拔4 500 m及以上的雷達(dá)回波數(shù)據(jù)，預(yù)報(bào)準(zhǔn)確率下降約為原來的60%左右。從圖5b可以看出雷達(dá)回波外推在水平方向變化，距離雷達(dá)中心越遠(yuǎn)，誤差越大。經(jīng)過綜合分析可知，雷達(dá)在本地區(qū)海拔4 000 m以下回波預(yù)測效果較好；在海拔3 000 m的預(yù)測誤差最小，表明回波與700 hPa的氣流有很強(qiáng)的相關(guān)性。

選取2020年6月25日的降水過程分析不同預(yù)測時(shí)長在水平方向上的誤差變化。模型的輸入時(shí)間恒為30 min，計(jì)算海拔3 000 m高度層各時(shí)刻第30分鐘，第60分鐘，第90分鐘及第120分鐘后的平均預(yù)測誤差，結(jié)果如圖9所示。

圖9 不同外推時(shí)間預(yù)測誤差變化Fig.9 Variation of prediction error of different extrapolation time

圖10 預(yù)測誤差隨時(shí)間變化趨勢Fig.10 Trend of forecast error over time

圖9a～圖9d分別為各時(shí)刻第30分鐘，第60分鐘，第90分鐘，第120分鐘的預(yù)測誤差累積在海拔3 000 m上的分布情況。從圖中可以直觀地看出，隨著外推時(shí)間的增加，誤差的累積也在逐漸增大。在圖9a～圖9d中分別有96.3%、95.5%、95.1%、81.7%的預(yù)測誤差小于5 dBz，相對誤差大于5 dBz的占比也由4%上升至18%左右。圖10中，在海拔3 000 m高度層上誤差隨著外推時(shí)間的增加逐步累積，預(yù)測誤差小于5 dBz的占比整體呈現(xiàn)波動(dòng)下降趨勢，相對誤差在數(shù)值上逐漸增大。整體來看，對于第30～第90分鐘的預(yù)測精度誤差的整體變化趨勢不大，但第60～第90分鐘的預(yù)測精度比第30分鐘略高，在第120分鐘時(shí)預(yù)測的誤差值有明顯的增加，預(yù)測精度較低。

4 討論與結(jié)論

雷達(dá)回波外推技術(shù)是臨近天氣預(yù)報(bào)的主要技術(shù)手段，目前深度學(xué)習(xí)網(wǎng)絡(luò)在預(yù)報(bào)領(lǐng)域的應(yīng)用逐漸增多。本文利用ConvLSTM模型進(jìn)行雷達(dá)回波外推，檢驗(yàn)該模型對青海大學(xué)X波段天氣雷達(dá)回波的預(yù)測能力。由于降水的形成與發(fā)展的過程較為復(fù)雜，預(yù)報(bào)存在較多的不確定因素。從本次實(shí)驗(yàn)結(jié)果綜合來看，該模型對西寧地區(qū)的雷達(dá)回波數(shù)據(jù)有著較好的預(yù)測精度，在海拔3 000 m高度層的預(yù)測精度較高，回波與700 hPa的氣流有很強(qiáng)的相關(guān)性，該結(jié)果與陸漢城等[18]所得結(jié)論一致。與傳統(tǒng)的光流法的預(yù)報(bào)結(jié)果相比，ConvLSTM模型通過訓(xùn)練能從大量的雷達(dá)數(shù)據(jù)中學(xué)習(xí)到更多的雷達(dá)回波變化時(shí)間和空間特征，有效地提高了數(shù)據(jù)的利用率，對于雷達(dá)回波數(shù)據(jù)有較強(qiáng)的適應(yīng)性。Shi等[14]通過對比ConvLSTM模型與光流法及其改進(jìn)方法的預(yù)測精度，該模型的精度高于其他對比方法，這與本文研究結(jié)果相同。從ConvLSTM模型的外推結(jié)果來看，在第1.5小時(shí)內(nèi)的表現(xiàn)較好，誤差較低，預(yù)報(bào)精度滿足應(yīng)用需求，能夠在總體上預(yù)報(bào)降水的變化趨勢，這與梁振清等[19]得出ConvLSTM 網(wǎng)絡(luò)在1h內(nèi)表現(xiàn)較為良好的結(jié)論稍有出入，這可能與西寧地區(qū)的降水類型、降水持續(xù)時(shí)間及云層的影響有關(guān)。相比于光流法，ConvLSTM模型能在一定程度學(xué)習(xí)到特征在空間及時(shí)間上的變化趨勢，并在一定程度上延長外推時(shí)長。ConvLSTM模型對有降水過程的雷達(dá)回波預(yù)測表現(xiàn)較好，在未來可以針對不同降水強(qiáng)度下雷達(dá)回波的變化對模型進(jìn)行適當(dāng)?shù)母倪M(jìn)，進(jìn)一步提高模型的準(zhǔn)確性與適用性。