CN109740209A - 高超声速飞行器参数在线辨识方法及使用其的力学模型 - Google Patents

高超声速飞行器参数在线辨识方法及使用其的力学模型 Download PDF

Info

Publication number
CN109740209A
CN109740209A CN201811560782.9A CN201811560782A CN109740209A CN 109740209 A CN109740209 A CN 109740209A CN 201811560782 A CN201811560782 A CN 201811560782A CN 109740209 A CN109740209 A CN 109740209A
Authority
CN
China
Prior art keywords
state
disturbance
estimate
priori
aircraft
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201811560782.9A
Other languages
English (en)
Inventor
杨志红
柳青
李瑜
徐宝华
佘智勇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Aerospace Technology Research Institute
Original Assignee
Beijing Aerospace Technology Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Aerospace Technology Research Institute filed Critical Beijing Aerospace Technology Research Institute
Priority to CN201811560782.9A priority Critical patent/CN109740209A/zh
Publication of CN109740209A publication Critical patent/CN109740209A/zh
Pending legal-status Critical Current

Links

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Traffic Control Systems (AREA)

Abstract

本发明提供了一种高超声速飞行器参数在线辨识方法及使用其的力学模型，该参数在线辨识方法包括：确定系统状态方程和观测方程，状态量包括不确定参数和确定参数；计算获取状态量先验估计值和观测量先验估计值；计算状态误差协方差矩阵先验估计值；计算获取卡尔曼增益更新方程；基于状态量先验估计值、观测量先验估计值和卡尔曼增益更新方程计算获取下一时刻的状态量估计值；步骤六，更新各变量，转至下一时刻，重复上述步骤，直至满足给定的计算次数，完成确定参数和不确定参数的辨识。应用本发明的技术方案，以解决现有技术中飞行器实际飞行过程中的气动、环境等不确定性对飞行器控制系统设计造成影响从而降低飞行器控制品质的技术问题。

Description

