黄土高原南北样带不同土层土壤容重变异分析与模拟

易小波 邵明安,2 赵春雷 张晨成,2

(1.西北农林科技大学资源环境学院， 陕西杨凌 712100；2.西北农林科技大学黄土高原土壤侵蚀与旱地农业国家重点实验室， 陕西杨凌 712100)

黄土高原南北样带不同土层土壤容重变异分析与模拟

易小波1邵明安1,2赵春雷1张晨成1,2

(1.西北农林科技大学资源环境学院， 陕西杨凌 712100；2.西北农林科技大学黄土高原土壤侵蚀与旱地农业国家重点实验室， 陕西杨凌 712100)

为探明黄土高原南北样带土壤容重空间分布特征，为土壤水文过程模拟与预测提供水力参数，采用经典统计学方法，分析了样带不同土层深度(0～10 cm、10～20 cm、20～40 cm)土壤容重的空间变异特征，并用多元逐步回归、传递函数方程和一阶自回归状态空间模型方法分别对土壤容重的空间分布进行了模拟。结果表明：样带0～20 cm深度土壤容重的变异为中等程度变异，20～40 cm为弱变异。状态空间方程转换系数表明，不同土层深度土壤容重的影响因素不同，0～10 cm主要为有机碳含量、黏粒和砂粒体积分数，10～20 cm为有机碳含量、黏粒和砂粒体积分数和降水量，20～40 cm为黏粒和砂粒体积分数、降水量和土地利用。状态空间模型的模拟效果均优于经典统计的多元逐步回归方程和传递函数方程，基于黏粒和砂粒体积分数、降水量和土地利用的状态空间模型可以解释样带20～40 cm容重92.3%的变异。一阶自回归状态空间模型可用于田间条件下土壤容重分布特征的预测。

容重； 状态空间； 传递函数； 模拟； 空间变异

引言

容重是土壤重要的基本物理性质，对土壤的透气性、入渗性、持水性、溶质迁移以及土壤抗蚀能力等特性都有重要影响，是衡量土壤环境好坏的重要指标之一[1-2]。土壤容重通过改变土壤硬度影响植物的生长，随着容重的提高，土壤硬度随之增大，土壤中根系伸长速度减慢且根变短变粗[3]。容重对于养分在土壤中的扩散也具有十分显著的影响，由于土壤固体颗粒的排列不同，造成了不同的曲折通路，单位土体内的电荷密度随容重的改变而改变，进而影响扩散系数[4]。此外，土壤容重也是准确估算土壤碳氮储量所必需的重要参数[1]。

土壤容重的传统测量方法——环刀法，操作耗时，且在野外大尺度实际应用中存在工作繁重，人力、物力耗费较大等问题[5]。近年来，土壤容重测量的新方法也不断涌现。研究人员利用土壤质量含水率、体积含水率和容重的关系，分别采用MP水分分析仪和TDR反射水分仪，同时结合土钻法间接测量土壤容重。但是，受限于土壤体积含水率获取的局限性，这种间接测量方法仍不能实现野外大面积的实时连续测量[6]。土壤容重预测模型作为简便、快捷的土壤容重获取方法，越来越多地受到国内外学者的关注和应用[7-10]。已有的土壤容重预测模型中，所需要的输入参数过多，而且大部分模型模拟精度相对较低，严重影响了模型的实用性[8-9]。WANG等[7]利用传递函数方程对黄土高原区域尺度表层土壤容重进行了模拟，并基于黏粒和粉粒体积分数、有机碳含量、坡度及其转换组合建立了传递函数预测模型。然而，其研究并未考虑土层深度对土壤容重模拟的影响。

由于成土母质、过程、气候、生物及耕作的影响，自然条件下的土壤容重具有较高的空间变异性[11]。传统的空间变异理论未考虑土壤参数与地理位置间的空间关系，对空间变异的解释不够充分，状态空间方程对土壤及植被的模拟效果要明显优于线性回归和多元回归等传统方法[12]。贾小旭等[13]对黄土高原北部典型坡面草地表层土壤水分进行了状态空间模拟，结果表明状态空间方程可用于环境因素复杂的黄土高原水蚀风蚀交错区，其模拟效果明显优于经典统计的线性回归模型。赵春雷等[14]对黄土高原坡面尺度土壤饱和导水率进行模拟，结果表明基于容重和土壤颗粒的状态空间方程可以解释坡面饱和导水率90%以上的变异。因此，与经典统计的回归方程相比，一阶自回归状态空间模型更能揭示空间上变量与变量之间的关系，是一种对变量进行模拟和预测的有效工具。

