閩江流域水質(zhì)時間序列變化趨勢識別及特征分析

王春曉，盧毅敏

(1.福州大學(xué) 空間數(shù)據(jù)挖掘與信息共享教育部重點實驗室，福建 福州 350116；2.福州大學(xué) 地理空間信息技術(shù)國家地方聯(lián)合工程研究中心，福建 福州 350002；3.福州大學(xué) 數(shù)字中國研究院(福建)，福建 福州 350002)

1 研究背景

地表水資源污染正嚴重威脅著人們的飲用水安全，同時制約著社會各個方面的發(fā)展，是發(fā)展過程中亟待解決的問題之一[1,2]。如何還原出水質(zhì)隨時間的變化趨勢，并根據(jù)不同水質(zhì)指標的變化趨勢確定流域內(nèi)不同污染需要治理的優(yōu)先程度，對流域內(nèi)水環(huán)境的針對性防治有重要的參考意義[3]。

對水質(zhì)指標變化趨勢進行識別有兩大難點，即季節(jié)性[4]和隨機性[5]特征。目前，常規(guī)的水質(zhì)時序研究方法可分為基于統(tǒng)計的方法[6-7]以及線性回歸法[8]，但這兩種方法的缺陷是難以反映水質(zhì)時間序列非線性和局部變化特征。STL法(seasonal and trend decomposition using Loess)是國外常用的時序數(shù)據(jù)研究方法，該方法最早應(yīng)用于大氣CO2濃度變化趨勢分析[9]，其采用局部加權(quán)回歸[10](Loess)策略擬合時間序列數(shù)據(jù)，最突出的特點是能夠有效處理非線性序列，并且能夠根據(jù)數(shù)據(jù)的特點反映出時間序列局部的變化特征[11]。Ndungu等[12]基于MODIS遙感數(shù)據(jù)并引入STL完成了肯尼亞奈瓦沙湖10年間葉綠素a的時空變化趨勢分析，表明STL方法對長期的時間序列的擬合具有良好的適應(yīng)性；Stow等[13]應(yīng)用STL季節(jié)性趨勢分解法來分析莫米河養(yǎng)分濃度和氨氮負荷的長期變化趨勢；梁中耀等[14]在滇池流域水質(zhì)指標時間序列數(shù)據(jù)分析中，利用STL方法準確評估了氨氮在總氮中隨時間變化越來越占據(jù)主導(dǎo)地位的趨勢，為滇池流域水環(huán)境污染評價、治理提供了一定的參考。盡管STL時序分解算法在數(shù)據(jù)趨勢分析中可有效克服統(tǒng)計模型法和線性回歸法存在的難以處理非線性序列數(shù)據(jù)、反映局部特征的問題，但STL方法不能直觀地體現(xiàn)出趨勢變化的顯著性[14]。

基于前人的研究成果和不足之處，本研究以閩江流域為研究對象，提出一種結(jié)合STL時序分解法和Mann-Kendall趨勢檢驗算法的水質(zhì)指標識別及變化趨勢分析策略，該策略在利用STL算法剔除人為和自然因素的干擾項及周期項后，保留能夠反映水質(zhì)指標變化的趨勢項，同時結(jié)合Mann-Kendall趨勢檢驗算法對分解出的趨勢項進行顯著性檢驗[15]，以此判斷其變化顯著性，為閩江流域水質(zhì)防治提供決策參考。

2 研究區(qū)概況

閩江是福建的母親河，發(fā)源于閩、贛邊界，長約541 km，流域總面積約為6.1×104km2。閩江上游包含沙溪、富屯溪、建溪三大支流，中下游支流眾多。流域內(nèi)共包含7個地市的36個縣，常駐人口近1 350萬，約占全省人口的1/3，流域GDP約占全省1/4。閩江上游以山地為主，高山峽谷密布，河道坡降比較大，上游三大支流在南平附近匯集，面積共約42 000 km2，約占流域總面積的68.7%。

