復(fù)雜地形區(qū)土壤有機(jī)質(zhì)空間變異性分析及制圖

張 歡,高小紅,2,3,4

(1.青海師范大學(xué) 地理科學(xué)學(xué)院,西寧 810000;2.青海省自然地理與環(huán)境過(guò)程重點(diǎn)實(shí)驗(yàn)室,西寧 810008;3.青藏高原地表過(guò)程與生態(tài)保育教育部重點(diǎn)實(shí)驗(yàn)室,西寧 810008;4.青海省人民政府—北京師范大學(xué) 高原科學(xué)與可持續(xù)發(fā)展研究院,西寧 810008)

有機(jī)質(zhì)在土壤中扮演著重要角色，是植物營(yíng)養(yǎng)中的重要來(lái)源，它是促進(jìn)植物生長(zhǎng)發(fā)育、土壤保肥性和緩沖性、改善土壤物理性質(zhì)的主要因素，目前土壤有機(jī)質(zhì)在改善農(nóng)業(yè)、管理耕地等方面具有重要意義，是土壤學(xué)、地球化學(xué)、農(nóng)學(xué)研究的重點(diǎn)內(nèi)容之一[1-2]。

傳統(tǒng)的土壤預(yù)測(cè)制圖是一個(gè)手工制圖的過(guò)程，需要對(duì)地形圖等進(jìn)行目視解譯以及對(duì)土壤類型斑塊進(jìn)行手工勾繪，除此之外還需要野外核查才能編輯成數(shù)據(jù)圖件[3]，需要耗費(fèi)大量的人力物力財(cái)力。數(shù)字土壤制圖是以土壤與景觀理論為基礎(chǔ)，借助數(shù)學(xué)、空間分析方法等技術(shù)手段的土壤制圖與調(diào)查方法[4]，制圖方法主要包括基于地統(tǒng)計(jì)的數(shù)字土壤制圖和基于專家知識(shí)的數(shù)字土壤制圖[3]。研究表明：地統(tǒng)計(jì)學(xué)方法在數(shù)字土壤制圖中不僅考慮空間上的鄰近點(diǎn)相關(guān)性，也考慮了該點(diǎn)的其他地理要素相關(guān)性，可以較好的用于數(shù)字土壤制圖[4]。許多學(xué)者陸續(xù)地使用GIS與地統(tǒng)計(jì)方法探討土壤屬性的空間特征[2]。趙永存等[5]發(fā)現(xiàn)回歸克里格的預(yù)測(cè)效果最好，泛克里格次之，而多元線性回歸方法最差，證明地統(tǒng)計(jì)學(xué)插值方法用于空間預(yù)測(cè)相對(duì)較好；席江勇等[6]在對(duì)土壤有機(jī)質(zhì)的插值研究中利用了確定性內(nèi)插和地統(tǒng)計(jì)插值進(jìn)行分析，使用確定性插值產(chǎn)生了牛眼，而地統(tǒng)計(jì)插值消除了這種現(xiàn)象，結(jié)果連續(xù)且光滑，能更好的反映其在空間上的變化；綜上所述，在研究土壤養(yǎng)分時(shí)考慮地統(tǒng)計(jì)插值方法能夠取得較好的預(yù)測(cè)精度。

