基于天气雷达反射率特征匹配的降水估测方法
【专利摘要】本发明公开了基于天气雷达反射率特征匹配的降水估测方法，首先构建雷达反射率特征库，然后确立实况雷达数据与雷达反射率特征库的匹配规则，并在此基础上动态确立Z?R关系中最优的参数A和参数b。本发明在短临降水预报业务应用中实现Z?R关系的动态最优，进一步提高雷达估测降水的准确性。
【专利说明】
基于天气雷达反射率特征匹配的降水估测方法
技术领域
[0001] 本发明属于大气科学领域，特别涉及了基于天气雷达反射率特征匹配的降水估测 方法。
【背景技术】
[0002] 天气雷达（以下简称雷达)是强对流天气探测的重要手段，它能够探测其扫描半径 内的雨强和雨量分布情况，进一步可以估测总降水量和降水率，尤其在短时临近预报等业 务中已得到广泛应用。
[0003] 基于Z-R关系的计算是雷达估测降水量的关键方法，而估测降水量的准确率很大 程度上依赖于Z-R关系（Z = ARb)中参数A和参数b的选择。很多学者通过对地面实测降水量 与雷达回波的分析研究，得出特定区域和(或)特定时期内相对最优的参数A和参数b，由此 作为未来对同区域、同气象条件下雷达估测降水Z-R关系计算时的重要参考。
[0004] 在短时临近预报业务软件中，其参数A和参数b或是针对当地气象条件和降水类型 设定为一组相对最优值，或是按时期、按降水类型(如大陆强对流降水型、热带降水型等)设 定多组相对最优值。
[0005] 然而，由于雷达回波的时空变化性很大，从雷达回波图像上看，几乎每次降水过程 的回波图像都是独一无二的，其Z-R关系的最优参数A和参数b也并非固定不变的。尽管部分 气象业务软件中Z-R关系的参数是可调的，但在某一降水事件发生前，无论人工或是软件自 动地确定适当的参数都是非常困难的。因此，如何在短临降水预报业务应用中，实现Z-R关 系的动态最优，进一步提高雷达估测降水的准确性需要特别关注。
【发明内容】
[0006] 为了解决上述【背景技术】提出的技术问题，本发明旨在提供基于天气雷达反射率特 征匹配的降水估测方法，通过构建雷达反射率特征库，以及将反射率特征与雷达反射率特 征库进行匹配，动态确立Z-R关系中最优的参数A和参数b，提高了降水估测的准确性。
[0007] 为了实现上述技术目的，本发明的技术方案为：
[0008] 基于天气雷达反射率特征匹配的降水估测方法，包括以下步骤：
[0009] 步骤一、构建天气雷达反射率特征库：
[0010] (1)提取天气雷达基数据中的基本反射率数据，并将基本反射率数据从原来的极 坐标形式转换为平面直角坐标系形式;分别计算P千米和q千米高度处的等高平面位置显示 图?△??1_口和0厶?？1_9,0厶?？1_口和0厶??1_9的结构相同，均为~\袖勺二维网格，求取0厶?？1_口 和CAPPI_q各网格点在同一水平位置的最大回波强度值CAPPI_MAX，即CAPPI_MAX(x，y) = Max(CAPPI_p(x，y)，CAPPI_q(x，y))，其中x、y分别表示平面直角坐标系下的横坐标和纵坐 标，Max为取最大值函数；
[0011] (2)根据得到的最大回波强度值CAPPI_MAX (x，y)，计算雷达回波的灰度图像Rad (x,y)：
[0012]
[0013] (3)计算灰度图像Rad(x，y)的一阶梯度值，从而强化回波图像的变化特征；
[0014] (4)新建一个MXM的雷达回波特征矩阵？68，]\1〈1'1/2，？63中每个单元的值计算方式 为：
[0015] Fea(x'，y '）= Hex(Fea'（X，，y '））
[0016]
[0017] 其中，STEP = |三|，1'￡[0，-1]，7'￡[0，-1]此叉为对括号中的值作16进制转换 Μ 的函数；
[0018] (5)将得到的雷达回波特征矩阵Fea，从Fea(0,0)到Fea(M-l，M-l)逐个添加到一个 字符串变量FeaStr中；
[0019] (6)选取雷达有效探测范围内的地面雨量观测设备，以当前雷达探测时间为基准， 提取这些地面雨量观测设备在基准后的X小时内的累积观测降水量，形成数据集0RF，其中 〇RF[c]表示第c个地面雨量观测设备的X小时降水量，xe[l，24]，将当前雷达有效探测范围 内地面雨量观测设备的总数记为MSum，即数据集0RF的大小为MSum，则ce [0，MSum-l];
[0020] (7)根据数据集0RF中每个地面雨量观测设备的经炜度信息，以及当前雷达的地理 坐标信息，提取各个地面雨量观测设备所对应地理位置的CAPPI_MAX值；定义CM为CAPPI_ MAX的一个子集，CM中的每个值CM[ c ]为与0RF[ c ]地理位置最相邻的网格点的回波强度值； 设置两个递增变量PA和Pb，其中PAe [ 100，400]，以10为步长递增，Pb ￡[1.0,2.0]，以0.1为 步长递增，分别计算0RF与CM中各记录项的误差平方和：
[0021]
[0022] (8)记录下所得D (PA，Pb)取值最小时的一组PA和Pb，分别记为;
[0023] (9)将步骤(5)中得到的FeaStr和步骤(8)得到的添加到数据库；
[0024] (10)依据上述步骤（1)-(9)，将有资料以来，雷达探测数据和地面雨量观测数据逐 一进行计算并入库，从而构建起天气雷达反射率特征库；
[0025] 步骤二:将实况天气雷达反射率数据与建立的天气雷达反射率特征库进行匹配：
[0026] (11)将需要匹配的天气雷达反射率数据按照上述步骤（1)-(5)，计算出能够反映 该雷达回波特征的字符串FeaStrCur;
[0027] (12)从天气雷达反射率特征库中遍历同一部天气雷达的FeaStr，形成数据集FS， FS中的每条记录FS[c]与FeaStrCur作相关性计算，计算方法为：
[0028]
[0029] 其中，FeaStrCur和FS[c]分别表示FeaStrCur与FS[c]中所有记录的算术平均值， FS [ c ] (k)表示FS [ c ]的第k项记录值；
[0030] (13)找到R(c)取值最大时，FS[c]所对应的数据库记录，该记录中对应的 Para_b即为当前天气雷达资料下Z-R关系的最优参数A和参数b，并据此估测降水。
[0031] 进一步地，在步骤(1)中，取p = l .5，q = 3.0。
[0032]进一步地，在步骤(3)中，灰度图像Rad(x，y)最边缘一圈的各像素的一阶梯度为0， 其余像素的一阶梯度采用Sobel检测算子来计算。
[0033] 进一步地，横向Sobel检测算弓
纵向Sobel检测算子
。
[0034] 进一步地，在步骤(5)中，字符串变量FeaStr中每2个字节代表1个Fea(x'，y'）。
[0035] 进一步地，在步骤(13)中，若找到的最大R( c)〈0.5时，则Z-R关系中的最有参数A和 参数b的取值改用预先设定值。
[0036] 采用上述技术方案带来的有益效果：
[0037] (1)本发明对雷达回波进行了一阶梯度信息提取，强化雷达回波强弱区间的边缘 特征，为不同时刻雷达反射率的比较、匹配提供了分析的数据基础。
[0038] (2)本发明将雷达回波特征信息以字符串形式加以存储，便于数据库的实现以及 特征匹配时的快速计算。
[0039] (3)本发明中基于字符串和16进制数据格式的特征值的相关性计算，具有较高的 计算性能，有利于本方法业务化实施时，提高雷达估测降水预报的时效性。
【附图说明】
[0040] 图1是本发明中构建雷达反射率特征库的流程图，对应步骤1-9,图中实线表示程 序流，虚线表示数据流，当程序流与数据流重叠时，虚线省略；
[0041] 图2是本发明中实况雷达数据与雷达反射率特征库的匹配的流程图，对应步骤11-13，图中FS. Count表示数据集FS中的总记录数。
【具体实施方式】
[0042] 以下将结合附图，对本发明的技术方案进行详细说明。
[0043]本实施例选用2010年至2015年期间，南京地区多普勒天气雷达9层体积扫的基本 反射率数据。该雷达数据格式为CINRAD-SA，体扫模式为VCP21，时间分辨率为6分钟，径向分 辨率为lkm，所处地理位置为东经118.698°、北炜32.191°。剔除没有明显降水回波的雷达数 据，并对剩余待计算的雷达数据作必要的质量控制。然后，对这些雷达数据逐一作如下处 理：
[0044] 步骤1、不妨设当前待处理的雷达基数据文件名为RF.bin，提取该文件中的基本反 射率数据，将原有极坐标形式的数据结构转换为平面直接坐标形式。分别计算1.5km和 3.01〇11高度的04??1，得到04??1_15和04??1_30。再求取这两个高度，各网格点在同一水平位 置上的最大回波强度值，即CAPPI_MAX( X，y)=MaX(CAPPI_p(X，y)，CAPPI_q(X，y))，其中 X、y 分别表示平面直角坐标系下的横坐标和纵坐标，Max表示取最大值的函数。根据该部雷达的 性能技术指标，定义N=230,因此xe[0,229]，ye[0,229]。
[0045] 步骤2、由CAPP I_MAX，输出雷达回波的灰度图像，记为Rad，Rad中的值满足：
[0046]
[0047]步骤3、采用Sobel检测算子计算Rad的一阶梯度值，计算方法为：
[0048]
[0049] 其中，*表示卷积计算，展开后的计算方法形如：
[0050]
[0051]
[0052] 步骤4、新建一个MXM的雷达回波特征矩阵Fea，M=46，SETP = 5，Fea中每个单元的 值计算方式为：
[0053] Fea(x'，y '）= Hex(Fea'（X，，y '））
[0054]
[0055] 其中，SETP = 5，x' e[0,45]，y' e[0,45]，Hex表示一个函数，功能是对其括号中的 值作16进制转换。
[0056]步骤5、将Fea 中的各个元素，从Fea [ 0，0 ]、Fea [ 0，1 ]、Fea [ 0，2 ]依次到Fea [ 45，45 ] 逐个添加到一个字符串变量FeaStr中。例如，假设?6&[0,0]、？6&[0，1]、？6 &[0,2]至?6&[0， 5]值分别为：0、61、18、42、59、29，则？6&3钍中的记录形如 ：003012243810。
[0057] 步骤6、选取以该雷达地理坐标(东经118.698°北炜32.191°)为中心，230km为半径 范围内的地面雨量观测设备(包括常规地面观测站、加密自动气象站、单要素雨量计），以该 雷达探测时间为基准，提取这些设备在此时间后的1小时内累积观测降水量，形成数据集 ORFARFk]则表示第c个地面雨量观测设备的1小时雨量观测记录。不妨假设当前0RF数据 集的记录总数为312个，即雷达有效覆盖区域内存在312个地面雨量观测设备。
[0058]步骤7、根据0RF中每个设备的经炜度信息，以及当前雷达的地理坐标信息，提取各 个设备所对应地理位置的CAPPI_MAX值。为便于表述，不妨定义CM为CAPPI_MAX的一个子集， CM中的每个值CM[ c ]为与0RF[ c ]最相邻地理位置的回波强度值。显然，CM数据集中的记录数 与0RF相同。设置两个递增变量PA和Pb，其中PAe [ 100，400]，以10为步长递增，Pb e [ 1. 0， 2.0]，以0.1为步长递增，分别计算0RF与CM中各记录项的误差平方和，计算方法为：
[0059]
[0060] 步骤8、对每一组PA和Pb，都可计算出一个D(PA，Pb)，记录下D(PA，Pb)最小时的一 组PA和Pb，记为。
[0061 ] 步骤9、将计算得到的？6&31：1'和?&抑_4和?&抑_13，添加到数据库的？6&1：11代111；1^〇数 据表中，该数据表的结构如表1所;
[0062] 表 1
[0063]
[0064] 上述步骤1-9的过程如图1所示。
[0065] 步骤10、依据上述步骤1-步骤9,将本实施例所选用的雷达和地面雨量观测的数据 逐一进行计算并入库。完成后，即得到该部雷达的反射率特征库。
[0066] 建立好雷达反射率特征库后，就可以在实际业务应用中，将实时、最新的雷达基数 据文件通过反射率特征提取、特征匹配，找到最佳的Z-R关系的参数A和参数b，从而进行雷 达降水估测。仍以这部南京的雷达为例，具体流程如下：
[0067] 步骤11、将最新接收到的雷达基数据文件按上述步骤1-步骤5,计算出回波特征字 符串，记为FeaStrCur。
[0068] 步骤12、从雷达反射率特征库中遍历具有相同RadarlD(即这部南京的雷达）的 FeaStr，形成数据集FS，FS [ c ]表示FS中的c条记录。将FS中的逐条记录FS [ c ]与FeaStrCur作 相关性计算，计算方法为：
[0069]
[0070] 其中，FeaStrCur和F'S.[c]分别表示FeaStrCur与FS[c]中所有记录的算术平均值。 由于FeaStr字符串的记录方式是每2个字节为一个特征值，用FS[c](k)表示第k项记录值， 且该值是16进制的数值，因此，在进行相关性R(c)计算前，先将数据FS[c](k)和FeaStrCur (k)转换为10进制，再进行计算。
[0071 ] 步骤13、找到R( c)值最大时FS [ c ]所对应的记录，该记录中对应的 即为当前雷达资料估测降水时Z-R关系的最优参数A和参数b，由此计算出的降水量具有更 高的准确性。需要补充说明的是，考虑到本方法投入实施前期，雷达反射率特征库中的特征 信息记录可能并不完善，又或者当前天气现象和雷达回波信息比较特殊，从雷达反射率特 征库中匹配计算出的相关系数R(c)都非常低，致使特征库中的参数A和参数b不适用于当前 Z-R关系。因此，限定当最大的R(c)〈0.5时，Z-R关系的参数A和参数b改为使用预先设定好的 参数值，如 2 = 3001^^2 = 2301^25 或 2 = 2001^·4 等。
[0072] 上述步骤11-13的过程如图2所示。
[0073]以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是 按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围 之内。
【主权项】
1.基于天气雷达反射率特征匹配的降水估测方法，其特征在于，包括W下步骤： 步骤一、构建天气雷达反射率特征库： (1) 提取天气雷达基数据中的基本反射率数据，并将基本反射率数据从原来的极坐标 形式转换为平面直角坐标系形式；分别计算P千米和q千米高度处的等高平面位置显示图 CAPPI_p和CAPPI_q，CAPPI_p和CAPPI_q的结构相同，均为NXN的二维网格，求取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):(3) 计算灰度图像Rad(x，y)的一阶梯度值，从而强化回波图像的变化特征； (4) 新建一个MXM的雷达回波特征矩阵化a，M<N/2Jea中每个单元的值计算方式为：其中，x'e阳;，Μ - 1]巧'￡[〇，1-1]，胎义为对括号中的值作16进制转换的 函数； (5) 将得到的雷达回波特征矩阵化曰，从Fea(0,0巧化ea(M-l，M-l)逐个添加到一个字符 串变量化aStr中； (6) 选取雷达有效探测范围内的地面雨量观测设备，W当前雷达探测时间为基准，提取 运些地面雨量观测设备在基准后的X小时内的累积观测降水量，形成数据集ORF，其中ORF k]表示第C个地面雨量观测设备的X小时降水量，xe [1,24]，将当前雷达有效探测范围内 地面雨量观测设备的总数记为MS皿，即数据集ORF的大小为MS皿，则C e [0 ,MS皿-1 ]; (7) 根据数据集ORF中每个地面雨量观测设备的经缔度信息，W及当前雷达的地理坐标 信息，提取各个地面雨量观测设备所对应地理位置的CAPPI_MAX值;定义CM为CAPPI_MAX的 一个子集，CM中的每个值CMk]为与ORFk]地理位置最相邻的网格点的回波强度值;设置两 个递增变量PA和Pb，其中?4￡[100,400]，^10为步长递增，？6￡[1.0,2.0]，^0.1为步长递 增，分别计算0RF与CM中各记录项的误差平方和：(8) 记录下所得D (PA，Pb)取值最小时的一组PA和Pb，分别记为?日拘_4和?日拘_6; (9) 将步骤(5)中得到的化aStr和步骤(8)得到的？日拘_4和?日拘_6添加到数据库； (10) 依据上述步骤(1)-(9),将有资料w来，雷达探测数据和地面雨量观测数据逐一进 行计算并入库，从而构建起天气雷达反射率特征库； 步骤二:将实况天气雷达反射率数据与建立的天气雷达反射率特征库进行匹配： (11) 将需要匹配的天气雷达反射率数据按照上述步骤（1)-(5)，计算出能够反映该雷 达回波特征的字符串化aStr化。 (12) 从天气雷达反射率特征库中遍历同一部天气雷达的FeaS化，形成数据集FS，!^S中 的每条记录FSk ]与化aStr化r作相关性计算，计算方法为：其中，晚a託蛇此和RS [C]分别表示FeaStrCur与FS [ C ]中所有记录的算术平均值，FS k ] 化)表示FSk]的第k项记录值； (13) 找到R(c)取值最大时，FSU]所对应的数据库记录，该记录中对应的化^_八和 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检测算子5. 根据权利要求1所述基于天气雷达反射率特征匹配的降水估测方法，其特征在于:在 步骤(5)中，字符串变量化aStr中每2个字节代表1个化a (X '，y '）。6. 根据权利要求1所述基于天气雷达反射率特征匹配的降水估测方法，其特征在于:在 步骤(13)中，若找到的最大R(cK〇.5时，则Z-R关系中的最有参数A和参数b的取值改用预先 设定值。
【文档编号】G01S13/95GK105974418SQ201610540931
【公开日】2016年9月28日
【申请日】2016年7月8日
【发明人】王兴, 苗春生, 王坚红, 钱代丽, 王丽娟
【申请人】南京信息工程大学