引洮工程受水区径流演变规律分析

Analysis of Runoff Variation Law in the Water Receiving Area of Yintao Water Transfer Project

摘要: 根据引洮工程受水区四个水文站的1956~2022年实测径流数据，运用趋势检验、小波分析、贝叶斯变点分析模型等方法分析径流年际变化，采用不均匀系数、集中度和集中期等表征径流年内分配情况并研究其变化趋势，结果表明：受水区年径流整体呈显著减小趋势，在年代际差异上存在较强跳跃性，最大概率突变年在1990年前后，此后减小趋势愈发显著；年径流大尺度主周期振幅在2005年以后显著减小，小尺度主周期振荡愈发剧烈，呈现更复杂的变化态势，预计未来2~3年内持续减小；汛期径流显著减小，径流年内分配不均匀程度和集中程度降低。整体来看，受水区径流显著减小，年内分配发生变化，周期振荡愈发复杂，水资源安全保障不确定性加大，亟需优化调控引洮工程供水与水资源配置。

Abstract: Based on observed runoff data from four hydrological stations in the receiving water area of the Yintao project from 1956 to 2022, a comprehensive analysis of annual runoff inter-annual variability was con-ducted using methods such as trend analysis, wavelet analysis, and Bayesian change-point models. Addi-tionally, characteristics of intra-annual runoff distribution, including the uneven coefficient, concentration degree, and concentration period, were employed to analyze their temporal changes. The results are as follows: The annual runoff in the receiving water area exhibits a significant decreasing trend, with pro-nounced decadal variability in runoff changes. The highest probability of a change-point occurring was around 1990, after which the declining trend became more prominent. On a larger scale, the amplitude of the main periodic oscillations in annual runoff showed a significant reduction after 2005, while small-er-scale oscillations became more intense. Consequently, the annual runoff displayed a more intricate changing pattern, with a projected continuation of reduction in the next 2~3 years. There is a significant decrease in flood season runoff in the receiving area, accompanied by a decrease in the unevenness and concentration degree of intra-annual distribution. Overall, the runoff in the water receiving area signifi-cantly reduced, the intra-annual distribution changes, the cycle oscillation is more complex, the uncer-tainty of water resources security increases, there’s an urgent need to optimize the regulation of the water supply of Yintao project and water resources allocation.

1. 引言

全球变暖加速了水循环及其时空变化 [1] [2] ，径流作为水循环的关键组分，具有高度不确定性 [3] 。同时，河川径流也是地表水资源的主要形式，由于干旱半干旱地区径流较湿润地区对气候变化和人类活动更敏感 [4] [5] [6] [7] ，全球变暖进一步加剧了水资源供需矛盾，水资源开发利用面临新的风险和挑战。调水工程是解决区域水资源短缺问题的重要途径之一，研究干旱半干旱地区调水工程受水区的径流演变规律对调水工程的规划运行与区域水资源开发利用具有重要意义。

引洮供水工程受水区位于甘肃省中部黄土丘陵区，气候干旱，严重缺水造成当地人畜饮水困难、经济发展缓慢。引洮工程作为甘肃有史以来规模最大的跨流域调水工程，是解决中部地区水资源短缺问题、保障和改善民生、促进经济社会可持续发展的重大战略性基础设施。研究引洮工程受水区径流演变规律，不仅是引洮工程优化调度运行的需要，而且也是保障甘肃省中部黄土丘陵区社会经济发展中水资源安全的重大现实需求。本文以引洮工程受水区1956~2022年站点径流数据为基础，运用趋势检验、小波分析、贝叶斯变点分析模型等方法分析径流年际变化，采用不均匀系数、集中度、集中期等表征径流年内分布特征并分析其变化情况，以期为引洮工程的效益发挥及受水区内的水资源配置提供科学依据。

2. 研究区概况及基础数据

2.1. 研究区概况

