CN110276160B - 一种平原区土壤有机碳三维空间分布模拟方法 - Google Patents

Abstract

本发明提供一种平原区土壤有机碳三维空间分布模拟方法，包括(1)土壤有机碳剖面变化最优深度函数的拟合及函数参数确定；(2)土壤有机碳剖面深度函数参数的影响因素分析；(3)三维空间分布模拟模型的构建与计算过程。本发明构建地形平缓区土壤性质空间分布模拟的环境辅助变量，利用土壤深度传递函数和回归克里格法，建立地形平缓区土壤有机碳三维空间分布模拟模型，开展该区域土壤有机碳的三维空间分布模拟方法研究，不仅可为准确掌握地形平缓区在高强度农业利用背景下土壤有机碳三维时空演变信息提供可靠的模拟方法，也可为准确估算区域有机碳储量和保障土壤安全及粮食安全提供方法参考。

Description

一种平原区土壤有机碳三维空间分布模拟方法

技术领域

本发明属于土壤环境技术领域，具体涉及一种平原区土壤有机碳三维空间分布模拟方法。

背景技术

土壤有机碳库是陆地生态系统中最大的碳库。据估算，全球1m深土壤有机碳储量大约有1500Pg(1Pg＝10¹⁵g)，是大气碳库的2倍多(750Pg)，是陆地植被碳库的2～3倍。土壤有机碳较小幅度的变化都可能引起大气二氧化碳浓度较大幅度的改变，从而影响全球气候变化。同时，土壤有机碳也是植物矿质营养及有机营养的源泉。土壤中有机碳的储量和动态变化不仅直接影响着全球碳循环，还与土壤肥力、农业生产和环境保护等密切相关。受多种因素综合影响，土壤有机碳具有高度的空间异质性，包括在垂直和水平方向的空间变化。准确获取土壤有机碳三维空间分布信息是区域土壤资源利用、土壤有机碳储量估算及全球碳循环研究的基础，越来越受到国内外土壤学乃至其他地球科学研究领域科研人员的关注。土壤有机碳三维空间分布模拟正在成为土壤性质时空演变研究的重点和热点。

利用相关环境因素作为辅助变量进行土壤有机碳空间分布模拟的方法，能较好地揭示土壤有机碳的空间异质性，进而大幅度提高模拟结果精度。二维水平空间分布模拟中，可直接建立土壤有机碳与环境辅助变量的映射关系。但在三维空间分布模拟中，分层建立不同土层土壤有机碳与环境因子间的对应关系，会割裂土壤有机碳三维空间变化与环境因素关系间的整体性与系统性。直接建立环境因子与土壤剖面有机碳整体变化特征间的映射关系则更符合景观背景下土体的发生发育过程。土壤有机碳深度函数的参数能反映土壤有机碳表层含量及其在剖面上随深度增加的变化速率，这些参数同时也具有明显的空间分布特征，并与各种环境因素密切相关。直接建立环境因素与土壤性质深度函数参数间的映射关系可以有效解决分层二维模拟存在的问题。

平原区地势平坦，人口密度大，农业集约化程度较高。高强度农业利用背景下，平原区土壤有机碳等土壤性质相较其他区域变化更为剧烈，土壤环境问题也更加突出。仅模拟表层土壤有机碳的时空演变，难以全面地掌握高强度农业活动背景下土壤有机碳变化的整体状况。同时，地形等常用环境辅助因子在地形平缓区差别不大，很难用作土壤性质空间分布模拟预测中的辅助变量。

发明内容

针对上述技术问题，本发明提供一种平原区土壤有机碳三维空间分布模拟方法，准确估算区域有机碳储量和保障土壤安全及粮食安全提供方法参考。

具体的技术方案为：

一种平原区土壤有机碳三维空间分布模拟方法，具体包括以下步骤：

1、土壤有机碳剖面变化最优深度函数的拟合及函数参数确定

分别采用对数函数、幂函数和指数函数对分层剖面采样获得的土壤有机碳含量(kg m^-3)的剖面变化进行拟合。利用拟合方程的决定系数(R²)以及误差评价指数(平均绝对误差(MAE)、平均相对误差(MRE)和均方根误差(RMSE))选择土壤有机碳剖面变化的最佳深度拟合函数，进行得到揭示土壤有碳含量剖面变化的最优深度函数类型及土壤有机碳深度函数参数。具体方法如下：

首先将土壤有机碳含量(g kg^-1)转换成有机碳密度(kg m^-3)，公式如下：

SOC_v＝SOC_m×B×(1-θ％) (1)

式中，SOC_v表示土壤有机碳含量(kg m^-3)；SOC_m表示土壤有机碳含量(g kg^-1)；B表示土壤容重(g cm^-3)；θ表示大于2mm的石砾所占百分比(％)。

拟合过程中以有机碳含量(g cm^-3)为纵坐标，以中心土层深度为横坐标，逐个拟合每个剖面土壤有机碳的剖面分布。考虑到农业表层(0～20cm)土壤受人为干扰较大，土壤有机碳含量趋于均一化，将土壤有机碳的剖面分布拟合分为两部分，一共包含土壤有机碳深度函数参数SOC_v1、C₀和k三个参数。

