一种大规模降雨监测雷达校验方法与流程

1.本发明申请为申请日2021年11月26日，申请号为：202111419496.2，名称为“一种大规模降雨监测异常站点筛查方法”的发明专利申请的分案申请。
2.本发明涉及一种大规模降雨监测异常站点筛查方法，属于气象水文领域，主要用于为暴雨洪水风险提示和预警提供准确、可靠的降雨监测信息。


背景技术：

3.降雨监测是水文监测的重要组成部分，是暴雨洪水灾害防御工作的耳目和参谋。21世纪以来，水利部门加大了对自动监测站建设的支持力度，特别是经过山洪灾害防治项目建设，全国山洪灾害自动监测站点达13.2万个，自动雨量站网平均密度为38km2/站，是2006年(6000站)的22倍，最小报汛时段缩短到10min，数据信息量增加100余倍，极大地缩小了降雨监测盲区，有力支撑了水旱灾害防御工作。但由于部分测站建设标准偏低，位于山丘区的测站运维困难，数据质量难以得到保障，常常出现冒大数、缺测、漏测等情况，且测站出现问题具有很强的随机性，完全抛弃某一测站不切实际。
4.为了站点的监测数据得到有效利用，需要从众多降雨监测站点中找出不同时段内监测数据准确的站点，剔除数据质量有问题的站点，但目前大多数识别方法具有很强的主观性，误判的可能性很大，业务当中应用时存在较大风险，一是难以识别出确有异常的站点，二是容易将正确站点错误认为是异常站点，与开展异常站点识别的初衷相悖。对于极个别异常识别准确性较高的方法，需要进行大量的距离计算与比较，对计算资源要求极高，而事实上业务应用对异常站点识别方法的稳定性、可靠性和高效性均提出很高的要求。提出一个能够业务化应用的大规模异常站点识别方法极为迫切。


技术实现要素：

