基于多时相Sentinel-1水稻种植范围提取

查东平， 蔡海生， 张学玲， 何庆港

(1.江西农业大学/江西省鄱阳湖流域农业资源与生态重点实验室，南昌 330045； 2.江西旅游商贸职业学院，南昌 330100)

0 引言

及时准确地监测和获取水稻种植的空间分布范围对保障全国粮食安全具有重要意义[1-2]。传统的统计方法有抽样调查和逐级上报，这两种方法获取数据成本巨大，且存在很严重的信息滞后[3]。遥感技术为大尺度农作物种植范围监测和提取提供了有效的技术支撑[4-5]。目前多采用Landsat系列、MODIS或Sentinel-2等光学卫星数据对农作物信息进行提取，这些数据免费且易于获取，相关学者基于此类数据开展了大量研究并取得了较好的成效[6]。Zhu等[7]利用Landsat系列数据对我国南方水稻产区的水稻种植面积进行了监测，结果认为2000—2019年水稻总种植面积减少了171 474 km2； Cao等[8]采用Landsat数据识别早稻、晚稻、双季稻和中稻，总体精度在71.83%～85.81%之间； Li等[9]通过融合MODIS和Landsat数据，采用随机森林模型和决策树对农作物进行提取，精度分别达到90%和89.7%。利用光学卫星影像数据对水稻田的提取方法多要依靠多时相数据，单时相数据提取的精度往往比较低，但是南方地区水稻生长季节往往为多云雨天气，光学数据质量难以保障，影响了实际应用成效[10-11]。

合成孔径雷达(synthetic aperture Radar,SAR)不受天气影响，可以全天时全天候对地观测，有效解决了光学数据获取难的问题[12-14]。基于SAR数据的水稻田分类虽然还没有大范围业务化应用，但是也开展了大量的研究，如Son等[15]利用Sentinel-1 SAR数据进行了水稻田提取，第一季和第二季水稻种植面积提取精度分别为84.2%和82.6%，用户精度分别为82.1%和85.3%； Graman等[16]利用Sentinel-1数据对Rabi地区水稻田进行提取，总体精度达到91.5%，Kappa系数为0.83。前人的研究证明了Sentinel-1数据VV(垂直)和VH(交叉)极化数据对提取水稻田的可行性。但是研究没有从水稻不同物候期数据去进行对比，筛选水稻田识别的最佳时相数据。目前成熟的，研究较多的分类方法包括最大似然法、支持向量机等，基于机器学习算法的分类方法应用在SAR图像水稻田提取的研究不多[17-18]。本研究拟通过分析不同物候期SAR数据水稻田与其他地类之间的可分离性，筛选出水稻田识别最佳时相组合，采用随机森林方法进行分类，以期得到较高精度的分类结果，为今后在水稻田识别和水稻生长监测中，因为云雨天气难以获得连续高质量的光学数据的情况下，采用SAR进行水稻监测提供方法上的参考。

1 研究区概况与数据源

1.1 研究区概况

研究区(图1)位于江西省南昌市南昌县北部、赣江下游蒋巷镇。地处E115°65′，N28°32′。蒋巷镇四面环水，东临鄱阳湖、南北傍赣江、西接南昌市城郊，赣江南、中支流环抱而过。全镇国土总面积约为240 km2，南北宽约8 km，东西长约45 km。受亚热带季风气候影响，研究区内气候温暖湿润，光照充足，主要种植作物为水稻，是南昌县重要的水稻种植区。

图1 研究区位置Fig.1 The location of research area

1.2 SAR数据获取与处理

