基于空间自相关的县域尺度土壤侵蚀抽样方法研究

陈 超，齐 斐，徐雁南*，黎家作，赵传普,苏新宇

(1.南京林业大学林学院，江苏 南京 210037；2.水利部淮河水利委员会淮河流域水土保持监测中心站,安徽 蚌埠 233001)

土壤侵蚀是一个受多种因素影响的生态环境过程。摸清土壤侵蚀状况，因地制宜地采取有效防治措施，对维护区域生态环境、实现社会经济的可持续发展具有重要的意义。目前，在土壤侵蚀状况调查中，通常采用抽样单元推算、经验模型、遥感调查以及多种手段结合应用等方法[1-5]。其中，美国主要采用抽样调查和模型法(USLE、WEQ)，澳大利亚及欧洲主要采用专家法和模型法(USLE、CEMS)[6-9],而我国多采用综合评判法以及CSLE模型法。近几年，中国土壤流失方程(Chinese soil loss equation, CSLE)在全国得到大范围应用，众多学者将其与GIS或抽样单元相结合，进行区域土壤侵蚀状况评估[10-12]。根据相关研究，在区域尺度上以GIS为基础采用栅格计算法/全域计算法可以获取详细的土壤侵蚀状况，但存在模型因子精度不一致、工作量大，且其结果受土地利用精度影响等问题。采用抽样调查与经验模型相结合的方法仍然是快速获取区域土壤侵蚀状况的重要方法之一。段倩等[13]以沂水县为例，采用分层系统抽样法，对比分析了1%和4%抽样密度下不同推算方法对土壤侵蚀状况监测结果的影响；邹丛荣等[14]探讨在分层系统抽样方法下，不同抽样密度和推算方法对估算县域土壤侵蚀状况监测结果的影响；李子轩等[15]对比全域覆盖计算与4%密度抽样单元推算方法对全县土壤侵蚀计算结果的差异；吴迪等[16]采用分系统层抽样法，提高山丘区抽样密度，评价土壤侵蚀状况；朱梦阳等[17]采用分层不等概系统空间抽样方法布设抽样单元，探究土壤侵蚀快速调查方法。以上研究针对抽样单元的土壤侵蚀预测进行研究，但仅讨论了基于系统抽样的抽样密度和推算方法对结果的影响，未对其他抽样方法和插值方法进行探索。因此，本研究以沂蒙山泰山国家级重点治理区中的典型县——蒙阴县为研究区域，基于空间自相关理论，综合应用遥感技术、野外调查和地统计方法，综合考虑抽样数量和精度以及空间预测精度，探讨适宜该区域的空间抽样和插值方法，以期为沂蒙山区县域尺度土壤侵蚀状况快速调查提供技术支撑。

1 材料与方法

1.1 研究区概况

山东省临沂市蒙阴县(117°45′～118°15′E， 35°27′～36°02′N)位于山东省中南部、沂蒙山区北部，面积1 602 km2，为典型的北方土石山区。境内海拔131～1 105 m，主要为低山丘陵,山丘区占97.77%；岩石多为石灰岩和页岩,土壤主要有粗骨土、棕壤、褐土、红黏土；多年平均气温12.8 ℃，多年平均降水量700 mm，降水主要集中在6—8月，属暖温带季风大陆性气候；境内有梓河、东汶河、蒙河等3条较大的河流，以及山东省第二大水库云蒙湖，属淮河流域沂河水系；植被属暖温带落叶阔叶林区，现多为人工植被，防护林树种主要有侧柏(Platycladusorientalis)、刺槐(Robiniapseudoacacia)、油松(Pinustabulaeformis)，经济林树种主要有桃(Amygdaluspersica)、苹果(Maluspumila)、核桃(Juglansregia)等，自然生长的灌木与草本植物主要有黄荆(Vitexnegundo)、胡枝子(Lespedezabicolor)、三裂绣线菊(Spiraeatrilobata)等。

1.2 数据与处理

