汾河上游径流序列突变点及前兆信号检测

秦超瑞，桑宇婷，赵雪花，祝雪萍

（太原理工大学 水利科学与工程学院，山西 太原 030024）

受自然因素及人类活动的影响，河流年径流量会发生显著变化［1-3］，径流序列发生突变［4］，导致旱涝灾害等一系列水文极值事件频发，研究径流序列的突变规律可以为流域治理与水资源合理开发利用提供科学依据。 目前，水文气象方面常用的突变检测方法主要有参数统计法和非参数秩检测统计法两大类。 其中参数统计法，如累积距平法［5］、滑动t 检测法［6-7］及线性回归法等，适用于满足一定分布的数据序列，但在实际应用中无法对样本序列的分布形态做出简单假设，且统计值会受到原始数据序列长度、数据异常及缺失的影响；非参数秩检测统计法，如Mann-Kendall 检测法［8］和Pettitt 突变检测法［9］等，无法去除样本数据所固有的相关属性，从而存在分析误差。 多尺度直线拟合法可避免原始数据属性及其固有相关属性的影响，很好地克服上述方法的缺点，已被应用于电液伺服阀的故障诊断［10］，并取得了很好的效果。 除径流突变规律外，突变前兆信号也很重要。 临界慢化现象由于渌等［11］于1984年提出，具有揭示复杂动力系统是否趋于临界灾变的优势。 2009年，Scheffer 等［12］研究表明，当系统趋近于分岔点时，临界慢化现象会导致动力学出现3 个可能的早期预警信号（扰动的恢复变慢、方差增大、自相关系数增大）。 由此为突变前兆信号研究提供了一种新途径，并被广泛用于多个领域，如晏锐等［13］将临界慢化现象用于地震前兆信号检测；颜鹏程等［14］和吴浩等［15］将临界慢化现象应用于气象突变前兆信号研究中，证明临界慢化现象是气候突变前的一个早期信号。 但临界慢化现象在径流序列中的应用仍较少，因此本文采用多尺度直线拟合法分析汾河上游水文站的径流序列突变点并根据系统恢复速率及恢复力、方差及自相关系数的演变规律，探讨前兆信号发生年份作为早期预警的可能性。

1 研究内容与研究方法

1.1 研究区域概况及数据来源

汾河位于山西省境内，是山西省第一大河流，流域面积为39471 km2。 汾河上游段从河源宁武县东寨镇管涔山到太原兰村，流域面积为7705 km2，是太原市主要的水源涵养区和水土保持生态功能区。 上游流域位于黄土高原地区，属于温带大陆性季风区，春冬干旱少雨，夏秋多雨，且大雨、暴雨多集中于7月、8月。 研究区地处山区，无大型城市群，城镇化程度较低，工农业处于落后水平，但矿产资源非常丰富。 20 世纪60年代开始，人类活动对径流影响加剧，主要体现在水利工程建设及水资源开采。 流域内有汾河水库和汾河二库两座大型水库，其中汾河水库是山西省最大的水库，位于汾河干流上游，地处娄烦县杜交曲镇，设计总库容为7.21 亿m3，是一座集防洪、灌溉、发电为一体的水利工程，于1956年拦洪，1961年竣工验收后正式投入运行。 本文以汾河上游4 个水文站（上静游站、汾河水库站、寨上站及兰村站）1956—2000年的年径流量序列（资料来自山西省水文局）为研究资料。 4 个水文站年降水量最大值均发生于1967年，最小值均发生于1972年。 其中上静游站位于汾河支流岚河上，控制流域面积为1140 km2。 汾河水库站、寨上站及兰村站位于汾河干流上，控制流域面积分别为5268、6819、7705 km2。 4 个水文站在汾河上游流域内的地理位置如图1所示。

图1 汾河上游流域水系及水文站分布

1.2 研究内容

