基于迭代滤波的高超声速飞行器气动热预测

Aerothermal heat distribution prediction of hypersonic vehicles based on iterative filtering

摘要

气动热是气体快速流动或被高度压缩过程中产生的热量，过量气动热可能严重损害飞行器的结构和性能。优化飞行器形状、减少阻力和热负荷是解决气动热问题的主要策略，而气动热分布的快速预测则是其实现前提。工程计算方法和数值计算方法是两种常用预测方法，二者各有优缺点：作为一种经验预测方法，工程计算方法预测速度虽快但计算误差大；数值计算方法的预测精度虽高，但计算速度较慢。目前，能够同时兼顾速度和精度的气动热预测方法还不多见。本文提出了一种新的气动热分布预测方法，方法的基本思想是对工程计算方法进行改进，基于飞行器飞行试验数据对其经验参数进行精准辨识，在不降低计算速度的前提下提高其预测精度。鉴于常见的容积Kalman滤波（Cubature Kalman Filter，CKF）、扩展Kalman滤波（Extended Kalman filter，EKF）及无际Kalman滤波（Unscented Kalman Filter，UKF）等参数辨识方法无法同时满足精度和稳定性要求，本文提出了一种迭代滤波参数辨识方法，引入不动点迭代策略对容积Kalman滤波方法进行了改进。仿真结果显示，相比CKF、EKF及UKF方法，本文的滤波参数识别方法有更高的精度、更强的稳定性。此外，本文的预测方法不但速度快，而且预测误差小，比改进前的工程计算方法低50%，也比数值计算方法低20%。

Abstract

Aerothermal heating refers to the heat generated in the process of rapid flow or high compression of gases.Due to the serious damage of excessive heat on the structure and performance of vehicle， aerothermal heating problem poses a significant challenge to the design of hypersonic vehicles.To address the aerothermal heating problem， the optimization strategy of vehicle shape aiming to reduce drag and thermal loads heavily relies on the rapid prediction of aerothermal heat distribution.Engineering calculation method and numerical method are main methods used in aerothermal heat prediction.The former is fast but inaccurate， while the latter is accurate but slow.Nowadays， prediction methods that can balance the calculation speed and prediction accuracy are highly demanded but still be unknown.In this paper， to supply this gap， a new prediction method based on the iterative filting is proposed by improving the engineering calculation method through precise identification of its empirical parameters based on the flight test data.Considering that the main parameter identification methods such as cubature Kalman filter （CKF）， extended Kalman filter （EKF） and unscented Kalman filter （UKF） fail to balance the accuracy and stability， a iterative filtering identification method is proposed by introducing the fixed-point iteration strategy into the CKF method.Simulation results show that the proposed identification method outperforms the CKF，EKF and UKF in accuracy and stability.Meanwhile， the proposed prediction method significantly increases the prediction accuracy by reducing the prediction error by 50% of the original engineering calculation method and by 20% of the numerical method， while keeping the same calculation speed of engineering calculation method.

Graphical abstract

关键词

容积Kalman滤波 / 气动热 / 参数辨识 / 迭代滤波

Key words

cubature Kalman filter / aerothermal heat / parameter identification / iteration filtering

气动热是高超声速飞行器飞行过程中因空气摩擦或压缩效应而在其表面集聚的大量热能，可能导致飞行器外壁面温度急剧升高，破坏飞行器结构及飞行稳定性。随着高超声速飞行器技术的快速发展，气动热问题逐渐成为飞行器设计的关键^［1］，对飞行器的设计优化、材料选择、散热及冷却系统设计提出了严峻挑战。

解决气动热问题的前提是必须准确预测飞行器在不同飞行状态下的热流分布和热流值，此即气动热预测。目前，实验方法和计算方法是业界普遍使用的两类气动热分布预测方法。其中，实验方法主要包括挂飞试验和风洞试验，虽然能够提供较为精确的飞行数据，但受限于高昂的实验成本和漫长的实验周期，难以在短时间内获取大量的热流分布数据，因而更适合作为精细化预测手段。计算方法又可分为工程计算方法和数值计算方法。数值计算方法主要利用网格划分建立飞行器热流模型，预测精度虽高但计算量大、耗时长，且模型缺乏显式表达式，可能给后续数学处理带来挑战。工程计算方法基于实验得到的经验参数和流体力学原理建立物理方程，计算量虽然小，但精度较低，仅适合在初始预测阶段应用。

