基于粒子滤波的双偏振雷达差分传播相移的估计方法

摘要

一种基于粒子滤波的双偏振雷达差分传播相移的估计方法。其包括利用雷达偏振参量间相互关系，建立状态和观测方程；在单个距离门中，进行初始化采样，利用状态方程计算重要性密度函数，将初始化采样数据和重要性密度函数结合得到状态预测值；在单个距离门中，利用观测方程求取似然函数，并实现重要性权值的迭代更新；将全部粒子的状态向量与重要性权值求取均值，由此实现对差分传播相移与差分传播相移率的估计。本发明方法不仅可准确估计差分传播相移，使滤波处理后的数据具有更好的连续性、平滑性和准确性，且由于以雷达偏振参量的不模糊范围为依据进行采样，还可以有效地抑制差分传播相移率的负值，保留真实的气象数据，使用的条件更加广泛。

说明书

技术领域

本发明属于气象雷达信号处理技术领域，特别是涉及一种基于粒子滤波的双偏振雷达差分传播相移的估计方法。

背景技术

为了解决气象灾害为我国经济和生活带来的严重影响，气象雷达已经广泛地应用于预防气象灾害、恶劣天气预报与人工影响天气等方面。气象雷达通过发射电磁波探测气象环境，根据回波的变化来评估气象目标的特性，当路径上存在降雨区时，会造成反射率的衰减，为了准确分析气象目标的真实特性，提高降水估测的精度，需要对反射率进行衰减订正。双偏振雷达通过发射水平和垂直极化电磁波不仅能探测到常规的多普勒参量，而且还能获取表征粒子相态和微物理特性的偏振参量，因此在识别粒子相态、定量估测降水等方面较常规多普勒雷达有很大的优势。对于双偏振多普勒雷达而言，差分传播相移率与降雨率之间不仅具有比较高的相关性，而且差分传播相移还具有不受波束传播阻碍效应、雷达校准、传播路径衰减影响的特性，因此可以使用差分传播相移与差分传播相移率进行反射率的衰减订正。

在实际检测中，气象环境的多样性、雷达系统的噪声以及由后向散射引起的差分散射相移都会影响差分传播相移的估计精度。差分传播相移率是由差分传播相移估算得到的，因此差分传播相移率的估算精度受差分传播相移测量值以及估算方法的影响。当差分传播相移估计不准时，会影响后续的雨衰订正结果的准确性，与真实的气象数据不符。因此，对受到污染的差分传播相移进行准确的估计对反射率的衰减订正尤其重要。

传统低通滤波器估计差分传播相移的方法是当连续多个距离门存在非零的差分散射相移时，并不能有效地对差分传播相移进行平滑处理，估计效果不好。通过迭代滤波既可以自动检测到差分散射相移，还能够达到剔除干扰的目的，但是迭代次数难以确定，因此数据处理时间较长。近期国内也开展了在差分传播相移方面的研究，引入卡尔曼滤波方法求取差分传播相移，该方法可以同步估计差分传播相移与差分传播相移率，有效地减小了差分传播相移的波动，但估计得到的差分传播相移率存在负值，与实际气象环境不符。通过小波分析估计得到的差分传播相移具有良好的平滑度，并且减少了差分传播相移率的负值，但是该方法在强降雨区容易受到衰减的偏振参量的影响，使估计结果不准确。综上，上述原因制约了差分传播相移在气象雷达信号处理中的应用与推广。

发明内容

为了解决上述问题，本发明的目的在于提供一种能够保证偏振参量的估计精度，同时在激励噪声的先验信息未知的情况下也能保留真实的气象信息的基于粒子滤波的双偏振雷达差分传播相移的估计方法。

为了达到上述目的，本发明提供的基于粒子滤波的双偏振雷达差分传播相移的估计方法包括按顺序进行的下列步骤：

1)利用雷达偏振参量之间的相互关系，建立状态方程和观测方程；

