气候因子对地表水资源量变化影响的定量分析

陈立华，王　焰，关昊鹏

(1.广西大学土木建筑工程学院，南宁 530004；2.广西防灾减灾与工程安全重点实验室，南宁 530004)

0　引　言

2017年《北部湾城市群发展规划》中钦州市作为“双轴”城市，区域内主要为源短流急中小入海河流，汛期降雨量占全年80%以上，全市各类水库共393座，总兴利库容仅为4.85 亿m3[1]。2012年《中国近海海洋调查评价》指出近90%的城市存在不同程度的缺水问题。在区域人口增长、经济发展与水资源短缺间矛盾日益凸显境况下，合理计算滨海区地表水资源量及其未来趋势变化，剖析在降雨、蒸发、气温及其他因子影响下，各因子对径流变化的贡献率大小，是分析水资源潜力和正确评价与配置的前提。

深入分析地表水资源量在各因子影响下趋势变化首要问题是径流还原。因滨海区不同于内陆，具有河网密布且水文站点不足等特点，导致计算难度偏大。因此在满足其精度条件下，有必要探究其他径流还原方法在滨海地区的适用性。第二次全国水资源综合规划的水资源调查评价中规定选用的径流还原方法为：降雨径流变化趋势法、双累积曲线(Double Mass Curve，简称DMC)模型、蒸发差值法、分项调查法等，其中DMC模型[2]和蒸发差值法[3]相比于其他方法，具有原理简单、操作简便、资料要求低、忽略要素少等特点，因此采用以上两种方法对径流进行还原计算，并结合钦州市已有简化分项调查法对其结果合理性、方法适用性方面进行评价。关于环境变化下各因子对径流影响的贡献率分析，目前在流域尺度上，可定量分析气候变化与其他因子对水文影响的研究方法主要有水文模型法[4-6]和定量评估法。其中定量评估方法主要有气候弹性系数[7-8]、敏感性分析法[9]、降雨-径流双累积曲线法等。水文模型法具有较好的物理基础，但参数敏感性存在一定的不确定性，若对模拟结果不进行验证，极可能使气候变化对径流的贡献值偏高[10]；尽管定量评估法所需数据较少，但所需较长数据序列的同时，长序列中噪音会对评估结果造成干扰[8]。王随继[11]等于2012年提出累积量斜率变化率比较法(Slope Changing Ratio of Cumulative Quantity，简称SCRCQ)法，可有效剔除噪音，较简便地分离出气候变化和其他因子对径流变化的影响比重，在洞庭湖三口河系[12]、长江荆南三口[13]、南水北调中线工程[14]、湟水川流域[15]、黑河流域中上游[16]等地区均得到较好的应用，2017年Wu[17]等人更是将其应用于岩溶区域的气候变化与人类活动研究。但上述针对SCRCQ法的研究，所假定的基准期均为实测径流序列，未对其进行还原计算以消除早期人类活动影响，未充分考虑气温对降雨与蒸发的影响，未深入探讨降雨、蒸发、气温对径流变化影响复杂的内在联系，这将导致气候因子贡献比重偏小，结果误差偏大。针对研究中存在的问题，有必要对SCRCQ法进行适当改进：基于SCRCQ法的基础上重新定义基准期，考虑气温对降雨与蒸发的影响，并结合水量平衡原理探讨各气候因子对径流变化影响的内在联系，并间接获取其他因子贡献率值，使SCRCQ法分析结果更贴近实际情况。

1　基础数据

钦州市地处桂南滨海，属水资源四级区，境内河网密布，多为中小入海河流。将研究区域划分为19个子流域如图1所示，钦州市境内河流四级分区如表1所示。研究所选用资料包括32个雨量站、5个水文站点、3个蒸发站(陆屋、黄屋屯、坡朗坪站)以及1个气象站(钦州站)分布如图1所示，研究所采用的降雨、径流、蒸发、气象资料统一采用1956-2013年月平均序列。钦江为钦州市最大流域，如图2所示，是钦南、钦北区及灵山县最重要的取水来源河流。钦江流域干流全长191 km，集水面积2 326 km2，干流坡降0.36%，地势东北部高、西南部低(钦江下游滨海平原)。

图1　钦州市子流域及站网图Fig.1 Station network and sub basin in Qinzhou

图2　钦江流域地理位置Fig.2 geographical location Qinjiangbasin

2　分析方法

2.1　气候变化对水文影响贡献率分析方法

