东北典型黑土区农田景观多尺度土壤养分时空分异研究进展

张少良，张海军，肖梓良，曲凤娟，王雪珊，霍纪平，张兴义，刘晓冰

（1.东北农业大学资源与环境学院，哈尔滨 150030；2.中国科学院东北地理与农业生态研究所，哈尔滨 150081）

土壤养分异质性具有普遍性，土壤养分时空分布格局和主要驱动机制是土壤养分管理重要依据[1-2]。研究表明土壤养分异质性给农田施肥带来困难，施肥不足降低作物产量和品质，施肥过量不利于作物生长，增加养分流失风险，特别是N、P流失对水环境质量构成威胁[3-4]。土壤养分时空异质性规律和主要驱动机制成为土壤学、生态学、环境科学和地理学等领域研究热点。土壤养分异质性程度，除与土壤养分自身性质有关，与外界自然环境和人文环境密切相关，特别是不同区域土壤类型、地形、气候和水文过程不同，异质性程度和主要驱动机制不同[5-6]。因此，不同区域土壤养分时空分布规律和主要驱动机制研究尤为必要。

黑土区以水蚀为主，长期面蚀和沟蚀严重影响区域农业可持续性。黑土土壤质量降低主要原因之一是土壤侵蚀，据报道每形成1 cm厚黑土至少需300～400年，而典型黑土区年均土壤流失量即达0.3～1.0 cm[7]。研究证实坡度、坡位、坡向均影响侵蚀过程，受侵蚀和沉积过程影响，部分地块露出黄土母质，土壤质量较差，土壤淤积较深改善或恶化土壤质量，改变土壤养分时空分布格局[8]。土地利用类型不同，地表覆盖和管理方式差异较大，影响土壤理化性质时空变化，如耕地比林地耕作频繁且施肥量高，林地和草地相对耕地有更好保土、保水效果[9]。耕作方式和作物种类影响坡地土壤理化性质异质性，如水保耕作方式和传统耕作方式下保土、保肥效果不同，如种植玉米和种植大豆其耕作方式、施肥量和施肥方式不同[10]。研究表明，不同尺度条件下景观异质性表现规律不同，异质性规律在各尺度之间相互推译关系存在不确定性[11]。因此，应充分考虑侵蚀、地形、土地利用、耕作方式、作物种类等因素作用，系统总结不同尺度条件下土壤养分时空分布规律和主要驱动机制，探索养分异质性尺度效应和影响机制。

本文系统分析当前黑土区不同尺度条件下土壤养分时空分异变化规律及其主要驱动机制研究进展，并展望未来黑土区土壤养分异质性研究需开展工作，旨在为黑土区土壤养分管理和区域生态环境保护等提供理论依据和技术支撑。

1 典型黑土养分异质性研究进展

东北典型黑土带跨越黑龙江、吉林和辽宁三省（见图1），总面积17万km2，开垦年限30～100年[12]。典型黑土区地势东北部高于西南部，海拔高度相差最高达200 m，典型地形为漫川漫岗；主要种植玉米和大豆；年均温度1.5℃，最高温度32℃，最低温度-37℃，年均降雨量500 mm，降雨集中在6～8月，因6月耕地植被覆盖度低，高强度降雨易发生强度土壤侵蚀；春季冻融频繁，积雪融化持续时间短（7～10 d），易发生融雪侵蚀[13]。开垦前典型黑土区黑土层厚度约80 cm，SOM含量60～120 g·kg-1，开垦后土层平均厚度约30 cm，耕层土壤SOM和土壤养分含量大幅下降[8,12]。

1.1 典型黑土带尺度土壤养分空间分异特征

典型黑土区土壤有机质（SOM）、全氮（TN）、全磷（TP）、速效氮（AN）、速效钾（AK）从南向北呈逐渐增高趋势；除全钾（TK）属于中等空间自相关强度外，其他指标均表现出较强空间自相关水平（表1）[14]；除TK属于弱空间变异程度外，其他养分均处于中度空间变异程度（0.1＜CV＜1.0）[15]。由于黑土区东西和南北海拔高度相差较小（200 m），南北跨度较大（900 km），垂直地带性可忽略。因此，典型黑土带SOM、TN、TP、AN、AK空间分布主要受太阳辐射影响，即呈现纬度地带性规律[16-17]。研究发现典型黑土带TK和速效磷（AP）未表现出明显纬度地带性规律，Zhang等认为主要受施肥历史影响，与土壤理化性质有关，如黑土成土母质中富含钾[14]，张少良等认为还与土地利用方式、种植作物种类和开垦年限等因素有关[18]，特别是黑土区钾肥用量较低，种植作物种类和年限差异较大，导致TK相对其他土壤养分表现出受随机因子影响更强特征。

1.2 市、县辖区尺度典型黑土养分空间分异特征

市、县辖区尺度，以黑龙江省哈尔滨市辖区（2×104km2）研究结果为例[19]（见表1），除AK属于强空间自相关强度外，SOM、TN、AN和AP均表现中等强度空间自相关水平[20]，相对黑土带尺度养分空间自相关程度减弱，受随机因子影响增强；SOM、TN、AN、AP和AK空间变异情况均处于中等强度水平（0.1＜CV＜1.0）[15]。分析认为除受太阳辐射影响（纬度地带性）外，还受季风等气候因素作用。研究区冬季刮西北风，夏季东南风，降水从西北向西南逐渐增加，导致哈尔滨市辖区土壤SOM和TN从西南向东北方向逐渐升高，呈现纬度地带性和海陆季风地带性交互特征。哈尔滨市辖区全量养分和速效养分（AN、AP、AK）表现明显非地带性规律，呈现养分“斑块化”特征，越靠近城市或乡镇，特别是靠近耕作历史较长城镇其速效养分含量较高，即“城市效应”，主要因为早期农民运输肥料移动半径有限，越靠近城市或乡镇土地开垦越早，施肥量和频率越高。海伦市土壤SOM和速效养分在部分乡镇附近也表现出明显“斑块化”特征[21]。因此，局部气候变化和土地耕作管理年限等改变土壤养分分布纬度地带性规律；养分分布在靠近城市或乡镇区域农田土壤N和P存在较高流失风险，是农业清洁生产及水体富营养化需要重点监管区域。

表1 不同尺度土壤养分含量、空间变异特征和空间自相关（所有样品均秋季采样）Table 1 Content,spatial variance,spatial autocorrelation of soil nutrients at various scales

1.3 村域尺度典型黑土养分空间分异特征

村域尺度，以黑龙江省绥化市前进乡光荣村（6.55 km2）研究结果为例[8]（见表1），研究表明，除耕层（0～20 cm）SOM呈强空间自相关水平外，TN和TP均表现中等强度空间自相关水平；村域尺度SOM和TN相对黑土带尺度养分空间自相关强度减弱，但相对市辖区尺度略有增强；SOM、TN和TP空间变异情况均处于中等程度变异水平（0.1＜CV＜1.0）[15]，低于黑土带尺度，却高于市辖区尺度。土壤SOM、TN和TP均表现出很强相关性，其空间分布规律相近，靠近流域地势较高区域含量相对较高。但由于村域尺度未考虑包含完全汇水单元（流域），地形（如坡向）和汇水通道等，难以从地形、水文和侵蚀沉积角度系统深入分析[22]。基于3个不完全汇水面积流域大量数据分析时发现，当同时考虑坡度和坡向时，研究区地形因子对土壤养分驱动作用更显著；相对坡度较陡区域，朝南坡向0～2%坡度区域和朝北坡向0～4%坡度区域养分含量相对较高，横坡垄作区域可显著提高土壤养分含量[8]。

