水体溶解氧(DO)是水生生命体赖以生存的核心元素和必要条件^[1-2]，也是反映河流、湖库等水体环境质量状况与自净能力的重要指标，参与水体到水体物质转化过程^[3-4]，对于维持水生生态系统健康状态有着重要的意义^[5-6]。水中的DO含量取决于耗氧及复氧过程及平衡条件^[7-11]。当水体的耗氧速率＞复氧速率时，DO值开始降低，进而威胁水生生态系统稳定^[12]。因此，研究水体DO的时空分布特征及其影响因素，对掌握水体生物多样性和水质保护具有重要意义^[13-14]。

国内关于DO的研究大多集中在富营养化的河流、水库、湖泊及近海海域^[15-23]，鲜见有关大型贫营养水体的研究。新丰江水库是全国第一个“中国优质饮用水资源开发基地”和东江水资源核心调节枢纽，也是东江饮用水源的“枢纽”，对粤港澳大湾区(简称“大湾区”)的社会经济发展起着至关重要的基础支撑作用^[24]。根据历年河源市环境质量公报，新丰江水库水质稳定保持在国家《地表水环境质量标准》(GB 3838—2002) Ⅰ类水质标准(DO≥7.5 mg/L)，其综合营养状态指数[TLI(∑)]评价处于贫营养状态[TLI(∑)≤26]^[25]。在全球变暖和人类活动不断加剧的大背景下，社会经济发展和人类生产生活方式转变，城镇生活源排放、畜禽养殖和农村面源污染等，都对新丰江水库的水质造成了严重威胁^[26-27]。目前，对水库的研究主要集中在浮游植物、浮游动物和沉积物等方面^[28-31]，对水库表层DO特征进行详细研究的鲜见报道。近年来，水库出现DO值低于地表水Ⅰ类水质标准限值状况，对水库的科学管理和正常运营带来了困扰。因此，亟需深入分析新丰江水库DO的变化特征及其影响因素，从而更加客观地评估贫营养水体的水生态环境质量。

本研究以新丰江水库为例，对DO的时空分布进行研究，并结合小波变换、最大信息系数和多元线性回归等多种方法量化分析各因素对DO的影响，以期为水库的科学管理与保护提供支撑，切实保障供水安全，为大湾区建设和可持续发展提供有力支持。

1 研究方法 1.1 研究区域概况

新丰江水库位于广东省河源市和韶关市(23°41′15″—24°7′45″N，114°19′30″—114°45′40″E)，始建于1958年，是全国第八大水库和华南第一大水库，流域面积5 813 km²，以丘陵和山地为主。水域面积370 km²，库容1.39×10¹⁰ m³，平均年入库水量0.6×10¹⁰ m³，多年平均水深28.7 m，水力滞留时间2 a，库区属南亚热带季风气候，年平均气温20.8 ℃^[24]。

1.2 样品采集及分析

在新丰江水库中心、左右两边和库区出入口共设置6个采样点位(a—f)，具体点位分布见图 1。

于2018年4月—2022年11月(共56个月)，以每月1次的频次进行样品采集，采用卡盖式水样采样器(青岛聚创环保设备有限公司)在水深0.5 m处采样，采用YSI EXO2 multiparameter sonde型深水多参数水质分析仪(美国赛莱默公司)现场测定水样的DO，水温(WT)，电导率(EC)和pH值，总氮(TN)、总磷(TP)、高锰酸盐指数(I_Mn)的实验室分析方法参照《国家地表水环境质量监测网监测任务作业指导书(试行)》^[32]。

1.3 分析方法

采用折线图、箱型图和柱状图分析DO值的月度变化、年际变化和空间变化。采用连续小波变换提取DO的周期，采用最大信息系数、交叉小波变换、相关系数和多元线性回归方法揭示DO时空变化的驱动因素。

1.3.1 连续小波变换

