黄河源区雨季降水与汛期径流量重建及其千年尺度下的演变特征

2022-02-13王文卓张建云贾本有吴淑君王国庆

水科学进展 2022年6期

王文卓，张建云，陈峰，贾本有，吴淑君，王国庆

(1.河海大学水文水资源学院,江苏南京 210098；2.河海大学长江保护与绿色发展研究院，江苏南京 210098；3.水利部应对气候变化研究中心,江苏南京 210029；4.云南大学国际河流与生态安全研究院,云南昆明 650500；5.南京水利科学研究院水文水资源与水利工程科学国家重点实验室,江苏南京 210029)

径流及其影响因子的样本长度及代表性，对准确揭示河川径流过去长期演变特征、发现大尺度下径流对气象因子的响应特征以及精确预测未来变化趋势有重要影响[1]。受限于水文气象设站较晚的现实情况，中国绝大多数测站累积的观测数据长度不足百年，缺乏古径流及古降水系列(百年及千年)限制了对降水及河川径流长期演变及其响应规律的深入理解。树木年轮(以下简称树轮)由于具有时间序列长、站点分布广、连续性强、分辨率高和野外测量简单等优势，成为量化过去千年气候水文变化最常用的代用指标之一，基于树轮的水文气象要素重建研究一直以来都是国内外研究的热点[2-5]。

目前，利用树轮资料重建气候水文变量多采用线性回归法，包括简单线性回归法、多元线性回归法、典型回归法和逐步线性回归法等，方法较为单一，且无法考虑不确定性。鉴于此，Devineni等[6]提出了一种分层贝叶斯回归模型，在考虑参数不确定性的情况下，给出了各重建序列的后验概率分布估计，降低了径流重建的不确定性；沈大军等[7]同样采用基于部分汇集的层次贝叶斯回归方法，仅利用少量观测数据重建了印度河上游流域3个水文站点的径流量，解决了该地区可用观测序列过短难题，从而证明层次贝叶斯回归考虑了重建的不确定性，可以在观测序列较短或有缺失值时重建径流。

黄河源区(本文指唐乃亥水文站以上区域)是黄河流域内最大的冰雪融水和降雨补水区[8]，河川径流量占全流域的35%，准确掌握黄河源区径流对气候因子的响应规律及其演变趋势有助于合理确定全流域的水资源利用及分配战略[9]。Gou等[10]和勾晓华等[11]基于祁连圆柏树，通过主成分线性回归分别重建了黄河源区过去593 a和1 234 a的年径流序列，重建结果显示15世纪晚期、1820—1830年和1480—1490年为低流量时段;Wang等[12]将嵌套主成分回归与最优子集选择法相结合，重建了黄河源区过去近2 000 a的年径流序列，并揭示了温暖气候条件下易出现高流量时期，而寒冷气候条件下则易出现低流量时期;Chen等[13]重建了黄河上游过去254 a(上年8月至当年6月)的降水序列，结果显示1873—1887年、1898—1923年和1989—2003年间降水较弱，而1769—1785年、1798—1833年、1924—1938年、1957—1968年和2004—2013年间降水较强。从当前见诸文献的成果来看，黄河源区径流重建研究主要集中在采用线性回归法的年径流量重建，尚未见到针对年内特定季节的水文气象变量重建研究。然而黄河源区雨季降水(5—9月)占全年的85%，汛期径流量(7—10月)占全年的59%，且宁蒙河段新的“地上悬河”正在形成，开展黄河源区雨季降水与汛期径流量重建及其非线性响应规律分析研究具有十分重要的理论意义和实用价值。相比于年尺度的水文气象变量，树轮生长对雨季降水和汛期径流量的响应更为复杂，利用树轮数据对其进行重建存在更大的不确定性。

本文建立考虑不确定性条件下适用于季节性水文气象变量的重建模型，展延黄河源区雨季降水，并在已有唐乃亥年径流量重建成果的基础上[12]辨识其年内分布规律，提出基于年径流量的汛期径流量重建模型以规避利用树轮资料重建季节性水文变量的不确定性，获得黄河源区汛期径流量重建序列，扩充黄河源区水文气象数据库，并揭示千年尺度下黄河源区雨季降水和汛期径流量的演变及响应规律，为流域径流预报预测提供有力支撑，推动黄河流域生态保护和高质量发展。

1 研究区概况与数据来源

1.1 研究区概况

黄河源区集水面积为12.20万km2，占黄河流域面积(不含闭流区)的16.2%，河长约300 km，为黄河主要产水区，被形象地称为“黄河的水塔”。黄河源区坐落于海拔高、气温低的青藏高原东北边缘(图1)，对气候变化敏感，是树轮气候学研究的理想场所。

