结合NCEP再分析资料优化的6S模式大气校正方法

以我国高光谱遥感卫星——环境1号卫星为例，开展结合NCEP再分析资料辅助优化的6S 大气校正方法的分析。首先，考虑到高光谱图像缺少标准反射率产品的问题，利用最优化估计方法构建高光谱反射率曲线，并作为标准曲线，用于大气校正结果的验证。其次，基于6S大气校正理论，开展了大气校正的敏感性分析，确定了气溶胶光学厚度的敏感因素以及气溶胶类型、大气模式和大气温湿度对大气校正系数的敏感性。在此基础上，提出了NCEP再分析资料辅助优化的6S大气校正方法，利用NCEP再分析资料提供的大气温湿度廓线、水平能见度反演的550 nm气溶胶光学厚度等数据资料，优化6S模式的输入参数，得到准确的大气校正系数X_a、X_b和X_c，获得大气校正后的不同地物反射光谱曲线。最后，选取西安作为试验区，以水体为例，进行波谱曲线对比，利用标准曲线对校正结果进行精度评价。对比分析结果表明，NCEP再分析资料辅助优化的6S模式校正的地面反射率结果明显优于6S的大气校正结果，与标准曲线具有一致的反射率变化趋势，二者的相关系数达到0.8596，标准差低于0.0685，各波段地面反射率逐像元误差的平均值和标准差接近0.02，反映了利用NCEP辅助数据优化的6S模式对大气校正有着明显的改善作用，可提高6S 大气校正的地物反射率反演精度。

Abstract

Objective

Unlike panchromatic and multispectral images, hyperspectral images have a high spectral resolution that causes more difficult atmospheric correction. Additionally, the atmospheric correction methods and correction accuracy based on the 6S model have been catching attention. China has launched several hyperspectral imaging sensor payloads, such as Environment 1 satellite, Gaofen 5, and resource satellites. Meanwhile, these satellites play an increasingly important position in agricultural remote sensing, forestry remote sensing, environmental monitoring and other fields, and the atmospheric correction technology for these hyperspectral remote sensing satellites is constantly developing. Taking China's hyperspectral remote sensing satellite Environment 1 satellite as an example, we carry out an optimized 6S atmospheric correction method combined with NCEP reanalysis data to improve 6S atmospheric correction accuracy.

Methods

First, considering the lack of standard reflectance products for hyperspectral images, the hyperspectral reflectance curve is constructed by the optimization estimation method and is regarded as the standard curve to verify the atmospheric correction results. Secondly, based on 6S atmospheric correction theory, we carry out sensitivity analysis and determine the most sensitive factors of aerosol optical thickness and the sensitivity of the aerosol model, atmospheric model, and atmospheric temperature and humidity to atmospheric correction coefficient. On this basis, an optimized 6S atmospheric correction method combined with NCEP reanalysis data is proposed. The aerosol optical thickness at 550 nm, atmospheric temperature and humidity profiles, and other data provided by NCEP are adopted to optimize the input parameters of the 6S model. Meanwhile, accurate atmospheric correction coefficients X_a, X_b, and X_c can be obtained, and the reflectance spectral curves of different ground objects are thus obtained after optimized atmospheric correction. Finally, by choosing Xi'an as the test area, the spectral curve of the water body is compared, and the accuracy of the correction results is evaluated via the standard curve.

Results and Discussions

The reflectance results by the 6S model and NCEP are significantly better than those by the 6S model. Compared with the standard curve, they have the same trend in spectral reflectance, and the correlation coefficient between them can reach 0.8596 with a standard deviation lower than 0.0685 (Fig. 13). The average and standard deviations of pixel-by-pixel error of ground reflectance in each band are close to 0.02, which demonstrates that the optimized 6S model with NCEP data has obvious improvement on the atmospheric correction.

Conclusions

The absolute error of the reflectance curve and the standard curve obtained by the atmospheric correction of the 6S model optimized by the NCEP data is much lower than that of the 6S model, and the average absolute error of each band is also less than that of the 6S model. The correlation coefficients of the three characteristic bands are higher than 0.85, the standard deviation is less than 0.07, and the mean and standard deviations of the ground reflectance per pixel error in each band are close to 0.02. Additionally, the determination coefficient between the ground object reflectance curve and the standard curve obtained by 6S+NCEP data reaches 0.78, which is higher than that by the 6S model. Meanwhile, the spectral angle of the optimized 6S model is reduced by 2.3565 and less than that of the 6S model, which indicates that the corrected spectral curve of the optimized 6S model is closer to the standard data. In conclusion, the atmospheric correction method in the 6S model of NCEP-assisted data optimization for HSI hyperspectral images can effectively improve the atmospheric correction effect.

