伴随人类活动加剧和流域经济发展，大量污染物随河流输送进入海洋，导致近岸海域富营养化、生态健康退化、生态系统服务功能下降。河流污染物输移通量是指单位时间内某污染物通过某断面的总量^[1-2]。污染物入海通量估算作为研究断面河流污染负荷，描述流域污染物产生和输送的主要方法，是入海污染物总量控制的关键，对近岸海域的污染防治及水生态环境保护决策具有重要意义。

在污染物通量估算中，河流水质监测数据源主要包括自动监测和人工监测两类，二者各有优缺点^[3-5]。基于自动监测的水质数据能以4 h 1次的监测频率及时、准确地记录水质变化；但受设备故障等技术问题和天气等客观条件影响，自动监测水质数据会出现短则缺失几个小时，长则中断数月的数据缺失情况，监测数据的时间覆盖范围有限，无法获取对应时间段内的污染物通量^[6-8]。基于人工监测的水质数据可提供时间跨度较为完整的月度水质信息，但每月1次的低频数据整体水质特征代表性不强，对应时段内河流污染物通量估算结果准确度不高^[9-10]。国内外学者对于河流污染物通量估算方法进行了多方面探索，无论是以通量定义为出发点进行统计分析^[11-14]，还是依托数值模型完成通量模拟及校核^[15-17]，现有的河流污染物通量估算研究大多固定一种监测数据来源，特别是单一水质数据源，鲜有考虑如何在单一源数据受限的情况下尽可能准确地获取污染物通量。而研究发现，采用数据融合方法或整合多个数据源的数据能够有效提高污染物浓度的准确度^[18-19]，弥补单一数据源在进行污染物通量估算时的不足，如数据质量参差不一，数据缺少一定的冗余和互补等，可用于获得更为准确的污染物通量。

因此，本研究以山东省小清河为示范河流，依托2018—2020年的流量与氮、磷浓度监测数据分别估算基于自动监测和人工监测水质数据的氮磷污染物入海通量，并在此基础上整合可利用的月通量结果，实现互为补充，相互补强，最终获得多源数据下的小清河氮磷污染物入海年通量及估算误差，以期为准确估算河流污染物入海通量提供方法和流程实证，为小清河流域污染物总量控制和水环境治理提供数据参考。

1 材料与方法 1.1 研究区概况

小清河发源于济南西郊睦里闸，干流全长229 km，流经济南、淄博、滨州、东营和潍坊5市12县(市、区)，最终于潍坊寿光汇入莱州湾(见图 1)。作为流入莱州湾西南岸的主要河流之一，小清河沿岸经由山东半岛蓝色经济区、黄河三角洲高效生态经济区和省会城市群经济圈，在山东省经济社会发展中扮演着重要角色。然而，城市化进程的加快让小清河受到了严重的人为污染，大量工业废水和生活污水排放其中，严重威胁河口及莱州湾的生态环境质量状况。

1.2 数据来源及预处理

为最大可能保证入海通量结果的准确性，依据就近原则选择流域中处于相同或相近位置的水文站与水质站进行通量计算。本文选取范李站监测结果代表小清河的入海水质数据，监测指标为总氮(TN)、总磷(TP)；该站点临近的水文站点为石村站。其中，人工监测水质数据为2018—2020年的每月监测结果，自动监测水质数据的时间跨度为2018年7月至2018年12月、2019年6月至2020年6月，监测频次为每4 h 1次。以上水质数据均来源于山东省生态环境厅。2018—2020年的逐日流量数据来源于山东省及各地市水文局。

受运行维护或仪器故障等客观因素影响，自动监测水质数据存在不同时间跨度上的缺失，因此需要对该部分氮磷浓度数据进行填充处理。整合2018—2020三年间TN、TP自动监测数据缺失情况如图 2所示。其中，深色部分表示该时间记录点内记有水质数据，白色部分则表示该时间记录点对应的水质数据缺失。由图 2可知，氮磷浓度数据的缺失情况有所差异，可被视作遵循随机缺失机制^[20]。

对于自动监测数据中的缺失部分，平均数填充、中位数填充及线性插值法是短时序缺失数据插补的常用方法^[21]，人为掩盖完整监测数据以模拟短时间序列缺失情景后调用上述三种方法进行插补。借助平均绝对百分比误差(Mean absolute percentage error, MAPE)评估插补效果，MAPE越小，填充值和原始值之间的差距越小，方法填充效果越好。填充效果评估指标MAPE计算公式如下：

$ M A P E=\frac{1}{n} \sum\limits_{i=1}^n\left|\frac{M_i-O_i}{O_i}\right| \times 100 \%。$

(1)

式中：M_i为模拟缺失情景下的填充值；O_i为对应的原始完整监测数据值。

