CN102901563B - 一种同时确定地表窄波段和宽波段比辐射率的方法及装置 - Google Patents

Abstract

一种同时确定地表窄波段和宽波段比辐射率的方法及装置，包括：A）对中分辨率成像光谱仪（MODIS）数据进行数据处理，包括辐射定标、大气纠正和云掩膜处理，获取晴空条件下MODIS数据位于0.4-2.1μm大气窗口区的可见光/近红外波段地表反射率数据；B）利用步骤A的反射率数据，结合发展的窄波段比辐射率反演算法，确定MODIS热红外光谱谱段第29（波谱范围8.4-8.7μm）、31（波谱范围10.78-11.28μm）和32（波谱范围11.77-12.27μm）波段的窄波段地表比辐射率；C）利用步骤B的窄波段地表比辐射率，结合发展的窄-宽波段比辐射率转换模型，得到3-14μm以及3-∞μm的宽波段地表比辐射率。

Description

一种同时确定地表窄波段和宽波段比辐射率的方法及装置

技术领域

本发明属遥感定量反演领域，特别涉及一种同时确定地表窄波段和宽波段比辐射率的方法及装置。

背景技术

地表比辐射率作为地表自然物质的固有属性，受土地覆盖类型，土壤水分，土壤有机质，植被浓密程度和结构、表面粗糙度和地表起伏等因素的影响，能够指示地表的物质组成，在地表分类、土壤形成和侵蚀，稀疏植被覆盖与变化估算，岩床制图和资源探测等方面都有重要的意义。除此之外，地表比辐射率在地表温度的反演中扮演着重要的角色，没有准确的比辐射率就很难得到精确的温度反演结果。当地表温度反演精度向着1K，甚至更高精度迈进时，地表比辐射率的反演水平也需要向前挺进。

地表比辐射率的反演是科学家公认的一个病态问题。从普朗克方程可以看出，某个波段接收到的辐射能是温度和比辐射率的函数。因此，N个波段观测的辐亮度，总有N+1个未知数（N个比辐射率和1个温度），温度和比辐射率始终耦合在一起，其中任何一个物理量的确定需要以另一个物理量的确定为前提。这种病态性成为了地表比辐射率反演的难点之一。虽然现有MODIS地表比辐射率产品已经业务化生产，但是该产品也仅仅是在大气窗口区间上6个波段在卫星观测方向上的窄波段地表比辐射率，而不是各种水热平衡模型和陆面过程模型等所要求的地表宽波段比辐射率。迄今为止，还没有一种切实可行的算法能够为气候系统模型提供地表宽波段比辐射率数据，这也是为什么目前各种气候系统模型都将地表宽波段比辐射率值设为1或接近于1的常数的原因。由于地表比辐射率决定了地表长波辐射大小，如果简单采用窄波段比辐射率替代宽波段比辐射率，长波辐射的反演误差可能达到100W/m²。由此可见，如何实现窄波段地表比辐射率向宽波段地表比辐射率的转换是目前遥感反演需要解决的难点和热点。

由于受测量仪器的限制，现有的地物波谱库中一般只测量了14μm以下的地物波谱，这在一定程度上阻碍了地表宽波段比辐射率的估算。为了实时高效的从卫星数据中直接获取地表3-14μm以及3-∞μm的宽波段比辐射率，解决地表比辐射率的反演技术本身涉及的病态核心问题，发展一种地表宽波段比辐射率估算方法和装置是本发明提出的初衷。

发明内容

本发明的技术解决问题：克服现有技术的不足，提供一种同时确定地表窄波段和宽波段比辐射率的方法及装置。

本发明的技术解决方案：地表窄波段和宽波段比辐射率确定方法，实现步骤如下：

（A）．对中分辨率成像光谱仪（MODIS）数据进行数据处理，包括辐射定标、大气纠正、以及云掩膜等处理，获取晴空条件下MODIS数据位于0.4-2.1μm大气窗口区的可见光/近红外波段地表反射率数据；

（B）．利用步骤A的反射率数据，结合建立的窄波段地表比辐射率反演算法，确定MODIS热红外光谱谱段第29（波谱范围8.4-8.7μm）、31（波谱范围10.78-11.28μm）和32（波谱范围11.77-12.27μm）波段的窄波段地表比辐射率；

（C）．利用步骤B的窄波段地表比辐射率，结合建立的窄-宽波段比辐射率转换模型，得到地表3-14μm以及3-∞μm的宽波段比辐射率。

