CN110940984A - 基于变分分析的双偏振雷达比差分相移快速估算方法 - Google Patents

Abstract

本发明公开了基于变分分析的双偏振雷达比差分相移快速估算方法，步骤为接收双偏振雷达一个平面扫描内，经过去系统相位、去折叠处理后的Φ_DP；将比差分相移K_DP作为变分分析的状态向量，将Φ_DP为K_DP的距离积分作为观测算子，将Φ_DP作为观测，构建变分分析的代价函数；在变分分析的状态向量和变分分析的代价函数中引入B样条插值算子；利用LBFGSB算法最小化代价函数，并在其中引入状态变量为非负的条件，求得非负的K_DP。该方法兼顾了业务天气雷达监测所需要的计算效率以及将K_DP应用到定量降雨估测时所需要的精度。

Description

基于变分分析的双偏振雷达比差分相移快速估算方法

技术领域

本发明涉及一种双偏振雷达比差分相移估算方法，尤其涉及一种基于变分分析的双偏振雷达比差分相移估算方法。

背景技术

双偏振天气雷达可以测量的反射率因子Z_H、差分反射率Z_DR、差分相位Φ_DP、以及零滞后相关系数ρ_hv等双偏振参量，相比常规天气雷达而言，可获取更多降水信息。其中，Φ_DP是相位特征，不受电磁波能量衰减的影响；其衍生产品比差分相移K_DP和降雨率之间的关系接近线性。Φ_DP和K_DP常被用到衰减订正和降水反演中去。K_DP不能被双偏振雷达所观测，需要从Φ_DP中计算。除了K_DP的距离积分的贡献以外，Φ_DP还包含了两部分，即随机误差ε和后向散射相位δ_hv的贡献。在实际计算K_DP的过程中，随机误差和后向散射相位会造成从Φ_DP中计算K_DP时的误差，甚至会导致负值的K_DP，造成误差最终降水估测的误差。因此构建K_DP的优化计算方法能够帮助提高双偏振雷达探测的准确度。

早期计算K_DP的方法是从Φ_DP的定义

出发的,如Istok[1]等使用的分段最小二乘拟合方法即对于Z_H大于(小于等于)40dBZ的观测库，选取该观测前后各1(3)千米的雷达库上的Φ_DP观测，进行最小二乘拟合得到该距离库的K_DP。然而，由于随机误差和后向散射相位的干扰，该方法计算的K_DP可能存在负值，误差较大，在后续应用中精度不高。而Giangrande等[2]提出基于最优化理论的线性规划方法，求取满足K_DP非负条件下的，和观测Φ_DP最接近的传播差分相位φ_DP解和对应的K_DP。这种方法保证了K_DP的非负性，但该方法一方面估计K_DP时受限于定义，只考虑了利用观测Φ_DP，在小雨的情况下，由于随机误差的作用精度较差。另外，该线性规划方法的求解过程较慢，很难应用到实时业务中去。

发明内容

本发明所要解决的技术问题是提出一种基于变分分析的快速估算比差分相移的方法。该方法避免了最小二乘拟合方法中随机误差和后向散射相位对比差分相移估算精度的影响；同时避免了线性规划方法中最终需求的计算量过大的问题。

为了解决上述技术问题，本发明基于变分分析的双偏振雷达比差分相移快速估算方法，包括如下步骤。

步骤1，接收双偏振雷达一个平面扫描内，经过去系统相位、去折叠处理后的Φ_DP。

步骤2，将比差分相移K_DP作为变分分析的状态向量，将Φ_DP为K_DP的距离积分作为观测算子，将Φ_DP作为观测，构建变分分析的代价函数。

步骤3，在变分分析的状态向量和变分分析的代价函数中引入B样条插值算子。

步骤4，利用LBFGSB算法最小化代价函数，并在其中引入状态变量为非负的条件，求得非负的K_DP。

上述技术方案中，步骤1中还需接收经过定标的Z_H、Z_DR、ρ_hv；并利用双偏振雷达参数的自洽关系，从Z_H和Z_DR观测中估计K_DP的背景值以及它的误差分布，并构建变分分析的代价函数中的背景项。

上述技术方案中，对Z_H和Z_DR进行径向滤波，所述双偏振雷达参数的自洽关系利用滴谱仪观测进行拟合。

上述技术方案中，K_DP的背景值为

其中a为10^-1至10^-10，b为0至10，c为0至-10，Z_h和Z_dr分别是Z_H和Z_DR的线性形式。

上述技术方案中，所述K_DP背景值的误差标准差理论值σ_t为

其中σ_DSD为雨滴谱不确定性导致的误差标准差，σ_m为观测误差导致的误差标准差。

