资源一号02D高光谱影像内陆水体叶绿素a浓度反演

刘瑶，李俊生，肖晨超，张方方，王胜蕾

1.自然资源部国土卫星遥感应用中心,北京100048;

2.中国科学院空天信息创新研究院 数字地球重点实验室,北京100094;

3.中国科学院大学 电子电气与通信工程学院,北京100049

1 引 言

叶绿素a（Chla）是浮游植物（藻类）中最常见的色素，叶绿素a浓度则是评价水体富营养化程度的核心水质参数之一（Cui 等，2020）。与传统地面测量方式相比，基于卫星遥感手段获取水体叶绿素a浓度具有空间连续观测和可长时序动态监测等优点。目前，面向大洋清洁水体的叶绿素a反演方法相对成熟，基于蓝绿波段比值（O’Reilly 和Werdell，2019）或海色指数OCI（Ocean Color Index）（Hu 等，2019）即可实现全球大洋水体叶绿素a 的高精度反演；但对于近岸和内陆浑浊水体，由于其光学特性复杂且随区域和季节变化大，实现精确的叶绿素a 浓度反演仍有较大难度（张兵等，2021；Neil等，2019）。

用于叶绿素a浓度反演的卫星遥感数据主要有3类：（1）波段宽度较大的多光谱数据，如Landsat TM/ETM+/OLI、Sentinel-2 MSI 等；（2）波段宽度较窄且数量较多但光谱离散的高光谱数据，如Envisat MERIS、Sentinel-3 OLCI 等；（3）波段宽度较窄、数量较多且连续的高光谱数据，如EO-1 Hyperion、GF-5 AHSI 等。然而，这3 类卫星数据应用于内陆水体的叶绿素a浓度监测都有着各自的局限性。首先，波段较宽的多光谱数据难以精确捕捉内陆水体的叶绿素a的光谱特征，反演叶绿素a 浓度的精度往往不高而且稳定性也较差（Gholizadeh 等，2016）。其次，MERIS、OLCI等波段离散的高光谱遥感器能够更精确地捕捉叶绿素a的光谱特征，进而提高叶绿素a的反演精度，但其空间分辨率往往较低，因此在中小型的内陆水体监测上能力有限（Kravitz 等，2020）。相比较而言，波段连续的高光谱数据在空间分辨率、光谱分辨率和波段数量上都具有优势，使其比前两类遥感器更适宜作为内陆水体叶绿素a 监测的数据源，不过其有限的幅宽和重访能力以及较低的信噪比（Giardino 等，2007）限制了其在内陆水体的叶绿素a反演方面的广泛应用。

2018-05，中国发射了高分五号卫星（GF-5），其搭载的可见短波红外高光谱相机AHSI（Advanced Hyperspectral Imager）有330 个波段，能够在可见近红外VN（Visible to Near-infrared）和短波红外波段SW（Short-Wave infrared）分别获取5 nm 和10 nm光谱分辨率的影像，且具有30 m的空间分辨率和60 km 的幅宽（Liu 等，2019）。紧接着在2019-09，中国又成功发射了资源一号02D 星（简称ZY-1 02D 卫星）。该卫星搭载了新一代AHSI 相机，具有与高分五号AHSI 相同的空间分辨率和幅宽，不同的是为了提高数据的信噪比，VN 和SW波段的光谱分辨率分别降低为10 nm 和20 nm，因此波段数也减少到166 个波段（刘银年等，2020a）。预计2021年还将发射ZY-1 02D 的后续星，双星组网将进一步提高其搭载的AHSI 的图像获取效率。而且，ZY-1 02D 高光谱数据的信噪比（刘银年 等，2020b）优于以往的EO-1 Hyperion、HJ-1 HSI 等高光谱数据，且幅宽更宽，覆盖能力更强。因此，在内陆水体的定量信息提取尤其是叶绿素a反演方面，ZY-1 02D高光谱数据具有很高的应用潜力，但是其具体应用效果还未见研究报道。

