1 引言(Introduction)

在地球陆地生态系统中，土壤储存的氮是整个生态系统氮的主体，而且最为活跃.土壤中氮素含量的大小和分布与陆地生态系统的养分循环、土地生产力以及生态系统的健康具有密切的相关性(程先富等，2007；徐国策等，2012；罗由林等，2015).目前在农业生产中氮素研究的主导思想已突破过去单一提高粮食单产的局限，向高产高效和生态环境保护协调发展的观念转变(朱兆良，2008；李林森等，2015).土壤氮素含量过高或过低都会引起作物生长、水体环境污染等一系列农业问题，特别是近些年，由于氮肥施用过量或施用方式不合理造成的生态环境问题日益受到人类重视(陈肖等，2007；闫建梅等，2014).过多的土壤氮已成为地表水和地下水污染的“源”，是农业非点源氮污染治理的关键，摸清区域内土壤氮素的“家底”及其空间规律对氮素的合理投入意义重大(朱兆良，2008；司涵等，2014；杨小林等，2013).土壤全氮(TN)是衡量土壤氮素供应状况的重要指标，碱解氮(AN)作为一种速效态氮，与农业生产有着直接关系.已有大量学者对不同区域土壤TN和AN的空间分异规律进行了研究(吴德勇等，2014；胡玉福等，2015；裴宏伟等，2015；Zhang et al.，2011；Zhang et al. ，2016)，但其研究多以行政区划或单一地块为研究单元，忽略了土壤中各养分元素累积形成的自然过程.流域是一种集社会经济和自然资源于一体的复合型系统，以流域为单元进行土壤TN和AN的空间分异规律研究能够更好地揭示社会经济及自然因素与土壤TN和AN形成及其空间分布特征的内在关系.中国东北黑土区以土壤肥沃、作物产量高著称，是我国重要的商品粮、畜牧基地.经过近几十年的高强度开发，黑土区土壤侵蚀现象日益严重，农业污染、土壤生态环境恶化等问题逐渐凸显(雷国平等，2009；王充等，2014；彭显龙等，2014).因此，以提升粮食产能和保护生态环境为共同目标，对我国东北黑土区小流域土壤中TN和AN含量进行深入研究在区域环境安全评价和改进施肥管理制度方面具有迫切的现实意义.

本文综合应用了地统计学和“3S”技术，对黑土区海沟河小流域土壤中TN和AN的空间分异特征及分布格局进行探讨，并对影响土壤TN和AN空间布局的地形、环境、土地利用等外部因素进行定性和定量的分析，尝试研究了环境等外部因素对土壤TN和AN的影响程度及其响应的区域范围，可为有针对性的土壤保护与培育及划定农业生产中环境保护的重点区域提供依据，并进一步促进黑土区粮食产能提升与生态环境保护的协同发展.

2 材料与方法(Materals and methods) 2.1 研究区概况

海沟河小流域(地理坐标为126°55′45″~127°10′05″E，45°34′18″~45°40′50″N)位于哈尔滨市东南，以料甸满族乡为主，兼有红星乡部分村屯，总面积为119.76 km²(图 1)，境内海沟河是松花江的二级支流，自西向东汇入阿什河并在哈尔滨市东南注入松花江.小流域内地形呈东高西低之势，海拔在54~366 m之间，东部为低山丘陵区，向西逐步过渡为平原耕作区.流域内土壤类型以黑土为主，兼有部分草甸土、白浆土和暗棕壤分布，属典型的黑土区.气候类型为中温带大陆性季风气候，冬季干旱漫长，夏季温热湿润，多年平均气温3.9℃，年降水量600~800 mm，多集中在每年的4—10月，占全年降雨量的90%.海沟河小流域土地资源开发历史较早，农业生产活动密集，西部平原以农业生产为主，耕地面积占总面积的70.97%，东部低山丘陵区多林业分布，森林覆盖面积占总面积的24%.研究区地处大城市边缘区，受城市社会经济辐射作用强，流域内交通发达，社会经济发展迅速，但同时也给流域内带来一定的负面影响：土壤污染、土壤肥力下降、水源污染等农业生态环境问题凸显.

2.2 数据采集与测试分析

