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

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

0 引言

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

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

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

1 雷電識別外推算法設(shè)計

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

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

圖1 卡爾曼濾波過程

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

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

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

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

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

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

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

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

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

1.1.2 卡爾曼預(yù)測

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

(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對應(yīng)B中每個雷暴j的代價函數(shù)Cij，其中0

Cij=dp+ds

(3)

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

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

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

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

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

1.1.4 卡爾曼更新

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

(4)

(5)

(6)

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

1.1.5 外推預(yù)測路徑

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

1.2 基于TITAN風(fēng)暴路徑的追蹤外推算法

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

1.2.1 閃電融合

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

圖3 閃電融合過程圖

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

1.2.2 聚類識別

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

1.2.3 外推預(yù)測

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

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

圖4 移動路徑外推圖

2 雷電預(yù)警模型和平臺設(shè)計

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

圖5 系統(tǒng)結(jié)構(gòu)圖

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

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

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

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

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

圖6 系統(tǒng)界面

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

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

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

3.2 檢驗(yàn)方法

定義命中率和提前量為兩個檢驗(yàn)指標(biāo)。

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

(7)

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

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

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

3.3 結(jié)果分析

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

圖7 檢驗(yàn)結(jié)果對比

4 結(jié)束語

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