CN108256186B - 一种在线计算查找表的逐像元大气校正方法 - Google Patents

Abstract

本发明公布了一种基于6S大气辐射传输软件包在线计算的遥感影像逐像元大气校正方法，本发明方法利用6S在线并行计算查找表，首先生成遥感影像在不同成像几何与大气模式下不同气溶胶光学厚度的查找表，基于该查找表获取气溶胶分布数据，然后生成大气校正系数的查找表，基于该查找表进行逐像元大气校正，该发明改善了传统方法离线建立查找表存在的适用范围有限、建表量大以及逐像元多维插值运算复杂的问题，对于高分四号地球同步轨道卫星在全天任意时刻获取的高分辨率遥感影像大气校正都具有较高的稳定性与适用性。

Description

一种在线计算查找表的逐像元大气校正方法

技术领域

本发明涉及一种遥感影像的大气校正技术，具体的说，涉及一种针对地球同步轨道高分辨率遥感卫星图像的在线计算查找表的逐像元大气校正技术。

背景技术

遥感图像的辐射预处理一直是遥感数据处理的主要课题之一。新的遥感卫星投入使用后，对新数据的预处理便是紧要的问题。遥感图像的预处理通常包含几何预处理于辐射预处理，其中的辐射预处理主要的目的是从遥感传感器观测到的辐射信号中去除大气的吸收、大气分子的瑞利散射和气溶胶的米散射影响，还原地表的真实反射率。大气校正是遥感图像辐射处理的关键步骤，目前常用的大气校正方法是基于辐射传输模型的方法，该类方法通过大气状况、卫星和太阳位置信息等输入参数，利用大气辐射传输模型计算大气透过率、程辐射等大气校正参数，以此为依据反演气溶胶光学厚度与计算地表反射率。目前已经产生了很多可选择的辐射传输模型算法软件包，不同作者的算法原理基本相同，差异主要体现在适用的假设条件和范围，常用的有6S(Second Simulation of the SatelliteSignal in the Solar Spectrum)、LOWTRAN(Low Resolution Transmission)、MORTRAN(Moderate Resolution Transmission)等。利用这些辐射传输模型算法软件包，开发了很多大气校正软件，比如ACORN(Atmospheric CORrection Now)、ATCOR(Aspatially-Adaptive Fast Atmospheric and Topographic Correction)和FLAASH(Fast Line ofsight Atmospheric Analysis of Spectral Hypercubes)。为了获取更高精度的大气校正结果，针对特定遥感卫星数据还出现了专门开发的算法，比如MODIS数据C5版本的地表反射率反演算法、Landsat TM/ETM+数据的LEDAPS算法和Landsat8 OLI数据的LASRC算法等。中国近年来连续发射了HJ系列、ZY系列和GF系列多颗卫星，针对搭载的CCD传感器获取的高分辨率多光谱数据，国内研究者也开发了多套针对国产卫星数据的大气校正算法软件用于地表反射率数据产品的生产。

高分四号卫星(简称GF-4)是中国于2015年12月发射的一颗地球同步轨道卫星，轨道高度35786千米，搭载空间分辨率为50米的全色、多光谱相机和400米分辨率的中波红外相机，采用面阵凝视方式成像，成像间隔快至20秒，具备高时间、高空间分辨率的优势，可在任意时刻对局部区域实现机动观测，在防灾减灾与应急响应、环境资源监测调查等方面具有重要意义。作为一颗地球静止卫星，GF-4卫星多光谱数据却具有50米的高分辨率，成像覆盖范围相对于以往几千米分辨率的静止卫星数据小很多，幅宽只有400千米，与传统的太阳同步轨道卫星获取的多光谱数据有很多相似的特性。针对GF-4数据生产大气校正产品是目前迫切的工程需求。

