适用于部分试验数据缺失的气动参数辨识方法
刘洋, 李春娜, 方远, 龚春林     
西北工业大学 航天学院, 陕西 西安 710072
摘要: 为解决试验中部分数据缺失或难以测量的问题, 提出了一种适用于部分试验数据缺失的气动参数联合辨识方法。该方法将气动参数辨识问题转为优化问题, 以飞行状态初值和气动参数插值表为设计变量, 使用包含全部气动参数的弹道模型, 构建包含多条数据的目标函数。优化中以气动参数数据库和现有方法辨识结果为先验知识, 拟合出未测量数据的初值作为基准值, 设计了可行样本选取方法, 利用差分进化算法进行求解。应用所提方法处理264条试验数据, 结果表明相比于现有气动参数辨识方法, 所提方法能辨识全部气动参数, 准确度更高, 且能反算出未测量的飞行试验数据, 具有实际工程意义。
关键词: 气动参数辨识    联合辨识方法    试验数据缺失    差分进化算法    全局优化    

气动参数辨识是提取飞行器气动参数、预测飞行状态的一种手段,在外形设计、飞行轨迹预测、弹道重构、射表编制等过程中起着重要作用。气动参数辨识可分为在线辨识和离线辨识。在线辨识是在飞行中进行参数辨识,目的是预报飞行状态和进行自适应控制。离线辨识是在飞行试验后进行参数辨识,目的是建立准确的“气动参数-运动状态”数学模型,同时包含模型辨识和参数辨识。模型辨识是在一系列候选模型集中选择出与系统输入-输出特性最匹配的数学模型[1],参数辨识是优化所选数学模型中的参数。

离线辨识中,给定试验数据并确定数学模型之后,气动参数辨识问题就转化为最优化问题,设计变量为飞行状态初值和气动参数,目标为最小化计算弹道和测量弹道之间的误差。根据试验数据数量的不同,可分为单条和多条试验数据气动参数辨识。

单条试验数据气动参数辨识发展成熟,重点研究不同飞行器的建模及算法匹配[2-4]。文献[5]对比了输出误差法(OEM)和两步法(2SM)在飞机气动参数辨识中的表现,结果表明OEM更准确,但2SM效率更高更灵活。文献[6]利用极大似然法(MLE)对飞机操稳特性大导数进行辨识,并与小扰动理论分析法进行对比,结果表明随机噪声对大导数的影响较小,但对长周期和滚转模态响应值影响较大。文献[2]提出了对飞行器扰流片参数辨识的元优化算法(MO)及初始样本选取方法,并对比了差分进化算法(DE)等其他9种算法,结果中仅有DE、MO和蚁群算法收敛至全局最优解。

多条试验数据参数辨识方法重点研究多条试验数据间的关系。因为针对飞行数据的气动参数辨识通常存在局限性、差异性和矛盾性。局限性表现为在单条试验数据的辨识中,所得结果仅对当前测量数据有效;差异性表现为受飞行器个体差异、环境、测量等影响,不同数据及其辨识值本身存差异;矛盾性表现为局限性和差异性会互相放大,从而无法将辨识值直接用到其他飞行试验上时。因此,按照这3种特性处理方式的不同可分为2类:①每条试验数据单独辨识,对多个辨识结果取平均值并进行符合修正,作为多条试验数据辨识结果;②同时对多条试验数据进行辨识,获得唯一的辨识结果。

第一类辨识方法原理简单,方法成熟,计算效率高,是射表编拟的标准方法,目前得到广泛应用。文献[7]使用极大似然法对8条导弹靶道试验数据进行辨识,利用辨识结果的均值分析非对称导弹的运动特征。文献[8]利用自适应扩展卡尔曼滤波和Modified Bryson Frazier Smoother(MBFS)对战略导弹的7条试验数据进行降噪,然后利用MLE和Broyden Fletcher Goldfarb Shanno(BFGS)对数据分别进行辨识,所得结果趋势一致,但在数值上仍有较大差异。此类辨识方法存在的主要问题有3个:①可能使用多个数学模型进行辨识[8],不同模型会有不同简化,增大了差异性;②辨识结果误差较大且可能有偏[7],没有完全解决矛盾性;③部分气动参数难以辨识[7-8]

在第二类方法中,根据飞行器及试验类型的不同,建模方法不完全相同。文献[9]将再入体的气动参数表示为马赫数和攻角的多项式,同时对9条测量数据进行辨识,所得结果与对相似飞行器的其他研究一致。文献[10]建立了联合雅克比矩阵,并使用梯度法对2种外形相似导弹的4条测量数据进行了辨识。结果表明,同时辨识4条测量数据比同时辨识1条或2条测量数据时目标函数值更小;使用标准模型获得的目标函数值比使用简化模型更小。文献[11]使用差分进化算法对炮弹的6条纸靶试验数据同时进行辨识,辨识结果的精度比第一类方法更好,更能反映实际飞行状态。相比第一类方法,第二类方法的计算成本更高,对试验数据的要求也更高,但能较好地解决辨识的局限性、差异性和矛盾性,所获的气动参数更能反映飞行器实际飞行状态。

