CN102778675B - 一种卫星遥感影像大气校正方法及其模块 - Google Patents

Abstract

本发明公开了一种卫星遥感影像大气校正方法及其模块，利用6S（Second Simulation of a Satellite Signal in the Solar Spectrum）辐射传输模型建立大气校正参数查找表LUT（Look Up Table），利用消息传递编程模型MPI（Message Passing Interface）的主从并行程序设计模式，设计数据块驱动与功能链驱动的并行大气校正模型，实现基于6S辐射传输查找表的并行大气校正，有效解决卫星遥感影像大气校正的时间和精度问题。本发明能够适用于各种卫星遥感影像的业务化大气校正。

Description

一种卫星遥感影像大气校正方法及其模块

技术领域

本发明涉及卫星遥感数据处理，特别是指一种大气校正方法。

背景技术

遥感影像大气校正是进行遥感定量反演的前提和关键，大气校正的精度与速度直接影响到遥感实际应用。

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

目前大气校正方法主要分为两种：绝对大气校正、相对大气校正。绝对大气校正方法有：辐射传输模型法，包括FLAASH（Fast Line of SightAtmospheric Analysis of Spectral Hypercube）模型、MORTRAN（ModerateResolution Transmission）模型、LOWTRAN（Low Resolution Transmission）模型、ATCOR（Atmospheric and Topographic Correction for Airborne ScannerData）模型、6S（Second Simulation of a Satellite Signal in the Solar Spectrum）模型等；统计学模型法（平场域、对数残差、内部平均、经验线性等）；黑暗像元法等。相对大气校正方法有：不变目标法，直方图匹配法，参考值大气校正法，大气抗阻植被指数法等。辐射传输模型法由于普适性好、精度较高，得到了广泛应用。其中，6S辐射传输模型考虑了气体吸收、分子和气溶胶散射，预先设置了50多种波段模型，可以很好地模拟太阳光在太阳-地面目标-传感器的传输过程中所受到的大气影响。然而，对于幅宽较小的遥感影像，大气影响因素相对均一，随着影像幅宽的增大、分辨率的提高，各个像元的差异比较大，此时需要考虑各个像元大气状况、观测角度等条件的差异，创建大气校正参数查找表，逐像元进行大气校正成为研究的热点。

随着面向遥感应用需求的拓展、深化，以及自主数据接受系统的建立，遥感信息应用所日常接收和处理的遥感影像数据量剧增，单机设备处理海量数据速度慢，难以满足业务化应用的需求。如何在保证大气校正精度的条件下，快速实现遥感影像的逐像元大气校正成为迫切需要解决的问题。

高性能集群技术的发展为遥感影像快速大气校正提供了有效的手段。集群计算机系统由一些相互独立的计算机（节点）高度紧密地协作完成计算工作，其处理能力、存储能力、管理能力、读写能力、可靠性都得到了大幅的提升，可有效缓解单机处理速度难以满足海量遥感影像数据处理时效要求的矛盾，实现海量遥感数据的高效处理。

因此，创建大气校正参数查找表，利用高性能计算技术和海量存储管理技术，快速实现遥感影像的逐像元精确大气校正，是当前遥感影像大气校正发展的一种重要趋势。

发明内容

有鉴于此，本发明的目的在于提出一种针对卫星遥感影像大气校正方法以及运用该方法完成大气校正的模块。通过本发明所提供的大气校正方法，能够快速实现卫星遥感影像的逐像元精确大气校正，克服传统大气校正精度相对低、速度慢等问题。

基于上述目的本发明提供的卫星遥感影像大气校正方法，至少包括以下步骤：

根据不同的气象条件参数，利用6S辐射传输模型，创建由大气校正参数组成的大气校正参数查找表；

根据各个波段的增益与偏移值将影像的像素值转为亮度值；

根据太阳天顶角、大气外部的太阳辐照度、地球太阳距离系数将亮度值转为大气顶层表观反射率；

根据像元所处的波段、太阳天顶角、卫星天顶角、相对方位角，查找大气校正参数查找表表中对应的大气校正参数，利用气溶胶光学厚度对该数据进行插值计算得到像元处的大气校正参数，根据此大气校正参数将大气层顶表观反射率转为地表反射率。

优选地，所述大气校正参数查找表的建立包括如下步骤，

根据研究区域的地理气象条件以及影像的获取时间，确定大气环境参数；利用6S辐射传输模型，创建不同大气校正输入参数等条件下的大气校正参数查找表；所述大气校正参数查找表的输出参数包括内部大气反射率、总透过率、球面反照率。

