基于坡度的黑土區(qū)切溝密度協(xié)同克里格插值方法研究

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

(1.黑龍江省水土保持科學(xué)研究所, 哈爾濱 150070;2.中國科學(xué)院 東北地理與農(nóng)業(yè)生態(tài)研究所, 哈爾濱 150081; 3.中國科學(xué)院大學(xué), 北京 100049)

基于坡度的黑土區(qū)切溝密度協(xié)同克里格插值方法研究

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

(1.黑龍江省水土保持科學(xué)研究所, 哈爾濱 150070;2.中國科學(xué)院 東北地理與農(nóng)業(yè)生態(tài)研究所, 哈爾濱 150081; 3.中國科學(xué)院大學(xué), 北京 100049)

切溝侵蝕已成為東北黑土區(qū)土壤侵蝕的重要組成部分。結(jié)合野外樣區(qū)實地調(diào)查切溝分布數(shù)據(jù)與空間插值方法是快速獲取大面積切溝密度值的有效手段。流域的平均坡度值是切溝形成的影響因素之一，為提高切溝密度值空間分布的插值精度，在黑龍江省海倫市內(nèi)，利用網(wǎng)格法均勻選取40個1 km2左右的小流域，在實地測量區(qū)域內(nèi)切溝密度數(shù)據(jù)的基礎(chǔ)上，應(yīng)用距離權(quán)重反比法、普通克立格和以小流域的平均坡度作為協(xié)同變量的協(xié)同克立格對其做空間插值。結(jié)果表明：切溝密度與協(xié)同區(qū)域化變量受結(jié)構(gòu)性因素的影響遠(yuǎn)大于隨機(jī)性因素，均為強(qiáng)空間自相關(guān)性；距離權(quán)重反比法與普通克里格法的預(yù)測精度相近；協(xié)同區(qū)域化變量的空間結(jié)構(gòu)性優(yōu)于單一變量，協(xié)同克立格法生成的空間分布圖精細(xì)度明顯提高，均方根誤差降低20%以上，預(yù)測值與實測值的相關(guān)系數(shù)提高89%以上，協(xié)同克里格可有效提高區(qū)域切溝插值精度。

黑土區(qū); 切溝密度; 協(xié)同克立格; 空間插值; 預(yù)測精度

水土流失導(dǎo)致耕層變薄、土壤退化、肥力降低,影響生產(chǎn)力的提高[1-2]。溝蝕在東北黑土區(qū)廣泛分布，在黑土水土流失中意義重大[3]。溝蝕不僅會加劇水土流失、影響農(nóng)業(yè)生產(chǎn)，而且會引發(fā)一系列社會問題[4]。朱顯謨[5]認(rèn)為切溝侵蝕是溝蝕的最后階段，以不能橫過耕作為其主要特征。切溝侵蝕是土地退化的主要過程之一[6]。野外證據(jù)表明，切溝侵蝕量可占到土壤侵蝕總量的30%～90%[7]。由此可見，確定切溝侵蝕強(qiáng)度對于開展東北黑土區(qū)土壤侵蝕及其防治的研究具有重要的實際意義。

