一种顾及时空分布非平稳特征的土壤有机碳空间抽样网络设计方法

摘要

本发明涉及一种顾及时空分布非平稳特征的土壤有机碳空间抽样网络设计方法，首先提取土壤有机碳空间抽样网络设计的基础数据并进行整合，然后应用时空回归方法建立土壤有机碳时空变化与多时间序列属性之间的回归关系，将土壤有机碳时空变化过程映射为一个由土壤有机碳多因素时空回归趋势和回归残差构成的时空随机过程，最后构建土壤有机碳空间抽样网络设计模型，以上述数据为模型输入数据，评价适应度函数，建立抽样方案设计问题与粒子群算法之间的映射关系，求解得出土壤有机碳空间抽样网络优化方案。该方法大幅度提高调查精度，同时充分发挥粒子群算法快速收敛和超强寻优能力，提升土壤有机碳空间抽样网络的合理性、适用性和设计效率。

说明书

技术领域

本发明涉及土地资源调查领域，尤其是涉及一种土壤有机碳空间抽样网络设计方法。

背景技术

土壤有机碳库是陆地生态系统中处于活跃状态的最大碳库，在减缓碳释放、调节全球碳平衡方面发挥着关键作用。准确测定土壤有机碳储量已经成为当前科学界和各国政府关注的焦点问题。在我国，随着“践行节能低碳、建设美丽中国”目标的提出，及时开展高精度的土壤有机碳抽样调查与监测工作，将为科学制定低碳发展战略提供有力的决策支持。面向科研与发展的多层面需求，研究高效的土壤有机碳时空抽样方法势在必行。

时空分布非平稳性是影响土壤有机碳抽样设计的关键因素。在生态系统类型、土壤生物活性、群落组成及土地利用等因素共同作用下，土壤有机碳时空变量的均值、方差具有高度的时空位置变化相关性，导致有机碳时空分布自相关结构(依赖性)呈现出阶跃、线性连续或非线性连续等不平稳特征，从而增加了时空抽样网络对其的拟合难度。准确刻画土壤有机碳时空非平稳特征有助于显著提升时空抽样方法的实用性。

现有时空抽样方法，包括完全重复抽样、完全替换抽样和适应性抽样方法，在描述时空非平稳性方面仍存在一定不足。具体表现在：(1)完全重复抽样与完全替换抽样仅能描述变量的空间非平稳性，忽略了其时间变化特征；(2)适应性抽样方法涉及了采用分层方法或几何比例法描述时空非平稳性及时空各向异性，但仍无法准确刻画时空连续非平稳特征。如何避免传统时空抽样网络设计方法的片面性，建立有效的土壤有机碳连续非平稳时空特征描述模型并用于指导抽样网络的优化布局，对于精准、高效开展土壤有机碳抽样与监测工作尤为重要。

发明内容

本发明主要是解决现有技术所存在的技术问题；提供了一种通过准确描述土壤有机碳时空分布特征，有效提高土壤有机碳空间抽样调查精度的一种土壤有机碳抽样网络优化设计方法。

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

一种顾及时空分布非平稳特征的土壤有机碳空间抽样网络设计方法，其特征在于，该方法包括如下步骤：

步骤1，提取土壤有机碳空间抽样网络设计的基础数据并进行整合，整合后的数据为地理空间栅格单元的多时间序列属性信息，其中，基础数据包括多时期的土壤有机碳空间分布数据、土地利用与覆被变化数据、数字高程模型数据；整合数据的具体方法是对空间数据进行配准并设置一致的空间参考坐标，将矢量数据转成精度相同的地理空间栅格数据，以及通过分析数字高程模型数据生成栅格单元的高程、坡度坡向、坡度位置、表面起伏度一系列地形指数信息。多时间序列属性信息包括地理空间栅格单元的多时期土壤有机碳含量、多时期土地利用与覆被类型以及高程、坡度坡向、坡度位置和表面起伏度信息。

步骤2，基于应用时空回归方法分析步骤1的数据，以地理空间栅格单元为数据处理单元，建立土壤有机碳变化与多时间序列属性的回归关系模型，将土壤有机碳时空变化过程映射为一个由土壤有机碳多因素时空回归趋势和回归残差构成的时空随机过程，包括以下子步骤：

步骤2.1，将土壤有机碳时空变化过程映射为一个由土壤有机碳多因素时空回归趋势m(x,t)和回归残差ε(x,t)构成的时空随机过程(式1)。

$>Z(x,t)=\underset{i=1}{\overset{N}{\Sigma }}{a}_i{b}_i(x,t)+ϵ(x,t),t=1,2,\mathrm{...},T,x=1,2,\mathrm{...},N---\left(1\right)$>

