金沙江下游遥测雨量站数据质量研究

刘秀林　张行南　方园皓　黄晴

摘要：随着水文自动测报技术的发展，越来越多的自记雨量站投入使用，而研究自记雨量站的数据质量对水利工程的运行具有重要意义。采用格拉布斯准则和K-medoids聚类方法对金沙江下游流域收集的2008～2015年遥测雨量站逐时降雨数据系列进行数据质量分析研究。结果显示，该方法在年累积雨量异常值的确定以及逐时降雨数据奇异点的寻找方面效果显著，能快速判别异常值和奇异点，且有统计学理论基础。所提出的方法为颗粒度越来越细的水文数据质量分析提供了一定参考依据。

关键词：遥测雨量;逐时降雨数据;格拉布斯准则;K-medoids聚类;奇异点;金沙江下游流域

中图法分类号：P33

文献标志码：A

DOI：10.16232/j.cnki.1001-4179.2019.03.023

随着我国社会经济水平的不断发展，水利基础工程设施建设的脚步也越来越快，截至2017年，长江流域的雨量站数量已达9959个，其中包括大量为水情自动测报系统服务的遥测雨量站。遥测雨量站的建设和应用加密了站网密度，满足了流域实时洪水预报的工程需求，在一定程度上为面雨量计算方法、流域洪水预报的研究提供了大量的数据支持。随着长江，上游水库工程的进一步规划建设，为满足水库的施工期洪水预报对预见期和预报精度以及水库建成后枢纽信息实时监测等计算要求，需进一步加密现有的雨量站网。由此可见，现阶段我国的雨量站数目仍将增加，同时，雨量数据的时序颗粒度也越来越细，目前的遥测雨量站数据储存皆是逐时雨量，一个雨量站每年约8760个值存储于数据库中。

雨量数据经遥测雨量站实时传输入数据库，大量的数据未经质量分析，给工程项目和科学研究带来了极大的不便。单站逐时降雨量时序数据不同于水位或者流量的时间系列，具有连续性的特点，在分析系列时可通过平均变化量（或者平均变化率）来确定异常点，逐时降雨量数据存在大量的零值，是离散的、随机的周期性过程。序列多会出现两类问题：模式异常和点异常。模式异常指在一条序列上与其他模式存在显著差异的、具有异常行为的模式;点异常指在某段时序区间内与其他序列点存在显著差异的、具有异常特征的序列点。这两类问题在单站逐时降雨量序列中具体表现形式分别为年累积雨量异常和区间奇异点。

本文以金沙江下游流域自建雨量站（指金沙江下游梯级水电站水情自动测报系统三期建设的站点，以及部分可直接或间接获取数据的信息共享站点）为例，提出面对大量数据时，降雨量数据质量分析的方法，以期为高效分析降雨量数据质量提供一定的参考依据。

1 研究区域与数据

1.1 研究区域

金沙江流域地处青藏高原、云贵高原和四川盆地西部边缘，呈东西短、南北长的狭长形状。金沙江下段自攀枝花站起，至四川省宜宾岷江口止，地理位置位于东经100°～105°、北纬24°～29°之间，河长782.4km，落差729m。该河段水量大、落差集中，是长江流域水能资源最丰富的河段，同时也是雨量站点较为密集的区域。

金沙江下游流域气候为亚热带季风湿润气候，多年平均降水量为893mm，水量丰沛稳定，暴雨多为两次以上的连阴雨天气形成，一次暴雨过程的历时约为3～6d，最大1d降雨量高达100～200mm。本次研究收集到金沙江下游流域95个自记雨量站点逐时雨量数据，站点分布情况见图1。

1.2 雨量站点信息

自记雨量站点信息及數据均为三峡集团公司梯调通信中心提供，因篇幅限制，各站点经纬度信息和高程信息均未列出，具体位置可参见图1。在95个雨量站中，大沙店、尼格、向家坝（专）和龙街（三）站点资料系列为2012～2015年，底坝和细沙站资料系列为2013～2015年，其余均为2008～2015年，数据总量为349651条。

2逐时降雨数据质量分析方法