在实际工程中,测量系统制约了气动参数辨识的发展。受限于传感器的安装空间、采样频率、测量精度等问题,试验中常有部分飞行数据无法测量或测量误差较大等问题;在部分试验中,测量数据仅能用于制导控制,无法用于气动参数辨识。以上问题存在于各类飞行器的气动参数辨识中,例如飞机[12]、导弹[13]等,增加了辨识难度,降低了辨识准确度,最终限制了飞行器设计。

针对测量系统中存在的问题,本文根据第二类多条试验数据辨识方法的基本原理,提出了一种适用于试验数据缺失的气动参数联合辨识方法。该方法先构建包含所有气动参数的飞行器运动模型,以气动参数插值表及状态初值为设计变量,以数值计算数据库、风洞试验数据库及第一类辨识方法所得结果为先验知识,获得一个或多个设计变量基准值,并生成设计空间,利用差分进化算法同时对多条试验数据进行辨识。针对部分飞行状态无测量数据的情况,将气动参数的先验知识与已有测量数据相结合,拟合出第一个测量点处的基准值,将其作为设计变量进行优化。该方法能解决并主动利用气动辨识的局限性、差异性和矛盾性,以获得符合全部飞行数据的全局最优解。将以上方法应用于某子弹多条仿真数据的气动参数辨识中,并与现有方法进行对比,以检验所提方法的正确性和有效性。

1 单条试验数据辨识方法

气动参数辨识本质上是选取一套气动参数C, 使飞行弹道计算值Ycal与测量值Yexp之差达到最小。由于飞行弹道的求解需要数值积分, 所以积分初值, 即t0时刻飞行状态初值Ycalt0的获取是气动参数辨识的重点。一般可将测量值Yexpt0作为Ycalt0, 但其测量误差会在积分中累积, 影响辨识的准确性。因此, 通常也将Ycalt0作为待辨识参数, 尽可能减小测量误差的影响。

将以上气动参数辨识问题转化为最优化问题, 设计变量为X={C, Ycalt0}, 约束为g(Ycal, X), 目标函数为f(X)=‖Ycal(X)-Yexp‖, 优化问题可表示为

(1)

图 1为优化流程, 具体介绍如下:

图 1 单条试验数据优化流程

1) 根据Yexpt0获得飞行状态初值Ycalt0的基准值, 将已有的气动参数数据库作为先验知识, 得到气动参数C的基准值。将Ycalt0C合并, 得到设计变量X的基准值;

2) 根据试验数据类型、长度、采样频率, 确定收敛指标ε和总迭代次数kt;

3) 将X代入数学模型中计算并得到飞行数据计算值Ycal。对于无控飞行器及开环辨识的有控飞行器, 该数学模型为飞行器运动方程组, 对于闭环辨识的有控飞行器, 该数学模型为飞行器运动方程组及反馈回路;

4) 根据计算值Ycal和测量值Yexp计算目标函数f(X);

5) 若f(X)>ε且当前迭代次数k<kt, 则计算未收敛, 根据辨识算法获得新的X, 并重复2)~5)。若f(X)≤εk=kt, 则达到收敛条件, 计算终止, X即为优化结果。

不同辨识问题的数学模型不同, 但模型的输入均为Ycalt0C, 输出为YcalYcal为速度矢量v、姿态角速度矢量ω、姿态角矢量Θ和位置矢量R组合构成的矢量, C为各气动参数构成的矢量。每个气动参数可表示为马赫数Ma、攻角α和侧滑角β的插值表或函数, 但在多数情况下, 仅能部分表示为αβ的函数。如对于轴对称飞行器, 阻力系数CD可表示为总攻角αt的多项式

(2)

式中: CD0为零升阻力系数;CD1CD2分别为阻力系数一次项和二次项;o(sin3αt)为更高阶项;CD0, CD1CD2仍为Ma的函数。由于Ma随时间不断变化, 需要分段辨识。

一条试验数据共包含nd个测量点, 将其分为ns段, 每段包含nf个测量点, 相邻两段的测量点可以重复。nf的选取不能过大或过小, 若nf小于设计变量个数, 噪声对优化的影响较大且优化结果不唯一, 若nf过大, 则该段Ma可能变化过大。根据图 1, 每一段试验数据单独构成一个优化问题, 其优化结果为该段平均Ma下的设计变量。

与一般优化问题不同的是, f(X)不能完全反映X与真实值的接近程度。将数据分段后, 优化受测量误差、模型简化等影响更大, 当f(X)很小但X与真实值相差较大时, 优化结果不能反映该飞行器的普遍飞行情况, 一旦工况发生改变, 所得飞行状态与实际可能相差很大。

