楊 飛，金 光 ，謝金華，邱振戈，曲宏松，賀小軍

(1.中國科學(xué)院長春光學(xué)精密機(jī)械與物理研究所，吉林長春130033;2.小衛(wèi)星技術(shù)國家地方聯(lián)合工程研究中心，吉林長春130033;3.中國科學(xué)院大學(xué)，北京100049;4.國家測繪地理信息局衛(wèi)星測繪應(yīng)用中心，北京101300;5.上海海洋大學(xué)，上海201306)

1 引言

由于航天任務(wù)具有高風(fēng)險(xiǎn)、高投入、研制周期長等特點(diǎn)，因此在地面實(shí)驗(yàn)環(huán)境下進(jìn)行光學(xué)成像建模與仿真是完成載荷和平臺大量測試工作的重要途徑。高精度、高分辨率立體測繪相機(jī)成像仿真在衛(wèi)星技術(shù)指標(biāo)論證、試驗(yàn)測試、在軌運(yùn)行評價(jià)以及故障模擬分析中都發(fā)揮著十分重要的作用。精確的成像仿真可以完善衛(wèi)星平臺及載荷指標(biāo)系統(tǒng)，使得在研衛(wèi)星性能達(dá)到最優(yōu)，從而保證衛(wèi)星在軌成像質(zhì)量。因此，基于在軌成像物理機(jī)理的全鏈路仿真在衛(wèi)星研制過程中占據(jù)著舉足輕重的地位。

對于可見光成像系統(tǒng)的建模與仿真，國內(nèi)外專家做了大量的工作，主要有物理仿真、半物理仿真和數(shù)學(xué)仿真3種技術(shù)途徑，其中數(shù)學(xué)仿真因其工程成本低而成為了研究的重點(diǎn)。1999年，德國宇航中心開發(fā)了光學(xué)遙感仿真軟件SENSOR［1-2］，首次實(shí)現(xiàn)了全鏈路成像的模擬流程，該軟件針對高光譜成像儀，通過光線追蹤算法完成了太陽、地物和傳感器之間的幾何建模，構(gòu)建某一成像時刻像元與觀測地物對應(yīng)的幾何關(guān)系，然后通過調(diào)用查找表計(jì)算到達(dá)傳感器入瞳處的輻亮度，最后經(jīng)過光電模型得到數(shù)字圖像。2010年，Cota［3］等人根據(jù)不同材質(zhì)的反射率將高質(zhì)量輸入場景轉(zhuǎn)化為反射系數(shù)，結(jié)合MODTRAN得到入瞳輻亮度，通過光學(xué)系統(tǒng)傳遞函數(shù)、CCD傳感器傳遞函數(shù)和信噪比模型建立了航天成像系統(tǒng)端到端仿真工具PICASSO。2014年，劉曉等人［4］針對高分辨率星載光學(xué)遙感器成像特點(diǎn)，提出一種基于低空遙感系統(tǒng)的成像仿真方法。2002年，王剛等人［5］通過分析成像鏈的光譜輻射響應(yīng)和空間響應(yīng)特性，提出了一種基于圖像仿真的對地遙感過程科學(xué)可視化方法。2007年，陳方等人［6］依靠波譜數(shù)據(jù)庫的數(shù)據(jù)支持，提出了利用寬光譜光學(xué)遙感圖像模擬細(xì)分光譜光學(xué)遙感圖像的模擬技術(shù)。

以帶波譜數(shù)據(jù)的高精度、高分辨率地表物理模型為輸入源，本文提出了一種基于在軌成像物理機(jī)理的立體測繪相機(jī)全鏈路數(shù)值仿真方法。高精度高分辨率地表物理模型由通過唯一ID關(guān)聯(lián)的三角網(wǎng)DEM和波譜屬性組成，本方法由建立的姿態(tài)軌道模型和探測器安裝位置生成需要跟蹤的光線，結(jié)合MODTRAN軟件模擬計(jì)算大氣輻射傳輸，求得相機(jī)入瞳處輻亮度［7-10］，進(jìn)而采用蒙特卡洛光線追跡算法仿真光線進(jìn)入相機(jī)后的變化，最后在相機(jī)輻射響應(yīng)模型中仿真CCD對光線的響應(yīng)，得到數(shù)字影像。輸入源為全鏈路成像仿真提供近似真實(shí)的影像，支撐高精度成像仿真技術(shù)實(shí)現(xiàn)，各模塊的精確性保證了本文全鏈路成像仿真技術(shù)的可行性。