在一定区域范围内，邻近雨量站之间的降雨量有一定的相关性，故在横向比较（单站各年份比较）的基础上，针对相应情况，还可进行纵向比较（邻近站点相应年份比较），以确定数据质量。WMO《水文实践指南》第1卷中提出：在温带、内陆热带山区，最合适的做法是按大约500m高差来规划高度带。针对金沙江下游流域雨量站分布位置高程差较大的特点，将与研究站点高程差500m作为搜索范围，寻找5个邻近站点进行对比分析。

2.1 累积雨量异常判别

累积雨量异常表现为相较于其他年份该年总量偏大或者偏小，即该年的逐时降雨量数据可能呈现系统性的偏大或者偏小，属于时间序列异常研究中的模式异常问题。这类问题需要特别警惕，应认真核实数据的正确性和可信度，因为它会影响水文研究中的水量平衡计算。对于该问题，考虑到本文单站年累积雨量统计值个数在3～8之间，故采用改进的格拉布斯准则来确定异常年份。

格拉布斯准则适用于测量次数较少的情况（3≤n<100），可一次性求出多个异常值。改进的格拉布斯准则是将原准则公式中的平均值用中位数代替，可有效消除同侧异常值的屏蔽效应，是更为稳健的处理方法。其判别方法如下。

先将样本从小至大排序为新的系列X=（x，x2，.x，），统计临界系数G（a，n）的值G，（查临界值表获得），然后计算G，G，：

公式

式中，a为显著性水平，n为测量次数，Xψ为样本中位数，σ为标准误差。

若G≥G，且G>Go，则x应予以剔除;若G，≥G且G>Go，则x。应予以剔除;若G

根据以上统计学方法，金沙江下游流域雨量站累积雨量异常分析总体思路为：先对单个站点计算各年的累积雨量，挑出异常年份;再用各年的累积雨量与邻近站点的同年累积雨量做格拉布斯准则分析。由于累积雨量在年际间相差可能较大，而邻近雨量站的雨量总体反映了丰水年、平水年或者枯水年，可作为较好的参考，故设定两者输出相同的异常年份则判定为需核查的数据。

2.2 奇异点识别

奇异点在本次研究中是指在某段时序区间内的极大值远远大于该段内的其他点，可能存在单点数据有误的随机误差。奇异点的存在并不是错误值，而是需要进一步核实是不是正确数据的极大值，以避免后续研究（如水文分析中的次洪参数率定等）中带来不必要的影响。

针对逐时降雨量数据存在大量的零值，且各个测量值之间离散相互独立的特点，采用时间序列异常研究中的基于特征空间的方法来识别奇异点：首先对时间序列进行分段，即分为非汛期和汛期，从分段中提取特征，然后在特征空间中应用无序数据集合中的异常点检测方法一聚类法来寻找奇异点，分段的思想考虑到了非汛期中的奇异点在汛期时段中显得平庸的特点，能够有效检测出各段奇异点。

聚类法是将数据集根据相似度划分成若干组的统计学方法，不同组中相似度低，相似度可用距离进行度量。K-medoids聚类法是系统聚类法中最为常用的一种，因其算法简单、收敛速度快、中心点明确以及局部搜索能力强的优点被应用到很多方面。K-medoids算法步骤如下：

（1） 针对数据集{y，，y2，.，yn}，适当选择k个样本作为初始聚类中心2，2，g};

（2） 对每个样本y;找到离它最近的聚类中心z。，并将其分配到z。所标明的类u;

（3） 更新每个类的中心：

公式

（5）如果D值收斂，则返回（z，z，*.，zn，U），并终止算法，否则转至步骤（2）。

K-medoids算法聚类的显著缺点是需提前指定分类数目，采用最优聚类数的评价指标Silhouette来确定分类数目。该指标反映了聚类结构的类内紧密性和类间分离性，既可用于估计最佳聚类数，也可用于评价聚类质量，Silhouette指标值在范围内变动，指标值越大表示聚类质量越好，最大值对应的类数为最佳聚类数。Silhouette指标值的计算公式为

公式

