CN102338871A - 地表反射率计算方法和装置 - Google Patents

Abstract

本发明涉及一种用于HJ-1A/B卫星的地表反射率计算方法和装置，其包括步骤：1)利用HJ-1B CCD及红外相机的中红外波段，基于暗目标的方法得到大气气溶胶光学厚度；利用HJ-1A大幅宽快速重访的特点，基于不变目标方法得到大气气溶胶光学厚度；2)基于HJ-1A/BXML文件和影像数据计算每日周期的太阳天顶角、太阳方位角、观测天顶角和观测方位角；3)基于上述步骤1)和2)中得到的参数，利用MODTRAN模拟星上幅亮度L_m；4)通过步骤3)建立地表反射率查找表。5)根据大气参数和待校正图像，利用查找表实现大气校正。利用本发明能有效的消除大气对遥感影像的吸收和散射，恢复地表目标的反射率，解决了行业应用的瓶颈，进一步扩展环境减灾小卫星数据的应用范围。

Description

地表反射率计算方法和装置

技术领域

本发明涉及一种地表反射率计算方法和装置，具体而言，涉及一种用于环境减灾小卫星HJ-1的地表反射率计算方法和装置。

背景技术

“环境与灾害监测预报小卫星星座”(以下简称“HJ-1”)是中国第一个专门用于环境与灾害监测预报的小卫星星座，其中HJ-1A/B星于2008年9月6日发射。卫星载荷如下表：

HJ-1A/B在以后的灾害监测和预报中起着越来越重要的作用。对图像进行大气校正处理，生成地表反射率产品，这对于后期产品(如植被指数、叶面积指数、土壤湿度/含水量、植被覆盖度)的生成尤为重要。如果大气校正的方便程度和校正精度能达到各行业对环境减灾小卫星数据的要求，将进一步推动环境减灾小卫星数据的使用，推动各个行业的发展。

由于大气的分子、气溶胶的散射以及臭氧、水汽等气体的吸收均可以削弱传感器接收到的信号，所以航空、航天遥感平台上的传感器接收到的地物信息不能真实地反映地表。因此准确获得地表信息非常有必要，它也是地表参数遥感定量反演的一个必备环节。而目前还没有针对HJ-1A/B载荷及成像特点的地表反射率反演的计算方法和装置。本方法和装置实现了HJ-1A/B的CCD和HSI数据的地表反射率的反演。

发明内容

对此，本发明提供了一种针对HJ-1A/B的地表反射率计算方法和装置。

本方法和装置首先理论上分析大气传输过程中的大气影响因子，然后利用MODTRAN对大气-地表-遥感器之间的辐射传输过程进行模拟，依据卫星同步反演得到的大气参数和辐亮度，针对HJ-1A/B传感器的特性，生成适用的查找表，并最终反演出地表反射率。

本发明的方法包括步骤：1)利用HJ-1B CCD及红外相机的中红外波段，基于暗目标的方法得到大气气溶胶光学厚度；利用HJ-1A大幅宽快速重访的特点，基于不变目标方法得到大气气溶胶光学厚度；2)基于HJ-1A/B XML文件和影像数据计算每日周期的太阳天顶角、太阳方位角、观测天顶角和观测方位角；3)基于上述步骤1)和2)中得到的参数，利用MODTRAN模拟星上辐亮度L_m；4)通过步骤3)建立地表反射率查找表。5)根据大气参数和待校正图像，利用查找表实现大气校正。

本发明的地表反射率计算装置包括：大气参数反演模块、太阳及观测几何参数反演模块和地表反射率计算模块；

大气参数反演模块利用暗目标法和不变目标法，基于经典的大气辐射传输理论，进行气溶胶光学厚度反演，得出地表反射率ρ与空间观测星载反射率(也称表观反射率)ρ’之间的关系。

太阳及观测几何参数反演模块基于HJ-1-A/B CCD 30米分辨率DN值数据和HJ-1-A/B CCD XML文件，计算得到具有地理参考的太阳天顶角数据、太阳方位角数据、观测天顶角数据和观测方位角数据。

地表反射率计算模块基于HJ-1-A/B CCD、HSI全波段Radiance数据或DN值和水平气象视距，太阳天顶角等大气参数，利用大气辐射传输模型MODTRAN来建立通用查找表，实现逐像元大气校正，得到CCD，HSI波段的地表反射率数据。

本发明的方法和装置具有较高的准确性，而且由于是在MODIS数据大气校正算法的基础上进行了改进，具有一定的创新性。

本发明的方法和装置对环境减灾小卫星数据进行大气校正。能有效的消除大气对遥感影像的吸收和散射，恢复地表目标的反射率，解决了行业应用的瓶颈，进一步扩展环境减灾小卫星数据的应用范围。

附图说明

以下结合附图详细说明本发明，其中：

图1为本发明的地表反射率计算方法的流程图；

图2为太阳天顶角(SZA)和MODTRAN运行的光谱幅亮度(mod_sr)的关系图；

图3为水平气象视距(VIS)和MODTRAN运行的光谱幅亮度(mod_sr)的关系图；

图4为发明的地表反射率计算方法的示意图；以及

图5为本发明的地表反射率计算装置的结构示意图。

具体实施方式

下面参照图1和图4详细描述本发明：

(一)大气参数反演

