CN105528500B - 一种分米级星载tdi ccd立体测绘相机成像仿真方法和系统 - Google Patents

Abstract

本发明公开了一种分米级星载TDI CCD立体测绘相机成像仿真方法和系统，以目前最高分辨率的商业卫星Worldview‑3立体影像作为数据源，通过立体匹配处理成“影像+地面点”的数据形式。以瞬时TDI CCD成像面辐照度的获取为核心步骤，将轨道、姿态参数，相机安装参数，TDI CCD几何参数，相机MTF等有机地串联起来。同时通过时间域积分和TDI CCD多级电荷累加实现TDI CCD多级动态积分的模拟，并以此为基础实现相机辐射响应和线阵影像成像的严格模拟。本发明可以避免传统成像仿真采用单一数据源而难以模拟高分辨率卫星立体成像中不同视角对应几何、辐射信息的差异性的缺点，提高了仿真立体测绘卫星影像的精度。

Description

一种分米级星载TDI CCD立体测绘相机成像仿真方法和系统

技术领域

本发明涉及卫星成像仿真技术领域，特别涉及一种分米级星载TDI CCD立体测绘相机成像仿真方法和系统。

背景技术

随着我国首颗民用传输型立体测绘卫星资源三号的发射成功，结束了我国以往完全依靠国外高分辨率影像数据的历史。大量试验证明，资源三号卫星无论是分辨能力还是几何精度均超过了国外同类卫星。资源三号卫星数据接收以后，地面处理能力也迅速跟进，保障了数据的有效利用。资源三号卫星的成功很大程度上得益于先前对卫星各项技术指标的详细论证，而卫星成像仿真在资源三号的技术指标论证和地面应用系统建设中都发挥了重要作用。卫星成像仿真已经成为现代卫星工程中的一个重要环节。

资源三号卫星正视相机分辨率最高为2.1米，前后视相机分辨率为3.6米，主要面向1∶5万立体测图，对于更大比例尺的测绘应用还不能满足需求。因此，面向1∶1万立体测图的分米级立体测绘卫星也很快被提上议事日程，因此需要研究分米级立体测绘卫星相机的成像仿真技术。

在可见、近红外光学成像卫星的成像仿真技术领域，国内外学者提出了物理仿真、半物理仿真和计算机数值仿真三种技术途径，其中计算机数值仿真的工程成本最低。采用各子系统的实验室标定数据结合设计数据建立各子系统的数学模型，也就是用数学公式描述信号传输和转换的物理过程，简称子系统建模。在子系统建模的基础上实现在轨动态成像过程的全链路的仿真流程，能比较全面地模拟卫星的设计技术指标对卫星在轨成像质量的影响，准确地模拟出延迟积分成像高分辨率光学卫星在轨成像过程中几何和辐射质量的下降与影像总体几何、辐射精度的对应关系，建立起卫星设计参数与在轨成像质量之间的直接关联。

真实的卫星相机成像是对真实的地面物理场景进行辐射信息采集和转化的过程，进行计算机数值仿真同样需要一个虚拟的数字物理场景，也就是一个数值化的地面模型，这个地面模型需要包括地面每一点的几何信息和辐射信息。换句话说，地面模型由几何模型和辐射模型构成，几何模型提供每一地面点的几何信息，辐射模型提供地面每一点的辐射信息。

现有的仿真软件和方法大都局限于卫星影像的辐射质量和几何分辨能力，不能解决分米级TDI CCD立体测绘相机最关注的立体几何精度仿真等问题，因此国产分米级立体测绘卫星的成像仿真研究十分迫切。卫星影像的分辨能力是采样间隔、MTF、噪声、平台运动、工作环境等一系列因素综合作用的结果，要达到分米级分辨率和1∶10000测绘精度仅仅实现CCD对应的地面采样大小达到分米级是远远不够的。需要综合考虑各种影响因素，设计合理的技术指标和工作模式才能实现。另外，卫星成像过程不是简单的几何位置转换或辐射信息转换过程，而是几何、辐射因素共同作用的结果，不能将二者分离开来。立体成像最显著的特点是不同视角获取的影像之间有一定的差异性，体现在几何和辐射两个方面。采用单一的地面几何模型和辐射模型很难对这两个方面进行高精度建模。而视角接近的其他甚高分辨率立体卫星影像本身就包含了这种不同视角间几何信息和辐射信息的差异性，因此可以作为立体成像仿真的数据源。

本发明所要解决的技术问题是：如何提高卫星成像仿真的模拟精度。

发明内容

为此，本发明提出一种分米级星载TDI CCD立体测绘相机成像仿真方法和系统，可充分地消除由于现有技术的限制和缺陷导致的一个或多个问题。

本发明另外的优点、目的和特性，一部分将在下面的说明书中得到阐明，而另一部分对于本领域的普通技术人员通过对下面的说明的考察将是明显的或从本发明的实施中学到。通过在文字的说明书和权利要求书及附图中特别地指出的结构可实现和获得本发明目的和优点。

本发明提供了一种分米级星载TDI CCD立体测绘相机成像仿真方法，其特征在于，包括以下步骤：

步骤S1，获取基础数据，并对所述基础数据进行预处理，生成卫星立体影像每一点的影像-几何模型；其中，所述基础数据为Worldview-3卫星立体影像和辅助RPC参数，所述预处理具体包括：

(S11)，通过多次前方交会和后方交会迭代计算，消除RPC参数的相对误差；

(S12)，通过准核线立体影像制作和半全局匹配算法计算左右影像的视差图，根据所述视差图计算同名像点坐标，进而前方交会计算准核线影像每一点的物方大地坐标(L，B，H)，其中，L表示经度、B表示纬度、H表示椭球高。

步骤S2，计算卫星轨道数据，所述卫星轨道数据包括地球固连坐标系下卫星的扫描时刻、位置、速度数据；

步骤S3，计算卫星姿态数据；

步骤S4，相机几何建模，其具体包括

步骤S401，计算相机安装模型；

步骤S402，计算TDI CCD畸变模型；

步骤S5，利用步骤S1-S4计算的数据，计算TDI CCD瞬时焦面能量影像的数据源范围；

步骤S6，计算滤波前的TDI CCD瞬时焦面辐照度影像；

步骤S7，用点扩散函数对TDI CCD瞬时焦面辐照度影像进行滤波，得到一个时刻的TDI CCD离散平均焦面辐照度；

步骤S8，计算一个积分时间区间内多个细分时刻的TDI辐照度影像，累加后取平均，得到一个积分区间内的TDI时间平均辐照度影像；

步骤S9，将一个积分区间内的TDI时间平均辐照度影像转换为光生电荷数影像，并根据霰粒噪声的生成机制加入噪声电荷，生成该积分区间的电荷数影像；

步骤S10，重复步骤S8和S9，计算下一个积分时间区间的电荷数影像，完成M级电荷数影像累加后，进行模数转换，输出一行灰度值影像；