2 立體測繪相機(jī)成像鏈模型

電磁輻射經(jīng)過大氣傳輸?shù)竭_(dá)地球表面，地表景物反射的電磁波再經(jīng)過大氣輻射傳輸進(jìn)入光學(xué)系統(tǒng)，光學(xué)系統(tǒng)對一定波段范圍的電磁波能量聚焦成像于探測器靶面［11］。遙感探測器大多采用TDI CCD固體探測器，利用光電轉(zhuǎn)換作用將到達(dá)靶面的光信號轉(zhuǎn)換為電信號，在信號處理電路中經(jīng)過放大處理和模數(shù)轉(zhuǎn)換后儲存。獲得的原始影像數(shù)據(jù)需經(jīng)過相應(yīng)圖像處理后用于圖像顯示和分析等。這些環(huán)節(jié)構(gòu)成了一個完整的物理成像模型。

基于在軌成像物理機(jī)理對立體測繪相機(jī)進(jìn)行建模與仿真時，需要分析平臺姿態(tài)軌道、光學(xué)系統(tǒng)、相機(jī)輻射響應(yīng)等主要成像鏈環(huán)節(jié)。

2.1 姿態(tài)軌道建模

平臺姿態(tài)及軌道動力學(xué)仿真模型是整個動態(tài)成像系統(tǒng)中的基本組成部分。軌道動力學(xué)模型主要模擬衛(wèi)星軌道運(yùn)動，計(jì)算軌道信息;姿態(tài)動力學(xué)模型主要模擬衛(wèi)星姿態(tài)運(yùn)動，計(jì)算姿態(tài)信息。該仿真模型基本計(jì)算過程如下:

(1)給定初始時刻衛(wèi)星平臺在J2000.0慣性坐標(biāo)系內(nèi)位置與速度;

(2)基于J2000.0慣性坐標(biāo)系內(nèi)的二體軌道動力學(xué)微分方程，利用數(shù)值積分算法，對二體軌道動力學(xué)微分方程積分，得到下一時刻衛(wèi)星在慣性坐標(biāo)系下的位置矢量和速度矢量;

(3)由衛(wèi)星在J2000.0慣性坐標(biāo)系下的位置速度矢量，計(jì)算相應(yīng)的軌道根數(shù)(半長軸、偏心率、軌道傾角、近心點(diǎn)角距、升交點(diǎn)赤經(jīng)、偏近點(diǎn)角及平近點(diǎn)角);

(4)給定初始時刻衛(wèi)星本體坐標(biāo)系相對于J2000.0慣性坐標(biāo)系的姿態(tài)四元數(shù)，下一時刻的姿態(tài)角信息由三次多項(xiàng)式得到;

(5)根據(jù)所建立的 J2000.0慣性坐標(biāo)系至WGS84坐標(biāo)系時的坐標(biāo)轉(zhuǎn)換矩陣模型，結(jié)合載荷在衛(wèi)星平臺的安裝矩陣及地面目標(biāo)點(diǎn)在WGS84系下的經(jīng)緯度坐標(biāo)，帶入下一時刻的軌道位置參數(shù)及衛(wèi)星平臺的姿態(tài)參數(shù)，可計(jì)算得到成像仿真所需的前視、正視相機(jī)與目標(biāo)區(qū)域間的相對幾何參數(shù)，為基于光線追跡的成像仿真模型提供衛(wèi)星姿態(tài)及軌道輸入?yún)?shù);

(6)重復(fù)上述過程，即可得到不同成像時刻的衛(wèi)星姿態(tài)及軌道仿真參數(shù)。

通過上述方法得到衛(wèi)星的外方位元素后，根據(jù)文獻(xiàn)［12］中提到的坐標(biāo)轉(zhuǎn)換方法，計(jì)算正視相機(jī)和前視相機(jī)像面上各點(diǎn)的光線在WGS84坐標(biāo)系下的矢量方向。

