CN115829879B - 面向敏捷卫星的姿态四元数处理方法、装置、设备 - Google Patents
面向敏捷卫星的姿态四元数处理方法、装置、设备 Download PDFInfo
- Publication number
- CN115829879B CN115829879B CN202211620133.XA CN202211620133A CN115829879B CN 115829879 B CN115829879 B CN 115829879B CN 202211620133 A CN202211620133 A CN 202211620133A CN 115829879 B CN115829879 B CN 115829879B
- Authority
- CN
- China
- Prior art keywords
- coordinate system
- rotation matrix
- gesture
- satellite
- quaternion
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000003672 processing method Methods 0.000 title claims description 5
- 239000011159 matrix material Substances 0.000 claims abstract description 155
- 238000003384 imaging method Methods 0.000 claims abstract description 83
- 238000000034 method Methods 0.000 claims abstract description 59
- 238000006243 chemical reaction Methods 0.000 claims abstract description 39
- 239000013598 vector Substances 0.000 claims description 43
- 238000004364 calculation method Methods 0.000 claims description 32
- 238000005070 sampling Methods 0.000 claims description 15
- 238000004590 computer program Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 description 14
- 230000006870 function Effects 0.000 description 8
- 230000008859 change Effects 0.000 description 7
- 238000007781 pre-processing Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 6
- 230000009466 transformation Effects 0.000 description 6
- 238000001914 filtration Methods 0.000 description 5
- 238000005259 measurement Methods 0.000 description 5
- PCTMTFRHKVHKIS-BMFZQQSSSA-N (1s,3r,4e,6e,8e,10e,12e,14e,16e,18s,19r,20r,21s,25r,27r,30r,31r,33s,35r,37s,38r)-3-[(2r,3s,4s,5s,6r)-4-amino-3,5-dihydroxy-6-methyloxan-2-yl]oxy-19,25,27,30,31,33,35,37-octahydroxy-18,20,21-trimethyl-23-oxo-22,39-dioxabicyclo[33.3.1]nonatriaconta-4,6,8,10 Chemical compound C1C=C2C[C@@H](OS(O)(=O)=O)CC[C@]2(C)[C@@H]2[C@@H]1[C@@H]1CC[C@H]([C@H](C)CCCC(C)C)[C@@]1(C)CC2.O[C@H]1[C@@H](N)[C@H](O)[C@@H](C)O[C@H]1O[C@H]1/C=C/C=C/C=C/C=C/C=C/C=C/C=C/[C@H](C)[C@@H](O)[C@@H](C)[C@H](C)OC(=O)C[C@H](O)C[C@H](O)CC[C@@H](O)[C@H](O)C[C@H](O)C[C@](O)(C[C@H](O)[C@H]2C(O)=O)O[C@H]2C1 PCTMTFRHKVHKIS-BMFZQQSSSA-N 0.000 description 4
- 244000007853 Sarothamnus scoparius Species 0.000 description 4
- 230000001133 acceleration Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 230000003287 optical effect Effects 0.000 description 4
- 230000036962 time dependent Effects 0.000 description 4
- 230000033001 locomotion Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000009499 grossing Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 230000002093 peripheral effect Effects 0.000 description 2
- 230000003068 static effect Effects 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 238000005295 random walk Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
Landscapes
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本申请涉及一种面向敏捷卫星的姿态四元数处理方法、装置、设备。获取姿态四元数、星历数据；根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系对应的第一姿态旋转矩阵；根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角；根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角；利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数。
Description
技术领域
本申请涉及敏捷卫星技术领域，尤其是涉及一种面向敏捷卫星的姿态四元数处理方法、装置、设备。
背景技术
敏捷卫星与传统卫星相比具有俯仰、滚动、偏航等姿态动作，能够展现出更大的对地观测范围，拥有更强的对地观测能力。因此通过敏捷卫星可以更灵活的获取对地观测数据。目前卫星地面预处理系统对卫星下传的姿态四元数通常采用球面线性内插的方法进行插值处理。这种方式将两个已知的四元数所代表的旋转角以及中间待插值的四元数旋转看作位于四维空间的单位球上，将四元数旋转在转化为单位球上的轴-角旋转模型，内插过程中保持旋转轴不变，对旋转角进行和时间相关的线性内插，插值四元数与已知四元数之间的旋转偏差通过四元数乘法表达，插值完成后将采用“轴-角”表示的旋转变换回四元数形式。这是一种数学意义上的处理。
但是球面线性内插的方法只针对四元数代表的旋转本身，无法解决敏捷卫星所产生的噪声问题，也无法对姿态数据进行外延处理。以至于数据处理精度较低。因此，对于具备复杂成像模式的敏捷卫星来说，球面线性内插的方法并不是理想的处理方法。
发明内容
本申请提供一种面向敏捷卫星的姿态四元数处理方法、装置、设备。
第一方面，本申请提供一种面向敏捷卫星的姿态四元数处理方法，包括：
获取姿态四元数、星历数据；
根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系对应的第一姿态旋转矩阵；
根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角；
根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角；
利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数。
本申请提供了一种面向敏捷卫星的姿态四元数处理方法、装置、设备。可以获取星上计算机传输的姿态四元数以及星历数据。在获取成功后，利用坐标系转换规则，确定卫星的本体坐标系到当地轨道坐标系所对应的第一姿态旋转矩阵。进而根据预设的三轴旋转顺序，将第一姿态旋转矩阵中的参数进行拆分，得到对应的欧拉角，使得具有数学意义的姿态四元数具有物理意义。进而通过卫星的推扫成像方式，将得到的欧拉角进行拟合，得到模型化后的欧拉角。通过这种方式将最初得到的姿态四元数中的噪声进行平滑处理，或者遵守姿态角、角速度和角加速度随时间的变化规律，提高数据内部相对几何精度。从而利用坐标系转换规则，将模型化后的欧拉角对应的序列所构成的旋转矩阵转换成处理后的姿态四元数。提高最初得到的姿态四元数的精度。
第二方面，本申请提供一种面向敏捷卫星的姿态四元数处理装置。
数据获取模块，用于获取姿态四元数、星历数据；
第一姿态旋转矩阵确定模块，用于根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系对应的第一姿态旋转矩阵；
欧拉角确定模块，用于根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角；
欧拉角拟合模块，用于根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角；
处理模块，用于利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数。
第三方面，本申请提供一种面向敏捷卫星的姿态四元数处理设备，包括：存储器和处理器，所述存储器上存储有能够被处理器加载并执行第一方面的方法的计算机程序。
第四方面，本申请提供一种计算机可读存储介质，存储有能够被处理器加载并执行第一方面的方法的计算机程序。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本申请一实施例提供的一种应用场景示意图；
图2为本申请一实施例提供的一种面向敏捷卫星的姿态四元数处理方法的流程图；
图3为本申请一实施例提供的另一种面向敏捷卫星的姿态四元数处理方法的流程图；
图4a-4c为本申请一实施例提供的一种被动推扫方式下三轴的原始欧拉角与拟合欧拉角曲线对比图；
图5a-5c为本申请一实施例提供的一种主动推扫方式下三轴的原始欧拉角与不同多项式外推的曲线之间的拟合情况对比图；
图6a-6b为本申请一实施例提供的一种使用二次多项式拟合的主动推扫与被动推扫姿态四元数的物方残差示意图；
图7a-7b为本申请一实施例提供的一种使用五次多项式拟合的主动推扫与被动推扫姿态四元的物方残差示意图；
图8为本申请一实施例提供的一种面向敏捷卫星的姿态四元数处理装置的结构示意图；
图9为本申请一实施例提供的一种面向敏捷卫星的姿态四元数处理设备的结构示意图。
具体实施方式
为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
另外，本文中术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况。另外，本文中字符“/”，如无特殊说明，一般表示前后关联对象是一种“或”的关系。
下面结合说明书附图对本申请实施例作进一步详细描述。
由于目前卫星地面预处理系统对卫星下传的姿态四元数通常采用球面线性内插的方法进行插值处理，来提高姿态四元数的精度。这种方式属于数学意义上的处理。由于敏捷卫星与传统卫星相比具有俯仰、滚动、偏航等姿态动作，因此在获取数据的过程中会产生噪声。但这种情况，通过球面线性内插的方法无法解决，因此使得数据处理的精度较低，并不能达到理想的精度。
基于此，本申请提供一种面向敏捷卫星的姿态四元数处理方法、装置、电子设备。在获取星上计算机系统下传的姿态四元数后，卫星地面预处理系统将姿态四元数根据坐标转换规则进行相应地转换后，得到对应的欧拉角，进而通过卫星的推扫成像模式，选择不同的多项式进行拟合得到模型化后的欧拉角，使得数据内部的相对精度以及绝对精度都能得到提高。进而将模型化后欧拉角重新利用坐标系转换规则转换回姿态四元数，完成姿态四元数的处理，提高数据的精度。
图1为本申请提供的一种应用场景示意图。在获取星上计算机系统下传的姿态四元数后，卫星地面预处理系统利用坐标系转换规则以及相应的计算规则，结合全球导航卫星系统进行相应的操作后，使得姿态四元数内部的相对精度以及绝对精度都能得到提高。
具体的实现方式可以参考以下实施例。
图2为本申请一实施例提供的一种数据传输方法的流程图，本实施例的方法可以应用于以上场景中的卫星地面预处理系统。如图2所示的，该方法包括：
S201、获取姿态四元数、星历数据。
姿态四元数是星上计算机系统对星敏感器感器获取的四元数以及陀螺角速度的数据进行联合滤波后得到的。姿态四元数是卫星本体坐标系到地心惯性坐标系的协议坐标系J2000.0的四元数。
星历数据是通过全球导航卫星系统时刻捕捉到的卫星在地心固定坐标系下的实时的位置和速度矢量。
具体的，卫星地面预处理系统通过星上计算机系统获取进行联合滤波后得到姿态四元数。并获取来自全球导航卫星系统时刻捕捉到的卫星在地心固定坐标系下的实时的位置和速度矢量。
S202、根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系对应的第一姿态旋转矩阵。
坐标系转换规则是根据不同坐标系之间的关系所设定的。比如，由于姿态四元数是卫星本体坐标系到地心惯性坐标系的协议坐标系J2000.0的四元数。所以本体坐标系到地心惯性坐标系的之间的坐标系转换可以通过姿态四元数对应分量的各个部分进行处理。又或者，由于瞬时坐标系是时刻变化的，所以地心固定坐标系与地心惯性坐标系之间实际上是两个协议坐标系之间的转换。地心固定坐标系的协议坐标系是WGS84，地心惯性坐标系的协议坐标系是J2000.0。两个协议坐标系之间通过参数可以得到岁差、章动、恒星时旋转和极移矩阵，通过上述矩阵可以得到地心固定坐标系与地心惯性坐标系之间转换。
本体坐标系是以卫星质心作为原点，以卫星横轴作为X轴；以卫星飞行方向作为Y轴；以右手规则确定Z轴所构成的坐标系；当地轨道坐标系是以卫星质心为原点；以地球质心作为Z轴；X轴位于包含Z轴与卫星速度矢量V构成的平面上，垂直于Z轴，和速度矢量的夹角小于90度；以右手规则确定Y轴所构成的坐标系。
第一姿态旋转矩阵是本体坐标系到当地轨道坐标系所形成的矩阵。
具体的，在通过上述步骤S201获取姿态四元数后，根据本体坐标系到当地轨道坐标系之间的关系，按照坐标系转换规则，得到本体坐标系到当地轨道坐标系对应的第一姿态旋转矩阵。
S203、根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角。
由于本体坐标系到当地轨道坐标系的三轴旋转顺序不同，欧拉角的表达式也不同。可能的坐标轴旋转顺序构成的旋转矩阵形式多达12种。最常见的坐标轴旋转顺序是非对称轴旋转：3-2-1。因此，在本实施例中使用常见的非对称轴旋转：3-2-1的顺序作为预设的三轴旋转顺序。
拆分第一姿态旋转矩阵即为将第一姿态旋转矩阵中的姿态四元数进行拆分，得到姿态四元数对应的矢量部分、标量部分。
具体的，利用拆分后的第一姿态旋转矩阵以及非对称轴旋转：3-2-1的顺序对应的计算公式求算对应的欧拉角。
S204、根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角。
卫星的推扫成像方式可以利用姿态四元数的前缀编码进行确认。
拟合是将离散的数值使用某种函数来表达。在本实施例中利用函数关系将离散的欧拉角连接起来，使其能够形成相对光滑的曲线。
通过不同的推扫成像方式选择拟合方式，将欧拉角进行拟合得到模型化后的欧拉角。
S205、利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数。
具体的，为了保持物理模型的简洁，将模型化后的欧拉角序列形成的旋转矩阵利用坐标系转换规则转换回四元数。得到处理后的姿态四元数。姿态四元数每个分量的符号与初始姿态四元数保持一致。
本实施例获取星上计算机传输的姿态四元数以及星历数据。在获取成功后，利用坐标系转换规则，确定卫星的本体坐标系到当地轨道坐标系所对应的第一姿态旋转矩阵。进而根据预设的三轴旋转顺序，将第一姿态旋转矩阵中的参数进行拆分，得到对应的欧拉角，使得具有数学意义的姿态四元数转换为具有物理意义的欧拉角。进而通过卫星的推扫成像方式，将得到的欧拉角进行拟合，得到模型化后的欧拉角。通过这种方式将最初得到的姿态四元数中的噪声进行平滑处理，或者遵守姿态角、角速度和角加速度随时间的变化规律，提高数据内部相对几何精度。然后利用坐标系转换规则，将模型化后的欧拉角对应的序列所构成的旋转矩阵转换成处理后的姿态四元数。提高最初得到的姿态四元数的精度。
在一些实施例中，根据预设的三轴旋转顺序，将第一姿态旋转矩阵进行拆分并根据相应的计算公式进行计算得到对应的欧拉角。具体的，可以包括：确定卫星的本体坐标系到当地轨道坐标系的三轴旋转顺序；根据三轴旋转顺序，确定欧拉角的计算公式；基于计算公式、第一姿态旋转矩阵，确定对应的欧拉角。
本体坐标系(body)到当地轨道坐标系(LVLH)的三轴旋转顺序不同，欧拉角的表达式也不同。可以使用常用的非对称轴旋转顺序：3-2-1，相应的欧拉角的计算公式为：
计算卫星本体坐标系到轨道坐标系的旋转侧摆角，如公式1所示：
其中，roll代表卫星本体坐标系到轨道坐标系的旋转侧摆角；Mbody2lvlh(2,3)代表第一姿态旋转矩阵中坐标(2,3)的分量；Mbody2lvlh(3,3)代表第一姿态旋转矩阵中坐标(3,3)的分量。
计算卫星本体坐标系到轨道坐标系的俯仰角，如公式2所示：
pitch＝asin(Mbody2lvlh(1,3)) (2)
其中，pitch代表卫星本体坐标系到轨道坐标系的俯仰角；Mbody2lvlh(1,3)代表第一姿态旋转矩阵中坐标(1,3)的分量。
计算卫星本体坐标系到轨道坐标系的偏航角，如公式3所示：
其中，yaw代表卫星本体坐标系到轨道坐标系的偏航角；Mbody2lvlh(1,2)代表第一姿态旋转矩阵中坐标(1,2)的分量；Mbody2lvlh(1,1)代表第一姿态旋转矩阵中坐标(1,1)的分量。
上述三轴旋转顺序以及对应的计算公式均为举例，在本实施例中应用此顺序及对应的计算公式使计算过程更加简便。但本领域技术人员应当知晓，其他三轴旋转顺序以及对应的计算公式同样可以得到欧拉角。通过本实施例将最初得到的姿态四元数转换为具有物理意义的欧拉角，使其可以进行后续的操作，使其可以提高数据的精度。
在一些实施例中，可以通过卫星的推扫成像方式，选择合适的多项式进行拟合，得到模型化后的欧拉角。具体的，可以包括：根据卫星拍摄模式确定卫星的推扫成像方式；根据所述推扫成像方式确定拟合所需多项式；根据确定的多项式将所述欧拉角进行拟合得到模型化后的欧拉角。
卫星的推扫成像方式包括被动推扫成像和主动推扫成像。
被动推扫成像是指卫星通过星上的动量轮等装置从前一个姿态达到目标姿态并恢复平稳状态后开始拍摄。推扫期间除了需要在偏航方向上实时调整偏流角以补偿像移，由于在滚动和俯仰方向的角速度变化非常小，星上的姿态测量设备——星敏感器和陀螺都在小角速度情况下正常开展工作，所以星敏感器能够有效捕获恒星，以固定频率输出绝对定姿结果。陀螺以较高频率测量三轴角速度并输出。星上计算机系统对星敏感器输出的四元数和陀螺输出的三轴角速度进行扩展卡尔曼滤波。
扩展卡尔曼滤波适用于非线性程度不高的动力学系统，尤其适合卫星的稳态成像过程。因此我们可以认为被动推扫成像的时候卫星姿态确定的绝对精度与相对精度都很高，短轨拍摄周期内欧拉角随时间的变化曲线应该是线性的，长轨的欧拉角时间变化规律应符合二次多项式的曲线特征。
陀螺的输出包括真实输出、陀螺漂移和测量噪声，其中陀螺漂移包括常值漂移和与时间相关的部分，与时间相关部分的陀螺漂移可以认为是随机游走，是典型的有色噪声。假定星敏感器感器的测量噪声和陀螺的测量噪声都是白噪声且互相独立，陀螺的常值漂移噪声的一阶导数是常值，陀螺输出由真实角速度、常值漂移和测量噪声组成。状态参量为误差四元数的矢量部分和陀螺的常值漂移噪声。
主动推扫成像指的是由于卫星的姿态角速度较大，星敏感器难以跟踪到恒星并输出有效观测数据，这种情况下姿态确定系统一般不使用星敏感器参与定姿，而是对陀螺进行积分得到姿态角的估计值。在大角速度情况下，陀螺的运动学模型与稳态时有很大的不同，在稳态时可以忽略的误差在姿态机动过程中变得突出，使用陀螺测量的角速度进行积分得到姿态角的方法导致卫星主动推扫成像时的姿态确定的绝对精度较低，但由于陀螺的输出频率较高(通常大于200hz)，动态响应较快，因此我们可以认为主动推扫成像过程定姿系统输出的姿态数据相对精度较高。卫星主动推扫成像时，卫星的姿态角、角速度、角加速度均为随时间连续变化的曲线，姿态确定误差的来源是随时间不断增大的陀螺漂移。
具体的，在根据上述特征确定卫星的推扫成像方式后，对稳态成像的欧拉角随时间变化的曲线使用二次多项式拟合，对非稳态成像的欧拉角随时间变化的曲线使用五次多项式拟合。即对被动推扫成像的欧拉角随时间变化的曲线使用二次多项式拟合，对主动推扫成像的欧拉角随时间变化的曲线使用五次多项式拟合。
使用低阶多项式(二次多项式)来拟合被动推扫成像模式的姿态四元数，能够在最初获取的姿态四元数有较多噪声的情况对其进行平滑处理；使用高次多项式(五次多项式)拟合主动推扫成像模式的姿态四元数，能够使得拟合曲线与姿态四元数观测值尽量保持一致，严格遵守姿态角、角速度和角加速度随时间的变化规律，提高主动推扫成像的姿态四元数内部的相对几何精度。
被动推扫成像模式使用二次多项式，如公式4所示：
其中，roll代表卫星本体坐标系到轨道坐标系的旋转侧摆角；pitch代表卫星本体坐标系到轨道坐标系的俯仰角；yaw代表卫星本体坐标系到轨道坐标系的偏航角；t代表采样历元；a1、a2、a3、b1、b2、b3、c1、c2、c3为多项式系数。
主动推扫成像模式使用五次多项式，如公式5所示：
其中，roll代表卫星本体坐标系到轨道坐标系的旋转侧摆角；pitch代表卫星本体坐标系到轨道坐标系的俯仰角；yaw代表卫星本体坐标系到轨道坐标系的偏航角；t代表采样历元；a1、a2、a3、a4、a5、a6、b1、b2、b3、b4、b5、b6、c1、c2、c3、c4、c5、c6为多项式系数。
模型化后的欧拉角时间序列为连续的光滑曲线，处处可导，既能反映卫星的真实姿态状态也有利于地面预处理，提高图像内部几何精度。在得到拟合多项式系数an、bn、cn之后可以根据实际需求对姿态角数据进行内插、外延、改变初始历元和采样间隔等操作。例如，从星务时间240385262.2578开始，获取频率为200hz的被动推扫成像模式的欧拉角的数值序列就可以使用an、bn、cn计算得到，如公式6所示：
t＝240385262.2578×(1:0.005:200)
其中，roll代表卫星本体坐标系到轨道坐标系的旋转侧摆角；pitch代表卫星本体坐标系到轨道坐标系的俯仰角；yaw代表卫星本体坐标系到轨道坐标系的偏航角；a1、a2、a3、b1、b2、b3、c1、c2、c3为多项式系数。
在本实施例中通过确定卫星的推扫成像方式，确定使用哪种多项式对欧拉角进行拟合。将卫星的拍摄模式作为区分姿态四元数拟合的标准，根据不同成像模式卫星采用的不同的姿控和姿态确定的特点，对姿态四元数的观测值进行分析，结合实际姿态状况确定不同的拟合模型，精准描述和拟合卫星的实际姿态状况，提高图像内部相对几何精度，达到更好地效果。
在一些实施例中，可以通过姿态四元数的前缀编码，确定推扫成像方式。具体的，可以包括：提取姿态四元数的前缀编码；将前缀编码与预设前缀编码进行匹配；根据匹配结果，确定推扫成像方式。
预设前缀编码是根据星上计算机获取数据过程中所形成的编码进行设置的。
具体的，由于卫星在成像过程中带有编码，所以将姿态四元数中的编码和预设前缀编码进行对比，确认编码对应的推扫成像方式。
本实施例通过提取的姿态四元数的前缀编码，并将前缀编码与预设前缀编码进行比对，确定所提取的姿态四元数是如何获取的。并根据匹配的结果，确定卫星的推扫成像方式。这种方法可以清晰的确定出姿态四元数的对应的卫星的推扫成像模式，在一定程度上减少了工作量，提高了工作效率。
在一些实施例中，通过获取的姿态四元数、星历数据，并利用坐标系转换规则，确定第一姿态旋转矩阵。具体的，可以包括：确定姿态四元数的标量部分、矢量部分；利用姿态四元数的标量部分、矢量部分，确定本体坐标系到地心惯性坐标系对应的第二姿态旋转矩阵；利用旋转参数，确定地心固定坐标系与地心惯性坐标系之间的第一旋转矩阵；通过第一旋转矩阵、星历数据，得到卫星在地心惯性坐标系的瞬时位置与速度矢量；根据瞬时位置与速度矢量，确定当地轨道坐标系到地心惯性坐标系对应的第二旋转矩阵；根据第二姿态旋转矩阵与所述第二旋转矩阵，确定所述第一姿态旋转矩阵。
第二姿态旋转矩阵是本体坐标系到地心惯性坐标系的旋转矩阵。
第一旋转矩阵是地心固定坐标系与地心惯性坐标系之间的旋转矩阵。地心固定坐标系与地心惯性坐标系之间的第一旋转矩阵事实上这是两个协议坐标系之间的转换，因为瞬时坐标系是时变的。ECEF的协议坐标系是WGS84，ECI的协议坐标系是J2000.0。
第二旋转矩阵是当地轨道坐标系到地心惯性坐标系的旋转矩阵。
旋转参数是通过国际地球自转服务IERS(International Earth Rotation andReference Systems Service)每个月发布的地球定向参数EOP(Earth OrientationParameter)提供的。
星历数据是通过全球导航卫星系统得到的卫星在地心固定坐标系下的位置和速度矢量，星历数据的获取频率通常为1hz。在实际实现方式中，该频率可以进行相应的调整。
姿态四元数中的数据可以分为两部分：矢量部分与标量部分。如公式7所示：
q＝[q1 q2 q3 q4] (7)
其中，q1 q2 q3 q4是四元数分量，q4是四元数的标量部分，q1 q2 q3是四元数的矢量部分。
由于接收到的姿态四元数本身是卫星本体坐标系到地心惯性坐标系(ECI)的协议坐标系J2000.0的四元数，所以可以通过姿态四元数得到本体坐标系到地心惯性坐标系的第二姿态旋转矩阵，如公式8所示：
其中，Mbody2ECI代表本体坐标系到地心惯性坐标系的第二姿态旋转矩阵，q1、q2、q3、q4是四元数分量。
获取星历数据后，首先利用拉格朗日多项式内插法将星历数据根据姿态数据的起始时间和采样间隔进行内插处理，保证星历数据与姿态数据的起始时间和采样间隔一致。
第一旋转矩阵如公式9所示：
MECEF2ECI＝ABCD (9)
其中，ABCD代表两个协议坐标系之间的姿态旋转矩阵，分别是岁差矩阵、章动矩阵、恒星时旋转矩阵和极移矩阵。
通过第一旋转矩阵、星历数据，得到卫星在地心惯性坐标系的瞬时位置与速度矢量后，利用下列公式得到当地轨道坐标系到地心惯性坐标系的第二旋转矩阵。如公式10所示：
其中，
其中，MLVLH2ECI代表当地轨道坐标系到地心惯性坐标系的第二旋转矩阵；代表卫星在地心惯性坐标系下的瞬时位置；代表卫星在地心惯性坐标系下的瞬时位置对应的速度矢量；XS YS ZS代表卫星在地心惯性坐标系不同坐标轴下的瞬时位置；XVs YVs ZVs代表卫星在地心惯性坐标系不同坐标轴下的速度矢量。
通过上述方式得到当地轨道坐标系到地心惯性坐标系的第二旋转矩阵后，即可利用第二旋转矩阵、第二姿态旋转矩阵，通过相应的计算得到卫星的本体坐标系到当地轨道坐标系的第一姿态旋转矩阵，如公式11所示：
Mbody2LVLH＝MLVLH2ECI′×Mbody2ECI (11)
其中，Mbody2LVLH代表卫星的本体坐标系到当地轨道坐标系的第一姿态旋转矩阵；MLVLH2ECI代表当地轨道坐标系到地心惯性坐标系的第二旋转矩阵；Mbody2ECI代表卫星的本体坐标系到地心惯性坐标系的第二姿态旋转矩阵。
在一些实施例中，将所述星历数据根据姿态四元数的起始时间与采样间隔进行内插得到内插后的星历数据；根据所述内插后的星历数据、所述第一旋转矩阵，得到卫星在所述地心惯性坐标系的瞬时位置与速度矢量。
具体的，通过上述方式得到第一旋转矩阵以及进行内插后的星历数据后，利用第一旋转矩阵，将星历数据进行处理得到卫星在地心惯性坐标系的瞬时位置与速度矢量，如公式12、公式13所示：
其中，代表卫星在地心惯性坐标系下的瞬时位置；代表全球导航卫星系统获取的卫星在地心固定坐标系下的tj时刻的位置；t代表成像时间。
其中，代表卫星在地心惯性坐标系下的瞬时位置对应的速度矢量；代表全球导航卫星系统获取的卫星在地心固定坐标系下的tj时刻的位置对应的速度；t代表成像时间。
将获取的星历数据按照姿态四元数的起始时间与采样间隔进行内插，使星历数据与姿态四元数的起始时间与采样间隔保持一致，得到内插后的星历数据。进而使用第一旋转矩阵，利用相应的计算公式得到卫星在地心惯性坐标系的瞬时位置与速度矢量。通过这种方式使星历数据与姿态四元数的起始时间和采样间隔保持一致，更有利于后续的计算以及其余操作。
在一些实施例中，利用坐标系转换规则，将模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的第二姿态旋转矩阵；根据预设分量计算规则、处理后的第二姿态旋转矩阵中的参数，计算四个分量的值；将四个分量的值进行排序，确定最大值对应的分量；根据预设四元数计算规则，利用所述最大值对应的分量计算四元数。
通过上述方式得到模型化后的欧拉角以后，确定所有模型化后的欧拉角的序列对应的旋转矩阵，并利用坐标系转换规则转换成处理后的第二姿态旋转矩阵，如公式14所示：
Mbody2_ECI＝MLVLH2ECI×Mbody2LVLH (14)
其中，Mbody2_ECI代表处理后的第二姿态旋转矩阵；MLVLH2ECI代表第二旋转矩阵；Mbody2LVLH代表第一姿态旋转矩阵。
使用姿态四元数的分量来表达姿态旋转矩阵：
其中，Mbody2_ECI代表处理后的第二姿态旋转矩阵；q1、q2、q3代表姿态四元数的矢量部分；q4代表姿态四元数的标量部分。
得到处理后的第二姿态旋转矩阵后，利用以下步骤得到本体坐标系到地心惯性坐标系的姿态四元数。如公式15所示，首先计算四个分量：
其中，d1、d2、d3、d4代表计算姿态四元数的中间变量； 代表处理后的第二姿态旋转矩阵对角线矩阵。
通过上述方式计算出中间变量后，确认哪一中间变量值最大，并根据确定的最大中间变量值计算姿态四元数的分量。若d1最大，则利用d1计算姿态四元数的分量，如公式16所示：
其中，q1、q2、q3、q4代表姿态四元数的分量； 代表处理后的第二姿态旋转矩阵对角线矩阵。
若d2最大，则利用d2计算姿态四元数的分量，如公式17所示：
其中，q1、q2、q3、q4代表姿态四元数的分量； 代表处理后的第二姿态旋转矩阵对角线矩阵。
若d3最大，则利用d3计算姿态四元数的分量，如公式18所示：
其中，q1、q2、q3、q4代表姿态四元数的分量； 代表处理后的第二姿态旋转矩阵对角线矩阵。
若d4最大，则利用d4计算姿态四元数的分量，如公式19所示：
其中，q1、q2、q3、q4为姿态四元数的分量； 为处理后的第二姿态旋转矩阵对角线矩阵。
通过上述方式确定姿态四元数的分量部分后，计算单位四元数。如公式20所示：
其中，q代表单位四元数。
本实施例利用相应的计算公式计算对应的分量的值，并进行排序，得到最大值对应的分量。进而按照预设四元数计算规则，利用最大值对应的分量计算四元数的标量部分以及矢量部分。通过得到最大值对应的分量可以避免了计算过程中分母过小的问题，可以得到有效的四元数分量，减小了误差。
本申请所设计的处理流程逻辑清晰，参量之间的转换严密，充分考虑目前敏捷型卫星在主动成像和被动成像时姿态方面的真实表现，最大限度的提升图像内部的相对几何精度。
图3为本申请的流程图。结合图3来描述本申请的内容。图中英文字母roll代表卫星本体坐标系到轨道坐标系的旋转侧摆角、pitch代表俯仰角，yaw代表偏航角。
本申请为顾及卫星实际姿态状况的四元数拟合方法，包括以下步骤：
步骤一、卫星下传的姿态数据为卫星本体坐标系到地心惯性坐标系的协议坐标系J2000.0的四元数，通常为星上计算机对星敏四元数和陀螺角速度数据进行联合滤波的结果，使用坐标系变换，获得本体坐标系到LVLH——local vertical local horizontal坐标系即瞬时轨道坐标系的欧拉角，欧拉角具有明确的物理意义，即本体坐标系三轴相对于LVLH的旋转角。下面详细介绍坐标系变换的过程和方法。
星上姿态确定系统提供的四元数为q，姿态四元数中的数据可以分为两部分：矢量部分与标量部分。如上述公式7所示。
1.本体坐标系到地心惯性坐标系的旋转矩阵，如上述公式8所示。
2.地心固定坐标系与地心惯性坐标系之间的旋转矩阵，如上述公式9所示。
事实上这是两个协议坐标系之间的转换，因为瞬时坐标系是时变的。ECEF的协议坐标系是WGS84，ECI的协议坐标系是J2000.0。
ABCD代表的是两个协议坐标系之间的旋转矩阵，分别是——岁差、章动、恒星时旋转和极移矩阵。每个旋转矩阵通过复杂的旋转参数完成，这些通过IERS每个月发布的EOP提供。
由GNSS(Global Navigation Satellite System)得到卫星在地心固定坐标系下的位置和速度矢量，星历数据的获取频率通常为1hz，在进行坐标系转换过程中，需要将星历数据根据姿态数据的起始时间和采样间隔进行内插处理，使用拉格朗日多项式内插法对星历数据进行内插，保证星历数据与姿态数据的起始时间和采样间隔一致。
拉格朗日(Lagrange)多项式内插法被广泛采用，因为这种内插法速度快且易于编程。拉格朗日插值法的原理十分简单：已知函数y＝f(x)的n+1个节点x1,x2,x3…xn及其对应的函数值y1,y2,y3…yn对插值区间内任一点x，可用下面的插值多项式来计算函数值，如公式21所示：
其中，自变量x和y之间也可以不存在函数关系，只是简单的对应关系，例如星历参数中，成像瞬时与位置矢量和速度矢量之间的对应关系。卫星成像瞬时的星历位置矢量和速度矢量通过时刻前后至少8个已知的星历参数进行插值计算得到就可以保证较高的内插精度。因此，如果我们要计算某个时刻卫星的位置和速度，可以通过下式进行差值得到。如上述公式12、公式13所示。
3.当地轨道坐标系(LVLH)到地心惯性坐标系(ECI)的旋转矩阵该旋转矩阵使用卫星在ECI坐标系下的每个瞬时的位置矢量和速度矢量表示。如上述公式10所示。
4.完成了上述坐标系转换之后，得到如上述公式11所示的旋转矩阵该矩阵有多种表达方式，例如四元数、欧拉角或方向余弦、Gibbs矢量等，反之，也可以通过姿态旋转矩阵来求算相关的姿态参数，在这个环节通过姿态矩阵直接求算欧拉角。需要注意的是本体坐标系(body)到当地轨道坐标系(LVLH)的三轴旋转顺序不同，欧拉角的表达式也不同，可能的坐标轴旋转顺序构成的旋转矩阵形式多达12种。最常见的坐标轴旋转顺序是非对称轴旋转：3-2-1，对应的欧拉角的求算公式如上述公式1、公式2、公式3所示。
步骤二、根据卫星拍摄模式判断卫星是被动推扫成像还是主动推扫成像。
姿态确定的绝对精度主要影响外方位元素的精度，进而影响卫星影像无控定位的精度；姿态确定的相对精度主要影响影像内部沿线阵推扫方向的相对几何精度。相对于影像的绝对几何精度来说，内部的相对精度更为重要。
步骤三、姿态角的初始形式为离散的采样时间和对应的欧拉角，使用最小二乘法求算以时间为自变量的多项式系数，被动推扫成像方式使用二次多项式，如上述公式4所示。
主动推扫成像方式使用五次多项式，如上述公式5所示。
t为采样历元，an、bn、cn为多项式系数。模型化后的欧拉角时间序列为连续的光滑曲线，处处可导，既能反映卫星的真实姿态状态也有利于地面预处理，提高图像内部几何精度。
步骤四、为了保持严密物理模型的简洁，最后将模型化后的欧拉角序列经过坐标系旋转转换回姿态四元数。姿态四元数每个分量的符号与初始姿态四元数保持一致，且为单位四元数。欧拉角构成的旋转矩阵经过坐标系转换可以得到从本体坐标系到地心惯性坐标系的转换矩阵，该矩阵即姿态四元数的姿态旋转矩阵，从该矩阵中求算姿态四元数每个分量，如上述公式14所示。
通过下面的步骤得到本体坐标系到地心惯性系的姿态四元数：
1.计算四个分量，如上述公式15所示。
2.寻找d1-d4中最大的值，计算四元数。
如果d1最大，则利用d1计算姿态四元数的分量，如上述公式16所示。
如果d2最大，则利用d2计算姿态四元数的分量，如上述公式17所示。
若d3最大，则利用d3计算姿态四元数的分量，如上述公式18所示。
若d4最大，则利用d4计算姿态四元数的分量，如上述公式19所示。
3.计算单位四元数，如上述公式20所示。
本申请关键创新点：将敏捷卫星的拍摄模式作为区分姿态数据拟合模型的标准，根据不同成像模式卫星采用的不同的姿控和姿态确定的特点，对姿态数据的观测值进行分析，结合实际姿态状况确定不同的拟合模型，精准描述和拟合卫星的实际姿态状况，提高图像内部相对几何精度。
本申请解决了目前地面预处理系统对星上下传的姿态观测值进行内插时只考虑其数学意义，不考虑其实际物理意义的问题。同时兼顾了姿态数据的数学意义和物理意义，能够精准拟合卫星实际姿态状况同时对姿态数据中的噪声进行平滑处理，提高了图像内部的相对几何精度。
分别对某敏捷卫星的被动成像和主动成像方式的姿态数据按照本发明的步骤进行处理并对图像内部几何精度进行检测。
1.被动推扫成像方式
得到某敏捷卫星被动推扫成像方式的轨道号为1721的前1秒的姿态数据，原始姿态数据是200hz的离散四元数，存在丢帧现象。使用本发明步骤一描述的方法将四元数转换为欧拉角，由于是稳态推扫的被动推扫成像方式，使用二次多项式拟合欧拉角。图4a-4c是被动推扫方式下三轴欧拉角原始姿态数据与拟合欧拉角曲线。
在对相机内方位参数进行在轨几何检校后，使用本发明的拟合方案处理的姿态数据建立严密物理模型，对L1图像内部几何精度进行检测。参考高精度DOM和DEM数据，对L1标准景数据进行稀疏控制点采集，使用仿射变换模型在像方优化rpc系数，对校正后图像的精度进行检查。表1为某卫星部分被动成像数据内部几何精度检测结果。
表1被动成像数据内部几何精度
从图4a-4c可以发现某敏捷成像卫星的稳态成像下传姿态信号有系统噪声，将其进行傅里叶变换会发现每轨信号抖动的频率和振幅都不一样，且在图像上没有发现任何由于姿态抖动造成的异常现象——低频噪声导致地物扭曲/高频噪声导致纹理模糊，因此判断这是姿态确定系统的某些失误导致的数据误差，并不是卫星平台姿态状态的真实反映。不能直接使用这样的姿态对图像进行几何处理，因此使用二次多项式对欧拉角进行拟合既能反映真实的平台状况又能对姿态数据进行噪声平滑。通过表1的内部几何精度检查可以证明姿态数据的拟合方式是合理有效的。
2.主动推扫成像方式
得到某敏捷卫星主动推扫成像方式的轨道号为7CD的原始姿态数据，将姿态四元数转换为欧拉角，为了覆盖图像拍摄时间需要将欧拉角沿当前历元向前外延，分别使用二次多项式和五次多项式拟合原始姿态数据进行外延。图5a-5c是首景数据需要外延的姿态，分别为不同角度状态下原始姿态离散点、二次多项式拟合曲线和五次多项式拟合曲线。
对首景数据采集检查点参考高精度DOM影像来验证使用不同拟合方法的姿态数据的物方定位情况。图6a-6b是使用二次多项式拟合的主动推扫与被动推扫姿态数据的物方残差；图7a-7b是使用五次多项式拟合的主动推扫与被动推扫姿态数据物方剩余残差。
图5a-5c表现了使用两种拟合得到的姿态角曲线的外延部分与原始姿态角数据的关系。通过图6a-6b、图7a-7b可以看出五次多项式拟合姿态的误差表现在行列方向均为线性函数，且精度较高，使用很少的控制点就可以实现对rpc系数的像方优化，而二次多项式拟合姿态的物方误差在行列方向的表现均为2次函数，在一景范围内误差变化剧烈，难以使用控制点进行rpc的像方优化。对于主动推扫得到的姿态数据来说，最重要的处理原则是尽量拟合原始姿态曲线的趋势，信任使用陀螺输出的角速度得到的积分数据。
对于敏捷成像卫星，其在被动成像和主动成像过程中所使用的姿态确定和姿态控制完全不同，将卫星成像模式作为区分的标志来对成像姿态角采用不同的拟合模型，通过上述实施例证明了本申请的严密和有效。
图8为本申请一实施例提供的一种面向敏捷卫星的姿态四元数处理装置的结构示意图，如图8所示的，本实施例的面向敏捷卫星的姿态四元数处理装置800包括：数据获取模块801、第一姿态旋转矩阵确定模块802、欧拉角确定模块803、欧拉角拟合模块804、处理模块805。
数据获取模块801，用于获取姿态四元数、星历数据；
第一姿态旋转矩阵确定模块802，用于根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系对应的第一姿态旋转矩阵；
欧拉角确定模块803，用于根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角；
欧拉角拟合模块804，用于根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角；
处理模块805，用于利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数。
可选的，所述欧拉角确定模块803在根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角时，具体用于：
确定卫星的本体坐标系到当地轨道坐标系的三轴旋转顺序；
根据所述三轴旋转顺序，确定欧拉角的计算公式；
基于所述计算公式、所述第一姿态旋转矩阵，确定对应的欧拉角。
可选的，所述欧拉角拟合模块804在根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角时，具体用于：
根据卫星拍摄模式确定卫星的推扫成像方式；
根据所述推扫成像方式确定拟合所需多项式；
根据确定的多项式将所述欧拉角进行拟合得到模型化后的欧拉角。
可选的，所述欧拉角拟合模块804在根据卫星拍摄模式确定卫星的推扫成像方式时，具体用于：
提取所述姿态四元数的前缀编码；
将所述前缀编码与预设前缀编码进行匹配；
根据匹配结果，确定所述推扫成像方式。
可选的，所述第一姿态旋转矩阵确定模块802在根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系的第一姿态旋转矩阵时，具体用于：
确定所述姿态四元数的标量部分、矢量部分；
利用所述姿态四元数的标量部分、矢量部分，确定本体坐标系到地心惯性坐标系对应的第二姿态旋转矩阵；
利用旋转参数，确定地心固定坐标系与地心惯性坐标系之间的第一旋转矩阵；
通过所述第一旋转矩阵、所述星历数据，得到卫星在所述地心惯性坐标系的瞬时位置与速度矢量；
根据所述瞬时位置与速度矢量，确定当地轨道坐标系到地心惯性坐标系对应的第二旋转矩阵；根据所述第二姿态旋转矩阵与所述第二旋转矩阵，确定所述第一姿态旋转矩阵。
可选的，所述第一姿态旋转矩阵确定模块802在通过所述第一旋转矩阵、所述星历数据，得到卫星在所述地心惯性坐标系的瞬时位置与速度矢量时具体用于：
将所述星历数据根据所述姿态四元数的起始时间与采样间隔进行内插得到内插后的星历数据；根据所述内插后的星历数据、所述第一旋转矩阵，得到卫星在所述地心惯性坐标系的瞬时位置与速度矢量。
可选的，所述处理模块805在利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数时，具体用于：
利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的第二姿态旋转矩阵；
根据预设分量计算规则、所述处理后的第二姿态旋转矩阵中的参数，计算四个分量的值；将所述四个分量的值进行排序，确定最大值对应的分量；
根据预设四元数计算规则，利用所述最大值对应的分量计算四元数。
本实施例的装置，可以用于执行上述任一实施例的方法，其实现原理和技术效果类似，此处不再赘述。
图9为本申请一实施例提供的一种面向敏捷卫星的姿态四元数处理设备的结构示意图，如图9所示，本实施例的面向敏捷卫星的姿态四元数处理设备900可以包括：存储器901和处理器902。
存储器901上存储有能够被处理器902加载并执行上述实施例中方法的计算机程序。
其中，处理器902和存储器901相连，如通过总线相连。
可选地，面向敏捷卫星的姿态四元数处理设备900还可以包括收发器。需要说明的是，实际应用中收发器不限于一个，该面向敏捷卫星的姿态四元数处理设备900的结构并不构成对本申请实施例的限定。
处理器902可以是CPU(Central Processing Unit，中央处理器)，通用处理器，DSP(Digital Signal Processor，数据信号处理器)，ASIC(Application SpecificIntegrated Circuit，专用集成电路)，FPGA(Field Programmable Gate Array，现场可编程门阵列)或者其他可编程逻辑器件、晶体管逻辑器件、硬件部件或者其任意组合。其可以实现或执行结合本申请公开内容所描述的各种示例性的逻辑方框，模块和电路。处理器902也可以是实现计算功能的组合，例如包含一个或多个微处理器组合，DSP和微处理器的组合等。
总线可包括一通路，在上述组件之间传送信息。总线可以是PCI(PeripheralComponent Interconnect，外设部件互连标准)总线或EISA(Extended Industry StandardArchitecture，扩展工业标准结构)总线等。总线可以分为地址总线、数据总线、控制总线等。为便于表示，图中仅用一条粗线表示，但并不表示仅有一根总线或一种类型的总线。
存储器901可以是ROM(Read Only Memory，只读存储器)或可存储静态信息和指令的其他类型的静态存储设备，RAM(Random Access Memory，随机存取存储器)或者可存储信息和指令的其他类型的动态存储设备，也可以是EEPROM(Electrically ErasableProgrammable Read Only Memory，电可擦可编程只读存储器)、CD-ROM(Compact DiscRead Only Memory，只读光盘)或其他光盘存储、光碟存储(包括压缩光碟、激光碟、光碟、数字通用光碟、蓝光光碟等)、磁盘存储介质或者其他磁存储设备、或者能够用于携带或存储具有指令或数据结构形式的期望的程序代码并能够由计算机存取的任何其他介质，但不限于此。
存储器901用于存储执行本申请方案的应用程序代码，并由处理器902来控制执行。处理器902用于执行存储器901中存储的应用程序代码，以实现前述方法实施例所示的内容。
其中，电子设备包括但不限于：移动电话、笔记本电脑、数字广播接收器、PDA(个人数字助理)、PAD(平板电脑)、PMP(便携式多媒体播放器)、车载终端(例如车载导航终端)等等的移动终端以及诸如数字TV、台式计算机等等的固定终端。还可以为服务器等。图9示出的电子设备仅仅是一个示例，不应对本申请实施例的功能和使用范围带来任何限制。
本实施例的电子设备，可以用于执行上述任一实施例的方法，其实现原理和技术效果类似，此处不再赘述。
本申请还提供一种计算机可读存储介质，存储有能够被处理器加载并执行如上实施例中的方法的计算机程序。
本领域普通技术人员可以理解：实现上述各方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成。前述的程序可以存储于一计算机可读取存储介质中。该程序在执行时，执行包括上述各方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
Claims (8)
1.一种面向敏捷卫星的姿态四元数处理方法，其特征在于，包括：
获取姿态四元数、星历数据；
根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系对应的第一姿态旋转矩阵；
根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角；
根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角；
利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数；
所述根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角包括：
根据卫星拍摄模式确定卫星的推扫成像方式；
根据所述推扫成像方式确定拟合所需多项式；
根据确定的多项式将所述欧拉角进行拟合得到模型化后的欧拉角；
所述利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数包括：
利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的第二姿态旋转矩阵；
根据预设分量计算规则、所述处理后的第二姿态旋转矩阵中的参数，计算四个分量的值；
将所述四个分量的值进行排序，确定最大值对应的分量；
根据预设四元数计算规则，利用所述最大值对应的分量计算四元数；
所述预设分量计算规则为在得到处理后的第二姿态旋转矩阵后，利用以下步骤计算四个分量的值：
；
其中，代表计算姿态四元数的中间变量；代表处理后的第二姿态旋转矩阵的对角线矩阵。
2.根据权利要求1所述的方法，其特征在于，所述根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角包括：
确定卫星的本体坐标系到当地轨道坐标系的三轴旋转顺序；
根据所述三轴旋转顺序，确定欧拉角的计算公式；
基于所述计算公式、所述第一姿态旋转矩阵，确定对应的欧拉角。
3.根据权利要求1所述的方法，其特征在于，所述根据卫星拍摄模式确定卫星的推扫成像方式包括：
提取所述姿态四元数的前缀编码；
将所述前缀编码与预设前缀编码进行匹配；
根据匹配结果，确定所述推扫成像方式。
4.根据权利要求1所述的方法，其特征在于，所述根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系的第一姿态旋转矩阵包括：
确定所述姿态四元数的标量部分、矢量部分；
利用所述姿态四元数的标量部分、矢量部分，确定本体坐标系到地心惯性坐标系对应的第二姿态旋转矩阵；
利用旋转参数，确定地心固定坐标系与地心惯性坐标系之间的第一旋转矩阵；
通过所述第一旋转矩阵、所述星历数据，得到卫星在所述地心惯性坐标系的瞬时位置与速度矢量；
根据所述瞬时位置与速度矢量，确定当地轨道坐标系到地心惯性坐标系对应的第二旋转矩阵；
根据所述第二姿态旋转矩阵与所述第二旋转矩阵，确定所述第一姿态旋转矩阵。
5.根据权利要求4所述的方法，其特征在于，所述通过所述第一旋转矩阵、所述星历数据，得到卫星在所述地心惯性坐标系的瞬时位置与速度矢量包括：
将所述星历数据根据所述姿态四元数的起始时间与采样间隔进行内插得到内插后的星历数据；
根据所述内插后的星历数据、所述第一旋转矩阵，得到卫星在所述地心惯性坐标系的瞬时位置与速度矢量。
6.一种面向敏捷卫星的姿态四元数处理装置，其特征在于，包括：
数据获取模块，用于获取姿态四元数、星历数据；
第一姿态旋转矩阵确定模块，用于根据所述姿态四元数、所述星历数据，利用坐标系转换规则，确定本体坐标系到当地轨道坐标系对应的第一姿态旋转矩阵；
欧拉角确定模块，用于根据预设的三轴旋转顺序，拆分所述第一姿态旋转矩阵，得到对应的欧拉角；
欧拉角拟合模块，用于根据确定出的卫星的推扫成像方式，将所述欧拉角进行拟合，得到模型化后的欧拉角；
处理模块，用于利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的姿态四元数；
所述欧拉角拟合模块具体用于：
根据卫星拍摄模式确定卫星的推扫成像方式；
根据所述推扫成像方式确定拟合所需多项式；
根据确定的多项式将所述欧拉角进行拟合得到模型化后的欧拉角；
所述处理模块具体用于：
利用所述坐标系转换规则，将所述模型化后的欧拉角的序列对应的旋转矩阵转换成处理后的第二姿态旋转矩阵；
根据预设分量计算规则、所述处理后的第二姿态旋转矩阵中的参数，计算四个分量的值；
将所述四个分量的值进行排序，确定最大值对应的分量；
根据预设四元数计算规则，利用所述最大值对应的分量计算四元数；
所述预设分量计算规则为在得到处理后的第二姿态旋转矩阵后，利用以下步骤计算四个分量的值：
；
其中代表计算姿态四元数的中间变量；代表处理后的第二姿态旋转矩阵的对角线矩阵。
7.一种面向敏捷卫星的姿态四元数处理设备，其特征在于，包括：存储器和处理器；
所述存储器，用于存储程序指令；
所述处理器，用于调用并执行所述存储器中的程序指令，执行如权利要求1-5任一项所述的方法。
8.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质中存储有计算机程序；所述计算机程序被处理器执行时，实现如权利要求1-5任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211620133.XA CN115829879B (zh) | 2022-12-15 | 2022-12-15 | 面向敏捷卫星的姿态四元数处理方法、装置、设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211620133.XA CN115829879B (zh) | 2022-12-15 | 2022-12-15 | 面向敏捷卫星的姿态四元数处理方法、装置、设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115829879A CN115829879A (zh) | 2023-03-21 |
| CN115829879B true CN115829879B (zh) | 2023-10-27 |
Family
ID=85545977
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211620133.XA Active CN115829879B (zh) | 2022-12-15 | 2022-12-15 | 面向敏捷卫星的姿态四元数处理方法、装置、设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115829879B (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104848860A (zh) * | 2015-05-19 | 2015-08-19 | 北京空间飞行器总体设计部 | 一种敏捷卫星成像过程姿态机动规划方法 |
| CN105160125A (zh) * | 2015-09-24 | 2015-12-16 | 航天东方红卫星有限公司 | 一种星敏感器四元数的仿真分析方法 |
| CN105444778A (zh) * | 2015-11-10 | 2016-03-30 | 北京空间飞行器总体设计部 | 一种基于成像几何反演的星敏感器在轨定姿误差获取方法 |
| WO2016077018A1 (en) * | 2014-11-14 | 2016-05-19 | Intel Corporation | Four-dimensional morton coordinate conversion processors, methods, systems, and instructions |
| CN108344396A (zh) * | 2018-01-24 | 2018-07-31 | 浙江大学 | 一种敏捷卫星斜条带成像模式姿态计算方法 |
| CN110211054A (zh) * | 2019-04-28 | 2019-09-06 | 张过 | 一种星载推扫式光学传感器无畸变影像制作方法 |
| CN113190028A (zh) * | 2021-03-31 | 2021-07-30 | 北京控制工程研究所 | 一种敏捷卫星指向控制方法及系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7874519B2 (en) * | 2006-02-25 | 2011-01-25 | Space Systems/Loral, Inc. | Spacecraft three-axis attitude acquisition from sun direction measurement |
-
2022
- 2022-12-15 CN CN202211620133.XA patent/CN115829879B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2016077018A1 (en) * | 2014-11-14 | 2016-05-19 | Intel Corporation | Four-dimensional morton coordinate conversion processors, methods, systems, and instructions |
| CN104848860A (zh) * | 2015-05-19 | 2015-08-19 | 北京空间飞行器总体设计部 | 一种敏捷卫星成像过程姿态机动规划方法 |
| CN105160125A (zh) * | 2015-09-24 | 2015-12-16 | 航天东方红卫星有限公司 | 一种星敏感器四元数的仿真分析方法 |
| CN105444778A (zh) * | 2015-11-10 | 2016-03-30 | 北京空间飞行器总体设计部 | 一种基于成像几何反演的星敏感器在轨定姿误差获取方法 |
| CN108344396A (zh) * | 2018-01-24 | 2018-07-31 | 浙江大学 | 一种敏捷卫星斜条带成像模式姿态计算方法 |
| CN110211054A (zh) * | 2019-04-28 | 2019-09-06 | 张过 | 一种星载推扫式光学传感器无畸变影像制作方法 |
| CN113190028A (zh) * | 2021-03-31 | 2021-07-30 | 北京控制工程研究所 | 一种敏捷卫星指向控制方法及系统 |
Non-Patent Citations (3)
| Title |
|---|
| Quaternion based linear time-varying model predictive attitude control for satellites with two reaction wheels;Ali Golzari等;《Aerospace Science and Technology》;第98卷 * |
| 大气折射对光学卫星遥感影像几何定位的影响分析;严明等;《测绘学报》;第44卷(第09期);第995-1002页 * |
| 超敏捷卫星姿态机动与稳定控制算法研究;曲直;《中国博士学位论文全文数据库 工程科技Ⅱ辑》(第08期);第C031-8页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115829879A (zh) | 2023-03-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11295472B2 (en) | Positioning method, positioning apparatus, positioning system, storage medium, and method for constructing offline map database | |
| Peng et al. | Pose measurement and motion estimation of space non-cooperative targets based on laser radar and stereo-vision fusion | |
| CN108592950B (zh) | 一种单目相机和惯性测量单元相对安装角标定方法 | |
| WO2020140431A1 (zh) | 相机位姿确定方法、装置、电子设备及存储介质 | |
| CN111156998A (zh) | 一种基于rgb-d相机与imu信息融合的移动机器人定位方法 | |
| US20170254906A1 (en) | Acceleration of real time computer vision processing on uavs through gps attitude estimation | |
| CN111380514A (zh) | 机器人位姿估计方法、装置、终端及计算机存储介质 | |
| CN110954134B (zh) | 陀螺仪偏差校正方法、校正系统、电子设备及存储介质 | |
| JP2007183138A (ja) | 小型姿勢センサ | |
| CN113551665B (zh) | 一种用于运动载体的高动态运动状态感知系统及感知方法 | |
| CN110824453A (zh) | 一种基于图像跟踪与激光测距的无人机目标运动估计方法 | |
| Sun et al. | Line‐of‐sight rate estimation based on UKF for strapdown seeker | |
| CN110645976A (zh) | 一种移动机器人的姿态估计方法及终端设备 | |
| CN111998870A (zh) | 一种相机惯导系统的标定方法和装置 | |
| CN115829879B (zh) | 面向敏捷卫星的姿态四元数处理方法、装置、设备 | |
| CN110160530B (zh) | 一种基于四元数的航天器姿态滤波方法 | |
| Delabie et al. | Highly efficient attitude-estimation algorithm for star trackers using optimal image matching | |
| US20240134033A1 (en) | Method for determining a movement state of a rigid body | |
| CN115655305A (zh) | 外参标定方法、装置、计算设备、存储介质和车辆 | |
| Xiaoqian et al. | Nonlinear Extended Kalman Filter for Attitude Estimation of the Fixed‐Wing UAV | |
| CN112230194B (zh) | 一种基于平移阵列的解模糊方法、设备及存储介质 | |
| Do et al. | DeRO: Dead Reckoning Based on Radar Odometry With Accelerometers Aided for Robot Localization | |
| Deans et al. | A sun tracker for planetary analog rovers | |
| Feetham et al. | Single camera absolute motion based digital elevation mapping for a next generation planetary lander | |
| Delabie et al. | Robustness and efficiency improvements for star tracker attitude estimation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |