基于Sentinel-1/2影像的水稻种植面积提取方法研究

doi:10.11873/j.issn.1004-0323.2023.1.0078

基于Sentinel-1/2影像的水稻种植面积提取方法研究

——以三江平原为例

乔树亭^,¹^,², 叶回春²^,³, 黄文江²^,³^,⁴, 黄珊瑜⁵, 刘荣豪^,¹, 郭安廷²^,³, 钱彬祥²^,⁴

1.太原理工大学水利科学与工程学院，山西太原 030024

2.中国科学院空天信息创新研究院数字地球重点实验室，北京 100094

3.中国科学院空天信息研究院海南研究院海南省地球观测重点实验室，海南三亚 572029

4.中国科学院大学，北京 100094

5.农业农村部规划设计研究院，北京 100125

Study on Extraction Method of Rice Planting Area based on Sentinel-1/2 Image: A Case Study of Sanjiang Plain

QIAO Shuting^,¹^,², YE Huichun²^,³, HUANG Wenjiang²^,³^,⁴, HUANG Shanyu⁵, LIU Ronghao^,¹, GUO Anting²^,³, QIAN Binxiang²^,⁴

1.College of Water Resources Science and Engineering，Taiyuan University of Technology，Taiyuan 030024，China

2.Key Laboratory of Digital Earth Science，Chinese Academy of Sciences，Beijing 100094，China

3.Hainan Institute，Institute of Aerospace Information Research，Chinese Academy of Sciences，Sanya 572029，China

4.University of Chinese Academy of Sciences，Beijing 100049，China

5.Academy of Agricultural Planning and Engineering，Beijing 100125，China

通讯作者: 刘荣豪（1985-），男，山东临沂人，副教授，主要从事植被定量遥感及应用研究。E⁃mail：liuronghao@tyut.edu.cn

收稿日期: 2022-02-16 修回日期: 2022-10-20

基金资助:

中国科学院战略性先导科技专项“水稻土和白浆土质量与产能提升三江示范区”.  XDA28100500
国家自然科学基金项目.  42001384
中国科学院青年创新促进会项目资助.  2021119

Received: 2022-02-16 Revised: 2022-10-20

作者简介 About authors

乔树亭（1998-），男，山西孝义人，硕士研究生，主要从事植被定量遥感及应用研究E⁃mail：tyutqiaoshuting@126.com , E-mail：tyutqiaoshuting@126.com

摘要

水稻是中国主要粮食作物之一，稻米产量关系到民生福祉。及时、准确地获取水稻种植面积信息及其空间分布状况对于区域农业发展规划和产量评估具有重要意义。针对水稻与其他农作物易混以及光学数据易受云雨天气影响等问题，以东北三江平原为例，利用中高分辨率Sentinel-1微波数据、Sentinel-2光学数据，分别构建时序水体指数SDWI和植被指数NDVI组成水稻完整的物候生长曲线，分析水稻移栽期、分蘖期、抽穗期、成熟期4个重要生长时期不同的光谱差异，通过阈值分割和组合不同时期的数据，来实现水稻不同物候时期种植面积的提取，并与传统的基于单一光学数据的方法进行对比。研究结果表明：经过地表样本点的验证，所构建方法可以精确提取三江平原水稻几个关键生育期的种植面积并且优于单一使用光学数据的方法。同时利用单生育期影像例如移栽期影像提取水稻面积也可使总体精度达到87.08%，随着生育期数据的完整，总体精度也不断提高，其中基于全生育期的面积提取总体精度也高达91.88%，Kappa系数为0.834，可以满足实际应用需求。因此这种的多源数据结合的水稻种植面积提取方法能够准确、高效地提取三江平原水稻不同物候时期种植面积，为短期内的农情调查管理和区域农业可持续发展提供依据。

关键词： 水稻面积提取 ; 物候特征 ; Sentinel⁃2 ; Sentinel⁃1 ; 三江平原

Abstract

Rice is one of the main grain crops in China， and rice yield is related to people's well-being. Timely and accurate acquisition of rice planting area information and its spatial distribution is of great significance for regional agricultural development planning and yield assessment. To solve the problems of rice mixing easily with other crops and optical data being susceptible to cloud and rain weather， taking the Sanjiang Plain in northeast China as an example， a complete rice phenological growth curve was constructed by using time-series water index SDWI and vegetation index NDVI， respectively， based on sentinel-1 microwave data and Sentinel-2 optical data. The spectral differences of four important growth stages of rice were analyzed， including transplanting stage， tillering stage， heading stage and maturity stage， and the planting area of rice in different phenological stages was extracted by threshold segmentation and combination of data of different stages， and compared with the traditional method based on single optical data. The results show that the proposed method can accurately extract the planting area of rice in several key growth stages in Sanjiang Plain， and is superior to the method using optical data alone. At the same time， the overall accuracy of rice area extraction from single growth period images such as transplanting period images can also reach 87.08%. With the completeness of growth period data， the overall accuracy of rice area extraction based on the whole growth period is also as high as 91.88%， and the Kappa coefficient is 0.834， which can meet the requirements of practical application. Therefore， the rice planting area extraction method combined with multi-source data can accurately and efficiently extract the rice planting area in different phenological periods in Sanjiang Plain， and provide a basis for short-term agricultural situation investigation and management and regional agricultural sustainable development.