步骤S11，计算一幅模拟影像的每一行灰度值影像，最终输出一幅模拟影像。

优选的，步骤S1中的子步骤(S11)具体包括：

首先通过SIFT匹配获取Worldview-3全色立体影像间的9个以上(例如，可以为9、16或25个)分布均匀的同名像点；通过同名像点前方交会获取物方地面点坐标；执行后方交会，所述后方交会为以地面点和对应像点作为控制点校正RPC参数的分子常数项和一次项，执行三次以上前方交会、后方交会迭代计算，消除RPC参数的相对误差。

优选的，步骤S2具体包括以下步骤：

步骤S201，计算惯性系轨道数据；

在惯性参考系的XOZ平面上，取一组卫星测量坐标系原点坐标，各原点坐标与惯性参考系原点的距离相同，均为R，R等于平均地球半径加上平均轨道高度；假设第一个原点对应时刻为0，相邻两个质心对应时刻相差0.1秒，由R可求得轨道周期T：

(单位：秒)，

其中G为万有引力常数，等于6.67*10^-11米³/(千克*秒²)，M为地球质量，等于5.98*10²⁴千克，

降交点时刻Tdsd(秒)减去T/4就是第一个质心对应的实际时间T₀，即

T₀＝T_dsd-T/4

以下取样时刻依次加0.1秒为相应的实际时间，dθ＝360/(T*10)为相邻两个原点的地心角，第一个质心在Z轴上(Z＝R)，第i(i＝0，1，…)个取样时刻的卫星原点位置为：

X＝R×cos(θ)

Y＝0

Z＝R×sin(θ) (1)

其中，θ＝(90-dθ×i)/180.0×π。

设轨道倾角为I，将取样坐标绕X轴顺时针旋转I-90度，然后绕Z轴逆时针旋转L度；

步骤S202，计算第i个时刻卫星原点(Xi，Yi，Zi)对应的初始大地坐标(Li，Bi，Hi)；

步骤S203，计算经度L与目标区域中心经度L_c的偏差ΔL；

步骤S204，卫星原点大地坐标的经度L_i(i＝0，1，…)加上ΔL，得到卫星原点的星下点通过目标区域中心点(L_c，B_c)的卫星原点大地坐标序列(L′_i，B_i，H_i)((i＝0，1，…)，L′_i＝L_i+ΔL)；

步骤S205，将卫星原点大地坐标序列(L′_i，B_i，H_i)转换成地心直角坐标序列(X′_i，Y′_i，Z′_i)(i＝0，1，…)；

步骤S206，返回步骤201，利用公式(1)计算卫星原点(Xi，Yi，Zi)时将角度θ加上0.00002度；

步骤S207，重复步骤S202至S206，得到新的卫星原点地心直角坐标序列(X″_i，Y″_i，Z″_i)；

步骤S208，计算卫星原点速度序列(V_xi，V_yi，V_zi)：

V_xi＝(X″_i-X′_i)/(T*0.00002/360)

V_yi＝(Y″_i-Y′_i)/(T*0.00002/360)

V_zi＝(Z″_i-Z′_i)/(T*0.00002/360)

其中第i组卫星原点位置速度对应时间为：T_i＝T₀+i*0.1；

步骤S209，保存卫星原点轨迹数据(T_i，X′_i，Y′_i，Z′_i，V_xi，V_yi，V_zi)。

优选的，步骤S3具体包括以下步骤：

步骤301，计算指向误差参数：

步骤302，利用功率谱分析法获取姿态参数：

步骤S303，将通过步骤S301得到的指向误差参数和通过步骤S302得到的姿态参数相加得到最终的姿态参数。

其中，步骤S303具体包括：

每隔1秒统计生成的姿态角序列f(t_i)的均方差σ，将姿态角序列的每个姿态角乘以σ/(3yaw₁)，然后加上偏航指向角系统量yaw′₀，得到偏航姿态角序列：

yaw(t_i)＝f(t_i)*σ/(3*yaw₁)+yaw′₀

将姿态角序列的每个姿态角乘以σ/(3pitch₁)，然后加上俯仰指向角系统量pitch′₀，得到俯仰姿态角序列：

pitch(t_i)＝f(t_i)*σ/(3*pitch₁)+pitch′₀

将姿态角序列的每个姿态角乘以σ/(3roll₁)，然后加上滚动指向角系统量roll′₀，得到滚动姿态角序列：

roll(t_i)＝f(t_i)*σ/(3*roll₁)+roll′₀。

优选的：

步骤S401具体包括以下子步骤：

S4011，将相机光学节点位置结合成像时刻转换到物方坐标系：

S4012，根据相机安装角计算相机测量坐标系相对卫星本体坐标系的旋转矩阵r_set；

其中，相机测量坐标系相对卫星本体坐标系的旋转矩阵r_set为：

其中，

S4013，计算相机光学节点在卫星本体坐标系的安装向量U在WGS84地球固连坐标系的向量U′：

其中，设相机光学节点在卫星本体坐标系的安装位置为(X_s，Y_s，Z_s)，则U为卫星本体坐标系原点到(X_s，Y_s，Z_s)的向量，为卫星本体坐标系相对t时刻的局部轨道坐标系的旋转矩阵，为卫星本体坐标系原点与WGS84地球固连坐标系对应旋转矩阵。

S4014，得到相机光学节点的运行位置为：

其中，为卫星位置归一化向量，为相机光学节点在卫星本体坐标系的安装位置向量在地球固连坐标系的形式。

S4015，设相机成像的积分时间值为Δt₀，t+Δt₀时刻的相机光学节点的运行位置为t时刻相机光学节点的速度矢量为：

该速度矢量在地面的投影为：

其中H_obt为轨道高度。

设CCD尺寸为a*a，a单位为米，相机主距为f，则地面分辨率GR＝a/f*Hobt。

实际积分时间为：

重新计算t时刻相机光学节点的速度矢量为：

优选的：

每片TDI CCD的安装误差用一个以CCD列坐标为变量的二阶多项式表示，根据相机光学系统的畸变模型将不同视场的畸变量换算为像面位置误差，与安装误差多项式的常数项相加，计算每个CCD在垂直TDI积分方向的畸变量σ_CSi和每个CCD在沿TDI积分方向的畸变量σ_ASi，两个CCD列之间的亚像元点的畸变量为：

σ_AS(c)＝σ_ASi-0.5＜c-i≤0.5

σ_CS(c)＝σ_CS(i+1)0.5＜c-i≤1

其中，i表示CCD序列号，C表示列坐标，i≤c≤i+1。

优选的：

首先将某一时刻的TDI CCD成像面细分为m*m个更小的子CCD，这时整个TDI CCD像面可以看作一幅分辨率为模拟影像分辨率m倍的框幅式影像。设TDI CCD积分级数为M，TDICCD单行探元数为N，则TDI CCD像面构成的框幅式影像(以下称为TDI影像)行数为M*m，列数为N*m；

