基于三明治估計的全國土壤墑情估計方法*

李淑華，郝星耀，周清波，高秉博，潘瑜春

(1.中國農業科學院農業資源與農業區劃研究所，北京 100081； 2.北京農業信息技術研究中心，北京 100097；3.國家農業信息化工程技術研究中心，北京 100097； 4.農業部農業信息技術重點實驗室，北京 100097；5.北京市農業物聯網工程技術研究中心，北京 100097)

?

·技術方法·

基于三明治估計的全國土壤墑情估計方法*

李淑華1， 3， 4， 5，郝星耀2， 3， 4， 5，周清波1※，高秉博2， 3， 4， 5，潘瑜春2， 3， 4， 5

(1.中國農業科學院農業資源與農業區劃研究所，北京 100081； 2.北京農業信息技術研究中心，北京 100097；3.國家農業信息化工程技術研究中心，北京 100097； 4.農業部農業信息技術重點實驗室，北京 100097；5.北京市農業物聯網工程技術研究中心，北京 100097)

土壤墑情監測對于科學灌溉、抗旱資源調配均意義重大，文章針對全國范圍墑情監測中由于監測點數據少、分布不均及數據點之間關系松散等導致的空間估計效果較差的問題，對全國范圍的土壤墑情空間分布規律與空間估計方法展開研究。采用適用于分層異質性表面進行空間估計的三明治估計方法，以全國縣域土壤墑情監測數據為采樣數據，以全國降水柵格數據與土壤分類數據作為分區數據，以縣行政區作為報告單元，實現土壤墑情監測點數據到面數據的上推。采用交叉驗證方法對三明治估計效果進行驗證，并與距離反比、克呂金、自然臨近點和樣條4種常規插值方法進行對比。結果表明：三明治估計的誤差均值為1.66%，誤差標準差為25.16%，墑情等級估計準確率為60.7%； 與其它對比方法相比，三明治估計的結果在誤差均值與誤差標準差上并沒有表現更好，三明治估計的誤差均值小于樣條插值方法但大于其它方法，誤差標準差與距離反比、克呂金及自然臨近點方法接近?？傮w來看，三明治估計能夠獲得與一般插值方法類似的估計效果，但是由于其更符合土壤墑情空間分布以及采樣數據情況，因此其結果更可靠和易于解釋，能更好地滿足全國土壤墑情監測的應用需要。

三明治估計 土壤墑情 空間插值 異質性表面 分區

0 引言

土壤水分的空間分布模式對于氣象、水土保持、農田保護和景觀管理等相關工作都具有重要的意義。在一個區域內進行土壤墑情監測，可以為農業生產中的科學灌溉提供決策基礎，在省或全國范圍內進行土壤墑情監測，對于水資源管理和抗旱政策制定及實施至關重要。但是受土壤性質、氣候條件等因素影響，土壤水分的時間空間變化很大，大面積連續監測比較困難[1]。

當前用于土壤水分測量的方法主要分為直接稱重與其他間接測量方法兩類，其中常用的間接測量方法有時域反射法(TDR)、頻域反射法(FDR)、探地雷達(GPR)及遙感(RS)等。在上述測量方法中，直接稱重、TDR、FDR方法是單點觀測，GPR與RS方法是區域觀測。探地雷達一天能夠觀測最多500m×500m的區域，在大范圍應用時經濟性不高[2]。遙感技術為大范圍地獲取土壤墑情相關參數提供了同步、時效性高且經濟的觀測方法，被認為是土壤墑情監測最為高效的一種方法。但土壤水分的遙感估計值容易受到地形與植被的影響，且特定的遙感技術與反演模型通常只能應用于一種特定的條件[3-5]，這些不足限制了遙感技術在大空間尺度下進行土壤墑情監測。土壤墑情的單點觀測值比區域觀測值更加精確，但是單點數據須采用空間估計方法上推到區域才能獲得更大空間范圍內的土壤墑情數據?？臻g插值是最常用的空間估計方法，其中在土壤墑情監測中應用最為廣泛的是克呂金(kriging)插值[6-8]。在較大的空間尺度下，土壤水分呈現較強的分區異質性，不能滿足普通克呂金插值方法的前提條件，因此其估計結果可靠性并不高。分區克呂金等改進方法則需要大量的采樣點，現實條件難以滿足。此外，土壤墑情常常需要按照不同的單元進行統計報告，簡單的分層采樣和推斷不能實現一次采樣多單元報告。

