一种适用于直升机自适应飞行控制的鲁棒辨识方法与流程

本发明涉及直升机飞行动力学与飞行控制技术领域，具体是一种适用于直升机自适应飞行控制的鲁棒辨识方法。

背景技术：

直升机具有悬停、垂直起降和低速机动能力，成为不可或缺的重要飞行器。然而，直升机固有的强耦合和不稳定性也导致了直升机飞行品质较差，不易操纵。为了从根本上解决这个问题，最有效的方法，同时也是当前最常用的方法就是设计一套高质量的飞行控制系统。由于直升机高度的非线性和非定常特性，经典控制方法的控制效果往往不够理想，而基于现代控制理论的自适应控制方法成为解决这一难题的有效手段。

自适应控制方法的两个主要关键问题就是自适应控制策略设计和控制对象的实时辨识技术。就直升机飞行动力学模型的实时辨识技术而言，由于其固有特性非常复杂，振动水平和噪声水平高导致测量数据污染严重，给准确辨识其动力学模型带来了较大困难。早期的直升机飞行动力学模型辨识技术以采用单纯的参数辨识方法辨识解耦模型为主，然而解耦模型与直升机真实动力学特性差异较大，对控制器设计的效果产生较大影响。之后，全耦合模型的辨识技术得到较大发展，出现了结构辨识和参数辨识结合的综合辨识方法，然而，相对于解耦模型，全耦合模型的复杂特性也导致各种辨识方法精度差异明显。从飞行器模型辨识领域的权威专家美国NASA Ames研究中心飞行力学与飞行控制团队主任Mark B.Tischler在AIAA出版社出版的最新专著《Aircraft and Rotorcraft System Identification:Engineering Methods with Flight Test Examples 2^nd Edition》可以看出，目前国际上直升机领域飞行动力学模型的主流辨识方法仍旧是以经典辨识理论为主，对于测量噪声等外界干扰对数据的影响采用忽略或简单工程修正的方式加以处理。

目前，在已发展的直升机自适应控制方法中，直升机飞行动力学模型的实时辨识主要采用递推最小二乘法，同时将噪声视为白噪声，并通过对误差方差进行简单修正以补偿噪声引起的估计误差。然而，直升机实际噪声为复杂的有色噪声，同时，在非主要通道，高噪声水平导致了测量数据的信噪比很低。从而简单的工程修正方法并不能很好的解决复杂噪声对模型辨识精度产生较大影响的困难。综上，从目前国内外现有的直升机飞行动力学模型实时辨识技术来看，尚未有从根本上解决实际飞行过程中包括测量噪声、突风等在内的外界干扰对模型辨识精度影响的方法，进而给直升机自适应控制系统的设计造成了一定的困难。

技术实现要素：

本发明为了解决当前基于递推最小二乘法的在线辨识方法鲁棒性不足的问题，提供了一种适用于直升机自适应飞行控制的鲁棒辨识方法，能够有效消除测量噪声等外界扰动对模型辨识精度的影响、显著提高动力学模型辨识的鲁棒性。

本发明包括以下步骤：

第一步，初始化直升机悬停状态的状态空间模型，如式(1)所示，计算出悬停状态下直升机状态空间模型中各个参数的取值，

式中，x为直升机状态向量，A为稳定性矩阵，B为操纵矩阵，y为直升机观测向量，C和D分别为观测矩阵以及操纵对观测向量的影响矩阵；

第二步，从接通飞控后开始，以第一步得到的参数值作为模型参数初值，基于四阶龙格-库塔法对式(1)所示的直升机飞行动力学状态空间方程进行单步求解，即得到t时刻的状态响应和观测向量x(t)以及y(t)，同时各个机载传感器记录下直升机实际观测向量y_m(t)，计算模型估计误差；

第三步，将实际观测向量通过式(3)所示的滤波器进行低通滤波处理，并记录滤波后的实际观测向量y_mf(t)，根据式(4)和式(5)分别计算N个观测量在时间窗宽度为L的噪声序列的均值和方差，并以此方差作为当前时刻N个观测量的噪声边界σ(t)，具体表达式如式(6)所示。与此同时，在直直升机悬停状态的状态空间模型的基础上建立如式(7)所示的模型参数灵敏度方程，并同样基于四阶龙格-库塔法求解直升机状态响应向量对各个模型参数的灵敏度

σ(t)＝[Var(n₁),Var(n₂),…,Var(n_N)]^T (6)，

其中，T为一阶滤波函数的时间常数，s为复频变量，为第i个观测量的噪声均值，Var(n_i)为第i个观测量的噪声方差。

第四步，基于最优边界椭球法对模型参数偏差进行估计，首先，定义模型参数偏差向量为Δθ，并设置其初值为0I，其中，I为单位矩阵，分别根据式(8)和式(9)计算中间输出向量和中间回归矩阵

