基于坡度的黑土区切沟密度协同克里格插值方法研究

王 平, 李 浩, 陈 帅, 徐金忠, 张兴义

(1.黑龙江省水土保持科学研究所, 哈尔滨 150070;2.中国科学院 东北地理与农业生态研究所, 哈尔滨 150081; 3.中国科学院大学, 北京 100049)

基于坡度的黑土区切沟密度协同克里格插值方法研究

王 平1, 李 浩2,3, 陈 帅2,3, 徐金忠1, 张兴义2

(1.黑龙江省水土保持科学研究所, 哈尔滨 150070;2.中国科学院 东北地理与农业生态研究所, 哈尔滨 150081; 3.中国科学院大学, 北京 100049)

切沟侵蚀已成为东北黑土区土壤侵蚀的重要组成部分。结合野外样区实地调查切沟分布数据与空间插值方法是快速获取大面积切沟密度值的有效手段。流域的平均坡度值是切沟形成的影响因素之一，为提高切沟密度值空间分布的插值精度，在黑龙江省海伦市内，利用网格法均匀选取40个1 km2左右的小流域，在实地测量区域内切沟密度数据的基础上，应用距离权重反比法、普通克立格和以小流域的平均坡度作为协同变量的协同克立格对其做空间插值。结果表明：切沟密度与协同区域化变量受结构性因素的影响远大于随机性因素，均为强空间自相关性；距离权重反比法与普通克里格法的预测精度相近；协同区域化变量的空间结构性优于单一变量，协同克立格法生成的空间分布图精细度明显提高，均方根误差降低20%以上，预测值与实测值的相关系数提高89%以上，协同克里格可有效提高区域切沟插值精度。

黑土区; 切沟密度; 协同克立格; 空间插值; 预测精度

水土流失导致耕层变薄、土壤退化、肥力降低,影响生产力的提高[1-2]。沟蚀在东北黑土区广泛分布，在黑土水土流失中意义重大[3]。沟蚀不仅会加剧水土流失、影响农业生产，而且会引发一系列社会问题[4]。朱显谟[5]认为切沟侵蚀是沟蚀的最后阶段，以不能横过耕作为其主要特征。切沟侵蚀是土地退化的主要过程之一[6]。野外证据表明，切沟侵蚀量可占到土壤侵蚀总量的30%～90%[7]。由此可见，确定切沟侵蚀强度对于开展东北黑土区土壤侵蚀及其防治的研究具有重要的实际意义。

沟蚀强度可以用沟蚀密度来表达[8-9]。目前切沟密度的获取多为依赖遥感、数字地形模型、计算机模拟等间接研究方法。遥感通过对航片或卫片的判读,获取切沟侵蚀的解译指标，其经济、快速、周期短的特点使其得到广泛应用[10]。然而受遥感影像分辨率的限制，尺度低于影像分辨率的切沟很难被解译，而未被评估的切沟可反映潜在沟蚀危害。野外调查发现黑土区沟蚀中相当比例的切沟宽度在5 m以下，而高分辨率影像的应用可提高切沟解译率[11]。数字地形模型利用多期航片生成的DEM估算切沟侵蚀量，效果较好，但难以获得切沟空间分布[12]。目前对切沟的野外观测主要集中在测量切沟的形态参数、发生部位，研究切沟侵蚀的演变过程、侵蚀量等[13]。然而，由于经济和人力的原因，大区域的切沟分布、沟蚀密度等参数的野外实测结果鲜见报道。因此，利用研究区域的采样样区的切沟密度值，通过空间插值来生成目标区域的沟蚀密度分布就成为获取切沟密度分布的一种解决方法。常用的空间插值方法有距离权重法(Inverse Distance)、多项式插值法(Interpolating Polynomials)和克立格法(Kriging)等。在这些方法中，距离权重法最为简便；多项式插值的物理意义不是很明确，容易得到一些难以解释的值。以上两种方法都是基于采样点与临近点的距离进行空间插值的，没有考虑变量的空间结构[14]。地统计学是一种研究地表变量空间变异特征的方法。其中协同克立格法(CoKriging)是对普通克立格法(Ordinary Kriging)的一种扩展应用，它允许使用辅变量对主变量进行预测，考虑了主辅变量之间统计相关性的同时，也考虑了主辅变量之间的互相关性，故较普通克立格法能更准确地预测结果[15-16]。

