一种面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法

摘要

本发明涉及一种面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法，通过从多时相亮温/陆面温度数据集中，任意选取两不同时序节点的亮温/陆面温度数据，构建多种指数（差异指数、比值指数、归一化差异指数）来反映目标土壤的空间变异特征，两不同时序节点的选择，并不局限于某一个或某一种具体的外界环境条件（如夏季或者降雨前后）的限制，而是首先计算所有可能时相组合下的指数，并利用目标土壤属性与所有时相亮温/陆面温度数据以及所有时相组合下指数之间的相关程度，并结合预设筛选规则，综合选择确定针对目标土壤属性预测与制图的输入变量。

说明书

技术领域

本发明涉及一种面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法，属于土壤属性协同变量构建技术领域。

背景技术

土壤是由固相、液相、气相三相共存的复杂异质体系，其形成受到气候、生物、地形、母质以及时间因素的影响，使得土壤在时空上具有强烈的变异特征。地统计学方法在土壤属性的空间预测与制图方面发挥着重要的作用，该方法通过空间插值技术(如普通Kriging)实现土壤的空间表达。但这种方法对土壤样品采集模式以及采样点的密度具有较高的要求，考虑到土壤强烈的时空变异特性，为了获得较高的土壤属性预测与制图精度，往往需要较多的土壤样品，其在大尺度土壤属性预测与制图方面的应用往往会受到采样密度的限制。随着数字地形以及遥感技术的发展，各种地形变量以及遥感手段获取的反映土壤和植被特征的参数，被广泛应用于土壤属性预测与制图；并且随着机器学习方法的不断发展，土壤属性与这些变量/参数间的关系得以更好的揭示，预测与制图的精度得到不断提高。

当前，用于土壤属性预测与制图的遥感参数主要来自于光学遥感影像，光学遥感影像多是记录地物的反射率特征，例如常用的归一化差异植被指数(NormalizedDifference Vegetation Index，NDVI)便是基于光学影像近红外与红色波段的反射率进行计算。与土壤反射率特征相比，土壤的温度属性对外界环境的变化更为敏感，除了受到土壤自身属性(如土壤质地、有机质等)以及土壤水分含量等因素的影响外，还受到大气温度的影响。因此，土壤温度是其自身特征与外部因素影响的综合反映。土壤温度信息以及相似外部条件下(如相同大气温度)土壤温度变化程度在不同空间位置处的差异，具有反映土壤空间变异特征的潜力。

发明内容

本发明所要解决的技术问题是提供一种面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法，能够有效、精准实现协同变量的构建与选择，提高实际工作效率。

本发明为了解决上述技术问题采用以下技术方案：本发明设计了一种面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法，用于获得目标土壤对应于面向土壤属性预测与制图的协同变量，其特征在于，包括如下步骤：

步骤A.根据目标土壤的土壤属性的时间变异特征，结合下垫面土壤以及植被状况，综合考虑目标土壤所处环境因素，确定热红外影像的选取时间区间，然后进入步骤B；

步骤B.获取目标土壤对应于选取时间区间的时间序列热红外影像，并计算或直接获取其中对应于各时序节点的亮温/陆面温度数据，构成目标土壤对应选取时间区间的亮温/陆面温度数据集T，然后进入步骤C；其中，所有各时序节点均分别对应于亮温数据，或者所有各时序节点均分别对应于陆面温度数据；

步骤C.采用时-空插值方法，针对目标土壤所对应选取时间区间内亮温/陆面温度数据集T中无法使用的区域进行填补，并以此针对目标土壤对应选取时间区间的亮温/陆面温度数据集T进行更新，然后进入步骤D；

步骤D.针对亮温/陆面温度数据集T中的各个数据，按两个不同时序节点所对应数据划分为一组，获得各个数据的所有分组；并针对各个分组，获得各分组分别所对应的差异指数、比值指数、归一化差异指数，构建差异指数集合DI、比值指数集合RI、归一化差异指数集合NDI，然后进入步骤E；