因此，本文基于ZY-1 02D 高光谱数据和同步获取的水面光谱和叶绿素a 浓度数据，开展ZY-1 02D 数据的叶绿素a 反演能力分析。开展基于不同波段遥感反射率组合的典型内陆水体叶绿素a浓度反演模型构建和反演精度分析，并评价ZY-1 02D数据在叶绿素a监测上的优势和不足。

2 研究区试验数据获取

2.1 研究区介绍

本文选择了营养状态和浑浊程度不同的太湖、小浪底水库和于桥水库作为研究区。太湖位于江苏省南部和浙江省北部交界处，面积约2300 km2，平均水深约2 m，属于浅水湖泊。近几十年，受到周围高强度的人类活动和环境变化影响，太湖的富营养化问题十分突出，叶绿素a浓度常年处于较高水平。小浪底水库位于河南省洛阳市与济源市之间，是黄河中下游库容量较大的峡谷型水库，水库面积约272 km2（武俐等，2020）。经过近些年来的治理和保护，小浪底水库水质显著好转，水体处于中营养级别（生态环境部，2020）。于桥水库位于天津市蓟州区，作为天津市的重要水源地，其主要功能是城市供水以及防洪。受水库上游及周边经济活动的影响，进入输水河道的污染物大量增加，导致于桥水库的富营养化趋势明显（岳昂等，2020）。

2.2 水面实验数据获取

2020-09—2020-11 期间，在太湖、小浪底水库和于桥水库分别开展了与ZY-1 02D 卫星的同步/准同步实验，共布设了46 个地面采样点（图1），在现场完成了光谱测量、透明度测量和表层水样的采集工作。表层水样于试验当天低温保存，运回实验室后依照热乙醇—分光光度计法（陈宇炜等，2006）测量叶绿素a浓度。各水面实验测量的叶绿素a浓度统计结果如表1所示。

图1 与ZY-1 02D卫星同步/准同步的研究区水面实验采样点分布图Fig.1 Sampling sites of in situ experiments in study areas during the ZY-1 02D satellite overpasses

水上实验时，按照“水面以上法”（唐军武等，2004）采用ASD Field Spec Pro 光谱仪在采样点进行光谱测量。采样点的遥感反射率通过测量参考板辐亮度（Lp（λ））、天空辐亮度（Lsky（λ））和水面辐亮度（Lu（λ））后计算得到。水面辐亮度测量时的观测天顶角为倾斜向下40°，观测方位角与太阳方向的夹角是135°；天空辐亮度观测天顶角为倾斜向上40°，观测方位角与测量水面辐亮度时一致。参考板辐亮度是垂直向下对准参考板测量。水面测量的遥感反射率（）计算公式如下：

式中，ρsky为天空光在气水界面的反射率，可以根据查找表确定（Mobley，1999），ρp为参考板反射率，由实验室内定标得到。3 个研究区水面测量的采样点遥感反射率，重采样到ZY-1 02D波段。

2.3 ZY-1 02D图像数据获取

太湖、小浪底水库和于桥水库的ZY-1 02D 高光谱数据获取日期分别为2020-09-06，2020-10-21 和2020-11-08。与表1 的试验时间比较可知，太湖和于桥水库的卫星影像和水面试验为同一天获取，小浪底水库开展水面试验比卫星影像获取早1 d。考虑到小浪底水库是深水水库，短期气象条件的改变对水质影响较小，可认为水面试验和影像获取时的水体光学特性基本不变。3个研究区的ZY-1 02D 卫星影像假彩色图像如图2 所示，其中，单景ZY-1 02D影像覆盖了太湖的大部分区域（图2（a））。影像的拍摄时间均为当地时上午11 点左右，影像质量良好，研究区水体上方无云。

表1 与ZY-1 02D卫星同步/准同步获取的典型研究区水面叶绿素a浓度数据统计Table 1 Minimum，maximum，and average values of measured Chla concentrations in study areas during the ZY-1 02D satellite overpasses

