解妍瓊 張云 楊波 王佳 王亞 張鵬

(1 國防科技大學(xué) 氣象海洋學(xué)院，南京 211101；2 中國氣象局 云霧物理環(huán)境重點實驗室；北京 100081；3 陸軍工程大學(xué) 教研保障中心，南京 210007)

引 言

天氣雷達(dá)雷暴自動識別與跟蹤是自動天氣分析和雷電臨近預(yù)報中的一個基本問題。雷暴一經(jīng)識別，其空間位置、移動路徑等發(fā)展特征可進(jìn)一步確定。利用天氣雷達(dá)回波數(shù)據(jù)進(jìn)行雷暴(或風(fēng)暴)識別、跟蹤和臨近預(yù)報的研究已經(jīng)進(jìn)行了幾十年[1-5]。天氣雷達(dá)反射率因子是雷暴識別中最常用的參數(shù)，識別風(fēng)暴的最簡單方法是設(shè)定一個合理的雷達(dá)反射率因子對圖像進(jìn)行閾值處理，由于該方法簡單實用，在許多風(fēng)暴識別方法中得到了應(yīng)用。兩種著名的方法是TITAN(雷暴識別、跟蹤、分析和臨近預(yù)報)[2]和SCIT(風(fēng)暴單體識別和跟蹤)[3]。還有很多關(guān)于雷暴雷達(dá)回波特征的研究[6-10]，主要以雷達(dá)的反射率因子為研究對象。如Trapp， et al[4]研究了1998—2000年發(fā)生的3 828次龍卷的雷達(dá)反射率因子圖像，以確定準(zhǔn)線性對流系統(tǒng)引發(fā)的龍卷風(fēng)的百分比。

在對天氣雷達(dá)回波雷暴特性的研究中，許多研究給出了與不同溫度層和反射率因子相關(guān)的統(tǒng)計結(jié)果[11-16]。Martinez， et al[17]發(fā)現(xiàn)，當(dāng)強(qiáng)度超過40 dBZ的回波對應(yīng)的高度大于7 km時才會發(fā)生閃電。LIU， et al[18]利用TRMM(Tropical Rainfall Mea￣suring Mission)1998—2010年共13 a的觀測資料研究了雷達(dá)反射率因子垂直剖面與閃電頻度(Lightning Flash Rates, LFR)之間的關(guān)系，檢驗了20、30和40 dBZ回波頂溫度與閃電活動的相關(guān)性。結(jié)果表明，回波頂溫度與閃電活動的相關(guān)性并不高，而混合相態(tài)的強(qiáng)回波區(qū)大小與閃電發(fā)生率及閃電發(fā)生區(qū)域有更高的相關(guān)性。文獻(xiàn)[19—21]都將-10 ℃高度層對應(yīng)的反射率因子值達(dá)到40 dBZ作為地閃活動的預(yù)報因子。

YANG， et al[22]提出了一種基于雙偏振雷達(dá)資料識別雷暴的算法，利用雙偏振雷達(dá)資料進(jìn)行粒子類型識別，然后綜合識別出的霰粒子面積與雷達(dá)反射率因子進(jìn)行雷暴識別。結(jié)果表明，在CAPPI 0 ℃高度層強(qiáng)回波區(qū)域反射率因子識別閾值為30 dBZ，霰粒子區(qū)篩選閾值為2 km2時，得到了最優(yōu)的雷暴識別結(jié)果，識別準(zhǔn)確度達(dá)到91%以上，識別虛警率(False Alarm Rate, FAR)僅為6.9%，臨界成功指數(shù)為85.3%，這一方法能夠達(dá)到很好的雷暴識別效果，但必須基于雙偏振雷達(dá)。雖然我國業(yè)務(wù)雷達(dá)正在進(jìn)行雙偏振的升級改造，但相當(dāng)長一段時間內(nèi)，單偏振雷達(dá)仍是業(yè)務(wù)主體。如何使單偏振天氣雷達(dá)具備相應(yīng)的雷暴自動識別能力一直是氣象業(yè)務(wù)的熱點和難點問題，本文在前期研究基礎(chǔ)上建立一種基于溫度層強(qiáng)回波區(qū)域面積的天氣雷達(dá)雷暴識別方法，解決單偏振天氣雷達(dá)雷暴自動識別問題，并與文獻(xiàn)[22]提出的算法進(jìn)行比較，分析本文算法的識別效果，為提高單偏振天氣雷達(dá)雷暴自動分析能力提供技術(shù)手段。

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

