王利民 杨玲波 刘 佳 杨福刚 姚保民

中国农业科学院农业资源与农业区划研究所, 北京 100081

作物长势指作物的生长状况和趋势[1], 可以为作物田间管理[2]、早期产量快速估算[3]、国家粮食政策制定[4]提供支持。由于遥感技术的时效性、便捷性、准确性等优势, 利用遥感技术进行农作物长势监测已成为当前的作物长势的主要获取手段[5-7]。20世纪80年代, 美国农业部建立了全球尺度的农情遥感监测系统, 实现全球主要农作物的长势监测和产量预报, 加拿大、澳大利亚、欧盟等发达国家也相继建立了作物长势遥感监测运行系统[8]; 中国对作物长势遥感监测的研究及应用稍晚于国外, 中国科学院、中国气象局、中国农业科学院、农业部规划设计研究院等建立了一系列的农情遥感监测系统并投入实际应用[9]。当前, 作物长势遥感监测的常用方法主要包括植被指数关联法和作物生长模型两大类[10]。遥感指数关联方法主要利用遥感影像获取与作物生长密切相关的植被指数, 并构建植被指数与长势之间的相互关系, 从而反演地面作物的生长情况或进行长势的分级[11-12], 常用到的遥感指数包括NDVI[13]、LAI[14]、VCI[15]、NDWI[16]等。Moriondo等[17]利用 NDVI作为作物长势指标, 统计分析多年的冬小麦长势规律, 并对意大利两个省份小麦的长势和产量进行预测, 取得了较高的精度; Zhao等[18]利用 NDVI构建了植被生长动力学指数, 通过将NDVI曲线坐标归一化, 根据有效积温划分作物不同生长阶段及对产量的影响因子, 对研究区进行大豆估产, 结果表明精度高达88.5%。作物生长模型同化方法则利用遥感数据与作物生长模型如WOFOST(world food studies)、DSSAT (decision support system for agrotechnology transfer)、SWAP (soil, water, atmosphere, and plant)等结合, 模拟作物整个生育期的生长环境, 并对长势和产量预估计算, 具有机理性强的优势。Maas[19]利用遥感影像获取叶面积指数LAI和水分胁迫系数, 并将其输入作物生长模型,获取了玉米地上生物量的模拟结果; 刘峰等[20]设计实现了基于极快速模拟退火算法的遥感数据与CERES-Wheat作物生长模型的同化原型系统, 以地面高光谱为遥感数据源, 同化LAI并与实测值对比,展现了良好的拟合度。

同时, 在当前的业务系统中, 利用 NDVI等指数进行作物生长状况分级评价, 或与多年均值对比分析, 获得监测区域作物长势[21-23]已成为主要方法。如王蕊等[24]通过对多年冬小麦 MODIS NDVI数据进行比较和分级, 评价了河北平原农田的生产力空间变异; 冯美臣等[25]对相邻年份冬小麦关键生育期 MODIS NDVI 值进行比值计算, 根据比值的大小将冬小麦长势分为 5个等级; 乌兰吐雅等[26]利用MODIS计算NDVI, 以NDVI值的大小对大兴安岭西麓春小麦苗青长势分析评价, 较好地反映了春小麦的长势状况。

MODIS数据具大幅宽(2330 km)、高重访(一日4次)、数据免费获取等优势, 是当前作物长势遥感监测的主要数据源, 但由于其分辨率较低(250 m), 主要应用于大范围区域监测, 如梁瀚月等[27]利用MODIS长时间序列NDVI, 通过构建植被条件指数、距平植被指数、与往年比较指数等, 进行西藏地区大面积作物的苗情监测; 而随着高分数据的增多、数据处理能力的增强, 利用高分辨率数据进行作物长势监测的研究日益增多, 监测范围也从小尺度[28]向大尺度转变[29], 方法从单纯的植被指数长势评价[30]到植株氮素、叶绿素含量等作物物理、生理参数建模反演方向发展[31], 如查海涅[32]基于 GF-1影像反演水稻的氮营养指数、谭昌伟等[33]基于HJ卫星影像进行了冬小麦开花期生物量、叶绿素含量、蛋白质含量等生长指标的遥感监测。