图2 与地面试验同步/准同步获取的ZY-1 02D高光谱影像假彩色合成图（红：860 nm，绿：654 nm，蓝：550 nm）Fig.2 False-color composite of concurrent ZY-1 02D hyperspectral images（Red：860 nm，Green：654 nm，Blue：550 nm）

3 方 法

3.1 ZY-1 02D水体图像大气校正及精度评价

（1）ZY-1 02D 图像预处理方法。ZY-1 02D 高光谱数据的预处理主要包括大气校正和水体范围提取。由于ZY-1 02D 的高光谱相机在VN 和SW 波段是分别成像，首先需要将VN 和SW 波段的影像数据进行合并，并通过头文件数据为合并后影像添加中心波长、半波宽和辐射定标参数信息。然后，基于辐射定标参数将DN 值转换为辐亮度图像，并转换为BIL（Band-interleaved-by-line）或BIP（Band-interleaved-by-pixel）的存储格式。最后，采用ENVI 的FLAASH 模块实施大气纠正。FLAASH 模块中，大部分参数参照官方文档进行设置；除此之外，对于ZY-1 02D 卫星数据，传感器高度设置为778 km，水汽反演选择940 nm 波段，spectral polishing的波段数设为3。

水体范围的提取采用光谱指数法。具体是基于FLAASH处理后的地表反射率图像进行归一化水体指数NDWI（Normalized Difference Water Index）计算，然后基于NDWI图像的直方图，设定阈值完成（Liu 和Xiao，2020）。其中，针对太湖所在的影像，还采用了基于浮游藻类指数FAI（Floating Algae Index）阈值法（Hu 等，2010）去除水华和水草区域，即，把FAI＞-0.004 的区域认为是水华或水草并去除。

（2）水面遥感反射率估算方法。FLAASH 处理后的地表反射率图像，需要去除天空光的影响，才能获得遥感反射率。面向ZY-1 02D 数据，本文采用了一种基于图像的遥感反射率近似估计方法（Wang等，2016），计算公式如下：

式中，(λ)代表图像上近似估计的遥感反射率，ρ表示地表反射率，min(ρNIR，ρSWIR)表示近红外和短波红外地表反射率的较小值。考虑到高光谱数据受噪声影响较大，ρNIR和ρSWIR分别取730—760 nm和1530—1630 nm 范围内ZY-1 02D 各波段反射率的平均值。

（3）水面遥感反射率精度评价方法。由于ZY-1 02D数据波段数量多，本文选取了叶绿素a浓度反演模型中使用的波段进行精度评价，具体评价结果见4.1节。精度评价方法是以地表实测的为真值，将图像估算的与其进行比较；采用的精度评价指标包括平均无偏相对误差AURE（Average Unbiased Relative Error），平均反射率比值和决定系数R2。其中，AURE和的计算公式为式中，N表示采样点的数量。AURE的计算中，考虑到图像和地表存在空间尺度和时间上存在差异，分母采用来减轻这种不一致性对误差的影响（Hu等，2012；Li等，2019）。

3.2 叶绿素a反演建模及其精度评价

（1）叶绿素a反演建模策略。目前，内陆水体的叶绿素a浓度反演主要有半经验模型、半分析模型、机器学习等方法。半经验模型在生物光学理论模型（Gordon 和Morel，1983）的基础上，基于特定假设提出与叶绿素a浓度相关的光谱指数，并建立统计关系实现叶绿素a 浓度反演（Gitelson 等，2007；Le 等，2013；Salem 等，2017）。半分析方法是从反射率信号中获取固有光学量，然后在辐射传输方程求解过程中实现叶绿素a 浓度估算（Liu 等，2020a）。机器学习模型则是将反射率和叶绿素a浓度数据输入网络训练，由网络自主学习二者之间的非线性关系，实现叶绿素a 浓度预测（Doerffer 和Schiller，2007；Pahlevan 等，2020）。其中，半分析模型的构建需要实测的吸收系数等固有光学量参数，而神经网络模型的有效训练又需要大量配套的反射率和叶绿素a浓度数据。相对地，半经验模型采用的光谱指数物理意义明确，且对实测数据的要求相对简单。因此，本文选择半经验模型进行ZY-1 02D 数据的叶绿素a 浓度反演模型构建。

