土壤属性空间插值方法及其精度验证

土壤属性空间插值方法及其精度验证

作为土壤变化的时空定量监测方法，土壤属性空间插值方法及其精度是土壤科学和“数字土壤”（数字土壤）领域的一项重要研究和方向。由于土壤母质、土壤类型、地形、气候、生物等环境条件的差异及人类活动的影响,土壤属性存在空间异质性。地理信息系统(GIS)、地统计学和曲面建模逐渐成为适合农业、环境管理等研究领域的常用工具。近些年来,由于精准农业的实施和土壤污染管理的需要使得这些技术进一步发展。例如,土壤肥力制图的质量影响定位土壤肥力管理的效果,土壤重金属分布图的质量影响区域环境保护措施的具体实施等。土壤属性空间插值能够为土壤的可持续利用及环境保护的宏观决策提供技术支持,但在一定程度上依赖于土壤属性空间插值的精度。因此,土壤属性的空间插值方法、精度及其验证是土壤属性空间插值中必须考虑的因素。本文系统总结了土壤连续属性空间插值及其精度的研究进展,并阐述了提高土壤属性插值精度的方法,最后探讨了未来的研究趋势。1土壤属性空间插值方法1.1变量的半方差函数拟合过程和正态分布检验Kriging插值法是一种无偏线性最优估值方法。对变量在点x处的估计值Z*(x),可以通过该点影响范围内的n个有效观测值Z(xi)的线性组合得到,即:式中:λi是赋予观测值Z(xi)的权重,表示各观测值对估计值Z*(x)的贡献,在保证估计值无偏性(即估值偏差的平均值为0)和最优性(即估值方差最小)条件下,可由变量的半方差函数计算得到。其半方差函数可用下式表示:其中:γ(h)是半方差,N(h)是距离等于h时的点对数,Z(xi)是在位置xi处的数值,Z(xi+h)是在距离(xi+h)处的数值。在进行半方差函数的拟合之前,需对数据进行正态分布检验,因为数据的非正态分布会影响半方差函数的精度,即当样品的平均值增加时,样品的方差也增加,使实际变异函数产生畸变、基台值和块金值增大,导致某些结构特征不明显。如不呈正态分布,需对其进行变换。正态变换的常用方法主要包括对数变换、Box-Cox变换和正态计分变换(NormalScoreTransformation)等。Kriging方法包括单变量Kriging方法和与辅助变量信息相结合的Kriging方法。不同Kriging方法的特点如表1所示。1.2距离待估点近的样点的权重反距离权重法(IDW)假设未采样点的属性值是其周围一定邻域内已知点的加权平均,距离待估点越近的样点,被赋予的权重越大。设Z*(x0)是变量Z在未采样点x0基于周围数据Z(xi)的估计值:其中wi是每个Z(xi)的权重,n是用来估计的邻域点的个数。反距离权重为:其中di是估计点和采样点之间的距离,p是指数参数。1.3曲面滑动地连接Spline样条函数以最小曲率面来充分逼近各采样点,包含的多项式可以将一系列线和曲面光滑地连接在一起。这种多项式拟合最大的缺陷在于除了通过采样点外,没有其他任何限制,导致所拟合的值在采样点与周围的点差距可能较大。因此,样条插值法对渐变的表面会有较好的插值结果,而对在很短的水平距离内曲面值有较大变化的情况通常不太合适。1.4bme方法的估计模型贝叶斯最大熵(Bayesianmaximumentropy,BME)最早由Christakos提出,是一种能够利用共区域化模型和可获得的软数据来重建空间变异性的方法。BME方法后来发展为可同时考虑不同类型和性质的数据,结合较少的硬数据和较易获得并且信息量较大的软数据(数据范围可以是精确的测量值,也可以是区间值或者概率密度分布,如专题地图、数值区间、概率密度函数、物理定律等)对目标变量的后验条件概率进行预测估计。BME作为一种非线性方法,可以提供每个估计点完整的概率分布。BME的估计模型为:式中:xmap=(xk,xhard,xsoft),xk、xhard、xsoft分别表示待估点的预测值、采样点的硬数据和软数据;fG(xmap)为先验概率,是在满足信息熵H最大的条件下利用拉格朗日乘数法求出的;fG(xhard,xsoft)是硬数据和软数据的联合概率密度函数,是一个归一化系数;fk(xk|xhard,xsoft)为待求的后验概率密度函数。对于连续变量,信息熵H的定义是:近10多年来BME方法的应用越来越广泛,它的主要特点是发挥了该方法在严密的逻辑下使用软数据的优势。如D’Or等基于贝叶斯最大熵(BME)方法仅利用软数据便达到了评价变量空间变异性的目的。Douaik等基于区间软数据,应用贝叶斯最大熵方法对土壤盐分进行了制图,结果表明BME方法可以结合较易获取的软数据信息,取得更准确的预测。Bogaert和D’Or将BME算法和蒙特卡罗过程结合(BME/MC),利用少量采样点的硬数据和详细的土壤图进行了土壤黏粒、粉粒和砂粒含量的空间预测,研究表明,BME的估计值不仅比简单Kriging准确,而且结果更能代表土壤的真实性。BME方法与传统的地统计学方法相比具有以下优势:(1)软数据无需“硬化”,对不同类型的软数据有相对灵活的表达方式;(2)能够以先验知识的方式有效利用不同来源的数据,从而提高空间插值精度;(3)是一种无需正态分布或空间均质假设的时空地统计学方法;(4)既可单独处理连续变量或范畴变量,也可同时处理以上两种变量。1.5gwr方法空间变量之间存在的空间非稳定性是简单的全局模型所不能解释的,Brunsdon等在1996年提出的地理加权回归(Geographicallyweightedregression,GWR)方法,是近10多年来国内外公认的模拟空间变异过程的有效方法。在GWR方法中,每个数据点根据回归点被赋予一个反距离权重,与回归点离得越近的点具有越高的权重。因此,该方法对权重函数的选择尤为敏感。有关GWR的理论在文献中有更多细节的描述,其参数可用GWR3软件进行预测。GWR模型可被描述为:其中:Yi是预测的在i点位置的土壤某种属性值,(ui,vi)是i点位置的坐标,β0到βk是回归系数,Xi1到Xik是i点位置的环境变量值,k是环境变量的个数。Mishra等应用GWR方法预测了区域土壤有机碳库,结果表明GWR的预测结果比多元线性回归和回归Kriging更准确,得到了土壤碳库预测的满意精度。1.6土壤属性的高精度曲面建模HASM是应用采样点数据通过一系列迭代计算对土壤属性曲面进行全局模拟,进而对未采样点的土壤属性进行插值的方法。假设{u珔i,j}是采样点{(xi,yj)}在曲面u上的采样值,uin,j(其中n≥0,0≤i≤I+1,0≤j≤J+1)是中心点为(xi,yj)的栅格值的第n次迭代值,其中ui0,j=u珘i,j,{u珘i,j}是基于采样点{u珔i,j}的插值结果。第n+1次迭代的有限差分基本表达式为:其中:;A1、A2和b1n、b2n分别是式(8)和(9)的系数矩阵和右端常数项矩阵。假设为了保证采样点所在的栅格模拟值与真实值相等或相近,HASM的表达式可转化为等式约束的最小二乘问题:其中C和d分别为采样点系数矩阵和采样点的值,它们的表达式分别为即第k个采样点(xi,yj)的值为为了求解最小二乘问题的方程组(10),引入一个足够大的参数λ(λ是赋予采样点的权重,它决定于采样点对模拟曲面的贡献),式(10)可被转化为无约束的最小二乘问题:也就是求解近年来,随着HASM方法的不断完善[37,38,39,40,41,42],土壤属性的高精度曲面建模从精度、速度、尺度各方面得到了长足发展。如Shi等将HASM应用在土壤pH的插值中,获得了比传统土壤属性插值方法更高的精度;在此基础上,Shi等提出了结合土地利用类型作为辅助信息的改进的HASM模型(HASM_LU),该方法不仅提高了土壤属性的插值精度,而且提供了更多的细节信息和更准确的空间结构信息。基于多重网格(MG)解法的HASM方法(HASM-MG),在保证精度的基础上,同时也提高了土壤属性插值的运算速度。李启权等应用HASM模拟了中国范围内土壤表层有机质和全氮的空间分布,得到了较为合理的插值图。2土壤服务验证方法精度提高土壤属性空间插值精度的验证受验证方法和验证指标影响,选择合适的验证数据集和验证指标,对准确评价土壤属性空间插值方法的精度至关重要。2.1验证方法的选择2.1.1应用交叉验证交叉验证包括两种,一是去一法交叉验证(cross-validationwithreplacement),该方法通过每次消去一个目标样点,用所有剩余样点来估计目标样点的方法;二是k折交叉验证(kfoldcrossvalidation),将训练集合分成k个大小相同的子集,其中一个子集用于测试,其他k-1个子集用于训练。这样,整个训练集中的每一个子集均被预测一次。因此,去一法可以说是k折交叉验证法的极端情况。交叉验证经常用来比较不同方法或同一方法不同参数选择的插值结果,所有数据都可用来预测,具有快速且经济的特点。例如,可用其从可供选择的半方差模型中选择最佳搜索半径和步长距离以使Kriging方差最小化,也可用来帮助测试IDW和Spline方法的最佳参数。另外,交叉验证也被许多学者应用到不同方法的插值精度对比中。但交叉验证在许多情况下不足以描述空间预测误差。例如,Laslett等在一个实例中证明了应用交叉验证可以导致两种方法精度比较的错误结论。以Kriging为例,每去掉一个点,半方差函数应该是发生变化的,而实际预测时,如果每去掉一个点就做一次半方差函数的计算费时费力。Mueller等在对比应用两种数据集计算均方根误差(RSME)时发现,应用所有数据交叉验证远比应用独立数据集验证的RSME大。2.1.2增加验证数据集的参独立数据集验证是估计空间不确定性更为直接和独立的方法,但需要额外采集验证数据集,预测成本较高。具体步骤包括三步:一是随机将数据集分成预测数据集和验证数据集;二是用预测数据集的值对验证数据集的值进行模拟;三是比较验证数据集每个点的模拟值和真实值。2.2克立格方差se、不精确,imsp,fmsp,akse和平均克立格标准表土壤属性空间插值精度的验证指标可分为绝对验证指标和相对验证指标。其中绝对验证指标主要包括平均误差(MeanError,ME)、平均绝对误差(MeanAbsoluteError,MAE)、均方根误差(RootMeanSquaredError,RMSE)、不精确度(Imprecision,IP)、平均克立格标准误差(AverageKrigingStandardError,AKSE)、均方根标准化预测误差(RootMeanSquareStandardizedPredictionError,RMSP)、平均标准化预测误差(MeanStandardizedPredictionError,MSPE)和被减少的克立格方差(ReducedKrigingVariance,RKV)等;相对验证指标包括预测效率(PredictionEfficiency,PE)、预测方法的级别(Meanrank,MR)和标准差(StandarddeviationoftheRank,SDR)、相关提高度(Relativeimprovement,RI)、散点图、残差图、半方差函数的残差等。2.2.1验证点位置预测的精度分析在对比不同插值方法的过程中,平均误差(ME)、平均绝对误差(MAE)和均方根误差(RMSE)被看作是三个最常用的衡量已知数据和预测数据差别的指标,根据每个验证点位置的观测值和预测值计算,公式如下:其中:N为验证点个数;z(xi)为验证点的测定值;z*(xi)为验证点的模拟值。ME是预测无偏的量度,越接近0越说明方法是无偏的;MAE和RMSE是预测精度的量度,越小则说明插值方法越精确。由RMSE和ME组成的“不精确度”(IP),也可以用来量测估计的质量,其中:IP(x0)是观测点x0位置的不精确度;RMSE(x0)是观测点x0位置的均方根误差;ME(x0)是观测点x0位置的平均误差。2.2.2克立格方差和被减少的克立格方差平均克立格标准误差(AKSE)、均方根标准化预测误差(RMSP)、平均标准化预测误差(MSPE)和被减少的克立格方差(RKV)4个指标用到了克立格方差,因此只适用于克立格方法。它们的计算公式分别是[18,20,21,53,58,59]:其中:N是数据集中点的个数,σ2是xi位置的克立格方差,MSE为均方差(MeanSquaredError),MKV是平均克立格方差。2.2.3少均方差mse的比选预测效率,又称降低误差指数,表明预测方法与田间平均相比提高的效率百分比,属于相对验证指标,表示某种方法与田间平均方法相比减少均方差(MSE)的百分比,可比较不同方法预测不同数据集、不同元素之间的效果,使预测方法在不同数据集、不同元素的插值效果比较成为可能[5,53,56,60,61,62]。定义如下:其中:MSEfieldavg是用田间采样点平均值作为所有验证点位置的估计值计算出的MSE;Zm是所有采样点的土壤属性的平均值;n为验证点个数;Z*(xi)为验证点的模拟值;Z(xi)为验证点的测定值。2.2.4各点中异常值的影响不同方法预测的准确度级别(MR)和级别的标准差(SDR)两个指标考虑到了离群值对插值方法精度的影响,能够适合多种方法插值精度的比较,排除验证点中异常值的影响,而ME、MAE和RMSE等指标并未考虑到以上影响。Ri和Si的计算公式如下:其中:第i种方法的平均级别(MR)是Ri,第i种方法的级别的标准差(SDR)用Si表示;n为验证点个数;l表示插值方法的种类;rij是第i种方法在第j个点的级别,rij的取值范围从1取到n,其中准确度最高的级别设为1,最低则设为n。最优方法应具备最小的平均级别和较小的级别标准差,可用平均级别(Ri)与级别的标准差(Si)的散点图来识别。2.2.5对数kring方法精度的提高Kravchenko和Bullock用相关提高度(RI)来表示对数Kriging比Kriging方法精度提高的程度,发现对数Kriging方法比Kriging方法是否提高精度与样本数有关。RI的计算公式为:其中:RKOlog和RKO是对数Kriging和Kriging的预测值和测量值的相关系数。2.2.6散点图的选取预测值和测定值的散点图经常被用来检查和评价预测质量,如果图的质量较好,数据的散点图会靠近1∶1的线。残差(预测值-观测值)图或半方差函数的残差可用来检查Kriging方法预测误差在空间上是否存在变异和存在怎样的变异。3提高土壤属性插值精度的方法3.1idw插值比对空间插值方法的选择是影响土壤属性预测图质量的最主要因素之一。前人已对Kriging、IDW和Spline三种常见的土壤属性插值方法的精度比较做了大量研究,但结论并不一致。一部分研究表明Kriging的效果好于IDW和Spline;而另一部分结论却恰恰相反;有时方法的精度优势也常常是互相混淆的[5,20,54,60,61,67]。例如,Panagopoulos等认为Kriging插值对于每种土壤属性来说均是最好的方法。但也有学者研究表明IDW插值法优于包括Kriging插值在内的其他方法。然而,Kravchenk认为具有已知半方差函数参数的Kriging插值的效果明显优于IDW插值,但如果从采样数据集中不能获得可信赖的半方差函数,精度就不如IDW。除了对Kriging、IDW、Spline等方法进行比较外,也有学者对不同Kriging插值方法进行了比较。Triantafilis等比较了普通Kriging、回归Kriging、3D-Kriging和协同Kriging方法对土壤盐分插值的结果,表明回归Kriging精度最好,而从平均值和标准差的等级来看,则协同Kriging表现最好。近几年的研究表明,在已有的应用实例中,HASM已表现出精度明显优于以上传统插值方法(Kriging、IDW和Spline)的特点。3.2插值精度分析同一种插值方法中,不同参数的设定对插值精度有显著影响。影响Kriging插值的因素包括:(1)数据的变异性和空间结构;(2)半方差函数模型;(3)搜索半径;(4)估计的邻域点数目等。影响IDW方法的参数有:(1)指数值的选择;(2)估计的邻域样点数目;(3)搜索半径。影响Spline精度的因素有:(1)样条方法;(2)权重值;(3)邻域点数目。仔细选择以上参数,均会不同程度地提高Kriging、IDW和Spline三种插值方法的土壤属性插值精度,但这些方法的绝大部分参数的确定没有太多可遵循的规律[20,53,60,65,66,68,69]。虽然已有研究表明数据的一些描述性统计特征与插值方法的部分参数确定有些联系,但尚未发现比较确定的规律,以上三种方法的大部分参数仍需要用交叉验证等方法做大量的尝试才可确定。例如,Kravchenko和Bullock用IDW方法对土壤属性进行插值时对幂指数的值作了仔细选择之后精度有了显著提高,发现大部分偏度大于2.5的数据选择4作为幂指数精度最高,偏度小于1.0的数据幂指数选择1精度最好;而Weber和Englund的研究结论恰恰相反,即较大的幂指数适用于较小偏度的数据。Gotway等发现变异系数小于25%的数据用IDW插值时选择较高的幂指数有较高的精度,反之幂指数选择1精度最好;而Kravchenko和Bullock则发现IDW插值中的幂指数选择与变异系数无显著相关性。3.3采样点分布的稳定性采样数目和密度的设计同样可以影响图的质量。较低密度的采样降低了空间结构的描述能力,增加了预测误差。一般地,采样点数越多,土壤属性图就越准确。对于Kriging插值来说,准确的半方差函数的构建是描绘准确的空间结构的前提,而获得稳定的半方差函数至少需要50~100个采样点;也有学者建议至少100~150个采样点数据才可达到一个平稳的变异函数;超过300个样点才适合分析各向异性。3.4土壤属性与空间自相关范围空间结构对于制图准确度的重要性已被广泛认识。数据的空间结构可用地统计学中的半方差函数描述,参数包括块金值/基台值(N/S)和空间自相关范围(变程)。其中N/S表明研究区不可被半方差函数描述的小范围变异的部分;空间自相关范围表示超过变程范围时土壤属性值之间没有相关性。空间自相关性越强,土壤属性制图越准确。较小的N/S和较大的空间自相关范围的土壤属性制图通常具有较高的准确度。3.5空间分布的变化实验分析顺序对土壤属性预测图所表现的空间分布准确性也有影响。按土壤采样的格网顺序或野外采样顺序分析与随机分析采样点,得出的土壤属性空间趋势会有差别。这是因为,这种采样顺序分析的相关性打乱了空间自相关的关系,导致较多异常值的出现。Mueller等以步长距离为横轴,实验半方差值减去真实的半方差值为纵轴作图,得到一些异常的土壤属性的空间分布趋势,这些异常有可能是由于数据的化学实验分析顺序的趋势所导致的。因此,在实验室分析之前,将采样编号随机编排对研究数据的空间结构和制作土壤属性图很重要。3.6土壤金属因素和土壤电导率的影响当土壤属性值在空间上相关性很差或变化比较剧烈时,用其他相关的连续变量或类别信息可以提高插值精度。土壤属性与很多连续变量存在重要的相关关系,利用DEM、遥感数据、其他与目标变量相关的土壤属性等作为辅助变量,均得到了较为满意的效果。O