基于在軌成像物理機理的立體測繪相機建模與仿真

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

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

1 引言

由于航天任務具有高風險、高投入、研制周期長等特點，因此在地面實驗環境下進行光學成像建模與仿真是完成載荷和平臺大量測試工作的重要途徑。高精度、高分辨率立體測繪相機成像仿真在衛星技術指標論證、試驗測試、在軌運行評價以及故障模擬分析中都發揮著十分重要的作用。精確的成像仿真可以完善衛星平臺及載荷指標系統，使得在研衛星性能達到最優，從而保證衛星在軌成像質量。因此，基于在軌成像物理機理的全鏈路仿真在衛星研制過程中占據著舉足輕重的地位。

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

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

2 立體測繪相機成像鏈模型

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

基于在軌成像物理機理對立體測繪相機進行建模與仿真時，需要分析平臺姿態軌道、光學系統、相機輻射響應等主要成像鏈環節。

2.1 姿態軌道建模

平臺姿態及軌道動力學仿真模型是整個動態成像系統中的基本組成部分。軌道動力學模型主要模擬衛星軌道運動，計算軌道信息;姿態動力學模型主要模擬衛星姿態運動，計算姿態信息。該仿真模型基本計算過程如下:

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

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

(3)由衛星在J2000.0慣性坐標系下的位置速度矢量，計算相應的軌道根數(半長軸、偏心率、軌道傾角、近心點角距、升交點赤經、偏近點角及平近點角);

(4)給定初始時刻衛星本體坐標系相對于J2000.0慣性坐標系的姿態四元數，下一時刻的姿態角信息由三次多項式得到;

(5)根據所建立的 J2000.0慣性坐標系至WGS84坐標系時的坐標轉換矩陣模型，結合載荷在衛星平臺的安裝矩陣及地面目標點在WGS84系下的經緯度坐標，帶入下一時刻的軌道位置參數及衛星平臺的姿態參數，可計算得到成像仿真所需的前視、正視相機與目標區域間的相對幾何參數，為基于光線追跡的成像仿真模型提供衛星姿態及軌道輸入參數;

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

通過上述方法得到衛星的外方位元素后，根據文獻［12］中提到的坐標轉換方法，計算正視相機和前視相機像面上各點的光線在WGS84坐標系下的矢量方向。

2.2 光學系統建模

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

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

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

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

式中，βsca為大氣散射系數。

根據CCD上光線出發點的密度將光學系統PSF離散化，取與CCD平均間隔相同，設PSF有效寬度為w，CCD光敏面尺寸為d×d，光線密度為m×m，則PSF的離散間隔為:

計算綜合PSF，以一維的點擴散函數為例，綜合PSF計算方法為:

式中，w為點擴散函數的有效范圍，pixsize為像元尺寸。

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

圖1 光學系統光路示意圖Fig.1 Light path diagram of optical system

2.3 相機輻射響應建模

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

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

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

電壓輸出的響應度R為:

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

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

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

模擬電壓信號VS經過模數轉換電路進行灰度量化，假設量化位數為 N，最大量化電壓為Vmax，最小量化電壓為Vmin，則得到系統的輸出圖像的灰度滿足如下關系:

2.3.2 焦平面光譜響應幅值分布特性模型

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

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

2.3.3 焦平面光譜響應SNR模型

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

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

3 仿真方法與分析

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

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

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

(2)根據衛星軌道高度、軌道傾角、偏心率、降交點時刻及成像目標位置，綜合考慮星載相機的視場角范圍，計算能夠對目標區域中心成像的衛星運行時刻，輸出這一時刻前后數組衛星位置、速度向量及對應時間。

(3)根據衛星姿態穩定度、姿態指向精度、偏流角修正精度等設計參數及不同顫振分量的頻率、幅值的統計信息［13-14］，生成瞬時軌道坐標系下的衛星姿態角。

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

(5)根據相機安裝矩陣、衛星姿態角、衛星位置速度向量，將觀測向量轉換到數字地面模型相同的坐標系，獲取地面元在相機入瞳處的輻亮度。

(6)根據基礎影像MTF曲線，目標影像MTF曲線(PSF離散二維矩陣)和相機相對孔徑、地面分辨率等參數，計算相同空間頻率下的MTF對應關系。通過快速傅立葉變換，在頻率域實現對MTF導致的輻亮度“影像”的分辨能力下降的模擬。

(7)對各個離散時刻CCD獲取目標區域的平均輻亮度和目標相機A/D轉換關系輸出目標影像對應行列的DN值。

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

4 成像仿真

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

表1 測試場景數據參數Tab.1 Data parameters of test scenario

表2 立體測繪相機系統參數Table 2 Parameters for stereo mapping camera system

圖4 點擴散函數計算模型Fig.4 Model for calculating the point spread function(PSF)

5 結論

圖5 仿真影像Fig.5 Simulated scene

基于在軌成像物理機理的高精度立體測繪相機建模與仿真在衛星的技術指標論證以及成像質量評價中起著關鍵作用。本文從立體測繪相機的成像鏈模型出發，考慮了能量傳輸路徑、衛星姿態軌道參數、光學系統成像以及光電轉換等環節，對成像鏈路中的每一物理環節進行精確建模與分析。仿真結果表明:正視相機幾何物理模型定位精度達124 m，前視相機定位精度達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］ 劉曉.基于低空遙感系統的星載光學遙感器成像仿真［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］ 王剛.基于圖像仿真的對地遙感過程科學可視化研究［J］.系統仿真學報，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］ 陳方，牛錚，覃馭楚，等.基于寬光譜光學遙感圖像的細分光譜光學遙感圖像的模擬［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］.測繪與空間地理信息，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］ 亓雪勇，田慶久.光學遙感大氣校正研究進展［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］ 江月松，李小路，趙一鳴，等.熱紅外對地遙感中的大氣散射效應［J］.光學學報，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］ 付強，相里斌，景娟娟，等.基于多光譜遙感成像鏈模型的系統信噪比分析［J］.光學學報，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］ 王家騏，于平，顏昌翔，等.航天光學遙感器像移速度矢計算數學模型［J］.光學學報，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］.中國光學，2013，6(5):767-772.XUE X CH，FU Y，HAN CH SH.Confirmation of satellite attitude stabilization for TDI CCD camera［J］.Chinese Optics，2013，6(5):767-772.(in Chinese)

［14］ 楊秀彬，常琳，金光.單框架控制力矩陀螺轉子動不平衡對遙感衛星成像的影響［J］.中國光學，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)