一种高精度的航天器辐射光压建模方法

1.本发明属于航空航天领域，尤其涉及一种高精度的航天器辐射光压建模方法。


背景技术：

2.辐射光压作为一种非保守力，是由光子与物体表面碰撞时的动量转移所产生的力学效应。对于地球轨道航天器而言，来自太阳、地球以及航天器自身的辐射光子流会与航天器的表面相互作用，从而产生辐射光压。虽然辐射光压的量级较小，但会对长期在轨运行的航天器轨道产生较大影响。与建模技术相对成熟且精度较高的地球引力、日月引力等相比，由于辐射强度的变化、航天器姿态控制的偏差以及表面材料的老化等因素的影响，高准确度的辐射光压建模较难实现，目前已经成为航天器动力学定轨的主要误差源。辐射光压的物理机制表明，不同航天器的辐射光压模型具有特异性，不可能建立一个适用于所有航天器的辐射光压模型，需要对不同的航天器建立不同的模型，从而大大增加了建模工作的复杂性。
3.随着越来越多航天任务对航天器在特定坐标框架下空间位置的高精度要求，建立高精确的航天器轨道力学模型成为工程需要。以全球导航卫星(gnss，global navigation satellite system)为例，利用其实现高精度定位需要高精度的卫星位置。辐射光压作为目前中高轨卫星最主要的误差源，其精确度直接影响卫星的精密定轨与轨道预报。对于执行高精度地球观测的低轨道卫星而言，辐射光压是其轨道力学的第二大误差来源，因此也需要对其进行高精度建模。
4.按照建模原理的不同，现有的辐射光压建模方法大致可以分为三类：(1)不使用任何物理信息的纯经验方法；(2)考虑部分几何结构和表面材料信息的半经验方法；(3)使用所有可用的几何形状、表面材料、姿态、空间环境等信息的纯物理方法。其中，以ecom(extended code orbit model)和gspm(gps solar pressure model)为代表的纯经验方法模型大都为截断傅里叶级数形式(常数项、三角函数项)，需要大量的观测数据和长时间的数据拟合才能达到较高精度，难以用于刚发射的卫星，而且模型参数与待估计轨道参数具有较高相关性。半经验方法虽然在一定程度上减少了所需数据，但是模型中参数的设定(参数化)也需要仔细斟酌，避免与所估计参数产生强相关性。纯物理方法模型利用一切可用的先验信息，更能反映物理本质，便于分析各物理效应对辐射光压的影响，而且只要模型输入参数的准确性得到保障，基于物理规律的辐射光压模型也具有高精度。
5.目前纯物理建模方法的实现方式主要为tpmc(test particle monte carlo)和射线追踪方法。但上述两种方法均存在计算量大，计算耗时长的缺陷，难以满足海量航天器的快速高精度辐射光压建模需求。因此，研究适用性强、精度高的通用航天器辐射光压物理建模方法是十分必要的。


技术实现要素：

