飞机飞行品质对飞行安全至关重要，目前低阶等效系统方法是评估电传飞机飞行品质的重要手段^[1-4]，但随着飞机进一步采用放宽稳定性设计和高增益飞控设计，飞机动响应阻尼比越来越大(纵向模态达到0.6~0.8，横航向超过0.5) 。传统的试验方法主要有两类：短时激励和长时激励。短时激励方法如脉冲/倍脉冲、阶跃/倍阶跃等方法，不容易充分激励飞机动态特性，而且飞机响应衰减非常快；而长时激励如“3-2-1-1”、扫频等方法虽能充分激励飞机动态特性，但是动作时间长，大气扰流激励等非飞行员输入激励累计误差大，加上测试引入的误差导致现有成熟辨识方法很难精确有效辨识低阶等效系统参数。辨识方法包括时域法和频域法两类，时域法发展相对成熟，试飞工程中一般采用最大似然辨识法，方法成熟可靠^[5-6]，但对激励/响应数据相关性要求高(试飞数据质量要求高)，对初值选取也有较高要求，特别是系统时间延迟直接参与辨识时，计算速度慢，不适合实时/准实时在线分析；频域法主要为基于最小二乘法的各类优化辨识方法^[7-15]，辨识收敛速度快，在工程中广泛应用，但对初值选取要求高，而且噪声干扰对辨识结果影响很大，同一试验点、不同激励方法，可能得到的辨识结果差别很大，有时甚至连品质等级评价都不一致。为此，有研究人员对最小二乘法进行了一些改进^[16]，也有学者采用子空间方法进行参数辨识^[17]。

本文提出了采用Chirp-Z法进行高精度时频变换，利用频域方程误差法(Equation error method, EEM)选取辨识参数初值，再利用输出误差法(Output error method, OEM)进行参数精确辨识的方法，无论何种试飞方法，只要能将飞行模态响应激励出来，均可精确辨识出系统参数，而且结果收敛性好(试验结果重复性好)，有效地解决了新一代电传飞机飞行品质评估技术难题。此外，本文还探讨了针对现代飞机开展低阶等效系统评价的试验方法及数据要求，保证了该方法的工程实用性。

1 基本原理和方法

频域法是低阶等效系统参数辨识的最基本的方法，其应用最为广泛，同时也是国军标GJB-2874、美军标MIL-STD-1797A所推荐的方法，这种方法反映了系统在主要频段内的基本特性。

1.1 高精度时频变换法(Chrip-Z变换)

离散傅里叶变换可定义为^[4]

${{X}_{k}}=\sum\limits_{i=0}^{N-1}{{{X}_{i}}}{{e}^{-j2\pi {{f}_{k}}{{t}_{i}}}}~k=0,1,\text{ }2,\ldots ,\text{ }M-1$

(1)

式中

${{f}_{k}}=\frac{k}{N\Delta t}=\frac{k}{T}k=0,\text{ }1,\text{ }2,\ldots ,\text{ }M-1$

(2)