当土层深度≤20cm时：

y＝SOC_v1 (2)

式中，y表示一定深度土壤有机碳含量(kg m^-3)；SOC_v1表示0～20cm范围内土壤有机碳含量(kg m^-3)。

当土层深度≥20cm时，分别采用指数函数、对数函数和幂函数拟合土壤有机碳含量剖面分布，其拟合公式如下：

y＝C₀e^-kx (3)

式中，y表示一定深度土壤有机碳含量(kg m^-3)；x表示中心土层深度(cm)；C₀表示地表土壤有机碳含量(kg m^-3)；k表示土壤有机碳剖面变化的斜率；

y＝-kln(x)+C₀ (4)

式中，y表示一定深度土壤有机碳含量(kg m^-3)；x表示中心土层深度(cm)；k表示土壤有机碳剖面变化的斜率；C₀表示地表土壤有机碳含量(kg m^-3)；

y＝C₀x^-k (5)

式中，y表示一定深度土壤有机碳含量(kg m^-3)；x表示中心土层深度(cm)；C₀表示土壤有机碳含量(kg m^-3)；k表示土壤有机碳剖面变化的斜率。

2、土壤有机碳剖面深度函数参数的影响因素分析

2.1影响因素的选择与处理

根据土壤发生学原理，从成土因素中选取和构建地形平缓区土壤有机碳剖面变化的影响因素，包括成土母质、地形与水文、植被指数、农用地密度、建设用地密度、人口密度和GDP。具体处理方法如下：

(1)成土母质

成土母质是土壤物质形成的基础，对土壤结构、孔隙度、土壤养分含量均具有重要影响，是制约陆地生态系统中植物养分吸收的驱动因素。不同土壤类型因成土过程和耕作管理措施的差异进而对土壤有机碳空间分布产生影响。将收集到的纸质成土母质图进行扫描，在ArcGIS中数字化，获取1:10000的SHP格式的成土母质和土壤类型的空间分布数据，并转化成30m分辨率的栅格数据。

(2)地形与水文

基于1:50000地形图在ArcGIS软件中生成30m分辨率的数字高程模型。利用研究区30m分辨率的数字高程模型进行水文分析，计算研究区汇流面积和地形湿度指数，其中汇流面积值较大，对其取自然对数。考虑到研究区河网密布，河流距离对土壤有机碳空间分布的影响与河流所在的流域有关。本文通过流域分析确定流域分水岭，再在每个流域范围内对河流进行缓冲距离分析，获取与河流距离的空间分布图。

(3)植被指数

采用归一化植被指数(NDVI)研究植被生长对研究区土壤有机碳空间分布的影响。根据土壤样品采集时间，选择与采样时间接近的16天合成的250m分辨率的MODIS数据，获取研究区NDVI信息。根据研究区植被生长情况，选取植被生长旺盛的的2月和7月的MODIS数据，经投影变换、边界裁剪后求取均值，获得研究区NDVI的空间分布信息

(4)人口密度

人口密度的栅格数据来自于中科院资源环境数据平台，空间分辨率为1km。选取与采样时间最为接近的年份的人口密度数据，经过投影变换、边界裁剪后，求取不同时段的均值作为研究区人口密度。将研究区人口密度重采样至30m分辨率。

(5)GDP

GDP的栅格数据来自于资源环境数据平台(徐新良,2017)，空间分辨率为1km。选取与采样时间最为接近的2010年和2015年的GDP数据，经过投影变换、边界裁剪后，求取两个时段的均值作为研究区GDP。将研究区GDP重采样至30m分辨率。

(6)建设用地密度

根据研究区土地利用现状图，结合土壤采样点密度，利用3×3km的空间格网计算单位格网内建设用地比例，将其赋值到3×3km的空间格网上，再转化成栅格数据，重采样至30m分辨率，获得研究区建设用地密度空间分信息。

(7)农用地密度

根据研究区土地利用现状图，结合土壤采样点密度，利用3×3km的空间格网计算单位格网内农用地所占比例，将其赋值到3×3km的空间格网上，再转化成栅格数据，重采样至30m，获得研究区农用地密度空间分信息。

2.2土壤有机碳部分函数参数与影响因素间相关关系分析

利用相关分析和回归分析方法分析土壤有机碳剖面变化深度函数的SOC_v1、C₀和k三个参数与成土母质、地形与水文、植被指数、农用地密度、建设用地密度、人口密度和GDP环境因子间的关系，从而确定建立模型所需的环境辅助因子。

3、三维空间分布模拟模型的构建与计算过程

(1)模型构建

结合所选的环境因子，利用回归克里格法预测土壤有机碳剖面变化深度函数三个参数(SOC_v1、C₀和k)的空间分布，从而实现对土壤有机碳三维空间分布的模拟。模型(RK_E)构建如下:

