一种基于卡尔曼滤波的攻角、侧滑角辨识方法与流程

本发明涉及航空飞行器飞行过程中的攻角和侧滑角测量技术领域，尤其涉及一种基于卡尔曼滤波的攻角、侧滑角辨识方法。

背景技术：

目前，对于飞航导弹或无人机等飞行器而言，普遍不装备攻角、侧滑角传感器，因此缺乏直接测量攻角、侧滑角的手段。而在飞行控制、发动机控制中，飞行器的攻角、侧滑角又是一个重要的参控量。

目前，控制专业常使用过载替代攻角、侧滑角。但是，影响过载的因素很多，除攻侧角外，还有动压、飞行速度等因素，因此精度较差，可靠性低。

技术实现要素：

本发明提供了一种基于卡尔曼滤波的攻角、侧滑角辨识方法，能够解决现有技术中在无攻角、侧滑角测量装置时无法精确测量攻角和侧滑角的技术问题。

本发明提供了一种基于卡尔曼滤波的攻角、侧滑角辨识方法，所述方法包括：

s1、建立攻角α与侧滑角β的递推模型；

s2、建立攻角α与侧滑角β的观测模型；

s3、基于攻角α和侧滑角β各自的卡尔曼滤波初值和所述递推模型，得到攻角α和侧滑角β的一步预测值xpre；

s4、基于误差方差阵p和系统方差阵q，得到一步预测均方差ppre；

s5、基于所述观测模型，得到测量矩阵h和测量量z；

s6、基于一步预测均方差ppre、测量矩阵h和测量方差阵r，得到滤波增益k；

s7、基于一步预测值xpre、滤波增益k、测量矩阵h、测量量z、当前时刻的攻角α(k)和当前时刻的侧滑角β(k)，得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1)；

s8、基于滤波增益k、测量矩阵h和一步预测均方差ppre，得到更新后的误差方差阵p′；

s9、基于卡尔曼滤波后的攻角α(k+1)、卡尔曼滤波后的侧滑角β(k+1)和所述递推模型，得到攻角α和侧滑角β更新后的一步预测值xpre′，并返回s4基于更新后的误差方差阵p′和更新后的一步预测值xpre′进行下一轮迭代。

优选的，建立的攻角α与侧滑角β的递推模型为：

式中，α(k-1)为上一时刻的攻角，β(k-1)为上一时刻的侧滑角。

优选的，建立的攻角α与侧滑角β的观测模型为：

式中，为俯仰角，θ为弹道倾角，γ为滚转角，ny为纵向过载，nz为侧向过载，k1为第一飞行相关系数，k2为第二飞行相关系数。

优选的，通过下式得到一步预测均方差ppre：

ppre＝a·p·a^t+q，

式中，a为模型方阵。

优选的，通过下式得到滤波增益k：

k＝ppre·h^t(h·ppre·h^t+r)^-1，

其中，

优选的，通过下式得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1)：

其中，

优选的，通过下式得到更新后的误差方差阵p′：

p′＝(i-kh)ppre，

式中，i为单位矩阵。

应用本发明的技术方案，通过建立攻角α与侧滑角β的递推模型和观测模型，并基于攻角α与侧滑角β的递推模型和观测模型进行卡尔曼滤波，可实现攻角α和侧滑角β的高精度在线辨识。本发明在无攻角、侧滑角测量装置的条件下使用，节省了一套传感采集设备，而且可实现对攻角、侧滑角的同步测量。本发明可广泛运用于各种巡航类飞行器上，参与飞控算法解算，达到提升系统可靠性目的。

附图说明

所包括的附图用来提供对本发明实施例的进一步的理解，其构成了说明书的一部分，用于例示本发明的实施例，并与文字描述一起来阐释本发明的原理。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1示出了根据本发明一种实施例的基于卡尔曼滤波的攻角、侧滑角辨识方法的流程图；

图2示出了根据本发明一种实施例的攻角辨识仿真结果；

图3示出了根据本发明一种实施例的侧滑角辨识仿真结果。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本发明的范围。同时，应当明白，为了便于描述，附图中所示出的各个部分的尺寸并不是按照实际的比例关系绘制的。对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为授权说明书的一部分。在这里示出和讨论的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它示例可以具有不同的值。应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步讨论。

如图1所示，示出了一种基于卡尔曼滤波的攻角、侧滑角辨识方法，所述方法包括：

s1、建立攻角α与侧滑角β的递推模型；

s2、建立攻角α与侧滑角β的观测模型；

s3、基于攻角α和侧滑角β各自的卡尔曼滤波初值和所述递推模型，得到攻角α和侧滑角β的一步预测值xpre；

s4、基于误差方差阵p和系统方差阵q，得到一步预测均方差ppre；

