張毓杉，周 曉

(武漢理工大學(xué) 機(jī)電工程學(xué)院，湖北 武漢 430070)

天氣預(yù)報(bào)分為長期預(yù)報(bào)，中期預(yù)報(bào)，短期預(yù)報(bào)，短時(shí)預(yù)報(bào)和臨近預(yù)報(bào)5種[1]。本文研究的臨近預(yù)報(bào)是指的2 h以內(nèi)的預(yù)報(bào)。臨近預(yù)報(bào)主要是通過雷達(dá)回波圖像的外推技術(shù)來實(shí)現(xiàn)的，對(duì)于雷達(dá)回波歷史圖像，經(jīng)過算法分析處理后，預(yù)測生成未來時(shí)刻的雷達(dá)回波圖像。雷達(dá)回波圖像是灰度圖像，其灰度值大小反映了當(dāng)前的降雨情況[2]，圖像上灰度值越大的區(qū)域，表示此處觀測到的回波值越大，也說明了此處降雨程度越嚴(yán)重。1978年Rinehart[3]提出了使用交叉相關(guān)法來跟蹤雷達(dá)回波。1993年，Dixon[4]提出了基于雷達(dá)回波數(shù)據(jù)對(duì)云團(tuán)進(jìn)行識(shí)別、追蹤和預(yù)報(bào)的TITAN(thunderstorm identification tracking analysis and nowcasting)算法。后來Johnson[5]提出了SCIT(storm cell identification and tracking)算法，對(duì)TITAN算法進(jìn)行了改進(jìn)。因?yàn)楣饬鞣ǖ呐d起，Bowler[6-7]等提出利用光流算法來外推雷達(dá)回波圖像來實(shí)現(xiàn)降雨預(yù)測的目的。

筆者以雷達(dá)回波降雨預(yù)測為研究目標(biāo)，結(jié)合傳統(tǒng)的光流法[8]和卷積神經(jīng)網(wǎng)絡(luò)提出了基于卷積神經(jīng)網(wǎng)絡(luò)的光流生成網(wǎng)絡(luò)(convolutional neural networks-flow, CNN-Flow)生成雷達(dá)回波圖像間的光流場。最后再利用光流場和半拉格朗日算法預(yù)測出未來時(shí)刻的圖像序列。

1 卷積光流生成網(wǎng)絡(luò)

卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural networks, CNN)擁有各種各樣的結(jié)構(gòu)，相對(duì)于全連接神經(jīng)網(wǎng)絡(luò)，由于其權(quán)重共享的特性，大大減少了訓(xùn)練時(shí)的參數(shù)數(shù)量，因此選擇CNN作為基礎(chǔ)結(jié)構(gòu)來生成光流場。假設(shè)雷達(dá)回波圖像的寬與長分別為w和h，且由于雷達(dá)回波圖像是灰度圖像，其通道數(shù)為1，則兩個(gè)輸入圖像的尺寸均為h×w×1。而光流場中每個(gè)點(diǎn)代表了圖像中每個(gè)像素點(diǎn)的橫縱移動(dòng)速度分量，因此網(wǎng)絡(luò)輸出的光流場尺寸為h×w×2，網(wǎng)絡(luò)的總體工作流程如圖1所示。

圖1 CNN-Flow網(wǎng)絡(luò)工作流程

從圖1可知，輸入網(wǎng)絡(luò)的兩幀圖像首先進(jìn)行通道融合，在通道維度上進(jìn)行疊加，再輸入網(wǎng)絡(luò)，經(jīng)過編碼解碼器的上下采樣U型結(jié)構(gòu)后，最后生成對(duì)應(yīng)的光流場。

1.1 網(wǎng)絡(luò)結(jié)構(gòu)