溝蝕強(qiáng)度可以用溝蝕密度來表達(dá)[8-9]。目前切溝密度的獲取多為依賴遙感、數(shù)字地形模型、計算機(jī)模擬等間接研究方法。遙感通過對航片或衛(wèi)片的判讀,獲取切溝侵蝕的解譯指標(biāo)，其經(jīng)濟(jì)、快速、周期短的特點使其得到廣泛應(yīng)用[10]。然而受遙感影像分辨率的限制，尺度低于影像分辨率的切溝很難被解譯，而未被評估的切溝可反映潛在溝蝕危害。野外調(diào)查發(fā)現(xiàn)黑土區(qū)溝蝕中相當(dāng)比例的切溝寬度在5 m以下，而高分辨率影像的應(yīng)用可提高切溝解譯率[11]。數(shù)字地形模型利用多期航片生成的DEM估算切溝侵蝕量，效果較好，但難以獲得切溝空間分布[12]。目前對切溝的野外觀測主要集中在測量切溝的形態(tài)參數(shù)、發(fā)生部位，研究切溝侵蝕的演變過程、侵蝕量等[13]。然而，由于經(jīng)濟(jì)和人力的原因，大區(qū)域的切溝分布、溝蝕密度等參數(shù)的野外實測結(jié)果鮮見報道。因此，利用研究區(qū)域的采樣樣區(qū)的切溝密度值，通過空間插值來生成目標(biāo)區(qū)域的溝蝕密度分布就成為獲取切溝密度分布的一種解決方法。常用的空間插值方法有距離權(quán)重法(Inverse Distance)、多項式插值法(Interpolating Polynomials)和克立格法(Kriging)等。在這些方法中，距離權(quán)重法最為簡便；多項式插值的物理意義不是很明確，容易得到一些難以解釋的值。以上兩種方法都是基于采樣點與臨近點的距離進(jìn)行空間插值的，沒有考慮變量的空間結(jié)構(gòu)[14]。地統(tǒng)計學(xué)是一種研究地表變量空間變異特征的方法。其中協(xié)同克立格法(CoKriging)是對普通克立格法(Ordinary Kriging)的一種擴(kuò)展應(yīng)用，它允許使用輔變量對主變量進(jìn)行預(yù)測，考慮了主輔變量之間統(tǒng)計相關(guān)性的同時，也考慮了主輔變量之間的互相關(guān)性，故較普通克立格法能更準(zhǔn)確地預(yù)測結(jié)果[15-16]。

土壤是一個時空連續(xù)的變異體，具有高度的空間異質(zhì)性[17]，而切溝的形成與發(fā)育受到土壤特性、環(huán)境及人為因素的共同影響，因此切溝的分布也應(yīng)當(dāng)具有一定的空間異質(zhì)性。前人研究表明，坡度決定侵蝕溝形成的位置和發(fā)育的大小。因此，能否通過將坡度信息融入到切溝密度的空間插值過程來提高切溝密度的插值精度是本文的研究目的。本文以黑龍江省海倫市為例，研究了切溝的空間變異性，選取平均坡度作為輔變量，應(yīng)用距離權(quán)重反比法、普通克立格法和協(xié)同克立格法空間插值方法對切溝密度做空間插值，并對該三種方法的插值結(jié)果進(jìn)行了驗證與選擇，繪制了切溝密度空間分布圖，旨在為農(nóng)業(yè)生產(chǎn)合理布局，以及綜合治理水土流失和土地退化等問題提供一定的理論參考。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)域位于黑龍江省海倫市，地理位置為南起46°59′N，北至47°38′N，西起126°16′E，東至127°07′E，地勢從東北到西南由低丘陵、高平原、河階地、河漫灘依次呈階梯形逐漸降低。境內(nèi)除少量殘丘外，大部分地勢為漫川漫崗平原，并有4條主要河流貫穿全境。海倫氣候條件屬溫帶大陸性季風(fēng)氣候區(qū)，冬季寒冷干燥，夏季高溫多雨，雨熱同季。極端最高氣溫37℃；極端最低氣溫-39.5℃?？偯娣e4.66×103km2，黑土占土地總面積的63.4%。2007年全市各種土地利用類型的占地面積比例分別為：旱田71.2%，水田6.4%，林地8.1%，草地0.9%，水域1.1%，城鎮(zhèn)建設(shè)用地0.3%，農(nóng)村建設(shè)用地4.8%，未利用土地7.1%[18]。

1.1 數(shù)據(jù)預(yù)處理

本文所論述的切溝均為在坡耕地上受人類活動影響(如耕作)而產(chǎn)生或加劇的，以不能橫過耕作為其主要特征的侵蝕溝，因此海倫市東北部的林區(qū)或未利用土地區(qū)域未包含在本文的研究區(qū)域中。研究區(qū)域面積4.01×103km2，將研究區(qū)域平均分成40塊方格，每個方格面積10 km×10 km。在每個方格的中心附近選取面積約為1 km2的小流域作為采樣樣方，共計40個采樣樣方。在每個采樣樣方內(nèi)，利用水平精度為5 m的手持Garmin60 CSX GPS測量并記錄切溝的溝頭、溝中、溝形折點及溝尾的點號和經(jīng)緯度，對切溝賦唯一編號，直至遍歷所有切溝。在獲取40個采樣樣方的切溝基本數(shù)據(jù)后，導(dǎo)入到ArcGIS9.3。應(yīng)用Et GeoWizard插件，以溝編號為關(guān)鍵詞，將每條切溝的溝頭、溝中溝形折點及溝尾連接，生成線狀切溝，統(tǒng)計每個采樣樣方內(nèi)的切溝總長度，并計算其相應(yīng)的切溝密度。其中，