z(x_i，k，y_j，k)＝f(x_i，k，y_j，k)+r(x_i，k，y_j，k) (6)

式中：Z(x_i,k,y_j,k)表示在点(x_i,y_j)处土壤有机碳深度函数参数(SOC_v1、C₀和k)的实测值，(x_i,y_j)为点位坐标；f(x_i,k,y_j,k)为点位(x_i,y_j)处由该点位成土母质、地形与水文、植被指数、农用地密度、建设用地密度、人口密度和GDP环境因子点位环境因素差异引起的土壤有机碳深度函数参数变异；r(x_i,k,y_j,k)为点位(x_i,y_j)处去除研究中所选因素的影响后剩余的残差项，由研究中未考虑到的其它因素或点位处其它因素引起。

假设2个部分相互独立，可分别由不同的方法获取：

(1)f(x_i,k,y_j,k)：反映点位(x_i,y_j)处土壤有机碳深度函数参数变异，由该点位处成土母质、地形与水文、植被指数、农用地密度、建设用地密度、人口密度和GDP环境因子点位环境因素差异决定，其空间变化采用多元逐步回归获得。

(2)r(x_i,k,y_j,k)：为点位(x_i,y_j)处去除研究中所选因素的影响后剩余的残差项，由样点实测值减去多元回归模型对f(x_i,k,y_j,k)的预测值获得。该部分的空间变化采用普通克里格法来进行模拟。

(2)计算过程

采用该方法模拟土壤有机碳三维空间分布的计算过程如下：

(1)首先利用建模样点，采用多元逐步回归分析方法建立土壤有机碳深度函数参数与环境因子之间的回归模型。然后在ArcGIS中运用所建立的多元回归模型在栅格计算器分别生成SOC_v1、C₀和k三个参数的f(x_i,k,y_j,k)曲面。

(2)以建模样点对应点位SOC_v1、C₀和k的值减去多元回归模型的模拟值，得到各样点去除所选环境因素后的残差项r(x_i,k,y_j,k)；然后在ArcGIS中以普通克里格法对SOC_v1、C₀和k三个参数多元回归模拟残差的空间分布进行模拟，分别生成残差SOC_v1、残差C₀和残差k3个参数的r(x_i,k,y_j,k)曲面。

(3)将多元回归模型模拟得到的SOC_v1、C₀和k三个参数的f(x_i,k,y_j,k)曲面分别与普通克里格法得到的三个参数残差项r(x_i,k,y_j,k)空间分布曲面相加，得到研究区土壤有机碳剖面变化深度函数3个参数的空间分布。

(4)在ArcGIS栅格计算器中，利用最优深度函数类型以及由回归克里格法预测得到的3个参数的空间分布数据，分别输入对应的土壤深度，从而模拟得到研究区0-20cm、20-40cm、40-60cm和60-100cm土层土壤有机碳空间分布。

本发明提供的一种平原区土壤有机碳三维空间分布模拟方法，构建地形平缓区土壤性质空间分布模拟的环境辅助变量，利用土壤深度传递函数和回归克里格法，建立地形平缓区土壤有机碳三维空间分布模拟模型，开展该区域土壤有机碳的三维空间分布模拟方法研究，不仅可为准确掌握地形平缓区在高强度农业利用背景下土壤有机碳三维时空演变信息提供可靠的模拟方法，也可为准确估算区域有机碳储量和保障土壤安全及粮食安全提供方法参考。

附图说明

图1a为实施例建模样点空间分布图；

图1b为实施例验证样点空间分布图；

图2不同预测方法土壤有机碳空间分布结果对比。

具体实施方式

结合实施例说明本发明的具体技术方案。

本实施例选择成都平原核心作为案例区，以上述方案开展该区域土壤有机碳三维空间分布模拟。具体如下：

1数据来源与研究方法

1.1研究区概况

研究区位于成都平原的核心区域，地理坐标介于北纬30°22′～31°50′、东经103°27′～104°14′之间。全区面积约为3170km²，主要包括成都市的彭州、新都、都江堰、郫都、温江、双流、新津、崇州、大邑和邛崃等10个区县的平原区。区内地势总体平缓，海拔在447～732m之间，由西北向东南降低。河网密布，主要有金马河、西河和斜江河三条岷江水系和沱江水系青白江流经境内。该区属亚热带湿润季风气候区，年平均气温在16℃左右，年平均降水量为900～1300mm，全年无霜期为278d。成土母质以岷江冲积物为主，包括灰色冲积物、灰棕冲积物、紫色冲积物以及主要分布于西南部的更新统老冲积物。由于长期的水耕熟化，区内主要土壤类型为水稻土和少量新积土，水稻土一共包含4个亚类。土地利用以耕地为主，包括水稻～小麦和水稻～油菜轮作两种典型种植制度。随着该区域经济发展，农业种植结构调整，形成了以耕地为主，以园地和农林混合用地等多种农业种植制度为辅的农业用地方式。

1.2数据来源与处理

1.2.1土壤有机碳数据来源

