基于GNSS绝对定位数据的卫星姿态导引角计算方法及系统与流程

基于gnss绝对定位数据的卫星姿态导引角计算方法及系统
技术领域
1.本发明涉及卫星姿态导引规划与计算领域，具体地，涉及一种基于gnss绝对定位数据的卫星姿态导引角计算方法及系统。


背景技术：

2.卫星在轨飞行受到地球自转和椭圆轨道的影响，卫星成像时会出现图像位置偏移，图像变形等问题。特别是对sar卫星来说，多普勒中心会由于地球曲率及自转的存在，导致偏移，典型偏移值高于系统脉冲重复频率，影响图像质量。因此需要对卫星进行二维姿态导引角计算以降低影响，提高成像质量、简化后期处理过程。
3.槐超，王文妍.insar编队卫星全零多普勒姿态导引研究[j].上海航天,2014，31(6).中推导了一种insar编队卫星全零多普勒姿态导引方法，由辅星多普勒中心频率，给出了偏航角和俯仰角的解析表达式。
[0004]
郭德明,徐华平,李景文.高分辨率星载斜视sar的姿态导引[j].宇航学报,2011(5).中根据星地空间几何关系和姿态导引原理，推导了实现全零多普勒导引的偏航角和俯仰角的解析表达式，并进一步给出了最小化多普勒中心频率随距离变化的姿态导引规律，并进行了仿真。
[0005]
但是上述文献中均未涉及工程实际中gnss绝对定位数据的处理与引入。
[0006]
专利cn106843249a“一种二维导引姿态控制方法”中介绍了本发明公开了一种高精度高稳定度二维导引姿态控制方法，可以实现卫星任意位置的快速二维导引接入控制，并且能够提高系统的控制精度。未涉及导引角计算。
[0007]
专利文献cn106915477b(申请号：cn201710128823.6)公开了一种姿态控制方法，包含以下过程：采用秒脉冲信号对影响卫星姿态精度的星上时间进行对准；实时对卫星姿态测量用的星敏感器误差进行修正；利用动力学耦合关系对斜飞卫星的重力梯度干扰力矩进行补偿；在姿态控制过程中利用输入成型控制方法对卫星挠性进行抑制；采用位置速度双回路和位置修正回路增加角速度前馈指令的姿态控制律实现高精度高稳定度姿态导引控制；采用饱和滑模变结构控制算法实现大惯量卫星的快速姿态机动。但该发明也未涉及工程实际中gnss绝对定位数据的处理与引入。


技术实现要素：