在SCRCQ法基础上重新定义基准期，为实测径流序列首个突变点前的还原径流序列，研究期为突变点后实测径流量序列。基准、研究期累积降水量～年份线性方程的斜率分别为Kpa(mm/a)、Kpb(mm/a)，斜率变化率为Sp；累积蒸散发量～年份线性方程的斜率分别为Kra(mm/a)、KEb(mm/a)，斜率变化率为SE；累积径流量～年份线性方程的斜率分别为KRa(亿m3/a)、KRb(亿m3/a)，斜率变化率为SR；累积气温～年份线性方程的斜率分别为KTa(℃/a)、KTb(℃/a)，斜率变化率为ST。

SP=(KPb-KPa)/KPa

(1)

SE=(KEb-KEa)/KEa

(2)

SR=(KRb-KRa)/KRa

(3)

ST=(KTb-KTa)/KTa

(4)

根据式(1)～式(4)，可求降雨、蒸发、气温对径流变化的贡献率，分别为CP、CE、CT：

CP=100×SP/SR

(5)

CE=100×SE/SR

(6)

CT=100×ST/SR

(7)

依据所得CP、CE、CT，我们可间接获取其他因子对径流变化的贡献率CO。根据水量平衡原理如式(9)，可知T为自变量，R为因变量，T能够直接引起R变化，贡献率为CT；不仅如此，T的改变可引起P、E变化，间接影响到R，贡献率同样为CT，但该部分CT对CP、CE进行了重复计算。因此，CO并非定值，而应处于是否考虑T影响下的值域即上下限内变化，如式(8)所示。

CO∈[100-CP-CE-CT, 100-CP-CE]

(8)

2.2　突变点检测方法

DMC模型与SCRCQ法均需确定基准期，其关键在于分析水文序列的突变点位置。水文序列突变检测方法众多，其中Mann-Kendall法[1]与累积距平法[20]原理简单，M-K法毋需要素样本遵从特定分布，不受异常值的干扰，检测范围宽以及理论基础较强等特点获得广泛应用；累积距平法的核心是判断离散数据对其均值的离散幅度，若累积距平增值大，表明离散数据大于其平均值，反之则小于其平均值[11]。但以上两种方法检测突变点存在检测结果不一致等问题，因此辅以滑动T检验法[21]检测时间序列突变点的显著性，若检测超过|t|≥tα=0.05=2.003显著性水平的时间节点即为突变点，否则为转折点。

2.3　地表水资源总量计算方法

DMC模型原理如图3所示，若径流序列在1986年发生突变，则1986年前序列为基准期①(实心方框)；1986年后为还原期②(空心圆框)。若处于天然条件下则k2=k1，而图中k1>k2，表明1986年后径流序列受外界环境改变影响，致使k1减小，因此以①期降雨-径流累积量线性关系式替代②期进行还原计算。

图3　DMC模型原理图Fig.3 Schematic diagram of DMC model

流域蒸发差值法遵循水量平衡原理[18]，地表水资源量R变化与降雨P、流域下垫面变化前的流域蒸发E前、土壤蓄水量W紧密相关，W又可表示为气温T与P的函数如式(9)。

R=P-E前-W=f(P,E前,T)

(9)

其中E前可表示为：

E前=[(A-A1)E陆+αA1E水]/A

(10)

式中：A为设计流域总面积，km2；A1为人类活动前流域内的水库、湖、塘等水面面积，km2；α为水面蒸发折算系数；E水为蒸发皿实测水面蒸发，mm；E陆为陆面蒸发。

研究区域属南方滨海区，因此陆面蒸发参考凯江公式[19]：

E陆=(αT+β)θ/L

(11)

θ=θ0[0.65-0.35(n-s)]

(12)

式中：θ为太阳总辐射，cal/cm2；T为日平均气温，℃；L为蒸发潜热；α和β为系数和常数；θ0为碧空无云时的太阳总辐射，cal/cm2；n为月平均云量；s为月平均日照量。

3　应用实例分析

3.1　基准期分析

对钦州市内有水文控制站流域径流量序列采用累积距平法及M-K法进行突变点检测，主要以钦江流域为例进行详细介绍。

由图4(a)可知，UF与UB曲线相交于1986、1993、1994、1996、1997、2001、2003年，累积距平曲线在1964、1986、1992、2003年为极值点，在以上时间节点处径流量极可能发生突变，但M-K法所检测出突变点均在α=0.05(±1.96)显著性水平，累积距平法显著性无法判断，因此需辅以滑动T对其进一步显著性水平检验。由滑动T结果可知，仅有1986年与2003年t值分别为2.47和2.55，均超过95%显著性水平，表明1986、2003年为陆屋站径流量序列突变年份；对于其他节点并未超过95%显著性水平，仅为转折点。采用滑动T对图4(b)、图4(c)以及图4(d)中降雨量、蒸发量、气温序列可能的突变节点进行显著性分析，降雨量序列在整个时间域内并未出现显著突变点；蒸发量、气温序列均仅有一个突变点超过95%显著性水平，分别为1973年和1996年。

