0 引言

汾河是黄河中游仅次于渭河的第二大支流，对黄河中下游的水沙演变过程有着重要影响。汾河发源于山西省宁武县管涔山和天池两处，流域面积39 471 km²，河长694 km，干流由北向南纵贯于云中山、吕梁山之间的谷地，经山西省河津市流入黄河。实测资料^[1]表明，近30年来，汾河径流呈锐减趋势，汾河入黄控制性水文站河津站1980年以来实测径流量相对于1920—1979年的平均流量减少高达66%。

汾河径流锐减原因非常复杂，其中包括水利工程建设、水土保持措施、煤炭资源开采等人类活动方面的原因，此外气候变化也是汾河径流锐减的一个重要原因。据1956—2000年资料^[2]统计，汾河流域的降水量呈递减趋势。如果以1956—1969年为基准期，基准期内年均降水量为528.8 mm，20世纪70、80、90年代年降水量较基准期分别减少了8.8%、13.6%、12.3%。汾河径流锐减带来了一系列后果，包括引发河道内水环境和水生态恶化、主河槽萎缩、入黄径流锐减等；因此，定量分析各个影响因素对汾河径流锐减形成的相对贡献，对提出具有针对性的区域旱、涝减灾措施有重要意义，这也是汾河流域水资源可持续利用与管理中迫切需要解决的问题。

气候变化使得流域降水、蒸发等水文要素发生改变，从而改变流域水资源的演变情势和河川径流过程。人类活动通过改变流域下垫面而影响产汇流机制和流域水循环过程，从而对河川径流产生显著影响。近20年来，变化环境下流域水循环过程演变过程问题及影响因素的量化问题，一直是国内外水文分析领域的研究热点。已有研究成果中的分析方法主要有：相似降水对比试验法^[3-5]、相似流域对比试验法^[6-8]、分项组合调查法^[9-11]、资料统计回归法^[12-14]、流域水文模拟法^[15-17]等。上述方法中，前三种方法在选择试验流域时存在较大的困难，并且对水文及相关资料的要求较高；因此，在实际运用中以后两种方法居多，尤其是流域水文模拟法，其近年来在评估气候变化和人类活动对河川径流过程影响问题中得到了大量的应用。但该方法的不足之处是水文模型的参数众多，率定参数所需的水文、气候、下垫面等实测资料量非常大，在准确率定水文模型参数方面存在一定困难；此外，尤其是1980年以前的土地利用/覆被变化(LUCC) 的历史资料难以获取，而我国大部分流域的下垫面在1980年以前已发生较大的变化，这使得该方法在实际应用过程中受到一定局限。

针对上述问题，本文以汾河流域作为研究区，提出构建流域降水径流多元动态回归分析方法，利用该方法定量分析降水变化和产流条件变化对汾河径流变化的相对贡献，以期为汾河流域水资源的合理开发利用、减轻旱涝灾害提供科学参考；同时在定量评价人类活动和气候变化对流域水文情势影响方法创新上也是一次探索和尝试。

1 研究思路和方法 1.1 降水径流双累积曲线法

双累积曲线法^[18-19]是检验两个参数间关系一致性及其变化的常用方法。所谓双累积曲线就是在直角坐标系中绘制的同期内一个变量的连续累积值与另一个变量连续累积值的关系线，它可用于水文所气象要素一致性的检验、缺值的插补或资料校正，以及水文气象要素的趋势性变化及其强度的分析。双累积曲线的理论基础为：在相同时段内只要给定的两个变量的两组样本数据成正比，那么一个变量的累积值与另一变量的累积值在直角坐标上可以表示为一条直线，其斜率为两个变量对应点的比例常数；如果双累积曲线的斜率发生突变，则意味着两个变量之间的比例常数发生了改变，或者其对应的累积值的比可能根本就不是常数；若接受两个变量累积值之间直线斜率已发生改变，那么斜率发生突变点所对应的年份就是两个变量累积关系出现突变的时间。

