基于等碰撞概率面法的航天器安全接近控制方法与流程

本发明涉及航天器近距离安全接近控制技术领域，具体涉及一种基于等碰撞概率面的航天器安全接近控制方法。

背景技术：

近年来随着在轨航天器失效案例日益增多，在轨服务受到了更多的关注。航天器在轨服务对提高航天器操作能力，增强航天器适应性和可伸缩性，延长航天器寿命，节省航天器研制费用等有积极作用。

近距离操作是在轨服务的重要关键技术，涉及到近距离操作的空间任务几乎都需要进行航天器相对控制，当服务航天器靠近目标时，需要保证航天器不与目标发生碰撞，因此航天器安全接近控制是近距离操作技术的关键技术之一。

目前，在实现航天器安全接近的方法中，人工势函数制导方法(artificialpotentialfunction,apf)由于具有在复杂环境下形式简单、计算量小等特点，但近年来在空间近距离操作中得到了广泛关注。

人工势函数制导方法最初源自机器人路径规划问题，ismaellopez和colinr.mcinnes首次将人工势函数制导方法应用到航天器自主避障的路径控制上来，在其文献中根据航天器相对运动c-w方程，将人工势函数制导方法分别应用于多个障碍物的交会和终端带约束的交会两个场景中。hyung-chullim等将人工势函数与滑模控制相结合，设计了针对卫星编队飞行的控制器。enderst.john-olcayto和colinr.mcinnes等人研究了在国际空间站的交会和近距离作业中，采用势函数制导的安全自主近距离机动，实现了v/r-bar方向的安全自主交会对接。shawnb.mccamish和marcelloromano等人针对圆轨道上多个航天器近距离交会问题提出了一种基于linearquadraticregulator(lqr)和人工势函数方法相结合的制导算法，并利用spheres进行地面试验，验证了所设计算法的可行性。lenelpalacios和matteoceriotti等利用lqr和人工势函数方法设计了一种适用于圆轨道和椭圆轨道的制导追踪算法，并利用并利用safe(spaceborneautonomousformationflyingexperiment)进行了验证。josued.munoz和normang.fitz-coy提出了一种自适应的人工势函数制导方法，针对静态障碍物，研究了航天器近距离操作中的轨道和姿态运动快速规划，比传统人工势函数制导方法更加节省燃料。高鹏和罗建军在josued.munoz和normang.fitz-coy工作的基础上，针对航天器近距离操作中的动态障碍物规避问题，研究了一种燃料较省、精度较高的规避动态障碍物的自适应人工势函数制导方法。张大伟针对非合作航天器，提出了将人工势函数与模糊控制相结合的制导方法，能够实现对动态障碍物的规避。针对apf局部极小点问题，可寻找更好的势场函数，使势场中出现极小值的概率减小；或将人工势场法和其他方法相融合，如采用沿墙跟踪方法及极限环法。鲁新军等提出了虚拟水流法，但效率不高。基于遗传算法的改进apf法实现最优路径搜索，收敛速度有待提高。为了解决实际应用中控制对象在全动态环境中的路径规划问题，殷路等提出了全动态环境中改进的势能函数构造方法，使研究对象时刻趋向与目标相同、与障碍物相悖的运动趋势。

但是目前人工势函数中的势函数未考虑近距离操作中的导航不确定性和控制不确定性，而目前国内外学者普遍采用碰撞概率函数作为不确定影响下对航天器可能发生的碰撞进行预警。1997年，chan针对空间中的两个物体，在两者位置误差互不关注的前提下，对误差协方差矩阵的可叠加性做出了验证，为碰撞概率研究奠定了基础。早期研究的航天器和空间碎片的碰撞问题主要针对典型的线性相对运动碰撞，其计算方法比较成熟。线性相对运动碰撞概率计算的基本假设为：

1)相遇时间非常短(通常为几秒量级)，相对运动可以简化为线性模型；

2)相对速度矢量在相遇期间保持不变，积分体近似为无限长圆柱体；

3)两空间目标位置误差pdf在相遇期间保持不变。

