闽江流域水质时间序列变化趋势识别及特征分析

王春晓，卢毅敏

(1.福州大学 空间数据挖掘与信息共享教育部重点实验室，福建 福州 350116；2.福州大学 地理空间信息技术国家地方联合工程研究中心，福建 福州 350002；3.福州大学 数字中国研究院(福建)，福建 福州 350002)

1 研究背景

地表水资源污染正严重威胁着人们的饮用水安全，同时制约着社会各个方面的发展，是发展过程中亟待解决的问题之一[1,2]。如何还原出水质随时间的变化趋势，并根据不同水质指标的变化趋势确定流域内不同污染需要治理的优先程度，对流域内水环境的针对性防治有重要的参考意义[3]。

对水质指标变化趋势进行识别有两大难点，即季节性[4]和随机性[5]特征。目前，常规的水质时序研究方法可分为基于统计的方法[6-7]以及线性回归法[8]，但这两种方法的缺陷是难以反映水质时间序列非线性和局部变化特征。STL法(seasonal and trend decomposition using Loess)是国外常用的时序数据研究方法，该方法最早应用于大气CO2浓度变化趋势分析[9]，其采用局部加权回归[10](Loess)策略拟合时间序列数据，最突出的特点是能够有效处理非线性序列，并且能够根据数据的特点反映出时间序列局部的变化特征[11]。Ndungu等[12]基于MODIS遥感数据并引入STL完成了肯尼亚奈瓦沙湖10年间叶绿素a的时空变化趋势分析，表明STL方法对长期的时间序列的拟合具有良好的适应性；Stow等[13]应用STL季节性趋势分解法来分析莫米河养分浓度和氨氮负荷的长期变化趋势；梁中耀等[14]在滇池流域水质指标时间序列数据分析中，利用STL方法准确评估了氨氮在总氮中随时间变化越来越占据主导地位的趋势，为滇池流域水环境污染评价、治理提供了一定的参考。尽管STL时序分解算法在数据趋势分析中可有效克服统计模型法和线性回归法存在的难以处理非线性序列数据、反映局部特征的问题，但STL方法不能直观地体现出趋势变化的显著性[14]。

基于前人的研究成果和不足之处，本研究以闽江流域为研究对象，提出一种结合STL时序分解法和Mann-Kendall趋势检验算法的水质指标识别及变化趋势分析策略，该策略在利用STL算法剔除人为和自然因素的干扰项及周期项后，保留能够反映水质指标变化的趋势项，同时结合Mann-Kendall趋势检验算法对分解出的趋势项进行显著性检验[15]，以此判断其变化显著性，为闽江流域水质防治提供决策参考。

2 研究区概况

闽江是福建的母亲河，发源于闽、赣边界，长约541 km，流域总面积约为6.1×104km2。闽江上游包含沙溪、富屯溪、建溪三大支流，中下游支流众多。流域内共包含7个地市的36个县，常驻人口近1 350万，约占全省人口的1/3，流域GDP约占全省1/4。闽江上游以山地为主，高山峡谷密布，河道坡降比较大，上游三大支流在南平附近汇集，面积共约42 000 km2，约占流域总面积的68.7%。

闽江流域水污染由于主要受人为因素的影响具有很强的地域性特征，其中闽江上游主要流经南平、三明市，污染主要来自于工业废水，年入江废水约4×108t，工业废水占比达90%；下游主要受到来自于福州市的点源、面源等污染，每年大约有2×108t废水排入闽江。

图1 闽江流域水系及水质监测站点分布图

3 数据来源与研究方法

3.1 数据来源

为全面识别和分析闽江流域水质状况及水质指标时序变化趋势，根据流域内主要干流、支流以及重点水库的位置，共选取12个水质自动监测站的水质监测指标时间序列数据作为实验数据源，其中8个站点位于闽江上游，4个站点位于闽江下游。水质自动监测站数据由福建省生态环境厅提供，共8个监测指标，分别为温度(TEMP，℃)、酸碱度(pH)、氨氮(NH3—N，mg/L)、电导率(CON，μS/cm)、浊度(TUR，NTU)、溶解氧(DO，mg/L)、高锰酸盐指数(CODMn，mg/L)、总磷(TP，mg/L)，时间跨度为2014年1月至2018年12月底，时间间隔为每天1次，12个站点共20 615条数据。由于某些不确定因素的影响(如网络波动、偶然误差等)，数据存在着一定的缺失和异常等情况。对于其中的异常数据，以《地表水环境质量标准》(GB 3838-2002)中的各水质指标能检测到的最值为范围，筛选出各水质时间序列中超出范围的数据(异常值)并赋值为零。对于缺失数据，使用多重插补法进行插补[16-17]，并以月时间尺度为长度，计算每个月的水质指标数据的平均值，得到共12×8×5组实验数据。

