基于卷積曲面的動(dòng)態(tài)實(shí)時(shí)星圖模擬

閆勁云， 劉慧， 趙偉強(qiáng)， 江潔

(1. 中國(guó)計(jì)量科學(xué)研究院, 北京 100029； 2. 北京航空航天大學(xué)儀器科學(xué)與光電工程學(xué)院, 北京 100083)

星敏感器是航天器中最重要的姿態(tài)敏感設(shè)備，在設(shè)計(jì)、調(diào)試和測(cè)試等環(huán)節(jié)中，需要提供輸入激勵(lì)以完善和評(píng)估其各項(xiàng)功能和性能，星模擬器是實(shí)現(xiàn)這一任務(wù)的最佳途徑[1-2]。星模擬器包括單星模擬器、靜態(tài)多星模擬器、動(dòng)態(tài)多星模擬器[3-4]。動(dòng)態(tài)多星模擬器可模擬全天區(qū)任意指向的星圖，然后由光學(xué)系統(tǒng)將星圖輸出為適合星敏感器接收的星光信號(hào)[5]。

隨著星敏感器的動(dòng)態(tài)性能和數(shù)據(jù)更新率等各種指標(biāo)的不斷提高，對(duì)多星模擬器的要求也隨之提高[6-7]。星敏感器的動(dòng)態(tài)性能是指當(dāng)其隨著載體運(yùn)動(dòng)時(shí)，星點(diǎn)成像發(fā)生拖尾的情況下，保持正常姿態(tài)輸出的能力。為了測(cè)試動(dòng)態(tài)性能，要求多星模擬器可以模擬拖尾的星點(diǎn)，目前多星模擬器對(duì)于動(dòng)態(tài)拖尾星點(diǎn)模擬的算法研究還較少[5,8]。當(dāng)前星敏感器的數(shù)據(jù)更新率普遍可達(dá)10 Hz，一些星敏感器可達(dá)30 Hz[6]，多星模擬器要求動(dòng)態(tài)刷新顯示變化的星空，其刷新率不低于星敏感器的數(shù)據(jù)更新率，通常大于30 Hz，因此對(duì)算法的實(shí)時(shí)性也提出了較高要求[9]。

多星模擬器算法的核心任務(wù)是高精度實(shí)時(shí)星圖生成，主要包括導(dǎo)航星搜索、星點(diǎn)坐標(biāo)變換、星點(diǎn)光斑彌散及像素離散化等步驟。導(dǎo)航星搜索是指根據(jù)星敏感器的視軸指向在全天區(qū)范圍內(nèi)搜索當(dāng)前視場(chǎng)內(nèi)的星點(diǎn)信息。星點(diǎn)坐標(biāo)變換將搜索到的星點(diǎn)信息變換為模擬星圖中的圖像坐標(biāo)和灰度信息。星點(diǎn)光斑彌散模擬了鏡頭的彌散作用，將星點(diǎn)變換為星點(diǎn)光斑[10]。像素離散化過(guò)程將連續(xù)的星點(diǎn)光斑變?yōu)殡x散的數(shù)字圖像。通過(guò)星表分區(qū)[3,5]和重新組織構(gòu)建星表[11]等方法，可以解決搜索導(dǎo)航星耗時(shí)較高的問(wèn)題，本文對(duì)此不再贅述。靜態(tài)下的星點(diǎn)坐標(biāo)變換模型較為成熟，但在動(dòng)態(tài)條件下，星點(diǎn)在模擬星圖中的坐標(biāo)構(gòu)成一條軌跡。劉海波[8]和劉朝山[12]等推導(dǎo)了星點(diǎn)在像面的速度與星敏感器角速度的關(guān)系，但未給出星點(diǎn)軌跡曲線的表達(dá)式。Sun等[7]在推導(dǎo)過(guò)程中使用了一階近似，因此僅適用于曝光時(shí)間較短的情況。本文在此基礎(chǔ)上，給出了完備的星點(diǎn)軌跡模型。對(duì)于模擬星點(diǎn)光斑彌散過(guò)程，Sun等[7]利用運(yùn)動(dòng)模糊圖像中的概念，將動(dòng)態(tài)拖尾星點(diǎn)在頻域中表示為鏡頭彌散函數(shù)和運(yùn)動(dòng)退化函數(shù)的乘積，但運(yùn)動(dòng)退化函數(shù)在時(shí)域中的物理含義仍有待進(jìn)一步明確。Yan等[6]利用線擴(kuò)散函數(shù)的概念，得出了拖尾星點(diǎn)光斑的解析表達(dá)式模型，但僅適用于星點(diǎn)在像面做勻速直線運(yùn)動(dòng)的情況。本文針對(duì)這些問(wèn)題，建立了拖尾星點(diǎn)的卷積曲面模型，更準(zhǔn)確地描述了動(dòng)態(tài)拖尾星點(diǎn)的成像過(guò)程。對(duì)于像素離散化過(guò)程，由于涉及對(duì)光斑曲面的二重積分，普遍存在計(jì)算復(fù)雜、實(shí)時(shí)性較差的問(wèn)題[13]。Yan等[6]的方法僅適用于勻速直線運(yùn)動(dòng)的情況，劉朝山等[12]的方法假定星點(diǎn)掃過(guò)每個(gè)像素的時(shí)間一致，因此不適合非勻速運(yùn)動(dòng)的情況。本文在卷積曲面模型的基礎(chǔ)上提出了一種快速像素離散算法，適用于星點(diǎn)做任意運(yùn)動(dòng)的情況，并且大大加快仿真速度，以滿足動(dòng)態(tài)多星模擬器的實(shí)時(shí)星圖仿真要求。

