1 引 言

土壤安全是农产品安全的重要基础，其直接影响着国民经济健康发展。当前我国土壤质量在工农业生产扰动作用下退化严重，开展土壤质量监测、掌握土壤质量变化规律势在必行。空间抽样是构建与优化土壤质量监测网络的关键技术^[1]。与传统全面调查方法相比，空间抽样技术能够通过布设较少的样点及时准确地获取抽样变量的动态变化信息，从而避免全面调查造成的巨大人力、物力和财力开销^{[2, 3, 4]}。

土壤空间抽样方法主要包括专家布样法、基于设计的经典统计抽样法和基于模型的抽样方法3类，其中专家布样法属于非概率抽样，依赖于专家主观知识，其结果在指定置信水平下不具有可比性，容易产生偏误^[5]；经典统计抽样包括随机抽样、系统抽样和分层抽样等方法，适用于静态调查(here and now)或者总量估计^{[2, 4]}；基于模型的抽样方法主要在地统计学理论方法基础上发展而来，适用于处理未知点预测、土壤制图和时序抽样设计，因此在土壤空间抽样与监测网络设计领域应用广泛^[6]。地统计背景下的模型抽样方法通常以减小协方差(变异函数)参数估计(estimation)误差和未知点插值预测(prediction)误差为准则，采用设计选择算法获取满足目标要求的最优抽样方案^[6]。此类方法通过在设计过程中引入抽样精度、成本、代表性等多目标以及可达性、空间自相关性和空间均质性等约束条件，能够显著提高抽样方案的精度与代表性^[7]。但是，上述机制同时也增加了优化问题的复杂程度，导致四叉树、序贯法等传统精确设计选择方法效率受限^{[8, 9]}。有鉴于此，一些学者将抽样设计映射为组合优化问题，广泛采用模拟退火、遗传算法等离散随机进化算法进行问题求解，并取得了显著效果^{[10, 11]}。

粒子群算法(particle swarm optimization,PSO)属于群智能进化算法，由粒子间相互协作表现出的智能行为控制寻优过程，具有搜索速度快、结构简单、易于实现的特点^{[12, 13]}，其应用范围广泛，优化性能相对高于模拟退火等离散优化算法^{[13, 14, 15]}。针对土壤抽样方案的变异函数参数估计精度和未知点预测精度需求，即样本方案需兼顾微观上拟合抽样变量变异规律的准确性和宏观上空间分布的均衡性^[16]，本文以最小克里金方差准则和极大熵准则为抽样目标，以空间可达性、最优样本量、最小成本和最优样点间距为约束条件进行样点布局优化。为了实现上述多目标优化以及样点布局微观调整等操作，本文结合多目标粒子群和邻域粒子群算法^{[17, 18]} 提出了多目标微观邻域粒子群方法(multi-objective micro-neighborhood PSO,MM-PSO)，并选择试验区进行方法有效性验证。

MM-PSO模型目的在于协同优化抽样方案的协方差函数(变异函数)参数拟合精度与未知点预测精度。对于网格化研究区域而言，抽样设计即是采用MM-PSO模型确定区域内抽样单元是否被选择为样点的0-1组合优化过程，其中样点的选取依赖于粒子速度更新策略和本文设置的粒子微观邻域操作策略。基于上述原理，MM-PSO模型主要包括3个关键步骤：① 定义粒子的空间化表征方式，建立粒子群与抽样问题的映射关系；② 结合土壤空间抽样特征，设计粒子群算法的多目标适应度函数和粒子的微观邻域操作策略；③ 应用多目标微观邻域粒子群算法寻找最优抽样方案。

存在离散化土壤抽样区域A=(s₁,s₂,…,s_N)，s_n∈A为网格单元，N为格网单元数量，MM-PSO模型主要通过建立二进制粒子及其维度与抽样区域及网格单元的映射关系实现粒子的空间化表征。设S表示种群规模为P的二进制粒子群，S^p(t)为粒子S^p在t时刻的位置，记为N维向量S^p(t)=[s₁^p(t) s₂^p(t) … s_N^p(t)]，则S^p(t)代表t时刻抽样区域的离散状态，s_n^p(t)=1表示格网单元i被选为样点，s_n^p(t)=0表示非样点，s_n^p(t)=－1表示采样不可达区域。相应的，t时刻粒子S_n^p代表的抽样方案可用集合表示为{s^p(t)|s_n^p(t)=1,n=1,2,…，N}。可见，粒子在土壤抽样目标和约束限制下，通过不断改变其维度值的组合方式，即可逐渐寻找到最优抽样方案。图 1描述了一个具有10个抽样单元的地理空间与粒子的映射关系，该粒子编码为(1001010101)，其对应抽样方案中包含5个样点。