所述步骤（A）中对中分辨率成像光谱仪（MODIS）数据进行数据处理的过程是：

（A1）根据MODIS提供的定标系数，对MODIS遥感数据进行绝对辐射定标，将影像灰度值转换为辐射亮度；

（A2）利用快速大气纠正模型（QUAC）和MODIS的云掩膜产品MOD35数据，完成MODIS遥感数据的大气纠正和有云像元的掩膜处理，获取MODIS数据位于0.4-2.1μm大气窗口区的可见光/近红外谱段的地表反射率数据。

所述步骤（B）中确定MODIS数据热红外光谱谱段29、31和32波段的窄波段地表比辐射率过程是：根据计算出的可见光/近红外地表反射率数据，计算归一化的积雪指数（NDSI）和归一化的植被指数（NDVI），并结合MODIS的陆地/海洋掩膜数据，对像元进行分类，分别分为水体像元，雪/冰像元和植被/裸土区像元三大类，分别得到所述比辐射率；

（B1）对于水体像元，利用ASTER波谱库和UCSB波谱库中测量的一些典型的水体波谱曲线，结合MODIS的第29、31和32波段的光谱响应函数分别进行卷积处理，获得每个波段上比辐射率的平均值为0.984、0.992和0.987，最后将该三个值分别作为MODIS热红外光谱谱段29、31和32波段的比辐射率；

（B3）对于雪/冰像元，利用ASTER波谱库和UCSB波谱库中测量的一些典型的雪/冰波谱曲线，结合MODIS的第29、31和32波段的光谱响应函数分别进行卷积处理，获得每个波段上比辐射率的平均值，最后将该三个值分别作为雪/冰情况下MODIS热红外光谱谱段29、31和32波段的比辐射率。

（B4）对于植被/裸土像元，根据计算出的NDVI值大小，进一步细分为裸土像元、植被像元和植被/裸土混合像元。

针对裸土像元，现有的方法主要是将热红外波段的比辐射率和红光波段的反射率之间建立一种简单的一元线性回归模型。可通过分析ASTER波谱库中不同裸土的波谱曲线发现，通过这样的简单一元回归模型获取的热红外波段比辐射率的精度并不高，误差会高达20%以上。本发明通过将MODIS热红外波段比辐射率与可见光/近红外波段多个波段的反射率建立多元线性回归模型，如公式（1），精度提高10%以上。

ε_si＝a_i+∑b_ijρ_j（i=29,31,32;j=1-7）              （1）

公式（1）中，ε_si表示为热红外第i波段裸土的比辐射率，ρ_j表示为是可见光/近红外第j波段的反射率，a_i和b_ij是模型转换系数，通过多元回归获取。

针对植被像元，现有的方法主要是通过赋常数值的方法来获取各个热红外波段的比辐射率，这种方法没有反映出植被的真正覆盖及分布结构信息。为了体现不同观测角度情况下植被的三维结构对热红外比辐射率的影响，本发明通过建立热红外比辐射率与NDVI的关系模型，来获取热红外光谱谱段的比辐射率值：

ε_vi＝c_0i+c_1iNDVI(i=29,31,32)               （2）

公式（2）中，ε_vi表示为热红外第i波段浓密植被的比辐射率，c_0i和c_1i是模型转换系数。

这样本发明通过上述获得的比辐率值能够反映出植被的真正覆盖及分布结构信息。

针对植被/裸土混合像元，通过建立热红外波段比辐射率与植被覆盖度P_v的关系模型，来获取热红外光谱谱段的比辐射率：

ε_i＝d_0i+d_1iP_v        （3）

公式（3）中，ε_i表示为热红外第i波段在植被/裸土混合区的比辐射率，d_0i和d_1i是模型转换系数，P_v表示为植被覆盖度，其值可以通过NDVI值来确定：

$P_v={\left[\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right]}^2---\left(4\right)$

式中，NDVI_min和NDVI_max分别为裸土和浓密植被的NDVI值。NDVI_mi_n和NDVI_max值可以从NDVI直方图获取。对于全球条件下，设定NDVI_min＝0.2和NDVI_max＝0.5。为了获取连续的P_v值，对于NDVI<NDVI_min的像元，设定P_v＝0；对于NDVI>NDVI_max的像元，设定P_v＝1。

所述步骤（C）中确定地表3-14μm以及3-∞μm的宽波段比辐射率过程是：

（C1）根据计算出的窄波段地表比辐射率，结合建立的窄-宽波段比辐射率转换模型（公式（5）），确定3-14μm地表宽波段比辐射率ε_3→14：

