一种随机分布型地理要素的多目标协同抽样方案设计方法

摘要

本发明涉及一种随机分布型地理要素的多目标协同抽样方案设计方法。1，选择多个具有随机分布特征的地理要素作为抽样指标，获取指标的预采样样点数据；2，利用普通克里金法探索预采样数据的空间变异结构，拟合生成各个指标的理论半变异函数，根据理论半变异函数和1中的预采样样点数据，采用高斯序贯随机模拟方法得到各个指标的空间分布总体；3，根据指标的精度要求和先验方差计算最小协同样本量，并设置协同微粒群算法参数；4，利用协同微粒群算法对现有预采样样点方案进行优化，得到多指标协同抽样方案。本发明能够节约实地抽样成本，提高了抽样样点布设过程的自组织性和优化性，提高了地理要素多指标协同抽样方案的设计效率。

说明书

技术领域

本发明涉及一种多目标协同抽样方案设计方法，尤其是涉及一种随机 分布型地理要素的多目标协同抽样方案设计方法。

背景技术

地理要素包括地形、气候、水文、生物和土壤等，对其实施准确的调 查是解释地理要素空间分布规律和演变过程的关键，也是资源管理、灾害 管理、生态环境治理以及全球变化研究的需要。在实际操作过程中，研究 人员为了节约抽样成本和提高抽样效率，通常仅设计一套抽样方案同时实 施多种地理要素指标的空间抽样调查任务，因此如何设计抽样方案保证多 指标的抽样精度要求成为学者们关注的焦点。

目前存在的两类空间抽样策略中，基于设计的抽样策略(Design based  sampling)适用于大尺度、大样本量的多抽样目标抽样，而在处理数据缺 失(Missing data)和数据分析等问题上存在不足；基于模型的抽样策略 (Model based sampling)通过构建总体的拟合模型能够避免数据缺失对 抽样结果的影响，在实际应用中较多关注单一指标。随机分布型自然地理 要素空间分布上同时具有随机性和空间结构性，因此适合应用地统计学相 关原理和方法探索其空间分布规律，构建空间分布的总体模型，辅助进行 空间抽样方案设计。可见，改进基于模型抽样策略，构建一种多要素协同 抽样的方法已经成为自然地理要素空间抽样方法的一种发展趋势。

从本质上看，多抽样目标协同抽样需要同时满足不同指标对抽样精度 的要求，因此属于多目标优化范畴内的问题。传统多目标协同方法将加权 平均法或最小阈值法与模拟退火求解算法进行结合，效果显著。加权平均 法(Weighted Average method，WA)的核心思想是采用加权求和的方式将 多目标问题转换为单目标问题，以多指标的联合克里金方差达到最小作为 最优抽样方案的判断依据。在加权求和时，为了保证不同指标的克里金方 差值具有相同的量纲，通过将克里金方差除以相应指标协半变异方差函数 的基台值方法方式对其进行标准化处理。最小阈值法(Minimum Threshold  Method，MT)的核心思想是为每个指标设定合适的最小克里金方差目标值， 并在每次样点位置变更时要保证所有指标的克里金方差低于或者尽可能接 近该目标值。如果某一指标的克里金方差等于其目标值表明当前样本方案 能够精确描述该指标的空间变异规律；而如果低于其目标值，则表明当前 样本方案存在一定的冗余。但是，上述传统方法存在一定不足，如对于复 杂的多目标问题加权平均法的权值难以获取，并且设定的权重通常会导致 某些目标达到最优而其他目标却较差的极端情况；最小阈值法要求所有目 标值均要达到预定的最优值，因此可能搜索不到最优解，并且对于一些目 标可能较难设定最优值；另外，已有更加高效的搜索算法可用于替代模拟 退火算法。

协同微粒群算法(Coorperative Particle Swarm Optimization，CPSO) 能够通过个体的相互影响和相互合作，自发地完成时间、空间和功能上有 序的群体行为，在多目标优化领域得到广泛的应用。该算法通过采用独立 与协同相结合的进化方式，大大地提高了算法的全局收敛能力。因此，引 入二进制协同微粒群算法进行多地理要素协同抽样具有较大的可操作性和 应用空间。

发明内容

本发明主要是解决现有技术所存在的技术问题；提供了一种能够节约 实地抽样成本，能提高地理要素多指标协同抽样方案设计效率，抽样样点 布设过程具有自组织性、优化性和智能性的一种随机分布型地理要素的多 目标协同抽样方案设计方法。