连续小波变换可以将时间序列分解为时域和频域2个维度，进而深入获取更多数据意义，已被广泛应用于数据的周期分析^[33]。计算公式如下：

$\mathrm{WT}(a, b)=\frac{1}{\sqrt{a}} \int x(t) \varphi\left(\frac{t-b}{a}\right) \mathrm{d} t$

(1)

式中：WT——小波系数；a——尺度参数，决定了信号的压缩程度；b——位移参数，决定了信号的平移程度；x(t)——原始信号序列；φ(t)——小波基函数；dt——积分变量。

基于现有研究成果^[34-35]，采用Morlet母小波(小波基函数)的实部分析DO的周期性^[7]，通过计算不同周期对应的小波方差，确定变化主周期。小波方差越大，代表其对应的周期越显著。小波方差的计算公式为：

$\mathrm{WV}=\int \mathrm{WT}(a, b)^2 \mathrm{~d} t$

(2)

式中：WV——小波方差；a——尺度参数；b——位移参数；dt——积分变量。

1.3.2 交叉小波变换

交叉小波变换能够同时分析时域和频域中2个序列之间多时间尺度的相关性关系^[36]，其计算公式为：

$W_n^{X Y}(s)=W_n^X(s) W_n^{Y *}(s)$

(3)

式中：W_n^XY(s)——交叉小波功率谱；W_n^X(s)——序列1的连续小波变换；W_n^Y*(s)——序列2的连续小波变换的复共轭；X——序列1{x₁，x₂，…，x_n}；Y——序列2{y₁，y₂，…，y_n}；n——参数变量；s——小波相位角的标准偏差。

交叉小波只能揭示序列2在高能量区的位相关系，对于低能量区位相关系解析不足，而相干小波可以弥补这一点，小波相干谱定义为^[37]：

$R_n^2(n)=\frac{\left|S\left(s^{-1} W_n^{X Y}(s)\right)\right|^2}{S\left(s^{-1}\left|W_n^X[s]\right|^2\right) \cdot S\left(s^{-1}\left|W_n^Y[s]\right|^2\right)}$

(4)

式中：R_n²(n)——相干性系数；S——平滑算子；s——小波相位角的标准偏差。

交叉小波能量谱和小波相干谱图中黑色细实线为小波边界效应影响锥，粗黑实线表示置信水平为95%以上区间，代表两者显著性。

1.3.3 最大信息系数

最大信息系数(MIC)可以有效地度量变量间的各种复杂关系，已经在许多领域得到了广泛应用^[38-39]，MIC若超过0.30，则可以认为变量之间存在相关关系^[7]。一般将最大信息系数和皮尔逊决定系数(R²)连用以量化变量间的线性关系和非线性关系，计算公式如下：

$\alpha=\mathrm{MIC}-R^2$

(5)

式中：α——非线性系数；MIC——最大信息系数；R²——皮尔逊决定系数。α＞0.2时，表示变量间的非线性关系较强，反之则表示变量间的线性关系较强。

1.3.4 多元线性回归

DO的时空变化受到多重因素的复合影响，选用多元线性回归来定量不同因素的贡献。计算公式如下：

$Y=a+b_1 X_1+\cdots+b_n X_n$

(6)

式中：Y——DO质量浓度，mg/L；X——水温等影响因素；a——截距；b——X和Y之间的回归系数；n——参数变量。

1.4 数据处理

以Excel 2019、SPSS 22.0、Matlab软件进行数据处理和统计学分析；以Matlab 2019软件进行最大信息系数、连续小波和交叉小波分析；以Origin pro 8.0、ArcGIS 10.6软件绘制统计图。

2 结果与分析 2.1 DO时空变化规律 2.1.1 DO月度变化