目前的一些研究表明，將土地利用類型結(jié)合到地統(tǒng)計(jì)學(xué)插值中能夠獲得較好的預(yù)測(cè)精度。文雯等[7]研究發(fā)現(xiàn)在黃土丘陵地區(qū)，利用土地利用修正克里格插值方法精度比其他插值方法精度高，能夠提高土壤屬性空間制圖的精度；顧成軍等[8]利用有限的采樣點(diǎn)將土地利用結(jié)合到克里格插值中，發(fā)現(xiàn)在該研究區(qū)土地利用對(duì)有機(jī)碳的空間分異有很好的解釋能力且不需要獲取其他因子；趙建華[9]在研究蘭州市土壤全氮和有機(jī)質(zhì)的空間變異時(shí)發(fā)現(xiàn)其空間分布特征與研究區(qū)的各種土地利用類型有重要的聯(lián)系。復(fù)雜地形區(qū)不同于地形比較簡(jiǎn)單的區(qū)域可以使用地統(tǒng)計(jì)分析的簡(jiǎn)單克里格進(jìn)行插值預(yù)測(cè)，需要考慮有關(guān)的輔助因子來(lái)進(jìn)行土壤屬性空間預(yù)測(cè)研究，以提高預(yù)測(cè)精度[10]。協(xié)同克里格法基于目標(biāo)變量與輔助變量，是一種考慮輔助變量的有效方法[11]，而在實(shí)際情況中，空間數(shù)據(jù)是隨著地理位置的變化而變化的,地理加權(quán)回歸模型正是考慮了空間數(shù)據(jù)的位置信息[12]，近年來(lái)在土壤屬性空間制圖中得到較好的應(yīng)用。地理加權(quán)回歸克里格是地理加權(quán)回歸與克里格法的有機(jī)結(jié)合，該方法不僅考慮了空間位置局部信息，還考慮了代表隨機(jī)性的殘差數(shù)據(jù)，能夠反映出更加真實(shí)的土壤屬性空間變異[13]。

湟水流域位于青藏高原與黃土高原過(guò)渡帶，其地形起伏大、地貌類型多樣，主要以農(nóng)業(yè)為主。目前對(duì)該研究區(qū)進(jìn)行土壤屬性空間制圖的研究還比較少見(jiàn)。代子俊[14]分別利用1985年和2015年的土壤采樣數(shù)據(jù)，對(duì)湟水流域土壤全氮進(jìn)行克里格插值制圖，但由于其土壤采樣點(diǎn)較少且在流域西北區(qū)域基本沒(méi)有采樣點(diǎn)，這在一定程度上將影響空間制圖的精度；同時(shí)雖然其在土壤采樣時(shí)考慮了土地利用類型因素，但在空間插值時(shí)并未將其結(jié)合到土壤全氮制圖中?；谝陨系目紤]，本文對(duì)湟水流域土壤有機(jī)質(zhì)含量進(jìn)行空間地統(tǒng)計(jì)學(xué)制圖，旨在研究所使用方法在復(fù)雜地形區(qū)的適用性以及比較協(xié)同克里格方法、地理加權(quán)回歸克里格法與結(jié)合土地利用類型后的協(xié)同克里格法在土壤有機(jī)質(zhì)制圖方面的精度，為湟水流域提供較為準(zhǔn)確的土壤有機(jī)質(zhì)數(shù)據(jù)支持。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

湟水流域位于青藏高原與黃土高原的交接地帶，位于36°02′—37°28′N，100°42′—103°04′E，地形在南北方向，具有“三山兩谷”獨(dú)特的構(gòu)造景觀，三山從南向北為拉脊山、達(dá)坂山、祁連山，兩谷為介于達(dá)坂山與拉脊山之間的湟水干流河盆谷地、達(dá)坂山與祁連山之間的大通河河谷。青海省境內(nèi)湟水流域面積約為16 120 km2，海拔高度在1 664～4 882 m，地勢(shì)呈現(xiàn)東南低西北高的趨勢(shì),具有多樣復(fù)雜的地形,黃土地貌在該區(qū)最為典型，區(qū)域內(nèi)土地利用類型多樣，在淺山區(qū)主要以農(nóng)田為主,川水區(qū)以農(nóng)田和城鄉(xiāng)工礦居住用地為主，海拔3 000 m以上主要是草地占多,湟水流域土壤類型主要是灰鈣土、黑鈣土和栗鈣土為主。

湟水流域年平均氣溫在2.5～7.8℃，該區(qū)冬季氣溫低、夏季氣溫適宜，最高氣溫可達(dá)30℃左右，而最低氣溫可達(dá)-28℃左右，日照長(zhǎng)、可吸收輻射大，且從東到西逐漸增大；區(qū)域平均降水量486mm，主要集中在6—9月份，10月到次年2月降水量最少，只占全年降水量的10%，而蒸發(fā)量最大可達(dá)1 000mm，且隨海拔高度的下降蒸發(fā)量是逐漸增加的，表現(xiàn)出西北向東南蒸發(fā)量增加的趨勢(shì)[15]。該區(qū)是青海省的主要農(nóng)業(yè)耕作區(qū)，也是青海省人口最密集、經(jīng)濟(jì)最發(fā)達(dá)的地區(qū)。