GF-4卫星成像的特点使得其多光谱数据的气溶胶反演与大气校正难以直接使用传统的基于离线查找表的方法框架。由于采用地球同步轨道，卫星可以在一天中的任意时间成像，使得卫星-地表-太阳三者之间的成像几何特性复杂，使得建立离线查找表需要考虑的角度组合数成倍增加，为了保证精度就需要建立更加复杂的离线查找表，同时逐像元多维插值的运算量也将大幅增加。实验表明复杂的查找表插值难以将单幅GF-4图像的处理时间同样控制在30分钟以内。在气溶胶估计环节，由于地球同步轨道成像时间的任意性与太阳同步轨道不同，无法利用相近时刻成像的其它卫星数据的地表反射率产品或气溶胶产品，并且由于缺少短波红外波段，通过图像数据估计气溶胶厚度更加困难。传统的基于查找表的大气校正方法需要进行结构上的调整，克服离线查找表存在的建表复杂与插值慢的问题，改善基于图像数据反演气溶胶的方法，针对高分辨率地球同步轨道卫星数据成像特性定制逐像元大气校正的专门算法流程框架，提供一种稳健的GF-4卫星逐像元大气校正算法程序，具有自动获取与分析大气校正所需的输入参数能力，具有批处理能力，以便满足工程化生产高精度GF-4大气校正数据产品，即气溶胶与地表反射率数据。

发明内容

本发明的目的是针对静止卫星高分辨率遥感图像辐射预处理应用，提供一种逐像元大气校正技术，特别是对于GF-4卫星50米空间分辨率的多光谱图像的L1级数据产品中的气溶胶光学厚度与地表反射率数据产品生产，提供一种自动稳健的大气校正算法流程。本技术基于成熟的6S大气辐射传输软件包及传统的基于查找表的气溶胶反演与大气校正方法，根据GF-4卫星图像的辐射预处理需求与图像成像特性定制的、在线生成查找表的快速大气校正算法流程。

本发明的基本思路为：考虑到GF-4卫星多光谱数据幅宽只有400千米，覆盖范围小，可以假设一幅图像内的太阳和传感器的观测角度一致，且具有相同的海拔高度，只有气溶胶光学厚度是不同的，于是仅对不同的气溶胶光学厚度建立气溶胶反演查找表与大气校正参数查找表，生成这两个查找表需要调用6S的次数不到一百次，可在线调用6S生成，同时由于查找表仅对气溶胶厚度插值，逐像元运算的速度大幅提高，单幅处理时间可提高到10分钟以内。基本处理流程是：首先是自动获取大气校正参数，通过分析图像辅助数据，获取大气校正所需成像几何角度信息与成像时间，并结合低分辨率全球DEM高程数据与全球地表分类数据选择对应大气模式与气溶胶模式等参数。然后针对不同的气溶胶光学厚度，在线调用6S生成针对当前图像的气溶胶光学厚度查找表与大气校正参数查找表，再利用气溶胶查找表结合图像波段特性反演气溶胶数据，最后利用大气校正参数查找表结合气溶胶数据完成逐像元大气校正，获得地表反射率数据。

所述的静止卫星高分辨率遥感图像，限定为地球同步轨道卫星获取的幅宽小于400千米的多光谱图像，或者单幅图像中幅宽小于400千米的局部区域，图像辅助数据包含大气校正所需的成像参数与辐射定标系数。

本发明技术方案提供的一种在线计算查找表的逐像元大气校正方法，其特征在于包括以下实施步骤：A数据准备，辐射定标与计算表观反射率，通过参数解析自动获得大气模式、气溶胶模式与平均高程，在6S软件包中添加卫星多光谱数据各波段光谱响应数据；

B在线生成查找表，利用自动获取的太阳天顶角、观测天顶角、相对方位角、成像时间、大气模式、气溶胶模式与平均高程，针对不同的气溶胶光学厚度，在线调用6S生成针对当前图像的气溶胶光学厚度查找表与大气校正参数查找表；

C气溶胶光学厚度反演，利用气溶胶查找表，结合图像波段特性与修正的DDV(darkdense vegetation)法反演暗植被像元处的气溶胶光学厚度，剔除非暗植被点，再通过克里金插值与大尺度中值滤波平滑获取整幅图像的气溶胶分布数据；

D逐像元大气校正，利用大气校正参数查找表与气溶胶分布数据，逐像元查找表插值出大气校正参数，并代入大气校正公式获得地表反射率数据。