1.4 流域尺度典型黑土养分空间分异特征

流域是相对完整汇水单元，几乎包含所有坡向和坡度，且可系统研究水文过程对养分异质性驱动作用[23]。以黑龙江省海伦市前进乡光荣小流域（1.86 km2）为例，利用随机采样法，分别研究0～20、20～30、30～40、40～50 和50～60 cm 土壤 SOM和土壤养分空间分布规律。研究发现相对黑土带、市辖区、村域尺度，流域内土壤养分空间分布规律和主要驱动机制改变（见表1）。0～60 cm土层土壤养分平均值在所有尺度中含量最低，其变异系数也处于较高水平。当仅考虑0～20 cm土层时，SOM、TP和AK平均值相对其他尺度略高，主要是因为研究区位于典型黑土带高纬度区域；TP和AK变异系数增加，因为流域尺度养分异质性主要驱动因子不完全同于其他尺度。P在土壤中移动性弱，施入土壤中P多数被固定[24]，黑土富含K素，当地旱作农业K肥使用量极低，K的有效性与土壤黏粒含量和水热分布密切相关[25]，因此可推断流域尺度TP和AK主要受侵蚀和沉积影响，说明流域尺度水文过程对土壤养分异质性影响权重较大。各土壤养分在不同土层深度空间分布情况见表2。

表2 流域尺度土壤养分含量、空间变异特征和空间自相关（N=122）Table 2 Content,spatial variance,spatial autocorrelation of soil nutrients at the watershed scale.

流域除50～60 cm土层SOM和0～30 cm土层TN处于中度自相关水平外，SOM、TN和AN在其他土层均处于强空间自相关水平，且SOM和AN空间自相关水平随土层深度增加而降低，而TN则呈相反趋势[25]，说明表土层TN受人为因素影响更大[26]。SOM、TN和AN空间变异强度多随土层深度增加而增加（见表2），主要受侵蚀和沉积影响[22-23]。研究区SOM和TN总体上处于较高水平，含量从表土层向深土层逐渐降低，表土层含量（0～20 cm）接近深土层含量（50～60 cm）2倍。受土壤侵蚀和沉积影响，水平方向上坡顶0～30 cm土层SOM和TN较高，坡中和部分坡底处较低，坡中处30～60 cm土层SOM和TN较低，而坡底，特别是靠近流域出口处最高。耕作方式和土地利用方式显著影响流域SOM和TN异质性，横坡垄作SOM和TN含量分别比顺坡垄作高34%和23%；退耕还林表层土壤SOM和TN小于耕地，但深土层高于耕地。道路和居民点也对SOM和TN分布影响显著，但影响权重小于地形因子。同样坡度是影响流域尺度SOM和TN主要因子，坡度对SOM和TN影响也与坡向有关。通常坡中处C和N储量最低，坡底处相对较高，特别是受沉积效应作用靠近流域出口处达到峰值[16]。C∶N影响土壤微生物活性进和土壤养分转化，改变土壤养分时空分布[24]。流域尺度C∶N空间异质性主要受结构因子影响，且C∶N与海拔高度呈正比，沿侵蚀沟或水道方向逐渐降低，在流域出口处达到谷值。C∶N同样受地形因子、耕作措施和土地利用等因素影响，通常C∶N在所有坡位均表现为随土层深度增加而增加，其中坡中处C∶N值相对较高，而坡顶和坡底较低。耕地深土层C∶N高于次生林地，大豆田深土层C∶N高于玉米田。横坡垄作同样因改变N流失过程，其C∶N通常高于顺坡垄作。研究区AN含量总体上处于较低水平，但在生态交错带附近（如，林地和耕地）及靠近侵蚀沟附近较高，横坡垄作区域AN显著增加。通常坡度与AN负相关；坡顶AN含量最高，其次为坡底；深土层AN受坡向影响较大。小流域坡中位置SOM和TN含量较低，而C∶N较高，需重点防控土壤退化，特别是土壤侵蚀的防控，还应结合水保措施综合增施N肥[16,22-23]。但在流域出口处养分含量较高，C∶N较低，特别是深土层SOM和TN含量高，C∶N较低，N存在较大流失风险[28]。