对大气参数的反演利用了暗目标法和不变目标法。

暗目标法：其基本原理就是在假定待校正的遥感图像存在黑暗像元区域、地表朗伯面反射、大气性质均一，忽略大气多次散射辐照作用和邻近像元漫反射作用的前提下，反射率很小的黑暗像元是由于大气的影响，而使得这些像元的反射率相对增加，可以认为这部分增加的反射率是由于大气程辐射的影响产生的。然而，如果区域内植被稀疏，该方法就不适用。由于要求存在水体或浓密植被的条件，使得该方法在应用于冬季获取的地表影像上受到很大的限制；对于环境星地表反射率反演的业务化、批量处理的要求而言。暗目标方法需要人工确定暗目标点等人及交互操作；而且对于不同植被条件下的中红外波段和可见光/近红外之间的经验关系是不同的，因此在应用该方法之前，必须有大量的实地影像的统计基础。

不变目标法：假定图像上存在具有较稳定反射辐射特性的像元，并且可确定这些像元的物理意义，那么就称这些像元为不变目标，这些不变目标在不同时相的遥感图像上的反射率将存在一种线性关系。当确定了不变目标以及它们在不同时相遥感图像中反射率的这种线性关系，就可以对遥感图像进行大气校正。该方法简单、直接，但是本质上是一种统计方法，对不同时间和区域的影像大气校正的应用程度值得考虑，也是一种相对校正方法，而很难校正异质大气气溶胶散射的影响是该方法的另一个先天性缺陷。

本发明基于HJ-1B 4谱段CCD相机和MODIS/Terra多光谱遥感数据，利用HJ-1B和MODIS/Terra双星协同反演陆地下垫面1000米高空间分辨率3个波段的大气气溶胶光学厚度遥感定量反演。其中，数据来源为：HJ-1B CCD 1-3波段DN数据；HJ-1B红外相机第2波段DN数据；MODIS/Terra第1、3、4、6波段DN数据；HJ-1B成像时间及传感器ID；以及MODIS成像太阳及观测几何数据。

本发明基于经典的大气辐射传输理论，推导得出地表反射率ρ与空间观测星载反射率(也称表观反射率)ρ′间的关系表达式：

${\mathrm{\ρ}}_{A,{\mathrm{\λ}}_i}=\frac{\left({{{\mathrm{\ρ}}_{A,{\mathrm{\λ}}_i}}^{\′}b-a}_A\right)+a_A(1-{{\mathrm{\ρ}}_{A,{\mathrm{\λ}}_i}}^{\′})e^{(a_A-b)\mathrm{\ϵ}(0.00879{\mathrm{\λ}}_i^{-4.09}+{\mathrm{\β}}_A{\mathrm{\λ}}_i^{-\mathrm{\α}})}}{({{\mathrm{\ρ}}_{A,{\mathrm{\λ}}_i}}^{\′}b-a_A)+b(1-{{\mathrm{\ρ}}_{A,{\mathrm{\λ}}_i}}^{\′})e^{(a_A-b)\mathrm{\ϵ}(0.00879{\mathrm{\λ}}_i^{-4.09}+{\mathrm{\β}}_A{\mathrm{\λ}}_i^{-\mathrm{\α}})}}---\left(1\right)$

${\mathrm{\ρ}}_{B,{\mathrm{\λ}}_i}=\frac{\left({{{\mathrm{\ρ}}_{B,{\mathrm{\λ}}_i}}^{\′}b-a}_B\right)+a_B(1-{{\mathrm{\ρ}}_{B,{\mathrm{\λ}}_i}}^{\′})e^{(a_B-b)\mathrm{\ϵ}(0.00879{\mathrm{\λ}}_i^{-4.09}+{\mathrm{\β}}_B{\mathrm{\λ}}_i^{-\mathrm{\α}})}}{({{\mathrm{\ρ}}_{B,{\mathrm{\λ}}_i}}^{\′}b-a_B)+b(1-{{\mathrm{\ρ}}_{B,{\mathrm{\λ}}_i}}^{\′})e^{(a_B-b)\mathrm{\ϵ}(0.00879{\mathrm{\λ}}_i^{-4.09}+{\mathrm{\β}}_B{\mathrm{\λ}}_i^{-\mathrm{\α}})}}---\left(2\right)$

其中：a＝secθ，b＝2，ε是后向散射系数，通常取值为0.1，α是波长指数，β是消光系数。太阳天顶角可由经度、纬度以及卫星过境时间计算得到。大气光学厚度取决于大气的浊度状态。

分别代表MODIS和HJ-1B影像数据的地表真实反射率，其中i＝1，2，3，分别代表HJ-1B影像数据的0.47，0.55，0.66μm等三个波段与MODIS对应波段。此时，我们将得到未知量远大于已知观测值的方程组(9个未知量，6个方程)，属于不定解问题。

然而，地表各波段反射率的变化与地表物质在分子尺度上对光的反射和吸收有关，一般可看为波长的强函数。气溶胶和大气分子的散射也是波长的函数。地表的二向反射特性主要受地表物质的宏观物理结构影响，随波长的变化较小。因此可以认为，在可见光和近红外波段，地表的二向反射率与波长无关，而与观测的空间几何位置关系有关。实验表明，不同角度的两次观测的地表真实反射率间存在如下关系：