设TDI影像四角像面坐标为[xi，yi](i＝1，2，3，4)，相机主距为f。则TDI影像四角与光学系统主点构成向量u0＝[xi，yi，-f](i＝1，2，3，4)；

根据相机安装角、卫星平台三轴姿态角、相机主点位置P[X，Y，Z]、速度向量V[vx，vy，vz]计算的相机安装矩阵Mset、姿态矩阵Matt、轨道-物方转换矩阵Mobt将像空间向量u0转换为物方向量u；

设u对应地面点高程为H，根据P、H、u可以计算地面点平面坐标(L，B)。将(L，B，H)通过Worldview-3对应视角的准核线影像RPC计算像面坐标，并获取该点的高程Htrue；

将H替换为Htrue，重复前面的计算过程，直到|H-Htrue|小于一定的限差dh，优选的，所述限差为0.1m；

这样就得到TDI影像四角坐标在Worldview-3对应视角的准核线影像上的四个坐标，也就确定了所需数据的概略范围。将该范围在TDI积分方向对应的Worldview-3准核线影像方向外扩一定距离确定参与计算的数据范围。

优选的：

将确定数据范围内的Worldview-3准核线影像像点对应地面点通过TDI影像的瞬时外方位元素投射到TDI影像上，加上相应的TDI CCD畸变改正量后，按子CCD大小取整后查找坐标对应子CCD的行、列号，每个CCD都有一个存储单元记录投射地面元在Worldview-3准核线影像的坐标和地面元与TDI影像投影中心的距离。所有地面点都完成投射后，有投射记录的子CCD将挑选出距离最近的地面元作为可见地面元，并将对应的Worldview-3准核线影像DN值转换为像面辐照度记录下来，没有投射记录的子CCD辐照度通过周围有记录的子CCD辐照度距离倒数加权法插值获得。

优选的：

S801，计算各细分时间段中心时刻的TDI CCD离散平均焦面辐照度；

S802，计算一个积分时间段内的TDI CCD时间平均焦面辐照度；具体为：将通过步骤S801得到的TDI CCD的每个CCD在各离散时刻的TDI CCD离散平均焦面辐照度取平均，得到一个积分时间段内的TDI CCD时间平均焦面辐照度。

本发明还提供了一种计算机可读介质，所述介质中存储有计算机指令，所述指令执行前述任意一种方法。

本发明以30厘米分辨率的Worldview-3全色立体影像作为仿真数据源，对50-80厘米分辨率的立体测绘卫星影像进行成像仿真，在仿真链路中把平台参数、相机参数有机地串联起来，实现了几何仿真和辐射仿真的结合。

本发明的有益效果是：本发明实现了基于甚高分辨率卫星Worldview-3立体影像的分米级星载TDI CCD立体测绘相机成像仿真，模拟的两个相机的影像源于视角接近的输入源立体影像。最大程度地避免了同一数据源条件下地面几何模型过度简化引起的立体影像逼真度不高的问题。通过输入甚高分辨率卫星影像的像点三维坐标向TDI像面投射的方法计算TDI影像子CCD的观测点，可以实现较快速的视点检测，避免了传统迭代方法在高程突变位置不收敛的问题。对TDI CCD相机多级动态积分的过程进行严格模拟，实现了卫星轨道、姿态数据，相机几何、辐射参数，TDI CCD多级动态积分这一融合了严格几何转换和辐射转换过程的高精度模拟，保证了模拟结果的可信度。

附图说明

图1为根据本发明实施例的、分米级星载TDI CCD立体测绘相机成像仿真方法的流程图；

图2为TDI CCD与细分CCD示意图。

具体实施方式

下面参照附图对本发明进行更全面的描述，其中说明本发明的示例性实施例。

如图1所示，一种分米级星载TDI CCD立体测绘相机成像仿真方法，其特征在于，包括以下步骤：

步骤S1，获取基础数据，并对所述基础数据进行预处理，生成卫星立体影像每一点的影像-几何模型。

根据本发明的实施例，通过导入数据的方式获取基础数据，导入的基础数据为Worldview-3卫星立体影像和辅助RPC(Rational Polynomial Coefficient)参数。本发明采用Worldview-3卫星立体影像数据，是因为Worldview-3卫星立体影像是目前商业应用立体测绘卫星影像中分辨率最高的。其他分米级立体测绘卫星的分辨率一般0.5米-0.8米左右。因此可以将其他设计分辨率低于Worldview-3的分米级立体测绘卫星影像仿真作为Worldview-3影像的退化处理。

对基础数据进行预处理包括：

(S11)，通过多次前方交会和后方交会迭代计算，消除RPC参数的相对误差；

具体为：首先通过SIFT匹配获取Worldview-3全色立体影像间的9个以上(9、16或25，即3*3、4*4或5*5均匀分布)分布均匀的同名像点；通过同名像点前方交会获取物方地面点坐标；执行后方交会，所述后方交会为以地面点和对应像点作为控制点校正RPC参数的分子常数项和一次项。执行三次以上前方交会、后方交会迭代计算，消除RPC参数的相对误差。

(S12)，通过准核线立体影像制作和半全局匹配算法计算左右影像的视差图，根据所述视差图计算同名像点坐标，进而前方交会计算准核线影像每一点的物方大地坐标(L，B，H)，其中，L表示经度、B表示纬度、H表示椭球高。

步骤S12具体可以为：通过投影轨迹法计算Worldview-3准核线立体影像和各自RPC参数；通过半全局匹配算法计算左右影像的视差图，通过左右一致性检测剔除误匹配点并通过距离倒数加权法拟合误匹配点的视差；根据视差图计算同名像点坐标，前方交会计算所述准核线立体影像每一点的物方大地坐标(L，B，H)，其中，L表示经度、B表示纬度、H表示椭球高。

通过以上预处理操作，Worldview-3准核线影像上的每个像点不仅包括灰度信息，还包括用于跟仿真影像建立映射关系的几何位置信息。

在立体影像的仿真过程中，每个视角的仿真相机对应的数据为观测角度接近的Worldview-3准核线立体影像、RPC参数和准核线立体影像每一点的物方大地坐标。这是保证立体影像仿真精度、体现立体仿真影像间几何、辐射差异性的数据基础，因此，本发明采用立体影像上每一点的“影像-几何模型”(即，每一点的物方大地坐标(L，B，H)和影像灰度)作为模拟数据源，以提高仿真精度，这是本发明的改进点之一。

步骤S2，计算卫星轨道数据，所述卫星轨道数据包括地球固连坐标系下卫星的扫描时刻、位置、速度数据。

本发明通过步骤201-210所述的计算流程来计算卫星轨道数据，步骤201-210所述的计算流程不是严格的动力学计算方法，而是一种简化的数学模型，这样可以更有效的计算卫星轨道数据。