图1 黄河源区及树轮站点地理信息Fig.1 Map of geographic information of the headwater catchment of Yellow River basin and tree ring stations

黄河源区气温和降水受孟加拉湾亚洲夏季风的影响，季节分布很不均匀，每年5—9月气温和降水都迅速上升和增加，气温为-4～2 ℃，降水为250～750 mm，同时黄河源区降水的空间分布存在显著差异，西北部年平均降水量仅为200 mm，但东南部可能超过700 mm[14]。

1.2 数据来源及其预处理

主要数据包括天然径流数据、树轮数据、气象网格数据以及大尺度气候数据，各数据站点地理位置如图1。黄河源区附近7个气象站1960—2016年的降水数据从中国地面气候资料日值数据集(V3.0)中获取；贵德水文站(101°24′E,36°02′N)1920—1976年以及唐乃亥水文站(100°09′E,35°30′N)1956—2016年的天然月径流数据从黄河水利委员会获取；唐乃亥159—1919年的径流序列为已有研究成果[12]；黄河源区附近16个树轮宽度年表从世界树轮银行及邵雪梅等[15]研究中获取，基本信息如表1所示。

表1 数据信息表Table 1Data information chart

采用泰森多边形法计算黄河源区降水面值。为得到更加稳定的模型，采用上下游相关法，以贵德水文站为参证站，根据其1920年以来的数据延长唐乃亥站径流序列。贵德水文站1954年建站以前的径流数据为黄河水利委员会利用兰州站1934—1953年的观测数据及兰州站1920—1933年的插补数据插补得到，该数据已用于龙羊峡水电站的设计与建设及国务院1987年批复实施的《黄河可供水量分配方案》的制定，具有较强的合理性及可靠性，其插补可靠性分析详见文献[12]。为与已有研究成果保持一致，选取上年11月至当年10月为1个水文年，延长后最终获得1921—2016年共96 a唐乃亥站水文年径流量。采用“Signal-free”法对树轮宽度原始数据进行去趋势，在尽可能保留低频气候信号的前提下去除树木的内在生长趋势，获得树轮标准年表[16]。该方法通过反复迭代得到平稳的年轮宽度指数序列，再对轮宽序列求平均进而建立零信号年表，可减少趋势扭曲问题，并保留较多的轮宽低频变化信息[17]。

1.3 研究区径流量年内分布特征

黄河源区是典型的高寒区，整体气温偏低，大面积的冰川和冻土不仅是径流量的重要补给来源，还影响径流的下渗作用，且受夏季西南及东南气流影响，年内降水集中。在这样复杂的水循环规律下，黄河源区径流量年内分布总体呈现夏秋多水、冬春少水的季节性分布特征。依据1920—2016年径流序列，绘制年内变化箱线图如图2(a)所示，4—10月的径流量占年径流量的83%，其中7—9月占比最大。7月径流量平均占年径流量的17.2%，平均径流量为35.03 亿m3，其中，最高占年径流量的28.7%，当月径流量为56.86 亿m3，当年径流量为198.4亿m3;8月径流量平均占年径流量的14.2%，平均径流量为28.65亿m3，其中，最高占年径流量的22.1%，当月径流量为45.61亿m3，当年径流量为206.8亿m3;9月径流量平均占年径流量的14.7%，平均径流量为30.20亿m3，其中，最高占年径流量的30.9%，当月径流量为92.11 亿m3，当年径流量为297.8亿m3。

唐乃亥汛期径流量(7—10月)与年径流总量关系如图2(b)。汛期与非汛期径流量在枯水年和丰水年的对比如图2(c)所示，枯水年汛期与非汛期多年平均径流量分别为93.56亿m3与70.06亿m3，而丰水年分别为151.4亿m3与90.08亿m3;与丰水年相比，枯水年汛期与非汛期径流量均较小，与年径流量变化一致。

唐乃亥汛期径流量年内占比受年径流量大小影响。汛期和非汛期径流量及其占比在枯水年和丰水年的对比如图2(d)所示，枯水年汛期径流量占比最小降至0.40，丰水年汛期径流量占比最高达0.74;枯水年汛期与非汛期径流量占比差距小，丰水年汛期与非汛期径流量占比差距大。

图2 黄河源区径流年内分布特征Fig.2 Intra-annual distribution features of streamflow of the headwater catchment of Yellow River basin

1.4 树轮年表的选择及相关性分析

