直升机相对于固定翼飞机具有垂直升降,全地形作战，近地飞行的独特优势，它的这些特点使之成为一种“万用”的交通工具^[1-2]。但是，由于旋翼气动特性与机身结构的复杂性，直升机又具有振动大、噪声水平高等缺点，这无疑对直升机飞行动力学的建模带来了很大的困难。

系统辨识的建模方法为提高直升机飞行动力学模型的精度提供了有效的技术手段^[3]，这种方法是根据飞行试验得到的输入输出数据建立系统的数学模型，其特点是不需要知道对象详细的物理细节。从20世纪50年代开始，飞行器飞行动力学就开始使用系统辨识的方法进行数学建模，然而，直到20世纪60年代末70年代初，系统辨识方法才开始在直升机飞行动力学建模中得到应用^[4]。文献^[5]将成功应用在固定翼飞机上的最小二乘法、输出误差法和极大似然法进行改造，并对MK-50海王直升机展开了相关的辨识工作，研究结果表明，经过适当改造的算法可适用于直升机飞行动力学的辨识。文献^[6]中提出了一种高效批量的子空间辨识算法并对V-22直升机的结构进行了辨识。文献^[7]中提出了一种基于径向基函数网络的参数估计方法，这种方法不需要考虑直升机的数学模型，旋翼的参数可以直接从飞行试验数据得到。文献^[8]中将两步方程误差/输出误差状态估计算法应用在了微型电动直升机上，得到了直升机在盘旋飞行时的线性模型。

中国国内直升机飞行动力学模型辨识技术的发展相对缓慢。文献^[9]中研究了求解直升机气动导数的拉格朗日法，通过引入拉格朗日乘子，减小了直升机非线性对辨识结果的影响，提高了辨识结果的精度，并将研究成果直接应用到了直6、直8、直9及直11等多种型号直升机的气动参数辨识当中。文献^[10]应用最小二乘一次完成算法与广义最小二乘法对直升机的自转着陆模型的气动参数辨识。文献^[11]中提出了一种渐消记忆的最小二乘逐状态辨识算法，对小型无人直升机悬停状态进行了辨识，该算法可以减少计算量并提高计算过程的稳定性。文献^[12]中提出了一种确定直升机垂直状态气动参数的辨识方法并且利用辨识结果确定了直升机的悬停升限。虽然系统辨识有了很大的发展，有如最小二乘法、输出误差法这些经典的辨识理论，也有如子空间辨识、集员辨识这些现代辨识理论，但是这些辨识算法都有各自的缺点。经典的辨识方法往往假定外在干扰是白噪声，现代辨识理论较为复杂，许多地方还不成熟，有待进一步的研究发展。直升机是一个有着非理想噪声干扰的对象，这就要求发展一套考虑非理想噪声干扰、简单易用、能进行在线辨识的算法。本文则针对有非理想噪声干扰的直升机飞行动力学模型问题，提出了这样一种辨识方法。

1 辨识模型

由于本文主要是研究非理想噪声对辨识结果的影响，为了消除耦合性带来的影响，选取UH-60直升机悬停状态下的纵向运动方程并线性化为状态空间的形式

$\dot{x}=Ax+B{{u}_{in}}$

(1)

式中：x=[Δu,Δw,Δq,Δθ]^T,为直升机的状态量的增量，为了表示方便，下文去掉增量符号Δ； u_in=Δδ_e ，表示直升机的纵向周期变距增量，为了表示方便，下文去掉增量符号Δ；A为气动导数矩阵；B为操纵导数矩阵。

$A=\left[ \begin{matrix} {{X}_{u}} & {{X}_{w}} & {{X}_{q}} & -gcos{{\theta }_{0}} \\ {{Z}_{u}} & {{Z}_{w}} & {{Z}_{q}} & -gsin{{\theta }_{0}}cos{{\theta }_{0}} \\ {{M}_{u}} & {{M}_{w}} & {{M}_{q}} & 0 \\ 0 & 0 & 1 & 0 \\ \end{matrix} \right]$

(2)