本发明的上述技术问题主要通过以下技术方案解决：

一种随机分布型地理要素的多目标协同抽样方案设计方法，其特征在 于，包括以下步骤：

步骤1，选择多个具有随机分布特征的地理要素作为抽样指标，获取指 标的预采样样点数据，将抽样区域进行格网划分，形成的网格单元作为抽 样单元，抽样区域内的抽样单元集合作为抽样框；

步骤2，利用普通克里金法探索预采样数据的空间变异结构，拟合生成 各个指标的理论半变异函数，根据理论半变异函数和步骤1中的预采样样 点数据，采用高斯序贯随机模拟方法得到各个指标的空间分布总体；

步骤3，根据指标的精度要求和先验方差计算最小协同样本量，并设置 协同微粒群算法参数；

步骤4，利用协同微粒群算法对现有预采样样点方案进行优化，得到多 指标协同抽样方案，所述现有预采样样点方案即步骤1所述的获取指标的 预采样样点数据。

本发明采用基于“环状模型”的二进制协同微粒群算法生成随机分布 型自然地理要素的多目标空间协同抽样方案，适用于在有效保持各个指标 样点的代表性和抽样精度的前提下通过构建一套抽样方案对多指标同时进 行空间抽样，能够节约实地抽样成本；以平均克里金方差和信息熵为微粒 适应度的衡量准则，采用最小协同样本量、抽样可达性和抽样费用作为样 点空间布设的条件约束，动态寻找可变样本量与可变样点布局的最优组合 方案，提高了抽样样点布设过程的自组织性和优化性；本发明方法继承了 协同微粒群算法的特性，自组织自学习能力强、智能化程度高，提高了自 然地理要素多指标协同抽样方案的设计效率。

在上述的一种随机分布型地理要素的多目标协同抽样方案设计方法， 所述的步骤3中，计算最小协同样本量的步骤如下：定义抽样样本空间A， 抽样单元数为N，抽样指标数量为V，设指标i存在一个最优的抽样设计，其 精度要求是p_i，对应的样本容量为n_i；

步骤3.1，设定指标1的精度要求p₁、总体均值μ₁、指标1抽样均值以及标准差s₁，采用传统随机抽样方法中的最小样本量计算公式确定指标1 的最小样本容量为n₁，形成的抽样样本为ξ₁；

其中，$n_1=\frac{{p_1}^2{s_1}^2}{{(\overline{x_{1n}}-{\mu }_1)}^2}$式一；

步骤3.2，依据指标2的抽样精度要求检验样本ξ₁，剔除不符合指标2 抽样设计要求的样点后，形成一组新的样本

步骤3.3，依据指标2的精度要求，推算出样本所能代表的最大子样 本空间A₂及其剩余样本空间

步骤3.4，在剩余样本空间中依据指标2的抽样设计要求采用随机 抽样方法确定容量为n₂的一个样本ξ₂，则容量为n₁+n₂的样本ξ₁∪ξ₂是同时满 足指标1和指标2抽样设计的一组最小样本；

步骤3.5，重复执行步骤3.1至步骤3.4，重复执行次数为V-1次，其 中，N₀＝n₁+n₂+…+n_v，并且ξ₁，ξ₂，Λ，ξ_v相互独立，则N₀就是满足多指标协 同抽样设计的最小样本容量。

在上述的一种随机分布型地理要素的多目标协同抽样方案设计方法，所 述的步骤4具体操作步骤如下：

步骤4.1，构建微粒与抽样区域的映射关系，依据映射关系初始化微粒 群S；

步骤4.2，根据指标数量V等分微粒群，当种群规模P被V整除时，则子 微粒群S_v的种群规模P_v＝P/V；当种群规模P不被V整除时，为了避免微粒群 微粒信息丢失，首先确定每个子群包含的基本微粒个数然后将剩余 的微粒依次分配给前个子群，对于未额外分配微粒的子群则从其 中随机选择一个微粒进行复制，保证每个子微粒群S_v的种群规模均为 其中为取整函数；

步骤4.3，构建子种群的环状协同结构，以一个随机的子种群为环状结 构的首端，采用随机选取的方式确定其邻域子种群，重复本步骤直到所有 子种群被组合进环状协同结构；

步骤4.4，在一个协同周期T_coor内，设定Flag标识子种群独立进化是否完成， t_coor标识协同周期内子群进化的代数，初始值为1，子微粒群S_v分别独立寻 找一个土壤属性的最佳抽样方案，其中，(1≤t_coor≤T_coor)；

