脉冲风洞天平短时振荡测力数据稳态值提取的优化识别方法

王 锋， 武 龙， 吴东升， 乐嘉陵

(中国空气动力研究与发展中心 超高速所 高超声速冲压发动机技术重点实验室， 四川 绵阳 621000)

激波风洞、炮风洞和燃烧风洞等脉冲式风洞是开展高超声速飞行器测力试验的重要设备。由于脉冲式风洞在极短时间内启动，气流以近似阶跃形式作用于试验模型，导致试验系统结构振动，在几十或几百毫秒的有效试验时间内，天平测力信号呈现大幅度振荡，无法像连续式风洞那样等衰减到稳定值再采集数据。因此，对于脉冲式风洞测力试验，从振荡信号中提取隐藏其中的稳态值是一个比较困难的问题。

对于这种短时振荡测力信号，可以采取在试验稳定段内取均值的方法来得到稳态测力值[1-2]。这种方法比较简单，但理论上只有当测力信号是缓慢衰减的正弦波，可以清晰准确地获取到整数个周期的数据，才能得到较准确的均值。然而，实际中的测力试验系统大多是复杂的结构振动系统，特别是对于较大尺度的飞行器模型，本身低频模态比较密集，在试验时很多个模态被激发，测力信号形态复杂，难以判断截取整数个周期，而且也几乎不可能在所取时间段内，使参与振动的主要模态正好都取到整数个周期。这时，对信号取平均值作为稳态载荷测量值就有可能带来较大误差。当然，系统的基频越高，稳态时段上的振动周期越多，误差就会越小。但是，对于接近真实飞行器的大尺度试验模型，质量很大(几百千克，甚至超过一吨)，在保证天平灵敏度的前提下，试验系统的基频难以做到很高，在有效时间段内，测力信号中低频成分的周期数可能只有3～5个甚至更少，这时平均法的误差就可能比较大。因此，有必要采取新的方法从振动信号中提取可靠的稳态值，而不必对试验系统的基频提出过高要求。

1 问题的求解模型

根据脉冲式风洞的工作特点，其理想的工作过程可以描述为：在发出启动信号后，在短暂时间内(几毫秒至几十毫秒)，风洞的喷管迅速建立起稳定的流场，并保持约几百毫秒，然后流场快速减弱消失，如图1所示。从模型的角度看，其在静止状态突然受到强大的瞬态气动载荷的冲击产生振动，然后持续受到稳定气动载荷的作用，最后有个突然的卸载过程。可把此过程区分为加载段、保持段和卸载段。通常只关心保持段。

图1 脉冲风洞理想工作过程 Fig.1 Ideal operation process of impulse tunnel

在加载段，试验系统结构经历复杂动态载荷作用下的受迫振动过程，其结果是为保持段提供了初始振动位移和速度。在保持段，试验系统在前段提供的初始扰动下产生自由衰减振动，同时在稳态气流作用下产生稳定的气动载荷。测力的目的就是得到在保持段被振动信号掩盖的稳态输出数值，最后利用标定的天平方程计算载荷。至于卸载段则可以直接忽略。

从结构上，试验系统一般由试验模型、天平和支撑件组成。试验中，系统仅允许微幅振动，因此可视为线性振动系统，可用模态叠加法来描述其振动运动，而每个模态就是一个单自由度振动系统。模型测力一般同时测量多个分量，因此天平有多路输出，但数据处理时可逐个处理，求得每一路的稳态值，最后再一起代入天平方程求模型载荷。在此仅考虑一路输出，基于前面的分析，根据单自由度系统的自由衰减振动响应和模态叠加原理[3]，在保持段，系统由初始扰动引起的理论输出可写为

(1)

式中：时间t以该段起点为零时刻，m是选取的系统振动模态数目，Ai、ξi、fi和φi分别是第i个模态的初振幅、阻尼比、频率和初相位，A0是稳态输出。这些参数均是未知的，但所关心的主要是A0，如果获得了每个输出通道的A0，利用天平方程就可计算出模型气动载荷的各个分量。因此，数据处理的任务就是识别参数A0，当然其它参数也会同时获得。

(2)

(3)

J中包含4m+1个未知参数，其显式的表达式如下

(4)

下面讨论如何求解该优化问题以获得未知参数。

2 参数识别算法

