基于溫度層強回波區(qū)域面積的天氣雷達雷暴自動識別方法

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

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

引 言

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

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

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

1 數據來源

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

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

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

2 雷暴識別算法

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

2.1 算法總體流程

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

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

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

(3) 通過設定的反射率因子閾值Zth1，利用組合反射率數據進行雷達回波強中心搜索；

(4) 通過設定的區(qū)域尺度閾值A1，對組合反射率數據雷達回波強中心進行區(qū)域尺度篩選，篩選出區(qū)域尺度大于區(qū)域尺度閾值A1的回波區(qū)域，篩選結果構成待定雷暴區(qū)域集合C1；

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

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

(7) 通過設定的區(qū)域尺度閾值A2，對CAPPI反射率因子雷達回波強中心進行區(qū)域尺度篩選，篩選出區(qū)域尺度大于區(qū)域尺度閾值A2的回波區(qū)域，構成待定雷暴區(qū)域集合C2；

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

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

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

2.2 雷暴強中心提取

(1)強回波點提取

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

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

(2)強回波段合成

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

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

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

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

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

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

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

2.3 區(qū)域尺度篩選

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

2.4 篩選結果匹配

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

3 結果檢驗與分析

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

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

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

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

(1)

(2)

(3)

(4)

其中：R1，R2和R3分別為識別準確數、漏警數和虛警數。

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

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

圖5 雷達回波雷暴識別結果檢驗(“-”代表負閃電，“+”代表正閃電)：(a) 系統(tǒng)性雷暴識別結果；(b)孤立單體雷暴識別結果；(c)颮線雷暴識別結果；(d)混合性雷暴識別結果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時的參數組合，研究HT為0、-10和-15 ℃層高度上區(qū)域尺度閾值A2分別是1、2和3 km2時的情況，此時，共有9種參數組合，利用這9種參數組合對表1所示天氣過程的雷達回波數據進行雷暴識別，識別統(tǒng)計結果如表2所示，統(tǒng)計的有閃電數據對應的雷暴單體共有312個。其中，識別結果較好的參數組合是HT為0 ℃高度的3種參數組合，第1種參數組合(HT為0 ℃高度，A2=1 km2)能夠得到最高的探測概率，POD達到89.7%，但虛警率略高，FAR為41.1%，臨界成功指數為55.2%；從CSI看，第5種參數組合(HT為-10 ℃，A2=2 km2)能夠得到最佳的識別效果，POD雖然低于第一種參數組合，僅達到了71.8%，但虛警率相對較低，FAR為20.6%，CSI為60.5%。此外，HT為-10 ℃時的三種組合CSI都優(yōu)于0 ℃和-15 ℃的幾種組合。

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

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

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

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

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

4 結論

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