红外成像制导中梯度温度气体介质的热辐射噪声获取方法

摘要

本发明公开了一种红外成像制导中梯度温度气体介质的热辐射噪声获取方法，包括以下几个步骤：步骤一：离散介质建立有限元分析环境；步骤二：光线追迹，寻找LOS路径；步骤三：采用离散传递法迭代获取谱带热辐射强度；步骤四：获取热辐射噪声。本发明克服辐射能量对空间方向的依赖性，使对于特定研究目的高温气体辐射问题可数值求解；引入LOS路径为计算梯度温度介质对红外成像形成扰动的热辐射效应的有效积分路径；发展温度梯度场中基于CFD网格的精确光线追迹模型，实现了包含温度函数的光线轨迹方程和辐射轨迹方程；利用辐射离散传递法，得到高温气体热辐射效应获取方法。

说明书

技术领域

本发明涉及红外成像制导中梯度温度气体介质的热辐射噪声获取方法，属于高速飞行器红外末制导技术领域。

背景技术

带有红外成像制导系统的导弹在大气层内高速飞行时，其光学头罩周围存在高温激波层。激波层中的气体密度、温度、组分高度不均匀，并产生强烈的红外辐射，对探测器成像产生热辐射噪声干扰，使得目标图像严重退化，称为气动热辐射效应。

气动光学热辐射效应是制约高速空空导弹应用红外成像制导技术实现高精度命中目标的关键问题，需要对其形成机理和校正方法进行深入的研究。

图1给出了气动热辐射效应的风洞测试试验图像。其中的a1-a3是透过不同温度ZnS窗口材料红外图像(温度T＝336K，T＝480K，T＝1600K)，b1-b3不同吹风时间红外图像(吹风时间t＝4秒，20秒，25秒)。

气动光学热辐射效应的研究，目前业界大多学者从分子热运动的角度研究大气热辐射机理，未结合光学传输，不能定量得到高温激波对红外目标成像造成具体波前畸变；而通过对气体介质求解辐射传递方程来研究此问题的学者，多采用统计方法(蒙特卡洛)来解决辐射积分路径问题，需要大量的样本，耗费巨大的计算量，并假定不同方向上的特征射线对此节点射出的谱带辐射强度贡献的权重相同，不能确定专用于红外成像噪声分析时辐射强度计算的积分路径。

发明内容

本发明的目的是是为了解决上述问题，提出一种红外成像制导中梯度温度气体介质的热辐射噪声获取方法。

红外成像制导中梯度温度气体介质的热辐射噪声获取方法，选定LOS路径为沿CFD网格计算时的辐射积分路径，根据特定研究目标，寻找到热辐射能量传递所依赖的空间方向；把高温气体介质划分为若干等温、等物性的体元，建立离散温度场并计算转化为离散折射率场，利用几何光学光线追迹的方法，寻找光线传播路径，运用动态规划理论证明其为LOS路径；运用离散传递法计算进入和离开每个单元格网控制体的谱带辐射强度，直至迭代得到光线穿过整体非均匀高温气体介质的总谱带辐射强度，包括以下几个步骤：

步骤一：离散介质建立有限元分析环境；

选取尺寸为a×b×c的长方体气体介质，单位为毫米，内温度、折射率分布分别为T(x，y，z)，n(x，y，z)，气体介质与环境的界面是镜反射半透明界面，对于在三维空间任意连续分布的温度T(x，y，z)，把三维气体介质离散为a×b×c个网格，每个网格的尺寸为1mm×1mm×1mm；

将每个网格的内温度设定为一阶函数，密度和组分各向均匀；当辐射能在网格之间的界面传递时，只考虑折射或全反射；

沿展向方向，即平行于Z轴，垂直于X轴的方向，对每个网格取界面，得到a个二维平面，通过介质离散，所瞄准目标辐射的红外光线穿过三维介质的热辐射噪声获取过程被分解为a个光线穿过二维平面介质的热辐射噪声获取过程；

