HY-1D COCTS 海表温度反演与印证

刘铭坤， 管磊， 刘凡莉， 刘建强

1. 中国海洋大学 信息科学与工程学部海洋技术学院, 青岛 266100;2. 中国海洋大学 三亚海洋研究院, 三亚 572024;3. 青岛海洋科学与技术试点国家实验室 区域海洋动力学与数值模拟功能实验室, 青岛 266071;4. 自然资源部 国家卫星海洋应用中心, 北京 100081;5. 自然资源部 空间海洋遥感与应用重点实验室, 北京 100081

1 引 言

海表温度SST（Sea Surface Temperature）是地球气候系统中一个非常重要的变量。 由于海表层位于海洋和大气的界面，SST对两者以及两者之间的热量、水分、动量和气体的交换都至关重要（Bentamy 等，2017；Wanninkhof 等，2009）。由于自然界的几乎所有动力过程都表现出对温度的依赖性，SST 的变化直接对气候变化造成影响。因此，SST广泛应用于气候和海洋模式的评估、气候变化的观测量化、过程理解和参数化、海洋生态学、海洋学和地球物理学等。

SST 的现场测量已经有150 多年的历史了，最初主要是船测，最近几十年补充了漂流浮标和固定浮标等测量方式（Kent等，2017）。卫星观测SST产品可提供更可靠和更完整的时空采样，主要观测波段包括红外和微波。红外SST观测相比微波观测具有更高的空间分辨率，通常使用中红外和热红外波段来测定其辐射，但这两个波段都不能穿透云层，并且受气溶胶和水蒸气的影响，所以红外SST的观测需要进行大气校正并且只能在无云区使用（Minnett等，2019）。自1981年搭载在美国国家海洋和大气管理局NOAA（National Oceanic and Atmospheric Administration）-7 卫星上的甚高分辨率辐射计AVHRR（Advanced Very High Resolution Radiometer）二代业务化运行以来，随着卫星传感器设计水平的提高和SST反演算法的改进，红外卫星业务化SST 观测的精度也得到了进一步的提高。目前国外主流的极轨卫星观测的红外SST业务化产品主要包括NOAA 发布的AVHRR Pathfinder SST 产品，美国宇航局NASA（National Aeronautics and Space Administration）发布的中分辨率成像光谱仪MODIS （Moderate-Resolution Imaging Spectroradiometer）SST 产品，NOAA 发布的可见光红外成像辐射计VIIRS（Visible Infrared Imaging Radiometer Suite）SST产品，欧洲航天局ESA（European Space Agency）气候变化专项CCI（Climate Change Initiative）发布的沿轨扫描辐射计ATSRs（Along Track Scanning Radiometers）SST 产品和欧洲气象卫星开发组织EUMETSAT（European Organisation for the Exploitation of Meteorological Satellites）开发的海洋陆地表面温度辐射计SLSTR（Sea and Land Surface Temperature Radiometer）SST 产品。其中AVHRR 和MODIS 的SST 产品精度在0.5—0.6 K 左右（Saha 等，2020；Kilpatrick 等，2015）；VIIRS 与现场实测SST 的比较结果表明二者的平均偏差为0.0 K（白天和夜间），标准偏差分别为0.466 K（白天）和0.359 K（夜间）（Petrenko 等，2014）；ATSR 系列SST 产品精 度 在0.2—0.3 K 左 右（Merchant 等，2019）；SLSTR SST 产品精度在0.3 K 左右（Luo等，2020）。中国也有多颗卫星可以用于红外SST观测，主要包括风云（Fengyun，FY） 系列和海洋（Haiyang，HY）系列卫星。FY-3C 卫星上搭载的可见光红外扫描辐射计VIRR（Visible and Infrared Radiometer）全球SST产品的精度在0.6 K左右（Li等，2021）。