$k_{\mathrm{\λ}}=\frac{{\mathrm{\ρ}}_{1,\mathrm{\λ}}}{{\mathrm{\ρ}}_{2,\mathrm{\λ}}}---\left(3\right)$

k为与波长无关的比例系数，A_1，λi是其中一个角度观测的地表真实反射率，A_2，λi是另外一个角度观测的地表真实反射率。

对于HJ-1B与MODIS双星观测而言，同样存在此关系：

$k_{\mathrm{\λi}}=\frac{{\mathrm{\ρ}}_{A,\mathrm{\λi}}}{{\mathrm{\ρ}}_{B,\mathrm{\λi}}}---\left(4\right)$

因为除了沙尘外，其他气溶胶的消光系数随着观测波长的增加而减小。当观测的波长为短波近红外波段时，气溶胶的光学厚度已经远远小于可见光通道，如果忽略大气分子对短波近红外波段的贡献，k_λ在短波近红外波段处的值可以近似地表达为卫星两次过境时刻的大气顶部表观反射率的比值。又因为k被认为与波长无关，所以短波近红外波段的k值与其他可见光通道(λ_i＝0.47，0.55，0.66μm)地表反射率的比例系数k相等。而搭载在HJ-1B上的红外相机正好有一通道位于1.6μm附近，因而利用CCD相机覆盖的0.47，0.55，0.66μm波段辐射亮度数据和红外相机的1.6μm辐射亮度数据，我们将得到由三个方程组成的闭合的非线形方程组，尽管未知量的个数与方程个数一样，还是很难直接获得非线形方程组的解析解。因此，本算法将采用数值计算方法：准牛顿迭代法，近似解算出大气气溶胶光学厚度数值所需要的波长指数α和

浊度系数β。至此，三个可见光波段上的气溶胶光学厚度值得以反演。

(二)太阳及观测的几何要素计算

对于传统的小幅宽高分辨率卫星影像，大气校正可以不考虑像元间太阳及观测几何的差异。而HJ-1A/B的CCD相机均为大幅宽、倾斜成像，整景影像上像元间太阳天顶角/方位角，观测天顶角/方位角较大，所以在地表反射率反演之前必须进行相应几何要素计算。以下的计算是基于HJ-1A/B XML文件(XML文件是该影像数据的说明信息，包括时间、轨道、投影、地理位置等信息。)

1、太阳天顶角计算

基于HJ-1A/B XML文件和影像数据实现每日周期的太阳天顶角计算。所使用的数据为HJ-1-A/B CCD 30米分辨率DN值数据和HJ-1-A/BCCD XML文件进行影像几何位置计算和太阳天顶角计算，得到具有地理参考的太阳天顶角数据。

计算方法如下：

地球上某点的切面的法线与此点和太阳连线的夹角称为此时此地的太阳天顶角。太阳高度角和太阳天顶角互为余角。依据下面的公式就可以计算出太阳天顶角的大小。

cosθ_s＝sinγ·sinδ+cosγ·cosδ·cosω

则

θ_s＝acos(sinγ·sinδ+cosγ·cosδ·cosω)

为了得到太阳天顶角θ_s，需要先求出赤纬δ和时角ω，他们的计算方法分别如下：

赤纬δ的计算：

(1)计算日地距离

日地距离＝1.000423+0.032359sinα+0.000086sin2α-0.008349cosα+0.000115cos2α

式中α称日角，即α＝2πt/365.2422

这里t又由两部分组成，即t＝N-N₀

式中N为积日，所谓积日，就是日期在年内的顺序号。

N₀＝79.6764+0.2422×(年份-1985)-0.25×(年份-1985)

(2)计算赤纬角

$\mathrm{\δ}=(0.006894-0.399512\cos \mathrm{\α}+0.072075\sin \mathrm{\α}-0.006799\cos 2\mathrm{\α}$

$+0.000896\sin 2\mathrm{\α}-0.002689\cos 3\mathrm{\α}+0.001516\sin 3\mathrm{\α})\×\frac{180}{\mathrm{\π}}$

时角ω的计算：

(1)计算当天的太阳角

自公历年1月1日零时起作为起始点，计算地球绕太阳公转的角度，该角度记为β.

β＝2.0*π*(nDayNumber-1)/365.0

其中nDayNumber指的是该天是一年中的第几天。

(2)计算时差

时差是平太阳时和真太阳时的差值，我们计算需要的是真太阳时。用TimeDifference表示时差，即

TimeDifference＝(0.000043+0.002061cosβ-0.032040sinβ-0.014974cos2β-0.040685sin2β)*229.183

(3)计算真实的时间

真实的时间差需要对前面计算得到的时间数加一个改正量RealTimeNumber＝TimeNumber+TimeDifference/60+(m_longn-120.0)*1/15

其中TimeNumber是指将平太阳时间换算成以小时为单位的时间，m_longn是观测点的经度。

(4)计算时角

时角TimeAngle＝(12-RealTimeNumber)*15*π/180

2、太阳方位角计算

基于HJ-1A/B XML文件和影像数据实现每日周期的太阳方位角计算，所使用的数据为HJ-1-A/B CCD 30米分辨率DN值数据和HJ-1-A/BCCD XML文件进行影像几何位置计算和太阳方位角计算，得到具有地理参考的太阳方位角数据。

计算方法如下：

