珠江口磨刀门水道盐度变化对水文气象要素的频率响应特征

陈子燊，高时友，2，陈玲舫

（1.中山大学 水资源与环境系，广东 广州，510275；2.珠江水利科学研究院，广东 广州，510611；3.三亚深海科学与工程研究所，海南 三亚，572031）

全球气候变化对珠江三角洲水文过程的影响十分明显，非汛期珠江流域降水少径流不足导致河口区咸潮上溯。近十几年来咸潮持续时间增长，上溯影响范围趋大已严重影响了珠江三角洲城镇饮水安全和农业灌溉用水的需求。其中，磨刀门水道是珠江流域的主要入海通道之一，是中山、珠海、澳门三地主要淡水源地，枯季持续的咸潮入侵使得沿岸许多取水口氯度超标，由此引起的淡水资源短缺对当地经济社会发展造成了十分不利的影响。

近年来，研究人员采用水动力模型分析磨刀门水道咸潮输移特征。其中，陈文龙等（2014） 分析了磨刀门水道盐淡水分层与混合特征、盐份物质输移机制，探讨了咸潮上溯强度时空分布差异的原因；Wang 等（2012） 和Gong 等（2011） 分别探讨了咸潮上溯规律及其与径流量和潮差的响应关系；Zhou 等（2012） 通过敏感性数值试验指出如海径流量是盐水入侵的控制因素，盛行冬季风可增强咸潮入侵。不少研究则从观测资料的动力学角度分析了枯季磨刀门河口盐度变化规律。陈荣力等（2012） 通过分析磨刀门水道咸潮上溯实测资料，分析了该水道咸潮运动规律及上溯机理，认为潮汐动力与咸潮运动密切相关，径流是抑制咸潮上溯的主要动力，咸淡水混合状态随潮差而变化。贾良文等(2006) 根据磨刀门枯季水文观测资料,对河口咸淡水混合等水文特征进行了初步研究，认为在枯季由于径流较弱,潮流成为主要动力。包芸等（2009）、刘杰斌等（2008） 利用磨刀门水道天河至挂定角之间8 个水厂水闸的表层盐度数据，通过拉格朗日插值方法作磨刀门水道近3 个月的纵向逐时咸界分布图，给出了0.5‰，2‰，5‰和8‰咸界的逐时变化情况，分析了其运动规律。吕爱琴等（2006） 通过分析认为磨刀门水道的咸潮上溯成因主要为枯季上游来水量减少、潮动力增强、河床演变、潮位的改变。陈水森等（2007） 导出含氯度与径流、潮流、河口地形等的关系式，建立了咸潮入侵的经验模型。宋晓飞等（2014） 分析指出珠江口磨刀门盐水入侵最直接的影响因素是径潮相互作用，磨刀门的盐水入侵活动规律复杂多变。

河口盐度变化过程也可视为河道地形制约下上游入海径流的冲淡作用、进入河口的潮波潮流涨落对咸潮进退的控制和风力对河口纵向净环流的扰动等共同作用下的多因素驱动复杂的过程—频率响应系统。入海径流和涨落潮对河口不同位置水体盐度变化的真实影响存在于不同的频率，盐度对二者驱动的真实响应频率需要在扣除潮周振荡作用后加以进一步识别。此外，枯季偏北风驱动表层冲淡水的向海输运可加强近底层高盐度海水的上溯作用，其对河口盐度的混合扩散也不可忽视。因此，枯季河口盐度变化和咸潮入侵范围是集径流冲淡、涨落潮与风应力驱动与混合作用等主要因素于一体的综合问题。本文拟对珠江口磨刀门水道沿岸盐度变化对水文气象环境要素的频率响应作经验诊断分析，此将有助于进一步深入认识河海气相互作用下的河口咸潮入侵过程与盐淡水混合作用。

1 研究环境与基本数据

1.1 研究环境水文特征

磨刀门是珠江流域西江水系主要入海口门，年径流量约923 亿m3，占珠江入海总径流量的28.3%，属于径强潮弱的河优型河口。但径流量年内分布不均，主要集中于每年的汛期，非汛期入海径流明显减小。马口水文站位于西北江三角洲顶部，是西江流入珠江三角洲的国家水文控制站，西江水主要经西江下游的磨刀门出海（图1）。监测显示，非汛期涨潮流可上溯至马口站上游的西江干流，是磨刀门水道盐水入侵的主要因素。根据F=（O1+K1） /M2划分潮汐类型，磨刀门潮性系数F 为1.54，属于不正规半日潮类型。三灶站平均潮差为1.11 m，大横琴站为0.94 m，竹银站0.80 m。进入磨刀门向河道上游传播过程的潮汐动力持续递减。

图1 西江下游磨刀门水道沿岸测站位置

1.2 基本数据