2 现有多条试验数据辨识方法

飞行试验往往是一类多组、一组多次地进行, 同组试验工况相同, 同类试验的测量系统相同, 不同类试验的测量系统、数学模型均可能不同。

在某飞行器的试验中, 共进行了nm次试验, 共有nm条试验数据, 第i条试验数据可分为nsi段。根据文献[1]中的多条试验数据优化结果取平均值的方法, 共有nt个优化问题

(3)

该方法的误差来源主要有以下3个方面:

1) 在每个优化问题中, 实际Ma和平均Ma有一定差异, 当采样频率较低或Ma变化率较大时, 可能导致低灵敏度设计变量辨识误差很大;

2) 不同试验的可测数据与选用数学模型可能不同, 对同一气动力、力矩的简化程度不同, 因此优化结果在总体上可能是有偏的[14];

3) 部分飞行试验会使用解耦模型, 耦合程度不同, 解耦带来的模型误差也不同, 使同一气动参数的辨识误差较大。

第一种误差可以通过提高采样频率或忽略低灵敏度的设计变量降低对优化的影响[15], 后2种误差无法规避, 误差大小受测量数据、设计变量、采样频率、测量噪声、气动非线性等因素的影响, 难以定量甚至定性分析。因此, 该方法常结合风洞试验数据进行修正[16], 或进行符合计算才能应用。

3 联合辨识方法 3.1 参数辨识方法

为解决现有多条数据气动参数辨识中的问题, 提出一种适用于试验数据缺失的气动参数联合辨识方法。为避免由模型简化引起的误差, 构建包含全部气动参数的数学模型, 以数值仿真、风洞试验数据库及现有气动参数辨识方法为先验知识, 获得C是对Ma的函数并作为基准值; 以Yexpt0作为Ycalt0的基准值, 对于缺失测量值的ω*Θ*, 可根据先验知识拟合出t0时刻的基准值。以最小化计算数据与测量数据的误差平方和为目标, 利用差分进化算法对该问题进行优化。该优化问题可表示为

(4)

式中: μi为权系数;ncμi所属试验类型的总数, 目的是f(X)对各类试验的灵敏度相差不会过大;nwnm条试验数据中的所有测量点

(5)

相比于现有多条试验数据辨识方法, 该方法主要做出3点改变:

1) 设计变量XC的形式与弹道仿真完全相同。在现有气动参数辨识方法中, 设计变量是对应平均Ma下的C, 但在弹道仿真中, C由当前Ma插值得到。由于C最终提供给弹道仿真使用, 为使气动参数更准确, 在联合辨识中将CMa的插值表作为设计变量, 使参数辨识和弹道仿真中飞行状态计算方法完全相同, 减小了因气动参数使用方式带来的误差。

2) 将nt条试验数据作为一个优化问题, 使用同一个运动方程组进行优化。在不考虑非典型力、力矩的情况下, 最能描述飞行状态的数学模型是包含全部气动力、力矩的六自由度模型。在联合辨识中, 使用该模型同时对所有测量点进行辨识, 所得X能反映试验的所有飞行状态, 消除了气动参数辨识的局限性和差异性。

3) 对于无法测量的Yexpt0, 可主动利用气动参数辨识的矛盾性, 将风洞试验值和数值计算值作为基准值开展优化。当Yexp中存在ω*Θ*时, f(X)仅计算已测量数据。由于v对力矩系数的灵敏度很低, 在现有辨识方法中, 可能会出现力矩系数极大地偏离真实值但f(X)很小的问题。在联合辨识中, 即使所有试验中均存在ω*Θ*, 但各试验Θ的初始相位不同, 若设计变量与真实值相差较大, f(X)会明显增大且不满足稳定性约束。从试验的角度出发, 在试验设计之初就应考虑飞行器的力矩系数辨识, 因此np次试验中至少有一部分包含Θ, 此时f(X)对力矩系数的灵敏度增大, 更有利于优化。

在优化前, 需给出X的基准值。C的基准值可结合先验知识及现有辨识结果给出, 而Ycalt0的基准值可由CYexpt0和现有辨识结果共同给出。

公式(4)中优化问题的维度很高, 且C的储存形式为插值表, 难以求解梯度信息, 在计算资源充足的条件下选用无梯度的差分进化算法作为优化算法。该算法适用于非线性问题, 全局搜索能力较强, 其种群的求解可以并行以提高计算效率, 且在气动参数辨识中的表现优于其他启发式算法[2]

3.2 可行样本选取方法

在联合辨识方法中, 初始设计空间对计算准确性和效率有显著的影响。在剔除数值计算数据库、风洞试验数据库和现有辨识结果中不合理的数据后, CmaxCmin分别为3个数据库中出现的最大值和最小值, 因此设计空间的边界[Clow, Cup]为