本文选取可能影响土壤容重的环境因素，包括黏粒和砂粒体积分数、有机碳含量、降水量、土地利用、海拔、坡度和坡向等作为潜在的预测变量，用一阶自回归状态空间模型及经典统计的多元逐步回归和传递函数方程分别对不同深度土层土壤容重的空间分布状况进行模拟，探索区域样带尺度不同土层深度土壤容重的空间分布特征及主要影响因素，以期为该区土壤水文过程模拟与预测提供水力参数。

1 材料与方法

1.1 研究区概况

黄土高原指日月山、贺兰山以东，太行山以西，秦岭以北，阴山以南的广大地区，面积64万km2。黄土高原年均降水量为150～800 mm，由东南向西北方向递减，多年平均降水量466 mm，多集中于6—9月份，占全年降水量的55%～78%。黄土高原属典型的温带大陆性季风气候，冬春季寒冷干燥，夏秋季温暖湿润，雨热同步，年平均气温3.6～14.3℃。该区光能资源丰富，年总辐射量为(5.02～6.70)×105J/cm2，呈东南向西北递增趋势。黄土高土壤主要为风成黄土，土层深厚，厚度多在50～200 m，部分地区深达300 m，土质疏松，易于侵蚀，且具有明显的地带性特征，由东南向西北，依次出现褐色土、塿土、黑垆土、黄绵土、灰钙土、栗钙土、棕钙土、棕漠土、风沙土等。黄土高原地貌类型主要有丘陵、塬、茆、沟、梁、阶地、平原、土石山地等[15]。

1.2 样品采集

2013年4月份在黄土高原南北方向布设一维样带，总长860 km。样带具有明显的环境梯度条件，年均降水量从南到北变化范围为620～400 mm，年均气温12.3～6.8℃，干燥度2.3～5.6，海拔高度383～1 529 m(图1)；土壤类型从南到北依次为塿土、黄绵土、风沙土；土地利用依次为农地(主要指一年生农作物地)、林地(包括天然/人工乔木林及灌木林地)、草地(包括人工草地和天然草地)[16]。根据地表景观和地形特征，沿样带每隔10 km布设一个样点，共86个样点，由南至北依次编号为1～86。在样带86个样点中，包含10个农地、54个林地和22个草地，为定量刻画土地利用对土壤容重的影响，将农地、林地和草地分别表示为0、1和2[7](图1)。王云强等[17]基于地统计分析发现，黄土高原区域尺度表层(0～5 cm)和底层(20～25 cm)土壤容重变程分别为448 km和300 km。根据变程可制定合理的采样策略，FLATMAN等[18]建议最佳的采样间距为变程的1/4～1/2之间。因此，研究黄土高原区域尺度土壤容重的空间变异性，表层土壤容重最低的采样间距应为112～224 km，底层应为75～150 km。本研究以典型黄土区为对象，结合黄土高原地区的数字地形图(DEM)和自然地理条件，采用高密度布点(样点间距约为10 km)策略，采样密度在其变程范围内，因此，样点设置合理，所得数据足以满足空间上的自相关要求。

图1 降水量、气温、相对高程和土地利用沿样带分布情况Fig.1 Distributions of precipitation, air temperature, relative elevation and land use types along south to north transect

利用RTK-GPS记录每个样点的经纬度坐标及海拔高度，利用罗盘仪测定每个样点的坡度和坡向，同时记录样点土地利用方式。在每个样点挖取40 cm深土壤剖面，用100 cm3的环刀分3层采集原状土壤样品(0～10 cm、10～20 cm、20～40 cm)。将原状土样带回实验室，利用烘干法测定土壤容重。同时采集不同土层深度的扰动土壤样品，风干后分别过0.25 mm和1.0 mm筛预处理，用于土壤黏粒、砂粒体积分数和有机碳含量测定[19]。

