胡 鵬,伍瑞林,伍光勝,張志堅

(1.廣州市突發(fā)事件預警信息發(fā)布中心，廣州 511430; 2.廣州市荔灣區(qū)氣象局，廣州 510140)

0 引言

雷電是一種常見的天氣現(xiàn)象，常伴隨著強對流等天氣過程出現(xiàn)，容易造成破壞性極大的雷災損失。廣東省地勢北高南低，珠三角三面環(huán)山，南朝南海，呈一個半封閉低谷，山地的向陽坡和迎風坡有利于海洋水汽抬升，冷空氣來臨時容易出現(xiàn)鋒面雷暴天氣，尤其廣州市地處珠江三角洲，城市范圍大，人口分布稠密，是廣東省雷擊風險高值區(qū)之一[1]。經(jīng)統(tǒng)計，2019年廣州市共發(fā)生云對地閃電約23萬次，平均地閃回擊密度值為30.64次/(平方公里·年)，平均地閃強度為10.44 kA，發(fā)生雷災事故37起，占全省雷災總數(shù)的14.5%，較2018年增加48%，因雷電災害造成的經(jīng)濟損失215.41萬元，其中電力核電行業(yè)是雷災的重災區(qū)，雷災起數(shù)最多，造成的經(jīng)濟損失最大。因此，加強對雷電災害產(chǎn)生的分析和研究，尤其是對雷暴運動發(fā)展過程的研究，從而實現(xiàn)高精度的雷電預警，對進一步完善城市雷電災害防御等方面有著重要意義。

目前，國內(nèi)外均主要利用閃電定位資料、大氣電場資料和雷達回波資料開展雷電預警相關技術研究。如美國、加拿大、巴西、法國、英國、日本等均建立了閃電監(jiān)測定位網(wǎng)，其中美國于20世紀七八十年代就已開始建設雷電預警系統(tǒng)[2]。在我國，中國氣象科學研究院近幾年也開發(fā)出了綜合利用雷達、衛(wèi)星、閃電定位、大氣電場及探空等多種資料的雷電臨近預警系統(tǒng)[3]，并投入業(yè)務運行；梁巧倩、林良勛根據(jù)雷電發(fā)生的天氣學分型和完全預報方法研究出了一種廣州地區(qū)雷電短期預報方案，對雷電潛勢預報具有較高參考價值[4]；王凱等通過利用黃山周邊地區(qū)雷暴過程大氣電場儀和閃電定位儀數(shù)據(jù)，選取預報因子建立了雷電預警的預報方程，獲得較好的應用效果[5]；路郁等對閃電定位數(shù)據(jù)進行聚類分析，利用“膨脹-腐蝕”算法還原雷暴云，并外推未來雷電落區(qū)[6]。以上雷電預警預報研究都取得了一定進展，但雷電發(fā)生時具有一定的隨機性、瞬間性、地域相關性等特征，增加了其規(guī)律研究的難度。近年，隨著計算機技術和模式識別技術的快速發(fā)展，又出現(xiàn)了利用閃電定位和雷達回波數(shù)據(jù)識別、追蹤雷暴云，從而采用外推法進行雷暴運動趨勢臨近預報的方法，并取得了一定成果[7]。如R.E.Rinchert提出利用模式識別方法識別和匹配相鄰強回波單體；J.T.Johnson提出改進的SCIT算法[8]，使雷暴跟蹤的準確率提高到90%；文獻[9]提出了利用閃電定位數(shù)據(jù)基于DBScan聚類算法進行雷暴云識別從而進行其運動趨勢預測的方法，可實現(xiàn)預測未來15 min內(nèi)雷暴位置和覆蓋區(qū)域預測。

本文提出了一種新型的基于clique聚類識別[10]和卡爾曼濾波[11-12]進行雷電識別及外推的方法，并在雷電預警監(jiān)測系統(tǒng)中分別利用本方法和傳統(tǒng)的TITAN風暴模型[13]路徑預測方法實現(xiàn)了未來1小時逐6分鐘的雷暴路徑外推，基于2020年廣州閃電定位數(shù)據(jù)對兩種方法進行了有效性檢驗評估。

1 雷電識別外推算法設計

1.1 基于卡爾曼濾波的追蹤外推算法

卡爾曼濾波是一種最優(yōu)化自回歸數(shù)據(jù)處理算法，是一種針對線性系統(tǒng)的高效率狀態(tài)濾波技術，在對時間序列分析、軌跡最優(yōu)化等方面都有較好的效果[14-15]。它是一種遞歸的濾波過程，即只要獲知上一時刻狀態(tài)的估計值以及當前狀態(tài)的觀測值就可計算出當前狀態(tài)的估計值，因此不需要記錄觀測或者估計的歷史信息(如圖1)。

