19882018年滇池氮磷比的时空演变特征与原因解析*

2021-01-12程国微严长安

湖泊科学 2021年1期

高伟，程国微，严长安，陈岩

(1：云南大学生态与环境学院，云南省高原山地生态与退化环境修复重点实验室，昆明 650091) (2：昆明市生态环境科学研究院，昆明市环境系统模拟与风险评估重点实验室，昆明 650032) (3：生态环境部环境规划院，国家环境保护环境规划与政策模拟重点实验室，北京 100012)

氮和磷是影响淡水生态系统初级生产力的两种主要营养元素，氮、磷富集导致的富营养化是全球面临的普遍水环境问题[1]. 氮磷比是水体氮、磷元素的质量或原子数量比值，是水体氮、磷相对水平的综合体现. 作为影响藻类生长的重要化学计量指标之一[2]，氮磷比变化对水域生态系统结构和功能具有重要影响，如低氮磷比与水体固氮蓝藻的出现具有关联性[3]，氮磷比也被用于判识水体氮磷限制类型[4-5]. 现有研究表明氮磷比对藻类的影响较为复杂，藻类生长受氮磷比和氮、磷在水体中绝对浓度的共同影响[6]，且在不同的水体中响应机制存在差异. 此外，水体氮磷比的变化不仅影响本地水体的生态系统健康，还可能对下游水体造成不利的生态影响[7]. 由于氮磷比的变化不仅取决于氮、磷的绝对浓度变化，还受到氮和磷相对浓度变化的直接影响，因此氮磷比的主要影响因素既是氮和磷浓度影响因素的集合，也与影响它们相对浓度的因素有关. 现有研究主要关注氮磷浓度的变化及其原因，但对它们相对浓度变化即氮磷比的驱动研究较少. 目前对湖泊氮磷浓度影响的定量研究主要包括相关性分析法、统计学建模和机理过程模型3种. 相关性分析法是基于相关系数(线性相关系数、非线性相关系数、灰色关联度等)判断因素之间的相关强度，如通过计算水质浓度与驱动因素之间的灰色关联度识别影响湖泊水质指标变化的因素[8]，该方法计算较为简单，但只能识别单因素的独立作用. 统计模型法是建立湖泊氮、磷浓度与多个驱动因子之间的回归关系，通过统计显著性筛选驱动因子，通过回归系数判断驱动强度，该方法在湖泊水质变化研究中得到广泛应用[9-10]，常用的统计方法有多元线性回归[9]、随机森林[11]和因素分解法. 机理模型法主要是采用湖泊模型建立不同输入条件与湖体水质的定量关系，通过运行情景模拟出湖体水质变化的主要因子，常用的有EFDC水动力-水质模型[12]. 然而，机理模型建模的难度较高，对数据量和建模者使用经验有较高要求，在实际应用中难度较高. 湖泊氮磷浓度的影响因素较为复杂，主要包括流域氮磷输入负荷[12]、气象条件[10]和氮磷削减工程措施[13]等. 统计模型能够有效地处理多因素综合影响效应问题，模型相对简单，本研究使用统计模型分析滇池氮磷浓度和氮磷比时空变化的影响因素及其效应.

滇池的富营养化持续历史较长，氮磷超标是滇池水环境治理长期面临的难题. 随着城市污水处理、环湖截污和生态补水等治理工程的实施，滇池的氮磷浓度发生了显著变化[11]，湖体氮磷比随之改变，并对湖体优势藻类产生了影响[14]. 目前针对滇池的氮磷比演变特征及其驱动力的研究相对较少，主要集中于氮磷比对藻类影响上[14]. 对滇池氮磷浓度变化与来源的研究相对较多[11-12, 15]，但主要集中于影响相对重要性的识别和污染来源角度. 滇池氮磷比长期演变趋势及其多种因素的综合驱动影响尚缺乏相关研究. 本研究基于滇池1988-2018年的多年水质观测数据，评估了滇池氮磷比在时间和空间上的变化趋势；基于氮磷比的影响因素，定量解析了不同因素对滇池氮磷比变化的贡献，从而识别滇池氮磷比的主控因子. 氮磷比变化与驱动效应研究对科学诊断湖泊生态系统状态变化方向并制定调控方法具有重要意义，本研究可为湖泊氮磷比变化趋势与驱动效应识别提供案例参考和方案借鉴.