轨道数据计算的目的是建立卫星本体坐标系的原点(X_o，Y_o，Z_o)在卫星运行的每一时刻，在WGS84地球固连坐标系中的运行轨迹模型，卫星原点轨迹数据用卫星相机成像时间段内每一时刻卫星本体坐标系原点的位置速度矢量(X_o，Y_o，Z_o，V_Xo，V_Yo，V_Zo)来表示。模拟卫星轨道数据的频率是卫星下传轨道数据频率的10倍，模拟成像任一时刻的卫星位置、速度矢量由该时刻前后各两组位置、速度矢量均通过拉格朗日插值获得。

本发明用规则圆轨道来描述卫星轨道，卫星轨道数据的计算具体包括以下步骤：

步骤S201，计算惯性系轨道数据。

在惯性参考系的XOZ平面上，取一组卫星测量坐标系原点坐标，各原点坐标与惯性参考系原点的距离相同，均为R，R等于平均地球半径加上平均轨道高度。假设第一个原点对应时刻为0，相邻两个质心对应时刻相差0.1秒，由R可求得轨道周期T：

(单位：秒)

其中G为万有引力常数，等于6.67*10^-11米³/(千克*秒²)，M为地球质量，等于5.98*10²⁴千克。

降交点时刻Tdsd(秒)减去T/4就是第一个质心对应的实际时间T₀，即

T₀＝T_dsd-T/4

以下取样时刻依次加0.1秒为相应的实际时间。dθ＝360/(T*10)为相邻两个原点 的地心角。第一个质心在Z轴上(Z＝R)，第i(i＝0，1，…)个取样时刻的卫星原点位置为：

其中，θ＝(90-dθ×i)/180.0×π。

设轨道倾角为1，将取样坐标绕X轴顺时针旋转I-90度，然后绕Z轴逆时针旋转L度。

步骤S202，计算第i个时刻卫星原点(Xi，Yi，Zi)对应的初始大地坐标(Li，Bi，Hi)。步骤S202具体包括：

对第i个取样时刻，将卫星原点(Xi，Yi，Zi)(第i个取样时刻的三轴坐标)位置绕Z轴旋转i×dL度。dL＝(1/2400)°，其为0.1秒内地球自转角度：

其中，

计算第i个时刻卫星原点(Xi，Yi，Zi)对应的初始大地坐标(Li，Bi，Hi)，计算公式为：

L＝atan(Y/X)

H＝P/cos(B)-N

其中，(L，B，H)为(L_i，B_i，H_i)序列的任一组数据，(X，Y，Z)为对应(X_i，Y_i，Z_i)序列的数据，(Xi，Yi，Zi)为第i个时刻卫星原点三维坐标。其中：

a＝6378137

θ＝atan(Z/P*a/b)

步骤S203，计算经度L与目标区域中心经度L_c的偏差ΔL，步骤S203具体包括：

用最小二乘法得到经度和纬度的对应关系式：

L＝a+b*B+c*B²+d*B³

其中a，b，c，d为星下点纬度B与经度L间的最小二乘拟合系数。

将模拟目标区域中心纬度Bc带入这个关系式，求得的经度L与目标区域中心经度L_c的偏差ΔL：

ΔL＝L_c-L

步骤S204，卫星原点大地坐标的经度L_i(i＝0，1，…)加上ΔL，得到卫星原点的星下点通过目标区域中心点(L_c，B_c)的卫星原点大地坐标序列(L′_i，B_i，H_i)((i＝0，1，…)，L′_i＝L_i+ΔL)。

步骤S205，将卫星原点大地坐标序列(L′_i，B_i，H_i)转换成地心直角坐标序列(X′_i，Y′_i，Z′_i)(i＝0，1，…)，计算公式为：

X＝(N+H)*cos(B)*cos(L)

Y＝(N+H)*cos(B)*sin(L)

Z＝(N*(1-e²)+H)*sin(B)

其中，(L，B，H)为(L′_i，B_i，H_i)序列的任一组数据，(X，Y，Z)为对应(X′_i，Y′_i，Z′_i)序 列的数据。

a＝6378137

步骤S206，返回步骤201，利用公式(1)计算卫星原点(Xi，Yi，Zi)时将角度θ加上0.00002度。

步骤S207，重复步骤S202至S206，得到新的卫星原点地心直角坐标序列(X″_i，Y″_i，Z″_i)。

步骤S208，计算卫星原点速度序列(V_xi，V_yi，V_zi)：

V_xi＝(X″_i-X′_i)/(T*0.00002/360)

V_yi＝(Y″_i-Y′_i)/(T*0.00002/360)

V_zi＝(Z″_i-Z′_i)/(T*0.00002/360)

其中第i组卫星原点位置速度对应时间为：T_i＝T₀+i*0.1。

步骤S209，保存卫星原点轨迹数据(T_i，X′_i，Y′_i，Z′_i，V_xi，V_yi，V_zi)。

即，生成了时刻Ti对应的卫星本体坐标系原点的位置、速度数据。

步骤S3，计算卫星姿态数据。

通过功率谱反推方法计算卫星的三轴高频姿态数据，与三轴指向精度数据相加得到卫星本体坐标系相对局部轨道坐标系的三轴姿态数据。

姿态建模的目的是获取符合卫星姿态设计指标的成像每一时刻卫星本体坐标系相对局部轨道坐标系的姿态角(yaw，pitch，roll)，其中yaw为偏航角，pitch为俯仰角，roll为滚动角。卫星的主要姿态设计指标包括：

三轴姿态指向精度(yaw₀，pitch₀，roll₀)，度(3σ)；

三轴姿态稳定度(yaw₁，pitch₁，roll₁)，度/秒(3σ)；

本发明分两步获取卫星姿态参数：首先生成符合卫星设计姿态指向精度的系统指向误差，然后，用功率谱分析法获取符合卫星设计姿态稳定度的姿态参数，将指向误差参数和姿态参数相加得到最终的姿态参数。

步骤S3具体包括：

步骤301，计算指向误差参数：

根据三轴姿态指向精度(yaw₀，pitch₀，roll₀)生成三轴指向角误差参数：

yaw₀′＝1.0·rand()/RAND_MAX·yaw₀

pitch₀′＝1.0·rand()/RAND_MAX·pitch₀

roll₀′＝1.0·rand()/RAND_MAX·roll₀

其中，rand()为C++的均匀分布随机整数生成器，RAND_MAX为最大随机整数。

步骤302，利用功率谱分析法获取姿态参数：

步骤S3021，根据模拟设置的多个频段中心频率v、频段宽度w及该频段的幅值A(w_i，v_i，A_i)(i表示第i组模拟频段参数)生成滤波器：

Q(v)为频率v处的滤波系数，n为设置的频段个数。

步骤S3022，生成随机振动功率谱

其中，F(v)为频率v处的功率谱数据。rand()为随机正整数，RAND_MAX为随机正整数的最大范围。

步骤S3023，对随机振动功率谱数据进行频率域滤波得到姿态振动功率谱：

F(v)′＝Q(v)*F(v)

步骤S3024，反傅立叶变换变回时间域得一个姿态角的姿态数据序列：

