基于光度序列同时估计空间失稳目标进动和自旋速率方法技术

本发明专利技术公开了一种基于光度序列同时估计空间失稳目标进动和自旋速率方法，包括：获取空间失稳目标的光度序列；将光度序列变分模态分解为多个本征模态函数；计算每个本征模态函数与光度序列的互信息；抽取最大的两个互信息值对应的本征模态函数频率作为空间失稳目标的自旋频率和进动频率，将空间失稳目标的自旋频率和进动频率转化为空间失稳目标的进动和自旋速率。本发明专利技术采用VMD对光度序列进行模态分解，基于互信息理论，对目标的进动和自旋频率进行有效的提取，可较高精度同时估计出目标进动和自旋速率，克服了频谱分析法、自相关和交叉残差法等估计方法对空间目标转动状态估计能力不足的问题。为空间目标转动状态的估计提供一种新的参考思路。

【技术实现步骤摘要】
基于光度序列同时估计空间失稳目标进动和自旋速率方法
本专利技术涉及的是一种基于光度序列同时估计空间失稳目标进动和自旋速率方法，属于计算领域。
技术介绍
随着人类开发和利用外层空间脚步的加快，近地空间运行着大量的人造物体。由于碎片碰撞、卫星突发故障、在轨失效等原因，产生了数以万计的废弃空间物体，给空间安全带来了极大威胁。低轨道废弃目标会受气动力矩、重力梯度力矩、地磁力矩等力的作用，几年至几十年后才会再入大气层烧毁；处于轨道空间拥挤的高轨废弃航天器，如不采取主动措施，其将会一直存在。以航天器燃料加注、在轨维修及装配、空间垃圾清除为主的在轨服务任务(On-OrbitServing，OOS)变得越来越重要，获取失稳非合作目标转动状态是对其进行服务的前提。失控或者失稳的空间目标通常处于高速旋转状态，旋转速率是衡量空间目标状态的重要参数，由于距离远无法通过雷达、激光探测等方法对目标进行状态探测，通过地基望远镜获取目标的光度数据并进行处理，对其转动状态进行估计，并判断其工作状态，是一种有效的态势感知手段。失稳旋转的目标，由于受到外界力矩的干扰以及质量分布不均匀等因素的影响，目标初始角速度矢量与惯量轴不重合，目标转动表现为自旋轴绕角动量轴做锥旋进动，自旋轴进动速率小于自旋速率。近年来，基于目标光度序列的频谱分析法、自相关法和交叉残差法等方法，对空间目标转动状态进行了反演，其中，通过对光度序列进行傅里叶变换，在长周期和短周期选取频谱中能量高的两个频率作为自旋和进动的频率，近似的估计出进动周期，自旋频率的估计频谱中出现了倍频，估计误差较大；交叉残差法相较于自相关法优势最为明显，但其受采样频率的影响较大。从反演的结果可以看出，上述方法只能进行进动或者自旋速率的独立反演，对于在某一时刻光度数列包含两种频率的失稳空间目标，由地基观测得到的光度序列，无法同时确定光度序列中包含的进动频率和自旋频率，原因如下：1)频谱分析法、交叉残差法和自相关法不能进行目标自旋和进动速率的同时反演，只能单独估计出自旋或者进动中的一个速率。2)经验模态分解(EMD)在分解过程中，频率混叠比较严重，且其模态数是自适应的，不可控制；频谱分析法倍频比较严重，対目标旋转速率的估计误差较大。3)交叉残差法和自相关法分析：对于包含周期性离散数据的信号，自相关计算公式为：式中，N为测量数据的总数，自相关系数L(Δn)表示为原数据x(n)与时延Δn后的数据乘积的均值。可以用自相关图来描述数据的周期性，自相关图峰值之间的间隔数值表示周期T或T的整数倍。交叉残差法表达式为：上式中，N为测量数据的总数，R(Δt)表示为原信号I(t)与延时Δt之后的信号I(t+Δt)两者之差平方的平均值，时延Δt的取值范围为[0,tmax/2]，且只能取采样间隔的整数倍。与自相关法相比，在数值上采取了平方的形式，使得峰谷之间的差异更大，能够考虑到目标光变曲线整体包络与细节性闪耀，与自相关方法相比结果较为准确。由于上述的两种方法时延的最小取值只能为采样间隔，这限制了上述方法的精度于局限数据的采样间隔。比如，对于不含任何噪声的数据，采样间隔为1秒，其精度只能达到秒，对目标旋转速率为4rves/min的数据，采取交叉残差法可以很精确的计算出峰值之间间距为15秒，其精度为1秒，而对于自旋速率6.5rves/min的数据，只能计算出间隔为9秒，实际约为9.231秒，受数据采样间隔的影响，计算误差较大。4)频率的选取。一是将频谱分析法得到的频率分为短周期谱和长周期谱，在短周期和长周期谱中分别选取能量最高的频率作为目标的转动频率，长短周期的区分没有明显的界限，且能量高的频率不一定是目标代表目标转动频率，由于该方法会出现倍频，倍频与倍频之后频率相近的频率混叠之后，其能量会出现比目标的旋转频率能量高的情况。二是由频谱分析法、自相关法和交叉残差法得到的不同频率的信号，通过峰值3σ准则进行筛选，3σ准则的基本思想是：假设随机误差服从正态分布，因此误差的绝对值主要集中在μ附近。3σ准则的公式为：P(|d|>3σ+μ)＝0.027P(|d|>2σ+μ)＝0.046P(|d|>σ+μ)＝0.317因此，当误差大于σ+μ时，则有68.3％的概率认为其是粗大误差；当误差大于2σ+μ时，则有95.4％的概率认为其为粗大误差；当误差大于3σ+μ时，则有97.3％的概率认为其为粗大误差。实际应用中，若某一频率处对应的峰值σ+μ，则有68.2％的概率认为其为主周期对应的尖峰，该频率则是旋转周期对应的频率。具体的选取方法为：(1)计算均方差与均值：均方差均值(2)求|di|(i＝1,2,3…)的最大值|d|max，若|d|max<σ+μ，则认为峰值不存在，否则进入下一步。(3)返回|d|max与其对应的频率i，回到步骤1。上述步骤中，di表示数据的幅值，N为数据的个数，为所有幅值的均方根，μ为幅值的均值。对数据进行频谱分析、自相关分析和交叉残差分析，根据所得结果峰值的横坐标以及峰值之间的距离对空间目标的旋转参数进行反演。对于峰值不明显、峰值间距离不等的处理结果，依据3σ准则进行峰值筛选。此种方法对于空间目标翻滚复杂的光度序列，频谱分析、自相关分析和交叉残差分析得到频率谱峰值差异均不明显，且峰值分布不集中，该方法适用性不强，同时间隔的选取依赖于峰值之间的距离，出现倍频的情况该方法将无法选取出目标的旋转速率。
技术实现思路
针对上述缺陷，本专利技术提供了一种基于光度序列同时估计空间失稳目标进动和自旋速率方法，该方法通过对光度序列进行变分模态分解(VariationalModeDecomposition，VMD)算法，将光度序列分解为不同频率的本征模态函数(IntrinsicModeFunction，IMF)，通过计算每个本征模态函数与光度序列的互信息(MutualInformation，MI)，进行目标进动和自旋频率的提取。本专利技术对目标的光度数据进行处理分析，可同时估计出目标的进动和自旋速率，从目标的光度数据中获取更多的目标状态信息。为达到上述目的，本专利技术通过以下技术方案来具体实现：本专利技术提供了一种基于光度序列同时估计空间失稳目标进动和自旋速率方法，该方法包括：步骤一、获取空间失稳目标的光度序列；步骤二、将光度序列变分模态分解为多个本征模态函数；步骤三、计算每个本征模态函数与光度序列的互信息；步骤四、抽取最大的两个互信息值对应的本征模态函数频率作为空间失稳目标的自旋频率和进动频率，将空间失稳目标的自旋频率和进动频率转化为空间失稳目标的进动和自旋速率。所述步骤一中，获取空间失稳目标的光度序列的方法包括：基于数值计算、实际观测和/或实验室模拟测量中的一种或多种方式获取空间失稳目标的光度序列。所述步骤二中，将光度序列变分模态分解为多个本征模态函数，包括：基于经典维纳滤波、希尔伯特变换和频率混合的方式将光度序列通过变分问题的构造和变分问题的求解两个过程变分模态分解为多个本征模态函数。所述变分问题的构造包括：假设每个模态是具有中心频率的有限带宽，变分问题描述为寻求K个本征模态函数uk(t)，使得每个模态的估计带宽之和最小，约束条件为各模态之和等于输入光度序列f，具体构造步骤如下：①通过希本文档来自技高网...

