马玉婷, 蔡华阳, 杨昊, 刘锋, 陈欧, 谢荣耀,欧素英, 杨清书

1. 中山大学, 海洋工程与技术学院, 河口海岸研究所, 广东 广州 510275;

2. 河口水利技术国家地方联合工程实验室, 广东 广州 510275;

3. 广东省海岸与岛礁工程技术研究中心, 广东 广州 510275;

4. 南方海洋科学与工程广东省实验室(珠海), 广东 珠海 519000

河口是海陆相互作用的过渡区域, 受径流和潮流的相互作用。余水位(即潮平均水位)是河口径流、潮汐、地形等多因素耦合的结果, 其空间形态变化及演变机制是探讨河口陆海相互作用及潮波衰减的重要基础及前沿问题, 对探讨河口区水资源高效开发利用具有重要科学意义。河口区潮平均状态下,其水位坡降(水位空间变化率)近似与底摩擦项(或底切应力项)相平衡(Buschman et al, 2009; Cai et al,2015; Hoitink et al, 2016), 因此平均态下三角洲河道及河口区余水面坡降梯度(即余水位曲率)与河床底切应力梯度相平衡, 可指示底切应力梯度的时空演变, 并在一定程度上指示水道泥沙净输运的方向和地形冲淤趋势。研究河口余水位的空间变化尤其是余水位曲率的时空变化, 可从水动力角度揭示自然因素和人类活动共同影响下河床冲淤演变趋势,为防洪纳潮、水沙调控和航道整治开发等提供重要理论支撑。

密度梯度、潮汐不对称和回水效应的共同作用使得河口的平均水位向陆地方向抬升, 余水位及其坡降随流量增大而增大(Cai et al, 2013, 2014)。但余水位坡降随流量的增大不是单调递增, 而是存在阈值效应(张先毅 等, 2019)。低流量时河口过渡区呈凹形水面轮廓, 形成空间减速和沉积区, 而高流量时形成水流加速和侵蚀区(Lamb et al, 2012)。近年来,强人类活动对河口的干预日趋严重, 超过自然恢复能力, 已成为河口三角洲演变过程中的第三驱动力(陈吉余 等, 2008)。联围筑闸、滩涂围垦、河道采沙等强人类活动已经对三角洲河床演变产生重要影响。作为珠江主要泄洪通道的磨刀门河口, 近年来在挖沙等人类活动影响下上游河道河床平均加深约0.7m, 中下游河道河床加深约1.4m, 河道逐渐向窄深化方向发展(钱挹清, 2004; 郑国栋, 2005; 刘锋等, 2011)。河床下切导致河口上段水位显著下降、潮汐动力增强、咸潮入侵加剧等(Cai et al, 2019; 洪鹏锋 等, 2019)。

水位作为河口三角洲地区观测站位最多、时间序列最长且最连续的水文要素, 其在自然和人类活动影响下的时间演变特征及机制都得到较深入的研究。例如杨昊等(2020)、张先毅等(2020)对磨刀门河口平均水位变化及机制进行研究, 指出余水位在时间演变上存在着明显的季节变化和人类活动影响下的阶段性异变特征, 并采用边界驱动的回归模型对挖沙等强人类活动的影响进行定量辨识, 指出高强度采沙导致的河床下切使河口上段水位与水面坡降显著减小, 特别是夏季洪水流量时水面坡降减小最大。但针对余水位的空间分布及其演变特征、水位曲率变化对强人类活动的响应以及其与地形冲淤演变的关联性等问题仍有待进一步深入研究。本文聚焦磨刀门河口余水位空间分布及其演变特征, 基于长序列水位观测数据, 在探究余水位空间分布的基础上, 采用双变量线性回归模型定量辨识人类活动对余水位坡降和曲率的影响, 初步从水动力演变角度探讨三角洲河道冲淤演变趋势。

1 研究区域概况