f(t_i)＝IFFT(F(v)′)

其中，i表示第i个姿态角，t_i表示第i个姿态角的时间，设频率域的最高频率为v_max，则t_i＝1.0/v_max*i。

步骤S303，将通过步骤S301得到的指向误差参数和通过步骤S302得到的姿态参数相加得到最终的姿态参数。

步骤S303具体包括：

每隔1秒统计生成的姿态角序列f(t_i)的均方差σ，将姿态角序列的每个姿态角乘以σ/(3yaw₁)，然后加上偏航指向角系统量yaw′₀，得到偏航姿态角序列：

yaw(t_i)＝f(t_i)*σ/(3*yaw₁)+yaw′₀

将姿态角序列的每个姿态角乘以σ/(3pitch₁)，然后加上俯仰指向角系统量pitch′₀，得到俯仰姿态角序列：

pitch(t_i)＝f(t_i)*σ/(3*pitch₁)+pitch′₀

将姿态角序列的每个姿态角乘以σ/(3roll₁)，然后加上滚动指向角系统量roll′₀，得到滚动姿态角序列：

roll(t_i)＝f(t_i)*σ/(3*roll₁)+roll′₀。

姿态数据是决定影像立体几何精度的关键数据。由于输入Pleiades数据范围远小于仿真影像的数据范围，最终的仿真影像只会覆盖整景仿真影像的一部分区域和高度。为了更好的评估整景影像的几何精度，需要进行多次实验，使有效仿真影像在任意姿态区间的立体几何精度。要实现这个目标需要每次仿真时在上述已获取姿态的起始时间上加入一个随机变量，本发明中将该随机变量的范围设为仿真姿态对应时间跨度的三分之一。

步骤S4，相机几何建模，其包括两个部分：第一部分是计算相机安装模型，所述相机安装模型计算所需的基本参数包括相机光学节点在卫星本体坐标系的安装位置(X_s，Y_s，Z_s)和相机测量坐标系相对卫星本体坐标系的安装角(yaw_s，pitch_s，roll_s)；第二部分是计算TDI CCD畸变模型。

步骤S4具体包括以下步骤：

步骤S401，计算相机安装模型。

相机安装模型计算所需的基本参数包括相机光学节点在卫星本体坐标系的安装位置(X_s，Y_s，Z_s)和相机测量坐标系相对卫星本体坐标系的安装角(yaw_s，pitch_s，roll_s)，由于相机安装的过程中会存在一定的偏差，实际的安装角并不等于理想值，本发明使用的相机安装角为理想安装值。步骤S401具体包括以下子步骤：

S4011，将相机光学节点位置结合成像时刻转换到物方坐标系：

S4012，根据相机安装角计算相机测量坐标系相对卫星本体坐标系的旋转矩阵r_set；

其中，相机测量坐标系相对卫星本体坐标系的旋转矩阵r_set为：

其中，

S4013，计算相机光学节点在卫星本体坐标系的安装向量U在WGS84地球固连坐标系的向量U′：

其中，设相机光学节点在卫星本体坐标系的安装位置为(X_s，Y_s，Z_s)，则U为卫星本体坐标系原点到(X_s，Y_s，Z_s)的向量，为卫星本体坐标系相对t时刻的局部轨道坐标系的旋转矩阵，为卫星本体坐标系原点与WGS84地球固连坐标系对应旋转矩阵。

下面作进一步说明：

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

上面已求得相机测量坐标系相对卫星本体坐标系的旋转矩阵为t_set。卫星本体坐标系相对t时刻的局部轨道坐标系的姿态角为(yaw，pitch，roll)，对应旋转矩阵为

其中，

卫星本体坐标系原点在WGS84地球固连坐标系的位置速度矢量为对应旋转矩阵为

为卫星位置归一化向量；

为卫星速度归一化向量；

[C_x，C_y，C_z]^T为和的叉积；

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

S4014，得到相机光学节点的运行位置为：

其中，为卫星位置归一化向量，为相机光学节点在卫星本体坐标系的安装位置向量在地球固连坐标系的形式。

S4015，设相机成像的积分时间值为Δt₀，t+Δt₀时刻的相机光学节点的运行位置为t时刻相机光学节点的速度矢量为：

该速度矢量在地面的投影为：

其中H_obt为轨道高度。

设CCD尺寸为a*a，a单位为米，相机主距为f，则地面分辨率GR＝a/f*Hobt。

实际积分时间为：

重新计算t时刻相机光学节点的速度矢量为：

将所有时间节点的速度矢量投影也保存下来用于计算仿真影像起算时间。

卫星成像时间区间内一定时间间隔(卫星下传轨道数据时间间隔的十分之一)的相机光学节点“时间-位置-速度”数据是仿真流程中光线向量计算的基础。

S402计算TDI CCD畸变模型。

每片TDI CCD的安装误差用一个以CCD列坐标为变量的二阶多项式表示，常数项表示平移，一阶项表示整体旋转，二阶项表示弯曲。根据相机光学系统的畸变模型将不同视场的畸变量换算为像面位置误差，与安装误差多项式的常数项相加。计算每个(用i表示CCD列号)CCD在垂直TDI积分方向的畸变量σ_CSi和每个CCD在沿TDI积分方向的畸变量σ_ASi。两个CCD列之间的亚像元点(列坐标为c，i≤c≤i+1)的畸变量为：

σ_AS(c)＝σ_ASi-0.5＜c-i≤0.5

σ_CS(c)＝σ_CS(i+1)0.5＜c-i≤1

设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焦面上任意一点的主光线对应的地面点。

通过以上两个几何模型的建立，获取了考虑畸变的像空间光线向量计算的数据基础和光学节点在物方坐标系的位置、速度向量形式。也就是说，加入畸变误差的像点，转到物方坐标系的光学节点和对应地面点符合摄影测量里面的共线方程。

需要指出的是，前面的4个步骤是为后面的仿真流程做准备，其执行顺序没有严格要求，仿真的主要工作就是建立仿真链路各个环节的数学模型供仿真链路调用。如第1步建立的地面模型，第2、3步是生成卫星的轨迹和姿态数据，第4步是生成后面光线向量计算所需的光学节点轨迹、旋转矩阵等数据。后面的仿真流程以辐射信息的获取和转换为主线，前面的四个步骤所做的几何数据准备是辐射信息获取的基础。

根据本发明的一优选实施例，在步骤4之后还可以包括计算相机的起始成像时间的步骤。

立体成像仿真就是模拟两个不同观测角度的影像。这一步骤的目的是确定立体测绘相机的起始成像时间，同轨立体成像模式下立体成像分两种方法：一是单相机沿轨道侧摆一定角度构建立体；二是两个以上不同安装角的相机分别成像构建立体。前者只需要在步骤3姿态建模的基础上对俯仰角pitch加上一个固定的观测角度Pset，作为姿态角的一个构成部分，而姿态角的数据形式不变。两者确定相机起始成像时间的方法相同，具体步骤如下：