土壤是一个时空连续的变异体，具有高度的空间异质性[17]，而切沟的形成与发育受到土壤特性、环境及人为因素的共同影响，因此切沟的分布也应当具有一定的空间异质性。前人研究表明，坡度决定侵蚀沟形成的位置和发育的大小。因此，能否通过将坡度信息融入到切沟密度的空间插值过程来提高切沟密度的插值精度是本文的研究目的。本文以黑龙江省海伦市为例，研究了切沟的空间变异性，选取平均坡度作为辅变量，应用距离权重反比法、普通克立格法和协同克立格法空间插值方法对切沟密度做空间插值，并对该三种方法的插值结果进行了验证与选择，绘制了切沟密度空间分布图，旨在为农业生产合理布局，以及综合治理水土流失和土地退化等问题提供一定的理论参考。

1 材料与方法

1.1 研究区概况

研究区域位于黑龙江省海伦市，地理位置为南起46°59′N，北至47°38′N，西起126°16′E，东至127°07′E，地势从东北到西南由低丘陵、高平原、河阶地、河漫滩依次呈阶梯形逐渐降低。境内除少量残丘外，大部分地势为漫川漫岗平原，并有4条主要河流贯穿全境。海伦气候条件属温带大陆性季风气候区，冬季寒冷干燥，夏季高温多雨，雨热同季。极端最高气温37℃；极端最低气温-39.5℃。总面积4.66×103km2，黑土占土地总面积的63.4%。2007年全市各种土地利用类型的占地面积比例分别为：旱田71.2%，水田6.4%，林地8.1%，草地0.9%，水域1.1%，城镇建设用地0.3%，农村建设用地4.8%，未利用土地7.1%[18]。

1.1 数据预处理

本文所论述的切沟均为在坡耕地上受人类活动影响(如耕作)而产生或加剧的，以不能横过耕作为其主要特征的侵蚀沟，因此海伦市东北部的林区或未利用土地区域未包含在本文的研究区域中。研究区域面积4.01×103km2，将研究区域平均分成40块方格，每个方格面积10 km×10 km。在每个方格的中心附近选取面积约为1 km2的小流域作为采样样方，共计40个采样样方。在每个采样样方内，利用水平精度为5 m的手持Garmin60 CSX GPS测量并记录切沟的沟头、沟中、沟形折点及沟尾的点号和经纬度，对切沟赋唯一编号，直至遍历所有切沟。在获取40个采样样方的切沟基本数据后，导入到ArcGIS9.3。应用Et GeoWizard插件，以沟编号为关键词，将每条切沟的沟头、沟中沟形折点及沟尾连接，生成线状切沟，统计每个采样样方内的切沟总长度，并计算其相应的切沟密度。其中，

切沟密度(km/km2)= 采样样方内切沟总长度(km)/

采样样方面积(km2)

(1)

采用研究区域的SRTM数据(分辨率90 m)生成的平均坡度作为协同克立格法的辅变量。将采样样方抽象成点，以样方质心点为坐标，并换算成米制单位。以每个样方的平面坐标X，Y值，以及切沟密度和平均坡度变量为属性，构建Excel文件。

1.2 空间插值方法

前人已对距离权重反比法(Inverse Distance Weighting，IDW)、普通克立格法(Ordinary Kriging，OK)做了详细描述[19]，在此不再赘述。协同克立格法(CoKriging)是对普通克立格法的一种扩展应用。它的优点是当主变量难以获得或者获取代价很高时，协同克立格法采用更易获取或样本分布密度更高的，且与主变量有一定相关性的辅变量对主变量进行预测，从而提高插值精度[20]。研究指出当主辅变量间的互相关性超过0.45的中等相关程度时，协同克立格法的插值结果精度即明显优于普通克立格法[21]。协同克立格法表达式为：