不同于SOM、TN和AN，流域尺度除0～20 cm土层AP处于中等强度空间自相关强度外，其他各土层TP和AP均处于强空间自相关水平，且TP空间自相关性强于AP[17]，由于P易被土壤固定、有效性低和移动性差[24]，区域以玉米和大豆种植为主，P肥施用量不高。耕地TP和AP通常从表土层向深土层逐渐降低，特别是耕层（0～20 cm）AP是深土层2倍；受土壤侵蚀影响TP空间变异程度从表土层向深土层逐渐增加，而AP空间变异程度则从30～40 cm土层分别向深土层和表土层降低，这与犁底层对土壤水和P迁移影响有关[17,26]。研究发现TP水平空间分布主要受结构因子影响，如受地形因子影响土壤侵蚀和沉积因子等，坡顶含量相对较高，坡中处最低，沿河道方向逐渐增加，在流域出口处达到峰值[17]。由于朝南坡向坡度较大，侵蚀程度相对较高[22-23]，导致朝南坡向TP含量总体低于朝北坡向。AP水平空间分布主要受土壤侵蚀、沉积、土壤水分、温度影响，南朝向坡面通常含量较高，沿坡向向下逐渐降低。同时由于施肥和入渗等在一定程度上影响AP纵向分布，必然也影响P水平分布[24]。土地利用和耕作措施均显著影响P空间分布，耕地表土层（0～20 cm）TP和AP均高于次生林地，而深土层则均低于林地；横坡垄作区域TP比顺坡垄作区域高22%。研究发现，当AP（Olsen-P）含量超过25～27 mg·kg-1时存在潜在流失风险，特别是靠近流域出口P流失风险最大[17,29]。

不同于SOM、TN、AN、TP和AP，受成土母质和地形因子等结构因子影响，流域尺度除30～40 cm土层AK处于中等自相关水平外，其他各土层均表现出强空间自相关水平，特别是靠近表土层和深土层处最高[25]。侵蚀和土地利用等因子显著影响流域AK空间分布，AK与侵蚀负相关，但与DEM呈正相关关系[25]。AK空间变异程度均处于中等水平，受作物吸收和土壤侵蚀沉积影响通常表土层和深土层空间变异程度相对较大。受地表水文过程影响[22-23]，流域内AK从林地向周边减小，沿侵蚀沟和水道方向减小，在流域出口处达最低。同时受太阳辐射影响[24]，北朝向坡面从坡顶向坡底处逐渐增加，南朝向坡面从坡顶向坡底逐渐降低。垂直方向上，AK通常从表土层向30 cm土层降低，然后又逐渐升高。由于区域施用K肥有限[8]，退耕还林地（20yr）0～60 cm土层AK含量均高于耕地。黑土区K含量丰富，根据其成土母质元素储量，土壤侵蚀损耗速率，作物生长消耗速率，K通常不作为多数作物生长限制因子[25]。

1.5 坡面尺度典型黑土养分时空分异特征