ε_3-14＝w₀+w₁*ε₂₉+w₂*ε₃₁+w₃*ε₃₂             (5)

式中，ε₂₉、ε₃₁和ε₃₂分别是MODIS热红外光谱谱段第29、31和32波段的窄波段比辐射率，w_i是模型转换系数。

（C2）根据计算出的窄波段地表比辐射率，结合建立的窄-宽波段比辐射率转换模型（公式（6）），确定3-∞μm地表宽波段比辐射率ε_3→∞：

ε_3-∞＝w′₀+w′₁*ε₂₉+w′₂*ε₃₁+w′₃*ε₃₂         （6）

式中，ε₂₉、ε₃₁和ε₃₂分别是MODIS热红外光谱谱段第29、31和32波段的窄波段比辐射率，w′_i是模型转换系数。

本发明提供的实现上述方法的装置，其包括：

数据预处理模块：对中分辨率成像光谱仪（MODIS）数据进行数据处理，包括辐射定标、大气纠正、以及云掩膜处理，获取晴空条件下MODIS数据位于0.4-2.1μm大气窗口区的可见光/近红外地表反射率数据；

窄波段比辐射率反演模块：利用数据预处理模块的反射率数据，结合窄波段地表比辐射率反演算法，确定MODIS热红外光谱谱段第29,即波谱范围8.4-8.7μm、第31，即波谱范围10.78-11.28μm和第32，即波谱范围11.77-12.27μm波段的窄波段地表比辐射率；具体为根据计算出的可见光/近红外地表反射率数据，计算归一化的积雪指数（NDSI）和归一化的植被指数（NDVI），并结合MODIS的陆地/海洋掩膜数据，对像元进行分类，分别分为水体像元，雪/冰像元和植被/裸土像元三大类；再结合获取的可见光/近红外地表反射率数据计算的归一化植被指数（NDVI）值的大小，进一步将植被/裸土像元分为裸土像元、植被像元和裸土/植被混合像元三类；针对每一类像元，分别利用自己建立的窄波段比辐射率反演算法，确定不同类别像元的MODIS热红外谱段29、31和32波段的窄波段比辐射率数据；

宽波段比辐射率反演模块：分别利用自己建立的3-14μm窄-宽波段地表比辐射率转换模型和3-∞μm窄-宽波段地表比辐射率转换模型，计算出每个像元3-14μm的地表宽波段比辐射率数据和3-∞μm的地表宽波段比辐射率数据。

本发明与现有技术相比的优点在于：

（1）通过本发明的步骤实现了同时确定地表窄波段和宽波段比辐射率的技术，具有实时和高效的优点。

（2）本发明提出的窄波段地表比辐射率确定方法中，针对裸土区，建立的MODIS热红外波段比辐射率与可见光/近红外波段反射率的多元线性回归模型，结果优于现有方法中利用热红外波段比辐射率仅与红光波段反射率建立线性关系模型的结果，精度提高10%以上；针对植被区，建立的MODIS热红外波段比辐射率与NDVI的线性回归模型，不仅克服了现有方法中对所有植被只赋一个常数值的缺陷，还保留了观测角度的信息。

（3）由于测量仪器性能的限制以及大气水汽对热红外波谱强烈吸收的影响，地面测量的光谱数据都局限在14μm以下，而且遥感传感器只在大气窗口区设置了几个窄的光谱波段，只能获取离散的窄波段比辐射率数据。因而，获取宽波段比辐射率数据（特别是3-∞μm）变得非常困难。目前国际上还没有现有的模型来反演3-∞μm宽波段地表比辐射率数据，而各种水热平衡模型和陆面过程模型等都需要3-∞μm宽波段的地表比辐射率数据作为模型输入。本发明建立的窄-宽波段地表比辐射率转换模型首次实现了3-∞μm的宽波段地表比辐射率反演，能满足现有各种水热平衡模型和陆面过程模型的需求。

（4）本发明建立的地表宽波段比辐射率确定装置是通过数据预处理模块、窄波段比辐射率反演模块和宽波段比辐射率反演模块来实现从MODIS原始数据中直接反演3-14μm以及3-∞μm宽波段地表比辐射率，具有操作简单、实用性强、及可扩展性强的特点。

附图说明

图1是本发明确定地表宽波段比辐射率的流程示意图；

图2是本发明快速大气校正处理流程示意图；