(2)

式中：Z*(x0)——预测点x0处的估计值；Z1(xi)，Z2(xj)——主变量Z1和辅变量Z2的实测值；λi，λj——分配给主变量Z1和Z2的实测值的权重，且∑λ1i=1，∑λ2j=0；n，p——参与x0点估计的主变量Z1和辅变量Z2的实测值的数目。

1.3 检验标准

每一个采样样方的沟蚀密度预测值都用周围采样样方的值来估算，然后计算所有采样样方预测值与实测值的误差，以此来评判插值方法的优劣。本文采用均方根误差(RMSE)和预测值与实测值的相关系数r2用来表征预测的精度[22]。RMSE越小、r2越大，则预测的精度越高。采用协同克立格法预测的RMSE与距离权重反比法、普通克立格法预测的RMSE减少的百分数(分别为RRMSEIDW，RRMSEOK)来表示预测精度的提高程度，用RIDW和ROK来表示协同克立格法对距离权重反比法和普通克立格法的预测值与实测值相关系数的提高程度。

(3)

式中：Z(xi)，Z*(xi)——实测值和预测值。

RRMSE1= (RMSE1-RMSECK)/

RMSE1×100%

(4)

R1=(rCK-r1)/rCK×100%

(5)

式中：RMSECK，RMSE1——协同克立格法和插值方法1预测的均方根误差；rCK和r1——协同克立格法和插值方法1的预测值和实测值之间的相关系数。

数据的经典统计分析在SPSS19.0中进行，地统计学的半方差分析、距离权重反比法与克立格法的插值与检验均在GammaDesignSoftware公司开发的GS+GeostatisticsforEnvironmentalSciences9.0中进行。

2 结果分析

2.1 描述性统计分析

首先对采样样区的沟蚀密度、平均坡度做平均值、标准差、变异系数分布类型等的经典统计分析(见表1)。

表1 切沟密度、平均坡度的经典统计结果

结果可见，切沟密度与平均坡度的分布服从对数正态分布，变异系数分别为90%与68%，有一定的差异，但均属中等程度的变异[23]。因此，首先对切沟密度与平均坡度的原始值做对数变换，使其符合标准正态分布以便于进一步的数据分析。

当主辅变量间的统计相关性超过0.45的中等相关程度时，协同克立格法的插值结果精度即明显优于普通克立格法。对切沟密度与平均坡度做相关性分析(见图1)发现，二者的Person相关性，即r值为0.54，达到了0.01(双侧)的极显著水平，因此适合选择平均坡度作为辅变量，应用协同克立格法对切沟密度进行插值。

图1 切沟密度与平均坡度的回归分析

2.2空间变异特征分析

图2是切沟密度和平均坡度的地统计学半方差函数拟合模型及其参数，及协同区域化变量交互半方差函数的拟合模型及其参数。切沟密度、平均坡度与交叉半方差函数的理论模型均符合球状模型。协同区域化变量可以是正相关，也可以是负相关，这与协同区域化现象发生的具体过程有关，本研究的结果都是正相关的。决定系数r2反映所选模型对半方差函数的拟合程度，协同区域化变量的交互半方差函数对试验变异函数的拟合程度从0.76提高到了0.90，说明交互半方差函数一方面很好的反映了切沟密度的空间结构特性，另一方面可以提高空间插值结果的精度。

图2 半变异函数图及其模型拟合结果

参数理论模型块金值基台值块金值/基台值/%变程/m决定系数残差切沟密度球状模型0.861.7250409370.140.139平均坡度球状模型0.722.0136520000.760.149协同区域化变量球状模型0.020.3633223000.900.024

