0 引 言

由于风洞自由飞试验数据采集的特殊性,即测量数据只包含模型位置和姿态角信息,无法获取线加速度、角速率及角加速度的测量量^[1],导致气动力系数的辨识准度受位置信息等测量精度影响很大^[2],故待估气动参数的散布较大,正确评价待估参数的辨识准度十分重要。目前,广泛采用Cramér-Rao界作为参数估计准度的度量,但通常情况下,Cramér-Rao界过于乐观^[3-4]。文献[4]的方法通过对残差进行频谱分析来对Cramér-Rao界进行修正,这需要分析者工程经验十足,且工作量较大。文献[5]首先提出一种方法来处理传统的最大似然估计残差,以便计算出有色残差情况下精确的Cramér-Rao下界,并将其应用于F-18飞行器的纵向机动飞行中^[6]。针对风洞自由飞试验中测量数据被有色噪声污染的实际情况，本文采用一种新的修正协方差方法替代传统方法计算Cramér-Rao界，通过对多次Monte Carlo仿真试验和风洞试验实测数据进行气动参数辨识结果的准度评价，验证了修正协方差方法计算Cramér-Rao界的有效性，其客观反映了参数辨识结果的精准度。

1 气动参数辨识准度评价方法 1.1 气动参数辨识原理

目前工程上应用最为广泛的气动力参数辨识方法是最大似然法(MLE：Maximum Likelihood Estimation)，该方法将参数辨识问题转化为一优化问题，通过优化选取气动力模型参数值，使模型输出与实测值间的偏差达极小^[7]，其形式如下^[8]：

式中：ν(i)为输出误差向量：

式中：为由状态方程和观测方程计算得到的观测量预测值，z(i)为测量向量。R为观测噪声协方差矩阵，当观测噪声的统计特性未知时，采用R的最优估计：

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

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

式中：

其中，定义灵敏度矩阵S，表示观测量对待辨识参数的灵敏度，其表达式为：

1.2 Cramér-Rao不等式

Cramér-Rao不等式是准度估计理论中最重要的结果，它给出参数估计准度可达到的理论极限。在某种意义上，Cramér-Rao不等式给出了数据中所含信息量的度量。对于有偏估计，Cramér-Rao不等式为：

式中：▽_θtrb(θ_tr)表示估计的偏差对参数真值的一阶梯度，M为信息矩阵。

对于无偏估计，b(θ_tr)=0，式(8)化为：

式(9)表明，参数估计值与真值的方差总是大于信息矩阵之逆，只有当信息矩阵无穷大时，才能充分接近真值，即当试验的数据足够多时，式(9)的等式成立。经验证明，只要数据含有系统的最大固有周期的几倍长的数据，式(9)的等号就接近成立^[9-11]。因此式(9)是最大似然参数估计准度的最好度量。在参数估计时，θ_tr是未知的，故以代替θ_tr。

对于飞行器动力学系统模型，当过程噪声分布矩阵为零，观测噪声分布矩阵已知时，信息矩阵M的表达式为：

定义信息矩阵M的逆矩阵为散布矩阵D：

那么，矩阵D的对角线元素d_jj给出第j个参数估计值方差的下界，Cramér-Rao不等式及Cramér-Rao下界表达式为：

Cramér-Rao界是预测准度的方法，这已被气动参数辨识的大量数字仿真所证实。但是对于飞行试验实测数据的参数辨识实践却表明，重复多次飞行试验所得参数估计结果的标准差比Cramér-Rao界高5~10倍^[12-14]。原因在于，采用式(13)估算Cramér-Rao界时，假定了系统噪声是平均分布于直至Nyquist频率的白噪声，而对气动参数辨识有重要影响的噪声却是局限于低频区域的有色噪声。噪声特性的差异导致的真实标准差与Cramér-Rao界的比值为Nyquist频率与系统截止频率之比，约为5~10。因此，采用Cramér-Rao界作为参数估计准度的度量往往偏小，过于乐观，在工程上参数估计的不确定度通常取Cramér-Rao界的5~10倍。

1.3 修正协方差方法

在有色测量噪声下，参数估计的协方差为^[15]：

当最大似然估计算法收敛时，离散噪声协方差阵和观测灵敏度独立于参数估计，于是上式可改写为：

对于有色残差，有：

