基于在轨成像物理机理的立体测绘相机建模与仿真

杨 飞，金 光 ，谢金华，邱振戈，曲宏松，贺小军

(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)