一种复杂外形航天器投影面积的确定方法及系统与流程

本发明涉及一种复杂外形航天器投影面积的确定方法及系统，属于飞行器气动特性设计和姿态轨道控制技术领域。

背景技术：

近年来，由于重力场和稳态海洋环流精确测量的需要，超低地球轨道(200-500km)航天器逐渐展现出广阔的应用前景，然而，实现重力梯度精确测量的先决条件是低轨卫星气动阻力的准确计算。气动阻力主要由当地大气密度、阻力系数和卫星沿运动方向的投影面积共同决定。当地大气密度一般由大气模型给出，气动阻力系数可以查阅数据库，还需要解决的问题就是计算其沿运动方向的投影面积。

最近，太阳帆探测器在深空探测领域显示出了巨大的潜力，因为其利用太阳光压力获取推动力，不需要携带大量推进剂。虽然太阳光压很小，但持续地加速可使探测器达到一个可观的速度，最终可达到传统航天器5-10倍的速度。太阳帆探测器的姿态和轨道控制系统需要太阳光压力的输入，太阳光压力则与太阳通量、太阳帆表面反射率和沿光线方向的投影面积有关。太阳通量为物理常数，反射率为材料属性，可通过查阅相关手册或实验测量获得，太阳光压力的计算就归结为太阳帆沿光线方向投影面积的计算。

此外，航天器防护系统设计和风险评估中重要的一环是微流星体或空间碎片与航天器部件撞击概率的预估，其前提是提供航天器及其各个部件在给定威胁攻击方向上的投影面积。同时，卫星沿某一方向的投影面积更直接是其姿态控制的必要输入参数，比如太阳瞄准(sun-pointing)姿态模式需要输入卫星沿太阳光线方向的投影面积，最小或最大阻力(minimum-dragormaximum-drag)姿态模式则需要识别并计算出卫星的最小或最大投影面积。

生成几何实体沿某方向的投影，在计算机图形学领域称为遮挡算法，是计算几何实体投影面积的一种直观数值方法，至今仍在不断的研究和完善之中。目前，计算机图形学广泛使用的投影算法是光线投射和光线追踪方法：光线投射方法计算从光源发射出的光线与几何实体表面的交点，通过最小距离识别出可见面(非可见面为背光面)，然后采用多边形投影算法对可见面进行投影即可生成几何实体的投影图像；光线追踪方法在光线投射方法的基础上，考虑了光线散射和反射效应。

然而，无论是低轨卫星大气阻力、太阳帆探测器光压的计算，还是航天器防护设计与风险评估，以及航天器姿态控制，都只需要航天器(部件)投影面积的输入，通过计算机图形学的方法生成投影图像再计算投影面积显然耗时耗力。国外有研究人员基于凸多边形理论，发展了一种计算航天器沿任意方向投影面积的解析方法。由于每个过程均采用解析解的形式来实现，该方法的计算效率很高。然而，复杂外形的航天器进行多边形离散后，其投影多边形难免存在两两相交甚至三个或三个以上多边形同时相交的情形，导致计算投影多边形相交部分的面积极端困难，甚至不可能。而且，解析的特性导致其鲁棒性和通用性较差，某些极端情况下，计算机固有的舍入误差会导致错误的计算结果。航天器外形越复杂，出现极端情况的可能性越大。因此，上述解析方法虽然高效，却很难应用于航天工程实践中。

技术实现要素：

本发明的技术解决问题是：克服现有技术的不足，提出了一种复杂外形航天器沿任意方向投影面积的确定方法及系统，解决了现有技术中无法计算复杂外形航天器投影面积如的问题。

本发明的技术解决方案是：

一种复杂外形航天器投影面积的确定方法，包括以下步骤：

(1)构建投影坐标系oxyz；

(2)构建一个包围航天器外表面的长方体盒子，长方体的3条相互垂直的棱分别平行于x、y和z轴；

(3)确定产生试验粒子的边界矩形；

(4)在步骤(3)给出的边界矩形中随机产生试验粒子，给出试验粒子在投影坐标系中初始坐标r0；

(5)判断试验粒子射线运动轨迹是否与航天器表面三角形单元相交：如果没有相交，直接转到步骤(6)；如果相交，则令相交试验粒子数q递增1，然后继续执行步骤(6)；

(6)判断产生的试验粒子总数m是否小于预设值，如果小于预设值，转到步骤(4)，否则继续执行步骤(7)；

(7)根据所述相交试验粒子数q和产生的试验粒子总数m，统计航天器的投影面积。

构建投影坐标系oxyz，具体为：

(1.1)首先构建机体坐标系oxyz，坐标原点即为航天器质心o，x、y和z轴与航天器惯性主轴重合；

