基于P-K濾波的X波段雷達(dá)雨衰補(bǔ)償研究

李 海, 羅 原, 馮興寰, 馮 青

(中國(guó)民航大學(xué)天津市智能信號(hào)與圖像處理重點(diǎn)實(shí)驗(yàn)室，天津 300300)

0 引 言

相比單偏振雷達(dá)，雙偏振雷達(dá)通過發(fā)射水平和垂直偏振電磁波不僅能探測(cè)到常規(guī)的多普勒參量，還能獲取表征粒子相態(tài)和微物理特性的其他偏振參量。天氣雷達(dá)的偏振參量對(duì)于確定雷達(dá)校準(zhǔn)、地面雜波的消除、降水分類和雨滴譜參數(shù)估計(jì)等方面具有更大的優(yōu)勢(shì)[1]。相比C波段和S波段天氣雷達(dá)，X波段天氣雷達(dá)對(duì)弱氣象目標(biāo)的探測(cè)更敏銳，并且具有天線尺寸小、易于移動(dòng)等特點(diǎn)。但由于降水區(qū)域會(huì)對(duì)天氣雷達(dá)的回波信號(hào)產(chǎn)生散射干擾以及吸收，會(huì)導(dǎo)致接收到的數(shù)據(jù)與實(shí)際值相比會(huì)有衰減，特別是波長(zhǎng)較短的X波段(中心波長(zhǎng)3 cm)，其回波信號(hào)的衰減現(xiàn)象更為嚴(yán)重，這也是X波段天氣雷達(dá)應(yīng)用推廣所面臨的主要問題。

X波段天氣雷達(dá)雨衰補(bǔ)償相比C波段和S波段天氣雷達(dá)的研究和發(fā)展較晚，因而對(duì)X波段天氣雷達(dá)雨衰補(bǔ)償研究的最開始階段是借鑒了C波段和S波段氣象雷達(dá)的相關(guān)方法以及研究結(jié)果[2]。2008年胡志群等人提出了差分傳播相移率-衰減率(Kdp-Adp)方法，并通過對(duì)Kdp設(shè)置閾值，對(duì)水平反射率因子(ZH)進(jìn)行衰減訂正，該方法有一定的效果，但由于其沒有考慮到Kdp的累積量在電磁波傳播路徑上會(huì)受到前向差分散射相移的影響，因而具有一定的局限性。2009年何宇翔等人引入卡爾曼濾波首先對(duì)差分傳播相移(Φdp)進(jìn)行濾波處理，進(jìn)而對(duì)ZH進(jìn)行衰減訂正，該方法有效地減小了水平反射率因子的誤差，但其處理過程是將降水區(qū)域視為線性關(guān)系，其具有一定的局限性。2014年魏慶等人采用線性滑窗、小波分析等方法對(duì)Φdp進(jìn)行了衰減訂正；2015年，孫躍、肖輝等人提出了線性擬合-遞推法(LFRM)對(duì)Φdp進(jìn)行質(zhì)量控制。這些方法都能對(duì)Φdp進(jìn)行一定程度的衰減訂正，但這些方法的共同點(diǎn)都是將降水區(qū)域視為一個(gè)連續(xù)的線性區(qū)域，但由于天氣雷達(dá)的回波功率只代表某一距離門的單位波束體積內(nèi)N個(gè)降水粒子的回波功率平均值，不能理想的將天氣雷達(dá)不同距離門的回波信號(hào)數(shù)據(jù)視為線性關(guān)系[3]，因而這些方法具有一定的局限性。