根据研究区成土母质图及土地利用现状图，利用3×3km格网法在室内进行采样点布设。在兼顾空间分布均匀性和样点代表性的基础上，参照成都市第二次土壤普查资料中的土壤发生层分类，按0～20、20～40、40～60、60～100cm，于2016年4月和2017年4月在研究区内采集土壤剖面，共采集171个土壤剖面，677个土壤样品。剖面采集过程中，在预先布设好的点位上，根据实际情况选择具有代表性的位置挖土壤剖面。首先去除地表枯落物，再划分土壤发生层并记录其厚度，最后按照0～20、20～40、40～60和60～100cm的间隔在剖面内利用多点混合法自下而上采集土壤样品。本实施例采集的土壤样品包括常规测试分析样品和土壤容重测定样品。常规测试样品每个土样大约采集1.5kg，装袋贴标签。样品采集过程中，利用手持GIS详细记录每个剖面的经纬度和海拔等地形信息，并同时记录每个剖面的成土母质、土壤类型、土地利用方式及灌排条件等环境信息。

土壤样品带回实验室，经自然风干后挑出动植物残体，碾磨过0.149mm筛待测。采用重铬酸钾氧化～外加热测定土壤有机碳含量。测定过程中采用三次重复、空白对照以及国家标准物质样进行数据质量控制。

为评价所建立方法对土壤有机碳三维空间分布模拟的精度，采用独立验证法对不同模拟方法预测结果进行精度评价。利用ArcGIS中的地统计模块从171个土壤剖面中随机抽取80％作为建模点，剩余20％作为验证点，即建模点137个，验证点34个。建模点和验证点位的统计结果、空间分布位置如图1a和图1b所示。

1.2.2环境辅助因素的选取与表达

根据土壤发生学原理，参照前人在该区域的研究结果，本实施例选取地形因子中的汇流面积和地形湿度指数、水文条件(与河流距离)、归一化植被指数(NDVI)和农用地密度、建设用地密度、人口密度和国民生产总值(GDP)作为建模的辅助环境变量。综合考虑研究区图件信息及样点空间采集信息，对所有环境辅助变量重采样到30m分辨率。

(1)地形与水文条件

地形通常是通过影响区域水热的再分配进而对土壤有机碳的空间分布产生影响。平原区内，水文状况通常是通过不同潮滩的淹水时间和水动力条件影响有机物质输入和分解，与河流距离的远近还反映了物质的沉积过程和成土时间的差异，进而影响土壤有机碳分布。研究区地形和水文因子均基于1:50000地形图在ArcGIS软件中生成的30m分辨率数字高程模型，包括汇流面积、地形湿度指数和与河流距离。利用研究区30m分辨率的数字高程模型进行水文分析，计算研究区汇流面积和地形湿度指数，其中汇流面积值较大，对其取自然对数。考虑到研究区河网密布，河流距离对土壤有机碳空间分布的影响与河流所在的流域有关。本文通过流域分析确定小流域分水岭，再在每个流域范围内对河流进行缓冲距离分析，获取与河流距离的空间分布图。

(2)植被指数

采用归一化植被指数(NDVI)研究植被生长对研究区土壤有机碳空间分布的影响。根据土壤样品采集时间和植被生长情况，选择2016年和2017年两年的2月和7月16天合成的250m分辨率MODIS NDVI数据，进行投影变换、边界裁剪、求取均值，进一步得到研究区植被状况。

(3)人口密度

人口密度的栅格数据来自于资源环境数据平台，空间分辨率为1km。选取与采样时间最为接近的2010年和2015年的人口密度数据，经过投影变换、边界裁剪后，求取两个时段的均值作为研究区人口密度。最后将研究区人口密度重采样至30m分辨率。

(4)GDP

GDP的栅格数据来自于资源环境数据平台，空间分辨率为1km。选取与采样时间最为接近的2010年和2015年的GDP数据，经过投影变换、边界裁剪后，求取两个时段的均值作为研究区GDP。最后将研究区GDP重采样至30m分辨率。

(5)建设用地密度

根据研究区2015年土地利用现状图，结合土壤采样点密度，利用3×3km的空间格网计算单位格网内建设用地比例并赋值对应格网上，然后转化成栅格数据，重采样至30m分辨率，从而获得研究区建设用地密度空间分信息。

(6)农用地密度

根据研究区2015年土地利用现状图，结合土壤采样点密度，利用3×3km的空间格网计算单位格网内农用地所占比例并赋值对应的格网上，然后转化成栅格数据，重采样至30m，从而获得研究区农用地密度空间分信息。

1.3空间分布预测方法

(1)普通克里格法

普通克里格法(Ordinary Kriging,OK)是一种最优线性无偏估计插值方法，该方法以区域化变量为基础，应用半变异函数理论，根据采样点实测值对其临近未采样区域内属性值进行无偏最优估计，对地形平缓区土壤性质空间分布模拟具有较好的预测能力。该方法采用半方差分析刻画土壤性质的空间结构，根据半方差最优理论模型和参数在ArcGIS中进行空间插值，被广泛应用于土壤性质的空间分布研究。

