基于卡爾曼濾波和TITAN路徑的雷電預警

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

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

0 引言

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

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

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

1 雷電識別外推算法設計

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

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

圖1 卡爾曼濾波過程

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

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

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

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

步驟1：閃電坐標網格化，將6分鐘內的所有閃電根據經緯度坐標放入網格化的平面直角坐標系；

步驟2：根據設定好的閾值參數，掃描符合閃電頻數閾值和雷暴面積閾值的網格并合并；

步驟3：記錄符合閾值要求的網格為閃電簇區域；

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

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

1.1.2 卡爾曼預測

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

(1)

(2)

其中:F是狀態轉移矩陣，代表系統變化的方向和幅度，Q是狀態轉移噪聲矩陣，代表該過程可能產生的誤差變量。

1.1.3 雷暴追蹤

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

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

Cij=dp+ds

(3)

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

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

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

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

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

1.1.4 卡爾曼更新

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

(4)

(5)

(6)

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

1.1.5 外推預測路徑

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

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

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

1.2.1 閃電融合

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

圖3 閃電融合過程圖

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

1.2.2 聚類識別

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

1.2.3 外推預測

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

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

圖4 移動路徑外推圖

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

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

圖5 系統結構圖

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

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

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

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

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

圖6 系統界面

3 實驗結果與分析

3.1 數據來源

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

3.2 檢驗方法

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

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

(7)

其中:N命中為命中次數，N空報為空報次數。

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

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

3.3 結果分析

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

圖7 檢驗結果對比

4 結束語

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