圖1 卡爾曼濾波過程

在雷暴預測應用中，可以首先利用Clique空間網(wǎng)格聚類算法對雷暴單體進行識別，然后利用卡爾曼濾波算法，對識別出的雷暴單體運動進行線性外推。其中每6分鐘設為一個時間段，在時間段結束時，采集融合該6分鐘內(nèi)的閃電數(shù)據(jù)，按圖2算法流程執(zhí)行。

圖2 卡爾曼濾波路徑外推流程圖

1.1.1 基于Clique聚類算法的雷暴識別

由于閃電基本上由雷暴云團生成，在空間上具有連續(xù)性，雷暴云團的移動意味著閃電位置的移動，因此在進行雷電外推時，可通過聚類算法，將在空間上具有一定連續(xù)性的離散閃電聚合成獨立的雷電單體，從而簡化雷電外推過程。此處選擇了運算速度較快，抗干擾性較好，可進行任意形狀聚合的Clique網(wǎng)格空間聚類算法[16-18]。其基本步驟如下。

步驟1：閃電坐標網(wǎng)格化，將6分鐘內(nèi)的所有閃電根據(jù)經(jīng)緯度坐標放入網(wǎng)格化的平面直角坐標系；

步驟2：根據(jù)設定好的閾值參數(shù)，掃描符合閃電頻數(shù)閾值和雷暴面積閾值的網(wǎng)格并合并；

步驟3：記錄符合閾值要求的網(wǎng)格為閃電簇區(qū)域；

步驟4：合并相鄰的閃電簇區(qū)域并編號；

步驟5：掃描所有閃電簇區(qū)域，利用橢圓法包絡閃電簇，返回識別的雷暴區(qū)域。

1.1.2 卡爾曼預測

把雷暴的移動過程看作是系統(tǒng)的變化，應用卡爾曼濾波中的時間更新方程對前一時間段的雷暴進行預測，從而得到現(xiàn)時間段的預測值。通過前一時間的系統(tǒng)狀態(tài)Xn-1和協(xié)方差矩陣Pn-1預測最新時間的系統(tǒng)狀態(tài)和協(xié)方差，并把預測結果作為預測值保存下來。該過程使用時間更新方程：

(1)

(2)

其中:F是狀態(tài)轉(zhuǎn)移矩陣，代表系統(tǒng)變化的方向和幅度，Q是狀態(tài)轉(zhuǎn)移噪聲矩陣，代表該過程可能產(chǎn)生的誤差變量。

1.1.3 雷暴追蹤

步驟1：獲取在前一時間段的N1個雷暴，在雷暴識別步驟中獲取現(xiàn)時間段識別的N2個雷暴，把兩個相鄰時間段的雷暴作為二分圖中的兩個集合A和B；

步驟2：計算A中每個雷暴i對應B中每個雷暴j的代價函數(shù)Cij，其中0

Cij=dp+ds

(3)

其中:dp是雷暴中心點的位置差；ds是雷暴面積的平方根的差。將Cij作為雷暴二分圖之間邊上的權值；

步驟3：針對帶權值的雷暴二分圖，利用二分圖最小權值匹配尋找到能使代價函數(shù)的和∑Cij的匹配方案；

步驟4：完成匹配后，把雷暴分為四個集合，分別為前一時間段的已匹配和未匹配雷暴集合Mn-1和UMm-1，以及現(xiàn)時間段的已匹配和未匹配雷暴集合Mn和UMn；

步驟5：結合卡爾曼預測所得的預測值，再分別在重合的基礎上匹配Mn-1和UMn、UMm-1和Mn，以此判斷是否產(chǎn)生了雷暴分裂和合并，并對分裂和合并的雷暴進行方向和速度上的校準；

步驟6：分析處理剩余雷暴，添加新生雷暴，舍棄追蹤失敗的雷暴，并把雷暴追蹤的結果作為觀察值。

1.1.4 卡爾曼更新

應用卡爾曼濾波中的狀態(tài)更新方程對雷暴狀態(tài)進行更新，從而得到現(xiàn)時間段的雷暴狀態(tài)值。首先利用預測過程得到的協(xié)方差預測值結合觀察矩陣計算卡爾曼增益Kn，再利用系統(tǒng)的預測值Xn-1和最新時間的觀察值Yn更新系統(tǒng)的最新狀態(tài)。該過程使用狀態(tài)更新方程：

(4)

(5)

(6)

