基于ERHIs法的雅鲁藏布江流域水文情势变化研究

杨 靖，蔺子琪，赵庆绪，谢 静，廖卫红

(1.华能西藏水电安全工程技术研究中心，四川 成都 610093；2.四川大学 水利水电学院，四川 成都 610065；3.华中科技大学 土木与水利工程学院，湖北 武汉 430072；4.长江勘测规划设计研究有限责任公司，湖北 武汉 430010；5.水利部水网工程与调度重点实验室，湖北 武汉 430010；6.中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室，北京 100038)

1 研究背景

雅鲁藏布江是世界上海拔最高、西藏最大的一条河流，也是我国西南地区非常重要的跨界河流之一。雅鲁藏布江流域径流主要由降水、冰雪融水和地下水补给组成[1]。由于流域面积广阔，流域内的降水量、蒸发量以及地形、冰雪覆盖面积、植被等条件差异很大，流域内的径流补给类型复杂多样[2-3]。精确描述径流特征及其变化的基础是获得大量高精度流量数据，特别是在水文站点缺乏的地区。分布式水文模型为流域径流数据的延展和加密提供了重要的支撑[4-5]，基于此，可量化评估人类活动、下垫面和气候变化等变化条件对水文特性影响程度，从而可为区域水资源合理开发利用提供重要的依据[6]。

针对气候变化和水资源开发利用而导致的水文改变度的分析一直为生态水文研究领域的热点问题[7-8]，Richter等[9]在1996年提出了一套能够评估水生态水文改变过程的指标体系(Indicators of Hydrologic Alteration，IHA)，目前被广泛应用于生态水文领域[10-14]。但Olden等[15]学者分析结果表明，众多水文变量之间存在十分突出的冗余性问题，而IHA指标与大多数水文变量之间存在较高的信息重叠现象。相对于早期研究中只考虑单一水文变量，尽管IHA体系较完善，但是水文变量仍然较多，IHA变量之间也存在着自相关性和严重的冗余性问题。IHA法目前被认为可以相对最全面、最系统的来评估水文情势变化特征。此外，雅鲁藏布江地区的水文情势研究较多，但具体对于环境流量组分的变化分析较少[16]，本文采用改进后的IHA法[17]降低生态水文变量中的冗余性，保留原始数据中的关键信息，分析雅鲁藏布江下游、尼洋河出口断面受上游水库影响的水文改变度，聚焦于对流量组分分级以及缺资料地区的水文情势评估。以期为分析雅鲁藏布江地区水资源管理和生态治理提供科学参考。

2 研究概况

2.1 研究区概况雅鲁藏布江(28°N—31°N、82°E—97°E)发源于西藏自治区西南部喜马拉雅山中段北麓的杰马央宗冰川，是世界上海拔最高、西藏最大的一条河流。雅鲁藏布江从源头～入海口全长约3162 km，其中，中国境内长约2057 km。入海口处总流域面积175万km2，在中国境内流域面积约33.20万km2，雅鲁藏布江径流总体以降水补给为主，流域内多年平均降水量300～500 mm，且年际变化较小，年内分配不均。流域多年平均净流量约1655亿m3，年径流深与年降水量的分布趋势具有较好的一致性[18]。测站分布位置如图1所示。

图1 研究区地理位置和观测站点示意图Fig.1 Schematic diagram of geographical location and observation stations in the study area

2.2 研究数据奴下水文站位于雅鲁藏布江流域干流，控制流域面积占其国内总流域面积的57.6%。尼洋河出口断面的更张测站为雅鲁藏布江的上游控制站，是青藏高原重要的水文站点之一[19]。本文根据雅鲁藏布江流域水文站点分布及各站基本资料条件，选用奴下、更张站作为径流计算的依据站，进行分析评价采用的是1965—2018年的径流日值数据。然而，由于资料获取难度问题，本文分析的几个断面的数据均存在不同程度的资料缺失，如奴下站点缺1958—1959年、1969年的实测资料。因此，本文采用分布式水文模型EasyDHM[20-22]进行主要站点逐日流量数据的插补还原。

