CN103955953A - 一种面向数字土壤制图的地形协同变量选取方法 - Google Patents

Abstract

本发明涉及一种面向数字土壤制图的地形协同变量选取方法，采用多地形因子多算法的选择方法，运用了功能测试策略对不同的地形因子变量进行预处理与选择，通过结合其与土壤属性的相关性机制实现了繁杂地形因子变量的快速准确选取，并采用了“评价分析为主，相关分析为辅”技术，实现了“不同地形因子变量，通用选取机制；不同依赖关系，动态因子筛选；评价控制策略，算法性能兼顾”的定量化数字土壤制图地形因子变量选取体系，具有广阔的工业化应用前景。

Description

一种面向数字土壤制图的地形协同变量选取方法

技术领域

本发明涉及一种面向数字土壤制图的地形协同变量选取方法。

背景技术

数字土壤制图（Digital Soil Mapping，简称DSM）是基于环境协变量，利用数学模型进行土壤观测、土壤知识推测土壤类型、属性时空演变分析的土壤信息系统，以优化土壤调查与制图技术、提高土壤信息服务质量为最终目标。传统土壤调查与制图已为各行业、各学科提供了大量的信息支持。目前，各行业迅速膨胀的应用需求、日益丰富的技术手段已对传统土壤制图技术方法体系提出了严峻的挑战。在近30年里，尤其是精准农业、环境管理、土地管理、生态水文模拟等应用对土壤图的精度及时效性提出更高的要求。土壤学理论认为土壤自身的性质与其地形序列有着某种较强的联系，因而土壤信息可以间接地通过地形因子得以表达。通过对有限点的采样，建立环境因素（如典型地形因子）与土壤类型和属性之间的关系模型，如模糊推理、分类树等，可以快速、准确、实时地预测该研究区的土壤类型和属性。数字土壤制图可以预测土壤性质、土壤种类和其他土壤实体，预测模型可以有机地与物理意义明确、易获取的景观属性（协变量）联系起来。按照性质，协变量的种类主要有：地貌、地质、植被群落、侵蚀模式、水文模式、气候状况以及土地利用方式。不同区域土壤形成本身就是一个漫长复杂的过程，期间土壤形成因素也在不断变化。因此，土壤属性的空间变异特征很难用标准的处理流程进行预测分析。

目前，面向数字土壤制图的常用协同变量获取途径主要有两大类：基于DEM的地形分析和遥感信息的自动解译。遥感影像提供了景观的大量有用信息，对遥感信息的自动解译也是土壤景观模型研究的重要方面。快速获取土壤各类理化参数的技术手段是土壤制图的重要研究方向，利用土壤本身的光谱特性，可以从空间连续的遥感观测来提取表层裸土土壤水分、属性分布信息。然而，利用高光谱监测土壤含水量的有效穿透深度有限，仅包含土壤表层几微米的湿度信息。时域反射仪（TDR）与频率反射仪（FDR）是目前最先进的土壤含水量测量仪器，其测量值与烘干法测量值尚有一定差异，应用时需要进行标定。大多数土壤制图方法都限制在小区域，认为决定土壤变化的主要景观因子是地形和水文状况。水文状况的实时数据难以获取与量测，DEM则可以为实际应用提供性价比较高的地形信息，从而成为相关部门使用频率最高的数据来源。

