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

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

0 引言

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

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

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

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

1 ConvLSTM模型

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

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

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

ConvLSTM內部結構如圖3所示。

其對應公式如下：

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

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

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

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

2 實驗過程

2.1 數(shù)據(jù)準備

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

（1）切割

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

（2）灰度轉換

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

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

（3）得到訓練集

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

2.2 ConvLSTM模型建立與分析

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

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

ConvLSTM模型樣本數(shù)據(jù)的輸入尺寸為如圖8所示的 5D 張量（samples，time, rows，cols，channels），要提前將訓練集和測試集reshape成如上形式的tensor張量，即樣本總數(shù)、各樣本幀數(shù)、圖像寬度、圖像高度和顏色通道數(shù)。輸出尺寸為，如果返回全部序列，則返回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的輸入。

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

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

（2）ConvLSTM模型建立

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

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

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

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

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

2.3 ConvLSTM模型訓練與預測

（1） 模型訓練

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

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

（2）預測

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

（3）可視化結果

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

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

2.4 算法評價

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

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

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

3 結語

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

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

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