2.3 研究方法论文基于Richter等[9，23-24]提出的IHA指标进行分析，通过重新构建IHA指标体系将降水时期划分为枯水期(1—5月和11—12月)和丰水期(6月—10月)两个时段，将IHA指标划分为不同时段内IHA指标，进一步优化IHA体系计算生态最相关水文指标 ERHIs，从而达到减少水文指标的同时尽可能多地保留原数据关键信息[17]。然后基于ERHIs进一步优选参数，采用6种突变检验诊断指标变化规律，最后应用水文改变度描述雅鲁藏布江径流变化情况。

2.3.1 水文改变指标(IHA) Richter提出的IHA体系，该方法定量分析了水文变异程度，目前该指标体系被普遍用来代表水文情势变化特征，该指标体系包含33个指标，目前被认为可以相对最全面、最系统的来评估水文情势变化特征。该体系共分为5组，详细指标划分可参照文献[23]。由于研究时段内雅鲁藏布江未出现断流现象，故本研究不考虑“断流天数”这一变量[25]。

2.3.2 突变诊断 采用贝叶斯推断[26]、曼-肯德尔(Mann-Kendall)法[27-28]、勒帕热(Le-Page)法[28]、BG(Bernaola-Galvan)分割算法[29]以及有序聚类法[30]对改进后的指标突变点进行识别，用以排除单一突变诊断方法的局限性以及数据漂移问题。

(1)贝叶斯推断

贝叶斯推断[26]是利用贝叶斯定理来修正特定假设的可能性，该方法可以充分利用先验信息和样本信息，得到统计结论。计算公式为：

p(X，α)=p(X|α)p(α)=p(α|X)p(X)

(1)

(2)

p(α|Y)∝p(α)p(Y|α)

(3)

(2)M-K法

M-K法[27-28]作为常见的非参数统计检验方法，无需样本遵从一定的分布，适用于类型变量和顺序变量，计算简便，可以明确突变开始的时间并指出突变区域。计算公式如下：

(4)

(5)

(6)

其中，ai是第i时刻数值个数的累计数，E(ai)，D(ai)分别是ai的均值和方差，统计量UFi为标准正态分布。按时间序列x逆序xn，xn-1，…，x1，再重复上述过程，同时使UBk=-UFk，其中k=n，n-1，…，1。

2）经营健康稳健，盈利质量和运营效率提升。公司过去三年每年经营活动现金流量净额均高于净利润，业绩质量较高。公司运营效率这两年也有很大改善，存货周转率不断优化，2016年和2017年存货周转率分别为1.69次和2.07次，今年上半年已经达到1.07次，较17年继续提升7%，运营效率的提升将会显著减少存货减值准备计提，对业绩提升有很大帮助。

(3)勒帕热法

勒帕热(Le Page)法[31]是一种具有两个样本不受分布影响的非参数判别法，计算公式如下：

(7)

(8)

U为统计量，式(8)是勒帕热统计量。注意到当样本量足够大时，WU渐近具有自由度为2的χ2分布。其余参数含义同上。

(4)BG分割算法

BG分割算法[29]主要针对非平稳时间序列进行突变检测，其核心思想是通过显著性检验找到突变点的左右两部分序列的差异最大的位置。计算公式如下：

(9)

(10)

(5)有序聚类法

有序聚类分析法[32]进行突变点的识别，该方法核心思想是在分类时不打乱次序的情况下，通过计算突变点前后点位的离差平方和寻求最优分割点。计算公式为：

2.3.3 水文改变度 基于IHA指标体系，Richter等[24]在1997年提出了用于判断IHA指标改变程度的变化范围法(Range of Variability Approach，RVA)，并定义了RVA阈值概念，即在IHA的基础上，将各指标发生几率的75%和25%作为RVA的上下阈值，用来作为评价各指标受影响程度的标准[33-35]。为了量化ERHIs受干扰后的变化程度，采用水文改变度来评估[36]，其定义如下：

式中：Di为第i个ERHIs指标的水文变化程度；Nio为第i个ERHIs指标变化在之后(之前)仍在RVA阈值范围内的实际观测年数；Nie为第i个ERHIs突变(前)后预计落在RVA阈值内的年数，可以表示为r×NT，其中r为突变前(后)第i个ERHIs落在RVA阈值内的比例(本文取50%)，NT是受突变(前)和后影响的时间序列中记录的总年数。D0为流域整体的水文改变度，规定水文变化程度在0～<33%之间为无变化或低程度变化，33%～<67%为中度变化，67%～≤100%之间是高度变化。

