1、滇池水質時間序列變化分析1引言湖泊富營養(yǎng)化是當前我國水環(huán)境領域面臨的重要問題之一，其中，滇池作為高原重污染湖泊的典型代表，自1980s以來受到人們的廣泛關注，研究人員也對此開展了大量的監(jiān)測、模擬、規(guī)劃和控制研究.在長期的研究中，如何評估滇池的水質變化趨勢、識別主要水質指標的演替特征 與規(guī)律，一直是人們廣為關注的熱點問題之一(萬能等，2007;鄒銳等，2011).在國務院發(fā)布的滇池流域水污染防治“十二五”規(guī)劃中，提出了全面推進及突出重點、兼顧全面的原則.為更好地推進“十二五”期間滇池富營養(yǎng)化控制和水質改善，需基于長時間序列的水質數(shù)據(jù)分析，識別滇池水質指標的變化趨勢和長期水平，進而區(qū)分不同水質指標

2、在滇池污染防治中的優(yōu)先程度，從而可以更具針對性地進行滇池污染防治水質變化趨勢的識別并非簡單的時間變化分析，而要考慮到水質變化過程中固有的周期性和隨機性特征，排除干擾誤差.在水質趨勢的時間序列分析中，統(tǒng)計模型是常用的方法.目前已有的研究多采用線性回歸或者基于次序統(tǒng)計量的非參數(shù)方法，但因其主要基于線性或者單調性假設， 不能反映局部變化.而水質由于受到人為活動干擾及其他自然因素影響，并不滿足線性、單調假 設.為解決這一問題，在前期的研究基礎上，STL(Seasonal-Trend Decomposition using LOESS)方法被應用于水質評價中，它采用局部加權回歸法(LOESS)進行擬合，

3、是一種可以處理非線性、局部趨勢的非參數(shù)統(tǒng)計方法.STL方法最早由Clevel and 等(1990)提出并應用于對大氣CO2濃度和美國失業(yè)人口數(shù)變化趨勢分析上.在水質變化分析中，Qian等(2000)最先采用STL方法對美國加州紐斯河口的氮(N)、磷(P)數(shù)據(jù)進行了趨勢識別.此后，STL方法在環(huán)境領域得到廣泛應用， 例如，Sellinger 等(2008)應用STL方法分析了密歇根-休倫湖水位的變化趨勢 ；Conrad等(2004) 應用STL方法結合動態(tài)線性模型 (DLM)分析了美國亞德金河懸沙濃度和水流量的變化趨勢及關 系.作為探索性數(shù)據(jù)分析的有效手段，STL方法亦有廣泛的應用(Lu et

4、 al., 2003; Carslaw et al.,2005; Jong et al., 2012).對于滇池而言，由于人為干擾的強度增大，水質指標變化具有很強的非線性和隨機性特征.因此，本文擬采用STL方法對水質數(shù)據(jù)進行時間序列分析，剔除干擾因素，從而可以更為準確地反映各個水質指標的變化趨勢.但STL方法的缺陷在于無法有效判定趨勢變化的顯著性，為此，本文采用穩(wěn)態(tài)轉換指數(shù)(Regime Shift Index , RSI)對趨勢的變化進行顯著性檢驗，從統(tǒng)計學意義上確定趨勢變化的顯著性，以期為進一步的滇池水質改善提供決策參考.2研究對象與方法研究對象本文的分析對象為滇池外海，選取昆明市環(huán)境監(jiān)測

