刘建阳，毛志华, ，陶邦一，马力，朱乾坤，黄海清，刘建强，丁静

(1.上海交通大学 海洋学院，上海 200240；2.自然资源部第二海洋研究所 卫星海洋环境动力学国家重点实验室，浙江杭州 310012；3.国家卫星海洋应用中心，北京 100081)

1 引言

20 世纪70 年代美国发射了第一代海色探测器海岸带水色扫描仪（CZCS），用于监测海洋水色和浮游生物情况。搭载中分辨率成像光谱仪（MODIS）的两颗卫星Terra 和Aqua 分别于1999 年12 月、2002 年5 月发射成功，至今仍然在轨运行。MODIS 共有36 个光谱波段，其中分辨率为250 m 的波段有2 个，分辨率为500 m 的波段有5 个，分辨率为1 km 的波段有29 个，其在110°的视场角内实现2 330 km 的扫描宽度[1]。Nishihama 等[2]在1997 年提出MODIS 几何定位算法，从焦平面位置开始，经过转换与计算得到探测器坐标系下视矢量，再经过卫星（SC）坐标系、轨道（ORB）坐标系、地心惯性（ECI）坐标系、地心旋转（ECR）坐标系转换，得到ECR 视矢量，结合卫星ECR 位置，计算ECR 视矢量与地球的交点，最后转换交点为地理经纬度。

HY-1A/B 卫星水色仪的几何定位方法利用下行数据中的GPS 数据和轨道拟合技术得到卫星轨道参数，计算卫星在扫描时刻的ECI 位置和速度，结合水色扫描仪工作原理和球面几何关系，求解扫描点在ECI 下的位置，过程中涉及SC 坐标系、ORB 坐标系、ECI 坐标系的转换，然后再转到ECR 坐标系，最后计算出扫描点的地理经纬度[3]。HY-1A/B 卫星水色仪与MODIS 几何定位的区别在于，一是利用GPS 数据和轨道拟合技术来计算卫星的位置和速度，二是未采用视矢量。

风云三号（FY-3）卫星微波成像仪采用天线绕轴旋转形成圆锥形跨轨的扫描方式，GPS 接收机提供卫星三维位置，利用两组以上卫星位置计算出卫星的ECI 实时速度。根据探测器观测几何、探测器空间位置和指向，建立观测像元与地面位置的关系模型。先将计算的天线视矢量经过坐标系转换，生成ECR 视矢量，然后结合卫星ECI 位置到ECR 位置的转换，即可计算出ECR 视矢量与地球的交点，最后再把交点转换成地理经纬度[4]。不同于MODIS，FY-3 微波成像仪根据GPS 数据中卫星位置计算速度。

以上几种探测器的几何定位方法一致，关键过程是计算出探测器视线与地球椭球体的交叉点，过程涉及坐标系转换与计算，最后转化为地理经纬度。该几何定位方法同样应用于陆地遥感卫星多光谱扫描仪（MSS）、可见光红外成像辐射仪（VIIRS）、甚高分辨率扫描辐射计（AVHRR）等探测器中[5]，应用过程中存在各种变形和改进，如传感器视线呈现形式不同，获得卫星位置速度的方式不同，初始视矢量所在坐标系不同，坐标系变换步骤不同，经纬度计算公式不同等。

潘德炉等[6]在1997 年提出了基于卫星轨道参数和探测器观测角的几何定位方法，用于卫星仿真，既不采用地面控制点，也不采用插值法近似逼近。根据先验知识获得的已知点作为地理参考点，利用卫星轨道参数和探测器观测角等参数把图像平面坐标依次转化为以星下点为球心的球面坐标、以星下点为极点的极坐标、以北极为极点的极坐标，最后转化为像元地理坐标，未采用常规坐标系。研究过程中通过逐步逼近法，推算出所有行的中心像元点地理坐标和卫星飞行方向，以此减少误差。

2018 年9 月7 日，我国在太原卫星发射中心成功发射HY-1C 卫星，2020 年6 月11 日又成功发射HY-1D卫星，两星组网，进一步提高了全球海洋水色、海洋带资源与生态环境的有效观测能力。作为HY-1A/HY-1B 卫星的后继者，HY-1C/D 卫星的正常运行标志着海洋一号系列卫星业务化运行的成功。HY-1C/D 卫星搭载中国海洋水色水温扫描仪（COCTS），COCTS几何定位是遥感资料预处理的核心要求，需要通过方法研究来提高定位精度，精确的地理定位数据可以提高COCTS 后续产品的质量。