基础数据：本研究采用Google Earth 2014年高分辨率遥感影像，通过选取地面控制点进行几何校正、配准，在ArcGIS10.0平台中采取人工目视判读法进行影像解译，获取海沟河小流域土地利用数据(图 2)，通过野外定点采样调研，遥感影像解译准确率达90%以上，满足本研究需要；地形数据主要从国际科学数据共享平台提供的30 m×30 m的DEM数据中获取(图 2)；通过对1：1万土壤类型图进行矢量化获取土壤类型数据，其他数据来源于相关政府部门.

采样过程：以研究区地形图、土地利用及土壤类型图等数据为辅助信息，综合考虑地形、土地利用方式及距离居民点、交通主干道、河流水系的距离等要素，在室内进行网格布点.为避免作物生长、耕作施肥等对本研究的影响，在农作物收获后，于2014年10月进行野外采样，主要集中于流域西部耕地区.采样过程中，通过手持GPS进行空间定位，与室内预布点位尽可能一致，每个样点按照不同方位等距离采集5次，采集土壤表层0~20 cm耕层土壤，用“四分法”取混合均匀的样品约1 kg，留足实验用量，并分别记录采样点坡位及周边环境等相关信息，共采集120个土壤样品(图 1).

测试分析：土壤TN含量分析采用半微量开氏法，土壤AN含量分析采用碱解扩散法(鲍士旦，2000).

2.3 数据处理

采用格拉布斯(Grubbs)法进行异常值识别，TN和AN分别有119个和117个样本点数据符合要求；通过随机抽取样本点，将样本点分为训练集(TN共89个样本点，AN共87个样本点)和验证集(TN和AN各30个样本点)，其中，训练集参与本文的地统计分析、空间插值及影响因子分析，验证集只参与对插值结果进行的独立验证，不参与后续分析过程.

2.4 研究方法 2.4.1 地统计模型

地统计学中的半方差函数是土壤属性空间变异特征研究中最有效的方法，是进行土壤属性空间插值的基础(史舟和李艳，2007；邓宝山等，2015).本研究在GS+9.0中主要对土壤TN含量、AN含量及其回归残差进行半方差函数拟合，其一般过程如下：

(1)

式中，γ(h)为变异函数；N(h)是分割距离为h时的样本点对总数；Z(x_i)为Z(x)在空间位置x_i处的实测值；Z(x_i+h)是Z(x)在x_i处距离偏离h的实测值.

2.4.2 空间插值方法

对海沟河小流域土壤TN和AN的分布进行空间插值，插值方法很大程度决定了插值精度的高低以及空间分布特征分析的准确性和结果的可靠性.回归协同克里格(Regression Co-Kriging)是一种考虑外部环境要素的插值方法，通过在SPSS中获取土壤TN和AN与指标因子的最佳线性回归关系，分别得到代表确定性部分的趋势项和随机部分的残差项，然后通过协同克里格法对残差项进行插值，最后将趋势项和残差插值结果相加得到回归协同克里格插值结果，本文采用回归协同克里格插值法来获取土壤TN和AN的空间分布情况，其过程如下(张素梅等，2010；杨顺华等，2015):

(2)

式中，f(x)为SOM在x处的模拟值；m(x)为逐步回归过程所得到的趋势项；ε(x)为协同克里格对残差项的插值结果.

对回归协同克里格插值结果通过预先布设的验证集样点采用独立验证法进行插值精度的检验，主要采取平均误差(ME)和均方根误差(RMSE)作为验证指标，已有研究对其相关计算过程及代表意义进行了描述(张素梅等，2010；杨顺华等，2015)，此处不再赘述.

2.4.3 主控因素分析

为尽可能全面地探究影响研究区土壤TN和AN空间分布的主控因素.一般认为，地形因子(高程、坡度等)能够从土壤自然形成的角度对土壤养分的空间分布进行解释，居民点分布、河流水系在一定程度上能够反映人类活动对土壤改造与演变过程的影响程度和范围，土地利用因子(土地利用类型、坡位)则反映了人类在农业生产过程中对土壤养分含量的直接影响过程.本研究对土壤TN和AN含量与基本地形因子进行相关性分析，提取与氮素含量显著性相关的地形因子；为克服样点所受外部性的影响，以空间距离接近的5个样点为1组，计算采样数据的均值及样点到河流水系、居民点间的平均欧氏距离，最后与土壤TN和AN含量做曲线估计，寻找土壤TN和AN含量在空间上发生变化的拐点及变化范围；对不同土地利用和坡位下样点的土壤TN和AN含量的差异性分析，确定其与TN和AN含量之间的定量关系.

