在国内，目前气动力参数的获取主要基于常规风洞试验，飞行力学特性数据主 要基于大气自由飞试验^[1]。风洞自由飞试验作为由常规风洞试验到全尺寸飞行器试飞试验之间的一个重要过渡环节，综合了以上两种试验手段，其拍摄记录了飞行器模型在风洞流 场中“自由飞行”的运动轨迹，既避免了风洞试验中模 型支撑装置对流场和气动参数测量结果的影响，又具有可控性、可重复性、费用低、周期短等优势。因此，风洞自由飞试验有利于获取飞行器气动力参数和飞行力学特性数据库，开展有关飞行器系统辨识的研究^[2]。

由于风洞自由飞试验技术在现代先进飞行器设计中的重要作用，欧美日等相关大学和研究机构均在开展有关该试验技术的研究。德国亚琛工业大 学(RWTH Aachen University)自2008年起开展了利用风洞自由飞实现飞行器系统辨识的研究^[3]；日本东海大学、九州大学自2000年起开展了利用风洞自由飞实现飞行器飞行力学特性研究^{[4, 5, 6]}；目前，我国也在该领域开展了相关研究工作，文献^[7]对有翼飞行器在脉冲风洞中自由飞试验的动导数测量结果进行了初步讨论；文献^[8]讨论了风洞自由飞实验结果精准度的考 核标准；中国空气动力研究与发展中心蔡金狮等人^[9]自1995年起开展了利用最大似然法对弹道靶自由飞试验进行气动参数辨识的研究。近年来，该中心正致力于风洞自由飞试验手段的各项研究工作^[10]。

风洞自由飞试验测量数据中只包含模型位置和姿态角信息，无法获取线加速度、角速率及角加速度的测量量，若对有限的含噪声的位置和姿态角数据进行二阶数值微分，将造成气动系数的严重误差。因此，本文将气动参数辨识技术应用于风洞自由飞试验，利用有限测量数据(x_d，y_d，z_d，θ，ψ，γ)获取飞行器的气动特性，并对此类试验提出几点建议，以提高气动参数辨识结果的可信度。

目前工程上应用最为广泛的气动力参数辨识方法是最大似然法(MLE：Maximum Likelihood Estimation)，该方法将参数辨识问题转化为一优化问题，通过优化选取气动力模型参数值，使模型输出与实测值间的偏差达极小^[11]。对于风洞自由飞试验，图 1给出了气动参数辨识最大似然法的基本原理。

图 1 飞行器气动参数辨识最大似然法的基本原理 Fig 1 MLE for aircraft aerodynamic parameter estimation

图选项

下文针对风洞自由飞试验，给出气动参数辨识最大似然法的数学模型、辨识准则、辨识算法以及辨识准度评价方法。

风洞自由飞试验气动力参数辨识问题的数学模 型由状态方程、观测方程、气动力模型构成。飞行器 六自由度动力学系统包括3个速度分量、3个角速率分量、3个欧拉角和3个位置坐标，共12个状态变量。参照航天工业标准QJ1028A中的弹体坐标系定义，原点o为模型质心，ox₁轴沿模型纵轴，指向锥顶；oy₁轴在纵对称平面内，垂直于ox₁轴，向上为正；z轴垂直于ox₁y₁平面，其方向按右手定则确定。状态方程可如下：

点击浏览公式

式中：D=1-J_xy²/(J_xJ_y)；V_x1，V_y1，V_z1为体轴系速度分量；ω_x，ω_y，ω_z为体轴系角速率分量；θ，ψ，γ为俯仰、偏航和滚转姿态角；x_d，y_d，z_d为大地系位置坐标分量；m为飞行器质量；g₀为重力加速度，取为常值；S为参考面积，l为参考长度；J_x，J_y，J_z，J_xy为惯性矩和惯性积；C_D，C_L，C_Z，m_x，m_y，m_z为速度系下气动力和力矩系数；q_∞为动压，q_∞=1/2ρV_r²，其中 V _r 为相对 速度，由模型速度 V 与来流风速 V _w合成，V _r= V － V _w；α，β分别为迎角和侧滑角，其计算公式如下：

点击浏览公式

点击浏览公式

式中：ν_i i=1,2,…,6 是观测噪声，为高斯白噪声。

在状态方程(1)中含有六分量气动系数。采用最大似然法进行气动参数辨识，需要给出参数化的空气动力数学模型。这里采用多项式模型，其中气动静导数取至3阶以反映非线性特性，动导数只取线性项。气动力系数的数学模型和相应气动参数的具体形式如下：

点击浏览公式

气动力模型(4)所包含的气动参数均为待辨识参数。由于状态变量的初值是未知的，它对最大似然参数辨识结果的影响又较显著，通常将状态变量的初值也作为未知参数参加辨识。因此，待辨识参数为