2 COCTS 观测特性

2.1 COCTS 介绍

HY-1C/D 卫星搭载5 个有效载荷，本文重点研究其中的水色仪COCTS。COCTS 重达50 kg，功率为70 W，数据传输速率为670 kb/s。COCTS 以1 km 的空间分辨率、114°的全视域角、2 900 km 的地面幅宽扫描监测中国邻海以及全球海洋，提供了海洋水色、植被产品以及海表温度等数据。沿着HY-1C/D 卫星飞行方向，COCTS 从右向左扫描，四元同步逐点摆扫，每元有1 664 个采样点，每个采样点有10 个波段（B1-B10），包括8 个可见波段和近红外波段、2 个热红外波段。水色仪地球区图像数据传输格式如表1 所示。

表1 COCTS 地球区图像数据传输格式Table 1 COCTS earth area image data transmission format

2.2 地面景物与COCTS 成像几何关系

COCTS 采用45°扫描镜和四元探测器扫描成像形式，其扫描方向表示为+Y，指向地心方向表示为+Z，飞行方向表示为+X，垂直于Y轴、Z轴所在平面指向外侧，建立的坐标系如图1 所示。COCTS 每360°扫描周期内都有一段对地球区进行扫描成像，图1 中地球区以星下点为中心，左右各转57°，即从成像开始逆时针旋转114°的区域。COCTS 安装在卫星载荷舱对地面，依靠45°扫描镜穿轨方向的扫描和卫星的飞行，获取地球区图像数据，10 波段焦平面与地面景物成像几何关系如图1，获取的波段数据将用于后续遥感图的绘制。

图1 水色仪COCTS 周期扫描成像几何关系示意Fig.1 Diagram of the geometric relation of COCTS periodic scanning imaging

3 几何定位方法

本文提出的COCTS 几何定位方法算法流程如图2所示。首先，构建ORB 到ECR 的转换矩阵TECR/ORB。从0 级数据中提取星下点计数值，计算出星下点时间，根据相邻采样时间关系计算出各点采样时间。在提取的星历数据中采用插值法获取各采样时间对应的ECR 卫星位置和速度，通过在ECR 坐标系中建立ORB 坐标系来构造ORB 到ECR 的转换矩阵TECR/ORB。

图2 COCTS 几何定位方法流程Fig.2 Flow chart of COCTS geometric positioning method

其次，计算扫描行各采样点的ORB 视矢量。探测器（INST）坐标系下原点视矢量为OZ=[0,0,1]，忽略仪器安装误差，理论上认为安装矩阵是单位矩阵。INST 到SC 坐标系转换矩阵TSC/INST为单位矩阵，热胀冷缩是影响该矩阵动态误差的因素，误差相对较小且难被分析，可忽略。SC 到ORB 坐标系转换矩阵TORB/SC由卫星姿态决定，使用HY-1C/D 卫星双敏感器与陀螺的搭配来确定卫星姿态，精度更高，稳定性更好，故本文采用基于扩展卡尔曼滤波（EKF）姿态确定算法，分成两个子系统进行研究[7]。由于此方法难度较大，本文暂时将TORB/SC简单视为单位矩阵，所以ORB 坐标系原点视矢量为OZ=[0,0,1]，扫描行的各采样点ORB 视矢量可以通过将原点视矢量绕X、Y轴旋转相应角度的方式来获得。

然后，计算采样点ECR 视矢量与地球椭球体的交叉点，并转换为地理经纬度。结合转换矩阵TECR/ORB，进一步计算出各采样点的ECR 视矢量，建立ECR 视矢量与地球椭球体交叉模型，从而求解出交叉点ECR 坐标，利用英国学者鲍林的研究思路导出的公式即可转换交叉点ECR 坐标为地理经纬度[8]。

最后，结合波段数据绘制遥感图。根据上述过程计算的各采样点经纬度以及提取的对应波段数据，设置像元尺度，投影各采样点，各投影像元匹配对应采样点的平均波段数据，据此可以绘制遥感图。

