基于普查样点的贵州省潜在水土流失空间插值方法比较

Comparison of Spatial Interpolation Methods for Potential Soil Erosion in Guizhou Province Based on Census Sampling Units

摘要

目的 为评估贵州省潜在水土流失状况并支撑石漠化与水土流失耦合研究，旨在筛选适用于该区域大尺度、低密度采样条件下的最优空间插值方法。 方法 基于全国第一次水利普查1 097个抽样单元，结合气候、土壤、地形等多源辅助环境变量，系统比较反距离加权法（IDW）、径向基函数（RBF）、多元线性回归（MLR）、普通克里金（OK）、回归克里金的B模型和C模型（RK-B和RK-C）及地理加权回归（GWRK）7种空间插值方法与地图代数法的差异，并通过独立验证和交叉验证法综合评价各方法性能。 结果 1）贵州省潜在水土流失量平均值为9 992.81 t/km²，空间分布呈“南高北低”趋势。主成分分析提取6个主成分，累积方差解释率达85.46%。2）地图代数法呈破碎化斑块，存在显著高估倾向；IDW、RBF和OK空间分布存在“牛眼”特征；RK-B和GWRK存在过度平滑和压缩低值情况；MLR和RK-C方法空间分布均衡性最优。3）独立检验中，RMSE排序为GWRK<RK-B<RK-C<MLR<IDW<RBF<OK，其中MLR的RI为0.30%，且RK-C的MAE最低为6 634.18 t/km²；交叉检验中，RMSE排序为IDW<OK<RK-C<RK-B<MLR<GWRK<RBF，其中RK-C在RMSE、ME和RI均优于RK-B和MLR。 结论 综合精度、空间分布合理性，RK-C被认为是大尺度、低密度条件下贵州省潜在水土流失空间插值的最优方法。

Abstract

Objective This study aims to assess the potential soil erosion status in Guizhou Province， China， and to support research on the coupling mechanisms between rocky desertification and soil erosion. The objective is to identify the optimal spatial interpolation method suitable for large-scale， low-density sampling conditions in this region. Methods Based on 1 097 sampling units from the First National Census for Water of China， and incorporating multiple auxiliary environmental variables such as climate， soil， and topography， the study systematically compared seven spatial interpolation approaches， including inverse distance weighting （IDW）， radial basis functions （RBF）， multiple linear regression （MLR）， ordinary kriging （OK）， regression-kriging B and C models （RK-B and RK-C）， and geographically weighted regression kriging （GWRK）， with the map algebra method. The performance of each method was comprehensively evaluated using independent validation and cross-validation. Results 1） The average potential soil erosion modulus in Guizhou Province was 9 992.81 t/km²， exhibiting a spatial pattern of higher in the south and lower in the north. Principal component analysis （PCA） extracted 6 principal components， with a cumulative variance explanation rate reaching 85.46%. 2） The map algebra method showed fragmented patches and exhibited a significant overestimation tendency. IDW， RBF， and OK exhibited "bull's-eye" effects in their spatial distribution. RK-B and GWRK demonstrated over-smoothing and low-value compression. MLR and RK-C methods achieved the most balanced spatial distributions. 3） In independent validation， the RMSE ranking was GWRK<RK-B<RK-C<MLR<IDW<RBF<OK. Among these， the RI of MLR was 0.30%， and RK-C achieved the lowest MAE of 6 634.18 t/km². In cross-validation， the RMSE ranking was IDW<OK<RK-C<RK-B<MLR<GWRK<RBF， where RK-C outperformed RK-B and MLR in terms of RMSE， ME， and RI. Conclusion Considering both accuracy and the rationality of spatial distribution， RK-C is identified as the optimal spatial interpolation method for estimating potential soil erosion under large-scale， low-density conditions in Guizhou Province.

Graphical abstract

关键词

贵州省 / 潜在水土流失 / 空间插值 / 低密度 / 普查采样单元

Key words

Guizhou Province / potential soil erosion / spatial interpolation / low-density / census sampling units

我国是世界上水土流失最为严重的国家之一。处于喀斯特中心地带的贵州省，因其较薄的土层厚度和持续的人为扰动，水土流失的危害更为严重。根据2024年中国水土保持公报^［1］公布结果，贵州省共有水土流失面积44 344 km²，占土地总面积的25.17%。持续发生的水土流失，导致土地石漠化。石漠化是指在亚热带脆弱的喀斯特环境背景下，受人类不合理社会经济活动的干扰破坏所造成的土壤严重侵蚀，基岩大面积出露，土地生产力严重下降，地表出现类似荒漠景观的土地退化过程^［2］。我国西南岩溶地区是喀斯特石漠化分布较为广泛且严重的区域，贵州省更是因其独特的地形、自然因素等原因，石漠化面积达247.0万hm²，居岩溶地区首位^［3］。

受易溶性石灰岩的控制，岩溶地区土壤层普遍较薄。降雨引起的水土流失持续发生的结果更容易导致土地石漠化。近年来，越来越多的学者致力于探究二者之间的耦合关系，以期进一步优化脆弱生态环境的治理方案。现有研究^［4-5］多聚焦于流域或区县等小尺度范围，主要利用遥感解译数据结合地图代数法分析二者的时空演变及相关性，揭示如石漠化与水土流失之间存在滞后性或不同岩性背景下的关联性等重要规律。

