一种基于太阳位置和双频组合的GNSS导航卫星PCO改正方法

技术领域

本发明涉及卫星定位技术领域，具体地说是一种基于太阳位置和双频组合的GNSS导航卫星PCO改正方法。

背景技术

GNSS是全球导航卫星系统，包含美国的GPS、中国的BDS、俄罗斯的GLONASS和欧洲的GALILEO。

GNSS卫星可向外广播发射GNSS导航信号和电文，低轨道地球卫星上搭载的GNSS接收机可接收来自GNSS卫星发射的信号，并进行捕获跟踪从而获取伪距、载波相位、信噪比等观测信息。

在GNSS遥感领域，比如GNSS无线电掩星探测、GNSS反射探测等技术中需要确定低轨卫星、导航星之间的几何关系或高度角关系，这就需要用到GNSS导航卫星的精确位置。

国际GNSS服务IGS官网或武大IGS网站公布了GNSS卫星的精密轨道和钟差信息，通常该轨道是卫星质心处的位置，一般GNSS卫星较大，其实际接收信号的点为接收天线相位中心，卫星质心离接收天线相位中心有3米左右的距离，故在精确计算GNSS卫星位置时需要对GNSS卫星的天线相位中心相对质心的偏差进行改正，同时也需要考虑到GNSS卫星姿态对其造成的影响。

现有技术中IGS官网或武大IGS中心公布的GNSS卫星轨道均为各卫星质心处的精确位置信息，不代表其接收信号时那一点的位置，使用卫星质心处位置代替GNSS信号接收天线处的位置会产生较大的误差，在GNSS掩星探测技术中，对GNSS掩星位置精度要求为20cm，该GNSS天线相位中心相对卫星质心偏差会导致计算的接收信号的位置偏差超过20cm，从而引起中性大气掩星事件反演误差过大。本专利提供了一种基于太阳位置和双频组合的GNSS导航卫星PCO改正方法。

因此，为了解决上述问题，本申请提出了一种基于太阳位置和双频组合的GNSS导航卫星PCO改正方法，解决了GNSS卫星计算不准的问题，提高计算GNSS卫星的精度，在GNSS掩星探测、GNSS反射探测以及其他低轨卫星载荷等遥感探测方面对GNSS卫星位置有较高精度要求的领域有重要应用价值。

发明内容

本发明的目的是克服现有技术的不足，提供了一种基于太阳位置和双频组合的GNSS导航卫星PCO改正方法，解决了GNSS卫星计算不准的问题，提高计算GNSS卫星的精度，在GNSS掩星探测、GNSS反射探测以及其他低轨卫星载荷等遥感探测方面对GNSS卫星位置有较高精度要求的领域有重要应用价值。

为了达到上述目的，本发明提供一种基于太阳位置和双频组合的GNSS导航卫星PCO改正方法，包括以下步骤：

S1，读取GNSS卫星PCO信息：

S1-1，下载GNSS卫星天线相位中心偏差文件；

S1-2，根据时间选择对应的相位中心偏差文件；

S1-3，相位中心偏差文件存储有各GNSS卫星不同频点上的天线相位中心偏差信息，存储方式依据导航系统、频点、三轴方向模式；

S2，通过解析法计算太阳在地心惯性系中的位置：

S2-1，计算TT时间：

式一：TT-TAI＝32.184；

式二：GPST-TAI＝-19.0；

TAI为世界原子时，单位秒；GPST是GPS时间，单位秒；通过式一和式二计算出已知GPST时间对应的TT时；

S2-2，计算太阳的位置：

将TT时转换成简化儒略日的形式，计算TT时下的旋转角θ，依据旋转角θ将全零的3×3矩阵对X轴进行旋转得到矩阵R，计算赤道面位置向量W，黄道坐标，并通过矩阵R将W旋转到赤道坐标上，进而得到太阳位置S0，如下式三：

式三：S0＝R·W；

S3，计算Z向量：

S3-1，采用拉格朗日插值算法获取当前时刻的GNSS导航卫星的位置P(x,y,z)向量；P向量是惯性坐标系ECI下的坐标向量，以地心为坐标原点；

S3-2，计算星固系到ECI的转换矩阵，转换矩阵在天球系统中进行计算；S3-3，通过Z＝0～P计算得到Z矩阵；

S3-4，将Z矩阵转化为单位向量，公式如下：

式四：

S4，计算S向量：

在S2-2和GNSS卫星位置已知的基础上，计算出太阳到GNSS卫星的位置向量S，通过两个向量作差的方式得到S向量，公式如下：

式五：S＝S0-P；

根据式四将S矩阵转化为单位向量；

S5，计算Y向量：

用Y向量表示Z向量和S向量组成平面的垂向量，利用向量叉乘计算可得Y向量，公式如下：

式六：Y＝Z×S；

