基于HY-1C卫星CZI影像的多类型内陆水体大气校正及水质参数反演

张方方，李俊生，3，王超，王胜蕾，王正，张兵，3

1.可持续发展大数据国际研究中心, 北京 100094;

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

3.中国科学院大学, 北京 100049;

4.河南省科学院地理研究所, 郑州 450052

1 引 言

海洋一号C 星（HY-1C）和D 星（HY-1D）是中国第三和第四颗海洋水色系列卫星，分别发射于2018年9月7日和2020年6月11日，两颗卫星采用上、下午卫星组网的观测形式，可增加观测次数，提高全球覆盖能力。HY-1C 和HY-1D 卫星上分别携带了一台海岸带成像仪（CZI），主要用于获取海陆交互作用区域的实时影像资料进行海岸带监测。CZI 传感器兼具宽幅盖（950 km）、短重访周期（3 d）、中等分辨率（50 m）的特点，可以兼顾内陆水体及陆地监测的需求（苗珊珊，2018；蔡婷，2000）。

从2018年发射以来，HY-1C CZI 影像在很多领域得到了应用。主要包括冬小麦识别（王利民等，2019）、珊瑚礁监测（Zou 等，2019）海面溢油识别（沈亚峰 等，2020）、自然灾害监测（刘建强 等，2020）、红树林提取及长势监测（梁超 等，2020；邹亚荣 等，2020）、火山喷发（刘建强 等，2021a）、冰架裂缝变化与断裂过程监测（刘建强等，2021b）等。HY-1C CZI 数据在近岸和内陆水色遥感方面也得到了初步应用，包括蓝藻水华监测（Cao 等，2021），叶绿素a（Chen 等，2019；Cao 等，2021）、悬浮物（Cai 等，2020；Cao 等，2021）、浊度（周屈 等，2020）等的反演。其中，内陆水体的应用区域仅仅有太湖和武汉知音湖、黄家湖，二者都是高度浑浊水体，还有很多类型的内陆水体有待检验与扩展。

水色遥感，特别是水质参数反演属于定量遥感的范畴，精确大气校正是卫星数据应用的瓶颈。大气校正的前提是准确的辐射定标，HY-1C 星上配置了具备高光谱分辨率与高辐射定标精度的定标光谱仪，实现对同一卫星平台载荷设备的交叉定标（张可立 等，2019）。HY-1C CZI 影像缺少短波红外波段用于水体大气校正，目前的研究主要利用6SV 和Flaash 做大气校正，6SV 需要的气溶胶数据多从MODIS 气溶胶产品中获取（Tong 等，2019；Cai 等，2020；Cao 等，2021）。6SV 和FLAASH 是针对陆地设计的，对于仅有4 个宽波段的HY-1C CZI影像来说，其水体大气校正结果在绝对精度和稳定性上仍然不够理想。因此，也有部分学者利用HY-1C CZI 的瑞利散射校正产品进行水质参数反演研究（周屈 等，2020）。综上，针对HY-1C CZI影像的水体大气校正技术还不成熟，仍然需要尝试更多的方法提高大气校正的精度和稳定性。

针对HY-1C CZI 影像大气校正算法不成熟，水质参数反演应用较少的问题，本研究将发展基于均匀不变地物的HY-1C CZI 影像相对大气校正算法，拓展HY-1C CZI影像水体叶绿素a浓度和透明度反演应用到多个类型水体（浑浊、较浑浊、较清洁）。通过实测数据验证HY-1C CZI 影像在多类型内陆水体水色遥感的应用潜力。

2 研究区与数据

2.1 研究区

本文的研究区为小浪底水库、官厅水库、丹江口水库、白龟山水库、白洋淀（图1）。5个湖库都分布在华北平原，但是湖泊类型多样，具有较好的代表性。其中，白洋淀水体较为破碎，由多块独立的水体组成，并由水道相连，水深较浅，水体较为浑浊；官厅水库两端为开阔水域，中间为狭窄河道，水体也较为浑浊；小浪底水库和白龟山水库都是饮用水源地，水体较为清洁；丹江口水库为南水北调中线水源地，是五个研究区中最为清洁的水体。以上5个研究区涵盖了较浑浊、较清洁、清洁等不同水体类型，有利于检验HY-1C CZI数据在不同类型内陆水体水质遥感监测的适用性。