(2)回归克里格法

回归克里格法(Regression kriging,RK)是一种结合多元线性回归(Multiplelinear regression,MLR)与普通克里格法(Ordinary kriging,OK)的空间预测方法。与OK方法相比，RK方法引入环境辅助变量，通过刻画土壤性质与环境因子之间的关系，进而预测土壤性质的空间分布，能较好的反应土壤性质与环境因子之间的关系，有助于提高土壤性质空间分布的预测精度。

1.4精度评价方法

采用样点实测值与预测值之间的平均绝对误差(MAE)、平均相对误差(MRE)和均方根误差(RMSE)对研究区土壤有机碳剖面分布拟合结果精度和土壤有机碳三维空间模拟结果精度进行评价。MRE、MAE和RMSE的计算公式如下：

式中，n为验证点的个数；SOC_obs(i)为第i个验证点的实测值；SOC_pred(i)为第i个验证点的模拟值。MAE、MRE和RMSE的值越小，精度越高。

2结果与分析

2.1最优土壤有机碳深度函数类型的确定

分别采用对数函数、幂函数和指数函数对研究区171个剖面水稻土壤有机碳(kgm^-3)的剖面分布进行拟合。三种函数的拟合R²均在0.85以上，能较好拟合研究区土壤有机碳的剖面分布。进一步综合考虑MAE、MRE及RMSE可知，研究区土壤有机碳剖面分布的最优拟合函数为对数函数。因此，本实施例最终选取对数函数作为土壤有机碳剖面变化的拟合方程。

表1不同土壤深度函数拟合精度评价

2.2土壤有机碳深度函数参数统计特征

对拟合研究区土壤有机碳剖面分布的对数函数参数SOC_v1、C₀及k进行统计，结果显示，参数SOC_v1、C₀及k的变异系数分别为20.46％、25.66％和34.37％，均在10％～100％之间，具有中等程度的空间变异性，但参数k的变异程度要高于C₀和SOC_v1。k-s检验结果表明，SOC_v1和k符合正态分布(P>0.05)，参数C₀符合对数正态分布。将参数C₀的原始数据进行自然对数转换，以满足后续数据分析要求。

表2土壤有机碳深度函数参数的描述性统计统计特征

注：*表示自然对数转换后的结果。

Note:*indicates data is transformed with natural logarithm.

2.3土壤有机碳深度函数参数空间结构特征

利用半方差模型刻画土壤有机碳空间变异的结构性和随机性。半方差分析结果显示(表3)，对数函数的3个参数均符合指数模型，模型拟合决定系数R²均在0.90以上，表明指数函数能较好拟合3个参数的空间结构特征。块金系数反映了区域化变量空间变异的随机性和结构性。若块金系数<25％，表明其空间变异主要受结构性因素影响，由气候、成土母质、地形等自然因素引起；若块金系数>25％，表明其空间变异主要受随机性因素影响，由耕作、施肥、农业管理措施等人为活动影响；若块金系数在25％～75％之间，则说明其空间变异由结构性因素和随机性因素共同影响。本实施例中，3个参数的块金系数在38.22％～40.16％之间，表明土壤有机碳的剖面分布由结构性因素和随机性因素共同影响。参数C₀的块金系数>SOC_v1>k，这一分析结果表明，相比于参数k，代表地表土壤有机碳含量的参数C₀和表示表层土壤有机碳含量的参数SOC_v1更倾向于受随机性因素影响。

表3土壤有机碳深度函数参数的半方差模型及参数

2.4影响因素分析

2.4.1土壤有机碳深度函数参数与环境因子间的关系

由相关分析可知(表4)，3个参数与8种定量因素之间均存在着一定的相关关系。其中，3个参数与NDVI和农用地密度表现出极显著的正相关关系，随NDVI和农用地密度的增加而增大；与建设用地密度、人口密度及GDP均表现出极显著的负相关关系(除SOC_v1与建设用地密度呈显著负相关外)，随影响因素的增加而减小；与河流距离表现出显著的正相关关系，随河流距离的增大而增加。尽管3个参数与汇流面积和地形湿度指数之间表现出一定的相关关系，但相关性并不显著。以上结果表明，与河流距离、NDVI、农用地密度、建设用地密度、人口密度及GDP是影响研究区土壤有机碳剖面分布的重要因素。

表4土壤有机碳深度函数参数与环境因子间的相关关系

注：*,**分别表示P<0.05和P<0.01。

Note:*,**indicate significant relationship is tested at the level of0.05and 0.01,respectively.

2.4.2影响程度分析