采用2020年4月2 m分辨率的GF-1和GF-6卫星遥感影像，以eCongnition 9.0软件为操作平台进行人机交互解译，参考GB/T 21010—2017《土地利用现状分类》[18]，获得土地利用类型和水土保持措施[19]，土地利用类型划分为：耕地(平缓地、坡耕地、梯田)，园地(平地果园、坡地果园、梯田果园)，林地(防护林、用材林、灌木林地、其他林地)，其他草地，建设用地(城镇建设用地、农村建设用地、采矿用地、其他建设用地)，交通运输用地(农村道路、其他交通用地)，水域及水利设施用地和其他土地。基于2017—2019年250 m空间分辨率MODIS遥感影像，采用参数转换法，计算24.5个月的平均植被覆盖度。降雨量采用沂蒙山区88个雨量站点1986—2017年逐日降水量≥10 mm的降雨数据；土壤数据采用山东省1∶50万土壤图，裁剪获取研究区土壤类型分布；地形数据采用1∶1万地形图，生成10 m栅格DEM；野外调查于2020年8月进行，采集土地利用、水土保持措施等数据[20]，用以修正土地利用数据并计算土壤侵蚀模数。土壤侵蚀模数根据CSLE模型[21]计算：

A=R·K·L·S·IBET。

(1)

式中：A为土壤流失量，t/(hm2·a)；R为降雨侵蚀力因子，(MJ·mm)/(hm2·h·a)，根据1986—2017年沂蒙山区88个雨量站点降雨数据采用普通克里金插值获取[22-23]；K为土壤可蚀性因子，其单位为(t·hm2·h)/(hm2·MJ·mm),通过野外采集土壤样品，采用Williams模型计算，并根据沂蒙山区径流小区数据修正；L和S为坡长因子和坡度因子(均无量纲)，根据1∶1万地形图，采用刘宝元修正算法[23]计算；IBET为水土保持因子(无量纲)，参考第一次全国水利普查水土保持情况普查措施因子值，采用径流小区数据修正[23]。

1.3 空间抽样方法

空间抽样基于空间自相关理论进行。空间自相关是指同一变量在不同空间位置上的相关性[24]。Moran’sI是描述空间相关性常用指标，其计算方法见文献[25-27]。

Moran’sI取值范围一般在-1～1之间，小于0表示负相关，大于0表示正相关，等于0表示不相关；绝对值越接近1代表单元间性质越相似，接近0则代表单元间不相关[25-27]。

与传统抽样方法相比，空间抽样是在样本单元间自相关显著的前提下计算样本容量，可有效减少样本冗余及减少抽样成本[25-26]。在传统抽样方法计算的样本容量的基础上，空间抽样样本容量计算公式[27]为:

n=nClassic·(1-c)。

(2)

式中：n为空间抽样方法样本量，nClassic为由传统简单随机抽样计算而得的样本容量[28]，c为总体相关系数。空间抽样中均值的方差估算参考文献[26]。

1.3.1 抽样单元尺寸确定方法

基于“抽样单元间相互独立”原则[26]，分别建立200 m×200 m、300 m×300 m、400 m×400 m、500 m×500 m和600 m×600 m等5种抽样单元，剔除面积小于50%的网格，并采用CSLE模型计算各个抽样单元土壤侵蚀模数；分别统计不同抽样单元的土地利用图斑数量特征和Moran’sI指数，利用标准化统计量Zscore对空间自相关的显著性水平进行检验，综合考虑图斑数量特征和空间自相关属性确定最优抽样单元尺寸。

1.3.2 空间抽样方法优选

研究选择空间随机抽样(spatial random sampling,SRS)、空间系统抽样(spatial systematic sampling,SSS系统)和空间分层抽样(spatial stratified sampling,SSS分层)3种抽样方法(图1)，为避免偶然性，减小因样本抽样过程产生的数据不稳定性，在可靠性95%、精度90%的水平下，分别抽取5套单元样本，分析比较3种方法的样本容量和抽样精度平均值，确定最优抽样方法。

图1 研究区野外调查单元及不同抽样方法的样点分布

1)空间随机抽样：从总体N个单元中，随机、独立抽取n个单元组成样本估计总体参数的方法[28]。样本容量计算方法参考文献[27-30]。为了降低野外调查难度，野外调查时考虑道路交通的可达性及调查难度，将样本容量扩大10%。本研究采用随机数法[31]，抽取5套样本单元，分别计算样本特征参数及其平均值。

2)空间系统抽样：又称等距离抽样，从总体N个单元中，随机确定起点后，按照间隔k抽取n个单元组成样本，用于估算总体。其样本容量n的确定、总体估算和误差估计公式与空间随机抽样一致[30]。取样时，首先随机选择一个样本作为起点，然后根据间隔大小在纵横双向上等间隔选择样本点，抽取5套抽样单元样本。

