土壤Hg和NDVI全局均值估计的样本点权重调整方法

张东悦 朱庆伟 董士伟 潘瑜春 吴亚楠 谭梦岩

(1.西安科技大学测绘科学与技术学院， 西安 710054； 2.北京市农林科学院信息技术研究中心， 北京 100097；3.国家农业信息化工程技术研究中心， 北京 100097； 4.肇东市农业技术推广中心， 肇东 151100)

0 引言

我国工业化进程的快速发展，导致土壤重金属在污染水平、污染范围等方面面临日益严重的趋势。为了应对这种趋势，国家采取相关措施，如《土壤污染防治行动计划》(土十条)等开展土壤污染调查，农用地分类管理，加强农业污染源监管。不同时期土壤整体污染情况可以基于长期定位监测样本点进行预测与评估[1-3]，而样本点空间分布的均匀性和代表性是土壤整体污染情况监测准确性的重要指标。目前野外作业采集的样本点数据大多数存在样本点均匀性和代表性差的问题，导致样本点全局均值估计数据有偏，因此有针对性地提出样本点权重调整方法进行相应样本点数据纠偏非常关键。

空间采样是监测土壤整体污染情况及空间信息变化的重要手段，而统计推断是采样点布局优化的目的。样本点全局均值估计是统计推断的重要组成部分，其目的是准确估计研究区内某一数据的整体质量水平[4]。通过土壤样本点重金属全局均值估计与土壤重金属背景值进行对比分析，可以判断土壤是否受到污染[5]及污染程度[6-7]。然而，样本点全局均值估计的准确性受样本点密度及分布合理性影响，样本点密度越大且分布越均匀其全局均值估计越准确[8]。由于受研究区范围、抽样成本、地形等因素限制，难以获得高密度代表性样本点[9-12]，样本点空间分布易产生地理空间分布不均匀和特征空间代表性差的问题[13-14]，导致样本点全局均值估计值与全局均值真实值之间产生偏差。样本点数据纠偏主要包括样本点增删和样本点权重调整。样本点增删主要基于样本点均匀性检测结果，对密集分布样本点进行删除处理，对稀疏分布样本点进行加密处理[15-16]。样本点增删可以有效提高样本点空间分布的均匀性和代表性，但由于样本点加密过程中需考虑加密样本点的土壤属性值，因此对加密样本点与原始样本点之间的时间间隔有很高的要求。样本点权重调整不需要加密或删除样本点，只需调整样本点属性值权重，实现对样本点数据纠偏。样本点权重调整方法主要包括层次分析法[17-18]、专家打分法、格网法[19]、自组织双分区法[20-21]以及分层权重调整法[9]等。综合分析现有样本点权重调整方法后发现主要存在两大问题：①格网法、自组织双分区法直接对全部样本点进行权重调整。这种样本点权重调整模式无法确保大部分样本点属性值的真实值，降低样本点全局均值估计的可信度和准确性。②层次分析法、专家打分法所构建的样本点权重调整规则通常是静态的，即样本点权重不随评价因子状态变化而变化[22]，这将造成样本点在全局均值估计过程中产生偏差。综上所述，当前研究主要集中在对全部样本点通过统一的权重调整规则进行权重调整，而通过划分样本点类型确定待调整样本点数量并根据待调整样本点类型构建相应权重调整规则的方法亟待研究。

样本点权重调整的数据纠偏效果评价通常采用研究区大数据量的样本点某一属性值数据模拟估计的参考值作为真实值，或借助可获得真实值的属性值数据进行评价。研究区采样点数量越多，估计的参考值越准确[2,23-24]；由于受采样成本的限制，野外作业只能获取有限的样本点，难以估计获取准确的样本点土壤重金属全局均值，可以借助条件模拟方法进行模拟估计。归一化植被指数(NDVI)全局均值真实值通过研究区内遥感影像计算获取，也可以作为样本点权重调整的数据纠偏效果评价的参数。因此，本文以北京市顺义区农业用地为例，构建样本点类型划分规则以实现样本点类型划分，对非均匀样本点进行变权调整，基于土壤重金属Hg含量和NDVI 2个参数，分别设置与原始样本点、常规面积比例权重调整法的对比实验，并利用全局均值估计相对误差定量评价样本点权重调整的数据纠偏效果，旨在改善样本点空间分布的均匀性和代表性，进而提高样本点统计推断精度。