式中，a（i）是样本i与类内所有其他样本的平均距离，b（i）为样本i到其他每个类中样本平均距离的最小值。

初步分析金沙江下游流域雨量站逐时降雨量（也称降雨强度，mm/h）数据，可以发现大量的零值数据（无降雨），故初步处理应将零值去掉，形成无零值的数据系列文件;对新的数据系列用K-medoids算法进行分类。实践中仅将分类数k取为2或3，采用Silhouette指标值进行比较确定最优分类。小时降雨强度r的等级划分标准为：r<2.5mm（小雨），2.5mm≤r<8mm（中雨），8.0mm≤r<15mm（中雨），15mm≤r（暴雨）。由此可看出，逐时降雨量的特点是存在大量的较小值较多的中间值和很少的较大值，奇异点存在于较大值的聚类中10]。本次研究对较大值所在的聚类需要做进一步的分析。资料显示，最大1h降雨强度极少超过100mm，故在该聚类中直接将100mm以上的数据标记为奇异点。另外，经人工逐一在excel中绘制降雨量柱状图发现，一年中各个分段出现奇异点的概率不大，故将该聚类中的阈值设定为3，若该区段中较大值的个数不大于3个，那么皆标记为奇异点，若大于3个且小于100个，则采用格拉布斯准则来寻找异常数据（奇异点），找到异常数据则输出奇异点，寻找不到则输出无奇异点，认为存在100个及以上的较大值在同一聚类中，属于正常降雨数据。

根据金沙江下游流域降雨特点，将降雨量系列分为3个区段，1～3月和11～12月降雨较少，分为第1区段;4，5月和10月降雨量中等，为第2区段;6～9月降雨量较多，为第3区段。根据数据特点，各区段正常降雨强度阈值设定为：≤15mm，≤25mm和≤35mm。由于相邻区域降雨存在相似性，即暴雨出现时间的近同步性，且存在每年汛期偏枯偏涝等不同的特点，故单站奇异点的分析采用同年邻近雨量站降雨对比寻找异常点，而不是该站各年份之间的比较。

2.3 实现方式

根据统计学方法搜寻金沙江下游流域雨量站逐时降雨量存在的累积量异常和奇异点问题，不仅将人工目估的定性方法转为定量，使得问题寻找有据可循，而且易于编程实现，节省了大量的人力时间。研究选用Matlab进行编程，利用内置的utmzone和mfwdtran函数将各雨量站的地理坐标转换为地图坐标，以便雨量站之间距离的测算，来搜寻邻近雨量站;利用内置的K-medoids和Silhouette函数实现聚类分析，完成累积量异常和奇异点问题年份和问题点的提取。应用格拉布斯准则搜寻异常累积雨量时，考虑到累积雨量值分布密集，显著性水平a取0.05即可，分析奇异点时，由于数据差异较大，离散程度高，显著性水平a取0.005。

3 结果及讨论

3.1 计算结果

利用本文提出的检测累积雨量异常和奇异点的方法，对金沙江下游95个雨量站点数据进行逐一核查，检测到2015年龙山村站、2012年大沙店站、2012年八家村水库站、2011年地索（二）站、2014年封过站、2011年后布列托站、2013年细沙站以及2012年龙街（三）站累积雨量异常。经核查原始数据，其中，大沙店站、八家村水库站、地索（二）站、细沙站和龙街（三）站均是由于该年建设站点投入使用时间较迟，数据系列大多从9月开始记录使得年累积雨量异常;龙山村站和淌塘站均是年内累积量异常偏小。这里仅详细分析龙山村站，在累积雨量表1中，先进行各年份比较，查找到2015年累积雨量（79mm）和2011年累积雨量（1031.5mm）异常，再计算同一年份邻近站点累积雨量，龙山村站2009，2015年的累积雨量异常，根据两者输出相同的异常年份为异常值，判定2015年数据需核查矫正;封过站2014年累积雨量达2555.6mm，逐时降雨数据系统性偏大，金沙江下游流域年雨量一般不超过2000mm，故该站点数据需核查校正。

