1. 引言

河流携沙入海是陆海物质循环的关键环节。近几十年来，在高强度人类活动和全球性气候变化的影响下，河流水沙通量正在发生重大变化 [1] 。河流水沙变化带来的河床淤积和侵蚀将引起流域地表形态的演变，进而给流域防洪、供水、航运、生态治理带来新的问题与挑战。因此，科学认识河流水沙情势演变特征及其驱动因素，对流域未来开展径流管理和水旱灾害防御具有重要意义。

珠江是我国华南地区第一大河，年径流量约3360亿m³，年入海泥沙量约为6680万t。珠江主要由西江、北江和东江组成，西、北江思贤滘以下及东江石龙以下约26,820 km²的珠江三角洲是我国经济最发达地区之一。自1970年代以来，随着珠江流域经济社会状况的剧变，珠江三角洲地区的水沙情势发生了显著变化 [2] 。针对珠江流域的水沙情势演变，国内已开展了相关研究，对西江、北江、东江的水沙变化特征及造成水沙变化的因素进行了分析。在水沙序列的分析方法上，目前普遍采用的是定性分析法和定量分析法两类 [3] 。定性分析主要通过对水文序列进行简单处理，直观地反映其变化趋势和突变特征，常用的方法包括累计距平法、Hurst系数法 [4] 、滑动平均法等。定量分析则是基于一定的理论及指标，对序列变异和突变时间进行判定，常用的方法包括Spearman秩次相关法 [5] 、Mann-Kendall检验法 [6] 、Bayesian检验法 [7] 、滑动T检验法、滑动F检验法、R/S检验法等。现有研究大多采用上述分析方法中的一种或两种来对水沙序列的趋势及变异特征进行诊断，但由于不同方法在检验原理上各有侧重，其得出的分析结果往往存在差异，甚至互相矛盾。因此，有必要构建一套水文变异诊断体系来对珠江流域的水沙情势演变规律进行综合性分析。

本文针对河流水沙情势演变给珠江河口治理带来了新的不确定性，基于马口、三水及博罗3个水文站1959~2019年共61年的水文数据，采用水文变异诊断系统开展珠江河口三角洲水沙情势演变分析，从年际变化趋势、突变特征、周期性规律及其驱动因素上深刻理解河流水沙运移的特征，从而为流域河川径流管理和水旱灾害防御提供科学依据。

2. 材料与方法

2.1. 研究区概况

珠江由西江、北江、东江以及珠江三角洲诸河组成，流经云南、贵州、广西、广东、湖南和江西6省(自治区) 46个地(州)市、215个县及香港、澳门特别行政区，全长约2400 km，流域面积约45.37万km²。西、北江于佛山市三水区思贤滘交汇，而后汇入珠江河口三角洲，东江于东莞市石龙镇汇入三角洲，三江汇流后经虎门、蕉门、洪奇沥、横门、磨刀门、鸡啼门、虎跳门和崖门八大口门注入南海，形成“八口出海”的格局(图1)。

2.2. 数据来源

本研究选取马口、三水和博罗站1959~2019年共60年的水文数据，来分析珠江河口三角洲的水沙情势演变。马口、三水和博罗站分别为西江、北江和东江进入河口三角洲的主要控制站(图1)。上述各水文站点1959~2019年的年、月径流量和输沙量数据来源于珠江流域水文年鉴和中国河流泥沙公报。

2.3. 研究方法

本研究采用水文变异诊断系统对马口、三水、博罗站径流量、输沙量序列进行变异分析。水文变异诊断系统 [8] 包括三个阶段，分别为初步诊断、详细诊断和综合诊断(图2)。首先采用Hurst系数法对水文序列进行初步诊断，判别序列是否存在变异及其变异程度；若存在变异，再利用多种检验方法对水文序列的趋势及跳跃特征进行详细诊断；根据详细诊断结果和效率系数判别水文序列是趋势变异还是跳跃变异；最后结合实际水文调查结果，从成因上对序列的变异形式进行分析和确认，并形成最终的诊断结果。