步骤二：光线追迹，寻找LOS路径；

设网格单元大小为d，即1mm×1mm；1-1表示第一行第一列处的节点，1-2表示第一行第二列处的节点，以此类推；S_i表示网格名称，θ_i表示光线在第S_i个网格中的偏移角；X_i表示光线在第S_i个网格中的偏移量；节点1-1处对应的温度为T₁₁，节点i-j处温度为T_ij；1-1节点与2-1节点间的界面以上网格温度由网格上、下左节点的温度值确定T₁₁(z)＝(T₂₁-T₁₁)z+T₁₁，界面下网格内的温度取T₂₁(z)＝(T₃₁-T₂₁)z+T₂₁，z为导弹行进反方向；

设携能光线初始入射角为θ₀，通过节点1-1后辐射传递轨迹切线方向与传递方向间的夹角为θ₁，第k次传递时对应切向角为θ_k，第k个拐变点处对应的坐标为(x，x(z))；其中，k为光传递轨迹z方向改变次数；

令Δx_k表示光线平移量，Δx₁＝线段X₁O，Δx₂＝线段X₂X₁，……Δx_k＝线段X_kX_k-1，用表示光线的总的实际偏移量；

其中x_i(z)为每个网格单元中的光线方程，是z的函数；对于三维(x，y，z)网格中的光线方程，变成沿z方向的一维方程；T(z)表示网格内温度，n(z)表示网格内折射率，θ₀表示进入这个网格的光线初始偏折角，θ表示在这个网格中发生的角度偏折，也是辐射轨迹上点(z_i，x(z_i))处的切向极角，s(z)表示辐射轨迹路径；

一维梯度温度介质，得到光线方程x＝x(z)为：

$x\left(z\right)=$

${\int }_0^z\frac{T\left(0\right)T\left(z\right)\sin {\theta }_0+{\delta }_{0}T\left(z\right)\sin {\theta }_0}{\sqrt{(T^2\left(0\right)+T^2\left(z\right))(1-{\sin }^2{\theta }_0)+2{\delta }_{0}T\left(0\right)T\left(z\right)(T\left(0\right)-T\left(z\right){\sin }^2{\theta }_0)+{\delta }_0^2(T^2\left(0\right)-T^2\left(z\right){\sin }^2{\theta }_0)}}\mathrm{dz}$

其边界条件为x(0)＝0，x′(0)＝tanθ₀；

式中δ₀＝T₀(n₀-1)，T₀和n₀为参考状态空气温度和折射率；T(0)为界面z＝0处介质的温度；

时，方程x(z)恒有解，若以点(0，x(0))处为原点，光线在x方向移动|x(z)|距离，通过判断|x(z)|是否大于1mm，确定下一个格网单元；

时，方程x(z)无解，光线在到达点(z，x(z))之前被全反射，光线不能到达下一网格的界面处；

时，从光线所在网格的初始界面出发的光线到达点(z，x(z))被全反射回此网格初始界面；

当规定网格内温度函数T_i(z)＝(T_i+1-T_i)z+T_i时，得到了光线偏移量x_i+1(z)的解析表达式；

$x_{i+1}\left(z\right)=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(z\right)}-\frac{m_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)z\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(z\right)}\end{array}\right)|-{\zeta }_i)$

$\sin {\theta }_{i+1}=\frac{x_{i+1}\left(z\right)}{\sqrt{1+{x^2}_{i+1}\left(z\right)}}$

式中：m₀＝T₀(n₀-1)，T₀和n₀分别为所研究激波外大气的平均温度和折射率，$m_i={(\frac{m_0}{T_i}+1)}^2\sin {\theta }_i^2;$

f(z)＝(1-m_i)(T_i+1-T_i)²z²+2(T_i+1-T_i)(T_i+m₀)z+(T_i+m₀)²

${\zeta }_i=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(T_i\right)}-\frac{m_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)T_i\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(T_i\right)}\end{array}\right)\left|\right);$

