0 引言

表面热流辨识属于一类热传导反问题，它是通过测量导热材料内壁的温度测点的温度历程，反演出外壁受热面的热流时间历程。热传导反问题在航空航天、机械制造、车辆工程及生物工程等领域有着广泛应用^[1]。高超声速飞行器在大气层内飞行时，面临严重的气动加热问题，而对飞行器服役过程中的温度、热流等关键参数的测量是评价热防护材料使用性能、验证气动热模型和算法、指导热防护设计的必要手段^[2]。但对于驻点等热流密度大的区域，通常不能直接安置温度传感器或热流传感器进行测量，一方面是由于结构开孔导致结构强度下降以及缝隙加热，引起烧蚀不同步等结构匹配问题^[3]；另一方面，有些传感器本体材料不能承受过高热载荷，而且传感器的嵌入带来壁温的不连续以及周向的干扰，导致热流测量结果并非真实的气动热^[4]，所以需要采用表面热流辨识手段监测热流。

国内外对热传导反问题进行了大量的研究，通常做法是选取合适的目标函数，将辨识问题转化为优化问题求解。美国在20世纪60年代的Reentry-F飞行试验项目中，曾在飞行器弹身上布置了21个温度测点，利用类似顺序函数法的反演算法估算飞行器表面热流，从而预测边界层转捩位置^[5]，与其它方法测得的转捩位置基本一致。钱炜祺^[6]分别用顺序函数法和共轭梯度法研究了一维表面热流辨识方法，并实现了二维和三维非稳态表面热流辨识^[7-8]；薛齐文^[9]应用Tikhonov方法研究了一维热传导反问题中，内热源强度、导温系数及边界条件的多宗量辨识；DENG^[10]和智会强等^[11]研究了神经网络和遗传算法求解热传导反问题，但只考虑了简单的热流加载情况；CUI Miao^[12]采用无量纲化目标方程，对热流模型进行参数辨识，但局限于已知的热流函数形式；钱炜祺^[13]考虑到材料烧蚀后退，利用简化后的热解面烧蚀模型，对一维烧蚀表面热流辨识进行了研究；张聪^[14]利用简化的一维和二维传热模型进行了高超声速燃烧室壁面热流的辨识，在轴对称模型下取得了较好的效果。

在实际工程应用中，待辨识的热流通常是随时间和空间复杂变化的，已有的辨识方法不仅在三维情况下辨识精度不高，而且计算量庞大，难以满足实时性要求。针对某些工程问题通常最关心重点区域驻点热流随时间的变化，本文提出神经网络和顺序函数法结合的方法。在导热材料内壁布置适量温度传感器，对于每个传感器测得的温度数据采用一维的顺序函数法得到对应的表面热流辨识数据，再以每个测点得到的热流序列作为神经网络的输入，以外壁受热面对应时刻的驻点热流序列作为网络输出，并引入粒子群算法优化神经网络的初始权值和阈值，通过数值仿真或风洞试验的方法获得训练样本。该方法保留了顺序函数法良好的抗噪性，同时神经网络本身有强大的非线性映射能力，适合用于三维效应的修正。算例验证表明经过训练的神经网络预测的驻点热流和实际值十分吻合。

1 非稳态表面热流辨识方法

一维、三维表面热流辨识与二维情况类似，以二维为例，现介绍顺序函数法的辨识方法。对于典型的二维矩形域，考虑非线性非稳态热传导控制方程为:

$ \begin{array}{c}{\frac{\partial}{\partial x}\left(k_{x}(T) \frac{\partial T}{\partial x}\right)+\frac{\partial}{\partial y}\left(k_{y}(T) \frac{\partial T}{\partial y}\right)=\rho C_{P}(T) \frac{\partial T}{\partial t}} \\ {x \in[0, a], y \in[0, b]}\end{array} $

(1)

边界条件为：

