基于哨兵SAR数据和多光谱数据的水稻识别研究

张征云,江文渊,张彦敏,罗 航 (.天津市生态环境科学研究院,天津 3009;.天津环科环境规划科技发展有限公司,天津 3009)

水稻是国家实施粮食安全战略的基础,我国粮食种植结构中,水稻以25.76%的占比居第2。掌握水稻的种植面积动态信息,可为监测水稻的生产状况、评估水稻产量、规划粮食生产布局、合理调整农业结构等提供依据[1]。同时,水稻的种植信息对于水资源的合理利用和监测以及评估人类活动对于大气环境效应的影响均有重要的意义。对于水稻提取的研究以使用光学遥感数据为主,由于部分地区在农作物生长时期内常被云雾覆盖,实时获取农作物种植与生长情况困难,而合成孔径雷达(synthetic aperture radar,SAR)与光学遥感相比,不依赖阳光、能够穿透云雾甚至土壤表层,实现全天时全天候观测地面,提供不同于光学遥感的信息,在水稻研究领域得到大量应用[2]。近年来国内利用遥感技术获取水稻等作物面积的相关研究较多,光学数据和SAR数据已成为作物种植信息获取的重要途径且方法趋于成熟[3]。

光学影像在近红外与短波红外波段对农作物有较好的区分。李杨等[4]和张薇薇[5]基于多特征波段构建决策树分类提取水稻信息。由于光谱分辨率的限制易产生“同物异谱”和“异物同谱”现象,为提高作物识别的精度,时序光谱和指数影像被广泛使用[6]。已有学者采用不同光学影像数据源,基于植被指数特征、波段组合特征等分类特征对水稻信息进行识别提取[7-12]。水稻的后向散射系数随时相的变化比其他作物明显,这是利用SAR进行水稻识别提取的主要机理[13]。SAR通过主动发射微波脉冲,并接收目标物体对微波信号的后向散射信号,根据其强度差异进行成像。SAR遥感技术在水稻信息识别中的应用始于20世纪八十年代。LE TOAN等[14]通过分析X波段SAR影像中地物的后向散射特征提取水稻信息。随着多极化、多入射角、更高分辨率以及全极化SAR的发展,已有学者基于SAR的交叉极化数据,融合时序统计参数与时序曲线相似性特征,对国内外不同水稻种植区进行了识别提取[15-18]。雷达遥感影像与光学遥感影像信息的结合能够更好地实现对水稻的识别和监测,也是目前的研究热点[19]。栗云峰等[20]提出了一种融合多时相的GF-1号卫星多光谱影像和Sentinel-1雷达影像的水稻田提取方法。

综上所述,对水稻的识别与提取无论是基于SAR数据还是基于光学数据的研究都比较多,但对基于这2种数据源的水稻识别精度进行对比的研究尚不多见,同时对华北单季稻稻作区的提取研究也较少。笔者以天津市为研究区域,基于Sentinel-1 SAR时序数据和Sentinel-2多光谱数据对华北单季稻稻作区天津市水稻种植区提取方法进行研究,对比分析了基于SAR数据和多光谱数据在天津水稻分类提取上的效率和精度,提出最利于华北单季稻稻作区水稻识别的方法,并对天津市水稻种植区2016、2018、2021年的变化趋势进行分析。

1 研究区概况

天津处于海河流域下游,上游各河道水流入渤海,其中北系有蓟运河、潮白新河、北运河和永定新河水系,西系有大清河和子牙河水系,南系有独流减河和南四河水系,境内河道沟渠纵横交错,洼淀较多,坑塘密布。天津历史上种稻较早,水稻栽培与育种研究历史久远,在北方稻区有较大的知名度。在全国综合农业区划中,天津属华北单季稻稻作区华北北部平原中早熟亚区[21]。根据天津市第三次农业普查数据,天津市水稻种植集中分布在宝坻和宁河区,除中心城区外在其他区均有零星种植[22]。

