基于CUDA的实时红外辐射传输及成像仿真设计

蔡方伟，卞春江，冯水春，曹景，毛博年

（中国科学院国家空间科学中心北京100190）

近几年来，科学研究对不同天候、温度、场景、干 扰情况下的红外图像的需求越来越大。但多数情况下，外场实拍不仅要投入大量的人力和物力而且还要耗费相当长的时间[1]；同时，该行为往往与军事目的有关，因而很可能根本无法进行，且无法穷尽所有天候条件。在该情况下，通过计算机仿真手段获取红外图像就显得方便而快捷，并可应用于军事、农业、地质、气象等多个领域[2]。

本文采用计算机数学仿真技术，利用OpenGL强大的图形算法，内容包括目标、背景的红外辐射传输及成像仿真以及二者的融合三部分。针对目标与背景在成像上的区别，我们将目标视为点辐射源的成像，而将背景分解为小面辐射源成像。对于目标，先利用辐射量计算公式直接计算成像辐射数据，再进行相应的成像显示。对于背景，在OpenGL中加载背景三维场景模型[3]，建立面元的红外成像链路，并进行场景模型绘制及渲染，通过合适的几何投影变换得到其二维图像。最终融合计算中，在背景成像上叠加目标成像影响，完成本文的设计任务。

需要注意的是，目前，不仅是在红外成像仿真方面，在计算机软件开发上的主流思想仍是CPU上的串行思想。随着计算机硬件，尤其是各厂商制造的GPU水平的提高，图形图像处理能力大幅进步，处理的面元与像元数据量也因此越来越大，但其计算、处理方法却又是类似的。如果采用按部就班的串行思想开发，不仅无法实现快速处理，也是对硬件资源的极大浪费。而GPU上的并行思想显然给这个问题的解决提供了新的思路。而CUDA便是目前实现GPU上并行处理最成熟、实用化程度最高的技术[4]。

1 目标辐射成像仿真及辐照度计算

上文提过，将目标视为“点辐射源”。其原因在于实际中，目标尺寸相对于背景过小。因而成像时，与广袤的背景相比，目标是像素级或亚像素级。那么，目标在最终成像中需要确定的，一是成像像素位置，二是通过计算得到的该像素位置上的辐照度值。

1.1 接口设计

目标仿真部分接口设计如表1所示。

1.2 确定目标成像位置

在OpenGL中进行投影时，我们设置相关的投影几何参数，包括视点位置、视场角大小、前后投影面位置。因此，我们能够快速求出投影面的解析形式。加之我们有了目标在场景中的具体三维坐标位置，那么在OpenGL中，我们能很容易的确定该目标三维位置对应于最终成像图像中的像素点坐标位置（xt，y）t。

表1 目标仿真接口设计

1.3 目标成像像素位置上辐照度Et计算

根据“点源”辐射关系，如式（1）所示[5]，由点的辐强度It得到其成像辐照度Et，我们需先计算目标成像角度θ及目标与成像位置的距离l。

对于θ的计算，我们首先计算投影面法矢。由于投影面解析形式已知，那么这便能很容易计算而出。连接目标最终的成像像素在三维坐标系中的位置与目标实际位置构成向量，其模即为l，也是二者之间的距离（也可用视点位置与目标实际位置之间距离近似），与投影面法矢的夹角便是目标的成像角度θ，利用单位向量点乘可直接求得cosθ。最终，利用式（1）便可求得目标对应最终成像像素位置上的辐照度Et。

2 背景辐射成像仿真及辐照度计算

在最终成像的画面上，显然背景应是成像画面的主要内容。背景成像借助导入的背景模型文件，是由三角形面元构成背景场景。因此，将背景分解为“小面辐射源”成像是十分恰当的，而下文中的所有计算也是针对每一个组成背景的面元的。

图1 背景辐射成像过程

本次仿真设计针对的是中红外波段的红外成像。因此，如图1所示[6]，探测器获取的红外辐射由两部分组成：一是场景中材质表面自身向外发出的辐射，图中标号为1，二是反射光源（太阳及其他普通光源）辐射，图中标号为2。在此我们暂忽略大气带来的影响（图中标号为3、4）。

2.1 接口设计

背景仿真部分接口设计如表2所示。

表2 背景仿真接口设计

2.2 成像链路的建立

2.2.1 背景材质自身辐射