珠江三角洲上承三江(西江、北江和东江), 下连八大口门(虎门、蕉门、洪奇门、横门、磨刀门、鸡啼门、虎跳门、崖门), 河网纵横, 洪潮叠加, 是世界上最为复杂的河网三角洲体系之一。西江干流进入珠江三角洲后于思贤 滘出现分流分沙, 其中思贤滘 西 滘口至百顷头河段称为西江干流水道(约80km),百顷头至口门河段称为磨刀门水道(约64km)。本文研究区域包括西江干流水道和磨刀门水道(以下简称为磨刀门河口), 即从上游的马口站延伸至下游口门附近的三灶站, 全程约144km (如图1 所示)。磨刀门河口为西北江河网的主要泄洪通道, 其多年输沙输水量居八大口门之首, 属于典型的河优型河口(吕海滨 等, 2006)。据马口水文站1965—2016 年月均流量资料统计, 其多年洪季平均流量为10688m3·s-1, 多年枯季平均流量为3583m3·s-1, 洪枯季差异显著。上游马口水文站多年月均水位为1.43m, 月平均潮差为0.34m, 下游三灶潮位站多年月均水位为-0.10m, 月平均潮差为1.09m, 口门至上游沿程水位抬升, 潮差减小。20 世纪80 年代中期开始, 珠江三角洲河网区出现大规模的河道挖沙现象。据前人统计, 西江干流水道1991—1997 年年均采沙量约为2.21×106m3, 采沙总深0.71m; 竹洲头至灯笼山1991—2000 年年均采沙量约为3.16×106m3,采沙总深1.06m; 年均采沙量约为年均来沙量的5倍, 引起河床高程不均匀下切0.7～2.3m, 河道逐渐向窄深化方向发展(钱挹清, 2004; 刘锋 等, 2011)。20 世纪60 年代中期起, 磨刀门口门地区开始大规模口门整治和围垦工程, 1966—1988 年建设白藤堵海大堤, 1984—1995 年建设横洲水道东西导堤和洪湾水道南北导堤, 口门区大片水域围垦成陆, 其水域面积减小约80%, 入海口门向外海延伸约16km,口门水位逐渐抬升(黎开志 等, 2005; 韩志远 等,2010; 胡煌昊 等, 2016)。磨刀门主干道疏浚工程于1999 年启动, 整治工程河段过洪断面面积加大, 附近水位略有下降(袁建国 等, 2009)。受地形边界(挖沙引起地形下切和口门围垦引起河道变窄以及疏浚引起断面面积加大等)和动力边界(海平面和上游流量变化等)的双重影响, 磨刀门河口水位呈现其独特的分布及阶段性变化。

图1 研究区域图及测站位置该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2019)4342 标准地图制作Fig. 1 Map of the study area and the stations used in this study

2 数据与方法

2.1 数据

本文收集马口、甘竹、江门、竹银、灯笼山和三灶6 个测站1965—2016 年月平均水位数据以及相应时段马口站的月均流量数据, 站位位置如图1 所示。所用数据来源于广东省水文局, 原始数据高程基面为冻结基面, 已统一校正至珠江基面高程, 基面数据缺失年份使用已知相邻年份的基面进行转换,可能导致水位有0～0.07m 的偏差, 相对于河段长度属于小量。

2.2 水面线沿程变化的数理统计

一维动量守恒方程中, 水位坡降项是河道水流运动的主要驱动力, 它是水位沿空间某个方向的变化率, 比如在一维河道中, 其水面坡降项是水位沿河道延伸方向, 即x方向(河道延伸方向)的一阶导数。但实际上由于河道地形、径潮动力的复杂性, 无法得到河口三角洲地区水位坡降(S)的理论解, 因此实际应用中主要采用其离散解来确定, 即:

由于在西江磨刀门河口144km 的范围内, 实测站点有限, 为更直观地显示水位空间分布, 假设三灶—马口之间的水位沿程变化可由此6 个站点(图2b*号)控制, 以三灶为起点, 设其x=0, 河道余水位(Z)沿河道距离x的变化可分别使用分段三次Hermite 插值、四阶多项式拟合、三次样条插值进行计算, 结果如图2b 所示(图中只为展示3 种方法的差异性, 其水位值非真实水位)。由图2b 可见, 分段三次Hermite 插值方法可避免过冲, 且可准确地连接平台区, 保证原始数据的准确性, 更能反映出相邻站位之间平均信息。因此, 本文主要采用分段三次Hermite 插值, 基于6 个站位的距离x与月平均余水位Z, 按空间间隔1km 对整个河口进行插值, 得出沿程 144 个断面的余水位Z, 并采用公式(1)和公式(2)分别计算沿程余水位坡降和余水位曲率, 用以定量分析月平均水位的空间分布及其变化规律。