$M=\left\{ \begin{matrix} N/2 & N为偶数 \\ \left( N+1 \right)\text{ }/\text{ }2 & N为奇数 \\ \end{matrix} \right.$

(3)

现定义

$\begin{align} & \theta ={{\omega }_{k}}\Delta t=2\pi {{f}_{k}}\Delta t=\frac{2\pi k}{N} \\ & \begin{array}{*{35}{l}} k=0,1,\text{ }2,\ldots ,\text{ }M-1 \\ {} \\ \end{array} \\ \end{align}$

(4)

可将式(1) 写为

${{X}_{k}}=\sum\limits_{i=0}^{N-1}{{{X}_{i}}}{{e}^{-j\theta i}}k=0,\text{ }1,\text{ }2,\ldots ,\text{ }M-1$

(5)

$X\left( \theta \right)=X(2\pi {{f}_{k}}\Delta t)\approx \Delta t{{X}_{k}}$

(6)

离散傅里叶变换的缺点是其频率间隔依赖于数据的长度T，对于较短数据，采用这种方法所得到的频率分辨率很低，这将导致在频域中重要细节的丢失，影响后续数据分析和建模的结果，在极端情况下，离散频率间隔过大将导致在频域中重要特征丢失。此外，在离散傅里叶变换中，频率点是从0至奈奎斯特频率f_n之间平均分布的(采样率一般为32，频率分布0~100 rad/s)，但飞行品质主要关心其中的一小部分频带(0.1~ 10 rad/s)，这意味着大多频率点在感兴趣频带之外，导致了可利用的数据点非常少，数据特征丢失，同样也降低了数据分析和建模的精度。而Chirp-Z变换方法可以在任意频带范围内进行频谱分析，从而改善了数据长度较短时导致的频率分辨率低的问题。

Chirp-Z变换假定从频带f₀ , f₁中选取M个离散频率

${{f}_{k}}={{f}_{0}}+k\Delta f\text{ }\Delta f=\frac{({{f}_{1}}-{{f}_{0}})}{M}$

(7)

因此式(1) 可以写为

${{X}_{k}}=\sum\limits_{i=0}^{N-1}{{{X}_{i}}}{{e}^{-j2\pi {{f}_{k}}{{t}_{i}}}}=\sum\limits_{i=0}^{N-1}{{{X}_{i}}}{{e}^{-j2\pi {{f}_{0}}i\Delta t}}{{e}^{-j2\pi k\Delta fi\Delta t}}$

(8)

定义

$\begin{array}{*{35}{l}} \begin{align} & {{\varphi }_{0}}=2\pi {{f}_{0}}\Delta t\text{ }\Delta \varphi =2\pi \Delta f\Delta t \\ & A={{e}^{j{{\varphi }_{0}}}}~Z={{e}^{j\Delta \varphi }} \\ \end{align} \\ \sum\limits_{i=0}^{N-1}{{{X}_{i}}}{{A}^{-i}}{{Z}^{-ki}}=\sum\limits_{i=0}^{N-1}{{{X}_{i}}}A{{Z}^{k}}^{-i} \\ {} \\ \end{array}$

(9)

式中：AZ^k代表了z平面单位圆半径；φ₀代表了从起始频率f₀处的相位；Δφ代表了沿z平面单位圆上频率增量Δf所引起的相位增量，它对应于k的增量，如图 1所示。

由式(6，9) ，有

$\theta ={{\varphi }_{0}}+k\Delta \varphi =2\pi \Delta t({{f}_{0}}+k\Delta f){{f}_{0}}=0$

(10)

假设${{f}_{0}}=0,{{\varphi }_{0}}=0,A=1,\Delta f=\frac{1}{2T}=\frac{1}{N\Delta t}$，则有$\Delta \varphi =\frac{2\pi }{N},Z=\text{ }{{e}^{j2\pi /N}}$。

由于Chirp-Z变换在所选频率范围内有非常好的频域分辨率，所以又称为可缩放傅里叶变换，可以将所有计算频率点控制在感兴趣的频带范围内。

1.2 低阶等效系统辨识方法 1.2.1 方程误差法(EEM)

EEM(图 2)包括了过程误差，但假定所有输入和状态变量测量都没有误差。以俯仰角速度对纵向杆力(杆位移)传递函数为例，有

${{q}_{s}}=\frac{as+b}{{{s}^{2}}+cs+d}{{e}^{-\tau s}}{{F}_{es}}$

(11)

式(11) 可改写为

${{q}_{s}}=\frac{{{e}^{-ts}}{{F}_{es}}sa+{{e}^{-ts}}{{F}_{es}}b-s{{q}_{s}}c-{{q}_{s}}d}{{{s}^{2}}}+\varepsilon ={{X}^{T}}\theta +\varepsilon \text{ }$

(12)

式中

$X={{\left[ {{e}^{-\tau s}}{{F}_{es}}/s{{e}^{-\tau s}}{{F}_{es}}/s\text{ }2-{{q}_{s}}/s-{{q}_{s}}/{{s}^{\text{ }2}} \right]}^{T}},\theta ={{\left[ abcd\text{ } \right]}^{T}}。$

代价函数为

$J=\frac{1}{2}{{\varepsilon }^{T}}\varepsilon =\frac{1}{2}{{({{q}_{s}}-{{X}^{T}}\theta )}^{T}}\cdot \text{ }({{q}_{s}}-{{X}^{T}}\theta )$

(13)

代价函数对参数的梯度必须满足

$\frac{\partial J}{\partial \theta }=-q_{s}^{T}{{X}^{T}}+\theta T\left( XXT \right)=0$

(14)

从而有

$\theta ={{(X{{X}^{T}})}^{-1}}Xq{{~}_{s}}$

(15)

即可得到待辨识参数的值。方程误差法属于一步迭代，其辨识结果可作为输出误差法迭代的初值。

1.2.2 输出误差法(OEM)

OEM(图 3)假定没有过程误差，其误差来源于观测的输出值。仍以俯仰角速度对纵向杆力(杆位移)传递函数为例，其输出误差法的代价函数为

$J=\frac{1}{2}{{({{q}_{s}}{{-}_{s}})}^{T}}({{q}_{s}}-{{q}_{s}})$

(16)

式中：q_s为模型输出值，q_s为观测值。

输出误差法的迭代算法为

${{\theta }_{k+1}}={{\theta }_{k}}-{{M}^{-1}}~\frac{\partial J}{\partial \theta }\left| _{k} \right.$

(17)

式中：$\frac{\partial J}{\partial \theta }=-{{\left( \frac{\partial {{q}_{s}}}{\partial \theta } \right)}^{T}}\left( {{q}_{s}}{{-}_{~}}{{q}_{s}} \right),\text{ }M=\frac{{{\partial }^{2}}J}{\partial \theta \partial {{\theta }^{T}}}\left| _{k} \right.$。对矩阵M可做如下近似处理

$M=\left( \frac{{{\partial }^{2}}J}{\partial \theta \partial {{\theta }^{T}}} \right)\left| _{k} \right.\approx \left( \frac{\partial {{q}_{s}}}{\partial \theta } \right){{~}^{T}}\left( \frac{\partial {{q}_{s}}}{\partial \theta } \right)$

(18)

从而可得到如下迭代算法

$\Delta \theta ={{\left( \left( \frac{\partial {{q}_{s}}}{\partial \theta } \right){{~}^{T}}\left( \frac{\partial {{q}_{s}}}{\partial \theta } \right) \right)}^{-1}}\left( \left( \frac{\partial {{q}_{s}}}{\partial \theta } \right){{~}^{T}}\left( {{q}_{s}}{{-}_{~}}{{q}_{s}} \right)~ \right)$

(19)

输出误差法的辨识精度高，结果一致性好，可作为辨识的最终结果。

2 算例

电传控制系统可有效增加飞机在全包线内的阻尼比，如某型飞机在某些状态下的纵向阻尼比高达0.6～0.8，抑制了飞机固有模态的充分激发。在传统的试飞输入方式下(如倍脉冲) ，不能在全频带范围内给飞机充分的激励，试飞数据所包含的信息比较有限，需要采用3-2 -1-1、扫频等激励方式才能获得频域辨识所需的充分的数据量。

以美国F-14(Tomcat)双发飞机陆用舰载超声速战斗机为例，在海平面以马赫数0.18飞行时(飞机在该状态下的飞行品质等级为1级)，增稳飞机高阶系统传递函数为

$\frac{q\left( s \right)}{{{F}_{e}}\left( s \right)}=\frac{1.034\left( 0.444 \right)\text{ }\left( 0.5 \right)\text{ }\left( 1.887 \right)\text{ }\left( 13.986 \right)}{\left[ 0.7,\text{ }1.05 \right]\left( 0.531 \right)\text{ }\left( 1.48 \right)\text{ }\left( 14.9 \right)\text{ }\left( 18.87 \right)}$

(20)

传递函数中(a)代表(s+a)，[ξ, ω]代表s²+2ξωs+ω²。

选用的LOES拟配模型为

$\frac{q\left( s \right)}{{{F}_{e}}\left( s \right)}=\frac{{{b}_{1}}s+{{b}_{0}}}{{{s}^{2}}+a{{~}_{1}}s+{{a}_{0}}{{e}^{-{{\tau }_{e}}s}}}$，拟配参数有：b₁，b₀，a₀，a₁，τ_e。飞机等效短周期频率ω_ne.sp、阻尼比ξ_e.sp 和时间常数T_θ2与拟配参数之间的关系如下

$\left\{ \begin{matrix} {{\omega }_{ne.sp}}=\sqrt{{{a}_{0}}} \\ {{\xi }_{e.sp}}={{a}_{1}}/\left( 2\sqrt{{{a}_{0}}}~ \right) \\ {{T}_{\theta 2}}={{b}_{0}}/{{b}_{1}} \\ \end{matrix} \right.$

(21)

本文采用数值仿真的方法，模拟“3-2-1-1”和扫频激励输入，使用本文方法对仿真结果进行辨识，结果如图 4～8和表 1所示。

由图 8和表 1可以看出，对于不同的试验方法，本文所用方法计算结果是一致的，品质评价结果与试飞结论相同。

为了对比说明本文的方法优于传统方法，还是采用“3-2-1-1”法激励和扫频法激励的输入输出信号，考虑实际采集数据中会存在干扰，在输出信号中加入了不同幅值的白噪声，分别利用本文方法和某基于最小二乘的优化辨识方法进行辨识，结果见表 2。可以发现超过5%的噪声后，后者可靠性越来越差，尤其是时间延迟，而本文方法辨识出的结果可靠性好，从表 2可以看出本文计算的阻尼比、频率、时延等参数结果稳定，甚至是加入了30%的噪声数据，辨识精度仍然很高。

工程实践表明，“3-2-1-1”、扫频等长时激励方法对试飞员驾驶技术较高，特别是扫频，对试飞员的体力要求也很高，因此，在飞机设计阶段可以在电传系统中增加信号发生器(系统)，专供试飞使用，减轻试飞员负担，有效提高试飞数据质量。对于无法使用信号激励系统的飞机，除了扫频等方法外，还可以采用连续倍脉冲方法(3次以上)、连续倍阶跃或组合来激励，但纵向与横航向操纵一定需要解耦，尽可能不引入操纵以外的响应。

3 结束语

本文提出利用Chirp-Z变换和EEM、OEM相结合频域辨识方法正确，辨识精度高，收敛速度快，能有效解决现代电传飞机大阻尼比状态飞行品质评估难题，在实践中，系统阻尼比达到0.95时，仍然能有效辨识。此外，该方法可以实现实时/准实时在线参数辨识，适用于中国各类军、民机试飞工程实践。