文章針對全國尺度的土壤墑情監測問題進行研究，從大空間尺度下的土壤水分空間分布特征出發，利用全國土壤墑情監測站點數據，采用三明治空間估計方法實現監測點數據到全國范圍的上推，并采用交叉驗證的方法對估計性能進行了檢驗。

1 數據與方法

1.1 土壤墑情數據

1.1.1 空間分布特點

土壤墑情由土壤中水分的輸入與輸出決定。對于農田土壤，水分輸入主要來源于灌溉與降水。灌溉對農田土壤墑情的影響范圍通常在一個較小的區域，在更大區域乃至全國范圍，降水是土壤水分輸入的主要來源[9]。土壤水分的主要輸出方式為蒸發與下滲， 2種方式的強度都極大地受到土壤質地的影響。因此在較大的空間范圍，降水與土壤類型形成了影響土壤墑情分布最為重要的兩個因素。

在全國范圍，較長時間內降水分布呈現很強的空間規律性，例如各地的年均降水量變化通常不大。但在較短的時期，降水分布則更多呈現斑塊狀，這就意味著即使空間分布相鄰的2個斑塊，其降水也可能存在很大差異。例如，某一縣域降水10mm，但其相鄰縣沒有降水，此時該相鄰縣之間的土壤墑情就會產生很大的差異。這種情況下采用一般的空間插值方法進行估計是非常困難的。土壤類型的空間分布同樣存在著斑塊性而非漸變或者同質，因此臨近監測單元之間的土壤類型與土壤水流動特性也可能存在很大差異。

根據前面的分析可知，土壤墑情的空間分布由降水與土壤類型的空間分布決定，在全國范圍較短時期內，其分布更接近于分層異質表面，即土壤墑情的分布表面能夠劃分為多個同質的分區。具有相同降水量和土壤類型的區域可以劃分為1個分區，在分區內土壤墑情基本一致，分區內的土壤墑情等級相同，并且能夠采用分區內采樣點的均值進行估計。劃分土壤墑情同質分區需要的包括降水與土壤類型等先驗知識通常都比較容易獲得。

空間差值是最常用的空間估計方法之一，分為確定性與地統計方法兩類。確定性插值方法(例如反距離插值)是基于與采樣信息之間的相似度，或者基于整個光滑表面來構造一個擬合表面。地統計方法(例如克呂金)在采樣點之間量化空間自相關，然后在預測點周圍建立采樣點的空間模型。采用非確定性的插值方法進行空間估計時，在采樣點處的估計值與測量值并不必然相等，有時甚至會差異較大。在較小的空間范圍內，這些差值方法對于土壤墑情分布的估計是有效的[10， 11]。在平穩隨機函數的假設之下，克呂金插值給出了中間值的線性無偏最優預測。但受土壤墑情空間分布的影響，在較大空間范圍內普通克呂金插值方法的估計性能會顯著下降。分區克呂金對于分層異質表面是有效的，但在每個同質分區內定義半變異函數與計算范圍都需要大量的采樣點(通常超過150個)，由于采集費用較高，通常采樣數很難達到要求，同時采樣過程耗時的增加對于監測與預測也非常不利。

1.1.2 數據采樣點

土壤墑情采樣數據來自農業部全國土壤墑情監測系統，采樣數據以縣為基本單元進行統計整理，數據采用烘干法測量，每月兩次，月值采用兩次測量的均值表示。每個縣的墑情數據通過計算若干個采樣點的均值獲得，這些采樣點代表了全縣的土壤墑情狀況。采樣點數據包含了4個深度的值： 20cm、40cm、60cm及80cm。為了降低影響因子的復雜度，該文僅采用20cm深度的土壤墑情數據進行計算，因為這個深度的土壤墑情對于降水與土壤類型的影響最為顯著，其它深度土壤墑情將在后續研究中分析。該文計算所用數據為2015年5月月均值數據，采樣數據分布在全國24個省市的112個不同的縣內，大部分西部非農產區沒有采樣分布(圖1)。