式中，x表示空间位置，t表示时间节点，Z(x,t)表示t时刻空间位置x的土壤有机碳含量实际值。b_i(x,t)为t时刻空间位置x处的属性b_i的值，a_i为回归系数。

步骤2.2，将每个地理空间栅格对应的土壤有机碳值、多时间序列属性值代入式1，采用Matlab软件的最小二乘法求取所有回归系数b_i。

步骤2.3，在所有b_i和多时间序列属性值已知情况下，计算每个栅格土壤有机碳的估计值，将估计值与真实值对比，其差值即为回归残差ε(x,t)。

步骤3，采用时空不可分协方差方法描述步骤2的土壤有机碳时空回归残差，建立土壤有机碳时空变化协方差模型。

步骤4，以随机生成的抽样网络为模型输入，应用二进制粒子群算法进行建模并对随机抽样网络进行优化设计。以地理空间栅格为数据处理单元，建立抽样网络布局问题与二进制粒子群算法之间的映射关系，求解得出土壤有机碳空间抽样网络的优化布局。所述地理空间栅格单元为二进制粒子群算法中的粒子维度，其值表征该单元是否作为抽样点；土壤有机碳空间抽样方案对应于粒子；抽样网络效能评价函数对应于粒子适应度函数，由土壤有机碳时空变化协方差模型推导出的最小预测误差表征。

在上述的一种顾及时空分布非平稳特征的土壤有机碳空间抽样网络设计方法，所述的步骤3中，采用时空不可分协方差方法建立土壤有机碳回归残差时空变化协方差模型的步骤包括以下子步骤：

步骤3.1，将步骤2的土壤有机碳回归残差的协方差看作空间维度协方差C_s与时间维度协方差C_t的组合：

C_st＝κ₁C_s+κ₂C_t-κ₃C_sC_t>

式中，κ₁,κ₂,κ₃为关联系数的组合(参数组合)。

步骤3.2，根据协方差与方差的关系，分别将土壤有机碳空间维度协方差C_s与时间维度协方差C_t描述为以下模型：

C_s(x_i,x_j|θ_s)＝σ_s(x_i|θ_s)σ_s(x_j|θ_s)ρ_s(|x_i-x_j|),i,j＝1,2,...,N>

C_t(x_i(t),x_j(t)|θ_t)＝σ_t(x_i(t)|θ_t)σ_t(x_j(t)|θ_t)ρ_t(|x_i(t)-x_j(t)|),t＝1,2,...,T>

式中，C_s,C_t分别为空间、时间维度上的协方差函数；σ_s,σ_t分别代表空间、时间维度上的方差分布函数，θ_s,θ_t为方差函数的未知参数，ρ_s(x_i,x_j),ρ_t(xi(t),x_j(t))分别表示空间维度和时间维度上两点之间的相关性，只与点间空间维度距离和时间维度距离h_s,h_t相关。

在上述的一种顾及时空分布非平稳特征的土壤有机碳空间抽样网络设计方法，所述的步骤4中，求解土壤有机碳空间抽样网络优化布局的步骤包括以下子步骤：

步骤4.1，设置抽样网络样点数量，在研究区范围内随机生成包含所设置的样点数量的抽样网络。

步骤4.2，粒子群种群规模，惯性权重，个体认知参数、社会认知参数、最大迭代次数，并以步骤4.1生成的随机抽样网络初始化每个粒子的维度值；

步骤4.3，设置粒子群算法的适应度函数；

步骤4.4，计算各个粒子的适应度，通过对比粒子适应度选择粒子自身所经历最优位置以及当前种群最优位置；

步骤4.5，根据粒子自身所经历最优位置和当前种群最优位置对粒子进行位置变更操作，生成新一批粒子；

步骤4.6，循环迭代，当达到最大迭代次数时搜索结束，将当前种群最优位置作为最优土壤有机碳空间抽样设计方案输出，否则返回步骤4.3。

在上述的一种顾及时空分布非平稳特征的土壤有机碳空间抽样网络设计方法，所述步骤4.4中，粒子位置更新操作在于对粒子的维度值进行变更，维度值取1代表对应地理栅格单元被选为抽样单元，取0代表为非抽样单元。粒子位置变更遵循采样可达性原则，即当粒子维度对应的地理栅格单元落在水域、建设用地或坡度大于60°区域内时，该栅格单元被抽样概率为0，同时粒子相应维度值保持0不变。

