GIS與地統計學的土壤水分空間插值方法

張 優, 王 娟, 張 杰, 彭文甫*, 朱亞蘭

(1. 四川師范大學 地理與資源科學學院, 四川 成都 610066；2. 四川師范大學 西南土地資源評價與監測教育部重點實驗室, 四川 成都 610066)

土壤是由多種因素相互作用而形成的,其屬性在空間分布上有一定的差異[1].實現精確農業養分管理和解決全球變化等環境問題的關鍵在于準確掌握土壤理化性狀的空間變異規律[2].土壤干濕狀況對監測土壤退化、土地干旱有著巨大作用,是監測土地狀況的重要指標[3].土壤水分的時空分布研究對農業干旱監測和生態環境監測具有重大意義.目前,揭示土壤屬性空間分布特征的主要手段仍是通過野外采樣和室內測定.但其問題在于,無論采樣密度多大,均無法得到所有區域的土壤屬性值.而空間插值技術可以將離散的數據點轉化為連續的數據曲面從而實現研究區土壤屬性的全覆蓋.作為土壤屬性變化的時空定量監測方法,土壤屬性空間插值方法及其精度是數字土壤研究的重點.姚雪玲等[4]以黃土丘陵溝壑區洋圈溝小流域為例,介紹了基于ArcGIS與地統計學的回歸克里格模型插值方法,詳細介紹了基于虛擬變量技術將分類變量引入模型的過程,同時將回歸克里格模型、回歸線性模型和普通克里格模型的插值精度進行對比分析.史文嬌等[5]介紹了土壤屬性空間插值的常用方法與精度驗證的方法和指標,并且總結了能夠提高土壤屬性插值精度的6種途徑.陳思萱等[6]針對廣東省某地區土壤中采樣樣品點的重金屬污染物砷(As),選取統計學中具有代表性的3種空間插值方法,對比分析不同空間插值方法下土壤重金屬濃度估算精度和污染格局空間區劃差異.龍軍等[1]根據不同地貌類型特征在福建省各地級市選取了9個典型縣,利用2008年采集的樣點數據評價不同空間插值方法對土壤有機質含量插值精度的影響.

本文以山地平原過渡地區四川省綿竹市為研究區,選擇了幾種常用的空間插值方法,包括反距離權重法(IDW)、徑向基函數空間插值方法、普通克里格插值方法以及結合線性回歸模型和普通克里格的回歸克里格插值方法,探究不同插值方法對山地平原過渡地帶土壤水分的影響.

1 研究區概況

綿竹市位于四川盆地西北部,位于山地平原過渡地區,地處30°09′—31°42′N,103°54′—104°20′E之間,面積1 245.3 km2,轄19鎮2鄉,總人口52萬.本文研究區為綿竹市部分區域(圖1).該區域屬四川盆地中亞熱帶濕潤氣候區,氣候溫和,年均溫15.7℃,年降水量為1 053.2 mm,雨量充沛,四季分明,大陸季風性氣候顯著.地貌形態區域分異十分明顯,西北部為山地,東南部為平原,地勢西北高,東南低,由西北向東南傾斜.森林占幅員面積49.8%.土地利用類型主要以林地、耕地、草地、建設用地、水體和未利用地6類為主.境內平原區土壤包括13個土屬,8個亞類,土地肥沃,物產豐富.

圖 1 研究區位置圖

2 數據準備及預處理

