李宗飛 李祥海 陳凱華 馮亮

(1 天津市氣象信息中心，天津 300000; 2 中國科學(xué)院大氣物理研究所云降水物理與強風(fēng)暴重點實驗室，北京 100029)

引言

隨著我國電子硬件和計算機軟件技術(shù)的發(fā)展，新一代多普勒天氣雷達(dá)正在向雙偏振天氣雷達(dá)升級。雙偏振天氣雷達(dá)不僅包含反射率因子(ZH)、徑向速度(V)、譜寬(W)，還包括了差分反射率(ZDR)、相關(guān)系數(shù)(ρhv)、差傳播相移(ΦDP)、差傳播相移率(KDP)等差分信息。雙偏振天氣雷達(dá)在冰雹、暴雨等災(zāi)害性天氣監(jiān)測預(yù)警中起到了重要作用[1]，能更好地計算降水量、雨滴譜等信息[2]，并且可以進(jìn)行水凝物相態(tài)識別，對于分析降水的物理過程有著重要意義。

雙偏振天氣雷達(dá)在相態(tài)識別方面表現(xiàn)較好，1993年Straka和Zrnic提出的布爾判決邏輯進(jìn)行水凝物分類[3]，由于觀測噪聲以及雷達(dá)參數(shù)之間的不完全排斥，所以該方法并不能完全分類降水相態(tài)。通過建立概率模型的方法也可以進(jìn)行相態(tài)識別，需要統(tǒng)計相應(yīng)的概率密度函數(shù)，并且需要先驗知識，實現(xiàn)較為困難[4]。模糊邏輯算法是目前較為常用的相態(tài)識別方法，模糊邏輯算法最早由Zadeh(1965)提出，曹俊武(2005)、何宇翔(2010)、劉亞男(2012)等先后使用該方法實現(xiàn)降水粒子的相態(tài)識別，模糊邏輯算法進(jìn)行相態(tài)識別首先需要建立成員函數(shù)，成員函數(shù)是指不同水凝物粒子的雷達(dá)參數(shù)取值范圍的函數(shù)，又稱為隸屬函數(shù)，成員函數(shù)的優(yōu)劣直接決定了相態(tài)識別效果[5-7]。

成員函數(shù)的建立有通過散射模擬獲得，也有根據(jù)觀測以及前人經(jīng)驗進(jìn)行修改得到，然而由于不同雷達(dá)接收機靈敏度、動態(tài)范圍、噪聲系數(shù)等差異[8]，該方法得到的成員函數(shù)對不同雷達(dá)的識別效果存在差異，尤其是當(dāng)雷達(dá)參數(shù)存在系統(tǒng)性偏差時，如雙偏振雷達(dá)水平和垂直通道增益和靈敏度不同，導(dǎo)致差分反射率存在系統(tǒng)性偏差，此時常用的成員函數(shù)與實際雷達(dá)觀測數(shù)據(jù)存在差異，不利于進(jìn)行相態(tài)識別[9]。為了得到更合適的成員函數(shù)，可使用雷達(dá)和機載粒子譜儀聯(lián)合觀測，然后調(diào)整相應(yīng)的的成員函數(shù)，該方法實驗難度較大，且觀測成本較高。針對此類情況，本文提出使用聚類分析的方法對水凝物粒子進(jìn)行分類，然后根據(jù)聚類中心判斷降水相態(tài)。

因為雷達(dá)老化出現(xiàn)系統(tǒng)性偏差是一個漸變過程，不能及時對雷達(dá)進(jìn)行定標(biāo)或系統(tǒng)性補償，該過程中采用模糊邏輯算法進(jìn)行相態(tài)識別將存在一定的誤差，所以本文使用聚類分析的方法完成相態(tài)識別，然后統(tǒng)計不同相態(tài)粒子的雷達(dá)參數(shù)概率分布特征，計算相應(yīng)的概率分布函數(shù)，對該函數(shù)進(jìn)行歸一化處理并構(gòu)建成員函數(shù)，可用于模糊邏輯算法進(jìn)行相態(tài)識別。因為該函數(shù)通過統(tǒng)計得到，所以不受系統(tǒng)偏差影響，更適合于該雷達(dá)使用。

1 資料來源