图1 研究区及采样点分布图Fig.1 Distribution of study area and sampling sites

2.2 数据获取

本研究在小浪底水库、官厅水库、丹江口水库、白龟山水库、白洋淀5 个研究区分别开展了1 次星地同步实验，一共获取了85 个采样点数据（如表1）。其中，56 个用于遥感反射率系统定标，29 个用于遥感反射率反演精度检验。由于部分采样点现场测量条件的限制、运输及室内测量存在问题，叶绿素a 和透明度数据少于85 个。其中，58 个采样点拥有叶绿素a 浓度数据（39 个用于建模，19个用于精度检验），72个采样点拥有透明度数据（53个用于建模，19个用于精度检验）。

表1 采样数据及影像数据明细表Table 1 Sampling data and image data

本文获取了5个研究区的HY-1C CZI影像，其中，白洋淀和丹江口水库的成像日期和采样日期为同一天，官厅水库、白龟山水库、小浪底水库的成像日期比采样日期晚一天。对HY-1C CZI 影像进行基于不变地物的相对大气校正需要Sentinel-2 MSI影像的辅助。每一景HY-1C CZI影像需配套邻近日期的两景Sentinel-2 MSI 影像，一共获取10 景影像。每个区域的两景Sentinel-2 MSI 影像的间隔时间为5—15天，可以方便提取不变地物。

另外，本文发展的基于均匀不变地物的相对大气校正算法的精度完全依赖Sentinel-2 MSI 影像的地表反射率产品（L2A）的精度。为了验证Sentinel-2 MSI 影像地表反射率产品的精度，本文在北运河（2019年10月15日，15 个采样点）、于桥水库（2018年11月1日，13个采样点）、玉渊潭（2019年8月21日，9 个采样点）、太湖（2020年9月5日，12 个采样点）等4 个研究区开展了与Sentinel-2 MSI 传感器的星地同步实验，获取了一共49 个采样点的实测遥感反射率数据。另外，本文用于HY-1C CZI建模和验证的5个研究区中的丹江口水库、白龟山水库、小浪底水库的49 个采样点也与Sentinel-2 MSI 影像同步，也用来检验Sentinel-2 MSI 影像地表反射率产品的精度。利用以上共计110个同步数据（野外实验和卫星成像日期均为同一天）评价Sentinel-2 MSI 影像地表反射率产品的精度以及在水色遥感研究中的可用性。

3 方 法

3.1 基于均匀不变地物的相对大气校正

地球表面时刻处于变化之中，有些地物变化剧烈，有些地物变化较为缓慢，变化缓慢的地物在一定的时间内可以视为不变地物（丁丽霞 等，2005）。分布均匀且面积较大的地物具有稳定的光学特性，且可以避免邻近效应的影响。变化缓慢、分布均匀且面积较大的地物即为均匀不变地物。可以通过建立均匀不变地物已知反射率与待校正影像像元之间的关系实现相对大气校正（Yuan 和Elvidge，1996；Du 等，2002；Canty 等，2004；郭丽峰 等，2009）。Sentinel-2 MSI L2A 数据是目前全球使用较为广泛的地表反射率产品，且时间分辨率较高。因此，本研究利用Sentinel-2 MSI L2A影像建立基于均匀不变地物的HY-1C CZI 影像大气校正方法。主要步骤如下（图2）：

图2 基于均匀不变地物的HY-1C CZI相对大气校正流程图Fig.2 Flow chart of HY-1C CZI relative atmospheric correction based on uniform and invariant ground pixels

（1）影像预处理：以Sentinel-2 MSI 影像为基准，对HY-1C CZI 影像进行几何配准；提取Sentinel-2 MSI影像与HY-1C CZI影像的对应波段：443 nm，560 nm，665 nm，842 nm。为了保持二者具有一致的空间分辨率，将Sentinel-2 MSI 影像重采样为50 m；并用式（1）将Sentinel-2 MSI L2A 的地表反射率（SR）数据转换为遥感反射率（Rrs）。

（2）基准影像均匀地物提取：分别求两景Sentinel-2 MSI 影像3×3 像元的变化率ρ，ρ＜5%的像元提取为均匀地物。ρ计算公式如下：

式中，SR(c)为中心像元地表反射率，SR(i)为非中心像元的地表反射率。