2.2 光學(xué)系統(tǒng)建模

用于立體測繪相機(jī)成像仿真的輸入場景是由一系列包含了幾何信息和輻射信息的三角面片組成。三角網(wǎng)數(shù)據(jù)包括三角面片ID標(biāo)識號、3個頂點(diǎn)的三維坐標(biāo)及輻射信息。光線求交算法計(jì)算2.1中求取的觀測光線矢量方向上三角面片與光線的交點(diǎn)，索引距離最近的三角面片幾何及輻射信息。

入瞳輻亮度可表示為［11］:

式中，Etop、Ed分別為總的上行輻照度和下行輻照度，ρ為目標(biāo)光譜反射率，τ為大氣透過率，σ為太陽天頂角，Lu為上行輻亮度。

式中，k為玻爾茲曼數(shù)，T為黑體的絕對溫度，h為普朗克常數(shù)。

式中，βsca為大氣散射系數(shù)。

根據(jù)CCD上光線出發(fā)點(diǎn)的密度將光學(xué)系統(tǒng)PSF離散化，取與CCD平均間隔相同，設(shè)PSF有效寬度為w，CCD光敏面尺寸為d×d，光線密度為m×m，則PSF的離散間隔為:

計(jì)算綜合PSF，以一維的點(diǎn)擴(kuò)散函數(shù)為例，綜合PSF計(jì)算方法為:

式中，w為點(diǎn)擴(kuò)散函數(shù)的有效范圍，pixsize為像元尺寸。

本文以離軸反射式光學(xué)系統(tǒng)為例，設(shè)計(jì)了一個焦距為10 m，F(xiàn)#為8.5的離軸反射式系統(tǒng)，系統(tǒng)的口徑為1 180 mm，體積為4 m×3 m×1.5 m，由四塊反射鏡形成，其中三塊非球面鏡，一塊為平面鏡。光學(xué)系統(tǒng)光路示意圖如圖1所示。

圖1 光學(xué)系統(tǒng)光路示意圖Fig.1 Light path diagram of optical system

2.3 相機(jī)輻射響應(yīng)建模

相機(jī)動態(tài)成像模擬實(shí)際是將連續(xù)物理的積分過程以足夠的精度離散化，用數(shù)值積分的方法計(jì)算空間TDI CCD相機(jī)動態(tài)成像的過程。高精度相機(jī)輻射響應(yīng)建模主要包括三部分內(nèi)容:TDI CCD延時積分特性的電子學(xué)采樣模型、焦平面光譜響應(yīng)幅值分布特性模型和焦平面光譜響應(yīng)SNR模型。

2.3.1 TDI CCD延時積分特性的電子學(xué)采樣模型

CCD相機(jī)在軌動態(tài)積分成像模擬模塊根據(jù)TDI CCD相機(jī)成像的工作原理，負(fù)責(zé)將輸入的目標(biāo)相機(jī)入瞳輻亮度場轉(zhuǎn)換為入瞳輻亮度影像。在動態(tài)積分成像模擬中要考慮量子效率、轉(zhuǎn)換效率、信噪比、MTF等，不考慮CCD器件工作的驅(qū)動時序和控制時序。TDI CCD成像由單個CCD像元曝光和多個CCD像元積分累加兩個步驟完成。

電壓輸出的響應(yīng)度R為:

式中:G為CCD片內(nèi)放大倍數(shù);q為電子電荷;η為量子效率;Ap為像元有效面積 ;K為TDI CCD電荷包總轉(zhuǎn)移效率;De單位時間內(nèi)暗信號電子數(shù)。

光譜范圍為λ1～λ2的CCD輸出電壓Vs為:

式中，R(λ)由式(8)得到;E(λ)為像面輻照度;tint為積分時間。

模擬電壓信號VS經(jīng)過模數(shù)轉(zhuǎn)換電路進(jìn)行灰度量化，假設(shè)量化位數(shù)為 N，最大量化電壓為Vmax，最小量化電壓為Vmin，則得到系統(tǒng)的輸出圖像的灰度滿足如下關(guān)系:

2.3.2 焦平面光譜響應(yīng)幅值分布特性模型