2)在单个距离门中，根据雷达偏振参量的不模糊范围进行初始化采样,然后利用上述状态方程计算重要性密度函数，之后将初始化采样数据和重要性密度函数结合起来共同完成状态预测，得到状态预测值；

3)在单个距离门中，利用步骤1)中得到的观测方程求取似然函数，并实现重要性权值的迭代更新；

4)根据步骤3)获得的重要性权值判定是否符合继续迭代的要求，若不符合要求，进行重采样并重复步骤3)，重新进行状态预测和重要性权值更新，否则进入下一步；

5)将符合要求的全部粒子的状态向量与重要性权值求取均值，由此实现对差分传播相移与差分传播相移率的估计。

在步骤1)中，所述的利用雷达偏振参量之间的相互关系，建立状态方程和观测方程的方法是：首先，将差分传播相移与差分传播相移率整合为状态向量以实现同步估计；然后，分析相邻距离门之间差分传播相移与差分传播相移率的相互关系，确定状态转移矩阵的具体形式；接着，定义观测向量为全差分相位，以避免衰减的偏振参量对估计结果的影响；最后，为了减少由于后向散射产生的估计误差，将后向散射的变化引入到观测方程中，最终确定状态方程和观测方程的具体形式。

在步骤2)中，所述的在单个距离门中，根据雷达偏振参量的不模糊范围进行初始化采样,然后利用上述状态方程计算重要性密度函数，之后将初始化采样数据和重要性密度函数结合起来共同完成状态预测，得到状态预测值的方法是：首先，根据已知的雷达偏振参量的不模糊范围进行初始化采样而得到；然后，根据步骤1)中建立的状态方程获得状态转移密度函数，并将其作为重要性密度函数；最后，根据重要性密度函数产生新的采样粒子，并与初始化采样数据结合得到状态预测值，由此实现状态的预测。

在步骤3)中，所述的在单个距离门中，利用步骤1)中得到的观测方程求取似然函数，并实现重要性权值的迭代更新的方法是：首先，利用步骤1)中得到的观测方程计算观测向量与状态预测值之间的冗余，并将其作为似然函数；然后，利用似然函数计算新的重要性权值；最后，循环处理当前距离门内的全部粒子，实现重要性权值的更新。

在步骤4)中，所述的根据步骤3)获得的重要性权值判定是否符合继续迭代的要求，若不符合要求，进行重采样并重复步骤3)，重新进行状态预测和重要性权值更新的方法是：为避免粒子退化造成的预测误差，利用步骤3)得到的重要性权值计算有效粒子数，当有效粒子数小于设定的阈值时，则进行重采样，重新进行状态预测并更新重要性权值，直至重要性权值符合继续迭代的要求，即有效粒子数大于设定的阈值时进行接下来的差分传播相移与差分传播相移率估计。

在步骤5)中，所述的将符合要求的全部粒子的状态向量与重要性权值求取均值，由此实现对差分传播相移与差分传播相移率的估计的方法是：将符合要求的全部粒子与更新得到的重要性权值相结合，求取各自的均值作为最终的状态估计值；最后，对全部距离门进行循环处理，最终计算出每个距离门的差分传播相移与差分传播相移率的估计值。

本发明提供的基于粒子滤波的双偏振雷达差分传播相移估计方法由于估计模型中的观测向量只依赖全差分相移，不受其他偏振参量的约束，而且以雷达偏振参量的不模糊范围为依据进行采样，可以有效地抑制差分传播相移率的负值，并在激励噪声先验信息未知的情况下也能保留真实的气象数据。本发明方法采用粒子滤波，利用观测到的偏振参量之间的关系建立状态与观测方程，并将方程应用于同步估计差分传播相移与差分传播相移率，还将X波段双偏振多普勒雷达X-SAPR的外场观测数据作为实验数据，从差分传播相移和差分传播相移率的估计效果与经过滤波处理后的反射率的衰减订正结果两方面验证了本方法的有效性。

附图说明

图1为本发明提供的基于粒子滤波的双偏振雷达差分传播相移估计方法流程图。