式(3)是一个非线性优化问题，可以通过适当的优化算法将未知参数辨识出来。目标函数比较复杂而且参数众多，为了能够得到全局最小点并提高求解效率，确定较好的参数初始值十分重要。

2.1 参数初值确定方法

(5)

取Y向量的前半部分，在其幅频曲线上取前m个最大峰值点，这些点对应的物理频率、幅值和相位分别作为m个模态的频率fi、幅值Ai和φi的初值

(6)

式中：fs为信号采样频率；k是峰值点对应的离散频率点序号；YR、YI分别是Y(k)的实部和虚部；sgn(·)是符号函数。

关于模态数量m的确定：取Y的幅值的最大值作为参考量，将各峰值点的幅值与该量的比值作为选择的依据，设定一个合适的阈值cr(比如0.1)，m就取比值大于cr的峰值点的个数。

这些初值的精度与信号本身以及DFT的频率分辨率有关，而频率分辨率由信号时长决定，为

(7)

式中：Δt是采样周期；n为所取测量数据点数。虽然对于确定参数的大致初值来讲，并不需要很精确，但是，如果试验系统的前m阶模态中存在两个模态的频率差小于Δf，则对应振动信号在DFT中将产生主瓣混叠，在幅频曲线上体现为一个峰值，无法将其直接区分，这时对m的确定就有问题。例如，信号时长为100 ms时，仅为10 Hz，此种情况很有可能发生。因此，当数据的时长较短时，为了准确的确定系统重要的模态数目，需要在数据长度有限的情况下，提高DFT的频率分辨率。对于无阻尼或弱阻尼信号，Sacchi等给出的高分辨率DFT方法[4-5]能提高短时数据DFT的频率分辨率，可用来解决该问题。

因为实际结构阻尼比通常很小，ξi的初值就取零。

2.2 优化求解算法

将所有待辨识参数列在一起组成参数向量，令

(8)

解问题式(3)即求满足

(9)

的参数θ。这是一个非线性方程组，在此采用拟牛顿法算法族[6]中的BFGS算法来求解。拟牛顿法避免了在迭代过程中求F(θ)关于θ的导数矩阵，也就是J(θ)关于θ的Hessian矩阵及其逆矩阵，在保持较快收敛速率的同时降低计算量。BFGS算法的迭代流程如下：

(0)准备初值θ0和一个(4m+1)×(4m+1)的正定对称矩阵B(通常取单位矩阵)，设定收敛判断参数εr，计算向量f0=F(θ0)；

(1)在迭代的第i步，计算fi=F(θi)，令

si=-Bifi

(10)

作为下一步的搜索方向；

(2)沿向量si执行一维搜索，求J(θi+αsi)的极小值，得到

(11)

更新参数

θi+1=θi+α*si

(12)

(3)如果‖F(θi+1)‖<εr，取θ*=θi+1，终止迭代，否则继续下一步；

(4)更新矩阵Bi至Bi+1

(13)

其中

di=θi+1-θi，gi=F(θi+1)-F(θi)

(5)执行下一迭代步：i→i+1，转到第(1)步。

3 应用算例

本节用4个算例来对方法的可行性进行检验。第一个算例是人为构造的仿真信号，用来检验辨识算法的准确性和鲁棒性，其余3个算例均针对实际风洞试验的测量信号，检验算法的实用性。

3.1 数值仿真算例

表1 信号参数真值Tab.1 True values of the parameters

本例中，已知一个由多模态振动叠加的信号，并考虑适当的噪声，然后，由此信号辨识出该信号的直流分量和各模态的振动参数。假设信号取为

(14)

式中：ε取标准差为0.2的高斯白噪声，其它参数任给合理的数值，如表1所示，其中第1、2模态的频率间隔仅0.6 Hz。信号时长取200 ms，采样频率取2 kHz，分别用普通DFT和高分辨率DFT对信号进行分析，其幅频曲线如图2所示，因信号时长太短，普通DFT不能区分第1、2模态，而误认信号只含两个模态，高分辨率DFT可以分辨出3个模态。分别用两种DFT结果给出的信号参数初值进行辨识计算，结果列于表 2，可见，尽管高分辨率DFT给出了3个模态的较好初值，但受噪声影响，优化算法并不能准确收敛到三个模态的参数值。但对于我们所关心的直流项，用两模态和三模态辨识结果精度都足够高。用识别参数计算得到的信号与原始信号的对比如图3所示，两模态与三模态的重建曲线几乎都与原始无噪声信号重合。因此，当存在接近重频的信号时，优化算法会收敛到一个合适的“假模态”，使得重建信号能很好拟合原始信号，而信号直流项的辨识精度并不受明显影响。