3.1 坐标系及转换矩阵 TECR/ORB的 构建

在水色仪COCTS 几何定位算法中，涉及到3 个主要坐标系，分别是ORB 坐标系、ECR 坐标系、GEO坐标系，坐标系及几何定位示意图如图3 所示。ORB坐标系基于惯性空间的卫星位置，以卫星质心为中心，Z轴从卫星质心指向地球质心，Y轴是Z轴与瞬时速度矢量的标准叉乘，X轴为Y轴与Z轴的标准叉乘。

图3 坐标系及几何定位示意Fig.3 Coordinate system and geometric positioning diagram

ECR 坐标系以地球质心为起点，X轴从地球质心指向格林威治子午线与赤道的交点，Z轴从地球质心指向北半球极点，Y轴为Z轴与X轴的叉乘，ECR 坐标系是一个空间直角坐标系，所以某点空间位置可以用该点在坐标轴上的投影值来表示。

GEO 坐标系是基于WGS-84，用经度、纬度表示坐标的坐标系。经度定义为本地子午线与格林威治子午线的夹角，格林威治子午线向东为正，称为东经，向西为负，称为西经，变化范围是环全球经度；纬度定义为椭球体法线与赤道的夹角，向北为正，称为北纬，向南为负，称为南纬，变化范围是90°S～90°N[8]。

传统的ORB 到ECR 坐标系的转换需要两个步骤，一是从ORB 转换到ECI，该过程需要用到卫星的ECI位置和速度；二是从ECI 转换到ECR，该过程需要考虑地球的自转、极移、章动、进动等因素的影响[9]。由于地球方位信息动态变化，所以需要实时从国际地球自转和参考系服务网站获取更新的信息。

传统方法较为复杂，考虑因素较多，并且需要实时更新信息，因此采用一种简捷、不需实时获取更新信息的直接转换方法显得尤为必要。本文根据卫星的ECR 位置和速度，建立了ORB 与ECR 的联系，本质上是ECR 坐标系分别绕着3 个坐标轴旋转一定角度，最终转换到ORB 坐标系。假设某一时刻卫星的ECR 位置和速度分别为PECR、VECR，则ORB 坐标系的单位矢量计算如下[10]：

则ORB 到ECR 的转换矩阵为TECR/ORB=[XORBYORBZORB]。

3.2 插值算法

星下点位置是地球采样区的中间位置，从0 级辅助数据中提取星下点时间计数值，计算出星下点时间，相邻采样点间隔124 μs，所以根据相邻采样时间关系可以推导出各点采样时间，本节插值算法即是针对采样时间进行的。从0 级GPS 定位广播数据中提取并计算出GPS 绝对定位时间以及对应的卫星ECR位置和速度，但是由于重复的数据较多，无法直接用于后续的插值算法，所以需要先进行重复数据的剔除处理。

卫星位置的每一个元素被当作独立的一维函数，速度作为其函数的一阶导数；每组位置和速度组合起来进行插值，从而保证插值位置和速度的一致性。从0 级辅助数据中提取的GPS 时间都是整数时间，相邻间隔为1 s，每个时间对应一组三维ECR 位置和速度，使用插值法计算各采样时间对应的卫星ECR 位置和速度，在GPS 时间范围外的采样时间无法进行插值计算。插值法示意如图4 所示。

图4 插值法示意Fig.4 Diagram of insertion method

为计算采样时间Ts对应的卫星ECR 位置和速度，在参考时间中找到采样时间Ts的两个最近相邻整数时间T1、T2，dT是相邻时间T1、T2之差，根据时间T1、T2及其对应的位置P1、P2和速度V1、V2，计算多项式系数a0、a1、a2、a3，公式如下[11]：

两个位置P1、P2是 连续相邻的，两个速度V1、V2也是连续相邻的。为了插入到合适的点，对采样时间Ts进行归一化处理，转化为0 与1 之间的数值：

最后采样时间对应的卫星ECR 位置、速度计算公式如下：

3.3 各采样点ECR 视矢量的构建