高超声速飞行器参数在线辨识方法及使用其的力学模型
技术领域
本发明涉及飞行器控制领域，尤其涉及一种高超声速飞行器参数在线辨识方法及使用其的力学模型。
背景技术
高超声速飞行器飞行包线大、飞行环境复杂、飞行高度及马赫数跨度范围大，气动特性有很剧烈的变化，气动参数难以准确预测，导致高超声速飞行器具有明显的“模型时变性”，这使得良好的制导系统和控制系统的设计成为一项重要难点。
目前，为解决高超声速飞行器气动参数剧烈变化对系统设计的影响，现有技术集中在研发增强控制系统和制导系统的鲁棒性上，即开发具有“快时变、强非线性、强耦合”的先进的控制与制导算法。但是，这类鲁棒性设计将表现出很大的保守性，牺牲了飞行器实现飞行任务的基本性能，不适用于高质量的飞行任务的完成。此外，在这样的先进控制律/制导律设计过程中，对飞行器相关参数的常值组合拉偏也难以覆盖可能随飞行状态参数快速变化的不确定因素，无法确保飞行器在极限恶劣情况的飞行状态。
发明内容
本发明提供了一种高超声速飞行器参数在线辨识方法及使用其的力学模型，解决了现有技术中难以对飞行器飞行过程中不确定参数进行在线辨识从而造成的飞行器制导系统和控制系统设计过于保守，难以完成高质量飞行任务的技术问题。
根据本发明的一方面，提供了一种高超声速飞行器参数在线辨识方法，高超声速飞行器参数在线辨识方法包括：步骤一，确定系统状态方程和观测方程，系统状态方程包括系统状态方程矩阵，观测方程包括观测方程矩阵，给定系统状态方程矩阵和观测方程矩阵的初值，给定系统状态方程中状态量及状态量的初值，状态量包括不确定参数和确定参数，根据状态量的初值计算获取状态量估计值初值以及状态误差协方差矩阵的初值，给定需要计算的计算次数；步骤二，基于当前时刻下的系统状态方程矩阵以及当前时刻下的状态量估计值计算获取状态量先验估计值，根据状态量先验估计值计算获取观测量先验估计值；步骤三，基于当前时刻下的系统状态方程矩阵以及当前时刻下的状态误差协方差矩阵，计算状态误差协方差矩阵先验估计值；步骤四，基于步骤三中的状态误差协方差矩阵先验估计值计算获取卡尔曼增益更新方程；步骤五，基于步骤二中的状态量先验估计值、观测量先验估计值和步骤四中的卡尔曼增益更新方程计算获取下一时刻的状态量估计值；步骤六，更新状态量估计值、系统状态方程矩阵、观测方程矩阵、观测量以及状态误差协方差矩阵，转至下一时刻，重复步骤二至步骤五，直至满足计算次数，完成对高超声速飞行器参数中确定参数和不确定参数的辨识。
进一步地，在步骤一中，不确定参数包括阻力系数扰动量、升力系数扰动量、密度扰动量、攻角扰动量或侧滑角扰动量。
进一步地，在步骤一中，观测量包括体轴轴向加速度、体轴法向加速度、体轴侧向加速度、动压或速度。
进一步地，在步骤一中，系统状态方程包括：
其中，r为地心到飞行器的距离，为r的导数，V为飞行器相对地球的速度，为V的导数，γ为航迹倾角，为γ的导数，θ为经度，为θ的导数，ψ为航迹偏角，为ψ的导数，φ为纬度，为φ的导数，T为发动机推力，α为攻角、β为侧滑角,L为无量纲阻力，D为无量纲阻力，σ为弹道倾角,C为侧向力,为实际飞行过程中阻力系数扰动量dCD的导数,为实际飞行过程中升力系数扰动量dCL的导数，为实际飞行过程中密度扰动量dρ的导数，Ω为弹体坐标系到速度坐标系的转换矩阵。
进一步地，在步骤一中，系统状态方程中的无量纲升力L、无量纲阻力D与升力系数扰动量dCL、阻力系数扰动量dCD以及密度扰动量dρ的关系为：
ρ＝ρ0+dρ
其中，Sref为参考面积，CL为升力系数，m为质量，CD为阻力系数，为阻力系数基准值，为升力系数基准值，ρ为大气密度，ρ0为大气密度基准值。
进一步地，在步骤一中，观测方程包括：
其中，z为观测量，axB为飞行器体轴的轴向加速度，ayB为飞行器体轴的法向加速度，azB为飞行器体轴的侧向加速度，Fx为阻力，Fy为升力，Fz侧向力，为轴向加速度的测量噪声，为法向加速度的测量噪声，为侧向加速度的测量噪声，νq为动压测量噪声，νV为速度测量噪声，q为动压，V真为真实速度，V为测量速度。
进一步地，在步骤二中，采用如下公式计算获取状态量先验估计值：
其中，为地心到飞行器的距离的先验估计值，为经度的先验估计值，为纬度的先验估计值，为速度的先验估计值，为航迹倾角的先验估计值，为航迹偏角的先验估计值，为阻力系数扰动量的先验估计值，为升力系数扰动量的先验估计值，为密度扰动量的先验估计值，为当前时刻的地心到飞行器的估计距离，为当前时刻的估计经度，为当前时刻的估计纬度，为当前时刻的估计速度，为当前时刻的估计航迹倾角，为当前时刻的估计航迹偏角，为当前时刻的估计阻力系数扰动量，为当前时刻的估计升力系数扰动量，为当前时刻的估计密度扰动量，Δt＝tk+1-tk为滤波步长，tk为当前时刻，tk+1为下一时刻。
进一步地，在步骤三中，采用如下公式计算状态误差协方差矩阵先验估计值：
P(k+1/k)＝Φ(k)P(k/k)ΦT(k)+Qk
其中，P(k+1/k)为状态误差协方差矩阵先验估计值，Φ(k)为当前时刻下的系统状态方程矩阵，P(k/k)为当前时刻下的状态误差协方差矩阵，Qk由公式E[w(k)wT(k)]＝Qk计算获取，w(k)为过程噪声。
进一步地，在步骤五中，采用如下公式计算下一时刻的状态量估计值：
其中，为下一时刻的状态量估计值，为下一时刻地心到飞行器的距离估计值，为下一时刻的经度估计值，为下一时刻的纬度的估计值，为下一时刻的速度估计值，为下一时刻的航迹倾角估计值，为下一时刻的航迹偏角估计值，为下一时刻的阻力系数扰动量估计值，为下一时刻的升力系数扰动量估计值，为下一时刻的密度扰动量估计值，为地心到飞行器的距离的先验估计值，为经度的先验估计值，为纬度的先验估计值，为速度的先验估计值，为航迹倾角的先验估计值，为航迹偏角的先验估计值，为阻力系数扰动量的先验估计值，为升力系数扰动量的先验估计值，为密度扰动量的先验估计值，K(k+1)为卡尔曼增益更新方程，axB(k+1/k+1)为tk+1时刻体轴轴向加速度，ayB(k+1/k+1)为tk+1时刻体轴法向加速度，azB(k+1/k+1)为tk+1时刻体轴侧向加速度，q(k+1/k+1)为tk+1时刻动压，V(k+1/k+1)为tk+1时刻速度，为体轴轴向加速度先验估计值，为体轴法向加速度先验估计值，为体轴侧向加速度先验估计值，为动压先验估计值，为速度先验估计值。
根据本发明的另一方面，提供了一种高超声速飞行器力学模型，高超声速飞行器力学模型使用了如上所述的高超声速飞行器参数辨识方法进行建模。
应用本发明的技术方案，通过基于EKF的高超声速飞行器在线建模方法，对飞行器实际飞行过程中的气动、环境等不确定性参数进行在线辨识，所得辨识结果可用于在线构建更为精准的飞行器力学模型和动力学模型，解决了飞行器实际飞行过程中的气动、环境等不确定性对控制系统设计影响的技术问题，避免了常见的鲁棒控制方法因考虑鲁棒性而导致的性能的牺牲，提高了发生较大模型偏差和干扰时飞行器自适应控制性能，确保控制品质不会过度降低。
附图说明
所包括的附图用来提供对本发明实施例的进一步的理解，其构成了说明书的一部分，用于例示本发明的实施例，并与文字描述一起来阐释本发明的原理。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1示出了根据本发明的具体实施例提供的高超声速飞行器参数在线辨识方法的流程图；
图2示出了根据本发明的具体实施例提供的升力系数扰动量dCL辨识结果示意图；
图3示出了根据本发明的具体实施例提供的阻力系数扰动量dCD辨识结果示意图；
图4示出了根据本发明的具体实施例提供的密度扰动量dρ辨识结果示意图。
具体实施方式
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本发明的范围。同时，应当明白，为了便于描述，附图中所示出的各个部分的尺寸并不是按照实际的比例关系绘制的。对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为授权说明书的一部分。在这里示出和讨论的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它示例可以具有不同的值。应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步讨论。
如图1至图4所示，根据本发明的具体实施例提供了一种高超声速飞行器参数在线辨识方法，高超声速飞行器参数在线辨识方法包括：步骤一，确定系统状态方程和观测方程，系统状态方程包括系统状态方程矩阵，观测方程包括观测方程矩阵，给定系统状态方程矩阵和观测方程矩阵的初值，给定系统状态方程中状态量及状态量的初值，状态量包括不确定参数和确定参数，根据状态量的初值计算获取状态量估计值初值以及状态误差协方差矩阵的初值，给定需要计算的计算次数；步骤二，基于当前时刻下的系统状态方程矩阵以及当前时刻下的状态量估计值计算获取状态量先验估计值，根据状态量先验估计值计算获取观测量先验估计值；步骤三，基于当前时刻下的系统状态方程矩阵以及当前时刻下的状态误差协方差矩阵，计算状态误差协方差矩阵先验估计值；步骤四，基于步骤三中的状态误差协方差矩阵先验估计值计算获取卡尔曼增益更新方程；步骤五，基于步骤二中的状态量先验估计值、观测量先验估计值和步骤四中的卡尔曼增益更新方程计算获取下一时刻的状态量估计值；步骤六，更新状态量估计值、系统状态方程矩阵、观测方程矩阵、观测量以及状态误差协方差矩阵，转至下一时刻，重复步骤二至步骤五，直至满足计算次数，完成对高超声速飞行器参数中确定参数和不确定参数的辨识。
应用此种配置方式，通过基于EKF的高超声速飞行器在线建模方法，对飞行器实际飞行过程中的气动、环境等不确定性参数进行在线辨识，所得辨识结果可用于在线构建更为精准的飞行器力学模型和动力学模型，解决了飞行器实际飞行过程中的气动、环境等不确定性对控制系统设计的影响，避免了常见的鲁棒控制方法因考虑鲁棒性而导致的性能的牺牲，提高了发生较大模型偏差和干扰时飞行器自适应控制性能，确保控制品质不会过度降低。
进一步地，为了对飞行器实际飞行过程中的升力、阻力和密度的不确定性参数进行在线辨识，在步骤一中，不确定参数包括阻力系数扰动量、升力系数扰动量、密度扰动量、攻角扰动量或侧滑角扰动量。
应用此种配置方式，通过对阻力系数扰动量、升力系数扰动量和密度扰动量进行在线辨识，可以有效地识别出飞行器在飞行过程中阻力、升力和大气密度的变化情况，将辨识结果应用于力学模型能够提高模型构建的精确性。
进一步地，为了方便传感器进行测量，在步骤一中，观测量包括体轴轴向加速度、体轴法向加速度、体轴侧向加速度、动压或速度。
应用此种配置方式，通过将观测量具体选取为体轴轴向加速度、体轴法向加速度、体轴侧向加速度、动压和速度，能够方便传感器对观测量进行测量，提高测量的精确程度，从而增加参数在线辨识的精确程度。
进一步地，为了对飞行器动力学模型进行具体描述，系统状态方程包括：
其中，r为地心到飞行器的距离，为r的导数，V为飞行器相对地球的速度，为V的导数，γ为航迹倾角，为γ的导数，θ为经度，为θ的导数，ψ为航迹偏角，为ψ的导数，φ为纬度，为φ的导数，T为发动机推力，α为攻角、β为侧滑角,L为无量纲阻力，D为无量纲阻力，σ为弹道倾角,C为侧向力,为实际飞行过程中阻力系数扰动量dCD的导数,为实际飞行过程中升力系数扰动量dCL的导数，为实际飞行过程中密度扰动量dρ的导数，Ω为弹体坐标系到速度坐标系的转换矩阵。
应用此种配置方程，通过采用上述系统状态方程引入阻力系数扰动量、升力系数扰动量和密度扰动量，对飞行器飞行过程中的扰动量进行辨识，将飞行器飞行过程中的扰动量纳入系统状态方程，从而能够更精确的在线描述飞行器的动力学规律。
进一步地，为了获取升力系数扰动量dCL、阻力系数扰动量dCD以及密度扰动量dρ，在步骤一中，系统状态方程中的无量纲升力L、无量纲阻力D与升力系数扰动量dCL、阻力系数扰动量dCD以及密度扰动量dρ的关系为：
ρ＝ρ0+dρ
其中，Sref为参考面积，CL为升力系数，m为质量，CD为阻力系数，为阻力系数基准值，为升力系数基准值，ρ为大气密度，ρ0为大气密度基准值。
应用此种配置方式，通过建立系统状态方程中的无量纲升力L、无量纲阻力D与升力系数扰动量dCL、阻力系数扰动量dCD以及密度扰动量dρ的关系，能够通过计算获取升力系数扰动量dCL、阻力系数扰动量dCD以及密度扰动量dρ，完成各扰动量的参数辨识。
进一步地，为了获取飞行器观测量和状态量的关联，观测方程包括：
其中，z为观测量，axB为飞行器体轴的轴向加速度，ayB为飞行器体轴的法向加速度，azB为飞行器体轴的侧向加速度，Fx为阻力，Fy为升力，Fz侧向力，为轴向加速度的测量噪声，为法向加速度的测量噪声，为侧向加速度的测量噪声，νq为动压测量噪声，νV为速度测量噪声，q为动压，V真为真实速度，V为测量速度。
应用此种配置方式，通过建立观测量和状态量的关系构建飞行器观测方程，通过系统状态方程和观测方程完成对高超声速飞行器不确定参数的辨识。
在完成步骤一构建系统状态方程和观测方程后，需要计算获取状态量先验估计值，为了完成对地心到飞行器的距离、经度、纬度、速度、航迹倾角、航迹偏角、阻力系数扰动量、升力系数扰动量和密度扰动量的在线辨识，步骤二中，采用如下公式计算获取状态量先验估计值：
其中，为地心到飞行器的距离的先验估计值，为经度的先验估计值，为纬度的先验估计值，为速度的先验估计值，为航迹倾角的先验估计值，为航迹偏角的先验估计值，为阻力系数扰动量的先验估计值，为升力系数扰动量的先验估计值，为密度扰动量的先验估计值，为当前时刻的地心到飞行器的估计距离，为当前时刻的估计经度，为当前时刻的估计纬度，为当前时刻的估计速度，为当前时刻的估计航迹倾角，为当前时刻的估计航迹偏角，为当前时刻的估计阻力系数扰动量，为当前时刻的估计升力系数扰动量，为当前时刻的估计密度扰动量，Δt＝tk+1-tk为滤波步长，tk为当前时刻，tk+1为下一时刻。
应用此种配置方式，通过上述状态方程进行计算以获取状态量先验估计值 完成对各状态量的先验估计。
在完成步骤二构建系统状态方程和观测方程后，需要计算获取状态误差协方差矩阵先验估计值，在步骤三中，采用如下公式计算获取状态误差协方差矩阵先验估计值：
P(k+1/k)＝Φ(k)P(k/k)ΦT(k)+Qk
其中，P(k+1/k)为状态误差协方差矩阵先验估计值，Φ(k)为当前时刻下的系统状态方程矩阵，P(k/k)为当前时刻下的状态误差协方差矩阵，Qk由公式E[w(k)wT(k)]＝Qk计算获取，w(k)为过程噪声。
应用此种配置方式，通过当前时刻下的系统状态方程矩阵以及当前时刻下的状态误差协方差矩阵，完成对状态误差协方差矩阵先验估计值的计算。作为本发明的一个具体实施例，采用白噪声为过程噪声。
在完成步骤四计算获取卡尔曼增益更新方程后，需要计算获取下一时刻的状态量估计值，为了完成基于状态量先验估计值、观测量先验估计值和卡尔曼增益更新方程计算获取下一时刻的状态量估计值，在步骤五中，采用如下公式计算下一时刻的状态量估计值：
其中，为下一时刻的状态量估计值，为下一时刻地心到飞行器的距离估计值，为下一时刻的经度估计值，为下一时刻的纬度的估计值，为下一时刻的速度估计值，为下一时刻的航迹倾角估计值，为下一时刻的航迹偏角估计值，为下一时刻的阻力系数扰动量估计值，为下一时刻的升力系数扰动量估计值，为下一时刻的密度扰动量估计值，为地心到飞行器的距离的先验估计值，为经度的先验估计值，为纬度的先验估计值，为速度的先验估计值，为航迹倾角的先验估计值，为航迹偏角的先验估计值，为阻力系数扰动量的先验估计值，为升力系数扰动量的先验估计值，为密度扰动量的先验估计值，K(k+1)为卡尔曼增益更新方程，axB(k+1/k+1)为tk+1时刻体轴轴向加速度，ayB(k+1/k+1)为tk+1时刻体轴法向加速度，azB(k+1/k+1)为tk+1时刻体轴侧向加速度，q(k+1/k+1)为tk+1时刻动压，V(k+1/k+1)为tk+1时刻速度，为体轴轴向加速度先验估计值，为体轴法向加速度先验估计值，为体轴侧向加速度先验估计值，为动压先验估计值，为速度先验估计值。
应用此种配置方式，根据观测量实际值与观测量先验估计值的误差，进一步结合状态量先验估计值，从而得到下一时刻状态量的估计值，完成对下一时刻状态量估计值的计算。
根据本发明的另一方面，提供了一种高超声速飞行器力学模型，高超声速飞行器力学模型使用了如上所述的高超声速飞行器参数辨识方法进行建模。应用此种配置方式，将不确定参数辨识结果应用于飞行器力学模型，避免因考虑不确定参数波动造成的飞行器性能的牺牲，有利于提高飞行器的控制品质。
为了对本发明有进一步的了解，下面结合图1至图4对本发明的高超声速飞行器参数在线辨识方法进行详细说明。
如图1所示，该方法具体包括以下步骤。
步骤一，确定系统状态方程和观测方程，给定状态量、系统状态方程矩阵以及观测方程矩阵的初值。
系统状态方程和观测方程如下：
x(k+1)＝Φ(k)x(k)+τ(k)w(k)
z(k)＝H(k)x(k)+v(k)
系统状态方程中，x(k+1)为下一时刻状态量，x(k)为当前时刻状态量，τ(k)为环境参数矩阵，在本实施例中，τ(k)为单位向量。观测方程中，z(k)为当前时刻观测量，H(k)为当前时刻观测方程系数矩阵，ν(k)为当前时刻量测噪声，本实施例中量测噪声为白噪声。
给定系统状态方程中的状态量的初值、系统状态方程初值以及观测方程系数矩阵初值，给定需要计算的计算次数。
根据状态量的初值x(0)计算得到初始时刻的状态量估计值
其中，m0为给定的初值期望。
根据状态量的初值x(0)以及第一次计算的状态量估计值计算获取状态误差协方差矩阵的初值P(0/0)：
本实施例中，选取上述初值以得到无偏估计。
具体的，在本实施例中，状态量包括确定参数r、θ、φ、V、γ和ψ，此外状态量还包括不确定参数dCD、dCL和dρ，系统状态方程的具体表达式为：
其中，Ω为弹体坐标系到速度坐标系的转换矩阵,其计算公式为：
系统状态方程中的升力L、阻力D与升力系数扰动量dCL、阻力系数扰动量dCD以及密度扰动量dρ的关系为：
ρ＝ρ0+dρ
在本实施例中，系统观测方程为：
步骤二，基于当前时刻下的系统状态方程矩阵Φ(k)以及当前时刻下的状态量估计值计算状态量先验估计值采用公式如下：
其中，为状态量先验估计值，Φ(k)＝I+F(k)Δt，I为单位矩阵，F(k)为状态方程的雅克比矩阵，在本实施例中，状态量先验估计值计算表达式为 具体表达式为：
其中，为当前时刻的估计速度，为当前时刻的估计航迹倾角，为当前时刻的估计航迹偏角，为当前时刻的地心到飞行器的估计距离，为当前时刻的估计纬度。
将上述带入公式得到扩展卡尔曼滤波的先验状态更新方程用以计算状态量的先验估计值，在本实施例中，先验状态更新方程为：
根据状态量的先验估计值计算得到观测量先验估计值
步骤三，基于当前时刻下的系统状态方程矩阵Φ(k)以及当前时刻下的状态误差协方差矩阵P(k/k)，计算状态误差协方差矩阵先验估计值P(k+1/k)，具体可采用下述公式进行计算。
P(k+1/k)＝Φ(k)P(k/k)ΦT(k)+τ(k)QkτT(k)
在本实施例中，τ(k)取为单位向量：
P(k+1/k)＝Φ(k)P(k/k)ΦT(k)+Qk。
步骤四，基于步骤三中的状态误差协方差矩阵先验估计值P(k+1/k)计算卡尔曼增益更新方程K(k+1)：
K(k+1)＝P(k+1/k)HT(k+1)[H(k+1)P(k+1/k)HT(k+1)+R(k+1)]-1
其中量测噪声矩阵R(k+1)＝E[ν(k+1)νT(k+1)]，ν(k+1)为tk+1时刻的量测噪声，H(k+1)为状态量先验估计值的观测方程系数矩阵的雅克比矩阵。
步骤五，基于步骤二中的状态量先验误差估计值观测量先验估计值以及步骤四中的卡尔曼增益更新方程K(k+1)计算获取下一时刻的状态量估计值采用下述公式进行计算：
其中，z(k+1/k+1)为tk+1时刻的观测量，其数值可由传感器测得。本实施例中，z(k+1/k+1)的表达式为：
其中，axB(k+1/k+1)为tk+1时刻体轴轴向加速度，ayB(k+1/k+1)为tk+1时刻体轴法向加速度，azB(k+1/k+1)为tk+1时刻体轴侧向加速度，q(k+1/k+1)为tk+1时刻动压，V(k+1/k+1)为tk+1时刻速度，V(k/k)为tk时刻速度。
扩展卡尔曼滤波中的后验状态更新方程为：
步骤六，更新状态量估计值、系统状态方程矩阵、观测方程矩阵、观测量以及状态误差协方差矩阵，转至下一时刻，重复步骤二至步骤五，直至满足计算次数，完成对高超声速飞行器参数中确定参数和不确定参数dCD、dCL和dρ的辨识。
其中tk+1时刻状态误差协方差矩阵采用下述公式：
P(k+1/k+1)＝P(k+1/k)-P(k+1/k)HT(k+1)
[H(k+1)P(k+1/k)HT(k+1)+R(k+1)]-1H(k+1)P(k+1/k)
在本实施例中，升力系数、阻力系数、密度分别存在标称值20％、-20％、10％的不确定扰动，如图2至图4所示。
本发明的高超声速飞行器参数在线辨识方法及使用其的力学模型，通过基于EKF(扩展卡尔曼滤波)的高超声速飞行器在线建模方法，实现对飞行器实际飞行过程中高度、经度、纬度、速度、航迹倾角或航迹偏角等确定参数以及升力系数扰动量、阻力系数扰动量或密度扰动量等不确定性参数进行在线联合辨识，该方法能够较好、较快地在线辨识出飞行器实际飞行过程中存在气动、环境和风场扰动所造成的飞行器气动系数、密度、攻角以及侧滑角中的不确定性，进一步可根据辨识结果完成高超声速飞行器的在线建模，构建更为精准的飞行器力学模型和动力学模型，为自适应控制律的设计奠定基础，避免了常见的鲁棒控制方法因考虑鲁棒性而导致的性能的牺牲，实现了发生较大模型偏差和干扰时飞行器自适应控制，确保控制品质不会过度降低。
为了便于描述，在这里可以使用空间相对术语，如“在……之上”、“在……上方”、“在……上表面”、“上面的”等，用来描述如在图中所示的一个器件或特征与其他器件或特征的空间位置关系。应当理解的是，空间相对术语旨在包含除了器件在图中所描述的方位之外的在使用或操作中的不同方位。例如，如果附图中的器件被倒置，则描述为“在其他器件或构造上方”或“在其他器件或构造之上”的器件之后将被定位为“在其他器件或构造下方”或“在其他器件或构造之下”。因而，示例性术语“在……上方”可以包括“在……上方”和“在……下方”两种方位。该器件也可以其他不同方式定位(旋转90度或处于其他方位)，并且对这里所使用的空间相对描述作出相应解释。
此外，需要说明的是，使用“第一”、“第二”等词语来限定零部件，仅仅是为了便于对相应零部件进行区别，如没有另行声明，上述词语并没有特殊含义，因此不能理解为对本发明保护范围的限制。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种高超声速飞行器参数在线辨识方法，其特征在于，所述高超声速飞行器参数在线辨识方法包括：
步骤一，确定系统状态方程和观测方程，所述系统状态方程包括系统状态方程矩阵，所述观测方程包括观测方程矩阵，给定系统状态方程矩阵和观测方程矩阵的初值，给定所述系统状态方程中状态量及所述状态量的初值，所述状态量包括不确定参数和确定参数，根据所述状态量的初值计算获取状态量估计值初值以及状态误差协方差矩阵的初值，给定需要计算的计算次数；
步骤二，基于当前时刻下的所述系统状态方程矩阵以及当前时刻下的状态量估计值计算获取状态量先验估计值，根据所述状态量先验估计值计算获取观测量先验估计值；
步骤三，基于当前时刻下的所述系统状态方程矩阵以及当前时刻下的状态误差协方差矩阵，计算状态误差协方差矩阵先验估计值；
步骤四，基于所述步骤三中的所述状态误差协方差矩阵先验估计值计算获取卡尔曼增益更新方程；
步骤五，基于所述步骤二中的所述状态量先验估计值、所述观测量先验估计值和所述步骤四中的所述卡尔曼增益更新方程计算获取下一时刻的状态量估计值；
步骤六，更新所述状态量估计值、所述系统状态方程矩阵、所述观测方程矩阵、所述观测量以及所述状态误差协方差矩阵，转至下一时刻，重复步骤二至步骤五，直至满足所述计算次数，完成对高超声速飞行器参数中确定参数和不确定参数的辨识。
2.根据权利要求1所述的高超声速飞行器参数在线辨识方法，其特征在于，在所述步骤一中，所述不确定参数包括阻力系数扰动量、升力系数扰动量、密度扰动量、攻角扰动量或侧滑角扰动量。
3.根据权利要求1所述的高超声速飞行器参数在线辨识方法，其特征在于，在所述步骤一中，所述观测方程中的所述观测量包括体轴轴向加速度、体轴法向加速度、体轴侧向加速度、动压或速度。
4.根据权利要求3所述的高超声速飞行器参数在线辨识方法，其特征在于，在所述步骤一中，所述系统状态方程包括：
其中，r为地心到飞行器的距离，为r的导数，V为飞行器相对地球的速度，为V的导数，γ为航迹倾角，为γ的导数，θ为经度，为θ的导数，ψ为航迹偏角，为ψ的导数，φ为纬度，为φ的导数，T为发动机推力，α为攻角、β为侧滑角,L为无量纲阻力，D为无量纲阻力，σ为弹道倾角,C为侧向力,为实际飞行过程中阻力系数扰动量dCD的导数,为实际飞行过程中升力系数扰动量dCL的导数，为实际飞行过程中密度扰动量dρ的导数，Ω为弹体坐标系到速度坐标系的转换矩阵。
5.根据权利要求4所述的高超声速飞行器参数在线辨识方法，其特征在于，在所述步骤一中，所述系统状态方程中的所述无量纲升力L、所述无量纲阻力D与所述升力系数扰动量dCL、所述阻力系数扰动量dCD以及所述密度扰动量dρ的关系为：
ρ＝ρ0+dρ
其中，Sref为参考面积，CL为升力系数，m为质量，CD为阻力系数，为阻力系数基准值，为升力系数基准值，ρ为大气密度，ρ0为大气密度基准值。
6.根据权利要求5所述的高超声速飞行器参数在线辨识方法，其特征在于，在所述步骤一中，所述观测方程包括：
其中，z为观测量，axB为飞行器体轴的轴向加速度，ayB为飞行器体轴的法向加速度，azB为飞行器体轴的侧向加速度，Fx为阻力，Fy为升力，Fz侧向力，为轴向加速度的测量噪声，为法向加速度的测量噪声，为侧向加速度的测量噪声，νq为动压测量噪声，νV为速度测量噪声，q为动压，V真为真实速度，V为测量速度。
7.根据权利要求6所述的高超声速飞行器参数在线辨识方法，其特征在于，在所述步骤二中，采用如下公式计算获取所述状态量先验估计值：
其中，为地心到飞行器的距离的先验估计值，为经度的先验估计值，为纬度的先验估计值，为速度的先验估计值，为航迹倾角的先验估计值，为航迹偏角的先验估计值，为阻力系数扰动量的先验估计值，为升力系数扰动量的先验估计值，为密度扰动量的先验估计值，为当前时刻的地心到飞行器的估计距离，为当前时刻的估计经度，为当前时刻的估计纬度，为当前时刻的估计速度，为当前时刻的估计航迹倾角，为当前时刻的估计航迹偏角，为当前时刻的估计阻力系数扰动量，为当前时刻的估计升力系数扰动量，为当前时刻的估计密度扰动量，Δt＝tk+1-tk为滤波步长，tk为当前时刻，tk+1为下一时刻。
8.根据权利要求1所述的高超声速飞行器参数在线辨识方法，其特征在于，在所述步骤三中，采用如下公式计算所述状态误差协方差矩阵先验估计值：
P(k+1/k)＝Φ(k)P(k/k)ΦT(k)+Qk
其中，P(k+1/k)为状态误差协方差矩阵先验估计值，Φ(k)为当前时刻下的系统状态方程矩阵，P(k/k)为当前时刻下的状态误差协方差矩阵，Qk由公式E[w(k)wT(k)]＝Qk计算获取，w(k)为过程噪声。
9.根据权利要求7所述的高超声速飞行器参数在线辨识方法，其特征在于，在所述步骤五中，采用如下公式计算下一时刻的状态量估计值：
其中，为下一时刻的状态量估计值，为下一时刻地心到飞行器的距离估计值，为下一时刻的经度估计值，为下一时刻的纬度的估计值，为下一时刻的速度估计值，为下一时刻的航迹倾角估计值，为下一时刻的航迹偏角估计值，为下一时刻的阻力系数扰动量估计值，为下一时刻的升力系数扰动量估计值，为下一时刻的密度扰动量估计值，为地心到飞行器的距离的先验估计值，为经度的先验估计值，为纬度的先验估计值，为速度的先验估计值，为航迹倾角的先验估计值，为航迹偏角的先验估计值，为阻力系数扰动量的先验估计值，为升力系数扰动量的先验估计值，为密度扰动量的先验估计值，K(k+1)为所述卡尔曼增益更新方程，axB(k+1/k+1)为tk+1时刻体轴轴向加速度，ayB(k+1/k+1)为tk+1时刻体轴法向加速度，azB(k+1/k+1)为tk+1时刻体轴侧向加速度，q(k+1/k+1)为tk+1时刻动压，V(k+1/k+1)为tk+1时刻速度，为体轴轴向加速度先验估计值，为体轴法向加速度先验估计值，为体轴侧向加速度先验估计值，为动压先验估计值，为速度先验估计值。
10.一种高超声速飞行器力学模型，其特征在于，所述高超声速飞行器力学模型使用如权利要求1至9中任一项所述的高超声速飞行器参数辨识方法进行建模。
CN201811560782.9A 2018-12-20 2018-12-20 高超声速飞行器参数在线辨识方法及使用其的力学模型 Pending CN109740209A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811560782.9A CN109740209A (zh) 2018-12-20 2018-12-20 高超声速飞行器参数在线辨识方法及使用其的力学模型

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811560782.9A CN109740209A (zh) 2018-12-20 2018-12-20 高超声速飞行器参数在线辨识方法及使用其的力学模型

