基于回归克里格的河北省土壤有机质时空变化特征研究

李 可，杨 勇，刘亚军，谭秀丽，杨 雪，李恩光，吴若景

（华中农业大学资源与环境学院，农业农村部长江中下游耕地保育重点实验室，湖北 武汉 430070）

土壤是人类进行各项生产和赖以生存的基地，土壤质量变化与质量高低关系到人类的生存状况，而土壤有机质（SOM）是构成土壤有机无机复合胶体的核心物质，是土壤的重要组成部分，对土壤结构、物理性状、化学性状和生态过程有着重要的影响与作用［1-2］。所以，研究SOM的时空变化特征，对协调耕地土壤肥力，提高土壤质量水平，揭示土壤有机质空间分布格局，了解土壤质量动态变化，优化农业种植结构和施肥管理具有科学指导作用［3-4］。

近年来，由于GIS技术和地统计学的快速发展，学者们开始利用两者相结合来研究各个尺度上SOM时空变化特性，为测土培肥等技术实施提供理论支持。如田块尺度上，李国良等［5］探究了广东荔枝园土壤养分肥力时空变化特征，提出该地区土壤养分肥力发生了很大变化，多种肥力属性有明显下降趋势，但有机质变化不大；县市尺度上，王少贺等［6］对台安县SOM的时空变化特征的研究表明：25年来，该县耕地土壤有机质平均含量下降了2.65 g/kg，整体已呈现出下降的趋势，这与人们长期的“重种轻养、只用不养”有关；鲍思屹等［7］研究了建瓯市近40年来竹林土壤有机质的变化特征，得出建瓯市毛竹林土壤有机质Ⅰ级与Ⅲ级土壤面积分别增加了4.68%和17.3%，而Ⅱ与Ⅳ级土壤面积分别减少了17.4%和4.63%，人类活动对竹林土壤有机质的影响程度有所增强；省级尺度上，赵明松等［8］以江苏省为研究对象，对比分析了近26年全省土壤有机质的时空变异特征，结果表明：SOM含量空间分布呈现出北增南减，沿江平原增，宁镇丘陵减，滨海平原基本持平的空间格局，增加幅度由北向南逐渐减小；全国尺度上，任意等［9］对我国东北、西北、华北、西南、中南和华东不同地区土壤养分的差异及变化趋势进行了研究，结果指出在1986～1997年中南、华北区有显著的上升趋势，华东、东北区有下降趋势，西北、西南区变化平稳，1998～2006年各个区域的土壤有机质变化基本平稳。以上研究均表明，SOM在各空间尺度上存在明显的时空分异，这对分析土壤肥力对人类活动的响应、指导农业生产、实现农业生态的可持续发展有重要意义。

从地理条件上讲，河北省地形地貌复杂，是全国唯一兼有高原、山地、丘陵、平原、湖泊和海滨的省份，土壤类型较多。从经济人口条件上讲，河北省总人口数全国排第6，约为7 185.4万（第六次全国人口普查数据）。而且河北省是农业大省，全省34.60%的面积为耕地，其粮食产量在全国位居前列。因此，本文选取了河北省作为研究区域，并采用GIS空间分析、地统计分析法以及回归克里格空间插值方法，对1980～2010年间河北省土壤有机质的时空变化进行了研究，为预测粮食产量、优化施肥管理、保护生态平衡和耕地的可持续发展等提供了科学依据。

1 材料与方法

1.1 研究区概况

河北省简称“冀”，位于东经113°27′～119°50′，北纬 36°05′～42°40′，地处华北地区，北靠燕山，南望黄海，西依太行，东临渤海，面积约18.8万km2。河北省地势自西北向东南倾斜，东部有太行山脉，全省依据地形地貌具体可分为河北平原区、坝上高原区、燕山和太行山地区3大部分。全省最高海拔2 855 m，平均海拔551 m，大部分地区还是属于平原丘陵。总的来说，河北省是全国唯一兼有高原、山地、丘陵、平原、湖泊和海滨的省份，其土壤类型按发生分类共分为21种土类、55种亚类，整体呈带状分布，变化较大，极具研究意义。

