基于地统计学和GIS的辽河上游区域土壤养分空间分异研究

汤 洁，林晓晟，侯克怡，李忠和

(1.吉林大学环境与资源学院，吉林 长春 130012；2.吉林省农业科学院，吉林 长春 130124)

基于地统计学和GIS的辽河上游区域土壤养分空间分异研究

汤 洁1，林晓晟1，侯克怡1，李忠和2

(1.吉林大学环境与资源学院，吉林 长春 130012；2.吉林省农业科学院，吉林 长春 130124)

利用3S技术，结合地统计学方法，对辽河吉林省段土壤总氮、总磷含量的空间变异特征和分布规律进行了研究.结果表明：研究区土壤总氮、总磷含量范围分别为0.33～3.21 g/kg，0.18～0.78 g/kg，变异系数分别为58.08%和41.5%，均为中等程度变异.土壤总氮、总磷的最优半变异函数模型分别为高斯模型和球状模型，空间预测模型以各向异性模型最佳.经Kriging插值分析，总氮、总磷含量为林地>旱田>水田>草地，呈纵向阶梯式变化趋势.土壤总氮密度变化主要受土地利用类型的影响，而总磷密度变化受农业活动的影响较大.区内总氮、总磷总储量分别为3 874.686 t，1 136.49 t，不同土地利用类型的总氮、总磷储量均表现为旱田>林地>水田>草地.

土壤养分；地统计学；半变异函数；Kriging插值

土壤受成土母质、地形、气候等自然因素和人类活动等人为因素的影响，使得土壤成为不均一变化的时空连续体，并具有高度的空间变异性[1].土壤空间异质性存在于多尺度且不同尺度上的主要影响因子有显著差异.研究土壤性质的空间异质性特征，对于了解土壤生态系统的过程具有重要意义.地统计学是基于区域化变量理论基础上的一种空间分析方法[2]，利用地统计学与GIS相结合的方法进行克里格(Kriging)插值，是研究土壤性质空间变异的一种有效方法.

土壤氮素和磷素对土壤理化性质[3]、植被生长[4]以及微生物活动[5]等有重要影响，被公认为最常见的限制植被生长和土壤生产力的关键指标[6].土壤全氮和全磷含量是土壤肥力和土壤质量的重要指标，通常用来表征土壤整体氮、磷含量的水平.其空间分布在水平和垂直方向都具有变异性，这种空间异质性受到气候因素[7]、土壤母质[8]、地形特征[9]、植被类型[10]和土地利用方式[11]等各种自然及人为因素的影响.相对于其他自然因素，土地利用类型在受到气候变化和人类活动影响时，会发生更快速和直接的变化[12].因此，土地利用对土壤全氮和全磷的空间分布具有更加积极的影响.

土壤氮素、磷素随水土流失进入水体可产生富营养化等问题[13-14].辽河吉林省段属于辽河的上游，随着农业和畜禽养殖业的迅猛发展，流域内含氮化肥、农药的施用及持续增加的畜禽粪便污水，致使农业面源污染日益严重，已经对流域生态安全和下游河流水环境产生了严重威胁.本文在前人研究的基础上，应用地统计学的基本原理和方法，利用GS+(9.0)软件和ArcGIS10软件中的地统计学分析模块(geostatistical analysis)，以吉林省内的辽河流域为例，对土壤总氮、总磷的空间变异情况和分布规律进行了研究，以为有效利用土壤、加强农业生产和流域环境管理工作提供依据.

1 研究方法

1.1 研究区概况

辽河流域涉及内蒙古自治区的东南部、吉林省西部和辽宁省的大部分地区.吉林省内的辽河流域位于吉林省西南部，地理位置为东经123°18′～125°36′、北纬42°36′～44°10′，流域面积15 746 km2，约占吉林省总土地面积的 8%.该流域地处松辽平原中部，低山、丘陵和平原兼备，地势由东南向西北缓降，海拔611～120 m，属温带大陆性季风气候区，四季分明；多年平均降水量为545 mm，大部分降水量集中在6—9月，约占全年降水量的80%；多年平均蒸发量为1 020 mm.流域内土壤以暗棕壤和黑土为主，垦殖指数高，耕地比重大，自然肥力高.土地覆盖类型以农地、林地和草地为主.农地以玉米、水稻、大豆为主.

1.2 样品采集与分析

