一种高阶大气紊流成型滤波器的设计方法及仿真方法与流程

本发明涉及飞行仿真领域，尤其涉及一种高阶大气紊流成型滤波器的设计方法及仿真方法。

背景技术：

各类飞行器在大气层内飞行时经常受到大气紊流的影响，大气紊流的干扰会导致飞行器飞行性能下降、操纵困难，甚至导致驾驶员诱发振荡，危及飞行安全。因此，发展一种高效精确的大气紊流仿真方法对飞行器飞行控制系统设计、飞行品质评估、结构载荷和疲劳预估以及驾驶员培训都有很重要的意义。

Dryden和Von Karman模型是航空航天领域常用的两种大气紊流模型。基于Taylor冻结场假设，Dryden模型由于其有理形式、可以通过简单的共轭分解获得大气紊流传递函数等优点而在各类飞行器飞行仿真中获得广泛的应用。然而Dryden模型的缺点在于其精度较低，Dryden谱函数的有效无量纲频率范围仅为0～10rad，导致基于该模型的飞行器大气紊流仿真在高空低速飞行状态的精度降低，使相关的飞行品质分析可能出现过大的误差。与Dryden模型相比，Von Karman模型在高频段满足Kolmogorov提出的-5/3衰减规律，与实验数据的吻合良好，具有极高的精度。因此，发展一种基于Von Karman模型的大气紊流仿真方法对飞行器飞行仿真和飞行品质分析是有必要的，尤其当飞行器处于高空低速的飞行状态时。

由于其无理形式，将Von Karman模型直接用于飞行器大气紊流仿真并不容易。其中一种思路是首先采用数值拟合的方式得到Von Karman谱函数的有理近似式，然后基于Taylor冻结场假设对有理近似谱函数进行共轭分解，获得有理形式的大气紊流传递函数。然而现有研究得到的大气紊流传递函数的有效无量纲频率范围不高，最大不超过50rad。另外一种思路是采用数学方法直接对Von Karman模型处理实现大气紊流仿真，这类方法包括傅里叶变换/滤波方法，卷积积分法和正弦函数叠加法。然而该种类型方法往往比较复杂，直接用于飞行器飞行仿真较为困难。

技术实现要素：

本发明提供一种高阶大气紊流成型滤波器的设计方法及仿真方法，具体包括一种高阶大气紊流传递函数、一种大气紊流离散递推数值算法及其计算机仿真实现步骤，以提升飞行器在高空低速飞行状态的大气紊流仿真和飞行品质分析精度。

一种高阶大气紊流成型滤波器的设计方法，包括以下步骤：

S1、采用最小二乘法对角频率形式的Von Karman谱函数进行数值拟合得到其高阶有理近似式，然后对有理形式的大气紊流谱函数进行共轭分解，最终得到有理形式的高阶大气紊流传递函数；

S2、采用Tustin变换将共轭分解后的大气紊流传递函数离散得到初步的离散递推数值算法，分析初步算法的频域响应特性，根据分析结果对初步算法进行修正，最终得到高精度的大气紊流离散递推数值算法；

S3、根据S2的高精度大气紊流离散递推数值算法设计计算机仿真实现步骤。

进一步的，所述高阶大气紊流传递函数为：

式中：H_u(ω)、H_v(ω)、H_w(ω)分别为纵向、横向和垂向的大气紊流速度传递函数，σ_u、σ_v、σ_w分别为纵向、横向和垂向大气紊流速度的强度，π为圆周率，γ_u,v,w＝V/L_u,v,w，其中V为飞行器飞行速度，L_u、L_v、L_w分别为纵向、横向和垂向大气紊流速度谱函数的特征长度。

进一步的，所述高效高精度大气紊流速度的离散递推数值算法包括纵向、横向和垂向三个方面，其中，纵向大气紊流速度的离散递推数值算法为：

式中：x_i代表离散Gauss白噪声随机信号，y_i代表离散纵向大气紊流速度，f_j(j＝1,…,7)为递推系数且其具体形式见表1；

横向和垂向大气紊流速度的离散递推数值算法为：

式中：y_i代表离散横向或垂向大气紊流速度，g_j(j＝1,…,7)为递推系数且其具体形式见表2。

本发明还提供了高阶大气紊流成型滤波器的仿真方法，包括以下步骤：

SS1、根据飞行器飞行状态，确定飞行器的飞行高度h、飞行速度V、垂向大气紊流速度的强度σ_w；