双累积曲线的绘制方法为：设有两个变量X(参考变量或基准变量) 及Y(被检验变量)，在N年的观测期内，有观测值X_i及Y_i，其中i=1，2，3，…，N。首先对变量X及变量Y按年计算各自的累积值，得到新的逐年累积序列X′_i及Y′_i，其中i=1，2，3，…，N，即

然后，在直角坐标第中绘制2个变量所对应点累积值的关系曲线。绘制的曲线图一般以被检验的变量为纵坐标(即Y轴)、参考变量或基准变量为横坐标(即X轴)。

1.2 动态回归分析法

受下垫面复杂状况的影响，不同降水要素下流域的径流系数会有所不同。以降水量不同时为例，当某年7月份流域的降水量为40 mm时，其与降水量为120 mm时的降水径流回归方程可能会有所不同。不妨假设在统计意义上存在某一阈值(假设为70 mm)：

当该年7月份降水量大于该阈值(70 mm) 时，其降水径流可以用回归方程(Ⅰ) 来表述；

当该年7月份降水量小于等于该阈值(70 mm) 时，其降水径流可以用回归方程(Ⅱ) 来表述。

当然，影响径流演变过程的因素非常多，例如蒸发、地下水位、水土保持面积等其他因素的状态改变时，径流的演变机制均可能发生改变；因此，在建立径流演变过程对主要影响因素的多元回归方程时，应考虑到各个主要影响因素状态发生改变时，上述多元回归方程会发生变化，即不同产汇流机制下需用不同的多元回归方程来描述。

基于上述考虑，本项研究将首先对降水进行丰枯划分，得到年降水的丰、平、枯三种状态下的年降水系列，然后与对应的年径流进行线性回归分析，从而建立起年降水年径流的动态回归方程：

式中：a₀、a₁、a₂、a₃、b₀、b₁、b₂、b₃、c₀、c₁、c₂、c₃为回归方程的回归系数；W_t为第t年的年径流量；P_t为第t年的年降水量；P_t－1为第t-1年的年降水量；P_t－2为第t-2年的年降水量；P_枯为年降水量丰枯划分的枯水阈值；P_丰为年降水量丰枯划分的丰水阈值。

2 汾河降水径流多元动态回归模型 2.1 流域产流条件变化分析

本项研究将首先运用双累积曲线法对河津站降水径流系列相关性发生显著改变前后的2个历史时期进行划分；然后分别建立降水径流相关关系发生改变前后的2个降水径流动态回归模型，运用这2个多元动态回归模型来推算不同降水条件在不同时期产流条件下的理论年径流量；并通过与实测年径流对比，来估算降水变化和产流条件变化对汾河径流减少的贡献。图 1是汾河入黄控制站河津站的降水径流双累积曲线。

由图 1可看出：河津站控制流域的产流条件前18个经验数据点可近似作为一个产流条件，后35个经验数据点可近似作为另一个产流条件；第18个数据点(1973年) 可以作为两个产流条件的分界点。因此，依据图 1对降水径流资料进行2个时期的划分：

时期1：1956—1973年，此阶段内汾河流域对应的气候和下垫面等水文要素形成的产流条件称之为产流条件1；

时期2：1974—2008年，此阶段内汾河流域对应的气候和下垫面等水文要素形成的产流条件称之为产流条件2。

2.2 年降水量的丰枯状态划分

在建立流域降水径流之间的多元回归方程过程中，应考虑到降水的丰枯状态发生改变时，降水径流回归方程可能会发生变化；因此，在研究流域降水变化对径流变化的影响工作中，需对流域的年降水量变化过程进行丰枯状态的划分。

丰枯评价因丰枯指标的不同而异，丰枯指标的确定有多种计算方法，如频率分析法、均值标准差法、模糊聚类法等。常用的丰枯划分方法主要是频率分析法和距平百分比法。本项研究采用距平百分比法对河津站年降水进行丰枯划分，距平百分比的计算公式为

式中：η_DP为距平百分比；P_i为某一年份的降水量；为降水量系列的多年平均值。根据《水文情报预报规范(SL2502000)》^[20]，具体的丰枯类型划分采用的标准如下。