图1 滇池流域及其气象、水文与水质监测站点分布 Fig.1 Location of metrological, hydrological and water quality gauging stations of Dianchi Basin

1 研究区概况

滇池是我国第六大淡水湖泊，是云贵高原面积最大的湖泊，湖泊总面积309.5 km2(1887.4 m)，湖容15.6亿m3，平均水深5.3 m(图1). 1996年滇池被海埂大坝分割为北部的草海和南部的外海两个子湖，外海是滇池的主体，面积为289.2 km2. 草海由人工出口西园隧道出流，外海由海口河出流，两者已没有水量联系. 滇池水质自1970s末开始恶化，成为我国水体富营养化严重的“三湖”之一. 为改善滇池的水环境，滇池流域自“八五”以来实施了一系列污染防治工程[16]，水质得到有效改善[11]. 2013年12月底牛栏江-滇池补水工程实施，牛栏江调水成为滇池的重要水源，对滇池草海和外海开展了生态补水. 相对于巢湖和太湖，滇池在“三湖”中的富营养化程度最高[15]. 根据云南省2018年环境状况公报[17]，滇池草海的总氮为劣V类，总磷为IV类，外海总氮和总磷均为IV类. 虽然近年来水质有所好转，特别是2009年以来滇池草海和外海的总氮与总磷浓度有显著下降[15]，但仍然未达到水环境功能区(草海IV类，外海III类)的要求，滇池仍然处于轻度富营养化状态.

2 材料与方法

2.1 研究方法

2.1.1 氮磷比时间变化趋势分析指标时间变化趋势判别是水环境质量评估的重要内容，目前主要采用主观判断和统计分析两种方法. 主观判断是根据指标数值随时间增长的变化方向判断数据的趋势，在湖泊污染指标变化趋势分析中获得广泛应用[18]，但存在主观性较强且在数据波动较大时不确定性极高的问题. 统计分析方法是时间序列分析的另一主流方法，常用的有参数检验法和非参数检验法. 参数检验法最常用的是线性趋势回归检验[19]；非参数检验法类型较多，常用的有Mann-Kendall检验法和Daniel检验法[8]. 本研究选择线性趋势检验和常用的Mann-Kendall非参数突变检验法对滇池氮磷比的时间变化趋势与突变进行判别，其中Mann-Kendall突变检验具体方法参见文献[20]，该方法不要求时间序列的数据满足正态分布.

2.1.2 氮磷比驱动效应模型驱动效应模型是用于解析影响因子对被影响变量变化贡献的数学模型. 本研究使用多元线性回归法，该方法是通过建立自变量和因变量的定量关系，基于回归系数和变量数值测算自变量对因变量的贡献. 该方法在污染源解析中应用较多[21]，近年来在水文变化驱动力分析中也有应用[22]. 由于多元线性回归模型具有结构清晰、操作简便和应用范围广的特点，本研究选择多元线性回归模型对氮磷比驱动效应进行解析. 多元线性回归模型对驱动因子的效应解析分为两步，第一步是构建多元回归方程，采用逐步回归法识别出具有统计显著性(P<0.05)的驱动因子及其回归系数；第二步，基于回归系数和变量均值测算各因子对因变量的贡献. 具体公式为：

(1)

(2)

(3)

为进一步分析湖泊氮和磷浓度变化的影响因素，建立了氮、磷浓度与降水、气温、流域氮磷输入等驱动因子的多元线性回归模型，方程的形式与公式1相同. 为了消除数据偏度分布的影响，对自变量和因变量进行了取自然对数处理[10]，对取值有0的数据进行Z值(归一化为均值为0、标准差为1的Z值)转换.