地球上某点的切面的法线与此点和太阳连线的夹角称为此时此地的太阳天顶角。太阳高度角和太阳天顶角互为余角。依据下面的公式就可以计算出太阳天顶角的大小。

则

为了得到太阳天顶角θ_s，需要先求出太阳高度角h和赤纬δ。

赤纬δ的计算：

(1)计算日地距离

日地距离＝1.000423+0.032359sinα+0.000086sin2α-0.008349cosα+0.000115cos2α

式中α称日角，即α＝2πt/365.2422

这里t又由两部分组成，即t＝N-N₀

式中N为积日，所谓积日，就是日期在年内的顺序号。

N₀＝79.6764+0.2422×(年份-1985)-0.25×(年份-1985)

(2)计算赤纬角

$\mathrm{\δ}=(0.006894-0.399512\cos \mathrm{\α}+0.072075\sin \mathrm{\α}-0.006799\cos 2\mathrm{\α}$

$+0.000896\sin 2\mathrm{\α}-0.002689\cos 3\mathrm{\α}+0.001516\sin 3\mathrm{\α})\×\frac{180}{\mathrm{\π}}$

3、观测天顶角计算

基于HJ-1A/B XML文件和影像数据实现每日周期的卫星观测天顶角计算，所使用的数据为HJ-1-A/B CCD 30米分辨率DN值数据和HJ-1-A/B CCD XML文件进行影像几何位置计算和观测天顶角计算，得到具有地理参考的卫星观测天顶角数据。

计算方法如下：

卫星的观测天顶角是指观测方向与法向方向的夹角，卫星是按照设计的轨道运行的，也就是说在某个时间点上的轨道位置是已知的，因此就可以利用卫星位置星下点的经纬度和地球半径、轨道高度和经纬差计算出卫星的观测天顶角。计算公式如下：

$\tan \mathrm{\θ}=\frac{R\·\sqrt{1-{\cos }^2{\mathrm{\θ}}_1\·{\cos }^2{\mathrm{\θ}}_2}}{H+(R-R\·\cos {\mathrm{\θ}}_1\·\cos {\mathrm{\θ}}_2)}$

则

$\mathrm{\θ}=a\tan \left(\frac{R\·\sqrt{1-{\cos }^2{\mathrm{\θ}}_1\·{\cos }^2{\mathrm{\θ}}_2}}{H+(R-R\·\cos {\mathrm{\θ}}_1\·\cos {\mathrm{\θ}}_2)}\right)$

其中，

R：地球半径

H：卫星高度

θ₁：卫星的星下点与当前点在经度方向上的弧度差

θ₂：是指卫星的星下点与当前点在纬度方向上的弧度差。

4、观测方位角计算

基于HJ-1A/B XML文件和影像数据实现每日周期的卫星观测方位角计算，所使用的数据为HJ-1-A/B CCD 30米分辨率DN值数据和HJ-1-A/B CCD XML文件进行影像几何位置计算和卫星观测方位角计算，得到具有地理参考的卫星观测方位角数据。

计算方法如下：

观测天顶角是在三维的球体中进行计算的，而观测方位角是在二维观测平面中计算的，因此省去了卫星高度和地球半径两个参数，具体公式如下：

则

其中，

θ₁：卫星的星下点与当前点在经度方向上的弧度差

θ₂：是指卫星的星下点与当前点在纬度方向上的弧度差。

(三)地表反射率反演

1、判定太阳天顶角(SZA)的变化间距

在判定太阳天顶角(SZA)的变化间距前，首先看下太阳天顶角(SZA)和MODTRAN运行的光谱幅亮度(mod_sr)的关系图，如图2。

下面开始判定太阳天顶角SZA的插值间距。

SZA最原始的点：0.0，6.0，12.0，24.0，35.2，48.0，54.0，60.0，66.0。SZA的第一次线性插值(取相邻点的中点)：3，9，18，29.6，41.6，51，57，63。

对上面自变量SZA的插值点，因变量光谱幅亮度可以通过线性插值计算和MODTRAN运行结果两种方法得到。根据前面介绍的插值方法和精度判别方法，精度在1％以内的插值点是可以接受的，也就是在其附近不需要再插值；对于大于1％的插值点，需要在其附近再次插值。

因此最终太阳天顶角(SZA)插入的值是：

0.0，6.0，12.0，24.0，35.2，48.0，54.0，60.0，66.0。

太阳天顶角(SZA)和MODTRAN运行的光谱幅亮度(mod_sr)的关系图如图2所示。

2、判定水平气象视距(VIS)的变化间距

在判定水平气象视距(VIS)的变化间距前，还是首先看下水平气象视距(VIS)和MODTRAN运行的光谱幅亮度(mod_sr)的关系图，如图3所示。

VIS最原始的点：0.4，1.4，2.6，6.1，13.1，27.2，80.0

VIS的第一次线性插值(取相邻点的中点)：0.9，2.0，4.35，9.6，20.15，53.6。

同SZA插值原理相同，插入点的精度在1％以内是可以接受的，对于大于1％的需要进行再次插值，后面插值过程同理。每次插值过程其插值结果如下：

第一次VIS中确定要插入的值是0.9，2，4.35，9.6，20.15，53.6，总的VIS就是0.4，0.9，1.4，2.0，2.6，4.35，6.1，9.6，13.1，20.15，27.2，53.6，80.0。(13个)