該網(wǎng)絡(luò)由編碼器和解碼器組成。編碼器的作用不僅是為了減輕訓(xùn)練壓力，也壓縮了圖像的尺寸，使得卷積層可以獲得更大的感受野，提取更加深層的特征。解碼器的作用是恢復(fù)圖像的尺寸，因?yàn)槔走_(dá)回波各個(gè)部分移動(dòng)的速度可能不一樣，有的區(qū)域位移過大，有的區(qū)域位移過小，為了能夠同時(shí)處理不同移動(dòng)幅度的回波，網(wǎng)絡(luò)在解碼過程的同時(shí)在不同尺寸的圖像上利用步長為1的卷積層提取不同尺度的光流場，然后將不同尺度的光流場同時(shí)與經(jīng)過縮放的真實(shí)光流場進(jìn)行損失值的計(jì)算，最后網(wǎng)絡(luò)輸出與原圖像尺寸一樣大小的光流場，CNN-Flow的具體結(jié)構(gòu)如圖2所示。

圖2 CNN-Flow網(wǎng)絡(luò)結(jié)構(gòu)

從圖2可知，由于網(wǎng)絡(luò)層數(shù)過少會(huì)使得提取的特征不足，層數(shù)過多對(duì)訓(xùn)練結(jié)果也不會(huì)有太大的提升，反而提高了訓(xùn)練難度。因此該網(wǎng)絡(luò)一共設(shè)有6層編碼層和6層解碼層，其中Conv代表卷積層，Deconv代表反卷積層，Concat代表連接層(用于將多個(gè)特征圖按通道維度進(jìn)行融合)，flow1-6代表不同尺寸大小的光流場，這些光流場均是在解碼過程中利用卷積層對(duì)不同分辨率的特征圖進(jìn)行提取而生成的，由于光流場的通道數(shù)是2(代表兩個(gè)速度分量)，因此提取光流場的卷積層的卷積核數(shù)量均為2。

CNN-Flow的編碼器由10個(gè)卷積層疊加而成，每一層的作用是向下一層或者同層解碼層傳遞特征，解碼器由6個(gè)反卷積層及6個(gè)連接層組成，解碼器的每一層輸入的特征都由上一層解碼層輸出的特征以及上一層解碼層提取的光流場和同一層的編碼層輸出的特征通過跳層連接融合而成的。網(wǎng)絡(luò)訓(xùn)練過程中將所有提取的不同尺寸的光流場與相對(duì)應(yīng)尺寸的真實(shí)光流場進(jìn)行損失函數(shù)計(jì)算，然后將所有損失值進(jìn)行加權(quán)和。由于LeakyReLU[9]激活函數(shù)解決了神經(jīng)元死亡的問題，因此網(wǎng)絡(luò)選取了LeakyReLU作為激活函數(shù)。

1.2 多權(quán)重組合損失函數(shù)

由于光流場的每一個(gè)格點(diǎn)均由兩個(gè)速度分量u和v組成，為了使得訓(xùn)練時(shí)輸出的光流場接近真實(shí)的光流場，使用如式(1)所示的指標(biāo)EPE(endpoint error,端點(diǎn)誤差)作為損失函數(shù)的基本計(jì)算公式。

(1)

網(wǎng)絡(luò)中一共輸出了6個(gè)不同尺度的光流場作為輸出，最后需要計(jì)算6個(gè)不同尺度的EPE，然后獲得EPE加權(quán)損失函數(shù)。尺寸越大的光流場所計(jì)算出的EPE的權(quán)重越大，加權(quán)損失函數(shù)計(jì)算公式如式(2)所示。

Loss=a1×EPE1+a2×EPE2+a3

×EPE3+a4×EPE4+a5*EPE5+

a6×EPE6

(2)

式中：EPE1為尺度最小的光流場；EPE6為尺度最大的光流場。經(jīng)過多次嘗試比較，最后權(quán)重設(shè)置如下：a1=0.004，a2=0.01，a3=0.02，a4=0.06，a5=0.1，a6=0.31。

2 雷達(dá)回波圖像預(yù)測

通過CNN光流網(wǎng)絡(luò)獲取光流場之后，還需要利用光流場以及計(jì)算出光流場的原始兩幀圖像進(jìn)行圖像外推計(jì)算。此處外推計(jì)算使用的是半拉格朗日外推法，一維的半拉格朗日算法為：現(xiàn)假設(shè)某一粒子的F值由其坐標(biāo)x和時(shí)間t所決定，即F(x,t)，則半拉格朗日的基本方程如式(3)和式(4)所示：

