長沙機場陣列天氣雷達地物識別算法

魏萬益 馬舒慶 楊 玲* 甄小瓊4) 呂寺煒

1)(成都信息工程大學電子工程學院， 成都 610225)2)(中國氣象局大氣探測重點開放實驗室， 成都 610225)3)(中國氣象局氣象探測中心， 北京 100081)4)(中國科學院大氣物理研究所， 北京 100029)5)(雷象科技(北京)有限公司， 北京 100089)

引 言

地物雜波識別和去除是天氣雷達數(shù)據(jù)質(zhì)量控制研究中的一項基本工作。地物會影響雷達回波，掩蓋天氣回波，甚至直接影響雷達產(chǎn)品的可信度。

為了消除地物雜波，國內(nèi)天氣雷達專家們提出了多種去除地物雜波的方法，這些方法可分為以下3類：①通過改變雷達的設計和布局位置減少地物雜波，②在信號處理端對雷達信號進行處理，③針對雷達基數(shù)據(jù)進行后期處理。Smith[1]和Mann等[2]通過選擇雷達的安放位置和采用更短的波長減少雜波及其相對強度。第2類方法通常是對收集到的所有數(shù)據(jù)采用時域[2-4]或頻域[5-7]濾波器去除速度接近0的回波數(shù)據(jù)。Li等[8]在頻譜上提取特征參數(shù)作為樸素貝葉斯算法的輸入，以判斷是否存在地物雜波并加以去除。針對雷達基數(shù)據(jù)的處理方法主要是從基數(shù)據(jù)中提取二維或三維地物識別特征參數(shù)。Steiner等[9]使用回波強度三維結(jié)構作為決策樹算法輸入?yún)?shù)構建自動化雜波抑制算法，發(fā)現(xiàn)最有效的3個參數(shù)是雷達回波的垂直范圍、回波強度的水平變化及回波強度的垂直梯度。Zhang等[10]使用類似的回波強度垂直梯度去除地物雜波，并在計算該參數(shù)時考慮仰角抬高對水平距離的影響。Lakshmanan等[11-12]采用神經(jīng)網(wǎng)絡方法，將回波強度的三維結(jié)構、速度和譜寬以及其他特征作為神經(jīng)網(wǎng)絡的輸入數(shù)據(jù)去除地物雜波。此外，Kessinger等[13]、Ferrer等[14]、Cho等[15]、劉黎平等[16]和江源等[17]在雷達基數(shù)據(jù)中提取各種特征，使用模糊邏輯算法識別并去除地物雜波。李豐等[18-19]采用類似的模糊邏輯算法識別非降水回波，并對比分析S波段和C波段多普勒天氣雷達各種特征的分布。

Ruiz等[20]使用相控陣天氣雷達探測數(shù)據(jù)，提取回波強度與時間的相關性紋理特征等參數(shù)，并將其用于樸素貝葉斯算法，有效識別和去除地物雜波。但該算法運算效率較低，不能滿足快速掃描相控陣天氣雷達的產(chǎn)品生成需求。文中同時指出，與傳統(tǒng)拋物線天線天氣雷達相比，日本大阪大學的相控陣天氣雷達具有較強的旁瓣和較寬的波束，會導致雷達數(shù)據(jù)中存在更多的地物雜波，對地物雜波的識別和去除算法要求更高。

為了進一步研究小尺度強對流天氣系統(tǒng)，在國內(nèi)外網(wǎng)絡化和相控陣技術發(fā)展的背景下[21-23]，2015年中國氣象局氣象探測中心設計并聯(lián)合有關廠家開始研制陣列天氣雷達，它結(jié)合網(wǎng)絡化雷達和分布式相控陣技術，具有高時空分辨率的特點。本文采用基于模糊邏輯的地物識別算法，參考日本大阪大學的相控陣天氣雷達地物雜波去除算法，提出了一個新的參數(shù)——回波強度時間變化量(time variability of reflectivity factor,TVR)，并在不同天氣情況下驗證了地物雜波去除算法的有效性。

1 陣列天氣雷達簡介

陣列天氣雷達是一種分布式相控陣天氣雷達，由控制處理中心和至少3個相控陣接收發(fā)射子陣(簡稱收發(fā)子陣)構成，通過增加收發(fā)子陣數(shù)量擴大探測區(qū)域[24]。每3個相鄰的收發(fā)子陣為一組進行協(xié)同掃描，保證3個相鄰收發(fā)子陣在空間同一點的數(shù)據(jù)時差小于2 s。

