基于空間插值的區域尺度土壤全鹽量時空變異規律研究

徐存東，谷豐佑，朱興林，王 燕,3，張 鵬，高德坤

(1.華北水利水電大學水利學院，鄭州 450046；2.水資源高效利用與保障工程河南省協同創新中心，鄭州 450046；3.河南省水工結構安全工程技術研究中心，鄭州 450046)

0 引 言

我國西北地區土地廣闊、光熱充足，制約大片宜耕土地利用的主要因素是可用水資源短缺，為了改善耕地環境，提高土地利用效率，發展區域農業經濟，必須大力發展引水灌溉工程[1]。同時，長期以來人工提灌工程的大面積實施，使得區域水文循環規律隨著灌溉方式和引水規模變化而不斷改變，也主導著區域尺度的土壤水鹽分異過程，改變并控制著區域的水鹽時空演化和分布特征[2]。由于該地區高蒸發低降雨的氣候條件，促使土壤深層的水分攜鹽分向土壤表層遷移，隨著水分的不斷蒸發，鹽分滯留在土壤表層，最終導致大面積的土地鹽漬化；為了灌區的可持續發展，有必要對區域內土壤鹽分的空間分布和演化態勢進行科學的評估和預測。當前，土壤全鹽量野外監測點數量較為稀少，無法形成對整個區域土壤鹽分的大規模監測，從中獲取的有限數據不能滿足對其時空變化的研究，且監測點在空間上分布不均勻，通過空間插值模型可以有效地對土壤鹽分空間分布特征進行預測[3]。如姚榮江、楊勁松等學者利用普通克里金法對蘇北海涂圍墾區土壤鹽分空間分布進行了模擬、插值和預測，為濱海地區鹽漬土壤的精確評估與高效改良提供了參考依據[4]；劉鑫等運用反距離權重法對河套灌區的土壤進行空間插值，實現了鹽化分布的空間可視化[5]；吳亞坤等運用反距離權重法對南疆綠洲區土壤鹽分進行了分布解析得到了較好的插值效果[6]；張飛以新疆精河綠洲為研究區，運用反距離權重法對表層土壤鹽分進行空間插值，得到了鹽漬化程度與分布情況[7]；姜涵等運用GIS 空間分析法對研究區土壤鹽分特征進行研究發現可以通過插值得到未檢測地區的數據，能夠有效的提高對土壤全鹽量的研究[8]。上述研究僅是使用了單一插值方法對研究區土壤全鹽量進行分析，沒有涉及不同插值方法下所得模擬結果的精度比較，不能驗證插值結果的準確性；同時由于研究區的土壤鹽分監測點分布較為稀疏，無法全面監測研究區，構建出的土壤鹽分變化空間插值模型，在進行土壤鹽分時空變化特征研究時有一定的局限性。本文以甘肅省景電一期灌區(以下稱“景電灌區”)為研究區，利用反距離權重法、趨勢面法、樣條函數法、普通克里金法四種不同的空間插值方法對研究區2015年土壤全鹽量數據進行插值，將插值結果與監測點的實際觀測數據進行對比分析，優選出區域土壤全鹽量最佳空間插值方法。并分別對1994年、2001年、2008年和2015年監測數據進行空間插值，對研究區土壤全鹽量的空間分布及變化規律進行分析，以期為灌區的土壤鹽漬化研究提供數據和理論支持。

1 材料與方法

1.1 研究區概況

景電灌區是地處西北干旱荒漠區的大型揚水灌區，東經103°20′～104°04′，北緯37°26′～38°41′，位于甘、寧、蒙的三省的交界地帶。北倚騰格里沙漠，東臨黃河, 南靠昌嶺山, 西連石羊河流域。地層區劃屬河西走廊六盤山分區, 武威中寧小區。總灌溉面積6.13 萬hm2，分布高程在1 596～1 906 m之間。屬于典型的內陸寒旱區氣候，少雨干旱、蒸發量大、植被稀疏、風積風蝕嚴重是其主要氣候特征。研究區地表水資源十分匱乏，多年平均年降雨量僅為186 mm，且集中在5-9 月份，蒸發量卻高達2 516 mm。研究區地表土層結構松散，為土壤鹽分溶質運移提供了良好的條件。加之低降雨高蒸發的氣候特征，鹽分因子極易留存于土層表面形成積鹽，日積月累造成土壤鹽漬化[9]。

圖1 研究區土壤全鹽量監測點位置圖

1.2 數據來源

根據《景泰川電力提灌工程第一期工程技施設計報告》《景泰川灌區歷年土地調查報告》和《景泰縣荒地資源及其開發利用報告》，結合相關資料可以得到1994、2001、2008和2015年的0～100 cm土層的土壤全鹽量和50個非插值點土壤全鹽量數據，表1為各監測點點位坐標及各個點位地下100 cm內土壤全鹽量數據。

1.3 研究原理與方法

1.3.1 插值方法