经分析，黄屋屯、坡朗坪、西塘、大江口站径流量均于1986年发生显著突变，而面降雨量序列均不存在明显突变点，因此将1956-1986年假定为DMC模型径流还原基准期，对1987-2013年序列进行径流还原并计算钦州市地表水资源总量。

3.2　地表水资源总量及其趋势分析

3.2.1钦州市地表水资源总量

DMC模型以1956-1986年为基准期，经Pearson相关分析及T统计量相关显著性分析可知，基准期内5个水文站降雨～径流相关系数r均超过99%显著性水平如表1所示，可对1987-2013年径流序列进行还原计算。蒸发差值法则对1956-2013年整个时段在降雨量的基础上逐年剔除流域年蒸发量，从而得到水文站天然径流量。在水文站天然径流量基础上采用代表站法分析不同保证率下各流域地表水资源量，结果如表1所示。DMC模型与蒸发差值法所得钦州市地表水资源总量分别为106.573、111.533 亿m3。

表1　各典型流域水文控制站降雨～径流Pearson相关分析Tab.1 Analysis of Pearson correlation between rainfall and runoff in typical watershed hydrological control station

图4　时间序列M-K法及累积距平法突变点检测Fig.4 Time series M-K method and cumulative anomaly detection

根据《广西钦州市水资源综合规划报告(2008年)》所采用《广西壮族自治区地表水资源(1984年)》中1956-1980年实测～天然径流相关系数，对1980-2000年实测径流进行还原，所得地表水资源量为104.41 亿m3。蒸发差值法严格遵循水量平衡原理，对整个序列均进行还原计算，从而消除了早期人类活动对径流所产生的影响，使地表水资源量更为接近于天然情况。综上，表明两种径流还原方法在钦州市内不仅能够很好的应用，而且计算简单、所需资料较少、结果与天然情况较为接近，在南方滨海地区水资源总量计算中值得推广。

从表2中可知，钦州市地表水资源量最为丰富的为独流入海河流水资源分区，其地表水资源量占全市总量62.7%～63.3%。由于独流入海河流源短流急，实际可利用水量不多，如钦江流域多年平均入海水量为17.40 亿m3，其实际可利用水量仅为3.082～4.601 亿m3。

3.2.2径流序列的趋势变化分析

采用3 a滑动平均与Morlet小波综合对钦州市地表水资源量进行分析如图5所示，58年序列中存在1960-1973年、1974-1987年、1988-2002年3个平均周期为14 a的完整丰-枯变化过程，1986年开始，降雨、径流波动幅度较之前明显增大。尽管整个水资源序列呈0.087 亿m3/a的上升趋势，但自20世纪70年代以来序列却呈-0.400 亿m3/a的下降趋势。进一步对不同年代均值进行比较可知，各年代均值呈丰-枯转换阶梯式下降趋势，降雨量序列趋势性与其相同。综上所述，可推测地表水资源量在2017-2023年左右将处于偏枯期，2024-2030年左右将处于偏丰期；根据天然与实测径流量差值曲线如图6可知，自20世纪60年代以来，各年代天然与实测径流量差值序列趋向率高达0.406 亿m3/a，总体呈明显梯级上升趋势，表明以人类活动为主的其他因子对实测径流减少的影响持续增大，因此有必要选取钦州市典型流域，进一步定量分析降雨、蒸发、气温及其他因子对径流减少影响的贡献率研究。

图5　钦州市水资源总量趋势分析Fig.5 Trend analysis of surfacewater resources in Qinzhou

表2　钦州市各水资源分区地表水资源量Tab.2 Surface water resources quantity of each water resources division in Qinzhou

图6　天然与实测径流量差值曲线趋势分析Fig.6 Analysis of difference between natural and measured runoff

3.3　气候变化对主要河流径流减少的贡献的定量分析

采用改进SCRCQ法对钦州市主要河流定量分析降雨、蒸发、气温变化对径流减少影响的贡献率的同时，间接分析其他因子对径流减少影响的贡献率，其中以钦江流域为例进行详细介绍，基准期天然径流量序列则采用蒸发差值法还原结果。

