一种黑土可蚀性的反演模型构建方法

1.本发明涉及一种黑土可蚀性的反演模型构建方法，属于土壤侵蚀技术领域。


背景技术：

2.水土流失对耕地质量建设与保护影响严重，如何科学、准确地评估耕地土壤侵蚀风险一直是耕地保护的研究热点。经典的通用土壤流失方程(usle)和其修订版(rusle)以及中国土壤流失方程(csle)均考虑了降雨侵蚀力(r)、土壤可蚀性(k)、地形(ls)、植被覆盖与管理(c)和水土保持措施(p)，为土壤侵蚀评估提供了简化的综合框架，是目前最适合在区域尺度上快速评估水蚀的模型。
3.土壤可蚀性被国内外学者公认为土壤侵蚀预报模型的核心因子，土壤可蚀性反演主要采用实地调查数据与经验模型计算结合的方法，主要依赖于1979年的第二次土壤普查历史资料，依据该资料难以开展大范围、高密度的土壤调查取样，原因是其不仅缺乏对土壤理化性质时空动态性、地表环境和资源分布情况的考虑，且忽视了近年来黑土地受高强度垦殖活动可能产生的基础理化属性重分布现象，测算结果难以真实有效反映当前漫川漫岗区核心产区的土壤可蚀性现状，评估结果的可信度被大大降低。


技术实现要素：