图2为X-SAPR雷达于2013年11月6日01:30 1.5°俯仰角、153°方位角差分传播相移的雷达观测数据以及经过不同滤波方法的径向距离廓线图。

图3为X-SAPR雷达于1.5°仰角观测数据全差分相位的PPI图和经过粒子滤波估计差分传播相移的PPI图。

图4为Kalman滤波和粒子滤波处理后的差分传播相移率径向距离廓线。

图5为衰减订正前后的反射率的径向距离廓线与PPI图。

图6为衰减订正前后的偏振参量散点图分析。

具体实施方法

下面结合附图和具体实施例对本发明提供的基于粒子滤波的双偏振雷达差分传播相移估计方法进行详细说明。

如图1所示，本发明提供的基于粒子滤波的双偏振雷达差分传播相移估计方法包括按顺序进行的下列步骤：

1)利用雷达偏振参量之间的相互关系，建立状态方程和观测方程；

基于粒子滤波的状态方程和观测方程分别表示如下：

其中x_k表示状态向量，T表示状态转移矩阵，表示激励噪声，y_k表示观测向量，F表示观测矩阵，表示观测噪声。

下面根据雷达偏振参量之间的关系说明上述状态方程的具体形式。为了同步估计差分传播相移与差分传播相移率，本发明定义状态向量x_k为：

其中，Φ_dp(k),k＝1,2,…,K表示差分传播相移，K_dp(k),k＝1,2,…,K表示差分传播相移率，为差分传播相移Φ_dp(k),k＝1,2,…,K随距离的变化率，k表示沿着传播路径电磁波到达的距离门，K为距离门的个数。将式(3)带入式(1)得到状态方程为：

表示激励噪声，为前向传播路径上由于气象环境、雷达系统等因素引起的不确定性参数，设定激励噪声服从正态分布。下面推导状态转移矩阵的具体形式，差分传播相移与差分传播相移率满足如下关系：

Φ_dp(k+1)＝Φ_dp(k)+2△rK_dp(k)(5)

其中，△r表示距离门长度。将式(5)带入式(4)，当后验状态估计的差分传播相移率K_dp(k)与先验状态估计的差分传播相移率K_dp(k+1)相等时，得到状态转移矩阵为：

为了避免衰减的雷达偏振参量对估计结果的影响，定义观测向量为：

y_k＝[Ψ_dp(k)-c](7)

其中Ψ_dp(k),k＝1,2,…,K为全差分相移，满足如下关系：

Ψ_dp(k)＝Φ_dp(k)+δ_hv(k)>

其中，差分传播相移Φ_dp(k),k＝1,2,…,K为有用信号，δ_hv(k)表示由后向散射引起的差分散射相移，为需要分离的高频噪声。在经典的估计方法中，认为差分传播相移率K_dp具有非负性，因此差分传播相移Φ_dp的距离廓线不可能出现下降的趋势。由于不同距离门的差分散射相移δ_hv的变化导致估计差分传播相移率K_dp时会存在不合理的负值，为了减少由于距离门的差分散射相移δ_hv产生的估计误差，将距离门的差分散射相移δ_hv的变化引入到观测方程中。根据Hubbert拟合得到的不同频率的雷达δ_hv-bK_dp的线性关系，可得到参数c为：

c＝δ_hv(k)-bK_dp(k)>

其中，参数b和c的取值依赖于差分传播相移率K_dp(k),k＝1,2,…,K的取值范围和雷达的频率。由式(8)、(9)相减，可得观测向量为：

y_k＝[Ψ_dp(k)-c]＝[Φ_dp(k)+bK_dp(k)]>

由式(2)、式(3)、式(10)得到观测方程为：

表示观测噪声，设定观测噪声服从正态分布。则观测矩阵为：

F＝[1 b](12)

参数b的选择依据式(9)中给出的线性拟合关系，参数c为人为引入用于衡量δ_hv(k),k＝1,2,…,K与bK_dp(k),k＝1,2,…,K之间冗余的测量值。