每年4月上旬至10月中旬是天津水稻生育期[23],生育期170～190 d。4月中旬播种、育秧(最早在3月末开始播种),5月上、中旬插秧(移栽),插秧前10～15 d灌水泡田,耙地找平。5月下旬至6月中上旬为水稻分蘖期,6月下旬至7月上旬为拔节期,7月中旬前后进入孕穗、抽穗期(各约15 d),8月底前进入齐穗期,10月进入成熟期,10月底前收获。

2 数据源与研究方法

2.1 数据源

SAR数据源为覆盖天津市水稻产区的Sentinel-1 GRDH数据,坐标采用UTM/WGS84投影系统。下载时结合天津水稻的生长时间,不同年份的4—11月每月下载2～3幅影像,涵盖了研究区水稻生长发育及收割后的整个时间过程,以便高精度地提取不同生长阶段的水稻信息。数据获取模式为干涉宽幅模式(IW),产品类型为Level-1级地距影像(GRD)。时间分辨率12 d,空间分辨率5 m×20 m,幅宽250 km。

光学影像数据源为覆盖天津市的Sentinel-2数据,覆盖全天津市需5景。Sentinel-2是全球环境与安全监测的第2颗卫星,卫星携带多光谱成像仪(MSI),涵盖了从可见光到短波红外共13个波段,按照空间分辨率可将波段分为10、20和60 m,不同分辨率的波段可应用于陆表监测、植被监测、地质监测、海洋监测等不同领域。Sentinel-2卫星幅宽达290 km,双星重复访问,时间分辨率为5 d。研究区2016和2018年的Sentinel-2影像数据为经过正射校正和亚像元级几何精校正后的L1C大气表观反射率产品,未进行大气校正。2021年光学影像数据为经大气校正后的L2A级产品。