随着国产高分系列卫星的不断发射, 目前中国农情遥感作物长势业务化监测的主要数据源由MODIS等国外中等分辨率数据逐渐转向GF-1等国产高分辨率卫星数据。本文以河北廊坊为主要研究区, 以冬小麦为目标作物, 分别使用GF-1和MODIS数据反演 NDVI作为作物遥感长势监测指标, 同步采集研究区冬小麦样本点, 以综合茎数、株高、叶绿素浓度等作为地面长势监测指标, 分析遥感与地面长势监测指标的关联性; 同时对比分析 GF-1和MODIS影像进行作物长势遥感监测的差异性, 评价作物长势不一致性、传感器差异、NDVI尺度效应等对作物长势的影响情况。本研究将有利于厘清遥感长势监测代表的实际含义及精度, 同时也将为中国当前农情遥感作物长势业务化监测的数据源由MODIS等国外中低分辨率数据向GF-1等国产高分数据转变提供依据。

1 研究区概况

研究区大部分位于河北省廊坊市(39°28′42″～39°32′54″N, 116°38′07″～116°44′06″E), 少部分位于天津市武清区和北京市大兴区, 总面积约 4500 km2(图1)。廊坊市地处中纬度地带, 属暖温带大陆性季风气候, 四季分明, 光热资源充足, 雨热同季, 有利于农作物生长, 但气象灾害较多, 干热风、雷雨冰雹大风、连阴雨、寒潮等灾害性天气常给农业生产造成不利影响。该区年平均气温 11.9℃, 年均无霜期183 d, 年平均降水量为554.9 mm, 年平均日照时数在2660 h左右。廊坊市除北部有76 km2的低山丘陵外, 98%的土地为平原, 地势开阔平坦、地层深厚,土壤类型多样, 主要作物为玉米、小麦、蔬菜、棉花等, 粮食播种面积30.62万公顷。

2 数据获取与处理

2.1 遥感数据获取及处理

使用当前农情遥感监测业务常用的 GF-1/WFV数据和 MODIS卫星遥感影像进行冬小麦长势评价分析。GF-1卫星发射于2013年4月26日, 是我国“高分专项”发射的高分系列卫星首星, 搭载了2台2 m分辨率全色/8 m分辨率多光谱相机和4台16 m分辨率多光谱相机(WFV1～WFV4)。WFV传感器4台相机组合可达 800 km幅宽, 标称重访周期 4 d。MODIS数据则是搭载在TERRA和AQUA卫星上的一个重要的传感器, 其数据通过X波段免费向全球广播, MODIS传感器最大分辨率可达250 m, 包含36个离散的波段, 光谱范围0.4～14.4 μm, 扫描宽度2330 km, 两颗卫星组合一天可过境同一地区4次(白天2次), 其中等分辨率、大幅宽、高重访周期、高光谱分辨率的优势使其在全球各项遥感监测工作中得到广泛应用, 同时也是当前农情遥感监测的重要数据源。

图1 研究区地理位置Fig. 1 Location of study area

选用2016年4月25日研究区的GF-1 WFV1影像(图 2-A), 进行几何校正、大气校正处理。几何校正使用卫星自带的 RPC参数, 并将坐标系转为UTM52, 分辨率为16 m。大气校正使用6S辐射传输模型, 能见度设置为40 km, 其他参数默认,获得GF-1地表反射率影像。MODIS数据选用2016年 4月 24日 AQUA卫星数据(图 2-B), 同样经过几何校正、大气校正处理, 并投影到UTM52坐标系上。

