不同插值方法在面源污染入河量计算中的适用性研究

章 晶，张新华，蔡育杰，金 涛，唐 玮

（四川大学水力学与山区河流开发保护国家重点实验室，四川 成都 610065）

随着全球经济发展，面源污染已经成为水污染的主要原因[1].其中污染入河量计算不仅有识别水体污染主要来源的作用，更能对污染源的特征进行明确分析，是一项必要的基础性工作[2]. 但面源污染所具有的随机性、广泛性和滞后性等特征使得在监测与计算中的难度较大[3]，而空间插值因其能够很好反映区域变量的空间真实变化等优点，使得其在面源污染计算的应用已经越来越广.目前在面源污染负荷计算应用最广泛的方法是1996 年Johnes[4]提出的将土地利用类型、畜禽养殖、农村生活排放等不同面源类型综合考虑的输出系数法. 蔡明等[5](2004)通过考虑流域损失改进了输出系数模型，探究了面源污染入河量的情况；胡昱欣等[6](2015)利用克里金插值方法对东辽河30 个子流域出口的入河系数进行空间分析，得到东辽河流域泥沙入河量总体较大，而面源入河量呈现东南小西北大的特征；胡富昶等[7](2019)采用反距离权重法对降雨影响因子进行空间插值，分析射洪县污染空间的分布状况.

然而不同的插值方法因其原理不同和对样本点的要求不一致，所反映的空间效果存在一定差异，这对面源污染入河量的计算存在很大影响.本文基于输出系数法对研究区3 种污染源进行负荷总量计算，分析该区2018 年TN 负荷总量空间分布特征；此外，针对反距离权重法受样本点密度的影响较大的特征，对该方法进行均匀增加样本的修正手段，并进一步对比普通克里金法和修正前后的反距离权重法在计算污染入河量的结果差异，利用渠江出口实际负荷量检验插值结果的精确性，探究不同插值方法在计算面源污染物入河量的适用性.

1 研究区域概况

渠江发源于川陕交界处米仓山系铁船山，是长江支流嘉陵江左岸最大支流. 自东北向西南经过平昌县、渠县、广安等县境，在重庆市合川区汇入嘉陵江.由于处在亚热带湿润季风区，受季风影响，降水充足，多年平均降水量在1 000 ～1 600 mm 之间.径流主要由降雨补给，5 ～10 月为汛期，11 月至次年4 月为枯季，径流年内分配极不均匀，汛期径流量占年径流量的86.1%，而枯期仅占13.9%[8].渠江流域作为川东北粮仓，全流域内耕地占比55.47%，特别是干流两岸及支流州河流域，农业生产水平较高，受面源污染影响严重.

研究区域为渠江流域主要涉及到的19 个区县，如图1 所示.涉及总人口1 634.24 万，其中农业人口有951.69 万人.本文以III 类水质为标准，依据渠江出口在2015 ～2018 年的监测断面资料可知，全流域主要污染物为总氮TN，其中2016 ～2018 年渠江出口处TN 浓度全面超标，2018 年渠江出口TN 的平均浓度为1.53 mg/L，根据地表水环境质量标准为V 类水质.

图1 渠江地理位置图Fig.1 The geographical position of Qujiang River

2 研究方法

2.1 输出系数法

输出系数法通过采用输出系数估算各类农业污染源(土地利用、畜禽养殖、农村生活)排污量，计算得到各区县污染物负荷总量.输出系数模型为[9]:

式中:Sj为污染物j在区域的总负荷量，t/a；j为污染物类型；Eij为污染物j在区域第i种土地利用类型的输出系数或牲畜、人口的输出系数，t/(km2·a)；i为区域土地利用类型的种类或牲畜、人口数量；Ai为第i种土地利用类型的面积(km2)，或牲畜数量(头)，或人口数量(人).

基于各区县污染物负荷总量，采用入河系数计算各区县面源污染物入水体负荷量. 入河量表达式为[10]:

式中:Lj为污染物j进入河流水体的负荷量，t/a；λ为入河系数.

2.1.1 输出系数的确定

本文根据国内相似区域已有的研究结果并结合研究区的实际情况[7，11-14]，确定出不同土地利用、畜禽以及农业生活的输出系数取值，如表1 所示.

表1 渠江流域面源污染物输出系数表Table 1 Qujiang river basin non-point source pollutant output coefficient table

2.1.2 入河系数的确定