选取涵盖2020年早稻整个物候期的8期Sentinel-1地距影像产品数据(https: //search.asf.alaska.edu/)。SAR过境时间分别对应水稻生长的出苗期(4月9日)、三叶期(4月21日)、移栽期(5月3日，直播稻为幼苗早期)、返青期(5月15日，直播稻为幼苗中期)、分蘖期(6月2日)、孕穗期(6月14日)、乳熟期(7月8日)、收割和二季晚稻移栽期(7月26日)。在ENVI SARscape 5.5软件中对SAR数据进行处理(图2)，包括读取、多视处理、配准和辐射定标等操作，采用“Adaptive Non-Local SAR Fil-tering”方法进行滤波处理，减少相干斑噪声。

图2 多时相Sentinel-1数据处理方法Fig.2 Multi-temporal Sentinel-1 data processing method

1.3 外业数据采集

为获取研究区的解译样本及对结果进行验证分析，对研究区开展了现场踏勘调查，获得了研究区主要地物类型分布，此外于2020年7月17日采用了大疆DJI Mavic 2无人机对研究区进行了飞行录像(无人机飞行轨迹如图3所示)。为节约成本，无人机采用“抬头”拍摄，飞行过程中调整拍摄角度，从受局限的俯视过程过渡到开阔的视角，“抬头”拍摄镜头前方可观测距离超过1 km，视野范围影像截图如图4所示。

图3 无人机飞行轨迹图Fig.3 Flight trajectory map of UAV

图4 无人机录像截图Fig.4 Video screenshot of UAV

本文研究目的是提取水稻种植范围，结合土地利用/覆被类型特点，参考邓刚等[19]、张悦琦等[20]对水稻种植范围提取辅助分类，依据多次野外调查结果，确定分类系统包括： 水稻田、林地、水体、莲田、人工表面和湿地等6种地类。根据外业调查和无人机飞行录像数据，在同时相GF-1影像图上勾绘地物标志样本。每类地物人工勾绘采选200个样本点，选择手工勾绘确保在影像图上获得纯净像元，便于后期分析。

在ArcGIS中对样本点矢量数据进行编号，奇数序号样本点导出并命名为训练样本，偶数样本点导出并命名为验证样本，各地类训练样本点和验证样本均100个，并且大致均匀分布于研究区。训练样本点和验证样本点分布情况如图5所示。

图5 样本点分布图Fig.5 Sample points distribution map

2 研究方法

2.1 后向散射的一般特征

SAR使用微波作为发射电磁波，电磁波照射到地物上，一部分能量通过反射进入雷达接收天线，另外一部分能量则透过地物表面继续向下[21]。雷达照射表面如果是粗糙地物发生散射，能量朝各个方向，如果是光滑表面则发生镜面反射[22]。雷达接收天线接收的能量取决于雷达波长、入射角和地面粗糙程度[21]。一般来说长波比短波具有更好的穿透性，对于探测冠层下垫面具有更好的效果。入射角越小后向散射回的能量越大，大角度的回波信息中包含了地物粗糙程度信息[21,23]。

通常，水面和光滑的地表可认为光滑表面，电磁波照射后易发生镜面反射，雷达天线接收不到或者只能接收少部分回波信号，在SAR图像上呈现黑色[24]。随着入射角的变化，接收到的信号也会有所不同，当入射角达到90°时，反射信号直接进入雷达天线，图像上呈现明亮点[24]。不同极化数据，水平极化和垂直极化反射和散射后到接收器上的能量不同，两种方式的差别可以达到5 dB，同向极化(HH和VV)与交叉极化(HV和VH)的后向散射系数也不同，二者可相差10～15 dB[25]。除水体，大部分地物都可认为是程度不相同的粗糙表面[16,25]。根据瑞利散射原则，地物散射特性跟地物表面和波长有关，地物的散射十分复杂。通常，粗糙表面的散射是一条平缓的曲线，曲线斜率变化与入射角相关性较小。农作物在不同的物候期，表现出来的散射特性不同。SAR接收天线接收到的地物后向散射能量受地物特性影响，如植被的密度、植被的含水量、植被高度、下垫面类型等。不同时相、作物处于不同的物候期，会产生不同的雷达回波，是雷达监测及农作物制图的基础。

2.2 随机森林分类法

本文随机森林训练流程如图6所示。随机森林是一个包含多个决策树的分类器[26]。1995年贝尔实验室华裔美国人何天琴(Tin Kam Ho)提出了随机决策森林(random decision forest)，2001年Leo Breiman和Adele Cutler发展并推出随机森林算法[27-28]。该算法通过构建多个决策树作为基本的分类器，每个决策树都是一个分类器，本质属于机器学习[26]。随机森林方法通过随机方式从原始训练集中选择K个训练集，对训练集构建分类决策树，集成多个分类投票的结果，将投票次数最多的分类作为最终地类并输出[26]。

图6 随机森林训练流程图Fig.6 Random forest training process

本研究随机森林分类通过对开源软件ENMAP-BOX v2.1.1源码封装后在ENVI 5.5软件中运行。

2.3 地物散射特征分析

2.3.1 水稻田散射特征分析

4月9日—7月26日水稻田的后向散射系数均值见表1，VV极化和VH极化后向散射系数变化趋势见图7。

表1 水稻田不同时相平均后向散射系数Tab.1 Average backscattering coefficient of paddy field in different time-phase (dB)

图7 水稻田不同时相平均后向散射系数Fig.7 Average backscattering coefficient ofpaddy field in different time-phase

水稻属于水生植物，大部分时间水稻根部被水体淹没，雷达波照射到水稻后大概可以分为3种散射过程：

1)水稻下垫面水体的镜面反射，水稻生长的早期，镜面反射是主要的散射过程，因此早期阶段水稻田的后向散射值较低。