SS2、根据飞行品质规范MIL-F-8785C的规定计算三个方向大气紊流速度谱函数的特征长度L_u,L_v,L_w和纵向大气紊流速度的强度σ_u、横向大气紊流速度的强度σ_v：

式中：f_uv＝0.177+0.0027h；

SS3、根据离散取样周期Δt及表1和表2计算纵向大气紊流速度离散递推系数f_i(i＝1,…,7)和横向、垂向大气紊流速度离散递推系数g_i(i＝1,…,9)；

SS4、获得均值为0方差为1的Gauss白噪声信号；

SS5、根据式(15)和(15)，计算纵向、横向和垂向大气紊流速度。

本发明有益效果是：

1)与现有大气紊流传递函数相比，本发明提供的一种高阶大气紊流传递函数具有更高的有效无量纲频率范围，因此基于该传递函数的大气紊流数值算法能够显著提升飞行器在高空低速飞行状态的大气紊流仿真精度；

2)本发明提供的高精度大气紊流离散递推数值算法具有良好的数值稳定性和较高的计算效率，更适合飞行器大气紊流飞行仿真和相关飞行品质分析的应用；

3)本发明提供的仿真方法步骤简单易实现，可方便地用于各种类型飞行器的大气紊流飞行仿真和飞行品质分析。

附图说明

图1为本发明提供的一种高阶大气紊流成型滤波器的原理图；

图2为本发明拟合近似谱函数与Dryden和Von Karman谱函数的对比图，其中图2.1为纵向大气紊流速度谱函数的对比图，图2.2横向大气紊流速度谱函数的对比图；

图3为大气紊流强度是1m/s时，飞行器飞行速度V与大气紊流速度谱函数特征长度L的比值γ＝1飞行状态下仿真的大气紊流速度功率谱密度与Von Karman谱函数的对比图，其中图3.1为纵向大气紊流速度功率谱密度的对比图，图3.2为横向大气紊流速度功率谱密度的对比图；

图4为大气紊流强度是1m/s时,飞行器飞行速度V与大气紊流速度谱函数特征长度L的比值γ＝0.2飞行状态下仿真的大气紊流速度功率谱密度与Von Karman谱函数的对比图，其中图4.1为纵向大气紊流速度功率谱密度的对比图，图4.2为横向大气紊流速度功率谱密度的对比图。

具体实施方式

下面结合附图，对本发明提供的一种高阶大气紊流成型滤波器的设计方法及仿真方法进行详细说明。

如图1所示，所述一种高阶大气紊流成型滤波器的原理为：

通过对无理形式的Von Karman谱函数进行有理近似获得有理形式的大气紊流谱函数，然后对有理近似谱函数进行共轭分解得到有理形式的大气紊流传递函数，最后以Gauss分布白噪声信号x通过大气紊流传递函数生成具有Von Karman谱特征的大气紊流速度y，即：

式中：Φ_X(ω)为输入Gauss白噪声信号的功率谱密度，Φ_Y(ω)为输出大气紊流速度的功率谱密度。

假如输入信号是理想白噪声，且其功率谱密度恒等于1，即

Φ_X(ω)＝1 (2)

则大气紊流速度的功率谱密度为

因此可以通过将功率谱密度为1的理想Gauss白噪声信号输入具有Von Karman谱特征的传递函数进行大气紊流仿真。

所述一种高阶大气紊流传递函数通过以下推导获得：

空间频率形式的Von Karman谱函数为：

式中：Φ_u(Ω),Φ_v(Ω),Φ_w(Ω)分别为纵向、横向和垂向大气紊流速度的谱函数，Ω为空间频率，π为圆周率，a＝1.339，σ_u,σ_v,σ_w分别为纵向、横向和垂向大气紊流速度的强度，L_u,L_v,L_w分别为纵向、横向和垂向大气紊流速度谱函数的特征长度。

基于Taylor冻结场假设，角频率和空间频率以及角频率形式和空间频率形式的大气紊流谱函数存在如下关系：

ω＝Ω·V (7)

Φ_u,v,w(ω)＝Φ_u,v,w(Ω)/V (8)

利用式(7)和(8)将空间形式的Von Karman谱函数转换为角频率形式，并无量纲化可得：

式中：为无量纲角频率，其中角频率无量纲参数γ_u,v,w＝V/L_u,v,w为飞行器前飞速度V与大气紊流速度谱函数特征长度L_u,v,w的比值。

采用最小二乘法对式(9)～(11)进行拟合，拟合所选取的无量纲频率范围为0～100rad，得到有理近似大气紊流谱函数，然后通过共轭分解得到三个方向大气紊流传递函数分别为：

