基于飞行数据的扰动风定制化建模和大气数据估计方法

1.本技术属于民航安全技术与飞行数据应用领域，具体涉及基于飞行数据的扰动风定制化建模和大气数据估计方法。


背景技术：

2.大气环境中除了缓慢变化的盛行风外，还有各种离散突风和连续湍流等扰动风。扰动风不仅严重影响飞机飞行安全，同时也造成空速、迎角和侧滑角的测量误差。民航飞机的飞行数据实时记录了多种飞行参数，在飞行异常事件和事故分析过程中起到关键作用。然而，在大气扰动诱发的事故分析中，受扰动风的影响，飞机大气数据系统采集的空速、迎角、侧滑角存在误差。
3.解决这一问题的有效途径是结合飞行动力学方程、以飞行数据提供的飞行参数为量测输入来构建空速、迎角和侧滑角的状态估计滤波算法，同时建立扰动风模型并集成在滤波算法内。然而，实际的复杂扰动风场不仅难以采集，且与扰动风理论模型存在显著差异。若在滤波算法中采用通用的扰动风经验模型，必然影响空速、迎角和侧滑角的估计精度。
4.基于飞行数据的飞行参数状态估计与一般的状态估计问题又有显著差异。表现在：在扰动风影响下，飞机动力学响应具有强非线性特征，采用传统的基于一阶线性化的扩展卡尔曼滤波算法将导致估计精度差；飞行数据记录参数具有不确定性，难以确定滤波算法中的过程噪声和量测噪声模型。


技术实现要素：