2013年10月，在研究区域内共设置28个采样点，采样采用随机布点并兼顾地形和土地利用差异，采集0～30 cm的表层土壤样品共28个.同时进行GPS定位，记录各个样点植被类型、坡度、高程等样点信息.采样点经ArcInfo投影转换，产生以米为单位的平面坐标系统(投影类型为：Transverse_Mercator投影，中央经线：123°)，结合研究区2012年土地利用遥感解译图，最终生成相应的采样点分布图(见图1).本文的土地利用类型图通过Landsat ETM遥感影像的人机交互解译获得，遥感数据来源于国家科学数据服务平台和美国地质调查局(USGS).

将土壤样品带回实验室风干、研磨，过筛后装入纸袋中.土壤总氮(TN)含量用Foss8400全自动凯氏定氮仪测定(型号：2300Ⅱ)；土壤总磷(TP)含量用分光光度计测定；土壤容重采用体积为100 cm3的环刀取样法测定，共取28个土壤剖面，水田、耕地、林地、草地分别为4，10，9和5个，带回实验室后在(105±2)℃烘箱内烘干10 h后测定容重.

图1 研究区土地利用及采样点空间分布图

1.3 数据分析与处理方法

描述性统计分析采用数理统计软件SPSS17.0进行，统计特征值包括平均值、标准差、变异系数等.半变异函数采用地统计学软件GS+(9.0)进行计算，空间分布图利用ArcGIS10.0中的地统计学模块下的Kriging插值完成.

半变异函数的理论模型可用来分析土壤理化性质空间变异的随机性和结构性，它是地统计学特有的工具和分析的基础[15].变异函数是研究空间变异的关键函数[16]，该函数为：

(1)

式中：γ(h)为变异函数；[z(xi)-z(xi+h)]为间隔为h的2个观测点的实测值，称为步长；N(h)为以h为步长的所有观测点的成对数目.由γ(h)对h作图可得到变异函数图.

实际工作中区域化变量的变异性往往很复杂，它可能在不同的方向上呈现不同的变异性，或者在同一方向上包含着不同尺度的多层次的变异性，为了准确分析区域化变量的主要空间结构，构建理论模型是必需的.本文依据半变异函数的决定系数R2和残差RSS，分别利用两种模型对土壤总氮、总磷的半变异函数进行拟合，得到二者半变异函数的理论模型[17].

半变异函数拟合模型为：

球状模型(spherical model)，其一般公式为：

(2)

式中：C0为块金常数；a为变程；C为拱高；C0+C为基台值.

高斯模型(Gaussian model)，一般公式为：

(3)

式中：C0为块金常数；C为拱高；C0+C为基台值.

土壤总氮、总磷储量的计算借鉴土壤有机碳储量的计算公式[18-19]，即

Ci=di×ρi×οi/100；

(4)

Si=Ai×Ci.

(5)

式中：i为土壤层次；C为土壤养分密度(kg/m2)；d为土壤厚度(cm)；ρ为土壤容重(g/cm3)；ο为土壤养分含量(g/kg)；A为各地类所占面积；S为土壤养分储量.

2 结果与讨论

2.1 土壤养分含量的统计特征

研究区表层土壤理化性质的统计结果(见表1)表明，土壤TN、TP含量最大值和最小值的差值较大，TN为2.88 g/kg，TP为0.33 g/kg，这主要是因为表层土壤受人为干扰较大，与耕作、施肥等农业活动有关，变异系数分别为58.08%，41.5%.根据Nielson和Bouma(1985)的分类系统(弱变异：变异系数≤10%；中等变异：10%≤变异系数≤100%；强变异：变异系数≤100%)，研究区土壤TN、TP含量均为中等变异.

表1 研究区土壤总氮、总磷的统计特征

使用Kriging插值法要求数据集满足正态分布才能正确、合理地对空间采样数据进行变异函数的建模，所以在进行地统计分析之前，需要检验土壤TN、TP的数据集是否满足正态分布，本文采用K-S检验来验证.检验结果表明，土壤TN、TP的K-S检验的P值分别为0.298，0.465.一般来讲，当P值大于0.05时，就说明数据服从正态分布.描述统计分析中的P-P图(见图2)是根据变量的累积比例与指定分布的累积比例之间的关系所绘制的图形，图中的斜线对应着一个均值为0的正态分布.如果图中的散点密切地散布在这条斜线附近，说明随机变量服从正态分布，从而证明样本确实是来自于正态总体.本研究中的28个散点大致分布于斜线附近(见图3)，因此可以认为样本数据服从正态分布.

图2 土壤总氮、总磷含量正态P-P分布图

