滇池水质时间序列变化分析

滇池水质时间序列变化分析

1 引言

湖泊富营养化是当前我国水环境领域面临的重要问题之一，其中，滇池作为高原重污染湖泊的典型代表，自1980s以来受到人们的广泛关注，研究人员也对此开展了大量的监测、模拟、规划和控制研究.在长期的研究中，如何评估滇池的水质变化趋势、识别主要水质指标的演替特征与规律，一直是人们广为关注的热点问题之一(万能等，2007;邹锐等，2011).在国务院发布的《滇池流域水污染防治“十二五”规划》中，提出了全面推进及突出重点、兼顾全面的原则.为更好地推进“十二五”期间滇池富营养化控制和水质改善，需基于长时间序列的水质数据分析，识别滇池水质指标的变化趋势和长期水平，进而区分不同水质指标在滇池污染防治中的优先程度，从而可以更具针对性地进行滇池污染防治.

水质变化趋势的识别并非简单的时间变化分析，而要考虑到水质变化过程中固有的周期性和随机性特征，排除干扰误差.在水质趋势的时间序列分析中，统计模型是常用的方法.目前已有的研究多采用线性回归或者基于次序统计量的非参数方法，但因其主要基于线性或者单调性假设，不能反映局部变化.而水质由于受到人为活动干扰及其他自然因素影响，并不满足线性、单调假设.为解决这一问题，在前期的研究基础上，STL(Seasonal-Trend Decomposition using LOESS)方法被应用于水质评价中，它采用局部加权回归法(LOESS)进行拟合，是一种可以处理非线性、局部趋势的非参数统计方法.STL方法最早由Clevel and 等 (1990)提出并应用于对大气CO2浓度和美国失业人口数变化趋势分析上.在水质变化分析中，Qian等(2000)最先采用STL方法对美国加州纽斯河口的氮(N)、磷(P)数据进行了趋势识别.此后，STL方法在环境领域得到广泛应用，例如，Sellinger等(2008)应用 STL方法分析了密歇根-休伦湖水位的变化趋势;Conrad等(2004)应用 STL方法结合动态线性模型(DLM)分析了美国亚德金河悬沙浓度和水流量的变化趋势及关系.作为探索性数据分析的有效手段，STL方法亦有广泛的应用(Lu et al., 2003; Carslaw et al., 2005; Jong et al., 2012).对于滇池而言，由于人为干扰的强度增大，水质指标变化具有很强的非线性和随机性特征.因此，本文拟采用STL方法对水质数据进行时间序列分析，剔除干扰因素，从而可以更为准确地反映各个水质指标的变化趋势.但STL方法的缺陷在于无法有效判定趋势变化的显著性，为此，本文采用稳态转换指数(Regime Shift Index，RSI)对趋势的变化进行显著性检验，从统计学意义上确定趋势变化的显著性，以期为进一步的滇池水质改善提供决策参考.

2 研究对象与方法

2.1 研究对象

本文的分析对象为滇池外海，选取昆明市环境监测中心在外海的8个常规监测点位(灰湾中、罗家营、观音山西、观音山中、观音山东、白鱼口、滇池南、海口西)为研究对象(图 1).根据数据的可得性，选取水温(T)、pH、透明度(SD)、溶解氧(DO)、BOD5、CODMn、氨氮(NH3-N)、总氮(TN)、总磷(TP)、叶绿素a(Chl a)10个水质指标进行分析，时间尺度为1998—2010年，时间分辨率为月.因此，每个监测点位的每个水质指标的数据样本为156个(Chl a时间尺度为1999—2009年，每个监测点位132个数据).数据缺失值比例为2.8%，采用中位数平滑方法进行插值;Q-Q 图(Q-Qplot)显示插值后数据与原始数据具有相同的分布，说明插值效果良好.本文对水质数据的分析均基于R 3.0.1版本(/).

图 1 滇池外海监测点位分布

首先采用STL方法对水质指标或水质指标的比例进行分解，并抽取分解后的趋势项，探究指标或其比值的变化趋势.STL方法可获取趋势项，但并不能对趋势项变化特征进行分析.为此，本文采用稳态转化指数(RSI)对趋势项变化的显著性进行定量分析，探究趋势项的突变和稳定区域，从而对趋势项的变化状态进行确认.

2.2 鲁棒局部加权回归法

鲁棒局部加权回归法(Robust LOESS)是一个迭代回归的过程，是STL方法采用的平滑方法，其主要步骤如下(Clevel and ，1979;1988).

2.2.1 LOESS过程

基于距离越近、相关性越强的假设，赋予不同位置的点不同的权重并进行局部加权回归.该

过程需要选定局部回归的窗口长度、回归方程阶数(d)及权重函数(W)，常采用立方权重函数：

假设一个正整数q≤n(n为时间序列长度)，令距离x点最近的q个点被选择参与回归，λq(x)为距离x点第q远的点与x点处的距离， xi-x 为xi点与x点之间的距离，则x的临近值权重

公式为vi(x)=W(

).选定回归阶数d后，根据最小二乘法得到回归结果(x).当q>n时，令λn(x)为离x最远点的距离，此时λq(x)=λn(x)q n .

2.2.2 鲁棒性过程

为了消除极端值对回归结果的影响，基于xi点处残差 Ri = g(xi)-yi大小，赋予xi权重，残差值大的点处被赋予小的权重.通常采用平方权重函数：

令h=6×median(Ri)，则各点处的Robust权重值为ρi=B(Ri /h)，此权重与vi(x)一起用于最小二乘法的参数估计.

2.2.3 迭代过程

重复LOESS过程和鲁棒性过程，直至收敛.

2.3 STL方法

STL是一种用LOESS作为平滑器，将时间序列分解为低频率的趋势项、高频率的周期项及不规则变化的残差项的非参数统计方法：

式中，Yv、Trendv、Seasonalv和Residualv分别为v时刻的观测值、趋势项、周期项和残差项.对于水质数据，趋势项可认为是低频率的变化趋势，周期项可认为是由于周期性稳定扰动造成的高频变化，而残差项可认为是随机扰动造成的不规则变化，因此，将周期项和残差项去除得到低频的趋势项，有利于准确认识水质变化趋势.STL方法是一个递归的过程，每一次递归要分别进行3次LOESS和滑动平均过程.鲁棒局部加权回归法方法的LOESS过程和鲁棒性过程分别在STL的内部环(图 2)和外部环中嵌套实现.