点击浏览公式

如1.1节所述，本文采用最大似然准则，其形式如下^[12]：

点击浏览公式

点击浏览公式

点击浏览公式

气动参数辨识问题就是寻求待辨识参数θ的最大似然估计值，使准则函数J( θ )达最小值，即：

点击浏览公式

本文采用修正Newton-Raphson迭代算法求解此优化问题。其迭代公式为：

点击浏览公式

本节以某10°半锥角尖锥模型为例，具体说明气动力参数辨识算法在风洞自由飞试验中的应用。模型及试验工况的基本参数为：m=0.2217kg，l=0.16896m，底部直径D=0.06m，相对质心位置 _cg=0.596，总压p₀=514520Pa，总温T₀=410K，来流马赫数Ma=5。

为验证算法的有效性，首先给定模型仿真初始条件，V=7m/s，θ=5°，ψ=5°，γ=0°，和1组气动力参数，由状态方程组(1)仿真得出各观测量的时间历程作为实测值。观测量取俯仰、偏航和滚转姿态角及3个位置坐标。选取以下14个气动参数进行辨识：C_D0,C^α2_D,C_L0,C^α_L,C_Z0,C^β_Z,m_x0,m^β_x,m_y0,m^β_y,m^ωy_y,m_z0,m^α_z,m^ωz_z。

主要气动力参数辨识值与真值(即仿真时用的参数值)的比较列于表 1，文中涉及到的角度单位采用弧度计。应用辨识所得参数对飞行器纵向位置坐标x_d、法向位置坐标y_d和迎角α进行估计，如图 2所示。由表 1和图 2可见，估计结果与仿真测量数据符合得较好，说明辨识结果是可信的。

图 2 估计结果与仿真值的比较 Fig 2 Comparison of estimation and simulation

图选项

在实际应用中，观测量的实测结果中不可避免地会存在测量噪声。将法向位置坐标的仿真计算结果分别叠加64组白噪声后作为测量值，此白噪声的标准差为对应测量值的1%，并计算辨识参数的平均值与标准差，列于表 1中。由表 1可见，法向位置坐标测量精度仅对法向气动参数影响较大，如零升系数C_L0、升力线斜率C^α_L和动导数m^ωz_z等，而对其它气动系数的辨识结果影响甚微，其中，零升系数和升力线 斜率的分布如图 3所示。由图 3可见，第4组试验 中，法向气动参数的辨识结果为C_L0=－0.00072,C^α_L=5.10605，与真值偏差最大，但应用该组辨识结果对主要观测量进行拟合，其与测量值差异很小，如图 4所示。由此可见，法向位置坐标的测量准度对法向力气动系数的影响很大。

图 3 法向气动参数辨识结果 Fig 3 Identified results of the normal aerodynamic parameters

图选项

图 4 估计结果与仿真值的比较(含噪声) Fig 4 Comparison of estimation and simulation (with noise)

图选项

为了进一步观察噪声对辨识结果的影响程度，加大噪声信号，将64组噪声幅值分别改为2%和5%叠加到法向位置坐标的仿真值中作为测量值，辨识结果同样列于表 1。由表 1可见，零升系数C_L0、升力线斜率C^α_L及动导数m^ωz_z等参数辨识结果误差进一步增大。在噪声幅值为1%、2%和5% 3种情况下，C^α_L的辨识误差分别为3.97%、17.14%和41.85%，动导数m^ωz_z的辨识误差分别为0.42%、1.74%和5.4%。因此，为了保证C_L0及C^α_L辨识结果的有效性，噪声幅值应控制在2%以内，即测量精度控制在10mm。在噪声幅值5%情况下，选取其中辨识结果与真值偏离最大的1组气动参数对观测量进行拟合，如图 4所示。由图 4可见，法向位置拟合结果与测量值偏差较大，与零升系数C_L0与升力线斜率C^α_L的辨识结果较差有关。因此，法向位置坐标的测量精度会影响法向力系数的可信度。此外，由于C_L0与C^α_L 2者均作用于升力系数，公式如下：

点击浏览公式

点击浏览公式

因此，风洞自由飞试验中的气动参数C_L0、C^α_L、C_Z0与C^β_Z的辨识结果受观测量精度影响较大，对风洞自由飞试验中观测量的测量精度有较高要求，位置坐标的精度需控制在10mm以内。相对而言，动导数m^ωy_y与m^ωz_z受观测量精度影响较小，位置坐标的测量精度控制在25mm的合理范围内即可。实际上，由于力矩导数主要影响短周期运动，力导数影响长周期运动^[13]，而风洞自由飞试验的空间较小，试验时间较短，长周期运动体现不充分，位置坐标随迎角、侧滑角的变化不明显，因此，力导数的辨识存在一定难度。只有在模型上安装过载测量设备的情况下，风洞自由飞试验才能对力导数进行有效辨识。