海洋一号（Haiyang-1，HY-1）卫星是中国第一代应用于海洋遥感的海洋卫星。海洋一号A（Haiyang-1A，HY-1A）卫星作为试验业务化卫星，于2002 年5 月发射成功。随后，海洋一号B（Haiyang-1B，HY-1B）卫星在HY-1A卫星的基础上研制，于2007 年4 月成功发射，其观测能力进一步增强和提高。海洋一号C（Haiyang-1C，HY-1C）和海洋一号D（Haiyang-1D，HY-1D）卫星分别于2018年9月和2020年6月成功发射，是中国HY-1 系列卫星在轨业务化运行的海洋遥感卫星。HY-1C/D 卫星采用上、下午卫星组网，提高了全球覆盖能力。HY-1D 卫星的轨道类型为太阳准同步近圆形极地轨道，轨道高度为782 km。卫星降交点地方时为下午13：30±30 min。星上搭载有10波段的海洋水色扫描仪COCTS（Chinese Ocean Color and Temperature Scanner）、4波段的海岸带成像仪、紫外成像仪和星上定标光谱仪。HY-1D 卫星的主要探测要素包括海水光学特性、叶绿素浓度、悬浮泥沙、可溶有机物及SST等。COCTS有8个可见光近红外通道（0.402—0.422 μm，0.433—0.453 μm，0.480—0.500 μm，0.510—0.530 μm，0.555—0.575 μm，0.660—0.680 μm， 0.730—0.770 μm， 和0.845—0.885 μm）用于海洋水色观测和2 个热红带通道（10.30—11.40 μm和11.40—12.50 μm）用于SST观测（http：//www.nsoas.org.cn/news/content/2018-11/23/44_5226.html［2021-11-02］）。星下点空间分辨率为1.1 km，刈幅宽度为2900 km。COCTS 的扫描方式为推帚式扫描，在垂直于扫描方向上共有8个探测器，每行扫描像素数为1656。Liu 等（2018）和Liu 等（2021）分别对HY-1B 和HY-1C COCTS 的两个热红外通道进行了交叉定标，交叉定标后，HY-1B 和HY-1C COCTS 两热红外通道的亮温精度在0.3 K左右（Liu等，2018，2021）。

本文基于HY-1D COCTS 的热红外通道观测亮温，开展云检测和SST反演。首先对HY-1C COCTS进行条纹噪声去除，接下来设计了贝叶斯云检测算法对COCTS 进行云检测，最后实现最优估计SST 反演，并采用现场实测和VIIRS SST 产品对反演得到的COCTS SST进行印证。

2 数据集

2.1 HY-1D COCTS L1A级数据

HY-1D COCTS L1A 级数据由国家卫星海洋应用中心（https：//oHDFsdds.nsoas.org.cn/［2021-11-02］）处理分发，以分层数据格式（Hierarchical Data Format）存储。

HY-1D COCTS L1A 级数据由全局属性信息和数据变量信息组成。全局属性信息包括传感器名称、总波段数、波段参数、白天/夜间标志信息、扫描起始时间、扫描终止时间、扫描4个角的经纬度信息、扫描行数、每行扫描数、扫描中心点地球位置信息等。数据变量信息分为传感器校准信息、辅助数据信息、地球物理信息、地理定位信息、质量控制信息和扫描属性信息等。其中地理定位信息包括定位点数目、首个定位点位置、定位点之间间隔、经度、纬度、太阳天顶角、太阳方位角、卫星天顶角、卫星方位角等。扫描属性信息包括起始扫描行每个点的经纬度、终止扫描行每个点的经纬度、每个扫描行中心点的经纬度、每个扫描行中心点的太阳天顶角、每个扫描行的探测器编号、红外通道每个扫描行的校正系数和偏差值、每行扫描时间等。地球物理数据中包括412 nm、443 nm、490 nm、520 nm、565 nm、670 nm、750 nm、865 nm、11 μm 和12 μm 波段大气层顶DN值。

2.2 NWP ERA5数据

本文使用NWP ERA5数据（Hersbach等，2020）作为大气辐射传输模型的输入数据。NWP ERA5数据由欧洲中期天气预报中心ECMWF（European Centre For Medium Range Weather Forecasts） 处 理分发，包括大气参数和表面参数，每隔一个小时提供数据，被存储为NetCDF4格式。