建模数据采用实测叶绿素a 浓度和对应ZY-1 02D 像元的遥感反射率。匹配采样点的对应像元时，除了保证像元与地面点的位置保持一致，还需要保证像元具有较好的均一性：以像元周围3×3 邻域像元值的标准差与均值的比值不超过0.4为准。此外，鉴于星地同步的点位数有限，如果将全部同步点分为建模和验证数据集，建模样本可能会缺乏代表性。因此，本文采用留一点法交叉验证LOOCV （Leave-One-Out Cross Validation）完成建模和精度评价（Feng 等，2015；Li 等，2019），即把N个样本中的1 个样本作为单独验证集，其余N-1个样本用于建模，然后循环N次。

（2）半经验模型中采用的光谱指数。本文选择了5 种典型的光谱指数（表2），开展基于ZY-1 02D 数据的叶绿素a 反演半经验模型构建。首先，通过循环优化（Le 等，2013；宋挺等，2017）确定各光谱指数经计算所用的ZY-1 02D 波段：即针对每种光谱指数进行反演模型构建和验证，建模时只改变模型中1个波段的波长，其他波段保持不变，选择验证时叶绿素a反演误差最小的模型所采用的波长，作为模型的最终波长。对于包含多个波段的模型，依次采用这种方式确定各波段的最佳波长。基于此方法，确定计算这5种光谱指数的λ1-λ4分别是中心波长为671 nm，705 nm，731 nm和748 nm 的ZY-1 02D 波段。表2 中前4 个模型与ZY-1 02D 的波长差异主要是由于传感器波段设置差异造成的，而Wynne 等（2013）与本文使用的波长差异是因为计算的基高位置不同。

表2 叶绿素a半经验模型采用的光谱指数Table 2 Spectral indices used in semi-empirical models of Chla retrieval

（3）叶绿素a 反演精度评价指标。为筛选在ZY-1 02D 数据上应用效果最佳的叶绿素a 反演模型，将基于ZY-1 02D 图像遥感反射率优化后的上述模型，进行叶绿素a浓度反演精度评价。对验证集进行精度评价的指标有3项，包括R2，AURE和均方根误差RMSE（Root Mean Square Error），AURE 和RMSE的计算公式如下：

式中，N为样本数量，Xi和Yi分别表示叶绿素a 浓度的实测值和反演值。

4 结 果

4.1 ZY-1 02D大气校正精度评价

按照3.2 的匹配标准，46 个采样点都符合建模对数据的要求，因此全部用于叶绿素a 模型构建。基于3.1中的方法，估算得到的46个采样点在图像中的水面遥感反射率光谱如图3所示。

图3 ZY-1 02D图像大气校正的采样点遥感反射率光谱Fig.3 Atmospherically corrected ZY-1 02D image remote sensing reflectance of in situ sampling sites

对建模使用的4 个波段开展46 个同步采样点所在图像遥感反射率的精度评价，各项精度评价指标数据见表3，各波段的散点图如图4所示。

从表3中各波段的精度指标可以看到，671 nm和705 nm波段的图像遥感反射率精度较高，AURE分别为22.1%和25.8%，都在0.8左右，与的R2也很高，说明这2个波段图像估算的遥感反射率与地表实测值一致性较高；731 nm 和748 nm波段的精度相对低一些，AURE 超过36%，为0.67 和0.74，R2也相对低一些。这主要是由于，671 nm 和705 nm 波段反射率相对高一些，而731 nm 和748 nm 波段的反射率相对比较低，因此731 nm 和748 nm 波段相对于671 nm 和705 nm波段更易受到数据噪声的影响，因而731 nm 和748 nm 波段图像大气校正得到的遥感反射率的相对误差比671 nm和705 nm波段更大，也就是671 nm和705 nm 波段的大气校正精度更高。从图4 也可以看到，671 nm 和705 nm 波段的图像遥感反射率与水面实测数据一致性更高。

