复杂地形区土壤有机质空间变异性分析及制图

张 欢,高小红,2,3,4

(1.青海师范大学 地理科学学院,西宁 810000;2.青海省自然地理与环境过程重点实验室,西宁 810008;3.青藏高原地表过程与生态保育教育部重点实验室,西宁 810008;4.青海省人民政府—北京师范大学 高原科学与可持续发展研究院,西宁 810008)

有机质在土壤中扮演着重要角色，是植物营养中的重要来源，它是促进植物生长发育、土壤保肥性和缓冲性、改善土壤物理性质的主要因素，目前土壤有机质在改善农业、管理耕地等方面具有重要意义，是土壤学、地球化学、农学研究的重点内容之一[1-2]。

传统的土壤预测制图是一个手工制图的过程，需要对地形图等进行目视解译以及对土壤类型斑块进行手工勾绘，除此之外还需要野外核查才能编辑成数据图件[3]，需要耗费大量的人力物力财力。数字土壤制图是以土壤与景观理论为基础，借助数学、空间分析方法等技术手段的土壤制图与调查方法[4]，制图方法主要包括基于地统计的数字土壤制图和基于专家知识的数字土壤制图[3]。研究表明：地统计学方法在数字土壤制图中不仅考虑空间上的邻近点相关性，也考虑了该点的其他地理要素相关性，可以较好的用于数字土壤制图[4]。许多学者陆续地使用GIS与地统计方法探讨土壤属性的空间特征[2]。赵永存等[5]发现回归克里格的预测效果最好，泛克里格次之，而多元线性回归方法最差，证明地统计学插值方法用于空间预测相对较好；席江勇等[6]在对土壤有机质的插值研究中利用了确定性内插和地统计插值进行分析，使用确定性插值产生了牛眼，而地统计插值消除了这种现象，结果连续且光滑，能更好的反映其在空间上的变化；综上所述，在研究土壤养分时考虑地统计插值方法能够取得较好的预测精度。

目前的一些研究表明，将土地利用类型结合到地统计学插值中能够获得较好的预测精度。文雯等[7]研究发现在黄土丘陵地区，利用土地利用修正克里格插值方法精度比其他插值方法精度高，能够提高土壤属性空间制图的精度；顾成军等[8]利用有限的采样点将土地利用结合到克里格插值中，发现在该研究区土地利用对有机碳的空间分异有很好的解释能力且不需要获取其他因子；赵建华[9]在研究兰州市土壤全氮和有机质的空间变异时发现其空间分布特征与研究区的各种土地利用类型有重要的联系。复杂地形区不同于地形比较简单的区域可以使用地统计分析的简单克里格进行插值预测，需要考虑有关的辅助因子来进行土壤属性空间预测研究，以提高预测精度[10]。协同克里格法基于目标变量与辅助变量，是一种考虑辅助变量的有效方法[11]，而在实际情况中，空间数据是随着地理位置的变化而变化的,地理加权回归模型正是考虑了空间数据的位置信息[12]，近年来在土壤属性空间制图中得到较好的应用。地理加权回归克里格是地理加权回归与克里格法的有机结合，该方法不仅考虑了空间位置局部信息，还考虑了代表随机性的残差数据，能够反映出更加真实的土壤属性空间变异[13]。

湟水流域位于青藏高原与黄土高原过渡带，其地形起伏大、地貌类型多样，主要以农业为主。目前对该研究区进行土壤属性空间制图的研究还比较少见。代子俊[14]分别利用1985年和2015年的土壤采样数据，对湟水流域土壤全氮进行克里格插值制图，但由于其土壤采样点较少且在流域西北区域基本没有采样点，这在一定程度上将影响空间制图的精度；同时虽然其在土壤采样时考虑了土地利用类型因素，但在空间插值时并未将其结合到土壤全氮制图中。基于以上的考虑，本文对湟水流域土壤有机质含量进行空间地统计学制图，旨在研究所使用方法在复杂地形区的适用性以及比较协同克里格方法、地理加权回归克里格法与结合土地利用类型后的协同克里格法在土壤有机质制图方面的精度，为湟水流域提供较为准确的土壤有机质数据支持。

1 数据与方法

1.1 研究区概况

