一种基于COCTS四元逐点摆扫方式的几何定位方法

一种基于cocts四元逐点摆扫方式的几何定位方法
技术领域
1.本发明涉及一种有效载荷的几何定位方法，具体涉及一种基于中国海洋水色水温扫描仪(chinese ocean color temperature scanner，简称水色仪cocts)四元逐点摆扫方式的几何定位方法。


背景技术：

2.2018年9月7日，我国在太原卫星发射中心由长征二号丙运载火箭成功发射hy-1c卫星，2020年6月11日又成功发射hy-1d卫星，两颗卫星组网，进一步提高了全球海洋水色、海洋带资源与生态环境的有效观测能力。作为hy-1a/hy-1b的后继者，hy-1c/d的正常运行标志着海洋一号系列卫星业务化运行的成功。hy-1c/d搭载空间分辨率约为1km的水色仪cocts，cocts几何定位是遥感资料预处理的核心要求，需要通过方法研究来提高定位精度，精确的地理定位数据可以提高cocts后续产品的质量。有效载荷的几何定位方法关键过程是计算出探测器视线与地球椭球体的交叉点，过程涉及坐标系转换与计算，最后转化为地理经纬度。各有效载荷在计算过程中存在各种变形和改进，如传感器视线呈现形式不同、获得卫星位置速度的方式不同、初始视矢量所在坐标系不同、坐标系变换步骤不同、经纬度计算公式不同等等。cocts的扫描方式是四元探头并列、逐点摆扫，目前国内外尚无关于cocts这种扫描方式的几何定位方法。


技术实现要素：

