CN102435586B - 地表反照率产品的生成方法及系统 - Google Patents

Abstract

本发明公开了一种地表反照率产品的生成方法及系统，涉及卫星遥感领域。所述方法包括步骤：根据经过分类的POLDER BRDF数据，得到POLDER反照率数据，转换POLDER反照率数据得到宽波段地表反照率数据；根据样本MODIS地表反射率数据和宽波段地表反照率数据建立第一转换模型，得到第一地表反照率数据；根据样本MODIS表观反射率数据和宽波段地表反照率数据建立第二转换模型，得到第二地表反照率数据；整合第一和第二地表反照率数据，及经网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。本发明生产出了长时间序列并且时空连续、高时空分辨率的地表反照率产品，提高了地表反照率产品反演的精度。

Description

地表反照率产品的生成方法及系统

技术领域

本发明涉及卫星遥感技术领域，特别涉及一种地表反照率产品的生成方法及系统。

背景技术

地表反照率反映了地球表面对太阳光辐射的反射能力，是地表辐射能量平衡以及地气相互作用中的驱动因子之一，其时空变化受到自然过程以及人类活动的影响，是全球环境变化的指示因子。

目前一些卫星反照率产品已经业务化运行和发布，空间分辨率从250m-20km，时间分辨率从日到月，其中极轨卫星产品包括MODIS，MISR，CERES，POLDER，MERIS。地球同步卫星产品包括Meteosat，and MSG。

目前针对在轨运行的卫星遥感数据以及部分航空遥感数据进行了大量地表反照率反演研究，并进行了有效的反照率产品的地面验证，其中基于线性核驱动模型的地表BRDF/反照率遥感反演模型算法是目前地表反照率遥感反演中应用最广泛的方法，已经在MODIS、MISR的地表反照率产品中得到了很好的应用。另外，POLDER/PARASOL系列传感器具有更好的多角度观测能力，空间分辨率稍低(6km)，也发布了非常有特色的全球长时间序列反照率产品。其他如静止轨道气象卫星MSG、METEOSET，极轨卫星传感器AVHRR、VEGETATION等传感器都有不同覆盖范围的反照率产品，我国气象卫星、环境减灾小卫星数据都可用于区域和全球范围的地表反照率反演，目前反照率产品生成算法正在研制过程中。

已有的反照率产品的生产系统环境，均是采用个人计算机少量生产，并且缺少长时间序列的全球陆表宽波段发射率数据集。迄今为止，针对全球陆面变化研究与陆面模型研发，国际陆地遥感领域仍然缺乏长时间序列并且时空连续、高时空分辨率的全球陆表特征参量产品。国内的遥感产品生产均采用个人计算机生产小批量产品，无法满足长时间序列、高时空分辨率和高质量的遥感生产需求。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：如何提供一种时空连续的地表反照率产品的生成方法及系统。

(二)技术方案

为解决上述技术问题，本发明提供一种地表反照率产品的生成方法，其包括步骤：

B：根据太阳天顶角、观测天顶角和相对方位角对MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据进行网格划分；

C：对POLDER BRDF数据进行优选，按照地表类型对经过优选的POLDER BRDF数据进行分类；

D：根据配套光谱数据模拟得到模拟POLDER BRDF数据和模拟MODIS地表反射率数据，建立所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系，根据所述多元线性转换关系和经过分类的POLDER BRDF数据，计算得到样本MODIS地表反射率数据；

E：根据经过分类的POLDER BRDF数据，得到POLDER反照率数据，转换所述POLDER反照率数据得到宽波段地表反照率数据；

F：根据所述样本MODIS地表反射率数据和所述宽波段地表反照率数据建立第一转换模型，根据所述第一转换模型和经过网格划分的MODIS地表反射率数据，得到第一地表反照率数据；

G：根据所述POLDER BRDF数据得到地表反射特性参数，根据所述地表反射特性参数和大气状态参数计算得到样本MODIS表观反射率数据；

H：根据所述样本MODIS表观反射率数据和所述宽波段地表反照率数据建立第二转换模型，根据所述第二转换模型和经过网格划分的MODIS表观反射率数据，得到第二地表反照率数据；

I：整合所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。

优选地，在所述步骤B之前还包括步骤A：输入POLDER BRDF数据、MODIS地表反射率数据、MODIS表观反射率数据、MODIS地表反照率数据和配套光谱数据。

优选地，所述步骤I具体包括步骤：

I1：对所述第一地表反照率数据和所述第二地表反照率数据进行插值弥补；

I2：对所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据进行归一化处理；

I3：对归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据进行数据融合，得到融合数据；

I4：对所述融合数据进行平滑处理；

I5：对平滑处理后的融合数据进行反归一化处理，得到所述时空连续的地表反照率产品。

优选地，所述反归一化处理的公式如下：