DEM包含的只有高程信息，需要借助数字地形分析(Digital Terrain Analysis，简称DTA）来提取地形协变量。数字地形分析是在数字高程模型(Digital Elevation Model，简称DEM）的基础上进行地形属性计算和特征提取的数字信息处理技术。数字高程模型DEM是用一组有序数值阵列形式表示地面高程的一种实体地面模型，是数字地形模型(DigitalTerrain Model，简称DTM)的一个分支，其它各种地形特征值均可由此派生。DEM的数字表示方法有格网、三角网、等高线三种，其中格网方式因其简单直观、存储方便的优点而被广泛使用。通常应用的协变量主要包括坡度、坡向、曲率、地形湿度指数等因子。

由于数字地形分析的理论研究迅速发展，DEM数据分辨率日益精细，面向不同地形特征的地形因子数量持续增长，从而导致在数字土壤制图过程中选择最优地形因子及最优算法十分困难，归纳起来有以下几点局限性：

（1）土壤协变量研究热点主要集中在土壤与环境关系知识提取方法、新的土壤协变量方面提取、土壤空间推理，其中对于土壤与环境关系知识提取方法还处在定性分析阶段，而对于众多地形协变量的建模分析以及量化研究尚未有所涉及。现有数字地形分析的处理模式多从单一视角出发，未从整体上对数据、物理意义、土壤属性空间分布进行综合整合，也未见有全面探讨三者之间相互关系及量化模型的研究成果，这就使得普通使用者很难把握如何根据地形因子在实际应用过程中建立量化模型并达到效率最大化的目的，在一定程度上制约了数字土壤制图技术的发展和应用。

（2）数字土壤制图的实施仍然依赖已有的土壤调查和知识。反之，历史数据也需要借助全新的地形协变量来提高制图效率和降低制图成本。如何借助已有的基础数据资源、丰富的地形因子提取方法、庞杂的历史土壤数据建立最优的土壤时空演变定量模型仍是一个新的挑战。

（3）地形协变量的选择对制图的结果影响非常大，而且单一的降维操作不能准确选择不同计算方法与不同地形因子的最优集合，预测不同的土壤属性需要选择不同的地形变量。地形分析结果的精度依赖于DEM的精细程度，较为精细的DEM能够提取较为精确的地形因子，然而，超高分辨率（0.5m）的DEM未必最优。因此，土壤制图过程中如何选择特定尺度最优地形因子组合也已彰显出其特有的重要性。

发明内容

本发明所要解决的技术问题是提供一种基于功能测试与量化处理进行地形协同变量选取的面向数字土壤制图的地形协同变量选取方法。

本发明为了解决上述技术问题采用以下技术方案：本发明设计了一种面向数字土壤制图的地形协同变量选取方法，包括如下步骤：

步骤01.针对各个土壤样本点对应的各个地形因子变量V_i,j进行预处理，使其符合正态分布；其中，1≤i≤m，1≤j≤n_i，m为地形因子的种类数，n_i为对应第i种地形因子的计算方法的数量；

步骤02.针对经步骤01处理后的各个地形因子变量V_i,j与各个土壤样本点的土壤属性数据S_long,lati分别进行单地形因子多算法的皮尔逊相关性分析，分别获得各种地形因子对应其各种计算方法与该土壤样本点的土壤属性数据ΣS_long,lati之间的相关性结果pr_i,j，并分别构成各种地形因子对应其各种计算方法与土壤属性数据ΣS_long,lati之间的相关性集合 ${\mathrm{pr}}_i=({\mathrm{pr}}_{i,1},\·\·\·,{\mathrm{pr}}_{i,j},\·\·\·,{\mathrm{pr}}_{{i,n}_i});$

步骤03.分别针对各种地形因子对应其各种计算方法与土壤属性数据ΣS_long,lati之间的相关性集合中的各个相关性结果pr_i,j，分别采用如下公式（1）进行评价：

${\mathrm{LDTASET}}_{i,j}=|\mathrm{\λ}\×\frac{{\mathrm{pr}}_{i,j}}{T_{i,j}}|\×\mathrm{\ξ}---\left(1\right)$

针对各个pr_i集合，将各个集合中的各LDTASET_i,j按大小降序排列，由大至小取出前u个LDTASET_i,j，分别获得各种地形因子分别对应其u种最优计算方法的地形因子变量V′_i,j，并分别构成各种地形因子对应其最优计算方法的地形因子变量集合V_i；其中，LDTASET_i,j为各种地形因子对应其各种计算方法的评价结果；T_i,j表示对应第i种地形因子的第j种计算方法的计算时间；ξ表示地形因子栅格值的分布是否符合正态分布，符合ξ为1，否则ξ为0；λ为膨胀系数，为预设计算精度控制系数；u为预设选取地形因子对应其最优算法的个数，且u≤min(n_i)；

步骤04.针对中各种地形因子对应其最优计算方法的地形因子变量集合V_i中的各个地形因子变量V′_i,j进行线性变化，并进行降维操作；

步骤05.针对步骤04中进行降维操作后获得的各种地形因子对应其最优计算方法的地形因子变量集合V_i中的各个地形因子变量V′_i,j与各个土壤样本点的土壤属性数据S_long,lati分别进行多地形因子单算法的皮尔逊相关性分析，分别获得各种地形因子对应其最优计算方法与土壤属性数据ΣS_long,lati之间的相关性结果pr′_i,j；

步骤06.针对步骤05中获得的各种地形因子对应其最优计算方法与土壤属性数据ΣS_long,lati之间的相关性结果pr′_i,j，分别采用如下公式（2）进行评价：

${\mathrm{MDTASET}}_{i,j}=\left|\frac{{\mathrm{pr}}_{i,j}^{\′}}{T_{i,j}}\right|---\left(2\right)$

其中，MDTASET_i,j为该土壤样本中各种地形因子对应其最优计算方法的评价结果；

步骤07.根据降维操作后获得的地形因子的种类数c，以及其各自对应的最优的计算方法的数量u，获得(c×u)种地形因子变量组合，针对该(c×u)种地形因子变量组合，根据如下公式（3）：

${\mathrm{MLDTA}}_k=\frac{\underset{k=1}{\overset{c\×m}{\mathrm{\Σ}}}|{\mathrm{pr}}_k^{\′}-{\mathrm{pr}}_{k-1}^{\′}|}{{\mathrm{MDTASET}}_k-{\mathrm{MDTASET}}_{k-1}}---\left(3\right)$

获得各种地形因子变量组合的评价结果MLDTA_k，1≤k≤(c×u)；取出MLDTA_k的最大值，根据对应地形因子变量组合获取地形协同变量。

作为本发明的一种优选技术方案：所述步骤02中按照如下公式：

${\mathrm{pr}}_{i,j}=\frac{\mathrm{b\Σ}{S}_{i,j}{S}_{\mathrm{long},\mathrm{lati}}-\mathrm{\Σ}{V}_{i,j}\mathrm{\Σ}{S}_{\mathrm{long},\mathrm{lati}}}{\sqrt{\mathrm{b\Σ}{V}_{i,j}^2-{\left(\mathrm{\Σ}{V}_{i,j}\right)}^2}\sqrt{\mathrm{b\Σ}{S}_{\mathrm{long},\mathrm{lati}}^2-{\left(\mathrm{\Σ}{S}_{\mathrm{long},\mathrm{lati}}\right)}^2}}$

针对所述经步骤01处理后的各个地形因子变量V_i,j与各个土壤样本点的土壤属性数据S_long,lati分别进行单地形因子多算法的皮尔逊相关性分析，分别获得各种地形因子对应其各种计算方法与土壤属性数据ΣS_long,lati之间的相关性结果pr_i,j，b为所述土壤样本点的数量。

作为本发明的一种优选技术方案：所述步骤04中.按如下公式（4）进行降维操作，

$\left\{\begin{array}{c}S=a_{11}{V}_1+a_{21}{V}_2+\·\·\·a_{m1}{V}_m\\ S=a_{12}{V}_1+a_{22}{V}_2+\·\·\·a_{m2}{V}_m\\ \·\\ \·\\ \·\\ S=a_{1m}{V}_1+a_{2m}{V}_2+\·\·\·a_{\mathrm{mm}}{V}_m\end{array}\right.---\left(4\right)$

其中，S为土壤属性数据, $a=\left[\begin{array}{cccc}{a}_{11}& a_{21}& \·\·\·& a_{m1}\\ a_{12}& a_{22}& \·\·\·& a_{m2}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ a_{1m}& a_{2m}& \·\·\·& a_{\mathrm{mm}}\end{array}\right]$ 为预设降维系数矩阵，(a_1h)²+(a_2h)²+…+(a_mh)²＝1，h∈{1、…、m}。

作为本发明的一种优选技术方案：所述膨胀系数λ通过如下步骤进行预设：

步骤a.获取所述各种地形因子的平均值mean，初始nmean＝0；

步骤b.迭代计算是否成立，若成立，则nmean的数值加1并赋值给参数nmean；mean的数值除以10并赋值给参数mean，继续执行本步骤，直到判断条件不成立；

步骤c.判断是否成立，若成立则计算结束，λ＝10^nmean。

本发明所述一种面向数字土壤制图的地形协同变量选取方法采用以上技术方案与现有技术相比，具有以下技术效果：

（1）本发明设计的面向数字土壤制图的地形协同变量选取方法，采用多地形因子多算法的选择方法，运用了功能测试策略对不同的地形因子变量进行预处理与选择，通过结合其与土壤属性的相关性机制实现了繁杂地形因子变量的快速准确选取，并采用了“评价分析为主，相关分析为辅”技术，实现了“不同地形因子变量，通用选取机制；不同依赖关系，动态因子筛选；评价控制策略，算法性能兼顾”的定量化数字土壤制图地形因子变量选取体系，具有广阔的工业化应用前景；

（2）本发明设计的面向数字土壤制图的地形协同变量选取方法中，提出的单地形因子多算法功能测试方法，结合算法复杂度的差异性及地形因子变量与土壤属性空间相关的特点，以皮尔逊相关性分析结果作为基数，这样使得地形因子变量在选取时充分考虑了其适用性，从而提高土壤制图的精度；

（3）本发明设计的面向数字土壤制图的地形协同变量选取方法中，提出多地形因子单算法功能测试方法，使用主成分分析的方法对原始数据变量进行降维操作，在保证地形因子变量数量合理的情况下，通过计算相关性分析与算法复杂度的加权评价指标，确定最优测试用例集合，避免了频繁执行土壤制图方法，大幅度提升了计算效率。

附图说明

图1是本发明设计的面向数字土壤制图的地形协同变量选取方法的流程示意图。

具体实施方式

下面结合说明书附图对本发明的具体实施方式作进一步详细的说明。

针对多种地形分析方法、多种计算方法，对土壤制图过程中的地形因子变量进行抽象，具体模型的数学表示为：

F＝(S_long,lati,V_i,j,R)

其中，F是土壤制图结果，S_long,lati是各个土壤样本点的土壤属性数据，long、lati分别是土壤样本点经纬度位置数据，V_i,j为地形因子变量，1≤i≤m，1≤j≤n_i，m为地形因子的种类数，n_i为对应第i种地形因子的计算方法的数量，R是其他的协同变量，如土地利用、年均气温等。V_i,j地形因子变量主要包括数据分辨率（Resolution），数据维度（Col列×Row行），起始坐标（RCol×RRow），每一个栅格的字节数（Bit），原始DEM生产方法（Type），水平分辨率误差（Ehorizon），垂直分辨率误差（Evertical），投影坐标（Projection）等。

地形因子变量之间的相关性包括：不相关、弱相关、相关三种。根据计算方法流程，地形因子间的关系分为：派生关系、无派生关系。

构建地形因子变量分布矩阵V，dem是二维高程栅格数据。

$V=\left(\begin{array}{ccccc}{V}_{\mathrm{1,1}}\left(\mathrm{dem}\right)& V_{\mathrm{1,2}}\left(\mathrm{dem}\right)& V_{\mathrm{1,3}}\left(\mathrm{dem}\right)& \·\·\·& V_{1,n1}\left(\mathrm{dem}\right)\\ V_{\mathrm{2,1}}\left(\mathrm{dem}\right)& V_{\mathrm{2,2}}\left(\mathrm{dem}\right)& V_{\mathrm{2,3}}\left(\mathrm{dem}\right)& \·\·\·& V_{2,n2}\left(\mathrm{dem}\right)\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ V_{m,1}\left(\mathrm{dem}\right)& V_{m,2}\left(\mathrm{dem}\right)& V_{m,3}\left(\mathrm{dem}\right)& \·\·\·& V_{m,\mathrm{nm}}\left(\mathrm{dem}\right)\end{array}\right)$

如图1所示，本发明设计的面向数字土壤制图的地形协同变量选取方法具体包括如下步骤：

步骤01.针对各个土壤样本点对应的各个地形因子变量V_i,j进行预处理，使其符合正态分布；其中，1≤i≤m，1≤j≤n_i，m为地形因子的种类数，n_i为对应第i种地形因子的计算方法的数量，其中，各个地形因子变量的正态分布可以通过Kolmogorov-Smirnov检验法进行检验。

步骤02.针对经步骤01处理后的各个地形因子变量V_i,j与各个土壤样本点的土壤属性数据S_long,lati按照如下公式：

${\mathrm{pr}}_{i,j}=\frac{\mathrm{b\Σ}{S}_{i,j}{S}_{\mathrm{long},\mathrm{lati}}-\mathrm{\Σ}{V}_{i,j}\mathrm{\Σ}{S}_{\mathrm{long},\mathrm{lati}}}{\sqrt{\mathrm{b\Σ}{V}_{i,j}^2-{\left(\mathrm{\Σ}{V}_{i,j}\right)}^2}\sqrt{\mathrm{b\Σ}{S}_{\mathrm{long},\mathrm{lati}}^2-{\left(\mathrm{\Σ}{S}_{\mathrm{long},\mathrm{lati}}\right)}^2}}$

分别进行单地形因子多算法的皮尔逊相关性分析，分别获得各种地形因子对应其各种计算方法与该土壤样本点的土壤属性数据ΣS_long,lati之间的相关性结果pr_i,j，并分别构成各种地形因子对应其各种计算方法与土壤属性数据ΣS_long,lati之间的相关性集合其中，b为所述土壤样本点的数量；其中，各个土壤样本点的土壤属性数据S_long,lati的读取，随机选取土壤属性数据的75%作为测试用例数据，剩余的25%作为验证数据集合。

其中基本的相关性分析评价如下表所示：

排序等级	相关系数大小	解释
			1	0～0.2	非常弱的相关或无关
2	0.2～0.4	弱相关
			3	0.4～0.6	中度相关
4	0.6～0.8	强相关
			5	0.8～1.0	非常强的相关

步骤03.分别针对各种地形因子对应其各种计算方法与土壤属性数据ΣS_long,lati之间的相关性集合中的各个相关性结果pr_i,j，分别采用如下公式（1）进行评价：

${\mathrm{LDTASET}}_{i,j}=|\mathrm{\λ}\×\frac{{\mathrm{pr}}_{i,j}}{T_{i,j}}|\×\mathrm{\ξ}---\left(1\right)$

针对各个pr_i集合，将各个集合中的各LDTASET_i,j按大小降序排列，由大至小取出前u个LDTASET_i,j，分别获得各种地形因子分别对应其u种最优计算方法的地形因子变量V′_i,j，并分别构成各种地形因子对应其最优计算方法的地形因子变量集合V_i；其中，LDTASET_i,j为各种地形因子对应其各种计算方法的评价结果；T_i,j表示对应第i种地形因子的第j种计算方法的计算时间；ξ表示地形因子栅格值的分布是否符合正态分布，符合ξ为1，否则ξ为0；u为预设选取地形因子对应其最优算法的个数，且u≤min(n_i)；λ为膨胀系数，λ≥1，为预设计算精度控制系数，取值为最接近当前地形因子平均值的10的n次方，需要根据地形因子的平均值来确定，所述膨胀系数λ通过如下步骤进行预设：

步骤a.获取所述各种地形因子的平均值mean，初始nmean＝0；

步骤b.迭代计算是否成立，若成立，则nmean的数值加1并赋值给参数nmean；mean的数值除以10并赋值给参数mean，继续执行本步骤，直到判断条件不成立；

步骤c.判断是否成立，若成立则计算结束，λ＝10^nmean。

步骤04.针对中各种地形因子对应其最优计算方法的地形因子变量集合V_i中的各个地形因子变量V′_i,j进行线性变化，并按如下公式（4）进行降维操作；

$\left\{\begin{array}{c}S=a_{11}{V}_1+a_{21}{V}_2+\·\·\·a_{m1}{V}_m\\ S=a_{12}{V}_1+a_{22}{V}_2+\·\·\·a_{m2}{V}_m\\ \·\\ \·\\ \·\\ S=a_{1m}{V}_1+a_{2m}{V}_2+\·\·\·a_{\mathrm{mm}}{V}_m\end{array}\right.---\left(4\right)$

其中，S为土壤属性数据, $a=\left[\begin{array}{cccc}{a}_{11}& a_{21}& \·\·\·& a_{m1}\\ a_{12}& a_{22}& \·\·\·& a_{m2}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ a_{1m}& a_{2m}& \·\·\·& a_{\mathrm{mm}}\end{array}\right]$ 为预设降维系数矩阵，(a_1h)²+(a_2h)²+…+(a_mh)²＝1，h∈{1、…、m}；

步骤05.针对步骤04中进行降维操作后获得的各种地形因子对应其最优计算方法的地形因子变量集合V_i中的各个地形因子变量V′_i,j与各个土壤样本点的土壤属性数据S_long,lati按照如下公式

${\mathrm{pr}}_{i,j}=\frac{\mathrm{b\Σ}{S}_{i,j}{S}_{\mathrm{long},\mathrm{lati}}-\mathrm{\Σ}{V}_{i,j}\mathrm{\Σ}{S}_{\mathrm{long},\mathrm{lati}}}{\sqrt{\mathrm{b\Σ}{V}_{i,j}^2-{\left(\mathrm{\Σ}{V}_{i,j}\right)}^2}\sqrt{\mathrm{b\Σ}{S}_{\mathrm{long},\mathrm{lati}}^2-{\left(\mathrm{\Σ}{S}_{\mathrm{long},\mathrm{lati}}\right)}^2}}$

分别进行多地形因子单算法的皮尔逊相关性分析，分别获得各种地形因子对应其最优计算方法与土壤属性数据ΣS_long,lati之间的相关性结果pr′_i,j；

步骤06.针对步骤05中获得的各种地形因子对应其最优计算方法与土壤属性数据ΣS_long,lati之间的相关性结果pr′_i,j，分别采用如下公式（2）进行评价：

${\mathrm{MDTASET}}_{i,j}=\left|\frac{{\mathrm{pr}}_{i,j}^{\′}}{T_{i,j}}\right|---\left(2\right)$

其中，MDTASET_i,j为该土壤样本中各种地形因子对应其最优计算方法的评价结果；

步骤07.根据降维操作后获得的地形因子的种类数c，以及其各自对应的最优的计算方法的数量u，获得(c×u)种地形因子变量组合，针对该(c×u)种地形因子变量组合，根据如下公式（3）：

${\mathrm{MLDTA}}_k=\frac{\underset{k=1}{\overset{c\×m}{\mathrm{\Σ}}}|{\mathrm{pr}}_k^{\′}-{\mathrm{pr}}_{k-1}^{\′}|}{{\mathrm{MDTASET}}_k-{\mathrm{MDTASET}}_{k-1}}---\left(3\right)$

获得各种地形因子变量组合的评价结果MLDTA_k，1≤k≤(c×u)；取出MLDTA_k的最大值，根据对应地形因子变量组合获取地形协同变量。

在获得所需土壤样本点的地形协同变量后，基于专家知识，对于具有特定含义的地形因子与经验最佳的算法进行选取，然后进行土壤制图与精度验证，其中，精度验证方式可以采用交叉验证的方式进行验证。

本发明设计的面向数字土壤制图的地形协同变量选取方法，采用多地形因子多算法的选择方法，运用了功能测试策略对不同的地形因子变量进行预处理与选择，通过结合其与土壤属性的相关性机制实现了繁杂地形因子变量的快速准确选取，并采用了“评价分析为主，相关分析为辅”技术，实现了“不同地形因子变量，通用选取机制；不同依赖关系，动态因子筛选；评价控制策略，算法性能兼顾”的定量化数字土壤制图地形因子变量选取体系，具有广阔的工业化应用前景；本发明设计的面向数字土壤制图的地形协同变量选取方法中，还提出了单地形因子多算法功能测试方法，结合算法复杂度的差异性及地形因子变量与土壤属性空间相关的特点，以皮尔逊相关性分析结果作为基数，这样使得地形因子变量在选取时充分考虑了其适用性，从而提高土壤制图的精度；以及在本发明方法中，提出了多地形因子单算法功能测试方法，使用主成分分析的方法对原始数据变量进行降维操作，在保证地形因子变量数量合理的情况下，通过计算相关性分析与算法复杂度的加权评价指标，确定最优测试用例集合，避免了频繁执行土壤制图方法，大幅度提升了计算效率。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (4)

1.一种面向数字土壤制图的地形协同变量选取方法，其特征在于，包括如下步骤：

步骤01.针对各个土壤样本点对应的各个地形因子变量V_i,j进行预处理，使其符合正态分布；其中，1≤i≤m，1≤j≤n_i，m为地形因子的种类数，n_i为对应第i种地形因子的计算方法的数量；

步骤02.针对经步骤01处理后的各个地形因子变量V_i,j与各个土壤样本点的土壤属性数据S_long,lati分别进行单地形因子多算法的皮尔逊相关性分析，分别获得各种地形因子对应其各种计算方法与该土壤样本点的土壤属性数据ΣS_long,lati之间的相关性结果pr_i,j，并分别构成各种地形因子对应其各种计算方法与土壤属性数据ΣS_long,lati之间的相关性集合 ${\mathrm{pr}}_i=({\mathrm{pr}}_{i,1},\·\·\·,{\mathrm{pr}}_{i,j},\·\·\·,{\mathrm{pr}}_{{i,n}_i});$

步骤03.分别针对各种地形因子对应其各种计算方法与土壤属性数据ΣS_long,lati之间的相关性集合中的各个相关性结果pr_i,j，分别采用如下公式（1）进行评价：

${\mathrm{LDTASET}}_{i,j}=|\mathrm{\λ}\×\frac{{\mathrm{pr}}_{i,j}}{T_{i,j}}|\×\mathrm{\ξ}---\left(1\right)$

针对各个pr_i集合，将各个集合中的各LDTASET_i,j按大小降序排列，由大至小取出前u个LDTASET_i,j，分别获得各种地形因子分别对应其u种最优计算方法的地形因子变量V′_i,j，并分别构成各种地形因子对应其最优计算方法的地形因子变量集合V_i；其中，LDTASET_i,j为各种地形因子对应其各种计算方法的评价结果；T_i,j表示对应第i种地形因子的第j种计算方法的计算时间；ξ表示地形因子栅格值的分布是否符合正态分布，符合ξ为1，否则ξ为0；λ为膨胀系数，为预设计算精度控制系数；u为预设选取地形因子对应其最优算法的个数，且u≤min(n_i)；

步骤04.针对中各种地形因子对应其最优计算方法的地形因子变量集合V_i中的各个地形因子变量V′_i,j进行线性变化，并进行降维操作；

步骤05.针对步骤04中进行降维操作后获得的各种地形因子对应其最优计算方法的地形因子变量集合V_i中的各个地形因子变量V′_i,j与各个土壤样本点的土壤属性数据S_long,lati分别进行多地形因子单算法的皮尔逊相关性分析，分别获得各种地形因子对应其最优计算方法与土壤属性数据ΣS_long,lati之间的相关性结果pr′_i,j；

步骤06.针对步骤05中获得的各种地形因子对应其最优计算方法与土壤属性数据ΣS_long,lati之间的相关性结果pr′_i,j，分别采用如下公式（2）进行评价：

${\mathrm{MDTASET}}_{i,j}=\left|\frac{{\mathrm{pr}}_{i,j}^{\′}}{T_{i,j}}\right|---\left(2\right)$

其中，MDTASET_i,j为该土壤样本中各种地形因子对应其最优计算方法的评价结果；

步骤07.根据降维操作后获得的地形因子的种类数c，以及其各自对应的最优的计算方法的数量u，获得(c×u)种地形因子变量组合，针对该(c×u)种地形因子变量组合，根据如下公式（3）：

${\mathrm{MLDTA}}_k=\frac{\underset{k=1}{\overset{c\×m}{\mathrm{\Σ}}}|{\mathrm{pr}}_k^{\′}-{\mathrm{pr}}_{k-1}^{\′}|}{{\mathrm{MDTASET}}_k-{\mathrm{MDTASET}}_{k-1}}---\left(3\right)$

获得各种地形因子变量组合的评价结果MLDTA_k，1≤k≤(c×u)；取出MLDTA_k的最大值，根据对应地形因子变量组合获取地形协同变量。

2.根据权利要求1所述一种面向数字土壤制图的地形协同变量选取方法，其特征在于：所述步骤02中按照如下公式：

${\mathrm{pr}}_{i,j}=\frac{\mathrm{b\Σ}{S}_{i,j}{S}_{\mathrm{long},\mathrm{lati}}-\mathrm{\Σ}{V}_{i,j}\mathrm{\Σ}{S}_{\mathrm{long},\mathrm{lati}}}{\sqrt{\mathrm{b\Σ}{V}_{i,j}^2-{\left(\mathrm{\Σ}{V}_{i,j}\right)}^2}\sqrt{\mathrm{b\Σ}{S}_{\mathrm{long},\mathrm{lati}}^2-{\left(\mathrm{\Σ}{S}_{\mathrm{long},\mathrm{lati}}\right)}^2}}$

针对所述经步骤01处理后的各个地形因子变量V_i,j与各个土壤样本点的土壤属性数据S_long,lati分别进行单地形因子多算法的皮尔逊相关性分析，分别获得各种地形因子对应其各种计算方法与土壤属性数据ΣS_long,lati之间的相关性结果pr_i,j，b为所述土壤样本点的数量。

3.根据权利要求1所述一种面向数字土壤制图的地形协同变量选取方法，其特征在于：所述步骤04中.按如下公式（4）进行降维操作，

$\left\{\begin{array}{c}S=a_{11}{V}_1+a_{21}{V}_2+\·\·\·a_{m1}{V}_m\\ S=a_{12}{V}_1+a_{22}{V}_2+\·\·\·a_{m2}{V}_m\\ \·\\ \·\\ \·\\ S=a_{1m}{V}_1+a_{2m}{V}_2+\·\·\·a_{\mathrm{mm}}{V}_m\end{array}\right.---\left(4\right)$

其中，S为土壤属性数据, $a=\left[\begin{array}{cccc}{a}_{11}& a_{21}& \·\·\·& a_{m1}\\ a_{12}& a_{22}& \·\·\·& a_{m2}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ a_{1m}& a_{2m}& \·\·\·& a_{\mathrm{mm}}\end{array}\right]$ 为预设降维系数矩阵，(a_1h)²+(a_2h)²+…+(a_mh)²＝1，h∈{1、…、m}。

4.根据权利要求1所述一种面向数字土壤制图的地形协同变量选取方法，其特征在于，所述膨胀系数λ通过如下步骤进行预设：

步骤a.获取所述各种地形因子的平均值mean，初始nmean＝0；

步骤b.迭代计算是否成立，若成立，则nmean的数值加1并赋值给参数nmean；mean的数值除以10并赋值给参数mean，继续执行本步骤，直到判断条件不成立；

步骤c.判断是否成立，若成立则计算结束，λ＝10^nmean。