步骤4.5，获取所有子群的当前全局历史最优位置作为最优解组成 Pareto解集，如果Flag＝true进入步骤4.6，否则将Pareto解集中的最优 解i分配给相邻子种群S_i+1作为S_i+1下一个协同周期的初始全局最优位置并返 回步骤4.4，当i＝V时，则将最优解i分配给子种群S₁；

步骤4.6，从Pareto解集中随机选择一个最优解作为随机分布型地理要 素的多目标协同抽样方案输出。

在上述的一种随机分布型地理要素的多目标协同抽样方案设计方法，所 述的步骤4.4中，独立进化步骤如下：

步骤4.4.1，计算S_v中每个微粒S_vi(t)的适应度；

步骤4.4.2，确定微粒个体历史最佳位置和全局历史最优位置；

步骤4.4.3，判断是否满足终止条件，如果满足则设定Flag＝true并执行 步骤4.5，同时将步骤4.4.1中生成的Sg_v(t)作为子群的最优抽样方案输出， 否则将进入步骤4.4.4；

步骤4.4.4，判断是否完成一个协同周期t_coor＝T_coor，如果完成则设定 Flag＝false跳至步骤4.5，并将步骤4.4.1中生成的Sg_v(t)作为子群的最优 抽样方案输出，否则进入步骤4.4.5；

步骤4.4.5，根据微粒当前个体历史最佳位置和当前全局历史最优位置 更新微粒的位置，t_coor自动加1，并返回步骤4.4.1。

在上述的一种随机分布型地理要素的多目标协同抽样方案设计方法，所 述的步骤3中，微粒群聚类算法的参数包括种群规模P、最大迭代次数、惯 性权重w、随机参数r₁和r₂、个体信息加速常数c₁、社会信息加速常数c₂、收 敛因子x和子种群协同周期T_coor，采用最大迭代次数作为算法终止条件。

在上述的一种随机分布型地理要素的多目标协同抽样方案设计方法， 所述的步骤4.1构建的映射关系中，一个二进制微粒为抽样单元的集合， 表示一种抽样方案，微粒中的每一个维度代表一个抽样单元，维度值为1 表示相应抽样单元被选为样点，维度值为0表示相应抽样单元未被选为样 点。

在上述的一种随机分布型地理要素的多目标协同抽样方案设计方法，所 述的步骤4.4.1中，微粒的适应度采用平均克里金方差和样点信息熵进行 衡量，确定微粒S_vi在t时刻个体历史最佳位置Sb_vi(t)与微粒群全局最佳位置 Sg_v(t)的步骤如下：

当t+1时刻微粒的适应度优于其历史最优的适应度时，则将微粒当前的 位置记录为个体历史最佳位置，否则保持不变；

当子群S_v在t+1时刻处于一个新的协同周期内的第一代时，其中 t_coor＝1，则子群全局最优位置等于子群S_v-1的上个协同周期结束时得到的 全局最优位置Sg_v-1(t)，从而用于保证子群间的信息交换，否则，对比当前 子群内的最优个体位置与子群历史全局最优位置，将较优者即为子群S_v的当 前全局最优位置：

${\mathrm{Sb}}_{\mathrm{vi}}(t+1)=\left(\begin{array}{cc}{\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)& \mathrm{if}[f_v\left(S_{\mathrm{vi}}(t+1)\right)\ge f_v\left({\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)\right)]\\ S_{\mathrm{vi}}(t+1)& \mathrm{otherwise}\end{array}\right)$式二

${\mathrm{Sg}}_v\left(t\right)=\left(\begin{array}{cc}{\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)& \mathrm{if}[f_v\left({\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)\right)={\min }_{1_S{i}_{S}N}{\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)]\mathrm{and}[t_{\mathrm{coor}}\ne 1]\\ {\mathrm{Sg}}_{v-1}\left(t\right)& \mathrm{if}[t_{\mathrm{coor}}=1]\end{array}\right)$式三。

在上述的一种随机分布型地理要素的多目标协同抽样方案设计方法，所 述的步骤4.4.5中，微粒位置更新即根据子群S_v中微粒S_vi(t)当前个体历史 最佳位置Sb_vi(t)和当前全局历史最优位置Sg_v(t)确定其t+1时刻维度j的取值 x_vi，j(t+1)：