現有的監測點數量相較于整個研究區面積略顯稀少，且分布不勻，有限的數據無法準確反映研究區域土壤全鹽量空間分布狀況及時空變化規律。因此根據現有的研究區監測點數據，利用ArcGIS中的空間插值模塊得出整個研究區域的土壤全鹽量插值圖，進而分析區域尺度土壤鹽分發展特征[10]。ArcGIS空間分析工具中有四種常用的空間插值方法：反距離權重法(IDW)、趨勢面法(TL)、樣條插值法(Splines)和克里金(Kriging)插值法[11]。反距離權重法(IDW)通過對已知點與預測點之間的距離進行賦權得到插值結果，賦權大小與距離遠近呈負相關[12]；趨勢面法(TL)以多項式回歸分析為基礎，遵循最小二乘法原理，對輸入的監測點數據進行擬合，反映出土壤鹽分的空間發展趨勢[13]；樣條插值法分為規則樣條插值(RS)和張力樣條插值(TS)，是由相鄰的兩個數據點通過數學函數確定一個多項式，多項式兩兩相連組成一條經過所有數據點的光滑樣條曲線[14]；克里金法是在區域化變量和半變異函數基礎上，對區域內的變量取值進行線性無偏最小二乘估計，以此反映區域化變量在空間的真實變化特性，其中又以普通克里金法(OK)的應用最為廣泛[15]。

表1 研究區監測點點位及其土壤全鹽量

1.3.2 精度評價方法

空間插值技術的發展為土壤全鹽量空間分析提供了新的方法，然而如何從多種插值方法中選取該研究區域的最優插值模型是空間插值在土壤全鹽量分析中的一個關鍵問題[16]。本文運用交叉驗證法對幾種插值方法進行精度評價。首先從16個監測點中選出一個點，設其數據為未知值，使用其他15個監測點實測數據構建插值模型，讀取該點的插值結果，計算插值結果與所選監測點實際數據之間的誤差值。依次從16個監測點中選取下一個點作為未知站點，重復進行操作，獲取這16個監測點誤差值。即將現有監測點數據分為已知值和驗證值，通過不同插值方法對已知值進行插值，讀取驗證值點位對應的插值結果，與驗證值進行對比，插值結果與實際數據之間的誤差越小則插值方法越精確[17]。

采用平均絕對誤差MAE，平均相對誤差MRE，均方根誤差RMSE來評價插值結果的精度[18]。

(1)

(2)

(3)

式中：Pei、Pai分別為第i個監測點的預測值和實測值；n為監測點數量。

2 結果與分析

2.1 不同插值方法比較與優選

為了更為準確地顯示研究區100 cm內土壤全鹽量在水平空間的分布情況，得出土壤全鹽量最優插值方法。運用ArcGIS中的Spatial Analyst模塊對2015年土壤全鹽量數據進行四種空間插值分析，繪制出各插值方法100 cm內土壤鹽分的空間分布圖。

樣條函數插值結果見圖2、圖3，從圖中可以看出運用樣條函數法進行插值時，規則樣條函數插值中權重為1的時候正負向外延較大，而權重為0時正負向外延較小，出現此現象主要是由于在插值時缺乏數據基礎，監測點數據分布不均勻且各點之間變化不夠平緩；而張力樣條曲線中權重越大越接近真實區間。但是兩種函數類型的插值結果負向數值始終為負數，與土壤全鹽量必須大于等于零的實際不符。

趨勢面法插值結果見圖4，趨勢面法通過數學函數(各階多項式)進行插值，可以看出在運用趨勢面插值法時，采用二階多項式插值時與真實數據擬合最好，一階多項式插值結果區間甚至小于實測數據真實區間。但是由于趨勢面法過于注重整體發展趨勢，導致插值結果負向數值為負數，不符合實際。

普通克里金插值結果見圖5，可以看出普通克里金插值結果綜合考慮了距離、方位的影響，對極值點的處理較好，所有插值結果都沒有出現負值，且區間與實驗數據真實區間基本吻合，僅有OK(指數函數)的區間略小于真實區間。

圖2 張力樣條函數插值圖

圖3 規則樣條函數插值圖

反距離權重法見圖6，綜合所有插值結果可得運用反距離權重法進行插值時，結果直觀，與真實區間符合。

運用誤差矩陣(Error Matrix)對3種反距離權重法和其他4種普通克里金法插值結果進行精度評價，以土壤全鹽量0.5%為一階段進行劃分，將誤差矩陣劃分為7個階段區。對50個非插值點的數據進行檢驗。用基于誤差矩陣的精度評價方法對土壤鹽分插值結果與非插值點數據進行評價，kappa系數高于0.60為高度一致性，低于0.60為中度一致性。由表2可知，插值總精度反距離權重法(冪3)最高(為76.00%)，其次是普通克里金法三角函數(為66.00%)，反距離權重法(冪1)最低(為56.00%)。以Kappa系數進行判定，僅有反距離權重法(冪3)屬于高度一致性，其余插值結果均屬于中等一致性，故反距離權重法(冪3)對土壤鹽分插值精度最高。