进一步采用回归分析定量揭示不同影响因素对3个参数的影响程度(表5)，结果显示2种定性因素中，亚类对参数C₀和k均无显著影响(P>0.05)，对SOC_v1有显著影响(P<0.05)，可解释其3.3％的空间变异。成土母质对3个参数均有极显著影响，对3个参数的解释能力均高于亚类，是影响研究区土壤有机碳剖面分布的主控定性因素。6种定量因素对3个参数均有显著影响或极显著影响，对参数SOC_v1空间变异的解释能力在2.3％～20.1％之间，对参数C₀空间变异的解释能力在3.2％～20.7％之间，对参数k空间变异的解释能力在2.1％～16.7％之间。可以看出，人口密度和GDP对参数SOC_v1空间变异的解释能力均为20.1％，对参数C₀空间变异的解释能力分别为20.7％和19.7％，对参数k空间变异的解释能力分别为16.7％和16.2％，均远高于其余4种定量因素，是影响研究区土壤有机碳剖面分布的主要定量因素。

表5与环境因子间的回归分析

2.5三维模拟模型构建与模拟结果

2.5.1建模样点土壤有机碳剖面函数参数与环境因子间的相关性分析

前述结果显示，研究区内成土母质、亚类、与河流距离、NDVI、农用地密度、建设用地密度、人口密度和GPD是影响研究区土壤有机碳三维空间分布的主要因素。定性因素中，成土母质的作用显著高于亚类，可作为影响研究区水稻土有机碳三维空间分布的主控定性因素。进一步分析建模样点与成土母质、与河流距离、NDVI、农用地密度、建设用地密度、人口密度和GPD的关系，结果显示，拟合土壤有机碳剖面分布的参数C₀、k和SOC_v1与各环境因子之间的相关性均要高于不同土层土壤有机碳(SOC_v1、SOC_v2、SOC_v3和SOC_v4)与环境因子之间的相关性，这一结果强调了环境因子对土壤有机碳剖面变异的影响作用。

表6建模样点与环境因子之间的相关性分析

2.5.2回归克里格法对深度函数参数空间分布的预测结果

采用逐步回归分析方法建立土壤有机碳深度函数参数与环境因子之间的回归模型(表7)。3个参数与环境因子之间的回归方程P值均小于0.01，表明回归方程是有效的。从决定矫正系数R²来看，环境因子对3个参数回归模型的解释能力在23.5％～27.7％之间。

表7建模数据集土壤有机碳深度函数参数回归模型

利用半方差函数分析回归模型预测残差的空间分布结构特征。结果表明，3个参数的块金系数在11.76％～21.3％之间，低于25％。这一结果表明，去除所选环境因子的影响后，3个参数仍然具有强烈空间自相关性，其空间变异主要受结构性因素影响。

表8建模数据集的半方差模型及参数

2.5.3融合深度和环境辅助变量的回归克里格法对土壤有机碳三维空间分布的预测结果

在ArcGIS栅格计算器中，利用最优深度函数类型以及由回归克里格法(RK_E)预测得到的3个参数的空间分布数据，分别输入对应的土壤深度，模拟得到研究区0-20cm、20-40cm、40-60cm和60-100cm土层土壤有机碳空间分布，并与分层普通克里格法(OK)、普通克里法结合深度函数的方法(OK_E)的预测结果进行比较。

(1)不同方法预测结果比较

图2为不同方法对土壤有机碳三维空间分布的预测结果。从图中可以看出，3种方法预测的土壤有机碳空间分布趋势基本一致。0～20cm土壤有机碳的高值区主要分布在研究区西南部；20～40cm土壤有机碳呈现出东北和西南高，中间低的空间分布格局；40cm以下有机碳高值区主要分布在研究区的东北部。从不同预测方法来看，普通克里格法(OK和OK_E)预测结果较为平滑，土壤有机碳的高值区和低值区呈明显的斑块状分布，仅能反应土壤有机碳的空间分布趋势。相比于OK，OK_E考虑了土壤有机碳在垂直剖面上的连续变化，弱化了OK方法中极大值和极小值区域的分布范围。相比于普通克里格方法(OK和OK_E)，回归克里格方法(RK_E)融合了定性和定量的环境辅助变量，反应了土壤有机碳空间分布随环境因子变化的关系，其预测结果体现了更多的细节信息，更符合土壤有机碳空间分布的实际情况。

(2)不同方法预测结果精度评价

采用平均绝对误差(MAE)、平均相对误差(MRE)和均方根误差(RMSE)对3种模拟方法的预测结果误差进行评价(表9)。对验证样点误差分析结果表明，结合土壤深度函数的三维空间分布模拟方法(RK_E和OK_E)对土壤有机碳空间分布的预测精度要高于传统的二维分层预测方法(OK)。基于传统二维分层模拟预测结果的平均绝对误差(MAE)为2.39kg m^-3，平均相对误差(MRE)为2.33％，均方根误差(RMSE)为3.08kg m^-3。相比于OK方法，RK_E和OK_E方法的平均绝对误差(MAE)下降了5.02％和和2.51％，平均相对误差(MRE)下降了2.85％和0.72％，均方根误差(RMSE)下降了6.10％和1.95％。结合了定性和定量环境辅助变量的RK_E方法的预测精度要高于仅考虑土壤有机碳空间自相关性的OK_E方法。相比于OK_E方法，RK方法预测结果的平均绝对误差(MAE)下降了2.58％，平均相对误差(MRE)下降了2.13％，均方根误差(RMSE)下降了4.24％。不同预测方法土壤有机碳空间分布结果对比如图2所示。