图2 余水位曲率示意图(a)及水位沿程分布拟合方法对比(b)Fig. 2 Sketch of the residual water level curvature (a) and comparison of fitting methods for water level profiles along the estuary (b)

2.3 定量辨识强人类活动影响的双变量线性回归模型

在三角洲地区, 强人类活动(如围垦、大规模挖沙、河道疏浚等)导致三角洲地形发生显著变化, 进而引起径潮动力结构的改变。Cai 等(2019)、Yang等(2020)采用上下动力边界驱动的线性回归模型,定量辨识河道下切等对三角洲径潮动力的影响。本文在此基础上, 同样采用基于数据驱动的双变量线性回归模型, 以探讨水位空间结构的影响因素, 通过输入流量(Q)和下游边界水位(Zdown), 重构出河口沿程各站点的平均水位(Zsim), 回归模型的效果用均方根误差(RMSE)和相关系数(R2)来评价:

式中:α1、α2和β分别为流量、下游水位边界的待定系数和背景值待定系数。n为输入数据的长度,Zobs和Zsim分别表示实测水位和计算水位值, 上划线表示平均值。

人类活动对河口动力格局的影响是一个不断累积的过程, 累积到一定程度动力格局会发生异变,根据Yang 等(2020)、张先毅等(2020)的研究, 磨刀门河口动力结构对人类活动的响应主要集中在1990—2000 年期间, 因此将1965 年至1990 年(共26年)定义为缓变阶段(即强人类活动前), 2001 年至2016 年(共16 年)定义为强人类活动影响后的调整阶段。基于各站点强人类活动前的动力边界驱动回归模型及其重构的月平均水位, 计算出强人类活动前的沿程余水位坡降和余水位曲率; 使用缓变阶段的待定系数, 假设地形边界条件不变, 采用调整阶段的流量及水位边界驱动, 预报并计算调整阶段的余水位坡降和曲率。本文聚焦强人类活动对余水位曲率变化的影响, 强人类活动前后曲率总变化( ΔTOT)可分解为地形边界影响( ΔGEO, 即挖沙等强人类活动改变河道地形边界后的影响)和动力边界影响( ΔBOU, 即上游径流量和外海潮汐动力改变后的影响)(Cai et al, 2019; 杨昊 等, 2020), 通过计算 ΔGEO来定量辨识强人类活动对曲率的影响程度:

式中:Cobs和Csim分别为实测曲率值和计算曲率值,下标pre 和post 分别表示缓变阶段(强人类活动前阶段)和调整阶段(强人类活动后阶段)

2.4 余水位曲率与底床冲淤趋势的关联性

式中:α为待定系数。

基于一维泥沙守恒方程, 河道底床冲淤变化模型可描述为:

3 水位分布演变特征

从实测水位数据出发, 对沿程6 个站点月平均余水位进行三次Hermite 插值(1998 至2000 年江门站无水位数据, 故不进行插值), 计算出沿程余水位坡降与余水位曲率的数据序列, 结果如图3 所示。由图3a 可见, 余水位从口门至上游沿程抬升, 上游洪枯季变化较下游显著。由图3b 可见, 在江门至甘竹段余水位坡降存在负向波动, 且这种波动多在枯季出现, 而在强人类活动干预后, 甘竹、马口潮汐动力增强, 潮区界、潮流界显著上移, 这种负向波动则逐渐消失。余水位曲率的时空演变特征如图3c 所示。空间上, 曲率沿程有明显的正负值波动, 沿程水面线呈现凹凸形态转变; 而时间上, 每段曲率呈现波动的阶段性逆转, 与杨昊等(2020)对磨刀门河口人类活动划分的阶段对应较好。沿程余水位曲率的正负波动较人类活动前阶段发生显著变化, 三灶—灯笼山段曲率由负变正, 河道由侵蚀趋势转为淤积趋势; 灯笼山—竹银段曲率由正转负, 即河道由淤积趋势转为侵蚀趋势; 竹银—江门段的曲率由负转正, 河道由侵蚀趋势转为淤积趋势; 江门—马口段曲率由正转负, 河道由人类活动前的淤积趋势转为侵蚀趋势。

