深圳市氣象局臨近預報技術進展

陳元昭 蘭紅平 劉琨

（深圳市氣象局，深圳 518030）

0 引言

基于雷達數據的雷暴識別追蹤和外推預報方法是最早出現的臨近預報技術[1]。20世紀50年代，Ligda[2]首先開始了雷達回波外推進行臨近預報研究，此后臨近預報技術得到不斷發展，主要體現在一些預報方法的提出。單體質心法的核心是首先將整個雷暴單體視為一個三維的單體，然后對單體進行識別、分析后，再進行追蹤，采用擬合外推法來做雷暴的臨近預報[3-5]。交叉相關法是通過計算整塊回波不同區域的連續時次的相關系數，采用最優相關法，來獲得回波的運動矢量場，然后通過這些運動矢量場來預報回波未來的位置和形狀。交叉相關法的優勢在于它的計算方法相對簡單，不但可以跟蹤對流降水，還可以對層狀云降水進行追蹤[6-8]，在氣象業務部門一度得到廣泛應用[9-10]。但交叉相關法對強度和形狀隨時間變化很快、局地生成的回波，常會出現跟蹤失敗的情況[11-12]。近年來，基于雷達的臨近預報方法發展迅速，中國氣象局在完成災害性天氣短時臨近預報系統（SWAN）[13]的基礎上進一步完善了強對流天氣的監測技術，研發了基于閃電和衛星產品的臨近預報產品[14]。光流法是計算機視覺領域中的重要方法[15]，近年來該方法在天氣臨近預報領域得到廣泛應用[16]，但對于西風帶系統，尤其是颮線的預報，光流法存在一定的局限性[17]；而未來幾年以雷達外推技術為主的臨近預報將繼續在業務中應用[18-19]。近年來，隨著計算能力、信息技術和智能算法技術的不斷突破，人工智能技術備受社會各界的高度重視，發展迅猛[20-21]，基于深度學習的智能天氣臨近預報方法也受到氣象領域的高度關注。

本文主要從技術角度對深圳市氣象局的臨近預報技術的發展歷程進行回顧總結。第一部分主要介紹交叉相關法、光流法、粒子濾波臨近預報方法以及雷達基數據質量控制方法，第二部分主要介紹基于生成對抗網絡（generative adversarial networks，GAN）[22]的人工智能臨近預報方法，最后對不同方法的優缺點進行總結和討論。

1 基于運動矢量場的臨近預報方法研究進展

深圳市氣象局臨近預報方法的研究始于2006年，在引進和吸收了國內外對流風暴追蹤技術算法的基礎上，2008年研發了交叉相關臨近預報方法，并應用到2011年深圳世界大學生運動會氣象服務中。隨后逐漸研發了局部約束光流法、粒子濾波融合法等臨近預報技術。以上述方法為核心搭建了臨近預報決策支持平臺（PONDS）系統，并在業務中實時運行，對廣東地區的強對流天氣過程進行實時監測和臨近預報，成為深圳市氣象局三大核心業務系統之一。

1.1 雷達資料處理

深圳臨近預報采用廣東省12部S波段天氣雷達（分別位于廣州、深圳、韶關、珠海、清遠、陽江、河源、汕尾、汕頭、梅州、湛江、肇慶）2.5 km高度的反射率因子雷達拼圖資料。雷達的虛假回波經常出現在低層，選取2.5 km高度的回波拼圖作為預報的初始場。這一層的雷達回波既能代表華南地區對流的水平分布特征，同時在很大程度上避免了虛假回波和地物雜波的影響。為計算方便，雷達資料從極坐標格式利用最近鄰居法和垂直方向的線性內插法相結合的插值法插值到三維直角坐標系中。拼圖的時間間隔為6 min。拼圖前對雷達基數據資料進行質量控制，剔除了超折射、地物等的影響。

1.2 交叉相關臨近預報方法