1.2 数据来源

1978～1980年（以下均表示为1980年）河北省表层土壤数据来自全国第二次土壤普查时期，《河北土种志》［10］中记录的334个典型土壤剖面，每个剖面的坐标位置通过《土种志》对典型剖面的地理位置描述与高清遥感影像结合得出。2009～2011年（以下均表示为2010年）表层土壤数据为国家科技基础性工作专项“我国土系调查与《中国土系质》［11］编制”中河北省162个典型土壤剖面数据。两个时期具体样点分布如图1所示，不同的土地利用类型上均有样点分布，1980年时耕地、林地、草地、水域、未利用地与城乡、工矿、居民用地6类土地利用类型上分别分布了195、30、31、10、5、63个土壤样点，2010年时则分别分布了68、41、35、6、4、8个土壤样点，其中土系调查72个典型土壤剖面均在第二次土壤普查时期的典型剖面附近采集。

本文采用的辅助变量为海拔高程、曲率、归一化植被指数、坡度、土地利用类型和土壤类型。因此，除了土壤剖面数据以外，还收集了研究区域1 km分辨率的数字高程图（DEM）、2010年植被指数（NDVI）空间分布图、1980和2010年土地利用数据（由Landsat TM/ETM遥感图像人工目视解译生成）等环境因子数据，来源于中国科学院资源环境科学数据中心（http：//www.resdc.cn）。河北省1∶100万土壤类型图为ArcGIS数字化数据，来源于国家科技基础条件平台—国家地球系统数据共享服务平台-土壤科学数据中心（http：//soil.geodata.cn）。坡度和曲率数据在ArcGIS 10.2软件中由DEM数据提取而来。结果分析时涉及统计数据来源于1995和2015年《河北农村统计年鉴》、1985和2011年《河北经济年鉴》，涉及部分土地利用数据来自上述土地利用遥感数据计算得出的结果。

图1 研究区域土壤样点分布图

1.3 研究方法

1.3.1 回归克里格插值（Regression Kriging，RK）

回归克里格的基本步骤原理：首先，对辅助环境变量和目标变量之间的回归关系进行计算；然后，对空间样点上的趋势项进行分离，用实际值减去趋势项值得到残差，对残差进行线性无偏估计处理；最后，将回归分析的趋势项同残差插值结果相加，从而得到目标变量的空间估测值［12-13］。假设目标变量为 z（s1），z（s2），…，z（sn），si=（xi，yi），i=1，2，…，n，si为实测点的坐标值，n为实测点的数目。则预测点s0的值z（s0）可表示为：

式中，趋势项m（s0）通常用一元或多元线性回归进行拟合，残差项ε（s0）用普通克里格插值。由上式可以看出，回归克里格方法能够综合多个确定性的辅助环境变量（如曲率、坡度、高程等）对目标变量进行解释预测［14］，在样点较稀疏的区域得到精度较高的预测结果［15］。鉴于本研究的区域面积较大，样点数量相对较少，因此，选用此方法对土壤养分进行空间插值。

1.3.2 哑变量引入

土地利用和土壤类型分布图均为范畴型变量，与DEM、NDVI等连续型变量不同的是，其Value值分别为标识代码和实际数值，两者无法同时直接用于进行回归分析，因此引入哑变量值［16］。哑变量为虚拟变量，常用于处理定性因子或分类变量，一般取值为0或1［17］。哑变量的定义为对于等级（定性）数据α，用变量 （α，β）表示为：