丰水年：η_DP>15%；平水年：－15%≤η_DP≤15%；枯水年：η_DP＜－15%

依据式(4) 和上述标准对河津站年降水量的划分结果如表 1所示。

采用上述丰枯划分指标对河津站1956—2008年间的丰枯年份进行统计，对应的丰、平、枯三种年份的统计结果如图 2所示。

2.3 降水径流多元动态回归分析

运用多元线性动态回归方法，分别对汾河流域时期1和时期2的降水径流进行动态回归分析。在进行回归分析之前，首先对当年径流量、当年降水量、前1年降水量、前2年降水量进行相关分析，表 2是上述2个时期内当年降水量为丰、平、枯3种状态下的降水径流多元相关系数矩阵。

对表 2进行分析可得以下结论：

1) 时期1内汾河流域在丰水年份的当年径流量，与当年降水量和前1年降水量的相关系数均较大，即丰水年份的当年径流量受当年降水量和前1年份降水量的影响。分析其原因可知，丰水年份的次降水量可能较大，从而出现超渗产流形成河道内的次洪水过程，甚至在局部地区形成蓄满产流而汇入河道，因而河津站的当年径流与当年降水呈较强的正相关关系；同时，汾河的河川基流来源于前1年份的流域内降水形成的深层地下水，因此又与前1年份的降水有正相关关系。在时期2内，丰水年份的当年径流量与前1年份降水量的相关系数为负值，违背水文降水径流成正相关的水文原理，说明这两个变量之间是伪相关。究其原因是因为选取的样本太短(1974—2008年间丰水年份样本只有4个)，缺乏代表性；如果系列足够长，两者应该是正相关的，即汾河流域丰水年份的当年径流量仍然受当年降水量和前1年降水量的影响。

2) 时期1内汾河河津站平水年份的当年径流量与当年降水量呈正相关，说明在平水年份的某几次降水以超渗产流或蓄满产流最终汇入到了汾河干流；与前1年份降水量、前2年份降水量均为负相关，依据流域产流条件原理可知此种相关关系为伪相关。需要注意的是，汾河基流由前1年份的降水形成，那么平水年份的径流量与前一年份的降水量为什么会出现负相关?这是因为选取的样本太短(时期1内平水年份样本只有6个),缺乏代表性；如果样本系列足够长，应该还是正相关的。观察表 2中时期2内平水年份的年径流量可知，其与当年降水量和前1年降水量均为正相关，且相关系数也较大，进一步说明了汾河流域平水年份的径流量受当年降水量和前1年降水量的影响。

3) 时期1和时期2内，汾河河津站枯水年份的当年径流量与当年降水量的相关系数均为负值且很小，而与前1年的降水量相关系数均大于0.900；说明枯水年份的径流主要由前1年的降水量影响而与当年降水量基本不相关。分析其原因可知，汾河流域地处干旱半干旱地区，枯水年份的次降水量一般较小，绝大部分耗于蒸发或被包气带、植被截留，很难汇集到河道形成次洪水，而汾河在枯水年份的径流是河川基流，河川基流一般由深层地下水补给，深层地下水流动非常缓慢，一般由前1年份的降水补给，即汾河流域枯水年份的径流量只受前1年降水量的影响，这一原理在文献^[21]中也有较为详尽的阐述。

基于上述分析可知，河津站的河川基流一般由前1年降水量形成，汾河的年径流量中剔除基流部分的剩余部分一般由当年降水量形成。因此，对汾河流域降水径流进行回归分析时，应根据汾河流域丰枯状态进行区别处理：当流域降水为丰水和平水年份时，应取当年降水量和前1年份降水量对年径流量进行回归；当流域降水为枯水年份时，应只选取前1年份降水量为自变量对年径流量进行回归。表 3分别为两个历史时期内3种丰枯状态下的多元动态回归分析的参数估计结果。