针对ORB 坐标系下视矢量的求解，以下降过程为例加以说明，如图5 所示，交点表示水色仪扫描地面所采样的4×1 664 个点。沿着卫星运行方向，COCTS探头从右向左扫描，图5 中最左边是第1 列，最右边是最后一列；每扫描行有四元即4 个探头，从下向上分别是第一元到第四元。坐标原点设置在正中间位置，卫星运行方向表示X轴正方向，扫描方向表示Y轴正方向，构建图5 所示的右正左负、下正上负坐标系。假设Z轴垂直X、Y轴所在平面指向地面，则水色仪在原点处ORB 视矢量可视为OZ=[0,0,1]。

图5 水色仪降COCTS 轨扫描示意Fig.5 Diagram of COCTS downward scanning

扫描行每元都有1 664 个采样点，各列采样通过COCTS 探头在太空中旋转一定角度来实现。以星下点为中心，COCTS 探头向左、向右最大旋转角度设计值都是57°，在空中的整个视域角为114°，实际在轨运行时视域角约为116°。COCTS 旋转扫描一圈的周期为640 ms，每个周期对应360°，相邻采样点扫描间隔为124 μs，那么相邻列之间旋转角roll=124/(640 000)×2π，原点ORB 视矢量需要绕X轴旋转角度XR=roll×[831.5,-1,-831.5]来获得各列视矢量；X轴方向上的最小旋转角度pitch=0.001 38，原点ORB 视矢量需要绕Y轴旋转角度YR=pitch×[1.5,0.5,-0.5,-1.5]来获得各元视矢量。所以各采样点ORB 视矢量根据原点视矢量OZ 分别绕X轴、Y轴旋转一定角度来获得[12]，每扫描行共计4×1 664 个视矢量，即为ORB 坐标系下的视矢量uORB。

利用3.1 节构建的ORB 到ECR 转换矩阵TECR/ORB，则ECR 坐标系下各采样点视矢量计算如下：

3.4 ECR 视矢量与地球交叉算法

利用地球椭球体长半轴a、短半轴b重新调节ECR 视矢量uECR和 卫星位置矢量pECR，下列等式未考虑光传输时间和因卫星移动或相对效应造成的异常，这将导致一点系统性的偏差。

使用二次方程式求解与单位球面交叉的视矢量u′的 缩放系数d，下式d是两个交点中靠近卫星的一个解，即较小的那个解，再使用d计算椭球体交叉矢量x。

所以，根据ECR 视矢量uECR、卫星位置矢量pECR以及推算出的缩放系数d，可以计算出ECR 视矢量与地球椭球体交叉点坐标pECR+duECR，即在ECR 空间直角坐标系下的坐标。

3.5 ECR 到GEO 坐标系转换

ECR 到GEO 坐标系转换是根据英国科学家鲍林研究思路推演而成的计算公式，把ECR 坐标系下的坐标（X，Y，Z）转化为大地坐标，即地理经纬度[8]。a、b分别为地球椭球体的长半轴、短半轴；f为椭球扁率；据此计算椭球第一偏心率平方e2、第二偏心率平方e′2：

则大地经度（L on ）、纬度（Lat）可表示为

本文未考虑地形对于几何定位精度的影响，根据后文对几何定位结果的验证，可以得知地形校正前该几何定位方法对经纬度误差的影响程度，并且在后续的研究中使用数字高程模型DEM 对地形加以校正，从而减少误差。

4 结果验证

在环全球经度范围，纬度范围90°S～90°N 的地图上，进行网格化处理，横向经度分成36 000 个格子，纵向纬度分成18 000 个格子，最终形成方格边长为0.01°的36 000×18 000 网格，表征经度范围为360°、纬度范围为180°的投影地图。上述几何定位方法计算出的经纬度，可以转化为在网格中的行列号，再转化行列号为在网格中的序号，网格按从上向下，从左向右顺序标序号。把投影在同一网格序号中的采样点的R、G、B 波段数据分别取平均，表示该网格序号对应的三波段数据，最后显示得到赤道区域分辨率约为1 km 的彩色遥感图。