1.2 數(shù)據(jù)來(lái)源與預(yù)處理

本文使用的數(shù)據(jù)主要有數(shù)字高程模型(Digital elevation model，DEM)(圖1A)、土壤采樣數(shù)據(jù)253個(gè)(圖1B)、土地利用數(shù)據(jù)(圖1C)、土壤類型數(shù)據(jù)(圖1D)、湟水流域邊界數(shù)據(jù)。其中土壤類型矢量數(shù)據(jù)通過(guò)1∶100萬(wàn)青海省土壤類型圖矢量數(shù)字化獲??；土壤有機(jī)質(zhì)來(lái)源于2016年10—11月份野外采集數(shù)據(jù)，采樣土層厚度為0—20 cm，土壤有機(jī)質(zhì)(Soil organic matter，SOM)的測(cè)試分析請(qǐng)參見(jiàn)李冠穩(wěn)[16]。將采樣數(shù)據(jù)的野外經(jīng)緯度坐標(biāo)(GPS獲取，WGS-84)轉(zhuǎn)換為UTM投影坐標(biāo)，DEM來(lái)源于USGS網(wǎng)站(http:∥www.usgs.gov/)，分辨率為30 m，對(duì)DEM進(jìn)行填洼、拼接、裁剪等預(yù)處理，最終得到湟水流域DEM數(shù)據(jù)，利用ArcGIS軟件計(jì)算出坡度(圖1E)、坡向(圖1F)以及高程；研究區(qū)土地利用類型數(shù)據(jù)根據(jù)2016年Landsat8/OLI影像通過(guò)隨機(jī)森林分類方法獲取，總體精度腦山區(qū)、淺山區(qū)、川水區(qū)分別為91.33%，92.09%，87.85%，見(jiàn)馬慧娟論文[17]，對(duì)該土地利用類型數(shù)據(jù)進(jìn)行合并后包括7類，即水澆地、旱地、草地、林地、水域、城鄉(xiāng)居住建設(shè)用地、未利用土地，并利用土壤有機(jī)質(zhì)采樣點(diǎn)提取相應(yīng)的土地利用類型數(shù)據(jù)。

圖1 研究區(qū)數(shù)據(jù)

1.3 研究方法

1.3.1 描述性統(tǒng)計(jì)與單方差分析 首先對(duì)土壤有機(jī)質(zhì)數(shù)據(jù)以及土地利用數(shù)據(jù)利用SPSS 17.0進(jìn)行描述性統(tǒng)計(jì)分析；統(tǒng)計(jì)出最大值、最小值以及峰度、偏度等數(shù)據(jù)值，并通過(guò)K-S檢驗(yàn)數(shù)據(jù)是否呈正態(tài)分布，由于有機(jī)質(zhì)數(shù)據(jù)不呈正態(tài)分布，對(duì)其進(jìn)行對(duì)數(shù)轉(zhuǎn)換或者剔除異常值，直到數(shù)據(jù)呈現(xiàn)正態(tài)分布。使用單因素方差分析比較主要土地利用類型的土壤有機(jī)質(zhì)是否具有顯著差異[18]，如果p>0.05，則土地利用類型來(lái)自于同一個(gè)正態(tài)總體，不用消除由地類差異引起的土壤有機(jī)質(zhì)差異，確保土地利用類型的有機(jī)質(zhì)可以看做是同一個(gè)區(qū)域化變量在空間上的分布[19]。