5、中心在外海的8個常規(guī)監(jiān)測點位(灰灣中、羅家營、觀音山西、觀音山中、觀音山東、白魚口、滇池南、?？谖?為研究對象(圖1).根據(jù)數(shù)據(jù)的可得性，選取水溫 (T)、pH、透明度(SD)、溶解氧(DO)、BOD5 CODMn氨氮(NH3-N)、總氮 (TN)、總磷(TP)、葉綠素a(Chl a)10個水質指標進行分析，時間尺度為1998 2010年，時間分辨率為月.因此，每個監(jiān)測點位的每個水質指標的數(shù)據(jù)樣本為156個(Chl a時間尺度為19992009年，每個監(jiān)測點位132個數(shù)據(jù)).數(shù)據(jù)缺失值比例為2.8%,采用中位數(shù)平滑方法進行插值；Q-Q圖(Q-Qplot)顯示插值后數(shù)據(jù)與原始數(shù)據(jù)具有相同的分布，

6、說明插值效果良好.本文對水質數(shù)據(jù)的分析均基于 R 3.0.1 版本( HYPERLINK / /).圖1滇池外海監(jiān)測點位分布首先采用STL方法對水質指標或水質指標的比例進行分解，并抽取分解后的趨勢項，探究指標或其比值的變化趨勢.STL方法可獲取趨勢項，但并不能對趨勢項變化特征進行分析.為此，本文采用穩(wěn)態(tài)轉化指數(shù)(RSI)對趨勢項變化的顯著性進行定量分析，探究趨勢項的突變和穩(wěn)定區(qū)域，從而對趨勢項的變化狀態(tài)進行確認.魯棒局部加權回歸法魯棒局部加權回歸法 (Robust LOESS)是一個迭代回歸的過程，是STL方法采用的平滑方法，其主要步驟如下(Clevel and , 1979;1988).LO

7、ESS 過程基于距離越近、相關性越強的假設，賦予不同位置的點不同的權重并進行局部加權回歸.該過程需要選定局部回歸的窗口長度、回歸方程階數(shù)(d)及權重函數(shù)(W),常采用立方權重函數(shù)：叭口二(1 一 W w W E ( |)京 I假設一個正整數(shù) qn時，令人n(x)為離x最遠點的距離，此時 q q(x)= n n(x)q n .魯棒性過程為了消除極端值對回歸結果的影響，基于 xi點處殘差Ri = g(xi)-yi 大小，賦予xi權重, 殘差值大的點處被賦予小的權重.通常采用平方權重函數(shù)：(2)H ) 士 J ( I 一尸.0 W u W(2)lo3 1令h=6x median(Ri)，則各點處的

8、Robust權重值為p i=B(Ri /h)，此權重與 vi(x) 一起用于 最小二乘法的參數(shù)估計.迭代過程重復LOESS過程和魯棒性過程，直至收斂 .STL方法STL是一種用LOESS乍為平滑器，將時間序列分解為低頻率的趨勢項、高頻率的周期項及不 規(guī)則變化的殘差項的非參數(shù)統(tǒng)計方法：Yr = Trend, + Seasonal + Residualr ( 3)式中，Yv、Trendv、Seasonalv和Residualv分別為v時刻的觀測值、趨勢項、周期項和殘 差項.對于水質數(shù)據(jù)，趨勢項可認為是低頻率的變化趨勢，周期項可認為是由于周期性穩(wěn)定擾動 造成的高頻變化， 而殘差項可認為是隨機擾動造成

9、的不規(guī)則變化，因此，將周期項和殘差項去除得到低頻的趨勢項，有利于準確認識水質變化趨勢.STL方法是一個遞歸的過程，每一次遞歸要分別進行3次LOESS滑動平均過程.魯棒局部加權回歸法方法的LOESSi程和魯棒性過程分別在STL的內部環(huán)(圖2)和外部環(huán)中嵌套實現(xiàn).滇池的水質數(shù)據(jù)是以年為周期的月時間序列，每個月份的數(shù)據(jù)組成1個子序列，共12個子序列.假設Yv、Skv、Tkv、Rkv分別表示v時刻的觀測值、k次迭代的周期項、趨勢項和殘差項， 則STL的內部環(huán)步驟如圖 2所示.對于波動劇烈的時間序列，還應加上外部環(huán)過程，即在內部環(huán) 完成n(i)次迭代前，根據(jù)每一時間點殘差的大小賦予魯棒權重(式(2),并

10、將權重應用于下一次迭代.外部環(huán)進行次數(shù)為 n(o),每進行一次外部環(huán)過程，都需要進行n(i)次內部環(huán)過程，初始Robust權重為1.平滑參數(shù)n(s)、n(l)、n(t)的選擇，內部環(huán)次數(shù)n(i)、外部環(huán)次數(shù) n(o)的選擇及收斂的規(guī)則參見文獻.穩(wěn)態(tài)轉換指數(shù)(RSI)RSI是一種滑動T檢驗(running T-test)方法，通過計算某一點的RSI值，來判斷該點前后M個值組成的樣本是否有顯著性的變化(Xu et al., 2013).t0 時刻的RSI值公式為：旦J 一(4) TOC o 1-5 h z RSI =。九)= “一)(4)F.Mi ，曠 1/訝 卜(*)-蒼+匯 22式中,10-1

11、1國二方 2 t(/)= w 2,x(t)表示t0時刻之前M個時間序列值的平均值，工2表示t0時刻之后M個時間序列值的平均值.如果RSI值不在對應自由度的T分布的置信區(qū)間中，則認為該點前后M個值組成的樣本有顯著的變化；如果某個區(qū)域都是顯著點，說明該區(qū)域變化顯著、迅速，而在顯著區(qū)域之間的區(qū) 域相對穩(wěn)定.RSI結果與STL趨勢項相比，在開始和末端各缺少M個值.在本文對滇池的水質數(shù)據(jù)分析中，選擇 M=24(Rosqvist et al., 2010),對應T分布的自由度為 (2M-2)即為46,選擇 95%置信區(qū)間作為顯著點的判斷標準.3結果與討論一般理化指標的趨勢分析與判定為便于分析并考慮到滇池的

12、富營養(yǎng)化特征，本文將水質指標分為 3類：一般理化指標、有機物指標和富營養(yǎng)化指標.一般理化指標包括 T、pH、SD和DO;有機物指標包含 CODM和BOD5由 于缺乏COD的數(shù)據(jù)，本文以 BOD5心ODM來近似表征水體的可生化性；富營養(yǎng)化指標包含 NH3-N、 TN TP和Chla.圖3展示了 4種理化指標的趨勢項及 SD和DO的RSI分析結果.圖3 一般理化水質指標趨勢及 RSI變化（指標趨勢圖中水平虛線為平均值，RSI圖中兩條水平虛線分別表示t檢馬金95洲信區(qū)間的端點值（下同）；豎直虛線為 RSI顯著區(qū)域的極值點，豎直 實線表示2003年1月時間點）由圖可知，各個站點水溫總體趨勢一致，呈現(xiàn)降

13、低-升高-降低-升高的趨勢，變動范圍在16.519.5 C之間，平均值在 15.517.5 C之間；水溫從2008年開始緩慢升高，至U 2010年底達 到13年平均水平，較 1998年水平低1 C左右；觀音山西和?？谖髌骄递^高.各個立點pH變化趨勢一致， 變動范圍在8.29.2之間，平土p pH在8.79.0之間;pH從1998 年1月至2001年5月下降，至2003年上升，波峰出現(xiàn)在 5月或9月，至2007年下降，波谷出 現(xiàn)在1月或5月，至2009年呈現(xiàn)上升趨勢，2009年之后維持在比較高的穩(wěn)定狀態(tài)，在9.2左右，為近13年來最高.已有的研究發(fā)現(xiàn)，湖泊中物種豐度在超過其最適pH范圍（6.0

14、8.5）時易于降低；劉春光等（2006）的研究也表明，藻類在pH=8.5時生長狀況最好.天然水體的pH主要決定于CO2 HCO-& CO2-3平衡體系中各組分的相對含量，Shapiro等（1990）提出高pH/低CO2條件有利于藍藻形成競爭優(yōu)勢的理論.滇池頻發(fā)的有毒藻類水華以微囊藻為優(yōu)勢種，陳建中等（2010）的研究表明：銅綠微囊藻的最適生長pH為8.09.0 ,當pH為9.09.5時，其生長量略有下降.由此可知，2009年之后滇池水體的 pH條件有利于銅綠微囊藻的競爭生長.從SD變化趨勢圖上看，各個站點的變化趨勢一致，在 2000年2月之前緩慢下降，至 2002 年8月（?？谖髟?003年1

15、月）升高并達到最大值，之后一直呈下降趨勢，并在 2010年底達到 13年來最低值.在空間分布上，以滇池南最高、灰灣中最低.分析RSI的結果可知，SD在2000年1月附近顯著上升，之后到 2004年1月附近顯著下降，期間為穩(wěn)定的高值狀態(tài)，之后又有顯 著下降的時間點，2008年1月之后為穩(wěn)定低值.從DO趨勢圖上看，各個站點DO值呈現(xiàn)先增加后減小的趨勢， 轉折點有所不同，觀音山西高 于其它站點.最低DO濃度大于6.0 mg - L-1 ,符合II類水質標準.水體中的DO濃度是水體復 氧過程和水體生物呼吸作用、水體植物光合作用的暫時平衡，受海拔、水溫、鹽度及耗氧有機物的分解速率等影響，DO濃度影響因子

16、的復雜性可能是導致各站點DO濃度差異的原因.從RSI圖上看，近幾年顯著下降區(qū)域超過上升區(qū)域， 結合趨勢圖可見滇池總體 DO濃度呈現(xiàn)顯著下降趨勢.另外，從2003年1月開始，DO和SD呈現(xiàn)同步下降趨勢，并且均達到近13來的最低水平.從SD和DO的RSI的結果可知，2003年以后存在大量的負值顯著點，說明這一階段趨勢下降迅 速且明顯.SD和DO的下降是湖泊富營養(yǎng)化的重要表征，這也說明2003 2010年間滇池富營養(yǎng)化有加重的趨勢.3.2有機物指標的趨勢分析與判定圖4展示了 CODMn BOD5及BOD5心ODM的趨勢和 RSI分析結果.從CODM趨勢圖上看，各 站點總體趨勢一致，平均值在8 mg

17、- L-1左右，站點間相差很小，2006年之前在低于平均值水平下先增加后減小再增加，2006年前后迅速增長，2008年后雖然各個站點的變化趨勢不同，但仍維持在近 13年的較高水平.1998 2010年的CODMAF均水平達到IV類水質標準，2010年底 水平為V類標準.從RSI圖上看，在2002年1月附近存在顯著下降的區(qū)域，在2006年37月附近存在顯著上升的區(qū)域，結合趨勢圖可見，CODM在2006年之前經歷了低于 13年平均值的較高水平和較低水平兩個穩(wěn)定狀態(tài)；2006年后，開始顯著上升，并維持在高于13年平均值的高值水平;13年總體趨勢顯著上升.圖4有機物指標趨勢及其 RSI變化BOD外各個

18、站點的總體變化趨勢一致， 2000年之前經歷了短暫的上升過程 （?？谖鹘档停?， 至2008年6月降低，之后上升，至 2010底接近13年的平均水平，達到 IV類水質標準.從RSI 圖上看，存在多處顯著下降的區(qū)域， 但2008年11月的上升較為顯著， 說明BOD5存在反彈趨勢.水體CO而一般均比 CODM值大（黃慧坤，2004）,以BOD5心ODM來衡量水體可生化性一般比BOD5心。皿到的數(shù)值大.從趨勢圖和RSI圖上看，除去滇池南存在顯著上升的趨勢外， 各個站 點的顯著變化點均為顯著降低，2008年之后的上升趨勢不顯著，說明該比例在顯著下降，2010年底其值大多在 0.4左右，可生化性差.從變化

19、趨勢分析，BOD5心ODM的下降是由于CODMn勺上 升和BOD5的下降引起的.點源是滇池流域的主要污染源，位于滇池北岸的昆明主城區(qū)是流域重污 染排水區(qū)（李躍勛等，2004）.昆明市自1991年開始建設污水處理廠，到“十一五”末已建成8座，點源處理能力不斷加強.污水處理廠主要采用 A2/O、ICEAS等工藝，對BOD5W良好的去除效 果，這是導致BOD5下降的主要原因.而對于COD的處理效果則較 BOD5差（余冬等，2008）;另外, 由于人口增長和城市擴張，產生大量生活污水，點源負荷增加，加之有些點源未經處理直接排入滇池（王紅梅等，2009;李躍勛等，2004）,導致湖體CODM憎加.由于滇

20、池外海水體可生化下降, 導致生物降解緩慢、水體污染物積累，也是導致CODMrfh高的原因.3.3富營養(yǎng)化狀態(tài)的趨勢變化識別圖5展示了 NH3-N、TN和NH3-N/TN的變化趨勢及 RSI結果.由圖可知，灰灣中 NH3-N平均 值顯著高于其它站點，且變化趨勢亦異于其它站點.灰灣中NH3-N濃度經歷1次顯著下降和3次顯著上升過程，至2009年達到最大值，符合III類水質標準；其他站點NH3-N濃度變化趨勢一致， 在平均值附近波動，2008年后呈現(xiàn)下降趨勢.圖5 NH3-N、TN及其比例年份趨勢及 RSI結果各個立點TN濃度總體變化趨勢一致， 并在2006年前后劃分為低值和高值 2部分，低值部分

21、除去觀音山東和滇池南兩個站點， 在2002年附近存在顯著下降區(qū)域， 高值部分在2006年附近存 在顯著上升區(qū)域，之后大部分站點在 2009年附近顯著下降，之后穩(wěn)定或者上升，未達到V類水質標準.各個站點的NH3-N/TN總體趨勢一致，呈現(xiàn)下降趨勢，且下降幅度大，說明滇池氮元素形態(tài) 分布發(fā)生了顯著變化，NH3-N在氮元素形態(tài)中的主導地位正逐漸降低.從RSI結果分析，部分站點在2002年附近存在短暫的顯著上升區(qū)域，之后都是顯著下降區(qū)域，說明下降趨勢明顯.根據(jù)對同期文獻（余冬等，2008）數(shù)據(jù)的分析計算，昆明市污水處理廠的NH3-N加權去除率為91.20%,而TN的加權去除率僅為 70.81%.圖6展

22、示了 TP、TN/TP和Chla的變化趨勢及 RSI結果.由圖可知，各個站點的 TP濃度總體 趨勢一致，經過短暫的上升（滇池南下降）后，開始迅速下降，至2003年開始在低值波動，從2009年開始上升，并接近或者超過13年的平均值，除?？谖骱突覟持?，其他站點均未達到V類水質標準.RSI圖的結果顯示，2002年附近區(qū)域 TP濃度顯著下降，有 4個站點2009年后的上升是顯 著的.IS * TP. TX. TP. Chi.時曜河 K、亮牝%* Eh 3 * Hhl d n。Mid q皿圖6 TP、TN/TP、Chla趨勢及其 RSI變化從TN/TP趨勢圖上看，各個站點總體變化趨勢一致，在 2009年

23、1月前后分為上升和下降2部分.從RSI結果看，上升階段存在多處顯著上升區(qū)域，而下降階段或僅有一個顯著區(qū)域，或不 顯著.根據(jù)前述分析，可以將TN/TP的變化趨勢分為 3個部分：2003年之前由于TP濃度迅速下降主導的迅速上升；至2009年由于TN濃度迅速上升主導的迅速上升；至2010年底由于TP濃度上升主導的下降.TN/TP值總體在518之間，從2001年開始在1018之間變動.Chla濃度總體呈現(xiàn)先下降后上升的趨勢，在2006年達到波谷，灰灣中 Chla濃度明顯高于其它站點.從RSI圖上看，下降過程中經歷了 2次顯著下降，在2007年附近存在顯著的升高趨勢 （除羅家營站）.Chla是表征水體中

24、藻類濃度的指標；有學者探討了氮磷比與藻類生長之間的關系，認為氮磷比是湖泊水華的限制因子，如 Smith提出適宜藍藻生長的 TN/TP值為1016,適合 真核藻類生長的比例為1623（胡鴻鈞，2011）;由此分析，理論上滇池的氮磷比應該適宜藍藻生長.但Chla的變化趨勢并不支持這種理論，說明該理論對于滇池并不適應.從TM TP和Chla的趨勢圖上可以看出，2003年之前TP和Chla的下降具有很好的協(xié)同性，而2006年之后TN和Chla的升高具有很好的一致性，說明滇池藻類的生長受到磷、氮的共同影響，但影響的時間階段不同，且近年來氮已成為滇池藻類爆發(fā)的限制性因素.這一結論與之前的研究相同（萬能等，2007;顏小 TOC o 1-5 h z 品等，2013）,但尚需細致的實驗數(shù)據(jù)進一步驗證綜上分析