$B=\left[ \begin{matrix} {{X}_{\delta e}} \\ {{Z}_{\delta e}} \\ {{M}_{\delta e}} \\ 0 \\ \end{matrix} \right]$

(3)

式中：θ₀为配平后俯仰角的值；X表示体轴系下X轴气动力；Z表示体轴系下Z轴气动力；M表示绕直升机的俯仰力矩。

2 辨识算法 2.1 最小二乘法

假设系统的输出输入关系可以描述成如下的最小二乘格式^[13]

$z\left( k \right)={{h}^{T}}\left( k \right)\theta \left( k \right)+n\left( k \right)$

(4)

式中n(k)为白噪声。

$\left\{ \begin{matrix} h\left( k \right)=[-z\left( k-1 \right),z\left( k-2 \right),\ldots ,-z\left( k-{{n}_{a}} \right), \\ u\left( k-1 \right),u\left( k-2 \right),\ldots ,u(k-{{n}_{b}}){{]}^{T}} \\ \theta \left( k \right)={{[{{a}_{1}},{{a}_{2}},\ldots ,{{a}_{{{n}_{a}}}},{{b}_{1}},{{b}_{2}},\ldots ,{{b}_{n}}_{b}]}^{T}} \\ \end{matrix} \right.$

(5)

为了使原模型能够使用最小二乘辨识算法，将式(2,3) 代入式(1) ，可得

$\left[ \begin{matrix} {\dot{u}} \\ {\dot{w}} \\ {\dot{q}} \\ {\dot{\theta }} \\ \end{matrix} \right]=\left[ \begin{matrix} {{X}_{u}} & {{X}_{w}} & {{X}_{q}} & -gcos{{\theta }_{0}} \\ {{Z}_{u}} & {{Z}_{w}} & {{Z}_{q}} & -gsin{{\theta }_{0}}cos{{\theta }_{0}} \\ {{M}_{u}} & {{M}_{w}} & {{M}_{q}} & 0 \\ 0 & 0 & 1 & 0 \\ \end{matrix} \right]\left[ \begin{matrix} u \\ w \\ q \\ \theta \\ \end{matrix} \right]+\left[ \begin{matrix} {{X}_{{{\delta }_{e}}}} \\ {{Z}_{{{\delta }_{e}}}} \\ {{M}_{{{\delta }_{e}}}} \\ 0 \\ \end{matrix} \right]{{\delta }_{e}}$

(6)

整理成最小二乘格式可得

$\left[ \begin{matrix} {{a}_{x}} \\ {{a}_{y}} \\ {\dot{q}} \\ \end{matrix} \right]=\left[ \begin{matrix} u & w & q & {{\delta }_{e}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & u & w & q & {{\delta }_{e}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u & w & q & {{\delta }_{e}} \\ \end{matrix} \right]{{\theta }_{1}}+\left[ \begin{matrix} {{v}_{1}} \\ {{v}_{2}} \\ {{v}_{3}} \\ \end{matrix} \right]$

(7)

式中：v₁，v₂，v₃为白噪声。

$\left\{ \begin{matrix} {{a}_{x}}=\dot{u}+gcos{{\theta }_{0}}\cdot \theta \\ {{a}_{y}}=\dot{w}+gsin{{\theta }_{0}}cos{{\theta }_{0}}\cdot \theta \\ {{\theta }_{1}}={{[{{X}_{u}},{{X}_{w}},{{X}_{q}},{{X}_{{{\delta }_{e}}}},{{Z}_{u}},{{Z}_{w}},{{Z}_{q}},{{Z}_{{{\delta }_{e}}}},{{M}_{u}},{{M}_{w}},{{M}_{q}},{{M}_{{{\delta }_{e}}}}]}^{T}} \\ \end{matrix} \right.$

(8)

将式(8) 与式(4) 进行比较，令

$\begin{align} & ~z\left( k \right)=\left[ \begin{matrix} {{a}_{x}}\left( k \right) \\ {{a}_{y}}\left( k \right) \\ \dot{q}\left( k \right) \\ \end{matrix} \right]{{h}^{T}}\left( k \right)= \\ & \left[ \begin{matrix} u\left( k \right) & w\left( k \right) & q\left( k \right) & {{\delta }_{e}}\left( k \right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & u\left( k \right) & w\left( k \right) & q\left( k \right) & {{\delta }_{e}}\left( k \right) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u\left( k \right) & w\left( k \right) & q\left( k \right) & {{\delta }_{e}}\left( k \right) \\ \end{matrix} \right] \\ \end{align}$

(9)

若k=1,2,…,L(L为数据长度)，则式(7) 的最小二乘估计为

${{\theta }_{LS}}={{({{H}^{T}}_{L}{{H}_{L}})}^{-1}}{{H}^{T}}_{L}{{z}_{L}}$

(10)

式中

$\begin{align} & {{z}_{L}}={{\left[ {{a}_{x}}\left( 1 \right),{{a}_{y}}\left( 1 \right),\dot{q}\left( 1 \right),\ldots ,{{a}_{x}}\left( L \right),{{a}_{y}}\left( L \right),\dot{q}\left( L \right) \right]}^{T}}{{H}_{L}}=\left[ \begin{matrix} {{h}^{T(1)}} \\ \vdots \\ {{h}^{T}}\left( L \right) \\ \end{matrix} \right]= \\ & \left[ \begin{matrix} u\left( 1 \right) & w\left( 1 \right) & q\left( 1 \right) & {{\delta }_{e}}\left( 1 \right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & u\left( 1 \right) & w\left( 1 \right) & q\left( 1 \right) & {{\delta }_{e}}\left( 1 \right) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u\left( 1 \right) & w\left( 1 \right) & q\left( 1 \right) & {{\delta }_{e}}\left( 1 \right) \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ u\left( L \right) & w\left( L \right) & q\left( L \right) & {{\delta }_{e}}\left( L \right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & u\left( L \right) & w\left( L \right) & q\left( L \right) & {{\delta }_{e}}\left( L \right) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u\left( L \right) & w\left( L \right) & q\left( L \right) & {{\delta }_{e}}\left( L \right) \\ \end{matrix} \right] \\ \end{align}$

(11)

2.2 本文算法

式(10) 是最小二乘的一次完成算法，虽然此算法对白噪声的模型辨识有一定的精度，且在理论研究中也有许多方便之处，但是计算方面需要计算矩阵的逆，特别是当H_L的维数比较大时，计算量比较大。而且本文的算法要求是能够进行在线辨识，显然一次完成算法是不符合要求的。直升机在飞行的过程当中会受到非理想噪声的干扰，所以用普通的最小二乘算法不能准确地估计出直升机飞行动力学的模型。本文基于增广最小二乘递推算法设计了一套同时考虑模型参数和噪声参数的综合辨识方法，将模型噪声作为辨识的参数加入到θ中去。

设非理想噪声的模型为

$e\left( k \right)=\left[ \begin{matrix} {{e}_{1}}\left( k \right) \\ {{e}_{2}}\left( k \right) \\ {{e}_{3}}\left( k \right) \\ \end{matrix} \right]=\left[ \begin{matrix} {{v}_{1}}\left( k \right)+{{d}_{1}}{{v}_{1}}\left( k-1 \right)+{{d}_{2}}{{v}_{1}}\left( k-2 \right) \\ {{v}_{2}}\left( k \right)+{{d}_{3}}{{v}_{2}}\left( k-1 \right)+{{d}_{4}}{{v}_{2}}\left( k-2 \right) \\ {{v}_{3}}\left( k \right)+{{d}_{5}}{{v}_{3}}\left( k-1 \right)+{{d}_{6}}{{v}_{3}}\left( k-2 \right) \\ \end{matrix} \right]$

(12)

将此噪声加入到模型中去，得到非理想噪声干扰的直升机飞行动力学纵向模型为

$\left[ \begin{matrix} {{a}_{x}} \\ {{a}_{y}} \\ {\dot{q}} \\ \end{matrix} \right]=\left[ \begin{matrix} u & w & q & {{\delta }_{e}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & u & w & q & {{\delta }_{e}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u & w & q & {{\delta }_{e}} \\ \end{matrix} \right]{{\theta }_{1}}+\left[ \begin{matrix} {{e}_{1}} \\ {{e}_{2}} \\ {{e}_{3}} \\ \end{matrix} \right]$

(13)

显然式(13) 是用普通最小二乘估计模型参数，将式(13) 变形写成增广最小二乘的格式为

$\begin{align} & \left[ \begin{matrix} {{a}_{x}}\left( k \right) \\ {{a}_{y}}\left( k \right) \\ \dot{q}\left( k \right) \\ \end{matrix} \right]= \\ & \left[ \begin{matrix} u\left( k \right) & w\left( k \right) & q\left( k \right) & {{\delta }_{e}}\left( k \right) & {{v}_{1}}\left( k-1 \right) & {{v}_{1}}\left( k-2 \right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & u\left( k \right) & w\left( k \right) & q\left( k \right) & {{\delta }_{e}}\left( k \right) & {{v}_{1}}\left( k-1 \right) & {{v}_{1}}\left( k-2 \right) & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u\left( k \right) & w\left( k \right) & q\left( k \right) & {{\delta }_{e}}\left( k \right) & {{v}_{3}}\left( k-1 \right) & {{v}_{3}}\left( k-2 \right) \\ \end{matrix} \right] \\ & \theta +\left[ \begin{matrix} {{v}_{1}}\left( k \right) \\ {{v}_{2}}\left( k \right) \\ {{v}_{3}}\left( k \right) \\ \end{matrix} \right] \\ \end{align}$

(14)

式中：θ=[X_u,X_w,X_q,X_δe,d₁,d₂,Z_u,Z_w,Z_q,Z_δe,d₃,d₄,M_u,M_w,M_q,M_δe,d₅,d₆]^T。

为了减小计算量且能够进行在线辨识，一次完成算法显然不适用，这时必须使用递推算法

$\left\{ \begin{matrix} K\left( k \right)=P\left( k-1 \right)h\left( k \right){{[{{h}^{T}}\left( k \right)P\left( k-1 \right)h\left( k \right)+I]}^{-1}} \\ P\left( k \right)=[I-K\left( k \right){{h}^{T}}\left( k \right)]P\left( k-1 \right) \\ \hat{\theta }\left( k \right)=\left( k-1 \right)+K\left( k \right)[z\left( k \right)-{{h}^{T}}\left( k \right)\hat{\theta }\left( k-1 \right)] \\ \end{matrix} \right.$

(15)

式中：K(k)为增益阵；P(k)为数据协方差阵；I为3维的单位矩阵。

3 模型的仿真验证与分析 3.1 仿真与辨识

本文的输入输出数据都是通过仿真来获得的，在MATLAB/SIMULINK平台建立如图 1所示的模型。

由于本文中辨识的模型为直升机的纵向通道模型，在获取仿真实验数据时，输入的信号为推杆扫频信号，最大幅值为1英寸。驾驶员在操纵驾驶杆时的频率不可能太大，所以扫频信号的频率范围为0.1~5 Hz，如图 2所示。

为了能够辨识噪声模型中的参数，本文中仅考虑加速度和角加速度被非理想噪声干扰，而速度和角速度的信号不被噪声干扰的情况。同时，为了使仿真实验的噪声数据接近真实的飞行试验，在参考国外的一些研究报告的基础上，确定了一般信噪比在50 dB左右，噪声较大的信噪比在25 dB以下。在推杆操纵的情况下，传感器所测得的纵向参数的灵敏度高，垂向参数的灵敏度低，所以仿真实验中，纵向输出数据的信噪比取为40 dB，垂向输出数据的信噪比为25 dB，则仿真输出信号如图 3,4所示。

从图 3，4可以看出，两条曲线的信噪比并不相同，主要是因为在驾驶员推杆的情况下，直升机在各通道的灵敏度不同，纵向通道灵敏度高，信噪比大，垂向通道灵敏度低，信噪比小，这和实际情况吻合。

根据式(10) ，使用普通最小二乘算法对直升机纵向通道模型参数进行辨识。

根据式(14,15) ,使用本文方法对直升机纵向通道的参数进行辨识。在进行递推之前，要选择初始的辨识参数${\hat{\theta }}$(0) 和P(0) ,一般情况下，取P(0) =a²I ，其中a为充分大的实数；${\hat{\theta }}$(0) =ε，ε为充分小的实向量。

由于在h(k)中含有不可测量的噪声${\hat{v}}$(k-1) ,…,${\hat{v}}$(k-n_d)，要用相应的估计值去替代真实的噪声。这时可以置

$\hat{v}\left( k \right)=z\left( k \right)-{{h}^{T}}\left( k \right)\hat{\theta }\left( k-1 \right)$

(16)

根据以上步骤，应用本文算法，可以辨识出气动导数矩阵A和控制矩阵B，从而得到UH-60直升机的纵向动力学模型为

$\begin{align} & {{\left[ \begin{matrix} {\dot{u}} & {\dot{w}} & {\dot{q}} & {\dot{\theta }} \\ \end{matrix} \right]}^{T}}=\left[ \begin{matrix} -0.0233 & 0.0232 & 2.8069 & -32.0272 \\ 0.0228 & -0.2935 & 0.3600 & -2.8283 \\ 0.0036 & 0.0018 & -0.8163 & 0 \\ 0 & 0 & 1 & 0 \\ \end{matrix}\text{ } \right] \\ & \left[ \begin{matrix} u \\ w \\ q \\ \theta \\ \end{matrix} \right]+\left[ \begin{matrix} -1.6584 \\ -0.1370 \\ 0.3347 \\ 0 \\ \end{matrix} \right]{{\delta }_{e}} \\ \end{align}$

(17)

3.2 模型验证

为了验证辨识参数的准确性，本文使用了3-2-1-1信号作为输入信号，输出拟合结果表示在图 5～7中。

为了验证本文方法的优越性，将本文方法所获得的辨识结果与普通最小二乘法的结果作对比。表 1中直接比较了辨识参数与真值的接近程度。表 2给出了3个输出的TIC(Theil inequality coefficient)值的比较。图 8，9为输出值的拟合情况对比。图 10为本文方法所估计的非理想噪声与仿真噪声的误差。

TIC值反映了辨识参数的预测能力。其值越小，预测能力越好。定义如下

$TIC=\frac{\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}{{{({{x}_{i}}-{{{\tilde{x}}}_{i}})}^{2}}}}}{\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}{{{x}^{2}}_{i}}}+\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}{{{{\tilde{x}}}^{2}}_{i}}}}$

(18)

式中:x_i表示仿真数据；x_i表示飞行试验数据，在本文中即为噪声数据；n表示数据的长度。

由表 1可知，本文方法的精度高于普通最小二乘法。特别是X_w,Z_w,M_w的值，其误差较大，甚至出现了反号的情况。这主要是因为给定的输入为推杆操纵，此时传感器测得的垂直方向的噪声较大，辨识的误差就会比较大。对于噪声模型的辨识参数，很明显一阶参数辨识精度较高，二阶参数的效果则不是很理想，随着阶数的增大，辨识精度会随之降低。表 2所示的是各个输出的TIC值，很显然，应用本文方法计算的值比普通最小二乘法的值要小，说明其预测能力更好。

图 8，9反映的是应用最小二乘方法与本文方法得到的计算值与仿真输出拟合对比，由于最小二乘法不能估计噪声而本文方法辨识出了噪声模型，且速度和角速度没有被噪声干扰，所以本文方法计算曲线与仿真曲线拟合的程度更好，最小二乘法的计算曲线则更光滑。图 10是估计的非理想噪声与仿真噪声的误差值，可以看出，在第1 s内，由于没有输入，噪声基本为零。在1 s后，随着递推步数的增加，其误差值越来越小，这也说明了利用本文方法估计噪声模型的可行性和有效性。

4 结束语

本文建立了UH-60直升机非理想噪声干扰的纵向飞行动力学模型，提出了一套适用于非理想噪声干扰下的直升机动力学模型的辨识方法。通过本文方法与最小二乘方法的对比,验证了本文方法建模的精确性以及相对于最小二乘方法的优越性。