聚類分析在X波段雙偏振雷達(dá)相態(tài)識(shí)別中的應(yīng)用

李宗飛 李祥海 陳凱華 馮亮

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

引言

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

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

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

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

1 資料來(lái)源

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

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

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

2.1 聚類分析

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

(1)

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

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

(2)

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

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

不同相態(tài)水凝物在雙偏振雷達(dá)參數(shù)上的表現(xiàn)是識(shí)別相態(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)識(shí)別[12]。模糊邏輯算法進(jìn)行相態(tài)識(shí)別是在此基礎(chǔ)上建立成員函數(shù)對(duì)雷達(dá)參數(shù)進(jìn)行模糊化，然后按照一定規(guī)則進(jìn)行數(shù)據(jù)合成，最后進(jìn)行退模糊化和識(shí)別相態(tài)。

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

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

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

3.1 資料預(yù)處理

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

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

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

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

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

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

在不考慮溫度分布的情況下，對(duì)雷達(dá)4種參數(shù)進(jìn)行聚類分析，使其充分分類，得到7個(gè)分類，聚類中心分別為小雨(含中雨){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)，主要因?yàn)樵搮^(qū)域反射率、差分反射率存在較大衰減，并且因?yàn)樾旁氡茸儾钏运胶痛怪狈较蚧夭ǖ南嚓P(guān)性變差，導(dǎo)致雷達(dá)參數(shù)不能用于相態(tài)識(shí)別。從圖2中看出在對(duì)流中心近地面位置主要是大雨、雹和雨夾雹，對(duì)流中心的高空為雹，外圍為霰。零度層亮帶識(shí)別為大雨和霰(此處無(wú)法正確分類)，零度層亮帶上方統(tǒng)一識(shí)別成雪(或有冰晶存在)，該區(qū)域的錯(cuò)誤識(shí)別皆因聚類分析算法自身局限性導(dǎo)致，不能對(duì)所有水凝物進(jìn)行完全分類，考慮引入溫度數(shù)據(jù)有可能解決該類問(wèn)題，但由于對(duì)流性降水溫度分布較為復(fù)雜，所以留待后續(xù)研究解決。

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

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

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

圖4為相態(tài)識(shí)別結(jié)果，聚類算法對(duì)各降水相態(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é)果可以看出對(duì)流中心位置5 km以上主要為冰雹，5 km以下主要為雨夾雹，而水平方向64 km以外存在因?yàn)樾盘?hào)衰減而無(wú)法觀測(cè)的冰雹區(qū)域。在5 km以上對(duì)流中心兩側(cè)主要為霰粒子，最外圍是雪或冰晶，其下方為降雨。圖中可以看出在降水回波邊緣位置，依然存在因回波信噪比較差而無(wú)法識(shí)別的未知降水類型。

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

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

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

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

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

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

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

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

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

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

4 結(jié)論與討論

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

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

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