Kalman滤波（Kalman Filter， KF）是一类重要的参数辨识方法，目前已被广泛应用于空气动力学研究与应用^［2-5］。Gao等利用向后滤波算法减少了风扰情况下瞬时风场和空气数据的估计误差^［6］。Hong等引入了扩展Kalman滤波（Extended Kalman Filter，EKF）方法，以便合理地设计运载火箭的热保护系统^［7］。Tang等用无迹Kalman滤波（Unscented Kalman Filter，UKF）方法实现了气动力系统在线参数的辨识^［1］。值得注意的是，由于气动热问题过于复杂，研究者主要将滤波方法用于实验数据修正，如张峻宾利用滤波方法对风洞试验所得数据进行了调整或校正^［8］，通过对工程计算气动热模型进行参数辨识来实现气动热预测的研究还很少见。

相比EKF和UKF，容积Kalman滤波算法（Cubature Kalman Filter，CKF）^［9］具有更高的滤波精度、更强的稳定性。CFK基于容积变换原理进行状态估计，同时运用Spherical-Radial准则来选择一组具有相同权值的容积点，以估算函数的均值和协方差，减少非线性函数在线性化过程中的误差^［10］。进一步，为了获得更好的滤波性能，研究者将迭代思想引入了CKF算法，提出了迭代容积Kalman （Iterative Cubature Kalman Filter，ICKF） 算法^［11-13］。值得注意的是，这些迭代方法都基于Banach不动点迭代原理，对迭代函数的要求较为苛刻，往往难以直接应用于气动热预测问题。

本文旨在发展一种计算速度快同时预测精度高的气动热预测方法。本文的基本思路是：对工程计算气动热预测方法进行改进，基于飞行器的实测飞行试验数据对工程计算方法的经验参数进行辨识，通过提高参数值的精度来提高预测的精度，同时保持其计算量小、易于实现的优点。本文在CFK算法基础上引入了收敛条件更弱的不动点迭代策略，利用得到的迭代容积Kalman滤波算法（后文简称迭代滤波算法）对工程计算预测方法进行参数辨识。

1 预备知识

1.1 气动热模型

本文首先给出工程计算预测方法用到的气动热模型。

已知来流参数M（马赫数）、T（静温）、P（压力）、\rho（空气密度）、T_w（壁面温度）以及X（母线上位置坐标）。热流q的气动热模型为（未写出含义的量为中间量）：

q={\rho }^{\mathrm{*}}{u}_e(h_r-h_w)S_t/\sqrt[]{X},

S_t=f\times \frac{\mathrm{1}}{\mathrm{2}}{C}_f^{\mathrm{*}}{p}_r^{-\mathrm{0.667}},

C_f^{\mathrm{*}}=\frac{\mathrm{0.664}}{\sqrt[]{(R_e^{\mathrm{*}}/L)}},

R_e^{\mathrm{*}}/L={\rho }^{\mathrm{*}}{u}_e/{\mu }^{\mathrm{*}},

{\rho }^{\mathrm{*}}=P_e/(R{T}^{\mathrm{*}}),

{\mu }^{\mathrm{*}}=(\mathrm{1.789}\mathrm{4}\mathrm{e}-\mathrm{5})\times {\left(\frac{T^{\mathrm{*}}}{\mathrm{288.15}}\right)}^{\mathrm{1.5}}\times \frac{\mathrm{288.15}+\mathrm{110.4}}{T^{\mathrm{*}}+\mathrm{110.4}},

T^{\mathrm{*}}=\mathrm{0.5}(T_s+T_W)+\mathrm{0.22}{P}_r^{\mathrm{0.5}}(T_{\mathrm{0}}-T_e),

u_e=M_e{a}_e,

M_e=\sqrt[]{\left(\frac{M^{\mathrm{2}}+\frac{\mathrm{2}}{r-\mathrm{1}}}{\frac{\mathrm{2}r{M}^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta }{r-\mathrm{1}}-\mathrm{1}}+\frac{\frac{\mathrm{2}{M}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\beta }{r-\mathrm{1}}}{M^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta +\frac{\mathrm{2}}{r-\mathrm{1}}}\right)},

P_e=P\left(\mathrm{1}+\frac{\mathrm{2}r}{r+\mathrm{1}}\times (M^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta -\mathrm{1})\right),

{\rho }_e=\frac{\rho (r+\mathrm{1})M^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta }{(r-\mathrm{1})M^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta +\mathrm{2}},

T_e=\frac{T{P}_e\rho }{{\rho }_{e}P},

a_e=\sqrt[]{rR{T}_e},

h_r=\mathrm{1004}{T}_r,

T_r=T_e+\sqrt[]{P_r}(T_{\mathrm{0}}-T_s),

h_w=\mathrm{1004}{T}_w,

\mathrm{t}\mathrm{a}\mathrm{n}\delta =\frac{M^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta -\mathrm{1}}{\left[M^{\mathrm{2}}\left(\frac{r+\mathrm{1}}{\mathrm{2}}-\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta \right)+\mathrm{1}\right]\mathrm{t}\mathrm{a}\mathrm{n}\beta },

T_{\mathrm{0}}=T\left(\mathrm{1}+\frac{r-\mathrm{1}}{\mathrm{2}}{M}^{\mathrm{2}}\right),

其中， R、r及壁面迎风角度\delta为已知参数，

P_r=\left\{\begin{array}{ll}\mathrm{0.7},& \mathrm{层流},\\ \mathrm{0.89},& \mathrm{湍流},\end{array}\right.

f=\left\{\begin{array}{ll}\mathrm{1.0},& \mathrm{f}\mathrm{l}\mathrm{a}\mathrm{t},\\ \sqrt[]{\mathrm{3}},& \mathrm{b}\mathrm{l}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{。}\end{array}\right.

以上气动热模型的优点在于计算速度快，缺点在于计算精度较差。作为经验参数，模型中的参数值是不准确的。本文考虑根据高超声速气动热实验获得的实测数据对模型中的部分经验参数值进行精准辨识，以便在不降低计算速度的前提下获得更加准确的气动热预测效果。

1.2 容积Kalman滤波

考虑如下的通用离散非线性状态与量测系统：

\left\{\begin{array}{l}{x}_{k+\mathrm{1}}=l(x_k)+u_k,\\ z_k=h(x_k)+v_k\end{array}\right.（1）

其中，h(x)和l(x)分别为系统非线性量测方程和非线性状态转移方程的传递函数，x_k、u_{k-\mathrm{1}}、z_k和v_k分别为系统的状态量、过程噪声、量测量和量测噪声，u_k和v_k相互独立且满足

u_k~\mathrm{𝒩}(\mathrm{0},{\mathit{Q}}_{\mathit{k}})，v_k~\mathrm{𝒩}(\mathrm{0},{\mathit{R}}_{\mathit{k}})，

{\mathit{Q}}_{\mathit{k}}和{\mathit{R}}_{\mathit{k}}分别表示u_k和v_k的协方差阵。

根据容积Kalman滤波算法中容积点的选取原则，设状态空间向量的维数为n，容积点采样点数为m=\mathrm{2}n。又设{\xi }_i为采样容积点，{\omega }_i为{\xi }_i所占权值，即

{\xi }_i=\sqrt[]{\frac{m}{\mathrm{2}}}\left\{\left[\begin{array}{c}\mathrm{1}\\ \mathrm{0}\\ ⋮\\ \mathrm{0}\end{array}\right],\cdots ,\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ ⋮\\ \mathrm{1}\end{array}\right],\left[\begin{array}{c}-\mathrm{1}\\ \mathrm{0}\\ ⋮\\ \mathrm{0}\end{array}\right],\cdots ,\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ ⋮\\ -\mathrm{1}\end{array}\right]\right\}，

{\omega }_i=\mathrm{1}/m（2）

假设k-\mathrm{1}时刻系统概率分布的均值为{\overline{x}}_{k-\mathrm{1}}，协方差{\mathit{P}}_{k-\mathrm{1}}正定。设{\mathit{S}}_{k-\mathrm{1}}=\sqrt[]{{\mathit{P}}_{k-\mathrm{1}}}为{\mathit{P}}_{k-\mathrm{1}}经过Cholesky分解后的结果。

算法1 容积Kalman滤波算法

已知量：{\widehat{x}}_{\mathrm{0}},{\mathit{P}}_{\mathrm{0}},{\mathit{Q}}_{\mathrm{0}},{\mathit{R}}_{\mathrm{0}}。

For kK do

预测步：计算容积点，

{\mathit{S}}_{k-\mathrm{1}}=\sqrt[]{{\mathit{P}}_{k-\mathrm{1}}}，

{\mathrm{𝒳}}_{k-\mathrm{1}}^i={\widehat{x}}_{k-\mathrm{1}}+{\mathit{S}}_{k-\mathrm{1}}{\xi }_i，

容积点状态预测，

{\widehat{\mathrm{𝒳}}}_{k|k-\mathrm{1}}^i=f\left({\mathrm{𝒳}}_{k-\mathrm{1}}^i\right),i=\mathrm{0},\cdots ,m，

目标状态预测，

{\widehat{x}}_{k|k-\mathrm{1}}=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}{\widehat{\mathrm{𝒳}}}_{k|k-\mathrm{1}}^i,

