一种基于试飞数据的飞行动力学模型频宽自适应辨识方法

摘要

本发明公开了一种基于试飞数据的飞行动力学模型频宽自适应辨识方法，属于飞机飞行试验和飞行动力学模型辨识技术领域。该方法首先采用平移机动动作至窗函数中心，和复制机动动作、补零圆整化数据预处理技术，以削减窗函数边缘缩减效应，提高频域辨识结果的精度。采用经典的Welch谱估计法计算输入、输出的功率谱密度、相干函数、频域响应，并采用多窗口综合技术，进一步提高频域辨识结果的精度。选定低阶等效的飞行动力学模型和等效拟配的频率点数；根据相干函数和功率谱密度，自适应的确定参数化模型辨识的频宽范围和频率节点位置；优化了最小二乘法频域低阶等效拟配，获得更准确的参数化的低阶等效模型。

说明书

技术领域

本发明涉及一种基于试飞数据的飞行动力学模型频宽自适应辨识方法，特别是低信噪比短时机动试飞数据的辨识方法，属于飞机飞行试验和飞行动力学模型辨识技术领域。

背景技术

对飞行器特定输入下的飞行试验数据进行飞行动力学系统辨识后，可获得应用于系统建模、风洞试验数据校准、飞行控制系统设计和飞行品质评价等多个方面的支持数据，是飞行器研制工程中的一项重要工作。由于飞行试验受到输入激励信号类型、测量设备精度、系统干扰、飞行大气环境等多种因素影响，从飞行试验数据中准确地辨识飞行器动力学模型十分困难。为确保基于试飞数据获得准确的飞行动力学模型辨识结果，需采用具有较强鲁棒性和自适应能力的系统辨识方法。

文献“飞机和旋翼机系统辨识：工程方法和飞行试验案例，北京：航空工业出版社，2012”已商业发布的频域辨识法CIFER软件，主要采用Chirp-z变换和组合窗处理技术来提升频域辨识结果的精度，为得到较准确的辨识结果，并特别建议采用历时数分钟的扫频信号作为飞行试验的输入激励。但扫频输入虽然提高了辨识结果的精度，却与试飞工程的实现性、安全性和经济性相矛盾。为此，在试飞工程中，尤其是试飞初期，期望采用较为简单的短时机动作，如倍脉冲，为优先考虑的激励输入信号。但是，短时激励信号频谱范围窄、信噪比低，即使CIFER等频域辨识软件对短时机动动作，往往也不能得到理想的动力学辨识结果。因此，针对飞行试验工程中，满足特定飞行动力学模型辨识工程应用需求、易于驾驶员实现，同时兼顾安全性与经济性的短时激励信号，获得准确的动力学模型辨识结果，成为飞行试验飞行动力学模型辨识技术领域一个重要研究课题。

发明内容

为了克服现有基于试飞数据的频域辨识法对长时扫频输入激励要求苛刻和短时脉冲、倍脉冲等辨识结果精度低的不足，本发明提供一种基于试飞数据的飞行动力学模型频宽自适应辨识方法。该方法首先采用平移机动动作至窗函数中心，和复制机动动作、补零圆整化数据预处理技术，以削减窗函数边缘缩减效应，提高频域辨识结果的精度。采用经典的Welch谱估计法计算输入、输出的功率谱密度、相干函数、频域响应，并采用多窗口综合技术，进一步提高频域辨识结果的精度。选定低阶等效的飞行动力学模型和等效拟配的频率点数；根据相干函数和功率谱密度，自适应的确定参数化模型辨识的频宽范围和频率节点位置；优化了最小二乘法频域低阶等效拟配，获得更准确的参数化的低阶等效模型。

本发明解决其技术问题所采用的技术方案是：一种基于试飞数据的飞行动力学模型频宽自适应辨识方法，其特点是包括以下步骤：

步骤一、