本章利用S波段雙偏振多普勒天氣雷達(dá)2014年和2015年雷暴探測實驗觀測資料結(jié)合江蘇省閃電定位資料、地面大氣電場資料和國家氣候基準(zhǔn)站南京站點的探空數(shù)據(jù)研究雷暴的雷達(dá)回波特征。觀測使用的S波段雙偏振多普勒天氣雷達(dá)，性能指標(biāo)可參考文獻(xiàn)[22-23]，不再贅述。江蘇省閃電定位系統(tǒng)為ADTD型雷電定位系統(tǒng)，該系統(tǒng)綜合使用磁定向法以及多站時差定位法探測地閃，定位精度約為500 m，探測效率約為90%～95%。

表1 樣本資料簡要說明Table 1 The brief description of sample data

為研究雷暴過程的雙偏振雷達(dá)回波特征，選取2014—2015年共17個雷暴天氣過程進(jìn)行分析，具體樣本日期見表1。為保證數(shù)據(jù)的有效性，選取大氣電場儀或閃電定位系統(tǒng)監(jiān)測區(qū)域內(nèi)的雷達(dá)回波進(jìn)行分析，采用了地面大氣電場資料、江蘇省閃電定位資料對對流單體進(jìn)行篩選，區(qū)分雷暴單體與非雷暴單體。

2 雷暴識別算法

由于雷暴是一種中小尺度強(qiáng)對流天氣系統(tǒng)，在水平范圍內(nèi)具有較強(qiáng)的雷達(dá)回波反射率因子，對流高度能夠達(dá)到幾公里甚至十幾公里，且對流越旺盛越有利于雷暴起電，因此本文雷暴識別算法的基本思想是利用雷達(dá)回波體掃數(shù)據(jù)首先在水平范圍內(nèi)搜索強(qiáng)回波區(qū)域；再根據(jù)某一溫度層上，達(dá)到設(shè)定閾值的雷達(dá)回波反射率因子積分面積，判斷水平范圍強(qiáng)回波區(qū)域是否符合雷暴特征；當(dāng)強(qiáng)回波區(qū)域符合雷暴特征時，判定該強(qiáng)回波區(qū)域為雷暴區(qū)域。根據(jù)雷暴的雷達(dá)回波特征，采用不同溫度層的雷達(dá)回波反射率因子積分面積實現(xiàn)雷暴識別。

2.1 算法總體流程

識別算法總體流程如圖1所示，主要包括：

(1) 將雷達(dá)體掃數(shù)據(jù)作為輸入量，完成數(shù)據(jù)解析；

(2) 利用體掃回波數(shù)據(jù)各層的反射率因子生成組合反射率數(shù)據(jù)；

(3) 通過設(shè)定的反射率因子閾值Zth1，利用組合反射率數(shù)據(jù)進(jìn)行雷達(dá)回波強(qiáng)中心搜索；

(4) 通過設(shè)定的區(qū)域尺度閾值A(chǔ)1，對組合反射率數(shù)據(jù)雷達(dá)回波強(qiáng)中心進(jìn)行區(qū)域尺度篩選，篩選出區(qū)域尺度大于區(qū)域尺度閾值A(chǔ)1的回波區(qū)域，篩選結(jié)果構(gòu)成待定雷暴區(qū)域集合C1；

(5) 利用步驟(1)解析得到的體掃數(shù)據(jù)，生成各溫度層CAPPI反射率因子數(shù)據(jù)；

(6) 利用某溫度層CAPPI反射率因子閾值Zth2進(jìn)行雷達(dá)回波強(qiáng)中心搜索；

(7) 通過設(shè)定的區(qū)域尺度閾值A(chǔ)2，對CAPPI反射率因子雷達(dá)回波強(qiáng)中心進(jìn)行區(qū)域尺度篩選，篩選出區(qū)域尺度大于區(qū)域尺度閾值A(chǔ)2的回波區(qū)域，構(gòu)成待定雷暴區(qū)域集合C2；