以MAPE最小的最优填充法填充自动监测水质数据中的短时序缺失部分。对于自动监测水质数据中的长时序缺失情景，类似地，选取线性插值法、K最近邻插补(K-Nearest neighbor, KNN)及链式方程的多重插补(Multiple imputation by chained equations, MICE)三种方法进行最优填充法评估^[22-23]，并用最优填充方法填充原始自动监测水质数据中的长时序缺失部分，以此迭代式地获得完整自动监测氮磷浓度数据。

1.3 污染物通量估算方法 1.3.1 Monte Carlo法

Monte Carlo法是基于概率估计不确定事件的一类方法^[24]，可通过使用大量的随机数组执行模拟得到一个问题的最优解，核心是随机抽样。在本文的Monta Carlo模拟中，借助随机抽样模拟每日氮磷浓度数据和流量数据的可能取值，并利用通量公式计算出对应取值下的河流污染物入海通量，不断重复该模拟过程，最大限度地获得河流污染物入海日通量的可能取值，累加后便可得到河流污染物入海月通量的估算结果和波动范围。

对于自动监测水质数据源，计算求得氮磷浓度的日均值。单日情况下，利用氮磷浓度日均值和日径流量根据公式(2)计算当日氮磷污染物通量，累加日通量结果后可得到基于自动监测水质数据的当月氮磷污染物通量估算值。

$ W_{\mathrm{d}}=0.086\;4 \times C_{\mathrm{d}} \times Q_{\mathrm{d} }。$

(2)

式中：W_d为当日污染物通量，t；C_d为当日氮磷浓度均值，mg/L；Q_d为日径流量，m³/s。

结合数据填充处理中相应MAPE结果可获得基于日均值的氮磷浓度波动范围。流量数据参照《河流流量测验规范》(GB50179—2015)中的流量断面单次测验允许误差(5%)，可获得基于日径流量的流量数据波动范围。借助Monte Carlo法分别在水质和流量数据范围内进行随机抽样10 000次，根据公式(2)计算当日污染物通量值，共计10 000组。累加当月10 000组污染物日通量结果可得到10 000个污染物月通量取值，以月通量计算结果的最小值和最大值为下限值和上限值，得到当月污染物通量的估算范围。

1.3.2 LOADEST模型

LOADEST模型利用给定的离散污染物浓度数据和对应的流量数据构建污染物通量与流量的回归方程，根据方程估算不同时间尺度下的污染物通量，同时给出估算通量的95%置信区间，增加结果的科学性^[25]。目前，该模型已被应用于国内外多个流域，是进行低频水质监测河流污染物通量的有效途径之一^[26-28]。该模型提供了3种参数拟合估算方法，分别是最大似然估计法(MLE)、渐进极大似然估计法(AMLE)和最小绝对偏差法(LAD)^[29-31]。对于残差符合正态分布的删失型数据，采用AMLE法；对于非删失型数据，MLE和AMLE均可使用；当数据残差不符合正态分布时，无论是否存在删失，均应采用LAD法。

模型运行时，根据基于最大熵原理的赤池信息准则(Akaike information criterion, AIC)和Schwarz后验概率准则(Schwarz posterior probability criteria, SPPC)，从备选的9个污染物通量回归方程中选取最小AIC值和最小SPPC值的方程为最优方程。对于通量回归方程拟合效果的检验，可借助判定残差序列相关系数(SCR)检验残差是否存在序列相关性。SCR值越小，残差间越独立。同时通过相关系数R²判定方程的总体拟合优度，R²越接近1，表示方程拟合效果越好。对于非删失型数据，运用概率曲线相关系数(PPCC)检验方程残差是否符合正态分布，PPCC越接近1，回归方程的残差越符合正态分布。

1.3.3 月单值法

月单值法利用每月一次的人工监测水质数据和日径流量数据借助算术方法估算污染物通量，具体计算公式如式(3)所示。月单值法误差较大，一般用于大致估计污染物通量及描述通量整体变化趋势。

$ W_{\mathrm{m}}=0.086\;4 \times C_{\mathrm{m}} \times Q_{\mathrm{m}} $

(3)

式中：W_m为当月污染物通量，t；C_m为月度污染物浓度值，mg/L；Q_m为月径流量，m³/s。

1.4 误差分析

污染物通量估算误差e的计算公式如公式(4)所示：

$ e=\left(\frac{W_{\max }-W_{\min }}{2 W}\right) \times 100 \%。$

(4)

式中：W为利用水质和流量数据直接得到的污染物通量估算值，t；W_max为污染物通量波动范围中的上限值，t；W_min为污染物通量波动范围中的下限值，t。

2 结果与讨论 2.1 基于自动监测水质数据的氮磷污染物入海通量估算

针对短时序缺失数据，TN、TP的最优插补法均为线性插值法，MAPE分别为5%和9.8%，小于10%这一《水污染源在线监测系统验收技术规范》(HJ/T354—2007)规定的质控样测定时的最大相对误差。因而，该部分自动监测数据的填充值可近似等同于原始值，数据波动的最值和日均值之间的相对误差也同原始完整数据保持一致，统一设为±10%。而对于长时序缺失情况，TN的最优插补法为MICE，MAPE为7.54%；TP的最优插补法为线性插值法，MAPE为17.76%。对比质控样测定的最大相对误差10%，该部分自动监测填充数据的波动范围分别设为±10%(TN)和±17.76%(TP)。