上述实施步骤的特征在于：

步骤A中所述数据准备，包括辐射定标与计算表观反射率，通过参数解析自动获得大气模式、气溶胶模式与平均高程，在6S软件包中添加GF-4卫星多光谱数据各波段光谱响应数据。所述辐射定标与计算表观反射率，目的是将遥感图像原始的DN(Digtial Number)值转换到表观反射率；所述通过参数解析自动获得大气模式、气溶胶模式与平均高程，是通过卫星数据XML辅助文件中获取的成像时间与地理范围，结合低分辨率的全球高程数据确定平均高程，结合低分辨率的全球地表分类数据确定大气模式与气溶胶模式，其中全球地表分类数据分类类别包含裸地、水体、城市、荒漠、高原积雪5类；所述在6S软件包中添加GF-4卫星多光谱数据各波段光谱响应数据，采用的6S版本为6SV2.1，通过修改6S源码的方式添加GF-4卫星多光谱数据各波段光谱响应数据并重新编译，使得6S软件包原生支持GF-4卫星数据，避免每次在线调用重新定义光谱条件；

步骤B中所述的在线生成查找表，包括利用自动获取的太阳天顶角、观测天顶角、相对方位角、成像时间、大气模式、气溶胶模式与平均高程，针对不同的气溶胶光学厚度，在线调用6S生成针对当前图像的气溶胶光学厚度查找表与大气校正参数查找表；所述针对不同的气溶胶光学厚度，是在常用气溶胶光学厚度取值范围[0,3.5]内，根据查找表建表插值精度确定气溶胶光学厚度的个数，气溶胶光学厚度查找表取值为(0,0.05,0.1,0.2,0.5,1,1.5,2)共8个值，大气校正参数查找表取值为(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,2,2.5,3,3.5)共20个值；所述气溶胶光学厚度查找表，内容为对于参与气溶胶估计的波段，每个气溶胶光学厚度值对应三个由在线调用6S获取的变量值T、ρ₀、S，采用的辐射传输方程为

其中ρ_s为地表反射率、ρ_TOA为表观反射率、T为大气透过率、ρ₀为大气的路径辐射项等效反射率、S为大气下界的半球反射率；所述大气校正参数查找表，内容为对于多光谱图像的所有波段，每个气溶胶光学厚度值对应三个由在线调用6S获取的变量值a、b、c，采用的大气校正公式为y＝aρ_TOA-b，

步骤C中所述序气溶胶光学厚度反演，包括利用气溶胶查找表，结合图像波段特性与修正的DDV法反演暗植被像元处的气溶胶光学厚度，剔除非暗植被点，再通过克里金插值与大尺度中值滤波平滑获取整幅图像的气溶胶分布数据；所述结合图像波段特性与修正的DDV法反演暗植被像元处的气溶胶光学厚度，是利用在暗植被处像元红波段地表反射率

与蓝波段地表反射率

存在固定的线性关系

红蓝地表反射率比率k取值为2，计算方法是利用表观反射率计算归一化植被指数

将NDVI值大于0.6的像元标记为暗植被，对于每一个暗植被像元，将红蓝波段气溶胶光学厚度查找表8个不同的气溶胶光学厚度值所对应的T、ρ₀、S参数带入辐射传输方程，由最满足该线性关系的两个气溶胶光学和厚度之间插值出该暗植被像元处的气溶胶光学厚度值；所谓剔除非暗植被点，是考虑到暗植被的选取是在有大气的影响下确定的，当气溶胶光学厚度过大时可能会出现误判，剔除误判的方法是将所有反演得到的气溶胶光学厚度代入红波段与近红外波段的大气校正参数查找表，插值出大气校正参数，代入大气校正公式得到地表反射率，根据地表反射率重新计算NDVI，将大气校正后NDVI小于0.6的像元反演的气溶胶光学厚度值剔除；所述通过克里金插值与大尺度中值滤波平滑获取整幅图像的气溶胶分布数据，是将只在暗植被处反演得到的气溶胶光学厚度值插值到整个图像的过程，克里金插值在固定插值半径的情况下被执行多次直到所有像素位置都获得气溶胶光学厚度，最后利用大尺度中值滤波进行平滑处理；所述大尺度中值滤波指滤波半径大于10个像素的中值滤波，目的是获得边缘过度更加平滑的、更符合实际气溶胶空间分布特性的气溶胶光学厚度图像；

