珠江流域东江(东莞段)溶解氧时空变化及其影响因素*

徐 闯,刘广州,陈晓宏

(1：中山大学土木工程学院，广州 510275) (2：湖北省生态环境科学研究院水环境研究所，武汉 430072)

氧是生命体赖以生存的核心元素，水体中的溶解氧是水生生态系统中生命体和生物地球化学过程的基础[1-2]. 当水体的耗氧速率大于复氧速率时，溶解氧浓度开始降低，并对水生生态系统形成威胁. 当溶解氧低于3 mg/L时，大量浮游植物的卵和幼虫会在几天之内死亡；严重的低氧甚至会在水体中形成“死区”，没有生命体可以在该区域内存活[3-5]. 然而，随着社会经济的不断发展，全球河流均出现了不同程度的污染，溶解氧浓度持续下降[6-7]. 因此，研究溶解氧的时空分布及其影响因素对理解水体生物多样性和水体生物地球化学循环具有重要意义.

溶解氧的分布具有时空特异性. 空间上，我国水体溶解氧可以被分为3个区域：北部高饱和溶解氧区、南部低饱和溶解氧区和中部饱和溶解氧过渡区[8]. 时间上，溶解氧的分布具有自相似性，例如，Li等[9]研究了珠江流域沙河站溶解氧的变化规律，发现其存在25 h的迟滞. 溶解氧的时空分布由耗氧过程(底泥耗氧作用和硝化作用等)和复氧过程(浮游植物光合作用产氧等)共同主导，受到自然要素(水温、降雨等)和人类活动(污染物排放、水库和丁坝工程建设等)的综合影响[10-11]. 例如，在Mexico湾北部，过量营养盐输入湾内，加上水体层化的作用，导致湾内出现了超过15000 km2的低氧区[12]. 在珠江口及其邻近陆架区，流量和西南风主导的动力过程以及营养盐主导的生物化学过程共同影响溶解氧的时空分布[13]. 可以看出，不同区域溶解氧时空分布的调控因素不同.

在高度城市化城市河网区，社会经济发达，工业化程度高，大量废污水被排放进入河流，导致水体受到严重污染[14]. 以往研究表明，包括汾河、渭河、西江在内的许多河流均面临着溶解氧偏低问题，严重影响了河流生态系统的健康发展[9,15]. 基于此，许多学者采用不同的方法对河流水质变化及污染成因进行研究. 例如，王昱等[16]发现黑河上中游河流的主要污染物为氨氮、五日生化需氧量和高锰酸盐指数，且这些污染物会严重影响河流溶解氧浓度. 苏春雅[17]则发现水温、藻类密度和水体底质是扬州城区11条内河溶解氧的主要调控因素. 可以看出，河网区溶解氧受到水文要素和水质要素的共同调控，且作用机理复杂. 然而，先前研究主要针对河流综合水质进行分析，在此基础上探讨溶解氧的变化情况，单一针对溶解氧时空分布及其影响因素的研究不多. 同时，先前研究多采用皮尔逊相关系数来度量变量间的相关关系[18]. 然而，溶解氧变化呈现出典型的非线性特征，导致该系数无法准确地捕捉到变量间的响应关系[19]. 近年来，一种新的变量关系度量方法最大信息系数在许多领域得到了广泛应用，该方法可以有效地度量变量间的各种复杂关系[20]. Fan等[21]将皮尔逊相关系数和最大信息系数连用，成功地将流量和各影响因素之间的线性关系和非线性关系区分开. 然而，最大信息系数在水环境领域的应用还比较少.

东江是珠江流域的重要组成部分，近年来随着水环境整治力度的加大，河流水质有了一定的提升，但溶解氧偏低问题依然突出[22-23]. 在全球变暖和人类活动不断加剧的大背景下，东江溶解氧的时空分布是否会改变，主要影响因素是什么，这些问题亟待回答. 基于此，本文以珠江流域东江(东莞段)为例，采用实测数据对该区域溶解氧的时空分布进行研究，并结合最大信息系数、交叉小波转换和多元线性回归等多种方法量化分析各因素对溶解氧分布的贡献，以期为河流治理提供借鉴.

1 研究区域和数据

1.1 研究区概况

图1 研究区域及点位示意图Fig.1 Schematic diagram of study area and stations

东江发源于江西省寻乌县桠髻钵山，自东北向西南流经广东省，最后汇入珠江口及其邻近陆架区. 东江是珠江的3条主要支流之一，长度约为562 km，集水面积约为35340 km2(图1). 东江流域位于亚热带季风区，年降雨量约为1800 mm. 降雨量的年内分布不均匀，4-9月为丰水期，降雨量约占全年降雨量的80%；10月至次年3月为枯水期，降雨量仅为全年降雨量的20%. 社会经济方面，东江承担着为粤港澳大湾区主要城市约4000万人口提供生产生活用水的重要任务[24-25]. 东莞段位于珠江流域东江下游，上起东莞市桥头镇，下至东莞市沙田镇，河长约70 km. 根据东莞市生态环保局发布的水环境质量信息[26]可知，①当市内的饮用水源地(万江水厂、第六水厂和樟村水厂)水质不达标时，溶解氧是主要的超标项目；②市内入海口沙田泗盛断面的溶解氧几乎全年均处于不达标状态；③市内水质自动监测站(桥头水站和虎门水站)的主要污染物为溶解氧. 综合来看，东江东莞段面临显著的溶解氧不达标问题. 因此，选择该区域进行溶解氧时空变化及其影响因素研究是合适的.

1.2 研究数据及来源

本文选取东江东莞段桥头、樟村和沙田泗盛3个水质站的监测数据进行实例研究. 如图1所示，3个站点分别位于上游干流区域、中游河网区域以及下游入海口区域，具有较好的代表性. 站点实测数据包括溶解氧、水温、电导率、浊度、高锰酸盐指数和氨氮. 桥头站的数据范围为2013-2019年，樟村站的数据范围为2011-2019年，沙田泗盛站的数据范围为2018-2019年.

2 研究方法

采用折线图和箱型图分析溶解氧的空间变化、季节变化和年际变化；采用连续小波转换提取溶解氧的周期；采用交叉小波转换、最大信息系数和多元线性回归揭示溶解氧时空变化的驱动因素.

2.1 连续小波转换

连续小波转换可以将原始时间序列分解为时域和频域两个维度，进而获取更多数据细节，已经被广泛应用于时间序列的周期分析[27]. 其计算公式为:

(1)

式中，WT代表小波系数；a代表尺度参数，决定了信号的压缩程度；b代表位移参数，决定了信号的平移程度；x(t)代表原始信号序列.

根据前人研究[28]，本文采用Morlet小波作为母小波，基于小波系数的实部来分析溶解氧的周期性. 为了确定主周期，同时计算了不同周期对应的小波方差. 小波方差越大，代表其对应的周期越显著. 小波方差的计算公式为:

(2)

2.2 交叉小波转换

交叉小波转换能够同时分析时域和频域中两个序列之间多时间尺度的相关关系[29]，其计算公式为:

(3)

2.3 最大信息系数

最大信息系数由Reshef等于2011年提出[30]. 与皮尔逊相关系数等只能获取变量间的线性关系不同，最大信息系数可以有效地度量变量间的各种复杂关系，已经许多领域得到了广泛应用[20]. 最大信息系数在0和1之间变化，数值越大代表变量之间的关系越强. 最大信息系数超过0.30，则可以认为变量之间存在相关关系. 一般将最大信息系数和皮尔逊相关系数连用以量化变量间的线性关系和非线性关系，如公式(4)所示.

α=MIC-R2

(4)

式中，α代表非线性系数，MIC代表最大信息系数，R代表皮尔逊相关系数.α>0.20表示变量间的非线性关系较强；反之则表示变量间的线性关系较强.