$\left(\begin{array}{c}{V}_{\mathrm{vi},j}(t+1)=\chi ({\mathrm{wV}}_{\mathrm{vi},j}\left(t\right)+c_1{r}_1({\mathrm{Sb}}_{\mathrm{vi},j}\left(t\right)-S_{\mathrm{vi},j}\left(t\right))+c_2{r}_2({\mathrm{Sg}}_{v,j}\left(t\right)-S_{\mathrm{vi},j}\left(t\right)))\\ x_{\mathrm{vi},j}(t+1)=\left(\begin{array}{cc}0& \mathrm{if}[r\ge \mathrm{sig}\left(V_{\mathrm{vi},j}(t+1)\right)]\\ 1& \mathrm{otherwise}\end{array}\right)\end{array}\right)$式四；

微粒位置更新受限于最小协同样本量、采样可达性和采样费用3个因素， 其中最小协同样本量由步骤3计算得到，限制抽样方案中样点的个数；采 样可达性限制包括抽样区域内的建筑物、水域范围或者陡坡区域，限制样 点的空间位置；采样费用限制样点的个数，包括总采样费用为FT、基础采 样费用为FB和单个样点采样费用三部分，单个样点采样费用按抽样区域分 为I级困难区域、II级困难区域和III级困难区域，设I级困难区域样点采样 费用为FS₁，样点数量为N₁，II级每个样点采样费用为FS₂，样点数量为N₂， III级区域为FS₃，样点数量为N₃，采样费用限制条件可以表示为：

FS₁×N₁+FS₂×N₂+FS₃×N₃≤FT-FB    式五。

因此，本发明具有如下优点：能够提高地理要素多指标协同抽样方案 设计效率，节约实地抽样成本，抽样样点布设过程具有自组织性、优化性 和智能性。

附图说明

图1是本发明工作流程示意图。

图2是本发明中微粒映射关系示意图。

图3是本发明中采样不可达区域示意图。

图4是本发明中采用费用与样本量的线性关系图。

图5是本发明实施例中试验区域位置图。

图6是本发明实施例中试验区数据图。

图7是本发明实施例中多土壤指标协同抽样方案设计结果图。

图8是本发明实施例中协同微粒群算法的收敛曲线图。

图9是本发明实施例中最小协同样本量对样本抽样精度的影响图.

图10是本发明实施例中最小协同样本量对算法收敛时间的影响图。

图11是本发明实施例中协同周期对微粒群收敛时间的影响图。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的 说明。

实施例：

本发明所涉及的理论基础如下：

多指标协同抽样方案设计需要满足两个条件：(1)抽样的样本量和样点 布局能够使任一指标的克里金插值结果达到最小精度要求；(2)抽样方法 必须是可量化的(Vasat et al.，2010)。由此可见，自然地理要素多指标 的空间协同抽样属于多目标优化问题的范畴。在抽样过程中，每个指标的 精度限制即为一个优化目标，由于指标的空间分布结构各异，这些优化目 标在同一套样点布局方案中是无法同时得到满足的。优化目标之间的矛盾 表现为当样本方案使一个指标的样本拟合精度达到最优时，另一个指标的 精度可能下降，另外在样本布设过程中还需要考虑最小样本量、采样费用 以及采样可达性等限制。多指标协同空间抽样设计过程就是在采样约束条 件下寻找一组样点布局方案使得所有优化目标通过适当妥协达到帕累托 (Pareto)最优，也即当该抽样方案使一个指标的拟合精度上升时，至少 不会导致其他指标的精度下降。

协同微粒群算法(Coorperative Particle Swarm Optimization，CPSO) 源于对群体智能理论的研究。群体智能理论本质是对自然界中群居型动物 行为进行模拟，探索群体在与外界有物质或能量交换的情况下，如何通过 个体的相互影响和相互合作，自发地完成时间、空间和功能上有序的群体 行为。该理论的研究促进了群体智能计算方法与协同理论的结合。基于协 同思想的群智能计算方法的核心思想为：若优化过程中存在N个优化目标， 则将整个群体划分为N个子群体，每个子群体对应于一个优化目标并以该 目标为个体适宜度判断准则展开并行的独立的最优解搜索，每次迭代过程 后子群体都会将其所获得的最优个体作为社会信息进行共享从而指导其他 子群体完成位置更新。通过这种独立与协同相结合的进化方式，算法的全 局收敛能力能够得到较大提高。

基于种群划分策略的协同微粒群算法对自然地理要素的多指标协同空 间抽样优化方案进行搜索。其基本思想为根据多目标优化问题的优化目标 个数将微粒群划分为V个子种群，对各个子种群进行独立的单目标优化，优 化过程中子种群之间可以交换关于各自最优解的信息，也可执行子种群重 组操作，优化完成后将各个子微粒群获取的Pareto最优解组合形成Pareto 最优解集并作为优化结果输出。