1.3.2 結(jié)合土地利用類型的協(xié)同克里格法 土地利用數(shù)據(jù)屬于分類變量或定性變量，不具備數(shù)學(xué)上的大小關(guān)系，因此不能直接用于地統(tǒng)計(jì)學(xué)分析。根據(jù)前人研究[19]，有兩種方法可將分類變量轉(zhuǎn)換為定量變量，第一種方法為空間分區(qū)，該方法只適合大面積連續(xù)的區(qū)域；第2種方法為對(duì)不同土地利用類型的土壤屬性值進(jìn)行均值、中值中心化，并利用其殘差進(jìn)行插值，本文研究區(qū)屬于復(fù)雜地形區(qū)，地形起伏大、土地類型斑塊破碎，不像平原地區(qū)具有大面積連續(xù)的單一類型區(qū)域，所以采用第2種方法將定性變量轉(zhuǎn)為定量變量，公式如下：

R(xi)=Z(xi)-Q(xi)

(1)

Z(i)=R(i)+Q(i)

(2)

式中：Z(xi)為有機(jī)質(zhì)實(shí)測(cè)值；Q(xi)為拆分趨勢(shì)項(xiàng)；R(xi)為拆分殘差項(xiàng)；Z(i)為未采樣點(diǎn)的土壤有機(jī)質(zhì)預(yù)測(cè)值。

1.3.3 地理加權(quán)回歸克里格法 地理加權(quán)回歸克里格是在地理加權(quán)回歸的基礎(chǔ)上進(jìn)行預(yù)測(cè)，所以本文首先對(duì)地理加權(quán)回歸模型進(jìn)行擬合：

(3)

式中：YGWRK(ui，vi)為因變量在點(diǎn)(ui，vi)處的估計(jì)值；xij為第i個(gè)解釋變量在(ui，vi)處的值；β0(ui，vi)為樣本點(diǎn)在(ui，vi)的特征彈性系數(shù)；βp(ui，vi)為第i個(gè)樣點(diǎn)上的第p個(gè)回歸系數(shù)。擬合的關(guān)鍵參數(shù)主要是核函數(shù)與帶寬，經(jīng)過(guò)測(cè)試，本文使用的核函數(shù)為自適應(yīng)高斯核函數(shù)，帶寬以AICc赤池信息量準(zhǔn)則為判斷標(biāo)準(zhǔn)，值越小說(shuō)明擬合后信息的不確定性越小[13]。地理加權(quán)回歸的模型擬合參數(shù)見(jiàn)表1。

表1 地理加權(quán)回歸模型擬合參數(shù)

1.3.4 半變異函數(shù) 半變異函數(shù)又稱半方差函數(shù)，是關(guān)于數(shù)據(jù)點(diǎn)間的變異性距離的函數(shù)，也是描述區(qū)域化變量結(jié)構(gòu)性、隨機(jī)性的基本手段。估算方程式如下：

(4)

式中：h是分隔距離；N(h)是(xi+h,xi)之間計(jì)算樣本變異函數(shù)值的樣本對(duì)數(shù)(如Z(xi+h)和Z(xi)為計(jì)算中的第一對(duì)數(shù)據(jù))[12]。

1.3.5 空間插值 本文創(chuàng)建訓(xùn)練集后首先對(duì)數(shù)據(jù)進(jìn)行正態(tài)分布檢驗(yàn)，均值與中值殘差、GWRK殘差、原值均符合正態(tài)分布，可以用于空間插值。

分別采用均值、中值消除趨勢(shì)項(xiàng)，并用得到的插值圖加上不同土地利用信息的趨勢(shì)項(xiàng)，最終得到結(jié)合土地利用類型的土壤有機(jī)質(zhì)預(yù)測(cè)圖。同時(shí)，進(jìn)行地理加權(quán)回歸模型擬合，得到的殘差項(xiàng)進(jìn)行空間插值，最后與協(xié)同克里格的插值結(jié)果比較，制作土壤有機(jī)質(zhì)分布圖；本文稱4種方法為均值修正協(xié)同克里格(Land Use Mean Modification Cokriging，MMCOK _LU)、中值修正協(xié)同克里格(Land Use Median Correction Cokriging，MCCOK_LU)、協(xié)同克里格(Cokriging，COK)、地理加權(quán)回歸克里格(Geographically Weighted Regression Kriging，GWRK)。