研究涉及數據主要包括：綿竹市行政區劃矢量圖,來源于四川省測繪地理信息局；綿竹市2016年土地利用類型數據為Landsat 8 OLI解譯數據,解譯為耕地、林地、草地、水域、建設用地和未利用地共6類[7],經實地樣點驗證,一級土地利用解譯精度滿足后續分析需要(總體精度達94.3%)；綿竹市數字地面模型(DEM)來自USGS(https://www.usgs.gov/),空間分辨率為30 m×30 m.

綿竹市土壤含水量樣點數據通過野外采集獲取,考慮到樣點空間分布對插值結果的影響,在樣點選取時按以下思路進行.首先在研究區行政區劃圖上創建2 km×2 km的網格,以各個格網中心作為預采樣點,然后疊合土地利用數據,對預采樣點進行適當調節以規避諸如水體、城市不透水面等區域.此外,為適應研究區西北部復雜地形條件下土壤水分空間分布的潛在異質性,采樣時對樣點進行了適當加密,最終確定在空間上盡可能均勻分布的樣點共184個.結合手持GPS,土壤采樣工作在2016年4～5月間進行,采樣時選取有代表性的新鮮土壤,采樣深度為0～20 cm.取樣時,先鏟出一個耕層斷面,再平行于斷面下鏟取土,混合土樣太多時用四分法將多余土壤棄去[8],采樣的同時記錄樣點的土地利用類型.采樣當天即采用烘干法進行土樣含水量測定,具體過程：將鮮土捏碎放入已知準確質量的鋁盒中,對其進行稱量,精確至0.001 g；然后,將其置于已預熱至(105±2) ℃的烘箱中烘烤12 h,取出冷卻至室溫后立即稱重,每個樣點做3個平行測定.若第一次測定結果中3個平行樣本之間的差距過大(超過0.05%),則進行二次烘干再次稱重[9].

利用ArcGIS對土壤含水量樣點數據進行數據子集提取,為滿足不同插值方法對數據的要求,本文按照7∶3進行插值模型輸入數據和插值精度驗證數據子集提取.插值模型輸入數據集(記為Ⅰ)共129個樣點,精度驗證數據集(記為Ⅱ)共55個樣點.表1為土壤含水量數據的主要統計量信息.從表1可看出,山地平原過渡地帶的土壤水分質量分數的變異系數為19%,屬于中等變異.

表 1 土壤含水量數據統計結果

3 研究方法

3.1 空間插值方法選擇空間插值方法的選擇是影響插值精度的主要因素之一[10].不同空間插值方法在數據要求、參數設置等方面均有較大差異[11].本文選擇了幾種常用的空間插值方法,包括：反距離權重法(IDW)、徑向基函數函數空間插值方法、普通克里格插值方法以及結合線性回歸模型和普通克里格的回歸克里格插值方法.

3.1.1反距離權重(IDW)和克里格插值法(Kriging) 反距離權重法是一種常用又簡單的空間確定性插值方法,待插值點的預測值為搜索半徑中已知點的值以與待插值點的距離為權重的加權平均,距離待插值點越近的已知樣點的權重越大[12].與此類似,克里格空間插值(Kriging)方法亦通過一定范圍內的實測樣點數據的線性加權預測待插值點,所不同的是,各點的權重由待插值變量的半方差函數得到.IDW和Kriging空間插值方法的一般模型可表示為

(1)

Wi=d-p

(2)

其中,dSi-S0為搜索半徑中的點Si到待插值點S0的距離,p為參數.因此,當通過搜索半徑或空間插值樣點數目確定參數N時,參數p成為IDW插值方法的主要影響因素.

Kriging方法中權重Wi的一般形式為

W

Z(xi+h)]2,

(3)

其中,N(h)為空間距離為h時點對數目,Z(xi)和Z(xi+h)分別為位置xi和(xi+h)處的屬性采樣值.實際上,Wi(h)可以看作是屬性空間自相關的度量函數,稱為半變異函數.

3.1.2徑向基函數插值方法(RBF) RBF徑向基函數插值方法是一種確定性空間插值方法,其基本思路是用一個通過所有屬性樣點的曲面來預測待插值點,這個擬合曲面擁有最小的空間曲率.徑向基函數包括5種常用的基本函數,包括規則樣條函數、張力樣條函數、高次曲面函數、反高次曲面樣條函數和平面樣條函數.本文以交叉驗證精度為標準選取最優的徑向基函數插值結果作為對比.

