基于理论干湿边与改进TVDI的麦田土壤水分估算研究

蔡庆空 陶亮亮 蒋瑞波 蒋金豹

(1.河南工程学院土木工程学院， 郑州 451191； 2.南京信息工程大学地理科学学院， 南京 210044；3.中国矿业大学(北京)地球科学与测绘工程学院， 北京 100083)

0 引言

土壤水分是地球系统的重要组成部分，影响着水圈、生物圈、大气圈的水热平衡，是农业、水文、生态、气候等领域的关键参数，对旱情监测及农作物生长发育有着十分重要的作用。

传统基于点的测量方法无法满足目前用于制作大面积、高精度土壤水分产品的应用需求，而遥感技术具有快速、大面积以及实时等优点，为估算大面积土壤水分时空信息、实时监测土壤水分动态变化提供了便利。光学遥感作为目前最常用、最成熟的土壤水分遥感监测方法而被广泛使用，其波段丰富、成像质量好、空间拓扑关系清晰。目前，利用光学遥感反演地表土壤水分含量最常用的监测方法包括归一化植被指数 (Normalized difference vegetation index，NDVI)[1]、表观热惯量 (Apparent thermal inertia，ATI)[2-3]、温度植被干旱指数 (Temperature vegetation drought index，TVDI)[4]、作物缺水指数 (Crop water stress index，CWSI)[5-6]以及高光谱遥感等。其中，TVDI方法通过将可见光波段与热红外波段结合构建特征空间来估算土壤水分含量，在土壤水分的监测中取得了较多的研究成果[7-12]。PATEL等[13]利用TVDI方法对印度地区的土壤水分进行估算，结果表明，TVDI与土壤水分之间存在非常好的负相关关系，尤其是在植被稀疏的情况下。ZHANG等[14]采用MODIS遥感产品数据构建地表温度及NDVI特征空间，对关中地区土壤水分空间变化进行研究，结果表明，TVDI与地表0～10 cm和10～20 cm土壤水分均有较高的相关性，因此TVDI可以反映并指示地表0～20 cm深度的土壤水分状况，可以用于关中地区干旱动态监测。

基于地表温度(Land surface temperature，LST)和植被指数(Vegetation index，VI)特征空间的TVDI在选取干湿边以及用于土壤水分监测时存在不足，理论上TVDI特征空间应包含3种植被覆盖类型，即裸土、部分植被覆盖以及全植被覆盖类型，但实际应用过程中，研究区域的覆盖状况往往无法完全满足植被覆盖条件，得到的干湿边实质上只是理论特征空间内部存在的边，而非理论干湿边。本文从地表能量平衡方程出发，根据能量交换过程中的极端状况，获取理论干湿边的端点方程，从而得到理论意义的特征空间，进而对试验区域内地表土壤水分进行反演研究，以期获得更高精度的土壤含水率估算值。

1 材料与方法

1.1 研究区概况

研究区位于关中平原中部地区，主要包括杨凌区、扶风县以及武功县(北纬34.1°～34.5°，东经107.8°～108.3°，图1)，地势平坦，气候温和，属于暖温带大陆性季风气候，年平均植被蒸发量为993.2 mm，平均日照时数2 163.8 h，主要种植作物为小麦和玉米。研究区内以水浇地(即农田)为主，部分区域有城镇用地、河流、山地等地表覆盖类型。

图1 研究区概况Fig.1 Overview of study area

以陕西杨凌区作为核心区，开展多次野外试验，主要获取地面土壤水分数据以及植株参数等，共设置3个小麦试验站点，每个试验站点根据农田观测面积设置不同的固定样点，同时以标杆标记并进行数据记录。观测时间为2014年3月29日和4月22日，观测数据主要包括地表土壤水分含量、地表温度、植株高度、叶面积指数以及光谱数据等。

所使用的气象科学数据从中国气象数据网(http:∥data.cma.cn/)下载，主要包括中国地面气候资料日值数据集和中国农作物生长发育和农田土壤湿度旬值数据集等。

1.2 卫星数据与预处理

中分辨率成像光谱辐射计(Moderate resolution imaging spectroradiometer，MODIS)共有36个离散波段，空间分辨率为250、500、1 000 m，扫描宽度为2 330 km。目前MODIS提供多种业务化运行产品，包括地表参量、海洋参量以及大气参量等从原始数据0级到应用模型开发5级及以上级别产品。本研究主要使用2014年3月29日(DOY088)和4月22日(DOY112)MYD09GA每日地表反射率数据以及MYD11A1每日地表温度数据(可从https:∥ladsweb.modaps.eosdis.nasa.gov/下载)。其中MYD09GA地表反射率数据共7个波段，空间分辨率为500 m。该产品为经过大气校正的2级地表反射率产品数据，可以得到星下点和标准太阳高度角，用于计算每日反射率及归一化植被指数(NDVI)。MYD11A1每日地表温度数据空间分辨率为1 000 m，利用MRT(MODIS reprojection tool)软件进行处理，投影类型为Albers Equal Area，并将影像分辨率从1 000 m重采样至500 m，用于计算空气温度、理论干湿边及观测干湿边。

1.3 研究方法

1.3.1温度植被干旱指数及Ts-Fv特征空间

研究发现地表温度与植被指数(VI)所构建的特征空间内存在多条土壤水分含量等值线，并且对于每条等值线上的点所对应的土壤水分含量均相等，该特征空间将地表覆盖区域分为裸土区、部分植被覆盖区以及全植被覆盖区，其所构建的干湿边与植被指数呈现显著的线性关系。基于此，SANDHOLT等[4]提出温度植被干旱指数(TVDI)，该指数与地表土壤水分密切相关，是表征植被受水分胁迫以及地表干湿状况的重要指标，其表达式为

TVDI=(Ts-Tmin)/(Tmax-Tmin)

(1)

其中

Tmax=a+bVI

(2)

Tmin=c+dVI

(3)

式中Ts——影像相应像元的地表温度

Tmin——Ts-VI特征空间的最低地表温度，即特征空间的“湿边”

Tmax——Ts-VI特征空间的最高地表温度，即特征空间的“干边”

VI——植被指数

a、b、c、d——经验参数，可线性拟合得到

目前最常用的植被指数为归一化植被指数，但是在植被高覆盖区域，NDVI易于饱和，同时其对于植被冠层背景影响较为敏感，使土壤水分估算不确定，因此本文将NDVI替换为改进植被覆盖度Fv[15]，其计算公式为

(4)

式中NDVIi——影像任一像元的NDVI

NDVImax——全植被覆盖对应的NDVI

NDVImin——全裸土覆盖对应的NDVI

构建的Ts-Fv特征空间(图2)不仅可以表示植被的生长状况，同时能够一定程度上避免TVDI特征空间估算时对于地表覆盖类型的限制，应用范围更广泛。

图2 理论干边与观测干边在梯形特征空间位置图Fig.2 Schematic of theoretical and observed dry edges in trapezoidal feature space

1.3.2地表能量平衡方程

无论是观测干湿边还是理论干湿边，实质上都受到地表能量平衡方程影响，而求导理论干边及理论湿边4个端点的方法均需要从地表能量平衡方程出发。地表能量平衡方程表示了太阳辐射能量除去被大气吸收以及被大气顶层反射回太空以外到达地表的能量守恒情况[16]。其可以表示为

Rn-G=H+LE

(5)

(6)

αs=0.160α1+0.291α2+0.243α3+0.116α4+
0.112α5+0.081α7-0.001 5

(7)

(8)

(9)

εs=0.261+0.314ε31+0.411ε32

(10)

G=Rn[Γv+(Γs-Γv)(1-fv)]

(11)

(12)

(13)

(14)

d0=0.65h

(15)

z0=0.13h

(16)

(17)

(18)

(19)

(20)

γ=0.665×10-3P

(21)

式中Rn——地表净辐射通量

G——地表土壤热通量

H——地表显热通量

LE——地表潜热通量

αs——地表反射率[17]

S0——太阳总辐射，利用6S辐射传输模型得到

σ——斯-玻常数，取5.67×10-8W/(m2·K4)

α1、α2、α3、α4、α5、α7——MODIS09反射率产品的1、2、3、4、5、7波段

εa——无云时大气有效发射率[18]

ea、Ta——参考温度水气压和近地表空气温度，可由半经验方程估算得到[19]

Tsky——天空等效温度[20]

εs——地表发射率[21]

ε31、ε32——MODIS产品第31、32波段发射率

Γv——植被覆盖下G和Rn比率，取0.05[22]

Γs——裸土覆盖下G和Rn比率，取0.315[22]

fv——植被覆盖度

ρ——空气密度，取1.29 kg/m3

CP——空气定压比热，取1.005 kJ/(kg·K)

ra——空气动力学阻抗[23]

z——地表参考高度，取2 m

k——卡曼常数，取0.41

u——参考高度2 m处的风速

d0——零平面位置高度[22]

z0——表面粗糙度[16]

h——植株平均高度

Δ——饱和水气压曲线斜率