然而，现有方法针对省以上大尺度的水土流失和石漠化数据，多依赖遥感数据解译结合地图代数方法，其精度易受到植被、水保措施等人为扰动因素影响，计算误差较大。相较而言，采用抽样调查方法能够获取更准确的调查单元水力侵蚀因子状况，有效避免遥感反演的系统误差，特别是在复杂地形区域，显著提高土壤水力侵蚀因子的计算精度^［6-7］。杨志成等^［8］、魏梦瑶等^［9］研究证实利用普查样点数据进行区域土壤侵蚀分布分析的可行性。为深入探究贵州省潜在水土流失与石漠化耦合关系，本研究聚焦于自然因素主导下的潜在水土流失，认为其与石漠化强度的耦合度能更直观地反映二者的关联性。若二者空间格局不耦合，则提示水土流失并非石漠化的唯一诱因。鉴于此，本文旨在通过比较不同插值方法的性能，力图找到基于抽样调查数据的贵州省潜在水土流失空间插值的最优方法。

在区域尺度制图中，地统计插值方法因其在处理空间数据上的优势而被广泛应用。李豪等^［10］、孙京禄等^［11］均基于不同调查对象布设抽样单元获取样本数据的基础上，对比不同空间插值方法所得结果的精度找到适宜各研究对象的区域性空间插值方法。鉴于此，本研究以贵州省为例，旨在系统比较不同空间插值方法在区域尺度潜在水土流失制图中的适用性。选取反距离加权法（inverse distance weighting， IDW）、径向基函数（radial basis functions， RBF）、多元线性回归（multiple linear regression， MLR）、普通克里金（ordinary kriging， OK）、回归克里金（regression kriging）的B模型和C模型（RK-B和RK-C）及地理加权回归（geographic weighted regression， GWRK）7种插值方法，并结合地图代数法共同进行对比。同时为提升模型的预测精度，本研究在第一次全国水利普查水土保持专项调查的1 097个抽样单元的基础上结合文献总结分析^［12-13］，引入11个与潜在水土流失过程密切相关的辅助环境变量，包括高程（digital elevation model， DEM）、坡长坡度因子（slope length and steepness factor， LS因子）、降雨侵蚀力因子（rainfall erosivity factor， R因子）、土壤可蚀性因子（soil erodibility factor， K因子）、归一化植被指数（normalized difference vegetation index， NDVI）、土体厚度、土壤体积质量、土壤有机碳含量、温度、砾石质量分数和降雨量；最后采用独立验证和交叉验证法综合评价各方法的插值精度，以筛选出最适宜于贵州省潜在水土流失空间插值的方法，为后续贵州省潜在水土流失和喀斯特石漠化耦合关系研究奠定基础。

1 材料与方法

1.1 研究区概况

贵州省（24°37′~29°13′N，103°36′~109°35′E）地处中国西南地区的东南部，全省总面积17.62万km²，常住人口约3 860万人，属于亚热带季风性气候，年平均降水量1 100 mm，年平均气温16 ℃。地势西高东低，以高原、山地和丘陵为主，岩溶地貌广布（图1）。土壤类型多为黄壤、石灰土等，植被覆盖度较低，在人类活动加剧和生态环境恶化的背景下水土流失和石漠化较为严重。

1.2 数据与预处理

本研究数据主要包括因变量（潜在水土流失量）和自变量（辅助环境变量）。

1）潜在水土流失量计算。核心数据来源于2010—2012年全国第一次水利普查水土保持专项调查在贵州省布设的1 097个抽样单元。潜在水土流失是指只考虑人为不能调控的自然因素情景下发生的水土流失^［14］。潜在水土流失量依据CSLE（Chinese soil loss equation）模型计算，即不考虑其中植被覆盖、耕作、工程等措施，以反映自然因素主导下的侵蚀潜力，具体计算因子及方法参照相关研究^［6，15］。每个抽样单元（不规则多边形）最终产出1个潜在水土流失量。为进行空间插值，提取每个调查单元的几何中心点作为样点，并将该单元潜在水土流失量赋予该样点。按8∶2的比例将1 097个样点随机划分为训练集（878个）和独立验证集（219个）（图1）。

2）辅助环境变量。为提高在低密度样点条件下的空间插值精度，并增强模型对地理过程的解释能力，本研究引入11个与潜在水土流失密切相关的辅助环境变量（表1）。以普查数据为核心时间基准，地形、土壤等相对稳定的因子，用该时期前后可获取的最佳数据集；对于气候、植被等动态因子，则选用最接近研究时段的年份或多年平均值，以代表研究区在目标时段的平均状况。为避免对较粗分辨率数据进行升采样产生伪精度，本研究采用双线性内插法对所有辅助环境变量重采样至30 m，坐标系统一为WGS_1984_Albers。

本研究对研究区抽样单元样点潜在水土流失量与11个辅助环境变量进行相关性分析发现，潜在水土流失量与K因子、LS因子、R因子、温度和降雨均呈正相关；与高程呈显著负相关。不同辅助环境变量之间也存在显著相关，如高程和LS因子呈显著负相关；K因子和土体厚度呈显著正相关。为避免辅助环境变量间存在多重共线性，对11个辅助环境变量标准化后进行主成分分析（principal component analysis， PCA），并选取累计方差贡献率超过85%的前6个主成分（PC1-PC6）作为最终参与回归类插值模型的自变量。

3）抽样单元几何中心点提取。提取各抽样单元中心点潜在水土流失量时，若该点非Nodata值则直接赋值，若该点是Nodata值则取选择该点9×9（像元）领域进行欧氏距离加权取均值，辅助环境变量值提取方法同上。本研究采用“样本平均值3倍的标准差”来界定特异值，在SPSS中进行检查并用正常值的最大值替代特异值。