圖1 采樣數據分布

從圖2中可以看出，監測數據主要分布在東部地區，且覆蓋區域較小。在經典空間估計方法中，采樣樣本需要足夠大且采樣點分布應該為簡單隨機或系統抽樣，才能保證結果的可靠性。在土壤墑情監測中，采樣規模大小受經濟與時效性限制，采樣點的空間分布也很難保證隨機性或者系統性。一個縣是否能夠為全國范圍的墑情監測提供有效的數據很大程度上取決于其經濟與技術發展水平。因此，全國墑情監測系統中的采樣數據不能滿足經典空間估計的數據要求。

1.2 三明治估計方法

1.2.1 架構與特點

三明治估計方法是王勁峰研究員提出的一種用于估計分層異質表面并進行多單元報告的空間估計方法，該方法的估計模型由報告層、分區層及采樣層組成，因此被稱為三明治估計[11]。分區層是將研究區分為多個空間屬性同質分區，采樣單元隨機分布在各區中，分區的值可以通過采樣均值進行估計，分區采樣誤差能夠通過誤差傳導計算得到。采樣數據通過分區層傳遞到了具有很多報告單元的報告層，從采樣層到分區層再到報告層，信息流攜帶著估計均值與采樣誤差[12]逐層傳遞。三明治估計的核心思想包括分層采樣與從相同屬性類中借力。三明治估計框架提供了一種非常直接且相對簡單的方法來解決多報告單元問題與在多邊形系統之間傳輸數據[12]。該方法能夠提供與分區克呂金、分層貝葉斯估計同樣質量的估計，且更容易在地理信息環境中實現應用。在三明治估計中，采樣數據通過分區層傳導，而分區層與報告層之間是相互獨立的，因此采樣層中采樣單元不必遵循特定的空間分布。在分區層中，每個分區被認為是一個同質區域，在分區內采樣數據的分布對分區估計也沒有任何影響。因此，該文所采用的全國土壤墑情采樣數據的樣本數量與分布對于三明治估計是合格的。需要特別注意的是，如果某一分區中不包含任何采樣數據，那么這個分區就無法獲得估計值。

1.2.2 分區層和報告層

如何在分區層劃分分區是實現三明治估計的關鍵。根據土壤墑情空間分布特征與主要影響因素，結合數據獲取難易程度，該文采用降水與土壤類型作為分區因素[13， 14]，可保證得到的分區斑塊破碎程度適中，有利于估計計算。根據分區同質性的要求，具有相同土壤類型與相似降水量的區域可認為具有相同的墑情水平。

該文采用的土壤分類數據是由中國科學院土壤研究所繪制的全國土壤分類圖，進行分區計算的基本單元為土壤類型二級分類(圖2)。

圖2 土壤分類

圖3 降水量分級

降水數據從中國氣象數據共享網獲得，為了與土壤墑情監測數據匹配，采用2015年5月累積降水量，原始數據為0.5度的經緯格網柵格數據，通過分級并轉為矢量圖層后作為分區的基礎數據(圖3)。

土壤墑情監測主要的應用包括灌溉決策與抗旱資源部署，在這些工作中基本組織單元為行政區劃，因此該文報告層采用縣行政單元作為數據分析和評價的單元。

1.2.3 計算流程

實現三明治估計的數據處理流程主要包括3個步驟：采樣數據處理、生成分區單元及計算報告單元(圖4)。

圖4 土壤墑情三明治估計的數據處理流程

首先需要對采樣數據進行處理，為進一步的計算做準備。每個縣的土壤墑情值是來自于多個監測點，其均值代表了整個縣的均值，將這些值通過統一的行政區代碼關聯到全國縣行政區圖上，然后將行政區矢量圖轉為柵格，柵格單元值為土壤墑情值。

然后使用降水數據與土壤類型數據生成分區單元。原始的降水數據為柵格數據，單位為毫米。為了獲得降水同質區，將原始降水量分為5個等級，相同降水量等級的區域可認為降水基本相同。將分級的降水柵格數據轉為矢量，與土壤分類矢量圖層進行疊加，生成的新圖層分區層，分區層中每個矢量單元內具有形同的土壤類型和降水等級。