使用數(shù)據(jù)來自于中國科學(xué)院大氣物理研究所的車載X波段雙偏振雷達(dá)，型號為714XDP-A，觀測地點為北京市順義區(qū)，經(jīng)緯度分別為116°40′48″E，40°10′48″N，雷達(dá)工作頻率9.370 GHz，距離分辨率150 m/250 m，最大探測距離150 km。該雷達(dá)工作年限較長，存在器件老化現(xiàn)象，其反射率ZH和差分反射率ZDR存在系統(tǒng)偏差，相關(guān)系數(shù)ρhv在信噪比較低時其相關(guān)性變差，為估計降水和相態(tài)識別增加了難度。

使用了2015年7月和2016年9月發(fā)生在北京的降雹過程資料，雷達(dá)掃描方式為垂直掃描(RHI)，用于相態(tài)識別的數(shù)據(jù)包括ZH、ZDR、KDP和ρhv。其中反射率ZH和差分反射率ZDR存在較大的降水衰減，所以本文使用了KDP綜合分類法進(jìn)行衰減訂正[10]。另外，為了減少數(shù)據(jù)中的噪聲，提高數(shù)據(jù)可用性，本文對以上四種數(shù)據(jù)進(jìn)行平滑處理。

2 聚類分析識別相態(tài)原理

2.1 聚類分析

聚類分析是一種群分析手段，是對樣品或指標(biāo)進(jìn)行多元統(tǒng)計的方法[11]。K-均值聚類方法通過計算樣本之間的距離進(jìn)行分類，可以將樣本或指標(biāo)看成n維空間上的點，例如樣本x1={x11x12x13},和x2={x21x22x23},在三維空上的歐氏距離d為

(1)

K-均值聚類方法由如下4個步驟完成：

j≠i},i=1,2,…,k

(2)

最終的各個聚點是每個分類的重心，具有各分類的鮮明特征，而且每個聚點之間有較大的歐氏距離，是區(qū)分各個分類的標(biāo)志性數(shù)據(jù)。另外，由于初始聚點為隨機數(shù)據(jù)，受算法自身的局限性影響，聚類分析可能會陷入局部穩(wěn)定的狀態(tài)，不能進(jìn)行充分的分類，需要重復(fù)計算，才能得到最佳結(jié)果。此為K均值聚類分析方法的自身局限性，本人通過軟件重復(fù)分類，選擇充分分類的結(jié)果做為最終分類，利用計算機可重復(fù)計算和不斷迭代的特點彌補算法的不足，表明該算法可以進(jìn)一步的優(yōu)化。

2.2 相態(tài)識別原理

不同相態(tài)水凝物在雙偏振雷達(dá)參數(shù)上的表現(xiàn)是識別相態(tài)的基礎(chǔ)，例如冰雹的反射率通常在35～50 dBz之間，差傳播相移率在-3～1 °/km，差分反射率在-3～0 dB之間，相關(guān)系數(shù)在0.94～0.98范圍之間，所以，能夠利用不同水凝物在雷達(dá)參數(shù)上的表現(xiàn)不同進(jìn)行相態(tài)識別[12]。模糊邏輯算法進(jìn)行相態(tài)識別是在此基礎(chǔ)上建立成員函數(shù)對雷達(dá)參數(shù)進(jìn)行模糊化，然后按照一定規(guī)則進(jìn)行數(shù)據(jù)合成，最后進(jìn)行退模糊化和識別相態(tài)。

根據(jù)不同水凝物雷達(dá)參數(shù)特征的表現(xiàn)，表明各水凝物的雷達(dá)參數(shù)存在歐氏距離，所以能夠使用聚類分析算法對不同相態(tài)的水凝物進(jìn)行區(qū)分。能夠反映水凝物相態(tài)的雷達(dá)參數(shù)包括反射率因子ZH(dBz)，差分反射率因子ZDR(dB)，差傳播相移率KDP(°/km)以及相關(guān)系數(shù)ρhv(無量綱)。首先，將相同位置的4個雷達(dá)參數(shù)構(gòu)成一個樣本(ZH，ZDR，KDP，ρhv)，然后對一定范圍的所有樣本進(jìn)行聚類分析，最后對分類結(jié)果的聚類中心進(jìn)行相態(tài)判斷，因為聚點與聚點之間有較大的歐氏距離，邊界清晰，所以很容易根據(jù)不同相態(tài)的降水特征進(jìn)行判別。

聚類分析完成相態(tài)識別后，可一進(jìn)步統(tǒng)計水凝物的雷達(dá)參數(shù)概率分布特征，并對其進(jìn)行函數(shù)擬合，從而建立不同水凝物的成員函數(shù)，作為該雷達(dá)專有的成員函數(shù)，可用于模糊邏輯算法進(jìn)行相態(tài)識別。