1.3.6 模型精度評(píng)價(jià) 本文通過(guò)內(nèi)部驗(yàn)證與外部驗(yàn)證進(jìn)行精度評(píng)價(jià)，內(nèi)部驗(yàn)證：標(biāo)準(zhǔn)均方根(Root-Mean-Square Standardized，RMSS)預(yù)測(cè)誤差和1最接近、標(biāo)準(zhǔn)平均值誤差(Mean Standardized Error，MSE)與0最接近、均方根(Root-Mean-Square，RMS)預(yù)測(cè)誤差達(dá)到最小、平均標(biāo)準(zhǔn)誤差(Average Standard Error，ASE)和均方根(Root-Mean-Square，RMS)預(yù)測(cè)誤差最相近。外部驗(yàn)證：平均絕對(duì)值誤差(Mean Absolute Error，MAE)越小模型精度越高、均方根誤差(Root Mean Squares Error，RMSE)越小模型越穩(wěn)定[20]，準(zhǔn)確度(Accuracy，AC)可評(píng)價(jià)預(yù)測(cè)的準(zhǔn)確性，越接近于1說(shuō)明預(yù)測(cè)越準(zhǔn)確，取值范圍為0～1[7]。r為預(yù)測(cè)集與驗(yàn)證集的相關(guān)系數(shù)，RI為比較GWRK相對(duì)于COK的提高度。

(5)

(6)

(7)

(8)

(9)

(10)

2 結(jié)果與分析

2.1 土壤有機(jī)質(zhì)描述性統(tǒng)計(jì)分析

從土壤有機(jī)質(zhì)描述性統(tǒng)計(jì)分析(表2)可以看出采樣點(diǎn)的土壤有機(jī)質(zhì)含量在7～150 g/kg，標(biāo)準(zhǔn)差為24.2，平均值為34.16 g/kg。

表2 土壤樣本有機(jī)質(zhì)描述性統(tǒng)計(jì)

CV為變異系數(shù)，根據(jù)CV<10%為弱的變異性,10%～100%屬于中等變異,CV>100%是強(qiáng)變異性，可知總體樣本CV為90.90%，屬于中等變異強(qiáng)度；偏度為2.047，屬于右偏態(tài)，峰度4.761，不是特別陡峭；在K-S檢驗(yàn)中，p的值等于0，說(shuō)明其不呈正態(tài)分布，不能直接用于空間插值研究，所以需要進(jìn)行對(duì)數(shù)轉(zhuǎn)換，轉(zhuǎn)換后其p值為0.055，大于0.05，可以用于地統(tǒng)計(jì)分析；除此之外根據(jù)表3對(duì)各土地利用類型的統(tǒng)計(jì)，土壤有機(jī)質(zhì)都處于中等變異水平，林地的變異程度大于其他3種土地利用類型。

表3 樣本所對(duì)應(yīng)的土地利用類型的描述性統(tǒng)計(jì)

不同土地利用類型的土壤有機(jī)質(zhì)單因素方差分析(表4)表明，由于p>0.05,四組土地利用的SOM并未存在很大的差異，故SOM值是屬于同一總體的，使水澆地、林地、旱地、草地的SOM值可以視為同一區(qū)域化變量，不用在插值前消除由這些類所引起的SOM差異，故其殘差值也無(wú)明顯的不同，都可以看成在空間上分布的同一個(gè)區(qū)域化變量[19]。

表4 不同土地利用類型SOM單因素方差分析

2.2 空間變異特征分析

首先利用土地利用類型，GWR模型得到殘差項(xiàng)，然后利用ArcGIS 10.2地統(tǒng)計(jì)模塊建立訓(xùn)練子集209個(gè)，驗(yàn)證子集31個(gè)，最后利用ArcGIS 10.2地統(tǒng)計(jì)模塊的協(xié)同克里格插值，插值過(guò)程中結(jié)合湟水流域DEM因子，擬合最優(yōu)的半變異函數(shù)模型(表5)。本文4種方法均為指數(shù)模型。

表5 SOM空間變異的理論模型與參數(shù)