3)空间分层抽样：指按照某种特征或标志将调查总体划分为若干层次或类型，然后在各层内独立、随机地进行抽样的抽样方法。①以坡度为分层标志，参考《坡度分级标准与编码表》[23]，结合研究区坡度特征，将研究区划分为<2°、≥2°～5°、≥5°～8°、≥8°～15°、≥15°～25°、>25°等6层(图1)。②确定样本容量时，各层最优样本数根据最优分配原则计算，总体估算和误差估计公式见文献[27-28]。③采用与简单随机抽样相同的方法在6个样本层中抽样，获得各个层内的样本单元。

1.4 土壤侵蚀预测的插值方法

在获得最优抽样方法的基础上，选择确定性方法和地统计方法进行水土流失空间分布预测，以栅格计算法(grid calculation,GC)为基准值进行精度评价。其中，确定性方法为反距离权重法(inverse distance weight, IDW)和张力样条函数法(spline function method,SFM)，地统计方法为普通克里金法(o-Kriging)和协同克里金法(co-Kriging,辅助变量为坡度)[24]。地统计插值时采用正态QQ图对样本数据进行正态分布检验，根据检验结果对其进行变换，并选择合适的变异函数模型；确定性插值数据则不进行正态变换处理。

栅格计算法(GC)根据野外调查单元数据面积加权平均，采用CSLE模型法进行土壤侵蚀计算。参考SL 665—2014《北方土石山区水土流失综合治理技术标准》[32]对土壤侵蚀强度分级，结合研究区土壤侵蚀模数特征，将土壤侵蚀强度分为微度[<200 t/(hm2·a)]、轻度[≥200～1 000 t/(hm2·a)]、中度[≥1 000～2 500 t/(hm2·a)]、强烈[≥2 500～4 000 t/(hm2·a)]和极强烈[>4 000 t/(hm2·a)]，其中，建设用地(不包括采矿用地)、交通运输用地、水域及水利设施用地、其他土地等土地利用类型不参与侵蚀预测，默认为微度侵蚀。

精度评价时采用交叉验证法，常用指标有平均误差(ME，式中记为δME)、均方根误差(RMSE，式中记为δRMSE)、平均标准误差(ASE，式中记为δASE)、平均标准化误差(MSE，式中记为δMSE)和均方根标准化误差(RMSSE，式中记为δRMSSE)；采用相对误差(σ)作为评价各级土壤侵蚀分布的指标，计算公式如下：

σi=|Vi-Ai|/Ai×100%。

(3)

式中：σi为第i级土壤侵蚀强度面积相对误差，Vi为空间插值方法计算的第i级土壤侵蚀强度面积，Ai为基准值。

2 结果与分析

2.1 土壤侵蚀抽样方案的确定

2.1.1 抽样单元尺寸的选择

Moran’sI是衡量空间自相关的常用指标，本研究抽样单元方案的空间自相关指数见图2。由图2可知5种抽样单元尺寸方案的变化趋势，即在200 m×200 m和500 m×500 m区间内随着抽样单元尺寸增大，对应的Moran’sI和Zscore则相应减小，单元尺寸大于500 m×500 m后Moran’sI呈上升趋势，Zscore减小。这说明，在200 m×200 m和500 m×500 m区间内，随着抽样尺寸的增加，空间相关性越弱；抽样单元尺寸在400 m×400 m和600 m×600 m之间，对应的Moran’sI斜率较小，说明其空间自相关性变化不大。同时，抽样单元尺寸越大，每个单元的平均图斑数和图斑种类越多，600 m×600 m的图斑数量是400 m×400 m图斑数量的1.26倍。因此，综合考虑后选取400 m×400 m作为最优的抽样单元尺寸，此时空间自相关性较弱，且图斑数量相对较少。

图2 不同抽样单元方案的空间自相关指数和平均图斑数

2.1.2 抽样方法的选择

在可靠性95%、精度90%的水平下，样本容量从高到低依次为空间系统抽样(278个)、空间随机抽样(252个)、空间分层抽样(150个),其中空间分层抽样各层样本量分别为6、16、19、35、53、17。抽样精度为空间分层抽样(96.69%)>空间系统抽样(94.33%)>空间随机抽样(90.10%)。总体来看，空间随机抽样样本容量居中，抽样精度最低；空间系统抽样样本量最多，抽样精度居中；空间分层抽样样本量最少，抽样精度最高。因此，空间分层抽样优势最明显。

2.2 基于普通克里金的不同抽样方法对土壤侵蚀计算结果的影响