研究所用的影像全部从欧洲航天局哥白尼开放数据访问中心(https:∥scihub.copernicus.eu/dhus/#/home)下载,数据源信息详见表1。

表1 研究区数据源Table 1 Study area data source

2.2 数据预处理

Sentinel-1和Sentinel-2数据预处理均采用欧洲航天局为处理Sentinel系列卫星数据而开发的配套软件SNAP(Sentinels Application Platform)。

Sentinel-1 IW GRD级数据是由SLC级经过多视处理、采用WGS84椭球投影至地距的聚焦数据,共含4个波段,Amplitude_VH、Amplitude_VV波段表示VH/VV通道的振幅,Intensity_VH、Intensity_VV波段表示VH/VV通道的强度(振幅的平方)。数据依次经热噪声去除、辐射定标、相干斑滤波(Refined Lee滤波器)、地形校正等处理后得到线性比例单位的后向散射系数σ0。以VV通道为例,σ0值在0～1之间,平均后向散点值分别为0.07(VV)和0.02(VH),数据视觉对比度较低。为便于可视化及数据分析,将雷达后向散射系数进行分贝化处理,公式为

(1)

式(1)中,P、P0分别为目标量和参考量。对于后向散射系数σ0的分贝化,按式(2)执行对数变换。变换后代表后向散射系数强度的分贝值大约在35～10 dB之间。

σ0,dB=10×lgσ0。

(2)

对Sentinel-2 L1C级多光谱数据进行大气校正后得到L2A级数据产品。将20 m分辨率处的6个波段使用最近邻法重采样到10 m分辨率,结合4个原始10 m分辨率的波段共得到10个10 m分辨率的波段,对5景影像进行裁剪后拼接形成天津市全域10 m分辨率的多光谱影像。

2.3 分类方法

2.3.1阈值法

在多数情景下,SAR图像中的目标区域与背景区域之间的灰度值存在较大差异[25]。阈值法分类的原理是利用不同通道对各类地物的敏感度差异进行单波段或多波段阈值的选定,从而将目标地物提取出来,调用一个或多个已提取地物的信息源,分层实现地物的再提取,最终完成影像分类。

2.3.2监督分类法

监督分类又称训练分类,是指用被确认类别的样本像元去识别其他未知类别像元的过程。通过Google高精度影像目视判读和野外调查,形成水稻、旱地、水体、林地、建设用地等用地类型的先验知识,对每类用地分别选取不低于100个训练样本,采用监督分类中最常用的最大似然分类法和支持向量机法进行分类提取。

2.4 分类精度验证

2.4.1样本与验证数据

样本质量直接关系到分类信息提取的精度。研究样本数据来自现场实地踏勘和依据高分辨率影像进行目视解译识别2种方式,选择具有典型性、代表性的纯净像元作为样本,选取水稻、旱地、水体、林地、建设用地各100处作为验证数据,用于分类结果的精度验证。

2.4.2精度验证

采用混淆矩阵(confusion matri)定量评价同一位置的实际像素地类和分类结果地类是否一致。判断指标有Kappa系数、总精度、生产者精度(某类)和使用者精度(某类),数值越接近1,分类结果越准确。其中Kappa系数和总精度是对分类结果的总体评价,生产者精度和使用者精度则反映某特定类别的分类精度。生产者精度指通过分类法将整个影像的像元正确分为某类的像元数与该类真实参考总数的比值,该值用于比较分类方法的好坏。用户精度是指正确分到某类的像元总数与通过分类法将整个影像的像元分为该类的像元总数比值,该值用来表示分类结果中各类别的可信度和可靠性。

3 结果与分析

3.1 后向散射特征分析

图1为不同地物在VV、VH极化下的后向散射系数时序变化特征。由图1可知,不同地物VH、VV极化的时序后向散射系数变化趋势基本一致。由于水平返回到传感器的垂直传输波的比例较小,VH波段的强度大多低于VV波段。在VV、VH极化通道下后向散射强度均表现为水体<水稻<旱地<林地<建筑物。

图1 不同地物在VV、VH极化下的后向散射系数时序特征Fig.1 Time series characteristic curves of backscattering coefficient for different surface features under VV and VH polarization

水稻VH、VV极化的时序后向散射系数分别在-28.25～-14.09、-20.32～-6.63 dB之间,波动幅度较大,波动趋势基本一致;建筑物VH、VV极化的时序后向散射系数在-6.38～-2.27、-0.24～3.05 dB之间波动,波动幅度较小,两者波动趋势类似,较为平稳;水体VH、VV极化的时序后向散射系数分别在-28.1～-24.96、-25.95～-16.25 dB之间波动,VV极化的时序后向散射系数波动趋势高于VH极化;林地VH、VV极化的时序后向散射系数分别在-16.87～-12.29、-10.80～-5.48 dB之间,波动趋势基本一致;旱地VH、VV极化的时序后向散射系数分别在-22.79～14.53、-15.04～-7.88 dB之间波动,波动趋势基本一致。综上可以看出,相对其他地物而言,水稻的后向散射系数年内变化幅度最大(波幅分别达13.9和14.0 dB),且有较为明显的2个波峰和波谷对应不同的生长期。

水稻不同生育期的时序后向散射系数为识别水稻提供了必要的信息。VH通道下,水稻与其他地物的时序后向散射曲线明显不同。在移栽期及生长初期,水稻植株个体小,稻田的后向散射接近于水面,散射强度较低。随着水稻的生长,后向散射以稻体散射和冠层面散射为主,后向散射系数不断增大,在拔节期后期、孕穗期前期达到高值,随后在抽穗期开始下降,在齐穗期降至低值,成熟期田中水全部排干,导致后向散射强度较高。利用这一特点,选取不同生育期的SAR影像进行合成,水稻与其他地物之间有明显区别而易于提取。

针对水稻不同于其他地物的后向散射系数时序变化特征,结合水稻关键生育期,选取3个特征时间,VV极化选择5月30日、6月23日、8月10日这3景影像进行RGB假彩色合成,VH极化选择5月30日、7月5日、8月10日这3景影像进行RGB假彩色合成(图2)。从视觉上突出水稻与其他地物的区别,进而使用监督分类方法来提取水稻面积。由图2可以看出,VH极化下的合成效果更加均一稳定,利于水稻的识别与提取。

绿色为水稻。图2 SAR数据波段的假彩色合成影像Fig.2 False color synthesis of SAR data in different phases

3.2 光谱特征分析

对2021年水稻生长各阶段及收割后Sentinel-2不同波段组合的对比分析发现,在水稻移栽期(5月30日)和成熟期(10月21日),稻田和周边其他作物光谱特征有明显区别,易于提取。以2021年影像为例,图3反映了在移栽期和成熟期近红外(B8)、短波红外(B11、B12)和可见光红光(B4)组合下,水稻光谱特征凸显。在B11+B8+B4、B12+B8+B4和B11+B8+B4波段组合下,水稻与其他地物的分离度较高,均在1.98～2.00之间。

B2、B3、B4—可见光;B8—近红外;B11、B12—短波红外。图3 Sentinel-2多光谱影像不同波段合成影像Fig.3 Sentinel-2 multispectral image synthesis of different bands

3.3 分类提取与精度验证

3.3.1分类提取

使用SAR和光学影像数据源所需波段及合成的影像,选取典型区(潮白新河以北宝坻区),采用阈值法、最大似然法、支持向量机法进行水稻分类提取(图4)。

图4 基于不同数据源和不同分类方法的分类结果Fig.4 Classification results based on different data sources and different classification methods

阈值法基于SAR数据源提取水稻信息。根据不同波段下二值法中色带量级确定水稻的像元散射系数范围。提取规则设置为同时满足VH通道下后向散射系数5月30日在-29.0～-21.0 dB区间,7月5日在-19.0～-13.5 dB区间,8月10日在-25.3～-16.2 dB之间,在ArcGIS软件的栅格计算器中运用con函数获取满足条件的水稻种植范围。

采用最大似然法和支持向量机法对2021年基于SA300R时序数据的影像、移栽期(5月30日)Sentinel-2影像、成熟期(10月21日)Sentinel-2影像进行分类。该研究所选取的合成影像均为水稻与其他地物的不同波段组合,水稻与其他地物区分度高(分离度均在1.98～2.00之间),易于识别提取。

3.3.2精度验证评价

使用混淆矩阵对分类精度进行评估。阈值法只针对水稻信息的提取,因此仅从生产者精度和用户精度进行分析评估。基于时序SAR数据、移栽期光学影像(B12+B8+B4)和成熟期光学影像(B12+B8+B4)的最大似然法和支持向量机法分类精度,则采用总分类精度、Kappa系数、各类地物生产者精度和用户精度4个指标进行分析评估(表2)。由于水稻明显有别于其他地物的时序SAR后向散射特征和光谱特征,无论是基于SAR数据还是光学遥感数据,3种分类法对水稻信息提取的生产者精度和用户精度均在90%以上。最大似然法分类结果的总体精度和Kappa系数优于支持向量机法,在水稻信息提取精度上,支持向量机法优于最大似然法。采用支持向量机法对水稻成熟期B12+B8+B4波段组合光学遥感数据进行分类提取效果最佳,水稻生产者精度为99.61%,用户精度为99.40%。

表2 分类精度验证Table 2 Classification accuracy verification

3.4 研究结果

根据以上分析,基于水稻成熟期的Sentinel-2多光谱光学影像B12+B8+B4波段组合影像,采用支持向量机法提取水稻种植区范围,并结合移栽期光学影像和SAR数据对水稻信息提取结果进行目视修正,水稻种植范围识别度接近100%。结果显示,天津水稻分布范围最北至青甸洼,最南至独流减河,独流减河以南仅在静海团泊洼东南有少量分布。天津水稻种植集中分布在潮白新河和蓟运河间的宝坻和宁河两区,两区在2016、2018和2021年的水稻种植面积分别占天津水稻种植面积的90.53%、86.32%和83.76%(表3、图5)。天津水稻种植面积从2016年的399.04 km2增长至2021年的764.55 km2,增加91.60%,符合天津市实际情况。水稻种植面积大幅度增加,得益于水资源环境持续改善,更重要的是天津市小站稻产业振兴政策的支持,地方种植水稻的积极性很高。

表3 2016—2021年天津市分区水稻种植情况Table 3 Rice planting in Tianjin from 2016 to 2021

4 讨论与结论

4.1 讨论

目前国内采用遥感手段识别水稻种植范围的研究较多,但对华北单季稻稻作区信息的提取以及对比基于哨兵SAR数据和光学数据的水稻识别精度的研究尚少。该研究以天津为研究范围,基于Sentinel SAR数据和光学数据分别提出一种水稻信息提取方法。对比分析后提出,基于水稻成熟期的Sentinel-2 B12+B8+B4波段组合,采用支持向量机法提取水稻信息是一种适用于华北单季稻识别的方法。

该方法分类总体精度为90.18%,Kappa系数为0.86,水稻生产者精度为99.61%,用户精度为99.40%,水稻提取精度高。现有基于多光谱数据提取水稻信息的方法多综合利用各类光谱特征、纹理特征、指数特征等,笔者提出的方法主要是基于水稻特定时期特定光谱波段组合特征来提取,提取方法所需数据量少且提取精度较好。

基于Sentinel-2光学遥感影像可以快速识别水稻,但是由于水稻种植时间不同会有所差异,仅依据基于生长期内某个时间点拍摄的遥感影像进行水稻信息提取会存在一些漏提现象。对遥感影像分析发现,部分水稻田种植时间较晚,在5月底尚未灌水移栽,在影像上表现为旱地或裸地,但在10月的收获季影像及实地调查中显示为水稻田,另有稻田早于区域成熟期进行了收割。这种情况涉及水稻面积小而分散,可将移栽期和成熟期两期影像结合起来,适度进行查漏补缺的人工干预,可使结果更加精确。

Sentinel-2光学影像数据具有来源免费、时间和空间分辨率高(5 d、10 m)、处理简易的特点,可以作为华北单季稻作区水稻种植面积监测的优选数据源。在北方地区,水稻移栽、分蘖期(4—5月)和成熟期(10月)获取低云量或无云量多光谱数据概率较高,基本可满足水稻种植区识别需求。SAR在水稻识别方面不受天气条件影响,在光学影像因云量太高而无法使用的情况下,可作为补充数据,使水稻的监测手段更加完善。笔者提出的方法在技术上相对简单易行、提取精度高,可为提高我国北方单季稻作区水稻监测精度与效率提供参考。

4.2 结论

(1)该研究提出一种基于Sentinel SAR的华北单季稻作区的水稻识别方法。时序后向散射系数变化特征是水稻识别与提取的关键所在。与其他地物不同,水稻在移栽期、拔节期、抽穗期的时序后向散射系数曲线形成一个显著波峰,组合移栽期、拔节期、抽穗期3个时期的SAR波段可识别水稻。采用阈值法、最大似然法和支持向量机法对水稻信息进行提取,其中支持向量机法较优,水稻生产者精度和用户精度分别为97.95%和92.92%。

(2)该研究还提出了一种基于Sentinel光学数据的华北单季稻作区水稻识别方法。在水稻移栽期和成熟期,Sentinel-2近红外(B8)、短波红外(B11、B12)和可见光红光(B4)等波段组合易识别水稻。分析结果显示,成熟期B12+B8+B4波段组合影像分类效果更佳,在总体分类精度上最大似然法优于支持向量机法,在水稻提取精度上支持向量机法优于最大似然法,水稻生产者精度和用户精度分别为99.61%和99.40%。最大似然法分类总体精度达93.93%,Kappa系数为0.92,分类效果好,推荐可用于土地覆盖类型解译。

(3)基于水稻成熟期的Sentinel-2 B12+B8+B4波段组合,采用支持向量机法进行分类提取水稻是一种适用于研究区水稻识别的推荐方法,基于SAR数据的方法可作为补充方案。利用推荐方法计算出研究区2016、2018和2021年水稻种植面积分别为399.04、586.67和764.55 km2,5 a增加365.51 km2,符合天津市实际情况。该方法在技术上简单易行,可为提高我国北方稻作区水稻监测精度与效率提供参考。