光譜響應(yīng)是指芯片對于不同波長光線的響應(yīng)能力，通常用光譜響應(yīng)曲線給出。通過光譜響應(yīng)曲線能夠直觀看出芯片對不同波長光線的響應(yīng)能力，與人眼相比，芯片的光譜響應(yīng)范圍要寬很多，對于波長小于截止波長的紅外、紫外和X-ray光子都能夠響應(yīng)。在選擇芯片時，要根據(jù)具體應(yīng)用的需求選擇光譜響應(yīng)合適的產(chǎn)品。同樣，一旦芯片確定后，其光譜響應(yīng)曲線也就隨之確定了，因此，在建模時，應(yīng)根據(jù)選擇芯片手冊提供的光譜響應(yīng)幅值分布曲線進(jìn)行相應(yīng)建模參數(shù)輸入。

圖2為SPRITE Sensor(IT-EB-4096)的三片TDI CCD的光譜響應(yīng)幅值與波長的關(guān)系曲線。每條響應(yīng)曲線都有一個峰值波長，峰值波長響應(yīng)度一半處對應(yīng)的波長稱為截止波長，只有波長小于截止波長的光線才能被芯片所感應(yīng)。

2.3.3 焦平面光譜響應(yīng)SNR模型

圖2 TDI CCD芯片典型光譜響應(yīng)曲線及有效量子效率圖Fig.2 Typical spectrum response curves and EQ chart

探測器成像過程中不可避免地引入噪聲，從而降低對感興趣目標(biāo)的辨識能力，從成像鏈中引入噪聲的來源看，噪聲主要分為霰粒噪聲Nshot、暗電流噪聲Ndark、復(fù)位噪聲NKTC、轉(zhuǎn)移噪聲NCTE等，噪聲沒有相關(guān)性，因此系統(tǒng)總噪聲可表示為:

3 仿真方法與分析

基于在軌成像物理機(jī)理的高精度立體測繪相機(jī)建模與仿真方法描述如下，如圖3所示。

圖3 全鏈路仿真流程圖Fig.3 Simulation flow chart of all link

(1)根據(jù)目標(biāo)相機(jī)成像時的大氣條件、太陽高度角、太陽方位角、衛(wèi)星高度角、衛(wèi)星方位角等參數(shù)，模擬大氣透過率、大氣后向散射、大氣鄰近效應(yīng)，利用MODTRAN軟件計(jì)算目標(biāo)相機(jī)觀測地面目標(biāo)在相機(jī)入瞳處的輻亮度。

(2)根據(jù)衛(wèi)星軌道高度、軌道傾角、偏心率、降交點(diǎn)時刻及成像目標(biāo)位置，綜合考慮星載相機(jī)的視場角范圍，計(jì)算能夠?qū)δ繕?biāo)區(qū)域中心成像的衛(wèi)星運(yùn)行時刻，輸出這一時刻前后數(shù)組衛(wèi)星位置、速度向量及對應(yīng)時間。

(3)根據(jù)衛(wèi)星姿態(tài)穩(wěn)定度、姿態(tài)指向精度、偏流角修正精度等設(shè)計(jì)參數(shù)及不同顫振分量的頻率、幅值的統(tǒng)計(jì)信息［13-14］，生成瞬時軌道坐標(biāo)系下的衛(wèi)星姿態(tài)角。

(4)根據(jù)目標(biāo)相機(jī)的探測器件安裝參數(shù)及相機(jī)的內(nèi)方位參數(shù)，計(jì)算CCD光敏面各亞像元區(qū)域中心的觀測向量，根據(jù)鏡頭畸變參數(shù)及亞像元中心觀測光線和主光軸夾角，計(jì)算該視場的鏡頭畸變量，并對先前計(jì)算的觀測向量進(jìn)行修正。

(5)根據(jù)相機(jī)安裝矩陣、衛(wèi)星姿態(tài)角、衛(wèi)星位置速度向量，將觀測向量轉(zhuǎn)換到數(shù)字地面模型相同的坐標(biāo)系，獲取地面元在相機(jī)入瞳處的輻亮度。