2.4.4 软件平台及数据分析

在SPSS19.0软件中对样本数据进行描述性统计分析及正态分布性检验，用K-S检验法进行非参数检验(α=0.05)，绘制频率分布图，并进行氮素、回归残差与指标因素间的相关性分析、回归分析及方差分析等；半方差函数及理论模型的拟合等计算过程均采用地统计学软件GS+9.0进行；空间插值、欧氏距离计算及相关绘图在ArcGIS10.0和Origin9.0中完成.

3 结果与分析(Results and analysis) 3.1 土壤中氮素的描述性统计特征

对海沟河小流域土壤TN和AN有效样点数据进行经典统计学分析，得到海沟河小流域土壤TN和AN有效样点集的描述性统计特征值(表 1).

从表 1可以看出，土壤TN与AN含量的范围分别是0.96~5.50 g·kg^-1和34.02~207.09 mg·kg^-1，均值分别为2.60 g·kg^-1和114.85 mg·kg^-1，参考全国第二次土壤普查的氮素含量分级标准及黑土区相关研究(黑龙江土地管理局，1992；韩秉进等，2007；彭显龙等，2014)，整体上研究区土壤中TN含量处于较高水平，AN含量为中等水平，变异系数分别为83.46%和25.77%，结合一般的CV评估标准(当CV < 0.1时，为弱变异程度，当0.1 < CV < 1时，为中等变异程度，CV > 1时为强变异程度)，研究区土壤TN和AN含量的空间变异程度均为中等变异程度.K-S检验表明TN含量需要进行对数转换后才符合正态分布，AN含量符合正态分布，最后可以进行样点数据的地统计与空间插值分析.

3.2 土壤中氮素空间变异特征分析

通过GS9.0+软件对研究区土壤TN和AN进行地统计分析，得到研究区土壤TN和AN各向同性条件下半方差函数模型(图 3)及相关拟合参数(表 2).结果表明土壤TN和AN的半方差理论模型均为高斯模型(Gaussian)时，拟合残差RSS分别为1.154×10^-4和5.257×10^-4，拟合程度较高，能够有效反映研究区土壤TN和AN的空间变异特征.