(6)

式中: c1, c2分别为上、下界放大系数, 不同气动参数的c1, c2不同, 需根据具体问题进行分析。

Yexp的测量误差根据测量系统的不同而不同, 需分别讨论。v可通过雷达、GPS等测量, 系统误差一般小于1 m/s, 随机误差由滤波器过滤, 因此vt0设计空间为[vexpt0-1, vexpt0+1]。R的测量方式和v相同, 但系统误差对优化结果的影响很小, 一般情况下t0 < 1 s, 根据v的系统误差, Rt0设计空间为[Rexpt0-1, Rexpt0+1]。

ωΘ的测量方式包括陀螺仪、地磁传感器、太阳方位传感器、闪光阴影靶道等, 但数据并不总能返回地面[10]。当测量数据包含ωt0Θt0时, 可根据测量仪器的系统误差和随机误差得到设计空间; 当测量数据有滚转角速度γ′, 但不包含其他姿态角数据时, 可以根据现有数据库拟合出αt0, αt1βt0, βt1, 再进行微分得到ωt0, 最后根据拟合精度和求导精度得到设计空间; 当没有滚转角γ时, 可根据飞行器类型及典型飞行弹道判断其γt0; 当无法判断时, 可认为γ=0°, γt0的设计空间为[0°, 359°]; 对于旋转弹来说, 可通过公式计算出口转速γ

(7)

式中: v0为表定初速; η为膛线缠度。该公式在工程中较准确, 但受膛线磨损、发射药温等影响, 可能会有一定偏差。

由于传统方法优化结果可能有偏, (6)式生成的初始设计空间不一定包含真实值, 因此在生成第一代种群后, 不对设计空间的边界进行约束, 使搜索可以在边界外进行[11]。在迭代中, 由于f(X)是多个‖Ycal-Yexp‖的加权和, 因此在计算Ycal的过程中就判断子代是否优于父代, 而不必等该个体完全计算完。即每计算完1条试验数据, 数据就计算1次f(X), 若子代f(X)大于父代f(X), 则直接跳出计算, 能提高计算效率。

3.3 具体流程

在该设计空间内, 以飞行稳定性为约束条件, 对公式(4)中的问题进行优化, 具体优化流程见图 2

图 2 联合辨识优化流程图

1) 根据现有多条试验数据辨识方法, 对nm条试验数据进行优化, 得到nt个优化结果;

2) 优化值作为基准值, 生成Ycalt0的初始设计空间, 根据先验知识拟合出t0时刻的Θ并微分出ω;

3) 根据气动参数的先验知识, 生成3个基准值CB1, CB2CB3, 相互补充并生成C的初始设计空间;

4) 将Ycalt0C合并为X, 利用随机拉丁超立方生成个体数为N的种群;

5) 根据试验的可信度和类型占比, 选取合适的μi;

6) 将X代入数学模型, 判断是否满足稳定性约束并计算出Ycal, 并根据(4)式计算f(X);

7) 判断优化是否满足收敛条件, 若满足, 跳转至8), 若不满足, 通过优化算法得到新的X, 并重复5)~7);

8) 计算结束, X即为优化结果。

图 2中, CB1, CB2CB3分别为数值仿真数据、风洞试验数据和现有方法辨识值。考虑到仅用一种方法难以获得全部的气动参数, 因此3个基准值相互补充。

4 仿真验证 4.1 参数设置

为验证本文所提方法的正确性, 对Sierra的7.62 mm竞技子弹仿真数据应用多条试验数据联合辨识方法进行优化。该公开数据包含0~2.5Ma范围内11个气动参数[17], 历经多次试验和比赛验证。该仿真数据有5个重要特点:

1) 仿真对象转速高, 马格努斯效应强, 提高了低灵敏度设计变量的灵敏度;

2) 在全弹道上, ωγ均不可测, αβ在长距离上难以测量, 但在短距离内可测;

3) 可在试验中调整αt0βt0, 使不同试验的飞行状态差异较大;

4) 不包含推力、科式力、闭环控制系统, 数学模型的输入-输出仅体现其气动特性;

5) 测量系统的系统误差容易标定;

该子弹速度衰减快、角运动频率高, 试验数据的采样频率相对较低, 在气动辨识上难度较大。该子弹的运动模型如下:

(8)

式中: vx, vy, vz分别为地面系下的速度3分量, x, y, z分别为地面系下的飞行器位置, ρ为大气密度, S为参考面积, l为参考长度, d为弹径, m为质量, IxIy分别为极转动惯量和赤道转动惯量, C0C2分别为升力系数导数常数项和二次项, C0C2分别为静力矩系数导数常数项和二次项, CMpα, CMpα0, CMpα2CMpαmax分别为马格努斯力矩系数导数、常数项、二次项和最大值, CMqα为赤道阻尼力矩系数导数, Clp为滚转阻尼力矩系数导数, wxwz分别为纵风和横风, vr为相对于风的速度, α′和β′分别为角速率。

