基于MODIS的河南省春旱遥感监测

2018-07-04 07:51刘英，岳辉，李遥，鲁杨

干旱地区农业研究 2018年3期

刘英，岳辉，李遥，鲁杨

(西安科技大学测绘科学与技术学院，陕西西安 710054)

干旱的频发己成为我国重要的自然灾害之一。遥感技术具高时效性和低成本性等特点，在农业旱情监测中广泛应用[1]。目前，干旱遥感监测方法主要有热惯量法[2]、微波遥感法[3]和基于植被指数、温度指数的综合监测方法[4-8]。其中，温度植被干旱指数(temperature vegetation drought index，TVDI)在国内外应用较为广泛[9]。Carlson等[10]发现当植被覆盖和土壤湿度变化范围较大时，以归一化植被指数(normalized difference vegetation index，NDVI)为横轴和以地表辐射温度(land surface temperature, Ts)为纵轴的散点图特征空间呈三角形；Moran等[11]指出NDVI和Ts的散点图呈现梯形分布特征；Sandholt[12]根据Ts和NDVI的关系提出了TVDI的计算原理和公式。刘英等[13]指出NDVI-Ts特征空间则呈双抛物线型特征。刘公英等[14]利用MODIS NDVI和Ts数据，计算TVDI并利用降水数据和20 cm 土壤湿度数据进行结果验证，指出TVDI与实测土壤湿度拟合方程可决系数R2在0.598至 0.722之间，满足精度要求。熊世为等[15]利用HJ-1B星CCD和IRS数据，建立NDVI-Ts特征空间并计算TVDI，对2012年3月26日宿迁市土壤水分信息进行遥感监测，指出TVDI与10 cm深和20 cm深土壤湿度的相关性达到极显著水平。王纯枝等[16]采用MODIS 数据产品MOD13A2 和MOD11A2数据，构建Ts-NDVI特征空间，监测黄淮海平原2004年冬小麦旱情，指出TVDI与地表10～20 cm土壤水分可决系数最高，TVDI适合黄淮海平原作物生长季的旱情动态监测。

将地表温度和植被指数相结合的干旱监测方法不仅可以弥补基于单一指数的监测方法的不足，还可以提高干旱遥感监测的精确度和实用性。本文以河南省2000年、2005年、2010年和2015年2月下旬-5月下旬的旬时间尺度MODIS NDVI和Ts数据为基础，建立双抛物线型NDVI-Ts特征空间，采用TVDI研究河南省旱情的时空分布特征和变化规律，以期为河南省干旱灾害监测等提供基础数据，为有关部门制定政策提供科学的决策依据。

1 研究区概况

河南省位于北纬31°23′～36°22′，东经110°21′～116°39′之间，地形西高东低(图1),见222页彩图。河南省是一个农业大省，棉、粮、油等主要农产品产量均居全国前列，是全国重要的优质农产品生产基地。河南耕地面积为7179.2万hm2，其中小麦种植面积约占河南省总种植面积的37%。立春后，随着研究区冬小麦陆续进入返青期其需水量明显增加；而气温不断回升使得田间蒸发量不断加大，旱情极有可能迅速发展，对冬小麦返青生长和分蘖造成危害，进而危及夏粮丰收。

2 数据与方法

2.1 数据来源与处理