2.2 土壤理化性质的空间变异性分析

利用地统计学软件GS+(9.0)中对土壤TN、TP含量进行半变异函数模拟，确定半变异函数的空间变异性，最终得到各自的半变异函数模型(见图3、图4).

图3 各向同性半变异函数模型

各向同性半变异函数模拟是在设置了有效距离间隔为102 360.22 m，步长间距为824.01 m后进行的.模拟结果(见图3)表明，TN、TP含量的空间变异函数随步长距离的增加，各向同性表现出一定的差异性，变异函数出现各向异性结构.

图4 各向异性半变异函数模型

土壤TN、TP含量各向异性半变异函数模型及参数值(见表2、图4)，拟合系数(R2)和残差平方和(RSS)的值较小，说明TN、TP的最优拟合模型分别为高斯模型和球状模型.模型的拟合精度较高，能够很好地反映研究区土壤理化性质的空间结构特征.

表2 不同深度土壤总氮、总磷的地统计学参数

块金值反映了区域变化量内部随机性的可能程度，其大小可用来反映区域土壤变化量的随机性.变程和基台值是描述变异函数曲线的两个重要指标，变程反映了土壤变化量影响范围的大小，或该变量自相关范围的大小；而基台值反映了区域土壤理化性质变化幅度的大小，即变异强度[15].土壤TN、TP含量的基台值分别为0.48和0.027，说明研究区TN的变化幅度较大.

块金系数，是块金值与基台值之比，反映土壤性质的空间依赖性，可表明系统变量空间相关性的强度.如果块金系数<25%，说明系统具有强烈的空间相关性；块金系数在25%～75%之间，表明系统具有中等的空间相关性；块金系数>75%，说明系统空间相关性很弱.区内土壤TN、TP含量的块金系数均大于75%，说明系统空间相关性很弱.土壤养分分布是结构性因素(气候、母质、地形、土壤类型)和随机性因素(施肥、耕作措施、种植制度)共同作用的结果.结构性因素可导致土壤变量的空间相关性增强；随机性因素会导致土壤变量的空间相关性减弱，朝均一化方向发展[20].

2.3 土壤总氮空间分布规律

本文的Kriging插值及空间分布图的生成是通过Arcgis10.0软件操作完成的.在选取最优半变异函数模型及其参数的基础上，采用简单克里格(simple Kringing)插值法进行最优插值.插值结果以ASCII数据的格式输出储存，然后生成栅格(grid)格式图形，在Spatial Analyst模块下利用研究区边界提取，最后生成相应的土壤TN、TP含量空间分布图(见图5).

图5 研究区土壤总氮、总磷含量空间分布

空间插值分布图可以更直观地展现土壤总氮含量的变化.从图5可以看出，土壤总氮含量自上游到下游逐渐降低，呈纵向梯度变化，这是因为流域上游为林地，生物固氮性较好，土壤养分含量较高；研究区中游地区为耕地，下游地区为耕地和草地，中游和下游地区受耕作等人为活动影响，使得总氮含量较上游地区要低.土壤总磷含量与总氮含量的空间变化趋势相同，表现为林地和耕地的总磷含量高于草地和水田；与土壤总氮含量不同的是，输入土壤中的磷主要是人为施肥累积所致，使得研究区内旱田的总磷含量较高.土壤中的磷主要以无机磷的形式存在，占土壤总磷含量的60%～80%[21].无机磷含量与土壤类型有关，其主要吸附在土壤黏粒上，导致质地黏重、有机质高的土壤总磷含量较高.研究区内林地土壤为黑土，属于有机质含量较高的黏土，所以总磷含量较高.

图6 不同土地利用类型的土壤总氮、总磷含量

土地利用是自然条件和人为活动的综合反映，土地利用与土壤总氮、总磷含量有密切关系[20].流域内土地利用类型主要有林地、旱田、水田和草地，在ArcGIS10.0软件下的ArcToolbox中，对研究区2012年土地利用遥感解译图和土壤总氮、总磷含量插值图进行叠加分析，提取了不同土地利用类型的总氮、总磷含量，结果见图6.

不同土地利用类型下的土壤TN、TP含量均表现为林地>旱田>水田>草地.各地类的总氮均值为：林地1.66 g/kg，旱田1.32 g/kg，水田1.21 g/kg，草地0.27 g/kg.总磷均值为：林地0.43 g/kg，旱田0.4 g/kg，水田0.39 g/kg，草地0.28 g/kg.总体来讲，东南部林地的总氮含量较高，旱田和草地的总氮含量较低.受人类活动的影响，多年连续耕作可导致土壤养分的下降，导致林地的总氮含量高于耕作地区.总磷的含量为林地和旱田较高，水田和草地较小，这是由人为施肥因素和土壤质地因素共同影响的结果.

