王彩云,趙煥玥,李曉飛,劉純勝,王立國(guó),王佳寧

(1.南京航空航天大學(xué)航天學(xué)院,江蘇 南京 210001;2.北京電子工程總體研究所,北京 100854)

1 引 言

由于航天技術(shù)的發(fā)展,各軍事強(qiáng)國(guó)發(fā)射進(jìn)入太空的空間目標(biāo)數(shù)量劇增。空間目標(biāo)的體積愈加小型化,結(jié)構(gòu)更加復(fù)雜化,形狀越來(lái)越多樣化,空間目標(biāo)探測(cè)被賦予新的挑戰(zhàn)。紅外探測(cè)作為一種重要的光學(xué)探測(cè)手段,因其對(duì)目標(biāo)背景溫差的高敏感性,探測(cè)的高隱蔽性以及可以進(jìn)行多光譜探測(cè)的優(yōu)點(diǎn)成為了理想的空間探測(cè)方式[1]。與外場(chǎng)試驗(yàn)相比較,紅外成像仿真可以獲得特定位置、時(shí)間和天氣下的紅外圖像,不僅能減少試驗(yàn)成本,還能提高試驗(yàn)效率,具有廣泛的應(yīng)用前景與研究?jī)r(jià)值。

國(guó)內(nèi)在空間目標(biāo)的紅外成像起步較晚,但是也取得了很大進(jìn)步。例如文獻(xiàn)[2]～[4]中闡述了空間目標(biāo)的紅外輻射特性分析建模方法,對(duì)其紅外輻射強(qiáng)度進(jìn)行了數(shù)值求解,并將計(jì)算結(jié)果轉(zhuǎn)化成灰度值,最后采用計(jì)算機(jī)圖形學(xué)算法生成紅外圖像,為空間目標(biāo)紅外探測(cè)提供了有力的工具。但是在紅外成像方面,已有的研究大都基于理論數(shù)值計(jì)算階段,尚未開(kāi)發(fā)一套從目標(biāo)表面溫度分布求解到目標(biāo)紅外輻射特性分析求解,再轉(zhuǎn)入目標(biāo)紅外成像的系統(tǒng)仿真軟件,而國(guó)外成熟的紅外成像商業(yè)軟件由于版權(quán)限制,價(jià)格昂貴,也不是最佳選擇[5]。

所以基于經(jīng)濟(jì)因素和可行性的綜合考慮,本文提出了一種基于三角剖分的空間目標(biāo)紅外成像仿真方法,建立了一套完整的空間目標(biāo)紅外成像仿真方法。此方法首先根據(jù)空間目標(biāo)熱平衡方程[6],利用節(jié)點(diǎn)網(wǎng)絡(luò)法求解出目標(biāo)的各節(jié)點(diǎn)溫度,然后利用STK得出探測(cè)器、目標(biāo)、地球與太陽(yáng)的位置關(guān)系,隨后由目標(biāo)節(jié)點(diǎn)溫度求得目標(biāo)紅外輻射強(qiáng)度,最后利用改進(jìn)Delaunay三角剖分的方法重建目標(biāo)得到紅外圖像。以某衛(wèi)星[7]為例進(jìn)行了成像仿真實(shí)驗(yàn),成像結(jié)果驗(yàn)證了本文方法的有效性,為進(jìn)一步的空間目標(biāo)的探測(cè)與識(shí)別提供了基礎(chǔ)。

2 空間目標(biāo)表面溫度求解模型

2.1 理論計(jì)算模型

空間目標(biāo)的紅外輻射特性主要取決于目標(biāo)的表面溫度分布。要對(duì)空間目標(biāo)進(jìn)行紅外成像,首先應(yīng)該建立目標(biāo)表面溫度求解模型,然后根據(jù)溫度分布求解其紅外輻射特性。

為了簡(jiǎn)化計(jì)算,忽略空間目標(biāo)自身工作的發(fā)熱,目標(biāo)表面溫度的分布主要由太陽(yáng)直接輻射[8]、地球反射太陽(yáng)輻射[9]、地球直接輻射和目標(biāo)表面向宇宙空間的輻射四個(gè)因素決定。