(1)：将相机起始成像时间设置为一个步骤S4获取相机光学节点数据序列总的时间区间的任一时刻ts。线性插值获取该时刻的相机安装角、卫星平台三轴姿态角、相机主点位置P[X，Y，Z]、速度向量V[vx，vy，vz]，计算相机安装矩阵Mset、姿态矩阵Matt、轨道-物方转换矩阵Mobt，将像空间向量u0＝[0，0，-f]转换为物方向量u。

u＝Mobt*(Matt*(Mset*u0))。

(2)：计算向量u与地球椭球的交点大地坐标(L，B，0)，并转换为地心直角坐标P(X，Y，Z)。

(3)：将目标区域中心大地坐标Pc(Lc，Bc，0)转换为地心直角坐标Pc1(Xc，Yc，Zc)。 计算P和Pc1间的距离D。根据该时刻相机光学节点在地面上的投影计算时间差：

(4)：将起始时间替换为：

(5)：重复步骤(1)到(4)，直至D小于一定门限(本发明设为0.1)。

(6)：设模拟影像行数为height，则起始成像时刻为：

ts＝ts-height/2*Δt。

其中Δt为步骤4计算的积分时间。

步骤S5，利用步骤S1-S4计算的数据，计算TDI CCD瞬时焦面能量影像的数据源范围。

首先将某一时刻的TDI CCD成像面细分为m*m个更小的子CCD，这时整个TDI CCD像面可以看作一幅分辨率为模拟影像分辨率m倍的框幅式影像。如图2所示，其中，21表示第i级CCD，22表示第j列CCD，23表示细分子CCD中心。

设TDI CCD积分级数为M，TDI CCD单行探元数为N，则TDI CCD像面构成的框幅式影像(以下称为TDI影像)行数为M*m，列数为N*m。

设TDI影像四角像面坐标为[xi，yi](i＝1，2，3，4)，相机主距为f。则TDI影像四角与光学系统主点构成向量u0＝[xi，yi，-f](i＝1，2，3，4)。

根据相机安装角、卫星平台三轴姿态角、相机主点位置P[X，Y，Z]、速度向量V[vx，vy，vz]计算的相机安装矩阵Mset、姿态矩阵Matt、轨道-物方转换矩阵Mobt将像空间向量u0转换为物方向量u。

u＝Mobt*(Matt*(Mset*u0))

设u对应地面点高程为H，根据P、H、u可以计算地面点平面坐标(L，B)。将(L，B，H)通过Worldview-3对应视角的准核线影像RPC计算像面坐标，并获取该点的高程Htrue。

将H替换为Htrue，重复前面的计算过程，直到|H-Htrue|小于一定的限差dh(如0.1m)。

这样就得到TDI影像四角坐标在Worldview-3对应视角的准核线影像上的四个坐标，也就确定了所需数据的概略范围。将该范围在TDI积分方向对应的Worldview-3准核线影像方向外扩一定距离确定参与计算的数据范围。

步骤S6，计算滤波前的TDI CCD瞬时焦面辐照度影像。

将步骤S5确定的数据范围内的Worldview-3准核线影像像点对应地面点通过TDI影像的瞬时外方位元素投射到TDI影像上，加上相应的TDI CCD畸变改正量后，按子CCD大小取整后查找坐标对应子CCD的行、列号，每个CCD都有一个存储单元记录投射地面元在Worldview-3准核线影像的坐标和地面元与TDI影像投影中心的距离。所有地面点都完成投射后，有投射记录的子CCD将挑选出距离最近的地面元作为可见地面元，并将对应的Worldview-3准核线影像DN值转换为像面辐照度记录下来。没有投射记录的子CCD辐照度通过周围有记录的子CCD辐照度距离倒数加权法插值获得。

步骤S7，用点扩散函数对TDI CCD瞬时焦面辐照度影像进行滤波，得到一个时刻的TDI CCD离散平均焦面辐照度。

以高斯点扩散函数模型作为光学点扩散函数模型，高斯点扩散函数为一个圆对称的二维函数，形式为：

PSF(x，y)＝exp(-(x²+y²)/(2σ²))

其中，exp为指数函数，σ为函数下降参数，表示点扩散函数值随(x，y)变化的速度，(x，y)为点扩散函数作用范围内(PSF(x，y)＞0)的一点相对点扩散函数原点的二维平面坐标。输入信号和输出信号都是连续信号，在模拟算法中通过获取CCD成像面上高密度的离散能量信号来近似代替连续的能量信号分布。设CCD成像面的离散采样间隔为aμm，将点扩散函数也进行离散化，这时的平面坐标变为整数型的离散序号。用离散化后的点扩散函数与离散正弦信号卷积，统计卷积前后的信号调制度计算MTF。计算一个σ序列(σ₁，σ₂，…，σ_n)对应的MTF序列(MTF₁，MTF₂，…，MTF_n)，模拟使用的MTF为MTFu(MTF_m＜MTF_u＜MTF_m+1)(n＞m≥1)，则MTFu对应σ_u通过σ_m和σ_m+1线性插值获得。最后通过σ_u计算离散PSF。

静态点扩散函数除了光学点扩散函数外，还包括CCD的矩形滤波作用，CCD的矩形滤波是一个将单个CCD焦面内的辐照度(或等效于CCD焦面内每一点对应入瞳辐亮度)取平均的过程。这里的辐照度(或入瞳辐亮度)是经光学点扩散函数滤波后的辐照度(或入瞳辐亮度)。

经过光学点扩散函数滤波和CCD矩形滤波后得到一个时刻的TDI CCD离散平均焦面辐照度。

步骤S8，计算一个积分时间区间内多个细分时刻的TDI辐照度影像，累加后取平均，得到一个积分区间内的TDI时间平均辐照度影像。

步骤S8具体包括以下步骤：

S801计算各细分时间段中心时刻的TDI CCD离散平均焦面辐照度。

步骤S7获取的是相机TDI CCD在一个时刻的TDI CCD离散平均焦面辐照度，而每一级CCD在一个积分时间段内积累的电荷是该积分时间段内CCD接收平均辐照度(或入瞳辐亮度)的函数。因此，需要获取一个积分时间段内CCD接收平均辐照度(或入瞳辐亮度)。这是一个在时间上连续的过程，

本发明将积分时间段细分(一般均匀细分为5个时间段)并获取每个细分时间段中心时刻的TDI CCD离散平均焦面辐照度。

S802计算一个积分时间段内的TDI CCD时间平均焦面辐照度

得到一个积分时间段内各离散时刻(即，每个细分时间段中心时刻)的TDI CCD离散平均焦面辐照度后，将TDI CCD的每个CCD在各离散时刻的TDI CCD离散平均焦面辐照度取平均，得到一个积分时间段内的TDI CCD时间平均焦面辐照度。

步骤S9，将一个积分区间内的TDI时间平均辐照度影像转换为光生电荷数影像，并根据霰粒噪声的生成机制加入噪声电荷，生成该积分区间的电荷数影像。

