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

（中國科學(xué)院國家空間科學(xué)中心北京100190）

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

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

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

1 目標(biāo)輻射成像仿真及輻照度計算

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

1.1 接口設(shè)計

目標(biāo)仿真部分接口設(shè)計如表1所示。

1.2 確定目標(biāo)成像位置

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

表1 目標(biāo)仿真接口設(shè)計

1.3 目標(biāo)成像像素位置上輻照度Et計算

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

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

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

在最終成像的畫面上，顯然背景應(yīng)是成像畫面的主要內(nèi)容。背景成像借助導(dǎo)入的背景模型文件，是由三角形面元構(gòu)成背景場景。因此，將背景分解為“小面輻射源”成像是十分恰當(dāng)?shù)模挛闹械乃杏嬎阋彩轻槍γ恳粋€組成背景的面元的。

圖1 背景輻射成像過程

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

2.1 接口設(shè)計

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

表2 背景仿真接口設(shè)計

2.2 成像鏈路的建立

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

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

式（2）中，ε為材質(zhì)在當(dāng)前波段上的發(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ā)出的輻亮度，顯然不應(yīng)該只有材質(zhì)自身輻射的貢獻，還應(yīng)有其反射光源輻射而帶來的貢獻。在此，我們采用的方法是利用BRDF數(shù)據(jù)直接計算材質(zhì)面元反射光源輻射的輻亮度。

圖2 面元M處BRDF反射模型

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

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

在理論長久發(fā)展和研究的基礎(chǔ)上，一些學(xué)者采用不同的方法對材質(zhì)的BRDF模型進行了研究，有些給出了BRDF具體的解析形式，例如文獻[9]所提出的5參數(shù)模型[9]，這對我們來說是最理想的。但更多情況下，BRDF的解析形式無法獲得，但不同波段下的BRDF實驗數(shù)據(jù)庫是可以得到的。該數(shù)據(jù)庫以類似于索引的方式給出了不同輸入?yún)?shù)組合情況下對應(yīng)的BRDF具體數(shù)值，而本設(shè)計單元反射部分的計算基于就是它。我們只需通過已知的入射波長λ，通過計算得到θi、θr及φ，就能通過索引的方式，確定當(dāng)前面元的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)化為對應(yīng)成像位置的輻照度Eout。

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

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

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

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

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

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

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

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

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

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

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

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

3 目標(biāo)與背景的融合

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

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

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

4 基于CUDA的并行計算

4.1 CUDA并行思想原理

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

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

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

從軟件上來看，一個CUDA程序會分成多個并行線程（thread）來執(zhí)行；多個線程構(gòu)成一個線程塊（block）；多個block構(gòu)成線程束（warp），是GPU執(zhí)行時的最基本調(diào)度和運行單元。當(dāng)一個kernel函數(shù)啟動后，thread會被分配到這些SM的SP中執(zhí)行，一個SP對應(yīng)一個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)三維圖像渲染，圖像和多媒體應(yīng)用處理、視頻編解碼、立體視覺和模式識別等過程[15]。

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

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

在此，我們所做的就是利用CUDA將面元映射到GPU上進行相關(guān)計算，實現(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 結(jié)束語

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

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

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

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

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

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

[6]季園園.典型地物紅外特性仿真關(guān)鍵技術(shù)研究[D].大連:大連海事大學(xué)，2014.

[7]胡海鶴.紅外視景仿真關(guān)鍵技術(shù)研究[D].北京:北京理工大學(xué)，2015.

[8]王霞，汪昊，徐超，等.紅外場景仿真技術(shù)發(fā)展綜述[J].紅外技術(shù)，2015，37（7）:537-543.

[9]謝東輝.目標(biāo)與地物背景光散射特性建模[D].西安:西安電子科技大學(xué)，2002.

[10]方寧，張曉林，高嘉，等.OpenGL在紅外輻射計算中的應(yīng)用[J].北京航空航天大學(xué)學(xué)報，2009（12）:1473-1476.

[11]婁樹理，楊詠東，周曉東.基于OpenGL紅外視景仿真研究[J].煙臺大學(xué)學(xué)報:自然科學(xué)與工程版，2004，17（2）:129-134.

[12]黃遷，張濤，呂銀環(huán)，等.基于 Visual C++和OpenGL的紅外目標(biāo)景象仿真研究[J].紅外技術(shù)，2010（2）:101-104.

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

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

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

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