0.   引言

作用距离是红外成像探测系统最重要的评估指标之一，它与目标的辐射强度、背景的辐射亮度、大气透过率、探测器的性能等多种因素有关。早期的作用距离评估方程主要有R. D. 小哈德逊作用距离方程^[1]、基于NETD（Noise Equivalent Temperature Difference，噪声等效温差）的作用距离方程^[2-4]，基于MDTD（Minimum Detectable Temperature Difference，最小可探测温差）的作用距离方程^[2]、基于对比度的作用距离方程等^[4]。当前的作用距离评估模型则同时考虑目标、探测设备和大气环境的相互作用，对于点目标探测来说，把以下3个要求作为点目标能否被探测到的基本条件^[5-6]：

1）目标在探测器上所成的像元必须具有一定大小，一般要求在单方向上的最小像元个数n在1~3之间；

2）目标在靶面上的照度产生的信噪比SNR（signal-to-noise ratio）要满足一定要求；

3）目标和背景在探测器上的照度必须具有一定的对比度，该对比度满足信号检测和处理所需要的最低对比度。

根据以上要求进行作用距离评估时，一般是先根据预想的大气环境条件和观测条件，设定目标到探测设备的距离R，利用LOWTRAN或MODTRAN软件计算平均大气透过率，再根据评估模型计算以上3个条件，如果满足，说明还不到最远作用距离，增大R进一步计算，直到其中有一个条件不满足为止，则此时得到的R值就是最远作用距离。如果初始设定的R值较大，不满足3个条件，则减小R直到3个条件都满足为止。以上方法可以手工计算，也可以通过编程方法自动计算，后者的效率大大提高。文献[7-8]利用编程方法调用MODTRAN软件逐步计算，得到作用距离。但在此方法中，采用通常的固定步长方法，步长太大，则最终的计算精度不够；步长太小，则计算速度较慢。下面采用的变步长方法有效地解决了计算精度和速度之间的矛盾，大大提高了计算效率。

1.   作用距离评估方程

1.1   目标像斑大小计算

不考虑弥散时，目标在焦平面上的像斑大小d₀可用下式表示：

式中：l是目标线度；f是光学系统焦距；R是目标到探测设备的距离。实际的目标成像大小与多种因素有关，主要包括系统对目标的夹角、仪器的振动、跟踪平稳情况、光学系统的弥散、大气湍流的弥散、探测器的分辨率等。

1）无弥散时目标的张角σ₀

2）大气抖动在积分时间内引起的角弥散σ₁

影响大气抖动造成的像点弥散的主要因素有：大气抖动频率和均方根角值、系统帧积分时间等。总的来说，随波长的增大，弥散斑越小。一般情况下，较好天气下的可见光波段的弥散取值为3"或2"。

3）探测器空间分辨率限制引起的角弥散σ₂

式中：α为单元探测器的尺度。

4）光学系统衍射极限引起的角弥散σ₃

式中：D为光学系统口径。

5）积分时间内跟踪系统抖动造成的角弥散σ₄

式中：ω为积分时间内目标在像面上的最大位移。

6）系统总的像斑大小d

考虑弥散情况下，系统总的像斑直径大小为：

7）像斑所占像元个数n

在实际应用中，主要考虑目标短边所成的像元个数。

1.2   信噪比计算

信噪比反映目标在探测器上的信号和探测器本身的噪声的相对强度，信噪比SNR可表示为^[1-2]：

式中：D^*是红外探测系统的探测率；E₀是点目标在探测器上的辐照度；A_d是单个像素的面积；Δf是系统的带宽。一般情况下，D^*、A_d和Δf作为系统参数给出，E₀则需要经计算给出，它的计算公式为：

式中：I₀是目标辐射强度；A_op是光学系统接收面积；τ_op是光学系统透过率；τ是大气平均透过率；A_m是目标在探测器上成像的面积。在这些参数中，A_op和τ_op是系统参数，τ可以由LOWTRAN或MODTRAN计算给出，A_m可以通过计算目标长边和短边在探测器上所成像斑大小的乘积得到，I₀可由下式计算得到：

式中：ε是目标的发射系数；A₀是目标的有效辐射面积；M₀是目标在工作波长λ₁~λ₂内的黑体辐射出射度，它可通过对普朗克公式积分得到：

式中：h是普朗克常数；c是光速；k是玻尔兹曼常数；T是目标的温度。对于气动加热的物体，表面温度可由下式计算^[1]：

式中：T₀是周围大气的温度；r是恢复系数，对于层流，r=0.82；对于紊流，r=0.87。γ取1.4，M是目标运动的马赫数。

1.3   调制对比度计算

调制对比度反映目标和天空背景在探测器上辐照度的相对比值，它可表示为：

式中：E_bg是天空背景在探测器上的辐照度，其计算公式为：

式中：L_bg是天空背景辐亮度，它可由LOWTRAN或MODTRAN计算给出，但两者给出的值是特定地区天空背景辐亮度的年平均值，无法反映L_bg随太阳高度角以及复杂多变的大气条件等的变化情况^[9-10]，所以有条件的话，可以用天光光学辐射特性测量设备进行实际测量。

2.   算法实现

在程序设计上，采用Matlab或VC进行界面设计，采用Matlab的优点是可以直接调用其积分函数integral，对普朗克公式积分得到目标的辐射出射度。在界面上输入目标、探测设备的参数，在LOWTRAN或MODTRAN的接口界面中输入大气模式等参数。LOWTRAN7的光谱分辨率为20 cm^-1，MODTRAN的光谱分辨率提高到1 cm^-1。下面以Matlab调用MODTRAN为例说明程序设计方法。