$\mathrm{\ξ}={\mathrm{\ξ}}^{\′}\stackrel{\‾}{\mathrm{\σ}}-\stackrel{\‾}{\mathrm{\μ}}$

其中，ξ是所述平滑处理后的融合数据，ξ是所述时空连续的地表反照率产品，

是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的均值，

是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的消除随机误差后的标准差的平均值。

优选地，所述步骤C具体包括步骤：

C1：从所述POLDER BRDF数据中剔除由于云和气溶胶影响以及地表状态发生变化而不满足二向反射模型的数据；

C2：对经过所述步骤C1处理后的POLDER BRDF数据，按照植被、冰雪和裸地三种地表类型进行分类。

优选地，所述步骤C2具体包括步骤：

C21：判断POLDER BRDF数据的NDVI值是否大于0.2，如果是，则判定所述POLDER BRDF数据对应的像元为植被，否则执行步骤C22；

C22：判断POLDER BRDF数据是否其蓝光波段反射率大于0.3或者红光波段反射率大于0.3，如果是，则判定所述POLDER BRDF数据对应的像元为冰雪，否则，判定所述POLDER BRDF数据对应的像元为裸地。

优选地，所述步骤E具体包括步骤：

E1：对所述经过分类的POLDER BRDF数据进行半球积分，得到各波段的POLDER反照率数据；

E2：根据窄波段向宽波段的转换公式，转换所述POLDER反照率数据得到所述宽波段地表反照率数据。

本发明还提供一种地表反照率产品的生成系统，其包括：

网格划分模块：用于根据太阳天顶角、观测天顶角和相对方位角对MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据进行网格划分；

POLDER BRDF数据预处理模块：用于对POLDER BRDF数据进行优选，按照地表类型对经过优选的POLDER BRDF数据进行分类；

样本MODIS地表模块：用于根据配套光谱数据模拟得到模拟POLDER BRDF数据和模拟MODIS地表反射率数据，建立所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系，根据所述多元线性转换关系和经过分类的POLDERBRDF数据，计算得到样本MODIS地表反射率数据；

POLDER转换模块：用于根据经过分类的POLDER BRDF数据，得到POLDER反照率数据，转换所述POLDER反照率数据得到宽波段地表反照率数据；

第一产品模块：用于根据所述样本MODIS地表反射率数据和所述宽波段地表反照率数据建立第一转换模型，根据所述第一转换模型和经过网格划分的MODIS地表反射率数据，得到第一地表反照率数据；

样本MODIS表观模块：用于根据所述POLDER BRDF数据得到地表反射特性参数，根据所述地表反射特性参数和大气状态参数计算得到样本MODIS表观反射率数据；

第二产品模块：用于根据所述样本MODIS表观反射率数据和所述宽波段地表反照率数据建立第二转换模型，根据所述第二转换模型和经过网格划分的MODIS表观反射率数据，得到第二地表反照率数据；

终产品模块：用于整合所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。

优选地，所述系统还包括：

数据接收模块：用于接收POLDER BRDF数据、MODIS地表反射率数据、MODIS表观反射率数据、MODIS地表反照率数据和配套光谱数据。

(三)有益效果

本发明所述地表反照率产品的生成方法及系统，以现有反照率产品为基础，通过数据融合技术，生产出长时间序列并且时空连续、高时空分辨率的地表反照率产品，对于全球陆面变化研究和陆面模型研发具有积极意义。

附图说明

图1是本发明实施例所述地表反照率产品的生成方法流程图；

图2是MCD43产品与第一地表反照率数据对比的散点图；

图3是Bondville站点2006年反照率地面测量结果与第二地表反照率数据的对比图；

图4是Flagstaff-Wildfire站点2006年反照率地面测量结果与第二地表反照率数据的对比图；

图5是Willow_Creek站点不同方法生成的地表反照率产品整合后的结果对比示意图。

具体实施方式

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

图1是本发明实施例所述地表反照率产品的生成方法流程图。如图1所示，所述方法包括：

步骤A：输入以下数据：

(1)POLDER BRDF数据

POLDER BRDF数据的空间分辨率为(6km×7km)，是目前能获得的最新多角度卫星遥感数据，有丰富的角度、光谱和极化信息，是多角度遥感的理想数据之一。输入的POLDER BRDF数据共包括13227个数据文件。

(2)MODIS地表反射率数据

MODIS地表反射率数据是指MOD09GA产品。MODIS地表反射率数据中含有MODIS的1-7波段的500米分辨率的反射率数据和1km分辨率的观测天顶角、观测方位角、太阳天顶角、太阳方位角和质量控制等信息。

(3)MODIS表观反射率数据

MODIS表观反射率数据是指MODIS大气层顶反射率数据，即MODIS的L1B数据，其中含有MODIS的36个波段的数据以及相应的辐射定标信息和地理定位信息。