图3是本发明植被和裸土混合区地表窄波段比辐射率确定的流程示意图；

图4是本发明采用的窄-宽波段比辐射率转换模型的流程示意图。

具体实施方式

如图1所示，本发明具体实现如下：

A，对中分辨率成像光谱仪（MODIS）数据进行处理，包括辐射定标、大气纠正、以及云掩膜等，获取MODIS数据位于0.4-2.1μm大气窗口区的可见光/近红外波段地表反射率数据。该步骤主要通过数据预处理模块来实现，具体实施方式如下：

（1）绝对辐射定标

MODIS1B存储的数据为影像灰度值，需要经过辐射定标将影像灰度值转换为辐射亮度：

L_sensor＝gain*DN+offset          （7）

式中，L_sensor为传感器观测的辐射亮度，单位为Wm^-2sr^-1um^-1；DN为影像灰度值；gain和offset分别为增益值和偏移值，可以从MODIS数据中直接获取。

（2）大气校正

快速大气校正（Quick Atmospheric Correction Algorithm-QUAC）是基于辐射传输的半经验大气校正算法，它能够利用很少的参数快速地对多光谱和高光谱数据(0.4~2.5μm)做大气校正处理获取地表反射率。如图2所示，QUAC首先根据波长位置自动从多光谱或者高光谱数据上收集不同物质的波谱信息，然后将暗目标选定为基准并确定基准光谱及反射率偏移量。接着对采集的波谱信息进行移除反射率偏移量的差值处理。通过基准波谱匹配，确定地表光谱端元。通过端元波谱匹配，获取端元平均光谱。在端元平均参考光谱的支持下，获取反射率增益值。最后利用反射率增益值和偏移值完成高光谱和多光谱的快速大气校正，获取地表的反射率光谱。相对于其他大气校正算法，QUAC不依赖先验知识，且大气校正速度快，适用于多种遥感卫星数据。

B，利用步骤A的反射率数据，结合建立的窄波段地表比辐射率反演算法，确定MODIS热红外光谱谱段第29、31和32波段的窄波段地表比辐射率。该步骤主要通过窄波段比辐射率反演模块来实现，具体实施方式如下：

（1）MODIS图像的像元分类

首先结合MODIS数据产品MOD03中陆地/海洋掩膜数据，将图像中的像元分为海洋像元、内陆水体像元和陆地像元三大类。然后，再对图像中的内陆水体和陆地像元进一步进行细分类。具体分类方法为：根据计算出的可见光/近红外地表反射率数据，计算归一化的积雪指数（NDSI）和归一化的植被指数（NDVI），将内陆水体像元和陆地像元进一步分为内陆水体像元、雪/冰像元、裸土像元、植被像元以及植被/裸土混合像元五类。具体实施方式如下：

建立归一化积雪指数模型（NDSI）：

$\mathrm{NDSI}=\frac{{\mathrm{\&rho;}}_4{-\mathrm{\&rho;}}_6}{{\mathrm{\&rho;}}_4+{\mathrm{\&rho;}}_6}---\left(8\right)$

式中，ρ₄和ρ₆分别是MODIS第4波段，即波谱范围0.545-0.565μm和第6波段，即波谱范围1.628-1.652μm的地表反射率数据；计算每个像元的NDSI值，当像元的NDSI值大于0.4时，则判别该像元为雪/冰像元，否则判别为非雪/冰像元，即植被/裸土像元；再进一步对植被/裸土像元进行分类，利用MODIS的第1波段，即波谱范围0.620-0.670μm和第2波段，即波谱范围0.841-0.876μm的地表反射率数据，建立归一化的植被指数模型（NDVI）：

$\mathrm{NDVI}=\frac{{\mathrm{\&rho;}}_2{-\mathrm{\&rho;}}_1}{{\mathrm{\&rho;}}_2+{\mathrm{\&rho;}}_1}---\left(9\right)$

式中，ρ₁和ρ₂分别是MODIS第1波段和第2波段的地表反射率，计算整景图像中每个像元的NDVI值，当NDVI<0.2时，则判别该像元为裸土像元；当0.2≤NDVI≤0.5时，则判别为裸土/植被混合像元；当NDVI>0.5时，则判别该像元为植被像元。这样，对于整景MODIS图像，就可根据上述分类方法将像元分为不同的类型。

（2）MODIS像元尺度红外窄波段比辐射率的确定