交叉相關法的核心就是通過計算連續時次雷達回波不同區域的最優空間相關系數，來確定回波在過去的移動矢量特征。交叉相關法把整個數據區域劃分成若干小區域，然后在相鄰時刻雷達回波圖像的小區域之間計算相關系數，通過最大相關系數確定相鄰時刻圖像中的區域對應關系，進而確定回波區域的平均運動。一個簡單的示意圖如圖1所示。將相鄰時刻（t1，t2）雷達反射率因子回波圖像（圖1），分別劃分成若干大小相同的子區域。對于t1時刻的一個子區域A，對應t2時刻回波圖像的候選區域（大小由平均風速決定），將所有可能的子區域分別與A做相關計算，相關系數R表示為：

圖1 交叉相關法示意圖 Fig. 1 Cross-correlation approach

式中，1和2分別是1和2時刻的反射率因子，是一個子區域內象素點的數量。t2時刻相關系數最大的子區域B即為與A匹配的子區域。將A和B的中心位置連接起來，即為回波運動矢量。

交叉相關法僅利用最近兩個時次的雷達回波獲得的運動矢量場是不穩定的，用這種風場外推的預報結果也具有較大誤差。雷達標定、地物及異常傳播的電磁波受到阻礙物遮擋、衰減、雷達波束的不完全充塞等隨機因素是出現這種誤差的主要原因。在進行運動矢量場計算時，需要抑制風場中的這些噪聲干擾，以提高回波預報穩定性和準確性。在交叉相關法中利用卡爾曼濾波法改進雷達回波交叉相關算法[23]，具體算法見[24-25]，計算中考慮了過去0.5～1h的回波趨勢。從效果上看，提升了回波預報的穩定性[23]。

在變化較為平緩的層狀云降水系統中，從雷達回波圖像上觀察到的雷暴的外形和移動速度變化都不是很劇烈，此時交叉相關法的效果較好[11]。但對強對流天氣過程，尤其是在雷暴的外形和移動速度隨時間變化很快的情況下，通過簡單的計算相關系數的方法難以保證追蹤的準確性，交叉相關法的預報效果就會明顯降低，跟蹤失敗的情況顯著增加，并最終影響預報的結果[11]。所以，對變化較為劇烈的強對流降水系統，交叉相關法有先天的缺陷，需要引入新的方法。

1.3 局部約束光流法臨近預報方法

1.3.1 Lucas-Kanade 法計算光流場

光流法的基本原理是以圖像亮度變化來識別運動目標和觀測器之間的相對運動產生的瞬時位移。圖像中所有像素點的瞬時位移就構成了圖像的光流場。而光流法的核心就是通過識別出的連續圖像系列計算光流場。簡單的來說，光流場就是類剛體物體的速度矢量場。光流方程如下：

計算光流場的約束方法通常有Lucas-Kanade局部約束[26]和Horn-Schunk全局約束[27]。局部約束方法是對某個點周圍給定小區域的光流給定限定條件，而全局約束方法是在整個圖像區域范圍內滿足一定的約束條件。臨近預報關注回波的局地變化，光流法中采用Lucas-Kanade局部約束法作為計算光流的約束條件。用獲得的光流場進行臨近預報。

與傳統的交叉相關法相比，光流法的優勢在于立足于變化，而不是選定不變特征再跟蹤不變特征移動的方式。

1.3.2 質量控制——中值濾波

在光流法中，對雷達數據的質量控制采用中值濾波法。中值濾波是一種非線性平滑技術，它將每一像素點的灰度值設置為該點某鄰域窗口內的所有像素點灰度值的中值[28]。

2013年4月30日，受鋒面低槽影響，一條颮線自西北向東南橫掃廣東全省。本次過程回波濾波前后對比可以看出（圖2），濾波前后颮線的長度、形狀、強回波中心等特征大體一致。但可以看出，濾波前（圖2a），颮線后部南端的層狀云回波比較零碎，雜亂。濾波后（圖2b）去噪效果好，圖像清晰，較好地保留了回波的邊緣特征，沒有出現明顯失真。

圖2 2013年4月30日14時（北京時，下同）廣東雷達拼圖光流法追蹤得到的回波：（a）中值濾波前；（b）中值濾波后 Fig. 2 Optical flow method based on the Guangdong multi-radar mosaic at 14:00 BT 30 April 2013: (a) before median filtering; (b) after median filtering

