CN102353478B - 半透明介质环境下非接触测温的校正方法 - Google Patents

Abstract

半透明介质环境下非接触测温的校正方法，属于高温测量技术领域。它解决了被测物体表面处于半透明介质覆盖的环境下时，采用传统方法探测到的辐射能量不能通过传统的材料表面发射率修正方法进行修正得到其真实温度的问题。首先判断半透明介质与被测材料表面是否接触，若接触，选择一维耦合换热模型，采用有限体积法进行正向模型的计算，得到测温设备能够获得的理论辐射能量值；若不接触，选择一维纯辐射换热模型，采用有限体积法进行正向模型的计算，得到测温设备能够获得的理论辐射能量值；然后测量被测材料表面的实际辐射能量值；采用智能微粒群优化算法反演被测材料表面的真实温度值。本发明适用于半透明介质环境下被测材料表面的温度测量。

Description

半透明介质环境下非接触测温的校正方法

技术领域

本发明涉及一种半透明介质环境下非接触测温的校正方法，属于高温测量技术领域。

背景技术

温度是确定物质状态的最重要参数之一，对温度进行测量与控制在国防、军事、科学实验及工农业生产中具有十分重要的作用，尤其是高温测量，在航天、材料、能源、冶金等领域占有极重要的地位。

温度测量大致可分为接触法测量和非接触法测量两大类。接触法测温包括热电偶测温和热电阻测温等，非接触法测温主要是以辐射测温为主。近二十年内，随着电子技术及计算机技术的飞速发展，辐射测温技术得到了长足的进步与发展。辐射测温具有无测量上限、响应速度快及不接触被测对象因而不影响测温场等优点，目前辐射温度计已经发展到使用硅电二极管作为检测器的光学测量和光电精密测量阶段。辐射法测量物体真温是各国学者一直关心的问题，已提出的辐射测温法如发射率修正法、逼近黑体法、反射率测量法及多光谱辐射测温法等。

对固体表面进行红外辐射测温的传统方法，主要是采用红外辐射热像仪得到两光谱(比色法)或多光谱下的辐射能量，再结合已知的固体材料表面发射率或假设表面发射率分布函数，通过最小二乘近似等方法得到其真实温度。但对于半透明材料而言，采用热像仪进行红外辐射测温时，由于辐射的沿程性，热像仪探测到的辐射能量来自半透明材料内部沿探测方向各点的辐射出射能量之和，这个辐射出射能量之和与材料内部的温度、物性等有关，不能通过传统的材料表面发射率修正方法得到其真实温度，而需要提供一种方法，综合考虑辐射物性的影响，再结合反问题算法反演来得到其真实温度。

发明内容

本发明的目的是解决被测物体表面处于半透明介质覆盖的环境下时，采用传统方法探测到的辐射能量不能通过传统的材料表面发射率修正方法进行修正得到其真实温度的问题，提供一种半透明介质环境下非接触测温的校正方法。

本发明所述半透明介质环境下非接触测温的校正方法，它包括以下步骤：

步骤一：判断半透明介质与被测材料表面是否接触，若接触，执行步骤二；否则，执行步骤三；

步骤二：选择一维耦合换热模型，采用有限体积法进行正向模型的计算，得到测温设备能够获得的理论辐射能量值，然后执行步骤四；

步骤三：选择一维纯辐射换热模型，采用有限体积法进行正向模型的计算，得到测温设备能够获得的理论辐射能量值，然后执行步骤四；

步骤四：采用测温设备实际测量被测材料表面，得到测温设备获得的实际辐射能量值；

步骤五：根据所述理论辐射能量值和实际辐射能量值，采用智能微粒群优化算法反演被测材料表面的真实温度值。

步骤二中获得测量设备的理论辐射能量值的具体方法为：

选择一维耦合换热模型，根据计算精度要求将半透明介质内部沿与被测材料表面垂直的方向进行网格划分，均分为多个网格单元，并采用天顶角和水平角均匀划分的方法将半透明介质内部空间立体角离散划分为N_Ω份，所述网格单元与被测材料表面平行；设定辐射源项q_R的初始值为0，利用该一维耦合换热模型的能量守恒方程和耦合换热边界条件求出半透明介质内部待求网格单元P中心节点的温度T_p：

能量守恒方程： $\mathrm{\λ}\frac{{\∂}^{2}T}{{\∂x}^2}-q_R=0,---\left(1\right)$

