CN102298150B - 全球陆表宽波段发射率反演方法及系统 - Google Patents

Abstract

本发明涉及卫星遥感技术领域，特别涉及一种全球陆表宽波段发射率反演方法及系统，该方法包括：S1：获取多种分辨率的地表反射率和反照率，并获取配套的土壤分类图数据；S2：对多种分辨率的地表反射率和反照率分别进行预处理；S3：对土壤分类图数据进行空间重采样，对经预处理后的多种分辨率的地表反射率和反照率进行空间匹配，并将经空间匹配后的地表反照率逐个像元进行标识；S4：判断陆地像元的地表类型；S5：建立窄波段发射率与宽波段发射率之间的换算关系，得到水体像元和冰雪像元宽波段发射率；S6：计算出陆地像元的宽波段发射率。本发明通过对反射率数据和反照率数据的处理，提高了宽波段发射率反演的精度。

Description

全球陆表宽波段发射率反演方法及系统

技术领域

本发明涉及卫星遥感技术领域，特别涉及一种全球陆表宽波段发射率反演方法及系统。

背景技术

地表发射率为地表固有特性，受地表的组成成分、含水量和粗糙度等影响，能够刻画地表热辐射能力的强弱。地表发射率对岩石中的矿物成分具有指示性，常用于岩性识别和地质填图等领域。地表发射率除了对地表温度估算有着重要的影响，它和地表温度决定了地表辐射能量平衡中的长波辐射，是气候、水文、生态和生物地球化学模式中的关键输入参数。

遥感是获得区域、全球尺度上陆表发射率的唯一手段。然而，地表温度和发射率耦合在一起，由传感器的辐射测量值进行反演属于病态问题，即使用N个观测值求解N+1个未知数，必须采用一些策略构造多余观测，使方程组完备。很多遥感科学家致力于解决这类问题，提出了很多算法，例如参考通道、发射率归一化、独立于温度的光谱指数、日夜算法、光谱比、alpha发射率、灰体发射率，TES方法、光谱迭代平滑、多像元、优化、多层感知器网络、相关性算法和逐步求精方法。虽然有这么多算法可以使用，但真正用于业务化反演发射率的却很少。究其原因主要是，地表类型复杂多变，对算法的普适性提出了严峻的挑战。用于业务化反演发射率方法主要有ASTER的TES算法、中分辨率成像光谱仪(Moderate-Resolution ImagingSpectroradiometer，MODIS)的日夜算法和归一化植被指数(Normalized Difference Vegetation Index，NDVI)阈值法。

此外，一些学者在反演大尺度陆表发射率方面进行有益的尝试，并取得了显著进展，例如：Wilber等将地表划分为18个类型，根据实验室反演的发射率光谱，获得全球10×10格网的12个窄波段(＞4.5μm)和宽波段发射率(5-100μm)，并将它们用于辐射传输模式和NASA云与地球辐射能量系统(CERES)。Seemann等提出了基线拟合方法(Baseline fit method)，由MODIS发射率反演产品(MODll)测量全球陆表0.05°空间分辨率、10个波段的发射率(3.6-14.3μm)，用于改善晴空时大气温湿度廓线回归反演的精度。Pequignot等提出了多光谱方法(Multispectral Method，MSM)，由三年的AIRS数据生成了全球南北纬30°之间、空间分辨率为1°×1°的月平均光谱发射率(3.7-14.0μm，0.05μm光谱分辨率)。这些陆表发射率计算方面的尝试同样存在一定缺陷，Wilber等的工作属于典型的根据地表类型赋予发射率，无法反映地表发射率的真实变化。Seemann等的工作使用单一传感器，无法生成长序列的全球陆表发射率。Pequignot等的工作和Seemann等的工作类似，且空间分辨率太粗。

当前的气候模式里面，对陆表发射率的处理很粗糙，经常当作常数处理。例如，在静止环境卫星(GOES)业务化的数据处理算法里面，地表被看作发射率为0.96的灰体；欧洲的中长期天气预报(ECMWF)模式在进行干旱监测时假设所有的地表类型的发射率为常数；NCAR的通用陆面模式(CLM2)里面由叶面积指数(LAI)计算冠层的发射率，分别将土壤和雪的发射率设置为0.96和0.97。具有一定时空分辨率、长时间序列的全球陆表宽波段发射率将会显著改善气候模式的精度，目前这样的数据集仍不具备。

现有的发射率反演方法存在明显的缺陷。例如，ASTER的TES算法对于农田区域的发射率反演误差很大；MODIS的昼夜算法得到的发射率时序变化不明显，未能表征植被覆盖对发射率的影响；NDVI阈值法无法反演裸土发射率的巨大变化。

