本发明涉及一种目标红外成像可探测性分析方法，具体而言是目标红外成像可探测时段分析方法，属于红外探测技术领域。
背景技术：
自然环境中，所有温度高于绝对零度的物体都在不停地向外界发射红外辐射。由于这种辐射包含有物体本身的特征信息，所以我们利用这种物体自发的辐射，发展出了成熟的红外成像技术。这种技术已经被广泛地应用于军事、工业、科研等各个领域。随着二代焦平面阵列红外成像技术的不断进步，以及其越来越广泛的应用，如何才能准确有效的评价红外成像系统对于特定目标的红外可探测性成为了目前的一个研究热点。
为了解决这个问题，我们需要建立一个标准的目标红外成像探测性能评估系统，根据系统的性能参数，去预测某种探测系统在一个特定环境中识别出目标的概率，避免设计出的探测系统不能满足既定要求，并造成极大的损失。这种红外可探测性分析模型理论上需要结合红外探测器本身的性能参数、当前的大气传输特性以及目标和背景的红外辐射特性进行建模和计算，并运用一定的探测准则来预测目标在特定环境中被探测的概率。并通过模拟出不同时间和大气状况下的大气透过率等参数，从而可以计算出目标最容易被探测的时间和地点。
目前大多数国内的评估方法依然是基于之前的数学模型，进行大量的现场实验取得实验数据，来预测不同作用距离处的探测识别概率，指导成像系统的设计和优化。这种传统的预测方法需要大量的外场试验，会耗费大量的人力物力，而且还必须考虑外场大气环境、探测器的各种典型物理效应以及背景杂波对建模的影响，测试MRTD曲线时，还存在许多人眼辨别的主观因素干扰实验结果的稳定性，使得传统的预测方法花费时间较长，且效果较差。
技术实现要素：
本发明的目的在于克服现有技术的不足，提供一种目标红外成像可探测性分析方法，该方法性能好、操作简洁、并且置信度高。
本发明目的通过如下技术方案予以实现：
提供一种目标红外成像可探测性分析方法，包括如下步骤：
S1.计算目标场景图像中每个像素点对应地理坐标点的温度，获得目标场景的温度场；
S2.计算目标场景图像中，用户指定的目标和背景的温差，以及目标的空间频率f；
S3.根据目标的空间频率f，计算最大可解周期数N；
S4.利用目标和背景的温差计算50％探测概率时的周期数N50；
S5.利用N/N50比值，计算探测概率Pinf。
同时提供一种根据目标选择红外成像系统的方法，包括如下步骤：
S1.计算目标场景图像中每个像素点对应地理坐标点的温度，获得目标场景的温度场；
S2.计算目标场景图像中，用户指定的目标和背景的温差，以及目标的空间频率f；
S3.根据目标的空间频率f，计算最大可解周期数N；
S4.利用目标和背景的温差计算50％探测概率时的周期数N50；
S5.利用N/N50比值，计算探测概率Pinf；
S6.如果红外成像系统探测概率Pinf小于设定阈值，则更换探测性能更高的红外成像系统，返回步骤S2；如果红外探测系统探测概率Pinf大于等于设定阈值，则确定该红外成像系统。
优选的，步骤S3中根据目标的空间频率f，计算最大可解周期数N的具体方法为：
其中Mag为红外系统的放大率，R为目标离光学系统的距离；Cd为目标的特性尺寸，Aw为用户指定目标所在区域面积。
优选的，步骤S4中利用目标和背景的温差计算50％探测概率时的周期数N50的具体方法为：
ΔTRSS为目标和背景的温差，C为背景杂波水平。
优选的，步骤S5中计算探测概率Pinf的具体方法为：
Ratio＝N/N50E＝2.7+0.7*Ratio。
优选的，步骤S1中计算目标场景图像中每个像素点对应地理坐标点的温度，获得目标场景的温度场的具体方法为：
S11.接受用户输入的目标场景的材质图像，根据目标场景的材质图像中的材质编号，从材质库中提取出每种材质的发射率ε和反射率α；
S12.计算目标场景中每个像素点的辐射能；
S13.建立目标场景的热平衡方程，求解每个像素点对应地理坐标点的材质表面温度ti,j，进而获得目标场景的温度场。
优选的，步骤S12中计算目标场景中每个像素点的辐射能的具体方法为：
1)对每个像素点计算接收到的太阳辐射能Qsun；
2)对每个像素点计算接收到的大气长波辐射能Qsky；
3)对每个像素点计算材质表面的辐射能Qrad；
4)对每个像素点计算显性热交换能QCV；
5)对每个像素点计算潜性热交换能QLE。
优选的，步骤S13中建立目标场景的热平衡方程为：
Qsun+Qsky＝Qrad+Qcv+Qle。
优选的，计算目标场景图像中，用户指定的目标和背景的温差的具体方法为：
S21.获取用户指定材质图像中的目标，其余为背景区域；
S22.计算目标所在区域和背景区域的温差ΔTRSS：
式中：ti,j代表目标区域像素点(i,j)处的材质表面温度，POT为目标像素点总个数，μbkg为环境背景温度，通过对背景区域各像素点温度取均值获得。
优选的，计算目标的空间频率f的具体方法为：
(1)计算红外成像系统的传递函数；
(2)计算显示器传递函数；
(3)计算人眼调制传递函数；
(4)将红外成像系统的传递函数、显示器传递函数以及人眼调制传递函数相乘，得到红外成像系统水平和竖直方向的传递函数MTFxs(f)，MTFys(f)；
(5)计算红外成像系统的最小可探测温差函数MRTD(f)，
其中Δfn红外成像系统的为噪声等效带宽，SNRDT为观察者能分辨线条的阈值视觉信噪比，NETD为红外系统的噪声等效温差，α为水平瞬时视场，β为垂直瞬时视场，τd为扫描驻留时间，fp为帧频；MTFs(f)分别取水平传递函数MTFxs(f)，竖直方向的MTFys(f)，计算或的水平和竖直方向的MRTD(f)；
(6)令水平和竖直方向的MRTD(f)分别等于目标背景温差ΔTRSS，求解出目标x、y方向上对应的空间频率fx、fy，则目标的空间频率为
优选的，计算红外成像系统的传递函数的方法为：
红外光经过光学系统进入探测器，转换为电信号后，进入信号处理电路进行处理；计算光学系统的传递函数Hopt(fx,fy)；计算探测器水平和竖直方向传递函数；计算信号处理电路传递函数MTFe；分别计算获得红外成像系统水平方向的传递函数和竖直方向的传递函数。
优选的，计算信号处理电路传递函数MTFe的具体方法为：计算低通滤波电路的传递函数MTFe1；计算高通滤波电路的传递函数MTFe2；计算高频提举电路的传递函数MTFe3；计算信号处理电路的传递函数为：
MTFe＝MTFe1MTFe2MTFe3。
优选的，显示器传递函数计算方法如下：
对于阴极射线管显示器传递函数MTFm计算方法如下：
MTFm＝exp(-2π2δ2f2)
式中δ为显示器光点分布的标准偏差；对于LED显示器水平和竖直方向传递函数MTFmx、MTFmy计算方法如下：
MTFmx＝sinc(πxfx)
MTFmy＝sinc(πyfy)
式中x、y分别为LED在归一化空间内的张角。
优选的，人眼调制传递函数计算方法如下：
式中M为红外成像系统角放大率；K为与显示屏亮度L有关的参量，
K＝1.272081-0.300182lgL+0.04261(lgL)2+0.00197(lgL)3。
本发明与现有技术相比具有如下优点：
(1)本发明通过红外成像性能评估模型将红外成像系统的硬件参数、大气环境和目标背景特性三部分结合起来，利用系统的综合性能参数数学模型，将影响红外探测系统性能的各因素联系起来，结合红外目标背景的特性和大气衰减，并使用特定探测准则来预测红外探测系统的现场性能和目标的可探测性，本发明对目标红外成像可探测性分析的性能好，操作简洁，并且置信度高。
(2)本发明通过数字计算的方式，确定目标可探测的概率，与外场试验相比，更加客观、效率更高、实现成本低。
(3)本发明的目标红外探测性分析方法适用于各类红外系统探测性能的评估以及目标可探测概率的评估，通用性好，可靠性高。
(4)本发明的目标红外可探测性分析方法能够准确评估红外系统探测性能，对于红外系统探测性能改进提供了理论指导依据。对于特定目标，通过对不同距离条件探测概率的计算，确定是否需要更换红外系统，选择适当的红外系统进行探测，保证目标探测概率和探测精度。
附图说明
图1为本发明的流程框图；
图2为本发明的目标提取图；
图3为本发明的分时段可探测性折线图。
具体实施方式
以下结合附图和具体实施例对本发明作进一步的详细描述，但该实施例不应理解为对本发明的限制。
目标红外成像可探测时段分析方法，结合图1，其步骤如下：
S1.计算目标场景图像中每个像素点对应地理坐标点的温度，获得目标场景的温度场：
S11.接受用户输入的目标场景的材质图像和保障数据，根据目标场景的材质图像中的材质编号，从材质库中提取出每种材质的发射率ε和反射率α；；
S12.计算目标场景中每个像素点的辐射能，具体步骤如下：
1)对每个像素点计算接收到的太阳辐射能太阳直射能可以按如下公式计算：
Qsund＝rIscPm cosi
式中，r为日地修正系数；Isc为太阳常数；P为目标场景所在区域的大气透明度，通过气象资料得到，非常好的晴天P取0.85，较好的晴天P取0.8，中等晴天时P取0.65，比较差的晴天P取0.532。m为目标场景所在区域的大气质量，i为该材质对应地物的太阳光入射角，根据保障数据中该像素点的法向量确定。
太阳散射能具体的计算公式如下：
式中，Isc为太阳常数，h为太阳高度角，β为该像素点对应地物的倾斜角，根据保障数据中该像素点的法向量确定。P为目标场景所在区域的大气透明度。
该像素点接收到的太阳辐射可以表示为：
Qsun＝(1-α)(Qsund+Qsuns)
式中，α为该像素点的材质的反射率，Qsund为太阳直射能量，Qsuns为太阳散射辐射能量。
2)对每个像素点计算接收到的大气长波辐射能
大气长波辐射一般和一些气象条件有关，天空无云条件下大气长波辐射为：
式中，a，b为经验常数，a取0.61，b取0.05；ε为材质发射率；σ为Stefan-Boltzmann常数，σ＝5.67×10-8W·m-2·K-4；eα为目标场景所在区域的近地面水汽压，单位是kPa，eα的计算公式如下：
式中，rh为目标场景所在区域的大气的相对湿度，Tα是目标场景所在区域的大气温度。
3)对每个像素点计算自身辐射能
任何温度高于绝对零度的物体都会不断地向空间辐射电磁波，将热能转化为辐射能。工程上可以将物体当作灰体处理，根据Stefan-Boltzmann定律，灰体的辐射值计算公式为：
式中，ε表示热源表面材质的发射率，σ为波尔兹曼常数，，ti,j为该像素点对应材质表面的温度，即像素点对应地理坐标点材质表面温度，为待求量。Qrad则为该像素点材质表面红外辐射量。
4)对每个像素点计算显性热交换能
当物体与周围的环境存在温度差时，就会有热量的交换，方向是从高温物体到低温物体，且温差越大，热量交换也就越大。红外成像仿真中主要涉及物体与周围大气间的热交换，称之为显性热交换。可以表示为：
QCV＝h(ti，j-Tα)
式中，h是与风速有关的参数，假定风速为v，单位是米每秒，接近物体表面的实效风速为V。对于迎风面，V等于0.25v，对于背风面，V等于0.3+0.05v。物体表面的热交换系统为h＝3.5+5.6V。
5)潜性热交换
对于一些含有水分的材质，由于水的蒸发作用，将会从材质中带走一部分的热量。这样由于水分散失引起的热交换称为潜性热交换，可用下式进行计算：
QLE＝ρaLCDvWs(qa-qc)(7-8)
式中，ρa为目标场景所在区域的空气密度；L为水的汽化潜热；CD为风的拖曳系数，取0.002；v为目标场景所在区域的风速；Ws为该像素点对应材质表面含水量；qa为目标场景所在区域空气中的比湿；qc为目标场景所在区域的环境温度下的空气的饱和比湿。
S13.将辐射量带入到热平衡方程，通过牛顿迭代法计算出目标场景的每种个像素点对应地理坐标点的材质表面温度ti,j，热平衡方程如下
Qsun+Qsky＝Qrad+Qcv+Qle
根据每个像素点的温度，获得目标场景的温度场。
S2.计算空间频率：
S21.获取用户指定材质图像中的目标，其余为背景区域；
S22.计算目标区域和背景区域的温差，本模型采用ΔTRSS统计温差，计算方法如下：
式中：ti,j代表目标区域像素点(i,j)处的材质表面温度，POT为目标像素点总个数，μbkg为环境背景温度，对背景区域各像素点温度取均值获得。这种方法较为全面地考虑到了目标温度的变化，比较准确地反映了真实的目标背景温差。
S23.探测器参数仿真
红外成像系统成像，红外光经过光学系统进入探测器，转换为电信号后，进入信号处理电路进行处理，经过显示器显示。红外成像系统的空间传递特性主要包括：光学系统的传递函数，探测器的传递函数，信号处理电路的传递函数，显示器的传递函数。这些物理效应可看作是线性传递效应利用各组成模块的调制传递函数(MTF)进行模拟，故总的系统传递函数可由各个模块的传递函数相乘得到。
1)光学系统的传递函数Hopt(fx,fy)计算：
光学系统的MTF取决于波长、焦距和孔径的形状。在成像系统中，光学系统的衍射限、几何像差也对其光学系统传递函数有很大的影响。通常衍射限对应的传递函数在低频下降较快，在高频拖尾较长；而几何像差所对应的传递函数在低频下降较慢，在高频下降较快且拖尾较短。
由线性理论可知，几何像差和衍射限像差无关，故光学系统总传递函数Hopt(fx,fy)为几何像差和衍射限像差函数之积：
Hopt(fx,fy)≈Hdiff(fx,fy)Haber(fx,fy)
衍射限制的光学系统的传递函数取决于波长和孔径的形状，fx,fy分别为目标水平、垂直方向的空间频率。
对于圆形孔径fx＝fy＝f：
fc＝D0/λ0为非相干光学系统的空间截止频率(cycle/mrad)；D0是光学系统的有效孔径(mm)；λ0为非相干光的中心波长(μm)，可取平均工作波长(λ1+λ2)/2；[λ1，λ2]为工作波长范围；f为空间频率。
对于方形孔径fx＝fy＝f:
fc＝D0/λ0为空间截止频率(cycle/mrad)；D0是光学系统的有效孔径(mm)。
非衍射限光学系统中，由像差引起的弥散斑的能量分布为高斯型，具有圆对称形式，其标准偏差为σc(mm)，其传递函数的形式为
Haber(f)＝exp(-2π2σ2f2)
σ是用角度度量的光学系统的标准方差，σ＝σc/f，σ单位为cycle/mrad，它与弥散圆内所占能量的百分比有直接的关系。假定所有的能量都集中在弥散斑内，那么标准方差σc就等于弥散圆直径的1/4。这种假设能很好地近似弥散斑的形状，而且易于控制。其中弥散斑的直径可由下式确定：
2)探测器水平和竖直方向传递函数计算：
探测器的类型主要有单元、多元(串联、并联、串并联)以及sprite探测器等多种类型，不同类型的探测器传递函数也有区别。其中，单元、多元探测器水平方向的空间传递函数为
式中：α为矩形探测器水平方向的空间张角。
对于sprite探测器，其水平方向的传递函数受两个因素影响，其一是读出区长度产生的影响，读出区影响造成的传递函数为MTFl：
MTFl＝sinc(πlfx)
式中：l为读出区长度对应的角距离(mrad)。
其二，是扩散和扫描速度失调造成的影响。在sprite探测器的扫积区，扩散和扫描速度失调造成的传递函数MTFdv为
MTFdv＝{(DτK2+1)2+[(μaE-vs)τK]2}-1/2
式中：D为载流子扩散系数(cm2·s-1)，为探测器的性能参数；τ为载流子寿命(s)，为探测器的性能参数；K＝2πfx/f′，f′为探测器的焦距(cm)；μa为双极迁移率(cm2·V-1·s-1)；E为偏置电场(V·cm-1)，为探测器的性能参数；vs为扫描速度(cm·s-1)，为探测器的性能参数。
在实际调试过程中，我们可以通过微调每一个条带的工作电压，同时观察其输出结果，从而确定最佳工作电压。这种情况下，可以忽略速度失调项的影响，只需要计算扩散造成的影响因素，此时，传递函数可简化为：
由上可得，sprite探测器的水平方向传递函数MTFdx为上述两个函数之积，即
MTFdx＝MTFlMTFdv
同时，探测器也可以等效为一个RC一阶低通滤波器，其传递函数为
式中：f0为空间频率域的3dB频率(c/mrad)。
通常作为基本参量给出的是时间频率域的3dB频率ft0(Hz)，可由时间频率域的3dB频率换算。有时ft0未给出，但给出载流子寿命τ，则可按下面的方法换算ft0：
ft0＝1/2πτ
MTFdt作为时间滤波，只作用于水平方向。
因此探测器的水平方向传递函数还应当与MTFdt相乘，获得水平方向最终的传递函数MTFdx·MTFdt。
在垂直方向上，探测器是通过离散的扫描线来完成空间采样成像的，因此具有相似的性质。严格讲，垂直方向不满足线性和空间不变性，但当空间频率小于扫描抽样的Nyquist频率时，仍可用传递函数来描述。
垂直方向传递函数形式如下：
式中：fy为垂直方向空间频率(c/mrad)；fNy为Nyquist频率。
3)信号处理电路传递函数MTFe的计算：
信号处理电路包括依次串联的低通滤波电路、高通滤波电路和高频提取电路。
低通滤波电路的传递函数为MTFe1：
低通滤波器在时间频率域的传递函数为
MTF＝[1+(ft/ft0)2]-1/2
式中：ft为目标的空间频率(Hz)；ft0为低通滤波器的3dB频率(Hz)，ft0＝1/2πRC，R、C分别为低通滤波器的电阻和电容。
在热成像系统中，电子线路接收信号的时间频率ft(Hz)和目标的空间频率f(c/mrad)是相关的，转换公式为
ω为红外系统的扫描角速度，τd为红外系统的扫描驻留时间，α为水平瞬时视场。将时间频率转换为空间频率，则可得空间频率域的传递函数为
MTFe1＝[1+(f/f0)2]-1/2
高通滤波电路的MTFe2：
在电路中完成类似微分处理的环节是高通R1C1滤波器，特征频率(即高通滤波器上升3dB频率)ft0为
ft0＝1/2πR1C1
用与低通滤波器类似的方法得：
式中：f0为高通滤波器在空间频率域的3dB频率。
高频提举电路的MTFe3：
高频提举电路是一种MTF增强电路，其MTF值可大于1，最大提举频率点为fmax，提举幅度为K(≥1)，传递函数可表示为
计算信号处理电路总的传递函数为MTFe＝MTFe1MTFe2MTFe3。
4)显示器传递函数的计算
在热成像系统中采用的显示器是阴极射线管(CRT)或者发光二极管(LED)阵列。
对于阴极射线管的MTFm：
通常认为CRT上光点亮度分布是高斯分布，所以传递函数为
MTFm＝exp(-2π2δ2f2)
式中：δ为显示器光点分布的标准偏差(mrad)，如果显示器x和y方向的标准偏差不同，则分别代入x和y方向的σx和δy计算x和y方向MTFm。
对于LED的MTFm：
通常假定LED的发光分布是矩形盒形函数，因此，LED的水平和竖直方向传递函数分别为
MTFmx＝sinc(πxfx)
MTFmy＝sinc(πyfy)
式中：x、y分别为LED在归一化空间内的张角(mrad)。
5)人眼调制传递函数
热成像系统探测的红外辐射图像需要在显示器上输出，最后由人眼观察并由人脑做出相应判断和决策，故热成像系统、人眼、大脑是一个有机结合体，在性能模型中必须考虑人眼的传递特性。
人眼传递函数为：
式中：M为红外系统角放大率，为红外系统的系统参数：K为与显示屏亮度L有关的参量，当L用cd/m2表示时，可表示为
K＝1.272081-0.300182lgL+0.04261(lgL)2+0.00197(lgL)3
S24.将S23步骤计算出的各部分MTF函数相乘，得到系统总的调至传递函数MTFs(f)，分别包括水平和竖直方向的传递函数MTFxs(f)，MTFys(f)。
根据MTFs(f)计算MRTD(f)，获得红外系统的MRTD曲线，计算公式为：
式中：MTFs(f)为系统的调制传递函数对于水平和竖直方向分别替换为MTFxs(f)，MTFys(f)，由各个模块的MTF曲线相乘可得；Δfn红外系统的为噪声等效带宽，SNRDT为观察者能分辨线条的阈值视觉信噪比，一般取2.25，f为目标空间频率，NETD为红外系统的噪声等效温差，α为水平瞬时视场(mrad)，β为垂直瞬时视场(mrad)，τd为扫描驻留时间(s)，fp为帧频(s-1)。
分别求解出x、y方向的MRTD曲线后，x和y方向的MRTD(f)分别等于目标背景温差ΔTRSS，求解出目标x、y方向上对应的空间频率fx、fy，则目标最终的空间频率为
S3.计算最大可解周期数N：
可分解的周期数N代表了目标被探测到的概率大小。一个很高的N值代表的很高的探测概率，定义式如下：
N–目标的可分解周期数(cycles/mrad)
f–空间频率(cycles/mrad)
Mag–红外系统的放大率
Cd–目标的特性尺寸(m)
R–目标离光学系统的距离(km)
特征尺寸为用户指定材质图像中的目标区域面积Aw(m2)开根号计算得出：
S4.计算N50(可探测概率为50％的可解周期数)：
计算N50的公式如下：
N50-50％探测概率时的周期数
ΔTRSS为目标/背景温差
C–背景杂波水平，根据经验人为指定，1.0低，1.5略低，2.0一般，2.7高。
S5.计算探测概率：
接下来需要计算出N/N50比值，以便下一步求解：
Ratio＝N/N50
Ratio–N与N50的比值
N–目标求解周期数
N50-50％探测概率时的周期数
Pinf是目标探测概率，通过一个经验主义的目标概率转换函数计算：
Pinf–目标探测概率
E-2.7+0.7*Ratio
Ratio–N与N50的比值。
可以人为设定阈值，如果Pinf大于等于阈值，表明该红外目标可探测，如果低于阈值表明该红外目标不可探测。如果低于阈值，可更换探测性能更高的红外成像系统进行探测，直至红外探测系统的Pinf大于等于阈值。
如图2所示为本发明的目标提取图；图3所示为本发明的分时段可探测性折线图，本发明通过用户输入的材质图像，利用热平衡方程计算出目标温度场，从温度场中提取出目标，并对探测器性能和目标特性进行建模，计算目标的可探测性。本发明效率高、误差小且执行度高。
以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。
本说明书中未作详细描述的内容，属于本专业技术人员公知的现有技术。