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版本(http://www.r-project.org/).
图 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)和外部环中嵌套实现.
图 2 STL方法内部环过程
滇池的水质数据是以年为周期的月时间序列，每个月份的数据组成1个子序列，共12个子序列.假设Yv、Skv、Tkv、Rkv分别表示v时刻的观测值、k次迭代的周期项、趋势项和残差项，则STL的内部环步骤如图 2所示.对于波动剧烈的时间序列，还应加上外部环过程，即在内部环完成n(i)次迭代前，根据每一时间点残差的大小赋予鲁棒权重(式(2))，并将权重应用于下一次迭代.外部环进行次数为n(o)，每进行一次外部环过程，都需要进行n(i)次内部环过程，初始Robust权重为1.平滑参数n(s)、n(l)、n(t)的选择，内部环次数n(i)、外部环次数n(o)的选择及收敛的规则参见文献.
2.4 稳态转换指数(RSI)
RSI是一种滑动T检验(running T-test)方法，通过计算某一点的RSI值，来判断该点前后M个值组成的样本是否有显著性的变化(Xu et al., 2013).t0时刻的RSI值公式为：
式中，
，x(t)表示t0时刻之前M个时间序列值的平均值，
2表示t0时刻之后M个时间序列值的平均值.如果RSI值不在对应自由度的T分布的置信区间中，则认为该点前后M个值组成的样本有显著的变化;如果某个区域都是显著点，说明该区域变化显著、迅速，而在显著区域之间的区域相对稳定.RSI结果与STL趋势项相比，在开始和末端各缺少M个值.在本文对滇池的水质数据分析中，选择M=24(Rosqvist et al., 2010)，对应T分布的自由度为(2M-2)即为46，选择95%置信区间作为显著点的判断标准.
3 结果与讨论
3.1 一般理化指标的趋势分析与判定
为便于分析并考虑到滇池的富营养化特征，本文将水质指标分为3类：一般理化指标、有机物指标和富营养化指标.一般理化指标包括T、pH、SD和DO;有机物指标包含CODMn和BOD5，由于缺乏COD的数据，本文以BOD5/CODMn来近似表征水体的可生化性;富营养化指标包含NH3-N、TN、TP和Chla.图 3展示了4种理化指标的趋势项及SD和DO的RSI分析结果.
图 3 一般理化水质指标趋势及RSI变化(指标趋势图中水平虚线为平均值，RSI图中两条水平虚线分别表示t检验95%置信区间的端点值(下同);竖直虚线为RSI显著区域的极值点，竖直实线表示2003年1月时间点)
由图可知，各个站点水温总体趋势一致，呈现降低-升高-降低-升高的趋势，变动范围在16.5~19.5 ℃之间，平均值在15.5~17.5 ℃之间;水温从2008年开始缓慢升高，到2010年底达到13年平均水平，较1998年水平低1 ℃左右;观音山西和海口西平均值较高.
各个站点pH变化趋势一致，变动范围在8.2~9.2之间，平均pH在8.7~9.0之间;pH从1998年1月至2001年5月下降，至2003年上升，波峰出现在5月或9月，至2007年下降，波谷出现在1月或5月，至2009年呈现上升趋势，2009年之后维持在比较高的稳定状态，在9.2左右，为近13年来最高.已有的研究发现，湖泊中物种丰度在超过其最适pH范围(6.0~8.5)时易于降低;刘春光等(2006)的研究也表明，藻类在pH=8.5时生长状况最好.天然水体的pH主要决定于CO2、HCO-3、CO2-3平衡体系中各组分的相对含量，Shapiro等(1990)提出高pH/低CO2条件有利于蓝藻形成竞争优势的理论.滇池频发的有毒藻类水华以微囊藻为优势种，陈建中等(2010)的研究表明：铜绿微囊藻的最适生长pH为8.0~9.0，当pH为9.0~9.5时，其生长量略有下降.由此可知，2009年之后滇池水体的pH条件有利于铜绿微囊藻的竞争生长.
从SD变化趋势图上看，各个站点的变化趋势一致，在2000年2月之前缓慢下降，至2002年8月(海口西在2003年1月)升高并达到最大值，之后一直呈下降趋势，并在2010年底达到13年来最低值.在空间分布上，以滇池南最高、灰湾中最低.分析RSI的结果可知，SD在2000年1月附近显著上升，之后到2004年1月附近显著下降，期间为稳定的高值状态，之后又有显著下降的时间点，2008年1月之后为稳定低值.
从DO趋势图上看，各个站点DO值呈现先增加后减小的趋势，转折点有所不同，观音山西高于其它站点.最低DO浓度大于6.0 mg · L-1，符合II类水质标准.水体中的DO浓度是水体复氧过程和水体生物呼吸作用、水体植物光合作用的暂时平衡，受海拔、水温、盐度及耗氧有机物的分解速率等影响，DO浓度影响因子的复杂性可能是导致各站点DO浓度差异的原因.从RSI图上看，近几年显著下降区域超过上升区域，结合趋势图可见滇池总体DO浓度呈现显著下降趋势.
另外，从2003年1月开始，DO和SD呈现同步下降趋势，并且均达到近13来的最低水平.从SD和DO的RSI的结果可知，2003年以后存在大量的负值显著点，说明这一阶段趋势下降迅速且明显.SD和DO的下降是湖泊富营养化的重要表征，这也说明2003—2010年间滇池富营养化有加重的趋势.
3.2 有机物指标的趋势分析与判定
图 4展示了CODMn、BOD5及BOD5/CODMn的趋势和RSI分析结果.从CODMn趋势图上看，各站点总体趋势一致，平均值在8 mg · L-1左右，站点间相差很小，2006年之前在低于平均值水平下先增加后减小再增加，2006年前后迅速增长，2008年后虽然各个站点的变化趋势不同，但仍维持在近13年的较高水平.1998—2010年的CODMn平均水平达到IV类水质标准，2010年底水平为V类标准.从RSI图上看，在2002年1月附近存在显著下降的区域，在2006年3—7月附近存在显著上升的区域，结合趋势图可见，CODMn在2006年之前经历了低于13年平均值的较高水平和较低水平两个稳定状态;2006年后，开始显著上升，并维持在高于13年平均值的高值水平;13年总体趋势显著上升.
图 4 有机物指标趋势及其RSI变化
BOD5在各个站点的总体变化趋势一致，2000年之前经历了短暂的上升过程(海口西降低)，至2008年6月降低，之后上升，至2010底接近13年的平均水平，达到IV类水质标准.从RSI图上看，存在多处显著下降的区域，但2008年11月的上升较为显著，说明BOD5存在反弹趋势.
水体COD值一般均比CODMn值大(黄慧坤，2004)，以BOD5/CODMn来衡量水体可生化性一般比BOD5/COD得到的数值大.从趋势图和RSI图上看，除去滇池南存在显著上升的趋势外，各个站点的显著变化点均为显著降低，2008年之后的上升趋势不显著，说明该比例在显著下降，2010年底其值大多在0.4左右，可生化性差.从变化趋势分析，BOD5/CODMn的下降是由于CODMn的上升和BOD5的下降引起的.点源是滇池流域的主要污染源，位于滇池北岸的昆明主城区是流域重污染排水区(李跃勋等，2004).昆明市自1991年开始建设污水处理厂，到“十一五”末已建成8座，点源处理能力不断加强.污水处理厂主要采用A2/O、ICEAS等工艺，对BOD5有良好的去除效果，这是导致BOD5下降的主要原因.而对于COD的处理效果则较BOD5差(余冬等，2008);另外，由于人口增长和城市扩张，产生大量生活污水，点源负荷增加，加之有些点源未经处理直接排入滇池(王红梅等，2009;李跃勋等，2004)，导致湖体CODMn增加.由于滇池外海水体可生化下降，导致生物降解缓慢、水体污染物积累，也是导致CODMn升高的原因.
3.3 富营养化状态的趋势变化识别
图 5展示了NH3-N、TN和NH3-N/TN的变化趋势及RSI结果.由图可知，灰湾中NH3-N平均值显著高于其它站点，且变化趋势亦异于其它站点.灰湾中NH3-N浓度经历1次显著下降和3次显著上升过程，至2009年达到最大值，符合III类水质标准;其他站点NH3-N浓度变化趋势一致，在平均值附近波动，2008年后呈现下降趋势.
图 5 NH3-N、TN及其比例年份趋势及RSI结果
各个站点TN浓度总体变化趋势一致，并在2006年前后划分为低值和高值2部分，低值部分除去观音山东和滇池南两个站点，在2002年附近存在显著下降区域，高值部分在2006年附近存在显著上升区域，之后大部分站点在2009年附近显著下降，之后稳定或者上升，未达到V类水质标准.
各个站点的NH3-N/TN总体趋势一致，呈现下降趋势，且下降幅度大，说明滇池氮元素形态分布发生了显著变化，NH3-N在氮元素形态中的主导地位正逐渐降低.从RSI结果分析，部分站点在2002年附近存在短暂的显著上升区域，之后都是显著下降区域，说明下降趋势明显.根据对同期文献(余冬等，2008)数据的分析计算，昆明市污水处理厂的NH3-N加权去除率为91.20%，而TN的加权去除率仅为70.81%.
图 6展示了TP、TN/TP和Chla的变化趋势及RSI结果.由图可知，各个站点的TP浓度总体趋势一致，经过短暂的上升(滇池南下降)后，开始迅速下降，至2003年开始在低值波动，从2009年开始上升，并接近或者超过13年的平均值，除海口西和灰湾中，其他站点均未达到V类水质标准.RSI图的结果显示，2002年附近区域TP浓度显著下降，有4个站点2009年后的上升是显著的.
图 6 TP、TN/TP、Chla趋势及其RSI变化
从TN/TP趋势图上看，各个站点总体变化趋势一致，在2009年1月前后分为上升和下降2部分.从RSI结果看，上升阶段存在多处显著上升区域，而下降阶段或仅有一个显著区域，或不显著.根据前述分析，可以将TN/TP的变化趋势分为3个部分：2003年之前由于TP浓度迅速下降主导的迅速上升;至2009年由于TN浓度迅速上升主导的迅速上升;至2010年底由于TP浓度上升主导的下降.TN/TP值总体在5~18之间，从2001年开始在10~18之间变动.
Chla浓度总体呈现先下降后上升的趋势，在2006年达到波谷，灰湾中Chla浓度明显高于其它站点.从RSI图上看，下降过程中经历了2次显著下降，在2007年附近存在显著的升高趋势(除罗家营站).Chla是表征水体中藻类浓度的指标;有学者探讨了氮磷比与藻类生长之间的关系，认为氮磷比是湖泊水华的限制因子，如Smith提出适宜蓝藻生长的TN/TP值为10~16，适合真核藻类生长的比例为16~23(胡鸿钧，2011);由此分析，理论上滇池的氮磷比应该适宜蓝藻生长.但Chla的变化趋势并不支持这种理论，说明该理论对于滇池并不适应.从TN、TP和Chla的趋势图上可以看出，2003年之前TP和Chla的下降具有很好的协同性，而2006年之后TN和Chla的升高具有很好的一致性，说明滇池藻类的生长受到磷、氮的共同影响，但影响的时间阶段不同，且近年来氮已成为滇池藻类爆发的限制性因素.这一结论与之前的研究相同(万能等，2007;颜小品等，2013)，但尚需细致的实验数据进一步验证.
综上分析，STL方法对于滇池水质数据时间序列趋势识别和特征分析具有很好的适应性.从分析结果看，在原始数据的基础上剥离了周期性干扰和随机扰动后，趋势项用平滑曲线表示，很好地描述了水质指标变化所具有的非线性、非单调性特征及局部性特征.此外，上述结果还表明，该方法适用于受人为干扰较为强烈的水质指标的时间序列分析，如在水体可生化性的分析中，就较好地揭示了污水处理等人为干扰对水质的影响.具体参见 污水处理技术资料或污水技术资料更多相关技术文档
4 结论
1)本文采用STL方法和RSI方法对滇池外海10个水质指标、13年的时间序列变化进行了趋势分析和判定，研究发现，滇池各污染物指标在外海各个监测点位的变化趋势基本一致，浓度有所不同，灰湾中水质最差，滇池南和海口西水质最好.
2)从各种物质的变化趋势及总体水平来看，TN浓度升高及SD和DO的协同下降表明1998—2010年间滇池富营养化状态加剧，而TP亦处于劣V类水质水平;BOD5和CODMn分别符合IV和V类水质标准，BOD5总体趋势下降，CODMn显著上升.根据各水质指标的变化趋势及长期水平，“十二五”期间滇池应重点防治氮、磷营养盐，兼顾有机物污染.
3)滇池湖体的部分污染物分布发生了显著变化.有机物可生化性下降应引起高度重视，且由此而引起的CODMn累积是滇池有机物污染治理需要考虑的问题;NH3-N在氮元素形态中的主导地位正逐渐降低，亟待对氮元素的形态变化状况和机理进行深入研究，提出更合适的防治方案.
4)本文提供了一种通过剥离扰动来分析周期性环境时间序列数据的方法，去除周期项和残差项对于深入分析数据在时间序列上的总体变化趋势具有重要意义，也具有延伸推广的价值.本文据此对滇池外海的水质变化趋势及部分指标间的定性关系进行了识别，下一步仍需加强对趋势项的定量关系分析;此外，STL方法在参数诊断方面尚缺乏定量化的方法，尚需开展深入的研究.