最后，得到基于粒子滤波估计的差分传播相移Φ_dp与差分传播相移率K_dp的状态方程与观测方程为：

2)在单个距离门中，根据雷达偏振参量的不模糊范围进行初始化采样,然后利用上述状态方程计算重要性密度函数，之后将初始化采样数据和重要性密度函数结合起来共同完成状态预测，得到状态预测值；

在本发明中，将雷达偏振参量的不模糊范围作为先验信息，进行初始化采样。

x_1:k＝{x₁,x₂,…,x_k}是从初始距离门到第k个距离门的状态集，用表示对第k个距离门的数据进行采样得到N个粒子，上标i表示采样得到的第i个粒子。为对x_1:k＝{x₁,x₂,…,x_k}进行采样得到的粒子集，y_1:k＝{y₁,y₂,…,y_k}是从初始距离门到第k个距离门的观测集。利用最容易得到的状态转移概率密度函数并作为重要性密度函数，并从中采样而产生新的采样粒子。可由下式表示：

则根据式(4)建立的状态方程进行状态预测而获得状态预测值：

3)在单个距离门中，利用步骤1)中得到的观测方程求取似然函数，并实现重要性权值的迭代更新；

通过重要性权值对状态进行更新，似然函数表征每个粒子的重要性权值的确定，由步骤1)得到的观测方程确定，可以用观测向量与状态预测值之间的冗余来表示：

利用似然函数计算新的重要性权值，然后循环处理当前距离门内的全部粒子，以实现重要性权值的迭代更新，更新公式为：

4)根据步骤3)获得的重要性权值判定是否符合继续迭代的要求，若不符合要求，进行重采样并重复步骤3)，重新进行状态预测和重要性权值更新，否则进入下一步；

在以上过程中，由于会经过多次迭代，有可能出现只有少量粒子具有较大的重要性权值的现象，称之为粒子退化现象，该现象会使得之后的预测结果产生较大误差，即所得到的重要性权值不能满足继续迭代的要求。为避免这一现象的影响，需要进行重采样，进行重采样的条件如下：

定义有效粒子数N_eff为设定阈值N_th＝0.5N。当N_eff<N_th时进行重采样，重新进行状态预测并更新重要性权值，直至N_eff≥N_th时，重要性权值符合继续迭代的要求，进行接下来的差分传播相移与差分传播相移率估计。

5)将符合要求的全部粒子的状态向量与重要性权值求取均值，由此实现对差分传播相移与差分传播相移率的估计。

对重要性权值进行归一化可得：

状态向量x_k的均值为：

即差分传播相移与差分传播相移率的估计值为：

最后，对全部距离门进行循环处理，最终计算出每个距离门的差分传播相移与差分传播相移率的估计值。

本发明提供的基于粒子滤波的双偏振雷达差分传播相移估计方法的效果可以通过以下仿真结果进一步说明。

仿真参数设置：利用ARM(Atmospheric Radiation Measurement ClimateResearch Facility)的X波段双偏振多普勒雷达X-SAPR的实测数据验证本发明方法的性能，该雷达在水平和垂直方向同步发射偏振波，差分传播相移不模糊的范围为0～180°。X-SAPR雷达的δ_hv—K_dp线性关系为：

由于K_dp(k),k＝1,2,…,K没有先验信息，所以参数^b和^c必须依赖于K_dp(k),k＝1,2,…,K的先验估计值。设定激励噪声服从均值为零的正态分布，观测噪声服从均值为零，方差为2的正态分布。

雷达观测地点位于纬度36°36'18.0"N，经度97°29'6.0"W。X-SAPR雷达于2013年11月6日探测到大平原南部俄克拉荷马州地区出现了范围较大、持续时间较长的降雨过程。选用11月6日01:30降水过程雷达PPI扫描资料进行分析。

图2为X-SAPR雷达于2013年11月6日01:301.5°俯仰角153°方位角Φ_dp的雷达观测数据以及经过不同滤波方法的径向距离廓线图。由图2可知经过Kalman滤波和粒子滤波处理后的距离廓线的波动和毛刺都得到了很好的抑制，保证了廓线的连续性和平滑度。