3.2 研究方法

3.2.1 STL时间序列分解算法 STL是一种非参数统计方法，它以具有鲁棒性的局部加权回归(robust loess)作为平滑方式，将时间序列数据以加法模型为原理分解为趋势项、周期项和不规则的余项[9]：

Yv=Tv+Sv+Rv

(1)

式中：Yv为时刻v的时间序列观测值;Tv，Sv，Rv分别为时刻v的趋势值、周期值和余项值。

局部加权回归(Loess)：根据点于点之间距离越小权重越大的原则，对各时间数据的离散点进行赋权值并采用局部回归法对时间序列数据拟合[18]。局部回归首先需要确定回归的阶数(d)、参与回归的序列长度(q)以及权重函数，设正整数q≤n(n为时序数据子序列个数)，并根据距离选取离x最近的q个点作为回归数据，对于每一个xi根据其与x的距离计算出权重，常采用的权重方法是如公式(2)所示的立方权重函数：

(2)

STL方法主要分为内循环和外循环两个部分：内循环的作用是进行趋势项和周期项的拟合；外循环的作用是调节时序数据点的稳健性权重，以消除异常值的影响，保证循环稳健性。

设Tv(k)、Sv(k)为第k-1次内循环结束时的趋势值和周期值，Tv(1)=0，且包含以下6个参数：

内循环次数n(i)；外循环稳健迭代次数n(o)；每一个周期时序数据样本点数量n(p)；分离周期项平滑参数n(s)；低通滤波平滑参数n(l)；趋势平滑参数n(t)，6个参数的选择规则参见参考文献[9]。

内循环具体步骤如下：

(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)；

外循环具体步骤如下：

(1)计算点i处的实际值与估计值之差：Bi=|g(xi)-yi|；

(2)计算xi点处健壮性权重大小(Bi越大，权重越小)，计算公式为：

(3)

设水质指标时序数据(x1，x2，…，xn)随机独立，并构造以下统计量S：

(4)

式中：n为水质指标时序数据数目;xj,xk分别为时序数据中第k个值和第k+1个值，Sign函数为：

(5)

当n>8时，其方差Var(S)为：

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

(6)

式中：tp为序列第p个数的值；q为具有相同数值序列的个数。

统计检验值Z的计算公式为：

(7)

Mann-Kendall法采用的是双边检验，当Z>0，表示序列有升高趋势；当Z<0，表示序列有降低趋势；若|Z|<1.96，表示变化不显著；若1.96<|Z|<2.56，变化显著；|Z|>2.56，变化极显著；且|Z|值越大，变化越显著。

例7：画面的色彩光线，还可以通过会声会影的各种滤镜来调节，比如利用“自动曝光”滤镜可以调节光线不足；综合调整“色彩平衡”、“色调和饱和度”、“亮度和对比度”等滤镜的自定义参数，改善画面存在的问题；运用“NewBlue色彩”滤镜，调节颜色、色彩、饱和、亮度、电影伽玛等参数，获得良好的视觉效果，增加作品的艺术魅力。比如校园的宣传视频中，春天的桃花，花朵更加色彩艳丽、娇艳欲滴；秋天的银杏，满树金黄的叶子在阳光的照耀下金光闪闪。

4 结果与分析

为了对闽江流域中各不同性质的水质指标时间序列进行准确地分析，本文依据水质指标的不同性质，将8个指标分为3组，即物理化学指标、有机物污染指标和富营养化指标。物理化学指标包括温度(TEMP，℃)、酸碱度(pH)、电导率(CON，μS/cm)、浊度(TUR，NTU)；有机物污染指标包括溶解氧(DO，mg/L)、高锰酸盐指数(CODMn，mg/L)；富营养化指标包括总磷(TP，mg/L)、氨氮(NH3—N，mg/L)。分别对3组数据进行STL时间序列分解与趋势项Mann-Kendall检验，结果如下：

4.1 物理化学指标实验结果