引洮供水工程受水区属黄土高原丘陵沟壑区，气候干旱，植被覆盖率低，生态环境脆弱，水土流失严重，区内产流少，水资源极度紧缺。受水区分属洮河、宛川河、祖厉河及渭河流域(北道以上)。引洮供水工程从洮河九甸峡水利枢纽有坝自流引水，受益区涉及兰州市榆中县，白银市会宁县，定西市临洮县、渭源县、安定区、陇西县、通渭县，天水市武山县、甘谷县、秦安县，以及平凉市静宁县共5个市辖的11个设计供水县(区)。引洮供水工程受水区地形、河流水系、工程界限等基本情况如图1(a)所示。

Figure 1. Schematic diagram of study region: (a) Topographic and administrative division map; (b) Basin river system and hydrological station map

图1. 研究区域示意图：(a) 地形与行政区划图；(b) 流域水系与站点分布图

2.2. 数据来源

选取引洮工程受水区内李家村、会宁、郭城驿、武山4个干流水文站1956~2022年逐月实测径流数据，根据各水文站集水面积统计得到径流深序列。水文站分布情况如图1(b)所示。

研究区流域水系与水文站集水区由数字高程模型(Digital Elevation Model, DEM)数据分析提取得到。DEM采用美国国家航空航天局2020年发布的全球30 m分辨率数据(NASA DEM)。

3. 研究方法

3.1. 预置白处理的Mann-Kendall检验法

Mann-Kendall (MK)统计检验 [8] 由Mann和Kendall提出的一种非参数趋势检验方法，广泛应用于诊断水文气象时间序列是否存在单调趋势和突变点，克服了水文数据偏态、非同分布、有异常值等问题 [9] 。TFPW-MK方法在进行MK检验前，采用趋势预置白法 [10] 消除序列正向自相关性，避免趋势显著性被放大，已广泛应用于径流方面 [11] [12] [13] 。采用TFPW-MK方法对年径流以及径流年分配特征指标进行趋势检验和变点识别。

3.2. 贝叶斯变点分析模型

MK检验法在识别时间序列的趋势变动和突变定性判断方面具有较强优势，而贝叶斯变点分析模型 [14] 则能进一步量化突变发生的后验概率 [15] 。贝叶斯变点分析模型假设长度为n水文序列的观测值x_i服从某一正态分布，在变点时刻k前后，x_i服从如下分布：

$x_i\sim N\left({\mu }_a,{\sigma }_a^2\right)\left(i=1,2,\cdots ,k\right)$ (1)

$x_i\sim N\left({\mu }_b,{\sigma }_b^2\right)\left(i=k+1,2,\cdots ,n\right)$ (2)

假设方差不变，只研究均值变化，即σ_a = σ_b = σ。假设σ_a、σ_b服从均值为μ、方差很大的正态分布，趋近于均匀分布(σ和μ由实测序列统计)。则可由实测值X_i推导出μ_a的后验分布：

$\begin{array}{l}{\mu }_a|X_i\sim N\left({\mu }_a^{\ast },{\sigma }_a^{\ast 2}\right)\left(i=1,\cdots ,k\right)\\ \begin{array}{c}{\mu }_a^{\ast }=\left(n^{\ast }\mu +{\displaystyle \underset{i=1}{\overset{k}{\sum }}{x}_i}\right)/\left(n^{*}+k\right)\\ {\sigma }_a^{\ast 2}={\sigma }^2/\left(n^{*}+k\right)\\ n^{*}={\sigma }^2/{\sigma }_a^2\end{array}\end{array}$ (3)

同理可得μ_b的后验分布，则观测值发生的联合分布函数为：

$p\left(X|k,{\mu }_a,{\mu }_b\right)={\displaystyle \underset{i=1}{\overset{k}{\prod }}\frac{1}{\sqrt{2\text{π}}\sigma }\exp \left[-\frac{{\left(x_i-{\mu }_a\right)}^2}{2{\sigma }^2}\right]}\cdot {\displaystyle \underset{i=k}{\overset{n}{\prod }}\frac{1}{\sqrt{2\text{π}}\sigma }\exp \left[-\frac{{\left(x_i-{\mu }_b\right)}^2}{2{\sigma }^2}\right]}$ (4)

