CN102103203A - 基于环境一号卫星的地表温度单窗反演方法 - Google Patents

Abstract

本发明公开了一种基于环境一号卫星的地表温度单窗反演方法，包括以下步骤：S1：获取待测区域地表的环境一号B星HJ-1B卫星的遥感数据，对该遥感数据进行预处理，得到影像数据，并将所述影像数据的像元亮度DN值转换为辐亮度；S2：采用基于影像的COST大气校正模型，对所述影像数据进行大气校正；S3：采用归一化植被指数NDVI阈值法估算待测区域地表的比辐射率ε；S4：根据步骤S1～S3的计算结果利用单窗算法反演地表温度。本发明把基于影像的COST模型大气校正方法和不需要大气水汽含量参数的地表温度单窗算法应用于HJ-1B卫星遥感数据，实现了地表温度的反演。

Description

基于环境一号卫星的地表温度单窗反演方法

技术领域

本发明涉及热红外遥感技术领域，特别是涉及一种基于环境一号卫星的地表温度单窗反演方法。

背景技术

地表温度是地表与大气之间辐射及能量交换和区域与全球尺度地表物理过程的关键参数，热红外遥感是探测地表温度的有效手段。环境一号卫星A、B(HJ-A/B)星是我国继气象、海洋、国土资源卫星之后2008年发射的一个全新的民用卫星。环境一号B(HJ-1B)卫星上的红外相机(IRS)的热红外波段空间分辨率为300m。对该星的热红外研究还不多，且主要集中在前期的模型模拟、在轨测评或反演过程相对复杂阶段，如段四波等(2008)利用HJ-1B数据模拟反演了地表温度。在卫星的最新定标参数发布后，目前急需结合地面观测数据或者其它卫星数据交叉标定，定量评价该星在环境和减灾领域中的应用能力。

地表温度反演的方法目前主要有单窗算法(也称单通道算法)、劈窗算法(也称双通道或分裂窗算法)和多通道算法。其中应用最多的是单窗算法和劈窗算法。单窗算法是利用卫星传感器上单独的一个热红外通道(一般在大气窗口内)获得的辐射能，借助于无线电探空或卫星遥感确定的大气廓线数据(温度、湿度、压力)，结合辐射传输方程来修正大气和比辐射率的影响，从而得出地表温度。

国内外利用单窗算法反演地表温度的研究已经非常成熟，其中较 具代表性的主要有辐射传输方程法、Jiménez-&Sobrino单窗算法(2003)和覃志豪单窗算法(Qin et al，2001)。辐射传输方程法较复杂，所需参数较多，在反演地表温度的实际中不便应用，而后两种虽然仅需要地表温度、大气水汽含量和地表比辐射率(覃志豪单窗算法还需要大气平均作用温度)，但大气水汽含量并不容易获取，尽管通过其它遥感数据，如MODIS(是搭载在terra和aqua卫星上的一个重要的传感器)的水汽产品或者波段比值法获得大气水汽分布，但不同传感器、不同分辨率的数据融合时，也易带来误差。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：如何在不获取大气水汽含量参数的情况下，仅利用遥感影像数据实现地表温度的反演。

(二)技术方案

为了解决上述技术问题，本发明提供了一种基于环境一号卫星的地表温度单窗反演方法，包括以下步骤：

S1：获取待测区域地表的环境一号B星HJ-1B卫星的遥感数据，对该遥感数据进行预处理，获得影像数据，并将所述影像数据的像元亮度值DN值转换为辐亮度；

S2：采用基于影像的COST大气校正模型，对所述影像数据进行大气校正；

S3：采用归一化植被指数NDVI阈值法估算待测区域地表的比辐射率ε；

S4：根据步骤S1～S3的计算结果利用单窗算法反演地表温度。

其中，所述影像数据包括电荷耦合器件CCD图像和热红外波段图像，所述步骤S1中将所述影像数据的像元亮度值转换为辐亮度的过程具体包括：

利用绝对定标系数将CCD图像的DN值转换为辐亮度的公式为：

$L=\frac{\mathrm{DN}}{A}+L_0$