(4)配套光谱数据

配套的光谱数据共有493个，分别是梁顺林教授2004年著作《定量遥感》所附光盘中提供的119条波谱、“我国典型地物标准波普数据库”提供的224条植被和土壤波谱、黑河综合遥感联合实验采集的103条典型地物波谱和格林兰采集的47条冰雪波谱数据。

(5)MODIS地表反照率数据

MODIS地表反照率数据(MCD43)是利用半经验的线性核驱动模型及16天合成的多角度多波段观测数据反演得到的地表反照率产品。是Terra和Aqua星数据共同合成的MODIS Level3标准数据产品之一。以MCD43B3产品为例，它具有1km的空间分辨率以及8天的时间分辨率，其数据产品中包含MODIS传感器1-7波段的7个窄波段以及可见光(0.3-0.7μm)、近红外(0.7-5.0μm)及短波波段(0.3-5.0μm)三个宽波段的白空反照率和黑空反照率。

步骤B：根据太阳天顶角、观测天顶角和相对方位角对所述MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据进行网格划分。

太阳天顶角以2度间隔进行划分，范围是0-80度，网格中心点分别是0，2，4.....度，共分为41个间隔。观测天顶角以2度间隔进行划分，范围是0-64度，网格中心点分别是0，2，4.....度，共分为33个间隔。相对方位角以5度间隔进行划分，范围是0-180度，网格中心点分别是0，5，10.....度，共分为37个间隔。因此根据太阳天顶角、观测天顶角和相对方位角将所述MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据分别分成41*33*37＝50061个网格。

步骤C：对所述POLDER BRDF数据进行优选，按照地表类型对经过优选的POLDER BRDF数据进行分类。

所述步骤C具体包括：

步骤C1：从所述POLDER BRDF数据中剔除由于云和气溶胶影响以及地表状态发生变化而不满足二向反射模型的数据。

POLDER BRDF数据总体质量非常好，但是部分数据由于云和气溶胶影响以及地表状态发生变化(比如降雨、降雪等过程)，不再满足二向反射模型的假设。因此需要对POLDER BRDF数据进行筛选，剔除不适合做训练数据的部分。

如果某个POLDER BRDF数据满足以下3个判别准则之一，则认为它是无效数据集：

490nm波段反射率的拟合均方根误差(RMSE)大于0.01，或者其除以490nm波段的平均反射率大于0.3；

6个波段反射率拟合均方根误差的和大于0.1，或者其除以6个波段的平均反射率的和大于0.2；

总观测数小于80个，或者观测的轨数小于4轨。

POLDER BRDF数据共有13227个数据文件，经过筛选，剔除4203个，剩余9024个数据文件成为优选的POLDER BRDF数据。

步骤C2：对经过所述步骤C1处理后的POLDER BRDF数据，按照植被、冰雪和裸地三种地表类型进行分类。

不同地表类型具有不同的双向反射率特征，有必要引入地物分类信息，进一步细分训练样本，减少后续步骤中线性回归模型的不确定性。

因此，我们选择直接根据遥感观测数据分类的策略，分类的具体步骤如下：

步骤C21：计算POLDER BRDF数据的NDVI(NormalizedDifference Vegetation Index，归一化植被指数)。判断NDVI是否大于0.2，如果是，则判定所述POLDER BRDF数据对应的像元为植被，否则执行步骤C22。其中，NDVI的具体计算公式为：

NDVI＝(R2-R1)/(R2+R1)    (1)

其中R2表示地物0.865微米的地表反射率值，R1表示地物0.67微米的地表反射率值。

步骤C22：判断POLDER BRDF数据是否蓝光波段(0.49微米)反射率大于0.3，或者红光波段(0.67微米)反射率大于0.3，如果是，则判定所述POLDER BRDF数据对应的像元为冰雪，否则，判断所述POLDER BRDF数据对应的像元为裸地。

步骤D：根据所述配套光谱数据模拟得到模拟POLDER BRDF数据和模拟MODIS地表反射率数据，建立所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系，根据所述多元线性转换关系和经过分类的POLDER BRDF数据，计算得到样本MODIS地表反射率数据。

所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系公式如下：

Mi＝offset+∑Ki*Pi    (2)

其中Mi是第i波段样本MODIS地表反射率数据，Ki为第i个波段的转换系数，Pi为第i个波段的经过分类的POLDER BRDF数据，offset是偏移值。生成的各波段样本MODIS地表反射率数据如表1所示。

表1各波段样本MODIS地表反射率数据表

步骤E：根据经过分类的POLDER BRDF数据，得到POLDER反照率数据，转换所述POLDER反照率数据得到宽波段地表反照率数据。

所述步骤E具体包括：

步骤E1：对所述经过分类的POLDER BRDF数据进行半球积分，得到各波段的POLDER反照率数据。