图3 磨刀门河口沿程余水位Z (a)、余水位坡降S (b)和余水位曲率C (c)的演变特征图中空白处由于1998 至2000 年江门站无水位数据, 故不进行插值。圆圈表示站位Fig. 3 Evolution of the residual water level (a), residual water level slope (b) and residual water level curvature (c) along the Modaomen estuary

为更直观的描述强人类活动前后水面坡降及曲率的季节变化及阶段性差异, 本文统计强人类活动前(1965—1990 年)与强人类活动后(2001—2016 年)多年月平均水位坡降与水位曲率, 结果如图4 所示。由图4 可见, 夏季水位坡降大, 冬季小, 上游坡降的季节性变化大于下游, 强人类活动后的调整期, 水位坡降均有不同程度的减小, 坡降的减小导致外海潮波传播所受的有效摩擦减弱、河口潮差增大、咸潮上溯加剧、余流结构改变等会进一步导致曲率发生变化甚至出现逆转。三灶至灯笼山年平均曲率减小0.13×10-4m·km-2, 河道倾向于全年淤积趋势; 灯笼山至竹银年平均曲率减小1.13×10-4m·km-2, 河道倾向于夏淤冬冲趋势; 竹银至江门段年平均曲率减小0.33×10-5m·km-2, 河道倾向于冲刷趋势; 江门至甘竹段年平均曲率增大1.30×10-4m·km-2, 河道倾向于夏冲冬淤趋势; 甘竹至马口段年平均曲率减小1.20×10-4m·km-2, 河道倾向于夏淤冬冲趋势。

图4 各区段缓变期(1965—1990 年)和调整期(2001— 2016 年)内多年月平均坡降S (a～e)与曲率C (f～j)变化Fig. 4 Monthly averaged variations in residual water level slope (a～e) and curvature (f～j) in the pre-human period(1965-1990) and post-human period (2001-2016)

4 水位空间分布影响因素

水位空间结构呈现阶段性变化, 主要由地形边界影响和动力边界影响共同作用, 本文采用上下动力边界驱动的线性回归模型, 计算得出沿程水位的预报结果, 从而定量辨识地形和动力边界对水位的影响程度, 进一步辨识出水位坡降和曲率的影响程度。模型率定计算结果(表1)和调整阶段预测结果(图 5)表明, 回归模型率定结果良好,缓变阶段(1965—1990 年)和调整阶段(2001—2016年)各站位RMSE 均在0.03～0.38m,R2均大于0.89,表明该模型能较好地模拟磨刀门河口的水位时空变化, 其误差主要取决于回归模型建立时的误差,如表1 所示, 误差范围在6%～11%, 江门站误差最大。利用缓变阶段的回归模型系数, 采用 2001—2016 年上下游动力边界来预估河道地形条件不变情况下的水位, 并计算实测与预测之差值, 此差值即为地形边界变化引起的水位变化量, 从图 5中可以看出, 2001—2016 年, 回归模型预测水位与实测水位有较大差异, 洪水期时, 其差值在0.2～2m, 远超出误差范围, 表明强人类活动对磨刀门河口水位变化有显著影响。

表1 磨刀门河口双变量三参数线性回归模型的率定结果Tab. 1 Calibration of the bivariate linear regression model with three parameters in the Modaomen estuary

图5 磨刀门河口缓变阶段回归模型的率定(a、c、e、g、i)和调整阶段水位变化过程的反演(b、d、f、h、j)a、b 为马口站; c、d 为甘竹站; e、f 为江门站; g、h 为竹银站; i、j 为灯笼山站。图中虚线为1:1 率定中线Fig. 5 Calibration of the regression model in the pre-human period (a, c, e, g, i) and reconstruction of the water level in the post-human period (b, d, f, h, j) in the Modaomen estuary

