MODIS地表反射率数据的预处理方法及系统

摘要

本发明公开了一种MODIS地表反射率数据的预处理方法及系统，所述方法包括以下步骤：S1：获取原始输入的遥感数据；S2：对所述原始输入的遥感数据进行缺失数据检测；S3：对所述原始输入的遥感数据进行云、雪初步检测；S4：得到的云和雪数据作为训练样本，对所有数据检测，标识异常数据；S5：对原始输入的遥感数据进行时空滤波与插值。所述系统包括分别用于实现上述步骤的：数据输入模块、数据缺失检测模块、云雪检测模块、异常数据检测模块以及时空滤波与插值模块。本发明通过对地表反射率数据中缺失和异常数据进行处理，可以生成长时间序列、时空连续一致的地表反射率数据，提高了后续应用和遥感反演的精度。

说明书

技术领域

本发明涉及数据预处理技术领域，尤其涉及一种MODIS地表反 射率数据的预处理方法及系统。

背景技术

预处理过程是保证全球陆表卫星数据产品质量的前提。由于遥感 数据在空间、波谱和时间上受到云、雪和云影等的影响，真实的地表 反射率往往受到干扰，从而很难精确地反映出地表特征参量产品的变 化规律。

从全球每天的遥感图像上分析，全球平均70％以上地区都被云覆 盖；在时间上，云覆盖还存在长期性、季节性、多变性；同时，云覆 盖带来了云影的存在；在中高纬度区，由于大量的可溶性雪的存在， 云和雪的相似性使得图像的判别又增加了很多困难。中分辨率成像光 谱仪MODIS(Moderate-Resolution Imaging Spectroradiometer，MODIS) 数据由于其空间分辨率低、覆盖范围广，受云、云影、雪、及其他异 常气候条件的影响十分明显，从而给数据使用带来不便影响。因此云 及云影的标识和去除成为影响地表特征参量产品精度的重要组成部 分。

现有技术中，辐射传输模型广泛应用于地表反射率和地表发射率 的计算，包括短波和长波辐射传输方程。

对于短波辐射，如果忽略极化效应，辐射传输方程可以表达为：

$\frac{\mathrm{dI}\left(s\right)}{\mathrm{ds}}=-K(I-J)$

其中I(s)表示特定方向s上的辐亮度；K称为体消光系数；J称为源 函数。

假设大气是水平均一的，而且介质粒子是各向同性的，那么辐射 传输方程又可以表达为：

其中τ是光学厚度，ω是单次散射反照率，P是相函数，J0是介质 的源函数。

地表是解辐射传输方程的下边界条件。云层的反射可以看作是影 响介质各参数的因素。

对于长波辐射，辐射传输方程可以写为：

同理，云层自身的温度和辐射影响了介质的各个参数。

直观上看，太阳辐射到达地面，再经地面反射回到传感器。由于 云的存在，太阳辐射有可能被云直接反射回到传感器。而太阳辐射被 云阻挡，或穿过薄云后到达地面，在地面形成云影，而传感器对云影 成像。同时，太阳、传感器角度的变化、地形影响，都可能造成传感 器输出辐射值的变化。只有考虑这些因素，全面应用于辐射传输方程， 才能获取正确的地表反射率。

从反射光谱来看，云具有可见光和短波红外高反射率的特性，雪 具有在可见光区域高反射、在短波红外高吸收的特性，如图3。从发 射光谱看，云和雪的温度较低。

同时，云具有高度的流动性和变化性，而地物往往保持了一定的 连续性。因此，识别云雪和地物既可以从其光谱属性，也可以从其温 度属性出发，也可以从时空变化属性出发。

从反射光谱来看，云具有可见光和短波红外高反射率的特性，雪 具有在可见光区域高反射、在短波红外高吸收的特性。同时，云具有 高度的流动性和变化性，而地物往往保持了一定的连续性。因此，识 别云雪和地物既可以从其光谱属性，也可以从其温度属性出发，或者 可以从时空变化属性出发。利用更多的信息，进一步改善云、雪、云 影的标识和去除成为影响地表特征参量产品精度的重要组成部分。

目前，对于云的识别主要分为两块：基于单个区域(Swath)的 识别和基于时间序列的识别。

