CN104061960B - 一种亚音速飞行器体上气压高度参数确定方法 - Google Patents

Abstract

一种亚音速飞行器体上气压高度参数确定方法，包括一套飞行器体上压力源的选取原则，用于压力源的选取定位；一套基于所选定压力源测得压力的数据处理方法，包括修正公式拟合及测量值回归计算方法，修正公式包括α、β修正及静压修正，回归计算方法为根据压力测量值迭代求出目标量的方法。试验结果表明，修正后的数据具有很高的还原度，结合回归计算方法可以得到精度较高的飞行器飞行姿态角及速度等目标量，为飞行器亚音速飞行状态测量提供了一种新方法。

Description

一种亚音速飞行器体上气压高度参数确定方法

技术领域

本发明涉及一种亚音速飞行器上用气压高度参数确定方法，用于飞行器实时解算部分飞行状态参数，提供飞控系统使用。

背景技术

任何飞行器都要求在一定的速度、高度和姿态角下飞行，因此适时确定飞行状态对于飞行器来说具有很重要的意义。测定飞行状态(亚音速)最常规的方法是测量静压、总压和姿态角感知压力(或风标)，采用空速管对来流总压和静压进行测量，利用等熵流关系式即可确定飞行速度，利用静压又可确定出飞行高度，姿态角可利用感知压力修正公式计算(或风标)得到，要想准确的确定飞行姿态常规空速管的等直段通常要做得很长，这样对于现代飞行器尤其是强调隐身性能的飞行器的设计很不利。在常规空速管的基础上发展出了气动补偿空速管，气动补偿空速管的长度和大小相对于常规空速管大大减小，但受飞行器其他各部件影响的程度也相应增加，为了消除亚音速下各部件对其的影响，补偿空速管外形需进行特殊的曲线处理，最重要的则是在测量时对特定状态下所测的压力进行补偿使其更接近真实静压，而要做到这一点则必须事先确定压力随飞行状态的变化规律找出补偿原则，但补偿空速管由于尺寸的限制测压孔数量及位置均存在很大的局限性，加之受干扰较大其补偿规律通常比较复杂，修正效果不佳。

发明内容

考虑到补偿空速管的工作原理，如果把静、总压测压源及姿态角感知压力源直接布在飞行器的特定位置上，并通过试验找出其压力补偿规律，从而准确的计算出飞行状态理论上也是可行的，为此特进行了本发明项目研究，结果表明本方法不仅可行而且具有较高的回归精度，可以替代空速管作为测量亚音速飞行状态参数的一种新手段。

本发明的技术解决问题是：克服现有技术的不足，提供一种亚音速飞行器上用气压高度参数确定方法，在理论计算、风洞试验、数据处理及回归算法推算的基础上提高飞行状态参数的测量精度，为飞控系统的准确控制提供条件。

本发明的技术解决方案是：

一种亚音速飞行器体上气压高度参数确定方法，步骤如下：

(1)对亚音速飞行器进行气动特性数值计算，得到亚音速飞行器体表面的压力分布；

(2)根据步骤(1)中计算得到的压力分布结果，在飞行器体上选择攻角α的测压点、侧滑角β的测压点、总压源以及静压源的位置；

(3)加工所述亚音速飞行器的缩比模型，对步骤(2)中选取的攻角α的测压点、侧滑角β的测压点、总压源以及静压源进行测压风洞试验，得到攻角α角度压力点、侧滑角β角度压力点、总压源和静压源的压力系数试验数据；

(4)将步骤(3)得到的试验数据进行飞行参数修正公式拟合；

(5)对步骤(4)得到的飞行参数修正公式的拟合结果进行相关性计算，根据各干扰项的偏相关系数进一步筛选干扰项，偏相关系数小于0.2的干扰项舍弃，得到各修正公式的复相关系数；

(6)对步骤(5)得到的所述各修正公式的复相关系数进行判断，如果复相关系数的精度大于99％，则保留该复相关系数对应的修正公式，之后进入步骤(7)，否则返回步骤(4)；

(7)如果保留的修正公式大于1个，则选取复相关系数的精度最高的修正公式，并装订到所述亚音速飞行器的计算机中，之后进入步骤(8)；如果保留的修正公式为1个，则直接装订到所述亚音速飞行器的计算机中，之后进入步骤(8)；