图2和表1为4个物理化学指标的时序分解和Mann-Kendall检验结果。结果显示pH值发生大幅度下降的有建阳将口站和福州原厝站，建阳将口站2015-2018年的pH值一直处于下降的趋势；各个站点的pH值变化范围总体在6～8之间，其中上游8个站点中pH值均值最小的站点为邵武和顺站，均值为6.84，其他7个站点水体pH值多大于7；下游4个站点的pH均值均小于7，除闽侯大樟溪外的3个站点均存在局部时间pH值小于6.5的情况。从表1检验结果可知，有3个站点的pH值发生了极显著上升，3个站点的pH值发生了极显著下降。

图2 2014-2018年闽江流域各监测站点水质物理化学指标STL分解结果

TEMP时间序列分解结果表明，10个站点的水温在2014-2018年均有一定的上升，各个站点的温度变化范围为19～24℃，其中，梅口悬索桥平均水温高于其他站点，顺昌富文站水温有较大幅度的上升；从变化趋势来看，12个站点中有10个站点的水温在2016-2018年间发生了具有相同趋势的显著上升，占比为83.3%，查询同期气温与日照强度资料得知，闽江流域2016-2018年间气温与日照强度较同期有明显的上升，并且降水量较往年有明显的减少。

CON时间序列分解结果表明，闽江上游各站点CON值除南平水汾桥站点外，平均值范围为46.76～82.81 μS/cm，均低于闽江下游4个站点；CON值除了受到带电荷物质浓度的影响以外，还受到温度的影响，从分解结果可以看到，2016-2018年间南平水汾桥及闽江下游共计5个站点的CON和TEMP具有一致的变化趋势，说明闽江下游CON值受温度影响较大。从表1检验结果可知，2个站点发生了显著或极显著的下降，7个站点发生了趋势极显著的上升。

表1 物理化学指标Mann-Kendall趋势检验结果

TUR值所衡量的是水体中所含颗粒物的浓度对光的散射程度的大小，通常TUR值越大，水体浑浊度越高。TUR时间序列分解结果表明闽侯大樟溪站点2014-2018年间TUR值持续升高，并且明显高于其他站点，说明水体浑浊状况持续加重。其中三明东牙溪站点TUR最低，并且保持稳定。从检验结果可知，5个站点出现了显著下降或者极显著下降的趋势，2个站点出现了极显著上升的趋势。

4.2 有机物污染指标实验结果

图3和表2分别为有机物污染指标DO、CODMn的时序分解结果和Mann-Kendall趋势检验结果。DO时间序列分解结果表明，清流安砂、邵武和顺、建阳将口3个站点的DO值自2016年开始一直处于升高的状态，闽江上游8个站点DO平均质量浓度在7.1～8.6 mg/L之间，各站点均符合Ⅱ类水质标准，下游4个站点DO平均质量浓度在6.2～8.0 mg/L之间，较上游值偏低，符合Ⅱ类水质标准。其中，闽侯竹岐站DO平均质量浓度最低，为6.2 mg/L，并且与福州原厝、古田黄田站在2017年以后有显著的下降趋势，有超过Ⅲ类水质的风险。除闽侯大樟溪外，闽江上游各站点DO平均浓度均大于闽江下游各站点。从表2趋势检验结果可知，12个站点中有7个站点的DO值呈现显著上升或者极显著上升，占总数的58.3%，2个站点的DO值发生了极显著下降，占总数的16.7%。

图3 2014-2018年闽江流域各监测站点水质有机物污染指标STL分解结果

CODMn时序分解和检验表明，12个站点中，有10个站点的CODMn质量浓度在2016-2017年附近有明显的下降，但是有8个站点在2018年附近有明显的上升；所有站点平均质量浓度均在1.67 ～2.52 mg/L之间，均符合Ⅱ类水质标准，闽江上游有6个站点的平均值高于闽江下游除闽侯大樟溪之外的3个站点。表2结果表明，6个站点的CODMn值存在着极显著的下降趋势。总体上，闽江流域CODMn污染2014-2018年整体上有好转的趋势，尤其是闽江下游，自2015年以来CODMn浓度一直维持在较低水平，但是2017-2018年间CODMn污染有加重的趋势，需要进一步加强治理。

表2 有机物污染指标Mann-Kendall趋势检验结果

4.3 富营养化指标实验结果

