梁欣阳，卢玉东，孙东永，张洪波，杨亚慧(1.长安大学环境科学与工程学院，西安 710054;2.长安大学 旱区地下水文与生态效应教育部重点实验室，西安 710054)

0 引 言

近些年，在全球气候变化与人类活动共同作用下，河川径流的演变规律受到各国学者的关注[1]。研究表明，变化环境下世界上很多河流的径流都表现出不同程度的下降趋势[2]，这一现象在我国也不例外，黄河上游表现尤为明显。近年来黄河上游径流量总体呈下降趋势，到1990年后平均水量只有270～290 亿m3，平均较常年偏少15.6%[3]。在黄河上游径流量的下降趋势过程中，针对水文变异情况作进一步分析对于具有操作性的河流生态水文调度具有深远的意义。如樊辉等[4]采用TFPW-MK法检验年径流与年输沙长期变化趋势，并利用贝叶斯变点分析方法诊断时序数据突变点，分析了黄河干、支流径流量年际变化特征；李剑锋等[5]在考虑水文变异的基础上，系统提出河道内生态需水的计算方法，并对黄河流域生态需水问题进行分析；于茜等[6]应用变化范围法(RVA)对黄河上游兰州站的生态水文特征变异进行分析，给出了不同水库运行模式对黄河兰州站生态水文特征变异的影响程度情况。然而以系统性的思路对黄河上游兰州断面径流序列进行变异诊断、分级和指标识别的研究还较少。本文基于黄河流域水文生态变异问题的重要性及其特殊意义，对研究方法及思路作了进一步的考虑：采用TFPW-MK突变检验法和秩和检验法2种方法分析黄河上游流域代表性断面径流变异，避免单一变异分析方法的缺陷，增加变异分析结果的可信度；进一步采用IHA水文指标改变法分析变异前后水文生态情势。

1 研究区概况

黄河上游区从源头至内蒙古托克托县河口镇，河道长3 472 km，流域面积约38.6 万km2，汇入的较大支流有43条，径流量占全河的60%，其水文生态的变化对区域的可持续发展具有重要的影响。根据已有成果[7-9]中水文分区与生态分区的特点，考虑不同位置的调蓄水库的影响控制范围以及干流控制断面的分布，将黄河上游划分为黄河源头-兰州河段、兰州-河口镇段、兰州-河口镇段。本文选取刘家峡水电站下游的兰州断面，其流量的变化可以反映龙羊峡、刘家峡联合调度对流量特征的影响，将对下游河段的供水灌溉、水沙关系以及生态保护等产生关键性的影响。

2 方 法

2.1 去趋势预置白Mann-Kendall检验法

基于秩的Mann-Kendall检验法[10,11]是一种非参数统计检验方法，除了可以进行趋势检验，也可用于变异点的检验分析。由于其不需要遵从一定的分布，不受少数异常值的干扰，且结构简单，计算方便，使其也成为水文序列的变异点分析中常用方法之一，但其检验受到趋势项的影响。章诞武等[12]通过对比分析普通Mann-Kendall检验法和去趋势预置白Mann-Kendall检验法得出：序列中存在的正(负) 自相关性将增大( 减少) 其方差，从而在某一显著性水平下，增大( 降低) 拒绝原假设的概率；同时，序列趋势项会增大其自相关系数的计算值，对于有明显趋势的序列，应先剔除趋势项后估算其自相关系数。

去趋势预置白Mann-Kendall检验 (TFPW-MK, Mann-Kendall test with trend-free pre-whitening)[13,14]： 考虑数据序列明显趋势对自相关系数估计的影响，通过去趋势预置白方法，剔出趋势对自相关系数估计的影响，更准确地对数据序列进行检验。步骤如下。

(1) 假定序列由线性趋势和AR(1)组成，采用Theil(1950年)和Sen(1968年)提出的方法计算新样本数据的线性趋势β：

β=median( ) ∀j

(1)

(2)去除趋势项，形成不含趋势项的序列Yt：

Yt=Xt-Tt=Xt-βt

(2)

(3)计算序列Yt的一阶自相关系数r1，剔除序列中的自相关项：