6.为解决上述现有技术存在的不足，本发明提出一种适用性强、精度高的航天器辐
射光压物理建模方法，该建模方法将辐射光子流量、航天器三维结构、以及辐射光子与航天器表面的相互作用作为基本要素，综合考虑了航天器形状和材料、多次反射等因素，并对太阳辐射光压、地球辐射光压、辐射光压和天线辐射光压进行了整合，降低了建模的整体复杂度，提高了模型的精度。
7.本发明的技术方案如下：
8.一种高精度的航天器辐射光压建模方法，包括以下步骤：
9.s1：建立辐射源模型
10.计算辐射源的量级大小和方向，所述辐射源包括太阳辐射、地球辐射和航天器热辐射和航天器天线辐射；
11.s2：构建航天器三维模型数据库
12.将航天器模型各组件对应于基本几何原型，将基本几何原型组合成复杂航天器三维模型，以参数化的形式描述各基本几何原型，结合航天器模型各组件的表面材料特性，建立航天器三维模型数据库；
13.s3：模拟辐射光子与航天器表面的相互作用
14.采用射线追踪算法，计算辐射源产生的光子与航天器表面作用产生的力的大小和方向，得到航天器辐射光压模型。
15.优选的，s1中，所述太阳辐射建模中考虑太阳到航天器的距离和航天器处于地影中的蚀因子；所述地球辐射建模中考虑地球的长波辐射和短波辐射；所述航天器热辐射建模考虑航天器表面吸收的太阳辐射和地球辐射转换的热辐射，以及航天器内部仪器运行所产生的热量。
16.优选的，s2中还包括，利用航天器三维模型数据库，构建一个kd树结构存储三维模型。
17.优选的，s1具体包括：
18.依据输入总太阳辐射、地球辐射数据，计算航天器所在位置的太阳辐射功率、地球辐射功率以及航天器热辐射功率；
19.其中太阳辐射功率的计算公式如下：
[0020][0021]
其中f
s
为根据阴影函数模型计算航天器当前位置的蚀因子函数，在0
‑
1之间取值，e
s
为输入的总太阳辐射(e
s
＝1361)，au为天文单位，r为航天器的位置，为航天器的速度，r
s
为根据行星星历计算太阳的惯性系坐标；
[0022]
其中地球辐射的计算方法如下：
[0023]
将地球表面进行三角格网化，对于长波地球辐射，如当前地球表面三角网格对于航天器可见，则长波地球辐射可以直接到达航天器所在位置；对于短波地球辐射，当前地球表面三角网格对于航天器和太阳同时可见，则短波地球辐射可以到达航天器所在位置，根据以上长波地球辐射和短波地球辐射的条件，遍历所有地球表面三角格网，计算航天器所在位置的长波辐射功率w
s，l
和短波辐射功率w
s，s
。
[0024]
其中航天器热辐射功率的计算方法如下：
[0025]
计算航天器表面的温度t
mli
：
[0026][0027]
其中，α为航天器表面材料的吸收率，θ为外部辐射方向与航天器表面法向的夹角，ε
eff
为航天器内部到外太空的多层隔热材料的有效热发射率，σ为史蒂芬
‑
玻尔兹曼常数，t
sc
为航天器内部温度，ε
mli
为航天器表面材料的发射率；
[0028]
进一步计算多层隔热材料的热辐射功率w
t
：
[0029][0030]
其中，n为航天器表面材料的法向量；
[0031]
在没有外部辐射的条件下，多层隔热材料的热辐射功率为：
[0032][0033]
优选的，所述s2中的基本几何原型包括平面、抛物面、圆、球、圆环、圆锥、圆柱和圆台；所述航天器表面材料特性包括：反射率、镜面反射率、发射率和透射率以及这些特性参数的测量标准差。
[0034]
优选的，s2中所述构建kd树结构存储三维模型的方法为：
[0035]
首先，构建一个轴对齐的矩形边界框(aabb)，即航天器的最大几何包络，将aabb作为树结构的根节点，然后设置划分平面将aabb在空间上划分为2个部分，并分别识别这2个部分里的航天器三维模型组件，将这2个部分作为树结构的2个叶子节点；重复上述过程直到空间划分里只有1个或者0个几何原型，则kd树结构构建结束。
[0036]
优选的，s3具体包括以下步骤：
[0037]
s31：对每一束光线执行射线追踪算法，利用kd树结构计算辐射光与航天器三维模型各组件之间的交点，计算交点上该束辐射光的太阳辐射光压、地球辐射光压、航天器热辐射光压以及航天器天线辐射光压；
[0038]
其中太阳辐射光压的计算方法如下：
[0039][0040]
其中，a为被辐射光照射的航天器组件表面面积，θ为入射辐射流量的入射角，c为光速，s为入射辐射单位向量，μ为航天器表面材料的镜面反射率；v为航天器表面材料的反射率，p为镜面反射向量；
[0041]
其中地球辐射光压的计算方法如下：
[0042][0043]
其中航天器表面热辐射光压的计算方法如下：
[0044][0045]
其中航天器天线辐射光压的计算方法如下：
[0046][0047]
其中，w
a
为航天器天线的辐射功率；
[0048]
s32：对于多次反射效应的计算，辐射光与航天器表面上次碰撞的交点作为下一次反射的起点，镜面反射方向作为下一次反射的光线方向，且反射的能量根据航天器表面材料的吸收率α进行调整，追踪每一束辐射光的反射直到该辐射光进入太空或者能量衰减为0；对于每一次反射，均采用s31中的方法进行辐射光压计算；
[0049]
s33：将所有辐射光的计算结果累加起来，即得当前辐射源在当前方向的太阳辐射光压、地球辐射光压、航天器热辐射光压以及航天器天线辐射光压；
[0050]
s34：将航天器周围的空间按照经度和纬度划分为181
×
361个格网，其中纬度有181个方向，经度有361个方向，对于每一个辐射源方向，均采用s31
‑
s33中的方法计算出其对应的辐射光压，并按照星固系的x，y和z分量生成3个格网文件，每个格网文件都存储181
×
361个数据。
[0051]
优选的，所述航天器辐射光压建模方法还包括s4：
[0052]
s41：用户根据实际情况计算所需要的辐射流量的经度、纬度和大小；
[0053]
s42：根据辐射流量的经度和纬度在s3所得格网文件中进行双线性插值得到相应的航天器姿态和辐射光压，然后根据实际的辐射流量大小进行缩放，得到真实的辐射光压。
[0054]
优选的，用户端使用双线性内插得到需要的航天器姿态数据，如下式：
[0055][0056]
其中f
b
为航天器所受辐射光压，m0为建模中航天器的质量，m为实际航天器质量；w为航天器实际受到的辐射功率；e
s
为建模中采用的太阳辐射功率；和λ分别为辐射源在航天器星固坐标系下的纬度和经度；分别为在当前经纬度下插值所得到的航天器星固坐标系3个方向的辐射光压加速度，即太阳辐射光压、地球辐射光压、航天器热辐射光压和航天器天线辐射光压。
[0057]
优选的，s2中航天器三维模型数据库还存储有航天器表面材料的光学和热学特性的不确定度，则s3中计算辐射源产生的光子与航天器表面作用产生的力的不确定度如下：
[0058][0059]
其中
[0060][0061]
其中，w
i
为入射辐射功率，f
i
(v，μ)表示第i束辐射光的射线追踪所计算的辐射光压力；和分别表示镜面反射率和反射率的测量方差；δ
xy
，δ
xz
，δ
yz
分别表示
该辐射光的辐射光压力在x，y和z方向的方差以及相应的协方差。
[0062]
相比于现有技术，本发明具有以下优势：
[0063]
1.本发明提出的高精度的航天器辐射光压建模方法能够建立高精度的辐射光压模型，且适用性强，可灵活用于不同的航天器建模。同时，作为一种非保守力建模方法，能够完善并提高航天器动力学模型，进一步提高精密定轨和轨道预报的精度。
[0064]
2.本发明通过对辐射源建模、航天器三维模型、以及辐射与航天器表面的相互作用这三个基本要素进行分步建模的方法，实现了对不同辐射源(太阳辐射光压、地球辐射光压、热辐射光压、天线推力)的统一整合建模，而并非是分开单独建模，降低了建模的复杂度，提高了精确度。
[0065]
3.本发明使用的参数化几何建模方法，通过几何原型(比如圆、多边形、圆锥、球、抛物面等)的拼接组合，更加精确地描述航天器的几何形状(相对于格网化方式)；综合几何形状与表面材料特性构建了统一的航天器三维模型数据库，规范化了用于航天器辐射光压建模的数据库格式，便于计算和管理。
[0066]
4.本发明使用射线跟踪法，可以得到辐射光任意多次反射后的结果；还可以结合gpu实现并行计算，大大提高计算速度。此外，本发明可以利用输入的表面材料光学特性参数计算辐射光压的不确定度。
[0067]
5.本发明采用kd树结构重新组织航天器三维模型，大大提高了射线追踪算法的计算效率；而且提出了基于输入表面材料光学特性参数的不确定性计算辐射光压不确定性的方法。
[0068]
6.本发明采用格网文件存储太阳辐射光压、地球辐射光压、热辐射光压等计算结果，能够适用于航天器的任何姿态。
附图说明
[0069]
图1为发明的高精度的航天器辐射光压建模方法的流程图；
[0070]
图2为本发明的地球辐射建模中长波和短波到达卫星的条件；
[0071]
图3为jason
‑
2海洋测高卫星的三维模型与aabb；
[0072]
图4为射线跟踪法基本流程示意图；
[0073]
图5为辐射光压产生的示意图。
具体实施方式
[0074]
下面结合附图，对本发明的具体实施方式进行详细描述，但应当理解本发明的保护范围并不受具体实施方式的限制。
[0075]
实施例1
[0076]
以jason
‑
2海洋测高卫星为例，对本发明提供的高精度的航天器辐射光压建模方法进行说明，具体包括如下步骤：
[0077]
s1：根据卫星和太阳在惯性系的位置，计算到达卫星位置的太阳辐射，以及地球短波辐射和长波辐射；其中长波辐射和短波辐射要分别计算，图2展示了地球表面的三角格网上计算长波辐射和短波辐射的过程，其中短波辐射只有在三角格网同时对卫星和航天器可见的时候才可以到达卫星，而长波辐射只要三角格网对卫星可见便可以到达卫星。
[0078]
s2：根据卫星实际情况构建三维数据库，需要包含卫星的三维几何原型以及其表面材料的反射率、镜面反射率、热发射率以及相应的标准差。图3展示了用本发明方法所构建的jason
‑
2卫星的星体结构。
[0079]
s3：根据卫星三维数据库中存储的卫星几何原型，提取出每个几何组件的8个顶点坐标，遍历所有的组件，从而构建整个卫星外形的aabb，并在此基础上构建整个卫星三维模型的kd树结构。
[0080]
s4：将整星的aabb按照辐射源的方向投影到平面上，其中投影平面的长度和宽度根据整星的aabb来确定，并在投影平面上建立像素阵列空间坐标系的x和y轴，辐射源到卫星质心的指向作为像素阵列空间坐标系的z轴。
[0081]
s5：在所建立的像素阵列平面上，按照所设置的空间分辨率要求，将像素平面划分为像素大小的格网。像素的尺寸可灵活设置，一般设置为1mm或者1cm。
[0082]
s6：在已建立好的像素阵列上，每一个像素的坐标作为一束辐射光的起点，辐射源到卫星质心的连线作为光的方向；对于每一束辐射光通过并行的方式启动射线追踪过程：
[0083]
s61：射线追踪计算中需要遍历卫星三维模型数据库中各个组件，由于已经将各个卫星组件放置于kd树结构中，对每个卫星组件的访问和遍历本质上是kd树的遍历，该遍历方法简洁快速，可以大大提高射线追踪计算效率。
[0084]
s62：对于某一束辐射光，利用kd树的遍历找到与这束光相交的交点，在交点上利用辐射光压的计算公式计算太阳辐射、星体热辐射、地球辐射光压以及卫星天线辐射，然后将镜面反射光线作为下一次反射的入射光，开始寻找下一次反射的交点；如此重复直到这一束辐射光不与卫星的任何部分相交或者光的能量在反射中消耗为0。图5展示了一束辐射光与平面的相互作用过程。整个射线追踪的计算过程如图4所示。
[0085]
s7：将所有辐射光的计算结果累加起来得到当前辐射源在当前方向的太阳辐射光压、地球辐射光压、星体热辐射光压以及卫星天线辐射光压。
[0086]
s8：在卫星的星体坐标系中，按照经纬度间隔为1度来设置辐射源的方向，共得181
×
361个方向；对于每一个方向重复步骤s4到s7计算，得到每一个辐射源方向辐射光压的x，y以及z分量，将x，y和z分量存放在3个格网文件中，每个格网文件都存储181
×
361个数据。
[0087]
根据以上步骤得到经过建模计算的x，y和z这3个方向的格网文件，用户可以根据格网文件计算任意卫星姿态、任意时刻的辐射光压加速度。用户算法步骤如下：
[0088]
s9：用户根据卫星的位置以及辐射源的种类(太阳辐射或者地球辐射)计算出辐射方向及量级；根据辐射方向利用提供的格网文件分别进行双线性插值得到x，y和z方向的辐射光压加速度；
[0089]
s10：根据实际的辐射功率，对以上插值得到的辐射光压及速度进行缩放；得到真实的辐射光压加速度。
[0090]
综上，本发明针对现有辐射光压建模中的不同辐射源独立建模、参数相关性高、建模工作复杂、难以拓展应用到各种航天器等问题，提出了一种适用性强、精度高的航天器辐射光压建模方法。该方法在提高了模型精度的同时，整合太阳辐射光压、地球辐射光压和热辐射压为一个统一模型，降低了模型的复杂度，采用几何原型精细化描述航天器模型，规范化了航天器辐射光压建模的数据库格式，便于统一管理；在射线跟踪法中采用并行计算大大提高了计算速率；输出结果为格网文件，可适用于不同类型的航天器的不同姿态，提高了
动力学定轨和轨道预报的精确度；而且可以通过表面材料的光学特性参数评估结果的不确定度。本发明提出的建模方法在航天器精密定轨以及导航卫星高精度定位中，有着广阔的应用前景。
[0091]
以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明。应当指出，对于本技术领域的普通技术人员来说，在不脱离发明原理的前提下，还可以对发明进行若干改进和修饰，这些改进和修饰也落入发明权利要求的保护范围内。因此，本发明的保护范围应所述以权利要求的保护范围为准。