预测协方差，

{\mathit{P}}_{k\left|k-\mathrm{1}\right.}=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}\left({\widehat{\mathrm{𝒳}}}_{k|k-\mathrm{1}}^i-{\widehat{x}}_{k|k-\mathrm{1}}\right){\left({\widehat{\mathrm{𝒳}}}_{k|k-\mathrm{1}}^i-{\widehat{x}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}}+{\mathit{Q}}_{k-\mathrm{1}};

更新步：计算容积点，

{\mathit{S}}_k=\sqrt[]{{\mathit{P}}_{k|k-\mathrm{1}}},

{\mathrm{𝒳}}_{k|k-\mathrm{1}}^i={\widehat{x}}_{k|k-\mathrm{1}}+{\mathit{S}}_{k|k-\mathrm{1}}{\xi }_i,

容积点量测预测，

{\widehat{Z}}_{k|k-\mathrm{1}}^i=h\left({\mathrm{𝒳}}_{k-\mathrm{1}}^i\right),i=\mathrm{0},\cdots ,m,

目标量测预测，

{\widehat{z}}_{k|k-\mathrm{1}}=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}{\widehat{Z}}_{k|k-\mathrm{1}}^i,

量测协方差，

{\mathit{S}}_{k\left|k-\mathrm{1}\right.}=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}\left({\widehat{Z}}_{k|k-\mathrm{1}}^i-{\widehat{z}}_{k|k-\mathrm{1}}\right){\left({\widehat{Z}}_{k|k-\mathrm{1}}^i-{\widehat{z}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}}+{\mathit{R}}_k,

状态量测协方差，

{\mathit{C}}_k=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}\left({\widehat{\mathrm{𝒳}}}_{k|k-\mathrm{1}}^i-{\widehat{x}}_{k|k-\mathrm{1}}\right){\left({\widehat{Z}}_{k|k-\mathrm{1}}^i-{\widehat{z}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}},

Kalman增益，

{\mathit{K}}_k={\mathit{C}}_k/{\mathit{S}}_{k\left|k-\mathrm{1}\right.};

估计：状态估计，

{\widehat{x}}_k={\widehat{x}}_{k|k-\mathrm{1}}+{\mathit{K}}_k\left(z_k-{\widehat{z}}_{k|k-\mathrm{1}}\right),

协方差估计，

{\mathit{P}}_k={\mathit{P}}_{k\left|k-\mathrm{1}\right.}-{\mathit{K}}_k{\mathit{S}}_{k\left|k-\mathrm{1}\right.}{\mathit{K}}_k^{\mathrm{T}};

end for。

1.3 迭代收敛性

不动点定理是设计迭代算法的基本工具。本节将介绍用到的不动点定理及推论。根据Banach不动点定理，在Banach空间内，压缩映射都有且仅有一个不动点^［14］。值得注意的是，已有的迭代容积Kalman滤波方法的收敛条件大都基于Banach不动点定理，并且对映射的连续可微性有较高要求^［12-13］。

对于容积Kalman滤波算子g(x)，若g(x)是压缩映射，则可证明迭代容积Kalman滤波方法存在不动点。根据容积Kalman滤波的定义，若状态方程为恒等方程

l(x+{\mathit{S}}_{k|k-\mathrm{1}}{\xi }_j)=x+{\mathit{S}}_{k|k-\mathrm{1}}{\xi }_j，

则有

\begin{array}{l}g(x)=x+\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}{\mathit{S}}_{k-\mathrm{1}}{\xi }_i+\\ \begin{array}{cc}& \end{array}{\mathit{K}}_k\left(z_k-\frac{\mathrm{1}}{m}{\displaystyle \sum _{j=\mathrm{1}}^m}h(x+{\mathit{S}}_{k|k-\mathrm{1}}{\xi }_j)\right)\end{array}（3）

从而，g\left(x\right)为压缩映射的条件是||\mathrm{I}-{\mathit{K}}_{k}h||\mathrm{1}，此时的容积Kalman滤波算法收敛。

在气动热模型中，量测函数并不一定连续，因而\mathrm{I}-{\mathit{K}}_{k}h也不一定是压缩映射。本文考虑逐步放宽对映射的要求，将原压缩映射放宽为非扩张映射甚至不连续映射。

首先，本文给出SKC条件作为进一步提出非扩张映射的基础。

定义1.1 任意x,y\in \mathbf{C},如果

\frac{\mathrm{1}}{\mathrm{2}}|Tx-x|\le |x-y|

蕴含

|Tx-Ty|\le N(x,y)，

则称映射T:\mathbf{C}\to \mathbf{C}满足SKC条件，其中

\begin{array}{l}N(x,y)=\mathrm{m}\mathrm{a}\mathrm{x}\left\{|x-\right.y|,\\ \begin{array}{cc}& \frac{|x-Tx|+|Ty-y|}{\mathrm{2}},\frac{|Ty-x|+|x-Ty|}{\mathrm{2}}\}\end{array}\end{array}（4）

压缩映射显然满足SKC条件。因此，SKC条件是对压缩映射的延拓。进而，根据上述定义，Karapinar和Tas提出了以下Suzuki 映射的推广。

定理1.2^［15］ 令T是 Banach 空间中紧凸子集K上的映射，且满足SKC条件。对于由x_{\mathrm{1}}\in K以及

x_{n+\mathrm{1}}=\lambda T{x}_n+(\mathrm{1}-\lambda )x_n

定义的序列\{x_n\}，假设

\underset{n\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}|T{x}_n-x_n|=\mathrm{0}成立，

则\{x_n\}强收敛到T的一个不动点，其中\lambda \in \left.\left[\mathrm{1}/\mathrm{2,1}\right.\right)\mathrm{。}

根据定理1.2，有如下关于紧凸集迭代序列的收敛结果。

推论1.3^［15］ 在紧凸集上，设x_{\mathrm{0}}生成的迭代序列为\left\{x_n\right\}\mathrm{。}令x_n=T^n{x}_{\mathrm{0}}\mathrm{。}如果x_n满足条件

\underset{n\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}|x_n-x_{n-\mathrm{1}}|=\mathrm{0}，

则\left\{x_n\right\}强收敛。

因

\left|x_n-x_{n-\mathrm{1}}\left|=\right|T{x}_{n-\mathrm{1}}-x_{n-\mathrm{1}}\right|\to \mathrm{0},n\to \mathrm{\infty },

若

|g^n(x)-g^{n-\mathrm{1}}(x)|\to \mathrm{0},n\to \mathrm{\infty }，

则\underset{n\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{g}^n(x)存在，且该极限就是g(x)的不动点。因此，推论1.3中的条件比Banach不动点定理中的压缩映射条件放宽许多。接下来本文将利用以上结果导出迭代容积Kalman滤波算法。

2 迭代容积Kalman滤波算法

对于系统（1），设j为迭代次数并设最大迭代次数为N\mathrm{。}本文在普通容积Kalman滤波算法中加入迭代策略，并另设k为当前时刻，时序数据长度指定最大时刻为K。

算法2 基于迭代容积Kalman滤波算法

已知量：{\mathit{Q}}_{\mathrm{0}},{\mathit{R}}_{\mathrm{0}},

{\widehat{\mathit{P}}}_{\mathrm{0}}^{\mathrm{0}}=E\left[\left(x_{\mathrm{0}}-{\widehat{x}}_{\mathrm{0}}^{\mathrm{0}}\right){\left(x_{\mathrm{0}}-{\widehat{x}}_{\mathrm{0}}^{\mathrm{0}}\right)}^{\mathrm{T}}\right],

\begin{array}{c}\begin{array}{l}{\widehat{x}}_{\mathrm{0}}^{\mathrm{0}}=E[x_{\mathrm{0}}^{\mathrm{0}}]={\int }_{{\mathrm{R}}^s}f(x_{\mathrm{0}}^{\mathrm{0}})p(x_{\mathrm{0}}^{\mathrm{0}})\mathrm{d}{x}_{\mathrm{0}}^{\mathrm{0}}=\\ {\int }_{{\mathrm{R}}^s}f(x_{\mathrm{0}}^{\mathrm{0}})\times N(x_{\mathrm{0}}^{\mathrm{0}};{\mathit{P}}_{\mathrm{0}})\mathrm{d}{x}_{\mathrm{0}}^{\mathrm{0}}\mathrm{。}\end{array}\end{array}

\mathrm{F}\mathrm{o}\mathrm{r}\begin{array}{cc}& jN\end{array}\begin{array}{cc}& \mathrm{d}\mathrm{o}\end{array}

\mathrm{F}\mathrm{o}\mathrm{r}\begin{array}{cc}& kK\end{array}\begin{array}{cc}& \mathrm{d}\mathrm{o}\end{array}

预测步：计算容积点，

{\mathit{S}}_{k-\mathrm{1}}=\sqrt[]{{\widehat{\mathit{P}}}_{k-\mathrm{1}}^j},

X_{k-\mathrm{1},i}^j={\widehat{x}}_{k-\mathrm{1}}^j+{\mathit{S}}_{k-\mathrm{1}}{\xi }_i,

容积点状态预测，

{\mathrm{𝒳}}_{k-\mathrm{1},i}^j=l({\mathrm{𝒳}}_{k-\mathrm{1},i}^j),i=\mathrm{0},\cdots ,m,

目标状态预测，

{\overline{x}}_k^j=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}{\mathrm{𝒳}}_{k,i}^j,

预测协方差，

{\mathit{P}}_k^j=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}\left({\mathrm{𝒳}}_{k,i}^j-{\overline{x}}_k^j\right)({\mathrm{𝒳}}_{k,i}^j-{\overline{x}}_k^j{)}^{\mathrm{T}}+{\mathit{Q}}_{k-\mathrm{1}};

更新步： 计算容积点，

{\mathit{S}}_k^j=\sqrt[]{{\mathit{P}}_k^j},

{\mathrm{𝒳}}_{k,i}^j={\overline{x}}_k^j+{\mathit{S}}_k^j{\xi }_i，

容积点量测预测，

{\mathit{Z}}_{k,i}^j=h\left({\mathit{X}}_{k,i}^j\right),i=\mathrm{0},\cdots ,m,

目标量测预测，

{\overline{z}}_k^j=\frac{\mathrm{1}}{m}{{\displaystyle \sum _{i=\mathrm{1}}^m}Z}_{k,i}^j，

量测协方差，

{\mathit{S}}_k^j=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}\left(Z_{k,i}^j-{\overline{z}}_k^j\right)(Z_{k,i}^j-{\overline{z}}_k^j{)}^{\mathrm{T}}+{\mathit{R}}_k,

状态量测协方差，

{\mathit{C}}_k^j=\frac{\mathrm{1}}{m}{\displaystyle \sum _{i=\mathrm{1}}^m}\left({\mathrm{𝒳}}_{k,i}^j-{\overline{x}}_k^j\right)(Z_{k,i}^j-{\overline{z}}_k^j{)}^{\mathrm{T}},

Kalman增益，

{\mathit{K}}_k^j={\mathit{C}}_k^j/{\mathit{S}}_k^j;

估计：根据定理1.2及推论1.3，设置t_j为每一轮的迭代因子，通过Matlab内置的随机种子函数给出，建立由基于过往信息的预测值和当前量测信息的状态估计及其协方差估计。

状态估计：

t_j\in [\mathrm{1}/\mathrm{2,1}]，

{\widehat{x}}_k^j=(\mathrm{1}-t_j){\overline{x}}_k^j+t_j\left[{\overline{x}}_k^j+{\mathit{K}}_k^j(z_k-{\overline{z}}_k^j)\right];

协方差估计：

{\widehat{\mathit{P}}}_k^j={\mathit{P}}_k^j-{\mathit{K}}_k^j{\mathit{S}}_k^j({\mathit{K}}_k^j{)}^{\mathrm{T}},

end for

end for。

接下来，本文利用以上的迭代滤波算法对工程计算预测方法中的气动热模型经验参数进行修正。

在飞行试验中，气流和飞行器参数P_r和f以及空气流动参数trans是不确定的，也是工程计算方法无法取得良好效果的重要原因。本文考虑将这些参数作为未知量，基于飞行实测数据进行精准辨识。

计算机验证发现，q=h\left(\begin{array}{c}{P}_r\\ f\\ \mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\end{array}\right)函数在某些点附近满足推论1.3中的不动点收敛条件，从而可以利用迭代滤波算法进行参数精准辨识。将P_r、f和trans视作状态量，则状态方程为

x_k=\left(\begin{array}{c}{P}_r^k\\ f^k\\ \mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}{\mathrm{s}}^k\end{array}\right)=l(x_{k-\mathrm{1}})\equiv I^{\mathrm{*}}\left(\begin{array}{c}\mathrm{0.89}\\ \sqrt[]{\mathrm{3}}\\ \mathrm{340}\end{array}\right)（5）

量测方程为气动热模型中的热流方程，则

z_k:=q_k=h(x_k)（6）

通过迭代滤波算法得到P_r、f和trans随时间的变化后，将其存储到离线文件中作为气动热模型参数辨识结果。鉴于该结果是一列时序数据，为了让模型具有泛化性，且适用于更多种传感器和流场类型（层流或湍流），本文分别对这些参数取均值，得到平滑后的参数，最后代入模型得到修正的气动热预测模型。

3 仿真分析

用Matlab 构建仿真平台，并利用一组气动热试验数据建立修正的气动热预测模型，然后用另一组试验数据验证改进的预测模型的效果。

本文分别进行参数辨识效果比较和气动热预测效果比较。在参数辨识效果比较中，本文分别用UKF算法、PF算法、CKF算法及本文提出的迭代滤波算法（本节中将其简记为ICKF）基于第一组实测数据进行参数辨识、调整模型参数，然后比较修正模型的预测结果与第一组实测数据的差异。在气动热预测效果比较部分，本文使用基于参数平滑后的ICKF预测算法与工程计算及数值计算预测方法对第二组数据进行预测，并比较其误差和运算时间。

3.1 参数辨识效果

在相同的初始状态和协方差下，利用第一组数据分别建立基于UKF、PF、CKF及ICKF参数辨识算法的气动热预测模型。选取不同流场下的550组数据（包括来流参数和对应飞行器上特定测点的热流值）作为各参数辨识算法的训练集，并取迭代次数为5。分别使用UKF算法、PF算法、CKF算法及ICKF算法（t_j取1）对热流模型进行参数辨识，并将辨识后的参数值按照时间序列代入模型，随后再将各气动热预测模型的预测结果与真实值进行比较，结果如图1~图4所示，图中的蓝色线表示参数辨识前气动热模型的预测结果与真实值的相对误差，红色线表示参数辨识后气动热模型的预测结果与真实值的相对误差。

从图1~图4可以发现，经CKF算法修正后的气动热预测模型的预测结果与真实值差距很大，经PF算法及UKF算法修正后的气动热预测模型的预测效果波动较大，而使用ICKF算法修正后的模型的预测结果与真实值的差距较小，且较为稳定。

3.2 气动热预测效果

首先使用UKF、PF、CKF及ICKF算法分别进行参数辨识，然后对所得参数取均值作为修正模型的新参数，代入原预测模型得到基于UKF、PF、CKF及ICKF算法的修正预测模型。对第一组实测数据进行预测并与真实值进行比较，结果如图5所示。可以发现，基于ICKF算法的预测模型确实表现最好。