以3种抽样方法获取的水土流失区域土壤侵蚀强度均以轻度为主，中度及以上侵蚀面积较小(表1，图3)。其中，基于空间随机抽样的水土流失面积为610.34 km2，占全县面积的38.02%；基于空间系统抽样的水土流失面积为578.22 km2，占全县面积的36.04%；空间分层抽样的水土流失面积548.58 km2，占全县面积的34.21%。与基准值比较，空间简单随机抽样计算水土流失面积与基准值误差较大，各级土壤侵蚀强度相对误差也较大，且随着侵蚀强度的增加而增大，尤其是强烈及以上侵蚀相对误差高于85.80%；基于空间系统抽样的水土流失面积与基准值相对误差为14.10%，中度侵蚀相对误差为17.07%，无强烈及其以上侵蚀；基于空间分层抽样的水土流失面积最接近基准值，相对误差8.25%。

表1 普通克里金法下不同空间抽样方法土壤侵蚀强度结果

图3 基于栅格计算法和空间插值法的研究区空间分层抽样土壤侵蚀强度分布

从空间分布来看(图3)，基于普通克里金法和空间分层抽样的轻度及以上土壤侵蚀区域主要在县域东北部和西南部，呈集中分布状态，而栅格计算法结果呈均匀异质性零散分布，存在明显差异。因此，综合考虑抽样样本数量以及空间预测结果，仍以空间分层抽样结果最优。

2.3 基于空间分层抽样的不同插值方法对土壤侵蚀计算结果的影响

由不同插值方法的交叉验证结果(表2)可知，普通克里金法和协同克里金法的ME、RMSE明显较小，说明二者预测误差较小；普通克里金法的ME、MSE略小于协同克里金法，而后者的σRMSSE和|δASE-δRMSE|略小于前者，说明前者拟合模型较优，后者拟合较有效。

表2 空间分层抽样下不同插值方法的交叉验证结果

基于空间分层抽样方法采用不同插值方法计算的蒙阴县水土流失状况的结果见表3。与栅格计算法结果相比，以反距离权重法计算的水土流失面积相差高达379.00 km2，轻度侵蚀面积相差最大，高达326.89 km2，其余各级土壤侵蚀强度面积相差0.56～78.87 km2；以样条函数法计算的水土流失面积相差高达126.89 km2，主要对中度和轻度侵蚀高估，分别相差65.75 km2和54.16 km2。普通克里金和协同克里金法插值结果差异较小，水土流失面积仅相差8.09 km2，占全县面积的0.50%，但各级土壤侵蚀面积差异较大。

表3 空间分层抽样下不同插值方法土壤侵蚀强度结果

与协同克里金法相比，普通克里金法结果更接近栅格计算法，水土流失面积计算结果相差41.82 km2，其中，轻度侵蚀相对误差最小为7.63%，强烈侵蚀面积为0.75 km2，相对误差最大为55.05%，其余各级相对误差约40%；基于协同克里金法的水土流失面积比基准值高49.90 km2，相对误差为9.85%，但中度侵蚀计算相对误差大于100%，达172.74%，其余各级相对误差为5.62%～28.81%。总体看，普通克里金法和协同克里金法较优。

从空间分布看(图3)，4种方法插值结果显示研究区水土流失区域均集中分布在蒙阴县东北部、中部和西北部，但协同克里金法轻度及以上土壤侵蚀分布区域与其他方法差异明显，土壤侵蚀区域分布相对更为细致，这是源于辅助变量坡度的影响。

2.4 空间分层抽样的样本典型性分析

根据土壤类型统计可知，蒙阴县土壤类型23种，空间分层抽样共涉及19种土壤类型，二者均以酸性粗骨土为主，其次为红黏土和石灰岩钙质粗骨土，但两者面积比例存在差异(图4)。空间分层抽样单元内酸性粗骨土和石灰岩钙质粗骨土所占比例均高于蒙阴县均值，绝对误差高于9.77%;红黏土所占比例低于全县，绝对误差为3.77%。