2)水稻叶面的奇次散射，电磁波照射到水稻叶片后一部分电磁波发生奇次散射，直接进入雷达接收天线。

3)水稻叶片和茎秆的体散射，主要表现在雷达波到达水面后，水面反射到水稻茎秆或者叶片发生的多次散射，雷达波到达叶片、叶片散射到其他茎秆叶片或者水面发生的散射。

研究区水稻为直播稻，但水稻种植时间略有不同，4月初部分地区水稻已播种，部分地区水田尚未耕种长满野草，野草叶面积指数要比水稻种植区高，其后向散射系数也更高。4月21日，水稻处于三叶期，但水稻叶面积指数仍小于野草，后向散射系数要低于4月9日。随着水稻的生长，四月下旬三叶期水稻冠层本身散射和冠层散射逐渐增强。从极化情况来看，VV极化的穿透性要高于VH，其后向散射能量包含更多的植被冠层和植被下垫面的回波能量。4月下旬至5月中旬，随着水稻的生长VV极化数据呈现上升趋势，但是VH极化数据变化情况不大。根据地面数据采集情况，拔节、孕穗期之后水稻含水量逐渐下降，水稻的介电常数与含水量存在正相关性，结果受水稻的介电常数和降雨等情况的综合影响。7月26日已种植第二季晚稻，种植方式为插播，正处于返青期，此时仍受水面镜面反射主导。

2.3.2 其他地物散射特征分析

对莲田、林地、人工表面、水体和湿地等样本点在SAR时相序列后向散射数据中进行取值计算，各地物VV极化后向散射系数均值结果如表2和图8所示，地物散射系数箱型图见图9； 各地物VH极化后向散射系数均值如表3和图10所示，箱型图见图11。从VV极化各地物后向散射系数均值来看，各地物后向散射系数在-21～-5 dB之间，其中水体最低，人工表面最高。从图9可知，仅从VV极化数据，水体、莲田两种地类均具有较好的分离度，最佳区分的时相分别为5月15日、6月2日和6月14日。莲田与其他地类进行区分，最佳区分时相为4月9日和5月3日。值得一提的是，随着时间推移，莲田的VV后向散射系数逐渐提升，在5月3日之前，莲田的后向散射系数要低于水稻田、湿地等地类，6月14日之后，莲田后向散射系数要高于这两种地类，在6月26日甚至要高于林地的后向散射系数。林地的后向散射系数一直比较稳定，而且标准差较小，各时相均能与水稻田和莲田区分，与人工表面的最佳区分时相为4月21日，与莲田在4月9日、4月21日、5月3日和7月26日数据中均能较好区分，与水稻田在4月21日、6月14日、7月8日和7月26日均能较好区分； 与湿地在7月8日和7月26日具有较好的分离度。水稻田在水稻生长后期与其他地类能够较好区分，如6月14日、7月8日和7月26日。7月8日和7月26日湿地与人工表面、林地、莲田和水稻田具有较好的分离度，但是这2期数据湿地VV极化后向散射系数与水体极为相似，但是其他几个时相湿地可以与水体进行分离。

表2 VV极化各地物不同时相后向散射系数Tab.2 Backscattering coefficients of VV polarized objects in different time-phase (dB)

图8 VV极化各地物不同时相后向散射系数图

图9 不同时相数据地类VV散射系数箱型图Fig.9 Box diagram of VV scattering coefficient for different time-phase data表3 VH极化各地物不同时相后向散射系数Tab.3 Backscattering coefficients of VH polarized objects in different time-phase (dB)

图10 VH极化各地物不同时相后向散射系数Fig.10 Backscattering coefficient diagrams of VH polarized objects in different time-phase

图11-1 不同时相数据地类VH散射系数箱型图

图11-2 不同时相数据地类VH散射系数箱型图

从VH极化各地物后向散射系数均值来看，各地物后向散射系数在-27～-14 dB之间，其中水体最低，林地和人工表面数值较为接近，且相对其他地类后向散射系数值较高。

从图11中可以看出，水体与其他地类具有较好的分离度，除7月8日和7月26日与湿地较为相似，其他几个时相数据均能较好分离。在VH极化数据中人工表面、林地的后向散射系数比较接近，难以分离。4月9日—5月3日数据莲田与水稻田分离性较低，原因是该时间段，水稻和莲都处于生长初期，地物主要受到水体的镜面反射影响，二者区分度不大，但是6月2日之后，莲田与水稻田的后向散射系数值逐渐拉大， 7月8日和7月26日二者具有较好的区分度。从单时相数据来看，6月14日、7月8日和7月26日水稻田与其他地类具有较好的区分度。

2.4 水稻田提取可分离性分析

2.4.1 Jeffries-Matusita距离

Jeffries-Matusita距离(J-M) 通过计算地类样本之间的特征距离判断地类的可分离程度，是分类过程中样本可分离程度评价的常用指数[29]。采用J-M距离计算地类之间的差异是一种具有通用性的定量方式，它不需要假定不同地物后向散射系数呈正态分布[29-30]。J-M距离计算结果值在0～2之间，值越大可分离性越好，若大于1.8则样本可分离性非常好，在1.4～1.8之间样本选择合格，样本可分离性小于1.4地类之间可分离性较差，建议重新选择样本，可分离性小于1则建议将地类样本进行合并[29-30]。

J-M距离j计算公式[29-30]为：

j=2(1-e-B)，

(1)

式中B为某一特征维上的巴氏距离。

不同地类所选取的样本巴氏距离B计算公式[29-30]为：

(2)

2.4.2 时相优选

为达到理想的分类精度，同时尽可能少地选用数据，降低数据获取和处理成本[31]，研究对时相数据进行优选分析，确定水稻田识别的最佳时相。统计了不同时相水稻田与其他地物的可分离性，结果如表4所示。从单时相上来看，除了4月9日和7月26日的湿地，5月15日和6月14日的水体，7月26日的莲田具有可分离性外，其他时相水稻田与其他地物之间的J-M距离小于1.4，可分离性不佳，证明双极化数据单时相在地类分类过程中存在一定的局限性。

表4 单时相数据水稻田与各地物J-M距离Tab.4 Single-phase data J-M value between paddy filed and other land types

4月9号水稻刚刚出苗，易与杂草相混淆，7月26号处于第二季水稻返青时期，第二季水稻与第一季水稻种植面积有所差异。为了探测第一季水稻的种植范围，根据表4，从4月9日—7月8日第一季水稻生长过程中的SAR影像数据可分离程度最高的时相为6月14日。通过在6月14日垂直极化数据的基础上加上其他极化条件，每次加入的数据以J-M距离作为标准。表5给出了不同时相条件组合条件下水稻与其他地物的J-M距离。结果可以看出，增加时相数据均能增加地类的可分离性。这表明采用多时相数据可有效提高地类的分离度。

表5 多时相数据水稻田与各地物J-M距离Tab.5 Multi-temporal data J-M value between paddy field and other land types

从表5中可以看出，加入7月26日影像数据水稻田与莲田可分离性最大达到1.675； 加入5月3日数据林地与水稻田可分离性最大，达到1.952； 加入4月21日数据人工表面与水稻田可分离性性最大，达到1.786； 加入7月26日数据水稻田与水体和湿地可分离性最大，分别达到1.954和1.998。

将4月21日、5月3日和7月26日数据与6月14日数据进行叠加生成一个包含不同时相的数据文件，计算水稻田与其他地类之间的J-M距离，结果如表6所示。

表6 水稻田与各地物J-M距离Tab.6 J-M value between paddyfield and other land types

采用以上数据进行组合的SAR数据水稻田与其他地类的可分离度均超过1.8，具有较好的分离度。研究区为实现各种地类有效区分，可选择6月14日、4月21日、5月3日和7月26日这4个时相数据。

2.4.3 水稻田识别精度评价

采用总体精度和F测度来分析和评价地物提取的精度，其中F测度包括准确率、召回率和F指数3个指标，准确率指在水稻田识别结果中水稻田被正确识别的比例，召回率指样本中水稻田被正确识别的比例，F指数指准确率P和召回率R的协调平均值[31]，总体精度A是所有验证样本中被正确识别的比例，计算公式为：

(3)

(4)

(5)

(6)

式中：TP为水稻田样本被正确识别的个数；FP为非水稻田样本被识别为水稻田的个数；FN为水稻田样本被识别为非水稻田的个数；TN为非水稻田样本被正确识别的个数[31]。

3 结果分析与验证

3.1 时相优选前的水稻田识别

时相优选前，采用水稻生长所有物候期影像数据进行分类。分类方法包括随机森林法、最大似然法、支持向量机和神经网络法4种。