1　引言

随着成像光谱技术的发展，高光谱成像作为一项新的遥感技术，续写和完善了光学成像从全色经多光谱到高光谱的全部图像信息链。因此，高光谱遥感成像技术在植被定量遥感、遥感智能解译、天空一体化综合监测等领域发挥了越来越重要作用^［1-4］。

定量遥感的应用是基于表面参数的定量反演，大气校正是定量分析地表参数的前提，定量反演必须从大气-陆地混合信号中提取有用的可见信息^［5］。利用高光谱遥感图像进行定量分析，首先需要获得地物真实的光谱特性，它需要消除辐射传输过程中大气对遥感图像的影响。但是，不同于全色和多光谱图像，高光谱图像具有高光谱分辨率和多光谱通道，高光谱遥感受到大气的影响更加严重，大气校正难度也更大。因此，针对高光谱遥感图像大气校正的研究变得越来越重要^［6］。

国内外的研究人员陆续提出了多种大气校正方法。早先，普遍采用基于统计学模式的方法，主要有经验线性法、内部平均法、平场域法和对数残差法等^［7-10］。1994年，Ben-Dor等^［11］采用统计学模式法对光谱通道的数据进行大气校正，使用3个参数即目标面积比（TAR）、相对于原始数据集的子集大小（SSR）和校正质量比（CQR）分析了校正效果。2002年，中国科学院遥感应用研究所提出一种平场域自动搜索（AFFF）方法，并将其用于对业务化运行航天高光谱数据进行大气校正^［12］。随着辐射传输模式的不断发展，辐射传输法已成为大气校正的主流方法，国内外众多科研工作者先后利用LOWTRAN、MODTRAN、ATREM、ACORN、FLAASH、6S等辐射传输模式开展大气校正方法研究^［13-15］。2010年，Proud等^［16］对SEVIRI数据进行大气校正，优化大气校正方法，并与6S模式校正结果进行对比。2015年，Bassani等^［17］利用气溶胶的局部微物理特性，从HICO大气数据中准确检索高光谱离水反射率，发现气溶胶微观物理特性可提高大气校正的准确性。2021年，Jha等^［18］开发了一种基于开源6S RT代码的柔性大气补偿技术（FACT），对高光谱（AVIRIS-NG和Hyperion）和多光谱（Landsat-8 OLI和WorldView-3）遥感数据进行了大气校正，显著提高了校正精度。2021年，Zhou等^［19］利用MODIS产品和30 m ASTER全球数字高程模式数据调整相关参数，搭建了利用这些辅助数据优化的6S模式，不同程度地提高了各波段的精度。2022年，Zhang等^［20］在6S辐射传输模式的基础上搭建了水陆边界临近效应修正模式WL-AE，提高了6S辐射传输模式在水陆边界的大气校正精度。可见，以6S模式为主的大气校正方法和校正精度一直备受关注。随着高光谱遥感技术的发展，高光谱遥感受到大气的影响更加严重，大气校正难度也更大^［21-22］。我国目前已发射多颗高光谱成像传感器载荷，如环境1号卫星、高分5号、资源卫星等，这些卫星在农业遥感、林业遥感、环境监测等领域发挥着越来越重要的作用，针对这些高光谱遥感卫星的大气校正技术在不断发展中。

本文针对高光谱遥感成像的大气校正技术，以我国高光谱遥感卫星——环境1号卫星（以下简称“HJ-1A星”）为例，开展结合NCEP再分析资料辅助优化的6S大气校正方法分析，以提高6S模式对高光谱遥感数据的大气校正精度。基于6S大气校正理论，对大气校正的敏感性进行分析，获得不同输入参数对大气校正系数的敏感性结果。在此基础上，采用NCEP作为辅助数据优化的6S模式进行大气校正，获得大气校正后不同地物的反射光谱曲线。为了验证结果，使用相同位置的相近日期的一景Landsat数据，利用最优化估计方法构建高光谱反射率曲线，对大气校正结果进行比较和分析，证明了NCEP辅助数据优化的6S模式对大气校正效果具有改善作用。

2　大气校正理论、方法和评价

2.1　大气校正理论

大气校正是一种将遥感影像数据转换为地面反射率的方法。由于遥感图像的成像过程具有复杂性，传感器接收的电磁波能量与目标本身辐射的能量是不一致的。传感器输出的能量包含由太阳位置和角度条件、大气条件、地形影响和传感器本身的性能所引起的失真，这些失真不是地面目标本身的辐射，因此需要对这些因素进行校正或者消除。在考虑大气作用下，传感器接收的反射辐射能量^［13］L_sensor表示为

