GWR模型下农用地土壤镍空间分布预测

王春帅，姚立伟，刘弋珲，牛瑞卿，任超

(1.河南省地质矿产勘察开发局 第一地质矿产调查院，河南 洛阳 471000；2.中国地质大学(武汉) 地球物理与空间信息研究所，武汉 430070)

0 引言

随着《中华人民共和国土壤污染防治法》的颁布与执行，土壤污染现状的调查工作变得越来越重要，而农用地土壤是食品安全的第一道关口，因而需要加强对农用地的监管工作。但是传统的土壤重金属监测方法需要耗费大量时间、人工成本，具有高成本、低效率、不能大范围进行快速检测的特点。利用光谱及各类其他因子进行重金属浓度快速反演的方法应运而生，这种方法具有高效、快速、无损、无害、可重复、成本低、尺度大的优点。利用光谱特征来拟合区域重金属污染空间分布的研究已经成为热点，常用的方法有最小二乘法回归、偏最小二乘回归、逐步回归、支持向量机、随机森林等，但是对于线性模型的使用主要以全局模型为主，局部模型拟合的研究很少[1-3]。

沈强等[4]通过对光谱曲线进行连续去除和一阶微分预处理，选择合适的高光谱数据进行逐步回归，建立了对矿业废弃地土壤重金属含量进行高光谱反演的模型；张秋霞等[5]利用偏最小二乘模型对高标准基本农田区域重金属含量进行高光谱反演，其模型决定系数达到了0.70～0.99，获得了极好的反演效果；陶超等[6]利用支持向量机模型对土壤重金属铅、锌进行高光谱反演，取得了较好的模型精度；徐良骥等[7]采用逐步回归、偏最小二乘回归、人工神经网络的方法分别对煤矸石充填复垦重构土壤中重金属含量进行高光谱反演，发现人工神经网络模型取得了最好的建模效果。可见，对重金属浓度反演所用到的主要为高光谱遥感影像，而Landsat-8因为数据源的开放性，使用更加普遍，利用Landsat-8进行反演重金属的研究已经取得了一定的研究成果。王海潇等[8]利用DEM数据对Landsat-8影像进行地形校正，然后利用校正后的反射率变换作为变量，利用PLSR模型进行建模反演重金属Cu含量，模型决定系数达到0.52，有较好建模效果；方媛[9]利用多时相Landsat-8影像，使用偏最小二乘回归、人工神经网络、支持向量机模型进行无偏评价与选择，发现偏最小二乘模型在土壤重金属反演中表现最佳，并利用该模型对Cu、As、Cr含量进行反演。

综上所述，对于土壤中重金属浓度的反演主要考虑光谱特征，很少考虑样品采集点的空间位置及地形信息。因此，本研究综合考虑光谱变换信息、采样点空间位置信息、地形信息等，选用逐步回归和地理加权回归，实现土壤中重金属镍含量的快速、精确反演，同时为以后的研究提供新的思路。

1 研究区域概况

研究区石宝沟位于栾川县石庙镇，地理范围为111°31′9.47″E至111°33′6.47″E，33°48′53.76″N至33°51′50.21″N之间，面积为3 km2。农用地污染土壤采样范围为0.50 km2，污染地块主要分布在石宝沟及银洞沟，呈条带状分布于石宝路两侧。调查区北至石宝沟三岔口，南至石宝沟沟口，与伊河河谷相接。研究区交通便利，有S328省道(栾川县城至石庙镇)，区内主干道路为石宝沟道路。研究区位于北至马超营断裂和南达栾川断裂之间的异常区，该异常区呈现为狭长的带状形态，属于卢氏—栾川钼、钨、铅、锌、银金地球化学异常带，其中广泛发育的元素有Mo、Ag、Pb、Zn、Au、W、Cu、Mn、Ni等。沿马超营断裂带，形成Au、Ag、Pb、Zn、Mo、Bi多元素的综合异常带，Ag、Pb、Zn元素异常套合较好。不同类型的岩体与岩体的不同部位分别形成钼、钨、银、铅、锌、铁、铜、锌、硫、金矿等，构成一个完整的成矿系列[10-16]。

