河床下切对珠江三角洲峰值水位演变的影响

2022-07-27郑金海鲍仕昱季小梅

河海大学学报（自然科学版） 2022年4期

郑金海，鲍仕昱, ，张蔚，季小梅

(1.河海大学海岸灾害及防护教育部重点实验室，江苏南京 210098；2.河海大学江苏省海岸海洋资源开发与环境安全重点实验室，江苏南京 210098)

河口三角洲是世界上经济最发达、人口最稠密的地区之一[1]。径流潮汐共同作用产生的非线性效应，给三角洲河网带来高度动态和复杂的环境[2]。随着社会和经济的高速发展，日益频繁的人类活动对河口三角洲地区长期动力环境造成的变化受到国内外学者的广泛关注[3]。人类活动，例如航道整治、河道采沙、通航疏浚等活动，直接改变了河口三角洲地貌形态，造成了河床下切。

水深的增加是造成河口水动力扰动的主要因素[4]。河床下切减小了河道底摩擦，从而可能导致潮汐过程的变化，并与风暴潮、盐水入侵等灾害有关。同时，河口地区的水位变化也与河床下切密切相关。已有的研究表明，河床下切后，在上游区域会出现同流量下水位下降的现象：Luo等[5]分析了珠江挖沙前后的数据资料，发现在相同流量下三水站2005年比1989年的水位下降1.28 m；Ralston等[4]通过数值模拟得到了哈德逊河口河床下切后洪水总体风险降低的结论。这一现象也在湄公河、哥伦比亚河[6]出现过。但珠江三角洲在经历河床下切后，中腹部地区却出现明显的壅水现象[7-8]。以往学者对三角洲水位的研究大多采用年极值水位，但是年极值水位通常包含风暴潮等极端事件[9]，同时会受潮汐调制周期的影响。因此，本文通过长时间序列诊断，剔除极端事件与潮汐调制的影响，深入分析珠江三角洲河网在河床下切背景下峰值水位的时空演变规律。

1 研究区域与方法

1.1 研究区域

珠江由3条主要支流组成：北江、东江和西江。这些支流入海，形成了珠江三角洲。该区域位于中国南海岸，北与广州接壤，西南与澳门接壤，东南与香港接壤，在21°40′N～23°N和112°E～113°20′E之间(图1)，面积约17 200 km2。珠江三角洲潮汐是以半日潮为主的混合型，平均潮差为1.0～1.7 m[10]。北江、东江和西江通过8个出口流入南海。西面有4个口门，分别为崖门、虎跳门、鸡啼门、磨刀门，其余4个口门(横门、洪奇门、角门、虎门)位于东侧。

图1 珠江河网和水文站的位置分布Fig.1 Channel network and locations of hydrologic stations in the Pearl River Delta

20世纪80年代中期以后，珠江三角洲的人类活动随着城市的发展逐渐加剧。20世纪90年代初到1995年，河道采沙达到顶峰，西江、北江和东江河床平均下切深度分别为0.59～1.73 m、0.34～4.43 m和1.77～6.48 m。这种不受控制的采沙加速了河床下切，导致整个河网的水动力发生改变。2000年开始，政府开始关注河床下切造成的不良影响，并制定法规和政策控制采沙，使之放缓。

1.2 数据来源

对珠江三角洲西江和北江干流的8个水文站1961—2012年的水文数据进行分析。由于马骝洲站数据长度不足，并且珠江西四口门动力特征相似，本文用黄金站替换马骝洲站，水文站的位置见图1。由于观测技术的限制，原始数据序列的每个全日潮周期(24.84 h)只有2个高水位和2个低水位。用三角插值法对水位数据进行1h间隔的填充，用于小波分解。水位数据来源于《中华人民共和国水文年鉴》。

1.3 概率密度函数法