式中A为绝对定标系数增益，L₀为绝对定标系数偏移量，转换后辐亮度的单位为W·m^-2·sr^-1·μm^-1，W表示瓦，sr表示单位立体角；

对于红外相机的第四波段的热红外波段图像，将其DN值转换为辐亮度的公式为：

$L=\frac{\mathrm{DN}-b}{g}$

其中g为绝对定标系数增益，b为绝对定标系数偏移量。

其中，所述步骤S2中，

COST大气校正模型的表达式如下：

ρ＝πd²(L_sat-L_p)/(E₀cosθ_zT_z)

上式中，ρ为地面相对反射率；d为日地天文单位距离；L_sat为入瞳辐射亮度，即大气顶层的辐射能量；L_p为大气层辐射值；E₀为大气顶层的太阳平均光谱辐射，即大气顶层太阳辐照度；θ_z为太阳天顶角，T_z是大气透过率。

其中，所述步骤S3具体为：

采用以下公式估算地表的比辐射率ε：

$\left\{\begin{array}{c}\mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_s,\mathrm{NDVI}<0.2\\ \mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_v{P}_v+{\mathrm{\&epsiv;}}_s(1-P_v)+\mathrm{d\&epsiv;},0.2\&le;\mathrm{NDVI}\&le;0.5\\ \mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_v,\mathrm{NDVI}>0.5\end{array}\right.$

其中NDVI、P_v和dε分别表示归一化植被指数、地表的植被比率和比辐射率误差，ε_s、ε_v分别是地表的裸土和植被的比辐射率。

其中，所述步骤S4具体为：

利用如下公式反演地表温度：

$T_s=\frac{T_B}{1+(\mathrm{\&lambda;}\&times;T_B/C)*\ln \mathrm{\&epsiv;}}$

式中，T_s是地表温度，T_B是亮度温度，λ是有效中心波长，C是常数，ε是地表的比辐射率。

其中，所述预处理包括几何校正、镶嵌、拼接、去云和辐射校正处理。

(三)有益效果

本发明把基于影像的COST模型大气校正方法和不需要大气水汽含量参数的地表温度单窗算法应用于HJ-1B卫星遥感数据，实现了地表温度的反演。

附图说明

图1为本发明实施例的方法流程图；

图2示出了HJ-1B卫星的热红外通道响应函数；

图3为基于本发明的方法反演的地表温度和MODIS温度产品的比较(300m)图；

图4为基于本发明的方法反演的地表温度和MODIS温度产品的比较(1000m)图。

具体实施方式

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

环境一号卫星只有一个合适的热红外波段(10.5-12.5um)，因此，本发明采用一种不需大气水汽含量参数的简便的单窗算法(是现有技术，由Artis和Carnahan在1982年提出；Weng等在2004年进行了应用)来反演地表温度。

图1示出了本发明实施例的方法流程图，如图1所示，该方法包括以下步骤：

S1、预处理

对宁夏研究区的遥感数据(即HJ-1B卫星的光学和热红外影像数据，所述光学图像数据为CCD(电荷耦合器件)图像，热红外影像数据为热红外波段图像)(数据见表1)进行预处理，包括几何校正、镶嵌、拼接、裁剪、去云和辐射定标。其中，HJ-1B卫星的热红外通道响应函数如图2。

表1环境卫星数据

表1中，CCD2表示HJ-1B上携带的两个CCD相机。按公式(1)利用绝对定标系数将CCD图像的DN值(是指遥感影像数据的像元亮度值)转换为辐亮度L：

$L=\frac{\mathrm{DN}}{A}+L_0---\left(1\right)$

式中A为绝对定标系数增益(单位是DN/(W·m^-2·sr^-1·μm^-1))，L₀为绝对定标系数偏移量，转换后的辐亮度单位为W·m^-2·sr^-1·μm^-1，其中，W表示瓦，sr表示单位立体角。

对于IRS-Band4(红外相机的第四波段)热红外波段图像，其辐亮度L的公式为：

$L=\frac{\mathrm{DN}-b}{g}---\left(2\right)$

其中g为绝对定标系数增益，b为绝对定标系数偏移量。这些绝对定标系数和大气层外太阳辐照度均来自环境卫星官方网站上公布的资料，详见表2。

表2HJ-1B波段的绝对定标系数和大气层外太阳辐照度

S2、大气校正

大气校正主要采用理论模型和基于影像的模型两种方法(可参考文献吾拉木等2004；彭光雄等，2007；Lu et al，2002)。理论模型如MODTRAN(Moderate Resolution Transmission)(可参考文献Berk et al，1983)或6S(Second Simulation of Satellite Signal in the Solar Spectrum)模型(可参考文献Vermote et al，1997)，二者都建立在辐射传输理论基础之上，但模型复杂，所需同步和实测参数较多。基于影像的如DOS模型、COST模型(可参考文献Moran et al，1992；Chavez，1996)等，这类模型的特点是所需参数少，不需要同步的实测资料。所以，本发明采用基于影像的COST模型。

该模型基于同一影像中大气影响相同的假设，因此每一波段均存在反射率为1％的黑体辐射，黑体辐射值取决于大气顶层的太阳辐照度；另外，大气性质是均一的，传感器每一波段的最小辐射亮度值除黑体辐射影响外，主要还有大气散射和反射作用的影响，需考虑大气透过率的影响。模型的表达式如下(可参考文献Chavez，1996；Sobrino et al，2004)：

ρ＝πd²(L_sat-L_p)/(E₀cosθ_zT_z)          (3)

式中ρ为地面相对反射率；d为日地天文单位距离；L_sat为入瞳辐射亮度，即大气顶层的辐射能量；L_p为大气层辐射值；E₀为大气顶层的太阳平均光谱辐射，即大气顶层太阳辐照度；θ_z为太阳天顶角，T_z是大气透过率。

这里，日地天文单位距离可以查相关手册获取或者通过儒略日的有关公式计算，本发明中其值取为1左右。L_sat即为公式(1)计算的辐亮度，大气层外太阳辐照度E₀见表2。太阳天顶角θ_z可以通过影像头文件的太阳高度角获取(二者互余)。因为HJ-1B星CCD2相机第三、四波段和TM波段范围相似，所以T_z的值参考TM波段设置，这里也分别取0.85和0.91(可参考文献Chavez，1996)。然而，在获取不到T_z值的时候，Chavez(1996)指出可以近似用太阳天顶角的余弦代替，不过易带来误差。

L_p的计算公式如下：

L_p＝L_min-L_1％         (4)

其中L_min为遥感器每一波段最小光谱辐射值，可由上文对应的辐亮度计算公式获取。L_1％为反射率为1％的黑体辐射值，其计算公式为：

$L_{1\%}=\frac{0.01\cos {\mathrm{\&theta;}}_z{T}_z{E}_0}{\mathrm{\&pi;}{d}^2}---\left(5\right)$

S3、比辐射率的估算

采用NDVI阈值法(可参考文献Sobrino et al，2001；Sobrino et al，2004)，估算地表比辐射率ε：

$\left\{\begin{array}{c}\mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_s,\mathrm{NDVI}<0.2\\ \mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_v{P}_v+{\mathrm{\&epsiv;}}_s(1-P_v)+\mathrm{d\&epsiv;},0.2\&le;\mathrm{NDVI}\&le;0.5\\ \mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_v,\mathrm{NDVI}>0.5\end{array}\right.---\left(11\right)$

其中NDVI(Normalized Difference Vegetation Index，归一化植被指数)、P_v和dε分别表示归一化植被指数、植被比率和比辐射率误差，ε_s、ε_v分别是裸土和植被的比辐射率，取值分别为0.973和0.99。

P_v和dε表示如下：

NDVI＝(ρ₄-ρ₃)/(ρ₄+ρ₃)       (12)

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

dε＝(1-ε_s)(1-P_v)Fε_v          (14)

ρ₃，ρ₄分别表示HJ-1B-CCD2的第3(红外)和第4(近红外)波段的反射率。NDVI_max，NDVI_min分别取0.2和0.5。Fε_v为常数0.55。

S4、无需大气水汽含量的地表温度单窗反演算法：

$T_s=\frac{T_B}{1+(\mathrm{\&lambda;}\&times;T_B/C)*\ln \mathrm{\&epsiv;}}---\left(6\right)$

式中，T_s是地表温度，T_B是亮度温度，λ是有效中心波长，C是常数(C＝hc/σ，h为普朗克常数，c为光速，σ为波尔兹曼常数)，C的值为1.43876869×10^-2m K，ε是地表比辐射率。

利用式(7)的通道响应函数(如)获取有效中心波长λ(Jim énez-&Sobrino，2003)：

$\mathrm{\&lambda;}=\frac{\&Integral;\mathrm{\&lambda;f}\left(\mathrm{\&lambda;}\right)\mathrm{d\&lambda;}}{\&Integral;f\left(\mathrm{\&lambda;}\right)\mathrm{d\&lambda;}}---\left(7\right)$

经过上式的积分比值，得到λ＝11.55μm。然而，在没有通道响应曲线时，可用中心波长代替，但也易带来误差。

亮度温度的计算可以从普朗克方程(式(8))推导(Jim énez-

&Sobrino，2003)：

$B(\mathrm{\&lambda;},T)=\frac{c_1{\mathrm{\&lambda;}}^{-5}}{\exp \left(\frac{c_2}{\mathrm{\&lambda;}{T}_s}\right)-1}---\left(8\right)$

式中，B(λ，T)是辐射亮度，B是λ、T(表示温度)的函数，λ是有效中心波长，T_s是地表温度，c₁和c₂是光谱常数(c₁＝1.19104×10⁸Wμm⁴m^-2sr^-1，c₂＝1.43877×10⁴μm K)。从上式反推出亮度温度T_B的计算公式：

$T_B=\frac{c_2}{\mathrm{\&lambda;}\ln [1+\frac{c_1}{{\mathrm{\&lambda;}}^{5}B(\mathrm{\&lambda;},T)}]}---\left(9\right)$

为方便计算，令K₁＝c₁/λ⁵和K₂＝c₂/λ，对上式简化得：

T_B＝K₂/ln(1+K₁/B(λ，T))           (10)

式中，K₁、K₂是常量，K₁＝579.20W m^-2sr^-1μm^-1，K₂＝1245.58K。

接下来可以对反演结果和精度验证。

对于地表温度反演的精度验证，由于同步的实测数据比较难于获取，况且实测值一般是点上的温度信息，而遥感数据代表面上的信息，所以本发明采用Terra/MODIS温度产品验证。MODIS温度产品的空间分辨率为1km，虽然较环境卫星的低，但MODIS温度产品比较成熟，具有小于1K的精度，并且Terra星和环境星的过境时间都是地方时上午10:30，因此具有同步的优势，也避免了由点到面转换带来的误差。

对一个反演算法的对比验证可根据实测值或其它标准值，由于实测数据较难获取，所以我们采用具有较高精度(＜1K)的MODIS温度产品作为标准值，来评价该反演算法。HJ-1B数据的热红外通道的空间分辨率为300m，而MODIS温度产品的空间分辨率为1km，为了全面有效地对该星热红外波段的地表温度探测能力进行评判，本文采用两个思路，一种是把MODIS的空间分辨率重采样为300m，另一种是把HJ-1B的数据重采样为1km，进而对同一种分辨率下的反演结果对比。反演结果的评价方法，本发明采用最小值、最大值和平均值这些统计指标。

利用上述的方法，对HJ-1B影像进行辐射校正、简单去云、大气校正等预处理，进而采用单窗算法反演地表温度，其反演的地表温度和MODIS温度产品的影像对比见图4和5，统计值比较结果如表3。从表3可以看出，对HJ-1B本身的分辨率来说(300m)，其反演的最大值误差为0.12K，平均值误差为-0.96K。重采样为1km分辨率后，其反演的最大值误差为-0.25K，平均值误差为-0.88K。与HJ-1B其星上的定标精度2K相比，两种思路下该单窗算法的平均误差均小于1K，表明该法具有可信的精度。图3和4的LST(Land Surface Temperature，地表温度)分布图也证明了该方法的有效性。

表3HJ-1B反演的地表温度和MODIS温度产品的统计值比较

由以上实施例可以看出，本发明实施例首先对影像进行几何校正、镶嵌、拼接、去云和辐射校正等预处理，然后通过基于影像的COST模型进行大气校正，再用NDVI阈值法估算地表比辐射率，最后把亮温、地表比辐射率等参数输入到无需大气水汽含量的单窗算法，进而反演地表温度。本发明可以应用于其它单通道热红外遥感数据，如TM等。其特点是简单、快速、准确，仅靠遥感影像而无需大气水汽参数即可利用单通道的热红外数据反演地表温度。

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

Claims (6)

1.一种基于环境一号卫星的地表温度单窗反演方法，其特征在于，包括以下步骤：

S1：获取待测区域地表的环境一号B星HJ-1B卫星的遥感数据，对该遥感数据进行预处理，得到影像数据，并将所述影像数据的像元亮度DN值转换为辐亮度；

S2：采用基于影像的COST大气校正模型，对所述影像数据进行大气校正；

S3：采用归一化植被指数NDVI阈值法估算待测区域地表的比辐射率ε；

S4：根据步骤S1～S3的计算结果利用单窗算法反演地表温度。

2.如权利要求1所述的基于环境一号卫星的地表温度单窗反演方法，其特征在于，所述影像数据包括电荷耦合器件CCD图像和热红外波段图像，所述步骤S1中将所述影像数据的像元亮度值转换为辐亮度的过程具体包括：

利用绝对定标系数将CCD图像的DN值转换为辐亮度的公式为：

$L=\frac{\mathrm{DN}}{A}+L_0$

式中A为绝对定标系数增益，L₀为绝对定标系数偏移量，转换后辐亮度的单位为W·m^-2·sr^-1·μm^-1，W表示瓦，sr表示单位立体角；

将热红外波段图像的DN值转换为辐亮度的公式为：

$L=\frac{\mathrm{DN}-b}{g}$

其中g为绝对定标系数增益，b为绝对定标系数偏移量。

3.如权利要求1所述的基于环境一号卫星的地表温度单窗反演方法，其特征在于，所述步骤S2中所使用的COST大气校正模型的表达式如下：

ρ＝πd²(L_sat一L_p)/(E₀cosθ_zT_z)

上式中，ρ为地面相对反射率；d为日地天文单位距离；L_sat为入瞳辐射亮度，即大气顶层的辐射能量；L_p为大气层辐射值；E₀为大气顶层的太阳平均光谱辐射，即大气顶层太阳辐照度；θ_z为太阳天顶角，T_z是大气透过率。

4.如权利要求1所述的基于环境一号卫星的地表温度单窗反演方法，其特征在于，所述步骤S3中，采用以下公式估算地表的比辐射率ε：

$\left\{\begin{array}{c}\mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_s,\mathrm{NDVI}<0.2\\ \mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_v{P}_v+{\mathrm{\&epsiv;}}_s(1-P_v)+\mathrm{d\&epsiv;},0.2\&le;\mathrm{NDVI}\&le;0.5\\ \mathrm{\&epsiv;}={\mathrm{\&epsiv;}}_v,\mathrm{NDVI}>0.5\end{array}\right.$

其中NDVI、P_v和dε分别表示归一化植被指数、地表的植被比率和比辐射率误差，ε_s、ε_v分别是地表的裸土和植被的比辐射率。

5.如权利要求4所述的基于环境一号卫星的地表温度单窗反演方法，其特征在于，所述步骤S4中，利用如下公式反演地表温度：

$T_s=\frac{T_B}{1+(\mathrm{\&lambda;}\&times;T_B/C)*\ln \mathrm{\&epsiv;}}$

式中，T_s是地表温度，T_B是亮度温度，λ是有效中心波长，C是常数，ε是地表的比辐射率。

6.如权利要求1～5任一项所述的基于环境一号卫星的地表温度单窗反演方法，其特征在于，所述预处理包括几何校正、镶嵌、拼接、去云和辐射校正处理。

Images