3 聚類分析相態(tài)識別實例

3.1 資料預(yù)處理

由于在強降水過程中X波段雷達(dá)存在較大的降水衰減，所以首先需要對雷達(dá)反射率及差分反射率進(jìn)行衰減訂正，本文使用KDP綜合分類法進(jìn)行衰減訂正，該方法分別對ZH和KDP進(jìn)行分類，然后采用自適應(yīng)算法計算不同分類中的衰減系數(shù)，并制定了衰減系數(shù)表格，最后使用表格數(shù)據(jù)計算衰減率，并進(jìn)行反射率和差分反射率的衰減訂正[10]。

聚類分析是按照各樣本之間的歐氏距離進(jìn)行分類的，為了解決不同參數(shù)在歐氏距離計算時的公平問題，需要對所有參數(shù)進(jìn)行歸一化，例如反射率因子取值范圍在-5～70 dBz之間，而相關(guān)系數(shù)取值范圍在0～1之間，如果不進(jìn)行歸一化，那相關(guān)系數(shù)對于距離的影響要遠(yuǎn)遠(yuǎn)小于反射率因子，不能充分發(fā)揮相關(guān)系數(shù)的作用。固態(tài)降水識別主要依賴相位數(shù)據(jù)，為提高相位參數(shù)的權(quán)重，可以對KDP和ρhv乘以大于1的系數(shù)，增加該參數(shù)的權(quán)重。另外，受接收機靈敏度和噪聲系數(shù)限制，當(dāng)回波較弱時雙偏振雷達(dá)的偏振參數(shù)誤差較大，如KDP和ρhv，所以需要剔除反射率小于5 dBz的回波，然后進(jìn)行相態(tài)識別。

3.2 相態(tài)識別實現(xiàn)

圖1為2016年9月7日18:48(北京時)冰雹降水的垂直掃描過程，所處方位角為212°，圖1a反射率和圖1b差分反射率均進(jìn)行了降水衰減訂正，回波強度在對流中心處得到一定恢復(fù)，但對流中心后方因電磁波無法穿透出現(xiàn)空洞，并且在回波邊緣處由于信噪比較低ZDR、KDP和ρhv均出現(xiàn)異常值，如ρhv在50 km后因小于0.7超出色標(biāo)范圍而無法顯示。

圖1中可以看出云團(tuán)頂高達(dá)到12 km，并且垂直3 km處可明顯看到零度層亮帶，距離雷達(dá)32 km處為對流中心，最強回波達(dá)到60 dBz左右。圖中顯示零度層亮帶以下反射率小于30 dBz， 差分反射率小于1.5 dB，差傳播相移率在0°/km左右，相關(guān)系統(tǒng)在0.98左右，表明該區(qū)域為小到中雨量級。在對流中心處反射率35～60 dBz之間，差分反射率從-3～3 dB之間，差傳播相移率從0～8 °/km，相關(guān)系統(tǒng)在0.88～0.96之間，且對流中心延伸到10 km高度，表明此處為強降雨和冰雹區(qū)域。零度層亮帶一般為融雪，其上多為雪和冰晶[13]。在近地面對流中心的后方區(qū)域，受降水衰減影響反射率為0～5 dBz。圖中反射率和差分反射率均已進(jìn)行衰減訂正，但差分反射仍存在小于-3 dB的數(shù)據(jù)，超出色標(biāo)顯示范圍，該區(qū)域差傳播相移率和相關(guān)系數(shù)信噪比較低，且有些數(shù)據(jù)低于色標(biāo)顯示范圍，不進(jìn)行相態(tài)識別。

圖1 2016年9月7日18:48北京X波段雙偏振雷達(dá)垂直剖面(方位212°):(a)訂正后反射率，(b)訂正后差分反射率，(c)差傳播相移率，(d)相關(guān)系數(shù)

