基于LOADEST模型和小波变换的河流氮磷污染动态分析

金亚楠，张柏发，郝 韵，邬建红，吕 军

(浙江大学 环境与资源学院，浙江 杭州310058)

非点源污染致使水体富营养化，是自然水体水质退化的主要原因[1]。尽管流域氮、磷等营养物质的非点源污染过程十分复杂，但最终都会反映在河流的水质变化上。水中氮、磷污染物浓度变化的直接监测和污染物通量估算，是进行河流水质评价和实施水污染控制的基础工作。目前，我国大多数河流都已建立了连续的流量监测系统，但水质的自动化监测尚未普及。一些已有的在线自动监测装置，也大多存在着受环境变化影响较大、监测阈值限制，以及安装、使用和维护的成本较高等问题[2]。传统的人工采样和水质测定方法虽然仍在水质监测中广泛应用，但由于耗时费力，很难获得连续的或高频的水质监测资料；离散的每月1～2次的水质采样分析结果，却又难以反映河流水质的连续变化。LOADEST模型是美国地质调查局(USGS)为此专门开发的，用于根据河流离散的水文水质监测资料建立营养物负荷量与水流通量之间关系的多元统计模型[3-5]。应用该模型就可以在每日流量资料的基础上估算氮、磷负荷量在河流中的连续变化。这一方法已被国内外广泛采用[6-8]。

LOADEST模型是优化的统计回归模型，而统计方法的基本假设是样本数据在时空上有独立性。但实际上，驱动河流水质变化的主要因素(水文气象因素和氮、磷的投、排放变化等)具有明显的时间上的周期性和延续性，以及空间上的传承性和加和性。小波变换是在序列数据的时间-频率分析领域迅速发展的一种新技术[9]，它可以通过解构序列数据，由粗到细地描述各种复杂信号的具体结构，从而实现对整个数据序列的局部化特征和多层次变化规律的分析。本研究在对浙江东部典型农业流域长乐江出口断面进行每月一次水质监测的基础上，分析其2003—2016年间河流水文水质的变化特征，建立LOADEST模型来核定河流总氮(TN)和总磷(TP)的每日负荷量，并采用小波分析方法，进一步探索河流水质和营养物负荷通量的动态变化，以期揭示我国东部地区农业流域河流非点源污染变化规律及其时间格局特征，为农业流域河流水质控制提供科学依据。

1 材料与方法

1.1 研究区概况

长乐江流域位于浙江省绍兴市嵊州市(120°35′56″～120°49′03″E，29°27′98″～29°35′12″N)，是我国东部地区典型的低山河谷平原农业流域，流域总面积约为653 km2，多年平均降水量1 307.4 mm，但全年约75%的降水集中在4—9月，多年平均气温17.48 ℃。土地利用类型主要有林地(46.14%)、水田(21.04%)、园地(18.30%)、旱地(5.87%)、人居地(5.71%)、水域(2.31%)和其他(0.62%)。流域内的主干河流——长乐江是曹娥江上游，全长70.5 km，平均河宽55 m，河道比降为3.6‰，多年平均流量为16.39 m3·s-1，输沙量和侵蚀模数分别为10.9万t·a-1和127 t·km-2·a-1。

1.2 水质监测方法与水文气象数据来源

本研究主要分析长乐江流域出口断面监测点(图1中S7)的氮、磷营养物变化。水样采集和监测分析时间为2003年8月—2016年12月，采样监测频率为每月1次。在监测断面河流主流道距水面30 cm左右位置采集水样，装入聚乙烯瓶，用稀硫酸酸化后，当天运回实验室就TN、TP浓度等水质指标进行分析测定。TN浓度参照GB 11894—89，采用碱性过硫酸钾消化-紫外分光光度法测定；TP浓度参照GB 11893—89，采用过硫酸钾消化-钼锑抗比色法测定。

图1 长乐江流域和水质监测点的位置示意图Fig.1 Location of Changle River watershed and water quality monitoring points