閩江流域水污染由于主要受人為因素的影響具有很強的地域性特征，其中閩江上游主要流經(jīng)南平、三明市，污染主要來自于工業(yè)廢水，年入江廢水約4×108t，工業(yè)廢水占比達90%；下游主要受到來自于福州市的點源、面源等污染，每年大約有2×108t廢水排入閩江。

圖1 閩江流域水系及水質(zhì)監(jiān)測站點分布圖

3 數(shù)據(jù)來源與研究方法

3.1 數(shù)據(jù)來源

為全面識別和分析閩江流域水質(zhì)狀況及水質(zhì)指標時序變化趨勢，根據(jù)流域內(nèi)主要干流、支流以及重點水庫的位置，共選取12個水質(zhì)自動監(jiān)測站的水質(zhì)監(jiān)測指標時間序列數(shù)據(jù)作為實驗數(shù)據(jù)源，其中8個站點位于閩江上游，4個站點位于閩江下游。水質(zhì)自動監(jiān)測站數(shù)據(jù)由福建省生態(tài)環(huán)境廳提供，共8個監(jiān)測指標，分別為溫度(TEMP，℃)、酸堿度(pH)、氨氮(NH3—N，mg/L)、電導(dǎo)率(CON，μS/cm)、濁度(TUR，NTU)、溶解氧(DO，mg/L)、高錳酸鹽指數(shù)(CODMn，mg/L)、總磷(TP，mg/L)，時間跨度為2014年1月至2018年12月底，時間間隔為每天1次，12個站點共20 615條數(shù)據(jù)。由于某些不確定因素的影響(如網(wǎng)絡(luò)波動、偶然誤差等)，數(shù)據(jù)存在著一定的缺失和異常等情況。對于其中的異常數(shù)據(jù)，以《地表水環(huán)境質(zhì)量標準》(GB 3838-2002)中的各水質(zhì)指標能檢測到的最值為范圍，篩選出各水質(zhì)時間序列中超出范圍的數(shù)據(jù)(異常值)并賦值為零。對于缺失數(shù)據(jù)，使用多重插補法進行插補[16-17]，并以月時間尺度為長度，計算每個月的水質(zhì)指標數(shù)據(jù)的平均值，得到共12×8×5組實驗數(shù)據(jù)。

3.2 研究方法

3.2.1 STL時間序列分解算法 STL是一種非參數(shù)統(tǒng)計方法，它以具有魯棒性的局部加權(quán)回歸(robust loess)作為平滑方式，將時間序列數(shù)據(jù)以加法模型為原理分解為趨勢項、周期項和不規(guī)則的余項[9]：

Yv=Tv+Sv+Rv

(1)

式中：Yv為時刻v的時間序列觀測值;Tv，Sv，Rv分別為時刻v的趨勢值、周期值和余項值。

局部加權(quán)回歸(Loess)：根據(jù)點于點之間距離越小權(quán)重越大的原則，對各時間數(shù)據(jù)的離散點進行賦權(quán)值并采用局部回歸法對時間序列數(shù)據(jù)擬合[18]。局部回歸首先需要確定回歸的階數(shù)(d)、參與回歸的序列長度(q)以及權(quán)重函數(shù)，設(shè)正整數(shù)q≤n(n為時序數(shù)據(jù)子序列個數(shù))，并根據(jù)距離選取離x最近的q個點作為回歸數(shù)據(jù)，對于每一個xi根據(jù)其與x的距離計算出權(quán)重，常采用的權(quán)重方法是如公式(2)所示的立方權(quán)重函數(shù)：

(2)

STL方法主要分為內(nèi)循環(huán)和外循環(huán)兩個部分：內(nèi)循環(huán)的作用是進行趨勢項和周期項的擬合；外循環(huán)的作用是調(diào)節(jié)時序數(shù)據(jù)點的穩(wěn)健性權(quán)重，以消除異常值的影響，保證循環(huán)穩(wěn)健性。