目前在大部分情况下，线性相对运动碰撞概率计算都适用。但是在近距离操作中，接近时相对速度较小，目标之间会长时间接近，接近期间相对速度的大小和方向都会发生变化，目标之间是非线性相对运动。此时线性相对运动假设的碰撞概率的描述方法和计算方法不再适用。对于非线性情况，chan分析了线性相对运动假设成立的条件，并讨论了非线性相对碰撞概率的计算方法。patera通过变化将相对位置误差协方差矩阵转换成三维空间中对称形式，此时平行于相对速度的一维就可以和其他两维互相解耦，然后得到某一小段事件内的碰撞概率分量，在整个事件间隔内对其积分就可以得到非线性碰撞概率值。alfano提出了两种计算非线性相对运动情况下碰撞概率计算方法。第一种方法对线性碰撞概率计算方法进行了扩展，将碰撞管道分割为一系列相邻的圆柱体，在圆柱体内运动是近似线性的。计算每段对应的碰撞概率，然后求和。第二种方法将碰撞管道在马氏空间内进行体积有限元划分，表示为在马氏空间中建立体积有限元网格，计算每个有限元网格内的概率并求和。mckinley通过建立积分坐标系，在此坐标系中积分得到碰撞概率。针对卫星与空间绳系卫星之间的相对运动，chobotovpatera采用一重曲线积分的方法对碰撞概率进行了计算，并对相对运动弱非线性运动时的碰撞概率进行了研究。但是，碰撞概率方法没有给出碰撞概率与碰撞参数之间的直接的和显式的联系，难以得到碰撞概率与影响因素之间的关系；而且碰撞概率法容易受到协方差不确定性的影响；而且其计算量相对于人工势函数法仍然很大。

技术实现要素：

针对现有技术中存在的缺陷，针对航天器近距离操作中的安全接近控制技术中的问题，本发明提供了基于等碰撞概率面的航天器安全接近控制方法。

为实现本发明之目的，采用以下技术方案予以实现：

基于等碰撞概率面的航天器安全接近控制方法，包括以下步骤：

(1)建立目标航天器的轨道坐标系；

(2)在目标航天器的轨道坐标系中，构建跟踪航天器的非线性相对运动方程；

(3)采用无迹变化方法来推导跟踪航天器的非线性相对运动方程的不确定性传播；

(4)计算目标航天器和跟踪航天器之间的碰撞概率密度；

(5)计算目标航天器和跟踪航天器两个卫星之间的碰撞概率梯度；

(6)确定最终作用在跟踪航天器上的控制力。

步骤(1)坐标系的建立方法如下：

以o-xyz表示j2000地球惯性坐标系，原点为地心，x轴指向历元j2000春分点，基本面为地球平赤道面，z轴指向地球北极，根据右手螺旋定律由x轴与z轴得出y轴的方向。

以o-xyz表示目标航天器的轨道坐标系，以目标航天器的质心为坐标原点，x轴由地心指向跟踪航天器的质心，z轴指向跟踪航天器轨道面法线方向，y轴由右手法则，垂直于x轴和z轴所组成的平面。

rtarget表示目标航天器在j2000地球惯性坐标系下的位置矢量，且rtarget＝[x1y1,z1]^t。

rtrack表示跟踪航天器在j2000地球惯性坐标系下的位置矢量，且rtrack＝[x2,y2,z2]^t。

rlvlh表示跟踪航天器在目标航天器的轨道坐标系下的位置矢量，且rlvlh＝[x,y,z]^t。

vlvlh表示跟踪航天器在目标航天器的轨道坐标系下的速度矢量，且vlvlh＝[vx,vy,vz]^t。

步骤(2)中，跟踪航天器的非线性相对运动方程，如下所示：

其中：x,y,z表示目标航天器在三个方向的位置参数；表示目标航天器在三个方向的速度参数；表示目标航天器在三个方向的加速度参数；r1表示目标航天器的质心相距地心的距离；r2表示跟踪航天器的质心相距目标航天器质心的距离；μ表示地球引力常数；u＝[ux,uy,uz]^t表示控制力；ω,分别表示目标航天器的转动速度和转动角速度，且如下表示：

其中：a是目标航天器轨道半长轴；e是目标航天器偏心率；f是目标航天器真近点角度；

假设状态向量x＝[rlvlh^t,vlvlh^t]^t，则方程(1)可改写成如下形式

xk+1＝f(xk,uk)(4)

其中xk+1表示第k+1步相对状态矢量；xk表示第k步相对状态矢量；uk表示第k步输入控制状态矢量。

步骤(3)中采用无迹变化方法推导方程(4)的不确定性传播，方法如下：

当l表示状态向量x的长度时，n＝2l+1表示需要的采样数。无迹变化方法是：取样一系列均值为m(t0)，方差为p(t0)的初始采样点r(tk)(k＝1,2,...,n)，然后代入方程(4)计算下一时刻的输出期望和方差。

其中：ωk表示采样点的权重因子。

步骤(4)中，目标航天器和跟踪航天器之间的碰撞概率密度如式(6)计算得到：

其中：表示目标航天器在三个方向的位置不确定性的方差；表示相对位置不确定性的协方差矩阵；将目标航天器简化成包络椭圆，v表示目标航天器的几何尺寸；r3表示在包络椭圆里的相对位置；r2表示跟踪航天器的相对位置。

步骤(5)的实现方法如下：

在时刻t1，跟踪航天器在目标航天器的轨道坐标系下的位置矢量可表示为其中相对位置矢量和相对速度矢量分别为则由式(6)可得目标航天器和跟踪航天器两个卫星之间的碰撞概率计算公式为