本发明具有如下优点：1、充分考虑了土壤有机碳时间维度和空间维度变化特征，在此基础上设计的空间抽样网络能够准确拟合土壤有机碳变化特征，提高调查精度；2、粒子群算法具备快速收敛和超强寻优能力，将土壤有机碳空间网络设计映射为一种二元组合优化问题，可以在适应度函数指导下优化空间抽样网络，以提高土壤有机碳空间抽样网络的合理性和适用性。

附图说明

图1是本发明的模型流程图。

图2是本发明的粒子抽样网络映射关系示意图。

图3是本发明的粒子位置更新操作示意图。

图4a是本发明的实施例不同样点数量下抽样网络优化结果示意图(样点数100)。

图4b是本发明的实施例不同样点数量下抽样网络优化结果示意图(样点数200)。

图4c是本发明的实施例不同样点数量下抽样网络优化结果示意图(样点数300)。

具体实施方式

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

本发明所采用的模型流程图如图1示。

该顾及时空分布非平稳特征的土壤有机碳空间抽样网络设计方法包括如下步骤：

步骤1，提取土壤有机碳空间抽样网络设计的基础数据并进行整合，整合后的数据为地理空间栅格单元的多时间序列属性信息。提取基础数据是指从其他数据库或其他系统获取抽样网络设计所需的数据，如多时期的土壤有机碳空间分布数据、土地利用与覆被变化数据、数字高程模型数据等。整合数据是指这些数据是多个时期、多种来源、不同格式的数据，要对这些时空数据进行统一规范化处理。使用ArcGIS 10.2.1的Spatial Adjustment工具对上述空间数据进行配准并设置一致的空间参考坐标，使用Feature to Raster工具将矢量数据转成精度相同的地理空间栅格数据，使用Surface analyst工具分析数字高程模型数据，生成栅格单元的高程、坡度坡向、坡度位置、表面起伏度一系列地形指数信息。多时间序列属性信息包括地理空间栅格单元的多时期土壤有机碳含量、多时期土地利用与覆被类型以及高程、坡度坡向、坡度位置和表面起伏度信息。

步骤2，以地理空间栅格单元为数据处理单元，应用时空回归方法分析抽样网络设计的基础数据，建立土壤有机碳变化与多时间序列属性的回归关系模型，将土壤有机碳时空变化过程映射为一个由土壤有机碳多因素时空回归趋势m(x,t)和回归残差ε(x,t)构成的时空随机过程(式1)。

$>Z(x,t)=\underset{i=1}{\overset{N}{\Sigma }}{a}_i{b}_i(x,t)+ϵ(x,t),t=1,2,\mathrm{...},T,x=1,2,\mathrm{...},N---\left(1\right)$>

式中，x表示空间位置，t表示时间节点，Z(x,t)表示t时刻空间位置x的土壤有机碳含量实际值。b_i(x,t)为t时刻空间位置x处的属性b_i的值，a_i为回归系数。

步骤3，计算土壤有机碳时空回归残差。将每个地理空间栅格对应的土壤有机碳值、多时间序列属性值代入式1，采用Matlab软件的最小二乘法求取所有回归系数b_i。在所有b_i和多时间序列属性值已知情况下，计算每个栅格土壤有机碳的估计值，将估计值与真实值对比，其差值即为回归残差ε(x,t)。

步骤4，采用时空不可分协方差方法描述土壤有机碳时空回归残差ε(x,t)，建立土壤有机碳时空变化协方差模型。具体将土壤有机碳回归残差的协方差看作空间维度协方差C_s与时间维度协方差C_t的组合：

C_st＝κ₁C_s+κ₂C_t-κ₃C_sC_t>

式中，κ₁,κ₂,κ₃为关联系数的组合，根据空间维度与时间维度土壤有机碳变异大小进行设置。

步骤5，对步骤4构建的土壤有机碳时空变化协方差模型进行实例化。首先根据协方差与方差的关系，分别将土壤有机碳空间维度协方差C_s与时间维度协方差C_t描述为以下模型：

C_s(x_i,x_j|θ_s)＝σ_s(x_i|θ_s)σ_s(x_j|θ_s)ρ_s(|x_i-x_j|),i,j＝1,2,…,N>

C_t(x_i(t),x_j(t)|θ_t)＝σ_t(x_i(t)|θ_t)σ_t(x_j(t)|θ_t)ρ_t(|x_i(t)-x_j(t)|),t＝1,2,...,T>

式中，C_s,C_t分别为空间、时间维度上的协方差函数；σ_s,σ_t分别代表空间、时间维度上的方差分布函数，θ_s,θ_t为方差函数的未知参数；ρ_s(x_i,x_j),ρ_t(x_i(t),x_j(t))分别表示空间维度和时间维度上两点之间的相关性，只与点间空间维度距离和时间维度距离h_s,h_t相关。