现考察算法对噪声的鲁棒性。改变信号的噪声水平，每个水平下随机生成50次信号，用两模态进行辨识，计算直流项辨识结果的相对误差绝对值的平均值，其随噪声水平的变化如图4所示，误差随噪声水平近似线性增长，即使在噪声标准差取0.5时，误差也仅有2.5%，而直接对原始无噪信号取平均得到的值则为5.21，其相对误差为4.2%。

图2 普通DFT与高分辨率DFT的幅频曲线 Fig.2 Amplitude-frequency plot of general DFT and high resolution DFT

参数i两模态12三模态123Ai3.21751.00400.02763.23390.9852ci0.62400.13990.83770.66000.4127ωi22.59437.37717.14322.59237.346ϕi0.65060.1193-1.60370.64970.1175A04.99434.9974

图3 重建信号与原始信号的对比 Fig.2 Comparison between reconstructed and original signals

图4 直流项的辨识误差随信号噪声水平的变化 Fig.4 Error of estimated DC value versus noise level of the signal

3.2 简单模型单分量测力试验数据

在脉冲燃烧风洞针对钝头锥模型进行了阻力测量试验[7]，因为钝头锥模型是接近实心的模型，同时支撑系统刚度也很高，而天平仅设计了一个应变测量通道，所以整个测力系统非常接近单自由度振动系统。图5给出了试验测得的风洞燃烧室总压曲线、来流皮托压力曲线和天平输出信号。天平信号在约0.55 s后的突变源于风洞扩压段反射气流的作用，是无效信号。总压和皮托压曲线反应了此风洞的典型运行特征。可见，在所谓稳定段(约0.1～0.55 s)，来流其实不很平稳，皮托压一直有明显波动。因此，在这一段，模型其实依然是受迫振动，只是由于载荷波动较小，系统振动主要由初始的冲击引起。图5中的天平信号可以定性地反映这点，其形态是带有畸变的、振幅逐渐衰减的正弦曲线。

图5 钝头锥吹风试验结果 Fig.5 Test results of a blunt cone model

取0.1～0.5 s的信号作为稳态段信号，辨识得到稳态值为29.744 mV，直接取平均为29.27 mV，二者差别不太大，因为信号的周期较多。用辨识参数重建该段信号，与原始信号对比于图6中。因为实际上该段载荷并不稳定，所以重建信号也不可能像上一例中那样很好地吻合原信号，其偏差主要源于载荷波动引起的振动响应。

图6 天平原始信号、重建信号与稳态值 Fig.6 The original and reconstructed balance signals and the estimated steady-state value

因为载荷不稳定，所以取不同时间段上的数据辨识得到的稳态值应该也会不同。为验证这一点，只取时长为0.1 s的数据，从第0.1 s开始，依次后移0.05 s，直到末端到达第0.5 s，共辨识得到7个稳态值，如图7所示，各段得到的稳态值大小是不同的，其趋势大致与皮托压的变化趋势一致，说明识别方法能够跟踪到一定时段上模型平均载荷的变化。

图7 天平信号分段辨识的稳态值 Fig.7 Piecewise estimated steady-state values of the balance signal

3.3 某升力体飞行器的测力试验数据

现在以某大尺度飞行器模型在脉冲燃烧风洞的测力试验数据来分析。仅取其升力通道的信号，如图8中的实线所示，尽管可以看出信号包含一优势频率，但显然不是单频的简谐振动信号。从风洞燃烧室总压曲线观察，大约0.8～1.1 s的时段，风洞进入稳定的试验状态。取0.825～1.1s时段的数据来辨识该段的稳态值，结果为3.148 mV，在图8中以水平点划线示出，而平均法的结果则为2.68 mV。同时辨识得到该段信号的3个主要频率成分，分别是20.521 Hz、41.646 Hz和72.317 Hz，其幅值分别为7.483 mV、22.917 mV和2.313 mV。用辨识参数重建的信号如图8中的点线所示，基本与原信号重合，说明参数的辨识是准确的。