步骤D中所述逐像元大气校正，包括利用大气校正参数查找表与气溶胶分布数据，逐像元查找表插值出大气校正参数，并代入大气校正公式获得地表反射率数据；所述逐像元查找表插值出大气校正参数，是对于GF-4每一个波段的表观反射率数据，逐像素在气溶胶光学厚度数据中找到对应的气溶胶光学厚度值，根据气溶胶光学厚度值在该波段的大气校正参数查找表中从该气溶胶光学厚度值所在的前后两个气溶胶光学厚度值对应的大气校正参数中线性插值出3个大气校正参数，与表观反射率一同代入大气校正公式计算得到该像元的地表反射率。

本发明与现有技术相比有如下特点：本发明提供了一种针对地球静止卫星高分辨率多光谱图像的气溶胶与地表反射率反演解决方案，相比较传统大气校正的离线建立查找表、在线多维复杂插值的方法，本发明根据GF-4卫星幅宽小的特点，利用在线生成针对当前图像的查找表的方式，避免了离线建立静止卫星多时刻成像导致的复杂观测几何所需要的庞大查找表，并且将多维查找表插值简化为只针对不同气溶胶光学厚度的单维线性插值，将一般单幅图像大气校正的时间由30分钟提高到5分钟以内。算法自动化程度高，利用全球地表分类数据与高程数据自动判断当前校正图像所对应的大气模式、气溶胶模式与平均高程，整个处理流程无需人机交互，工程上可实现批处理。涉及的关键步骤采用成熟的算法实现，具有较高的稳定性与适用性。为GF-4卫星数据预处理中的气溶胶与地表反射率数据产品的生产提供了关键的技术支撑。

附图说明：

图1是在线计算查找表的逐像元大气校正方法流程图

图2是大气校正参数估计用的高程数据示意图

图3是大气校正参数估计用的全球分类数据示意图

图4是GF-4卫星蓝光、绿光、红光、近红外四个波段的光谱响应函数示意图

图5是气溶胶光学厚度数据克里金插值与大尺度中值滤波平滑处理示意图

具体实施方式：

本发明实现一种在线计算查找表的逐像元大气校正方法流程如图1所示，现结合附图对其进行描述。

图1中“输入”部分的“图像DN值数据”针对地球同步轨道卫星高分辨率多光谱数据，限制图像幅宽不超过400千米，图像波段至少包含蓝波段、绿波段、红波段、近红外波段4个。本发明实例数据是GF-4卫星多光谱图像标准数据产品。

图1中“输入”部分的“XML文件”是图像标准数据产品包含的、用于保存图像成像基本信息的文件，目前卫星数据的该类文件通常保存为XML格式，通过该文件通常可直接获取、或者通过简单的数据转换获取大气校正所需要的参数，包括用于辐射定标的定标系数、太阳天顶角、观测天顶角、图像获取时间与覆盖区域地理坐标范围等，相对方位角可通过太阳方位角与观测方位角计算得到。

图1中“输入”部分的“全球高程数据”是一幅单波段低分辨率全球遥感图像，像素值为地表高度。本发明实例采用的全球高程数据采用WGS-84经纬度投影，像素分辨率为0.0083度，图像尺寸为43200×20880像素，见示意图2。图中亮度越大表示地表高度越高，另外根据大气校正的需要，将高程值为负值的像素、以及海洋区域的像素都赋为0值。

图1中“输入”部分的“全球地表类别”是一幅单波段低分辨率全球遥感分类图像，分类类别包含裸地、水体、城市、荒漠、高原积雪5类，像素值为地表分类类别。本发明实例采用的全球地表分类数据采用WGS-84经纬度投影，像素分辨率为0.0083度，图像尺寸为43200×21600像素，像素值为1为裸地、2为水体、3为城市、5为荒漠、254为高原积雪、255为未分类别，其中1至5的编号与6S软件包中的气溶胶类型编号是对应的。全球地表类别见示意图3。图中左边为全球范围示例，右边为北京与天津局部区域放大示例。