利用黄土高原68个国家气象观测站1951—2012年月值气象数据，包括月均降水量和气温，使用反距离加权插值方法将气象要素插值生成空间上连续分布的气象数据(空间分辨率为100 m×100 m)，然后提取每一样点的降水量和气温数据[16]。

1.3 分析方法

1.3.1 经典统计学空间变异分析理论

经典统计学进行空间变异分析时要求分析的空间变量在研究区域内符合随机分布特征，对变异特征的分析主要以变异系数Cv进行，其定义为

(1)

根据变异程度分级，Cv≤10%属于弱变异性，10%

1.3.2 状态空间方程

状态空间方程是基于马尔可夫(Markovian)过程的特征建立起来的，主要用于描述在同一地点采集到的j组数据Yj(xi)的随机过程。根据给定现有状态的数据，可以依据过去的过程，预测未来的状况。在这些系统中，过程状态包括用来预测将来所需的过去所有信息[20-21]。状态空间模型包含2个方程：状态方程和观测方程。系统在位置i处的状态通过状态空间方程与位置i-h(h=1，2，…，n)的状态方程发生联系。当h=1时，状态方程可表示为

Xi=φXi-1+wi

(2)

式中Xi——在位置i处的状态向量，通常指若干个变量在i处的值

φ——状态转化系数矩阵

wi——模型误差

φ是一个表示状态在位置i与i-1处的空间线性回归情况的p′×p′矩阵，wi也是一个普通的p′×p′协方差矩阵，且在空间位置上不相关、方差恒定呈正态分布且均值为0。状态向量以及转化系数矩阵中的参数可通过迭代回归求得。状态向量Xi通过与观测向量相联系，从而产生观测方程。观测方程为

Yi=AiXi+vi

(3)

式中vi——在空间位置上不相关，方差恒定呈正态分布，且均值为0的协方差矩阵

在观测方程中，通过观测矩阵Ai以及观测误差vi将系统的观测向量和系统的真实状态向量联系起来。状态方程和观测方程中2个误差项wi和vi相互独立。在应用过程中，Xi为系统的真实状态，Yi为测量值，但不一定能全部反映系统的真实状态，而只是系统状态的间接反映，它等于系统的真实状态再加上一个不能确定的误差。状态方程和观测方程的求解同时进行，状态向量以及转化系数矩阵φ中的参数可通过迭代回归求得[22]。

为消除不同量纲、变量自身变异和数值对模拟结果的影响，状态空间模拟要求对数据进行标准化，标准化公式为

zi=[Zi-(m-2s)]/(4s)

(4)

式中Zi——标准化前的变量zi——标准化后的变量m——标准化前变量的平均值

标准化后，各变量的均值为0.5，标准差为0.25。标准化后可以根据转换系数来反映变量之间的相对贡献率[23]。本文各种模拟分析中均采用标准化后的值。在评价状态空间模型的优劣时，统一使用决定系数(R2)和均方根误差(RMSE)来评价模型的优劣。

决定系数R2可揭示模型对样带土壤容重空间变异的解释量。R2值越大，表示模拟值所能解释土壤容重的变异越多。R2值越大，RMSE值越小，表示模拟效果越好。

1.4 数据处理

本研究中不同土层深度下土壤容重的统计特征值、Pearson相关分析、土壤容重的传递函数方程和多元逐步回归方程的建立均利用SPSS 13.0进行。土壤容重的状态空间模型在ASTSA软件下进行。

2 结果与分析

2.1 土壤容重、有机碳含量、黏粒和砂粒体积分数沿样带分布情况

