[发明专利]航天TDI CCD相机全链路数值仿真方法

摘要

基于在轨成像物理机理的航天TDICCD相机全链路数值仿真方法，属于光电成像技术领域，为了解决卫星全链路全物理硬件仿真存在的技术问题，本发明方法包括以下步骤：步骤一，建立包含波谱属性的观测场景；步骤二，轨道建模；步骤三，姿态建模；步骤四，相机光学系统建模；步骤五，星体观测矢量建模；步骤六，目标相机入瞳辐亮度场建模；步骤七，相机辐射响应建模；步骤八，MTF退化模拟；步骤九，加噪及压缩模拟；本发明节约了仿真成本；排除了全物理仿真平台自身引入的误差对成像的影响，准确计算出不同成像时刻传感器像面各点与相应地物点的投影关系；能够反映成像过程中的几何畸变问题；每个模块都预备接口，便于根据实际科研需求进行变更升级。

主权项

航天TDI CCD相机全链路数值仿真方法，其特征是，包括以下步骤：步骤一，建立包含波谱属性的观测场景；首先几何模型的建立采用三角网DEM的组织形式，三角网DEM包括三角形顶点位置、顶点纹理坐标、顶点法向量、三角形顶点在顶点数组中的索引，整个观测场景由一系列三角网DEM组成；其次将三角网DEM和波谱属性通过唯一ID进行关联，波谱数据库包含的数据组成：观测地点信息描述观测地的经纬度、海拔高度以及面积信息；观测目标信息描述观测的典型地物的名称、所属类别；反射率值描述了此典型地物在不同波长上的反射率，通过索引ID读取相应的波谱数据；步骤二，轨道建模；首先给定初始时刻卫星平台在J2000.0惯性坐标系内的位置和速度；然后基于J2000.0惯性坐标系内的二体轨道动力学微分方程，利用数值积分算法，对二体轨道动力学微分方程积分，得到下一时刻卫星在惯性坐标系下的位置矢量和速度矢量；最后根据卫星在J2000.0惯性坐标系下的位置速度矢量计算相应的轨道六根数包括半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角及真近点角；步骤三，姿态建模；首先给定初始时刻卫星本体坐标系相对J2000.0惯性坐标系的姿态四元数，卫星本体坐标系相对于J2000.0惯性坐标系的姿态角速度在体系下的分量形式；然后对卫星姿态动力学模型和卫星姿态四元数运动学方程数值积分，得到下一时刻体系下卫星本体坐标系相对J2000.0惯性坐标系的姿态角速度及卫星本体坐标系相对J2000.0惯性坐标系的姿态四元数；步骤四，相机光学系统建模；相机内方位元素由焦距、主点位置、图像畸变组成，通过平行光管和星点靶标对相机进行标定得到；将光学系统各光学元件的加工误差，装调误差输入到ZEMAX光学设计软件中，获得相应的二维点扩散函数矩阵；基于光线追踪法建模中光谱响应特性的模型：设光线入射到光学反射表面的入射角为i，可通过光线追踪法算出；入射波长为λ，则光线经过光学反射表面后的光谱响应特性：R(i,λ)＝r(i)∫σ(λ)dλ其中r(i)为镀膜材料对不同角度入射光线的反射率；σ(λ)为镀膜材料的光谱响应函数；设光学系统有n个反射表面，则该光线由入瞳到像面的光谱响应特性为：$R={\mathrm{\Σ}}_j^n{R}_j(i,\mathrm{\λ})$采用菲德尔法对光学系统光路进行几何建模，给定光学系统的结构参数，再给出入射到光学系统的光线方向和空间位置，最后求出光线通过该系统后的方向和空间位置；步骤五，星体观测矢量建模；首先根据目标相机的探测器件安装参数及步骤四中建立的相机内方位元素和光学系统光线几何模型，计算CCD光敏面各亚像元区域中心的观测向量，根据镜头畸变参数及亚像元中心观测光线和主光轴夹角，计算该视场的镜头畸变量，并对先前计算的观测向量进行修正；然后根据J2000.0惯性坐标系至WGS84坐标系时的坐标转换矩阵模型，结合载荷在卫星平台的安装矩阵及地面目标点在WGS84系下的经纬度坐标，带入步骤二中下一时刻的轨道位置参数及步骤三中卫星平台下一时刻的姿态参数，可计算得到成像仿真所需的相机与目标区域间的相对几何参数，即光线矢量数据；在光线追踪算法中，将光线矢量数据和步骤一建立的观测场景数据库进行求交运算，也就是计算观测光线方向上离视点最近的三角网与光线的交点，求得交点后通过索引ID号得到该三角网的波谱信息；步骤六，目标相机入瞳辐亮度场建模；根据目标相机成像时的大气条件、太阳高度角、太阳方位角、卫星高度角和卫星方位角参数，通过大气辐射传输模型模拟大气透过率、大气后向散射和大气邻近效应，获取地面元在相机入瞳处的辐亮度；步骤七，相机辐射响应建模；对各个离散时刻各级CCD获取目标区域的平均辐亮度和目标相机A/D转换关系输出目标影像对应行列的像元亮度值，得到数字影像；最终到达探测器靶面上的辐通量可表示为$\mathrm{\Φ}=\frac{S_{\det ector}{\mathrm{\β}}_{aperture}{\cos }^n\mathrm{\θ}}{f^2{(M+1)}^2}\×\underset{{\mathrm{\λ}}_{\min }}{\overset{{\mathrm{\λ}}_{\max }}{\∫}}{\mathrm{\τ}}_{optics}\left(\mathrm{\λ}\right)\left(L_{pupil}\right(x,y,\mathrm{\λ})*h(x,y,\mathrm{\λ}\left)\right)d\mathrm{\λ}$式中：S_detector为像元有效面积，S_aperture为光学系统孔径面积，θ为入射角，n为光学系统反射面个数，f为相机焦距，M为垂轴放大率，τ_optics(λ)为光谱透过率，h(x,y,λ)为系统点扩散函数，L_pupil(x,y,λ)为相机入瞳辐亮度，由步骤六得到；积分时间t内，探测器接收到的电子数为$N_p=t\frac{S_{\det ector}{\mathrm{\β}}_{aperture}{\cos }^n\mathrm{\θ}}{f^2{(M+1)}^2}\×\underset{{\mathrm{\λ}}_{\min }}{\overset{{\mathrm{\λ}}_{\max }}{\∫}}\frac{\mathrm{\λ}}{hc}{\mathrm{\η\τ}}_{optics}\left(\mathrm{\λ}\right)\left(L_{pupil}\right(x,y,\mathrm{\λ})*h(x,y,\mathrm{\λ}\left)\right)d\mathrm{\λ}$式中：h为普朗克常量，c为光速，η为量子效率；步骤八，MTF退化模拟；根据基础影像MTF曲线、由步骤四得到的PSF离散二维矩阵、相机的相对孔径和地面分辨率参数，计算相同空间频率下的MTF对应关系；通过快速傅立叶变换，在频率域实现对MTF导致的步骤七中辐亮度影像的分辨能力下降的模拟；步骤九，加噪及压缩模拟；将模拟噪声信号加到步骤八得到的无噪声信号上，判断是否超出最大灰度值，若超出设为最大值；霰粒噪声近似用离散泊松分布函数表示；暗电流噪声的概率密度分布为泊松分布；转移噪声可表示为平均转移损失电荷数的均方根值；复位噪声可表示为C为扩散电容，k为波尔兹曼常数，q＝1.602×10^‑19C，T为TDI CCD工作温度；量化噪声可表示为Num_well为满阱电荷数，N为量化位数。