如果∑Δx_k＞l×d，l∈N为计数值，用于标定所追击光线是否超出给定的CFD数据网格范围，超出即表示计算无效；

步骤三：采用离散传递法迭代获取谱带热辐射强度；

(1)确定辐射传递轨迹；

以点(0，x(0))处为原点，则由光线方程式得出点(z，x(z))在辐射轨迹上的坐标为：

$s\left(z\right)={\int }_0^z\frac{(T\left(0\right)T\left(z\right)+{\delta }_{0}T\left(0\right))\mathrm{dz}}{\sqrt{(T^2\left(0\right)+T^2\left(z\right))(1-{\sin }^2{\theta }_0)+2{\delta }_{0}T\left(0\right)T\left(z\right)(T\left(0\right)-T\left(z\right){\sin }^2{\theta }_0)+{\delta }_0^2(T^2\left(0\right)-T^2\left(z\right){\sin }^2{\theta }_0)}}$

由点(z，x(z))在辐射轨迹上的坐标可得，从辐射轨迹上点(z_i，x(z_i))到达另一点(z_i+1，x(z_i+1))的光线行程长度为：

${\delta }_{s(i+1)}\left(z\right)=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(z\right)}-\frac{{m_{i}m}_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)z\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(z\right)}\end{array}\right)|-{\zeta }_i)$

式中

${\zeta }_i=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(T_i\right)}-\frac{m_i{m}_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)T_i\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(T_i\right)}\end{array}\right)\left|\right),$

θ_i为辐射轨迹上点(z_i，x(z_i))处的切向极角；

因此确定得到了辐射特征射线在控制体单元网格中的行程δ_s(i+1)(z)的解析表达式；

(2)运用离散传递法

设特征射线s(z)进入某单元格网(i，j)时的谱带辐射强度为I_i，j，离开该网格的谱带辐射强度为I_i，j+1；

沿射线在控制体(i，j)上积分辐射传递方程直至最终离开此网格，忽略散射，射线穿过控制体时谱带辐射强度的变化为：

$I_{i,j+1}=I_{i,j}\exp (-{\kappa }_k{\delta }_{s\left(j\right)})+\frac{B_{k.T_{i,j}}\sigma (T_{i,j}+{\delta }_0)T_{i,j}^4}{\pi }[1-\exp (-\kappa {\delta }_{s\left(j\right)})]$

其中δ_s(j)：特征射线在控制体(i，j)中的行程，T_i，j：控制体(i，j)的温度；温度下Δλ_k谱带内黑体辐射能占总辐射能的份额；

$B_{k,T_{i,j}}=\underset{{\mathrm{\Delta \lambda }}_k}{\int }{E}_{\mathrm{b\lambda }}\left(T_{i,j}\right)\mathrm{d\lambda }/\underset{0}{\overset{\infty }{\int }}{E}_{\mathrm{b\lambda }}\left(T_{i,j}\right)\mathrm{d\lambda }=\underset{{\mathrm{\Delta \lambda }}_k}{\int }\frac{c_1{\lambda }^{-5}}{\exp {[c_2/\left(\lambda T_{i,j}\right)]}^{-1}}\mathrm{d\lambda }/\underset{0}{\overset{\infty }{\int }}\frac{c_1{\lambda }^{-5}}{\exp {[c_2/\left(\lambda T_{i,j}\right)]}^{-1}}\mathrm{d\lambda }$

其中Δλ_k为高速导弹制导红外成像涉及的3～8μm中远红外波段；

步骤四：获取热辐射噪声；

联立方程

