微小卫星飞行姿态动环境规划算法创新机制

　0 引 言
　　随着微小卫星在侦察、导航、数据中继及早期预警等领域得到越来越广泛的应用，微小卫星的安全性能也越来越得到重视。空间探测设备，反卫星武器的快速发展越来越多地限制了微小卫星的使用[1?3]。为提高在轨微小卫星的使用效能及生存能力，对其进行低可观测性设计是一个有效的方案。目前，对中、低轨卫星的探测主要由雷达完成，因此低可观测性设计中低雷达探测特性是关键。除了对微小卫星的外形进行低可观测外形设计外，对其进行适当的飞行姿态规划以提高其低可观测特性也很有必要[4?5]。
　　对微小卫星的飞行姿态规划可从静环境和动环境两个方面进行规划。静环境主要是指环境中的威胁分布和各威胁设备的威胁性均已知，可以精确地进行规划，有效地提高卫星的低可观测特性；动环境规划则是规划前并不知道威胁分布，也不知道有哪些威胁设备及威胁设备的威胁性，只有当卫星处于当前位置时才获得该位置所对应的威胁设备的分布情况，所以需要对卫星在飞行过程中进行实时调整来降低卫星的低可观测性能。本文主要对微小卫星动环境规划进行研究。   []
　　微小卫星姿态规划与无人机等的航迹规划的基本思想都是在一定的约束条件下寻找一种飞行方式，使无人机等所受到威胁最低，所以对微小卫星的姿态规划可以参考航迹规划方面的算法。目前国内外对路径规划算法的研究很多，常见的有Dijkstra算法、A*算法及相关的改进A*算法[6]、概率地图方法（Probabilistic Roadmap Method，PRM）[7]、Voronoi图法[8]、快速扩展随机树算法（Rapidly Random Exploring Trees，RRT）[9]等以及智能优化算法如粒子群算法[10]，遗传算法[11]，混沌算法[12]等，前几种方法都属于图搜索算法，算法的实现需要知道威胁分布情况，对于动环境规划并不适用，所以本文采用智能优化算法对微小卫星低可观测飞行姿态进行动环境规划。
　　1 混沌优化算法
　　混沌优化算法（Chaos Optimization Algorithm）的基本思想是把混沌变量从混沌空间映射到解空间，然后利用混沌变量具有遍历性、随机性和规律性的特点进行搜索。混沌优化算法具有对初值敏感、易跳出局部极小、搜索速度快、计算精度高、全局渐近收敛的特点。
　　混沌算法首先应使优化变量在空间中处于混沌状态，一般选用式（1）所示的Logisic映射来产生混沌变量，其中[μ]是控制参量。设[0≤xn≤1，n=0，1，2，…，][μ=4]时，Logisic映射完全处于混沌状态[12]。混沌优化算法需要在混沌区间任意设定[i]个初值（不能为方程（1）的不动点0.25，0.5，0.75）作为[i]个初始混沌变量。
　　[xn+1=μxn（1-xn）， 0≤xn≤1] （1）
　　设一类连续对象的优化问题为求函数最小值，如式（2）所示：
　　[minf（xi）， i=1，2，…，n ； a≤xi≤b] （2）
　　采用混沌优化算法的基本步骤如下：
　　（1） 算法初始化：给式（1）中的[xn]赋予[i]个大小在[0，1]之间的初值，可以得到[i]个变量[xi.n+1，]将这[i]个变量作为初始混沌变量；
　　（2） 首先初始化相应的性能指标，随机设定一个最优解[xopt，][fopt]作为初始最优解，因为优化的是最小值，所以先初始化一个较大的[fopt，]保证[fopt]在接下来的迭代搜索中能取到当前优化问题的最优解；
　　（3） 用混沌变量进行迭代搜索；
　　Step1：因为优化问题的变量区间为[[a，b]，]所以需要将[i]个混沌变量[xi.n+1，]按式（3）分别转化为[i]个混沌变量[x′i.n+1，][x′i.n+1]在[[a，b]]之间：
　　[x′i.n+1=a+（b-a）xi.n+1] （3）
　　Step2：用混沌变量[x′i.n+1]进行迭代搜索，[k]是迭代次数，[fi（k）]表示变量为[x′i.n+1]时的函数值。若[fi（k）≤fopt，] 则[fopt=fi（k），][xopt=x′i.n+1，]否则不执行任何操作；
　　Step3：[k=k+1；]
　　Step4：将[x′i.n+1]通过式（3）的逆运算转换回混沌变量区间[0，1]，即将[x′i.n+1]转换回[xi.n+1；]
　　Step5：将[xi.n+1]代入式（1）中的[xn，]继续进行混沌映射，得到新的混沌变量；
　　Step6：重复步骤（3），当[fopt]保持不变时，迭代搜索完成，输出[xopt，][fopt。]
　　图1为初始值分别为0.501和0.502两个点的混沌演化轨道，初始距离仅为0.001，采用Logistic映射迭代后，两个点逐渐分离。可以看出，前6次迭代两个点的距离还比较近，在图中表现为两个点基本重合，当迭代次数超过7次后，两个点迅速分离，分别按照各自的混沌轨道运行。由此例可以看出混沌运动的初值敏感性和遍历性。
　图1 0.501和0.502两点的混沌演化轨道
　　2 粒子群优化算法
　　粒子群优化算法（Particle Swarm Optimization，PSO）是基于群体智能的一种进化计算方法。PSO算法中每个粒子就是解空间中的一个解，它根据自己的经历和整个粒子群的经历来调整自己。每个粒子在飞行过程的最好位置就是该粒子本身所找到的最优解。整个粒子群在飞行过程中经历的最好位置就是整个种群目前所找到的最优解。粒子本身的最优解叫做个体最优解[pbest，]整个粒子群体找到的最优解叫做全局最优解[gbest。]每个粒子根据上述两个最优解结合更新公式不断更新自己，从而产生新一代粒子群体，粒子的“好坏”程度由适应度函数来评价。与一般的进化算法相比，PSO概念简单、容易实现并且需要调整的参数少，目前广泛应用于各种优化领域。
　　在PSO算法中，每个粒子可以看作是解空间中的一个点，假设粒子种群规模为[N，]则第[i（i=1，2，…，N）]个粒子的位置可表示为[Xi。]粒子的位置即是适应度函数的变量，根据粒子位置及适应度函数可以计算出粒子的适应度，根据适应度判断粒子的“好坏”程度。粒子所经历过的“最好”位置记为[pbest（i）]，该粒子的更新速度用[Vi]来表示，粒子群体所经历的“最好”的位置的用[gbest]表示，第[i]个粒子的速度和位置更新公式为：
　　[Vi（k+1）=ω?Vi（k）+c1?r1?（pbest（i）-Xi（k））+c2?r2?（gbest-Xi（k）） ] （4）
　　[Xi（k+1）=Xi（k）+Vi（k+1） ] （5）
　　式中：[c1，][c2]为常数，称为学习因 子；[r1，][r2]是[0，1]上的随机数；[ω]称为惯性权重，同时粒子在更新自己的速度和位置的时候还受最大更新速度[Vmax]和最小更新速度[Vmin]的约束[13?15]，即[Vi∈[Vmax，Vmin]]。
　　粒子群算法的步骤如下：
　　（1） 首先在搜索空间里初始化粒子种群。假设粒子种群规模为[M，]种群中的粒子记为[Xi，]粒子维数为[n，]表示为[（ai1，ai2，…，ain），]即粒子在搜索空间中位置的坐标。同时初始化每个粒子的飞行速度[V，]也是一个[n]维向量[（Vi1，Vi2，…，Vin），]其中[i=1，2，…，M；]
　　（2） 根据每个粒子的位置和适应度函数，计算出每个粒子的适应度；
　　（3） 比较适应度的大小，将每个粒子的当前适应度与该粒子的个体最优解相比较，若当前适应度优于个体最优解，则用当前位置取代个体最优解的位置；否则，个体最优解保持不变；
　　（4） 首先比较得出最优粒子，然后将该粒子的适应度与种群最优解比较，若当前适应度优于种群最优解，则用当前位置取代个体最优解的位置；否则，种群最优解保持不变；
　　（5） 根据更新公式（4），（5）来更新粒子的速度和位置；
　　（6）重复步骤（2）～（5），直到适应度函数达到最优或者满足迭代条件。
　　3 动环境规划算法设计   []
　　3.1 雷达分布及威胁水平计算模型
　　由于雷达是低可观测微小卫星最主要的威胁，所以在威胁水平评估建模中，主要考虑雷达作为威胁设备。当卫星处于某一位置时，并不是地球上所有的雷达都对卫星具有威胁性，而只是在某一角度内的雷达才对卫星具有威胁性，如图2所示。为了计算出卫星在某一位置时的威胁性水平，首先应分辨出哪些雷达对卫星具有威胁性[16]。
　　
　　图2 威胁区、威胁设备及威胁方向定义
　　建立如图2所示的直角坐标系。取微小卫星位置[S]为坐标系原点，[x]轴指向为星体飞行方向，[z]轴指向为地球质心方向。图2中，[Oe]为地球质心，▲表示威胁雷达或激光设备，△表示无威胁雷达或激光设备。
　　自卫星位置[S]向威胁分布球面引切线，由切点可确定一个平面[AOcB]垂直于[z]轴，平面[AOcB]球面划分曲面[ACB、]曲面[ADB]两个部分。由雷达的工作原理可知仅分布在曲面[ACB]上的雷达能够探测到微小卫星，因此，可将该曲面定义为威胁区，威胁区内的设备定义为威胁设备，威胁设备相对于微小卫星的方向称为威胁方向。
　　威胁区俯仰角[θz]为切线[SA]与[z]轴的夹角，[（θT，φT）]为威胁方向的俯仰角及方位角，[（θs，φs，γs）]为微小卫星低可观测外形的姿态，[θs，][φs，][γs]分别表示微小卫星相应的俯仰角、方位角及横滚角[4]。
　　威胁区对应俯仰角[θz]可根据卫星飞行高度[h]及地球半径[Re]按公式（6）求得：
　　[θz=arctanRe（Re+h）2-R2e] （6）
　　3.2 雷达威胁性定义
　　根据雷达方程及二元假设检验理论[4?5]，可得到雷达检测概率[p]随[σR4]单调递增。因此雷达的威胁性可由雷达监测概率来表示，本文将雷达威胁性[T]定义为[T=kσR4] （7）
　　式中：[R]为雷达作用距离；[σ]为微小卫星的雷达散射截面积（Radar Cross?Section，RCS）；[k]为雷达威胁等级加权系数。当一个区域内有多部雷达时，只要被一部雷达探测到，则微小卫星就被发现，所以只要保证威胁性最大的雷达探测不到微小卫星即可。将威胁性最大的雷达对微小卫星的威胁性定义为卫星在该区域面临的威胁水平[T，]即 ：
　　[T=maxT1（t），T2（t），…，Tn（t）] （8）
　　式中：[n]为威胁雷达数；[Ti]（[i=]1，2，[…，n]）为第[i]部雷达的威胁性。
　3.3 规划空间压缩
　　（1） 对卫星产生威胁的雷达是威胁区内的雷达，所以微小卫星俯仰角[θ]的规划只需要在威胁区俯仰角[θz]内即可，即：
　　[θ∈[0，θz]] （9）
　　（2） 根据姿控最大速度，前一时刻卫星的飞行姿态为[（θ0，φ0，γ0）]，则[t]时刻卫星外形姿态[（θ，φ，ν）]应满足：
　　[θ∈[θ0-ωθt，θ0+ωθt]] （10）
　　[φ∈[φ0-ωφt，φ0+ωφt]] （11）
　　[γ∈[γ0-ωγτ，γ0+ωγτ]] （12）
　　式中：[ωθ，ωφ，ωγ]分别为卫星俯仰角[θ，]方位角[φ，]横滚角[γ]单位时间内的最大调整角度。
　　（3） 目前，卫星低可观测外形多为锥形轴对称设计，所以对横滚角的调整对规划性能的影响很小，因此，可将4维规划空间[{θ，φ，γ，t}]简化为3维[{θ，φ，t}]，既可降低规划空间规模和规划计算复杂度，又可以降低规划算法设计难度[16]。
　　3.4 规划代价
　　卫星规划性能主要考虑两个方面：威胁水平和姿控能耗。提高微小卫星的低可观测性能首先要保证卫星拥有较低的威胁水平，其次，由于微小卫星的工作特点，对微小卫星的姿态的控制需要在有限的姿控能耗下完成。
　　假设威胁雷达的威胁等级相同，则可取威胁等级系数[k=1。]根据式（8）则可以得到，微小卫星在[t]时刻面临的威胁水平[T（t）]为：
　　[T（t）=max（σ1R41，σ2R42，…，σnR4n）] （13）
　　在[t]时刻，定义卫星的姿控能耗代价[C（t）]为：
　　[C（t）=aθ?Δθ，aφ?Δφ，aγ?Δγ] （14）
　　式中：[aθ，aφ，aγ]分别为[θ，φ，γ]的能耗加权系数；[Δθ，Δφ，][Δγ]分别为前一时刻到[t]时刻[θ，φ，γ]的变化角度，且满足：
　　[ΔθΔt≤ωθ] （15）
　　[ΔφΔt≤ωφ] （16）
　　[ΔγΔt≤ωγ] （17）
　　卫星的规划代价则为：
　　[f（t）=T（t）+C（t）] （18）
　　3.5 混沌粒子群（chaos PSO）规划算法
　　上述分析可知，[t]时刻雷达对微小卫星的威胁性可由[T]及[C]决定。对微小卫星飞行姿态规划就是寻找一组最优飞行姿态[（θt，φt，γt），]使式（18）所确定的[f（t）]的值最小 ，因此，可以采取智能优化算法进行规划。在粒子群算法的基础上引入混沌运动，增加了粒子在规划空间的遍历性，从而提高了规划算法的精确性。
　　根据粒子群优化算法特征[13?15]，惯性权值[ω]可定义为：
　　 （19）
　　式中：分别为惯性权值的初始值和终止值。粒子的适应度函数即为卫星的规划代价，粒子的位置即为卫星飞行姿态。
　　结合混沌优化算法和粒子群优化算法的基于局部近似最优解的混沌粒子群优化规划算法的步骤如下[17?18]：
　　（1） 以微小卫星之前的飞行姿态为基础，根据式（20）～（25） 初始化规模为[N]的初始粒子种群：
　　[pθ（i）=θinitial+rθ?ωθ?tc] （20）
　　[pφ（i）=φinitial+rφ?ωφ?tc] （21）
　　[pγ（i）=γinitial+rγ?ωγ?tc] （22）
　　[vθ，φ，γ（i）=0] （23）
　　[fopt（i）=∞] （24）
　　[gopt=∞] （25）
　　式中：[θinitial，][φinitial，][γinitial]分别为动环境规划中的初始规划姿态（俯仰角、方位角、横滚角），即微小卫星前一时刻的飞行姿态，[rθ，][rφ，][rγ]为[-1，1]内的随机数；[tc]为规划时长，[pθ，][pφ，][pγ]为初始化的粒子飞行姿态，[i]表示第[i]个粒子，[vθ，φ，γ（i）]表示第[i]个粒子的初始更新速度；[fopt]为个体最优解；[gopt]为种群最优解。
　　（2） 计算出当前位置内的威胁雷达数目，根据威胁雷达分布计算出卫星规划代价。
　　（3） 将当前的规划代价与个体最优解相比较，若小于个体最优解，则用当前姿态代替粒子最优姿态，保存粒子个体最优姿态，同样的，将当前种群中的最优解与先前保存的种群最优解相比较，若小于先前种群最优解，则用当前种群最优个体的信息取代先前的种群最优个体的信息，并保存种群最优个体。
　　（4） 更新粒子信息
　　速度更新公式：
　　[vθ（k+1）=ω?vθ（k）+c1?r1?（θfoptn（i）-θ（i））+c2?r2?（θgopt-θ（i））] （26）
　　[vφ（k+1）=ω?vφ（k）+c1?r1?（φfoptn（i）-φ（i））+c2?r2?（φgopt-φ（i））] （27）
　　[vγ（k+1）=ω?vγ（k）+c1?r1?（γfoptn（i）-γ（i））+c2?r2?（γgopt-γ（i））] （28）
　　姿态更新公式：
　　[θ（k+1）=θ（k）+vθ（k+1）] （29）
　　[φ（k+1）=φ（k）+vφ（k+1）] （30）
　　[γ（k+1）=γ（k）+vγ（k+1）] （31）
　　式中：[c1，][c2]是常数，称为学习因子；[r1，][r2]是[0，1]上的随机数，[ω]也为一个常数，称为惯性权值；粒子移动速度[v∈（vmin，vmax），][vmin，][vmax]是常数，用来设定粒子移动速度。粒子移动速度与惯性权值都是用来维护全局最优解与局部最优解的平衡，本文引入动态权值[ω，]即减小了粒子群算法对粒子移动速度的依赖，所以本文中粒子移动速度的设定对粒子群算法性能影响较小。
　若更新后的姿态不满足式（10）～（12）的范围，若 [θ<θ0-ωθt，]则[θ=θ0-ωθt；]若[θ>θ0+ωθt，]则[θ=θ0+ωθt。]同理判断出[φ，][γ]的取值。
　　（5） 如果粒子有重叠现象，即粒子聚集在某个极值附近的情况，则保留其中一个粒子不变，其他粒子赋予混沌运动，首先将粒子逆运算到混沌区间即[0，1]，若粒子区间为[[a，b]，]根据逆运算公式[x=（x-a）（b-a）]转换，然后再采用Logistic映射进行迭代，再根据式（3）变换到优化变量空间中，其中需要判断当更新后的姿态不满足式（10）～（12）的范围时，若[θ<θ0-ωθt，]则[θ=θ0-ωθt；]若[θ>θ0+ωθt，]则[θ=θ0+ωθt。]同理判断出[φ，][γ]的取值。
　　（6） 重复算法步骤（2）～（5）直到满足迭代数或者最优解稳定，输出最优粒子的信息。
　　4 算法仿真结果及分析   []
　　本文的仿真条件如下：地球半径取6 371 km，卫星的飞行高度为500 km，俯仰角。
　　首先对卫星在处于某一位置时，分别采用粒子群算法和混沌粒子群算法进行规划仿真。当两种算法的迭代次数均为1～50次时，规划得到的最小威胁代价水平如图3所示。
　　
　　图3 迭代次数为1～50次时chaos PSO和PSO的规划结果
　　当算法迭代次数均为50次时，将两种算法均重复运行50次，规划得到的最小威胁水平如图4所示。
　　从仿真结果可以看出，随着迭代次数的增加，混沌粒子群算 法规划处的结果趋于稳定，而粒子群算法则有一定的波动性，因为粒子群算法本身容易陷入局部最优解，混沌运动的遍历性和随机性可以使其跳出局部最优解，增强了粒子群算法的搜索性能；同时迭代次数均为50次时，混沌粒子群算法也比粒子群算法更稳定，所以可以看出混沌粒子群算法性能更优。
　　
　　图4 迭代次数均为50次时chaos PSO和PSO
　　重复运行50次的规划结果
　　接下来对卫星绕地球运行一周进行规划，每个位置算法规划迭代次数均为50次，将两种算法进行比较的同时与无规划（[θ，][φ]均为[0°]）情况进行比较。卫星绕球心每飞行[1°]进行飞行姿态规划，即[tc=1°]。
　　图5所示为雷达在300 MHz工作频率时对卫星进行飞行姿态规划后卫星在飞行过程中的俯仰角和方位角的变化，实线和点线分别表示采用chaos PSO规划后的卫星俯仰角[θ]，方位角[φ]的变化，点划线和虚线分别表示采用PSO规划后的卫星俯仰角[θ，]方位角[φ]的变化。由于设定的俯仰角的能耗系数比方位角的大，所以在图中俯仰角的变化比较小，主要通过调整方位角来减小卫星威胁代价。
　　
　　图5 300 MHz频率采用chaos PSO和PSO规划时
　　卫星飞行俯仰角[θ，]方位角[φ]变化
　　图6所示是进行飞行姿态规划后的卫星在飞行过程中的威胁代价图。实线（chaos PSO）采用chaos PSO算法规划后的威胁代价水平分布。点线（PSO）表示的是采用PSO算法进行规划后的威胁代价水平分布，虚线（noplan）表示的则是无规划时的卫星威胁代价水平分布。
　　表1 300 MHz频率卫星无规划时，PSO规划时
　　及chaos PSO规划时的仿真数据
　　[＼&威胁代价 /dBsm＼&姿控能耗代价＼&规划代价＼&无规划＼&516.600 0＼&0＼&516.600 0＼&PSO＼&-466.263 1＼&24.040 0＼&-442.223 1＼&chaos PSO＼&-658.921 3＼&18.640 0＼&-640.281 3＼&]
　　由表1和图5，图6可以看出不管是PSO算法还是chaos PSO算法，都可以在很大程度上降低卫星的威胁代价，图中在某些地方规划后的威胁代价大于未规划时的威胁代价，是因为图中的威胁水平曲线是微小卫星在当前位置面临的威胁水平，而本文的规划代价为初始位置到当前位置的总体威胁水平，所以会出现在某些位置的代价值偏大的情况，但是总体上的代价值会远远小于未规划时的代价。由于规划是实时规划，因为最大调整角度的约束，所以卫星某一时刻某个位置的规划结果很依赖于前一时刻的卫星飞行姿态。如图中所示，有些地方PSO的规划结果会优于chaos PSO就是因为前一时刻的飞行姿态不同。但是从整体规划来看，chaos PSO的规划结果还是优于PSO，因为chaos PSO对最优解的搜索能力比PSO更强更稳定。
　图6 300 MHz频率采用chaos PSO和PSO规划时
　　微小卫星的威胁代价水平分布
　　5 结 论
　　本文以粒子群算法和混沌算法为基础，对微小卫星进行了动环境的飞行姿态规划的算法设计。采用粒子群算法和混沌粒子群算法分别对微小卫星在某一位置和飞行一个周期两种情况进行规划，规划结果显示混沌粒子群算法比粒子群算法的规划性能更好，同时通过与无规划时微小卫星所面临的威胁水平相比较，可以看出两种算法均能有效提高微小卫星的低可观测性能。
　　微小卫星飞行姿态规划是个复杂的多学科课题，对卫星进行姿控规划的同时也对能源系统设计提出新的要求， 同时进行姿态规划可能会对卫星功能及寿命产生一定影响。这些问题需要通过多学科的交流合作进行进一步研究。
　　参考文献
　　[1] 黄汉文.卫星隐身概念研究[J].航天电子对抗，2010，26（6）：22?24.
　　. Avitio Week&Space Technology， 2001， 8： 82?85.
　　[3] 张建民，孙健.国外隐身技术的应用与发展分析[J].舰船电子工程，2012，32（4）：18?21.
　　[4] 苏抗，周建江.有限姿控能力的低可观测微小卫星姿态实时规划[J].航空学报，2010，31（9）：1841?1846.
　　[5] 苏抗.微小卫星低可观测关键技术研究[D].南京：南京航空航天大学，2011.
　　[6] 史辉，曹闻，朱述龙，等.A*算法的改进及其在路径规划中的应用[J].测绘与空间地理信息，2009，32（6）：208?211.
　　[7] 刘贤敏，王琪，许书诚.改进型概率地图航迹规划方法[J].火力与指挥控制，2012，37（4）：121?124.
　　[8] 吴海彬，林宜.基于改进Voronoi图的移动机器人在线路径规划[J].中国工程机械学报，2007，5（1）：117?121.
　　[9] 崔挺，李俨，张明庄.基于改进RRT算法的无人机航迹规划[J].电子设计工程，2013，21（12）：50?53.
　　[10] 王新增，慈林林，李俊山，等.基于改进粒子群优化算法的无人机实时航迹规划[J].微电子学与计算机，2011，28（4）：87?90.   []
　　[11] 鲁艺，吕跃，罗燕，等.基于改进遗传算法的UAV航迹规划[J].电光与控制，2012，19（1）：29?33.
　　[12] 李兵，蒋慰孙.混沌优化方法及其应用[J].控制理论与应用，1997，14（4）：613?615.
　　[13] MASEHIAN E， SEDIGHIZADEH D. A multi?objective PSO?based algorithm for robot path planning [C]// 2010 IEEE International Conference on Industrial Technology （ICIT）. [S.l.]： IEEE， 2010： 465?470.
　　[14] 吴晓军，李峰，马悦，等.均匀搜索粒子群算法的收敛性分析[J].电子学报，2012，40（6）：1115?1120.
　　. Applied Soft Computing， 2011， 11（4）： 3658?3670.
　　[16] 苏抗，周建江.微小卫星低可观测外形飞行姿态规划[J].航空学报，2011，32（4）：720?728.