步骤E2：根据窄波段向宽波段的转换公式，转换所述POLDER反照率数据得到所述宽波段地表反照率数据。

其中，POLDER反照率数据的窄波段向宽波段转换公式为：

A＝C0+∑CiBi    (3)

其中，A是宽波段地表反照率数据，Ci为POLDER反照率数据第i个波段的转换系数，Bi为POLDER反照率数据第i个波段的反照率数据。因为雪的光谱特征与植被-土壤体系显著不同，雪的反照率向宽波段转换需要不同的转换系数，对雪覆盖和无雪覆盖地表给出不同的转换系数，见表2：

表2转换系数表

步骤F：根据所述样本MODIS地表反射率数据和所述宽波段地表反照率数据建立第一转换模型，根据所述第一转换模型和经过网格划分的MODIS地表反射率数据，得到第一地表反照率数据。

本步骤的整体思路是根据所述样本MODIS地表反射率数据和所述宽波段地表反照率数据，建立各网格内不同地物的宽波段地表反照率转换模型；再使用输入的经过网格划分的MODIS地表反射率数据，生成第一地表反照率数据。

反演第一地表反照率数据的算法可以表述为如下公式

其中α_ws为宽波段地表白空反照率，α_bs(θ_k)为宽波段地表黑空反照率，θ_k为黑空反照率的太阳天顶角，从0到80度，5度为步长共计17个间隔(θ_k＝0，5，...，75，80)，k＝1，2，3...，16，17，i为传感器波段(i＝1，2，3，4)，c_i为回归得到的回归系数，

为MODIS发布的地表反射率产品MOD09中提供的经过大气校正的方向反射率。

公式4和公式5代表了地表多波段方向反射率与宽波段地表反照率之间的关系，这个关系是由回归系数决定，对于3种地表类型以及每一个太阳/观测角度网格，就有一组回归系数。

为了查看回归效果，我们统计了使用训练数据(即由公式2计算得到的样本MODIS地表反射率数据)时的WSA(白空反照率)和45度太阳角BSA(黑空反照率)反演误差，因为共有50061个网格，我们统计了RMSE(均方根误差)的平均值，如下表3。

表3作用于训练数据时的误差统计

训练数据类别	WSA的平均RMSE	45度太阳角BSA的平均RMSE
			植被	0.0099	0.0082
裸地	0.0132	0.0120
			冰雪	0.0214	0.0170

为了用地面观测数据定量评价第一地表反照率数据的精度，我们使用计算得到的结果与地面观测的反照率数据进行比较。

直接用每日的第一地表反照率数据与地面观测值进行比较，无雪条件下共统计了5240个有效观测，第一地表反照率数据与地面观测值的系统偏差为0.002775，第一地表反照率数据略小于台站观测，均方根误差为0.04751，总体分析表明，第一地表反照率数据的平均值与地面观测平均值符合较好。

图2是MCD43产品与第一地表反照率数据(AB1产品)对比的散点图。如图2所示，我们选取北美10个地表均匀的观测站点的第一地表反照率数据和MCD43产品进行对比，时间跨度根据各站点收集数据的情况不同，最大从2001年至2009年。结果表明二者之间的一致性比较好，在总共1067个观测值中，二者的系统偏差(bias)为-0.0096，MCD43产品的数值比第一地表反照率数据的数值偏低，均方根误差(RMSE)为0.0206，相关系数为0.9538。

步骤G：根据所述POLDER BRDF数据得到地表反射特性参数，根据所述地表反射特性参数和大气状态参数计算得到样本MODIS表观反射率数据。

本步骤的本质是通过MODIS传感器大气层顶方向反射率(TOAreflectance)与地表反照率(Albedo)的训练数据集，通过划分格网回归的方法建立二者之间的统计回归关系。为了获得具有普适性的训练结果，需要建立一个涵盖全球各种二向性反射特性的地表在不同大气状况下(大气类型、气溶胶类型、气溶胶光学厚度和目标海拔高度等)的大气层顶方向反射率及其对应的地表宽波段反照率的数据集。由于直接采用6S大气辐射传输模型模拟非朗伯地表的大气层顶方向反射率计算量非常巨大，因此我们采用具有较高精度的近似公式来计算大气层顶方向反射率，所用公式如下：

${\mathrm{\ρ}}^{*}(i,v)=t_g\{{\mathrm{\ρ}}_0(i,v)+\frac{T\left(i\right)\·R(i,v)\·T\left(v\right)-t_{\mathrm{dd}}\left(i\right)\·t_{\mathrm{dd}}\left(v\right)\·|R(i,v)\·\stackrel{\‾}{\mathrm{\ρ}}|}{1-r_{\mathrm{hh}}\stackrel{\‾}{\mathrm{\ρ}}}\}---\left(6\right)$

其中矩阵T(i)，R(i，v)，T(v)分别定义为