这种方法叫做定性因子二值化展开，变量（α，β）就称为哑变量。因定性变量为取0或1的数值向量，便于用数值方法进行处理［18］。本文将土地利用数据按照一级分类分为耕地、林地、草地、水域、未利用地与城乡、工矿、居民用地6类，土壤类型数据按照发生分类中的土纲分为9类，然后按上述方法进行哑变量换算，符合相应类型即为1，不符合即为0，最后将哑变量和其它连续变量结合做出回归分析。另外，本研究是直接利用全国第二次土壤普查图件河北省土壤图提取土纲进行哑变量转换，通过地图套合与2010年土系调查的土壤样点对应获得样点发生分类名称，而非将2010年土壤样点的系统分类的土壤类型直接对照为发生分类土壤类型。

1.4 分析工具

回归克里格的计算与插值采用Matlab R2014a软件和ArcGIS 10.2的地统计模块来完成；空间数据的处理与制图依靠ArcMap 10.2来完成；基本的统计分析由SPSS 19和Excel 2010软件完成。

2 结果与分析

2.1 1980～2010年全省土壤有机质的时间变化

表1为第二次全国土壤普查土壤有机质分级标准［19］，表2为1980和2010年土壤典型剖面点描述性统计结果。由表2可知，1980和2010年土壤有机质含量范围分别为0.80～164.50和2.60～208.50 g/kg，均值分别为22.02和18.38 g/kg，属于III和IV级标准（对比表1），这说明从总体上看，近30年来SOM含量有一定程度的下降。而且，1980和2010年SOM含量变异系数分别为1.10和1.07，均大于1，属于高度变异，这表示在该时期内河北省的该类土壤有机质含量的空间分布差异较大，受人为影响严重。同时，对比两时期内SOM含量的变异系数可发现，近30年来土壤有机质含量的变异系数有所减小，这表示随着耕作制度、施肥量变化、施肥方式和田间管理措施等人为因素的影响，土壤有机质含量的空间变异特征表现为逐渐减弱的趋势。

表1 第二次全国土壤普查土壤有机质分级标准 （g/kg）

表2 1980和2010年土壤样点SOM描述性统计数据

表3为河北省1980和2010年各地市SOM含量变化结果。由表3来看，1980年时石家庄、保定、张家口、承德地区的SOM含量均高于全省平均值，2010年时只有承德地区SOM含量高于全省平均值。从两个时期对比变化来看，保定和张家口地区SOM含量在1980年时明显高于其他地区，但随着近30年来SOM含量变化，其整体呈现降低趋势，且差值很大，均从II级标准直跌至IV级标准。石家庄、承德、邢台3地SOM含量同样呈降低趋势。其中，石家庄和承德SOM含量降低量适中，邢台地区SOM含量降低量较小。其他地市中，除邯郸地区SOM含量无明显变化外，剩余地区SOM含量均呈增长状态，且增量都较小。

造成上述土壤有机质时间变化的原因：张家口和承德均位于北部，此地以棕壤和栗钙土为主，腐殖质层较厚，所以两地区在1980年时土壤有机质较高。而在1990～2008年间，由于大量退耕还林还草导致张家口和承德两地常用耕地面积分别减少了220 362、79 709 hm2，但近30年中两地区总耕地面积增量百分比分别只有0.65%、-0.11%。这表明至2010年，张家口和承德地区在总耕地面积变化极小的前提下，大量非耕地转变为耕地，所以SOM呈下降趋势。同时，张家口地区耕地压力指数增大［20］，至2010年，其耕地肥料使用量虽增加至113.50 kg/hm2，但相对其它市区来说其施肥量最低，所以土壤养分无法得到及时补充，也会导致SOM含量降低。承德地区耕地压力指数也在增加，加之近些年各类土地利用强度加大，农业技术及管理水平相对较低，因此土壤有机质在大幅度消耗后无法及时补足，SOM含量降低［20-21］。而保定、石家庄、邢台地区，1980年SOM含量较高样点所处的区域不属于耕地范围，所以不会对其进行施肥等处理，随着土壤养分流失，其SOM含量值同样也就会降低。而其它地区可能由于土壤类型、20世纪80年代以前本省耕地管理问题等的影响［22］，在1980年时土壤有机质相较于上述地市偏低，但经过近30年施肥、秸秆还田、耕作制度调整等措施，促进了土壤有机质的提升。