河流日流量数据由浙江省水文局提供；气象数据来自于中国气象科学数据共享服务网(http://www.escience.gov.cn/metdatapageindex.html)。

1.3 LOADEST模型

USGS专门开发并发布了基于河流流量变化与水质变化关系的河流污染物负荷量估算和优化模型，即LOADEST模型和相应的软件[3]。LOADEST模型的主要思路是河流营养物的负荷量是河流流量的函数，并在不同的条件下伴随着不同的时间周期性变化。LOADEST模型包含11个河流营养物负荷量与流量关系的多元回归方程，基于水质监测得到的营养物负荷量与流量数据进行各方程的拟合，通过赤池信息量准则(AIC)[10-11]和Schwarz后验概率准则(SPPC)[12]比较各方程的拟合效果，最后采用具有最小AIC值和SPPC值的方程作为最优的描述河流污染物负荷量与流量关系的回归方程。

LOADEST模型软件提供了3种方程参数的拟合估计方法，即最小方差无偏估计(MVUE)[13]、渐进极大似然估计(AMLE)[14]和最小绝对偏差估计(LAD)[15]，可对非正态分布数据资料、删失型数据资料(具有不完整的或异常的序列数据资料)和变量共线性问题进行规范、方便的处理。

1.4 小波分析方法

小波分析的基本思想是，一个由变化的序列数据组成的函数(母函数)，可以被分解为多组不同变频的子函数来表示或逼近，从而可以通过对子函数的分析来表征母函数变化的细节。比如河流水体中的营养物负荷量的变化(母函数)，可以通过正交小波变换分解为不同变频的子函数组合，从而通过对子函数变化的特殊纹理和突变信息的分析，揭示负荷量变化的细节和规律。因此，小波分析又被称为数字信号的“数学放大镜”。小波分析的基本原理是通过增加或减小伸缩尺度(a)来得到信号的低频或高频信息，然后分析信号的概貌或细节，从而实现对信号不同时间尺度和空间局部特征的分析。在实际研究中，最主要的就是要由小波变换方程得到小波系数，然后通过这些系数来分析时间序列的时频变化特征[16]。通过小波变换分析得到的小波系数，可以描述不同时期的信号做出周期性变化(信号振荡)所需要的时间跨度(时间尺度)。如果要更深入地了解各时期信号振荡周期的整体性规律，就需要借助方差分析的方法[17]。将小波系数的平方差在b域上进行积分，就可得到小波方差。小波方差随尺度a的变化过程，称为小波方差图。小波方差图可用来描述信号在不同时间尺度扰动的相对强度，或者是不同振荡强度所对应的主要时间跨度(尺度)，即信号变化(振荡)的主次周期。

2 结果与分析

2.1 河流水文水质特征

2003—2016年间，长乐江出口断面的多年平均流量为16.39 m3·s-1，河流年平均流量呈显著增加趋势。河流水质监测结果表明，实测水体中的总氮和总磷浓度平均值分别为(3.74±0.91)mg·L-1和(0.136±0.086)mg·L-1，远高于该河流总氮和总磷浓度的参照标准(即未受人为污染时的河流水质指标参考标准：总氮1.70 mg·L-1，总磷0.055 mg·L-1)[18]。从年内的分布来看，月平均流量的变异较大，流量最大的6月份的多年平均流量为枯水期(10月—次年2月)月平均流量的2～3倍(图2-a)。

图2-b是研究期内河流监测断面上TN浓度的月平均分布状况，可以看出，枯水季节河流TN浓度具有高于丰水季节的趋势，这与我国南方其他地区的研究结果一致[19]，意味着多雨季节径流对河流的稀释作用强于对土壤的冲刷和对养分的淋溶作用。但在图2-c中，河流中TP浓度的月平均变化在枯水期和丰水期之间未呈现出显著差异，推测是在本研究区良好的植被和以水稻为主的种植制度下[20]，多雨季节径流对土壤冲刷造成的河流中磷的增加与稀释作用的效果基本相当。

图2 2003—2016年长乐江出口断面流量(a)、TN浓度(b)、TP浓度(c)各月份实测平均值Fig.2 River flow (a)， TN concentration (b) and TP concentration (c) measured at outlet of Changle River in each month from 2003 to 2016

河流中的养分通量主要取决于2个因素——河流流量与养分浓度。但相关分析表明，河流中TN和TP的负荷量都只与河流的流量呈极显著(P<0.01)相关，而与它们的浓度相关性很弱(图3)。实测河流TN负荷量与河流流量呈极显著线性相关(R2=0.732 9，P<0.01)，而与TN浓度相关性较弱(R2=0.247 5，P>0.05)；实测河流TP负荷量与河流流量也呈极显著线性相关(R2=0.676 1，P<0.01)，而与TP浓度相关性较弱(R2=0.191 4，P>0.05)。可见，虽然河流养分的负荷量来自于河流流量与养分浓度的乘积，但在研究期间长乐江TN和TP的负荷量变化都主要取决于河流流量的变化，养分浓度变化的影响比较弱。本研究显示，河流养分负荷量与流量之间具有显著的相关关系，这与许多前人的研究结果一致[1，21-22]；而污染物负荷量与浓度的相关关系则取决于河流的状态，通常流量较稳定的河流，其河流污染物负荷量与浓度的相关性会高于流量涨落起伏较大的源头地区的河流。

图3 长乐江出口断面实测TN、TP负荷量与河流流量，以及TN、TP浓度的相关分析Fig.3 Correlation of measured TN， TP load with river flow and TN， TP concentration at outlet section of Changle River

2.2 河流氮、磷负荷量估算

为了更好地对河流流量和养分浓度的每日连续性变化进行表征，采用LOADEST模型，先建立河流营养物负荷量与河流流量之间的关系，再利用每日的流量资料来估算河流中每日的营养物负荷量。

首先，把2003—2016年的水质监测数据和水文资料分成奇数年和偶数年2组，基于奇数年各月份的水质监测资料与监测当日的流量资料，通过LOADEST模型来建立和优化长乐江出口断面上的TN、TP负荷量与河流流量之间的回归方程。LOADEST模型运行结果表明，通过AIC和SPPC优选，2003—2016年间TN、TP负荷量与出口断面流量之间的最优模型均为7参数模型(即LOADEST系列模型中的第9个方程)，模型各参数的优化标定结果如表1所示。然后，用偶数年水质监测当日的平均流量来模拟估算所有偶数年水质监测当日的TN和TP负荷量，并将计算结果与实测负荷量进行比较，来验证回归方程模拟结果的准确性和可靠性(图4)。结果表明，在参数标定期(奇数年)，TN拟合方程的决定系数R2=0.89，纳氏模拟效率系数Nse=0.86；验证期(偶数年)模拟验证结果的决定系数R2=0.78，纳氏模拟效率系数Nse=0.76，达到了“非常好”的水平(R2和Nse均大于0.75)[23]。TP的方程拟合和验证效果与TN模型相比稍差，但也达到了“好”(R2和Nse均大于0.60)的水平：标定期(奇数年)模型拟合的决定系数R2=0.76，纳氏模拟效率系数Nse=0.69；验证期(偶数年)模型拟合的决定系数R2=0.67，纳氏模拟效率系数Nse=0.66。

图4 长乐江出口断面TN、TP负荷量的LOADEST模型拟合效果和验证结果Fig.4 Calibrations and validations of LOADEST model for TN and TP load in outlet section of Changle River

表1 长乐江2003—2016奇数年TN和TP负荷量的LOADEST模型率定结果Table 1 Calibration results for LOADEST model of TN and TP loads in Changle River using 2003-2016 odd years’ water quality mornitoring data

验证结果表明，所优化的LOADEST模型可以直接用于所有研究期间长乐江每日TN和TP负荷量的模拟。为了进一步提高模拟的精度和可靠性，采用这14 a的全部实测资料，重新拟合优化LOADEST模型，所得到的TN拟合模型的R2=0.82，Nse=0.80；TP拟合模型的R2=0.76，Nse=0.67。模型具体形式如下：

lnLTN=7.483 4+1.141 8lnQ-0.004 1lnQ2+0.035 1sin(2πT)+0.221 4cos(2πT)-0.004 1T+0.003 4T2；

(1)

lnLTP=4.300 2+1.053 2lnQ-0.048 4lnQ2-0.048 5sin(2πT)+0.056 8cos(2πT)+0.038 9T-0.005 1T2。

(2)

式(1)、(2)中，LTN和LTP分别为通过监测断面的TN和TP的日负荷量(kg·d-1)，Q为日平均流量(m3·s-1)，T为分数形式的日期(如2016年2月5日记为36/366)。基于(1)式和(2)式，应用2003—2016年的每日流量资料，估算长乐江出口断面研究期间的每日TN和TP负荷量，进而获得各个时间尺度(月、年等)的TN、TP通量(图5)。结果表明，TN的年负荷量为989.03～3 578.46 t·a-1，平均年负荷量为(2 207.71±862.48)t·a-1。山区河谷地区河流流量(特别是洪峰流量)变化巨大，导致TN负荷量的日平均值变化很大，为(6.04±11.81)t·d-1。TP的年负荷量为25.80～122.21 t·a-1，平均年负荷量为(71.43±31.56)t·d-1，平均日负荷量为(0.19±0.26)t·d-1。长乐江TN和TP的负荷量在2003—2016年间总体呈上升趋势。从年内变化来看，多年平均的TN月负荷量在12月份最低(1 395.94 t)，而最高值出现在6月份(4 316.65 t)，与流量最大值的出现时间一致；多年平均的TP月负荷量最低值出现在1月份(49.18 t)，最高值与TN一样出现在6月份，为143.53 t，变化趋势与在其他河流上的研究结果类似[21-22]。需要注意的是，在特别高的流量条件下(如洪峰期间)，TN、TP负荷量估算可能存在低估的现象[1]。

图5 2003—2016年长乐江出口断面TN和TP年负荷总量和平均月负荷量分布Fig.5 Annual loadsand mean monthly loads of TN and TP at toutlet section of Changle River from 2003 to 2016

2.3 河流氮磷负荷量变化小波分析

采用Morlet小波对2003—2016年长乐江出口断面的月TN和TP负荷进行连续小波变换，获得它们在研究期间的小波系数图谱(图6)。图6-a展示了长乐江出口断面以月为单位的TN负荷量小波系数变化(振荡)，大约以12～24个月的时间跨度(时间尺度)为周期，发生明显的高低循环变化；特别是自2009年起变化加剧，反映了TN负荷量振荡强度加剧，振荡强度的正负值差加大。此外，在时间跨度为3～12个月的尺度上，也有一些较明显的周期性变化。TP负荷量小波系数的变化周期与TN负荷量的变化类似，但变化强度相对较小，并且变化强度在研究期间一直呈增强趋势，最大振荡的周期出现在2014年以后(图6-b)。

分析研究期间监测断面上的河流流量(图6-c)、TN浓度(图6-d)和TP浓度(图6-e)的小波系数变化频谱，可以发现，河流流量的小波系数变化与TN负荷量小波系数的变化最为相似，包括高低值振荡的时间跨度(时间尺度)和振荡发生的日期等；TN浓度小波系数的变化在时间尺度上与TN负荷量的变化类似，但变化强度大约仅为河流流量变化的1/200。TP浓度小波系数的变化，在振荡的主要时间尺度和振荡强度上均与TP负荷量小波系数的变化图谱相差较大。因此，虽然TN和TP的负荷量是河流流量与TN浓度的乘积，但小波系数频谱分析表明，由于流量的变化远大于TN和TP浓度的变化，所以长乐江TN和TP负荷量的变化主要取决于流量的变化，而TN和TP浓度变化的影响相对较小。由于河流非点源污染过程主要受水文循环驱动，因此尽管气候和物候条件并不相同，不同地区河流氮、磷负荷量的年度变化周期性却既有类似的规律，也有区域性的特征。例如在松花江哈尔滨段上的类似分析表明，污染物负荷量在1 a(12个月)和2 a(24个月)时间尺度上的2个变化主周期基本与本研究一致。这是因为，无论是在浙江还是在黑龙江，河流污染物负荷量变化都与地表径流过程的变化基本一致；但同时它们也有各自其他时间变化尺度上的特点[24]。

黄色和绿色表示接近小波系数的平均值，红色和蓝色代表小波系数偏高和偏低的状态，线条为偏高(正值)或偏低(负值)相对量的等值线。Yellow and green represented data close to the average of wavelet coefficients. Red and blue represented high and low wavelet coefficients. The line represented the isoline of relative quantity of high (positive) or low (negative).图6 长乐江出口断面2003—2016年TN(a)和TP(b)负荷量、河流流量(c)、TN(d)和TP(e)浓度小波系数频谱等值线图Fig.6 Contour of wavelet coefficient for TN (a) and TP (b) loads， river flow (c)， TN (d) and TP (e) concentrations at outlet section of Changle River from 2003 to 2016

利用小波方差可以更明确地定量判别不同时间尺度下所发生的TN和TP负荷量振荡的相对强弱和各种振荡发生的时间周期。如图7所示，TN和TP的小波方差的峰值分别出现在6、12和18个月的时间尺度上(图7-a和7-b)，意味着它们负荷量从低值到高值的轮回变化平均分别有6、12和18个月3个周期，其中18个月时间尺度上的小波方差最大，是变化最明显的主周期。TN和TP负荷的这种变化规律首先是由河流流量的变化周期决定的(图7-c)，而TN浓度(主周期为19个月)和TP浓度(24个月)的变化影响较小，因为它们的变化强度远远小于流量的变化强度(图6)。

图7 长乐江出口断面2003—2016年TN(a)和TP(b)负荷量、流量(c)、TN(d)和TP(e)浓度的小波方差Fig.7 Wavelet variance of TN (a) and TP (b) loads， flow (c)， TN (d) and TP (E) concentrations at outlet section of Changle River from 2003 to 2016

在一些大型流域，地下径流的水龄可长达数十年[25]，而通过地下径流产生的河流非点源污染也具有数月至数年的滞后作用[26]。因此，长乐江TN、TP负荷量的变化，除了受河川径流的影响而形成以18个月为主周期的重要变化规律以外，还可能与当地农业种植制度、TN和TP的投/排放，以及上游水库的蓄/放水操作等有关，形成了其他次要的多尺度周期性变化。这些多尺度的变化机制更为复杂，尚待进一步研究。

3 结论

(1)2003—2016年间长乐江出口断面水质监测表明，该河流水质明显受到氮、磷的污染，水体中的TN和TP浓度平均值分别为(3.74±0.91)mg·L-1和(0.136±0.086)mg·L-1，均高于该河流TN和TP浓度的参照标准。河流中TN和TP的实测负荷量都与河流的流量呈极显著(P<0.01)相关，决定系数R2分别为0.732 9和0.676 1；但与TN和TP浓度的相关性较弱，决定系数R2分别为0.247 5和0.191 4。

(2)低频率的人工采样水质监测不能反映河流水体中负荷量的每日连续变化。利用LOADEST模型和2003—2016年间每月1次的水质监测资料，建立了长乐江出口断面上TN和TP的负荷量估算模型。模型验证结果较好，其中TN负荷量的模拟结果与实测值对比的决定系数和纳氏模拟效率系数分别为0.78和0.76；TP负荷量的模拟结果与实测值对比的决定系数和纳氏模拟效率系数分别为0.67和0.66。

(3)利用LOADEST模型核定了长乐江TN和TP的每日负荷量，进而获得各个不同时间尺度上TN和TP的负荷量。结果表明，在2003—2016年间，长乐江出口断面的TN平均年负荷量为2 207.71 t·a-1，TP平均年负荷量为71.43 t·a-1；它们在研究期间均总体呈上升趋势。多年平均的TN月负荷量在12月份最低(1 395.94 t)，在6月份最高(4 316.65 t)。

(4)2003—2016年以月为单位的TN和TP负荷量小波系数分析结果表明，长乐江TN负荷量以18个月的时间跨度为主周期发生着周期性的变化；TP负荷量的变化主周期与TN负荷量的变化类似，但变化强度相对较小。TN和TP负荷量变化的小波方差分析还发现，长乐江TN和TP负荷量的变化主要取决于河流流量的变化，而受TN和TP浓度变化的影响相对较小，因为TN和TP浓度的变化强度远小于流量的变化强度。