土壤光谱发射率变化大，其发射率确定比较困难。专门针对土壤，生成1公里或5公里空间分辨率发射率的算法，未见文献报道。Sobrino等人基于NDVI阈值法将地表区分土壤，部分植被覆盖和完全植被覆盖，分别获得其发射率，并给出了针对AVHRR、MODIS、SEVIRS、AASTR和TM等数据的数学表达式。Sobrino等人基于NDVI阈值法将地表区分土壤，部分植被覆盖和完全植被覆盖，分别获得其发射率，并给出了针对AVHRR、MODIS、SEVIRS、AASTR和TM等数据的数学表达式。事实上地表土壤发射率变化范围约在0.65-1.0之间，使用NDVI阈值法在全球尺度上用一个值代替或者用红光波段的反射率来表达，显然不合理。

ASTER的TES算法对于农田区域的发射率反演误差很大；MODIS的昼夜算法得到的发射率时序变化不明显，未能表征植被覆盖对发射率的影响；NDVI阈值法无法反演裸土发射率的巨大变化。

基于植被指数估算地表发射率是遥感反演发射率的常见方法，Sobrino等人基于NDVI阈值法来区分植被与非植被区，通过植被覆盖度获得像元的发射率，并给出了针对AVHRR、MODIS、SEVIRS、AASTR和TM等数据的数学表达式。然而，他们的结果多来自于对地面测量数据或者有限遥感数据的统计，并且难以处理高植被区和低植被区的情况，发射率产品的空间连续性差，将其运用于生成全球地表发射率产品必然存在较大的误差。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：如何提供一种适用于全球陆表大尺度范围的宽波段发射率反演算法，以提高宽波段发射率反演的精度。

(二)技术方案

为解决上述技术问题，本发明提供了一种全球陆表宽波段发射率反演方法，包括以下步骤：

S1：通过MODIS或AVHRR卫星遥感传感器获取多种分辨率的地表反射率和反照率，并根据所述地表反射率和反照率获取多种分辨率的以及配套的ASTER宽波段发射率数据、ASTER波谱库数据、MODISUCSB发射率波谱库数据和土壤分类图数据；

S2：对所述多种分辨率的地表反射率和反照率分别进行预处理；

S3：对所述土壤分类图数据进行空间重采样，对经预处理后的多种分辨率的地表反射率和反照率进行空间匹配，并将进行空间匹配后的地表反照率逐个像元标示出水体像元、冰雪像元和陆地像元，通过进行空间匹配后的地表反射率的红波段和近红外波段对所标示的陆地像元计算归一化植被指数；

S4：通过所述归一化植被指数判断出所述陆地像元的地表类型，所述陆地像元的地表类型分为：土壤类型、植被类型和过渡区类型，并根据经空间重采样后的土壤分类图对所述陆地像元进行土壤类型的划分；

S5：通过所述ASTER波谱库数据、以及MODIS UCSB发射率波谱库数据，建立窄波段发射率与宽波段发射率之间的换算关系，得到水体像元和冰雪像元宽波段发射率；

S6：基于ASTER宽波段发射率数据，分别计算出陆地上土壤类型、植被类型和过渡区类型三个类型像元的宽波段发射率。

其中，步骤S2具体包括：

S21：对所述多种分辨率的地表反射率和反照率进行判断，初步筛选出正常的像元和异常的像元；

S22：从正常的像元中选取已被卫星遥感传感器标识好的像元作为训练样本，分别计算所述异常的像元和正常的像元与训练样本之间的相关系数或相似系数，并判断计算出的相关系数或相似系数是否大于或等于相关系数阈值或相似系数阈值，若是，则判定为正常的像元，否则判定为异常的像元；

S23：基于地理位置、时间和归一化雪被指数判断云和雪；

S24：通过时空插值滤波的方式来填充在长时间序列中缺失的像元和空间上异常的像元。

其中，所述相关系数阈值为0.9，相似系数阈值为0.98。

本发明还公开了一种全球陆表宽波段发射率反演系统，包括：

获取模块，用于通过MODIS或AVHRR卫星遥感传感器获取多种分辨率的地表反射率和反照率，并根据所述地表反射率和反照率获取多种分辨率的以及配套的ASTER宽波段发射率数据、ASTER波谱库数据、MODIS UCSB发射率波谱库数据和土壤分类图数据；

预处理模块，用于对所述多种分辨率的地表反射率和反照率分别进行预处理；

匹配标示模块，用于对所述土壤分类图数据进行空间重采样，对经预处理后的多种分辨率的地表反射率和反照率进行空间匹配，并将进行空间匹配后的地表反照率逐个像元标示出水体像元、冰雪像元和陆地像元，通过进行空间匹配后的地表反射率的红波段和近红外波段对所标示的陆地像元计算归一化植被指数；