可选的，所述大气环境参数包括大气模式、气溶胶模式、气压、水汽含量、臭氧含量、目标海拔高度；所述大气校正输入参数包括气溶胶光学厚度、太阳天顶角、卫星天顶角、太阳与卫星相对方位角。

所述影像亮度值的计算使用如下公式：

$L_{\mathrm{TOA}}^k=\frac{X^k}{{\mathrm{\α}}_k}+{\mathrm{\β}}_k---\left(5\right)$

其中：

为亮度（单位：Wm^-2sr^-1um^-1），X^k为输入影像的像素值，α_k为每个波段的增益，β_k为每个波段的偏移；

所述大气层顶表观反射率的计算使用如下公式：

${\mathrm{\ρ}}_{\mathrm{TOA}}^k=\frac{\mathrm{\π}{L}_{\mathrm{TOA}}^k}{E_S^k\·\cos \left({\mathrm{\θ}}_s\right).d/d_0}---\left(6\right)$

其中，

为大气层顶表观反射率，θ_s为太阳天顶角，

为离地球d₀量测的大气外部的太阳辐照度，d/d₀为地球太阳距离系数；

所述地表反射率的计算依据如下公式：

${\mathrm{\ρ}}_S^{\mathrm{unif}}=\frac{A}{1+\mathrm{SA}}---\left(7\right)$

$A=\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}-{\mathrm{\ρ}}_{\mathrm{atm}}}{T\left(\mathrm{\μS}\right).T\left(\mathrm{\μV}\right).t_g^{\mathrm{al\; lg\; as}}}---\left(8\right)$

其中，

为假设的朗伯面和均一环境下的地表反射率，ρ_TOA为大气顶部的反射率，ρ_atm为内部大气反射率，为大气球面反照率（总大气传输），T(μS)为下行散射，T(μV)为上行散射，S为球面反照率。

优选地，所述卫星遥感影像大气校正方法采用数据块驱动与功能链驱动的并行大气校正模型，且基于消息传递编程模型的主从并行程序设计模式来执行大气校正。所述主从并行的模式包括：

主进程分配数据及参数，所述主进程根据总进程数和影像大小，计算得到每个从进程分配到的影像块在整幅影像中的位置，将所述影像块及配置参数传递给所述从进程；

所述从进程读取大气校正参数，并进行参数配置，所读取的大气参数为：增益偏移值、辐照度、大气校正参数查找表；

所述从进程根据主进程传递过来的位置信息，读取相应数据块，进行大气校正；

主进程不断接收从进程发送的结果影像块，并将其写入输出文件。

进一步，本发明提出一种卫星遥感影像大气校正模块，包括如下子模块：

查找模块、计算模块、输入模块和输出模块；

所述查找模块根据不同的气象条件参数，利用6S辐射传输模型，创建由大气校正参数组成的大气校正参数查找表；

所述输入模块用于输入大气校正所需参数；

所述计算模块根据各个波段的增益与偏移值将影像的像素值转为亮度值；根据太阳天顶角、大气外部的太阳辐照度、地球太阳距离系数将亮度值转为大气顶层表观反射率；根据像元所处的波段、太阳天顶角、卫星天顶角、相对方位角，查找大气校正参数查找表表中对应的大气校正参数，利用气溶胶光学厚度对该数据进行插值计算得到像元处的大气校正参数，根据此大气校正参数将大气顶层表观反射率转为地表反射率；

所述输出模块输出大气校正结果。

所述计算模块计算影像亮度值时使用公式（5）；

所述计算模块计算大气层顶表观反射率时使用公式(6)；

所述计算模块计算地表反射率时使用公式（8）。

从上面所述可以看出，对比现有利用均一输入参数校正，本发明提供的卫星遥感影像大气校正方法具有更高的精度和实用性；本发明所采用的环境与单机处理环境相比，速度大幅度提高，能够有效解决卫星遥感影像大气校正的时间和精度问题，适用于各种卫星遥感影像的业务化大气校正。

具体而言，本发明所提供的卫星遥感影像大气校正方法针对各个像元大气状况、观测角度等条件的不同，利用6S辐射传输模型，创建大气校正参数查找表，基于查找表逐像元进行大气校正，相对于均一输入参数而言，大气校正精度明显提高。

同时，本发明面向遥感业务化应用需求，利用高性能计算技术和海量存储管理技术，建立块驱动与功能链驱动的并行大气校正模型，在保证精度的条件下，快速实现逐像元大气校正，大气校正速度显著提高。