表3 ZY-1 02D图像大气校正的遥感反射率精度Table 3 Accuracies of atmospherically corrected ZY-1 02D image remote sensing reflectance

图4 同步采样点（N=46）在ZY-1 02D图像校正的遥感反射率与水面实测遥感反射率散点图Fig.4 Scatterplots of ZY-1 02D image-corrected Rrs versus in situ measured Rrs at sampling sites（N=46）

4.2 叶绿素a反演模型精度评价

基于地面测量的Chla 浓度和对应像元的，分别采用3.2 节中介绍的5 种模型进行反演模型构建和精度验证，获得的模型公式、R2和验证集上的误差指标如表4 所示。由表4 的建模R2可知，BR、NDCI、TBI、ETBI 与实测Chla 浓度均具有较好的相关性；BH 与实测Chla 浓度的相关性较差。但是，基于ETBI 的模型在应用时不稳定，出现少数采样点的Chla浓度反演值为负。

表4 基于同步ZY-1 02D图像遥感反射率和实测叶绿素a数据建立的叶绿素a反演模型和精度评价结果Table 4 Chla retrieval models calibrated using ZY-1 02 image and in situ measured Chla concentrations and accuracy analyses for these models

表4 基于同步ZY-1 02D图像遥感反射率和实测叶绿素a数据建立的叶绿素a反演模型和精度评价结果Table 4 Chla retrieval models calibrated using ZY-1 02 image and in situ measured Chla concentrations and accuracy analyses for these models

序号1 2 3 4 5模型简称BR NDCI TBI ETBI BH模型参数x Rrs(705)Rrs(671)Rrs(705) - Rrs(671)Rrs(705) + Rrs(671)(1 Rrs(671) -1 Rrs(705) )·Rrs(731)(1 Rrs(671) -1 Rrs(705) )/(1 Rrs(748) -1 Rrs(705) )Rrs(705) - Rrs(671) -705- 671 731- 671 (Rrs(731) - Rrs(671))建模模型y y=45.34x-32.04 y=87.06x2+92.77x+13.35 y=55.25x2+108.9x+13.36 y=80.03x+14.63 y=124.54x0.35 R2 0.78 0.78 0.77 0.72 0.50检验R2 0.76 0.75 0.75 0.00 0.35 AURE/%13.5 13.7 15.9—22.5 RMSE/（mg/m3）4.5 4.7 4.7—9.0

其中，基于BR 的经验模型在验证集取得了最小的AURE 和RMSE 误差，其建模的散点图和在3 个研究区各自的反演结果和精度如图5 所示。从图5 中的反演结果和反演精度可以看到，BR 模型在3个区域都具有较好的适用性，反演结果与实测值一致性较好。

图5 基于波段比值的叶绿素a浓度反演模型和反演结果散点图Fig.5 Chla retrieval model based on the band ratio and the scatterplot of estimated Chla concentrations

4.3 研究区叶绿素a反演结果图像

通过上述对不同模型建模后的精度分析可以看到，BR模型在ZY-1 02D高光谱数据上具有最佳的叶绿素a 反演效果。因此，基于该模型，对太湖、小浪底水库和于桥水库的ZY-1 02D 图像进行叶绿素a浓度反演，得到的叶绿素a浓度分布如图6所示。