步骤E.分别计算亮温/陆面温度数据集T中各数据相对目标土壤的土壤属性的相关系数，进而构建针对目标土壤中土壤属性所对应的亮温/陆面温度相关系数集RT；

同时，分别针对差异指数集合DI、比值指数集合RI、归一化差异指数集合NDI，计算获得集合中各数据相对目标土壤的土壤属性的相关系数，进而构建针对目标土壤中土壤属性所对应的差异指数相关系数集RDI、比值指数相关系数集RRI、归一化差异指数相关系数集RNDI；然后进入步骤F；

步骤F.根据目标土壤所对应的亮温/陆面温度相关系数集RT、差异指数相关系数集RDI、比值指数相关系数集RRI、归一化差异指数相关系数集RNDI，确定获得目标土壤对应于面向土壤属性预测与制图的协同变量。

作为本发明的一种优选技术方案，所述步骤B包括如下步骤：

步骤B1.获取目标土壤对应于选取时间区间的时间序列热红外影像，然后进入步骤B2；

步骤B2.针对目标土壤所对应时间序列热红外影像，进行辐射定标、几何校正、大气校正预处理更新，然后进入步骤B3；

步骤B3.针对目标土壤所对应时间序列热红外影像中的各个时序节点，计算时序节点所对应的亮温数据，当陆面发射率准确获取时，则计算时序节点所对应的陆面温度数据，并对各亮温/陆面温度数据按时间先后顺序进行排序，构成目标土壤对应选取时间区间的亮温/陆面温度数据集T＝{T₁、T₂、…、T_n}，n为目标土壤对应于选取时间区间的时间序列热红外影像的时序节点数目，然后进入步骤C；其中，所有各时序节点均分别对应于亮温数据，或者所有各时序节点均分别对应于陆面温度数据。

作为本发明的一种优选技术方案，所述步骤B3中，针对目标土壤所对应时间序列热红外影像中的各个时序节点，基于普朗克定律，计算影像亮温数据，并且基于陆面发射率准确获取时，计算陆面温度数据。

作为本发明的一种优选技术方案，所述步骤D包括如下步骤：

步骤D1.针对亮温/陆面温度数据集T中的各个数据，按两个不同时序节点所对应数据划分为一组，获得各个数据的所有分组，分组的个数为M，M＝C(n,2)＝n！/(2！(n-2)！)，然后进入步骤D2；

步骤D2.分别针对各个分组，按如下公式：

DI_m＝T_i-T_j

RI_m＝T_i/T_j

NDI_m＝(T_i-T_j)/(T_i+T_j)

获得各分组分别所对应的差异指数、比值指数、归一化差异指数；其中，T_i∈T，T_j∈T，i＞j，1≤m≤M，DI_m表示第m个分组所对应的差异指数，RI_m表示第m个分组所对应的比值指数，NDI_m表示第m个分组所对应的归一化差异指数，然后进入步骤D3；

步骤D3.根据各分组分别所对应的差异指数、比值指数、归一化差异指数，构建差异指数集合DI＝{DI₁、DI₂、…、DI_M}、比值指数集合RI＝{RI₁、RI₂、…、RI_M}、归一化差异指数集合NDI＝{NDI₁、NDI₂、…、NDI_M}，然后进入步骤E；

作为本发明的一种优选技术方案，所述步骤F包括如下步骤：

步骤F1.按照R＝RT∪RDI∪RRI∪RNDI，针对目标土壤所对应的亮温/陆面温度相关系数集RT、差异指数相关系数集RDI、比值指数相关系数集RRI、归一化差异指数相关系数集RNDI进行合并，获得相关系数数据集合R，然后进入步骤F2；

步骤F2.获取相关系数数据集合R中的最大值Max(R)和最小值Min(R)，然后进入步骤F3；