表9不同预测方法误差比较

3结论

本实施例以成都平原为案例区，以与河流距离、NDVI、农用地密度、建设用地密度、人口密度及GDP作为地形平缓区土壤有机碳空间分布模拟模型的辅助变量，采用土壤深度函数和回归克里格方法(RK_E)相结合，构建了融合定性和定量辅助变量的土壤有机碳三维空间分布预测方法。与传统二维分层模拟方法(OK)，结合土壤深度函数的三维空间预测方法(RK_E和OK_E)对模拟精度更高。与仅考虑土壤有机碳空间自相关性的OK_E方法相比，融合了环境辅助变量的RK_E方法反映了更多的细节信息，模拟结果精度明显提高。结合地形平缓区土壤有机碳的影响因素建立的辅助变量，以及基于土壤深度函数和回归克里格方法构建的土壤有机碳三维空间分布预测方法(RK_E)，为准确获得平原区土壤有机碳三维空间分布特征提供了方法参考。

Claims (3)

1.一种平原区土壤有机碳三维空间分布模拟方法，其特征在于，包括以下步骤：

(1)土壤有机碳剖面变化最优深度函数的拟合及函数参数确定

分别采用对数函数、幂函数和指数函数对分层剖面采样获得的土壤有机碳含量的剖面变化进行拟合；利用拟合方程的决定系数R²以及误差评价指数即平均绝对误差MAE、平均相对误差MRE和均方根误差RMSE，选择土壤有机碳剖面变化的最佳深度拟合函数，进行得到揭示土壤有碳含量剖面变化的最优深度函数类型及其土壤有机碳剖面变化深度函数参数即SOC_v1、C₀和k三个参数；

步骤(1)具体方法为：

首先将土壤有机碳含量转换成有机碳密度，公式如下：

SOC_v＝SOC_m×B×(1-θ％) (1)

式中，SOC_v表示土壤有机碳含量，kg m^-3；SOC_m表示土壤有机碳含量，gkg^-1；B表示土壤容重，g cm^-3；θ表示大于2mm的石砾所占百分比，％；

拟合过程中以单位为g cm^-3的土壤有机碳含量为纵坐标，以中心土层深度为横坐标，逐个拟合每个剖面土壤有机碳的剖面分布；

考虑到农业表层即0～20cm土壤受人为干扰较大，土壤有机碳含量趋于均一化，将土壤有机碳的剖面分布拟合分为两部分，一共包含SOC_v1、C₀和k三个参数；

当土层深度≤20cm时：

y＝SOC_v1 (2)

式中，y表示一定深度土壤有机碳含量，kg m^-3；SOC_v1表示0～20cm范围内土壤有机碳含量，kg m^-3；

当土层深度≥20cm时，分别采用指数函数、对数函数和幂函数拟合土壤有机碳含量剖面分布，其拟合公式如下：

y＝C₀e^-kx (3)

式中，y表示一定深度土壤有机碳含量，kg m^-3；x表示中心土层深度，cm；C₀表示地表土壤有机碳含量，kg m^-3；k表示土壤有机碳剖面变化的斜率；

y＝-kln(x)+C₀ (4)

式中，y表示一定深度土壤有机碳含量，kg m^-3；x表示中心土层深度，cm；k表示土壤有机碳剖面变化的斜率；C₀表示地表土壤有机碳含量，kg m^-3；

y＝C₀x^-k (5)

式中，y表示一定深度土壤有机碳含量，kg m^-3；x表示中心土层深度，cm；C₀表示土壤有机碳含量，kg m^-3；k表示土壤有机碳剖面变化的斜率；

(2)土壤有机碳剖面深度函数参数的影响因素分析

(2.1)影响因素的选择与处理

根据土壤发生学原理，从成土因素中选取和构建地形平缓区土壤有机碳剖面变化的环境因子，包括成土母质、地形与水文、植被指数、农用地密度、建设用地密度、人口密度和GDP；

(2.2)土壤有机碳部分函数参数与影响因素间相关关系分析

利用相关分析和回归分析方法分析土壤有机碳剖面变化深度函数的SOCv1、C0和k三个参数与成土母质、地形与水文、植被指数、农用地密度、建设用地密度、人口密度和GDP环境因子间的关系，从而确定建立模型所需的环境辅助因子；

(3)三维空间分布模拟模型的构建与计算过程

(3.1)模型构建

结合所选的环境因子，利用回归克里格法预测土壤有机碳剖面变化深度函数参数即SOC_v1、C₀和k三个参数的空间分布，从而实现对土壤有机碳三维空间分布的模拟；模型RK_E构建如下：

z(x_i，k，y_j，k)＝f(x_i，k，y_j，k)+r(x_i，k，y_j，k) (6)