$\left(\begin{array}{c}\mathrm{\Sigma \Delta }{x}_k=\underset{i=1}{\overset{k}{\Sigma }}{x}_i\left(z\right)\\ \left(\begin{array}{cc}x\left(z\right)=& \end{array}\right)\\ {\int }_0^z\frac{T\left(0\right)T\left(z\right)\sin {\theta }_0+{\delta }_{0}T\left(z\right)\sin {\theta }_0}{\sqrt{(T^2\left(0\right)+T^2\left(z\right))(1-{\sin }^2{\theta }_0)+2{\delta }_{0}T\left(0\right)T\left(z\right)(T\left(0\right)-T\left(z\right){\sin }^2{\theta }_0)+{\delta }_0^2(T^2\left(0\right)-T^2\left(z\right){\sin }^2{\theta }_0)}}\mathrm{dz}\\ x_{i+1}\left(z\right)=\frac{1}{T_{i+1}-T_i}(\frac{1}{1{-}_{m}i}\sqrt{f\left(z\right)}-\frac{m_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)z\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(z\right)}\end{array}\right)|-{\zeta }_i)\\ \sin {\theta }_{i+1}=\frac{x_{i+1}\left(z\right)}{\sqrt{1+{x^2}_{i+1}\left(z\right)}}\end{array}\right)$

通过X_(i+1)(z)得到光线穿过介质后因热辐射造成的光线所经历的实际几何路程，据此得到成像畸变分析需要的光程差；

联立方程

$\left(\begin{array}{c}s\left(z\right)={\int }_0^z\frac{(T\left(0\right)T\left(z\right)+{\delta }_{0}T\left(0\right))\mathrm{dz}}{\sqrt{(T^2\left(0\right)+T^2\left(z\right))(1-{\sin }^2{\theta }_0)+2{\delta }_{0}T\left(0\right)T\left(z\right)(T\left(0\right)-T\left(z\right){\sin }^2{\theta }_0)+{\delta }_0^2(T^2\left(0\right)-T^2\left(z\right){\sin }^2{\theta }_0)}}\\ {\delta }_{s(i+1)}\left(z\right)=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(z\right)}-\frac{m_i{m}_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)z\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(z\right)}\end{array}\right)|-{\zeta }_i)\\ I_{i,j+1}=I_{i,j}\exp (-{\kappa }_k{\delta }_{s\left(j\right)})+\frac{B_{k.T_{i,j}}\sigma (T_{i,j}+{\delta }_0)T_{i,j}^4}{\pi }[1-\exp (-{\mathrm{\kappa \delta }}_{s\left(j\right)})]\\ B_{{k,T}_{i,j}}=\underset{\Delta {\lambda }_k}{\int }{E}_{\mathrm{b\lambda }}\left(T_{i,j}\right)\mathrm{d\lambda }/\underset{0}{\overset{\infty }{\int }}{E}_{\mathrm{b\lambda }}\left(T_{i,j}\right)\mathrm{d\lambda }=\underset{\Delta {\lambda }_k}{\int }\frac{c_1{\lambda }^{-5}}{\exp {[c_2/\left(\lambda T_{i,j}\right)]}^{-1}}\mathrm{d\lambda }/\underset{0}{\overset{\infty }{\int }}\frac{c_1{\lambda }^{-5}}{\exp {[c_2/\left(\lambda T_{i,j}\right)]}^{-1}}\mathrm{d\lambda }\end{array}\right)$

通过I_i，j+1得到光线穿过所分析的介质后因热辐射造成的辐射强度，可据此得到光电探测仪成像时的能量畸变；

光线穿过介质后热辐射造成光程差和能量改变组成热辐射噪声，最后获取基于红外成像制导时高速导弹光学窗口外的梯度温度气体介质造成的热辐射噪声。

本发明的优点在于：

(1)克服辐射能量对空间方向的依赖性，使对于特定研究目的高温气体辐射问题可数值求解；

(2)引入LOS路径为计算梯度温度介质对红外成像形成扰动的热辐射效应的有效积分路径；

(3)发展温度梯度场中基于CFD网格的精确光线追迹计算模型，实现了包含温度函数的光线轨迹方程和辐射轨迹方程；

(4)利用辐射离散传递法，得到高温气体热辐射效应获取方法；

(5)由于热辐射效应计算模型的建立，图像复原和识别的约束条件和先验知识的丰富和具体，建立更加符合实际物理空间结构的图像退化模型。

附图说明