我们知道，一切温度高于绝对零度的物体无时无刻不在产生着辐射能量，从而被红外探测器发现与识别。该能量与地物表面性质、温度有关[7]。在此，我们将地物视为灰体，那么其自身辐射计算便较为直观，在已知温度、材质、发射率、探测器探测波段的情况下，对普朗克公式进行波段积分[8]，便可以直接求得背景地物材质本身在一定波长范围内向外的辐射量，用辐出度M表示，单位为w·m-2。

式（2）中，ε为材质在当前波段上的发射率，λ1、λ2为探测器波段范围，c1为第一辐射常数，c1=3.741 8×10（8W·cm-2·μm4），c2为第二辐射常数，c2=1.438 8×104（μm·K），T是材质温度（K），λ为波长（μm）。由此，求材质自身向外辐射的辐亮度Le，在半球条件2π下，可如式（3）所示。

2.2.2 基于BRDF的反射光源辐射

而在实际的中波红外有光源照射的情况下，探测器能够拍摄到的场景所发出的辐亮度，显然不应该只有材质自身辐射的贡献，还应有其反射光源辐射而带来的贡献。在此，我们采用的方法是利用BRDF数据直接计算材质面元反射光源辐射的辐亮度。

图2 面元M处BRDF反射模型

如图2所示，BRDF被称作双向反射分布函数，用于描述表面入射光和反射光的关系，表征特定反射方向上离开表面的光辐亮度与入射光辐照度的比例，如式（4）所示。

其中，M表示入射的面元，Lr表示从（θr，φr）方向反射而出的辐亮度；E表示从（θi，φ）i方向入射的辐照度；θi为光源方向与面元法矢的夹角（即入射天顶角），θr为某一出射方向（即探测器方向与面元法矢夹角），φi为光源方向在面元上的投影与在面元上建立的局部坐标系某一坐标轴的夹角（即入射方位角），同理，φr为探测器方位角，λ为入射波长。由于局部坐标系建立的不确定性，φi与φr具体的大小我们不得而知，但无论怎么建立坐标系，二者之间的夹角不变。因此，定义φ=φr-φi，范围为[0，π]。背景材质在本仿真设计中视为各向同性，因而其BRDF具有旋转不变性，即无论φi与φr的大小如何改变，只要二者之差不变，对应的BRDF值就不变，从而式（4）中的fr便可以降维表示为f（rM，θi，θr，φi，λ）。

在理论长久发展和研究的基础上，一些学者采用不同的方法对材质的BRDF模型进行了研究，有些给出了BRDF具体的解析形式，例如文献[9]所提出的5参数模型[9]，这对我们来说是最理想的。但更多情况下，BRDF的解析形式无法获得，但不同波段下的BRDF实验数据库是可以得到的。该数据库以类似于索引的方式给出了不同输入参数组合情况下对应的BRDF具体数值，而本设计单元反射部分的计算基于就是它。我们只需通过已知的入射波长λ，通过计算得到θi、θr及φ，就能通过索引的方式，确定当前面元的BRDF值，从而计算得到其反射光源辐射的辐亮度贡献L（rM），如式（5）所示：

E(M,θi,φi,λ)为以（θi，φ）i角度入射的光源辐照度，fr(M,θi,θr,φ,λ)为入射波长是λ的情况下，以3个角度为输入索引而出的 BRDF 数值；Lr(M,θi,θr,φ,λ)即为在光源以（θi，φ）i角入射，探测器以（θr，φr）拍摄情况下，面元反射光源辐射进入到探测器中辐亮度。

2.2.3 成像链路模型建立

基于上两个小节，就可以完整地通过建立探测器入瞳时，背景的成像仿真链路模型，综合式（3）与式（5），同时，如果存在多个光源，得：

那么，在获取了探测器波段范围、光源入射波长、光源辐照度参数、面元BRDF参数后，我们便可以利用公式（6）的成像链路模型得到背景材质每一个面元上的辐亮度值Lout。

2.2.4 背景面元辐照度Eout计算

最终利用的是面元辐照度进行成像仿真，并用于度量最终所有成像位置上的辐射量。那么，对于作为“小面源”的背景面元在探测器入瞳镜头表面的成像，我们利用式（7）进行辐照度的计算。

其中，As为每一个背景面元的面积，单位m2，由于模型文件中给出了面元3个顶点位置，那么利用海伦公式便可以直接求得As具体数值，海伦公式如式（8）所示，a、b、c为三角形面元三边长，p为三角形半周长；ls为每一个背景面元中心点与视点连线向量的模。那么，每一个背景面元的辐亮度Lout便可通过式（7）转化为对应成像位置的辐照度Eout。

2.3 背景仿真成像及辐照度计算

