基于ConvLSTM的雷達回波外推預報

高云，郭艷萍，周建慧，張葉娥，楊澤民

山西大同大學計算機與網絡工程學院，山西大同，037009

0 引言

隨著全球氣候變暖，工業污染加重，全球生態遭到嚴重破壞，引發的氣候異常問題也越來越嚴重，中國的自然災害是世界上最嚴重的國家之一。在各種重大自然災害中，冰雹、暴雨、雷暴等突發性強的氣象災害，發生頻率最高；又因其程度強而持續時間短而難以提前預知和防范，導致危害程度較大。如2005年07月，印度西部連日特大暴雨引發多起洪水和泥石流等災害，造成至少900人死亡。2005年06月，黑龍江省寧安市洪災死亡人數達106人，給人民的生命財產帶來很大危害。

強對流天氣是具有重大破壞性的災害性天氣之一。目前，通常使用雷達回波外推的方法對強對流天氣進行監測和預警。雷達回波外推是指對回波運行的速度、方向、強度和形態變化等特征進行跟蹤和預測[1-2]。傳統雷達回波外推方法存在預測回波變化較快的天氣過程時準確率大幅下降，預測時效較短，隨著預測時長的增加準確性快速下降等不足。

早期的外推方法最常使用TREC技術，即交叉相關法。Hilst等先后在文獻[3-6]中應用TREC技術對整體移動的風暴簇進行了跟蹤，結果表明該技術能夠對風暴簇的內部運動進行較好的反演。后經過多位學者的改進，徐亞欽等在文獻[7]中提出了多重動態區域TREC算法，對TREC方法進行改進。單體質心法是較常用的一種回波外推法。喬春貴等在文獻[8]中對整幅雷達回波圖通過單體質心法進行線性外推，得出對于穩定的層狀云該方法雖然有較好的外推能力，但對其他類型的云則外推結果較差。1950年Gibson在文獻[9]中首先提出了光流法，張蕾等在文獻[10]中對傳統光流法進行了改進。陳家慧等在文獻[11]中使用BP神經網絡進行混合型降水回波的外推，表現出了較好的外推能力。

雷達回波外推方法近年來雖然取得了較大的進步，外推預報的準確性逐漸提高，但目前的雷達回波外推算法，對于1小時以上的外推預報，準確性下降很快。雷達回波外推方法還需要進一步地探索和改進。

1 ConvLSTM模型

ConvLSTM是LSTM結構的變體，通過將W的權值計算變成了卷積運算來提取出圖像的特征。圖1所示即為其單元結構圖。

ConvLSTM的思路就是使用卷積來代替矩陣乘法。它在輸入到狀態以及狀態到狀態轉換中都采用了卷積結構，通過堆疊多個ConvLSTM層并形成編碼預測結構來建立更一般的時空序列預測模型。

該設計的一個顯著特點是所有輸入X1,...,Xt，細胞輸出C1,...,Ct，隱藏狀態H1,...,Ht，和ConvLSTM的3個門it，ft，ot都是3維張量，2維網格向3維張量的轉換如圖2所示。

ConvLSTM內部結構如圖3所示。

其對應公式如下：

其中星號*表示卷積，小圓圈○表示哈達瑪乘積。可以看出，這里在輸入門、遺忘門和輸出門這3個門的輸入上，都考慮了細胞狀態Ct-1。

為了確保細胞狀態Ct-1具有與輸入相同的行數和列數，在應用卷積運算之前需要先進行padding。將邊界點上隱藏狀態的填充視為使用外部世界的狀態進行計算。通常，在第一個輸入到來之前，將LSTM的所有狀態初始化為零，這對應于對于未來的“完全無知”。

ConvLSTM模型用于預測網絡時其編解碼結構如圖4所示。

ConvLSTM也可以作為更復雜結構的構建塊。對于時空序列預測問題，使用圖 4所示的結構。它包括編碼網絡和預測網絡2個網絡，將編碼網絡的最后狀態復制到預測網絡的初始狀態和單元輸出，2個網絡都是通過堆疊多個ConvLSTM層形成的。預測目標與輸入具有相同的維度，將預測網絡中的所有狀態連接起來并發送到1×1卷積層以生成最終預測。

2 實驗過程

2.1 數據準備

本算法進行分析的原始數據序列為山西省大同市2019-6-1～2019-6-6共1191幅雷達站基本反射率圖。雷達站名：大同；雷達掃描間隔時間為5′47"；數據范圍200km；顯示仰角：0.5/1.5/2.4。以原始基數據的命名方式進行命名，方便后期進行按順序索引。其中2019-6-5 16:14:55BJT基本反射率圖如圖5所示。