（3）基准影像不变地物提取：将计算（2）步中两景影像根据地理坐标信息进行叠加，计算两景图像重叠像元的变化率，变化率ρ′＜10%的像元提取为不变地物；

（4）基准影像均匀不变地物提取：将（2）提取的均匀地物和（3）提取的不变地物进行叠加，二者的重合区域即为均匀不变地物；

（5）目标影像均匀地物提取：参考（2）的方法计算目标影像HY-1C CZI上的均匀地物；

（6）均匀不变地物提取：将（4）的基准影像均匀不变地物和（5）的目标影像均匀地物叠加，求得最终的均匀不变地物；

（7）相对大气校正模型构建：分别提取（6）得到的均匀不变地物的HY-1C CZI影像像元值和对应的两景Sentinel-2 MSI 影像像元值的平均值，并建立二者之间的线性关系，得到相对大气校正模型，根据模型计算得到相对大气校正反射率Rrs'。

（8）基于实测数据的HY-1C CZI系统定标：由于HY-1C CZI 和Sentinel-2 MSI 传感器的波段设置和光谱响应函数存在一定差异（表2、图3），将会导致相对大气校正后HY-1C CZI 影像的遥感反射率和Sentinel-2 MSI 的相似度较高，但是和HY-1C CZI原本应有的数值存在系统性偏差。因此，本文将根据5个研究区的56个采样点的实测数据对相对大气校正后结果进行系统定标，根据二者的线性关系，建立系统定标模型（公式：Rrs = Rrs’*Gain），其中，Rrs’为7）得到的相对大气校正结果，Gain为校正系数，Rrs为最终的水体遥感反射率。

表2 HY-1C CZI和Sentinel-2 MSI传感器光谱通道参数比较Table 2 Comparison of spectral channels parameters between HY-1C CZI and Sentinel-2 MSI

图3 HY-1C CZI和Sentinel-2 MSI传感器光谱响应函数（蓝色、绿色、红色、黑色线条分别表示蓝、绿、红、近红外波段）Fig.3 Spectral response function of HY-1C CZI and Sentinel-2 MSI（blue，green，red and black lines represent blue，green，red and near infrared bands，respectively）

3.2 水质参数反演模型构建

本文尝试了多种波段比值和差值模型进行水体叶绿素a 浓度和透明度反演建模（Duan 等，2007；Tebbs 等，2013；Zhang 等，2021）。首先将实测遥感反射率使用HY-1C CZI 光谱响应函数进行波段等效，等效后实测光谱也是4波段数据，利用等效后数据进行水体叶绿素a 和透明度反演建模。

3.3 精度评价

本研究在5 个研究区获取了85 个采样点的水面遥感反射率光谱和水体叶绿素a浓度、透明度数据。遥感反射率数据较多，其中，29 个采样点数据用于遥感反射率精度验证，部分采样点未测量叶绿素和透明度，所以19 个采样点数据用于叶绿素a和透明度反演精度验证。采用平均无偏相对误差（AURE）、均方根误差（RMSE）验证大气校正得到的每个波段的遥感反射率光谱的精度；采用光谱相关系数（r）和光谱角度距离（SAD）验证大气校正得到的遥感反射率光谱形状与实测光谱形状的相似性。水体叶绿素a 浓度和透明度采用AURE和RMSE进行反演精度评价。公式如下：

式中，VE，i和VM，i分别指反演（或大气校正后光谱）数据和实测数据，N指验证数据量，Cov（VE，VM）指大气校正后光谱和实测光谱的协方差，Var［VE］和Var［VM］分别为大气校正后光谱和实测光谱的方差。

本文采用FLAASH（Fast Line-of-sight Atmospheric Analysis of Hypercubes）大气校正结果和本研究发展的基于均匀不变地物的大气校正算法作比较。FLAASH 是ENVI（The Environment for Visualizing Images）遥感图像处理软件集成的遥感影像大气校正程序，其内核算法为MODTRAN（MODerate resolution atmospheric TRANsmission）辐射传输模型。FLAASH 被广泛应用于各类高光谱和多光谱影像的大气校正。本研究将利用FLAASH 对HY-1C CZI影像进行大气校正，生成地表反射率数据，忽略水面反射的情况下，地表反射率可以粗略的用式（1）转换为遥感反射率，方便和本研究发展的基于均匀不变地物的相对大气校正结果做精度比较和分析。