当研究幅度进一步降低时发现坡面尺度土壤TN处于中等和强度空间自相关水平，AN主要处于强空间自相关水平[30]。TN和AN多属于低等和中等程度空间变异水平。当仅施用底肥N，不施用追肥N时，作物生育期内NH4-N与NO3-N在不同坡位呈相反变化趋势；施用底肥和追肥时，作物生育期内NH4-N与NO3-N在不同坡位呈相近变化趋势；无论春季施足量N（超过120 kg·hm-1）底肥还是夏季施足量N（超过120 kg·hm-1）追肥，经历降水后，通常坡顶AN含量显著增加；生育期内AN峰值从坡顶向坡底逐渐运动，秋季后坡底AN含量最高，坡背最低[26]。生育期内AN异质性程度较高，特别是夏季降雨频繁季节最高。生育期内AN主要受结构因子影响，其分布主要受侵蚀和沉积影响，而侵蚀沉积主要受地形因子影响[24]，其次是不同坡位水热差异，以及作物生长对养分吸收作用差异[24,31]。考虑到夏季降雨频繁，特别是作物覆盖度较低季节（6～7月），AN含量通常也处于较高水平，是水土流失防治关键时期[9]，作物生育期内坡顶和坡底通常需施用足量N追肥，坡背处应该采用合理水土保持措施（如免耕、秸秆还田、生物黑炭还田、等高垄作等）后增加N肥投入，可有效提高黑土区坡耕地产量[30]。由于种植大豆施用N肥量较少，相对玉米连作，玉米-大豆轮作通常降低坡面土壤AN（低于3.6%）；受微生物对土壤C、N利用关系影响，秸秆还田通常降低坡面土壤AN含量[31]。

不同于TN和AN，作物生长期内坡面尺度土壤TP处于中等空间自相关水平，除施肥后AP异质性主要受随机因子影响外，其他时期AP空间分布主要受结构因子影响，空间自相关性较强[33]。TP和AP多处于低等和中等空间变异程度。作物生长期AP空间分布受施肥、作物种类、作物生长过程、降雨、侵蚀等多因子影响[31]，特别是6～7月主要受侵蚀影响，8～10月主要受作物生长影响，通常在作物生长末期坡底AP含量较低。春季坡顶AP含量较低，受土壤水热影响，作物生长中期（7～8月）各坡位和处理土壤AP含量显著高于苗期和成熟期，但各坡位和各处理间差异较小[32]。作物生长期内，秸秆还田和轮作措施对坡面AP时空分布影响较弱或在短时期内影响不显著。种植玉米施用N肥较多，受微生物活性影响，比种植大豆能更有效提高坡面土壤AP含量[33]。

综上所述，水土流失是作物生长期坡面尺度N、P流失主要原因，7月份是土壤N、P流失重点防治时期，坡面均存在较大N、P流失风险。

2 冻融驱动下典型黑土养分时空分异特征

2.1 冻融驱动下流域尺度典型黑土养分时空分异特征

冻融增强流域尺度0～20 cm土壤TN空间自相关性，但对空间变异影响较小，增加AN空间自相关性与空间变异性[27]。研究发现冻融后流域TN主要在3个区域明显增加:①靠近流域出口，特别是0～5 cm和5～10 cm土层，②流域部分坡顶坡度较缓和位置，③农林交错带部分位置。AN在3个区域明显增加:①沿侵蚀沟两侧位置至流域出口，②林地种植区附近，③ 坡底位置。流域尺度平均值分析表明，冻融后，0～20 cm土壤TN减少5.7%，而AN减少7.6%且达到显著水平（P＜0.05）。冻融后，除坡下0～10 cm土壤AN上升1.1-8%外，土壤TN和AN元素在坡上、坡中、坡下均表现为下降趋势。冻融后，TN和AN在面南和面北坡面0～20 cm均表现为下降趋势。与面南坡面对比，冻融后面北坡向土壤TN下降1.8%，而AN则相反；面南坡向AN相对较高，相比于面北坡高5.5%。与冻融前相比较，冻融后0～5、5～10和10～20 cm土壤中TN和AN平均值表现为下降趋势，分别下降4.6%～6.4%和4.7%～10.1%。冻融后林地0～20 cm土壤TN下降幅度高于农田7.5%，农田AN减少10.3%，而林地AN却上升3.7%。冻融后，大豆田和玉米田TN和AN均表现为下降趋势，除玉米田0～5 cm外，两者AN下降幅度均高于TN。冻融后，无论横坡还是顺坡耕作方式下，土壤TN和AN平均值均表现为下降趋势，其中AN下降幅度高于TN。相比于顺坡，横坡垄作在0～10 cm可减少2.7%～8%AN和1.1%～3.5%TN流失或转化[27]。总体来说，冻融减少流域土壤N含量，其中融雪侵蚀、径流输运以及氮转化在冻融过程中具有重要作用。