上一节算例是将各观测量的仿真时间历程作为实测值，并未考虑到真实风洞试验环境下不可避免的传感器常值误差、传感器安装位置误差、传感器时间延长误差、测量噪声等一系列干扰的影响，即便对仿真值叠加白噪声也并不能完全反映真实试验飞行历程，所以本节针对10°半锥角尖锥模型的风洞自由飞试验实测数据开展气动力参数辨识研究。

试验是在中国空气动力研究与发展中心超高速空气动力研究所Φ1m高超声速风洞中进行的。试验工况、模型构型、观测量与待辨识气动参数的选取与2.1节一致。所得气动参数的辨识结果列于表 2。利用气动参数辨识结果对状态变量进行拟合，与实测数据的比较如图 5所示。

图 5 拟合结果与测量值的比较 Fig 5 Comparison of fitting results and measurement

图选项

同时，本文给出了m^β_y、m^α_z、m^ωy_y、m^ωz_z等气动参数的近似解析解。考虑一种比较理想的简单情况，即在自旋速度很小时，对于轴对称模型，较为简单准确 的线性力学模型是俯仰角θ与偏航角ψ组成的常系数微分动力学方程组，公式如下^[14]：

点击浏览公式

点击浏览公式

由于模型为轴对称体，因此方程组(13)有解析解：

点击浏览公式

由方程组(14)出发，分别通过观测量 t_i,θ_{i i=0,1,…,n}和 t_i,ψ_{i i=0,1,…,n}，采用遗传算法(GA：Genetic Algorithms^[15])识别出方程中m^β_y、m^α_z、m^ωy_y、m^ωz_z等气动参数和姿态角的拟合曲线，如表 2和图 5所示。图中，“Identified”和“Analytical”分别表示利用辨识结果和近似解析解得到的姿态角拟合曲线。此外，表 2中也列出部分气动参数近似牛顿解，由气动力工程计算软件DATCOM^[16]得出。

综合比较表 2与图 5结果，下面对模型的气动参数辨识结果进行简要分析。

(1) 静导数m^α_z与m^β_y辨识结果近似相等，与DATCOM计算结果比较，相对误差分别为5.51%和28.41%，基本满足工程精度要求。造成m^β_y误差较大的主要原因是偏航角实测数据质量较差，如图 5(b)中偏航角多次出现尖点，对航向静导数辨识结果造成较大影响。

(2) 俯仰阻尼导数m^ωz_z与DATCOM计算结果相比偏小26.09%，文献[17]、[18]中辨识结果也比近似牛顿解计算结果偏小。

(3) 偏航阻尼导数m^ωy_y辨识结果比DATCOM计算结果偏小约1倍。造成大误差的原因主要有2方面：一方面，DATCOM等工程计算软件算法精度有限，下一步考虑开展计算流体力学研究，并应用其计算结果进行比对；另一方面，偏航角实测数据质量较差，且航向静导数误差已较大，进而对航向动导数辨识结果造成很大影响。

(4) 在不考虑模型自旋情况下，最大似然辨识方法与近似解析方法的理论基础相同，两者得到的气动参数结果一致。但当考虑模型自旋时，将计入模型进动和章动的影响，解析方法将成为三周期法，而最大似然辨识方法仅增加一个待辨识参数，方法仍适用。

(5) 利用气动参数辨识结果对姿态角进行拟合，与实测数据一致性较好，说明辨识结果具有一定可信度。 4 初步结论

通过仿真辨识算例和实际风洞试验辨识算例结果可得出以下初步结论：

(1) 静导数辨识结果m^α_z与m^β_y近似相等，且与DATCOM计算结果接近，动导数m^ωz_z变化趋势与国外文献发表结果一致。说明将基于最大似然法的气动参数辨识技术应用于风洞自由飞试验中，是获取飞行器气动特性的有效途径之一。

(2) 位置坐标的测量误差将造成长周期气动参数C^α_L、C^β_Z等结果的可信度较低，甚至不可辨识。因此，可提高测量数据精准度(通过提高高速摄像设备的性能实现，如提升拍摄速度与增多像元素)，如法向、侧向位置坐标的测量精度控制在10mm以内，有效提高参数辨识结果的可信度。

(3) 动导数m^ωy_y与m^ωz_z的可辨识性较高，观测量的测量精度控制在25mm的合理范围内即可。此外，在模型上安装过载测量设备，增加加速度等测量量，亦可提高待辨识气动参数的可信度。