根据贝叶斯法则，突变发生位置k的后验分布密度函数为：

(5)

最终得到后验概率的期望值为：

$p\left(k|X\right)={\displaystyle \iint p\left(k|X,{\mu }_a,{\mu }_b\right)d}{\mu }_{a}d{\mu }_b$ (6)

3.3. 小波分析

小波函数φ(t)是指具有振荡特性、有限长或能迅速衰减到零的一类函数。小波分析是将某一时间序列信号分解为基本小波函数φ_a,b(t)经伸缩或平移形成的一簇函数。a为尺度(周期)因子，反映小波的周期长度；b为时间因子，反映在时间上的平移。

对于二维时间序列f(t)，其连续小波变换为：

$W_f\left(a,b\right)={\left|a\right|}^{-1/2}{\displaystyle {\int }_{-\infty }^{+\infty }f\left(t\right)\bar{\phi }\left(\frac{t-b}{a}\right)\text{d}t}$ (7)

式中， $\bar{\phi }\left(t\right)$ 为φ(t)的复共轭函数，W_f(a,b)为小波转换系数。

小波变换能同时反映时间序列f(t)映射到时域和频域特性。W_f(a,b)随参数a和b变化，可作出时频域一定范围内W_f(a,b)实部的二维等值线图。W_f(a,b)实部时频图反映了不同时间和尺度下小波系数的变化特征：小波系数绝对值越大，表明该时间尺度变化越显著 [16] 。

小波方差反映信号波动能量强弱，表示为：

$Var\left(a\right)={\displaystyle {\int }_{-\infty }^{+\infty }{\left|W_f\left(a,b\right)\right|}^2}\text{d}b$ (8)

3.4. 不均匀系数

采用径流年内分配不均匀系数(或称完全调节系数) C_v来量化径流年内分配的不均匀性，其计算公式为：

(9)

式中，R(t)为年内各月径流量， $\bar{R}$ 为年内平均月径流量。

3.5. 径流集中度与集中期

把一年内各月径流看作矢量，合矢量模与月径流均值之比为年径流集中度，合矢量方向对应年径流集中期 [17] 。集中度与集中期量化表征径流年内分配的集中程度和集中重心出现的时间。

4. 结果与讨论

4.1. 径流年际变化规律

洮河流域李家村站、渭河流域武山站、祖厉河流域郭城驿站和会宁站1956~2022年年径流深的统计特征如表1所示，图2(a)~(d)给出了年际变化过程。TFPW-MK趋势检验统计量Z与线性趋势显示四站年径流序列整体均呈波动减小趋势，且Z值均通过了1%显著性水平检验。变差系数介于0.34~0.83之间，年际波动较显著。

Table 1. Statistical characteristics of station-measured runoff

表1. 年径流深统计特征

Figure 2. Station-measured runoff trend indicators and their inter-annual variation process: (a)~(d) Trends of annual runoff depth; (e)~(h) Statistics of mutation test; (i)~(l) Cumulative anomalies and Bayesian posterior probabilities

图2. 年径流深相关趋势指标及其变化过程：(a)~(d) 年径流深变化趋势；(e)~(h) 突变检验统计量曲线；(i)~(l) 累积距平和贝叶斯后验概率

位于洮河流域干流的李家村站集水面积和多年平均径流深最大，极值比和变差系数最小，年际振荡较小，Z的绝对值也最小，径流深年均减小速率为1.32 mm/a；会宁站位于祖厉河流域河源区的支流河口处，集水面积和多年平均径流深最小，最小值仅0.27 mm (2015年)，极值比高达161.85，变差系数和Z的绝对值都明显大于其他水文站，年际振荡最大，减小趋势也最为显著。除会宁站外，其余三个干流水文站年径流深的最大观测值均出现在1967年，1997年达到一定低值后均有所回升，会宁站年径流深则持续减小，乃至2015年后出现全年断流现象。

