CN106289568A - 一种爆震火焰温度测量***及三维温度场的重建方法 - Google Patents

Abstract

本发明涉及一种爆震火焰温度测量***及三维温度场的重建方法，利用一定数量声学传感器和扬声器对相同位置处的爆震燃烧室温度场进行测量。对数据进行分析处理，然后利用互相关技术在信噪比低的情况下仍然能够得到准确的时延估计。针对燃烧室环境混响干扰，研究基于倒谱处理的相关算法，提高时延估计精度。再在MATLAB中利用基于LSQR算法(Least Squares QR‑Decomposition)的区域离散法重建脉冲爆震发动机燃烧室内的温度场。

Description

一种爆震火焰温度测量***及三维温度场的重建方法

技术领域

本文涉及到一种燃烧温度场分布的测量方法，尤其是爆震燃烧室温度场的测量方法，属于燃烧诊断领域。

背景技术

文献“基于声波传感器航空发动机燃烧室出口温度分布测量研究，沈阳航空工业学院报，Vol.26，No.2”公开了一种基于声波传感器测量航空发动机燃烧室出口温度分布的方法。该方法先利用多个声学传感器测量声波在该温度场飞渡所用的时间，再由声波的传播速度与介质温度的关系，得出该传播路径的平均温度，然后利用基于最小二乘法的计算机断层技术(Computed Tomography)实现二维温度场的仿真重建。文献所述方法，该方法较精确的测出了燃烧室出口温度分布，较好的反映了火焰的真实温度。但是由于燃烧室内环境恶劣，出口噪声很大，而且声波会有反射，所以不能准确测量燃烧室出口温度分布。而且，该法只测量了燃烧室出口温度，不能研究燃烧室内温度场分别的特点。因此，文献所述方法很难对恶劣环境中的燃烧室整体温度分布特点进行准确测量。

发明内容

要解决的技术问题

脉冲爆震发动机(PDE，Pulse Detonation Engine)是一种基于爆震燃烧的新概念发动机。爆震燃烧产生的爆震波使可爆燃料的压力、温度迅速升高。为了克服恶劣环境中火焰温度难以测量的难题，本发明提出了一种爆震火焰温度测量***及三维温度场的重建方法。

技术方案

本发明利用一定数量声学传感器和扬声器对相同位置处的爆震燃烧室温度场进行测量。对数据进行分析处理，然后利用互相关技术在信噪比低的情况下仍然能够得到准确的时延估计。针对燃烧室环境混响干扰，研究基于倒谱处理的相关算法，提高时延估计精度。再在MATLAB中利用基于LSQR算法(Least Squares QR-Decomposition)的区域离散法重建脉冲爆震发动机燃烧室内的温度场。

一种爆震火焰温度测量***，其特征在于包括带声卡和显示器的主机、功率放大器、继电器、数采卡、端子板、信号调理器、扬声器和麦克风；扬声器和麦克风位于在爆震燃烧室内，带声卡和显示器的主机、功率放大器、继电器、数采卡、端子板和信号调理器位于爆震燃烧室外，扬声器与继电器连接，麦克风与信号调理器连接；主机中的声波测温软件***发出数字声波信号，经声卡转换为模拟信号后，再经功率放大器放大输出到继电器开关的一端；与此同时，主机上的声波测温软件***向数采卡上的数字I/O输出高电平，高电平由与数采卡相连的端子板引出，经信号控制电缆去启动继电器，将功率放大器输出的模拟信号传输到扬声器，扬声器将模拟信号转化为声波信号，麦克风将接收到的声波信号转换为电压信号，经信号调理器放大后，被数采卡采集，然后输入到声波测温软件***进行三维温度场的重建。

扬声器和麦克风的数量分别为20个，一个扬声器和一个麦克风组成一组，20组均匀的分布在立方体空间的各个边上。

一种对爆震火焰温度测量***测得的电压信号进行三维温度场的重建方法，其特征在于步骤如下：

步骤1：采用倒谱法对两麦克接收到的电压信号进行倒谱变换：