4）K因子值修正。本研究抽样单元中土壤类型主要为黄壤和石灰土，二者性状相近，普查数据调查单元的K因子平均值为0.004 584（t·hm²·h）/（hm²·MJ·mm），而刘宝元等^［19］在国家地球系统科学数据中心发布的数据得到的贵州省K因子值平均值为0.004 607 54（t·hm²·h）/（hm²·MJ·mm），2个K值相近，故取2个预测值平均值为0.004 595 5（t·hm²·h）/（hm²·MJ·mm）。然而，宏观预测值可能无法准确反映研究区的局部特性。为提升参数本地化精度，本研究采纳KE等^［25］在贵州省毕节市通过径流小区观测获得的黄壤的K因子实测值0.002 9（t·hm²·h）/（hm²·MJ·mm）作为更可靠的基准。因此，本研究K因子对KE等^［25］的实测数据进行修正，修正系数为0.630 983。故本研究抽样单元K因子值均乘以系数0.630 983进行修正。

5）0值数据处理。本研究抽样单元中存在近300个0值，该值代表该点不存在水土流失。统计效果具有偏态性，直接影响数据正态性与后续克里金系列方法的适用性。借鉴CASTRO等^［26］在“复杂地形日降水空间插值”中的处理思路，本研究首先剔除0值样本，并对剩余非0值进行平方根转换以服从正态分布，即转换后偏度与峰度均接近0，同时Kolmogorov-Smirnov（KS检验）显著性>0.05。随后，利用二元逻辑回归构建0值发生概率图，即若像元值为1则代表有侵蚀，值为0代表无侵蚀。用各插值方法得到的非0表面逐像元乘以与该概率图，从而把0值空间信息回代至最终结果。

1.3 插值方法介绍

本研究选用常见且具有代表性的反距离加权法（IDW）、径向基函数（RBF）、多元线性回归（MLR）、普通克里金（OK）、回归克里金的B模型和C模型（RK-B和RK-C）及地理加权回归（GWRK）7种空间插值方法，其中IDW、RBF、MLR和OK较为常见，不作赘述，在ArcGis10.2软件中直接插值即可。

RK-B法是将预测项进行插值后与残差项用OK插值的结果进行相加，得到RK-B的插值结果^［27］；而RK-C则是用多元线性回归的线性方程拟合预测值的趋势面，再与残差数据的OK插值结果相加，得到RK-C的插值结果^［28］；GWRK则是将RK-B中预测值的拟合换成GWR局部拟合，后续操作与RK-B方法一致^［29］。

OK、RK-B、RK-C和GWRK中涉及的Kriging插值均需保证数据服从正态分布，再在GS⁺9.0进行半方差分析找到最优半方差模型参数，最后在ArcGis 10.2软件中进行插值；其中GWRK中的GWR局部拟合在GWR4中进行，而OK、RK-B和RK-C在ArcGis10.2软件中拟合趋势面。同时为了比较不同插值方法和地图代数方法插值的性能差异，地图代数法的结果用表1中获取到的贵州省R因子、K因子和LS因子于ArcGis 10.2软件中使用栅格计算器直接计算。

1.4 检验精度方法

本研究选取独立验证和交叉验证法中平均误差（mean error， ME）、平均绝对误差（mean absolute error， MAE）、均方根误差（root mean square error， RMSE）和相关系数（correlation coefficient， R）4个指标评价插值方法性能，其中ME、MAE、RMSE的值均越小，R值越接近1，模型性能越好。

考虑到数据量和计算速度，本研究采用K=5的交叉验证方法，即将数据集随机划分为5份，每次选取其中1份作为测试集，其余4份作为训练集；在训练集上完成模型构建后，对测试集进行预测，重复此过程5次，最终以5次各指标计算结果的平均值作为该方法交叉验证的精度指标。

鉴于已有研究^{［27，30］}指出，大尺度、低密度采样下RK-C精度普遍优于其他方法，本文以RK-C为基准，计算相对提高指数（relative improvement， RI），直观对比其他方法与RK-C的性能差异，计算公式为：

