CN106872981A - 雨量雷达的降水强中心跟踪与预报方法 - Google Patents

Abstract

本发明提出了一种雨量雷达的降水强中心跟踪与预报方法，包括：跟踪并获得TREC矢量场的步骤，对移动矢量场进行平滑处理的步骤，降水强中心回波位置预报和回波强度预报的步骤。本发明利用TREC方法对X波段中频相参雨量雷达进行降水中心跟踪，消除了因地物杂波存在产生的错误的移动方向或零矢量，提高了TREC矢量场的连续性，此外针对多个数据移入的格点和无数据格点的特殊情况进行了额外处理，使得预报准确性大大提高。

Description

雨量雷达的降水强中心跟踪与预报方法

技术领域

本发明属于雨量雷达数据处理技术领域，具体涉及雨量雷达的降水强中心相关数据分析方法。

背景技术

对雷达回波的跟踪和外推是进行自动化强对流天气临近预报的关键。新一代多普勒天气雷达的降水强度信息有较高的时空分辨率，是监测暴雨、冰雹，进行临近预报的重要手段之一。

我们可以利用雨量雷达对暴雨进行实时监测，分析暴雨的一些重要特征，对其进行预警和临近预报，以减少灾害的发生。以往的监测方法有质心跟踪法，是对某块回波相继两时刻的距心位置进行最小二乘拟合后外推预报的。这种方法在算法上是把每块回波简化为一点，所以比较适合对孤立的大而强的回波块进行跟踪，而对于零散的回波，以及回波发生合并或分列时，跟踪成功率较低。后来广泛发展的是TREC方法，TREC方法是交叉相关法的一种发展，这种方法较之质心跟踪法预报准确性上有所提高，但仍需改善。因此，目前的几种方法均存在缺点。

发明内容

为解决上述问题，本发明基于TREC算法进行改进，提出了一种预报准确度更为精确的雨量雷达的降水强中心跟踪与预报方法。

为了达到上述目的，本发明提供如下技术方案：

雨量雷达的降水强中心跟踪与预报方法，包括如下步骤：

步骤一，跟踪并获得TREC矢量场

基于以下方程获得强度Z：

其中，P_t为发射机功率，G为天线增益，θ和为水平和垂直波束宽度，λ为波长，R为距离，Z为雷达反射率，为目标物的散射特征，ψ为波束冲塞系数，k为衰减系数；

将回波区域划分成若干个矩形区域，跟踪每个区域的移动，将第一时刻的回波图像中的一个矩形区域，在搜索半径内向任一方向移过一定的距离，通过下式计算此矩形区域与第二时刻相同大小的矩形区域之间的交叉相关系数R：

式中Z₁，Z₂分别为T时刻和T+Δt时刻反射率因子(雨强)的矩阵，n为矩阵的数据点数，求出间隔Δt时间的两个矩阵的相关系数，重复该过程，直到找到最大的相关系数，找到最大相关系数时，得到该T时刻矩形区域的中心位置指向T+Δt时刻矩形区域的中心位置的矢量为TREC矢量；

步骤二，对移动矢量场进行平滑处理

步骤二-1，当一个TREC矢量是零矢量或与其周围矢量的平均方向相差大于或等于20度，则用其周围8个点的平均矢量来替代这个TREC矢量；

步骤二-2，对TREC矢量的u，v分量分别进行客观分析，采用如下公式得到连续的位移矢量场：

α^*(i,j)＝α₀(i,j)+Δα(i,j)

其中，α₀(i,j)为变量α在格点(i，j)处的估测值；α^*(i,j)为变量α在格点(i，j)处的校正值；Δα(k)为目标点k处的差值，这里为该点TREC矢量的分量与估测值之差；w(i,j,k)为s目标点k在格点(i，j)处的加权函数；d_m为分析点与格点之间的距离及TREC矢量的函数，R是影响半径；

步骤三：降水强中心预报

步骤三-1，回波位置预报