图1显示了样点降水量、气温、相对高程和土地利用沿样带的分布状况。样点降水量和气温由南向北逐渐降低，这由研究区气候特征决定。1～10号样点相对高程较低，11～86号样点相对高程较高，这是因为1～10号样点分布在关中平原地区，而11～86号样点分布于黄土高原丘陵沟壑区，海拔高度较高，地形地貌复杂[24]。南北样带0～10 cm、10～20 cm和20～40 cm深度土壤容重的平均值分别为1.24、1.33、1.37 g/cm3(表1和图2)，表明随土层深度的增加，土壤容重逐渐增大。这可能与土壤有机碳含量有关，本研究中0～10 cm土壤有机碳含量高于10～20 cm土层，20～40 cm土壤有机碳含量最低(图2)，这与刘志鹏等[12]在区域尺度下的研究结果一致，有机碳含量在一定程度上影响着容重沿剖面的分布。WANG等[7]分析研究了黄土高原地区0～25 cm土壤容重分布状况，结果表明容重与有机碳含量显著负相关。BENITES等[25]通过研究巴西区域尺度0～100 cm土壤容重状况，也发现有机碳含量与容重存在显著的负相关关系。样带0～10 cm和10～20 cm土层深度土壤容重的变异系数分别为11.01%和11.07%，为中等程度变异；20～40 cm土层容重变异系数仅为8.82%，为弱变异。这主要是由于上层土壤相较于下层土壤更容易受到环境因素的影响[7]，因此，空间变异性较强。样带砂粒和黏粒体积分数沿样带分别呈现出逐渐增加和减小的趋势，这主要是由于黄土高原形成过程中风搬运颗粒逐级分选造成的[26]。

表1 不同土层深度下土壤容重统计特征值

注：N为服从正态分布。

2.2 容重与环境变量的Pearson相关分析

表2为所选变量与不同土层深度土壤容重(B)的相关系数。

结果表明，0～10 cm土壤容重与有机碳含量、黏粒体积分数、降水量和坡度显著负相关，与土地利用、海拔高度和坡向无显著相关关系；10～20 cm土壤容重与有机碳含量、黏粒体积分数、降水量和土地利用显著负相关，与海拔高度、坡度和坡向无显著相关关系；20～40 cm土壤容重与黏粒体积分数、降水量、土地利用和海拔高度显著负相关，与有机碳含量、坡度和坡向无显著相关关系；此外，各土层容重均与砂粒体积分数呈极显著正相关关系，表明砂粒体积分数越高，土壤容重越大，这与WANG等[7]在黄土高原的研究结果一致。本研究只选择与不同土层深度土壤容重显著相关的变量分别建立状态空间方程、传递函数方程和多元逐步回归方程，筛选出不同方法最优的预测方程，并比较不同方程对样带不同深度土壤容重的模拟效果。因此，对于0～10 cm土层，以有机碳含量、黏粒和砂粒体积分数、降水量和坡度作为预测变量；对于10～20 cm土层，以有机碳含量、黏粒和砂粒体积分数、降水量和土地利用作为预测变量；对于20～40 cm土层，以黏粒和砂粒体积分数、降水量、土地利用和海拔高度作为预测变量。

2.3 土壤容重的多元逐步回归方程和传递函数方程

多元逐步回归方程作为经典的统计预测方法已得到广泛应用。本文采用经典统计的多元逐步回归方法，将与土壤容重显著相关的变量作为自变量对因变量土壤容重进行线性回归模拟。不同土层深度多元逐步回归方程均经过显著性检验，且均在p<0.01 水平上显著。在3个多元逐步回归方程中均包含了黏粒和砂粒体积分数，表明在研究区，黏粒和砂粒体积分数是影响土壤容重的重要因素(表3)。一般来说，随着砂粒体积分数的增加和黏粒体积分数的减小，由于土壤总孔隙度的减小，容重呈现增加的趋势。 WANG等[7]运用人工神经网络和线性回归的方法模拟预测了黄土高原区域尺度表层土壤容重状况，研究结果表明，包含黏粒和砂粒体积分数等因素在内的线性回归模型可以较好地预测区域尺度下0～25 cm深度土壤容重。本研究中，不同土层深度的逐步回归方程决定系数R2为0.444～0.493，均方根误差RMSE为0.177～0.185 g/cm3，表明在当前研究尺度下，模型变量所能解释的土壤容重空间变异仅为44.4%～49.3%，即模拟精度较低。

图2 不同土层深度土壤容重、有机碳含量、黏粒和砂粒体积分数沿样带分布情况Fig.2 Distributions of soil bulk density, soil organic carbon content and clay and sand contents in different soil layers along south to north transect

表2 不同土层深度土壤容重与各变量的相关系数