切溝密度(km/km2)= 采樣樣方內(nèi)切溝總長度(km)/

采樣樣方面積(km2)

(1)

采用研究區(qū)域的SRTM數(shù)據(jù)(分辨率90 m)生成的平均坡度作為協(xié)同克立格法的輔變量。將采樣樣方抽象成點，以樣方質(zhì)心點為坐標(biāo)，并換算成米制單位。以每個樣方的平面坐標(biāo)X，Y值，以及切溝密度和平均坡度變量為屬性，構(gòu)建Excel文件。

1.2 空間插值方法

前人已對距離權(quán)重反比法(Inverse Distance Weighting，IDW)、普通克立格法(Ordinary Kriging，OK)做了詳細(xì)描述[19]，在此不再贅述。協(xié)同克立格法(CoKriging)是對普通克立格法的一種擴(kuò)展應(yīng)用。它的優(yōu)點是當(dāng)主變量難以獲得或者獲取代價很高時，協(xié)同克立格法采用更易獲取或樣本分布密度更高的，且與主變量有一定相關(guān)性的輔變量對主變量進(jìn)行預(yù)測，從而提高插值精度[20]。研究指出當(dāng)主輔變量間的互相關(guān)性超過0.45的中等相關(guān)程度時，協(xié)同克立格法的插值結(jié)果精度即明顯優(yōu)于普通克立格法[21]。協(xié)同克立格法表達(dá)式為：

(2)

式中：Z*(x0)——預(yù)測點x0處的估計值；Z1(xi)，Z2(xj)——主變量Z1和輔變量Z2的實測值；λi，λj——分配給主變量Z1和Z2的實測值的權(quán)重，且∑λ1i=1，∑λ2j=0；n，p——參與x0點估計的主變量Z1和輔變量Z2的實測值的數(shù)目。

1.3 檢驗標(biāo)準(zhǔn)

每一個采樣樣方的溝蝕密度預(yù)測值都用周圍采樣樣方的值來估算，然后計算所有采樣樣方預(yù)測值與實測值的誤差，以此來評判插值方法的優(yōu)劣。本文采用均方根誤差(RMSE)和預(yù)測值與實測值的相關(guān)系數(shù)r2用來表征預(yù)測的精度[22]。RMSE越小、r2越大，則預(yù)測的精度越高。采用協(xié)同克立格法預(yù)測的RMSE與距離權(quán)重反比法、普通克立格法預(yù)測的RMSE減少的百分?jǐn)?shù)(分別為RRMSEIDW，RRMSEOK)來表示預(yù)測精度的提高程度，用RIDW和ROK來表示協(xié)同克立格法對距離權(quán)重反比法和普通克立格法的預(yù)測值與實測值相關(guān)系數(shù)的提高程度。

(3)

式中：Z(xi)，Z*(xi)——實測值和預(yù)測值。

RRMSE1= (RMSE1-RMSECK)/

RMSE1×100%

(4)

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

(5)

式中：RMSECK，RMSE1——協(xié)同克立格法和插值方法1預(yù)測的均方根誤差；rCK和r1——協(xié)同克立格法和插值方法1的預(yù)測值和實測值之間的相關(guān)系數(shù)。

數(shù)據(jù)的經(jīng)典統(tǒng)計分析在SPSS19.0中進(jìn)行，地統(tǒng)計學(xué)的半方差分析、距離權(quán)重反比法與克立格法的插值與檢驗均在GammaDesignSoftware公司開發(fā)的GS+GeostatisticsforEnvironmentalSciences9.0中進(jìn)行。

2 結(jié)果分析

2.1 描述性統(tǒng)計分析

首先對采樣樣區(qū)的溝蝕密度、平均坡度做平均值、標(biāo)準(zhǔn)差、變異系數(shù)分布類型等的經(jīng)典統(tǒng)計分析(見表1)。