进一步，本文比较基于ICKF算法的预测模型与原工程计算方法及数值计算方法对于第二组实测数据的预测效果，结果如图6所示。可以看到，基于ICKF算法的预测模型的误差比工程计算方法低50%，比数值计算方法也低20%。

此外，表1给出了基于ICKF算法的预测模型、工程计算预测方法及数值计算预测方法的计算时间。可以看到，前者的运算速度相对工程计算预测方法并没有明显降低。

4 结论

本文提出了基于迭代容积Kalman滤波的气动热预测方法，通过对气动热工程计算预测方法中的部分经验参数进行辨识和修正来提高其预测精度，同时保持其计算速度。本文的方法有望为高空长航时高超声速飞行器设计提供重要帮助。

参考文献

[1] Tang X， Jiang Z， Chen H.Online identification of aerodynamic parameters of experimental rockets based on unscented Kalman filtering ［J］.Int J Aerospece Eng， 2024， 2024： 4541120.

[2] Gontumukkala S S T， Godavarthi Y S V， Godavarthi B R R T， et al.Kalman filter and proportional navigation based missile guidance system ［C］//2022 8th ICACCS.Piscataway： IEEE， 2022： 1731-1736.

[3] Xu Z， Mao B Q， Xu L C， et al.Study of missile radiator predicting and tracking technology on the basis of Kalman filter and mean-shift algorithm ［J］.Adv Mater Res， 2011， 383-390： 1584-1589.

[4] Siouris G M， Chen G， Wang J.Tracking an incoming ballistic missile using an extended interval Kalman filter ［J］.IEEE T Aero Elec Sys， 1997， 33（1）： 232-240.

[5] Stoffel T D， Karlgaard C D， White T R， et al.Fusion of in-flight aerothermodynamic heating sensor measurements using Kalman filtering ［J］.J Spacecraft Rockets， 2024， 61（2）： 599-610.

[6] Gao Z， Wang H， Xiang Z， et al.Flight data-based wind disturbance and air data estimation ［J］.Atmosphere， 2021， 12（4）： 470.

[7] Hong Y， Ma Y， Wen S， et al.A reconstructed approach for online prediction of transient heat flux and interior temperature distribution in thermal protect system ［J］.Int Commun Heat Mass， 2023， 148：107055.

[8] Zhang J B， Xu X B， Wang X， et al.Data processing technology of balanced dynamic characteristics based on wavelet reconstruction ［J］.J BUAA， 2023， 49（6）： 1362-1371.

[9] 张峻宾， 许晓斌， 王雄， 等.基于小波重构的天平动态特性数据处理技术 ［J］.北京航空航天大学学报， 2023， 49 （6）： 1362-1371.

[10] Arasaratnam I， Haykin S.Cubature Kalman filters ［J］.IEEE T Autom Control， 2009， 54（6）： 1254-1269.

[11] Hu Z T， Yuan G Y， Hu Y M， et al.Training method of neural network based on cubature Kalman filter ［J］.Control and Decision， 2016， 31（2）： 355-360.

[12] 胡振涛， 袁光耀， 胡玉梅， 等.基于容积卡尔曼滤波的神经网络训练算法［J］.控制与决策， 2016， 31（2）： 355-360.

[13] Wang J， Chen X， Shi C， et al.Robust M-estimation-based ICKF for GNSS outlier mitigation in GNSS/SINS navigation applications ［J］.IEEE T Instrum Meas， 2023， 72： 8505617.

[14] Cui B， Chen X， Xu Y， et al.Performance analysis of improved iterated cubature Kalman filter and its application to GNSS/INS ［J］.ISA T， 2017， 66： 460-468.

[15] Mu J， Cai Y L.Iterated cubature Kalman filter and its application ［C］//2011 IEEE CYBER.Piscataway： IEEE， 2011： 33-37.

[16] Shi J.Iterative algorithms for fixed points of several generalized nonexpansive mappings in Banach spaces ［D］.Wuhan： Wuhan University， 2020.

[17] 史杰.Banach 空间中几类广义非扩张映射不动点的迭代算法［D］.武汉：武汉大学， 2020.

[18] Karapinar E， Tas K.Generalized （C）-conditions and related fixed point theorems ［J］.Comput Math Appl， 2011， 61（11）： 3370-3380.

基金资助

四川省自然科学基金(2023NSFSC0061)