本文使用FLAASH 对HY-1C CZI 影像做大气校正时，各参数设置如下：（1）输入为大气层顶辐亮度图像；（2）经纬度为ENVI 软件自动读取；（3）传感器类型为UNKNOWN-MSI，Sensor Altitude为782 km，Ground Elevation 为Google Earth 上获取的海拔数据（单位转换为km），Pixel Size 为50 m；（4）Flight Date/Time 从图像参数文件获取（格林尼治时间）；（5）大气模式选择了中纬度夏季，因为本文所用影像都是中国北方夏半年的数据；（6）气溶胶类型选择了Rural 乡村型；（7）因为HY-1C CZI没有水汽通道，因此Water Retrieval选择了No，不反演水汽参数；（8）HY-1C CZI 的光谱响应函数从卫星海洋应用中心获取；（9）其它参数均采用FLAASH默认参数。

4 结果与分析

4.1 实测光谱及水质参数

本文在5个研究区开展了水体遥感反射率（如图4）和透明度现场测量（表3），并采集水样，在实验室测量了水体叶绿素a浓度（表3）。从遥感反射率光谱曲线来看，蓝波段反射率较低，550 nm附近的绿波段反射率较高，白洋淀、官厅水库、白龟山水库绿波段反射峰出现了不同程度往长波方向移动的现象。590—670 nm 的黄红波段反射率下降较快，675和700 nm附近有明显的叶绿素反射谷和反射峰，其中，丹江口水库在该波段的峰谷特征不明显。750 nm 以后的近红外波段反射率普遍较低。因此，以上5个研究区是比较典型的内陆水体。

图4 实测水体光谱曲线图（光谱曲线来自85个采样点，每个研究区光谱为该区域所有采样点的平均值）Fig.4 In situ measured water spectrums of remote sensing reflectance（The spectrums come from 85 sampling points，and the spectrum of each study area is the average of all sampling points in the area）

表3 研究区实测水质参数统计表Table 3 Statistical of measured water quality parameters

本文统计了5 个研究区的叶绿素a 浓度和透明度的最小值、最大值、平均值和标准差（表3）。结果表明，白洋淀、官厅水库的叶绿素a 浓度较高，其次是小浪底和白龟山水库，丹江口水库叶绿素a浓度最低；丹江口水库透明度最高，其次是小浪底水库，白龟山水库、白洋淀、官厅水库透明度较低。叶绿素a浓度和透明度的统计结果和遥感反射率光谱特征具有较好的一致性。

4.2 Sentinel-2 MSI L2A地表反射率产品精度验证

Sentinel-2 MSI L2A 的地表反射率产品是基于均匀不变地物的HY-1C CZI 数据大气校正的基准，Sentinel-2 MSI L2A 的地表反射率产品的精度直接决定了HY-1C CZI数据大气校正的精度。图5是利用98 个实测数据对Sentinel-2 MSI L2A 的地表反射率产品精度检验的结果。结果表明，Sentinel-2 MSI L2A 的地表反射率产品和实测数据相比各波段都存在一定的高估现象，原因主要是气溶胶散射校正不足、缺乏水面天空光校正。但是高估具有很强的规律性，检验所用的7个区域不同，成像时间不同，气溶胶类型和光学厚度不同，但是表现出了非常高的相关性和一致性，绿、红波段的R2高达0.85 以上。蓝和近红外波段效果稍差，R2小于绿、红波段，但也超过了0.7。综上，Sentinel-2 MSI L2A 地表反射率产品尤其是可见光的蓝、绿、红波段）可以用于HY-1C CZI 数据大气校正的基准数据，但是需要对残余误差和系统性偏差做进一步校正。

图5 Sentinel-2 MSI L2A地表反射率产品精度（纵轴为地表反射率SR/π，约等于遥感反射率。B1、B3、B4、B8分别为Sentinel-2 MSI和HY-1C CZI对应的蓝、绿、红、近红外波段）Fig.5 The accuracy of Sentinel-2 MSI L2A surface reflectance product（where the vertical axis is the surface reflectance SR / π，which is about equal to the remote sensing reflectance.B1，B3，B4 and B8 are the blue，green，red and near-infrared bands of Sentinel-2 MSI corresponding to HY-1C CZI）

