基于LSWI和NDVI时间序列的水田信息提取研究

刘 元 亮，李 艳,吴 剑 亮

(南京大学国际地球系统科学研究所，江苏省地理信息技术重点实验室，江苏 南京 210023)



基于LSWI和NDVI时间序列的水田信息提取研究

刘 元 亮，李 艳*,吴 剑 亮

(南京大学国际地球系统科学研究所，江苏省地理信息技术重点实验室，江苏 南京 210023)

利用多时相的Landsat TM/ETM+和DEM数据，分别采用直方图匹配(HM)和准不变特征点(PIFs)方法对影像序列进行相对辐射校正，减少了影像的光照和大气条件在时间上的不确定性，提高了归一化植被指数(NDVI)和地表水分指数(LSWI)的计算精度。根据水稻在不同生育期表现出的生理特征，基于LSWI和NDVI时间序列及高程特征，采用二叉树方法提取了浙江省金华市水田信息。经过验证，在空间上水田信息的提取精度达到92.3%，在县域尺度上提取面积与统计年鉴具有0.97的相关度。

水田；LSWI；NDVI时间序列；金华

0 引言

水稻是世界主要的粮食作物之一，水田面积约占世界农田总面积的11%，其灌溉用水约占全球淡水资源的70%，而且水田是温室气体甲烷排放变化的重要影响因素，也是碳汇的重要途径，因此，及时、准确地了解水田的面积及空间分布状况，对粮食安全、温室气体排放监测、水资源管理等方面具有重要意义[1-4]。目前，部分学者对水田信息的提取进行了相关研究：如杨艳昭等基于TM数据，采用水稻移栽期的地表水分指数(Land Surface Water Index，LSWI)等指标获取了吉泰盆地的水田分布信息[5]；钟仕全等基于HJ-1B卫星遥感数据，通过对HJ-1B星CCD数据水稻的光谱反射特性进行分析，建立了水稻遥感信息识别模型，应用决策树分类方法提取了广西宾阳县的水稻作物信息[6]；王力凡等利用水稻成熟期的CBERS-02B卫星遥感影像对比分析了水稻与其他地物的归一化植被指数(Normalized Difference Vegetation Index，NDVI)值和各个波段DN值的差异，并通过运算扩大数值差异来建立水稻像元提取模型，进而提取了南京市溧水县的水稻信息[7]。但以上研究均是面向像元的信息提取方法，容易产生由“同物异谱”或“同谱异物”引起的“椒盐效应”；另外，以上研究多是基于单一时相的遥感影像数据进行水田信息提取，由于水稻和其他植被同一季节的光谱特征类似而易造成误分现象。

本文以浙江省金华市的Landsat TM/ETM+和数字高程模型(DEM)为基本数据，利用面向对象的分类方法，定量提取了研究区的耕地信息，根据水稻在不同生育时期表现出的光谱特征，综合水稻的生长发育规律，基于LSWI和NDVI时间序列，在耕地信息的基础上进一步提取水田信息，为我国作物面积的信息提取提供了新的思路。

1 研究区和数据

金华市位于浙江省中部(东经119°14′～120°46′30″，北纬28°32′～29°41′)，地处金衢盆地东段，为浙中丘陵盆地地区，地势南北高、中间低。年平均气温为16.7℃(浦江)～18.2℃(永康)，年降水量为1 839.9 mm(永康)～2 052.3 mm(浦江)，四季分明，热量雨量丰富，干湿两季明显，属于亚热带季风气候。金华作为重要的粮食生产基地，素有“浙中粮仓”的美誉。2011年，金华市粮食作物总产量达90.7万t，目前，水稻仍是其主要种植作物。

