蔡方偉，卞春江，馮水春，曹景，毛博年

（中國科學院國家空間科學中心北京100190）

近幾年來，科學研究對不同天候、溫度、場景、干 擾情況下的紅外圖像的需求越來越大。但多數(shù)情況下，外場實拍不僅要投入大量的人力和物力而且還要耗費相當長的時間[1]；同時，該行為往往與軍事目的有關，因而很可能根本無法進行，且無法窮盡所有天候條件。在該情況下，通過計算機仿真手段獲取紅外圖像就顯得方便而快捷，并可應用于軍事、農(nóng)業(yè)、地質(zhì)、氣象等多個領域[2]。

本文采用計算機數(shù)學仿真技術，利用OpenGL強大的圖形算法，內(nèi)容包括目標、背景的紅外輻射傳輸及成像仿真以及二者的融合三部分。針對目標與背景在成像上的區(qū)別，我們將目標視為點輻射源的成像，而將背景分解為小面輻射源成像。對于目標，先利用輻射量計算公式直接計算成像輻射數(shù)據(jù)，再進行相應的成像顯示。對于背景，在OpenGL中加載背景三維場景模型[3]，建立面元的紅外成像鏈路，并進行場景模型繪制及渲染，通過合適的幾何投影變換得到其二維圖像。最終融合計算中，在背景成像上疊加目標成像影響，完成本文的設計任務。

需要注意的是，目前，不僅是在紅外成像仿真方面，在計算機軟件開發(fā)上的主流思想仍是CPU上的串行思想。隨著計算機硬件，尤其是各廠商制造的GPU水平的提高，圖形圖像處理能力大幅進步，處理的面元與像元數(shù)據(jù)量也因此越來越大，但其計算、處理方法卻又是類似的。如果采用按部就班的串行思想開發(fā)，不僅無法實現(xiàn)快速處理，也是對硬件資源的極大浪費。而GPU上的并行思想顯然給這個問題的解決提供了新的思路。而CUDA便是目前實現(xiàn)GPU上并行處理最成熟、實用化程度最高的技術[4]。

1 目標輻射成像仿真及輻照度計算

上文提過，將目標視為“點輻射源”。其原因在于實際中，目標尺寸相對于背景過小。因而成像時，與廣袤的背景相比，目標是像素級或亞像素級。那么，目標在最終成像中需要確定的，一是成像像素位置，二是通過計算得到的該像素位置上的輻照度值。

1.1 接口設計

目標仿真部分接口設計如表1所示。

1.2 確定目標成像位置

在OpenGL中進行投影時，我們設置相關的投影幾何參數(shù)，包括視點位置、視場角大小、前后投影面位置。因此，我們能夠快速求出投影面的解析形式。加之我們有了目標在場景中的具體三維坐標位置，那么在OpenGL中，我們能很容易的確定該目標三維位置對應于最終成像圖像中的像素點坐標位置（xt，y）t。

表1 目標仿真接口設計

1.3 目標成像像素位置上輻照度Et計算

根據(jù)“點源”輻射關系，如式（1）所示[5]，由點的輻強度It得到其成像輻照度Et，我們需先計算目標成像角度θ及目標與成像位置的距離l。

對于θ的計算，我們首先計算投影面法矢。由于投影面解析形式已知，那么這便能很容易計算而出。連接目標最終的成像像素在三維坐標系中的位置與目標實際位置構成向量，其模即為l，也是二者之間的距離（也可用視點位置與目標實際位置之間距離近似），與投影面法矢的夾角便是目標的成像角度θ，利用單位向量點乘可直接求得cosθ。最終，利用式（1）便可求得目標對應最終成像像素位置上的輻照度Et。

2 背景輻射成像仿真及輻照度計算

在最終成像的畫面上，顯然背景應是成像畫面的主要內(nèi)容。背景成像借助導入的背景模型文件，是由三角形面元構成背景場景。因此，將背景分解為“小面輻射源”成像是十分恰當?shù)?，而下文中的所有計算也是針對每一個組成背景的面元的。

