基于粒子滤波的差分传播相移估计方法

(中国民航大学天津市智能信号与图像处理重点实验室,天津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.