大沽河水沙演化规律及驱动因素研究

胡萌　宋委　刘克金　吴振　雷炳霄　王毅

［关键词］Mann-Kendall检验；突变检验；小波分析；气候变化；径流；泥沙；大沽河

［摘要］基于1956—2016年长系列水沙资料，利用MK、滑动T检验、累积距平曲线、小波分析等方法分析了大沽河水沙变化规律，在此基础上，计算了气候变化和人类活动对径流量、输沙量变化影响的贡献率。结果表明：大沽河径流量、输沙量呈显著减少趋势，降水量呈不显著减少趋势；径流量在1976年、输沙量在1965年和1979年均发生了由多到少的突变；天然径流量存在8 a、12 a、21 a三个变化周期，其中21 a是第一主周期；输沙量周期变化存在5～8 a、9～17 a和18～30 a三类较明显的时间尺度，变化的第一主周期是12 a；人类活动是造成径流量、输沙量减少的主要原因，人类活动对径流量变化的贡献率为76.8%～97.4%，对输沙量变化的贡献率为76.2%～98.4%。

［中图分类号］ TV213.9［文献标识码］ A［文章编号］ 1000－0941（2023）08－0032－06

流域系统是对气候变化与人类活动异常敏感的动态系统，水沙演化是其中最为活跃的部分[1]。径流量和输沙量的变化不仅关系到河流自身的发展演变，也反映了流域的环境状况、水土流失程度及人類活动的影响[2-3]。

已有的研究主要关注大型流域系统水沙运移规律及其影响因素，而中小型流域系统的重要性往往被忽视或低估［4-6］。大沽河是胶东半岛最大的河流，地理位置与气候环境独特，流域面积占青岛市域总面积的42%，被誉为青岛市的“母亲河”。

目前，对大沽河水沙研究多集中于单一要素方面［7-10］，但尚未见针对大沽河流域开展气候变化与人类活动双重驱动下水沙耦合演化方面的研究。本研究以大沽河流域降水、径流、泥沙为研究对象，分析水沙演化规律，定量研究气候变化和人类活动对径流量和输沙量变化影响的贡献率，以期为流域环境保护与生态文明建设提供参考。

1研究范围与方法

1.1研究区域及数据

大沽河地处山东省胶东半岛西部，位于120°07′～120°34′E、36°02′～37°05′N之间。河流发源于招远市东北部的阜山，自北往南流经招远市、莱西市、平度市、即墨区、胶州市、城阳区6个市（区），于胶州市营海街道办事处东营村入胶州湾，干流全长199 km，流域面积6 205 km2。

选取大沽河干流南村水文站1956—2016年逐月降水量、实测径流量、实测输沙量资料，采用青岛市第三次水资源调查评价成果，对南村水文站控制流域内受水利工程和其他因素的影响而损耗或增加的河川径流量逐年进行还原计算，求得历年天然月、年径流量。该数据系列完整，无缺失。南村水文站位于青岛市平度市南村镇南村，控制流域面积3 724 km2，距入海口62 km，在大沽河流域具有较好的代表性。

1.2研究方法

1.2.1变化趋势研究

本研究采用线性趋势法[11]、Mann-Kendall秩相关检验法[12]（以下简称“MK法”）和改进MK法（包括MMK法[13-14]、PW-MK法[15]、TFPW-MK法[16]）进行趋势分析。

MK法是世界气象组织推荐并已被广泛使用的一种非参数检验方法，该方法的优点是不需要待检序列遵从一定的分布。构建检验统计量Z，当Z>0时，存在上升的趋势；当Z<0时，存在下降的趋势。当样本数量n比较大时，Z近似服从标准化正态分布，给定显著性水平α，则可根据Z与临界值Z（1-α/2）的比较结果，判定序列趋势的统计显著性。

MK法基于序列独立性假设，而水文序列的自相关性将影响检验结果的显著性水平。改进的MMK法、PW-MK法和TFPW-MK法均对时间序列的相关性进行了处理，减小自相关性对检验结果的影响，详见表1。

1.2.2突变检验

本研究采用滑动T检验法、累积距平曲线法进行突变检验。

滑动T检验是通过考察两个样本平均值的差异是否显著来检验突变点。对于具有n个样本的时间序列x，人为设置某一时刻为基准点，基准点前后两端子序列x1和x2的均值差异如果超过了一定的显著性水平，可以认为均值发生了质变，有突变发生[17]。该方法的缺点是子序列时段的选择具有人为性。为避免任意选择子序列长度造成突变点的漂移，具体使用这一方法时，可以反复变动子序列长度进行比较，提高计算结果的可靠性[18]。