图1是背景技术中气动热辐射效应的风洞测试试验图像示意图；

图2是本发明的方法流程图；

图3是三维介质格网划分图；

图4是截取展向的二维平面的方法示意图；

图5是介质离散最终得到的二维平面介质图

图6为不均匀温度介质的光线追迹图；

图7为光线沿z方向传递轨迹示意图；

图8为离散传递的计算模型图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是红外成像制导中梯度温度气体介质的热辐射噪声获取方法，考虑研究目标是最终服务于气动热辐射退化图像的复原，选定LOS(Line Of Sight)路径为沿CFD(Computational Fluid Dynamics)网格计算时的辐射积分路径；建立了满足精度的网格(64×80×64)个，每个网格大小为1mm×1mm×1mm)，依托湍流中光传输机理建模技术，完成辐射计算积分路径的寻找；采用热辐射计算中经典离散传递法，通过气体总体加权和模型，最终可得到热辐射效应的精确数值解。

流程如图1所示，具体包括以下步骤：

步骤一：离散介质建立有限元分析环境；

选取尺寸为64×80×64(单位为毫米)的长方体气体介质，气体介质为高温且物性分布高度不均匀，内温度、折射率分布分别为T(x，y，z)，n(x，y，z)，气体介质与环境的界面是镜反射半透明界面，对于在三维空间任意连续分布的温度T(x，y，z)，把三维气体介质离散为64×80×64个网格，每个网格的尺寸为1mm×1mm×1mm，如图3所示。图3中坐标轴Z代表流向，Y坐标为法向，X坐标为展向，S_i表示网格的标号，下标i表示第i个网格。

将每个网格的内温度设定为一阶函数，密度和组分各向均匀；当辐射能在网格之间的界面传递时，只考虑它的折射或全反射，虽然与热辐射在阶跃界面处的传递不相符，但符合热辐射在梯度温度介质内的传递特点。

沿展向方向，即平行于Z轴，垂直于X轴的方向，在所分析的64×80×64长方形介质的中心，即X轴方向64个格的一半，第32格处取截面。如图4所示，P₀P₁代表穿过截取后平面的光线。得到了二维平面，如图5所示。入射光表示导弹所瞄准目标辐射的红外光线，经过位于导弹外部的、被离散化的梯度温度气体介质，最终用本方法得到的热辐射噪声将反映在光学探测仪所捕获的图像上。

通过介质离散，所瞄准目标辐射的红外光线穿过三维介质的热辐射噪声获取过程被分解为64个光线穿过二维平面介质的热辐射噪声获取过程，本发明解决二维介质平面内热辐射噪声获取方法；同时将辐射传递方程在微元体半球空间的热辐射求解简化为垂直于每个边界网格面元的均匀谱带辐射强度的求解。

步骤二：光线追迹，寻找(LOS(Line Of Sight))路径；

图6中，设网格单元(正方网格)大小为d(1mm×1mm)；1-1表示第一行第一列处的节点，1-2表示第一行第二列处的节点，以此类推；S_i表示网格名称，θ_i表示光线在第S_i个网格中的偏移角；X_i表示光线在第S_i个网格中的偏移量。节点1-1处对应的温度为T₁₁，节点i-j处温度为T_ij。1-1节点与2-1节点间的界面以上网格温度由网格上、下左节点的温度值确定T₁₁(z)＝(T₂₁-T₁₁)z+T₁₁，界面下网格内的温度取T₂₁(z)＝(T₃₁-T₂₁)z+T₂₁，z为导弹行进反方向。

设携能光线初始入射角为θ₀，通过节点1-1后辐射传递轨迹切线方向与传递方向间的夹角为θ₁，第k次(k定义为光(辐射)传递轨迹z方向改变次数)传递时对应切向角为θ_k，第k个拐变点处对应的坐标为(x，x(z))。

令Δx_k表示光线平移量，Δx₁＝线段X₁O，Δx₂＝线段X₂X₁，……Δx_k＝线段X_kX_k-1，用表示光线的总的实际偏移量。