2.4 土壤总氮、总磷储量分布规律

本文根据公式(4)计算研究区土壤TN、TP密度，并在ArcGIS10.0软件中进行Kriging插值后，生成相应的空间分布等值线图(见图7).

区内总氮密度的变化范围为0.06～0.56 kg/m2.林地总氮密度高于其他地类，西北部草地的总氮密度最小，这一地区土壤类型为沙地和盐碱地，总氮含量低，导致总氮密度最小.从空间分布图还可以发现，研究区总氮密度较小的地区包括居民区，说明土壤总氮密度变化主要受土地利用类型的影响.总磷密度的变化范围为0.03～0.1 kg/m2，空间变化为区域中心位置密度值较高，并由该中心向区域边界逐渐降低.旱田的总磷密度明显高于其他地类，说明总磷密度变化受人为干扰和农业活动的影响较大.

对研究区TN、TP储量的计算结果(见图8)表明，区内TN、TP总储量分别为3 874.686 t和1 136.49 t.不同地类的总氮储量为：水田151.378 t，旱田2 908.499 t，林地797.627 t，草地17.181 t.不同地类的总磷储量为：水田51.9 t，旱田847.5 t，林地231.174 t，草地5.913 t.不同土地利用类型的总氮、总磷储量均表现为旱田>林地>水田>草地.

图8 不同土地利用类型的土壤总氮、总磷储量

3 结论

(1) 数据统计分析表明，采样点数据服从正态分布.研究区土壤总氮、总磷含量范围分别为0.33～3.21 g/kg，0.18～0.78 g/kg，均为中等程度变异.土壤总氮含量的最优半变异函数模型为高斯模型(Gaussian model)，总磷含量的最优半变异函数模型为球状模型(spherical model).土壤总氮、总磷的基台值较小，空间变异性较弱.土壤总氮、总磷的空间变异性均体现为各向异性结构.

(2) Kriging插值的土壤总氮、总磷含量空间分布表明，土壤总氮含量由东南向西北呈纵向梯度变化.东南部的林地土壤肥沃、养分高，富集了较多的氮；中部为农耕地，土壤肥力相对较丰富，且人为施肥量较高；而西北部以沙壤和盐碱土为主，土壤含氮量很低.土壤总磷含量与总氮含量的空间变化趋势相同，表现为林地和耕地的总磷含量高于草地和水田，这是由于人为施肥因素和土壤质地因素共同影响的结果.土壤总氮、总磷含量均表现为林地>旱田>水田>草地.

(3) 流域总氮密度的变化范围为0.06～0.56 kg/m2，主要受土地利用类型的影响.总磷密度的变化范围为0.03～0.1 kg/m2，受人为干扰和农业活动的影响较大.区内TN、TP总储量分别为3 874.686 t和1 136.49 t，不同土地利用类型的总氮、总磷储量均表现为旱田>林地>水田>草地.

氮、磷元素是植物生长的营养元素，同时也是水体富营养化的重要污染物，为此，土壤总氮、总磷空间分异和分布规律的研究，不仅可为区域土壤环境背景值研究、土地利用规划和管理提供数据支撑，也可为流域水体面源污染的防治提供科学依据.

[1] 秦松，樊燕，刘洪斌，等.地形因子与土壤养分空间分布的相关性研究[J]. 水土保持研究，2007，14(4)：275-279.

[2] 严连香，黄标，邵学新，等.长江三角洲典型地区土壤有效铜和锌的时空变化及其影响因素研究[J]. 土壤通报，2007，38(5)：971-977.

[3] HATI K M，SWARUP A，MISHRA B，et al. Impact of long-term application offertilizer，manure and lime under intensive cropping on physical properties andorganic carbon content of an Alfisol[J]. Geoderma，2008，148：173-179.

[4] QUILCHANO C，MARAN T，PÉREZ-RAMOS I M，et al. Patterns and ecologicalconsequences of abiotic heterogeneity in managed cork oak forests of Southern Spain[J]. Ecological Research，2008，23：127-139.

[5] LIU E，YAN C，MEI X，et al. Long-term effect of chemical fertilizer，straw，andmanure on soil chemical and biological properties in northwest China[J]. Geoderma，2010，158：173-180.