基于以上原理，我们便可以进行相应的仿真设计，以下的内容也体现了我们本次设计为实时性红外仿真成像做出的努力之一。

首先，在OpenGL中导入模型文件，文件中包含顶点、面元、材质等信息，同时建立面元数据结构，求各面元单位法矢。

第二，利用给出的探测器方向，遍历各面元，求视线天顶角θr。如果小于90°，那么该面元就对探测器可见，即在成像后能够被我们所看到，我们只针对这些面元进行辐射计算。否则，不必对其进行计算。

第三，利用给出的光源方向，求可见面元的光源入射天顶角θi。同理，如果不小于90°，那么当前面元未被光源照到，只进行自身辐射Le的计算，成像链路只由自身辐射辐亮度构成；否则，面元能被光源照到，还要进行反射光源辐射环节Lr的计算，成像链路Lout由二者叠加组成，最后利用式（7）转换为对应的Eout。

第四，灰度量化，实现变面元辐亮度值为OpenGL中可利用的“颜色”（实为灰度）信息。OpenGL本身没有红外场景绘制功能，其只能通过面元为每个像元赋予一定的颜色灰度[10]来绘制可见光场景。因此，要想在OpenGL的成像方法上将面元辐亮度对像元灰度的影响（红外成像原理）体现出来，就需要建立辐亮度与颜色灰度之间的映射关系，由面元的辐亮度信息得到其灰度信息[11]，进而才能借助OpenGL实现红外场景的绘制及渲染。

由于最终要实现目标与背景的成像融合，因此，我们必须将目标成像位置上辐亮度数据Et与背景面元的辐亮度数据Eout同时量化，如式（9）所示。

将面元颜色的RGB 3个通道分量都赋为同一个值，就能得到相应的灰度图像[12]。

第五，利用OpenGL强大的绘制功能，以三角形面元为最小单元进行背景的场景绘制及渲染，得到整个红外背景场景。

最后，利用OpenGL的几何投影方法，将绘制好的场景转化为最终二维屏幕平面上显示的图像，完成最终的背景红外仿真成像，投影参数可按需设定。

2.4 成像辐照度数据获取

利用OpenGL自身功能，我们可以直接获取背景二维成像中每一个像素位置上的像素信息。那么将得到的灰度（红外灰度值），进行反量化操作，便可以获得每一个像素位置上的辐照度Efin数据，反量化公式是式（9）的逆过程，如式（10）所示。但由灰度还原回辐照度的过程中，存在精度损失，这也是本仿真中所存在的一个不足。最终，按照“（Px，Py）:Efin”的格式，将结果输出。

3 目标与背景的融合

至此，对目标和背景的红外仿真成像设计方法均已给出，但需要将二者融合，实现最终完成的成像及辐射数据输出，包含成像和数据两方面。

首先，成像融合。我们现已求得目标所在成像像素位置（xt，y）t、该位置上的目标成像灰度值gt，r，g，b，以及背景的红外仿真成像图。由于未加入目标的成像，因而，该背景画面中不含目标：画面中，目标所在成像像素位置上的灰度值还是背景的。那么，要在背景成像中体现目标成像情况，只需在背景的成像仿真图中，将目标所在像素位置上的灰度值（背景的灰度值）替换为gt，r，g，b即可。

其次，辐射数据融合。同理，要在辐照度数据中体现目标情况，只需在背景辐照度数据文件中，将像素（xt，y）t位置上的辐照度值替换为目标的辐照度Et即可。

4 基于CUDA的并行计算

4.1 CUDA并行思想原理

统一设备架构（CUDA）是NVIDIA公司提出的一个基于GPU通用计算的开发环境，它针对GPU多处理单元的特性，通过并行计算提高大规模运算的速度[13]。它包含了GPU的指令集及并行计算引擎，可用C语言来调用GPU硬件以并行的方式来解决数据巨大但结构类似的复杂计算问题[14]。

在CUDA架构下，一个程序分为Host端和Device端。Host端在CPU上执行，而Device端在GPU上执行，又称为“kernel”函数。通常Host端程序会将数据准备好后，复制到GPU，再由Device端并行执行，完成后再返回至CPU的Host端。

CUDA利用GPU实现并行。从硬件上来看，GPU中的SP（Streaming Processor）也称为CUDA核，是GPU最基本的处理单元，具体的指令和任务都是在SP上进行的。GPU进行并行计算，也就是多个SP同时工作。多个SP及其他资源构成一个SM（Streaming Multiprocessor），也称为 GPU 大核，是GPU的心脏，多个SM和其他资源构成整个GPU。