(8) 將第4步得到的待定雷暴區(qū)域集合C1與第7步得到的待定雷暴區(qū)域集合C2進(jìn)行匹配，當(dāng)集合C2中某區(qū)域中心的水平投影落在集合C1中某區(qū)域的范圍內(nèi)時，判斷為匹配成功，將C1對應(yīng)的待定雷暴區(qū)域確定為雷暴單體。

其中，步驟(5)中的CAPPI的高度根據(jù)探空數(shù)據(jù)的溫度確定，本文主要采用了0 ℃、-10 ℃和-15 ℃對應(yīng)的高度，由于探空數(shù)據(jù)一般情況不會出現(xiàn)0 ℃、-10 ℃和-15 ℃這幾個高度，因此，溫度層對應(yīng)的高度在實際處理過程中采用插值方法得到。

圖1 雷暴識別算法流程Fig.1 The overall process of thunderstorm identification algorithm

2.2 雷暴強(qiáng)中心提取

(1)強(qiáng)回波點提取

圖2 強(qiáng)回波點提取效果：(a)原始雷達(dá)回波；(b)強(qiáng)回波點提取結(jié)果Fig.2 Effect of intense echo point extraction: (a) the original radar echo image;(b) the intense echo point extraction results

設(shè)組合反射率雷達(dá)回波圖和CAPPI反射率因子回波圖為H×V像素的圖像，則圖像共包含H×V個像素點，由ZH(H，V)表示各像素點的雷達(dá)回波反射率因子值，當(dāng)ZH(i，j)>ZTh時，ZH(i，j) 即為強(qiáng)回波點，其中，i為回波圖橫坐標(biāo)，j為回波圖縱坐標(biāo)，0

(2)強(qiáng)回波段合成

在橫向或縱向上，逐點對上一步驟得到的強(qiáng)回波點集合ZHC進(jìn)行相鄰點合并，構(gòu)成橫向或縱向強(qiáng)回波段集合，表示為ZHS(N)，N為強(qiáng)回波段集合大小。以縱向強(qiáng)回波段為例，強(qiáng)回波段ZHS(x)是由多個橫坐標(biāo)相同，在縱向坐標(biāo)位置連續(xù)的點構(gòu)成的一維回波數(shù)據(jù)點集合。當(dāng)某強(qiáng)回波點ZHC(i，j)為孤立點時，則判定該回波點為無效點，直接刪除。

(3)強(qiáng)回波區(qū)域合成

對強(qiáng)回波段合成后得到的強(qiáng)回波段集合ZHS(N)進(jìn)行相鄰強(qiáng)回波段合成，即可構(gòu)成強(qiáng)回波區(qū)域集合ZHA(M)，M為強(qiáng)回波區(qū)域集合大小。當(dāng)某強(qiáng)回波段ZHS(x)為孤立強(qiáng)回波段時，則判定該回波段為無效回波段，直接刪除。

強(qiáng)回波區(qū)域合成后，可利用邊界點識別算法提取各強(qiáng)回波區(qū)域的邊界點，以便于CAPPI識別出的強(qiáng)回波區(qū)域和組合反射率識別出的強(qiáng)回波區(qū)域進(jìn)行水平投影位置匹配。本文采用遍歷各強(qiáng)回波點的算法提取邊界點，如圖3所示。紅色區(qū)域表示強(qiáng)回波區(qū)域，每一個格點代表一個回波點，遍歷各格點數(shù)據(jù)，當(dāng)某數(shù)據(jù)點的相臨點個數(shù)小于4時，該數(shù)據(jù)點為邊界點，圖中A、B兩點的相鄰點的個數(shù)為4，判定為非邊界點，除這兩點外的各數(shù)據(jù)點的相鄰點個數(shù)都小于4，判定為邊界點。采用該算法即可提取出強(qiáng)回波區(qū)域的邊界。

圖3 強(qiáng)回波區(qū)域邊界點提取算法示意Fig.3 Schematic diagram of boundary point extractionalgorithm in intense echo region