1.3.3 預報效果個例對比

2014年5月20日中午，一條近東西向的弓形回波從珠江口西側向東移動，追上原位于其東側的正在減弱的回波。這兩條回波合并后，迅速增加，在珠江口附近發展成為倒“Y”字形的強回波。給珠江口一帶帶來了局地強降水和短時雷雨大風。

光流法對這次強降水過程30 min外推預報效果較好。從圖3可知，光流法預報出了回波合并后強度加強，預報的30 min后回波的位置與實況基本相符，也基本預報出了回波的倒“Y”字形結構，只是西南部的回波預報比實況偏弱。

對華南地區回波預報評估結果表明[16]，光流法整體預報效果優于交叉相關法，尤其對移動型局地生成的回波及強度和形狀隨時間變化很快的雷暴。

光流法對西風帶系統，尤其是颮線的預報仍存在一定的局限性[19]，需要研制新方法來彌補雷達臨近預報技術在0～1 h的預報能力方面的缺陷。

1.4 粒子濾波臨近預報方法

1.4.1 粒子濾波法基本原理

粒子濾波（Particle filter）由Carpenter等[29]首次提出，主要是通過非參數化的蒙特卡洛（Monte Carlo）模擬來實現遞推貝葉斯（Bayes）濾波估計。

在進行雷達外推預報前需獲得雷達回波的運動矢量場。和獲得唯一運動矢量場不同，粒子濾波算法是利用同一時次的雷達回波采用兩種不同的估測方法得到的相同時次的兩組回波運動矢量場，然后采用粒子濾波算法對兩組運動矢量場進行融合，再利用得到的運動矢量場來進行回波的外推預報。在粒子濾波算法中，回波矢量場可表示為：

圖3 光流法對2014年5月20日珠江口附近迅速加強回波預報：（a）14時實況；（b）14時的30 min外推預報；（c）14：30雷達回波實況 Fig. 3 Optical flow method based on the Guangdong multi-radar mosaic on 20 May 2014:(a) real radar echo at 14:00 BT; (b) 30 min extrapolation forecast at 14:00 BT; (c) real radar echo at 14:30 BT

式中，f(x)、h(x)分別為狀態方程和觀測方程。x(t)、y(t)、m(t)、n(t)分別表示系統的狀態、觀測值、過程噪聲和觀測噪聲。其中f(x)為采用Lucas-Kanade約束的光流法得到的運動矢量場，而h(x)采用基于Harris角點算法[30]追蹤回波得到運動矢量場。

圖4 粒子濾波法流程圖 Fig. 4 Flow chat of particle filtering fusion algorithm

圖4為粒子濾波法流程圖。粒子濾波計算方法的核心是用一組加權隨機樣本（粒子）來近似表示后驗概率密度函數。遞推貝葉斯濾波理論給出了該問題的嚴格求解。在進行求解時引入了優化算法對粒子濾波參數進行尋優，使算法對整個參數空間進行高效搜索以獲得最優解。相比其他方法，粒子濾波可以得到更優化的回波運動矢量場。粒子濾波利用一系列帶權值的空間隨機采樣，來逼近后驗概率密度函數。

1.4.2 質量控制—雙邊濾波

1998年，Overton等[31]最先提出了雙邊濾波（bilateral filtering）的概念，它是在基于空間分布的高斯濾波函數基礎上提出的。傳統高斯濾波僅僅關注像素的空間距離而忽略了像素值的變化，因此會在濾除噪聲的同時造成邊緣的模糊。雙邊濾波在傳統高斯濾波器的基礎上添加了考慮像素值變化程度的權重。因此，在濾除圖像噪聲的同時考慮到灰度相似度的信息，使得權重系數隨著圖像灰度的變化而改變，有效地保持了圖像邊緣結構，最終能夠進行自適應的濾波。

