基于粒子濾波的差分傳播相移估計方法

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

0 引言

為了解決氣象災害為我國經濟和生活帶來的嚴重影響,氣象雷達已經廣泛地應用于預防氣象災害、惡劣天氣預報與人工影響天氣等方面[1]。氣象雷達通過發射電磁波探測氣象環境,根據回波的變化來評估氣象目標的特性,當路徑上存在降雨區時,會造成反射率的衰減,為了準確分析氣象目標的真實特性,提高降水估測的精度,需要對反射率進行衰減訂正[2]。雙偏振雷達通過發射水平和垂直極化電磁波不僅能探測到常規的多普勒參量,而且還能獲取表征粒子相態和微物理特性的偏振參量,因此在識別粒子相態、定量估測降水等方面較常規多普勒雷達有很大的優勢[3-4]。對于雙偏振多普勒雷達而言,差分傳播相移率與降雨率之間不僅具有比較高的相關性,而且差分傳播相移還具有不受波束傳播阻礙效應、雷達校準、傳播路徑衰減影響的特性,因此可以使用差分傳播相移與差分傳播相移率進行反射率的衰減訂正[4]。

在實際檢測中,氣象環境的多樣性、雷達系統的噪聲以及由后向散射引起的差分散射相移都會影響差分傳播相移的估計精度[5]。差分傳播相移率是由差分傳播相移估算得到的,因此差分傳播相移率的估算精度受差分傳播相移測量值以及估算方法的影響[6]。當差分傳播相移估計不準時,會影響后續雨衰訂正結果的準確性,與真實的氣象數據不符。因此,對受到污染的差分傳播相移進行準確估計對反射率的衰減訂正尤其重要[1]。Hubbert提出采用有限沖激響應(FIR)和無限沖激響應(IIR)低通濾波器估計差分傳播相移的方法,但在連續多個距離門存在非零的差分散射相移時,該方法并不能有效地對差分傳播相移進行平滑處理,估計效果不好。Hubbert和Bringi[7]提出了迭代濾波方法,通過迭代濾波既可以自動檢測到差分散射相移,還能夠達到剔除干擾的目的,但是迭代次數難以確定,數據處理時間較長。近期國內在差分傳播相移方面也開展了研究,曹俊武等[8]采用了多點平滑的方法處理雷達數據,只能粗略地濾除高頻分量,濾波的效果不明顯。何宇翔等[9]提出卡爾曼濾波方法求取差分傳播相移,該方法可以同步估計差分傳播相移與差分傳播相移率,有效地減小了差分傳播相移的波動,但估計得到的差分傳播相移率存在負值,與實際氣象環境不符[10]。胡志群等[11]提出了小波濾波法,通過小波濾波估計得到的差分傳播相移具有良好的平滑度,并且減少了差分傳播相移率的負值。但是該方法通過變距離法對濾波處理后的差分傳播相移進行最小二乘擬合得到差分傳播相移率,再進行后續的衰減訂正,小波濾波中引入的誤差會傳遞到差分傳播相移率的估計以及反射率的衰減訂正中。

本文提出的粒子濾波方法的估計模型只依賴全差分相移,不受其他偏振參量的約束,而且以雷達偏振參量的不模糊范圍為依據進行采樣,能夠有效地抑制差分傳播相移率的負值,在低信噪比的情況下也能保留真實的氣象信息,使用的條件更加廣泛。該方法首先利用偏振參量之間的關系建立狀態與觀測方程,然后利用粒子濾波方法同步估計差分傳播相移與差分傳播相移率,并以X波段雙偏振多普勒雷達X-SAPR的外場觀測數據作為試驗數據,從差分傳播相移和差分傳播相移率的估計效果與反射率的衰減訂正結果兩方面進行了驗證和分析。

1 粒子濾波的狀態方程和觀測方程的建立

粒子濾波的狀態方程和觀測方程[12-13]表示如下:

式中,xk為系統的狀態向量,T為狀態轉移矩陣,εxk為激勵噪聲,yk為觀測向量,F為觀測矩陣,εyk為觀測噪聲。

下面根據偏振參量之間的關系說明狀態方程的具體形式。為了同步估計差分傳播相移與差分傳播相移率,定義狀態向量xk為

式中,Φdp(k)(k=1,…,K)表示差分傳播相移,Kdp(k)(k=1,…,K)表示差分傳播相移率,為差分傳播相移Φdp(k)(k=1,…,K)隨距離的變化率,k表示沿著傳播路徑電磁波到達的距離門,K表示距離門的個數。將式(3)代入式(1)得到狀態方程為

式中,εxk表示前向傳播路徑上由于氣象環境、雷達系統等引起的不確定性,設定εxk服從正態分布。下面推導狀態轉移矩陣的具體形式,由文獻[14]可知差分傳播相移與差分傳播相移率滿足如下關系:

式中,Δr表示距離門長度。將式(5)代入式(4),當后驗狀態估計Kdp(k)與狀態先驗估計Kdp(k+1)相等[14]時,得到狀態轉移矩陣為

為了避免衰減的偏振參量對估計結果的影響,定義觀測向量為

式中,Ψdp(k)(k=1,…,K)為全差分相移,滿足如下關系:

式中,差分傳播相移Φdp(k)為有用信號,δhv(k)表示由后向散射引起的差分散射相移,為需要分離的高頻噪聲。在經典的估計方法中,認為Kdp具有非負性,因此Φdp的距離廓線不可能出現下降的趨勢[15]。由于不同距離門δhv的變化導致估計Kdp時會存在不合理的負值[16]。為了減少由于δhv產生的估計誤差,將δhv的變化引入到估計模型中。根據文獻[14]中,Hub bert擬合得到的不同頻率的雷達δhv-Kdp的線性關系,可得到c為

式中,b和c的取值依賴于Kdp(k)(k=1,…,K)的取值范圍和雷達的頻率。由式(8)、式(9)相減,可得觀測向量為

由式(3)、式(10)得到觀測方程為

式中,εyk表示觀測引起的誤差,設定εyk服從正態分布。則觀測矩陣為

參數b的選擇依據式(9)中給出的線性擬合關系,c為人為引入用于衡量δhv(k)與b Kdp(k)(k=1,…,K)之間冗余的測量值。根據文獻[14]的方法確定參數b與c的取值。

最后,得到基于粒子濾波估計Φdp與Kdp的狀態方程與觀測方程為

2 基于粒子濾波的差分傳播相移與差分傳播相移率的估計

粒子濾波的思想基于蒙特卡羅實驗方法,旨在通過尋求一組在狀態空間中的隨機樣本對條件后驗概率密度函數進行近似,用樣本均值來替代積分運算,以求得狀態的最小均方誤差估計[13]。

x1∶k={x1,x2,…,xk}是從初始距離門到第k個距離門的狀態集,用表示對第k個距離門的數據進行采樣得到N個粒子,上標i表示采樣得到的第i個粒子。為對x1∶k={x1,x2,…,xk}進行采樣得到的粒子集,y1∶k={y1,y2,…,yk}是從初始距離門到第k個距離門的觀測集。利用最容易得到的狀態轉移概率密度函數作為重要性密度函數,并從中采樣產生粒子。可由下式表示[13]:

則根據式(5)建立的狀態方程進行預測:

可以通過觀測方程迭代更新重要性權值,更新為

權值進行歸一化可得

狀態xk的估計為

采用Smith等[17]提出的多項式采樣方法進行重采樣。根據重要性權值重新采樣得到新的粒子集,并更新粒子的對應權值。最后計算出差分傳播相移與差分傳播相移率的估計值。粒子濾波算法可總結歸納如表1所示。

表1 基于粒子濾波的差分傳播相移估計的算法流程

3 仿真分析

利用ARM(Atmospheric Radiation Measurement Climate Research Facility)的X波段雙偏振多普勒雷達X-SAPR的實測數據驗證算法性能,該雷達在水平和垂直方向同步發射偏振波,差分傳播相移不模糊的范圍為0°～180°[4]。根據文獻[14]得到符合X-SAPR雷達的δhv-Kdp線性關系為

由于Kdp(k)(k=1,…,K)沒有先驗信息,所以b和c必須依賴于Kdp(k)(k=1,…,K)的先驗估計值。根據文獻[9]設定激勵噪聲εxk服從均值為零、方差為10的正態分布,觀測噪聲εyk服從均值為零、方差為2的正態分布。

雷達觀測地點位于緯度36°36′18.0″北、經度97°29′6.0″西。X-SAPR雷達于2013年11月6日探測到大平原南部俄克拉荷馬州地區出現了范圍較大、持續時間較長的降雨過程。選用2013年11月6日1時30分降水過程雷達PPI掃描資料進行分析。

下面通過仿真試驗分析不同濾波方法的估計效果,以下均用“Kalman濾波”來表示何宇翔在文獻[9]中提出的濾波方法。用“滑動平均”表示滑動平均的方法,用“迭代濾波”表示迭代濾波的方法。應用魏慶在文獻[10]定義的FIX參數描述濾波后的性能。

3.1 差分傳播相移濾波效果分析

圖1為X-SAPR雷達于2013年11月6日1時30分1.5°俯仰角、153°方位角Φdp的雷達觀測數據,以及經過不同濾波方法的徑向距離廓線圖。由圖1可知,經過Kalman濾波和粒子濾波處理后的距離廓線的波動和毛刺都得到了很好的抑制,保證了廓線的連續性和平滑度。