5.本发明提出了一种大规模降雨监测异常站点筛查方法，其解决的技术问题是基于hampel法、格拉布斯准则、周边测站分析法和雷达辅助校验等方法，建立递进式异常站点筛查体系，并通过k-d tree(k-dimension tree)高级数据结构和并行计算方法提高计算效率，为大规模降雨监测数据异常识别、充分利用雨量监测站有效信息提供可靠的方法。
6.为了解决上述存在的技术问题，本发明采用了以下方案：
7.一种大规模降雨监测异常站点筛查方法，包括以下步骤：
8.步骤1、采用hample法通过某一测站降雨时间序列查找出测站异常年份，在时间维度上初步判定基准站；
9.步骤2、采用改进的格拉布斯准则进行空间维度判别，通过异常年份测站的周边测站雨量确定是否为异常，再次判别基准站；
10.步骤3、确定基准站后采用周边测站分析法基于逐小时降雨监测数据判别异常站点；
11.步骤4、雷达辅助校验异常站点。
12.建立了“年尺度初筛—小时尺度细判—小时尺度验证”的多时空异常站点筛查体系。
13.进一步，步骤1、2提出了以年降雨序列为数据基础的“时空维度 先后规则”的基准站判别方法，即从大规模测站中，先通过某一测站降雨时间序列查找出测站的异常年份，再通过周边测站雨量确定其是否为异常，能够适应我国大多数测站年观测序列短的现状，提高了准确度。
14.进一步，步骤1中的采用hampel法初步筛选出雨量监测较为稳定、数据质量相对较高的测站，首先是基于长序列测站的年雨量值，利用hampel法对单一测站监测数据的异常年份进行识别。hampel法可用于异常极端值的判别，其基本原理是对给定的数据集假设一个分布和概率模型，然后根据假设采用不一致检验对数据系列进行处理，其公式如下：
[0015][0016]
其中，xi为数据序列x中的某一值，median为x的中位数，mad(median absolute deviation)为数据集y的中位数，x＝{x1,x2,
…
xn}，是测站年降雨量数据序列，y＝{y1,y2,
…
,yn}＝{x
1-median,x
2-median,
…
,x
n-median}。
[0017]
进一步，步骤1采用hampel法进行时间维度判别，对于年序列小于10的测站，阈值设定为2.24。当zi(i＝1,2,
…
,n)值大于2.24时，则判定xi为异常点，i为该测站的异常年份。
[0018]
步骤1可以根据研究区面积、测站分布情况因素，在考虑待测站点与周边站点距离的基础上，通过经纬度将目标分区，计算机程序在查找20km范围内测站时，仅在其中一个区域内查找即可,目标分区的区域同时进行搜索对比,已经判断为异常的测站在后续判别待测站时不作为周边测站参与比较，减少部分异常值对判别效果的不良影响。
[0019]
考虑到大量测站建设年代较近，数据序列较短，仅通过hampel法从时间维度上判定测站监测数据的异常年份，可靠性还不够。
[0020]
进一步，步骤2采用改进的格拉布斯准则进行空间维度判别，以监测出异常年份的站点为中心，20km为半径划定影响区域，对测站的异常年份进行验证。
[0021]
进一步，步骤2中格拉布斯准则适用于测量次数较少的情况(3≤n《100)，可一次性求出多个异常值。其基本原理是判断可疑值取舍的过程中，将正态分布中平均值和方差这两个最重要的参数引进来，有利于提高判断的准确性。考虑到本文单站年累积雨量的时间序列较短，为了提高异常值判定的准确性，需借助周边测站从空间维度再做判断，格拉布斯准则判别法较为适用。
[0022]
改进的格拉布斯准则是将原准则公式中的平均值用中位数代替，可有效消除同侧异常值的屏蔽效应，是更为稳健的处理方法，其方法如下:
[0023]
先将通过hampel法初步判定为某年异常的测站选出，在其周围以20km为半径划定区域，区域内测站约50个左右，区域内所有站点的年雨量值构成样本，通过从小至大排序为样本序列x＝(x1,x2,xn)，设定临界系数g(a,n)的值g0，g0通过查临界值表获得，a为显著性水平，取a为0.05，并计算g1，gn：
[0024]
g1＝(x
中-x1)/σ
[0025]
gn＝(x
n-x
中
)/σ
[0026]
其中，n为测站数量，x
中
为样本中位数，σ为标准差。
[0027]
若g1≥gn且g1＞g0，则判定x1为异常值并予以剔除；若gn≥g1且gn＞g0，则xn为异常值并予以剔除；若g1＜g0且gn＜g0，则不存在异常值。若存在异常值，则剔除之后用剩余站点的年雨量值重新计算，重复上述步骤，直到无异常值为止。
[0028]
为提高计算效率，缩短程序运行时对待判断测站与周边站点距离关系的判断所用时长，步骤2采用改进的格拉布斯准则进行空间维度判别时，根据研究区面积、测站分布情况因素，在考虑待测站点与周边站点距离的基础上，通过经纬度将目标分区，计算机程序在查找20km范围内测站时，仅在其中一个区域内查找,目标分区的区域同时进行搜索对比,已经判断为异常的测站在后续判别待测站时不作为周边测站参与比较，减少部分异常值对判别效果的不良影响。
[0029]
比如：通过经纬度将目标区福建省分为7个区域。程序在查找20km范围内测站时，仅在其中一个区域内查找，全省范围内7个区域同时进行搜索对比。已经判断为异常的测站在后续判别待测站时不作为周边测站参与比较，减少部分异常值对判别效果的不良影响。
[0030]
通过hampel法和格拉布斯准则共同判断后，则完成基准站初步判定。
[0031]
进一步，步骤3中采用周边测站分析法对逐小时的降雨监测数据进行异常识别。采用周边测站分析时，优先选用基准站与待评估测站进行同时段雨量比较，基准站距离较远的(超过某一阈值)经基准站评估后，确定为正常站点的非基准站与待评估测站进行同时段雨量比较。为了避免初步确定的基准站在某一时刻降雨监测出现问题，而依然被当做监测正确的测站，在采用基准站评估待评估测站时，待评估测站也包括基准站。通过比较待评估测站雨量与同时段周边基准站(或已评估合格的非基准站)的平均雨量，判断待评估测站是否异常。
[0032]
步骤3逐小时降雨监测数据判别时，采用周边测站分析法，仅对1h或3h、6h雨量超过10mm、12h雨量超过15mm或24小时雨量超过25mm的测站进行评估。评估时分别比较待测站不同时间段1h、3h、6h、12h、24h雨量值与周边站点相应时段的平均雨量值，当雨量值相差超过一个等级时，则判定为异常站点。其中，雨量等级划分如下表所示：
[0033][0034][0035]
为选取合适的评估范围，分别对距离待评估测站5～30km的范围进行了测试，当周边距离为5km、10km、15km、20km、25km、30km时，待测站的周边站点平均数量分别为4个、13个、30个、50个、78个、115个。对异常站点识别准确率及计算时长进行了统计。
[0036]
步骤3逐小时降雨监测数据判别时，周边分析法分析范围的最佳半径r与监测站网密度p关系为p＝0.0267r2 0.4667r 12，半径取值范围为[5km,30km]，p取值范围为[15km2/站，50km2/站]。考虑到福建省雨量监测站网密度约为25km2/站，且准确率在15km时达到最
大，计算时长相对适中，故采用周边测站分析时，分析范围的半径设为15km，分析范围内站点平均数量30个。
[0037]
福建省雨量监测站数量多，密度大，降雨监测数据量庞大。为了能够达到实时计算并判断异常站点的目的，研究采用k-dtree(k-dimension tree)高级数据结构和并行计算方法，大幅提高了计算效率。经测试，全省站点进行一次异常识别，计算时长约5～8min。
[0038]
步骤1-步骤3，以年降雨序列为数据基础的“时空维度 先后规则”的基准站判别方法，即从大规模测站中，先通过某一测站降雨时间序列查找出测站的异常年份，再通过周边测站雨量确定其是否为异常，
[0039]
通过对目标区域进行分区，在单个区域内采用k-d tree(k-dimension tree)高级数据结构和并行计算方法，大幅提高了计算效率。
[0040]
步骤4中通过雷达辅助校验，对筛选出的异常站点做进一步验证，确定异常站点。经过基准站初步判定和周边测站分析后，已初步完成异常站点的筛查，但人工校验时发现，处于雨区与非雨区边界、雨强差异较大的雨区边界报讯正常的测站也容易被判定为异常站点，因此仍需要进一步验证筛查结果。尽管雷达降雨反演的精度受到反演算法等因素的影响，但依据反射率仍能判定雷达探测覆盖范围是否降雨或判断降雨量级，能够充分反映某一时段降雨的空间分布特征，故采用雷达辅助校验异常站点初步筛查结果是合理的。
[0041]
进一步，步骤4采用雷达辅助校验的方法，是对处于雨区与非雨区边界、雨强差异较大的雨区边界报讯正常的测站进行筛查，筛查条件有三，缺一不可，一是用雷达低层仰角反射率是否超过20dbz阈值，验证雨区与非雨区边界处测站判定结果；二是通过雷达基数据反演降雨强度，并与测站降雨量级进行对比，验证测站是否异常；三是通过对比反射率空间变化梯度值与雨量站观测梯度值，验证雨强差异较大的雨区边界处测站判定结果，当反射率空间变化梯度值超过雨量站观测梯度值1倍判定异常，当反射率空间变化梯度值小于雨量站观测梯度值1倍判定异常。
[0042]
基于hampel法、格拉布斯准则及周边测站法异常识别等方法进行联合判断，异常站点识别准确率高。本发明通过建立递进式异常站点筛查体系，提高了大规模异常站点识别准确率，对充分利用雨量监测站有效信息、剔除无效信息具有重要意义。
[0043]
本发明一种大规模降雨监测异常站点筛查方法，具有以下有益效果：
[0044]
(1)本发明中hampel法从时间维度上判定站点监测数据的异常年份，改进的格拉布斯准则借助周边测站从空间维度再做判断，以年降雨序列为数据基础的“时空维度 先后规则”，能够适应我国大多数测站年观测序列短的现状，提高了准确度。
[0045]
(2)本发明中采用的改进的格拉布斯准则进行空间维度判别时，根据面积、测站分布情况，对检测区域进行分区，查找周边测站时，仅在其中一个分区内查找，已判定为异常的测站在后续判别待测站时，不再作为周边测站参与比较，提高了计算效率及准确性。
[0046]
(3)本发明中基于hampel法、格拉布斯准则、周边测站分析法和雷达辅助校验等方法，建立了递进式异常站点筛查体系，可以提高大规模雨量监测站数据异常识别的准确性和稳定性，为暴雨洪水风险提示和预警提供准确、可靠的依据。
[0047]
(4)本发明中采用的雷达辅助校验方法，是对被误判为异常站的正常测站，处于雨区与非雨区边界、雨强差异较大的雨区边界报讯正常的测站进行筛查，筛查条件有三，缺一不可，提高了异常识别的准确性。
[0048]
附图(表)说明
[0049]
图1：本发明中一种大规模降雨监测异常站点筛查方法流程图。
[0050]
图2：本发明中距离与计算时长/准确率的对应关系。
[0051]
图3：本发明中2015～2020年异常站点数量及占比情况。
[0052]
图4：本发明中6月27、28日4个时刻异常站点数量及占比情况。
具体实施方式
[0053]
本发明所采用的技术方案是基于hampel法、格拉布斯准则、周边测站分析法和雷达辅助校验等方法，建立递进式异常站点筛查体系，并通过k-d tree(k-dimension tree)高级数据结构和并行计算方法提高计算效率，为大规模降雨监测数据异常识别、充分利用雨量监测站有效信息提供可靠的方法。
[0054]
一种大规模降雨监测异常站点筛查方法，包括以下步骤：
[0055]
步骤1、采用hample法通过某一测站降雨时间序列查找出测站异常年份，在时间维度上初步判定基准站；
[0056]
步骤2、采用改进的格拉布斯准则进行空间维度判别，通过异常年份测站的周边测站雨量确定是否为异常，再次判别基准站；
[0057]
步骤3、基准站确定后采用周边测站分析法基于逐小时降雨监测数据判别异常站点；
[0058]
步骤4、雷达辅助校验异常站点。
[0059]
建立了“年尺度初筛—小时尺度细判—小时尺度验证”的多时空异常站点筛查体系。
[0060]
进一步，步骤1、2提出了以年降雨序列为数据基础的“时空维度 先后规则”的基准站判别方法，即从大规模测站中，先通过某一测站降雨时间序列查找出测站的异常年份，再通过周边测站雨量确定其是否为异常，能够适应我国大多数测站年观测序列短的现状，提高了准确度。
[0061]
进一步，步骤1中的采用hampel法初步筛选出雨量监测较为稳定、数据质量相对较高的测站，首先是基于长序列测站的年雨量值，利用hampel法对单一测站监测数据的异常年份进行识别。hampel法可用于异常极端值的判别，其基本原理是对给定的数据集假设一个分布和概率模型，然后根据假设采用不一致检验对数据系列进行处理，其公式如下：
[0062][0063]
其中，xi为数据序列x中的某一值，median为x的中位数，mad(median absolute deviation)为数据集y的中位数，x＝{x1,x2,
…
xn}，是测站年降雨量数据序列，y＝{y1,y2,
…
,yn}＝{x
1-median,x
2-median,
…
,x
n-median}。
[0064]
步骤1采用hampel法进行时间维度判别，对于年序列小于10的测站，阈值设定为2.24。当zi(i＝1,2,
…
,n)值大于2.24时，则判定xi为异常点，i为该测站的异常年份。
[0065]
步骤1可以根据研究区面积、测站分布情况因素，在考虑待测站点与周边站点距离的基础上，通过经纬度将目标分区，计算机程序在查找20km范围内测站时，仅在其中一个区
域内查找即可,目标分区的区域同时进行搜索对比,已经判断为异常的测站在后续判别待测站时不作为周边测站参与比较，减少部分异常值对判别效果的不良影响。
[0066]
考虑到大量测站建设年代较近，数据序列较短，仅通过hampel法从时间维度上判定测站监测数据的异常年份，可靠性还不够。
[0067]
进一步，步骤2采用改进的格拉布斯准则进行空间维度判别，以监测出异常年份的站点为中心，20km为半径划定影响区域，对测站的异常年份进行验证。
[0068]
进一步，步骤2中格拉布斯准则适用于测量次数较少的情况(3≤n《100)，可一次性求出多个异常值。其基本原理是判断可疑值取舍的过程中，将正态分布中平均值和方差这两个最重要的参数引进来，有利于提高判断的准确性。考虑到本文单站年累积雨量的时间序列较短，为了提高异常值判定的准确性，需借助周边测站从空间维度再做判断，格拉布斯准则判别法较为适用。改进的格拉布斯准则是将原准则公式中的平均值用中位数代替，可有效消除同侧异常值的屏蔽效应，是更为稳健的处理方法，其方法如下:
[0069]
先将通过hampel法初步判定为某年异常的测站选出，在其周围以20km为半径划定区域，区域内测站约50个左右，区域内所有站点的年雨量值构成样本，通过从小至大排序为样本序列x＝(x1,x2,xn)，设定临界系数g(a,n)的值g0，g0通过查临界值表获得，a为显著性水平，取a为0.05，并计算g1，gn：
[0070]
g1＝(x
中-x1)/σ
[0071]
gn＝(x
n-x
中
)/σ
[0072]
其中，n为测站数量，x
中
为样本中位数，σ为标准差。
[0073]
若g1≥gn且g1＞g0，则判定x1为异常值并予以剔除；若gn≥g1且gn＞g0，则xn为异常值并予以剔除；若g1＜g0且gn＜g0，则不存在异常值。若存在异常值，则剔除之后用剩余站点的年雨量值重新计算，重复上述步骤，直到无异常值为止。
[0074]
为提高计算效率，缩短程序运行时对待判断测站与周边站点距离关系的判断所用时长，步骤2采用改进的格拉布斯准则进行空间维度判别时，根据研究区面积、测站分布情况因素，在考虑待测站点与周边站点距离的基础上，通过经纬度将目标分区，计算机程序在查找20km范围内测站时，仅在其中一个区域内查找,目标分区的区域同时进行搜索对比,已经判断为异常的测站在后续判别待测站时不作为周边测站参与比较，减少部分异常值对判别效果的不良影响。
[0075]
比如：根据经纬度将目标区域福建省分为7个区域。程序在查找20km范围内测站时，仅在其中一个区域内查找，全省范围内7个区域同时进行搜索对比。已经判断为异常的测站在后续判别待测站时不作为周边测站参与比较，减少部分异常值对判别效果的不良影响。
[0076]
通过hampel法和格拉布斯准则共同判断后，则完成基准站初步判定。
[0077]
进一步，步骤3中采用周边测站分析法对逐小时的降雨监测数据进行异常识别。采用周边测站分析时，优先选用基准站与待评估测站进行同时段雨量比较，基准站距离较远的(超过某一阈值)经基准站评估后，确定为正常站点的非基准站与待评估测站进行同时段雨量比较。为了避免初步确定的基准站在某一时刻降雨监测出现问题，而依然被当做监测正确的测站，在采用基准站评估待评估测站时，待评估测站也包括基准站。通过比较待评估测站雨量与同时段周边基准站(或已评估合格的非基准站)的平均雨量，判断待评估测站是
否异常。
[0078]
步骤3逐小时降雨监测数据判别时，采用周边测站分析法，仅对1h或3h、6h雨量超过10mm、12h雨量超过15mm或24小时雨量超过25mm的测站进行评估。评估时分别比较待测站不同时间段1h、3h、6h、12h、24h雨量值与周边站点相应时段的平均雨量值，当雨量值相差超过一个等级时，则判定为异常站点。其中，雨量等级划分如下表所示：
[0079][0080]
为选取合适的评估范围，分别对距离待评估测站5～30km的范围进行了测试，当周边距离为5km、10km、15km、20km、25km、30km时，待测站的周边站点平均数量分别为4个、13个、30个、50个、78个、115个。对异常站点识别准确率及计算时长进行了统计。
[0081]
步骤3逐小时降雨监测数据判别时，周边分析法分析范围的最佳半径r与监测站网密度p关系为p＝0.0267r2 0.4667r 12，半径取值范围为[5km,30km]，p取值范围为[15km2/站，50km2/站]。考虑到福建省雨量监测站网密度约为25km2/站，且准确率在15km时达到最大，计算时长相对适中(图2)，故采用周边测站分析时，分析范围的半径设为15km，分析范围内站点平均数量30个。
[0082]
福建省雨量监测站数量多，密度大，降雨监测数据量庞大。为了能够达到实时计算并判断异常站点的目的，研究采用k-dtree(k-dimension tree)高级数据结构和并行计算方法，大幅提高了计算效率。经测试，全省站点进行一次异常识别，计算时长约5～8min。
[0083]
进一步，步骤4中通过雷达辅助校验，对筛选出的异常站点做进一步验证，确定异常站点。经过基准站初步判定和周边测站分析后，已初步完成异常站点的筛查，但人工校验时发现，处于雨区与非雨区边界、雨强差异较大的雨区边界报讯正常的测站也容易被判定为异常站点，因此仍需要进一步验证筛查结果。尽管雷达降雨反演的精度受到反演算法等因素的影响，但依据反射率仍能判定雷达探测覆盖范围是否降雨或判断降雨量级，能够充分反映某一时段降雨的空间分布特征，故采用雷达辅助校验异常站点初步筛查结果是合理的。
[0084]
进一步，步骤4采用雷达辅助校验的方法，是对处于雨区与非雨区边界、雨强差异较大的雨区边界报讯正常的测站进行筛查，筛查条件有三，缺一不可，一是用雷达低层仰角反射率是否超过20dbz阈值，验证雨区与非雨区边界处测站判定结果；二是通过雷达基数据反演降雨强度，并与测站降雨量级进行对比，验证测站是否异常；三是通过对比反射率空间变化梯度值与雨量站观测梯度值，验证雨强差异较大的雨区边界处测站判定结果，当反射率空间变化梯度值超过雨量站观测梯度值1倍判定异常；当反射率空间变化梯度值小于雨量站观测梯度值1倍判定异常
[0085]
本研究利用福建省5234个雨量站2015-2021年实测降雨资料，基于hampel法、格拉
布斯准则、周边测站分析法和雷达辅助校验等方法构建了递进式异常站点筛查体系，对降雨监测数据进行了异常识别。经过hampel法和格拉布斯准则判别法对2015—2020年异常站点进行判定发现，2015年的异常站点数量最多，占比11.5％，之后异常站点逐年减少，至2020年，异常站点数量占比仅5.18％(图3)。
[0086]
从通过周边测站法得出的异常结果中分别选取2016—2020年每年的6月、7月和8月的1日8点、10日14点、30日20点的异常值识别结果进行人工验证，异常识别正确率均为90％(表1)。
[0087]
表1基于周边测站分析的雨量异常值识别正确率(％)
[0088][0089]
为进一步验证方法的可用性，以2021年6月27—28日福建全省的雨量监测站进行了异常站点实时判别(图4)。考虑到福建全省在6月27日19:00-21:00和6月28日09:00-15:00降雨范围较广，报讯站较多，所以选择如下6个时刻的异常结果进行雷达回波的校验。通过雷达辅助校验法对被误判为异常站的正常测站进行筛查，雷达校验前异常识别结果的平均准确率为89％，经雷达校验后平均准确率提升至95％(表2)，这表明雷达辅助校验方法非常适用于应对处于雨区与非雨区边界、雨强差异较大的雨区边界的正常站点被错误判断为异常站点的情况。
[0090]
表2 6个时刻雷达校验前后福建省雨量监测站异常识别结果准确率(％)
[0091][0092]
基于hampel法、格拉布斯准则及周边测站法异常识别等方法进行联合判断，异常站点识别准确率高。本发明通过建立递进式异常站点筛查体系，提高了大规模异常站点识别准确率，对充分利用雨量监测站有效信息、剔除无效信息具有重要意义。
[0093]
上面结合附图对本发明进行了的描述，显然本发明的实现并不受上述方式的限制，只要采用了本发明的方法构思和技术方案进行的各种改进，或未经改进将本发明的构思和技术方案直接应用于其它场合的，均在本发明的保护范围内。