针对汾河上游流域4 个水文站的年径流量序列进行突变点及突变前兆信号检测。 具体分为如下3 个部分：采用多尺度直线拟合法检测年径流量序列突变点，根据系统的恢复速率及恢复力演变规律检测突变前兆信号发生年份，根据序列的方差及自相关系数演变规律检测突变前兆信号发生年份。

1.3 研究方法

1.3.1 多尺度直线拟合法

根据多尺度直线拟合法原理［10］，选取合适的尺度，用多条线段代替原序列，可利用多条拟合线段提取原序列的主要信息，以达到快速有效地反映原序列信息的目的。 初始尺度r1计算公式如下：

式中：n为径流序列长度。

采用式（1）计算得到初始尺度r1后，可将原序列划分成m段，每段用di（i取1 ～m）来表示，如最后一段长度小于r1，则将其并入前一段即可。 对于每段序列，均可用最小二乘法进行直线拟合。 根据最小二乘法原理可求得每段拟合线段的斜率ki，并依次求出相邻两段线段斜率差值的绝对值，找出最大值记为emax，说明原序列在该相邻两段线段间的变化最大，则突变点一定落在该区间内，此时将r1变为r2，r2=E（r1／2）。在r2尺度下对新区间采用上述方法逼近突变点，直至尺度rn=2 时，此时的区间长度为3，区间的中点即为原径流序列的突变点。

1.3.2 基于系统的恢复速率及恢复力演变规律检测突变前兆信号

径流系统是一个复杂的动力学系统，径流量在变化环境中演变时往往会发生突变，而径流系统在发生突变前总是先演变到某个临界阈值附近，一段时间后系统才发生突变。 因此，径流系统在临界阈值附近的一些特征可作为判断径流系统临近突变的前兆信号。由临界慢化现象［11］（动力系统由一种相态转变为另一种相态之前，系统在趋近于临界点附近，表现出涨落幅度增大、扰动的恢复速率变慢以及恢复到旧相位能力变小的现象）可知，系统的恢复速率及恢复力可作为表征径流系统突变前兆信号的特征量。 构建径流序列混沌模型，假设径流量变化是从一个平衡态经历一次突变后向另一个平衡态转移的过程，系统的恢复速率及恢复力是指径流序列在受到突变干扰后重新恢复到平衡态的速率和能力。 系统的恢复速率可表示为

式中：x为径流状态变量；λ为径流量的增长系数；λx为径流量的增长速率；k为径流量的减少系数；kx2为径流量的减少速率；t为时间。

式（2）的物理意义是径流状态变量x随着时间t增大（减小）的速率。 该模型中的初始值和参数λ、k只能取正值，为了便于讨论，令λ／k=μ，由式（2）可知，当dx／dt=0 时，x=0 和x=μ为式（2）的两个解，此时系统的状态变量不变，即径流量达到平衡状态。 式（2）等号两边同时对t求导，可以得到系统的恢复力：

由式（3）可知，当x=0、x=μ／2 和x=μ时，系统的恢复力为零。

该模型描述的是径流状态变量由初始平衡态（x=x0）经历一次突变后趋向于另一个平衡态（x=μ）。 当径流状态变量开始发生演变时，系统的恢复速率和恢复力也发生演变，且变化幅度大于径流状态变量的变化幅度。 根据临界慢化现象，在径流状态变量趋于突变点之前会出现系统的恢复速率及恢复力变小的现象，表现为存在某一时刻系统的恢复速率及恢复力同时为零，并发生相位转变，该时刻能很好地反映径流状态变量偏离平衡态的演变过程，可作为径流系统发生突变的前兆信号。

1.3.3 基于方差及自相关系数演变规律检测突变前兆信号

临界慢化现象除了导致系统的恢复速率及恢复力出现上述演变规律外，还将导致序列的方差和自相关系数增大［11］。 因此，径流系统趋近阈值附近时序列的方差和自相关系数同时达到极小值并开始逐渐增大，表明径流系统即将发生突变，也可作为径流系统发生突变的前兆信号。

