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方法在參數診斷方面尚缺乏定量化的方法，尚需開展深入的研究.