树轮年表主要依据地理位置、年表时间跨度及与研究区降水的相关性等来进行选择。由于降水的空间一致性不受限于流域范围，没有严格的界限[18]，因此挑选黄河源区内及附近的树轮作为备选。鉴于大部分树轮年表的截止年份晚于2001年，选取截止年份晚于2001年的树轮年表以将尽可能多的年表用于重建。考虑树轮生长对降水响应的滞后性，同时分析同期树轮和滞后1 a树轮年表与降水的相关性，选取相关性显著性水平达0.05的年表[19]，但为了避免过拟合，仅选取其中相关性较强的1个。最终共筛选出16个树轮年表用于降水重建，基本信息如表2所示。筛选结果显示少部分源区以外树轮同样可以反映黄河源区降水特征，且其地理位置相对集中，表明该区域与黄河源区降水特征相对一致。相关性分析表明树轮与黄河源区雨季降水呈正相关，表明所选树木生长情况受当地气候干湿条件影响显著，且雨季降水对树轮生长存在1 a的滞后影响期，即当年雨季降水对次年树木生长影响显著。由于树轮采样区地处青藏高原和黄土高原的“干旱-半干旱”过渡地带，植被对于地表水和与之紧密联系的地下水变化异常敏感，雨季降水可以有效补充地表及地下水量，是树木生长的主要影响因子之一。

表2 树轮年表基本信息Table 2Information of tree-ring chronologies

续表

2 研究方法

2.1 嵌套主成分层次贝叶斯回归模型

为了尽可能挖掘不同时间跨度的树轮信息并延长重建时间，本文根据树轮年表开始的年份划分出一系列不断向过去延伸的嵌套时段。在划分嵌套时段时，依据可用树轮年表开始年份的分布情况，选取2个相邻嵌套时间段开始年份间隔的最小值，以避免大量冗杂的重建过程。由于可用树轮年表的增加，开始年较晚的时段重建通常比开始年较早的时段重建更可靠。因此，通过不断将开始年更为早远的时段重建值补充至开始年较晚的时段重建序列获得最终的重建值，具体示意如图3所示。

图3 嵌套重建方法示意Fig.3 Schematic diagram of nested reconstruction method

主成分层次贝叶斯回归模型首先对具有自相关性的树轮年表集做主成分分析，将年表集转换成一系列独立的主成分，并采用双侧t检验方法，挑选满足置信水平为α=0.05的主成分作为模型的回归因子，以提高模型效率、降低计算维度及过滤无关信息[20]。

第二步建立层次贝叶斯回归模型。层次贝叶斯回归模型具有双层次结构，其回归参数为一系列超参数描述的先验分布函数[6]。为使先验信息较弱，选取均值为0、方差为104的正态分布作为回归系数的先验分布[21]。嵌套层i内，t年雨季降水Pi(t)的主成分层次贝叶斯回归模型为：

(1)

以及先验分布函数：

αi～N(0,104)
βi～MVN(0,104I)
σi～C(0,2)

(2)

式中：Uij(j=1,2,…,n)为嵌套层i中与雨季降水Pi(t)显著相关的主成分；ni为嵌套层i中相关性显著的主成分个数；εi(t)为嵌套层i中回归模型的残差,符合均值为0、σi为方差的正态分布N(0,σi);αi为嵌套层i中回归模型的截距；βi中为嵌套层i中主成分的回归系数,其先验分布为均值为0、方差矩阵为104I的多元正态分布MVN(0,104I),I为单位矩阵;C(0,2)为位置参数为0、尺度参数为2的柯西分布。

2.2 分类占比回归法

分类占比回归法以不同丰枯水平年汛期径流量占比特征为基础，按丰枯水平分类，分别建立丰水年和枯水年汛期径流量与年径流量的线性回归方程。由于年径流量是依据树轮年表重建获得，建立汛期径流量与年径流量的一元回归函数，以避免重复利用树轮信息。考虑到当年径流量为0时，汛期径流量也为0，将截距项设为0。建立汛期径流量与年径流量的一元回归方程如下：

(3)

2.3 检验统计量

选取误差缩减值(ER)和有效系数(EC)检验重建结果，其值大于0代表重建结果精度高于均值，且越大代表重建结果与记录数据越相似[22]。设验证期年数为m，ER和EC的公式如下[23-24]：

(4)

(5)

3 结果及分析

3.1 雨季降水重建

3.1.1 嵌套时段划分及其参数估算