图1中“输入”部分的“各波段光谱响应数据”是传感器在每个波长处接收的辐射亮度与入射的辐射亮度的比值，每个传感器的光谱响应数据都不同，本发明实例使用的GF-4卫星蓝光、绿光、红光、近红外四个波段的光谱响应函数见示意图4，可以看出前三个波段的波段宽度约为100纳米，近红外波段宽度约为150纳米。该数据按照6S软件包的光谱分辨率重采样后，通过扩展编写6S的Fortran源代码并重新编译获得6S对GF-4卫星数据的光谱支持。

图1中“数据准备”部分的“辐射定标与表观反射率计算”在遥感数据辐射处理中是相对固定的。辐射定标的目的是消除传感器本身产生的误差，利用一个线性公式将数据的原始DN值转化到表观辐射亮度L(单位：W·m^-2·sr^-1·μm^-1)。本发明实例使用的GF-4卫星多光谱数据的定标系数记录在XML文件中，该定标系数是经过外场同步定标，由中国资源卫星应用中心2016年10月中旬对外发布星上外场辐射定标系数。由于GF-4需根据地表亮度调整积分时间，各个波段的积分时间不一致。业务化运行采用了五种组合方式。下表给出了5种状态的定标系数，例如状态6-4-6-6，指的是蓝、绿、红、近红波段的积分时间分别是6、4、6和6ms。根据定标系数，可利用公式L＝DN/gain实现观测数据的辐射定标，获得各个波段的辐亮度。其中gain为表中的数据，GF-4定标数据的偏移项offset为0。

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

其中π为常量，D为日地之间距离(天文单位)，该距离和一年里的天数之间存在固定关系，可通过程序内建天数和D之间的查找表，根据天数确定。ESUN为大气层顶的平均太阳光谱辐照度(单位：W·m^-2·μm^-1)，该值与波段有关，本发明实例使用的GF-4多光谱各波段数据分别为蓝波段1967.207，绿波段1858.935，红波段1575.866，近红外波段1074.335。θ为太阳天顶角，从XML文件中获取。

考虑到单幅GF-4图像数据尺寸大，难以整个读入内存，辐射定标与表观反射率计算过程采用图像分块处理的方式进行，计算得到的表观反射率数据保存为float数据类型的tif图像磁盘文件，用户可选择作为中间结果保存，或者在大气校正结束后删除。

图1中“数据准备”部分的“参数自动解析”。目的是获取运行6S生成查找表的输入参数，其中数据获取时间、数据覆盖地理范围、太阳天顶角、观测天顶角、相对方位角这些参数通过解析XML文件并结合简单计算都能获取，平均高程获取需要根据数据覆盖地理范围，读取全球高程对应范围的高程值，计算排除前后2％极大值与极小值之后的、剩余像素的高程平均值。气溶胶模式获取需要根据数据覆盖地理范围，读取全球地表类别对应范围的各种类别对应的像素数量，计算像素数量最多的地表类别作为判断气溶胶模式的依据。地表类别与气溶胶模式的对应关系是裸地对应大陆气溶胶模式、水体对应海洋气溶胶模式、城市对应城市气溶胶模式、荒漠对应荒漠气溶胶模式、高原积雪对应无气溶胶。大气的水汽含量与气温含量很难获取，卫星数据辅助文件中也不会有记载，所以通过固定的大气模式替代输入6S，大气模式包含热带模式T(Tropical)、中纬度夏季模式MLS(Mid-LatitudeSummer)、中纬度冬季模式MLW(Mid-Latitude Winter)、亚极地夏季模式SAS(Sub-AricticSummer)、亚极地冬季模式SAW(Sub-Arictic Winter)五种，算法通过数据获取的纬度与时间判断所采用的大气模式，下表是大气模式与纬度、时间的关系表：