方差是概率论中用来度量随机变量与其均值之间偏离程度的指标。 自相关系数度量的是同一变量在两个不同时期的相关程度，即自身过去行为对现在的影响，自相关系数c的定义如下：

式中：j为滞后长度；xi为序列中第i个样本值； ˉx为序列均值；s为序列均方差。

2 结果与分析

2.1 年径流量序列突变点检测

将多尺度直线拟合法应用于汾河上游4 个水文站年径流量序列突变点检测，结果如图2 所示。

图2 年径流量序列突变点检测结果

图2（a）为上静游站年径流量序列的突变点检测结果，经计算可得分割点3、4 之间线段的斜率和分割点4、5 之间线段的斜率差值绝对值最大，由多尺度直线拟合法原理可知突变点一定落在分割点3 与分割点5 之间，因此将寻求突变点的范围缩小至分割点3 与分割点5 之间。 采用同上的方法不断缩小范围和尺度，直至最后尺度rn=2 时，此时区间的中点即为原始序列的突变点。 图中A点为最终检测到突变点的位置，由此可知上静游站年径流量序列突变年份为1975年，年径流量均值在1975年前后分别为0.57 亿m3和0.44 亿m3，突变后较突变前减少22.81%。 图2（b）～（d）为汾河水库站、寨上站和兰村站年径流量序列突变点检测结果，经计算可得3 个水文站初始尺度下拟合线段斜率差值绝对值最大值均为分割点2、3 之间的线段与分割点3、4 之间的线段之差，因此将寻求突变点的范围缩小至分割点2 与分割点4 之间，采用同上的方法逐渐逼近突变点，最终A点为检测到的突变点位置，由此可知汾河水库站、寨上站和兰村站年径流量序列的突变年份均为1966年。 其中汾河水库站年径流量均值在1966年前后分别为4.00 亿m3和3.21亿m3，突变后较突变前减少19.75%；寨上站年径流量均值在1966年前后分别为4.74 亿m3和3.57 亿m3，突变后较突变前减少24.68%；兰村站年径流量均值在1966年前后分别为4.72 亿m3和3.25 亿m3，突变后较突变前减少31.14%。 干流上3 个水文站径流量的减少程度自上而下依次递增。 通过降水径流双累积曲线斜率变化率计算［16］，可得降水和人类活动分别对4 个水文站径流量序列突变的贡献率，结果见表1。

表1 降水和人类活动对突变的贡献率

由表1可知，汾河干流上的水文站年径流量序列突变一致性较好，且人类活动是引起干流径流量序列突变的主要因素。 而支流的年径流量序列突变一致性较差，降水是引起上静游站径流量序列突变的主要因素。 降水是径流产生的源泉，因此降水变化是影响径流变化最敏感的因素。 4 个水文站年降水量均于1967年达到最大值，于1972年减到最小值。 由于短时期内降水量发生较大幅度的波动且支流控制的流域面积远小于干流，因此上静游站1975年发生突变，而干流沿程不断有支流汇入导致干流站控制的流域面积不断增大以缓冲降水波动带来的影响，汾河水库站、寨上站和兰村站受汾河水库的调节作用使得这次降水量大幅度波动未对干流上的水文站径流产生突变影响。