【技术保护点】
1.一种基于光度序列同时估计空间失稳目标进动和自旋速率方法，其特征在于，该方法包括：步骤一、获取空间失稳目标的光度序列；步骤二、将光度序列变分模态分解为多个所述本征模态函数；步骤三、计算每个本征模态函数与光度序列的互信息；步骤四、抽取最大的两个互信息值对应的本征模态函数频率作为空间失稳目标的自旋频率和进动频率，将空间失稳目标的自旋频率和进动频率转化为空间失稳目标的进动和自旋速率。

【技术特征摘要】
1.一种基于光度序列同时估计空间失稳目标进动和自旋速率方法，其特征在于，该方法包括：步骤一、获取空间失稳目标的光度序列；步骤二、将光度序列变分模态分解为多个所述本征模态函数；步骤三、计算每个本征模态函数与光度序列的互信息；步骤四、抽取最大的两个互信息值对应的本征模态函数频率作为空间失稳目标的自旋频率和进动频率，将空间失稳目标的自旋频率和进动频率转化为空间失稳目标的进动和自旋速率。2.如权利要求1所述的方法，其特征在于，所述步骤一中，获取空间失稳目标的光度序列的方法包括：基于数值计算、实际观测和/或实验室模拟测量中的一种或多种方式获取空间失稳目标的光度序列。3.如权利要求1或2所述的方法，其特征在于，所述步骤二中，将光度序列变分模态分解为多个本征模态函数，包括：基于经典维纳滤波、希尔伯特变换和频率混合的方式将光度序列通过变分问题的构造和变分问题的求解两个过程变分模态分解为多个本征模态函数。4.如权利要求3所述的方法，其特征在于，所述变分问题的构造包括：假设每个模态是具有中心频率的有限带宽，变分问题描述为寻求K个本征模态函数uk(t)，使得每个模态的估计带宽之和最小，约束条件为各模态之和等于输入光度序列f，具体构造步骤如下：①通过希尔伯特变换，得到每个本征模态函数uk(t)的解析信号，并得到其单边谱：②将各模态解析信号混合，预估中心频率将每个模态的频谱调制到相应的基频带：③计算解调信号梯度的平方L2范数，估计出各模态信号带宽，受约束的变分问题如下：其中，{uk}＝{u1,...,uk}、{ωk}＝{ω1,...,ωk}为变分模态分解到K个本征模态函数分量及本征模态函数中心频率的集合，空间失稳目标的转动频率包含在中心频率的集合中，f为输入的光度序列，δ(t)为冲击函数，*表示卷积，j为虚数单位，光度序列f长度为T，则t∈[0,T]，s.t.表示限制条件，符号表示变量对t求偏导数。5.如权利要求3所述的方法，其特征在于，所述变分问题的求解包括：①引入二次惩罚因子α(也成为平衡约束系数)和拉格朗日算子λ(t)，将约束性变分问题变为非约束性变分问题，扩展的拉格朗日表达式如下：②采用乘法算子交替方向法解决以上变分问题，通过交替更新和寻求扩展拉格朗日表达式的鞍点，其中的取值问题可表述为：其中，*表示卷积，n表示迭代...

【专利技术属性】
技术研发人员：李智，徐灿，张峰，谢杨峻，李鹏，周驰，汉京滨，张学文，周童，
申请(专利权)人：中国人民解放军战略支援部队航天工程大学，
类型：发明
国别省市：北京,11