收發(fā)子陣是陣列天氣雷達的前端，由于缺少獨立的控制和數(shù)據(jù)處理部分，不是完整的雷達，故被稱為收發(fā)子陣而非相控陣雷達。收發(fā)子陣由相控陣天線陣列、收發(fā)模塊陣列、信號處理陣列、方位旋轉(zhuǎn)伺服、同步和通信模塊等部分組成。水平方向采用機械掃描方式覆蓋0°～360°，垂直方向發(fā)射4個寬度為22.5°的波束覆蓋0°～90°，接收到回波信號后，通過數(shù)字波束得到16個平均波束寬度為1.6°的波束。2018年3月在長沙機場布設的陣列天氣雷達，天線旋轉(zhuǎn)速度為30 °·s-1，完成1個體掃僅需要12 s，約為現(xiàn)在業(yè)務中使用的多普勒天氣雷達常用體掃模式1個體掃時間的1/30。每個收發(fā)子陣的最大探測距離為20.28 km，徑向分辨率為30 m。

2 基于模糊邏輯的地物識別算法

2.1 數(shù)據(jù)預處理

由于雷達基數(shù)據(jù)中含有部分孤立的點狀或線狀回波(圖1a)，本文采用Zhang等[10]提出的孤立點濾波算法對雷達基數(shù)據(jù)進行預處理。具體算法如下：以空間中點X為中心位置，在方位、徑向和俯仰上5×5×1個距離庫中有效回波點的百分比為PX，如果PX小于預設閾值(默認為75%)，則認為該點為孤立點并將其剔除。圖1為孤立點濾波前后的回波對比。由圖1b可知，通過孤立點濾波可以消除大部分孤立回波點。

圖1 孤立點濾波前(a)、濾波后(b)回波強度Fig.1 Intensity echo of isolated point before filtering(a) and after filtering(b)

2.2 識別算法介紹

陣列天氣雷達的時間分辨率空間分辨率均高于常規(guī)多普勒天氣雷達。本文在已有地物特征參數(shù)基礎上加入表征回波強度隨時間變化的參數(shù),回波強度時間變化量(TVR)，作為模糊邏輯算法的輸入?yún)?shù)來識別地物雜波[25-27]。其他特征參數(shù)包括回波強度紋理(TDBZ)、回波強度垂直梯度(GDBZ)、徑向庫間回波強度變化程度(SPIN)、徑向速度平均值(MDVE)和速度譜寬平均值(MDSW)。這些參數(shù)的定義見式(1)～(7)[16]：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

其中，NA，NR表示在方位和徑向定義的計算范圍；i，j為范圍內(nèi)的索引；Z，V和SW分別為回波強度、徑向速度和速度譜寬；Zlow和Zup為對應的本層和上層PPI的回波強度；Zthresh為徑向上相鄰距離庫強度變化的閾值；ZN和ZL分別表示當前時刻和上一時刻的強度。在實際特征值計算過程中，設定NA=5，NR=5，Zthresh=2 dBZ，ZN和ZL的時間差為1 min。

除少數(shù)特殊地物外，大部分地物位置固定不動且形狀大小不隨時間改變，所以地物雜波的徑向速度主要分布在0附近，譜寬相對較小，回波強度時間變化量(TVR)接近0。相比降水回波，地物雜波大部分出現(xiàn)在低仰角數(shù)據(jù)中，強度在空間上變化大，且隨著仰角的抬高，強度迅速減小。因此地物雜波的強度紋理(TDBZ)和徑向庫間強度變化程度(SPIN)較大，強度的垂直梯度(GDBZ)為負值且絕對值比較大。

本文收集了2019年3—7月不同天氣情況下長沙機場陣列天氣雷達基數(shù)據(jù)，根據(jù)地物雜波和天氣回波的特征預先確定回波類型。判斷依據(jù)：地物主要出現(xiàn)在較低的仰角，且隨著仰角的抬高強度迅速減小，地物的徑向速度接近零，位置不隨時間改變；對于降水回波，可根據(jù)回波形狀、垂直結(jié)構、演變和移動等方面確定[16]。使用部分回波數(shù)據(jù)，統(tǒng)計得到地物和降水各個參數(shù)的概率分布情況(圖2)，其余回波數(shù)據(jù)驗證，檢驗算法的識別效果。