2018—2022年各采样点位DO值月度变化趋势见图 2(a)—(f)。从图 2可见，2018—2022年DO月平均值均呈现出先降低后升高的变化趋势，DO月平均值为6.94~9.48 mg/L，均能达到《GB 3838—2002》Ⅰ类或Ⅱ类水质标准值(6.0 mg/L)。DO值一般从2—4月开始降低，直到9—10月达到最小值，之后逐渐上升，年内整体变化呈现“V”字型。以2021年为例，4个点位最大值出现在3月，2个点位最大值出现2月。5个点位最小值出现在9月。时间上，水库56组DO月平均值监测数据中，有8组数值低于《GB 3838—2002》 Ⅰ类水质标准(DO＜7.5 mg/L)，占比14.3%。空间上，a、b、c、d、e和f点位的56组DO月平均值数据中，分别有4，9，11，8，17，14次低于Ⅰ类水质标准，分别占比7.1%，16.1%，19.6%，14.3%，25%和30.4%。其中，a点位占比最低，f点位占比最高。

2.1.2 DO年际变化规律

各点位DO值年际变化趋势见图 3。

由图 3可见，除a点位外，其他各点位的DO年均值变化趋势基本一致，年际间存在一定变化幅度，均可达到Ⅰ类标准限值(DO≥7.5 mg/L)。a点位DO值在2018年为7.91 mg/L，2019年上升至8.71 mg/L，然后逐年下降，在2022年下降至7.87 mg/L。而其他点位DO年平均值分别从2018年的7.88(b点位)，7.81(c点位)，8.01(d点位)，7.84(e点位)和7.95 mg/L(f点位)，在2019年上升为8.53(b点位)，8.028(c点位)，8.38(d点位)，8.43(e点位)和8.36 mg/L(f点位)，然后在2020和2021年持续下降至7.74(b点位)，7.8(c点位)，7.85(d点位)，7.56(e点位)和7.71 mg/L(f点位)，在2022年DO值再次上升至8.48(b点位)，8.34(c点位)，8.45(d点位)，8.07(e点位)和8.15 mg/L(f点位)。

2.1.3 DO空间变化规律

2018—2022年各采样点位DO值箱型图见图 4。由图 4可见，2018—2022年a, b, c, d, e, f点位的DO年平均值分别为8.16，8.14，8.07，8.15，7.95和8.02 mg/L，均能达到《GB 3838—2002》Ⅰ类水质标准。

2.1.4 DO周期性

为了进一步揭示DO在不同时间尺度上的变化规律，采用连续小波转换提取了DO的周期及其对应的小波方差，具体见图 5(a)—(f)和表 1。由图 5可见，水库的DO在整个时域范围内出现了11~12和35~37个月2个震荡周期(以颜色的变化表示周期性)，a、c、d和f点位对应的震荡周期均为12和36个月，b点位对应的震荡周期为12和37个月，e点位对应的震荡周期为11和35个月，说明水库DO存在显著的年周期。此外，在部分时域内，水库点位DO还存在3~6个月的微弱震荡周期，尽管其对应的小波方差较小(表 1)。

2.2 DO时空变化与环境因子关系分析 2.2.1 最大信息系数结果

调查期间，各点位DO与WT之间的MIC均最高(表 2)，MIC为0.41~0.53，均＞0.30，说明WT是水库DO时空变化最重要的影响因素。其中，仅c点位的DO与WT之间为线性负相关关系(α=0.18)，其他点位的DO与WT均为非线性关系(α＞0.2)。f和e点位pH值的MIC均为0.35，说明pH值在一定程度上影响这2个点位DO值的变化，但这2个点位的DO值与pH值均为非线性关系(α＞0.2)。c和f点位EC的MIC分别为0.34和0.38，说明EC在一定程度上影响这2个点位DO值的变化，但其DO值与EC均为非线性关系(α＞0.2)。b点位的I_Mn的MIC为0.32，说明I_Mn在一定程度上影响该点位DO值的变化，但其点位DO值与I_Mn呈非线性关系(α＞0.2)。