附图说明

图1为本发明所采用的基于消息传递编程模型的并行大气校正流程示意图；

图2为本发明的卫星遥感影像大气校正方法流程图；

图3为本发明实施例的大气校正参数查找表建立流程示意图；

图4为本发明实施例的并行大气校正流程图示意图；

图5为本发明实施例的光谱响应函数曲线图之一；

图6为不同大气校正方法植被样点平均反射率的对比线图；

图7为不同的大气校正方法随机点归一化差分植被指数对比线图；

图8为本发明实施例的大气校正并行执行时间线图；

图9为本发明实施例的大气校正并行执行加速比线图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

为了提高大气校正的精度，本发明所提供的卫星遥感影像大气校正方法，在执行大气校正之前，基于6S辐射传输模型，建立大气校正参数查找表LUT（Look Up Table），实现基于6S辐射传输查找表的逐像元大气校正。为了提高大气校正速度，本发明利用消息传递编程模型MPI（Message PassingInterface）主从并行程序设计模式，提供数据块驱动与功能链驱动的并行大气校正模型，根据数据块驱动策略将输入影像进行分块，主进程负责向从进程分发分块数据，从进程负责根据大气校正功能链处理不同的影像块，处理完毕后再将结果返回给主进程。其包括以下过程：

步骤0：建立大气校正参数查找表。该步骤是进行大气校正的前提，为卫星影像大气校正提供了大气校正参数查找表。

步骤1：影像数据分块。主进程分配数据及参数，主进程根据总进程数和影像大小，计算得到每个从进程分配到的影像块在整幅影像中的位置，将影像块及配置参数传递给各从进程。

步骤2：数据输入。从进程分别读取大气校正参数，并进行参数配置。所要读取的参数包括波段增益偏移值、大气外部太阳辐照度、LUT表。所述波段增益偏移值、大气外部太阳辐照度可以从卫星数据官方得到。

步骤3：执行并行大气校正。从进程根据主进程传递过来的位置信息，读取相应数据块，根据“影像亮度值计算-表观反射率计算-地表反射率计算”大气校正功能链进行大气校正；具体如下：

根据影像块的像素值与每个波段的增益、偏移的关系，将影像的像素值转为亮度值；

根据影像亮度值与太阳天顶角、大气外部的太阳辐照度、地球太阳距离系数的关系，将亮度值转为大气层顶表观反射率；

结合大气校正参数查找表，将表观反射率转换为地表反射率。根据像元所处的波段、太阳天顶角、卫星天顶角、相对方位角，查找LUT表中对应的大气校正参数，利用气溶胶光学厚度对该数据进行线性插值得到像元处的大气校正参数，利用该大气校正参数将大气层顶表观反射率转为地表反射率。

步骤4：主进程接收并写入处理结果。从进程的各个模块分别输出转换后的上述数据，主进程接收所述转换后的数据并写入处理结果。

步骤5：合并处理结果。主进程不断接收从进程发送的结果影像块，并将其写入输出文件。

本发明根据数据块驱动策略将输入影像进行分块，并行计算时节省了内存资源，同时，并行计算时间与影像大小相关。衡量并行计算效率的指标包括如下三个，

串行计算时间复杂度，计算公式为：

T_S=O(M²)          （1）

并行计算时间复杂度，计算公式为：

T_P=O(M²/(P-1))    （2）

加速比，计算公式为：

$S_P=\frac{T_S}{T_P}---\left(3\right)$

效率，计算公式为：

$E_P=\frac{S_P}{P}---\left(4\right)$

其中，M为影像规模大小，P为进程个数，O表示数量级的符号，加速比与进程数成线性关系，效率用于衡量多个处理器计算时的使用情况，通常小于等于1，越接近1，说明P个处理器的使用效率越高。

所述步骤0的建立大气校正参数查找表包括以下具体过程。该准备步骤在主从进程进行大气校正之前完成。具体如下：

根据研究区域的地理、气象条件以及影像的获取时间，确定大气模式、气溶胶模式、气压、水汽含量、臭氧含量、目标海拔高度等输入参数；利用6S辐射传输模型，创建不同气溶胶光学厚度、太阳天顶角、卫星天顶角、太阳与卫星相对方位角等条件下的LUT表。所述LUT表由内部大气反射率ρ_atm、总透过率

球面反照率S三个参数组成。