(3)

(4)

式(3)和式(4)表明F在沿著粒子運(yùn)動(dòng)軌跡時(shí)是一個(gè)恒定不變的值。U(x,t)代表著當(dāng)前點(diǎn)的速度。

圖3為半拉格朗日外推方法的大致模型。圖3表示了粒子的移動(dòng)情況，AC這條粗線表示了粒子從t-Δt時(shí)刻到t+Δt時(shí)刻的真實(shí)移動(dòng)路徑；aC這條細(xì)線表示了粒子從t-Δt時(shí)刻到t+Δt時(shí)刻的近似移動(dòng)路徑。A和a作為粒子的起始點(diǎn)；B是粒子的中間點(diǎn)；C是粒子的終點(diǎn)。一般終點(diǎn)在圖像格點(diǎn)上，而起點(diǎn)和中間點(diǎn)都不會(huì)剛好在網(wǎng)格點(diǎn)上。Δt是粒子一次移動(dòng)的時(shí)間間隔；2α是粒子從t-Δt時(shí)刻到t+Δt時(shí)刻的移動(dòng)距離。

圖3 半拉格朗日外推方法的大致模型

由圖3可以得到式(5)和式(6)：

(5)

α=ΔtU(x-α,t)

(6)

式(6)中α值可由式(7)進(jìn)行迭代求出。

(7)

這樣就能在t-Δt時(shí)刻找到t+Δt時(shí)刻對(duì)應(yīng)的F值。

對(duì)于雷達(dá)回波的外推問題，因?yàn)槭菆D像外推，所以需要將一維的半拉格朗日擴(kuò)展成二維，現(xiàn)假設(shè)粒子(圖像中的像素點(diǎn))的F值(即灰度值)由其坐標(biāo)x和y以及時(shí)間t決定，二維半拉格朗日的方程變?yōu)槿缦滦问剑?/p>

(8)

(9)

(10)

式(8)～式(10)中U為像素點(diǎn)x方向的速度；V為像素點(diǎn)y方向上的速度。于是像素點(diǎn)的運(yùn)動(dòng)軌跡可由式(11)～式(13)決定：

(11)

α=ΔtU(x-α,y-β,t)

(12)

β=ΔtV(x-α,y-β,t)

(13)

式(12)和式(13)依然可用迭代方式求出。由于起始點(diǎn)和中間點(diǎn)不可能會(huì)剛好落在圖像整數(shù)坐標(biāo)上，因此F(x-2α,y-2β,t-Δt)，U(x-α,y-β,t)，V(x-α,y-β,t)需要通過圖像插值的方式求出。假設(shè)(i,j)，(i,j+1)，(i+1,j)，(i+1,j+1)4個(gè)點(diǎn)分別為鄰近的4個(gè)格點(diǎn)，而起始點(diǎn)和中間點(diǎn)會(huì)落在(i+u,j+v)處，u和v均為浮點(diǎn)數(shù)，可采用雙線性插值法，即通過式(14)求出。

f(i+u,j+v)=(1-u)×(1-v)×f(i,j)+

u×(1-v)×f(i+1,j)+(1-u)×v×

f(i,j+1)+u×v×f(i+1,j+1)

(14)

通過上述方法就能利用光流場以及原始的兩幀圖像推出第三幀的圖像。由于CNN-Flow結(jié)合半拉格朗日外推算法后只能通過輸入兩幀圖像預(yù)測出第三幀圖像，筆者采用遞歸調(diào)用的方式生成10幀預(yù)測圖像，通過輸入最新生成的圖像以及上一幀圖像，重新計(jì)算光流場并推算下一幀圖像。

3.1 數(shù)據(jù)集

采用的數(shù)據(jù)集來自于深圳市氣象局提供的粵港澳大灣區(qū)的雷達(dá)回波圖像數(shù)據(jù)集(SRAD)，SRAD共包含20 000個(gè)樣本，其中每個(gè)雷達(dá)數(shù)據(jù)個(gè)案樣本覆蓋時(shí)長為4小時(shí)，樣本中每幀圖像采集時(shí)間間隔6分鐘，圖像長和寬均為256(即大約覆蓋255 km×255 km的區(qū)域)。

