应用粒子滤波同化条件植被温度指数的土壤水分估测

解 毅，王鹏新，李 俐，荀 兰，张树誉

(1. 中国农业大学信息与电气工程学院，北京 100083； 2. 农业部农业灾害遥感重点实验室，北京 100083；3. 陕西气象局，陕西 西安 710014)

大面积的土壤水分监测是农业水管理、灌溉制度确定、农作物旱情预报以及农业增产的重要内容[1-2]。水文、气象站点仅能提供有限点的土壤水分信息，虽然精度较高，但不能从空间上反映土壤水分分布与变化规律，而遥感观测技术能够获取地表时空信息，达到全面监测土壤水分的目的[3-4]。

遥感监测土壤水分的方法主要包括：直接反演土壤水分，通过作物生长状况反映土壤水分，以及通过下垫面蒸散亏缺反映陆地生态系统总体水分状况[5]。直接反演土壤水分的方法适用于裸露地表，如有植被存在就会受到干扰[5]。利用遥感数据反演的植被指数能够反映作物生长状况，进而用于土壤水分的估测，然而，植被指数对土壤水分的响应有一定的滞后性。为了综合利用土壤和植被两方面信息监测土壤水分状况，Moran等[6]从能量平衡方程出发，提出了基于陆气温差-植被指数梯形特征空间的水分亏缺指数，用于判断陆面整体水分亏缺。王鹏新等[7]基于遥感反演归一化植被指数(NDVI)和地表温度(LST)的散点图呈三角形区域分布的特征，提出了条件植被温度指数(CVTI)的干旱监测方法，并用于区域干旱监测和预测。Sun等[8]将关中平原农田土壤浅层含水量和CVTI进行线性相关分析，结果表明，CVTI和0～20 cm层土壤含水量间的线性相关性显著。Patel等[9]的研究表明，CVTI和作物水分指数间的相关性显著。因此，CVTI适用于区域土壤水分的估测[10]。

基于遥感数据能够获取大面积的土壤水分信息，但仅能估测卫星过境时刻的瞬时状态，在时间上是不连续的。作物生长模型能够连续模拟单点尺度的土壤含水量，但当其应用到区域尺度时，模型参数的区域化方面存在困难，此外，模拟误差在模型运行过程中不断积累，导致模拟精度无法保证[11]。利用数据同化方法将作物生长模型和遥感数据相结合，能够发挥作物模型和遥感数据各自的优势，同时，在卫星过境时刻实时修正模型模拟的土壤水分，从而提高土壤水分的估测精度。陈鹤等[11]将集合卡尔曼滤波同化方法集成到水文强化陆面过程模型HELP(hydrologically-enhanced land process)中，对模型中的土壤水分和表面温度等状态变量进行优化，结果表明，数据同化方法能够提高土壤水分的模拟精度。Nagarajan等[12]应用粒子滤波算法同化田间实测土壤水分和SVAT-植被耦合模型，并模拟甜玉米生长季的根区土壤水分，结果表明，同化后的土壤水分模拟精度比模型直接模拟土壤水分的精度得到提高。

本文利用Landsat-8遥感数据反演CVTI，基于CVTI和实测土壤含水量间的线性相关性构建土壤水分反演模型。利用粒子滤波(PF)算法同化遥感反演的和CERES-Wheat模型模拟的土壤水分，得到以天为步长的土壤水分同化值，基于田间实测土壤水分分别检验模拟的、反演的土壤水分和同化的土壤水分的精度，分析数据同化方法对土壤水分估测精度的影响，以期获取准确的农田土壤水分信息，为大面积农业干旱监测和农田水分管理奠定基础。

1 土壤水分的反演

1.1 遥感数据预处理

研究区域关中平原覆盖了3个Landsat卫星轨道：126/036(关中平原东部)、127/036(关中平原中部)和128/036(关中平原西部)。本研究获取了2013—2015年关中平原冬小麦主要生育期(3月上旬至5月下旬，包括返青期、拔节期、抽穗～灌浆期和乳熟期)的Landsat-8遥感影像，所获取影像的云覆盖面积均小于5%。对Landsat-8影像进行辐射定标和大气校正，获取红光波段(第4波段)和近红外波段(第5波段)的反射率，以及热红外波段(第10、11波段)的辐射亮度。

1.2 地表温度的反演

根据Rozenstein等[13]提出的适用于Landsat-8数据的劈窗算法反演LST：

LST=A0+A1·BT10-A2·BT11

(1)

式中，BT10和BT11分别为第10、11波段的亮度温度(BT)，BTi=K2/ln1+K1/Li，i=10, 11，其中，Li为第i波段的辐射亮度，K1和K2为发射前预设的常量，在影像头文件中获取；A0、A1和A2是系数：

(2)

(3)

(4)

式中，ai和bi分别为第10、11波段根据不同温度范围确定的回归系数[13]；Ci、Di是由地表比辐射率(εi)和大气透射率(τi)所确定的参数：

Ci=εi·τi

(5)

Di=1-τi·1+(1-εi)·τi

(6)

根据覃志豪等[14]和宋挺等[15]对自然表面的比辐射率估计方法进行εi的估算：

εi=Pv·Rv·εiv+1-Pv·Rs·εis+dε

(7)

式中，εiv和εis分别为植被和裸土在第i波段的地表比辐射率，分别取ε10v=0.98672，ε11v=0.98990，ε10s=0.96767，ε11s=0.97790；在地表相对较平整的情况下，一般取dε=0；Pv为植被覆盖度；Rv和Rs分别为植被和裸土的温度比率。Pv的计算方法为：

Pv=(NDVI-NDVIs)/(NDVIv-NDVIs)

(8)

式中，NDVIv和NDVIs分别为完全植被和完全裸土覆盖像素的NDVI，取NDVIv为0.5，NDVIs为0.2。当NDVI>NDVIv时，取Pv=1；当NDVI

利用中纬度夏季大气模式估算Landsat-8 TIRS第10、11波段的大气透射率[13]，即τ10=1.0335-0.1134ω，τ11=1.0078-0.1546ω，其中，ω为大气水分含量。采用与Landsat-8影像获取日期相同且过境时刻接近的MODIS L1B Calibrated Radiances产品，通过其第2、19波段的比值计算大气水分含量[16]：

ω=α-lnρ19/ρ2/β2
(α=0.02;β=0.6321)

(9)

式中，ρ2和ρ19分别为MODIS数据第2和第19波段的表观反射率。

1.3 CVTI的反演

应用Landsat-8数据第4、5波段的反射率计算NDVI。基于NDVI、LST的散点图呈三角形区域分布的基础上计算CVTI：

(10)

其中

LSTmax(NDVIi)=a+b·NDVIi

(11)

LSTmin(NDVIi)=a′+b′·NDVIi

(12)

式中，LSTmax(NDVIi)和LSTmin(NDVIi)分别表示在研究区域内，当NDVIi等于某一特定值时所有像素的LST最大值和最小值，被称作热边界和冷边界；LST(NDVIi)表示像素的NDVI为NDVIi时的LST；a、b、a′和b′为待定系数，由研究区域的NDVI和LST的散点图近似获得。

1.4 土壤水分反演模型的构建

基于研究样点的地理坐标信息，获取样点所在像素的CVTI。将样点的实测0～20 cm层土壤含水量(θ)和反演CVTI进行线性相关分析，并构建CVTI和θ间的线性回归模型，进而反演区域土壤含水量。

2 土壤水分的同化

2.1 田间实测数据

在2013—2015年冬小麦生长季，课题组在陕西省关中平原选取了12～15个典型的冬小麦种植区域作为研究样点，包括灌溉样点和旱作样点(图1)。分别在冬小麦的播种期、拔节期、抽穗期以及成熟期实测和调查研究样点的田间数据，调查日期根据Landsat卫星的过境时间确定。

在冬小麦播种期实测样点土壤剖面各层的理化参数和土壤含水量。土壤剖面分为7层(0～12 cm、12～20 cm、20～50 cm、50～80 cm、80～120 cm、120～160 cm、160～200 cm)，采集土壤剖面各层的土壤样品，通过烘干法得到土壤剖面各层的体积含水量；在小麦拔节期和抽穗期，实测样点的土壤耕层(0～20 cm)含水量，同时调查作物田间管理数据和观测小麦生长状况；在小麦收获时测量土壤剖面各层的体积含水量，并实测籽粒单产。

2.2 CERES-Wheat模型

CERES-Wheat模型的输入参数包括气象数据、土壤参数、田间管理数据和作物遗传参数[17]，其中，气象数据通过分布在关中平原的气象站点观测获得；土壤参数和田间管理数据通过田间试验和地面调查得到。基于田间观测数据对冬小麦遗传参数进行标定，实现模型的“本地化”[18]，应用标定的CERES-Wheat模型模拟冬小麦整个生长季以天为步长的土壤水分，并将模拟土壤水分和遥感反演土壤水分进行同化，得到小麦主要生育期的土壤水分同化值。

图1 研究区域位置和典型样点分布Fig.1 Location of the research areas and distribution of the sampling sites

2.3 PF算法

(13)

(14)

再次进入预测阶段，重复上述过程，得到第k+2时刻的土壤水分同化值，依次重复，直至将所有的反演θ代入同化过程为止。

2.4 区域土壤水分估测

将研究样点的θ反演值和θ同化值进行线性相关分析，并构建小麦不同生育期的Landsat数据反演θ和同化θ间的线性回归模型。将遥感反演的区域土壤水分代入回归模型，估算关中平原冬小麦主要生育期的区域土壤含水量[19]。基于不同区域的累计降水量对估测土壤水分的空间分布进行分析。

3 结果与分析

3.1 CVTI和土壤水分的相关性分析

在冬小麦各主要生育期分别选取和田间试验日期相对应的Landsat遥感数据，将所有样点的遥感反演CVTI和田间实测θ进行线性回归分析，构建反演CVTI和实测θ间的线性回归模型。由于2013—2014年在126/036轨道未选取试验样点，以及2014—2015年在126/036轨道仅获取了1景云覆盖面积小于5%的影像且仅选取了3个试验样点，因此，以2013—2015年127/036和128/036轨道的分析结果为例(表1)，将不同生育期的反演CVTI和实测θ间的相关性进行对比。结果表明，在小麦主要生育期，反演CVTI和实测θ间的线性相关性显著。随着冬小麦进入返青期(3月15日)，至拔节期(3月24日)，CVTI和θ间的相关性逐渐增大；在拔节期和抽穗～灌浆期(4月28日)，CVTI和θ间的相关性达到极显著水平(拔节期：r=0.78，P<0.01；抽穗～灌浆期：r=0.79，P<0.01)；进入乳熟期后(5月11日和5月18日)，小麦叶片逐渐变黄和枯萎，CVTI和θ间的相关性逐渐减小，其中，5月18日的CVTI和θ间的相关性(r=0.65，P<0.05)相对较低。在冬小麦生长旺盛的拔节期和抽穗～灌浆期，CVTI和θ间的相关性大于其它生育期CVTI和θ间的相关性，说明CVTI适用于高植被覆盖条件下的土壤水分监测。

表1 CVTI和实测θ 间的线性回归分析

注： *表示0.05显著水平；**表示0.01显著水平。下同。

Note: * means significant at 0.05 level; ** means significant at 0.01 level. The same as below.

3.2 土壤水分的估测

3.2.1 土壤水分的同化 基于CVTI和实测土壤水分间的线性回归模型反演冬小麦不同生育期的土壤含水量。利用PF算法同化CERES-Wheat模型模拟的和遥感反演的θ，获取冬小麦主要生育期以天为步长的θ同化值，以2013—2014年灌溉样点扶风县杏林镇、武功县代家乡和旱作样点武功县武功镇、永寿县甘井镇的同化结果为例(图2)，将同化θ分别和模拟θ、反演θ进行对比。总体上，同化θ保持了CERES-Wheat模型模拟θ的变化趋势，同时，在遥感反演θ的影响下，θ同化值得到增大或减小，且更接近θ实测值。

在扶风县杏林镇播种后第161天，CERES-Wheat模型模拟的θ偏低，在遥感反演θ的修正下，θ同化值得到提高；在播种后第164天，θ同化值(θ=0.22 mm3·mm-3)比θ模拟值(θ=0.16 mm3·mm-3)更接近θ实测值(θ=0.23 mm3·mm-3)。同样，在武功县代家乡播种后第165天，在遥感反演θ的修正下，θ同化值比θ模拟值增大；在播种后第166天，θ同化值和θ实测值间的绝对误差(AE)(AE=0.025 mm3·mm-3)低于θ模拟值和θ实测值间的AE(AE=0.060 mm3·mm-3)，说明θ同化值的精度大于θ模拟值的精度，同化过程降低了CERES-Wheat模型模拟θ的累积误差。

分析所有样点的θ模拟值、θ反演值和θ同化值分别与θ实测值间的线性相关性，并计算其与θ实测值间的平均相对误差(MRE)和均方根误差(RMSE)(表2)。结果表明，θ同化值与θ实测值间的线性相关性(r=0.96，P<0.001)大于θ模拟值与θ实测值间的相关性(r=0.71，P<0.01)；θ同化值与θ实测值间的相关性大于θ反演值与θ实测值间的相关性(r=0.89，P<0.001)。θ同化值的RMSE和MRE比θ模拟值的RMSE和MRE分别降低了0.025 mm3·mm-3和2.70%，比θ反演值的RMSE和MRE分别降低了0.016 mm3·mm-3和4.15%。综上所述，基于PF算法同化过程能够充分结合遥感反演θ和CERES-Wheat模型模拟θ各自的优势，得到的θ同化值的精度高于θ模拟值和θ反演值的精度，因此，同化θ比模拟θ和反演θ更准确地反映了土壤水分的时间变化特征。

3.2.2 区域土壤水分的估测 将小麦不同生育期的Landsat数据反演θ和研究样点的θ同化值进行线性回归分析，得到反演θ和同化θ间的线性回归模型(以2014年3月15日的回归模型为例，代表返青期)：

PFθ1=1.28θ1-0.05R2=0.87P<0.001

(15)

式中，PFθ1为2014年3月15日的θ同化值；θ1为遥感反演θ。

表2 模拟θ、反演θ 和同化θ 与实测θ 间的RMSE、MRE和r

注：***代表0.001显著水平。

Note: *** means significant at the 0.001 level.

将遥感反演区域θ代入反演θ和同化θ间的线性回归模型，估测关中平原冬小麦不同生育期的区域土壤含水量，以2014年3月15日(返青期)、5月18日(乳熟期)和2015年4月28日(抽穗～灌浆期)的θ估测结果为例，基于不同区域的累计降水量对估测土壤含水量的空间分布进行分析(图3)。在2014年3月15日，凤翔县和乾县北部的土壤含水量低于岐山县东南部、扶风县西南部和眉县以北等区域(渭河以北区域)的土壤含水量。原因可能为，在2014年3月1日至3月15日期间，关中平原西部仅有两日出现降水，其中，凤翔县和乾县的累计降水量相对较低，分别为4.7 mm和6.4 mm；渭河以北的区域多为灌溉小麦地，虽然降水量较低，但3月10日左右的田间灌溉对土壤含水量影响较大，灌溉量一般为75～90 mm，灌溉区域的土壤含水量相对较高。

整体上，2015年4月28日的研究区域土壤含水量相对高于2014年3月15日和5月18日的土壤含水量。原因可能为，关中平原每年4月份的降水量较多，在2015年4月1日至4月28日期间，有10 d左右出现降水且降水量较大，例如，乾县4月1日至4月28日的累计降水量为106 mm，武功县为102 mm，临潼区为108 mm，蓝田县为114 mm；而在2014年5月1日至5月18日期间，关中平原西部仅有4 d左右出现降水且降水量相对较少，例如，凤翔县的累计降水量为29 mm，岐山县为26 mm，扶风县为25 mm。在2015年4月28日，白水县以及临潼区东部、临渭区南部和华县西南部等区域的土壤含水量高于其它区域的土壤含水量，原因可能为，白水县在4月26日和27日均出现降水，临潼区、临渭区和华县在4月26日出现降水，且临潼区、临渭区和华县在4月18日至28日的累计降水量分别为57.5 mm、49.8 mm和40.0 mm，降水量相对较高；而相邻的富平县、泾阳县和高陵县在4月26日和27日均无降水，其在4月18日至28日的累计降水量分别为20.1 mm、14.5 mm和24.7 mm，降水量相对较低。综上所述，不同研究区域的估测土壤含水量和累计降水量较吻合，准确反映了累计降水量对土壤水分的影响。

图2 2013—2014年冬小麦主要生育期的土壤水分同化曲线Fig. 2 The soil moisture assimilated curve in the main growing period of winter wheat during 2013 to 2014

图3 区域尺度土壤水分估测结果Fig.3 Estimated results of soil moisture at the regional scale

4 讨 论

土壤有效水分的减少会导致作物生长受限制，NDVI减小，通过计算同一时期NDVI和多年NDVI平均值间的差距，能够反映土壤水分状况，然而，NDVI对土壤水分的响应有一定的滞后性，不能用于农业干旱的实时监测。LST比NDVI更早地响应土壤水分的变化，在土壤水分降低时，植株叶片气孔关闭以减少蒸腾造成的水分损失，导致潜热通量减小，由于能量通量必须保持平衡，显热通量增大，LST增大。在NDVI、LST的散点图呈三角形区域分布的基础上反演得到CVTI，其不仅与NDVI的变化相关，而且在不同像素的NDVI相同时，与LST的变化相关。以往研究表明，CVTI是一种近实时的干旱监测方法，被广泛应用于干旱监测、预测和干旱影响评估等。本文将CVTI和0～20 cm层土壤含水量进行线性相关分析，在冬小麦主要生育期，CVTI和土壤含水量间的相关性显著，说明CVTI适用于高植被覆盖条件下的土壤水分监测。

CERES-Wheat模型能够连续模拟土壤含水量的时间变化信息[20]，然而，由于模型结构的不确定性，模拟误差在模型运行过程中不断积累，导致θ模拟值的精度降低。基于Landsat数据反演θ具有实时监测土壤含水量的优势，缺点为数量少，不能获取土壤水分的连续变化信息。利用数据同化算法将作物生长模型和遥感数据相结合，能够实现作物生长模型和遥感数据的优势互补。目前得到广泛应用的集合卡尔曼滤波(EnKF)同化算法，是在卡尔曼滤波基础上发展而来的，EnKF算法由于采用集合的思想近似地估计状态变量的概率密度分布，能用于非线性系统的数据同化，然而，EnKF算法仍没能摆脱状态变量服从高斯分布的假设，相比于卡尔曼滤波系列算法，粒子滤波算法不受模型状态变量和误差高斯分布假设的约束，适用于任意非线性非高斯分布的动态系统。本文基于PF同化算法，通过离散的遥感反演θ实时修正CERES-Wheat模型连续模拟的θ，得到以天为步长的θ同化值，且θ同化值的精度大于θ模拟值和θ反演值的精度，说明PF算法同化过程能够提高时间序列土壤含水量的估测精度。

由于Landsat卫星重访周期的限制和云对影像的干扰，难以获取频率较高的大面积土壤水分信息。MODIS-CVTI的空间分辨率较低，但其具有较高的时间分辨率。基于降尺度算法和MODIS-CVTI数据，得到空间尺度为30 m且周期较短的CVTI，进而反演土壤水分信息，将是今后研究工作的重点。

5 结 论

利用Landsat-8数据反演CVTI，基于CVTI和田间实测θ间的线性相关性构建土壤水分反演模型，并反演区域θ，通过PF算法同化反演的和CERES-Wheat模型模拟的θ，得到冬小麦主要生育期的θ同化值。主要结论为：从小麦的返青期至乳熟期，CVTI和实测θ间的线性相关性显著，尤其在小麦拔节期和抽穗-灌浆期，其相关性达到极显著水平，表明CVTI适用于高植被覆盖条件下的土壤水分监测。基于PF同化算法实现了遥感反演θ和CERES-Wheat模型模拟θ的优势互补，在反演θ的修正下，同化θ和实测θ间的线性相关性大于模拟θ和实测θ间的相关性以及反演θ和实测θ间的相关性，且同化θ的RMSE和MRE分别低于模拟θ和反演θ的RMSE和MRE，同化方法提高了θ的估测精度。此外，基于CVTI和PF算法估测的区域θ和累计降水量较吻合，反映了降水量对土壤水分的影响。因此，应用数据同化方法能够准确估测冬小麦不同生育期的土壤耕层含水量，为实现区域尺度作物干旱的实时监测提供了方法和依据。

