王道生，杜克平，陈树果，薛程，叶小敏，李忠平
1.厦门大学 海洋与地球学院, 厦门 361102;
2.北京师范大学 地理科学学部, 北京 100875;
3.中国海洋大学 海洋技术学院, 青岛 266100;
4.中国海洋大学 三亚海洋研究院, 三亚 572024;
5.自然资源部国家卫星海洋应用中心, 北京 100081;
6.麻省州立大学波士顿分校 环境学院, 美国波士顿 02125
自1970年代美国发射CZCS（Coastal Zone Color Scanner Experiment）水色卫星，国际上发射运行了多颗关注全球海洋生态环境的水色遥感卫星，特别是20世纪90年代末发射、运行先进的SeaWiFS（Sea-viewing Wide Field-of-view Sensor，生命周期：1997年—2010年）、MODIS（Moderate-resolution Imaging Spectroradiometer，生命周期：2002年至今）等水色卫星，在轨运行期间积累大量丰富的海洋水色卫星观测数据。同时，NASA 由1993年开始组织开发相应的水色卫星数据处理系统SeaWiFS Data Analysis System（SeaDAS），最终成为一个集处理、显示、分析和质量控制为一体的水色卫星数据综合处理系统（Fu等（1998））。起初，系统仅支持SeaWiFS 的数据处理和产品制作，包括从L0（经过数据解包、元数据校核，未做任何纠正的相机载荷原始数据及辅助信息）到L3（统一时空网格尺度的全球地图投影后图像产品）所有等级数据进行显示、处理、分析和质量控制。经过几十年的不断改进和扩展，除了保持对SeaWiFS的支持，现已扩展到CZCS、MODIS、VIIRS（Visible and Infrared Imagere/Radiometer Suite）、 MERIS（Medium Resolution Imaging Spectrometer）、GOCI（Geostationary Ocean Color Imager）、Landsat8、Sentinel2/3、HICO（Hyperspectral Imager for the Coastal Ocean）、HawkEye 等卫星传感器。目前最新版本为SeaDAS8.1.0，于2021年6月24日发布（https：//seadas.gsfc.nasa.gov ［2021-10-04］），其中的Ocean Color Science SoftWare package（OCSSW）科学数据处理包需要在Linux 环境下运行。它通过组件形式，可实现对自选卫星传感器数据的处理，包括数据辐射定标、几何定标、大气校正、产品反演和地图投影融合等，并且开放源代码，支持用户添加自定义卫星传感器及要素产品。国内外学者Wang 等（2014）、Franz 等（2015）、Wang 等（2016）、Chen 等（2021）也已在这方面做过尝试，实现了新版SeaDAS 对自选卫星传感器的新算法植入和遥感产品制作。比如Frouin 等基于NASA 的6.4 版本SeaDAS（不支持VIIRS 传感器的处理），开发NOAA-SeaDAS 版本，对VIIRS 等卫星数据提供大气校正功能，实现了对当时新型卫星传感器的大气校正处理。Bryan 等针对Landsat 8 的OLI（Operational Land Imager）传感器，同样扩展了当时SeaDAS 版本对其的支持，并植入多种大气校正算法可供程序运行处理时选择。Wang等在NOAASeaDAS 版本基础上，进一步拓展其功能，利用近红外波段和短波红外波段对VIIRS传感器进行在轨替代定标工作。Chen等利用SeaDAS-C版本实现对HY-1C的COCTS（Chinese Ocean Color and Temperature Scanner）传感器卫星数据的处理，并评估COCTS传感器在全球水色遥感反演中的性能。
根据中国民用空间基础设施规划，正全力发展海洋水色、海洋动力环境和海洋监视监测3大系列的海洋卫星（蒋兴伟 等，2016）。其中，海洋水色观测卫星主要通过获取、分解海洋的水体光谱，反演得到海洋水体的固有光学量以及叶绿素、悬浮颗粒物和有色溶解有机物浓度等水质信息。这些产品可应用于海洋资源调查、水质监测、赤潮监测和海岸带植被监测等方面。在2002、2007、2018 和2020年，国家卫星海洋应用中心已先后发射、运行HY-1A/1B/1C/1D 水色卫星。特别是，HY1C/1D 实现双星上下午组网运行，不但提高探测效率，也有助于分析浮游植物的日内动态变化。星上搭载海洋水色水温仪、海岸带成像仪、紫外成像仪、星上定标光谱仪和一套船舶监测系统（AIS）。其中海洋水色水温仪（COCTS）覆盖周期单颗为1 d，双星组网为0.5 d，星下点地面像元小于1100 m，幅宽大于2900 km，可实现对西北太平洋区域的实时观测和西北太平洋区域之外的全球其他区域的非实时观测，因此能够对广袤的海洋进行高效率的测量（http：//www.nsoas.org.cn/［2021-11-04］）。中国海洋水色卫星星座的发展历经从HY-1A/1B 试验型卫星到HY-1C/1D 业务型卫星转变，未来还将继续发射新一代海洋水色HY1E卫星和正在规划的静止轨道卫星（林明森 等，2019），进一步提升对海洋生态环境的观测能力。
随着海洋卫星观测星座的完备，需与之建立相匹配的通用卫星数据地面处理系统，避免出现一颗卫星一套系统的情况。由于中国海洋卫星发展起步较晚，科学研究长期依赖国外卫星数据及其处理系统。因此，虽然有了水色卫星及其标准产品，但过去没有一套针对HY-1C/D 的离线数据处理系统，用户则不能够根据自己的需求在HY-1C/D 的标准产品之外产生自定义水色产品，妨碍HY-1C/D 的广泛应用。另外，前人在国内外水色卫星的大气校正和水色要素反演方面积累了大量算法，若能在类似SeaDAS 平台上实现各自算法的集成，将算法代码共享发布，有助于算法的进一步推广应用及适应性改进。因此基于国内外学者的工作基础，一个行之有效的方式是借鉴SeaDAS科学处理包软件框架，根据HY-1C/D 的载荷特征和数据格式特点，针对性开发组件实现从L1B 级别数据到L2 产品及后续L3 级产品的制作。本文则描述以HY-1C/D 的水色水温仪L1B 数据为数据源，开发的适用于国产水色卫星大气校正及水色要素反演、卫星原始分辨率下的离线处理系统OffLine-COCPS（OffLine-Chinese Ocean Color satellite Processing System）。旨在通过本文研究工作，阐明对多颗水色卫星数据集成处理的系统工作方式，该系统不仅保持对各种卫星数据格式的识别支持，同时也支持对前人已有算法产品进行集成，通过一套系统，用户即可对国产水色卫星数据处理实现自主扩展应用。
2.1 数据源
国家卫星海洋应用中心已经通过中国海洋卫星数据服务系统将HY-1C/D 卫星数据进行全面的共享（https：//osdds.nsoas.org.cn［2021-11-04］），通过该网站可获取到经过几何定位和辐射定标后的COCTS L1B 级别卫星数据。OffLine-COCPS 将以此数据作为输入数据，开发多项数据模块，从而自主实现HY1C/1D 数据的大气校正和产品要素反演。
L1B 数据是以HDF5 格式进行存储，内容包含Global Attributes 全局属性、QC Attributes QC 属性、Calibration 定标属性、Scan_Line Attributes 扫描行属性、Geophysical Data 地球物理数据和Navigation Data图像定位数据等。
本文采用2019年9月8日HY-1C白天共119景L1B数据及2020年12月8日HY-1D白天共140景L1B数据作为数据源，部分数据的覆盖情况如图1所示。
图1 HY-1C/D卫星COCTS数据全球覆盖示意图Fig.1 Global coverage of HY-1C/D COCTS data
验证数据采用NASA 官方提供的MODIS 和VIIRS 对应9 km 叶绿素产品和遥感反射比产品（https：//oceancolor.gsfc.nasa.gov/l3/［2022-03-18］），对HY1C/1D COCTS 传感器反演得到的叶绿素产品和遥感反射比产品进行评价。
2.2 COCPS系统架构
为实现多个传感器的自由扩展，OCSSW 处理包采用如图2的接口设计方式，将传感器观测到的辐亮度值和其他辅助文件、查找表、实测值，以及公用的辅助文件和每个传感器独有的配置文件分离处置，因此，在添加新传感器时，只需按照相应格式配置专用传感器的配置文件。
图2 OCSSW传感器独立设计方式Fig.2 OCSSW sensor-independent design approach
考虑到HY-1C/D COCTS L1B 数据格式均以HDF5 方式存储，仅在全局属性Global Attributes 中以title或Satellite Name字段内容存在区别，在两个传感器集成方面，以title属性作为关键字段，分别对HY-1C 和HY-1D 传感器数据识别分类处理，以实现在一套系统上对多传感器的支持。
2.3 COCPS系统开发
根据NASA 官方OBPG（Ocean Biology Processing Group）组发布的说明文档，为了使对OCSSW科学处理包的修改生效，需在对应的开发环境中对OCSSW 处理包进行重新编译。以SeaDAS-7.3.2 版本为例，包含组件l2gen、l2bin、l3bin、l2mapgen、l3mapgen、smigen 等，第三方插件hdf5、levmar、jpeg、lapack等。
本文以CentOS7 系统为例，详细说明系统搭建所需环境配置，对于OCSSW 处理包则需要Python3.6 及以上版本。OffLine-COCPS 编译所需的步骤可参考https：//oceancolor.gsfc.nasa.gov/docs/ocssw［2021-10-10］说明，根据使用到的组件，这里主要以l2gen 为例，其功能为实现L1B 数据到L2 数据的处理，分别在头文件中增加HY-1C 传感器的定义以及在原有源文件中新增及修改部分代码，HY-1D 的添加方式同HY-1C；
最后编译生成可执行程序供命令行或第三方调用，该程序既可分别对应识别HY-1C/D 的COCTS 传感器，也可支持原有SeaWiFS 传感器，可供验证程序修改正确性。该模式除了支持增加新的传感器，同样也支持新的算法植入及生产新的水色要素产品，其方法同上。
2.4 COCPS处理流程
OffLine-COCPS 将以上L1B 数据和相应的辅助配置文件作为输入，包括考虑COCTS 波段特性，采用基于逐次散射法开发的矢量海气耦合的辐射传输模型自主生成的HY-1C/1D 专用的瑞利、气溶胶、大气漫射透过率查找表，调用Gordon 和Wang（1994）以及Bailey 等（2010）的大气校正算法对可见光波段进行大气校正，得出海面遥感反射率（Rrs）；
其中气溶胶查找表的生成方式基于Ahmad等（2010）。进而根据水色要素反演算法得出期望的水色要素产品。
3.1 大气校正
通常将从水色遥感传感器接收到的总辐射中扣除大气的贡献和影响，得到携带水体性质的海面向上的出射辐亮度这一过程称为大气校正。典型情况下，水色遥感器接收总辐射的80%左右来自大气层散射和海洋表面的反射，而真正携带海洋水色信息的离水辐射只占其中很小的一部分。大气校正算法的过程就是准确地剔除大气影响，进而根据光辐射与水体中物质含量之间的经验或理论关系，推导出水体中物质含量的信息。
处理过程中，为了剔除太阳辐射变化的贡献，通常以天顶的反射率为输入开始计算。该反射率（ρ）通常定义为
式中，Lt为进入卫星传感器的辐亮度，F0为大气顶端的太阳辐照度，θ0为太阳天顶角。在此基础上，传感器在某一波长λ 接收的ρ可表达为（Gordon 和Wang（1994））
式中，ρr代表大气分子瑞利散射反射率，ρa代表气溶胶散射反射率，ρcoupl代表太阳耀斑、分子散射和气溶胶散射之间多次散射的反射率，ρg代表太阳直射反射率。ρw代表离水辐射的反射率，toz为臭氧漫透射率，T为大气直射透射率，t为目标到传感器的透射率，t0为太阳到目标的透射率。其中ρw（=πRrs，Rrs表示水面之上离水辐亮度与水面之上下行辐照度的比值）则是需要从卫星测量中获取的有关水体信息的表观光学量。以下分别对各个分量的计算进行简单介绍：
大气分子瑞利散射反射率：首先，利用矢量辐射传输模型提前生成COCTS 的专用瑞利查找表，其目的是在保证计算精度的前提下，尽可能提高计算效率。查找表计算中假设大气为一层结构（海洋部分无散射），由纯大气分子组成，其光学厚度为传感器各波段等效瑞利散射光学厚度（Bodhaine 等，1999），单次散射反照率为1，退偏系数跟波段相关（Bates，1984）。查找表中太阳天顶角变化范围为0°—88°，变化步长为2°（45 个角度）；
传感器天顶角为勒让德多项式在（0，1）区间上的零点，变化范围大致为0°—84.3°，变化步长大约为2°（41个角度）；
风速变化范围为0—25 m/s，变化步长为2.5 m/s（11 个），考虑阴影效应（Cox和Munk，1954；
Mishchenko 和Travis，1997）。查找表中存储了Stokes矢量对相对方位角傅立叶展开的3 个系数。方位在COCTS 输入的观测几何和海表风速条件下，利用COCTS 瑞利散射查找表插值出给定太阳方位角、观测方位角、相对方位角和海表风速下的瑞利反射率系数；
其次，根据地球表面的大气压P对瑞利光学厚度进行校正；
最后利用瑞利散射系数和瑞利光学厚度计算出大气分子瑞利散射反射率数值。
气溶胶散射反射率：采用Ahmad 等（2010）分析AERONET 中3 个浑浊近岸站点和8 个清澈海洋站点的多年长期现场观测数据，模拟出的包括8种相对湿度和10种细模态体积浓度分数下共80组非吸收性和弱吸收性气溶胶模型，利用Mie散射模型（Mishchenko 等，2002）计算气溶胶单次散射反照率及散射相矩阵，进一步通过矢量辐射传输模拟构建包括气溶胶单次散射反射率和气溶胶反射率单多次变换系数的查找表；
查找表计算中假设大气为两层结构（海洋部分无散射），上层为大气分子（计算方式同瑞利散射查找表），下层为气溶胶，假定为平静水面（风速效应在瑞利散射查找表中考虑）。查找表中太阳天顶角变化范围为0°—80°，变化步长为2.5°（33个角度）；
传感器天顶角为勒让德多项式在（0，1）区间上的零点，变化范围大致为0°—76°，变化步长大约为2°（35 个角度）；
相对方位角为0°—180°，变化步长为10°（19 个角度）；
气溶胶光学厚度分别为0.01、0.05、0.1、0.15、0.2、0.25、0.3、0.4、0.5和0.7（10种）。查找表中存储不同气溶胶模式、太阳天顶角、传感器天顶角及相对方位角条件下，气溶胶多次散射反射率（包含气溶胶与分子相互作用部分）与单次散射反射率拟合二次多项式的3 个系数（Gordon 和Wang，1994）。基于Gordon 和Wang（1994）提出近红外波段离水辐射率为零的暗像元假设理论，在考虑相对湿度和气溶胶模型权重影响下，应用于气溶胶校正（Bailey等，2010）。
太阳直射反射率：根据Cox 和Munk（1954）的计算方法，在给定太阳方位角、卫星观测角和相对方位角，海面风速的情况下，在有轻微太阳耀斑的影响下，可计算出太阳直射反射率。
大气漫射透过率：利用矢量辐射传输模型计算不同天顶角、水面风速、以及气溶胶模式、光学厚度组合情况下大气漫射透过率查找表。查找表计算中假设大气为两层结构，上层为大气分子（计算方式同瑞利散射查找表），下层为气溶胶，假定为平静水面（风速效应在瑞利散射查找表中考虑），海洋部分为纯吸收一层结构，假定了一个单位强度各向同性上行辐射（Yang 和Gordon，1997）。查找表中天顶角变化范围为0°—80°，变化步长为2.5°（33 个角度）；
气溶胶模式及光学厚度与气溶胶查找表一致。查找表中存储不同条件下的大气漫射透过率。
考虑到暗像元假设成立的条件及不同传感器波段配置，系统中预留了不同的气溶胶反射率计算模式可供选择，用户可根据传感器及感兴趣反演区域选择性地使用对应计算模式，具体可见Wang 和Shi（2007），Vanhellemont 和Ruddick（2015）。
3.2 叶绿素与透明度反演
OCSSW 科学处理包中l2gen组件集成了多种水色要素产品及多种反演算法，用户可根据需求采用备用算法，可在配置文件中指定或程序调用时在参数中指定。若没有指定，水色要素产品大部分算法采用l2gen 组件中默认自带的算法，如叶绿素-a（Chlora），可根据需要调用OCX（O’Reilly 等（1998））、OCI（Hu 等（2012））或其他已经提出并开发完成的模块。OCI算法是将经验波段比值OCX 算法和基于波段差值的CI 算法相结合，其中CI 算法主要应用于开阔大洋的清洁水体中，OCX算法主要应用于近岸和内陆水体，二者过渡区域采用CI 和OCX 相结合的线性插值方法。系统根据COCTS 的波段设置自动将Rrs插值到反演算法所需的波段上，其中的反演系数用户可根据需要在配置文件中进行调整。此外，如上所述，系统可增加自定义反演算法。比如，本文根据Lee 等（2015）提出的新的透明度算法（Zsd），将透明度和水体漫衰减系数相关联，将其算法植入重新编译后的l2gen中。该算法的核心在于，基于辐射传输和人眼目标探测原理，简化后的Zsd可表达为
式中，Kdtr为漫衰减系数Kd在可见光波段（400—700 nm）之间的最小值，为对应的遥感反射比。Kd可由Lee 等（2005）提出的半分析算法从Rrs 遥感反射比中反演得出。
首先，针对遥感反射比产品做出评价，将NASA 官方提供的MODIS 卫星9 km 各波段Rrs产品作为参考产品，与VIIRS对应波段进行比较，验证作为参考产品的准确度，比较结果如图3 所示。可以看出MODIS 载荷生成的产品保持较好的稳定性，各波段的散点图均匀地分布在1∶1 线两侧，相关性均达到0.9以上，考虑到MODIS波段配置与HY-1C/1D相近，本文选取MODIS产品作为参考产品，对反演生成的HY-1C/D产品进行评价。
图3 2020年12月8日VIIRS与MODIS反演全球各波段遥感反射比比较散点频数图Fig.3 Comparison between VIIRS and MODIS global ocean Rrs for December 8，2020
利用以上系统，参照SeaWiFS、MODIS等卫星的L2 级别数据，输出HY-1C/D 的L2 级nc 文件，大气校正后的遥感反射比产品及水色要素产品均存储在Geophysical Data 属性组中。统计计算得到的2020年12月8日HY-1D的遥感反射比（Rrs）产品和同期MODIS 卫星9 km 的各波段Rrs 产品进行比较，结果如以下散点图所示。412 nm 波段和443 nm 波段的Rrs 高密度线性分布在1∶1 线上，表明整个波段的空间分布趋势与MODIS获得的分布基本一致，进一步说明利用自主生成的HY-1C/D大气校正相关查找表基本满足了在全球范围内的应用，而490 nm，520 nm，565 nm，670 nm 波段散点对比结果略低于412 nm 波段和443 nm 波段，这一方面可能与两个传感器在这些波段设置的偏差有关，另一方面在红光波段Rrs值本身就低，轻微的误差就容易引起整个散点图偏离1∶1 线。这与Chen 等（2021）得到的结果一致，在490 nm，520 nm，565 nm，670 nm 波段上整体的离散度和不确定度要高于412 nm波段和443 nm波段。
图4 2020年12月8日HY-1D卫星COCTS与MODIS反演全球各波段遥感反射比比较散点频数图Fig.4 Comparison between HY-1D COCTS and MODIS global ocean Rrs for December 8，2020
其次，利用叶绿素和透明度反演算法，绘制2019年9月8日1 点55 分（UTC 时间）HY1C 拍摄到的澳大利亚西北沿岸的叶绿素和水体透明度产品。整体分布情况如图5、图6 所示，整体分布趋势合理，表明调用系统原有叶绿素算法抑或是重新植入的透明度算法，均可满足生产需求。
图5 HY-1C卫星COCTS 2级产品叶绿素-a分布图Fig.5 HY-1C COCTS level-2 chlorophyll-a product map
图6 HY-1C卫星COCTS水体透明度产品分布图Fig.6 HY-1C COCTS Secchi disk depth product map
基于2 级产品质量标识，进一步对2 级产品进行统计和投影分析，得到9 km 分辨率的单日全球分布图；
图7、图8 分别展示HY-1C 载荷COCTS的2019年9月8日和HY1D 载荷COCTS 的2020年12月8日的单天叶绿素产品分布图。
图7 2019年9月8日HY-1C卫星COCTS反演全球海洋叶绿素-a产品图Fig.7 HY-1C COCTS global chlorophyll-a product map for September 8，2019
图8 2020年12月8日HY-1D卫星COCTS反演全球海洋叶绿素a产品图Fig.8 HY-1D COCTS global chlorophyll-a product map for December 8，2020
将HY-1D与同一日期的MODIS对应9 km叶绿素-a 产品对比，可得如下结果，采用最小二乘法得到拟合方程式，可以看出散点整体分布在1∶1左右，颜色越深表示该处散点越密集，且高密度散点线性分布在1∶1 线上，说明反演得到的HY-1D 叶绿素结果同MODIS的整体空间分布趋势基本一致。但HY-1D 的结果略偏高于MODIS 结果，这可能与大气校正后得到的Rrs 遥感反射比存在偏差有关，也与叶绿素反演算法中的经验系数有关。本文主要阐述系统建立过程的关键技术，旨在提升国产海洋水色卫星数据的实际应用，大气校正后得到的遥感反射比和反演产品精度的提升，不在本文研究范围内，后续研究中可进一步对数据产品的精度进行深入探究。
图9 2020年12月8日HY1D卫星COCTS与MODIS全球海洋叶绿素-a比较散点图Fig.9 Comparison between HY1D COCTS and MODIS global ocean chlorophyll-a product for December 8，2020
基于SeaDAS 的科学处理包框架OCSSW 开发的OffLine-COCPS，不仅保留对原有传感器的处理功能，同时增加了对国产水色卫星HY1C/1D 载荷COCTS 的支持，此外支持水色反演要素产品（如透明度）的定制开发，实现了多星传感器在一套系统下的集成，并兼容原有多种大气校正方法及水色要素反演模块，保持程序使用的灵活性。基于此框架，研究人员可扩展到对应的国产水色卫星处理系统，不仅限于HY-1C/1D 的COCTS 水色水温仪传感器，为后续研究提供便利。通过此平台，可以实现算法和传感器的集成；
同时，开源代码也可使用户在使用过程中更方便提出建议，本研究仅对水色要素反演算法进行扩展，后续亦可根据需要对大气校正算法进行开发，如加入光谱优化算法或神经网络算法等。关于产品的精度分析，涉及到辐射定标、大气校正、算法反演等诸多方面，后续可针对具体问题，对产品精度做进一步分析与提升。
志 谢感谢国家卫星海洋应用中心提供开放HY-1C/1D 数据，SeaDAS 数据处理框架来自https：//seadas.gsfc.nasa.gov/installers/［2021-11-04］。
猜你喜欢水色气溶胶反射率影响Mini LED板油墨层反射率的因素印制电路信息(2022年11期)2022-11-30近岸水体异源遥感反射率产品的融合方法研究海洋通报(2022年4期)2022-10-10具有颜色恒常性的光谱反射率重建光谱学与光谱分析(2022年4期)2022-04-06水色小资CHIC！ELEGANCE(2021年28期)2021-08-18雨花·艺术 徐华翎作品雨花(2020年9期)2020-09-10CF-901型放射性气溶胶取样泵计算公式修正辐射防护通讯(2019年3期)2019-04-26气溶胶中210Po测定的不确定度评定四川环境(2019年6期)2019-03-04基于地面边缘反射率网格地图的自动驾驶车辆定位技术汽车文摘(2018年2期)2018-11-27Aerosol absorption optical depth of fine-mode mineral dust in eastern ChinaAtmospheric and Oceanic Science Letters(2016年1期)2016-11-23四川盆地秋季气溶胶与云的相关分析高原山地气象研究(2016年2期)2016-11-10