地表类型判断模块，用于通过所述归一化植被指数判断出所述陆地像元的地表类型，所述陆地像元的地表类型分为：土壤类型、植被类型和过渡区类型，并根据经空间重采样后的土壤分类图对所述陆地像元进行土壤类型的划分；

建立换算关系模块，用于通过所述ASTER波谱库数据、以及MODIS UCSB发射率波谱库数据，建立窄波段发射率与宽波段发射率之间的换算关系，得到水体像元和冰雪像元宽波段发射率；

计算模块，用于基于ASTER宽波段发射率数据，分别计算出陆地上土壤类型、植被类型和过渡区类型三个类型像元的宽波段发射率。

(三)有益效果

本发明提供了一种适用于全球陆表大尺度范围的宽波段发射率反演算法，通过对反射率数据和反照率数据的处理，提高了宽波段发射率反演的精度。

附图说明

图1是按照本发明一种实施方式的宽波段发射率反演方法的流程图；

图2a是波段为8-13.5μm时，通过图1所示的方法计算的宽波段发射率和ASTER宽波段发射率的散点图；

图2b是波段为8-13.5μm时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的偏差直方图；

图3a是波段为3-14μm时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的散点图；

图3b是波段为3-14μm时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的偏差直方图；

图4a是波段为8-13.5μm，且土壤类型为Andisols时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的散点图；

图4b是波段为8-13.5μm，且土壤类型为Andisols时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的偏差直方图；

图5a是波段为3-14μm，且土壤类型为Andisols时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的散点图；

图5b是波段为3-14μm，且土壤类型为Andisols时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的偏差直方图；

图6a是波段为8-13.5μm，且土壤类型为Inceptisols时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的散点图；

图6b是波段为8-13.5μm，且土壤类型为Inceptisols时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的偏差直方图；

图7a是波段为3-14μm，且土壤类型为Inceptisols时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的散点图；

图7b是波段为3-14μm，且土壤类型为Inceptisols时，ASTER宽波段发射率和通过图1所示的方法计算的宽波段发射率的偏差直方图；

图8a是MODIS波段1反射率和AVHRR波段1反射率之间的关系示意图；

图8b是MODIS波段2反射率和AVHRR波段2反射率之间的关系示意图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

图1是按照本发明一种实施方式的宽波段发射率反演方法的流程图，包括以下步骤：

S1：通过中分辨率成像光谱仪(Moderate Resolution ImagingSpectroradiometer，MODIS)或改进型甚高分辨率辐射仪(AdvancedVery High Resolution Radiometer，AVHRR)卫星遥感传感器获取多种分辨率的地表反射率和反照率，并根据所述地表反射率和反照率获取多种分辨率的以及配套的ASTER宽波段发射率数据、ASTER波谱库数据、MODIS UCSB发射率波谱库数据和土壤分类图数据；

本实施方式中通过MODIS卫星遥感传感器获取2000年至2010年的遥感数据，通过所述AVHRR传感器获取1985年至2000年的遥感数据。

分析全球每天的遥感图像，全球平均70％以上地区都被云覆盖；云覆盖存在时间上的长期性、季节性、多变性；同时，云覆盖带来了云影的存在；中高纬度地区存在大量的可溶性雪，云和雪又存在较大相似性。因此，遥感数据所反映的地表反射率和反照率往往受到干扰，从而很难精确地反映出地表特征参量产品的变化规律。因此，需对发射率算法的输入数据进行预处理。

S2：对所述多种分辨率的地表反射率和反照率分别进行预处理。

步骤S2具体包括以下步骤：

S21：对所述多种分辨率的地表反射率和反照率进行判断，初步筛选出正常的像元和异常的像元；

S22：从正常的像元中选取已被卫星遥感传感器标识好的像元作为训练样本，分别计算所述异常的像元和正常的像元与训练样本之间的相关系数或相似系数，并判断计算出的相关系数或相似系数是否大于或等于相关系数阈值或相似系数阈值，若是，则判定为正常的像元，否则判定为异常的像元，优选地，相关系数阈值为0.9，相似系数阈值为0.98；

S23：基于地理位置、时间和归一化雪被指数(NormalizedDifference Snow Index，NDSI)判断云和雪，具体为：

若NDSI＞0.5，且地理位置和时间符合下雪条件的，标识为纯雪，其中，NDSI＝(R4-R6)/(R4+R6)，R4为波长为0.555微米的地表反射率，R6为波长为1.64微米的地表反射率；

若0.4＜NDSI＜0.5，则根据已有云雪标识的像元作为训练样本，建立云和雪两类。根据训练结果，利用最大似然法对无法识别的云和雪分类；

