基于时间序列的鄱阳湖Chl-a预测模型优化构建

钱春龙,曾一川,袁伟皓,吴 怡

(1.江苏中瑞咨询有限公司,南京 210036; 2.河海大学 环境学院,南京 210098)

0 引 言

水体富营养化问题是全球水环境由于人口数量的增加及工业化进程加快所必须面临的长期治理挑战[1]。全球范围内的很多重要湿地、湖泊、河流,如拉姆萨尔湿地[2]、肯尼亚维多利亚湖[3]、加拿大安大略湖[4]、巴西米纳斯·格莱斯多斯河[5]等,水体均面临承载不同程度的营养负荷压力。水质恶化会为饮用水安全、水生态系统健康带来系列严重危害。

随着江西省及环鄱阳湖经济带的快速发展,鄱阳湖湖区污染负荷日益加重。2018年江西省生态环境厅组织江西省环境保护科学研究院对鄱阳湖水体总磷(TP)污染的成因和来源进行了分析,形成了初步结论:近5 a来,鄱阳湖水体总体达到《地表水环境质量标准》(GB 3838—2002)中Ⅳ类标准要求,Ⅰ—Ⅲ类水质优良点位比例总体呈下降趋势。由于湖体北部连通长江,湖流形态与水沙交换与其他湖泊有着较大差异,体现出显著季节性[6],吸引多种迁徙候鸟[7]及包括江豚[8]在内多种珍稀物种中转或栖息。因此全面考虑多项水质指标,综合评估水质状况、提高藻类增殖等可能保护区水生态安全威胁因子模拟预测水平,对鄱阳湖地区物种资源、水资源、生态资源良性化发展利用具有现实意义。

近年来,关于鄱阳湖水质评价方法多样,包括聚类分析[9]、水质指数法(Water Quality Index,WQI)[10]、综合营养状态指数法[11]等,分别从各自角度对鄱阳湖水质保护提供重要参考。但这些方法存在评价体系复杂,当存在劣Ⅴ类监测值时无法完全评价等情况。考虑到总磷(TP)、氨氮(NH3-N)是鄱阳湖主要水污染指标[12],且总氮(TN)多处于《地表水环境质量标准》(GB 3838—2002)中Ⅴ类标准或劣Ⅴ类;高锰酸盐指数(CODMn)作为表征水体中还原性物质总量的一项指标,对于衡量水体含氧量因生物降解消耗而变化具有重要意义。

本文选择2012—2020年间对鄱阳湖涵盖三大重点保护区的4处代表性采样点的NH3-N、TN、TP、CODMn等4项水质因子进行月度监测,并采用综合水质标识指数法(Water Quality Identification Index,WQII)进行综合分析。依托单因子指数法,该法可以量化水质参数监测值在地表水5级水质标限(包括劣V类标准)之中的具体位置[13],并综合考虑各参数的影响。依托上述8 a时间序列数据同步构建以叶绿素a(Chl-a)为表征藻类的因变量的多元线性逐步回归方程(Multiple Linear Stepwise Regression,MLSR)、季节性自回归求和移动平均(Seasonal Autoregressive Integrated Moving Average,SARIMA)模型。鉴于水温是藻类增殖的主要限制因子,选择2020年6—8月份水温增幅明显时期的Chl-a浓度预测值与实测值偏差比较,检验2种模型对于各监测点位的适用程度与抗干扰性。该套组合方法对更全面掌握水质综合状态,保障环鄱阳湖地区饮用水、水生态安全,增强蓝藻水华爆发风险预测能力具有现实意义,同时以期为各类水体关键生态因子评价预测提供参考。

1 材料与方法

1.1 研究区域与监测点位

鄱阳湖流域气候类型属于亚热带季风气候,位于江西省的北部、长江中下游的南岸(115°4′13″E—116°24′6″E,26°44′40″N—29°44′48″N)[14],是中国最大淡水湖,也是世界第23大淡水湖,同时还是世界上重要的生态湿地,在维护生物多样性和保障防洪安全方面发挥着重要的作用[15]。鉴于鄱阳湖面积广、保护区众多及监测工作开展情况,选择包括长江江豚保护区、国家种质资源保护区、鄱阳湖国家级自然保护区在内的三大保护区作为分析对象;代表性监测点选择康山、蚌湖、都昌、蛇山4个点位,获取自2012年起的月度监测数据对历史水质变化进行综合研判分析。监测过程执行《水环境监测规范》(SL 219—2013)规定,选择《地表水环境质量标准》(GB 3838—2002)的Ⅲ类标准作为水质要求标准。研究区域保护区分布及点位从属关系如图1所示。

1.2 综合水质标识指数法(WQII)

就TN指标而言,康山、蚌湖、都昌、蛇山4个点位均有属于劣Ⅴ类监测值,分别占检测时段总体样本的8.33%、18.52%、17.59%、21.30%。WQII评价体系可以将这种超过Ⅴ类标准的情况囊括,同时对所有监测值进行归一化、定量化判断。

(1)单因子水质指数。首先将所选择的4个水质参数指标依照GB 3838—2002对于各参数标准限值计算出各自单因子水质指数,将不同浓度限值因子标准归一化,便于进一步处理。单因子水指数质公式为

式中:Qi为第i项水质因子的水质指数;Ci为第i项待评水质参数的水质类别,取值为0,1,2,…,Ci≥5时表示为劣Ⅴ类;Pi为待评价数据在Ci类水与Ci+1类水之间与中所处分位,即待评价数据与所在水质类别的标准下限值之比,保留一位小数(四舍五入);ρi表示水质因子实际浓度;ρCi和ρCi+1分别表示该水质因子Ci类水与Ci+1类水在GB3838—2002中的实际浓度。

(2)综合水质标识指数。该指数为单因子水质指数的综合体现,其表达式为

WQII=W1.W2W3W4。

(2)

式中:W1.W2为综合水质指数,由各项CiPi取均值确定(W2为四舍五入保留一位的小数);W3、W4为标识码,W3表示低于所评价的水功能区水质要求标准(此次评价选择Ⅲ类)的水质因子个数;W4为所评价的水功能区水质要求标准(此次评价选择Ⅲ类)与综合水质指数的差距[16],如果综合水质指数更好则为0。WQII法评价流程实例示意如图2所示。

图2 WQII法的评价流程Fig.2 Assessment process of the WQII method

1.3 预测模型体系

1.3.1 MLSR方程

(3)

式中b0,b1,…,bn表示线性回归系数。

依照最小二乘法原理,计算因变量的观测值与估计值之差的平方和,选出该值最小的一组作为方程的线性回归系数。进一步通过逐个引入并剔除偏回归平方和不显著的备选自变量[18],排除严重多重共线性,以得到最优MLSR方程。该方程的建立通过Minitab 18实现。

1.3.2 SARIMA模型

1.3.2.1 模型公式

时间序列方法是通过综合考虑研究对象的现状及过去以及各因素间相关关系所建立的预测模型,在评估研究对象时需要考察其历史演化规律[19]。综合水质标识指数法则以一种直观且更深入的方式呈现水质指标在一段时间的观测值(即时间序列)上的变化情况,并综合评估水质情况。依托该时间序列,Chl-a浓度的季节性自回归求和移动平均模型(SARIMA)通过时间序列中所包含规律的建立模型[20]并通过运行R语言可视化平台Rstudio实现。如果时间序列在经历s个时间间隔之后出现相似的变化趋势,则该序列周期为s,阶数用(p,d,q)×(P,D,Q)s表示[21]。其中(p,d,q)和(P,D,Q)s分别是模型的非季节性和季节性部分[22],公式为

Φp(L)AP(Ls)ΔdΔDsyt=Θq(L)BQ(Ls)εt。(4)

式中:Φp(L)和Θq(L)分别是平稳的自回归算子和可逆的移动平均算子[23];AP(Ls)和BQ(Ls)分别为非季节性的自回归算子和移动平均算子;Δd、ΔDs分别为非季节性差分算子和季节性差分算子;yt为时间序列的观测值(t为时间序列长度);L为滞后算子;p、d、q分别是非季节性的自回归阶、差分次数、移动平均项的阶;P、D、Q分别是季节性周期的自回归阶、差分次数、移动平均项的阶;εt为随机模型具有独立同分布和共同方差的随机误差[24]。

1.3.2.2 最优SARIMA模型选取

作为SARIMA模型中的一个重要假设,监测数据时间序列必须是平稳的,因此运用单位根检验(Augmented Dickey-Fuller Test,ADF)判定序列数据的平稳性。具体采用R语言可视化平台Rstudio内置tseries程序包中的adf.test函数对序列进行检验,显著性水平为0.01,对立假设为“序列平稳”。当差分n次时,若出现p<0.01,则否定原假设,接受对立假设,即此序列d=n。本次模型预测的4个监测点除都昌差分次数为1以外,其余序列均差分次数均为0。