式中:R_νν(i－j)为输出残差的自相关函数^[16]。其无偏估计为：

将式(16)代入式(15)，得：

式(18)给出了有色测量噪声情况下参数估计的协方差阵，其中，R_νν(i－j)和R均采用式(17)进行估计。

参数估计的标准差为：

2 仿真辨识算例 2.1 白噪声对辨识结果的影响

本节以某10°半锥角尖锥模型为例，具体说明气动参数辨识算法在高超声速风洞自由飞试验中的应用。模型及试验工况的基本参数为：m=0.2217kg，l=0.168 96m，底部直径D=0.06m，总压p₀=514 520Pa，总温T₀=410K，来流马赫数Ma=5，采样频率f=2000Hz。给定模型仿真初始条件:V=7m/s、θ=5°、ψ=5°、γ=0°和1组气动参数，在6个观测量上添加白噪声，其中，位置坐标、俯仰角、偏航角和滚转角的噪声标准差分别为1mm、0.1°、0.1°和0.02°。进行300次Monte Carlo仿真试验，采用最大似然法对其进行气动参数辨识，待辨识参数为。

表 1列出了Monte Carlo试验气动参数估计的统计结果，文中涉及到的角度单位以弧度计。表 1中，第2列“真值”为仿真时所用参数值，第3~4列为300次试验的参数估计值的均值和标准差s，第5~6列为传统的Cramér-Rao界方法给出的标准差的均值σ及其与s的比值s/σ，第7~8列为修正协方差方法给出的标准差的均值σ_c及其与s的比值s/σ_c。

由表 1可见，参数估计的统计标准差s与Cramér-Rao界方法、修正协方差方法给出的标准差均值一致。结果表明，在测量噪声为白噪声情况下，2种方法均较好地反映了辨识结果的分布，均可以作为参数估计准度的有效度量。

2.2 有色噪声对辨识结果的影响

由于对气动参数辨识有重要影响的噪声主要是低频区域的有色噪声，故在6个观测量上添加不同水平的有色噪声。下面简要介绍以气动参数辨识为目的添加的有色测量噪声的生成过程。由于高超声速风洞自由飞试验空间较小，试验时间较短等固有特点，故本文采用的轴对称尖锥模型的纵向短周期运动为主要运动模态，有色噪声的主要功率谱应分布于[0，ω_n]上，大于Nyquist频率的区间上功率谱为0，其余部分采用高斯白噪声描述。ω_n为纵向短周期模态近似频率^[17]，表达式如下：

式中：

根据式(21)计算得出ω_n约为25Hz。有色噪声的窄带部分先由零均值高斯白噪声通过一个如图 1所示的通带边界频率为25Hz、阻带截止频率为50Hz的6阶Chebyshev II型低通滤波器生成，再与一个独立的宽带零均值高斯白噪声相叠加。通过对噪声幅值进行调节，达到所需信噪比。在进行Monte Carlo仿真时，对所有仿真输出量均采用此种方式独立地添加噪声。图 2给出了第300次仿真中叠加到θ上的有色噪声时间序列和频谱图。采用此种方法生成的有色噪声较好地模拟了真实飞行试验数据残差序列的特点，这也是本文采用此有色噪声模型的原因所在。

在6个观测量上添加上文所述形式的有色噪声，噪声水平与2.1节白噪声水平一致，进行300次Monte Carlo仿真试验，并对其进行参数辨识。图 3示出了C_L^α、C_z^β、m_y^ωy和m_z^ωz的估计结果分布，可见对于最大似然估计，当试验次数足够多时，参数估计值的概率分布趋近于以真值为中心的高斯正态分布。

表 2列出了Monte Carlo试验气动参数估计的统计结果。由表 2可见，参数估计的统计标准差s约为传统的Cramér-Rao界方法给出的标准差均值σ-的5倍，表明传统的Cramér-Rao界方法过于乐观。而参数估计的统计标准差s与修正协方差方法给出的标准差均值σ_c大致相当，表明修正协方差方法较为合理，可以作为参数估计准度的度量。

3 实测数据辨识算例 3.1 实测数据算例

上文中仿真算例是将各观测量的仿真时间历程作为实测值，即使对仿真值叠加白噪声或有色噪声，也并不能完全反映真实试验飞行历程，所以本节针对10°半锥角尖锥模型的风洞自由飞试验实测数据开展气动参数辨识研究。