(8)进行飞行试验或风洞验证试验，得到静压源压力和总压源压力数据，解算出马赫数，并将该马赫数带入步骤(7)中确定的修正公式计算得到攻角α、侧滑角β、飞行高度以及静压源修正系数，视为一次计算值；

(9)根据步骤(8)中得到的攻角α、侧滑角β、飞行高度以及静压源修正系数的一次计算值再次解算出马赫数，并带入步骤(7)中确定的修正公式解算得到迭代后的攻角α、侧滑角β、飞行高度以及静压源修正系数的目标值即为二次迭代值；

(10)将步骤(9)中二次迭代后的目标值与一次计算值进行误差分析，如果误差满足条件A或者条件B，则回归算法固化为一次迭代；否则返回步骤(9)继续进行下一级迭代直至误差满足要求，届时完成的迭代次数即为回归算法的迭代级数，将该迭代级数装订到所述亚音速飞行器的计算机中，之后进入步骤(11)；

条件A：平均标准偏差马赫数小于0.01、攻角和侧滑角均小于0.1度、飞行高度小于10米；

条件B：误差变化量在5％以内；

(11)将所述步骤(2)中确定的攻角α的测压点、侧滑角β的测压点、总压源以及静压源的位置布置在真实亚音速飞行器上，重复执行步骤(8)-(9)，直到达到装订的迭代级数，所述亚音速飞行器的计算机实时解算出真实亚音速飞行器的气压高度参数，所述气压高度参数包括马赫数、攻角、侧滑角和飞行高度。

所述在飞行器体上选择攻角α的测压点以及侧滑角β的测压点的选取原则为a、b、c三个条件至少满足两个：

(a)压力差随角度变化斜率在最大斜率90％以上；

(b)压力差随角度变化线性度在90％以上；

(c)压力差变化受马赫数影响小于20％。

所述选取静压源位置的原则为a、b、c三个条件至少满足两个：

(a)静压源压力系数绝对值在0.1以内；

(b)静压源压力受弹头、弹翼及其他弹体附件影响小于10％；

(c)静压源压力随攻角、侧滑角及马赫数变化规律符合多项式规律60％以上。

所述多项式是指未知数个数不大于5个，次数不大于3次的方程。

所述步骤(4)中飞行参数修正公式选取原则为：选择所述飞行参数修正公式中干扰项的数量及形式：干扰项对待修正值影响大于30％的保留，次数不大于3次；单干扰项影响大于20％的组合形成组合干扰项，该组合干扰项中单干扰项个数不多于2个；所述干扰项包括单干扰项和组合干扰项，单干扰项包括马赫数、攻角和侧滑角；组合干扰项是指马赫数、攻角、侧滑角、马赫数多次方、攻角多次方和侧滑角多次方之间的组合。

所述步骤(4)中飞行参数修正公式包括：

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M$

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M^2$

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×\mathrm{\β}$

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×\mathrm{\β}+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M$

$\mathrm{\β}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×M+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×M^2$

$\mathrm{\β}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×\mathrm{\α}$

$\mathrm{\β}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×\mathrm{\α}+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×M$

Cp_∞＝A(1)+A(2)×β+A(3)×β²+A(4)×α+A(5)×α²

Cp_∞＝A(1)+A(2)×α+A(3)×β+A(4)×α×β+A(5)×α×M+A(6)×β×M

Cp_∞＝A(1)×α+A(2)×α²+A(3)×β+A(4)×β²+A(5)×α×β×M

Cp_∞＝(A(1)+A(2)×α)×(A(3)+A(4)×(M-A(5)))

其中，α为攻角；β为侧滑角；M为马赫数；A为系数；，A(1)是第一个干扰项的系数，A(2)是第二个干扰项的系数，A(3)是第三个干扰项的系数，A(4)是第四个干扰项的系数，A(5)是第五个干扰项的系数，ΔCp_α为攻角压力源压力系数；ΔCp_β为侧滑角压力源压力系数；Cp_∞为静压源压力系数。

本发明与现有技术相比的有益效果是：

(1)测压点及压力源直接布置在飞行器体上，减少凸出物，更利于现代飞行器的隐身设计；

(2)原有空速管技术受机体干扰严重修正量较大且修正公式中干扰项较多，本发明方法减少了干扰因素，可选择更接近目标量的一次测量值，减少修正量提高精度；