根据式四将Y矩阵转化为单位向量；

S6，计算X向量：

用X向量表示Y向量和Z向量组成平面的垂向量，利用向量叉乘计算可得X向量，公式如下：

式七：X＝Y×Z；

根据式四将X矩阵转化为单位向量；

S7，计算旋转矩阵S2C：

用S2C表示星固系到ECI下的转换矩阵，公式如下：

式八：

S8，读取对应卫星的PCO信息：

根据当前GNSS导航卫星PRN号查找处对应的卫星两个频点上的PCO，分别用pco1(x,y,z)和pco2(x,y,z)表示；

S9，通过S2C将PCO旋转到ECI坐标系下：

通过S2C矩阵分别将pco1(x,y,z)和pco2(x,y,z)旋转到ECI坐标系下，公式如下：

C1[0]＝S2C[0]*pco1[0]+S2C[1]*pco1[1]+S2C[2]*pco1[2]

C1[1]＝S2C[3]*pco1[0]+S2C[4]*pco1[1]+S2C[5]*pco1[2]

式九：C1[2]＝S2C[6]*pco1[0]+S2C[7]*pco1[1]+S2C[8]*pco1[2]；

C2[0]＝S2C[0]*pco2[0]+S2C[1]*pco2[1]+S2C[2]*pco2[2]

C2[1]＝S2C[3]*pco2[0]+S2C[4]*pco2[1]+S2C[5]*pco2[2]

式十：C2[2]＝S2C[6]*pco2[0]+S2C[7]*pco2[1]+S2C[8]*pco2[2]；

C1和C2分别为旋转后两个频点上的PCO；

S10，双频组合计算GNSS卫星的组合PCO：

采用双频组合的方式计算出一个天线相位中心偏差C(x,y,z)代表该GNSS导航星接收信号的点，公式如下：

式十一：

C[0]＝(Freq_1*Freq_1*C1[0]-Freq_2*Freq_2*C2[0])/(Freq_1*Freq_1-Freq_2*Freq_2)

C[1]＝(Freq_1*Freq_1*C1[1]-Freq_2*Freq_2*C2[1])/(Freq_1*Freq_1-Freq_2*Freq_2)

C[2]＝(Freq_1*Freq_1*C1[2]-Freq_2*Freq_2*C2[2])/(Freq_1*Freq_1-Freq_2*Freq_2)；

Freq_1和Freq_2分别表示GNSS导航星两个频点上的频率，L1和L2分别为GPS卫星的两个频率，B1和B3分别为BDS卫星两个频率，E1和E5分别为GALILEO卫星两个频率，G1和G2分别为GLONASS卫星两个频率；

S11，将组合后的PCO改正到GNSS卫星位置上：

通过向量相加的方式将C改正到GNSS导航卫星位置P上得到改正天线相位中心偏差后的位置P2(x,y,z)，公式如下：

P2[0]＝P[0]+C[0]

P2[1]＝P[1]+C[1]

P2[2]＝P[2]+C[2]。

本发明同现有技术相比，具备以下有益效果：

本发明提高了计算GNSS卫星的精度，在GNSS掩星探测、GNSS反射探测以及其他低轨卫星载荷等遥感探测方面对GNSS卫星位置有较高精度要求的领域有重要应用价值。

附图说明

图1为本发明方法的技术路线示意图。

具体实施方式

现结合附图对本发明做进一步描述。

参见图1，一种基于太阳位置和双频组合的GNSS导航卫星PCO改正方法：

1.读取GNSS卫星PCO信息：

从外网下载GNSS卫星天线相位中心偏差文件igs14.atx，依据时间不同选择对应的文件，2022年年积日344日以前选择igs14.atx，之后选择igs20.atx。该文件中存储有各GNSS卫星不同频点上的天线相位中心偏差信息，依据导航系统、频点、三轴方向模式存储该信息。

2.计算太阳的位置：

在GNSS导航卫星PCO改正过程中用到了太阳的位置，需要对太阳的位置进行计算，本发明使用解析法计算太阳在地心惯性系中的位置。

太阳位置是在地心动力学时TT时间系统下进行的计算，通过以下公式：

TT-TAI＝32.184

GPST-TAI＝-19.0

计算出TT时间，其中，TAI为世界原子时，单位是秒，GPST是GPS时间，单位是秒，通过两个公式可以计算出已知GPST时间对应的TT时。

将TT时转换成简化儒略日的形式。计算TT时下的旋转角θ，依据该旋转角θ将全零的3×3矩阵对X轴进行旋转得到矩阵R，计算赤道面位置向量W，黄道坐标，并通过矩阵R将W旋转到赤道坐标上，进而得到太阳位置S0，如以下公式所示：

S0＝R·W

3.计算Z向量：