1 動(dòng)態(tài)拖尾星點(diǎn)模擬

1.1 星點(diǎn)運(yùn)動(dòng)軌跡

當(dāng)星敏感器隨載體發(fā)生角位移時(shí)，在曝光時(shí)間內(nèi)，星點(diǎn)在星敏感器像面移動(dòng)成像，星點(diǎn)中心的運(yùn)動(dòng)軌跡為一條曲線段。

設(shè)星表中恒星i在天球坐標(biāo)系下的映像為Vi。在t時(shí)刻，此恒星在星敏感器坐標(biāo)系下的映像為Wi(t)，則Vi和Wi(t)可分別表示為

(1)

(2)

式中：αi和δi分別為赤經(jīng)和赤緯；(xi(t),yi(t))為星點(diǎn)在星敏感器像面的投影坐標(biāo)；F為星敏感器的鏡頭焦距。

Wi(t)和Vi是同一個(gè)矢量在不同坐標(biāo)系下的映像，二者關(guān)系為

Wi(t)=R(t)Vi

(3)

式中：R(t)為星敏感器的姿態(tài)矩陣。

設(shè)星敏感器相對(duì)于天球坐標(biāo)系的旋轉(zhuǎn)角速度為ω=[ωx,ωy,ωz]T，式(3)對(duì)時(shí)間t求導(dǎo)得[8]

-[ω×]Wi(t)

(4)

(5)

若曝光起始時(shí)刻t0星點(diǎn)在像面坐標(biāo)為(xi0,yi0)，則動(dòng)態(tài)下星點(diǎn)在像面的運(yùn)動(dòng)軌跡C為

(6)

1.2 基于卷積曲面的光斑彌散模型

星敏感器為了獲得亞像素級(jí)的定位精度，無(wú)窮遠(yuǎn)處的星光經(jīng)過(guò)星敏感器鏡頭在圖像傳感器像面并非匯聚于單個(gè)像素，而是成像為一個(gè)彌散光斑。靜態(tài)下星點(diǎn)光斑能量分布可表示為[6-7]

f(x,y)=ΦTh(x-xc,y-yc)

(7)

式中：Ф為星光被鏡頭收集并入射到像面的輻射通量；T為曝光時(shí)間；(xc,yc)表示星點(diǎn)光斑中心位置；h(x,y)為鏡頭的點(diǎn)擴(kuò)散函數(shù)，通常以二維高斯函數(shù)來(lái)描述[6]。

(8)

(9)

