我国首颗民用传输型立体测绘卫星资源三号的发射成功结束了我国以往完全依靠国外高分辨率影像数据的历史。资源三号卫星的成功很大程度上得益于先前对卫星各项技术指标的详细论证，而卫星成像仿真在资源三号的技术指标论证和地面应用系统建设中都发挥了重要作用。在卫星成像仿真方面，2001年，德国宇航中心(DLR)开发了新一代的光学遥感仿真软件SENSOR，首次实现全链路的模拟流程。从几何和辐射两方面，精确模拟链路中每一成像环节，系统地仿真光学、电学成像模型^{[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]}；2003年，法国图卢兹试验室开发基于“三维体”的辐射传输模型DART，成功地模拟出不同大气条件及不同传感器响应条件下的图像^[11]；美国罗彻斯特理工学院的数字成像与遥感实验室开发的遥感成像仿真软件DIRSIG，已经应用于NASA众多项目中，如OrbView-4商业卫星等项目^{[12, 13]}。Murali、Stephane较早进行了CCD相机成像仿真的研究^{[19, 20]}。Jensen通过光线追踪方法进行了影像合成方面的研究^[21]。Guanter实现了EnMAP的成像仿真研究^[16]。Xiao实现了CCD成像的建模与仿真^[17]。Fukue等实现了ALOS卫星的数据仿真^[18]。我国针对成像系统仿真工作起步较晚，长春光机所编制了遥感图像模拟软件RSIS1.0。安徽光机所开发了光学遥感图像仿真软件ORSIS。中国空间技术研究院508所也在进行遥感器成像仿真与评价方面的研究工作，编制图像模拟软件SOSRA2.0。董龙等做了线阵CCD成像仿真方面的研究^[2]。颜昌祥等对一定频率的方波信号产生的静态、动态MTF下降进行了仿真研究^[4]。岳庆兴对卫星相机动态成像造成的几何质量下降进行了仿真研究^{[8, 10]}。这些软件和方法大都局限于卫星影像的辐射质量和几何分辨能力，不能解决亚m级TDI CCD(time delay and iintegration charge coupled devices)立体测绘相机最关注的立体几何精度仿真等问题,因此，国产亚m级立体测绘卫星的成像仿真研究十分迫切。卫星影像的分辨能力是采样间隔、MTF、噪声、平台运动、工作环境等一系列因素综合作用的结果^{[3, 6, 9, 14]}，要达到亚m级分辨率需要综合考虑各种影响因素，设计合理的技术指标和工作模式才能实现。另外，卫星成像过程是几何、辐射因素共同作用的结果，不能将两者分离开来。本文基于精细的地面场景数据，提出了亚m级TDI CCD立体测绘相机的几何、辐射一体化成像仿真技术。 1 精细地面模型与入瞳辐亮度获取

实现亚m级的成像仿真需要建立高分辨率的地面模型，这里称为精细地面模型，包括精细几何模型、纹理模型及两者之间的对应关系。同时需要针对精细地物模型开发新的仿真算法。本文使用的精细地面模型为三角网加贴面纹理影像的数据结构。

卫星相机对地面目标成像实际上是对地面目标辐射亮度的响应或模数转换。本文将贴面纹理影像每个像素的RGB值各乘以一个比例系数来近似代替RGB三个波段内的反射率，同样对RGB波段设置一个波段范围。对某一波长处的反射率乘以辐照度，再除以半球立体角就得到该点的辐射亮度。本文中使用的辐照度查找表从国家质量技术监督局于1999年发布的《太阳能在不同接收条件下的太阳光谱辐照度标准》文件中获取。

得到地面目标的辐亮度后，要得到相应的入瞳辐亮度还要获得大气透过率和大气后向散射辐亮度两个关键物理量。大气透过率受气象视距和观测天底角的影响较大，大气后向散射辐亮度受气象视距、太阳天顶角、观测天底角影响较大，而且随目标反射率变化而变化。同时，两个物理量都是波长的函数，因此也将两个物理量称作波谱透过率和波谱后向散射辐亮度，最终获取的入瞳辐亮度是波谱辐亮度在设置波长范围内的积分辐亮度。设目标波谱辐亮度 为L(λ)，波谱大气后向散射辐亮度为L_p(λ)，这里两个波谱辐亮度单位为W/(m²·sr·μm)波谱大气透过率为τ(λ)，则波段(λ₁、λ₂)(单位为μm)内的入瞳积分辐亮度为：