首先对飞行试验机动动作输入x(t)、输出y(t)数据进行采样、去稳态值、平滑滤波等常规的系统辨识数据预处理；针对Welch功率谱估计窗函数的非线性边缘衰减，对短时间的脉冲、倍脉冲、211和3211输入激励，通过平移时间段至窗函数中心附近抵消边缘衰减效应，进一步的，复制1～2个机动，并在整个时间历程内大致均匀分布，以提高机动动作出现在窗函数的中心的概率；通过补零，进行数据圆整化，得到离散化的输入x(n)、输出y(n)。数据预处理，主要消除了数据测量误差，和窗函数引起的有效信号的非线性畸变，可获得更准确频域辨识结果。

步骤二、

采用Welch谱估计法计算输入的自功率谱密度P_xx(ω)、输出的自功率谱密度P_yy(ω)。确定自功率谱密度P_xx(ω)和P_yy(ω)的最小值P_xx(ω)_min、P_yy(ω)_min及其频率点位置P_xx(ω_Ind)_min、P_yy(ω_Ind)_min，进行最小值归零，即和同时归零频率点的位置确定和的最大值频率点位置然后采用去线性趋势项的方法，将自功率谱密度最大值频率点位置相对于另一侧的非零最小值归零，得到预处理后的输入自功率谱密度输出自功率谱密度在每个频率点上对输入、输出的自功率谱密度和进行数值积分，即和计算面积归一化的输入、输出综合功率谱密度消除单位量纲不同引起的谱密度不可比性。

步骤三、

定义一组w个(通常w为3～5)不同时间宽度的窗函数(哈明窗、海宁窗等)。利用Welch谱估计法，分别对第i个时间宽度的窗函数，计算输入x(n)自功率谱密度输出y(n)自功率谱密度和输入输出的互谱密度

飞行动力学系统g(x)，第i个窗函数的频率响应函数G_i(ω)为：

$G_i\left(\omega \right)=\frac{P_{xy}^i\left(\omega \right)}{P_{xx}^i\left(\omega \right)}$

其中，系统频域响应G_i(ω)的为复数形式：

$G_i\left(\omega \right)=G_R^i\left(\omega \right)+{\mathrm{jG}}_I^i\left(\omega \right)$

其中，为复数的实部，为复数的虚部。

G_i(ω)的幅值|G_i(ω)|和相位∠G_i(ω)可表示为：

$\left|G_i\left(\omega \right)\right|=\sqrt{G_R^{i2}\left(\omega \right)+{\mathrm{jG}}_I^{i2}\left(\omega \right)}$

$\angle G_i\left(\omega \right)={\tan }^{-1}\left[\frac{G_I^i\left(\omega \right)}{G_R^i\left(\omega \right)}\right]$

相干函数为：

$r_i^2\left(\omega \right)=\frac{{\left|P_{xy}^i\left(\omega \right)\right|}^2}{P_{xx}^i\left(\omega \right)P_{yy}^i\left(\omega \right)}$

其中，表示互谱密度模的平方，根据相干函数，对不同窗函数i辨识的高阶频域特性幅值|G_i(ω)|、相位∠G_i(ω)和相干函数按公式

$\left|\bar{G}\left(\omega \right)\right|=\frac{\underset{i=1}{\overset{w}{\Sigma }}{G}_i\left(\omega \right)*{\left(r_i^2\right(\omega \left)\right)}^2}{\underset{i=1}{\overset{w}{\Sigma }}{\left(r_i^2\right(\omega \left)\right)}^2}$

$\angle \bar{G}\left(\omega \right)=\frac{\underset{i=1}{\overset{w}{\Sigma }}\angle G_i\left(\omega \right)*{\left(r_i^2\right(\omega \left)\right)}^2}{\underset{i=1}{\overset{w}{\Sigma }}{\left(r_i^2\right(\omega \left)\right)}^2}$

${\bar{r}}^2\left(\omega \right)=\frac{\underset{i=1}{\overset{w}{\Sigma }}{r}_i^2\left(\omega \right)*{\left(r_i^2\right(\omega \left)\right)}^2}{\underset{i=1}{\overset{w}{\Sigma }}{\left(r_i^2\right(\omega \left)\right)}^2}$