設(shè)Tv(k)、Sv(k)為第k-1次內(nèi)循環(huán)結(jié)束時的趨勢值和周期值，Tv(1)=0，且包含以下6個參數(shù)：

內(nèi)循環(huán)次數(shù)n(i)；外循環(huán)穩(wěn)健迭代次數(shù)n(o)；每一個周期時序數(shù)據(jù)樣本點數(shù)量n(p)；分離周期項平滑參數(shù)n(s)；低通濾波平滑參數(shù)n(l)；趨勢平滑參數(shù)n(t)，6個參數(shù)的選擇規(guī)則參見參考文獻[9]。

內(nèi)循環(huán)具體步驟如下：

(1)去除趨勢項，v時刻的時序觀測值減去上一輪的趨勢值(Yv-Tv(k))；

(2)時序子序列平滑，利用Loess對子序列進行擬合，并且向子序列前后擴展一個周期，記為Cv(k+1)；

(3)平滑時序子序列Cv(k+1)低通濾波。低通濾波由n(p)組成，對Cv(k+1)分別進行窗口長度為n(p)、n(p)、3的滑動平滑，然后做d=1、q=n(l)的Loess回歸，得到Lv(k+1)；

(4)時序子序列去趨勢，Sv(k+1)=Cv(k+1)-Lv(k+1)；

外循環(huán)具體步驟如下：

(1)計算點i處的實際值與估計值之差：Bi=|g(xi)-yi|；

(2)計算xi點處健壯性權(quán)重大小(Bi越大，權(quán)重越小)，計算公式為：

(3)

設(shè)水質(zhì)指標時序數(shù)據(jù)(x1，x2，…，xn)隨機獨立，并構(gòu)造以下統(tǒng)計量S：

(4)

式中：n為水質(zhì)指標時序數(shù)據(jù)數(shù)目;xj,xk分別為時序數(shù)據(jù)中第k個值和第k+1個值，Sign函數(shù)為：

(5)

當n>8時，其方差Var(S)為：

Var(S)=(n(n-1)(2n+5)-

(6)

式中：tp為序列第p個數(shù)的值；q為具有相同數(shù)值序列的個數(shù)。

統(tǒng)計檢驗值Z的計算公式為：

(7)

Mann-Kendall法采用的是雙邊檢驗，當Z>0，表示序列有升高趨勢；當Z<0，表示序列有降低趨勢；若|Z|<1.96，表示變化不顯著；若1.96<|Z|<2.56，變化顯著；|Z|>2.56，變化極顯著；且|Z|值越大，變化越顯著。

例7：畫面的色彩光線，還可以通過會聲會影的各種濾鏡來調(diào)節(jié)，比如利用“自動曝光”濾鏡可以調(diào)節(jié)光線不足；綜合調(diào)整“色彩平衡”、“色調(diào)和飽和度”、“亮度和對比度”等濾鏡的自定義參數(shù)，改善畫面存在的問題；運用“NewBlue色彩”濾鏡，調(diào)節(jié)顏色、色彩、飽和、亮度、電影伽瑪?shù)葏?shù)，獲得良好的視覺效果，增加作品的藝術(shù)魅力。比如校園的宣傳視頻中，春天的桃花，花朵更加色彩艷麗、嬌艷欲滴；秋天的銀杏，滿樹金黃的葉子在陽光的照耀下金光閃閃。

4 結(jié)果與分析