若NDSI＜0.4，由于有可能是云，也有可能是云和雪的混合，则都判为云；

S24：通过时空插值滤波的方式来填充在长时间序列中缺失的像元(由于接收处理失败或数据质量可能存在问题，会缺少某一天或几天的数据)和空间上异常的像元，具体为：利用一年时间序列内的数据，根据同一类地物光谱在时间和空间上的连续性和相关性特性，采用多项式拟合的方法进行填充插值。

S3：对所述土壤分类图数据进行空间重采样，对经预处理后的多种分辨率的地表反射率和反照率进行空间匹配，并将进行空间匹配后的地表反照率逐个像元标示出水体像元、冰雪像元和陆地像元，通过进行空间匹配后的地表反射率的红波段和近红外波段对所标示的陆地像元计算归一化植被指数；

对0.03°×0.03°土壤分类图重采样和重投影，使之与MODIS反照率产品投影相同，均为正弦(sinusoidal)投影，得到lkm的土壤分类图，本实施方式中，针对于MODIS数据，利用MODIS反照率质量控制产品MCD43B2的数据集Snow_BRDF_Albedo和BRDF_Albedo_Ancillary对MODIS反照率逐个像元区分水体、冰/雪和陆地像元，对于陆地像元利用MODIS反射率的红光波段和近红外波段的值，计算NDVI，针对AVHRR数据，利用AVHRR反照率质量控制产品AVH09C1中的水体标示对AVHRR反照率逐个像元区分水体、冰/雪和陆地像元，对于陆地像元利用AVHRR反射率的红光波段、近红外波段的值，计算其NDVI。

NDVI的计算公式为

NDVI＝(R(nir)-R(red))/R(nir)+R(red))

其中R(nir)表示地物近红外波段的反射率，R(red)表示地物红光波段的反射率。

S4：通过所述归一化植被指数判断出所述陆地像元的地表类型，所述陆地像元的地表类型分为：土壤类型、植被类型和过渡区类型，并根据经空间重采样后的土壤分类图对所述陆地像元进行土壤类型的划分。

对陆地像元利用NDVI阈值法划分土壤类型(NDVI≤0.1)、植被覆盖类型(NDVI≥0.2)和过渡区类型(0.1＜NDVI＜0.2)，再分别针对MODIS数据和AVHRR数据，根据经空间重采样后的土壤分类图等数据对所述陆地像元逐像元进行土壤类型的划分，以标识出不同的土壤类型，如表1为陆地像元标识。

表1

S5：通过所述ASTER波谱库数据、以及MODIS UCSB发射率波谱库数据，建立窄波段发射率与宽波段发射率之间的换算关系，得到水体像元和冰雪像元宽波段发射率；

由于MODIS和AVHRR热红外传感器的波段及波段的带宽有限，无法直接得到3-14μm和8-13.5μm两个宽波段发射率，因此需要使用ASTER波谱库和MODIS UCSB波谱库中的波谱数据，建立窄波段发射率和宽波段发射率之间的换算关系，由此得到MODIS和AVHRR两个传感器数据的水体像元和冰雪像元宽波段发射率。

由于ASTER波谱数据库中的数据为半球方向光谱反射率，利用基尔霍夫定律(ε_λ＝1-ρ_λ)将其转化为窄波段发射率ε(λ)，则宽波段发射率可定义如下式：

$\stackrel{\‾}{\mathrm{\ϵ}}=\frac{{\∫}_{{\mathrm{\λ}}_1}^{{\mathrm{\λ}}_2}\mathrm{\ϵ}\left(\mathrm{\λ}\right)B(\mathrm{\λ},T)\mathrm{d\λ}}{{\∫}_{{\mathrm{\λ}}_1}^{{\mathrm{\λ}}_2}B(\mathrm{\λ},T)\mathrm{d\λ}}---\left(1\right)$

其中，λ₁和λ₂分别为宽波段发射率的波长上限和下限，B(λ，T)为普朗克函数，T为地物温度。

对于水体像元，考虑到发射率光谱变化很小，可以近似为常数。水体像元的发射率值可以根据ASTER波谱库和MODIS UCSB波谱库中水体波谱计算获得。基于水体发射率光谱数据，根据上式(1)计算得到波段为8-13.5μm和3-14μm的宽波段发射率。ASTER光谱中样本seawater和tapwater的宽波段发射率皆为0.984；MODIS UCSB光谱中样本distd_wn，diswat06，seagol01，seagol02，seawat10的宽波段发射率皆为0.985。因此，水体的宽波段发射率可统一赋值为0.985。

冰雪像元的发射率光谱变化与水体像元类似，采取同样的方法计算其宽波段发射率。