使用式(6)～(8)分别计算出水位、水位坡降、水位曲率的 ΔTOT、ΔGEO、ΔBOU沿程变化, 结果如图6和表2 所示。由图6a 可见, 强人类活动后, 全区段水位下降, 且上游较下游下降程度更明显, 上游流量输入减小和底床下切的叠加效应导致整体水位下降, 而下游灯笼山站挖沙活动引起的水位下降效应与口门围垦、海平面变化引起的水位抬升效应相抵消。不同程度的水位下降将反映在水位坡降的沿程变化上。由图6b 可见, 强人类活动后, 全区段水位坡降减小。三灶至竹银段水位坡降略有减小, 平均坡降减小0.22×10-2m·km-1; 竹银站附近及河口上游人为采沙剧烈, 导致水位坡降急速减小, 竹银至江门、甘竹至马口段平均水位坡降减小 0.9 6×10-2m·km-1; 江门至甘竹段坡降略有增大, 平均坡降增大0.001×10-2m·km-1, 整体上地形边界影响量ΔGEO与总变化值 ΔTOT趋势一致, 动力边界影响量ΔBOU沿程变化不明显, 对坡降的影响为负, 仅在口门处占总变化的比值较高。图6c 显示, 强人类活动后, 曲率的变化较为复杂, 口门段(三灶至竹银)和河口上段(甘竹至马口)水位曲率减小, 分别减小0.41×10-4m·km-2和1.04×10-4m·km-2; 河口中段(竹银至甘竹) 水位曲率增大, 平均曲率增大0.21×10-4m·km-2。地形边界影响量 ΔGEO与总变化值ΔTOT基本一致, 动力边界对曲率基本无影响, 表明水位曲率变化主要由强人类活动导致的河床地形变化引起, 而围垦、挖沙、疏浚等人类活动导致的水面坡降减缓、曲率正负转变是不可逆的。上述结果表明调整阶段后珠江磨刀门河口的冲淤演变、潮波传播等已经处于新的格局。

图6 磨刀门河口水位Z (a)、水位坡降S (b)、水位曲率C (c)变化值的沿程变化Fig. 6 Longitudinal variations in the residual water level (a), residual water level slope (b) and the residual water level curvature (c) along the Modaomen estuary

表2 强人类活动对沿程水位、坡降、曲率的影响程度Tab. 2 Influence of strong human interventions on the residual water level, slope and curvature along the estuary

通过前文对水位空间分布的分析, 可知径流向下游传播过程中存在一个坡降变缓的区域, 枯季时甚至出现倒坡降的现象。撇开地形变化对水面形态的影响, 探究仅由流量、外海海平面(以三灶站月平均水位作为代表)两个因子驱动对水位空间分布的影响。由图7 可见, 流量对上游水位的影响较大, 而外海海平面仅对下游附近水位的影响较大。正是由于流量、海平面高度与沿程水位的非线性耦合关系,引起坡降在空间上的非线性变化, 即从口门至上游,水位坡降沿程增大, 在距口门 80～100km 处减缓,100km 处后再增大, 其中流量越大, 水位坡降也随之增大。在下游水位项设置为0m 的条件下, 当流量低于4500m3·s–1时, 在中部河段(江门—甘竹)甚至会出现坡降为负的情况, 这是由于河流下泄受到潮流的顶托作用, 形成回水带, 表现为局部壅高现象,故坡降在此处会有所减缓。特别是在低流量条件下,全河道整体水位均较低, 潮汐动力相对强一些, 使得部分下游水位高于上游水位, 故形成倒坡降。坡降的沿程分布反映到曲率的变化上, 从口门至上游,曲率沿程变化为正-负-正, 随着流量的增大和海平面的上升, 曲率均没有发生逆转, 表明上下游动力边界驱动无法引起曲率的正负逆转。但在高强度的人类活动后, 水面坡降减缓, 潮汐动力增强, 其负向波动消失, 曲率出现正负逆转, 已然超出上下游动力边界对水面形态的驱动, 说明曲率的正负逆转是由地形变化反馈到动力结构上导致的。