在不考慮溫度分布的情況下，對雷達(dá)4種參數(shù)進(jìn)行聚類分析，使其充分分類，得到7個分類，聚類中心分別為小雨(含中雨){17.42 dBz，-0.18 dB，-0.99°/km，0.970}、大雨{36.49 dBz，0.64 dB，2.04°/km，0.951}、雹{36.81 dBz，-1.47 dB，-0.83°/km，0.963}、雨夾雹{47.97 dBz，1.53 dB，7.50°/km，0.923}、雪或冰晶{11.90 dBz，-1.02 dB，-1.26°/km，0.977}、霰{21.48 dBz，-1.28 dB，-1.130°/km，0.979}、未知{7.84 dBz，-1.18 dB，-1.84°/km，0.778}。圖2中深灰色為未知相態(tài)，主要因為該區(qū)域反射率、差分反射率存在較大衰減，并且因為信噪比變差所以水平和垂直方向回波的相關(guān)性變差，導(dǎo)致雷達(dá)參數(shù)不能用于相態(tài)識別。從圖2中看出在對流中心近地面位置主要是大雨、雹和雨夾雹，對流中心的高空為雹，外圍為霰。零度層亮帶識別為大雨和霰(此處無法正確分類)，零度層亮帶上方統(tǒng)一識別成雪(或有冰晶存在)，該區(qū)域的錯誤識別皆因聚類分析算法自身局限性導(dǎo)致，不能對所有水凝物進(jìn)行完全分類，考慮引入溫度數(shù)據(jù)有可能解決該類問題，但由于對流性降水溫度分布較為復(fù)雜，所以留待后續(xù)研究解決。

圖2 2016年9月7日18:48北京X雙偏振雷達(dá)相態(tài)識別

圖3為2015年7月29日22：57(北京時)觀測數(shù)據(jù)，該數(shù)據(jù)為冰雹過程的垂直掃描數(shù)據(jù)，方位角為190°。圖3中分別為訂正后的反射率，訂正后差分反射率，差傳播相移率和相關(guān)系數(shù)，圖中顯示本過程為強對流過程，對流高度達(dá)到12 km以上，最大強度為50 dBz。數(shù)據(jù)在水平方向的64 km之后因電磁波無法穿透而完全衰減。圖中在水平方向60～64 km處對流高度及強度均為最大區(qū)域，是該過程的對流中心，該區(qū)域的差分反射率和相關(guān)系數(shù)均表現(xiàn)出冰雹特征。圖中反射率和差分反射率均已進(jìn)行降雨衰減訂正，差分反射率低于-3 dB的區(qū)域未能顯示，不進(jìn)行相態(tài)識別。

圖3 2015年7月29日22:57北京X波段雙偏振雷達(dá)垂直剖面(方位190°):(a)訂正后反射率，(b)訂正后差分反射率，(c)差傳播相移率，(d)相關(guān)系數(shù)

圖4為相態(tài)識別結(jié)果，聚類算法對各降水相態(tài)分類中心為小雨(含中雨){14.15 dBz，-0.13 dB，-1.14°/km，0.968}、大雨{35.17 dBz，0.9 dB，1.04°/km，0.976}、雹{38.02dBz，-1.22dB，-0.91°/km，0.956}、雨夾雹{41.76 dBz，0.74 dB，3.47°/km，0.942}、雪或冰晶{13.16 dBz，-1.23 dB，-1.20°/km，0.943}、霰{25.99 dBz，-1.01 dB，-0.93°/km，0.969}、未知{11.59 dBz，-0.42 dB，1.1°/km，0.794}。從分類結(jié)果可以看出對流中心位置5 km以上主要為冰雹，5 km以下主要為雨夾雹，而水平方向64 km以外存在因為信號衰減而無法觀測的冰雹區(qū)域。在5 km以上對流中心兩側(cè)主要為霰粒子，最外圍是雪或冰晶，其下方為降雨。圖中可以看出在降水回波邊緣位置，依然存在因回波信噪比較差而無法識別的未知降水類型。

圖4 2015年7月29日22:57北京X雙偏振雷達(dá)相態(tài)識別

聚類中心即為該相態(tài)粒子雷達(dá)參數(shù)的平均值，表明該相態(tài)粒子的平均特征，從數(shù)據(jù)上看與以往經(jīng)驗數(shù)據(jù)略有差異，各類水凝物雷達(dá)參數(shù)均值偏低，表明雷達(dá)可能存在老化現(xiàn)象，回波接收功率及回波相關(guān)性變差。

3.3 統(tǒng)計參數(shù)特征