$T\left(i\right)=\left[t_{\mathrm{dd}}\left(i\right)t_{\mathrm{dh}}\left(i\right)\right],T\left(v\right)=\left[\begin{array}{c}{t}_{\mathrm{dd}}\left(v\right)\\ t_{\mathrm{hd}}\left(v\right)\end{array}\right],$ $R(i,v)=\left[\begin{array}{cc}{r}_{\mathrm{dd}}(i,v)& r_{\mathrm{dh}}\left(i\right)\\ r_{\mathrm{hd}}\left(v\right)& r_{\mathrm{hh}}\end{array}\right]$

其中i，v分别代表太阳入射方向和传感器观测方向，ρ(i，v)为大气层顶方向反射率，ρ₀(i，v)为大气程反射率，t_g为气体吸收透过率，

为大气球形反照率，t，r分别代表大气透过率和地表反射率，其中下标h，d分别代表散射(半球)和直射(方向)。t_dd(i)，t_dh(i)，t_dd(v)，t_hd(v)分别代表大气下行直射透过率、大气下行方向半球透过率、大气上行直射透过率和大气上行半球方向透过率，当i＝j时，有t_dd(i)＝t_dd(v)，t_dh(i)＝t_hd(v)。r_dd(i，v)，r_dh(i)，r_hd(v)，r_hh分别为地物的二向性反射因子、方向半球反射率、半球方向反射率和双半球反照率，在考虑地表互易的情况下有，当i＝j时，有r_dh(i)＝r_hd(v)。其中地表反射特性参数r_dd(i，v)，r_dh(i)，r_hd(v)，r_hh为通过POLDER-3/PARASOL BRDF数据集经过波段转换和半球积分得到，作为参数输入，而大气状态参数ρ₀(i，v)，t_dd(i)，t_dh(i)，t_dd(v)，t_hd(v)，

则通过大气辐射传输模型模拟建立的大气参数查找表(Look Up Table)获得。

该大气参数查找表是一个以大气类型、气溶胶类型、气溶胶光学厚度、目标海拔高度、太阳天顶角、观测天顶角和相对方位角为维度的7维大气参数查找表，具体如下表4。在大气辐射传输模型模拟过程中，通过改变模型的输入参数来模拟各种大气状况：其中大气类型设置为热带、中纬度夏季、中纬度冬季、副极地夏季、副极地冬季和US62标准大气6种；气溶胶类型设置为大陆型气溶胶、海洋型气溶胶、城市型气溶胶、沙漠型气溶胶、生物燃烧型气溶胶和霾型气溶胶6种，其中霾型气溶胶假定气溶胶中沙尘、水溶性、烟尘和海洋粒子的组成百分比分别为15％，75％，10％和0％；550nm的气溶胶光学厚度设置为0.01，0.05，0.1，0.2共4个梯度，包含了从清洁大气到较浑浊大气的情况；目标海拔高度设置为0到3.5km，0.5km为步长共计8个梯度；太阳天顶角设置为0到80度，4度为步长共计21个间隔；观测天顶角为0到64度，4度为步长共计17个间隔；相对方位角为0到180度，20度为步长共计10个间隔。因此入射/出射几何空间(太阳天顶角、观测天顶角和相对方位角3个维度)共被划分为3570个格网(angular bin)，通过6S大气辐射传输模型计算格网中央点的大气参数代表整个格网的大气参数。经过参数敏感性分析，该格网划分方案可以满足大气辐射传输模拟精度的要求，且极大地减少了进行大气辐射传输模拟的计算量。其中水汽含量可以通过MODIS大气水汽含量产品错误！未找到引用源。获得，在模拟过程中不改变该参数。输入数据集包括POLDER-3/PARASOL BRDF数据集中的地表方向反射率共计9024条，其中植被+缓冲1的数据为5873条，裸地+缓冲1+缓冲2为3660条，冰雪+缓冲2为750条。经过大气辐射传输模拟，在每个格网(angular bin)中获得MODIS大气层顶方向反射率及其对应的地表宽波段反照率数据10395648条，其中植被+缓冲1的数据为6765696条，裸地+缓冲1+缓冲2为4216320条，冰雪+缓冲2为864000条。

表4大气参数查找表

6S大气参数	参数设置
		大气类型	热带、中纬度夏季、中纬度冬季、副极地夏季、副极地冬季、US62标准大气
气溶胶类型	大陆型、海洋型、城市型、沙漠型、生物燃烧型、霾型
		气溶胶光学厚度	0.01，0.05，0.1，0.2
目标海拔高度	0，0.5，1.0，1.5，2.0，2.5，3，3.5(km)
		太阳天顶角	0，4，8，……，76，80(度)
观测天顶角	0，4，8，……，60，64(度)
		相对方位角	0，20，40，……，160，180(度)