所述6S辐射传输模型的输入参数主要类型包括：时间参数、几何参数、大气模式、气溶胶参数、海拔高度、光谱条件、地表参数等。其中气溶胶光学厚度、太阳天顶角、卫星天顶角、太阳与卫星相对方位角对大气校正的影响较大。由此输入参数设置为，气溶胶光学厚度从0-1以0.1步长变化，共11个值；太阳天顶角从0-90以10步长变化，共有10个值：卫星天顶角从0-90以10步长变化，共有10个值；太阳与卫星相对方位角从0-180以10步长变化，共19个值；根据各个传感器特征设置光谱参数，其他参数设为常量，本实施例中，传感器设有4个波段。将上述参数组合代入6S辐射传输模型循环计算，得到大气校正输入参数查找表，其流程如图2所示，包括如下步骤：

步骤21：设置气溶胶光学厚度取值个数i为0，所述气溶胶光学厚度从0开始取值，每次气溶胶光学厚度值增加0.1，且每增加一次，i值增加1。

步骤22：当i小于11时，执行下一步骤；否则，结束。

步骤23：设置太阳天顶角的取值个数j为0，所述太阳天顶角从0开始取值，每次太阳天顶角的角度值增加10，且每增加一次，j值增加1。

步骤24：当j小于10时，执行下一步骤；否则，结束。

步骤25：设置卫星天顶角的取值个数k为0，所述卫星天顶角从0开始取值，每次卫星天顶角的角度值增加10，且每增加一次，k值增加1。

步骤26：当k小于10时，执行下一步骤；否则，结束。

步骤27：设置相对方位角取值个数l为0，所述相对方位角从0开始取值，每次相对方位角的角度值增加10，且每增加一次，l值增加1。

步骤28：当l小于19时，执行下一步骤；否则，结束。

步骤29：设置传感器波段中心波长取值个数m为0，对所有传感器逐一取值，每取值一次，m增加1。

步骤210：m小于4时，执行下一个步骤；否则，结束。

步骤211：将上述所取的各参数值代入6s辐射传输模型进行计算。

依据上述程序，6S查找表输入参数设置如表1。

表1  6S查找表输入参数设置

上述6S查找表输入参数设置并不是固定的，可以根据研究区域的实际情况，细化输入参数，得到更精细的输入参数设置表。

步骤3中所述的大气校正执行过程，包括图3所示的以下步骤：

步骤3 1：计算输入影像亮度值。根据像输入影像的像素值与输入参数的每个波段的增益、偏移的关系，计算输入影像亮度值。

$L_{\mathrm{TOA}}^k=\frac{X^k}{{\mathrm{\α}}_k}+{\mathrm{\β}}_k---\left(5\right)$

其中：

为亮度（单位：Wm^-2sr^-1um^-1），X^k为输入影像的像素值，α_k为每个波段的增益，β_k为每个波段的偏移。

步骤32：计算大气层顶表观反射率。根据影像亮度值与太阳天顶角、大气外部的太阳辐照度、地球太阳距离系数的关系，计算大气层顶表观反射率。

${\mathrm{\ρ}}_{\mathrm{TOA}}^k=\frac{\mathrm{\π}{L}_{\mathrm{TOA}}^k}{E_S^k\·\cos \left({\mathrm{\θ}}_s\right).d/d_0}---\left(6\right)$

其中，

为大气层顶表观反射率，θ_s为太阳天顶角，为离地球d₀量测的大气外部的太阳辐照度，d/d₀为地球太阳距离系数，π取值为3.1415926，为步骤31所求得到的影像亮度值。

步骤33：计算地表反射率。根据像元所处的波段、太阳天顶角、卫星天顶角、相对方位角，查找LUT表中对应的大气校正参数，利用气溶胶光学厚度对该数据进行线性插值得到像元处的大气校正参数，利用该大气校正参数将大气层顶表观反射率转为地表反射率。

${\mathrm{\ρ}}_S^{\mathrm{unif}}=\frac{A}{1+\mathrm{SA}}---\left(7\right)$

$A=\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}-{\mathrm{\ρ}}_{\mathrm{atm}}}{T\left(\mathrm{\μS}\right).T\left(\mathrm{\μV}\right).t_g^{\mathrm{al\; lg\; as}}}---\left(8\right)$

其中，

为假设的朗伯面和均一环境下的地表反射率，即为大气校正数据结果，ρ_TOA为大气顶部的反射率，ρ_atm为内部大气反射率，

为大气球面反照率（总大气传输），T(μS)为下行散射，T(μV)为上行散射，S为球面反照率。