目前比较成熟的单Swath云掩模算法的特点是基于单一像元，根 据云的强反射和低温度特性，采用反射光谱和亮温的阈值来判断该像 元是否被云污染。

MODIS云掩模目前采用的算法基本可分为如下几类：

1)基于热红外波段亮温(尤其以11μm为主)：单通道亮温(BT)、 通道间亮温差(BTD)；

2)基于反射波段反射率：单通道反射率、通道间反射率比值；

3)基于水汽吸收或通透波段(以1.38μm、11μm为主)；

4)CO2切片法；

5)时空不一致性检测。

这些算法的原理来自于可见光、红外的大气辐射传输模型和地物 反射、发射特性。检测时的阈值设置最终来源于时间、气象、大气辐 射传输模型、地物生态特性。

为了判定合适的阈值，一般有很多先验知识输入，MODIS云掩 模采用的先验知识有：

1)白天/黑夜；

2)太阳耀斑：计算太阳散射角；

3)雪/冰：NSDI，或每日雪/冰图；

4)每日海洋冰产品；

5)水/陆图(1km)：水、陆、沙漠、海岸；

6)太阳高度角、方位角和视角；

7)地形；

8)生态系统；

9)其它尽可能获取的知识。

目前，由于已有多样化的其它卫星产品，更多的先验知识可以用 于检测云，如：

1)气象学、臭氧、气溶胶数据集：温度、湿度、风、臭氧、气 溶胶的垂直分布

2)地表发射率图

3)地表类型图

4)晴空反照率图

在最近几年，关于以上方法和阈值改进的文章也有不少。但是， 由于这些方法在时间上和同一地类内采用了固定阈值，而其阈值来源 于经验，因此当其应用于全球范围内时，在地物和云皆千变万化的情 况下，容易产生各种误判和漏判。即使后来的阈值和应用地区进一步 细化和改进，但是由于本质没有改变，仍然存在同样的问题。

CAPPELLUTI于2006年提出了局部云的自动判别方法。

而刘荣高等于2007年以来提出了很多新的云检测方法，包括动态 阈值法，概率表达法、针对低云的波段3、9、31检测法。

这些新方法很大程度上改进和扩展了单Swath的云检测。

另一方面，迄今为止，像MODIS、AVHRR之类的极轨卫星，已 经形成了多年的每天时间序列数据，这使得利用时间序列信息判断云 也逐渐成为一个新的手段。在高频率的观测下，相对于变化较快的云， 陆地表面可以被看作是静态或缓变的背景。

最初，Rossow，W.B等提出的ISCCP云掩模算法根据以前测量的 数据建立清晰天空合成图，对比当前测量值和清晰天空参考值判断每 个像元是否为云。而参考值的不确定度由地物的变化和传感器的噪声 确定，根据时间序列数据计算。

A.Lyapustin等根据这个思路改进了方法，提出了利用单波段空间 图像的协方差变化分析识别清晰天空图像和多云图像的方法。

O.Hagolle等同样也沿用了这个思路，提出了某像元蓝光通道， 不同时间反射率和参考清晰天空反射率之差关于相应天数间隔dD的 模型，认为反射率之差大于0.03*(1+dD/30)时为云。

关于云影方面，一般方法都以太阳高度角和云位置的几何投影关 系为主，结合少量的波段阈值判断。

其余关于这两类方法的文献还有很多，上述所有这些方法在地物 的光学属性和空间的时间变化统计特征上做出了很多成果，但是没有 使用光谱的时间变化统计特征。

在冰雪的识别方面，由于冰雪在可见光和近红外高反射，在短波 红外高吸收的特征，使得NDSI(归一化冰雪指数)成为识别云和冰 雪的一个指标。但是当冰雪和土壤混合时，或者卷云时，云和冰雪的 区分出现了困难。

MODIS和NASA团队采用了NDSI以及几个波段的亮温判别指标 来判别冰雪。

尽管云在可见光和近红外表现为高反射性，但是在实际处理时， 存在以下几个难点：

云的多样性、混合像元的存在、前端处理的干扰(时间序列)、 云和冰雪等的区分、长期云覆盖下的数据填充、考虑BRDF下的数据 填充。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：提供一种MODIS地表反射率数据 的预处理方法及系统，以生成长时间序列、时空连续一致的地表反射 率数据，为后续数据处理提供方便。