Keywords： Rice area extraction ; Phenological characteristics ; Sentinel-2 ; Sentinel-1 ; Sanjiang plain

PDF (5909KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

乔树亭, 叶回春, 黄文江, 黄珊瑜, 刘荣豪, 郭安廷, 钱彬祥. 基于Sentinel-1/2影像的水稻种植面积提取方法研究. 遥感技术与应用[J], 2023, 38(1): 78-89 doi:10.11873/j.issn.1004-0323.2023.1.0078

QIAO Shuting, YE Huichun, HUANG Wenjiang, HUANG Shanyu, LIU Ronghao, GUO Anting, QIAN Binxiang. Study on Extraction Method of Rice Planting Area based on Sentinel-1/2 Image: A Case Study of Sanjiang Plain. Remote Sensing Technology and Application[J], 2023, 38(1): 78-89 doi:10.11873/j.issn.1004-0323.2023.1.0078

1 引言

中国既是粮食生产大国又是消费大国，粮食安全不仅影响社会稳定和经济发展，更关系国家长治久安^［1］。三江平原位于黑龙江省东北部，不仅是我国黑土地宝贵资源的重要组成部分，还是我国粮食安全关键地区和优质商品粮基地^［2］。区内水资源丰富，人均耕地面积大致相当于全国平均水平的5倍^［3］，水田面积266.7万hm²，旱田346.7万hm²，生态环境条件好，被誉为“北大仓”^［4］。水稻是三江平原的主要农作物之一，其播种面积在所有粮食作物中约占33%^［5］。因此，及时、准确地获取水稻种植面积信息及其空间分布状况对于对区域农业精细化生产管理与产量评估具有重要意义^［6］。

传统方法获取水稻种植面积依赖于统计上报与地方实测，缺少空间分布的刻画且耗资巨大^［7］。遥感技术因覆盖范围大、探测周期短、成本低等优势，已被广泛应用于大范围的种植类型识别和植被监测等相关研究，也为高效获取水稻种植面积提供了方法^［8］。总体上，基于遥感数据的水稻种植面积提取方法根据输入数据的不同，可以分为基于光学数据和微波数据两种方法。基于单一光学数据的遥感方法主要是分析某一光谱特征的时间序列剖面，找出从中提取出水稻的关键物候期（如灌水移栽期），然后确定区分特定作物的光谱特征阈值来识别水稻与其他作物^［9］。随着免费的中高分辨率卫星的迅速发展，欧空局发射的Sentinel-2多光谱卫星以其10 m的空间分辨率和5 d的重访周期被广泛应用到水稻面积提取中来。如姜依兰等^［10］以2019~2020年杞县地区多时相Sentinel-2影像为数据源利用不同作物在时间序列NDVI曲线上的差异构建NDVI指数来提取水稻种植区；Amir等^［11］利用Sentinel-2影像的时间序列检测水稻，研究了水稻栽培时期近红外波段反射率的变化规律以及时间归一化差异植被指数数据检测水稻的方法，并将其与其他作物区分开来。研究表明：在光学数据能有效覆盖整个生长季的前提下，基于上述的分类方法能保证一定的精度^［12］。但是，在一定程度上云雨等因素会严重影响光学遥感数据质量，对于重访周期较长的中、高分辨率数据的影像更为显著，可能造成1~2个月的数据缺失^［7］。因此，单纯地利用光学遥感数据难以实现对水稻关键物候信息（特别是灌水移栽期）的准确提取。

基于单一微波数据的遥感方法主要是根据物候特征分析水稻的散射机制进而提取水稻信息。Sentinel-1是欧空局发射的地球观测卫星，由两颗卫星组成，是目前唯一能够免费获取数据，并且定期发布影像的中高空间分辨率的合成孔径雷达，未来将成为微波遥感在各个研究领域的重要数据源。其中Son等^［13］通过利用时间序列Sentinel-1合成孔径雷达（SAR），构建平滑的12 d时间序列VH偏振复合数据，然后利用物候指标对水稻种植区进行基于对象的图像分类并评价结果的精度评估；黄翀等^［14］基于Sentinel-1 SAR时间序列数据，提出一种融合时序统计参数与时序曲线相似性特征的热带地区水稻种植结构提取方法。基于上述研究不难得出基于单一微波数据提取水稻面积最大的问题在于需要多时期、长时间序列的影像，对于短期内的农情调查则很难发挥其作用。

综上可知，单一遥感数据往往不足以满足作物分类的要求，现阶段光学与遥感数据结合来提取水稻种植面积的方法逐渐开始发展。Cai等^［15］基于Sentinel-2 NDVI时间序列、物候数据和时间序列基于JBh距离的Sentinel-1 SAR后向散射图像，使用随机森林分类器利用最优特征组合对水稻进行提取；Fiorillo等^［16］基于Sentinel-1和Sentinel-2影像使用随机森林算法绘制Sédhiou地区的低地水稻作物区域，并对2017~2019年的生长季节进行分析。因此，选择合适的数据组合和方法来减少云污染对水稻关键时期的影响是当前分类研究的重点^［17］。

以三江平原为研究区域，基于Google Earth Engine平台，耦合Sentinel-1微波数据和Sentinel-2光学数据来构建作物完整物候曲线，利用阈值分割法来提取水稻种植区面积，主要目标为：①与基于单一光学数据的水稻面积提取方法进行对比，降低在提取水稻移栽期面积时光学数据容易受到云雨天气的影响。②基于水稻4个关键生育期进行不同的数据组合来提取相应时期的水稻面积，为短期内的农情调查和管理提供依据。

2 数据与研究区

2.1　研究区概况

三江平原位于中国东北黑龙江省的东部，土地总面积为10.9×10⁴km^2［2］。包括黑龙江、松花江、乌苏里江的冲积平原以及完达山以南由乌苏里江和兴凯湖冲积而成的低地平原^［18］。三江平原属于温带湿润、半湿润大陆性季风气候，平均海拔50~60 m，年平均气温1~4 ℃，7月的平均气温为21~22 ℃，适合作物生长，年降水量为500~650 mm，降水主要集中在6~10月^［19］。同时也是中国重要的商品粮基地，年产量约1 500万t^［20］。人均耕地面积和人均粮食产量是全国水平的5倍^［21］。三江平原主要农作物为水稻、玉米和大豆，其中玉米和大豆种植面积各占其农作物总种植面积的51.3%和23.4%，而水稻种植面积约占其所种植农作物总面积的24.7%^［5］，是世界上最大的优质粳稻种植区之一，种植方式以一年一熟为主^［22］，同时3种作物的生育期都约为4月底至10月底，且生长曲线都较为接近。图1为三江平原土地利用数据^［23］、位置及样本点分布图。表1为三江平原主要农作物生育期情况汇总。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 三江平原土地利用数据、位置及样本点分布

Fig.1 Land use data， location and sample point distribution in Sanjiang Plain

表1 三江平原主要农作物生育期情况汇总

Table 1 Summary of growth period of main crops in Sanjiang Plain

月份		4			5			6			7			8			9			10
旬	上	中	下	上	中	下	上	中	下	上	中	下	上	中	下	上	中	下	上	中	下
水稻		播种		灌水移栽			苗期		分蘖			抽穗				成熟			收获
玉米			播种		出苗				拔节			抽穗		开花		灌浆		成熟		收获
大豆			播种		幼苗			分化			开花			结荚			成熟		收获

新窗口打开| 下载CSV

2.2　数据获取及预处理

2.2.1　地面验证数据

为了保证足够的数量和质量用于分类及精度验证，本研究的样本点来源于野外调查数据和Google Earth二者相结合的总体数据。其中地面样本点来源于2019年7月份作物生长季的实地调查工作，样本数据主要记录了作物的详细信息，主要包括农作物类型、物候期、生长状况和半径5 m定位精度的GPS采集的经纬度信息。通过整体野外调查数据与Google Earth结合的总体数据，得到有效的样本数据813个。其中水稻样本点330个，玉米样本点375个，大豆样本点108个。按照3∶1的比例将所有样本点随机分为用于本文水稻、玉米和大豆各指数时间序列曲线的样方点以及水稻面积提取、精度评价的分类精度验证点。各样本采样点位置如图1所示。

2.2.2　Sentinel-1/2卫星数据及处理

“Sentinel”系列卫星是欧洲哥白尼（Copernicus）计划空间部分（GSC）的专用卫星系列，由欧洲委员会（EC）投资，欧洲航天局（ESA）研制^［24］。而其中的Sentinel-1、Sentinel-2卫星数据数据用于本研究并在GEE平台中存档，二者相结合，可以获取大范围、高重访周期的数据。提供高分辨率的环境监测能力^［25］。

Sentinel-1卫星由两颗卫星组成，分别为Sentinel-1A和Sentinel-1B，载有C波段合成孔径雷达，可提供连续图像（白天、夜晚和各种天气），分辨率为10 m，每6 d对整个地球进行一次成像，提供的极化方式有单极化和双极化^［27-28］。垂直极化后Sentinel-1数据调用GEE平台中的预处理模块完成标定轨道参数、去除热噪声、辐射定标、地表几何纠正等操作。

Sentinel-2由两颗卫星组成，重访周期为5 d，携带多光谱成像仪共13个波段，其中4个波段的分辨率为10 m，6个波段分辨率为20 m，3个波段的分辨率为60 m，光谱覆盖范围从可见光到近红外至短波红外^［26］。根据GEE平台上官方发布的Sentinel-2影像数据说明，Sentinel-2遥感影像产品等级包括1C级和2A级，其中2A级数据是已经经过大气校正的大气底部反射数据可以直接运用到本研究中。同时考虑到云雾对数据的影响，可以采用GEE平台自带的QA60波段去除云像元^［10］。

2.2.3　统计数据

实际水稻面积根据《黑龙江省2020年统计年鉴》、《黑龙江垦区2020年统计年鉴》和依兰县、穆棱市农业统计数据来获取2019年三江平原水稻种植面积统计数据，以此对农作物提取面积进行精度验证。

3 研究方法

本研究中基于多源遥感数据进行水稻面积提取的原理是根据物候生长曲线，将水稻按主要生长时期分为移栽期、分蘖期、灌浆期和成熟期4个时期，首先基于水体指数提取灌水移栽期水稻面积，然后通过时序特征NDVI等值进行阈值分割，选择不同的数据组合来获取相应生育期的水稻面积，最终通过对比单一光学数据提取水稻面积和微波与光学数据相结合来提取水稻面积两种方法并验证其精度来进行比较，主要流程如图2所示。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 提取水稻种植区面积技术流程图

Fig.2 Extract technical flow chart of rice planting area

3.1　光学数据提取水稻种植面积方法

移栽蓄水期的水稻生长在被淹没的土壤中，稻田是开放水域和绿色水稻植物的混合物，具备更易识别的地物特征，对于水稻的面积提取是最重要的一个时期^［29］。所以本研究提取水稻面积最重要的部分就是基于不同的水体指数提取灌水移栽期水稻面积。其中基于单一光学数据的方法主要是采用Sentinel-2卫星数据通过不同时期的数据组合来提取水稻相应生育期的种植面积。其中，在提取水稻灌水移栽期面积时本研究采用了归一化差异水体指数（MNDWI）进行阈值分割来提取水田部分作为水稻移栽期面积的基础，因为该指数采用了归一化的比值运算，可以一定程度上消除地形差异的影响，从而解决了水体信息中有夹有阴影的问题，更适合提取水体信息^［30］。其中基于训练样本绘制的三江平原3种主要作物的MNDWI时序曲线如图3所示，且MNDWI指数公式如下：

M N D W I = (B 3 - B 11) / (B 3 + B 11)

(1)

其中：B3和B11分别为Sentinel-2卫星的绿波段和短波红外波段。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 时序MNDWI曲线

（已经过去云筛选及SG滤波处理，标注为显著置信区间）

Fig.3 Time-series MNDWI curve

归一化植被指数（NDVI）作为反映植被生长状态和覆盖度的指数，其值的增长或降低与作物生长变化过程相对应^［10］。其中基于训练样本绘制的三江平原3种主要作物的NDVI时序曲线如图4所示，且NDVI指数公式如下：

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 时序NDVI曲线

（已经过去云筛选及SG滤波处理，标注为显著置信区间）

Fig.4 Temporal NDVI curve

N D V I = (B 8 - B 4) / (B 8 + B 4)

(2)

其中：B4和B8分别为Sentinel-2卫星的红波段和近红外波段。

由平滑后3种作物的MNDWI及NDVI时序曲线可知，3种作物基本都于4月至5月之间种植，水稻由于其灌水移栽期的特性，NDVI值在种植前期一直小于-0.1，但是与其他作物区别并不大，随着灌水移栽期的结束，水稻进入分蘖期，生长旺盛，其NDVI值在逐步增大到0.4以上，然后于8月份之后达到峰值0.7左右，此后NDVI值逐步下降，到10月份前后水稻完成收割，NDVI变为最小值（<0.3）。而其MNDWI值在4月份开始不断增大，5月份前后达到峰值于0.1~0.3之间浮动，明显大于另外两种作物，然后随着分蘖期的到来逐渐降低和其他作物逐渐接近。

由上述分析可以得出，水稻的NDVI曲线与其他作物并没有明显区别，但是可以反映水稻生长发育的全过程。而MNDWI指数在灌水移栽期可以很好地将水体部分与其他作物区分开。所以研究所涉及的使用传统的根据光学数据提取水稻种植区面积的方案①为使用MNDWI指数来区分水稻与玉米、大豆，将水田面积提取出来。然后在方案①的基础上通过计算不同时期的NDVI指数进行阈值分割来提取另外不同时期的面积，其中在方案①基础上加入分蘖期的NDVI数据为方案②，在方案②基础上加入抽穗期的NDVI数据为方案③，在方案③基础上加入成熟期的NDVI数据为方法方案④。具体数据组合如表2所示。

表2 基于Sentinel-2数据的不同生育期水稻种植面积提取方案

Table 2 Extraction scheme of rice planting area at different growth stages based on Sentinel-2 data

日期	方案
日期	①	②	③	④
4月20日~5月20日	Sentinel-2 MNDWI (移栽期)	Sentinel-2 MNDWI (移栽期)	Sentinel-2 MNDWI (移栽期)	Sentinel-2 MNDWI (移栽期)
6月20日~7月20日		Sentinel-2 NDVI (分蘖期)	Sentinel-2 NDVI (分蘖期)	Sentinel-2 NDVI (分蘖期)
7月20日~8月20日			Sentinel-2 NDVI (抽穗期)	Sentinel-2 NDV I(抽穗期)
9月20日~10月20日				Sentinel-2 NDVI (成熟期)

注：方案①②③④为单一光学数据组合，①为移栽期数据提取水稻面积；②为移栽期+分蘖期数据提取水稻面积；③为移栽期+分蘖期+抽穗期数据提取水稻面积；④为4个生育期全数据提取水稻面积

新窗口打开| 下载CSV

3.2　微波数据与光学数据相结合提取水稻种植面积方法

微波数据与光学遥感数据结合提取水稻的方法主要使用Sentinel-1微波数据、DEM高程数据与Sentinel-2光学数据等。SAR通过发射电磁波照射到地物表面，再接受地物的后向散射，获得地物的电磁特性，平静的水面会发生镜面反射，当水面有波浪或不平整时，表面粗糙度就会增大，会发生漫反射^［31］，所以我们选用在水稻灌水移栽期结束之后10 d左右的出苗期影像，会使水体表面发生漫反射更容易凸显出水体。首先在提取移栽期水稻面积时，使用Sentinel-1 SAR数据计算得出的SDWI指数（Sentinel-1 Dual-Polarized Water Index）^［27］来提取水稻种植区面积，将其命名为方案⑤。其中基于训练样本绘制的三江平原3种作物的SDWI时序曲线如图5所示，且SDWI指数计算公式如下：

S D W I = l n (10 ∙ V V ∙ V H) - 8

(3)

其中：SDWI是水体提取指数，当其值大于0时为水体，小于0时为非水体；VV和VH为双极化方式，8是本文所选择的最初阈值^［27］。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 时序SDWI曲线

（已经过SG滤波平滑处理且标注为显著置信区间）

Fig.5 SDWI curve of timing sequence

由于SAR数据质量很大程度受观测角度影响，即距成像中心距离越大的像元噪声越强烈。因此，预处理后的SAR数据以月为单位，逐像元计算其信号在本时间段内的中位数，以排除信号对成像角度及地表状况过敏感带来的噪声^［7］。

在使用水体指数来提取水稻灌水移栽期种植面积时，光学数据一定程度会受到云雨天气的影响，需要进行筛选，而微波数据则可以避免这一问题，由图5可知，当4月水稻灌水移栽期开始时，其SDWI值开始突然上升（>0.1），与其他两种作物发生明显差异，在五月份开始到达最大值（>0.3），然后在6月份开始逐渐降低，与其他两种作物差别并不是十分明显。所以方案⑤就是利用微波水体指数，通过阈值分割来区分水稻与玉米、大豆，然后提取水稻灌水移栽期面积。然后在移栽期水稻面积的基础上使用Sentinel-2数据得出的NDVI值继续进行阈值分割来提取剩余生育期面积，与前一种方法相同并以此类推可得出为方案⑥、⑦、⑧。具体数据组合如表3所示。

表3 基于Sentinel-1和Sentinel-2数据的不同生育期水稻种植面积提取方案

Table 3 Extraction scheme of rice planting area at different growth stages based on Sentinel-1 and Sentinel-2 data

日期	方案
日期	⑤	⑥	⑦	⑧
5月20日~6月20日	Sentinel-1 SDWI(移栽期)、DEM	Sentinel-1 SDWI(移栽期)、DEM	Sentinel-1 SDWI(移栽期)、DEM	Sentinel-1 SDWI(移栽期)、DEM
6月20日~7月20日		Sentinel-2 NDVI(分蘖期)	Sentinel-2 NDVI(分蘖期)	Sentinel-2 NDVI(分蘖期)
7月20日~8月20日			Sentinel-2 NDVI(抽穗期)	Sentinel-2 NDVI(抽穗期)
9月20日~10月20日				Sentinel-2 NDVI(成熟期)

注：方案⑤~⑧为微波数据与光学数据相结合；⑤为利用移栽期微波数据提取水稻面积；⑥为移栽期+分蘖期数据提取水稻面积；⑦为移栽期+分蘖期+抽穗期数据提取水稻面积；⑧为4个生育期全数据提取水稻面积。

新窗口打开| 下载CSV

3.3　验证方法

在利用水稻面积提取方法获得三江平原水稻种植面积和空间分布信息的基础上，利用地面实测数据和统计年鉴对2019年三江平原水稻提取结果进行精度验证。在定量统计方面，通过地面验证点建立混淆矩阵，计算其用户精度（User Accuracy，UA，%）、生产者精度（Producer Accuracy，PA，%）、总精度（Overall Accuracy，OA，%）和Kappa系数，其中Kappa系数分为5组来表示不同级别的一致性：0.0~0.20为极低、0.21~0.40为一般、0.41~0.60为中等、0.61~0.80为高度的和0.81~1为几乎完全一致。其计算如式（4）~（7）所示。

U A = \frac{n_{i i}}{n_{i ∙}} \times

100%(4)

P A = \frac{n_{i i}}{n_{∙ i}} \times

100%(5)

O A = \frac{\sum_{i = 1}^{q} n_{i i}}{n} \times

100%(6)

K a p p a = \frac{[n ∙ \sum_{i = 1}^{q} n_{i i} - \sum_{i = 1}^{q} (n_{i ∙} ∙ n_{∙ i})]}{[n^{2} - \sum_{i = 1}^{q} (n_{i ∙} ∙ n_{∙ i})]}

(7)

其中：n_ii 为混淆矩阵中第i行第i列的值；n_i· 为混淆矩阵中第i行的和；n_·i 为混淆矩阵中第i列的和；n为验证样本总数；q为混淆矩阵行列数。

4 结果与分析

4.1　水稻提取特征敏感性分析

经过上文所述，本研究根据水稻生育期内与其他作物的物候生长曲线来设置判别规则，考虑到不同下垫面之间的差异，本文将作物类型主要分为水稻和非水稻（玉米、大豆），基于训练样本绘制的SDWI、MNDWI和NDVI这3个主要分类参数的时序曲线，如图3~图5所示。从图中可以看出，在水稻生长期，几种典型的农作物的NDVI值与水稻NDVI值时间序列曲线有着一定差异但不明显，其中在4月底至5月初水稻灌水移栽期的时候，水稻NDVI值略小于其他两种作物且小于0，约为-0.1，因为之间的差异并不明显，所以不能很好地区分水田和其他作物。同时由于稻田的特殊性，所以使用水体指数来提取稻田面积更有优势。其中在单一使用常规的高精度光学影像Sentinel-2来计算水体指数提取水稻移栽期面积时，根据图3（b）设置MNDWI_①值范围为（0.1，0.3）可以将永久性水体和其他地物排除在外。根据NDVI生长曲线，在使用方案②时，NDVI_②阈值设置为大于0.4；在使用方案③时，作物的NDVI_③值随着抽穗期的结束，此时已达到峰值，其阈值为大于0.6，同时设置NDVI_③大于NDVI_②。最后在使用方案④计算水稻成熟期面积时，在方案③提取面积的基础上计算成熟期水稻NDVI平均值，因为随着水稻被收割之后，NDVI值又趋于裸地，所以设置判别条件NDVI_④小于NDVI_③，同时NDVI_④小于0.3。基于上述条件，分别使用方案①、②、③、④可以得出4个时期的水稻种植面积，且方法中的各指数划分区间都尽量使用30 d为单位，来减少同种作物物候的影响。

光学与微波数据相结合的方法与单一使用光学数据的方法最大的不同在于使用水体指数作为阈值提取的基础，所以本研究使用Sentinel-1 SAR微波数据计算得出的SDWI水体指数来代替方案①中的MNDWI指数。提取水稻灌水移栽期面积主要选取Sentinel-1卫星5月20日至6月20日的数据，计算得出水体指数SDWI之后，根据图3（b）选取阈值为0.1~0.4剔除掉非水体和永久性水体，根据地方标准《水稻种植技术规范》再选用DEM高程数据选择小于280 m高程地形来减少SAR数据在山体阴影等地的影响。在此基础上可以得出移栽期水稻的面积。在使用方案⑥、⑦、⑧面积时与基于光学数据的方案完全相同，均采用Sentinel-2数据并计算NDVI，采用相同时间的影像、处理和阈值分割。基于上述方案可以分别得出4个时期的水稻种植面积。

4.2　不同遥感方法水稻种植区提取结果分析及精度评价

根据上述方法进行三江平原水稻面积提取，其中基于单一光学数据的4种方案提取的水稻种植面积和基于微波和光学数据相结合的4种方案的水稻种植面积分别为图6所示。从整体水稻空间分布方面来看，46°~48° N是三江平原水稻的主要分布区域，西南部长白山和东部的完达山地势高，土质较差，农田面积较小、收割较困难、不易大规模种植作物。大规模的水稻种植区分布与北部黑龙江和松花江相邻的大面积平原处，该区域气温高降水足，沿河水塘密集，灌溉充分，腐殖质土壤营养丰富，适宜种植水稻。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 不同方案提取水稻各生育期面积结果

Fig.6 Extraction results of rice area at different growth stages by different methods

基于2.2.1节获取的地面验证数据（包括83个水稻样本点，94个玉米样本点和27个大豆样本点）计算混淆矩阵，并根据统计年鉴进行定量统计，结果如表4所示。从基于光学数据和基于微波和光学数据相结合的水稻面积提取结果来分析可得出：

表4 不同方案水稻种植面积验证精度对比

Table 4 Comparison of rice planting area validation accuracy among different schemes

方案	面积精度 /%	生产者精度/%	用户精度 /%	总体精度 /%	Kappa 系数
①	87.67	77.06	90.61	85.24	0.702
②	92.03	79.05	90.30	86.35	0.723
③	94.30	79.89	90.30	86.84	0.733
④	89.49	85.38	88.48	89.18	0.776
⑤	83.21	77.64	95.76	87.08	0.742
⑥	88.83	80.00	95.76	88.56	0.770
⑦	96.46	81.61	95.45	89.42	0.786
⑧	90.45	86.46	94.85	91.88	0.834

注：表4中的方案①~④、⑤~⑧指上文3.1和3.2中的不同方法组合，其中方案①~④为单一光学数据组合，①为移栽期数据提取水稻面积；②为移栽期+分蘖期数据提取水稻面积；③为移栽期+分蘖期+抽穗期数据提取水稻面积；④为4个生育期全数据提取水稻面积；方案⑤~⑧为微波数据与光学数据相结合；⑤为利用移栽期微波数据提取水稻面积；⑥为移栽期+分蘖期数据提取水稻面积；⑦为移栽期+分蘖期+抽穗期数据提取水稻面积；⑧为4个生育期全数据提取水稻面积。

新窗口打开| 下载CSV

（1）从方案①和方案⑤两个提取灌水移栽期的水稻种植面积到方案④和方案⑧两个提取全生育期的水稻种植面积都保持一定的相关性，水稻空间分布基本一致，并且随着生育期的完善水稻种植面积精确度也越来越高。其中两种基于全生育期的方案④和⑧的总体精度高达89.18%和91.88%，Kappa系数为0.776和0.834。说明两种方法分类精度均很高，但是微波与光学数据相结合的方法提取的水稻种植面积与地面验证数据之间的一致性更强，更适合于实际应用。

（2）基于光学数据的方法还是容易受到云雨天气的影响，尤其是对于水稻来说，影响到灌水移栽时期的面积提取就直接影响到了整个生育期的结果。如图7（b）、图7（c）所示，在研究区的右下角及兴凯湖的上方明显可见由于影像的缺失导致面积提取的不完整，一大片区域因为云污染导致空白，导致每个生育期都缺失这一部分的水稻面积且非常显著，一定程度上影响到了该方法的面积精度。同时，由图7（e）、图7（f）可知，在同样的位置微波与光学相结合的方法可以一档程度上避免云雨天气的干扰（主要是灌水移栽时期的水稻种植面积不会受到影响。）

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 基于全生育期数据的方案④与方案⑧的水稻面积提取结果细节对比图

Fig. 7 Detailed comparison diagram of rice area extraction results based on whole growth period

4.3　水稻不同生育期组合方法面积提取精度评价

两种方法的生产者精度、总体精度和Kappa系数都随着生育期的完善越来越高，其中根据年鉴所得三江平原水稻种植面积为3 731 977.4 m²，各方案与年鉴相对比的面积精度也均高于80%。其中：

（1）在移栽期提取水稻种植面积时，方案⑤的总体精度高达87.08%，说明单生育期影像在提取水稻面积时也达到了一定的精度，证明数据提取移栽期水稻面积也具有一定的优势。

（2）在水稻分蘖期、抽穗期提取水稻种植面积时，方案⑥、⑦的总体精度都高达88.56%和89.42%，且Kappa系数也分别为0.770和0.786，属于与地面验证数据高度一致性。证明在短时间内没有全生育期影像的情况下使用该方法提取水稻种植面积也可以保证相对较高的精度。

（3）但基于4个主要生育期提取水稻面积时，微波与光学相结合的方法中的用户精度也都明显高于基于单一光学数据的方法，说明误判误差比其要更低，但是二者生产者精度差距并不大，说明漏判误差都很接近，证明微波数据和光学数据相结合提取水稻不同生育期的种植面积比使用单一光学数据更具有优势。

5 讨论

当前，遥感提取水稻种植面积的方法有很多种，前人的研究中并没有普遍适用且精度高且精度最高的一种方法，水稻种植信息的提取本质上是结合了物候信息的遥感图像分类问题，常规的基于光学数据的水稻分类方法在实际的水稻主产区的分类工作中极易受到多云天气的影响尤其是对于最重要的灌水移栽期，无法形成有效的光谱时间序列，影响进一步提取水稻关键生长期的特征^［32］。Guo等^［33］利用全极化Radarsat-2数据模拟了多时相CP SAR数据然后利用height-sensitive指数（RHSI）并建立了决策树方法基于RHSI的最优CP参数，在此基础上进行了水稻分类取得了较高的精度；Chang等^［34］建立一种基于利用完整时间序列SAR后向散射函数拟合出水稻生长季节的模型提出了提高水稻检测精度的方法，总体精度可达91.3%；谭深等^［7］在线调用中分辨率的光学、微波遥感数据，创新性地采用了按月提取、按直方图大小提取特征的方式，采用随机森林分类器，绘制海南省2016年10 m分辨率水稻种植分布图，总体精度达到了93.2%，可以满足实际应用需求。但是这些方法大多数都需要完整生育期的影像才能达到一定的精度并且工作量巨大。

而本研究采用微波和光学数据结合的方式，基于GEE平台，利用不受云层影响的Sentinel-1微波数据并结合Sentinel-2数据合成NDVI时间序列数据集组成水稻完整的物候生长曲线。使用阈值分割法，在使用微波水体指数来提取移栽期稻田面积的基础上，利用不同生育期的NDVI数据进行组合来得出不同生育时期的水稻种植面积，并与单一使用光学数据的方法进行对比。最后利用地面验证数据与统计年鉴进行精度验证，得出总体精度高达91.88%，Kappa系数为0.834。证明该方法可以结合实际应用于北方单季稻区并能保证相对能以较高的精度来提取水稻种植区面积。

并且本研究在水稻重要的4个生育期（移栽期、分蘖期、抽穗期和成熟期）提取种植面积时都可以保证高于80%的总体精度。同时，研究最大的优势就在于在保证精度的前提下，在水稻生育期的某一时间段，不需要完整的生育期影像也可以很好地提取出水稻种植区的面积。比如利用微波数据提取移栽期水稻种植面积时，精度就可以达到87.08%，Kappa系数也达到了0.742；在7月初的分蘖期提取水稻面积时，总体精度可以达到88.56%，Kappa系数为0.770。这些都可以证明本研究可以在水稻的关键生育期提供种植区提取结果，为区域农情管理提供第一手资料，支撑三江平原规模化、智能化、现代化农业生产。

本研究还存在一定的不足之处，在Sentinel-1微波数据处理时尽管使用DEM数据来去除一部分的阴影，但是因为DEM数据分辨率为30 m，不能很好地与本文10 m的分辨率叠加，其次微波数据的处理还有进一步提升的空间。在除移栽期水稻面积部分外，其余生育期的面积提取也一定程度上需要光学数据的辅助，需要通过后续研究来弥补云雨天气的影响和进一步提高水稻种植面积的提取精度。而且近年来深度学习在自然场景分类的成就吸引了遥感领域的广泛关注^［35］。该方法表达、自学习及容错能力突出，能够分层、高效且准确地从遥感图像中提取具有代表和区分性的特征^［36］。由于图像质量和数据集规模等限制，深度学习针对遥感影像尤其SAR的水稻等作物提取仍处于挖掘阶段，未来的研究重点将放在这方面^［37］。

6 结论

本研究将Sentinel-1微波数据与Sentinel-2光学数据相结合，通过水稻物候生长曲线进行不同的组合，使用阈值分割方法来提取水稻不同时期的种植区面积。得出以下结论：

（1）与单一使用Sentinel-2时间序列NDVI数据进行面积提取的方法来进行对比，证明Sentinel-1雷达数据在提取水稻种植面积有很好的应用潜力并且微波水体指数在提高水稻识别提取精度方面存在一定的优势。

（2）证明了使用微波数据与光学数据相结合提取水稻种植区面积的可行性程度，虽然在完整生育期下提取水稻种植区面积的精度固然最高，但是在水稻的某一生长期需要面积数据时本方法也可以提供较高精度的结果，对于区域农情管理可以提供第一时间的资料，在北方单季稻种植区可以普遍推广应用。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

LIN

Zhihui

， LIU

Xianfeng

， CHEN

Ying

， et al.

Water-food-energy nexus：Progress，chanllenges and prospect

［J］. Acta Geographica Sinica， 2021，76（7）：1591-1604.