$ \left\{\begin{array}{l}{x=0 :-k_{x}(T) \frac{\partial T}{\partial x}=q(t, y)} \\ {x=a :-k_{x}(T) \frac{\partial T}{\partial x}=0} \\ {y=0 :-k_{y}(T) \frac{\partial T}{\partial y}=0} \\ {y=b :-k_{y}(T) \frac{\partial T}{\partial y}=0}\end{array}\right. $

(2)

初始条件：

$ t=0 : T(x, y, 0)=T_{0} $

(3)

非线性偏微分方程可以采用拟线性隐式差分或有限体积法求解。

假设温度测点数量为M个，由于导热过程的时间延迟性，测点t时刻温度受到前r个时刻热流密度的影响，已知t_i前各时刻的热流值，所以要辨识t_i时刻的热流q(t_i, y)，即要求如下目标函数极小：

$ J = \sum\limits_{m = 1}^M {\sum\limits_{k = i}^{i + r - 1} {{{\left[ {\tilde T\left( {{x_m}, {y_m}, {t_k}} \right) - T\left( {{x_m}, {y_m}, {t_k}} \right)} \right]}^2}} } $

(4)

其中: $ \tilde{T}\left(x_{m}, y_{m}, t_{k}\right)$是测点m温度测量值，T(x_m, y_m, t_k)是对应位置温度计算值。由于计算温度场需要知道后续时刻热流值，所以假设q(t_i, y)到q(t_i+r-1, y)为线性变化，即：

$ \begin{array}{l}{q\left(t_{i+s}, y\right)=q\left(t_{i}, y\right)+s \times\left[q\left(t_{i}, y\right)-q\left(t_{i-1}, y\right)\right]}, \\ {s=1, 2, \cdots, r-1}\end{array} $

(5)

顺序函数法引入灵敏度函数的概念，即：

$ X(x, y, t)=\frac{\partial T(x, y, t)}{\partial q\left(t_{i}, y_{j}\right)} $

(6)

X(x, y, t)表示温度场对热流密度的灵敏度，灵敏度函数可以通过对热传导方程求偏导数得到：

$ \frac{\partial}{\partial x}\left(k_{x}(T) \frac{\partial X}{\partial x}\right) +\frac{\partial}{\partial y}\left(k_{y}(T) \frac{\partial X}{\partial y}\right)=\rho C_{P}(T) \frac{\partial X}{\partial t} \\ x \in[0, a], y \in[0, b], t \in\left[t_{i-1}, t_{i+r-1}\right] $

(7)

边界条件：

$ \left\{\begin{array}{l} x=0 :-k_{x}(T) \frac{\partial X}{\partial x}=n+1, t \in\left[t_{i+n-1}, t_{i+n}\right], \\~~~~~~~~~~~~~~~~~~~~~~ n=0, 1, \ldots, r-1 \\ x=a :-k_{x}(T) \frac{\partial X}{\partial x}=0 \\ y=0 :-k_{y}(T) \frac{\partial X}{\partial y}=0 \\ y=b :-k_{y}(T) \frac{\partial X}{\partial y}=0 \end{array}\right. $

(8)

初始条件：

$ X\left(x, y, t_{i-1}\right)=0, x \in[0, a], y \in[0, b] $

(9)

式(7)是一个近似方程，忽略了导热系数和比热对热流的导数项，由于二者相对较小，出于简化的目的可以忽略。

然后采用牛顿-拉夫逊算法得到待辨识热流分量的迭代修正公式：

$ \begin{array}{*{20}{l}} {{q^{k + 1}}\left( {{t_i}, {y_j}} \right) = {q^k}\left( {{t_i}, {y_j}} \right) + }\\ {\gamma \frac{{\sum\limits_{m = 1}^M {\sum\limits_{n = 0}^{r - 1} {\left[ {\tilde T\left( {{x_m}, {y_m}, {t_{i + n}}} \right) - T\left( {{x_m}, {y_m}, {t_{i + n}}} \right)} \right]} } X\left( {{x_m}, {y_m}, {t_{i + n}}} \right)}}{{\sum\limits_{m = 1}^M {\sum\limits_{n = 0}^N {{{\left[ {X\left( {{x_m}, {y_m}, {t_{i + n}}} \right)} \right]}^2}} } }}} \end{array} $

(10)

k为当前迭代次数，γ为防止迭代发散的收敛因子。每步迭代中，由修正后的热流可以计算得到新的目标函数值，当满足：

$ \frac{\left|J^{k+1}-J^{k}\right|}{J^{k}} \leqslant \varepsilon $

(11)

可以认为迭代已经收敛，ε为很小的正数，此时的q(t_i, y_j)作为辨识值。沿时间推进，即可得到连续的热流辨识序列。

顺序函数法中未来时间步r的选择直接影响到热流辨识的精度和鲁棒性，目前r的选取具有一定的经验性。为了选取最优的未来时间步，文献[15]提出了一种基于残差原理的最优未来时间步的确定方法，取不同的r，当满足方程(12)时，此时的r则为最优未来时间步。

$ {S_T}(r) = {\left[ {\frac{1}{{K - 1}}\sum\limits_{k = 1}^K {{{\left( {Y_j^k - T_j^k} \right)}^2}} } \right]^{1/2}} = \sigma $

(12)

其中:Y_j^k代表测点j在k时刻的测量温度，T_j^k则为求解正问题测点温度的计算值，K为时间离散点数，σ为温度测量误差的标准差。

为了避免r选择不当带来的影响，可以采用共轭梯度法(CGM)进行表面热流的辨识，它将辨识问题分解为热传导求解、灵敏度求解和伴随变量求解这三个问题进行求解，具体方法可参考文献[6]。

2 顺序函数法非稳态热流辨识算例

本文采用的算例计算域为：厚度L_x=0.01 m; 宽度L_y=0.1 m; 材料比热c_p=500 J/(kg·K); 密度ρ=8000 kg/m³; 导热系数k=80-0.02 T, W/(m·K); 初始温度300 K。测温点数量为21个，时间步长0.2 s，时间步r取25，实际加载的热流函数为：

$ \begin{aligned} q(t, y)=& 10000\left[1+4 \sin \left(\frac{t}{t_{\_\text {total }}} \pi\right)\right] \cdot \\ &\left[1+\sin \left(\frac{y}{L_{y}} \pi\right)\right] \end{aligned} $

(13)

总加热时间为t_{_total}=20 s，采用数值仿真的温度结果作为实测值，为检验算法的抗噪性能，在温度实测值基础上分别加入方差D=0.04和D=0.25的零均值均匀分布白噪声，真实热流，辨识热流以及加入噪声之后的辨识热流结果分别如图 1-4所示，Q表示热流密度。无噪声辨识相对误差和加入噪声后的辨识相对误差分布分别如图 5~图 7所示。

可见，无噪声时热流辨识结果最大相对误差约3.6%，加入D=0.04和0.25的噪声时，相对误差最大值约为5.7%和8.0%。由于噪声的增大，放大了反问题的不适定性，所以辨识结果相对误差也会增大，但总体精度仍然保持在较高水平，证明了顺序函数法用于非稳态表面热流辨识具有很好的准确性和抗噪性。

3 人工神经网络三维效应修正模型

根据国际著名神经网络专家Hecht Nielsen的观点，人工神经网络(Artificial Neural Networks, ANN)是由人工建立的、以有向图为拓扑结构的动态系统，它通过对连续或断续的输入作状态响应而进行信息处理^[16]。为了使神经网络修正模型尽可能简单稳定，本文选择目前应用广泛的BP(Back Propagation)神经网络。它是一种多层前馈神经网络，该网络拓扑结构如图 8所示。

BP网络主要特点是信号前向传递，误差反向传递。输入层和隐含层、隐含层和输出层神经元之前通过权值连接，隐含层和输出层的每个神经元分别有相应的阈值，隐含层和输出层还要设置相应的激励函数类型。一个典型的隐含层神经元输入输出关系为：

$ {H_j} = f\left( {\sum\limits_{i = 1}^n {{\omega _{ij}}} {x_i} - {a_j}} \right), j = 1, 2, \cdots , l $

(14)

其中，l为隐含层节点数，n为输入层节点数，f、ω、a分别为激励函数、权值和阈值。输出层输入输出关系类似。

神经网络初始化后，需要输入训练数据样本，若输出不满足误差要求，则误差转入反向传递，调整相应的权值和阈值，直到输出满足期望输出要求。经过大量的数据训练成熟的神经网络，将具有强大的非线性映射和泛化能力。

BP神经网络的特点，使其能够用于表面热流辨识三维效应的修正。基本思想是在飞行器表面峰值热流出现的重点区域内壁面，布置M个温度传感器，根据每个传感器实时测量的温度数据，用一维顺序函数法得到M个实时辨识的热流结果，数据归一化之后，将其作为神经网络的输入，神经网络的输出再反归一化，就是此刻该区域对应的峰值热流。

4 粒子群算法优化神经网络

遗传算法(GA)和粒子群算法(PSO)都是已被广泛应用的仿生全局优化方法。遗传算法模拟生物种群进化机制筛选最优个体，粒子群算法模拟鸟群觅食行为逼近食物位置。鉴于粒子群算法相比于遗传算法，不需要编码且没有交叉操作，使得模型更简单，大多数情况下，粒子群的单项信息共享机制通常能更快地收敛于最优解，故本文采用PSO算法优化神经网络的初始权值和阈值，得到最优的网络修正模型。

PSO算法首先在问题的可解空间中将一群粒子初始化，用位置、速度和适应度来表示该粒子特征，适应度值是根据适应度函数计算得到的，代表着粒子的优劣，每个粒子都代表问题的一个潜在解。粒子需要在解空间不停运动，其位置和速度的更新是根据个体最优值和群体最优值选择的，通过不断的更新群体，使群体最优值逼近全局最优解。

假设在一个D维的搜索空间中，n个粒子组成的种群为X=(X₁, X₂, …, X_n)，X_i=(X_i1, X_i2, …, X_iD)^T代表第i个粒子在搜索空间中的的位置，粒子速度表示为V_i=(V_i1, V_i2, …, V_iD)^T，个体极值表示为P_i=(P_i1, P_i2, …, P_iD)^T，全局极值表示为P_g=(P_g1, P_g2, …, P_gD)^T。则粒子速度和位置更新公式为：

$ \begin{aligned} V_{i d}^{k+1}=& \omega V_{i d}^{k}+c_{1} r_{1}\left(P_{i d}^{k}-X_{i d}^{k}\right)+\\ & c_{2} r_{2}\left(P_{g d}^{k}-X_{i d}^{k}\right) \end{aligned} $

(15)

$ X_{i d}^{k+1}=X_{i d}^{k}+V_{i d}^{k+1} $

(16)

其中：ω为惯性权重；c₁和c₂为加速度因子，为非负常数；r₁和r₂是介于0和1之间的随机数。通常要将粒子的位置和速度限制在一定区间，避免出现盲目搜索。为防止PSO早熟收敛，借鉴遗传算法的变异思想，在粒子每次更新之后，以10%的概率重新初始化粒子。

在实际应用中，将神经网络权值和阈值排列为行向量，对应粒子位置。这样每个粒子都对应一个神经网络，以训练样本的输出误差平方和作为粒子的适应度值，显然适应度值越低，代表该粒子越优。PSO算法最后得到的群体最优值作为神经网络模型初始权值和阈值。

5 三维效应修正模型算例验证

算例模型采用外径0.2m，厚度0.03m的空心半球壳，如图 9所示。材料比热c_p=500 J/(kg·K)，密度ρ=8000 kg/m³，导热系数k=80 W/(m·K)，初始温度300 K，时间步长0.2 s，总加热时间200 s，时间步r取30，测温点数量为9个，呈放射状均布于内壁面，布局在底面投影示意图如图 10所示。测点1位于球壳端部，内侧2-5测点距离测点1约0.03 m，外侧6-9测点距离测点1约0.06 m。

在球壳外表面施加随时间和空间变化的热流，热流函数形式为：

$ \begin{aligned} Q=& 1.5 \times 10^{6} \times\left[0.01+\frac{\mathrm{e}^{-50 R}}{\left(1+\mathrm{e}^{-50 R}\right)^{2}}\right]\cdot \\ &\left[\sin \left(\frac{t \pi}{200}\right)+0.1\right] \end{aligned} $

(17)

R为任意点与对称轴之间的距离，当t=100 s时，沿对称轴所在平面任意半圆周热流分布如图 11所示。

可见，该热流形式随空间变化明显，且与飞行器表面热流分布相似，可以较好地模拟气动加热。采用数值仿真软件计算出9个测点的时变温度数据，并叠加标准差0.5K的零均值白噪声作为后续输入数据，峰值热流为：

$ Q_{\max }=3.9 \times 10^{5} \times\left[\sin \left(\frac{t \pi}{200}\right)+0.1\right] $

(18)

在ANSYS APDL中加载(17)形式的热流边界条件，端部测点1的温度随时间变化如图 12所示。测点温度先升高后降低，是因为材料内部的热扩散是三维的，表面热流空间分布极不均匀，且随时间先增大后减小，测点周围单元的热扩散从以径向扩散升温为主，逐步转变为以周向扩散降温为主。

先仅采用端部测点1的温度数据进行一维辨识，结果如图 13所示，可见三维效应是很显著的，一维辨识误差很大，需要对其进行修正。

将9个测点的一维热流辨识结果作为神经网络输入，峰值热流序列作为网络输出。实际得到1000组输入输出数据，随机选取100个时刻测试数据，剩下900组数据用于网络训练。网络输入层为9个神经元，输出层为1个神经元，对于最佳隐含层节点数的选取，通常有如下经验公式：

$ l=n-1 $

(19)

$ l=\sqrt{m+n}+a $

(20)

$ l=\log _{2} n $

(21)

n，l，m分别为输入层，隐含层和输出层节点数，a为0~10之间的常数，本文分别设置隐含层神经元数目5~17。网络隐含层和输出层激发函数分别采用logsig和purelin函数，最大迭代训练次数为100，粒子群算法中，加速度因子设为1.49，惯性权重为0.5，粒子群规模为20个粒子，寻优迭代次数为100。当隐含层节点数为16时，粒子群群体最佳适应度值随寻优次数变化如图 14所示，图 15为训练后的网络峰值热流预测值与真实值比较结果，图 16为神经网络预测误差。将测试样本输出按时间顺序排列后与真实热流对比结果如图 17所示。

可见经过PSO算法优化的BP网络预测精度极高，隐含层节点16时，网络最大预测相对误差仅有0.20 %，表 1是不同隐含层节点数时，网络预测最大相对误差。

可以看出，隐含层节点数在较大范围变化时，神经网络预测精度仍然能保持在较高水平，隐含层节点数9时预测精度最优。

在工程应用中，通常希望在达到预期辨识精度前提下尽量减少测点数量。本文只选取1、4、6三个测点温度数据，设计神经网络结构为3输入，单输出，隐层节点数为7，相同数量样本训练后的预测结果和相对误差如图 18和图 19所示。

当测点数目减少到三个时，峰值热流预测最大相对误差为2.3%，虽然误差略有增大，但预测精度仍然较好，所以工程应用中可以适当减少测点数量以降低成本。

由于热流的时空分布形式不一，特别是随时间变化形式多样，需要验证基于一种热流时空分布建立的修正模型，在另外一种区别较大的热流加载下的表现。因此，本文根据半球壳在方程(17)的热流工况，建立并训练了BP网络修正模型(12个隐含层节点)，然后设计了具有不同时空分布形式的热流，其函数形式如式(22)：

$ \begin{aligned} Q &=10^{6} \times\left[0.02+\frac{\mathrm{e}^{-30 R}}{\left(1+\mathrm{e}^{-30 R}\right)^{2}}\right] \times \\ &\left(8 \times 10^{-5} t^{2}-1.6 \times 10^{-2} t+1\right) \end{aligned} $

(22)

其中:R为任意点与对称轴之间的距离，t为时间。

同样采用有限元软件计算出9个测点的时变温度数据，并叠加标准差0.5K的零均值白噪声作为后续输入数据，将一维的热流辨识结果输入前文已经建立的BP网络修正模型，驻点热流预测结果与实际值比较如图 20所示。

可见当应用于不同的热流工况时，修正模型的驻点热流预测值可以较好地跟随实际热流的变化趋势，但精度会显著降低，但相对于一维的辨识结果，仍然有了极大的提高。本文算例采用的材料较厚，测点离表面较远，同时模型头部半径较小，这些导致传热的三维效应非常明显，热流辨识的不适定性增加，这些也是图 20对比结果区别明显的原因。

6 结论

本文提出神经网络和一维顺序函数法结合的方法，构建表面热流辨识三维效应修正模型，由内壁测点的温度测量数据直接获得重点区域表面热流峰值实时辨识结果。通过数值仿真的算例测试结果可以看出，本文提出的方法对于峰值热流的辨识结果准确度高，模型训练可在线下进行，避免了三维辨识耗时大、精度低、无法实现实时在线辨识的缺点，同时具有良好的抗噪性和稳定性，在航天器表面峰值热流在线实时辨识中有广阔的应用前景。