人类活动影响的极值水位的演变一直受到学者们的关注。同时估算极值水位的方法也随之发展，包括经典的年极值法、r最大值法、蒙特卡罗法、经验模拟法、联合概率法。然而，这些方法都不能分离人类活动对水位趋势的影响。潮汐、径流、风暴潮是影响水位趋势最主要的因素[11]。风暴潮由于其出现频率较低影响较大，影响了河床下切对极值水位的长周期趋势分析。本文根据Vellinga等[12]提出的概率密度函数法定义峰值水位，以便进行趋势分析。这种方法可以有效地降低风暴潮和热带气旋等极端事件对峰值水位的影响，适用于受强烈人类活动影响的河口三角洲。

概率密度函数的方法是根据连续一年的高低潮位绘制水位的概率密度曲线(图2)。曲线用2个峰值和1个极小值分为4个区间。在每个区间内，将频率与相应的水位乘积积分得到4个质心(从左到右为形心1、形心2、形心3和形心4)。形心1为概率密度函数法确定的峰值低水位(PLWL)，形心4为概率密度函数法确定的峰值高水位(PHWL)。

图2 水位概率密度曲线Fig.2 Probability density function curve of water level

1.4 S-Tide方法

潮汐是海水在太阳和月球的引力作用下产生的涨落运动，表现为垂直周期性运动。潮汐的涨落周期约为12 h或24 h，同时潮汐还有半月、月、年以及18.61 a的长周期变化。月球绕地轨道与赤道成一定夹角，月球的近地点和远地点的连线每8.85 a转一周，上升节点与下降节点位置每18.61a向西绕一周，因此主要的天文潮的振幅和相位在长变化中处于18.61a的节点调制[13]。

节点调制是潮汐对水位的主要长周期调制。为了避免节点调制对水位趋势分析的影响，本文采用S-Tide方法[14]。S-Tide方法是一种将经验模态分解(empirical mode decomposition, EMD)引入河流潮汐分解的非稳态潮汐分析方法。该方法无须预先设定基函数，而是根据自身的时间尺度进行信号分解。这种特征使经验模态分解方法在理论上可以应用于所有类型信号的分解。

S-Tide方法与NS-Tide方法类似，均假定分潮的振幅和迟角随时间而变化：

(1)

式中：Z(t)为t时刻观测到的水位；S(t)为t时刻的平均水位；J为分析的分潮总个数；aj和bj的平方和是第j个分潮对应振幅的平方；σj是第j个分潮的角速度。

S-Tide通过基于样条函数方法的独立点方案求解式(1)。独立点的个数和使用的分潮是S-Tide最核心的2个输入参数。独立点个数越多，S-Tide得到的分潮振幅和迟角的变化就越复杂；独立点越少, 得到的振幅和迟角的变化就越简单。当独立点个数为1时，S-Tide得到的分潮振幅和迟角都是不随时间变化的。基于此将S-Tide用于水位和潮汐振幅数据中18.61 a节点调制的提取，并去除峰值水位时间序列中18.61a节点调制的影响。

2 结果与讨论

2.1 18.61 a节点调制对峰值水位的影响

图3展示了从水位及分潮簇振幅时间序列中提取18.61 a的节点调制。周期为全日的分潮簇(D1)和周期为半日的分潮簇(D2)潮汐振幅通过连续小波分解(CWT)得到。结果表明，越靠近上游的站点受到潮汐节点调制越小，越靠近海洋的站点潮汐节点调制越大。这是因为潮汐在三角洲河道网中向上游传播的过程中受到河流流量、浅水地形等多种混杂因素的影响，节点调制逐渐减小。总体来看，D1和D2分潮簇都存在18.61a的节点调制，这是由天文潮引起的。由于高频率的潮汐对径流等造成的阻尼更加敏感，在传递过程中衰减更快，因此D1分潮簇的节点调制比D2分潮簇更加明显。分潮簇振幅的节点变化要比水位中大得多，这是因为每个分潮簇都有不同的相位[15]。当不同分潮簇叠加时，会削弱整体节点调制的变化。因此，虽然节点调制对潮汐振幅的变化有明显影响，但对基于概率密度函数方法的水位指标变化的趋势影响有限，且这种影响越向上游越小。