块金值(C₀)可以揭示区域变量在小于抽样尺度上发生突然变异的程度(史舟和李艳，2007；李启权等，2010).海沟河小流域土壤TN和AN含量的块金值分别为0.001和0.002，表明该区域内存在某种不可忽略的生态过程对土壤TN和AN的空间分布产生影响.基台值(C₀+C)通常用来表征系统内总的变异程度，块金值与基台值之比(C/(C₀+C)则表征了系统内随机部分引起的变异占系统总变异的比例(刘璐等，2010；苏松锦等，2012)，研究区土壤TN和AN的C/(C₀+C)值分别为0.03和0.114，表明海沟河小流域土壤TN空间变异更容易受到成土母质、地形等结构性要素影响，而受人类活动(耕作制度、施肥措施等)随机性要素影响小(连纲等，2008；Zhang et al.，2011；胡玉福等，2015)；同时，土壤AN的C/(C₀+C)值明显大于TN，在该研究区内，人类的农业生产活动对土壤AN空间变异的影响作用更大，在一定程度上增加了研究区土壤AN空间变异及分布研究中由人类活动所产生的不确定性.变程是半方差达到基台值的样本间距，也称空间最大相关有效距离，通常采样距离超过该范围，则没有空间相关性(史舟和李艳，2007)，海沟河小流域土壤TN和AN的变程分别为900 m和1282 m，表明土壤TN和AN含量均存在较强的空间自相关性，该结果可以作为以后在相似尺度上进行土壤TN和AN空间分异规律研究时野外采样间隔距离的参考值.

研究区地形地貌复杂，基于以上分析进行土壤TN和AN的各向异性分析，绘制0°、45°、90°和135°方向的半方差函数图(图 4).结果表明：土壤TN含量在间隔距离小于3 km时，南-北(N-S)向变异程度较大，东-西(E-W)向变异较小，之后趋于同性，整体上，西北-东南(WN-ES)、西南-东北(WS-EN)两个方向随着间隔距离的变化可近似为各向同性，但当间隔距离大于8 km时各向变异均逐渐减小，各项变异的差异也逐渐变小.土壤AN含量在间隔距离小于3 km时，东西(E-W)变异小，其余各向变异相对一致，间隔距离大于3 km时，各向异性逐步增加，其中东-西(E-W)向增长幅度较大，当间隔距离大于8 km时各向均趋于同性，但呈各向异性增加趋势.综合土壤TN和AN含量的各向异性分析，可以得到土壤TN含量与AN含量随间隔距离的变化表现出十分相似的趋势，即东-西(E-W)变异程度占据主导作用，各向变异在达到一定间隔距离后又趋于同性，这与研究区地形从东向西依次为高山丘陵、漫川漫岗、平地，而且土地利用依次为林地、旱田、水田，各土地利用类型之间跨度较大有关，在较小采样尺度下变异程度小于较大采样尺度的变异程度.该结果可以进一步指导该尺度上未来研究过程中采样点的合理布设.

3.3 土壤中氮素的回归协同克里格插值

本文采用回归协同克里格法进行所有样点数据的空间插值，并对其空间分布状况和分布规律进行分析.在SPSS19.0中对所选外部指标与TN和AN含量进行相关性分析，以提取参与回归分析的指标(表 2，图 7)，该过程优先选取易于提取的地形指标因子，然后采用逐步回归法进行回归拟合，得到最终的回归方程(方程3，4)及回归趋势项，其中，TN和AN拟合方程的决定系数分别达到0.356和0.362.

(3)

(4)

式中：m(x)_TN为TN回归趋势项；SD为坡度；K_h为平面曲率；m(x)_AN为AN回归趋势项；H为高程；RD为距离居民点距离；K_v为剖面曲率.

在回归分析的基础上，对土壤TN和AN的回归残差进行描述性统计分析，TN的残差范围为-3.79～14.98 g·kg^-1，偏度和峰度分别为4.054、29.414，AN的残差范围为-79.86~160.02 mg·kg^-1，偏度和峰度分别为0.728、2.320，两者均通过K-S检验，呈现正态分布特征(图 5)，适合进行地统计插值，同时，半方差函数拟合模型分析表明TN和AN残差分别符合指数模型和高斯模型，且拟合效果很好，可以进一步指导克里格空间插值.

我们认为，回归残差的插值过程制约了空间插值的精度，尝试对回归残差与地形指标进行相关性分析显示，发现TN的回归残差与坡度呈显著性水平，AN的回归残差与高程呈显著性水平，故分别选取坡度和高程作为辅助变量对TN和AN回归残差进行协同克里格插值，交叉验证结果显示，TN和AN的残差插值的标准平均值误差(MSE)分别为-0.018和-0.024，标准均方根误差(RMSSE)分别为0.9750和0.9368，表明残差插值效果很好，将回归方程所得趋势项与残差插值结果叠加得到研究区土壤TN和AN的空间分布图(图 6).最后，对回归协同克里格插值进行的独立性检验结果表明，土壤TN和AN回归克里格插值的平均误差(ME)分别为0.024和-0.465，均方根误差(RMSE)分别为1.754和3.731，插值效果很好，能够反映和揭示研究区土壤TN和AN的空间分布特征.

从图 6可以看出：海沟河小流域土壤TN和AN含量的范围分别为0.023~5.706 g·kg^-1和10.19~274.51 mg·kg^-1.其中，土壤TN的高值区出现在研究区中部城镇周边，低值区出现在研究区西部水田集中分布及与旱地交接的区域，高值与低值区之间过度比较平缓，整体上流域内东部TN含量较高，向西依次呈现“高-低”交替的带状趋势，与土地利用方式在东西方向上的交替情形一致，最大值分布在中心乡镇附近.土壤AN空间分布与研究区地形地势表现关系密切，高值区出现在东部山区，低值区出现在中部旱地分布集中的区域，最小值则出现在北部海沟河上游的山谷中，整体上呈现出“两边高，中间低”和“山地高于平原，下游高于上游”的趋势.

3.4 土壤中氮素空间分布的主控因素分析

地形因素被认为是影响土壤理化性质空间分布的重要因素，其作用贯穿土壤的整个形成与发育过程，而随着人类对土地资源的不断开发，人类活动在土壤养分空间分布及再配置中起到越来越明显的作用(朱兆良，2008)，本文尝试按照基本地形指标因子(高程、坡度等)、外部环境因子(距居民点距离、距水系距离)与土地利用因子(坡位、土地利用类型)分类，对与土壤TN和AN分布关系密切的地形、外部环境、土地利用因子做进一步的分析，以量化土壤TN和AN的响应范围.

3.4.1 土壤氮素与基本地形因子的相关性分析

通过ArcGIS10.0平台直接从研究区DEM数据中提取主要地形指标因子，相关指标的计算过程参考已有研究内容(连纲等，2008；连纲等，2009；王充等，2014；李启权等，2014；黄魏等，2015)，利用SPSS19.0对土壤TN和AN含量与地形指标进行相关性分析(表 3)，可以看出：海沟河小流域耕地土壤TN与坡度(SD)呈显著正相关(r=0.449)水平，与平面曲率(K_h)呈显著负相关(r=-0.267)水平，表明：在坡度越大的地方土壤TN含量越高，坡度越小的地方TN含量越低，在平面曲率大的地方TN含量低，反之则较高；另一方面，土壤AN与高程(H)呈显著正相关(r=0.497)，与剖面曲率(K_V)呈显著负相关(r=-0.227).海拔高的地方土壤AN含量高，这可能与高海拔地区多森林覆盖，植被种类复杂，地表枯枝落叶较多，易形成有利于土壤TN和AN累积的地带性的气候条件有关，从而加速了土壤TN和AN在该区域的转换与形成，其次，地势高的地方土壤岩层在长期演变过程中易受风化作用的影响，促进了土壤TN和AN的长期累积过程；剖面曲率反映了地表坡度的变化率，其值越大的地方，一般受降雨冲刷越严重，造成水土流失的同时也导致了该地段土壤AN含量显著减少，进一步向低洼处累积.

此外，研究区土壤TN含量与高程(H)、剖面曲率(K_V)呈比较弱的负相关关系，与复合地形指数(CTI)、汇流动力指数(SPI)、沉积物运移指数(STI)呈比较弱的正相关关系；土壤AN含量与坡度(SD)、复合地形指数(CTI)、平面曲率(K_h)呈正相关关系，与汇流动力指数(SPI)、沉积物运移指数(STI)呈比较弱的负相关关系，其相关性表现较弱，影响作用有限，不再进行深入探讨.

3.4.2 土壤氮素与外部环境因子的函数关系分析

在SPSS19.0中进行回归分析中的曲线估计，获得海沟河小流域外部环境因子(距居民点距离、距水系距离)与TN和AN含量间的函数关系(图 6)，除水系距离与土壤TN和AN曲线拟合的R²值较低外，居民点距离与土壤TN和AN曲线拟合的R²值均很高，分别为0.966和0.982.进一步分析表明：居民点的分布与土壤TN和AN呈现明显的函数关系(图 6和图 6)，距离居民点近的地方土壤TN和AN含量明显较高，并且向远处逐渐减小，其中，距居民点距离与AN之间呈现出“倒数”的函数关系，该情形可能与农村居民点作为人类活动的聚集地，其日常的生产生活(燃烧取暖、废物堆积、污水排放)、畜牧生产等行为所产生的氮素在居民点周围产生累积密切相关.研究区土壤TN和AN与距河流距离间规律性较差(图 6a和图 6c)，但整体上也呈现“先增后减”的趋势，距离河流越近土壤TN和AN含量越低，但在距离河流1.2 km和2 km处分别出现最大值，这可能与河流作为流域内地表物质的“汇”有关，其附近地表物质运移的强度和频度均较大，小范围内的TN和AN受土壤侵蚀和地表水运动的影响，更多地进入水体，产生了一系列的环境问题，并且随着距离的增加，这种影响逐渐弱化，出现TN和AN的高值区.

综上，河流水系、居民点等外部环境因子与土壤TN和AN的空间分布之间均存在比较明显的临界范围，并且在一定范围内表现出较强的规律性，该结果在宏观层面上指导区域性施肥模式，制定污染防治措施服务，定量并定向地指导农业生产和生态环境协调发展方面具有一定的指导意义.

3.4.3 土地利用因子对土壤TN和AN的影响比较分析

研究区土地利用复杂，坡度变化较大，本文对不同土地利用的土壤TN和AN进行了均值比较分析(图 7):不同土地利用方式下，土壤TN含量平均值表现为：林地(FL)>玉米地(CL)>水稻地(RL)，3种土地利用方式下的土壤AN含量与TN含量表现特征差别较大，为玉米地(CL)>水稻地(RL)>林地(FL)，其中，林地AN含量明显低于水稻和玉米地，可能与其特殊的地表覆盖、地理环境及其形成的地带性气候有关，由于本次采样时间为10月份，研究区内玉米地、水稻地进行秸秆还田、焚烧的行为较多，并且研究区内畜牧养殖发达，更多的人畜粪便等有机肥进入田间，这些可能是造成水稻地和玉米地AN含量普遍较高的主要原因，且其作用过程远大于对土壤TN含量的影响.研究区不同坡位TN含量的平均值表现为坡顶(TS)>坡底(BS)>坡中(CS)>平地(FS)，AN含量的平均值表现为坡顶(TS)>坡底(BS)>平地(FS)>坡中(CS)，坡顶和坡底的TN和AN含量均高于坡中和平地，但差异不显著.这可能与土壤TN和AN的形成过程和随水土的运移结果有关，即坡顶土壤岩层在长期演变过程中易受风化作用的影响，促进了TN和AN的长期累积过程，同时又受雨水冲刷、地表径流等作用的影响，大量的TN和AN随水土运移并在坡底累积，出现富集现象(胡玉福等，2008；杨玲等，2012).总体上，研究区土地利用及坡位对土壤TN和AN含量的影响作用比较显著，并且存在较大差异.

4 结论(Conclusions)

1) 海沟河小流域土壤表层(0~20 cm)中TN和AN的含量分别在0.96~5.50 g·kg^-1和34.02~207.09 mg·kg^-1之间，处于较高和中等水平，空间变异程度均为中等变异程度；半方差函数分析显示土壤TN和AN的变程分别为900 m和1282 m，存在较强的空间自相关性，TN和AN的空间变异整体上受人类活动(耕作制度、施肥措施、秸秆处理)等随机因素影响小，受地形、成土母质等结构性要素的影响大，东-西(E-W)方向变异程度在空间变异中均占主导作用.

2) 回归克里格插值法可以有效提高土壤TN和AN的空间插值精度；插值结果表明海沟河小流域内土壤TN的高值区在该区中部城镇周边，低值区出现在西部水田集中分布的区域，高值与低值区之间过度比较平缓，向西依次呈现“高-低”交替的带状趋势，与土地利用方式在东西方向上的演替趋势相近.土壤AN空间分布与区域地形地势关系密切，高值区出现在东部山区，低值区出现在中部旱地分布较为集中的区域，整体上呈现出“两边高，中间低”和“山地高于平原，下游高于上游”的状况.

3) 坡度和平面曲率与土壤TN含量在空间上分别极显著性正相关和显著性负相关，高程、坡位与土壤AN含量在空间上显著正相关，剖面曲率与AN含量显著负相关；河流水系、居民点等非面状要素与土壤TN和AN的空间分布之间均存在比较明显的作用距离，并且在一定范围内表现出较强的规律性(图 6)；土地利用及坡位对土壤TN和AN含量影响作用显著，并且存在较大差异.以上研究结果可以在宏观层面上制定区域性的施肥模式，对定量定向地指导农业生产与环境治理具有重大的实际意义.