年径流深突变检验统计量变化曲线如图2(e)~(h)所示。1970年以前，李家村和武山站年径流深呈现上升趋势但变化趋势并不显著。郭城驿和会宁站则没有明显的增大或减小趋势；1975年左右，武山和会宁站年径流深开始减小，李家村和郭城驿站径流深的减小趋势则在1985年左右出现；1985~1990年间李家村、武山和郭城驿站年径流序列发生明显突变，突变发生后年径流减小趋势愈发显著，并逐渐超过1%显著性水平，图2(c)中的累积距平曲线直观反映了这一变化。会宁站年径流深变化虽然也在年代际差异上存在较强跳跃性，但顺逆序统计量交叉点位于5%显著性水平界线外，未发生明显突变。

Figure 3. Morlet wavelet analysis of station-measured runoff: (a)~(d) Time-frequency distribution from the real part of wavelet transform coefficients; (e)~(h) Time-frequency distribution from the modular square of wavelet transform coefficients; (i)~(l) Curves of wavelet variance

图3. 年径流深Morlet小波分析结果：(a)~(d) 小波变换系数实部时频图；(e)~(h) 小波变换系数模平方时频图；(i)~(l) 小波方差曲线

采用贝叶斯变点分析模型进一步量化年径流发生突变的后验概率，结果如图2(i)~(l)所示。李家村和武山站最大概率突变年分别是1985和1986年，郭城驿站则是1992年。TFPW-MK突变检验结果表明会宁站年径流序列突变特征不明显，但1978突变发生的贝叶斯后验概率异常偏高，由于贝叶斯后验概率与变点前后序列统计特征(均值)的差异大小呈正相关，这可能是此前出现的几次峰值扰动所致。

图3(a)~(d)中变换系数实部的分布反映了不同特征时间尺度下的径流丰枯变化特性，图3(i)~(l)的小波方差曲线进一步刻画了不同时间尺度的波动能量强弱。由图可知，祖厉河流域两水文站实测年径流的波动能量较弱，周期特征较不明显，主周期尺度也较小，位于支流的会宁站年径流3年和12年小尺度主周期振荡最显著，位于干流的郭城驿站也存在3年尺度的主周期振荡，但波动能量不大，第一和第二主周期分别为14年和8年，两水文站在30年处均表现出较明显的正负位相交替现象，但等值线未闭合，还需更长序列的数据验证。洮河流域李家村和渭河流域武山站年径流存在相似的周期特征，第一主周期分别为20年和22年，大尺度的振荡更强烈，并在其主导下呈现出较明显的“枯–丰–枯–丰–枯–丰–枯”演变特征。

小波变换后的模平方时频分布如图3(e)~(h)所示，刻画了不同特征时间尺度周期信号的强弱分布。李家村站年径流15~20年尺度的周期信号最强，振荡中心在1970年左右，但1985年前后年径流发生突变后逐渐减弱，5~10年的周期信号逐渐变强，振荡中心在2017年左右；武山站以15~25年的波动信号为主，兼有5年以下的周期信号，但振荡范围也集中在1985年之前，此后大尺度周期信号显著减弱，2010年以后虽有5~10年周期信号产生但并不强烈；1985年以前，郭城驿和会宁站均以5年以下的周期信号为主，此后郭城驿站5~10年的周期信号渐强，是影响目前年径流波动的主要尺度，会宁站则不再表现出特定尺度的周期振荡。

为进一步刻画连续小波转换后主周期尺度下年径流的波动特性，在小波变换系数实部时频图上取主周期做平行于时间轴的切割线，采用双调和样条插值法沿切割线作Morlet小波系数实部随时间变化的过程线，叠加在年径流深距平线上，如图4所示。

Figure 4. Anomalies of station-measured runoff depth and curves of the real part of Morlet wavelet transform coefficients at the dominant scales

图4. 年径流深距平与主周期尺度下Morlet小波变换系数实部的变化过程