在TREC矢量的u、v分量上分别应用半拉格朗日一维平流公式，分别计算出u、v分量上的位移，每间隔一定时间步长输出一次预报结果，该时刻每个格点上的雨强值及移向、移速，用作下一连续时刻的积分，所述半拉格朗日一维平流公式为：

其中，U(x,t)为给定的函数，在这里U(x,t)为TREC矢量u、v分量中一个，x方向与TREC矢量u、v分量的方向相同，公式(7)说明标量Z沿着轨迹是不变的，在这里Z为反射率因子值；

根据公式(7)，有

x_*＝x_j-U(x,t)Δt

其中，为n+1时刻在点j处的反射率因子值；为n时刻在点*处的反射率因子值；x_*为点*的x坐标；x_j为点j的x坐标；U(x,t)为TREC矢量的一个分量；Δt为n+1时刻与n时刻的时间差；

步骤三-2，回波强度预报

采用线性外推的方法，通过连续两个时刻的雷达反射率因子进行线性拟合，得到线性外推的回波强度。

进一步的，所述步骤一中搜索半径根据最大期望速度计算：

S_半径＝V_max×Δt

式中V_max——雨带移动的最大期望速度(这里采用60公里/小时)，Δt为时间间隔。

进一步的，步骤三-2中，对于有多个数据移入的格点，本发明采用加权平均的方法，具体公式如下：

上式中，w_i通过如下公式进行计算：

其中w_i为第i个移入格点回波强度Z_i对应的权重值，N为移入同一格点的个数；Zmin、Zmax分别为移入格点数据中回波强度最小值、最大值。

进一步的，当无数据移入的格点周围超过一半的数值有效时，进行3×3均值替代该无格点。

与现有技术相比，本发明具有如下优点和有益效果：

本发明提供的雨量雷达的降水强中心跟踪与预报方法，利用TREC方法对X波段中频相参雨量雷达进行降水中心跟踪，消除了因地物杂波存在产生的错误的移动方向或零矢量，提高了TREC矢量场的连续性，此外针对多个数据移入的格点和无数据格点的特殊情况进行了额外处理，使得预报准确性大大提高。

附图说明

图1为TREC方法示意图。

图2为两个时间层的半拉格朗日积分方法。

图3为个例1回波风场结构外推和雨强外推图。

图4为个例2回波风场结构外推和雨强外推图。

图5为个例3回波风场结构外推和雨强外推图。

具体实施方式

以下将结合具体实施例对本发明提供的技术方案进行详细说明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

雨量雷达系统发射的射频微波经过天线馈线分系统的馈线部分到达天线，向空间定向辐射。天线定向辐射的电磁波能量遇到云、雨等降水目标时，便会产生散射，其中后向散射的一部分，形成气象目标的射频回波信号被天线接收。天线接收到的射频回波信号，经过馈线部分被送往接收分系统，由数字中频转换器进行A/D变换后成为数字中频回波信号送往信号处理分系统。信号处理分系统对来自接收分系统的数字中频回波信号，做DSU处理后在数字域内处理形成I、Q正交信号，并对其作平均处理，地物对消滤波处理，得到反射率的估测值即强度Z；通过脉冲对处理(PPP)或快速富里叶变换(FFT)处理，从而得到散射粒子群的平均径向速度V和速度的平均起伏即速度谱宽W。

具体地说，本发明方法包括如下步骤：

步骤一，跟踪并获得TREC矢量场

强度Z通过Z-R关系得到，具体采用以下方程获得：

其中，P_t为发射机功率，G为天线增益，θ和为水平和垂直波束宽度，λ为波长，R为距离，Z为雷达反射率，为目标物的散射特征，ψ为波束冲塞系数，k为衰减系数。

将回波区域划分成若干个矩形区域，如图1所示，跟踪每个区域的移动，将第一时刻的回波图像中的一个矩形区域，在搜索半径内向任一方向移过一定的距离。

上述搜索半径根据最大期望速度计算：

S_半径＝V_max×Δt (2)

式中V_max——雨带移动的最大期望速度(这里采用60公里/小时)，Δt为时间间隔。

然后通过下式计算此矩形区域与第二时刻相同大小的矩形区域之间的交叉相关系数R：