使用六自由度刚体弹道方程组生成264条仿真数据, 其中长距离飞行射程为150 m, 有4个初速v0, 每个v0有4个射角θ0, 每个θ0有6个起始攻角和角速度, 共96条数据; 短距离飞行射程为10 m, 共7个v0和4个θ0, 每个θ0有6个起始攻角和角速度, 共168条试验数据。对于长距离飞行数据, 其测量方式为弹道跟踪雷达, 仅有vR的测量数据。对于短距离飞行数据, 为方便研究, 选择每0.5 m得到1个测量数据。测量噪声服从正态分布(μ, σ), 具体参数设置如表 1所示。

表 1 仿真参数设置
x/m vt0/(m·s-1) θ0/(°) 噪声 可测数据
150 (800, 700, 650, 590) (1, 3, 5, 10) N(0.01, 0.002) vx, vy, vz, x, y, z
10 (783, 730, 680, 620, 580, 530, 490) (1, 5, 10, 15) N(0.01, 0.002) vx, vy, vz, α, β, x, y, z

在长距离飞行试验中, 采样间隔为1 ms, 对于小口径子弹, 更大的采样间隔可能会导致计算发散; 在短距离飞行试验中, 由于照相机的位置固定, 采样间隔为0.5 m, 其数学模型转换遵从(9)式

(9)

(9) 式中将测量数据对t的微分转为对x的微分, 使其积分步长为整数。

4.2 现有方法优化结果

采用现有多条试验数据气动参数辨识方法对全部试验进行辨识, 优化算法为极大似然法。

长距离飞行试验中无攻角数据, 因此设计变量C中仅包含CD0, 优化结果如图 3所示。优化得到的CD0误差较大。但是, 优化值均大于真实值, 且攻角越大误差越大, 相对误差最大的点αt=14.63°。根据文献[1]中的方法求CD0的均值并与真实值对比, 结果如表 2所示。

图 3 长距离飞行试验优化结果
表 2 CD0均值
Ma 均值 真实值 误差/%
1.5 0.411 2 0.397 5 3.45
1.7 0.405 0 0.375 0 7.99
1.9 0.390 8 0.357 5 9.32
2.1 0.358 8 0.344 5 4.16
2.3 0.373 0 0.332 7 12.12

表 2中最大误差为12.12%, 其余Ma下的误差均小于10%。需要注意的是, 由于飞行数据是仿真数据, 为充分研究不同攻角的影响, 起始攻角从大到小均有分布, 但在实际飞行中, 很可能所有的起始攻角都较大, 这会导致均值与真实值相差很大。

短距离飞行试验中包含αβ的数据, 设计变量C={CD0, CD2, C0, C0, CMqα, CMpα}。C2C2的敏感度较低, 若将其作为变量, 会使优化结果严重偏离真实值[14]。由于xt0=1 m, nd=19, 因此无需分段, 设置ns=1, nf=19, nt=1 008。优化结果如图 4所示。

图 4 短距离飞行数据优化结果

图 4中可以看出, 优化结果在趋势上正确, 但数值相差较大, 难以直接使用。CD0均值的最大相对误差为14.17%, 最小为8.14%, 均大于表 2中的优化结果, 但未出现误差非常大的解; CD2的优化结果很集中, 仅1个离群点出现在1.986MaC0的分布与CD0相似, 均值的最大相对误差为1.423 2Ma下的21.63%, 其余Ma下误差均小于16%。

C0是所有优化结果中误差最小的, 其均值的误差均在5%以下。但需要注意的是, 系统误差对C0的影响很大, 在实际工程中要尽可能地消除此误差。CMqα的灵敏度并不高, 均值的最大误差为25.57%, 其余误差均小于20%。由于CMqα是非定常气动力系数, 由赤道阻尼力矩和下洗力矩构成, 难以通过其他方式获取准确的值, 因此气动辨识是获得准确CMqα最重要的手段。

在弹道中, CMpα基本达到限幅, CMpα的优化结果相比于CMpα0更接近CMpαmax, 但CMpαmax难以使用极大似然法辨识, 一般通过拟合得到。综合来看, 辨识值均是有偏的, 主要有以下2点原因:

1) 数学模型为简化模型。在长距离飞行试验中, 仅对CD0进行辨识, 会使结果中包含其他气动参数的影响, 该影响与攻角呈正相关。由于攻角不为0, 因此虽然理论上CD0与攻角无关, 但辨识值有偏且取平均无法消除。

2) 气动参数相互影响。如CD0CD2共同构成了阻力系数, 因此容易出现一个增大一个减小的情况, 因此结果是有偏的。

总体上, 现有辨识方法在缺少部分试验数据的条件下, 所得结果误差较大且无法获得全部的气动参数, 在实际工程中直接使用会有较大的误差。

4.3 联合辨识方法优化结果

根据飞行数据的Ma覆盖范围, C的插值表Ma点为1.4, 1.6, 1.8, 2.0, 2.2和2.4。理论设计变量个数为2 970。Yexpt0的测量精度很高, γ′可通过公式估算出较准确的值, 可缩小对应设计变量的设计空间, 重要设计变量个数为986。

使用CFD方法对的关键气动参数进行计算。网格为结构化网格, 第一层附面层高度为10-6 mm, 求解器为Fluent的密度基求解器, 湍流模型为k-ω的SST两方程模型, 气体为标准理想气体模型。分别对3种网格进行网格无关性验证, 计算工况为Ma=1.6、α=10°, 所得结果如表 3所示。

表 3 网格无关性验证
网格数 CD CL CM
103万 0.562 0.591 0.126
156万 0.594 0.613 0.137
204万 0.598 0.615 0.139

表 3中结果来看, 156万网格和204万网格计算结果相差很小, 而103万网格与156万网格相差较大, 尤其是俯仰力矩系数相差8%。综合计算精度和计算效率, 选取156万网格进行CFD计算, 其表面网格如图 5所示, 计算结果如表 4所示。

图 5 7.62 mm子弹表面网格
表 4 CFD计算结果
Ma CD0 CD2 C0 C2 C0 C2
1.6 0.419 5.783 3.13 13.3 0.816 -1.05
2.0 0.382 5.442 3.14 16.9 0.734 -0.766
2.4 0.344 5.046 3.05 18.2 0.683 -0.802

结合试验数据的优化结果及CFD计算结果, 根据设计空间生成方法, 在剔除远离其他优化结果的几个解后, 获得C的初始设计空间。C2C2仅有CFD数据, 因此取C2C2的最大值的1.5倍为上界, 最小值的0.5倍为下界, 获得初始设计空间。CMpα0CMpαmax的初始设计空间相同, 而CMpα2通过图 4c)中同一Ma下不同CMpα的拟合, 给出一个较大的设计空间。该优化问题的约束为

(10)

式中: Sg为陀螺稳定因子, P为转速项, M0M2为静力矩项, 具体计算方法见(8)式。由于C的形式是插值表, 且γ是时变的, 因此无法在获得新的样本时检验是否满足约束, 仅能在弹道计算中进行检验。若不满足则使用罚函数法构建目标函数并跳出当前计算

(11)

式中: c3为惩罚因子。

对于差分进化算法, 设置代数kt=8 000, 种群数目N=4 900, 交叉概率和变异概率因子均为0.8, 种群初始化方法为随机拉丁超立方。优化过程的收敛性曲线如图 5所示。

图 6 目标函数收敛曲线

图 5中可以看出, f(X)在k=8时比k=0时减少了1个数量级, 但之后f(X)的下降速度逐渐变慢。由于优化问题的非线性很强, 在优化前期, 敏感度较高的设计变量最先接近真实值, 使f(X)大幅度下降, 但在优化中后期, 数百个敏感度较低的设计变量是优化的主要对象, 因此优化效率降低。

由于αβα′和β′的灵敏度为1, 对αβ的优化实际上就是对α′和β′的优化。在所有Ycalt0中, 长距离飞行试验的αt0βt0最难优化, 该变量决定了飞行攻角大小, 且缺少测量数据。因此, 根据表 3中的CFD计算值, 对vt0, vt1vt2进行拟合得到αt0, βt0αt1, βt1, 并进行微分得到αt0, βt0αt0βt0的优化速度近似为f(X)的下降速度, αt0βt0的优化结果如7图所示。

图 7 长距离飞行试验αt0βt0优化结果

图 7中可以看出, αt0βt0基准曲线在趋势上较为准确, 这是f(X)能收敛的重要原因。优化结果的误差大部分在3%~8%, 最大误差为10.4%, 由于该参数的灵敏度较低且缺少测量数据, 该优化结果可以接受。气动参数C优化结果如图 8所示。从图 8中可以看出, 联合辨识方法的优化结果远好于文献[1]中优化方法所得结果。CD0C0最大误差仅为0.43%和2.8%, 理论上, vCD0的灵敏度最高, αβC0的灵敏度最高, 且数学模型精度很高, 因此CD0C0的优化结果最准确。

图 8 气动参数优化结果

CD2最大误差为10.4%, C0C2的最大误差分别为5.1%和8.5%。由于在v的计算中这3个气动参数与αβ有关, 因此准确度比CD0低。