将N维向量微粒群P用于优化具有V个目标的多目标优化问题F(p)时，种 群P首先被划分V个子群体，设S_v为第v个子群体，S_vi为子群S_v中第i个微粒， V_vi(t)为微粒S_vi在t时刻的飞行速度，S_vi(t)为其t时刻的位置，Sb_vi(t)为该微粒 在t时刻经历的历史最优位置，Sg_v(t)为该微粒所属的子微粒群在t时刻的全 局最优位置，Sg(t)为整个微粒群的全局最优位置。

下面是本发明的方法流程步骤：

本发明的实施过程可以划分为宏观过程、中观过程和微观过程三个部 分：宏观过程主要是指子种群划分、子群间协同操作，其最小操作单位为 子种群；中观过程主要包括微粒评价、子种群更新等操作，其最小操作单 位为微粒；微观过程主要包括微粒初始化、微粒位置更新等操作，其最小 单位为微粒维度。具体实施步骤如附图1所示：

步骤1，选择具有随机分布特征的多个地理要素作为抽样指标，获取各 个指标的预采样样点数据，将抽样区域进行格网划分，形成的网格单元作 为抽样单元，抽样区域内的抽样单元集合作为抽样框；

步骤2，利用ESRI公司的ARCGIS 10软件中的普通克里金工具(Ordinary  Kriging)探索预采样数据的空间变异结构并拟合生成各个指标的理论半变 异函数，确定各个指标的空间变异结构拟合函数以及函数参数基台值、块 金值和变程，根据理论半变异函数和指标的预采样样点数据，采用ARCGIS 10 软件中的高斯序贯随机模拟方法(Gaussian sequence randomly simulating  method)得到各个指标的空间分布总体；

步骤3，设置协同微粒群算法参数，包括种群规模p、最大迭代次数、 惯性权重w、随机参数r₁和r₂、个体信息加速常数c₁、社会信息加速常数c₂、 收敛因子x和子种群协同周期T_coor，采用最大迭代次数作为算法终止条件； 根据抽样指标的精度要求和先验方差计算最小协同样本量，具体步骤如下：

对于抽样样本空间A，抽样单元数为N，抽样指标数量为V，设指标i存 在一个最优的抽样设计，其精度要求是p_i，对应的样本容量为n_i。

步骤(3.1)，采用单指标最小样本量计算方法和随机抽样方法确定指 标1的最小样本容量为n₁、最小抽样样本为ξ₁，具体公式如下：

$n_1=\frac{{p_1}^2{s_1}^2}{{(\overline{x_{1n}}-{\mu }_1)}^2}$式一；

式中指标1的精度要求p₁、总体均值μ₁、指标1抽样均值以及标准差 s₁。

步骤(3.2)，依据指标2的抽样精度要求检验样本ξ₁，剔除不符合指 标2抽样设计要求的样点后，形成一组新的样本

步骤(3.3)，依据指标2的精度要求，推算出样本所能代表的最大 子样本空间A₂及其剩余样本空间

步骤(3.4)，在剩余样本空间中依据指标2的抽样设计要求采用 随机抽样方法确定容量为n₂的一个样本ξ₂，则容量为n₁+n₂的样本ξ₁∪ξ₂是同 时满足指标1和指标2抽样设计的一组最小样本。

步骤(3.5)，如此反复，N₀＝n₁+n₂+…+n_v，并且ξ₁，ξ₂，Λ，ξ_v相互独立， 则N₀就是满足多指标抽样设计的最小样本容量。

步骤4，利用协同微粒群算法对现有预采样样点方案进行优化，得到多 指标协同抽样方案，其步骤如下：

步骤(4.1)，构建微粒与抽样区域的映射关系(附图2)，一个二进制 微粒为抽样单元的集合，表示一种抽样方案，微粒中的每一个维度代表一 个抽样单元，维度值为1表示相应抽样单元被选为样点，维度值为0表示 相应抽样单元未被选为样点；依据映射关系和预采样样点数据随机初始化 微粒群S，即将微粒中与预采样样点对应的维度值设置为1；