C0为块金值，表示由采样误差、短距离的变异、随机和固有变异引起的基底效应；C0+C为基台值，即Sill；块金值与基台值的比值[C0/(C0+C)]表示空间自相关性程度[24]。切沟密度、平均坡度的C0/(C0+C)值分别为0.50与0.36，表现为中等强度的空间自相关性，说明获取的切沟密度、平均坡度分布是由结构性因素(如气候、母质、地形、土壤类型等)和随机性因素(如耕作措施、土地利用等各种人为活动)共同作用的结果[25]。协同区域化变量的C0/(C0+C)值为0.056，小于0.10，表明二者具有强烈的空间自相关性，结构性因素远大于随机性因素引起的空间异质性。另一方面，协同区域化变量的C0/(C0+C)比切沟密度的值小，表明辅变量的应用增强了结构性因素造成的互相关作用。

RSS(Residue Sums of Squares)为残差平方和，表示试验模型与理论模型的拟合程度，简称残差。残差最小是回归模型所追求的目标。从残差来看，尽管辅变量平均坡度的残差(0.149)大于主变量切沟密度的残差(0.139)，但协同区域化变量的残差显著(0.024)小于主变量残差(0.139))，说明交互半方差函数的拟合效果要明显优于半方差函数，这也表明协同区域化变量的空间结构性要优于单一变量。

本研究中辅变量平均坡度的C0/(C0+C)值、RSS值均差于主变量切沟密度的值，但协同区域化变量的r2值、C0/(C0+C)值、RSS值均优于主变量，说明了即使在主辅变量的空间结构性中等情况下，如果主辅变量间互相关性达到中等相关程度(相关性0.76)，协同区域化变量的的空间结构性就可优于单一变量。

2.3 切沟密度的空间分布特征分析

图3是距离权重反比法、普通克立格法、协同克立格法对切沟密度进行空间插值形成的空间分布图。可以看出，3种方法得到的切沟密度分布特征是一致的，均为条带状和斑块状格局。总体来说，海伦市切沟密度基本走向是北低南高，西低东高，并在东南部与南部有两块高值区域。西部切沟密度值最低，低于0.6 km/km2，东部与东南部值最高，超过了2.4 km/km2。导致这种空间分布相似性的原因是多方面的。一方面，坡度是决定侵蚀沟形成的位置和发育的大小的主要因素，而海伦市东北为小兴安岭西南麓，坡度较大；西面与松嫩平原接壤，坡度较缓。另一方面，耕作是加速切沟侵蚀的重要人为因素[26]。海伦市耕地的开垦顺序是从南向北，南部地区的开垦时间最长，开垦强度大。

2.4 不同插值方法的评价

尽管三种插值方法反映的切沟密度总体分布特征相近，但协同克立格法得到的结果更为精细，准确。比较图3a，3b，3c可知，距离权重反比法与普通克立格法插值形成的切沟密度空间分布具有较强的相似性，均无密度大于2.7 km/km2的高值区域，且斑块化现象严重。而在协同克立格法得到的结果中，受辅变量的协助，对局部变异细节的描述更为详细，层次感鲜明；得到的信息更为丰富，局部多中心聚集现象更显著。表3列出的是距离权重反比法、普通克立格法与协同克立格法对切沟密度的预测精度的比较。从检验结果来看，距离权重反比法的预测值与实测值之间的相关系数比普通克立格法的稍高，从0.01提高到0.04，但总体而言是很低的。这可能是切沟密度的取样密度不高导致的。二者的均方根误差相近，均为0.94左右。协同克立格法预测所产生的均方根误差明显低于距离权重反比法和普通克立格法，从0.94左右降低到了0.74，减少了20%以上。但预测值与实测值之间的相关系数有了明显的提高，增加了89%以上。总之协同克立格法相对于距离权重反比法与普通克立格法的预测精度有较大的提高，预测的结果明显。

图3 不同插值方法下的切沟密度分布

方法RMSEr2RRMSE/%R/%IDW0.940.0421.389.8OK0.950.0122.297.9COK0.740.42——

3 结论与讨论

本文比较了距离权重反比法、普通克里格法与协同克里格法对切沟密度做空间插值的预测精度，结论如下：