3.1.3回歸克里格插值方法 地物屬性受多自然因素甚至社會人文因素的影響,且這些影響因影響因素的空間異質、時間變率不定而對不同區域具有不同的影響深度.一些影響因素空間分布廣、影響時間尺度大、空間影響較為均一,使得受其影響的地物屬性具有明顯的空間趨勢.而一些區域性或短時間尺度的因子則多體現為局部影響,這就在由大尺度因素構建的“趨勢”上疊加了短程變異.回歸克里格插值方法將這種多尺度因素影響下的地物屬性進行分離,其概念模型為

(4)

研究表明,太陽輻射、地形、土地利用類型、坡度、坡向、海拔高程等是影響土壤水分的重要因子[13].結合數據的可獲取性和研究區特點,本文將土地利用類型、太陽輻射和坡度因子作為土壤含水量“趨勢”擬合的考慮因子.需要注意的是,土地利用類型數據是一種“類別集合”形式的空間“軟數據”,利用其進行土壤含水量多因素線性回歸時需要對其進行“硬化”[14].本文引入虛擬變量對土地利用類型軟數據進行“硬化“.在研究區中,結合采樣時記錄的土地利用現狀,最終確定6種土地利用類型：耕地、河灘、未利用地、草地、林地、濕地.根據顯著性差異,分為4類:林地歸為第一類(C1),耕地歸為第二類(C2),荒地歸為第三類(C3),河灘地和草地歸為第四類(C4).對各個土地利用類型進行土壤含水量方差分析,以此為參考將土地利用類型做聚合歸并,并進行虛擬變量處理(表2)和回歸因子顯著性檢驗(表3).

因此多元回歸預測結果為

μ(fi)=17.780+[landusetype].weight+

0.006*[solar]+0.05*[slope].

(5)

多元回歸預測結果和樣點實測結果的差值即為局部、短時間尺度因素引起的變異部分(ε(gj)),對其進行空間插值即得到待差值點的變異值,本文采用普通克里格插值方法作為土壤含水量變異值的空間分布.最后根據(5)式得到最終的空間插值結果.

3.2 數據探索和插值參數確定為合理設定空間插值模型參數、滿足空間插值模型對輸入數據的要求(比如克里格插值要求樣點數據滿足正態分布),對樣點數據的統計特征、空間分布特征及異常值進行探索是空間插值的必要過程.正態QQ圖是用于反映樣點數據與標準正態分布的接近程度,樣點數據的分布越接近一條直線,則越接近于服從正態分布.同時,圖中遠離標準直線的點為潛在的異常值點,因而QQ正態圖還可用于異常值剔除.圖2為本文土壤含水量樣點數據經異常值剔除后的正態QQ圖,可以看出,樣點數據具有很好的正太分布特點.如前所述,IDW和Kriging插值方法對待插值點的預測是基于空間距離的,對IDW插值模型中距離參數(或點數量)N和Kriging方法中的權重參數Wi(h)的確定均要求對樣點的空間分布特征進行探索.圖3為本文土壤含水量樣點數據的半變異函數擬合結果.從圖3可以看出,研究區土壤含水量最大空間自相關距離(變程)約為1 600 m,因此以1 600 m為IDW空間插值方法中的搜索半徑閾值.

圖 2 土壤含水量QQ正態分布

圖 3 土壤含水量半變異函數擬合

3.3 精度評價指標空間插值精度評價指標可分為相對驗證指標和絕對驗證指標.本文選擇絕對驗證指標有均方根誤差和平均絕對百分比誤差,均方根誤差(RMSE)表示預測值與實測值偏差的平方和觀測次數比值的平方根,均方根誤差對預測中的最大值和最小值比較敏感；平均絕對百分比誤差(MAPE)表示預測值和實測值的偏差占實測值的百分比的平均根.相對驗證指標選取模擬優度指數(G)它是衡量模型模擬效果是否優于僅對實測值取平均值得的指標.指標選取的優劣對于空間插值精度評價有重要的影響.