本研究区成矿地质条件极为优越，存在丰富的矿产资源，但是由于历史原因，曾经有多达24个采矿权，无序开采严重，矿产开采过程中长期存在大矿小开、一矿多家分块开采现象，资源利用效率低，并且不可避免地对研究区域内农用地产生污染，导致该区域农用地重金属含量异常。

2 数据来源及研究方法

2.1 数据来源及处理

本研究根据石宝沟和银洞沟农用地土壤分布特点，以自然分割的地块划分调查采样区[17-18]。采样田块共计60个，每个地块构成一个单独的采样区，每个采样区划分若干个采样单元，采用1∶1 000田块测量成果图对采样区进行网格化剖分。在石宝沟主沟、支沟银洞沟、支沟倒回沟沟口农用地土壤污染区共布设表层土壤样点277个。将采集的样品依据土壤环境监测技术规范样品流转的要求，运输至实验室进行镍浓度的化验分析。采用的分析方法为石墨炉原子吸收分光光度法。地块划分及检测结果如图1所示。

图1 研究区农用地Ni空间分布

本研究采样时间为2018年7月末至8月中旬，从美国地质调查局网站(www.usgs.gov)下载了研究区2018年8月的Landsat-8影像，并利用ENVI 5.3对图像进行预处理，包括辐射定标与大气校正。利用BIGMAP下载了研究区1 m×1 m的DEM数据，该数据精度较高，用于计算坡度、坡向等地形因子。

2.2 研究方法

土壤的反射光谱包含了丰富的土壤信息，是土壤形状的综合反映，在近红外波段范围(700～1 000 nm)内，土壤光谱的吸收特性主要是由于金属离子的电子跃迁形成；而在短波红外区域(1 000～3 000 nm)，土壤的吸收主要是由于有机质、碳酸盐、层状硅酸盐等矿物的各类分子团中化学键的弯曲、伸展、变形等振动造成。重金属元素在土壤中通常存在多种化学形态，不同的形态具有不同的敏感光谱特征。通过多波段的组合运算，将不同波段的光谱信息进行相互补充，使得组合后结果与重金属元素浓度的相关系数显著高于单一波段[19-23]。

本研究通过提取研究区域农用地的可见光-近红外-短波红外的波段信息，并进行一定的波段组合运算，结合地形因子、与污染企业的空间位置距离因子等，选择不造成明显多重共线性的解释因子，使用地理加权回归的方式进行建模，并对所建模型进行评估，最终对整个研究区域进行重金属镍反演，评价其空间分布特征。技术流程见图2。

图2 技术流程图

1)变量选取。首先，利用Landsat-8遥感影像提取研究区农用地的各种植被指数作为因子，包括归一化植被指数(NDVI)、绿色归一化植被指数(GNDVI)、比值植被指数(RVI)、增强型植被指数(RVI)等共计9种；其次，对Landsat-8影像的7个波段进行不同方式的波段运算，包括对波段取对数，开根号，任意2个波段之间的加、减、乘、除运算，以及任意3个波段之间的加、乘运算等；然后，考虑到研究区多年矿产背景，农用地中重金属可能来源于企业排放废气的大气沉降作用、尾矿库泄露等，利用ArcGIS及DEM数据提取采样点的高程、坡度、坡向、地面曲率因子、剖面曲率因子、平面曲率因子、坡度变率因子以及坡向变率因子等地形因子；最后，利用ArcGIS的近邻分析工具，分别计算每一个采样点与尾矿库、矿产加工企业以及主干线(石宝路)的最近距离，将结果作为分析因子。

综上，本研究共获取的自变量因子有243个，计算每一种因子与镍浓度的皮尔逊相关性系数，筛选出显著性水平α=0.01的变量因子33个。

2)OLS模型。OLS模型利用n组观察值求得p个独立变量与因变量之间残差平方和最小的拟合方法是一种利用n组观测值，使得p个独立变量与因变量之间的残差平方和取得最小的拟合方法。计算方法如式(1)所示。

(1)

式中：Y是农用地土壤中重金属镍的浓度；βi为未知的模型拟合参数；p是用于建模的变量个数；ε为模型的残差项，服从正态分布N(0，σ2)。通过矩阵方程的方式进行估计，如式(2)所示。

