基于坡度的黑土區切溝密度協同克里格插值方法研究

王 平, 李 浩, 陳 帥, 徐金忠, 張興義

(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