大气校正算法对高光谱反演水体叶绿素a 浓度的影响

孟祥亮，冯建飞，付萍杰，∗，张家威，张雨煊，孟飞

(1.山东建筑大学测绘地理信息学院，山东 济南 250101；2.山东省生态环境监测中心，山东 济南 250101；3.山东科技大学测绘与空间信息学院，山东 青岛 266000)

0 引言

水体叶绿素a(Chla)是水体富营养化程度的重要指标之一，对湖库叶绿素a 浓度的广域定量监测和实时定点分析是监测湖库水体营养状态的基本需求。水体叶绿素a 浓度的监测难点之一在于其空间异质性，地面监测虽然精度较高，但局限于某一特定时刻，难以满足动态变化的监测。 遥感技术具有监测范围广、监测时序长的优势[1]，但大气中水汽分子、气溶胶等物质会消耗掉大部分的离水辐射信号(90%～95%)，且内陆水体空间布局分散，地理状况多变，其光学特性表现为复杂性和区域性[2－4]。 大气校正能够消除部分大气影响，还原目标地物的光谱信号，高精度的大气校正可以降低不同空间和不同时相影像大气校正误差对模型构建和应用带来的不确定性，是提升水质遥感监测模型时空移植性的前提，对于提高水色反演模型的准确性和普适性具有重要意义。

学者们针对水质反演的大气校正方法开展了多项研究。 杜挺等[5]采用HJ－1B 多光谱影像，对太湖流域进行的反演中发现光谱超立方体的快速视线大气分析(Fast Line－of－Fight Atmospheric Analysis of Spectral Hypercubes， FLAASH)和对太阳光谱中卫星信号的二次分析(Second Simulation of the Satellite Signal in the Solar Spectrum， 6S)算法效果较好，但快速大气校正( Quick Atmospheric Correction，QUAC)算法下的典型地物光谱出现失真现象。 商子健等[6]以归一化植被指数(Normalized Difference Vegetation Index，NDVI)为指标，验证了6S 大气校正模型在GF－1 多光谱影像水体反演方面的适用性，其效果优于FLAASH 大气校正。 曾群等[7]将不同的波段组合因子运用到大气校正模型中，发现对于浑浊、高动态特性的水体，FLAASH 优于QUAC。 潘岑岑等[8]利用Hyspex 高光谱数据的研究指出，基于统计学的QUAC 和经验线性法两种大气校正方法效果不稳定，导致其系数变化较大。 崔成岭等[9]指出6S 大气校正查找表精度较低的现状，通过实时创建查找表的方式对资源1 号02D(ZY－1 02D)高光谱影像做大气校正处理，总体精度与FLAASH 相近，且发现6S 对蓝光波段的校正不完全。 刘瑶等[10]测试了资源1 号02D 高光谱影像在内陆水体叶绿素a 浓度反演方面的适用性，其波段比值模型取得了最好的效果，并指出针对于ZY－1 02D 水体图像的降噪和大气校正方法是未来的发展方向。 高光谱影像的大气校正在针对不同地物的反演中亦表现出不同的效果[11－12]，了解地物反演的参数贡献有利于提高大气校正的精度和效率。 目前，针对高光谱影像在二类水体叶绿素a 浓度反演方面的大气校正方法相关研究较少，由于波段连续细微的高光谱数据对大气干扰极为敏感，需要更精准的大气校正处理，才能满足水体叶绿素a 浓度反演的数据需求。

文章结合内陆水环境研究对影像光谱高精度的需求及高光谱卫星遥感器的特点，以南四湖为研究区，以(ZY－1 02D)高光谱遥感影像为数据源，按照6S、FLAASH、QUAC 大气校正的特点，从影像外部大气产品校正、影像内部大气补偿参数校正和影像的光谱均值校正3 个角度分别对影像进行大气校正，进而提取多维光谱指数，并利用半分析和机器学习方法验证反演模型的反演精度。

1 数据源与实验区选取

1.1 研究区和水体样本采集

南四湖位于微山县境内，由微山、昭阳、独山和南阳等4 个湖泊南北相连组成，是南水北调东线工程重要的水源地[13－14]，湖内水生生物与水禽种类众多[15]，其水质状况对东线工程以及生物多样性影响巨大。 南四湖面积约为1 266 km2，平均水深为1.5 m[16]，其流域多年平均年降水量在700 mm 以上，入湖河流有50 多条[17]，且来水河流近九成集中于上级湖，外流入湖和湖内运作均对南四湖水质产生不同程度影响[18]。

根据卫星过境时间，按照均匀布点原则，在独山湖周边采集38 个表层水样，分布如图1 所示。 文章所用地图审图号为鲁SG(2023)031 号。

图1 南四湖位置及实测点分布示意图

1.2 数据源及其预处理