利用公式(1), 分别计算GF-1和MODIS影像的NDVI指数, 作为作物长势遥感监测的评价指标。为计算方便, 本文的NDVI值扩大10 000倍并取整。

图2 研究区GF-1卫星影像(A)和MODIS卫星影像(B)Fig. 2 GF-1 (A) and MODIS (B) images of the study area

2.2 冬小麦长势地面调查

2016年4月24日在河北廊坊及天津武清, 共采集 39个有效的冬小麦长势监测地面样本点(图 3),此时正值冬小麦拔节末期, 是生长速度最快、生长量最大的时期, 遥感识别特征明显。选取位于冬小麦生长密集区域的完整地块, 调查样本点的地理坐标、现场照片、冬小麦综合茎数、株高、叶绿素浓度。使用手持式GPS采集地理坐标。冬小麦综合茎数为每平方米总茎数, 使用皮尺测量结合人工计数方式获取; 用卷尺测量株高; 用SPAD-502叶绿素仪测定叶片叶绿素浓度, 以 SPAD值表现叶绿素相对含量。

3 算法描述

3.1 技术思路

使用同一时间段获取的 GF-1和 MODIS影像,并同步采集地面实测的综合茎数、株高和叶绿素浓度 3个指标, 分析 NDVI指数与冬小麦作物长势地面监测指标的关系, 并评价GF-1和MODIS影像作物长势遥感监测的精度。采用GF-1影像提取的冬小麦作物种植区域, 获取MODIS影像和GF-1影像冬小麦作物长势遥感监测结果, 从传感器光谱响应函数差异、NDVI尺度效应、地面长势状况等方面, 分析造成两类影像 NDVI长势表达不一致性的原因,比较两类影像作为长势监测数据源的优劣(图4)。

3.2 冬小麦种植范围获取

从高分辨率的 GF-1影像中提取冬小麦种植区域。首先使用最大似然分类方法, 选取冬小麦样本区域, 进行样本训练, 并得到研究区冬小麦初步分类结果, 再结合人工目视解译初步精修正分类结果,获取较为精确的研究区冬小麦种植空间分布图(图5)。利用 2016年 4月 21日的覆盖部分研究区的GF-2/PMS 4 m分辨率卫星影像, 通过监督分类结合目视解译进行冬小麦面积提取, 对GF-1获得的冬小麦种植区域结果精度验证表明, 分类总体精度为92.5%, 冬小麦的用户精度87.9%, 制图精度89.4%。

图3 冬小麦地面调查样本点位置分布图Fig. 3 Distribution of winter wheat ground sample points in the investigation area

3.3 遥感NDVI与冬小麦长势地面监测指标分析方法

根据地面调查时获取的冬小麦地面样本点经纬度坐标, 提取样本点位置处影像 NDVI值, 分别与冬小麦综合茎数、株高、叶绿素浓度地面监测指标进行线性拟合, 评价 NDVI指数与各作物长势地面监测指标的关系, 并分别评价GF-1和MODIS影像的冬小麦长势遥感监测精度, 为遥感影像作物长势监测分级提供理论依据。

3.4 GF和MODIS数据NDVI指数冬小麦长势表达差异性分析方法

由于 MODIS的地面分辨率较低(250 m), 对于冬小麦的识别能力相比GF-1卫星16 m分辨率影像较差, 同时由于混合像元的存在, 其对于冬小麦长势监测的可靠性也不如GF-1。但是, MODIS具有重访周期短(每日2次), 覆盖范围大(幅宽2330 km)等优势, 在实际业务工作中经常使用其进行作物面积获取、作物长势监测、产量评估等。为详细研究MODIS卫星影像长势监测的有效性及与GF-1影像长势监测的差异, 对比冬小麦区域 MODIS与 GF-1影像长势指标的一致性。该部分主要包括 MODIS影像重采样及与 GF-1像元对齐、渔网分割统计、MODIS冬小麦纯像元提取、MODIS与GF-1长势对比、MODIS与GF-1传感器差异及NDVI尺度效应对作物长势遥感监测的影响分析等步骤。