本文按照以上方法對冰雹過程的多次RHI掃描進(jìn)行相態(tài)識別，并將雷達(dá)各類參數(shù)根據(jù)水凝物進(jìn)行分類，然后統(tǒng)計各參數(shù)的概率密度分布，并進(jìn)行歸一化， 最后得到屬于該雷達(dá)在本次過程的相態(tài)識別成員函數(shù)，可進(jìn)一步用于模糊邏輯算法進(jìn)行相態(tài)識別，圖5、圖6給出了2015年和2016年的大雨和冰雹的統(tǒng)計數(shù)據(jù)。

圖5為2016年9月7日的大雨和冰雹雷達(dá)參數(shù)的概率分布情況，從圖中可以看出，本次降水的冰雹和大雨反射率較為接近，主要分布在25～45 dBz之間，差分反射率和差傳播相移率有較大差異，其中冰雹的差分反射率和差傳播相移率主要為負(fù)值，而大雨為正值，兩者相關(guān)系數(shù)主要分布在0.97以上，差別不大。該統(tǒng)計顯示了該過程中冰雹及大雨的參數(shù)分布差異。

圖6為2015年7月29日的大雨和冰雹雷達(dá)參數(shù)的概率分布特征，與圖5相比本次過程的反射率、差分反射率和差傳播相移率的中心位置均小于圖6，表明該過程弱于2016年的過程。從兩次過程反映出冰雹的差分反射率和差傳播相移率小于大雨過程，并且其分布范圍較窄。

圖5 2016年9月7日過程大雨和冰雹雷達(dá)參數(shù)概率分布(a)反射率,(b)差分反射率,(c)差傳播相移率,(d)相關(guān)系數(shù)

圖6 2015年7月29日過程大雨和冰雹雷達(dá)參數(shù)據(jù)概率分布(a)反射率,(b)差分反射率,(c)差傳播相移率,(d)相關(guān)系數(shù)

經(jīng)與劉亞男2012年使用模糊邏輯算法相比[7]，雷達(dá)參數(shù)的分布特征存在一定的差異，首先是反射率和差分反射率明顯偏弱，其中馮亮在2018年論文中提到該雷達(dá)差分反射率需要一定的補償[9]，表明該雷達(dá)確實存在老化現(xiàn)象。另外，差傳播相移率分布寬度窄于模糊邏輯算法中給出的成員函數(shù)，表明在單次過程中的參數(shù)分布包含于平均特征分布中，考慮增加更多降水個例，可能會趨于一致。

從圖5、圖6數(shù)據(jù)可以看出，在本次過程中，冰雹與大雨的差分反射率和差傳播相移率的分布特征有較大的差異，以此分布曲線擬合成員函數(shù)，在模糊邏輯算法中將明確分辨出兩種降水類型，所以使用聚類分析方法構(gòu)造成員函數(shù)用于模糊邏輯算法是可行的。

4 結(jié)論與討論

K-均值聚類分析方法在雙偏振雷達(dá)相態(tài)識別中應(yīng)用效果較好，該方法區(qū)別于傳統(tǒng)算法，是對大數(shù)據(jù)算法的一種典型應(yīng)用，彌補了傳統(tǒng)算法在相態(tài)識別中對先驗數(shù)據(jù)的依賴，如模糊邏輯算法完全依賴于成員函數(shù)。但該算法同樣存在自身的局限性，如在計算時存在局部穩(wěn)定狀態(tài)，需要重復(fù)聚類才能解決。該問題可采用其他聚類分析方法，擇優(yōu)處理。另外，該方法是在當(dāng)前大數(shù)據(jù)應(yīng)用環(huán)境下的一種嘗試，是模糊邏輯算法的一種補充方案。

經(jīng)聚類分析進(jìn)行相態(tài)識別后的數(shù)據(jù)，可分別進(jìn)行概率分布統(tǒng)計，完成成員函數(shù)的構(gòu)建，通過該方法構(gòu)建的成員函數(shù)更為適合該雷達(dá)進(jìn)行水凝物分類，并且不受雷達(dá)系統(tǒng)誤差的影響。因為不同降水類型的數(shù)據(jù)分布特征不一致，所以需要對多類型降水進(jìn)行完全統(tǒng)計，才能構(gòu)建合理的各成員函數(shù)。

溫度分布在相態(tài)識別中起著重要作用，在不考慮溫度的情況下，在雷達(dá)參數(shù)上無法識別融雪與霰和雨、雪和冰晶等的區(qū)別，因為它們的雷達(dá)參數(shù)不存在排斥，所以加入溫度信息是利用雷達(dá)參數(shù)進(jìn)行相態(tài)識別的重要手段。