結(jié)構(gòu)、隨機(jī)性變異是空間變異的兩大模塊，而C0，C0+C，C0/(C0+C)常常用來(lái)表示半變異函數(shù)建模的空間變異程度，其中C0為塊金值，表示可能的隨機(jī)性，這種隨機(jī)因素可能為土地的耕作、施肥、管理措施等；C結(jié)構(gòu)方差由自然因素，即地形、地質(zhì)、土壤母質(zhì)、土壤類型等引起的土壤性狀的結(jié)構(gòu)變異；C0+C為基臺(tái)值，空間總的變異一般由基臺(tái)值表示；C0/(C0+C)為塊金與基臺(tái)值之比，稱為塊金系數(shù)，其大小表示空間的變異程度，值小則表示空間變異程度大多是結(jié)構(gòu)性因素引起的，反之值大則為隨機(jī)性因素引起[2]。

研究區(qū)MMCOK_LU，MCCOK_LU，COK模型的塊金系數(shù)分別為17.8%，19.9%，19.3%，屬于強(qiáng)的空間相關(guān)性，GWRK的塊金系數(shù)為47.5%，屬于中等的空間相關(guān)性；前3個(gè)模型表明土壤有機(jī)質(zhì)空間變異主要由結(jié)構(gòu)性因素引起，而最后一個(gè)模型表明土壤有機(jī)質(zhì)空間變異主要由結(jié)構(gòu)性和隨機(jī)性因素引起，原因可能是在考慮土地利用類型后，空間細(xì)節(jié)信息更豐富，人為的影響更突出[C0/(C0+C)>75%:弱空間相關(guān)；C0/(C0+C)<25%:強(qiáng)空間相關(guān)；之間屬于中等空間相關(guān)]。變程表示變量在空間上的自相關(guān)性，MMCOK_LU和MCCOK_LU的變程均為6 899 m，COK為6 790 m，GWRK為5 977 m，表明在這個(gè)范圍內(nèi)空間自相關(guān)是存在的，均選擇指數(shù)模型作為最佳模型。

2.3 模型精度分析

由表6可知RMSS標(biāo)準(zhǔn)均方根預(yù)測(cè)誤差MMCOK_LU最接近于1，4種方法的MSE相差不大，COK，GWRK僅比MMCOK_LU與MCCOK_LU小0.001，COK的ASE與RMS接近程度為0.003 9，大于MMCOK_LU與MCCOK_LU；MAE平均絕對(duì)誤差MCCOK _LU最小，GWRK最大，說(shuō)明MCCOK _LU的模型精度最高，其次分別為MMCOK_ LU，COK；RMSE均方根誤差GWRK最小，MCCOK_LU最大，說(shuō)明GWRK比COK，MMCOK_LU，MCCOK_LU穩(wěn)定性高；而對(duì)于模型的準(zhǔn)確度，由表6中的AC可知預(yù)測(cè)準(zhǔn)確度由大到小依次為MCCOK_LU>GWRK>MMCOK_LU>COK；與COK相比，GWRK相對(duì)提高3.3%。綜上所述，可以發(fā)現(xiàn)結(jié)合土地利用中值殘差的MCCOK_LU模型精度最高。

表6 模型預(yù)測(cè)精度評(píng)價(jià)

2.4 有機(jī)質(zhì)空間插值及結(jié)果與分析

圖2中可以看出4種方法的SOM預(yù)測(cè)值分布范圍相差不大，COK(協(xié)同克里格)與GWRK(地理加權(quán)回歸克里格)的土壤有機(jī)質(zhì)空間分布較為連續(xù)，而MMCOK_LU(均值修正協(xié)同克里格)、MCCOK_LU(中值修正協(xié)同克里格)分布比較破碎，但是反應(yīng)的地類細(xì)節(jié)特征以及信息的豐富程度較好；MCCOK_LU高值區(qū)相較于MMCOK_LU少，GWRK的低值區(qū)最少，四者的分布中間值最多。