其中x_i(z)为每个网格单元中的光线方程，是z的函数，以x(z)表示。对于三维(x，y，z)网格中的光线方程，变成沿z方向的一维方程。如图7所示，T(z)表示网格内温度，n(z)表示网格内折射率，θ₀表示进入这个网格的光线初始偏折角，θ表示在这个网格中发生的角度偏折，也是辐射轨迹上点(z_i，x(z_i))处的切向极角，s(z)表示辐射轨迹路径。

对于图7所示的一维梯度温度介质，得到光线方程x＝x(z)为：

$x\left(z\right)=$

${\int }_0^z\frac{T\left(0\right)T\left(z\right)\sin {\theta }_0+{\delta }_{0}T\left(z\right)\sin {\theta }_0}{\sqrt{(T^2\left(0\right)+T^2\left(z\right))(1-{\sin }^2{\theta }_0)+2{\delta }_{0}T\left(0\right)T\left(z\right)(T\left(0\right)-T\left(z\right){\sin }^2{\theta }_0)+{\delta }_0^2(T^2\left(0\right)-T^2\left(z\right){\sin }^2{\theta }_0)}}\mathrm{dz}$

其边界条件为x(0)＝0，x′(0)＝tanθ₀

式中δ₀＝T₀(n₀-1)，T₀和n₀为参考状态空气温度和折射率。T(0)为界面1(z＝0)处介质的温度。

时，方程x(z)恒有解，若以点(0，x(0))处为原点，光线在x方向移动|x(z)|距离，通过判断|x(z)|是否大于1mm，可确定下一个格网单元；

时，方程x(z)无解，光线在到达点(z，x(z))之前被全反射，光线不能到达下一网格的界面处。。

时，从光线所在网格的初始界面出发的光线到达点(z，x(z))被全反射回此网格初始界面。

当规定网格内温度函数T_i(z)＝(T_i+1-T_i)z+T_i时，得到了光线偏移量x_i+1(z)的解析表达式。

$x_{i+1}\left(z\right)=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(z\right)}-\frac{m_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)z\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(z\right)}\end{array}\right)|-{\zeta }_i)$

$\sin {\theta }_{i+1}=\frac{x_{i+1}\left(z\right)}{\sqrt{1+{x^2}_{i+1}\left(z\right)}}$

式中：m₀＝T₀(n₀-1)，T₀和n₀分别为所研究激波外大气的平均温度和折射率，$m_i={(\frac{m_0}{T_i}+1)}^2\sin {\theta }_i^2;$

f(z)＝(1-m_i)(T_i+1-T_i)²z²+2(T_i+1-T_i)(T_i+m₀)z+(T_i+m₀)²

${\zeta }_i=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(T_i\right)}-\frac{m_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)T_i\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(T_i\right)}\end{array}\right)\left|\right).$

如果∑Δx_k＞l×d，l∈N为计数值，用于标定所追击光线是否超出给定的CFD(Computational Fluid Dynamics)数据网格范围，超出即表示计算无效。

步骤三：采用离散传递法迭代获取谱带热辐射强度；

(1)确定辐射传递轨迹；

以点(0，x(0))处为原点，则由光线方程式可得出点(z，x(z))在辐射轨迹上的坐标为：

$s\left(z\right)={\int }_0^z\frac{(T\left(0\right)T\left(z\right)+{\delta }_{0}T\left(0\right))\mathrm{dz}}{\sqrt{(T^2\left(0\right)+T^2\left(z\right))(1-{\sin }^2{\theta }_0)+2{\delta }_{0}T\left(0\right)T\left(z\right)(T\left(0\right)-T\left(z\right){\sin }^2{\theta }_0)+{\delta }_0^2(T^2\left(0\right)-T^2\left(z\right){\sin }^2{\theta }_0)}}$

由点(z，x(z))在辐射轨迹上的坐标可得，从辐射轨迹上点(z_i，x(z_i))到达另一点(z_i+1，x(z_i+1))的光线行程长度为：