為選取土壤全鹽量最優插值方法，本研究采用交叉驗證法通過平均相對誤差(MAE)、平均絕對誤差(MRE)、均方根誤差(RMSE)3個指標來評價不同插值方法插值結果精度。從表3中可以看出，從MRE來看，TS>RS>TL>OK>IDW，說明IDW的插值相對誤差最小，張力樣條函數法最大；從MAE來看，RS>TS>TL>OK>IDW，即規則樣條函數法插值誤差最大，IDW的絕對誤差最小；從RMSE來看，TS>RS>TL>IDW>OK，說明OK的靈敏度最高，張力樣條函數法最低，但IDW靈敏度僅低于OK，且在數值上相差不大。綜合3個誤差評價指標以及誤差矩陣分析結果，可得IDW(冪3)相比較于其他方法在插值精度和預測準確度上均為最佳，所以選擇IDW(冪3)對土壤全鹽量進行空間插值。

圖4 趨勢面法插值圖

表2 土壤鹽分插值精度評價表

2.2 區域土壤鹽分插值結果驗證

運用ArcGIS軟件使用反距離權重法(冪3)對2015年土壤鹽分進行空間插值，提取插值圖中相應監測點的預測值，選取最近的年份2015年，將所得預測值與實測值進行對比，結果如圖7所示。可以看出，運用反距離權重法(冪3)對研究區土壤全鹽量進行空間插值得到的預測值與所選年份的各個測點的實際值基本一致。說明利用反距離權重法(冪3)得到的結果圖可以準確的反映研究區各個年份的土壤全鹽量空間分布特征，并可以此分析土壤全鹽量的年際動態變化規律。

圖5 普通克里金插值圖

圖6 反距離權重插值圖

表3 5種插值方法插值土壤鹽分結果的誤差統計

注：MRE為平均相對誤差；MAE為平均絕對誤差；RMSE為均方根誤差。

圖7 2015年各監測點預測值與實際值對比圖

2.3 區域土壤全鹽量時空進程分析

運用反距離權重法(冪3)對1994、2001、2008和2015年土壤全鹽量進行空間插值，插值結果見圖8。可以看出，土壤全鹽量呈現西南低東北高分布狀況，從1994年到2015年南部低鹽分含量地區不斷向北發展至邊界，西部土壤全鹽量在0.5%～1.0%的地區逐年減少為一較小區域，東北部高鹽分含量地區面積雖無明顯增長，但其土壤全鹽量一直在增大，從1994年的2.0%～2.5%增大到2015年的2.5%～3.0%，并且局部區域的土壤全鹽量已增大到了3.0%～3.5%。

圖8 各年份土壤全鹽量空間插值結果圖(單位：%)

通過ArcGIS對圖8中各年份的土壤全鹽量特征值進行提取。由表4可知，研究區土壤全鹽量最小值從1994年到2015年前期增長后期下降，截止2015年相比于1994年僅提高了0.000 8%，幾乎沒有發生變化。最大值由1994至2015年持續增長，到2015年較1994年提高了0.964 6%。土壤全鹽量平均值總體上發生了增長，2015年較1994年總體增長了0.116 7%。

表4 研究區各年份土壤鹽分特征值 %

3 結 論

(1)文中所述的5種插值方法中，規則樣條函數、張力樣條函數和趨勢面法插值結果負向數值為負數，與實際情況不符；普通克里金法除了以指數函數為插值函數的插值結果外，都與實驗數據基本吻合，但是反距離權重法結果更加直觀，且與真實區間符合，經交叉驗證可知，反距離權重法插值結果與實際的監測結果吻合度較高，結果最為精確，可較準確預測灌區土壤全鹽量的變化趨勢，為灌區區域尺度鹽分空間處理分析提供指導。

(2)使用反距離權重法對研究區1994、2001、2008和2015年土壤鹽分含量進行空間插值，插值結果表明，研究區土壤鹽分含量平均值由1994年的1.006 1%增大為2015年的1.122 8%，除少數區域由于長期以來的灌溉洗鹽，土壤鹽分有所降低外，其他區域土壤鹽分總體上呈現增大趨勢，且土壤鹽漬化程度也越來越嚴重。

(3)根據研究區土壤鹽漬化發展趨勢，針對輕度鹽堿化區域，引進先進的灌溉技術如噴灌、滴灌等，加強灌溉水的利用率，控制地下水位的同時加速土壤表層脫鹽；針對中度鹽堿化區域，在淋溶洗鹽的基礎上，大量種植鹽生植物擴大地表植被覆蓋，降低土體與大氣間的水熱交換速率，以此來減小土壤表面積鹽的動力；針對重度鹽堿化區域，前期通過土地平整、農田水利工程等措施有效降低土壤內含鹽量，后期通過深耕曬垡結合有機肥、土壤翻耕結合生物菌肥等綜合措施改善土壤生態結構，降低土壤鹽堿化程度。

□