湟水流域位于青藏高原与黄土高原的交接地带，位于36°02′—37°28′N，100°42′—103°04′E，地形在南北方向，具有“三山两谷”独特的构造景观，三山从南向北为拉脊山、达坂山、祁连山，两谷为介于达坂山与拉脊山之间的湟水干流河盆谷地、达坂山与祁连山之间的大通河河谷。青海省境内湟水流域面积约为16 120 km2，海拔高度在1 664～4 882 m，地势呈现东南低西北高的趋势,具有多样复杂的地形,黄土地貌在该区最为典型，区域内土地利用类型多样，在浅山区主要以农田为主,川水区以农田和城乡工矿居住用地为主，海拔3 000 m以上主要是草地占多,湟水流域土壤类型主要是灰钙土、黑钙土和栗钙土为主。

湟水流域年平均气温在2.5～7.8℃，该区冬季气温低、夏季气温适宜，最高气温可达30℃左右，而最低气温可达-28℃左右，日照长、可吸收辐射大，且从东到西逐渐增大；区域平均降水量486mm，主要集中在6—9月份，10月到次年2月降水量最少，只占全年降水量的10%，而蒸发量最大可达1 000mm，且随海拔高度的下降蒸发量是逐渐增加的，表现出西北向东南蒸发量增加的趋势[15]。该区是青海省的主要农业耕作区，也是青海省人口最密集、经济最发达的地区。

1.2 数据来源与预处理

本文使用的数据主要有数字高程模型(Digital elevation model，DEM)(图1A)、土壤采样数据253个(图1B)、土地利用数据(图1C)、土壤类型数据(图1D)、湟水流域边界数据。其中土壤类型矢量数据通过1∶100万青海省土壤类型图矢量数字化获取；土壤有机质来源于2016年10—11月份野外采集数据，采样土层厚度为0—20 cm，土壤有机质(Soil organic matter，SOM)的测试分析请参见李冠稳[16]。将采样数据的野外经纬度坐标(GPS获取，WGS-84)转换为UTM投影坐标，DEM来源于USGS网站(http:∥www.usgs.gov/)，分辨率为30 m，对DEM进行填洼、拼接、裁剪等预处理，最终得到湟水流域DEM数据，利用ArcGIS软件计算出坡度(图1E)、坡向(图1F)以及高程；研究区土地利用类型数据根据2016年Landsat8/OLI影像通过随机森林分类方法获取，总体精度脑山区、浅山区、川水区分别为91.33%，92.09%，87.85%，见马慧娟论文[17]，对该土地利用类型数据进行合并后包括7类，即水浇地、旱地、草地、林地、水域、城乡居住建设用地、未利用土地，并利用土壤有机质采样点提取相应的土地利用类型数据。

图1 研究区数据

1.3 研究方法

1.3.1 描述性统计与单方差分析 首先对土壤有机质数据以及土地利用数据利用SPSS 17.0进行描述性统计分析；统计出最大值、最小值以及峰度、偏度等数据值，并通过K-S检验数据是否呈正态分布，由于有机质数据不呈正态分布，对其进行对数转换或者剔除异常值，直到数据呈现正态分布。使用单因素方差分析比较主要土地利用类型的土壤有机质是否具有显著差异[18]，如果p>0.05，则土地利用类型来自于同一个正态总体，不用消除由地类差异引起的土壤有机质差异，确保土地利用类型的有机质可以看做是同一个区域化变量在空间上的分布[19]。

1.3.2 结合土地利用类型的协同克里格法 土地利用数据属于分类变量或定性变量，不具备数学上的大小关系，因此不能直接用于地统计学分析。根据前人研究[19]，有两种方法可将分类变量转换为定量变量，第一种方法为空间分区，该方法只适合大面积连续的区域；第2种方法为对不同土地利用类型的土壤属性值进行均值、中值中心化，并利用其残差进行插值，本文研究区属于复杂地形区，地形起伏大、土地类型斑块破碎，不像平原地区具有大面积连续的单一类型区域，所以采用第2种方法将定性变量转为定量变量，公式如下：

R(xi)=Z(xi)-Q(xi)

(1)

Z(i)=R(i)+Q(i)

(2)

式中：Z(xi)为有机质实测值；Q(xi)为拆分趋势项；R(xi)为拆分残差项；Z(i)为未采样点的土壤有机质预测值。

1.3.3 地理加权回归克里格法 地理加权回归克里格是在地理加权回归的基础上进行预测，所以本文首先对地理加权回归模型进行拟合：

(3)