图3 PHWL和潮汐振幅移除18.61 a节点调制的影响过程Fig.3 Comparison of PHWL and tidal amplitudes before and after removing the nodal cycle of 18.61 a

2.2 峰值水位趋势

Zhang等[16]基于Pettitt检验，指出在河床下切影响下，珠江三角洲西江与北江水位趋势的突变点主要集中于20世纪80年代中期。这表明河床下切对水位的影响大约始于1985年。本文将以此为分界点，将数据年限分为采沙前期(1961—1985年)和采沙期(1986—2012年)，并比较采沙前和采沙期的水位变化趋势。基于概率密度函数法定义峰值水位，并将该峰值水位应用于西江与北江干流的8个主要水文站。对采沙前后各站点的年峰值水位进行线性趋势分析(图4)，图4中细线表示95%的置信区间,MWL(mean water level)表示平均水位。

图4 西北江干流去除节点调制的水位参数的变化Fig.4 Variation of modulated water level metrics by removing nodal influence at the main channel of the West River and the North River

图4比较了1961—1985年(采砂前)和1986—2012年(采砂期)的去掉节点变化的水位变化趋势，结果表明由概率密度函数计算得到的质心参数的变化趋势因地区而异。对靠海水文站而言，采沙前北江入海口门(万顷沙西站)PHWL的增加趋势比西江入海口门(黄金站)明显，这可能与伶仃湾向陆方向的宽度收敛有关。对于内陆站，潮差的减小将双峰概率密度曲线转化为单峰概率密度曲线，概率密度函数可以仅用2个质心(PLWL和PHWL)作为参数描述。

珠江三角洲水位趋势呈现时间和空间上的变化(图4)。采沙前，由于联围筑闸的影响，三角洲上游的MWL呈逐渐上升趋势[17]。其他水位参数也会随着MWL的增加而增加。MWL在靠海的水文站以每年大约0.24 cm的速度上升，与南海平均海平面上升速率(0.27 cm/a)相似[18]。

在河床下切的影响下，不同区域的峰值水位趋势变化不同。在1980年后的一段时间内，由于频繁的河床采沙与航道疏浚，珠江三角洲上游站的MWL、PLWL和PHWL呈显著下降趋势，尤其是三水站和马口站水位下降最为明显。值得注意的是，PHWL在三角洲中腹部仍然显示出上升的趋势，这意味着不仅高水位在上升，而且极端水位也在上升。竹银站MWL和PLWL呈下降趋势，三善滘站MWL呈上升趋势，PLWL呈下降趋势，表明MWL的变化并没有直接转化为高水位和低水位的变化。在河床下切后，下游站点的水位参数均呈现上升趋势，这与三角洲口门处的土地围垦等人类活动有关[19]。

图5为西江和北江水文站的峰值水位累计频率曲线。为了确定河床下切后的峰值水位变化，基于概率密度函数法计算了采沙前和采沙期PHWL的重现期。为了消除潮汐节点调制对峰值水位的影响，本次分析取采沙前1961—1980年和采沙期1993—2012年的水位序列。结果表明，在上游(马口和三水站)，河床下切后10 a重现期和50 a重现期的峰值水位比河床下切前减少了1～2 m，说明河床下切导致三角洲上游峰值水位降低。与上游相反，河床下切后中腹部区域10 a和50 a重现期的峰值水位比河床下切前高0.1～0.2 m，表明河床下切增加了中腹部区域的洪水风险。该区域峰值水位上升可能是潮汐顶托径流，阻碍了径流下泄。虽然挖沙活动造成了河道河床下切，但同时也增大了河道体积，使得纳潮量增加。增大的纳潮量不仅可以防止中腹部地区的水位随河床降低而下降，还可以减缓径流速度，而导致了峰值水位的升高[20]。在内陆河交汇口，由于干、支流的交汇存在水体间的相互顶托作用，造成局部水流壅塞[21]。在下游，峰值水位上升，这是由于河床下切导致潮汐振幅增大，增大的潮汐振幅和海平面上升的综合效应导致了潮控区峰值水位的增加。