基于反距离加权插值法评价海域水质类别空间分布

2019-02-20 06:23李翔宇包艳英李玉璞
中国环境监测 2019年6期
关键词:插值法栅格点位

李翔宇,李 曌,包艳英,李 利,李玉璞

1.大连市环境监测中心,辽宁 大连 116023 2.中国环境监测总站,国家环境保护环境监测质量控制重点实验室,北京 100012 3.海军大连舰艇学院,信息系统系,辽宁 大连 116018

《水污染防治行动计划》[1]中规定了近岸海域水质优良(一、二类)比例目标,在计算海水优良比例时,一般按照点位法或面积法进行计算。点位比例计算简单直接,责任明确,但评价结果不连续,与环境要素的管理目标不相吻合。面积比例计算一般采用代表面积法或者插值面积法,代表面积法(即通过核算,给每个监测点位赋以固定的面积)与污染物实际扩散规律不相适应,相邻点位评价结果不连续[2];插值面积法可以反映近岸海域水质空间分布情况,评价结果相对连续,可以作为点位比例法的有效补充。插值方法的选择十分重要,目前常用的空间插值方法有反距离加权插值法(IDW)、克里金插值法(Kringing)、样条函数法(Spline)和趋势面法(Trend Surface)等[3-6]。多项研究表明[7-9],IDW插值方法在海水污染物分布空间表达方面有一定优势。付瑞全等[7]基于ArcGIS软件空间插值模块采用IDW插值法实现了海水污染面积计算,但并未给出具体的插值参数建议;李俊龙等[8]分别采用IDW法、Spline法和克里格法等3种方法13种模型进行了全国近岸海域海水水质模拟,并通过交叉验证得出IDW(p=2)最适宜模拟全国近岸海域(299个点位)水质分布的结论,提出不同区域、不同要素有待进一步研究;李冕等[9]对IDW插值模型进行改进,给出一种自适应凸包选点的插值方法,认为该法能较好地反映海水水质质量趋势。对于空间插值来说,没有绝对最佳的插值方法[3],要针对研究区的情况,对数据样本进行反复试验,得到理想的插值效果。

笔者重点研究了传统IDW插值具体参数(幂参数、参与插值的周围点的数量、插值网格大小)的选取,比较了IDW方法在点位密度不同的开放海域和近封闭内海海域的适用性,并将IDW插值法水质类别综合评价结果与常用的点位法和代表面积法评价结果进行比较分析,指出了方法应用的注意事项及局限性,为有类似评价需求的用户提供参考和借鉴。

1 研究区与数据

以全国近岸海域和渤海海域为研究区(分别代表开放海域与近封闭内海海域)。全国近岸海域布设413个国控监测点位(图1),平均点位密度约为13.8个/万km2(统计数据不包含港、澳、台地区,下同);渤海海域布设163个监测点位(图2),平均点位密度约为21.2个/万km2。2个研究区监测点位分布特点均为临岸较密,离岸较疏,渤海海域的监测点位比全国近岸海域的监测点位更为密集。

海水水质数据来源于全国近岸海域环境监测网某水期监测数据,共有29个监测项目:pH、溶解氧、化学需氧量、生化需氧量、无机氮、非离子氨、活性磷酸盐、汞、铅、镉、六价铬、总铬、砷、铜、锌、硒、镍、氰化物、硫化物、挥发性酚、石油类、六六六、滴滴涕、马拉硫磷、甲基对硫磷、苯并(a)芘、阴离子表面活性剂、粪大肠菌群和大肠菌群等。

图1 全国近岸海域评价区域监测点位分布Fig.1 Assessment range of China coastal area and the distribution of monitoring points

图2 渤海海域监测点位分布Fig.2 Assessment range of the Bohai Sea area and the distribution of monitoring points

2 研究方法

2.1 IDW空间插值

IDW插值算法是基于地理学第一定律相近相似的原理,认为每个采样点都对邻近区域的插值点有一定的影响,且影响的大小随距离的增大而减小。一般公式为

(1)

(2)

(3)

IDW法属于精确插值方法,即插值后表面通过采样点,该特征可以保证插值结果在采样点位处与实际监测值一致。

从公式可以看出,在样本点位固定的情况下,对一个点进行IDW插值时涉及的参数主要是幂参数(p)和参与插值的采样点数量(N)。对一个区域进行插值并展示时,还应涉及插值网格的大小。

2.2 插值法水质类别综合评价