(6)根據(jù)基礎(chǔ)影像MTF曲線，目標(biāo)影像MTF曲線(PSF離散二維矩陣)和相機(jī)相對孔徑、地面分辨率等參數(shù)，計(jì)算相同空間頻率下的MTF對應(yīng)關(guān)系。通過快速傅立葉變換，在頻率域?qū)崿F(xiàn)對MTF導(dǎo)致的輻亮度“影像”的分辨能力下降的模擬。

(7)對各個離散時刻CCD獲取目標(biāo)區(qū)域的平均輻亮度和目標(biāo)相機(jī)A/D轉(zhuǎn)換關(guān)系輸出目標(biāo)影像對應(yīng)行列的DN值。

(8)根據(jù)不同輻亮度對應(yīng)的信噪比計(jì)算不同亮度下的等效噪聲亮度值并根據(jù)目標(biāo)相機(jī)的A/D轉(zhuǎn)換關(guān)系將其轉(zhuǎn)換為噪聲DN值，將噪聲DN值乘以－1～1之間符合噪聲分布特點(diǎn)的隨機(jī)數(shù)作為模擬噪聲附加在第(7)步的輸出影像上，并作溢出判斷。

4 成像仿真

成像鏈模型對構(gòu)成立體測繪相機(jī)成像的各個物理環(huán)節(jié)建模，因此利用上述計(jì)算結(jié)果對系統(tǒng)成像進(jìn)行仿真。采用高精度高分辨率地表物理模型作為仿真目標(biāo)圖像，測試場景數(shù)據(jù)如表1所示，系統(tǒng)參數(shù)如表2所示。圖4中將本文算法計(jì)算的點(diǎn)擴(kuò)散函數(shù)與光學(xué)分析Zemax軟件計(jì)算值進(jìn)行對比，兩者之間偏差 RMS=0.009 6，誤差不超過1%。圖5(a)為26°前視相機(jī)仿真影像，圖5(b)為5°后視相機(jī)仿真影像，圖5(c)為立體影像。將測試場景中心作為地面控制點(diǎn)，根據(jù)光線矢量與地表物理模型求交后交點(diǎn)幾何信息求取定位精度。正視相機(jī)幾何物理模型定位精度達(dá)124 m，前視相機(jī)定位精度達(dá)193 m。建立的立體測繪相機(jī)全鏈路模型滿足幾何精度，能夠合成立體影像，仿真方法可行，能夠滿足立體測繪相機(jī)技術(shù)指標(biāo)論證及成像質(zhì)量評價(jià)的要求。

表1 測試場景數(shù)據(jù)參數(shù)Tab.1 Data parameters of test scenario

表2 立體測繪相機(jī)系統(tǒng)參數(shù)Table 2 Parameters for stereo mapping camera system

圖4 點(diǎn)擴(kuò)散函數(shù)計(jì)算模型Fig.4 Model for calculating the point spread function(PSF)

5 結(jié)論

圖5 仿真影像Fig.5 Simulated scene

基于在軌成像物理機(jī)理的高精度立體測繪相機(jī)建模與仿真在衛(wèi)星的技術(shù)指標(biāo)論證以及成像質(zhì)量評價(jià)中起著關(guān)鍵作用。本文從立體測繪相機(jī)的成像鏈模型出發(fā)，考慮了能量傳輸路徑、衛(wèi)星姿態(tài)軌道參數(shù)、光學(xué)系統(tǒng)成像以及光電轉(zhuǎn)換等環(huán)節(jié)，對成像鏈路中的每一物理環(huán)節(jié)進(jìn)行精確建模與分析。仿真結(jié)果表明:正視相機(jī)幾何物理模型定位精度達(dá)124 m，前視相機(jī)定位精度達(dá)193 m，證明了本文全鏈路模型的可行性及可靠性。

［1］ BORNER A，WIEST L，REULKE R，et al..Sensor:a tool for the simulation of hyperspectral remote sensing systems［J］.ISPRS J.Photogrammetry and Remote Sensing，2001，55(5-6):299-312.

［2］ BORNER A，WIEST L.The simulation of APEX data:the SENSOR approach［J］.SPIE，1999，3780:23-32.

［3］ COTA S A，BELL J T，BOUCHER R H，et al..PICASSO:an end-to-end image simulation tool for space and airborne imaging systems［J］.J.Appl.Remote Sensing，2010，4(1):043535.