ASTER波谱库中样本fine snow，medium snow，coarse snow，ice，frost的宽波段发射率分别为0.99，0.986，0.982，0.983，0.984。

MODIS UCSB光谱中的不同冰样本的宽波段发射率分别为0.982，0.983，0.983；雪的宽波段发射率分别为0.986，0.99。

本实施方式反演了某地不同粒径雪、不同观测角度的光谱发射率，采取与水体相同的方法将其转化为宽波段发射率，可得到不同粒径雪不同观测角度发射率光谱对应的宽波段发射率，如表2所示。

观测角度	0°	15°	30°	45°	60°	75°
							Ice	0.983	0.982	0.981	0.978	0.953	0.861
Coarse	0.986	0.986	0.987	0.986	0.978	0.971
							Crust	0.983	0.984	0.984	0.981	0.974	0.954
Fine	0.988	0.989	0.989	0.988	0.987	0.986
							medium	0.987	0.987	0.986	0.987	0.984	0.984

表2

显然，冰/雪宽波段发射率随观测角度和粒径的变化比较明显。虽然给定雪的等效粒径，可以使用辐射传输模型来模拟雪的发射率光谱，冰的发射率光谱也可由菲涅耳定律模拟。现实中，很难获取雪的粒径。此外，冰表面含有水或风蚀特征时，无法进行模拟。观测视角通常小于45度，冰/雪的宽波段发射率赋值为0.985，带来的误差小于0.005。

S6：基于ASTER宽波段发射率数据，分别计算出陆地上土壤类型、植被类型和过渡区类型三个类型像元的宽波段发射率；

(1)像元为土壤类型

若所述陆地像元的地表类型判定为土壤类型，则通过下列方式进行计算：

使用的数据源来自于实际的卫星产品，即ASTER宽波段发射率数据，并提取出同时期的MODIS的7个波段的反照率数据。根据土壤分类图，选择试验区，表3列出了试验区的地理位置、土壤类型以及数据获取的时间范围。

土壤类型	地理位置	时间范围
			Alfisols	Lat：38.67°-41.77°；Lon：-86.7°--80.03°	2009.01-2009.03
Andisols	Lat：18.3°-20.7°；Lon：-103.3°--96.7°	2009.02-2009.03
			Aridisols	Lat：26.7°-43.4°；Lon：-120°--103.3°	2009.01-2009.02
Entisols	Lat：13.3°-26.7°；Lon：-3.3°-30°	2009.02-2009.02
			Gelisols	Lat：30.03°-36.7°；Lon：86.63°-96.63°	2008.01-2008.01
Inceptisols	Lat：37.61°-41.36°；Lon：-83.24°-74.51°	2009.01-2009.03
			Mollisols	Lat：30°-53.4°；Lon：-113.3°--93.4°	2009.03-2009.03
Oxidols	Lat：-27.5°--22°；lon：-55°--50°	2009.07-2009.09
			Ultisols	Lat：30°-38.4°；Lon：-96.7°--75°	2008.01-2009.02
Vertisols	Lat：13.3°-26.7°；Lon：73.3°-80°	2009.02-2009.03

表3

对于每个试验区提取土壤像元的反照率和ASTER宽波段发射率，进行线性最小二乘拟合计算，得到如表4所示的函数关系，表达式中a是反照率albedo的缩写，a1至a7分布表示MODIS的7个可见光-近红外波段的黑天空反照率，具体的波段如下表：

反照率	波段(nm)
		a1	620-670
a2	841-876
		a3	459-479
a4	545-565
		a5	1230-1250
a6	1628-1652
		a7	2105-2155

其中，波段的选择是通过多次试验确定的，例如对于Andisols来说，8-13.5宽波段发射率与MODIS的第1，3和7波段的反照率之间的关系比其他的组合关系更好。

8种土壤类型(Alfisols，Aridisols，Gelisols，Mollisols，Oxisols，Ultisols，Histosols，Spodosols)的宽波段发射率Emis计算表达式非常相近，将其合并在一起称为Ensemble-8，总共得到表4所示的3套公式，其余未分土壤类型均采用Ensemble-8的表达式计算。

表4

根据图2a、2b、4a、4b、6a和6b可观察出：波段为8-13.5μm时，通过图1所示方法计算的宽波段发射率和ASTER波谱库的宽波段发射率之间的相关系数大于0.60，均方根误差小于0.01，偏差很小基本上可以忽略。

根据图3a、3b、5a、5b、7a和7b可观察出：波段为3-14μm时，通过图1所示方法计算的宽波段发射率和ASTER波谱库的宽波段发射率之间的相关系数大于0.64，均方根误差小于0.008，偏差很小基本上可以忽略。表5中列出了图2a和2b中每种土壤类型对应的均方根误差，均方根误差小于0.013。