2.2.2 交叉小波转换结果

为了进一步揭示DO和影响因素的响应关系，采用交叉小波转换对原始数据进行分析，结果见图 6(a)—(r)和图 7(a)—(r)。图中颜色代表相关性，黄色越深，相关性越高；向左表示负相关关系，即变量之间变化趋势相反，向右表示正相关关系，即变量之间变化趋势相同，向上代表自变量超前因变量的变化，向下代表自变量滞后因变量的变化。水库各点位的WT与DO在全时间尺度的8—16个月周期上均存在显著负相关关系。a，b，c，d和f点位的pH值与DO在10~14个月周期上存在显著相关关系。其中，a点位pH值与DO呈负相关，f点位pH值与DO呈正相关，b、c和d点位的pH值变化滞后DO。e点位的pH值与DO在2~5个月周期上存在显著正相关关系。b，c，d，f和XH点位的EC与DO在9~15个月周期上存在显著负相关关系，所有点位在2~5个月周期上均存在一定的正相关关系。I_Mn与DO的相关性与EC类似。各点位TP与DO的周期相关性差异较大，其中a点位TP与DO在10~13个月周期上存在显著负相关关系。c，d，f和e点位的TP与DO分别在10~12，9~14，7~9和11~12个月周期上呈显著正相关，同时f点位的TP与DO在11~15个月周期上也存在显著正相关关系。水库b，c，d和f点位的TN与DO在8~11个月周期上存在显著相关关系，其中b和c点位的TN值变化滞后DO，d和f点位呈正相关。e点位的TN与DO在6~8个月周期上存在显著正相关关系；各点位的DO和TN在1~4个月周期上存在一定相关关系，相关关系有所差异。

2.2.3 多元线性回归

为了量化主要影响因素对DO时空变化的贡献，采用多元线性回归方法来对DO进行拟合。根据最大信息系数和交叉小波转换的结果，本研究首先采用单一WT指标拟合DO，在此基础上叠加MIC在0.30以上的其他因素拟合DO，最后输入所有影响因素拟合DO，结果见表 3。水库各点位的WT可解释DO整体变化的量并不高，最高仅为26%(c点位)，最低为14%(a和e点位)，叠加其他因素后相关性系数提升明显，最高提升了30%(e点位)，说明WT并不是影响DO的唯一因素。其中，叠加pH值相关性系数提升最为明显，最高提升了22%(b点位)，说明pH值是影响DO的次要主因；而叠加I_Mn，TN和TP后，相关性系数提升并不明显，提升幅度为2%~11%，e点位相对较高，说明该点位除WT和pH值外，I_Mn、TN和TP的综合作用对DO变化贡献也很明显。

3 讨论 3.1 新丰江水库DO时空变化特征

水体DO水平受有机物降解、生物呼吸的耗氧及空气中氧溶解、水生植物光合作用复氧的双重作用，这2种作用的消长使水中的DO呈现出时空变化^[40]。新丰江水库DO总体呈“V”字型变化，夏末秋初最低(9~10月)，DO是决定水库月度水质达Ⅰ类或Ⅱ类标准的关键因子。DO和WT呈现极显著负相关(P＜0.01)，说明水库DO值符合随水温的升高而降低的物理规律，与泸沽湖^{[14, 41]}、千岛湖^[42]的研究结果相似。水库点位间的DO值月度变化存在一定差异，两侧位置(f和e点位)的DO值为Ⅱ类标准比例明显高于水库出口点位(a点位)，f和e点位的DO与pH值和TP呈正相关，a点位与pH值无相关性，与TP呈负相关，说明两侧区域受雨季有机物降解耗氧影响大，与千岛湖^[42]研究成果相似。