服務管理平臺的開發、建設和應用，不僅有效提升了信息中心IT基礎環境的管理水平，也顯著改善了信息中心在開展服務外包時的管理能力，主要功能和作用如下：

圖2為X-SAPR雷達于1.5°仰角觀測數據Ψdp的PPI圖和經過粒子濾波估計Φdp的PPI圖。從圖2(a)可見,由于雷達遠端的信噪比比較低,信號受噪聲影響比較嚴重,導致Ψdp原始數據的PPI圖存在很多波動數據點。圖2(b)為經過粒子濾波處理后的PPI圖,呈現出數據良好的平滑度,有效地剔除了遠端低信噪比區域的干擾以及后向散射相位的影響。

圖1 不同濾波方法處理后的Φdp徑向距離廓線

圖2 2013年11月6日1時30分1.5°俯仰角濾波處理前后Φdp PPI

為了進一步對不同濾波方法的效果進行對比,通過平均波動指數(FIX)來比較距離廓線的波動情況。FIX的定義[10]如下:

FIX越大說明距離廓線的波動就越大。觀測數據Ψdp、滑動平均、迭代濾波、Kalman濾波、粒子濾波的計算結果如表2所示,可見粒子濾波與Kalman濾波都具有一定的濾波效果,使得距離廓線的波動變小,但粒子濾波的波動更小。由此可見,粒子濾波的效果更好。

表2 Φdp徑向距離廓線波動指數統計

3.2 差分傳播相移率的估計分析

圖3為滑動平均、迭代濾波、Kalman濾波和粒子濾波處理后的Kdp徑向距離廓線。結果表明,經過滑動平均、迭代濾波、Kalman濾波和粒子濾波處理后估計的Kdp的負值數量分別為124,92,85和56。說明粒子濾波同步估計Φdp與Kdp的效果比較好,能夠有效地減少Kdp的負值,保留數據的真實信息。

圖3 2013年11月6日1時30分1.5°俯仰角、153°方位角濾波處理后Kdp的距離廓線

3.3 X-SAPR雷達反射率的衰減訂正分析與結果驗證

采用自適應約束算法對反射率Zh進行衰減訂正[4]。

由于Zh在S波段的衰減很小,可以作為真值用來進行Zh訂正前后的對比[6]。S波段雷達KVNX位于緯度36°44′26.9″北、經度98°7′39.0″西,距離庫長為250 m,掃描開始的時間為01:29:41。兩部雷達之間的直線距離為59 km。由于距離雨區的相對距離以及掃描時間的不同,導致X波段雷達與S波段雷達的Zh觀測值會有所偏移,但并不影響Zh訂正效果的驗證。

首先對雷達近端降雨區的衰減訂正效果進行仿真分析。圖4(a)為衰減訂正前后Zh的徑向距離廓線,圖4(b)為衰減訂正前的ZhPPI圖,圖4(c)為同一時段S波段KVNX雷達的ZhPPI圖,圖4(d)～圖4(g)分別為采用滑動平均、迭代濾波、Kalman濾波與粒子濾波處理后進行衰減訂正后的ZhPPI圖。從圖中可以明顯看出,X-SAPR雷達經過迭代濾波、Kalman濾波和粒子濾波處理后訂正的Zh都得到了衰減補償的效果,但圖4(f)中黑色方塊所示的區域Kalman濾波訂正的Zh超過了Zh的真值,出現了過訂正的情況,這也與圖4(a)中Kalman濾波比粒子濾波的Zh的取值高出2～8 d B相對應,因此經過粒子濾波處理訂正后的Zh與Zh的真值更加接近。

通過Park由散射模擬建立的偏振參量的經驗關系驗證衰減訂正的效果[4],比較了X波段訂正前后的Ah～Zh和Zh～Kdp之間的散點圖特性,Ah表示水平方向的衰減率。圖6(a)和圖6(b)分別為訂正前后的Zh～Kdp的散點圖,實線為Park通過散射模擬建立的Zh～Kdp的經驗關系。由圖6(a)可以發現,訂正前的散點圖比較分散,Zh大約分布在10～30 dBz,Kdp分布在0～6°/km,與Park的模擬曲線有很大偏移。經過訂正,Zh～Kdp的散點分布與Park曲線比較接近。圖6(c)和圖6(d)分別為訂正前后Ah～Zh的散點圖,實線是Park依據公式Ah=a Zβh經過散射模擬得到的曲線。通過對比發現,訂正后的散點圖分布與Park的模擬曲線比較相似,而訂正前的偏移較大。由此可見,訂正后的偏振參量與Park的散射模擬結果基本一致,進一步驗證了本文估計方法的有效性。

圖4 雷達近端降雨區訂正前后Zh的比較

圖5 雷達遠端降雨區訂正前后Zh的比較

圖6 訂正前后的偏振參量散點圖分析

4 結束語