(1) 切沟密度、平均坡度均呈中等变异强度；二者表现为中等程度的统计相关性，达到了极显著的水平，因此适合选取平均坡度作为辅变量，应用协同克立格法对切沟密度进行插值。

(2) 切沟密度、平均坡度与协同区域化变量的半方差函数模型均符合球状模型；切沟密度与协同区域化变量受结构性因素的影响远大于随机性因素，均为强空间自相关性；而平均坡度受结构性因素与随机性因素的影响相近，表现为中等空间自相关性；协同区域化变量表现为正相关，且交互半方差函数的块金值、基台值都有所降低，协同区域化变量的空间结构性要优于单一变量，对试验变异函数的拟合程度有较大的提高。

(3) 距离权重反比法与普通克里格法的预测精度相近，这可能与采样样方的密度不高有关；协同克立格法对于切沟密度的局部变异细节的描述更为详细，更接近切沟密度的真实分布情况；与二者相比，协同克立格法的均方根误差均降低20%强，相关系数提高89%以上，表明利用平均坡度采用协同克立格法可以提高切沟密度的预测精度。

本文结合切沟密度与平均坡度的关系，应用不同空间插值方法对切沟密度做插值，得到了以流域平均坡度为辅助变量，应用协同克里格对切沟密度做空间插值的精度大于距离权重反比法与普通克里格法的结果。将来的工作应当放在影响切沟形成与发展其他的环境因素，如距水系的距离，以及人为影响因素，如耕作措施等的影响。

[1] 唐克丽,史立人,史德明,等.中国水土保持[M].北京:科学出版社,2004.

[2] 张兴义,刘晓冰,隋跃宇,等.人为剥离黑土层对大豆干物质积累及产量的影响[J].大豆科学,2006,25(2):123-126.

[3] 水利部,中国科学院,中国工程院.中国水土流失防治与生态安全:东北黑土区卷[M].北京:科学出版社,2010.

[4] 孟令钦,李勇.东北黑土区沟蚀研究与防治[J].中国水土保持,2009(12):40-42.

[5] 朱显谟.黄土区土壤侵蚀的分类[J].土壤学报,1956,4(2):99-115.

[6] Oostwoud Wijdenes D J, Poesen J, Vandekerckhove L, et al. Spatial distribution of gully head activity and sediment supply along an ephemeral channel in a Mediterranean environment[J]. Catena,2000,39(3):147-167.

[7] Poesen J, Nachtergaele J, Verstraeten G, et al. Gully erosion and environmental change: importance and research needs[J]. Catena,2003,50(2):91-133.

[8] 刘秉正,吴发启.土壤侵蚀[M].陕西:陕西人民出版社,1997.

[9] 中华人民共和国水利部.SL446—2009黑土区水土流失综合防治技术标准[S].北京:中国水利水电出版社,2009.

[10] 赵英时,陈冬梅,李小文,等.遥感应用分析原理与方法[M].北京:科学出版社,2003.

[11] 李浩,张兴义,刘爽,等.典型黑土区村级尺度侵蚀沟演变[J].中国水土保持科学,2012,10(2):21-28.

[12] Martinez-Casasnovas J A, Ramos M C, Ribes-Dasi M. Soil erosion caused by extreme rainfall events: mapping and quantification in agricultural plots from very detailed digital elevation models[J]. Geoderma,2002,105(1):125-140.

[13] Vandekerckhove L, Poesen J, Govers G. Medium-term gully headcut retreat rates in Southeast Spain determined from aerial photographs and ground measurements[J]. Catena,2003,50(2):329-352.

[14] 李 新,程国栋,卢玲.空间内插方法比较[J].地球科学进展,2000,15(3):260-265.

[15] Stein A, Van Dooremolen W, Bouma J, et al. Cokriging point data on moisture deficit[J]. Soil Science Society of America Journal,1988,52(5):1418-1423.