参考文献：

[1] 吴春雷, 秦其明, 李梅, 等. 基于光谱特征空间的农田植被区土壤湿度遥感监测[J]. 农业工程学报, 2014, 30(16): 106-112.

[2] 周鹏, 丁建丽, 王飞, 等. 植被覆盖地表土壤水分遥感反演[J]. 遥感学报, 2010, 14(5): 959-973.

[3] 郑小坡, 孙越君, 秦其明, 等. 基于可见光-短波红外波谱反射率的裸土土壤含水量反演建模[J]. 光谱学与光谱分析, 2015, 35(8): 2113-2118.

[4] 胡猛, 冯起, 席海洋. 遥感技术监测干旱区土壤水分研究进展[J]. 土壤通报, 2013, 44(5): 1270-1275.

[5] 王敏政, 周广胜. 基于地面遥感信息与气温的夏玉米土壤水分估算方法[J]. 应用生态学报, 2016, 27(6): 1804-1810.

[6] Moran M S, Clarke T R, Inoue Y, et al. Estimating crop water deficit using the relation between surface-air temperature and spectral vegetation index.[J]. Remote Sensing of Environment, 1994, 49(3): 246-263.

[7] 王鹏新, 孙威. 基于植被指数和地表温度的干旱监测方法的对比分析[J]. 北京师范大学学报(自然科学版), 2007, 43(3): 319-323.

[8] Sun W, Wang P X, Zhang S Y, et al. Using the vegetation temperature condition index for time series drought occurrence monitoring in the Guanzhong Plain, PR China[J]. International Journal of Remote Sensing, 2008, 29(17-18): 5133-5144.

[9] Patel N R, Parida B R, Venus V, et al. Analysis of agricultural drought using vegetation temperature condition index (CVTI) from Terra/MODIS satellite data.[J]. Environmental Monitoringand Assessment, 2012, 184(12): 7153-7163.

[10] Peng J, Loew A, Zhang S Q, et al. Spatial downscaling of satellite soil moisture data using a vegetation temperature condition index[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 54(1): 558-566.

[11] 陈鹤, 杨大文, 刘钰, 等. 集合卡尔曼滤波数据同化方法改进土壤水分模拟效果[J]. 农业工程学报, 2016, 32(2): 99-104.

[12] Nagarajan K, Judge J, Graham W D, et al. Particlefilter-based assimilation algorithms for improved estimation of root-zone soil moisture under dynamic vegetation conditions[J]. Advances in Water Resources, 2011, 34(4): 433-447.

[13] Rozenstein O, Qin Z, Derimian Y, et al. Derivation of land surface temperature for Landsat-8 TIRS using a split window algorithm[J]. Sensors, 2014, 14(4): 5768-5780.

[14] 覃志豪, 李文娟, 徐斌, 等. 陆地卫星TM6波段范围内地表比辐射率的估计[J]. 国土资源遥感, 2004, 16(3): 28-32.

[15] 宋挺, 段峥, 刘军志, 等. Landsat 8数据地表温度反演算法对比[J]. 遥感学报, 2015, 19(3): 451-464.

[16] 胡德勇, 乔琨, 王兴玲, 等. 单窗算法结合Landsat8热红外数据反演地表温度[J]. 遥感学报, 2015, 19(6): 964-976.

[17] 王文佳, 冯浩, 宋献方. 基于DSSAT模型陕西杨凌不同降水年型冬小麦灌溉制度研究[J]. 干旱地区农业研究, 2013, 31(4): 1-10.

[18] 解毅, 王鹏新, 刘峻明, 等. 基于四维变分和集合卡尔曼滤波同化方法的冬小麦单产估测[J]. 农业工程学报, 2015, 31(1): 187-195.

[19] 解毅, 王鹏新, 王蕾, 等. 基于作物及遥感同化模型的小麦产量估测[J]. 农业工程学报, 2016, 32(20): 179-186.

[20] 赵海燕, 侯美亭, 王志伟. 利用CERES-Wheat模型分析冬小麦所需灌溉量的时空变化—以河南省为例[J]. 干旱地区农业研究, 2015, 33(4): 125-133.