采用2010年11月1日0 点-2011年3月31日23 点马口水文站流量、三灶站潮位、竹银站的风速风向与挂定角、联石湾、竹排沙、平岗、竹银5 个氯度测站长度为5 个月151 天共3 624 h 的逐时同步观测数据。各测站的位置见图1。对氯度（Cl） 数据作了以下处理：按S‰ = 1.805 0Cl +0.030 的关系式把氯度转换为盐度（S）。针对磨刀门水道走向，对竹银站的风速风向做了坐标轴旋转，采用西北—东南向作为磨刀门水道的主流向风速。分别建立5 个测站盐度对三灶站潮位、马口站流量、竹银站主流向风速的多变量频率响应系统。5 个测站的盐度、马口站流量、三灶站潮汐和竹银站主流向风速过程见图2、图3。

2 研究方法

为了确定多变量驱动系统（x1，x2，…xn） 中任意两个变量间的相关是直接或间接的，需要把双变量常相干分析扩展到多变量偏相干分析。偏相干分析理论就是计算模型中的条件输入和偏相干函数。由于模型中的条件输入已经去除了其他输入的相干部分信息,所以偏相干函数能有效表示各输入对输出的独立影响。采用Marcos（2007） 提出的由功率谱密度矩阵求解偏相干函数的计算方法，简述如下：

图2 马口流量、三灶潮位和竹银风速过程

图3 磨刀门水道5 个测站盐度

设有多维时间序列{x1t，x2t，…xnt；t∈T}，计算的交叉谱矩阵可写成分块矩阵形式：

式中：

当i≠j 时，上式中Crij（ω）是{xit} 和{xjt} 之间的交叉谱，即：Crij（ω）= Pij（ω）+ iQij（ω）。式中，Pij（ω）为协谱；Qij（ω）为正交谱；ω 为角频率。当i =j 时，Crij（ω）是{xjt} 的功率谱。

由式（1） 计算偏交叉谱阵：

阵中元素分别表示{x1t} 和{x2t} 之间的偏谱和偏交叉谱。由式（2） 计算{x1t} 和{x2t} 之间的偏相干谱为：由式（2） 计算{x1t} 和{x2t} 之间的偏相干谱和偏位相谱分别为：

式中：Cr12k（ω）、φ12k（ω）和PCoh12k（ω）代表了时间{x1t}， {x2t} 序列，在排除了其它序列影响下的偏交叉谱、偏位相谱和偏相干谱。采用Halliday 等（1995） 提出的信度水平α 下的临界偏相干函数值：

式中：α=0.05；υ =2L；L：样本分段数；L0：表示余变量的维数。常相干谱的估计方法可参见黄嘉佑（1984） 而不再列出。

3 计算结果与经验诊断分析

采用MATLAB 编程计算各测站氯度对三灶站潮汐、马口站径流和竹银站主流向风速的频率响应。首先对各输入输出变量去除趋势和标准化处理，随之分别采用常相干谱和偏相干谱估计方法计算盐度对三灶站潮汐、马口站径流和竹银站风速的频率响应特征。按式（5） 计算的临界偏相干系数为0.488。三灶站潮汐、马口站流量和竹银站主流向风速功率谱见表1 和图4，表2 为流量-潮位的常相干谱估计提取的最大4 个峰值（临界常相干系数为0.327），各测站盐度对潮汐、流量、风速响应的最大4 个偏相干系数及相应的频率见表3 和图5，表2 中加粗的数字表示为大于临界值的偏相干系数（为省篇幅，略去常相干谱计算结果）。

图4 三灶潮汐、马口流量、竹银主流向风速的功率谱

功率谱计算结果表明，潮汐与径流普遍存在准1日、半日、1/3日、1/4日和1/6日等振荡周期，而竹银站的主流向风速则显示出典型的随机振荡特征。

流量与潮位的常相干谱的最大4 个峰值分别代表口门外三灶站潮汐与上游马口站流量之间存在的准1日、半日、1/4日和1/3日耦合振荡的显著周期，不同频带的位相差正负值反映流量对潮汐调制的响应存在时间超前或滞后。

挂定角、联石湾、灯笼山、平岗、竹银5 个测站的盐度分别对三灶站潮汐与马口站流量的常相干谱分析显示出盐度对潮汐或流量的响应都显示出相同的潮周振荡特征（图5）。参考表2 和图4 易见，磨刀门沿岸各测站盐度对上游马口站流量和口门外三灶站的潮周响应实际上只是向上游传播的潮汐调制结果，并未排除入海径流对潮汐和盐度的直接影响。

表1 三灶站潮汐、马口站径流和竹银站主流向风速功率谱计算结果

表2 流量-潮位常相干谱最大4 个峰值

表3 各测站盐度对三灶潮汐、马口流量、竹银主流向风速响应的偏相干估计结果

5 个盐度测站的盐度变化与三灶站潮汐、上游马口站径流和竹银站主流向风速的偏相干响应具有以下特征：

（1） 盐度变化与潮汐主要振荡周期不完全一致。沿磨刀门水道向上游，无论是显著偏相干响应频带数或显著偏相干系数值都反映潮汐对盐度变化的影响递减，竹银站盐度已不存在对潮汐作用的显著响应频带。换言之，竹银站盐度可视为潮汐作用影响下咸潮上溯的平均上限位置，此对于在竹银修建的抽淡与应急供水水库提供了一定的科学依据。此外，盐度变化响应潮汐过程存在的多个非潮周期是否与亚潮振荡相关有待进一步深入研究。