其中，接收到的电压信号为：

x₁(t)＝h₁(t)*s(t-τ₁)+n₁(t)

x₂(t)＝h₂(t)*s(t-τ₂)+n₂(t)

其中，h₁(t)、h₂(t)为燃烧室的单位冲激响应；n₁(t)、n₂(t)为麦克风接收的噪声信号；s(t-τ₁)、s(t-τ₂)为两麦克风接收扬声器发出的信号；

进行倒谱变换后的电压信号为：

x₁(t)＝s(t-τ₁)+n₁(t)

x₂(t)＝s(t-τ₂)+n₂(t)

步骤2：对倒谱变换后的x₁(t)和x₂(t)进行互相关运算：

$R_{x_1{x}_2}\left(\mathrm{\τ}\right)=\frac{1}{2\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}{x}_1(t-\mathrm{\τ})x_2\left(t\right)dt=\frac{1}{2\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\[s_1(t-\mathrm{\τ})+n_1(t-\mathrm{\τ})\]\[{\mathrm{As}}_2\left(t\right)+n_2\left(t\right)\]dt$

化简可得：

$R_{x_1{x}_2}\left(\mathrm{\τ}\right)=R_{s_1{s}_2}\left(\mathrm{\τ}\right)+R_{s_1{n}_2}\left(\mathrm{\τ}\right)+R_{n_1{s}_2}\left(\mathrm{\τ}\right)+R_{n_1{n}_2}\left(\mathrm{\τ}\right)$

其中，A为声波传播过程中的衰减系数，为信号的互相关运算结果，分别为信号和噪声的互相关运算结果，是噪声互相关运算结果，易知为零，可简化为：

$R_{x_1{x}_2}\left(\mathrm{\τ}\right)=R_{s_1{s}_2}\left(\mathrm{\τ}\right)=\frac{1}{2\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}{s}_1(t-\mathrm{\τ})s_2\left(t\right)dt=\frac{1}{2\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}s(t+\mathrm{\Δ}t-\mathrm{\τ})s\left(t\right)dt=R_s(\mathrm{\τ}-\mathrm{\Δ}t)$

其中，Δt为时间延迟，因为R_s(τ-Δt)≤R_s(0)，当τ＝Δt时，取得最大值，则通过计算，互相关函数峰值出现的时刻即为时间延迟Δt即飞渡时间τ；

步骤3：采用LSQR算法重建温度场：

将三维燃烧室测量空间均分为3×3×3＝27个空间网格，则

${\mathrm{\τ}}_i=\underset{j=1}{\overset{N}{\Σ}}{w}_{ij}{f}_j,i=1,2...,58$

式中：w_ij为第i条路径通过第j个空间网格的长度；τ_i为测得的第i条声波路径的声波飞渡时间；f_j为待测区域内声波传播速度的倒数；

控制20个扬声器依次发射声波信号，组成一个测量周期，则得到58条声波路径，相应的线性方程组为：

$\begin{array}{c}{w}_{11}{f}_1+w_{12}{f}_2+\mathrm{...}+w_{1,27}{f}_{27}={\mathrm{\τ}}_1\\ w_{21}{f}_1+w_{22}{f}_2+\mathrm{...}+w_{2,27}{f}_{27}={\mathrm{\τ}}_2\\ .\\ .\\ .\\ w_{58,1}{f}_1+w_{58,2}{f}_2+\mathrm{...}{w}_{58,27}{f}_{27}={\mathrm{\τ}}_{58}\end{array}$

该方程组用矩阵表示为

Ax＝b，即为三维温度场重建的数学模型

式中：为58×27维权因子矩阵，反映了第j空间网格对第i条声波路径的贡献：x＝(f₁,f₂,…,f₂₇)^T，为27维温度场向量；b＝(τ₁,τ₂,…,τ₅₈)^T，为步骤2得到的58维飞渡时间测量数据向量；

求出矩阵x＝(f₁,f₂,…,f₂₇)^T，将x＝(f₁,f₂,…,f₂₇)^T带入其中Z为声学常数，即求解出各个网格内的平均温度，将此温度作为每个网格空间几何中心的温度，再在MATLAB中调用三次样条插值方法即可得出整个待测空间区域的三维温度场。

有益效果

本发明提出的一种爆震火焰温度测量***及三维温度场的重建方法，使用声学测温具有响应速率快，不受外界壁温影响辐射等优点。在获得时间时，由于采用了互相关技术，可以在低信噪比的环境中得到较为准确的时延估计。另外，采用倒谱的方法对混响进行处理，减少了声波反射对测量的影响。在温度场重建中，采用了基于LSQR算法的区域离散法，可以实现温度场的三维构建，连续测量和实时监测。最终实现高温、高传播速度的爆震火焰瞬态温度场分析。

附图说明

图1是本发明所述的方法对爆震燃烧室温度场测量的实验示意图

图2是传感器安装位置示意图

图3时本发明提出方法的流程图

具体实施方式

现结合实施例、附图对本发明作进一步描述：

所述爆震火焰图像测温方法的步骤如下：

步骤1：爆震实验台的搭建。参照图1，针对某一台脉冲爆震发动机样机，在爆震燃烧室的某一个空间内分别安装20个麦克风和扬声器(扬声器和麦克风共同组成了声波传感器)，通过数据采集卡连接计算机；首先由主机中的声波测温软件***发出数字声波信号，经声卡转换为模拟信号后，再经功率放大器放大输出，到达继电器开关的一端，与此同时，主机上的声波测温软件***向数采卡上的数字I/O输出高电平，此高电平由与数采卡相连的端子板引出，经控制信号电缆去启动继电器，将功率放大器输出的模拟信号与扬声器相连通，使扬声器发射出声波信号。麦克风将到达的声波信号转换为电压信号，经信号调理器放大后，被数采卡采集，电压信号被主机上的声波测温软件调入。

步骤2：时间获取。

数采卡上有20根数字I/O输出线对应控制继电器开关板上的20个继电器的闭合，从而实现声波信号发射的程序控制化。将扬声器和麦克风安装在燃烧室的一个空间内，参照图2，按顺时针方向进行编号，分别为1-20号，每次只有一个扬声器发射声波信号，当1号扬声器发射声波信号时，用1号的麦克风采集的数据与其他19个编号的麦克风采集的数据进行互相关运算，对应得出19条路径上声波飞渡时间，1-20扬声器依次发射声波信号，组成一个测量周期。将测得的电压信号调入声波测温软件，再用互相关算法，可得到声波飞渡时间。计算过程为：

步骤2.1：混响处理－倒谱法。

当燃烧室内存在混响时，假设两麦克接收到的声波信号转化为的电压信号可表示为：

x₁(t)＝h₁(t)*s(t-τ₁)+n₁(t)

x₂(t)＝h₂(t)*s(t-τ₂)+n₂(t)

其中，h₁(t)和h₂(t)燃烧室的单位冲激响应；n₁(t)，n₂(t)为麦克风接收的噪声信号；s(t-τ₁)，s(t-τ₂)为两麦克风接收到的扬声器发出的信号。

对上式做倒谱变换：

$\stackrel{\‾}{X}\left(k\right)=\stackrel{\‾}{h}\left(k\right)+\stackrel{\‾}{s}\left(k\right)+\stackrel{\‾}{\mathrm{\ϵ}}\left(k\right)$

其中，分别是x(t)，h(t)，s)t)的倒谱。