1 研究区概况与数据来源

1.1 研究区概况

顺义区位于北京市东北部，地理坐标为东经116°28′～116°58′，北纬40°00′～40°18′，境域东西长45 km，南北宽30 km，总面积为1 020 km2。其中农业用地面积为494.34 km2，是北京市重要的粮食和蔬菜生产基地。顺义区气候属暖温带半湿润大陆性季风性气候，年平均气温为11.5℃，年均降雨量约625 mm。土壤为河流洪水携带沉积物质造成，表面堆积物主要是砂土、亚砂土，北部山地最高点海拔为637 m，境内最低点海拔为24 m。选取高丽营镇、赵全营镇、牛栏山镇和北石槽镇的农业用地为研究区域，其中农用地主要是菜地、水浇地、苗圃等，农业用地面积为114.379 km2。由于该区域农业用地所占面积较大，对该区域农业用地土壤重金属污染状况进行监测有利于更好地了解顺义区土壤质量和农产品质量。

1.2 数据来源及预处理

采用2017年北京市顺义区样本点数据和NDVI数据。

(1)样本点数据：样本点数据来源于北京市农林科学院农产品质量安全管理平台。该平台样本点数据采集工作于2007年开展，样本点的初始数量为125个，并且根据土地利用类型及面积的变化，定期对样本点进行增删处理。本研究样本点数据采集于2017年，样本点布局和数量根据田块的土地利用方式和面积确定，采样时利用GPS获得样本点的地理坐标。原始采样中，每个样本点在10 m×10 m正方形网格4个顶点和中心点处选择0～20 cm耕层土壤混合，按四分法选取分析样品1.0 kg，最终得到66个农业用地样本点。所有土样在室内自然风干，碾压磨碎后，过100目尼龙网筛，按照国家标准分析测定土壤重金属Hg含量。

(2)NDVI数据：基于欧洲空间局空间分辨率为10 m的Sentinel-2卫星数据计算获取。影像获取时间为2017年9月5日。首先利用欧洲空间局发布的针对哨兵卫星数据的Sen2Cor软件对其进行辐射定标及大气校正，然后利用红波段(Sentinel-2A的第4波段)和近红外波段(Sentinel-2A的第8波段)的反射率数据计算顺义区NDVI[25]，并利用ENVI快速统计工具得到顺义区NDVI全局均值。

研究区内的样本点主要是作为样本点权重调整的实验数据，空间分布如图1所示。由于研究区内所布设的66个样本点无法准确估计土壤重金属全局均值的真实值，因此选择土壤重金属Hg的序贯高斯条件模拟全局均值与NDVI全局均值作为参数，分别用于定量评价样本点权重调整的数据纠偏效果。土壤重金属Hg含量侧重反映难以获取全局均值真实值的数据情景，NDVI侧重反映可以获取全局均值真实值的数据情景。2个参数的设置能够更好地反映样本点权重调整效果评价的客观性和科学性。

图1 农业用地样本点空间分布Fig.1 Spatial distribution of sampling sites of agricultural land

基于ArcGIS 10.5软件实现遥感影像裁剪、泰森多边形计算和序贯高斯条件模拟；基于ENVI 5.3软件实现NDVI提取及相应计算。

2 研究方法

技术流程如图2所示，3种权重调整方法包括泰森多边形权重调整法、原始样本点(无权重调整)、常规面积比例权重调整法。基于土壤重金属Hg含量和NDVI 2个参数，分别通过样本点全局均值估计相对误差对不同样本点权重调整的数据纠偏效果进行评价。其中常规面积比例权重调整法没有划分样本点类型，主要是将样本点所在泰森多边形面积与样本点泰森多边形平均面积的比值作为该样本点的调整权重，实现样本点权重调整。

图2 技术路线图Fig.2 Technical process