2.3.1. 初步诊断

初步诊断旨在从整体上判断某一水文序列是否存在包含趋势、跳跃、周期等变异特征，以及初步评估变异序列的变异程度。Hurst系数是一种用于描述时间序列的长期相关性和自相似性的统计量 [4] ，已被广泛应用于水文、气象、信号处理等领域。Hurst系数h的取值范围为[0, 1]，其取值代表的含义为：当 $h=0.5$ 时，表明时间序列呈随机特性，即未来的变化趋势与过去无关；当 $h>0.5$ 时，表明时间序列的趋势具有持续性，即未来的变化趋势与过去的趋势一致；当 $h<0.5$ 时，表明时间序列的趋势呈反回归性，即未来的变化趋势与过去的趋势相反。Hurst系数可通过时间序列极差R与序列标准偏差S的比值R/S计算得到，见式(5)：

$R/S={\left(c\tau \right)}^h$ (1)

基于实测水文资料，采用最小二乘法可求得参数c和Hurst系数h。序列的变异程度(无变异、弱变异、中变异、强变异、巨变异)可结合Hurst系数h和某一显著性水平下的相关系数最大临界值 $r_{\alpha }$ 进行初步判断 [4] ，变异程度划分见表1。

2.3.2. 趋势检验

为定量判断某一水文序列是否具有增大、减小等趋势性变异规律，以及这种趋势变异是否显著，本文所建立的综合诊断系统将在初步诊断结果的基础上，对水文序列进行趋势变异分析。趋势变异分析的方法主要包括：相关系数检验法、Spearman检验法 [5] 、Mann-Kendal检验法 [6] 等。其中，相关系数检验法是评估两组变量之间的线性关系强度的常用方法之一；Spearman检验法主要针对水文数据的非线性特征，来描述两组变量在变大或变小的趋势上有多大程度保持一致；Mann-Kendall检验法首先假设获取的序列数据不具有趋势性，并构造呈标准正态分布的统计量Z，结合累计正态分布表可判断系列是否具有趋势性。

2.3.3. 跳跃检验

除趋势变异外，水文序列还可能存在序列值剧烈变化的情况，即跳跃变异。针对跳跃变异的分析方法主要包括滑动F检验、滑动T检验、滑动秩和检验、滑动游程检验法、有序聚类法、Lee-Heghinan检验法 [9] 、Brown-Forsythe检验法 [10] 、Mann-Kendall检验法 [6] 、Bayes方法 [7] 等。其中，Mann-Kendall检验法为具有代表性的非参数检验方法，该方法无需样本遵从一定的分布，并较少受到异常值的干扰 [4] 。对于样本序列 $X\left\{x_1,x_2,\cdots ,x_n\right\}$ ，Mann-Kendall检验法首先构造秩序列 $S_k$ ：

$S_k={\displaystyle \underset{i=1}{\overset{k}{\sum }}{r}_i}k=2,3,\cdots ,n$ (2)

其中：