4 結果及精度對比分析

為了明確不同插值方法對山地平原過渡帶土壤水分空間數據的影響,通過幾種插值方法得到結果如圖4,并將不同的插值精度進行比較,篩選出適宜過渡帶的土壤水分空間插值方法.

由表4可知,對于RMSE,預測誤差表現為：

表 4 4種插值方法的精度評估

RBF>RK>IDW>Kriging,說明Kriging的預測誤差最小,RBF的預測誤差最大；對于MAPE,預測誤差表現為：RBF>RK>IDW>Kriging,說明Kriging的精度最高,RBF的預測精度最低；對于G,預測誤差表現為：RK>Kriging>IDW>RBF,說明回歸克里格插值的精度(0.36)高于普通克里格法(0.27)和反距離權重插值法(0.19)與樣條插樣法(0.08),其預測模型效果最優.不同檢驗指標精度評價不同,絕對驗證指標顯示,克里格插值精度最高,其次為反距離權重法,回歸克里格插值次之,徑向基函數效果最差.而相對驗證指標顯示,回歸克里格插值高于普通克里格與其他2種方法,說明驗證指標的選取也會導致評價精度的不同.

5 討論

山地平原過渡帶具有明顯地域分異,特殊的地形地貌造成土壤水分含量區域范圍內的差異.不同的插值方法插值精度不同,在眾多插值方法中,克里格插值仍然是一種比較適用的方法.在目標變量空間自相關性較弱,與環境因子相關性較強的時候,回歸克里格是較為優越的插值方法.本文的回歸克里格考慮了山地平原過渡帶的土壤水分受多種因素的影響,嘗試將坡度、年太陽輻射與土地利用類型等環境因子引入模型.反距離權重法是一種空間確定性插值方法,插值表面中的最小值和最大值出現在采樣點處,輸出表面對異常值十分敏感.目前對此模型使用的合理性證明仍較困難,但反距離權重法假定數據中存在空間自相關,當目標變量出現局部變異性時,此方法較為適宜.徑向基函數法可為平緩變化的表面生成很好的結果,但在表面值在短距離內出現劇烈變化或懷疑樣本值很可能有測量誤差或不確定性時,這種方法不適用.從其插值效果來看,山地平原過渡帶的土壤水分空間分布特征為西北山區的土壤水分含量明顯高于東南平原,城鎮周邊土壤水分含量明顯低于周圍耕地.

土壤水分對農作物的生長發育有著重要的影響,研究區的土壤水分質量分數為13.22%～72.52%,對于中生性作物而言,對土壤含水量的要求各不相同.馬鈴薯等豆類作物等最適土壤含水量相當于田間持水量的70%～80%,禾谷類作物為60%～70%.土壤含水量低于最適值時,光合作用降低.各種作物光合作用開始降低時的土壤含水量(占田間持水量之百分數)分別為：水稻57%,大豆45%,小麥41%,花生32%.因此,對于綿竹地區而言,馬鈴薯等豆類植物適宜種植在西北山區,這里土壤水分充足,利于作物的生長；對于東南平原而言,適宜種植水稻、大豆、小麥、花生等糧食作物和經濟作物.土壤水分含量的探討對于農事活動的安排具有重要的作用.

土壤屬性的空間連續數據在理論與實踐方面有著重要的作用,如何提高其精度是未來發展的一個重要研究方向.影響空間插值精度的原因除了插值方法的選擇外,采樣點的密度和數目[15]、插值方法中參數的設定以及空間自相關的范圍和程度都是影響其精度的因素[16].對于克里格插值方法而言,獲得穩定的半方差函數需要50～100個采樣點,采樣點位100～150個才可達到一個平穩的變異函數[17].插值參數的設定對其精度有一定的影響,但是目前絕大多數參數的確定還未有太多的規律可循.