注：*表示相关关系在p<0.05水平上具有显著性(双尾检验)；** 表示相关关系在p<0.01水平上具有显著性(双尾检验)。

为了提高统计方程的模拟精度，进一步将与土壤容重显著相关的变量以及其4种常见转换(指数、对数、倒数和平方转换)作为自变量再次进行逐步回归分析，以期得到对土壤容重空间变异具有更强解释能力的传递函数方程，并选出各土层模拟精度最优的传递函数方程(表3)。不同土层深度传递函数方程均经过显著性检验，且均在p<0.01水平上显著。结果表明，不同土层深度包含了变量及其常见转换的传递函数方程在模拟精度上得到了不同程度的提高。传递函数方程分别可以解释0～10 cm、10～20 cm和20～40 cm土层深度土壤容重空间变异的50.4%、58.7%和54.0%，RMSE分别为0.174、0.154、0.165 g/cm3，均低于多元逐步回归方程的RMSE。与多元逐步回归方程相比，0～10 cm、10～20 cm和20～40 cm土层深度传递函数方程模拟精度分别提高了6.3%、9.5%和32.2%。这主要是因为传递函数方程包含了不同变量的转换，可以更好地刻画变量与容重之间的关系[7]。

表3 不同土层深度下土壤容重模拟的状态空间方程、传递函数方程和多元逐步回归方程

2.4 土壤容重的状态空间模拟

将与不同土层深度土壤容重显著相关的变量及其不同因素组合分别引入状态空间模型中，以期找出能够模拟样带不同土层深度土壤容重的最佳因素组合。所有状态空间方程均经过显著性检验，且均在p<0.01水平上显著。表3为不同土层深度土壤容重模拟效果最好的状态空间方程。结果表明，不同土层深度最优状态空间方程包含变量不同，0～10 cm为有机碳含量、黏粒和砂粒体积分数；10～20 cm为有机碳含量、黏粒和砂粒体积分数和降水量；20～40 cm为黏粒和砂粒体积分数、降水量和土地利用。由决定系数R2可知，最优状态空间方程分别可以解释样带0～10 cm、10～20 cm和20～40 cm土层深度土壤容重空间变异信息的69.1%、72.7%和92.3%。因此，基于黏粒和砂粒体积分数、降水量和土地利用的状态空间方程对样带20～40 cm容重模拟效果最好，且95%置信区间最窄。图3显示了不同土层深度最优状态空间方程以及土壤容重实测值和模拟值的分布情况。从图3可以看出，20～40 cm土层土壤容重的模拟值与实测值重合性很好，且置信区间相对较窄，置信区间内包括了所有的实测值，而0～10 cm和10～20 cm土层模型的置信区间相对较宽且有部分实测值在置信区间以外，其模拟效果逊于20～40 cm，这可能与不同土层深度土壤容重的变异强弱有关。0～20 cm深度土壤容重为中等程度变异，20～40 cm为弱变异(表 1)，即表层土壤容重比深层土壤容重可能受外界因素影响更加复杂，而依据有限的变量表层模拟精度要低于深层。

图3 不同土层深度模拟土壤容重最优的状态空间模型Fig.3 Optimal state-space models for simulating soil bulk density in different soil layers

图4 不同土层深度最优状态空间模型中各变量对土壤容重的相对贡献率Fig.4 Relative contribution of each variable in optimal state-space model for soil bulk density in different depths

根据模拟结果的转换系数，可计算各状态空间方程中i-1处变量对相邻位置处土壤容重的相对贡献率(图4)。在0～10 cm土层，有机碳含量、黏粒和砂粒体积分数对容重的相对贡献率分别为16.1%、9.4%和14.8%；在10～20 cm土层，有机碳含量、黏粒和砂粒体积分数和降水量对容重的相对贡献率分别为34.2%、5.1%、9.1%和23.7%；在20～40 cm土层，黏粒和砂粒体积分数、降水量和土地利用对容重的相对贡献率分别为28.8%、40.4%、12.6%和17.3%。因此，在0～20 cm土层，有机碳含量对土壤容重空间分布的相对贡献率最大，是影响土壤容重空间分布的主要因素。然而，10～20 cm土层中有机碳含量的相对贡献率高于0～10 cm，这主要是由于表层土壤容重的影响因素相对复杂，有机碳含量对其变异的刻画能力相对较弱[25]。在20～40 cm土层，有机碳含量没有被引入最优的状态空间方程，这主要是由于深层土壤中有机碳含量较低且变异较小(图2)，这与KAUR等[27]的结论一致。值得注意的是，在表层0～10 cm，影响容重空间分布主要是土壤因素，而在10～20 cm和20～40 cm土层，气候因素和土地利用应当被考虑。