研究選取2016年4月13日颮線過程的回波進行濾波試驗。從濾波前的回波（圖5a）可以看出，強回波的邊界、強回波區的對流回波及部分層狀云回波比較零碎、雜亂，部分回波呈鋸齒狀。從濾波后回波對比可以看出，中值濾波和雙邊濾波后（圖5b和5c）回波的形狀、強回波中心以及強回波后的層狀云回波和實況大致一樣，強回波的邊緣、強回波中心，后部的層狀云回波變得更平滑、有序。對比濾波后颮線前沿回波，中值濾波后颮線前沿的回波被削弱，而雙邊濾波較好保留了前沿回波。中值濾波可以改善回波的質量，但是會削弱回波的邊沿，尤其是颮線的邊沿，在進行回波運動矢量場計算時，邊沿的風矢量場只能由估計分析得來，從而影響運動估計的質量。而雙邊濾波在提高回波質量的同時，也很好地保留了回波邊沿，在進行運動矢量計算時可以直接獲得該處的風矢量場，質量相對來說更可靠。

圖5 2016年4月13日雷達回波濾波前后對比：（a）濾波前；（b）中值濾波后；（c）雙邊濾波后 Fig. 5 Comparison between the Guangdong multi-radar mosaic on 13 April 2016: (a) before filtering; (b) after median filtering; (c) after bilateral filtering

1.4.3 運動矢量場對比

獲得真實、平滑的運動矢量場是進行外推預報的關鍵。分別用交叉相關法、光流法、粒子濾波法追蹤2016年4月13日06：30颮線過程回波的運動矢量場，其中交叉相關法和光流法均采用中值濾波法對雷達基數據進行質量控制，粒子濾波法采用雙邊濾波進行質量控制。

該颮線在06：30—07：30的1 h內大致呈西北—東南方向移動，粒子濾波融合法預報的移動路徑與實況大致相同。在實況中颮線的不同位置移動方向有所差別，颮線南段向偏東南方向移動，北段向偏東移方向移動，東南段介于兩者之間。粒子濾波融合算法追蹤的颮線運動矢量風場（圖6a）中，颮線南段大致為西北氣流控制，北段偏西氣流控制為主，東南段為西北氣流控制，與颮線相應位置的移動方向基本一致，對矢量場的刻畫也更精細。交叉相關法追蹤的矢量場大致為偏西風（圖6b），預報了颮線在1 h內朝偏東方向移動。而光流法追蹤的颮線南段為偏西風，北段為西南風（圖6c），預報颮線在1 h內朝偏東北方向移動。交叉相關法和光流法在進行運動矢量追蹤時，由于質量控制的原因，強回波中心或者旋轉處運動矢量會被周邊的運動矢量所平滑，導致預報的移動路徑出現偏差。由此可知，粒子濾波融合法獲得了更接近實況、更真實的、更精細的回波運動矢量。

圖6 2016年4月13日06：30運動矢量場對比：（a）粒子濾波融合法；（b）交叉相關法；（c）光流法 Fig. 6 Comparison between motion vector field at 06:30 BT 13 April 2016: (a) particle filtering fusion algorithm;(b) cross-correlation approach; (c) optical flow method

檢驗結果表明[19]，粒子濾波法能對各中天氣類型的回波進行較好的預報，尤其是颮線預報效果好于光流法和交叉相關法。

2 人工智能臨近預報方法研究

2.1 人工智能臨近預報方法研究進展

20世紀80年代開始，人工智能技術已經在氣象領域得到探索性的應用，從氣候資料的處理分析到預報產品的制作都有涉足，其中主要有[32]基于人工神經網絡、基于遺傳計算及其融合算法、基于支持向量機方法、基于貝葉斯方法等氣象預報方法。由于計算能力、數據存儲分析能力、以及人工智能算法等基礎條件的限制，盡管當時的一些成果已顯示出較好的應用前景，人工智能技術在氣象領域的進展仍受到了制約[20]。