第二次VIS中确定要插入的值是0.65，1.15，1.7，2.3，3.475，5.225，7.85，11.35，16.625，23.675，40.4，66.8；对数据进行四舍五入，VIS中所有的值就是0.4，0.7，0.9，1.2，1.4，1.7，2.0，2.3，2.6，3.5，4.4，5.2，6.1，7.9，9.6，11.4，13.1，16.6，20.2，23.7，27.2，40.4，53.6，66.8，80.0。(25个)

第三次VIS中确定需要插入的点为：3.05，3.95，7.0，33.8；对数据进行四舍五入，VIS中所有的值就是0.4，0.7，0.9，1.2，1.4，1.7，2，2.3，2.6，3.1，3.5，4.0，4.4，5.2，6.1，7.0，7.9，9.6，11.4，13.1，16.6，20.2，23.7，27.2，33.8，40.4，53.6，66.8，80.0。(29个)

第四次VIS中确定需要插入的点为：无

因此最终的水平气象视距(VIS)的所有值是：0.4，0.7，0.9，1.2，1.4，1.7，2，2.3，2.6，3.1，3.5，4.0，4.4，5.2，6.1，7.0，7.9，9.6，11.4，13.1，16.6，20.2，23.7，27.2，33.8，40.4，53.6，66.8，80.0。(29个)

水平气象视距(VIS)和MODTRAN运行的光谱幅亮度(mod_sr)的关系图如图3所示。

3、气溶胶光学厚度与水平气象视距的转换

气溶胶光学厚度是描述气溶胶光学特性的重要参数之一，在辐射传输和大气校正中被广泛应用；但是，在大气辐射传输软件MODTRAN中却采用水平气象视距(VIS)作为其输入参数，因此在反演反射率的时候，需要将气溶胶光学厚度转化为水平气象视距。

在输入参数时，MODTRAN要求输入VIS，根据VIS和用户选择的季节、水汽量等，可以计算出整层大气的透过率，它同时也给出气溶胶的透过率。根据Lambert-Beer定律，可以得到

${\mathrm{\τ}}_a={\∫}_{Z1}^{Z2}N(z,\mathrm{VIS})\·\mathrm{EXT}\left(\mathrm{\λ}\right)m_a\left(z\right)\mathrm{dz}=-\mathrm{InTa}---\left(5\right)$

式中N(z，VIS)为与高度(z)和能见度(VIS)有关的气溶胶粒子密度系数，EXT(λ)为消光系数，EXT(λ)为消光系数，λ为波长，m_a(z)为大气质量，Ta为气溶胶引起的透过率。

这样通过得到不同VIS，不同季节下的气溶胶光学厚度，我们就可以得到它们之间的关系。由于水汽对两者关系的影响非常小，所以就没考虑水汽的变化。

最终的拟合结果如下：

春夏： $\frac{1}{t\left(550\right)}=0.1202V+0.29735---\left(6\right)$

秋冬： $\frac{1}{t\left(550\right)}=0.1419V+0.1377---\left(7\right)$

公式中Vt(550)表示550nm处的气溶胶光学厚度，V代表水平气象视距(km)。

4、表观反射率计算

通过表观辐亮度及反射率的反演，得到大气上界的表观反射率。

步骤一：表观辐亮度计算：

用定标系数将图像上原始DN值转换为大气层外反射的表观辐亮度L，其关系为：L_λ＝DN·gain+bias

式中L_λ为测量的光谱辐亮度；DN为记录的电信号值；gain为响应函数的斜率(通道增益)；bias为响应函数的截距。对于环境与灾害监测预报小卫星影像来说，定标系数也就是响应函数的斜率和截距。

步骤二：表观反射率算法

表观反射率的计算公式为：

式中ρ_λ为波段λ的表观反射率，L_λ为波段λ的光谱辐亮度，d为日地天文单位距离，ESUN_λ是波段λ处的大气上界太阳光谱辐照度，θ_s是太阳天顶角。

5、地表反射率算法

基于辐射传输理论，假设地面为朗伯面，则在传感器处接受到的表观辐亮度可以表示为：

$L_m=L_0+\frac{\mathrm{\ρ}}{1-\mathrm{s\ρ}}\·\frac{T{F}_d}{\mathrm{\π}}---\left(8\right)$

式中L₀表x示零地表反射时大气引起的程辐射，T表示地表到传感器的透过率，s为大气球形反照率，ρ为地表反射率，F_d地表下行通量。所有变量都与波段范围有关；其中，s是波长、大气光学性质、以及一系列位置参数，如海拔高度、传感器高度、传感器方位角、视角以及太阳高度角方位角等的函数；F_d、T由大气条件及下垫面几何条件所决定的系数，与地表反射率无关。

6、利用MODTRAN建立通用查找表

在利用MODTRAN建立通用查找表的过程中，自变量是大气模式、水平气象视距、太阳天顶角、卫星天顶角、地表反射率，因变量是光谱幅亮度。通过地表反射率算法知道，在一种大气参数(大气模式、水平气象视距、太阳天顶角、卫星天顶角)下，利用三个地表反射率和对应的三个光谱幅亮度可以求出公式(8)的三个参数；也就是一种大气参数对应一组参数。因此可以将查找表设置成自变量是大气模式、水平气象视距、太阳天顶角，因变量是L₀，s及TF_d/π。