表3 1980～2010年河北省各地市土壤有机质含量对比

2.2 回归克里格插值结果

利用回归克里格插值方法对两时期样点数据进行插值。在插值过程中分别随机选取20%样点作为验证点进行精度评价，不参与插值，回归分析及精度评价结果见表4。表4中R2、F、P为回归方程评价指标，RMSE为插值精度评价指标。由表4可得，1980和2010年SOM的分析P值均小于0.01，属于极显著水平。虽然SOM在两时期内的决定系数和统计量较低，但其最终两时期SOM含量的结果精度值RMSE均小于其对应均值的1/2（对比表2），所以整体插值精度也属于误差接受范围。特别的是，两时期中1980年SOM含量的RMSE值更小，可能是由于此时期用于插值的样点数目为2010年的2倍，从而使得插值结果的精度更高。

2.3 1980～2010年全省土壤有机质的空间变化

图2为1980～2010年SOM含量空间分布图。由图2可知，1980和2010年SOM含量整体上空间分布均是西北高东南低，并且整体呈现由西北向东南降低的趋势，至2010年SOM含量东南与西北差异减小。由图2a可得，1980年时SOM含量达到I、II级标准的高值区域主要分布在承德西北部、张家口东北部和南部、保定西北部，V、VI级的低值区域主要在衡水、邢台东部、唐山北部。且SOM值大多集中在III级和IV级标准区间内，分别占全省面积的16.35%和33.85%，主要分布在河北省南部平原地区。由图2b分析得，2010年SOM含量高值区域主要在承德、张家口西北部、保定北部等地，而低值区域面积较小，分布零散。此时，SOM含量同样集中在III级和IV级标准区间内，分别占全省面积的47.82%和20.73%，主要分布在南部平原地区、唐山地区。

图3为两时期SOM含量的空间变化分布图，白色和浅灰色区域代表SOM含量减少，深灰色和黑色区域代表其含量增加。从整体上看，SOM含量整体上呈现西北减东南增的趋势，增幅总体上由西北向东南依次增加；SOM含量增加的面积占全省面积的65.54%，主要分布在南部平原、唐山、秦皇岛地区，其中增幅在0～10 g/kg间占比最高，约占全省面积的36.29%。从局部看，SOM减少区域主要位于1980年SOM高值分布地区，SOM增值较高的区域则对应处于1980年SOM低值分布区域。同时，SOM含量增幅最大（＞20 g/kg）的区域在承德东部、唐山北部和秦皇岛，降幅最大（＜-20 g/kg）的区域在承德西北部、石家庄和保定西部。造成上述变化的原因是由于降幅较大的地区在1980年时SOM含量就较高，但多为草地、林地，基本不属于耕地。其中，位于承德西北部的围场、丰宁地处河北坝上地区，至2005年时，此处放牧、樵采过度，坝上草场面积已由20世纪50年代初的86.7万hm2减少到现在的22.6万hm2，草场覆盖度由90%降低到44%，森林面积也呈明显减少趋势，土地沙漠化严重，导致该区域土壤肥力流失严重，SOM含量大幅度降低［23］。而石家庄和保定西部地处太行山区，属于石质山区，由于过度追求开垦种植、砍伐森林，导致森林覆盖率低，生态环境脆弱，土层养分涵养差，SOM含量显著降低［24］。而原来的低值区域大多处于耕地范围，自20世纪90年代起，河北省开始大力推广绿肥种植、秸秆还田等措施，目前河北省东南部平原区域基本实现100%秸秆还田率，使得该区域土壤有机质状况得到较好改善［21，25］。

表4 1980和2010年河北省SOM回归分析及插值精度结果

图2 1980～2010年SOM空间分布图

图3 1980～2010年SOM空间变化分布图

3 讨论