(1)太陽(yáng)直接輻射

太陽(yáng)的半徑為6.960×105km,可以被看作是溫度為5900 K的黑體。把太陽(yáng)與地球的距離當(dāng)作其平均距離(1.4961×108km),地球大氣層外法線(xiàn)方向的太陽(yáng)總輻照度Es為:

(1)

式中,Eλ是太陽(yáng)的光譜輻射照度。

另外,因?yàn)榈厍虬霃竭h(yuǎn)遠(yuǎn)小于太陽(yáng)與地球距離,可以將太陽(yáng)光視為輻射度均勻的平行光,則在目標(biāo)表面的第i個(gè)單元吸收的太陽(yáng)直接輻射qi1為:

qi1=εsEsSiFs,i

(2)

式中,εs表示材料對(duì)太陽(yáng)輻射的吸收率;Si為表面面積;Fs,i是外表面對(duì)于太陽(yáng)直接輻射的角系數(shù)。

(2)地球反射太陽(yáng)輻射

除了太陽(yáng)直接輻射,目標(biāo)表面還會(huì)接收地球反射太陽(yáng)輻射,被視為地球表面的漫反射,則該表面單元所吸收的地球反射太陽(yáng)輻射qi2為:

qi2=εrEsρFse,i

(3)

式中,εr為吸收率;ρ為地球反射率;Fse,i為角系數(shù)。

(3)地球直接輻射

地球可以被視為一個(gè)溫度280 K的黑體,其空間某處的總輻射照度Ee為:

(4)

式中,Re為 6450 km,是地球及其大氣系統(tǒng)總半徑;h為目標(biāo)和大氣層的距離,則該表面單元吸收的地球直接輻射qi3為:

qi3=εeEeSiFe,i

(5)

式中,εe為吸收率；Fe,i為角系數(shù)。

(4)目標(biāo)向宇宙空間的輻射

若目標(biāo)的溫度不低于熱力學(xué)零度,則目標(biāo)表面一定會(huì)向外輻射能量。第i個(gè)單元表面向宇宙空間的輻射能qi4為:

qi4=εσTi4(t)Si

(6)

式中,ε為發(fā)射率;σ為斯蒂芬-玻爾茲曼常數(shù);Ti(t) 為該部分在時(shí)間t的溫度。

(5)內(nèi)能的增加

空間目標(biāo)的內(nèi)能會(huì)根據(jù)溫度的變化而發(fā)生變化,在任一瞬間,第i個(gè)單元內(nèi)能的增加量為:

(7)

式中,c表示表面的比熱容;mi為質(zhì)量。

2.2 熱平衡方程的求解

根據(jù)上文的描述,熱平衡方程如下:

qi1+qi2+qi3-qi4=qi5

(8)

(9)

(2)將高次項(xiàng)變?yōu)橐淮雾?xiàng):

Ti4(t)=4Ti3(t-Δt)Ti(t)-3Ti4(t-Δt)

(10)

如果給定太陽(yáng)、目標(biāo)、地球的位置關(guān)系,就能計(jì)算出式(8)的角系數(shù)Fs,i、Fse,i和Fe,i的值,然后由式(8)～(10),可知對(duì)目標(biāo)任意表面單元,只要知道其在t-Δt時(shí)刻的溫度值Ti(t-Δt),其在t時(shí)刻的溫度值Ti(t)就可以推算得到。所以,若知道了目標(biāo)表面的初始溫度分布以及目標(biāo)與太陽(yáng)、地球的方位,目標(biāo)的表面溫度分布就能求出。

本文采用節(jié)點(diǎn)網(wǎng)絡(luò)法對(duì)空間目標(biāo)的表面溫度分布進(jìn)行計(jì)算。節(jié)點(diǎn)網(wǎng)絡(luò)法是把空間表面劃分為許多熱力學(xué)特性基本一致的單元,把單元作為節(jié)點(diǎn),對(duì)節(jié)點(diǎn)和背景的熱交換進(jìn)行分析,建立各節(jié)點(diǎn)的瞬時(shí)熱平衡方程。在給出的初始條件下,遞推計(jì)算這些方程組,便可以計(jì)算出目標(biāo)表面節(jié)點(diǎn)的溫度值和其變化規(guī)律[9]。