［4］ 劉曉.基于低空遙感系統(tǒng)的星載光學(xué)遙感器成像仿真［J］.紅外與激光工程，2014，43(1):217-225.LIU X.Satellite borne optical remote sensor imaging simulation based on low-altitude remote sensing system［J］.Infrared and Laser Engineering，2014，43(1):217-225.(in Chinese)

［5］ 王剛.基于圖像仿真的對地遙感過程科學(xué)可視化研究［J］.系統(tǒng)仿真學(xué)報(bào)，2002，14(6):756-760.WANG G.Study on scientific visualization of earth remote sensing based on imagery simulation［J］.J.System Simulation，2002，14(6):756-760.(in Chinese)

［6］ 陳方，牛錚，覃馭楚，等.基于寬光譜光學(xué)遙感圖像的細(xì)分光譜光學(xué)遙感圖像的模擬［J］.光電工程，2007，34(5):89-96.CHEN F，NIU ZH，QIN Y CH，et al..Simulation of an image with a subsection of spectral band using an image with a wider spectral band［J］.Opto-Electronic Engineering，2007，34(5):89-96.(in Chinese)

［7］ 毛克彪，覃志豪.大氣輻射傳輸模型及MODTRAN中透過率計(jì)算［J］.測繪與空間地理信息，2004，27(4):1-3.MAO K B，QIN ZH H.The transmission model of atmospheric radiation and the computation transmittance of MODTRAN［J］.Geomatics＆ Spatial Information Techonology，2004，27(4):1-3.(in Chinese)

［8］ 亓雪勇，田慶久.光學(xué)遙感大氣校正研究進(jìn)展［J］.國土資源遙感，2005，66(4):1-6.QI X Y，TIAN Q J.The advances in the study of atmospheric correction for optical remote sensing［J］.Remote Sensing for Land ＆ Resources，2005，66(4):1-6.(in Chinese)

［9］ 江月松，李小路，趙一鳴，等.熱紅外對地遙感中的大氣散射效應(yīng)［J］.光學(xué)學(xué)報(bào)，2007，26(12):1766-1771.JIANG Y S，LI X L，ZHAO Y M，et al..Effects of atmospheric scattering on thermal infrared remote sensing of the earth［J］.Acta Optical Sinica，2007，26(12):1766-1771.(in Chinese)

［10］ BERK A，ANDERSON G P，BERNSTEIN L S，et al..MODTRAN4 radiative transfer modeling for atmospheric correction［J］.SPIE，1999，3866:24-193.

［11］ 付強(qiáng)，相里斌，景娟娟，等.基于多光譜遙感成像鏈模型的系統(tǒng)信噪比分析［J］.光學(xué)學(xué)報(bào)，2012，32(2):0211001.FU Q，XIANG L B，JING J J，et al..System signal-to-noise ratio analysis based on image chain model in multispectral remote sensing［J］.Acta Optical Sinica，2012，32(2):0211001.

［12］ 王家騏，于平，顏昌翔，等.航天光學(xué)遙感器像移速度矢計(jì)算數(shù)學(xué)模型［J］.光學(xué)學(xué)報(bào)，2004，24(12):1585-1589.WANG J Q，YU P，YAN C X，et al..Space optical remote sensor image motion velocity vector computational modeling［J］.Acta Optical Sinica，2004，24(12):1585-1589.(in Chinese)

［13］ 薛旭成 ，傅瑤，韓誠山.TDI CCD相機(jī)的衛(wèi)星姿態(tài)穩(wěn)定度確定［J］.中國光學(xué)，2013，6(5):767-772.XUE X CH，F(xiàn)U Y，HAN CH SH.Confirmation of satellite attitude stabilization for TDI CCD camera［J］.Chinese Optics，2013，6(5):767-772.(in Chinese)

［14］ 楊秀彬，常琳，金光.單框架控制力矩陀螺轉(zhuǎn)子動不平衡對遙感衛(wèi)星成像的影響［J］.中國光學(xué)，2012，5(4):358-364.YANG X B，CHANG L，JIN G.Influence of dynamic imbalance of SGCMG rotor on remote sensing satellite imaging［J］.Chinese Optics，2012，5(4):358-364.(in Chinese)