（1）切割

從圖 5可以看出，雷達站大同掃描范圍為200km，除覆蓋大同地區以外，還覆蓋了山西省其他地區以及內蒙古自治區和河北省等部分地區。本實驗只對大同地區的強對流氣候進行預測，排除其他區域的干擾，先設定圖片輸出分辨率和初始化最終數組，處理的第一步是對原始雷達回波圖像進行切割，只切割原始圖像（200,150, 310, 350）的區域作為本次實驗的數據。切割后的圖像如圖 6所示，只保留了大同市周邊地區。同時減少比色卡、坐標、圖片標題等干擾項，增加模型訓練的精度。

（2）灰度轉換

第二步是將三通道的彩色圖片轉換為單通道灰度圖片。由于計算機硬件限制，原始繪制的300dpi高分辨率圖像不適合輸入到神經網絡中直接進行訓練，首先要將裁切后的圖像降低分辨率（壓縮）至100*100。原始雷達回波圖像為RGB三通道圖像，為減少訓練時間可將原始圖片轉換為灰度圖像，每個像素用8個bit表示，RGB轉換為灰度圖像的公式為：

將圖像抗鋸齒，取消字體平滑，抗混疊。每個雷達站每個時刻的體掃數據均可得到一個100*100分辨率大小的灰度圖像，轉換為灰度的圖像如圖7所示。

（3）得到訓練集

增加一個維度備用，每幅圖像載入內存后即為一個100*100*1大小的三維數組。處理完單個體掃數據之后再依次遍歷所有圖片，將每個圖片得到的數組連接成一個四維的數組序列。為了排除雜點的干擾，再次將數組轉化為2維（圖像數量，100*100），依次遍歷數組中每個點的值，如果該點的值＜50，則置為0。得到所需的圖像數組序列。最后將所得數組保存為.npz格式的numpy數組方便后續取用，這就是最終可以輸入到網絡中進行訓練的樣本集。

2.2 ConvLSTM模型建立與分析

使用ConvLSTM模型進行雷達回波圖像外推步驟如下：

（1）訓練集和測試集劃分

ConvLSTM模型樣本數據的輸入尺寸為如圖8所示的 5D 張量（samples，time, rows，cols，channels），要提前將訓練集和測試集reshape成如上形式的tensor張量，即樣本總數、各樣本幀數、圖像寬度、圖像高度和顏色通道數。輸出尺寸為，如果返回全部序列，則返回5D張量，即（samples, timesteps, output_row, output_col, filters）；否則，返回4D張量，即（samples,output_row, output_col, filters）。o_row和o_col取決于filter和padding的尺寸。

假如上一層是ConvLSTM2D layer，那么其輸出為以上形式的4D張量或5D張量，當后面再接另外一個layer時，就要考慮該layer是否能接受4D張量或5D張量，即要考慮ConvLSTM2D的輸出能否作為該layer的輸入。

數據準備生成的.npz格式的numpy數組是2維的，首先將其讀入并reshape為4維數據，即（NUMBER, WIDTH, HEIGHT, 1），樣本集總數為NUMBER，數據集為灰度圖像，將顏色通道數設為1。將數據集設置為2個5D張量的數組，BASIC_SEQUENCE和NEXT_SEQUENCE，分別存放當前時刻圖像和下一時刻圖像作為訓練集和結果集，數據為5維(NUMBER - FRAMES,FRAMES, WIDTH, HEIGHT, 1)，樣本總量=原樣本總數-各樣本幀數。

訓練時將整個樣本集切割成長度統一的樣本，16（FRAMES）幀圖像為一個單元。共分為（NUMBER –FRAMES）個單元，即可直接訓練。

（2）ConvLSTM模型建立

卷積長短時記憶神經網絡自提出以來逐步完善已經成為了較為成熟的圖像序列預測模型，目前Python中也有不少對其模型的封裝，方便研究人員直接構建網絡對自己的數據進行訓練和預測。

本文所建模型為如圖 9所示的Keras中的序貫模型，即逐層嵌套依次連接，數據流由輸入端輸入，逐層流動，反向傳播時更新參數，逐步降低損失函數。該網絡由四層ConvLSTM2D網絡層堆疊而成，最后為一個三維卷積層用以格式化輸出數組以便求取損失函數或獲得預測結果。