目前中国采用的海水点位水质类别评价方法为单因子指数法[10]。插值法海域水质类别综合评价与单因子指数法理论一致,不同的是将“点”扩展到“面”。步骤如下:①单指标空间插值,即针对每项监测指标,对所有监测点位的监测结果进行空间插值,获得每项指标浓度的连续空间分布栅格数据集;②单指标水质分类,依据《海水水质标准》(GB 3097—1997)[11]对单指标插值结果进行水质类别分类赋值,以1、2、3、4、5分别代表一类、二类、三类、四类、劣四类海水,获得每项指标水质类别的空间分布栅格数据集;③多指标水质类别判定,通过栅格数据叠加分析,取某处最差水质类别(多层栅格数据同一位置的最大值)作为该处最终水质类别,获得评价区域内的海水水质类别空间分布情况;④通过统计各水质类别栅格数量计算其面积及所占比例,输出近岸海域水质类别空间分布图。

2.3 方法的实现

单指标空间插值、单指标水质分类和多指标水质类别判定采用ArcGIS 10.2;数据的统计及相关分析采用SPSS 19.0。矢量数据采用GCS_WGS_1984坐标系统,Albers正轴割圆锥等面积投影,中央经线为105°,标准纬线为25°和47°。

2.3.1 幂参数的选择

采样点在插值计算过程中所占权重的大小受p(正实数)的影响,随着采样点与插值点之间距离的增加,采样点对插值点影响的权重按指数规律减小。幂值大,较近点的影响较大;幂值小,较远点的影响增大。

2008—2017年的全国近岸海域环境质量公报显示,无机氮和活性磷酸盐是中国近岸海域海水中的主要污染物,选用代表性较强的无机氮测值进行不同幂参数插值实验,并采用10折交叉验证法[12]验证各插值模型的预测效果。采用模拟值与实测值的平均误差(ME)、平均绝对误差(MAE)、均方根误差(RMSE)和Pearson相关系数等指标评价插值模型的精度。ME反映总体模拟值误差的大小,MAE反映模拟值存在的可能误差范围,RMSE反映模拟值的灵敏度和极值效应,相关系数反映预测值和模拟值之间的吻合程度[8]。

10折交叉验证法的原理是使用系统抽样法把数据集平均分为10份(做法是将所有点位以纬度为第一字段、经度为第二字段进行升序排列,并循环进行1~10编码,该法基本可以保证编号一致的点位均匀分布,将编号一致的数据集视为同一份数据)。每次从10份数据集中拿出1份作为验证集,剩下的9份数据集作为训练集,对每份训练集进行IDW插值实验。

2.3.2 参与插值的周围采样点的数量选择

参与插值的采样点数量与相邻点位的数量相关,通过监测点位所构成的泰森多边形判断一个离散点与哪些离散点相邻,每个泰森多边形的边数代表了相邻点位的个数。选取相邻点位个数出现频率最高的作为插值采样点数量参数,基本可以保证待插值点周围各个方向的近距离点位均参与插值。

2.3.3 插值栅格大小的选择

栅格的大小决定了插值结果对水质类别反映的精细程度,也影响插值运算效率。通过一系列栅格大小梯度实验,对比选取能较好反映水质类别分布情况且运算效率较高的栅格值。

3 结果与讨论

3.1 参数选择结果

3.1.1 幂参数

幂参数p依次选取0.5、1、2、3、4、5、6进行2个研究区无机氮插值实验,并进行误差分析,详见表1和表2。

表1 全国近岸海域10折交叉验证误差结果Table 1 The precision analysis of China coastal area with 10% off cross-certification

表2 渤海海域10折交叉验证误差结果Table 2 The precision analysis of the Bohai Sea area with 10% off cross-certification

对全国近岸海域来说:

ME:p(1)

MAE:p(3)

RMSE:p(3)

相关系数:p(3) >p(2) >p(4) >p(5) >p(1) >p(6) >p(0.5)。

对渤海海域来说:

ME:p(4)

MAE:p(1)

RMSE:p(2)

相关系数:p(3) >p(2) >p(4) >p(5) >p(6) >p(1) >p(0.5)。

2个研究区各模型的预测值与实测值在置信度0.01水平上均显著性相关,说明IDW插值结果与实测结果相吻合,能够反映海水污染物的水平分布情况。