上述技术方案中，根据K_DP背景值的误差标准差理论值σ_t、ρ_hv以及Z_H，得到变分分析代价函数中使用的的K_DP的背景值的误差标准差σ′_t。当Z_H大于50dBZ时，σ′_t取值为10的5次方以上任意一个数；当Z_H小于等于50dBZ时，如果σ_t小于0.05且ρ_hv大于等于0.9，σ′_t取值为0.05，如果σ_t小于0.05且ρ_hv小于0.9时，σ′_t取值为0.2；当Z_H小于等于50dBZ时，如果σ_t大于等于0.05且ρ_hv小于0.9，σ′_t取值为4倍的σ_t，如果σ_t大于等于0.05且ρ_hv大于等于0.9，σ′_t取值和σ_t相等。

上述技术方案中，所述变分分析的代价函数为

其中x_b是K_DP的背景值，B是K_DP背景值的误差，W是B样条插值算子，

代表需要估算的K_DP；y是Φ_DP观测，

是Φ_DP观测误差,

是从K_DP预测Φ_DP的观测算子。

上述技术方案中，所述LBFGSB算法最小化代价函数的梯度公式为

其中x_b是K_DP的背景值，B是K_DP背景值的误差，W是B样条插值算子，

代表需要估算的K_DP；y是Φ_DP观测，

是Φ_DP观测误差,

是从K_DP预测Φ_DP的观测算子。

上述技术方案中，B样条插值算子为3阶B样条插值算子。

本发明的基于变分分析的双偏振雷达比差分相移快速估算方法，兼顾了业务天气雷达监测所需要的计算效率以及将K_DP应用到定量降雨估测时所需要的精度。在计算精度方面，K_DP变分分析求解过程中引入了状态变量非负的条件，避免了常规的最小二乘拟合过程中经常存在的误差大，精度低缺陷，另一方面，K_DP变分分析的代价函数引入了利用Z_H和Z_DR观测估计K_DP值以及它的误差，降低了降雨较弱的情况下K_DP估算的精度。在计算效率方面，业务天气雷达扫描为6分钟一次，考虑到一般业务台站的机器性能，K_DP估算所需要的计算资源不能过大；发明中使用了线性算子作为K_DP变分分析的观测算子，具有十分高的计算效率。

附图说明

图1世界时2016年5月9日18:00广州S波段双偏振雷达1.5度仰角扫描，(a-d)分别为K_DP估算所用到的双偏振雷达观测的Z_H、Z_DR、ρ_hv和Φ_DP，接收到的这些参数已经经过偏差订正、非气象回波去除预处理、以及Φ_DP的去折叠。

图2广东地区2016-2017年二维视频滴谱仪观测模拟结果拟合的双偏振量自洽关系(a)，以及观测误差、雨滴谱不确定性统计得到双偏振量自洽性计算的K_DP的误差标准差(b)，模拟的雷达的频率为2.89GHz，温度为20℃。

图3基于图1(a-c)并利用双偏振量自洽性估计的K_DP(a)及其使用在变分分析中的误差标准差σ′_t(b)。

图4基于世界时2016年5月9日18:00广州S波段双偏振雷达1.5度仰角扫描(见图1a-d)，利用发明的变分法估算的K_DP(a)和Φ_DP(b)。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明。图1～4是以广州S波段双偏振多普勒天气雷达为例，展示本发明算法的实现过程和算法性能。

具体步骤如下:

步骤1，接收双偏振雷达一个平面扫描内，经过定标的Z_H、Z_DR和ρ_hv，经过去系统相位、去折叠处理后的Φ_DP，如图1所示。

步骤2，对获取的Z_H和Z_DR进行径向7点滤波；同时，利用双偏振雷达参数之间的自洽性，从Z_H和Z_DR中计算K_DP背景值，双偏振雷达参数之间的自洽关系可以利用滴谱仪观测进行拟合。如图2中的a，利用广东地区二维视频滴谱仪在2016-2017年观测拟合的关系为例，K_DP的背景值

其中Z_h和Z_dr分别是Z_H和Z_DR的线性形式

图1中的结果计算的K_DP背景值如图3中a。

步骤3，计算K_DP背景值

的误差的理论值。K_DP背景值的误差标准差的理论值为σ_t，主要包括两部分：(1)由分段统计散点的标准差，参见图2a；雨滴谱不确定性导致的误差标准差σ_DSD，参见图2b的划线；(2)由观测误差导致的误差标准差σ_m，它和Z_H、Z_DR误差之间的关系为

参见图2b中的点线；最终

的误差标准差σ_t为

参见图2b的实线。

步骤4，根据K_DP背景值的误差标准差理论值σ_t、ρ_hv以及Z_H，得到变分分析代价函数中使用的的K_DP的背景值的误差标准差σ′_t：(1)当Z_H大于50dBZ时，σ′_t取值为10的5次方以上任意一个数，本发明中使用了10的20次方。结果如图3b所示，其中当Z_H大于50dBZ时的σ_t为空白，也就未填色。(2)当Z_H小于等于50dBZ时，如果σ_t小于0.05且ρ_hv大于等于0.9，σ′_t取值为0.05，如果σ_t小于0.05且ρ_hv小于0.9时，σ′_t取值为0.2；(3)当Z_H小于等于50dBZ时，如果σ_t大于等于0.05且ρ_hv小于0.9，σ′_t取值为4倍的σ_t，如果σ_t大于等于0.05且ρ_hvd阿玉等于0.9，σ′_t取值和σ_t相等。