L_{s e n s o r} = L_{p a t h} + L_{r e f l e c t e d} + L_{a d j a c e n c y}

，（1）

式中：L_path表示太阳光经大气衰减后照射到地面，经地面反射后，又经大气第二次衰减进入传感器的能量；L_reflected为大气散射、反射和辐射的能量；L_adjacency为地面本身辐射的能量经大气后进入传感器的能量。

大气校正的主要步骤包括：辐射定标、大气辐射传输过程模拟和地面反射率估算。首先，进行辐射定标。将原始遥感数据转化为能量单位，确定每个波段的实际辐射值，转化公式^［13］为

R_{r a d i a n c e} = G_{g a i n} \cdot g_{D N} + d_{o f f s e t}

（2）

式中：G_gain为传感器增益；g_DN为原始的图像灰度值；d_offset为传感器偏移。然后，进行大气校正，即根据地球大气的特性和观测条件，模拟大气辐射传输过程，通过现有的大气校正方法去除大气散射和吸收产生的影响，估算出地面反射率。

一般情况下，大气校正的精度受到很多因素的影响，包括大气参数的准确性、气溶胶物理特性和地理环境等，因此，需要根据具体场景和需求选择合适的校正方法和参数，以获得准确的地面反射率信息。

2.2　基于辐射传输的大气校正方法

利用基于电磁波在大气中的辐射传输原理所建立的模式对遥感图像进行大气校正的方法，被称为基于辐射传输的大气校正方法。其中，6S辐射传输模式可用来模拟光束在太阳-地面目标-传感器的传输过程中所受到的大气影响，因此，6S辐射传输模式常被用来进行大气校正。基于6S辐射传输模式的大气校正原理^［13］可表示为

ρ_{T O A} (θ_{s}, θ_{v}, Δ ϕ) = t_{g} (θ_{s}, θ_{v}) \{ρ_{a} (θ_{s}, θ_{v}, Δ ϕ) + [e x p (\frac{- τ}{μ s}) + t_{d} (θ_{s})] \frac{ρ_{s u} t_{d} (θ_{v})}{1 - ρ_{e} S}\}

，（3）

式中： $ρ_{T O A}$ 、 $ρ_{a}$ 、 $ρ_{s u}$ 、 $ρ_{e}$ 分别表示大气顶反射率、大气反射率、地面反射率和均匀反射率； $Δ ϕ$ 为太阳方位与传感器方位的夹角；t_g为气体的总透射率，其考虑了不同气体的吸收； $τ$ 表示大气光学厚度； $t_{d} (θ_{s})$ 和 $t_{d} (θ_{v})$ 表示大气漫透射率；S为大气的球形反照率； $1 - ρ_{e} S$ 表示考虑了地表与大气之间的多次散射。

经6S模式模拟运行后，可得到辐照度、透射率等一系列大气参数，通过对这些参数的处理可以得到大气校正系数^［12］X_a、X_b和X_c，相应的表达式为

X_{a} = \frac{π \cdot d \cdot e \cdot f}{c o s (\frac{z \cdot π}{180}) \cdot g \cdot h}

，（4）

X_{b} = \frac{i \cdot h}{g \cdot f}

，（5）

X_{c} = j

，（6）

式中：d为函数滤波器；e为太阳光谱；f为向上的总散射率；z为太阳天顶角；g为总透射率；h为向下的总散射率；i为垂直辐照度；j为总的球面反照率。

根据大气校正系数可计算得到大气校正后的地面反射率^［13］R_ACR，表示为

R_{A C R} = \frac{k \cdot X_{a} + X_{b}}{1 + (k \cdot X_{a} + X_{b}) \cdot X_{c}}

，（7）

式中： $k$ 为原始图像经过辐射定标后得到的辐射值。

2.3　评价方法

校正后的地面反射率结果，可通过其与标准反射率产品的相关系数、相对误差、光谱角、决定系数、均方根误差等指标进行评价。相关系数用来度量两个变量间的线性关系，一般用字母 $r$ 表示，即

r (X, Y) = \frac{E_{C o v} (X, Y)}{\sqrt[]{E_{V a r} [X] \cdot E_{V a r} [Y]}}

，（8）

式中： $E_{C o v} (X, Y)$ 为X与Y的协方差； $E_{V a r} [X]$ 为X的方差； $E_{V a r} [Y]$ 为Y的方差。

光谱角的定义为两条光谱之间的夹角，通过计算两个向量间的夹角来度量光谱间的相似性，夹角越小，两条光谱越相似，其计算公式为