${\delta }_{s(i+1)}\left(z\right)=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(z\right)}-\frac{{m_{i}m}_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)z\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(z\right)}\end{array}\right)|-{\zeta }_i)$

式中

${\zeta }_i=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(T_i\right)}-\frac{m_i{m}_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)T_i\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(T_i\right)}\end{array}\right)\left|\right),$

θ_i为辐射轨迹上点(z_i，x(z_i))处的切向极角。

因此确定得到了辐射特征射线在控制体单元网格中的行程δ_s(i+1)(z)的解析表达式。

(2)运用离散传递法

离散传递法的基本思想是沿着某根特征射线求解辐射传递方程。从每个边界网格面元的节点向半球空间引出若干条特征射线，每条特征线都由引出点经过半透明介质到达另一边界面，如图8所示，I_i，j和I_i，j+1表示所在网格的谱带辐射强度，A_i，j表示所在网格的辐射出发点，下标i，j表示网格标号。

辐射强度对角度的依赖性是使辐射问题复杂化的关键因素。因为研究目标最终服务于气动热辐射退化图像的复原，选定LOS(Line Of Sight)路径为沿CFD(Computational Fluid Dynamics)网格计算时的辐射积分路径。

设特征射线s(z)进入某单元格网(i，j)时的谱带辐射强度为I_i，j，离开该网格的谱带辐射强度为I_i，j+1，如图8。

沿射线在控制体(i，j)上积分辐射传递方程直至最终离开此网格，忽略散射，射线穿过控制体时谱带辐射强度的变化为：

$I_{i,j+1}=I_{i,j}\exp (-{\kappa }_k{\delta }_{s\left(j\right)})+\frac{B_{k.T_{i,j}}\sigma (T_{i,j}+{\delta }_0)T_{i,j}^4}{\pi }[1-\exp (-\kappa {\delta }_{s\left(j\right)})]$

其中δ_s(j)：特征射线在控制体(i，j)中的行程，T_i，j：控制体(i，j)的温度；温度下Δλ_k谱带内黑体辐射能占总辐射能的份额。

$B_{k,T_{i,j}}=\underset{{\mathrm{\Delta \lambda }}_k}{\int }{E}_{\mathrm{b\lambda }}\left(T_{i,j}\right)\mathrm{d\lambda }/\underset{0}{\overset{\infty }{\int }}{E}_{\mathrm{b\lambda }}\left(T_{i,j}\right)\mathrm{d\lambda }=\underset{{\mathrm{\Delta \lambda }}_k}{\int }\frac{c_1{\lambda }^{-5}}{\exp {[c_2/\left(\lambda T_{i,j}\right)]}^{-1}}\mathrm{d\lambda }/\underset{0}{\overset{\infty }{\int }}\frac{c_1{\lambda }^{-5}}{\exp {[c_2/\left(\lambda T_{i,j}\right)]}^{-1}}\mathrm{d\lambda }$

其中Δλ_k为高速导弹制导红外成像涉及的3～8μm中远红外波段。

步骤四：获取热辐射噪声；

联立方程

$\left(\begin{array}{c}\mathrm{\Sigma \Delta }{x}_k=\underset{i=1}{\overset{k}{\Sigma }}{x}_i\left(z\right)\\ \left(\begin{array}{cc}x\left(z\right)=& \end{array}\right)\\ {\int }_0^z\frac{T\left(0\right)T\left(z\right)\sin {\theta }_0+{\delta }_{0}T\left(z\right)\sin {\theta }_0}{\sqrt{(T^2\left(0\right)+T^2\left(z\right))(1-{\sin }^2{\theta }_0)+2{\delta }_{0}T\left(0\right)T\left(z\right)(T\left(0\right)-T\left(z\right){\sin }^2{\theta }_0)+{\delta }_0^2(T^2\left(0\right)-T^2\left(z\right){\sin }^2{\theta }_0)}}\mathrm{dz}\\ x_{i+1}\left(z\right)=\frac{1}{T_{i+1}-T_i}(\frac{1}{1{-}_{m}i}\sqrt{f\left(z\right)}-\frac{m_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)z\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(z\right)}\end{array}\right)|-{\zeta }_i)\\ \sin {\theta }_{i+1}=\frac{x_{i+1}\left(z\right)}{\sqrt{1+{x^2}_{i+1}\left(z\right)}}\end{array}\right)$

