基于HJ 卫星数据的土壤含水量反演及其旱情预测

熊世为， 李卫国， 贾天山， 胡姗姗， 景元书

(1.滁州市气象局，安徽 滁州239000;2.江苏省农业科学院农业经济与信息研究所，江苏 南京210014;3.南京信息工程大学应用气象学院，江苏 南京210044)

农业上将长期无雨或少雨，使土壤水分不足、作物水分平衡遭到破坏而减产的灾害定义为干旱。干旱是中国常发的灾害之一，严重影响社会经济发展。以2010 年西南旱灾为例，此次灾害共造成5.03 ×106hm2面积的农作物受灾，其中绝收面积达1.12 ×106hm2，直接经济损失达2.366 ×1010元［1］。苏北地区是中国重要的粮食生产基地，由于其特殊的地理、气候条件，冬春季常常遭受干旱影响，因此，实施及时的干旱监测对保护该地区农业生产以及发展精细化农业具有积极意义。传统的干旱监测是基于台站式的“点上”监测，不仅耗时、费力，而且难以获得灾情发展的“面上”信息，在时空上具有局限性。作为一种能够快速获取大面积地物光谱信息的技术，遥感在干旱监测方面已经得到了广泛应用［2-3］。目前常用的基于卫星遥感的干旱监测模式主要有热惯量法、植被指数法和地表温度法等，每种方法都存在各自的优缺点，如热惯量法物理机制明确，但仅适用于裸土或者植被生长早期［4］;植被指数法在高植被覆盖区具有较好的应用效果，却容易受其他环境因素干扰［5］;地表温度对土壤水分变化比较敏感，但反演过程复杂，且精度有待提高［6］。耦合多种模式的综合模型已成为该领域的发展趋势，其中NDVI-Ts特征空间法耦合了归一化植被指数(Normalized difference vegetation index，NDVI)和地表温度(Temperature of surface，Ts)，物理意义明确，该方法衍生出来的植被供水指数(Vegetation supply water index，VSWI)、作物水分亏缺指数(Crop water stress index，CWSI)、植被温度状态指数(Vegetation temperature condition index ，VTCI)、温度植被干旱指数(Temperature vegetation dryness index，TVDI)［7］都已经在干旱监测研究中得到了较好的应用。

HJ-1A/B 系列卫星是中国自主研发的用于环境与灾害监测的小卫星，兼具MODIS 数据重返周期短及TM 数据空间分辨率较高的优点，已广泛应用于作物长势监测［8］、产量估算［9］等方面，但在大范围、业务化监测(尤其是干旱监测)的应用仍显匮乏［10-11］。本研究通过构建基于HJ 卫星数据的NDVI-Ts 特征空间，提取温度植被干旱指数并结合实测数据将遥感指数转化为土壤相对湿度，并探讨利用HJ 数据的TVDI 法在苏北地区干旱监测中的适用性，以期为国产卫星数据在资源与环境监测及评估方面的深入应用奠定基础。

1 材料与方法

1.1 试验数据的获取

本试验研究区为江苏省宿迁市，介于33°80' ～34°25' N、117°56' ～19°10' E 之间，年均降水量为892.3 mm，受季风影响年际变化不大，但分布不均，易形成春旱、夏涝、秋冬干燥天气。本研究采用手持式GPS 定位仪在宿迁市区、沭阳县、泗阳县以及泗洪县各选取10 个具有代表性的样点并进行取土采样。样点的选定原则是大面积成片的田块，地表植被单一(调查期间种植小麦)，在距离道路、水体等至少100 m 的位置取土，样点间距5 km。调查期间发现部分田块已出现龟裂、小麦叶片萎蔫等旱象。

采用称质量法获取土壤含水量，即土壤水分占土壤干质量的百分率。用0.01g 精度的天平称取土样的质量，记作土样的湿质量(M)，在105 ℃的烘箱内将土样烘6 ～8 h 至恒质量，然后测定烘干土样，记作土样的干质量(Ms)，则土壤含水量(W)可表示为:

田间持水量室内测定采用环刀法。土壤相对湿度即土壤含水量与田间最大持水量的比值:

式中R为相对湿度，W为含水量，W0为田间最大持水量。

取土深度为0 ～10 cm、10 ～20 cm 和40 ～50 cm。为方便讨论，将各层次的含水量记为10 cm、20 cm、50 cm 处的土壤含水量。本研究共进行4 次地面取土试验，并选取4 景与之同步的HJ 卫星数据(表1)。