步骤S10，重复步骤S8和S9，计算下一个积分时间区间的电荷数影像，完成M级电荷影像累加后，进行模数转换，输出一行灰度值影像。

具体的，设积分级数为M，依次计算下M-1个积分时间区间的光生电荷数影像，前后两幅光生电荷数影像错开一行CCD累加，M级积分计算后，各级CCD的累加次数分别为M-1，M-2，……，0，这时，第一级CCD完成电荷积累过程。将辐射分辨率参数转换为以光子数度量的形式，以此为参量将总的电荷数进行截断处理并通过模数转换公式转换为影像DN值，以一行影像的形式输出。

步骤S11，成像时刻增加积分时间Δt，重复步骤S5-S10，计算一幅模拟影像的每一行灰度值影像，最终输出一幅模拟影像。其中，Δt通过步骤401计算。

具体的，S11在S10的基础上进行下一积分时间区间的电荷数影像，与之前的累加后又一行完成M级电荷数累加，模数转换输出又一行灰度值影像，以此类推直至输出一景影像的所有行。

通过上述说明，需要再次指出的是：

本发明以30厘米分辨率的Worldview-3全色立体影像作为仿真数据源，对50-80厘米分辨率的立体测绘卫星影像进行成像仿真，在仿真链路中把平台参数、相机参数有机地串联起来，实现了几何仿真和辐射仿真的结合。

本发明的有益效果是：本发明实现了基于甚高分辨率卫星Worldview-3立体影像的分米级星载TDI CCD立体测绘相机成像仿真，模拟的两个相机的影像源于视角接近的输入源立体影像。

另外，本发明还提出一种计算机可读介质，所述介质中存储有计算机指令，所述指令执行上文所述的方法。

另外，本发明还提出一种分米级星载TDI CCD立体测绘相机成像仿真系统，该系统中可以设置相应的功能模块来执行上文所述的方法，如本领域技术人员可以理解的，该系统的功能模块可以由硬件或者软件来实现。

在本说明书的描述中，参考术语“实施例一”、“实施例二”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体方法、装置或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、方法、装置或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种分米级星载TDI CCD立体测绘相机成像仿真方法，其特征在于，包括以下步骤：

步骤S1，获取基础数据，并对所述基础数据进行预处理，生成卫星立体影像每一点的影像-几何模型；其中，所述基础数据为Worldview-3卫星立体影像和辅助RPC参数，所述预处理具体包括：

(S11)，通过多次前方交会和后方交会迭代计算，消除RPC参数的相对误差；

(S12)，通过准核线立体影像制作和半全局匹配算法计算左右影像的视差图，根据所述视差图计算同名像点坐标，进而前方交会计算准核线影像每一点的物方大地坐标(L，B，H)，其中，L表示经度、B表示纬度、H表示椭球高；

步骤S2，计算卫星轨道数据，所述卫星轨道数据包括地球固连坐标系下卫星的扫描时刻、位置、速度数据；

步骤S3，计算卫星姿态数据；

步骤S4，相机几何建模，其具体包括

步骤S401，计算相机安装模型；

步骤S402，计算TDI CCD畸变模型；

步骤S5，利用步骤S1-S4计算的数据，计算TDI CCD瞬时焦面能量影像的数据源范围；

步骤S6，计算滤波前的TDI CCD瞬时焦面辐照度影像；

步骤S7，用点扩散函数对TDI CCD瞬时焦面辐照度影像进行滤波，得到一个时刻的TDICCD离散平均焦面辐照度；

步骤S8，计算一个积分时间区间内多个细分时刻的TDI辐照度影像，累加后取平均，得到一个积分区间内的TDI时间平均辐照度影像；

步骤S9，将一个积分区间内的TDI时间平均辐照度影像转换为光生电荷数影像，并根据霰粒噪声的生成机制加入噪声电荷，生成该积分区间的电荷数影像；

步骤S10，重复步骤S8和S9，计算下一个积分时间区间的电荷数影像，完成M级电荷数影像累加后，进行模数转换，输出一行灰度值影像；

步骤S11，计算一幅模拟影像的每一行灰度值影像，最终输出一幅模拟影像；

其中，步骤S2具体包括以下步骤：

步骤S201，计算惯性系轨道数据；

在惯性参考系的XOZ平面上，取一组卫星测量坐标系原点坐标，各原点坐标与惯性参考系原点的距离相同，均为R，R等于平均地球半径加上平均轨道高度；假设第一个原点对应时刻为0，相邻两个质心对应时刻相差0.1秒，由R可求得轨道周期T：

单位：秒，

其中G为万有引力常数，等于6.67*10^-11米³/(千克*秒²)，M为地球质量，等于5.98*10²⁴千克，

降交点时刻Tdsd秒减去T/4就是第一个质心对应的实际时间T₀，即

T₀＝T_dsd-T/4

以下取样时刻依次加0.1秒为相应的实际时间，dθ＝360/(T*10)为相邻两个原点的地心角，第一个质心在Z轴上，Z＝R，第i个取样时刻的卫星原点位置为：

X＝R×cos(θ)

Y＝0

Z＝R×sin(θ) (1)

其中，θ＝(90-dθ×i)/180.0×π,其中，i＝0,1,…；

设轨道倾角为I，将取样坐标绕X轴顺时针旋转I-90度，然后绕Z轴逆时针旋转L度；

步骤S202，计算第i个时刻卫星原点(Xi，Yi，Zi)对应的初始大地坐标(Li，Bi，Hi)；

步骤S203，计算经度L与目标区域中心经度L_c的偏差ΔL；

步骤S204，卫星原点大地坐标的经度L_i(i＝0,1,…)加上ΔL，得到卫星原点的星下点通过目标区域中心点(L_c,B_c)的卫星原点大地坐标序列(L′_i,B′_i,H′_i),其中，i＝0,1,…,L′_i＝L_i+ΔL；

步骤S205，将卫星原点大地坐标序列(L′_i,B_i,H_i)转换成地心直角坐标序列(X′_i,Y′_i,Z′_i)；

步骤S206，返回步骤201，利用公式(1)计算卫星原点(Xi，Yi，Zi)时将角度θ加上0.00002度；

步骤S207，重复步骤S202至S206，得到新的卫星原点地心直角坐标序列(X”_i,Y”_i,Z”_i)；

步骤S208，计算卫星原点速度序列(V_xi,V_yi,V_zi)：

V_xi＝(X”_i-X′_i)/(T*0.00002/360)

V_yi＝(Y”_i-Y′_i)/(T*0.00002/360)

V_zi＝(Z”_i-Z′_i)/(T*0.00002/360)

其中第i组卫星原点位置速度对应时间为：T_i＝T₀+i*0.1；

步骤S209，保存卫星原点轨迹数据(T_i,X′_i,Y′_i,Z′_i,V_xi,V_yi,V_zi)。

2.根据权利要求1所述的方法，其特征在于，步骤S1中的子步骤(S11)具体包括：