Y′t=Yt-r1Yt-1

(3)

(4)补还趋势项Tt，得到不含自相关影响的新序列：

Y″t=Y″t+Tt

(4)

该序列将不再受自相关性的影响。

对经过去趋势预白化后形成的新序列使用Mann-Kendall检验法进一步检验序列变异，定义统计量为：

(5)

(6)

(7)

E(sk)=k(k+1)/4

(8)

(9)

将时间序列x按降序排列，再按上式计算，同时使：

UFk=-UB′k,k′=n+1-k,k=1,2,…,n

(10)

通过分析统计序列UFk和UBk，不仅可以进一步分析序列x的趋势变化，还可以明确突变的时间，指出突变的区域。若UFk值大于0，则表明序列呈上升趋势，小于0则呈下降趋势。2统计序列构成的曲线分别标记为UF和UB，如果2条曲线超过临界直线，则表明上升或下降趋势显著；如果2条曲线出现交点，且交点在临界直线之间，那么交点对应的时刻就是突变开始的时刻。

此方法与别的检验方法不同的是，在检验时，其2条曲线的交点通常为与2个年份之间，即突变点不指向某一个确定的点。我们需要采取一定的方法从2个点中选取一点作为突变点。

2.2 滑动秩和检验(Mann-Whitney U检验)

秩和检验法，也称Mann-WhitneyU检验法[15,16]，是非参数统计方法，它不依赖于总体分布类型，不以推断总体参数为目的，而旨在检验2种或2种以上观察变量的分布有无显著差异。由于秩和检验法对变异点诊断有适应性强且效率高(均值和Cv、Cs)的特点，本文采用该法来验证TFPW-MK检验的结果是否合理。滑动秩和检验法计算步骤如下[15-17]。

(1)考察点前后2序列总体的分布函数Fpre(x)和Fpost(x)，从总体Fpre(x)和Fpost(x)中分别抽取容量为npre和npost的2个样本，要求检验原假设Fpre(x)=Fpost(x)。

(2)将2组数据分别由小到大排列，再将2组数据由小到大统一编秩。如有原始数据相同时，可按相同数据全在同组内不用求平均秩次，不同组间有相同数据需求平均秩次的原则处理。

(3)求秩和并确定检验统计量。当2样本容量不同时，以容量小者容量为n1，容量大者容量为n2，小容量样本的秩和为统计检验量T1，大容量样本的秩和为统计检验量T2。若2样本容量相同，则n1=n2，任取一组的秩和检验统计量为T1，另一组为T2。

(4)计算和作出推断结论。统计量U值计算公式如下：

(11)

或

(12)

式中：U为统计量；n1为较小容量样本容量；n2为较大容量样本容量；T1为小容量样本的秩和；T2为大容量样本的秩和。

如果n2≤8，则分别用式(11)和式(12)计算，取较小值作为U值。然后查阅Mann-WhitneyU检验中U的相伴概率表得相伴概率p。

如果n2>8，无法查阅相伴概率标，可把正态近似用于U的抽样分布来检验：

(13)

式中：z为统计检验量；其他参数意义同前。

2.3 水文变异的影响域识别方法

水文生态变异指标的定义：在有关生态水文联系的研究中，人们通常用指标体系来表征和评估一个流域或区域的生态水文系统健康状态。指标体系一般由一系列具有生态意义的水文要素或特征所组成，并作为生态水文系统对水需求的外在表现形式。

水文变异法(RVA, Range of Variability Approach)[18]是对扰动河流进行水文情势变异评估的一种方法，其基本原理在于以变异点位置为依据，进行水文阶段划分，结合水文变异指标(IHA, Indicators of Hydrologic Alteration)，识别变异前后水文序列的变异领域(变异点前称为参考序列，变异点后称为变异序列)，对水文变异进行细节解析，找出发生强烈变异的水文指标。其以与生态相关的流量特征的统计分析为基础，从量、时间、频率、延时和变化率5个方面对32个河流水文特征(指标)进行提取(见表1)。指标的选取应遵循以下原则。

(1)立足流域河情和生态水文特点，充分考虑河流的水流特性，选取能够反映生态水文特征的指标。