（2） 经分离潮汐直接影响后的盐度与马口站流量之间依然存在准1日、半日和1/4日等振荡周期的显著偏相干性，反映了受到潮汐调制后的流量对盐度混合变化的间接影响。

图5 各测站盐度对潮汐、流量、风速响应的相干性

（3） 尽管风的变化随机性很强，但枯季偏北风可进一步驱使表层落潮阶段的冲淡水向海流动和底层入侵海水的上溯作用，从而加强河口区纵向净环流作用。对比偏相干谱和常相干谱分析结果可看出，在消除潮汐和径流对盐度直接影响后，盐度变化对风的响应明显加强。其中，竹银站和平岗站的盐度存在对竹银站的主流向风速的两个显著偏相干性频带。其余3 个监测站距离竹银站加大，不存在盐度对竹银站主流向风速的显著频率响应，此反映出磨刀门沿岸盐度监测站位置受风的影响程度的不一致性，要真实反映不同测站盐度变化对风扰动的实际响应需要采用测站所在地的测风资料。

（4） 偏位相谱估计反映盐度变化对潮汐、流量与风的响应过程存在时间超前或滞后，但不存在普遍的规律。此进一步说明研究期间磨刀门水道不同测站位置的盐度混合过程对潮汐、流量以及风的响应的不一致性和盐度混合作用的复杂性，河口咸淡水混合过程是多个水文气象要素共同作用的结果。

4 主要结论

本文对2010年11月-2011年3月磨刀门水道挂定角、联石湾、竹排沙、平岗、竹银5 个测站的盐度对上游马口水文站流量、口门外三灶站潮位和竹银站的主流向风速的频率响应的经验诊断分析，获得以下认识：

（1） 入海径流存在准日潮、半日、1/4日和1/3日等显著潮致振荡周期；

（2） 各测站盐度变化与潮汐主要振荡周期并不一致，沿磨刀门水道向上游，潮汐对盐度变化的影响递减，竹银站盐度可视为潮汐直接作用下咸潮上溯的平均上限位置。

（3） 消除潮汐和径流的影响后风对盐度变化的作用明显加强；

（4） 各测站盐度变化对潮汐、流量与风的响应过程存在时间超前或滞后但无特定规律，说明磨刀门水道盐度变化过程十分复杂。

Gong W,Shen J,2011. The response of salt intrusion to Changes in river discharge and tidal mixing during the dryseason in the Modaomen Estuary,China.Continental Shelf Research,31(7/8）：769-788.

Halliday D M, Rosenberg J R, Amjad A M, et al. 1995. A framework for the analysis of mixed time series/point process data- Theory and application to the study of physiological tremor, single motor unit discharges and electromyograms, Prog. Biophys. Molec. Biol. 64,237-278.

Underwood M A,KELL ER T,2007.Using the Spectral Density Matrix to Determine Ordinary, Partial, and Multiple Coherence. Spectral Dynamics,9(2）：11-20.

Wang B,Zhu J,Wu H,et al,2012.Dynamics of saltwater intrusion in the Modaomen Waterway of the Pearl River Estuary. Science China :Earth Sciences,55(11）：1 901-1 918.

Zhou W, Wang D, Lou L, 2012. Investigation of saltwater intrusion and salinity stratification in winter of2007/2008 in the Zhujiang River Estuary in China.Scientia Geographica Sinica,31(3）：31-46.

包芸，刘杰斌，任杰，等，2009.磨刀门水道盐水强烈上溯规律和动力机制研究.中国科学（G 辑） .39（10）：1 527-1 534.

陈荣力，刘诚，高时友，2012.磨刀门水道枯季咸潮上溯规律分析.水动力学研究与进展A 辑.26（3）：312-317.

陈水森，方立刚，李宏丽，等，2007.珠江口咸潮入侵分析与经验模型-以磨刀门水道为例.水科学进展，18（5）：751-755.

陈文龙，邹华志，董延军，2014.磨刀门水道咸潮上溯动力特性分析.水科学进展， 5. http://www.cnki.net/kcms/ detail/ 32.1309.P.20140806.0837.018.html.

黄嘉佑，1984.气象中的谱分析.气象出版社.

贾良文，吴超羽，任杰，等，2006.珠江口磨刀门枯季水文特征及河口动力过程.水科学进展.17（1）：82-88.

刘杰斌，包芸，2008.磨刀门水道枯季盐水入侵咸界运动规律研究.中山大学学报（自然科学版） .47（增刊2）：122-125.

吕爱琴，杜文印，2006.磨刀门水道咸潮上溯成因分析.广东水利水电，5：50-53.

宋晓飞，石荣贵，孙羚晏，等，2014.珠江口磨刀门盐水入侵的现状与成因分析.海洋通报，33（7）：7-15.