C2的最大误差为6.7%且趋势与真实值完全一致。CMqα曲线规律较弱但误差较小, 这是因为CMqα起抑制角运动的作用。理论上αβ的衰减越明显, 优化结果越好, 短距离飞行试验的αβ衰减量较小, 长距离飞行试验中, vCMqα的灵敏度很低, 因此优化结果规律性较差。

由于子弹转速非常高, 马格努斯效应很强, 因此CMpα0, CMpα2CMpαmax的辨识结果较好。在大部分的飞行数据中, CMpα均达到限幅, 即CMpα=CMpαmax, 因此CMpαmax的优化精度最高, 最大误差仅为3.3%, CMpα0CMpα2需小攻角飞行数据才能得到准确的优化结果, 但攻角越小, f(X)对CMpα的灵敏度越低。理论上, γ′越大或α′和β′越小, CMpα的灵敏度越高, 因此想得到更精准的CMpα0CMpα2非常困难, 需专门设计针对性试验。

由于Clp对应的飞行状态γ′不作为可测数据,现有方法在该情况下无法进行辨识。Clp仅能通过γ′的衰减获得, 短距离飞行试验中γ′几乎没有变化, 而长距离飞行数据vγ′的灵敏度非常低, 因此Clp辨识结果最大误差达到了12.7%, 但规律正确。

图 8b)8h)8i)中初始设计空间没有完全包含真实值, 这说明联合辨识方法有能力在初始设计空间不包含真实值的情况下, 得到接近真实值的优化结果。

5 结论

本文提出了一种适用于试验数据缺失的气动参数联合辨识方法, 并给出基准值及初始设计空间的获取准则。利用7.62 mm子弹仿真数据对该辨识方法进行验证, 并与现有多条试验数据辨识方法进行对比, 得出以下结论:

1) 在11个气动参数中, 极大似然法可以获得6个气动参数。相比于现有方法, 联合辨识方法在满足约束的条件下, 使用同一数学模型对多条不同类型试验数据进行辨识, 能较为准确地获得全部气动参数, 且结果具有唯一性。其中CD0的最大误差仅为0.43%, 比极大似然法误差小了1个量级; 与极大似然法误差最接近的是C0的2.8%, 仍小于极大似然法。在极大似然法难以辨识的5个气动参数中, 联合辨识方法所得结果的误差仍然较小, 对于完全没有对应数据的Clp, 辨识误差仍也仅有12.7%。以上结果表明联合辨识方法消除了气动参数辨识的局限性和差异性, 且能减少测量噪声的影响。

2) 联合辨识方法在部分试验数据缺失的条件下, 通过主动利用气动辨识的矛盾性, 能得到对应气动参数及试验数据辨识值。对于无法测量的攻角, 其t0时刻的辨识值αt0βt0误差大多在3%~8%, 最大误差为10.4%, 表明联合辨识方法具有实际工程意义。

3) 以数值仿真数据、风洞试验数据及现有方法辨识值作为先验知识, 确定设计变量的基准值和设计空间, 可以有效减少弹道计算的发散, 提高优化效率, 确保获得准确的优化结果。

对于导弹等包含三通道舵偏角导数的飞行器, 过强的非线性运动、过多的气动参数和敏感的闭环控制系统都会增大联合辨识的难度, 且该类飞行器通常可以有针对性地进行程控飞行, 较少面临没有测量数据的情况。因此, 联合辨识方法更适用于不包含舵面和闭环控制系统的无控段飞行, 此类对象仅需替换相关的运动方程即可。