流域尺度，冻融后除5～10 cm土壤AP空间自相关性呈微弱下降趋势（5.7%）外，其他土层TP和AP空间自相关性均存在不同程度增加[34]。冻融过程增加TP空间变异，且变化幅度随深度增加而增加；而冻融后除10～20 cm土层外，AP空间变异呈下降趋势。冻融过程减少流域尺度下0～20 cm耕层TP（9.52%）和AP（3.79%）含量；在垂直方向上，通常冻融前TP随土层深度增加而升高，而冻融后则相反；冻融前后AP含量均随土层深度增加而递减。研究发现在空间上冻融后流域部分区域TP和AP均有不同程度增加或降低。冻融后TP主要在流域3个区域明显增加:①靠近流域侵蚀沟附近区域，②农林交错带，③流域南部部分横坡垄作区域。AP在3个区域明显增加:①流域沿侵蚀沟附近位置，②林地种植区附近，③居民点附近区域。这一变化主要受融雪侵蚀、土地利用类型、地形因子、土壤性质以及冻融环境等共同影响。冻融后不同坡位P均存在不同程度减少，其中坡上和坡下是冻融过程P主要流失风险区，流失量分别为0.087和0.065 g·kg-1。不同土地利用类型下，农田TP减少0.09 g·kg-1，而林地增加0.04 g·kg-1。此外，冻融过程增加林地AP释放和沉积，冻融后0～20 cm AP增加5.98 mg·kg-1。流域尺度下，顺坡耕作，玉米种植区、南向坡、低植被覆盖区以及高积雪覆盖区是冻融过程P流失重点防控区域[25]。

2.2 冻融驱动下坡面尺度典型黑土养分时空分异特征

坡面尺度冻融前TN处于中度或强空间自相关水平，除个别年份无覆盖10～20 cm土层变化较大外，冻融后有覆盖（土壤表层被带有透气、防降水、侵蚀圆柱型小桶覆盖，排除大气沉降，排除水土流失，验证冻融过程养分二次迁移结果）或无覆盖处理土壤TN通常仍处于中度或强空间自相关水平，且两者间因年份不同规律也不同[23]。冻融前AN处于中度或强空间自相关水平，冻融后无覆盖或有覆盖处理土壤AN通常仍处于中度或强空间自相关水平，同样两者间因年份不同规律也不同。冻融前，典型坡面土壤TN和AN主要受到地形因子影响，旋翻起垄降低土壤0～10 cm土层中TN和AN空间自相关性。受侵蚀和硝化反硝化作用，冻融通常降低坡面AN含量，但TN因年降雪量和气温等因素影响年际间差异较大，各坡位间差异也较大。覆盖条件通常可有效减少N流失，增加土壤N有效性，径流、土壤侵蚀和水分垂直方向上二次迁移是冻融过程影响坡面N增减主要因子[23,26]。