图 1 抽样区域,粒子与抽样方案之间的映射关系 Fig. 1 Mapping relationships between the sampling area, the particle and the sampling network

图选项

MM-PSO算法的多目标适应度函数用于计算粒子适应度值并评价粒子的优化程度，也即用于评判粒子对应的抽样方案是否满足预定抽样目标。适应度函数由最小克里金方差和极大熵准则构成。

误差大小是评价土壤空间抽样方案精度的主要指标，最小克里金方差准则能够使抽样点的分布处于合适的位置以使克里金插值后生成的表面精度最高或由已知抽样点对所有未抽样点的估计精度最高。在已知土壤变量的协方差函数的情况下，克里金方差主要取决于样点位置信息，而与样点具体取值无关，因此成为目前应用最为广泛的抽样设计准则之一^[16]。设粒子S^p(t)中任意一个未采样维度(抽样单元)为s_x^p(t)，则该点克里金方差取决于抽样变量的先验方差与单元s_x^p(t)到其邻域范围内预采样样点协方差系数的差值，具体公式如下

点击浏览公式

式中，Cov_z(s_i^p(t),s_j^p(t))为样点s_i^p(t)和其邻域内样点s_j^p(t)的协方差，则λ_j与μ由下式计算得到

点击浏览公式

式中，s_i^p(t),s_j^p(t)=1;Cov_z(0)为抽样变量的先验方差;λ_j为克里格差值权重系数;μ为拉格朗日系数。

根据信息熵理论，在已知部分知识的前提下，关于未知分布最合理的推断就是符合已知知识最不确定或最随机的推断，亦即熵值最大的情况下作出的选择不包含任何主观假设与约束。土壤空间抽样总体具有随机性，因此抽样方案对应的信息系统熵值越大，越能客观地揭示土壤变量的总体空间变异程度。该准则对于多目标抽样尤为实用，它能够使抽样方案尽可能包含每个抽样目标的相关信息^[19]。极大熵准则用于衡量粒子S^p(t)熵值的数学表达式如下

点击浏览公式

式中，z(s_i^p(t))为粒子Sp(t)维度i对应样点的观测值；p为该观测值在抽样变量总体中出现的概率。</sup>

模型采用动态加权求和法(dynamic weighted aggregation,DWA)将上述准则进行整合，具体如式所示

点击浏览公式

式中，权重系数w₁和w₂在一定循环周期R内动态变化(式(4))，从而避免了传统静态权重过于主观的不足^[20]。

点击浏览公式

设V^p(t)=(v₁^p(t),v₂^p(t),…,v_N^p(t))为粒子S^p在t时刻的飞行速度，sb^p(t)为粒子S^p在t时刻所经历过的个体最好位置，Sg(t)为种群S在t时刻的全局最好位置，则粒子t+1时刻位置受限于粒子t时刻位置、个体最优位置、种群全局最优位置以及微观邻域规则，具体公式如下所示

点击浏览公式

式中，r₁、r₂、r∈[0,1]为随机生成数；c₁、c₂分别为个体信息参数和社会信息参数；Sig()为Sigmoid函数，对t+1时刻的速度v_n^p(t+1)进行0-1转换；I[]为取值函数，当括号内条件满足时值为1，否则为0；c_i(S^p(t))为粒子微观邻域操作规则取值。

(1) 样点可达性是影响样本代表性的重要因素。当土壤抽样单元位于建筑物或者水域范围内，或者样点位于坡度较大的山地上时，采样人员通常无法实地获取所需数据，则该备选单元无法成为样点，具体可表示为

点击浏览公式

(2) 在顾及样点间关联情况时，所需样点数量要小于样点独立情况时的样本量N_max，则样本量约束取值范围如下式所示

点击浏览公式

式中

点击浏览公式

式中，z_a/2表示置信度；σ表示样点观测值的标准差；μ代表样点观测值的均值；p代表空间抽样的相对误差。

(3) 采样成本限制条件是指根据《测绘生产成本费用定额》将抽样区域分为Ⅰ级困难区域(平原地区为主)、Ⅱ级困难区域(丘陵地区为主)和Ⅲ级困难区域(山地地区为主)，设总采样费用为FT，基础采样费用为FB，Ⅰ级困难区域样点采样费用为FS₁，样点数量为N₁，Ⅱ级每个样点采样费用为FS₂，样点数量为N₂，Ⅲ级区域为FS₃，样点数量为N₃，则存在