3 模型验证与指标确定

3.1 测站验证依据流域和数据资料情况，论文对奴下、更张断面的1999到2012年水文日模拟精度进行了参数率定，将优化得到的参数对流域各主要断面1999—2018年的逐日水文过程进行了模拟，各站Nash系数评价见表1，说明论文采用的EasyDHM模型用于日径流资料的延展和加密是可行的。依据流域和数据资料情况，本研究对奴下、更张断面挑选洪水场次进行了参数率定和验证，并对精度进行了评价，针对流域特点，主要针对洪峰误差精度对方案合格程度进行了评定，其中水库入库数据主要通过出库和蓄量反推，本身就存在误差，因此对其洪峰许可误差做了适当放大。选取1998—2018年间的20组径流过程进行率定，径流过程拟合效果较好，洪峰流量误差和峰现时间误差基本符合要求。对1998—2018年间的20组径流过程样本进行评定，洪水场次合格率为90%。EasyDHM模型采用SCE-UA全局参数优化算法进行参数率定[37-40]。

表1 EasyDHM模型模拟径流Nash系数表Table 1 Nash coefficient of simulated runoff by EasyDHM model

3.2 IHA的冗余性分析IHA体系包含33个指标，由于研究时段内雅鲁藏布江未出现断流现象，故本研究不考虑“断流天数”这一变量[25]。论文对所选择的32个IHA指标进行冗余性分析，主要针对指标间的相关系数，认为相关系数小于0.4为冗余度较小。从图2可以看出，部分水文变量之前展现了较强的相关性，10月流量，11月流量，12月流量，年最小 1 日、3 日、7 日、30 日流量之间的相关性较强，高流量年内发生次数，高流量年内平均历时，上升率，下降率之间相关性较强。从图3可以看出，年最小 3 日、7 日、30 日流量，年最大 1 日、3 日、7 日、30 日、90日流量之间也具有较强的相关性，其相关系数甚至超过了0.9。结果表明，雅鲁藏布江的站点断面的IHA指标变量相互之间显然存有高度自关性，即指标冗余性问题非常明显。

图2 奴下断面IHA指标相关关系热力图Fig.2 Correlation heat map of IHA index under Nuxia section

图3 更张断面IHA指标相关关系热力图Fig.3 Correlation heat map of IHA index under Gengzhang section

3.3 IHA指标体系重建鉴于IHA变量之间存在的冗余性问题比较突出，本文通过重新构建IHA指标体系将降水时期划分为枯水期(1—5月和11—12月)和丰水期(6月—10月)两个时段，将IHA指标划分为不同时段内IHA指标，进一步优化IHA体系计算ERHIs，在缓解IHA变量冗余性时，减少水文指标，在保证全面描述流域水文情势变化时尽可能多地保留原数据关键信息，通过涵盖丰枯时段以及指标所在组别信息来确定。

流域的径流量以及径流特征受到季节变化的影响，由于不同季节降水量的不同，可将降水时期分为丰水期和枯水期。由于一年内枯水期历史长短与流域自然地理及气象条件相关，根据丰水期和枯水期的特点将雅鲁藏布江奴下断面6月至10月划分为丰水期，将1—5月、11—12月划分为枯水期(表2)。

表2 重建后IHA指标Table 2 IHA indicator after reconstruction

为使筛选出的ERHIs指标具有较高的可靠性，本文通过选取样本年并挑选各样本年所对应的相似年，计算相似年各IHA指标的贡献率，选取贡献率排在前六位且重复出现次数最多的IHA指标为ERHIs指标。

相似年通常指在表征同一流域有两年或多个年份具有相同的水文指标或径流特征或其他流域特征时，都可以称这两年或多个年份之间在某一特征前提下存在相似性，即这两年或多个年互为相似年[41-42]。相似年在各个领域的应用十分广泛，在天文、气象、水利等各行业的预报都有应用。本文从奴下、更张断面57年数据中分别随机选取十年作为相应断面的样本年，将能够反应径流特征的水文统计值作为选取样本年所对应相似年的基准，运用二维欧氏距离筛选方法结合统计值对选取好的样本年进行相似年份的选取[41-42]。