(2)

逐步回归是指将所有自变量按其对因变量的影响程度从大到小排序，把通过了偏F统计量检验的自变量逐个引入回归方程；同时，对已引入的变量进行检验，若未通过检验则剔除，直到既不能引入新变量也不能剔除已有变量为止。

3)空间自相关。空间自相关是进行地理加权回归分析的前提条件，一般用于描述某一区域的样本观测数据与其他区域的样本观测数据之间的相互依赖程度，常用来确定数据间是否具有依赖性的指数为全局莫兰指数(global Moran’S I)。

4)地理加权回归。地理加权回归(geographic weighted regression，GWR)模型是对传统线性回归模型的扩展，通过把数据的位置信息通过权重的方式嵌入到回归参数来反映数据在空间上的关系差异，其实质是一种局部加权最小二乘法。基本的GWR模型表达如式(3)所示。

(3)

式中：i=1，2，3，…，n；(μi，υi)为第i个样本数据的位置；γn(μi，vi)为第i个样本数据中的第n个估计参数，是能够反映不同地理位置的函数项；γ0(μi，υi)为估计方程的常数项；δin为第i个样本数据的第n项解释变量；δi～N(0，σ2)；Cov(εi，εj)=0(i≠j)。选择一个合适的最优带宽是构建高质量GWR模型的关键。通过对比验证，本文使用固定距离法构成核函数，并选择最小信息准则AIC指标来决定最优带宽(式(4))。

(4)

5)模型拟合度评估。建立模型后，对测试集数据应用模型，计算应用平均误差(mean error，ME)、平均绝对误差(mean absolute error，MAE)、均方根误差(root mean square error，RMSE)和预测精度(accuracy，Acc)对模型进行精度评价。ME衡量模型预测偏移程度，该值越接近于0，偏移程度越小；MAE和RMSE是衡量模型估算精度的指标，该值越小，模型精度越高；Acc是描述模型预测精度更直观的指标，越接近于1，模型预测能力越好。

3 模型构建与评估

3.1 OLS模型

从所有277个检测数据中随机抽取75%(207组)的数据作为建模组，将剩余25%(70组)的数据作为验证组。

设定自变量进入模型的显著性水平SLE=0.1，保留在模型中的显著性水平SLS=0.05，通过逐步回归的方法最终选出采样点与厂区的最短距离、ln(band2/band3)、band3-band5为基础变量的模型。

模型计算的方程分析见表1，回归模型的系数、显著性及共线性诊断见表2。可见，农用地土壤中镍含量与尾矿库、工矿企业等有关，所建立的模型中所有的变量因子显著性检验均小于0.01，并且VIF值小于10，说明变量因子之间不存在严重的多重共线性。此外，模型决定系数(R2)达到0.325，在一定程度上已达到统计模型所需的基本要求，但是仍然偏低。

表1 方差分析表

表2 OLS回归模型系数、显著性及共线性诊断

3.2 空间自相关分析

利用ArcGIS的空间自相关分析工具，计算研究区用于建模样本的局部莫兰指数，计算结果见表3。由莫兰指数可知，本研究区重金属镍浓度呈空间正相关，根据Z得分可认为该区域空间显著相关。此外，通过显著性水平α=1的显著性检验得知，本研究区呈现显著空间自相关，GWR模型适用。

表3 全局空间自相关分析结果表

3.3 GWR模型

采用与矿产加工企业厂区最短距离、ln(band2/band3)、band3-band5作为模型的解释变量，以高斯函数为核函数，选择使AIC值取得最小值的带宽为模型的最终带宽，带宽为460.97 m。在这个带宽范围内，每个样点平均有25个相邻样点参与建模，模型统计数据如表4所示，模型决定系数达到0.64，取得较好的拟合效果。

表4 GWR模型统计结果

GWR模型获取了每一个样品点的回归模型，其系数的统计结果见表5。根据ArcGIS官方介绍，所计算模型的条件数如果大于30，表示存在多重共线性。由表5可知，本研究建立的模型不存在大于30的条件数，故本模型不存在明显的多重共线性。

表5 GWR模型条件数及各变量系数的统计特征

3.4 模型评估

