CN103913148B - 航天tdi ccd相机全链路数值仿真方法 - Google Patents

Abstract

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

Description

航天TDI CCD相机全链路数值仿真方法

技术领域

本发明涉及光电成像技术领域，具体涉及一种基于在轨成像物理机理的航天TDICCD相机全链路数值仿真方法。

背景技术

随着卫星遥感影像在各领域的广泛应用，对航天TDI CCD相机成像质量要求不断提高。考虑到航天产品成本很高的情况，准确的地面仿真显得尤为必要。高精度高分辨率卫星成像仿真在卫星的技术指标论证、试验测试、在轨运行评价以及故障模拟分析中都能够发挥十分重要的作用。目前，卫星对地成像过程仿真大多采用地面全物理成像试验手段：根据卫星在轨飞行状态和相机参数计算光学遥感器和成像目标之间的相对位置及相对速度关系，利用三轴气浮转台模拟卫星平台，动态靶标模拟地面景物，并使用光纤陀螺、磁强计等实时测量姿态角和姿态角速度信息。然而，这种全物理硬件仿真具有如下不足：

1、全物理仿真使用硬件实物模拟卫星在轨成像状态，成本较高、调试周期长；

2、全物理仿真平台大多采用平面靶标模拟地面景物，或是带有立体起伏地形地物的成像靶标，很难反映出地物各点位置和像面各点的精确投影关系，尤其是难以模拟地球曲率对光学遥感器几何定位及成像畸变的影响；

3、不同类型的传感器所需采用的硬件设备也不尽相同，因此单机复用性较差，带来资源浪费。

发明内容

为了解决上述航天TDI CCD相机全链路全物理硬件仿真存在的技术问题，本发明提出一种针对航天TDI CCD相机在轨成像过程，采用数值模拟技术路线对成像中的每一环节进行精确建模。

航天TDI CCD相机全链路数值仿真方法，包括以下步骤：

步骤一，建立包含波谱属性的观测场景；

首先几何模型的建立采用三角网DEM的组织形式，三角网DEM包括三角形顶点位置、顶点纹理坐标、顶点法向量、三角形顶点在顶点数组中的索引，整个观测场景由一系列三角网DEM组成；

其次将三角网DEM和波谱属性通过唯一ID进行关联，波谱数据库包含的数据组成：观测地点信息描述观测地的经纬度、海拔高度以及面积信息；观测目标信息描述观测的典型地物的名称、所属类别；反射率值描述了此典型地物在不同波长上的反射率，通过索引ID读取相应的波谱数据；

步骤二，轨道建模；

首先给定初始时刻卫星平台在J2000惯性坐标系内的位置和速度；然后基于J2000惯性坐标系内的二体轨道动力学微分方程，利用数值积分算法，对二体轨道动力学微分方程积分，得到下一时刻卫星在惯性坐标系下的位置矢量和速度矢量；最后根据卫星在J2000惯性坐标系下的位置速度矢量计算相应的轨道六根数包括半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角及真近点角；

步骤三，姿态建模；

首先给定初始时刻卫星本体坐标系相对J2000.0惯性坐标系的姿态四元数及体坐标系相对于J2000.0惯性坐标系的姿态角速度在体系下的分量形式；然后对卫星姿态动力学模型和卫星姿态四元数运动学方程数值积分，得到下一时刻体系下卫星本体坐标系相对惯性坐标系的姿态角速度及体坐标系相对惯性坐标系的姿态四元数；

步骤四，相机光学系统建模；

相机内方位元素由焦距、主点位置、图像畸变组成，通过平行光管和星点靶标对相机进行标定得到；将光学系统各光学元件的加工误差，装调误差输入到ZEMAX光学设计软件中，获得相应的二维点扩散函数矩阵；

基于光线追踪法建模中光谱响应特性的模型；设光线入射到光学反射表面的入射角为i，可通过光线追踪法算出；入射波长为λ，则光线经过光学反射表面后的光谱响应特性：

R(i,λ)＝r(i)∫σ(λ)dλ

其中r(i)为镀膜材料对不同角度入射光线的反射率；σ(λ)为镀膜材料的光谱响应函数；设光学系统有n个反射表面，则该光线由入瞳到像面的光谱响应特性为：

$R={\mathrm{\Σ}}_j^n{R}_j(i,\mathrm{\λ})$

采用菲德尔法对光学系统光路进行几何建模，给定光学系统的结构参数，再给出入射到光学系统的光线方向和空间位置，最后求出光线通过该系统后的方向和空间位置；

步骤五，星体观测矢量建模；

首先根据目标相机的探测器件安装参数及步骤四中建立的相机内方位参数和光学系统光线几何模型，计算TDI CCD光敏面各亚像元区域中心的观测向量，根据镜头畸变参数及亚像元中心观测光线和主光轴夹角，计算该视场的镜头畸变量，并对先前计算的观测向量进行修正；然后根据J2000.0惯性坐标系至WGS84坐标系时的坐标转换矩阵模型，结合载荷在卫星平台的安装矩阵及地面目标点在WGS84系下的经纬度坐标，带入步骤二中下一时刻的轨道位置参数及步骤三中卫星平台的姿态参数，可计算得到成像仿真所需的相机与目标区域间的相对几何参数，即光线矢量数据；

在光线追踪算法中，将光线矢量数据和步骤一建立的观测场景数据库进行求交运算，也就是计算观测光线方向上离视点最近的三角网与光线的交点，求得交点后通过索引ID号得到该三角网的波谱信息；

步骤六，目标相机入瞳辐亮度场建模；

根据目标相机成像时的大气条件、太阳高度角、太阳方位角、卫星高度角和卫星方位角参数，通过大气辐射传输模型模拟大气透过率、大气后向散射和大气邻近效应，获取地面元在相机入瞳处的辐亮度；

步骤七，相机辐射响应建模；

对各个离散时刻各级CCD获取目标区域的平均辐亮度和目标相机A/D转换关系输出目标影像对应行列的DN值，得到数字影像；

最终到达探测器靶面上的辐通量可表示为