本文提出了基于粒子濾波的X波段雙偏振氣象雷達差分傳播相移與差分傳播相移率估計的新方法,利用ARM的X波段雙偏振雷達X-SAPR的實測數據驗證了算法的性能。該方法在低信噪比的情況下能夠有效剔除Φdp存在的波動數據點和毛刺,使數據具有良好的濾波效果,Φdp距離廓線體現出較好的收斂性,更加符合實際降水過程距離廓線的變化。其次該方法能夠有效地保持Kdp的非負性,保留數據的真實信息。最后X-SAPR雷達經過粒子濾波處理后訂正的Zh有了明顯的增強,強回波位置Zh值與S波段KVNX雷達更加接近。通過分析訂正前后X波段雙偏振雷達參量之間的散點圖發現,經過訂正后的散點圖與偏振參量之間的經驗公式具有更強的一致性,證明了本文方法的有效性。

[1]LIM S,CHANDRASEKAR V.A Robust Attenuation Correction System for Reflectivity and Differential Reflectivity in Weather Radars[J].IEEE Trans on Geoscience and Remote Sensing,2016,54(3):1727-1737.

[2]曹俊武,胡志群.X波段測雨雷達強度資料評估及改進方法[J].雷達科學與技術,2016,14(3):237-243.CAO Junwu,HU Zhiqun.Evaluation and Its Improved Method for Reflectivity Data Quality of Rain Radar[J].Radar Science and Technology,2016,14(3):237-243.(in Chinese)

[3]宋新景.基于極化特征的雷達目標識別技術[J].雷達科學與技術,2016,14(1):39-44,53.SONG Xinjing.Radar Target Recognition Based on Polarization Feature[J].Radar Science and Technology,2016,14(1):39-44,53.(in Chinese)

[4]BRINGI V N,KEENAN T D,CHANDRASEKAR V.Correcting C-Band Radar Reflectivity and Differential Reflectivity Data for Rain Attenuation:A Self-Consistent Method with Constraints[J].IEEE Trans on Geoscience and Remote Sensing,2001,39(9):1906-1915.

[5]汪旭東,胡志群,劉浩,等.雙線偏振天氣雷達性能測試方法研究[J].雷達科學與技術,2015,13(4):395-401,409.WANG Xudong,HU Zhiqun,LIU Hao,et al.Research on Measurement Method for Dual-Polarization Weather Radar[J].Radar Science and Technology,2015,13(4):395-401,409.(in Chinese)

[6]魏慶,胡志群,劉黎平,等.C波段偏振雷達數據預處理及在降水估計中的應用[J].高原氣象,2016,35(1):231-243.

[7]HUBBERT J,BRINGI V N.An Iterative Filtering Technique for the Analysis of Copolar Differential Phase and Dual-Frequency Radar Measurements[J].Journal of Atmospheric and Oceanic Technology,1995,12(3):643-648.

[8]曹俊武,胡志群,陳曉輝,等.影響雙線偏振雷達相位探測精度的分析[J].高原氣象,2011,30(3):817-822.

[9]何宇翔,呂達仁,肖輝,等.X波段雙線極化雷達反射率的衰減訂正[J].大氣科學,2009,33(5):1027-1037.

[10]魏慶,胡志群,劉黎平.雙偏振雷達差分傳播相移的五種濾波方法對比分析[J].成都信息工程學院學報,2014,29(6):596-602.

[11]杜牧云,劉黎平,胡志群,等.雙線偏振雷達差分傳播相移的小波濾波初探[J].暴雨災害,2012,31(3):248-254.

[12]ALA-LUHTALA J,WHITELEY N,HEINE K,et al.An Introduction to Twisted Particle Filters and Parameter Estimation in Non-Linear State-Space Models[J].IEEE Trans on Signal Processing,2016,64(18):4875-4890.

[13]李天成,范紅旗,孫樹棟.粒子濾波理論、方法及其在多目標跟蹤中的應用[J].自動化學報,2015,41(12):1981-2002.

[14]SCHNEEBELI M,GRAZIOLI J,BERNE A.Improved Estimation of the Specific Differential Phase Shift Using a Compilation of Kalman Filter Ensembles[J].IEEE Trans on Geoscience and Remote Sensing,2014,52(8):5137-5149.

[15]HUANG H,ZHANG G,ZHAO K,et al.A Hybrid Method to Estimate Specific Differential Phase and Rainfall with Linear Programming and Physics Constraints[J].IEEE Trans on Geoscience and Remote Sensing,2017,55(1):96-111.

[16]曹楊,蘇德斌,周筠珺,等.C波段雙線偏振多普勒雷達差分相位質量分析[J].高原氣象,2016,35(2):548-559.

[17]SMITH A F M,GELFAND A E.Bayesian Statistics Without Tears:a Sampling-Resampling Perspective[J].The American Statistician,1992,46(2):84-88.