(3)原有空速管技术由于尺寸限制不利于压力测孔的布置，本发明方法可选择多位置测点形成压力源，有利于提高修正精度及减少干扰项的相关系数；

(4)本发明方法可减少飞行状态目标量回归计算的迭代次数。

附图说明

图1是本发明方法流程图；

图2是某型号气压高度试验模型示意图；其中，图2(a)为测压点轴向位置示意图；图2(b)为攻角、侧滑角测压点周向示意图；图2(c)为静压源测点周向示意图；

图3是体纵截面压力系数随攻角变化曲线；

图4是体纵截面压力系数随马赫数变化曲线；

图5是体纵截面上、下表面压力系数差值随攻角变化曲线；

图6是体横截面周向4点压力系数平均值随攻角变化曲线；

图7是体纵截面上下测点压差随攻角变化曲线；

图8是体纵截面上下测点压差随马赫数、侧滑角变化曲线；

图9是体静压源压力系数随攻角、侧滑角变化曲线；

图10是体静压源压力系数随马赫数、侧滑角变化曲线；

图11是静压回归计算绝对误差分布图；

图12是静压回归计算相对误差分布图；

图13是马赫数回归计算绝对误差分布图；

图14是马赫数回归计算相对误差分布图；

图15是攻角回归计算绝对误差分布图；

图16是侧滑角回归计算绝对误差分布图；

具体实施方式

下面结合附图对本发明的具体实施方式进行进一步的详细描述。

本发明提供了一种亚音速飞行器体上气压高度参数确定方法，是一种通过理论计算、风洞试验、公式拟合及回归算法推算提供一种飞行器飞行参数测定的新的系统设计方法，为提高飞行状态参数的测量精度，为飞控系统的准确控制提供条件。

如图1所示，本发明包括步骤如下：

(1)利用EULER算法或者N-S算法等数值计算方法对亚音速飞行器进行气动特性数值计算，得到亚音速飞行器体表面的完整压力分布；

(2)根据步骤(1)中计算得到的压力分布结果，在飞行器体上选择攻角α的测压点、侧滑角β的测压点、总压源以及静压源的位置；

所述在飞行器体上选择攻角α的测压点以及侧滑角β的测压点的选取原则为：

(2.1)压力差随角度变化斜率在最大斜率90％以上；

(2.2)压力差随角度变化线性度在90％以上；

(2.3)压力差变化受马赫数影响小于20％。

以上三个条件至少同时满足两个；

所述选取静压源位置通过如下原则进行：

(2.4)静压源压力系数绝对值在0.1以内；

(2.5)静压源压力受弹头、弹翼及其他弹体附件影响小于10％；

(2.6)静压源压力随攻角、侧滑角及马赫数变化规律符合多项式规律60％以上。

以上三个条件至少同时满足两个；

所述多项式是指未知数个数不大于5个，次数不大于3次的方程。

(3)加工所述亚音速飞行器的缩比模型，对步骤(2)中选取的攻角α的测压点、侧滑角β的测压点、总压源以及静压源进行测压风洞试验，得到攻角α角度压力点、侧滑角β角度压力点、总压源和静压源的压力系数试验数据；

(4)将步骤(3)得到的试验数据进行飞行参数修正公式拟合；

飞行参数修正公式选取原则为：选择所述飞行参数修正公式中干扰项的数量及形式：干扰项对待修正值影响大于30％的保留，次数不大于3次；单干扰项影响大于20％的组合形成组合干扰项，该组合干扰项中单干扰项个数不多于2个；所述干扰项包括单干扰项和组合干扰项，单干扰项包括马赫数、攻角和侧滑角；组合干扰项是指马赫数、攻角、侧滑角、马赫数多次方、攻角多次方和侧滑角多次方之间的组合；

本发明选取的飞行参数修正公式包括：

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M^2$

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M$

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×\mathrm{\β}$

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×\mathrm{\β}+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M$

$\mathrm{\β}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×M+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×M^2$

$\mathrm{\β}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×\mathrm{\α}$

$\mathrm{\β}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×\mathrm{\α}+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\β}}\×M$

Cp_∞＝A(1)+A(2)×β+A(3)×β²+A(4)×α+A(5)×α²

Cp_∞＝A(1)+A(2)×α+A(3)×β+A(4)×α×β+A(5)×α×M+A(6)×β×M