本文选用的遥感影像由2008年、2009年和2010年4-11月份时相组合而成，空间分辨率为30 m，轨道号为119/40和118/40(http://www.gscloud.cn/)。由于Landsat 7 ETM+机载扫描校正器(SLC)发生了故障，造成2003年6月以后的图像出现条带而丢失了25%左右的数据，因此，在数据下载时利用地理空间数据云的在线条带修复模型(多影像局部自适应回归分析模型)进行了条带修复，最大限度减小了由条带缺失对信息提取结果造成的影响。

2 研究方法

本文首先基于决策树的分类方法定量提取研究区的耕地信息(即先区分原始影像中的植被与非植被地类，再区分植被地类中的林地与耕地)，然后利用特征时间段遥感影像的LSWI以及水稻生长期的NDVI时间序列，提取研究区水田的空间分布信息。

2.1 数据预处理

首先对遥感影像进行辐射定标和大气校正，将DN值转换为地表反射率；其次对各期遥感影像之间进行几何校正，使其误差在0.5个像元以内，以减小由影像之间的位置偏差对结果造成的影响；然后对遥感影像进行拼接裁剪，得到研究区的遥感影像；最后，因为遥感影像在不同季节的成像条件存在差异，从而会使各期遥感影像的成像效果有所不同，所以需要对各期影像进行相对辐射校正，确保各期影像之间的成像效果差异最小。研究中以7月份影像为参考影像，对同季相的8月、9月份的遥感影像分别采用直方图匹配的方法进行相对辐射校正[8]，对异季相的其他月份的遥感影像则分别采用“准”不变特征点(PIFs)法[9，10]，即选取影像之间相同地区地物类型没有变化的像元集合，构造两幅影像间灰度统计值的线性变换函数，进行相对辐射校正。图1展示了从各月份影像中选取的不变地物在相对辐射校正前后LSWI值和NDVI值与参考影像值的对比状况，数据点越靠近对角线说明与参考值越接近，即精确度越好。从图1可以看出，经过相对辐射校正，各月份同名不变地物点LSWI值和NDVI值与参考数据的差异得到缩小。

图1 相对辐射校正前后对比

Fig.1 Contrast before and after relative radiometric correction

2.2 研究区水稻信息物候历

旱地作物的光谱特征和水稻类似，会对水田信息的提取产生影响。基于遥感的水田信息提取主要依据水稻与其他旱地作物的物候期差异进行。本文利用水稻生长期影像的植被指数，通过其在不同生长期表现的物候特征对水田进行提取。通常，水田的动态变化一般包括灌水与秧苗移栽期、秧苗生长期、成熟收割期、休耕期。根据研究区农业气象观测数据可知水稻生长发育的时间分布特征，单季稻一般6月中下旬进行灌水移栽，10月上旬进入成熟收割期，其中8月为水稻的生长旺季；双季早稻则在4月下旬至5月上旬进行移栽，7月下旬至8月上旬成熟，其中6月份进入生长旺季；双季晚稻则在8月上旬灌水移栽，10月下旬收割，其中9月份生长最旺盛。

在水稻移栽期间，需要大量灌溉，其水深通常在2～15 cm，影像中水田表现出水体和秧苗的混合光谱特征，图2展示了典型水体、秧苗及移栽期含水体秧苗的光谱特征。数据分别采样自7月、8月、6月的样本。

旱地农作物(如棉花)不具备在第5波段反射率下降这一特征[11]。因此本研究首先采用对土壤湿度和植被含水量较为敏感的地表水分指数(LSWI)提取大部分水田信息，计算公式为：

LSWI=(ρnir-ρswir)/(ρnir+ρswir)

(1)

图2 典型地物光谱特征

Fig.2Thespectralcharacteristicsoftypicalobjects

式中:ρnir是近红外波段(TM或ETM+的第4波段)的反射率，ρswir则是短波红外波段(TM或ETM+的第5波段)的反射率。由于6月和8月分别是单季稻和双季晚稻的灌水移栽期，与生育期相同的旱地作物相比，地表湿度较大，此时相水田的LSWI值相对于这些旱地的LSWI值高，而在其他时相会由于植被含水量的原因使得旱地和水田不易区分，因此，可以通过设置6月份和8月份恰当的LSWI阈值将大部分的水田信息提取出来[5]。图3(见封2)是研究区耕地样本6月份和8月份LSWI变换图像，图3a中蓝色地物为单季稻移栽期水田，绿色地物为旺盛期双季早稻田，粉红色地物为旱地；图3b中蓝色地物为双季晚稻移栽期水田，绿色地物为旺盛期单季稻田，粉红色地物为旱地；图3c和图3d分别为图3a和图3b的LSWI变换图像，从图中可以看出，旱地的LSWI值较低，因此可以通过设置6月份和8月份的LSWI阈值剔除旱地信息，保留水田信息。但是，某些旱地作物在6月和8月生长较为茂盛，受植被含水量的影响其LSWI值较大；加之6月和8月恰逢我国南方多雨时节，降水会导致某些旱地表面被雨水覆盖，表现出和水田类似的光谱特征，LSWI值相对于普通旱地作物会偏高。因此，基于LSWI提取的水田可能含有部分旱地信息。

NDVI时间序列能够反映植被随季节变化的规律，由于不同植被在各季节表现出的生理特征不同，其植被指数时间序列差异较大，因此，利用植被指数时间序列进行植被信息提取，能够区分同一时段内无法区分的植被，减少因为季节变化而造成的误分现象[12]。因此,本文在利用LSWI提取水田的基础上，又根据水稻在不同生育期表现出的生理特征，利用NDVI在整个水稻生育期的变化状况，进一步提取水稻信息。NDVI的计算公式为：

NDVI=(ρnir-ρred)/(ρnir+ρred)

(2)

式中:ρnir是近红外波段(TM或ETM+的第4波段)的反射率值，ρred是红光波段(TM或ETM+的第3波段)的反射率值。

根据水稻在各个生育期的物候特征可知，单季稻在6月份进行移栽，由于地表水体的影响，NDVI值较低，随后稻苗开始返青，NDVI值逐渐升高，到8月份进入生长旺季，NDVI值达最大，随后进入成熟收割期，NDVI值开始下降；双季早稻则在4月、5月份进行移栽，NDVI值较低，随后开始升高，到6月达最大；双季晚稻则在7月下旬进行移栽，9月进入生长旺季，NDVI值达最大。根据实测数据和谷歌地球高清数据，各选取6个双季稻和单季稻样本点，分析相对辐射校正前后水稻在其生育期内NDVI值的动态变化，如图4和图5所示。

由图4可以看出，相对辐射校正后双季稻在生育期内(5-10月份)具有两个NDVI峰值，分别出现在6月份和9月份，单季稻在生育期内(6-10月份)具有单个峰值，出现在8月份，校正后水稻的NDVI值的动态变化状况可以基本反映水稻在生育期的生长状况，从而可以作为区分水田和旱地的依据。

图4 单季稻(左)和双季稻(右)NDVI曲线(校正后)

Fig.4 The NDVI value of single cropping rice(left) and double cropping rice(right) in different months (after correction)

图5 单季稻(左)和双季稻(右)NDVI曲线(校正前)

Fig.5 The NDVI value of single cropping rice(left) and double cropping rice(right) in different months (before correction)

由图5可以看出，相对辐射校正前双季稻在生育期内虽然也具有两个NDVI峰值，但在其他月份的NDVI值变化不尽一致，单季稻NDVI峰值出现的时间具有不确定性，因此，相对辐射校正前水稻的NDVI值随时间的变化与实际情况不相符，由此提取的水田信息必然存在较大偏差。

2.3 分类方法

利用eCognition软件，将分割尺度设为30，形状参数设为0.15，紧凑度设为0.5，对相对辐射校正后的遥感影像进行分割，得到含有语义信息的影像对象。采用二叉树分类体系:二叉树第一层两个叶片为植被和非植被，根据7月份NDVI值，通过设定恰当的阈值T1将植被提取出来。二叉树第二层叶片为林地和耕地，其光谱特征类似，但是植被指数随时间变化的规律差异较大，林地的NDVI值长期较高且保持稳定，而耕地的NDVI值变化相对较大，因此，利用5-10月份的NDVI值大于某一数值T2可剔除林地；此外，耕地一般分布在海拔较低的平原或山谷地带，所以通过设定恰当的高程阈值T3可进一步提取耕地信息。经过试验，用于分类提取的阈值最终设置如下：T1为0.25，T2为0.65，T3为330。经过影像分割、分类规则集构建、多阈值分类及后处理等过程，提取耕地信息(图6a)。二叉树第三层叶片为水田和旱地，以提取得到的耕地信息为基础，将LSWI的阈值设置为T4=0.06，T5=0.08，提取规则如下:

其中：g(x,y)=1代表此对象为水田，g(x,y)=0则代表此对象为旱地，LSWI06和LSWI08分别为6月份和8月份的LSWI。

根据以上规则提取了混有旱地的水田，然后计算提取的影像对象在各个时相的NDVI值。通过对单、双季稻样点各时相NDVI统计分析，采用以下规则提取水田信息(图6b)：1)在水稻生长期内NDVI有两个峰值，分别出现在6月份和9月份；2)在水稻生长期内NDVI有一个峰值，出现在8月份。本文主要采用基于对象水平的二次差分方法对单峰值和双峰值进行识别，具体过程如下，