式中：YGWRK(ui，vi)为因变量在点(ui，vi)处的估计值；xij为第i个解释变量在(ui，vi)处的值；β0(ui，vi)为样本点在(ui，vi)的特征弹性系数；βp(ui，vi)为第i个样点上的第p个回归系数。拟合的关键参数主要是核函数与带宽，经过测试，本文使用的核函数为自适应高斯核函数，带宽以AICc赤池信息量准则为判断标准，值越小说明拟合后信息的不确定性越小[13]。地理加权回归的模型拟合参数见表1。

表1 地理加权回归模型拟合参数

1.3.4 半变异函数 半变异函数又称半方差函数，是关于数据点间的变异性距离的函数，也是描述区域化变量结构性、随机性的基本手段。估算方程式如下：

(4)

式中：h是分隔距离；N(h)是(xi+h,xi)之间计算样本变异函数值的样本对数(如Z(xi+h)和Z(xi)为计算中的第一对数据)[12]。

1.3.5 空间插值 本文创建训练集后首先对数据进行正态分布检验，均值与中值残差、GWRK残差、原值均符合正态分布，可以用于空间插值。

分别采用均值、中值消除趋势项，并用得到的插值图加上不同土地利用信息的趋势项，最终得到结合土地利用类型的土壤有机质预测图。同时，进行地理加权回归模型拟合，得到的残差项进行空间插值，最后与协同克里格的插值结果比较，制作土壤有机质分布图；本文称4种方法为均值修正协同克里格(Land Use Mean Modification Cokriging，MMCOK _LU)、中值修正协同克里格(Land Use Median Correction Cokriging，MCCOK_LU)、协同克里格(Cokriging，COK)、地理加权回归克里格(Geographically Weighted Regression Kriging，GWRK)。

1.3.6 模型精度评价 本文通过内部验证与外部验证进行精度评价，内部验证：标准均方根(Root-Mean-Square Standardized，RMSS)预测误差和1最接近、标准平均值误差(Mean Standardized Error，MSE)与0最接近、均方根(Root-Mean-Square，RMS)预测误差达到最小、平均标准误差(Average Standard Error，ASE)和均方根(Root-Mean-Square，RMS)预测误差最相近。外部验证：平均绝对值误差(Mean Absolute Error，MAE)越小模型精度越高、均方根误差(Root Mean Squares Error，RMSE)越小模型越稳定[20]，准确度(Accuracy，AC)可评价预测的准确性，越接近于1说明预测越准确，取值范围为0～1[7]。r为预测集与验证集的相关系数，RI为比较GWRK相对于COK的提高度。

(5)

(6)

(7)

(8)

(9)

(10)

2 结果与分析

2.1 土壤有机质描述性统计分析

从土壤有机质描述性统计分析(表2)可以看出采样点的土壤有机质含量在7～150 g/kg，标准差为24.2，平均值为34.16 g/kg。

表2 土壤样本有机质描述性统计

CV为变异系数，根据CV<10%为弱的变异性,10%～100%属于中等变异,CV>100%是强变异性，可知总体样本CV为90.90%，属于中等变异强度；偏度为2.047，属于右偏态，峰度4.761，不是特别陡峭；在K-S检验中，p的值等于0，说明其不呈正态分布，不能直接用于空间插值研究，所以需要进行对数转换，转换后其p值为0.055，大于0.05，可以用于地统计分析；除此之外根据表3对各土地利用类型的统计，土壤有机质都处于中等变异水平，林地的变异程度大于其他3种土地利用类型。

表3 样本所对应的土地利用类型的描述性统计

不同土地利用类型的土壤有机质单因素方差分析(表4)表明，由于p>0.05,四组土地利用的SOM并未存在很大的差异，故SOM值是属于同一总体的，使水浇地、林地、旱地、草地的SOM值可以视为同一区域化变量，不用在插值前消除由这些类所引起的SOM差异，故其残差值也无明显的不同，都可以看成在空间上分布的同一个区域化变量[19]。

表4 不同土地利用类型SOM单因素方差分析

2.2 空间变异特征分析

首先利用土地利用类型，GWR模型得到残差项，然后利用ArcGIS 10.2地统计模块建立训练子集209个，验证子集31个，最后利用ArcGIS 10.2地统计模块的协同克里格插值，插值过程中结合湟水流域DEM因子，拟合最优的半变异函数模型(表5)。本文4种方法均为指数模型。