步骤(4.2)，根据抽样指标数量V等分微粒群，当P被V整除时，其中 子微粒群S_v的种群规模P_v＝P/V；当P不被V整除时，为了避免微粒群微粒信 息丢失，首先确定每个子群包含的基本微粒个数然后将剩余的微粒 依次分配给前个子群，对于未额外分配微粒的子群则从其中随机 选择一个微粒进行复制，保证每个子微粒群S_v的种群规模均为 $P_v=\mathrm{int}\left(\frac{p}{v}\right)+1;$

步骤(4.3)，构建子种群的环状协同结构，以一个随机的子种群为环 状结构的首端，采用随机选取的方式确定其邻域子种群，以此类推直到所 有子种群被组合进环状协同结构。

步骤(4.4)，在一个协同周期T_coor内，子微粒群S_v分别独立寻找一个抽 样指标的最佳抽样方案。设定Flag标识子种群进化是否完成， t_coor(1≤t_coor≤T_coor)用于标识协同周期内子群进化的代数，初始值为1，独 立进化步骤如下：

(4.4.1)以样点的平均克里金方差f₁和信息熵f₂作为S_v中每个微粒S_vi(t) 的适应度衡量准则，计算二准则的加权求和值作为每个微粒的适应度值， 其中权重由用户在实际实施过程中指定：

$f_1=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}(\mathrm{cov}\left(0\right)-W^T\mathrm{cov}\left(x_i{x}_{\mathrm{ij}}\right))$

式中，cov(0)为抽样指标预采样数据的方差，W^T为样点x_i邻近区域内样 点x_ij的权重矩阵，cov(x_ix_ij)为样点x_i邻近区域内样点x_ij的协方差。根据普通 克里金法，协同差为样点x_i到邻近区域内样点x_ij距离的函数，因此f₁的取值 仅与样点的空间分布有关，与样点的具体取值无关。

f₂＝E[-log{f(x₁，x₂，...，x_n)}]

式中，f(x₁，x₂，...，x_n)为其概率密度函数，x_n为步骤2得到的指标空间分 布总体中样点n处的取值。

(4.4.2)确定微粒S_vi在t时刻个体历史最佳位置Sb_vi(t)与微粒群全局最 佳位置Sg_v(t)，具体步骤如下：

当t+1时刻微粒的适应度优于其历史最优的适应度时，则将微粒当前的 位置记录为个体历史最佳位置，否则保持不变；当子群S_v在t+1时刻处于一 个新的协同周期内的第一代时(t_coor＝1)，则子群全局最优位置等于子群 S_v-1的上个协同周期结束时得到的全局最优位置Sg_v-1(t)，从而用于保证子群 间的信息交换，否则，对比当前子群内的最优个体位置与子群历史全局最 优位置，将较优者即为子群S_v的当前全局最优位置：

${\mathrm{Sb}}_{\mathrm{vi}}(t+1)=\left(\begin{array}{cc}{\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)& \mathrm{if}[f_v\left(S_{\mathrm{vi}}(t+1)\right)\ge f_v\left({\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)\right)]\\ S_{\mathrm{vi}}(t+1)& \mathrm{otherwise}\end{array}\right)$式二；

${\mathrm{Sg}}_v\left(t\right)=\left(\begin{array}{cc}{\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)& \mathrm{if}[f_v\left({\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)\right)={\min }_{1_S{i}_{S}N}{\mathrm{Sb}}_{\mathrm{vi}}\left(t\right)]\mathrm{and}[t_{\mathrm{coor}}\ne 1]\\ {\mathrm{Sg}}_{v-1}\left(t\right)& \mathrm{if}[t_{\mathrm{coor}}=1]\end{array}\right)$式三；

(4.4.3)判断是否满足终止条件，如果满足则设定Flag＝true并跳至 (4.5)，并将(4.4.1)中生成的Sg_v(t)作为子群的最优抽样方案输出，否 则将进入(4.4.4)。

(4.4.4)判断是否完成一个协同周期(t_coor＝T_coor)，如果完成则设定 Flag＝false跳至(4.5)，并将(4.4.1)中生成的Sg_v(t)作为子群的最优抽 样方案输出，否则进入(4.4.5)。

(4.4.5)根据微粒当前个体历史最佳位置Sb_vi(t)和当前全局历史最优 位置Sg_vj(t)更新微粒的位置，t_coor自动加1，并返回步骤(4.4.1)，微粒位 置更新即确定t+1时刻子群S_v中微粒S_vi(t)的维度j的取值x_vi，j(t+1)：