其次，根据式(10)—(12)计算用于优化边界椭球的三个系数C₁、C₂和C₃，再利用式(13)计算得到最优边界椭球加权系数λ(t)，

C₀＝m[tr[σ(t)σ(t)^T]-tr[e(t)e(t)^T]]-κ(t-1)x(t)^TP(t-1)x(t) (12)，

其中，tr表示矩阵的迹，m为回归矩阵的行数，e(t)＝y(t)-y_mf(t)为模型计算的响应与实际测量的响应之间的误差，P为协方差矩阵，κ为椭球形状参数，协方差矩阵P和椭球形状参数κ在第一次计算中取初值为P(0)＝10-⁶I和κ(0)＝1，I为单位矩阵；

第五步，建立式(14)—式(16)所示的迭代算法，并利用第四步中得到的最优边界椭球加权系数修正参数的辨识结果，得到待辨识参数的最优可行解集，并以椭球中心点作为模型参数估计偏差Δθ的辨识结果，

Δθ(t)＝Δθ(t-1)+λ(t)P(t)x(t)e(t)^T (14)，

第六步，根据当前模型参数估计偏差的辨识结果利用式(17)修正模型参数并更新式(1)中的A和B矩阵，使用更新后的直升机悬停状态的状态空间模型回到第二步继续进行辨识计算。

θ_i＝θ_i-1+Δθ (17)

其中，θ_i为更新后的模型参数向量，θ_i-1为更新前的上一迭代过程的模型参数向量。

第一步中，基于非线性机理建模，然后通过配平和线化计算出式(1)所示的模型中各个参数的取值。

第一步中，通过一次单独的悬停状态扫频飞行试验，利用最小二乘法计算出式(1)中各个参数的取值。

所述的最小二乘法具体步骤如下，首先，试飞员将直升机在悬停状态进行配平，然后依次施加总距、纵向周期变距、横向周期变距和尾桨总距的正弦扫频激励信号，每组操纵持续10—30秒，将得到的飞行试验数据进行滤波、传感器位置校正和数据相容性检查，之后可根据式(2)所示的标准最小二乘法估计出式(1)中各个参数的取值，

θ＝(X^TX)^-1X^TY (2)，

其中，θ为待辨识参数，由式(1)中矩阵A和矩阵B中的参数组成，Y是输出向量，由式(1)中的y向量组成，X是回归矩阵，由式(1)中的x向量组成。

本发明有益效果在于：

1)以参数估计偏差向量作为中间辨识参数解决了边界椭球法无法用于微分系统的难题，实现了边界椭球法在直升机飞行动力学模型辨识中的应用。

2)通过实时估计噪声边界有效提高边界椭球的优化效果，增强了参数辨识精度。

3)现有递推最小二乘方法需要知道噪声的具体统计特性，由于实际噪声的统计学特性无法获取导致方法精度尤其是鲁棒性的降低。而本发明以边界椭球参数加权的参数递推辨识算法仅需要考虑噪声边界，所以适用于各种非理想噪声干扰环境下的模型辨识，相比现有方法能够提高参数辨识的精度，而鲁棒性则更得到显著改善。

附图说明

图1是本发明的鲁棒辨识方法实施流程。

图2是直升机飞行动力学模型初始化流程。

图3是实施例中本发明方法辨识结果的验证。

图4是实施例中高信噪比条件下本发明方法与递推最小二乘法的辨识效果对比。

图5是实施例中低信噪比条件下本发明方法与递推最小二乘法的辨识效果对比。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步说明。

本发明提供了一种适用于直升机自适应飞行控制的鲁棒辨识方法，流程如图1所示，具体包括以下几个步骤：

第一步，如图2所示，针对对象直升机进行飞行动力学模型初始化，即初始化悬停状态的状态空间模型，如式(1)所示。其中，x为直升机状态向量，包括体轴系下的速度、角速度和姿态角，A为稳定性矩阵，B为操纵矩阵，y为直升机观测向量，如空速、姿态角、角速度、加速度等等，C和D分别为观测矩阵以及操纵对观测向量的影响矩阵。一般情况，D为零矩阵，C由所选观测量和直升机状态变量的转换关系直接确定。需要进行初始化工作的主要就是A和B两个矩阵。由于对于直升机，其每次飞行总是从垂直起飞到悬停开始，所以对飞行动力学模型的初始化只需要针对悬停状态进行即可。此外，对于同一直升机，尽管可能每次飞行的起飞重量有所区别，但在同一飞行状态下对于A和B矩阵中的无量纲参数影响不大，所以只需进行1次初始化计算即可。当初始化完成后，将得到悬停状态模型参数记录到模型数据库中，之后的每次飞行可直接导入模型数据即可完成初始化工作，无需重新计算。

