CN105974418A - 基于天气雷达反射率特征匹配的降水估测方法 - Google Patents

Abstract

本发明公开了基于天气雷达反射率特征匹配的降水估测方法，首先构建雷达反射率特征库，然后确立实况雷达数据与雷达反射率特征库的匹配规则，并在此基础上动态确立Z‑R关系中最优的参数A和参数b。本发明在短临降水预报业务应用中实现Z‑R关系的动态最优，进一步提高雷达估测降水的准确性。

Description

基于天气雷达反射率特征匹配的降水估测方法

技术领域

本发明属于大气科学领域，特别涉及了基于天气雷达反射率特征匹配的降水估测方法。

背景技术

天气雷达(以下简称雷达)是强对流天气探测的重要手段，它能够探测其扫描半径内的雨强和雨量分布情况，进一步可以估测总降水量和降水率，尤其在短时临近预报等业务中已得到广泛应用。

基于Z-R关系的计算是雷达估测降水量的关键方法，而估测降水量的准确率很大程度上依赖于Z-R关系(Z＝AR^b)中参数A和参数b的选择。很多学者通过对地面实测降水量与雷达回波的分析研究，得出特定区域和(或)特定时期内相对最优的参数A和参数b，由此作为未来对同区域、同气象条件下雷达估测降水Z-R关系计算时的重要参考。

在短时临近预报业务软件中，其参数A和参数b或是针对当地气象条件和降水类型设定为一组相对最优值，或是按时期、按降水类型(如大陆强对流降水型、热带降水型等)设定多组相对最优值。

然而，由于雷达回波的时空变化性很大，从雷达回波图像上看，几乎每次降水过程的回波图像都是独一无二的，其Z-R关系的最优参数A和参数b也并非固定不变的。尽管部分气象业务软件中Z-R关系的参数是可调的，但在某一降水事件发生前，无论人工或是软件自动地确定适当的参数都是非常困难的。因此，如何在短临降水预报业务应用中，实现Z-R关系的动态最优，进一步提高雷达估测降水的准确性需要特别关注。

发明内容

为了解决上述背景技术提出的技术问题，本发明旨在提供基于天气雷达反射率特征匹配的降水估测方法，通过构建雷达反射率特征库，以及将反射率特征与雷达反射率特征库进行匹配，动态确立Z-R关系中最优的参数A和参数b，提高了降水估测的准确性。

为了实现上述技术目的，本发明的技术方案为：

基于天气雷达反射率特征匹配的降水估测方法，包括以下步骤：

步骤一、构建天气雷达反射率特征库：

(1)提取天气雷达基数据中的基本反射率数据，并将基本反射率数据从原来的极坐标形式转换为平面直角坐标系形式；分别计算p千米和q千米高度处的等高平面位置显示图CAPPI_p和CAPPI_q，CAPPI_p和CAPPI_q的结构相同，均为N×N的二维网格，求取CAPPI_p和CAPPI_q各网格点在同一水平位置的最大回波强度值CAPPI_MAX，即CAPPI_MAX(x,y)＝Max(CAPPI_p(x,y),CAPPI_q(x,y))，其中x、y分别表示平面直角坐标系下的横坐标和纵坐标，Max为取最大值函数；

(2)根据得到的最大回波强度值CAPPI_MAX(x,y)，计算雷达回波的灰度图像Rad(x,y)：