累积距平曲线法是一种直观判断趋势变化的方法，核心是离散数据大于平均值，累积距平值增大，曲线呈现上升趋势，反之则呈现下降趋势；根据累积距平值曲线的波动起伏，可以判断长期的演变趋势以及变化趋势发生突变的时间[19]。

1.2.3周期分析

本研究应用Morlet连续复小波分析法[20-21]研究径流量和输沙量的周期性变化规律，提取变化的主周期，揭示多时间尺度变化的复杂结构。小波分析是一种窗口大小固定但时间窗和频率窗都可改变的时频局部化方法，可以同时揭示时频细节信息。小波方差随尺度的变化过程称小波方差图，它反映了波动的能量随尺度的分布情况。通过小波方差图，可以确定一个水文序列中存在的主要时间尺度，即主周期。

1.2.4影响因素的贡献率计算

本研究采用双累积曲线法[22]分析不同时期气候和下垫面变化对径流量影响的贡献率，采用累积量斜率变化率法[23-24]分析气候变化和人类活动对输沙量影响的贡献率。

1）双累积曲线法。通过突变检验识别径流量时间序列的突变点，以径流量变异点前作为受人类活动影响较小的基础期，采用双累积曲线建立径流量变异前累积天然径流量与累积降水量关系式，并用人类活动影响显著时期（变异期）累积降水量计算得到变异期的模拟累积天然径流量，认为模拟径流量与基准期人类活动影响条件近似相同，对比变异期与基础期的实测径流量、天然径流量、模拟径流量等变化，可分析不同时期降水量和人类活动对径流量影响的贡献率。

2）累积量斜率变化率法。设累积输沙量—年份线性关系式的斜率在变异点前后两个时期分别为SR1和SR2，则累积输沙量斜率变化率计算公式[24]为

RS=（SR2-SR1）/SR1×100%（1）

RS为正表示斜率增大，为负表示斜率减少。同理，求得累积降水量斜率变化率PS。

降水量变化对输沙量变化的贡献率Cp的计算公式为

CP=PS/RS×100%（2）

人类活动对输沙量变化的贡献率CH计算公式为

CH=1-CP（3）

2水沙演化特征

2.1变化趋势

大沽河南村水文站降水量、径流量与输沙量年际变化线性趋势分析见图1。由图1可知，1956—2016年大沽河降水量、实测径流量、天然径流量、输沙量均呈下降趋势，平均减少速率分别为2.146 3 mm/a、0.118 3亿m3/a、0.123 1亿m3/a、1.927 3万t/a；其中径流量、输沙量均达到显著性水平（p＜0.05），而降水量未通过显著性检验。

人类活动对大沽河年输沙量的影响主要表现在修建水利工程和进行水土流失治理。大沽河流域共12座大中型水库，除北墅水库（1974年）、庙埠河水库（1971年）外，其余10座均建成于1958—1960年。此后，20世紀60年代也是小型水库、塘坝建设的高峰期。由图1（b）可知，大沽河输沙量在1965年以后急剧减少。分布在流域范围内的水库、塘坝一方面拦蓄了坡面汇入河道的径流和泥沙，导致大沽河干流特别是中下游输沙量大幅减少；另一方面，河道内的径流被拦蓄后，上游来水失去了冲刷下游河道的动力，从而达到减蚀的效果，这在很大程度上改变了流域的水沙状况。

采用MK法、MMK法、PW-MK法、TFPW-MK法进行趋势检验（见表2），检验统计量Z均小于0，说明4个序列都呈下降趋势；实测径流量、输沙量的MK、MMK统计量与天然径流量的MK、MMK、PW-MK、TFPW-MK统计量|Z|>1.96，下降趋势显著（p＜0.05），其他项目均未通过显著性检验。

2.2突变分析

图2为大沽河南村水文站径流量、输沙量累积距平曲线。由图2可知，1956—1958年径流量累积距平值呈下降趋势，1958—1965年出现明显的上升趋势，1965—1976年小幅波动，1976—2016年呈现逐渐减少趋势。由此可见，实测径流量、天然径流量的突变点均为1976年，均发生了由多到少的突变。而输沙量累积距平曲线大致表现为3个变化阶段：1956—1965年呈明显上升趋势，1965—1979年缓慢下降，1979—2016年呈线性减少趋势。综合分析，输沙量的突变点为1965年、1979年，均发生了由多到少的突变。