图4 本研究技术流程框图Fig. 4 Structure of technical flow in this study

图5 研究区冬小麦分类结果Fig. 5 Winter wheat classification in study area

3.4.1 GF-1影像重采样及与MODIS像元对齐

由于MODIS的分辨率为250 m, 与GF-1卫星16 m分辨率影像无法缜密套合, 因此将MODIS影像重采样至256 m, 这样每个MODIS像元都将包含16×16共256个GF像元, 便于进行统计分析。像元对齐主要为防止MODIS像元与GF-1影像像元之间中心点位差导致两者无法完全套合。最后获得的是范围完全一致的MODIS和GF-1卫星NDVI影像。

3.4.2 渔网分割统计 对整个研究区的 MODIS影像构建大小为256 m×256 m的渔网, 每个渔网都完整包含一个MODIS的像元和256个GF-1像元(图6), 渔网的总个数为MODIS影像的总像元数, 渔网的ID则唯一标示MODIS的像元序号。

3.4.3 MODIS冬小麦纯像元提取 由于 MODIS的分辨率较低, 存在较多的混合像元, 不利于冬小麦长势的定量分析。因此, 本文利用高分辨率的GF-1影像冬小麦分类结果(图5)来提取MODIS的冬小麦纯像元, 剔除非冬小麦及混合像元区域。使用渔网分区统计进行提取, 统计各个渔网内的GF-1影像冬小麦像元个数, 个数为 256则表示为冬小麦纯像元。

3.4.4 MODIS与GF-1长势对比分析 针对提取出来的 MODIS冬小麦纯像元, 统计其像元内的GF-1像元 NDVI的平均值及标准差, 建立 MODIS像元NDVI与GF-1像元NDVI均值之间的关系, 并统计像元内部冬小麦长势不一致性(NDVI标准差)对MODIS长势遥感监测的影响。

3.4.5 传感器光谱响应函数差异、NDVI尺度效应对作物长势遥感监测的影响 由于 MODIS和GF-1卫星传感器光谱响应函数不一致, 导致地表反射率存在一定差异, 通过统计研究区内 MODIS和GF-1卫星地表反射率的关系, 使 MODIS地表反射率转与 GF-1地表反射率一致, 并计算 NDVI; 同时考虑到 NDVI随分辨率尺度的非线性变化特性, 重新计算降分辨率后的 GF-1影像 NDVI; 最后, 对比考虑 NDVI尺度效应和传感器差异后的 MODIS和GF-1影像NDVI指数, 评价传感器差异及NDVI尺度效应对作物长势遥感监测的影响情况。

图6 GF-1 NDVI影像及按照MODIS像元大小划分渔网示意图Fig. 6 GF-1 NDVI image and the fishnet based on MODIS pixels

4 结果与分析

4.1 GF-1和MODIS的NDVI指数与作物长势地面监测指标的相关性

4.1.1 GF-1 NDVI指数与作物长势地面监测指标的相关性 针对全部 39个冬小麦长势地面实测点,根据调查时同时记录的GPS坐标信息, 在GF-1影像上获取相应位置的像元 NDVI值作为遥感监测长势指标, 与地面获取的冬小麦综合茎数、株高、叶绿素浓度对比分析。结果, GF-1影像冬小麦像元的NDVI值与综合茎数有很强的相关性(图 7-A), 而与冬小麦的株高(图 7-C)和叶绿素浓度(图 7-E)之间的线性相关性相对较弱。冬小麦长势包括物理性指标(如综合茎数和株高)和生理性指标(如叶片叶绿素浓度), 与NDVI指数相关性最好的是综合茎数。传统冬小麦长势地面监测中划分 1、2、3类苗的主要指标是单位面积的总茎数[34], 而利用 GF-1卫星影像进行冬小麦长势监测可以与传统地面长势监测相结合, 具有较高的精度和可行性。