由于水体、雪/冰的地表窄波段比辐射率通常较大，且光谱变异性较小，所以水体和雪/冰不同样本的波段比辐射率的均值即可以代表MODIS像元尺度红外窄波段比辐射率。因此水体29、31和32波段的窄波段比辐射率可以近似估算为0.984、0.992和0.987。相应的雪/冰29、31和32波段的窄波段比辐射率可以近似估算为0.988、0.987和0.966。

对于植被/裸土下垫面，由于植被长势以及植被与裸土比辐射率迥异的差别导致不能简单利用不同样本的波段比辐射率均值来代表MODIS像元尺度红外窄波段比辐射率。较为合理的方法是利用组分比例加权纯植被和纯裸土的比辐射率。考虑到归一化植被指数NDVI能够反映植被的长势即组分比例变化信息，因此基于NDVI值，可以区分不同的地表长势情况，如裸土，植被、裸土/植被的混合。针对不同的长势情况，MODIS像元尺度红外窄波段比辐射率的确定方法略有不同，如图3所示。

（i）裸土像元

对ASTER波谱库（此为美国宇航局公开的数据库，且为本领域技术人员常采用）中的裸土比辐射率和红光波段的反射率值比较分析发现，两者之间很难建立一种线性关系，而MODIS 29，31和32波段的比辐射率与MODIS 1—7波段的反射率之间则存在较好的线性关系。通过分析，建立了如下利用可见光/近红外反射率来估算比辐射率的多元线性模型：

ε_si＝a_i+∑b_ijρ_j          （10）

其中系数a和b如下表所示：

表1.裸土像元MODIS 29,31和32波段比辐射率线性模型系数和统计结果

（ii）植被和裸土混合像元

ε_i＝ε_v,iP_v+ε_s,i(1-P_v)+C_i          （11）

式中，ε_v,i和ε_s,i分别为第i波段的植被和裸土比辐射率；P_v为植被覆盖度；C_i为取决于地表特性的并考虑腔体效应的增加项。对于水平均一地表，C_i≈0；对于粗糙非均一地表，C_i通过下式获得：

C_i＝(1-ε_s,i)ε_v,iF(1-P_v)            （12）

式中，F为取决于地表几何分布的形状因子，其值范围为0-1。由于从遥感可见光-近红外或者热红外数据不能估算F，通常假定不同的几何分布取其平均值，即0.55。

P_v可以根据NVDI值由下式获得：

$P_v={\left[\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right]}^2---\left(13\right)$

式中，NDVI_min和NDVI_max分别为裸土和浓密植被的NDVI值。NDVI_min和NDVI_max值可以从NDVI直方图获取。对于全球条件下，设定NDVI_min＝0.2和NDVI_max＝0.5。为了获取连续的P_v值，对于NDVI<NDVI_min的像元，设定P_v＝0；对于NDVI>NDVI_max的像元，设定P_v＝1。按照以上的分析，MODIS 29，31和32波段的比辐射率与植被覆盖度（Pv）之间的关系可以表示为：

ε₂₉＝0.971+0.010P_v

ε₃₁＝0.987-0.005P_v          （14）

ε₃₂＝0.990-0.006P_v

（iii）植被像元

对于植被覆盖区像元，即NDVI>NDVI_max时，通常可以将比辐射率设定为一个常数（即0.99）。考虑到NDVI能够在一定程度上反映出植物长势以及观测方向，因此可以建立比辐射率与NDVI的关系模型，结合MODIS 29，31和32波段的光谱响应函数，对植被覆盖区的比辐射率和NDVI值进行了分析，建立了如下的关系模型：

ε_vi＝c_0i+c_1iNDVI           （15）

其中系数c₀和c₁如下表所示：

表2.植被像元MODIS 29,31和32波段比辐射率线性模型系数和统计结果

C，利用步骤B的窄波段地表比辐射率，结合建立的窄-宽波段比辐射率转换模型，得到3-14μm以及3-∞μm的宽波段地表比辐射率。该步骤主要通过宽波段比辐射率反演模块来实现，具体实施方式如下：

（1）3-14μm宽波段地表比辐射率

地表比辐射率被定义为地表向外发射的热辐射与同温度黑体地表向外发射的热辐射的比率。MODIS传感器为了增加观测的信噪比，波段被设计成具有一定宽度的光谱谱段，因此MODIS红外波段比辐射率被定义为：