[16] Ersahin S. Comparing ordinary kriging and cokriging to estimate infiltration rate[J]. Soil Science Society of America Journal,2003,67(6):1848-1855.

[17] 王政权.地质统计学及在生态学中的应用[M].北京:科学出版社,1999.

[18] 谢叶伟,刘兆刚,赵军,等.基于RS与GIS的典型黑土区土地利用变化分析：以海伦市为例[J].地理科学,2010(3):428-434.

[19] 林忠辉,莫兴国,李宏轩.中国陆地区域气象要素空间插值[J].地理学报,2002,57(1):47-56.

[20] Goovaerts,P. Geostatistics for Natural Resources Evaluation[M]. New York:Oxford University Press,1997.

[21] Asli M, Marcotte D. Comparison of approaches to spatial estimation in a bivariate context[J]. Mathematical Geology,1995,27(5):641-658.

[22] Goovaerts P. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall[J]. Journal of Hydrology,2000,228(1):113-129.

[23] Zhang X Y, Sui Y Y, Zhang X D, et al. Spatial variability of nutrient properties in black soil of northeast China[J]. Pedosphere,2007,17(1):19-29.

[24] Chien Y J, Lee D Y, Guo H Y, et al. Geostatistical analysis of soil properties of mid-west Taiwan soils[J]. Soil Science,1997,162(4):291-298.

[25] Cambardella C A, Moorman T B, Parkin T B, et al. Field-scale variability of soil properties in central Iowa soils[J]. Soil Science Society of America Journal,1994,58(5):1501-1511.

[26] Zhang S, Zhang X, Huffman T, et al. Soil loss, crop growth, and economic margins under different management systems on a sloping field in the Black soil area of Northeast China[J]. Journal of Sustainable Agriculture,2011,35(3):293-311.

InterpolationofPermanentGullyDensityBasedonSlopeSteepnessinBlackSoilArea

WANG Ping1, LI Hao2,3, CHEN Shuai2,3, XU Jin-zhong1, ZHANG Xing-yi2

(1.HeilongjiangInstituteofSoilandWaterConservationScience,Harbin150070,China; 2.KeyLaboratoryofMollisolsAgroecology,NortheastInstituteofGeographyandAgroecology,ChineseAcademyofSciences,Harbin150081,China; 3.UniversityofChineseAcademyofSciences,Beijing100039,China)

Gully erosion is serious in northeast China. Combining field-measured gully length density and spatial interpolation is an efficient method to identify gully erosion wizard in large area. Steepness is an important factor affecting gully development. As an auxiliary variable whether it could improve the spatial interpolation performance of gully length density was investigated in this research. The spatial variability of permanent gully density was interpolated by Inverse Distance Weighting, Ordinary Kriging and CoKriging with mean slope steepness of the field sample area from the 40 field-measured sampling data in Hailun county, Heilongjiang Province, located in the black soil area, northeastern of China, and their prediction accuracies were compared. The results indicated that the permanent gully density was strong spatial autocorrelation. The permanent gully density and its coregionalized variables were much more affected by structure factors than stochastic factors. Inverse Distance Weighting got similar prediction accuracy with Ordinary Kriging. Compared with Inverse Distance Weighting and Ordinary Kriging, the accuracy of permanent gully density interpolated by CoKriging was much improved, the root-mean-square error decreased more than 20%, and the determination coefficient between the observed and the predicted values increased more than 89%. Hence, CoKriging is a high accuracy method for the permanent gully density interpolation in northeastern China.

black soil area; permanent gully density; CoKriging; interpolation; prediction accuracy

2014-03-14

：2014-03-27

国家自然科学基金(41171230，41071201)

王平(1965—),女,黑龙江省哈尔滨人,学士,高级工程师,主要从事土壤侵蚀机理及水土流失监测技术方面研究。E-mail:wp626588@163.com

张兴义(1966—),男,黑龙江省密山人,博士,研究员,主要从事黑土侵蚀和保护研究。E-mail:zhangxy@iga.ac.cn

P941

：A

：1005-3409(2014)06-0312-06