圖1 背景輻射成像過程

本次仿真設計針對的是中紅外波段的紅外成像。因此，如圖1所示[6]，探測器獲取的紅外輻射由兩部分組成：一是場景中材質(zhì)表面自身向外發(fā)出的輻射，圖中標號為1，二是反射光源（太陽及其他普通光源）輻射，圖中標號為2。在此我們暫忽略大氣帶來的影響（圖中標號為3、4）。

2.1 接口設計

背景仿真部分接口設計如表2所示。

表2 背景仿真接口設計

2.2 成像鏈路的建立

2.2.1 背景材質(zhì)自身輻射

我們知道，一切溫度高于絕對零度的物體無時無刻不在產(chǎn)生著輻射能量，從而被紅外探測器發(fā)現(xiàn)與識別。該能量與地物表面性質(zhì)、溫度有關[7]。在此，我們將地物視為灰體，那么其自身輻射計算便較為直觀，在已知溫度、材質(zhì)、發(fā)射率、探測器探測波段的情況下，對普朗克公式進行波段積分[8]，便可以直接求得背景地物材質(zhì)本身在一定波長范圍內(nèi)向外的輻射量，用輻出度M表示，單位為w?m-2。

式（2）中，ε為材質(zhì)在當前波段上的發(fā)射率，λ1、λ2為探測器波段范圍，c1為第一輻射常數(shù)，c1=3.741 8×10（8W·cm-2·μm4），c2為第二輻射常數(shù)，c2=1.438 8×104（μm·K），T是材質(zhì)溫度（K），λ為波長（μm）。由此，求材質(zhì)自身向外輻射的輻亮度Le，在半球條件2π下，可如式（3）所示。

2.2.2 基于BRDF的反射光源輻射

而在實際的中波紅外有光源照射的情況下，探測器能夠拍攝到的場景所發(fā)出的輻亮度，顯然不應該只有材質(zhì)自身輻射的貢獻，還應有其反射光源輻射而帶來的貢獻。在此，我們采用的方法是利用BRDF數(shù)據(jù)直接計算材質(zhì)面元反射光源輻射的輻亮度。

圖2 面元M處BRDF反射模型

如圖2所示，BRDF被稱作雙向反射分布函數(shù)，用于描述表面入射光和反射光的關系，表征特定反射方向上離開表面的光輻亮度與入射光輻照度的比例，如式（4）所示。

其中，M表示入射的面元，Lr表示從（θr，φr）方向反射而出的輻亮度；E表示從（θi，φ）i方向入射的輻照度；θi為光源方向與面元法矢的夾角（即入射天頂角），θr為某一出射方向（即探測器方向與面元法矢夾角），φi為光源方向在面元上的投影與在面元上建立的局部坐標系某一坐標軸的夾角（即入射方位角），同理，φr為探測器方位角，λ為入射波長。由于局部坐標系建立的不確定性，φi與φr具體的大小我們不得而知，但無論怎么建立坐標系，二者之間的夾角不變。因此，定義φ=φr-φi，范圍為[0，π]。背景材質(zhì)在本仿真設計中視為各向同性，因而其BRDF具有旋轉(zhuǎn)不變性，即無論φi與φr的大小如何改變，只要二者之差不變，對應的BRDF值就不變，從而式（4）中的fr便可以降維表示為f（rM，θi，θr，φi，λ）。

在理論長久發(fā)展和研究的基礎上，一些學者采用不同的方法對材質(zhì)的BRDF模型進行了研究，有些給出了BRDF具體的解析形式，例如文獻[9]所提出的5參數(shù)模型[9]，這對我們來說是最理想的。但更多情況下，BRDF的解析形式無法獲得，但不同波段下的BRDF實驗數(shù)據(jù)庫是可以得到的。該數(shù)據(jù)庫以類似于索引的方式給出了不同輸入?yún)?shù)組合情況下對應的BRDF具體數(shù)值，而本設計單元反射部分的計算基于就是它。我們只需通過已知的入射波長λ，通過計算得到θi、θr及φ，就能通過索引的方式，確定當前面元的BRDF值，從而計算得到其反射光源輻射的輻亮度貢獻L（rM），如式（5）所示：