ERA5 数据被存储在0.25°×0.25°经纬度网格中，包括大气廓线数据和表面地球物理参数。ERA5数据提供了37层大气的廓线参数，包括大气气压、温度、相对湿度、臭氧质量混合比、云液态水量和云覆盖量等。表面参数包括海陆标识、总云量、表面压力、表层温度、SST、水汽总量、10 m 风速水平分量、10 m 风速垂直分量、2 m 温度、2 m露点温度等。

2.3 iQuam数据

本文选择由NOAA 开发的iQuam 数据用于SST的评估和印证（Xu 和Ignatov，2014）。iQuam 数据存储为NetCDF4 格式，每月存储为一个数据文件，每24 h 更新一次。数据文件分为全局属性和数据变量两个部分。全局属性包括iQuam数据的生成机构、数据版本、处理级别、起始时间、终止时间、地球位置范围等。数据变量包括年、月、日、小时、分钟、秒等时间信息，经度、纬度等地球位置信息，大气温度、大气压力、风向、风速、云覆盖、SST等观测信息，以及质量控制标识信息和观测平台类型标识信息等。质量控制信息提供了SST 数据的整体质量级别，分为0—5 共6 个级别，5表示质量最高级别。本文选用最高质量级别以及观测类型为船、漂流浮标、热带系泊浮标、海岸系泊浮标和高分辨率漂流浮标的数据。

2.4 Suomi NPP VIIRS L2级数据

本文选择NOAA 提供的Suomi NPP VIIRS L2 级SST 数据用于与COCTS SST 的对比（ftp.nodc.noaa.gov［2021-11-02］）。VIIRS L2级数据是未经几何校正的卫星轨道SST数据，存储为NetCDF4格式。数据分为全局属性和数据变量两部分。全局属性包括数据版本、处理级别、传感器信息、数据时间覆盖范围、地理空间范围等。数据变量包括经度、纬度、卫星天顶角等地球位置信息，SST以及10 m风速、海冰分数等地球物理变量。质量控制信息提供了SST数据的整体质量级别，分为0—5共6个级别，0 表示没有数据，1 表示无效数据，5 表示质量最高级别。本文选用最高质量级别数据。

3 HY-1D COCTS 亮温数据条纹噪声去除

HY-1D COCTS 是扫帚式扫描辐射计，共有8 个探测器并行扫描。由于不同探测器的光谱响应存在差异，因此其辐射量图像在沿轨方向出现了明显的条纹噪声。该条纹噪声导致辐射量数据出现较大的误差，从而会对后续反演得到的SST精度造成影响。因此，在进行SST反演之前，本文先对HY-1D COCTS的亮温数据进行条纹噪声去除。

本文采用单方向变分迭代算法进行条纹噪声的去除，该算法曾成功地应用于MODIS 条纹噪声去除中（Bouali 和Ignatov，2013）。条纹噪声不影响图像的水平梯度，所以可以被假定为单方向噪声。在Bouali和Ignatov（2013）的研究中，作者提出条纹噪声的去除可以看作是一个基于单向变分模型最小化的优化问题。基于Gauss-Seidel 不动点迭代算法，可以得到Euler-Lagrange 方程的解，从而实现条纹噪声的去除（Bouali 和Ignatov，2013）。图1 和图2 展示了COCTS 11μm 通道亮温图像条纹噪声去除前后的变化情况。图1 为11μm 通道区域亮温图像，图2 为图1 黑色线位置处的亮温在沿轨方向的变化，其中图2（a）表示原始亮温，图2（b）表示条纹噪声去除后的亮温。类似的，图3 和图4展示了COCTS 12 μm 通道亮温图像条纹噪声去除前后的变化情况。通过两通道亮温图像和沿轨方向亮温的波动情况的对比，可以看出，在采用了该算法后，COCTS 亮温数据的条纹噪声被明显减弱，这对于后续反演得到高精度的SST 非常重要。

图1 HY-1D COCTS 11 μm 通道亮温图Fig. 1 The demo of HY-1D COCTS 11 μm brightness temperature