3 空間目標(biāo)的紅外輻射特性計(jì)算

空間目標(biāo)對(duì)空間某處的紅外輻射包含兩部分:目標(biāo)自身紅外輻射和目標(biāo)反射背景輻射。

(1)目標(biāo)自身紅外輻射

在空間目標(biāo)表面溫度分布被計(jì)算出以后,可以采用普朗克公式,在其紅外波段范圍的積分得到:

(11)

式中,λ1,λ2為紅外波段范圍的上下限;C1為第一輻射常數(shù),其值為3.742×108W·μm4/m2;T為單元表面溫度;C2為第二輻射常數(shù),其值為1.439×104μm·K。

目標(biāo)自身紅外輻射強(qiáng)度采用漫輻射源的輻射特性得出:

(12)

式中,ΔA為目標(biāo)表面單元面積;θi為目標(biāo)表面單元法線(xiàn)與輻射方向的夾角。

(2)目標(biāo)反射輻射

目標(biāo)對(duì)背景的反射輻射亮度為:

(13)

式中,ρe為等效反射系數(shù);E為外界對(duì)目標(biāo)輻射量的總和。目標(biāo)反射強(qiáng)度可采用漫反射源的輻射特性得到:

Ir=∑LrΔAcosθi

(14)

式中,ΔA為目標(biāo)表面單元面積;θi為目標(biāo)表面單元面法線(xiàn)與反射方向的夾角。

因此,空間目標(biāo)的總輻射強(qiáng)度為:

I=Ifa+Ir

(15)

4 空間目標(biāo)的紅外輻射特性成像

4.1 圖像灰度值計(jì)算

前一節(jié)計(jì)算出的輻射強(qiáng)度并不是紅外成像仿真需要的最終結(jié)果,需要把經(jīng)過(guò)成像系統(tǒng)效應(yīng)作用后的輻射強(qiáng)度轉(zhuǎn)化成灰度值。

現(xiàn)在通常使用的是均勻量化,輻射強(qiáng)度和灰度值之間是一個(gè)線(xiàn)性關(guān)系。一幅圖像中像面微元對(duì)應(yīng)的最高和最低輻射強(qiáng)度分別為Imin和Imax,設(shè)圖像灰度等級(jí)為0～255級(jí),則對(duì)應(yīng)于輻射強(qiáng)度為I(j)的像面微元的灰度值為:

(16)

4.2 基于改進(jìn)Delaunay三角剖分的紅外輻射特性成像方法

前面已經(jīng)計(jì)算出紅外圖像中節(jié)點(diǎn)的灰度值,下一步是將這些離散的節(jié)點(diǎn)導(dǎo)出,重建為空間目標(biāo),并且選取探測(cè)器所在角度進(jìn)行成像,這一步要采用Delaunay三角剖分[11]方法。

Delaunay三角剖分算法常用的分為兩種:基于投影的三角剖分算法[12],此算法容易出現(xiàn)拼接錯(cuò)誤并且不利于復(fù)雜曲面的重構(gòu)；另一種為基于Delaunay的生長(zhǎng)算法,即在三維離散點(diǎn)中構(gòu)造一個(gè)種子三角形,隨后根據(jù)Delaunay準(zhǔn)則尋找最優(yōu)點(diǎn)擴(kuò)展形成一個(gè)新的三角形,以此三角形為基礎(chǔ)再次向外搜索最優(yōu)點(diǎn),直到所有的散射點(diǎn)都加入到三角網(wǎng)中構(gòu)成三角形[13]。此算法的優(yōu)點(diǎn)是思想簡(jiǎn)單,能對(duì)復(fù)雜三維曲面進(jìn)行重構(gòu),但是此方法在搜索滿(mǎn)足條件的第三點(diǎn)時(shí)計(jì)算復(fù)雜,耗費(fèi)大量時(shí)間,且易出現(xiàn)曲面重疊、空洞等現(xiàn)象[14]。所以本文將采用三維柵格法[15]和自適應(yīng)空間外接球方法[16]對(duì)其進(jìn)行改進(jìn),提出一種改進(jìn)Delaunay三角剖分算法。

