不同水文序列突变检测方法在漳河观台站径流分析中的对比研究

鲍振鑫，张建云，王国庆，刘翠善，严小林，刘 晶，刘 悦

(1.南京水利科学研究院 水文水资源与水利工程科学国家重点实验室，南京 210029； 2.水利部应对气候变化研究中心，南京 210029；3.河海大学，南京 210098)

0 引 言

随机性是水文序列的重要特性之一。受多种因素的影响和制约，长时间尺度水文序列的演变规律表现出显著的随机性，很难用确定的数学物理方程来描述[1]。统计方法是研究水文序列随机性的重要工具之一。利用统计方法分析水文序列的均值、方差等统计特征，有一个基本前提，即统计样本服从独立同分布[2]。但是近些年来由于环境变化的影响，水文序列的一致性受到破坏，对传统统计方法研究随机水文特征带来了挑战[3]。因此，识别水文序列的突变特征，寻找可用于统计分析的水文序列段，对于科学认识水文要素的演变规律，分析环境变化对水文过程的影响，具有重要的科学意义和实用价值。

突变检测方法包括参数检验法和非参数检验法等，被广泛应用于相关流域水文序列的突变点分析，取得了大量的研究成果[4-6]。例如，田小靖等[7]利用有序聚类分析法、Mann-Kendall检验法、Pettitt法等9种突变检测方法，研究了黄河中游头道拐站和龙门站的输沙变化。袁满[8]等同时考虑类内部的离差较小和类间的离差较大原则，提出了改进的有序聚类分析法，识别了年平均流量序列的突变点，结果更合理。梁欣阳等[9]利用秩和检验法和TFPW-MK检验法分析了黄河上游兰州站的水文序列突变特征，结果表明1985年前后水文指标发生了显著变化。黄丽娜[10]利用空间差异突变诊断模型研究了淮河上游降水和径流的突变特征，结果表明降水和径流突变发生在20世纪80年代中期。张敬平[11]等利用漳泽水库的降水和径流序列，对比分析了有序聚类法和启发式分割算法在水文序列突变检测中的应用，结果表明启发式分割法能更好地适用于非线性、非平稳水文序列的变异诊断。

以往主要侧重于利用突变检测方法分析水文序列的突变点，关于多种突变检测方法有何异同，特别是不同方法的适用范围，及其检测结果的物理意义有何异同等方面的研究较少。本文主要对比分析不同水文序列突变检测方法，研究其适用范围及突变点的物理意义，为水文序列演变规律分析提供技术支撑。

1 研究区与资料

以海河流域漳河水系观台水文站以上流域为研究对象，流域概况见图1。流域地处中国北方，位于北纬35.9°～37.6°，东经112.4°～114°，地跨山西、河北和河南三省，流域面积17 800 km2。流域以山区地貌为主，处于半湿润区，属于温带大陆性季风气候，四季分明，多年平均气温13.2 ℃，降水560 mm，径流深53 mm，降水-径流系数小于0.1。流域内支流众多，水系呈扇形分布，有北部的清漳河与南部的浊漳河两大水系。流域内有漳泽等3座大型水库，12座中型水库和93座小型水库，总库容约14 亿m3，并建有4个大型跨流域引水工程[12]。利用观台水文站观测的1951-2015年径流序列，分析漳河流域水文特征的突变性。年径流数据来源于水利部刊印的中国水文年鉴，该资料经过水文资料整编等多重手续，结果可靠。

图1 观台水文站以上漳河流域概况Fig.1 The basic information of the Zhanghe River Basin above the Guantai hydrologic station

2 研究方法

水文序列的突变检测方法种类繁多，包括经验方法、统计学方法和水文模拟法等，其中统计学方法种类最多，应用较广泛[13]。本文选用常用的经验方法中的降水-径流双累积曲线，统计学方法中的有序聚类突变检测方法和Mann-Kendall突变检测方法，以及水文模拟法中的基于水文模拟的突变检测方法等多种方法对比分析水文序列突变检测结果的异同。

