CN101915914A - 一种基于查找表的遥感影像逐像元大气校正方法 - Google Patents

Abstract

本发明公布了一种基于查找表的遥感影像逐像元大气校正方法，本发明方法利用6S大气辐射传输模型离线计算建立了不同气溶胶光学厚度、太阳天顶角、传感器天顶角以及地表海拔条件下大气校正系数的查找表，基于该查找表对MODIS遥感影像进行逐像元大气校正。本发明可以有效改善由于大气条件、传感器位置等空间分布差异对MODIS图像大气校正的影响。

Description

一种基于查找表的遥感影像逐像元大气校正方法

技术领域

本发明涉及一种遥感影像的大气校正方法，具体涉及一种基于查找表的遥感影像逐像元大气校正方法。

背景技术

随着遥感理论和技术的不断发展及其在应用中的不断完善，遥感已经从开始的定性解释发展到了现在的定量分析。大多数地表参数的定量遥感模型如反照率、叶面积指数等等都是建立在地表反射率的基础上，通过大气校正获取精确的地表反射率信息是遥感定量分析的关键问题，也是主要难点之一。随着定量遥感技术迅速发展，特别是利用多传感器、多时相遥感数据进行土地利用和土地覆盖变化监测、全球资源环境分析、气候变化监测等的需要，使得遥感图像大气校正方法的研究越来越受到重视，对遥感图像大气校正的研究具有重大意义。

太阳辐射从穿过大气层到达地表然后被反射再次穿过大气层到达遥感传感器的过程中，受到大气分子、气溶胶和云粒子等大气成份吸收与散射的影响，使其获取的遥感信息中带有一定的非目标地物的成像信息，因此需要进行大气校正去除遥感图像中的大气影响。大气对太阳辐射的影响主要来源于大气的吸收、大气分子的瑞利散射和气溶胶的米散射作用，大气校正的目标就在于从遥感传感器观测到的辐射信号中去除这些作用的影响，得到地面目标物的真实反射率。

目前大气校正的方法有暗目标像元法、地表实测线性回归法、大气辐射传输模型法等等。其中辐射传输模型法由于普适性好、精度较高在近些年来得到广泛应用。辐射传输模型法的基本原理是通过大气状况、卫星和太阳位置信息等输入参数，利用合适的大气辐射传输模型计算出大气校正参数的数值解(如大气透过率、程辐射等)，以此为依据计算地表真实反射率。目前已有多位学者利用此方法对遥感影像进行了大气校正处理。由于前人的研究主要针对范围较小的TM、ASTER、CBERS等遥感影像，考虑的变化因素相对较少，或者把整幅影像的输入参数看成均一，或者仅仅考虑气溶胶光学厚度的变化。也有学者利用大气辐射传输模型对MODIS影像进行了大气校正，但是所选研究区范围都不大，并非针对全景影像，也都是采用了统一的输入参数进行大气校正。但是对于大范围的遥感影像而言，各个像元的差异比较大，用统一的输入参数进行校正无法考虑各个像元大气状况、观测角度等条件的差异，但是如果每个像元都在线运行辐射传输模型计算将耗费大量的计算机资源和时间，如何在保证精度的条件下快速实现遥感影像的逐像元大气校正成为迫切需要解决的问题。

发明内容

本发明目的是针对现有技术存在的缺陷，在分析大气校正结果对于大气气溶胶光学厚度、卫星和太阳位置、海拔高度输入参数变化敏感性的基础上，通过6S大气辐射传输模型离线计算建立了不同条件下大气校正系数的查找表，基于查找表对MODIS影像进行逐像元的大气校正。本发明的校正结果与利用统一输入参数校正结果进行比较分析表明，本发明的方法具有更高的精度和实用性。

本发明为实现上述目的，采用如下技术方案：

本发明一种基于查找表的遥感影像逐像元大气校正方法包括如下步骤：

(一)离线模拟

(1)采用6S模型进行大气校正，模拟太阳辐射在大气中的传输状况得到地表真实反射率ρ_s，并给出大气校正系数xa、xb和xc：

$\left\{\begin{array}{c}{{\mathrm{\ρ}}^{\′}}_s=\frac{\mathrm{\πL}}{T_g{E}_s\cos \left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)}-\frac{{\mathrm{\ρ}}_{R+A}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)}{T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)}\\ {\mathrm{\ρ}}_s=\frac{{{\mathrm{\ρ}}^{\′}}_s}{1+S{{\mathrm{\ρ}}^{\′}}_s}\end{array}\right.,$

$\mathrm{xa}=\frac{\mathrm{\π}}{T_g{E}_s\cos \left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)},$ $\mathrm{xb}=\frac{{\mathrm{\ρ}}_{R+A}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)}{T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)},$ xc＝S；