5.本技术提出了基于飞行数据的扰动风定制化建模和大气数据估计方法，基于飞行数据，能够实现在大气扰动飞行中空速、迎角和侧滑角的准确估计，克服以往的状态估计方法缺少定制化的精确扰动风模型、动力学模型需进行一阶线化、噪声模型不确定而导致的估计误差等缺点，并利用平滑估计实现提高空速、迎角、侧滑角估计精度的发明目的。
6.为实现上述目的，本技术提供了如下方案：
7.基于飞行数据的扰动风定制化建模和大气数据估计方法，所述大气数据估计包括空速、迎角、侧滑角估计，包括如下步骤：
8.s1.基于选取的飞行记录数据，获得扰动风时间序列，基于所述扰动风时间序列，获得满足平稳性和正态性的湍流序列；
9.s2.基于所述湍流序列，建立所述湍流序列中每个分量的高斯-马尔可夫模型，完成扰动风定制化建模；
10.s3.基于所述湍流序列中每个分量的高斯-马尔可夫模型，并根据所述飞行力学原理，建立包含扰动风影响的空速、迎角和侧滑角的状态方程和量测方程；
11.s4.基于所述状态方程和所述量测方程，建立自适应平方根无迹卡尔曼平滑器，并基于所述自适应平方根无迹卡尔曼平滑器，得到任意时刻的空速、迎角、侧滑角的最优平滑
估计，完成大气数据估计。
12.可选的，所述s1中，获取所述湍流序列的方法包括：
13.根据空间扰动风与地速、空速的矢量三角形关系，得到飞机沿飞行轨迹所经历的所述扰动风三轴分量，并获得扰动风时间序列；
14.对所述扰动风时间序列进行平稳性检验和正态性检验，若检验通过，则将所述扰动风时间序列转换为湍流序列；
15.其中，x轴湍流序列为w
tx
＝[w
tx1
，w
tx2
，...，w
txn
]
t
，y轴湍流序列为w
ty
＝[w
ty1
，w
ty2
，...，w
tyn
]
t
，z轴湍流序列为w
tz
＝[w
tz1
，w
tz2
，...，w
tzn
]
t
。
[0016]
可选的，所述平稳性检验采用单位根检验方法；
[0017]
所述正态性检验采用分位点q-q图检验方法。
[0018]
可选的，如果所述扰动风时间序列不能通过所述平稳性检验和所述正态性检验，则建立高斯过程回归模型，并设定长度尺度进行模型预测，得到预测值序列；
[0019]
使用所述扰动风时间序列减去所述预测值序列得到残差序列；
[0020]
对所述残差序列进行所述平稳性检验和所述正态性检验，若检验通过，获得所述湍流序列。
[0021]
可选的，对所述残差序列进行所述平稳性检验和所述正态性检验不通过，则将所述长度尺度乘以2，获得新的预测值序列，重新获取残差序列并进行所述平稳性检验和所述正态性检验，重复本步骤，直到通过所述平稳性检验和所述正态性检验为止，得到所述湍流序列。
[0022]
可选的，所述s2中，所述湍流序列中x轴的高斯-马尔可夫模型为
[0023][0024]
所述湍流序列中y轴的高斯-马尔可夫模型为
[0025][0026]
所述湍流序列中z轴的高斯-马尔可夫模型为
[0027][0028]
可选的，所述s3中，建立所述状态方程和所述量测方程的方法包括：
[0029]
根据飞行力学原理和所述湍流序列中每个分量的高斯-马尔可夫模型，建立以x＝[v
t
，α，β，w
tx
，w
ty
，w
tz
]
t
为状态量的所述状态方程，所述状态方程为
[0030][0031]
[0032][0033][0034]
基于所述飞行力学原理，建立以y＝[v
gx
，v
gy
，v
gz
]
t
为量测参数的所述量测方程；
[0035]
所述量测方程为
[0036][0037]
可选的，所述s4包括：
[0038]
对所述自适应平方根无迹卡尔曼平滑器的前向滤波初始化，所述初始化方法为
[0039][0040]
其中，是初始状态估计协方差，通过cholesky分解，获得所述初始状态估计协方差的平方根s0，将所述初始状态估计协方差的平方根s0作为初始过程噪声协方差的平方根是初始量测估计协方差，通过cholesky分解，获得初始量测噪声协方差的平方根
[0041]
对所述自适应平方根无迹卡尔曼平滑器的前向时间更新，所述时间更新为
[0042][0043]
对所述自适应平方根无迹卡尔曼平滑器的前向量测更新及噪声协方差阵的自适应更新，包括量测噪声协方差阵平方根的自适应更新、量测更新和过程噪声协方差阵平方根的自适应更新；
[0044]
所述量测噪声协方差阵平方根的自适应更新为
[0045][0046]
所述量测更新为
[0047][0048]
所述过程噪声协方差阵平方根的自适应更新为：
[0049][0050]
使用更新后的所述自适应平方根无迹卡尔曼平滑器得到任意时刻的空速、迎角、侧滑角的最优平滑估计，包括：
[0051]
选取预设的时间间隔k，并初始化所述自适应平方根无迹卡尔曼平滑器
[0052][0053]
通过定义状态组合，得到后验估计及相应的协方差阵，根据贝叶斯方法，获得获得
k时刻状态估计及其协方差：
[0054]
基于所述k时刻状态估计及其协方差，根据贝叶斯方法，获得组合状态的平滑估计及其协方差，并得到k时刻的平滑估计及其协方差；
[0055]
对所述k时刻的平滑估计及其协方差进行qr分解和cholesky秩1更新，得到k时刻的平滑估计及其协方差阵的平方根为
[0056][0057]
本技术的有益效果为：
[0058]
本技术公开了基于飞行数据的扰动风定制化建模和大气数据估计方法，通过飞行数据推算获得飞机沿飞行轨迹的扰动风，基于扰动风时间序列建立扰动风的定制化模型，该模型具有通用性，能够准确还原飞机飞行过程中经历的扰动风，从而有助于提升空速、迎角和侧滑角估计精度；利用平方根无迹卡尔曼滤波提高滤波精度；通过过程噪声和量测噪声协方差阵的自适应调整，有效降低飞行参数噪声不确定性对滤波精度和稳定性的影响；考虑到飞行数据离线事后分析的特点，对有限长飞行数据进行平滑估计以进一步提高估计精度。本技术实现了基于飞行数据准确估计扰动风影响下的飞机空速、迎角、侧滑角的目的，从而辅助进行大气扰动诱发的飞行事故分析。
附图说明
[0059]
为了更清楚地说明本技术的技术方案，下面对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本技术的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0060]
图1为本技术实施例的基于飞行数据的扰动风定制化建模和大气数据估计方法流程示意图；
[0061]
图2为本技术实施例中获取x轴湍流序列的流程示意图。
具体实施方式
[0062]
下面将结合本技术实施例中的附图，对本技术实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本技术一部分实施例，而不是全部的实施例。基于本技术中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本技术保护的范围。
[0063]
为使本技术的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本技术作进一步详细的说明。
[0064]
如图1所示，为本技术实施例的基于飞行数据的扰动风定制化建模和大气数据估计方法流程示意图，主要包括如下步骤：
[0065]
s1.基于选取的飞行记录数据，获得扰动风时间序列，基于所述扰动风时间序列，
获得满足平稳性和正态性的湍流序列。
[0066]
在本实施例中，首先根据空间扰动风与地速、空速的矢量三角形关系，可得到飞机沿飞行轨迹所经历的扰动风三轴分量w，表达式如式(1)所示：
[0067][0068]
其中w＝[w
x
，wy，wz]
t
为扰动风三轴分量，vg＝[v
gx
，v
gy
，v
gz
]
t
为飞行数据记录的地速值，为机体系到地面系的转移矩阵，其中的φ为飞机滚转角，θ为俯仰角，ψ为偏航角，c和s分别为cos和sin的缩写(下同)，为气流系到机体系的转移矩阵，其中α为迎角，β为侧滑角。为空速矩阵。从而根据当前一段时间内采集到的有限的飞行数据记录的地速vg、空速v
t
、滚转角φ、俯仰角θ、偏航角ψ、迎角α，侧滑角β，获得扰动风时间序列
[0069]
在本实施例中，以扰动风的x轴分量序列w
x
＝[w
x1
，w
x2
，...，w
xn
]
t
为例，如图2所示，首先对其进行平稳性和正态性检验，在本实施例中，序列的平稳性采用单位根检验方法，序列的正态性采用分位点q-q图检验方法。如果检验通过，则认为是满足平稳性和正态性的湍流序列。否则建立高斯过程回归模型，用扰动风序列w
x
＝[w
x1
，w
x2
，...，w
xn
]
t
进行模型训练。此处设定长度尺度n进行模型预测，获得预测值序列将原始扰动风序列减去预测值序列后形成残差序列，再次进行平稳性和正态性检验。若通过检验，则获得满足平稳性和正态性的湍流序列w
tx
＝[w
tx1
，w
tx2
，...，w
txn
]
t
；否则，令n＝n
×
2，再次获得预测值序列，执行后续过程，直至通过平稳性和正态性检验、得到湍流序列为止。
[0070]
在本实施例中，获得满足平稳性和正态性的湍流序列为
[0071]
s2.基于所述湍流序列，建立所述湍流序列中每个分量的高斯-马尔可夫模型。
[0072]
在本实施例中，仍以湍流序列的x轴为例，高斯-马尔可夫模型可表示为
[0073][0074]
其中，ts为序列采样周期，由飞行参数的记录频率确定。τ
x
是高斯-马尔可夫过程的时间相关系数。
[0075]
设为驱动白噪声的强度，将序列样本的自相关系数函数(1)代入1阶自回归模型的yule-walker方程：
[0076][0077]
其中为序列w
tx
的自相关函数。由式(3)计算获得时间相关系数τ
x
。进一步由
[0078][0079]
估计出模型白噪声w
x
的强度其中，为序列w
tx
的自协方差函数。从而建立湍流分量w
tx
序列的高斯-马尔可夫模型。
[0080]
同理，本实施例建立的湍流分量w
ty
序列和w
tz
序列的高斯-马尔可夫模型分别为：
[0081][0082][0083]
s3.基于所述湍流序列中每个分量的高斯-马尔可夫模型，并根据所述飞行力学原理，建立包含扰动风影响的空速、迎角和侧滑角的状态方程和量测方程。
[0084]
建立x＝[v
t
，α，β，w
tx
，w
ty
，w
tz
]
t
为状态参数的状态方程的过程：
[0085]
在本实施例中，将式(1)作变换，可获得在扰动气流影响下飞机相对地面系的地速矢量vg为
[0086][0087]
对上式求导，并代入牛顿第二定律可得
[0088][0089]
对式(8)进行转换，有由于则代入式(8)的右侧并经整理可得
[0090][0091]
从而有
[0092][0093]
上式中，为斜对称阵，[p，q，r]
t
为三轴角速度，由飞行数据记录值提供。由于飞机在空中的受力f包括了重力、气动力和发动机推力，难以准确求解。本实施例根据飞行数据记录的加速度值a＝[a
x
，ay，az]
t
取代受力的计算，从而将式(10)转化为
[0094][0095]
式(11)作为含空速、迎角和侧滑角的状态方程。
[0096]
再结合(2)、(5)、(6)式提供的湍流序列模型，建立以x＝[v
t
，α，β，w
tx
，w
ty
，w
tz
]
t
为状态量的状态方程。
[0097]
建立以y＝[v
gx
，v
gy
，v
gz
]
t
为量测参数的量测方程的过程：
[0098]
将式(1)直接转换为
[0099][0100]
从而建立以y＝[v
gx
，v
gy
，v
gz
]
t
为量测参数的量测方程。
[0101]
基于式(11)、(12)获得离散化的非线性系统状态方程和量测方程为，
[0102][0103]
其中xk、yk分别为k时刻系统的状态向量和量测向量；wk为过程噪声，vk为量测噪
声。
[0104]
s4.基于所述状态方程和所述量测方程，建立自适应平方根无迹卡尔曼平滑器，并基于所述自适应平方根无迹卡尔曼平滑器，得到任意时刻的空速、迎角、侧滑角的最优平滑估计。
[0105]
与传统的基于一阶线性化近似的扩展卡尔曼滤波算法相比，无迹卡尔曼滤波算法无需对非线性系统线性化，具有更高的估计精度。进一步采用协方差阵的平方根进行算法递推，可以避免协方差阵递推过程中的奇异问题。本技术在基本平方根无迹卡尔曼滤波算法的基础上，考虑飞行数据记录噪声的不确定性，对过程噪声和量测噪声的协方差阵进行自适应更新；利用有限长飞行数据的全部观测值，获得任意时刻空速、迎角、侧滑角的最优平滑估计。
[0106]
具体的，本步骤包括了四个小步骤：
[0107]
s4.1：前向自适应平方根无迹卡尔曼滤波初始化。
[0108]
在本实施例中，初步给定过程噪声和量测噪声的统计特性为
[0109][0110]
基本平方根无迹卡尔曼滤波算法的初始化方法为：
[0111][0112]
其中，是初始状态估计协方差，通过cholesky分解，获得所述初始状态估计协方差的平方根s0，将所述初始状态估计协方差的平方根s0作为初始过程噪声协方差的平方根是初始量测估计协方差，通过cholesky分解，获得初始量测噪声协方差的平方根
[0113]
在任意k时刻(k≥1)，选取2n+1个sigma点：
[0114][0115]
其中s
i,k-1
表示平方根阵的第i行。sigma点其对应的权值分别为
[0116][0117]
其中，w
i(c)
与w
i(m)
分别表示状态估计及其协方差的权值；pa描述了sigma点在均值周围的延伸程度，pb是先验分布描述参数。本算法中取pa＝1，pb＝2。
[0118]
s4.2：自适应平方根无迹卡尔曼平滑器的前向时间更新。其过程为
[0119][0120]
在时间更新过程中，通过qr分解和cholesky秩1更新获得平方根阵(18)式中，qr表示qr分解，用于求取协方差阵的平方根。cholupdate表示cholesky秩1更新，其中的sgn表示取符号。
[0121]
s4.3：前向平方根无迹卡尔曼滤波的量测噪声协方差阵平方根的自适应更新、量测更新及过程噪声协方差阵平方根的自适应更新。
[0122]
根据先验估计值再次选取2n+1个sigma点：
[0123][0124]
量测值的先验估计为
[0125]
[0126]
为进行量测噪声协方差阵的自适应更新，首先计算量测值yk与其先验估计的差在此基础上，量测噪声协方差阵按下式进行调整
[0127][0128]
其中，dk∈[0,1]是自适应滤波的加权系数，其数值决定了窗口长度。若为等加权方法，则当噪声参数不变时，噪声参数估计会越来越准确。但针对实际飞行数据记录噪声可能发生缓慢变化的情况，采用指数渐消记忆的加权方法
[0129]dk
＝(1-b)/(1-b
k+1
)
ꢀꢀ
(22)，
[0130]
其中b∈(0,1)为遗忘因子。从而，量测噪声协方差的平方根为：
[0131][0132]
在算法递推实现时，首先对其做qr分解，并进行两次cholesky秩1更新，量测噪声协方差阵平方根的自适应更新为
[0133][0134]
其中，diag是指取矩阵对角线元素组成对角阵。由此获得量测噪声协方差阵平方根的自适应更新。在此基础上，平方根无迹卡尔曼滤波的量测更新为
[0135][0136]
其中，仍采用qr分解和cholesky秩1更新获得量测估计协方差阵的平方根s
y,k
。p
xy,k
为先验状态估计和先验量测估计之间的互协方差阵。kk为卡尔曼滤波增益。最终取得状态估计协方差阵的平方根sk。
[0137]
按下式进行过程噪声协方差阵的自适应更新
[0138][0139]
对其作qr分解，并进行cholesky秩1更新
[0140][0141]
从而获得噪声协方差的平方根的自适应更新为：
[0142][0143]
s4.4：基于有限长飞行数据的自适应平方根无迹卡尔曼平滑估计。
[0144]
根据飞行数据离线事后分析的特点，利用固定时间间隔的全部量测值yn＝{y1,
y2,
…
,yn}，来获得该间隔内全部时刻的最优平滑估计。从而在基于步骤s4.1-s4.3获得每个时刻的状态估计后，再经最优平滑获得
[0145]
首先选取固定时间间隔k＝n-1,n-2,
…
,1,0,初始化滤波器
[0146][0147]
其中，上标s代表平滑估计。定义状态组合zk＝[xk，x
k+1
]
t
，在已知yk＝{y1,y2,
…
,yk}条件下，其后验估计为相应的协方差阵为
[0148][0149]
其中，为的协方差阵，为的协方差阵。为和的互协方差阵。从而，在已知测量值yk的情况下，的估计误差服从分布为x
k+1
的估计误差服从分布为根据贝叶斯公式，在已知测量值yk和的条件下，获得k时刻状态估计及其协方差为：
[0150][0151]
在通过平滑获得x
k+1
的基础上，可以得到p(xk|x
k+1
，yk)＝p(xk|x
k+1
，yk，y
k+1
，y
k+2
，
……
，yn)，即p(xk|x
k+1
，yk)＝p(xk|x
k+1
，yn)。这说明在平滑递推过程中，量测值yn,y
n-1
,
…
,y
k+1
对x
k+1
的平滑估计体现了条件概率推理中的马尔可夫性。从而在已知的基础上，根据贝叶斯公式获得组合状态zk＝[xk，x
k+1
]
t
的平滑估计及其协方差为：
[0152]
[0153]
上式中，为滤波器增益。从而得到k时刻的平滑估计及其协方差为：
[0154][0155]
在算法递推实现时，对平滑估计协方差阵的平方根进行qr分解和cholesky秩1更新，最终得到k时刻的平滑估计及其协方差阵的平方根为
[0156][0157]
综上所述，基于有限长飞行数据，通过(15)、(16)、(17)、(18)、(19)、(20)、(24)、(25)、(28)、(29)、(34)获得空速、迎角、侧滑角的最优平滑。
[0158]
以上的实施例仅是对本技术优选方式进行的描述，并非对本技术的范围进行限定，在不脱离本技术设计精神的前提下，本领域普通技术人员对本技术的技术方案做出的各种变形和改进，均应落入本技术权利要求书确定的保护范围内。