2.4 多元线性回归

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

Y=a+b1X1+…+bnXn

(5)

式中，Y表示溶解氧，X表示水温等影响因素，a表示截距，b表示X和Y之间的回归系数.

图2 溶解氧月变化趋势Fig.2 Monthly trend of dissolved oxygen

3 结果与分析

3.1 溶解氧时空变化规律

3.1.1 溶解氧月变化规律 从图2可见，2011-2019年各年月均溶解氧浓度均呈现出先降低后升高的变化趋势. 一般而言，溶解氧浓度从1-2月开始降低，直到6-8月达到最小值，之后逐渐上升，年内整体变化呈现V字型. 以2018年为例，3个站点溶解氧最大值均出现在2月，分别为9.03、7.74和5.25 mg/L. 桥头站和沙田泗盛站溶解氧最小值出现在8月，分别为2.95和1.56 mg/L；樟村溶解氧最小值出现在6月，为3.17 mg/L.

3.1.2 溶解氧丰-枯水期变化规律 由图3a可知，桥头、樟村和沙田泗盛3个站点多年溶解氧平均值分别为6.66、5.68和3.11 mg/L，分别对应地表水Ⅱ类、Ⅲ类和Ⅳ类标准. 由此可知，从东江上游干流至下游入海口，溶解氧浓度呈现出逐渐递减的趋势. 尽管统计的年份并不完全一致，但从图2可以看出，3个站点的溶解氧在2018-2019年也符合空间递减规律. 由图3b～c可知，丰水期3个站点溶解氧均值分别为5.82、4.64和2.23 mg/L，对应地表水Ⅲ类、Ⅳ类和Ⅴ类标准；枯水期3个站点溶解氧均值分别为7.50、6.67和4.00 mg/L，对应地表水Ⅰ类、Ⅱ类和Ⅳ类水标准. 丰水期溶解氧浓度比枯水期分别低1.68、2.03和1.77 mg/L. 整体而言，东江东莞段溶解氧展现出了明显的丰-枯季变化.

图3 多年所有季节(a)、丰水期(b)和枯水期(c)溶解氧箱型图Fig.3 Dissolved oxygen box diagrams of all seasons (a), wet season (b) and dry season (c)

图4 溶解氧年际变化趋势Fig.4 Inter-annual variation of dissolved oxygen

3.1.3 溶解氧年际变化规律 图4为溶解氧浓度的年际变化趋势. 因沙田泗盛站仅有两年数据，因而未展示在图中. 由图4可知，桥头站溶解氧浓度在2013-2017年变化不大，约为7.00 mg/L；2018和2019年溶解氧浓度降低，分别为5.68和6.40 mg/L. 和桥头站相比，樟村站溶解氧浓度的年际变化较大. 2011年溶解氧浓度为5.50 mg/L，到2012年下降为4.55 mg/L，之后逐渐增加，至2017年达到最大值(6.50 mg/L)后又逐渐降低. 同时，丰水期和枯水期溶解氧的年际变化趋势和年均溶解氧的变化趋势类似. 另外，由图4也可以看出，桥头站溶解氧浓度高于樟村站，且丰水期溶解氧浓度低于枯水期，这和图2～3中得出的结论一致.

3.1.4 溶解氧周期性 为了进一步揭示溶解氧在不同时间尺度上的变化规律，采用连续小波转换提取了溶解氧的周期及其对应的小波方差，如图5和表1所示. 可以看出，桥头站和樟村站的溶解氧均在整个时域范围内呈现出8～16个月的震荡周期，对应的震荡中心均为11.64个月，说明这两个站点溶解氧存在显著的年周期. 此外，在部分时域内，桥头站溶解氧还存在6.24和18.91个月的周期，樟村站溶解氧还存在4.12和6.69个月的周期，尽管其对应的小波方差较小. 沙田泗盛站仅有两年数据，因此溶解氧未表现出年周期. 图5c显示2018和2019年溶解氧分别在0～4个月和5～7个月的周期上震荡. 表1显示该站点溶解氧的主周期为3.58 和5.43个月. 总的来说，东江东莞段溶解氧最显著的周期为年周期，此外还存在一些半年和季节等周期.