(2)考虑已有的开发既定事实。指标的选取不仅仅以满足河流生态需求为标准，同时还要兼顾生防洪、供电的需要。

(3)重视指标的代表性。不同断面的水文特征不同，选取的指标要具有独立性和代表性，同时要避免不必要的重复。

(4)提高指标的可操控性，为后续的生态调控奠定基础。

(5)考虑指标的获取难度，保证在一定时期能够获取全部指标用以评价或表征。

将参考序列的提取结果作为扰动阈值， 通过对比评价序列的河流水文特征落入扰动阈值的数量， 来反映河流流量受扰动的影响程度，一般以32个IHA 指标以及综合影响程度的等级来表示，其中各个IHA指标代表的生态环境影响不同。水文指标是通过其参数的改变来反映影响生态变化的程度，水文指标的高程度改变无疑会导致相应的生态体系发生相对应的高强度扰动。表1为各个指标具体含义及其对河流生态系统的影响。

表1 河流的主要生态水文改变指标(IHA)及其对河流生态系统的影响Tab.1 Summary of IHA parameters in rivers and their ecosystem influences

2.4 方法耦合

本文采用TFPW-MK检验法为主要变异检测方法，结合滑动秩和检验方法验证检测结果，得到变异年份后作为RVA法的间隔年份对IHA变异指标进行分析，见图1。

图1 方法耦合Fig.1 Method coupling

3 实例分析

选取兰州断面实测日径流资料，资料来源于黄委水文资料汇编，序列长度为1950-01-01-2004-12-31，共55 a。本文将对黄河上游的生态水文特征变化包括径流年际变化、汛期变化、年内分配、极端流量变化、水文极值发生时间、高流量和低流量、水流条件变化频率与速率、大小洪水进行分析。图2为兰州断面年径流变化趋势。

3.1 变异点的识别分析

鉴于变异点检验的特殊性和复杂性，采用 TFPW-MK变异点检验对兰州断面实测年径流序列进行变异点研究。检验结果如图2所示，可以看到，所得趋势线UFk数值基本都为负值，表示径流序列整体呈现连续的下降趋势，且在1985年径流序列有显著下降趋势。根据文献分析[19]，从20世纪80年代后期到90年代末，黄河上游气温逐步升高，降水量也呈现明显减少趋势， 降水量减少自然会导致径流减少，且径流变化比例大于降水量变化比例。因此，黄河上游径流趋势和突变特征在一定程度上是对流域内降水变化的响应，特别是年径流序列20世纪80年代中期与流域降水突变时间相吻合。

图2 兰州站年天然径流量变化趋势线Fig.2 Annual natural runoff change trend line of Lanzhou station

依据TFPW-MK突变检验结果，UFk和UBk曲线交点位于1985年，即初步认为兰州断面径流变异年份为1985年。

图3 TFPW-MK变异检验Fig.3 Mann-Kendall test with trend-free pre-whitening

3.2 滑动秩和检验

依照时间顺序，对水文序列逐点进行秩和检验，即逐点作为分割点，以此把水文序列分成2个样本。设定2个假设，即：①H0，2个样本的分布相同，即无显著差异；②H1，2个样本的分布不同，即具有显著差异。

取置信水平α=0.05 ，采用双侧分位数检验，查得α=0.05 时，Z1-α/2=1.96 ，若|z|≤Z1-α/2，则接受H0，2个样本无显著差异；若|z|>Z1-α/2，则接受H1，2个样本有显著差异，即发生变异，并以|z|最大处作为水文变异点。表2为兰州断面滑动秩和检验成果。

表2 兰州断面首个变异点滑动秩和检验成果Tab.2 The moving rank sum test results of the first change point of Lanzhou

步骤1：对兰州整个序列逐年进行秩和检验，计算出首个变异点。序列分段为1950-1985、1986-2004 时，统计检验量|z|最大，且|z|>1.96，检验结果为拒绝H0，2样本分布有显著差异。

步骤2：接着找出所有可能变异点。对1950-1985年进行滑动秩和检验，序列分段为1950-1963和1964-1985时，统计检验量最大，检验结果为拒绝H0，2样本分布有显著差异；对1986-2004 年进行滑动秩和检验，该时段内不存在变异点。现在，兰州序列分段为1950-1963，1964-1985，1986-2004。

