作品编号:D737(三等奖)
作品名称:面向江汉平原农情监测的多源遥感数据时空融合应用
作品单位:中国地质大学(武汉)李四光学院
小组成员:杜恩宇,赵赛赛,郑铠沅,周泉
指导老师:沈永林,沈国铃
作品简介
一、作品简介
由于技术和经费的限制,很难获取同时具有高空间分辨率和高时间分辨率的遥感数据。而时空融合算法可以融合高空间分辨率遥感数据和高时间分辨率遥感数据,从而生成具有高时空分辨率的遥感数据,对遥感数据的应用具有重要意义。本项目在已有的STARFM(Spatial and Temporal Adaptive Reflection Fusion Model)和ESTARFM(Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model)两种时空融合模型基础上,提出一种基于高斯函数的时空融合算法,用IDL语言编程实现三种模型并进行对比,证明提出的高斯时空融合算法用于时空融合的可行性,并以Landsat 8和MODIS MOD09GQ为主要数据源,以江汉平原的早稻、中晚稻、棉花为研究对象,结合农作物的物候规律,获取高空间分辨率时间序列的NDVI图像,最终应用于江汉平原的农作物分类及农情监测。
二、应用背景
农业遥感是遥感应用领域的一个重要应用方向,长时间连续的高空间分辨率观测数据对于作物生长状况监测十分重要,然而由于多种因素的限制,使得监测具有很大程度上的“时空数据缺失”,即缺少同时具有高空间和高时间分辨率的遥感数据。
遥感数据时空融合技术是一种融合中分辨率影像(Landsat)的高空间分辨率特征和低空间分辨率影像(MODIS)的高时间分辨率特征从而生成高时空分辨率遥感数据的方法。近年来,国内外学者提出了几种遥感影像时空融合方法。STARFM和ESTARFM是目前已有的应用较多的时空融合算法,在此基础上,为满足田块破碎、地形复杂区特定作物农情监测的应用需求,本项目提出了一种基于高斯函数模型的时空融合算法:根据NDVI指数在作物生长期内的变化规律,选取高斯函数对时空融合模型进行表征。通过中分辨率NDVI和低分辨率NDVI的函数映射关系构建时空融合模型,重建密集NDVI时序影像。
三、应用目标
1、运用高斯时空融合算法,融合Landsat8与MODIS MOD09GQ遥感数据,“预测”出同时具有高时间和高空间分辨率的NDVI数据。
2、利用融合后的30m时序NDVI影像的空间特征、时间特征以及农作物的物候特征对江汉平原的农作物进行马氏距离分类和光谱角分类(Spectral Angle Mapper SAM)两种分类方法,并进行种植面积估算并验证,证明高斯时空融合算法可提高分类精度。
3、应用融合后的30m时序NDVI影像,对江汉平原进行农情监测。
四、主要技术流程
图1 主要技术流程图
1.数据获取
(1)下载江汉平原2017年4-11月可用的Landsat8和MODIS MOD09GQ遥感数据。
(2)从指导老师处获得用于监督分类的江汉平原早稻、晚稻和棉花三种作物类别样本。
(3)下载USGS发布的全球农田30米高分辨率互动地图,利用ARCGIS裁剪出江汉平原地区范围内的农田分布数据,其中黄色为农田。
图2 江汉平原矢量 | 图3 USGS农田30m互动地图江汉平原内类别 |
2. 数据处理
(1)利用ENVI Desktop对Landsat 8影像进行辐射定标、大气校正、重投影、裁剪拼接和计算NDVI处理。
(2)利用MODIS Reprojection Tool(MRT)工具以及IDL实现MODIS影像数据的批处理,包括重投影、裁剪拼接、计算NDVI及重采样。
(3)分类处理
利用ENVI Desktop完成Landsat 8数据的分类和精度评价。ROI选取物候变化区域,采用最大似然或马氏距离监督分类,计算混淆矩阵。分类结果作为是否为物候变化区的判断准则应用于高斯时空融合算法中。
3.算法实现、优化和验证
(1)使用IDL开发语言编程实现STARFM、ESTARFM和高斯时空融合算法,并使用GPU进行算法加速,以减小时间复杂度。下图为算法界面:
图4 STARFM算法界面 | 图5 ESTARFM算法界面 | 图6 高斯算法界面 |
(2)利用模拟数据,以均方根误差为指标对ESTARFM和高斯时空融合算法进行精度评价并对比。优化高斯时空融合算法,包括参数设置、窗口大小调整和减少算法时间复杂度等。模拟出的Landsat8和MODIS的两条高斯函数,共分图7、8、9、10四种情况。图12、13、14为两高斯函数均值方差都不同时,其中3种输入情况中,所有输出影像的平均RMSE。
图7 同均值同方差 | 图8 同均值不同方差 | 图9 不同均值同方差 | |||
图10 不同均值不同方差 | 图11 ESTARFM与高斯优化算法各区域平均RMSE | ||||
图12 输入图像时间位于高斯函数对称轴左 | 图13输入图像时间位于高斯函数对称轴右 | 图14 第1幅输入图像时间位于高斯函数左,第2幅右 |
(3)在预处理后的真实数据上选取样本,拟合NDVI变化曲线,进行算法实验,通过比较均方根误差对ESTARFM和高斯时空融合算法进行精度评价并对比,如图11。
4.应用算法进行时空融合
合理安排数据的输入顺序及其对应的预测时间范围,融合出214张30m的NDVI图像,将这些图像合为一张,得到30m空间分辨率的NDVI时序图像。图15为融合结果。
图15 融合结果 | 图16 各作物时序NDVI曲线 |
5.江汉平原应用
(1)分类
用Savitzky-Golay滤波对融合后的30m空间分辨率NDVI时序图像做平滑处理,再利用中晚稻样方,分别进行马氏距离监督分类和光谱角分类。如图16,光谱角分类可充分利用各农作物时序NDVI曲线的差异。接着利用美国的30m分辨率的农田影像对分类结果进行分类后处理,只有在美国30m分辨率农田影像中属于农田范围的最终判定为农田。最后计算混淆矩阵评价分类精度并提取中晚稻种植面积,将两个分类结果和用大气校正后的Landsat8数据分类结果共三种进行综合对比,再与国家统计局网站获取的中晚稻种植面积进行对比,分析精度。表1是真实数据及各分类方法得到的种植面积和像元数;图17、18、19是中晚稻的三种分类结果图。
表1 中晚稻分类结果对比表
真实数据 | 单幅Landsat影像马氏距离分类 | NDVI时序影像马氏距离分类 | NDVI时序影像光谱角分类 | |
像元数(个) | 9121638 | 7984975 | 8851928 | 8988151 |
种植面积(千公顷) | 820 | 718 | 796 | 808 |
图17 中晚稻单幅Landsat影像 马氏距离分类 | 图18 中晚稻NDVI时序影像马氏距离分类 | 图19 中晚稻NDVI时序影像光谱角分类 |
由上表可知,对单幅Landsat影像进行马氏距离分类的结果精度最低;对融合后的NDVI时序影像进行分类,可以提高分类精度;基于融合后获取的NDVI时序曲线进行的光谱角分类精度则最高。采用本算法融合结果进行分类,可显著提高分类精度,因此本算法可广泛用于农作物分类相关的应用。
(2)农情监测
从时序NDVI图像中选择3个波段,做NDVI的分级,分析长势监测农情。如下图,颜色越深,NDVI值越大,农作物生长越旺盛。4月1日,左上角及左下角颜色较深,该处农作物最早进入旺盛生长;7月9日,整个江汉平原颜色较深,正值盛夏,各农作物均生长旺盛;10月27日,图中绿色较上一幅变浅,农作物生长变得缓慢,如:水稻已变黄结穗。
图20 2017年4月1日江汉平原农作物长势 | 图21 2017年7月9日江汉平原农作物长势 | 图22 2017年10月27日江汉平原农作物长势 |
五、关键技术
1.高斯时空融合算法的实现
本项目最重要的环节就是高斯优化算法的实现,在实现已有ESTARFM算法的基础上,有针对性的提出并实现基于高斯函数模型的时空融合算法,通过进行模拟数据实验和真实数据实验验证分析算法,最终证明高斯时空融合算法用于农情监测的可行性。
该函数模型分为线性部分和非线性部分。其中,线性部分主要反映因不同传感器的轨道参数、波段宽度、获取时间、光谱响应函数等差异,造成的对应 NDVI 的不一致;非线性部分反映作物物候变化,用高斯函数表征。具体思路为:首先获取研究区范围内低分辨率(MODIS)和中分辨率(Landsat 8)遥感影像,从中挑选出时刻对应、数据质量良好、无云覆盖的影像对;进行数据预处理;融合时,用物候不变化像元求解其线性部分映射关系;用物候变化像元,求解中分辨率NDVI和低分辨率NDVI的非线性映射关系,即高斯函数的映射关系,根据该关系构建时空融合模型;重建密集 NDVI 时序影像。
图23 算法思路图
例如,
时刻中分辨率的NDVI值
和其对应重采样后的低分辨NDVI值
,可表示为:
式中,
表示每个中空间分辨率像元在低空间分辨率像元中的位置;
为线性部分的参数;
为非线性部分的函数。
2.算法的加速
由于算法的时间复杂度很高,经计算正常得到一景Landsat8影像融合结果需10天左右,故本项目使用CUDA进行GPU硬件加速,可将分类时间缩短到40分钟,达到时间上的优化。CUDA(Compute Unified Device Architecture)是一种通用并行计算架构,CUDA在架构上采用了一种全新的计算体系结构来使用GPU提供的硬件资源,从而给大规模的数据计算应用提供了一种比CPU更加强大的计算能力。
3.MODIS数据的批处理
本项目在数据源上选用了250米MODIS MOD09GQ数据,由于其较高的时间分辨率,若在研究区域和时间内对其按顺序进行预处理会耗费大量时间,故使用MRT工具进行MODIS数据的批处理。MRT ( MODIS Reprojection Tool,MODIS重投影工具)被开发用于支持MODIS陆地高级产品,这些MODIS产品为HDF-EOS栅格文件,基于分幅的Sinusoidal(一种投影方式)投影子。MRT通过提供地图投影、格式转换和光谱与空间子集选项等功能方便了MODIS陆地产品的应用。
4.利用S-G滤波器进行滤波
将融合得到的时序NDVI影像进行Savitzky-Golay滤波。S-G滤波器被广泛地运用于数据流平滑除噪,是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变。
5.使用马氏距离和光谱角匹配对时序NDVI影像分类
30m时序NDVI图像包含了丰富的空间和时间信息,可提高单时相Landsat8的分类精度。故在对S-G滤波后的时序NDVI影像进行分类操作时,采用两种分类方法,分别是马氏距离监督分类和光谱角匹配分类。
光谱角分类方法是一种光谱的匹配技术,这种技术基于估计像元光谱与样本光谱或是混合像元中亚像元组分(endmember)光谱的相似性来区分各像元点的光谱曲线,以运算影像像元光谱与样本参考光谱之间的夹角来进行分类。