本发明涉及飞行器飞行试验实测数据处理技术领域，特别是一种轴对称无翼无舵导弹主动段零升阻力系数辨识方法。
背景技术：
为有效拦截高超声速侵袭目标，高速高加速防御性战术导弹的弹道和控制系统设计对最大速度、平均速度预示精准度提出了非常高的要求。当该类型导弹采用轴对称无翼无舵气动布局以及趋零攻角飞行弹道设计时，导弹的零升阻力特性，也即小攻角轴向力特性，直接影响导弹的速度特性。空气动力学上零升阻力一般表达为零升阻力系数与飞行动压的乘积。因此，该类型导弹对零升阻力系数预示精准度提出了非常高的要求。获取零升阻力系数有cfd数值模拟、风洞试验和飞行试验三种方法。cfd数值模拟对底部由高温燃气射流、剪切流、旋涡、分离流以及激波干扰流构成的复杂流动模拟精准度不足，难以准确预示阻力系数；因轴对称无翼无舵导弹零升阻力小，试验成本高和技术难度大等原因，地面风洞试验也难以给出宽空域宽速域高精准度的阻力系数。飞行试验是获取飞行器在真实飞行条件下气动力参数的重要途径，基于导弹多次飞行试验实测数据开展气动参数辨识可获取高置信度零升阻力系数(即小攻角轴向力系数)，其结果对理论计算和风洞实验气动数据的评估以及促进气动、弹道、控制改进设计具有重要意义。
基于飞行试验实测数据的气动参数辨识方法主要有以最大似然估计法、卡尔曼滤波法为代表系统辨识方法和直接理论求解的反算方法等。最大似然法和卡尔曼滤波法实质上是以观测数据和未知参数为自变量的统计估计方法，由于建立的状态及观测方程非线性强、耦合程度高，同时没有针对无翼无舵轴对称导弹进行模型简化导致需估计的参数多，具有耗时长、会发散等缺点。同时由于涉及的理论多，这两类方法对技术人员的要求高，都不易普及。基于直接理论求解的反算方法，从牛顿第二定律出发直接求解给出零升阻力系数，具有计算量小且直观、易于理解，便于在工程实际中推广应用等优点。
当前工程实际中应用的零升阻力传统反算辨识法中存在许多不足：一是，受各种因素影响飞行试验实测数据会在某一范围内剧烈跳动，因此以需求辨识点为基准在实测参数中选取若干点进行反算辨识(称为点辨识)存在随机偏差大的缺点。二是，气动参数辨识常由气动专业技术人员开展，但主动段零升阻力系数辨识与发动机推力密切相关，由于气动、动力分属于两个不同专业，容易割裂，从而导致主动段零升阻力系数辨识结果存在置信度偏低的问题。三是，受弹道飞行参数范围限制主动段零升阻力系数辨识结果只能作为小攻角轴向力系数构成窄幅速域和空域范围轴向力系数模型，存在与现有宽速域和宽空域数据表轴向力系数模型不相容的问题。
技术实现要素：
本发明目的在于针对轴对称无翼无舵导弹主动段零升阻力系数预示与建模需求，提供一种基于飞行试验实测数据直接理论求解的反算辨识方法以及基于多次飞行试验辨识结果和其它数据源融合的具有较高置信度并与现有气动模型相容的主动段阻力系数建模方法，以解决使用系统辨识方法计算量大、会发散、技术要求高、不易普及的问题，解决传统反算辨识方法中存在的点辨识随机偏差大、气动与动力割裂导致辨识结果置信度低以及窄域辨识结果与现有宽域气动模型不相容等问题。
对此，本发明提出一种轴对称无翼无舵导弹主动段零升阻力系数辨识方法，其特征在于，其步骤为：第一步、预处理实测数据，对试验数据进行各种预处理，数据预处理包括采样时刻校正、异常值剔除；第二步、预处理理论数据，针对每次飞行试验获取质量、室压、推力数据，并对上述数据进行插值加密处理；第三步、获取推导参数，按照海拔高度获取动压、马赫数推导参数；第四步、导弹飞行推力修正，给出导弹主动段飞行过程中发动机推力修正值；第五步、数据截取，按照如下两个准则截取数据：1)截取主动段数据，以发动机点火时刻开始到轴向视加速度出现负值结束作为主动段，从实测数据中截取该时间段加速度；2)截取小攻角飞行段数据，截取弹体攻角αk≤2°时间段的数据；第六步、计算零升阻力系数，选取导弹底部面积作为参考面积，重力加速度g，计算零升阻力系数；第七步、零升阻力系数平滑与数据表构建；第八步、多次飞行试验辨识结果融合；第九步、获取其他数据源；获取地面计算或风洞试验给出的主动段宽域ma-h-ca数据表。采用同时模拟外流与底部喷流的计算流体力学方法获取由不同攻角、不同高度、不同马赫数构成的宽域工况下的主动段轴向力系数表ma-h-αk-ca；第十步、与其他数据源融合，形成易用的气动参数模型：在ma-h-αk-ca数据表插值获取与
其中，所述第一步、预处理实测数据包括：
1)轴向视加速度ax，单位m/s2；
2)飞行海拔高度h，单位取km，如果没有该数据，需由遥测数据中导航基准系测量的导弹位置x、y、z并结合发射点海拔高度h0、地球半径r按照公式(1)计算给出，
3)飞行速率v，单位取m/s。如果没有该数据，可由导航基准系导弹速度v的分量vx、vy、vz按照公式(2)计算给出，
4)弹体攻角αk，单位取“°”，由导航基准系导弹速度分量按照公式(3)计算给出，
5)摆动喷管θ，单位取“°”，由喷管y、z两个方向的摆角θy、θz按照公式(4)计算给出，
6)发动机贮室压强pzf，单位取mpa。
其中，所述第二步、预处理理论数据包括：
1)导弹飞行理论质量m，对理论质量进行线性插值加密获得导弹飞行过程每一采样时刻的导弹质量；
2)导弹飞行理论室压pzt，对燃烧室压力理论值进行线性插值加密获得导弹飞行过程每一采样时刻的室压；
3)导弹飞行理论推力ft，对理论推力进行线性插值加密获得导弹飞行过程每一采样时刻的导弹推力；
其中，所述第三步、获取推导参数包括：
1)导弹主动段飞行动压q，所述动压用于计算气动系数，标准大气参数表中大气密度与海拔高度对应，按照第一步中的高度h可查表获取当地密度ρ，按照公式(5)计算给出动压，
2)飞行马赫数ma，气动模型按照飞行马赫数给出，由高度h在大气参数表中查询获取当地音速a后按照公式(6)计算给出，
其中，所述第四步导弹飞行推力修正包括：
按照公式(7)计算给出导弹主动段飞行过程中发动机推力修正值ff，
其中，所述第六步、计算零升阻力系数包括：
取导弹底部面积作为参考面积s，重力加速度g＝9.8m/s2，按照公式(8)计算零升阻力系数ca0，
其中，所述第七步、零升阻力系数平滑与数据表构建包括：
1)采用savitzky-golay滤波器方法对第六步获取的零升阻力系数进行平滑处理，每一滤波窗口选取200点，获取由[ma，h，ca0]构成的数据表；
2)按照δma＝0.1间隔在上述数据表中插值获取飞行马赫数ma＝0和本次飞行试验最大飞行马赫数ma之间ma-h-ca0数据表。
其中，所述第八步、多次飞行试验辨识结果融合包括：
1)重复上述第一步到第七步，获取多条飞行试验弹道反算辨识结果数据；
2)计算每一马赫数下所有辨识结果高度平均值
3)剔除rsd≥0.05的数据条
本发明具有以下优点：
1)采用的基于飞行试验实测数据直接理论求解的反算辨识方法相比系统辨识方法具有直观、计算量小、不发散、技术门槛低、易普及等优点。
2)采用的气动、动力专业耦合辨识方法可将单次飞行试验辨识结果偏差由10％～30％提高至5％～10％，同时采用室压遥测结果与理论室压比对的方法实现飞行试验推力修正具有操作简便、精度满足需求的优点。
3)相比“点”辨识方法，采用的基于滤波光顺平滑方法的整个飞行弹道零升阻力系数“线”辨识方法，可降低随机误差10％以上；线辨识参数结果可与理论结果绘制曲线进行直观比较，利于发现偏差明显区域。同时本专利采用的s-g滤波光滑方法在origin、matlab等常用工具软件中都有集成，具有应用简便的优点。
4)采用的基于数据平均和相对标准差阈值的多条弹道辨识结果数据融合以及取舍评判方法，可将辨识结果置信度提高5％，具有意义明确、易操作、易调整等优点。
5)提出的反算辨识结果与其它计算或试验数据源融合的方法，解决了窄域辨识结果与宽域气动模型不相容的问题，利于工程实际使用，同时给出的采用喷流与外流耦合cfd仿真的数值模拟方法获取宽速域宽空域下的主动段零升阻力系数的建议具有成本低、周期短等优点。
6)本发明已用于某实际导弹型号，最大速度预示误差减小了20％～30％，通过了飞行试验考核。
具体实例方式
以下对本发明的具体实施方式做出详细说明。
本发明提出一种轴对称无翼无舵导弹主动段零升阻力系数辨识方法的技术方案为：
第一步预处理实测数据。
对试验数据进行各种预处理，以尽可能地消除各种确定性误差和随机误差。数据预处理包括采样时刻校正、异常值剔除等。需处理的实测数据包括：
1)轴向视加速度ax，单位m/s2。视加速度不包括引力加速度，由惯导测量，数据精度高。
2)飞行海拔高度h，单位取km。如果没有该数据，需由遥测数据中导航基准系测量的导弹位置x、y、z并结合发射点海拔高度h0、地球半径r按照公式(1)计算给出。
3)飞行速率v，单位取m/s。如果没有该数据，可由导航基准系导弹速度分量vx、vy、vz按照公式(2)计算给出。
4)弹体攻角αk(ti)，单位取(°)。由导航基准系导弹速度分量按照公式(3)计算给出。
5)摆动喷管θ，单位取(°)。由喷管y、z两个方向的摆角θy、θz按照公式(4)计算给出。
6)发动机贮室压强pzf，单位取mpa。
第二步预处理理论数据。
针对每次飞行试验获取质量、室压、推力等理论数据，并插值加密。
1)导弹飞行理论质量m。对理论质量进行线性插值加密获得导弹飞行过程每一采样时刻的导弹质量。因飞行过程质量无法实测，用理论值近似代替。
2)导弹飞行理论室压pzt。对燃烧室压力理论值进行线性插值加密获得导弹飞行过程每一采样时刻的室压。
3)导弹飞行理论推力ft。对理论推力进行线性插值加密获得导弹飞行过程每一采样时刻的导弹推力。
第三步获取推导参数。
按照海拔高度获取动压、马赫数等推导参数。
1)导弹主动段飞行动压q。动压用于计算气动系数。标准大气参数表中大气密度与海拔高度对应，因此按照第一步中的高度h可查表获取当地密度ρ，按照公式(5)计算给出动压。
2)飞行马赫数ma。气动模型按照飞行马赫数给出，可由高度h在大气参数表中查询获取当地音速a后按照公式(6)计算给出。
第四步导弹飞行推力修正。
按照公式(7)计算给出导弹主动段飞行过程中发动机推力修正值ff。
第五步数据截取。
按照如下两个准则截取数据。
1)截取主动段数据。以发动机点火时刻到轴向视加速度出现负值结束作为主动段，从实测数据中截取该时间段加速度以及下文涉及的相关数据。
2)截取小攻角飞行段数据。只有小攻角飞行段可用于零升阻力系数辨识。截取弹体攻角αk≤2°时间段的数据。
第六步计算零升阻力系数。
选取导弹底部面积作为参考面积s，重力加速度g＝9.8m/s2，按照公式(8)计算零升阻力系数ca0。
第七步零升阻力系数平滑与数据表构建。
1)采用savitzky-golay滤波器(s-g滤波器)方法对第六步获取的零升阻力系数进行平滑处理，每一滤波窗口选取200点，获取由[ma，h，ca0]构成的数据表；
2)按照δma＝0.1间隔在上述数据表中插值获取ma＝0和本次飞行试验最大ma之间ma-h-ca0数据表。
第八步多次飞行试验辨识结果融合。
1)重复上述第一步到第七步，获取多条飞行试验弹道反算辨识结果数据。
2)计算每一马赫数下所有辨识结果高度平均值
3)剔除rsd≥0.05的数据条
第九步获取其他数据源。
获取地面计算或风洞试验等其他方式给出的主动段宽域ma-h-ca数据表。一般风洞试验成本较高，可采用同时模拟外流与底部喷流的计算流体力学方法(cfd)获取由不同攻角、不同高度、不同马赫数构成的宽域工况下的主动段轴向力系数表ma-h-αk-ca。
第十步与其他数据源融合，形成易用的气动参数模型。
在ma-h-αk-ca数据表插值获取与
本发明针对轴对称无翼无舵导弹主动段零升阻力系数预示与建模需求，提供了一种基于飞行试验实测数据直接理论求解的反算辨识方法以及基于多次飞行试验辨识结果和其它数据源融合的具有较高置信度并与现有气动模型相容的主动段阻力系数建模方法，以解决使用系统辨识方法计算量大、会发散、技术要求高、不易普及的问题，解决传统反算辨识方法中存在的点辨识随机偏差大、气动与动力割裂导致辨识结果置信度低以及窄域辨识结果与现有宽域气动模型不相容等问题。拟保护的创新点主要包括：采用气动、动力专业耦合的方法，提高了辨识结果置信度，同时采用室压遥测结果与理论室压比对的方法较为简便的实现了飞行试验推力修正；提出了基于滤波光顺平滑方法的整个飞行弹道零升阻力系数“线”辨识方法，解决了“点”辨识方法随机误差偏大的问题，并给出了滤波平滑关键控制参数；针对多条弹道辨识结果数据融合及取舍需求，提出了数据平均以及相对标准差阈值简化评判方法，并结合实际应用给出了阈值建议值；提出了反算辨识结果与其它计算或试验数据源融合的方法，既提高了辨识结果置信度，又解决了气动模型相容问题，利于工程实际使用，同时还给出了采用喷流与外流耦合cfd仿真的数值模拟方法低成本短周期获取宽速域宽空域下的主动段零升阻力系数的建议。