面對信息技術和智能預報技術的快速發展，國內氣象界也在積極探索將智能預報技術與氣象業務科技的發展相結合，并開展了積極的嘗試。如中央氣象臺、福建省氣象臺等氣象部門通過合作開展基于卷積神經網絡、決策樹模型等算法的智能臨近預報方法研究。香港天文臺從2014年開始和香港科技大學合作開發基于卷積神經網絡預測模型的深度學習降水臨近預報方法[33]。通過機器深度學習，模擬雷達回波未來2 h的移動路徑，較傳統基于光流法預測雷達回波移動具有了新的優勢。

2.2 基于生成對抗網絡GAN 的智能預報技術

基于生成對抗網絡GAN的智能預報技術是從一系列雷達觀測資料中，運用圖像識別等方法提取知識和信息建立預測模型，并通過損失函數訓練模型，輸出產品，進行預報的人工智能臨近預報新技術。該方法可快速和智能化地識別回波，從而預報未來一段時間內對流天氣發生和影響區域（圖7）。

圖7 基于GAN的人工智能臨近預報方法示意圖 Fig. 7 Artificial intelligent nowcasting approach based on GAN

GAN由谷歌大腦團隊的Goodfellow等[22]提出，近幾年在學術界受到廣泛的關注，但目前研究主要集中在傳統的視覺領域，即利用GAN模型生成傳統視覺領域圖片，如鳥、貓、狗等。深圳市氣象局首次使用GAN方法來做回波臨近預報。

2.2.1建 模模型時的，建需立對及原訓始練雷達圖片進行5次卷積[34]。其中第一到第四次卷積分別對輸入圖片采用卷積層和糾正線性函數ReLU[35]進行卷積，分別獲取圖片特征，第五次卷積采用卷積層和雙曲正切函數[36]進行卷積，獲得第五特征圖，再對5次獲取的卷積進行堆疊完成模型。建模應用了2015—2017年共3年的雷達數據。

通過堆疊5個卷積層從原始圖像中提取的特征，建立預報模型，輸出預報目標，通過計算預報目標與真實目標圖像之間的損失函數，再通過損失函數對模型進行訓練，促使模型向正確的方向優化。

嚴格來說，一個GAN框架，最少擁有兩個組成部分，一個是生成器G（Generator），一個是判別器D（Discriminator）。生成器G基于模型生成預報回波圖像序列。判別器則需要對生成器生成的圖片進行真假判別。在訓練過程中，會把生成器生成的樣本和真實樣本信息傳遞給判別器D。判別器D的目標是盡可能正確地識別出真實樣本，和盡可能正確地揪出生成的樣本，也就是假樣本。而生成器的目標則和判別器相反，就是盡可能最小化判別器揪出它的概率。

生成器G 和判別器D 組成一個最小最大游戲（min-max game），在訓練過程中雙方都不斷優化自己，直到達到平衡——雙方都無法變得更好，也就是假樣本與真樣本完全不可區分，從而完成模型的優化訓練。

GAN方法本身沒有明確的“預報”過程。在模型的訓練過程中，使用歷史序列雷達圖片和1 h內預測目標的數據對作為訓練樣本。這樣模型在訓練中學到的輸出就是未來1 h的預報。對抗學習模型見圖8。

圖8 基于GAN的智能臨近預報對抗學習外推模型 Fig. 8 Adversarial learning model based on GAN

2.2.2 模型的效果

受南海西北部熱帶低壓和西南氣流的共同影響，2018年6月5日廣東大部分地區出現了暴雨降水。降水主要出現在5日上午后。低壓外圍的對流云團不斷涌現，給廣東帶來持續的間歇性降水。

選取當日12時的回波進行60 min的外推預報試驗。對比這次降水過程雷達反射率因子拼圖13時的實況（圖9c）和智能預報方法在12時起報的13時外推預報（圖9b），智能預報方法預報出了主要的四塊回波：粵西一帶的東北—西南向的條狀回波、珠江口西側的偏西北—東南向的條狀回波、珠江口東側的回波及粵東沿海的條狀回波。預報的回波位置與實況大致相當，但強度偏弱，強回波的范圍偏大。