根据所筛选出的16个树轮年表起始年份，使嵌套时段开始年份至少相差30 a，划分出10个嵌套时段，时段1—时段10起始年份依次为1300年、1230年、1163年、1130年、943年、900年、857年、681年、470年及159年，分别对每个嵌套时段建立主成分层次贝叶斯回归模型。各时段的起始年份、可用树轮个数、相关性强的主成分个数及第一主成分相关系数如图4所示，起始年份最早的时段为159—2000年，其可用树轮数量最少，只有1个。随着起始年份不断后延，各嵌套时段内可用的树轮个数不断增加，但最终用于回归的主成分个数与树轮个数没有明显的相关关系；各嵌套时段第一主成分与汛期径流量的相关性均满足95%的置信区间，且除第10嵌套时段外其余嵌套时段第一主成分与汛期径流量的相关系数绝对值均在0.45以上，相关性较强。

图4 各嵌套时段基本信息Fig.4 Information of different nested periods

将降水记录数据时段(1960—2000年)分为校准期(1971—2000年)和验证期(1960—1970年)，利用校准期数据校准模型参数。各嵌套时段层次贝叶斯回归模型的参数后验分布计算结果如图5所示。其中时段9和时段10中回归因子为单因子，回归系数平均值分别为18.02和20.52，表明时段内树轮年表集主成分与雨季降水呈正相关，这与图4相关性分析结果一致。

图5 各嵌套时段主成分层次贝叶斯回归模型参数后验分布估算Fig.5 Posterior distribution estimations of parameters in principal component hierarchical Bayesian models in each nested periods

3.1.2 模型检验及雨季降水重建

利用建立的嵌套主成分层次贝叶斯回归模型重建验证期雨季降水。重建的验证期雨季降水与记录值对比如图6所示，重建序列与记录序列波动特征较为一致，除1967年雨季降水记录值演替特征为增长而重建值为减少外，其余年份重建序列均能够较为真实地反映记录序列的演变规律；重建序列验证期的均值为75.4 mm、中位数为75.9 mm、最小值64.7 mm，记录序列的均值为77.7 mm、中位数为79.7 mm、最小值66.1 mm，均值、中位数及最小值统计值相似，说明模型能够较好地重建降水总体情势，且对干旱年份降水模拟效果较好。但与记录序列相比，重建序列标准差略小，最大值偏小，1967年和1968年的强降水年重建精度不足，这表明树轮由于生长所需水分有限，且当年降水过程及气温、气压等综合气候条件不利，无法精确反映特大降水年份的降水总量。

图6 验证期雨季降水重建值与记录值对比Fig.6 Comparison of reconstruction and observation of rainy season precipitation in verification period

同时利用EC和ER进一步检验模型效果。检验结果如表3所示，各嵌套层EC和ER均显著高于0。除嵌套时段8—时段10外，各嵌套时段EC和ER值均高于0.43，EC值最高达0.58，ER值最高达0.59。为保证重建结果的可靠性，删去嵌套时段8—时段10。

表3 各嵌套时段EC和ER检验值Table 3Values of EC and ER evaluation indices of different nested periods

基于筛选出的16个树轮年表，利用建立的嵌套主成分层次贝叶斯回归模型，展延黄河源区雨季降水至857年，并采用5 a滑动平均和11 a滑动平均方法对重建序列进行一定平滑处理以突出低频变化。重建结果及其95%和5%频率线如图7所示。

图7 黄河源区857—2016年雨季降水重建序列Fig.7 Reconstruction of rainy season precipitation of the headwater catchment of Yellow River basin in period of 857—2016 C.E.

分析5 a滑动平均值可知，近期没有出现千年尺度上的雨季降水极端时期，而在896—901年、907—914年、1141—1145年和1491—1495年分别发生了雨季弱降水时期，870—874年、969—974年、982—988年、1343—1349年和1681—1685年分别发生了雨季强降水时期。分析11 a滑动平均值可知，9世纪晚期到10世纪初期和15世纪中期到后期整体雨季降水较弱，10世纪中期到晚期、14世纪中期到晚期及16世纪早期到晚期整体雨季降水较强，而21世纪前后也出现了较为干旱的气候。其中，10世纪中期到晚期发生了最为强烈的持续强降水时期，期间共13 a的雨季降水量超过5%频率线;14世纪中期到晚期和16世纪早期到晚期出现了持续时间最长的湿润时期，期间分别有44 a和51 a雨季降水量远超多年平均值，3个持续强降水时期成果均与Yang等[25]提出的过去3 500 a祁连山降水重建序列一致;而15世纪中期到后期出现了持续时间最久的干旱时期，雨季降水持续偏弱，期间1464—1498年共35 a的多年平均雨季降水量仅66.9 mm，这一结果与Gou等[26]重建的祁连山6—7月降水序列一致。