采用滑动T检验法对径流量、输沙量进行突变分析，步长k分别取2，3，…，20，所得结果见图3。由图3可知，k≥7时，实测径流量通过α=0.01显著性检验，

有一个突变点，突变年份为1976年（图3取k=7为代表曲线）；k≥16时，天然径流量通过α=0.01显著性检验，有一个突变点，突变年份为1976年（取k=18作为代表曲线）；k≥8时，输沙量通过α=0.01显著性检验，有两个突变点，突变年份分别为1965年、1979年（取k=9作为代表曲线）。降水量未检出突变点。

综上可知，累积距平法与滑动T检验法的结论基本一致，实测径流量、天然径流量突变年份为1976年左右，输沙量突变年份为1965年左右、1979年左右，均发生了由多到少的突变。径流量、输沙量突变年份前、后10 a均值及其变幅见表3。

2.3变化周期

2.3.1天然径流量周期及震荡情况分析

图4（a）为天然径流量小波系数实部图，正值表示丰水阶段，负值表示枯水阶段。1956—2016年间，天然径流量在年际和年代际上都存在明显周期变化，包括5～9 a和10～15 a的小尺度信号以及16～30 a的大尺度信号。其中，在5～9 a时间尺度有周期表现，但不明显，其尺度中心在8 a左右；10～15 a时间尺度上正负相位明显且较稳定，其尺度中心在12 a左右；16～30 a时间尺度内小波曲线闭合完整，正负相位交替出现，说明其周期表现明显，其尺度中心在21 a左右。从图4（a）还可以看出，以上3个尺度的周期变化在整个分析时段前期（1976年以前）较为稳定，1976年以后天然径流量年际间变幅较小，震荡周期表现不够稳定。1976年之前小波系数实部等值线密集，1976年之后变稀，说明进入20世纪80年代以来，大沽河径流量变幅趋缓，处于低流量期，与前述突变检验的结论一致。

图4（b）为天然径流量小波方差图，图中共有3个峰值，第一峰值为21 a，该时间尺度波动能量最强，是丰枯变化的主周期；其次还存在12 a和8 a两个峰值，分别对应着第二、第三主周期。上述3个周期共同对天然径流量变化起作用，3个周期的波动控制着天然径流量在整个时间域内的变化特征。

2.3.2输沙量周期及震荡情况分析

图5为大沽河南村水文站输沙量小波系数实部、小波方差图。由图5可知，输沙量呈明显的周期性变化，包括5～8 a和9～17 a的小尺度信号以及18～30 a的大尺度信号。以上3个尺度的周期变化在整个分析时段前期（1979年以前）表现比较稳定。1979年之前小波系数实部等值线密集，1979年之后变稀，说明进入20世纪80年代以来，大沽河输沙量变幅趋缓，处于低输沙量期，与前述突变检验的结论一致。由图5（b）可知，输沙量变化的第一主周期为12 a，第二、第三主周期分别为21 a和7 a。

3气候变化与人类活动对大沽河水沙演化的贡献率

3.1对径流量变化影响的贡献率

由上述分析可知，大沽河径流量在1976年有1个突变点。因此，近似假定1956—1976年间是受人类活动影响相对较小、下垫面变化较小的基准期，建立这一时期的累积降水量—天然径流量双累积曲线关系，将1977—2016年累积降水量数据代入双累积曲线关系式，计算该时段在假设下垫面条件与1956—1976年相近条件下的累积径流量，从而可得到各年的拟合径流量R拟合。将1980—1989年、1990—1999年、2000—2009年、2010—2016年各时段的拟合径流量、天然径流量、实测径流量分别取均值，记为R拟合、R天然、R实测。R拟合与R天然、R实测对比关系见表4。

人类活动对地表径流的影响，可分为水资源开发利用活动的直接影响（主要是从河道内取水、修建水库闸坝等拦蓄工程）和下垫面渐变累积的间接影响。人类活动对流域径流量的直接影响可表示为实测径流量R实测与天然径流量R天然的差值，即对应表4中的“耗水量（直接影响）”一列；人类活动对径流量的间接影响可表示为1976年后R天然与R拟合差值，即下垫面的影响，对应表4中的“下垫面影响（间接影响）”一列；人类活动对径流量的综合影响为1976年后实测径流量R实测与R拟合的差值，即表4中“综合影响”一列。