本文首先通过6S模型计算典型观测条件下的大气透过率和大气后向散射辐亮度查找表。在仿真运算时通过插值获取需要的物理量，保证后面获取的焦面能量在合理的范围之内。 2 建立“视线目标索引”图

成像仿真的一个重要步骤是确定影像点对应的地面目标位置及辐射亮度。其中确定像点对应地面目标的主要工作是解决视线遮挡问题，确定目标的辐射亮度的主要工作是分析该点与太阳之间是否存在遮挡。像点和投影中心连线方向可能存在多个地面点，需要在这多个地面点中找出离视点最近的点，即该像点对应的地面目标点。地面目标点与太阳连线之间同样可能存在多个地面点，只要有一个点遮挡，则太阳无法照射到该目标地面点，该地面目标点就处在阴影中。本文通过建立“视线目标索引”的方法，可以快速确定视线对应的最近目标点并判断该目标点是否对太阳可见。具体算法如下。

1) 不考虑TDI CCD的级数，将每个三角网各三角形顶点坐标根据定向参数投射到模拟影像上，得到像方三角形。对于线阵成像来说，这是一个以轨道和姿态数据为基础的迭代过程^[1]。

2) 计算统计落在像方三角形内的像点坐标，记录三角形内像点的“视线目标索引，即对应的三角网和三角形序号。

3) 以太阳中心作为框幅式成像的投影中心，根据太阳高度角和方位角计算该框幅式相机的姿态角。

4) 将每个三角网的各三角形顶点地面坐标根据定向参数投射到步骤3)的框幅式影像上，得到框幅式影像像方三角形。

5) 统计落在框幅式影像像方三角形内的像点坐标，记录三角形内像点的“视线目标索引”，即对应的三角网和三角形序号。

6) 计算模拟成像时刻t对应的卫星轨道(位置和速度)、姿态，计算TDI CCD上的像点(r，c)与投影中心构成的视线向量 u ，根据相机安装角和卫星轨道、姿态将 u 转换到物方空间得到 u ′。

7) 设时刻t最近的整数影像行对应时刻为t₀、t₀，对应整数行为r′(r′=r或r′=r+1)，单行CCD积分时间为dt，步骤6)中r对应的整数TDI CCD行为r₀，c对应整数TDI CCD列为c₀，则对应的“视线目标索引”对应的整数影像行列(R,C)为：

式中，int表示取整。

8) 在 以(R,C)为中心的3×3窗口内计算 u ′与该窗口内所有“视线目标索引”对应地面三角形的交点序列，统计离t时刻投影中心最近的地面点，记录相应的地面点P、三角网W和三角形T。

9) 将步骤8)求得的地面点P投影到太阳方向的框幅式影像上得到像点(r_s,c_s)，对应视线向量为u_s。以(r_s,c_s) 为中心的3×3窗口内计算u_s与该窗口内所有框幅式影像“视线目标索引”对应地面三角形的交点序列，统计框幅式影像投影中心最近的地面点，记录相应的三角网W′和三角形T′。

10) 如果W=W′且T=T′，则地面点P对相机和太阳同时可见，否则，地面点P对相机可见但对太阳不可见。

11) 设框幅式影像投影中心为(X_s,Y_s,Z_s)，对应物距为d,日地距离为D，太阳直径为d_s。将框幅式影像投影中心在与框幅式影像主光轴垂直且过(X_s,Y_s,Z_s)点的直径为d_s*d/D的圆 平面内作随机移动(可预先确定一定间距doff的均匀分布点，然后在这些点的基础上加入范围在(-doff,-doff)-(doff,doff)之间的随机量),对每个投影中心根据步骤9)~步骤10)判断地面点是否对太阳可见(对应每个投影中心都把太阳看作点光 源)。统计可视光线比例,模拟软阴影效果。 3 静态MTF

由于光学系统像差、装调误差及光的波动性，源于物方同一点的光线到达像方焦面时会形成一个弥散斑，弥散斑的形状和尺度表征光学系统的质量。弥散斑中的光能量分布函数即点扩散函数(Point Spread Function,PSF)。点扩散函数的存在使得像方任一点的辐射能量既包含来自“目标”点的辐射能量，同时包含来自“目标点”周围点的辐射能量。由于每个点的点扩散函数都有一定的作用范围，造成影像的模糊程度也是与空间尺度或频率有关的。通常采用正弦信号作为输入，不考虑点扩散函数的空间和尺度的变化，输出信号(影像)也是正弦信号。点扩散函数的存在造成信号对比度的下降，通常采用调制度传递函数(Modulation Transfer Function,MTF)来描述。设输入正弦信号的最大值为max_i，最小值为min_i，输出正弦信号的最大值为max_o，最小值为min_o，则MTF计算方法为：