2000年、2005年、2010年和2015年2月26日到6月1日MODIS数据来源于美国宇航局NASA(http://revert. echo.nnsn.gov/revert/)。NDVI数据来自16d合成分辨率为250m的MOD13Q1数据集 (MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid)，Ts数据来自分辨率为1000 m的MOD11A2数据集(MODIS/Terra Land Surface Temperature/Emissivity 8-Day L3 Global/1 km SIN Grid)。16天数据选取时间段为：2000年、2005年、2010年和2015年的2月26日～3月13日、3月14日～3月29日、3月30日～4月14日、 4月15日～4月30日、5月1日～5月16日和5月17日～6月1日，涵盖冬小麦生长周期的冬灌～返青、返青～拔节、拔节～抽穗、抽穗、乳熟和成熟六个阶段。地面观测数据包括土壤相对湿度和降雨量数据，其中前者来自中国农作物生长发育和农田土壤湿度数据集，观测站点个数为17个(图1),见222页彩图；后者来自于中国逐日网格降水量实时分析系统(1.0版)数据集，两类数据下载网址为http://data.cma.cn/data/。

利用MRT (modis reprojection tool) 软件结合Cygwin，实现对2000～2015年MODIS NDVI和Ts数据的投影转换(将integerized sinusoidal grid 转换为Geographic 投影)和拼接。由于NDVI影像是16d合成的，Ts影像是8d合成的，为了保证时相一致，利用最大值合成法将两期8d合成Ts影像合并为一期16d合成影像，并将其重采样为250m。以0.01为步长，利用JAVA语言编程求取相同NDVI下对应的地表温度最大值和最小值并在Excel工具下制作得到干湿边散点图。将每日降水量栅格数据累加合成与MODIS数据时相相同的16d合成数据。借助ENVI 5.1和ArcGis 10.1，完成数据剪裁、公式计算以及旱情等级图的制作等。

2.2 双抛物线型特征空间的建立

Sandholt等[12]基于NDVI-Ts三角形或梯形特征空间，提出了TVDI，计算公式为：

(1)

式中，TVDI为温度植被干旱指数；Ts为任意像元地表温度(K)；Tsmax为相同NDVI的最小地表温度，对应NDVI-Ts特征空间的干边；Tsmin为相同NDVI的最小地表温度，对应NDVI-Ts特征空间的湿边；Tsmax和Tsmin可通过NDVI-Ts特征空间干、湿边模拟得到，即

Tsmax=a1+b1NDVITsmin=b2+b2NDVI

(2)

式中，a1、a2、b1、b2为干湿边方程系数，由NDVI-Ts特征空间散点图线性拟合得到；NDVI为归一化差异植被指数。

我国学者刘英等[17]将NDVI扩展到0.15以内，提出了NDVI-Ts双抛物线型特征空间及其TVDI计算方法，其干湿边方程为：

Tsmax=a1NDVI2+b1NDVI+c1

Tsmin=a2NDVI2+b2NDVI+c2

(3)

式中，c1、c2为干湿边方程拟合系数；NDVI为归一化植被指数。

TVDI的取值范围是0到1。TVDI越大，Ts越接近干边，越干旱；反之，TVDI 越小，Ts越接近湿边，越湿润。

3 旱情时空分布特征分析

3.1 干旱指数的构建和适用性

经处理得到NDVI-Ts特征空间干湿边散点图(图2)。由图2可知，NDVI-Ts特征空间干边Tsmax均为极大值向下开口的抛物线，干边可决系数基本接近0.9；湿边Tsmin为极小值向上开口的抛物线，干湿边呈明显双抛物线型特征。由表1进一步分析可知，干湿边拟合方程可决系数均通过了显著性检验(P≤0.05)。

3.2 TVDI与实测土壤湿度相关性分析

相关研究表明NDVI与表层10cm土壤相对湿度的相关性较好，与绿色植物的活力和密度关系密切[18-19]。以河南省17个气象站点观测的10cm深土壤湿度数据(soil moisture, SM)为横坐标，以双抛物线型NDVI-Ts特征空间反演得到的3×3像元的TVDI均值为纵坐标，构建SM-TVDI散点图(图3)，并计算其相关系数。

由图3可知，SM与TVDI总体呈现负相关线性关系，土壤湿度越大，TVDI越小，反之亦然。从实测土壤湿度和TVDI的线性拟合结果来看，拟合方程可决系数R2在0.3～0.5左右，均通过P≤0.05显著性检验，表明TVDI可作为干旱监测指标。

3.3 旱情遥感分析

采用TVDI作为干旱分级指标，将干旱划分为5级[16]，即：极湿润(TVDI为0～0.2)；湿润(TVDI为0.2～0.4)；正常(TVDI为0.4～0.6 )；干旱(TVDI为0.6～0.8 )；极干旱(TVDI为0.8～1)。因为研究区存在水体、建筑物及数据存在缺失情况，因此添加以下三种分类：无数据(TVDI为-1.5～1.49)；水体(TVDI为-1.49～1.2)；建设用地(TVDI为-1.2～0)。利用ArcGis制作河南省2000、2005、2010、2015年2月26日到6月01日旱情等级时空分布图(图4)。

图2 河南省2000年2月26日～6月1日 NDVI-Ts特征空间干、湿边散点图Fig.2 Scatter diagram of dry and wet edge in NDVI-Ts space of Henan from February 26 to June 1 in 2000

时间 Time干边(R2为拟合方程的可决系数)Dry edge (R2 is the decisive coefficient of the fitting equation)湿边(R2为拟合方程的可决系数)Wet edge (R2 is the decisive coefficient of the fitting equation)2005-02-26～2005-03-13Tsmax=-36.193NDVI2+19.604NDVI+298.6R2=0.9038,P≤0.05Tsmin=18.829NDVI2-6.2479NDVI+280.54R2=0.6463,P≤0.052005-03-14～2005-03-29Tsmax=-43.784NDVI2+30.12NDVI+298.91R2=0.9225,V≤0.05Tsmin=21.856NDVI2-17.632NDVI+279.35R2=0.5951,P≤0.052005-03-30～2005-04-14Tsmax=-49.713NDVI2+39.997NDVI+304.63R2=0.7914,P≤0.05Tsmin=14.901NDVI2-11.993NDVI+286.7R2=0.5099,P≤0.052005-04-15～2005-04-30Tsmax=-42.462NDVI2+30.703NDVI+306.67R2=0.9378,P≤0.05Tsmin=9.2086NDVI2-5.6821NDVI+291.59R2=0.4371,P≤0.052005-05-01～2005-05-16Tsmax=-37.707NDVI2+30.334NDVI+305.56R2=0.9295,P≤0.05Tsmin=21.82NDVI2-26.262NDVI+287.51R2=0.6857,P≤0.052005-05-17～2005-06-01Tsmax=-45.436NDVI2+31.361NDVI+314.22R2=0.8898,P≤0.05Tsmin=32.747NDVI2-27.05NDVI+289.64R2=0.6696,P≤0.052010-02-26～2010-03-13Tsmax=-37.69NDVI2+23.619NDVI+292.33R2=0.7011,P≤0.05Tsmin=52.369NDVI2-34.387NDVI+273.08R2=0.7745,P≤0.052010-03-14～2010-03-29Tsmax=-44.551NDVI2+29.063NDVI+298.4R2=0.918,P≤0.05Tsmin=26.586NDVI2-19.199NDVI+278.96R2=0.6617,P≤0.052010-03-30～2010-04-14Tsmax=-33.58NDVI2+22.204NDVI+305.25R2=0.774,P≤0.05Tsmin=7.1429NDVI2-2.0879NDVI+287.28R2=0.65,P≤0.052010-04-15～2010-04-30Tsmax=-33.522NDVI2+22.848NDVI+298.71R2=0.8367,P≤0.05Tsmin=4.2102NDVI2-1.6872NDVI+284.54R2=0.32082010-05-01～2010-05-16Tsmax=-28.946NDVI2+18.682NDVI+304.82R2=0.8426,P≤0.05Tsmin=14.462NDVI2-12.478NDVI+288.48R2=0.2977,P≤0.052010-05-17～2010-06-01Tsmax=-33.936NDVI2+24.526NDVI+309.14R2=0.9232,P≤0.05Tsmin=11.296NDVI2-14.937NDVI+291.39R2=0.7197,P≤0.052015-02-26～2015-03-13Tsmax=-43.419NDVI2+27.773NDVI+300.03R2=0.899,P≤0.05Tsmin=0.0717NDVI2+10.786NDVI+274.27R2=0.4313,P≤0.052015-03-14～2015-03-29Tsmax=-38.718NDVI2+26.753NDVI+302.43R2=0.9086,P≤0.05Tsmin=7.699NDVI2-3.9548NDVI+283.37R2=0.4255,P≤0.052015-03-30～2015-04-14Tsmax=-34.49NDVI2+24.773NDVI+299.83R2=0.8959,P≤0.05Tsmin=7.5417NDVI2-1.7161NDVI+281.23R2=0.6559,P≤0.052015-04-15～2015-04-30Tsmax=-26.491NDVI2+19.428NDVI+305.37R2=0.9404,P≤0.05Tsmin=16.194NDVI2-11.503NDVI+290.56R2=0.5822,P≤0.052015-05-01～2015-05-16Tsmax=-30.881NDVI2+21.194NDVI+304.87R2=0.8666,P≤0.05Tsmin=8.506NDVI2-6.6436NDVI+290.12R2=0.2371,P≤0.052015-05-17～2015-06-01Tsmax=-28.616NDVI2+22.155NDVI+308.23R2=0.853,P≤0.05Tsmin=4.1778NDVI2-12.399NDVI+291.99R2=0.608,P≤0.05

图3 2010年2月26日～6月1日10cm深度土壤湿度与TVDI的拟合(图中R2均为拟合方程的可决系数)Fig.3 Correlation between TVDI and soil moisture in 10cm depth from February 26 to June 1 in 2010

以2010年为例，结合2010年日降水栅格数据对2010年河南省旱情时空分布特征进行分析。

(a)2010-02-26～2010-03-13：旱情主要分布在西部、东南部和北部局地，此时段为冬小麦冬灌～返青阶段，这个阶段是决定成穗率高低的关键时期。此时间段内河南省累加降雨量为0～59.75mm，研究区西部、中部、北部降水量少，仅为0～5mm；降雨量主要集中在河南省南部区域。因此，旱情主要发生在西部、东南部和北部局地。

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

(b)2010-03-14～2010-03-29：西部和北部局地干旱加重，此时段为冬小麦返青～拔节阶段，此时植株生长快，需要大量水分和养料，但此时间段几乎无降雨，仅为0～3.7mm。为了促进冬小麦的生长，采取地面灌溉方式，使得研究区东部、南部及北部部分区域旱情有所缓解。

(d)2010-04-15～2010-04-30：全省旱情得到缓解，此时段为冬小麦抽穗阶段。2010-04-16日降水覆盖全省，研究区东部和南部2010-04-17、2010-04-20、2010-04-21日分别有降水，研究区西部2010-04-23、2010-04-28日分别降水，此时间段累加降雨量为9.15～60.54mm，但北部局地降雨量较小。

(e)2010-05-01～2010-05-16：全省旱情普遍加重，除东部局地外，其余地方旱情极为严重，此时段为冬小麦乳熟阶段。此时间全省绝大部分区域的降雨量为0.64-13.78mm，降水很快被蒸发和吸收，未能有效缓解旱情。

图4 2010年2月26日～6月1日河南省旱情等级时空变化图Fig.4 Temporal and spatial change map of Henan Province drought grade from February 26 to June 1 in 2010

(f)2010-05-17～2010-06-01：全省2010-05-28和2010-06-1日分别出现了均值为22.51mm和15.39mm的降水，旱情得到全面减轻，但西南部、西北部旱情仍然持续，此时段为冬小麦成熟阶段。

4 讨论与结论

文章基于双抛物线型NDVI-Ts特征空间，利用TVDI对河南省旱情等级时空变化进行分析，并结合同期降雨数据，探讨分析河南省旱情时空变化规律，主要结论为：

(1)利用MODIS NDVI和Ts数据，发现NDVI-Ts特征空间呈双抛物线型，干边方程可决系数R2在0.9以上；TVDI与实测土壤湿度呈现负相关关系，两者之间的线性拟合方程通过了P≤0.05的显著性检验；TVDI监测河南省旱情空间分布基本与降雨量空间分布一致。因而，基于双抛物线型NDVI-Ts特征空间的TVDI可以用于研究区旱情监测。

(2)由2000年、2005年、2010年和2015年的冬小麦旱情时空分布特征可知，研究区中南部冬小麦主产区土壤比较湿润，能够满足冬小麦需水量。

2.3 两组身长净增值比较干预组和对照组在身长净增值上的对比显示出生至42 d、42 d～2个月、2～3个月，3～6个月4个阶段，且干预组婴儿的身长净增值较对照组多，差异有统计学意义(P<0.05)。见表3。

利用气象站点实测土壤湿度与TVDI进行精度验证过程中，两者相关性较低，这是由于一方面TVDI是利用多天合成数据所得，而地面实测数据特定时间的数据，二者在时间尺度方面存有一定的误差；另一方面，气象站数据为点数据，而TVDI数据分辨率为250m，二者在空间上难以精确对应[20]。于敏等[21]、孔令寅等[22]、姜亚珍等[23]、李红军等[24]、闫峰等[25]研究也表明，NDVI/EVI-Ts特征空间方法结合了地表反射率和热红外辐射信息，能反映出一定时空尺度的土壤湿度状况，但容易受到植被覆盖状态和卫星观测质量的影响，导致特征空间中湿边的拟合度一般。因此，在后续工作中利用TVDI监测旱情需要从空间和时间尺度进行完善。

参考文献：

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

[2] Yu T, Tian G L. The application of thermal Inertia method the monitoring of soil moisture of North China plain based on NOAA-AVHRR data[J]. Journal of Remote Sensing,1997,1(1): 24-32.

[3] Sun L, Chen H W, Zhao L J, et al. The advances of drought monitoring by remote sensing [J]. Journal of Agro-Environment Science, 2004, 23(1): 202-206.

[4] 杨曦,武建军,闫峰,等.基于地表温度-植被指数特征空间的区域土壤干湿状况[J].生态学报,2009,29(3):1205-1216.

[5] 江东,王乃斌,杨小唤,等.植被指数-地面温度特征空间的生态学内涵及其应用[J].地理科学进展,2001,20(2):146-152.

[6] 齐述华,王长耀,牛铮.利用温度植被旱情指数(TVDI)进行全国旱情监测研究[J].遥感学报,2003,7(5):420-427.

[8] Wang C, Qi S, Niu Z, et al. Evaluating soil moisture status in China using temperature-vegetation dryness index (TVDI)[J]. Canadian Journal of Remote Sensing, 2004, 30(5): 671-679.

[9] 李春强,李红军. TVDI 在冬小麦春季干旱监测中的应用[J].遥感技术与应用,2008,23(2): 161-165.

[10] Carlson T N, Gillies R R, Perry E M. A method to make use of thermal infrared temperature and NDVI measurements to infer surface soil water content and fractional vegetation cover[J]. Remote Sensing Review,1994,9(1-2) :161-173.

[11] Moran M, Clarke T, 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.

[12] Sandholt I, Rasmussen K, Andersen J. A simple interpretation of the surface temperature/vegetation index space for assessment of surface moisture status[J]. Remote Sensing of Environment, 2002, 79(2):213-224.

[13] 刘英,马保东,吴立新,等.基于NDVI-ST双抛物线特征空间的冬小麦旱情遥感监测[J].农业机械学报,2012,43(5):55-63.

[14] 刘公英,申海凤,胡佳,等.基于TVDI 指数的冬小麦旱情动态研究—以河北省邢台市为例[J].干旱地区农业研究,2015,33(4):227-232.

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

[16] 王纯枝,毛留喜,何延波,等.温度植被干旱指数法(TVDI)在黄淮海平原土壤湿度反演中的应用研究[J].土壤通报,40(5):998-1005.

[17] Liu Y, Wu L X, Yue H. Bi-parabolic NDVI-Tsspace and soil moisture remote sensing in an arid and semi-arid area[J]. Canadian Journal of Remote Sensing, 2015, 41(3), 159-169.

[18] 张树誉,赵杰明,袁亚社,等. NOAA/AVHRR资料在陕西省干旱动态监测中的应用[J].中国农业气象,1998,19(05):27-29+33.

[19] 詹志明,冯兆东. 区域遥感土壤水分模型的方法初探[J].水土保持研究,2002,9(3):227-230.

[20] 孙振蓉. 京津冀地区冬小麦面积估算和旱情遥感监测研究[D]. 北京: 北京林业大学, 2015.

[21] 于敏,程明虎,刘辉.地表温度-归一化植被指数特征空间干旱监测方法的改进及应用研究[J].气象学报,2011,69(5):922-931.

[22] 孔令寅,延昊,鲍艳松,等.基于关键发育期的冬小麦长势遥感监测方法[J].中国农业气象,2012,33(3):424-430.

[23] 姜亚珍,张瑜洁,孙琛,等.基于TVDI河北省干热风同期土壤湿度监测研究[J].遥感技术与应用,2014,29(3):442-450.

[24] 李红军,雷玉平,李春强,等.地表温度-植被指数特征空间时空尺度效应分析[J].中国生态农业学报, 2014, 22(10):1252-1258.

[25] 闫峰,艳姣,吴波.多时间尺度温度-植被指数特征空间旱情监测的差异性[J].地理科学,2014,34(8):987-993.