步骤5，对步骤1中获取的雷达数据逐个径向数据进行处理，假设一个雷达径向有n个有降水的距离库，这些距离库上的K_DP组成向量x＝[K_DP,1,K_DP,2,…,K_DP,n]^T，将x表示成矩阵W和长度为m的向量

的乘积

其中(m<n)。W中的元素是由3次B样条插值算子组成(Prenter 1975)；在此基础上构建变分估算K_DP的代价函数

其中x_b是该雷达径向上获得的

B是对角矩阵，对角元素分别是

的误差方差σ′_t×σ′_t；对角矩阵

的对角元素对应该径向上的各距离库Φ_DP的观测误差方差，假定Φ_DP的观测的标准差为2度为例；n×n的雅可比矩阵H中第i行和第j列元素为：当i≤j时等于0，当i>j时等于2δr，其中δr为雷达的距离库长。

步骤6，利用LBFGS-B算法(limited-memory Broyden-Fletcher-Goldfarb-Shannoalgorithm for bound optimization)计算步骤5中的代价函数最小时的

同时在最小化时设置

的元素大于0；在最小化时用到的代价函数的梯度的公式为

步骤7，利用

从步骤6)中获得的

计算每个步骤5中选取的径向的各距离库上的K_DP。

步骤8，选取下一个雷达径向，重复步骤5-步骤8直至该雷达扫描结束。

在变分分析的状态向量以及代价函数中引入B样条插值算子，这一步的目的是在估计K_DP时保证它的空间连续性，降低随机误差的影响。

在变分分析中使用了固定的雅可比矩阵H，不随着代价函数的最小化过程而变化。

Claims (9)

1.基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：

步骤1，接收双偏振雷达一个平面扫描内，经过去系统相位、去折叠处理后的Φ_DP；

步骤2，将比差分相移K_DP作为变分分析的状态向量，将Φ_DP为K_DP的距离积分作为观测算子，将Φ_DP作为观测，构建变分分析的代价函数；

步骤3，在变分分析的状态向量和变分分析的代价函数中引入B样条插值算子；

步骤4，利用LBFGSB算法最小化代价函数，并在其中引入状态变量为非负的条件，求得非负的K_DP。

2.如权利要求1所述的基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：还需接收经过定标的Z_H、Z_DR、ρ_hv；并利用双偏振雷达参数的自洽关系，从Z_H和Z_DR观测中估计K_DP的背景值以及它的误差分布，并构建变分分析的代价函数中的背景项。

3.如权利要求2所述的基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：对Z_H和Z_DR进行径向滤波，所述双偏振雷达参数的自洽关系利用滴谱仪观测进行拟合。

4.如权利要求3所述的基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：K_DP的背景值为

其中a为10^-1至10^-10，b为0至10，c为0至-10，Z_h和Z_dr分别是Z_H和Z_DR的线性形式。

5.如权利要求4所述的基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：获取K_DP背景值的误差标准差理论值σ_t为

其中σ_DSD为雨滴谱不确定性导致的误差标准差，σ_m为观测误差导致的误差标准差。

6.如权利要求5所述的基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：根据K_DP背景值的误差标准差理论值σ_t、ρ_hv以及Z_H，得到变分分析代价函数中使用的的K_DP的背景值的误差标准差σ′_t。当Z_H大于50dBZ时，σ′_t取值为10的5次方以上任意一个数；当Z_H小于等于50dBZ时，如果σ_t小于0.05且ρ_hv大于等于0.9，σ′_t取值为0.05，如果σ_t小于0.05且ρ_hv小于0.9时，σ′_t取值为0.2；当Z_H小于等于50dBZ时，如果σ_t大于等于0.05且ρ_hv小于0.9，σ′_t取值为4倍的σ_t，如果σ_t大于等于0.05且ρ_hv大于等于0.9，σ′_t取值和σ_t相等。

7.如权利要求6所述的基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：所述变分分析的代价函数为

其中x_b是K_DP的背景值，B是K_DP背景值的误差，W是B样条插值算子，

代表需要估算的K_DP；y是Φ_DP观测，

是Φ_DP观测误差，

是从K_DP预测Φ_DP的观测算子。

8.如权利要求7所述的基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：所述LBFGSB算法最小化代价函数的梯度公式为

其中x_b是K_DP的背景值，B是K_DP背景值的误差，W是B样条插值算子，

代表需要估算的K_DP；y是Φ_DP观测，

是Φ_DP观测误差，

是从K_DP预测Φ_DP的观测算子。

9.如权利要求8所述的基于变分分析的双偏振雷达比差分相移快速估算方法，其特征在于：B样条插值算子为3阶B样条插值算子。

Images