(二)技术方案

为解决上述问题，一方面，本发明提供了一种MODIS地表反射 率数据的预处理方法，包括以下步骤：

S1：通过MODIS星载传感器，获取原始输入的遥感数据；

S2：对所述原始输入的遥感数据进行异常数据检测；

S3：对所述原始输入的遥感数据，基于云、雪先验知识进行云、 雪初步检测；

S4：得到的云和雪数据作为训练样本，利用所述训练样本，对所 有数据检测，标识出所有云、雪数据为异常数据；

S5：对原始输入的遥感数据进行时空滤波与插值，以填充在长时 间序列中缺失的像元和空间上异常的像元。

优选地，所述步骤S2中进行异常数据检查的包括：

填充像元检查：标识提取所述原始输入的遥感数据中存在的填充 像元，用于后续插值和滤波；

错误像元检查：识别所述原始输入的遥感数据中存在的与真实数 据相差很远的数据以及饱和数据，用于后续插值和滤波；

缺失像元数据检查：利用数据的时间标识检查所述原始输入的遥 感数据中存在的时间缺失数据，用于后续插值和滤波。

优选地，所述错误像元的检查方法包括：

先验知识：地表反射率数据应该介于0和1之间，凡是低于0或 高于1的值均认为是错误像元；

长时间序列像元相关系数：以一段时间的数据为处理单元，同一 像元在不同的时间会形成一个时间序列数据；根据两个像元时间序列 数据可以计算出所述长时间序列像元相关系数，公式如下：

$r_{\mathrm{xtyt}}=\frac{{(\mathrm{xt}-{\mu }_{\mathrm{xt}})}^T(\mathrm{yt}-{\mu }_{\mathrm{yt}})}{\left|\right|\mathrm{xt}-{\mu }_{\mathrm{xt}}\left|\right|\left|\right|\mathrm{yt}-{\mu }_{\mathrm{yt}}\left|\right|}$

xt和yt分别为所述原始输入的遥感数据同一波段在不同时间上 的反射光谱值，μ_xt和μ_yt分别为xt和yt在时间上的均值；

其中，长时间序列像元相关系数的绝对值R位于0和1之间，R 越大，误差越小，变量之间的线性相关程度越高。

优选地，所述步骤S3中进行云、雪初步检测的方法为：

S31：根据所述原始输入的遥感数据中相邻像元在不同波段的反 射光谱值，计算所述相邻像元在不同波段反射谱值的光谱角度和光谱 相关系数；

S33：通过所述光谱角度和光谱相关系数对所有时间点数据进行 云、雪检测。

优选地，在所述步骤S31和S33之间还包括步骤S32：设置相关 系数阈值，通过所述相关系数阈值对所述原始输入的遥感数据进行筛 选，去除异常数据，留下正常数据。

优选地，所述光谱角度和光谱相关系数的计算公式为：

光谱角度θ：$\theta =\mathrm{ar}\cos \frac{x^{T}y}{\left|\right|x\left|\right|\left|\right|y\left|\right|}$

像元光谱相关系数r_xy：$r_{\mathrm{xy}}=\frac{{(x-{\mu }_x)}^T(y-{\mu }_y)}{\left|\right|x-{\mu }_x\left|\right|\left|\right|y-{\mu }_y\left|\right|}$

其中，x和y分别为所述原始输入的遥感数据相邻像元的不同波 段反射光谱值，μ_x和μ_y分别为x和y在波段上的均值。

优选地，步骤S4具体包括以下步骤：

S41：计算数据的归一化雪被指数NDSI：NDSI＝(R4-R3)/(R4+R3)

其中R4为波长为0.555微米的地表反射率，R3为波长为1.64 微米的地表反射率；

S42：根据计算的数据的NDSI数值对进行判断：

如果NDSI＞0.5，且地理位置和时间符合下雪条件的，则该数据 识别为纯雪；

如果NDSI＜0.4，则该数据识别为云；

S43：将已识别为云和雪的数据分别作为云训练样本和雪训练样 本，根据训练结果，利用最大似然法对0.4＜NDSI＜0.5的异常数据进 行云和雪的分类。