步骤F3.按照L＝T∪DI∪RI∪NDI，针对亮温/陆面温度数据集T、差异指数集合DI、比值指数集合RI、归一化差异指数集合NDI进行合并，获得数据集合L，然后查找获得Max(R)、Min(R)分别对应数据集合L中的数据，然后进入步骤F4；

步骤F4.按预设面向土壤属性预测与制图的协同变量的筛选规则，针对步骤F3所查找获得的数据进行筛选，获得目标土壤对应于面向土壤属性预测与制图的协同变量。

本发明所述一种面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法采用以上技术方案与现有技术相比，具有以下技术效果：

本发明所设计一种面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法，与现有技术相比，通过从多时相亮温/陆面温度数据集中，任意选取两不同时序节点的亮温/陆面温度数据，构建多种指数(差异指数、比值指数、归一化差异指数)来反映目标土壤的空间变异特征，两不同时序节点的选择，并不局限于某一个或某一种具体的外界环境条件(如夏季或者降雨前后)的限制，而是首先计算所有可能时相组合下的指数，并利用目标土壤属性与所有时相亮温/陆面温度数据以及所有时相组合下指数之间的相关程度，并结合预设筛选规则，综合选择确定针对目标土壤属性预测与制图的输入变量。

附图说明

图1是本发明所设计面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法的流程图；

图2是实施案例中两期陆面温度数据时-空插值前后的对比图；

图3a至图3d是实施案例中土壤粘粒含量与陆面温度以及各指数间的相关系数直方图。

具体实施方式

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

基于热红外影像获取的亮温或陆面温度，能够有效反映土壤及其它地物的温度信息，而且作为一种被动式获取方式，不会受到太阳辐射的限制，夜晚仍能够进行影像的获取。除卫星遥感影像平台外，随着航空、无人机等技术的快速发展，基于航空、无人机等快速部署平台的热红外遥感，能够更加方便地获取多时相(时间序列)高频率的亮温与陆面温度数据。通过获取较高时间频率的热红外影像，能够更加有效地揭示土壤温度随时间变化的规律与程度，不同空间位置处土壤温度的差异以及土壤温度时间变化规律的差异能够反映土壤的空间变异特征。因此，基于不同时相亮温与陆面温度数据构建的反映土壤温度变化的指数在土壤属性预测与制图方面具有重要价值。

本发明设计了一种面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法，用于获得目标土壤对应于面向土壤属性预测与制图的协同变量，如图1所示，实际应用中，具体包括如下步骤：

步骤A.根据目标土壤的土壤属性的时间变异特征，结合下垫面土壤以及植被状况，综合考虑目标土壤所处环境因素，确定热红外影像的选取时间区间，然后进入步骤B。

步骤B.获取目标土壤对应于选取时间区间的时间序列热红外影像，并计算或直接获取其中对应于各时序节点的亮温/陆面温度数据，构成目标土壤对应选取时间区间的亮温/陆面温度数据集T，然后进入步骤C；其中，所有各时序节点均分别对应于亮温数据，或者所有各时序节点均分别对应于陆面温度数据。

上述步骤B在实际应用中，具体执行如下步骤B1至步骤B3。

步骤B1.获取目标土壤对应于选取时间区间的时间序列热红外影像，然后进入步骤B2。

步骤B2.针对目标土壤所对应时间序列热红外影像，进行辐射定标、几何校正、大气校正预处理更新，然后进入步骤B3。

步骤B3.针对目标土壤所对应时间序列热红外影像中的各个时序节点，基于普朗克定律(Planck’s law)，计算时序节点所对应的亮温数据，当陆面发射率准确获取时，则计算时序节点所对应的陆面温度数据，并对各亮温/陆面温度数据按时间先后顺序进行排序，构成目标土壤对应选取时间区间的亮温/陆面温度数据集T＝{T₁、T₂、…、T_n}，n为目标土壤对应于选取时间区间的时间序列热红外影像的时序节点数目，然后进入步骤C；其中，所有各时序节点均分别对应于亮温数据，或者所有各时序节点均分别对应于陆面温度数据。