圖9 2018年6月5日反射率因子拼圖實況和預報：（a）16：12實況；（b）16：12的60 min預報；（c）17：12實況 Fig. 9 Comparison between the Guangdong multi-radar mosaic on June 5, 2018: (a) real radar echo at 16:12 BT;(b) 60 min forecast at 16:12 BT; (c) real radar echo at 17:12 BT

相對于天氣過程的非線性、多因子作用的特點，只用雷達數據做預報顯然是不夠的。在進行模型訓練時，除了雷達數據，還應該包含大氣動力信息、熱力信息，甚至地形信息等，以提高模型的適應能力，提高智能預報水平。對算法的改進，可以考慮在模型中加入Dense Net結構、Batch Normalization批歸一化層等，以提高模型的性能。

3 預報效果評估

多年臨近預報方法的研究，極大提升了深圳市氣象局臨近預報預警水平。對流天氣預警時間提前量由30 min左右提高到2018年的56 min。為了對臨近預報效果進行定量檢驗，分別計算粒子濾波法、光流法和GAN智能臨近預報方法的擊中率（POD）、空報率（FAR）、臨界成功指數（CSI），具體算法見文獻[16]。選取2018年的18個降水過程，包括西風帶、西南季風、臺風降水等，將回波強度分為25～80，35～80，45～80 dBz三個等級。結果如圖10。

圖10 粒子濾波法、光流法及GAN方法1 h預報試驗結果擊中率、空報率和臨界成功指數對比 Fig. 10 Comparation of hit rate, false forecasting rate and critical success index (CSI) in 1 hour later forecasting experiment via particle filtering method, optical flow and GAN

由圖10可知，三種方法對25 dBz以上中等強度回波的POD都在55%以上，其中光流法最高，GAN方法次之，略低于光流法。而成功臨近指數CSI中，光流法和GAN方法大致相當。由此可知，總體優勢上，GAN和傳統的光流法大致相當。對光流法，三種強度回波的預報效果均好于GAN方法和粒子濾波法，目前業務應用以光流法為主。但對不同的過程，三種方法預報效果不盡相同，光流法對局地生成及強度隨時間快速變化的回波預報效果較好，粒子濾波法對颮線的預報效果較好，而GAN方法對颮線和臺風的預報也具有較大的參考價值。

4 結論和討論

近些年來，深圳市氣象局逐漸開發了基于運動矢量場的交叉相關法、光流法、粒子濾波融合法臨近預報技術，并與哈爾濱工業大學深圳研究生院合作開展了基于生成對抗網絡（GAN）的人工智能臨近預報方法研究。

基于運動矢量場的臨近預報方法中，三種預報方法對于不同天氣過程各有優勢，業務中很難有一種方法能很好地做出所有天氣過程的臨近預報，對于不同的天氣過程，只能通過預報員的判斷來選擇采用哪種預報方法?；谏蓪咕W絡（GAN）的人工智能臨近預報方法大致可以預報1 h的回波，預報結果有一定的參考價值。

目前，基于雷達回波的外推臨近預報技術是基于雷達反射率因子單個因子，并沒有考慮影響風暴生消的動力學、熱力學等因素，幾乎沒有能力預報回波強度的變化，對于生命史小于預報時效的雷暴，預報意義不大。而對于風暴生命史較長的過程，比如超級單體風暴、持續幾個小時的強颮線、由鋒面造成的強降水雨帶以及臺風等大尺度天氣過程，外推臨近預報具有一定的現實意義。

對于人工智能臨近預報方法，災害性天氣事件如冰雹、龍卷、臺風等的氣象觀測數據目前尚未達到深度學習需要的樣本數量，所以基于深度學習的智能臨近預報方法的研究是一個持續的過程，需要持續的投入研究。

對于臨近預報未來發展主要依賴以下幾方面：一是數值模式的發展，將先進的外推預報方法同快速更新循環的高時空分辨率數值模式預報以及二者的融合，是未來強對流天氣短時臨近預報的重要發展方向[13]；二是大數據的挖掘，發展基于人工智能的臨近預報方法，希望可以尋找到一些規則并應用到業務服務中；三是風暴追蹤算法的不斷改進；四是預報員能力的逐漸提高。