图5 桥头站(a)、樟村站(b)和沙田泗盛站(c)溶解氧周期性Fig.5 Periodicity of dissolved oxygen at Qiaotou (a), Zhangcun (b) and Shatiansisheng (c) stations

表1 3个站点溶解氧的主周期及 其对应的小波方差

3.2 溶解氧时空变化的影响因素

3.2.1 最大信息系数结果 2013-2019年桥头站溶解氧和水温之间的最大信息系数最高(0.72)，说明水温是溶解氧时空变化最重要的影响因素；溶解氧和高锰酸盐指数以及氨氮之间的最大系数值超过0.30，说明这两个因素也会在一定程度上影响溶解氧的变化. 不同之处在于，溶解氧和水温之间为线性负相关关系(R=-0.84)，和高锰酸盐指数以及氨氮之间为非线性关系(表2).

2011-2019年樟村站溶解氧和水温以及高锰酸盐指数之间的最大信息系数达到0.30以上，分别为0.61和0.30，说明水温仍是该站点溶解氧时空变化的主要影响因素. 其中，溶解氧和水温之间为线性负相关关系(R=-0.76)，和高锰酸盐指数之间为非线性关系(表2).

表2 溶解氧和各驱动因素之间的最大信息系数、相关性系数和非线性系数值

2018-2019年沙田泗盛站溶解氧和所有因素之间的最大信息系数值均达到0.30以上，分别为0.56(水温)、0.69(电导率)、0.35(浊度)、0.31(高锰酸盐指数)和0.48(氨氮). 其中，溶解氧和水温之间为线性负相关关系(R=-0.60)，和电导率之间为线性正相关关系(R=0.71)，和浊度、高锰酸盐指数以及氨氮之间为非线性关系(表2).

3.2.2 交叉小波转换结果 为了进一步揭示溶解氧和影响因素在不同周期下的响应关系，采用交叉小波转换对原始数据进行分析，结果如图6所示. 图中细实线代表影响锥边界，落在影响锥内部的为有效谱值区域；粗实线为显著性95%以上的置信区间. 颜色代表相关性：黄色越深，相关性越高. 箭头代表相位:箭头向左代表负相位关系，即变量之间变化趋势相反；箭头向右代表正相位关系，即变量之间变化趋势相同；箭头向上代表自变量超前因变量的变化；箭头向下代表自变量滞后因变量的变化.

图6 桥头站(a～e)、樟村站(f～j)和沙田泗盛站(k～o)溶解氧交叉小波转换结果Fig.6 Cross wavelet transform results of dissolved oxygen at Qiaotou (a-e), Zhangcun (f-j) and Shatiansisheng (k-o) stations

在桥头站，溶解氧和水温在全时间尺度的8～16个月周期上存在显著相关关系，和电导率、浊度、高锰酸盐指数以及氨氮仅在部分时域内存在相关性. 以电导率为例，在2014-2016年，其和溶解氧在3～15个月的周期上存在显著相关区域. 和桥头站类似，樟村站溶解氧和水温在全时间尺度的8～16个月周期上也存在显著的相关关系，而溶解氧和其他因素多在2014-2017年的8～16个月周期内存在显著相关关系，但其相位有所差异. 沙田泗盛站的数据较少，溶解氧主要和高锰酸盐指数在6～7个月的周期上存在显著相关关系.

3.2.3 多元线性回归结果 为了量化主要影响因素对溶解氧时空变化的贡献，采用多元线性回归来对溶解氧进行拟合. 根据最大信息系数和交叉小波转换的结果，本文首先采用单一水温指标拟合溶解氧，在此基础上叠加最大信息系数在0.30以上的其他因素拟合溶解氧，最后输入所有影响因素拟合溶解氧，结果如表3所示.