$Rad(x,y)=\left\{\begin{array}{cc}0,& CAPPI\_MAX(x,y)<20\\ CAPPI\_MAX(x,y),& 20\&le;CAPPI\_MAX(x,y)\&le;64\\ 64,& CAPPI\_MAX(x,y)>64\end{array}\right.;$

(3)计算灰度图像Rad(x,y)的一阶梯度值，从而强化回波图像的变化特征；

(4)新建一个M×M的雷达回波特征矩阵Fea，M<N/2，Fea中每个单元的值计算方式为：

Fea(x',y')＝Hex(Fea'(x',y'))

${\mathrm{Fea}}^{,}(x^{,},y^{,})=\frac{\underset{k_0=0}{\overset{STEP-1}{\&Sigma;}}\underset{k_1=0}{\overset{STEP-1}{\&Sigma;}}Deg(x^{,}\&times;STEP+k_0,y^{,}\&times;STEP+k_1)}{STEP\&times;STEP\&times;4}$

其中，x'∈[0,M-1]，y'∈[0,M-1]，Hex为对括号中的值作16进制转换的函数；

(5)将得到的雷达回波特征矩阵Fea，从Fea(0,0)到Fea(M-1,M-1)逐个添加到一个字符串变量FeaStr中；

(6)选取雷达有效探测范围内的地面雨量观测设备，以当前雷达探测时间为基准，提取这些地面雨量观测设备在基准后的x小时内的累积观测降水量，形成数据集ORF，其中ORF[c]表示第c个地面雨量观测设备的x小时降水量，x∈[1,24]，将当前雷达有效探测范围内地面雨量观测设备的总数记为MSum，即数据集ORF的大小为MSum，则c∈[0,MSum-1]；

(7)根据数据集ORF中每个地面雨量观测设备的经纬度信息，以及当前雷达的地理坐标信息，提取各个地面雨量观测设备所对应地理位置的CAPPI_MAX值；定义CM为CAPPI_MAX的一个子集，CM中的每个值CM[c]为与ORF[c]地理位置最相邻的网格点的回波强度值；设置两个递增变量PA和Pb，其中PA∈[100,400]，以10为步长递增，Pb∈[1.0,2.0]，以0.1为步长递增，分别计算ORF与CM中各记录项的误差平方和：

$D(PA,Pb)=\underset{c=1}{\overset{MSum}{\&Sigma;}}{(ORF\&lsqb;c\&rsqb;-\sqrt[Pb]{\frac{{10}^{\frac{CM\&lsqb;c\&rsqb;}{10}}}{PA}})}^2;$

(8)记录下所得D(PA,Pb)取值最小时的一组PA和Pb，分别记为Para_A和Para_b；

(9)将步骤(5)中得到的FeaStr和步骤(8)得到的Para_A和Para_b添加到数据库；

(10)依据上述步骤(1)-(9)，将有资料以来，雷达探测数据和地面雨量观测数据逐一进行计算并入库，从而构建起天气雷达反射率特征库；

步骤二：将实况天气雷达反射率数据与建立的天气雷达反射率特征库进行匹配：

(11)将需要匹配的天气雷达反射率数据按照上述步骤(1)-(5)，计算出能够反映该雷达回波特征的字符串FeaStrCur；

(12)从天气雷达反射率特征库中遍历同一部天气雷达的FeaStr，形成数据集FS，FS中的每条记录FS[c]与FeaStrCur作相关性计算，计算方法为：

$R\left(c\right)=\frac{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}\left(FeaStrCur\right(k)-\stackrel{\&OverBar;}{FeaStrxCur})\&CenterDot;(FS\&lsqb;c\&rsqb;(k)-\stackrel{\&OverBar;}{FS\&lsqb;c\&rsqb;})}{\sqrt{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}{\left(FeaStrCur\right(k)-\stackrel{\&OverBar;}{FeaStrCur})}^2}\&CenterDot;\sqrt{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}{(FS\&lsqb;c\&rsqb;(k)-\stackrel{\&OverBar;}{FS\&lsqb;c\&rsqb;})}^2}}$

其中，和分别表示FeaStrCur与FS[c]中所有记录的算术平均值，FS[c](k)表示FS[c]的第k项记录值；

(13)找到R(c)取值最大时，FS[c]所对应的数据库记录，该记录中对应的Para_A和Para_b即为当前天气雷达资料下Z-R关系的最优参数A和参数b，并据此估测降水。

进一步地，在步骤(1)中，取p＝1.5，q＝3.0。

进一步地，在步骤(3)中，灰度图像Rad(x,y)最边缘一圈的各像素的一阶梯度为0，其余像素的一阶梯度采用Sobel检测算子来计算。