$\left(\begin{array}{c}{V}_{\mathrm{vi},j}(t+1)=\chi ({\mathrm{wV}}_{\mathrm{vi},j}\left(t\right)+c_1{r}_1({\mathrm{Sb}}_{\mathrm{vi},j}\left(t\right)-S_{\mathrm{vi},j}\left(t\right))+c_2{r}_2({\mathrm{Sg}}_{v,j}\left(t\right)-S_{\mathrm{vi},j}\left(t\right)))\\ x_{\mathrm{vi},j}(t+1)=\left(\begin{array}{cc}0& \mathrm{if}[r\ge \mathrm{sig}\left(V_{\mathrm{vi},j}(t+1)\right)]\\ 1& \mathrm{otherwise}\end{array}\right)\end{array}\right)$式四；

微粒位置更新受限于最小协同样本量、采样可达性和采样费用3个因 素，其中最小协同样本量由步骤3计算得到，限制抽样方案中样点的个数； 采样可达性限制包括抽样区域内的建筑物、水域范围或者陡坡区域，限制 样点的空间位置(附图3)；采样费用限制样点的个数，包括总采样费用为FT、 基础采样费用为FB和单个样点采样费用三部分(附图4)，根据《测绘生产 成本费用定额》的规定单个样点采样费用按抽样区域可分为I级困难区域 (平原地区为主)、II级困难区域(丘陵地区为主)和III级困难区域(山 地地区为主)三种标准，设I级困难区域样点采样费用为FS₁，样点数量为N₁， II级每个样点采样费用为FS₂，样点数量为N₂，III级区域为FS₃，样点数量 为N₃，采样费用限制条件可以表示为：

FS₁×N₁+FS₂×N₂+FS₃×N₃≤FT-FB    式五；

步骤(4.5)，获取所有子群的当前全局历史最优位置作为最优解组成 Pareto解集，如果Flag＝true进入(4.6)，否则将Pareto解集中的最优解i 分配给相邻子种群S_i+1作为S_i+1下一个协同周期的初始全局最优位置并返回 (4.4)，当i＝V时，则将最优解i分配给子种群S₁。

步骤(4.6)，从Pareto解集中随机选择一个最优解作为自然地理要素 的多指标的协同抽样方案输出。

以下是采用上述方法流程对陕西省横山县进行的实验：

选择陕西省横山县作为试验区(附图5)，依据本发明编写计算机程序， 对试验区的土壤有机质含量(SOM)、土壤含水量(WATER)和速效钾含量(AK) 三个指标进行协同空间抽样方案设计。

(1)实验数据与预处理(附图6)。实验数据包括土地利用现状数据、 1∶5万数字高程模型(Digital Elevation Model，DEM)数据和2004年的 土地质量监测预采样数据。其中，土地利用现状数据和DEM数据作为样点 布设的约束条件，预采样数据作为实验区的先验知识用于对总体均值、方 差等参数以及各个变量空间变异结构的估计。

预采样数据共包括251个样点，有机质含量(SOM)、土壤含水量(WATER)、 速效钾(AK)的统计特征如表1所示：

表1土壤指标的统计特征

对实验数据进行了如下处理：以行政区为抽样底图，采用750m*750m 格网划分底图生成抽样单元，划分后抽样框总体规模为7864；

(2)模拟各个指标的实际总体。对预采样数据中不服从正态分布或近 似正态分布的变量进行正态变换(Log/Normal score变换)；采用普通克里 金法分别对上述三个土壤指标的半变异函数进行拟合，得到拟合结果如表 所示：

表2 SOM、WATER和AK的拟合半变异函数参数

  变量  拟合函数   块金值   基台值   变异程度   SOM  Gaussian   0.029   0.104   0.283   WATER  Spherical   11.235   85.213   0.132   AK  Gaussian   1882.2   1920.5   0.980

以预采样数据为条件约束，在拟合的半变异函数基础上采用序贯高斯 模拟方法分别对上述变量的总体模拟100次，取其平均结果作为各自指标 的实际总体，用于指标数据的重复采样和不同方法获取的抽样结果的对比 分析。

(3)计算最小协同样本量。①、根据土壤有机质含量(SOM)、土壤含 水量(WATER)和速效钾含量(AK)的总体方差和总体规模，在不考虑样点 之间空间关联性和指标之间协同关系的情况下应用下式计算得到了置信水 平为95％、期望精度为90％时各个指标的最小样本量并采用随机抽样方式分 别确定其样点布局(表3)；②、应用4.4.4.1节的流程依次确定协同抽样 时各个指标的最小样本量分别为90、198和79，则最终确定的多指标协同 下的最小样本量为367：