${\mathrm{\&epsiv;}}_i=\frac{\underset{\mathrm{\&lambda;i}1}{\overset{\mathrm{\&lambda;i}2}{\&Integral;}}{f}_i\left(\mathrm{\&lambda;}\right)\mathrm{\&epsiv;}\left(\mathrm{\&lambda;}\right)B(\mathrm{\&lambda;},T)\mathrm{d\&lambda;}}{\underset{\mathrm{\&lambda;i}1}{\overset{\mathrm{\&lambda;i}2}{\&Integral;}}{f}_i\left(\mathrm{\&lambda;}\right)B(\mathrm{\&lambda;},T)\mathrm{d\&lambda;}}---\left(16\right)$

其中f_i是红外波段的光谱响应函数，ε(λ)是比辐射率曲线，B(λ,T)是对应波长λ和温度T时的辐射亮度。

基于斯蒂芬-波尔兹曼（Stefan-Boltzmann）定律，对于一个自然物体，它的总发射辐射通量φ可以通过普朗克函数在整个电磁波谱区间积分而获得，即：

$\mathrm{\&phi;}=\underset{0}{\overset{\&infin;}{\&Integral;}}{\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}={\mathrm{\&epsiv;}}_b\mathrm{\&sigma;}{T}^4---\left(17\right)$

式中，ε_λ是在波长为λ处，物体的发射辐射与黑体的发射辐射B_λ(T)的一个比值，T是温度（K），ε_b是整个电磁波谱区间的比辐射率，σ是Stefan-Boltzmann’s常数(5.67×10^-8Wm^-2K^-4)。根据上式，可以推导出对应波长区间内的宽波段比辐射率为：

${\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}1\&RightArrow;\mathrm{\&lambda;}2}=\frac{\underset{\mathrm{\&lambda;}1}{\overset{\mathrm{\&lambda;}2}{\&Integral;}}{\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}}{\underset{\mathrm{\&lambda;}1}{\overset{\mathrm{\&lambda;}2}{\&Integral;}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}},---\left(18\right)$

式中，λ₁和λ₂分别对应于波谱区间的下边界和上边界。对于3-14μm电磁波谱段而言，λ₁等于3μm，而λ₂等于14μm。

结合窄波段比辐射率的定义，即方程式（8），并假定一个类似于矩形的光谱响应函数，于是宽波段比辐射率可以重新写为：