为便于查看,放大2个研究区典型部分来对比不同幂参数的无机氮插值效果(图3)。可以看出,p取值为0.5时,过于强调较远点的影响,各类水质边界不平滑,部分区域插值结果失真(如辽东湾东北部区域二类海水与四类海水直接相邻,缺乏过渡);p取值为1时,上述情况略有改善,海水类别跳跃分布现象仍然存在,并且出现围绕监测点的牛眼现象;p取值为2或3时,各类水质边界趋于平滑,能较好反映各类水质过渡情况;p≥4时,插值效果较之前变化越来越不显著,但仍可以看出,随着p取值增大,各类水质边界更加平滑,但过于强调较近点的影响,当一个区域缺少点位时,评价结果开始从周围向这些区域过度挤压(如长江口邻近海域三、四类海水面积几乎被挤压殆尽,渤海普兰店湾附近二、三、四类海水也存在同样的情况)。

综合考虑各模型误差结果及图上效果,全国近岸海域幂参数取值为3,渤海海域幂参数取值为2。从表1和表2可以看出,虽然渤海海域点位密度较大,但全国近岸海域插值结果误差分析多项指标均优于渤海海域,这说明IDW法在开放海域的模拟效果整体优于近封闭海域。究其原因,海域水质变化虽然都呈现近岸到远海逐渐变好的趋势,但是对于开放海域来说,基本是单一方向的变化,对于近封闭海域来说则是从四周向中间的多方向变化,点位监测值规律性较差,影响了插值效果。

3.1.2 参与插值的周围采样点数量

根据监测点位的分布情况,建立泰森多边形。2个研究区每个监测点位相邻点的数量(N)从1个到8个不等,全国近岸海域以4个为最多,渤海海域以5个为最多,见表3。可将全国近岸海域参与插值的周围采样点数量设置为N=4,渤海海域设置为N=5。

图3 不同p值所对应的无机氮水质类别插值结果Fig.3 Spatial interpolation results of dissolved inorganic nitrogen(DIN) with different power

表3 监测点位构成的泰森多边形相邻点位统计Table 3 The neighbor point statistics using Thiessen polygons created by monitoring points

3.1.3 插值栅格大小

分别将栅格大小(C)设置为1 000、500、200、100 m,比较评价范围吻合程度、统计结果的稳定性及运行效率等(表4)。

栅格大小为1 000 m时,锯齿化情况严重,一些海湾和岛屿的陆上部分被侵占,随着栅格取值的减小,锯齿化情况不断改善,直至栅格大小为200、100 m时,各类别之间的界线平滑,插值结果与评价范围吻合较好。统计结果的稳定性方面,从表4可以看出,随着栅格变小,插值后各类别海水面积占比之间的差异也逐渐变小,栅格大小为200、100 m时,差异基本可忽略不计,这种差异对于评价范围及其子集之间统计结果的稳定有重要意义,如全国与各省(区、市),各省(区、市)与各市,差异越小,越有利。运行效率方面,随着栅格的变小,运行效率不断降低,但均在可接受的范围内。

因此,对全国或海区进行评价时,为能较好反映小水域详细的水质状况,保证插值结果的稳定性,同时确保较高的处理效率,栅格大小设置为200 m左右为宜。

当对河口、海湾等小区域进行单独评价或监测点位密度更大时,也可将栅格设置的更小。同一次评价单个指标图层栅格大小必须保持一致,以便进行多图层叠加。

表4 不同栅格大小对应的各类别海水面积比例差值Table 4 The interpolation results comparison to different cell sizes %

3.2 与常用评价方法的对比

利用IDW插值法评价全国近岸海域和渤海海域水质类别空间分布(29个监测项目)。其中全国近岸海域插值参数设置为p=3,N=4,C=200 m;渤海海域插值参数设置为p=2,N=5,C=200 m。统计海水优良面积比例,并与常用的点位法和代表面积法的统计结果进行比较,结果见表5。

表5 不同方法获得近岸海域海水优良比例对比Table 5 The percentage of the first and second class seawater by different methods %

3种方法中,IDW插值法优良比例值高于点位法,低于代表面积法(泰森多边形)。点位法优良比例相对较低的原因是该法认为每个点位权重是一样的,但监测点位分布不均匀,对评价区域而言,密集区域的单个点位影响范围小于稀疏区域,而水质较差点位主要分布在靠近陆地的点位的密集区域,因此监测点位分布均匀性越差,两者差异越大。

2种面积法水质类别统计结果非常接近,但从空间分布直观展示效果来看(以渤海海域为例,图4),IDW插值法更为流畅,更符合污染物扩散规律。

图4 水质类别空间分布效果对比Fig.4 The seawater quality spatial distribution using IDW interpolation and thiessen polygon