表3 桥头站、樟村站和沙田泗盛站溶解氧多元线性回归结果

在桥头站，水温可以解释溶解氧整体变化的70%，叠加其他因素后仅可将相关性系数提升2%，说明水温是该站点溶解氧时空变化最重要的影响因素，这和表2所得出的结论一致. 在樟村站，水温可以解释溶解氧时空变化的57%，叠加高锰酸盐指数后，相关性系数并未明显提升，叠加其他因素后，相关性系数提升至65%，这说明除水温外的其他因素贡献了溶解氧时空变化的8%. 在沙田泗盛站，水温仅可以解释溶解氧时空变化的33%，叠加电导率、浊度、高锰酸盐指数和氨氮后，可以解释整体变化的64%. 这说明其他因素的综合作用可以与水温的作用相当.

4 讨论

本文基于桥头、樟村和沙田泗盛3个水质站点的监测数据，在月、季节和年际3个时间尺度下研究了溶解氧的时空分布规律；同时采用最大信息系数、交叉小波转换和多元线性回归等多种研究方法对溶解氧的影响因素进行解析. 可以看出，东江东莞段溶解氧展现出了显著的时空变化，且该变化受到多重因素的综合调控.

溶解氧浓度在1-12月大体呈现出先递减后递增的趋势，具有明显的丰、枯季变化规律. 实际上，这一规律在珠江口、香港大鹏湾等区域也有发现[13, 31]. 如Zhang等[31]研究了珠江口溶解氧的年内变化，发现磨刀门附近枯水期溶解氧浓度可以达到8 mg/L，然而丰水期仅约3 mg/L. 在本文中，桥头、樟村和沙田泗盛3个站点丰水期溶解氧浓度比枯水期分别降低1.68、2.03和1.77 mg/L. 这一差异的原因可以被归结为以下两点：丰水期水温较高，使得氧气在水体中的溶解度变低；根据前人研究，珠江流域养殖业、种植业和农村污水等带来的氮、磷面源污染较多[32]，丰水期降雨量远大于枯水期，降雨冲刷作用使得水体中耗氧物质增加，溶解氧浓度降低.

年际上，桥头和樟村2个站点溶解氧均以11.64个月为周期进行波动. 相比而言，樟村站溶解氧的年际差异较大. 其中，2012和2013年溶解氧浓度较低(4.55和4.99 mg/L)，2017年溶解氧浓度最大(6.50 mg/L). 这3年樟村站的水温、电导率、高锰酸盐指数和氨氮浓度较为接近，但浊度差别较大. 2012和2013年的浊度分别为26.75和31.37 NTU，远高于2017年的7.18 NTU. 溶解氧和浊度之间存在负相关关系(表2)，水体浊度增加不利于水生生物的光合作用，限制了溶解氧的源，使其浓度降低. 这一结果与Hu等[33]在珠江口的最新研究一致.

就空间变化而言，在桥头和樟村站，水温是溶解氧时空变化的首要影响因子，两者之间存在显著的线性负相关关系. 此外，表2显示这两个站点的高锰酸盐指数和氨氮也会调控溶解氧的变化. 高锰酸盐指数和氨氮实际上反映了水体中各有机物和无机物的耗氧过程，其中氨氮引发的硝化反应会消耗溶解氧进而使其转化为亚硝态氮和硝态氮，而高锰酸盐指数则反映了各有机物和无机物的氧化过程. 在下游沙田泗盛站，溶解氧仅达到地表IV类水标准，比上游降低了3.55 mg/L. 该站点的水温和上游2个站点类似. 但是高锰酸盐指数为3.67 mg/L，氨氮浓度为0.96 mg/L，均高于上游站点，说明该站点水体中各污染物的耗氧过程更加显著. 同时，沙田泗盛站的浊度为144.75 NTU，比上游站点高了约100 NTU，说明浊度对水生生物光合作用的限制也影响了该站点溶解氧浓度. 值得注意的是，沙田泗盛站的溶解氧和电导率之间呈现出了显著的正相关关系，这与上游站点不同. 沙田泗盛站距离入海口最近，受到的潮汐作用最为显著，潮汐携带的大量氯离子导致该站点电导率增加. 实际上，潮汐可以促进水体混合，进而缓解水体低氧[34]. 但是沙田泗盛站位于多支流交汇处，接纳了来自上游各支流的大量污染物，水体耗氧过程占据了主导地位，导致溶解氧浓度下降.