表1 切溝密度、平均坡度的經(jīng)典統(tǒng)計結(jié)果

結(jié)果可見，切溝密度與平均坡度的分布服從對數(shù)正態(tài)分布，變異系數(shù)分別為90%與68%，有一定的差異，但均屬中等程度的變異[23]。因此，首先對切溝密度與平均坡度的原始值做對數(shù)變換，使其符合標(biāo)準(zhǔn)正態(tài)分布以便于進(jìn)一步的數(shù)據(jù)分析。

當(dāng)主輔變量間的統(tǒng)計相關(guān)性超過0.45的中等相關(guān)程度時，協(xié)同克立格法的插值結(jié)果精度即明顯優(yōu)于普通克立格法。對切溝密度與平均坡度做相關(guān)性分析(見圖1)發(fā)現(xiàn)，二者的Person相關(guān)性，即r值為0.54，達(dá)到了0.01(雙側(cè))的極顯著水平，因此適合選擇平均坡度作為輔變量，應(yīng)用協(xié)同克立格法對切溝密度進(jìn)行插值。

圖1 切溝密度與平均坡度的回歸分析

2.2空間變異特征分析

圖2是切溝密度和平均坡度的地統(tǒng)計學(xué)半方差函數(shù)擬合模型及其參數(shù)，及協(xié)同區(qū)域化變量交互半方差函數(shù)的擬合模型及其參數(shù)。切溝密度、平均坡度與交叉半方差函數(shù)的理論模型均符合球狀模型。協(xié)同區(qū)域化變量可以是正相關(guān)，也可以是負(fù)相關(guān)，這與協(xié)同區(qū)域化現(xiàn)象發(fā)生的具體過程有關(guān)，本研究的結(jié)果都是正相關(guān)的。決定系數(shù)r2反映所選模型對半方差函數(shù)的擬合程度，協(xié)同區(qū)域化變量的交互半方差函數(shù)對試驗變異函數(shù)的擬合程度從0.76提高到了0.90，說明交互半方差函數(shù)一方面很好的反映了切溝密度的空間結(jié)構(gòu)特性，另一方面可以提高空間插值結(jié)果的精度。

圖2 半變異函數(shù)圖及其模型擬合結(jié)果

參數(shù)理論模型塊金值基臺值塊金值/基臺值/%變程/m決定系數(shù)殘差切溝密度球狀模型0.861.7250409370.140.139平均坡度球狀模型0.722.0136520000.760.149協(xié)同區(qū)域化變量球狀模型0.020.3633223000.900.024

C0為塊金值，表示由采樣誤差、短距離的變異、隨機(jī)和固有變異引起的基底效應(yīng)；C0+C為基臺值，即Sill；塊金值與基臺值的比值[C0/(C0+C)]表示空間自相關(guān)性程度[24]。切溝密度、平均坡度的C0/(C0+C)值分別為0.50與0.36，表現(xiàn)為中等強(qiáng)度的空間自相關(guān)性，說明獲取的切溝密度、平均坡度分布是由結(jié)構(gòu)性因素(如氣候、母質(zhì)、地形、土壤類型等)和隨機(jī)性因素(如耕作措施、土地利用等各種人為活動)共同作用的結(jié)果[25]。協(xié)同區(qū)域化變量的C0/(C0+C)值為0.056，小于0.10，表明二者具有強(qiáng)烈的空間自相關(guān)性，結(jié)構(gòu)性因素遠(yuǎn)大于隨機(jī)性因素引起的空間異質(zhì)性。另一方面，協(xié)同區(qū)域化變量的C0/(C0+C)比切溝密度的值小，表明輔變量的應(yīng)用增強(qiáng)了結(jié)構(gòu)性因素造成的互相關(guān)作用。

RSS(Residue Sums of Squares)為殘差平方和，表示試驗?zāi)Ｐ团c理論模型的擬合程度，簡稱殘差。殘差最小是回歸模型所追求的目標(biāo)。從殘差來看，盡管輔變量平均坡度的殘差(0.149)大于主變量切溝密度的殘差(0.139)，但協(xié)同區(qū)域化變量的殘差顯著(0.024)小于主變量殘差(0.139))，說明交互半方差函數(shù)的擬合效果要明顯優(yōu)于半方差函數(shù)，這也表明協(xié)同區(qū)域化變量的空間結(jié)構(gòu)性要優(yōu)于單一變量。