最后一步是通過分區統計方法計算報告單元值。將采樣數據傳輸到分區層的過程中，輸出數據采用柵格，因為柵格數據可以用作報告層計算的輸入數據。在數據從分區層傳輸到報告層的過程中，數據輸出格式采用統計表，以方便報告單元結果的統計。

整個數據處理過程可以采用ArcGIS中的分析工具來進行建模(圖5)，實現自動化的數據處理流程，其中標①圖框表示外部輸入數據，標②圖框表示數據處理過程，標③圖框表示數據處理結果，并且可作為后續數據處理過程的輸入數據。

圖5 數據處理流程建模

1.2.4 誤差評價方法

采樣數據通過2次分區統計傳導至報告層，三明治估計中的3層相互獨立，因此從采樣層傳輸到報告層的誤差能夠通過誤差傳導規則進行計算。

假設，采樣點數據值為ysi，i=1， 2，…，N，采樣數據方差為σs，分區層中每個分區的估計值與方差可表示為：

(1)

(2)

(3)

類似于分區單元的估計量，報告層中每個單元的估計量能夠用以下方程進行表達：

(4)

(5)

(6)

根據上述的推導，可以根據樣本數據計算報告層的方差。但是，在實際的土壤墑情估計中，更多關心的是一次估計的準確性，因此該文案例采用交叉驗證的方式來計算估計誤差，并采用誤差均值與誤差方差來對結果可靠性進行了評價。

2 結果與討論

2.1 三明治估計結果

使用降水數據和土壤類型數據計算的到三明治估計的分區層，結果見圖6。

圖6 生成分區單元的數據輸入輸出

按照三明治估計方法將采樣數據通過分區層傳遞到報告層，結果見圖7。

圖7 土壤墑情數據在不同層間的傳遞

在報告層中，共獲得2519個縣報告單元的估計值，有效報告單元的分布以及墑情級別與采樣數據的分布基本一致。采用交叉驗證方法對三明治估計結果準確性進行驗證，即從采樣數據中依次剔除1個縣域采樣數據后進行計算，計算未參與計算縣的估計值與采樣值之差作為誤差，最終得到誤差的均值為1.66，誤差標準差為25.16。在全國范圍的墑情監測中常使用土壤墑情等級而非土壤含水量，因此將交叉驗證結果進行分級后對比，估計值等級與采樣值等級相同的數量占全部采樣數據的比例為60.7%，等級差異小于等于1級的比例為76.9%。

從報告單元范圍來看，三明治估計有效地實現了對采樣覆蓋范圍內無值區域的估計，在采樣未覆蓋區域沒有估計值，符合空間估計規律。從誤差評價結果來看，誤差均值不為0表明估計結果存在一定程度的偏差。三明治估計也屬于統計學方法，因此采樣點的估計值與測量值之間存在一定差異，加之價值估計方法本身的誤差，導致等級預測準確率并不高，但等級差超過1級的僅為少數。全國范圍的墑情監測重點關注墑情的分布趨勢和數據覆蓋，因此三明治估計的總體效果基本能夠滿足監測需要。

2.2 對比分析

為了驗證三明治估計的效果，同時采用了距離反比、克呂金、自然臨近點和樣條4種常用插值方法進行估計計算，同樣采用交叉驗證方法來獲取誤差均值與誤差標準差兩個評價指標，誤差均值表示了系統誤差的大小，誤差標準差表示了估計的穩定性，結果見圖8。

圖8 不同估計方法的交叉驗證結果

從結果可以看出，與距離反比、克呂金、自然臨近點及樣條插值方法的結果相比，三明治估計的結果在誤差均值與誤差標準差上并沒有表現更好，三明治估計的誤差均值小于樣條插值方法，但大于其它方法，誤差標準差與距離反比、克呂金及自然臨近點接近?？傮w上，三明治估計方法并沒有明顯優于其它4種方法，但從空間估計的理論上分析，三明治估計更符合土壤墑情的分布規律與采樣數據的現實條件，同時三明治估計方法更加直接與易于理解，因此后續通過進一步改進三明治估計，應該能獲得明顯優于一般空間插值的估計結果。