泰森多边形权重调整法主要步骤分为：①样本点类型划分：构建样本点泰森多边形，根据样本点所在泰森多边形面积确定面积置信区间，实现样本点类型划分及不同类型样本点空间位置确定。②样本点权重调整：分别构建不同类型样本点对应的权重调整规则，实现研究区内非均匀样本点变权调整。③样本点权重调整的数据纠偏效果评价：通过全局均值估计相对误差对样本点权重调整的数据纠偏效果进行评价。

2.1 样本点类型划分规则

样本点类型划分是确定非均匀样本点位置及数量的关键，合理的样本点类型划分规则可以提高样本点全局均值估计的准确性[25]。通过面积置信区间来划分样本点的类型，将泰森多边形面积小于面积置信区间下限的样本点划分为聚集样本点，将泰森多边形面积大于面积置信区间上限的样本点划分为稀疏样本点，将泰森多边形面积在面积置信区间之间的样本点划分为均匀样本点。

(1)通过样本点平均采样面积划分样本点所在泰森多边形类型。首先构建样本点泰森多边形，将样本点泰森多边形总面积与总采样点数的比值确定为样本点平均采样面积S0；将样本点所在泰森多边形面积按从小到大顺序排列，当样本点所在泰森多边形面积小于S0时，则划分为Ⅰ类泰森多边形(V1)，反之则划分为Ⅱ类泰森多边形(V2)。

(2)根据划分的两种类型泰森多边形面积确定面积置信区间。计算V1的面积均值μ1及面积标准差σ1来确定面积置信区间下限，计算V2的面积均值μ2及面积标准差σ2来确定面积置信区间上限，面积置信区间计算式为

(1)

2.2 样本点权重调整规则

样本点权重调整规则基于划分样本点类型确定，聚集样本点减小权重，稀疏样本点增大权重，均匀样本点无需进行权重调整。

(1)聚集样本点权重调整

在特定区域范围内，由于聚集样本点比理想均匀分布的样本点数量多，应该减小单个样本点的权重。聚集样本点权重调整公式为

Wi=1-ΔWi

(2)

(3)

式中 ΔWi——第i个聚集样本点权重调整量

Si——第i个聚集样本点所在泰森多边形面积，km2

(2)稀疏样本点权重调整

在特定区域范围内，由于稀疏样本点比理想均匀分布的样本点数量少，应该增大单个样本点的权重。稀疏样本点权重调整公式为

Wj=1+ΔWj

(4)

(5)

式中 ΔWj——第j个稀疏样本点权重调整量

Sj——第j个稀疏样本点所在泰森多边形面积，km2

(3)均匀样本点权重调整

均匀样本点不影响全局均值估计的准确性，则无需对这类样本点进行权重调整，对应的调整权重为1。

2.3 样本点权重调整的数据纠偏效果评价

利用全局均值估计相对误差来对样本点权重调整的数据纠偏效果进行评价。全局均值估计相对误差表示样本点全局均值估计值与全局均值真实值之间的偏离程度，计算公式为

(6)

式中Eg——全局均值估计相对误差，%

Mr——全局均值真实值

Me——全局均值估计值

当计算得到的全局均值估计相对误差越大时，全局均值估计值与全局均值真实值越偏离；当计算得到的全局均值估计相对误差越小时，全局均值估计值与全局均值真实值越接近。

2.4 序贯高斯条件模拟

序贯高斯条件模拟(Sequential Gaussian conditional simulation, SGCS)是一种对连续变量进行随机模拟的地统计学方法，该模拟方法结果在现有条件样本点处的模拟值为该点实测的真实值[26]。SGCS将现有数据作为条件数据计算待模拟点值的条件概率分布，从该分布中随机抽取一个值作为模拟值，连同原始数据、模拟数据一起作为条件数据代入下一点的模拟。该模型能产生多个可能的空间分布格局，有效避免普通克里格插值结果的平滑效应[27]。SGCS将研究数据作为一个整体，最大程度模拟研究区土壤重金属真实空间分布，模拟的结果与实测数据的直方图相同，并且模拟结果与实测结果具有相同的空间自相关关系，因此SGCS可以实现研究区土壤重金属含量模拟。

