一种分米级星载tdi ccd立体测绘相机成像仿真方法和系统的制作方法
【技术领域】
[0001] 本发明设及卫星成像仿真技术领域，特别设及一种分米级星载TDI CCD立体测绘 相机成像仿真方法和系统。
【背景技术】
[0002] 随着我国首颗民用传输型立体测绘卫星资源=号的发射成功，结束了我国W往完 全依靠国外高分辨率影像数据的历史。大量试验证明，资源=号卫星无论是分辨能力还是 几何精度均超过了国外同类卫星。资源=号卫星数据接收W后，地面处理能力也迅速跟进， 保障了数据的有效利用。资源=号卫星的成功很大程度上得益于先前对卫星各项技术指标 的详细论证，而卫星成像仿真在资源=号的技术指标论证和地面应用系统建设中都发挥了 重要作用。卫星成像仿真已经成为现代卫星工程中的一个重要环节。
[0003] 资源=号卫星正视相机分辨率最高为2.1米，前后视相机分辨率为3.6米，主要面 向1:5万立体测图，对于更大比例尺的测绘应用还不能满足需求。因此，面向1:1万立体测图 的分米级立体测绘卫星也很快被提上议事日程，因此需要研究分米级立体测绘卫星相机的 成像仿真技术。
[0004] 在可见、近红外光学成像卫星的成像仿真技术领域，国内外学者提出了物理仿真、 半物理仿真和计算机数值仿真=种技术途径，其中计算机数值仿真的工程成本最低。采用 各子系统的实验室标定数据结合设计数据建立各子系统的数学模型，也就是用数学公式描 述信号传输和转换的物理过程，简称子系统建模。在子系统建模的基础上实现在轨动态成 像过程的全链路的仿真流程，能比较全面地模拟卫星的设计技术指标对卫星在轨成像质量 的影响，准确地模拟出延迟积分成像高分辨率光学卫星在轨成像过程中几何和福射质量的 下降与影像总体几何、福射精度的对应关系，建立起卫星设计参数与在轨成像质量之间的 直接关联。
[0005] 真实的卫星相机成像是对真实的地面物理场景进行福射信息采集和转化的过程， 进行计算机数值仿真同样需要一个虚拟的数字物理场景，也就是一个数值化的地面模型， 运个地面模型需要包括地面每一点的几何信息和福射信息。换句话说，地面模型由几何模 型和福射模型构成，几何模型提供每一地面点的几何信息，福射模型提供地面每一点的福 射信息。
[0006] 现有的仿真软件和方法大都局限于卫星影像的福射质量和几何分辨能力，不能解 决分米级TDI CCD立体测绘相机最关注的立体几何精度仿真等问题，因此国产分米级立体 测绘卫星的成像仿真研究十分迫切。卫星影像的分辨能力是采样间隔、MTF、噪声、平台运 动、工作环境等一系列因素综合作用的结果，要达到分米级分辨率和1:10000测绘精度仅仅 实现CCD对应的地面采样大小达到分米级是远远不够的。需要综合考虑各种影响因素，设计 合理的技术指标和工作模式才能实现。另外，卫星成像过程不是简单的几何位置转换或福 射信息转换过程，而是几何、福射因素共同作用的结果，不能将二者分离开来。立体成像最 显著的特点是不同视角获取的影像之间有一定的差异性，体现在几何和福射两个方面。采 用单一的地面几何模型和福射模型很难对运两个方面进行高精度建模。而视角接近的其他 甚高分辨率立体卫星影像本身就包含了运种不同视角间几何信息和福射信息的差异性，因 此可W作为立体成像仿真的数据源。
[0007] 本发明所要解决的技术问题是:如何提高卫星成像仿真的模拟精度。
【发明内容】
[0008] 为此，本发明提出一种分米级星载TDI CCD立体测绘相机成像仿真方法和系统，可 充分地消除由于现有技术的限制和缺陷导致的一个或多个问题。
[0009] 本发明另外的优点、目的和特性，一部分将在下面的说明书中得到阐明，而另一部 分对于本领域的普通技术人员通过对下面的说明的考察将是明显的或从本发明的实施中 学到。通过在文字的说明书和权利要求书及附图中特别地指出的结构可实现和获得本发明 目的和优点。
[0010] 本发明提供了一种分米级星载TDI CCD立体测绘相机成像仿真方法，其特征在于， 包括W下步骤：
[0011] 步骤SI,获取基础数据，并对所述基础数据进行预处理，生成卫星立体影像每一点 的影像-几何模型；其中，所述基础数据为Worldview-3卫星立体影像和辅助RPC参数，所述 预处理具体包括：
[0012] (Sl 1 )，通过多次前方交会和后方交会迭代计算，消除RPC参数的相对误差；
[0013] (S12)，通过准核线立体影像制作和半全局匹配算法计算左右影像的视差图，根据 所述视差图计算同名像点坐标，进而前方交会计算准核线影像每一点的物方大地坐标化， B，H)，其中，L表示经度、B表示缔度、H表示楠球高。
[0014] 步骤S2,计算卫星轨道数据，所述卫星轨道数据包括地球固连坐标系下卫星的扫 描时刻、位置、速度数据；
[0015] 步骤S3，计算卫星姿态数据；
[0016] 步骤S4，相机几何建模，其具体包括
[0017] 步骤S401，计算相机安装模型；
[0018] 步骤S402,计算TDI CCD崎变模型；
[0019] 步骤S5,利用步骤S1-S4计算的数据，计算TDI CCD瞬时焦面能量影像的数据源范 围；
[0020] 步骤S6,计算滤波前的TDI CO)瞬时焦面福照度影像；
[0021] 步骤S7,用点扩散函数对TDI CO)瞬时焦面福照度影像进行滤波，得到一个时刻的 TDI CCD离散平均焦面福照度；
[0022] 步骤S8,计算一个积分时间区间内多个细分时刻的TDI福照度影像，累加后取平 均，得到一个积分区间内的TDI时间平均福照度影像；
[0023] 步骤S9,将一个积分区间内的TDI时间平均福照度影像转换为光生电荷数影像，并 根据靈粒噪声的生成机制加入噪声电荷，生成该积分区间的电荷数影像；
[0024] 步骤S10,重复步骤S8和S9,计算下一个积分时间区间的电荷数影像，完成M级电荷 数影像累加后，进行模数转换，输出一行灰度值影像；
[0025] 步骤Sll，计算一幅模拟影像的每一行灰度值影像，最终输出一幅模拟影像。
[00%]优选的，步骤SI中的子步骤(Sll)具体包括：
[0027] 首先通过SIFT匹配获取Worldview-3全色立体影像间的9个W上（例如，可W为9、 16或25个)分布均匀的同名像点;通过同名像点前方交会获取物方地面点坐标;执行后方交 会，所述后方交会为W地面点和对应像点作为控制点校正RPC参数的分子常数项和一次项， 执行=次W上前方交会、后方交会迭代计算，消除RPC参数的相对误差。
[00%]优选的，步骤S2具体包括W下步骤：
[0029] 步骤S201，计算惯性系轨道数据；
[0030] 在惯性参考系的XOZ平面上，取一组卫星测量坐标系原点坐标，各原点坐标与惯性 参考系原点的距离相同，均为R，R等于平均地球半径加上平均轨道高度;假设第一个原点对 应时刻为0，相邻两个质屯、对应时刻相差0.1秒，由R可求得轨道周期T:
[0032] 其中G为万有引力常数，等于6.67*1〇-11米^(千克*秒2)，M为地球质量，等于5.98* 1〇24千克，
[0033] 降交点时刻Tdsd(秒)减去T/4就是第一个质屯、对应的实际时间To,即
[0034] To = Tds<i-T/4
[0035] W下取样时刻依次加0.1秒为相应的实际时间，d0 = 36O/(T*lO)为相邻两个原点 的地屯、角，第一个质屯、在巧由上(Z = R),第i(i = 0，l，一)个取样时刻的卫星原点位置为：
[0036] X = RXcos(O)
[0037] Y = O
[003引 Z = RXsin(目） (1)
[0039] 其中，目= (90-d 目 Xi)/180.0XJi。
[0040] 设轨道倾角为I，将取样坐标绕X轴顺时针旋转1-90度，然后绕巧自逆时针旋转L度； [0041 ] 步骤S202,计算第i个时刻卫星原点(Xi，Yi，Zi)对应的初始大地坐标化i ,Bi，Hi);
[0042] 步骤S203,计算经度L与目标区域中屯、经度Lc的偏差A L
[0043] 步骤S204,卫星原点大地坐标的经度^(i = 0，l，???）加上A L，得到卫星原点的星 下点通过目标区域中屯、点化c，Bc)的卫星原点大地坐标序列化/ i，Bi，出）（（i = 0,1，…），!/ i 二以+A L);
[0044] 步骤S205,将卫星原点大地坐标序列化^i,Bi,出)转换成地屯、直角坐标序列i， Y'i，Z'i)(i = 0，l，...）；
[0045] 步骤S206,返回步骤201，利用公式（I)计算卫星原点^1，￥1，21)时将角度目加上 0.00002度；
[0046] 步骤S207,重复步骤S202至S206,得到新的卫星原点地屯、直角坐标序列(X"i，Y"i， Z"i);
[0047] 步骤S208,计算卫星原点速度序列(Vxi，Vyi，Vzi):
[004引 Vxi = (X"i-X' i)/(T*0.00002/360)
[0049] Vyi= (Y"i-Y'i)/(T*0.0000 2/360)
[0050] Vzi= U"i-Z'i)/(T*0.00002/360)
[(K)引]其中第i组卫星原点位置速度对应时间为:Ti = To+i*0.1;
[0化2] 步骤S209，保存卫星原点轨迹数据(Ti，X' i，Y' i，Z' i，Vxi，Vyi，Vzi)。
[0053] 优选的，步骤S3具体包括W下步骤：
[0054] 步骤301，计算指向误差参数：
[0055] 步骤302，利用功率谱分析法获取姿态参数：
[0056] 步骤S303,将通过步骤S301得到的指向误差参数和通过步骤S302得到的姿态参数 相加得到最终的姿态参数。
[0化7] 其中，步骤S303具体包括：
[0058] 每隔1秒统计生成的姿态角序列f(ti)的均方差0,将姿态角序列的每个姿态角乘 W〇/(3yawi)，然后加上偏航指向角系统量yaw^o,得到偏航姿态角序列：
[0059] yaw(ti) = f (ti)*〇/(3*yawi)+yaw' O
[0060] 将姿态角序列的每个姿态角乘^〇/(391*(：111)，然后加上俯仰指向角系统量 pitch/日，得到俯仰姿态角序列：
[0061 ] pitch(ti) = f (ti)*〇/(3*pitchi)+pitch' O
[0062] 将姿态角序列的每个姿态角乘W日/(3rolh)，然后加上滚动指向角系统量roll%， 得到滚动姿态角序列：
[0063] roll(ti) = f (ti)*〇/(3*rolli)+roll^ 〇〇
[0064] 优选的；
[00化]步骤S401具体包括W下子步骤：
[0066] S4011，将相机光学节点位置结合成像时刻转换到物方坐标系：
[0067] S4012,根据相机安装角计算相机测量坐标系相对卫星本体坐标系的旋转矩阵 Tset;
[006引其中，相机测量坐标系相对卫星本体坐标系的旋转矩阵rset为：
[0069] Mr* M,'
[0070] 其中， -
O
[0074] S4013,计算相机光学节点在卫星本体坐标系的安装向量U在WGS84地球固连坐标 系的向量:
[007引 =
[0076]其中，设相机光学节点在卫星本体坐标系的安装位置为(Xs，Ys，Zs)，则U为卫星本 体坐标系原点到(Xs，Ys，Zs)的向量，r。'为卫星本体坐标系相对t时刻的局部轨道坐标系的旋 转矩阵，r。'为卫星本体坐标系原点与WGS84地球固连坐标系对应旋转矩阵