土壤类型	8-13.5μm均方根误差	3-14μm均方根误差
			Alfisols	0.009	0.013
Aridisols	0.009	0.008
			Entisols	0.011	0.010
Gelisols	0.010	0.008
			Mollisols	0.006	0.006
Oxidols	0.015	0.010
			Ultisols	0.009	0.007
Vertisols	0.007	0.006

表5

对于AVHRR数据，其在可见/近红外波段只有2个波段，没有反照率。针对MODIS的土壤发射率算法不适用于AVHRR。AVHRR的第1、2波段和MODIS的第1、2波段中心波长相近，可建立MODIS的第1、2波段反射率和ASTER宽波段发射率之间的非线性关系，如果关系不显著，则借鉴NDVI阈值法，对土壤宽波段发射率进行赋值；如果关系显著，建立MODIS的第1、2波段反射率与AVHRR的第1、2波段反射率的转换关系，得到AVHRR反射率和ASTER宽波段发射率之间的联系，用于AVHRR的宽波段发射率反演。

使用JHU波谱库中的岩石、土壤、植被光谱，分别与AVHRR和MODIS的第1、2波段的响应函数卷积，得到相应的波段反射率，进行回归，得到如图8a和8b所示的结果，两者呈非常好的线性相关，相关性达到0.999以上。因此，认为MODIS反射率和AVHRR反射率是可以近似等价的，直接将MODIS反射率和发射率之间的关系用于AVHRR，不加修正。直接使用AVHRR数据反演宽波段发射率时，使用AVHRR的NDVI划分地表为土壤或者植被覆盖。

8种土壤类型(Alfisols，Aridisols，Gelisols，Mollisols，Oxisols，Ultisols，Histosols，Spodosols)的宽波段发射率计算表达式非常相近，将其合并在一起称为Ensemble-8，表6给出了不同土壤类型宽波段发射率和反射率的非线性关系，表6中，R1为AVHRR的红光波段的反射率，R2为近红外波段的反射率。

表6

表7给出了六种土壤类型的宽波段发射率均值、偏差及均方根误差，宽波段发射率均值位于0.959和0.972之间，变化幅度达0.13，根据土壤类型分别确定其均值要优于NDVI阈值法中对所有的类型指定统一的值。

土壤类型	均值	均方根误差
			Alfisols	0.965	0.006
Aridisols	0.964	0.010
			Gelisols	0.959	0.012
Mollisols	0.967	0.013
			Oxidols	0.968	0.001
Ultisols	0.972	0.003

表7

(2)像元为植被类型

若所述陆地像元的地表类型判定为植被类型，则通过下列方式进行计算：

建立植被覆盖地区ASTER宽波段发射率与MODIS NDVI、反照率的经验关系，用于宽波段发射率反演。在全球范围内选择40个站点，下载2000-2010年晴空下的ASTER瞬时发射率产品与同期MODIS的NDVI产品、7个波段的反照率产品MCD43B3及质量控制产品MCD43B2，建立宽波段发射率与MODIS的NDVI、反照率的经验关系。本实施方式设定1km尺度上宽波段发射率与NDVI、反照率之间存在线性关系，如下式：

$\stackrel{\‾}{\mathrm{\ϵ}}=C+A_0\·\mathrm{NDVI}+\underset{i=1}{\overset{7}{\mathrm{\Σ}}}{A}_i{\mathrm{\α}}_i$

其中，A₀和A_i分别为NDVI和反照率的系数，C为常数。

针对MODIS数据，将植被类型分为了5类，即森林区(Forest)、灌木区(Shrublands)、草原区(Savannas)、草地区(Grassland)和耕地区(Cropland)，其中森林区包括落叶林、阔叶林和混交林，灌木区包括浓密灌木和稀疏灌木。从40个站点的ASTER与MODIS数据中选择出60656个有效样本数据，各种植被类型占有的样本数目和比例如表8所示。理论上，这些样本数据对全球植被的发射率具有很强的代表性。

植被类型	像元个数	所占比例(％)
			Shrublands	13954	23.0
Savannas	11661	19.2
			Grassland	12593	20.8
Cropland	13373	22.0
			Forest	9075	15.0

表8

为了高效准确的反演的全球地表宽波段发射率，不给出针对于各类植被类型的经验关系，而是基于所有样本数据建立一套普适性的求取宽波段发射率的公式。表9给出了利用MODIS的NDVI、7个波段反照率求取宽波段发射率的经验关系系数和相关统计信息。

表9

经计算，利用经验关系计算的宽波段发射率与样本实测值之间的均方根误差(RMSE)集中在0.01左右，其中3-14μm宽波段发射率的RMSE等于0.008，8-13.5μm宽波段发射率的RMSE等于0.01。

