基于Copula分析生态恢复对黄土高原水沙频率耦合关系的影响

薛 帆，易海杰，宋松柏，张晓萍,

(1.西北农林科技大学水土保持研究所，黄土高原土壤侵蚀与旱地农业国家重点实验室，陕西 杨凌 712100；2.中国科学院教育部水土保持与生态环境研究中心，陕西 杨凌 712100；3.西北农林科技大学水利与建筑工程学院，陕西 杨凌 712100)

黄土高原自然环境复杂，地形破碎，土质疏松，植被破坏后极易发生水土流失，黄河年均输沙量曾达到16亿t，其中90%来源于黄土高原。自1970s，黄土高原开展了大规模的水土流失治理，1999年后又大力实施退耕还林(草)措施，土地利用和地表覆被发生剧烈变化，流域水沙数量、水沙关系成为研究热点。

水沙关系，其内涵是指某时空尺度上流域径流量与输沙或泥沙含量的关系，受自然因素，如降水、土壤、地形、植被等的影响，而大规模人类活动导致的下垫面条件变化，如大规模植被恢复工程、水土流失治理、农业灌溉等，对流域径流和泥沙过程均产生深刻影响。一方面表现在径流量和泥沙数量的协同增加或减少效应及密切程度；另一方面，还表现在其发生频率的协同性。目前，已有许多学者探究了黄土高原典型流域生态建设对水沙数量关系及其密切程度的调控效应，如刘宝元、Zhang等研究表明，生态环境治理下，黄河径流量由1919—1959年的年均426亿m减少到2010年以来的近300亿m，同时期输沙量由年均15.8亿t减少至不到2亿t。针对黄土高原流域水沙关系变化的研究，多是基于水、沙数量的角度进行理解与分析。Zhao等、Wang等研究表明，2000年以来黄土高原高含沙量明显降低；Zheng等研究表明，在流域尺度上，植被和坡面水保措施仅通过减水来减沙；Zhang等、Gao等基于经验回归关系式探究水沙关系的变化，2000年后水、沙数量相关性降低。而水沙关系在发生频率协同性上的研究，如郭爱军等基于Copula函数建立泾河流域水沙联合分布模型，结果表明，受生态建设影响，水沙丰枯各种遭遇组合频率分布更为均匀。以上研究主要从数量角度分析黄土高原典型流域水沙关系对生态恢复的响应特征，而对于其频率耦合关系的研究开展得还不够。

要深入剖析黄土高原流域水沙频率耦合关系的演变特征，需全面考虑径流、输沙2个相关变量的联合分布状况。Copula函数是定义域为[0，1]均匀分布的多维联合分布函数，它可以根据多个随机变量和其边缘分布函数的相关关系构建联合分布模型，具有计算方便、形式简单的优越性。目前，已经有学者将Copula函数应用于多变量水文分析计算的研究中，如水文随机变量概率分布、洪水频率分析、降雨频率分析、水文干旱特征分析、水文随机模拟等，均表明根据Copula理论建立的联合分布模型模拟精度较高。

黄土高原丘陵沟壑区，不仅是黄河下游多沙粗砂的主要策源地，也是黄土高原植被恢复最为显著的区域，开展该区生态环境变化下流域水沙特征演变规律具有很好的典型性。本文基于典型丘陵沟壑区北洛河上游刘家河水文站1960—2019年的日径流量和输沙量数据，采用Copula函数和经验模态分解(empirical mode decomposition)，建立水沙联合分布模型，分析多时间尺度下流域水沙频率耦合关系对水土保持与生态恢复的响应，为开展黄土高原水土流失治理规划、生态恢复等提供参考依据。

1 研究区概况