4.1.2 MODIS NDVI指数与作物长势地面监测指标的相关性 同样使用 MODIS影像 NDVI值与上述 3个冬小麦长势监测指标进行相关分析, 发现相关度均不高(图 7-B, D, F), 说明基于 MODIS影像NDVI指数的冬小麦长势遥感监测精度比GF-1低。其原因可能是地面实地测量的冬小麦长势样本点范围较小, 而MODIS单个像元覆盖面积较大, 大量混合像元及作物长势不一致性情况的存在, 使得MODIS NDVI指数进行冬小麦作物长势遥感监测精度相对较低, 仅适合用于获取大范围内冬小麦趋势性的长势情况。

4.2 冬小麦长势一致性对MODIS数据NDVI长势表达的影响

将MODIS影像重采样至256 m分辨率, 使每个MODIS影像都包含256个完整的GF-1像元。利用GF-1影像获取的冬小麦空间分布, 剔除混合像元,只保留纯冬小麦 MODIS像元, 通过比较 MODIS NDVI与像元内GF-1 NDVI平均值之间的关系, 分析 MODIS冬小麦长势遥感监测精度与像元内冬小麦平均长势的关系。共获取3588个MODIS冬小麦像元, 统计其与像元内GF-1像元NDVI之间的关系,发现MODIS像元NDVI值与像元内部GF-1 NDVI值之间整体趋势一致(图8-A),R2为0.3944。

统计各MODIS冬小麦纯像元内的GF-1 NDVI的标准差, 并按标准差大小排序, 分别统计全部像元及标准差＜1100、＜1000、＜900、＜800、＜700、＜600、＜500、＜400、＜300、＜200 的 MODIS NDVI与像元内GF-1像元NDVI均值的相关性, 发现随着标准差的递减, MODIS与 GF-1的 NDVI相关性越来越强(图8-C),R2从0.3944增加到0.8284。据此认为, 只有当地面冬小麦长势较为均一时, MODIS的 NDVI值才能有效代表其像元内部的冬小麦平均长势。

4.3 传感器差异和 NDVI尺度效应对 GF-1和MODIS影像作物长势监测精度的影响

由于MODIS传感器和GF-1传感器的光谱响应函数之间存在差异, 导致MODIS影像和GF-1影像的反射率不一致。根据研究区MODIS和GF-1影像地表反射率, 利用线性最小二乘法拟合研究区MODIS和GF-1影像红光波段和近红外波段的转换公式(2)和(3), 并利用转换后的 MODIS地表反射率计算 NDVI指数, 减少传感器不一致造成的植被指数差异。

根据公式(1), NDVI值随分辨率尺度非线性变化, 当分辨率变化时, 不能直接用原始 NDVI值重采样获取变化后的NDVI。依据地表反射率随尺度线性变化的特点, 可用公式(4)计算分辨率变化后的NDVI值。

式中,n为分辨率变化后影像每个像元对应原始影像像元的个数, NDVIs为转换后NDVI理论值。

利用该公式重新计算 GF-1影像重采样到MODIS分辨率后的 NDVI值。将重新调整后的MODIS和GF-1影像中冬小麦像元NDVI指数对比,在重新考虑MODIS与GF-1传感器的差异, 以及不同分辨率尺度NDVI的尺度效应后, MODIS与GF-1影像对于冬小麦的长势监测结果一致性显著提高(图8-B),R2从0.3944提高到0.4633。

图7 GF-1影像(A, C, E)和MODIS影像(B, D, F)的NDVI与冬小麦长势地面监测指标对比Fig. 7 Relationship between NDVI of GF-1 image (panels A, C, and E) or MODIS image (panels B, D, and F) and ground monitoring growth index of winter wheat