3.2 汛期径流量重建

3.2.1 模型构建

以1943—2001年共58 a的天然径流数据作为校准期数据，建立汛期径流量与年径流量的一元回归方程如下：

(6)

汛期径流量与年径流量的一元回归方程回归系数在丰水年和枯水年分别为0.625和0.553，表明丰水年汛期径流量占全年径流量的比例约为62.5%，而枯水年汛期径流量约占55.3%。

3.2.2 模型检验及汛期径流量重建

利用构建的基于年径流的一元回归模型重建验证期汛期径流量。重建的验证期汛期径流量与记录值对比如图8所示。对比结果表明，基于年径流量的一元回归模型精度高于基于树轮的主成分层次贝叶斯回归模型，其重建序列呈现与记录序列几乎一致的波动特征，重建精度在各个年份相对稳定，较好地规避了利用树轮重建的不确定性；重建序列验证期的径流量均值为106.2亿m3、中位数为95.50亿m3、标准差31.38亿m3，记录序列的径流量均值为111.0亿m3、中位数为105.5亿m3、标准差33.22亿m3，各统计值相似。可见，利用汛期径流量的年内占比重建汛期径流量可以较好地避免利用树轮信息重建古径流系列的不确定性，并与年径流量重建结果保持一致。计算ER和EC值可得，重建结果ER值达0.90、EC值达0.88，重建精度高。

图8 验证期汛期径流量重建值与记录值对比Fig.8 Comparison of reconstruction and observation of flood season streamflow in verification period

利用构建的一元占比回归模型重建汛期径流量，重建黄河源区汛期径流量序列至公元159年，并采用5 a 滑动平均和11 a滑动平均方法对重建序列进行一定平滑处理以突出低频变化。重建结果及其95%和5%频率线如图9所示。

图9 黄河源区公元159—2016年汛期径流量重建序列Fig.9 Reconstruction of flood season streamflow of the headwater catchment of Yellow River basin in period of 159—2016 C.E.

分析5 a滑动平均值可知，即使在千年时间尺度下，20世纪七八十年代(1979—1985年)亦是较为不寻常的汛期高径流量时期。除此之外，1107—1113年、1141—1146年、1477—1482年和1713—1717年为汛期低径流量时期，196—201年、230—236年、569—574年、1229—1234年、1243—1249年、1900—1906年和1979—1985年为汛期高径流量时期。分析11 a滑动平均值可知，5世纪中晚期、7世纪早期到8世纪中期、15世纪早期到晚期及1700年前后整体汛期径流偏低，3世纪晚期到4世纪早中期、6世纪晚期、12世纪晚期到13世纪早期、16世纪中期整体汛期径流偏高。其中15世纪早期到晚期发生了较为强烈的持续低径流量时期，这与大量现有研究一致，该时期中国北方大部分地区天然来水偏枯[10，27-28]。

采用50 a低通滤波器对雨季降水和汛期径流量重建序列进行处理，其Z值序列如图10所示。雨季降水与汛期径流量呈现相似的演变特征，雨季强降水时期与汛期高径流时期、雨季弱降水时期与汛期低径流时期几乎一致。雨季降水较为稀少的9世纪晚期到10世纪初期和较为充沛的10世纪中晚期即使没有发生千年尺度上显著的汛期低径流或高径流，但在百年尺度上亦发生了显著的极端汛期径流。同时Z值图显示黄河源区雨季降水和汛期径流量在15世纪以后出现先减小后增长的趋势。结合全球气候变化可知，雨季降水与径流量减少的15世纪到19世纪属于气候寒冷的小冰期时期，20世纪暖期开始后雨季降水与径流量逐渐开始回升[29]。该结果表明黄河源区雨季降水和径流量受大尺度气候影响，近现代温暖气候下黄河源区雨季降水和径流量呈增长趋势。未来全球气温不断升高的趋势下，黄河源区气候较大概率向总体暖湿化发展[30]，极端降水和极端干旱事件发生频率增高，流域来水非一致演变日渐强烈，流域防洪工程的运行安全和运行参数需进一步校核，水利工程调度运行规则需进行必要的调整，进行流域径流预报、水量配置等水资源管理工作时迫切需要进一步考虑来水非一致性演变的影响。

图10 黄河源区雨季降水和汛期径流量重建序列Z值对比Fig.10 Z-scores of reconstructions of rainy season precipitation and flood season streamflow of the headwater catchment of Yellow River basin