选 择COCTS 从2020 年10 月27 日 02:33:00 UTC开始扫描到12:35:00 UTC 结束的0 级数据，卫星在该段时间内共绕地运行6 轨，扫描范围可以覆盖欧亚非区域。按上述方法计算该区域的经纬度以及提取波段数据，作图6a 所示欧亚非大陆遥感图，图6a 为1 km分辨率遥感图，1 km 高分辨率体现地图的精细化信息，根据遥感图估算的经纬度数据也更加准确。MODIS 1 级产品光谱分辨率有250 m、500 m、1 000 m 3 种模式，图6b 为相同区域的MODIS 遥感图。从图6 可以看出，COCTS 与MODIS 遥感图轮廓大体一致。

图6 欧亚非大陆遥感图Fig.6 Remote sensing image of Europe,Asia and Africa

4.1 HY-1C 卫星COCTS 与MODIS 海岸线定性比较

选择阿拉伯半岛和渤海湾两区域进行海岸线定性验证。阿拉伯半岛海陆分界线轮廓清晰，经过Canny 边缘检测、人工优化处理即可提取海岸线，然后进行后续的轮廓匹配；而渤海湾区域具有光谱反射率低、上空多云多雾等特征，造成陆海分界线在视觉上不明显，需要先对波段数据进行归一化植被指数（NDVI）处理，以此增强海陆色度差，提高区分度，再实施Canny 边缘检测、人工优化来提取海岸线[13]，最后对提取的海岸线进行轮廓匹配。

1）NDVI 处理

NDVI 是衡量植被覆盖率的一个参数，根据NIR和R 两个波段反射值之差与之和的比来计算，范围处于-1 与1 之间。根据R、NIR 波段范围，从COCTS的10 个波段中分别选择第6、8 波段来计算NDVI，正值表征植被，植被覆盖率越高，NDVI 越大。通过对NDVI 进行适合阈值设定，从而剔除绝大部分白云、水、岩石、裸土以及绿色植被覆盖较少的地方，剩下的深绿区域植被覆盖率较高[14-15]。

2）Canny 边缘检测

首先，高斯滤波确定高斯核的尺寸为5×5，标准差为0.5，与图像进行离散卷积；其次，考虑到Sobel 算子边缘强弱分明，抗噪声好，选定该算子计算X方向和Y方向像素梯度；然后，将当前像素梯度强度与沿正负梯度方向上的相邻像素梯度强度进行比较，抑制非极大值像素梯度，消除边缘检测带来的杂散响应；最后，设置高低阈值，定义真实边缘点[16-18]。

3）人工优化

经过前面两步处理，海陆交界处像素梯度显著，所以海岸线被提取得很好。然而海岸图仍有尚未剔除的干扰线段，所以最后需要对海岸图进行人工优化，用干扰线段附近的图像数据取代干扰区域，从而屏蔽干扰线段，得到最终清晰的海岸线。

4.1.1 区域1—阿拉伯半岛

根据MODIS 在2020 年10 月27 日的1 级产品画出1 km 空间分辨率的阿拉伯半岛地区遥感图，如图7a所示。对其进行Canny 边缘检测，调整阈值参数，得到图7b，该海岸线还存在干扰线段，因此还需进行人工优化，剔除噪声干扰，优化如图7c 或图7d所示，该图即为根据MODIS 数据最终提取的阿拉伯半岛地区海岸线。从2020 年10 月27 日 05:54:23 UTC 至09:15:17 UTC，COCTS 绕地运行两轨，覆盖阿拉伯半岛，该时间段的0 级数据按上述方法计算出1 级数据，据此作遥感图7e，与图7a 同样经纬度范围。利用CorelDRAW 软件把图7d 海岸线叠加到图7e 遥感图中，通过调整透明度进而验证海岸线吻合情况，因此从图7f 中可以定性判断MODIS 白色海岸线与HY-1C 卫星COCTS 遥感图海岸线匹配一致，二者变化趋势吻合。

图7 阿拉伯半岛海岸线提取与匹配Fig.7 Extraction and match of coastline in the Arabian Peninsula

4.1.2 区域2—渤海湾

COCTS 从2020 年10 月19 日 02:12:25 UTC 至02:22:37 UTC 扫描探测渤海湾区域，该时间段对应的0 级数据按上述几何定位方法计算出经纬度，提取波段信息，据此作1 km 分辨率遥感图，依次进行图8 所示处理。图8e 为根据同一天的MODIS 1 级产品作同一经纬区域的遥感图，图8f 为COCTS 提取的海岸线与MODIS 遥感图海岸线轮廓匹配结果，因此可以定性判断COCTS 白色海岸线与MODIS 遥感图海岸线匹配一致，二者变化趋势吻合。