2.1   Matlab与MODTRAN的接口

主程序与MODTRAN软件的数据交换是通过MODTRAN软件的tape5和tape6文件实现的。在MODTRAN软件的输入界面上选择相应的选项后，MODTRAN会生成一个tape5文件，其中包含了用户输入的各种选项的格式化数据。在执行MODTRAN软件时，MODTRAN调用tape5中的数据进行计算，得到大气透过率和天空辐亮度等数据，保存在tape6文件中。

主程序与MODTRAN软件的接口之一是读取tape6中的大气平均透过率和天空背景辐亮度数据，为此，需要在MODTRAN软件的运行模式（Execution mode）中选择“Thermal Radiance”选项，这样，MODTRAN会同时计算大气透过率和天空背景辐亮度。这两个数据在tape6文件的最后部分，如图 1所示，需要想办法从文本文件中把它们提取出来。

图  1  Tape6文件的部分数据结构

Figure  1.  Part of the data structure of file Tape6

下载: 全尺寸图片 幻灯片

程序设计方法是把tape6中从文件结束前的1000个字符读入一个字符串数组中，再利用字符串比较命令strcmp来逐一比较字符串数组中的元素，如果连续找到3个字符串“AVERAGE”、“TRANSMITTANCE”和“=”，则其后的一个字符串就是大气平均透过率的值。同样，如果连续找到4个字符串“INTEGRATED”、“TOTAL”、“RADIANCE”和“=”，则其后的一个字符串就是天空背景辐亮度的值。

在读取tape6文件时要注意的是，一定要等到tape6中的数据生成后再读取。由于MODTRAN在运行时，是先删除旧的tape6文件，再生成新的tape6文件，所以程序有可能在MODTRAN运行时读到空字符串。为此，程序设计了两重防读空措施，一是在读取tape6时，先判断tape6文件是否存在，如果不存在，则一直等待。二是读取数据后，先判断字符串数组是否为空，如果为空则继续等待。程序代码如下：

while(~exist('tape6', 'file'))

end

fp=fopen ('tape6', 'r');

fseek (fp, -1000, 'eof');

while(1)

str=textscan(fp, '%s');

if (~isempty(str))

break;

end

end

fclose(fp);

主程序与MODTRAN软件的另一个接口是改写tape5文件中的R值，即目标与探测设备的距离。改写的方法是把tape5文件读入到一个字符串数组中，删去原tape5文件，把对应R值的元素值改写后再按格式写入到新的tape5文件中。另一个更简单的方法是先用试探法找到R值在文件中的位置，然后直接在此处写入新的R值即可。例如，假设R值在距文件结尾处前112个字符处，则可用以下代码实现：

fp=fopen('tape5', 'r+');

fseek(fp, -112, 'eof');

fprintf(fp, '%7.3f', R/1000);

fclose(fp);

由于tape5中的距离以km为单位，所以R值在写入前要先转换一下。

2.2   变步长的循环算法设计

程序运行时，先设定一个较小的R值，计算3个约束条件，如果满足，则每次加一个较大的步长。例如初次步长设为100 km，等到其中一个约束条件不满足时，则反向减小步长，第二次的步长可设为25 km。等到再次满足条件时，完成了程序的一次循环，此时作用距离的精度是25 km。然后再执行一次循环，正向步长设为5 km，反向步长设为1 km，作用距离的精度即为1 km。如果想进一步提高精度，多执行几次循环即可，每次循环的步长都逐步减小，直到满足需要的精度要求为止。程序设计流程图如图 2所示。

图  2  变步长的一次循环计算程序框图

Figure  2.  Program diagram for variable step length loop calculation

下载: 全尺寸图片 幻灯片

要注意的是，Matlab程序在用while语句进行循环运行时，GUI界面中编辑框的数据不更新显示。而为了监视程序中间的运行结果，又需要观察每次循环计算时像元个数、信噪比和调制对比度的变化情况，因此需要在输出显示代码后，加一条延时程序，延时时长为1 ms即可，以不影响程序运行速度，具体代码如下：

set(handles.edit14, 'string', R/1000);

set(handles.edit15, 'string', n);

set(handles.edit16, 'string', SNR);

set(handles.edit17, 'string', Cm);

pause(0.001);

其中的pause语句起到延时作用，这样，GUI编辑框中的数据即可正常更新显示。此外，Matlab在调用MODTRAN程序时，为了不对GUI的显示界面造成干扰，需要MODTRAN程序以后台方式运行，可执行以下代码实现：

system('mod4.exe', '-echo');

3.   运行结果

软件在Matlab中的设计界面如图 3所示。在主程序界面中输入目标参数和设备参数，点击“作用距离计算”按纽，程序就自动运行，并在输出框中显示每次循环计算的中间及最终结果。

图  3  主程序界面

Figure  3.  Interface of main program

下载: 全尺寸图片 幻灯片

变步长循环计算的效率是很高的。假定初始R值为50 km，最终的作用距离值为300 km。采用固定步长计算时，为了达到1 km的精度，需要设定步长为1 km，这样程序需要循环计算250次。采用前述的变步长方法，第一次循环需要3+3=6次，此时R=275 km。第二次循环需要5+5=10次。这样总共只需要16次就可达到所需的精度。程序耗时主要在MODTRAN程序的运行上。假定程序循环计算一次耗时0.5 s，固定步长计算需要耗时125 s，变步长计算只需要8 s。

4.   结论

本文在Matlab环境下结合MODTRAN软件，根据红外成像探测系统评估模型实现了作用距离的自动计算，为红外系统性能评估提供了便利手段。针对固定步长计算耗时较长的不足，采用了变步长的循环计算方法，大大提高了计算效率。