进行加权综合，获得多窗口综合的高阶频域响应特性和

步骤四、

根据参数最小化原则选择低阶等效的飞行动力学模型(如品质评价时，品质规范推荐的飞行动力学模型)，设定等效拟配的频率点个数N_ω(如品质评价时，在0.1～10rad/s的频率范围内，拟配频率节点个数一般取为20)和拟配频率范围[ω_min,ω_max]。根据相干函数准则，在[ω_min，ω_max]内，缩减高阶频域响应可用于低阶等效的频宽范围

对面积归一化的综合功率谱密度在频率范围内进行数值积分，即：

$\left(\begin{array}{cc}{S}_P\left({\omega }_k\right)=\underset{i=1}{\overset{k}{\Sigma }}\bar{P}\left({\omega }_i\right)\Delta \omega & {\omega }_i\in [{\omega }_{\min }^{LOES},{\omega }_{\max }^{LOES}]\end{array}\right)$

其中S_P(ω_k)为积分值序列，ω、ω_i、ω_k表示离散的频率点，Δω为按对数等间隔频率。

按S_P(ω_k)和积分值等分序列采用一维线性插值法确定对应的频率点位置，获得内与功率谱密度关联的等效拟配的频率节点序列ω_P：

${\omega }_P=[{\omega }_{P1}^{LOES},{\omega }_{P2}^{LOES},\mathrm{...},{\omega }_{{\mathrm{PN}}_{\omega }}^{LOES}]$

对范围内的采用综合功率谱密度同样的积分值等分的方法，确定相干函数关联的等效拟配的频率节点序列ω_r：

${\omega }_r=[{\omega }_{r1}^{LOES},{\omega }_{r2}^{LOES},\mathrm{...},{\omega }_{{\mathrm{rN}}_{\omega }}^{LOES}]$

取ω_P序列的权重为其中ω_r序列的权重为得到综合加权的等效拟配频率节点序列为：

$\bar{\omega }={\omega }_P*W_{{\omega }_P}+{\omega }_r*(1-W_{{\omega }_P})$

步骤五、

采用带阻尼的最小二乘优化算法，对综合加权的高阶频域响应幅值和相位在频率节点序列上，进行频域低阶等效拟配，获得失配参数、低阶等效拟配参数，和参数化的低阶等效模型。

步骤六、

以飞行试验时域输入序列为模型输入信号，对低阶等效拟配的参数化模型进行时域仿真，获得同输入信号下的低阶模型输出响应时间序列与飞行试验的输出时间序列y_h(i)对比验证。飞行动力学模型辨识结果的可靠性评估，采用如下式定义的曲线拟合重合度指标来衡量，评估辨识模型失配的显著性水平。

$fit=100*(1-\frac{\left|\right|y_h\left(i\right)-{\hat{y}}_l\left(i\right)|{|}_2}{\left|\right|y_h\left(i\right)-\underset{i=1}{\overset{n}{\Sigma }}{y}_h\left(i\right)|{|}_2})$

其中n为时间序列的数值个数。

本发明的有益效果是：由于该方法采用平移机动动作至窗函数中心，和复制机动动作、补零圆整化数据预处理技术，以削减窗函数边缘缩减效应，并采用多窗口综合技术，提高了频域辨识结果的精度。在参数化模型辨识过程中，通过引入表征输入机动特征的功率谱密度，和表征输入、输出线性相关性的相干函数的综合加权，可以自适应的选择拟配频宽和频率节点，使得低阶等效拟配过程与输入、输出信号高度相关，优化了最小二乘低阶等效拟配算法，提高参数化模型辨识的精度。该方法能同时兼顾长时间尺度的扫频机动动作，和脉冲、倍脉冲、211、3211等短时机动飞行试验数据的辨识结果的精度，是一种具有强鲁棒性和自适应性的辨识方法。

附图说明