图1中“在线生成查找表”部分的“在线生成气溶胶查找表”。利用自动获取的太阳天顶角、观测天顶角、相对方位角、成像时间、大气模式、气溶胶模式与平均高程，针对不同的气溶胶光学厚度，在线调用6S生成针对当前图像的气溶胶光学厚度查找表与大气校正参数查找表；所述针对不同的气溶胶光学厚度，是在常用气溶胶光学厚度取值范围[0,3.5]内，根据查找表建表插值精度确定气溶胶光学厚度的个数，气溶胶光学厚度查找表取值为(0,0.05,0.1,0.2,0.5,1,1.5,2)共8个值。气溶胶光学厚度查找表内容为对于参与气溶胶估计的波段，每个气溶胶光学厚度值对应三个由在线调用6S获取的变量值T、ρ₀、S，采用的辐射传输方程为

其中ρ_s为地表反射率，ρ_TOA为表观反射率，T为大气透过率，对应6S输出结果total sca.trans.项的total值，ρ₀为大气的路径辐射项等效反射率，对应6S输出结果reflectance I项的total值，S为大气下界的半球反射率，对应6S输出结果spherical albedo项的total值。气溶胶光学厚度查找表的个数是参与气溶胶反演的卫星图像波段数量，GF-4卫星利用红光与蓝光2个波段进行气溶胶反演，则需要对于这两个波段建立2个气溶胶光学厚度查找表。

图1中“在线生成查找表”部分的“在线生成大气校正查找表”。利用自动获取的太阳天顶角、观测天顶角、相对方位角、成像时间、大气模式、气溶胶模式与平均高程，针对不同的气溶胶光学厚度，在线调用6S生成针对当前图像的气溶胶光学厚度查找表与大气校正参数查找表；所述针对不同的气溶胶光学厚度，是在常用气溶胶光学厚度取值范围[0,3.5]内，根据查找表建表插值精度确定气溶胶光学厚度的个数，大气校正参数查找表取值为(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,2,2.5,3,3.5)共20个值。大气校正参数查找表内容为对于多光谱图像的所有波段，每个气溶胶光学厚度值对应三个由在线调用6S获取的变量值a、b、c，对应6S输出结果atmospheric correction result项的coefficients xap xb xc值。采用的大气校正公式为y＝aρ_TOA-b，

气溶胶光学厚度查找表的个数是参与地表反射率反演的卫星图像波段数量，GF-4卫星为蓝光、绿光、红光、近红外全部波段建立4个大气校正查找表。

图1中“在线生成查找表”部分调用外部6S的程序次数有限，以GF-4为例，在线生成气溶胶查找表需要调用6S程序的次数为8个气溶胶厚度值与2个波段的乘积，共16次，在线生成大气校正查找表需要调用6S程序的次数为20个气溶胶厚度值与4个波段的乘积，共80次，总计耗时在2分钟左右，编程时采用简单的CPU多核加速技术将耗时提高到30秒以内，整个线生成查找表的算法耗时不高。图1中“在线生成查找表”部分生成的“气溶胶查找表”与“大气校正查找表”直接保存在内存的结构数组中。

图1中“气溶胶反演”部分的“气溶胶反演”。考虑到气溶胶的空间分布特性，气溶胶反演时对参与运算的波段数据进行降采样以减小运算量是常用的做法，对于GF-4数据气溶胶反演涉及的蓝光、红光与近红外3个波段的表观反射率数据，先降采样4倍，即200米分辨率，然后将采样后的波段数据全部读入内存，计算NDVI，将NDVI值大于0.6的像元标记为暗植被。对于暗植被像元，将气溶胶查找表中8个不同气溶胶光学厚度值所对应的T、ρ₀、S参数值代入辐射传输方程：

计算得到8组暗植被点在蓝光、红光波段的地表反射率，为了线性插值出满足

关系的气溶胶厚度值，计算8组地表反射率的

寻找c_i值出现正负号过渡的位置，比如假若出现c₁>0与c₂<0，c₁对应的气溶胶光学厚度为x₁，c₂对应的气溶胶光学厚度为x₂，则说明气溶胶光学厚度满足