从上面所述可以看出，本实施例通过建立6S辐射传输大气校正参数查找表，解决时间和精度问题，通过建立块驱动与功能链驱动的并行大气校正模型，提高大气校正的计算速度，能够有效处理海量遥感影像数据，满足业务化大气校正的需求。

以环境星CCD（Charge Coupled Device，电荷耦合器件）数据大气校正为实施例，采用本发明所提供的卫星遥感影像大气校正模块进行大气校正。其中，CCD幅宽为360公里，各个像元的大气条件、太阳位置、卫星位置等均不相同，甚至有较大差异。本实施例详细描述环境星CCD数据的并行大气校正过程，试验数据为2009年6月25日太湖区域的HJ1A CCD1（环境1号A卫星搭载的CCD传感器1号）卫星影像。图4为本实施例的卫星遥感影像大气校正模块所执行的校正步骤，具体如下：

步骤41：建立大气校正参数查找表。

利用6S模型构建大气校正参数查找表。输入参数如下表2所示，

表2  环境星CCD 6S查找表输入参数设置

利用6S模型计算得到的大气校正参数查找表示例如表3，依据表2的数据取值个数，所建立的查找表由13×13×16×11×4=118976组大气校正参数（内部大气反射率、总透过率、球面反照率）构成。

表3  大气校正参数查找表示例

步骤42a，获取大气校正输入参数，包括像元角度信息、像元气溶胶光学厚度值、定标系数、波段光谱曲线、大气层外太阳辐照度信息。每个参数的获取方式如下，

像元气溶胶光学厚度值。实测该区域典型样点的气溶胶光学厚度值，或者从“AERONET”站点获取，或者通过反演的方法获取，通过插值方法得到每个像元的气溶胶光学厚度值。

定标系数。从HJ CCD（环境卫星上的CCD传感器）影像的元数据XML文件中获取HJ CCD相机定标系数，即各个波段的增益与偏移值。

波段光谱曲线。从中国资源卫星应用中心官方网站上获取环境星CCD光谱响应函数值，本实施例中所获取的曲线如图5所示。

大气层外太阳辐照度信息。从中国资源卫星应用中心官方网站上获取大气层外太阳辐照度信息。

传感器	波段1	波段2	波段3	波段4
					HJ-1-A CCD1	1914.324	1825.419	1542.664	1073.826
HJ-1-A CCD2	1929.810	1831.144	1549.824	1078.317
					HJ-1-B CCD1	1902.188	1833.626	1566.714	1077.085
HJ-1-B CCD2	1922.897	1823.985	1553.201	1074.544

表4  HJ1 A/B星CCD大气层外太阳辐照度（单位：w/m2）

表4  为从中国资源卫星应用中心官方网站上获取的大气外太阳辐照度信息。其中HJ-1-A CCD2为环境1号A卫星搭载的CCD传感器2号，HJ-1-BCCD1为环境1号B卫星搭载的CCD传感器1号，HJ-1-B CCD2为环境1号B卫星搭载的CCD传感器2号。所述波段1、波段2、波段3、波段4分别为蓝波段、绿波段、红波段、近红外波段。

步骤42b，配置并行计算环境。

基于遥感数据集群处理平台，配置执行大气校正的资源数，如运行队列、节点数、进程数、所需内存大小等。并行环境的配置如下表所示：

表5  遥感数据集群处理平台配置

步骤42c，配置并行大气校正参数。主进程配置大气校正参数查找表、像元气溶胶光学厚度值、像元角度信息、定标系数、大气层外太阳辐照度信息等参数，将配置参数传递给从进程。同时，主进程根据总进程数和影像大小，计算得到每个从进程分配到的影像块在整幅影像中的位置，将数据块传递给各从进程。

步骤43，执行并行大气校正。根据主进程传递过来的位置信息，读取相应的影像数据块，根据大气校正功能链进行大气校正，进程之间互不干预，完全处于独立状态。各从进程分别从本地读取大气校正参数，所读取的参数包括波段增益偏移值、辐照度、大气校正参数查找表，所述影像数据块为并行处理影像时，将影像分成的方块；所述大气校正功能链为大气校正流程的自动化表达。

根据波段增益偏移值进行影像亮度值计算；根据影像亮度值与太阳天顶角、大气外部的太阳辐照度、地球太阳距离系数的关系，进行大气顶部表观反射率计算；根据本数据块中像元所处的波段、太阳天顶角、卫星天顶角、相对方位角，查找LUT表中对应的大气校正参数，利用气溶胶光学厚度对内部大气反射率、总透过率、球面反照率这三个大气校正参数进行线性插值，得到各个像元的大气校正参数，利用该大气校正参数进行地表反射率计算。

步骤44，主进程接收并写入处理结果。从进程的各个模块分别输出转换后的上述数据，主进程接收所述转换后的数据并写入处理结果。

步骤45，合并并行执行结果，从进程将大气校正结果即地表反射率传递给主进程，主进程不断接收从进程发送的结果影像块，并将其写入输出影像。

为了验证本实施例的卫星遥感影像大气校正模块执行大气校正的精度，本发明与均一输入参数的6S模型大气校正进行比较。

参数	值
		获取时间	2009年6月25日
大气模式	中纬度夏季
		传感器高度	-1000
太阳天顶角	26.45

太阳方位角	301.3
		卫星天顶角	32.5
卫星方位角	119.9
		气溶胶模式	大陆型
气溶胶光学厚度	0.2

表6  HJ-1A/B星CCD均一输入参数设置

表6为HJ-1A/B星CCD均一输入参数设置。将本实施例所提供的卫星遥感影像大气校正模块所获得的大气校正结果与均一参数大气校正模块所获得的结果进行比较，植被样点平均反射率如表7所示。

波段	校正前反射率	均一参数反射率	本发明反射率
				1	0.1276	0.0771	0.0729
2	0.1123	0.0657	0.0604
				3	0.0891	0.0623	0.0539
4	0.1989	0.2162	0.2044

表7  不同校正方法植被样点平均反射率

表7中所述的波段1、2、3、4分别为环境卫星CCD影像的蓝波段、绿波段、红波段、近红外波段。表7中所反映的植被样点平均反射率的对比如图6所示。从表7和图6中可以看出，大气校正后可见光波段的平均反射率明显比校正前小，蓝波段平均反射率下降明显；近红外波段的平均反射率比校正前稍稍增大，但变化程度较小。

进一步，利用NDVI值（归一化差分植被指数，Normalized DifferenceVegetation Index）验证大气校正的效果，在所述影像块中随机选取45个像素点进行分析，所选取的45个点覆盖不同地表覆盖类型，分析结果如图7所示，大气校正使NDVI值有所增加；因为大气校正后红光波段反射率降低，而近红外波段反射率增加，大气传播过程中受气溶胶、水气、空气分子等因素的影响，使彼此之间的对比度降低，造成了NDVI信号的衰减，而本发明所提供的大气校正的方法正好能够弥补了这种衰减。从上面所述可以看出，本发明所提供的卫星遥感影像大气校正后NDVI增加幅度大于均一参数大气校正后的NDVI，证明了本发明所提供的大气校正方法的有效性和具有较高的精确性。

为了验证本发明所提供的卫星遥感影像大气校正速度，将其校正速度数据与串行大气校正进行比较。CCD影像（15811×13974）并行大气校正执行时间、加速比和效率见表8。

进程数	执行时间（秒）	加速比	效率
				1	1223	1	1
2	1165	1.05	0.52
				3	608	2.01	0.67
4	388	3.15	0.79
				5	306	4	0.8
6	251	4.87	0.81
				7	208	5.88	0.84
8	181	6.77	0.85

表8  CCD影像并行大气校正执行时间、加速比和效率

图8为表8中所述的大气校正并行执行时间随着进程数的变化图。从所述的并行执行时间可以看出，进程数从1增加到2时，运算时间有小幅度的减小，2个进程时只有1个从进程参与运算。进程数小于7时，运算时间随着进程数的增加而大幅减少，大于7时减少缓慢，原因是随着进程数的增加，计算时间不断减少，通信时间和文件读写时间比重逐渐上升。

图9为表8中所述的大气校正并行执行加速比与进程数的线性关系图。从所述加速比可以看出，并行计算的加速比与进程数成线性关系，随着处理进程数的增加，并行算法中由于通信时间占计算时间的比例上升，导致处理进程空闲和通信的增多，因而加速比略有减小。

通过对比分析大气校正精度与速度，证明本发明所提供的卫星遥感影像大气校正模块的大气校正精度优于均一参数大气校正模块，并行大气校正速度明显优于串行大气校正。因此，本发明是一种适合于卫星遥感影像业务化运行的并行大气校正方法。

进一步，本发明提出一种卫星遥感影像大气校正模块，采用上面所述的卫星遥感影像大气校正方法，进行大气校正；包括至少一个子模块，所述子模块执行卫星遥感影像的并行大气校正。