单位冲激响应h(t)的傅立叶变换为H(ω)，它可以分为最小相位部分H_ap(ω)，和全通部分H_min(ω)，如下式：

H(ω)＝H_ap(ω)H_min(ω)

其中，|H_ap(ω)|＝1。

单位冲激响应h(t)的倒谱为：

$\stackrel{\‾}{h}\left(k\right)={\stackrel{\‾}{h}}_{ap}\left(k\right)+{\stackrel{\‾}{h}}_{min}\left(k\right)$

其中分别对应全通部分和最小相位部分的倒谱，由最小相位的性质可得倒谱的最小相位部分为：

${\stackrel{\&OverBar;}{h}}_{min}\left(k\right)=\left\{\begin{array}{c}0,k<0\\ \stackrel{\&OverBar;}{h}\left(k\right),k=0\\ \stackrel{\&OverBar;}{h}\left(k\right)+\stackrel{\&OverBar;}{h}(-k),k>0\end{array}\right.$

${\stackrel{\&OverBar;}{h}}_{ap}\left(k\right)=\left\{\begin{array}{c}\stackrel{\&OverBar;}{h}\left(k\right),k<0\\ \stackrel{\&OverBar;}{h}\left(0\right),k=0\\ -\stackrel{\&OverBar;}{h}(-k),k>0\end{array}\right.$

因为，输入信号x(t)的倒谱可以表示为：

$\stackrel{\&OverBar;}{X}\left(k\right)={\stackrel{\&OverBar;}{h}}_{ap}\left(k\right)+{\stackrel{\&OverBar;}{h}}_{min}\left(k\right)+\stackrel{\&OverBar;}{s}\left(k\right)+\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}\left(k\right)$

时延τ在频域转换成e^jωr的形式，它将只影响全通部分，而不对最小相位部分产生影响，从信号倒谱中减去部分得：

$X\left(k\right)={\stackrel{\&OverBar;}{h}}_{ap}\left(k\right)+\stackrel{\&OverBar;}{s}\left(k\right)+\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}\left(k\right)$