由土壤有机质质量的演化过程可知，1980～2010 年土壤有机质含量的时空变化主要与自然和社会经济发展有关，自然过程的影响因素主要涉及气候、温度、降水等，但这些因素的变化从近30年来看，总体波动幅度较小，对SOM的累积与变化影响不大，甚至可忽略不计。因此，本文主要侧重分析社会经济因素对土壤有机质变化的影响。由李继明等［26］的研究可知，绿肥与化肥长期配合施用会使得土壤有机质、全氮和全磷均有所积累，其积累的量与肥料施用量及有机肥种类有关。所以，接下来本研究主要分析土壤有机质含量的时空变化同肥料施用量间的关系。

1980～2010年河北省化肥施用量（折纯量）由73.10万t增加至322.86万t，耕地面积由664.79万hm2减少至598.89万hm2，平均用量由109.96 kg/hm2增加至539.10 kg/hm2，其中氮肥、磷肥、钾肥和复合肥的用量有不同程度的增加（图4）。1985～2010年间，全省氮肥和复合肥用量增加较大，分别由109.79、24.38 kg/hm2增加至255.64、159.63 kg/hm2；磷肥和钾肥的用量增加较平稳。同时，1980～2010年期间，全省粮食（稻谷、小麦、玉米、大豆、薯类）总产量由1 575万t增加至2 976万t，粮食平均产量由2.15 t/hm2增加至4.88 t/hm2，这主要得益于种植高产作物及使用更多的化肥、农药、地膜等［20］。同时，由于化肥使用量的大幅度提升，全省耕地土壤肥力质量也有一定提升［22］。而对比SOM空间变化分布图也可发现，绝大多数SOM含量增长的区域属于耕地范围，与上述分析结果相符。

图4 1980～2010年肥料施用总量和平均量

由图5知，近些年来唐山、邯郸、秦皇岛等地化肥平均施用量提升较高，SOM含量增幅也较高；张家口耕地肥料施用提升最少，SOM含量增幅则较小。这表明耕地化肥施用量、耕作制度及种植结构等社会经济活动与土壤肥力质量变化总体上呈正相关，对SOM水平的提升具有一定正效应。其中石家庄、保定、承德地区化肥施用增量居中，但其SOM增幅较小，可能是由于肥料利用率低、耕地管理差导致的，也可能是由于1980年时此处高值采样点较多且过于集中，导致其代表性较差引起的。曾招兵等［27］分析得出广东省随着化肥的投入增加了作物产量，使得作物以根系或秸秆废弃物形式归还土壤的数量相应增长，耕层土壤有机质含量也相应增加；赵吉霞等［28］研究表明1985～2006年云南墨江县土壤全氮含量随着氮肥施用量增加而产生了较小幅度的增加。上述对不同地区的研究结果表明，河北省不同地市SOM含量总体上可能随化肥施用量增加而增加，但由于各地市消耗不同，具体增量不一样，部分地市增量甚至为负值，导致河北省有机质总量降低。

图5 河北省各地市化肥平均施用增量与SOM平均增量关系

4 结论

从时间上看，1980～2010年间河北省SOM含量总体呈降低趋势，平均降低3.6 g/kg。从空间上看，近30年全省SOM含量整体上均呈现西北减东南增的趋势，增幅总体上由西北向东南依次增加；且其增加的面积占全省面积的65.54%，主要分布在南部平原地区和唐山、秦皇岛等地。上述土壤有机质变化格局的形成主要与施用化肥、实施秸秆还田等人类活动有关。因此，要想改善土壤肥力下降的问题，就应坚持合理施用化肥、实施秸秆还田、制定合理的耕作制度、树立用地养地相结合的耕作意识。

值得注意的是，本研究所采用的1980和2010年土壤样点数据差距较大，而且空间分布位置和密度不一，所以通过回归克里格插值模拟出的空间分布图精度差距较大，其产生的误差是否会对研究结果造成一定影响，也是值得探索的方面。