点击浏览公式

(4) 样点间距过大容易导致样点布局无法准确描述地理空间的变化趋势，而间距过小又将造成采样经费超支。最小空间相关性约束旨在限制样点间距离，使得样点能够相对均衡地分布于抽样区域内。根据文献^[21]，设α为土壤变量变异函数描述的变程，则合理的样点间距d应满足

点击浏览公式

MM-PSO优化抽样方案的流程如下：

(1) 随机初始化粒子群S，根据式(9)确定粒子包含的初始样点数量不大于N_max，样点位置随机确定。

(2) 根据循环周期R和公式(5)计算优化目标权重；根据公式(4)计算每个粒子S^p(t)的优化目标函数值(粒子适应度)，并通过比较粒子适应度确定单个粒子的历史最佳位置sb^p(t)和整个种群的全局最优位置Sg(t)。

(3) 采用公式(6)更新粒子的位置，粒子维度值变化受4类粒子微观邻域规则约束。更新后，当满足粒子微观领域规则的样点数量少于初始样点数量时，则随机补全样点，从而维持样点数量不变并保证粒子群种群多样性，此随机补全样点过程同样受粒子微观邻域规则约束。

(4) 判断算法是否满足终止条件，如果满足则进入步骤(5)，否则返回步骤(2)。算法终止条件包括：①当在一定迭代次数T⁰内适应值的变化阈值小于某一设定值ρ时算法终止；②最大迭代次数T_max。

(5) 输出全局最优方案Sg(t)为最优抽样方案。

试验区横山县位于陕西省北部，毛乌素沙漠南缘(108°65′E～110°02′E、37°22′N～38°74′N)，属于陕北黄土高原与风沙高原的过渡区(如图 2)。历史悠久的农业开垦和长期的土地不合理利用行为，使得自然植被破坏殆尽，加之相对干旱的气候和疏松的土质，导致土壤在水蚀和风蚀作用发生严重退化，从而对当地粮食安全产生影响，一定程度阻碍了社会经济发展。

图 2 研究区域位置图 Fig. 2 The location of the study area

图选项

收集了试验区1∶5万土地利用现状数据、坡度数据以及土壤质量预采样数据。土地利用现状数据与坡度数据作为样点空间布设的约束条件。土壤质量预采样数据来源于2004年野外实地调查，共计251个样点，主要作为先验知识揭示区域土壤指标空间变异结构。本试验拟对土壤有机质含量(soil organic matter,SOM)的抽样方案进行优化，从而为该区域土壤肥力变化、土壤碳效应等研究提供技术支撑。

为了满足空间分析的要求，试验采用ArcGIS10软件对数据进行LOG变换，使得原始数据服从偏正态分布。以稳态协方差理论模型拟合土壤有机质空间变异规律，得到块金值(nugget)为0.029 3，基台值(sill)为0.278 8，变程α为2 683.95 m。采用普通克里金法和序贯高斯随机法模拟土壤有机质总体100次，用于提取所设计抽样方案中各点的估计值。