本发明涉及的土壤有机碳时空维度上的方差分布函数σ_s,σ_t表示一阶(均值)平稳条件下可能存在的二阶过程(方差)变化，可以采用K级阶跃函数表示如下：

$>{\sigma }_s\left(x_i\right|{\theta }_s)={\theta }_{sl}I(x_i\in D_l)+...+{\theta }_{sK}I(x_i\in D_K),\forall x_i\in D_k,D_k⋐D---\left(5\right)$>

$>{\sigma }_t\left(x_t\right(t\left)\right|{\theta }_t)={\theta }_{tl}I\left(x_i\right(l)\in D_i(l\left)\right)+...+{\theta }_{tK}I\left(x_i\right(K)\in D_i(K\left)\right),\forall x_i\left(t\right)\in D_i\left(K\right),D_i\left(K\right)⋐D---\left(6\right)$>

式中，θ_s,θ_t分别为空间维度和时间维度K级阶跃函数的参数组合，分别由K个θ_sk,θ_tk组成；θ_sk,θ_tk分别表示空间维度和时间维度上x落在第k个二阶平稳子域(D_k或D_i(k))内的方差。

本发明涉及的时间维度和空间维度上二阶平稳子域采用ARCGIS10.2软件的自然断裂法(Natural break)工具进行划分，划分后求取属于同一子域的地理栅格单元的土壤有机碳实际值的方差(式7)即为相应子域对应的K级阶跃函数的参数θ_sk,θ_tk。

$>{\theta }_{sk}=\frac{1}{N_k}\underset{i=1}{\overset{N_k}{\Sigma }}{\left(z\left(x_i\right)-\overline{z\left(x_i\right)}\right)}^2,x_i\in D_k,$>

$>{\theta }_{tk}=\frac{1}{N_k}\underset{t=1}{\overset{N_k}{\Sigma }}{\left(z\left(x_i\left(t\right)\right)-\overline{z\left(x_i\left(t\right)\right)}\right)}^2,x_i\left(t\right)\in D_i\left(k\right)---\left(7\right)$>

式中，Z(x_i)表示空间位置x_i的土壤有机碳含量所有时刻实际值的均值，Z(x_i(t))表示t时刻空间位置x_i的土壤有机碳含量实际值，表示所有属于子域D_k范围内的栅格单元对应的土壤有机碳含量在所有时刻实际值的均值，表示空间位置x_i所有属于子域D_i(k)时间序列内的栅格单元对应的土壤有机碳含量真实值的均值。

本发明涉及的ρ_s(x_i,x_j),ρ_t(x_i(t),x_j(t))分别表示空间维度和时间维度上两点之间的相关性，可以采用理论模型进行拟合，如拟合为指数模型ρ_s(x_si,x_sj)＝exp(-h_s)，ρ_t(x_ti,x_tj)＝exp(-h_t)。

步骤6，设置空间抽样网络样点数量，采用ARCGIS10.2软件generate randompoints工具生成随机抽样网络，应用二进制粒子群算法进行建模并对随机抽样网络进行优化设计。以地理空间栅格为数据处理单元，建立抽样网络布局问题与二进制粒子群算法之间的映射关系，求解得出土壤有机碳空间抽样网络的优化布局。所述地理空间栅格单元为二进制粒子群算法中的粒子维度，其值表征该单元是否作为抽样点；土壤有机碳空间抽样方案对应于粒子；抽样网络效能评价函数对应于粒子适应度函数，由土壤有机碳时空变化协方差模型推导出的最小克里金误差表征。

步骤7，设置粒子群种群规模，惯性权重，个体认知参数、社会认知参数、最大迭代次数，并以步骤6生成的随机抽样网络初始化每个粒子的维度值。

步骤8，根据适应度函数，计算每个粒子的适应度。具体方法是在土壤有机碳回归残差的协方差模型基础上，推导最小克里金方差准则作为土壤有机碳空间抽样网络的评价函数，也即粒子群算法的适应度函数，计算公式如下：

$>\begin{array}{cc}Minimize& F_{mkv}=C_{st}\left(0\right)-{\Sigma }_{j-1}^N{\lambda }_j·C_{st}\left(x_i\left(t\right),x_j\left(t\right)\right)+\mu \end{array}---\left(5\right)$>

式中，F_MKV为适应度函数；C_st(0)为土壤有机碳时空分布的先验方差；为克里格差值权重系数，μ为拉格朗日系数；C_st(x_i(t),x_j(t))为待添加样点x_i(t)与已存在样点x_j(t)的协方差，根据公式2、3、4计算得到。