由于奇异点数量较多，这里仅列出表现异常的第1区段（1月至3月和11月至12月）奇异点数据表和数量最多的第3区段（6～9月）奇异点分布直方图。从表2中可以看出，检测到的奇异点，100%的降雨强度大于15mm（暴雨等级），30.3%的降雨强度大于50mm，18.2%大于90mm。可以认为，这些奇异点均属异常数据，需详细核实。图2为第3区段奇异点降雨量级直方图，从95个雨量站共734个年逐时雨量数据中找出164个异常点，可以清晰地看到奇异点逐时降雨量范围在30～50mm居多，有107个，占比65.2%。实践中可核查该区域降雨强度在该段内是否属于正常现象，若属正常，那么采用该方法时正常降雨阈值设定可放宽到50mm;地索（二）站2011年5月22日11时降雨数据为235.9mm，考虑到直方图间距未加入图中，该数据需校正。

3.2 讨论

（1） 格拉布斯准则。剔除粗大误差的统计学方法有很多。虽然经证明，格拉布斯准则适用于测量次数n大于3小于100的情况，但是异常值判别标准G（a，n）与显著性水平a的选取同样相关。在累积雨量异常判别时，为严格起见，取a=0.05。对于该例而言，敏感程度较高，从表2中可看到寻找到了较多的异常值，而该异常值在纵向对比时又为正常值，所以综合考量，实验判定横向纵向对比均为异常数据的值为可疑值。

在奇异点判别时，显著性水平a=0.005，这主要是考虑了逐时降雨数据离差系数大的特点。但从图2中可看到挑选出的异常值65.2%聚集在30～50mm之间。这种现象可能预示着正常数据被当作可疑数据被挑选出来，表明较小的显著性水平依旧无法正确判别，而必须加入阈值设定环节才能使判别更加高效。

（2） 聚类分析。在逐时雨量异常数据判别时，聚类分析契合了逐时降雨数据量大且类间差别大，类内差距小的特点，可以高效地将较大值分成一类，相当于剪枝的思想，而将实验研究对象快速聚焦到奇异点的分析上。但从表2可以看到，区段划分具有极强的主观性，3月中下旬在第1区段中表现异常的值有一部分小于25mm，若是将该时段划分到第2区段，可能属于正常数据而被淹没。在实践中，还需根据研究流域实际降雨特征来划分区段以及选取阈值。

（3） 问题核实。用年鉴资料来核查数据，以2013年向家坝（专）站为例，本文方法搜寻到奇异点2013年2月25日21：00降雨36mm，当日降雨36mm，而年鉴当日显示未降雨;9月3日13：00降雨75.5mm，当日总降雨99.5mm，而年鉴资料当日仅降雨11mm。可见，为严谨使用数据，前期降雨数据质量分析工作非常重要。

本文仅研究极大值出现异常的情况，而未对降雨量小的值进行研究，一方面是方法限制，无法获得正确的参照值;另一方面是因为仪器测量本身存在一定误差，较小的降雨量产生异常一般在误差允许范围内。

查找到的问题数据除了与年鉴进行比较外，也可与其他可靠的气象产品比较，将错误数据矫正，矫正的方法可直接采用正确数据，无准确参考值的情况下可采用邻近站点降雨插值校正等。

4 结语

（1） 本文详细讨论了金沙江下游流域逐时降雨数据质量分析方法，结合格拉布斯准则对累积雨量异常年份进行确定，利用K-medoids聚类方法挑选可疑极大值数据，确定奇异点。年累积雨量挑选方法科学易行且准确率高，奇异点分析方法由于数据离散程度较大，需结合阈值选取以避免将正常值作为异常值的错误。

（2） 基于特征空间的方法来识别奇异点，关键在于分段，需综合考虑流域降雨特点，将降雨特性相近的时段合并，差异大的时段分开，可提高奇异点提取准确性。

（3） 遥测雨量站逐时降雨数据产生的误差分系统误差和粗大误差两类，即累积雨量异常和奇异点两类问题。本文提出的方法是对初步处理大量该类数据的一次探索，以期为越来越精细化的水文数据质量分析提供一定参考依据。