基于赤池信息量准则(Akaike Information Criterion,AIC)[25],forecast程序包auto.arima函数可以实现AIC值最小的SARIMA作为最优模型的自动选取。此时还需要将d值替换成由adf.test函数确定好的差分次数以消除模型的趋势性(auto.arima函数输出最优模型后带有“with drift”字样),成为各监测点最优SARIMA模型的最终形态。模型采用的训练样本为2012年1月—2020年5月监测数据,预测2020年6—8月Chl-a浓度,并与实测数据比较。

2 结果与分析

2.1 基于综合水质标识指数法的水质评价

2.1.1 鄱阳湖重点保护区监测点单因子水质指数分析

2.1.1.1 康山监测点单因子水质指数分析

如图3所示,对于南矶湿地国家级自然保护区代表监测点康山,CODMn、NH3-N、TP、TN这4种水质参数单因子水质指数均值从最优开始排布为:CODMn(Qi=2.2)>NH3-N (Qi=2.3)>TP(Qi=2.5)> TN(Qi=4.7)。就时间序列方差来看,TN波动程度最高,其次为NH3-N;TP和CODMn均在0.2以下。除TN以外,其他参数的Qi在2015年以后多保持在3.0以内,达到Ⅲ类水体标准;水质类别最低时Qi分别接近Ⅳ类(CODMn,Qi=3.9)、Ⅴ类(NH3-N,Qi=4.6;TP,Qi=4.8),TN为劣Ⅴ类(Qi=8.0)。自2016年6月之后CODMn一直保持在Ⅱ类标准附近;氨氮则2015年出现Ⅳ类监测值后浓度出现明显下降,在2020年最新监测中接近Ⅰ类水体标准。

图3 康山2012—2020年单因子水质指数月际变化Fig.3 Monthly variation of single-factor water quality index at Kangshan from 2012 to 2020

2.1.1.2 蚌湖监测点单因子水质指数分析

对于代表长江江豚保护区的蚌湖监测点而言,平均水质略低于康山。NH3-N成为水质最优因子, 其次是TP和CODMn,TN最差。将标准差作为判断波动程度标准,CODMn(0.52)是波动幅度最小的参数,其次分别为TP(0.59)、NH3-N(0.98)、TN(1.51)。各参数最低类别监测值均出现在2016年1月份之前,其中TN单因子水质指数虽然仍在5.0左右波动,但自2018年1月起出现较为明显的浓度降低趋势;NH3-N自2016年1月份出现一次极值后,Qi一直保持在较为优良的水平,2016年2月—2020年12月均值为2.0;CODMn、TP均值分别为2.3、2.4,均满足Ⅲ类水体要求,如图4所示。

图4 蚌湖2012—2020年单因子水质指数月际变化Fig.4 Monthly variation of single-factor water quality index at Banghu from 2012 to 2020

2.1.1.3 都昌监测点单因子水质指数分析

如图5所示,对于代表长江江豚保护区的都昌监测点而言,将2012—2020年监测时段内的4种水质参数单因子水质指数均值从表现最优开始排布,依次为:CODMn(Qi=2.2)>TP(Qi=2.4)>NH3-N(Qi=2.6),均满足Ⅲ类水体标准限值;排在最后的依然是TN指标(Qi=4.9)氮素指标的单因子水质指数的波动范围较大,最优时接近或超过在Ⅱ类水体标准(NH3-N,Qi=1.2;TN(Qi=2.2);CODMn(Qi=1.5)、TP(Qi=1.5))最优值达到了Ⅱ类水体标准。时间趋势角度而言,而Qi最高的2014年2月时NH3-N、TN分别为7.7、8.4。NH3-N和CODMn在2016年后出现了一定程度的降幅;而TP和TN在2020年依旧保持前8 a的浓度水平。

图5 都昌2012—2020年单因子水质指数月际变化Fig.5 Monthly variation of single-factor water quality index at Duchang from 2012 to 2020

2.1.1.4 蛇山监测点单因子水质指数分析

如图6所示,对于代表鄱阳湖国家种质资源保护区的蛇山监测点而言,各水质参数单因子水质指数均值从表现最优开始排布,依次为:CODMn(Qi=2.2)> TP(Qi=2.7)> NH3-N(Qi=2.8),均满足Ⅲ类水体标准限值;与都昌监测点类似,排在最后的依然是TN指标(Qi=5.2,劣Ⅴ类)。将标准差作为判断波动程度标准,从大到小依次为:NH3-N(1.93)>TN(1.69)>TP(0.71)>CODMn(0.33)。NH3-N、CODMn、TP的Qi指数最好时优于Ⅱ类水质标准;而TN最优值也满足Ⅲ类水质标准限值(Qi=2.5)。NH3-N和TN在2015年之前浓度监测值较为同步,2014年2月均出现最大值(NH3-N,Qi=10.1;TN,Qi=11.0);进入2016年,NH3-N浓度出现明显下降,2016—2020年单因子指数为2.0,同期TN则为5.0。