圖4 強(qiáng)回波區(qū)域邊界提取效果Fig.4 Effect of intense echo region boundary extraction

根據(jù)上述算法，采用合適的參數(shù)即可提取出強(qiáng)回波區(qū)域，提取效果如圖4所示。強(qiáng)回波區(qū)域包含強(qiáng)回波段，強(qiáng)回波段包含強(qiáng)回波點，因此強(qiáng)回波區(qū)域包含區(qū)域內(nèi)各強(qiáng)回波點的坐標(biāo)信息和反射率因子信息。各強(qiáng)回波區(qū)域的邊界可以從回波反射率因子圖中提取出來。實際處理算法中可對提取出來的每個強(qiáng)回波區(qū)域進(jìn)行編號，再通過遍歷的方式對這些強(qiáng)回波區(qū)域進(jìn)行區(qū)域尺度篩選和位置匹配。

2.3 區(qū)域尺度篩選

由于雷暴單體的尺度一般在幾公里到上百公里，對應(yīng)尺度過小的區(qū)域，可通過區(qū)域尺度篩選將其判斷為非雷暴區(qū)域。雷暴識別算法總體流程中，組合反射率回波圖和CAPPI回波圖中都采用了區(qū)域尺度篩選，對尺度過小的區(qū)域進(jìn)行篩選。在組合反射率回波圖中以區(qū)域尺度閾值A(chǔ)1，CAPPI回波圖中以區(qū)域尺度閾值A(chǔ)2進(jìn)行篩選，當(dāng)強(qiáng)回波區(qū)域的面積小于對應(yīng)的區(qū)域尺度閾值時，強(qiáng)回波區(qū)域集合中刪除該區(qū)域。將組合反射率篩選得到的待定雷暴區(qū)域集合C1，CAPPI在某高度層上篩選得到的待定雷暴區(qū)域集合C2。

2.4 篩選結(jié)果匹配

最后將待定雷暴區(qū)域集合C1和C2進(jìn)行水平投影位置匹配，若匹配成功，即判定對應(yīng)的組合反射率區(qū)域為雷暴區(qū)域。水平投影匹配基本思路是首先利用2.2節(jié)中提取的各強(qiáng)回波區(qū)域邊界，求取CAPPI各強(qiáng)回波區(qū)域的幾何中心點Pt，當(dāng)Pt的水平投影落在組合反射率的某個待定雷暴區(qū)域時，則判定匹配成功。

3 結(jié)果檢驗與分析

為檢驗上述雷暴識別算法的有效性，利用表1所示天氣過程的雷達(dá)回波數(shù)據(jù)和對應(yīng)的閃電定位數(shù)據(jù)進(jìn)行驗證，為便于雷暴識別準(zhǔn)確率統(tǒng)計，將閃電定位數(shù)據(jù)通過經(jīng)緯度映射到雷達(dá)回波圖像中。

根據(jù)上述雷暴識別算法流程，需要在組合反射率和CAPPI中分布搜索反射率因子強(qiáng)回波中心。由于本算法雷暴識別的關(guān)鍵在于CAPPI不同高度層上反射率因子閾值參數(shù)選擇，為減小統(tǒng)計工作量，將組合反射率強(qiáng)回波中心搜索閾值設(shè)為固定值40 dBZ，區(qū)域尺度閾值A(chǔ)1設(shè)為固定值1 km2。在此基礎(chǔ)上，識別過程中還需設(shè)置CAPPI圖像的強(qiáng)回波中心搜索反射率因子閾值Z2，溫度層HT以及反射率因子面積閾值A(chǔ)2，構(gòu)成識別算法參數(shù)組合(Z2，HT，A2)。

在溫度層HT的處理上本文算法做了近似，由于雷暴是強(qiáng)對流系統(tǒng)，其內(nèi)部存在強(qiáng)烈的上升氣流和下沉氣流，云內(nèi)某一溫度層不可能在固定的高度，考慮到環(huán)境溫度是雷暴起電的重要影響因素，在大量相關(guān)的研究中都采用探空數(shù)據(jù)給出的溫度層(固定高度層)作為宏觀特征量對雷暴進(jìn)行研究，因此本文算法也將溫度層近似認(rèn)為是一個固定的高度，利用最近的探空站點最接近的時間點探空數(shù)據(jù)判定溫度層的高度，然后在CAPPI中得到不同高度的反射率因子圖像。