(1.2)根据机体坐标系oxyz，构建投影坐标系oxyz，投影坐标系原点位于航天器质心o，z轴正方向为投影方向，令oz在平面的投影为ozp，则oy位于xoy平面内且垂直于ozp，ox通过右手定则确定。

所述步骤(2)构建一个包围航天器外表面的长方体盒子，该长方体的6个面在投影坐标系表示为

航天器表面由三角形单元离散，n为航天器表面三角形单元数量，(xi,yi,zi),i∈[1,n]为三角形单元节点在投影坐标系中的坐标。

所述步骤(3)边界矩形垂直于投影方向，且位于所述长方体盒子的上游，其坐标满足

步骤(4)中初始坐标r0(x0,y0,z0)具体为：

式中，r1、r2均为(0,1)区间均匀分布的随机数。

所述步骤(5)中判断试验粒子射线运动轨迹是否与航天器表面三角形单元相交，具体包括以下步骤：

(5.1)计算试验粒子射线运动轨迹：根据试验粒子初始位置坐标r0(x0,y0,z0)和投影方向，得射线参数方程为

其中，t是时间参数；

(5.2)计算三角形单元所在平面的方程：令三角形abc任意一个顶点的坐标为(x1,y1,z1)，表面外法向矢量为n＝(n1,n2,n3)，则其所在平面的方程为

n·(x-x1,y-y1,z-z1)＝0；

(5.3)判断试验粒子射线轨迹是否与三角形单元所在平面相交：联立射线参数方程和三角形所在平面的方程

若n3＝0，则平面平行于射线，两者不相交，否则求得

若t≤0，则射线与平面不相交，进而不会与三角形单元相交，若t＞0，则进入

步骤(5.4)进一步判断交点是否在三角形abc内部；

(5.4)代入参数t的值到射线参数方程，得交点p的坐标(xp,yp,zp)为

令

若同时满足

则交点p在三角形abc内部，试验粒子射线运动轨迹与航天器表面三角形单元相交，否则交点在三角形abc外部，射线运动轨迹与三角形单元不相交。

所述步骤(6)中的预设值大于等于10⁶。

所述步骤(7)中统计航天器的投影面积ap，具体公式为

ap＝(xmax-xmin)(ymax-ymin)(q/m)

式中，m为产生并跟踪的试验粒子总数，q为运动轨迹与航天器表面相交的试验粒子数。

所述预设值的取值范围10⁶～10⁸。

一种复杂外形航天器投影面积的确定系统，包括：

坐标系构建模块：用于构建投影坐标系oxyz；

盒子构建模块：用于构建一个包围航天器外表面的长方体盒子，长方体的3条棱分别平行于x、y和z轴；

边界矩形产生模块：用于确定产生试验粒子的边界矩形；

试验粒子产生模块：用于在边界矩形产生模块确定的边界矩形中随机产生试验粒子，给出试验粒子在投影坐标系中初始坐标r0；

相交判断模块：用于判断试验粒子射线运动轨迹是否与航天器表面三角形单元相交，并统计相交试验粒子数q；

投影面积计算模块：用于根据相交判断模块统计出的相交试验粒子数q和产生的试验粒子总数m，统计航天器的投影面积。

本发明与现有技术相比的优点：

(1)当前复杂外形航天器投影面积的计算方法，要么通过计算机图形学的方法生成投影图像再计算投影面积显得耗时耗力，要么基于凸多边形理论采用解析的方法导致通用性和鲁棒性较差而无法应用于航天工程。基于此，本发明借鉴统计学中montecarlo技术的思想，结合计算几何学中基本的射线与三角形相交理论，发展了一种复杂外形航天器沿任意方向投影面积的计算方法。

(2)相对于计算机图形学的投影算法，本发明方法计算速度快，存储要求低，相对于基于凸多边形理论的解析方法而言，本发明方法鲁棒性和通用性好，具备处理任意复杂工程外形的能力。此外，本发明方法计算时间随着计算规模、问题复杂程度和空间维数的增大呈线性关系，是复杂外形航天器投影面积计算的理想方法。

附图说明：

图1：本发明方法流程图；

图2：机体坐标系和投影坐标系示意图；

图3：单面圆形平板示意图(a)及其三角形单元离散化(b)；

图4：单面圆形平板计算结果与理论结果对比图；

图5：samson卫星模型示意图；

图6：samson卫星投影面积随着其姿态的变化，其中，图6(a)为三维分布图，图6(b)为等值曲线；

图7：goce卫星模型示意图，其中图7(a)为头部方向视图，图7(b)为尾部方向视图；

图8：goce卫星投影面积随着其姿态的变化；其中，图8(a)为三维分布图，图8(b)为等值曲线。

具体实施方式

