技术领域
[0001] 本
发明涉及一种基于
土壤发生层厚度预测的土壤有机碳储量估算方法，属于计量土壤技术领域。
背景技术
[0002] 土壤碳主要包括土壤有机碳和土壤无机碳，土壤无机碳库(碳酸盐碳)较为稳定。土壤有机碳主要分布在1米深度土体内，对土壤的物理、化学性质有重要影响,并直接影响着土壤
质量。作为土壤肥
力评价的重要指标之一，土壤有机碳的不断降低会直接导致农业耕作土壤的贫瘠问题。土壤有机碳储量不仅是生态公益调查重点关注的社会问题之一，也是全球性的
基础课题。国内外已有众多应用生产、环境保护部
门开始对人类生存发展面临的全球
气候增暖问题给予广泛关注。相关国际条例也对国家级碳收支计算技术提出了不同需求，土壤有机碳储量估算已经成为影响国家经济、外交的重大生态与环境生产技术问题。
[0003] 土壤有机碳
密度指代一定面积固定土壤深度的碳储量，单位是kg/m2，一定区域的土壤有机碳储量是有机碳密度与区域面积的乘积，单位是kg。常规的碳储量估算技术主要包括土壤类型法、植被类型法、生命带类型法、模型法、相关关系统计方法、GIS(
地理信息系统)空间预测法。基于GIS的空间预测方法是现代数字土壤制图广泛采用的一种空间预测模式。有别于传统的土壤调查和制图技术，该方法采用国内外土壤科学家广泛肯定的量化土壤-景观模型，能够有效地集成遥感
图像处理技术、数字地形分析技术、GIS空间分析技术和土壤调查技术，通过对景观信息的分析来预测土壤有机碳密度的空间分布。GIS空间预测方法的主要作业流程是通过不同地理
位置的野外
采样、实验室分析土壤理化属性、基于环境变量与已知样点土壤属性数据建立碳储量估算模型，进而估算目标区域的有机碳储量。
[0004] 目前，常用的有机碳储量估算技术主要基于以下两种生产计算模式：先计算再预测(CTM)、先预测再计算(MTC)，这里的“计算”是计算单个土体的土壤有机碳密度(基于土壤有机碳含量、土壤容重、土壤砾石含量、土壤厚度)，“预测”是空间插值预测目标土壤属性的空间分布。
现有技术应用主要采用一种计算模式，对比应用模式尚鲜有资料可参阅。
[0005] 近年来，随着计算机、遥感、土壤调查、测绘技术的发展，特别是高
分辨率遥感、数字测绘技术的快速发展，区域级别的环境变量信息包含被测更多异质性细节的数据获取成为可能，由空间离散的土壤样点数据，定量、客观、实时、准确地模拟复杂景观区域土壤属性空间分布给土壤信息相关部门的具体生产提出了更高的要求，这也发展成为国际高
精度土壤有机碳储量估算的重要发展方向。
[0006] 发生层是
鉴别土壤类型的重要依据，并在性质上有一系列的定量说明。作为碳储量估算的预测基础，土壤-景观模型假设土壤属性的实际分布状况(包括
水平维与垂直维)是与景观属性密切相关的。尤其是在土壤属性的垂直维分布方面，发生层往往具有均质的土壤属性特征。然而，现有土壤有机碳储量估算模型均面向固定深度的土壤属性预测。全球土壤数字制图计划(GSM)联合协议规定专题土壤图的生产土壤层次厚度均采用0-5，5-15，15-30,30-60,60-100,100-200cm的固定方式。
[0007] 目前，土壤有机碳储量估算在复杂景观区域和生产技术上存在一定的局限性，具体归纳起来有以下几点：
[0008] (1)固定层次厚度的计算模式从某种程度忽略了土壤发生的理论模型，在生产加工上具有一定的局限性。由于地形、气候、人为干扰的原因土壤垂直维空间变异的复杂程度远远超出了现有计算技术的模拟能力。如果土壤采样点的
属性信息缺失，尤其是样点总量较少且样点的空间代表性较差时，实际应用中很难准确估算区域的土壤有机碳储量。事实上，这也是造成现有数字土壤制图表达存在严重区域性不确定现象的主要原因之一。中国科学院南京土壤研究所与ISRIC——世界土壤信息参比中心(荷兰)均已指出，现代数字土壤制图技术的应用
框架不应忽略发生层的土壤信息特征。
[0009] (2)在实际土壤信息生产应用过程中，现有生产技术很难综合考虑土壤发生层厚度。现有土壤调查工程仍采用传统的野外调查模式，不仅无法实时获取大规模区域的土壤理化属性数据与描述性信息，也往往受限于调查经费支持与人员的调查作业的技术能力。土壤系统分类涉及诸多的土壤诊断特性，不同的诊断特性具有定量与定性的土壤发生模型理论，给技术人员的实际生产带来了很大困难。因此，现有基础性土壤有机碳储量估算工程往往缺乏综合考虑土壤发生层特征的归并技术。
[0010] (3)现有技术多集中于一种固定计算模式，缺乏强化对比研究。国内外大量应用案例已表明，不同的生产模式(CTM、MTC)往往在不同的模型假设、景观突变区域表现出迥异的区域性精度问题。因此，强化对比不同的计算模式对于具体的生产环节具有重要的指导意义。同时，土壤数据非结构化的数据结构与发生层的计量土壤特征是导致现有生产技术不完善的重要原因。
[0011] 以上所述现有土壤有机碳储量估算技术的不足，在不同应用部门生产加工土壤信息产品和工程应用中带来较大困难，在诸如偏远山区、人为影响较复杂区域、高度变异林区等区域的碳储量估算工程应用上会带来潜在的错误决策支持，进而给国家经济规划直接造成损失。
发明内容
[0012] 本发明所要解决的技术问题是提供一种采用全新架构设计，能够有效提高土壤有机碳储量估算精度与估算工作效率的基于土壤发生层厚度预测的土壤有机碳储量估算方法。
[0013] 本发明为了解决上述技术问题采用以下技术方案：本发明设计了一种基于土壤发生层厚度预测的土壤有机碳储量估算方法，包括如下步骤：
[0014] 步骤001.针对目标土壤区域，设置各个样点位置，并统计样点位置的数量为I，然后进入步骤002；
[0015] 步骤002.在目标土壤区域中，分别针对各个样点位置，获得样点位置的环境信息，同时若该样点位置竖直向下土壤区域的深度大于等于预设深度L，则获得该样点位置竖直向下预设深度L的剖面区域；若该样点位置竖直向下土壤区域的深度小于预设深度L，则获得该样点位置竖直向下土壤的剖面区域；并根据中国土壤系统分类的作业要求，针对该剖面区域划分土壤发生层，获得该样点位置所对应的各个土壤发生层；再分别获得该样点位置所对应各个土壤发生层的厚度信息和土壤形态信息；进而分别获得目标土壤区域中各个样点位置的环境信息，以及各个样点位置分别所对应各个土壤发生层的厚度信息和土壤形态信息，接着进入步骤003；
[0016] 步骤003.分别针对目标土壤区域中的各个样点位置，在样点位置所对应的各个土壤发生层中，分别采集预设质量的土壤样品，并分别测定获得各个土壤发生层的土壤属性信息，进而获得目标土壤区域中各个样点位置分别所对应各个土壤发生层的土壤属性信息，然后进入步骤004；其中，土壤属性信息包括土壤有机碳含量、大于预设直径的砾石含量和土壤容重；
[0017] 步骤004.分别针对目标土壤区域中的各个样点位置，根据样点位置所对应各个土壤发生层的土壤属性信息、厚度信息，获得各个样点位置的土壤有机碳密度实测值SOCDi，i＝{1,…,I}，然后进入步骤005；
[0018] 步骤005.分别针对目标土壤区域中的各个样点位置，针对样点位置所对应的各个土壤发生层进行归并处理，获得该样点位置所对应的各个归并层，然后根据该样点位置各个归并层所分别对应的各个土壤发生层，针对该样点位置各个土壤发生层的土壤属性信息进行加权计算，获得该样点位置所对应各个归并层的土壤属性信息，以及针对该样点位置各个土壤发生层的厚度信息进行求和计算，获得该样点位置所对应各个归并层的厚度；进而分别获得目标土壤区域中各个样点位置分别对应的各个归并层，以及各个归并层的土壤属性信息和厚度信息，再进入步骤006；
[0019] 步骤006.分别针对目标土壤区域中的各个样点位置，判断样点位置所对应各个归并层的厚度之和是否小于预设深度L，是则在该样点位置所对应各个归并层之下，以基岩层设置归并层，使得该样点位置所对应各个归并层的厚度之和等于预设深度L，并进入步骤007；否则直接进入步骤007；
[0020] 步骤007.获取目标土壤区域中所有样点位置所对应归并层的种类，构成目标土壤区域归并层种类集合；根据目标土壤区域归并层种类集合，针对目标土壤区域中各个样点位置所对应的归并层进行统一的操作，使得目标土壤区域中各个样点位置所对应归并层的种类彼此相同，并进入步骤008；
[0021] 步骤008.采用线性同余
算法，将目标土壤区域中的样点位置划分为预测样点位置集合和验证样点位置集合，并根据目标土壤区域中各个样点位置的土壤有机碳密度实测值SOCDi，获得验证样点位置集合中各个验证样点位置的土壤有机碳密度实测值 构成验证样点位置集合中各个验证样点位置的土壤有机碳密度实测值集合V；并进入步骤009，其中，i2＝{1,…,I2}，其中，I2为验证样点位置集合中验证样点位置的数量，预测样点位置集合中预测样点位置的数量为I1，I1>I2；
[0022] 步骤009.根据预测样点位置集合中各个预测样点位置的环境信息，以及各个预测样点位置分别所对应各个归并层的厚度信息，采用
随机森林方法，训练获得分别以各个归并层厚度信息为目标的各个
预测模型，构成第一预测模型集合；
[0023] 同时，根据预测样点位置集合中各个预测样点位置的环境信息，以及各个预测样点位置分别所对应各个归并层的土壤属性信息，采用随机森林方法，训练获得分别以各个归并层土壤属性信息为目标的各个预测模型，构成第二预测模型集合；然后进入步骤010；
[0024] 步骤010.分别针对验证样点位置集合中的各个验证样点位置，根据验证样点位置的环境信息，通过第一预测模型集合中的各个预测模型，分别获得该验证样点位置分别所对应各个归并层的第一预测厚度信息；同时，分别针对验证样点位置集合中的各个验证样点位置，根据验证样点位置的环境信息，通过第二预测模型集合中的各个预测模型，分别获得该验证样点位置所分别对应各个归并层的第一预测土壤属性信息；进而获得验证样点位置集合中各个验证样点位置分别所对应各个归并层的第一预测厚度信息和第一预测土壤属性信息；然后进入步骤011；其中，第一预测土壤属性信息包括第一有机碳含量预测信息、第一土壤容重预测信息和大于预设直径的第一砾石预测含量；
[0025] 步骤011.分别针对验证样点位置集合中的各个验证样点位置，获得L相对验证样点位置所对应各个归并层第一预测厚度信息之和的拉伸系数Rc；然后根据验证样点位置集合中各个验证样点位置分别所对应各个归并层的第一预测土壤属性信息、第一预测厚度信息，获得验证样点位置集合中各个验证样点位置的第一土壤有机碳密度预测值 构成验证样点位置集合中各个验证样点位置的第一土壤有机碳密度预测值集合PMTC-D；然后进入步骤012；
[0026] 步骤012.根据预测样点位置集合中各个预测样点位置的环境信息，以及各个预测样点位置的土壤有机碳密度实测值，采用随机森林方法，训练获得以土壤有机碳密度实测值为目标的第三预测模型；然后分别针对验证样点位置集合中的各个验证样点位置，根据验证样点位置的环境信息，通过第三预测模型，获得该验证样点位置的第二土壤有机碳密度预测值，进而获得验证样点位置集合中各个验证样点位置的第二土壤有机碳密度预测值，构成验证样点位置集合中各个验证样点位置的第二土壤有机碳密度预测值集合PCTM；然后进入步骤013；
[0027] 步骤013.针对目标土壤区域中所有样点位置，统一按预设划分规则，基于L将样点位置所对应的各个归并层划分为各个拟合层，各个样点位置所对应的拟合层、以及拟合层的数量相同，并获得各个拟合层的厚度信息；然后分别针对预测样点位置集合中的各个预测样点位置，根据预测样点位置各个拟合层所分别对应的各个归并层，针对该预测样点位置各个归并层采样发生层的土壤属性信息进行拟合，获得预测样点位置各拟合层的土壤属性信息；进而获得预测样点位置集合中各个预测样点位置分别所对应各个拟合层的土壤属性信息，并进入步骤014；
[0028] 步骤014.根据预测样点位置集合中各个预测样点位置的环境信息，以及各个预测样点位置分别所对应各个拟合层的土壤属性信息，采用随机森林方法，训练获得以各个拟合层土壤属性信息为目标的各个预测模型，构成第四预测模型集合；然后分别针对验证样点位置集合中的各个验证样点位置，根据验证样点位置的环境信息，通过第四预测模型集合中的各个预测模型，获得该验证样点位置分别所对应各个拟合层的第三土壤属性预测信息，进而获得验证样点位置集合中各个验证样点位置分别所对应各个拟合层的第三土壤属性预测信息，然后进入步骤015；其中，第三预测土壤属性信息包括第三有机碳含量预测信息、第三土壤容重预测信息和大于预设直径的第三砾石预测含量；
[0029] 步骤015.根据验证样点位置集合中各个验证样点位置所对应各个拟合层的第三土壤属性预测信息，以及厚度信息，获得验证样点位置集合中各个验证样点位置的第三土壤有机碳密度预测值 构成验证样点位置集合中各个验证样点位置的第三土壤有机碳密度预测值集合PMTC-F；然后进入步骤016；
[0030] 步骤016.根据土壤有机碳密度实测值集合V进行精度检验，获得第一土壤有机碳密度预测值集合PMTC-D、第二土壤有机碳密度预测值集合PCTM和第三土壤有机碳密度预测值集合PMTC-F中的最优土壤有机碳密度预测值集合，并进入步骤017；
[0031] 步骤017.将目标土壤区域离散化空间栅格数据，将目标土壤区域中所有样点位置分别所对应发生层的采样数据作为预测数据集合，若最优土壤有机碳密度预测值集合为第一土壤有机碳密度预测值集合PMTC-D，则采用步骤009至步骤011的方法获得目标土壤区域的土壤有机碳密度空间分布栅格数据；若最优土壤有机碳密度预测值集合为第二土壤有机碳密度预测值集合PCTM，则采用步骤012的方法获得目标土壤区域的土壤有机碳密度空间分布栅格数据；若最优土壤有机碳密度预测值集合为第三土壤有机碳密度预测值集合PMTC-F，则采用步骤013至步骤015的方法获得目标土壤区域的土壤有机碳密度空间分布栅格数据；然后进入步骤018；
[0032] 步骤018.根据目标土壤区域的土壤有机碳密度空间分布栅格数据，获得目标土壤区域的土壤有机碳储量。
[0033] 作为本发明的一种优选技术方案，所述步骤001具体包括如下步骤：
[0034] 步骤00101.获得目标土壤区域的土壤类型分布图、土壤利用分布图和土壤地质分布图，并将土壤类型分布图、土壤利用分布图和土壤地质分布图进行空间
叠加，获得目标土壤区多层叠加图，然后进入步骤00102；
[0035] 步骤00102.分别获得目标土壤区多层叠加图中各个图斑区域的面积，并统计获得面积比例超过图斑面积
阈值的各个图斑区域；然后分别针对该各个图斑区域，根据实地可达性检测要求，设置样点位置，统计样点位置的数量为I，再进入步骤002。
[0036] 作为本发明的一种优选技术方案：所述环境信息包括植被
覆盖率、
岩石露头面积比、地形、重要标志物、地表粗碎
块大小、地表裂隙状况、地表盐斑信息。
[0037] 作为本发明的一种优选技术方案：所述土壤形态信息表示土壤干湿状况、土壤
颜色、根系信息、孔隙信息、样品结构、斑纹组成物质、瘤状结核物质、胶结程度、石灰反应信息。
[0038] 作为本发明的一种优选技术方案，所述步骤004具体包括如下操作：
[0039] 分别针对目标土壤区域中的各个样点位置，根据样点位置所对应各个土壤发生层的土壤属性信息、厚度信息，按如下公式：
[0040]
[0041] 获得各个样点位置的土壤有机碳密度实测值SOCDi，其中，i＝{1,…,I}，ni＝{1,…Ni}，ni表示第i个样点位置所对应的第n个土壤发生层，Ni表示第i个样点位置所对应土壤发生层的总数； 表示第i个样点位置所对应第n个土壤发生层的有机碳含量， 表示第i个样点位置所对应第n个土壤发生层的土壤容重， 表示第i个样点位置所对应第n个土壤发生层中大于预设直径的砾石含量， 表示第i个样点位置所对应第n个土壤发生层的厚度信息，然后进入步骤005。
[0042] 作为本发明的一种优选技术方案，所述步骤005中，分别针对目标土壤区域中的各个样点位置，针对样点位置所对应的各个土壤发生层进行归并处理，获得该样点位置所对应的各个归并层的过程，具体包括如下步骤：
[0043] 步骤00501.分别针对目标土壤区域中的各个样点位置，根据样点位置所对应各个土壤发生层的土壤形态信息，以及中国土壤系统分类，获得该样点位置的各个土壤发生层分别所对应的特性表达符号；进而获得目标土壤区域中各个样点位置的各个土壤发生层分别所对应的特性表达符号，然后进入步骤00502；
[0044] 步骤00502.根据中国土壤系统分类，建立发生层特性归并表，如下表1所示：
[0045]
[0046]
[0047] 表1
[0048] 其中，阿拉伯数字为
土壤特性表达符号相同，且存在进一步划分的各个土壤发生层自土壤剖面由自上至下顺序的排列；
[0049] 然后分别针对目标土壤区域中的各个样点位置，针对样点位置所对应各个土壤发生层，按
腐殖质表层、淀积层、母质层进行划分，并按发生层特性归并表中各归并层次所对应的土壤发生层特性表达符号，自该样点位置剖面区域由上而下，依次划分该样点位置所对应的各个土壤发生层，获得该样点位置所对应的各个归并层；
[0050] 与此同时，在分别针对目标土壤区域中各个样点位置进行上述归并操作的过程中，若样点位置所对应两个归并层次中出现相同的土壤发生层特性表达符号时，则进入步骤00503；
[0051] 步骤00503.判断该样点位置所对应的各个土壤发生层特性表达符号中，是否存在位于发生层特性归并表中由上向下第一个归并层次中的土壤发生层特性表达符号，是则将该相同土壤发生层特性表达符号分别所对应的土壤发生层进行合并，然后再执行步骤00502进行归并处理，获得该样点位置所分别对应的各个归并层；否则进入步骤00504；
[0052] 步骤00504.将该相同土壤发生层特性表达符号中，位于下层归并层中的土壤发生层特性表达符号，划分至与其相同的土壤发生层特性表达符号所在的归并层中，并将该样点位置所对应的各个土壤发生层特性表达符号，按发生层特性归并表中的归并层次顺序向上移动一个归并层次，获得该样点位置所分别对应的各个归并层。
[0053] 作为本发明的一种优选技术方案，所述步骤007中，根据目标土壤区域归并层种类集合，针对目标土壤区域中各个样点位置所对应的归并层进行统一的操作，使得目标土壤区域中各个样点位置所对应归并层的种类彼此相同，具体包括如下操作：
[0054] 分别针对目标土壤区域中的各个样点位置，判断样点位置所对应归并层是否等于目标土壤区域归并层种类集合，是则不做任何操作；否则向该样点位置设置其相对目标土壤区域归并层种类集合缺少的归并层，同时设置该归并层的厚度信息为0.00001cm，以及设置该归并层的土壤属性信息为0.00001；进而使得目标土壤区域中各个样点位置分别所对应的归并层均等于目标土壤区域归并层种类集合。
[0055] 作为本发明的一种优选技术方案，所述步骤011具体包括如下操作：
[0056] 分别针对验证样点位置集合中的各个验证样点位置，获得L相对验证样点位置所对应各个归并层第一预测厚度信息之和的拉伸系数Rc；然后根据验证样点位置集合中各个验证样点位置分别所对应各个归并层的第一预测土壤属性信息、第一预测厚度信息，按如下公式：
[0057]
[0058] 获得验证样点位置集合中各个验证样点位置的第一土壤有机碳密度预测值构成验证样点位置集合中各个验证样点位置的第一土壤有机碳密度预测值集合PMTC-D；其中， 表示验证样点位置集合中第i2个验证样点位置所对应归并层的总数， 表示验证样点位置集合中第i2个验证样点位置所对应的第n个归并层，表示验证样点位置集合中第i2个验证样点位置所对应第n个归并层的第一有机碳含量预测信息， 表示验证样点位置集合中第i2个验证样点位置所对应第n个归并层的第一土壤容重预测信息， 表示验证样点位置集合中第i2个验证样点位置所对应第n个归并层中大于预设直径的第一砾石预测含量， 表示验证样点位置集合中第i2个验证样点位置所对应第n个归并层的第一预测厚度信息；然后进入步骤012。
[0059] 作为本发明的一种优选技术方案，所述步骤013中，获得预测样点位置集合中各个预测样点位置分别所对应各个拟合层的土壤属性信息，具体包括如下操作：
[0060] 分别针对预测样点位置集合中的各个预测样点位置，根据等面积Spline函数，采用如下公式：
[0061]
[0062] 获得该预测样点位置各个拟合层的土壤属性信息；进而获得预测样点位置集合中各个预测样点位置分别所对应各个拟合层的土壤属性信息，其中，i1＝{1,…,I1}，I1表示预测样点位置集合中预测样点位置的数量， 表示预测样点位置集合中第i1个预测样点位置所对应拟合层的层数， 表示预测样点位置集合中第i1个预测样点位置所对应第k层拟合层的土壤属性信息， 表示 函数在 层与 层拟合结果的平均值， 表示预测样点位置集合中第i1个预测样点位置所对应第k层拟合层的实验室测定分析误差， 函数是等面积Spline函数。
[0063] 作为本发明的一种优选技术方案，所述步骤015具体包括如下操作：
[0064] 根据验证样点位置集合中各个验证样点位置所对应各个拟合层的第三土壤属性预测信息，以及厚度信息，按如下公式：
[0065]
[0066] 获得验证样点位置集合中各个验证样点位置的第三土壤有机碳密度预测值构成验证样点位置集合中各个验证样点位置的第三土壤有机碳密度预测值集合PMTC-F；其中， 表示验证样点位置集合中第i2个验证样点位置所对应拟合层的总数， 表示验证样点位置集合中第i2个验证样点位置所对应的第k个拟合层，表示验证样点位置集合中第i2个验证样点位置所对应第k个拟合层的第三有机碳含量预测信息， 表示验证样点位置集合中第i2个验证样点位置所对应第k个拟合层的第三土壤容重预测信息， 表示验证样点位置集合中第i2个验证样点位置所对应第k个拟合层中大于预设直径的第三砾石预测含量， 表示验证样点位置集合中第i2个验证样点位置所对应第k个拟合层的厚度信息；然后进入步骤016。
[0067] 本发明所述一种基于土壤发生层厚度预测的土壤有机碳储量估算方法采用以上技术方案与现有技术相比，具有以下技术效果：
[0068] (1)本发明所设计一种基于土壤发生层厚度预测的土壤有机碳储量估算方法，通过针对发生层非结构化信息的封装，为考虑土壤属性水平维空间分布的土体连续性预测提供了较好的技术思路；其中，采用“发生层归并，预测再计算”技术，在保证土壤发生层特性信息不缺失的同时，修正了传统预测方法在土体连续性描述的局限性，实现了“描述规范，预测准确”的通用土壤有机碳储量估算技术，在农业应用、环境保护、国土资源等相关部门的工程调查方面具有广阔的工业化应用前景；
[0069] (2)本发明所设计一种基于土壤发生层厚度预测的土壤有机碳储量估算方法中，所提出的发生层厚度预测及其拉伸系数计算具有一定的普适性，其技术方案不仅面向土壤有机碳密度预测，也可以与发生层类型预测相结合，构成土壤类型预测的技术流程，因此，本发明提出的预测技术还具有较好的
稳定性，动态构建的发生层符号归并机制有望在保证预测精度的同时，也降低了土壤发生层厚度预测过程的误差；
[0070] (3)本发明所设计一种基于土壤发生层厚度预测的土壤有机碳储量估算方法中，所提出的发生层归并技术能够最大程度上简化土壤信息的结构化应用，有望为其他的工程应用，如为三维土壤制图、全球数字土壤制图及区域碳储量
时空演化分析提供技术指导，并且所使用的强化对比模式兼顾了不同计算模式的优点，确保了能够使用最优的计算模式对土壤有机碳密度进行空间预测，进而能够更加定量、客观评估目标区域的有机碳储量。
附图说明
[0071] 图1是本发明设计基于土壤发生层厚度预测的土壤有机碳储量估算方法的流程示意图；
[0072] 图2a是本发明所应用
实施例中当土壤深度大于等于预设深度1m时，采集1m深度的土壤剖面区域的示意图；
[0073] 图2b是本发明所应用实施例中当土壤深度小于预设深度1m时，采集竖直向下土壤的剖面区域的示意图；
[0074] 图3是本发明实施例中不同土壤发生层厚度空间变化的示意图；
[0075] 图4是本发明实施例中针对1m深度剖面区域划分获得固定层次厚度归并层的示意图；
[0076] 图5a至图5i是本发明实施例中各个归并层的预测厚度的空间分布图；
[0077] 图6是本发明实施例中1m深度剖面区域相对预测厚度的拉伸系数分布示意图；
[0078] 图7a至图7c是本发明实施例中不同方法预测结果的对比分析示意图。
具体实施方式
[0079] 下面结合
说明书附图对本发明的具体实施方式作进一步详细的说明。
[0080] 针对上述技术
缺陷，本发明以土壤演化的专家知识为指导，设计了一套新的基于土壤发生层厚度预测的土壤有机碳储量估算工程技术方案。该技术系统地给出了有机碳储量估算的预测的技术流程。同时，考虑土壤有机碳密度的高度空间变异特性，本发明创新地提出对比使用不同的碳储量估算计算模式。强化对比策略有望提供传统基于GIS估算技术所忽略的最优计算模式。
[0081] 本发明的基本思想是在土壤调查采样过程中，系统记录不同土体的发生层特性与描述信息。在实验室土壤属性化验分析后，使用结构化的土壤属性数据结构输入非结构化的土壤信息。以发生层特性表达符号为操作单元，构建具有相似发生层与相等发生层数量的标准化土壤剖面数据。使用
机器学习算法预测目标点位的发生层厚度与土壤属性，基于预测发生层厚度的拉伸系数，更新土壤有机碳密度的计算。通过对比分析其他的计算模式，遴选出最优的计算模式，最终估算目标区域的有机碳储量。
[0082] 图1所示，本发明所设计的一种基于土壤发生层厚度预测的土壤有机碳储量估算方法，在实际应用过程当中，具体包括如下步骤：
[0083] 步骤001.针对目标土壤区域，设置各个样点位置，并统计样点位置的数量为I，然后进入步骤002。
[0084] 其中，步骤001具体包括如下步骤：
[0085] 步骤00101.获得目标土壤区域的土壤类型分布图、土壤利用分布图和土壤地质分布图，并将土壤类型分布图、土壤利用分布图和土壤地质分布图进行空间叠加，获得目标土壤区多层叠加图，然后进入步骤00102。
[0086] 步骤00102.分别获得目标土壤区多层叠加图中各个图斑区域的面积，并统计获得面积比例超过图斑面积阈值的各个图斑区域；然后分别针对该各个图斑区域，根据实地可达性检测要求，设置样点位置，统计样点位置的数量为I，再进入步骤002。
[0087] 步骤002.在目标土壤区域中，分别针对各个样点位置，获得样点位置的环境信息，环境信息包括植被覆盖率、岩石露头面积比、地形、重要标志物、地表粗碎块大小、地表裂隙状况、地表盐斑信息；同时，如图2a所示，若该样点位置竖直向下土壤区域的深度大于等于预设深度L，则获得该样点位置竖直向下预设深度L的剖面区域；如图2b所示，若该样点位置竖直向下土壤区域的深度小于预设深度L，则获得该样点位置竖直向下土壤的剖面区域；并根据中国土壤系统分类的作业要求，针对该剖面区域划分土壤发生层，获得该样点位置所对应的各个土壤发生层；再分别获得该样点位置所对应各个土壤发生层的厚度信息和土壤形态信息，土壤形态信息表示土壤干湿状况、土壤颜色、根系信息、孔隙信息、样品结构、斑纹组成物质、瘤状结核物质、胶结程度、石灰反应信息；进而分别获得目标土壤区域中各个样点位置的环境信息，以及各个样点位置分别所对应各个土壤发生层的厚度信息和土壤形态信息，接着进入步骤003。
[0088] 其中，中国土壤系统分类表如下表2所示：
[0089]发生层符号 发生层指代信息描述 发生层符号 发生层指代信息描述
A 腐殖质表层 l 网纹
B 淀积层 m 强胶结
C 母质层 n 钠聚积
R 基岩 o 根系盘结
a 高分解有机物质 p 耕作影响
b 埋藏层 q 次生
硅聚积
c 结皮 r
氧化还原
d 冻融性片状结构 s
铁锰聚积
e 半分解有机物质 t 黏粒聚积
f 永冻层 u 人为堆积影响
g 潜育斑 v 变性特征
h 腐殖质聚积 w 就地
风化形成的显色、有结构层
i 低分解和未分解有机物质 x 固态坚硬的胶结，未形成磐
j 黄
钾铁矾 y
石膏聚积
k 碳酸盐聚积 z 可溶盐聚积
[0090] 表2
[0091] 步骤003.分别针对目标土壤区域中的各个样点位置，在样点位置所对应的各个土壤发生层中，分别采集预设质量的土壤样品，并分别测定获得各个土壤发生层的土壤属性信息，进而获得目标土壤区域中各个样点位置分别所对应各个土壤发生层的土壤属性信息，然后进入步骤004；其中，土壤属性信息包括土壤有机碳含量(单位：g/kg)、大于预设直径0.2mm的砾石含量(体积百分比，％)和土壤容重(单位：g/cm3)。
[0092] 步骤004.分别针对目标土壤区域中的各个样点位置，根据样点位置所对应各个土壤发生层的土壤属性信息、厚度信息，获得各个样点位置的土壤有机碳密度实测值SOCDi，i＝{1,…,I}，然后进入步骤005。
[0093] 所述步骤004具体包括如下操作：
[0094] 分别针对目标土壤区域中的各个样点位置，根据样点位置所对应各个土壤发生层的土壤属性信息、厚度信息，按如下公式：
[0095]
[0096] 获得各个样点位置的土壤有机碳密度实测值SOCDi，其中，i＝{1,…,I}，ni＝{1,…Ni}，ni表示第i个样点位置所对应的第n个土壤发生层，Ni表示第i个样点位置所对应土壤发生层的总数； 表示第i个样点位置所对应第n个土壤发生层的有机碳含量， 表示第i个样点位置所对应第n个土壤发生层的土壤容重， 表示第i个样点位置所对应第n个土壤发生层中大于预设直径的砾石含量， 表示第i个样点位置所对应第n个土壤发生层的厚度信息，然后进入步骤005。
[0097] 步骤005.分别针对目标土壤区域中的各个样点位置，针对样点位置所对应的各个土壤发生层进行归并处理，获得该样点位置所对应的各个归并层；然后根据该样点位置各个归并层所分别对应的各个土壤发生层，针对该样点位置各个土壤发生层的土壤属性信息进行加权计算，获得该样点位置所对应各个归并层的土壤属性信息，以及针对该样点位置各个土壤发生层的厚度信息进行求和计算，获得该样点位置所对应各个归并层的厚度；进而分别获得目标土壤区域中各个样点位置分别对应的各个归并层，以及各个归并层的土壤属性信息和厚度信息，再进入步骤006。
[0098] 步骤005中，分别针对目标土壤区域中的各个样点位置，针对样点位置所对应的各个土壤发生层进行归并处理，获得该样点位置所对应的各个归并层的过程，具体包括如下步骤：
[0099] 步骤00501.分别针对目标土壤区域中的各个样点位置，根据样点位置所对应各个土壤发生层的土壤形态信息，以及中国土壤系统分类，获得该样点位置的各个土壤发生层分别所对应的特性表达符号；进而获得目标土壤区域中各个样点位置的各个土壤发生层分别所对应的特性表达符号，然后进入步骤00502。
[0100] 步骤00502.根据中国土壤系统分类，建立发生层特性归并表，如下表1所示：
[0101]
[0102]
[0103] 表1
[0104] 其中，阿拉伯数字为土壤特性表达符号相同，且存在进一步划分的各个土壤发生层自土壤剖面由自上至下顺序的排列；注：(1)主要发生层或特性发生层可按其发生程度上的差异进一步细分为若干亚层。均以阿拉伯数字与大写字母并列表示，例如C1，C2，Bt1，Bt2，Bt3；(2)异元母质土层表示：用阿拉伯数字置于发生层符号前表示，如2C。
[0105] 然后分别针对目标土壤区域中的各个样点位置，针对样点位置所对应各个土壤发生层，按腐殖质表层、淀积层、母质层进行划分，并按发生层特性归并表中各归并层次所对应的土壤发生层特性表达符号，自该样点位置剖面区域由上而下，依次划分该样点位置所对应的各个土壤发生层，获得该样点位置所对应的各个归并层。
[0106] 与此同时，在分别针对目标土壤区域中各个样点位置进行上述归并操作的过程中，若样点位置所对应两个归并层次中出现相同的土壤发生层特性表达符号时，则进入步骤00503。
[0107] 步骤00503.判断该样点位置所对应的各个土壤发生层特性表达符号中，是否存在位于发生层特性归并表中由上向下第一个归并层次中的土壤发生层特性表达符号，是则将该相同土壤发生层特性表达符号分别所对应的土壤发生层进行合并，然后再执行步骤00502进行归并处理，获得该样点位置所分别对应的各个归并层；否则进入步骤00504。
[0108] 步骤00504.将该相同土壤发生层特性表达符号中，位于下层归并层中的土壤发生层特性表达符号，划分至与其相同的土壤发生层特性表达符号所在的归并层中，并将该样点位置所对应的各个土壤发生层特性表达符号，按发生层特性归并表中的归并层次顺序向上移动一个归并层次，获得该样点位置所分别对应的各个归并层。
[0109] 步骤006.如图3所示，分别针对目标土壤区域中的各个样点位置，判断样点位置所对应各个归并层的厚度之和是否小于预设深度L，是则在该样点位置所对应各个归并层之下，以基岩层设置归并层，使得该样点位置所对应各个归并层的厚度之和等于预设深度L，并进入步骤007；否则直接进入步骤007。
[0110] 步骤007.获取目标土壤区域中所有样点位置所对应归并层的种类，构成目标土壤区域归并层种类集合；根据目标土壤区域归并层种类集合，针对目标土壤区域中各个样点位置所对应的归并层进行统一的操作，使得目标土壤区域中各个样点位置所对应归并层的种类彼此相同，并进入步骤008。
[0111] 步骤007中，根据目标土壤区域归并层种类集合，针对目标土壤区域中各个样点位置所对应的归并层进行统一的操作，使得目标土壤区域中各个样点位置所对应归并层的种类彼此相同，具体包括如下操作：
[0112] 分别针对目标土壤区域中的各个样点位置，判断样点位置所对应归并层是否等于目标土壤区域归并层种类集合，是则不做任何操作；否则向该样点位置设置其相对目标土壤区域归并层种类集合缺少的归并层，同时设置该归并层的厚度信息为0.00001cm，以及设置该归并层的土壤属性信息为0.00001；进而使得目标土壤区域中各个样点位置分别所对应的归并层均等于目标土壤区域归并层种类集合。
[0113] 步骤008.采用线性同余算法，将目标土壤区域中的样点位置划分为预测样点位置集合和验证样点位置集合，并根据目标土壤区域中各个样点位置的土壤有机碳密度实测值SOCDi，获得验证样点位置集合中各个验证样点位置的土壤有机碳密度实测值 构成验证样点位置集合中各个验证样点位置的土壤有机碳密度实测值集合V；并进入步骤009，其中，i2＝{1,…,I2}，其中，I2为验证样点位置集合中验证样点位置的数量，等于25％I，预测样点位置集合中预测样点位置的数量为I1，等于75％I。
[0114] 步骤009.根据预测样点位置集合中各个预测样点位置的环境信息，以及各个预测样点位置分别所对应各个归并层的厚度信息，采用随机森林方法，训练获得分别以各个归并层厚度信息为目标的各个预测模型，构成第一预测模型集合。
[0115] 同时，根据预测样点位置集合中各个预测样点位置的环境信息，以及各个预测样点位置分别所对应各个归并层的土壤属性信息，采用随机森林方法，训练获得分别以各个归并层土壤属性信息为目标的各个预测模型，构成第二预测模型集合；然后进入步骤010。
[0116] 步骤010.采用MTC-D方法，分别针对验证样点位置集合中的各个验证样点位置，根据验证样点位置的环境信息，通过第一预测模型集合中的各个预测模型，分别获得该验证样点位置分别所对应各个归并层的第一预测厚度信息，如图5a至图5i所示；同时，分别针对验证样点位置集合中的各个验证样点位置，根据验证样点位置的环境信息，通过第二预测模型集合中的各个预测模型，分别获得该验证样点位置所分别对应各个归并层的第一预测土壤属性信息；进而获得验证样点位置集合中各个验证样点位置分别所对应各个归并层的第一预测厚度信息和第一预测土壤属性信息；然后进入步骤011；其中，第一预测土壤属性信息包括第一有机碳含量预测信息、第一土壤容重预测信息和大于预设直径的第一砾石预测含量。
[0117] 步骤011.分别针对验证样点位置集合中的各个验证样点位置，获得L相对验证样点位置所对应各个归并层第一预测厚度信息之和的拉伸系数Rc，如图6所示；然后根据验证样点位置集合中各个验证样点位置分别所对应各个归并层的第一预测土壤属性信息、第一预测厚度信息，获得验证样点位置集合中各个验证样点位置的第一土壤有机碳密度预测值构成验证样点位置集合中各个验证样点位置的第一土壤有机碳密度预测值集合PMTC-D；然后进入步骤012。
[0118] 其中，步骤011具体包括如下操作：
[0119] 分别针对验证样点位置集合中的各个验证样点位置，获得L相对验证样点位置所对应各个归并层第一预测厚度信息之和的拉伸系数Rc；然后根据验证样点位置集合中各个验证样点位置分别所对应各个归并层的第一预测土壤属性信息、第一预测厚度信息，按如下公式：
[0120]
[0121] 获得验证样点位置集合中各个验证样点位置的第一土壤有机碳密度预测值构成验证样点位置集合中各个验证样点位置的第一土壤有机碳密度预测值集合PMTC-D；其中， 表示验证样点位置集合中第i2个验证样点位置所对应归并层的总数， 表示验证样点位置集合中第i2个验证样点位置所对应的第n个归并层，表示验证样点位置集合中第i2个验证样点位置所对应第n个归并层的第一有机碳含量预测信息， 表示验证样点位置集合中第i2个验证样点位置所对应第n个归并层的第一土壤容重预测信息， 表示验证样点位置集合中第i2个验证样点位置所对应第n个归并层中大于预设直径的第一砾石预测含量， 表示验证样点位置集合中第i2个验证样点位置所对应第n个归并层的第一预测厚度信息；然后进入步骤012。
[0122] 步骤012.采用CTM方法，根据预测样点位置集合中各个预测样点位置的环境信息，以及各个预测样点位置的土壤有机碳密度实测值，采用随机森林方法，训练获得以土壤有机碳密度实测值为目标的第三预测模型；然后分别针对验证样点位置集合中的各个验证样点位置，根据验证样点位置的环境信息，通过第三预测模型，获得该验证样点位置的第二土壤有机碳密度预测值，进而获得验证样点位置集合中各个验证样点位置的第二土壤有机碳密度预测值，构成验证样点位置集合中各个验证样点位置的第二土壤有机碳密度预测值集合PCTM；然后进入步骤013。
[0123] 步骤013.如图4所示，针对目标土壤区域中所有样点位置，统一按预设划分规则，基于L将样点位置所对应的各个归并层划分为各个拟合层，各个样点位置所对应的拟合层、以及拟合层的数量相同，并获得各个拟合层的厚度信息，在本实施例例中，对于1m深度的剖面区域，具体划分为0-5cm、5-15cm、15-30cm、30-60cm与60-100cm的固定深度的各个拟合层，共划分为5个拟合层；然后分别针对预测样点位置集合中的各个预测样点位置，根据预测样点位置各个拟合层所分别对应的各个归并层，针对该预测样点位置各个归并层采样发生层的土壤属性信息进行拟合，获得预测样点位置各拟合层的土壤属性信息；进而获得预测样点位置集合中各个预测样点位置分别所对应各个拟合层的土壤属性信息，并进入步骤014。
[0124] 其中，所述步骤013中，获得预测样点位置集合中各个预测样点位置分别所对应各个拟合层的土壤属性信息，具体包括如下操作：
[0125] 分别针对预测样点位置集合中的各个预测样点位置，根据等面积Spline函数，采用如下公式：
[0126]
[0127] 获得该预测样点位置各个拟合层的土壤属性信息；进而获得预测样点位置集合中各个预测样点位置分别所对应各个拟合层的土壤属性信息，其中，i1＝{1,…,I1}，I1表示预测样点位置集合中预测样点位置的数量， 表示预测样点位置集合中第i1个预测样点位置所对应拟合层的层数， 表示预测样点位置集合中第i1个预测样点位置所对应第k层拟合层的土壤属性信息， 表示 函数在 层与 层拟合结果的平均值， 表示预测样点位置集合中第i1个预测样点位置所对应第k层拟合层的实验室测定分析误差， 函数是等面积Spline函数。
[0128] 步骤014.采用MTC-F方法，根据预测样点位置集合中各个预测样点位置的环境信息，以及各个预测样点位置分别所对应各个拟合层的土壤属性信息，采用随机森林方法，训练获得以各个拟合层土壤属性信息为目标的各个预测模型，构成第四预测模型集合；然后分别针对验证样点位置集合中的各个验证样点位置，根据验证样点位置的环境信息，通过第四预测模型集合中的各个预测模型，获得该验证样点位置分别所对应各个拟合层的第三土壤属性预测信息，进而获得验证样点位置集合中各个验证样点位置分别所对应各个拟合层的第三土壤属性预测信息，然后进入步骤015；其中，第三预测土壤属性信息包括第三有机碳含量预测信息、第三土壤容重预测信息和大于预设直径的第三砾石预测含量。
[0129] 步骤015.根据验证样点位置集合中各个验证样点位置所对应各个拟合层的第三土壤属性预测信息，以及厚度信息，获得验证样点位置集合中各个验证样点位置的第三土壤有机碳密度预测值 构成验证样点位置集合中各个验证样点位置的第三土壤有机碳密度预测值集合PMTC-F；然后进入步骤016。
[0130] 所述步骤015具体包括如下操作：
[0131] 根据验证样点位置集合中各个验证样点位置所对应各个拟合层的第三土壤属性预测信息，以及厚度信息，按如下公式：
[0132]
[0133] 获得验证样点位置集合中各个验证样点位置的第三土壤有机碳密度预测值构成验证样点位置集合中各个验证样点位置的第三土壤有机碳密度预测值集合PMTC-F；其中， 表示验证样点位置集合中第i2个验证样点位置所对应拟合层的总数， 表示验证样点位置集合中第i2个验证样点位置所对应的第k个拟合层，表示验证样点位置集合中第i2个验证样点位置所对应第k个拟合层的第三有机碳含量预测信息， 表示验证样点位置集合中第i2个验证样点位置所对应第k个拟合层的第三土壤容重预测信息， 表示验证样点位置集合中第i2个验证样点位置所对应第k个拟合层中大于预设直径的第三砾石预测含量， 表示验证样点位置集合中第i2个验证样点位置所对应第k个拟合层的厚度信息；然后进入步骤016。
[0134] 步骤016.使用Lins’s一致性相关系数(Lins's concordance correlationcoefficient)作为精度验证指标，进行预测结果的精度评价：
[0135]
[0136] 其中，表示验证样点位置集合中第i2个验证样点位置的土壤有机碳密度实测值，即土壤有机碳密度实测值 表示验证样点位置集合中第i2个验证样点位置的土壤有机碳密度预测值，即分别对应第一土壤有机碳密度预测值、第二土壤有机碳密度预测值、第三土壤有机碳密度预测值；表示验证样点位置集合中所有验证样点位置土壤有机碳密度实测值的平均值；表示验证样点位置集合中所有验证样点位置土壤有机碳密度预测值 的平均值；Lins’s一致性相关系数与精度检验的目标数据量级无关，值域是[-1,1],最大值1表示最好的预测结果，与观测数据完全一致。实际应用过程当中，即分别针对第一土壤有机碳密度预测值集合PMTC-D、第二土壤有机碳密度预测值集合PCTM和第三土壤有机碳密度预测值集合PMTC-F，分别采用上述精度评价方法，其中，采用第一土壤有机碳密度预测值集合PMTC-D时，则将第一土壤有机碳密度预测值集合PMTC-D中的各个值对应于 采用第二土壤有机碳密度预测值集合PCTM时，则将第二土壤有机碳密度预测值集合PCTM中的各个值对应于 采用第三土壤有机碳密度预测值集合PMTC-F时，则将第三土壤有机碳密度预测值集合PMTC-F中的各个值对应于
[0137] 根据土壤有机碳密度实测值集合V进行精度检验，获得第一土壤有机碳密度预测值集合PMTC-D、第二土壤有机碳密度预测值集合PCTM和第三土壤有机碳密度预测值集合PMTC-F中的最优土壤有机碳密度预测值集合，并进入步骤017。
[0138] 步骤017.将目标土壤区域离散化空间栅格数据，将目标土壤区域中所有样点位置分别所对应发生层的采样数据作为预测数据集合，若最优土壤有机碳密度预测值集合为第一土壤有机碳密度预测值集合PMTC-D，则采用步骤009至步骤011的方法获得目标土壤区域的土壤有机碳密度空间分布栅格数据；若最优土壤有机碳密度预测值集合为第二土壤有机碳密度预测值集合PCTM，则采用步骤012的方法获得目标土壤区域的土壤有机碳密度空间分布栅格数据；若最优土壤有机碳密度预测值集合为第三土壤有机碳密度预测值集合PMTC-F，则采用步骤013至步骤015的方法获得目标土壤区域的土壤有机碳密度空间分布栅格数据；然后进入步骤018。
[0139] 步骤018.根据目标土壤区域的土壤有机碳密度空间分布栅格数据，获得目标土壤区域的土壤有机碳储量。
[0140] 将本发明上述设计基于土壤发生层厚度预测的土壤有机碳储量估算方法应用于实际中，如以辽宁省土壤有机碳储量估算为例。
[0141] 辽宁省位于中国东北地区的南部，东北与吉林省接壤，西北与内蒙古自治区为邻。辽宁省土地资源不足，耕地资源更少，土地利用类型较多。以2009年的土壤调查数据为输入，在本发明的基础上，使用不同的预测方法与发生层归并策略，可以获取到基于最优方法的有机碳储量估算值。以辽宁省的实际应用中，最终如图7a至图7c所示，对于不同预测方法PMTC-D、PCTM、PMTC-F的预测结果，使用Lins’s一致性相关系数(Lins'sconcordance correlation coefficient)作为精度验证指标，进行预测结果的精度评价。精度评价结果为CTM最高(ρc＝0.21)，MTC-D其次(ρc＝0.18)，MTC-F最差(ρc＝0.13)。因此，对于辽宁省的有机碳储量估算来说，推荐使用CTM方法预测该区域的有机碳储量；则将目标土壤区域离散化空间栅格数据，将目标土壤区域中所有样点位置分别所对应发生层的采样数据作为预测数据集合，则采用步骤012的方法获得目标土壤区域的土壤有机碳密度空间分布栅格数据；
最后根据目标土壤区域的土壤有机碳密度空间分布栅格数据，获得目标土壤区域的土壤有机碳储量。
[0142] 有别于常规基于GIS的土壤有机碳储量估算方法，本发明充分考虑土壤属性在发生层基础上的均质性与水平维的空间连续性。交叉验证的使用为目标研究区域碳储量估算工程提供了预测精度与最优计算机制保证，有效地客服了一种空间预测技术存在局限性的技术
瓶颈。该技术对于科学制定有效的管理机制、实现资源的可持续发展及充分发挥
生态系统的生态效益具有十分重要的指导意义。
[0143] 上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。
