[摘 要]
水稻是中国65%以上人口的主粮。及时准确地掌握水稻种植面积和空间分布,对于制定农业政策、保障粮食安全具有重要意义。本文提出了一种基于水稻种植模式的小型水体识别与阈值分割相结合的方法进行省域水稻田识别。该方法数据采用哨兵2号光学遥感影像,利用谷歌地球引擎(GEE)平台在广西壮族自治区进行了试验和精度验证。本文方法在水稻田识别方面表现出可靠、快速的性能,Kappa 系数为0.76,总体精度为92%,生产者精度为70%,用户精度为95%。这种计算开销节约的方法实现了提取精度与计算效率的平衡,可有效应用在省域或更大区域的单季稻水稻田识别中。
[关键词]
遥感指数;农业遥感;水稻田识别;哨兵2号;谷歌地球引擎(GEE)
引言
中国是世界上最大的水稻生产国,水稻产量占世界总产量的30.7%。水稻是65%以上中国人口的主粮。2006 年,中国水稻年总种植面积为2 940 万hm
2
,占世界种植面积的19%,居世界第二位
[1]
。及时准确地掌握水稻种植面积和空间分布情况,对于制定农业政策和保障粮食安全具有重要意义
[2]
。
由于新品种的开发以及氮肥和灌溉等的改进,中国水稻产量在过去50 年间增长了2 倍多。然而,中国水稻生产体系中的几个关键问题阻碍了水稻总产量的可持续增长,其中包括过度施肥和使用农药、水资源短缺等。特别是水稻田面积减少情况下,一些农民为了短期经济利益,将原来种植主粮的土地改种经济作物。这种土地利用管理方式的转变如果得不到及时监测,将影响中国的食品安全。然而,传统的水稻田种植面积信息的获取方法主要是人工统计,这种方法费时费力,准确率低。而遥感技术作为一种有效的工具,被用来分析和提供从局部到全球范围内不同尺度的数据,为监测作物生长提供精确、及时的植被物候和发育信息。
在遥感图像的分析中,使用光谱数据对地表特征进行分类是常见的方法
[3]
。然而,不同类型的农作物在某些生长期可能具有相似的植被光谱特征,导致分类错误。此外,水稻田内部的条件(如水位变化、病虫害等)也可以显著改变光谱反射特性,导致使用光谱数据进行准确分类的难度增加。因此,可以通过分析时间序列上的信息,如植物物候和种植模式来提高水稻田识别精度
[4]
。目前有多种应用于遥感时间序列图像的耕地制图方法。文献[5]使用陆地卫星(Landsat)时间序列复合影像半自动地探测农业生长季节的作物田块。文献[6]提出了一种综合使用高空间分辨率(20 m)和中分辨率(300 m)土地覆被信息进行年度耕地测绘的自动化方法。这类数据源主要是中低空间分辨率的图像,适用于种植面积较大、云量较少的地区。然而,对于地形复杂和多云天气较多的地区,稻田呈现出小尺度、零碎和分散的特点
[7]
,其识别需要空间分辨率更高的数据源。文献[8]于2017 年开发了一种方法,将密集的哨兵1号一级地距多视影像数据与精确的作物物候信息结合,以识别特定作物类型的物候时间序列模式。然而,作为雷达数据,哨兵1号不像光学影像那样能够直接识别地物的外观或形状,需要额外的算法和处理来进行地物分类和识别
[9]
。文献[10]提出了一种提取小水体的方法,采用指数假彩色合成和色彩空间变换,利用哨兵2 号影像精确快速地进行小型水体提取。然而,该算法并不针对水稻田提取。稻田在注水期具有明显的水体植被混合特征,故可在小水体研究基础上,结合稻田植被期呈现出的异于其他水体的特征,筛选出水稻田。综上,本研究采用哨兵2号光学影像数据,在保证高空间分辨率的同时实现了快速计算。在利用六角锥体(hue saturation value,HSV)色彩系统增强显示稻田与水体的基础上,利用多时相数据,结合单季稻的物候信息,判别本区域内水稻的种植耕作方式,利用该信息进一步精确识别水稻。
1 研究区和数据
1.1 研究区域
广西壮族自治区位于中国南部,地处东经104°28′—112°04′、北纬20°54′—26°24′。广西壮族自治区由14 个地级市组成,总占地面积23.76 万 km
2
。全区属亚热带至热带季风气候,年平均气温17.5℃~23.5℃。大部分地区雨热充沛,日照适中,干湿季分明,为水稻生长提供了有利条件。平均年降水量为1 653 mm,海拔在0~2 800 m之间。在平坦的平原和低山丘陵上,种植着多种农作物,其中水稻是主要的农作物类型。
1.2 卫星数据
哨兵2 号是欧洲航天局(European Space Agency,ESA)开发的一组遥感卫星,其目标是提供宽波段、高分辨率、多光谱遥感影像,以支持包括植被、土壤和水覆盖在内的地球环境监测。所有哨兵2 号数据在谷歌地球引擎(Google Earth Engine,GEE)平台上免费开放,时间分辨率为5 d
[11]
。本研究使用的哨兵2 号数据时间范围为2016-01-01 至2019-05-01。选取了空间分辨率为10 m的波段2、波段4和波段8作为提取水稻田的基础数据,其中心波长分别为490 nm,665 nm和842 nm。在阿尔克·吉斯(ArcGIS)10.2 中对研究区内14个地级市进行了边界裁剪,并将生成的矢量数据文件导入GEE平台。
2 方法
本研究方法的操作流程包括以下步骤。
1)第一步,潜在水稻田与水体识别。计算研究区域的归一化差异植被指数(normalized difference vegetation index,NDVI),将指定波段和NDVI 指数合成红、绿、蓝(red,green and blue,RGB)伪彩色图像。通过色彩转换系统增强显示,识别灌溉期可能成为潜在水田的水体。
2)第二步,精确水田识别。计算研究区域的地表水指数(land surface water index,LSWI)、增强植被指数(enhanced vegetation index,EVI)和归一化差异建筑指数(normalized difference built-up index,NDBI)。基于研究区域灌溉管理现状,根据水田水体的时空分布特征精确识别水田。
3)第三步,精度评估。利用目视判读法评估本水稻田识别方法的准确性。
总体技术路线图如图1 所示,图中L1C 表示哨兵2号影像数据的级别。
图1 技术路线
2.1 潜在水稻田识别
由于研究区云量较多,本研究在GEE 平台导入了哨兵2 号多光谱影像(Level-1C 级别,属于经过几何精校正的正射影像)的3 年数据(2017—2019 年)。在每一年的遥感影像数据中筛选出一季稻生长期的数据,组成图像集,并基于哨兵2号表示云掩码信息的QA60 波段,去除被云遮挡的像素。
单季稻的栽培周期包括播种期、灌水期、插秧期、抽穗期和收割期。插秧期通常从5 月中旬开始,到6月中旬结束。在插秧期和抽穗期之间,稻株会经历返青和分蘖。返青通常从6月下旬开始,到8月下旬结束。从抽穗到收割,稻粒会从乳粒发育到成熟粒,成熟粒出现在9 月。水稻田自4 月份灌水期起,稻田里就充满了水,并混杂着水稻植株,田地呈现出植被和水体混合的地物特征,直至10月份完熟收割。进一步从选取的原始遥感图像集中筛选出水稻田灌满水的时间段,并计算这些图像中每个像素的NDVI值,公式为
式中,
b
4
和
b
8
分别是遥感卫星图像中红光波段和近红外波段;
I
NDVI
是NDVI值。
在由B2、B8和NDVI构建的RGB 假彩色合成图像中,充满水并混有水稻植株的水田被增强显示为蓝紫色。水体显示为蓝色,森林显示为深绿色,而建筑物、裸地等地物显示为灰白色或土黄色。使用这种RGB 假彩色合成图像可以在地图上突出显示水稻田,从而将其与其他地物类型区分开来。
传统方法是使用RGB 作为图像的色彩系统,然而,在RGB 色彩空间中,光照条件的变化可能导致不同地物在同一颜色通道上的数值变化,给阈值选择带来困难
[12]
。HSV 色彩空间将颜色信息分为色调(hue,H)、饱和度(saturation,S)和亮度(value,V)三个分量。在遥感影像处理中,色调对应地物的颜色信息,可以用于区分不同的地物类型。饱和度反映了颜色的纯度和鲜艳程度,可以用于区分不同地物的材质和特征。亮度则表示图像的亮度信息,可以提供关于地物反射率和辐射度的信息。在HSV 色彩系统中,V 通道对光照的变化相对较不敏感,使得阈值选择更加稳定。将亮度信息与色彩信息分离,可以更好地适应光照条件的变化,提高识别的鲁棒性
[13]
。将遥感影像中的RGB 通道分别映射到不同的遥感指标,并转换到HSV 色彩系统进行识别,可以充分利用图像中包含的遥感信息,将遥感指数信息特别是水体信息转换为更稳定的遥感指标。
将由B2 波段、NDVI 和B8 波段合成的RGB假彩色图像转换为HSV 色彩系统后,基于多种阈值的效果比对,将H 值的范围设定为0.57~1.00,以涵盖各种水体类型。由于该方法基于变换后的HSV 色彩系统进行提取,故将该方法称为色彩变换法。筛选出的这部分混有稻田的水体将在下一步使用遥感指数阈值分割法进行区分,以提高识别的准确性。
2.2 基于时序信息精确提取水稻田
不同的土地利用类型在不同的时间段会表现出不同范围的遥感指标,在进行识别时,单一的遥感指标并不能反映全面的地物特征,也就无法实现地物的精确识别。
2.2.1 代表性指数选择
LSWI 以近红外波段(near infrared,NIR)和短波红外波段(short wave infrared,SWIR)的归一化比率计算,公式为
式中,
ρ
swir
和
ρ
nir
分别代表短波红外反射率和近红外 反射 率
[14]
;
I
LSWI
是地表水指数(land surface water index,LSWI)。
EVI是一种优化的植被指数,用于识别稻田、评估季节变化和跟踪作物生长的物候事件
[15]
。EVI的计算公式为
式中,
ρ
nir
、
ρ
red
和
ρ
blue
分别代表近红外、红色和蓝色反射率;
I
EVI
是增强植被指数(enhanced vegetation index,EVI),可以明显区分植被和非植被,但水稻、草地和耕地这三类植被,EVI 的差异很小,难以区分。
NDBI的计算公式为
式中,
I
NDBI
是归一化差异建筑指数(normalized difference built-up index,NDBI);
b
5
和
b
4
分别代表遥感卫星图像中红外波段和近红外波段
[16]
。在春季和冬季,水田和水体的NDBI 值差异明显,可以此区分水田、森林和水体。
2.2.2 单一作物的灌溉模式和阈值的确定
本研究方法基于水稻的独特灌溉模式:在作物生长之前,主田首先要在潮湿的条件下浸泡、翻耕和积水。在此期间,水田呈现出土壤和水体混合的状态
[17]
。然后,将育好的秧苗移栽到主田,并在主田灌入5~10 cm 深的水。随着水稻的生长,稻田呈现出植被和水体混合的特征
[18]
。本研究选取了六大地物类型的样本点,记录了样本点在一年内的LSWI、EVI、NDBI 指数的变化并作时序曲线图,如图2 所示。根据图2(a),从1 月到4月,耕地的LSWI值较低,而林地的LSWI值相对较高,高于0.2。通过将该时间段的阈值设置为-0.2~0.2,可以将林地和部分其他土地利用类型筛选掉。从6 月到10 月,通过将阈值设为0.12~0.52,可以筛选掉建筑用地、草地和部分其他土地利用类型。根据图2(b),从6 月到10 月,水稻田的EVI 值在0.2 以上,而水体的EVI 值则低于0.1。将阈值设置为0.1~0.5,可将水体筛选掉。根据图2(c),从1 月到4 月,耕地的NDBI 值相对较高,高于0,而水体和林地的NDBI值低于-0.1,其他作物的NDBI 值低于0。将阈值设置为0~0.2,可将水体、林地以及其他作物筛选掉。图上的红色框与文字标示出了所使用的阈值以及使用该阈值范围筛选掉的地物类型。根据上述选定的阈值区间对图像进行掩膜处理,将其他五种土地利用类型全部从地图中遮罩掉,得到水稻田提取结果。
图2 6种地物利用类型的时间序列图
2.3 精度评估
本研究提出的无监督分类方法具有快速的优点。其准确性通过目视判读来验证。目视判读操作简单,且能确保较高的准确性。在对高分辨率遥感图像进行解译时,视觉解译是基本且必不可少的。通过目视判读,可以验证计算机自动分类的结果。
为了验证本方法的分类准确性,在广西壮族自治区柳州市选取了一个面积为224 km
2
的小型研究区域,随机生成了1 000 个样本点。使用当地最新0.3 m 的高分辨率遥感影像作为视觉判读的基础。对于每个样本点,分类算法结果和目视判读结果分别记录在误差矩阵中(如表2 所示),以评估总体分类的准确性。由误差矩阵计算出生产者准确度、用户准确度、卡伯(Kappa)系数和总体准确度
[19]
。
3 结果
3.1 研究区域的水稻田识别结果
如图3 所示,在利用阈值分割法对水稻田进行精细化提取后,提取精度有了明显的提升。图3(a)显示了使用单一的基于色彩系统转换法的小水体识别结果。在小水体识别的基础上,图3(b)加入了LSWI 指数的阈值分割筛选,进一步将部分小水体像素筛除。图3(c)进一步加入了EVI 指数的阈值分割筛选,图3(d)为使用LSWI、EVI、NDBI 三个遥感指数阈值分割的提取结果。与仅使用LSWI 和EVI 指数相比,加入了NDBI 对精度提升的效果最为明显,不仅去除了非水稻田的水体,也去除了建筑物、林地、草地等地物类型。
图3 初筛与加入阈值分割法后的提取结果对比
广西壮族自治区的水稻田提取结果如图4所示。
图4 广西壮族自治区各市水稻地图
(审图号:GS(2019)1822号)
广西壮族自治区14 个研究城市的稻田提取结果如表1 所示,表中列出了广西壮族自治区每个城市提取出的水稻田面积及其在城市总面积中的占比。
表1 广西壮族自治区各市水稻面积与城市总面积
续表1
3.2 精度评定
基于目视判读法,首先,对仅使用色彩变换法初筛出的水稻田和小水体混合地块的结果进行了精度验证,结果记录在表2 中。阈值分割法从初筛出的水体和潜在稻田中,挑选出了37%的像素判别为水稻田。经阈值分割法进一步精细化提取的水稻田识别结果的精度验证信息记录在表3中。
表2 初筛提取结果的误差矩阵 单位:个
表3 水稻田识别结果的误差矩阵 单位:个
所有分类结果和每个样本点的实际类型都记录在误差矩阵中,如表3所示。根据矩阵,计算出Kappa 系数为0.764,总体精度为92.2%,用户精度为95.6%,生产者精度为70.7%。其中,Kappa 系数较低,可能是水稻田与非水稻田样本分布较为不均衡造成的。
本文方法与其他方法相比,在单季稻田识别精度上有所提高,达到了92.2%的总体精度。文献[20]采用哨兵1 号和哨兵2 号数据对江汉平原水稻种植区进行了提取,达到了79%以上的总体提取精度。文献[21]基于哨兵2号数据分别采用最大似然法、支持向量机法和面向对象分类法对研究区水稻进行分类识别,其中最大似然法总体分类精度为89.35%,支持向量机法总体分类精度为84.75%,面向对象分类法总体分类精度为76.9%。
4 结束语
本研究的整个过程都基于GEE 平台,与传统的遥感图像处理软件相比,基于云的计算平台表现出了较好的计算性能,实现了提取精度和计算效率的平衡,更加适用于省级大面积水稻田识别。本文所提出的基于多种遥感指数的水稻识别研究方法为监测其他作物的种植面积提供了重要参考。该方法考虑了作物物候的详细信息,为作物测绘提供了一种有前景的方法。虽然使用小水体识别方法无法提高识别旱地作物的准确性,但仍可通过对多种遥感指数设置不同的阈值来区分目标作物和其他土地类型。需要说明的是,由于不同地区的气候不同,水稻种植模式也会不同。在其他地区进行水稻识别时,需要重新选择种植模式、设置适当时期的阈值来完成分类。
参考文献
[1]FAIRHURST T,DOBERMANN A.Rice in the global food supply[J].World,2002,5(7):454,349-511,675.
[2]何泽,李世华.水稻雷达遥感监测研究进展[J].遥感学报,2023,27(10):2363-2382.
[3]柴颖,阮仁宗,柴国武,等.基于光谱特征的湿地植物种类识别[J].国土资源遥感,2016,28(3):86-90.
[4]柳文杰,曾永年,张猛.融合时间序列环境卫星数据与物候特征的水稻种植区提取[J].遥感学报,2018,22(3):381-391.
[5]GRAESSER J,RAMANKUTTY N.Detection of cropland field parcels from Landsat imagery[J].Remote Sensing of Environment,2017,201:165-180.
[6]MATTON N,CANTO G,WALDNER F,et al.An Automated Method for annual cropland mapping along the season for various globally-distributed agrosystems using high spatial and temporal resolution time series[J].Remote Sensing,2015,7(10):13208-13232.
[7]蔡耀通,刘书彤,林辉,等.基于多源遥感数据的CNN 水稻提取研究[J].国土资源遥感,2020,32(4):97-104.
[8]BARGIEL D.A new method for crop classification combining time series of radar images and crop phenology information[J].Remote Sensing of Environment,2017,198:369-383.
[9]查东平,蔡海生,张学玲,等.基于多时相Sentinel-1 水稻种植范围提取[J].自然资源遥感,2022,34(3):184-195.
[10]BIE W,FEI T,LIU X,et al.Small water bodies mapped from Sentinel-2 MSI(MultiSpectral Imager)imagery with higher accuracy[J].International Journal of Remote Sensing,2020,41(20):7912-7930.
[11]SON N T,CHEN C F,CHEN C R,et al.Classification of multitemporal Sentinel-2 data or field-level monitoring of rice cropping practices in Taiwan[J].Advances in Space Research,2020,65(8):1910-1921.
[12]杨凯,刘如飞,崔立军,等.全景影像中道路交通标志的自动识别[J].遥感信息,2020,35(4):56-61.
[13]王蜜蜂,缪剑,李星全,等.基于RGB 和HSI 色彩空间的遥感影像阴影补偿算法[J].地理空间信息,2014,12(6):107-109+4.