2.1.3 氮磷比驱动力计算模型在本研究中，根据文献调研和数据可得性，主要评估流域氮磷输入、气象条件和水文条件3类驱动因子的驱动效应. 其中流域氮和磷输入分别采用NANI[23]和NAPI[24]模型计算，两者均是基于物料平衡原理计算流域的人为氮、磷输入量，NANI估算流域化肥施用输入、大气沉降输入、食品与饲料净输入和固氮农作物种植输入等来源的氮输入量，NAPI估算流域化肥施用、食品与饲料净输入和非食品输入3项磷输入量，NANI和NAPI的具体模型与参数参见文献[9, 25]. 气象条件是基于逐日降水和逐日平均气温计算了滇池流域的平均降水、汛期降水和平均气温. 水文指标主要使用湖泊平均水位作为水文驱动因子. 污染物削减和调水数据分别采用污水处理厂削减量和牛栏江滇池补水工程调水量，表征污染控制措施.

2.2 数据来源

根据研究内容，使用的数据分为4类：滇池氮磷监测数据、气象水文数据、人为氮磷输入数据和污染控制措施(污水处理厂削减和外流域调水). 氮、磷数据包括滇池流域草海和外海10个点位(图1)1988-2018年逐月总氮和总磷浓度数据(表1)，数据来源于昆明市生态环境监测中心. 气象数据主要收集了降水、湖泊水位和气温，其中降水数据收集了1988-2018年流域27个降水站(图1)的逐日降水量，水文数据收集了滇池大观楼、海埂、中滩站3个水位站的同期逐月平均水位(图1)，降水和水文数据来源于云南省水文水资源局；同期昆明、晋宁、太华山和呈贡4个国家气象站(图1)的逐日平均气温数据来源于资源环境数据云平台(http://www.resdc.cn/Default.aspx). 滇池流域1988-2015年的NANI和NAPI数据来源于文献[13]，2016-2018年基于经济社会统计数据和文献相同的方法计算，统计数据来源于昆明市统计年鉴. 污染物削减和调水数据来源于昆明市生态环境保护局和昆明市水务局.

表1 1988-2018年滇池氮、磷浓度基本统计特征

3 结果与讨论

3.1 滇池氮磷比的时间演变趋势

基于滇池草海和外海的多年氮磷监测数据，得到两个子湖泊氮磷比的逐年变化. 1988-1997年草海的氮磷比平均为21.3±11.2，到2009-2018年增长到51.3±13.7；外海从1988-1997年的20.9±5.10上升到近10年的32.8±4.9. 线性趋势检验结果表明，1988-2018年滇池草海和外海的氮磷比均呈现显著线性上升趋势，其中草海年均上升1.3(P<0.01)，外海年均上升0.7(P<0.01)，草海氮磷比上升速度是外海的2倍(图2). 草海和外海氮磷比的Mann-Kendall趋势检验的统计值分别为4.62和3.98，均呈显著上升趋势(P<0.01)，与线性趋势检验的结果一致. 然而，值得注意的是，虽然草海和外海的氮磷比均呈现上升趋势，但草海、外海的氮磷比线性相关系数仅为0.32，低于0.05的显著性水平，说明两者的环境驱动因素可能存在差异. 从Mann-Kendall的突变检验结果看，研究期间滇池草海和外海的氮磷比均存在显著突变. 其中，草海在2008年发生了上升突变，外海在2004年也发生了上升突变，两湖的突变方向一致，但时间上草海滞后于外海(图3).

图2 1988-2018年滇池氮磷比变化Fig.2 Change of nitrogen to phosphorus ratio of Lake Dianchi from 1988 to 2018

3.2 滇池氮、磷浓度的变化特征

根据氮磷比的突变时间，对突变前后的湖体氮、磷浓度变化进行了分析(图4). 对于草海，突变前后的氮、磷浓度存在相反的变化趋势，1988-2008年，草海的TN和TP浓度均呈现显著的上升趋势(P<0.01)，2008-2018年则转变为显著下降(P<0.01). 从两个阶段的斜率看，突变前的TN和TP年均浓度分别上升0.49 和0.04 mg/L，相对于突变前期(1988-2008年)均值分别提高0.05和0.04倍，TN浓度的相对增长速度较TP快，因此导致前期氮磷比呈现上升趋势；突变后期(2008-2018年)草海TN和TP年均浓度分别下降

图3 1988-2018年滇池(a：草海；b：外海)氮磷比Mann-Kendall突变检测(UFk和UBk分别为Mann-Kendall正向和逆向序列的统计检验值)Fig.3 Abrupt change of nitrogen to phosphorus ratio of Lake Dianchi (a: Lake Caohai; b: Lake Waihai) from 1988 to 2018 based on Mann-Kendall test (UFk and UBk are statistics of Mann-Kendall)

图4 1988-2018年滇池氮、磷浓度变化(a：草海；b：外海)Fig.4 Temporal change of nitrogen and phosphorus concentrations of Lake Dianchi from 1988 to 2018 (a: Lake Caohai; b: Lake Waihai)

1.31 和0.11 mg/L，相对于突变后期均值分别下降0.18和0.26倍，TN浓度的下降速度慢于TP，从而导致氮磷比上升. 对于外海，突变前期(1988-2004年)仅TN浓度有显著上升趋势(P<0.01)，因此突变前期氮磷比上升是TN浓度升高导致的结果；突变后期(2004-2018年)的TN和TP浓度均呈现显著下降趋势(P<0.01)，由于TN浓度较TP的相对下降速度慢(年均下降量相对于突变后期均值分别为0.03和0.04)，导致氮磷比呈现上升趋势. 综上，滇池草海和外海的氮磷比在突变前期与后期均呈现上升趋势，但原因有所差异. 滇池草海在突变前是由于TN浓度上升相对速度较TP快，突变后是由于TN浓度下降速度较TP慢；外海在突变期是TN浓度上升单一因素导致的，突变后是TN浓度下降速度慢于TP导致的.

综合现有文献[8,10-11,13,25]，可将湖体氮、磷水平的影响因素归纳为流域气象、水文、污染负荷和治理水平等. 气温可通过影响氮磷相关化学过程和微生物活动改变湖体的氮、磷水平[26]，进而对氮磷比可能造成影响；降水是湖泊水量的直接来源之一，同时也是陆域非点源入河的关键影响因素，因此可通过改变湖体水文和污染入湖量影响氮、磷浓度水平[10]；研究表明，风速是湖泊浮游植物生物量的重要影响因素，成为部分湖泊蓝藻水华暴发的重要驱动力[26]；水位是湖泊水量的直接表征指标，水位变化可通过影响水环境容量改变湖体氮、磷浓度[13]；人为氮、磷输入是流域最大的污染来源，可直接改变湖体的氮、磷水平[11, 13]；调水量对湖体氮、磷的影响主要是通过水文和氮、磷浓度变化两个途径[13]；此外，流域的污染减排措施可降低氮、磷的入湖水平，对湖体的氮、磷浓度及其相对水平产生影响[13].

以氮、磷浓度为因变量，以流域年均气温、年降水、汛期降水、年均风速、年均水位、NANI、NAPI、年调水量、流域TN削减量、流域TP削减量10个变量为自变量，对1988-2018年滇池氮、磷浓度进行多元线性回归，结果表明，草海的TN浓度主要受人类活动净氮输入量、流域TN削减量和滇池水位控制，其中人类活动净氮输入与草海TN浓度呈正相关，而流域TN削减量和滇池水位与草海TN浓度呈负相关. 草海的TP浓度受到人类活动净磷输入量、流域TP削减量、牛栏江调水和风速的综合影响，其中调水量和流域TP削减量表现出显著的负相关，说明近年来的调水和污染治理工程对滇池草海TP浓度的降低有显著作用. 对于外海，影响TN浓度的主要变量是流域人类活动净氮输入和外流域调水. 从可决系数看，外海TN浓度对这些变量的线性响应关系不如草海，可决系数低于草海，但仍然能够达到0.613；外海TP浓度回归方程的可决系数较低，仅0.356，只有调水量一个指标进入方程(表2)，说明牛栏江调水对滇池外海的TP浓度产生了较显著的抑制作用，但是滇池外海TP浓度变化过程较为复杂，可能存在复杂的非线性过程和底泥的源汇作用[27]. 总体来看，滇池氮磷浓度受陆域氮磷输入负荷、跨流域调水、流域氮磷削减和气象水文条件影响，其中流域氮输入负荷对草海和外海的总氮有显著影响，磷输入负荷对草海总磷有影响，这与基于随机森林统计模型得

表2 滇池氮、磷浓度的多元线性回归模型

出的结论较为一致[11]，流域氮、磷削减对草海氮、磷浓度影响较大，跨流域调水量对外海的氮、磷浓度和草海的总磷浓度影响较大.

3.3 滇池氮磷比的外在驱动因子及其变化特征

湖泊氮磷比是湖体总氮和总磷浓度的相对数量关系，因此直接取决于湖体的氮、磷水平. 基于氮、磷驱动的主要影响因素，本研究使用11项指标，分别表征流域气象条件、水文条件、污染负荷和污染治理水平. 滇池流域的氮、磷驱动因子变化如图5所示.

图5 1988-2018年滇池氮磷比驱动因子变化Fig.5 Driving forces of nitrogen to phosphorus ratio in Lake Dianchi from 1988 to 2018

3.4 滇池氮磷比变化的驱动因子及其效应解析

以湖体氮磷比为因变量，以驱动因子为自变量，分别构建草海和外海的多元线性回归模型. 基于逐步回归法，得到草海和外海的氮磷比回归方程(表3). 影响草海氮磷比的主要驱动因素是流域总磷削减量、人为输入氮磷比和平均气温，3个因子可解释草海氮磷比变化的80.4%，且均为正向驱动效应. 相对于草海，影响外海氮磷比具有显著性的驱动因子只有NANI和平均气温，回归模型的解释力度为0.676，低于草海，表明外海氮磷比的影响因素和氮、磷变化过程更为复杂. 平均气温对外海氮磷比的影响是负向驱动，即气温越高，氮磷比越低，这与草海的效应相反，说明温度在外海和草海的氮、磷循环中的作用可能存在差异. 已有研究结果表明[28]，气象因子(包括温度)对外海TN浓度的影响高于TP，且主要表现为负向效应，这与本研究较为一致，气温越高，氮浓度越低，氮磷比呈现下降趋势. 根据驱动因子贡献的解析公式，得到滇池草海和外海各驱动因子对氮磷比的影响(图6). 对草海，年平均气温是氮磷比变化的最主要影响因素，其次是流域的人为输入氮磷比，两者对草海氮磷比的贡献为53%，流域的磷削减量对氮磷比的贡献较小，仅为2%，此外还有45%的氮、磷是由其他未识别因素贡献的，可能为草海存量；温度同样是外海氮磷比的最主要影响因素，占46%，其次是流域人为氮输入，占14%，剩余40%的贡献来自于其他未识别因子. 由于底泥是影响滇池氮、磷平衡的重要因素，未识别因子可能主要来自于底泥. 在已有的滇池污染源解析研究中，根据未识别源的贡献特征将其定义为底泥内源[28]. 湖泊的氮、磷浓度受控因素较多，生物因素也是影响氮、磷的重要因素. 本研究在文献综述基础上选择了这些非生物因素，未考虑生物因素的可能影响，主要原因是缺乏长时间序列生物观测资料. 此外，生物因素的影响也会与非生物因素耦合且有时互为因果，如氮、磷浓度变化可改变浮游植物数量，浮游植物数量反过来又改变氮、磷浓度. 现有研究一般将氮、磷浓度变化作为浮游动植物变化的驱动，而不是结果. 综上，本研究主要考虑非生物因素的影响.