不同土层深度最优的状态空间方程的RMSE分别为0.118、0.114、0.028 g/cm3，远低于多元逐步回归方程和传递函数方程。因此，由R2和RMSE可知，在3个土层中，状态空间方程对样带土壤容重空间分布的模拟效果均优于经典统计的多元逐步回归方程和传递函数方程，可解释土壤容重空间变异的69.1%～92.3%。这主要是因为经典统计的回归分析不考虑采样点的空间位置，假设样点与样点之间是相互独立的，即在整个研究区域上把土壤容重与其他因素的相关性看作是单一的，而忽略容重与其他因素的关系在空间位置上的变异性[28]，而状态空间模型则考虑了变量之间相关性的空间变异。此外，状态空间方程考虑了模型的误差和观测误差，进而通过迭代算法调整模型参数，最终获取预测精度较高的模型参数[12]。因此，相对于经典统计的回归分析，状态空间模型更能揭示其他因素与土壤容重之间的关系，是一种对土壤容重进行预测和模拟的有效工具。

3 结论

(1)黄土高原南北样带0～20 cm深度土壤容重变异程度为中等变异，20～40 cm为弱变异。

(2)状态空间方程的模拟效果优于经典统计的多元逐步回归方程和传递函数方程，基于黏粒和砂粒体积分数、降水量和土地利用的状态空间方程对20～40 cm深度土壤容重的模拟效果最好，决定系数R2在0.9以上，可用于田间条件下区域土壤容重分布特征的预测。

(3)样带不同土层深度土壤容重的影响因素不同，0～10 cm主要为有机碳含量、黏粒和砂粒体积分数，10～20 cm为有机碳含量、黏粒和砂粒体积分数和降水量，20～40 cm为黏粒和砂粒体积分数、降水量和土地利用。

1 郑纪勇, 邵明安, 张兴昌. 黄土区坡面表层土壤容重和饱和导水率空间变异特征[J]. 水土保持学报, 2004, 18(3):53-56. ZHENG Jiyong, SHAO Ming’an, ZHANG Xingchang. Spatial variation of surface soils bulk density and saturated hydraulic conductivity on slope in loess region[J]. Journal of Soil and Water Conservation, 2004, 18(3):53-56. (in Chinese)

2 ALEXANDER E B. Bulk densities of California soils in relation to other soil properties[J]. Soil Science Society of America Journal, 1980, 44(4):689-692.

3 刘晚苟, 山仑, 邓西平. 不同土壤水分条件下土壤容重对玉米根系生长的影响[J]. 西北植物学报, 2002, 22(4):831-838. LIU Wan’gou, SHAN Lun, DENG Xiping. Effects of soil bulk density on the growth of maize root system under different water condition[J]. Acta Botanica Boreali-occidentalia Sinica, 2002, 22(4):831-838. (in Chinese)

4 王辉, 王全九, 邵明安. 表层土壤容重对黄土坡面养分随径流迁移的影响[J]. 水土保持学报, 2007, 21(3):10-13. WANG Hui, WANG Quanjiu, SHAO Ming’an. Effect of soil bulk density on soil nutrient in runoff from loess slope[J]. Journal of Soil and Water Conservation, 2007, 21(3):10-13. (in Chinese)

5 韩光中, 王德彩, 谢贤健. 中国主要土壤类型的土壤容重传递函数研究[J]. 土壤学报, 2016, 53(1):93-102. HAN Guangzhong,WANG Decai, XIE Xianjian. Pedotransfer functions for prediction of soil bulk density for major types of soils in China[J]. Acta Pedologica Sinica, 2016, 53(1):93-102. (in Chinese)