李家村站年径流第一主周期(20年)振幅在1956~1990年间相对稳定，随后明显减小，其余主周期振幅呈现先减小后增大的趋势，7年尺度的周期振荡逐渐占据主导；武山站第一主周期(22年)振幅变化与李家村类似，其余主周期振幅虽有所减小但在1990年后保持相对稳定，第四主周期(3年)在2015年后有增大趋势；郭城驿站第一主周期(20年)振幅相对稳定，第二主周期(8年)和第三主周期(3年)振幅也呈现先减小后增大的变化趋势；会宁站所有主周期振幅从1980年前后开始均明显减小。

总体来看，除郭城驿站年径流大尺度主周期振幅相对稳定外，其余三站都呈现减小趋势，2005年以后波动不明显；小尺度主周期振荡愈发剧烈，年径流呈现更加复杂的变化趋势。根据主周期尺度下小波系数实部的波动情况，预计未来2~3年内受水区内李家村、武山和郭城驿站年径流将仍呈现减小趋势；由于大尺度主周期振荡趋缓和小尺度主周期振荡加剧，大小周期的交替变化又隐藏嵌套在更大尺度的宏观结构中 [16] ，从不同时间和不同尺度上研究径流丰枯可能会得出不同的结论。

4.2. 径流年内分配特性及其变化趋势

各水文站实测月径流深的年内分布情况如图5所示。总体来看，四个水文站多年平均年内分布均呈“单峰型”。祖厉河流域内，郭城驿和会宁站年径流分配整体表现为较“尖瘦”的峰型，具有显著的季节性差异，年径流最低值基本出现在冬季，年内月径流峰值则基本出现在8月；洮河流域的李家村站与渭河流域的武山站则呈现相对“矮胖”的峰型，年内月径流峰值大多集中在7~10月，最低值则在1~2月，相对变幅较郭城驿和会宁两站明显更小。1960~2020年间，各水文站径流明显减小，在汛期尤为显著，径流年内分配出现“坦化”现象。

Figure 5. Inter-decadal variations of intra-annual distribution of station-measured runoff depth

图5. 不同年代际径流深年内分布

表2给出了各水文站汛期(5~10月)径流年内分配统计特征。总体来看，各水文站汛期月径流深下降趋势明显，Z值均通过1%显著性水平检验，其中会宁站最为显著，呈现明显的年代际差异。祖厉河流域内，郭城驿和会宁站汛期径流年内占比较其余两站明显更高。趋势检验显示，李家村和郭城驿站汛期径流年内占比的Z值分别为−3.45和−3.50，均通过1%显著性水平检验，呈显著减小趋势；武山站为−1.30，未通过5%显著性水平检验，减小趋势不显著；会宁站则为3.00，虽然会宁站汛期月径流表现出显著的下降趋势，但汛期径流年内占比却显著增大了，径流年内分配更加集中，这与会宁站近20年来非汛期出现断流现象有关。

Table 2. Distribution characteristics of station-measured runoff depth across inter-decadal flood seasons

表2. 不同年代际汛期径流深年内分配特性

图6进一步给出了年内分配的趋势检验结果。除郭城驿站外，其余三站年、季、月尺度径流深均呈减小趋势且大部分通过5%显著性水平检验，其中会宁站最为显著，会宁站又以冬季最为显著。虽然郭城驿站与会宁站同处祖厉河流域，且年径流深序列整体均呈现显著减小趋势，但位于干流的郭城驿站的河川径流年内分配出现夏秋两季更干、冬春两季更湿的变化特征，夏季径流深显著减小的同时，冬季径流深显著增大，径流年内变化渐趋平缓，图5(c)也映证了这一规律。李家村与武山站类似，除冬季径流深减小趋势不显著以外，其余三季径流深显著减小。

Figure 6. Trend test results of annual, seasonal, and monthly station-measured runoff depth: * indicates failure to pass the 5% notability test

图6. 年、季、月尺度径流深趋势检验结果：*表示未通过5%显著性水平检验