針對(duì)X波段天氣雷達(dá)的偏振參量衰減訂正，本文首先采用粒子濾波算法對(duì)Φdp進(jìn)行濾波處理，在經(jīng)過粒子濾波算法處理后，對(duì)數(shù)據(jù)利用卡爾曼濾波進(jìn)行最優(yōu)化估計(jì)，以達(dá)到對(duì)Φdp數(shù)據(jù)的精確估計(jì)訂正。在得到Φdp濾波數(shù)據(jù)后，將Φdp數(shù)據(jù)作為參照量采用自適應(yīng)算法，對(duì)ZH進(jìn)行衰減訂正處理。天氣雷達(dá)的偏振參量回波數(shù)據(jù)是一個(gè)距離門內(nèi)所有采樣點(diǎn)的疊加值，常規(guī)的窗函數(shù)濾波等方法對(duì)所包含的噪聲信號(hào)的濾除效果有限。粒子濾波能夠利用偏振參量之間的關(guān)系建立合適的狀態(tài)方程和觀測(cè)方程，可以有效地進(jìn)行數(shù)據(jù)的濾波處理，但由于粒子濾波的時(shí)間復(fù)雜度隨著粒子數(shù)的增加而增加，于是引入卡爾曼濾波，能夠?qū)?jīng)過較少粒子數(shù)的粒子濾波處理后數(shù)據(jù)進(jìn)行最優(yōu)化估計(jì)，從而達(dá)到時(shí)間復(fù)雜度的平衡[4]。采用P-K算法能夠得到更平滑和準(zhǔn)確的Φdp，進(jìn)而對(duì)水平反射率因子進(jìn)行更精確的衰減訂正。

1 基于P-K濾波的差分傳播相移衰減訂正

X波段雙偏振雷達(dá)回波數(shù)據(jù)包含有ZH,Kdp，Φdp等偏振參量信息。在對(duì)ZH進(jìn)行衰減訂正之前需要先對(duì)Φdp進(jìn)行衰減訂正。

1.1 建立粒子濾波方程

直接從天氣雷達(dá)回波數(shù)據(jù)中提取到的并不是差分傳播相移而是包含了前向差分散射相移和后向差分散射相移的全差分相移(Ψdp)，其中后向差分散射相移也就是差分傳播相移(Φdp)。單個(gè)距離門的全差分相移(Ψdp(k))定義為

Ψdp(k)=Φdp(k)+δhv(k)，k=1,…，K

(1)

式中：k表示傳播路徑電磁波到達(dá)的距離門，總共有K個(gè)距離單元；Φdp(k)表示第k個(gè)距離門的差分傳播相移；δhv(k)表示第k個(gè)距離門的前向差分散射相移，它是單基地天氣雷達(dá)回波信號(hào)中的有害高頻噪聲成分。X波段電磁波在小雨到中雨區(qū)域偏離瑞利散射很小，δhv(k)可以忽略，但對(duì)于較強(qiáng)的雨區(qū)，就可能產(chǎn)生較強(qiáng)的瑞利散射，因而不能被忽略，所以雨衰補(bǔ)償?shù)闹攸c(diǎn)就是從全差分相移中將δhv(k)濾除。

從全差分相移中無法直接通過濾波等方式將前向差分散射相移濾除，Hubbert擬合得到的不同頻率的雷達(dá)的δhv(k)與單個(gè)距離門的差分傳播相移率Kdp(k)有如下關(guān)系[5]：

δhv(k)=c+bKdp(k)，k=1,…，K

(2)

式中，b和c的取值依賴于Kdp(k)(k=1,…，K)的取值[6]，當(dāng)Kdp(k)≤2.5°/km時(shí)，b取2.37，c取0.054；當(dāng)Kdp(k)>2.5°/km時(shí)，b取0.27，c取6.16。

將式(2)代入式(1)可得到Ψdp(k)與Φdp(k)及Kdp(k)之間的關(guān)系：

Ψdp(k)-c=Φdp(k)+bKdp(k)，k=1,…,K

(3)

而式(1)中的差分傳播相移與差分傳播相移率之間又有如下關(guān)系[7]：

Φdp(k+1)=Φdp(k)+2Δr·Kdp(k)，

k=1,…，K

(4)

式中，Δr表示相鄰距離門之間的距離。