针对AVHRR数据，由于其在可见光/近红外波段并未向MODIS那样设置了7个波段，而只是设置了红光和近红外。因此在植被区，仍以上文提及的40个站点的数据为基础，首先将MODIS的红光、近红外反射率转化为AVHRR的红光和近红外波段的反射率，然后建立起宽波段发射率与AVHRR的红光、近红外以及NDVI的经验关系，用于生成AVHRR的全球宽波段发射率产品。研究发现，宽波段发射率与AVHRR的红光、近红外反射率之间的非线性关系能够使结果具有较小的误差，表10给出了AVHRR的NDVI、红光和近红外求取3-14μm和8-13.5μm宽波段发射率的经验关系系数和相关统计信息。

波段	表达式
		3-14μm	Emis＝0.976-0.055*R1-0.032*R2+0.495*R1*R1+0.068*R2*R2-0.519*R1*R2
8-135μm	Emis＝0.983-0.066*R1-0.049*R2+0.66*R1*R1+0.116*R2*R2-0.759*R1*R2

表10

其中，经计算，受波段数目的限制，仅仅从AVHRR红光、近红外反射率求取宽波段发射率的精度整体上小于MODIS的算法，但是两个波长范围内的宽波段发射率的RMSE均小于0.015。

(3)过渡区类型

若所述陆地像元的地表类型判定为过渡区类型，则根据下列方式进行计算宽波段发射率：

根据NDVI值对地表进行分类，具有不确定性，尤其对于稀疏植被地表。为此，根据NDVI值，确定土壤和植被的过渡区域。将NDVI值位于0.1和0.2之间的地区，标识为过渡区域。最终的过渡区域宽波段发射率值根据NDVI的归属确定。过渡区域的像元会同时属于两个地表类型，土壤和土壤过渡区域，或者植被与植被过渡区域，对两者分别计算的宽波段发射率进行平均，达到最终的宽波段发射率。

使用MODIS数据时，建立如表11所示的函数关系。8种土壤类型(Alfisols，Aridisols，Gelisols，Mollisols，Oxisols，Ultisols，Histosols，Spodosols)的宽波段发射率计算表达式非常相近，将其合并在一起称为Ensemble-8，总共得到表12所示的3套公式。其余未分土壤类型时采用Ensemble-8的公式。

表11

使用AVHRR数据时，建立过渡区域AVHRR宽波段发射率确定方法。表12给出了不同土壤类型宽波段发射率和反射率的非线性关系。表13给出了6种土壤类型的宽波段发射率的均值和均方根误差，宽波段发射率均值位于0.959和0.972之间，变化幅度达0.13，根据土壤类型分别确定其均值要好于NDVI阈值法中对所有的类型指定统一的值，表14为土壤类型对应的均方根误差，均方根误差小于0.011，总体上偏差绝对值很小，可以忽略。

表12

土壤类型	均值	均方根误差
			Alfisols	0.965	0.006
Aridisols	0.964	0.010
			Gelisols	0.959	0.012
Mollisols	0.967	0.013
			Oxidols	0.968	0.001
Ultisols	0.972	0.003

表13

土壤类型	8-13.5μm均方根误差	3-14μm均方根误差
			Ensemble-8	0.010	0.010
Incepti sol s	0.004	0.011
			Andisols	0.008	0.006

表14

步骤S6之后还可基于不变目标，对本实施方式中26年长时间序列发射率产品进行一致性检验与优化，最终获得全球陆表高精度的宽波段发射率。

针对1985年至1999年使用AVHRR数据获得的宽波段发射率和2000年至2010上使用MODIS数据获得的宽波段发射率，在全球选择大的沙漠作为不变目标，分析26年长时间序列宽波段发射率变化趋势，进行一致性校正，确保1999年至2000年数据的连续性，如果发现数据存在异常，则进行线性优化处理，最终得到高精度、并且具有时间一致性的长时间序列宽波段发射率。

本发明还公开了一种全球陆表宽波段发射率反演系统，包括：

获取模块，用于通过MODIS或AVHRR卫星遥感传感器获取多种分辨率的地表反射率和反照率，并根据所述地表反射率和反照率获取多种分辨率的以及配套的ASTER宽波段发射率数据、ASTER波谱库数据、MODIS UCSB发射率波谱库数据和土壤分类图数据；

预处理模块，用于对所述多种分辨率的地表反射率和反照率分别进行预处理；

匹配标示模块，用于对所述土壤分类图数据进行空间重采样，对经预处理后的多种分辨率的地表反射率和反照率进行空间匹配，并将进行空间匹配后的地表反照率逐个像元标示出水体像元、冰雪像元和陆地像元，通过进行空间匹配后的地表反射率的红波段和近红外波段对所标示的陆地像元计算归一化植被指数；