采用不均匀系数、集中度和集中期表征1957~2022年各水文站年内分配特征，变化过程如图7所示。

就年内分配的不均匀系数和集中度而言，祖厉河流域郭城驿和会宁两站比李家村和武山站更大，多年平均不均匀系数与集中度分别超过0.4和0.6 (表3)，且在2000年以前表现出较好的同步性，2000年以后会宁站的显著增大并表现出更强的年际波动，郭城驿站则相反。洮河流域李家村站和渭河流域武山站不均匀系数和集中度的大小和变化趋势在计算期内基本一致，多年均值约为0.24和0.37。图中不均匀系数与集中度的变化过程具有较高的一致性，但集中度有比不均匀系数更大的敏感性和分辨力，能更充分地表征径流年内分配的特点 [12] ；就集中期而言，各水文站最大月径流出现时间集中分布在6~8月之间，7月居多。以合成向量方向表示的集中期变化，显示李家村和武山站年内径流峰值出现的时间比郭城驿和会宁站稍晚一些，尤其是在2000~2010年间。2010年以后，除会宁站外，其余水文站径流年内分配的不均匀程度和集中程度减小，集中期波动范围更大。

Figure 7. Intra-annual distribution characteristics of station-measured runoff depth

图7. 径流深年内分配特征的变化过程

Table 3. Trend analysis results of intra-annual distribution characteristics of station-measured runoff depth

表3. 径流深年内分配特征的趋势分析结果

^*表示未通过5%显著性水平检验。

表3进一步给出了不均匀系数、集中度和集中期的趋势分析结果。1957~2022年间，除武山站外其余水文站不均匀系数和集中度的减小趋势显著，均通过5%显著性水平检验，郭城驿和会宁站通过了1%显著性水平检验，其中郭城驿站更为显著。集中期Z值均未通过5%显著性水平检验，没有明显的变化趋势。

5. 结论

1) 受水区内河川径流年际波动大，洮河流域干流李家村站集水面积和年径流深，明显大于其余水文站，位于祖厉河流域上游支流处的会宁站则面临断流风险；李家村、郭城驿、武山三个干流水文站年径流深的整体变化趋势一致，2000年以前受水区年径流深波动减小，2000年以后有所回升，会宁站则持续减小；各水文站实测年径流量均呈显著减小趋势，其中以会宁站最为显著；各水文站年径流深变化在年代际差异上存在较强跳跃性，最大概率突变年集中在1990年前后，突变发生后衰减趋势愈发显著。

2) 祖厉河流域两水文站实测年径流的波动能量较弱，周期特征较不明显，主周期尺度也较小；洮河流域李家村和渭河流域武山站年径流存在相似的周期特征，大尺度的振荡更强烈，并分别在20年和2年的第一主周期振荡下呈现较明显的“枯–丰–枯–丰–枯–丰–枯”演变特征；受水区年径流大尺度主周期振幅在2005年以后普遍减小，小尺度主周期振荡愈发剧烈，预计未来2~3年内各水文站年径流仍将继续减小；大小周期的交替变化及其与更大尺度宏观结构的隐藏嵌套，使得受水区径流年际变化不确定性加剧，从而呈现更复杂的变化态势。

3) 各水文站径流年内分配整体呈“单峰型”，祖厉河流域年内分配不均匀性和集中程度更大，集中期更早，月径流峰值则基本出现在8月，其余两站则大多集中在7~10月；1956年以来受水区汛期径流显著减小，除会宁站外各水文站汛期径流年内占比普遍减小；受水区径流年内分配不均匀程度和集中程度显著降低，集中期没有明显变化。

整体来看，引洮工程受水区径流显著减小，周期振荡愈发复杂，年内分配发生变化，水资源安全保障不确定性加大，亟需优化调控引洮工程供水与水资源配置。

基金项目

国家自然科学基金资助项目(52379022)；甘肃省重点研发计划(23YFFA0018, 23YFFA0017)；甘肃省水利厅2022年前期费项目(甘水规计发[2022] 331号)。