D1i=NDVIi+1-NDVIi

(3)

D2i=1ifD1i>0; D2i=-1ifD1i≤0

(4)

D3i=D2i+1-D2i

(5)

根据式(3)计算每个对象前后两个时相的NDVI差值，得到D1序列；根据式(4)判断规则得到D2序列；最后由式(5)计算D2前后两元素之差，得到D3序列。在D3序列中，若有两个元素的值为-2，且位置分别在2和5，则表示在水稻生育期有两个峰值且分别在6月份和9月份；若只有一个元素的值为2，且位置为4，则表示在水稻生育期内有一个峰值且出现在8月份。

图6 耕地、水田和旱地信息提取

Fig.6 The spatial distribution of cultivated land,paddy field and dry land

3 结果与分析

基于以上水稻识别规则提取了水稻信息，经过适当的修正处理得到了研究区水稻信息的空间分布状况(图6b)。结合实验区野外实测数据和谷歌地球高清数据目视判读相结合的方法，对研究区选取的样点进行相关指标计算，对耕地和水田提取结果进行了精度分析和评价。从表1和表2可见，耕地信息提取的总体精度达93.7%，Kappa系数为0.874；水田信息提取的总体精度达92.3%，Kappa系数为0.84。此外，还用相同方法对相对辐射校正前影像信息提取的结果进行了精度验证，发现校正前耕地信息的总体精度为86.9%，Kappa系数为0.737；水田信息提取的总体精度为85.6%，Kappa系数为0.712。