表 3中：R为模型的可决系数，它可以反映多个自变量对因变量的解释程度，其数值越大，回归模型中自变量对因变量的解释程度越好；F为模型进行卡方检验的检验统计量；P值为模型通过方差检验的概率，即模型计算值与样本实测值之间的相关性是由偶然性造成的概率，当P＜0.05时，说明模型通过显著性检验。由表 3可知：

1) 时期1内，丰水年份回归模型的可决系数R为0.728，可知丰水年份的回归模型中的2个自变量也可以很好地解释因变量(解释的程度为72.8%)，二者的线性关系良好；模型的总体显著性检验P值为0.039，即只有3.9%的概率拒绝该模型成立，因此，接受该模型；平水年份回归模型的可决系数R为0.571，可知平水年份的回归模型中的一个自变量勉强可以解释因变量(解释的程度为57.1%)，二者的线性关系一般，模型拟合度一般；模型的总体显著性检验P值为0.082，即只有8.2%的概率拒绝该模型成立，因此，可以勉强接受该模型。枯水年份回归模型的可决系数R为0.983，接近于1.000，可知枯水年份的回归模型中的1个自变量可以很好地解释因变量(解释程度为98.3%)，二者的线性关系很高，模型拟合度很高：模型的总体显著性检验P值为0.084，即只有8.4%的概率拒绝该模型成立，因此，接受该模型。

2) 时期2内，丰水年份回归模型的可决系数R为0.028，可知丰水年份的回归模型中的一个自变量解释因变量的效果很差，解释的程度仅为0.028%，二者的线性关系很差；模型的总体显著性检验P值为0.834，即有83.4%的概率拒绝该模型成立，因此，不能接受该模型。平水年份回归模型的可决系数R为0.694，可知平水年份的回归模型中的1个自变量可以较好地解释因变量(解释程度为69.4%)，二者的线性关系较好，模型拟合度也较高；模型的总体显著性检验P值为0.000，即拒绝该模型成立的概率基本为0.000，因此，接受该模型。枯水年份回归模型的可决系数R为0.995，接近于1.000，可知枯水年份的回归模型中的1个自变量也可以很好地解释因变量(解释程度为99.5%)，二者的线性关系很高，模型拟合度很高；模型的总体显著性检验P值为0.005，即只有0.5%的概率拒绝该模型成立，因此，接受该模型。

由此，可列出汾河河津站1956—1973年间的降水径流动态回归方程组：

同理，可建立汾河河津站1974—2008年间的降水径流动态回归方程组：

3 降水和产流条件对汾河径流影响的定量分析 3.1 降水变化对汾河径流的影响分析

1) 降水变化在产流条件1下对径流的影响程度

时期1实测降水量以P₁表示，实测径流量(气候变化前的径流量) 以W₁表示，把时期2实测降水量P₂代入回归方程组(5)，得到产流条件1下的推算径流量(气候变化后的径流量)，以W′₁表示，将降水丰枯变化在产流条件1下引起的径流变化量以ΔW₁^P表示，相应的变化率以η₁^P表示，则有

将汾河流域的实测降水、径流资料代入式(5) 和式(7)，可得结果如表 4所示。

由表 4可知：相对于时期1的实测径流量，时期2的丰水年份的推算年径流量平均减少6.51%，平水年份的推算年径流量平均减少4.70%，枯水年份的推算年径流量平均减少50.68%；3种水平年加权平均后表明，与时期1的降水条件相比，时期2的降水条件在产流条件1下将使得河津站年径流量平均减少22.71%。

2) 降水变化在产流条件2下对径流的影响程度

时期1的实测年降水量以P₁表示，时期2的实测年径流量(气候变化后的径流量) 以W₂表示，把P₁代入回归方程组(6) 得到产流条件2下的推算径流量(气候变化前的径流量) 以W′₂表示，将降水丰枯变化在产流条件2下引起的径流变化量以ΔW₂^P表示，相应的变化率以η₂^P表示，则有

将降水、径流资料代入式(6) 和式(8)，计算结果如表 5所示。