图4和表3分别展示了NH3—N、TP的时序分解结果及Mann-Kendall趋势检验结果。结果表明，闽江上游8个站点NH3—N浓度在平均值附近波动并符合水质Ⅰ、Ⅱ类标准，南平水汾桥NH3—N浓度自2014年开始出现明显的下降趋势，但仍然保持较高的浓度；闽江下游4个站点平均值比上游8个站点略高，但在2014-2015年附近，下游4个站点的NH3—N浓度有趋势相同的显著性下降，并且在2015年以后一直保持着较低的水准。从表3检验结果可得，闽江流域1个站点的NH3—N浓度检验结果为极显著上升，但依然维持在较低的水准，8个站点NH3—N浓度为极显著或显著下降。

图4 2014-2018年闽江流域各监测站点水质富营养化指标STL分解结果

TP时序分解结果表明，闽江流域12个监测站点的TP平均质量浓度范围在0.02～0.12 mg/L之间，共 11个站点水质符合Ⅰ、Ⅱ类水质标准，回垅大桥站2014-2018年之间平均TP浓度符合Ⅲ类水质标准，但下降趋势极为显著；下游4个站点水质TP浓度符合Ⅱ类水质标准，闽侯竹岐站TP浓度下降趋势明显。从表3检验结果可知，闽江流域12个站点中，TP浓度出现极显著性下降趋势的有7个站点，占比为58.3%，有2个站点有极显著上升或显著上升趋势；闽江下游有3个站点出现极显著的下降趋势，作为闽江下游最大的支流，闽侯大樟溪在2017年以后TP浓度有显著性的上升，有超过Ⅲ类水质的风险。总体来说，闽江流域TP指标有着明显的好转，但是仍然存在着局部水质TP浓度超过Ⅲ类水质的危险。

表3 富营养化指标Mann-Kendall趋势检验结果图

4.4 实验结果分析

实验结果表明，2014-2018年间，闽江流域12个监测站点的8个监测指标中，除TEMP以外，其他7个指标的平均值及变化趋势在闽江上下游都呈现出不同的态势。闽江上游的7个水质指标中，除CODMn以外，其他6个指标平均值均优于闽江下游，闽江流域以高污染、低附加值工业为主，共有重点污染企业2000余家，其中又以造纸、塑料、化工等高污染企业为主，闽江上游三明、南平两市“三高”产业在产业结构中占比重较大，大部分县区工业万元产值CODMn排放量均大于闽江流域平均值。2014年，三明市和南平市工业CODMn排放量分别占整个流域工业CODMn排放总量的45%和41.1%，工业废水排放量占上游废水排放总量的90%，过多的工业废水排放是导致闽江上游CODMn浓度大于闽江下游的主要原因。闽江下游主要流经福州市区，沿江两岸人口密集，工商业发达，GDP总量连年升高，污染主要来源于城镇废水污染和工业污染等，2014年闽江福州段城镇生活产生的CODMn和NH3—N分别占整个市区排放总量的47.9%和48.8%，大量未经处理的城镇污水排放是导致闽江下游6个指标劣与上游的主要原因。“十二五”与“十三五”时期，闽江流域三市积极推进水环境的治理与改善工作，优化产业结构，关闭了一大批严重污染水环境的“十小”项目与企业，并对现有的污水处理厂进行效率提升改造，污水日处理量显著提升，福州市污水处理率由2014年的87.7%提高到2018年的90.3%，CODMn与NH3—N、TP减排明显，并且各指标有显著的改善趋势，整体上水质有明显的改善，但DO浓度依然保持在较低的水平，尤其是闽江下游，并逐渐成为影响闽江下游水质等级的主导因素。

5 结 论

(1)采用STL时间序列分解算法与Mann-Kendall检验算法相结合对闽江流域12个水质监测站点2014-2018年的8个水质指标时间序列进行分析，研究发现，2014-2018年间，闽江流域整体上水质较好。

(2)闽江流域各污染指标都存在着显著变好的趋势，其中闽江下游NNH3—N与TP浓度发生了显著性的下降；闽江上游水体CODMn浓度较下游偏高，但闽江下游水体DO浓度较上游偏低，存在着下降的趋势，有超过Ⅲ类水质标准的风险，并逐渐成为影响闽江下游水质的主要因素。

(3)利用STL方法和Mann-Kendall方法相结合对闽江流域水质时间序列数据进行趋势识别和分析，可以直观地对剥离自然和人为干扰的水质时间序列进行观察，对受人为、自然影响因素较大的水质时间序列具有很好的适应性，能较好地揭示水质时间序列的非线性、非单调性和局部特征的规律。