综合来看，东江东莞段溶解氧偏低问题突出，尤其是沙田泗盛站. 针对这一问题，提出以下几点建议:①源头管控. 当前对于养殖业，农村生活污水等面源污染的摸排是不到位的，应该加强排查，摸清本底；另外，考虑到水动力、浮游植物和底泥等因素均会对溶解氧浓度产生较大影响，应该完善监测网络，在常规监测中同步测量流量、潮位、叶绿素a和底质相关指标；同时，应该在沙田泗盛上游各支流布设监测站，以明确不同支流的贡献. 在此基础上，可以采取提高排污标准，科学施肥等方式在源头上减少污染物的排放. ②过程管控. 要加快管网建设，推进雨污分流；同时，要加大巡查力度，快速发现违法排污问题并及时进行管控；③推进智能化平台建设，整合预报、预警和溯源等功能，以便未雨绸缪.

河流溶解氧的时空变化受到自然要素和人类活动的综合影响. 自然要素中，本文主要考虑了水温的变化，研究结果也显示水温是溶解氧时空变化的主要影响因素. 就人类活动而言，工厂、船舶和农田等点、面源污染会带来大量的耗氧物质，影响水体中溶解氧浓度，这一点在本文中已有体现. 除此之外，人类活动主导的水力调控也会影响溶解氧浓度，例如水库泄洪和水坝构建等水力调控措施可以影响河流的水动力条件，进一步影响溶解氧浓度[35]. 在后续研究中应进一步考虑水力调控的作用. 另一方面，溶解氧浓度的变化也会对水生生态系统造成严重影响，例如鱼类产卵场对溶解氧的变化非常敏感[36]. 因此，在后续研究中要进一步关注溶解氧引发的生态效应. 总的来说，本文基于实测数据分析了珠江流域东江东莞段溶解氧的时空变化规律及其影响因素，提出了对应的管控措施，并探讨了后续的研究方向，相关结果可以为本流域及其他流域的水环境综合整治和水生态管控等提供参考.

5 结论

1)东江东莞段溶解氧存在显著的时空变化. 空间上，从干流桥头站至入海口沙田泗盛站多年平均溶解氧浓度下降了3.55 mg/L. 季节上，桥头、樟村和沙田泗盛3个站点丰水期溶解氧浓度分别比枯水期低1.68、2.03和1.77 mg/L. 年际上，不同年份的溶解氧呈现出波动特征，主要波动周期为11.64个月.

2)最大信息系数、交叉小波转换和多元线性回归能够较好地揭示溶解氧和各影响因素之间的响应关系. 在桥头站和樟村站，水温是溶解氧时空变化的首要影响因子，可以解释其整体变化的57%以上. 此外，高锰酸盐指数主导的氧化作用和氨氮主导的硝化作用是溶解氧重要的汇，所有因素叠加后可以解释溶解氧整体变化的65%以上.

3)沙田泗盛站受上游多条支流汇入的影响，整体水质较差. 电导率、浊度、高锰酸盐指数和氨氮浓度均高于上游站点. 增强的水体耗氧作用及减弱的光合作用导致该站点溶解氧浓度低于上游站点. 总的来说，所有因素叠加后可以解释该站点溶解氧整体变化的64%.