由表 5可知：相对于由时期1降水量代入式(7) 得到的推算年径流量，时期2平水年份的实测径流量平均减少24.45%，枯水年份的实测径流量平均减少34.73%；两种水平年份加权平均后表明，时期1降水量演变为时期2降水量后，在产流条件2下将使河津站年径流量减少约18.95%。

3.2 产流条件改变对汾河径流的影响分析

1)1956—1973年降水条件下产流条件变化对径流的影响程度

1956—1973年实测降水量以P₁表示，实测年径流量(产流条件变化前的径流量) 以W₁表示，由P₁代入回归方程组(7) 得到的径流机制2下的推算径流量(产流条件变化后的径流量) 以W′₂表示，那么在1956—1973年降水条件下，汾河的产流条件由产流条件1变化为产流条件2后引起的径流变化量以ΔW₁^m表示，相应的变化率以η₁^m表示，则有

将降水、径流资料代入式(7) 和式(10)，计算结果如表 6所示。

由表 6可知，在1956—1973年降水条件下(平水、枯水共9个年份)，汾河流域的由产流条件1演变为产流条件2后，河津站的年径流量在此9年间累计减少53.937×10⁸ m³，平均每年减少5.993×10⁸ m³；相对于产流条件1下的实测年均径流量12.721×10⁸ m³，产流条件2下的年均径流量减少了47.11%。

2)1974—2008年降水条件下产流条件变化对径流的影响程度

1974—2008年实测降水量以P₂表示，实测径流量(产流条件变化后的径流量) 以W₂表示，由实测降水量P₂代入回归方程组(6) 得到的径流条件1下的推算径流量(产流条件变化前的径流量) 以W′₁表示，那么在1974—2008年间的降水条件下，汾河的产流条件由产流条件1变化为产流条件2后引起的径流变化量以ΔW₂^m表示，相应的变化率以η₂^m表示，则有

将降水、径流资料代入式(6) 和式(11)，计算结果如表 7所示。

由表 7可知：在1974—2008年间降水条件下，汾河流域的由产流条件1演变为产流条件2后，河津站的年径流量将在此35 a间累计减少252.758×10⁸ m³，平均每年减少7.222×10⁸ m³；相对与产流条件1下的年均径流量13.084×10⁸ m³，产流条件2上的年均径流量减少了55.20%。

结合表 4可知，时期2的降水条件在产流条件1下将使得河津站年径流量平均减少22.71%；又由表 7可知，时期2的降水条件下，相对于产流条件1，产流条件2使得河津站年径流量平均减少55.20%。即，相对于时期1降水条件和产流条件，时期2的降水条件使径流减少22.71%，时期2的产流条件使径流减少55.20%；那么降水变化的相对贡献为22.71÷(22.71+55.20)×100%=29.15%，产流条件变化的相对贡献为55.20÷(22.71+55.20)×100%=70.85%。

4 结论

1) 在1956—1973年间的产流条件1下，降水变化使得河津站年均径流量平均减少约22.71%；在1974—2008年间的产流条件2下，降水变化使得河津站年均径流量在平、枯两种年份平均减少了约18.95%。

2) 在1956—1973年间的年降水为平水或枯水条件下(仅包括平、枯2种年份)，产流条件改变使得河津站年均径流量减少了约47.11%。在1974—2008年间的降水条件下(包括丰、平、枯3种年份)，产流条件改变使得河津站年均径流量减少了约55.20%。

3) 产流条件变化是汾河流域径流锐减的主要原因(70.85%)，同时，1974—2008年间丰水年份的大幅减少也是汾流径流锐减的重要影响因素(29.15%)。

本项研究的主要创新是明确提出了不同降雨条件下产流条件改变对河川径流的影响，以及降水丰枯变化在不同产流条件上对河川径流的影响，并通过统计相关分析和产流条件分析发现了汾河在丰水、平水年份的径流量主要与当年径流量和前一年径流量相关，枯水年份的径流量主要与前一年径流量相关。研究思路和方法可以供缺乏下垫面历史资料地区径流变化原因的定量评估所参考。