为了更好反映区域污染物的实际入河程度，通常以入河系数作为重要参数，利用空间插值方法得到全区域污染物的入河情况[15]. 本文利用两种插值方法进行空间分析，获取全区域TN 的入河情况，其中采用的入河系数点来自对渠江流域构建的SWAT 模型，根据TN 模拟结果，计算出渠江流域50 个子流域的入河系数，得到入河系数的范围在0.005 ～0.99 之间.

2.2 普通克里金法(OK)

普通克里金法是随着地质统计学的发展而衍生出的一种无偏估计的空间插值方法.其原理是利用已知的插值样本加权平均值估计其余待预测点值，该方法能使预测值等于实际值的数学期望值并保持方差最小，能保证局部误差最小等优点，但也会产生平滑效应，即较小值易被夸大，较大值易被低估.其计算公式如下[16]:

式中:Z为预测值；λ表示克里金法权重系数；Z(Xi)表示实测点Xi处的值.

2.3 反距离权重法(IDW)

反距离权重插值是基于相近相似原理，其内涵为离得越近的样本之间有着相似的性质；相反，样本之间离得越远其相似性就越小.具体计算公式如下[16]:

式中:Z是预测值，n为总样本数，Zi为第i个样本点值，di为预测点与样本点的分离距离，p为加权幂.

这种方法在计算中相对容易，但在入河系数的应用上有一定局限性:1)与地质统计学方法相比，IDW无法在无样本处估计预测值方差；2)当样本点空间分布不均匀时，IDW 可能存在问题；3)IDW 的计算精度受样本点密度的影响较大[17].

针对以上问题，本文在利用IDW 对入河系数进行空间插值时，因样本点分布不均匀，在样点未覆盖到的岳池县、仪陇县和大竹县区域根据相近相似原则增加了样本点，并且对污染入河量较大的渠县区域的样本点密度进行了加密，并进行了修正前后的对比.

3 结果与分析

3.1 面源污染物负荷总量空间分析

通过输出系数法计算得到研究区2018 年面源污染负荷总量Sj(表2)，由表2 可知研究区产生的TN负荷量为67 512.46 t/a，土地利用产生的污染负荷量最大，其占总负荷量的84.4%；其中耕地对土地利用产生的污染负荷总量的贡献为51 191.22 t/a，占整个土地利用产生的75.8%. 从行政单元来看，巴中市和达州市由于耕地面积较大，且农村人口较多，故产生的面源污染负荷总量也较大.

表2 研究区面源污染负荷总量计算结果Table 2 Calculation results of total non-point source pollution load in the study area

通过Arcgis 软件对TN 负荷量进行空间分析，得到研究区2018 年各区县TN 负荷量分布图(图2)和各区县3 类污染源贡献情况(图3). 由图2 看出，通江县和宣汉县产生的TN 负荷量最大，而通川区、开江县、前锋区和华蓥市产生的TN 负荷量较小；由图3 可以看出，对每个区县来说土地利用产生的污染负荷量比重均最大，农业生活产生污染量次之，畜禽养殖污染产生量最小.可以推断是因为研究区人口密度大且农业人口占一半以上，农民以农业生产为主，加上耕地面积占比大，因而土地利用产生的污染负荷量偏大，另外由于农村居民分布不集中，集水管网和污水处理设施的建设远远不够，导致生活污染物以及生活废水任意排放.

图2 渠江各区县TN 负荷总量空间分布Fig.2 Spatial distribution of total TN load in each district and county of the Qujiang River

图3 渠江各区县不同污染源占比情况Fig.3 The proportion of different pollution sources in each district and county of the Qujiang River

3.2 不同插值方法下TN 入河系数与入河量空间分布

1)不同插值方法对TN 入河系数的空间分布影响

由图4 可以看出，两种插值方法一定程度上均能反映研究区入河系数的空间分布特征，整体上大致呈东北部向西南逐渐增加后减少的趋势.研究区东北部位于渠江源头处，植被覆盖面积大，对污染物入河有一定程度的阻挡作用；又因上游巴河、州河在渠县开始汇入渠江干流，污染物入河系数也逐渐增加至渠县境内达到最大；最后随着地势变缓，研究区西南部的污染物入河系数渐渐减小.

图4 不同插值方法下的入河系数空间分布图Fig.4 Spatial distribution of river entry coefficients under different interpolation methods