式中，MLT_i和MLT_o分别称为输入信号和输出信号的调制度。

本文采用高斯分布的点扩散函数，计算方法为：

式中，(x,y)表示以像方弥散斑中心为原点的像面位置。对像面进行细分，用细分后的离散平面位置近似代替连续的平面位置，细分程度越高，近似程度也越高。首先通过一维正弦信号建立PSF系数σ、细分倍数、极限频率MTF之间的关系，这里的极限频率指两个CCD大小代表的频率，如CCD尺寸为10 um，则极限频率为1/(2×10/1 000)=50线对/mm。

用设计极限频率MTF对应的点扩散函数对理想“像点”构成的焦面影像进行滤波即实现了光学系统MTF的仿真。但对于每个CCD器件来说，进行光电转换的能量是经过光学系统滤波后“落在”CCD矩形面元内部的平均辐射能量。因此在进行PSF空间滤波后还需要计算每个CCD内部的平均辐射照度，这是一个矩形滤波过程。设PSF作用半径为r，如果不进行PSF滤波，该平均能量是由以矩形CCD中心为原点，往矩形CCD外部四个方向扩展r的矩形范围内辐射能量的加权函数。可以通过预先计算权矩阵一步实现PSF滤波和矩形滤波的综合效果。根据点扩散函数作用半径r和矩形CCD尺寸d(设长宽相等)的相对大小，权矩阵计算公式如下：

式中，x为像点离矩形CCD中心的距离。

4 SNR与TDI动态积分

TDI CCD是解决微光摄像或高分辨力、高速CCD相机应用的首选器件。本文重点对TDI CCD的信噪比和多级动态积分进行仿真。CCD器件的噪声包括霰粒噪声、暗电流噪声、复位噪声、放大器噪声、1/f噪声、光子响应非均匀性、量化噪声等^[5]。其中霰粒噪声是CCD本身的物理特性决定的，霰粒噪声对应噪声电荷数与信号电荷数的开平方成正比，如信号电荷数为 n_e，则霰粒噪声的幅值为n_n= 。即实际电荷数在n_e的基础上作-n_n到n_n的上下浮动 。TDI CCD将多极CCD对同一目标产生的电荷累加输出，在积分级数为M的条件下，总电荷数相比单线阵CCD增加了 M倍，而噪声电荷增加倍，两者比值相应增加M/= 倍，这就是TDI CCD能提高信噪比的原因。设相机的相对孔径为F，某一波长范围内的入瞳积分辐射亮度为L(W/m²·sr)，光学透过率为 τ，则视场角α对应的焦面辐照度为：

对视场角较小的卫星相机，cos⁴α≈1，这时可近似地认为：

CCD光敏面面积为 A，像面在积分时间T内接受的光能量为：

根据量子理论，照射到光敏面上的光子数为：

式中，λ为光的平均波长；h为量子效率；c为真空中的光速。

假设光电转换效率为η，产生电荷数为n_e=ηn，霰粒噪声的幅值为；在积分级数为M的条件下，n_e=Mηn，相应地，霰粒噪声幅值为。

电荷最后经增益、模数转换等过程输出为数字信号，即DN值为：

式中,GRD为量化等级(如8 bits，10 bits等)； n_ful为饱和电子数；G为增益；n_e为总的光生电荷数。

最终的仿真数字影像DN值与光生电荷数成正比，也就是与入瞳辐亮度成正比。在本文中，仿真数字影像DN值与对应贴面纹理的像素灰度值成正比,因此，在实际的仿真运算之前首先根据式(1)及式(6)~式(10)描述的严密物理转换过程计算仿真数字影像DN值与对应地面模型贴面纹理的像素灰度值之间的比例系数,而不需要对每一个点都进行严密的转换。

将行积分时间dt细分成多个离散时刻，计算各个时刻的“静态焦面能量影像”，然后将它们相加得到“行积分时间平均焦面能量影像”。每得到一个行积分时间段内的“行积分时间平均焦面能量影像”。根据信号的大小对“行积分时间平均焦面能量影像”的每个“像点”加入霰粒噪声。将第一级CCD与后面依次推扫过地面目标的各级CCD的“行积分时间平均焦面能量影像”累加作为一行输出影像，实现TDI积分的仿真。

5 多视角立体成像仿真