图6 蛇山2012—2020年单因子水质指数月际变化Fig.6 Monthly variation of single-factor water quality index at Sheshan from 2012 to 2020

2.1.2 各监测点综合水质标识指数分析

鄱阳湖重点保护区各监测点WQII月际变化如图7所示。

图7 2012—2020年鄱阳湖重点保护区代表性监测点WQII月际变化Fig.7 Monthly variation of WQII at representative monitoring sites in key reserves of Poyang Lake from 2012 to 2020

从平均值来看,WQII值主要在3.01附近分布,分别为康山(2.91)>都昌(3.01)>蚌湖(3.11)>蛇山(3.31)。WQII高值集中出现在11、12、1、2月份共4个枯季月份。从时间趋势分析,都昌和蛇山站点在2014年第一季度以及2015年第一季度出现标识指数极大值,而后出现明显的下降,并在3.01左右波动;康山在2018年出现近8 a历史最高值之后下降到3.00以下,2020年末降至2.50左右。若将Qi≤3.0作为合格标准,各监测点达标率从最高开始排序,顺次如下:康山(69.44%)、都昌(62.04%)、蛇山(57.41%)、蚌湖(55.56%)。由于各保护区水质指数有部分样本超过3.00的情况(即高于Ⅲ类水质),可以将3.0都昌(2.60)>康山(2.70)>蛇山(2.80)。

2.2 SARIMA模型可靠性评价

如图8所示,qqnorm和qqline函数显示4个监测点SARIMA模型的残差均满足均值为0的正态分布,即正态性假设。同时,模型还需满足残差的自相关系数都应为0。在此选择Box.test函数,对立假设为“至少存在一个自相关系数不为0”,显著性水平设置为0.01。经残差检验,4个点位模型的p值分别为康山(0.85)、蚌湖(0.97)、都昌(0.33)、蛇山(0.95),均大于显著性水平。则否定对立假设,模型残差“自相关系数均为0”成立,即满足独立性。总体而言,SARIMA模型对各点位监测序列的拟合效果良好、可靠性高。

图8 各监测点SARIMA模型残差的正态分布Q-Q图Fig.8 Normal distribution Q-Q plots to determine whether the residuals of SARIMA model for each moni-toring point series satisfy the normality assumption

2.3 Chl-a月际浓度的预测优度比较

如表1所示,同步对鄱阳湖重点保护区各监测点建立以Chl-a为因变量的多元线性逐步回归方程(MLSR),并计算预测Chl-a与真实监测值的偏差,以比较2种模型拟合的精度。结合江西省水文局水质监测资料类别,考虑到藻类生长受水温、营养盐等要素影响明显,因此选择水温(Tw)、透明度(Ds)、

表1 鄱阳湖重点保护区各监测点位Chl-a的SARIMA模型与MLSR方程平均预测误差对比Table 1 Comparison of average prediction errors of Chl-a concentration between SARIMA model and MLSR equation for monitoring sites of key reserves in Poyang Lake

CODMn、NH3-N、TP、TN 6个因子作为备选自变量。为进一步比较逐步回归方程与SARIMA两种模型预测的可靠性,选择模型预测数据与监测数据之间的绝对误差(Absolute Percentage Error,AEP)、均方根误差(Root Mean Square Error,RMSE)、一致性指数(Index of Agreement, IoA)作为判断2种模型的精准度依据。AEP、RMSE是误差评估参数,参数值越小代表预测值与实测值差异小,模拟精准度更高;IoA综合体现了模式与观测的一致程度,其值越大代表一致度越高。各站点取误差平均值比较,如表2所示。

各监测点Chl-a浓度及最佳SARIMA模型如图9所示,中括号内数字代表周期s=12,即将12个月视为该序列的周期。蚌湖点位Chl-a浓度范围为0.5～51.0 μg/L,多年均值为11.01 μg/L,最高浓度出现在2018年10月,标准差为9.19 μg/L,波动范围与浓度最高。康山除2016年7月出现极值32.70 μg/L以外,Chl-a浓度保持在0.10～19.20 μg/L,标准差相较蚌湖降低5.17 μg/L,波动程度明显趋缓。都昌与蛇山点位Chl-a浓度与蚌湖相比大幅降低,浓度范围分别为0.20～18.50、0.57～23.40 μg/L;多年均值分别比蚌湖低5.12、5.86 μg/L;标准差分别为3.50、3.37 μg/L。