[0008]
针对现有技术中的缺陷，本发明的目的是提供一种基于gnss绝对定位数据的卫星姿态导引角计算方法及系统。
[0009]
根据本发明提供的一种基于gnss绝对定位数据的卫星姿态导引角计算方法，包括：
[0010]
步骤s1：根据gnss实时绝对定位数据中的时间、位置、速度转换到惯性坐标系下，得到惯性坐标系下的位置、速度；
[0011]
步骤s2：将gnss实时绝对定位的时间下位置速度数值积分递推至当前星时，得到
当前星时时刻的位置和速度；
[0012]
步骤s3：将当前时刻的位置速度信息转换得到对应时刻的轨道参数；
[0013]
步骤s4：根据当前星时时刻轨道参数，计算当前时刻偏航导引角和俯仰导引角。
[0014]
优选地，在所述步骤s1中：
[0015]
步骤s1.1：
[0016]
根据时间t0计算j2000地心惯性系到地固坐标系的坐标变换矩阵：
[0017]meci2ecf
＝ep
·
er
·
nr
·
pr
[0018]
其中，ep为极移矩阵，简化为3
×
3单位阵；er为地球自转矩阵，nr为章动矩阵，pr为岁差矩阵；
[0019]
步骤s1.2：
[0020]
将gnss绝对定位数据位置速度信息由wgs84坐标系转换至j2000.0惯性坐标系，得到位置r
eci0
、速度v
eci0
为：
[0021]reci0
＝(m
eci2ecf
)-1
·recf0
[0022]veci0
＝(m
eci2ecf
)-1vecf0
+ωe×recf0
[0023]
其中，r
ecfo
为t0时刻wgs84坐标系下的位置，v
ecfo
为t0时刻wgs84坐标系下的速度；上表t表示转置，上表
·
表示对时间求导。
[0024]
优选地，在所述步骤s2中：
[0025]
根据t0时刻位置速度，采用四阶龙格库塔公式，数值积分递推至当前星时t，得到t时刻的位置r
eci
、速度v
eci
。
[0026]
优选地，在所述步骤s3中：
[0027]
步骤s3.1：计算过程参数：
[0028]
位置标量r＝|r
eci
|，速度标量v＝|v
eci
|
[0029]
轨道坐标系到j2000.0惯性坐标系转换矩阵：
[0030][0031]
其中，y＝v
eci
×reci
，x＝r
eci
×
y，均为3
×
1列向量；|x|表示x的模；
[0032]
轨道能量常量：
[0033][0034]
其中，地球引力常数μ＝398600km3/s2；
[0035]
动量矩：
[0036]
h＝r
eci
×veci
[0037]
动量矩的模：
[0038]
h＝|h|
[0039]
半通径：
[0040]
p＝h2/μ
[0041]
拉普拉斯矢量：
[0042][0043]
其中，l
x
为、ly为、lz分别拉普拉斯矢量l在惯性坐标系中的三个坐标轴的分量；
[0044]
步骤s3.2：计算t时刻轨道根数：
[0045]
半长轴：
[0046][0047]
偏心率：
[0048][0049]
轨道倾角：
[0050]
i＝arccos(h
zi
/h)
[0051]
其中，h
zi
为动量矩h在z向分量；
[0052]
纬度幅角：
[0053]
u＝arctan2(-m
orb2eci33
,m
orb2eci31
)
[0054]
其中，m
orb2eci33
为矩阵m
orb2eci
第3行第3列的数值，m
orb2eci31
为矩阵m
orb2eci
第3行第1列的数值；
[0055]
升交点赤经：
[0056]
ω＝arctan2(-m
orb2eci12
,m
orb2eci22
)
[0057]
其中，m
orb2eci12
为矩阵m
orb2eci
第1行第2列的数值，m
orb2eci22
为矩阵m
orb2eci
第2行第2列的数值；
[0058]
近地点幅角：
[0059][0060]
真近点角：
[0061]
f＝u-ω。
[0062]
优选地，在所述步骤s4中：
[0063]
偏航导引角ψ为：
[0064][0065]
其中：ωe为地球自转角速度，a为半长轴，r为位置标量，u为纬度幅角，i为轨道倾角，e为偏心率，f为真近点角；
[0066][0067]
μ＝398600km3/s2；
[0068]
俯仰导引角θ为：
[0069][0070]
根据本发明提供的一种基于gnss绝对定位数据的卫星姿态导引角计算系统，包括：
[0071]
模块m1：根据gnss实时绝对定位数据中的时间、位置、速度转换到惯性坐标系下，得到惯性坐标系下的位置、速度；
[0072]
模块m2：将gnss实时绝对定位的时间下位置速度数值积分递推至当前星时，得到当前星时时刻的位置和速度；
[0073]
模块m3：将当前时刻的位置速度信息转换得到对应时刻的轨道参数；
[0074]
模块m4：根据当前星时时刻轨道参数，计算当前时刻偏航导引角和俯仰导引角。
[0075]
优选地，在所述模块m1中：
[0076]
模块m1.1：
[0077]
根据时间t0计算j2000地心惯性系到地固坐标系的坐标变换矩阵：
[0078]meci2ecf
＝ep
·
er
·
nr
·
pr
[0079]
其中，ep为极移矩阵，简化为3
×
3单位阵；er为地球自转矩阵，nr为章动矩阵，pr为岁差矩阵；
[0080]
模块m1.2：
[0081]
将gnss绝对定位数据位置速度信息由wgs84坐标系转换至j2000.0惯性坐标系，得到位置r
eci0
、速度v
eci0
为：
[0082]reci0
＝(m
eci2ecf
)-1
·recf0
[0083]veci0
＝(m
eci2ecf
)-1vecf0
+ωe×recf0
[0084]
其中，r
ecfo
为t0时刻wgs84坐标系下的位置，v
ecfo
为t0时刻wgs84坐标系下的速度；上表t表示转置，上表
·
表示对时间求导。
[0085]
优选地，在所述模块m2中：
[0086]
根据t0时刻位置速度，采用四阶龙格库塔公式，数值积分递推至当前星时t，得到t时刻的位置r
eci
、速度v
eci
。
[0087]
优选地，在所述模块m3中：
[0088]
模块m3.1：计算过程参数：
[0089]
位置标量r＝|r
eci
|，速度标量v＝|v
eci
|
[0090]
轨道坐标系到j2000.0惯性坐标系转换矩阵：
[0091][0092]
其中，y＝v
eci
×reci
，x＝r
eci
×
y，均为3
×
1列向量；|x|表示x的模；
[0093]
轨道能量常量：
[0094][0095]
其中，地球引力常数μ＝398600km3/s2；
[0096]
动量矩：
[0097]
h＝r
eci
×veci
[0098]
动量矩的模：
[0099]
h＝|h|
[0100]
半通径：
[0101]
p＝h2/μ
[0102]
拉普拉斯矢量：
[0103][0104]
其中，l
x
为、ly为、lz分别拉普拉斯矢量l在惯性坐标系中的三个坐标轴的分量；
[0105]
模块m3.2：计算t时刻轨道根数：
[0106]
半长轴：
[0107][0108]
偏心率：
[0109][0110]
轨道倾角：
[0111]
i＝arccos(h
zi
/h)
[0112]
其中，h
zi
为动量矩h在z向分量；
[0113]
纬度幅角：
[0114]
u＝arctan2(-m
orb2eci33
,m
orb2eci31
)
[0115]
其中，m
orb2eci33
为矩阵m
orb2eci
第3行第3列的数值，m
orb2eci31
为矩阵m
orb2eci
第3行第1列的数值；
[0116]
升交点赤经：
[0117]
ω＝arctan2(-m
orb2eci12
,m
orb2eci22
)
[0118]
其中，m
orb2eci12
为矩阵m
orb2eci
第1行第2列的数值，m
orb2eci22
为矩阵m
orb2eci
第2行第2列的数值；
[0119]
近地点幅角：
[0120][0121]
真近点角：
[0122]
f＝u-ω。
[0123]
优选地，在所述模块m4中：
[0124]
偏航导引角ψ为：
[0125]
[0126]
其中：ωe为地球自转角速度，a为半长轴，r为位置标量，u为纬度幅角，i为轨道倾角，e为偏心率，f为真近点角；
[0127][0128]
μ＝398600km3/s2；
[0129]
俯仰导引角θ为：
[0130][0131]
与现有技术相比，本发明具有如下的有益效果：
[0132]
1、本计算方法实效性高；
[0133]
2、本发明结合工程实际，可以高精度实时计算卫星姿态导引角；
[0134]
3、本方法可以实时高精度计算卫星姿态导引角，有效支持高精度的姿态控制，有效补偿地球自转引起的图像形变，提高成像质量。
附图说明
[0135]
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
[0136]
图1为本发明基于gnss绝对定位数据的卫星姿态导引角计算流程图；
[0137]
图2为根据本发明方法计算某型号卫星得到的姿态导引角，其中滚动方向导引角为0。
具体实施方式
[0138]
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
[0139]
实施例1：
[0140]
卫星在轨依靠gnss系统实时提供一定频率更新的高精度的绝对定位数据，本发明利用该数据实时高精度计算当前时刻的姿态导引角。有效支持高精度的姿态控制，进而提高成像质量。
[0141]
根据本发明提供的一种基于gnss绝对定位数据的卫星姿态导引角计算方法，如图1-图2所示，包括：
[0142]
步骤s1：根据gnss实时绝对定位数据中的时间、位置、速度转换到惯性坐标系下，得到惯性坐标系下的位置、速度；
[0143]
具体地，在所述步骤s1中：
[0144]
步骤s1.1：
[0145]
根据时间t0计算j2000地心惯性系到地固坐标系的坐标变换矩阵：
[0146]meci2ecf
＝ep
·
er
·
nr
·
pr
[0147]
其中，ep为极移矩阵，简化为3
×
3单位阵；er为地球自转矩阵，nr为章动矩阵，pr为岁差矩阵；
[0148]
步骤s1.2：
[0149]
将gnss绝对定位数据位置速度信息由wgs84坐标系转换至j2000.0惯性坐标系，得到位置r
eci0
、速度v
eci0
为：
[0150]reci0
＝(m
eci2ecf
)-1
·recf0
[0151]veci0
＝(m
eci2ecf
)-1vecf0
+ωe×recf0
[0152]
其中，r
ecfo
为t0时刻wgs84坐标系下的位置，v
ecfo
为t0时刻wgs84坐标系下的速度；上表t表示转置，上表
·
表示对时间求导。
[0153]
步骤s2：将gnss实时绝对定位的时间下位置速度数值积分递推至当前星时，得到当前星时时刻的位置和速度；
[0154]
具体地，在所述步骤s2中：
[0155]
根据t0时刻位置速度，采用四阶龙格库塔公式，数值积分递推至当前星时t，得到t时刻的位置r
eci
、速度v
eci
。
[0156]
步骤s3：将当前时刻的位置速度信息转换得到对应时刻的轨道参数；
[0157]
具体地，在所述步骤s3中：
[0158]
步骤s3.1：计算过程参数：
[0159]
位置标量r＝|r
eci
|，速度标量v＝|v
eci
|
[0160]
轨道坐标系到j2000.0惯性坐标系转换矩阵：
[0161][0162]
其中，y＝v
eci
×reci
，x＝r
eci
×
y，均为3
×
1列向量；|x|表示x的模；
[0163]
轨道能量常量：
[0164][0165]
其中，地球引力常数μ＝398600km3/s2；
[0166]
动量矩：
[0167]
h＝r
eci
×veci
[0168]
动量矩的模：
[0169]
h＝|h|
[0170]
半通径：
[0171]
p＝h2/μ
[0172]
拉普拉斯矢量：
[0173][0174]
其中，l
x
为、ly为、lz分别拉普拉斯矢量l在惯性坐标系中的三个坐标轴的分量；
[0175]
步骤s3.2：计算t时刻轨道根数：
[0176]
半长轴：
[0177][0178]
偏心率：
[0179][0180]
轨道倾角：
[0181]
i＝arccos(h
zi
/h)
[0182]
其中，h
zi
为动量矩h在z向分量；
[0183]
纬度幅角：
[0184]
u＝arctan2(-m
orb2eci33
,m
orb2eci31
)
[0185]
其中，m
orb2eci33
为矩阵m
orb2eci
第3行第3列的数值，m
orb2eci31
为矩阵m
orb2eci
第3行第1列的数值；
[0186]
升交点赤经：
[0187]
ω＝arctan2(-m
orb2eci12
,m
orb2eci22
)
[0188]
其中，m
orb2eci12
为矩阵m
orb2eci
第1行第2列的数值，m
orb2eci22
为矩阵m
orb2eci
第2行第2列的数值；
[0189]
近地点幅角：
[0190][0191]
真近点角：
[0192]
f＝u-ω。
[0193]
步骤s4：根据当前星时时刻轨道参数，计算当前时刻偏航导引角和俯仰导引角。
[0194]
具体地，在所述步骤s4中：
[0195]
偏航导引角ψ为：
[0196][0197]
其中：ωe为地球自转角速度，a为半长轴，r为位置标量，u为纬度幅角，i为轨道倾角，e为偏心率，f为真近点角；
[0198][0199]
μ＝398600km3/s2；
[0200]
俯仰导引角θ为：
[0201][0202]
实施例2：
[0203]
实施例2为实施例1的优选例，以更为具体地对本发明进行说明。
[0204]
本领域技术人员可以将本发明提供的一种基于gnss绝对定位数据的卫星姿态导引角计算方法，理解为基于gnss绝对定位数据的卫星姿态导引角计算系统的具体实施方式，即所述基于gnss绝对定位数据的卫星姿态导引角计算系统可以通过执行所述基于gnss绝对定位数据的卫星姿态导引角计算方法的步骤流程予以实现。
[0205]
根据本发明提供的一种基于gnss绝对定位数据的卫星姿态导引角计算系统，包括：
[0206]
模块m1：根据gnss实时绝对定位数据中的时间、位置、速度转换到惯性坐标系下，得到惯性坐标系下的位置、速度；
[0207]
具体地，在所述模块m1中：
[0208]
模块m1.1：
[0209]
根据时间t0计算j2000地心惯性系到地固坐标系的坐标变换矩阵：
[0210]meci2ecf
＝ep
·
er
·
ar
·
pr
[0211]
其中，ep为极移矩阵，简化为3
×
3单位阵；er为地球自转矩阵，nr为章动矩阵，pr为岁差矩阵；
[0212]
模块m1.2：
[0213]
将gnss绝对定位数据位置速度信息由wgs84坐标系转换至j2000.0惯性坐标系，得到位置r
eci0
、速度v
eci0
为：
[0214]reci0
＝(m
eci2ecf
)-1
·recf0
[0215]veci0
＝(m
eci2ecf
)-1vecf0
+ωe×recf0
[0216]
其中，r
ecfo
为t0时刻wgs84坐标系下的位置，v
ecfo
为t0时刻wgs84坐标系下的速度；上表t表示转置，上表
·
表示对时间求导。
[0217]
模块m2：将gnss实时绝对定位的时间下位置速度数值积分递推至当前星时，得到当前星时时刻的位置和速度；
[0218]
具体地，在所述模块m2中：
[0219]
根据t0时刻位置速度，采用四阶龙格库塔公式，数值积分递推至当前星时t，得到t时刻的位置r
eci
、速度v
eci
。
[0220]
模块m3：将当前时刻的位置速度信息转换得到对应时刻的轨道参数；
[0221]
具体地，在所述模块m3中：
[0222]
模块m3.1：计算过程参数：
[0223]
位置标量r＝|r
eci
|，速度标量v＝|v
eci
|
[0224]
轨道坐标系到j2000.0惯性坐标系转换矩阵：
[0225][0226]
其中，y＝v
eci
×reci
，x＝r
eci
×
y，均为3
×
1列向量；|x|表示x的模；
[0227]
轨道能量常量：
[0228][0229]
其中，地球引力常数μ＝398600km3/s2；
[0230]
动量矩：
[0231]
h＝r
eci
×veci
[0232]
动量矩的模：
[0233]
h＝|h|
[0234]
半通径：
[0235]
p＝h2/μ
[0236]
拉普拉斯矢量：
[0237][0238]
其中，l
x
为、ly为、lz分别拉普拉斯矢量l在惯性坐标系中的三个坐标轴的分量；
[0239]
模块m3.2：计算t时刻轨道根数：
[0240]
半长轴：
[0241][0242]
偏心率：
[0243][0244]
轨道倾角：
[0245]
i＝arccos(h
zi
/h)
[0246]
其中，h
zi
为动量矩h在z向分量；
[0247]
纬度幅角：
[0248]
u＝arctan2(-m
orb2eci33
,m
orb2eci31
)
[0249]
其中，m
orb2eci33
为矩阵m
orb2eci
第3行第3列的数值，m
orb2eci31
为矩阵m
orb2eci
第3行第1列的数值；
[0250]
升交点赤经：
[0251]
ω＝arctan2(-m
orb2eci12
,m
orb2eci22
)
[0252]
其中，m
orb2eci12
为矩阵m
orb2eci
第1行第2列的数值，m
orb2eci22
为矩阵m
orb2eci
第2行第2列的数值；
[0253]
近地点幅角：
[0254][0255]
真近点角：
[0256]
f＝u-ω。
[0257]
模块m4：根据当前星时时刻轨道参数，计算当前时刻偏航导引角和俯仰导引角。
[0258]
具体地，在所述模块m4中：
[0259]
偏航导引角ψ为：
[0260][0261]
其中：ωe为地球自转角速度，a为半长轴，r为位置标量，u为纬度幅角，i为轨道倾角，e为偏心率，f为真近点角；
[0262][0263]
μ＝398600km3/s2；
[0264]
俯仰导引角θ为：
[0265][0266]
实施例3：
[0267]
实施例3为实施例1的优选例，以更为具体地对本发明进行说明。
[0268]
针对卫星在轨实时计算二维导引角以提高卫星载荷成像质量的需求，本发明公开了一种基于gnss绝对定位数据计算卫星姿态导引角的方法。通过对gnss实时输出的绝对定位数据进行转换处理，得到卫星当前时刻轨道根数，计算卫星偏航和俯仰导引角。该方面结合工程实际，可以高精度实时计算卫星姿态导引角。为了实现上述目的，本发明通过以下技术方案实现，具体包括以下环节。
[0269]
一种基于gnss绝对定位数据的卫星姿态导引角计算方法，包括以下步骤：
[0270]
步骤1：根据gnss实时绝对定位数据中的时间t0、位置r
ecf0
、速度v
ecf0
转换到惯性坐标系下描述得到位置r
eci0
、速度v
eci0
。
[0271]
步骤2：将t0时刻位置速度，数值积分递推至当前星时t，得到t时刻的位置r
eci
、速度v
eci
。
[0272]
步骤3：将当前时刻的位置速度信息转换得到对应时刻的轨道参数。
[0273]
步骤4：根据t时刻轨道参数，分别计算当前时刻偏航导引角和俯仰导引角。
[0274]
所述的步骤1包括以下内容：
[0275]
步骤1-1：根据时间t0计算j2000地心惯性系到地固坐标系的坐标变换矩阵m
eci2ecf
＝ep
·
er
·
nr
·
pr，ep为极移矩阵，简化为3
×
3单位阵；er为地球自转矩阵，nr为章动矩阵，pr为岁差矩阵。
[0276]
步骤1-2：将gnss绝对定位数据位置速度信息由wgs84坐标系转换至j2000.0惯性坐标系，得到位置r
eci0
、速度v
eci0
为
[0277]reci0
＝(m
eci2ecf
)-1
·recf0
[0278]veci0
＝(m
eci2ecf
)-1vecf0
+ωe×recf0
[0279]
其中，上表t表示转置，上表
·
表示对时间求导。
[0280]
根据t0时刻位置速度，采用四阶龙格－库塔公式，数值积分递推至当前星时t，得到t时刻的位置r
eci
、速度v
eci
。
[0281]
所述的步骤3包括以下内容
[0282]
步骤3-1：计算过程参数
[0283]
位置标量r＝|r
eci
|，速度标量v＝|v
eci
|
[0284]
轨道坐标系到j2000.0惯性坐标系转换矩阵其中，y＝v
eci
×reci
，x＝r
eci
×
y，均为3
×
1列向量。，|x|表示x的模。
[0285]
轨道能量常量动量矩h＝r
eci
×veci
，动量矩的模h＝|h|，半通径p＝h2/μ，拉普拉斯矢量地球引力常数μ＝398600km3/s2。
[0286]
步骤3-2：计算t时刻轨道根数
[0287]
半长轴
[0288]
偏心率
[0289]
轨道倾角i＝arccos(h
zi
/h)，h
zi
为动量矩h在z向分量。
[0290]
纬度幅角u＝arctan2(-m
orb2eci33
,m
orb2eci31
)。
[0291]
升交点赤经ω＝arctan2(-m
orb2eci12
,m
orb2eci22
)。
[0292]
近地点幅角
[0293]
真近点角f＝u-ω
[0294]
偏航导引角ψ为：
[0295][0296]
其中：ωe为地球自转角速度
[0297]
俯仰导引角θ为：
[0298][0299]
本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统、装置及其各个模块以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统、装置及其各个模块以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同程序。所以，本发明提供的系统、装置及其各个模块可以被认为是一种硬件部件，而对其内包括的用于实现各种程序的模块也可以视为硬件部件内的结构；也可以将用于实现各种功能的模块视为既可以是实现方法的软件程序又可以是硬件部件内的结构。
[0300]
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本技术的实施例和实施例中的特征可以任意相互组合。