由于此信号的稳态值与信号振幅相比较小，因此直接取平均的方法容易引起较大的相对误差。

图8 某飞行器测力信号及辨识结果 Fig.8 The original balance signal and the estimated results of a aircraft model

3.4 JF12风洞尖锥模型测力数据

JF12激波风洞[8]由于有效试验时间稍短，约100 ms，则模型的振动周期数较少。以该风洞上10°尖锥模型三分量测力试验数据为例，取频率最低的俯仰通道的数据，如图9中的实线所示。根据皮托管信号，在约0.13～0.23 s上是比较稳定的时段，取该段数据进行辨识，得到稳态值为-0.231 4 mV，以点划线示于图9之中，而直接平均的结果为-0.152 1 mV，由于时间短、周期数少，与辨识结果的差别较大。根据辨识参数重建的信号用点线绘于图9中，与原信号吻合良好，间接说明辨识结果的可信性。另外，天平输出信号中虽含有大量毛刺，但并未对参数辨识产生明显的影响，表明方法有较好的鲁棒性。

图9 尖锥模型测力信号及辨识结果 Fig.9 The original balance signal and the estimated results of a cone model

4 结 论

本文将脉冲风洞平稳段的测力天平信号看作是试验系统在初始条件下引起的围绕未知稳态值的多模态振动信号的叠加，利用优化算法辨识信号主要成分的频率、幅值、阻尼、相位等未知参数，其中零频信号的幅值即为所要提取的天平稳态输出值。这种方法比直接取数据均值的方法更客观和准确，特别是对于风洞稳定段时间短、信号周期数较少的情况，优势更加明显。从另一方面讲，采用此方法后，可适当降低对试验系统基频的要求，即适当降低天平测力元件的刚度，提高天平灵敏度，从而有助于提高测力精度。

但实际中，即使在平稳段，模型载荷也并非恒值，作为进一步的深化研究，可以将该段载荷也进行合理的参数化近似建模，与振动信号参数同时进行辨识，以得到更接近实际的时变载荷。

致谢

感谢中国空气动力研究与发展中心吸气式高超声速技术研究中心贺伟研究员和中国科学院力学研究所刘云峰副研究员提供试验数据。

[ 1 ] 贺伟, 童泽润, 李宏斌. 单模块超燃发动机推力测量天平研制[J]. 航空动力学报, 2010,25(10): 2285-2289.

HE Wei, TONG Zerun, LI Hongbin. Investigation of thrust balance for the single module scramjet[J]. Journal of Aerospace Power, 2010, 25(10): 2285-2289.

[ 2 ] 刘云峰, 汪运鹏, 苑朝凯, 等. JF12激波风洞10°尖锥标模气动力特性试验研究[C]∥ 第十七届全国激波与激波管学术会议, 成都, 2016.

[ 3 ] THORBY D. Structural dynamics and vibration in practice[M]. Oxford, UK: Elsevier, 2008.

[ 4 ] SACCHI M D, ULRYCHlrych T J, WALKER C J. Interpolation and extrapolation using a high-resolution discrete fourier transform[J]. IEEE Transactions on Signal Processing, 1998, 46(1): 31-38.

[ 5 ] SACCHI M D, ULRYCH T J. Estimation of the discrete Fourier transform, a linear inversion approach[J]. Geophysics, 1996,61(4): 1128-1136.

[ 6 ] RAO S S. Engineering optimization: theory and practice[M]. 4th Ed. Hoboken, New Jersey: John Wiley & Sons, 2009.

[ 7 ] 王锋, 任虎, 周正, 等. 载荷辨识方法用于脉冲风洞模型阻力测量研究[J]. 振动与冲击, 2015, 34(24): 202-208.

WANG Feng, REN Hu, ZHOU Zheng, et al. Drag force measurement in impulse facilities by identification method[J]. Journal of Vibration and Shock, 2015, 34(24): 202-208.

[ 8 ] 俞鸿儒. 大幅度延长激波风洞试验时间[J]. 中国科学：物理学 力学 天文学, 2015,45(9): 94-701.

YU Hongru. A big increase in shock tunnel test times[J]. Scientia Sinica Physica, Mechanica & Astronomica, 2015, 45(9): 94-701