图6（a）中，西太湖的红色区域为水华，而东太湖主要是水草。从图6（a）可以看到，基于ZY-1 02D图像反演的太湖叶绿素a结果空间分布状况与已有研究（Liu 等，2020b；Zhang 等，2019；Shi等，2017；宋挺等，2017）一致性很高：西太湖区域藻华严重，靠近藻华的区域叶绿素a浓度也比较高，这是由于太湖西部是入湖河流集中的地方，随河流汇入带来了大量营养盐，导致藻类生长旺盛，从而引发蓝藻水华和叶绿素a 浓度升高；湖心和东太湖区域受到的影响相对小一些，因此水质状况好于西太湖（Li等，2019），东部有水草生长，且叶绿素a浓度较低。图6（b）的反演结果显示小浪底水库主要区域叶绿素a 浓度都在10 mg/m3以下，与小浪底水库处于中营养状态的相关研究一致（黄新莹等，2020；琚会艳和李琳，2017）。图6（c）中于桥水库的西北岸叶绿素a 浓度最高，这与岳昂等（2019）发现的西北岸水质较差且经常发生藻华的情况相符，也与透明度北高南低（殷子瑶等，2021）的空间分布情况一致。上述对比分析表明，基于ZY-1 02D 高光谱数据反演的研究区叶绿素a浓度空间分布比较合理。

图6 基于ZY-1 02D图像反演的研究区叶绿素a浓度分布图Fig.6 Chla maps retrieved based on ZY-1 02D images in study areas

5 讨 论

综合试验结果与分析可以看到，ZY-1 02D 高光谱数据波段丰富，能够比多光谱数据提供更佳的叶绿素a反演波段，且空间分辨率较高，可用于中小型湖库的叶绿素a监测。然而，由于ZY-1 02D单星的重访周期长（大约55 d），大范围和多频次的水体监测需求难以得到满足。考虑到后续还将发射高光谱卫星，未来有望通过多星组网提升监测能力。另外，从反演结果中发现，ZY-1 02D 数据的叶绿素a反演结果会受到影像噪声和遥感反射率精度的影响。图6 中，小浪底水库的叶绿素a 反演结果中有明显噪声；这是由于小浪底水库的水体相对清洁，反射信号低，受到的高光谱数据噪声影响比较突出，因此需要发展针对于ZY-1 02D水体图像的降噪方法。此外，尽管本文中采用基于图像的水面遥感反射率估算方法取得了不错的效果，但更准确的遥感反射率获取，能进一步提高叶绿素a浓度反演的精度。因此，还需要进一步发展针对于ZY-1 02D 水体图像的高精度大气校正方法，比如检验和优化适用于清洁水体的近红外暗目标法和适用于浑浊水体的短波红外暗目标法等。

6 结 论

本文基于地面实测数据和同步获取的ZY-1 02D卫星高光谱图像，开展了面向ZY-1 02D 数据的内陆湖库叶绿素a浓度反演试验，分析了典型叶绿素a反演模型在ZY-1 02D 图像上的适用性和反演精度，并评价了ZY-1 02D 在叶绿素a 反演的优势和局限性，主要结论如下：

（1）基于ZY-1 02D 图像估算的遥感反射率可进行叶绿素a 浓度反演。构建的叶绿素a 反演模型中，波段比值模型（Rrs（705）/Rrs（671））的反演效果最佳。

（2）ZY-1 02D 高光谱数据波段丰富，能够比多光谱数据提供更佳的叶绿素a反演波段，且空间分辨率较高，可用于中小型湖库的叶绿素a 监测。但ZY-1 02D 的重返周期较长，未来需要通过多星组网提升监测能力。未来还需要发展针对于ZY-1 02D水体图像的降噪和大气校正方法。

（3）考虑到ZY-1 02D 卫星运行时间较短，积累的地面同步数据相对有限，且大气校正精度也有待提升，ZY-1 02D 高光谱数据谱段丰富这一优势尚未完全展现。下一步的研究计划是积累更多不同光学特性的内陆水体地面试验数据，通过这些水体进行分别建模和验证，或者发展基于软分类的叶绿素a 反演模型，进一步提高ZY-1 02D 高光谱数据应用于叶绿素a浓度监测的能力。