A.土壤类型soil type：Ⅰ.酸性粗骨土acid skeleton soil;Ⅱ.砂质非灰性河潮土sandy non-gray moisture soil;Ⅲ.棕壤brown soil;Ⅳ.麻砂棕壤性土sandy brown soil;Ⅴ.红黏土red clay;Ⅵ.灰质淋溶褐土lime leached cinnamon soil;Ⅶ.黏质中层灰质淋溶褐土clay middle calcareous leached cinnamon soi;Ⅷ.石灰岩钙质粗骨土limestone calcareous skeleton soil;Ⅸ.红土red soil;Ⅹ.洪冲积潮棕壤flood alluvial moisture brown soil。B.土地利用类型landuse type：1.平缓地flat land;2.坡耕地sloping land;3.梯田terrace;4.平地果园flat orchard;5.坡地果园sloping orchard;6.梯田果园terrace orchard;7.防护林shelter belt;8.用材林timber forest;9.灌木林地shrubbery;10.其他林地other woodland;11.其他草地mining land;12.建设用地construction。C. 坡度slope: ⅰ.0°～2°;ⅱ.≥2°～5°;ⅲ.≥5°～8°;ⅳ.≥8°～15°;ⅴ.≥15°～25°;ⅵ.≥25°。

在CSLE模型因子中，研究区K均值为0.009 4 (t·hm2·h)/(hm2·MJ·mm)，空间分层抽样K均值为0.009 3 (t·hm2·h)/(hm2·MJ·mm)，相差0.000 1 (t·hm2·h)/(hm2·MJ·mm)。因此，虽然二者涵盖的土壤类型存在差异，但对K因子影响不大。

从土地利用看，空间分层抽样的样本单元内土地利用类型与蒙阴县全县整体情况存在差异，全县土地利用类型以平缓地、旱梯田、梯田果园和防护林为主，分别占14.07%、18.99%、17.02%和14.31%，空间分层抽样结果则以旱梯田、梯田果园和防护林为主，分别占20.59%、20.15%和22.15%，空间分层抽样平缓地和平地果园比例分别比全县低6.91%和2.76%，梯田果园和防护林面积比例分别比全县高3.13%和7.84%，其他土地利用类型面积与全县相似，绝对误差小于1.9%。空间分层抽样主要土地利用类型与全县基本保持一致，但面积比具有明显差异，这可能是造成水土流失面积估算值偏大的原因。

根据坡度统计，研究区坡度以0°～2°和≥8°～15°为主，分别占比23.94%和20.79%，平均坡度10.20°。空间分层抽样样本主要分布在≥15°～25°和≥8°～15°的区域，0°～2°的区域分布最少；样点平均坡度为14.76°，比研究区高4.56°，这也可能会造成水土流失面积估算偏大。

3 讨 论

1)在抽样方案设计中，研究设计抽样单元的形状为正方形网格，但是在实际野外调查过程中，样本的可达性、样本空间位置以及样本边界等因素对抽样精度、抽样效率以及土壤侵蚀模数的计算都会产生影响。因此，以小流域为抽样单元对水土流失面积进行空间抽样设计有待于进一步研究。

2)邹丛荣等[14]在沂蒙山区蒙阴县探索了1%和4%抽样密度下，不同推算方法对县域土壤侵蚀状况监测的影响，发现抽样密度能够对单元直接外推法和单元插值外推法的精度造成较大的影响，对栅格计算法影响较小，但是，该研究未考虑土壤侵蚀分布具有空间自相关性。本研究则探究了不同尺度下土壤侵蚀的空间自相关性，并采用3种空间抽样方法和普通克里金法对土壤侵蚀状况进行预测，比较发现400 m×400 m的土壤侵蚀空间自相关性较弱，且不同抽样方法对县域尺度土壤侵蚀状况估计结果差异较大，空间分层抽样的预测结果明显优于空间随机抽样和空间系统抽样结果，分析其原因发现空间分层抽样中层的划分较为合理，层间差异大，层内差异小，样本布局更加合理。

3)本研究中，以确定性方法(反距离权重法和样条函数法)计算水土流失面积明显高于栅格计算法结果，相对误差为74.79%和25.04%，地统计方法(普通克里金法和协同克里金法)计算水土流失面积相对误差较小(8.25%和9.85%)。这可能是因为研究区土壤侵蚀存在空间异质性，确定性方法在计算过程中未考虑方向，而地统计方法能够通过实测值建立有效且合理变异函数模型，分析空间自相关性的变化情况[24]。在地统计插值中，协同克里金法与普通克里金法计算结果相近，但前者测算的水土流失面积相对误差较高，尤其是对中度侵蚀面积的计算结果明显高于栅格计法，主要原因可能是受单一辅助因素坡度的影响。研究表明，土壤侵蚀与多种因素有关，受地形、土壤及降雨等因素交互作用的影响，侵蚀过程比较复杂[31]。因此，可以对区域土壤侵蚀因子敏感性进行分析，以助于选择合适的辅助变量，从而提高空间预测精度。