本研究中輔變量平均坡度的C0/(C0+C)值、RSS值均差于主變量切溝密度的值，但協(xié)同區(qū)域化變量的r2值、C0/(C0+C)值、RSS值均優(yōu)于主變量，說明了即使在主輔變量的空間結(jié)構(gòu)性中等情況下，如果主輔變量間互相關(guān)性達(dá)到中等相關(guān)程度(相關(guān)性0.76)，協(xié)同區(qū)域化變量的的空間結(jié)構(gòu)性就可優(yōu)于單一變量。

2.3 切溝密度的空間分布特征分析

圖3是距離權(quán)重反比法、普通克立格法、協(xié)同克立格法對切溝密度進(jìn)行空間插值形成的空間分布圖?？梢钥闯觯?種方法得到的切溝密度分布特征是一致的，均為條帶狀和斑塊狀格局?？傮w來說，海倫市切溝密度基本走向是北低南高，西低東高，并在東南部與南部有兩塊高值區(qū)域。西部切溝密度值最低，低于0.6 km/km2，東部與東南部值最高，超過了2.4 km/km2。導(dǎo)致這種空間分布相似性的原因是多方面的。一方面，坡度是決定侵蝕溝形成的位置和發(fā)育的大小的主要因素，而海倫市東北為小興安嶺西南麓，坡度較大；西面與松嫩平原接壤，坡度較緩。另一方面，耕作是加速切溝侵蝕的重要人為因素[26]。海倫市耕地的開墾順序是從南向北，南部地區(qū)的開墾時間最長，開墾強(qiáng)度大。

2.4 不同插值方法的評價

盡管三種插值方法反映的切溝密度總體分布特征相近，但協(xié)同克立格法得到的結(jié)果更為精細(xì)，準(zhǔn)確。比較圖3a，3b，3c可知，距離權(quán)重反比法與普通克立格法插值形成的切溝密度空間分布具有較強(qiáng)的相似性，均無密度大于2.7 km/km2的高值區(qū)域，且斑塊化現(xiàn)象嚴(yán)重。而在協(xié)同克立格法得到的結(jié)果中，受輔變量的協(xié)助，對局部變異細(xì)節(jié)的描述更為詳細(xì)，層次感鮮明；得到的信息更為豐富，局部多中心聚集現(xiàn)象更顯著。表3列出的是距離權(quán)重反比法、普通克立格法與協(xié)同克立格法對切溝密度的預(yù)測精度的比較。從檢驗結(jié)果來看，距離權(quán)重反比法的預(yù)測值與實測值之間的相關(guān)系數(shù)比普通克立格法的稍高，從0.01提高到0.04，但總體而言是很低的。這可能是切溝密度的取樣密度不高導(dǎo)致的。二者的均方根誤差相近，均為0.94左右。協(xié)同克立格法預(yù)測所產(chǎn)生的均方根誤差明顯低于距離權(quán)重反比法和普通克立格法，從0.94左右降低到了0.74，減少了20%以上。但預(yù)測值與實測值之間的相關(guān)系數(shù)有了明顯的提高，增加了89%以上。總之協(xié)同克立格法相對于距離權(quán)重反比法與普通克立格法的預(yù)測精度有較大的提高，預(yù)測的結(jié)果明顯。

圖3 不同插值方法下的切溝密度分布

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

3 結(jié)論與討論

本文比較了距離權(quán)重反比法、普通克里格法與協(xié)同克里格法對切溝密度做空間插值的預(yù)測精度，結(jié)論如下：

(1) 切溝密度、平均坡度均呈中等變異強(qiáng)度；二者表現(xiàn)為中等程度的統(tǒng)計相關(guān)性，達(dá)到了極顯著的水平，因此適合選取平均坡度作為輔變量，應(yīng)用協(xié)同克立格法對切溝密度進(jìn)行插值。