表5 SOM空间变异的理论模型与参数

结构、随机性变异是空间变异的两大模块，而C0，C0+C，C0/(C0+C)常常用来表示半变异函数建模的空间变异程度，其中C0为块金值，表示可能的随机性，这种随机因素可能为土地的耕作、施肥、管理措施等；C结构方差由自然因素，即地形、地质、土壤母质、土壤类型等引起的土壤性状的结构变异；C0+C为基台值，空间总的变异一般由基台值表示；C0/(C0+C)为块金与基台值之比，称为块金系数，其大小表示空间的变异程度，值小则表示空间变异程度大多是结构性因素引起的，反之值大则为随机性因素引起[2]。

研究区MMCOK_LU，MCCOK_LU，COK模型的块金系数分别为17.8%，19.9%，19.3%，属于强的空间相关性，GWRK的块金系数为47.5%，属于中等的空间相关性；前3个模型表明土壤有机质空间变异主要由结构性因素引起，而最后一个模型表明土壤有机质空间变异主要由结构性和随机性因素引起，原因可能是在考虑土地利用类型后，空间细节信息更丰富，人为的影响更突出[C0/(C0+C)>75%:弱空间相关；C0/(C0+C)<25%:强空间相关；之间属于中等空间相关]。变程表示变量在空间上的自相关性，MMCOK_LU和MCCOK_LU的变程均为6 899 m，COK为6 790 m，GWRK为5 977 m，表明在这个范围内空间自相关是存在的，均选择指数模型作为最佳模型。

2.3 模型精度分析

由表6可知RMSS标准均方根预测误差MMCOK_LU最接近于1，4种方法的MSE相差不大，COK，GWRK仅比MMCOK_LU与MCCOK_LU小0.001，COK的ASE与RMS接近程度为0.003 9，大于MMCOK_LU与MCCOK_LU；MAE平均绝对误差MCCOK _LU最小，GWRK最大，说明MCCOK _LU的模型精度最高，其次分别为MMCOK_ LU，COK；RMSE均方根误差GWRK最小，MCCOK_LU最大，说明GWRK比COK，MMCOK_LU，MCCOK_LU稳定性高；而对于模型的准确度，由表6中的AC可知预测准确度由大到小依次为MCCOK_LU>GWRK>MMCOK_LU>COK；与COK相比，GWRK相对提高3.3%。综上所述，可以发现结合土地利用中值残差的MCCOK_LU模型精度最高。

表6 模型预测精度评价

2.4 有机质空间插值及结果与分析

图2中可以看出4种方法的SOM预测值分布范围相差不大，COK(协同克里格)与GWRK(地理加权回归克里格)的土壤有机质空间分布较为连续，而MMCOK_LU(均值修正协同克里格)、MCCOK_LU(中值修正协同克里格)分布比较破碎，但是反应的地类细节特征以及信息的丰富程度较好；MCCOK_LU高值区相较于MMCOK_LU少，GWRK的低值区最少，四者的分布中间值最多。

4种方法所产生的高值均分别集中在大通县西北部、海晏县南部和湟源县西北部等地区，互助县中西部以及湟中县北部也有零星高值出现；其中大通县西北部、大通县与海晏县交界处的高值可以结合DEM图、土地利用图以及土壤类型图(图1)进行分析，该区海拔相对较高，薛晓娟等[21]的研究发现，海拔3 400 m以上土壤有机质含量开始急剧上升；且该区受人为因素影响小，多为林地和草地(动植物残体腐烂分解为腐殖质),所以该地区主要为SOM丰富的地区，且该区主要发育高山草甸土，通过青海土壤[22]高山草甸土有机质含量在不同的区域差异较大但是普遍含量较高，这是因为高山草甸土的土壤形成过程慢，表层粗的有机质易于积累，所以有机质含量较高，该区土壤有机质含量大约在38.8 g/kg左右，预测结果与该值相符；在海晏县南部地区和湟源县西北部地区、互助县西部、西纳川、石崖庄主要分布栗钙土，部分山地草甸土与高山草甸土，该区栗钙土的土壤有机质含量较高，分解微生物的能力较强[14]，根据青海土壤[22]该区水热条件较好，土壤较为肥沃，土壤有机质在40～55 g/kg左右，预测结果在50 g/kg左右，两者大致相符，且土地利用类型主要是草地、林地、部分旱地和水浇地，草地、林地有机质含量丰富，而旱地作物生长周期比较短，施肥较多，所以SOM值会偏高；民和县西部、西南部、乐都县东南部主要土地利用类型为耕地与草地，栗钙土分布广，根据代子俊[14]的研究，发现该区农作物产出率有所提高，离不开当地农民对农田的精细管理和施肥，所以增加了土壤有机质的含量，该区SOM含量大致在55 g/kg左右。有机质低值区主要出现在大通回族土族自治县西南部、民和县东南部、乐都县南部地区；民和县东南部，海拔最低，坡度相对较大且处于流域下游，如果遇有雨水，则对土壤的冲刷能力较大，而该区的栗钙土热量有余，水分不足，水热矛盾突出，SOM含量低，不易积累，加之水土流失强[22]，所以会呈现出较低值；大通回族土族自治县西南部主要分布黑钙土与部分栗钙土，土地利用类型主要是坡旱地，海拔在2 500～3 000 m左右，土温较高，土壤的水分较适中，土地耕作后透气性加快，使SOM矿化加快，除此之外该地耕种时间长，耕层经常烧灰，土壤颜色变浅，表层的SOM含量有明显下降趋势[22]，所以该地有机质含量不是很高；其他地区有机质含量基本处于预测值中间水平，即27 g/kg左右。