图9 基于SARIMA模型的各监测点Chl-a浓度预测及其放大Fig.9 SARIMA-based Chl-a predictions for monitoring sites and its zoomed-in plot

从模型精度角度比较,都昌最优模型ARIMA(0,1,1)(0,0,1)[12]的RMSE为2.09,比MLSR方程少3.01,预测优势明显。同样,蚌湖和蛇山SARIMA模型预测值的RMSE相较MLSR方程的预测值更小,IoA比MLSR方程预测值更大。8月份Chl-a在都昌和蛇山监测点的监测浓度出现明显的离群值,分别降至1.76、0.57 μg/L,均为监测序列中最小值,因此2020年8月2种预测模型均出现了较高的误差。虽然康山、都昌的IoA表现出MLSR方程略高于SARIMA模型,但从整体来看,SARIMA模型AEP、RMSE和IoA均值较MLSR方程分别低1.38、1.43和高0.05;除去2020年8月都昌和蛇山Chl-a浓度极小值后,SARIMA模型平均相对误差为17.14%,MLSR方程为34.09%。因此可以认为SARIMA模型精度更高,拟合效果更好。

3 讨 论

湖泊水体Chl-a的变化与工业与人类生活[26]、土地利用[27]、内源再释放[28]、水体局部理化条件[29]等多重作用因素有关。基于已有数据建立的MLSR方程拟合效果好,可以根据系数给出每个自变量的解释并消除多重共线性[30]。但逐步回归的线性预测结果往往不能全面考虑多种因素,且叶绿素与各要素驱动的水质、物化因子存在不可忽视的非线性关系[31-32]。SARIMA模型对含有季节因素的时间序列适用性较好,在考虑历史数据的优势以外还考虑了自变量的随机干扰[33]。从模型对各点位Chl-a浓度的预测精度来看,除都昌和蛇山点位于2020年8月的预测以外,SARIMA模型均优于MLSR方程,这也与Khodakhah等[34]和Chen等[35]在预测过程中所得出SARIMA模型精准度更高的结论一致。

本文选取的鄱阳湖区4个代表性水质监测点分布于鄱阳湖中部主体湖区的三大保护区内,除TN因子明显超标以外,自2016年开始总体水质改善明显。将各监测点之间进行比较,位于南部的康山最优,其次是北部的都昌和蚌湖,排在最后的蛇山,且水质较差的情况多集中在鄱阳湖“低水似河”的枯水期,这可能由于较长的水力停留水平有关[36]。但藻类生长并不完全与氮磷营养盐浓度呈正比,如水温这一控制藻类生长的重要因素在康山和蛇山体现较为显著。受风浪扰动和航运采砂影响,作为悬移质泥沙颗粒浓度的主要表征参数之一的SD是都昌点位藻类增殖主要驱动要素,悬移质泥沙颗粒因其特殊的组成与性质成为氮磷等营养的重要载体。

综合水质标识指数超过5.002的总计15次监测值全部来自于3个监测点,分别为蚌湖(2次)、都昌(3次)、蛇山(10次)。从标准差作为波动判断标准来看,蛇山同样也是波动幅度最大的监测点,在TN仍处于高位的前提下,NH3-N、TP指标在2020年1月的监测中单因子指数分别为7.8、3.5,进一步增加了蛇山这一点位的水质保护压力。考虑到其位于国家种质资源保护区的敏感性,有必要对其周边污染源溯源及优化识别、强化管控工作的接续开展。

4 结 论

基于2012—2020年监测数据,以综合水质标识指数法较为全面地分析鄱阳湖重点保护区代表性监测点水质;构建以Chl-a为因变量的MLSR方程与SARIMA模型,选取2020年6—8月进行预测并与实测值比较,以判断模型适用性与精度,结论具体如下:

(1)鄱阳湖重点保护区代表性监测点综合水质良好,平均值自最优开始排布,依次为:康山(2.91)>都昌(3.01)>蚌湖(3.11)>蛇山(3.31)。综合水质标识指数高值多集中出现在水力停留时间较长的枯季。

(2)SARIMA模型相比于MLSR方程能更好地预测保护区各点位Chl-a浓度,同时建议提高监测频次、强化序列数据可靠性。

(3)蚌湖点位Chl-a浓度最高、蛇山综合水质波动最为显著,建议开展两点位周边污染源溯源及优化识别、强化管控工作。