北洛河是渭河的一级支流，黄河的二级支流，流经陕西省榆林、延安、铜川、渭南及甘肃省庆阳等5个地(市)的18个县(区)，于大荔县东南注入渭河，地处34°42′—37°19′N，107°33′—110°10′E。刘家河水文站位于北洛河流域上游，其控制面积为7 325 km，约占流域总面积(26 905 km)的1/4，地貌类型为典型的黄土丘陵沟壑区，属暖温带半干旱气候，年平均降水量452 mm，多集中在7—9月；黄绵土为主要土壤类型，抗侵蚀能力差，极易发生严重的水土流失。自1970s以来，北洛河上游大规模实施生态建设和水土保持措施，如造林种草、修建淤地坝和梯田，2006底，上游地区各项水保措施面积达26.64万hm(表1)。1999年又推行退耕还林(草)政策，林草植被盖度由1970年的59.8%增加到2019年的81.1%。经过持续的生态建设，与1960s相比，刘家河水文站2010s径流量、输沙量分别减少约1/2和9/10，人类活动对径流量、输沙量变化贡献分别为94.7%，91.6%；不同发生频率的径流、输沙和含沙量存在不同程度的显著减少趋势，丰水期(1%～5%)的减少尤其显著。

表1 北洛河上游水土保持措施累积面积[24]

2 数据与方法

2.1 数据来源

本文收集了北洛河流域上游刘家河水文站1960—2019年实测日径流量、输沙量，并构建月径流量()和输沙量()序列进行研究分析，所用数据来源于黄河中游水文资料及年鉴。

2.2 数据处理方法

2.2.1 水沙数量时间变化趋势诊断 Pettitt非参数检验法广泛应用于水文序列突变检验。首先建立Mann-Whitney统计量(,)：

(1)

式中：和=2，…，。

令-=，则sgn()值为：

(2)

检验统计量()为：

=Max1≤≤|，|

(3)

统计量相关概率的显著性检验公式为：

(4)

2.2.2 二维Copula联合分布的建立

(1)确定边缘分布函数。通过皮尔逊III型曲线确定边缘函数：

(5)

式中：为皮尔逊III型曲线的形状参数；为尺度参数；为位置参数；Γ()为的伽马函数。

、和的矩法估计计算公式为：

(6)

(2)边缘分布函数检验。采用检验方法对边缘分布函数进行假设检验。设样本总体的分布函数为()且已知。将原假设为：总体的分布函数为(：，，…，)

检验统计量可表示为：

(7)

当为真时近似地有：

(8)

其拒绝域为：

≥(--1)

(9)

式中：为显著水平。

(3)Copula函数参数估计与选择。本研究选取Frank Copula、Gumbel-Hougaard(GH)Copula、Clayton Copula 3种函数构建二维水沙联合分布模型：

①Frank Copula

(10)

②Gumbel-Hougaard(GH)Copula

(11)

③Clayton Copula

(12)

式中：为Copula函数的参数。

为确定3种Copula函数参数，首先需得到径流量和输沙量的Kendall秩相关系数():

(13)

式中：、、、为观测数据；sgn为符号函数。

Copula函数的参数与τ存在的关系为：

①Frank Copula

(14)

②Gumbel-Hougaard(GH)Copula

(15)

③Clayton Copula

(16)

(4)Copula函数的拟合优度评价

本研究采用3种拟合优度评价指标：RMSE(均方根误差)、AIC(赤池信息量准则法)、BIC(贝叶斯信息准则法)，计算的评价指标值越小，则模型模拟越好。

①RMSE：

(17)

式中：(1,2,…,)为经验频率；(1,2,…,)为理论频率；为联合分布函数位数；为观测样本长度。

②AIC:

(18)

AIC=ln(MSE)+2

(19)

③BIC:

BIC=ln(MSE)+ln()

(20)

式中：为模型参数的个数。

(5)两变量联合概率与重现期。设边缘分布函数为()=，边缘分布函数为()=，则两变量联合概率分布为：

(,)=(<,<)=(,)

(21)

联合重现期为：

(22)

同现重现期为：

(23)

2.2.3 经验模态分解 经验模态分解EMD适用于非线性、非平稳性的信号处理，可将复杂的原始信号分解为有限个本征模函数IMF(intrinsic mode function)，分解出来的各IMF分量包含了原始信号在不同时间尺度的局部特征信号。

