GIS與地統(tǒng)計(jì)學(xué)的土壤水分空間插值方法

張 優(yōu), 王 娟, 張 杰, 彭文甫*, 朱亞蘭

(1. 四川師范大學(xué) 地理與資源科學(xué)學(xué)院, 四川 成都 610066；2. 四川師范大學(xué) 西南土地資源評價(jià)與監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室, 四川 成都 610066)

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

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

1 研究區(qū)概況

綿竹市位于四川盆地西北部,位于山地平原過渡地區(qū),地處30°09′—31°42′N,103°54′—104°20′E之間,面積1 245.3 km2,轄19鎮(zhèn)2鄉(xiāng),總?cè)丝?2萬.本文研究區(qū)為綿竹市部分區(qū)域(圖1).該區(qū)域?qū)偎拇ㄅ璧刂衼啛釒駶櫄夂騾^(qū),氣候溫和,年均溫15.7℃,年降水量為1 053.2 mm,雨量充沛,四季分明,大陸季風(fēng)性氣候顯著.地貌形態(tài)區(qū)域分異十分明顯,西北部為山地,東南部為平原,地勢西北高,東南低,由西北向東南傾斜.森林占幅員面積49.8%.土地利用類型主要以林地、耕地、草地、建設(shè)用地、水體和未利用地6類為主.境內(nèi)平原區(qū)土壤包括13個(gè)土屬,8個(gè)亞類,土地肥沃,物產(chǎn)豐富.

圖 1 研究區(qū)位置圖

2 數(shù)據(jù)準(zhǔn)備及預(yù)處理

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

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

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

表 1 土壤含水量數(shù)據(jù)統(tǒng)計(jì)結(jié)果

3 研究方法

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

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

(1)

Wi=d-p

(2)

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

Kriging方法中權(quán)重Wi的一般形式為

W

Z(xi+h)]2,

(3)

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

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

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

(4)

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

因此多元回歸預(yù)測結(jié)果為

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

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

(5)

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

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

圖 2 土壤含水量QQ正態(tài)分布

圖 3 土壤含水量半變異函數(shù)擬合

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

4 結(jié)果及精度對比分析

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

由表4可知,對于RMSE,預(yù)測誤差表現(xiàn)為：

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

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

5 討論

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

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

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