优选地，所述步骤S5中的插值方法为：利用一年时间序列内的 数据，根据同一类地物光谱在时间和空间上的连续性和相关性特性， 采用多项式拟合的方法进行填充插值。

另一方面，本发明还提供了一种MODIS地表反射率数据的预处 理系统，所述系统包括：

数据输入模块，用于通过MODIS星载传感器，获取原始输入的 遥感数据；

数据缺失检测模块，用于对所述原始输入的遥感数据进行缺失数 据检测；

云雪检测模块，用于对所述原始输入的遥感数据，基于云、雪先 验知识进行云、雪初步检测；

异常数据检测模块，用于得到的云和雪数据作为训练样本，并利 用所述训练样本，对所有数据检测，标识出所有云、雪数据为异常数 据；

时空滤波与插值模块，用于对原始输入的遥感数据进行时空滤波 与插值，以填充在长时间序列中缺失的像元和空间上异常的像元。

(三)有益效果

本发明利用时空序列数据的物理特性，判别其误判和漏判的数 据，并加以滤波和缺失填充，最终去除云、云影等的影响，生成长时 间序列、时空连续一致的地表反射率数据，能够有效提高这些地表特 征参量数据的精度；使得经过本发明预处理后的数据能够直接服务于 GLASS发射率、反照率、叶面积指数和辐射产品的生产。

附图说明

图1为根据本发明实施例预处理方法的步骤流程图；

图2为根据本发明实施例预处理方法步骤3的具体步骤流程图；

图3为根据本发明实施例预处理方法步骤4的具体步骤流程图。

图4为2003年MODIS地表反射率数据h24v05上某一像元预处理 前后时间序列数据图；

图5为2003年MODIS地表反射率数据h20v03上某一像元预处理 前后时间序列数据图。

具体实施方式

下面结合附图及实施例对本发明进行详细说明如下。

实施例一：

如图1所示，一种MODIS地表反射率数据的预处理方法，包括 以下步骤：

S1：通过MODIS星载传感器，获取原始输入的遥感数据；

在本实施例中，所述原始输入的遥感数据为2000年至2010年中某 一年的MODIS地表反射率数据，可以是下面4类数据中的一种：

(1)MOD09A1数据集；

(2)MYD09A1数据集；

(3)MOD09GA数据集；

(4)MYD09GA数据集。

其中，MOD09A1和MYD09A1数据是MODIS 8天地表反射率； MOD09GA和MYD09GA是MODIS 1天地表反射率。这些数据从2000 年开始收集，数据中含有MODIS的1-7波段500米分辨率的反射率数据 和1km分辨率的观测天顶角、观测方位角、太阳天顶角、太阳方位角 和质量控制等信息。

S2：对所述原始输入的遥感数据进行异常数据检测；

具体包括以下几种方法：

填充像元检查：由于MODIS09产品几乎每个Tile都存在一些填 充数据，因此需要标识提取所述原始输入的遥感数据中存在的填充像 元，用于后续插值和滤波；

错误像元检查：由于探测器错误或饱和、电子线路问题等，在 MODIS09产品中存在一些与真实数据相差很远的数据，比如冰雪的 高端饱和与水面的低端饱和，因此需要识别所述原始输入的遥感数据 中存在的与真实数据相差很远的数据以及饱和数据，用于后续插值和 滤波；这些异常数据特征往往是某个波段出现异常，如负值或毛刺， 因此识别这些错误像元的识别方法包括：

先验知识：地表反射率数据应该介于0和1之间，凡是低于0或 高于1的值均认为是错误像元；

长时间序列像元相关系数：以一段时间的数据为处理单元，同一 像元在不同的时间会形成一个时间序列数据；根据两个像元时间序列 数据可以计算出所述长时间序列像元相关系数r_xtyt，公式如下：

$r_{\mathrm{xtyt}}=\frac{{(\mathrm{xt}-{\mu }_{\mathrm{xt}})}^T(\mathrm{yt}-{\mu }_{\mathrm{yt}})}{\left|\right|\mathrm{xt}-{\mu }_{\mathrm{xt}}\left|\right|\left|\right|\mathrm{yt}-{\mu }_{\mathrm{yt}}\left|\right|}$

xt和yt分别为所述原始输入的遥感数据同一波段在不同时间上 的反射光谱值，μ_xt和μ_yt分别为xt和yt在时间上的均值；