EMD经验模态分解的步骤为：

(1)输入需要处理的原始信号()；

(24)

(25)

(4)得到IMF1分量后，()减分量IMF1作为新的原始信号，进行(1～3)，得到IMF2分量，以此类推，最终完成分解。

EMD分解原始序列()还得到1个趋势项()：

(26)

式中：()为原始序列的各IMF分量。

3 结果与分析

3.1 水沙数量变化的阶段性

由图1可知，流域上游5，7月径流量均于1980年和2000年发生明显转折，且均达到极显著水平(<0.01)，而月输沙量突变时间与月径流量保持高度一致性。上游水沙3阶段的变化特性，与1980s以来大面积开展水土流失综合治理，以及1999年来退耕还林(草)措施实施的阶段性事件相吻合。因此，在分析中，将水文序列分为P1(1960—1979年)、P2(1980—1999年)和P3(2000—2019年)3个阶段。

注：Z表示不同月份P1(1960—1979年)、P2(1980—1999年)和P3(2000—2019年)的Mann-Kendall趋势检验统计量值；NS表示不显著。

3.2 水沙联合分布模型

假设P1、P2和P3的径流量()和输沙量()均服从皮尔逊III型分布，并采用检验对其进行假设检验，当显著水平为0.05时，自由度为2的临界值为5.99，而3个阶段、的边缘分布检验结果均小于该临界值，所以通过原检验。表明各阶段的、均服从皮尔逊III型分布(表2)。

表2 皮尔逊III型分布参数及检验结果

本文选用Frank、Gumbel-Houggard(GH)和Clayton 3种Copula函数分别构建刘家河站P1、P2和P3的水沙联合分布模型。从表3可以看出，P1、P2、P3阶段GH Copula函数的AIC、BIC和RMSE均小于Frank和Clayton Copula函数的计算结果。表明P1、P2和P3的水沙联合分布模型均采用GHCopula函数拟合精度最优。

表3 二维Copula函数参数估计

1(1960—1979年)：

2(1980—1999年)：

3(2000—2019年)：

由图2可知，P1、P2和P3径流量()和输沙量()的经验联合频率与GH Copula得到的理论联合频率拟合系数分别为0.999 0，0.998 1和0.993 5，表明各阶段所选的Copula函数与经验值拟合程度较好，选择合理。

图2 P1、P2和P3经验累积频率与理论累积频率的一致性比较

3.3 水沙丰枯遭遇频率分析

由表4可知，1960—2019年重现期为5～100年对应的径流量()和输沙量()均呈减少趋势。与P1相比，P2和P3各重现期的减小幅度分别为2.5%～26.9%和21.5%～55.7%，的减小幅度分别为31.4%～59.5%和72.7%～88.7%，的减小程度更为剧烈。造成上游和大幅降低的主要原因可以归结为2个方面：一是P2阶段的水土流失综合治理；二是P3退耕还林(草)工程实施以来林草植被覆盖率大幅提升(表1)。P2重现期为100年的径流量相较于P1表现为增加趋势，造成这一现象的原因可能是上游1994年遭遇了千年一遇的极端降水。

由表4还可看出，单变量的重现期大于两者的联合重现期(即发生大水或大沙事件)，小于两者的同现重现期(即大水大沙同时发生的事件)。如P1时段5年一遇的和的联合重现期为4.6年，同现重现期为8.1年，即较大的径流量可能伴随着较大的输沙量，也可能携带着较小的输沙量；径流量和输沙量同时超过实测最大值的事件发生概率减小。同时Coupla函数对3时期相应径流量和输沙量重现期的估算结果表明，随着时期进展，重现期相同的情况下，5～20年重现期的水沙事件，其联合重现期和同现重现期大多表现出小幅增加趋势，而重现期为50～100年的极端事件，在P1～P3阶段联合重现期和同现重现期均大幅增加。表明水土流失治理和生态环境建设等下垫面条件变化，对20年重现期以上的水沙极端事件发生概率的影响程度，要远大于20年重现期以内事件发生概率的影响。

