一种基于回波强度和回波顶高外推的短时定量降水预报方法
阅读说明：本技术 一种基于回波强度和回波顶高外推的短时定量降水预报方法 (Short-time quantitative rainfall forecasting method based on echo intensity and echo top height extrapolation ) 是由 邹海波 吴珊珊 单九生 吴文心 易雪婷 于 2021-09-01 设计创作，主要内容包括：本发明涉及大气科学领域,具体是一种基于回波强度和回波顶高外推的短时定量降水预报方法,用于基于雷达的短时定量降水预报。为了达到上述目的,本发明所采用的技术方案是：雷达基数据处理,剔除雷达回波观测中的非气象回波,并将数据从极坐标格式转换至笛卡尔直角坐标格式；雷达回波外推,对雷达回波强度和雷达回波顶高同时进行外推,获得两者的预报值；雷达降水反演,基于实时雷达回波强度、雷达回波顶高和降水观测资料,拟合出不同回波顶高的Z-R关系表达式；雷达降水预报,基于雷达回波强度和回波顶高的预报值,利用实时拟合出的不同回波顶高的Z-R关系表达式,对外推预报的回波强度进行降水反演,获得短时定量降水预报。(The invention relates to the field of atmospheric science, in particular to a short-time quantitative rainfall forecasting method based on extrapolation of echo intensity and echo top height, which is used for short-time quantitative rainfall forecasting based on radar. In order to achieve the purpose, the invention adopts the technical scheme that: radar base data processing, namely eliminating non-meteorological echoes in radar echo observation, and converting data from a polar coordinate format to a Cartesian rectangular coordinate format; extrapolating the radar echo, namely extrapolating the intensity of the radar echo and the height of the radar echo simultaneously to obtain predicted values of the intensity and the height of the radar echo; radar rainfall inversion, namely fitting Z-R relational expressions with different echo top heights based on real-time radar echo intensity, radar echo top height and rainfall observation data; and radar rainfall forecast, namely carrying out rainfall inversion on the echo intensity of the extrapolation forecast by utilizing Z-R relational expressions with different echo top heights, which are fit in real time, based on the forecast values of the radar echo intensity and the echo top heights, so as to obtain the short-time quantitative rainfall forecast.)
技术领域
本发明涉及大气科学领域，尤其是一种基于回波强度和回波顶高外推的短时定量降水预报方法。
背景技术
多普勒天气雷达资料有很高的时空分辨率，在中尺度气象业务和研究中具有其他观测资料无可比拟的优势。由于数值模式存在spin-up问题(即在非平衡初值或扰动条件下，模式需要进行调整而达到平衡态的过程)，短时(2小时内)降水的预报效果不佳。目前，短时降水的业务预报主要依赖于雷达回波的外推预报。
目前，基于雷达回波外推的短时降水预报方法主要都是先采用不同技术(如交叉相关跟踪法TREC、基于连续方程约束的COTREC、通过Barnes空间滤波处理的BTREC、光流法等)获得雷达回波的移动矢量。再在假定回波强度短时间内无明显变化的条件下，对雷达回波进行外推，获得未来2小时内不同时刻的回波预报值。最后再利用降水强度R与雷达回波强度Z的经验公式Z＝aRb(即Z-R关系)反算(反演)出未来2小时内不同时刻的雷达降水量，即预报的降水。
可见，基于雷达回波外推的短时降水预报精度，不仅受雷达回波外推的精度影响，还要受外推回波降水反演的精度影响。以往的雷达降水预报方法多集中在雷达回波外推精度的改进上，如TREC、COTREC、BTREC和光流法等，对外推回波的降水反演则多采用固定Z-R关系(即Z＝300R1.4)算法。然而，由于Z-R关系中的系数a和b主要由降水雨滴普的特征决定，且随天气系统、降水类型、地理位置、季节等因素的变化而变化。因此，即使在雷达回波外推预报精度较高的情况下，用固定Z-R关系也很难准确地获得不同地区不同季节的雷达降水预报值。近年来，有学者也基于Z-R关系中的系数a和b随季节、降水类型、地理位置等变化的特性，发展一些精度更高的Z-R关系算法，如动态Z-R关系算法、基于回波顶高分级的动态Z-R关系算法等，但这些算法不仅需要用最优法进行实时拟合，还需要额外的预报参数(如回波顶高等)，因此，很少被用于基于雷达的短时降水预报上。
发明内容
(一)解决的技术问题
本发明的目的是克服现有技术的不足，本发明提供了一种基于回波强度和回波顶高外推的短时定量降水预报方法，通过对雷达回波强度(可反映云团中水凝物的浓度)和雷达回波顶高(可反映云团内部的上升运动)的共同外推，获得两者的预报值，解决基于回波高分级动态Z-R关系算法在雷达降水预报上的应用问题，有利于提升短时降水的定量预报精度。
(二)技术方案
本发明的技术方案：一种基于回波强度和回波顶高外推的短时定量降水预报方法，用于基于雷达的短时定量降水预报。为了达到上述目的，本发明所采用的技术方案是：雷达基数据处理，剔除雷达回波观测中的非气象回波，并将数据从极坐标格式转换至笛卡尔直角坐标格式；雷达回波外推，对雷达回波强度和雷达回波顶高同时进行外推，获得两者的预报值；雷达降水反演，基于实时雷达回波强度、雷达回波顶高和降水观测资料，拟合出不同回波顶高的Z-R关系表达式；雷达降水预报，基于雷达回波强度和回波顶高的预报值，利用实时拟合出的不同回波顶高的Z-R关系表达式，对外推预报的回波强度进行降水反演，获得短时定量降水预报。
具体包括如下步骤：
S1：雷达基数据处理：对雷达基数据进行解码、质量控制和坐标转换，获得笛卡尔坐标的雷达回波资料。
步骤S2：雷达回波外推：基于网格化的雷达组合反射率数据，通过交叉相关跟踪(TREC)法获得雷达回波移动矢量(即TREC矢量)，再用Barnes空间滤波更重矢量中的噪声和不协调，获得BTREC矢量，最后再用后向外推方案对雷达回组合反射率和雷达回波顶高进行外推，获得未来2小时内雷达回波强度和回波顶高的预报值。
步骤S3：雷达降水反演(基于回波顶高分级的Z-R关系)：将经过质控处理和网格化处理后的雷达组合反射率和雷达回波顶高数据，用双线性插值法插值至降水观测站处。根据回波顶高，按1km的间距将雷达组合反射率和观测降水分成不同类，如回波顶高在0-1km、1-2km、2-3km、3-4km、4-5km、5-6km、6-7km、7-8km、8-9km、9-10km、10-11km、11-12km、12-13km、13-14km、14-15km、＞15km范围内的降水与组合反射率，再用最优法拟合出不同类的最优Z-R关系。最优法的公式为其中：n为观测降水的样本数，Gi为第i个观测站的观测降水量，为第i个观测站处的雷达反演降水量，Zi为第i个观测站的雷达组合反射率，aj＝1,2，……，1200，bk＝100,101，……，300，不断调整aj和bk，使CTF值最小，此时aj和bk对应的Z-R关系即为最优Z-R关系。
步骤S4：雷达降水预报。基于外推的未来2小时内不同时次的雷达组合反射率和回波顶高数据，利用当前时次(动态)拟合的不同回波顶高的Z-R关系式，计算出未来2小时内不同时刻不同地点的定量降水预报值。
本发明技术方案中，所述S1具体内容包括：
S11：雷达基数据解码：根据雷达的型号(如CINRADSA/SB/SC/CB/CC/CD等)及其基数据的存储格式，对雷达基数据进行解码，获得不同仰角的极坐标格式的雷达回波资料。
S12：雷达回波质量控制：先利用回波(反射率)垂直梯度RGDZ＝W(Z-Zup)≥20dBZ剔除雷达回波观测中的异常地物杂波，其中：权重系数W在雷达斜距为0-40km时等于1，在40-200km间线性递减至0，Z为当前仰角的雷达回波(dBZ)，Zup为上一层仰角的雷达回波(dBZ)。再利用模糊逻辑算法Px＝N/25≤0.75剔除孤立非气象回波，其中N为方位径向平面(即雷达在某仰角上的扫描平面)上x点周围5×5窗口中雷达有效观测的数量。
S13：组合反射率和回波顶高的计算：基于质量控制处理后的不同仰角的极坐标格式的雷达反射率资料，计算出组合反射率(不同仰角层的最大值)和雷达回波顶高(即，≥18dBZ反射率因子的最高高度)，雷达回波顶高(相对雷达站)的计算公式为h＝Rsinθ+R2cos2θ/(2Re)，其中R为雷达探测的斜距，Re＝8500km为标准大气下的等效地球半径，θ为雷达探测的≥18dBZ反射率因子的最高仰角。
S14：雷达资料的坐标转换：先以雷达测站为中心，构建一个600×600个格点、水平分辨率为1km×1km的笛卡尔直角坐标网格(如图2虚线网格所示)，再根据雷达波束水平距离计算公式其中hr为雷达站的海拔高度，计算出雷达斜距在水平面上的投影，即图2中的L。最后，再通过x＝L·sinβ和y＝L·cosβ将极坐标系(L，β)的雷达资料转至笛卡尔直角坐标系，其中β为雷达观测的方位角。
本发明技术方案中，所述S2具体内容包括：
S21：TREC矢量的计算。将t时刻网格化(格点数600×600，水平分辨率1km×1km)的雷达组合反射率数据分成40×40个没有重叠的小区域，其中每个小区域的格点数为15×15，并将每个小区域(初始小区域)的雷达反射率数据按先x方向系后y方向的方式转换成一个一维序列Z1。将Z1序列与t+Δt(默认为12分钟)时刻的雷达组合反射率数据中任意15×15个格点组成的Z2序列进行空间相关分析，其相关系数计算公式为并将相关系数最大的Z2序列对应的小区域作为初始小区域的目标区域，从初始小区域中心至目标区域中心的矢量便是TREC矢量，反映雷达回波的移动情况，如图3所示。
S22：BTREC矢量的计算。由于回波小尺度的变化及跟踪失败等原因，TREC矢量中会有一些噪声或不协调的矢量，显著地影响着外推回波的精度。将TREC矢量在x方向和y方向的分量分别作为变量F(x,y)输入公式和获得BTREC矢量(剔除了噪声和更正了不协调)的分量F1(x,y)。其中，x和y是雷达回波在笛卡尔直角坐标系中的位置，高斯参数和M为距离(x,y)点100km范围内TREC矢量的格点数，rk为点(x,y)和点(xk,yk)之间的距离，滤波参数G＝0.35和C＝300。
S23：雷达回波外推。先利用双线性插值法将15km×15km水平分辨率的BTREC矢量插值至雷达组合反射率和雷达回波顶高的网格(1km×1km)上，再用雷达回波后向外推法进行雷达组合反射率和回波顶高的外推，获得未来2小时内任意时刻(如逐6分钟)的雷达组合反射率和雷达回波顶高的外推预报值。雷达回波后向外推法为Zi,j(t+Δt)＝Zi-Δi,j-Δj(t)，其中Z为雷达组合反射率或回波顶高，i和j分别为某一格点在x和y方向的序号，Δt为外推预报时间(如6分钟，半小时或2小时等)，Δi＝int(u·Δt+0.5)为预报时间内该格点在x方向上的移动距离，Δj＝int(v·Δt+0.5)格为预报时间内该格点在y方向上的移动距离，u和v分别为BTREC矢量在x和y方向上的分量。
(三)有益效果
本发明的优点在于：基于雷达回波强度和回波顶高的共同外推，实现了在雷达降水预报中引入了能反映云团内部上升运动的回波顶高，有利于提升短时降水定量预报的精度。
附图说明
图1是本发明方法流程图；
图2是雷达极坐标资料网格化示意图；
图3是TREC矢量计算示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1
请参阅图1-图3，一种基于回波强度和回波顶高外推的短时定量降水预报方法，如图1所示，包括以下步骤：
S1：雷达基数据处理。对雷达基数据进行解码、质量控制和坐标转换，获得笛卡尔坐标的雷达回波资料，具体内容包括：
S11：雷达基数据解码。根据雷达的型号(如CINRAD SA/SB/SC/CB/CC/CD等)及其基数据的存储格式，对雷达基数据进行解码，获得不同仰角的极坐标格式的雷达回波资料。
S12：雷达回波质量控制。先利用回波(反射率)垂直梯度RGDZ＝W(Z-Zup)≥20dBZ剔除雷达回波观测中的异常地物杂波，其中权重系数W在雷达斜距为0-40km时等于1，在40-200km间线性递减至0，Z为当前仰角的雷达回波(dBZ)，Zup为上一层仰角的雷达回波(dBZ)。再利用模糊逻辑算法Px＝N/25≤0.75剔除孤立非气象回波，其中N为方位径向平面(即雷达在某仰角上的扫描平面)上x点周围5×5窗口中雷达有效观测的数量。
S13：组合反射率和回波顶高的计算。基于质量控制处理后的不同仰角的极坐标格式的雷达反射率资料，计算出组合反射率(不同仰角层的最大值)和雷达回波顶高(即，≥18dBZ反射率因子的最高高度)，雷达回波顶高(相对雷达站)的计算公式为h＝Rsinθ+R2cos2θ/(2Re)，其中R为雷达探测的斜距，Re＝8500km为标准大气下的等效地球半径，θ为雷达探测的≥18dBZ反射率因子的最高仰角。
S14：雷达资料的坐标转换。先以雷达测站为中心，构建一个600×600个格点、水平分辨率为1km×1km的笛卡尔直角坐标网格(如图2虚线网格所示，其中虚线网格为笛卡尔坐标系；椭圆粗实线为雷达极坐标；x、y分别为雷达极坐标观测在笛卡尔坐标的投影；A为雷达观测点；L为雷达斜距在平面上的投影；β为雷达观测的方位角。)，再根据雷达波束水平距离计算公式其中hr为雷达站的海拔高度，计算出雷达斜距在水平面上的投影，即图2中的L。最后，再通过x＝L·sinβ和y＝L·cosβ将极坐标系(L，β)的雷达资料转至笛卡尔直角坐标系，其中β为雷达观测的方位角。
S2：雷达回波外推。基于网格化的雷达组合反射率数据，通过交叉相关跟踪(TREC)法获得雷达回波移动矢量(即TREC矢量)，再用Barnes空间滤波更重矢量中的噪声和不协调，获得BTREC矢量，最后再用后向外推方案对雷达回组合反射率和雷达回波顶高进行外推，获得未来2小时内雷达回波强度和回波顶高的预报值，具体内容包括：
S21：TREC矢量的计算。将t时刻网格化(格点数600×600，水平分辨率1km×1km)的雷达组合反射率数据分成40×40个没有重叠的小区域，其中每个小区域的格点数为15×15，并将每个小区域(初始小区域)的雷达反射率数据按先x方向系后y方向的方式转换成一个一维序列Z1。将Z1序列与t+Δt(默认为12分钟)时刻的雷达组合反射率数据中任意15×15个格点组成的Z2序列进行空间相关分析，其相关系数计算公式为并将相关系数最大的Z2序列对应的小区域作为初始小区域的目标区域，从初始小区域中心至目标区域中心的矢量便是TREC矢量，反映雷达回波的移动情况，如图3所示。
S22：BTREC矢量的计算。由于回波小尺度的变化及跟踪失败等原因，TREC矢量中会有一些噪声或不协调的矢量，显著地影响着外推回波的精度。将TREC矢量在x方向和y方向的分量分别作为变量F(x,y)输入公式和获得BTREC矢量(剔除了噪声和更正了不协调)的分量F1(x,y)。其中，x和y是雷达回波在笛卡尔直角坐标系中的位置，高斯参数和分，M为距离(x,y)点100km范围内TREC矢量的格点数，rk为点(x,y)和点(xk,yk)之间的距离，滤波参数G＝0.35和C＝300。
S23：雷达回波外推。先利用双线性插值法将15km×15km水平分辨率的BTREC矢量插值至雷达组合反射率和雷达回波顶高的网格(1km×1km)上，再用雷达回波后向外推法进行雷达组合反射率和回波顶高的外推，获得未来2小时内任意时刻(如逐6分钟)的雷达组合反射率和雷达回波顶高的外推预报值。雷达回波后向外推法为Zi,j(t+Δt)＝Zi-Δi,j-Δj(t)，其中Z为雷达组合反射率或回波顶高，i和j分别为某一格点在x和y方向的序号，Δt为外推预报时间(如6分钟，半小时或2小时等)，Δi＝int(u·Δt+0.5)为预报时间内该格点在x方向上的移动距离，Δj＝int(v·Δt+0.5)格为预报时间内该格点在y方向上的移动距离，u和v分别为BTREC矢量在x和y方向上的分量。
S3：雷达降水反演(基于回波顶高分级的Z-R关系)。将经过质控处理和网格化处理后的雷达组合反射率和雷达回波顶高数据，用双线性插值法插值至降水观测站处。根据回波顶高，按1km的间距将雷达组合反射率和观测降水分成不同类，如回波顶高在0-1km、1-2km、2-3km、3-4km、4-5km、5-6km、6-7km、7-8km、8-9km、9-10km、10-11km、11-12km、12-13km、13-14km、14-15km、＞15km范围内的降水与组合反射率，再用最优法拟合出不同类的最优Z-R关系。最优法的公式为其中n为观测降水的样本数，Gi为第i个观测站的观测降水量，为第i个观测站处的雷达反演降水量，Zi为第i个观测站的雷达组合反射率，aj＝1,2，……，1200，bk＝100,101，……，300，不断调整aj和bk，使CTF值最小，此时aj和bk对应的Z-R关系即为最优Z-R关系。
S4：雷达降水预报。基于外推的未来2小时内不同时次的雷达组合反射率和回波顶高数据，利用前一个时次(动态)拟合的不同回波顶高的Z-R关系式，计算出未来2小时内不同时刻不同地点的定量降水预报值。
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的保护范围由所附权利要求及其等同物限定。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：测距方法和装置