图2 HY-1D COCTS 11 μm 通道亮温在沿轨方向的变化图Fig. 2 Variation of HY-1D COCTS 11 μm brightness temperature along track

图3 HY-1D COCTS 12 μm 通道亮温图Fig. 3 The demo of HY-1D COCTS 12 μm brightness temperature

图4 HY-1D COCTS 12 μm 通道亮温在沿轨方向的变化图Fig. 4 Variation of HY-1D COCTS 12 μm brightness temperature along track

4 贝叶斯云检测

在进行COCTS SST 反演之前，首先需要对COCTS 热红外通道亮温图像进行云检测。考虑到贝叶斯云检测算法可以克服传统阈值云检测算法需要根据传感器性能、云类型或大气状态的改变而调整阈值的缺陷，本文采用贝叶斯云检测算法对COCTS 热红外通道进行云检测（Merchant 等，2019）。贝叶斯云检测算法就是基于大气辐射传输模型和NWP 海表温度和大气参数背景场信息，采用贝叶斯理论实现云检测（Merchant等，2019）。

在给定卫星观测值yo和背景场信息xb的情况下，根据贝叶斯理论，晴空概率P(c|yo，xb)可以用式（1）计算得到（Merchant等，2019）：

式中，P代表概率或者概率密度函数，c代表晴空海洋，y是卫星观测值，x是状态向量，上标o 代表观测信息，上标b 代表背景信息（先验参数）。本文用到的卫星观测数据包括11 μm和12 μm通道观测亮温，以及11 μm 通道亮温的区域标准偏差，即以该象元为中心的3 × 3 矩阵的标准偏差。背景场状态向量包括背景场海表温度和背景场大气水汽总量，二者均来自于ERA5 数据。根据式（1），在给定先验晴空概率的估计值P(c)以及在晴空和有云状态下的背景场信息的情况下，可以计算得到卫星观测的晴空概率。

由于卫星观测向量yo包括光谱信息（即通道亮温）和纹理信息（即区域标准偏差），且假设区域标准偏差的概率密度函数与观测亮温是不相关的，因此晴空概率和有云概率均可以表示为光谱部分与纹理部分晴空概率密度函数的乘积。晴空概率的光谱部分可以用高斯分布的概率密度分布函数计算得到，而晴空概率的纹理部分、有云概率的光谱部分以及有云概率的纹理部分都不满足高斯分布，因此采用基于经验统计的查找表计算得到。ESA CCI 项目提供了基于AATSR 观测的概率密度查找表（Bulgin，2019），可用于COCTS 的云检测。但是由于COCTS 的噪声水平和AATSR 的不同，因此，纹理部分的概率密度查找表需要进行调整后才能使用。将AATSR 的概率密度分布函数与COCTS的噪声进行卷积计算，得到了适用于COCTS的概率密度分布查找表。图5（a）和（b）分别为调整后的COCTS 纹理部分晴空和有云情况下的概率密度分布。对于光谱部分的概率密度查找表，主要包括表1 所示的参数。图6 展示了在给定大气路径长度（1—1.35）和NWP SST（280—281 K）区间时，概率密度在不同“11 μm 通道亮温与NWP SST 差”和“11 μm 与12 μm 通道亮温差”情况下的分布情况。因此在给定11 μm 通道亮温、12 μm 通道亮温、背景场SST、卫星天顶角和太阳天顶角的情况下，就可以根据该查找表计算得到晴空条件下光谱部分的概率。P(cˉ)可以通过NWP 数据中包括的云覆盖量比例TCC（Total Cloudy Cover）来确定（Merchant 等，2019）。如果TCC 的值在0.5—0.95，则将TCC 值赋给P(cˉ)；如果TCC 小 于0.5，则P(cˉ) = 0.5；如 果TCC 大 于0.95，P(cˉ) = 0.95（Merchant等，2019）。

图5 大气路径长度在1—1.35时纹理部分晴空概率密度Fig. 5 The clear-sky probability density of texture part when the length of atmospheric path is between 1—1.35