式中：E0為0等星在地球大氣外表面處的輻照度[1]，E0= 2.96×10-14W/mm2；m為入射星光的星等；D為星敏感器的光學(xué)鏡頭孔徑；τ為鏡頭的透射率；ρ為鏡頭的彌散半徑。

動(dòng)態(tài)情況下，星點(diǎn)在像面發(fā)生移動(dòng)，星點(diǎn)光斑的能量分布構(gòu)成一個(gè)二維卷積曲面[14-15]，如圖1所示。圖1中，星點(diǎn)的軌跡曲線為C，r為軌跡曲線 上的點(diǎn)，r=(xc(t),yc(t))，其參數(shù)方程見(jiàn)式(6)。

圖1 拖尾星點(diǎn)光斑成像過(guò)程Fig.1 Imaging process of smeared star spot

p為像平面內(nèi)的任意一點(diǎn)，p=(x,y)。h(p)為星敏感器光學(xué)成像系統(tǒng)的點(diǎn)擴(kuò)散函數(shù)，h(|p-r|)表示此函數(shù)的中心位于點(diǎn)r。記權(quán)重函數(shù)g(r)表示星點(diǎn)經(jīng)過(guò)軌跡曲線上的弧微元|dr|所需的時(shí)間，則有

(10)

在點(diǎn)r附近取弧微元|dr|，星點(diǎn)掃過(guò)此微元期間的能量為Φg(r)|dr|。顯然，當(dāng)星點(diǎn)的速度越快，星光累積的能量越少；反之，當(dāng)星點(diǎn)的速度越慢，能量累積越多。由于光學(xué)鏡頭的彌散作用，此星點(diǎn)能量按照點(diǎn)擴(kuò)散函數(shù)描述的方式彌散到了像面，能量分布可表示為Φg(r)h(|p-r|)|dr|，那么沿曲線積分后可得星點(diǎn)在像面的總能量分布為

Φg(p)?h(p)

(11)

式中：?為卷積運(yùn)算符。星點(diǎn)能量分布為一卷積曲面，即f(p)為點(diǎn)擴(kuò)散函數(shù)h(p)與權(quán)重函數(shù)g(r)的卷積。拖尾星點(diǎn)光斑的卷積曲面模型具有清晰的物理含義，但不便于運(yùn)算，為此進(jìn)一步導(dǎo)出動(dòng)態(tài)下星點(diǎn)能量分布積分模型為

(12)

式(11)和式(12)所表示的模型精確描述了動(dòng)態(tài)條件下的拖尾星點(diǎn)成像過(guò)程。

1.3 卷積曲面的快速像素離散算法

入射到像面的光斑能量，在像元內(nèi)轉(zhuǎn)化為電荷信號(hào)，再經(jīng)過(guò)讀出和放大，轉(zhuǎn)換為數(shù)字圖像。那么位于第i行和第j列的像素的圖像灰度值為[1,5-8,12-13]

(13)

式中：K為轉(zhuǎn)換系數(shù)，是圖像傳感器的量子效率和轉(zhuǎn)換增益的乘積。圖像傳感器對(duì)光斑能量的作用不僅僅體現(xiàn)為離散采樣，還包括積分作用，此過(guò)程稱為像素離散。顯然，像素灰度值為一個(gè)三重積分表達(dá)式，因此仿真一幅星圖需要對(duì)每個(gè)星點(diǎn)光斑范圍內(nèi)的每一個(gè)像素計(jì)算三重積分，運(yùn)算的時(shí)間開(kāi)銷較大，不能滿足動(dòng)態(tài)多星模擬器實(shí)時(shí)性的要求。

由拖尾星點(diǎn)的卷積曲面模型的物理含義出發(fā)，提出一種快速算法，以解決運(yùn)算量較大的問(wèn)題。具體思路為：將曝光時(shí)間T等間隔地分為N段，每段時(shí)間間隔ΔT=T/N。相應(yīng)的，軌跡曲線也被分為了N段。當(dāng)ΔT足夠短時(shí)，地軌跡上的每一段可近似為一點(diǎn)，此時(shí)的星點(diǎn)光斑成像近似為靜態(tài)下的情況。靜態(tài)下對(duì)二維高斯函數(shù)的二重積分可推導(dǎo)得其原函數(shù)，將耗時(shí)的數(shù)值積分轉(zhuǎn)化為函數(shù)求值。最后將這N點(diǎn)沿軌跡累加，即可得到