钦江流域各因子突变年份如图4所示，径流序列(1986、2003)，降雨量序列无显著突变点(为便于计算，假定其突变点与径流序列一致)，蒸发量序列(1973、2003年)，气温序列(1996年)。若仅考虑径流序列突变点情况，将1956-2013年序列统一划分为3个时期(1956-1986年为AR期，即基准期；1987-2003年为BR期，2004-2013年为CR期，BR与CR期为研究期)，则将使结果偏差较大。因此根据突变点情况，蒸发量、气温序列的时期划分需额外考虑：蒸发量序列AR期仅为1956-1973年，气温序列BR期仅为1997-2003年，如图7所示。

BR与AR期相比，累积径流深线性关系式斜率减少229.39 mm/a，减少率为23.04%；同期相比累积降雨量线性关系式斜率增加48.11 mm/a，增加率为2.73%；累积蒸发量线性关系式斜率减少74.73 mm/a，减少率为7.56%；累积气温线性关系式斜率增加0.97 ℃/a，增加率为4.41%。根据式(9)水量平衡可知地表水资源量变化是降雨量、蒸发量、气温的函数，均为气候因子，而实测径流量则存在其他因子影响，由其他因子影响较弱的AR期的累积径流量与年份之关系(图7)可知，假设BR期不存在其他因素影响，则上述气候因子斜率变化率之和与径流斜率变化率相等，而实际上径流深变化是气候及其他因子综合作用的结果。根据式(1)～(7)计算结果可知，BR与AR期相比，降水、蒸发以及气温对径流深减少的贡献率为11.84%、32.80%和19.14%，由式(8)可知其他因子贡献率并非定值，而应处于是否考虑气温影响下的值域，即[36.21%,55.35%]的上下限内变化。若依据原方法，在基准期仍为实测径流序列、不考虑气温因素的影响，则其他因子对径流减少的贡献率为28.24%，气候因子贡献率明显偏大。同理可知，CR与AR期相比，降水量、蒸发量的减少以及气温的升高对径流减少的贡献率分别为19.02%、5.32%和10.65%，其他因子对径流减少的贡献率为65.01%～75.66%的值域内，钦州市主要河流气候因子对径流减少的贡献率如表4所示。

图7　累积径流深、降雨量、蒸发量、气温与年份之关系Fig.7 Relationship between cumulative runoff depth, rainfall, evaporation, temperature and year

表3　钦江流域气候因子对径流变化的贡献率分析Tab.3 Analysis on contribution rate of climatic factors to runoff variation in Qinjiang Basin

表4　钦州市主要河流气候因子对径流变化的贡献率　%Tab.4 The contribution of climatic factors to runoff change in Qinzhou

由表4可知，BR与AR期相比，位于钦州西南、西、西北部区域的钦江、茅岭江，径流减少的贡献因子主要以其他因子与蒸发为主，降雨、气温的影响贡献率相对较小，而钦州东、北、南以及中部区域的大风江、小江、武思江，径流减少的贡献因子则以蒸发为主导，降雨与其他因子贡献率次之，不同流域受其他因子的贡献率排序：钦江>茅岭江>小江>武思江>大风江；CR与AR期相比，各流域径流减少的主导贡献为其他因子，具体为：其他因子>降雨>蒸发>气温，不同流域受其他因子的贡献率为：钦江>茅岭江>小江>大风江>武思江。

综上所述，各流域径流减少的主导贡献因子由蒸发向其他因子转移，以人类活动为主的其他因子从BR期的钦江流域逐步扩展至CR期的钦州市全境，而降雨、气温因子对径流减少影响的贡献率在整个研究阶段(1956-2013年)内无明显变化，且居于次要地位。

4　结　语

(1)为确定钦州市主要流域径流序列基准期，采用M-K、累积距平法并辅以滑动T检测各时间序列突变点。钦江、茅岭江、大风江径流量序列存在1986、2003年两个突变点显著，小江、武思江流域为1986年；各流域面雨量序列无明显突变点，陆屋站蒸发量序列为1973年突变显著，面蒸发量序列在此基础上增加一个突变点于2003年。相对于上述突变点而言其他可能的突变点突变特征不显著，仅在特定年份发生转折或处于丰枯变化节点上，可视为转折点。

(2)采用DMC模型及蒸发差值两种简化方法进行径流还原，径流还原结果偏差控制在9.2%以内，钦州市地表水资源量处于106.573～111.533 亿m3范围内，相对偏差仅为4.4%；自20世纪70年代以来地表水资源总量总体呈-0.400 亿m3/a的趋势下降，而天然与实测径流量差值曲线呈明显阶梯上升趋势；地表水资源量存在14 a左右的丰-枯转换周期，且在2017和2024年左右将分别步入偏枯、偏丰阶段。