(2) 切溝密度、平均坡度與協(xié)同區(qū)域化變量的半方差函數(shù)模型均符合球狀模型；切溝密度與協(xié)同區(qū)域化變量受結(jié)構(gòu)性因素的影響遠(yuǎn)大于隨機(jī)性因素，均為強(qiáng)空間自相關(guān)性；而平均坡度受結(jié)構(gòu)性因素與隨機(jī)性因素的影響相近，表現(xiàn)為中等空間自相關(guān)性；協(xié)同區(qū)域化變量表現(xiàn)為正相關(guān)，且交互半方差函數(shù)的塊金值、基臺值都有所降低，協(xié)同區(qū)域化變量的空間結(jié)構(gòu)性要優(yōu)于單一變量，對試驗變異函數(shù)的擬合程度有較大的提高。

(3) 距離權(quán)重反比法與普通克里格法的預(yù)測精度相近，這可能與采樣樣方的密度不高有關(guān)；協(xié)同克立格法對于切溝密度的局部變異細(xì)節(jié)的描述更為詳細(xì)，更接近切溝密度的真實分布情況；與二者相比，協(xié)同克立格法的均方根誤差均降低20%強(qiáng)，相關(guān)系數(shù)提高89%以上，表明利用平均坡度采用協(xié)同克立格法可以提高切溝密度的預(yù)測精度。

本文結(jié)合切溝密度與平均坡度的關(guān)系，應(yīng)用不同空間插值方法對切溝密度做插值，得到了以流域平均坡度為輔助變量，應(yīng)用協(xié)同克里格對切溝密度做空間插值的精度大于距離權(quán)重反比法與普通克里格法的結(jié)果。將來的工作應(yīng)當(dāng)放在影響切溝形成與發(fā)展其他的環(huán)境因素，如距水系的距離，以及人為影響因素，如耕作措施等的影響。

[1] 唐克麗,史立人,史德明,等.中國水土保持[M].北京:科學(xué)出版社,2004.

[2] 張興義,劉曉冰,隋躍宇,等.人為剝離黑土層對大豆干物質(zhì)積累及產(chǎn)量的影響[J].大豆科學(xué),2006,25(2):123-126.

[3] 水利部,中國科學(xué)院,中國工程院.中國水土流失防治與生態(tài)安全:東北黑土區(qū)卷[M].北京:科學(xué)出版社,2010.

[4] 孟令欽,李勇.東北黑土區(qū)溝蝕研究與防治[J].中國水土保持,2009(12):40-42.

[5] 朱顯謨.黃土區(qū)土壤侵蝕的分類[J].土壤學(xué)報,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] 劉秉正,吳發(fā)啟.土壤侵蝕[M].陜西:陜西人民出版社,1997.

[9] 中華人民共和國水利部.SL446—2009黑土區(qū)水土流失綜合防治技術(shù)標(biāo)準(zhǔn)[S].北京:中國水利水電出版社,2009.

[10] 趙英時,陳冬梅,李小文,等.遙感應(yīng)用分析原理與方法[M].北京:科學(xué)出版社,2003.

[11] 李浩,張興義,劉爽,等.典型黑土區(qū)村級尺度侵蝕溝演變[J].中國水土保持科學(xué),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] 李 新,程國棟,盧玲.空間內(nèi)插方法比較[J].地球科學(xué)進(jìn)展,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] 王政權(quán).地質(zhì)統(tǒng)計學(xué)及在生態(tài)學(xué)中的應(yīng)用[M].北京:科學(xué)出版社,1999.

[18] 謝葉偉,劉兆剛,趙軍,等.基于RS與GIS的典型黑土區(qū)土地利用變化分析：以海倫市為例[J].地理科學(xué),2010(3):428-434.

[19] 林忠輝,莫興國,李宏軒.中國陸地區(qū)域氣象要素空間插值[J].地理學(xué)報,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

國家自然科學(xué)基金(41171230，41071201)

王平(1965—),女,黑龍江省哈爾濱人,學(xué)士,高級工程師,主要從事土壤侵蝕機(jī)理及水土流失監(jiān)測技術(shù)方面研究。E-mail:wp626588@163.com

張興義(1966—),男,黑龍江省密山人,博士,研究員,主要從事黑土侵蝕和保護(hù)研究。E-mail:zhangxy@iga.ac.cn

P941

：A

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