遥感影像预处理过程包括:首先采用校正过的Landsat/TM 影像对HJ 卫星的CCD 和IRS 影像进行几何精纠正，保证平均误差在1 个像元以内;再对CCD 和IRS 影像进行辐射定标，将灰度值图像转化为具有物理意义的辐亮度图像，定标公式及系数均来自影像头文件;最后采用具有较高波谱还原精度的中光谱分辨率大气辐射传输模式(MODTRAN)对CCD 影像进行大气校正，以获取真实地表反射率。

1.2 归一化植被指数(NDVI)和地表温度(Ts)提取

NDVI被定义为近红外波段与红光反射率之差与两者反射率之和的比值，表示为［12］:

式中，ρR、ρNIR分别表示为红光、近红外通道的反射率值。

由于IRS 数据只有一个热红外通道(10.5 ～12.5 μm)，与TM 数据相似，因此可以采用单窗算法进行地表温度的反演［13］:

式中，TS是地表温度，TB为地表辐射亮度对应的亮温，有效中心波长λe=11.511 μm［14］，常数C =1.438 768 69 ×10-2mK，ε 为地表比辐射率。其中，TB可以根据Plank 公式推导［15］:

式中，B为IRS 数据热红外通道辐射亮度，K1、K2均 为 常 数，对 于 IRS 数 据:K1= 579.20 W/(m2·sr·μm)，K2=1 245.58 K［13］。地表比辐射率ε 根据覃志豪等［16］提出的NDVI阈值法计算。反演的地表温度参照文献［17］进行详细的验证和正确性讨论。

1.3 NDVI-Ts 特征空间与温度植被干旱指数(TVDI)

土壤湿度与地表温度(TS)密切相关，根据地表能量守恒，土壤蒸发、冠层叶片蒸腾作用越小，带走潜热能量就越小，地表感热能量大，TS就高;反之，蒸发、蒸腾作用大，潜热增加，感热能量降低，导致TS低。而土壤湿度是影响土壤蒸发及冠层蒸散阻抗的重要因素。研究结果表明，在不同土壤表层含水量和地表覆盖条件下，TS和NDVI的散点图呈现出一种三角形空间，在此基础上定义TVDI为［17］:

式中，TS是特征空间内某给定像素的地表温度;NDVI为该像素的归一化植被指数;TSmax和TSmin分别表示该NDVI对应的干、湿边上的地表温度的最大值和最小值;a1、b1和a2、b2分别为干、湿边方程的斜率和截距，可通过线性拟合获取。因此TVDI是一个归一化值，介于0 ～1，其值越高表示土壤湿度越低，反之土壤湿度越高。

2 结果与分析

2.1 特征空间构建

根据上述方法进行归一化植被指数(NDVI)和地表温度(TS)计算，并构建NDVI-TS特征空间。由于NDVI＜0 对应地表类型为水体(可以认定含水量为100%)，在特征空间构建时不予讨论。因此，本试验在NDVI大于等于0 的范围内，以0.01 为步长，找出各NDVI值对应的最大地表温度TSmax和最小地表温度值TSmin，得到研究区各时相的NDVI-TS特征空间(图1)。可以看出，TSmax随着NDVI的增大(NDVI＞0)先有所升高，再逐渐降低，但总体上看是随着NDVI增大呈减小趋势;TSmin随NDVI的增大逐渐升高，但升高的幅度没有TSmax降低的幅度大。总体而言，NDVI与TSmax和TSmin基本呈线性关系(表2)，随着植被覆盖程度增加，地表温度最大值逐渐减小同时地表温度最小值逐渐升高并交汇于一点，也证明了TS与NDVI之间确实存在一种三角形的特征空间关系。其中TSmin的变化幅度虽然没有TSmax的大，但并不平行于NDVI值。不同时相特征空间的干、湿边的离散程度有所区别，分析原因，应该是不同NDVI值对应的象元量差别导致的。

2.2 特征空间干、湿边确定

前人在进行类似研究时，大都是直接将特征空间所有非水体象元(NDVI值为0 ～1)对应的地表温度最大值与最小值进行干、湿边拟合的。但在NDVI的某些范围内，其对应的地表温度最大值、最小值并不符合理论模型的变化趋势，如本研究中在NDVI处于0 ～0.2时，其对应的地表温度最大值并没有随NDVI的升高而降低。地表温度最大值应随NDVI增大而线性减小，这是建立在假设植被覆盖度与NDVI呈线性相关的基础上的，但实际情况并非如此。有研究结果表明，当植被覆盖度不足15%时，NDVI对植被覆盖度的指示作用很低［18］，也就是说，在植被覆盖低的区域，NDVI夸大了植被信息，因此本试验在进行特征空间干、湿边确定时，对NDVI的范围进行优化处理，选择NDVI≥0.2 的植被覆盖地像元，舍弃NDVI介于0 ～0.2 裸地岩石等低植被覆盖区域。该方法有别于常规的考虑所有非水体象元，这样做既合理舍弃了NDVI指示性低的区间，又有效提高了高植被区的反演精度。