冻融前坡面TP表现为中等强度空间自相关水平，覆盖改变各土层TP空间自相关程度[24]。冻融后，覆盖可增加或维持TP空间自相关性，而无覆盖则相反。冻融前TP主要受结构因子影响，冻融后有覆盖主要受结构因子影响，而无覆盖表土层0～5和10～20 cm TP主要受随机因子影响。冻融前AP表现出纯块金效应，覆盖降低AP空间自相关，而无覆盖增加AP空间自相关；冻融后除有覆盖0～5 cm土层处于强空间自相关水平外，其他土层无变化。因为冻融前耕地作旋翻起垄，此时AP主要受随机因子影响，但冻融后无覆盖0～20 cm和有覆盖0～5 cm土层主要受结构因子影响，但有覆盖5～10 cm仍主要受随机因子影响。冻融后除轮作不还田坡顶外，土壤AP含量普遍高于冻融前土壤AP含量，特别是轮作不还田-无覆盖-坡脚处理处高出1倍，冻融增加耕层P有效性。冻融结束后所有处理无覆盖土壤坡顶TP降低，其他所有处理各坡位均增加，有覆盖土壤TP显著增加（37.1%）。这是由于冻融过程中，增大蒸发量以及合适土壤含水量有利于土壤P向上运动，但受侵蚀和沉积影响表层土壤TP含量在部分区域降低或增加。坡背处土壤P含量增加，因土壤容重较大，适当增加容重有利于P向表土层运移。覆盖（秸秆还田和小桶覆盖）虽可减少蒸发（减少空气流动，但增加土壤温度），防止大气沉降，但因无侵蚀或侵蚀量少而降低土壤P流失。无覆盖蒸发量虽大，P向上迁移量大，但侵蚀严重降低表土层P含量。耕种方式即玉米大豆轮作和玉米连作对冻融期间耕层土壤TP动态变化影响较小[24]。

野外观测发现，冻融过程多数坡位细菌和真菌生物多样性降低，特别是坡顶和坡底细菌多样性，坡顶和部分坡背真菌多样性均降低；冻融过程优势科细菌数量降低，如Pseudomonadaceae和Bacillaceae数量均显著降低，优势科真菌数量降低，如Cystofilobasidiaceae，冻融过程降低微生物分解土壤有机物能力，特别是坡顶细菌对N、P循环作用，为冻融后坡顶N、P通常向上迁移量较少原因[35]。综合不同尺度野外观测以及室内模拟结果，初步结论证实冻融后通常降低大部分区域土壤TN、TP、AN和AP含量，但增加局部区域TN、TP、AN和AP含量，如横坡垄作、林地和部分坡顶位置，蒸发和冻融过程毛管力作用下土壤水势发生变化，引起养分向上迁移，通常水分含量越高养分向上运动越剧烈，而水土流失严重降低土壤养分，特别是坡中位置。土壤N硝化-反消化过程也增加冻融过程N损失[36]。

3 典型黑土养分空间分异特征研究方法进展

以往土壤养分空间异质性分析多是针对某个区域先采样，测定养分指标如满足正态分布和二阶平稳假设便可作克里格插值和空间分析，当数据无法满足正态分布，即使通过不同方法转换也很难得到准确空间插值或分析结果，受时间限制也很难再次补充采样。因此，空间土壤采样原则上应增加采样点，但采集土样过多，增加劳动强度、采样和分析成本。为此，有必要优化特定区域土壤采样方法，即明确采样最佳样本量、粒度或密度阈值，优化采样路径；有必要改进空间插值方法，借助相关技术和算法提高有限样本插值精度。

3.1 典型黑土养分异质性调查采样方法

研究表明，黑土区采样过程结合网格采样法和随机采样法可有效减少采样过程损失空间信息。当研究区面积较小，地形复杂，土地利用类型多，异质性程度高时，可采用标准和分层网格采样法；当采样幅度较大时，异质性程度较低，土地利用类型单一，地形起伏变化较小，可采用随机采样法，需注意样点均匀度。黑土区采样粒度通常不应低于0.025 km2，样本容量多于120个；当研究区面积增大时采样粒度可适当降低，但样本容量不应减少，值得注意的是，当采样幅度和粒度增加时反映异质性结果通常也发生变化[37]。关于采样路线设计和布点方位，要根据土地利用、地形、特殊环境位置等实际情况综合设计；因典型黑土主要分布在漫川漫岗地区，不同区域坡度、坡长、海拔、水系密度、水流方向等均不同，特别是坡度较大区域不同坡位侵蚀和沉积效果不同，设计采样路线时在缩短采样路线距离、降低风险、降低难度、提高效率前提下，还应兼顾坡顶、坡背和坡底3个关键位置作系统采样[37]（见图2）；采样深度通常依据调查目的而定，如调查耕层土壤养分空间异质性时，要依据犁底层位置确定耕层深度，多数区域为0～20 cm，部分区域为0～15或0～25 cm[36-37]。