进一步地，横向Sobel检测算子纵向Sobel检测算子

进一步地，在步骤(5)中，字符串变量FeaStr中每2个字节代表1个Fea(x’,y’)。

进一步地，在步骤(13)中，若找到的最大R(c)<0.5时，则Z-R关系中的最有参数A和参数b的取值改用预先设定值。

采用上述技术方案带来的有益效果：

(1)本发明对雷达回波进行了一阶梯度信息提取，强化雷达回波强弱区间的边缘特征，为不同时刻雷达反射率的比较、匹配提供了分析的数据基础。

(2)本发明将雷达回波特征信息以字符串形式加以存储，便于数据库的实现以及特征匹配时的快速计算。

(3)本发明中基于字符串和16进制数据格式的特征值的相关性计算，具有较高的计算性能，有利于本方法业务化实施时，提高雷达估测降水预报的时效性。

附图说明

图1是本发明中构建雷达反射率特征库的流程图，对应步骤1-9，图中实线表示程序流，虚线表示数据流，当程序流与数据流重叠时，虚线省略；

图2是本发明中实况雷达数据与雷达反射率特征库的匹配的流程图，对应步骤11-13，图中FS.Count表示数据集FS中的总记录数。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

本实施例选用2010年至2015年期间，南京地区多普勒天气雷达9层体积扫的基本反射率数据。该雷达数据格式为CINRAD-SA，体扫模式为VCP21，时间分辨率为6分钟，径向分辨率为1km，所处地理位置为东经118.698°、北纬32.191°。剔除没有明显降水回波的雷达数据，并对剩余待计算的雷达数据作必要的质量控制。然后，对这些雷达数据逐一作如下处理：

步骤1、不妨设当前待处理的雷达基数据文件名为RF.bin，提取该文件中的基本反射率数据，将原有极坐标形式的数据结构转换为平面直接坐标形式。分别计算1.5km和3.0km高度的CAPPI，得到CAPPI_15和CAPPI_30。再求取这两个高度，各网格点在同一水平位置上的最大回波强度值，即CAPPI_MAX(x,y)＝Max(CAPPI_p(x,y),CAPPI_q(x,y))，其中x、y分别表示平面直角坐标系下的横坐标和纵坐标，Max表示取最大值的函数。根据该部雷达的性能技术指标，定义N＝230，因此x∈[0,229]，y∈[0,229]。

步骤2、由CAPPI_MAX，输出雷达回波的灰度图像，记为Rad，Rad中的值满足：

$Rad(x,y)=\left\{\begin{array}{cc}0,& CAPPI\_MAX(x,y)<20\\ CAPPI\_MAX(x,y),& 20\&le;CAPPI\_MAX(x,y)\&le;64\\ 64,& CAPPI\_MAX(x,y)>64\end{array}\right..$

步骤3、采用Sobel检测算子计算Rad的一阶梯度值，计算方法为：

$Deg(x,y)=\sqrt{{\left(\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]*Rad\right)}^2+{\left(\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]*Rad\right)}^2}$

其中，*表示卷积计算，展开后的计算方法形如：

$\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]*Rad\&DoubleRightArrow;\begin{array}{c}-Rad(x-1,y-1)+Rad(x+1,y-1)\\ -2ad(x-1,y)+2\&times;Rad(x+1,y)\\ -Rad(x-1,y+1)+Rad(x+1,y+1)\end{array}$

$\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]*Rad\&DoubleRightArrow;\begin{array}{c}-Rad(x-1,y-1)-2\&times;Rad(x,y-1)-Rad(x+1,y-1)\\ Rad(x-1,y+1)+2\&times;Rad(x,y+1)+Rad(x+1,y+1)\end{array}$

步骤4、新建一个M×M的雷达回波特征矩阵Fea，M＝46，SETP＝5，Fea中每个单元的值计算方式为：