步骤H：根据所述样本MODIS表观反射率数据和所述宽波段地表反照率数据建立第二转换模型，根据所述第二转换模型和经过网格划分的MODIS表观反射率数据，得到第二地表反照率数据。

在大气层顶表观反射率及其对应地表反照率数据集创建后，采用分格网回归分析的方法建立二者之间的经验关系。考虑到POLDER-3/PARASOL BRDF数据集中的地表方向反射率数据波段(b1-490nm，b2-565nm，b3-670nm，b4-765nm，b5-865nm，b6-1020nm)与MODIS前4个光学波段(b1-648nm，b2-859nm，b3-466nm，b4-554nm)在波长范围上较为一致，相对于MODIS后3个光学波段(b5-1244nm，b6-1631nm，b7-2119nm)波段转换的RMSE更小，转换精度更高。并且MODIS后3个光学波段容易受到大气中水汽吸收的影响，因此最终我们选择采用MODIS前4个波段(b1-648nm，b2-859nm，b3-466nm，b4-554nm)作为回归分析中大气层顶方向反射率的输入数据。利用大气辐射传输模型模拟获得的MODIS 4个波段(b1-648nm，b2-859nm，b3-466nm，b4-554nm)的大气层顶反射率与宽波段白空反照率、黑空反照率建立多元线性回归关系，如公式7、公式8所示。

其中α_ws为地表宽波段白空反照率，α_bs(θ_k)为地表宽波段黑空反照率，θ_k为黑空反照率的太阳天顶角，从0到80度，5度为步长共计17个间隔(θ_k＝0，5，…，75，80)，k＝1，2，3…，16，17，i为传感器波段(i＝1，2，3，4)，c_i为回归得到的回归系数，

为经过气体吸收透过率归一化公式得到的反射率，所述归一化公式如下：

其中

为大气层顶反射率，t_g为大气吸收透过率。

选取FluxNet，GC-Net观测站点结果对生成的第二地表反照率数据进行比较。图3是Bondville站点(40.0061°N，88.2919°W，农作物类型)2006年反照率地面测量结果与第二地表反照率数据的对比图。其中，圆点表示地面实测的反照率，叉点为第二地表反照率数据，从图中可以发现第二地表反照率数据与地面实测结果具有较好的一致性，在农作物播种、生长和收割的整个耕作周期内，包含了从植被、裸地以及冰雪类型的地表场景，在这一过程中第二地表反照率数据与地面站点测量结果比较一致，与地面观测结果对比RMSE为0.034。

图4是Flagstaff-Wildfire站点(35.446°N，111.727°W，草地类型)2006年反照率地面测量结果与第二地表反照率数据的对比图。其中，圆点表示地面实测的反照率，叉点为第二地表反照率数据，从图中可以发现第二地表反照率数据与地面实测结果在均值上非常接近(RMSE＝0.046)，但第二地表反照率数据较地面测试结果扰动较大，经过平滑处理后会得到更好的结果。从图中可以发现第二地表反照率数据对10天，80天和350天左右的短期降雪过程有较好的响应，在不被云完全覆盖的时段都能较好地反映降雪和融雪过程导致的地表反射率突然大幅度升高和降低的现象。

步骤I：整合所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。

所述第一地表反照率数据(以下简称AB1)、第二地表反照率数据(以下简称AB2)，再加上MODIS地表反照率产品(以下简称MOD43)，同一个地区具有3种地表反照率数据。本步骤的目的是整合这3种数据，生成时间上和空间上连续的地表反照率产品。

无论AB1、AB2还是MCD43，都存在个别数据缺失的现像，因此通过数据整合，可以填补缺失的数据，生成时间上和空间上连续的数据。

所述步骤I具体包括以下：

步骤I1：对所述第一地表反照率数据和所述第二地表反照率数据进行插值弥补。首先选用MCD43统计全球地表反照率时空分布的先验知识，称之为全球地表反照率背景场。考虑数据存储量以及MCD43的时间、空间代表性，最终确定背景场时间分辨率8天，空间分辨率5km。背景场中统计量有：全球逐像元的反照率均值、反照率方差以及同一像元不同时间的反照率之间的相关性。根据先验知识，判断AB1和AB2产品中对应像元的值是不是缺失，如果缺失，则使用前后时间像元值进行插值弥补，否则，无需插值。

步骤I2：对所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据进行归一化处理。AB1、AB2和MCD43之间可能存在系统偏差，在融合不同数据产品之前应该首先纠正它们之间的系统偏差。我们采用简单的归一化方法，先用线性变换消除不同数据产品之间的偏差，进行融合等处理之后再反变换回来。

步骤I3：对归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据进行数据融合，得到融合数据。

步骤I4：对所述融合数据进行平滑处理。