参考文献
[1] 蔡金狮. 飞行器系统辨识[M]. 北京: 中国宇航出版社, 1995: 2-5.
CAI Jinshi. Aircraft system identification[M]. Beijing: China Astronautic Publishing House, 1995: 2-5. (in Chinese)
[2] MATTHEW G, MARK C. Smart projectile parameter estimation using meta-optimization[J]. Journal of Spacecraft & Rockets, 2019, 56(5): 1508-1519.
[3] 邵干, 张曙光, 唐鹏. 小型无人机气动参数辨识的新型HGAPSO算法[J]. 航空学报, 2017, 38(4): 120365.
SHAO Gan, ZHANG Shuguang, TANG Peng. HGAPSO: a new aerodynamic parameters identification algorithm for small unmanned aerial vehicles[J]. Acta Aeronauticaet Astronautica Sinica, 2017, 38(4): 120365. (in Chinese)
[4] GAO B B, GAO S S, HU G G, et al. Maximum likelihood principle and moving horizon estimation based adaptive unscented Kalman filter[J]. Aerospace Science and Technology, 2018, 73: 184-196. DOI:10.1016/j.ast.2017.12.007
[5] JOAO O, CHU Q P, MULDER J A, et al. Output error method and two step method for aerodynamic model identification[C]//AIAA Guidance, Navigation, and Control Conference and Exhibit, San Francisco, 2005
[6] 丁娣, 钱炜祺, 汪清. 飞机操稳特性大导数辨识及随机噪声影响分析[J]. 航空学报, 2015, 36(7): 2177-2185.
DING Di, QIAN Weiqi, WANG Qing. Identification of aircraft stability and control characteristics derivatives and analysis of random noise[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(7): 2177-2185. (in Chinese)
[7] FRANK F, BERNARD G, ILMARS C, et al. Flight behavior of an asymmetric body through spark range experiments using roll-yaw resonance for yaw enhancement[J]. Journal of Spacecraft and Rocket, 2015, 54(1): 266-277.
[8] KAR P K, SARKAR A K, UMAKANT J. Aerodynamic coefficients estimation of a flight vehicle from different flight trials under limited measurements[C]//AIAA Atmospheric Flight Mechanics Conference, Portland, 2011
[9] MARIE A, SIMONA D, CLAUDE B, et al. Aerodynamic coefficient identification of a space vehicle from multiple free-flight tests[J]. Journal of Spacecraft and Rockets, 2017, 54(2): 1-10.
[10] BRADLEY T B. Aerodynamic parameter identification for symmetric projectiles: an improved gradient based method[J]. Aerospace Science and Technology, 2013, 30: 119-127. DOI:10.1016/j.ast.2013.07.010
[11] 刘洋, 常思江, 魏伟. 多发弹气动参数联合辨识方法研究[J]. 兵工学报, 2020, 41(5): 890-901.
LIU Yang, CHANG Sijiang, WEI Wei. Study on combined aerodynamic parameters identification using flight data of multiple projectiles[J]. Acta Armamentarii, 2020, 41(5): 890-901. (in Chinese) DOI:10.3969/j.issn.1000-1093.2020.05.008
[12] EUGENE A M. Real-time aerodynamic parameter estimation without air flow angle measurements[J]//Journal of Aircraft, 2012, 49(4): 1064-1076
[13] ARDA A. Aerodynamic parameter estimation of a missile without wind angle measurements[C]//AIAA Atmospheric Flight Mechanics Conference, Atlanta, 2014
[14] 刘洋. 高速旋转弹丸气动参数辨识方法研究[D]. 南京: 南京理工大学, 2020
LIU Yang. Study on Aerodynamic parameters identification of high-spin-stabilized projectile[D]. Nanjing: Nanjing University of Science and Technology, 2020 (in Chinese)
[15] DAWID M, MARIE A, SIMONA D. Global sensitivity analysis for modeling the free-flight behavior of an artillery projectile[J]. AIAA Journal, 2020, 58(7): 3139-3148.
[16] 李金晟, 庄凌, 宋加洪, 等. 适用于工程数据的飞行器气动特性修正框架[J]. 航空学报, 2021, 42(12): 125157.
LI Jinsheng, ZHUANG Ling, SONG Jiahong, et al. Aerocraft aerodynamic characteristic correction framework for engineering data[J]. Acta Aeronautica et Astronautica Sinica, 2021, 42(12): 125157. (in Chinese)
[17] MCCOY R L. Modern exterior ballistics[M]. Atglen, PA: Schiffer Publishing Ltd, 2012.
A combined aerodynamic parameter identification method for missing test data
LIU Yang, LI Chunna, FANG Yuan, GONG Chunlin     
School of Astronautics, Northwestern Polytechnical University, Xi'an 710072, China
Abstract: In order to solve the problem that some test data cannot be measured or the measurement is difficult, a combined aerodynamic parameters identification method for missing test data is proposed. In this method, the aerodynamic parameter identification problem is modified into an optimization problem. The initial value of the flight state and the aerodynamic parameter interpolation table are used as design variables, and the motion equation of the aircraft including all aerodynamic parameters is used as a model to construct an objective function containing multiple pieces of test data information. In the optimization, the aerodynamic parameter database and the identification results of the existing methods are used as the prior knowledge. The initial value of the unmeasured data is fitted as the reference value. Then, the feasible sample selection method is designed. Finally, the differential evolution algorithm is used to solve the problem. The proposed method is used to process 264 pieces of test data, and the results show that compared with the existing aerodynamic parameter identification methods, the proposed identification method can obtain all aerodynamic parameters with higher accuracy and can inversely calculate and obtain unmeasured flight test data practical engineering significance.
Keywords: aerodynamic parameter identification    combined identification method    missing test data    differential evolution algorithm    global optimization    
西北工业大学主办。
0

文章信息

刘洋, 李春娜, 方远, 龚春林
LIU Yang, LI Chunna, FANG Yuan, GONG Chunlin
适用于部分试验数据缺失的气动参数辨识方法
A combined aerodynamic parameter identification method for missing test data
西北工业大学学报, 2023, 41(2): 282-292.
Journal of Northwestern Polytechnical University, 2023, 41(2): 282-292.

文章历史

收稿日期: 2022-06-10

相关文章

工作空间