图8 渤海湾海岸线提取与匹配Fig.8 Extraction and match of coastline in the Bohai Gulf

4.2 HY-1C 卫星COCTS 海岸线经纬度与GSHHG 数据库定量比较

GSHHG 是全球一致、分层、高分辨率地理数据库，由3 个数据库构成：世界矢量海岸线、世界数据库II 和冰冻圈地图集，分别是海岸线、湖泊和南极海岸线的基础。GSHHG 地理数据分辨率有5 个级别，选取最高分辨率“full”模式。因为“high”模式空间分辨率为200 m，“full”模式所占磁盘空间是“high”模式的4 倍多，所以“full”分辨率是优于50 m 的[19]。为了检验几何定位的精度，根据几何定位方法计算的结果作1 km 分辨率遥感图，从GSHHG 数据库中提取相同区域的海岸线坐标，在该遥感图上叠加相应的海陆分界线，该海路分界线即为参考，通过其与实际遥感图海岸线的吻合程度来判断几何定位的精度。

选择上述COCTS 在2020 年10 月27 日扫描覆盖阿拉伯半岛的数据，作1 km 分辨率遥感图（图9a）。选定阿拉伯半岛区域重点研究，用海岸线提取工具，选定始末坐标，开始于35.91°N，30.78°E，结束于11.34°N，60.84°E。该始末坐标范围内的海岸线坐标都可以提取出来，叠加到图9a 中，红色即为参考海岸线。选取中间区域的波斯湾北部，进一步放大如图9b 所示；选取两侧区域的波斯湾东南部，进一步放大如图9c 所示，由此二图可以清晰看到像元级别。由于COCTS的扫描特征，越靠近扫描区域的两侧，采样越稀疏，采样不足导致存在无效数据点，在遥感图上体现为黑色像元点。红色参考海岸线与实际作的遥感图海岸线吻合程度反映了几何定位精度，如图9 所示，COCTS几何定位方法计算的经纬度误差在2 个像元内，即2 km 内。

图9 阿拉伯半岛遥感图参考海岸线叠加（a），星下点附近中间放大区域（b），星下点两侧边缘放大区域（c）Fig.9 Reference coastline overlaying Arabian remote sensing image (a),middle enlargement area near Nadir (b),marginal enlargement area away from Nadir (c)

4.3 HY-1C 卫星COCTS 星下点经纬度与卫星工具软件结果定量比较

成熟的卫星工具软件（Satellite Tool Kit，STK）是航天领域先进的系统分析软件，用于分析复杂的陆地、海洋等任务以及提供精确的分析报告。高精度轨道生成模块考虑点重力模型、日月重力影响、大气阻力、章动、自旋、质心变化等很多因素，以此确保生成高精度轨道。

选择某一区域数据，从0 级文件中提取GPS 时间及对应的位置和速度，利用STK 软件的高精度轨道生成模块模拟卫星轨道，然后根据输入的星下点时间，输出对应的卫星位置和速度。利用获得的星下点时间、位置、速度在STK 软件上再次模拟出卫星轨道，星下点作为特征点，一个星下点时间对应一组星下点经纬度。然后导出LLA 文件，文件中包括星下点经纬度[20-21]，该数据作为参考标准，与本文几何定位方法计算的星下点经纬度进行比较和分析。

首先选择COCTS 从2020 年10 月19 日 02:12:25 UTC 至02:22:37 UTC 扫描覆盖渤海湾区域的0 级数据，利用上述定位方法计算957 组星下点经纬度，STK输出对应的参考值。计算值与参考值对比结果为，误差在0.001°～0.01°之间的星下点占22.46%，误差在0.01°～0.02°之间的星下点占77.54%，误差在2 个像元内。再选取长时间大范围的数据，选择COCTS 从2020 年10 月27 日 07:34:50 UTC 至09:15:17 UTC 扫描经过阿拉伯半岛的一整轨0 级数据，同样方法计算出9 418 组星下点经纬度以及输出STK 对应参考值。为了进一步对比分析，把数据分成南北纬0°～30°、南北纬30°～60°、南北纬60°～90°共低、中、高纬3 部分，计算值与参考值对比，误差占比统计结果如表2 所示。分析星下点误差可知，本文几何定位方法计算的星下点误差都在0.02°内，即2 个像元内；从低纬到高纬误差变化趋势可以判断，赤道误差最小，越往两极，误差越大。