水浇地区域有机质含量受水量的影响较大，丰水期水储量大时会加速土壤的分解，会产生更多的有机质，此时有机质比较丰富，但是研究区处于干旱半干旱区，降水变化大，水量不是很充足，且日照大蒸发快，所以水浇地的SOM含量随气候的不同而发生变化；结合图2与图1C发现草地的有机质含量在大部分地区高于林地，通过文献[23]发现主要原因是该地温性草原分布广，由放牧等人为因素影响较大，动物粪便等都会对SOM产生影响。

图2 土壤有机质含量分布

3 讨论与结论

本研究表明复杂地形区的湟水流域2016年SOM平均含量为34.16 g/kg，属于中等水平。变异系数为90.90%，存在中等空间变异性，说明作为重要的农业区域，其土壤有机质同时受土壤母质、土壤类型以及土地利用管理方式、施肥等人类活动的共同影响。

对于地形复杂区，利用土地类型信息辅助研究土壤属性空间分布制图的文献还较少，本文除了考虑高程外，还将土地利用信息对SOM的影响结合到有机质空间变异性分析与空间制图中，研究表明，利用较易获取的土地利用类型研究SOM的空间变异性，能有效的降低空间异常值的影响，提升模型的预测精度，其插值精度均较高。所使用的4种方法中结合土地利用类型的中值修正协同克里格精度最高，其次为地理加权回归克里格、均值修正协同克里格、协同克里格,模型的准确度从高到低分别为0.923，0.909，0.905，0.883，GWRK相对COK提高3.3%。本研究中的土地利用中值修正协同克里格制图精度要好于均值修正，这与吴子豪等[19]研究土壤有机碳密度的空间异质性结果基本一致；文雯等[7]在研究黄土丘陵小流域时，中值修正克里格相对于未修正克里格准确度提升0.042 4，本文中的中值修正克里格相对于未修正克里格准确度提升0.040，稍低于前者的研究，原因可能是研究区相对于黄土丘陵区地形更加复杂，且研究中只考虑了4种土地利用类型，差异相对较小；本文利用协同克里格与辅助变量相结合，提高了模型的预测精度，与前人的研究结果一致[24]。

地理加权回归克里格、协同克里格均能有效地对复杂地形区的湟水流域土壤有机质进行制图，但前者精度高于后者。协同克里格与土地利用类型结合后，中值协同克里格与均值协同克里格预测精度均稍高于普通协同克里格。因此，考虑土地利用类型信息的地统计学在复杂地形区土壤有机质制图是非常有效的。研究区的土地利用类型为7类，对SOM的贡献程度分别为：草地为0.121，水浇地为0.063，旱地为0.035，林地为参照(将土地利用类型从名义变量转换为定量变量，做出土地利用类型的回归方程，得出土地利用方式的贡献程度)，贡献程度最高的草地，采样点数仅占总样本数的17%，小于一些学者[7]的研究，采样点空间分布不均匀以及个别样本较少是本文的不足。除此之外本文未将定性变量考虑进地理加权回归模型，目前的研究未能很好解决的这个问题，今后将对其进一步的研究。