α = a r c c o s [\frac{\sum_{i = 1}^{n_{b}} t_{i}}{{(\sum_{i = 1}^{n_{b}} t_{i}^{2})}^{\frac{1}{2}} {(\sum_{i = 1}^{n_{b}} u_{i}^{2})}^{\frac{1}{2}}}]

，（9）

式中：n_b表示波段数目；t_i和u_i分别表示第i条测试光谱和参考光谱。

决定系数是用来衡量回归模型对数据拟合程度的指标，表征趋势线的可靠性，愈接近1表示拟合程度愈高，其计算公式为

R^{2} = 1 - \frac{\sum_{i = 1}^{n} (y_{i} - {\hat{y}}_{i})^{2}}{\sum_{i = 1}^{n} (y_{i} - {\bar{y}}_{i})^{2}}

，（10）

式中： $y_{i}$ 表示真实值； ${\hat{y}}_{i}$ 表示观测值； ${\bar{y}}_{i}$ 表示平均值。

均方根误差是预测值与真实值偏差的平方与观测次数n比值的平方根，用来衡量观测值与真实值之间的误差，其计算公式为

E_{R M S} = \sqrt[]{\frac{1}{n} \sum_{i = 1}^{n} ({\hat{y}}_{i} - y_{i})^{2}}

。（11）

3　最优化估计获得标准反射率曲线

一般来说，在反射率数据验证方面，可与卫星标准产品或者标准波谱库进行比较。但是，在标准波谱库中，水体的光谱曲线易受到地区、海拔、河水含沙量、叶绿素含量等因素的影响，植被的光谱曲线又与品种有关，同时还存在背景条件不匹配等问题；在卫星标准产品方面，Landsat、MODIS等卫星虽可提供多光谱的数据产品，但存在光谱数据少的局限性。本文主要以我国高光谱遥感HJ-1A星为例进行大气校正方法的研究。HJ-1A星由两颗光学卫星和一颗合成孔径雷达卫星组成，配备了宽覆盖CCD相机、红外多光谱扫描仪、高光谱成像光谱仪、合成孔径雷达等多种类型的先进传感器。其中，高光谱成像仪共设置有115个谱段，谱段范围为450~950 nm，平均光谱分辨率为5 nm，空间分辨率为100 m。

目前，还没有高光谱遥感卫星的标准光谱曲线以及标准产品。本文以HJ-1A星和Landsat标准产品为基础，使用最优化估计与最优化模拟结合的方法来构建高光谱曲线。最优化估计是指通过寻找最优解来确定一个或多个参数的值，使得目标函数达到最小值或最大值^［23］。其原理是通过优化算法迭代地调整参数值，使目标函数尽可能地接近真实观测值。在本文中，使用最优化估计确定大气校正系数X_a、X_b和X_c。由于HJ-1A星的高光谱图像光谱范围为450~950 nm，Landsat卫星的光谱范围为430~2290 nm，其C1标准产品仅能提供7个波段的反射率值，这7个波段分别为433~453 nm、450~515 nm、525~600 nm、630~680 nm、845~885 nm、1560~1651 nm和2100~2300 nm。因此，考虑到HJ-1A星和Landsat卫星的光谱范围和分辨率，首先提取4个波段（482.6、561.3、654.6、864.6 nm）的高光谱图像辐亮度和Landsat C1产品数据的地面反射率，然后通过二者的重采样和波段匹配，使用拟牛顿法求解4组点的无约束最小值。

在上述基础上，使用最优化估计模拟方法来构建450~950 nm谱段的高光谱反射率标准曲线。最优化模拟通过对系统进行建模和仿真，获得不同参数组合下的系统性能指标，并通过对指标的评估和优化，找到最优化的参数组合。通常采用随机搜索、遗传算法等方法来实现最优化模拟。本文使用最优化估计得到的大气校正系数，结合光谱响应函数来构建最优化模拟系统，并采用随机搜索方法进行模拟，从而得到模拟曲线。

最优化估计和最优化模拟的具体步骤如下：1）对HJ-1A星高光谱图像进行垂直条纹去除等预处理；2）基于Landsat和HJ-1A星的波段匹配，对高光谱图像进行重采样，使其与Landsat波段对应；3）进行最优化估计，提取4个波段（482.6、561.3、654.6、864.6 nm）的高光谱图像辐亮度和Landsat C1产品数据的地面反射率，根据式（7）通过最优化估计方法求解4组点的无约束最小值；4）对4组最优化的最小值进行对比，选取误差最小的一组作为大气校正系数（X_a2、X_b2、X_c2）；5）基于大气校正系数对原始高光谱图像的辐亮度进行运算，利用最优化模拟方法得到模拟曲线，并将其作为标准反射率曲线。最优化估计和最优化模拟的流程如图1所示。