遥感影像数据为2021 年9 月12 日的ZY－1 02D 高光谱数据，其搭载的可见短波红外高光谱相机(Advanced Hyperspectral Imager， AHSI)传感器比高分5 号的信噪比更高，但波段数量减少到166 个，光谱分辨率也随之降低[19]。 原始影像数据经辐射定标处理，得到大气校正需要的表观反射率和辐射亮度值。

外业数据采集利用美国ASD 公司生产的FieldSpec 4 Hi－Res 便携式地物光谱仪，按照白板、水体、天空、白板的顺序，实测采样点水体高光谱数据，每点重复5 次。 瓶装湖面表层水样[20]，放入黑色不透光袋中保存，并用GPS 定位仪记录采集点位置信息。 返回室内避光处， 采用美国安诺牌ChloroTech 121A 手持式叶绿素测定仪测定水体叶绿素a 浓度。

2 理论和方法

2.1 大气校正模型

大气校正可以保留由影像地物重要组成成分差异导致的反射率的微小差别信息[8]。 针对不同地物的光谱异质性及周边环境的差异，衍生出的大气校正方法也具有各自的侧重性。 6S 模型侧重影像过境时的观测及大气条件参数，并基于此参数建立的查找表进行逐像元插值得到监测点的大气校正系数；FLAASH 模型则侧重于影像像素光谱特征，考虑周边像元对目标像元造成的“邻近效应”、光谱噪声等因素对目标地物校正；QUAC 模型则是整幅影像的光谱信息收集，监测点的光谱特征即为目标地物的光谱均值。

6S 辐射传输模型由5S 模型发展而来[21]，波段处理范围为0.25 ～4.0 μm，通过模拟信号传输过程的太阳辐射，并计算信号在进入传感器前的辐射能量，得到校正参数进行大气校正[7]。 ZY－1 02D 高光谱影像需进行可见光近红外波段(Visible and Near-Infrared，VNIR)和短波红外波段(Short Wave Infrared，SWIR)的数据整合，将表观辐亮度转换为表观反射率，由式(1)表示为

式中ρTOA为表观反射率； Lλ为波段λ 处的中心波长，nm；d 为日地距离，AU；Eλ为波段λ 处的太阳光谱辐照度，W/(m2·μm)；θ 为太阳天顶角，(°)。

6S 模型通过输入表观反射率和影像元文件，得到各波段大气校正系数xa、xb、xc。 将表观反射率转化为真实的地表反射率ρr，可由式(2)和(3)表示为

FLAASH 模型基于MODTRAN4 ＋辐射传输模型，波段区间0.4 ～2.5 μm。 通过设定模型参数，逐像元反演校正参数。 针对不同的波段区间，进行特定的水汽反演，采用暗目标法进行气溶胶厚度的反演。 影像像元光谱辐射亮度由式(4)表示为

式中L 为像元辐射亮度；ρ 和ρe分别为像元与相邻像元地表反射率均值；S 为大气球面反照率；Lα为大气程辐射； A、B 是依赖于大气(透过率) 和几何状况的系数。S、Lα、A、B 可由辐射传输模型MODTRAN4＋计算得到。

QUAC 是基于影像本身的统计方法，采集像元内的地物光谱值，取同一地物的光谱平均值作为判别地物的经验值，利用端元平均反射率与参考物进行对比确定大气的影响。 QUAC 大气校正算法并不依赖影像获取过程中的各类参数信息，大气参数和仪器标定的误差对校正结果的影响较小[7]。 端元平均反射率ρ′由式(5)表示为

式中ρ1，ρ2，…，ρn为视场内各物质的端元光谱反射率；n 为端元个数。

2.2 水体叶绿素a 浓度反演

(1) 三波段光谱指数

光谱指数是一种基于光谱反射率构建的指标函数，依据地物的光谱特性，通过波段组合度量地表参量[22]。 在尝试了多种三波段组合后，选择遍历后相关系数较高的公式反演三波段指数(Three-Band Index，TBI)，由式(6)表示为

式中Rλ1、Rλ2、Rλ3分别为395 ～900 nm 范围内波长为λ1、λ2和λ3处的遥感反射率。

(2) 四波段光谱指数

按照三波段光谱指数法选取波段组合的方法，运用水体反演叶绿素a 浓度的四波段公式反演四波段指数(Four-Band Index，FBI)，由式(7)[23]表示为

式中Rλ4为395～900 nm 范围内波长为λ4处的遥感反射率。

(3) CatBoost 算法

CatBoost 算法是梯度提升决策树(Gradient Boosting Decision Tree，GBDT)算法框架的改进算法，运用排序提升方法构建模型，在不同的迭代中会采用不同的排列顺序，训练集中的噪声点被消除，梯度估计与浓度预测误差得到了解决。 CatBoost 算法通过减少超参数调优来抵抗模型过度拟合，增强了算法的准确性和泛化能力，其计算过程添加了先验值和权重参数，一些低频率类型的数据和噪声点对数据整体趋势的影响得到有效控制。 CatBoost 算法可由式(8)表示为