通用查找表中自变量大气模式只有中纬度夏季大气和中纬度冬季大气两种情况，而水平气象视距、太阳天顶角只有部分数据，因此在计算过程中，有时需要对这两个变量进行线性插值，也就是二元变量的线性插值。

二元插值是对两个变量的函数z＝f(x，y)进行插值。x，y是两个独立变量，代表太阳天顶角和水平气象视距；z是一个因变量，代表L₀，s或TF_d/π。在插值过程中，如果x是查找表中已有的值，则只需要把z看成是y的函数，进行一元插值。如果x不等于查找表中已有的值，则首先在查找表中查找到与x值相邻的两个值，然后分别求得在这两个值处关于y变量一元插值的z值，最后对x值进行线性插值即可得到所需的z值。

插值判别精度的方法：

其中mod_sr是MODTRAN运行的光谱幅亮度，inter_sr是插值计算得到的光谱幅亮度。误差在1％以内是可以接受的，对于大于1％的需要进行再次的插值。

以下为MODTRAN简介：

MODTRAN(Moderate resolution atmospheric Transmittace andRadiance code)模型主要是对LOWTRAN 7模型的光谱分辨率进行了改进，它把光谱分辨率从20cm^-1减少到2cm^-1，发展了一种2cm^-1光谱分辨率的分子吸收的算法和更新了对分子吸收的气压温度关系的处理，同时维持LOWTRAN 7的基本程序和使用结构。重新处理的分子有水汽、二氧化碳、臭氧、一氧化二氮、一氧化碳、甲烷和氧气、一氧化氮、二氧化氮、氨气和硝酸。新的带模式参数仍是从HITRAN谱线参数汇编计算的，范围覆盖了0～17900cm^-1。而可见在紫外线这些较短的波长上，仍使用LOWTRAN7的20cm^-1的分辨率。在MODTRAN中，分子透过率的带参数在1cm^-1的光谱间隔上计算。

在程序处理上。MODTRAN的构造保持了对LOWTRAN的改动最小，MODTRAN作为附加中分辨率光谱计算能力的选择项而不干扰原LOWTRAN的执行。气溶胶模式、多次散射计算、用户定义的大气模式等都不改变。MODTRAN使用的分子带模式参数数据作为外部数据文件读入。

LOWTRAN 7模型是以20cm^-1的光谱分辨率的单参数带模式计算0cm^-1到50000cm^-1的大气透过率、大气背景辐射、单次散射的光谱辐射亮度、太阳直射辐射度。LOWTRAN7增加了多次散射的计算及新的带模式、臭氧和氧气在紫外波段的吸收参数。它提供了6种参考大气模式的温度、气压、密度的垂直廓线，H₂O、O₃、O₂、CO₂、CH₄、N₂O的混合比垂直廓线及其他13种微量气体的垂直廓线，城乡大气气溶胶、雾、沙尘、火山喷发物、云、雨廓线和辐射参量如消光系数、吸收系数、非对称因子的光谱分布、还包括地外太阳光谱。

MODTRAN和LOWTRAN在计算分子透过率上的区别有几点：LOWTRAN使用单参数带模式及分子密度的标度函数。MODTRAN使用3个与温度有关的参数：吸收系数、线密度和平均线宽。对于每种分子，其光谱区分为1cm^-1的间隔。将线中心在间隔内和间隔外的邻近谱线分别建模。间隔内的谱线采用Voigt线型积分得到。采用Curtis-Godson近似将多层的分层路径近似为等价的均匀路径。

MODTRAN的输入参数从类型上讲，共分五类：

第一为控制运行参数，如采用何种辐射传输程序，是否进行多次散射计算等。这些主要在CARD1中完成，CARD5提供了多重复计算的选项。

第二类参数为传感器的参数，如传感器的通道参数，观测的波束(波长范围)。其中CARD1A中有是否输入传感器通道响应函数的选项。在CARD1A3中输入通道响应函数的文件名。在CARD4中输入模拟计算的波长范围(涵盖所有波段，稍长)。

第三类参数为大气参数。其中大气类型通过CARD1中的选项确定。其它具体参数包括气溶胶主要通过CARD2来进行选择。

第四类为观测几何条件。在CARD1中有关于几何条件的选项，另外在CARD3中主要为几何参数的输入选项。他通过多种方式组合来实现几何参数的输入。可根据计算的方便进行选择。

第五类为地表参量。在CARD1提供了地表参数设定的初步选项，只好在CARD4根据CARD1中设定的参数对地表的参数进行具体设定。

所有的输入参数都通过CARD1进行控制，然后在由后续的CARD进行具体设置。

7、反演每个像元的地表反射率

对遥感影像进行大气校正，第一步是获取相应的大气参数，包括大气模式、水平气象视距、卫星天顶角、太阳天顶角，第二步是通过辐射传输模型模拟计算出L₀，s及

三个参数，然后将三个参数带入公式(8)即可推出地表反射率。

根据输入的大气参数，MODTRAN运行三次，这三次分别假设三个地表反射率(ρ＝0，0.5，0.8)，每次运行的结果是得到三个模拟的表观辐亮度L_m。