图6 大气路径长度在1—1.35、NWP SST在280—281 K，有云情况下的光谱部分晴空概率密度Fig. 6 The clear-sky probability density when the length of atmospheric path is between 1—1.35 and NWP SST is between 280—281 K under cloudy condition

表1 AATSR两通道概率密度函数查找表参数（Bulgin，2019）Table 1 The parameters of AATSR two channels’ probability density function lookup table（Bulgin，2019）

本文将最终计算得到的P(c|yo，xb)大于0.9 的像元视为晴空像元。图7（a）和（b）分别为2011 年5 月7 日COCTS 11 μm 通道云检测前后的条带亮温。其中，灰色区域表示为陆地，图7（b）中的白色表示贝叶斯云检测算法检测到的云像元。可以看出，云检测后的结果没有明显漏检的情况，云检测效果较为理想。另外，反演得到的SST的印证结果也可以对云检测效果进行间接评估。本文将在第6部分进行介绍。

图7 HY-1D COCTS 11 μm 通道亮温Fig. 7 HY-1D COCTS 11 μm brightness temperature

5 最优估计海表温度反演

最优估计SST反演算法是在给定大气和海洋先验状态值的前提下，结合卫星观测亮温和大气辐射传输模拟亮温以及其不确定性信息，对先验背景状态值和真实状态值的偏差进行估计，从而实现对真实状态值的估计（Merchant 等，2008，2009，2013）。最优估计SST 可以通过式（2）得到（Merchant等，2008，2009，2013）：

式中，xa为先验状态值，包括先验SST 和大气水汽总量TCWV（Total Column Water Vapor）。yo为COCTS 11 μm 和12 μm 的观测亮温，F(xa)为在先验状态为xa时，大气辐射传输模拟计算得到的亮温。本文采用中分辨率大气辐射模型MODTRAN（MODerate resolution atmospheric TRANsmission）用于模拟亮温的计算。K为偏导数矩阵，Sε为协方差矩阵，表示在忽略背景场信息误差的假设下，模拟亮温与观测亮温差的协方差总和。本文将11 μm 和12 μm 通道的模拟误差设置为0.15 K，并假设两个通道误差之间不存在相关性。根据COCTS 的传感器参数，11 μm 和12 μm 通道的等效噪声温差均为0.2 K。Sa为先验状态向量的协方差矩阵，包括先验SST 背景场的误差标准偏差eSSTa和先验大气水汽含量的误差标准偏差ewa。对于ewa的设置，本文参考Merchant 等（2013）中的参考值。eSSTa的设置根据HY-1B COCTS 最优估计SST 反演算法，确定为1.2 K。对于最优估计SST 反演算法，卫星观测值与模拟值的一致性可以作为评估SST置信指数的一个有效标准。Rodgers（2000）提出了最优估计SST置信指数的计算方法，如式（3）所示：

图8 HY-1D COCTS 2021年1月27日SST （审图号：GS鲁（2022）0019号）Fig. 8 HY-1D COCTS SST on January 27， 2021

6 HY-1D COCTS海表温度印证

本文采用iQuam 现场实测数据和VIIRS SST 对反演得到的COCTS SST 进行印证。对HY-1D COCTS SST 和VIIRS SST 均进行投影，投影方式为0.01°×0.01°。COCTS 与iQuam 以及COCTS 与VIIRS SST 分别进行匹配，匹配的空间窗口为0.01°，时间窗口为1 h。表2 给出了COCTS 与iQuam SST 以及COCTS 与VIIRS SST 的比较结果。可以看出，COCTS SST 与VIIRS SST 偏差的平均值为-0.05℃，标准偏差为0.49℃。COCTS SST 与iQuam SST 的比较结果表明，二者偏差的平均值为-0.04℃，标准偏差为0.45℃。

表2 HY-1D COCTS SST与VIIRS SST/iQuam SST 的比较结果Table 2 Comparison results of HY-1D COCTS SST with VIIRS SST / iQuam SST