所述子模块包括查找模块、计算模块、输入模块和输出模块。

所述查找模块根据不同的气象条件参数，利用6S辐射传输模型，创建由大气校正参数组成的大气校正参数查找表。

所述输入模块用于输入大气校正所需参数。

所述计算模块根据各个波段的增益与偏移值将影像的像素值转为亮度值；计算影像亮度值时使用如下公式：

$L_{\mathrm{TOA}}^k=\frac{X^k}{{\mathrm{\α}}_k}+{\mathrm{\β}}_k---\left(5\right)$

其中：

为亮度（单位：Wm^-2sr^-1um^-1），X^k为输入影像的像素值，α_k为每个波段的增益，β_k为每个波段的偏移。

所述计算模块根据太阳天顶角、大气外部的太阳辐照度、地球太阳距离系数将亮度值转为大气层顶表观反射率；大气层顶表观反射率计算时使用如下公式：

${\mathrm{\ρ}}_{\mathrm{TOA}}^k=\frac{\mathrm{\π}{L}_{\mathrm{TOA}}^k}{E_S^k\·\cos \left({\mathrm{\θ}}_s\right).d/d_0}---\left(6\right)$

其中，

为大气层顶表观反射率，θ_s为太阳天顶角，

为离地球d₀量测的大气外部的太阳辐照度，d/d₀为地球太阳距离系数。

根据像元所处的波段、太阳天顶角、卫星天顶角、相对方位角，查找大气校正参数查找表表中对应的大气校正参数，利用气溶胶光学厚度对该数据进行插值计算得到像元处的大气校正参数，根据此大气校正参数将大气层顶表观反射率转为地表反射率，地表反射率计算时依据如下公式：

${\mathrm{\ρ}}_S^{\mathrm{unif}}=\frac{A}{1+\mathrm{SA}}---\left(7\right)$

$A=\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}-{\mathrm{\ρ}}_{\mathrm{atm}}}{T\left(\mathrm{\μS}\right).T\left(\mathrm{\μV}\right).t_g^{\mathrm{al\; lg\; as}}}---\left(8\right)$

其中，

为假设的朗伯面和均一环境下的地表反射率，ρ_TOA为大气顶部的反射率，ρ_atm为内部大气反射率，

为大气球面反照率（总大气传输），T(μS)为下行散射，T(μV)为上行散射，S为球面反照率。

所述输出模块输出大气校正结果。

所属领域的普通技术人员应当理解：以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种卫星遥感影像大气校正方法，其特征在于，至少包括以下步骤：

根据不同的气象条件参数，利用6S辐射传输模型，创建由大气校正参数组成的大气校正参数查找表；

根据各个波段的增益与偏移值将影像的像素值转为亮度值；

根据太阳天顶角、大气外部的太阳辐照度、地球太阳距离系数将所述亮度值转为大气层顶表观反射率；

根据像元所处的波段、太阳天顶角、卫星天顶角、相对方位角，查找大气校正参数查找表中对应的大气校正参数，利用气溶胶光学厚度对查找到的大气校正参数进行插值计算得到像元处的大气校正参数，根据该像元处的大气校正参数将大气顶层表观反射率转为地表反射率；

所述影像的亮度值的计算使用如下公式：

$L_{\mathrm{TOA}}^k=\frac{X^k}{{\mathrm{\α}}^k}+{\mathrm{\β}}_k---\left(5\right)$

其中，

为亮度，X^k为输入影像的像素值，α_k为每个波段的增益，β_k为每个波段的偏移；

所述大气层顶表观反射率的计算使用如下公式：

${\mathrm{\ρ}}_{\mathrm{TOA}}^k=\frac{\mathrm{\π}{L}_{\mathrm{TOA}}^k}{E_S^k.\cos \left({\mathrm{\θ}}_s\right).d/d_0}---\left(6\right)$

其中，

为大气层顶表观反射率，θ_s为太阳天顶角，

为离地球d₀量测的大气外部的太阳辐照度，d/d₀为地球太阳距离系数；

所述地表反射率的计算依据如下公式：

${\mathrm{\ρ}}_S^{\mathrm{unif}}=\frac{A}{1+\mathrm{SA}}---\left(7\right)$

$A=\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}-{\mathrm{\ρ}}_{\mathrm{atm}}}{T\left(\mathrm{\μS}\right).T\left(\mathrm{\μV}\right).t_g^{\mathrm{allgas}}}---\left(8\right)$

