基于三明治估计的全国土壤墒情估计方法*

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

(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%； 与其它对比方法相比，三明治估计的结果在误差均值与误差标准差上并没有表现更好，三明治估计的误差均值小于样条插值方法但大于其它方法，误差标准差与距离反比、克吕金及自然临近点方法接近。总体来看，三明治估计能够获得与一般插值方法类似的估计效果，但是由于其更符合土壤墒情空间分布以及采样数据情况，因此其结果更可靠和易于解释，能更好地满足全国土壤墒情监测的应用需要。

三明治估计 土壤墒情 空间插值 异质性表面 分区

0 引言

土壤水分的空间分布模式对于气象、水土保持、农田保护和景观管理等相关工作都具有重要的意义。在一个区域内进行土壤墒情监测，可以为农业生产中的科学灌溉提供决策基础，在省或全国范围内进行土壤墒情监测，对于水资源管理和抗旱政策制定及实施至关重要。但是受土壤性质、气候条件等因素影响，土壤水分的时间空间变化很大，大面积连续监测比较困难[1]。

当前用于土壤水分测量的方法主要分为直接称重与其他间接测量方法两类，其中常用的间接测量方法有时域反射法(TDR)、频域反射法(FDR)、探地雷达(GPR)及遥感(RS)等。在上述测量方法中，直接称重、TDR、FDR方法是单点观测，GPR与RS方法是区域观测。探地雷达一天能够观测最多500m×500m的区域，在大范围应用时经济性不高[2]。遥感技术为大范围地获取土壤墒情相关参数提供了同步、时效性高且经济的观测方法，被认为是土壤墒情监测最为高效的一种方法。但土壤水分的遥感估计值容易受到地形与植被的影响，且特定的遥感技术与反演模型通常只能应用于一种特定的条件[3-5]，这些不足限制了遥感技术在大空间尺度下进行土壤墒情监测。土壤墒情的单点观测值比区域观测值更加精确，但是单点数据须采用空间估计方法上推到区域才能获得更大空间范围内的土壤墒情数据。空间插值是最常用的空间估计方法，其中在土壤墒情监测中应用最为广泛的是克吕金(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 不同估计方法的交叉验证结果

从结果可以看出，与距离反比、克吕金、自然临近点及样条插值方法的结果相比，三明治估计的结果在误差均值与误差标准差上并没有表现更好，三明治估计的误差均值小于样条插值方法，但大于其它方法，误差标准差与距离反比、克吕金及自然临近点接近。总体上，三明治估计方法并没有明显优于其它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