表4 径流量(Q)与输沙量(QS)不同组合下的重现期

采用频率法将水沙频率分为丰、平、枯3种状态，本文以=25%，=75%作为水沙丰枯划分的频率，丰枯遭遇组合：

水丰沙丰：=(≥，≥)

水丰沙平：=(≥，<<)

水丰沙枯：=(≥，≤)

水平沙丰：=(<<，≥)

水平沙平：=(<<，<<)

水平沙枯：=(<<，≤)

水枯沙丰：=(≤，≥)

水枯沙平：=(≤，<<)

水枯沙枯：=(≤，≤)

从表5可以看出：(1)P1、P2、P3的水沙丰枯遭遇频率变化范围分别为0.6%～41.3%，1.1%～38.6%和0.6%～28.4%。(2)P1、P2、P3的水沙丰枯遭遇频率均表现为同步遭遇频率均大于异步遭遇频率；同步遭遇频率中，同平表现最大，同丰和同枯相近；异步遭遇频率中，水丰沙平与水平沙丰、水丰沙枯与水枯沙丰、水平沙枯与水枯沙平分别相等，水丰沙枯与水枯沙丰最小。表明北洛河上游水沙序列之间具有较强的相关性，出现极端情况的概率较小。(3)与P1相比，P2水沙丰枯各异步组合频率均增加；同步组合中同枯频率增加，而同丰和同平频率均减小。(4)与P1相比，P3水丰沙枯和水枯沙丰频率降低，其余异步组合频率均增大；各同步组合频率均减小。(5)总体来看，与P1相比，P2、P3水沙丰枯异步频率均增加，同步频率均降低，反映水土流失综合治理和生态建设等人类活动对流域水沙频率一致性的影响在逐渐扩大。

表5 不同时段水沙丰枯遭遇频率

3.4 水沙联合概率的多时间尺度特征

从表6和图3可以看出，北洛河上游水沙联合概率可分解为5或6个具有不同波动周期的本征模态分量和1个趋势项，表明北洛河上游存在复杂的水沙过程。P1具有平均周期为3.9，9.0，15.9，27.2，47.1，60.0月6个不同尺度的本征模态分量，P2具有平均周期为3.4，8.7，15.3，25.6，43.5，57.0月6个不同尺度的分量，P3则有5个本征模态分量，平均周期分别为3.3，8.1，13.8，24.8，39.7月。相同本征模态分量的平均波动周期，随时间推移均表现出减小的趋势。

注：图(a)为上游Copula水沙联合概率分布，图(b～e)为基于上游水沙联合分布概率，使用经验模态分解方法分解概率IMF1、IMF2、IMF3、IMF6/7；图(f)为经验模态分解概率以后的残差图(趋势项)。

表6 不同时段经验模态分解后各IMF分量的平均周期 单位：月

IMF1表明，无论是年、月尺度，P1水沙联合概率波动总体较稳定。与P1相比，P2水沙联合概率的年内波动频率略微增加，振幅呈减小趋势；然而1994—1996年波动频率明显减小，振幅显著增加，可能是受到1994年刘家河千年一遇暴雨的影响。P3水沙联合概率年内呈现2～4次波动；较大的波动出现在2000—2002年，之后波动逐渐放缓，趋于平稳。

IMF2表明，P1水沙联合概率年内基本为1次波动，年际间波动总体较稳定。P2水沙联合概率年内主要呈现1次波动，伴随着少数的2次年内波动，且振幅不明显。P3水沙联合概率年内为1～2次波动，主要表现在年内丰枯月份的交替，年际间波动较稳定。

IMF3表明，P1水沙联合概率主要存在1～3年的周期，到P1后期，周期明显变短，振幅显著增加。P2水沙联合概率主要表现为1～3年的周期，振幅总体较稳定。P3水沙联合概率初期波动较大，后期波动比较稳定。值得注意的是，如2012—2013年，振幅明显减小，年际间变化极不明显，可能是因为该时期径流量明显增大，但是输沙量却无显著变化，表现为水丰沙枯，且出现该现象的概率较小(表5)。