Publications (1)

Publication Number Publication Date
CN109740209A true CN109740209A (zh) 2019-05-10

Family

ID=66360777

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811560782.9A Pending CN109740209A (zh) 2018-12-20 2018-12-20 高超声速飞行器参数在线辨识方法及使用其的力学模型

Country Status (1)

Country Link
CN (1) CN109740209A (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110162071A (zh) * 2019-05-24 2019-08-23 北京控制工程研究所 一种高超声速飞行器再入末段姿态控制方法及系统
CN111125971A (zh) * 2019-12-26 2020-05-08 北京航空航天大学 一种吸气式高超声速飞行器推力不确定性确定方法
CN111425304A (zh) * 2020-04-23 2020-07-17 南京航空航天大学 基于复合模型预测控制的航空发动机直接推力控制方法
CN111443726A (zh) * 2020-03-02 2020-07-24 北京空天技术研究所 基于飞行试验数据的弹道重构方法
CN112668104A (zh) * 2021-01-04 2021-04-16 中国人民解放军96901部队22分队 一种高超声速飞行器气动参数在线辨识方法
CN113805602A (zh) * 2021-10-23 2021-12-17 北京航空航天大学 一种考虑阵风影响的无人机飞行高度控制方法
CN113919194A (zh) * 2021-09-07 2022-01-11 中国空气动力研究与发展中心计算空气动力研究所 一种基于滤波误差法的非线性飞行动力学系统辨识方法
CN117521561A (zh) * 2024-01-03 2024-02-06 中国人民解放军国防科技大学 巡航飞行器的气动力与推力在线预示方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015033503A1 (ja) * 2013-09-05 2015-03-12 カルソニックカンセイ株式会社 推定装置及び推定方法
CN105973238A (zh) * 2016-05-09 2016-09-28 郑州轻工业学院 一种基于范数约束容积卡尔曼滤波的飞行器姿态估计方法
CN109033493A (zh) * 2018-06-01 2018-12-18 南京理工大学 基于无迹卡尔曼滤波的辨识高速旋转弹气动参数滤波方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015033503A1 (ja) * 2013-09-05 2015-03-12 カルソニックカンセイ株式会社 推定装置及び推定方法
CN105973238A (zh) * 2016-05-09 2016-09-28 郑州轻工业学院 一种基于范数约束容积卡尔曼滤波的飞行器姿态估计方法
CN109033493A (zh) * 2018-06-01 2018-12-18 南京理工大学 基于无迹卡尔曼滤波的辨识高速旋转弹气动参数滤波方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
吴楠等: "高超声速滑翔再入飞行器弹道估计的自适应卡尔曼滤波", 《航空学报》 *
戚磊等: "基于扩展卡尔曼滤波的超燃冲压发动机空气质量流量估计方法", 《战术导弹技术》 *
赵志刚等: "基于迭代滤波的上升段不确定参数辨识", 《飞行力学》 *

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110162071A (zh) * 2019-05-24 2019-08-23 北京控制工程研究所 一种高超声速飞行器再入末段姿态控制方法及系统
CN110162071B (zh) * 2019-05-24 2022-04-22 北京控制工程研究所 一种高超声速飞行器再入末段姿态控制方法及系统
CN111125971A (zh) * 2019-12-26 2020-05-08 北京航空航天大学 一种吸气式高超声速飞行器推力不确定性确定方法
CN111443726A (zh) * 2020-03-02 2020-07-24 北京空天技术研究所 基于飞行试验数据的弹道重构方法
CN111443726B (zh) * 2020-03-02 2023-08-15 北京空天技术研究所 基于飞行试验数据的弹道重构方法
CN111425304B (zh) * 2020-04-23 2021-01-12 南京航空航天大学 基于复合模型预测控制的航空发动机直接推力控制方法
CN111425304A (zh) * 2020-04-23 2020-07-17 南京航空航天大学 基于复合模型预测控制的航空发动机直接推力控制方法
CN112668104A (zh) * 2021-01-04 2021-04-16 中国人民解放军96901部队22分队 一种高超声速飞行器气动参数在线辨识方法
CN112668104B (zh) * 2021-01-04 2022-11-29 中国人民解放军96901部队22分队 一种高超声速飞行器气动参数在线辨识方法
CN113919194A (zh) * 2021-09-07 2022-01-11 中国空气动力研究与发展中心计算空气动力研究所 一种基于滤波误差法的非线性飞行动力学系统辨识方法
CN113805602A (zh) * 2021-10-23 2021-12-17 北京航空航天大学 一种考虑阵风影响的无人机飞行高度控制方法
CN113805602B (zh) * 2021-10-23 2022-04-08 北京航空航天大学 一种考虑阵风影响的无人机飞行高度控制方法
CN117521561A (zh) * 2024-01-03 2024-02-06 中国人民解放军国防科技大学 巡航飞行器的气动力与推力在线预示方法
CN117521561B (zh) * 2024-01-03 2024-03-19 中国人民解放军国防科技大学 巡航飞行器的气动力与推力在线预示方法

Similar Documents

Publication Publication Date Title
CN109740209A (zh) 高超声速飞行器参数在线辨识方法及使用其的力学模型
CN108827299B (zh) 一种基于改进四元数二阶互补滤波的飞行器姿态解算方法
Sabet et al. A low-cost dead reckoning navigation system for an AUV using a robust AHRS: Design and experimental analysis
Zhong et al. Sensor fault detection and diagnosis for an unmanned quadrotor helicopter
CN106354901B (zh) 一种运载火箭质量特性及动力学关键参数在线辨识方法
CN109708663B (zh) 基于空天飞机sins辅助的星敏感器在线标定方法
CN105783943A (zh) 一种基于无迹卡尔曼滤波的极区环境下舰船大方位失准角传递对准方法
Morelli et al. Real-time dynamic modeling: data information requirements and flight-test results
CN103852081B (zh) 用于大气数据/捷联惯导组合导航系统的真空速解算方法
CN108613674A (zh) 一种基于自适应差分进化bp神经网络的姿态误差抑制方法
Gao et al. Rapid alignment method based on local observability analysis for strapdown inertial navigation system
Sabet et al. Experimental analysis of a low-cost dead reckoning navigation system for a land vehicle using a robust AHRS
CN101122637A (zh) 一种sar运动补偿用sins/gps组合导航自适应降维滤波方法
CN108871319B (zh) 一种基于地球重力场与地磁场序贯修正的姿态解算方法
CN110132287A (zh) 一种基于极限学习机网络补偿的卫星高精度联合定姿方法
CN111189442A (zh) 基于cepf的无人机多源导航信息状态预测方法
CN109000639B (zh) 乘性误差四元数地磁张量场辅助陀螺的姿态估计方法及装置
Yu et al. Stochastic observability-based analytic optimization of SINS multiposition alignment
Cho et al. Airflow angle and wind estimation using GPS/INS navigation data and airspeed
CN117875090B (zh) 一种考虑风干扰的固定翼无人机增量元飞行气动建模方法
CN113919194B (zh) 一种基于滤波误差法的非线性飞行动力学系统辨识方法
CN103954288B (zh) 一种卫星姿态确定系统精度响应关系确定方法
CN107063300A (zh) 一种基于反演的水下导航系统动力学模型中扰动估计方法
CN111412887B (zh) 一种基于卡尔曼滤波的攻角、侧滑角辨识方法
Hao et al. Rapid transfer alignment based on unscented Kalman filter

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20190510

RJ01 Rejection of invention patent application after publication