三維柵格法首先讀取三維點(diǎn)云坐標(biāo),求得點(diǎn)云在三個(gè)坐標(biāo)軸X、Y和Z軸上的最大值和最小值,確定能將所有數(shù)據(jù)包含在內(nèi)的最小長(zhǎng)方體,然后選區(qū)合適的邊長(zhǎng)值,將長(zhǎng)方體分割成若干個(gè)相同大小的立方體。最后依據(jù)各點(diǎn)的三位坐標(biāo)值將點(diǎn)存儲(chǔ)在相對(duì)應(yīng)的小立方體編號(hào)數(shù)組中,完成了點(diǎn)云數(shù)據(jù)的空間劃分。該方法可以根據(jù)點(diǎn)云的疏密程度自由調(diào)整子區(qū)域的大小,在領(lǐng)域點(diǎn)搜索的過(guò)程中,可以減小搜索范圍,提高效率。

自適應(yīng)空間外接球方法的思想是找到點(diǎn)集中距離生長(zhǎng)邊的中心點(diǎn)最近的點(diǎn)組成新三角形,確定此三角形外接球并找到球內(nèi)部的點(diǎn)云數(shù)據(jù),最后通過(guò)迭代找到最優(yōu)擴(kuò)展點(diǎn)。

綜上,本文的改進(jìn)Delaunay三角剖分算法的步驟為:

(1)利用上述三維柵格法完成某個(gè)點(diǎn)的k鄰域?qū)ふ遥?/p>

(2)構(gòu)造種子三角形,利用種子三角形來(lái)生長(zhǎng)三角網(wǎng)；

(3)采用自適應(yīng)空間外接球方法尋找三角形生長(zhǎng)的最優(yōu)點(diǎn)；

(4)對(duì)未使用的點(diǎn)集進(jìn)行迭代生長(zhǎng)。

本文采用這種方法先對(duì)空間目標(biāo)進(jìn)行重建,將紅外輻射強(qiáng)度轉(zhuǎn)化的灰度值分布于空間目標(biāo)的表面,隨后采用STK計(jì)算出探測(cè)器方位,根據(jù)其方位矢量得出二維圖像。

5 空間目標(biāo)的紅外輻射成像仿真實(shí)驗(yàn)

本文以某衛(wèi)星為例,進(jìn)行紅外輻射成像仿真實(shí)驗(yàn),首先進(jìn)行目標(biāo)三維建模并計(jì)算衛(wèi)星表面溫度分布,隨后利用第二節(jié)中的公式求出目標(biāo)紅外輻射特性大小,最后采用改進(jìn)的Delaunay三角剖分成像方法得到目標(biāo)的紅外圖像。

5.1 實(shí)驗(yàn)參數(shù)

首先表1中列出了探測(cè)目標(biāo)的物性參數(shù)。表2列出了目標(biāo)軌道與探測(cè)軌道參數(shù)。

表1 目標(biāo)物性參數(shù)表

表2 目標(biāo)與探測(cè)軌道參數(shù)

目標(biāo)和探測(cè)器軌道均為太陽(yáng)同步軌道,其中ad為軌道半長(zhǎng)軸;ed為軌道偏心率;id為軌道傾角;Ωd為升交點(diǎn)赤經(jīng);ωd為近地點(diǎn)輻角;τd為升交點(diǎn)地方時(shí)。

5.2 結(jié)果及分析

為了驗(yàn)證本文計(jì)算空間目標(biāo)紅外輻射強(qiáng)度的模型與成像方法是否合理有效。圖1、圖2首先比較了在3～5μm與8～12μm兩個(gè)波段下目標(biāo)本體面元與太陽(yáng)帆板正面面元總紅外輻射強(qiáng)度變化,從2009 年 6 月 5 日 15∶30開(kāi)始,橫坐標(biāo)為距開(kāi)始時(shí)間的秒數(shù),縱坐標(biāo)為紅外輻射強(qiáng)度大小。由于紅外輻射強(qiáng)度隨目標(biāo)周期性運(yùn)動(dòng)也呈現(xiàn)周期性變化,圖1展現(xiàn)了紅外輻射強(qiáng)度的兩個(gè)周期變化曲線(xiàn)。