从软件上来看，一个CUDA程序会分成多个并行线程（thread）来执行；多个线程构成一个线程块（block）；多个block构成线程束（warp），是GPU执行时的最基本调度和运行单元。当一个kernel函数启动后，thread会被分配到这些SM的SP中执行，一个SP对应一个thread。大量的thread可能会被分到不同的SM中，但同一block中的thread必然在同一SM中执行。实际上并不是所有的thread能够在同一时刻执行。NVIDIA把32个thread组成一个warp。一个warp中所有的thread并行的执行相同的指令后，才转到下一warp。一个warp需要占用一个SM，多个warp轮流进入SM。

GPU并行执行将数据元素映射到并行的线程上，实现加速处理及计算。目前，利用GPU已实现三维图像渲染，图像和多媒体应用处理、视频编解码、立体视觉和模式识别等过程[15]。

4.2 利用CUDA实现红外仿真成像

在上述成像仿真的数据计算过程中，我们利用CUDA实现了CPU上的串行到GPU上并行编程的移植，充分利用了现代GPU的强大数据处理能力[16]，也是为实时性计算及仿真做出的努力之二。

在此，我们所做的就是利用CUDA将面元映射到GPU上进行相关计算，实现实时红外仿真成像。具体体现在，一是在遍历面元求单位法矢时；二是为确定面元BRDF值而对光源天顶角θi、探测器天顶角θr、方位角φ进行计算时；三是求每个面元面积As及与投影点距离ls时；四是通过以上面元数据并由3个角度索引出其BRDF值后，由成像链路模型计算得到面元入瞳总辐亮度时——均是在GPU上实现的。

将这些数据的计算由CPU上的串行转移到GPU上利用核函数进行并行计算，只需在开发环境中，按需求定义多个kernel函数，同时做好CPU、GPU中内存空间的开辟及二者之间接口数据的传递即可。

5 结束语

文中提出了一种利用CUDA并行思想，在GPU上实现的实时中波红外辐射传输计算及成像仿真方法，对目标和背景的成像仿真进行了原理分析及具体设计。对目标采用“点源”成像思路进行计算及设计；对于背景采用分解为“小面源”的思路。同时建立不同于可见光成像原理的红外入瞳成像仿真链路模型，包含背景材质自身辐射及材质反射光源辐射两大部分。其中在反射光源辐射部分，特别提出了利用BRDF模型数据计算反射辐射量的方法。本设计对红外仿真成像的准确性和实时性均会带来很大的提高。

[1]王长胜.典型地物表面温度特性研究[D].大连:大连海事大学，2013.

[2]周方方，司斌，王晖，等.红外图像仿真的发展综述[J].航空兵器，2013（2）:29-32.

[3]陈珊，孙继银，李琳琳，等.建筑物目标三维红外场景仿真[J].激光与红外，2010，40（5）:500-505.

[4]张军华，臧胜涛，单联瑜，等.高性能计算的发展现状及趋势[J].石油地球物理勘探，2010，45（6）:918-925.

[5]杨风暴.红外物理与技术[M].北京:电子工业出版社，2014:36-39.

[6]季园园.典型地物红外特性仿真关键技术研究[D].大连:大连海事大学，2014.

[7]胡海鹤.红外视景仿真关键技术研究[D].北京:北京理工大学，2015.

[8]王霞，汪昊，徐超，等.红外场景仿真技术发展综述[J].红外技术，2015，37（7）:537-543.

[9]谢东辉.目标与地物背景光散射特性建模[D].西安:西安电子科技大学，2002.

[10]方宁，张晓林，高嘉，等.OpenGL在红外辐射计算中的应用[J].北京航空航天大学学报，2009（12）:1473-1476.

[11]娄树理，杨咏东，周晓东.基于OpenGL红外视景仿真研究[J].烟台大学学报:自然科学与工程版，2004，17（2）:129-134.

[12]黄迁，张涛，吕银环，等.基于 Visual C++和OpenGL的红外目标景象仿真研究[J].红外技术，2010（2）:101-104.

[13]董荦，葛万成，陈康力.CUDA并行计算的应用研究[J].信息技术，2010（4）:11-15.

[14]库克，苏统华.CUDA并行程序设计:GPU编程指南[M].北京:机械工业出版社，2014.

[15]Nicholas Wilt.CUDA专家手册:GPU编程权威指南:The CUDA handbook:a comprehensive guide to GPU programming[M].北京:机械工业出版社，2014.

[16]吴永香.基于CUDA的红外成像系统虚拟样机研究[D].西安:西安电子科技大学，2012.