${\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}1\&RightArrow;\mathrm{\&lambda;}2}=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&lambda;}\left(i\right)}{\overset{\mathrm{\&lambda;}(i+1)}{\&Integral;}}{\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}}{\underset{\mathrm{\&lambda;}1}{\overset{\mathrm{\&lambda;}2}{\&Integral;}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{g}_i{\mathrm{\&epsiv;}}_i^{\&prime;}\&ap;\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{g}_i{\mathrm{\&epsiv;}}_i,---\left(19\right)$

其中，

${\mathrm{\&epsiv;}}_i^{\&prime;}=\frac{\underset{\mathrm{\&lambda;}\left(i\right)}{\overset{\mathrm{\&lambda;}(i+1)}{\&Integral;}}{\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}}{\underset{\mathrm{\&lambda;}\left(i\right)}{\overset{\mathrm{\&lambda;}(i+1)}{\&Integral;}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}},---\left(20\right)$

$g_i=\frac{\underset{\mathrm{\&lambda;}\left(i\right)}{\overset{\mathrm{\&lambda;}(i+1)}{\text{\&Integral;}}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}}{\underset{\mathrm{\&lambda;}1}{\overset{\mathrm{\&lambda;}2}{\&Integral;}}{B}_{\mathrm{\&lambda;}}\left(T\right)\mathrm{d\&lambda;}}.---\left(21\right)$

式中，ε_i是热红外第i波段的比辐射率，n表示波段的个数。

从上式不难发现，宽波段比辐射率与窄波段波段比辐射率可以近似于一种线性关系，线性关系中的系数几乎与温度T无关，即：

ε_3-14＝w₀+w₁*ε₂₉+w₂*ε₃₁+w₃*ε₃₂         （22）

由此可见，要想获取3-14μm的宽波段比辐射率，必须确定拟合系数a₀-a₃。拟合系数主要通过模拟数据来确定，如图4所示。首先，从ASTER地物光谱数据库中提取出岩石、土壤、植被和水体在3-14μm的地表比辐射率谱。然后，结合常温的普朗克函数以及MODIS的波段响应函数，利用方程式（16）和方程式（18）沿着波长方向卷积，分别获取窄波段和宽波段的地表比辐射率。最后，利用最小二乘法获取所需的宽波段比辐射率拟合系数。地表3-14μm宽波段比辐射率拟合系数如下表所示：

表3.宽波段比辐射率拟合系数

（2）地表3-∞μm宽波段比辐射率

由于受测量仪器的限制，现有的地物波谱库中一般只测量了14μm以下的地物波谱，因此简单采用地表3-14μm宽波段比辐射率的确定方法无法获取地表3-∞μm宽波段比辐射率。为了获得14μm以上的波谱信息，首先分析计算了整个ASTER波谱中不同地物类型波谱曲线的平均值和标准差。从中可以发现以下的规律：

1）对于一个给定的波长，比辐射率值变化越大，平均比辐射率值就相对越小；

2）波长越长，波谱比辐射率值的变化相对越小。

这样，对于波长大于14μm的波谱区间，可以用12.5-13.5μm波谱区间的比辐射率平均值来代替计算。除此之外，当温度从260K到340K变化时，对于一个黑体，3-100μm波谱区间发射辐射能在整个电磁波谱区间上所占辐射能的百分比为99.52%。由此可见，3-100μm的宽波段比辐射率可以近似代替3-∞μm的宽波段比辐射率。与确定3-14μm的宽波段比辐射率估算方法类似，3-∞μm的宽波段比辐射率可以表示为：

ε_3-∞＝w′₀+w′₁*ε₂₉+w′₂*ε₃₁+w′₃*ε₃₂            （23）

其中，3-∞μm宽波段比辐射率拟合系数如下表所示：

表3.宽波段比辐射率拟合系数

本发明未详细阐述部分属于本领域公知技术。

以上所述，仅为本发明部分具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本领域的人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (5)

1.一种同时确定地表窄波段和宽波段比辐射率的方法，其特征在于实现步骤如下：

（A）．对中分辨率成像光谱仪（MODIS）数据进行数据处理，包括辐射定标、大气纠正、以及云掩膜处理，获取晴空条件下MODIS数据位于0.4-2.1μm大气窗口区的可见光/近红外波段地表反射率数据；

（B）．利用步骤A的反射率数据，结合窄波段地表比辐射率反演算法，确定MODIS热红外光谱谱段第29,即波谱范围8.4-8.7μm、第31，即波谱范围10.78-11.28μm和第32，即波谱范围11.77-12.27μm波段的窄波段地表比辐射率；

（C）．利用步骤B的窄波段地表比辐射率，结合窄-宽波段比辐射率转换模型，得到地表3-14μm以及3-∞μm的宽波段比辐射率。

2.根据权利要求1所述的一种同时确定地表窄波段和宽波段比辐射率的方法，其特征于所述步骤（A）中对中分辨率成像光谱仪（MODIS）数据进行数据处理的过程是：

（A1）根据MODIS提供的定标系数，对MODIS遥感数据进行绝对辐射定标，将影像灰度值转换为辐射亮度；

（A2）利用快速大气纠正模型（QUAC）和MODIS的云掩膜产品MOD35数据，完成MODIS遥感数据的大气纠正和有云像元的掩膜处理，获取MODIS数据位于0.4-2.1μm大气窗口区的可见光/近红外谱段的地表反射率数据。

3.根据权利要求1所述的一种同时确定地表窄波段和宽波段比辐射率的方法，其特征于所述步骤（B）确定MODIS数据热红外光谱谱段29、31和32波段的窄波段地表比辐射率过程如下：根据计算出的可见光/近红外地表反射率数据，计算归一化的积雪指数（NDSI）和归一化的植被指数（NDVI），并结合MODIS的陆地/海洋掩膜数据，对像元进行分类，分别分为水体像元，雪/冰像元和植被/裸土像元三大类分别得到所述比辐射率；

（B1）对于水体像元，利用ASTER波谱库和UCSB波谱库中测量的一些典型的水体波谱曲线，结合MODIS的第29、31和32波段的光谱响应函数分别进行卷积处理，获得每个波段上比辐射率的平均值为0.984、0.992和0.987，最后将该三个值分别作为MODIS热红外光谱谱段29、31和32波段的比辐射率；

（B2）对于雪/冰像元，利用ASTER波谱库和UCSB波谱库中测量的一些典型的雪/冰波谱曲线，结合MODIS的第29、31和32波段的光谱响应函数分别进行卷积处理，获得每个波段上比辐射率的平均值，最后将该三个值分别作为雪/冰情况下MODIS热红外光谱谱段29、31和32波段的比辐射率；

（B3）对于植被/裸土像元，根据计算出的NDVI值大小，进一步细分为裸土像元、植被像元和植被/裸土混合像元；

针对裸土像元，将MODIS热红外波段比辐射率与可见光/近红外波段多个波段的反射率建立多元线性回归模型，如公式（1）：

ε_si＝a_i+∑b_ijρ_j，其中i=29,31,32;j=1-7            （1）

公式（1）中，ε_si表示为热红外第i波段裸土的比辐射率，p_j表示为可见光/近红外第j波段的反射率；a_i和b_ij是模型转换系数，通过多元回归获取；

针对植被橡元，建立热红外比辐射率与NDVI的关系模型获取热红外光谱谱段的比辐射率值：

ε_vi＝c_0i+c_1iNDVI  ，其中i＝29,31,32              （2）

公式（2）中，ε_vi表示为热红外第i波段植被的比辐射率，c_0i和c_1i是模型转换系数；

针对植被/裸土混合像元，通过建立热红外波段比辐射率与植被覆盖度P_v的关系模型来获取热红外光谱谱段的比辐射率：

ε_i＝d_0i+d_1iP_v                   （3）

公式（3）中，ε_i表示为热红外第i波段在植被/裸土混合区的比辐射率，d_0i和d_1i是模型转换系数，P_v表示为植被覆盖度，其值通过NDVI值来确定：

$P_v={\left[\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right]}^2---\left(4\right)$

式中，NDVI_min和NDVI_max分别为裸土和植被的NDVI值；NDVI_min和NDVI_max值能够从NDVI直方图获取。

4.根据权利要求1所述的一种同时确定地表窄波段和宽波段比辐射率的方法，其特征于所述步骤（C）中确定地表3-14μm以及3-∞μm的宽波段比辐射率过程是：

（C1）根据计算出的窄波段地表比辐射率，结合建立的窄-宽波段比辐射率转换模型（公式（5）），确定3-14μm地表宽波段比辐射率ε_3→14：

ε_3-14＝w₀+w₁*ε₂₉+w₂*ε₃₁+w₃*ε₃₂              (5)

式中，ε₂₉、ε₃₁和ε₃₂分别是MODIS热红外光谱谱段第29、31和32波段的窄波段比辐射率，w_i是模型转换系数；

（C2）根据计算出的窄波段地表比辐射率，结合建立的窄-宽波段比辐射率转换模型（公式（6）），确定3-∞μm地表宽波段比辐射率ε_3→∞：

ε_3-∞＝w′₀+w′₁*ε₂₉+w′₂*ε₃₁+w′₃*ε₃₂           （6）

式中，ε₂₉、ε₃₁和ε₃₂分别是MODIS热红外光谱谱段第29、31和32波段的窄波段比辐射率，w′_i是模型转换系数。

5.一种实现权利要求1所述方法的装置，其特征在于包括：数据预处理模块、窄波段比辐射率反演模块和宽波段比辐射率反演模块，其中：

数据预处理模块：对中分辨率成像光谱仪（MODIS）数据进行数据处理，包括辐射定标、大气纠正、以及云掩膜处理，获取晴空条件下MODIS数据位于0.4-2.1μm大气窗口区的可见光/近红外地表反射率数据；

窄波段比辐射率反演模块：利用数据预处理模块的反射率数据，结合窄波段地表比辐射率反演算法，确定MODIS热红外光谱谱段第29,即波谱范围8.4-8.7μm、第31，即波谱范围10.78-11.28μm和第32，即波谱范围11.77-12.27μm波段的窄波段地表比辐射率；具体为根据计算出的可见光/近红外地表反射率数据，计算归一化的积雪指数（NDSI）和归一化的植被指数（NDVI），并结合MODIS的陆地/海洋掩膜数据，对像元进行分类，分别分为水体像元，雪/冰像元和植被/裸土像元三大类；再结合获取的可见光/近红外地表反射率数据计算的归一化植被指数（NDVI）值的大小，进一步将植被/裸土像元分为裸土像元、植被像元和裸土/植被混合像元三类；针对每一类像元，分别利用自己建立的窄波段比辐射率反演算法，确定不同类别像元的MODIS热红外谱段29、31和32波段的窄波段比辐射率数据；

宽波段比辐射率反演模块：分别利用自己建立的3-14μm窄-宽波段地表比辐射率转换模型和3-∞μm窄-宽波段地表比辐射率转换模型，计算出每个像元3-14μm的地表宽波段比辐射率数据和3-∞μm的地表宽波段比辐射率数据。