3 結論

在全國范圍空間尺度下土壤墑情的分布表現為分層異質性，因此基于空間相關的空間插值方法就失去了利用采樣數據估計未知區域的理論基礎。另外受限于采樣代價采樣數據量常常不能滿足分區抽樣的要求，且數據分布不均勻，這就導致了在全國范圍進行土壤墑情監測非常困難。對于實現利用采樣數據來估計未知區域的目標，三明治估計是一種簡單可行的實現方式，其主要優勢在于：采樣數據通過分區層傳輸到報告層而非直接使用，這種機制提供了有效的從其它屬性相似單元借力的可能而不僅是空間相似區域，同質分區內通過簡單的加權平均方法就能實現； 三明治估計在不進行重采樣的條件下就能夠實現改變報告單元，在只能夠進行一次采樣，但需要進行不同報告單元數據的情況下非常經濟有效。

實例計算結果表明，在土壤墑情月值估計中三明治估計與常用的插值方法表現類似，但三明治估計在計算簡便性與多單元報告上的優勢明顯，同時三明治估計在實際應用中更容易被理解與接受。因此總體而言，三明治估計方法在全國范圍中長期土壤墑情估計中具有很大的應用潛力。

在先驗知識與數據的支持下，分區單元將估計限定到了一個相對穩定的分區內進行，但分區層的引入給估計過程帶來了兩個可預見的問題：首先，因為分區值是通過采樣數據均值進行估計，采樣數據中的極值有可能被平均，分區單元會消除這些特征值，盡管這些極值本身有時對于土壤墑情監測有著特殊意義； 其次，三明治估計的效果受分區質量的影響較大，土壤墑情是一個動態指標，因此分區也應該動態且與采樣匹配，另外氣象因子對土壤墑情影響可能存在的滯后效應也應被納入考慮以實現更合理的分區層。

[1] 劉明超， 秦其明，吳春雷，等.一種基于誤差分級的土壤水分數據同化方法.干旱區資源與環境, 2012，(11)： 139～144

[2] Dobriyal Pariva，Qureshi Ashi，Badola Ruchi，et al.A review of the methods available for estimating soil moisture and its implications for water resource management.Journal of Hydrology, 2012，(10)： 110～117

[3] Tao Yang，Huili Gong，Xiaojuan Li，et al.Progress of soil moisture monitoring by remote sensing.Acta Ecologica Sinica， 2010, 30(22): 6264～6277

[4] 鄧輝， 周清波.土壤水分遙感監測方法進展.中國農業資源與區劃, 2004，25(3)： 46～49

[5] 王利民， 劉佳，鄧輝，等.我國農業干旱遙感監測的現狀與展望.中國農業資源與區劃, 2008，29(6)： 4～8

[6] 王紅全， 李超.汾西灌區冬小麥田間土壤含水量空間分布研究.太原理工大學學報， 2009, 40(6): 636～639

[7] 楊紹輝， 王一鳴，孫凱，等.土壤墑情(旱情)監測與預測預報系統的設計與開發.中國農業大學學報， 2007, 12(4): 75～79

[8] Zhang Man，Li Minzan，Wang Weizhen，et al.Temporal and spatial variability of soil moisture based on WSN.Mathematical and Computer Modelling， 2013, 58(3): 826～833

[9] Tuttle Samuel E.，Salvucci Guido D.A new approach for validating satellite estimates of soil moisture using large-scale precipitation：Comparing AMSR-E products.Remote Sensing of Environment， 2014， 142： 207～222

[10]Nasta P.，Sica B.，Chirico G.B.，et al.Analysis of Near-surface Soil Moisture Spatial and Temporal Dynamics in an Experimental Catchment in Southern Italy.Procedia Environmental Sciences， 2013, 142(2): 188～197

[11]劉述強， 張長利，高君峰，等.土壤水分空間變異性的研究.東北農業大學學報， 2008, 39(8): 122～126