2.1 降水-径流双累积曲线

降水-径流双累积曲线是一种广泛用于分析流域降水-径流特性是否发生变化的经验方法[14]，其横坐标为累积年降水量(PSt)，纵坐标为累积年径流量(RSt)，计算公式如下：

(1)

式中：P和R分别为年降水量和年径流量；n为序列长度。

通过查看降水-径流双累积曲线是否发生转折来判断序列的突变点。一般为了排除偶然因素的干扰，只有在双累积曲线的坡度变化显著，且转折后有连续5年以上的观测资料时，才判定降水-径流序列发生了突变。该方法主要依赖于专家经验判断，未通过统计变量检测。

2.2 有序聚类突变检测方法

有序聚类突变检测方法最早由丁晶等人[15]提出，广泛应用于水文序列的突变分析。该方法的基本思想是，有一个样本长度为n的水文序列，首先任意假定某一点(τ)为突变点，将水文序列分成两段，则这两段水文序列分别代表两类水文样本，统计两段样本的离差平方和(Snτ)，其计算公式为：

(2)

根据同类水文样本间的方差较小，类之间的方差较大这一原则，找到Snτ的最小值，则相应的为水文样本的突变点。

2.3 Mann-Kendall突变检测方法

Mann-Kendall检测方法是一种非参数检验方法，最早由H B Mann和M G Kendall提出[16]。首先根据样本序列计算xi>xj(i>j)的个数(Sk)，其计算公式为：

(3)

(4)

式中：xi是样本长度为n的水文序列。

再统计Si的期望值[E(Si)]和方差[Var(Si)]：

(5)

则Mann-Kendall统计量UFi可由下式估计：

(6)

对于某一个显著性水平α，查得正态分布临界值Uα/2。如果|UFi|

2.4 基于水文模拟的突变检测方法

考虑到传统基于统计方法的突变检测技术只针对单一的水文序列如何改变，即将人类活动和气候变化对水文过程的影响混合在一起，为了更科学地分析人类活动改变降水-径流特性，进而引起水文过程改变的机理，Wang等[17]提出了基于水文模拟的突变检测技术。首先根据研究流域的实际情况，选择未受人类活动影响的天然时期，利用天然时期的水文数据率定水文模型参数，满足天然时期的观测值和模拟值没有系统偏差，即模拟值围绕观测值波动。天然时期率定的模型参数反映流域在天然时期未受人类活动干扰的下垫面特性。利用天然时期率定的模型参数，模拟整个时期的径流过程。如果模拟径流与实测径流出现系统偏差，则将系统偏差的起始点认为是人类活动对流域水文过程干扰的起始时刻，即流域降水-径流特性发生改变的突变点。定义模拟径流和实测径流的偏差(Ki)与累积偏差(Km)如下：

(7)

式中：QSi和QOi分别为模拟径流和实测径流。

在天然时期，累积偏差(Km)在0值附近波动，将其持续大于0或者小于0的起始点判定为流域降水-径流特性发生改变的突变点。

目前水文模型有很多种，可以根据流域特征和径流模拟的需要来选择合适的水文模型。本文选用VIC(Variable Infiltration Capacity)模型来模拟观台水文站的径流过程。VIC模型是一个基于正交网格的分布式水文模型，它能够同时模拟地表间的能量平衡和水量平衡，广泛应用于全球很多流域的径流模拟[18, 19]。VIC模型将土层划分为3层，总径流由上两层土壤产生的地表径流和第3层土壤产生的基流两部分组成。借用了新安江模型中的蓄水容量曲线来描述土壤含水量的空间分布不均匀性，计算地表径流过程[20]。基流部分采用Arno模型中的原理来计算，即当土壤含水量在某一阈值以下时，基流是线性消退的；而高于此阈值时，基流过程是非线性的[21]。在模型计算时，首先分别在每个网格内独立进行降水-蒸发-产流过程的计算，然后再统一汇流到流域出口断面形成流量过程。VIC模型中共有6个参数需要利用观测的流量资料来进行率定。包括第2层和第3层的土层厚度，描述土壤蓄水容量曲线的形状参数，以及控制基流的3个参数。模型参数的率定主要考虑2个目标函数：Nash-Sutcliffe效率系数和相对误差。其中，Nash-Sutcliffe效率系数反映了模拟的流量过程和观测的流量过程之间的吻合程度，其值越接近1，表示模拟效果越好；而相对误差是一个水量平衡指标，它反映了模拟总径流量和观测总径流量之间的相对误差，其值越接近0，则表示模拟效果越好。