通过X_(i+1)(z)可得到光线穿过介质后因热辐射造成的光线所经历的实际几何路程，可据此得到成像畸变分析需要的光程差。

光程差解释为：1、光程：光在媒质中所经历的实际几何路程与该种媒质的折射率的乘积。

2、光程差：两同频、同振动方向的光源发出的光，经不同媒质，在相遇点处任意时刻的振动相位差唯一地决定于差值n₂r₂-n₁r₁，其中n表示折射率，r表示光线经历的实际几何路程，下标表示第一个和第二个媒质。我们得到了光线经历的实际几何路程，已知气体介质的每个网格的温度，通过温度可以算出每个网格对应的折射率，就可以得到这种光穿过所分析的特殊介质的光程，减去普通大气的折射率与普通大气的实际光线几何路程(即直线)的乘积，就可以得到计算成像畸变所需的光程差。

联立方程

$\left(\begin{array}{c}s\left(z\right)={\int }_0^z\frac{(T\left(0\right)T\left(z\right)+{\delta }_{0}T\left(0\right))\mathrm{dz}}{\sqrt{(T^2\left(0\right)+T^2\left(z\right))(1-{\sin }^2{\theta }_0)+2{\delta }_{0}T\left(0\right)T\left(z\right)(T\left(0\right)-T\left(z\right){\sin }^2{\theta }_0)+{\delta }_0^2(T^2\left(0\right)-T^2\left(z\right){\sin }^2{\theta }_0)}}\\ {\delta }_{s(i+1)}\left(z\right)=\frac{1}{T_{i+1}-T_i}(\frac{1}{1-m_i}\sqrt{f\left(z\right)}-\frac{m_i{m}_0}{\sqrt{{(1-m_i)}^3}}\ln |\left(\begin{array}{c}2(1-m_i)(T_{i+1}-T_i)z\\ +2(T_i+m_0-T_i{m}_i)+2\sqrt{1-m_i}\sqrt{f\left(z\right)}\end{array}\right)|-{\zeta }_i)\\ I_{i,j+1}=I_{i,j}\exp (-{\kappa }_k{\delta }_{s\left(j\right)})+\frac{B_{k.T_{i,j}}\sigma (T_{i,j}+{\delta }_0)T_{i,j}^4}{\pi }[1-\exp (-{\mathrm{\kappa \delta }}_{s\left(j\right)})]\\ B_{{k,T}_{i,j}}=\underset{\Delta {\lambda }_k}{\int }{E}_{\mathrm{b\lambda }}\left(T_{i,j}\right)\mathrm{d\lambda }/\underset{0}{\overset{\infty }{\int }}{E}_{\mathrm{b\lambda }}\left(T_{i,j}\right)\mathrm{d\lambda }=\underset{\Delta {\lambda }_k}{\int }\frac{c_1{\lambda }^{-5}}{\exp {[c_2/\left(\lambda T_{i,j}\right)]}^{-1}}\mathrm{d\lambda }/\underset{0}{\overset{\infty }{\int }}\frac{c_1{\lambda }^{-5}}{\exp {[c_2/\left(\lambda T_{i,j}\right)]}^{-1}}\mathrm{d\lambda }\end{array}\right)$

通过I_i，j+1得到光线穿过所分析的介质后因热辐射造成的辐射强度，可据此得到光电探测仪成像时的能量畸变。

光线穿过介质后热辐射造成光程差和能量改变组成热辐射噪声，最后获取基于红外成像制导时高速导弹光学窗口外的梯度温度气体介质造成的热辐射噪声。