為了對閩江流域中各不同性質(zhì)的水質(zhì)指標時間序列進行準確地分析，本文依據(jù)水質(zhì)指標的不同性質(zhì)，將8個指標分為3組，即物理化學(xué)指標、有機物污染指標和富營養(yǎng)化指標。物理化學(xué)指標包括溫度(TEMP，℃)、酸堿度(pH)、電導(dǎo)率(CON，μS/cm)、濁度(TUR，NTU)；有機物污染指標包括溶解氧(DO，mg/L)、高錳酸鹽指數(shù)(CODMn，mg/L)；富營養(yǎng)化指標包括總磷(TP，mg/L)、氨氮(NH3—N，mg/L)。分別對3組數(shù)據(jù)進行STL時間序列分解與趨勢項Mann-Kendall檢驗，結(jié)果如下：

4.1 物理化學(xué)指標實驗結(jié)果

圖2和表1為4個物理化學(xué)指標的時序分解和Mann-Kendall檢驗結(jié)果。結(jié)果顯示pH值發(fā)生大幅度下降的有建陽將口站和福州原厝站，建陽將口站2015-2018年的pH值一直處于下降的趨勢；各個站點的pH值變化范圍總體在6～8之間，其中上游8個站點中pH值均值最小的站點為邵武和順站，均值為6.84，其他7個站點水體pH值多大于7；下游4個站點的pH均值均小于7，除閩侯大樟溪外的3個站點均存在局部時間pH值小于6.5的情況。從表1檢驗結(jié)果可知，有3個站點的pH值發(fā)生了極顯著上升，3個站點的pH值發(fā)生了極顯著下降。

圖2 2014-2018年閩江流域各監(jiān)測站點水質(zhì)物理化學(xué)指標STL分解結(jié)果

TEMP時間序列分解結(jié)果表明，10個站點的水溫在2014-2018年均有一定的上升，各個站點的溫度變化范圍為19～24℃，其中，梅口懸索橋平均水溫高于其他站點，順昌富文站水溫有較大幅度的上升；從變化趨勢來看，12個站點中有10個站點的水溫在2016-2018年間發(fā)生了具有相同趨勢的顯著上升，占比為83.3%，查詢同期氣溫與日照強度資料得知，閩江流域2016-2018年間氣溫與日照強度較同期有明顯的上升，并且降水量較往年有明顯的減少。

CON時間序列分解結(jié)果表明，閩江上游各站點CON值除南平水汾橋站點外，平均值范圍為46.76～82.81 μS/cm，均低于閩江下游4個站點；CON值除了受到帶電荷物質(zhì)濃度的影響以外，還受到溫度的影響，從分解結(jié)果可以看到，2016-2018年間南平水汾橋及閩江下游共計5個站點的CON和TEMP具有一致的變化趨勢，說明閩江下游CON值受溫度影響較大。從表1檢驗結(jié)果可知，2個站點發(fā)生了顯著或極顯著的下降，7個站點發(fā)生了趨勢極顯著的上升。

表1 物理化學(xué)指標Mann-Kendall趨勢檢驗結(jié)果

TUR值所衡量的是水體中所含顆粒物的濃度對光的散射程度的大小，通常TUR值越大，水體渾濁度越高。TUR時間序列分解結(jié)果表明閩侯大樟溪站點2014-2018年間TUR值持續(xù)升高，并且明顯高于其他站點，說明水體渾濁狀況持續(xù)加重。其中三明東牙溪站點TUR最低，并且保持穩(wěn)定。從檢驗結(jié)果可知，5個站點出現(xiàn)了顯著下降或者極顯著下降的趨勢，2個站點出現(xiàn)了極顯著上升的趨勢。

4.2 有機物污染指標實驗結(jié)果