参考文献：

[1]胡海洪.遥测雨量站在青海省中小河流监测站网中的应用[J].中国高新技术企业，2013（21）：59-61.

[2]吴骞，吴绍春.基于离群分析的水位异常识别研究[J].硅谷，2010（24）：45.

[3]詹艳艳，徐荣聪.时间序列异常模式的K-均距异常因子检测[J].计算机工程与应用，2009，45（9）：141-145.

[4]张睿，周建中，肖舸，等.金沙江下游梯级和三峡梯级水电站群联合调度补偿效益分析[J].电网技术，2013，37（10）：2738-2744.

[5]高俊刚，吴雪，张镱锂，等.基于等级层次分析法的金沙江下游地区生态功能分区[J].生态学报，2016，36（1）：134-137.

[6]世界气象组织（WMO）著，赵珂经等译.水文实践指南[M].北京：水利电力出版社，1987.，

[7]熊艳艳，吴先球.粗大误差四种判别准则的比较和应用[J].大学物理实验，2010，23（1）：66-68.

[8]夏宁霞，苏一丹，覃希.一种高效的K-medoids聚类算法[J].计算机应用研究，2010，27（12）：4517-4519.

[9]周世兵，徐振源，唐旭清.新的K-均值算法最佳聚类数确定方法[J].计算机工程与应用，2010，46（16）：27-31.

[10]杨茂林，卢炎生.基于剪枝的海量数据离群点挖掘[J].计算机科学，2012，39（10）：152-156.

[1]RoySS，RouaultM.Spatialpattermsofseasonalscaletrendsinex-

tremehourlyprecipitationinSouthAfrica[J].AppliedGeography，2013（39）：151-157.

[12]Gwo-FongLin，Ming-JuiChang，Chian-FuWang.ANovel

SpatiotemporalStatisticalDownscalingMethodforHourlyRainfall[J].WaterResourManage，2017（31）：3465-3489.

[13]周國良，高唯清，黄昌兴.2016年我国极端暴雨事件浅析[J].中国防汛抗旱，2017，27（1）：75-78，87.

[14]赵超，包为民，瞿思敏，等.遥测系统降雨观测粗差修正研究[J].人民长江，2003，34（2）：4-5，55.

[15]李朋军.遥测与虹吸雨量计降水数据对比分析[J].水科学与工程技术，2016（2）：51-53.

[16]楊旭，刘志武，李波.多源降水数据在长江上游流城比较研究[J].长江流域资源与环境，2016，25（1）：131-139.

引用本文：刘秀林，张行南，方园皓，黄晴.金沙江下游遥测雨量站数据质量研究[J].人民长江，2019，50（3）：131-135.

Study on data quality of hourly rainfall of telemetry rainfall stations in lower reaches of Jinsha River

LIU Xiulin'，ZHANG Xingnan' 1，2，3，FANG Yuanhao'，HUANG Qing*

（1. College of Hydrology and Water Resources，Hohai University，Nanjing 210098，China;2. National Cooperative Innovation Center for Water Safety & Hydro-Science，Hohai University，Nanjing 210098，China;3. National Engineering Research Center of Water Resources Efficient Utilization and Engineering Safety，Hohai University，Nanjing 210098，China;4. Geographic and Oceano-graphic Sciences College，Nanjing University，Nanjing 210098，China）

Abstract：Along with the development of automatic measurement and forecast technology，more and more automatic rainfallstations are put into operation. Analyzing the quality of data measured by automatic rainfall station is of great significance to the operation of water project. The analysis on the hourly rainfall series from telemetry rainfall stations in the lower reaches of the Jinsha River from 2008 to 2015 was carried out by using the Grubbs criterion and k-medoids clustering algorithm. Based on the statistical theory，the results showed that the methods have a remarkable effect in determining the abnormal value of annual accumulated rainfall and the singular points of hourly rainfall data. The method provides some references for hydrological data qualityanalysis that is changed by increasingly fine particles.

Key words：telemetry rainfall data;hourly rainfall series;Grubbs criterion;k-medoids clustering;singular points;lower reaches of Jinsha River