图9 为COCTS 最优 估计SST 与iQuam SST 偏差随iQuam SST 的变化。其中紫色曲线表示平均偏差的变化，黑色竖线表示在iQuam SST 在每隔1 ℃区间内，匹配点的标准偏差。背景颜色表示匹配点的数目。可以看出，多数匹配点位于暖水区，且在这个SST 范围区间内，二者SST 平均偏差的浮动范围在±0.5℃之间，没有明显的与SST 的依赖关系。图10 为COCTS 最优估计SST 与iQuam SST 偏差随大气水汽总量的关系，可以看出，在整个大气水汽总量的范围区间内，二者SST的平均偏差保持平稳，标准偏差在低水汽含量范围内较高水汽含量范围内有少许偏大，整体来讲，二者SST偏差没有明显的与大气水汽含量的依赖关系，表明该反演算法可以较好地进行大气校正。图11 为COCTS 最优估计SST 与iQuam SST 所有匹配点的SST 偏差的空间分布图（审图号为“GS 鲁（2022）0019号”）。83.98%的匹配点的SST偏差在-0.5℃—0.5℃。在30°N 以北的近岸海域，COCTS SST 较iQuam SST 偏高；而在近赤道海域，二者负偏差的匹配点较多，约在-0.3℃左右。图12 为COCTS 最优估计SST 与VIIRS SST 所有匹配点的SST 偏差的直方图。77.70%的匹配点的SST 偏差在-0.5℃—0.5℃。SST 偏差基本呈现正态分布，在-0.1℃左右的匹配点最多，约为4100000 个。SST 偏差的冷尾不明显，表明COCTS云检测不存在明显的漏检。

图9 HY-1D COCTS SST与iQuam SST的偏差随iQuam SST的变化Fig. 9 The variation of HY-1D COCTS SST minus iQuam SST difference against iQuam SST

图10 HY-1D COCTS SST与iQuam SST的偏差随大气水汽总量的变化Fig. 10 The variation of HY-1D COCTS SST minus iQuam SST difference against TCWV

图11 HY-1D COCTS SST与iQuam SST的偏差（审图号：GS鲁（2022）0019号）Fig. 11 The geolocation distribution of HY-1D COCTS SST minus iQuam SST difference

图12 HY-1D COCTS SST与VIIRS SST的偏差的直方图Fig. 12 The histogram of HY-1D COCTS SST minus VIIRS SST difference

总的来说，采用贝叶斯云检测算法和最优估计SST 反演算法得到的HY-1D COCTS SST 的精度较高，与国际业务化海表温度产品达到相当水平。

7 结 论

本文从HY-1D COCTS L1级数据出发，经过条纹噪声去除、云检测和最优估计SST 反演3 个步骤，反演得到了HY-1D COCTS SST。选择的单方向变分迭代算法有效地去除了由于探测器响应不一致造成的条纹噪声，对后续获取高精度的SST数据打下了基础。结合大气辐射传输模拟计算和背景场信息，采用贝叶斯理论对COCTS 的亮温数据进行云检测。云检测后的图像中未发现明显的云漏检或者过检情况，并且后续SST的精度评估结果也表明本文中的云检测算法效果理想。最后选择最优估计算法反演得到SST，并采用iQuam SST 和VIIRS SST 对COCTS SST 进行印证。印证结果表明，COCTS 与iQuam 的平均偏差为-0.04℃，标准偏差为0.45℃；COCTS 与VIIRS 的平均偏差为-0.05℃，标准偏差为0.49℃。以上印证结果表明，HY-1D COCTS 最优估计SST 的精度较高，已到达国际业务化海表温度产品的水平。由于MODTRAN模拟计算耗时较长，因此本文将西北太平洋海域作为研究测试海域。在未来的研究中会将考虑使用快速模式进行大气辐射传输模拟，将研究区域扩展到全球海域，对HY1D COCTS进行全球SST反演与印证。

志 谢HY-1D COCTS 辐射量数据由国家卫星海洋应用中心提供；iQuam 现场实测SST数据由NOAA 提供；VIIRS SST 数据由NOAA 提供；ERA5数值天气预报数据由ECMWF 提供。在此对数据提供方表示感谢！