(14)

(15)

式中：(xck,yck)為第k時(shí)間段時(shí)星點(diǎn)光斑中心在軌跡上的坐標(biāo)；erf(x)為誤差函數(shù)。式(15)將復(fù)雜的三重積分化簡(jiǎn)為函數(shù)求值、乘法和加法運(yùn)算，大大提高了拖尾星點(diǎn)的仿真速度。

2 仿真及分析

仿真計(jì)算機(jī)的配置如下。CPU：Inter Core i7，內(nèi)存8 GB，仿真環(huán)境：MATLAB 2015b。仿真星敏感器的主要參數(shù)為：視場(chǎng)20°×20°，鏡頭直徑27 mm，焦距37.92 mm，像素尺寸5.5 μm × 5.5 μm，圖像分辨率2 048像素×2 048像素。采用SAO星表，選取亮于6.0等星的星作為仿真星表。

2.1 精度對(duì)比

本節(jié)首先對(duì)比3種星點(diǎn)仿真算法的精度，再研究對(duì)軌跡劃分到何種精度，本文算法的精度可以達(dá)到滿意結(jié)果。式(13)表示的三重積分法是未優(yōu)化的原始算法，速度慢，但精度高，得到的星點(diǎn)光斑可認(rèn)為是真值。本文算法的精度依賴于對(duì)積分劃分的細(xì)度，隨著N增大，精度逐漸逼近真值。第3種算法為運(yùn)動(dòng)模糊圖像卷積法，拖尾星點(diǎn)由原始星點(diǎn)圖像與運(yùn)動(dòng)模糊核函數(shù)圖像卷積而得[7,10]。定義星點(diǎn)圖像灰度誤差為

ferr=|E(j,i)-I(j,i)|

(16)

進(jìn)一步定義最大灰度誤差為星點(diǎn)區(qū)域內(nèi)誤差絕對(duì)值最大的灰度值；定義灰度RMS誤差為星點(diǎn)區(qū)域內(nèi)每一個(gè)像素的誤差灰度值的均方差。

設(shè)[ωx,ωy,ωz]=[5(°)/s, 5(°)/s,0]，曝光時(shí)間T=20 ms，鏡頭的彌散半徑為0.7像素，星點(diǎn)初始位置為像面中心，則由式(6)可得，星點(diǎn)拖尾長(zhǎng)度為14.33像素。采用三重積分法得到的拖尾星點(diǎn)如圖2(a)所示；采用本文算法，將星點(diǎn)軌跡劃分為143段，仿真得到的拖尾星點(diǎn)如圖2(b)所示；采用圖像卷積法得到拖尾星點(diǎn)如圖2(c)所示；后兩者的星點(diǎn)圖像灰度誤差分別如圖2(d)、(e)所示。為便于比較，圖中星點(diǎn)灰度值在最大值處進(jìn)行了歸一化。

對(duì)比圖2(a)、(b)可知，本文仿真結(jié)果與真值較為一致。由圖2(d)可知，僅起始點(diǎn)和終點(diǎn)附近存在一定誤差，星點(diǎn)拖尾中間部位的誤差可忽略，其灰度RMS誤差為0.005 4。由圖2(a)、(c)可知，采用圖像卷積法得到的拖尾星點(diǎn)，其星點(diǎn)能量彌散程度較大，圖2(e)所示的星點(diǎn)誤差也較大，其灰度RMS誤差為0.102 0。

由式(14)可知，當(dāng)軌跡劃分越細(xì)，則仿真精度越高。本文算法在不同劃分間隔下，星點(diǎn)灰度誤差如表1所示。

由表1可知，隨著軌跡劃分間隔增加，星點(diǎn)灰度誤差逐漸增大。由于成像過(guò)程中存在的背景噪聲及其他噪聲，星點(diǎn)圖像灰度存在起伏，若仿真星點(diǎn)最大誤差小于噪聲引起的灰度起伏，則仿真誤差對(duì)最終結(jié)果可忽略。若以最大灰度誤差為限，軌跡劃分間隔為0.05像素時(shí)，仿真星點(diǎn)精度滿足要求；若以RMS誤差為限，軌跡劃分為0.10像素即可滿足要求。