s5、基于所述观测模型，得到测量矩阵h和测量量z；

s6、基于一步预测均方差ppre、测量矩阵h和测量方差阵r，得到滤波增益k；

s7、基于一步预测值xpre、滤波增益k、测量矩阵h、测量量z、当前时刻的攻角α(k)和当前时刻的侧滑角β(k)，得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1)；

s8、基于滤波增益k、测量矩阵h和一步预测均方差ppre，得到更新后的误差方差阵p′；

s9、基于卡尔曼滤波后的攻角α(k+1)、卡尔曼滤波后的侧滑角β(k+1)和所述递推模型，得到攻角α和侧滑角β更新后的一步预测值xpre′，并返回s4基于更新后的误差方差阵p′和更新后的一步预测值xpre′进行下一轮迭代。

本发明通过建立攻角α与侧滑角β的递推模型和观测模型，并基于攻角α与侧滑角β的递推模型和观测模型进行卡尔曼滤波，可实现攻角α和侧滑角β的高精度在线辨识。本发明在无攻角、侧滑角测量装置的条件下使用，节省了一套传感采集设备，而且可实现对攻角、侧滑角的同步测量。本发明可广泛运用于各种巡航类飞行器上，参与飞控算法解算，达到提升系统可靠性目的。

根据本发明的一种实施例，在飞行动力学中，可以默认在极短时间内，攻角α与侧滑角β不会发生变化，从而建立攻角α与侧滑角β的递推模型为：

式中，α(k)为当前时刻的攻角，β(k)为当前时刻的侧滑角，α(k-1)为上一时刻的攻角，β(k-1)为上一时刻的侧滑角。

在工程应用中，常使用纵向过载ny和侧向过载nz来估计攻角α和侧滑角β，以用于飞行控制及发动机控制中。目前使用的观测模型如下式所示：

式中，ny为纵向过载，nz为侧向过载，k1为第一飞行相关系数，k2为第二飞行相关系数。

其中，

式中，fy为纵向气动力，fz为侧向气动力，cy^α为纵向气动力系数，cy^β为侧向气动力系数，s为飞行器特征面积，q为动压。其中，纵向气动力系数cy^α和侧向气动力系数cy^β由风洞试验吹风获得，飞行器特征面积s是飞行器的固有状态变量，动压q通过下式获得：

q＝ρv²/2，

式中，ρ为空气密度，v为空速。其中，空气密度ρ依据标准大气并通过惯导测量的高度获得。由于很多飞行器未装备空速测量装置，因此在空速不可测的情况下，则可使用惯导测量的地速代替，从而获得近似动压q。

显然，受气动参数误差、实际空气密度与标准大气偏移、风场变化等因素影响，仅使用纵向过载ny和侧向过载nz来估计攻角α和侧滑角β存在较大误差。

在本发明中，攻角α和侧滑角β的观测模型增加了额外的约束，以提高攻角α和侧滑角β测量的精度和可靠性。

在本发明中，攻角α、侧滑角β与欧拉动力学关系如下式所示：

式中，为俯仰角，θ为弹道倾角，γ为滚转角，ψ为航向角，ψv为航迹角。其中，弹道倾角θ指空速弹道倾角，由于一般情况下，空速弹道倾角无法直接测量获得，因此，可采用地速弹道倾角近似代替。

在飞行过程中，由于攻角α、侧滑角β、航向角ψ与航迹角ψv之差一般为小量，上述的攻角α、侧滑角β与欧拉动力学关系可简化为下式：

设ψ-ψv为第一中间变量g，为第二中间变量l，则可进一步简化为下式：

将上式中的第一中间变量g消除，可得到下式：

对上式进行整理可得到攻角α、侧滑角β、俯仰角弹道倾角θ和滚转角γ之间的关系如下式所示：

在本发明中，对攻角α和侧滑角β的观测模型增加了俯仰角弹道倾角θ和滚转角γ的约束，从而提高攻角α和侧滑角β测量的精度和可靠性。建立的攻角α与侧滑角β的观测模型为：

其中，俯仰角滚转角γ、纵向过载ny和侧向过载nz均可由惯导直接获得，弹道倾角θ可由惯导测量的速度计算获得。

根据本发明的一种实施例，对于卡尔曼滤波，精确的初值可以减短收敛时间。因此，选取攻角α和侧滑角β各自的卡尔曼滤波初值为：

式中，α(0)为初始时刻的攻角，β(0)为初始时刻的侧滑角，为初始时刻的俯仰角，θ(0)为初始时刻的弹道倾角。

根据本发明的一种实施例，通过下式得到攻角α和侧滑角β的一步预测值xpre：