3 结果与分析

3.1 样本点类型划分结果

以北京市顺义区农业用地样本点为例，构建农业用地样本点泰森多边形并计算泰森多边形面积，进行农业用地样本点类型划分，结果如图3所示。图中样本点所在泰森多边形面积按照从小到大顺序排列。通过样本点划分规则判断存在2个聚集样本点，所在泰森多边形面积分别为0.265、0.423 km2；1个稀疏样本点，所在泰森多边形面积为3.439 km2；63个均匀样本点，所在泰森多边形面积在0.653～3.109 km2之间。不同类型样本点的空间分布如图4所示，图中2个聚集样本点分别是聚集样本点A和聚集样本点B。

图3 农业用地样本点类型划分结果Fig.3 Classification of sampling sites of agricultural land

图4 不同类型样本点的空间分布Fig.4 Spatial distribution of different types of sampling sites

3.2 样本点权重调整结果

利用泰森多边形权重调整法实现非均匀样本点权重调整，通过式(3)计算得到聚集样本点A和聚集样本点B权重调整量分别为0.391和0.117，通过式(2)计算得到对应的调整权重分别为0.609和0.883；通过式(4)、(5)分别计算得到稀疏样本点权重调整量为0.068，调整权重为1.068；对于所有均匀样本点不需进行权重调整，调整权重均为1。利用常规面积比例权重调整法计算得到全部样本点的调整权重，其对应范围为0.176～2.289。不同权重调整方法的样本点权重调整结果如表1所示。

表1 不同权重调整方法的样本点权重调整结果Tab.1 Weight adjustment results of sampling sites for different weight adjustment methods

3.3 数据纠偏效果评价与对比分析

基于土壤重金属Hg含量和NDVI 2个参数，分别设置与原始样本点及常规面积比例权重调整法的对比实验，计算不同方法对于不同评价参数的样本点全局均值估计相对误差，定量评估泰森多边形权重调整法数据纠偏效果。

3.3.1土壤重金属Hg全局均值相对误差

基于泰森多边形权重调整法、原始样本点和常规面积比例权重调整法计算样本点的调整权重，顺义区农业用地土壤重金属Hg含量的全局均值估计结果分别为0.069 9、0.070 9、0.068 1 mg/kg。基于原始土壤重金属Hg含量数据序贯高斯条件模拟得到顺义区农用地土壤重金属Hg含量全局均值参考值。首先，基于ArcGIS 10.5软件的探索数据工具探索最优变换方式及最优参数设置，通过地统计向导工具得到简单克里金插值结果；其次，基于插值结果进行序贯高斯条件模拟，为了避免由于模拟随机性带来的模拟误差，实验分别进行了10、15、20、25、30、35、40、45、50、55、60等不同次数模拟，直至模拟均值结果趋于平稳，获取最适的条件模拟次数为50；最后，对样本点进行50次序贯高斯条件模拟得到土壤重金属Hg含量全局均值为0.070 2 mg/kg。通过式(6)计算得到泰森多边形权重调整法、原始样本点以及常规面积比例权重法的土壤重金属Hg含量的全局均值相对误差分别为0.413%、1.056%和2.933%，如表2所示。结果表明，泰森多边形权重调整法具有较好的数据纠偏效果，可以有效地提高土壤重金属Hg含量全局均值估计结果的准确性。

表2 不同方法样本点重金属Hg含量全局均值估计相对误差Tab.2 Relative errors of Hg mean estimation of sampling sites for different methods

3.3.2农业用地NDVI全局均值相对误差

基于泰森多边形权重调整法、原始样本点和常规面积比例权重调整法计算样本点的调整权重，顺义区农业用地NDVI的样本点全局均值估计结果分别为0.691、0.697和0.700。基于研究区农用地全部996 148个像元计算获得NDVI全局均值为0.680，以此作为全局均值真实值。通过式(6)计算得到泰森多边形法、原始样本点和常规面积比例权重法的NDVI全局均值相对误差分别为1.617%、2.500%及2.941%，如表3所示。结果表明，与原始样本点和常规面积比例权重法相比，泰森多边形权重调整法NDVI全局均值相对误差最小，说明该方法对样本点的数据纠偏效果最好。

