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

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

0 引言

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

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

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

1 粒子濾波的狀態(tài)方程和觀測方程的建立

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

式中,xk為系統(tǒng)的狀態(tài)向量,T為狀態(tài)轉(zhuǎn)移矩陣,εxk為激勵噪聲,yk為觀測向量,F為觀測矩陣,εyk為觀測噪聲。

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

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

式中,εxk表示前向傳播路徑上由于氣象環(huán)境、雷達系統(tǒng)等引起的不確定性,設(shè)定εxk服從正態(tài)分布。下面推導(dǎo)狀態(tài)轉(zhuǎn)移矩陣的具體形式,由文獻[14]可知差分傳播相移與差分傳播相移率滿足如下關(guān)系:

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

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

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

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

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

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

式中,εyk表示觀測引起的誤差,設(shè)定εyk服從正態(tài)分布。則觀測矩陣為

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

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

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

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

x1∶k={x1,x2,…,xk}是從初始距離門到第k個距離門的狀態(tài)集,用表示對第k個距離門的數(shù)據(jù)進行采樣得到N個粒子,上標i表示采樣得到的第i個粒子。為對x1∶k={x1,x2,…,xk}進行采樣得到的粒子集,y1∶k={y1,y2,…,yk}是從初始距離門到第k個距離門的觀測集。利用最容易得到的狀態(tài)轉(zhuǎn)移概率密度函數(shù)作為重要性密度函數(shù),并從中采樣產(chǎn)生粒子?？捎上率奖硎綶13]:

則根據(jù)式(5)建立的狀態(tài)方程進行預(yù)測:

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

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

狀態(tài)xk的估計為

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

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

3 仿真分析

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

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

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

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

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

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

服務(wù)管理平臺的開發(fā)、建設(shè)和應(yīng)用，不僅有效提升了信息中心IT基礎(chǔ)環(huán)境的管理水平，也顯著改善了信息中心在開展服務(wù)外包時的管理能力，主要功能和作用如下：

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

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

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

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

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

表2 Φdp徑向距離廓線波動指數(shù)統(tǒng)計

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

[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].雷達科學(xué)與技術(shù),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]宋新景.基于極化特征的雷達目標識別技術(shù)[J].雷達科學(xué)與技術(shù),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].雷達科學(xué)與技術(shù),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波段偏振雷達數(shù)據(jù)預(yù)處理及在降水估計中的應(yīng)用[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].大氣科學(xué),2009,33(5):1027-1037.

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

[11]杜牧云,劉黎平,胡志群,等.雙線偏振雷達差分傳播相移的小波濾波初探[J].暴雨災(zāi)害,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]李天成,范紅旗,孫樹棟.粒子濾波理論、方法及其在多目標跟蹤中的應(yīng)用[J].自動化學(xué)報,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波段雙線偏振多普勒雷達差分相位質(zhì)量分析[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.