Fea(x',y')＝Hex(Fea'(x',y'))

${\mathrm{Fea}}^{,}(x^{,},y^{,})=\frac{\underset{k_0=0}{\overset{STEP-1}{\&Sigma;}}\underset{k_1=0}{\overset{STEP-1}{\&Sigma;}}Deg(x^{,}\&times;STEP+k_0,y^{,}\&times;STEP+k_1)}{STEP\&times;STEP\&times;4}$

其中，SETP＝5，x’∈[0,45]，y’∈[0,45]，Hex表示一个函数，功能是对其括号中的值作16进制转换。

步骤5、将Fea中的各个元素，从Fea[0,0]、Fea[0,1]、Fea[0,2]依次到Fea[45,45]逐个添加到一个字符串变量FeaStr中。例如，假设Fea[0,0]、Fea[0,1]、Fea[0,2]至Fea[0,5]值分别为：0、61、18、42、59、29，则FeaStr中的记录形如：003D122A3B1D。

步骤6、选取以该雷达地理坐标(东经118.698°北纬32.191°)为中心，230km为半径范围内的地面雨量观测设备(包括常规地面观测站、加密自动气象站、单要素雨量计)，以该雷达探测时间为基准，提取这些设备在此时间后的1小时内累积观测降水量，形成数据集ORF。ORF[c]则表示第c个地面雨量观测设备的1小时雨量观测记录。不妨假设当前ORF数据集的记录总数为312个，即雷达有效覆盖区域内存在312个地面雨量观测设备。

步骤7、根据ORF中每个设备的经纬度信息，以及当前雷达的地理坐标信息，提取各个设备所对应地理位置的CAPPI_MAX值。为便于表述，不妨定义CM为CAPPI_MAX的一个子集，CM中的每个值CM[c]为与ORF[c]最相邻地理位置的回波强度值。显然，CM数据集中的记录数与ORF相同。设置两个递增变量PA和Pb，其中PA∈[100,400]，以10为步长递增，Pb∈[1.0,2.0]，以0.1为步长递增，分别计算ORF与CM中各记录项的误差平方和，计算方法为：

$D(PA,Pb)=\underset{c=1}{\overset{MSum}{\&Sigma;}}{(ORF\&lsqb;c\&rsqb;-\sqrt[Pb]{\frac{{10}^{\frac{CM\&lsqb;c\&rsqb;}{10}}}{PA}})}^2$

步骤8、对每一组PA和Pb，都可计算出一个D(PA,Pb)，记录下D(PA,Pb)最小时的一组PA和Pb，记为Para_A和Para_b。

步骤9、将计算得到的FeaStr和Para_A和Para_b，添加到数据库的FeatureInfo数据表中，该数据表的结构如表1所示。

表1

上述步骤1-9的过程如图1所示。

步骤10、依据上述步骤1-步骤9，将本实施例所选用的雷达和地面雨量观测的数据逐一进行计算并入库。完成后，即得到该部雷达的反射率特征库。

建立好雷达反射率特征库后，就可以在实际业务应用中，将实时、最新的雷达基数据文件通过反射率特征提取、特征匹配，找到最佳的Z-R关系的参数A和参数b，从而进行雷达降水估测。仍以这部南京的雷达为例，具体流程如下：

步骤11、将最新接收到的雷达基数据文件按上述步骤1-步骤5，计算出回波特征字符串，记为FeaStrCur。

步骤12、从雷达反射率特征库中遍历具有相同RadarID(即这部南京的雷达)的FeaStr，形成数据集FS，FS[c]表示FS中的c条记录。将FS中的逐条记录FS[c]与FeaStrCur作相关性计算，计算方法为：

$R\left(c\right)=\frac{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}\left(FeaStrCur\right(k)-\stackrel{\&OverBar;}{FeaStrxCur})\&CenterDot;(FS\&lsqb;c\&rsqb;(k)-\stackrel{\&OverBar;}{FS\&lsqb;c\&rsqb;})}{\sqrt{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}{\left(FeaStrCur\right(k)-\stackrel{\&OverBar;}{FeaStrCur})}^2}\&CenterDot;\sqrt{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}{(FS\&lsqb;c\&rsqb;(k)-\stackrel{\&OverBar;}{FS\&lsqb;c\&rsqb;})}^2}}$