其中，长时间序列像元相关系数的绝对值R位于0和1之间，R 越大，误差越小，变量之间的线性相关程度越高。

时间缺失数据检查：遥感数据应该是时间上连续的，例如每年的 MODIS09A1产品都是应该从001编号开始，361编号结束，每隔8 天一个产品，但是也会存在时间上不连续、丢失或无产品的情况，因 此需要利用数据的时间标识检查所述原始输入的遥感数据中存在的 时间缺失数据，用于后续插值和滤波。

S3：对所述原始输入的遥感数据，基于云、雪先验知识进行云、 雪初步检测；

如图2所示，所述步骤3具体包括以下步骤：

S31：根据所述原始输入的遥感数据中相邻像元在不同波段的反 射光谱值，计算所述相邻像元在不同波段反射谱值的光谱角度和光谱 相关系数；

S32：设置相关系数阈值，通过所述相关系数阈值对所述原始输 入的遥感数据进行筛选，去除异常数据，留下正常数据；

S33：通过所述光谱角度和光谱相关系数对所有时间点数据进行 云、雪检测。

所述光谱角度和光谱相关系数的计算公式为：

光谱角度θ：$\theta =\mathrm{ar}\cos \frac{x^{T}y}{\left|\right|x\left|\right|\left|\right|y\left|\right|}$

像元光谱相关系数r_xy：$r_{\mathrm{xy}}=\frac{{(x-{\mu }_x)}^T(y-{\mu }_y)}{\left|\right|x-{\mu }_x\left|\right|\left|\right|y-{\mu }_y\left|\right|}$

其中，x和y分别为所述原始输入的遥感数据相邻像元在不同波 段反射光谱值，μ_x和μ_y分别为x和y在波段上的均值。

S4：识别云和雪两类的数据，对得到的云和雪数据作为训练样本， 利用所述训练样本，对所有异常数据分类；

如图3所示，所述步骤S4具体包括以下步骤：

S41：计算数据的归一化雪被指数NDSI：NDSI＝(R4-R3)/(R4+R3)

其中R4为波长为0.555微米的地表反射率，R3为波长为1.64 微米的地表反射率；

S42：根据计算的数据的NDSI数值对进行判断：

如果NDSI＞0.5，且地理位置和时间符合下雪条件的，则该数据识 别为纯雪；

如果NDSI＜0.4，则该数据识别为云；

S43：将已识别为云和雪的数据分别作为云训练样本和雪训练样 本，根据训练结果，利用最大似然法对0.4＜NDSI＜0.5的异常数据进 行云和雪的分类。

S5：对原始输入的遥感数据进行时空滤波与插值，以填充在长时 间序列中缺失的像元和空间上异常的像元。其中的插值方法具体为： 利用一年时间序列内的数据，根据同一类地物光谱在时间和空间上的 连续性和相关性特性，采用多项式拟合的方法进行填充插值。

图4和图5分别为两段实际的数据中某一像元在预处理前后时间 序列数据图，可以看出经过本发明的预处理之后，有效地去除了云雪 的影响，使数据更符合地表实际状况。

实施例二：

本实施例记载了一种MODIS地表反射率数据的预处理系统，所 述系统包括：

数据输入模块，用于通过MODIS星载传感器，获取原始输入的 遥感数据；

数据缺失检测模块，用于对所述原始输入的遥感数据进行缺失数 据检测；

云雪检测模块，用于对所述原始输入的遥感数据，基于云、雪先 验知识进行云、雪初步检测；

异常数据检测模块，用于得到的云和雪数据作为训练样本，并利 用所述训练样本，对所有数据检测，标识出所有云、雪数据为异常数 据；

时空滤波与插值模块，用于对原始输入的遥感数据进行时空滤波 与插值，以填充在长时间序列中缺失的像元和空间上异常的像元。

本发明对遥感数据进行预处理，从而生成长时间序列、时空连续 一致的地表反射率数据，能够有效提高这些地表特征参量数据的精 度；使得经过本发明预处理后的数据能够直接服务于GLASS发射率、 反照率、叶面积指数和辐射产品的生产。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关 技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下， 还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明 的范畴，本发明的专利保护范围应由权利要求限定。