风力搅动会对上层水体DO的水平分布造成影响^[43]。新丰江水库DO平均值的空间差异性小，各点位的DO变化主周期介于11~12个月和35~37个月，点位间DO最大均值和最小均值差值仅为0.2 mg/L，但从水库上游入库口至下游出库口DO值呈现缓慢递减的趋势，在长潭水库也有类似的现象出现^[44]。a点位DO年均值有别于其他点位，尤其在2022年DO值不升反降，比最高点位低了7.2%。其原因在于a点位于水库出库口，空间细长狭窄，营养盐浓度长期处于极低水平，藻类活性受磷浓度水平限制^[29]，而其他点位处于相对开阔水域，风力、风速与出库口点位明显不同，且雨季营养盐浓度相对较高，说明水库营养盐状况和环境特征决定了DO空间分布，与千岛湖^[42]、抚仙湖^[43]的研究结果相似。

3.2 新丰江水库DO的影响因素

淡水中氧的溶解度主要取决于水温^[45]，水体中氧的饱合度随水温升高而降低是基本物理规律。最大信息系数、交叉小波转换和多元线性回归等多种研究方法结果表明，水温是水库DO时空变化的主要影响因素，与千岛湖^[42]、抚仙湖^[43]的研究结果相似。浮游植物光合作用释放氧气的同时吸收二氧化碳，使得水中碳酸氢根大量积累，引起pH值升高^[46]。水库各点位pH值的MIC均相对较高，其中f和e点位pH值的MIC达到0.35；交叉小波转换结果显示，pH值与DO在2—5个月和10—14个月周期上存在显著相关关系；WT和pH值组合与DO的多元线性回归结果明显提升了相关性系数，说明pH值也是新丰江水库DO值变化的重要因素。氮磷等营养盐浓度作为影响浮游藻类季节变化的关键因子^[29]，控制着藻类增殖强度，从而影响水体DO水平。新丰江水库I_Mn、TN和TP与DO均呈线性关系，但MIC均＜0.3，且WT、pH值、EC、I_Mn、TN和TP组合与DO的相关性系数提升并不明显(介于2%~8%之间)，说明I_Mn、TN和TP并不是影响水库DO的重要因素，贫营养水体营养盐对DO变化贡献相对较弱，与洱海^[3]有所不同。

总体而言，DO作为影响新丰江水库水质状况的关键指标，水库DO的时空变化主要受自然要素WT影响，其次为人类活动产生的营养物质。就人类活动而言，过去几十年来新丰江水库流域经济社会发展迅速，增加了水库的污染负荷，尤其是生活、农田等点源和面源会带来大量的耗氧物质，影响水库DO水平，进而对水库水生态系统造成一定影响。因此，首先应加强水库集水区污染源头管控，有效控制农村生活污水、畜禽养殖、农业面源等污染强度，降低入库负荷。其次改进生产方式，推广科学施肥、种养结合、循环利用等方式，从源头减少污染物的排放。最后建立完善的监测网络，构建涵盖水平和垂直方向的DO及水生态等相关指标相结合的综合监测体系，持续观测DO变化可能引发的生态效应。

4 结论

(1) 新丰江水库作为典型的南亚热带贫营养型水库，2018—2022年水库水体的DO年均值为7.95~8.16 mg/L，均能达到《GB 3838—2002》Ⅰ类水质标准，DO月均值为6.94~9.48 mg/L，均能达到Ⅰ类或Ⅱ类水质标准。水库的DO空间差异性较小，时间上总体呈“V”字型变化，夏末秋初(9—10月)DO值最低。

(2) 新丰江水库DO变化主周期时间相对较长。连续小波转换分析结果表明，水库DO主要在11~12个月和35~37个月的年周期下循环变化。部分时域内，水库DO还存在3—6个月的变化周期。

(3) 基于最大信息系数、交叉小波变换、相关系数和多元线性回归分析结果，新丰江水库的DO受WT、pH值、EC、I_Mn、TN和TP等多个因素的综合影响，其中WT是水库DO时空分布的最关键因素，其次为pH值。