圖3和表2分別為有機物污染指標DO、CODMn的時序分解結(jié)果和Mann-Kendall趨勢檢驗結(jié)果。DO時間序列分解結(jié)果表明，清流安砂、邵武和順、建陽將口3個站點的DO值自2016年開始一直處于升高的狀態(tài)，閩江上游8個站點DO平均質(zhì)量濃度在7.1～8.6 mg/L之間，各站點均符合Ⅱ類水質(zhì)標準，下游4個站點DO平均質(zhì)量濃度在6.2～8.0 mg/L之間，較上游值偏低，符合Ⅱ類水質(zhì)標準。其中，閩侯竹岐站DO平均質(zhì)量濃度最低，為6.2 mg/L，并且與福州原厝、古田黃田站在2017年以后有顯著的下降趨勢，有超過Ⅲ類水質(zhì)的風(fēng)險。除閩侯大樟溪外，閩江上游各站點DO平均濃度均大于閩江下游各站點。從表2趨勢檢驗結(jié)果可知，12個站點中有7個站點的DO值呈現(xiàn)顯著上升或者極顯著上升，占總數(shù)的58.3%，2個站點的DO值發(fā)生了極顯著下降，占總數(shù)的16.7%。

圖3 2014-2018年閩江流域各監(jiān)測站點水質(zhì)有機物污染指標STL分解結(jié)果

CODMn時序分解和檢驗表明，12個站點中，有10個站點的CODMn質(zhì)量濃度在2016-2017年附近有明顯的下降，但是有8個站點在2018年附近有明顯的上升；所有站點平均質(zhì)量濃度均在1.67 ～2.52 mg/L之間，均符合Ⅱ類水質(zhì)標準，閩江上游有6個站點的平均值高于閩江下游除閩侯大樟溪之外的3個站點。表2結(jié)果表明，6個站點的CODMn值存在著極顯著的下降趨勢。總體上，閩江流域CODMn污染2014-2018年整體上有好轉(zhuǎn)的趨勢，尤其是閩江下游，自2015年以來CODMn濃度一直維持在較低水平，但是2017-2018年間CODMn污染有加重的趨勢，需要進一步加強治理。

表2 有機物污染指標Mann-Kendall趨勢檢驗結(jié)果

4.3 富營養(yǎng)化指標實驗結(jié)果

圖4和表3分別展示了NH3—N、TP的時序分解結(jié)果及Mann-Kendall趨勢檢驗結(jié)果。結(jié)果表明，閩江上游8個站點NH3—N濃度在平均值附近波動并符合水質(zhì)Ⅰ、Ⅱ類標準，南平水汾橋NH3—N濃度自2014年開始出現(xiàn)明顯的下降趨勢，但仍然保持較高的濃度；閩江下游4個站點平均值比上游8個站點略高，但在2014-2015年附近，下游4個站點的NH3—N濃度有趨勢相同的顯著性下降，并且在2015年以后一直保持著較低的水準。從表3檢驗結(jié)果可得，閩江流域1個站點的NH3—N濃度檢驗結(jié)果為極顯著上升，但依然維持在較低的水準，8個站點NH3—N濃度為極顯著或顯著下降。

圖4 2014-2018年閩江流域各監(jiān)測站點水質(zhì)富營養(yǎng)化指標STL分解結(jié)果

TP時序分解結(jié)果表明，閩江流域12個監(jiān)測站點的TP平均質(zhì)量濃度范圍在0.02～0.12 mg/L之間，共 11個站點水質(zhì)符合Ⅰ、Ⅱ類水質(zhì)標準，回垅大橋站2014-2018年之間平均TP濃度符合Ⅲ類水質(zhì)標準，但下降趨勢極為顯著；下游4個站點水質(zhì)TP濃度符合Ⅱ類水質(zhì)標準，閩侯竹岐站TP濃度下降趨勢明顯。從表3檢驗結(jié)果可知，閩江流域12個站點中，TP濃度出現(xiàn)極顯著性下降趨勢的有7個站點，占比為58.3%，有2個站點有極顯著上升或顯著上升趨勢；閩江下游有3個站點出現(xiàn)極顯著的下降趨勢，作為閩江下游最大的支流，閩侯大樟溪在2017年以后TP濃度有顯著性的上升，有超過Ⅲ類水質(zhì)的風(fēng)險。總體來說，閩江流域TP指標有著明顯的好轉(zhuǎn)，但是仍然存在著局部水質(zhì)TP濃度超過Ⅲ類水質(zhì)的危險。