将X(k)转换至时域：

x(t)＝s(t-τ)+n(t)

步骤2.2：互相关时延估计法

将步骤2.1处理过的信号，即为x(t)，进行互相关运算：

$R_{x_1{x}_2}\left(\mathrm{\&tau;}\right)=\frac{1}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}{x}_1(t-\mathrm{\&tau;})x_2\left(t\right)dt=\frac{1}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}\&lsqb;s_1(t-\mathrm{\&tau;})+n_1(t-\mathrm{\&tau;})\&rsqb;\&lsqb;{\mathrm{As}}_2\left(t\right)+n_2\left(t\right)\&rsqb;dt$

化简可得：

$R_{x_1{x}_2}\left(\mathrm{\&tau;}\right)=R_{s_1{s}_2}\left(\mathrm{\&tau;}\right)+R_{s_1{n}_2}\left(\mathrm{\&tau;}\right)+R_{n_1{s}_2}\left(\mathrm{\&tau;}\right)+R_{n_1{n}_2}\left(\mathrm{\&tau;}\right)$

其中，为信号的互相关运算结果，分别为信号和噪声的互相关运算结果，是噪声互相关运算结果。易知为零。可简化为：

$R_{x_1{x}_2}\left(\mathrm{\&tau;}\right)=R_{s_1{s}_2}\left(\mathrm{\&tau;}\right)=\frac{1}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}{s}_1(t-\mathrm{\&tau;})s_2\left(t\right)dt=\frac{1}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}s(t+\mathrm{\&Delta;}t-\mathrm{\&tau;})s\left(t\right)dt=R_s(\mathrm{\&tau;}-\mathrm{\&Delta;}t)$

其中，Δt为时间延迟，因为R_s(τ-Δt)≤R_s(0)，当τ＝Δt时，取得最大值，则通过计算，互相关函数峰值出现的时刻即为时间延迟Δt。

相关函数最大值点对应的横坐标就是时延。

通过互相关技术可有效的抑制燃烧室噪声的干扰，在低信噪比的环境中获得较为准确的时延估计。

步骤3：实验标定。

设定D是扬声器和麦克风之间的距离，单位为m；τ是声波飞渡时间，单位为ms，T为温度，单位是摄氏度。则温度公式：

$T=\frac{D^2}{{\mathrm{Z\&tau;}}^2}\&times;{10}^6-273.16---\left(1\right)$

其中，Z为声学常数，单位m²/(s²K)，由温度公式反推出

$Z=\frac{D^2}{(T+273.16){\mathrm{\&tau;}}^2}\&times;{10}^6---\left(2\right)$