式中Z₁，Z₂分别为T时刻和T+Δt时刻反射率因子(雨强)的矩阵。n为矩阵的数据点数。那么，通过上式就可以求出间隔Δt时间的两个矩阵的相关系数，重复这个过程，直到找到最大的相关系数。

对于不同的移动位置，会得到不同的相关系数值，直到找到极大值R_max为止，具有最大相关系数的移动就是TREC矢量，即用作外推预报的位移矢量。通过最大相关系数法找到当前回波位置和下时刻回波，即可确定该回波移动的方向和距离，得到TREC移动矢量。通过公式(3)找到最大相关系数时，从该T时刻矩形区域的中心位置指向T+Δt时刻矩形区域的中心位置的矢量即为TREC矢量。据此外推T+N*Δt时刻暴雨的位置。

步骤二，对移动矢量场进行平滑处理

由TREC方法得到的移动矢量场，由于地物杂波等因素影响，会存在一些错误的移动方向或零矢量，为了消除这些错误的移动方向，提高TREC矢量场的连续性，我们采取了以下两个步骤进行连续性检查：

步骤二-1，如果一个TREC矢量与其周围矢量的平均方向相差超过(大于或等于)20度，则用其周围8个点的平均矢量来替代这个TREC矢量，零矢量也采用这种办法。

步骤二-2，对TREC矢量进行客观分析，得到连续的位移矢量场：

对TREC矢量的u，v分量分别进行客观分析，客观分析采用如下公式：

α^*(i,j)＝α₀(i,j)+Δα(i,j) (4)

其中，α₀(i,j)为变量α在格点(i，j)处的估测值(格点是雷达测量雨量等到的最小范围，是雷达扫描投影到地面的一个正方型区域。不同型号雷达差异较大，X波段雷达最小可以到60米*60米。)，这里取TREC矢量场分量的平均值即用周围8个点的u、v矢量平均值；α^*(i,j)为变量α在格点(i，j)处的校正值；Δα(k)为目标点k处的差值，这里为该点TREC矢量的分量与估测值之差；w(i,j,k)为s目标点k在格点(i，j)处的加权函数；d_m为分析点与格点之间的距离及TREC矢量的函数，R是影响半径。

步骤三：降水强中心预报。预报是以运动估算、反射率因子的增长和减弱信息以及当前数据为基础，进行外推。对于每一个矩形区域来说，通过计算其均值的变化可以得到云团反射率因子增长减弱的信息，以此预报反射率因子的趋势变化。通过TREC矢量得到回波的移动方向和速度，预报回波的位置。也可以拟合连续5个时刻的TREC矢量场，作为回波的移动方向和速度。本步骤主要包括两个分步骤——回波位置预报和回波强度预报。

步骤三-1，回波位置预报

在对台风暴雨预报的时候，由于台风雨带作旋转移动，如果预报时效较长，线性外推必然引起较大的误差，因此本发明采用两个时间层的半拉格朗日积分方法(图2所示)，半拉格朗日一维平流公式为：

U(x,t)为给定的函数，在这里U(x,t)为TREC矢量u、v分量中一个，x方向与TREC矢量u、v分量的方向相同，公式(7)说明标量Z沿着轨迹是不变的，在这里Z为反射率因子值。根据公式(7)，在图2中有

x_*＝x_j-U(x,t)Δt (9)

其中，为n+1时刻在点j处的反射率因子值；为n时刻在点*处的反射率因子值；x_*为点*的x坐标；x_j为点j的x坐标；U(x,t)为TREC矢量的一个分量；Δt为n+1时刻与n时刻的时间差。

在对台风暴雨的外推预报中，在TREC矢量的u、v分量上分别应用半拉格朗日一维平流公式，时间步长为5分钟，分别计算出u、v分量上的位移，每5分钟输出一次预报结果，该时刻每个格点上的雨强值及移向、移速，用作下一连续时刻的积分。

步骤三-2，回波强度预报

回波强度的变化是一个复杂的过程，连续两个时刻的雷达反射率因子的差异可以作为预报回波强度变化的信息。在预报回波强度时，我们采用线性外推的方法，即通过连续两个时刻的雷达反射率因子进行线性拟合，得到线性外推的回波强度。由于TREC技术得到回波的多个移动矢量在外推的过程中，难免有些格点的位置会有多个数据移入，同时又有些格点会没有数据移入，对于有多个数据移入的格点，本发明采用加权平均的方法，具体公式如下：

上式中，w_i通过如下公式进行计算：

其中w_i为第i个移入格点回波强度Z_i对应的权重值，N为移入同一格点的个数；Zmin、Zmax分别为移入格点数据中回波强度最小值、最大值；权重的大小与数据值的大小有关，数值大的权重大，数值小的权重小。为了消除那些无数据格点造成的影响，如果当无数据移入的格点周围超过一半的数值有效，即进行3×3均值替代该无格点，以达到整体数据平滑处理效果。

采用本发明方法进行强中心跟踪与外推的个例分析数据如下：

个例1

江西赣州雨量雷达在2016年03月20日02：08时刻观测到一次混合降水天气过程，采用地形匹配模式进行扫描。如图3为该时刻回波风场结构外推和雨强外推反演。图3a为02:08时刻雨量雷达实时采集的回波强度；图3b为下一时刻回波强度和风场结构外推效果图；图3c为02:18时刻雨量雷达反演的雨量强度图；图3d为02:18时刻雨强外推效果。由此可知，强中心回波外推效果较好(椭圆区)，和实况基本吻合，降水边缘出现虚假弱回波，分析原因是因为降水的发展是生消过程，且是非线性过程，上一时刻的弱降水发展趋势较强降水回波更难以控制。

个例2

江西省南昌市雨量雷达在2015年06月01日18:38时刻观测到一次强对流天气过程，采用地形匹配模式进行扫描。如图4为该时刻回波风场结构外推和雨强外推反演。图4a为18:38时刻雨量雷达实时采集的回波强度；图4b为18:43时刻回波强度和风场结构外推效果图；图4c为18:43时刻雨量雷达反演的雨量强度图；图4d为18:43时刻雨强外推效果。由此可知，强中心回波外推效果较好(椭圆区)，和实况基本吻合。

个例3

江西省赣州市雨量雷达在2016年08月13日15:28时刻观测到一次局地强对流天气过程，采用地形匹配模式进行扫描。如图5为该时刻回波风场结构外推和雨强外推反演。图5a为15:28时刻雨量雷达实时采集的回波强度；图5b为15:33时刻回波强度和风场结构外推效果图；图5c为15:38时刻雨量雷达反演的雨量强度图；图5d为15:38时刻雨强外推效果。由此可知，强中心回波外推效果较好(椭圆区)，和实况基本吻合。

本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (4)

1.雨量雷达的降水强中心跟踪与预报方法，其特征在于，包括如下步骤：

步骤一，跟踪并获得TREC矢量场

基于以下方程获得强度Z：

其中，P_t为发射机功率，G为天线增益，θ和为水平和垂直波束宽度，λ为波长，R为距离，Z为雷达反射率，为目标物的散射特征，ψ为波束冲塞系数，k为衰减系数；

将回波区域划分成若干个矩形区域，跟踪每个区域的移动，将第一时刻的回波图像中的一个矩形区域，在搜索半径内向任一方向移过一定的距离，通过下式计算此矩形区域与第二时刻相同大小的矩形区域之间的交叉相关系数R：

$R=\frac{{\mathrm{\ΣZ}}_1\left(i\right)Z_2\left(i\right)-n^{-1}{\mathrm{\ΣZ}}_1\left(i\right){\mathrm{\ΣZ}}_2\left(i\right)}{{\[\left({\mathrm{\ΣZ}}_1^2\right(i)-n{\stackrel{\‾}{Z}}_1^2)\left({\mathrm{\ΣZ}}_2^2\right(i)-n{\stackrel{\‾}{Z}}_2^2)\]}^{1/2}}$