4種方法所產(chǎn)生的高值均分別集中在大通縣西北部、海晏縣南部和湟源縣西北部等地區(qū)，互助縣中西部以及湟中縣北部也有零星高值出現(xiàn)；其中大通縣西北部、大通縣與海晏縣交界處的高值可以結(jié)合DEM圖、土地利用圖以及土壤類型圖(圖1)進(jìn)行分析，該區(qū)海拔相對(duì)較高，薛曉娟等[21]的研究發(fā)現(xiàn)，海拔3 400 m以上土壤有機(jī)質(zhì)含量開(kāi)始急劇上升；且該區(qū)受人為因素影響小，多為林地和草地(動(dòng)植物殘?bào)w腐爛分解為腐殖質(zhì)),所以該地區(qū)主要為SOM豐富的地區(qū)，且該區(qū)主要發(fā)育高山草甸土，通過(guò)青海土壤[22]高山草甸土有機(jī)質(zhì)含量在不同的區(qū)域差異較大但是普遍含量較高，這是因?yàn)楦呱讲莸橥恋耐寥佬纬蛇^(guò)程慢，表層粗的有機(jī)質(zhì)易于積累，所以有機(jī)質(zhì)含量較高，該區(qū)土壤有機(jī)質(zhì)含量大約在38.8 g/kg左右，預(yù)測(cè)結(jié)果與該值相符；在海晏縣南部地區(qū)和湟源縣西北部地區(qū)、互助縣西部、西納川、石崖莊主要分布栗鈣土，部分山地草甸土與高山草甸土，該區(qū)栗鈣土的土壤有機(jī)質(zhì)含量較高，分解微生物的能力較強(qiáng)[14]，根據(jù)青海土壤[22]該區(qū)水熱條件較好，土壤較為肥沃，土壤有機(jī)質(zhì)在40～55 g/kg左右，預(yù)測(cè)結(jié)果在50 g/kg左右，兩者大致相符，且土地利用類型主要是草地、林地、部分旱地和水澆地，草地、林地有機(jī)質(zhì)含量豐富，而旱地作物生長(zhǎng)周期比較短，施肥較多，所以SOM值會(huì)偏高；民和縣西部、西南部、樂(lè)都縣東南部主要土地利用類型為耕地與草地，栗鈣土分布廣，根據(jù)代子俊[14]的研究，發(fā)現(xiàn)該區(qū)農(nóng)作物產(chǎn)出率有所提高，離不開(kāi)當(dāng)?shù)剞r(nóng)民對(duì)農(nóng)田的精細(xì)管理和施肥，所以增加了土壤有機(jī)質(zhì)的含量，該區(qū)SOM含量大致在55 g/kg左右。有機(jī)質(zhì)低值區(qū)主要出現(xiàn)在大通回族土族自治縣西南部、民和縣東南部、樂(lè)都縣南部地區(qū)；民和縣東南部，海拔最低，坡度相對(duì)較大且處于流域下游，如果遇有雨水，則對(duì)土壤的沖刷能力較大，而該區(qū)的栗鈣土熱量有余，水分不足，水熱矛盾突出，SOM含量低，不易積累，加之水土流失強(qiáng)[22]，所以會(huì)呈現(xiàn)出較低值；大通回族土族自治縣西南部主要分布黑鈣土與部分栗鈣土，土地利用類型主要是坡旱地，海拔在2 500～3 000 m左右，土溫較高，土壤的水分較適中，土地耕作后透氣性加快，使SOM礦化加快，除此之外該地耕種時(shí)間長(zhǎng)，耕層經(jīng)常燒灰，土壤顏色變淺，表層的SOM含量有明顯下降趨勢(shì)[22]，所以該地有機(jī)質(zhì)含量不是很高；其他地區(qū)有機(jī)質(zhì)含量基本處于預(yù)測(cè)值中間水平，即27 g/kg左右。

水澆地區(qū)域有機(jī)質(zhì)含量受水量的影響較大，豐水期水儲(chǔ)量大時(shí)會(huì)加速土壤的分解，會(huì)產(chǎn)生更多的有機(jī)質(zhì)，此時(shí)有機(jī)質(zhì)比較豐富，但是研究區(qū)處于干旱半干旱區(qū)，降水變化大，水量不是很充足，且日照大蒸發(fā)快，所以水澆地的SOM含量隨氣候的不同而發(fā)生變化；結(jié)合圖2與圖1C發(fā)現(xiàn)草地的有機(jī)質(zhì)含量在大部分地區(qū)高于林地，通過(guò)文獻(xiàn)[23]發(fā)現(xiàn)主要原因是該地溫性草原分布廣，由放牧等人為因素影響較大，動(dòng)物糞便等都會(huì)對(duì)SOM產(chǎn)生影響。