其中，和分别表示FeaStrCur与FS[c]中所有记录的算术平均值。由于FeaStr字符串的记录方式是每2个字节为一个特征值，用FS[c](k)表示第k项记录值，且该值是16进制的数值，因此，在进行相关性R(c)计算前，先将数据FS[c](k)和FeaStrCur(k)转换为10进制，再进行计算。

步骤13、找到R(c)值最大时FS[c]所对应的记录，该记录中对应的Para_A和Para_b即为当前雷达资料估测降水时Z-R关系的最优参数A和参数b，由此计算出的降水量具有更高的准确性。需要补充说明的是，考虑到本方法投入实施前期，雷达反射率特征库中的特征信息记录可能并不完善，又或者当前天气现象和雷达回波信息比较特殊，从雷达反射率特征库中匹配计算出的相关系数R(c)都非常低，致使特征库中的参数A和参数b不适用于当前Z-R关系。因此，限定当最大的R(c)<0.5时，Z-R关系的参数A和参数b改为使用预先设定好的参数值，如Z＝300R^1.4、Z＝230R^1.25或Z＝200R^1.4等。

上述步骤11-13的过程如图2所示。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (6)

1.基于天气雷达反射率特征匹配的降水估测方法，其特征在于，包括以下步骤：

步骤一、构建天气雷达反射率特征库：

(1)提取天气雷达基数据中的基本反射率数据，并将基本反射率数据从原来的极坐标形式转换为平面直角坐标系形式；分别计算p千米和q千米高度处的等高平面位置显示图CAPPI_p和CAPPI_q，CAPPI_p和CAPPI_q的结构相同，均为N×N的二维网格，求取CAPPI_p和CAPPI_q各网格点在同一水平位置的最大回波强度值CAPPI_MAX，即CAPPI_MAX(x,y)＝Max(CAPPI_p(x,y),CAPPI_q(x,y))，其中x、y分别表示平面直角坐标系下的横坐标和纵坐标，Max为取最大值函数；

(2)根据得到的最大回波强度值CAPPI_MAX(x,y)，计算雷达回波的灰度图像Rad(x,y)：