式中，αpre为攻角α的一步预测值，βpre为侧滑角β的一步预测值，α(k)为当前时刻的攻角，β(k)为当前时刻的侧滑角。

根据本发明的一种实施例，通过下式得到一步预测均方差ppre：

ppre＝a·p·a^t+q，

式中，a为模型方阵，a^t为模型方阵的转置。举例来说，因此

其中，误差方差阵p在滤波开始前根据系统特性直接赋值，然后在每一步迭代计算后进行更新。在本实施例中，系统方差阵q在滤波开始前根据系统特性直接赋值，但是在每一步迭代计算后不更新。在本实施例中，

根据本发明的一种实施例，通过下式得到滤波增益k：

k＝ppre·h^t(h·ppre·h^t+r)^-1，

其中，测量方差阵r在滤波开始前根据系统特性直接赋值，但是在每一步迭代计算后不更新。在本实施例中，

根据本发明的一种实施例，通过下式得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1)：

其中，

根据本发明的一种实施例，通过下式得到更新后的误差方差阵p′：

p′＝(i-kh)ppre，

其中，i为单位矩阵，

下面以某型号飞行器为例，对本发明提供的基于卡尔曼滤波的攻角、侧滑角辨识方法进行具体说明。

第一步，建立的攻角α与侧滑角β的递推模型为：

建立的攻角α与侧滑角β的观测模型为：

第二步，选取攻角α与侧滑角β的观测模型中的第一飞行相关系数k1和第二飞行相关系数k2。

由飞行动力学可知，第一飞行相关系数k1和第二飞行相关系数k2与气动系数、特征面积、动压、飞行器质量及当前飞行重力加速度相关。其中重力加速度、特征面积数值较为固定；动压与气动系数随飞行状态变化较大；飞行器质量是可估计的换变量。因此，若想提高辨识精度，可以依据气动数据、大气模型、耗油率，制定随飞行时间、高度、速度变化的差值表。以本实施例中的飞行器为例，设定k1＝1170，k2＝-0.1。

第三步，选取卡尔曼滤波的三个参数，误差方差阵p、系统方差阵q和测量方差阵r。

卡尔曼滤波参数关系卡尔曼滤波器性能。误差方差阵p、系统方差阵q和测量方差阵r的选取与飞行器性能等多因素相关。以本实施例中的飞行器为例，卡尔曼滤波的三个参数分别为：

第四步，基于攻角α与侧滑角β的递推模型和观测模型进行卡尔曼滤波，并得到攻角α和侧滑角β的仿真结果，如图2、图3所示。由图2、图3可知，若飞行器处于遇到弹道机动阶段，攻角α与侧滑角β的辨识精度稍差，攻角α误差为0.29°，侧滑角β误差为0.26°；若飞行器处于巡航平飞阶段，攻角α与侧滑角β的辨识精度较高，攻角α误差为0.0033°，侧滑角β误差为0.0077°，可以忽略不计；若飞行器处于倾斜转弯(btt)阶段，攻角α误差为0.026°，侧滑角β误差为0.02°。

综上所述，飞行器在稳定巡航飞行过程中，攻角α与侧滑角β的辨识误差均不超过0.1°，在弹道机动的过程中，攻角α与侧滑角β的辨识误差均不超过0.3°，满足测量的精度要求。

在本发明的描述中，需要理解的是，方位词如“前、后、上、下、左、右”、“横向、竖向、垂直、水平”和“顶、底”等所指示的方位或位置关系通常是基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，在未作相反说明的情况下，这些方位词并不指示和暗示所指的装置或元件必须具有特定的方位或者以特定的方位构造和操作，因此不能理解为对本发明保护范围的限制；方位词“内、外”是指相对于各部件本身的轮廓的内外。

为了便于描述，在这里可以使用空间相对术语，如“在……之上”、“在……上方”、“在……上表面”、“上面的”等，用来描述如在图中所示的一个器件或特征与其他器件或特征的空间位置关系。应当理解的是，空间相对术语旨在包含除了器件在图中所描述的方位之外的在使用或操作中的不同方位。例如，如果附图中的器件被倒置，则描述为“在其他器件或构造上方”或“在其他器件或构造之上”的器件之后将被定位为“在其他器件或构造下方”或“在其他器件或构造之下”。因而，示例性术语“在……上方”可以包括“在……上方”和“在……下方”两种方位。该器件也可以其他不同方式定位(旋转90度或处于其他方位)，并且对这里所使用的空间相对描述作出相应解释。

此外，需要说明的是，使用“第一”、“第二”等词语来限定零部件，仅仅是为了便于对相应零部件进行区别，如没有另行声明，上述词语并没有特殊含义，因此不能理解为对本发明保护范围的限制。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。