為檢驗各種參數(shù)的識別效果，采用探測概率 (Probability of Detection，POD)、漏警率 (Miss Rate，MR)、虛警率 (False Alarm Rate，F(xiàn)AR)和臨界成功指數(shù) (Critical Success Index，CSI)指數(shù)來衡量檢測結(jié)果，分別表示為：

(1)

(2)

(3)

(4)

其中：R1，R2和R3分別為識別準(zhǔn)確數(shù)、漏警數(shù)和虛警數(shù)。

檢驗過程中，利用上述識別雷暴算法結(jié)合識別參數(shù)組合，生成識別結(jié)果回波圖像，再將雷達(dá)回波數(shù)據(jù)對應(yīng)時段的閃電定位結(jié)果疊加到識別結(jié)果回波圖像上，如圖5所示。其中橢圓內(nèi)的強(qiáng)回波區(qū)為通過上述識別算法自動計算得到的雷暴區(qū)域；當(dāng)閃電定位結(jié)果落在雷暴識別的區(qū)域時，判定雷暴識別結(jié)果準(zhǔn)確，相應(yīng)R1增加1；當(dāng)某區(qū)域有閃電發(fā)生，但未識別為雷暴區(qū)域時，判定為漏警，相應(yīng)R2增加1；當(dāng)某區(qū)域被識別為雷暴區(qū)域，但該區(qū)域沒有閃電發(fā)生時，判定為虛警，相應(yīng)R3增加1。

現(xiàn)有的研究表明，將不同的溫度層的反射率因子用于雷暴預(yù)測時，會得到不同的預(yù)測效果，且多項研究表明在-10 ℃高度層回波強(qiáng)度達(dá)到40 dBZ以上，能夠較好地區(qū)分雷暴和非雷暴單體[17-19]。文獻(xiàn)[21]研究了0、-10和-15 ℃高度層霰粒子分布對雷暴識別結(jié)果的影響，為對比有無雙偏振信息對結(jié)果的影響，本文對應(yīng)地也采用這三個高度層進(jìn)行研究。

圖5 雷達(dá)回波雷暴識別結(jié)果檢驗(“-”代表負(fù)閃電，“+”代表正閃電)：(a) 系統(tǒng)性雷暴識別結(jié)果；(b)孤立單體雷暴識別結(jié)果；(c)颮線雷暴識別結(jié)果；(d)混合性雷暴識別結(jié)果Fig.5 Radar echo thunderstorm identification result test chart( “-” represents negative lightning and “+” represents positive lightning): (a) identification results of systematic thunderstorm; (b) identification results of isolated thunderstorm cell;(c) identification results of squall line thunderstorm;(d) identification results of mixed thunderstorm

檢驗過程中，首先考慮Z2為30 dBZ時的參數(shù)組合，研究HT為0、-10和-15 ℃層高度上區(qū)域尺度閾值A(chǔ)2分別是1、2和3 km2時的情況，此時，共有9種參數(shù)組合，利用這9種參數(shù)組合對表1所示天氣過程的雷達(dá)回波數(shù)據(jù)進(jìn)行雷暴識別，識別統(tǒng)計結(jié)果如表2所示，統(tǒng)計的有閃電數(shù)據(jù)對應(yīng)的雷暴單體共有312個。其中，識別結(jié)果較好的參數(shù)組合是HT為0 ℃高度的3種參數(shù)組合，第1種參數(shù)組合(HT為0 ℃高度，A2=1 km2)能夠得到最高的探測概率，POD達(dá)到89.7%，但虛警率略高，F(xiàn)AR為41.1%，臨界成功指數(shù)為55.2%；從CSI看，第5種參數(shù)組合(HT為-10 ℃，A2=2 km2)能夠得到最佳的識別效果，POD雖然低于第一種參數(shù)組合，僅達(dá)到了71.8%，但虛警率相對較低，F(xiàn)AR為20.6%，CSI為60.5%。此外，HT為-10 ℃時的三種組合CSI都優(yōu)于0 ℃和-15 ℃的幾種組合。