IMF5表明，在更大时间尺度上，水沙联合概率呈现3～6年的波动周期，且随着时间推移，周期逐渐变短。图3(f-1)表明上游P1水沙联合概率总体呈上升趋势，图3(f-2)、(f-3)表明P2、P3水沙联合概率总体呈下降趋势。

4 讨 论

滥砍滥伐和陡坡开垦等不合理的土地利用导致黄土高原生态环境恶化，侵蚀加剧。为有效治理水土流失，1970s以来，国家相继在黄河中游开展了大规模的水土流失治理和生态建设工程。一系列的生态恢复等人类活动使流域土壤侵蚀明显减弱，河流输沙量显著降低。与P1(1960—1979年)相比，P2(1980—1999年)和P3(2000—2019年)相同重现期的径流量和输沙量均明显减小，且输沙量的减少程度更加显著(表4)。表明生态恢复措施在相同重现期表现出减水更减沙的水土保持功能。

水土保持工程措施可以通过减小含沙量，进而实现水土流失的治理。北洛河上游1980s开展的水土流失综合治理，以修建淤地坝、梯田为主(表1)。其中，淤地坝的拦截效应使得排出淤地坝的洪水含沙量显著变小。与P1相比，P2大水伴随大沙的频率降低，水沙丰枯遭遇的异步频率略有增加；径流量和输沙量同丰、同平的频率明显降低(表4和表5)。而且同时期水沙联合概率的波动周期降低，振幅呈减小趋势(图3)。表明水土保持工程措施能够有效减小水沙联合概率和同现概率，降低高径流含沙量，调节水沙频率耦合关系。

研究表明，植被可以在一定程度上削减洪水的洪峰流量，有效降低降水及地表径流对土壤的侵蚀，改变水沙极端事件的发生概率。1999以来，黄土高原大力实施退耕还林(草)工程，北洛河上游植被盖度由2001年的30%增加到2019年的60%左右。与P1相比，2000年以来重现期为50～100年的水沙极端事件，联合重现期和同现重现期均大幅增加，水沙丰枯遭遇的同步组合频率均明显降低。然而，重现期为5年的水沙联合重现期减小，水丰沙平和水平沙枯的频率增加，这是由于植被能显著增加河流枯季径流量。林草植被盖度在小于40%～50%时，减水减沙、调节水沙极端事件发生概率的效果非常明显，但植被盖度达到60%左右时，产流、产沙系数均趋于稳定。P3水沙联合概率初期波动较大，后期无论是波动周期、振幅，均比较稳定(图3)。表明P3退耕还林(草)措施显著影响水沙联合概率，尤其是在2000年初期，而后期水沙联合概率逐渐稳定，这可能与植被盖度持续增加有关。

流域水沙频率耦合关系的变化，不仅受剧烈的人类活动影响，还与极端天气的发生息息相关。1994年刘家河发生了千年一遇的暴雨，使得P2重现期为100年时的径流量相比于P1表现为增加趋势(表4)，而且导致1994—1996年水沙联合概率的波动周期和振幅明显增加(图3)。说明极端暴雨对流域水沙联合概率有显著的影响，其内在影响机制还需进一步探讨。

5 结 论

(1)Gumbel-Houggard Copula能很好地模拟上游水沙联合概率分布。

(2)随时段推移，5～100年重现期的径流量、输沙量均明显减小，各设计频率的水沙数量在发生概率上的异步性增强。

(3)水土流失治理和生态环境建设等下垫面条件变化，大大降低20年以上重现期水沙极端事件发生的概率。

(4)水沙联合概率在不同时间尺度上具有非常复杂的变化特征。

通过Copula函数建立的水沙联合分布模型，可以进一步了解生态恢复对黄土高原水沙频率耦合关系的影响，为黄土高原的水沙调控提供参考依据。