卫星立体成像主要通过卫星整星在沿轨或穿轨方向摆动一定角度和多个不同安装角的相机沿轨推扫两种方式实现，因此，立体成像仿真需要综合考虑轨道、姿态、相机安装等多个因素。

5.1 轨道姿态模型

本文轨道数据采用文献^[7]的简化轨道模型获取。通过两步获取卫星姿态参数：首先用功率谱分析法获取符合卫星设计姿态稳定度的姿态曲线，然后生成符合卫星设计姿态指向精度的系统指向误差，将指向误差参数和姿态参数相加得到最终的姿态参数。功率谱分析法获取姿态曲线的方法参考文献^[8],在此不再赘述。

5.2 TDICCD畸变模型

相机垂直TDI积分方向的畸变参数为“CCD列－畸变量(垂直TDI积分方向σ_CS，沿TDI积分方向σ_AS)”数据序列，含输入参数的CCD之间的CCD畸变通过线性插值的方法获得。设TDI CCD焦面上有一 点p(c,r)(c为列坐标，r为行坐标单位像素)，过p的主光线对应物方点为P，光学系统入瞳中心为O。点p加入畸变量后变为p′(c′,r′)(c′=c+σ_CS,r′=r+σ_AS)，过p′并与主光轴平行的光线与相机光学系统物方节面的交点为C。则C、O、P三点 共线，根据这一关系可以得到过CCD焦面上任意一点的主光线对应的地面点。

5.3 相机安装模型

相机安装模型包括相机光学节点在卫星本体坐标系的安装位置 (X_s,Y_s,Z_s)和相机测量坐标系相对卫星本体坐标系的安装角(yaw_s,pitch_s,roll_s)。相机测量坐标系相对卫星本体坐标系的旋转矩阵为 r _set= r _ps× r _rs× r _ys。其中，r _ps、 r _rs、 r _ys分别为安装角(yaw_s,pitch_s,roll_s)对应的旋转矩阵。

由于地面模型所在的坐标系是WGS84地心直角坐标系，因此,需要将相机节点的运行轨迹转换到WGS84地心直角坐标系。设卫星本体坐标系原点到( X_s,Y_s,Z_s)的向量为 U ，则有：

卫星本体坐标系相对t时刻的局部轨道坐标系的姿态角为(yaw,pitch,roll)，对应旋转矩阵为 r ^t_a：

式中，r _pitch、 r _roll、 r _yaw分 别为姿态角为(yaw,pitch,roll)对应的旋转矩阵。

卫星本体坐标系原点在WGS84地球固连坐标系的位置速度矢量为 (X_o^t,Y_o^t,Z_o^t,V_Xo^t,V_Yo^t,V_Zo^t)，对应旋转矩阵为 r _o^t：

式中，[X_o^t,Y_o^t,Z_o^t]^T为卫星位置归一化向量； [V_Xo^t,V_Yo^t,V_Zo^t]^T为卫星速度归一化向量； [C_x,C_y,C_z^T为[X_o^t,Y_o^t,Z_o^t]^T和[V_Xo^t,V_Yo^t,V_Zo^t]^T的叉积。

则相机光学节点在卫星本体坐标系的安装向量 U 在WGS84地心直角坐标系的向量形式为：

相应地，相机光学节点的运行位置为：

式中，[X_s^t,Y_s^t,Z_s^t]^T为相机光学节点在卫星本体坐标系的安装位置向量在地球固连坐标系的形式。

设相机成像的行积分时间为Δt，t+Δt时刻的相机光学节点的运行位置为[X_opt^t+Δt,Y_opt^t+Δt,Z_opt^t+Δt]^T，t时刻相机光学节点的速度矢量计算方法为：

已知光学节点位置和光线向量，光线与地面模型的三角面片交点采用将直线方程转换为参数方程形式，然后代入三角面片所在平面方程的方法计算获得。

6 结 语

本文建立了亚m级TDI CCD立体测绘相机的成像仿真技术流程。针对精细地面模型，通过建立“视线目标索引”实现了阴影的仿真；建立了点扩散函数系数、CCD细分倍数和极限频率MTF之间的对应关系，并提出了静态MTF的快速仿真算法；通过分析入瞳辐亮度到最终DN值的转换流程结合霰粒噪声的生成机制，严格计算目标影像灰度值和仿真影像DN值间的转换系数，实现了霰粒噪声仿真。通过获取“行积分时间平均焦面能量影像”实现TDI积分的仿真；提出在获取卫星姿态数据、CCD畸变参数的条件下建立相机安装模型的详细方法，实现了任意侧摆角度和相机安装角下的成像仿真。