其中，T_g为大气中O₃、H₂O、N₂、CO₂、O₂、CH₄气体对太阳辐射的吸收透过率，ρ_R+A为大气层辐射反射率即气溶胶程辐射与瑞利散射程辐射之和，θ_s、θ_v分别为太阳天顶角和传感器天顶角，T(θ_s)、T(θ_v)分别为太阳至地面、地面至传感器的大气路径透过率，S为大气底层向下的半球反射率，L为表观辐亮度，E_s为对应波长的大气层顶太阳辐照度；

(2)敏感性分析遥感影像，只考虑气溶胶光学厚度以及传感器天顶角的变化，太阳天顶角及地表海拔设为常量；

(3)利用6S模型根据步骤(1)所述的大气校正系数xa、xb、xc公式离线计算不同气溶胶光学厚度、太阳天顶角、传感器天顶角以及地表海拔情况下的大气校正系数xa、xb、xc，并建立针对气溶胶光学厚度、太阳天顶角、传感器天顶角以及地表海拔相对于大气校正系数xa、xb、xc的查找表；

(二)在线校正

(4)对遥感影像进行大气校正时，从遥感影像数据中导出表观辐亮度、太阳天顶角、传感器天顶角数据，从遥感影像的气溶胶光学厚度中读出光学厚度数据，结合1km分辨率的SRTM DEM数据，逐像元根据步骤(3)所述的查找表进行线性插值，计算出当前像元对应波段的大气校正系数xa、xb、xc。

本发明具有如下优点：

(1)本发明依据MODIS图像自身的信息，不需要地表实测数据，对于输入参数的依赖性较小。

(2)本发明的校正方法可以有效改善由于大气条件、传感器位置等空间分布差异对MODIS影像大气校正的影响。

(3)本发明相对于6S辐射传输模型在线校正法而言，精度相近，速度要快捷很多；相对于统一输入参数法而言，精度有明显提高，能够体现不同像元大气性质、观测条件和地表海拔的空间分布差异。

附图说明

图1：输入参数变化对地表反射率校正结果的影响；

图2：基于查找表的MODIS逐像元大气校正流程图；

图3：扫描线上各像元的气溶胶光学厚度、太阳天顶角、传感器天顶角、地表海拔；

图4：MODIS第1～4波段三种不同方法校正结果的比较；

具体实施方式

下面结合附图对本发明进一步详述：

本发明一种基于查找表的遥感影像逐像元大气校正方法，包括如下步骤：

1)对辐射传输方程进行分析推导

6S(Second Simulation of the Satellite Signal in the Solar Spectrum radiative code)模型是目前世界上发展得比较完善的大气校正模型之一，由Vermote等人在5S(the Simulationofthe Satellite Signal in the Solar Spectrum radiative code)模型的基础上改进而来，适用于0.25～4μm波长范围内电磁波的大气辐射传输模拟。利用6S模型进行大气校正的工作流程是：将大气参数、几何参数、观测波段和海拔高度等条件输入6S模型之后，通过计算模拟太阳辐射在大气中的传输状况，并给出大气校正系数xa、xb和xc。根据下面的公式就可以计算得到经过大气校正的地表真实反射率。

y＝xa*(measured radiance)-xb，acr＝y/(1.+xc*y)        (1)

公式1虽然形式简单，但是并不是一个单纯的经验公式。通过对辐射传输方程进行推导和变换可知，系数xa、xb、xc和方程是有着具体物理意义的。

在6S模型中，表观反射率可以描述为：

${\mathrm{\ρ}}_{\mathrm{TOA}}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)=T_g[{\mathrm{\ρ}}_{R+A}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)+T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)\frac{{\mathrm{\ρ}}_s}{1-S{\mathrm{\ρ}}_s}]---\left(2\right)$

式中，ρ_TOA为表观反射率，T_g为大气中O₃、H₂O等气体对太阳辐射的吸收透过率，ρ_R+A为大气程辐射反射率(气溶胶程辐射与瑞利散射程辐射之和)，θ_s、θ_v分别为太阳天顶角和卫星天顶角，T(θ_s)、T(θ_v)分别为太阳-地面、地面-传感器的大气路径透过率，S为大气底层向下的半球反射率，ρ_s为地表真实反射率。

上式可改写为：

${\mathrm{\ρ}}_s=\frac{\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)}{T_g}-{\mathrm{\ρ}}_{R+A({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)}}{T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)+[\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)}{T_g}-{\mathrm{\ρ}}_{R+A}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)]S}---\left(3\right)$

其中表观反射率ρ_TOA＝πL/E_scos(θ_s)，式中，L为表观辐亮度，E_s为对应波长的大气层顶太阳辐照度。另外，再设

则上式可改写为