$r_i=\{\begin{array}{l}+1x_i>x_j\\ 0x_i\le x_j\end{array}j=1,2,\cdots ,i$ (3)

基于上述时间序列定义统计量 $U{F}_k$ 如下：

$U{F}_k=\frac{S_k-E\left(S_k\right)}{\sqrt{Var\left(S_k\right)}}k=1,2,\cdots ,n$ (4)

式中， $U{F}_1=0$ ； $E\left(S_k\right)$ 为序列 $S_k$ 的均值； $Var\left(S_k\right)$ 为序列 $S_k$ 的方差。 $U{F}_k$ 为标准正态分布，对于给定显著性水平 $\alpha$ ，若 $\left|U{F}_k\right|>U_{\alpha }$ ，则表明序列存在明显的趋势变化。再对样本序列X的逆序 $\left\{x_n,x_{n-1},\cdots ,x_1\right\}$ 重复上述计算过程，得到统计量 $U{B}_k$ 。对于显著性水平 $\alpha$ 可得到相应的临界线，曲线 $U{F}_k$ 和 $U{B}_k$ 在临界线之间的交点即为突变开始的时刻。

2.3.4. 综合诊断分析

在趋势检验和跳跃检验的基础上，结合效率系数来最终判定水文序列是趋势变异，还是跳跃变异。若水文序列仅表现为单一的趋势变异或跳跃变异，则可以直接判定序列的变异形式，反之则需要计算序列趋势成分和跳跃成分的效率系数，以两者中的较大值来确定序列的变异形式。效率系数的计算式如下：

$R^2=1-\frac{{\displaystyle {\sum }_{i=1}^n{\left(Q_{obs,i}-Q_{sim,i}\right)}^2}}{{\displaystyle {\sum }_{i=1}^n{\left(Q_{obs,i}-{\bar{Q}}_{obs}\right)}^2}}$ (5)

式中， $Q_{obs,i}$ 为实测水文序列； ${\bar{Q}}_{obs}$ 为实测水文序列的均值； $Q_{sim,i}$ 为拟合水文序列：对于趋势变异， $Q_{sim,i}$ 为一元回归线上各点的值；对于跳跃变异，则有：

$Q_{sim,i}=\{\begin{array}{l}\frac{1}{k}{\displaystyle \underset{i=1}{\overset{k}{\sum }}{Q}_{obs,i}}i=1,2,\cdots ,k\\ \frac{1}{n-k}{\displaystyle \underset{i=k+1}{\overset{n}{\sum }}{Q}_{obs,i}}i=k+1,\cdots ,n\end{array}$ (6)

式中，k为跳跃点。

3. 珠江河口三角洲水沙情势演变规律

3.1. 水沙序列的多时间尺度变化特征

3.1.1. 水沙年际变化特征

根据1959~2019年径流及泥沙序列资料，绘制马口、三水、博罗三站的年径流量及输沙量变化曲线，如图3所示。在年际变化规律上，输沙量与径流量表现出一定的相关性，输沙量的波峰与波谷与径流量的丰枯趋势基本一致。总体来看，近60年珠江上游来水量的整体变化趋势不明显，但来沙量自1990年代起呈现显著下降趋势(图3(d))。1963年为有气象记录以来的广东省特大干旱年，其年径流量较多年平均偏少53.3%，进而引起年输沙量较多年平均偏少73.1%。近十年间(2010~2019年)，马口、三水、博罗三站的最大输沙量分别为2218万t、258万t、201万t，较多年平均输沙量分别减少了65.45%、65.53%、11.54%。

3.1.2. 水沙年段变化特征

为了更直观地反映珠江三角洲的水沙情势变化，图4统计了不同年代三水、马口、博罗三站的平均径流量和输沙量情况。由图可知，珠江三角洲不同年代间的径流总量相对稳定，而输沙量则从1990年代起迅速下降。

如图4(b)所示，1960~1969年间，珠江三角洲上游三站的年均来沙总量为8162万t；1970~1979年间，年均来沙总量达次大，为9130万t；1980~1989年间年均来沙总量微幅增加至最大，达9225万t，此后，便进入持续下降阶段。1990~1999年开始，珠江三角洲年均来沙总量大幅下降，减小至7169万t；2000~2009年间，珠江三角洲年均来沙总量仅为3089万t，约为统计年段最大值的33%；2010~2019年间，年均来沙总量继续下降至2271万t，仅为统计年段最大值的25%。

3.2. 水沙序列变异诊断结果

3.2.1. 初步诊断

本文采用水文变异诊断系统对马口、三水、博罗站径流量、输沙量序列进行变异分析，以了解三角洲的水

沙变异规律。在初步诊断阶段，马口、三水、博罗三站的年径流量序列Hurst系数h分别为0.679、0.878、0.612 (见表2)。在显著性水平 $\alpha =0.05$ 和 $\beta =0.01$ 下，资料长度为60年时(自由度为58)的相关系数临界值分别为0.254和0.330，对应的Hurst系数 $h_{\alpha }=0.663$ ， $h_{\beta }=0.706$ 。根据表1可判断，马口、三水、博罗的年径流量序列分别弱变异、强变异，无变异。同样地，由马口、三水、博罗三站年输沙量序列Hurst系数可判别其变异程度分别为中变异、弱变异、无变异。以上初步诊断的结果仅表明水沙序列可能存在趋势或跳跃变异，但具体变异形式和变异发生时间还需结合详细诊断进一步分析。

注：表中+表示跳跃或趋势显著，−表示跳跃或趋势不显著，/表示不能进行显著性检验；第一显著性水平α = 0.05，第二显著性水平β = 0.01；↑↓分别表示趋势上升、下降。

3.2.2. 详细诊断

采用水文变异诊断系统中3种趋势检验法和7种跳跃检验法对马口站、三水站和博罗站的年径流量、年输沙量序列做进一步分析。在显著性水平 $\alpha =0.05$ 下，三站的水沙序列详细诊断结果如表2所示。可以看出，在趋势变异上，3种诊断方法得到的检验结果一致，其中，三水站水沙序列趋势变异显著；马口站年径流序列无显著趋势变异，而年输沙序列存在一定趋势变异；博罗站的水沙序列均无显著趋势变异。但在跳跃变异上，不同诊断方法得出的跳跃变异点并不一致，因此，还需结合实际的水沙序列判断跳跃变异的发生年份。

以MK检验为例，图5~7分别为马口、三水、博罗三站的MK检验曲线，在置信度区间内，UF曲线和UB曲线的交点即为发生变异的年份。如图5(a)所示，马口站年径流量的UF、UB曲线相交于1985年，且该点之后UF值持续小于0，表明1985年之后，马口站的年径流量为下降趋势，但由于UF值处于置信水平负区间[−1.96, 0]内，说明马口站年径流量在1985年后的下降趋势不显著。马口站年输沙量的MK检验结果(图5(b))显示，1999年为输沙序列的跳跃变异点，且该点以后的UF值小于0并超过0.01置信水平区间，表明马口站的年输沙量在变异点后呈显著的持续下降趋势。

三水站年径流序列的UF、UB曲线在置信水平区间内相交于1989年(图6(a))，从该站的年径流变化情况(图3(b))来看，1990年代以后，三水站的年径流量明显提高，即MK检验得到的跳跃变异点基本正确。同理，结合MK检验结果(图6(b))和三水站的实际年输沙量变化情况可知，三水站的年输沙量大致在2000年左右发生跳跃变异。

博罗站年径流量序列和年输沙量序列的UF曲线和UB曲线在置信水平区间内均有多个交点(图7)，表明其水沙序列可能存在多个变异点。但由于博罗站年径流量序列的MK检验曲线并未超过0.05置信水平区间，因此可认为博罗站年径流量序列无跳跃变异点。

3.2.3. 综合诊断

由上述详细诊断结果(表2)可知，不同诊断方法得到的跳跃变异点可能不同，即使同种诊断方法也可能得到多个跳跃变异点。但哪些变异点才是真正的变异点，哪个变异点的变异程度更大？还需要通过综合诊断做进一步判断。此外，对于同时存在趋势变异和跳跃变异的水沙序列，需计算序列中趋势成分和跳跃成分的效率系数，以判定序列最终的变异形式。

水沙序列的跳跃变异点根据该点被诊断为变异点的次数来确定，即马口站年径流序列发生跳跃变异的年份为1983年，年输沙序列发生跳跃变异的年份为1997年；三水站年径流序列发生跳跃变异的年份为1992年，年输沙序列发生跳跃变异的年份为1998年；博罗站年径流序列无明显跳跃变异点，年输沙序列发生跳跃变异的年份为1987年。

基于确定的跳跃变异点，进一步可得各水沙序列趋势变异和跳跃变异的效率系数(见表2)，最终确定马口站年径流和年输沙序列的变异形式分别为1983年跳跃变异和1997年跳跃变异；三水站年径流和年输沙序列的变异形式分别为上升趋势变异和下降趋势变异；博罗站年径流序列无明显变异特征，年输沙序列为1987年跳跃变异。

3.3. 西、北江水沙分配比变化

西、北江来水在思贤滘连通，经天然平衡调节后重新分配进入珠江河口三角洲。在西、北江水沙情势发生变化后，两江交汇后的水沙分配比也随之改变。

在分流比方面，马口站年径流量占思贤滘总径流量的比例在1990年代以后呈现出下降趋势(图8(a))，该站在1989年前后的年均分流比分别为85.8%和78.9%，后期分流比下降约8.1%，而三水站在此期间的分流比相应增加了8.1%。三水站分流比增加意味着进入三角洲腹地的洪水总量将会增加，若遭遇潮汐顶托，则可能导致三角洲腹地滞洪瓶颈河段上游洪水位异常壅高，使三角洲地区的防洪安全面临严峻压力。

在分沙比方面，自1990年代起，马口、三水两站的分沙比波动明显(图8(b))，但整体来看，随着三水站分流比增大、河流挟沙量增加，三水站的分沙比由1989年以前的9.9%上升至12.7%，分沙比增幅为28.3%，而马口站的分沙比相应降低28.3%。三水站分沙比增加使得北江的输沙量增大，进而加剧北江河道淤积，导致航道宽度变窄，这将加大广州出海航道的维护难度。

3.4. 八大口门水沙分配比变化

西、北江水沙分配比的改变将引起珠江河口八大口门水沙分配比的变化。总体来看，东四口门的分流比高于西四口门，但近年来东四口门的分流比有所减小，而西四口门的分流比增大。如表3所示，2010~2018年间，注入伶仃洋东四口门的年径流量占珠江河口总年径流量的58.9%，较1985~2000年间降低了3.6%。从八大口门的分沙比来看(表4)，东四口门的分沙比呈增加趋势，由1990年代的54%上升至近期的56.8%；而西四口门的分沙比则由1990年代的46%减小至近期的43.3%。

八大口门水沙分配比的变化将导致河口水动力和泥沙特征发生改变，尤其是东四口门分沙比增加后，伶仃洋悬沙浓度整体升高，平均增加了0.02~0.05 kg/m³，珠江河口区悬沙高浓度范围将明显向东扩展。

4. 珠江河口水沙情势演变成因分析

珠江河口水沙情势的变化与流域大规模的人类活动密切相关，本节将结合马口、三水、博罗三站水沙序列的变异诊断结果与上游水利工程建设、流域水土保持工作等，对珠江河口三角洲水沙情势演变的成因展开分析。

4.1. 水利工程

1960年代起，珠江流域陆续进入了大规模的水库建设时期，水库蓄水直接导致泥沙被拦截在库，从而引起下游输沙量下降。图9~11分别展示了马口站、三水站、博罗站的水沙相关关系，结合三站所在的西江、北江、东江流域的重点水库建成时间可以看出：对于西江流域，早期建成的大王滩、西津等水库的拦水拦沙作用较为有限，而1990年代以后建成的天生桥、龙滩等水库因其规模较大，对西江径流量和输沙量都造成了显著影响，使得西江年输沙量在1990年代中后期出现变异；对于北江流域，1999年建成的飞来峡水库是导致北江输沙量下降的重要原因；对于东江流域，其水沙情势演变相对平稳，由图11可以看出，新丰江、枫树坝、白盆珠三座水库的修建并未明显改变水沙关系。

4.2. 流域采砂及水土保持工作

珠江河口区水沙的变化除了受河道水利工程的影响外，还与流域内采砂活动及水土保持工作的开展密切相关。1990年代初至2006年，北江流域采砂活动频繁，大规模采砂导致北江河网河床下切明显 [11] ，进而引起北江径流规律发生较大改变(见图3(b))，1990年以后的北江年径流量较1990年以前增加约53.9%；与年径流量相比，北江流域的年输沙量在1990年代末期才发生明显变异，1998年后三水站的年输沙量较此前显著降低。从时间节点上看，北江流域输沙量骤降除了与兴建水库、大规模采砂有关外，还与广东省于2000年左右开展的水土保持工作有关 [11] ；此外，2003~2004年间，北江流域年均降雨量仅为1369 mm左右，较邻近年份偏少约30%，降雨的减少也进一步加剧了北江上游来沙量的降低。

1970年代至1980年代，西江流域森林乱砍乱伐导致流域水土流失严重，进而造成这一阶段的河流输沙量呈增加趋势(见图4(b))。据统计，西江流域在20世纪七八十年代的总体水土流失面积近2000 km²，对此，西江流域自1990年代起开始实施水土保持系列措施，使广西、广东两省(区)的森林覆盖率分别由1980年代的22%、27.7%提升至2005年的52.71%、55.9%。水土保持在一定程度上减少了河流的输沙量，因此自1990年代中后期，西江的年输沙量呈下降趋势(见图4(b))。

1980年代至1990年代，在国家水土保持政策的引导下，东江流域植被覆盖率较治理前提高约1.2倍，年入河泥沙减少约525万t，受此影响，东江博罗站的年输沙序列在1987年出现变异(见表2)，年输沙量呈现出下降趋势。另外，2000年前后博罗站附近的大量采砂活动也在一定程度上造成博罗站年输沙量下降，但其对博罗输沙序列的整体演变趋势的影响较为有限。

5. 结论

本文基于珠江河口区马口、三水、博罗3个水文站的实测水沙资料，采用水文变异诊断系统，耦合多种检验方法，对珠江河口三角洲的水沙情势变化进行分析，得出主要结论如下：

1) 水沙年际变化趋势显示，1959~2019年间，珠江上游来水量基本保持不变，但珠江上游输沙量自1990年代起则显著减少，特别是马口站和三水站，2010~2019年间的最大输沙量较多年平均值分别减少了65.45%和65.53%，输沙量的变异程度明显高于径流量的变异程度。

2) 根据水文序列变异诊断结果，马口站的年径流和年输沙序列分别在1983年和1997年发生跳跃变异；三水站的年径流和年输沙序列的主要变异形式分别为上升趋势变异和下降趋势变异；博罗站年径流序列无明显变异特征，年输沙序列在1987年发生跳跃变异。

3) 在西、北江分流比方面，1990年代以来，经北江三水站进入三角洲河网区的年径流量和年输沙量比例上升，增幅分别为8.1%和28.3%。三水站分流比、分沙比增大导致进入三角洲腹地的洪水总量增加，并使得北江输沙量增加，这将加大三角洲地区的防洪压力及广州出海航道的维护难度。

4) 西、北江分流比和分沙比改变，使珠江河口八大口门水沙分配比出现相应的调整。近年来注入伶仃洋东四口门的分流比略有降低，但东四口门的分沙比由1990年代的54%增加至近期的56.8%。

5) 流域上游拦河筑坝、大规模采砂、水土保持等人类活动是导致珠江河口输沙量变化的主要原因，河道水沙不平衡将可能引起三角洲河道沿程冲刷下切、河口海岸带侵蚀后退等不良现象。

基金项目

国家重点研发计划项目(2021YFC3001000)。

参考文献