3.3 IDW插值法应用注意事项及方法局限性

3.3.1 方法应用注意事项

3.3.1.1 基础数据监测不全的处理

当某点位存在部分评价项目未监测的情况时,该点位不应参与该缺失项目插值。需要在数据资料准备阶段,将表中相关单元格设置为<空>值(注意并非赋为0值)。

3.3.1.2 插值范围及掩膜范围

插值范围的设置是为了划定数据分析范围,应不小于插值点位外接多边形范围(保证相关点位均参与插值)及评价区域范围(保证评价区域内全部被插值)。一般可设置为评价区域范围图层,也可根据上述原则另行绘制。

若不设置插值范围,则默认输出区域为所有点位的外接矩形,当评价区域涉及海湾时,由于监测点位位于海湾内部,不能将评价区域完全包含,会出现插值不完全现象,部分评价区域存在空白,设置插值范围后避免了此类问题。

掩膜范围的设置主要是因为插值输出数据集的范围通常会大于评价区域范围,为获得以评价区域为边界的插值结果,需要对其进行掩膜处理,也就是使用评价区域对插值结果进行裁剪。掩膜范围应设置为评价范围图层。

3.3.1.3 跨岛问题的处理

为去除半岛及岛屿两侧点位之间的相互影响,需要对其进行阻断,使其无法参与对方区域的插值运算,这个目标可通过绘制限制采样点搜索的隔断线,并指定为障碍参数来实现。隔断线绘制的原则是插值点与用来进行插值的采样点之间地理上应相互通视。图5展示的是广东东海岛附近海域插值结果。东海岛南侧2个点位均为一类水质,北侧1个点位为四类水质,未设置隔断线时,东海岛南侧海域受北侧四类水质点位的影响,临岸有一部分区域为二类海水,这与实际情况相悖;设置隔断线后,南侧海域均为一类水质,使插值结果更接近实际。需要注意的是,使用隔断线会显著增加系统处理时间,因此,隔断线不是越多越好,可根据实际需要进行设置。

图5 设置隔断线对评价结果的影响对比Fig.5 The comparison with or without barrier

3.3.2 方法局限性

IDW插值法的特征之一是会产生围绕观测点位置的“牛眼”,这不符合海水污染物的局部扩散规律。

IDW插值法的特征之二是其预测值处于所输入的采样数据的最小值与最大值之间。这个特点虽然保证了插值结果不会出现过高或过低的超失真情况,但同时也抑制了污染物空间分布预测向高值或低值推进的趋势。研究中对各验证集预测值及对应实测值的水质类别进行对比分析,发现分布在内侧或者外侧、测值明显高于或者明显低于周围点位的测值容易产生预测错误,这主要是由于该特征所决定的。

IDW插值法公式与任何物理过程都不关联,但海水有海流、潮汐、潮流、波浪等流动特性,且具有很强的季节性,污染物的扩散极为复杂,该法未将水动力条件考虑在内。

4 结论

关于IDW插值参数的选取,应随不同区域、不同点位分布而不同。幂参数(p)可通过多幂值实验,利用交叉验证误差分析选取适宜的值;参与插值的周围采样点数量(N)可通过构建泰森多边形统计相邻点位个数出现频率最高的值确定,点位密度小,N较小,点位密度大,N较大;插值栅格大小(C)需综合考虑评价范围大小、评价结果稳定性和运行效率选定。针对目前研究区范围及监测点位的分布情况,全国近岸海域参数设置为p=3,N=4,C=200 m较为适宜;渤海海域设置为p=2,N=5,C=200 m较为适宜。

IDW插值方法在开放海域和近封闭内海海域均适用,但在开放海域的模拟结果好于近封闭海域。

海域水质面积评价的IDW插值法与代表面积法统计结果接近,但前者在空间分布展示上更接近实际;点位分布的均匀性决定点位法与面积法评价结果的接近程度,点位分布越均匀,两者越接近。

猜你喜欢
插值法栅格点位
InSAR形变场最佳插值算法对比研究
基于邻域栅格筛选的点云边缘点提取方法*
基于结构光视觉的钻孔点位法矢检测技术研究
基于A*算法在蜂巢栅格地图中的路径规划研究
《计算方法》关于插值法的教学方法研讨
《计算方法》关于插值法的教学方法研讨
浅谈舞台灯光工程配电回路设计
大盘仍在强烈下跌趋势中
克里金插值法内插IGS电离层图精度分析
不同剖面形状的栅格壁对栅格翼气动特性的影响