步骤3：每相邻两段序列整合成一个序列，进行变异点分析，对1964-2004年进行滑动秩和检验，序列分段为1964-1985、1986-2004时，统计检验量最大，检验结果为拒绝H0，2样本分布有显著差异；对1950-1985年进行滑动秩和检验，序列分段为1950-1963、1964-1985时，统计检验量最大，检验结果为拒绝H0，2样本分布有显著差异。每相邻两段序列合并成为一个序列进行检验，其变异点均没有发生变化。需要再次用步骤2对1950-1963年进行滑动秩和检验，该时段内不存在变异点；对1964-1985年进行滑动秩和检验，该时段内不存在变异点；对1986-2004年进行滑动秩和检验，该时段内不存在变异点，各序列均没有变异点。已经找出所有可能的变异点，最后得出唐乃亥的水文变异分段为：1950-1963、1964-1985、1986-2004。

根据图3所示TFPW-MK检验结果与秩和检验结果对照分析，在1963年，UFk和UBk曲线没有交点，且趋势线UFk变化较为平缓，则认为在1963年未发生变异。而在1985年，TFPW-MK检验结果与秩和检验结果吻合，均显示发生水文变异，则确定1985年为黄河上游兰州断面径流序列的突变年份。

3.3 IHA指标分析

通过以上分析，确定了水文序列的变异时间、变异强度，但是各水文指标所发生的变异程度仍旧未知，故本文通过RVA法评价了变异点前后序列的变异范围(指标)以及强度，结果如图4所示。

图4 指标改变度分析Fig.4 Alteration analyzing

通过分析可知，兰州断面的生态水文特征变化较大，整体改变度为0.61，属高度改变。32个指标中发生高度改变的近1/2，并且大多集中在7、8、9、10月平均流量、年最大流量平均值、每年高流量次数和延时以及每年流量逆转次数。高度改变的水文指标改变度如表3所示。

表3 高度改变的水文指标Tab.3 Highly changed alterations

以水文指标的改变为依据，来分析其对生态系统的影响。

(1)7、8、9、10月属于汛期，植物生长茂盛，蒸发强度大，因此对水的需求量大，1985年后7、8、9、10月平均流量高度改变，尤其是9月份从1 583 m3/s减少到863 m3/s，改变度达到-0.1，这对下游依赖洪水流量及其创造的栖息地环境而生存的生物物种而言，影响是破坏性的。

(2)最小化指标整体趋势趋于平缓，最大化指标则受其削减洪峰的影响，在1985年后出现了锐减，说明大洪水的洪峰一般只持续1～5 d， 最大流量平均值的高度改变会在一定程度上影响河流渠道地貌和自然栖息地的构建。

(3)高脉冲流量通过与漫滩和高地的连通,可大量地输送营养物质并塑造漫滩多样化形态，维系河道并育食河岸生物，从而提高河流生物的生物量和多样性。1985年后高脉冲流量的次数和延时均有明显下降趋势，将会导致生物多样性的下降。

(4)每年流量逆转次数明显增多，而水流条件的变化频率对干旱地区河流生物物种影响显著，陡涨陡落的流量变化会导致水生生物被冲刷或搁浅，洪水的暴落导致生物幼苗种群不能建立。由以上分析可以看出1985年后兰州断面流量的主要变化可以概括为汛期流量下降、洪峰削减、高脉冲流量减少和水流条件变化频率加快，这些改变与兰州断面受龙羊峡、刘家峡水库不同组合调节联系密切。1985年以后，龙羊峡水库建成，龙刘联合调度使得防洪能力和调蓄能力进一步增强,蓄丰补枯，汛期水流被削平，枯水流量被填补，峰谷缺失，以致水流情势改变严重。7-10月份流量过程趋于平缓，最大化流量平均值发生较大改变。同时水库兼有供水灌溉、发电、防洪任务,相机供水的次数增多，水流流量起落变化频繁。

4 结 语

通过对黄河上游兰州断面的实测径流水文序列变异点的识别和分析可得出以下结论。