3 结 果

3.1 观台水文站年径流序列变化特征

1951-2015年观台水文站观测的年径流序列见图2。从图2可以看出，1951年以来，观台水文站观测的年径流序列呈现出持续性的下降趋势，其中最大值发生在1963年，为259 mm，最小值发生在1999年，仅为1.67 mm，最大值是最小值的155倍。20世纪50和60年代处于丰水期，年径流在100 mm左右波动；从70年代初开始，年径流持续下降，到70年代末80年代初，年径流降低到10 mm以下；随后该流域进入连续枯水期，年径流在16 mm左右波动。从观台水文站实测年径流序列的演变特征来看，径流序列的突变点发生在20世纪70年代。

图2 观台水文站年径流序列Fig.2 The observed annual streamflow at the Guantai hydrologic station

3.2 观台水文站年径流序列突变检测分析

观台水文站以上流域年降水量和径流量双累积曲线见图3。从图3可以看出，降水-径流关系的突变点发生在1978年。在同样的降水条件下，1978年以后的年径流量明显小于1977年之前的径流量。例如1951-1977年，年降水径流系数为0.158，但是1978-2015年的降水径流系数仅为0.034。除了1978年这一明显的转折点以外，在1971年降水量和径流量双累积曲线开始有所转折，但是不是很显著。这表明降水-径流关系从1971年开始有所变化，但是变化程度不大；从1978年开始降水-径流关系发生了突变。

图3 观台水文站年降水量和径流量双累积曲线Fig.3 The double-accumulation curve of precipitation and runoff at the Guantai hydrologic station

观台水文站实测年径流序列的有序聚类突变检测结果见图4。从图4可以看出，年径流序列离差平方和的最小值位于1977年，则可以将整个序列分成1951-1977和1978-2015两段。这一结果和降水-径流双累积曲线的分析结果相一致。离差平方和的次小值是1973年，其值与1977年的值相差很小，而且从20世纪60年代后期到70年代末，年径流序列离差平方和的值都比较小，在谷值附近波动。这表明年径流序列的突变不一定是在某一个年份突然发生的，而是在70年代缓慢进行的。这一结论和年径流的演变特征相吻合，即在20世纪70年代，径流持续下降，而不是在某一年突然跃变。

图4 观台水文站年径流序列有序聚类突变检测Fig.4 The breakpoint of annual streamflow at the Guantai hydrologic station detected by the orderly clustering breakpoint detection method

观台水文站实测年径流序列的Mann-Kendall突变检测结果见图5。从图5可以看出，Mann-Kendall统计量为-6，表明年径流呈现出显著的下降趋势，达到了1%的显著性水平。UF和UB曲线的交点位于1973和1974年之间，表明年径流序列的突变点发生在1974年，则可以将整个序列分成1951-1973和1974-2015两段。这一结果略早于前两种方法检测的突变点，但是仍处于年径流持续下降的20世纪70年代之中。

图5 -观台水文站年径流序列Mann-Kendall突变检测Fig.5 The breakpoint of annual streamflow at the Guantai hydrologic station detected by the Mann-Kendall’s test methodology