(3)采用SCRCQ法定量分析了钦州市主要流域地表水资源量受各因子影响下径流减少的贡献率，其中钦江流域BR与AR期相比，降雨、蒸发和气温对径流减少的贡献率分别为11.84%、32.80%、19.14%，其他因子的贡献率为36.21%～55.35%范围内；CR与AR时期相比，降雨、蒸发和气温对径流减少的贡献率为19.02%、5.32%、10.65%，其他因子的贡献率有所增加为65.01%～75.66%。钦州市主要流域各因子对径流减少的贡献研究表明，1986-2003年导致径流减少的主要贡献为气候变化，其中又以蒸发因子最为显著，而在2004-2013年间，其他因子对径流减少的贡献明显高于气候变化。

□

参考文献：

[1]陈立华,王焰,易凯,等.钦州市降雨及入海河流径流演变规律与趋势分析[J].水文,2016,36(6):89-96.

[2]申红彬,吴保生,郑珊,等.黄河上游兰州至头道拐河段流量频率变化分析[J].水力发电学报,2014,33(1):7-14.

[3]魏茹生.径流还原计算技术方法及其应用研究[D].西安:西安理工大学,2008.

[4]王国庆,张建云,贺瑞敏.环境变化对黄河中游汾河径流情势的影响研究[J].水科学进展, 2006,17(6):853-858.

[5]夏智宏,刘敏,王苗.1990s以来气候变化和人类活动对洪湖流域径流影响的定量辨识[J].湖泊科学,2014,26(1):515-521.

[6]贺瑞敏,张建云,鲍振鑫.海河流域河川径流对气候变化的响应机理[J].水科学进展,2015,26(1):1-9.

[7]Fu Guobin, Charles Stephen P, Chiew Francis H S.A two-parameter climate elasticity of streamflow index to assessclimate change effects on annual streamflow[J].Water Resources Research, 2007,43(11):2 578-2 584.

[8]Sankarasubramanian A, Vogel R M, Limbrunner J F.Climate elasticity of streamflow in the United States[J].Water Resources Research, 2001,37(6):1 771-1 781.

[9]Jones R N, Chiew F H S, Boughton W C, et al.Estimating the sensitivity of mean annual runoff to climate change using selected hydrological models[J].Advances in Water Resources,2006,29(10):1 419.

[10]Legesse D，Vallet C C, Gasse F.Hydrological response of acatchment to climate and land use changes in Tropical Africa: Case study South Central Ethiopia[J].Journal of Hydrology,2003, 275(1-2):67-85.

[11]王随继,闫云霞,颜明,等.皇甫川流域降水和人类活动对径流量变化的贡献率分析-累积量斜率变化率比较方法的提出及应用[J].地理学报,2012,67(3):388-397.

[12]李杨,杨文尧.洞庭湖三口河系地区径流演变情势与农业用水面临的挑战[J].中国农村水利水电,2016,(11):86-92.

[13]帅红,李景保,何霞,等.环境变化下长江荆南三口径流变化特征检测与归因分析[J].水土保持学报,2016,30(1):83-88.

[14]李凌程,张利平,夏军,等.气候变化和人类活动对南水北调中线工程典型流域径流影响的定量评估[J].气候变化研究进展,2014,10(2):118-126.

[15]张调风,朱西德,王永剑,等.气候变化和人类活动对湟水河流域径流量影响的定量评估[J].资源科学,2014,36(11):2 256-2 262.

[16]何旭强,张勃,孙力炜,等.气候变化和人类活动对黑河上中游径流量变化的贡献率[J].生态学杂志,2012,31(11):2 884-2 890.

[17]WuL H, Wang S J, Bai X Y, et al.Quantitative assessment of the impacts of climate change and humanactivities on runoff change in a typical karstwatershed, SWChina[J].Science of the Total Environment,2017, 601:1 449-1 465.

[18]李林,汪青春,张国胜,等.黄河上游气候变化对地表水的影响[J].地理学报,2004,59(5):716- 722.

[19]金栋梁,杨世才.用普通气象资料计算土壤蒸发量的方法(及凯江蒸发公式)[J].人民长江,1981,(4):47-52.

[20]Ran L S, Wang S J, Fan X L.Channel change at Toudaoguai station and itsresponses to the operation of upstream reservoirs in the upper Yellow River[J].Journal of Geographical Sciences, 2010,20(2):231-247.

[21]余予,孟晓艳,张欣.1980-2011年北京城区能见度变化趋势及突变分析[J].环境科学研究,2013,26(2):129-136.