图1为本发明所述的一种基于试飞数据的飞行动力学模型频宽自适应辨识方法流程图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

本发明提供一种基于试飞数据的飞行动力学模型频宽自适应辨识方法，如图1所示，以某飞行器飞行试验数据短时倍脉冲机动为例，给出其飞行动力学模型频宽自适应辨识方法的实施过程。

步骤一、数据预处理。

对飞行试验机动动作数据输入x(t)、输出y(t)进行采样、去稳态值、平滑滤波等常规的系统辨识数据预处理；针对Welch功率谱估计法估计窗函数的非线性边缘衰减，对长时间的扫频输入激励，通过窗口重叠抵消边缘缩减；对短时间的脉冲、倍脉冲、211和3211输入激励，通过平移时间段至窗函数中心附近抵消边缘衰减效应，进一步的，复制1～2个机动，并在整个时间历程内大致均匀分布，以提高机动动作出现在窗函数的中心的概率；通过补零，进行数据圆整化，得到离散化的输入x(n)、输出y(n)。数据预处理，主要消除了数据测量误差和窗函数引起的有效信号的非线性畸变，可获得更准确的频域辨识结果。

步骤二、计算面积归一化的综合功率谱密度。

采用Welch谱估计法计算输入的自功率谱密度P_xx(ω)、输出的自功率谱密度P_yy(ω)。确定自功率谱密度P_xx(ω)和P_yy(ω)的最小值P_xx(ω)_min、P_yy(ω)_min及其频率点位置P_xx(ω_Ind)_min、P_yy(ω_Ind)_min，进行最小值归零，即和同时归零频率点的位置确定和的最大值频率点位置然后采用去线性趋势项的方法，将自功率谱密度最大值频率点位置相对于另一侧的非零最小值归零，得到预处理后的输入自功率谱密度输出自功率谱密度在每个频率点上对输入、输出的自功率谱密度和进行数值积分，即和计算面积归一化的输入、输出综合功率谱密度消除单位量纲不同引起的谱密度不可比性。

步骤三、计算多窗口综合的高阶频域响应和相干函数。

定义一组w个(通常w为3～5)不同时间宽度的窗函数(哈明窗、海宁窗等)。利用Welch谱估计法，分别对第i个时间宽度的窗函数，计算输入x(n)自功率谱密度输出y(n)自功率谱密度和输入输出的互谱密度

飞行动力学系统g(x)，第i个窗函数的频率响应函数G_i(ω)为：

$G_i\left(\omega \right)=\frac{P_{xy}^i\left(\omega \right)}{P_{xx}^i\left(\omega \right)}$

其中，频率响应函数G_i(ω)的复数形式为：

$G_i\left(\omega \right)=G_R^i\left(\omega \right)+{\mathrm{jG}}_I^i\left(\omega \right)$

其中，为复数的实部，为复数的虚部。

G_i(ω)的幅值|G_i(ω)|和相位∠G_i(ω)可表示为：

$\left|G_i\left(\omega \right)\right|=\sqrt{G_R^{i2}\left(\omega \right)+{\mathrm{jG}}_I^{i2}\left(\omega \right)}$

$\angle G_i\left(\omega \right)={\tan }^{-1}\left[\frac{G_I^i\left(\omega \right)}{G_R^i\left(\omega \right)}\right]$

相干函数为：

$r_i^2\left(\omega \right)=\frac{{\left|P_{xy}^i\left(\omega \right)\right|}^2}{P_{xx}^i\left(\omega \right)P_{yy}^i\left(\omega \right)}$

其中，表示互谱密度模的平方，根据相干函数，对不同窗函数i辨识的高阶频域特性幅值|G_i(ω)|、相位∠G_i(ω)和相干函数按公式：

$\left|\bar{G}\left(\omega \right)\right|=\frac{\underset{i=1}{\overset{w}{\Sigma }}{G}_i\left(\omega \right)*{\left(r_i^2\right(\omega \left)\right)}^2}{\underset{i=1}{\overset{w}{\Sigma }}{\left(r_i^2\right(\omega \left)\right)}^2}$