观台水文站实测与VIC模型模拟的1955-1970年月径流过程见图6。从图6可见VIC模型具有较好的径流模型效果，Nash-Sutcliffe效率系数超过了0.7，相对误差小于5%，可用于观台站的径流模拟。观台水文站实测年径流序列基于水文模拟的突变检测结果见图7。从图7可以看出，模拟年径流和实测年径流的累积差值在20世纪50年代和60年代围绕0值附近波动，在1971年以后持续偏大。这表明从70年代初开始，由于人类活动的影响，实测径流小于模拟径流，即降水-径流关系开始发生改变。这与实测径流下降的起始时刻相吻合。

图6 观台水文站实测与模拟流量过程Fig.6 The observed and simulated streamflow at the Guantai hydrologic station

图7 观台水文站年径流序列基于水文模拟的突变检测Fig.7 The breakpoint of annual streamflow at the Guantai hydrologi c station detected by the hydrological simulation based method

3.3 多种突变检测方法对比分析

根据上文分析结果，利用降水-径流双累积曲线、有序聚类突变检测方法、Mann-Kendall突变检测方法、基于水文模拟的突变检测方法等多种方法，分析得到观台水文站年径流序列的突变点分别是1978、1978、1974和1971年。可见基于水文模拟的突变检测方法分析得到的突变点最早，Mann-Kendall突变检测方法结果其次，而降水-径流双累积曲线和有序聚类突变检测方法的结果最晚。

气候要素的分析结果表明，漳河流域降水在20世纪50和60年代较大，70-90年代持续下降，但是在21世纪有所回升。从60年代初开始，流域内修建了大量的水库和引调水工程。1978年中国开始实行改革开放，经济快速发展，对水资源的需求迅速增加，极大地改变了流域的降水-径流特性。结合径流的演变特征和流域内的气候变化与人类活动情势，4种突变检测方法检测的径流突变点物理意义有所差别。降水-径流双累积曲线检测的突变点反映的是流域降水-径流关系发生改变的时刻，在本案例中从1971年开始流域降水径流关系开始发生变化但是程度不大，从1978年开始发生了显著变化。有序聚类突变检测方法检测的突变点反映的是径流序列的聚类特征，在本案例中以1978年为界分成前后两类，但是从20世纪60年代后期到70年代末的离差平方和都较小，表明以其中任意一年为突变点都具有一定的意义。Mann-Kendall突变检测方法检测的突变点反映的是径流序列的突变特征，在本案例中其分析结果和有序聚类突变检测方法的次小值相一致。基于水文模拟的突变检测方法检测的突变点反映的是人类活动对流域降水径流关系干扰的起始点，即本案例中的1971年，因此其结果早于其他方法的检验结果。综合分析四种突变检测方法计算结果的物理意义和流域实际气候、水文、人类活动的状况相符合。

4 主要结论

科学诊断水文序列的突变特性，对于认识水文循环的演变规律，研究环境变化对水文过程的影响具有重要的科学意义和实用价值。基于漳河观台水文站观测的1951-2015年径流序列，利用降水-径流双累积曲线、有序聚类突变检测方法、Mann-Kendall突变检测方法、基于水文模拟的突变检测方法等多种方法分析了其突变特征。对比分析了各种突变检测方法的异同及适用范围。主要结论如下。

(1)1951年以来，观台水文站实测年径流序列呈现出持续性的下降趋势，20世纪50和60年代处于丰水期，从70年代初开始年径流持续下降，80年代进入连续枯水期。

(2)降水-径流双累积曲线、有序聚类突变检测方法、Mann-Kendall突变检测方法和基于水文模拟的突变检测方法等多种方法，检测得到观台水文站年径流序列的突变点分别是1978、1978、1974和1971年。

(3)相比而言，降水-径流双累积曲线检测的突变点反映的是流域降水-径流关系发生明显改变的时刻；有序聚类突变检测方法检测的突变点反映的是径流序列的聚类特征；Mann-Kendall突变检测方法检测的突变点反映的是径流序列的综合突变特征；基于水文模拟的突变检测方法检测的突变点反映的是人类活动对流域降水径流关系干扰的起始点。

在分析流域水文序列的突变特性时，应根据研究的需要合理选用相应的突变检测方法，科学识别水文序列的变异情势。

□