步骤C.采用时-空插值方法，针对目标土壤所对应选取时间区间内亮温/陆面温度数据集T中无法使用的区域进行填补，并以此针对目标土壤对应选取时间区间的亮温/陆面温度数据集T进行更新，然后进入步骤D。此步骤，即填补了由于大气中水、云等对影像产生干扰而无法有效使用的区域。

步骤D.针对亮温/陆面温度数据集T中的各个数据，按两个不同时序节点所对应数据划分为一组，获得各个数据的所有分组；并针对各个分组，获得各分组分别所对应的差异指数、比值指数、归一化差异指数，构建差异指数集合DI(Difference index)、比值指数集合RI(Ratio index)、归一化差异指数集合NDI(Normalized difference index)，然后进入步骤E。

步骤D在实际应用中，具体执行如下步骤D1至步骤D3。

步骤D1.针对亮温/陆面温度数据集T中的各个数据，按两个不同时序节点所对应数据划分为一组，获得各个数据的所有分组，分组的个数为M，M＝C(n,2)＝n！/(2！(n-2)！)，然后进入步骤D2。

步骤D2.分别针对各个分组，按如下公式：

DI_m＝T_i-T_j

RI_m＝T_i/T_j

NDI_m＝(T_i-T_j)/(T_i+T_j)

获得各分组分别所对应的差异指数、比值指数、归一化差异指数；其中，T_i∈T，T_j∈T，i＞j，1≤m≤M，DI_m表示第m个分组所对应的差异指数，RI_m表示第m个分组所对应的比值指数，NDI_m表示第m个分组所对应的归一化差异指数，然后进入步骤D3。

步骤D3.根据各分组分别所对应的差异指数、比值指数、归一化差异指数，构建差异指数集合DI＝{DI₁、DI₂、…、DI_M}、比值指数集合RI＝{RI₁、RI₂、…、RI_M}、归一化差异指数集合NDI＝{NDI₁、NDI₂、…、NDI_M}，然后进入步骤E。

步骤E.分别计算亮温/陆面温度数据集T中各数据相对目标土壤的土壤属性的相关系数，进而构建针对目标土壤中土壤属性所对应的亮温/陆面温度相关系数集RT；

同时，分别针对差异指数集合DI、比值指数集合RI、归一化差异指数集合NDI，计算获得集合中各数据相对目标土壤的土壤属性的相关系数，进而构建针对目标土壤中土壤属性所对应的差异指数相关系数集RDI、比值指数相关系数集RRI、归一化差异指数相关系数集RNDI；然后进入步骤F。

步骤F.根据目标土壤所对应的亮温/陆面温度相关系数集RT、差异指数相关系数集RDI、比值指数相关系数集RRI、归一化差异指数相关系数集RNDI，确定获得目标土壤对应于面向土壤属性预测与制图的协同变量。

上述步骤F在实际应用中，具体执行如下步骤F1至步骤F3。

步骤F1.按照R＝RT∪RDI∪RRI∪RNDI，针对目标土壤所对应的亮温/陆面温度相关系数集RT、差异指数相关系数集RDI、比值指数相关系数集RRI、归一化差异指数相关系数集RNDI进行合并，获得相关系数数据集合R，然后进入步骤F2。

步骤F2.获取相关系数数据集合R中的最大值Max(R)和最小值Min(R)，然后进入步骤F3。

步骤F3.按照L＝T∪DI∪RI∪NDI，针对亮温/陆面温度数据集T、差异指数集合DI、比值指数集合RI、归一化差异指数集合NDI进行合并，获得数据集合L，然后查找获得Max(R)、Min(R)分别对应数据集合L中的数据，然后进入步骤F4。