表3 不同方法样本点NDVI全局均值估计相对误差Tab.3 Relative errors of NDVI mean estimation of sampling sites for different methods

由表2、3可以看出，泰森多边形权重调整法对于不同类型的评价参数均能有效地进行数据纠偏，可以提高样本点在全局均值估计的代表性，使其统计推断结果更加准确。分析其原因可能包括：①泰森多边形权重调整法仅对非均匀样本点进行权重调整，能确保大部分样本点属性值的真实性，使样本点全局均值估计结果更加可靠。②泰森多边形权重调整法构建不同权重调整规则可实现对非均匀样本点变权调整，针对性和目标性更强，调整后的样本点更具代表性，进而提高样本点全局均值估计的准确性。

4 讨论

泰森多边形权重调整法是一种面向全局均值统计推断的样本点数据纠偏方法，通过对样本点类型划分来确定不同类型样本点数量及其空间位置，并基于样本点所在泰森多边形面积构建不同类型样本点权重调整规则及其权重调整量，进而实现对聚集样本点和稀疏样本点的变权调整。该方法能够确保大部分样本点属性值的真实性，并且对非均匀样本点变权调整，提高样本点的代表性，最终获得更加准确、可靠的全局均值估计结果。与常规样本点权重调整方法相比，设置的土壤重金属Hg含量及NDVI数据纠偏效果评价对比实验，均可验证泰森多边形权重调整法对全局均值估计的样本点数据纠偏效果更好。面向全局均值估计的样本点权重调整方法包括两个具体的应用场景：①用于国家及各级地方政府监测管理所属管辖行政单元的土壤重金属整体污染评价，辅助宏观政策决策与管理。②应用于某个行政单元内不同空间分区的土壤重金属污染评价，辅助区域土壤重金属污染治理与防治。然而，该方法在研究过程中存在一定的不确定性。

(1)不同类型样本点数量划分存在不确定性。不同类型样本点数量通过面积置信区间确定，在不同置信水平下，面积置信区间的范围会随之发生变化，导致不同类型样本点的数量也会随之改变，权重调整样本点数量较多时会对样本点数据纠偏效果产生影响。

(2)尽管序贯高斯条件模拟方法能够获取比较准确的土壤重金属含量模拟参考值[22-23]，但是参考值不一定等于真实值，未来可以通过高密度代表性采样来进一步模拟逼近真实值的土壤重金属元素全局均值。

(3)泰森多边形权重调整法调整样本点并未改变样本点空间位置，而是通过改变样本点某个属性值的权重实现对样本点的数据纠偏，这仅仅是通过对样本点特征空间来进行调整，但并未对样本点的地理空间进行调整，未来将进一步研究耦合地理空间和特征空间的样本点权重调整方法，以兼顾样本点在地理空间的均匀性和特征空间的代表性。

5 结论

(1)泰森多边形权重调整法划分得到2个聚集样本点、1个稀疏样本点和63个均匀样本点。2个聚集样本点调整权重分别为0.609和0.883，稀疏样本点调整权重为1.068，均匀样本点调整权重均为1。

(2)泰森多边形权重调整法、原始样本点、常规面积比例权重调整法进行顺义区农业用地土壤重金属Hg含量和NDVI的全局均值估计相对误差分别为0.413%、1.056%、2.933%和1.617%、2.500%、2.941%。泰森多边形权重调整法对于不同评价参数的全局均值估计相对误差均为最小，与全局均值真实值更为接近，表明泰森多边形权重调整法数据纠偏效果更好。

(3)泰森多边形权重调整法准确定位到非均匀样本点位置并实现非均匀样本点的变权调整，可以提高样本点的代表性，更加准确地进行全局均值估计。但该方法仅在特征空间对样本点进行权重调整，且在面积置信区间确定、序贯高斯条件模拟土壤重金属含量等方面存在一定的不确定性。