6 王巧利, 林剑辉, 许彦峰. 基于BP神经网络的土壤容重预测模型[J]. 中国农学通报, 2014, 30(24):237-245. WANG Qiaoli, LIN Jianhui, XU Yanfeng. The model for predicting soil bulk density based on the BP neural network[J]. Chinese Agricultural Science Bulletin, 2014, 30(24):237-245. (in Chinese)

7 WANG Yunqiang, SHAO Ming’an, LIU Zhipeng, et al. Prediction of bulk density of soils in the Loess Plateau region of China[J]. Surveys in Geophysics, 2014, 35(2):395-413.

8 BALLAND V, POLLACCO J A P, ARP P A. Modeling soil hydraulic properties for a wide range of soil conditions[J]. Ecological Modelling, 2008, 219(3-4):300-316.

9 KAUR R, KUMAR S, GURUNG H P. A pedo-transfer function (PTF) for estimating soil bulk density from basic soil data and its comparison with existing PTFs[J]. Australian Journal of Soil Research, 2002, 40(5):847-858.

10 HAN Guangzhong, ZHANG Ganlin, GONG Zitong, et al. Pedotransfer functions for estimating soil bulk density in China[J]. Soil Science, 2012, 177(3):158-164.

11 汪懋华. “精细农业”发展与工程技术创新[J]. 农业工程学报, 1999, 15(1):1-8. WANG Maohua. Development of precision agriculture and innovation of engineering technologies[J]. Transactions of the CSAE, 1999, 15(1):1-8. (in Chinese)

12 刘志鹏, 邵明安, 王云强. 区域尺度下黄土高原土壤全钾含量的空间模拟[J]. 农业工程学报,2012,28(22): 132-140． LIU Zhipeng, SHAO Ming’an, WANG Yunqiang. Spatial simulation of soil total potassium in regional scale for Loess Plateau region[J]. Transactions of the CSAE, 2012, 28(22):132-140. (in Chinese)

13 贾小旭, 邵明安, 魏孝荣, 等. 黄土高原北部草地表层土壤水分状态空间模拟[J]. 农业工程学报,2010,26(10):38-44． JIA Xiaoxu, SHAO Ming’an，WEI Xiaorong, et al. State-space simulation of soil surface water content in grassland of northern Loess Plateau[J]. Transactions of the CSAE, 2010, 26(10):38-44. (in Chinese)

14 赵春雷, 邵明安, 贾小旭. 黄土高原北部坡面尺度土壤饱和导水率分布与模拟[J]. 水科学进展, 2014, 25(6):806-815. ZHAO Chunlei, SHAO Ming’an, JIA Xiaoxu. Distribution and simulation of saturated soil hydraulic conductivity at a slope of northern Loess Plateau[J].Advances in Water Science, 2014, 25(6):806-815. (in Chinese)

15 杨勤业, 张伯平, 郑度. 关于黄土高原空间范围的讨论[J]. 自然资源学报, 1988, 3(1):9-15. YANG Qinye, ZHANG Boping, ZHENG Du. On the boundary of the Loess Plateau[J]. Journal of Natural Resources, 1988,3(1):9-15. (in Chinese)

16 JIA X, SHAO M, ZHANG C, et al. Regional temporal persistence of dried soil layer along south-north transect of the Loess Plateau, China[J]. Journal of Hydrology, 2015, 528: 152-160.

17 王云强，邵明安. 黄土高原地区土壤干层的空间分布与影响因素[D]. 杨凌：中国科学院研究生院, 2010. WANG Yunqiang, SHAO Ming’an.Spatial distribution and influence factors of dried soil layers across the Loess Plateau[D]. Yangling: University of Chinese Academy of Sciences, 2010. (in Chinese)

18 FLATMAN G T, YFANTIS A A. Geostatistical strategies for soil sampling: the survey and the census [J]. Environmental Monitoring and Assessment, 1984, 4(4): 335-349.

19 ZHAO C, SHAO M, JIA X, et al. Using pedotransfer functions to estimate soil hydraulic conductivity in the Loess Plateau of China[J]. CATENA, 2016, 143:1-6.

20 WENDROTH O, REUTER H I, CHRISTIAN K K. Predicting yield of barley across a landscape: a state-space modeling approach[J]. Journal of Hydrology, 2003, 272(1): 250-263.