式中：Z(x_i，k，y_j，k)表示在点(x_i，y_j)处土壤有机碳深度函数参数的实测值，(x_i，y_j)为点位坐标；f(x_i，k，y_j，k)为点位(x_i，y_j)处由该点位环境因子点位环境因素差异引起的土壤有机碳深度函数参数变异；r(x_i，k，y_j，k)为点位(x_i，y_j)处去除研究中所选因素的影响后剩余的残差项，由研究中未考虑到的其它因素或点位处其它因素引起；

(3.2)计算过程

采用该方法模拟土壤有机碳三维空间分布的计算过程如下：

(1)首先利用建模样点，采用多元逐步回归分析方法建立土壤有机碳深度函数参数与环境因子之间的回归模型；然后在ArcGIS中运用所建立的多元回归模型在栅格计算器分别生成SOC_v1、C₀和k三个参数的f(x_i，k，y_j，k)曲面；

(2)以建模样点对应点位SOC_v1、C₀和k的值减去多元回归模型的模拟值，得到各样点去除所选环境因素后的残差项r(x_i，k，y_j，k)；然后在ArcGIS中以普通克里格法对SOC_v1、C₀和k三个参数多元回归模拟残差的空间分布进行模拟，分别生成残差SOC_v1、残差C₀和残差k 3个参数的r(x_i，k，y_j，k)曲面；

(3)将多元回归模型模拟得到的SOC_v1、C₀和k三个参数的f(x_i，k，y_j，k)曲面分别与普通克里格法得到的三个参数残差项r(x_i,k,y_j,k)空间分布曲面相加，得到研究区土壤有机碳深度函数参数的空间分布；

(4)在ArcGIS栅格计算器中，利用最优深度函数类型以及由回归克里格法预测得到的土壤有机碳深度函数参数的空间分布数据，分别输入对应的土壤深度，从而模拟得到研究区0-20cm、20-40cm、40-60cm和60-100cm土层土壤有机碳空间分布。

2.根据权利要求1所述的一种平原区土壤有机碳三维空间分布模拟方法，其特征在于，步骤(2.1)具体方法为：

具体处理方法如下：

(1)成土母质

获取1:10000的SHP格式的成土母质和土壤类型的空间分布数据，并转化成30m分辨率的栅格数据；

(2)地形与水文

基于1:50000地形图在ArcGIS软件中生成30m分辨率的数字高程模型；利用研究区30m分辨率的数字高程模型进行水文分析，计算研究区汇流面积和地形湿度指数，其中汇流面积值较大，对其取自然对数；考虑到研究区河网密布，河流距离对土壤有机碳空间分布的影响与河流所在的流域有关；通过流域分析确定流域分水岭，再在每个流域范围内对河流进行缓冲距离分析，获取与河流距离的空间分布图；

(3)植被指数

采用归一化植被指数NDVI研究植被生长对研究区土壤有机碳空间分布的影响；根据土壤样品采集时间，选择与采样时间接近的16天合成的250m分辨率的MODIS数据，获取研究区NDVI信息；根据研究区植被生长情况，选取植被生长旺盛的的2月和7月的MODIS数据，经投影变换、边界裁剪后求取均值，获得研究区NDVI的空间分布信息；

(4)人口密度

人口密度的栅格数据来自于中科院资源环境数据平台，空间分辨率为1km；选取与采样时间最为接近的年份的人口密度数据，经过投影变换、边界裁剪后，求取不同时段的均值作为研究区人口密度；将研究区人口密度重采样至30m分辨率；

(5)GDP

GDP的栅格数据经过投影变换、边界裁剪后，求取两个时段的均值作为研究区GDP；将研究区GDP重采样至30m分辨率；

(6)建设用地密度

根据研究区土地利用现状图，结合土壤采样点密度，利用3×3km的空间格网计算单位格网内建设用地比例，将其赋值到3×3km的空间格网上，再转化成栅格数据，重采样至30m分辨率，获得研究区建设用地密度空间分信息；

(7)农用地密度

根据研究区土地利用现状图，结合土壤采样点密度，利用3×3km的空间格网计算单位格网内农用地所占比例，将其赋值到3×3km的空间格网上，再转化成栅格数据，重采样至30m，获得研究区农用地密度空间分信息。

3.根据权利要求1所述的一种平原区土壤有机碳三维空间分布模拟方法，其特征在于，步骤(3.1)模型构建中，f(x_i,k,y_j,k)和)r(x_i,k,y_j,k)的获得方法为：

假设2个部分相互独立，可分别由不同的方法获取：

(1)f(x_i,k,y_j,k)：反映点位(x_i,y_j)处土壤有机碳深度函数参数变异，由该点位处环境因子点位环境因素差异决定，其空间变化采用多元逐步回归获得；

(2)r(x_i,k,y_j,k)：为点位(x_i,y_j)处去除研究中所选因素的影响后剩余的残差项，由样点实测值减去多元回归模型对f(x_i,k,y_j,k)的预测值获得；该部分的空间变化采用普通克里格法来进行模拟。

Images