3.针对现有研究方法的不足，本发明提供了一种基于cocts四元逐点摆扫方式的几何定位方法，结合卫星星历和姿态数据、cocts的参数以及卫星空间几何关系，使用地球椭球体和地形表面信息计算采样像元的地理位置。整套算法涉及有效载荷几何关系、视矢量计算、插值算法、坐标系的定义和转换、交叉点算法等，方法的核心是计算cocts各视矢量与地球椭球体的交叉点。
4.本发明为解决上述问题所提供的技术方案为：
5.一种基于cocts四元逐点摆扫方式的几何定位方法，cocts在地面的采样点构成4x1664的矩阵，沿着卫星运行方向，cocts探头从右向左扫描，每扫描行有四元即四个探头，从下向上分别是第一元到第四元，坐标原点设置在扫描区域的正中心位置，卫星运行方向表示x轴正方向，扫描方向表示y轴正方向，z轴垂直x、y轴所在平面，且z轴正方向指向地面，构建坐标系，cocts在卫星轨道坐标系orb下的原点视矢量为oz＝[0 0 1]；
[0006]
该方法包括如下步骤：
[0007]
s1：从cocts 0级数据中提取星下点计数值，计算出星下点时间，根据相邻采样时间关系计算出所有采样点的采样时间；
[0008]
s2：根据提取的卫星轨道数据，采用插值法获取各采样时间对应的地心旋转坐标系ecr下的卫星位置和速度；然后根据ecr卫星位置和速度计算orb坐标系到ecr坐标系的转换矩阵；
[0009]
s3：将cocts在orb下的原点视矢量oz绕x、y轴分别旋转相应角度，获得各采样点的orb视矢量；
[0010]
s4：根据orb坐标系到ecr坐标系的转换矩阵，将各采样点的orb视矢量转换为ecr视矢量；
[0011]
s5：计算ecr视矢量与地球椭球体的交叉点坐标(x，y，z)；
[0012]
s6：把ecr坐标系下的交叉点坐标(x，y，z)转化为大地坐标，即经纬度；
[0013]
s7：从cocts 0级数据中提取波段数据，据此绘制等经纬度遥感图。
[0014]
进一步地，所述s2中，
[0015]
(1)根据采样时间ts找到与其相邻的两个整数时间t1、t2及其对应的位置p1、p2和速度v1、v2，计算多项式系数a0、a1、a2、a3和归一化时间间隔t：
[0016]
公式如下：
[0017]
a0＝p1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0018]
a1＝dt*v1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0019]
a2＝3*(p
2-p1)-dt*(2*v1+v2)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0020]
a3＝2*(p
1-p2)+dt*(v1+v2)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0021]
其中，dt是相邻时间t1、t2之差；
[0022]
(2)对采样时间ts进行归一化处理，转化为0与1之间的数值：
[0023]
t＝(t
s-t1)/dt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0024]
(3)根据下式计算采样时间ts对应的卫星ecr位置、速度：
[0025]
p
ecr
＝a0+a1*t+a2*t2+a3*t3ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0026]vecr
＝(α1+2*a2*t+3*a3*t2)/dt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0027]
进一步地，所述s2中，orb坐标系到ecr坐标系的转换矩阵t
ecr/orb
的计算公式如下：
[0028]
t
ecr/orb
＝[x
orbyorbzorb
]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0029]zorb
＝-p
ecr
/|p
ecr
|
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0030]yorb
＝z
orb
×vecr
/|z
orb
×vecr
|
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0031]
x
orb
＝y
orb
×zorb
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0032]
进一步地，所述s3中，
[0033]
各采样点的orb的各列视矢量vecto
rco
l的计算公式如下：
[0034][0035]
xr＝roll*[831.5:-1:-831.5]
ꢀꢀꢀꢀ
(13)
[0036]
roll＝124/(640000)
×
2π
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0037]
其中，xr为原点orb视矢量oz需要绕x轴旋转的角度，roll为相邻列之间的旋转角；[831.5:-1:-831.5]表示从831.5依次递减-1至-831.5形成的1x1664的矩阵；
[0038]
各采样点的orb的四元视矢量vector
row
的计算公式如下：
[0039]
1]。
[0057]
如图2所示，本发明的基于cocts四元逐点摆扫方式的几何定位方法包括如下步骤：
[0058]
s1：从cocts 0级数据中提取星下点计数值，计算出星下点时间，根据相邻采样时间关系计算出所有采样点的采样时间；
[0059]
星下点位置是地球采样区的中间位置，从0级辅助数据中提取星下点时间计数值，计算出星下点时间，相邻采样点间隔124us，所以根据相邻采样时间关系可以推导出各点采样时间，这里的插值算法即是针对采样时间进行的。从0级gps定位广播数据中提取并计算出gps绝对定位时间以及对应的卫星ecr位置和速度，但是由于重复的数据较多，无法直接用于后续的插值算法，所以需要先进行重复数据的剔除处理。
[0060]
s2：根据提取的卫星轨道数据，采用插值法获取各采样时间对应的地心旋转坐标系ecr下的卫星位置和速度；然后根据ecr卫星位置和速度计算orb坐标系到ecr坐标系的转换矩阵；
[0061]
卫星位置的每一个元素被当作独立的一维函数，速度作为其函数的一阶导数；每组位置和速度组合起来进行插值，从而保证插值位置和速度的一致性。从0级辅助数据中提取的gps时间都是整数时间，相邻间隔1s，每个时间对应一组三维ecr位置和速度，使用插值法计算各采样时间对应的卫星ecr位置和速度，在gps时间范围外的采样时间无法进行插值计算。插值法示意如图3所示。
[0062]
为计算采样时间ts对应的卫星ecr位置和速度，在参考时间中找到采样时间ts的两个相邻整数时间t1、t2，dt是相邻时间t1、t2之差，根据时间t1、t2及其对应的位置p1、p2和速度v1、v2，计算多项式系数a0、a1、a2、a3，公式如下：
[0063]
a0＝p1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0064]
a1＝dt*v1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0065]
a2＝3*(p
2-p1)-dt*(2*v1+v2)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0066]
a3＝2*(p
1-p2)+dt*(v1+v2)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0067]
两个位置p1、p2是连续相邻的，两个速度v1、v2也是连续相邻的。为了插入到合适的点，对采样时间ts进行归一化处理，转化为0与1之间的数值：
[0068]
t＝(t
s-t1)/dt(5)
[0069]
最后采样时间对应的卫星ecr位置、速度计算公式如下：
[0070]
p
ecr
＝a0+a1*t+a2*t2+a3*t3ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0071]vecr
＝(a1+2*a2*t+3*a3*t2)/dt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0072]
然后通过在ecr坐标系中建立orb坐标系来构造orb到ecr的转换矩阵t
ecr/orb
。
[0073]
传统的orb到ecr坐标系的转换需要两个步骤，一是从orb转换到eci，该过程需要用到卫星的eci位置和速度；二是从eci转换到ecr，该过程需要考虑地球的自转、极移、章动、进动等因素的影响。由于地球方位信息动态变化，所以需要实时从国际地球自转和参考系服务网站获取更新的信息。传统方法较为复杂，考虑因素较多，并且需要实时更新信息，因此采用一种简捷、不需实时获取更新信息的直接转换方法显得尤为必要。
[0074]
本发明根据卫星的ecr位置和速度，建立了orb与ecr的联系，本质上是ecr坐标系分别绕着三个坐标轴旋转一定角度，最终转换到orb坐标系。假设某一时刻卫星的ecr位置
和速度分别为p
ecr
，v
ecr
，则orb坐标系的单位矢量计算如下：
[0075]zorb
＝-p
ecr
/|p
ecr
|
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0076]yorb
＝z
orb
×vecr
/|z
orb
×vecr
|
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0077]
x
orb
＝y
orb
×zorb
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0078]
则orb到ecr的转换矩阵t
ecr/orb
＝[x
orbyorbzorb
]。
[0079]
s3：将cocts在orb下的原点视矢量oz绕x、y轴分别旋转相应角度，获得各采样点的orb视矢量；
[0080]
如图1所示，交点表示水色仪扫描地面所采样的4x1664个点。扫描行每元都有1664个采样点，各列采样通过cocts探头在太空中旋转一定角度来实现。以星下点为中心，cocts探头向左、向右最大旋转角度设计值都是57
°
，在空中的整个视域角为114
°
，实际在轨运行时视域角约为116
°
。cocts旋转扫描一圈的周期为640ms，每个周期对应360
°
，相邻采样点扫描间隔124us，那么相邻列之间旋转角roll＝124/(640000)x2π，原点orb视矢量需要绕x轴旋转角度xr＝roll*[831.5:-1:-831.5]来获得各列视矢量；x轴方向上的最小旋转角度pitch＝0.00138，原点orb视矢量需要绕y轴旋转角度yr＝pitch*[1.5 0.5-0.5-1.5]来获得各元视矢量。所以各采样点orb视矢量根据原点视矢量oz分别绕x轴、y轴旋转一定角度来获得，每扫描行共计4x1664个视矢量，即为orb坐标系下的视矢量u
orb
。
[0081]
s4：根据orb坐标系到ecr坐标系的转换矩阵t
ecr/orb
，将各采样点的orb视矢量转换为ecr视矢量u
ecr
；
[0082]uecr
＝t
ecr/orb
×uorb
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0083]
＝[x
orbyorbzorb
]
×uorb
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0084]
s5：计算ecr视矢量与地球椭球体的交叉点坐标(x，y，z)；
[0085]
s6：把ecr坐标系下的交叉点坐标(x，y，z)转化为大地坐标，即经纬度；
[0086]
利用地球椭球体长半轴a、短半轴b重新调节ecr视矢量u
ecr
和卫星位置矢量p
ecr
，下列等式未考虑光传输时间和因卫星移动或相对效应造成的异常，这将导致一点系统性的偏差。
[0087][0088]
使用二次方程式求解与单位球面交叉的视矢量u
′
的缩放系数d，下式d是两个交点中靠近卫星的一个解，即较小的那个解，再使用d计算椭球体交叉矢量x。
[0089][0090]
x＝p
ecr
+du
ecr
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0091]
所以，根据ecr视矢量u
ecr
、卫星位置矢量p
ecr
以及推算出的缩放系数d，可以计算出ecr视矢量与地球椭球体交叉点坐标p
ecr
+du
ecr
，即在ecr空间直角坐标系下的坐标。ecr到geo(大地坐标系)转换是根据英国科学家鲍林研究思路推演而成的计算公式，把ecr坐标系下的坐标(x，y，z)转化为大地坐标，即地理经纬度。a、b分别为地球椭球体的长半轴、短半轴，f为椭球扁率，据此计算椭球第一偏心率平方e2、第二偏心率平方e
′2：
[0092]
e2＝2f-f2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0093][0094]
再根据e2、e
′2，计算辅助参数u：
[0095][0096]
则大地经度lon、纬度lat可表示为：
[0097][0098][0099]
s7：从cocts 0级数据中提取波段数据，据此绘制等经纬度遥感图。
[0100]
在经度范围-180
°
～180
°
，纬度范围-90
°
～90
°
的地图上，进行网格化处理，横向经度分成36000个格子，纵向纬度分成18000个格子，最终形成方格边长为0.01
°
的36000x18000网格，表征经度范围为360
°
、纬度范围为180
°
的投影地图。利用几何定位方法计算出的经纬度，可以转化为在网格中的行列号，再转化行列号为在网格中的序号，网格按从上向下，从左向右顺序标序号。把投影在同一网格序号中的采样点的r、g、b波段数据分别取平均，表示该网格序号对应的三波段数据，最后显示得到分辨率约为1km的彩色遥感图。
[0101]
下面通过具体的试验数据证明本发明的方法的优点。
[0102]
根据modis在2020-10-27的1级产品画出1km空间分辨率的阿拉伯半岛地区等经纬度遥感图像，如图4所示。对其进行canny边缘检测，调整阈值参数，得到图4(b)，该海岸线还存在干扰线段，因此还需进行人工优化，剔除噪声干扰，优化如图4(c)所示，该图即为根据modis数据最终提取的阿拉伯半岛地区海岸线。从2020-10-27 05:54:23utc至09:15:17utc，cocts绕地运行两轨，覆盖阿拉伯半岛，该时间段的0级数据按上述方法计算出1级数据，据此作遥感图4(d)，与图4(a)同样经纬度范围。利用coreldraw软件把图4(c)海岸线叠加到遥感图4(d)中，通过调整透明度进而验证海岸线吻合情况。因此，从图4(e)中可以定性判断，modis白色海岸线与hy-1c cocts遥感图海岸线匹配一致，二者变化趋势吻合。
[0103]
gshhg是全球一致、结构分层、分辨率较高的地理数据库，由三个数据库构成：世界矢量海岸线、世界数据库ii和冰冻圈地图集，分别是海岸线、湖泊和南极海岸线的基础。gshhg地理数据分辨率有五个级别，选取最高分辨率“full”模式。因为“high”模式空间分辨率为200m，“full”模式所占磁盘空间是“high”模式的四倍多，所以“full”分辨率是优于50m的。为了检验像元定位的精度，根据像元定位方法计算的结果作1km分辨率遥感图，从gshhg数据库中提取相同区域的海岸线坐标，在该遥感图上叠加相应的海陆分界线，该海路分界线即为参考，通过其与实际遥感图海岸线的吻合程度来判断像元定位的精度。
[0104]
选择上述cocts在2020-10-27扫描覆盖阿拉伯半岛的数据，作1km分辨率等经纬度遥感图。选定阿拉伯半岛区域重点研究，用海岸线提取工具，选定始末坐标，开始于东经30.78
°
、北纬35.91
°
，结束于东经60.84
°
、北纬11.34
°
。该始末坐标范围内的海岸线坐标都可以提取出来，叠加到等经纬度遥感图中。参考海岸线与实际作的遥感图海岸线吻合程度反映了像元定位精度，cocts像元定位方法计算的经纬度误差在2个像元内，即0.02
°
内。
[0105]
成熟的卫星工具软件stk(satellite tool kit)是航天领域先进的系统分析软
件，用于分析复杂的陆地、海洋等任务以及提供精确的分析报告。高精度轨道生成模块考虑点重力模型、日月重力影响、大气阻力、章动、自旋、质心变化等很多因素，以此确保生成高精度轨道。
[0106]
选择某一区域数据，从0级文件中提取gps时间及对应的位置和速度，利用stk软件的高精度轨道生成模块模拟卫星轨道，然后根据输入的星下点时间，输出对应的卫星位置和速度。利用获得的星下点时间、位置、速度在stk软件上再次模拟出卫星轨道,星下点作为特征点，一个星下点时间对应一组星下点经纬度。然后导出lla文件，文件中包括星下点经纬度，该数据作为参考标准，与本文几何定位方法计算的星下点经纬度进行比较和分析。
[0107]
首先选择cocts从2020-10-19 02:12:25utc至02:22:37utc扫描覆盖渤海湾区域的0级数据，利用上述定位方法计算957组星下点经纬度，stk输出对应的参考值。计算值与参考值对比结果为，误差在0.001-0.01
°
之间的星下点占22.46％，误差在0.01-0.02
°
之间的星下点占77.54％，误差在2个像元内。再选取长时间大范围的数据，选择cocts从2020-10-2707:34:50utc至09:15:17utc扫描经过阿拉伯半岛的一整轨0级数据，同样方法计算出9418组星下点经纬度以及输出stk对应参考值。为了进一步对比分析，把数据分成南北纬0-30
°
、南北纬30-60
°
、南北纬60-90
°
共低、中、高纬三部分，计算值与参考值对比，误差占比统计结果如表1所示。分析星下点误差可知，本发明的几何定位方法计算的星下点误差都在0.02
°
内，即2个像元内；从低纬到高纬误差变化趋势可以判断，赤道误差最小，越往两极，误差越大。
[0108]
表1计算值与参考值误差统计
[0109][0110]
综上所述，经过定性、定量验证与分析，以卫星轨道参数和水色仪cocts参数为基础的几何定位方法是可行的，满足一定的定位精度要求，可以用于cocts遥感图像预处理的几何定位。
[0111]
本领域普通技术人员可以理解，以上所述仅为发明的优选实例而已，并不用于限制发明，尽管参照前述实例对发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在发明的精神和原则之内，所做的修改、等同替换等均应包含在发明的保护范围之内。