其中:Xn是更新后的后驗狀態(tài)矩陣，Pn是更新后的后驗協(xié)方差矩陣，H為觀察模型，R是觀察噪聲矩陣。

1.1.5 外推預測路徑

根據(jù)更新后得到的雷暴狀態(tài)，通過線性外推得到雷暴預測路徑，以此外推未來1小時逐6分鐘(即6、12、……、60 min)的雷暴狀態(tài)，得到一條雷暴外推路徑；在每次時間段結束都可以進行迭代外推，從而不斷更新且精確化外推路徑。

1.2 基于TITAN風暴路徑的追蹤外推算法

基于TITAN風暴路徑進行雷電追蹤和外推則是目前應用較多的另一種雷電臨近預報方法。其主要過程也包括：閃電融合、聚類識別、追蹤外推[19-20]。

1.2.1 閃電融合

由于TITAN風暴移動路徑時間間隔和雷達回波探測間隔均為6 min，首先需要將每次采集到的零散閃電定位數(shù)據(jù)融合成逐6 min過程，與雷達的探測間隔相對應，融合過程如圖3所示。

圖3 閃電融合過程圖

當前時間采集到閃電，計算出起始和截止時間，將其對應06/12/18/24/30/36/42/48/54/00分鐘時段內(nèi)，按照時段的起止時間，從數(shù)據(jù)表中檢索出已存在的閃電進行融合，從而形成一個逐6 min的閃電集合，作為后續(xù)進行聚類識別和外推的雷電數(shù)據(jù)。

1.2.2 聚類識別

此處采用DBSCAN (density-based spatial clustering of applications with noise)密度聚類算法，用距離作為閃電密度的聚合因素，將多個離散的閃電聚合成若干個雷電單體，離散的閃電聚合成雷電單體后，利用雷電單體內(nèi)的閃電位置，加權計算出雷電單體的質(zhì)心位置，再以質(zhì)心位置與風暴移動路徑進行空間關聯(lián)和外推。

1.2.3 外推預測

利用成熟的TITAN雷暴識別、追蹤、分析和臨近預報算法，從前后兩個時次的雷達回波圖中識別出雷暴單體，利用兩個單體之間的位置偏移，計算出雷暴的移動方向和移動速度，從而得到風暴的移動路徑，將其作為雷電的移動路徑。由于TITAN數(shù)據(jù)接口提供的是未來1小時逐10 min的風暴移動路徑，還需要進行差值運算生成未來逐6 min的風暴移動路徑以匹配雷電預報路徑。

通常雷電單體會隨著風暴移動，但雷電單體和風暴位置往往并不重疊，同時由于閃電定位儀和雷達探測設備之間數(shù)據(jù)處理存在時間差，因此需利用時間和空間關聯(lián)方法，根據(jù)聚類出來的雷電單體的時間和質(zhì)心位置，匹配到最近時間和空間距離最近的風暴移動路徑，將風暴預報路徑平移后，作為雷電單體的預報路徑，對雷電單體進行逐6 min的位置外推(圖4中虛線部分)，生成未來1小時逐6 min的雷電預報產(chǎn)品。

圖4 移動路徑外推圖

2 雷電預警模型和平臺設計

基于上述兩種雷電識別外推方法和廣州范圍內(nèi)閃電定位儀、大氣電場儀、雷達等雷電實況監(jiān)測數(shù)據(jù)，廣州市氣象局開發(fā)了雷電監(jiān)測預警系統(tǒng)，可在GIS平臺上有效提供當前雷電實況數(shù)據(jù)和未來1小時雷電預報產(chǎn)品，并對兩種識別外推方法進行有效性檢驗評估，可廣泛應用于防雷減災、雷災調(diào)查、預警預報服務等領域。其整體方案如圖5所示。

圖5 系統(tǒng)結構圖

1)雷電實況監(jiān)測資料：包括EN全閃、粵港澳閃電、風暴路徑和雷達回波及廣州新建大氣電場儀等。

2)雷電預警模型：包括采集、融合、聚類分析、外推、預警和網(wǎng)格化等過程。①采集：從廣東省氣象探測數(shù)據(jù)中心的通用氣象數(shù)據(jù)訪問接口采集雷達回波、EN全閃、粵港澳閃電、風暴路徑等雷電實況資料，入庫存儲；②融合：將逐分鐘的閃電定位融合成逐6 min的數(shù)據(jù)，與雷達觀測時間進行匹配；③聚類分析：將離散的閃電定位聚合成雷電單體；④外推：利用風暴移動路徑，對雷電單體進行未來1小時的外推；⑤預警和網(wǎng)格化：將地圖劃分成1 * 1 km分辨率的網(wǎng)格，設置每個格點的雷電紅色和黃色預警標準，利用未來1小時逐6 min的雷電單體位置，按照距離權重法，插值成未來1小時逐6分鐘的雷電精細化格點預警產(chǎn)品。