MM-PSO算法参数参照文献^[12]及100次独立重复试验结果综合确定：种群规模为50，社会和个体参数c₁=2.25、c₂=2.05，惯性权重w初始值设为0.9并逐渐匀速降为0.6，算法终止条件T⁰=10、ρ=0.001、T_max=600，权重系数循环周期R=5。微观邻域算子参数设定如下：采样可达性约束中坡度上限为60°，采样成本限制中Ⅰ级困难区域坡度范围[0°,5°]、Ⅱ级困难区域坡度范围(5°,25°］、Ⅲ级困难区域坡度范围(60°,90°］，95%置信度、3%相对精度下的最大样本量N_max=421。试验对比分析了不同目标权重方案(MKV/ME/MKV+ME)、不同样本量(100/200/300)及不同算法(多目标空间模拟退火算法MO-SSA与MM-PSO)设计的抽样方案。为了确保可比性，不同方案所涉及的粒子群算法参数保持不变，MO-SSA算法采用与MM-PSO相同的抽样目标和微观邻域操作算子，并设定关键参数初始温度、终止温度和扰动系数分别为2.0、0.1和0.9。试验环境为单机PC，处理器为Intel E5800、内存2 G，算法采用C#.NET2005+GDAL组件库编程实现。

图 3分别给出了不同目标权重方案所获取的最优抽样方案。从不同样本量方案的样点空间分布情况看，MKV+ME多目标生成的样点方案相对合理：① 从整个区域来看，MKV方案样点分布偏向区域中心，ME方案样点分布偏向与区域边界，而MKV+ME方案样点空间分布相对比较均衡，合理兼顾了边界和区域内部的样点密度^{[10, 22]}；② 在区域内部，单一目标对应的样点配置方案存在较多间隙，而MM-PSO在微观邻域规则的作用下样点分布更加均衡和合理。统计分析中均值、标准差、偏度和峰度是描述统计变量空间分布情况的基本统计量，可以衡量不同抽样方案与预采样样点属性空间分布模式的相似程度^[16]。表 1结果显示MM-PSO获取的样本与预采样数据的分布模式更为接近，表明在多抽样目标情况下获取的抽样方案更加具有代表性。

图 3 不同权重方案与样本量条件下土壤有机质含量最优抽样方案 Fig. 3 Optimal sampling networks of SOM for 100,200 and 300 samples by different weighted scenarios

图选项

表 1中，假设均值x，标准差s，偏度skewness=，衡量SOM 采样数据分布非对称程度；峰度kurtosis=，表征概率密度分布曲线在平均值处峰值高低的特征数。

为了进一步验证上述结论，图 4对比分析了不同方案的抽样目标适应度值与均方根误差(RMSE)，均方根误差由预采样样点估计值与观测值计算得到^[24]，其中所有适应度值均被归一化到^{[0, 1]}。结果显示随着样本量的增大，抽样方案目标适应度与精度不断提升，提升幅度则略有下降，说明MM-PSO能够准确揭示抽样方案代表性与样本量之间的幂函数正比关系。在适应度提升方面，MKV+ME多目标所获取优化方案的MKV_A目标分量和ME_A目标分量稍劣于MKV和ME两个单目标结果，一定程度说明MM-PSO算法的多目标寻优能力源于MKV和ME单目标适应度的适当降低，从而寻求达到二者之间的帕累托(Pareto)平衡，表明了MM-PSO具有较强的多目标协同寻优能力。从抽样精度(RMSE)来看，多目标(MKV+ME)的精度明显高于单目标方案，与多目标样点空间布局相对均衡的结论保持一致。

图 4 不同权重、样本量时SOM抽样方案的目标适应度与均方根误差 Fig. 4 Fitness values and RMSE of the sampling networks of SOM

图选项

收敛性能与稳定性是衡量算法效率的关键指标。试验从最优收敛值、收敛代数、收敛时间和收敛曲线4个方面对比分析了不同算法的运行效率。由表 2可知，MM-PSO算法明显具有更强的收敛能力，其收敛代数和时间仅是MO-SSA算法的3.8%和43.9%。从收敛过程来看，MM-PSO的收敛过程均呈现了前段陡降、后端平缓的趋势，收敛曲线均比较平滑，表明粒子群具有较强的全局和局部搜索能力和收敛性能^[23]，而MO-SSA收敛曲线前段坡度平缓，收敛过程波动较大，表明其在本问题中收敛效果较差(图 5)。

图 5 MM-PSO与MO-SSA算法收敛曲线对比 Fig. 5 The convergence curves of the MM-PSO and the spatial simulated annealing algorithm in 100 times of simulations

图选项

本文提出的多目标微观邻域粒子群算法(MM-PSO)，设计了粒子的多目标适应度函数和微观邻域操作策略，用于高效协调土壤空间抽样过程中的多目标冲突，并完成了陕西省衡山县土壤有机质空间抽样方案设计。试验结果表明：① MM-PSO算法综合考虑了土壤抽样最小克里金方差和最大信息熵目标以及样点数量、布局等微观准则，对多抽样目标的良好协同能力使其能够获得较优的收敛值，样点空间分布更加均衡，统计特性更加接近观测值特征；② 相比SPSO和MO-SSA算法，MM-PSO稳定性、收敛能力和寻优效率更高。模型的应用能够较好地协调土壤空间抽样过程的拟合精度与插值预测精度之间的矛盾，大幅度提升抽样设计的效率，对设计与构建中、大尺度的土壤空间优化抽样方案与土壤环境监测网络具有重要意义。此外，由于地理要素空间抽样技术框架具有一致性，本研究提出的多目标微观邻域粒子群抽样方法亦可扩展用于支持大气、水体和植被等地理要素抽样网络的空间优化决策，从而为空间抽样模型的研究提供了一种思路。

在实际应用中，模型仍有待改进和完善，例如探讨模型不确定性对抽样方案的影响，针对空间变异异质性总体进行抽样方案设计，探讨抽样先验知识对优化结果的影响作用等^[24]。

致谢：中国土地勘测规划院国土资源部土地利用重点实验室为本文提供了研究数据，在此表示衷心感谢。