$\mathrm{\Φ}=\frac{S_{\det ector}{S}_{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为光学系统孔径面积，f为相机焦距，M为垂轴放大率，τ_optics(λ)为光谱透过率，h(x,y,λ)为系统点扩散函数，L_pupil(x,y,λ)为相机入瞳辐亮度，由步骤六得到；

积分时间t内，探测器接收到的电子数为

$N_p=t\frac{S_{\det ector}{S}_{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导致的步骤七中辐亮度影像的分辨能力下降的模拟；

步骤九，加噪及压缩模拟；

将模拟噪声信号加到步骤八得到的无噪声信号上，判断是否超出最大灰度值，若超出设为最大值；霰粒噪声近似用离散泊松分布函数表示；暗电流噪声的概率密度分布为泊松分布；转移噪声可表示为平均转移损失电荷数的均方根值；复位噪声可表示为k为波尔兹曼常数，q＝1.602×10^-19C，T为TDI CCD工作温度；量化噪声可表示为Num_well为满阱电荷数，N为量化位数。

本发明的有益效果：

1、节约了仿真成本，并适用于不同单机，解决了全物理仿真应用单一的问题，提高了资源利用率。

2、排除了全物理仿真平台自身引入的误差对成像的影响，准确计算出不同成像时刻传感器像面各点与相应地物点的投影关系。

3、利用包含实际地物信息的高分辨率高精度地面场模型作为输入源，能够反映成像过程中的几何畸变问题。

4、本发明包含了光学系统模块、轨道姿态模拟单元及相机辐射响应模块等，每个模块都预备接口，便于根据实际科研需求进行变更升级。

附图说明

图1为本发明航天TDI CCD相机全链路数值仿真方法流程图。

图2为空间相机成像几何关系示意图。

图3为本发明建立的星体观测矢量软件。

具体实施方式

本发明的特征在于采用数值模拟技术路线对航天TDI CCD相机在轨成像过程的每一个环节进行精确建模，技术实施方法如图1所示，具体通过以下步骤实现：

步骤一：建立包含波谱属性的观测场景

观测场景的建模分为几何和辐射特性两部分，首先几何模型的建立采用三角网数字高程模型(DEM)的组织形式，三角网DEM包括三角形顶点位置、顶点纹理坐标、顶点法向量、三角形顶点在顶点数组中的索引，整个观测场景由一系列三角网DEM组成。

其次将三角网DEM和波谱属性通过唯一ID进行关联，波谱数据库包含的数据组成：观测地点信息描述观测地的经纬度、海拔高度以及面积信息；观测目标信息描述观测的典型地物的名称、所属类别；反射率值描述了此典型地物在不同波长上的反射率，通过索引ID读取相应的波谱数据。

步骤二：轨道建模

首先给定初始时刻卫星平台在J2000惯性坐标系内的位置和速度；然后基于J2000惯性坐标系内的二体轨道动力学微分方程，利用数值积分算法，对二体轨道动力学微分方程积分，得到下一时刻卫星在惯性坐标系下的位置矢量和速度矢量；最后根据卫星在J2000惯性坐标系下的位置速度矢量计算相应的轨道六根数(半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角及真近点角)。

步骤三：姿态建模

首先给定初始时刻卫星本体坐标系相对J2000.0惯性坐标系的姿态四元数及体坐标系相对于J2000.0惯性坐标系的姿态角速度在体系下的分量形式；然后对卫星姿态动力学模型和卫星姿态四元数运动学方程数值积分，得到下一时刻体系下卫星本体坐标系相对惯性坐标系的姿态角速度及体坐标系相对惯性坐标系的姿态四元数。

步骤四：相机光学系统建模

相机光学系统建模包含两大类：第一类是采用理论公式推导的方法对相机光学系统基础理论建模，包括相机内方位元素建模、光学系统点扩散函数建模、入瞳到像面的光谱辐射响应特性建模；第二类是采用光线追迹法对光学系统光路几何建模。

相机内方位元素由焦距、主点位置、图像畸变组成，通过平行光管和星点靶标对相机进行标定得到；将光学系统各光学元件的加工误差，装调误差输入到ZEMAX光学设计软件中，获得相应的二维点扩散函数矩阵；本发明给出基于光线追踪法建模中光谱响应特性的模型。设光线入射到光学反射表面的入射角为i，可通过光线追踪法算出；入射波长为λ，则光线经过光学反射表面后的光谱响应特性：

R(i,λ)＝r(i)∫σ(λ)dλ

其中r(i)为镀膜材料对不同角度入射光线的反射率；σ(λ)为镀膜材料的光谱响应函数。设光学系统有n个反射表面，则该光线由入瞳到像面的光谱响应特性为：

$R={\mathrm{\Σ}}_j^n{R}_j(i,\mathrm{\λ})$

本发明中采用菲德尔(D.Feder)法对光学系统光路进行几何建模，给定光学系统的结构参数，如半径、厚度或间隔、折射率等结构参数，再给出入射到光学系统的光线方向和空间位置(也就是目标的位置)，最后求出光线通过该系统后的方向和空间位置。

步骤五：星体观测矢量建模

空间相机成像几何关系如图2所示。首先根据目标相机的探测器件安装参数及步骤四中建立的相机内方位参数和光学系统光线几何模型，计算CCD光敏面各亚像元区域中心的观测向量，根据镜头畸变参数及亚像元中心观测光线和主光轴夹角，计算该视场的镜头畸变量，并对先前计算的观测向量进行修正；然后根据J2000.0惯性坐标系至WGS84坐标系时的坐标转换矩阵模型，结合载荷在卫星平台的安装矩阵及地面目标点在WGS84系下的经纬度坐标，带入步骤二中下一时刻的轨道位置参数及步骤三中卫星平台的姿态参数，可计算得到成像仿真所需的相机与目标区域间的相对几何参数，即光线矢量数据。

在光线追踪算法中，将光线矢量数据和步骤一建立的观测场景数据库进行求交运算，也就是计算观测光线方向上离视点最近的三角网与光线的交点，求得交点后通过索引ID号得到该三角网的波谱信息。本发明中建立的星体观测矢量软件如图3所示。

步骤六：目标相机入瞳辐亮度场建模

根据目标相机成像时的大气条件、太阳高度角、太阳方位角、卫星高度角、卫星方位角等参数，通过大气辐射传输模型(MODTRAN4.0)模拟大气透过率、大气后向散射、大气邻近效应，获取地面元在相机入瞳处的辐亮度。

步骤七：相机辐射响应建模

对各个离散时刻各级CCD获取目标区域的平均辐亮度和目标相机A/D转换关系输出目标影像对应行列的DN值，得到数字影像。

最终到达探测器靶面上的辐通量可表示为

$\mathrm{\Φ}=\frac{S_{\det ector}{S}_{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为光学系统孔径面积，f为相机焦距，M为垂轴放大率，τ_optics(λ)为光谱透过率，h(x,y,λ)为系统点扩散函数，L_pupil(x,y,λ)为相机入瞳辐亮度，由步骤六得到。

积分时间t内，探测器接收到的电子数为

$N_p=t\frac{S_{\det ector}{S}_{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曲线，目标影像MTF曲线(PSF离散二维矩阵，由步骤四得到)和相机的相对孔径、地面分辨率等参数，计算相同空间频率下的MTF对应关系。通过快速傅立叶变换，在频率域实现对MTF导致的步骤七中辐亮度“影像”的分辨能力下降的模拟。

步骤九：加噪及压缩模拟

TDI CCD成像过程中信号电荷会受到霰粒噪声、非均匀性噪声、暗电流噪声、转移噪声、复位噪声、量化噪声的影响。将模拟噪声信号加到步骤八得到的无噪声信号上，判断是否超出最大灰度值，若超出设为最大值。霰粒噪声近似用离散泊松分布函数表示；暗电流噪声的概率密度分布为泊松分布；转移噪声可表示为平均转移损失电荷数的均方根值；复位噪声可表示为k为波尔兹曼常数，q＝1.602×10^-19C，T为TDI CCD工作温度；量化噪声可表示为Num_well为满阱电荷数，N为量化位数。

采用JPEG2000压缩方法对影像进行压缩处理。

Claims (1)

1.航天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为量化位数。