根据现有复杂外形航天器投影面积计算的需要，借鉴统计学中montecarlo技术的思想，结合计算几何学中基本的射线与三角形相交理论，本发明提出了一种复杂外形航天器沿任意方向投影面积的计算方法。

首先，构建一个恰好包围航天器外表面的长方体盒子，该长方体有一个表面垂直于投影方向；然后，确定产生试验粒子的边界矩形，边界矩形所在平面垂直于投影方向，且位于包围盒子的上游；再次，在边界矩形产生试验粒子，给出其位置坐标，并跟踪该试验粒子的运动轨迹，判断轨迹是否与航天器表面三角形单元相交，若相交则使相交试验粒子数递增1，否则进行下一步；最后，重复上一步过程直至试验粒子总数足够大，以保证计算结果收敛，根据相交试验粒子数和产生的试验粒子总数统计出航天器的投影面积。

如图1所示，具体步骤为：

(1)构建机体坐标系b和投影坐标系p，具体实现过程如下：

机体坐标系b的坐标原点位于航天器质心o，x、y和z轴与航天器惯性主轴重合。投影坐标系p的原点同样位于航天器质心o，其z轴与投影方向重合，即投影方向为oz轴，令oz在平面的投影为ozp，则oy位于xoy平面内且垂直于ozp，ox通过右手定则确定，

如图2所示。采用航天器轨道力学或天文学的术语，oz与xoy平面的倾角定义为赤纬δ，其范围为δ∈[-π/2,π/2]，ozp与ox的夹角定义为赤经α，其范围为α∈[0,2π)。若非特殊说明，下文中的坐标或者矢量均为投影坐标系p中的量。

(2)构建一个恰好包围航天器外表面的长方体盒子，长方体的3条相互垂直的棱分别平行于投影坐标系的x、y和z轴。通常，复杂外形航天器表面由三角形单元离散，设表面三角形单元数量为n，单元节点在投影坐标系中的坐标为(xi,yi,zi),i∈[1,n]，则该长方体盒子的6个面由下式确定

(3)确定产生试验粒子的边界矩形。边界矩形所在平面垂直于投影方向，且位于包围盒子的上游，其坐标满足

(4)在步骤(3)给出的边界矩形随机产生试验粒子，给出试验粒子在投影坐标系中位置，初始坐标r0(x0,y0,z0)为

r1、r2均为(0,1)区间均匀分布的随机数。

(5)判断试验粒子射线运动轨迹是否与航天器表面三角形单元相交：如果没有相交，直接转到步骤(6)；如果相交，则令相交试验粒子数q递增1，然后继续执行步骤(6)。q的初始值为0；具体实现过程如下：

(5.1)计算试验粒子射线运动轨迹。根据试验粒子初始位置坐标r0(x0,y0,z0)和投影方向，容易得其参数方程为

(5.2)计算三角形单元所在平面的方程。令三角形abc任意一个顶点的坐标为(x1,y1,z1)，表面外法向矢量为n＝(n1,n2,n3)，则其所在平面的方程为

n·(x-x1,y-y1,z-z1)＝0

(5.3)判断试验粒子射线轨迹是否与三角形单元所在平面相交。联立射线参数方程和三角形所在平面的方程

若n3＝0，则平面平行于射线，两者显然不相交，否则容易求得

若t≤0，则射线与平面不相交，自然不会与三角形单元相交，否则需要进一步判断交点是否在三角形abc内部。

(5.4)若试验粒子射线运动轨迹与三角形abc所在平面相交，需要判断交点p是否在三角形abc内部。代入参数t的值到射线参数方程，容易得交点p的坐标(xp,yp,zp)为

令

若同时满足

则交点p在三角形abc内部，试验粒子射线运动轨迹与航天器表面三角形单元相交，否则交点在三角形abc外部，射线运动轨迹与三角形单元不相交。

(6)判断产生和跟踪的试验粒子总数是否小于预设值，如果小于预设值，转到步骤(4)，否则继续执行步骤(7)。

一般跟踪的试验粒子的数量大于等于10⁶时，开始统计投影面积。根据统计误差和计算时间情况综合考虑，可以在10⁶～10⁸内取值。

(7)统计航天器的投影面积。设产生并跟踪的试验粒子数为m，其中，运动轨迹与航天器表面相交的试验粒子数为q，则航天器沿z轴的投影面积ap为

ap＝(xmax-xmin)(ymax-ymin)(q/m)

统计误差ε为单一误差分布标准差σ与试验粒子总数平方根的商，即

所以，产生并跟踪足够多的试验粒子之后，投影面积一定会收敛，且统计误差为

复杂外形航天器投影面积计算的具体求解实例如下：