其中，Z₁，Z₂分别为T时刻和T+Δt时刻反射率因子(雨强)的矩阵，n为矩阵的数据点数，求出间隔Δt时间的两个矩阵的相关系数，重复该过程，直到找到最大的相关系数，找到最大相关系数时，得到该T时刻矩形区域的中心位置指向T+Δt时刻矩形区域的中心位置的矢量为TREC矢量；

步骤二，对移动矢量场进行平滑处理

步骤二-1，当一个TREC矢量是零矢量或与其周围矢量的平均方向相差大于或等于20度，则用其周围8个点的平均矢量来替代这个TREC矢量；

步骤二-2，对TREC矢量的u，v分量分别进行客观分析，采用如下公式得到连续的位移矢量场：

α^*(i,j)＝α₀(i,j)+Δα(i,j)

$\mathrm{\Δ}\mathrm{\α}(i,j)=\frac{\underset{k}{\Σ}w{(i,j,k)}^2\×\mathrm{\Δ}\mathrm{\α}\left(k\right)}{\underset{k}{\Σ}w(i,j,k)}$

$w(i,j,k)=\left\{\begin{array}{cc}\frac{R^2-d_m^2}{R^2+d_m^2},& d_m^2<R^2\\ 0,& d_m^2>R^2\end{array}\right.$

其中，α₀(i,j)为变量α在格点(i，j)处的估测值；α^*(i,j)为变量α在格点(i，j)处的校正值；Δα(k)为目标点k处的差值，这里为该点TREC矢量的分量与估测值之差；w(i,j,k)为s目标点k在格点(i，j)处的加权函数；d_m为分析点与格点之间的距离及TREC矢量的函数，R是影响半径；

步骤三：降水强中心预报

步骤三-1，回波位置预报

在TREC矢量的u、v分量上分别应用半拉格朗日一维平流公式，分别计算出u、v分量上的位移，每间隔一定时间步长输出一次预报结果，该时刻每个格点上的雨强值及移向、移速，用作下一连续时刻的积分，所述半拉格朗日一维平流公式为：

$\frac{dZ}{dt}=\frac{\&part;Z}{\&part;t}+U(x,t)\frac{\&part;Z}{\&part;x}=0$

其中，U(x,t)为给定的函数，在这里U(x,t)为TREC矢量u、v分量中一个，x方向与TREC矢量u、v分量的方向相同，公式(7)说明标量Z沿着轨迹是不变的，在这里Z为反射率因子值；

根据公式(7)，有

$Z_j^{n+1}=Z_{*}^n$

x_*＝x_j-U(x,t)Δt

其中，为n+1时刻在点j处的反射率因子值；为n时刻在点*处的反射率因子值；x_*为点*的x坐标；x_j为点j的x坐标；U(x,t)为TREC矢量的一个分量；Δt为n+1时刻与n时刻的时间差；

步骤三-2，回波强度预报

采用线性外推的方法，通过连续两个时刻的雷达反射率因子进行线性拟合，得到线性外推的回波强度。

2.根据权利要求1所述的雨量雷达的降水强中心跟踪与预报方法，其特征在于，所述步骤一中搜索半径根据最大期望速度计算：

S_半径＝V_max×Δt

式中V_max——雨带移动的最大期望速度(这里采用60公里/小时)，Δt为时间间隔。

3.根据权利要求1所述的雨量雷达的降水强中心跟踪与预报方法，其特征在于：对于有多个数据移入的格点，采用下式进行加权平均：

$Z=\frac{\underset{i=1}{\overset{N}{\&Sigma;}}{w}_i\&times;Z_i}{N}$

上式中，w_i通过如下公式进行计算：

$w_i=\frac{Z_i-Z_{min}}{Z_{\max }-Z_{min}}$

其中，w_i为第i个移入格点回波强度Z_i对应的权重值，N为移入同一格点的个数；Zmin、Zmax分别为移入格点数据中回波强度最小值、最大值。

4.根据权利要求1所述的雨量雷达的降水强中心跟踪与预报方法，其特征在于：当无数据移入的格点周围超过一半的数值有效时，进行3×3均值替代该无格点。