式中：H_u(s),H_v(s),H_w(s)分别为纵向、横向和垂向大气紊流速度的传递函数。

所述大气紊流离散递推数值算法通过以下推导获得：

采用Tustin变换将s域的大气紊流传递函数转换到z域，即

：

将式(14)代入式(12)离散后得到纵向大气紊流速度的离散递推数值算法为：

式中：x_i为离散Gauss白噪声信号，y_i表示离散的纵向大气紊流速度，f_j(j＝1,…,7)为递推系数。

将式(14)代入式(13)离散后得到横向和垂向大气紊流速度的离散递推数值算法为：

式中：y_i表示离散的横向或垂向大气紊流速度，g_j(j＝1,…,7)为递推系数。

由式(15)或(16)的离散递推数值算法进行计算得到的离散大气紊流速度的功率谱密度为：

Φ_DY(ω_D)＝|H_D(ω_D)|²Φ_DX(ω_D) (17)

式中：Φ_DX(ω_D),Φ_DY(ω_D)分别为输入离散Gauss白噪声信号和输出离散大气紊流速度的功率谱密度，H_D(ω_D)表示离散递推数值算法的频域响应特性，ω_D为离散角频率。

常用窄带高斯随机信号的强度为1，且其功率均匀分布于-ω_Dmax～ω_Dmax区间内。其中，ω_Dmax＝π/Δt为离散取样的Nyquist频率，Δt为离散取样周期。因此，其单侧功率谱密度为：

同时，在远小于Nyquist频率的频段，由Tustin变换得到的离散递推数值算法的频域特性H_D(ω_D)与离散前传递函数的频域特性H(ω)吻合良好，因此：

|H_D(ω_D)|²≈|H(ω_D)|²＝Φ_Y(ω_D) (19)

将式(18)、(19)代入式(17)可得：

由上式可以看出，由于常用的窄带Gauss白噪声信号强度为1，功率谱密度为Δt/π，导致由式(15)和(16)离散递推数值算法计算的大气紊流速度的功率谱密度缩小了Δt/π倍。因此将输入Gauss随机信号的系数放大倍再进行计算即可得到具有理想谱特征的大气紊流速度。修正后式(15)和(16)的递推系数分别见表1和表2。

表1

表2

本发明还提供了高阶大气紊流成型滤波器的仿真方法，包括如下步骤：

SS1、根据飞行器飞行状态，确定飞行器的飞行高度h、飞行速度V、垂向大气紊流速度的强度σ_w；

SS2、根据飞行品质规范MIL-F-8785C的规定计算三个方向大气紊流速度谱函数的特征长度L_u,L_v,L_w和纵向大气紊流速度的强度σ_u、横向大气紊流速度的强度σ_v：

式中：f_uv＝0.177+0.0027h；

SS3、根据离散取样周期Δt及表1和表2计算纵向大气紊流速度离散递推系数f_i(i＝1,…,7)和横向、垂向大气紊流速度离散递推系数g_i(i＝1,…,9)；

SS4、获得均值为0方差为1的Gauss白噪声信号；

SS5、根据式(15)和(16)计算纵向、横向和垂向大气紊流速度。

图2所示为本发明拟合近似谱函数与Dryden和Von Karman谱函数的对比图。从图中可以看出，与Dryden谱函数相比，本发明拟合的近似大气紊流谱函数在更高的有效无量纲频率范围内与Von Karman谱函数吻合良好。有效无量纲频率范围的定义为满足下式的无量纲角频率：

式中：为理论Von Karman谱函数，为其他近似大气紊流谱函数。

根据上式的定义，Dryden谱函数的有效频率范围为0～10rad，而本发明近似大气紊流谱函数的有效无量纲频率范围高达0～400rad。

使用本发明提供的一种高阶大气紊流成型滤波器的仿真方法进行FORTRAN编程仿真，并分别对飞行器飞行速度V与大气紊流速度谱函数特征长度L的比值γ＝1和0.2两种飞行状态进行仿真获得大气紊流速度时间历程，然后计算大气紊流速度时间历程的功率谱密度并与Von Karman谱函数进行对比，如图3和图4所示。其中两种飞行状态的大气紊流强度均为1m/s。从图中可以看出，即使对于的飞行状态，本发明提供的一种高阶大气紊流成型滤波器仿真得到的大气紊流速度的功率谱密度仍然在飞行品质相关频率范围内(1～10rad/s)与Von Karman谱函数吻合良好。

本发明方法的具体应用途径很多，以上所述仅是本发明的优选实施方式。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些改进也应视为本发明的保护范围。