由表4可见，与气候变化相比，人类活动是造成大沽河径流量变化的主要因素，影响贡献率在76.8%～97.4%之间。人类活动的直接影响和间接影响都促使径流量减小，人类活动的直接影响在20世纪90年代之前逐渐增大，2000—2016年趋于减小；人类活动的间接影响在1980—2016年期间呈现“增大—减小—增大”的趋势；人类活动的综合影响在1990—1999年较大，2000年以后比较稳定。

发生这种情况的主要原因是，20世纪80年代初期胶东半岛遭遇大旱，青岛市供水频频告急，大沽河水资源被大量开采。此后，大沽河流域与周边区域经济社会发展加速，河道内取水量不断增加，径流量不断减少，强烈的人类活动促使大沽河基本演变为季节性河流[7]。

3.2对输沙量变化影响的贡献率

由上述分析可知，输沙量有2个突变点，分别为1965年、1979年。据此将研究阶段划分为3个时期：A（1956—1965年）、B（1966—1979年）、C（1980—2016年）。由于大沽河流域在20世纪60年代以前受人类活动的影响很小，因此可将1956—1965年这一时期作为基准期，1966—1979年和1980—2016年两个时期为计算期。3个不同时期累积输沙量与年份、累积降水量与年份关系曲线见图6、图7。表5给出了累积输沙量和累积降水量的斜率及其变化情况。

气候变化对输沙量变化的贡献率为累积降水量-年份关系斜率变化率与累积输沙量-年份关系斜率变化率的比值[23-24]。因此，与基准期（A）相比，计算期（B）的降水量和人类活动对输沙量减小的贡献率分别为23.8%和76.2%，计算期（C）的降水量和人类活动对输沙量减小的贡献率分别为21.3%和78.7%；与计算期（B）相比，计算期（C）的降水量和人类活动对输沙量减小的贡献率分别为1.6%和98.4%。

除水利工程外，水土保持治理措施也是影响大沽河年输沙量的重要因素。大沽河上游主要为青岛平度市与莱西市。新中国成立后，该地区先后开展了农村“四旁”植树、荒山荒滩造林、山区小流域治理、农田林网建设等，尤其在1962年、1981年两次林木确权发证后，农民植树的积极性明显高涨。1985年，平度市森林覆盖率为10.7%，20世纪末达到20.78%，2005年为26.55%。1983年，莱西市森林覆盖率为7.97%，1985年为21.5%，2012年达到34.4%。增加植被覆盖率，增强下渗能力和蓄水能力，可有效减少流域侵蚀产沙，保持水土，从而减少河道内输沙量[25]。

4结论

1）1956—2016年大沽河实测径流量、天然径流量、输沙量均呈显著减少趋势，平均减少速率分别为0.118 3亿m3/a、0.123 1亿m3/a、1.927 3万t/a；降水量呈不显著減少趋势，平均减少速率为2.146 3 mm/a。

2）1956—2016年大沽河实测径流量、天然径流量有1个突变点（1976年），发生了由多到少的突变，突变年份前、后10 a均值分别减少了3.22亿m3、2.56亿m3，变幅分别为73%、40%。输沙量有两个突变点（1965年、1979年），均发生了由多到少的突变，突变年份前、后10 a均值分别减少了113.00万t、22.94万t，变幅分别为84%、91%。

3）大沽河径流量、输沙量均呈现明显的周期性震荡，天然径流量第一主周期为21 a，第二、第三主周期分别为12 a和8 a；输沙量第一主周期为12 a，第二、第三主周期分别为21 a和7 a。

4）人类活动是造成大沽河径流量、输沙量减少的主要原因，人类活动引起的下垫面变化对天然径流量变化的贡献率为76.8%～97.4%；人类活动对输沙量变化影响的贡献率为76.2%～98.4%。

[参考文献]

［1］ 胡春宏，王延贵，张燕菁，等.中国江河水沙变化趋势与主要影响因素[J].水科学进展，2010，21（4）：524-532.

［2］ 王赫，朱秀端，林敬兰，等.汀江流域长汀段水沙变化特征及其影响因素[J].水土保持学报，2022，36（2）：122-129，137.