\mathrm{R}\mathrm{I}=\frac{\mathrm{R}\mathrm{M}\mathrm{S}{\mathrm{E}}_i^{\text{'}}-\mathrm{R}\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{R}\mathrm{K}-\mathrm{C}}^{\text{'}}}{\mathrm{R}\mathrm{M}\mathrm{S}{\mathrm{E}}_i^{\text{'}}}\times \mathrm{100}\mathrm{\%}（1）

式中：\mathrm{R}\mathrm{M}\mathrm{S}{\mathrm{E}}_i^{\text{'}}为第i个方法对应的RMSE，\mathrm{R}\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{R}\mathrm{K}-\mathrm{C}}^{\text{'}}为RK-C方法下的RMSE；i为IDW、RBF、MLR、OK、RK-B、GWRK和地图代数方法。

RI值为正且越大，说明RK-C方法与第i个方法相比精度提高更显著，RI为负，则说明RK-C方法不如与之相比方法精度提高的程度。

2 结果与分析

2.1 潜在水土流失量描述性统计分析

对研究区内训练集的潜在水土流失量进行描述性统计见表2。从统计结果可以看出，原始训练集（868个样点）潜在水土流失量为0~39 562.74 t/km²，平均值为9 992.81 t/km²，标准差为8 630.91 t/km²，数据呈现正偏态。剔除0值并经平方根转换后，数据服从正态分布。训练集（701个样点）潜在水土流失量为22.77~195.35 t/km²，平均值为105.14 t/km²，标准差为36.34 t/km²，适用于后续克里金系列插值。

2.2 辅助环境变量主成分分析

为消除辅助环境变量间的多重共线性，对11个辅助环境变量进行主成分分析（PCA）。由表3可知，前6个主成分累积方差解释率为85.46%，有效压缩原始数据信息。因此，本研究选取这6个主成分作为后续回归克里金等模型的核心自变量，以提升模型的稳健性和解释能力。

2.3 不同插值方法的插值结果分析

由图2可知，总体上所有结果的空间分布均呈现西南和东南值高、西北值低的“南高北低”趋势。南部黔南州、黔中安顺市及黔西南部分区域为高值集中区，北部遵义市、毕节市大部则以低值区为主，中等侵蚀区在中部地区形成连续过渡带，与南部喀斯特山区地形破碎、降雨集中且土壤抗蚀性弱，北部高原丘陵区地形相对平缓、降雨强度较低的自然背景一致。然而，不同方法在空间细节的表达和数值分布的合理性上表现出显著差异。

地图代数法的结果在空间上呈现破碎化的斑块，其数值分布极不均匀。由表4可知，最大值达758 717.00 t/km²，平均值38 689.10 t/km²，极端值占比高达0.74%，与其他插值方法相比，整体数值分布较大且未能形成清晰的地理分异规律。IDW、RBF和OK作为不依赖辅助环境变量的插值方法，其结果在样点稀疏区域出现明显的“牛眼”特征，即围绕少数高值样点生成不自然的独立高值圈，正是由于空间平滑性不足使空间变异被过度夸大。而IDW插值结果的标准差高达3 840.11 t/km²，与区域侵蚀过程的连续性不符。

引入辅助环境变量后的回归类模型在一定程度上克服上述问题，但各自也存在局限。RK-B和GWRK方法存在明显的过度平滑和低值压缩问题。RK-B和GWRK的标准差显著偏小分别为1 807.03、1 831.42 t/km²，且预测的最小潜在水土流失量分别高达4 093.23、4 041.38 t/km²，与样点数据实际存在大量0值和低值区的现实情况严重不符，表明2种方法未能准确表现出区域内存在侵蚀较弱的区域。

相比之下，MLR和RK-C的插值结果在空间分布上表现出更好的均衡性，不仅全局趋势一致，还在细节上与地理环境相契合，即南边的高值区域与该地喀斯特峰林地貌及高强度降雨相匹配，而北边的低值区与该地乌蒙山区相较周边较为平坦的地势和较少的降雨相符合。同时RK-C的最大值和最小值分别为30 683.80 、839.44 t/km²；MLR的最大值和最小值分别为29 943.50、790.84 t/km²，极端值面积占比均仅0.02%，没有产生虚假的局部热点避免“牛眼”的特征出现，更重要的是克服了RK-B和GWRK模型的低值压缩缺陷，其数值范围能够更真实地反映从无侵蚀到强烈侵蚀的完整过渡。因此，从空间和数值的分布来看，MLR和RK-C差异不显著，二者均能有效刻画研究区潜在水土流失量的宏观格局。

2.4 不同插值方法预测精度分析

为定量评价各方法的预测性能，通过独立验证法用219个测试集数据与该点预测值进行比较，各方法指标值见表5。首先，地图代数法与其他空间插值法的精度差异巨大，RMSE高达55 861.27 t/km²，且ME为32 842.70 t/km²，显示出极为严重的系统性高估倾向，其预测结果基本不具备可靠性。

其他插值方法的RMSE排序为GWRK<RK-B<RK-C<MLR<IDW<RBF<OK，MAE也呈相似规律，即引入辅助环境变量方法的预测精度普遍高于未引入的，说明结合地理环境变量能够显著提升模型对潜在水土流失空间异质性的解释能力。

在引入辅助环境变量的4个模型中来看，ME均为负值，存在轻微的低估倾向，但与地图代数法的巨大偏差相比可以忽略。RK-C的MAE为6 634.18 t/km²，是所有方法中的最小值，表明RK-C的平均预测误差最小。尽管其RMSE略高于GWRK和RK-B，且RI略差于RK-B和GWRK，但考虑到差异极小，而MAE对异常值不敏感，因此RK-C的整体预测稳健性更优；同时其R值与MLR均为0.17，高于RK-B和GWRK，表明其预测值与测试集实际值的线性相关性更优，故从独立验证检验精度来看RK-C方法最优。

交叉验证结果进一步揭示方法的稳定性差异见表6。RMSE的排序为IDW<OK<RK-C<RK-B<MLR<GWRK< RBF，MAE的排序也类似。其中RBF的RMSE高达14 001.59 t/km²，且R值为0.09，模型预测值与实测值几乎无线性相关性，稳定性较差；IDW和OK在交叉验证和独立验证结果中的RMSE表现差异显著，说明二者性能受数据集划分的影响，泛化能力不足。余下引入辅助环境变量的方法中，GWRK的RMSE高达13 105.24 t/km²，R为0.03，与RBF类似，存在严重过拟合问题，实际运用价值较低。

在核心对比MLR、RK-B和RK-C 3个方法时，RK-C优势显著。RK-C的RMSE和ME较之RK-B和MR均较小，尽管三者均存在对潜在水土流失量的低估，但RK-C的ME相较于其他2个偏差较小；MLR和RK-B的RI分别为0.88%和0.62%，也体现出RK-C更优的预测能力，因此在交叉验证方法下，RK-C是数据层面表现最优的方法。

综上所述，引入辅助环境变量的回归类模型在独立验证中，预测精度普遍高于未引入辅助环境变量的方法。在回归类方法内部比较中，尽管RK-B和GWRK的RMSE值最小且RI为负数，但其空间分布中存在明显的低值压缩和过度平滑问题，导致空间格局失真。而MLR的RI和MAE略逊于RK-C，故从独立验证中RK-C法最优。交叉验证中RMSE存在的显著差异进一步显示，IDW和OK性能更依赖于数据集的划分，泛化能力较弱，且RK-C相较于其他方法在各个指标中均有优势，尽管相较于MLR和RK-B差异较小，但结合空间分布合理性的前提下，本研究选择RK-C作为贵州省大尺度、低密度条件下潜在水土流失空间插值的最优方法。

3 讨 论

3.1 主成分对潜在水土流失量的影响分析

本研究根据主成分分析得到的6个主成分进行线性回归，得到回归方程Y=\mathrm{103.834}-\mathrm{4.637}\mathrm{P}\mathrm{C}\mathrm{1}-\mathrm{3.573}\mathrm{P}\mathrm{C}\mathrm{2}+\mathrm{1.206}\mathrm{P}\mathrm{C}\mathrm{3}+\mathrm{3.766}\mathrm{P}\mathrm{C}\mathrm{4}+\mathrm{12.29}\mathrm{P}\mathrm{C}\mathrm{5}+\mathrm{1.306}\mathrm{P}\mathrm{C}\mathrm{6}（R²=0.172），参考SUMFLETH等^［31］的研究，该R²值处于合理范围。从回归方程可以看出，PC5对潜在水土流失量影响最显著，且完全由K因子组成，表明土壤自身的抗侵蚀能力是决定贵州省潜在水土流失强度的首要因素。贵州省广泛分布的黄壤和石灰土，其土壤结构相对薄弱，抗蚀性差，因此K因子的微小变化显著加剧侵蚀潜力；PC1、PC2负向作用明显，土壤有机碳提升可增强团聚性，厚土体、高砾石能削弱侵蚀动力。高海拔区域，温度较低，生物活动减缓，有机质分解较慢，有利于有机碳的积累，从而间接增强土壤的抗侵蚀能力；PC4次之，温度、R因子及LS因子协同，驱动侵蚀势能累积；PC3与PC6影响相对微弱，区域降雨格局与植被覆盖的侵蚀调控效应，受其他因素叠加影响。综合来看，贵州省的潜在水土流失空间分异，在宏观上是由土壤内在的抗蚀性和土层、砾石所主导，并受到有机碳的显著调控；而侵蚀的驱动力则主要来源于降雨和地形。

3.2 区域尺度潜在水土流失图插值方法的选择

区域水土流失空间表达的常见方法是地图代数法和空间插值法，而地图代数法因其计算方便而经常被用于水土流失图的制作方法，但其受限于遥感数据解译的精度，较难得到与工程措施精度10 m相匹配的解译精度，使其计算的精度大大降低^［32］。而空间插值法获取到的抽样单元具体信息，数据本身的精度很高，但可能在插值过程中受限于抽样单元的密度和插值方法的选择等因素而产生不同的效果。本研究通过比较在贵州省区域尺度下7种不同空间插值方法和地图代数法时发现，地图代数方法RMSE为55 861.27 t/km²，显著高于其他插值方法的RMSE，即在区域尺度下地图代数得到的结果在精度上远不及空间插值法的结果，与朱浩楠等^［33］、DONG等^［34］的研究结果一致。

本研究发现，引入辅助环境变量插值方法的预测精度普遍高于未引入的方法，即使本文所选的辅助环境变量数量有限，且与潜在水土流失量的相关性不高（回归模型R²=0.172），但其为插值方法提供关键的地理约束，从而显著提高模型在低密度样点条件下的预测稳健性，与CARRERA-HERNANDEZ等^［35］得到的观点一致。其中RK-C和MLR 2个方法得到的空间分布图趋势大体一致，RK-C方法较之MLR的差异在于对MLR拟合得到的回归残差进行OK插值，旨在捕捉MLR模型未能解释的空间自相关信息。然而，本研究中MLR方法相较于RK-C不论是空间分布还是精度检验的数值上，都不存在显著差异，主要归因于当前大尺度、低密度采样条件的限制。RK-C的优势高度依赖于对其残差进行有效的空间结构拟合，即构建稳定的半方差函数模型，但需要足够密集的样点来捕捉短程变异。辅助环境变量对潜在水土流失量的全局解释力有限（R²=0.172），意味着MLR模型留下大量未被解释的方差。此时稀疏的样点难以支撑其空间自相关结构的稳健拟合，因此克里金修正的有效性受限，从而掩盖RK-C本应具有较MLR的精度优势，与张小美等^［36］的研究结果一致。尽管如此，RK-C在精度指标上仍表现出最稳定的综合性能，如独立验证中RK-C的MAE和RI均略胜MLR，交叉验证中RK-C的RMSE和ME也略优于MLR，且其理论框架在引入辅助变量的基础上，进一步考虑空间自相关性。因此，本研究选择RK-C作为大尺度、低密度下空间插值的最优方法，与陈昌华等^［30］在东北进行土壤全氮空间插值方法选择得到的结果一致。但需注意的是，RK-C方法的最优性并非普适于所有大尺度、低密度地形场景，其适用性高度依赖地形特征与辅助环境变量的匹配度，以及对空间变异规律的契合性。

在高原地区，若地形条件较为特殊，且局部气候存在差异，传统RK-C的全局回归则难以捕捉局部趋势，如聂磊等^［37］在四川省复杂地形的月平均气温插值中发现，混合地理加权回归克里金（（mixed） geographically weighted regression，（m） GWRK）通过将海拔设为全局变量，坡向设为局部变量，此时（m） GWRK法在四川西北高原的RMSE比传统RK-C降低33.7%，而RK-C无法适应局部非平稳性。同时在丘陵地区也存在，CHHIPA等^［38］在印度丘陵区土壤pH插值中发现，若目标变量与地形因子存在强相关性，RK-C仍能保持优势，且边缘区域方差更均匀；但GIA等^［39］在越南丘陵区插值土壤有机碳时发现，当辅助环境变量与目标变量相关性弱时，RK-C的精度略低于OK，此时OK因无需依赖辅助环境变量，计算更简洁，更具有优势，表明大尺度、低密度下RK-C是否为最优插值方法，影响因素很多。一方面是辅助环境变量与目标变量的相关性，相关性越强，RK-C的趋势面拟合越精准，残差修正的空间越小，优势越显著；另一方面，还需要考虑地形的空间平稳性。若地形整体趋势清晰，则RK-C可以捕捉趋势后发挥其优势；但若局部变异剧烈，GWRK能够通过局部回归更好地掌握不同区域的局部特征从而提高预测精度。如果辅助环境变量相关性不足，但抽样单元数量充足、分布均匀的情况下，选择OK、MLR等操作更简洁的方法也不失为一种好的选择。因此，不论是哪种空间插值方法，其适用性均需适配特定的场景。

3.3 不足与展望

本研究所用普查抽样单元相对于贵州省17.62万km²广阔的面积而言密度偏低，且分布并非完全均匀，这种稀疏的采样条件不仅限制插值方法的精度，更制约RK-C等混合方法理论优势的发挥，即稀疏采样点难以支撑对复杂残差场进行稳健的半方差函数拟合。数据源的局限性会进一步加剧模型的不确定性。一方面部分辅助环境变量的获取年份与2010—2012年的核心调查时段未能完全同步，这种时间上的不匹配可能引入误差；另一方面，辅助环境变量与潜在水土流失量的相关性偏低，不仅削弱MLR趋势面的预测能力，也导致RK-C需要拟合一个方差更大、结构更复杂的残差场，在低密度采样条件下尤为困难。

鉴于上述限制，本研究只能基于精度评价指标进行相对优劣的筛选。未来希望可以在考虑不确定性分析的前提下用到随机森林（RF）和随机森林加残差克里格（RFRK）^［40］等混合机器学习方法，结合更高精度、更全面的辅助环境变量，以期更充分地挖掘高维辅助环境变量中的非线性信息，从而进一步提高大尺度、低密度插值的精度。本研究对K因子的修正是基于宏观预测数据与本地实测数据间的尺度差异所采取的必要措施。鉴于研究区以性状相近的黄壤和石灰土为主，当前的统一校正是提升模型本地化精度的合理方案。但贵州省实际土壤类型丰富，未来工作时应获取更多土壤类型的实测数据以实现更精细的参数校正。本研究对于0值的处理选择二元逻辑回归法生成0概率图，尽管该方法优于IDW插值得到的0概率图性能，但R²仅0.006。因此，后续研究亟需借鉴并探索更稳健的0值处理方法^{［26，35］}。

4 结 论

1）基于第一次全国水利普查的1 097个抽样单元数据统计，贵州省潜在水土流失量为0~39 562.74 t/km²，平均值为9 992.81 t/km²，经主成分分析，提取6个主成分，累积方差解释率达85.46%。有效融合气候、土壤、地形等辅助环境变量信息。

2）所有插值结果均揭示贵州省潜在水土流失呈“南高北低”的总体分布格局。南部黔南州、黔中安顺市及黔西南部分区域为高值区，北部遵义市、毕节市大部为低值区，中部为中等侵蚀过渡带。IDW、RBF和OK存在显著局部高值，产生“牛眼”特征；RK-B和GWRK最小值分别为4 093.23、4 041.38 t/km²，较之其余方法值过大，存在过度平滑和压缩低值的情况；相比之下，MLR和RK-C的空间分布最为合理，高值区与喀斯特地貌及强降雨匹配，低值区与北部平缓地形及少雨特征一致。

3）从独立验证精度来看，地图代数方法精度最低RMSE高达55 861.27 t/km²，存在高估倾向。RMSE排序为GWRK<RK-B<RK-C<MLR<IDW<RBF<OK，RK-C的MAE最低为6 634.18 t/km²，且RI略高于MR。从交叉验证精度来看，RMSE排序为IDW<OK<RK-C<RK-B<MLR<GWRK<RBF，IDW和OK较独立验证存在显著差异，说明其性能高度依赖于数据集的划分。RK-C在ME和RI方面均优于MLR和RK-B，故综合考虑精度和空间分布合理性，RK-C被确认为大尺度、低密度条件下贵州省潜在水土流失空间插值的最优方法。

参考文献

[1] 中华人民共和国水利部.中国水土保持公报（2024年）［R/OL］.2025-03-21.

[2] Ministry of Water Resources of the People's Republic of China. 2024 China soil and water conservation bulletin［R］.2025-03-21. http：//www.mwr.gov.cn/sj/tjgb/zgstbcgb/202510/t20251001_2067996.html.

[3] 王世杰.喀斯特石漠化概念演绎及其科学内涵的探讨［J］.中国岩溶，2002，21（2）：31-35.

[4] WANG S J. Concept deduction and its connotation of karst rocky desertification［J］.Carsologica Sinica，2002，21（2）：31-35.

[5] 国家林业和草原局，国家公园管理局.中国·岩溶地区石漠化状况公报［R/OL］.2018-12-14.

[6] National Forestry and Grassland Administration， Natinal Park Administration. Bulletin on the status of rocky desertification in karst areas of China［R/OL］.2018-12-14.

[7] 何茂林，张玉珊，高家勇，等.喀斯特区土壤侵蚀与石漠化协同演变及交互关系［J］.水土保持学报，2023，37（1）：140-150.

[8] HE M L， ZHANG Y S， GAO J Y， et al. Synergistic evolution and interaction of soil erosion and rocky desertification in karst area［J］.Journal of Soil and Water Conservation，2023，37（1）：140-150.

[9] 罗旭玲，白晓永，谭秋，等.不同岩性背景下土壤侵蚀与石漠化关联性分析［J］.生态学报，2018，38（24）：8717-8725.

[10] LUO X L， BAI X Y， TAN Q， et al. Effect of lithology background on the correlation between soil erosion and rock desertification［J］.Acta Ecologica Sinica，2018，38（24）：8717-8725.

[11] 刘宝元，郭索彦，李智广，等.中国水力侵蚀抽样调查［J］.中国水土保持，2013（10）：26-34.

[12] LIU B Y， GUO S Y， LI Z G， et al. China's hydraulic erosion sampling survey［J］.Soil and Water Conservation in China，2013（10）：26-34.

[13] YANG Q K， ZHU M Y， WANG C M， et al. Study on a soil erosion sampling survey in the pan-third pole region based on higher-resolution images［J］.International Soil and Water Conservation Research，2020，8（4）：440-451.

[14] 杨志成，张卓栋，张科利，等.基于普查样点的贵州省耕地土壤侵蚀空间分布［J］.中国水土保持科学，2020，18（2）：62-68.

[15] YANG Z C， ZHANG Z D， ZHANG K L， et al. Spatial distribution of cropland soil erosion in Guizhou Province based on sample data of the first national water census［J］.Science of Soil and Water Conservation，2020，18（2）：62-68.

[16] 魏梦瑶，张卓栋，刘瑛娜，等.基于CSLE模型的广西土壤侵蚀规律［J］.水土保持研究，2020，27（1）：15-20.

[17] WEI M Y， ZHANG Z D， LIU Y N， et al. Characteristics of soil erosion in Guangxi based on CSLE［J］.Research of Soil and Water Conservation，2020，27（1）：15-20.

[18] 李豪，文安邦，刘涛，等.基于¹³⁷Cs本底值计算模型和地理加权回归克里格对¹³⁷Cs本底值空间分布的预测［J］.中国水土保持科学，2018，16（5）：57-66.

[19] LI H， WEN A B， LIU T， et al. Spatial distribution of ¹³⁷Cs reference inventory in Sichuan province using geographically weighted regression Kriging combined with ¹³⁷Cs reference inventory mathematical model［J］.Science of Soil and Water Conservation，2018，16（5）：57-66.

[20] 孙京禄，张贝尔，王淑文，等.基于BP神经网络的陕北黄土高原沟壑密度空间插值研究［J］.中国水土保持，2023（8）：24-27.

[21] SUN J L， ZHANG B E， WANG S W， et al. Research on spatial interpolation of gully density in northern Shaanxi of the Loess Plateau based on the back propagation neural network［J］.Soil and Water Conservation in China，2023（8）：24-27.

[22] 朱青，国佳欣，郭熙，等.基于随机森林算法的土壤侵蚀影响因子研究：以赣江上游流域为例［J］.水土保持通报，2020，40（2）：59-68.

[23] ZHU Q， GUO J X， GUO X， et al. Research on influencing factors of soil erosion based on random forest algorithm：A case study in upper reaches of Ganjiang River basin［J］.Bulletin of Soil and Water Conservation，2020，40（2）：59-68.

[24] ISLAM F， AHMAD M N， JANJUHAH H T， et al. Modelling and mapping of soil erosion susceptibility of murree， sub-Himalayas using GIS and RS-based models［J］.Applied Sciences，2022，12（23）：e12211.

[25] 马晓微，杨勤科.基于GIS的中国潜在水土流失评价指标研究［J］.水土保持通报，2001，21（2）：41-44.

[26] MA X W， YANG Q K. A study on indexes choice and extraction of China potential soil and water loss based on GIS［J］.Bulletin of Soil and Water Conservation，2001，21（2）：41-44.

[27] 刘子玥，章文波，殷兵，等.基于CSLE的泛第三极地区土壤侵蚀强度空间分布特征［J］.水土保持研究，2024，31（3）：1-9.

[28] LIU Z Y， ZHANG W B， YIN B， et al. Spatial distribution pattern of soil erosion intensity in the pan-third pole region based on CSLE［J］.Research of Soil and Water Conservation，2024，31（3）：1-9.

[29] HAWKER L， UHE P， PAULO L， et al. A 30 m global map of elevation with forests and buildings removed［J］. Environmental Research Letters，2022，17（2）：e24016.

[30] SUN Y W， ZHANG H M， YANY Q K， et al. A new high-resolution global topographic factor dataset calculated based on SRTM［J］.Scientific data，2024，11：e101.

[31] 谢云，刘宝元.中国30 m分辨率的降雨侵蚀力图（1981—2010年）［DS/OL］.（2024-05-07）.国家地球系统科学数据中心-黄土高原分中心. DOI：10.11866/db.loess.2015.02 .

[32] XIE Yun， LIU Baoyuan. Graphs of rainfall erosivity with 30 m resolution on the Loess Plateau （1981—2010）. ［DS/OL］.（2024-05-07）.National Earth System Science Data Center，Sub-center of Loess Plateau. DOI：10.11866/db.loess.2015.02 .

[33] 刘宝元，梁音，曹龙熹，等.中国土壤可蚀性因子K 30 m格网数据［DS/OL］.（2021-03-30）.北京师范大学地理科学学部地理数据与应用分析中心.

[34] LIU B Y， LIANG Y， CAO L X， et al. Grid data on soil erodibility factor K 30 meter in China［DS/OL］. （2021-03-30）Center for Geographic Data and Application Analysis， Department of Geographic Sciences， Beijing Normal University.

[35] 高吉喜，史园莉，张宏伟，等.中国区域250 m归一化植被指数数据集（2000—2023）［DS/OL］.（2025-04-27）.国家青藏高原数据中心.

[36] GAO J X， SHI Y L， ZHANG H W， et al. China regional 250 m normalized difference vegetation index data set （2000—2023）. ［DS/OL］.（2025-04-27）. National Tibetan Plateau / Third Pole Environment Data Center.

[37] LIU F， ZHANG G L， SONG X D，et al. High-resolution and three-dimensional mapping of soil texture of China［J］.Geoderma，2020，361：e114061.

[38] 贵州省 250 m分辨率土壤容重数据集［DS/OL］.国家地球系统科学数据中心，黄土高原科学数据中心2023-06-13.

[39] meter resolution soil bulk density dataset in Guizhou Province［DS/OL］.National Earth System Science Data Center， Loess Plateau Science Data Center（2023-06-13）.

[40] 贵州省 250 m分辨率土壤有机碳含量数据集［DS/OL］.国家地球系统科学数据中心，黄土高原科学数据中心2023-06-13.

[41] meter resolution soil organic carbon content dataset in Guizhou Province［DS/OL］.National Earth System Science Data Center， Loess Plateau Science Data Center（2023-06-13）.

[42] PENG S Z， DING Y X， LIU W Z， et al. 1 km monthly temperature and precipitation dataset for China from 1901 to 2017［J］.Earth System Science Data，2019，11：1931-1946.

[43] KE Q H， ZHANG K L. Interaction effects of rainfall and soil factors on runoff， erosion， and their predictions in different geographic regions［J］.Journal of Hydrology，2022，605：e127291.

[44] CASTRO L M， GIRONÁS J， FERNÁNDEZ B. Spatial estimation of daily precipitation in regions with complex relief and scarce data using terrain orientation［J］.Journal of Hydrology，2014，517：481-492.

[45] ODEH I O A， MCBRATNEY A B， CHITTLEBOROUGH D J. Further results on prediction of soil properties from terrain attributes： Heterotopic cokriging and regression-Kriging［J］.Geoderma，1995，67（3）：215-226.

[46] HENGL T， HEUVELINK G B M， STEIN A. A generic framework for spatial prediction of soil variables based on regression-Kriging［J］.Geoderma，2004，120（1/2）：75-93.

[47] FOTHERINGHAM A S， CHARLTON M E， BRUNSDON C. Geographically weighted regression： A natural evolution of the expansion method for spatial data analysis［J］.Environment and Planning A： Economy and Space，1998，30（11）：1905-1927.

[48] 陈昌华，陈锡云，徐英.基于大尺度低密度样点的东北土壤全氮空间插值方法比较［J］.水土保持通报，2014，34（6）：153-161.

[49] CHEN C H， CHEN X Y， XU Y. Comparison of spatial interpolation methods for soil total nitrogen content at large scale using low density soil survey data in northeast China［J］.Bulletin of Soil and Water Conservation，2014，34（6）：153-161.

[50] SUMFLETH K， DUTTMANN R. Prediction of soil property distribution in paddy soil landscapes using terrain data and satellite information as indicators［J］.Ecological Indicators，2008，8（5）：485-501.

[51] TENG H F， HU J， ZHOU Y， et al. Modelling and mapping soil erosion potential in China［J］.Journal of Integrative Agriculture，2019，18（2）：251-264.

[52] 朱浩楠，杨勤科，杜朝正，等.基于CSLE模型的巴基斯坦土壤侵蚀评价［J］.水土保持学报，2021，35（2）：22-30.

[53] ZHU H N， YANG Q K， DU C Z， et al. Assessment of soil erosion in Pakistan based on CSLE model［J］.Journal of Soil and Water Conservation，2021，35（2）：22-30.

[54] DONG L X， FU S H， LIU B Y， et al. Comparison and quantitative assessment of two regional soil erosion survey approaches［J］.International Soil and Water Conservation Research，2023，11（4）：660-668.

[55] CARRERA-HERNÁNDEZ J J， GASKIN S J. Spatio temporal analysis of daily precipitation and temperature in the basin of Mexico［J］.Journal of Hydrology，2007，336（3/4）：231-249.

[56] 张小美，高春瑞，闫晓斌，等.典型山西黄土高原区土壤有机质的空间异质性及空间插值预测［J］.山西农业科学，2023，51（7）：785-792.

[57] ZHANG X M， GAO C R， YAN X B， et al. Spatial heterogeneity and spatial interpolation prediction of soil organic matter in the typical Loess Plateau area of Shanxi Province［J］.Journal of Shanxi Agricultural Sciences，2023，51（7）：785-792.

[58] 聂磊，舒红，刘艳.复杂地形地区月平均气温（混合）地理加权回归克里格插值［J］.武汉大学学报（信息科学版），2018，43（10）：1553-1559.

[59] NIE L， SHU H， LIU Y. Interpolation of monthly average temperature by using （mixed） geographically weighted regression Kriging in the complex terrain region［J］.Geomatics and Information Science of Wuhan University，2018，43（10）：1553-1559.

[60] CHHIPA V， STEIN A， SHANKAR H， et al. Assessing and transferring soil health information in a hilly terrain［J］.Geoderma，2019，343：130-138.

[61] PHAM T GIA， KAPPAS M， VAN HUYNH C， et al. Application of ordinary Kriging and regression Kriging method for soil properties mapping in hilly region of central Vietnam［J］.ISPRS International Journal of Geo-Information，2019，8（3）：e147.

[62] GUO P T， LI M F， LUO W， et al. Digital mapping of soil organic matter for rubber plantation at regional scale： An application of random forest plus residuals Kriging approach［J］.Geoderma，2015，237：49-59.

基金资助

国家自然科学基金重点项目(41730748)