根據(jù)式(3)和式(4)即可建立如下的粒子濾波方程組[8]：

k=1,…,K

(5)

k=1,…,K

(6)

(7)

以觀測(cè)量與預(yù)測(cè)狀態(tài)之間的差值作為似然函數(shù)，計(jì)算公式如式(8)表示：

(8)

可以通過觀測(cè)方程迭代更新重要性權(quán)值，更新步為

(9)

權(quán)值進(jìn)行歸一化可得

(10)

進(jìn)而可以得到狀態(tài)xk的估計(jì)為

(11)

即

(12)

1.2 建立卡爾曼濾波方程

粒子濾波通過尋找一組在狀態(tài)空間的隨機(jī)樣本對(duì)概率密度函數(shù)進(jìn)行近似，以后驗(yàn)概率密度所產(chǎn)生的N個(gè)獨(dú)立同分布樣本的均值代替積分運(yùn)算，從而獲得狀態(tài)最小方差分布，但由于粒子數(shù)的增加會(huì)導(dǎo)致其運(yùn)算量急劇增加，而較少的粒子又不能得到最優(yōu)的最小方差狀態(tài)值，因而建立卡爾曼濾波對(duì)數(shù)據(jù)進(jìn)行再處理，從而得到結(jié)果值。

(13)

(14)

卡爾曼濾波主要包括預(yù)測(cè)和更新兩個(gè)過程，預(yù)測(cè)過程如式(15)、式(16)所示[10]：

(15)

(16)

更新過程是結(jié)合先驗(yàn)狀態(tài)估計(jì)與當(dāng)前距離門觀測(cè)向量的數(shù)據(jù)進(jìn)行后驗(yàn)估計(jì)過程，更新過程為

(17)

進(jìn)而得到卡爾曼濾波估計(jì)方程為

(18)

而協(xié)方差更新方程為

(19)

(20)

2 基于自適應(yīng)算法的水平反射率因子衰減訂正

對(duì)差分傳播相移的衰減訂正是對(duì)水平反射率因子訂正的基礎(chǔ)，在利用P-K濾波算法得到誤差更小的差分傳播相移數(shù)據(jù)后，利用自適應(yīng)算法，對(duì)反射率因子(ZH)進(jìn)行衰減訂正。

利用Bringi提出的自適應(yīng)約束(self-consistent method with constraints)算法對(duì)ZH進(jìn)行衰減訂正[11]。ZH的衰減訂正的關(guān)鍵就是確定雨區(qū)的衰減率AH，設(shè)定雨區(qū)的距離門范圍為k0～k1，可得到第k個(gè)距離門的衰減率AH(k)的計(jì)算公式為[12]

AH(k)=

(21)

(22)

(23)

3 實(shí)驗(yàn)結(jié)果及分析

3.1 P-K濾波結(jié)果分析

本仿真實(shí)驗(yàn)的X波段氣象回波數(shù)據(jù)來源于大氣輻射測(cè)量氣候研究中心(ARM)的網(wǎng)站。其數(shù)據(jù)包括了Φdp，Kdp，ZH等重要回波參數(shù)，每個(gè)徑向上有400個(gè)距離單元，每個(gè)距離單元代表100 m。實(shí)驗(yàn)選用的X波段雷達(dá)于2018年7月3日11:00的回波數(shù)據(jù)。

一個(gè)徑向上的Φdp數(shù)據(jù)經(jīng)過不同方法的濾波效果與源數(shù)據(jù)效果對(duì)比如圖1所示。

圖1 同一徑向不同方法的Φdp濾波效果圖

由圖1可以看出，相對(duì)比原始數(shù)據(jù)，卡爾曼濾波以及P-K濾波均能對(duì)反射率因子的波動(dòng)和高頻噪聲進(jìn)行抑制，保證了廓線的連續(xù)性和平滑度。但相比卡爾曼濾波，P-K濾波能夠?qū)?shù)據(jù)進(jìn)行更好的平滑處理。