为对应不同时期IHA指标体系，将1962—2018年流量数据按枯水期(1—5月和11—12月)和丰水期(6月—10月)分开，分别计算57年数据的枯水期统计值和丰水期统计值。统计值包括：均值、方差、极差、中位数、25分位数、75分位数、变异系数、偏度和峰度。

3.4 ERHIs指标确定通过计算各样本年所对应十个相似年枯水期和丰水期IHA指标的变异系数，将十个相似年中各年枯水期和丰水期变异系数前六个指标挑选出来并记录各指标在十年中出现的重叠次数，挑选出枯水期和丰水期出现次数排在前六的指标，将这十二个指标再次按照出现次数排序，查看前六个指标在一级指标中的覆盖度是否合理。通过计算各样本年所对应十个相似年枯水期和丰水期IHA指标的变异系数，将十个相似年中各年枯水期和丰水期变异系数前六个指标挑选出来并记录各指标在十年中出现的重叠次数，挑选出枯水期和丰水期出现次数排在前六的指标，将这十二个指标再次按照出现次数排序，查看前六个指标在一级指标中的覆盖度是否合理。经合理性确认，将前六个指标确定为ERHIs。在合理性判定的基础上计算各指标的相关性，经计算后发现只有两个指标存在相关性，其他指标之间的相关系数都小于0.4为不相关，将这些指标确定为筛选指标。具体步骤见图4。

图4 ERHIs指标确定步骤Fig.4 ERHIs indicator determination steps

指标的变异系数反应该指标在所有IHA指标中的鉴别能力。变异系数越大，表明该指标的信息分辨能力就越强，从而保证了筛选出的指标对流域影响最大[43]。初步筛选出的ERHIs指标见表3、表4。

表3 更张断面ERHIs指标Table 3 ERHIs index of Gengzhang section

表4 奴下断面ERHIs指标Table 4 ERHIs index of Nuxia section

3.4.1 ERHIs指标的合理性分析 为验证2个断面选出的ERHIs指标相互独立且具有较高的代表性，现对选出的ERHIs指标进行合理性分析，本文中主要通过ERHIs指标的覆盖性分析和冗余性分析来进行合理性分析。在减少IHA变量冗余性的同时，以保证ERHIs指标可以均匀覆盖不同时期以及不同组别，有效保留原始信息。

(1)ERHIs覆盖性分析

覆盖性分析通过对各ERHIs指标在组别中的覆盖范围以及各指标在不同时期的覆盖范围进行指标的合理性分析，研究断面指标覆盖范围分析结果见表5、表6。

表5 更张断面覆盖性分析Table 5 Coverage analysis of Gengzhang

表6 奴下断面覆盖性分析Table 6 Coverage analysis of Nuxia

由表可以得出，更张断面六个ERHIs指标涉及到了第一组、第三组至第五组且达到丰枯水期全覆盖，奴下断面的六个ERHIs指标涉及到了第一组至第五组指标且达到丰枯水期全覆盖。可以看出雅鲁藏布江下游的2个断面挑选出的ERHIs指标在组别中的覆盖率以及各指标在不同时期的覆盖率都很高，得出挑选出的ERHIs指标具有一定的合理性。

(2)ERHIs冗余性分析

为验证5个断面选出的ERHIs指标相互独立且具有较高的代表性，现对选出的ERHIs指标进行相关性分析，确定ERHIs指标各指标间的相关性，图5展示了各指标间的相关系数。结果表明：大部分指标之间存在较小的相关性，即信息并不存在冗余性问题，则可最终确定为ERHIs指标。

图5 更张断面和奴下断面ERHIs相关关系热力图Fig.5 Heatmap of the correlation relationship of ERHIs in Gengzhang and Nuxia section

4 计算结果分析

4.1 初步诊断初步诊断是对已经挑选出的ERHIs指标进行指标变异强度的确定，将强变异指标确定为最终用于识别变异年份的指标。变异系统的初步诊断采用Hurst系数法[44]分别对雅鲁藏布江的更张水文站、奴下水文站的ERHIs指标进行变异检验。经检验得出尼洋河站的10月指标Hurst系数为0.90属于强变异，奴下站的4月指标Hurst系数为0.86属于强变异。计算结果见表7、表8。