具体来看， IDW 法在修正前的插值结果整体小于其他两个，且呈现出入河系数大的位点分布不均的状态，这是由于插值点的不均匀性导致的. 而OK 法和修正后的IDW 法的插值结果趋势上差距不大，但OK 法相比IDW 法存在平滑效应，使得总体插值范围在0.134 ～0.849 之间，与样本点实际的入河系数范围0.005 ～0.99 相比缩小了许多. 另外，IDW 法受样本距离远近的影响，对上游南江县、万源市、宣汉县的入河系数插值结果偏小. 综上，由于插值方法原理的不同和插值样本的处理方式不一致，两种插值结果在各个区县的空间分布上存在差异.

2)不同插值方法对TN 入河量的空间分布影响

不同插值方法对TN 入河量的空间分布如图5 所示，研究区的TN 入河量在空间上分布不均，但总体趋势为中间高并依次向两边递减. 其中通江县、宣汉县和平昌县虽然位于研究区的渠江上游处，但因其较其他区县的耕地面积大，导致入河量也大；渠县的TN 入河量贡献为研究区最大，这是因为渠县的入河系数大小为上游城市的两倍左右；而广安区虽然入河系数大，但其地域面积比渠县小，导致TN 负荷总量小，入河量也减小.

图5 渠江各区县TN 入河量空间分布Fig.5 Qujiang River spatial distribution of TN intake in each district and county

总体来看，OK 法计算出的研究区TN 入河量在空间上更能反映各区县污染的实际入河情况，而IDW法因在南江县、仪陇县境内的插值点较少，使得这两个区县的入河量偏小. 其中IDW 法在修正前的计算结果整体上小于其余两个计算结果，这是由于整个研究区的入河系数在空间插值后的结果偏小导致；而修正后的IDW 法与OK 法相比，入河量结果稍微偏大，但相差大为减小，满足误差要求.

3.3 结果验证

本文采用两种插值方法对入河系数进行空间分析，并对IDW 方法修正前后的结果也进行了对比，结合输出系数法计算出2018 年研究区TN 入河量，并选取渠江出口断面的实际负荷量来对三种计算结果进行验证.结果见表3，可以看出OK 法在对TN 入河量计算上有很好的效果，其相对误差仅为0.73%；而修正后的IDW 法相对于修正前相对误差降低了28.59%.由此表明IDW 法采用的相近相似原理在对TN 入河量计算上有一定的偏差，在使用时需要对插值点进行修正，修正后的插值结果可以达到较好精度.

表3 结合不同插值方法的ECM 模拟结果对比Table 3 Comparison of ECM simulation results combining different interpolation methods

4 讨论

从污染源角度分析，土地利用产生的污染负荷量达到总量的84.4%，这与研究区以农业生产为主的情况相符合；从污染入河情况分析，研究区TN 平均入河系数为0.39，其中渠县的污染流失程度最大，应当受到重视；从插值方法分析，不同插值方法计算出的污染入河系数空间分布略有不同，OK 法的插值结果不可避免存在平滑效应，即夸大较小值，低估较大值，但OK 法可以保证局部误差达到最小，其在面源污染入河计算结果上也相对最优，而IDW 法的计算精度受插值点数量和密度影响较大，要得到良好的插值结果对插值点的空间分布情况要求会更高.

5 结论

通过输出系数法对研究区2018 年TN 负荷总量情况进行空间分析，并采用OK 法和修正前后的IDW法对入河系数进行空间插值，探究了不同插值方法对TN 入河量的影响，得到结论如下:

1)在对研究区的TN 入河情况进行空间分析时，两种插值法均能反映研究区TN 入河情况. 采用OK法能够更全面真实的反映各个区县的污染入河情况，而IDW 法在使用时要注意对插值点的空间分布形态、均匀程度等进行分析，对分布不均，密度不够应适当修定.本文采用加密方法进行适当修正，修正后的IDW 法和OK 法在入河量计算上都较为精确.

2)从污染源角度来看，研究区农业人口数量较大，耕地面积广，在耕种过程中的各类化肥以及农药的使用造成耕地的TN 负荷量高，其中耕地产生的TN负荷量占总量的75.8%；从污染流失程度来看，渠县作为污染流失最严重的地区，面临的面源污染风险最大，是今后污染治理的首要对象.其次，要大力普及测土配方施肥，推广病虫害绿色防控技术，着力提升农业科技水平，方能达到有效控制面源污染的目的.