关系的值介于x₁与x₂之间，使用简单的线性插值即可，公式为：

在所有暗植被点都插值出气溶胶光学厚度后，进行非暗植被点的剔除。剔除的方法是将所有反演得到的气溶胶光学厚度代入红波段与近红外波段的大气校正参数查找表，插值出大气校正参数，代入大气校正公式得到地表反射率，根据地表反射率计算NDVI，将大气校正后NDVI小于0.6的像元反演的气溶胶光学厚度值剔除。剔除非暗植被点对于最终气溶胶反演结果精度非常重要，实验发现很多极端情况，比如估算的气溶胶厚度超过2，容易出现暗植被的误判，并且使得最终反演的气溶胶光学厚度值出现局部范围的偏高。

在剔除非暗植被像元点之后，得到的气溶胶光学厚度数据还只是离散的点，为了获得整幅图像分布的气溶胶光学厚度数据，需要进行克里金插值与大尺度中值滤波平滑处理。克里金插值在固定插值半径的情况下被执行多次直到所有像素位置都获得气溶胶光学厚度，大尺度中值滤波指滤波半径大于10个像素的中值滤波，目的是获得边缘过度更加平滑的、更符合实际气溶胶空间分布特性的气溶胶光学厚度图像。克里金插值与大尺度中值滤波平滑处理过程见示意图5。气溶胶光学厚度图像重采样到原始图像尺寸并分块保存为float数据类型的tif图像磁盘文件，用户可选择作为中间结果保存，或者在大气校正结束后删除。

图1中“大气校正”部分的“大气校正”。是利用大气校正参数查找表与气溶胶分布数据，逐像元查找表插值出大气校正参数，并代入大气校正公式获得地表反射率数据。该过程是逐像元大气校正，是对于GF-4每一个波段的表观反射率数据，逐像素在气溶胶光学厚度数据中找到对应的气溶胶光学厚度值，根据气溶胶光学厚度值从该波段的大气校正参数查找表中从该气溶胶光学厚度值所在的前后两个气溶胶光学厚度值对应的大气校正参数中线性插值出3个大气校正参数，与表观反射率一同代入大气校正公式计算得到该像元的地表反射率。相比较传统大气校正的查找表多维插值的方法，该过程简化为只针对不同气溶胶光学厚度的单维线性插值，将一般单幅图像大气校正的时间由30分钟提高到5分钟以内。考虑到GF-4图像数据尺寸大，难以整个读入内存，表观反射率与气溶胶光学厚度的读入、地表反射率的输出都是采用分块处理的方式进行。

图1中“输出”部分的“地表反射率”。最终结果保存为tif图像磁盘文件，用户可选择输出float类型，或者乘以缩放系数后保存为int16类型节省一半存储空间。

本发明在青年科学基金项目的支持下(批准号41601384)，程序实例已经在PC平台上实现，目前已经交付用户方进行测试与使用，作为GF-4数据辐射预处理中气溶胶与地表反射率反演关键技术。

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。

Claims (1)

1.一种在线计算查找表的逐像元大气校正方法，该方法针对静止卫星高分辨率遥感图像辐射预处理应用，其特征在于包括以下实施步骤：

A数据准备，辐射定标与计算表观反射率，通过参数解析自动获得大气模式、气溶胶模式与平均高程，在6S软件包中添加卫星多光谱数据各波段光谱响应数据；所述辐射定标与计算表观反射率，目的是将遥感图像原始的Digtial Number值转换到表观反射率；所述通过参数解析自动获得大气模式、气溶胶模式与平均高程，是通过卫星数据XML辅助文件中获取的成像时间与地理范围，结合低分辨率的全球高程数据确定平均高程，结合低分辨率的全球地表分类数据确定大气模式与气溶胶模式，其中全球地表分类数据分类类别包含裸地、水体、城市、荒漠、高原积雪5类；所述在6S软件包中添加GF-4卫星多光谱数据各波段光谱响应数据，采用的6S版本为6SV2.1，通过修改6S源码的方式添加GF-4卫星多光谱数据各波段光谱响应数据并重新编译，使得6S软件包原生支持GF-4卫星数据，避免每次在线调用重新定义光谱条件；