由圖2可知，地物的徑向速度平均值(MDVE)集中分布在0附近，與降水的區(qū)別較大(圖2a)。地物和降水回波的速度譜寬平均值(MDSW)重合面積較大，區(qū)別較小，未達到預期效果，因此，在識別算法中剔除該參數(shù)(圖2b)。與降水回波相比，地物的回波強度紋理(TDBZ)分布更加離散，符合地物雜波空間分布差異較大的特征(圖2c)。地物雜波的強度垂直梯度(GDBZ)存在少數(shù)大于0的部分，與預期的效果存在差異，但與降水回波明顯不同，不影響作為識別地物的輸入?yún)?shù)(圖2d)。地物和降水回波的強度徑向庫間變化程度(SPIN)分布的波峰位置區(qū)別明顯，地物雜波的波峰位置更靠右，即在計算區(qū)域內(nèi)，地物雜波強度波動超過2 dB的數(shù)據(jù)點更多(圖2e)?；夭◤姸葧r間變化量(TVR)在地物和降水回波中有明顯差異,地物雜波主要分布在0～2 dB范圍內(nèi)，而降水回波的分布則相對離散(圖2f)。綜上所述，地物和降水回波的大部分參數(shù)存在較明顯差異，但也有部分重合區(qū)域，如在回波強度紋理(TDBZ)靠近0的區(qū)域也有部分地物雜波，徑向速度平均值(MDVE)較小的區(qū)域也有降水回波。因此，使用單一特征參數(shù)難以識別地物雜波，需綜合使用多種特征參數(shù)才能較好地區(qū)分地物和降水回波。

為了簡化計算，本文根據(jù)地物和降水特征分布的范圍、交點和波峰位置確定采用梯形隸屬函數(shù)(圖3)。根據(jù)隸屬函數(shù)，計算每個點各參數(shù)的判據(jù)M。然后，對每個點的所有判據(jù)(N個)進行累加平均(式(8))，可得到每個點地物的判斷值EM(0～1)。該值越大，代表該點是地物的可能性越大，反之，代表降水回波的可能性越大。最后將得到的判斷值EM與預先設定的閾值進行比較，大于閾值判斷為地物，反之判斷為降水。

(8)

圖2 地物和降水回波各特征參數(shù)的概率分布(a)MDVE，(b)MDSW，(c)TDBZ，(d)GDBZ，(e)SPIN，(f)TVRFig.2 Probability distributions of ground clutter and precipitation echo(a)MDVE，(b)MDSW，(c)TDBZ，(d)GDBZ，(e)SPIN，(f)TVR

圖3 地物識別的隸屬函數(shù)(a)MDVE，(b)TDBZ，(c)GDBZ，(d)SPIN，(e)TVRFig.3 Membership functions of ground clutter detection(a)MDVE，(b)TDBZ，(c)GDBZ，(d)SPIN，(e)TVR

3 地物識別效果

不同閾值的回波強度時間變化量(TVR)參數(shù)對識別效果影響見表1。隨著閾值的提高，地物識別準確率和降水誤判率降低，在閾值不大于0.55時，TVR對改善地物識別準確率和降水識別誤判率有貢獻。當閾值大于0.55后，TVR不再對識別算法產(chǎn)生影響。

雖然閾值為0.4時地物識別的準確率最高，但降水識別的誤判率也最大。綜合考慮地物識別的準確率和降水識別誤判率，最終選取閾值為0.45，并用于分析3個子陣在無降水、混合性降水和對流性降水天氣情況的地物識別算法效果。

表1 地物識別準確率和降水識別誤判率Table 1 Accuracy of ground clutter detection and error rate of precipitation detection

3.1 無降水情況地物識別效果

為了檢驗算法在無降水情況下的地物識別效果，采用2019年3—7月無降水的雷達數(shù)據(jù)，分析了3個子陣在無降水情況下算法識別效果(表2)。結(jié)果表明：地物識別算法在無降水情況下表現(xiàn)良好，地物識別的準確率高達96%。圖4為子陣1在2019年7月31日10：21地物識別前后的對比，經(jīng)過地物識別算法，大部分地物雜波可被有效識別并去除。

表2 無降水情況下地物識別算法效果Table 2 Accuracy of ground clutter detection algorithm under no precipitation condition

圖4 2019年7月31日10：21子陣1在1.4°仰角的回波強度(相鄰距離圈間距為5 km)(a)地物識別前，(b)地物識別后Fig.4 The echo intensity of subarray 1 at 1.4° elevation angle at 1027 BT 31 Jul 2019(the distance between adjacent rang rings is 5 km)(a)before ground clutter detection,(b)after ground clutter detection

3.2 混合性降水情況地物識別效果

分析混合性降水情況下地物識別算法效果(表3)。在此類降水情況下，各個子陣地物雜波識別的準確率均達到91%以上，而降水回波識別的誤判率最高僅10%。圖5為子陣1在2019年6月21日15：17地物識別前后的對比。由圖5可知，降水識別誤判率較高主要由降水回波邊緣部分的誤判引起。相對于無降水情況，地物識別準確率略有下降，從速度對比可以看到，與降水回波重合的地物雜波部分未被準確識別。