如圖2所示，為一個(gè)俯仰角的全徑向PPI圖對(duì)比，與卡爾曼濾波結(jié)果相比較，P-K濾波能更有效地將數(shù)據(jù)中的高頻噪聲部分剔除，且數(shù)據(jù)能呈現(xiàn)較好的非負(fù)性，更加符合真實(shí)氣象特征[14]。

(a) 卡爾曼濾波

3.2 水平反射率因子訂正結(jié)果及分析

圖3為衰減訂正前后的ZH的一個(gè)徑向距離廓線，由圖可以看到相比與卡爾曼濾波，P-K濾波在衰減訂正的同時(shí)，保留了更多的原始數(shù)據(jù)的輪廓細(xì)節(jié)。

圖3 一個(gè)徑向ZH訂正結(jié)果顯示

圖4(a)是相同時(shí)間，相近地點(diǎn)的KVNX雷達(dá)的S波段的ZH圖；圖4(b)是衰減訂正前的ZH數(shù)據(jù)的PPI圖；圖4(c)、(d)是應(yīng)用卡爾曼濾波、P-K濾波處理后進(jìn)行衰減訂正后的ZH的PPI圖。由于S波段的天氣雷達(dá)其雨衰相對(duì)較小，因而可以當(dāng)作參考，對(duì)比黑色方框區(qū)域可看到，P-K濾波的結(jié)果相對(duì)卡爾曼濾波結(jié)果，其雨衰補(bǔ)償訂正效果更好，與S波段KVNX 雷達(dá)的參考值更加接近。并且在雷達(dá)遠(yuǎn)端低SNR的條件下，依然具有良好的訂正效果。

(a) S波段ZHPPI掃描圖

通過散射模擬建立的偏振參量的經(jīng)驗(yàn)關(guān)系驗(yàn)證衰減訂正的效果，可以比較X波段訂正前后的AH～ZH和ZH～Kdp之間的散點(diǎn)圖特性。圖5(a)、(b)分別為訂正前后的ZH～Kdp的散點(diǎn)圖，實(shí)線為Park通過散射模擬建立的ZH～Kdp的經(jīng)驗(yàn)關(guān)系，由圖5(a)可以發(fā)現(xiàn)，經(jīng)過P-K濾波后訂正的ZH～Kdp的散點(diǎn)分布與Park曲線更加接近。圖5(c)、(d)分別為訂正前后AH～ZH的散點(diǎn)圖，由于全掃描的數(shù)據(jù)點(diǎn)過多，此處只顯示50個(gè)距離門的數(shù)據(jù)。通過對(duì)比發(fā)現(xiàn)，P-K濾波后訂正的散點(diǎn)圖分布與Park的模擬曲線更加相似。由此可見，訂正后的偏振參量與Park的散射模擬結(jié)果基本一致，進(jìn)一步驗(yàn)證了本訂正方法的有效性。

(a) 訂正前Kdp與ZH的散點(diǎn)圖

4 結(jié)束語

本文根據(jù)雙偏振天氣雷達(dá)偏振參量之間的相互關(guān)系，建立了適當(dāng)?shù)牟罘謧鞑ハ嘁坪筒罘謧鞑ハ嘁坡实臓顟B(tài)方程和過程方程，利用粒子濾波-卡爾曼濾波的綜合算法，首先對(duì)差分傳播相移數(shù)據(jù)進(jìn)行了效果顯著的消噪和估計(jì)過程，并利用自適應(yīng)算法，進(jìn)行了反射率因子的衰減訂正。對(duì)比了衰減前后以及只用卡爾曼濾波算法的訂正結(jié)果。結(jié)果表明，本文采用的方法能夠有效地進(jìn)行衰減訂正，能夠使X波段的雷達(dá)偏振參量訂正到一個(gè)良好的水平，這對(duì)降水量估計(jì)以及降水分類等有重要的作用，所以本方法具有一定的實(shí)際應(yīng)用價(jià)值。