将277个样本点数据中的验证组数据进行验证，并计算ME、MAE、RMSE、Acc值，计算结果见表6，2种模型验证结果的残差分布见图3。通过对比可以发现，利用地理加权回归模型GWR，明显可以获得更好的拟合效果，其MAE和RMSE计算结果均较低，而且Acc计算结果显示模型预测精度达到了96.51%。另一方面，从残差分布图(图4)也可以看出，GWR的残差分布更加集中，并更接近于0。

表6 模型评估结果

图3 OLS模型反演效果图

图4 2种模型拟合研究区验证样本镍浓度的残差分布

4 反演

利用所建立的模型，对整个研究区农用地进行了反演。利用OLS模型反演的结果见图3，利用GWR模型反演结果见图5，2种模型反演结果的统计特征见表7。

由表7可知，利用OLS模型反演的结果，重金属镍含量范围在20.80～62.03 mg/kg，GWR模型反演的结果中金属镍含量范围在19.61～66.99 mg/kg之间。实际采样检测数据中，重金属镍含量范围在12.36～124.00 mg/kg之间，范围大的主要原因在于研究区内所采集样品中存在少数偏高的离群点。例如研究区偏南侧2个地块，其周围重金属镍含量均属正常状态，只有它们的浓度存在异常(其镍含量分别为124 mg/kg、60.9 mg/kg)，原因可能在于这2个地块位于企业周边(栾川县天罡矿业有限公司)，受到了该企业的影响，或者是由于该地块采集样品的例外性导致了值的异常。而OLS模型和GWR模型均无法对孤立的异常点进行反演。分别比较OLS模型分位数与实测值分位数之间的差异、GWR模型分位数与实测值分位数之间的差异，可见GWR分位数之间的差距更加小，说明了GWR模型的可靠性。

图5 GWR模型反演效果图

表7 全局实测值及各模型计算结果统计特征

由图1、图3(a)、图5(a)可见，研究区中部偏北区域重金属含量均较高。细节见图3(b)、图5(b)。该区域重金属浓度的统计信息见表8。该区域实测值浓度的范围为12.36～84.47 mg/kg，OLS模型模拟计算的结果范围为36.52～60.13 mg/kg，而GWR模型模拟计算的结果范围为39.46～66.99 mg/kg，对于最大值和最小值的模拟效果较差；中间各分位数与平均值相比较而言，GWR取得更好的模拟效果。总之，GWR模型的计算结果可以真实反映该区域重金属浓度的异常状况，该区域重金属浓度异常的原因可能与上游存在的大型矿产加工企业鑫川化工的生产经营活动有关。

表8 研究区中部偏北区域实测值及各模型计算结果统计特征

综上，本研究在研究区范围内选取了一些未进行实际采样的区域进行2种模型的反演，结果发现倒回沟西侧区域重金属镍含量较高(图3(c)、图5(c))。通过对该区域模型反演结果进行统计特征分析(表9)发现，GWR模型反演浓度范围为29.73～49.57 mg/kg，该区域平均值为41.66 mg/kg，含量高于整体的平均值，其主要原因与该区域历史上存在尾矿库有关。该区域农田为尾矿库复填重构土壤，土壤中重金属镍浓度与其他区域相比存在异常。

表9 倒回沟西侧区域各模型反演结果统计特征

5 结束语

本研究通过将采样点与矿产加工企业的最短距离、ln(band2/band3)、band3-band5作为基础变量，以栾川县石宝沟农用地土壤实际检测重金属镍浓度作为目标变量，分别使用逐步回归与地理加权回归的方法进行模型拟合(模型决定因子分别达到了0.325和0.64)，之后利用建立模型前随机抽取的验证样本对模型进行了验证。经过计算，Acc值分别为72.88%和96.51%。综合比较，使用GWR模型获得了更加优秀的拟合效果。

通过对整个研究区所有农用地镍含量的空间分布应用2种模型进行反演，并与实测值进行对比，进一步发现GWR模型的有效性较佳，估算结果与实际情况吻合良好。

对于无实测值区域，通过反演发现，倒回沟西侧区域重金属镍含量异常，其原因和该区域用地历史有关。该区域曾经存在矿产加工企业，同时部分农用地为尾矿库复垦重构所得。