$\left\{\begin{array}{c}{{\mathrm{\ρ}}^{\′}}_s=\frac{\mathrm{\πL}}{T_g{E}_s\cos \left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)}-\frac{{\mathrm{\ρ}}_{R+A}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)}{T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)}\\ {\mathrm{\ρ}}_s=\frac{{{\mathrm{\ρ}}^{\′}}_s}{1+S{{\mathrm{\ρ}}^{\′}}_s}\end{array}\right.---\left(4\right)$

将公式1和4对照，发现

xc＝S。xa、xb、xc三个系数与大气状况、太阳与传感器的位置、地面的海拔、观测波段等有关，而与地表反射特性无关。在其它条件相同的情况下改变目标的表观反射率或者辐亮度，xa、xb、xc三个校正系数不发生改变。

2)敏感性分析

通常在利用6S模型进行大气校正时，往往假设整个研究区大气条件均一、太阳和传感器的观测角度一致、地表具有相同的海拔高度，然后利用6S模型计算出大气校正系数，并应用于整个图像，得到整个研究区的地表真实反射率。对于小范围的遥感影像(如TM或者裁切为小块的MODIS数据)而言，这种方法是可行的，但是对于大范围的遥感影像(如整幅MODIS或者AVHRR数据)而言，这样的简单假设就显得不合理了，势必会影响校正得到的地表反射率精度。

以MODIS数据为例，一景影像范围超过400万平方公里，各个像元的大气条件、太阳位置、传感器位置、海拔高度等均不相同，甚至有较大差异。本发明首先分析了6S模型地表反射率校正结果对于这些参数的敏感性，从而更深层次的了解各种参数对于大气校正结果的影响程度，为后面在大气校正过程中的参数选取提供参考。本发明选择了大气气溶胶光学厚度、太阳天顶角、传感器天顶角以及地表海拔这4个参数来分析这些参数的变化对6S模型大气校正结果的影响。

6S模型的初始输入参数为：太阳天顶角设为30度，方位角为0度，卫星天顶角设为30度，方位角为90度，成像日期7月1日，大气模式为中纬度夏季，气溶胶类型选择大陆型气溶胶，550nm光学厚度为0.5，地表海拔设为1km，波段选用MODIS的第1、2、3、4波段，表观反射率设为0.3。然后分别计算气溶胶光学厚度在0～1之间变化、太阳天顶角在5～35度之间变化、传感器天顶角在0～65度之间变化、海拔高度在0～4000m之间变化对于校正得到的地表反射率的影响，结果如图1所示。

从图1中可以看出，气溶胶光学厚度对于大气校正的影响明显要高于其它因素，传感器天顶角对于大气校正的影响也比较大，太阳天顶角以及海拔高度的变化对于6S模型大气校正的影响相对较小。如果对于大气校正精度要求不高的话，可以只考虑气溶胶光学厚度以及传感器天顶角的变化，其它要素设为常量，这样查找表比较简单而且运算速度大大提高。

3)建立查找表

利用6S模型离线计算不同气溶胶光学厚度、太阳天顶角、传感器天顶角以及地表海拔情况下的校正系数xa、xb、xc，建立针对这4个因素的查找表。6S模型的初始参数设置基本同前面敏感性分析时的输入参数，不过气溶胶光学厚度、太阳天顶角、传感器天顶角和海拔的取值范围有所变化。在敏感性分析中4个参数取的都是比较常见的值，而在建立查找表时还要考虑一些不常见的情况，所以此处参数取值范围要大于前面敏感性分析时的取值。光学厚度取值0～2，在0到1之间以0.1步长变化，在1～2之间以0.2步长变化共16个值；太阳天顶角从0度至80度以10度的步长变化，共有9个值；感器天顶角从0度至70度以10度的步长变化，共有8个值；地表海拔选取0m、200m、500m、1000m、2000m和4000m共6个值。将上述参数组合代入6S模型循环计算，得到MODIS前7个波段在不同气溶胶光学厚度、太阳方位角、传感器方位角、地表海拔下的大气校正系数xa、xb、xc查找表。

4)逐像元校正方法

对MODIS影像进行大气校正时，从MODIS L1B数据中导出表观辐亮度、太阳天顶角、传感器天顶角数据，从MODIS的气溶胶光学厚度产品中读出光学厚度数据，结合1km分辨率的SRTM DEM数据，逐像元根据查找表进行线性插值，计算出当前像元对应波段的大气校正系数xa、xb、xc，并应用公式1计算该像元的地表反射率。整个校正流程见图2。

5)实例验证

MODIS影像范围很大，很难像TM大气校正那样在卫星过境同时进行地表同步观测来对校正结果进行验证。前人的研究成果已经表明，在大气参数、几何位置信息等输入参数较为准确的情况下，小范围的遥感影像经过6S模型的大气校正之后得到的地表反射率信息比较准确。本发明主要验证的是在比较大的范围内，利用查找表逐像元进行大气校正与整景影像利用统一输入参数进行大气校正的差别。