图7 水位Z (a、b)、坡降S (c、d)、曲率C (e、f)随流量和下游水位变化的等值线图图c 中红色实线为对应坡降为0 的等值线Fig. 7 Contour plot of residual water level (a, b), slope (c, d) and curvature (e, f) as a function of river discharge and downstream water level. The red line is the contour with zero residual water level slope

5 讨论

基于韩志远等(2010)、刘锋等(2011)、张子昊(2018)、Liu 等(2019)和Yan 等(2020)对磨刀门河口地形演变的统计结果, 归纳曲率的正负与河道冲淤的对应关系, 结果如表3 所示。其中, 平均水深变化为正表示水深增大、河床呈冲刷趋势; 反之, 平均水深变化为负表示水深减小、河床呈淤积趋势。1962至1977 年间, 河道呈现中上段冲刷、口门段淤积,其中灯笼山至百顷头平均水深减小, 河槽容积因淤积减小15.2×106m3, 此河段C>0; 百顷头至马口段平均水深增大, 河槽容积因冲刷增大33.7×106m3,水面曲率C<0。从水位反演结果来看, 中下段通过余水位曲率可以较好的反演其冲淤趋势。1977 至1999 年间, 河道冲淤趋势较1977 年前发生空间上的逆转, 中下段及口门段(灯笼山至百顷头)受人类活动影响转为冲刷, 但水位曲率的调整并没有及时反馈, 此时水位曲率C>0, 指示河道仍呈现淤积趋势。1999 至2005 年间, 河道呈现中上段冲刷、口门段淤积, 其中灯笼山至竹银段平均水深减小0.06m,竹银至马口段平均水深增大0.59m; 全河道中下段河床演变趋势与曲率对应良好。2005 至2013 年间,河道呈现中段淤积、上下段冲刷, 其中百顷头至甘竹段平均水深减小0.02m, 其他河段平均水深均增大, 全段冲淤趋势与曲率反演结果对应良好。综上可见, 实测河道不同段的平均冲淤趋势与水位曲率变化基本一致。与曲率对应较差的河段基本是曲率为正, 而实测地形显示为冲刷趋势, 平均冲淤趋势与水位曲率对应不上的原因可能是: 1) 河床的冲淤受多种因素影响, 人类高强度挖沙活动使得整体河床向窄深化发展, 人为活动对河床的不均匀改造比水位反演来的更直接, 且水位反演有一定的滞后效应; 2) 公式(13)将平均水深作为检测冲淤变化的标准, 试想, 河宽变窄平均水深变深, 而过水断面面积变小的情况下, 实质上是淤积的情况, 可是平均水深变深被判定为冲刷, 冲刷和淤积还是不能仅单一的从平均水深变化来进行判断, 需结合过水断面面积; 3) 公式(13)的推导基于均匀流和平衡输沙假设得到, 而实际输沙并不一定是处于平衡输沙, 比如珠江洪水期来沙量变化、人类活动挖沙引起水体含沙量变化、地形不均匀下切引起的局部环流等因素导致河道非平衡输沙, 会对冲淤判断造成一定的误差, 在实际应用中有一定的局限性; 同时人类活动河道挖沙所形成的深坑, 尽管水深变大, 但坑内的水流是趋于静止, 余水位曲率不会有所影响。

表3 磨刀门河口不同河段曲率与冲淤趋势的对应关系Tab. 3 Correlation between curvature and tendency of erosion and deposition in different reaches of the Modaomen estuary