耦合换热边界条件： $\left\{\begin{array}{cc}x=0,& T=T_w\\ x=L,& \mathrm{\λ}\frac{\∂T}{\∂x}=h_f(T_w^{\′}-T_f)\end{array}\right.,---\left(2\right)$

采用有限体积法离散上述方程后，求解获得离散后半透明介质内部待求网格单元P的中心点温度T_P：

$T_P=(T_Z+T_Y-\frac{q_{R,P}\·{\mathrm{\Δx}}^2}{\mathrm{\λ}})/2;---\left(3\right)$

式中：λ表示半透明介质的导热系数，T表示半透明介质内部温度，x表示半透明介质水平坐标，T_w表示半透明介质与被测材料表面的边界温度，T′_w表示半透明介质环境表面的边界温度，T_f表示半透明介质环境表面侧的环境温度，h_f表示半透明介质的环境表面的对流换热系数，L表示半透明介质厚度，T_Y表示与待求网格单元P相邻的一侧的网格单元Y的中心点温度，T_Z表示与待求网格单元P相邻的另一侧的网格单元Z的中心点温度，q_R，P为离散后的半透明介质内部网格单元的辐射源项，Δx表示相邻网格单元的中心点之间的距离；

根据半透明介质内部温度T，利用半透明介质辐射传输模型中的辐射传输方程和边界条件采用有限体积法求得半透明介质内部网格单元在k谱带内方向辐射强度

所述辐射传输方程为：

$\frac{{\mathrm{dI}}_k^m\left(s\right)}{\mathrm{ds}}=-{\mathrm{\κ}}_{e,k}{I}_k^m\left(s\right)+{\mathrm{\κ}}_{a,k}{I}_{b,k}^m\left(s\right)+\frac{{\mathrm{\κ}}_{s,k}}{4\mathrm{\π}}{\∫}_{{\mathrm{\Ω}}^{m^{\′}}=4\mathrm{\π}}{I}_k^{m^{\′}}\left(s\right){\mathrm{\Φ}}_k({\mathrm{\Ω}}^{m^{\′}},{\mathrm{\Ω}}^m){\mathrm{d\Ω}}^{m^{\′}},---\left(4\right)$

所述边界条件：

$I_{w,k}^m={\left(\frac{n}{n_0}\right)}^2(1-{\mathrm{\&rho;}}_{0,k})I_{0,k}^m+\frac{f_0{\mathrm{\&rho;}}_k}{\mathrm{\&pi;}}{\&Integral;}_{n_w\&CenterDot;{\mathrm{\&Omega;}}^{m^{\&prime;}}<0}{I}_{w,k}^{m^{\&prime;}}|n_w\&CenterDot;{\mathrm{\&Omega;}}^{m^{\&prime;}}|{\mathrm{d\&Omega;}}^{m^{\&prime;}}+(1-f_0)I_{w,k}^{m^{\&prime;\&prime;}},---\left(5\right)$

采用有限体积法求得半透明介质内部待求网格单元P在k谱带内方向辐射强度

为：

$I_{k,P}^m=(a_{k,Y}^m{I}_{k,Y}^m+a_{k,Z}^m{I}_{k,Z}^m+c_{k,P}^m)/a_{k,P}^m,---\left(6\right)$

其中： $a_{k,P}^m=\underset{j=z,y}{\mathrm{\&Sigma;}}\max (A_j{D}_j^m,0)+{\mathrm{\&kappa;}}_{e,k,P}{V}_P{\mathrm{\&Delta;\&Omega;}}^m,---(6-1)$

$a_{k,J}^m=\max ({-A}_j{D}_j^m,0)(J=Z,Y,j=z,y),---(6-2)$

$C_{k,P}^m={\mathrm{\&kappa;}}_{e,k,P}\&CenterDot;U_{k,P}^m\&CenterDot;V_P\&CenterDot;{\mathrm{\&Delta;\&Omega;}}^m,---(6-3)$

$U_{k,P}^m=\frac{{\mathrm{\&kappa;}}_{a,k,P}}{{\mathrm{\&kappa;}}_{e,k,P}}\frac{{\mathrm{\&sigma;T}}_P^4}{\mathrm{\&pi;}}+\frac{{\mathrm{\&kappa;}}_{s,k,P}}{{\mathrm{\&kappa;}}_{e,k,P}}\&CenterDot;\frac{1}{4\mathrm{\&pi;}}\underset{m^{\&prime;}}{\mathrm{\&Sigma;}}{I}_{k,P}{\mathrm{\&Phi;}}_k({\mathrm{\&Omega;}}^m,{\mathrm{\&Omega;}}^{m^{\&prime;}}){\mathrm{\&Delta;\&Omega;}}^{m^{\&prime;}},---(6-4)$

式中，待求网格单元P在k谱带内的方向辐射强度中，m表示第m个立体角方向Ω^m，k表示谱带，

表示在k谱带内Ω^m方向的半透明介质壁面辐射强度，

表示在k谱带内Ω^m′方向的半透明介质壁面辐射强度，

表示在k谱带内方向黑体辐射强度，

表示与离散后的待求网格单元P相邻的靠近被测材料表面侧的网格单元Z在k谱带内的Ω^m方向辐射强度，表示与离散后的待求网格单元P相邻的靠近环境表面侧的网格单元Y在k谱带内的Ω^m方向辐射强度；

κ_e，k表示谱带衰减系数，κ_e，k表示谱带吸收系数，κ_s，k表示谱带散射系数，κ_e，k，P表示离散的待求网格单元P的谱带衰减系数；

表示线性方程组的系数、

表示线性方程组的系数、

表示线性方程组的系数，

表示线性方程组的常数项，

表示线性方程组常数项

的一部分，为中间变量，Ω^m′表示第m′个立体角方向，ΔΩ^m′表示第m′个立体角的大小，Ω^m″表示Ω^m的镜反射方向，Φ_k表示k谱带内Ω^m′方向在Ω^m方向上的介质散射相函数，其中Ω^m′表示除了Ω^m以外的其他立体角方向；

n表示半透明介质的折射率，n₀表示环境折射率，ρ_0，k表示环境的谱带反射率，ρ_k表示半透明介质的谱带反射率，f₀表示半透明介质壁面上的漫反射占总反射能量的比例，n_w表示半透明介质壁面的法向量，

j表示半透明介质待求网格单元的表面序号，A_j表示半透明介质第j个待求网格单元表面的面积，

表示j表面法向在Ω^m方向的权重，V_P表示离散后的待求网格单元的体积，σ表示斯蒂芬-波尔兹曼常数：σ＝5.67×10^-8[W/(m²·K⁴)]；

根据半透明介质内部待求网格单元P的温度T_P和半透明介质内部网格单元在k谱带内方向辐射强度

计算出半透明介质内部待求网格单元P的辐射源项q_R，P：

$q_{R,P}=\underset{k=1}{\overset{M_b}{\mathrm{\&Sigma;}}}{\mathrm{\&kappa;}}_{a,k,P}[{4B}_{k,T_P}{\mathrm{\&sigma;T}}_P^4-\underset{m=1}{\overset{N_{\mathrm{\&Omega;}}}{\mathrm{\&Sigma;}}}{I}_{k,P}^m{\mathrm{\&Omega;}}^m];---\left(7\right)$

式中M_b表示谱带份数，κ_a，k，P表示待求网格单元P在k谱带内的吸收系数，

表示在温度T_P下，谱带模型k谱带内的辐射能量占总辐射能量的比例；

将半透明介质内部待求网格单元P的辐射源项q_R，P代入方程(1)，重复前述过程，直至求得收敛的半透明介质内部每个网格单元的温度T_P，即前后两次迭代的相对误差小于预设精度阈值，再根据收敛的半透明介质内部每个网格单元的温度T_P计算出边界出射辐射热流密度q_w，P：

$q_{w,P}=\underset{k=1}{\overset{M_b}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_{w,k}[B_{k,T_w}{\mathrm{\&sigma;T}}_w^4-\underset{m=1}{\overset{N_{\mathrm{\&Omega;}}/2}{\mathrm{\&Sigma;}}}{I}_{w,k}^m{D}_w^m],---\left(8\right)$

式中ε_w，k表示半透明介质壁面的谱带发射率，

表示在半透明介质与被测材料表面的边界温度T_w下，谱带模型k谱带内的辐射能量占总辐射能量的比例，

表示半透明介质壁面法向在Ω^m方向的权重；

最后根据下述公式计算获得测温设备能够获得的理论辐射能量值Q，完成正算：

$Q={\mathrm{\&tau;q}}_{w,P}{A}_P+(1-\mathrm{\&tau;}+\mathrm{\&rho;\&tau;})\underset{k=1}{\overset{M_n}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_k{\mathrm{\&sigma;T}}_g^4,---\left(9\right)$

式中，τ表示环境大气透过率，A_P表示离散后的半透明介质边界单元面积，ε_k表示环境气体光谱发射率，T_g表示环境大气温度。

步骤三中获得测温设备的理论辐射能量值的具体方法为：

假设被测材料表面的温度为T_w，根据被测材料表面的光谱发射率ε_k计算被测材料表面的出射辐射强度

$I_{0,k}^m=\frac{{\mathrm{\&epsiv;}}_k{\mathrm{\&sigma;T}}_w^4}{\mathrm{\&pi;}},---\left(10\right)$

将出射辐射强度

作为半透明介质一维纯辐射换热模型的边界条件，根据一维纯辐射换热模型，采用有限体积法求得半透明介质内部的方向谱带辐射强度

辐射传输方程为：

$\frac{{\mathrm{dI}}_k^m\left(s\right)}{\mathrm{ds}}=-{\mathrm{\&kappa;}}_{e,k}{I}_k^m\left(s\right)+{\mathrm{\&kappa;}}_{a,k}{I}_{b,k}^m\left(s\right)+\frac{{\mathrm{\&kappa;}}_{s,k}}{4\mathrm{\&pi;}}{\&Integral;}_{{\mathrm{\&Omega;}}^{m^{\&prime;}}=4\mathrm{\&pi;}}{I}_k^{m^{\&prime;}}\left(s\right){\mathrm{\&Phi;}}_k({\mathrm{\&Omega;}}^{m^{\&prime;}},{\mathrm{\&Omega;}}^m){\mathrm{d\&Omega;}}^{m^{\&prime;}},---\left(11\right)$

边界条件为：

$I_{w,k}^m={\left(\frac{n}{n_0}\right)}^2(1-{\mathrm{\&rho;}}_{0,k})I_{0,k}^m+\frac{f_0{\mathrm{\&rho;}}_k}{\mathrm{\&pi;}}{\&Integral;}_{n_w\&CenterDot;{\mathrm{\&Omega;}}^{m^{\&prime;}}<0}{I}_{w,k}^{m^{\&prime;}}|n_w\&CenterDot;{\mathrm{\&Omega;}}^{m^{\&prime;}}|{\mathrm{d\&Omega;}}^{m^{\&prime;}}+(1-f_0)I_{w,k}^{m^{\&prime;\&prime;}},---\left(12\right)$

采用有限体积法求解半透明介质内部的方向谱带辐射强度

$I_{k,P}^m=(a_{k,Y}^m{I}_{k,Y}^m+a_{k,Z}^m{I}_{k,Z}^m+c_{k,P}^m)/a_{k,P}^m,---\left(13\right)$

其中： $a_{k,P}^m=\underset{j=z,y}{\mathrm{\&Sigma;}}\max (A_j{D}_j^m,0)+{\mathrm{\&kappa;}}_{e,k,P}{V}_P{\mathrm{\&Delta;\&Omega;}}^m,---(13-1)$

$a_{k,J}^m=\max ({-A}_j{D}_j^m,0)(J=Z,Y,j=z,y),---(13-2)$

$C_{k,P}^m={\mathrm{\&kappa;}}_{e,k,P}\&CenterDot;U_{k,P}^m\&CenterDot;V_P\&CenterDot;{\mathrm{\&Delta;\&Omega;}}^m,---(13-3)$

$U_{k,P}^m=\frac{{\mathrm{\&kappa;}}_{a,k,P}}{{\mathrm{\&kappa;}}_{e,k,P}}\frac{{\mathrm{\&sigma;T}}_P^4}{\mathrm{\&pi;}}+\frac{{\mathrm{\&kappa;}}_{s,k,P}}{{\mathrm{\&kappa;}}_{e,k,P}}\&CenterDot;\frac{1}{4\mathrm{\&pi;}}\underset{m^{\&prime;}}{\mathrm{\&Sigma;}}{I}_{k,P}{\mathrm{\&Phi;}}_k({\mathrm{\&Omega;}}^m,{\mathrm{\&Omega;}}^{m^{\&prime;}}){\mathrm{\&Delta;\&Omega;}}^{m^{\&prime;}},---(13-4)$

重复前述过程，直至求得收敛的半透明介质内部的方向谱带辐射强度

即前后两次迭代的相对误差小于预设精度阈值；

将半透明介质内部沿与被测材料表面垂直的方向进行网格划分，均分为多个网格单元，所述网格单元与被测材料表面平行；根据半透明介质内部每个网格单元的中心点温度TP计算出边界出射辐射热流密度q_w，P：

$q_{w,P}=\underset{k=1}{\overset{M_b}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_{w,k}[B_{k,T_w}{\mathrm{\&sigma;T}}_w^4-\underset{m=1}{\overset{N_{\mathrm{\&Omega;}}/2}{\mathrm{\&Sigma;}}}{I}_{w,k}^m{D}_w^m],---\left(14\right)$

最后根据下述公式计算获得测温设备能够获得的理论辐射能量值Q，完成正算：

$Q={\mathrm{\&tau;q}}_{w,P}{A}_P+(1-\mathrm{\&tau;}+\mathrm{\&rho;\&tau;})\underset{k=1}{\overset{M_n}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_k{\mathrm{\&sigma;T}}_g^4.---\left(15\right)$

本发明的优点是：本发明方法提供了半透明介质环境下非接触测温的校正技术，由于采用非接触式辐射测温方法测量半透明介质环境下的被测材料表面温度，辐射测温设备会受到半透明介质本身辐射的影响，因此会出现测量结果不准确的问题。本发明方法将半透明介质与被测材料表面是否接触两种状况，分别建立了正算模型，同时结合智能微粒群算法PSO反演被测材料表面的温度，同时能够获得半透明介质内部的温度分布情况，通过正算和反算的循环得到更为准确的温度值，实现了对传统辐射测温结果的校正。当被测材料表面受到半透明参与性介质本身的辐射影响时，通过测温设备测量获得的热流密度值，采用PSO算法结合正算模型能够准确反演出被测材料表面温度，同时可获得半透明介质内部的温度分布。

本发明方法适用于航天、材料、能源及冶金等高温测量的工程领域，对半透明介质环境下的非接触红外测温技术具有直接的理论指导意义。

附图说明

图1为将半透明介质均分为多个网格单元的局部示意图；

图2为半透明介质与被测材料表面不接触的辐射能量传输示意图；

图3为半透明介质与被测材料表面接触的辐射能量传输示意图。

具体实施方式

具体实施方式一：下面结合图1至图3说明本实施方式，本实施方式所述半透明介质环境下非接触测温的校正方法，它包括以下步骤：

步骤一：判断半透明介质与被测材料表面是否接触，若接触，执行步骤二；否则，执行步骤三；

步骤二：选择一维耦合换热模型，采用有限体积法进行正向模型的计算，得到测温设备能够获得的理论辐射能量值，然后执行步骤四；

步骤三：选择一维纯辐射换热模型，采用有限体积法进行正向模型的计算，得到测温设备能够获得的理论辐射能量值，然后执行步骤四；

步骤四：采用测温设备实际测量被测材料表面，得到测温设备获得的实际辐射能量值；

步骤五：根据所述理论辐射能量值和实际辐射能量值，采用智能微粒群优化算法反演被测材料表面的真实温度值。

本实施方式中所述的测温设备指能测出其所接受能量的辐射测温设备，测温设备需要沿着被测材料表面的法线方向进行测量，假设测温表面是一个平面，这样就可将计算模型转化为一维无限大平板模型。当半透明介质与被测材料表面接触时，不能忽略导热效应，这时正算模型为一个一维耦合换热模型，在准稳态时用辐射测温设备测量出接收到的能量；如果半透明介质不与被测材料表面接触，则正算模型为纯辐射传输问题，在半透明介质遮蔽被测材料表面时采用辐射测温设备测出接收到的能量，通过正算和反演得出更准确的温度值。

具体实施方式二：下面结合图1和图3说明本实施方式，本实施方式为对实施方式一的进一步说明，步骤二中获得测量设备的理论辐射能量值的具体方法为：

选择一维耦合换热模型，根据计算精度要求将半透明介质内部沿与被测材料表面垂直的方向进行网格划分，均分为多个网格单元，并采用天顶角和水平角均匀划分的方法将半透明介质内部空间立体角离散划分为N_Ω份，所述网格单元与被测材料表面平行；设定辐射源项q_R的初始值为0，利用该一维耦合换热模型的能量守恒方程和耦合换热边界条件求出半透明介质内部待求网格单元P中心节点的温度T_p：

能量守恒方程： $\mathrm{\&lambda;}\frac{{\&PartialD;}^{2}T}{{\&PartialD;x}^2}-q_R=0,---\left(1\right)$

耦合换热边界条件： $\left\{\begin{array}{cc}x=0,& T=T_w\\ x=L,& \mathrm{\&lambda;}\frac{\&PartialD;T}{\&PartialD;x}=h_f(T_w^{\&prime;}-T_f)\end{array}\right.,---\left(2\right)$

采用有限体积法离散上述方程后，求解获得离散后半透明介质内部待求网格单元P的中心点温度T_P：

$T_P=(T_Z+T_Y-\frac{q_{R,P}\&CenterDot;{\mathrm{\&Delta;x}}^2}{\mathrm{\&lambda;}})/2;---\left(3\right)$

式中：λ表示半透明介质的导热系数，T表示半透明介质内部温度，x表示半透明介质水平坐标，T_w表示半透明介质与被测材料表面的边界温度，T′_w表示半透明介质环境表面的边界温度，T_f表示半透明介质环境表面侧的环境温度，h_f表示半透明介质的环境表面的对流换热系数，L表示半透明介质厚度，T_Y表示与待求网格单元P相邻的一侧的网格单元Y的中心点温度，T_Z表示与待求网格单元P相邻的另一侧的网格单元Z的中心点温度，q_R，P为离散后的半透明介质内部网格单元的辐射源项，Δx表示相邻网格单元的中心点之间的距离；

根据半透明介质内部温度T，利用半透明介质辐射传输模型中的辐射传输方程和边界条件采用有限体积法求得半透明介质内部网格单元在k谱带内方向辐射强度

所述辐射传输方程为：

$\frac{{\mathrm{dI}}_k^m\left(s\right)}{\mathrm{ds}}=-{\mathrm{\&kappa;}}_{e,k}{I}_k^m\left(s\right)+{\mathrm{\&kappa;}}_{a,k}{I}_{b,k}^m\left(s\right)+\frac{{\mathrm{\&kappa;}}_{s,k}}{4\mathrm{\&pi;}}{\&Integral;}_{{\mathrm{\&Omega;}}^{m^{\&prime;}}=4\mathrm{\&pi;}}{I}_k^{m^{\&prime;}}\left(s\right){\mathrm{\&Phi;}}_k({\mathrm{\&Omega;}}^{m^{\&prime;}},{\mathrm{\&Omega;}}^m){\mathrm{d\&Omega;}}^{m^{\&prime;}},---\left(4\right)$

所述边界条件：

$I_{w,k}^m={\left(\frac{n}{n_0}\right)}^2(1-{\mathrm{\&rho;}}_{0,k})I_{0,k}^m+\frac{f_0{\mathrm{\&rho;}}_k}{\mathrm{\&pi;}}{\&Integral;}_{n_w\&CenterDot;{\mathrm{\&Omega;}}^{m^{\&prime;}}<0}{I}_{w,k}^{m^{\&prime;}}|n_w\&CenterDot;{\mathrm{\&Omega;}}^{m^{\&prime;}}|{\mathrm{d\&Omega;}}^{m^{\&prime;}}+(1-f_0)I_{w,k}^{m^{\&prime;\&prime;}},---\left(5\right)$

采用有限体积法求得半透明介质内部待求网格单元P在k谱带内方向辐射强度

为：

$I_{k,P}^m=(a_{k,Y}^m{I}_{k,Y}^m+a_{k,Z}^m{I}_{k,Z}^m+c_{k,P}^m)/a_{k,P}^m,---\left(6\right)$

其中： $a_{k,P}^m=\underset{j=z,y}{\mathrm{\&Sigma;}}\max (A_j{D}_j^m,0)+{\mathrm{\&kappa;}}_{e,k,P}{V}_P{\mathrm{\&Delta;\&Omega;}}^m,---(6-1)$

$a_{k,J}^m=\max ({-A}_j{D}_j^m,0)(J=Z,Y,j=z,y),---(6-2)$

$C_{k,P}^m={\mathrm{\&kappa;}}_{e,k,P}\&CenterDot;U_{k,P}^m\&CenterDot;V_P\&CenterDot;{\mathrm{\&Delta;\&Omega;}}^m,---(6-3)$

$U_{k,P}^m=\frac{{\mathrm{\&kappa;}}_{a,k,P}}{{\mathrm{\&kappa;}}_{e,k,P}}\frac{{\mathrm{\&sigma;T}}_P^4}{\mathrm{\&pi;}}+\frac{{\mathrm{\&kappa;}}_{s,k,P}}{{\mathrm{\&kappa;}}_{e,k,P}}\&CenterDot;\frac{1}{4\mathrm{\&pi;}}\underset{m^{\&prime;}}{\mathrm{\&Sigma;}}{I}_{k,P}{\mathrm{\&Phi;}}_k({\mathrm{\&Omega;}}^m,{\mathrm{\&Omega;}}^{m^{\&prime;}}){\mathrm{\&Delta;\&Omega;}}^{m^{\&prime;}},---(6-4)$

式中，待求网格单元P在k谱带内的方向辐射强度

中，m表示第m个立体角方向Ω^m，k表示谱带，

表示在k谱带内Ω^m方向的半透明介质壁面辐射强度，

表示在k谱带内Ω^m′方向的半透明介质壁面辐射强度，

表示在k谱带内方向黑体辐射强度，

表示与离散后的待求网格单元P相邻的靠近被测材料表面侧的网格单元Z在k谱带内的Ω^m方向辐射强度，

表示与离散后的待求网格单元P相邻的靠近环境表面侧的网格单元Y在k谱带内的Ω^m方向辐射强度；

κ_e，k表示谱带衰减系数，κ_a，k表示谱带吸收系数，κ_s，k表示谱带散射系数，κ_e，k，P表示离散的待求网格单元P的谱带衰减系数；

表示线性方程组的系数、

表示线性方程组的系数、

表示线性方程组的系数，

表示线性方程组的常数项，

表示线性方程组常数项

的一部分，为中间变量，Ω^m′表示第m′个立体角方向，ΔΩ^m′表示第m′个立体角的大小，Ω^m″表示Ω^m的镜反射方向，Φ_k表示k谱带内Ω^m′方向在Ω^m方向上的介质散射相函数，其中Ω^m′表示除了Ω^m以外的其他立体角方向；

n表示半透明介质的折射率，n₀表示环境折射率，ρ_0，k表示环境的谱带反射率，ρ_k表示半透明介质的谱带反射率，f₀表示半透明介质壁面上的漫反射占总反射能量的比例，n_w表示半透明介质壁面的法向量，

j表示半透明介质待求网格单元的表面序号，A_j表示半透明介质第j个待求网格单元表面的面积，

表示j表面法向在Ω^m方向的权重，V_P表示离散后的待求网格单元的体积，σ表示斯蒂芬-波尔兹曼常数：σ＝5.67×10^-8[W/(m²·K⁴)]；

根据半透明介质内部待求网格单元P的温度T_P和半透明介质内部网格单元在k谱带内方向辐射强度计算出半透明介质内部待求网格单元P的辐射源项q_R，P：

$q_{R,P}=\underset{k=1}{\overset{M_b}{\mathrm{\&Sigma;}}}{\mathrm{\&kappa;}}_{a,k,P}[{4B}_{k,T_P}{\mathrm{\&sigma;T}}_P^4-\underset{m=1}{\overset{N_{\mathrm{\&Omega;}}}{\mathrm{\&Sigma;}}}{I}_{k,P}^m{\mathrm{\&Omega;}}^m];---\left(7\right)$

式中M_b表示谱带份数，κ_a，k，P表示待求网格单元P在k谱带内的吸收系数，表示在温度T_P下，谱带模型k谱带内的辐射能量占总辐射能量的比例；

将半透明介质内部待求网格单元P的辐射源项q_R，P代入方程(1)，重复前述过程，直至求得收敛的半透明介质内部每个网格单元的温度T_P，即前后两次迭代的相对误差小于预设精度阈值，再根据收敛的半透明介质内部每个网格单元的温度T_P计算出边界出射辐射热流密度q_w，P：

$q_{w,P}=\underset{k=1}{\overset{M_b}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_{w,k}[B_{k,T_w}{\mathrm{\&sigma;T}}_w^4-\underset{m=1}{\overset{N_{\mathrm{\&Omega;}}/2}{\mathrm{\&Sigma;}}}{I}_{w,k}^m{D}_w^m],---\left(8\right)$

式中ε_w，k表示半透明介质壁面的谱带发射率，

表示在半透明介质与被测材料表面的边界温度T_w下，谱带模型k谱带内的辐射能量占总辐射能量的比例，

表示半透明介质壁面法向在Ω^m方向的权重；

最后根据下述公式计算获得测温设备能够获得的理论辐射能量值Q，完成正算：

$Q={\mathrm{\&tau;q}}_{w,P}{A}_P+(1-\mathrm{\&tau;}+\mathrm{\&rho;\&tau;})\underset{k=1}{\overset{M_n}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_k{\mathrm{\&sigma;T}}_g^4,---\left(9\right)$

式中，τ表示环境大气透过率，A_P表示离散后的半透明介质边界单元面积，ε_k表示环境气体光谱发射率，T_g表示环境大气温度。

具体实施方式三：下面结合图1和图2说明本实施方式，本实施方式为对实施方式一的进一步说明，步骤三中获得测温设备的理论辐射能量值的具体方法为：

假设被测材料表面的温度为T_w，根据被测材料表面的光谱发射率ε_k计算被测材料表面的出射辐射强度

$I_{0,k}^m=\frac{{\mathrm{\&epsiv;}}_k{\mathrm{\&sigma;T}}_w^4}{\mathrm{\&pi;}},---\left(10\right)$

将出射辐射强度

作为半透明介质一维纯辐射换热模型的边界条件，根据一维纯辐射换热模型，采用有限体积法求得半透明介质内部的方向谱带辐射强度

辐射传输方程为：

$\frac{{\mathrm{dI}}_k^m\left(s\right)}{\mathrm{ds}}=-{\mathrm{\&kappa;}}_{e,k}{I}_k^m\left(s\right)+{\mathrm{\&kappa;}}_{a,k}{I}_{b,k}^m\left(s\right)+\frac{{\mathrm{\&kappa;}}_{s,k}}{4\mathrm{\&pi;}}{\&Integral;}_{{\mathrm{\&Omega;}}^{m^{\&prime;}}=4\mathrm{\&pi;}}{I}_k^{m^{\&prime;}}\left(s\right){\mathrm{\&Phi;}}_k({\mathrm{\&Omega;}}^{m^{\&prime;}},{\mathrm{\&Omega;}}^m){\mathrm{d\&Omega;}}^{m^{\&prime;}},---\left(11\right)$

边界条件为：

$I_{w,k}^m={\left(\frac{n}{n_0}\right)}^2(1-{\mathrm{\&rho;}}_{0,k})I_{0,k}^m+\frac{f_0{\mathrm{\&rho;}}_k}{\mathrm{\&pi;}}{\&Integral;}_{n_w\&CenterDot;{\mathrm{\&Omega;}}^{m^{\&prime;}}<0}{I}_{w,k}^{m^{\&prime;}}|n_w\&CenterDot;{\mathrm{\&Omega;}}^{m^{\&prime;}}|{\mathrm{d\&Omega;}}^{m^{\&prime;}}+(1-f_0)I_{w,k}^{m^{\&prime;\&prime;}},---\left(12\right)$

采用有限体积法求解半透明介质内部的方向谱带辐射强度

$I_{k,P}^m=(a_{k,Y}^m{I}_{k,Y}^m+a_{k,Z}^m{I}_{k,Z}^m+c_{k,P}^m)/a_{k,P}^m,---\left(13\right)$

其中： $a_{k,P}^m=\underset{j=z,y}{\mathrm{\&Sigma;}}\max (A_j{D}_j^m,0)+{\mathrm{\&kappa;}}_{e,k,P}{V}_P{\mathrm{\&Delta;\&Omega;}}^m,---(13-1)$

$a_{k,J}^m=\max ({-A}_j{D}_j^m,0)(J=Z,Y,j=z,y),---(13-2)$

$C_{k,P}^m={\mathrm{\&kappa;}}_{e,k,P}\&CenterDot;U_{k,P}^m\&CenterDot;V_P\&CenterDot;{\mathrm{\&Delta;\&Omega;}}^m,---(13-3)$

$U_{k,P}^m=\frac{{\mathrm{\&kappa;}}_{a,k,P}}{{\mathrm{\&kappa;}}_{e,k,P}}\frac{{\mathrm{\&sigma;T}}_P^4}{\mathrm{\&pi;}}+\frac{{\mathrm{\&kappa;}}_{s,k,P}}{{\mathrm{\&kappa;}}_{e,k,P}}\&CenterDot;\frac{1}{4\mathrm{\&pi;}}\underset{m^{\&prime;}}{\mathrm{\&Sigma;}}{I}_{k,P}{\mathrm{\&Phi;}}_k({\mathrm{\&Omega;}}^m,{\mathrm{\&Omega;}}^{m^{\&prime;}}){\mathrm{\&Delta;\&Omega;}}^{m^{\&prime;}},---(13-4)$

重复前述过程，直至求得收敛的半透明介质内部的方向谱带辐射强度

即前后两次迭代的相对误差小于预设精度阈值；

将半透明介质内部沿与被测材料表面垂直的方向进行网格划分，均分为多个网格单元，所述网格单元与被测材料表面平行；根据半透明介质内部每个网格单元的中心点温度T_P计算出边界出射辐射热流密度q_w，P：

$q_{w,P}=\underset{k=1}{\overset{M_b}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_{w,k}[B_{k,T_w}{\mathrm{\&sigma;T}}_w^4-\underset{m=1}{\overset{N_{\mathrm{\&Omega;}}/2}{\mathrm{\&Sigma;}}}{I}_{w,k}^m{D}_w^m],---\left(14\right)$

最后根据下述公式计算获得测温设备能够获得的理论辐射能量值Q，完成正算：

$Q={\mathrm{\&tau;q}}_{w,P}{A}_P+(1-\mathrm{\&tau;}+\mathrm{\&rho;\&tau;})\underset{k=1}{\overset{M_n}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_k{\mathrm{\&sigma;T}}_g^4.---\left(15\right)$

本实施方式中，由于半透明介质与被测材料表面不接触，因此属于辐射非平衡问题，只需求解辐射传输方程，而不需要求解能量守恒方程。

具体实施方式四：本实施方式为对实施方式二或三的进一步说明，步骤五中采用智能微粒群优化算法反演被测材料表面的真实温度值的具体方法为：

将步骤二或步骤三中获得的理论辐射能量值Q与步骤四中获得的实际辐射能量值Q_m作差，以该差值的最小二乘平方为目标函数，目标函数OF表示如下：

$\mathrm{OF}=\frac{1}{2}{[Q_m-Q]}^2,---\left(16\right)$

对被测材料表面和半透明介质内部温度的反演，采用以目标函数OF值为适应度值的智能微粒群优化算法进行迭代求解；

首先：假设在n维搜索空间内，由M个微粒组成一个微粒群，其中第i个微粒在n维搜索空间内的解为X_i，当半透明介质与被测材料表面接触时，X_i表示半透明介质内和被测材料表面的边界温度；当半透明介质与被测材料表面不接触时，X_i表示被测材料表面的温度，将X_i代入目标函数OF中，计算出目标函数OF的适应度值，并根据该适应度值衡量X_i的优劣；

采用P_i表示第i个微粒在n维搜索空间内所经历的最优解；同时，每个微粒的飞行速度为V_i，所有微粒经历过的位置中全局最优解为P_g，与P_g相应的全局适应度值为F_g，对于每一个微粒在n维搜索空间内的解，其迭代方程如下：

V_i(t+1)＝V_i+C₁·R₁·[P_i(t)-X_i(t)]+C₂·R₂·[P_g(t)-X_i(t)]，    (17)

X_i(t+1)＝X_i(t)+V_i(t+1)，    (18)

其中，C₁和C₂均为迭代系数，为正常数，C₁用来调节微粒飞向自身最优解方向的步长，C₂用来调节微粒飞向全局最优解方向的步长；R₁和R₂均为在[0，1]范围内变化的随机数；

其次：在对微粒进行搜索时，微粒X_i的值被最大值x_max和最小值x_min限制，

当某一微粒的值大于最大值x_max时，被强制赋值为x_max；

当某一微粒的值小宇最小值x_min时，被强制赋值为x_min；

最后：更新所有微粒X_i的值时，同时更新P_g和P_i的相应值，循环迭代计算，迭代终止的条件是达到最大迭代次数或达到设定的最小适应度值，当适应度值小于某一预设精度值时，对应的微粒群中的全局最优解P_g即为被测材料表面的真实温度值，同时获得半透明介质内部真实温度值。

本实施方式中，目标函数的值即为智能微粒群算法的适应度值，适应度值越小表示反演的参数值越准确，最小适应度值对应的温度值X_i即视为被测材料表面的真实温度值。

在智能微粒群优化算法的过程中，微粒最大位置x_max和最小位置x_min由本发明方法中要反演的参数的物理意义决定。

具体实施方式五：本实施方式为对实施方式四的进一步说明，所述智能微粒群优化算法的参数选择如下：

微粒的个数范围为20～50；C₁＝C₂∈[0，2]。

本实施方式对微粒群中微粒个数的选择在20～50之间，既能保证PSO算法的计算精度又可保证其计算效率。

具体实施方式六：本实施方式为对实施方式五的进一步说明，所述C₁＝C₂＝1。

对于标准PSO算法而言，C₁＝C₂＝1可使得微粒群算法达到最佳收敛速度，有助于快速搜索到温度最优值。

本发明方法中的测温设备可以采用红外热像仪。

本发明方法首先对测温模型进行了理论建模，做出一些合理的简化和假设，然后根据不同的工况选择不同的数学物理模型，再根据物性参数和几何参数，通过正算和反演，最终得到被测材料表面更为准确的温度以及半透明介质中的温度分布。

Claims (3)

1.一种半透明介质环境下非接触测温的校正方法，其特征在于：它包括以下步骤： 

步骤一：判断半透明介质与被测材料表面是否接触，若接触，执行步骤二；否则，执行步骤三； 

步骤二：选择一维耦合换热模型，采用有限体积法进行正向模型的计算，得到测温设备能够获得的理论辐射能量值，然后执行步骤四； 

步骤三：选择一维纯辐射换热模型，采用有限体积法进行正向模型的计算，得到测温设备能够获得的理论辐射能量值，然后执行步骤四； 

步骤四：采用测温设备实际测量被测材料表面，得到测温设备获得的实际辐射能量值； 

步骤五：根据所述理论辐射能量值和实际辐射能量值，采用智能微粒群优化算法反演被测材料表面的真实温度值； 

将步骤二或步骤三中获得的理论辐射能量值Q与步骤四中获得的实际辐射能量值Q_m作差，以该差值的最小二乘平方为目标函数，目标函数OF表示如下： 

对被测材料表面和半透明介质内部温度的反演，采用以目标函数OF值为适应度值的智能微粒群优化算法进行迭代求解； 

首先：假设在n维搜索空间内，由M个微粒组成一个微粒群，其中第i个微粒在n维搜索空间内的解为X_i，当半透明介质与被测材料表面接触时，X_i表示半透明介质内和被测材料表面的边界温度；当半透明介质与被测材料表面不接触时，X_i表示被测材料表面的温度，将X_i代入目标函数OF中，计算出目标函数OF的适应度值，并根据该适应度值衡量X_i的优劣； 

采用P_i表示第i个微粒在n维搜索空间内所经历的最优解；同时，每个微粒的飞行速度为V_i，所有微粒经历过的位置中全局最优解为P_g，与P_g相应的全局适应度值为F_g，对于每一个微粒在n维搜索空间内的解，其迭代方程如下： 

V_i(t+1)＝V_i+C₁·R₁·[P_i(t)-X_i(t)]+C₂·R₂·[P_g(t)-X_i(t)]，     (17) 

X_i(t+1)＝X_i(t)+V_i(t+1)，         (18) 

其中，V_i(t+1)表示t+1时刻下第i个微粒的飞行速度；P_i(t)表示第t时刻下第i个微粒在n维搜索空间内所经历的最优解；X_i(t)表示第t时刻下第i个微粒在n维搜索空间内的解；P_g(t)表示第t时刻下所有微粒经历过的位置中全局最优解；X_i(t+1)表示第t+1时刻下第i 个微粒在n维搜索空间内的解，C₁和C₂均为迭代系数，为正常数，C₁用来调节微粒飞向自身最优解方向的步长，C₂用来调节微粒飞向全局最优解方向的步长；R₁和R₂均为在[0，1]范围内变化的随机数； 

其次：在对微粒进行搜索时，微粒X_i的值被最大值x_max和最小值x_min限制， 

当某一微粒的值大于最大值x_max时，被强制赋值为x_max； 

当某一微粒的值小于最小值x_min时，被强制赋值为x_min； 

最后：更新所有微粒X_i的值时，同时更新P_g和P_i的相应值，循环迭代计算，迭代终止的条件是达到最大迭代次数或达到设定的最小适应度值，当适应度值小于某一预设精度值时，对应的微粒群中的全局最优解P_g即为被测材料表面的真实温度值，同时获得半透明介质内部真实温度值。 

2.根据权利要求1所述的半透明介质环境下非接触测温的校正方法，其特征在于：所述智能微粒群优化算法的参数选择如下： 

微粒的个数范围为20～50；C₁＝C₂∈[0，2]。 

3.根据权利要求2所述的半透明介质环境下非接触测温的校正方法，其特征在于：所述C₁＝C₂＝1。 

Images