地表类型判断模块，用于通过所述归一化植被指数判断出所述陆地像元的地表类型，所述陆地像元的地表类型分为：土壤类型、植被类型和过渡区类型，并根据经空间重采样后的土壤分类图对所述陆地像元进行土壤类型的划分；

建立换算关系模块，用于通过所述ASTER波谱库数据、以及MODIS UCSB发射率波谱库数据，建立窄波段发射率与宽波段发射率之间的换算关系，得到水体像元和冰雪像元宽波段发射率；

计算模块，用于基于ASTER宽波段发射率数据，分别计算出陆地上土壤类型、植被类型和过渡区类型三个类型像元的宽波段发射率。

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

Claims (4)

1.一种全球陆表宽波段发射率反演方法，其特征在于，包括以下步骤：

S1：通过MODIS或AVHRR卫星遥感传感器获取多种分辨率的地表反射率和反照率，并根据所述地表反射率和反照率获取多种分辨率的以及配套的ASTER宽波段发射率数据、ASTER波谱库数据、MODIS UCSB发射率波谱库数据和土壤分类图数据；

S2：对所述多种分辨率的地表反射率和反照率分别进行预处理；

S3：对所述土壤分类图数据进行空间重采样，对经预处理后的多种分辨率的地表反射率和反照率进行空间匹配，并将进行空间匹配后的地表反照率逐个像元标示出水体像元、冰雪像元和陆地像元，通过进行空间匹配后的地表反射率的红波段和近红外波段对所标示的陆地像元计算归一化植被指数；

S4：通过所述归一化植被指数判断出所述陆地像元的地表类型，所述陆地像元的地表类型分为：土壤类型、植被类型和过渡区类型，并根据经空间重采样后的土壤分类图对所述陆地像元进行土壤类型的划分；

S5：通过所述ASTER波谱库数据、以及MODIS UCSB发射率波谱库数据，建立窄波段发射率与宽波段发射率之间的换算关系，得到水体像元和冰雪像元宽波段发射率；

S6：基于ASTER宽波段发射率数据，分别计算出陆地上土壤类型、植被类型和过渡区类型三个类型像元的宽波段发射率。

2.如权利要求1所述的宽波段发射率反演方法，其特征在于，步骤S2具体包括：

S21：对所述多种分辨率的地表反射率和反照率进行判断，初步筛选出正常的像元和异常的像元；

S22：从正常的像元中选取已被卫星遥感传感器标识好的像元作 为训练样本，分别计算所述异常的像元和正常的像元与训练样本之间的相关系数或相似系数，并判断计算出的相关系数或相似系数是否大于或等于相关系数阈值或相似系数阈值，若是，则判定为正常的像元，否则判定为异常的像元；

S23：基于地理位置、时间和归一化雪被指数判断云和雪；

S24：通过时空插值滤波的方式来填充在长时间序列中缺失的像元和空间上异常的像元。

3.如权利要求2所述的宽波段发射率反演方法，其特征在于，所述相关系数阈值为0.9，相似系数阈值为0.98。

4.一种全球陆表宽波段发射率反演系统，其特征在于，包括：

获取模块，用于通过MODIS或AVHRR卫星遥感传感器获取多种分辨率的地表反射率和反照率，并根据所述地表反射率和反照率获取多种分辨率的以及配套的ASTER宽波段发射率数据、ASTER波谱库数据、MODIS UCSB发射率波谱库数据和土壤分类图数据；

预处理模块，用于对所述多种分辨率的地表反射率和反照率分别进行预处理；

匹配标示模块，用于对所述土壤分类图数据进行空间重采样，对经预处理后的多种分辨率的地表反射率和反照率进行空间匹配，并将进行空间匹配后的地表反照率逐个像元标示出水体像元、冰雪像元和陆地像元，通过进行空间匹配后的地表反射率的红波段和近红外波段对所标示的陆地像元计算归一化植被指数；

地表类型判断模块，用于通过所述归一化植被指数判断出所述陆地像元的地表类型，所述陆地像元的地表类型分为：土壤类型、植被类型和过渡区类型，并根据经空间重采样后的土壤分类图对所述陆地像元进行土壤类型的划分；

建立换算关系模块，用于通过所述ASTER波谱库数据、以及MODIS UCSB发射率波谱库数据，建立窄波段发射率与宽波段发射率 之间的换算关系，得到水体像元和冰雪像元宽波段发射率；

计算模块，用于基于ASTER宽波段发射率数据，分别计算出陆地上土壤类型、植被类型和过渡区类型三个类型像元的宽波段发射率。 

Images