1.本发明属于航天测量与控制领域，涉及一种卫星受力的动量分析方法。
背景技术：2.在轨卫星运行期间面临着复杂的空间环境威胁，空间碎片、微流星体在空间高速运行并具有极高的动能，如果与航天器发生碰撞，会给后者造成严重危害。我国具备对较大空间目标监视及卫星规避控制能力，但受技术手段限制，微小空间碎片、微流星体仍对我国在轨运行航天器产生很大威胁。在受到无法监控的微小碎片冲击的小概率碰撞情况下，卫星受瞬时冲力影响引发后果，一是导致星上器部件受损或者整星致命性损伤，甚至整星解体；二是卫星姿态可能失稳，造成异常喷气控制，并导致卫星轨道变化；三是瞬时受力冲击会对整星带来速度增量，从而也会引起轨道变化。
3.针对此类异常，如何分析卫星受到的瞬时受力冲击是由内因导致还是外物撞击所致，对于异常的准确定位、规范处置、国家对此类空间事件的正确研判和对外发声非常重要。目前通常通过卫星遥测数据变化和事后定轨方式对卫星受力情况进行分析，通过异常发生前后的遥测数据分析法难以准确地追因溯源，无法直接证明卫星异常为外力撞击还是内力导致；通过事后定轨方法仅能评估已知力模型的影响，不能评估未知力的大小、受力时刻和受力方向问题。
技术实现要素：4.为了克服现有技术的不足，本发明提供一种识别卫星空间瞬时受力冲击的动量分析方法，基于卫星的姿轨控设计特点，通过分析确定喷气脉宽、喷气时长、喷气推力大小、喷气方向、轨道系三轴姿态等，计算异常发生后多个姿控发动机喷气产生的速度增量导致卫星半长轴的变化，与精密定轨结果进行比较，识别导致姿控异常喷气是否由于外力冲击产生，能够为卫星姿态异常定位、影响域分析、溯源分析及对撞击物质量与大小评估提供直接数据支撑。
5.本发明解决其技术问题所采用的技术方案包括以下步骤：
6.(1)确定当前条件下姿控推力大小曲线，确定不同喷气时长条件下姿控推力器推力大小；
7.(2)令卫星推力器喷气脉宽序列为{tk}，相应脉宽的推力序列为{f
tk
}，各推力器的安装矩阵为m
bt
，计算各推力器在各控制节拍中的冲量ik＝f
tk
tk；
8.如果推力器非正向安装，根据各推力器在卫星本体坐标系下的安装方向，计算卫星本体坐标系下各推力器在各控制节拍中冲量矢量在卫星本体坐标系下的矢量表示i
bk
＝m
bt
ik；
9.(3)计算卫星本体系到卫星轨道坐标系的转换矩阵m
ob
＝m
oi
·
(m
bi
)
t
，其中，m
bi
表示惯性系到本体系的转换矩阵，m
oi
表示2000.0惯性坐标系至轨道坐标系的转换矩阵；
10.(4)根据第j个推力器产生的冲量以及卫星在运动方向的冲量，计算卫星半长轴变
化量δag；
11.(5)利用异常发生前后的精密轨道确定数据，计算异常前后实际的半长轴变化量δa，得到差值δas＝δa-δag，该差值如果大于设定阈值，即判定是受到空间外力冲击产生。
12.所述的步骤(3)包括以下步骤：以卫星异常前一拍姿态测量数据为初值，令j2000惯性系下的姿态用四元数为q
gu0
，将陀螺角速度向量ω
bk
扩展为4维向量后，积分计算惯性系姿态四元数q
guk
＝[q0q1q2q3]；根据惯性姿态四元数，计算得到惯性系到本体系转换矩阵m
bi
；由卫星轨道坐标系定义可知空间矢量由2000.0惯性坐标系至轨道坐标系的转换矩阵m
oi
；计算得到卫星本体系到卫星轨道坐标系的转换矩阵m
ob
＝m
oi
·
(m
bi
)
t
。
[0013]
所述的惯性系到本体系转换矩阵
[0014][0015]
所述的2000.0惯性坐标系至轨道坐标系的转换矩阵m
oi
为：
[0016][0017][0018]moi
(1,m)＝m
oi
(2,p)
×moi
(3,k)
[0019]
其中，i＝1,2,3对应于转换矩阵m
oi
中每个行向量的三个分量，m
oi
(2,p)和m
oi
(3,k)为矩阵m
oi
中第二行和第三行的行向量，表示为空间矢量r的单位矢量，空间矢量r的速度矢量。
[0020]
所述的步骤(4)中，第j个推力器在轨道系中产生的冲量n表示第j个推力器在喷气时间内共产生n个喷气脉宽，i
bji
表示本体系下第j个推力器的第i个脉宽产生的冲量。
[0021]
所述的步骤(4)中，假定某卫星有m个姿控推力器，则各方向推力为m个姿控推力器在各方向产生的分量之和，在轨道系中前进方向t方向的冲量
[0022]
所述的步骤(4)中，令卫星质量为m，卫星轨道运行速率为n，轨道偏心率e，真近点角为f，则卫星姿控喷气导致的半长轴变化量i
ox
表示卫星在x方向的冲量。
[0023]
所述的步骤(5)中，差值δas＝δa-δag如果大于计算误差的3σ，即判定是受到空间外力冲击产生。
[0024]
本发明的有益效果是：完整描述了利用卫星姿控推力器异常喷气冲量分析估算卫星轨道半长轴变化量，通过与精密轨道确定数据差值的方法，实现卫星受外力冲击分析研判。本发明介绍了流程和实现的方法，后续还可根据情况变化，通过增加或减少计算过程环
节、改变过程方法、优化提高计算精度等方式对本发明不断完善和优化。一是标定姿控推力器推力环节，在推力器无需标定或推力器具备推力测定装置时，该过程可省略或直接利用已知值计算；二是确定卫星姿态矩阵环节，该过程还可通过改变计算姿态转换顺序、改变计算姿态中间矩阵、姿态描述方式等建立不同的姿态确定方法；三是提高姿态估计精度，通过其他姿态测量方法、误差修正方法、优化计算方法等方式提高姿态估计精度，比如通过对陀螺积分的修正提高姿态精度。此外在分析航天器受外力可能性研判方面，除了对航天器受力大小进行估算外，可在受力结果判定中根据半长轴的变化量级给出研判标准，将受力大小和轨道变化量级结合起来，全面地对航天器受外力情况准确判断，以提高方法应用的精准性。
[0025]
本发明突破了传统方法中按照航天器已知受力情况进行受力结果估算的思维，采取逆向推算方式，设计的基于瞬时受力冲击动量分析模型算法，从结果和过程出发，逆向溯源分析，为航天器异常分析和判断提供了新的思路，为在轨卫星受外力冲击后的结果认定提供研判方案，为分析撞击物的质量和大小提供了数据源，完善了航天器在轨状态感知方法，不仅对航天领域空间外力分析有重大创新，且本发明的计算结果为此类事件的定位与定性提供了有力证据。
附图说明
[0026]
图1是本发明算法流程图。
[0027]
图2是本发明算例姿控喷气平均脉宽示意图。
[0028]
图3是本发明算例本体系下各喷气节拍产生冲量示意图。
[0029]
图4是本发明算例惯性系下姿态四元数示意图。
[0030]
图5是本发明算例轨道系下各节拍产生冲量示意图。
具体实施方式
[0031]
下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。
[0032]
如图1所示，本发明在标定姿控推力器推力的基础上，计算喷气产生的冲量，将本体系下卫星喷气冲量与轨道系姿态结合后，计算轨道系产生的冲量及轨道变化量，与事后精密定轨结果差值以估算受外力冲击可能性。具体计算步骤如下：
[0033]
1.标定姿控推力器推力
[0034]
卫星姿控推力器大小由入口压力与喷气时长共同决定，因此应先确定当前条件下姿控推力大小曲线，确定不同喷气时长条件下姿控推力器推力大小。
[0035]
2.计算姿控异常喷气冲量
[0036]
令卫星推力器喷气脉宽序列为{tk}，相应脉宽的推力序列为{f
tk
}，各推力器的安装矩阵为m
bt
。
[0037]
计算各推力器在各控制节拍中的冲量ik＝f
tk
tk。
[0038]
如果推力器非正向安装时，根据各推力器在卫星本体坐标系下的安装方向，计算本体系下各推力器在各控制节拍中冲量矢量在本体坐标系下的矢量表示i
bk
＝m
bt
ik。
[0039]
3.确定卫星姿态矩阵
[0040]
卫星姿态异常喷气期间，恒星敏感器姿态角信息输出无效，为获得真实姿态信息，需要利用陀螺输出的角速度信息，积分计算卫星姿态。
[0041]
·
计算卫星瞬时姿态四元数。
[0042]
以卫星异常前一拍姿态测量数据为初值，令j2000惯性系下的姿态用四元数为q
gu0
＝[q
0 q
1 q
2 q3]，陀螺角速度向量表示为ω
bk
＝[ω
x ω
y ωz]
t
，扩展为4维向量为ω
bk
＝[ω
x ω
y ω
z 0]
t
。积分计算惯性系姿态四元数q
guk
＝[q
0 q
1 q
2 q3]。
[0043]
·
计算惯性系到本体系转换矩阵。
[0044]
根据惯性姿态四元数，计算得到惯性系到本体系转换矩阵m
bi
。
[0045][0046]
·
计算惯性系到轨道系转换矩阵。
[0047]
由卫星轨道坐标系定义可知，空间矢量由2000.0惯性坐标系至轨道坐标系的转换矩阵m
oi
为：
[0048][0049][0050]moi
(1,m)＝m
oi
(2,p)
×moi
(3,k)
[0051]
其中：i＝1,2,3对应于转换矩阵m
oi
中每个行向量的三个分量。m
oi
(2,p)和m
oi
(3,k)为矩阵m
oi
中第二行和第三行的行向量。表示为空间矢量r的单位矢量，空间矢量r的速度矢量。
[0052]
·
计算卫星本体系到卫星轨道坐标系的转换矩阵。
[0053]mob
＝m
oi
·
(m
bi
)
t
。
[0054]
4.计算姿控喷气轨道变化量
[0055]
·
计算第j个推力器产生的冲量
[0056]
以第j个推力器为例，喷气冲量计算为例，设第j个推力器在喷气时间内，共产生n个喷气脉宽。
[0057][0058]
·
计算卫星在运动方向的冲量
[0059]
假定按照有m个姿控推力器，则卫星在x方向的冲量i
ox
为：
[0060]
·
计算卫星半长轴变化量
[0061]
令卫星质量为m，卫星轨道运行速率为n，轨道偏心率e，真近点角为f，则卫星姿控喷气导致的半长轴变化量δag为：
[0062]
[0063]
5.分析卫星受外力冲击可能性
[0064]
利用异常发生前后的精密轨道确定数据，可以计算异常前后实际的半长轴变化量δa，得到差值δas＝δa-δag，该差值如果大于计算误差的3σ，即能够判断是受到空间外力冲击产生。
[0065]
下面以具体实例说明，基于瞬时受力冲击动量分析模型算法，计算某在轨卫星(以下称y星)的受力影响，通过计算结果数据对比，结合姿控喷气冲量分析后，可确定受外力冲击，并估算出外力冲击导致轨道半长轴变化量，为y星受外力冲击影响分析提供了有效分析方法。
[0066]
假定y星共在-x面安装6个推力器，姿控喷气会对卫星产生+x向速度增量，使半长轴抬高。
[0067]
通过建立姿控喷气模型，可计算喷气对卫星半长轴影响，此过程需考虑以下几个问题：
[0068]
(1)由于控制分系统在姿控喷气控制时喷气脉宽是不同的，而推力器在不同推力器入口压力与不同喷气脉宽条件下对应的推力也是不同的，因此应先确定当前条件下姿控控制推力的大小；
[0069]
(2)由于遥测参数设计及采样周期限制，没有直接参数表征推力器喷气时长及喷气脉宽，因此应考虑如何确定当前条件下推力器喷气时长及喷气脉宽；
[0070]
(3)用于姿控喷气的推力器存在安装角度区别；
[0071]
(4)卫星喷气期间，星敏输出无效，星上没有真实姿态信息。
[0072]
1.推力器推力确定
[0073]
(1)第一步：某压力下条件下，测量多台同型号推力器在不同工作脉宽下推力大小，以平均值作为基准。
[0074]
(2)第二步：考虑到当前卫星的实际推力器入口压力大小，插值计算出卫星当前推力大小与喷气脉宽的对应关系。结果如表1所示：
[0075]
表1推力器推力大小随喷气脉宽的变化
[0076][0077]
2.异常喷气冲量计算
[0078]
(1)第一步：根据卫星遥测数据分析，假定异常姿控喷气时长为119秒，参与喷气的推力器工作情况如表2所示。
[0079]
表2推力器喷气统计结果
[0080]
[0081]
(2)第二步：根据每帧遥测数据中喷气时长增量与喷气次数计算平均脉宽，其结果如图2所示。
[0082]
(3)第三步：根据推力标定结果得到实际推力大小，可得到各推力器在各控制节拍中的冲量。根据推力器安装布局，可计算得到各推力器在各控制节拍中冲量矢量在本体坐标系下冲量。计算结果如图3所示。
[0083]
3.基于陀螺数据的姿态矩阵确定
[0084]
根据陀螺测量数据，积分计算卫星姿态四元数，结果如图4所示。根据卫星姿态四元数即可得到卫星本体系到轨道系的转换矩阵。
[0085]
4.喷气导致轨道变化的估计
[0086]
根据上节计算的姿态转换矩阵，计算卫星在轨道坐标系下的各节拍喷气冲量，如图5所示。
[0087]
将上述计算得到的各推力器在各控制节拍中产生的冲量矢量在轨道坐标系x轴方向的分量i
x
进行累加，即可得到所有工作的推力器在整个控制过程中所产生的沿轨道坐标系x轴方向的冲量分量。
[0088]
经计算可得：i
x
＝-508.2ns。
[0089]
由冲量是负值，即冲量方向为沿-x轴方向，因此产生的轨道高度变化应为增大。根据目前卫星质量及当前轨道，计算可得喷气导致的半长轴变化量1034m。
[0090]
5.外力冲击可能性估计
[0091]
姿控喷气结束后，根据轨道确定结果，实际轨道抬高972米，低于计算理论半长轴62米。在不考虑推力器推力标定误差的情况下，其他误差分析主要包括：
[0092]
(1)卫星姿态确定存在误差。星上遥测没有真实姿态信息，地面通过陀螺积分计算卫星三轴姿态，受陀螺常值漂移、采样周期影响，姿态误差约1度左右，对轨道影响可以忽略。
[0093]
(2)遥测数据中，姿态采样周期为0.5秒，而喷气时长采样周期有两种，分别为2秒、10秒，对10秒采样周期数据进行均值处理，喷气时采用的姿态信息仍会存在偏差，姿态误差约在3度以内，对轨道影响影响较小。
[0094]
(3)考虑了喷气脉宽对推力大小的影响，对姿控喷气脉宽进行了计算，并引入对整星冲量计算，在不考虑推力模型误差，推力误差很小，可以忽略。
[0095]
(4)喷气初期，催化床温度偏低，对喷气效率有一定影响，计算中对前10秒推力效率按80％预估，理论抬高轨道值降低约16米。
[0096]
综上所述，实际轨道比理论计算轨道降低62米以上，超过理论误差的3σ，可判断卫星受空间瞬时外力冲击的可能性较大。
技术特征：1.一种识别卫星空间瞬时受力冲击的动量分析方法，其特征在于，包括以下步骤：(1)确定当前条件下姿控推力大小曲线，确定不同喷气时长条件下姿控推力器推力大小；(2)令卫星推力器喷气脉宽序列为{t
k
}，相应脉宽的推力序列为{f
tk
}，各推力器的安装矩阵为m
bt
，计算各推力器在各控制节拍中的冲量i
k
＝f
tk
t
k
；如果推力器非正向安装，根据各推力器在卫星本体坐标系下的安装方向，计算卫星本体坐标系下各推力器在各控制节拍中冲量矢量在卫星本体坐标系下的矢量表示i
bk
＝m
bt
i
k
；(3)计算卫星本体系到卫星轨道坐标系的转换矩阵m
ob
＝m
oi
·
(m
bi
)
t
，其中，m
bi
表示惯性系到本体系的转换矩阵，m
oi
表示2000.0惯性坐标系至轨道坐标系的转换矩阵；(4)根据第j个推力器产生的冲量以及卫星在运动方向的冲量，计算卫星半长轴变化量δa
g
；(5)利用异常发生前后的精密轨道确定数据，计算异常前后实际的半长轴变化量δa，得到差值δa
s
＝δa-δa
g
，该差值如果大于设定阈值，即判定是受到空间外力冲击产生。2.根据权利要求1所述的识别卫星空间瞬时受力冲击的动量分析方法，其特征在于，所述的步骤(3)包括以下步骤：以卫星异常前一拍姿态测量数据为初值，令j2000惯性系下的姿态用四元数为q
gu0
，将陀螺角速度向量ω
bk
扩展为4维向量后，积分计算惯性系姿态四元数q
guk
＝[q
0 q
1 q
2 q3]；根据惯性姿态四元数，计算得到惯性系到本体系转换矩阵m
bi
；由卫星轨道坐标系定义可知空间矢量由2000.0惯性坐标系至轨道坐标系的转换矩阵m
oi
；计算得到卫星本体系到卫星轨道坐标系的转换矩阵m
ob
＝m
oi
·
(m
bi
)
t
。3.根据权利要求1所述的识别卫星空间瞬时受力冲击的动量分析方法，其特征在于，所述的惯性系到本体系转换矩阵4.根据权利要求1所述的识别卫星空间瞬时受力冲击的动量分析方法，其特征在于，所述的2000.0惯性坐标系至轨道坐标系的转换矩阵m
oi
为：为：m
oi
(1,m)＝m
oi
(2,p)
×
m
oi
(3,k)其中，i＝1,2,3对应于转换矩阵m
oi
中每个行向量的三个分量，m
oi
(2,p)和m
oi
(3,k)为矩阵m
oi
中第二行和第三行的行向量，表示为空间矢量r的单位矢量，空间矢量r的速度矢量。5.根据权利要求1所述的识别卫星空间瞬时受力冲击的动量分析方法，其特征在于，所述的步骤(4)中，第j个推力器在轨道系中产生的冲量n表示第j个推力器在喷气时间内共产生n个喷气脉宽，i
bji
表示本体系下第j个推力器的第i个脉宽产生的冲
量。6.根据权利要求1所述的识别卫星空间瞬时受力冲击的动量分析方法，其特征在于，所述的步骤(4)中，假定某卫星有m个姿控推力器，则各方向推力为m个姿控推力器在各方向产生的分量之和，在轨道系中前进方向t方向的冲量7.根据权利要求1所述的识别卫星空间瞬时受力冲击的动量分析方法，其特征在于，所述的步骤(4)中，令卫星质量为m，卫星轨道运行速率为n，轨道偏心率e，真近点角为f，则卫星姿控喷气导致的半长轴变化量i
ox
表示卫星在x方向的冲量。8.根据权利要求1所述的识别卫星空间瞬时受力冲击的动量分析方法，其特征在于，所述的步骤(5)中，差值δa
s
＝δa-δa
g
如果大于计算误差的3σ，即判定是受到空间外力冲击产生。
技术总结本发明提供了一种识别卫星空间瞬时受力冲击的动量分析方法，确定当前条件下姿控推力大小曲线，确定不同喷气时长条件下姿控推力器推力大小；依次计算各推力器在各控制节拍中的冲量、卫星本体系到卫星轨道坐标系的转换矩阵和卫星半长轴变化量；利用异常发生前后的精密轨道确定数据，计算异常前后实际的半长轴变化量，该差值如果大于设定阈值，即判定是受到空间外力冲击产生。本发明能够为卫星姿态异常定位、影响域分析、溯源分析及对撞击物质量与大小评估提供直接数据支撑。小评估提供直接数据支撑。小评估提供直接数据支撑。
技术研发人员：刘凯 孙先伟 胡兴 陈军 蔡立锋 饶斌 宁小娟 高燕 罗亮 刘伟 梁欣
受保护的技术使用者：中国西安卫星测控中心
技术研发日：2022.03.21
技术公布日：2022/7/5