設(shè)計(jì)的CNN-Flow需要兩幀連續(xù)的雷達(dá)回波圖像作為輸入，而輸出是一幀兩通道的光流場。網(wǎng)絡(luò)要求的數(shù)據(jù)集的每一個(gè)樣本應(yīng)該包含兩張雷達(dá)回波圖像以及一個(gè)光流場文件，但是由于目前還沒有能夠直接檢測到光流的傳感器，無法建立標(biāo)準(zhǔn)的雷達(dá)回波圖像的光流數(shù)據(jù)集，因此筆者首先將SRAD進(jìn)行滑動(dòng)窗口切分，保持每個(gè)樣本中包含有兩幀連續(xù)的雷達(dá)回波圖像，再對(duì)其打亂并隨機(jī)篩選切分成了8 000個(gè)訓(xùn)練樣本及2 000個(gè)測試樣本，最后用Farneback算法處理每一個(gè)樣本的兩幀圖像，生成對(duì)應(yīng)的光流場文件加入對(duì)應(yīng)的樣本中，然后將新的數(shù)據(jù)集稱為FSRAD(flow SRAD)。

3.2 評(píng)價(jià)指標(biāo)

實(shí)驗(yàn)采取的是降水臨近預(yù)報(bào)評(píng)估方法，其中包含有CSI(critical success index,臨界成功指數(shù))，F(xiàn)AR(false alarm rate,誤報(bào)率)，POD(probebility of detection,檢測概率)[10]3個(gè)指標(biāo)來評(píng)判預(yù)報(bào)的準(zhǔn)確性。首先通過設(shè)定一個(gè)降雨閾值，將預(yù)測圖像與真實(shí)圖像中在閾值內(nèi)的格點(diǎn)設(shè)置為下雨點(diǎn)(即標(biāo)記為1)，閾值外的格點(diǎn)設(shè)置為無雨點(diǎn)(即標(biāo)記為0)，然后按像素遍歷預(yù)測圖像及真實(shí)圖像，當(dāng)預(yù)測值為1，真實(shí)值為1時(shí)記為一個(gè)hits；當(dāng)預(yù)測值為1，真實(shí)值為0時(shí)記為一個(gè)alarms；當(dāng)預(yù)測值為0，真實(shí)值為1時(shí)記為一個(gè)misses。最后統(tǒng)計(jì)hits，alarms以及misses的個(gè)數(shù)。式(15)為3個(gè)指標(biāo)的計(jì)算公式。

(15)

式中：nhits為hits點(diǎn)的個(gè)數(shù)；nalarms為alarms點(diǎn)的個(gè)數(shù)；nmisses為misses點(diǎn)的個(gè)數(shù)。其中CSI與POD的值越高代表預(yù)測越準(zhǔn)確，F(xiàn)AR的值越低代表預(yù)測誤報(bào)率越低。

3.3 模型訓(xùn)練參數(shù)設(shè)置

網(wǎng)絡(luò)訓(xùn)練超參數(shù)設(shè)置如表1所示。

表1 網(wǎng)絡(luò)的超參數(shù)

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

由于EPE指標(biāo)可用于驗(yàn)證生成的光流場及真實(shí)光流場之間的差距，且EPE值越小代表生成的光流場越接近真實(shí)值，因此該網(wǎng)絡(luò)在訓(xùn)練過程中記錄了迭代過程中在FSRAD訓(xùn)練集及測試集上的平均EPE數(shù)值變化，結(jié)果如圖4所示。

圖4 EPE變化曲線

從圖4可知，迭代了160次之后，訓(xùn)練集和測試集的平均EPE已經(jīng)開始收斂，并且訓(xùn)練集的EPE收斂于0.9201，測試集的EPE收斂于0.724 4。表2為CNN-Flow和傳統(tǒng)光流法在FSRAD測試集上生成一幀光流場所需的平均時(shí)間對(duì)比。

表2 生成光流場的平均時(shí)間消耗

從表2可知，CNN-Flow速度明顯快于傳統(tǒng)的光流法，因此基于深度學(xué)習(xí)的光流法在運(yùn)行效率上有了較大的提升。為了驗(yàn)證CNN-Flow在預(yù)測圖像上的降雨預(yù)測精度，將生成的光流場結(jié)合半拉格朗日外推法，通過兩幀雷達(dá)回波圖像預(yù)測未來10幀的圖像，實(shí)驗(yàn)中將回波閾值設(shè)為30 dB，以此閾值在CNN-Flow和傳統(tǒng)光流法預(yù)測的10幀圖像以及實(shí)際的10幀圖像上用3個(gè)降雨指標(biāo)進(jìn)行對(duì)比分析。表3為兩種算法在測試集上所計(jì)算的3種指標(biāo)的均值，可以看出CNN-Flow在CSI指標(biāo)和FAR指標(biāo)上均略差于傳統(tǒng)光流法，在POD指標(biāo)上強(qiáng)于傳統(tǒng)光流法，結(jié)合表2和表3分析可得，CNN-Flow網(wǎng)絡(luò)在總體精度上和傳統(tǒng)光流法差不多的情況下，運(yùn)行速度遠(yuǎn)遠(yuǎn)超過了傳統(tǒng)光流法。

表3 算法在測試集上的平均指標(biāo)

為了測試算法在個(gè)例樣本上的預(yù)測效果，選取了某個(gè)時(shí)刻的個(gè)例樣本繪制了兩種算法的CSI，F(xiàn)AR，POD指標(biāo)隨著預(yù)測時(shí)間增加的變化曲線，如圖5所示，其中橫坐標(biāo)代表輸入前兩幀圖像，預(yù)測3至12幀圖像。從圖5可知，CNN-Flow和傳統(tǒng)光流法在預(yù)測時(shí)間較短時(shí)，精度都較高，但是當(dāng)預(yù)測時(shí)間較長時(shí)，精度都嚴(yán)重下降，這是由于CNN-Flow和傳統(tǒng)光流法都是僅僅分析了兩幀圖像而預(yù)測后續(xù)圖像，隨著時(shí)間增加，不確定性會(huì)逐漸增加，例如回波的突然消失和誕生，這種情況導(dǎo)致了精度的嚴(yán)重下降。CNN-Flow在CSI指標(biāo)上和傳統(tǒng)光流法在預(yù)測前期幾乎相等，隨著時(shí)間增加，逐漸略強(qiáng)于傳統(tǒng)光流法；在FAR指標(biāo)上，略差于傳統(tǒng)光流法；在POD指標(biāo)上始終略強(qiáng)于傳統(tǒng)光流法。

圖5 評(píng)價(jià)指標(biāo)曲線圖

為了在視覺角度分析本文算法在預(yù)測圖像上的效果，圖6是CNN-Flow和傳統(tǒng)光流法在不同樣本上預(yù)測的10幀圖像中第1、4、7、10幀圖像的對(duì)比。從圖6可知，兩個(gè)算法均能成功預(yù)測出回波圖像的正確移動(dòng)趨勢，并且在預(yù)測時(shí)間較短時(shí)，預(yù)測的回波和真實(shí)的回波幾乎一樣，但是當(dāng)預(yù)測時(shí)間增加時(shí)，兩種算法對(duì)于圖像細(xì)節(jié)的預(yù)測能力均有減弱。

圖6 算法預(yù)測圖像對(duì)比

4 結(jié)論

筆者在卷積神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上結(jié)合傳統(tǒng)的光流法提出了使用CNN-Flow網(wǎng)絡(luò)及半拉格朗日算法進(jìn)行雷達(dá)回波圖像預(yù)測，該算法采用了SRAD數(shù)據(jù)集作為原數(shù)據(jù)進(jìn)行實(shí)驗(yàn)分析，并在時(shí)間、降雨指標(biāo)以及圖像上對(duì)算法進(jìn)行了效果的驗(yàn)證，在實(shí)驗(yàn)對(duì)比中發(fā)現(xiàn)該算法能夠較準(zhǔn)確地預(yù)測雷達(dá)回波圖像，并在計(jì)算速度上遠(yuǎn)遠(yuǎn)超過了傳統(tǒng)光流法。