其中：表示t1时刻目标航天器在三个方向的位置不确定性的方差；表示t1时刻相对位置不确定性的协方差矩阵；

时刻t1的碰撞概率梯度即此时目标航天器和跟踪航天器之间的碰撞概率对跟踪航天器对位置矢量求导，如式(8)所示：

其中：表示在时刻t1时跟踪航天器的相对位置矢量，表示沿碰撞概率梯度在跟踪航天器上施加避撞作用力后跟踪航天器在轨道坐标系中的运动矢量；

由式(8)可得到碰撞概率梯度为

按照式(9)进行碰撞概率的梯度计算时，计算量比较大。为了减少计算量，本发明步骤(5)也可采用以下方法实现：

采用辅助函数即式(10)近似计算碰撞概率梯度，

其中：表示t1时刻目标航天器在三个方向的位置不确定性的方差；表示t1时刻相对位置不确定性的协方差矩阵；λ0是用来确定推力数量级的正数；d0是用来描述目标航天器影响区域的半径表示相对平行速度，且

其中：表示由追踪航天器指向目标航天器的单位向量；当时，表示跟踪航天器远离目标航天器，因此不需要在跟踪航天器上施加避撞作用力；反之，跟踪航天器飞向目标航天器，且同时在影响区域范围内时，需要在跟踪航天器上施加避撞作用力。

目标航天器的影响区域范围被定义为

d0＝d0(r0+ds)(12)

其中：d0表示一个正常数；r0表示目标航天器外包络半径；ds表示最小停止距离，且

其中：amax表示跟踪航天器上控制器的最大推力加速度；

由式(10)可得，施加在跟踪航天器上的避撞作用力uavoidance为

步骤(6)中，最终作用在跟踪航天器上的控制力是复合控制力，其由保持轨迹正常的跟踪力utrack和避撞作用力uavoidance组成，表示为

utotal＝uavoidance+utrack(15)

相对于现有技术，本发明有效地解决了目前安全接近控制方法中存在的计算量大、实时性不好以及未能考虑不确定性影响等问题。本发明产生了以下有益技术效果：

(1)考虑到两目标交会时的位置速度几何关系以及航天器与危险目标位置速度的不确定性；

(2)复杂度不高，能够获得解析表达式的制导方法；

(3)能够在理论上证明其稳定性，便于理论分析；

(4)计算效率高，实时性好，便于航天器自主安全接近。

附图说明

图1是本发明第一步建立的坐标系示意图。

具体实施方式

下面将结合本发明实施例图中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，做进一步详细说明，但本发明的实施方式不仅限于此。

基于等碰撞概率面的航天器安全接近控制方法，包括以下步骤：

(1)建立目标航天器的轨道坐标系；

参照图1，以o-xyz表示j2000地球惯性坐标系，原点为地心，x轴指向历元j2000春分点，基本面为地球平赤道面，z轴指向地球北极，根据右手螺旋定律由x轴与z轴得出y轴的方向；

以o-xyz表示目标航天器的轨道坐标系，以目标航天器(target)的质心为坐标原点，x轴由地心指向跟踪航天器(track)的质心，z轴指向跟踪航天器轨道面法线方向，y轴由右手法则，垂直于x轴和z轴所组成的平面。

j2000地球惯性坐标系o-xyz和轨道坐标系o-xyz如图1所示。同时如图1所示：

rtarget表示目标航天器在j2000地球惯性坐标系下的位置矢量，且rtarget＝[x1y1,z1]^t；

rtrack表示跟踪航天器在j2000地球惯性坐标系下的位置矢量，且rtrack＝[x2,y2,z2]^t；

rlvlh表示跟踪航天器在目标航天器的轨道坐标系下的位置矢量，且rlvlh＝[x,y,z]^t；

vlvlh表示跟踪航天器在目标航天器的轨道坐标系下的速度矢量，且vlvlh＝[vx,vy,vz]^t。

(2)在目标航天器的轨道坐标系中，构建跟踪航天器的非线性相对运动方程，如下所示：

其中：x,y,z表示目标航天器在三个方向的位置参数；表示目标航天器在三个方向的速度参数；表示目标航天器在三个方向的加速度参数。r1表示目标航天器的质心相距地心的距离；r2表示跟踪航天器的质心相距目标航天器质心的距离；μ表示地球引力常数；u＝[ux,uy,uz]^t表示控制力；ω,分别表示目标航天器的转动速度和转动角速度，且如下表示：

其中：a是目标航天器轨道半长轴；e是目标航天器偏心率；f是目标航天器真近点角度。

假设状态向量x＝[rlvlh^t,vlvlh^t]^t，则方程(1)可改写成如下形式

xk+1＝f(xk,uk)(4)