步骤F4.按预设面向土壤属性预测与制图的协同变量的筛选规则，针对步骤F3所查找获得的数据进行筛选，获得目标土壤对应于面向土壤属性预测与制图的协同变量。

本发明的基本思想是利用土壤的温度属性来反映土壤的空间变异特征，而影像亮温/陆面温度是土壤温度属性的有效表现形式，因此利用亮温/陆面温度能够反映土壤的空间变异特征；多时相亮温/陆面温度能够反映土壤温度的时间变化规律，选择两期不同时相亮温/陆面温度构建的差值指数、比值指数、归一化差值指数是反映土壤温度变化规律的一种简单表示方式，这种不同空间位置处土壤温度的变化程度同样能够反映土壤的空间变异特征；通过建立多时相较高频率的亮温/陆面温度数据集，并计算任意两期不同时相亮温/陆面温度数据构建的指数，基于目标土壤属性与亮温/陆面温度数据以及两期亮温/陆面温度构建指数之间的相关系数，并结合预设筛选规则，实现面向土壤属性预测与制图过程中的变量选择。

以MODIS(Moderate Resolution Imaging Spectroradiometer)八天合成的陆面温度数据(MOD11A2)为例，该数据集包括白天和夜晚的陆面温度产品，空间分辨率为1km，实施例中仅选用白天陆面温度数据；以河南省封丘县为研究区区域，研究区域农作物种植方式为冬小麦-夏玉米，以土壤粘粒含量为目标土壤属性，土壤样品数量为117，采自2014年6月夏小麦收获之后。如图1所示，本发明所设计面向土壤属性预测与制图的亮温/陆面温度协同变量构建选择方法在实际应用过程当中，具体包括如下步骤：

步骤1.依据土壤粘粒含量短时间变异较小的特征，结合下垫面土壤以及作物生长规律，综合考虑气温、降水等因素，同时考虑到土壤样品采集时间是2014年6月，本实施案例中遥感热红外影像的时间区间设定为[2013年1月1日，2015年12月31]，涵盖三年的时间跨度。选择这一较长时间段内的影像，既能够考虑土壤粘粒含量时间变异较小的特点，又能够避免短时间段内可能存在的异常因素对亮温/陆面温度数据的干扰。

步骤A.对于热红外影像的选取时间区间的选择，在本实施例中，依据土壤粘粒含量短时间变异较小的特征，结合下垫面土壤以及作物生长规律，综合考虑气温、降水等因素，同时考虑到土壤样品采集时间是2014年6月，本实施案例中热红外影像的选取时间区间设定为[2013年1月1日，2015年12月31]，涵盖三年的时间跨度，选择这一较长时间段内的影像，既能够考虑土壤粘粒含量时间变异较小的特点，又能够避免短时间段内可能存在的异常因素对亮温/陆面温度数据的干扰，然后进入步骤B。

接下来针对实施例，通过步骤B1至步骤B3的执行，实现设计方案中步骤B的操作，然后进入步骤C。

步骤B1.获取实施例中目标土壤对应于选取时间区间[2013年1月1日，2015年12月31]的时间序列热红外影像，然后进入步骤B2。

然后依次执行步骤B2、步骤B3。

MODIS MOD11A2数据集能够直接提供陆面温度产品，而且具有较高的精度，因此，本实施案例中直接使用MOD11A2数据集中提供的陆面温度数据，不再进行步骤b和步骤c的计算。本实施案例中仅选用MOD11A2数据集中的白天(卫星过境时间为当地时间10：30AM)陆面温度数据，经过步骤B1至步骤B3，形成实施例中目标土壤对应选取时间区间的亮温/陆面温度数据集T＝{T₁、…、T₁₃₈}，共计138期陆面温度数据。