[12]Wang jin-Feng，Haining Robert，Liu Tie-Jun，et al.Sandwich spatial estimation for multi-unit reporting on a stratified heterogeneous surface.Environment and Planning， 2013， 45： 2515～2534

[13]Cho Eunsang，Choi Minha.Regional scale spatio-temporal variability of soil moisture and its relationship with meteorological factors over the Korean peninsula.Journal of Hydrology， 2014, 516(8): 317～329

[14]Letian Zhu，Yuanxin Liu，Feng Jiao.Time Series Analysis of Spatial Variability of Soil Moisture in Loess Hilly Region.Procedia Earth and Planetary Science, 2012，(5)： 346～353

NATIONWIDE SOIL MOISTURE ESTIMATION BASED ON SANDWICH SPATIAL ESTIMATION METHOD IN CHINA*

Li Shuhua1， 3， 4， 5，Hao Xingyao2， 3， 4， 5，Zhou Qingbo1※，Gao Bingbo2， 3， 4， 5，Pan Yuchun2， 3， 4， 5

(1.Institute of Agricultural Resources and Regional Planning，Chinese Academy of Agricultural Sciences，Beijing 100081，China；2.Beijing Research Center for Information Technology in Agriculture，Beijing 100097，China；3.National Engineering Research Center for Information Technology in Agriculture，Beijing 100097，China；4.Key Laboratory of Agri-informatics，Ministry of Agriculture，Beijing100097，China；5.Beijing Engineering Research Center of Agricultural Internet of Things，Beijing 100097，China)

Soil moisture monitoring is essential for scientific irrigation and drought resisting resource deployment. But in the current technology circumstance, it is difficult to survey soil moisture in large scale. The indirect estimation and derivation are main approaches to obtain soil moisture information. Remote sensing and spatial interpolation provide effective schemes to estimate soil moisture with affordable sampling quantity. But these two methods have some flaws in large scale application. Remote sensing estimation is susceptible to the effects of topography and vegetation. Further, specific remote sensing techniques and retrieval model can only be applied to a particular condition. On the other hand, ordinary interpolation methods are not applicable in large scale because the sampling points are insufficient and relationship between known points and unknown points is loose. In this paper, the spatial distribution characteristic of soil moisture in large scale was analyzed usingthe sandwich estimation method and the soil moisture surface was zoned into homogeneous subareas based on relevant knowledge and data. The soil type zones and precipitation data were used to partition the estimation region into some homogeneous zones, in which the soil moisture could be regarded asthe same. The reporting units were independent to the zoning and sampling layer, and could be changed without resampling or recalculate homogeneous zones. In order to evaluate the performance of sandwich estimation, cross validation method was used and compared with the other four empirical interpolation methods, including IDW, kriging, natural neighbor and spline. The results showed that the mean error was 1.66%, the standard deviation was 25.16%, the accuracy rate of moisture level estimation was 60.7%; the sandwich method effectively estimated the no-value area and was correspond with sampling data. The comparison results showed that the sandwich estimation had same performance asother methods. The mean error of sandwich estimation was less than that of the spine interpolation. The standard deviation of sandwich method was close to that of IDW, kriging, and natural neighbor. Overall, the sandwich estimation performed as well as other interpolation methods which frequently used in soil moisture monitoring. But the sandwich estimation was far superior for computational simplicity and multi-unit reporting. The sandwich estimation method can be better applied in the soil moisture estimation in nation-wide.

sandwich estimation； soil moisture； spatial interpolation； heterogeneous surface； zoning

10.7621/cjarrp.1005-9121.20161001

2015-10-08

李淑華(1979—)，女，河北衡水人，博士、助理研究員。研究方向：農業物聯網與云服務技術應用。

※通訊作者：周清波(1965—)，男，湖南沅江人，研究員、博士生導師。研究方向：農業遙感。Email：zhouqingbo@caas.cn

*資助項目：國家高技術研究發展計劃(863 計劃)“作物需水信息采集與智能控制灌溉技術”(2011AA100509)； 國家農業信息化工程技術研究中心開放課題“基于點對數頻率分布函數的空間變異函數構建方法研究”(KF2012N08-055)

S152.7

A

1005-9121[2016]10-0001-08