4.3 大气校正模型及精度验证

由于Sentinel-2 MSI 和HY-1C CZI 存在系统性偏差，经过基于均匀不变地物的相对大气校正后，需要进一步做系统定标。本研究建立了基于实测等效遥感反射率的HY-1C CZI 系统定标模型（图6）。结果表明：第1、2、3波段定标模型的R2较高，都超过了0.8。第4 波段数值较小，拟合趋势较差，部分点偏离较大，该波段需谨慎使用。

为了和FLAASH 大气校正结果做公平对比，本研究也建立了基于FLAASH 的大气校正结果的系统定标模型（图6）。结果表明：第1、3 波段定标模型的R2较低；第2、4波段定标模型的R2稍高；除第4 波段外，其他3 个波段定标模型的拟合趋势全面低于相对大气校正后定标模型。

图6 HY-1C CZI影像系统定标模型Fig.6 System calibration models of HY-1C CZI images

经过基于实测数据的系统定标以后获得了最终的HY-1C CZI 遥感反射率数据。采用同步实测的29 个独立采样点对相对大气校正得到的遥感反射率数据进行精度验证（图7）。结果表明：HY-1C CZI影像的第1、2、3波段误差都较小，数据点分布于1∶1 线附近；第4 波段误差较大，且数据点分布较散，数据的可靠性较低。综上，HY-1C CZI 遥感反射率数据的第1-3 波段精度较高，可以用于后续水质参数反演，第4波段精度较差，水质参数建模中应谨慎使用，尽量不用。

图7 HY-1C CZI影像大气校正精度评价Fig.7 The accuracies of HY-1C CZI image atmospheric correction

基于FLAASH 的大气校正精度如图5 所示。结果表明：第1、3、4 波段的AURE 和RMSE 都远大于相对大气校正的结果，特别是第1波段出现了较多负值，过校正现象严重。第2波段在数值较大区间表现较好，但是在数值较小的区间出现了数值低估现象，拉大了整体误差。

评价大气校正遥感反射率精度的另一个指标是光谱形状，用光谱相关系数和光谱角度距离表示。和实测光谱相比，全部29 个采样点的平均光谱相关系数r为0.98，平均光谱角度距离SAD 为0.11。表明遥感反射率光谱形状与实测光谱一致性较高（图8）。如图8（a）中光谱形状相似度极高，且每个波段的绝对数值也差异很小；图8（b）中光谱形状相似度很高，但是绝对数值上有一定差异，水质参数反演建模时采用比值和差值模型可以降低对反演结果的影响。与之形成对比，FLAASH 大气校正结果与实测光谱的平均相关系数r为0.68，平均光谱角度距离SAD 为0.56，校正精度远低于相对大气校正算法。

图8 大气校正后光谱形状对比Fig.8 Comparison of spectral shapes after atmospheric correction

4.4 水质参数反演建模及精度验证

利用实测水面遥感反射率，等效为HY-1C CZI波段，分别利用39 个和53 个实测等效数据建立水体叶绿素a和透明度反演模型。由于HY-1C CZI的第4 波段误差较大，因此，主要对比了前面3 个波段比值和差值模型。结果表明，B1/B2 和B1/B3 分别为叶绿素a 浓度和透明度反演的最优光谱指数。幂指数模型和线性模型分别应用于叶绿素a浓度和透明度反演可以得到最佳的拟合趋势（图9）。

图9 HY-1C CZI影像水体叶绿素a浓度和透明度反演模型Fig.9 Estimation models of chlorophyll-a concentration and secchi disk depth for HY-1C CZI images

本研究利用19 个独立采样点的实测数据评价了水体叶绿素a 浓度和透明度反演精度（图10）。结果表明，叶绿素a 浓度反演的AURE 为33.8%，RMSE 为4.8 μg/L；透明度反演的AURE 为25.0%，RMSE 为34.9 cm。散点图上，数据点总体在1∶1线附近，反演趋势合理。在光学复杂的内陆水体，对于波段较宽的4谱段多光谱影像来说，该精度达到了较好水平，具有应用于内陆水体水质遥感监测的潜力。