B在线生成查找表，利用自动获取的太阳天顶角、观测天顶角、相对方位角、成像时间、大气模式、气溶胶模式与平均高程，针对不同的气溶胶光学厚度，在线调用6S生成针对当前图像的气溶胶光学厚度查找表与大气校正参数查找表；所述针对不同的气溶胶光学厚度，是在常用气溶胶光学厚度取值范围[0,3.5]内，根据查找表建表插值精度确定气溶胶光学厚度的个数，气溶胶光学厚度查找表取值为(0,0.05,0.1,0.2,0.5,1,1.5,2)共8个值，大气校正参数查找表取值为(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,2,2.5,3,3.5)共20个值；所述气溶胶光学厚度查找表，内容为对于参与气溶胶估计的波段，每个气溶胶光学厚度值对应三个由在线调用6S获取的变量值T、ρ₀、S，采用的辐射传输方程为

其中ρ_s为地表反射率、ρ_TOA为表观反射率、T为大气透过率、ρ₀为大气的路径辐射项等效反射率、S为大气下界的半球反射率；所述大气校正参数查找表，内容为对于多光谱图像的所有波段，每个气溶胶光学厚度值对应三个由在线调用6S获取的变量值a、b、c，采用的大气校正公式为y＝aρ_TOA-b，

C气溶胶光学厚度反演，利用气溶胶查找表，结合图像波段特性与修正的DDV(darkdense vegetation)法反演暗植被像元处的气溶胶光学厚度，剔除非暗植被点，再通过克里金插值与大尺度中值滤波平滑获取整幅图像的气溶胶分布数据；所述结合图像波段特性与修正的DDV(dark dense vegetation)法反演暗植被像元处的气溶胶光学厚度，是利用在暗植被处像元红波段地表反射率

与蓝波段地表反射率

存在固定的线性关系

红蓝地表反射率比率k取值为2，计算方法是利用表观反射率计算归一化植被指数

将NDVI值大于0.6的像元标记为暗植被，对于每一个暗植被像元，将红蓝波段气溶胶光学厚度查找表8个不同的气溶胶光学厚度值所对应的T、ρ₀、S参数带入辐射传输方程，由最满足该线性关系的两个气溶胶光学厚度之间插值出该暗植被像元处的气溶胶光学厚度值；所谓剔除非暗植被点，是考虑到暗植被的选取是在有大气的影响下确定的，当气溶胶光学厚度过大时可能会出现误判，剔除误判的方法是将所有反演得到的气溶胶光学厚度代入红波段与近红外波段的大气校正参数查找表，插值出大气校正参数，代入大气校正公式得到地表反射率，根据地表反射率重新计算NDVI，将大气校正后NDVI小于0.6的像元反演的气溶胶光学厚度值剔除；所述通过克里金插值与大尺度中值滤波平滑获取整幅图像的气溶胶分布数据，是将只在暗植被处反演得到的气溶胶光学厚度值插值到整个图像的过程，克里金插值在固定插值半径的情况下被执行多次直到所有像素位置都获得气溶胶光学厚度，最后利用大尺度中值滤波进行平滑处理；所述大尺度中值滤波指滤波半径大于10个像素的中值滤波，目的是获得边缘过度更加平滑的、更符合实际气溶胶空间分布特性的气溶胶光学厚度图像；

D逐像元大气校正，利用大气校正参数查找表与气溶胶分布数据，逐像元查找表插值出大气校正参数，并代入大气校正公式获得地表反射率数据；所述逐像元查找表插值出大气校正参数，是对于GF-4每一个波段的表观反射率数据，逐像素在气溶胶光学厚度数据中找到对应的气溶胶光学厚度值，根据气溶胶光学厚度值在该波段的大气校正参数查找表中从该气溶胶光学厚度值所在的前后两个气溶胶光学厚度值对应的大气校正参数中线性插值出3个大气校正参数，与表观反射率一同代入大气校正公式计算得到该像元的地表反射率。

Images