专利名称：地表反照率产品的生成方法及系统的制作方法
技术领域：
本发明涉及卫星遥感技术领域，特别涉及一种地表反照率产品的生成方法及系统。
背景技术：
地表反照率反映了地球表面对太阳光辐射的反射能力，是地表辐射能量平衡以及地气相互作用中的驱动因子之一，其时空变化受到自然过程以及人类活动的影响，是全球环境变化的指示因子。目前一些卫星反照率产品已经业务化运行和发布，空间分辨率从250m_20km，时间分辨率从日到月，其中极轨卫星产品包括MODIS，MISR，CERES, POLDER, MERIS0地球同步卫星产品包括 Meteosat, and MSG。目前针对在轨运行的卫星遥感数据以及部分航空遥感数据进行了大量地表反照率反演研究，并进行了有效的反照率产品的地面验证，其中基于线性核驱动模型的地表 BRDF/反照率遥感反演模型算法是目前地表反照率遥感反演中应用最广泛的方法，已经在 M0DIS、MISR的地表反照率产品中得到了很好的应用。另外，POLDER/PARASOL系列传感器具有更好的多角度观测能力，空间分辨率稍低(6km)，也发布了非常有特色的全球长时间序列反照率产品。其他如静止轨道气象卫星MSG、METE0SET，极轨卫星传感器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具体包括步骤Il 对所述第一地表反照率数据和所述第二地表反照率数据进行插值弥补；12:对所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的 MODIS地表反照率数据进行归一化处理；13 对归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据进行数据融合，得到融合数据；14 对所述融合数据进行平滑处理；15 对平滑处理后的融合数据进行反归一化处理，得到所述时空连续的地表反照
率产品。优选地，所述反归一化处理的公式如下ξ = ξσ-^ι其中，ξ是所述平滑处理后的融合数据，ξ是所述时空连续的地表反照率产品，；; 是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的均值，G是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的消除随机误差后的标准差的平均值。优选地，所述步骤C具体包括步骤Cl 从所述POLDER BRDF数据中剔除由于云和气溶胶影响以及地表状态发生变化而不满足二向反射模型的数据；C2 对经过所述步骤Cl处理后的POLDER BRDF数据，按照植被、冰雪和裸地三种地
表类型进行分类。 优选地，所述步骤C2具体包括步骤
C21 判断POLDER BRDF数据的NDVI值是否大于0. 2，如果是，则判定所述POLDER BRDF数据对应的像元为植被，否则执行步骤C22 ；C22 判断POLDER BRDF数据是否其蓝光波段反射率大于0. 3或者红光波段反射率大于0. 3，如果是，则判定所述POLDER BRDF数据对应的像元为冰雪，否则，判定所述POLDER BRDF数据对应的像元为裸地。优选地，所述步骤E具体包括步骤El 对所述经过分类的POLDER BRDF数据进行半球积分，得到各波段的POLDER反照率数据；E2:根据窄波段向宽波段的转换公式，转换所述POLDER反照率数据得到所述宽波段地表反照率数据。本发明还提供一种地表反照率产品的生成系统，其包括网格划分模块用于根据太阳天顶角、观测天顶角和相对方位角对MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据进行网格划分；POLDER BRDF数据预处理模块用于对POLDER BRDF数据进行优选，按照地表类型对经过优选的POLDER BRDF数据进行分类；样本MODIS地表模块用于根据配套光谱数据模拟得到模拟POLDER BRDF数据和模拟MODIS地表反射率数据，建立所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系，根据所述多元线性转换关系和经过分类的P0LDERBRDF数据，计算得到样本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数据的空间分辨率为(6kmX 7km)，是目前能获得的最新多角度卫星遥感数据，有丰富的角度、光谱和极化信息，是多角度遥感的理想数据之一。输入的POLDER BRDF数据共包括13227个数据文件。(2)MODIS地表反射率数据MODIS地表反射率数据是指M0D09GA产品。MODIS地表反射率数据中含有MODIS 的1-7波段的500米分辨率的反射率数据和Ikm分辨率的观测天顶角、观测方位角、太阳天顶角、太阳方位角和质量控制等信息。(3)MODIS表观反射率数据MODIS表观反射率数据是指MODIS大气层顶反射率数据，即MODIS的LlB数据，其中含有MODIS的36个波段的数据以及相应的辐射定标信息和地理定位信息。(4)配套光谱数据配套的光谱数据共有493个，分别是梁顺林教授2004年著作《定量遥感》所附光盘中提供的119条波谱、“我国典型地物标准波普数据库”提供的2M条植被和土壤波谱、黑河综合遥感联合实验采集的103条典型地物波谱和格林兰采集的47条冰雪波谱数据。(5)MODIS地表反照率数据MODIS地表反照率数据(MCD43)是利用半经验的线性核驱动模型及16天合成的多角度多波段观测数据反演得到的地表反照率产品。是Terra和Aqua星数据共同合成的MODIS Level3标准数据产品之一。以MCD43B3产品为例，它具有Ikm的空间分辨率以及8天的时间分辨率，其数据产品中包含MODIS传感器1-7波段的7个窄波段以及可见光 (0. 3-0. 7 μ m)、近红外(0. 7-5. Oym)及短波波段(0. 3-5. Oym)三个宽波段的白空反照率和黑空反照率。步骤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具体包括步骤Cl 从所述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个，剩余90 个数据文件成为优选的POLDER BRDF数据。步骤C2 对经过所述步骤Cl处理后的POLDER BRDF数据，按照植被、冰雪和裸地
三种地表类型进行分类。不同地表类型具有不同的双向反射率特征，有必要引入地物分类信息，进一步细分训练样本，减少后续步骤中线性回归模型的不确定性。因此，我们选择直接根据遥感观测数据分类的策略，分类的具体步骤如下步骤C21 计算 POLDER BRDF 数据的 NDVI (Normalized Difference Vegetation Index,归一化植被指数)。判断NDVI是否大于0. 2，如果是，则判定所述POLDER BRDF数据对应的像元为植被，否则执行步骤C22。其中，NDVI的具体计算公式为NDVI = (R2-R1)/(R2+R1) (1)其中R2表示地物0. 865微米的地表反射率值，Rl表示地物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地表反射率数据表
权利要求
1.一种地表反照率产品的生成方法，其特征在于，包括步骤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所述的方法，其特征在于，在所述步骤B之前还包括步骤A输入 POLDER BRDF数据、MODIS地表反射率数据、MODIS表观反射率数据、MODIS地表反照率数据和配套光谱数据。
3.如权利要求1所述的方法，其特征在于，所述步骤I具体包括步骤II对所述第一地表反照率数据和所述第二地表反照率数据进行插值弥补；12 对所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS 地表反照率数据进行归一化处理；13:对归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的 MODIS地表反照率数据进行数据融合，得到融合数据； 14 对所述融合数据进行平滑处理；15:对平滑处理后的融合数据进行反归一化处理，得到所述时空连续的地表反照率产PΡΠ O
4.如权利要求3所述的方法，其特征在于，所述反归一化处理的公式如下 ξ=ξσ-μ其中，ξ是所述平滑处理后的融合数据，ξ是所述时空连续的地表反照率产品，G是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的均值，G是归一化处理后的第一地表反照率数据、第二地表反照率数据和经过网格划分的MODIS地表反照率数据的消除随机误差后的标准差的平均值。
5.如权利要求1所述的方法，其特征在于，所述步骤C具体包括步骤Cl 从所述POLDER BRDF数据中剔除由于云和气溶胶影响以及地表状态发生变化而不满足二向反射模型的数据；C2 对经过所述步骤Cl处理后的POLDER BRDF数据，按照植被、冰雪和裸地三种地表类型进行分类。
6.如权利要求5所述的方法，其特征在于，所述步骤C2具体包括步骤C21 判断POLDER BRDF数据的NDVI值是否大于0. 2，如果是，则判定所述POLDER BRDF 数据对应的像元为植被，否则执行步骤C22 ；C22 判断POLDER BRDF数据是否其蓝光波段反射率大于0. 3或者红光波段反射率大于 0. 3，如果是，则判定所述POLDER BRDF数据对应的像元为冰雪，否则，判定所述POLDER BRDF 数据对应的像元为裸地。
7.如权利要求1所述的方法，其特征在于，所述步骤E具体包括步骤El 对所述经过分类的POLDER BRDF数据进行半球积分，得到各波段的POLDER反照率数据；E2 根据窄波段向宽波段的转换公式，转换所述POLDER反照率数据得到所述宽波段地表反照率数据。
8.—种地表反照率产品的生成系统，其特征在于，包括网格划分模块用于根据太阳天顶角、观测天顶角和相对方位角对MODIS地表反射率数据、MODIS表观反射率数据和MODIS地表反照率数据进行网格划分；POLDER BRDF数据预处理模块用于对POLDER BRDF数据进行优选，按照地表类型对经过优选的POLDER BRDF数据进行分类；样本MODIS地表模块用于根据配套光谱数据模拟得到模拟POLDER BRDF数据和模拟 MODIS地表反射率数据，建立所述模拟POLDER BRDF数据和所述模拟MODIS地表反射率数据之间的多元线性转换关系，根据所述多元线性转换关系和经过分类的P0LDERBRDF数据，计算得到样本MODIS地表反射率数据；POLDER转换模块用于根据经过分类的POLDER BRDF数据，得到POLDER反照率数据， 转换所述POLDER反照率数据得到宽波段地表反照率数据；第一产品模块用于根据所述样本MODIS地表反射率数据和所述宽波段地表反照率数据建立第一转换模型，根据所述第一转换模型和经过网格划分的MODIS地表反射率数据， 得到第一地表反照率数据；样本MODIS表观模块用于根据所述POLDER BRDF数据得到地表反射特性参数，根据所述地表反射特性参数和大气状态参数计算得到样本MODIS表观反射率数据；第二产品模块用于根据所述样本MODIS表观反射率数据和所述宽波段地表反照率数据建立第二转换模型，根据所述第二转换模型和经过网格划分的MODIS表观反射率数据， 得到第二地表反照率数据；终产品模块用于整合所述第一地表反照率数据、所述第二地表反照率数据和经过网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。
9.如权利要求8所述的系统，其特征在于，所述系统还包括数据接收模块用于接收POLDER BRDF数据、MODIS地表反射率数据、MODIS表观反射率数据、MODIS地表反照率数据和配套光谱数据。
全文摘要
本发明公开了一种地表反照率产品的生成方法及系统，涉及卫星遥感领域。所述方法包括步骤根据经过分类的POLDER BRDF数据，得到POLDER反照率数据，转换POLDER反照率数据得到宽波段地表反照率数据；根据样本MODIS地表反射率数据和宽波段地表反照率数据建立第一转换模型，得到第一地表反照率数据；根据样本MODIS表观反射率数据和宽波段地表反照率数据建立第二转换模型，得到第二地表反照率数据；整合第一和第二地表反照率数据，及经网格划分的MODIS地表反照率数据，生成时空连续的地表反照率产品。本发明生产出了长时间序列并且时空连续、高时空分辨率的地表反照率产品，提高了地表反照率产品反演的精度。
文档编号G01N21/55GK102435586SQ20111027610
公开日2012年5月2日 申请日期2011年9月16日 优先权日2011年9月16日
发明者刘强, 张鑫, 梁顺林, 王立钊, 瞿瑛, 赵祥 申请人:北京师范大学