式中xσj，k、xσp，k分别为第k 个训练样本的第j、p 个类别特征；xi，k为类别特征平均值；[]为指示函数运算，即括号内两个量相等时取1，否则取0；Yj为第j 个样本的标签； P 为添加的先验项；a 为先验值的权重。

(4) 反演模型构建

对大气校正的影像提取像元光谱，并基于叶绿素a 实测浓度，提取最优TBI 和FBI 用于CatBoost模型的特征组建。 CatBoost 模型采取相同的随机种子进行训练集(70%)与验证集(30%)的划分，网格搜索法在指定范围内进行枚举[24]，将效果最好的参数用于模型构建。

一元线性模型的建模参数分别为TBI 和FBI 的最优遍历结果[25]，训练集与验证集的划分同CatBoost 模型。

2.3 大气校正及叶绿素a 反演精度评价

为定量比较大气校正模型的校正效果，使用光谱角制图(Spectral Angle Mapper，SAM)[26]、均方根误差(Root Mean Square Error，RMSE)和相关系数来衡量不同模型校正后的影像光谱反射率与实测反射率之间的接近程度[21]。 SAM 光谱角α 的余弦值由式(9)表示为

式中N 为样本总数；Ai、Bi分别为第i 个像元向量的光谱值。

相关系数r 由式(10)表示为

式中Xi、Yi分别为像元光谱值和实测光谱值；和分别为像元、实测的光谱均值。

叶绿素a 浓度反演精度采用决定系数R2和均方根误差RMSE 评价。 RMSE 和R2分别由式(11)和(12)表示为

式中Xai、Yai分别为叶绿素a 浓度的实测值和反演值；yi、分别为实测值和模型反演值；为N 个yi的均值。

3 结果与分析

3.1 大气校正

3.1.1 大气校正下的光谱曲线对比

由于实测点位分布在独山湖区域，因此仅针对南四湖的独山湖区域展开分析。 不同大气校正参数见表1，观测日期为2021 年10 月22 日，过境时间为世界标准时间(Universal Time Coordinated， UTC)03:19:42，中心经纬度为35.28°N、116.80°E。 其中，FLAASH、QUAC 大气校正参数均在ENVI 中完成，由于没有相应的传感器类型，所以统一选择“unknown sensor”。

表1 大气校正参数表

为保证光谱信息的可信度，任取一点(14 号点)的实测值和全部监测点光谱均值，分别与大气校正光谱对比，结果如图2 所示。 由于叶绿素和类胡萝卜素的吸收，在440 和490 nm 附近形成两个小反射谷；藻类色素的低吸收以及水中悬浮物的散射，使得光谱曲线在570 nm 附近形成了大反射峰；4 种光谱曲线均在675 和700 nm 附近分别形成了大反射谷和大反射峰，这与前人的研究结果一致[27]。 对比图2(a)和(b)发现，除反射率略有升高外，光谱均值与14 号点光谱曲线的整体趋势基本一致。 由于受到水面天空光等多种信号的影响，大气校正曲线普遍高于实测光谱曲线。 校正方式的不同使得光谱反射率的值高于实测值的程度也不一样，6S 借助大气校正产品逐一计算得到每波段每像素的地表反射率值，光谱反射率最接近实测值；而FLAASH 是基于整张图像[4]，且重采样减少了光谱的波段数量，对于蓝光波段的不完全校正使得FLAASH 曲线出现部分负值。 QUAC 则是基于影像本身，且水体占比较大，有利于进行水体的光谱校正。

图2 不同大气校正光谱对比

3.1.2 精度评价

38 个采样点的实测光谱数据与大气校正数据的指标值如图3 所示。 6S 和QUAC 模型的光谱角余弦值都在0.9 以上，而FLAASH 模型的余弦值略低一些，且各点之间的差异较大。 同时6S 模型的r在3 种大气校正模型中最高，这也说明6S 模型对地物光谱信息的保持度高。 6S 与FLAASH 模型的RMSE 值均约为0.15，QUAC 的值与实测值的离散程度较大。 综上所述，6S 模型的评估结果均表现良好，FLAASH 模型的光谱信息保持较弱，QUAC 模型校正后的影像反射率缺乏代表性，校正效果不稳定。

图3 大气校正精度指标图

3.2 叶绿素a 浓度反演

3.2.1 多维光谱指数提取

在文章研究的光谱波段范围内(395～900 nm)，通过MATLAB 软件实现基于TBI 和FBI 公式的波段反射率循环迭代。 基于实测叶绿素a 浓度，利用相关系数r 提取最优的波段组合，计算模型参数。选取相关性最高的前5 个光谱指数组合作为CatBoost 反演模型的参数，选取相关性最优的光谱指数进行一元线性回归模型的反演。 TBI 和FBI 的相关结果分别见表2 和3。