$\angle \bar{G}\left(\omega \right)=\frac{\underset{i=1}{\overset{w}{\Sigma }}\angle G_i\left(\omega \right)*{\left(r_i^2\right(\omega \left)\right)}^2}{\underset{i=1}{\overset{w}{\Sigma }}{\left(r_i^2\right(\omega \left)\right)}^2}$

${\bar{r}}^2\left(\omega \right)=\frac{\underset{i=1}{\overset{w}{\Sigma }}{r}_i^2\left(\omega \right)*{\left(r_i^2\right(\omega \left)\right)}^2}{\underset{i=1}{\overset{w}{\Sigma }}{\left(r_i^2\right(\omega \left)\right)}^2}$

进行加权综合，获得多窗口综合的高阶频域响应特性的幅值相位和相干函数

步骤四、自适应计算等效拟配频率范围。

根据参数最小化原则选择低阶等效的飞行动力学模型(如品质评价时，品质规范推荐的飞行动力学模型)，设定等效拟配的频率点个数N_ω(如品质评价时，在0.1～10rad/s的频率范围内，拟配频率节点个数一般取为20)和拟配频率范围[ω_min,ω_max]。根据相干函数准则，在[ω_min，ω_max]内，缩减高阶频域响应可用于低阶等效的频宽范围

对面积归一化的综合功率谱密度在频率范围内进行数值积分，即：

$\left(\begin{array}{cc}{S}_P\left({\omega }_k\right)=\underset{i=1}{\overset{k}{\Sigma }}\bar{P}\left({\omega }_i\right)\Delta \omega & {\omega }_i\in [{\omega }_{\min }^{LOES},{\omega }_{\max }^{LOES}]\end{array}\right)$

其中S_P(ω_k)为积分值序列，ω、ω_i、ω_k表示离散的频率点，Δω为按对数等间隔频率。

按S_P(ω_k)和积分值等分序列采用一维线性插值法确定对应的频率点位置，获得内与综合功率谱密度关联的等效拟配的频率节点序列ω_P：

${\omega }_P=[{\omega }_{P1}^{LOES},{\omega }_{P2}^{LOES},\mathrm{...},{\omega }_{{\mathrm{PN}}_{\omega }}^{LOES}]$

对范围内的采用综合功率谱密度同样的积分值等分的方法，确定相干函数关联的等效拟配的频率节点序列ω_r：

${\omega }_r=[{\omega }_{r1}^{LOES},{\omega }_{r2}^{LOES},\mathrm{...},{\omega }_{{\mathrm{rN}}_{\omega }}^{LOES}]$

取等效拟配的频率节点序列ω_P的权重为其中ω_r序列的权重为得到综合加权的等效拟配频率节点序列为，

$\bar{\omega }={\omega }_P*W_{{\omega }_P}+{\omega }_r*(1-W_{{\omega }_P})$

步骤五、最小二乘低价等效拟配。

采用带阻尼的最小二乘优化算法，对综合加权的高阶频域响应幅值和相位在频率节点序列上，进行频域低阶等效拟配，获得失配参数、低阶等效拟配参数和参数化的低阶等效模型。

步骤六、高低阶时域响应对比验证。

以飞行试验时域输入序列为模型输入信号，对低阶等效拟配的参数化模型进行时域仿真，获得同输入信号下的低阶模型输出响应时间序列与飞行试验的输出时间序列y_h(i)对比验证。飞行动力学模型辨识结果的可靠性评估，采用如下式定义的曲线拟合重合度指标来衡量，评估辨识模型失配的显著性水平fit。

$fit=100*(1-\frac{\left|\right|y_h\left(i\right)-{\hat{y}}_l\left(i\right)|{|}_2}{\left|\right|y_h\left(i\right)-\underset{i=1}{\overset{n}{\Sigma }}{y}_h\left(i\right)|{|}_2})$

其中n为时间序列的数值个数。