表3 富營養(yǎng)化指標Mann-Kendall趨勢檢驗結(jié)果圖

4.4 實驗結(jié)果分析

實驗結(jié)果表明，2014-2018年間，閩江流域12個監(jiān)測站點的8個監(jiān)測指標中，除TEMP以外，其他7個指標的平均值及變化趨勢在閩江上下游都呈現(xiàn)出不同的態(tài)勢。閩江上游的7個水質(zhì)指標中，除CODMn以外，其他6個指標平均值均優(yōu)于閩江下游，閩江流域以高污染、低附加值工業(yè)為主，共有重點污染企業(yè)2000余家，其中又以造紙、塑料、化工等高污染企業(yè)為主，閩江上游三明、南平兩市“三高”產(chǎn)業(yè)在產(chǎn)業(yè)結(jié)構(gòu)中占比重較大，大部分縣區(qū)工業(yè)萬元產(chǎn)值CODMn排放量均大于閩江流域平均值。2014年，三明市和南平市工業(yè)CODMn排放量分別占整個流域工業(yè)CODMn排放總量的45%和41.1%，工業(yè)廢水排放量占上游廢水排放總量的90%，過多的工業(yè)廢水排放是導(dǎo)致閩江上游CODMn濃度大于閩江下游的主要原因。閩江下游主要流經(jīng)福州市區(qū)，沿江兩岸人口密集，工商業(yè)發(fā)達，GDP總量連年升高，污染主要來源于城鎮(zhèn)廢水污染和工業(yè)污染等，2014年閩江福州段城鎮(zhèn)生活產(chǎn)生的CODMn和NH3—N分別占整個市區(qū)排放總量的47.9%和48.8%，大量未經(jīng)處理的城鎮(zhèn)污水排放是導(dǎo)致閩江下游6個指標劣與上游的主要原因。“十二五”與“十三五”時期，閩江流域三市積極推進水環(huán)境的治理與改善工作，優(yōu)化產(chǎn)業(yè)結(jié)構(gòu)，關(guān)閉了一大批嚴重污染水環(huán)境的“十小”項目與企業(yè)，并對現(xiàn)有的污水處理廠進行效率提升改造，污水日處理量顯著提升，福州市污水處理率由2014年的87.7%提高到2018年的90.3%，CODMn與NH3—N、TP減排明顯，并且各指標有顯著的改善趨勢，整體上水質(zhì)有明顯的改善，但DO濃度依然保持在較低的水平，尤其是閩江下游，并逐漸成為影響閩江下游水質(zhì)等級的主導(dǎo)因素。

5 結(jié) 論

(1)采用STL時間序列分解算法與Mann-Kendall檢驗算法相結(jié)合對閩江流域12個水質(zhì)監(jiān)測站點2014-2018年的8個水質(zhì)指標時間序列進行分析，研究發(fā)現(xiàn)，2014-2018年間，閩江流域整體上水質(zhì)較好。

(2)閩江流域各污染指標都存在著顯著變好的趨勢，其中閩江下游NNH3—N與TP濃度發(fā)生了顯著性的下降；閩江上游水體CODMn濃度較下游偏高，但閩江下游水體DO濃度較上游偏低，存在著下降的趨勢，有超過Ⅲ類水質(zhì)標準的風(fēng)險，并逐漸成為影響閩江下游水質(zhì)的主要因素。

(3)利用STL方法和Mann-Kendall方法相結(jié)合對閩江流域水質(zhì)時間序列數(shù)據(jù)進行趨勢識別和分析，可以直觀地對剝離自然和人為干擾的水質(zhì)時間序列進行觀察，對受人為、自然影響因素較大的水質(zhì)時間序列具有很好的適應(yīng)性，能較好地揭示水質(zhì)時間序列的非線性、非單調(diào)性和局部特征的規(guī)律。