$n_1=\frac{{p_1}^2{s_1}^2}{{(\overline{x_{1n}}-{\mu }_1)}^2}$式六

式中指标1的精度要求p₁、总体均值μ₁、指标1抽样均值以及标准差 s₁。

表3多指标协同下的最小样本量

(4)设置参数，根据预处理结果设置最小协同样本量为367；根据坡度 图和区域实际情况设定采样可达性中的坡度条件为小于等于60°。由于采 样费用与实地采样难易度和采样实施单位的财政能力密切相关，因此较难 确定确切的单个样本费用，实验中未考虑采样费用的约束。协同微粒群算 法参数如表4所示：

表4二进制协同微粒群算法参数设定

其中将惯性权因子w看作是微粒群迭代次数的线性递减函数，并设置其 范围为0.98到0.48，则其计算公式如下：

$w\left(t\right)=0.98-\frac{t}{T_{\mathrm{Max}}}\times 0.5$式七；

其中T_Max为最大迭代次数。

WA方法与模拟退火算法参数：根据变量对土壤质量分类能力高低，加 权平均法(WA)中SOM、WATER和AK的权重分别设为0.408、0.316和0.276 (具体计算过程参见4.5.3.1)；另外，WA中空间模拟退火算法的初始温度T₀设为1，等温系数L为300次，降温系数ρ＝0.25，能量函数采用最小克里金 方差准则。

(5)多土壤指标抽样方案设计。针对SOM、WATER和AK三个指标，采 用协同微粒群方法和加权求和与模拟退火结合法进行多指标协同抽样方案 设计，得到附图7所示的优化结果。对上述优化结果进行统计分析，得到 不同样本的优化目标值和样本量等统计结果。从表5中可以看出：(1)协 同微粒群算法的协同能力整体上要优于WA方法，其对应三个变量的样本MKV 值平均增量小于WA，尤其是它提升了空间变异程度较大的AK变量样本的精 度，并且尽量保持了样本量较低的涨幅。但是对于其中个别指标(SOM)， WA方法在抽样方案设计过程中会偏向于权重最大的指标，从而优先满足其 抽样精度要求，因此WA方法中SOM样本抽样精度的降幅较协同微粒群算法 中该指标抽样精度小。

表5多指标协同样本方案的抽样精度对比

另外从算法效率上，协同微粒群算法也要优于加权求和方法中的模拟 退火算法(SSA)(表6，附图8)。一方面，协同微粒群方法收敛效果相对 较好，其中WATER和AK两个抽样方案的克里金方差目标值均要小于SSA的 收敛结果；另一方面，协同微粒群算法收敛速度较快，显示其收敛所消耗 的时间仅是模拟退火算法的76％，且收敛曲线相对稳定，三个指标的克里金 方差目标值呈稳定下降趋势。

表6协同微粒群算法的收敛结果

(6)抽样方案中最小协同样本量c₁(S)和协同微粒群算法的协同周期 T_coor的影响分析。

①最小协同样本量的影响。以计算的最小协同样本量为中心点，将样 本量分别放大和缩小20％分析其对抽样样本相对精度(标准化到[0，100]范 围内的平均克里金精度)和微粒群效率的影响(附图9和图10)。可见，随 着最小样本量的增加，样本抽样精度逐渐提高，这与传统抽样理论的结论 是一致的，并且上升趋势的显著性以最优协同样本量为界线，当样本量小 于该值时，抽样精度上升较快，而大于该值时则抽样精度提升不再显著； 另外随着样本量增加，微粒群收敛时间不断增加，这主要是由于微粒群优 化算法的搜索空间规模与样本量、样点布局所有可能的组合数量呈正比， 该数量随着最小样本量的增加呈几何趋势增长。但是，由于最小协同样本 量仅是算法优化过程中的一个弹性约束条件，因此其对抽样精度和算法效 率的影响并非完全拟合上述趋势。

②协同周期影响。将协同微粒群算法的协同周期从1调整到10，分析 其对微粒群算法收敛时间的影响，结果如附图11所示。结果显示，协同周 期越小，子微粒群更加趋向于独立进化，子群之间信息交流次数增多，算 法收敛时间变长；随着协同周期增大，算法收敛先呈现线性递减，而当达 到一定周期后(T_coor＝5)，算法收敛时间基本保持不变。但是可以预见，过 大的协同周期将导致子微粒群之间共享信息更新速度边慢，子微粒群之间 协同进化能力退化，从而使得微粒群收敛能力下降。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明 所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或 补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权 利要求书所定义的范围。