4.本发明针对以上问题，提出一种黑土可蚀性的反演模型构建方法，该方法能够快速且稳健地获取空间连续的土壤可蚀性数据，为科学评估大面积长时间序列条件下的土壤侵蚀强度奠定精度和效率基础。
5.本发明的技术方案为：一种黑土可蚀性的反演模型构建方法，包括以下步骤：
6.(1)野外黑土土壤数据采集；
7.(2)对采集到的黑土土壤数据进行处理，获取其理化数据和光谱数据；
8.(3)获取多种环境因子，并基于多时相卫星遥感影像识别稳定耕地；
9.(4)基于土壤理化性质和多因子算法构建黑土土壤可蚀性评估模型。
10.进一步地，步骤(1)的具体过程如下：
11.(1.1)依据soilgrids中土壤有机碳以及土壤质地性质的值域，设置多个土壤采样点，并使用手持gps记录每个采样点的经纬度坐标；
12.(1.2)记录农业管理方式，包括记录种植结构、轮作方式、农田灌溉条件、秸秆覆盖程度信息；
13.(1.3)在各采样点进行地表土样的采集，并将每个采样点采集到的土样一分为二，一份保持原状，用于后期样本留档和测试对照，另一份再分为多份用于理化分析和光谱测试。
14.进一步地，步骤(2)土样的理化数据和光谱数据获取如下：
15.(2.1)将步骤(1.3)得到的于各采样点采集到的用于理化分析和光谱测试的其中一份土样分别进行土壤有机碳含量测定：具体为使用variomaxcn分析仪通过干烧原理得到
总碳含量，然后通过使用压力钙计法测量出10％hcl处理下出现明显反应样品的无机碳含量，最后从总碳中减去无机碳含量即可得到土壤有机碳含量；
16.(2.2)对于各采样点采集到的用于理化分析和光谱测试的其他份土样分别进行土壤ph、土壤质地、容重和含水量的测定；
17.(2.3)采用湿筛法对土样的土壤团聚体稳定性中的平均质量直径指标进行测定；
18.(2.4)采用fieldspec4地物光谱仪和vertex70傅里叶变换红外光谱仪对土样在350～2500nm波段范围进行数据采集，重复扫描30次，并取其平均值作为输出数据；
19.(2.5)删除步骤(2.4)所得的波段范围两端的低信噪比数据，仅保留400～2450nm波段范围，采用二阶多项式拟合和11个平滑点对光谱进行savitzky-golay平滑处理，采用主成分分析法结合马氏距离剔除离散数据。
20.进一步地，步骤(3)具体为：
21.(3.1)收集包括气候、地形、植被和土壤相关的研究区环境因子空间数据；
22.(3.2)对于通过多时相卫星遥感影像识别得到的稳定耕地图像使用snap软件哨兵二号l1c多光谱影像完成大气校正，生成boa影像；然后使用sen2cor插件进一步去除识别为薄云、云阴影、深色羽毛阴影和薄卷云像素坏点；
23.(3.3)基于联合国粮农组织作物日历确定裸土时间窗口，构建单时相遥感影像裸土像元数据集，利用几何平均值、中位数方法，通过多线性函数合成5年期多时相裸土像元复合数据；
24.(3.4)利用深度可分离卷积u-net算法分离耕地，对照上述单时相裸土像元数据和第三次国土调查中的不稳定耕地，剔除不稳定耕地图斑。
25.进一步地，步骤(4)的具体过程如下：
26.(4.1)基于传统指标和土壤可蚀性综合指数计算土壤可蚀性，包括：
27.平均重量直径:
28.几何平均直径:
29.土壤结构稳定指数：
30.sco胶结剂指数：
31.粘粒比：
32.k因子：k因子：
33.土壤可蚀性综合指数:
34.式中，xi为第i粒级团聚体粒度的平均粒径，yi为第i粒级团聚体粒度重量；som＝c
×
1.724，sn1＝1-sand/100，clay是土壤黏粒含量，土壤黏粒的直径小于0.002mm；silt是土
壤粉粒含量，土壤粉粒的直径为0.002～0.05mm；c是土壤有机碳含量；sand是土壤砂粒含量，土壤砂粒的直径是0.05～2.0mm；n为指标个数，wei和sei分别为第i个土壤可蚀性指数的权重和得分；
35.(4.2)基于主成分分析法得到土壤可蚀性综合指数的权重值，scai、cr、k的得分通过如下公式取得：
[0036][0037]
mwd、gmd和sssi的得分通过如下公式取得：
[0038][0039]
(4.3)利用snap软件提取采样点的裸土像元数据，选取b2-b8、b8a、b11和b12波段作为反演光谱模型的自变量，选取偏最小二乘法回归、随机森林、支持向量机、人工神经网络方法构建土壤属性反演模型；
[0040]
(4.4)采用独立的验证数据集和10倍交叉验证测试模型性能，利用决定系数r2、均方根误差rmse、相对分析误差rpd对预测结果进行评价和验证，选取评估土壤可蚀性的最优反演模型，公式分别如下：
[0041][0042][0043][0044]
式中，yi和分别为检验样本的实测值和预测值，为样本实测值的平均值，n为预测样本数；sd为实测值标准偏差。
[0045]
进一步地，所述的地表土样的采集在春耕前进行，且地表土样为距离表层0～20cm的土壤样品。
[0046]
进一步地，所述的土壤ph、土壤质地、容重和含水量的测定分别使用使用phs-3e仪、ls13320激光衍射仪、环刀法和恒温烘箱进行操作。
[0047]
本发明利用有限的地面采样和卫星影像数据，进行大范围多因子的黑土可蚀性反演，解决当前难以开展大范围、高密度土壤调查取样的限制，有效提升了黑土土壤可蚀性反演的时空分辨率和准确度，为科学评估大面积长时间序列条件下的土壤侵蚀强度奠定了精度和效率基础。
附图说明
[0048]
图1是本发明的工作流程图；
[0049]
图2是本发明实施例研究区域多时相裸土合成稳定耕地技术示意图；
[0050]
图3是本发明实施例研究区域多时相裸土合成稳定耕地效果图；
[0051]
图4是本发明实施例的土壤平均重量直径反演模型验证结果图；
[0052]
图5是本发明实施例的土壤粒径分级反演模型验证结果图；
[0053]
图6是本发明实施例的土壤有机碳反演模型验证结果图；
[0054]
图7是本发明的土壤可蚀性因子空间分布示意图。
具体实施方式
[0055]
为能清楚说明本专利的技术特点，下面通过具体实施方式，并结合其附图，对本专利进行详细阐述。
[0056]
本实施例的研究区位于吉林省白城市，如图1所示，黑土可蚀性的反演模型构建方法的步骤如下：
[0057]
(1)根据白城市行政区划大小，在arcgis等gis工具中基于公里格网选取14个近似均匀分布的地点，进行野外土壤数据采集，具体步骤为：
[0058]
(1.1)在野外进行地表土样本采集，关键过程包括依据soilgrids中土壤有机碳以及土壤质地性质的值域，在春耕前进行表层土壤(0～20cm)样品的采集，并将采集到的各采样点的土样一分为二，一份保持原状，用于后期样本留档和测试对照，另一份再分为多份用于理化分析和光谱测试；
[0059]
(1.2)记录采样点坐标，使用手持gps记录每个采样点的经纬度坐标，所选用的gps设备没有特殊要求，采用本领域技术人员熟知的操作进行即可；
[0060]
(1.3)记录农业管理方式，包括记录种植结构、轮作方式、农田灌溉条件、秸秆覆盖程度信息；
[0061]
(2)对各采集点采集到的土样分别进行内业数据处理得到土壤理化数据，包括测定总碳含量、土壤ph、土壤质地、容重、含水量、平均质量直径，并获取所采集土壤400～2450nm波段的光谱数据；具体的步骤包括：
[0062]
(2.1)将土样的总碳含量测定使用variomaxcn分析仪通过干烧原理进行，使用压力钙计法测量10％hcl处理下出现明显反应样品的无机碳含量，然后从总碳中减去无机碳得到土壤有机碳含量；
[0063]
(2.2)分别使用使用phs-3e仪、ls13320激光衍射仪、环刀法和恒温烘箱对上述土样进行土壤ph、土壤质地、容重、含水量的测定；
[0064]
(2.3)采用湿筛法对土壤团聚体稳定性中的平均质量直径(mwd)进行测定；
[0065]
(2.4)采用fieldspec4地物光谱仪和vertex70傅里叶变换红外光谱仪对供试土壤在350～2500nm波段范围进行数据采集；每个采样点的土样重复扫描30次，并取其平均值作为输出数据；
[0066]
(2.5)删除波段范围两端的低信噪比数据，仅保留400～2450nm波段范围，采用二阶多项式拟合和11个平滑点对光谱进行savitzky-golay平滑处理，采用主成分分析法结合马氏距离剔除离散数据。
[0067]
(3)获取包括降水、地温、地形、植被、土壤研究区环境因子数据，并基于哨兵2号l1c多时相卫星遥感影像、联合国粮农组织作物日历、深度可分离卷积u-net算法获取稳定耕地数据；如图2和图3所示，具体步骤包括：
[0068]
(3.1)收集包括气候(trmm降水、modis地温)、地形(alosdem)、植被(哨兵2ndvi、modisnpp)、土壤相关的研究区环境因子空间数据；
[0069]
(3.2)对于通过多时相卫星遥感影像识别得到的稳定耕地图像使用snap软件哨兵二号l1c多光谱影像完成大气校正，生成boa影像；然后使用sen2cor插件进一步去除识别为薄云、云阴影、深色羽毛阴影和薄卷云像素坏点；
[0070]
(3.3)基于联合国粮农组织作物日历确定裸土时间窗口，构建单时相遥感影像裸土像元数据集，利用几何平均值、中位数方法，通过多线性函数合成5年期多时相裸土像元复合数据；
[0071]
(3.4)利用深度可分离卷积u-net算法分离耕地，对照上述单时相裸土像元数据和第三次国土调查中的不稳定耕地，剔除不稳定耕地图斑。
[0072]
(4)基于土壤理化数据和多因子算法构建吉林省白城市土壤可蚀性评估模型，并得到白城市土壤可蚀性分布图；具体步骤包括：
[0073]
(4.1)基于6种传统方式和土壤可蚀性综合指数csei计算土壤可蚀性，包括：平均重量直径:
[0074]
几何平均直径:
[0075]
土壤结构稳定指数：
[0076]
sco胶结剂指数：
[0077]
粘粒比：
[0078]
k因子：k因子：
[0079]
土壤可蚀性综合指数:
[0080]
式中，xi为第i粒级团聚体粒度的平均粒径，yi为第i粒级团聚体粒度重量；som＝c
×
1.724，sn1＝1-sand/100，clay是土壤黏粒含量，土壤黏粒的直径小于0.002mm；silt是土壤粉粒含量，土壤粉粒的直径为0.02～0.05mm；c是土壤有机碳含量；sand是土壤砂粒含量，土壤砂粒的直径是0.05～2.0mm；n为指标个数，wei和sei分别为第i个土壤可蚀性指数的权重和得分；
[0081]
(4.2)基于主成分分析法得到土壤可蚀性综合指数的权重值，scai、cr、k的得分通过如下公式取得：
[0082][0083]
mwd、gmd和sssi的得分通过如下公式取得：
[0084][0085]
(4.3)利用snap软件提取采样点的裸土像元数据，选取b2-b8、b8a、b11和b12波段作为反演光谱模型的自变量，选取偏最小二乘法回归、随机森林、支持向量机、人工神经网络方法构建土壤属性反演模型；
[0086]
(4.4)采用独立的验证数据集和10倍交叉验证测试模型性能，利用利用决定系数r2、均方根误差rmse、相对分析误差rpd对预测结果进行评价和验证，选取评估土壤可蚀性的最优反演模型，公式分别如下：
[0087][0088][0089][0090]
式中，yi和分别为检验样本的实测值和预测值，为样本实测值的平均值，n为预测样本数；sd为实测值标准偏差。
[0091]
步骤(4)中的指标反演和验证结果如图4、图5和图6所示，其反演结果具有较好的统计学意义；其中，图4为土壤平均重量直径反演模型验证结果图，图中的ptf方法为通过多元线性回归方法从测量土壤性质中导出土壤传递函数来预测指标，spf方法为通过建立光谱传递函数，将土壤可见-近红外光谱与预测指标直接联系起来，采用偏最小二乘回归(plsr)建立spf模型，spf+ptf方法为将ptf方法和spf方法结合起来，即首先利用光谱数据预测土壤基本性质，然后根据预测的土壤性质建立土壤传递函数；
[0092]
图5为土壤粒径分级反演模型验证结果图，分别为粒径(predicted)小于63μm、粒径处于63～250μm之间以及粒径大于250μm时，分别使用ptf方法、spf方法以及spf+ptf方法的反演模型验证结果；
[0093]
图6为土壤有机碳反演模型验证结果图；
[0094]
步骤(4.4)最终所得模型即为利用有限地面采样和卫星遥感影像的大范围多因子黑土可蚀性反演模型，基于所得反演模型，通过常见的栅格计算工具如arcgis栅格计算器，计算并输出白城市区域的土壤可蚀性分布图，如图7所示。