在热态实验结束，烟气温度冷却到室温时，测得室温。根据公式(2)，计算10组数据，求得Z，取平均值，得到Z为20.19m²/(s²K)。

步骤4：爆震燃烧室温度场的建立。LSQR算法重建温度场如下：

令L_i为第i条声波路径，τ_i为步骤2中测得的声波飞渡时间，有公式：

${\mathrm{\&tau;}}_i={\&Integral;}_{L_i}\frac{1}{v(x,y,z)}dl={\&Integral;}_{L_i}f(x,y,z)dl---\left(3\right)$

式中：v(x,y,z)为待测区域内声波传播速度，f(x,y,z)为v(x,y,z)的倒数。

爆震燃烧室尺寸一般为10×10×100cm，选取从燃烧室出口到10cm处的一个测量空间(10×10×10cm)。声学测温空间中有58条声波路径(除去重合和壁面上的)，将三维燃烧室测量空间均分为3×3×3＝27个空间网格，则根据公式(3)可得：

${\mathrm{\&tau;}}_i=\underset{j=1}{\overset{N}{\&Sigma;}}{w}_{ij}{f}_j,i=1,2...,58$

式中：w_ij为第i条路径通过第j个空间网格的长度。

经过一个周期的测量，可以得到一个线性方程组：

$\begin{array}{c}{w}_{11}{f}_1+w_{12}{f}_2+\mathrm{...}+w_{1,27}{f}_{27}={\mathrm{\&tau;}}_1\\ w_{21}{f}_1+w_{22}{f}_2+\mathrm{...}+w_{2,27}{f}_{27}={\mathrm{\&tau;}}_2\\ .\\ .\\ .\\ w_{58,1}{f}_1+w_{58,2}{f}_2+\mathrm{...}{w}_{58,27}{f}_{27}={\mathrm{\&tau;}}_{58}\end{array}$

该方程组用矩阵可表示为

Ax＝b，即为三维温度场重建的数学模型

式中：为58×27维权因子矩阵，反映了第j空间网格对第i条声波路径的贡献：x＝(f₁,f₂,…,f₂₇)^T，x为27维速度倒数向量：b＝(τ₁,τ₂,…,τ₅₈)^T，为58维飞渡时间测量数据向量。

求出矩阵x，即将其带入公式(1)，利用步骤三得出得Z,可得温度即求解出各个网格内的平均温度，将此温度作为每个网格空间几何中心的温度，再在MATLAB中调用三次样条插值方法即可得出整个待测空间区域的三维温度场。

Claims (3)

1.一种爆震火焰温度测量***，其特征在于包括带声卡和显示器的主机、功率放大器、继电器、数采卡、端子板、信号调理器、扬声器和麦克风；扬声器和麦克风位于在爆震燃烧室内，带声卡和显示器的主机、功率放大器、继电器、数采卡、端子板和信号调理器位于爆震燃烧室外，扬声器与继电器连接，麦克风与信号调理器连接；主机中的声波测温软件***发出数字声波信号，经声卡转换为模拟信号后，再经功率放大器放大输出到继电器开关的一端；与此同时，主机上的声波测温软件***向数采卡上的数字I/O输出高电平，高电平由与数采卡相连的端子板引出，经信号控制电缆去启动继电器，将功率放大器输出的模拟信号传输到扬声器，扬声器将模拟信号转化为声波信号，麦克风将接收到的声波信号转换为电压信号，经信号调理器放大后，被数采卡采集，然后输入到声波测温软件***进行三维温度场的重建。

2.根据权利要求1所述的一种爆震火焰温度测量***，其特征在于扬声器和麦克风的数量分别为20个，一个扬声器和一个麦克风组成一组，20组均匀的分布在立方体空间的各个边上。

3.一种对权利要求2所述的爆震火焰温度测量***测得的电压信号进行三维温度场的重建方法，其特征在于步骤如下：

步骤1：采用倒谱法对两麦克接收到的电压信号进行倒谱变换：

其中，接收到的电压信号为：