20 世纪60 年代中期起, 磨刀门口门地区先后建设白藤堵海大堤、横洲水道东西导堤和洪湾水道南北导堤, 口门区大片水域围垦成陆, 形成横洲水道和洪湾水道“一主一支”的入海格局, 加之联围筑闸和航道整治, 使得主流集中、流速增大、水流挟沙能力提高, 导致河床产生冲刷。20 世纪80 年代至21 世纪初, 大规模且高强度的河道采沙导致河道逐渐向窄深化发展。地形变化反馈到河道整体动力结构上, 口门围垦导致下游过水断面变窄, 且口门外海平面呈上升趋势, 水位空间分布为适应缩窄断面调整, 形成上凸形态, 下游曲率转为负值, 此空间动力结构的调整又驱使下游冲淤趋势发生逆转, 口门段由原来的淤积转为冲刷; 采沙使整体河床下切,航道整治扩宽过水断面, 水位空间分布为适应地形调整, 竹银至江门段水位坡降减小0.88×10-2m·km-1,而江门至甘竹段水位坡降增大0.17×10-2m·km-1, 形成江门下凹甘竹上凸形态, 高阻点向上游迁移, 中间段曲率转为正值, 上游段曲率转为负值, 此空间动力结构的调整又驱使冲淤趋势发生逆转, 中间段由原来的冲刷转为淤积趋势, 上段由原来的淤积转为冲刷趋势。人类活动导致地形变化反馈到水位空间分布上(即曲率), 水位空间分布又能进一步反演河道冲淤趋势, 基于长序列水位数据, 可提供一种简易预测河床冲淤演变趋势的方法, 为河口治理与水沙调控提供技术支撑。

6 结论

河口区域受径流、潮汐、地形等多因素的耦合影响, 其水面线形态的时空变化复杂, 余水位曲率既是刻画水面线沿程形态变化的参数, 亦是河床地形冲淤变化的指示因子。本文基于珠江磨刀门河口6 个站点1965—2016 年的月均余水位数据及相应时段马口站的月均流量数据, 通过数理统计方法和分段三次 Hermite 插值方法分析水位空间结构变化特征, 采用双变量线性回归模型, 定量辨识地形和动力边界对水位坡降和曲率的影响程度, 并结合地形数据, 初步从动力角度探讨潮平均态下曲率与底床冲淤趋势的关联性, 结果表明:

1) 在径潮动力影响下, 水位沿程分布具有明显的季节性波动。水位坡降与水位曲率洪枯季差异明显, 且枯季在江门—甘竹段水位坡降存在负向波动。强人类活动影响下余水位沿程分布呈现阶段性异变特征, 其中上段水位显著下降, 水面线坡降变缓, 江门—甘竹段负向波动逐渐消失, 各河段曲率出现正负逆转, 即从口门至上游, 水面线形态由下凹(C>0)-上凸(C<0)-下凹转变为上凸-下凹-上凸,河口冲淤趋势也随之改变。

2) 强人类活动引起河床地形下切是水面空间形态发生逆转的主要原因。河床下切引起余水位不均匀下降0.015～0.450m; 引起中上段(竹银至马口段)和下段( 三灶至竹银段) 余水位坡降分别下降0.46×10-2m·km-1和0.128×10-2m·km-1; 口门段(三灶至竹银段)和河口上段(甘竹至马口段)水位曲率分别减小0.41×10-4m·km-2和1.04×10-4m·km-2, 河口中段(竹银至甘竹段)水位曲率增大0.21×10-4m·km-2。上游流量和海平面高度的驱动仅能改变整体坡降的大小, 而地形异变是导致坡降沿程变化趋势改变的主要原因。

3) 水面线的凹凸性质(即曲率的正负)能较好地指示河床冲淤变化趋势, 河段平均曲率为正, 指示河床趋于淤积, 曲率为负, 河床趋于冲刷侵蚀。珠江磨刀门河口大规模且高强度人为活动导致的河床地形下切, 驱使水面坡降不均匀变缓, 水面曲率空间上发生逆转; 动力结构的调整驱使冲淤趋势发生逆转, 下段和上段由原来的淤积转为冲刷趋势, 中间段由原来的冲刷转为淤积趋势。基于长序列水位数据, 分析河口沿程水位分布特征及其阶段性变化,可探究人类活动、动力结构、地形变化之间的耦合反馈过程, 并提供一种简易预测河床冲淤演变趋势的方法, 为河口治理与水沙调控提供技术支撑。

需要指出的是, 本文发现强人类活动前枯季江门段月平均水位出现逆坡降, 但由于影响因素(比如潮非线性相互作用、江门水闸、冬季风)的不确定性,未对其形成机理进行深入分析, 这是本文的不足之处, 也是接下来需要通过珠江河网数值模拟来解决的重点问题。