3)雷電數(shù)據(jù)庫：采用關系型數(shù)據(jù)庫和NetCDF文件庫相結合的方式來存儲雷電實況監(jiān)測和預警數(shù)據(jù)。其中關系型數(shù)據(jù)庫用于保存閃電和風暴路徑等結構化雷電數(shù)據(jù)；NetCDF文件庫用于存儲雷達回波和雷電精細化格點預警產(chǎn)品等網(wǎng)格數(shù)據(jù)。

4)雷電預警模型檢驗系統(tǒng)：在WebGIS地理信息系統(tǒng)上，實現(xiàn)雷達回波、閃電、風暴移動路徑和雷電精細化格點預警產(chǎn)品的綜合顯示，通過對雷電數(shù)據(jù)的可視化分析，檢驗雷電預警模型的輸出結果。

系統(tǒng)界面如圖6，其中黑色路徑為基于TITAN風暴路徑生成的未來1小時雷電路徑預報，灰色路徑為基于卡爾曼濾波算法生成的未來1小時雷電路徑預報。

圖6 系統(tǒng)界面

3 實驗結果與分析

3.1 數(shù)據(jù)來源

為了對兩種外推算法結果的有效性進行檢驗評估，本文使用了2020年5月1日至2020年10月27日廣東省氣象局數(shù)據(jù)接口所提供的廣東EN全閃觀測網(wǎng)的全省閃電定位數(shù)據(jù)作為檢驗結果的數(shù)據(jù)來源。

3.2 檢驗方法

定義命中率和提前量為兩個檢驗指標。

其中命中率是針對每個預報路徑點，讀取對應時刻發(fā)生的閃電，如閃電定位數(shù)據(jù)顯示半徑5 km范圍內(nèi)有地閃，則認為命中，否則認為空報，即：

(7)

其中:N命中為命中次數(shù)，N空報為空報次數(shù)。

提前量則是針對每個預報路徑點，其對應時刻的未來1小時內(nèi)發(fā)生閃電，則認為命中，計算第一次閃電和當前預報之間的時間差，作為預報的提前量。

統(tǒng)計時，以間隔6分鐘為一個時間片，從起始時間至結束時間整個時間段內(nèi)逐時間片計算每個預測路徑的命中率和提前量。

3.3 結果分析

統(tǒng)計后的檢驗結果如圖7可見，兩種外推算法的命中率均隨預測時間的增加而衰減，越是臨近當前時刻，預測命中率越高，未來6分鐘預測命中率最高可達到87.5%，未來60分鐘預測命中率最高27.2%。且由于2020年8月份對基于卡爾曼濾波的外推算法進行了參數(shù)優(yōu)化，將聚類識別過程中的網(wǎng)格密度增大后，其預測命中率在2020年9月至10月的檢驗中已明顯優(yōu)于基于TITAN路徑的外推算法。而兩種外推算法在各時次的預測提前量上差別不大，平均相差1分鐘左右，整體表現(xiàn)比較接近。但在命中次數(shù)上，基于TITAN路徑的外推算法要明顯高于基于卡爾曼濾波的外推算法，其主要原因還是因為后者所采用的Clique聚類算法對分散的較小型雷暴單體識別率低于前者采用的DBSCAN密度聚類算法，需進一步進行優(yōu)化改進。

圖7 檢驗結果對比

4 結束語

本文提出了一種基于clique聚類識別和卡爾曼濾波算法進行雷電識別及外推的方法，并利用該方法和傳統(tǒng)的基于TITAN風暴路徑外推的兩種算法，研究開發(fā)了廣州雷電監(jiān)測預警系統(tǒng)，該系統(tǒng)接入實時雷電氣象數(shù)據(jù)，可生成未來1小時逐6分鐘1*1 km分辨率的雷電精細化格點預警產(chǎn)品，并對預報產(chǎn)品進行檢驗，可為雷電預警業(yè)務化提供有效平臺。實際數(shù)據(jù)檢驗表明，兩種算法均能有效識別追蹤雷暴并進行未來1小時路徑預測，經(jīng)過參數(shù)優(yōu)化后，基于卡爾曼濾波的算法在預測命中率上已優(yōu)于傳統(tǒng)的基于TITAN風暴路徑的外推算法，下一步還需要在聚類識別算法上進行研究優(yōu)化，進一步提高雷電預警的準確性、及時性。