技术特征：
1.一种黑土可蚀性的反演模型构建方法，其特征在于，包括以下步骤：(1)野外黑土土壤数据采集；(2)对采集到的黑土土壤数据进行处理，获取其理化数据和光谱数据；(3)获取多种环境因子，并基于多时相卫星遥感影像识别稳定耕地；(4)基于土壤理化性质和多因子算法构建黑土土壤可蚀性评估模型。2.根据权利要求1所述的黑土可蚀性的反演模型构建方法，其特征在于，步骤(1)的具体过程如下：(1.1)依据soilgrids中土壤有机碳以及土壤质地性质的值域，设置多个土壤采样点，并使用手持gps记录每个采样点的经纬度坐标；(1.2)记录农业管理方式，包括记录种植结构、轮作方式、农田灌溉条件、秸秆覆盖程度信息；(1.3)在各采样点进行地表土样的采集，并将每个采样点采集到的土样一分为二，一份保持原状，用于后期样本留档和测试对照，另一份再分为多份用于理化分析和光谱测试。3.根据权利要求1或2所述的黑土可蚀性的反演模型构建方法，其特征在于，步骤(2)土样的理化数据和光谱数据获取如下：(2.1)将步骤(1.3)得到的于各采样点采集到的用于理化分析和光谱测试的其中一份土样分别进行土壤有机碳含量测定：具体为使用variomaxcn分析仪通过干烧原理得到总碳含量，然后通过使用压力钙计法测量出10％hcl处理下出现明显反应样品的无机碳含量，最后从总碳中减去无机碳含量即可得到土壤有机碳含量；(2.2)对于各采样点采集到的用于理化分析和光谱测试的其他份土样分别进行土壤ph、土壤质地、容重和含水量的测定；(2.3)采用湿筛法对土样的土壤团聚体稳定性中的平均质量直径指标进行测定；(2.4)采用fieldspec4地物光谱仪和vertex70傅里叶变换红外光谱仪对土样在350～2500nm波段范围进行数据采集，重复扫描30次，并取其平均值作为输出数据；(2.5)删除步骤(2.4)所得的波段范围两端的低信噪比数据，仅保留400～2450nm波段范围，采用二阶多项式拟合和11个平滑点对光谱进行savitzky-golay平滑处理，采用主成分分析法结合马氏距离剔除离散数据。4.根据权利要求3所述的黑土可蚀性的反演模型构建方法，其特征在于，步骤(3)具体为：(3.1)收集包括气候、地形、植被和土壤相关的研究区环境因子空间数据；(3.2)对于通过多时相卫星遥感影像识别得到的稳定耕地图像使用snap软件哨兵二号l1c多光谱影像完成大气校正，生成boa影像；然后使用sen2cor插件进一步去除识别为薄云、云阴影、深色羽毛阴影和薄卷云像素坏点；(3.3)基于联合国粮农组织作物日历确定裸土时间窗口，构建单时相遥感影像裸土像元数据集，利用几何平均值、中位数方法，通过多线性函数合成5年期多时相裸土像元复合数据；(3.4)利用深度可分离卷积u-net算法分离耕地，对照上述单时相裸土像元数据和第三次国土调查中的不稳定耕地，剔除不稳定耕地图斑。5.根据权利要求4所述的黑土可蚀性的反演模型构建方法，其特征在于，步骤(4)的具
体过程如下：(4.1)基于传统指标和土壤可蚀性综合指数计算土壤可蚀性，包括：平均重量直径:几何平均直径:土壤结构稳定指数：sco胶结剂指数：粘粒比：k因子：k因子：土壤可蚀性综合指数:式中，xi为第i粒级团聚体粒度的平均粒径，yi为第i粒级团聚体粒度重量；som＝c
×
1.724，sn1＝1-sand/100，clay是土壤黏粒含量，土壤黏粒的直径小于0.002mm；silt是土壤粉粒含量，土壤粉粒的直径为0.002～0.05mm；c是土壤有机碳含量；sand是土壤砂粒含量，土壤砂粒的直径是0.05～2.0mm；n为指标个数，wei和sei分别为第i个土壤可蚀性指数的权重和得分；(4.2)基于主成分分析法得到土壤可蚀性综合指数的权重值，scai、cr、k的得分通过如下公式取得：mwd、gmd和sssi的得分通过如下公式取得：(4.3)利用snap软件提取采样点的裸土像元数据，选取b2-b8、b8a、b11和b12波段作为反演光谱模型的自变量，选取偏最小二乘法回归、随机森林、支持向量机、人工神经网络方法构建土壤属性反演模型；(4.4)采用独立的验证数据集和10倍交叉验证测试模型性能，利用决定系数r2、均方根误差rmse、相对分析误差rpd对预测结果进行评价和验证，选取评估土壤可蚀性的最优反演模型，公式分别如下：
式中，y
i
和分别为检验样本的实测值和预测值，为样本实测值的平均值，n为预测样本数；sd为实测值标准偏差。6.根据权利要求2所述的黑土可蚀性的反演模型构建方法，其特征在于，所述的地表土样的采集在春耕前进行，且地表土样为距离表层0～20cm的土壤样品。7.根据权利要求3所述的黑土可蚀性的反演模型构建方法，其特征在于，所述的土壤ph、土壤质地、容重和含水量的测定分别使用使用phs-3e仪、ls13320激光衍射仪、环刀法和恒温烘箱进行操作。

技术总结
本发明公开了一种黑土可蚀性的反演模型构建方法，步骤包括：野外黑土土壤数据采集；对采集到的黑土土壤数据进行处理，获取其理化数据和光谱数据；获取多种环境因子，并基于多时相卫星遥感影像识别稳定耕地；基于土壤理化性质和多因子算法构建黑土土壤可蚀性评估模型。本发明利用有限的地面采样和卫星影像数据，进行大范围多因子的黑土可蚀性反演，解决了当前难以开展大范围、高密度土壤调查取样的限制，有效提升了黑土土壤可蚀性反演的时空分辨率和准确度，为科学评估大面积长时间序列条件下的土壤侵蚀强度奠定了精度和效率基础。的土壤侵蚀强度奠定了精度和效率基础。的土壤侵蚀强度奠定了精度和效率基础。


技术研发人员：祝元丽 桑学佳 吴子豪 冯向阳
受保护的技术使用者：成都信息工程大学
技术研发日：2023.05.17
技术公布日：2023/8/9