表3 混合性降水算法識別效果Table 3 Accuracy of ground clutter detection and error rate of precipitation detection under mixed prcipitation condition

圖5 2019年6月21日15：17子陣1在1.4°仰角的回波強度和徑向速度(相鄰距離圈間距為5 km)(a)識別前的回波強度，(b)識別后的回波強度，(c)識別前的徑向速度，(d)識別后的徑向速度Fig.5 The echo intensity and radial velocity of subarray 1 at 1.4°elevation angle at 1517 BT 21 Jun 2019(the distance between adjacent range rings is 5 km)(a)echo intensity before ground clutter detection,(b)echo intensity after ground clutter detection,(c)radial velocity before ground clutter detection,(d)radial velocity after ground clutter detection

3.3 對流性降水情況地物識別效果

分析對流性降水情況下地物識別算法效果見表4。由于子陣1和子陣3的個例中地物雜波和降水回波無重合部分，地物識別的準確率與無降水情況下差別不大，分別達到92%和94%。子陣2的大部分個例中降水和地物雜波重合區(qū)域中存在不能被準確識別部分，準確率僅91%。圖6和圖7分別是2019年7月21日13：57子陣2的1.4°和2.8°仰角地物識別前后的對比，2.8°仰角地物明顯減少，降水區(qū)域被更完整保留。因此，該類型降水情況下，算法識別處理后，大部分降水回波被較好保留，僅降水回波邊緣部分和西南方向一小塊孤立的降水回波被錯誤識別。

表4 對流性降水算法識別效果Table 4 Accuracy of ground clutter detection and error rate of precipitation detection under convective precipitation condition

圖6 2019年7月21日13：57子陣2在1.4°仰角的回波強度和徑向速度(相鄰距離圈間距為5 km)(a)識別前的回波強度，(b)識別后的回波強度，(c)識別前的徑向速度，(d)識別后的徑向速度Fig.6 The echo intensity and radial velocity of subarray 2 at 1.4°elevation angle at 1357 BT 21 Jul 2019(the distance between adjacent range rings is 5 km)(a)echo intensity before ground clutter detection,(b)echo intensity after ground clutter detection,(c)radial velocity before ground clutter detection,(d)radial velocity after ground clutter detection

圖7 2019年7月21日13：57子陣2在2.8°仰角的回波強度和徑向速度(相鄰距離圈間距為5 km)(a)識別前的回波強度，(b)識別后的回波強度，(c)識別前的徑向速度，(d)識別后的徑向速度Fig.7 The echo intensity and radial velocity of subarray 2 at 2.8° elevation angle at 1357 BT 21 Jul 2019(the distance between adjacent range rings is 5 km)(a)echo intensity before ground clutter detection,(b)echo intensity after ground clutter detection,(c)radial velocity before ground clutter detection,(d)radial velocity after ground clutter detection

續(xù)圖7

4 小 結(jié)

本文利用陣列天氣雷達高時空分辨率的特點，提出回波強度時間變化量(TVR)這一新參數(shù)。采用2019年3—7月長沙黃花機場X波段陣列天氣雷達數(shù)據(jù)，分析地物識別算法中6個特征參數(shù)的統(tǒng)計特征，改進基于模糊邏輯的地物識別算法，并利用雷達基數(shù)據(jù)分析TVR對算法的貢獻和不同天氣情況下識別算法的效果，得到以下結(jié)論：

1) 降水回波和地物雜波回波強度紋理(TDBZ)、回波強度垂直梯度(GDBZ)、徑向庫間回波強度變化程度(SPIN)、回波強度時間變化量(TVR)和徑向速度平均值(MDVE)5個參數(shù)特征分布差異較大，但速度譜寬平均值(MDSW)差異較小，因此，只采用前5個參數(shù)作為模糊邏輯算法的輸入特征。

2) 對回波強度時間變化量(TVR)的貢獻程度進行定量評價，通過選擇最優(yōu)的閾值，采用該參數(shù)可將地物識別準確率提高4%，降水識別誤判率降低2%。

3) 在不同天氣情況下分別對地物識別算法的表現(xiàn)進行驗證。無降水時，地物識別的準確率為94%；有降水時，地物識別準確率為92%，降水識別誤判率穩(wěn)定在10%左右。

本算法在陣列天氣雷達地物識別中取得了較好的效果，為生成更加準確的氣象產(chǎn)品提供了保障。但還存在一些問題：降水回波的邊緣部分可能被誤判，且當降水回波與地物回波混合在一起時，去除地物回波也損失了降水回波信息。在后續(xù)的研究中，將考慮采用多子陣聯(lián)合的地物去除方法，進一步提高識別準確率，以及補償與地物一起去除的降水回波。