Cp_∞＝A(1)×α+A(2)×α²+A(3)×β+A(4)×β²+A(5)×α×β×M

Cp_∞＝(A(1)+A(2)×α)×(A(3)+A(4)×(M-A(5)))

其中，α为攻角；β为侧滑角；M为马赫数；A为系数；，A(1)是第一个干扰项的系数，A(2)是第二个干扰项的系数，A(3)是第三个干扰项的系数，A(4)是第四个干扰项的系数，A(5)是第五个干扰项的系数，ΔCp_α为攻角压力源压力系数；ΔCp_β为侧滑角压力源压力系数；Cp_∞为静压源压力系数；

上述11个公式是本发明给出的，实际拟合的时候，从里面任选其一即可，即：α的公式中选一个，β的公式中选一个，Cp_∞的公式中选一个。

(5)对步骤(4)得到的飞行参数修正公式的拟合结果进行基于最小二乘法的相关性计算，根据各干扰项的偏相关系数进一步筛选干扰项，偏相关系数小于0.2的干扰项舍弃，得到各修正公式的复相关系数；

(6)对步骤(5)得到的所述各修正公式的复相关系数进行判断，如果复相关系数的精度大于99％，则保留该复相关系数对应的修正公式，之后进入步骤(7)，否则返回步骤(4)；

(7)如果保留的修正公式大于1个，则选取复相关系数的精度最高的修正公式，并装订到所述亚音速飞行器的计算机中，之后进入步骤(8)；如果保留的修正公式为1个，则直接装订到所述亚音速飞行器的计算机中，之后进入步骤(8)；

(8)进行飞行器实物飞行试验或缩比模型风洞验证试验，得到静压源压力和总压源压力数据，解算出马赫数，并将该马赫数带入步骤(7)中确定的修正公式计算得到攻角α、侧滑角β、飞行高度以及静压源修正系数，视为一次计算值；

(9)根据步骤(8)中得到的攻角α、侧滑角β、飞行高度以及静压源修正系数的一次计算值修正静压源压力后再次解算出马赫数，并带入步骤(7)中确定的修正公式解算得到迭代后的攻角α、侧滑角β、飞行高度以及静压源修正系数的目标值即为二次迭代值；

(10)将步骤(9)中二次迭代后的目标值与一次计算值进行误差分析，如果误差满足条件A或者条件B，则回归算法固化为一次迭代；否则返回步骤(9)继续进行下一级迭代直至误差满足要求，届时完成的迭代次数即为回归算法的迭代级数，将该迭代级数装订到所述亚音速飞行器的计算机中，之后进入步骤(11)；

条件A：平均标准偏差马赫数小于0.01、攻角和侧滑角均小于0.1度、飞行高度小于10米；

条件B：收敛误差变化量在5％以内；

(11)将所述步骤(2)中确定的攻角α的测压点、侧滑角β的测压点、总压源以及静压源的位置布置在真实亚音速飞行器上，重复执行步骤(8)-(9)，直到达到装订的迭代级数，所述亚音速飞行器的计算机实时解算出真实亚音速飞行器的气压高度参数，所述气压高度参数包括马赫数、攻角、侧滑角和飞行高度。

本方法中修正拟合公式的计算精度和适用范围与试验所得数据有关，在技术要求范围内马赫数、攻角及侧滑角分布越密、试验重复次数越多，计算精度也就越高；其次，由于目标方程的多样性，正确的分析试验数据以确定几个较合理的修正公式目标方程加以比较不仅可以大幅减少计算量，对提高计算精度也有帮助。

下面以一个具体实例进一步说明本发明的工作过程。