步骤I5：对平滑处理后的融合数据进行反归一化处理，得到所述时空连续的地表反照率产品。其中，反归一化处理的公式如下：

$\mathrm{\ξ}={\mathrm{\ξ}}^{\′}\stackrel{\‾}{\mathrm{\σ}}-\stackrel{\‾}{\mathrm{\μ}}---\left(10\right)$

其中，ξ是所述平滑处理后的融合数据，ξ是所述时空连续的地表反照率产品，

是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的均值，

是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的消除随机误差后的标准差的平均值。

图5是Willow_Creek站点不同方法生成的地表反照率产品整合后的结果对比示意图。如图5所示，从整合结果可以看出，通过本发明方法整合后的地表反照率数据可以弥补其他方法的缺点，得到长时间序列并且时空连续、高时空分辨率的地表反照率产品。

本发明实施例所述地表反照率产品的生成系统包括：

数据接收模块：用于接收POLDER BRDF数据、MODIS地表反射率数据、MODIS表观反射率数据、MODIS地表反照率数据和配套光谱数据。

网格划分模块：用于根据太阳天顶角、观测天顶角和相对方位角对MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据进行网格划分；

POLDER BRDF数据预处理模块：用于对POLDER BRDF数据进行优选，按照地表类型对经过优选的POLDER BRDF数据进行分类；

样本MODIS地表模块：用于根据配套光谱数据模拟得到模拟POLDER BRDF数据和模拟MODIS地表反射率数据，建立所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系，根据所述多元线性转换关系和经过分类的POLDERBRDF数据，计算得到样本MODIS地表反射率数据；

POLDER转换模块：用于根据经过分类的POLDER BRDF数据，得到POLDER反照率数据，转换所述POLDER反照率数据得到宽波段地表反照率数据；

第一产品模块：用于根据所述样本MODIS地表反射率数据和所述宽波段地表反照率数据建立第一转换模型，根据所述第一转换模型和经过网格划分的MODIS地表反射率数据，得到第一地表反照率数据；

样本MODIS表观模块：用于根据所述POLDER BRDF数据得到地表反射特性参数，根据所述地表反射特性参数和大气状态参数计算得到样本MODIS表观反射率数据；

第二产品模块：用于根据所述样本MODIS表观反射率数据和所述宽波段地表反照率数据建立第二转换模型，根据所述第二转换模型和经过网格划分的MODIS表观反射率数据，得到第二地表反照率数据；

终产品模块：用于整合所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。

本发明实施例所述地表反照率产品的生成方法及系统，以现有反照率产品为基础，通过数据融合技术，生产出长时间序列并且时空连续、高时空分辨率的地表反照率产品，对于全球陆面变化研究和陆面模型研发具有积极意义。

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

Claims (7)

1.一种地表反照率产品的生成方法，其特征在于，包括步骤：

A：输入POLDER BRDF数据、MODIS地表反射率数据、MODIS表观反射率数据、MODIS地表反照率数据和配套光谱数据；

其中，POLDER BRDF数据的空间分辨率为6km×7km，是多角度卫星遥感数据，有丰富的角度、光谱和极化信息，是多角度遥感的理想数据之一，所述POLDER BRDF数据共包括13227个数据文件；

MODIS地表反射率数据是指MOD09GA产品，MODIS地表反射率数据中含有MODIS的1-7波段的500米分辨率的反射率数据和1km分辨率的观测天顶角、观测方位角、太阳天顶角、太阳方位角和质量控制信息；

MODIS表观反射率数据是指MODIS大气层顶反射率数据，其中含有MODIS的36个波段的数据以及相应的辐射定标信息和地理定位信息；

MODIS地表反照率数据是利用半经验的线性核驱动模型及16天合成的多角度多波段观测数据反演得到的地表反照率产品，是Terra和Aqua星数据共同合成的MODIS Level3标准数据产品之一；

B：根据太阳天顶角、观测天顶角和相对方位角对MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据进行网格划分；

C：对POLDER BRDF数据进行优选，按照地表类型对经过优选的POLDER BRDF数据进行分类；

D：根据配套光谱数据模拟得到模拟POLDER BRDF数据和模拟MODIS地表反射率数据，建立所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系，根据所述多元线性转换关系和经过分类的POLDER BRDF数据，计算得到样本MODIS地表反射率数据；

E：根据经过分类的POLDER BRDF数据，得到POLDER反照率数据，转换所述POLDER反照率数据得到宽波段地表反照率数据；

F：根据所述样本MODIS地表反射率数据和所述宽波段地表反照率数据建立第一转换模型，根据所述第一转换模型和经过网格划分的MODIS地表反射率数据，得到第一地表反照率数据；

G:根据所述POLDER BRDF数据得到地表反射特性参数，根据所述地表反射特性参数和大气状态参数计算得到样本MODIS表观反射率数据；