表2 三波段光谱指数及相关系数

表3 四波段光谱指数及相关系数

3.2.2 叶绿素a 浓度反演

对影像光谱数据和实测叶绿素a 浓度，数据测试集和训练集分别占30%、70%。 将基于TBI 和FBI 构建的一元线性反演模型应用到遥感影像上，分别得到12 个点的反演值。 将不同算法的反演值与实测值进行拟合，结果如图4 所示。 6S 大气校正后得到的叶绿素a 浓度反演值，在以三波段和四波段作为特征参数的线性模型中，均具有良好表现。QUAC 大气校正与FLAASH 大气校正后的线性模型反演精度略有差异，其中QUAC 四波段线性模型精度最高。 整体来说，反演值与实测值的R2最高能达到0.74，模型效果较好。

图4 不同大气校正算法反演结果与实测结果拟合图

选择与一元线性模型相同的训练集，由于CatBoost 不需要过多调参，因此主要调节模型决策树的数量(n＿estimators)、学习率(learning＿rate)和决策树的深度(depth)，其余参数保持默认，按照参数调优标准，得到CatBoost 模型参数见表4。

表4 CatBoost 模型反演参数设置

AHSI 影像反演模型和精度评价结果见表5。CatBoost 模型的建模精度总体良好，测试集最高相关系数R2＝0.80、RMSE ＝4.97 μg/L。 由于CatBoost 模型在不同参数组合中取得的精度最高，且在6S 和QUAC 大气校正后的参数反演结果均较好，因此选用CatBoost 模型进行南四湖水体叶绿素a 浓度反演。

表5 AHSI 影像反演模型和精度评价结果

对比不同大气校正的模型建模结果，6S 大气校正的模型拟合效果最好，其中四波段参数CatBoost模型的拟合精度最高。 基于QUAC 大气校正的三波段参数和四波段参数组合的CatBoost 模型拟合精度均优于FLAASH 大气校正。 QUAC 大气校正效果缺乏稳定性，模型最低拟合精度R2＝0.35。

将模型应用到影像，得到南四湖叶绿素a 浓度反演结果如图5 所示。 可以看出，质量浓度空间差异明显，高值主要分布在南阳湖的东部沿岸和独山湖的湖心，流动性较弱区域，这与前人的研究结果相符[28]。 然而，该组合方法下模型的泛化能力较差，选取的波段对该地区的水体叶绿素a 缺乏敏感性，影像模糊是由于反射信号包含高光谱噪声。

图5 不同大气校正下的CatBoost 模型AHSI 影像南四湖叶绿素a 浓度

对比不同光谱指数的反演，四波段参数选择的波段反射率更接近实测反射率，参数信息更能体现叶绿素a 浓度变化。 这一点在6S 与QUAC 大气校正的影像中表现明显。

6S 大气校正的三波段与四波段参数模型反演，R2均稳定在较高水平，表明6S 大气校正对于高光谱反演水体叶绿素a 浓度具有良好的适用性，可作为影像预处理时的优先选择大气校正方法。

考虑到QUAC 模型以影像本身作为对象，对同类地物光谱采集取均值的校正特性，对比叶绿素a浓度反演结果分析，特征影像受地物异物同谱特性的影响较大。 由于受大气等光学条件的影响，不同水体叶绿素a 浓度监测点的光谱可能存在较大差异[5]，在求取均值后，不同叶绿素a 浓度对应的光谱值存在相近或相等的情况，这将对反演造成困难，三波段和四波段参数的光谱值出现误差，造成最终反演精度出现不确定性。

FLAASH 大气校正后的反演效果均不理想，可能是选取的TBI 与FBI 所用波段对叶绿素a 并不敏感[29]，FLAASH 模型以大气校正产品作为系数，不同波段参数的组合对影像的泛化能力也不尽相同，导致构建的模型不能很好地预测水体中的叶绿素a浓度。

4 结论

基于南四湖ZY－1 02D 影像，分别进行了6S、FLAASH 和QUAC 的大气校正处理，并对提取的光谱反射率进行多种组合，分别利用线性回归模型和CatBoost 模型中反演水体中叶绿素a 浓度，得到主要结论如下:

(1) 6S 模型校正效果最佳，QUAC 次之，FLAASH 模型效果最差；

(2) CatBoost 模型能在一定程度上提高反演精度，排序提升的模型构建以及先验值的添加消除了误差噪声，在不需要过多调参的情况下即可达到较好的回归预测精度；

(3) 6S 模型的四波段组合CatBoost 模型反演结果的R2最高为0.80，更有利于提高光谱数据在南四湖中部的独山湖叶绿素a 浓度反演中的真实性和反演精度。