21 TIMM L C, REICHARDT K, OLIVEIRA J C M, et al. State-space approach for evaluating the soil-plant-atmosphere system[R]. Lectures Given at the College on Soil Physics Trieste, 2003： 387-428.

22 SHUMWAY R H, STOFFER D S. An approach to time series smoothing and forecasting using the EM algorithm[J]. Journal of Time Series Analysis, 1982, 3(4): 253-264.

23 WENDROTH O, JURSCHIK P, KERSEBAUM K C, et al. Identifying, understanding, and describing spatial processes in agricultural landscapes—four case studies[J]. Soil & Tillage Research, 2001, 58(3-4): 113-127.

24 WANG Y, SHAO M, LIU Z, et al. Investigation of factors controlling the regional-scale distribution of dried soil layers under forestland on the Loess Plateau, China[J]. Surveys in Geophysics, 2012, 33(2): 311-330.

25 BENITES V M, MACHADO P L O A, FIDALGO E C C, et al. Pedotransfer functions for estimating soil bulk density from existing soil survey reports in Brazil[J]. Geoderma, 2007, 139(1-2):90-97.

26 ZHAO C, SHAO M, JIA X, et al. Particle size distribution of soils (0～500 cm) in the Loess Plateau, China[J]. Geoderma Regional, 2016, 7(3):251-258.

27 KAUR R, KUMAR S, GURUNG H P. A pedo-transfer function (PTF) for estimating soil bulk density from basic soil data and its comparison with existing PTFs[J]. Australian Journal of Soil Research, 2002, 40(5):847-858.

28 JIA X, SHAO M，WEI X, et al. Estimating total net primary productivity of managed grasslands by a state-space modeling approach in a small catchment on the Loess Plateau，China[J]. Geoderma, 2011, 160(3-4):281-291.

Variation Analysis and Simulation of Soil Bulk Density within Different Soil Depths along South to North Transect of Loess Plateau

YI Xiaobo1SHAO Ming’an1,2ZHAO Chunlei1ZHANG Chencheng1,2

(1.CollegeofResourcesandEnvironment,NorthwestA&FUniversity,Yangling，Shaanxi712100,China2.StateKeyLaboratoryofSoilErosionandDrylandFarmingontheLoessPlateau,NorthwestA&FUniversity,Yangling，Shaanxi712100,China)

To provide soil hydraulic parameters for simulating soil hydrological processes, classical statistics were used to characterize the spatial distribution of soil bulk density (BD) within different soil depths (0～10 cm, 10～20 cm and 20～40 cm) along south to north transect of the Loess Plateau. Furthermore, the multiple stepwise regressions, pedotransfer functions and the first order autoregressive state-space models were applied to simulate the distribution of BD. The magnitude of BD variability in the 0～20 cm soil layer was moderate according to the coefficient of variations, while the BD variability in the 20～40 cm soil layer was weak. The key factors affecting the spatial distribution of BD differed in different soil depths. Soil organic carbon, clay and sand contents were the key factors in the 0～10 cm soil layer; soil organic carbon content, clay and sand contents and precipitation were the key factors in the 10～20 cm soil layer; while in the 20～40 cm soil layer, clay and sand contents, precipitation and land use were the key factors to affect the spatial distribution of BD along south to north transect. State-space models were consistently more effective than multiple stepwise regression functions and pedotransfer functions for estimating spatial distribution of BD. State-space model that included clay and sand contents, precipitation and land use factors showed the best simulation result, and the combination of such variables explained 92.3% of the total variation of BD. State-space models were recommended for studying spatial relations between soil bulk density and other variables on the Loess Plateau.

bulk density; state-space; pedotransfer function; simulation; spatial variation

10.6041/j.issn.1000-1298.2017.04.026

2016-08-14

2016-09-28

国家自然科学基金项目(41530854、41501233)

易小波(1986—)，男，博士生，主要从事环境生态研究，E-mail： abobobi@163.com

邵明安(1956—)，男，研究员，博士生导师，主要从事土壤物理研究，E-mail： shaoma@igsnrr.ac.cn

S152.7

A

1000-1298(2017)04-0198-08