E(M,θi,φi,λ)為以（θi，φ）i角度入射的光源輻照度，fr(M,θi,θr,φ,λ)為入射波長是λ的情況下，以3個角度為輸入索引而出的 BRDF 數(shù)值；Lr(M,θi,θr,φ,λ)即為在光源以（θi，φ）i角入射，探測器以（θr，φr）拍攝情況下，面元反射光源輻射進入到探測器中輻亮度。

2.2.3 成像鏈路模型建立

基于上兩個小節(jié)，就可以完整地通過建立探測器入瞳時，背景的成像仿真鏈路模型，綜合式（3）與式（5），同時，如果存在多個光源，得：

那么，在獲取了探測器波段范圍、光源入射波長、光源輻照度參數(shù)、面元BRDF參數(shù)后，我們便可以利用公式（6）的成像鏈路模型得到背景材質(zhì)每一個面元上的輻亮度值Lout。

2.2.4 背景面元輻照度Eout計算

最終利用的是面元輻照度進行成像仿真，并用于度量最終所有成像位置上的輻射量。那么，對于作為“小面源”的背景面元在探測器入瞳鏡頭表面的成像，我們利用式（7）進行輻照度的計算。

其中，As為每一個背景面元的面積，單位m2，由于模型文件中給出了面元3個頂點位置，那么利用海倫公式便可以直接求得As具體數(shù)值，海倫公式如式（8）所示，a、b、c為三角形面元三邊長，p為三角形半周長；ls為每一個背景面元中心點與視點連線向量的模。那么，每一個背景面元的輻亮度Lout便可通過式（7）轉(zhuǎn)化為對應成像位置的輻照度Eout。

2.3 背景仿真成像及輻照度計算

基于以上原理，我們便可以進行相應的仿真設計，以下的內(nèi)容也體現(xiàn)了我們本次設計為實時性紅外仿真成像做出的努力之一。

首先，在OpenGL中導入模型文件，文件中包含頂點、面元、材質(zhì)等信息，同時建立面元數(shù)據(jù)結構，求各面元單位法矢。

第二，利用給出的探測器方向，遍歷各面元，求視線天頂角θr。如果小于90°，那么該面元就對探測器可見，即在成像后能夠被我們所看到，我們只針對這些面元進行輻射計算。否則，不必對其進行計算。

第三，利用給出的光源方向，求可見面元的光源入射天頂角θi。同理，如果不小于90°，那么當前面元未被光源照到，只進行自身輻射Le的計算，成像鏈路只由自身輻射輻亮度構成；否則，面元能被光源照到，還要進行反射光源輻射環(huán)節(jié)Lr的計算，成像鏈路Lout由二者疊加組成，最后利用式（7）轉(zhuǎn)換為對應的Eout。

第四，灰度量化，實現(xiàn)變面元輻亮度值為OpenGL中可利用的“顏色”（實為灰度）信息。OpenGL本身沒有紅外場景繪制功能，其只能通過面元為每個像元賦予一定的顏色灰度[10]來繪制可見光場景。因此，要想在OpenGL的成像方法上將面元輻亮度對像元灰度的影響（紅外成像原理）體現(xiàn)出來，就需要建立輻亮度與顏色灰度之間的映射關系，由面元的輻亮度信息得到其灰度信息[11]，進而才能借助OpenGL實現(xiàn)紅外場景的繪制及渲染。

由于最終要實現(xiàn)目標與背景的成像融合，因此，我們必須將目標成像位置上輻亮度數(shù)據(jù)Et與背景面元的輻亮度數(shù)據(jù)Eout同時量化，如式（9）所示。

將面元顏色的RGB 3個通道分量都賦為同一個值，就能得到相應的灰度圖像[12]。

第五，利用OpenGL強大的繪制功能，以三角形面元為最小單元進行背景的場景繪制及渲染，得到整個紅外背景場景。