圖2 土壤有機(jī)質(zhì)含量分布

3 討論與結(jié)論

本研究表明復(fù)雜地形區(qū)的湟水流域2016年SOM平均含量為34.16 g/kg，屬于中等水平。變異系數(shù)為90.90%，存在中等空間變異性，說(shuō)明作為重要的農(nóng)業(yè)區(qū)域，其土壤有機(jī)質(zhì)同時(shí)受土壤母質(zhì)、土壤類型以及土地利用管理方式、施肥等人類活動(dòng)的共同影響。

對(duì)于地形復(fù)雜區(qū)，利用土地類型信息輔助研究土壤屬性空間分布制圖的文獻(xiàn)還較少，本文除了考慮高程外，還將土地利用信息對(duì)SOM的影響結(jié)合到有機(jī)質(zhì)空間變異性分析與空間制圖中，研究表明，利用較易獲取的土地利用類型研究SOM的空間變異性，能有效的降低空間異常值的影響，提升模型的預(yù)測(cè)精度，其插值精度均較高。所使用的4種方法中結(jié)合土地利用類型的中值修正協(xié)同克里格精度最高，其次為地理加權(quán)回歸克里格、均值修正協(xié)同克里格、協(xié)同克里格,模型的準(zhǔn)確度從高到低分別為0.923，0.909，0.905，0.883，GWRK相對(duì)COK提高3.3%。本研究中的土地利用中值修正協(xié)同克里格制圖精度要好于均值修正，這與吳子豪等[19]研究土壤有機(jī)碳密度的空間異質(zhì)性結(jié)果基本一致；文雯等[7]在研究黃土丘陵小流域時(shí)，中值修正克里格相對(duì)于未修正克里格準(zhǔn)確度提升0.042 4，本文中的中值修正克里格相對(duì)于未修正克里格準(zhǔn)確度提升0.040，稍低于前者的研究，原因可能是研究區(qū)相對(duì)于黃土丘陵區(qū)地形更加復(fù)雜，且研究中只考慮了4種土地利用類型，差異相對(duì)較?。槐疚睦脜f(xié)同克里格與輔助變量相結(jié)合，提高了模型的預(yù)測(cè)精度，與前人的研究結(jié)果一致[24]。

地理加權(quán)回歸克里格、協(xié)同克里格均能有效地對(duì)復(fù)雜地形區(qū)的湟水流域土壤有機(jī)質(zhì)進(jìn)行制圖，但前者精度高于后者。協(xié)同克里格與土地利用類型結(jié)合后，中值協(xié)同克里格與均值協(xié)同克里格預(yù)測(cè)精度均稍高于普通協(xié)同克里格。因此，考慮土地利用類型信息的地統(tǒng)計(jì)學(xué)在復(fù)雜地形區(qū)土壤有機(jī)質(zhì)制圖是非常有效的。研究區(qū)的土地利用類型為7類，對(duì)SOM的貢獻(xiàn)程度分別為：草地為0.121，水澆地為0.063，旱地為0.035，林地為參照(將土地利用類型從名義變量轉(zhuǎn)換為定量變量，做出土地利用類型的回歸方程，得出土地利用方式的貢獻(xiàn)程度)，貢獻(xiàn)程度最高的草地，采樣點(diǎn)數(shù)僅占總樣本數(shù)的17%，小于一些學(xué)者[7]的研究，采樣點(diǎn)空間分布不均勻以及個(gè)別樣本較少是本文的不足。除此之外本文未將定性變量考慮進(jìn)地理加權(quán)回歸模型，目前的研究未能很好解決的這個(gè)問(wèn)題，今后將對(duì)其進(jìn)一步的研究。