运用Monte Carlo法计算小清河氮磷污染物的入海通量(见表 1)。结果表明：基于自动监测水质数据估算得到的各月份小清河氮磷污染物入海通量整体分布较为集中，污染物通量估算误差大致呈5%分布，其中TN为4.77%，TP为5.37%。

2.2 基于人工监测水质数据的氮磷污染物入海通量估算

根据LOADEST模型，优选出2018—2020年间可用于估算氮磷污染物入海通量的最佳回归方程如表 2所示。其中，Load为污染物通量；Q为监测断面的流量；dtime为模型根据监测数据的频率复核得到的时间累积。

以2019年为例，模型各参数的优化标定结果如表 3所示。TN的通量回归方程检验结果为R²=0.990 3，PPCC=0.978 4，SCR=-0.107 6；TP的通量回归方程检验结果为R²=0.896 0，PPCC=0.987 7，SCR=-0.323 2。检验结果表明：TN和TP的通量回归方程拟合效果较好，残差不存在序列相关性，且服从正态分布。因而，基于LOADEST模型优选出的污染物通量回归方程能够用于该站人工监测水质数据下的氮磷污染物入海通量估算，通量结果见表 4。其中上下限值指95%置信区间。

由模型95%置信区间内的通量估算值、上限值、下限值计算月通量估算误差，得到TN的e为30.87%、TP的e为34.66%。对比根据自动监测水质数据估算得到的污染物通量，2018—2020三年间小清河基于人工监测水质数据的氮磷污染物入海月通量估算结果区间分布相对较广。结果表明：尽管LOADEST模型的模拟结果可靠，但针对月通量模拟值的通量浮动范围较宽，模拟准确度不高，这也印证了该模型本身的不足之处^[1]。

2.3 基于数据融合的小清河氮磷污染物入海通量估算

为进一步探讨不同数据源污染物入海通量结果间的可比性，引入月单值法结果对上述小清河氮磷污染物入海通量的估算值进行综合分析。以2020年为例，将相同月份的自动监测水质数据污染物通量估算值、人工监测水质数据污染物通量估算值及月单值结果进行比较(见图 3)。结果显示：小清河各月份基于自动监测水质数据的污染物通量估算结果、基于人工监测水质数据的污染物通量估算结果和月单值法的计算值三者间数据相差不大，大部分情况下的通量值较为一致，变动趋势也基本同步。

综合上述分析，对于存在自动监测水质数据记录的月份，以基于自动监测水质数据计算得到的污染物通量结果及其上下限值为该月份的通量取值，其余月份选取基于人工监测水质数据的LOADES模型模拟结果及其限值范围。整合数据，累加各月通量估算值可得当年污染物年通量估算值。而为获取入海年通量的波动情况，再次借助Monte Carlo法对2018—2020年TN、TP的月通量范围进行抽样，重复10 000次后累加得到当年污染物年通量的上下限值，汇总数据于表 5。结果表明：研究期间小清河TN入海输送量整体呈现先上升后下降趋势，TP入海通量总体呈现下降趋势，且TN、TP入海年通量的整体估算误差均低于25%，实现了基于多源数据的河流污染物入海通量估算。

3 结论

基于河流监测断面的水量数据、自动监测和人工监测氮磷浓度数据，本文提出一种基于多源数据的河流污染物入海通量估算方法，用以补足单一数据源，特别是单一水质数据源无法计算完整河流污染物入海通量或通量估算准确度不高的现实情况，并利用示范河流小清河2018—2020年的监测数据进行分析与研究。主要结论如下：

(1) 结合小清河水文水质特征，对高频自动监测水质数据进行预处理后利用Monte Carlo法估算得到2018—2020年小清河TN和TP入海月通量值及估算误差。结果表明：基于自动监测水质数据得到的各月份小清河氮磷污染物入海通量整体分布较为集中，上下限值围绕估算值5%左右分布，其中TN为4.77%，TP为5.37%。

(2) 利用LOADEST模型建立TN、TP与河流流量的回归方程，总体拟合效果较好。借助模型估算得到2018—2020年间小清河入海月通量值及其估算误差。结果显示：基于人工监测水质数据的氮磷污染物入海月通量估算结果相对分布较广，TN估算误差为30.87%、TP估算误差为34.66%。

(3) 综合自动监测和人工监测水质数据，借助Monte Carlo法计算得到多源数据下的2018—2020年小清河TN、TP入海年通量及变化范围，估算误差整体在25%以内，为准确估算河流氮磷污染物入海通量提供方法借鉴，为小清河氮磷污染物总量控制及污染负荷优化提供数据支持。