针对一如图2(a)、(b)、(c)所示的旋成体弹身和弹翼组合体外形，图2中1点为总压测点，2、3点为攻角压力测点，4、5点为侧滑角压力测点，6、7、8、9点组成静压源；首先按步骤(1)计算弹身压力分布，图3、4为弹身表面压力分布计算结果曲线；根据此压力分布结果进行步骤(2)对数据进行分析用以选择弹身攻角α、侧滑角β测压点及总、静压源位置，α角测压点、β角测压点选取原则：1，压力差随角度变化斜率在最大斜率90％以上；2，压力差随角度变化线性度在90以上；3，压力差变化受马赫数影响小于20％。从图2、3中可初步判断在X＝0.1-0.5之间可作为选择区域，如图5中所示的曲线尖点处为压力差随攻角变化的曲线斜率最大位置，也是弹身攻角测压点的最佳轴向位置(见图2中2、3点)，由于弹身为旋成体侧滑角最佳测压点与攻角测压点在同一轴向剖面内左右布置即可(图2中4、5点)。静压源是多个测压点共同组成的静压感受系统，测点可位于同一弹身轴向剖面内也可位于不同剖面，静压源感受的压力为多个测压点压力平均值，静压源选取原则：1，静压源压力系数绝对值在0.1以内；2，静压源压力受弹头、弹翼及其他弹体附件影响小于10％；3，静压源压力随攻角、侧滑角及马赫数变化规律符合多项式规律60％以上，所述多项式是指未知数个数不大于5个，次数不大于3次的方程。根据上述规则实例模型中选取了图2(a)中的B截面为静压源截面，在此截面上按与空速管相类似的原则上下各开一对称于弹身对称平面的测压孔，4个测压点形成静压源，该静压源压力系数随攻角绝对值呈现单调变化趋势，如图6所示在弹体轴线X＝2.1-2.3区域内4点压力系数平均值均满足静压源选取原则，静压源应在该区域内选取。

