李 海, 羅 原, 馮興寰, 馮 青
(中國民航大學天津市智能信號與圖像處理重點實驗室,天津 300300)
相比單偏振雷達,雙偏振雷達通過發射水平和垂直偏振電磁波不僅能探測到常規的多普勒參量,還能獲取表征粒子相態和微物理特性的其他偏振參量。天氣雷達的偏振參量對于確定雷達校準、地面雜波的消除、降水分類和雨滴譜參數估計等方面具有更大的優勢[1]。相比C波段和S波段天氣雷達,X波段天氣雷達對弱氣象目標的探測更敏銳,并且具有天線尺寸小、易于移動等特點。但由于降水區域會對天氣雷達的回波信號產生散射干擾以及吸收,會導致接收到的數據與實際值相比會有衰減,特別是波長較短的X波段(中心波長3 cm),其回波信號的衰減現象更為嚴重,這也是X波段天氣雷達應用推廣所面臨的主要問題。
X波段天氣雷達雨衰補償相比C波段和S波段天氣雷達的研究和發展較晚,因而對X波段天氣雷達雨衰補償研究的最開始階段是借鑒了C波段和S波段氣象雷達的相關方法以及研究結果[2]。2008年胡志群等人提出了差分傳播相移率-衰減率(Kdp-Adp)方法,并通過對Kdp設置閾值,對水平反射率因子(ZH)進行衰減訂正,該方法有一定的效果,但由于其沒有考慮到Kdp的累積量在電磁波傳播路徑上會受到前向差分散射相移的影響,因而具有一定的局限性。2009年何宇翔等人引入卡爾曼濾波首先對差分傳播相移(Φdp)進行濾波處理,進而對ZH進行衰減訂正,該方法有效地減小了水平反射率因子的誤差,但其處理過程是將降水區域視為線性關系,其具有一定的局限性。2014年魏慶等人采用線性滑窗、小波分析等方法對Φdp進行了衰減訂正;2015年,孫躍、肖輝等人提出了線性擬合-遞推法(LFRM)對Φdp進行質量控制。這些方法都能對Φdp進行一定程度的衰減訂正,但這些方法的共同點都是將降水區域視為一個連續的線性區域,但由于天氣雷達的回波功率只代表某一距離門的單位波束體積內N個降水粒子的回波功率平均值,不能理想的將天氣雷達不同距離門的回波信號數據視為線性關系[3],因而這些方法具有一定的局限性。
針對X波段天氣雷達的偏振參量衰減訂正,本文首先采用粒子濾波算法對Φdp進行濾波處理,在經過粒子濾波算法處理后,對數據利用卡爾曼濾波進行最優化估計,以達到對Φdp數據的精確估計訂正。在得到Φdp濾波數據后,將Φdp數據作為參照量采用自適應算法,對ZH進行衰減訂正處理。天氣雷達的偏振參量回波數據是一個距離門內所有采樣點的疊加值,常規的窗函數濾波等方法對所包含的噪聲信號的濾除效果有限。粒子濾波能夠利用偏振參量之間的關系建立合適的狀態方程和觀測方程,可以有效地進行數據的濾波處理,但由于粒子濾波的時間復雜度隨著粒子數的增加而增加,于是引入卡爾曼濾波,能夠對經過較少粒子數的粒子濾波處理后數據進行最優化估計,從而達到時間復雜度的平衡[4]。采用P-K算法能夠得到更平滑和準確的Φdp,進而對水平反射率因子進行更精確的衰減訂正。
X波段雙偏振雷達回波數據包含有ZH,Kdp,Φdp等偏振參量信息。在對ZH進行衰減訂正之前需要先對Φdp進行衰減訂正。
直接從天氣雷達回波數據中提取到的并不是差分傳播相移而是包含了前向差分散射相移和后向差分散射相移的全差分相移(Ψdp),其中后向差分散射相移也就是差分傳播相移(Φdp)。單個距離門的全差分相移(Ψdp(k))定義為
Ψdp(k)=Φdp(k)+δhv(k),k=1,…,K
(1)
式中:k表示傳播路徑電磁波到達的距離門,總共有K個距離單元;Φdp(k)表示第k個距離門的差分傳播相移;δhv(k)表示第k個距離門的前向差分散射相移,它是單基地天氣雷達回波信號中的有害高頻噪聲成分。X波段電磁波在小雨到中雨區域偏離瑞利散射很小,δhv(k)可以忽略,但對于較強的雨區,就可能產生較強的瑞利散射,因而不能被忽略,所以雨衰補償的重點就是從全差分相移中將δhv(k)濾除。
從全差分相移中無法直接通過濾波等方式將前向差分散射相移濾除,Hubbert擬合得到的不同頻率的雷達的δhv(k)與單個距離門的差分傳播相移率Kdp(k)有如下關系[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)之間的關系:
Ψdp(k)-c=Φdp(k)+bKdp(k),k=1,…,K
(3)
而式(1)中的差分傳播相移與差分傳播相移率之間又有如下關系[7]:
Φdp(k+1)=Φdp(k)+2Δr·Kdp(k),
k=1,…,K
(4)
式中,Δr表示相鄰距離門之間的距離。
根據式(3)和式(4)即可建立如下的粒子濾波方程組[8]:

k=1,…,K
(5)

k=1,…,K
(6)

(7)
以觀測量與預測狀態之間的差值作為似然函數,計算公式如式(8)表示:
(8)
可以通過觀測方程迭代更新重要性權值,更新步為
(9)
權值進行歸一化可得
(10)
進而可以得到狀態xk的估計為
(11)
即
(12)
粒子濾波通過尋找一組在狀態空間的隨機樣本對概率密度函數進行近似,以后驗概率密度所產生的N個獨立同分布樣本的均值代替積分運算,從而獲得狀態最小方差分布,但由于粒子數的增加會導致其運算量急劇增加,而較少的粒子又不能得到最優的最小方差狀態值,因而建立卡爾曼濾波對數據進行再處理,從而得到結果值。
(13)
(14)

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

更新過程是結合先驗狀態估計與當前距離門觀測向量的數據進行后驗估計過程,更新過程為
(17)
進而得到卡爾曼濾波估計方程為
(18)
而協方差更新方程為
(19)

(20)
對差分傳播相移的衰減訂正是對水平反射率因子訂正的基礎,在利用P-K濾波算法得到誤差更小的差分傳播相移數據后,利用自適應算法,對反射率因子(ZH)進行衰減訂正。
利用Bringi提出的自適應約束(self-consistent method with constraints)算法對ZH進行衰減訂正[11]。ZH的衰減訂正的關鍵就是確定雨區的衰減率AH,設定雨區的距離門范圍為k0~k1,可得到第k個距離門的衰減率AH(k)的計算公式為[12]
AH(k)=
(21)

(22)

(23)
本仿真實驗的X波段氣象回波數據來源于大氣輻射測量氣候研究中心(ARM)的網站。其數據包括了Φdp,Kdp,ZH等重要回波參數,每個徑向上有400個距離單元,每個距離單元代表100 m。實驗選用的X波段雷達于2018年7月3日11:00的回波數據。
一個徑向上的Φdp數據經過不同方法的濾波效果與源數據效果對比如圖1所示。

圖1 同一徑向不同方法的Φdp濾波效果圖
由圖1可以看出,相對比原始數據,卡爾曼濾波以及P-K濾波均能對反射率因子的波動和高頻噪聲進行抑制,保證了廓線的連續性和平滑度。但相比卡爾曼濾波,P-K濾波能夠對數據進行更好的平滑處理。
如圖2所示,為一個俯仰角的全徑向PPI圖對比,與卡爾曼濾波結果相比較,P-K濾波能更有效地將數據中的高頻噪聲部分剔除,且數據能呈現較好的非負性,更加符合真實氣象特征[14]。

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

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

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

(a) 訂正前Kdp與ZH的散點圖
本文根據雙偏振天氣雷達偏振參量之間的相互關系,建立了適當的差分傳播相移和差分傳播相移率的狀態方程和過程方程,利用粒子濾波-卡爾曼濾波的綜合算法,首先對差分傳播相移數據進行了效果顯著的消噪和估計過程,并利用自適應算法,進行了反射率因子的衰減訂正。對比了衰減前后以及只用卡爾曼濾波算法的訂正結果。結果表明,本文采用的方法能夠有效地進行衰減訂正,能夠使X波段的雷達偏振參量訂正到一個良好的水平,這對降水量估計以及降水分類等有重要的作用,所以本方法具有一定的實際應用價值。