再研究Z2為35 dBZ時的參數(shù)組合，與Z2為30 dBZ時的分析方法相同，HT設(shè)為0、-10和-15 ℃，A2設(shè)為1、2和3 km2，共9種參數(shù)組合，識別統(tǒng)計結(jié)果如表3所示，可發(fā)現(xiàn)在0 ℃高度層上的識別效果最好，(35 dBZ， 0 ℃， 2 km2)參數(shù)組合能夠獲得最好的識別效果，POD達(dá)到87.5%，虛警率為32.9%，臨界成功指數(shù)為61.2%。0 ℃高度層上的3種參數(shù)組合POD和CSI相對于其余兩種高度層的參數(shù)組合，有較大的優(yōu)勢，但FAR則相對較高。

最后考慮ZTh為40 dBZ時的參數(shù)組合，采用相同的(HT，A2)參數(shù)組合進(jìn)行檢驗，識別統(tǒng)計結(jié)果如表4所示，可以發(fā)現(xiàn)對應(yīng)的9種參數(shù)組合都僅能得到很低的POD，最高僅達(dá)到了43.2%，CSI最高僅為42.2%。

表2 Z2=30 dBZ時的雷暴識別結(jié)果Table 2 Thunderstorm identification results (Z2=30 dBZ)

表3 Z2=35 dBZ時的雷暴識別結(jié)果Table 3 Thunderstorm identification results (Z2=35 dBZ)

由以上統(tǒng)計結(jié)果，在某溫度層上，對反射率因子大于閾值的區(qū)域進(jìn)行面積求和，將總面積作為雷暴識別的基本指標(biāo)合理且能夠得到較好的識別結(jié)果。又由于不同雷暴的對流發(fā)展不同，導(dǎo)致不同雷暴在不同的溫度層的超過反射率因子閾值的面積不同，對流越旺盛，在較高高度上的反射率因子越大，越符合雷暴的特征，雷暴的閃電活動也越劇烈。在將-10 ℃或-15 ℃高度層的反射率因子閾值作為識別參數(shù)時，盡管識別POD相對較低，但識別到的雷暴區(qū)域是所有雷暴區(qū)域中閃電活動相對較為劇烈的區(qū)域，因此，以強(qiáng)雷暴篩選為目的時，可以選擇相應(yīng)的參數(shù)組合。從CSI結(jié)果來看，27種參數(shù)組合中，(35 dBZ， 0 ℃， 2 km2)參數(shù)組合能夠達(dá)到最優(yōu)的識別效果。

4 結(jié)論

本文提出了基于溫度層強(qiáng)回波區(qū)域面積的天氣雷達(dá)雷暴自動識別方法，該方法將溫度層、反射率因子閾值和強(qiáng)回波區(qū)域面積相結(jié)合，包括組合反射率強(qiáng)反射率因子區(qū)域篩選、尺度篩選，某溫度層上CAPPI強(qiáng)反射率因子區(qū)域篩選，組合反射率識別區(qū)域和CAPPI識別區(qū)域匹配幾個步驟，該方法中的關(guān)鍵在于選擇恰當(dāng)?shù)腃APPI篩選溫度層、強(qiáng)回波區(qū)反射率因子閾值和強(qiáng)回波區(qū)域面積尺度閾值。研究中，采用了17個雷暴天氣過程的312個雷暴區(qū)域檢驗識別算法，結(jié)果表明：當(dāng) CAPPI數(shù)據(jù)處理中采用(35 dBZ， 0 ℃， 2 km2)參數(shù)組合能夠達(dá)到最優(yōu)的識別效果，識別概率可達(dá)87.5%，虛警率為32.9%，臨界成功指數(shù)CSI為61.2%。該方法可以用于天氣雷達(dá)的雷暴計算機(jī)自動識別業(yè)務(wù)。此外，雖然對比文獻(xiàn)[22]將雙偏振參數(shù)用于雷暴識別，本文的雷暴識別能力要弱一些，但本方法實現(xiàn)了只有單偏振雷達(dá)時的雷暴自動識別，為目前氣象業(yè)務(wù)提供重要指導(dǎo)。