(1)本文所采用的方法经过耦合从多方面对水文序列进行了检验，较全面地反映了时间序列的变异特性，因而检验分辨率较高。

(2)通过变异诊断分析，发现兰州断面水文生态指标在龙、刘两库运行后发生较高程度改变，这无疑非常不利于区域水文生态安全，建议在黄河生态水文指标体系的构建和水库生态调度中给予充分考虑。

(3)本研究以黄河流域兰州断面为例，研究了黄河上游生态水文变异情况，取得了一定的成果，但在实践应用层面仍存在问题，比如如何使用生态水文变异分析结果指导区域的水资源管理实践，在后续研究中需重点考量。作者认为生态水文变异分析应该成为河流水资源管理和水利工程规划中的重要限制性因素。因此，建立生态水文联系变异指标与水事活动间的关系，并依据其关系设置限建规则，将成为生态水文变异分析应用于流域水资源管理实践的重要手段。

□

[1] IPCC.The physical science basic[C]∥Contribution of working group I to the fourth assessment report of the intergovernmental panel on climate change(IPCC). Paris, 2007.

[2] Kundzewicz Z W,Graczyk D. Trend detection in river flow series:1. Annual maximum flow [J]. Hydrological Sciences Journal, 2005,50(5):797-810.

[3] 王金花， 康玲玲，余 辉，等. 气候变化对黄河上游天然径流量影响分析[J]. 干旱区地理, 2005,8(3):88-291.

[4] 樊 辉，杨晓阳. 黄河干、支流径流量与输沙量年际变化特征[J]. 泥沙研究, 2011,(4):11-15.

[5] 李剑锋，张 强，陈晓宏，等. 考虑水文变异的黄河干流河道内生态需水研究[J]. 地理学报, 2011,66(1):99-110.

[6] 于 茜，黄 强，张洪波. 不同水库运行模式对黄河兰州站生态水文特征变异的影响[J]. 干旱区地理, 2010,33(5):747-755.

[7] 水利部黄河水利委员会.黄河流域及西北内陆河水功能区划[R]. 黄河水利委员会,2001.

[8] 黄河流域气候编写组.黄河流域气候[R].黄河水利委员会水文局,1987.

[9] 沈珍瑶，杨志峰.黄河流域水资源可再生性评价指标体系与评价方法[J].自然资源学报,2002,17(2):188-197.

[10] Kendall M G.Rank correlation measures[M]. London:Charles Griffin,1976:110-120.

[11] Mann H B.Non-parametric tests against trend[J].Econometrica,1945,13:245-259.

[12] 章诞武，丛振涛，倪广恒. 基于中国气象资料的趋势检验方法对比分析[J]. 水科学进展, 2013,24(4):490-496.

[13] Sheng Yue,Wang Chunyuan.The mann-kendall test modified by effective sample size to detect trend in serially correlated hydrological series[J].Water Resources Management,2004,18:208-218.

[14] Sheng Yue,Paul P,Bob P,et al.The influence of autocorrelation on the ability to detecttrend in hydrological series[J].Hydrological Processes,2002,16(1):1 807-1 829.

[15] Seigel S. Nonparametric statistics for the behavioral sciences[Z]. New York: McGraw-Hill Book Company, 1956.

[16] Mann H B, Whitney D R. On a test of whether one of two random variables is stochastically larger than the other[J]. The Annals of Mathematical Statistics, 1947,18(1):50-60.

[17] 拜存有. 渭河流域关中段径流过程变异点诊断研究[D]. 陕西杨凌: 西北农林科技大学,2008.

[18] Brian D Richter,Jeffrey V Baumgartner,Jennifer Powell,et al.A method for assessing hydrologic alteration within ecosystems[J].Conservation Biology,1996,10(4):1 163-1 174.

[19] 赵芳芳，徐宗学. 黄河兰州以上气候要素长期变化趋势和突变特征分析[J].气象学报，2006，64(2)：246-255.

[20] 冯 波，闫佰忠，章光新. 松花江流域水文干旱联合概率分布特征研究[J]. 节水灌溉，2014,(5):38-42.