利用四分位数方法对 GF-1和 MODIS影像的NDVI值进行分级划分, 分为冬小麦长势较差、长势普通、长势较好、长势好4个等级。可以看出, GF-1影像和 MODIS影像长势分级结果的空间分布趋势基本一致, 在冬小麦分布密集的研究区东北部和西部长势较好, 其他区域相对较差(图 9), 表明利用MODIS影像进行宏观尺度冬小麦长势遥感监测是可行的。

5 讨论

由于MODIS影像高重访、免费接收且发射时间较早, 基于其影像的长势监测研究和业务应用都已较为成熟[35-37], 如在中国农业科学院农业资源与农业区划研究所参与建设的中国农业遥感监测业务系统(China agricultural remote sensing monitoring system, CHARMS)[38-39]中, 即以 MODIS为主要数据源, 实现对全国主要作物信息的快速获取及分析, 定期上报我国农业监测决策部门,为保护我国粮食安全、提高农民收益等提供科学的决策依据。

图8 MODIS与GF-1像元NDVI对冬小麦长势监测精度比较Fig. 8 Comparison of accuracy on monitoring winter wheat growth between MODIS NDVI and GF-1 NDVI

图9 基于四分位数划分的GF-1 (A)和MODIS影像(B)冬小麦长势分级Fig. 9 Quartile grading of winter wheat growth in GF-1 (A) and MODIS-image (B)

然而, 本研究发现, MODIS的分辨率较低, 其长势监测结果的精度不足, 仅能进行大尺度区域的作物长势趋势性分析。相比之下, GF-1具有更高的分辨率, 获取的长势信息更加精确, 且农业部门作为该卫星的主要用户单位, 利用GF-1影像进行农作物长势监测, 可以有效提高长势监测产品的可用性,为农业部门精细化决策、作物田间管理、早期估产等提供及时可靠的信息。

值得注意的是, 作物长势的内涵丰富, 除了本文选取的 NDVI指数、综合茎数、株高、叶绿素浓度等外, 还包括干物质重量、千粒重、分蘖数、叶面积指数、EVI (enhanced vegetation index, 增强植被指数)等一系列特征[1]。由于遥感数据的宏观性,当前开展的研究主要着重于植被指数等反应宏观生长状态的指标, 对于作物植株个体、群体长势特征的反演较少[40-43], 如姜城等[44]发现 TM 影像 NDVI值在表现小麦千粒重指标上有一定的效果, 李卫国等[45]发现TM影像NDVI值与小麦植株氮素含量之间存在较高相关性。本文则初步研究了不同分辨率卫星影像NDVI指数与地面长势监测指标综合茎数、株高、叶绿素浓度之间的关系, 表明GF-1 NDVI值等综合长势特征与作物个体(如株高)、群体(如综合茎数)长势特征之间存在一定的关联性, 将有助于提高遥感长势监测结果的直观性和客观性。但由于NDVI值在植被生长旺盛时会出现饱和现象[46], 探索其他改进指数如 LAI[47]、EVI[48]等作为长势遥感监测指标反演不同作物长势, 将是今后研究的一个趋势。

6 结论

对GF-1和MODIS数据的NDVI长势监测影响最大的是冬小麦的综合茎数, 这与地面观测结果的分析一致。进行作物长势监测的卫星影像分辨率越高, 则其 NDVI值越能反映实际的作物长势。与GF-1 NDVI结果相比较, MODIS像元内冬小麦长势一致性越高, 基于NDVI的MODIS与GF数据冬小麦长势监测结果越一致。从区域长势监测角度来看,尽管MODIS与GF-1数据的监测结果趋势较为一致,并且通过光谱响应函数调整、NDVI尺度归一化后能够进一步提高监测结果的一致性, 但MODIS NDVI长势精度较低, 难以满足作物长势情况复杂地区长势监测的需求。为实现作物长势遥感监测精细化的需要, 应当逐步使用高分辨率的遥感数据替代中低分辨率遥感数据, 并重点研究利用高分辨率卫星影像进行作物长势遥感监测。