最后，利用OpenGL的幾何投影方法，將繪制好的場景轉(zhuǎn)化為最終二維屏幕平面上顯示的圖像，完成最終的背景紅外仿真成像，投影參數(shù)可按需設定。

2.4 成像輻照度數(shù)據(jù)獲取

利用OpenGL自身功能，我們可以直接獲取背景二維成像中每一個像素位置上的像素信息。那么將得到的灰度（紅外灰度值），進行反量化操作，便可以獲得每一個像素位置上的輻照度Efin數(shù)據(jù)，反量化公式是式（9）的逆過程，如式（10）所示。但由灰度還原回輻照度的過程中，存在精度損失，這也是本仿真中所存在的一個不足。最終，按照“（Px，Py）:Efin”的格式，將結果輸出。

3 目標與背景的融合

至此，對目標和背景的紅外仿真成像設計方法均已給出，但需要將二者融合，實現(xiàn)最終完成的成像及輻射數(shù)據(jù)輸出，包含成像和數(shù)據(jù)兩方面。

首先，成像融合。我們現(xiàn)已求得目標所在成像像素位置（xt，y）t、該位置上的目標成像灰度值gt，r，g，b，以及背景的紅外仿真成像圖。由于未加入目標的成像，因而，該背景畫面中不含目標：畫面中，目標所在成像像素位置上的灰度值還是背景的。那么，要在背景成像中體現(xiàn)目標成像情況，只需在背景的成像仿真圖中，將目標所在像素位置上的灰度值（背景的灰度值）替換為gt，r，g，b即可。

其次，輻射數(shù)據(jù)融合。同理，要在輻照度數(shù)據(jù)中體現(xiàn)目標情況，只需在背景輻照度數(shù)據(jù)文件中，將像素（xt，y）t位置上的輻照度值替換為目標的輻照度Et即可。

4 基于CUDA的并行計算

4.1 CUDA并行思想原理

統(tǒng)一設備架構（CUDA）是NVIDIA公司提出的一個基于GPU通用計算的開發(fā)環(huán)境，它針對GPU多處理單元的特性，通過并行計算提高大規(guī)模運算的速度[13]。它包含了GPU的指令集及并行計算引擎，可用C語言來調(diào)用GPU硬件以并行的方式來解決數(shù)據(jù)巨大但結構類似的復雜計算問題[14]。

在CUDA架構下，一個程序分為Host端和Device端。Host端在CPU上執(zhí)行，而Device端在GPU上執(zhí)行，又稱為“kernel”函數(shù)。通常Host端程序會將數(shù)據(jù)準備好后，復制到GPU，再由Device端并行執(zhí)行，完成后再返回至CPU的Host端。

CUDA利用GPU實現(xiàn)并行。從硬件上來看，GPU中的SP（Streaming Processor）也稱為CUDA核，是GPU最基本的處理單元，具體的指令和任務都是在SP上進行的。GPU進行并行計算，也就是多個SP同時工作。多個SP及其他資源構成一個SM（Streaming Multiprocessor），也稱為 GPU 大核，是GPU的心臟，多個SM和其他資源構成整個GPU。

從軟件上來看，一個CUDA程序會分成多個并行線程（thread）來執(zhí)行；多個線程構成一個線程塊（block）；多個block構成線程束（warp），是GPU執(zhí)行時的最基本調(diào)度和運行單元。當一個kernel函數(shù)啟動后，thread會被分配到這些SM的SP中執(zhí)行，一個SP對應一個thread。大量的thread可能會被分到不同的SM中，但同一block中的thread必然在同一SM中執(zhí)行。實際上并不是所有的thread能夠在同一時刻執(zhí)行。NVIDIA把32個thread組成一個warp。一個warp中所有的thread并行的執(zhí)行相同的指令后，才轉(zhuǎn)到下一warp。一個warp需要占用一個SM，多個warp輪流進入SM。

