復雜地形區土壤有機質空間變異性分析及制圖

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

(1.青海師范大學 地理科學學院,西寧 810000;2.青海省自然地理與環境過程重點實驗室,西寧 810008;3.青藏高原地表過程與生態保育教育部重點實驗室,西寧 810008;4.青海省人民政府—北京師范大學 高原科學與可持續發展研究院,西寧 810008)

有機質在土壤中扮演著重要角色，是植物營養中的重要來源，它是促進植物生長發育、土壤保肥性和緩沖性、改善土壤物理性質的主要因素，目前土壤有機質在改善農業、管理耕地等方面具有重要意義，是土壤學、地球化學、農學研究的重點內容之一[1-2]。

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

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

湟水流域位于青藏高原與黃土高原過渡帶，其地形起伏大、地貌類型多樣，主要以農業為主。目前對該研究區進行土壤屬性空間制圖的研究還比較少見。代子俊[14]分別利用1985年和2015年的土壤采樣數據，對湟水流域土壤全氮進行克里格插值制圖，但由于其土壤采樣點較少且在流域西北區域基本沒有采樣點，這在一定程度上將影響空間制圖的精度；同時雖然其在土壤采樣時考慮了土地利用類型因素，但在空間插值時并未將其結合到土壤全氮制圖中。基于以上的考慮，本文對湟水流域土壤有機質含量進行空間地統計學制圖，旨在研究所使用方法在復雜地形區的適用性以及比較協同克里格方法、地理加權回歸克里格法與結合土地利用類型后的協同克里格法在土壤有機質制圖方面的精度，為湟水流域提供較為準確的土壤有機質數據支持。

1 數據與方法

1.1 研究區概況

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

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

1.2 數據來源與預處理

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

圖1 研究區數據

1.3 研究方法

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

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

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

(1)

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

(2)

式中：Z(xi)為有機質實測值；Q(xi)為拆分趨勢項；R(xi)為拆分殘差項；Z(i)為未采樣點的土壤有機質預測值。

1.3.3 地理加權回歸克里格法 地理加權回歸克里格是在地理加權回歸的基礎上進行預測，所以本文首先對地理加權回歸模型進行擬合：

(3)

式中：YGWRK(ui，vi)為因變量在點(ui，vi)處的估計值；xij為第i個解釋變量在(ui，vi)處的值；β0(ui，vi)為樣本點在(ui，vi)的特征彈性系數；βp(ui，vi)為第i個樣點上的第p個回歸系數。擬合的關鍵參數主要是核函數與帶寬，經過測試，本文使用的核函數為自適應高斯核函數，帶寬以AICc赤池信息量準則為判斷標準，值越小說明擬合后信息的不確定性越小[13]。地理加權回歸的模型擬合參數見表1。

表1 地理加權回歸模型擬合參數

1.3.4 半變異函數 半變異函數又稱半方差函數，是關于數據點間的變異性距離的函數，也是描述區域化變量結構性、隨機性的基本手段。估算方程式如下：

(4)

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

1.3.5 空間插值 本文創建訓練集后首先對數據進行正態分布檢驗，均值與中值殘差、GWRK殘差、原值均符合正態分布，可以用于空間插值。

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

1.3.6 模型精度評價 本文通過內部驗證與外部驗證進行精度評價，內部驗證：標準均方根(Root-Mean-Square Standardized，RMSS)預測誤差和1最接近、標準平均值誤差(Mean Standardized Error，MSE)與0最接近、均方根(Root-Mean-Square，RMS)預測誤差達到最小、平均標準誤差(Average Standard Error，ASE)和均方根(Root-Mean-Square，RMS)預測誤差最相近。外部驗證：平均絕對值誤差(Mean Absolute Error，MAE)越小模型精度越高、均方根誤差(Root Mean Squares Error，RMSE)越小模型越穩定[20]，準確度(Accuracy，AC)可評價預測的準確性，越接近于1說明預測越準確，取值范圍為0～1[7]。r為預測集與驗證集的相關系數，RI為比較GWRK相對于COK的提高度。

(5)

(6)

(7)

(8)

(9)

(10)

2 結果與分析

2.1 土壤有機質描述性統計分析

從土壤有機質描述性統計分析(表2)可以看出采樣點的土壤有機質含量在7～150 g/kg，標準差為24.2，平均值為34.16 g/kg。

表2 土壤樣本有機質描述性統計

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

表3 樣本所對應的土地利用類型的描述性統計

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

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

2.2 空間變異特征分析

首先利用土地利用類型，GWR模型得到殘差項，然后利用ArcGIS 10.2地統計模塊建立訓練子集209個，驗證子集31個，最后利用ArcGIS 10.2地統計模塊的協同克里格插值，插值過程中結合湟水流域DEM因子，擬合最優的半變異函數模型(表5)。本文4種方法均為指數模型。

表5 SOM空間變異的理論模型與參數

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

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

2.3 模型精度分析

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

表6 模型預測精度評價

2.4 有機質空間插值及結果與分析

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

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

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

圖2 土壤有機質含量分布

3 討論與結論

本研究表明復雜地形區的湟水流域2016年SOM平均含量為34.16 g/kg，屬于中等水平。變異系數為90.90%，存在中等空間變異性，說明作為重要的農業區域，其土壤有機質同時受土壤母質、土壤類型以及土地利用管理方式、施肥等人類活動的共同影響。

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

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