采用随机森林方法对影像进行分类，需要设置决策树个数n和节点分裂过程中抽取的特征个数m。特征变量m抽取时，默认选择总特征个数的算术平方根作为不同方案的特征个数； 理论上，决策树n的个数越多分类精度越高，但时间成本也越高。有研究发现，决策树n≥100时，所有方案的袋外(out of bag,OOB)误差逐渐收敛并趋于稳定[26,32]。因此选取n=100作为生成决策树的数量。最大似然法核函数选择径向基核函数(radial basis function)，惩罚参数设置100，其他参数默认； 最大似然法概率阈值设置单值(single value)，其他参数默认； 神经网络法训练阈值贡献0.9，训练率0.2，其他参数默认。各种方法执行分类后在ENVI软件中通过主要分析(majority)处理后导出，结果如图12所示。

图12 时相优选前分类结果

随机森林、最大似然法和支持向量机对水稻田识别的总体精度分别为0.937，0.927和0.933，但最大似然法解译的人工表面比实际情况偏大，验证样本中3个水稻田样本被错分为人工表面，2个样本被错分为莲田，支持向量机2个样本被错分为人工表面，1个样本被错分为林地和2个样本被错分为莲田。神经网络法分类法被错分为莲田的面积最大，其水稻田验证样本中1个样本被错分为人工表面，6个样本被错分为莲田。

精度验证结果如表7所示。随机森林的准确率、召回率、F指数和总体精度均为最高，而神经网络分类结果各项指标都要比其他3种方法差，从水稻田识别的初步结果来看，在时相优选前，采用水稻生长各物候期的影像进行分类，随机森林分类方法所得到的结果总体最高，达到0.937。

表7 时相优选前结果精度比较Tab.7 Accuracy comparison of classification resultsbefore time phase optimization

3.2 时相优选结果分析

针对时相优选结果，同时使用随机森林、最大似然法、支持向量机和神经网络4种分类器进行分类并对比结果，解译结果如图13所示。解译结果总体上与图12保持一致。不过由于时相减少，优选后的解译结果有所变化。4种不同分类器对水稻田识别的结果见表8，随机森林、最大似然法和神经网络法解译结果的总体精度比时相优选前的结果均有所提高，说明时相优选后保留的影像数据在水稻田识别冗余信息减少。时相优选后，准确率最高的是神经网络法，达到0.936，随机森林法略低，为0.913； 从召回率上来看最大似然法、支持向量机和神经网络3种分类器在时相优选后均有提升，解译结果的漏检率减少，也说明时相优选后对水稻田识别有用的特征增多。从总体精度上来看，随机森林方法进行水稻田识别效果最佳，达到0.943，Kappa系数0.932。

图13 时相优选后分类结果

表8 时相优选后分类结果精度比较Tab.8 Accuracy comparison of classificationresults after time phase optimization

4 结论与讨论

本文采用多时相Sentinel -1数据，多种分类方法提取南昌县蒋巷镇水稻种植范围，得到如下结论：

1)从J-M距离上来看，单时相Sentinel-1数据用于地物可分离度不高，水稻田容易与其他地类相混淆。Sentinel-1的VV和VH双极化单时相数据在不同的物候期，水稻田与其他地类J-M距离无法全部达到1.4以上。

2)通过增加时相数据可以增加地物的可分离性，选择最优时相数据不但可以降低成本、提升效率，还能降低数据冗余、提升数据解译精度。早稻物候期中采用三叶期、移栽期、孕穗期和二季晚稻移栽期SAR数据组合水稻田与莲田、林地、人工表面和水体的J-M距离分别可达到1.832，1.979，1.890和1.975。解译结果上来看，时相优选后，解译精度普遍得到了提升。

3)随机森林法、最大似然法、支持向量机和神经网络法在水稻田识别的潜力不同，其中无论是时相优选前还是时相优选后，随机森林分类方法的总体精度最高，时相优选后采用随机森林方法解译总体精度达到0.943，Kappa系数为0.932。

本文对比研究了随机森林、最大似然法、支持向量机和神经网络方法对水稻田进行提取，结果表明，随机森林方法提取异质性较强的水稻田具有较高的精度。但是在实际使用过程中不可避免地存在一些局限，如Sentinel-1 数据处理后的空间分辨率为20 m，在小尺度上可能难以满足应用需求，因此后续如果需要提取更加精细的地类分布范围，可以考虑将SAR数据与高空间分辨率数据进行融合再对数据进行提取。