技术特征：
1.一种基于gnss绝对定位数据的卫星姿态导引角计算方法，，其特征在于，包括：步骤s1：根据gnss实时绝对定位数据中的时间、位置、速度转换到惯性坐标系下，得到惯性坐标系下的位置、速度；步骤s2：将gnss实时绝对定位的时间下位置速度数值积分递推至当前星时，得到当前星时时刻的位置和速度；步骤s3：将当前时刻的位置速度信息转换得到对应时刻的轨道参数；步骤s4：根据当前星时时刻轨道参数，计算当前时刻偏航导引角和俯仰导引角。2.根据权利要求1所述的基于gnss绝对定位数据的卫星姿态导引角计算方法，其特征在于，在所述步骤s1中：步骤s1.1：根据时间t0计算j2000地心惯性系到地固坐标系的坐标变换矩阵：m
eci2ecf
＝ep
·
er
·
nr
·
pr其中，ep为极移矩阵，简化为3
×
3单位阵；er为地球自转矩阵，nr为章动矩阵，pr为岁差矩阵；步骤s1.2：将gnss绝对定位数据位置速度信息由wgs84坐标系转换至j2000.0惯性坐标系，得到位置r
eci0
、速度v
eci0
为：r
eci0
＝(m
eci2ecf
)-1
·
r
ecf0
v
eci0
＝(m
eci2ecf
)-1
v
ecf0
+ω
e
×
r
ecf0
其中，r
ecfo
为t0时刻wgs84坐标系下的位置，v
ecfo
为t0时刻wgs84坐标系下的速度；上表t表示转置，上表
·
表示对时间求导。3.根据权利要求1所述的基于gnss绝对定位数据的卫星姿态导引角计算方法，其特征在于，在所述步骤s2中：根据t0时刻位置速度，采用四阶龙格库塔公式，数值积分递推至当前星时t，得到t时刻的位置r
eci
、速度v
eci
。4.根据权利要求1所述的基于gnss绝对定位数据的卫星姿态导引角计算方法，其特征在于，在所述步骤s3中：步骤s3.1：计算过程参数：位置标量r＝r
eci
，速度标量v＝v
eci
轨道坐标系到j2000.0惯性坐标系转换矩阵：其中，y＝v
eci
×
r
eci
，x＝r
eci
×
y，均为3
×
1列向量；x表示x的模；轨道能量常量：其中，地球引力常数μ＝398600km3/s2；动量矩：
h＝r
eci
×
v
eci
动量矩的模：h＝h半通径：p＝h2μ拉普拉斯矢量：其中，l
x
为、l
y
为、l
z
分别拉普拉斯矢量l在惯性坐标系中的三个坐标轴的分量；步骤s3.2：计算t时刻轨道根数：半长轴：偏心率：轨道倾角：i＝arccos(h
zi
/h)其中，h
zi
为动量矩h在z向分量；纬度幅角：u＝arctan2(-m
orb2eci33
,m
orb2eci31
)其中，m
orb2eci33
为矩阵m
orb2eci
第3行第3列的数值，m
orb2eci31
为矩阵m
orb2eci
第3行第1列的数值；升交点赤经：ω＝arctan2(-m
orb2eci12
,m
orb2eci22
)其中，m
orb2eci12
为矩阵m
orb2eci
第1行第2列的数值，m
orb2eci22
为矩阵m
orb2eci
第2行第2列的数值；近地点幅角：真近点角：f＝u-ω。5.根据权利要求1所述的基于gnss绝对定位数据的卫星姿态导引角计算方法，其特征在于，在所述步骤s4中：偏航导引角ψ为：
其中：ω
e
为地球自转角速度，a为半长轴，r为位置标量，u为纬度幅角，i为轨道倾角，e为偏心率，f为真近点角；μ＝398600km3/s2；俯仰导引角θ为：6.一种基于gnss绝对定位数据的卫星姿态导引角计算系统，其特征在于，包括：模块m1：根据gnss实时绝对定位数据中的时间、位置、速度转换到惯性坐标系下，得到惯性坐标系下的位置、速度；模块m2：将gnss实时绝对定位的时间下位置速度数值积分递推至当前星时，得到当前星时时刻的位置和速度；模块m3：将当前时刻的位置速度信息转换得到对应时刻的轨道参数；模块m4：根据当前星时时刻轨道参数，计算当前时刻偏航导引角和俯仰导引角。7.根据权利要求6所述的基于gnss绝对定位数据的卫星姿态导引角计算系统，其特征在于，在所述模块m1中：模块m1.1：根据时间t0计算j2000地心惯性系到地固坐标系的坐标变换矩阵：m
eci2ecf
＝ep
·
er
·
nr
·
pr其中，ep为极移矩阵，简化为3
×
3单位阵；er为地球自转矩阵，nr为章动矩阵，pr为岁差矩阵；模块m1.2：将gnss绝对定位数据位置速度信息由wgs84坐标系转换至j2000.0惯性坐标系，得到位置r
eci0
、速度v
eci0
为：r
eci0
＝(m
eci2ecf
)-1
·
r
ecf0
v
eci0
＝(m
eci2ecf
)-1
v
ecf0
+ω
e
×
r
ecf0
其中，r
ecfo
为t0时刻wgs84坐标系下的位置，v
ecfo
为t0时刻wgs84坐标系下的速度；上表t表示转置，上表
·
表示对时间求导。8.根据权利要求6所述的基于gnss绝对定位数据的卫星姿态导引角计算系统，其特征在于，在所述模块m2中：根据t0时刻位置速度，采用四阶龙格库塔公式，数值积分递推至当前星时t，得到t时刻的位置r
eci
、速度v
eci
。9.根据权利要求6所述的基于gnss绝对定位数据的卫星姿态导引角计算系统，其特征在于，在所述模块m3中：模块m3.1：计算过程参数：位置标量r＝|r
eci
|，速度标量v＝|v
eci
|
轨道坐标系到j2000.0惯性坐标系转换矩阵：其中，y＝v
eci
×
r
eci
，x＝r
eci
×
y，均为3
×
1列向量；x表示x的模；轨道能量常量：其中，地球引力常数μ＝398600km3/s2；动量矩：h＝r
eci
×
v
eci
动量矩的模：h＝h半通径：p＝h2μ拉普拉斯矢量：其中，l
x
为、l
y
为、l
z
分别拉普拉斯矢量l在惯性坐标系中的三个坐标轴的分量；模块m3.2：计算t时刻轨道根数：半长轴：偏心率：轨道倾角：i＝arccos(h
zi
/h)其中，h
zi
为动量矩h在z向分量；纬度幅角：u＝arctan2(-m
orb2eci33
,m
orb2eci31
)其中，m
orb2eci33
为矩阵m
orb2eci
第3行第3列的数值，m
orb2eci31
为矩阵m
orb2eci
第3行第1列的数值；升交点赤经：ω＝arctan2(-m
orb2eci12
,m
orb2eci22
)其中，m
orb2eci12
为矩阵m
orb2eci
第1行第2列的数值，m
orb2eci22
为矩阵m
orb2eci
第2行第2列的数值；近地点幅角：
真近点角：f＝u-ω。10.根据权利要求6所述的基于gnss绝对定位数据的卫星姿态导引角计算系统，其特征在于，在所述模块m4中：偏航导引角ψ为：其中：ω
e
为地球自转角速度，a为半长轴，r为位置标量，u为纬度幅角，i为轨道倾角，e为偏心率，f为真近点角；μ＝398600km3/s2；俯仰导引角θ为：

技术总结
本发明提供了一种基于GNSS绝对定位数据的卫星姿态导引角计算方法及系统，包括：步骤S1：根据GNSS实时绝对定位数据中的时间、位置、速度转换到惯性坐标系下，得到惯性坐标系下的位置、速度；步骤S2：将GNSS实时绝对定位的时间下位置速度数值积分递推至当前星时，得到当前星时时刻的位置和速度；步骤S3：将当前时刻的位置速度信息转换得到对应时刻的轨道参数；步骤S4：根据当前星时时刻轨道参数，计算当前时刻偏航导引角和俯仰导引角。本方法可以实时高精度计算卫星姿态导引角，有效支持高精度的姿态控制，有效补偿地球自转引起的图像形变，提高成像质量。高成像质量。高成像质量。


技术研发人员：栗双岭 温俊健 王赟 李楠 孙英梅 韩东升 仲兆宇
受保护的技术使用者：上海卫星工程研究所
技术研发日：2022.11.24
技术公布日：2023/4/5