表2 计算值与参考值误差统计Table 2 Error statistics between calculated and reference values

4.4 HY-1D 卫星与HY-1C 卫星COCTS 经纬度定量比较

HY-1D 卫星与HY-1C 卫星功能一致，共同监测全球水色和海表温度，两颗卫星组网，提高全球扫描覆盖能力，共同构成中国首个海洋民用卫星星座。两星搭载同样的水色仪COCTS，COCTS 的硬件结构、运行模式和工作原理等完全一致。HY-1D 卫星与HY-1C 卫星的区别是，在遥感图上HY-1D 卫星从右下往左上运行成像，而HY-1C 是从右上往左下运行成像。

将2020 年9 月23 日 HY-1D COCTS 的遥感图与上述2020 年10 月27 日 HY-1C COCTS 的遥感图进行比较。选择波斯湾西部一个特征小区域，遥感图放大6 倍，如图10 所示，图10a 为HY-1C 特征区域，图10b 为HY-1D 卫星特征区域，两图中红色像元点为用于定量比较的选取点，把选取点在全球范围遥感图中的行列号记录在表3。经统计，在1 km 分辨率前提下，HY-1C 卫星与HY-1D 卫星的纬度、经度之差都为0 像元，所以COCTS 几何定位方法对HY-1D 卫星与HY-1C 卫星都适用，二者定位结果一致。

图10 HY-1C 卫星（a）与HY-1D 卫星（b）对应特征点采样Fig.10 Corresponding feature points sampled from HY-1C satellite (a) and HY-1D satellite (b)

表3 HY-1C/D 卫星特征区域采样点比较Table 3 Sampling points comparison of HY-1C/D satellite feature areas

5 像元尺度效应对几何定位误差影响分析

COCTS 几何定位方法的地理定位结果存在误差，本节分析像元尺度效应对几何定位误差的影响。地球是个不规则的三维椭球体，两极略尖，赤道略鼓，从而半轨成像具有两极宽、赤道窄的特点。选取COCTS 在2020 年10 月27 日扫描覆盖阿拉伯半岛的半轨经纬度数据，对于1 664 列数据，从中选择7 对相邻样本列，分别是1 和2 列、208 和209 列、416 和417 列、832 和833 列、1 247 和1 248 列、1 454 和1 455 列、1 663 和1 664 列。为了详细分析半轨经纬度数据，把半轨扫描区域分成3 部分：靠近北极的前1/3 区域、中间1/3 区域、靠近南极的后1/3 区域。分别计算各区域每对样本列的平均距离，记录如表4。根据表4中数据作7 对样本列距离变化趋势图，便于直观看出经纬度数据的特征趋势，如图11 所示，黑色、绿色、蓝色虚线分别表示前1/3 区域、中间1/3 区域、后1/3 区域，红色实线表示半轨区域，同时用柱形图表示该半轨区域的7 对样本列距离。结合水色仪COCTS的扫描原理，分析表4 数据和图11 的变化趋势，可以得出，从星下点到两侧边缘，从赤道到两极，采样像元尺度呈增大趋势。

图11 7 对样本列距离变化趋势Fig.11 Distance change trends of seven pairs of adjacent sample arrays

表4 相邻样本列距离（单位：（°））Table 4 Distance of adjacent sample arrays (unit:(°))

半轨扫描区域如图12a 中黑色线框范围，为了具体分析两侧和星下点附近中间区域的误差特征，从两侧和中间区域各选取40 个海岸特征点，分别用红色、黄色标志，进一步放大如图12b 所示。两侧和中间区域的特征点在1 km 分辨率遥感图中的行列号以及与GSHHG 数据库海岸线坐标的像元误差分别记录在表5 和表6 中，经误差统计，两侧区域的纬度方向、经度方向平均误差分别是 0.95、1.25 个像元，即误差为1.57 个像元；中间区域的纬度方向、经度方向平均误差分别是 0.30、0.65 个像元，即误差为0.72 个像元，所以扫描区域两侧的几何定位误差比中间大。理论分析结合实际结果验证可知，星下点误差最小，越接近两侧边缘，采样像元畸变越大，定位误差随之增大。因此，从星下点到两侧边缘，定位误差随采样像元尺度增大而呈增大趋势，二者呈正相关。