表2 各时相特征空间干、湿边方程Table 2 Dry and wet edge equations of feature space of each time phase

2.3 TVDI 与实测土壤水分含量的关系

根据实测土壤含水量数据，分析TVDI与不同深度土壤相对湿度之间的关系。图2 反映了2012年3 月26 日各土层实测结果与提取的TVDI线性关系，可以看出，各层土壤相对湿度与TVDI均存在一定的线性负相关关系(其他时相结果类似，图略)，即随着TVDI指数增加，土壤相对湿度呈逐渐减小趋势，这与TVDI的原理一致。此外，图中数据点相对比较离散，经分析有两类因素导致误差，第一类是空间性因素，包含两点:(1)地面试验每个样方采集的五点并不能完全代表样方的平均水平;(2)几何校正误差导致像元与地面样方不能完全匹配。第二类是时间性误差，本研究地面调查数据都是两天内完成的，而遥感影像获取的是某一时刻的土壤水分信息，二者在时间尺度上不能完全匹配。进一步分析各层实测土壤相对湿度与TVDI指数的相关性，并进行显著性检验(表3)，可以看出，TVDI与每个时相各层实测土壤相对湿度的相关性大多达到显著或极显著相关的水平(2012-05-17 和2012-05-28 2 个时相的40 ～50 cm土层除外)，说明采用TVDI作为该区域土壤水分含量的监测指标是合理的。

图2 TVDI 与不同深度土壤相对湿度的关系Fig.2 The relationship between TVDI and relative soil water content in different depths

从不同土层深度上看，0 ～10 cm 深度土层实测含水量与TVDI的相关系数在-0.45 ～-0.60，10 ～20 cm 土层深度实测含水量与TVDI的相关性则更高，为-0.650 ～-0.800，而40 ～50 cm 土层相对湿度与TVDI的相关性最差。这是由于光谱辐射的反射面在浅层地表，可见光、近红外波段的辐射能量很难到达50 cm 的土层;由于大气、植被、人类活动等各种因素的影响，表层土壤(0 ～10 cm 土层)一般都是以疏松、充满间隙的物理形态出现，光谱透射率高，所以TVDI对10 ～20 cm 土层水分含量的指示性优于0 ～10 cm 土层。

2.4 旱情监测效果

尽管TVDI与实际土壤水分含量不能完全匹配，但二者的相关程度较高，尤其是与10 ～20 cm 土层相对湿度的相关系数达到0.8，因此，以TVDI作为浅层(10 ～20 cm)土壤水分含量的监测指标是可行的。根据前文10 ～20 cm 土层相对湿度与TVDI的拟合方程，将基于HJ 遥感影像的TVDI 图转化为相对湿度分布图，在此基础上根据中国气象局农业气候中心提出的旱情等级划分标准，得到研究区冬小麦浅层土壤(10 ～20 cm)旱情等级分布信息图(图3)。

4 个时相发生最重的干旱等级为中旱，虽然没有出现重旱，但地面试验期间发现部分麦田土壤出现干裂、植物叶片萎蔫的现象。统计各等级旱情占全市面积的比率(表4)，4 个时相中旱情发生最严重的是2012 年5 月28 日，该时相中旱情发生面积占全市总面积的16.57%，轻旱情、中旱情总的发生面积接近75.00%;其次为2012 年5 月17 日，该时相轻旱情、中旱情总发生面积接近50.00%，但中旱情发生面积比率仅为3.78%;旱情程度最轻的是2012 年4 月28 日，全区域没有干旱发生，且湿润地区超过60.00%;2012 年3 月26 日出现了轻度旱情，但发生面积较小，不足20.00%。

同期降水数据显示，宿迁市2012 年2 月～5 月降水比历年同期少70.00%，春旱旱情明显，虽然期间发生过几次降水过程，但主要集中在2 个时段内:第一个时段为2 月中旬至3 月中旬出现了低温连阴雨天气，这次降水过程虽然持续时间长，但是雨量不大，且前期的旱情较重，所以3 月26 日的遥感影像仍然监测出部分地区有轻度干旱发生;第二个时段为4 月下旬至5 月初，分别为4 月20 日～21 日、4月25 日、4 月29 日～5 月1 日，这次的寡照连阴雨天气持续时间长且降水量大，对应本试验4 月28 日反演结果全市绝大部分土壤呈湿润状态。之后天气晴好，农田水分蒸散导致干旱程度逐渐增强。

3 讨论

利用HJ 数据提取研究区的NDVI、TS影像，构建NDVI-TS特征空间并确定其干、湿边参数，分析了TVDI对不同土壤深度含水量的估算能力，主要结论如下:

表4 各等级旱情占全区面积的比率Table 4 The pecentages of areas with different soil moistures to the whole area

构建NDVI-Ts特征空间时发现当NDVI处于0 ～0.2 区间时，发现由于NDVI对植被覆盖程度低的地类指示性不足，导致地表温度最大值并没有如理论模型一般随NDVI增大而减小，因此，在进行干、湿边方程拟合时，只考虑NDVI＞0.2 的区间，这样做既合理舍弃了NDVI指示性低的区间，又有效提高了高植被区的反演精度。

TVDI与各层土壤实测含水量数据都有一定程度的相关性，尤其是与10 ～20 cm 土层相关系数为-0.65 ～-0.80，说明TVDI对于土壤含水量尤其是浅层土壤具有较高的指示意义。

根据TVDI与10 ～20 cm 土层相对湿度的线性关系，得到10 ～20 cm 土层相对湿度分布图，并进行旱情等级划分，同期降水数据证明反演结果是合理的。

TVDI耦合了植被信息与地表温度信息，在干旱监测领域有着较大优势。本试验初步探讨了利用HJ 数据TVDI法在苏北地区干旱监测中的适用性。然而干旱不仅与土壤含水量有关，还受土壤性质、植被类型、种植制度等因素影响，为了更准确地描述受旱程度，耦合多种要素建立更具机理性的监测指数是该领域的一个趋势。此外，本研究中选取的HJ-1B 影像的热红外通道空间分辨率为300 m，相较MODIS 数据的1 000 m 已有较大进步，但仍然达不到田块尺度要求，今后的研究可以尝试中、高分辨率数据融合来提高反演精度。

［1］ 张元元.应用FY-2 地表蒸散产品监测西南特大干旱［J］.气象，2011，37(8):999-1005.

［2］ 郭 巍，魏 明.2010 年春西南酷旱的监测与机理分析［J］.自然资源学报，2011，26(9):1628-1636.

［3］ 孙 丽，王 飞，吴 全.干旱遥感监测模型在中国冬小麦区的应用［J］.农业工程学报，2010，26 (1):243-249.

［4］ 李卫国.作物旱涝灾情遥感监测进展与思考［J］.江苏农业学报，2013，29(6):1503-1506.

［5］ 仝兆远，张万昌.土壤水分遥感监测的研究进展［J］.水土保持学报，2007，21(24):107-113.

［6］ ALDERFASI A，NIELSEN D C.Use of crop water stress index for monitoring water status and scheduling irrigation in wheat［J］.Agricultural Water Management，2001，47(1):69-75.

［7］ 熊世为，李卫国，景元书.基于HJ-1B 卫星遥感影像的小麦冠层温度的反演［J］.江苏农业学报，2012，28(6):1466-1470.

［8］ 李卫国，赵丽花.中高分辨率遥感影像在小麦监测中的比较［J］.江苏农业学报，2011，27(4):736-739.

［9］ 李卫国，王纪华，李存军，等.冬小麦花期生理形态指标与卫星遥感光谱特征的相关性分析［J］.麦类作物学报，2009，29(1):79-82.

［10］ 冯海霞，秦其明，蒋洪波，等.基于HJ-1A/1B CCD 数据的干旱监测［J］.农业工程学报，2011，27(增1):358-365.

［11］ WAN Z M，WANG P X，LI X W.Using MODIS land surface temperature and normalized difference vegetation index products for monitoring drought in the southern great plains，USA［J］.International Journal of Remote Sensing，2003，24:1-12.

［12］ 李卫国，赵春江，王纪华，等.基于卫星遥感的冬小麦拔节期长势监测［J］.麦类作物学报，2007，27(3):523-527.

［13］ 赵少华，秦其明，张 峰，等.基于环境减灾小卫星(HJ-1B)的地表温度单窗反演研究［J］.光谱学与光谱分析，2011，31(6):1152-1156.

［14］ JIMENEZ-MUNOZ J C，SOBRINO J A.A generalized singlechannel method for retrieving land surface temperature from remote sensing data［J］.Journal of Geophysical Research，2003，108:4688-4696.

［15］ 杨沈斌，赵小燕，申双和.基于Landsat TM/ETM+数据的北京城市热岛季节特征研究［J］.大气科学学报，2010，33(4):427-435.

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

［17］ 熊世为，景元书，李卫国.基于HJ-1B 遥感数据的冬小麦旱情监测研究［J］.麦类作物学报，2013，33(1):84-88.

［18］ 姚春生，张增祥，汪 潇.使用温度植被干旱指数(TVDI)法反演新疆土壤湿度［J］.遥感技术与应用，2004，19(6):473-478.