［3］ ELHAM S，MAJID M，ARTEMI C.Sediment transport mechanisms and selective removal of soil particles under unsteadystate conditions in a sheet erosion system[J].International Journal of Sediment Research，2022，37（2）：151-161.

［4］ 李昌文，黄艳，严凌志.变化环境下长江流域超标准洪水灾害特点研究[J].人民长江，2022，53（3）：29-43.

［5］ 王亚迪，权全，薛涛涛，等.气候变化对黄河源区的水文影响分析[J].水文水资源，2018，7（2）：135-143.

［6］ 韩沙沙，谈广鸣，赵连军，等.黄河入海水沙情势变化分析[J].中国农村水利水电，2022（2）：1-5，19.

［7］ 姜德娟，王晓利.胶东半岛大沽河流域径流变化特征[J].干旱区研究，2013，30（6）：965-972.

［8］ 盛茂刚，黄修东，左林远，等.大沽河流域近60年降水量时空变化特征分析[J].水资源与水工程学报，2016，27（2）：65-68.

［9］ 金秉福，宫立新，宋键.大沽河泥沙来源的重矿物分析及其环境意义[J].海洋科学，2010，34（10）：71-76.

［10］ 薛伟，徐映雪，羊艳.青岛市大沽河地下水水源地应急能力分析[C].//水利部水文局，国际水文计划中国国家委员会，国际水文科学协会中国国家委员会，等.中国水文科技新发展：2012中国水文学术讨论会论文集.南京：河海大学出版社，2012：833-836.

［11］ 赵嘉阳，郭福涛，梁慧玲，等.福建长汀红壤区1965—2013年气温和降水量的变化趋势[J].福建农林大学学报（自然科学版），2016，45（1）：77-83.

［12］ 张润润.香港地区降水趋势及其演变过程分析[J].河海大学学报（自然科学版），2010，38（5）：505-510.

［13］ AHAMAD I，TANG D S，WANG T F，et al.Precipitation trends over time using Mann-Kendall and Spearmans rho Tests in Swat River Basin，Pakistan[J].Advances in Meteorology，2015（2015）：1-15.

［14］ 白勇.变化环境下半干旱流域径流演化特征及驱动因素分析[D].呼和浩特：内蒙古农业大学，2018：8-9.

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

［16］ 许继军，杨大文，雷志栋，等.长江流域降水量和径流量长期变化趋势检验[J].人民长江，2006，37（9）：63-67.

［17］ 马耘秀，夏露，范毅，等.土地利用和气候变化对黄土区典型流域水沙变化的影响[J].水土保持学报，2021，35（2）：38-46.

［18］ 魏凤英.现代气候统计诊断与预测技术[Ｍ].2版.北京：气象出版社，2007：63-65.

［19］ 夏伟，周维博，李文溢，等.气候变化和人类活动对沣河流域径流量影響的定量评估[J].水资源与水工程学报，2018，29（6）：47-52.

［20］ 魏晓婷，黄生志，黄强，等.定量分解气候变化与人类活动对季节径流变异的贡献率[J].水土保持学报，2019，33（6）：182-189.

［21］ 祁文燕，钱鞠，葛雷，等.湟水干流近60年径流变化特征分析[J].水资源与水工程学报，2018，29（3）：45-49.

［22］ 师忱，袁士保，史常青，等.滦河流域气候变化与人类活动对径流的影响[J].水土保持学报，2018，32（2）：264-269.

［23］ 郭文献，赵瑞超，付体江，等.乌江流域水沙变化趋势及驱动力分析[J].人民长江，2021，52（9）：71-78.

［24］ 王随继，李玲，颜明.气候和人类活动对黄河中游区间产流量变化的贡献率[J].地理研究，2013，32（3）：395-402.

［25］ 彭涛，田慧，秦振雄，等.气候变化和人类活动对长江径流泥沙的影响研究[J].泥沙研究，2018，43（6）：54-60.

收稿日期：2022-06-22

基金项目：山东省地矿局海岸带地质环境保护重点实验室2021年度开放基金项目（SYS202110）；西北农林科技大学旱区农业水土工程教育部重点实验室访问学者基金项目

第一作者：胡萌（1981—），男，山东日照人，高级工程师，博士，主要从事水资源方面的研究。

通信作者：雷炳霄（1984—），男，陕西渭南人，高级工程师，学士，副院长，主要从事水工环地质方面的研究。

E-mail：bingxiao8426@163.com

（责任编辑杨傲秋）