图3为X-SAPR雷达于1.5°仰角的观测数据Ψ_dp的PPI图和经过粒子滤波估计的差分传播相移Φ_dp的PPI图。从图3a中可见，由于雷达远端的信噪比比较低，信号受噪声影响比较严重，导致观测数据Ψ_dp的PPI图存在很多波动数据点。图3b为经过粒子滤波处理后的PPI图，呈现出数据良好的平滑度，有效地剔除了远端低信噪比区域的干扰以及后向散射相位的影响。

为了进一步对Kalman滤波和粒子滤波的滤波效果进行对比，通过平均波动指数FIX来比较距离廓线的波动情况。FIX的定义如下

FIX越大说明距离廓线的波动就越大。观测数据Ψ_dp、Kalman滤波、粒子滤波的计算结果如表1所示，可见粒子滤波与Kalman滤波都具有一定的滤波效果，使得距离廓线的波动变小，但粒子滤波的波动更小，由此可见，粒子滤波的效果更好。

图4为Kalman滤波和粒子滤波处理后的差分传播相移率K_dp径向距离廓线。结果表明，经过Kalman滤波和粒子滤波处理后估计的差分传播相移率K_dp的负值数量分别为85和56。说明粒子滤波同步估计差分传播相移Φ_dp与差分传播相移率K_dp的效果比较好，能够有效地减少差分传播相移率K_dp的负值，保留数据的真实信息。

通过将估计结果应用于衰减订正进一步对本发明方法的有效性进行验证和分析。本发明方法采用自适应约束(self-consistent method with constraints)算法对反射率Z_h进行衰减订正。由于反射率Z_h在S波段的衰减很小，可以作为真值用来进行反射率Z_h订正前后的对比。S波段雷达KVNX位于纬度36°44'26.9"N，经度98°7'39.0"W，距离库长为250m，扫描开始的时间为01:29:41。两部雷达之间的直线距离为59km。由于距离雨区的相对距离以及扫描时间的不同，导致X波段雷达与S波段雷达的反射率Z_h观测值会有所偏移，但并不影响反射率Z_h订正效果的验证。

图5a为衰减订正前后的反射率Z_h的径向距离廓线。图5b是衰减订正前的反射率Z_hPPI图，图5c是同一时段S波段KVNX雷达的Z_hPPI图，图5d、5e是应用Kalman滤波与粒子滤波处理后进行衰减订正后的反射率Z_hPPI图。从图中可以明显看出，X-SAPR雷达经过Kalman滤波和粒子滤波处理后订正的反射率Z_h都得到了衰减补偿的效果，但图5d中黑色方块所示的区域Kalman滤波订正的反射率Z_h超过了反射率Z_h的真值，出现了过订正的情况，这也与图5a中Kalman滤波比粒子滤波的反射率Z_h的取值高出2～8dB，相对应，因此经过粒子滤波处理订正后的反射率Z_h与反射率Z_h的真值更加接近。

通过Park由散射模拟建立的偏振参量的经验关系验证衰减订正的效果，比较了X波段订正前后的A_h～Z_h和Z_h～K_dp之间的散点图特性，A_h表示水平方向的衰减率。图6a、6b分别为订正前后的Z_h～K_dp的散点图，实线为Park通过散射模拟建立的Z_h～K_dp的经验关系，由图6a可以发现，订正前的散点图比较分散，反射率Z_h大约分布在10～30dBZ，差分传播相移率K_dp分布在0～6°/km，与Park的模拟曲线有很大偏移；经过订正，Z_h～K_dp的散点分布与Park曲线比较接近。图6c、6d分别为订正前后A_h～Z_h的散点图，实线是Park依据公式经过散射模拟得到的曲线。通过对比发现，订正后的散点图分布与Park的模拟曲线比较相似，而订正前的偏移较大。由此可见，订正后的偏振参量与Park的散射模拟结果基本一致，进一步验证了本发明方法的有效性。

表1