[6] GIESLER R，PETTERSSON T，HGBERG P. Phosphorus limitation in boreal forests：effects of aluminum and iron accumulation in the humus layer[J]. Ecosystems，2002(5)：300-314.

[7] PATIL R H，LAEGDSMAND M，OLESEN J E，et al. Effect of soil warming and rainfallpatterns on soil N cycling in Northern Europe[J]. Agriculture Ecosystems andEnvironment，2010，139：195-205.

[8] LIN J S，SHI X Z，LU X X，et al. Storage and spatial variation of phosphorus inpaddy soils of China[J]. Pedosphere，2009，19：790-798.

[9] REZAEI S A，GILKES R J. The effects of landscape attributes and plant communityon soil chemical properties in rangelands[J]. Geoderma，2005，125：167-176.

[11] ROSS D J，TATE K R，SCOTT N A，et al. Land-use change：effects on soil carbon，nitrogen and phosphorus pools and fluxes in three adjacent ecosystems[J]. Soil Biology and Biochemistry，1999，31：803-813.

[12] OSTWALD M，CHEN D. Land-use change：impacts of climate variation and policiesamong small-scale farmers in the Loess Plateau，China[J]. Land Use Policy，2006，23：361-371.

[13] SCHINDLER D W. Recent advances in the understanding and management ofeutropication[J]. Limnology and Oceanography，2006(5)：356-363.

[14] CHEN M，CHEN J，SUN F. Agricultural phosphorus flow and its environmentalimpacts in China[J]. Science of the Total Environment，2008，405：140-152.

[15] 刘爱丽，王培法，丁园圆. 地统计学概论[M]. 北京：科学出版社，2012：5.

[16] 徐国策，李占斌，李鹏，等. 丹江中游典型小流域土壤总氮的空间分布[J].地理学报，2012，67(11)：1547-1555.

[17] 王政权.地统计学及其在生态学中的应用[M]. 北京：科学出版社，1999：35-149.

[18] SUN WENYI，GUO SHENGLI. The spatial distribution of soil organic carbon and its influencing factors in hilly region of the Loess Plateau[J]. Acta Ecologica Sinica，2011，31(6)：1604-1616.

[19] XU WENQIANG，CHEN XI，LUO GEPING，et al. Soil organic carbon storage and its spatial distribution characteristic in the Sangong River watershed of arid region[J]. Journal of Natural Resources，2009，24(10)：1740-1747.

[20] WANG ZHENGQUA，WANG QINGCHENG. The spatial heterogeneity of soil physical properties in forests[J]. Acta Ecologica Sinica，2000，20(6)：945-950.

[21] 沈善敏. 中国土壤肥力[M ].北京：中国农业出版社，1997.

(责任编辑：方 林)

Study on spatial variation and distribution characteristics of soilnutrient in Liao River Basin based on geostatistical methods and GIS

TANG Jie1，LIN Xiao-sheng1，HOU Ke-yi1，LI Zhong-he2

(1.College of Environment and Resource，Jilin University，Changchun 130012，China；2.Jilin Academy of Agricultural Sciences，Changchun 130124，China)

Based on 3S technology，using geostatistical methods to study the spatial variation characteristics and distribution law of soil total nitrogen and phosphorus content in Liao River Basin of Jilin Province. The results indicated that the total nitrogen and phosphorus content was 0.33～3.21，0.18～0.78 g/kg，and its coefficient of variation were 58.08% and 41.5% which belongs to moderate variation. The best fitted model in different soil type were Gaussian models in TN and Spherical models in TP. A proper model for the Spatial prediction of TN and TP was the Anisotropy model. The Kriging interpolation showed that the spatial distribution of soil total nitrogen and phosphorus content exhibited a trend of longitudinal gradient change，TN and TP storage of different land use types is forest>crop field>paddy field>grassland. TN and TP total reserves are divided into 3 874.686，1 136.49 t，different land use types of total nitrogen，total phosphorus reserves are characterized by dry field>forest land>paddy field>grassland.

soil nutrient；geostatistical method；semivariagram model；Kriging interpolation

1000-1832(2014)04-0139-08

10.11672/dbsdzk2014-04-025

2014-05-22

国家水体污染控制与治理科技重大专项课题(2012ZX07202-009).

汤洁(1957—)，女，教授，博士研究生导师，主要从事生态环境系统理论与技术研究.

S 153.6 [学科代码] 210·5015

A