首先通过SIFT匹配获取Worldview-3全色立体影像间的9个以上分布均匀的同名像点；通过同名像点前方交会获取物方地面点坐标；执行后方交会，所述后方交会为以地面点和对应像点作为控制点校正RPC参数的分子常数项和一次项，执行三次以上前方交会、后方交会迭代计算，消除RPC参数的相对误差。

3.根据权利要求1所述的方法，其特征在于，步骤S3具体包括以下步骤：

步骤301，计算指向误差参数；

步骤302，利用功率谱分析法获取姿态参数；

步骤S303，将通过步骤S301得到的指向误差参数和通过步骤S302得到的姿态参数相加得到最终的姿态参数；

其中，步骤S303具体包括：

每隔1秒统计生成的姿态角序列f(t_i)的均方差σ，将姿态角序列的每个姿态角乘以σ/(3yaw₁)，然后加上偏航指向角系统量yaw′₀，得到偏航姿态角序列：

yaw(t_i)＝f(t_i)*σ/(3*yaw₁)+yaw′₀

将姿态角序列的每个姿态角乘以σ/(3pitch₁)，然后加上俯仰指向角系统量pitch′₀，得到俯仰姿态角序列：

pitch(t_i)＝f(t_i)*σ/(3*pitch₁)+pitch′₀

将姿态角序列的每个姿态角乘以σ/(3roll₁)，然后加上滚动指向角系统量roll′₀，得到滚动姿态角序列：

roll(t_i)＝f(t_i)*σ/(3*roll₁)+roll′₀，其中，(yaw₁,pitch₁,roll₁)表示三轴姿态稳定度。

4.根据权利要求1所述的方法，其特征在于：

步骤S401具体包括以下子步骤：

S4011，将相机光学节点位置结合成像时刻转换到物方坐标系：

S4012，根据相机安装角计算相机测量坐标系相对卫星本体坐标系的旋转矩阵r_set；

其中，相机测量坐标系相对卫星本体坐标系的旋转矩阵r_set为：

其中，

其中，(yaw_s,pitch_s,roll_s)表示相机测量坐标系相对卫星本体坐标系的安装角；

S4013，计算相机光学节点在卫星本体坐标系的安装向量U在WGS84地球固连坐标系的向量U′：

其中，设相机光学节点在卫星本体坐标系的安装位置为(X_s,Y_s,Z_s)，则U为卫星本体坐标系原点到(X_s,Y_s,Z_s)的向量，为卫星本体坐标系相对t时刻的局部轨道坐标系的旋转矩阵，为卫星本体坐标系原点与WGS84地球固连坐标系对应旋转矩阵；

S4014，得到相机光学节点的运行位置为：

其中，为卫星位置归一化向量，为相机光学节点在卫星本体坐标系的安装位置向量在地球固连坐标系的形式；

S4015，设相机成像的积分时间值为Δt₀，t+Δt₀时刻的相机光学节点的运行位置为t时刻相机光学节点的速度矢量为：

该速度矢量在地面的投影为：

其中H_obt为轨道高度；

设CCD尺寸为a*a，a单位为米，相机主距为f，则地面分辨率GR＝a/f*H_obt；

实际积分时间为：

重新计算t时刻相机光学节点的速度矢量为：

5.根据权利要求1所述的方法，其特征在于，步骤S402具体包括：

每片TDI CCD的安装误差用一个以CCD列坐标为变量的二阶多项式表示，根据相机光学系统的畸变模型将不同视场的畸变量换算为像面位置误差，与安装误差多项式的常数项相加，计算每个CCD在垂直TDI积分方向的畸变量σ_CSi和每个CCD在沿TDI积分方向的畸变量σ_ASi，两个CCD列之间的亚像元点的畸变量为：

σ_AS(c)＝σ_ASi-0.5＜c-i≤0.5

σ_CS(c)＝σ_CS(i+1)0.5＜c-i≤1

其中，i表示CCD序列号，C表示列坐标，i≤c≤i+1。

6.根据权利要求1所述的方法，其特征在于，步骤S5具体包括以下步骤：

首先将某一时刻的TDI CCD成像面细分为m*m个更小的子CCD，这时整个TDI CCD像面可以看作一幅分辨率为模拟影像分辨率m倍的框幅式影像；设TDI CCD积分级数为M，TDI CCD单行探元数为N，则TDI CCD像面构成的框幅式影像行数为M*m，列数为N*m，以下将TDI CCD像面构成的框幅式影像称为TDI影像；

设TDI影像四角像面坐标为[xi,yi](i＝1，2，3，4)，相机主距为f，则TDI影像四角与光学系统主点构成向量u0＝[xi,yi,-f](i＝1，2，3，4)；

根据相机安装角、卫星平台三轴姿态角、相机主点位置P[X,Y,Z]、速度向量V[vx,vy,vz]计算的相机安装矩阵Mset、姿态矩阵Matt、轨道-物方转换矩阵Mobt将像空间向量u0转换为物方向量u；

设u对应地面点高程为H，根据P、H、u可以计算地面点平面坐标(L,B)，将(L,B,H)通过Worldview-3对应视角的准核线影像RPC计算像面坐标，并获取该点的高程Htrue；

将H替换为Htrue，重复前面的计算过程，直到|H-Htrue|小于一定的限差dh，所述限差为0.1m；

这样就得到TDI影像四角坐标在Worldview-3对应视角的准核线影像上的四个坐标，也就确定了所需数据的概略范围，将该范围在TDI积分方向对应的Worldview-3准核线影像方向外扩一定距离确定参与计算的数据范围。

7.根据权利要求1所述的方法，其特征在于，步骤S6具体包括：

将将步骤S5确定的数据范围内的Worldview-3准核线影像像点对应地面点通过TDI影像的瞬时外方位元素投射到TDI影像上，加上相应的TDI CCD畸变改正量后，按子CCD大小取整后查找坐标对应子CCD的行、列号，每个CCD都有一个存储单元记录投射地面元在Worldview-3准核线影像的坐标和地面元与TDI影像投影中心的距离，所有地面点都完成投射后，有投射记录的子CCD将挑选出距离最近的地面元作为可见地面元，并将对应的Worldview-3准核线影像DN值转换为像面辐照度记录下来，没有投射记录的子CCD辐照度通过周围有记录的子CCD辐照度距离倒数加权法插值获得。

8.根据权利要求1所述的方法，其特征在于，步骤S8具体包括以下步骤：

S801，计算各细分时间段中心时刻的TDI CCD离散平均焦面辐照度；

S802，计算一个积分时间段内的TDI CCD时间平均焦面辐照度；具体为：将通过步骤S801得到的TDI CCD的每个CCD在各离散时刻的TDI CCD离散平均焦面辐照度取平均，得到一个积分时间段内的TDI CCD时间平均焦面辐照度。

9.一种计算机可读介质，所述介质中存储有计算机指令，所述指令执行根据权利要求1所述的方法。