3.2 典型黑土养分异质性统计分析方法和空间插值方法

普通克里格、简单克里格、反距离权重等插值方法是早期常见插值方法，特别是普通克里格（OK，Ordinary kriging）方法，因其建立在半方差函数基础上，提高土壤养分空间模拟和预测精度，应用广泛。随着统计学方法革新胶土壤养分空间插值精度要求不断提高，出现一些新土壤养分空间插值方法，如协克里格（CK）、多线性回归（MLR）、地理回归（GWR）、回归克里格（RK）、地理加权回归克里格（GWRK）等方法，但在不同区域，不同插值方法效果存在较大差异。以AN、AK和AP为例，利用其与SOM和TN显著回归关系，借助CK插值方法在一定程度上提高速效养分空间插值精度[20]；以TP空间插值为例，在典型黑土区4个不同区域（土地利用结构、地形复杂程度、样点分布等均不同），借助坡度或亮度指数作为辅助变量，利用RK和GRK插值方法，显著提高土壤TP空间插值精度。在插值过程中还发现地形起伏、采样点均匀度、采样点密度等均不同程度影响插值结果。研究结果同时证实，在典型黑土区开展土壤养分空间异质性分析过程中，当样点数量足够多，样点分布足够均匀，通常应用OK插值方法也可获得较高模拟精度，不需辅助变量，分析过程相对简单[38]。

4 展望

在长期人为因素作用下，典型黑土区耕地土壤养分时空异质性与尺度密切相关，不同尺度条件下土壤养分异质性格局、变异和主要驱动过程不一致，证实黑土区土壤养分模拟过程中不同尺度之间具有不可直接推绎性。黑土带尺度主要受太阳辐射影响。市县辖区尺度受太阳辐射和海陆季风影响。村域尺度主要受土地利用、坡度、坡向和坡位等地形因子影响。流域尺度除受土地利用、地形因子影响外，还受与耕作、侵蚀、沉积、水文过程等因子影响。作物生长期，受施肥、降雨、侵蚀、水热、作物生长等因素综合作用，坡面尺度土壤速效N和P呈周期性变化。冻融过程，受耕作措施、地形、侵蚀、沉积、土壤水分和容重等多因素影响，坡面尺度土壤N、P时空格局成规律性变化，流域尺度上表现出局部增加或降低现象。阐明不同尺度土壤养分空间分布格局和主要驱动机制，是正确认识土壤养时空分异格局变化本质需求，也是合理施肥和环境质量管理重要依据。

当前典型黑土土壤养分时空分异规律和主要驱动机制研究取得重要进展，但仍有待深入研究:①较大时间尺度上，特别是跨越几年或几十年尺度土壤养分空间分异规律和主要驱动机制；②较小空间尺度上，特别是不同作物根际范围土壤养分时空变化规律和主要驱动机制研究；③特定环境下或复杂环境背景下土壤养分时空分异规律和主要驱动机研究，如特殊地理位置:农林交错带、水陆交错带、侵蚀沟影响环境条件下等，以及气候变化背景下土壤养分时空变化特征和驱动机制；④多尺度土壤养分时空分布特征和主要驱动机制研究取得一系列进展，但如何精准实现空间异质性时间序列模拟和预测还需深入研究；⑤土壤养分空间异质性变化空间统计方法和模型还需进一步改善，应重点研究降低采样样本数前提下构建更简单、高效、精准、通用的空间插值方法。