圖2 拖尾星點(diǎn)光斑仿真精度對(duì)比Fig.2 Simulation accuracy comparison of smeared star spots

劃分間隔/像素0.020.050.100.20最大灰度誤差0.00530.01290.02600.0525灰度RMS誤差0.00110.00270.00540.0109

2.2 仿真速度對(duì)比

令多星模擬器產(chǎn)生1 000個(gè)隨機(jī)視軸指向，分別利用三重積分法、運(yùn)動(dòng)模糊圖像卷積法和本文算法仿真生成拖尾星圖，統(tǒng)計(jì)每幅仿真星圖的仿真時(shí)間，3種算法仿真耗時(shí)的直方圖統(tǒng)計(jì)如圖3所示。

由圖3可知，三重積分法仿真一幅星圖所需時(shí)間為幾秒，運(yùn)動(dòng)模糊圖像卷積法所需時(shí)間為幾百毫秒，本文算法耗時(shí)為幾十毫秒，3種算法仿真一幅星圖的平均耗時(shí)分別為3 030、309和22 ms。顯然，本文算法的仿真速度大大優(yōu)于前兩種算法。

圖3(c)所示的個(gè)別星圖的仿真耗時(shí)大于33 ms，不滿足多星模擬器刷新率的要求。為此，繪制模擬星圖仿真耗時(shí)與星點(diǎn)數(shù)量的關(guān)系圖(見(jiàn)圖4)，以進(jìn)一步分析超時(shí)原因。

由圖4可知，星圖仿真耗時(shí)與星點(diǎn)數(shù)量是線性關(guān)系。星點(diǎn)數(shù)量較多時(shí)，仿真耗時(shí)也越多。對(duì)于單個(gè)拖尾星點(diǎn)，平均仿真耗時(shí)為0.45 ms，那么為保證單幅星圖仿真時(shí)間小于33 ms，只要星圖中 星點(diǎn)個(gè)數(shù)小于70即可。為此，對(duì)于模擬星圖中星點(diǎn)數(shù)量大于70顆星的情況，根據(jù)星點(diǎn)的亮度信息，篩選其中較亮的70顆星，其余亮度較暗的星予以剔除。這樣的處理對(duì)實(shí)際使用影響較小，原因在于星敏感器為獲得較高的姿態(tài)精度，僅采用較亮的若干顆星進(jìn)行姿態(tài)解算。較暗的星形成的星點(diǎn)光斑圖像的信噪比低，因此對(duì)其定位的精度也較低，應(yīng)予以剔除。增加對(duì)多余暗星剔除的步驟，得到的仿真耗時(shí)直方圖如圖5所示。

圖3 不同算法仿真耗時(shí)直方圖Fig.3 Histogram of simulation time consumption for different algorithm

圖4 星圖仿真耗時(shí)與星圖中星點(diǎn)數(shù)量的關(guān)系Fig.4 Star map simulation time consumption vs. star number of star image

由圖5可知，剔除多余暗星后，本文的快速星點(diǎn)仿真算法可滿足動(dòng)態(tài)多星模擬器的實(shí)時(shí)星圖仿真的要求。

圖5 剔除暗星后仿真耗時(shí)直方圖Fig. 5 Histogram of simulation time consumption after dark star is excluded

3 結(jié) 論

針對(duì)動(dòng)態(tài)多星模擬器對(duì)拖尾星點(diǎn)仿真難以滿足實(shí)時(shí)性的問(wèn)題，本文首先建立了動(dòng)態(tài)拖尾星點(diǎn)的卷積曲面成像模型，然后提出了一種基于卷積曲面的快速像素離散算法。得到的仿真星點(diǎn)誤差低于圖像噪聲水平，仿真精度高；同時(shí)仿真一幅星圖耗時(shí)小于33 ms，比現(xiàn)有算法速度提高14倍，滿足了多星模擬器刷新率30 Hz的要求。