步骤C.采用时-空插值方法，针对实施例中目标土壤对应选取时间区间的亮温/陆面温度数据集T＝{T₁、…、T₁₃₈}中无法使用的区域进行填补，并以此针对目标土壤对应选取时间区间的亮温/陆面温度数据集T进行更新，然后进入步骤D。此步骤，即填补了由于大气中水、云等对影像产生干扰而无法有效使用的区域。

时-空插值方法既能够考虑地物在水平地理空间的依赖性，又能够考虑其在时序上的变化，因此在填补陆面温度数据中的无效数据区域时具有更好的效果。本实施例中采用的时-空插值方法，公开在下列文件中：Gerber,F.,Jong,R.d.,Schaepman,M.E.,Schaepman-Strub,G.,Furrer,R.,2018.Predicting Missing Values in Spatio-Temporal Remote Sensing Data.IEEE Transactions on Geoscience and RemoteSensing,56(5),2841-2853.doi:10.1109/TGRS.2017.2785240。实施案例中陆面温度数据时-空插值前后的对比结果如图2所示。

接下来针对实施例，通过步骤D1至步骤D3的执行，实现设计方案中步骤D的操作，然后进入步骤E。

步骤D1.针对实施例中目标土壤对应选取时间区间的亮温/陆面温度数据集T＝{T₁、…、T₁₃₈}中的各个数据，按两个不同时序节点所对应数据划分为一组，获得各个数据的所有分组，分组的个数为M，M＝C(n,2)＝n！/(2！(n-2)！)，然后进入步骤D2。本实施例中，M＝C(138,2)＝138！/(2！(138-2)！)＝9453。

然后执行步骤D2的操作后，继续执行如下步骤D3。

步骤D3.根据各分组分别所对应的差异指数、比值指数、归一化差异指数，构建差异指数集合DI＝{DI₁、DI₂、…、DI₉₄₅₃}、比值指数集合RI＝{RI₁、RI₂、…、RI₉₄₅₃}、归一化差异指数集合NDI＝{NDI₁、NDI₂、…、NDI₉₄₅₃}，然后进入步骤E。

按上述所设计步骤E进行执行，本实施案例中，土壤粘粒含量与陆面温度及各指数间相关系数直方图如图3a至图3d所示。之后进入步骤F，并具体执行步骤F1至步骤F4，其中，对于步骤F2，其中对应于本实施例，土壤粘粒含量与陆面温度及各指数间相关系数的最大值为0.61，最小值为-0.57，均能够通过显著性水平p＝0.01的显著性检验；最后执行步骤F3，对应于本实施例，土壤粘粒含量与各指数间相关系数的最大值均大于与陆面温度之间的相关系数，最小值均小于与陆面温度之间的相关系数；而且各指数间具有相等的相关系数最大值与最小值。

最后执行步骤F4，对应于本实施例，以所查找各亮温/陆面温度数据或指数为土壤属性预测与制图模型的输入，结合预设筛选规则，确定面向土壤属性预测与制图模型的输入变量(亮温/陆面温度或指数)。本实施案例中，土壤粘粒含量与各指数间具有相等的相关系数最大值与最小值，而且所对应的陆面温度日期也相同，考虑到NDI指数具有范围固定([-1,1])的优点，在作为土壤属性预测与制图模型的输入变量时具有一定的优势，本实施案例中选用NDI指数；土壤粘粒含量与NDI指数间相关系数的最大值对应的两期陆面温度日期分别为2014年12月19日、2014年10月8日，最小值对应的两期陆面温度日期分别为2014年10月8日、2014年1月17日。在2014年1月17日及12月19日附近，冬小麦并未完全覆盖土壤表面，2014年10月8日前后，处于夏玉米收获及冬小麦播种期，此时土壤表面小麦覆盖程度最小。因此，可以将基于2014年12月19日与2014年10月8日，以及2014年10月8日与2014年1月17日间陆面温度构建的NDI指数作为土壤粘粒含量预测与制图模型的输入变量。

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