为了对本方法的大气校正结果进行验证，本发明对比了利用6S模型逐像元在线计算、基于查找表逐像元校正以及整景影像利用统一输入参数的6S模型校正方法对于MODISL1B影像的校正结果。为了陈述方便，三种方法分别简称为Online法(6S onlinecorrection)、LUT法(Look-up table method)和PU法(parameter-uniform method)。因为6S模型在线计算耗时太长，本文取了一景2005年10月17日的AQUA/MODIS L1B数据的一个扫描行(2137个像元)来进行验证。

图3给出了扫描线上自西向东各个像元的气溶胶光学厚度、太阳天顶角、传感器天顶角和地表海拔。从图上可以看出，气溶胶变化比较剧烈，有些像元低于0.2，有些像元高于0.6；太阳天顶角自西向东逐渐增大；而传感器天顶角则从星下点开始向两边逐渐增大；海拔高度总体上来说西高东低，西部总体上在1500m至3000m之间，东部海拔一般低于50m，为平原区域。

6S在线逐像元校正过程中，各个像元的气溶胶光学厚度数据、太阳天顶角、传感器天顶角和地表海拔分别来自各像元对应的数据。整景影像利用统一输入参数的6S模型校正过程中，气溶胶光学厚度取中心像元值0.15，太阳天顶角取值49，传感器天顶角取值0，地表海拔取值30m。

三种不同校正方法得到的MODIS第1、2、3、4波段地表反射率见图4，从图中可以看出，第1、3和4波段的LUT法和Online法校正结果很接近，而P-U法到的地表反射率与前两者的差别要大一些，第2波段(近红外波段)受到气溶胶散射影响相对较小，LUT法、Online法和P-U法这三种方法计算得到的地表反射率比较接近。

以Online法的校正计算结果为基准计算了LUT法和P-U法计算得到的第1～4波段反射率的均方根误差RMSE(表1)来定量评价这两种方法。从表中也看出，除了第2波段两种方法的误差相差较小之外，其余波段查找表法的校正精度要显著高于统一参数法，表明了本发明大气校正结果相对于常用的统一参数法校正结果的优越性。

表1LUT法和P-U法校正结果的均方根误差

Claims (1)

1.一种基于查找表的遥感影像逐像元大气校正方法，其特征在于包括如下步骤：

(一)离线模拟

(1)采用6S模型进行大气校正，模拟太阳辐射在大气中的传输状况得到地表真实反射率ρ_s，并给出大气校正系数xa、xb和xc：

$\left\{\begin{array}{c}{{\mathrm{\ρ}}^{\′}}_s=\frac{\mathrm{\πL}}{T_g{E}_s\cos \left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)}-\frac{{\mathrm{\ρ}}_{R+A}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)}{T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)}\\ {\mathrm{\ρ}}_s=\frac{{{\mathrm{\ρ}}^{\′}}_s}{1+S{{\mathrm{\ρ}}^{\′}}_s}\end{array}\right.,$

$\mathrm{xa}=\frac{\mathrm{\π}}{T_g{E}_s\cos \left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)},$ $\mathrm{xb}=\frac{{\mathrm{\ρ}}_{R+A}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)}{T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)},$ xc＝S；

其中，T_g为大气中O₃、H₂O、N₂、CO₂、O₂、CH₄气体对太阳辐射的吸收透过率，ρ_R+A为大气层辐射反射率即气溶胶程辐射与瑞利散射程辐射之和，θ_s、θ_v分别为太阳天顶角和传感器天顶角，T(θ_s)、T(θ_v)分别为太阳至地面、地面至传感器的大气路径透过率，S为大气底层向下的半球反射率，L为表观辐亮度，E_s为对应波长的大气层顶太阳辐照度；

(2)敏感性分析遥感影像，只考虑气溶胶光学厚度以及传感器天顶角的变化，太阳天顶角及地表海拔设为常量；

(3)利用6S模型根据步骤(1)所述的大气校正系数xa、xb、xc公式离线计算不同气溶胶光学厚度、太阳天顶角、传感器天顶角以及地表海拔情况下的大气校正系数xa、xb、xc，并建立针对气溶胶光学厚度、太阳天顶角、传感器天顶角以及地表海拔相对于大气校正系数xa、xb、xc的查找表；

(二)在线校正

(4)对遥感影像进行大气校正时，从遥感影像数据中导出表观辐亮度、太阳天顶角、传感器天顶角数据，从遥感影像的气溶胶光学厚度中读出光学厚度数据，结合1km分辨率的SRTM DEM数据，逐像元根据步骤(3)所述的查找表进行线性插值，计算出当前像元对应波段的大气校正系数xa、xb、xc。

Images