H：根据所述样本MODIS表观反射率数据和所述宽波段地表反照率数据建立第二转换模型，根据所述第二转换模型和经过网格划分的MODIS表观反射率数据，得到第二地表反照率数据；

I：整合所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。

2.如权利要求1所述的方法，其特征在于，所述步骤I具体包括步骤：

I1：对所述第一地表反照率数据和所述第二地表反照率数据进行插值弥补；

I2：对所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据进行归一化处理；

I3：对归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据进行数据融合，得到融合数据；

I4：对所述融合数据进行平滑处理；

I5：对平滑处理后的融合数据进行反归一化处理，得到所述时空连续的地表反照率产品。

3.如权利要求2所述的方法，其特征在于，所述反归一化处理的公式如下：

$\mathrm{\ξ}={\mathrm{\ξ}}^{\′}\stackrel{\‾}{\mathrm{\σ}}-\stackrel{\‾}{\mathrm{\μ}}$

其中，ξ'是所述平滑处理后的融合数据，ξ是所述时空连续的地表反照率产品，是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的均值，

是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的消除随机误差后的标准差的平均值。

4.如权利要求1所述的方法，其特征在于，所述步骤C具体包括步骤：

C1：从所述POLDER BRDF数据中剔除由于云和气溶胶影响以及地表状态发生变化而不满足二向反射模型的数据；

C2：对经过所述步骤C1处理后的POLDER BRDF数据，按照植被、冰雪和裸地三种地表类型进行分类。

5.如权利要求4所述的方法，其特征在于，所述步骤C2具体包括步骤:

C21：判断POLDER BRDF数据的NDVI值是否大于0.2，如果是，则判定所述POLDER BRDF数据对应的像元为植被，否则执行步骤C22；

C22：判断POLDER BRDF数据是否其蓝光波段反射率大于0.3或者红光波段反射率大于0.3，如果是，则判定所述POLDER BRDF数据对应的像元为冰雪，否则，判定所述POLDER BRDF数据对应的像元为裸地。

6.如权利要求1所述的方法，其特征在于，所述步骤E具体包括步骤：

E1：对所述经过分类的POLDER BRDF数据进行半球积分，得到各波段的POLDER反照率数据；

E2：根据窄波段向宽波段的转换公式，转换所述POLDER反照率数据得到所述宽波段地表反照率数据。

7.一种地表反照率产品的生成系统，其特征在于，包括：

数据接收模块：用于接收POLDER BRDF数据、MODIS地表反射率数据、MODIS表观反射率数据、MODIS地表反照率数据和配套光谱数据；

其中，POLDER BRDF数据的空间分辨率为6km×7km，是多角度卫星遥感数据，有丰富的角度、光谱和极化信息，是多角度遥感的理想数据之一，所述POLDER BRDF数据共包括13227个数据文件；

MODIS地表反射率数据是指MOD09GA产品，MODIS地表反射率数据中含有MODIS的1-7波段的500米分辨率的反射率数据和1km分辨率的观测天顶角、观测方位角、太阳天顶角、太阳方位角和质量控制信息；

MODIS表观反射率数据是指MODIS大气层顶反射率数据，其中含有MODIS的36个波段的数据以及相应的辐射定标信息和地理定位信息；

MODIS地表反照率数据是利用半经验的线性核驱动模型及16天合成的多角度多波段观测数据反演得到的地表反照率产品，是Terra和Aqua星数据共同合成的MODIS Level3标准数据产品之一；

网格划分模块：用于根据太阳天顶角、观测天顶角和相对方位角对MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据进行网格划分；

POLDER BRDF数据预处理模块：用于对POLDER BRDF数据进行优选，按照地表类型对经过优选的POLDER BRDF数据进行分类；

样本MODIS地表模块：用于根据配套光谱数据模拟得到模拟POLDER BRDF数据和模拟MODIS地表反射率数据，建立所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系，根据所述多元线性转换关系和经过分类的POLDERBRDF数据，计算得到样本MODIS地表反射率数据；

POLDER转换模块：用于根据经过分类的POLDER BRDF数据，得到POLDER反照率数据，转换所述POLDER反照率数据得到宽波段地表反照率数据；

第一产品模块：用于根据所述样本MODIS地表反射率数据和所述宽波段地表反照率数据建立第一转换模型，根据所述第一转换模型和经过网格划分的MODIS地表反射率数据，得到第一地表反照率数据；

样本MODIS表观模块：用于根据所述POLDER BRDF数据得到地表反射特性参数，根据所述地表反射特性参数和大气状态参数计算得到样本MODIS表观反射率数据；

第二产品模块：用于根据所述样本MODIS表观反射率数据和所述宽波段地表反照率数据建立第二转换模型，根据所述第二转换模型和经过网格划分的MODIS表观反射率数据，得到第二地表反照率数据；

终产品模块：用于整合所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。

Images