首先考虑考虑外形简单、投影面积存在解析解的单面圆形平板，便于对比分析。单面圆形平板示意图、机体坐标系和投影方向如图3a所示，轴向、展向和法向分别为x、y和z轴，投影轴为z轴。圆形平板半径为1m，由于轴对称特性，只需要考虑赤经α＝0°的情形。因此，z轴位于xoz平面，与x轴的夹角为赤纬δ，与z轴的夹角为赤纬的余角。采用608个三角形单元进行离散，如图3b所示，试验粒子总数为2×10⁷，此时统计误差为

图4给出了该单面平板投影面积随赤纬的变化，赤纬的变化范围是δ∈[0°,90°]，即投影方向从平行于平板表面变化到垂直于平板表面。投影方向平行于平板表面时(δ＝90°)，投影面积为0；投影方向垂直于平板表面时(δ＝90°)，投影面积为平板面积π；随着投影方向从平行于平板表面变化到垂直于平板表面时，投影面积呈正弦函数增大。

对于该简单外形，投影面积的理论解为平板面积与赤纬正弦的乘积，即

ap＝πcos(π/2-δ)＝πsinδ

显然，本发明montecarlo统计结果与理论解完全一致，验证了方法的可靠性。

然后，考虑欧洲空间局的samson小卫星，如图5所示。表1列出了samson卫星沿几个典型方向(赤纬δ和赤经α)的投影面积，给出了本发明计算结果、可靠的文献结果和两者的相对误差。赤纬δ＝90°时，投影方向与机体坐标系的z轴重合，故投影面积与赤经无关。试验粒子总数为2×10⁷，对应的统计误差为满足工程需求。从表中可以看到，本发明计算结果和文献报道的数据一致，且相对误差的量级为与统计误差量级相同。因此，本发明提出的montecarlo技术适用于复杂工程应用外形，具备准确计算复杂外形航天器沿任意方向投影面积的能力，且相对误差满足统计分布规律。

图6是samson卫星沿不同方向的投影面积，图6(a)为投影面积三维分布图，图6(b)为对应的二维等值线图，赤纬的变化范围是δ∈[-90°,90°]，赤经的变化范围为α∈[0°,360°]。赤纬δ＝-90°时，投影方向为z轴反方向；赤纬δ＝0°时，投影方向平行于xoy平面；赤纬δ＝90°时，投影方向为z轴正方向。显然，samson卫星沿不同方向的投影面积变化范围跨度较大，从最小值apmin＝0.03m²变化到最大值apmax＝0.194m²。此外，当太阳帆展开时，赤纬，即投影方向与xoy平面的夹角，对投影面积的影响更为明显。相对而言，赤经的影响更小一些，这是samson卫星的几何外形所决定的。在赤纬δ＝0°附近，投影面积最小，此时赤经对投影面积的影响最大；赤纬δ＝±90°时，投影面积最大，此时赤经对投影面积没有影响。图6揭示了卫星姿态与其投影面积的关系，是卫星姿态控制的必要输入条件。

表1samson卫星典型姿态下的投影面积

最后，作为投影面积montecarlo模拟方法在航天工程的应用，考虑欧洲空间局的重力场和稳态海洋环流探测卫星(gravityfieldandsteadystateoceancirculationexplorer,goce)，如图7(a)、7(b)所示。

图8是goce卫星沿不同方向的投影面积，图8(a)为投影面积三维分布图，图8(b)为对应的二维等值线图，赤纬的变化范围是δ∈[-90°,90°]，赤经的变化范围为α∈[0°,360°]。与samon卫星类似，在赤纬δ＝0°附近，goce投影面积最小，此时赤经对投影面积的影响最大；赤纬δ＝±90°时，投影面积最大，此时赤经对投影面积没有影响。然而，goce卫星细长体的外形，外加太阳翼较大的面积，导致其沿不同方向的投影面积跨度比samson卫星更大，从最小值apmin＝0.817m²变化到最大值apmax＝10.273m²。随着其姿态的改变，goce卫星投影面积如此之大的变化范围导致其所受大气阻力和太阳光压变化巨大，同时也会引起太阳翼发电效率的剧烈变化，对其稳定性、姿态和轨道控制提出了巨大的挑战。图8揭示了goce卫星姿态与其投影面积的关系，是太阳翼发电量计算与监测、大气阻力和太阳光压预测的输入参数，也是卫星姿态和轨道控制的必要条件。

综合以上单面圆形平板模型、samson卫星和goce卫星投影面积的计算可以得到以下结论：本发明航天器投影面积计算方法具备准确计算任意复杂航天器外形沿任意方向投影面积的能力，鲁棒性和通用性好，能够给出卫星姿态与其投影面积的关系，是复杂外形航天器投影面积计算的理想方法。

本发明未公开技术属本领域技术人员公知常识。