其中，

为假设的朗伯面和均一环境下的地表反射率，ρ_TOA为大气顶部的反射率，ρ_atm为内部大气反射率，

为大气球面反照率，T(μS)为下行散射，T(μV)为上行散射，S为球面反照率。

2.根据权利要求1所述的卫星遥感影像大气校正方法，其特征在于，所述大气校正参数查找表的建立包括如下步骤：

根据研究区域的地理气象条件以及影像的获取时间，确定大气环境参数；

利用6S辐射传输模型，创建不同大气校正输入参数条件下的大气校正参数查找表；

所述大气校正参数查找表输出的大气校正参数包括内部大气反射率、总透过率、球面反照率。

3.根据权利要求2所述的卫星遥感影像大气校正方法，其特征在于，所述大气环境参数包括大气模式、气溶胶模式、气压、水汽含量、臭氧含量、目标海拔高度；所述大气校正输入参数包括气溶胶光学厚度、太阳天顶角、卫星天顶角、太阳与卫星相对方位角。

4.根据权利要求1所述的卫星遥感影像大气校正方法，其特征在于，所述卫星遥感影像大气校正方法采用数据块驱动与功能链驱动的并行大气校正模型，且基于消息传递编程模型的主从并行程序设计模式来执行大气校正。

5.根据权利要求4所述的卫星遥感影像大气校正方法，其特征在于，所述卫星遥感影像大气校正方法采用主从并行计算的模式，包括：

主进程分配数据及参数，所述主进程根据总进程数和影像大小，计算得到每个从进程分配到的影像块在整幅影像中的位置，将所述影像块及主进程所分配的参数传递给所述从进程；

所述从进程读取大气校正参数，并进行参数配置，所读取的大气参数为：增益偏移值、辐照度、大气校正参数查找表；

所述从进程根据主进程传递过来的所述影像块在整幅影像中的位置信息，读取相应数据块，进行大气校正；

主进程不断接收从进程发送的结果影像块，并将其写入输出文件。

6.一种卫星遥感影像大气校正模块，其特征在于，包括如下子模块：

查找模块、计算模块、输入模块和输出模块；

所述查找模块根据不同的气象条件参数，利用6S辐射传输模型，创建由大气校正参数组成的大气校正参数查找表；

所述输入模块用于输入大气校正所需参数；

所述计算模块根据各个波段的增益与偏移值将影像的像素值转为亮度值；根据太阳天顶角、大气外部的太阳辐照度、地球太阳距离系数将亮度值转为大气顶层表观反射率；根据像元所处的波段、太阳天顶角、卫星天顶角、相对方位角，查找大气校正参数查找表表中对应的大气校正参数，利用气溶胶光学厚度对查找到的大气校正参数进行插值计算得到像元处的大气校正参数，根据此大气校正参数将大气顶层表观反射率转为地表反射率；

所述输出模块输出大气校正结果；

所述计算模块计算影像亮度值时使用如下公式：

$L_{\mathrm{TOA}}^k=\frac{X^k}{{\mathrm{\α}}^k}+{\mathrm{\β}}_k---\left(5\right)$

其中，

为亮度，X^k为输入影像的像素值，α_k为每个波段的增益，β_k为每个波段的偏移；

所述计算模块计算大气层顶表观反射率时使用如下公式：

${\mathrm{\ρ}}_{\mathrm{TOA}}^k=\frac{\mathrm{\π}{L}_{\mathrm{TOA}}^k}{E_S^k.\cos \left({\mathrm{\θ}}_s\right).d/d_0}---\left(6\right)$

其中，

为大气层顶表观反射率，θ_s为太阳天顶角，

为离地球d₀量测的大气外部的太阳辐照度，d/d₀为地球太阳距离系数；

所述计算模块计算地表反射率时依据如下公式：

${\mathrm{\ρ}}_S^{\mathrm{unif}}=\frac{A}{1+\mathrm{SA}}---\left(7\right)$

$A=\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}-{\mathrm{\ρ}}_{\mathrm{atm}}}{T\left(\mathrm{\μS}\right).T\left(\mathrm{\μV}\right).t_g^{\mathrm{allgas}}}---\left(8\right)$

其中，

为假设的朗伯面和均一环境下的地表反射率，ρ_TOA为大气顶部的反射率，ρ_atm为内部大气反射率，为大气球面反照率，T(μS)为下行散射，T(μV)为上行散射，S为球面反照率。

Images