根据第一步得到的结果ρ和L_m，利用公式(8)建立三个以L₀，s及TF_d/π为变量的三元一次方程组；解方程组，求出一组与输入大气参数对应的L₀，s及TF_d/π。

输入不同的大气参数，根据第二步和第三步的方法，建立通用查找表。查找表自变量是大气参数(大气模式、太阳天顶角、卫星天顶角、水平气象视距)，因变量是地表反射率公式中的L₀，s及TF_d/π。

读取环境与灾害监测预报小卫星的CCD数据，如果是DN值，则通过定标系数将DN值图像转化为表观幅亮度图像。

读取气溶胶光学厚度(AOT)数据的每个像元，根据每个像元的气溶胶光学厚度值，找到与该大气参数对应的三个参数L₀，s及TF_d/π。

逐像元校正CCD或HSI数据，得到每个像元的地表反射率。

本发明还提供一种用于HJ-1A/B的地表反射率计算装置，其基于HJ-1A/B双星数据实现地表反射率反演。

对HJ-1A/B数据的地表反射率反演首先通过表观辐亮度及反射率的反演，得到大气上界的表观反射率；然后通过MODTRAN模拟得到不同气溶胶条件下的大气辐射传输特性并生成相应查找表，通过逐像元表观反射率及相应的气溶胶光学厚度条件查找得到对应地表反射率产品。

在生态参数反演模块输入HJ-1-A/B CCD影像数据及相关参数，可以得到植被指数、叶面积指数和植被覆盖度三个生态参数。

该地表反射率计算装置包括：表观反射率计算模块、大气参数反演模块、太阳及观测几何参数反演模块和地表反射率反演模块。

表观反射率计算模块用来通过表观辐亮度及反射率的反演，得到大气上界的表观反射率。

大气参数反演模块基于HJ-1-A/B CCD 30米分辨率DN值数据和HJ-1-A/B CCD XML文件，利用暗目标法和不变目标法，基于经典的大气辐射传输理论，进行气溶胶光学厚度反演，得出地表反射率ρ与空间观测星载反射率(也称表观反射率)ρ’之间的关系。具体过程如本发明的方法中所述，在此不再赘述。

太阳及观测几何参数反演模块基于HJ-1-A/B CCD 30米分辨率DN值数据和HJ-1-A/B CCD XML文件，计算得到具有地理参考的太阳天顶角数据、太阳方位角数据、观测天顶角数据和观测方位角数据。具体过程如前面本发明的方法中所述，在此不再赘述。

地表反射率反演模块基于HJ-1-A/B CCD、HSI全波段Radiance数据或DN值和水平气象视距，太阳天顶角等大气参数，利用大气辐射传输模型MODTRAN来建立通用查找表，实现逐像元大气校正，得到CCD，HSI波段的地表反射率数据。具体过程如前面本发明的方法中所述，在此不再赘述。

所述地表反射率计算装置还可以包括生态参数反演模块，其实现植被指数、叶面积指数及植被覆盖度三个生态参数的反演：(1)通过输入HJ-1A/B CCD蓝光、红光和近红外地表反射率数据和高光谱反射率数据计算得到植被指数；(2)通过土地利用/土地覆盖数据和HJ-1A/HJ-B CCD第3、第4地表反射率数据计算得到叶面积指数LAI；(3)运用叶面积指数和植被指数反演得到植被覆盖度。

植被指数反演算法：对于环境与灾害监测预报小卫星，可以使用宽覆盖CCD相机的第一、三、四波段，来计算归一化及增强植被指数，公式如下：

NDVI＝(CCD4-CCD3)/(CCD4+CCD3)

EVI＝2.5×(CCD4-CCD3)/(CCD4+6.0×CCD3-7.5×CCD1+1.0)

对于环境与灾害监测预报小卫星高光谱载荷，选取中心波长分别位于近红外和红光的谱段进行归一化植被指数计算：

Hyp_NDVI＝(Hyp_Nir-Hyp_Red)/(Hyp_Nir+Hyp_Red)

EVI计算所需三个经验系数，即：大气修正参数C1(6.0)、C2(7.5)，背景调节参数L(1.0)，分别用于减少土壤背景噪声的影响及通过蓝光波段来修正大气对红光波段的影响。

叶面积指数反演算法：对叶面积指数反演采用了基于查找表和基于经验公式的反演方法，其中基于查找表方法是主反演程序，而基于经验公式的反演方法是备用反演程序。

基于查找表的方法首先需要根据遥感机理模型生成一个查找表，然后用如下代价函数反演LAI：

cost＝∑[(NIR_LUP-NIR)²+(RED_LUP-RED)²]

其中NIR_LUP和RED_LUP分别为查找表中的近红外和红色波段反射率，NIR和RED分别是HJ-1A/B卫星的近红外和红色波段反射率。cost是一个代价函数。反演LAI是在查找表中寻找使cost最小的反射率对应的LAI数值，然后将所有LAI取平均得到LAI的反演值。

对不同的地物类型选用不同的查找表，如果查找表方法反演失败，则启动备用反演程序，用基于经验的方法反演LAI，即：

LAI＝f(NDVI，LandCover)

其中f为一个根据遥感机理模型和地物类型拟合的经验公式。在具体实现时同样采用查找表实现反演，即建立包括一系列生态系统类型、NDVI和LAI对应数据的数据表。

使用基于反射率查找表反演失败的标准主要针对以下情况：

1.对于不在6类生态系统类型(草地和谷类作物、灌木、阔叶作物、稀疏草原、阔叶森林、针叶森林)的地物类型，需要用上述经验公式处理。

2.对于由于云或其它原因造成像素反演质量控制值为2和3的情况，采用经验公式处理。

植被覆盖度反演：用来估算植被覆盖度的公式表示如下：

$f=\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{\mathrm{NDV}{I}_g-{\mathrm{NDVI}}_{\min }}$

NDVI_g＝NDVI_max-(NDVI_max-NDVI_min)exp(-K·LAI)

式中，NDVI_max代表完全植被覆盖区域的NDVI值，即纯植被像元的NDVI值。NDVI_min代表裸露地表覆盖区域的NDVI值，即纯裸土像元的NDVI值。K为消光系数，这里K＝1。

通过本发明的方法和装置，利用MODTRAN对大气-地表-遥感器之间的辐射传输过程进行模拟，依据卫星同步反演得到的大气参数和幅亮度、针对环境减灾小卫星传感器的特性，能够生成适用于环境减灾小卫星数据的查找表，并最终反演出地表反射率。构建一个完整的环境减灾小卫星数据的大气校正模块。用本算法形成的大气校正模块能有效的消除大气对遥感影像的吸收和散射，恢复地表目标的反射率，解决了行业应用的瓶颈，进一步扩展环境减灾小卫星数据的应用范围。

Claims (14)

1.一种用于HJ-1A/B卫星的地表反射率计算方法，其特征在于，包括步骤：

1)基于暗目标法或不变目标方法得到大气气溶胶光学厚度；

2)基于HJ-1A/B XML文件和影像数据计算每日周期的太阳天顶角、太阳方位角、观测天顶角和观测方位角；

3)基于上述步骤1)和2)中得到的参数，利用MODTRAN模拟星上幅亮度L_m；

4)通过步骤3)建立地表反射率查找表；

5)根据大气参数和待校正图像，利用查找表实现大气校正。

2.根据权利要求1所述的地表反射率计算方法，其特征在于，在步骤3)之前，还包括将气溶胶光学厚度转换为水平气象视距。

3.根据权利要求1所述的地表反射率计算方法，其特征在于，在步骤3)之前，还包括通过表观辐亮度及反射率的反演，得到大气上界的表观反射率。

4.根据权利要求1所述的地表反射率计算方法，其特征在于，在步骤4)中，查找表的自变量是大气模式、水平气象视距、太阳天顶角、卫星天顶角，因变量是求地表反射率的三个参数。

5.根据权利要求4所述的地表反射率计算方法，其特征在于，还包括步骤6)：参考MODIS大气校正算法中水平气象视距、太阳天顶角这两个变量的变化间距，然后根据MODTRAN运行结果来判定水平气象视距、太阳天顶角的变化间距上是否符合精度，若不符合精度，则进行插值。

6.根据权利要求1所述的地表反射率计算方法，其特征在于，通过输入HJ-1A/B CCD蓝光、红光和近红外地表反射率数据和高光谱反射率数据，实现植被指数的反演。 

7.根据权利要求1所述的地表反射率计算方法，其特征在于，所述生态参数反演模块通过植被指数和叶面积指数LAI，实现植被覆盖度的反演。

8.一种用于HJ-1A/B卫星的地表反射率计算装置，其特征在于，包括：

大气参数反演模块，其基于暗目标法或不变目标方法得到大气气溶胶光学厚度；

太阳及观测几何参数反演模块，其基于HJ-1A/B XML文件和影像数据计算每日周期的太阳天顶角、太阳方位角、观测天顶角和观测方位角；

地表反射率反演模块，其利用大气参数反演模块和太阳及观测几何参数反演模块计算得到的参数，利用MODTRAN模拟星上幅亮度L_m；并建立地表反射率查找表；然后根据大气参数和待校正图像，利用查找表实现大气校正。

9.根据权利要求8所述的地表反射率计算装置，其特征在于，还包括生态参数反演模块，其通过输入HJ-1A/B CCD蓝光、红光和近红外地表反射率数据和高光谱反射率数据，实现植被指数的反演。

10.根据权利要求9所述的地表反射率计算装置，其特征在于，所述生态参数反演模块通过土地利用/土地覆盖数据和HJ-1A/HJ-BCCD第3、第4地表反射率数据，实现叶面积指数LAI的反演。

11.根据权利要求10所述的地表反射率计算装置，其特征在于，所述生态参数反演模块通过植被指数和叶面积指数LAI，实现植被覆盖度的反演。

12.根据权利要求8所述的地表反射率计算装置，其特征在于，所述地表反射率反演模块还将气溶胶光学厚度转换为水平气象视距。 

13.根据权利要求12所述的地表反射率计算装置，其特征在于，所述查找表的自变量是大气模式、水平气象视距、太阳天顶角、卫星天顶角，因变量是求地表反射率的三个参数。

14.根据权利要求13所述的地表反射率计算装置，其特征在于，所述地表反射率反演模块参考MODIS大气校正算法中水平气象视距、太阳天顶角这两个变量的变化间距，然后根据MODTRAN运行结果来判定水平气象视距、太阳天顶角的变化间距上是否符合精度，若不符合精度，则进行插值。 

Images