GPU并行執(zhí)行將數(shù)據(jù)元素映射到并行的線程上，實現(xiàn)加速處理及計算。目前，利用GPU已實現(xiàn)三維圖像渲染，圖像和多媒體應用處理、視頻編解碼、立體視覺和模式識別等過程[15]。

4.2 利用CUDA實現(xiàn)紅外仿真成像

在上述成像仿真的數(shù)據(jù)計算過程中，我們利用CUDA實現(xiàn)了CPU上的串行到GPU上并行編程的移植，充分利用了現(xiàn)代GPU的強大數(shù)據(jù)處理能力[16]，也是為實時性計算及仿真做出的努力之二。

在此，我們所做的就是利用CUDA將面元映射到GPU上進行相關計算，實現(xiàn)實時紅外仿真成像。具體體現(xiàn)在，一是在遍歷面元求單位法矢時；二是為確定面元BRDF值而對光源天頂角θi、探測器天頂角θr、方位角φ進行計算時；三是求每個面元面積As及與投影點距離ls時；四是通過以上面元數(shù)據(jù)并由3個角度索引出其BRDF值后，由成像鏈路模型計算得到面元入瞳總輻亮度時——均是在GPU上實現(xiàn)的。

將這些數(shù)據(jù)的計算由CPU上的串行轉(zhuǎn)移到GPU上利用核函數(shù)進行并行計算，只需在開發(fā)環(huán)境中，按需求定義多個kernel函數(shù)，同時做好CPU、GPU中內(nèi)存空間的開辟及二者之間接口數(shù)據(jù)的傳遞即可。

5 結束語

文中提出了一種利用CUDA并行思想，在GPU上實現(xiàn)的實時中波紅外輻射傳輸計算及成像仿真方法，對目標和背景的成像仿真進行了原理分析及具體設計。對目標采用“點源”成像思路進行計算及設計；對于背景采用分解為“小面源”的思路。同時建立不同于可見光成像原理的紅外入瞳成像仿真鏈路模型，包含背景材質(zhì)自身輻射及材質(zhì)反射光源輻射兩大部分。其中在反射光源輻射部分，特別提出了利用BRDF模型數(shù)據(jù)計算反射輻射量的方法。本設計對紅外仿真成像的準確性和實時性均會帶來很大的提高。

[1]王長勝.典型地物表面溫度特性研究[D].大連:大連海事大學，2013.

[2]周方方，司斌，王暉，等.紅外圖像仿真的發(fā)展綜述[J].航空兵器，2013（2）:29-32.

[3]陳珊，孫繼銀，李琳琳，等.建筑物目標三維紅外場景仿真[J].激光與紅外，2010，40（5）:500-505.

[4]張軍華，臧勝濤，單聯(lián)瑜，等.高性能計算的發(fā)展現(xiàn)狀及趨勢[J].石油地球物理勘探，2010，45（6）:918-925.

[5]楊風暴.紅外物理與技術[M].北京:電子工業(yè)出版社，2014:36-39.

[6]季園園.典型地物紅外特性仿真關鍵技術研究[D].大連:大連海事大學，2014.

[7]胡海鶴.紅外視景仿真關鍵技術研究[D].北京:北京理工大學，2015.

[8]王霞，汪昊，徐超，等.紅外場景仿真技術發(fā)展綜述[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]黃遷，張濤，呂銀環(huán)，等.基于 Visual C++和OpenGL的紅外目標景象仿真研究[J].紅外技術，2010（2）:101-104.

[13]董犖，葛萬成，陳康力.CUDA并行計算的應用研究[J].信息技術，2010（4）:11-15.

[14]庫克，蘇統(tǒng)華.CUDA并行程序設計:GPU編程指南[M].北京:機械工業(yè)出版社，2014.

[15]Nicholas Wilt.CUDA專家手冊:GPU編程權威指南:The CUDA handbook:a comprehensive guide to GPU programming[M].北京:機械工業(yè)出版社，2014.

[16]吳永香.基于CUDA的紅外成像系統(tǒng)虛擬樣機研究[D].西安:西安電子科技大學，2012.