图10 HY-1C CZI影像水体叶绿素a浓度和透明度反演精度评价Fig.10 Estimation accuracy evaluation of chlorophyll-a concentration and secchi disk depth for HY-1C CZI images

4.5 研究区水质参数遥感反演结果

应用以上建立的HY-1C CZI 影像水体叶绿素a浓度反演模型，反演了官厅水库、小浪底水库、白洋淀、白龟山水库、丹江口水库的叶绿素a 浓度，生成了5 个研究区的叶绿素a 浓度分布图（图11）。结果表明：官厅水库、白洋淀和小浪底水库下游的西霞院水库的叶绿素a浓度较高，大部分区域的叶绿素a 浓度在20 μg/L 左右；小浪底水库（主库区）、白龟山水库的叶绿素a 浓度较低，大部分区域在10 μg/L左右；丹江口水库的叶绿素a浓度最低，中心库区在5 μg/L 左右，水库边缘及水库东部沿岸有少量区域叶绿素a浓度稍高，但大部分区域也低于10 μg/L。

图11 基于HY-1C CZI影像反演的典型内陆水体叶绿素a浓度分布图Fig.11 Distribution of chlorophyll-a concentration in several typical inland waters based on HY-1C CZI image

5 个研究区的透明度反演结果如图12 所示。结果表明：白洋淀和官厅水库水体透明度最低，大部分区域透明度在100 cm 左右；白龟山水库透明度空间分布差异较大，北部和东部沿岸透明度较低，大概在150 cm 左右，中西部透明度较高，达到200 cm，部分区域超过250 cm。小浪底水库透明度较高，且空间差异很小，全库区基本都在250 cm 左右。丹江口水库透明度最高，中心库区超过350 cm，库区边缘也超过了300 cm。

图12 基于HY-1C CZI影像反演的典型内陆水体透明度分布图Fig.12 Distribution of secchi disk depth in several typical inland waters based on HY-1C CZI image

5 结 论

本研究在华北平原的不同浑浊程度的小浪底水库、官厅水库、丹江口水库、白龟山水库、白洋淀开展了5 次星地同步实验，获取了85 个采样点的水面光谱、水体叶绿素a浓度、透明度及其它辅助数据。根据实测数据开展了HY-1C CZI 影像在内水体大气校正、水体叶绿素a浓度、透明度反演的建模及评价，验证了HY-1C CZI 影像在内陆水体水色遥感研究的应用潜力。

针对HY-1C CZI 影像缺少短波红外大气校正波段，精确大气校正较为困难的问题，发展了基于均匀不变地物的相对大气校正算法。该算法利用最近的两景Sentinel-2 MSI L2A 地表反射率产品提取均匀不变地物，建立HY-1C CZI 影像上均匀不变地物像元的辐射值与相应的Sentinel-2 MSI L2A地表反射率之间的关系，即为相对大气校正模型。最后根据星地同步实测光谱进行系统定标，校正HY-1C CZI 和Sentinel-2 MSI 波段设置差异等导致的系统误差。HY-1C CZI第1-4波段相对大气校正的无偏相对误差分别为：14.7%、11.2%、28.9%、41.7%，说明前3 个波段的精度相对较高，第4 波段精度较低，第4 波段不适宜参与后续水质参数反演。相对大气校正后光谱和实测等效光谱的平均相关系数为0.98，平均光谱角度距离为0.11，结果表明：相对大气校正后HY-1C CZI 影像的光谱形状和实测光谱一致性较高。

应用基于均匀不变地物大气校正后的遥感反射率影像建立了水体叶绿素a浓度和透明度遥感反演模型，通过比较认为B1/B2 和B1/B3 分别为叶绿素a 浓度和透明度反演的最优光谱指数，据此建立的叶绿素a浓度和透明度反演模型的平均无偏相对误差和均方根误差分别为：33.8%和25.0%、4.8 μg/L和34.9 cm。将反演模型应用于HY-1C CZI影像上，生产了5 个研究区的叶绿素a 浓度和透明度分布图。结果表明，叶绿素a浓度和透明度反演结果和实测数据一致，反演结果较为可靠。

随着HY-1D 卫星的发射，和HY-1C 形成上下午组网观测，可以实现更高的观测频率，对内陆水体水色遥感研究具有重要意义。但是，其应用效果仍需要在更多类型的研究区进行检验和应用。