试验是在中国空气动力研究与发展中心超高速空气动力研究所Φ1m高超声速风洞中进行的。共进行5次试验，试验模型构型及工况与2.1节近似，各具体参数如表 3所示。观测量与待辨识气动参数的选取与2.1节一致，主要气动参数的估计结果及其Cramér-Rao界如表 4和图 4所示。图 4中符号“○”表 示参数估计值，分别将采用传统方法和修正协方差方法计算的待估参数的标准差σ和σ_c作为误差带，以红色(左侧)和蓝色(右侧)实线表示。

3.2 结果分析

由图 4可见，各参数采用修正协方差方法计算的标准差σ_c远大于传统方法计算结果σ，说明采用传统方法计算的Cramér-Rao界过于乐观。图 4中，C_D^α2、C_L^α和C_Z^β的散布较大，相应地，表 4中对应气动参数的标准差亦很大，这是由力系数对位置坐标的测量准度十分敏感造成的，与文献[18]中的计算结果一致。实际上，由于力矩导数主要影响短周期运动，力导数影响长周期运动^[19]，而风洞自由飞试验的空间较小，试验时间较短，长周期运动体现不充分，位置坐标随迎角、侧滑角的变化不明显，因此，力导数的辨识存在一定难度。只有在模型上安装过载测量设备的情况下，风洞自由飞试验才能对力导数进行有效辨识。

由图 4可见，第3、4次风洞试验中，动导数m_z^ωz的辨识结果均大于0，与文献^[18]中近似解析解和近似牛顿解不符。造成的原因可能是流场中产生不稳定绕流，使气动参数发散，进而产生动不稳定锥运动，俯仰阻尼导数m_z^ωz的估计值大于0，这种简单细长锥体在小迎角情况下可能会产生锥动失稳的现象，与文献[20-22]一致。因此，在进行辨识结果的统计计算时，不予采纳。则第1、2和5次风洞试验得到的动导数m_z^ωz的估计结果为：，σ=0.01961，s/σ=12.15，σ_c=0.05033，s/σ_c=4.73，该统计结果与文献[2]、[18]和[24]中结果一致。

从表 4可见，采用传统方法计算的Cramér-Rao界过于乐观，远小于由5次飞行试验得到的待估气动参数的标准差，而采用修正协方差方法计算的Cramér-Rao界是传统方法计算结果的3~5倍，且与5次飞行试验得到参数的标准差接近。因此，在测量数据被有色噪声污染情况下，采用修正协方差方法计算得到的Cramér-Rao界是待估气动参数准度的一个有效度量。

可注意到，表 4中部分气动参数采用修正协方差方法计算的Cramér-Rao界虽然优于传统方法计算结果，但仍与参数估计的统计标准差略有差别，这是由于本文中高超声速风洞自由飞试验仅进行了5次，样本量较少，且每次试验所用模型及试验工况等各项条件不可能做到完全一致造成的，故可通过增加风洞试验次数和提高样本量来修正。

4 结 论

本文将基于最大似然理论的气动参数辨识方法推广应用于高超声速风洞自由飞试验中。多次Monte Carlo仿真试验和风洞实测数据的辨识结果表明：

(1) 在风洞试验测量数据被白噪声污染情况下，传统方法和修正协方差方法计算的Cramér-Rao界接近，与多次重复试验的参数估计的统计标准差一致，2种方法均较好地反映了辨识结果的散布程度，均可以作为参数估计准度的有效度量。

(2) 当试验测量数据被有色噪声污染情况下，传统方法计算的Cramér-Rao界远小于参数估计的统计标准差，将其作为参数估计准度的度量过于乐观，而采用修正协方差方法计算的Cramér-Rao界与参数估计的统计标准差十分接近，将其作为辨识参数结果的准度评价，更能准确地反映试验中待估气动参数的散布程度。

(3) 采用Chebyshev II型低通滤波器生成的有色噪声信号能量主要集中于低频段，而在高频段能量较低，此种有色噪声信号较好地仿真了真实的高超声速风洞试验测量噪声，并将其应用于Monte Carlo仿真试验中。

(4) 采用修正协方差方法的Cramér-Rao界的全部计算在时域上进行，避免了对有色残差进行频域分析，因此该方法适用范围广，对使用者的工程经验不作要求。