图12 中间及两侧边缘特征点采样（a）及放大（b）Fig.12 Feature points sampled in the middle and marginal sides (a) and enlargement (b)

表5 靠近边缘两侧区域特征点误差Table 5 Errors of feature points on both edge-nearing sides of the track

表6 中间区域特征点误差Table 6 Errors of feature points in the middle area of the track

续表 6

为了具体分析从赤道到两极的误差特征，在扫描范围内的赤道附近区域和南北纬35°附近区域分别选取20 个海岸特征点，并且特征点都分布在星下点附近区域，如图13a 所示，进一步放大如图13b 所示。赤道附近、35°N 附近和35°S 附近的特征点在遥感图中的行列号，以及特征点与GSHHG 数据库海岸线坐标的像元误差分别记录在表7、表8、表9 中，经误差统计，赤道附近特征点的纬度方向、经度方向平均误差分别是0.20、0.40 个像元，即误差为0.45 个像元；35°N 附近特征点的纬度方向、经度方向平均误差分别是0.55、0.80 个像元，即误差为0.97 个像元；35°S附近特征点的纬度方向、经度方向平均误差分别是0.50、0.70 个像元，即误差为0.86 个像元，因此扫描区域赤道附近定位误差比两侧小。理论分析结合实际结果验证可知，赤道定位误差最小，越靠近两极，采样像元畸变越大，定位误差随之增大。因此，从赤道到两极，定位误差随采样像元尺度增大而逐渐增大，二者呈正相关。

图13 赤道及两侧特征点采样（a）及放大（b）Fig.13 Feature points sampled near the equator and its both sides (a) and enlargement (b)

表7 赤道附近特征点误差Table 7 Errors of feature points near the equator

表8 35°N 附近特征点误差Table 8 Errors of feature points near 35°N

表9 35°S 附近特征点误差Table 9 Errors of feature points near 35°S

6 结论

HY-1C/D 卫星水色仪COCTS 几何定位方法是遥感资料预处理系统的重要组成部分，精确的地理定位数据即1 级产品可以提高HY-1C/D 卫星后续产品的质量，因此建立一套将水色仪0 级产品处理成1 级产品的系统，包括数据解包、几何定位，显得尤为重要。本文几何定位方法取代了传统需要根据6 个轨道参数来计算卫星位置的方法，也简化了传统ORB到ECR 的复杂转换过程，在卫星星历中采用插值法获得采样时间对应的卫星位置和速度，进而直接计算ORB 到ECR 的转换矩阵。基于COCTS 逐点摆扫、四元并扫的扫描方式以及相关参数，计算采样点ORB 视矢量，建立采样点视矢量与地球交叉点关系模型，从而对COCTS 采样像元进行地理定位。通过对HY-1C/D 卫星水色仪COCTS 几何定位结果的验证与分析，得出以下结论：

（1）针对海岸线定性验证，COCTS 与MODIS 遥感图海岸线匹配一致，二者变化趋势吻合。

（2）针对海岸线定量验证，COCTS 几何定位方法计算的海岸线经纬度与GSHHG 参考海岸线坐标相比，误差在2 个像元内，即2 km 内。

（3）针对星下点定量验证，星下点坐标计算值与STK 生成的参考值相比，误差在2 个像元内且从赤道至两极，星下点误差呈递增趋势。

（4）COCTS 几何定位方法对HY-1D 卫星与HY-1C 卫星COCTS 的定位结果一致。

（5）像元尺度效应与几何定位误差呈正相关。从星下点到两侧边缘，定位误差随采样像元尺度增大而呈增大趋势，星下点误差最小；从赤道到两极，定位误差也随采样像元尺度增大而逐渐增大，赤道定位误差最小。

综上所述，经过定性、定量验证与分析，以卫星轨道参数和水色仪COCTS 参数为基础的几何定位方法是可行的，满足一定的定位精度要求，可以用于COCTS 遥感图像预处理的几何定位。