本发明涉及的克里格插值权重系数和拉格朗日系数可以根据以下公式求得：

$>\left(\begin{array}{cccccc}0& r_{12}& r_{13}& \mathrm{...}& r_{1n}& 1\\ r_{21}& 0& \mathrm{...}& \mathrm{...}& r_{2n}& 1\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ r_{n1}& r_{n2}& \mathrm{...}& \mathrm{...}& 0& 1\\ 1& 1& \mathrm{...}& \mathrm{...}& \mathrm{...}& 0\end{array}\right)\left(\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ \mathrm{...}\\ {\lambda }_n\\ \mu \end{array}\right)=\left(\begin{array}{c}{r}_{10}\\ r_{20}\\ r_{n0}\\ 1\end{array}\right)---\left(6\right)$>

式中，r_nn代表已存在样点之间的协方差C_st，r_n0代表待添加样点与已存在样点之间的协方差C_st，协方差值均根据公式2、3、4计算得到。

步骤9，通过对比粒子适应度选择粒子自身所经历最优位置以及当前种群最优位置；根据粒子自身所经历最优位置和当前种群最优位置对粒子进行位置变更操作，生成新一批粒子；

$>\begin{array}{c}{v}_m(t+1)={\mathrm{wv}}_m\left(t\right)+c_1{r}_1\left({\mathrm{Pb}}_m\right(t)-x_m(t\left)\right)+c_2{r}_2\left(P_g\right(t)-x_m(t\left)\right)\\ x_m(t+1)=x_m\left(t\right)+v_m(t+1)\end{array}---\left(7\right)$>

式中，v_m代表粒子m维度的速度，x_m为粒子m维度的取值，Pb_i(t)和P_g(t)分别为应用适应度函数确定的t时刻粒子所经过的最佳位置和当前种群全局最佳位置，惯性权重w主要用于调整粒子群的全局搜索能力和局部搜索能力从而加速算法收敛。个体认知参数c₁和社会认知参数c₂能够起到调整粒子群收敛速度的作用。

本发明涉及的粒子位置更新操作在于应用粒子群速度更新公式，对粒子的维度值进行变更，变更后的x_m为连续值，采用Sigmoid函数将x_m映射到[0,1]之间，当映射值大于0.5时，粒子维度值取1，表示该维度对应的栅格单元被选为抽样单元，否则维度值取0代表为非抽样单元。粒子位置变更同时遵循采样可达性原则，即当粒子维度对应的地理栅格单元落在水域、建设用地或坡度大于60°区域内时，该栅格单元被抽样概率为0，同时粒子相应维度值保持0不变。

步骤10，循环迭代，当达到最大迭代次数时搜索结束，否则返回步骤8。

步骤11，选择出粒子群中适应度最高的粒子，将其转换为土壤有机碳空间抽样网络，得到抽样网络的优化设计结果。

实施例1：

1、提取某一区域的土壤有机碳抽样网络设计的基础数据并进行整合，整合后的数据为各个地理空间栅格单元的10年序列属性信息。

2、以地理栅格单元为基本单元，应用时空回归方法分析每个单元土壤有机碳值与不同时间节点土地利用与覆被类型、高程、坡度坡向、坡度位置和表面起伏度等属性的回归关系，建立时空回归方程。应用回归方程计算地理栅格单元土壤有机碳值估计量，降之与真实值求差得到土壤有条件时空回归的回归残差。

3、建立土壤有机碳时空变化协方差模型，根据研究区土壤有机碳时空维度变异系数设置κ₁,κ₂,κ₃分别为0.5，0.8和0.2。采用ARCGIS10.2软件natural>

4、分别设置抽样网络样点数量为100、200和300，应用ARCGIS10.2软件的Generaterandom points随机生成抽样网络方案。

5、设置粒子群种群规模为20个，惯性权重0.85，个体认知参数1.5，社会认知参数0.8，最大迭代次数600次，利用随机生成的抽样网络对所有粒子进行初始化操作。

6、根据适应度函数，计算每个粒子的适应度。

7、根据适应度函数值，通过对比粒子适应度选择粒子自身所经历最优位置以及当前种群最优位置；根据粒子自身所经历最优位置和当前种群最优位置对粒子进行位置变更操作，生成新一批粒子。

8、循环迭代，当达到最大迭代次数时搜索结束，否则返回步骤6。

9、选择出粒子群中适应度最高的粒子，将其转换为土壤有机碳空间抽样网络，得到不同样点数量情况下抽样网络的优化设计结果。

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