采用拉格朗日插值算法获取当前时刻的GNSS导航卫星的位置P(x,y,z)向量(下面简称P)，P向量是惯性坐标系ECI下的坐标向量，以地心为坐标原点。GNSS卫星的PCO信息均在星固系下定义的，GNSS导航星位置是ECI坐标系下，故需要计算星固系到ECI的转换矩阵，该转换矩阵在天球系统中进行计算，故通过Z＝0-P计算得到Z矩阵。由以下公式将Z矩阵转化为单位向量：

4.计算S向量：

太阳到GNSS卫星的位置用向量S表示。在上述太阳位置和GNSS卫星位置已知的基础上，计算出太阳到GNSS卫星的位置向量S，直接由两个向量作差的方式得到S向量：S＝S0-P。并由3中的公式将S矩阵转化为单位向量。

5.计算Y向量：

用Y向量表示Z向量和S向量组成平面的垂向量，利用向量叉乘计算可得Y向量：Y＝Z×S。并由3中的公式将Y矩阵转化为单位向量。

6.计算X向量：

用X向量表示Y向量和Z向量组成平面的垂向量，利用向量叉乘计算可得X向量：X＝Y×Z。并由3中的公式将X矩阵转化为单位向量。

7.计算旋转矩阵S2C：

GNSS导航卫星的天线PCO需要从星固系下旋转到ECI惯性系下才能改正到GNSS导航卫星的位置上，用S2C表示星固系到ECI下的转换矩阵，计算如以下公式所示。

8.读取对应卫星的PCO信息：

GNSS卫星天线PCO依据导航系统、频点、三轴方向模式存储,根据当前GNSS导航卫星PRN号查找处对应的卫星两个频点上的PCO分别用pco1(x,y,z)和pco2(x,y,z)表示；

9.通过S2C将PCO旋转到ECI坐标系下：

通过S2C矩阵分别将pco1(x,y,z)和pco2(x,y,z)旋转到ECI坐标系下，具体如以下公式所示，其中C1和C2分别为旋转后两个频点上的PCO；

C1[0]＝S2C[0]*pco1[0]+S2C[1]*pco1[1]+S2C[2]*pco1[2]

C1[1]＝S2C[3]*pco1[0]+S2C[4]*pco1[1]+S2C[5]*pco1[2]

C1[2]＝S2C[6]*pco1[0]+S2C[7]*pco1[1]+S2C[8]*pco1[2]

C2[0]＝S2C[0]*pco2[0]+S2C[1]*pco2[1]+S2C[2]*pco2[2]

C2[1]＝S2C[3]*pco2[0]+S2C[4]*pco2[1]+S2C[5]*pco2[2]

C2[2]＝S2C[6]*pco2[0]+S2C[7]*pco2[1]+S2C[8]*pco2[2]

10.双频组合计算GNSS卫星的组合PCO：

对GNSS遥感技术处理中通常使用一个GNSS位置，但是i)中计算出两个频点上的旋转后的C1和C2，采用双频组合的方式计算出一个天线相位中心偏差C(x,y,z)代表该GNSS导航星接收信号的点，具体计算如以下公式所示，其中Freq_1和Freq_2分别表示GNSS导航星两个频点上的频率，对于GPS卫星这两个频率分别为L1和L2，对于BDS卫星这两个频率分别为B1和B3，对于GALILEO卫星这两个频率分别为E1和E5，对于GLONASS卫星这两个频率分别为G1和G2。

C[0]＝(Freq_1*Freq_1*C1[0]-Freq_2*Freq_2*C2[0])/(Freq_1*Freq_1-Freq_2*Freq_2)

C[1]＝(Freq_1*Freq_1*C1[1]-Freq_2*Freq_2*C2[1])/(Freq_1*Freq_1-Freq_2*Freq_2)

C[2]＝(Freq_1*Freq_1*C1[2]-Freq_2*Freq_2*C2[2])/(Freq_1*Freq_1-Freq_2*Freq_2)；

11.将组合后的PCO改正到GNSS卫星位置上：

通过向量相加的方式将C改正到GNSS导航卫星位置P上得到改正天线相位中心偏差后的位置P2(x,y,z)，如公式所示：

P2[0]＝P[0]+C[0]

P2[1]＝P[1]+C[1]

P2[2]＝P[2]+C[2]。

以上仅是本发明的优选实施方式，只是用于帮助理解本申请的方法及其核心思想，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

本发明从整体上解决了现有技术中由于全球GNSS服务披露的轨道为卫星质心处的位置，与接收天线处位置偏差导致的计算偏差，进一步影响到对于掩星事件反演的误差过大的问题。通过引入一种基于太阳位置和双频组合的GNSS导航卫星PCO改正方法，提高了计算GNSS卫星的精度，在GNSS掩星探测、GNSS反射探测以及其他低轨卫星载荷等遥感探测方面对GNSS卫星位置有较高精度要求的领域有重要应用价值。