表1 耕地精度评价

Table 1 Accuracy evaluation of cultivated land

地物类别验证样本耕地其他总计用户精度(%)耕地其他总计生产精度(%)总精度(%)102510795．393．789210090Kappa系数110972070．87492．794．8

表2 水田精度评价

Table 2 Accuracy evaluation of paddy land

地物类别验证样本水田其他总计用户精度(%)水田其他总计生产精度(%)总精度(%)97710493．392．399610591．4Kappa系数1061032090．8491．593．2

分别对研究区各个县区进行水田面积验证，由于研究区缺少2010年土地利用面积统计年鉴，但土地利用的年际变化较小，因此本文利用2007年统计年鉴的土地利用面积数据与本研究结果进行对比。由统计年鉴可知，研究区水田总面积为1 420.41 km2,本文提取的水田总面积为1 404.91 km2，面积一致性达R=0.989。在县区尺度上，由统计年鉴得到的水田面积和本实验提取的水田面积具有较好的相关性(图7)。由图6b可知，金华地区水田呈现西多东少的态势，主要集中在河流两侧的河谷地带及水库周边等灌溉条件较好的地区，在灌溉条件较差的山地丘陵地区，水田分布较少。

图7 面积验证

Fig.7 The validation of area in county scale

4 结论

本文基于LSWI和NDVI时间序列，利用2010年前后水稻生长期的多时相遥感影像数据，依据水稻在各个生长期与其他作物表现出的生理差异特征，通过一定的识别规则提取了浙江省金华市水田的空间分布信息，得出以下结论：1)6月和8月LSWI值较高；双季稻NDVI峰值出现在6月和9月，单季稻NDVI峰值出现在8月。2)相对辐射校正使得耕地及水田的提取总体精度分别比校正前提高了6.8%和6.7%。3)二叉树水稻提取方法综合利用了遥感及空间数据在时间、空间两方面的信息优势，获得了92.3%的精度，在县域尺度上与统计年鉴具有0.97的相关度。

但是，本研究中未对未耕种水田信息进行处理,数据成像时间受TM数据的时间分辨率和天气状况影响，不容易获取满足作物关键生长期的影像，今后可以考虑利用更高时间分辨率数据进一步提高作物信息提取的精度。