圖1 兩種波段下目標(biāo)本體面元總輻射強(qiáng)度隨時(shí)間變化曲線(xiàn)

由圖1(a)與圖1(b)對(duì)比可知,8～12 μm波段探測(cè)的目標(biāo)本體面元輻射強(qiáng)度比3～5 μm更強(qiáng),采用8～12 μm波段更易于探測(cè)目標(biāo)。在某一時(shí)間段,目標(biāo)本體面元輻射強(qiáng)度達(dá)到最小保持不變,這時(shí)面元處于陰影區(qū)。

圖2 兩種波段下目標(biāo)太陽(yáng)帆板面元總輻射強(qiáng)度隨時(shí)間變化曲線(xiàn)

對(duì)比圖1、圖2,可以看出在光照區(qū),目標(biāo)太陽(yáng)帆板正面面元總輻射強(qiáng)度在8～12 μm波段大于目標(biāo)本體面元的總輻射強(qiáng)度,在3～5 μm波段并不明顯,這是因?yàn)樘?yáng)輻射能量主要集中在中短波,目標(biāo)本體反射率較大,所以在3～5 μm波段反射了更多的太陽(yáng)能量。

圖3為目標(biāo)紅外輻射成像仿真結(jié)果,仿真從5400 s(17∶00)開(kāi)始,直到114000 s(18∶40)結(jié)束,每隔1200 s成一次像。

圖3 不同時(shí)間目標(biāo)雙波段成像結(jié)果

由圖3可以看出,同一時(shí)間8～12 μm波段輻射強(qiáng)度大于3～5 μm波段。目標(biāo)的成像角度隨時(shí)間不斷變化,在光照區(qū),目標(biāo)的太陽(yáng)帆板正面的亮度大于目標(biāo)本體,滿(mǎn)足理論實(shí)際。綜上所述,本文提出的基于改進(jìn)Delaunay三角剖分的空間目標(biāo)紅外輻射成像方法有效。

為了比較傳統(tǒng)的Delaunay三角剖分算法與本文改進(jìn)后的Delaunay三角剖分算法的重建效率,進(jìn)行了兩種算法對(duì)同一目標(biāo)的重建時(shí)間對(duì)比實(shí)驗(yàn),實(shí)驗(yàn)平臺(tái)均為Matlab R2017a,且均在安裝 64 位Windows 10 操作系統(tǒng)的同一臺(tái)計(jì)算機(jī)上完成,表3是運(yùn)行時(shí)間比較結(jié)果。

表3 算法運(yùn)行時(shí)間比較

由表3可以看出,本文改進(jìn)后的算法重建時(shí)間比傳統(tǒng)Delaunay三角剖分算法重建時(shí)間短,效率更高。

6 結(jié) 論

本文提出了一種基于改進(jìn)Delaunay三角剖分的空間目標(biāo)紅外輻射成像方法,建立了一套完整的空間目標(biāo)紅外輻射成像方法。此方法首先根據(jù)熱平衡方程,利用節(jié)點(diǎn)網(wǎng)絡(luò)法求解出空間目標(biāo)的各節(jié)點(diǎn)溫度,然后利用STK得出探測(cè)器、目標(biāo)、地球與太陽(yáng)的位置關(guān)系,隨后由目標(biāo)節(jié)點(diǎn)溫度求得目標(biāo)紅外輻射強(qiáng)度,最后利用改進(jìn)的Delaunay三角剖分的方法重建目標(biāo)得到紅外圖像。本文以某衛(wèi)星為例進(jìn)行了成像仿真實(shí)驗(yàn),并對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行了分析,驗(yàn)證了本文方法的有效性,此方法對(duì)空間目標(biāo)的紅外探測(cè)以及識(shí)別具有重要意義。