本文所建模型由4層ConvLSTM2D網絡層+1層ConvLSTM3D網絡層堆疊而成，每層卷積核數目不同。4層ConvLSTM2D網絡層卷積核數目分別為40、40、60、40；卷積核大小均為3*3；補0策略為′same′，即輸出shape與輸入shape相同；使用線性激活函數；返回輸出序列中的全部序列。ConvLSTM3D網絡層卷積核數目分別為1；卷積核大小為3*3*3；補0策略為′same′，即輸出shape與輸入shape相同；使用′sigmoid′激活函數；輸出中維度的順序為通道在后面。

為了保證在線性向非線性轉變時，權重的尺度不變，所以使用BatchNormalization在激活函數前對輸入進行了批標準化操作，優化神經網絡，將分散的數據統一。

由于sigmod函數具有梯度飽和現象的缺點，所以本實驗使用的是學習率為0.01的隨機梯度下降優化器optimizers.SGD，將所有參數梯度裁剪到數值范圍為[-0.5,0.5]的區間內。

2.3 ConvLSTM模型訓練與預測

（1） 模型訓練

在訓練模型之前，需要通過compile來對學習過程進行配置。compile接收三個參數，分別為優化器、損失函數和指標列表。優化器optimizer為一個已預定義的優化器名，或者一個類型為Optimizer的對象；損失函數loss為目標函數，目的模型優化，可以使用預定義的損失函數名或一個自定義的損失函數；指標列表metrics的設置是針對分類問題，一般設metrics值為′accuracy′。指標可以是一個預定義指標的名字，也可以是一個用戶定制的函數。指標函數應該返回單個張量，或一個完成metric_name - ＞ metric_value映射的字典。

本模型的損失函數為對數損失函數′binary_crossentropy′，以adadelta為優化方法進行編譯。為了降低訓練時間，在Keras中創建一個多GPU模型，設定4GPU并行處理。輸入的數據即為2.1中處理的BASIC_SEQUENCE和NEXT_SEQUENCE作為訓練集和結果集序列；batch_size值為10，即計算梯度所需的樣本數量為10；epochs值為10，即代表樣本集內所有的數據經過了10次訓練；validation_split=0.05，即樣本集中5%的數據為驗證集。訓練樣本集，取其中一回合的部分訓練結果如圖10所示。

（2）預測

模型訓練完成后存入“nice_model.h5”文件中，取樣本集中第600組圖片中的12幀進行預測，預測16次，訓練結果連接形成預測集。

（3）可視化結果

將預測結果與訓練所用樣本集及其對應的結果集對比形成可視化結果。結果分為16組圖片，前8組圖片顯示初始軌跡與標記數據的對比結果，如圖11所示為第4幅圖片。

后8組圖片顯示預測結果與標記數據對應結果的對比結果，如圖 12所示為第12幅圖片。

2.4 算法評價

從圖10所示的模型訓練結果來看，損失函數趨于減小的趨勢，最終穩定在一個較小的值，說明模型的擬合程度較好。從圖11的初始軌跡對比來看，軌跡運動與標記數據吻合度很好，從圖12的預測結果來看，預測圖像與標記數據的結果集吻合度基本一致，說明預測結果正確性好。

值得注意的是，模型圖像的處理對后期結果的影響很大，由于不同時段雷達回波反射率因子分布區間不同，在進行灰度轉換時，若使用原始數據繪制圖像可能使比色卡（color map）的范圍不同，這將導致不同圖像中的相同強度的反射率因子值對應的RGB色彩不同，若不進行處理將極大影響訓練效果。

另外，為了得到更為準確的外推結果，可在模型訓練完成后，進行3次以上的外推，可以提高預測精度。

3 結語

本算法使用大同市歷史氣象雷達數據，重點介紹了ConvLSTM模型在雷達回波外推短期降水預測的應用方法，并詳細介紹了ConvLSTM模型在雷達回波外推即時空序列分析建模的整個過程，同時對相應的算法和整個建模流程給出了實現及結果說明。本算法可以從以下方面進行拓展思考：

（1）本實驗在建模時使用的損失函數是對數損失函數，損失函數有很多種，實際應用中均方誤差函數也是常用的一種損失函數。后續實驗中應考慮如何優化損失函數，選擇使模型擬合性更好的損失函數。

（2）ConvLSTM是Xingjian Shi等提出的，后面他針對雷達圖的天氣預測又提出了TrajGRU，基于運行軌跡對圖像做更精準的捕捉。本實驗如何使用TrajGRU模型進行更好的預測，這些問題都需要進一步地進行探討。