x₁(t)＝h₁(t)*s(t-τ₁)+n₁(t)

x₂(t)＝h₂(t)*s(t-τ₂)+n₂(t)

其中，h₁(t)、h₂(t)为燃烧室的单位冲激响应；n₁(t)、n₂(t)为麦克风接收的噪声信号；s(t-τ₁)、s(t-τ₂)为两麦克风接收扬声器发出的信号；

进行倒谱变换后的电压信号为：

x₁(t)＝s(t-τ₁)+n₁(t)

x₂(t)＝s(t-τ₂)+n₂(t)

步骤2：对倒谱变换后的x₁(t)和x₂(t)进行互相关运算：

$R_{x_1{x}_2}\left(\mathrm{\&tau;}\right)=\frac{1}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}{x}_1(t-\mathrm{\&tau;})x_2\left(t\right)\mathrm{\&omega;}=\frac{1}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}\&lsqb;s_1(t-\mathrm{\&tau;})+n_1(t-\mathrm{\&tau;})\&rsqb;\&lsqb;{\mathrm{As}}_2\left(t\right)+n_2\left(t\right)\&rsqb;dt$

化简可得：

$R_{x_1{x}_2}\left(\mathrm{\&tau;}\right)=R_{s_1{s}_2}\left(\mathrm{\&tau;}\right)+R_{s_1{n}_2}\left(\mathrm{\&tau;}\right)+R_{n_1{s}_2}\left(\mathrm{\&tau;}\right)+R_{n_1{n}_2}\left(\mathrm{\&tau;}\right)$

其中，A为声波传播过程中的衰减系数，为信号的互相关运算结果，分别为信号和噪声的互相关运算结果，是噪声互相关运算结果，易知为零，可简化为：

$R_{x_1{x}_2}\left(\mathrm{\&tau;}\right)=R_{s_1{s}_2}\left(\mathrm{\&tau;}\right)=\frac{1}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}{s}_1(t-\mathrm{\&tau;})s_2\left(t\right)dt=\frac{1}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}s(t+\mathrm{\&Delta;}t-\mathrm{\&tau;})s\left(t\right)dt=R_s(\mathrm{\&tau;}-\mathrm{\&Delta;}t)$

其中，Δt为时间延迟，因为R_s(τ-Δt)≤R_s(0)，当τ＝Δt时，取得最大值，则通过计算，互相关函数峰值出现的时刻即为时间延迟Δt即飞渡时间τ；

步骤3：采用LSQR算法重建温度场：

将三维燃烧室测量空间均分为3×3×3＝27个空间网格，则

${\mathrm{\&tau;}}_i=\underset{j=1}{\overset{N}{\&Sigma;}}{w}_{ij}{f}_j,i=1,2...,58$

式中：w_ij为第i条路径通过第j个空间网格的长度；τ_i为测得的第i条声波路径的声波飞渡时间；f_j为待测区域内声波传播速度的倒数；

控制20个扬声器依次发射声波信号，组成一个测量周期，则得到58条声波路径，相应的线性方程组为：

w₁₁f₁+w₁₂f₂+…w_1,27f₂₇＝τ₁

w₂₁f₁+w₂₂f₂+…w_2，27f₂₇＝τ₂

w_58，1f₁+w_58，2f₂+…w_58，27f₂₇＝τ₅₈

该方程组用矩阵表示为

Ax＝b，即为三维温度场重建的数学模型

式中：为58×27维权因子矩阵，反映了第j空间网格对第i条声波路径的贡献：x＝(f₁,f₂,...,f₂₇)^T，为27维温度场向量；b＝(τ₁,τ₂,…,τ₅₈)^T，为步骤2得到的58维飞渡时间测量数据向量；

求出矩阵x＝(f₁,f₂,…,f₂₇)^T，将x＝(f₁,f₂,…,f₂₇)^T带入其中Z为声学常数，即求解出各个网格内的平均温度，将此温度作为每个网格空间几何中心的温度，再在MATLAB中调用三次样条插值方法即可得出整个待测空间区域的三维温度场。