初始化悬停状态的直升机状态空间方程有两种途径，一种是基于非线性机理建模，然后通过配平和线化得到式(1)所示的模型。另外一种途径也是本发明采用和推荐的方法，即通过一次单独的悬停状态扫频飞行试验，利用最小二乘法得到式(1)中A和B矩阵中的模型参数，具体流程如图2所示。首先，试飞员将直升机在悬停状态进行配平，然后依次施加总距、纵向周期变距、横向周期变距和尾桨总距的正弦扫频激励信号，每组操纵持续10—30秒。将得到的飞行试验数据进行滤波、传感器位置校正和数据相容性检查，确保数据可信。之后可根据如式(2)所示的标准最小二乘法估计出模型参数。其中，θ为待辨识参数，由矩阵A和B中的参数组成，Y是输出向量，由式(1)中的y向量组成，X是回归矩阵，由式(1)中的x向量组成。

θ＝(X^TX)^-1X^TY (2)。

第二步，从接通飞控后开始，进入模型实时辨识环节，即本发明的主体部分。首先基于四阶龙格-库塔法对式(1)所示的直升机飞行动力学状态空间方程进行单步求解，即得到t时刻的状态响应和观测向量x(t)以及y(t)，同时各个机载传感器记录下直升机实际观测向量y_m(t)。

第三步，将实际观测向量通过式(3)所示的滤波器进行低通滤波处理，并记录滤波后的实际观测向量y_mf(t)。根据式(4)和式(5)分别计算N个观测量在时间窗宽度为L的噪声序列的均值和方差，并以此方差作为当前时刻N个观测量的噪声边界σ(t)，具体表达式如式(6)所示。与此同时，在直升机状态空间模型(1)的基础上建立如式(7)所示的模型参数灵敏度方程，并同样基于四阶龙格-库塔法求解直升机状态响应向量对各个模型参数的灵敏度

σ(t)＝[Var(n₁),Var(n₂),…,Var(n_N)]^T (6)，

其中，T为一阶滤波函数的时间常数，s为复频变量，为第i个观测量的噪声均值，Var(n_i)为第i个观测量的噪声方差。

第四步，基于最优边界椭球法对模型参数偏差进行估计。首先，定义模型参数偏差向量为Δθ，并设置其初值为0I，其中，I为单位矩阵。分别根据式(8)和式(9)计算中间输出向量和中间回归矩阵

其次，根据式(10)—(12)计算用于优化边界椭球的三个系数。其中，tr表示矩阵的迹，m为回归矩阵的行数，e(t)＝y(t)-y_mf(t)为模型计算的响应与实际测量的响应之间的误差，协方差矩阵P和椭球形状参数κ在第一次计算中可取初值为P(0)＝10^-6I和κ(0)＝1，I为单位矩阵,此后这两个参数都通过迭代计算得到，具体算法在本发明的第五步中介绍。利用式(13)计算得到最优边界椭球加权系数λ(t)。

C₀＝m[tr[σ(t)σ(t)^T]-tr[e(t)e(t)^T]]-κ(t-1)x(t)^TP(t-1)x(t) (12)，

第五步，本发明建立式(14)—式(16)所示的迭代算法，并利用第四步中得到的最优边界椭球加权系数修正参数的辨识结果，得到待辨识参数的最优可行解集，并以椭球中心点作为模型参数估计偏差Δθ的辨识结果。

Δθ(t)＝Δθ(t-1)+λ(t)P(t)x(t)e(t)^T (14)，

第六步，根据当前模型参数估计偏差的辨识结果利用式(17)修正模型参数并更新式(1)中的A和B矩阵。由于本发明的实时辨识算法主要面向飞行控制系统使用，所以只要飞行控制系统没有断开，则该辨识流程持续执行，即完成第六步后，使用更新后的状态空间方程(1)，并回到第二步继续进行辨识计算。

θ_i＝θ_i-1+Δθ (17)

其中，θ_i为更新后的模型参数向量，θ_i-1为更新前的上一迭代过程的模型参数向量。

本发明的一种具体实施例如下：

在本实施例中，将本发明的直升机飞行动力学模型鲁棒辨识方法用于实际直升机的在线实时辨识，采用的是某2吨级的轻型直升机作为辨识对象，进行了大约30秒左右的飞行试验，并基于实时飞行测试数据完成了模型辨识，图3展示了本发明方法的辨识效果，可以看出，辨识精度是比较高的。为了显示本发明的优越性，同时采用了递推最小二乘法也进行了实时辨识，并与本发明方法的辨识效果进行了对比。从图4和图5可以看出，在信噪比较高的条件下(图4)，本发明方法比递推最小二乘法精度略高，但在噪声水平提高，信噪比较低的情况下(图5)，递推最小二乘法的精度显著下降，而本发明方法的精度仍然保持较好，具有较高的鲁棒性。

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。