e0(Ta)——饱和水气压[24]

dz——参考高度水汽饱和差

Tdmax、Tdmin——日最高气温、日最低气温

γ——干湿球常数P——大气压

1.3.3理论干湿边端点选取及方程推算

传统选取干湿边的方法即为选择植被指数最大和最小时地表温度的最大与最小值，构建的Ts-NDVI特征空间一般将湿边处理为与NDVI轴相平行的直线，干边处理为与NDVI呈线性关系，但是多数情况下，干边的最高温度往往随着NDVI的增大呈现先增大后减小的趋势，并且随着时间、季节以及纬度的变化而有所差异。因此利用地表能量平衡方程获取理论干湿边的端点显得尤为重要，理论湿边是由蒸腾和蒸发最大时裸土及植被覆盖下最低地表温度确定，理论干边则是由无蒸发和蒸腾时裸土和植被覆盖下最高地表温度确定，由此得到的干湿边具有明确的物理意义，更加适合植被覆盖下的水分估算研究。

在确定理论干边的端点时，地表处于没有土壤水分蒸发以及没有植被蒸腾的状态，此时的热量交换只有显热交换，能量平衡方程变为

Rn-G=H

(22)

对于裸土或植被分布较少的区域，如果仅仅选择NDVI最小时所对应的地表温度为理论最高温度，可能会出现理论最高温度低于实际温度的情况，因此选择地表温度随着NDVI增大而升到最高的点作为干边的端点，此时地表温度记为Tsd。植被完全覆盖情况下，随着NDVI升高，地表温度表现为稳定下降趋势，因此选择地表覆盖度fv趋近于1时地表温度最高值的点作为干边的另一个端点，地表温度记为Tvd，其计算公式分别为

Tsd=

(23)

Tvd=

(24)

式中Tave——地表平均温度

理论湿边上的点具有最大蒸发和蒸腾的特点，热交换形式为潜热交换，因此能量平衡方程可表示为

Rn-G=LE

(25)

在选择湿边端点时，对于裸土最大蒸发的点选择fv接近0时所对应的最低地表温度，记为Tsw，对于全植被覆盖下最大蒸腾的点则选择fv接近1时所对应的最低地表温度，记为Tvw，其计算公式分别为

(26)

(27)

1.4 评价指标

主要采用相关系数r和均方根误差RMSE 2个评价参数进行结果分析和验证，r越大，RMSE越小，表明土壤水分的估算精度越高，地表干湿状况监测

更精确。

另外为了更好地验证与评价结果，将计算得到的温度植被干旱指数转换为土壤水分估算值，其转换公式为

SMestimated=(1-TVDI)(SMmax-SMmin)+SMmin

(28)

式中SMmax——实测土壤水分含量的最大值

SMmin——实测土壤水分含量的最小值

2 结果与分析

2.1 特征空间构建及理论干湿边端点选取

结合MODIS遥感影像数据及地面观测数据，对计算得到的NDVI以0.01的步长提取对应地表温度影像上的最大值和最小值，同时对所得到的散点图进行线性拟合得到观测干湿边，如图3所示。

图3 观测、理论干湿边及其特征空间分布图Fig.3 Observation, theoretical wet and dry edges and their feature space

由图3可知，DOY112的观测干湿边上的点分布较为分散，尤其是在[0.3，0.7]区间内的点，而DOY088上的点分布较为集中均匀，同时研究区内当日的温度较高，最高地表温度达到307.44 K，最低地表温度为295.22 K，而DOY112的当日温度相对较低，最高地表温度为298.96 K，最低地表温度为290.16 K，另外DOY088及DOY112上观测干湿边线性拟合结果的相关系数均较高(表1)，同时图3也表示出研究区内理论干湿边的端点位置(黑色方形实点)及特征空间分布情况，其中绿色实线为理论干边，绿色虚线为理论湿边。由此可见，观测干湿边形成的特征空间仅为理论干湿边特征空间的一部分，根据理论干湿边进行干旱指数构建以及土壤水分估算将更加贴近实际情况。

表1 观测及理论干湿边相关参数Tab.1 Observation and theoretical parameters of wet and dry edges

2.2 土壤含水率估算结果与分析

根据式(1)和表1各干湿边方程分别计算DOY088和DOY112以观测干湿边以及理论干湿边为特征空间的温度植被干旱指数，分别记为TVDIob和TVDIth。本次试验DOY088和DOY112的样点数据分别为37、35个，利用其中的23、22个样点数据(约2/3)分别对计算得到的温度植被干旱指数进行拟合分析，其拟合结果如图4和表2所示。

图4 TVDIob和TVDIth与实测土壤含水率拟合结果Fig.4 Fitting results between TVDIob, TVDIth and measured soil moisture content

表2 TVDIob和TVDIth与实测土壤含水率拟合参数Tab.2 Fitting parameters between TVDIob, TVDIth and measured soil moisture content

图4a、4b为地表温度Ts与改进植被覆盖度Fv构建的特征空间所计算得到的温度植被干旱指数与实测土壤含水率拟合结果，而图4c、4d为地表温度Ts与归一化植被指数(NDVI)构建的特征空间所计算得到的温度植被干旱指数与实测土壤含水率拟合结果。从图中可以看出，不论是DOY088还是DOY112，利用观测干湿边及理论干湿边计算得到的温度植被干旱指数TVDIob和TVDIth与实测土壤含水率均有比较好的负相关性，同时TVDIth与实测土壤含水率之间的相关性比TVDIob高，相关系数均在-0.700左右，同时均方根误差也比TVDIob小，均不大于0.060 cm3/cm3。说明利用理论干湿边得到的温度植被干旱指数对土壤水分有着更高的敏感性，更适合植被覆盖区域的土壤水分监测，主要是因为理论干湿边更能体现水热平衡的实际状况。同时在DOY088和DOY112，Ts和Fv所计算得到的温度植被干旱指数均比由Ts和NDVI的计算值与实测土壤含水率的拟合效果更好，说明改进植被覆盖度的引入和使用能够拓宽干旱指数的应用范围，避免了传统特征空间干湿边估算时必须包含裸土、部分植被及全植被覆盖的地表植被覆盖状况的限制，为不完全包含所有的地表覆盖类型的区域及影像应用提供了技术支持。

另外，比较DOY088和DOY112的评价结果可以发现，DOY088时的拟合结果均比DOY112相应的结果好，DOY088时TVDIth与实测土壤含水率拟合结果相关性最高，相关系数达到-0.715，均方根误差也最小，为0.029 cm3/cm3。同时DOY112时散点分布相比较于DOY088分散得多，主要是因为小麦在4月处于生长状态，其长势旺盛，同时降水量远远达不到小麦的生理需求量，因此当时大部分小麦种植区域进行农业灌溉，而有些区域未能及时进行灌溉，使得DOY112时的土壤含水率的分布范围较大，其拟合精度小于DOY088。

利用另外独立的DOY088和DOY112时采集的14、13个样点数据分别对土壤含水率估算结果进行验证评价。利用式(28)将DOY088和DOY112时计算得到的TVDIth均转换为土壤含水率估算值。其验证结果如图5所示。从图5可以看出，由TVDIth计算得到的土壤含水率估算值与实测土壤含水率均有着较好的线性相关性，其相关系数r在DOY088时为0.805，在DOY112时为0.637，均方根误差分别为0.002、0.018 cm3/cm3。而对于DOY088，Ts-Fv的拟合结果比Ts-NDVI模式下的好，但是在DOY112时期其差距较小，主要是因为在DOY088时地表覆盖类型未能完全达到具备完全植被覆盖的情况，而DOY112时的地表覆盖情况对于温度植被干旱指数的计算更加适合，使得不管利用Fv还是NDVI，对土壤水分的估算结果影响不大。

图5 基于TVDIth计算的土壤含水率估算值与实测值验证结果Fig.5 Validation results between estimated soil moisture from TVDIth and measured soil moisture content

3 结论

(1)根据地表能量平衡方程，提出一种理论干湿边端点获取方法及以Ts-Fv特征空间构建的改进TVDI模型，拓宽了TVDI在干旱监测及土壤水分估算中的应用范围，引入改进植被覆盖度参数，在一定程度上避免了TVDI参数对地表覆盖类型的限制。

(2)理论干湿边所构建的特征空间更能反映自然现实的水热转换状况，更适合真实土壤水分的遥感估算。

(3)改进TVDIth与实测土壤含水率有着较好的线性相关性，相关系数在-0.700左右，均方根误差不大于0.060 cm3/cm3。DOY088和DOY112的土壤含水率估算值均与实测土壤含水率有着较好的拟合关系，尤其是DOY088的反演结果更加贴近实际地表干湿状况。

(4)理论干湿边端点的选取方法需要多种地表参数，使得该方法更加复杂，对地面观测数据的完备性要求较高。因此，在后续的工作中，将针对小麦不同阶段的物候状态及不同地表覆盖类型的研究区域开展研究和应用，进一步验证该方法的有效性。