表7 更张站ERHIs指标变异强度Table 7 Intensity of variation of ERHIs indicators in Gengzhang section

表8 奴下站ERHIs指标变异强度Table 8 Intensity of variation of ERHIs indicators in Nuxia section

4.2 详细诊断根据初步诊断结果针对各断面挑选出的强变异指标进行详细诊断。本文利用水文变异诊断系统的跳跃诊断法和气候突变检测方法包括Mann-Kendall方法、有序聚类法、贝叶斯推断方法、Le-Page方法和BG(Bernaola-Galvan)分割法等方法进行突变年份的确定。更张断面有4种方法诊断10月指标在2000年发生显著突变，有2种方法诊断其在2014年发生显著突变。奴下断面有3种方法诊断4月指标在2013年发生显著突变。具体诊断结果见表9。

表9 雅鲁藏布江各断面强变异指标的突变诊断结果Table 9 Mutation diagnosis results of strong variation indicators in each section of the Yarlung Zangbo River

结果表明，不同突变方法对突变点的识别有所差异，由于LePage方法需要人为确定子序列长度，可能会由于突变点的漂移造成误差，只用一种突变方法难以全面的评判指标突变的结果。综上所述，运用多种突变检验方法分析更张、奴下站57年径流数据得出以下结论：①对于更张断面：10月指标突变时段为1990—2014年。②对于奴下断面：4月指标突变时段综合序列为1970—2016年。造成奴下断面指标综合突变时间较长的原因可能为上游支流拉萨河直孔及旁多水库对下游的调蓄作用，2010年冲久水库投入运行后，改变了下游河段的天然径流过程。因而，奴下断面通过3种突变方法中体现出这种突变状态，其中Hurst系数指标显示年最小流量出现时间、年最小1日流量指标均有中等程度的改变。

4.3 强变异指标变化程度本文通过判断变异前后时间序列的长短来确定稳定时段，通过以上对奴下和更张断面水文情势突变点的研究，认为时间序列较长一段的水文情势更为稳定[24]。将更张断面1965—1989年作为稳定时段，将奴下断面1965—1969年作为稳定时段。以Richter等[26]定义的阈值概

念即将稳定时段强变异指标发生几率的75%作为RVA的上阈值，发生几率的25%作为RVA的下阈值作为标准。

论文分别计算各断面强变异指标包括更张断面10月指标、奴下断面4月指标的水文改变度(表10)。结果显示，两个断面的强变异指标都发生了高度改变。综上得出，更张断面和奴下断面整体水文改变度为高度改变。

表10 雅鲁藏布江2个断面强变异指标改变度计算结果Table 10 Calculation results of the change degree of strong variation index of two sections of the Brahmaputra River

5 结论

(1)通过EasyDHM插补雅江流域站点数据是可行的，其中，参数优化后的Nash系数最大为0.97，最小为0.67，均在合理范围内，在缺资料地区运用较好。此外，ERHIs方法有效的通过减少变量简化水文响应推求过程，以获得水文情势变化情况，提高水资源调配效率，对于水资源管理和生态系统保护都具有重要的参考价值和借鉴意义。

(2)ERHIs不仅有效缓解了IHA的冗余性问题，还有利于抓住最关键的生态水文变量。从IHA指标相关性热图可以看出，IHA参数之间的相关性较高，说明指标间的冗余问题十分突出。ERHIs指标突变结果分析表明，运用欧式距离结合相似年筛选IHA指标与多种突变结果的横向对比以及ERHIs自身的纵向对比，都表明所选取的生态最相关指标是合理的。

(3)气候变化和人类活动对环境组分流量和水文状况的构成都有影响，而大规模水利建设是其中的主导因素。针对更张断面和奴下断面的水文改变度分析可以得出，在受到多布水电站建设的影响，尼洋河下游水文改变度差异显著，这和张冬艳、刘海平等人的研究结果相一致[45-46]。

(4)上游水电站建设运行对水文情势有不同程度的影响，水文情势变化总体属高度改变，非汛期所受影响程度比汛期更大。干旱季节的径流变化较大，除受到上游水电站影响外，气温的增温效应导致冰川融雪的补给量增加也在影响雅鲁藏布江的年内径流分配，但这种变化还需相关的资料和分析去验证。