$Rad\left(x,y\right)=\left\{\begin{array}{cc}0,& CAPPI\_MAX\left(x,y\right)<20\\ GAPPI\_MAX\left(x,y\right),& 20\&le;CAPPI\_MAX\left(x,y\right)\&le;64\\ 64,& CAPPI\_MAX\left(x,y\right)>64\end{array}\right.;$

(3)计算灰度图像Rad(x,y)的一阶梯度值，从而强化回波图像的变化特征；

(4)新建一个M×M的雷达回波特征矩阵Fea，M<N/2，Fea中每个单元的值计算方式为：

Fea(x',y')＝Hex(Fea'(x',y'))

${\mathrm{Fea}}^{,}(x^{,},y^{,})=\frac{\underset{k_0=0}{\overset{STEP-1}{\&Sigma;}}\underset{k_1=0}{\overset{STEP-1}{\&Sigma;}}Deg(x^{,}\&times;STEP+k_0,y^{,}\&times;STEP+k_1)}{STEP\&times;STEP\&times;4}$

其中，y'∈[0,M-1]，Hex为对括号中的值作16进制转换的函数；

(5)将得到的雷达回波特征矩阵Fea，从Fea(0,0)到Fea(M-1,M-1)逐个添加到一个字符串变量FeaStr中；

(6)选取雷达有效探测范围内的地面雨量观测设备，以当前雷达探测时间为基准，提取这些地面雨量观测设备在基准后的x小时内的累积观测降水量，形成数据集ORF，其中ORF[c]表示第c个地面雨量观测设备的x小时降水量，x∈[1,24]，将当前雷达有效探测范围内地面雨量观测设备的总数记为MSum，即数据集ORF的大小为MSum，则c∈[0,MSum-1]；

(7)根据数据集ORF中每个地面雨量观测设备的经纬度信息，以及当前雷达的地理坐标信息，提取各个地面雨量观测设备所对应地理位置的CAPPI_MAX值；定义CM为CAPPI_MAX的一个子集，CM中的每个值CM[c]为与ORF[c]地理位置最相邻的网格点的回波强度值；设置两个递增变量PA和Pb，其中PA∈[100,400]，以10为步长递增，Pb∈[1.0,2.0]，以0.1为步长递增，分别计算ORF与CM中各记录项的误差平方和：

$D\left(PA,Pb\right)=\underset{c=1}{\overset{MSum}{\mathrm{\&Sigma;}}}{\left(ORF\&lsqb;c\&rsqb;-\sqrt[Pb]{\frac{{10}^{\frac{CM\&lsqb;c\&rsqb;}{10}}}{PA}}\right)}^2;$

(8)记录下所得D(PA,Pb)取值最小时的一组PA和Pb，分别记为Para_A和Para_b；

(9)将步骤(5)中得到的FeaStr和步骤(8)得到的Para_A和Para_b添加到数据库；

(10)依据上述步骤(1)-(9)，将有资料以来，雷达探测数据和地面雨量观测数据逐一进行计算并入库，从而构建起天气雷达反射率特征库；

步骤二：将实况天气雷达反射率数据与建立的天气雷达反射率特征库进行匹配：

(11)将需要匹配的天气雷达反射率数据按照上述步骤(1)-(5)，计算出能够反映该雷达回波特征的字符串FeaStrCur；

(12)从天气雷达反射率特征库中遍历同一部天气雷达的FeaStr，形成数据集FS，FS中的每条记录FS[c]与FeaStrCur作相关性计算，计算方法为：

$R\left(c\right)=\frac{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}\left(FeaStrCur\right(k)-\stackrel{\&OverBar;}{FeaStrCur})\&CenterDot;(FS\&lsqb;c\&rsqb;(k)-\stackrel{\&OverBar;}{FS\&lsqb;c\&rsqb;})}{\sqrt{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}{\left(FeaStrCur\right(k)-\stackrel{\&OverBar;}{FeaStrCur})}^2}\&CenterDot;\sqrt{\underset{k=0}{\overset{M^2-1}{\&Sigma;}}{(FS\&lsqb;c\&rsqb;(k)-\stackrel{\&OverBar;}{FS\&lsqb;c\&rsqb;})}^2}}$

其中，和分别表示FeaStrCur与FS[c]中所有记录的算术平均值，FS[c](k)表示FS[c]的第k项记录值；

(13)找到R(c)取值最大时，FS[c]所对应的数据库记录，该记录中对应的Para_A和Para_b即为当前天气雷达资料下Z-R关系的最优参数A和参数b，并据此估测降水。

2.根据权利要求1所述基于天气雷达反射率特征匹配的降水估测方法，其特征在于：在步骤(1)中，取p＝1.5，q＝3.0。

3.根据权利要求1所述基于天气雷达反射率特征匹配的降水估测方法，其特征在于：在步骤(3)中，灰度图像Rad(x,y)最边缘一圈的各像素的一阶梯度为0，其余像素的一阶梯度采用Sobel检测算子来计算。

4.根据权利要求3所述基于天气雷达反射率特征匹配的降水估测方法，其特征在于：横向Sobel检测算子纵向Sobel检测算子

5.根据权利要求1所述基于天气雷达反射率特征匹配的降水估测方法，其特征在于：在步骤(5)中，字符串变量FeaStr中每2个字节代表1个Fea(x’,y’)。

6.根据权利要求1所述基于天气雷达反射率特征匹配的降水估测方法，其特征在于：在步骤(13)中，若找到的最大R(c)<0.5时，则Z-R关系中的最有参数A和参数b的取值改用预先设定值。