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

(中國民航大學天津市智能信號與圖像處理重點實驗室，天津 300300)

0 引 言

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

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

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

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

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

1.1 建立粒子濾波方程

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

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

(1)

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

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

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

(2)

式中，b和c的取值依賴于Kdp(k)(k=1,…，K)的取值[6]，當Kdp(k)≤2.5°/km時，b取2.37，c取0.054；當Kdp(k)>2.5°/km時，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)

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

(8)

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

(9)

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

(10)

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

(11)

即

(12)

1.2 建立卡爾曼濾波方程

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

(13)

(14)

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

(15)

(16)

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

(17)

進而得到卡爾曼濾波估計方程為

(18)

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

(19)

(20)

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

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

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

AH(k)=

(21)

(22)

(23)

3 實驗結(jié)果及分析

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

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

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

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

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

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

(a) 卡爾曼濾波

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

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

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

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

(a) S波段ZHPPI掃描圖

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

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

4 結(jié)束語

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