其中xk+1表示第k+1步相对状态矢量；xk表示第k步相对状态矢量；uk表示第k步输入控制状态矢量。

(3)采用无迹变化方法(unscentedtransformation,ut)来推导方程(4)的不确定性传播。

当l表示状态向量x的长度时，n＝2l+1表示需要的采样数。无迹变化方法是：取样一系列均值为m(t0)，方差为p(t0)的初始采样点r(tk)(k＝1,2,...,n)，然后代入方程(4)计算下一时刻的输出期望和方差。

其中：ωk表示采样点的权重因子。

采用无迹变化方法(unscentedtransformation,ut)来推导方程(4)的不确定性传播，具体计算步骤如下所示：

令l表示状态向量x的长度，n＝2l+1表示需要的采样数，给定初始采样点r(tk)(k＝1,2,...,n)，初始采样点的均值为m(t0)，方差为p(t0)；产生n个西格玛采样点和权重如下所示：

其中：α,β,κ,λ是尺度参数，α>0,β>0,κ>0,λ>0；分别是均值和方差的权重矢量；

得到最终时刻的期望值m(tf)和方差p(tf)

χ⁽ⁱ⁾(tf)＝f(χ⁽ⁱ⁾(t0),uk-1),i＝0,1,…,2l(18)

其中：χ⁽ⁱ⁾(t0)是矩阵χ(t0)的第i列；χ⁽ⁱ⁾(tf)表示tf时刻的相对状态矢量；χ⁽ⁱ⁾(t0)表示t0时刻的相对状态矢量；uk-1表示第k-1步输入控制状态矢量。

(4)计算目标航天器和跟踪航天器之间的碰撞概率密度；

通过静态轨道确定和不确定性推导理论，目标航天器和跟踪航天器之间的碰撞概率密度可以如下计算得到：

其中：表示目标航天器在三个方向的位置不确定性的方差，由式(18)至式(20)计算得到；表示相对位置不确定性的协方差矩阵；v表示目标航天器的几何尺寸，在此将目标航天器简化成包络椭圆；r3表示在包络椭圆里的相对位置；r2表示跟踪航天器的相对位置。

由式(6)计算出某一时刻目标航天器周边点的碰撞概率，并且将碰撞概率相同的点连接成曲面形成等碰撞概率面。该方法与等压面类似，是指在空间中碰撞概率相等的点所组成的曲面。

(5)计算目标航天器和跟踪航天器两个卫星之间的碰撞概率梯度。对于该步骤，本发明给出了两种实现方案，具体如下：

方案一：在时刻t1，跟踪航天器在目标航天器的轨道坐标系下的位置矢量可表示为其中相对位置矢量和相对速度矢量分别为则由式(6)可得目标航天器和跟踪航天器两个卫星之间的碰撞概率计算公式为

其中：表示t1时刻目标航天器在三个方向的位置不确定性的方差；表示t1时刻相对位置不确定性的协方差矩阵；

由式(7)可得时刻t1的碰撞概率的梯度，即此时目标航天器和跟踪航天器之间的碰撞概率对跟踪航天器对位置矢量求导，

其中：表示在时刻t1时跟踪航天器的相对位置矢量，表示沿碰撞概率梯度在跟踪航天器上施加避撞作用力后跟踪航天器在轨道坐标系中的运动矢量。

由式(8)可得到碰撞概率的梯度为

方案二：如果按照式(9)进行碰撞概率的梯度计算时，计算量比较大。为了减少计算量，采用一个辅助函数，当不确定足够大或者运行时间足够长时，该辅助函数能够近似计算碰撞概率的梯度，且设置一个实施避撞机动的条件，则该辅助函数可表示如下

其中：λ0是用来确定推力数量级的正数；d0是用来描述目标航天器影响区域的半径表示相对平行速度，且

其中：表示由追踪航天器指向目标航天器的单位向量。当时，表示跟踪航天器远离目标航天器，因此不需要在跟踪航天器上施加避撞作用力。反之，跟踪航天器飞向目标航天器，且同时在影响区域范围内时，需要在跟踪航天器上施加避撞作用力。

目标航天器的影响区域范围被定义为

d0＝d0(r0+ds)(12)

其中：d0表示一个正常数；r0表示目标航天器外包络半径；ds表示最小停止距离，且

其中：amax表示跟踪航天器上控制器的最大推力加速度。

由式(10)可得，施加在跟踪航天器上的避撞作用力uavoidance为

(6)确定最终作用在跟踪航天器上的控制力

最终作用在跟踪航天器上的控制力是复合控制力，其由保持轨迹正常的跟踪力utrack和避撞作用力uavoidance组成，可表示为

utotal＝uavoidance+utrack(15)

综上所述，虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明，任何本领域普通技术人员，在不脱离本发明的精神和范围内，当可作各种更动与润饰，因此本发明的保护范围当视权利要求书界定的范围为准。