人类活动对支流的上静游站影响较小，对干流径流变化的影响主要体现在水利工程和水资源开采两个方面。 水利工程主要为汾河水库的修建，汾河水库投入使用后防洪效益显著，1965年全年来水量为2.01 亿m3，而水库给下游放水6.29 亿m3，充分发挥了水库的调节作用，保证了下游的用水需求［16］。 工业用水量和农业用水量急剧增加，工业用水量在1949年约为180万m3，从汾河水库1960年投入运行后增至3051.11万m3，1956年汾河灌区面积为8.32 万hm2，在汾河水库1960年投入运行后灌溉面积增至8.46 万hm2［17］。因此，汾河水库下游工农业用水量急剧增加是干流水文站径流量序列发生突变的原因之一。 汾河上游水土流失严重，汾河水库的修建使上游来沙大量被拦蓄，下游输沙量急剧减少（1965年水库泥沙淤积为0.0125 亿m3，1966年为0.1270 亿m3）［18］，因此汾河水库拦蓄泥沙也是造成干流水文站径流量序列发生突变的原因之一。 随着城镇化程度提高，汾河上游流域内地下水过量开采（1949年以前每日不超过1 万m3，到1965年增至每日32.93 万m3），导致地下水位下降，深层地下水漏斗从1965年起迅速扩大［19］，而取水口一般设在支流和干流汇合后的河段上，干流上3 个水文站均处于支流和干流汇合后的河段上，因此人类对地下水资源的过度开采也是导致干流径流量序列发生突变的原因之一。 干流上3 个水文站径流量减少程度自上而下依次递增的原因有两个：一是汾河水库修建于汾河干流上游，自上而下对汾河水库站、寨上站和兰村站的调节作用递减；二是沿着干流不断有支流汇入，导致径流量减少程度增大。

综上，位于支流的上静游站径流量序列1975年发生突变的主要影响因素为降水，而位于干流的汾河水库站、寨上站和兰村站径流量序列1966年发生突变的主要影响因素为人类活动。

2.2 基于系统的恢复速率及恢复力演变规律检测突变前兆信号结果分析

为了更好地检测径流突变年份，将系统恢复速率及恢复力应用于汾河上游水文站径流量序列突变前兆信号检测中，如图3 所示。 由图3（a）和（b）可知上静游站年径流量序列的系统恢复速率在1973年为零并由正相位向负相位转变（见图3（a）中A点），恢复力在1973年为零并由正相位转变为负相位（见图3（b）中B点），恢复速率和恢复力同时为零且发生相位转变的点只有1973年，则上静游站年径流量序列的突变前兆信号发生于1973年。 由图3（c）～（h）可知，汾河水库站、寨上站和兰村站年径流量序列的系统恢复速率均在1958年附近为零并由负相位转变为正相位（见图3（c）（e）（g）中A点），恢复力均在1958年附近为零并由正相位转变为负相位（见图3（d）（f）（h）中B点），恢复速率和恢复力同时为零且发生相位转变的点只有1958年，则干流水文站的突变前兆信号均发生于1958年。 综上，根据径流量序列的系统恢复速率及恢复力演变规律，得到上静游站的突变前兆信号发生于1973年，汾河水库站、寨上站和兰村站的突变前兆信号均发生于1958年。

图3 年径流量序列恢复速率及恢复力演变过程

在适当的时间范围内前兆信号发生年份较突变年份越早，系统的恢复速率及恢复力对引起突变的因子越敏感，可为预防突变发生提供充足的准备时间，因此预警效果越优。 根据径流状态变量恢复速率及恢复力的演变规律，得到上静游站年径流量序列的突变前兆信号发生于1973年，比实际突变年份1975年提前2 a。 汾河水库站、寨上站和兰村站年径流量序列的突变前兆信号均发生在1958年，比实际突变年份1966年提前8 a。 突变前兆信号均早于突变年份10 a 之内，说明径流状态变量的恢复速率及恢复力的演变规律能为径流量序列提供可靠的突变前兆信号。 由于前兆信号发生年份在支流上只比突变年份提前2 a，和突变年份较为接近，因此预警效果不如干流。 前兆信号发生年份在干流上比突变年份提前8 a，因此针对汾河上游流域，基于系统恢复速率及恢复力的演变规律更适合求取干流径流量序列的突变前兆信号。 由于干流径流量序列的突变受人类活动的影响较大，受降水影响较小，因此径流状态变量的恢复速率及恢复力对人类活动干扰更敏感。

2.3 基于方差及自相关系数演变规律检测突变前兆信号结果分析

基于临界慢化现象，将方差和自相关系数应用于汾河上游水文站年径流量序列突变前兆信号检测中，结果如图4 所示。 由图4（a）和（b）可知，上静游站径流量序列的方差和自相关系数同时达到极小值的点只有1968年（见图4（a）和（b）中的A点和B点），之后开始逐渐增大，则上静游站径流量序列的突变前兆信号发生于1968年。 由图4（c）～（h）可知，汾河水库站、寨上站和兰村站的方差和自相关系数同时达到极小值的点只有1960年（见图4（c）～（h）中A点和B点），之后开始逐渐增大，则3 个水文站径流量序列的突变前兆信号均发生于1960年。 综上，根据径流量序列的方差和自相关系数演变规律得到上静游站径流量序列的突变前兆信号发生于1968年，干流3 个水文站径流量序列的突变前兆信号均发生于1960年。

图4 年径流量序列方差及自相关系数演变过程

根据临界状态时径流量序列的方差及自相关系数的演变规律，得到上静游站径流量序列突变的前兆信号发生于1968年，比实际突变年份1975年提前7 a。汾河水库站、寨上站和兰村站3 个水文站径流量序列突变的前兆信号均发生于1960年，比实际突变年份1966年提前6 a。 突变前兆信号均早于突变年份10 a之内，说明径流量序列的方差及自相关系数的演变规律能为径流量序列提供可靠的突变前兆信号。 前兆信号发生年份在支流上较突变年份出现更早，预警效果优于干流。 因此，基于方差及自相关系数的演变规律检测突变前兆信号更适用于支流径流量序列。 由于支流径流量序列的突变受降水影响较大，受人类活动的影响较少，因此径流量序列的方差和自相关系数对降水量变化更敏感。

3 讨论与结论

3.1 讨论

采用多尺度直线拟合法对汾河上游径流量序列进行突变点检测，计算降水和人类活动分别对突变的贡献率并分析其原因。 根据临界慢化现象检测突变前兆信号发生年份，探讨了前兆信号作为早期预警的可能性。 需要指出的是多尺度直线拟合法只能检测出序列中突变程度最显著的一个突变点，在今后的研究工作中，需改进检测方法，使其可以检测出序列中所有的突变点。 同时，对临界慢化现象是否受突变类型及突变强度影响仍然需要进一步研究。

3.2 结论

（1）位于支流的上静游站年径流量序列的突变年份为1975年，年径流量均值在1975年前后分别为0.57 万m3和0. 44 万m3， 突 变 后 较 突 变 前 减 少22.81%，降水对突变的贡献率为79.67%，是突变的主要影响因素。 位于干流的汾河水库站、寨上站和兰村站年径流量序列的突变年份均为1966年，其中汾河水库站年径流量均值在1966年前后分别为4.00 万m3和3.21 万m3，突变后较突变前减少19.75%，人类活动对突变的贡献率为62.46%；寨上站年径流量均值在1966年前后分别为4.74 万m3和3.57 万m3，突变后较突变前减少24.68%，人类活动对突变的贡献率为64.97%；兰村站年径流量均值在1966年前后分别为4.72 万m3和3.25 万m3，突变后较突变前减少31.14%，人类活动对突变的贡献率为68.94%，干流径流量序列突变的主要影响因素为人类活动。

（2）根据临界慢化现象，径流量序列在发生突变前径流状态变量的恢复速率变慢和恢复力变小，得到上静游站的突变前兆信号发生于1973年，汾河水库站、寨上站和兰村站的突变前兆信号均发生于1958年，且径流状态变量的恢复速率及恢复力对人类活动干扰更敏感，更适用于求取汾河上游干流的突变前兆年份。

（3）根据临界慢化现象，径流量序列在发生突变前将导致序列的方差和自相关系数增大，得到上静游站年径流量序列的突变前兆信号发生于1968年，汾河水库站、寨上站和兰村站的年径流量序列的突变前兆信号均发生于1960年，且径流量序列的方差和自相关系数对降水量变化更敏感，更适用于求取汾河上游支流的突变前兆年份。