本文工作受软件新技术与产业化协同创新中心部分资助，此致谢忱！

[1] XIAO X M,STEPHEN B,STEVE F,et al.Mapping paddy rice agriculture in south and southeast Asia using multi-temporal MODIS images[J].Remote Sensing of Environment,2006,100:95-113.

[2] 张友水,原立峰,姚永慧.多时相水田信息提取研究[J],遥感学报,2007,11(2):282-288.

[3] DENIER VAN DER GON H.Changes in CH4 emission from rice fields from 1960 to 1990s:Impacts of modern rice technology[J].Global Biogeochemical Cycles,2000,14(1):61-72.

[4] 郑长春,王秀珍,黄敬峰.多时相MODIS影像的浙江省水稻种植面积信息提取方法研究[J].浙江大学学报(农业与生命科学版),2009,35(1):98-104.

[5] 杨艳昭,郑亚南,姜鲁光,等.吉泰盆地水稻熟制遥感识别研究[J].地理与地理信息科学,2014,30(3):16-20.

[6] 钟仕全,莫建飞,陈燕丽,等.基于HJ-1B卫星遥感数据的水稻识别技术研究[J].遥感技术与应用,2010,25(4):464-468.

[7] 王力凡,潘剑君.基于CBERS-02B卫星影像光谱系想你的水稻种植面积提取方法——以南京市溧水县为例[J].南京农业大学学报,2013,36(1):87-91.

[8] 潘志强,顾行发,刘国栋,等.基于探元直方图匹配的CBERS-01星CCD 数据相对辐射校正方法[J].武汉大学学报(信息科学版),2005,30(10):925-927.

[9] 张友水,冯学智,周成虎.多时相TM影像相对辐射校正研究[J].测绘学报,2006,35(2):122-127.

[10] 张鹏强,余旭初,刘 智,等.多时相遥感图像相对辐射校正[J].遥感学报,2006,10(3):339-344.

[11] 魏新彩.基于HJ卫星影像的水稻种植面积遥感信息提取方法研究[D].武汉:湖北大学,2013.

[12] 朱满,胡光宇,于之峰.基于融合NDVI和EVI时间序列的遥感影像分类研究[J].遥感信息,2009(5):44-46.

Study on Extraction of Paddy Fields Based on LSWI and Time-Series NDVI

LIU Yuan-liang,LI Yan,WU Jian-liang

(InternationalInstituteforEarthSystemScience,NanjingUniversity,JiangsuProvincialKeyLaboratoryofGeographicInformationScienceandTechnology,Nanjing210023,China)

Rice is one of the most important food crops in China.Timely and accurately acquiring the area and spatial distribution of paddy fields is significant to the crop yield estimates,food security and water resources management and so on.In this paper,the multi-temporal Landsat TM/ETM+ images and DEM are used and histogram matching and pseudo invariant feature spot methods are employed respectively to make relative radiometric correction for the TM/ETM+ image sequence.It reduces the uncertainty of the lighting and atmospheric conditions in the temporal phase for the images,and improves the accuracy of the NDVI and LSWI.Firstly,the cultivated land of study area is extracted quantificationally according to multi-temporal NDVI,then based on the physiological characteristics of rice in different growth stages,the paddy field containing partial dry land is extracted by defining the threshold of LSWI in June and August.On the basis,according to the feature that the NDVI of double cropping rice has two peaks which appear in June and September respectively and the NDVI of single cropping rice has one peak in August,the paddy field is extracted more accurately.The results show that the precisions of paddy field extraction are as high as 92.3%.Besides,the correlation coefficient between extraction area and statistical yearbook is 0.97 in county scale.

paddy field;LSWI;time-series NDVI;Jinhua

2014-11-09；

2015-01-16

国家自然科学基金项目“遥感影像小支撑规范化去卷积快速算法研究”(41371331)

刘元亮(1990-)，男，硕士研究生，主要从事遥感与GIS应用方面的研究。*通讯作者E-mai:liyan@nju.edu.cn

10.3969/j.issn.1672-0504.2015.03.007

TP751

A

1672-0504(2015)03-0032-06