根据上述测压点选择结果加工缩比模型进行地面风洞试验，即步骤(3)，获得飞行范围(M、α、β)内尽可能多的各测压点、压力源压力测量值，为提高精度可进行多次重复性试验。之后进行步骤(4)利用试验得到的压力系数数据进行飞行参数修正公式拟合，图7、8为攻角测点压力差随飞行参数变化试验曲线，可见其与攻角线性关系很好且受马赫数、侧滑角影响较小，修正公式中压力差为主要变量；图9、10为静压源压力系数随马赫数和攻角变化曲线，攻角的变化对静压影响大于30％且规律性好，由于模型为旋成体，静压源压力受侧滑角的影响与攻角一致，至此可以确定静压修正公式应包含攻角α和侧滑角β；马赫数的变化对静压也有影响，但由于静压源测压孔布置的原因，对称攻角下静压随马赫数的变化规律并不一致，且静压源压力受马赫数影响变化的幅度小于受攻角和侧滑角影响所变化的幅度，因此静压修正公式可包含马赫数M项也可不考虑马赫数的影响。找出对攻角压差及静压差值CP∞影响较大的参数，可作出多个修正公式的目标方程，例如：

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A{\left(3\right)}_{10}\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M+A\left(4\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M^2$

$\mathrm{\α}=A\left(1\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}+A\left(2\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}^2+A\left(3\right)\×\mathrm{\ΔC}{p}_{\mathrm{\α}}\×M$

Cp_∞＝A(1)+A(2)×β+A(3)×β²+A(4)×α+A(5)×α²

Cp_∞＝A(1)+A(2)×α+A(3)×β+A(4)×α×β+A(5)×α×M+A(6)×β×M

Cp_∞＝A(1)×α+A(2)×α²+A(3)×β+A(4)×β²+A(5)×α×β×M

进行步骤(5)、(6)利用基于最小二乘法的多元回归法确定各组方程的系数、各系数的偏相关系数及各修正公式的复相关系数，根据偏相关系数进行干扰项筛选，根据复相关系数进行修正公式筛选，选取其中一组精度最好(复相关系数最高，大于99％)的方程为最终修正公式。

以第二个公式为例，其拟合结果如下：

回归系数

A(1)＝15.170530

A(2)＝-0.418769

A(3)＝-3.011092

偏差平方和Q，平均标准偏差S，复相关系数R

Q＝19.629020 S＝.172456 R＝.999628

偏相关系数V

V(1)＝.999702

V(2)＝.573959

V(3)＝.974956

该公式复相关系数0.999628表明其具有很高的精度，可以作为最终修正公式，变量1、3偏相关系数均在0.97以上表明该变量对数据修正贡献度很大，选择合理；变量2偏相关系数低于20％说明该变量可以不带入修正公式中。根据以上原则最终可确定攻角、侧滑角和静压修正公式并装载到弹载计算机中。

对缩比模型进行风洞验证试验，可回代入试验压力值进行计算，计算结果与真值比较以确定计算精度。回代计算分为两步：一次计算和多次迭代计算。如步骤(8)所示首先采用测得的总压和静压直接算出来流马赫数M，之后利用步骤(7)中确定的修正公式计算出攻角α、侧滑角β以及静压修正系数，至此一次计算结束；二次迭代计算时利用一次计算中求出的α、β和静压修正系数代入静压修正公式修正静压值求出马赫数，并再次带入步骤(7)中确定的修正公式解算得到攻角α、侧滑角β以及静压修正系数，将其误差与上级计算误差比较分析如满足精度要求则回代计算完成，如不满足要求则利用新求出的静压计算马赫数M，之后解算攻角α、侧滑角β再次代入静压修正公式求出静压并与上级计算精度比较直到误差满足要求或误差收敛为止，届时完成的迭代次数即为回归算法的迭代级数，固化后即可作为亚音速气压高度系统修正算法装载到弹载计算机中。

本实例一次计算和迭代计算各计算项平均标准偏差如下：

图11至图16为误差图谱，其中包括静压P、马赫数M的绝对误差和相对误差，攻角α、侧滑角β的绝对误差。计算结果表明，二次迭代计算结果已经具有较高的精度，其相对误差最高不超过0.5％，迭代次数固化为二次。但需要注意的是，并不是迭代计算次数越多计算结果精度越高，如本实例二次迭代计算结果中攻角α的计算精度就低于一次计算，这与多方面因素有关，如静压源开孔位置、修正公式选择得是否合理等，具体回代次数的选择应根据实际情况而定。

将步骤(2)中确定的攻角α的测压点、侧滑角β的测压点、总压源以及静压源的位置布置在真实弹体上，重复执行步骤(8)-(9)，直到达到装载的迭代级数，弹载计算机可实时解算出真实亚音速飞行器的气压高度参数，所述气压高度参数包括马赫数、攻角、侧滑角和飞行高度。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (6)

1.一种亚音速飞行器体上气压高度参数确定方法，其特征在于步骤如下：

(1)对亚音速飞行器进行气动特性数值计算，得到亚音速飞行器体表面的压力分布；

(2)根据步骤(1)中计算得到的压力分布结果，在飞行器体上选择攻角α的测压点、侧滑角β的测压点、总压源以及静压源的位置；

(3)加工所述亚音速飞行器的缩比模型，对步骤(2)中选取的攻角α的测压点、侧滑角β的测压点、总压源以及静压源进行测压风洞试验，得到攻角α角度压力点、侧滑角β角度压力点、总压源和静压源的压力系数试验数据；

(4)将步骤(3)得到的试验数据进行飞行参数修正公式拟合；

(5)对步骤(4)得到的飞行参数修正公式的拟合结果进行相关性计算，根据各干扰项的偏相关系数进一步筛选干扰项，偏相关系数小于0.2的干扰项舍弃，得到各修正公式的复相关系数；

(6)对步骤(5)得到的所述各修正公式的复相关系数进行判断，如果复相关系数的精度大于99％，则保留该复相关系数对应的修正公式，之后进入步骤(7)，否则返回步骤(4)；

(7)如果保留的修正公式大于1个，则选取复相关系数的精度最高的修正公式，并装订到所述亚音速飞行器的计算机中，之后进入步骤(8)；如果保留的修正公式为1个，则直接装订到所述亚音速飞行器的计算机中，之后进入步骤(8)；

(8)进行飞行试验或风洞验证试验，得到静压源压力和总压源压力数据，解算出马赫数，并将该马赫数带入步骤(7)中确定的修正公式计算得到攻角α、侧滑角β、飞行高度以及静压源修正系数，视为一次计算值；

(9)根据步骤(8)中得到的攻角α、侧滑角β、飞行高度以及静压源修正系数的一次计算值再次解算出马赫数，并带入步骤(7)中确定的修正公式解算得到迭代后的攻角α、侧滑角β、飞行高度以及静压源修正系数的目标值即为二次迭代值；

(10)将步骤(9)中二次迭代后的目标值与一次计算值进行误差分析，如果误差满足条件A或者条件B，则回归算法固化为一次迭代；否则返回步骤(9)继续进行下一级迭代直至误差满足要求，届时完成的迭代次数即为回归算法的迭代级数，将该迭代级数装订到所述亚音速飞行器的计算机中，之后进入步骤(11)；

条件A：平均标准偏差马赫数小于0.01、攻角和侧滑角均小于0.1度、飞行高度小于10米；

条件B：误差变化量在5％以内；

(11)将所述步骤(2)中确定的攻角α的测压点、侧滑角β的测压点、总压源以及静压源的位置布置在真实亚音速飞行器上，重复执行步骤(8)-(9)，直到达到装订的迭代级数，所述亚音速飞行器的计算机实时解算出真实亚音速飞行器的气压高度参数，所述气压高度参数包括马赫数、攻角、侧滑角和飞行高度。

2.根据权利要求1所述的一种亚音速飞行器体上气压高度参数确定方法，其特征在于：所述在飞行器体上选择攻角α的测压点以及侧滑角β的测压点的选取原则为a、b、c三个条件至少满足两个：

(a)压力差随角度变化斜率在最大斜率90％以上；

(b)压力差随角度变化线性度在90％以上；

(c)压力差变化受马赫数影响小于20％。

3.根据权利要求1所述的一种亚音速飞行器体上气压高度参数确定方法，其特征在于：所述选取静压源位置的原则为a、b、c三个条件至少满足两个：

(a)静压源压力系数绝对值在0.1以内；

(b)静压源压力受弹头、弹翼及其他弹体附件影响小于10％；

(c)静压源压力随攻角、侧滑角及马赫数变化规律符合多项式规律60％以上。

4.根据权利要求3所述的一种亚音速飞行器体上气压高度参数确定方法，其特征在于：所述多项式是指未知数个数不大于5个，次数不大于3次的方程。

5.根据权利要求1所述的一种亚音速飞行器体上气压高度参数确定方法，其特征在于：所述步骤(4)中飞行参数修正公式选取原则为：选择所述飞行参数修正公式中干扰项的数量及形式：干扰项对待修正值影响大于30％的保留，次数不大于3次；单干扰项影响大于20％的组合形成组合干扰项，该组合干扰项中单干扰项个数不多于2个；所述干扰项包括单干扰项和组合干扰项，单干扰项包括马赫数、攻角和侧滑角；组合干扰项是指马赫数、攻角、侧滑角、马赫数多次方、攻角多次方和侧滑角多次方之间的组合。

6.根据权利要求5所述的一种亚音速飞行器体上气压高度参数确定方法，其特征在于：所述步骤(4)中飞行参数修正公式包括：

$\mathrm{\α}=A\left(1\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}+A\left(2\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}^2+A\left(3\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}\×M+A\left(4\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}\×M^2$

$\mathrm{\α}=A\left(1\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}+A\left(2\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}^2+A\left(3\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}\×M$

$\mathrm{\α}=A\left(1\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}+A\left(2\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}^2+A\left(3\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}\×\mathrm{\β}$

$\mathrm{\α}=A\left(1\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}+A\left(2\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}^2+A\left(3\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}\×\mathrm{\β}+A\left(4\right)\×{\mathrm{\ΔCp}}_{\mathrm{\α}}\×M$

$\mathrm{\β}=A\left(1\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}+A\left(2\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}^2+A\left(3\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}\×M+A\left(4\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}\×M^2$

$\mathrm{\β}=A\left(1\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}+A\left(2\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}^2+A\left(3\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}\×\mathrm{\α}$

$\mathrm{\β}=A\left(1\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}+A\left(2\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}^2+A\left(3\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}\×\mathrm{\α}+A\left(4\right)\×{\mathrm{\ΔCp}}_{\mathrm{\β}}\×M$

Cp_∞＝A(1)+A(2)×β+A(3)×β²+A(4)×α+A(5)×α²

Cp_∞＝A(1)+A(2)×α+A(3)×β+A(4)×α×β+A(5)×α×M+A(6)×β×M

Cp_∞＝A(1)×α+A(2)×α²+A(3)×β+A(4)×β²+A(5)×α×β×M

Cp_∞＝(A(1)+A(2)×α)×(A(3)+A(4)×(M-A(5)))

其中，α为攻角；β为侧滑角；M为马赫数；A为系数；A(1)是第一个干扰项的系数，A(2)是第二个干扰项的系数，A(3)是第三个干扰项的系数，A(4)是第四个干扰项的系数，A(5)是第五个干扰项的系数，ΔCp_α为攻角压力源压力系数；ΔCp_β为侧滑角压力源压力系数；Cp_∞为静压源压力系数。