CN106126839B - 一种三线阵立体测绘卫星成像仿真方法和系统 - Google Patents

一种三线阵立体测绘卫星成像仿真方法和系统 Download PDF

Info

Publication number
CN106126839B
CN106126839B CN201610500333.XA CN201610500333A CN106126839B CN 106126839 B CN106126839 B CN 106126839B CN 201610500333 A CN201610500333 A CN 201610500333A CN 106126839 B CN106126839 B CN 106126839B
Authority
CN
China
Prior art keywords
satellite
origin
image
ccd
sequence
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201610500333.XA
Other languages
English (en)
Other versions
CN106126839A (zh
Inventor
岳庆兴
唐新明
高小明
范大昭
刑帅
耿则勋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
SATELLITE SURVEYING AND MAPPING APPLICATION CENTER NASG
Original Assignee
SATELLITE SURVEYING AND MAPPING APPLICATION CENTER NASG
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by SATELLITE SURVEYING AND MAPPING APPLICATION CENTER NASG filed Critical SATELLITE SURVEYING AND MAPPING APPLICATION CENTER NASG
Priority to CN201610500333.XA priority Critical patent/CN106126839B/zh
Publication of CN106126839A publication Critical patent/CN106126839A/zh
Application granted granted Critical
Publication of CN106126839B publication Critical patent/CN106126839B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/30—Circuit design
    • G06F30/36—Circuit design at the analogue level
    • G06F30/367—Design verification, e.g. using simulation, simulation program with integrated circuit emphasis [SPICE], direct methods or relaxation methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00—Indexing scheme for image data processing or generation, in general
    • G06T2200/08—Indexing scheme for image data processing or generation, in general involving all processing steps from image acquisition to 3D model generation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30168—Image quality inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Remote Sensing (AREA)
  • Computer Graphics (AREA)
  • Microelectronics & Electronic Packaging (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及一种三线阵立体测绘卫星成像仿真方法和系统，其中，所述方法包括：对获取的CCD相机在积分时间段内各个离散时刻的焦面辐照度在时间和级数上积分，得到CCD相机的焦面辐照度；根据通过CCD相机的焦面辐照度转换的入瞳辐亮度计算行影像的影像灰度值；根据通过CCD相机的入瞳辐亮度求取的信噪比计算行影像的噪声影像灰度值；将噪声影像灰度值与影像灰度值相加得到积分时间段对应的行影像的最终影像灰度值；继续获取CCD相机在下一个积分时间段内各个离散时刻的焦面辐照度的步骤，直到计算完所有积分时间段对应的行影像的最终影像灰度值，输出所有行组成的整幅影像。本发明能够提高立体影像仿真的精度。

Description

一种三线阵立体测绘卫星成像仿真方法和系统
技术领域
本发明涉及卫星成像仿真，具体涉及一种三线阵立体测绘卫星成像仿真方法和系统。
背景技术
我国第一颗民用测绘卫星资源三号已于2012年1月成功发射。资源三号卫星集测绘和资源调查功能于一体，用于1:5万立体测图及更大比例尺基础地理信息产品的生产和更新，开展国土资源调查与监测。该卫星工程的建设，将大力增强我国独立获取地理空间信息的能力，提升我国测绘服务保障水平，提高国土资源调查与监测的数据保障能力，加快空间数据基础设施建设，推动地理信息产业发展。除了搭载一台多光谱相机外，资源三号的主要载荷为由正视、前视、后视三台相机构成的三线阵立体测绘相机，其中正视相机分辨率约为2.1米，前视、后视相机分辨率约为3.6米。三线阵立体测绘相机可以在同一轨道获取高质量的三视立体影像，进而可以通过立体匹配获取高质量的数字表面模型(Digital SurfaceModel，简称DSM)。除了资源三号卫星以外，日本的ALOS卫星也是典型的三线阵立体测绘卫星，三个相机的分辨率均约为2.5米。法国的Pleiades卫星可以在同一轨道通过单相机沿轨道方向摆动实现三视立体成像，三个视角的影像分辨率均约为0.5米。
卫星工程属于国家重大工程，每颗高分辨率光学卫星的研制费用都很高。用户需要有一套可靠的技术方法对卫星获取的影像质量进行尽量准确的评估，即需要对卫星在轨拍摄影像的过程进行精确的仿真，用仿真影像对卫星的成像质量进行评估，甚至反馈到卫星的设计中，优化后续星的某些设计指标以期能够获取更高质量的影像。同时，可以用仿真影像对应用系统进行测试，验证应用系统关键技术效果，保障整个卫星目标的顺利实现，保障国家重大投资的成功率。
在可见、近红外光学成像卫星的成像仿真技术领域，国内外学者提出了物理仿真、半物理仿真和计算机数值仿真三种技术途径，其中计算机数值仿真的工程成本最低。采用各子系统的实验室标定数据结合设计数据建立各子系统的数学模型，也就是用数学公式描述信号传输和转换的物理过程，简称子系统建模。在子系统建模的基础上实现在轨动态成像过程的全链路的仿真流程，能比较全面地模拟卫星的设计技术指标对卫星在轨成像质量的影响，准确地模拟出延迟积分成像高分辨率光学卫星在轨成像过程中几何和辐射质量的下降与影像总体几何、辐射精度的对应关系，建立起卫星设计参数与在轨成像质量之间的直接关联。
真实的卫星相机成像是对真实的地面物理场景进行辐射信息采集和转化的过程，进行计算机数值仿真同样需要一个虚拟的数字物理场景，也就是一个数值化的地面模型，这个地面模型需要包括地面每一点的几何信息和辐射信息。换句话说，地面模型由几何模型和辐射模型构成，几何模型提供每一地面点的几何信息，辐射模型提供地面每一点的辐射信息。
以DEM和DOM为输入数据的卫星成像仿真比较容易实现，但这种输入几何模型中仅包含地形信息，而房屋、树木等地物的几何信息都被过滤掉了。DOM虽然可以采用较高分辨率的数据从而包含精细的纹理信息，但这些纹理信息与几何信息在地物丰富的区域往往不能统一。无法满足三线阵立体测绘卫星的主要评估目标之一DSM精度的预评估。
发明内容
本发明所要解决的技术问题是提供一种三线阵立体测绘卫星成像仿真方法和系统，能够提高立体影像仿真的精度。
本发明解决上述技术问题的技术方案如下：一种三线阵立体测绘卫星成像仿真方法，包括：
S1，获取CCD相机在积分时间段内各个离散时刻的焦面辐照度；
S2，对各个离散时刻的焦面辐照度在积分时间段上积分，和在级数上积分，得到CCD相机的焦面辐照度；
S3，将CCD相机的焦面辐照度转换为CCD相机的入瞳辐亮度；
S4，根据CCD相机的入瞳辐亮度计算所述积分时间段对应的行影像的影像灰度值；
S5，根据CCD相机的入瞳辐亮度求取信噪比，并根据所述信噪比计算所述积分时间段对应的行影像的噪声影像灰度值；
S6，将所述噪声影像灰度值与所述影像灰度值相加得到所述积分时间段对应的行影像的最终影像灰度值；
S7，返回S1继续执行，计算下一积分时间段对应的行影像的最终影像灰度值，直到计算完所有积分时间段对应的行影像的最终影像灰度值，根据所有行影像的最终影像灰度值输出所有行组成的整幅影像。
本发明的有益效果是：通过对获取CCD相机在积分时间段内各个离散时刻的焦面辐照度在积分时间段上积分，和在级数上积分，得到CCD相机的焦面辐照度；将CCD相机的焦面辐照度转换为CCD相机的入瞳辐亮度；根据CCD相机的入瞳辐亮度计算所述积分时间段对应的行影像的影像灰度值；根据CCD相机的入瞳辐亮度求取信噪比，并根据所述信噪比计算所述积分时间段对应的行影像的噪声影像灰度值；将所述噪声影像灰度值与所述影像灰度值相加得到所述积分时间段对应的行影像的最终影像灰度值；继续获取CCD相机在下一个积分时间段内各个离散时刻的焦面辐照度的步骤，直到计算完所有积分时间段对应的行影像的最终影像灰度值，根据所有行影像的最终影像灰度值输出所有行组成的整幅影像，从而提高立体影像仿真的精度。
进一步，所述获取CCD相机在积分时间段内各个离散时刻的焦面辐照度，具体为：
将CCD相机中每个CCD划分为多个子CCD，得到若干个子CCD；
求取每个子CCD在离散时刻t的入瞳辐亮度；其中，离散时刻t为CCD相机在一个积分时间段中的某一时刻；
根据子CCD在t时刻的入瞳辐亮度求取子CCD在t时刻的焦面辐照度并滤波；
将每个CCD内部的所有子CCD的滤波后的焦面辐照度取平均值，得到每个CCD相机在t时刻的焦面辐照度。
进一步，所述求取每个子CCD在离散时刻t的入瞳辐亮度，具体为：
将输入的几何模型剖分得到多个物方三角形，并对所述多个物方三角形编号；
将每个所述物方三角形的三个顶点坐标分别投射到CCD相机上，得到多个像方三角形；
统计落在每个像方三角形内的像点坐标，并记录所述像方三角形内每个像点对应的多个物方三角形的编号；
根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与所述物方三角形的交点；
将距离相机光学节点最近的交点作为实际交点；
获取预先输入的三视Pleiades立体卫星影像中与CCD相机视角对应的Pleiades立体卫星影像上与实际交点对应位置处的影像灰度值，并反算为实际交点入瞳辐亮度，得到子CCD的入瞳辐亮度。
采用上述进一步方案的有益效果是：通过以三视Pleiades立体影像为纹理信息数据源，以高精度Lidar数据为几何信息数据源，将Lidar数据剖分为三角面片并向模拟相机投射建立索引，并在获取实际交点后，从预先输入的三视Pleiades立体卫星影像中与CCD相机视角对应的Pleiades立体卫星影像上获取与实际交点对应位置处的影像灰度值，并反算为实际交点入瞳辐亮度，得到子CCD的入瞳辐亮度，进而达到快速计算模拟相机光线对应的地面目标的效果，同时进一步提高立体影像仿真的精度。
进一步，所述根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与所述物方三角形的交点的步骤中，具体为根据卫星轨道运行参数、卫星姿态参数和相机安装模型计算每个像点对应的主光线向量，并转换到物方坐标系。
进一步，所述卫星轨道运行参数按照以下方法获取：
计算卫星轨道周期T；具体为：
式中，G为万有引力常数；M为地球质量；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
计算第1个离散时刻卫星原点对应的实际时间T0；具体为：
T0＝Tdsd-T/4；
式中，Tdsd为降交点时刻；
计算第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)；具体为：
θ＝(90-dθ×i)/180×π；
X＝R×cos(θ)；
Y＝0；
Z＝R×sin(θ)；
式中，θ为地标原点的地心角；(X,Y,Z)为第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)中任一组数据；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
将第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)绕X顺时针旋转I-90度后绕Z轴逆时针旋转L度，得到第一次旋转后的第i个离散时刻的卫星原点位置序列(X′i,Y′i,Z′i)；I为卫星轨道倾角，L为待拍摄目标区域中心的经度；
对第一次旋转后的第i个离散时刻的卫星原点位置序列(X′i,Y′i,Z′i)绕Z轴旋转i×dL度，得到第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi″,Yi″,Zi″)；dL为0.1秒内地球的自转角度；dL＝(1/2400)°；
根据第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi″,Yi″,Zi″)中任一个卫星原点坐标(Xi,Yi,Zi)对应的初始大地坐标序列中(Li,Bi,Hi)的初始大地坐标(L,B,H)计算第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)；具体为：
L＝atan(Y/X)
H＝P/cos(B)-N；
式中，(L,B,H)为(Li,Bi,Hi)序列的任一组数据，(X,Y,Z)为对应(L,B,H)的(Xi,Yi,Zi)序列中的数据；a＝6378137；
将第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)的经度Li加上ΔL，得到卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L′i,Bi,Hi)；其中，Li为第i个离散时刻的卫星原点的经度，Bi为第i个离散时刻的卫星原点的纬度，Hi为第i个离散时刻的卫星原点的高度；i为离散时刻的编号，i＝0,1,L N；Lc为卫星原点的星下点通过目标区域中心点的经度，Bc为卫星原点的星下点通过目标区域中心点的纬度；ΔL＝Lc-L+RandL；
RandL＝1.0*rand()/RAND_MAX*(SW-GW)/2/6371000/(2π)*360；
式中，GW为输入的几何模型宽度；SW为Pleiades立体卫星影像的地面幅宽；rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值；
将卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L′i,Bi,Hi)转换成卫星原点地心直角坐标序列(X″′i,Y″′i,Z″′i)，i＝0,1,L N，具体为：
式中，(L′,B,H)为(L′i,Bi,Hi)序列的任一组数据，(X″′,Y″′,Z″′)为对应(X″′i,Y″′i,Z″′i)序列的任一组数据；
将地标原点的地心角θ加上阈值，求取新的卫星原点地心直角坐标序列(X″″i,Y″″i,Z″″i)；
根据卫星原点地心直角坐标序列(X″′i,Y″′i,Z″′i)和新的卫星原点地心直角坐标序列(X″″i,Y″″i,Z″″i)计算卫星原点速度序列(Vxi,Vyi,Vzi)：
Vxi＝(X″″i-X″′i)/(T*0.00002/360)
Vyi＝(Y″″i-Y″′i)/(T*0.00002/360)
Vzi＝(Z″″i-Z″′i)/(T*0.00002/360)
式中，T为轨道周期；
根据第i个离散时刻卫星原点位置、速度和两者对应的实际时间Ti，得到卫星轨道运行参数(Ti,X′i,Y′i,Z′i,Vxi,Vyi,Vzi)；式中，Ti＝T0+i*0.1，T0为第1个离散时刻对应的实际时间。
采用上述进一步方案的有益效果是：通过轨道的设置使小的目标区域随机成像在整景大区域影像的任意位置。通过多个小区域的成像仿真和精度评估来解决大区域影像测绘精度评估的问题。
进一步，卫星姿态参数的获取包括以下步骤：
将t+Randt时刻的卫星姿态参数作为t时刻的卫星姿态参数，其中t表示第i个离散时刻对应的时间；式中，Randt表示时间随机变量，Randt＝1.0*rand()/RAND_MAX*(SW-GH)/2/6371000/(2π)*T；
式中，T为轨道周期，GH为高度，SW为仿真影像地面幅宽，rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值。
采用上述进一步方案的有益效果是：通过成像时间的设置使小的目标区域随机成像在整景大区域影像的任意位置。通过多个小区域的成像仿真和精度评估来解决大区域影像测绘精度评估的问题。
本发明解决上述技术问题的另一种技术方案是：一种三线阵立体测绘卫星成像仿真系统，包括：
焦面辐照度获取模块，用于获取CCD相机在积分时间段内各个离散时刻的焦面辐照度；
积分模块，用于对各个离散时刻的焦面辐照度在积分时间段上积分，和在级数上积分，得到CCD相机的焦面辐照度；
转换模块，用于将CCD相机的焦面辐照度转换为CCD相机的入瞳辐亮度；
影像灰度值计算模块，用于根据CCD相机的入瞳辐亮度计算所述积分时间段对应的行影像的影像灰度值；
噪声影像灰度值获取模块，用于根据CCD相机的入瞳辐亮度求取信噪比，并根据所述信噪比计算所述积分时间段对应的行影像的噪声影像灰度值；
最终影像灰度值获取模块，用于将所述噪声影像灰度值与所述影像灰度值相加得到所述积分时间段对应的行影像的最终影像灰度值；
循环及输出模块，循环调用焦面辐照度获取模块、积分模块、转换模块、影像灰度值计算模块、噪声影像灰度值获取模块和最终影像灰度值获取模块，计算下一积分时间段对应的行影像的最终影像灰度值，直到计算完所有积分时间段对应的行影像的最终影像灰度值，根据所有行影像的最终影像灰度值输出所有行组成的整幅影像。
本发明的有益效果是：通过积分模块对获取CCD相机在积分时间段内各个离散时刻的焦面辐照度在积分时间段上积分，和在级数上积分，得到CCD相机的焦面辐照度；转换模块将CCD相机的焦面辐照度转换为CCD相机的入瞳辐亮度；影像灰度值计算模块根据CCD相机的入瞳辐亮度计算所述积分时间段对应的行影像的影像灰度值；噪声影像灰度值获取模块根据CCD相机的入瞳辐亮度求取信噪比，并根据所述信噪比计算所述积分时间段对应的行影像的噪声影像灰度值；最终影像灰度值获取模块将所述噪声影像灰度值与所述影像灰度值相加得到所述积分时间段对应的行影像的最终影像灰度值；循环及输出模块继续获取CCD相机在下一个积分时间段内各个离散时刻的焦面辐照度的步骤，直到计算完所有积分时间段对应的行影像的最终影像灰度值，根据所有行影像的最终影像灰度值输出所有行组成的整幅影像，从而提高立体影像仿真的精度。
本发明在解决上述技术问题的基础上，还能做出如下改进：
进一步，所述焦面辐照度获取模块包括：
划分子模块，用于将CCD相机中每个CCD划分为多个子CCD，得到若干个子CCD；
入瞳辐亮度获取子模块，用于求取每个子CCD在离散时刻t的入瞳辐亮度；其中，离散时刻t为CCD相机在一个积分时间段中的某一时刻；
子CCD焦面辐照度获取子模块，用于根据子CCD在t时刻的入瞳辐亮度求取子CCD在t时刻的焦面辐照度并滤波；
焦面辐照度获取子模块，将每个CCD内部的所有子CCD的滤波后的焦面辐照度取平均值，得到每个CCD相机在t时刻的焦面辐照度。
进一步，所述入瞳辐亮度获取子模块包括：
剖分单元，用于将输入的几何模型剖分得到多个物方三角形，并对所述多个物方三角形编号；
投射单元，用于将每个所述物方三角形的三个顶点坐标分别投射到CCD相机上，得到多个像方三角形；
统计及记录单元，用于统计落在每个像方三角形内的像点坐标，并记录所述像方三角形内每个像点对应的多个物方三角形的编号；
转换单元，用于根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与所述物方三角形的交点；
实际交点获取单元，用于将距离相机光学节点最近的交点作为实际交点；
子CCD入瞳辐亮度获取单元，用于获取预先输入的三视Pleiades立体卫星影像中与CCD相机视角对应的影像上与实际交点对应位置处的影像灰度值，并反算为实际交点入瞳辐亮度，得到子CCD的入瞳辐亮度。
进一步，所述转换单元在根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与所述物方三角形的交点的步骤中，具体为根据卫星轨道运行参数、卫星姿态参数和相机安装模型计算每个像点对应的主光线向量，并转换到物方坐标系。
进一步，所述卫星轨道运行参数按照以下方法获取：
计算卫星轨道周期T；具体为：
式中，G为万有引力常数；M为地球质量；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
计算第1个离散时刻卫星原点对应的实际时间T0；具体为：
T0＝Tdsd-T/4；
式中，Tdsd为降交点时刻；
计算第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)；具体为：
θ＝(90-dθ×i)/180×π；
X＝R×cos(θ)；
Y＝0；
Z＝R×sin(θ)；
式中，θ为地标原点的地心角；(X,Y,Z)为第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)中任一组数据；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
将第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)绕X顺时针旋转I-90度后绕Z轴逆时针旋转L度，得到第一次旋转后的第i个离散时刻的卫星原点位置序列(X′i,Y′i,Z′i)；I为卫星轨道倾角，L为待拍摄目标区域中心的经度；
对第一次旋转后的第i个离散时刻的卫星原点位置序列(Xi′,Yi′,Zi′)绕Z轴旋转i×dL度，得到第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi″,Yi″,Zi″)；dL为0.1秒内地球的自转角度；dL＝(1/2400)°；
根据第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi″,Yi″,Zi″)中任一个卫星原点坐标(Xi,Yi,Zi)对应的初始大地坐标序列中(Li,Bi,Hi)的初始大地坐标(L,B,H)计算第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)；具体为：
L＝atan(Y/X)
H＝P/cos(B)-N；
式中，(L,B,H)为(Li,Bi,Hi)序列的任一组数据，(X,Y,Z)为对应(L,B,H)的(Xi,Yi,Zi)序列中的数据；a＝6378137；
将第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)的经度Li加上ΔL，得到卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L′i,Bi,Hi)；其中，Li为第i个离散时刻的卫星原点的经度，Bi为第i个离散时刻的卫星原点的纬度，Hi为第i个离散时刻的卫星原点的高度；i为离散时刻的编号，i＝0,1,L N；Lc为卫星原点的星下点通过目标区域中心点的经度，Bc为卫星原点的星下点通过目标区域中心点的纬度；ΔL＝Lc-L+RandL；
RandL＝1.0*rand()/RAND_MAX*(SW-GW)/2/6371000/(2π)*360；
式中，GW为输入的几何模型宽度；SW为Pleiades立体卫星影像的地面幅宽；rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值；
将卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L′i,Bi,Hi)转换成卫星原点地心直角坐标序列(X″′i,Y″′i,Z″′i)，i＝0,1,L N，具体为：
式中，(L′,B,H)为(L′i,Bi,Hi)序列的任一组数据，(X″′,Y″′,Z″′)为对应(X″′i,Y″′i,Z″′i)序列的任一组数据；
将地标原点的地心角θ加上阈值，求取新的卫星原点地心直角坐标序列(X″″i,Y″″i,Z″″i)；
根据卫星原点地心直角坐标序列(X″′i,Y″′i,Z″′i)和新的卫星原点地心直角坐标序列(X″″i,Y″″i,Z″″i)计算卫星原点速度序列(Vxi,Vyi,Vzi)：
Vxi＝(X″″i-X″′i)/(T*0.00002/360)
Vyi＝(Y″″i-Y″′i)/(T*0.00002/360)
Vzi＝(Z″″i-Z″′i)/(T*0.00002/360)
式中，T为轨道周期；
根据第i个离散时刻卫星原点位置、速度和两者对应的实际时间Ti，得到卫星轨道运行参数(Ti,X′i,Y′i,Z′i,Vxi,Vyi,Vzi)；式中，Ti＝T0+i*0.1，T0为第1个离散时刻对应的实际时间。
采用上述进一步方案的有益效果是：通过轨道的设置使小的目标区域随机成像在整景大区域影像的任意位置。通过多个小区域的成像仿真和精度评估来解决大区域影像测绘精度评估的问题。
进一步，卫星姿态参数的获取包括以下步骤：
将t+Randt时刻的卫星姿态参数作为t时刻的卫星姿态参数，其中t表示第i个离散时刻对应的时间；式中，Randt表示时间随机变量，Randt＝1.0*rand()/RAND_MAX*(SW-GH)/2/6371000/(2π)*T；
式中，T为轨道周期，GH为高度，SW为仿真影像地面幅宽，rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值。
采用上述进一步方案的有益效果是：通过成像时间的设置使小的目标区域随机成像在整景大区域影像的任意位置。通过多个小区域的成像仿真和精度评估来解决大区域影像测绘精度评估的问题。
附图说明
图1为本发明一种三线阵立体测绘卫星成像仿真方法的流程示意图；
图2为仿真影像像点与地面三角面片对应关系示意图；
图3为本发明一种三线阵立体测绘卫星成像仿真系统的结构示意图。
具体实施方式
以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
图1为本发明一种三线阵立体测绘卫星成像仿真方法的流程示意图。
如图1所示，一种三线阵立体测绘卫星成像仿真方法，包括：
S1，获取CCD相机在积分时间段内各个离散时刻的焦面辐照度；具体为：
S11，将CCD相机中每个CCD划分为多个子CCD，得到若干个子CCD；
S12，求取每个子CCD在离散时刻t的入瞳辐亮度；其中，离散时刻t为CCD相机在一个积分时间段中第i个离散时刻对应的时刻；
求取每个子CCD的入瞳辐亮度；具体如下：
S121，将输入的几何模型剖分得到多个物方三角形，并对多个物方三角形编号；其中，几何模型为Lidar数据；
S122，将每个物方三角形的三个顶点坐标分别投射到CCD相机上，得到多个像方三角形；
S123，统计落在每个像方三角形内的像点坐标，并记录像方三角形内每个像点对应的多个物方三角形的编号；
如图2所示，例如：Lidar数据的平面坐标构成的矩形网格被对角线剖分为两个物方三角形，两个物方三角形根据上下顺序用0和1区分，总体位置用左上角矩形格网点行列坐标(R_1d，C_1d)表示，投射到CCD相机得到两个像方三角形，像点P将位于两个像方三角形内部，分别用P1和P2表示，则该像点P对应的物方三角形的编号为(R_1d，C_1d，0)和(R_1d，C_1d，1)。
S124，根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与物方三角形的交点；
上述步骤S124根据卫星轨道运行参数、卫星姿态参数和CCD相机安装模型将每个像点对应的主光线向量转换到物方坐标系；具体包括如下步骤：
S1241，设CCD相机的TDI CCD焦平面上任意一个点p的像平面坐标为(x,y)，x＝(c-c0)*d；y＝r*d，单位为米；将p点的像平面坐标加入畸变量得到加入畸变量后的像平面坐标p′(x′,y′)，如下：
式中，c为p点在对应视角CCD相机的TDI CCD焦平面的列坐标，单位为(像素，浮点数)，c＝i-0.5+m/k，i为CCD的级数；r为p点在对应视角CCD相机的TDI CCD焦平面的行坐标，单位为(像素，浮点数)，r＝j-0.5+n/k，j为CCD的列数；σAS和σCS分别为p点沿TDI CCD焦平面积分方向和垂直TDI CCD焦平面积分方向的畸变量；c0为像方主点在p点所在CCD相机中的列坐标，在每个CCD相机中，单个CCD长、宽均为d，对于正视CCD相机，d＝0.000007米；对于前视CCD相机和后视CCD相机，d＝0.00001米；
其中，对于p点来说，畸变量σAS和σCS分别用σASj和σCSj来表示；j为CCD序号；
则沿垂直TDI CCD焦平面积分方向的畸变量σCSj按照以下方法获取：
垂直TDI CCD焦平面积分方向的畸变量格式为“CCD列坐标－畸变量”数据序列，
两个CCD之间的CCD畸变量通过线性插值的方法获得：
通过插值可以得到每个CCD在垂直TDI CCD焦平面积分方向的畸变量σCSj。
沿TDI CCD焦平面积分方向的CCD畸变量σASj按照以下方法获取：
沿TDI CCD焦平面积分方向的畸变用CCD阵列拼接误差、CCD非直线性等指标来衡量。在三线阵相机TDI CCD焦平面积分方向的CCD畸变模型建模中主要考虑CCD阵列拼接误差σc和CCD非直线性σr两个指标，并分两步实现：
首先，以第一片CCD(即第一个CCD阵列)为基准，将第二片CCD中每个CCD相对理想安装位置均偏移σc1(σc1为绝对值小于σc的随机数)，再以第二片CCD为基准，将第三片CCD中每个CCD相对理想安装位置均偏移σc2(σc2为绝对值小于σc的随机数)；然后，将每片CCD以各自的片CCD的中心CCD为中心，旋转σrj(σrj为绝对值小于σr的随机数，i＝1,2,3，为CCD阵列号)。其他片CCD以此类推，最终可以得到每个CCD在沿TDI积分方向的畸变量σASj＝σcj+σrj。
在模拟算法中认为该畸变是光学畸变和CCD畸变的综合值；
对于两个CCD列之间的亚像元点，畸变量σAS和σCS分别用σAS(c)和σCS(c)表示，具体如下：
σAS(c)＝σASj*(1-temp)+σAS(j+1)*temp
σCS(c)＝σCSj*(1-temp)+σCS(j+1)*temp
temp＝c-int(c)
式中，c为亚像元点列坐标，j≤c≤j+1；int(c)表示取c的整数部分；
S1242，计算p′在对应视角CCD相机的像方坐标系的主光线向量u：
式中，-f是相机的焦距；x′是加入了畸变量的点p的行坐标；y′是加入了畸变量的点p的列坐标；
S1243，将主光线向量u转换到卫星本体坐标系得向量u1：具体按照以下方法转换：
u1＝rset*u
其中，rset为CCD相机的安装矩阵，具体按照以下方法获取：
CCD相机的安装矩阵由CCD相机安装模型获取，CCD相机安装模型包括CCD相机光学节点在卫星本体坐标系的安装位置(Xs,Ys,Zs)和CCD相机测量坐标系相对卫星本体坐标系的安装角(yaws,pitchs,rolls)，yaws为CCD相机测量坐标系相对卫星本体坐标系的安装角中的俯仰角；pitchs为CCD相机测量坐标系相对卫星本体坐标系的安装角中的滚动角；rolls为CCD相机测量坐标系相对卫星本体坐标系的安装角中的偏航角；理想情况下，正视相机的安装角为(0,0,0)，前视相机的安装角为(0,pitchF,0)，后视相机的安装角为(0,pitchB,0)。由于相机安装的过程中会存在一定的偏差，实际的安装角并不等于理想值，因此，相机安装角可以通过实验室标定方法获得，本发明使用的相机安装角为实验室标定结果。相机测量坐标系相对卫星本体坐标系的旋转矩阵rset为：
其中，
式中，表示CCD相机测量坐标系相对卫星本体坐标系的安装角中pitchs角的余弦值；表示CCD相机测量坐标系相对卫星本体坐标系的安装角中pitchs角的正弦值；为的相反数；表示CCD相机测量坐标系相对卫星本体坐标系的安装角中rolls角的余弦值；表示CCD相机测量坐标系相对卫星本体坐标系的安装角中rolls角的正弦值；为的相反数；表示CCD相机测量坐标系相对卫星本体坐标系的安装角中yaws角的余弦值；表示CCD相机测量坐标系相对卫星本体坐标系的安装角中yaws角的正弦值；为的相反数；
S1244，将光线向量u1转换到局部轨道坐标系得到光线向量u2：
u2＝rplm2obt*u1
其中，rplm2obt为卫星本体坐标系到局部轨道坐标系的转换矩阵，计算方法为：
rplm2obt＝rpitch*rroll*ryaw
其中，pitch、roll、yaw为卫星本体坐标系相对局部轨道坐标系的三个姿态角，yaw为偏航角，pitch为俯仰角，roll为滚动角；三个姿态角pitch、roll、yaw的获取具体为：
首先用功率谱分析法获取符合卫星设计姿态稳定度的姿态曲线，然后生成符合卫星设计姿态指向精度的系统指向误差，将指向误差参数和卫星姿态参数相加得到最终的卫星姿态参数。
其中，使用功率谱分析法获取姿态曲线的方法为：
步骤S401，根据卫星主要姿态设计指标中的三轴姿态指向精度(yaw0,pitch0,roll0)生成三轴指向角系统量：
yaw0′＝1.0·rand()/RAND_MAX·yaw0
pitch0′＝1.0·rand()/RAND_MAX·pitch0
roll0′＝1.0·rand()/RAND_MAX·roll0
式中，rand()为C++的均匀分布随机整数生成器，RAND_MAX为最大随机整数；
步骤S402，根据模拟设置的多个频段v、频段宽度w及该频段的幅值A(wi′,vi′,Ai′)(i′表示第i′组模拟频段参数)生成滤波器：
式中，n为设置的频段个数，Q(v)为频率v对应的滤波系数，是一个归一化比例系数。
步骤S403，生成随机振动功率谱；
式中，F(v)为频率v对应的功率值；
步骤S404，频率域滤波得到姿态振动功率谱：
F(v)′＝Q(v)*F(v)
式中，F(v)′为频率v对应的滤波后的功率值；
步骤S405，反傅立叶变换变回时间域得到一个姿态角的姿态数据序列：
f(ti″)＝IFFT(F(v)′)
其中，i″表示第i″个姿态角，ti″表示第i″个姿态角的时间，设频率域的最高频率为vmax，则ti″＝1.0/vmax*i″。
步骤S406，每隔1秒统计生成的姿态角序列f(ti″)的均方差σ，将姿态角序列的每个姿态角乘以σ/(3yaw1)，然后加上偏航指向角系统量yaw′0，得到偏航姿态角序列：
yaw(ti″)＝f(ti″)*σ/(3*yaw1)+yaw′0
将姿态角序列的每个姿态角乘以σ/(3pitch1)，然后加上偏航指向角系统量pitch′0，得到偏航姿态角序列：
pitch(ti″)＝f(ti″)*σ/(3*pitch1)+pitch′0
将姿态角序列的每个姿态角乘以σ/(3roll1)，然后加上偏航指向角系统量roll′0，得到偏航姿态角序列roll(ti″)。
根据离散时刻对应的时间t和第i″个姿态角的时间ti″就可以从获取的姿态角序列{yaw(ti″),roll(ti″),pitch(ti″)}得到离散时刻对应的时间的三个姿态角pitch、roll、yaw。
由于输入的Pleiades立体卫星影像的数据范围远小于仿真三线阵影像的数据范围，最终的仿真影像只会覆盖整景仿真影像的一部分区域和高度。为了更好的评估整景三线阵影像的几何精度，需要进行多次实验，使有效仿真影像在任意姿态区间的立体几何精度。要实现这个目标需要每次仿真时将t+Randt时刻的卫星姿态参数作为t时刻的卫星姿态参数，其中t表示第i个离散时刻对应的时间；式中，Randt表示时间随机变量，Randt＝1.0*rand()/RAND_MAX*(SW-GH)/2/6371000/(2π)*T；其中，卫星姿态参数指卫星本体坐标系相对局部轨道坐标系的三个姿态角；
式中，T为轨道周期，GH为高度，SW为仿真影像地面幅宽，rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值。
S1245，将光线向量u2转换到WGS84地心直角坐标系得到像点p的主光线向量u3：
u3＝ropt2wgs84*u2
其中，
为t时刻CCD相机光学节点位置的归一化向量；
为t时刻CCD相机光学节点速度的归一化向量；
[Cx,Cy,Cz]T为和的叉积。
t时刻CCD相机光学节点位置的归一化向量根据以下方法获取：
式中，是CCD相机光学节点在卫星本体坐标系的安装位置向量在地球固连坐标系的形式；为相机光学节点在卫星本体坐标系的安装位置向量在地球固连坐标系的形式；
式中，卫星本体坐标系原点在WGS84地球固连坐标系的位置和速度矢量为对应的旋转矩阵：
为卫星原点位置归一化向量；
为卫星原点速度归一化向量；
[Cx,Cy,Cz]T为和的叉积；
U为相机光学节点在卫星本体坐标系的安装向量，由已知的卫星本体坐标系原点到相机光学节点在卫星本体坐标系的安装位置[Xs,Ys,Zs]的向量获取；
t时刻CCD相机光学节点速度根据以下公式获取：
式中，Δt为CCD相机成像的积分时间；为t+Δt时刻的相机光学节点的运行位置；
t时刻位置速度向量由t时刻前后各四组卫星轨迹运行数据经拉格朗日插值获得；其中，t为第i个离散时刻对应的时间；
t时刻卫星轨迹运行数据(Ti,X′i,Y′i,Z′i,Vxi,Vyi,Vzi)按照以下方法求取：
设卫星轨迹运行数据用CCD相机成像时间段内某一时刻卫星本体坐标系原点的位置及速度矢量(Xo t,Yo t,Zo t,VXo t,VYo t,VZo t)来表示，具体按照以下方法获取：
若仿真卫星轨迹为太阳同步可重访圆轨迹，轨迹平均高度为H，轨迹倾角为I，降交点时间为Tdsd。本发明用规则圆轨迹来描述卫星轨迹，卫星轨迹数据的获取具体包括以下步骤：
卫星轨迹运行参数按照以下方法获取：
在惯性参考系的XOZ平面上，取一组卫星测量坐标系原点坐标，各原点坐标与惯性参考系原点的距离均相等，为R，R等于平均地球半径加上平均轨迹高度，计算卫星轨迹周期T；具体为：
式中，G为万有引力常数；M为地球质量；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
计算第1个离散时刻卫星原点对应的实际时间T0；具体为：
T0＝Tdsd-T/4；
式中，Tdsd为降交点时刻；
计算第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)；具体为：
θ＝(90-dθ×i)/180×π；
X＝R×cos(θ)；
Y＝0；
Z＝R×sin(θ)；
式中，θ为地标原点的地心角；(X,Y,Z)为第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)中任一组数据；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
将第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)绕X顺时针旋转I-90度后绕Z轴逆时针旋转L度，得到第一次旋转后的第i个离散时刻的卫星原点位置序列(X′i,Y′i,Z′i)；I为卫星轨迹倾角，L为待拍摄目标区域中心的经度；
对第一次旋转后的第i个离散时刻的卫星原点位置序列(Xi′,Yi′,Zi′)绕Z轴旋转i×dL度，得到第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi″,Yi″,Zi″)；dL为0.1秒内地球的自转角度；dL＝(1/2400)°；
根据第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi″,Yi″,Zi″)中任一个卫星原点坐标(Xi,Yi,Zi)对应的初始大地坐标序列中(Li,Bi,Hi)的初始大地坐标(L,B,H)计算第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)；具体为：
L＝atan(Y/X)
H＝P/cos(B)-N；
式中，(L,B,H)为(Li,Bi,Hi)序列的任一组数据，(X,Y,Z)为对应(L,B,H)的(Xi,Yi,Zi)序列中的数据；a＝6378137；
将第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)的经度Li加上ΔL，得到卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L′i,Bi,Hi)；其中，Li为第i个离散时刻的卫星原点的经度，Bi为第i个离散时刻的卫星原点的纬度，Hi为第i个离散时刻的卫星原点的高度；i为离散时刻的编号，i＝0,1,L N；Lc为卫星原点的星下点通过目标区域中心点的经度，Bc为卫星原点的星下点通过目标区域中心点的纬度；ΔL＝Lc-L+RandL；
RandL＝1.0*rand()/RAND_MAX*(SW-GW)/2/6371000/(2π)*360；
式中，GW为输入的几何模型宽度；SW为Pleiades立体卫星影像的地面幅宽；rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值；
将卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L′i,Bi,Hi)转换成卫星原点地心直角坐标序列(X″′i,Y″′i,Z″′i)，i＝0,1,L N，具体为：
式中，(L′,B,H)为(L′i,Bi,Hi)序列的任一组数据，(X″′,Y″′,Z″′)为对应(X″′i,Y″′i,Z″′i)序列的任一组数据；
将地标原点的地心角θ加上阈值，求取新的卫星原点地心直角坐标序列(X″″i,Y″″i,Z″″i)；
根据卫星原点地心直角坐标序列(X″′i,Y″′i,Z″′i)和新的卫星原点地心直角坐标序列(X″″i,Y″″i,Z″″i)计算卫星原点速度序列(Vxi,Vyi,Vzi)：
Vxi＝(X″″i-X″′i)/(T*0.00002/360)
Vyi＝(Y″″i-Y″′i)/(T*0.00002/360)
Vzi＝(Z″″i-Z″′i)/(T*0.00002/360)
式中，T为轨迹周期；
根据第i个离散时刻卫星原点位置、速度和两者对应的实际时间Ti，得到卫星轨迹运行参数(Ti,X′i,Y′i,Z′i,Vxi,Vyi,Vzi)；式中，Ti＝T0+i*0.1，T0为第1个离散时刻对应的实际时间。
S125，将距离相机光学节点最近的交点作为实际交点；具体为：
计算S1214中得到的主光线向量u3与输入的几何模型的交点，选取距离相机光学节点最近的交点作为实际交点；具体如下：
设距离时刻t最近的整数行影像对应时刻为t0，t0对应整数行为r′(r′＝r或r′＝r+1)，单行CCD积分时间为dt，TDI CCD行坐标r对应的整数行为rI，TDI CCD列坐标c对应整数列为cI，则像点对应的物方三角形对应的整数影像行列(R,C)为：
R＝rI+INT(r-rI+(t-t0)/dt+0.5)
C＝cI+INT(c-cI+0.5)
式中，INT表示取整。
在以(R,C)为中心的子CCD窗口内计算u3与该窗口内所有三角面片索引对应物方三角形的交点序列，统计离t时刻投影中心最近的几何模型上的点，作为地面点，并将地面点转换为大地坐标，取最小的大地坐标，作为实际交点的坐标。
S126，获取预先输入的三视Pleiades立体卫星影像中与CCD相机视角对应的影像上与实际交点对应位置处的影像灰度值，并反算为实际交点入瞳辐亮度，得到子CCD的入瞳辐亮度。三视Pleiades立体卫星影像与输入的几何模型同区域，包含三个视角的影像信息，是分别与三个CCD相机视角对应的三个视角的Pleiades立体卫星影像，Pleiades立体卫星影像为三视Pleiades立体卫星影像中的一幅，与其中一个CCD相机的视角相同，具体为：
根据实际交点和对应视角的Pleiades立体卫星影像的RPC参数计算得到该实际交点在Pleiades立体卫星影像上影像像点坐标，通过双线性内插得到该影像像点的灰度值，以及根据Pleiades立体卫星影像的辐射校正参数将影像灰度值反算为入瞳辐亮度，得到子CCD的入瞳辐亮度；RPC参数为几何模型上的物点与对应像点的几何转换参数。
S13，根据子CCD在t时刻的入瞳辐亮度求取子CCD在t时刻的焦面辐照度并滤波，得到子CCD在t时刻的实际焦面辐照度；具体按照以下方法求取：
其中，根据子CCD在t时刻的入瞳辐亮度求取子CCD在t时刻的焦面辐照度具体为：
设仿真三线阵相机的某一台CCD相机相对孔径为F，光学透过率为τ，则实际交点坐标的入瞳辐亮度L对应该CCD相机上子CCD焦面辐照度E为：
其中，α为仿真相机光学向量与相机主光轴的夹角；
对子CCD在t时刻的焦面辐照度滤波具体为通过光学点扩散函数对子CCD焦面辐照度滤波；
光学点扩散函数为：
PSF(x,y)＝exp(-(x2+y2)/(2σ2))；
式中，exp为指数函数；σ为函数下降参数，表示点扩散函数值随(x,y)变化的速度；(x,y)为点扩散函数作用范围内，即PSF(x,y)＞0内的一点相对点扩散函数原点的二维平面坐标；点扩散函数原点为每个子CCD中心点；
S14，将每个CCD内部的所有子CCD的滤波后的焦面辐照度取平均值，得到每个CCD相机在t时刻的焦面辐照度。
S2，对各个离散时刻的焦面辐照度在积分时间段上积分，和在级数上积分，得到CCD相机的焦面辐照度；
由于经过步骤S1获取的是CCD相机的一级CCD在一个时刻的离散平均焦面辐照度，因此，需要将积分时间段细分，获取每个细分时间段的中心时刻对应的离散平均焦面辐照度，将每一级CCD在各个离散时刻的离散平均焦面辐照度取平均值，作为CCD相机在一个积分时间段内的平均焦面辐照度；
接下来，对一个积分时间段内的平均焦面辐照度在级数上积分，设CCD相机的积分级数为M，则将每一级CCD在上述积分时间段内的平均焦面辐照度取平均值，作为仿真相机中多级CCD的焦面辐照度。
S3，将CCD相机的焦面辐照度转换为CCD相机的入瞳辐亮度；具体为根据以下公式：
S4，根据CCD相机的入瞳辐亮度计算积分时间段对应的行影像的影像灰度值；
S5，根据CCD相机的入瞳辐亮度求取信噪比，并根据信噪比计算积分时间段对应的行影像的噪声影像灰度值；
S6，将噪声影像灰度值与影像灰度值相加得到积分时间段对应的行影像的最终影像灰度值；
S7，返回S1继续执行，计算下一积分时间段对应的行影像的最终影像灰度值，直到计算完所有积分时间段对应的行影像的最终影像灰度值，根据所有行影像的最终影像灰度值输出所有行组成的整幅影像。
图3是本发明一种三线阵立体测绘卫星成像仿真系统的结构示意图。
如图3所示，一种三线阵立体测绘卫星成像仿真系统，包括：
焦面辐照度获取模块，用于获取CCD相机在积分时间段内各个离散时刻的焦面辐照度；
积分模块，用于对各个离散时刻的焦面辐照度在积分时间段上积分，和在级数上积分，得到CCD相机的焦面辐照度；
转换模块，用于将CCD相机的焦面辐照度转换为CCD相机的入瞳辐亮度；
影像灰度值计算模块，用于根据CCD相机的入瞳辐亮度计算积分时间段对应的行影像的影像灰度值；
噪声影像灰度值获取模块，用于根据CCD相机的入瞳辐亮度求取信噪比，并根据信噪比计算积分时间段对应的行影像的噪声影像灰度值；
最终影像灰度值获取模块，用于将噪声影像灰度值与影像灰度值相加得到积分时间段对应的行影像的最终影像灰度值；
循环及输出模块，用于循环调用焦面辐照度获取模块、积分模块、转换模块、影像灰度值计算模块、噪声影像灰度值获取模块和最终影像灰度值获取模块，计算下一积分时间段对应的行影像的最终影像灰度值，直到计算完所有积分时间段对应的行影像的最终影像灰度值，根据所有行影像的最终影像灰度值输出所有行组成的整幅影像。
焦面辐照度获取模块包括：
划分子模块，用于将CCD相机中每个CCD划分为多个子CCD，得到若干个子CCD；
入瞳辐亮度获取子模块，用于求取每个子CCD在离散时刻t的入瞳辐亮度；其中，离散时刻t为CCD相机在一个积分时间段中的某一时刻；
子CCD焦面辐照度获取子模块，用于根据子CCD在t时刻的入瞳辐亮度求取子CCD在t时刻的焦面辐照度并滤波；
焦面辐照度获取子模块，将每个CCD内部的所有子CCD的滤波后的焦面辐照度取平均值，得到每个CCD相机在t时刻的焦面辐照度。
入瞳辐亮度获取子模块包括：
剖分单元，用于将输入的几何模型剖分得到多个物方三角形，并对多个物方三角形编号；
投射单元，用于将每个物方三角形的三个顶点坐标分别投射到CCD相机上，得到多个像方三角形；
统计及记录单元，用于统计落在每个像方三角形内的像点坐标，并记录像方三角形内每个像点对应的多个物方三角形的编号；
转换单元，用于根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与物方三角形的交点；
实际交点获取单元，用于将距离相机光学节点最近的交点作为实际交点；
子CCD入瞳辐亮度获取单元，用于获取预先输入的三视Pleiades立体卫星影像中与CCD相机视角对应的影像上与实际交点对应位置处的影像灰度值，并反算为实际交点入瞳辐亮度，得到子CCD的入瞳辐亮度。
转换单元在根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与物方三角形的交点的步骤中，具体为根据卫星轨道运行参数、卫星姿态参数和相机安装模型计算每个像点对应的主光线向量，并转换到物方坐标系。
卫星轨道运行参数按照以下方法获取：
计算卫星轨道周期T；具体为：
式中，G为万有引力常数；M为地球质量；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
计算第1个离散时刻卫星原点对应的实际时间T0；具体为：
T0＝Tdsd-T/4；
式中，Tdsd为降交点时刻；
计算第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)；具体为：
θ＝(90-dθ×i)/180×π；
X＝R×cos(θ)；
Y＝0；
Z＝R×sin(θ)；
式中，θ为地标原点的地心角；(X,Y,Z)为第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)中任一组数据；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
将第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)绕X顺时针旋转I-90度后绕Z轴逆时针旋转L度，得到第一次旋转后的第i个离散时刻的卫星原点位置序列(X′i,Y′i,Z′i)；I为卫星轨道倾角，L为待拍摄目标区域中心的经度；
对第一次旋转后的第i个离散时刻的卫星原点位置序列(Xi′,Yi′,Zi′)绕Z轴旋转i×dL度，得到第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi″,Yi″,Zi″)；dL为0.1秒内地球的自转角度；dL＝(1/2400)°；
根据第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi″,Yi″,Zi″)中任一个卫星原点坐标(Xi,Yi,Zi)对应的初始大地坐标序列中(Li,Bi,Hi)的初始大地坐标(L,B,H)计算第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)；具体为：
L＝atan(Y/X)
H＝P/cos(B)-N；
式中，(L,B,H)为(Li,Bi,Hi)序列的任一组数据，(X,Y,Z)为对应(L,B,H)的(Xi,Yi,Zi)序列中的数据；
将第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)的经度Li加上ΔL，得到卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L′i,Bi,Hi)；其中，Li为第i个离散时刻的卫星原点的经度，Bi为第i个离散时刻的卫星原点的纬度，Hi为第i个离散时刻的卫星原点的高度；i为离散时刻的编号，i＝0,1,L N；Lc为卫星原点的星下点通过目标区域中心点的经度，Bc为卫星原点的星下点通过目标区域中心点的纬度；ΔL＝Lc-L+RandL；
RandL＝1.0*rand()/RAND_MAX*(SW-GW)/2/6371000/(2π)*360；
式中，GW为输入的几何模型宽度；SW为Pleiades立体卫星影像的地面幅宽；rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值；
将卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L′i,Bi,Hi)转换成卫星原点地心直角坐标序列(X″′i,Y″′i,Z″′i)，i＝0,1,L N，具体为：
式中，(L′,B,H)为(L′i,Bi,Hi)序列的任一组数据，(X″′,Y″′,Z″′)为对应(X″′i,Y″′i,Z″′i)序列的任一组数据；
将地标原点的地心角θ加上阈值，求取新的卫星原点地心直角坐标序列(X″″i,Y″″i,Z″″i)；
根据卫星原点地心直角坐标序列(X″′i,Y″′i,Z″′i)和新的卫星原点地心直角坐标序列(X″″i,Y″″i,Z″″i)计算卫星原点速度序列(Vxi,Vyi,Vzi)：
Vxi＝(X″″i-X″′i)/(T*0.00002/360)
Vyi＝(Y″″i-Y″′i)/(T*0.00002/360)
Vzi＝(Z″″i-Z″′i)/(T*0.00002/360)
式中，T为轨道周期；
根据第i个离散时刻卫星原点位置、速度和两者对应的实际时间Ti，得到卫星轨道运行参数(Ti,X′i,Y′i,Z′i,Vxi,Vyi,Vzi)；式中，Ti＝T0+i*0.1，T0为第1个离散时刻对应的实际时间。
卫星姿态参数的获取包括以下步骤：将t+Randt时刻的卫星姿态参数作为t时刻的卫星姿态参数，其中t表示第i个离散时刻对应的时间；式中，Randt表示时间随机变量，Randt＝1.0*rand()/RAND_MAX*(SW-GH)/2/6371000/(2π)*T；式中，T为轨道周期，GH为高度，SW为仿真影像地面幅宽，rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值。
以上仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (12)

1.一种三线阵立体测绘卫星成像仿真方法，其特征在于，包括：
S1，获取CCD相机在积分时间段内各个离散时刻的焦面辐照度；
S2，对各个离散时刻的焦面辐照度在积分时间段上积分，和在级数上积分，得到CCD相机的焦面辐照度；
S3，将CCD相机的焦面辐照度转换为CCD相机的入瞳辐亮度；
S4，根据CCD相机的入瞳辐亮度计算所述积分时间段对应的行影像的影像灰度值；
S5，根据CCD相机的入瞳辐亮度求取信噪比，并根据所述信噪比计算所述积分时间段对应的行影像的噪声影像灰度值；
S6，将所述噪声影像灰度值与所述影像灰度值相加得到所述积分时间段对应的行影像的最终影像灰度值；
S7，返回S1继续执行，计算下一积分时间段对应的行影像的最终影像灰度值，直到计算完所有积分时间段对应的行影像的最终影像灰度值，根据所有行影像的最终影像灰度值输出所有行组成的整幅影像。
2.根据权利要求1所述一种三线阵立体测绘卫星成像仿真方法，其特征在于，所述获取CCD相机在积分时间段内各个离散时刻的焦面辐照度，具体为：
将CCD相机中每个CCD划分为多个子CCD，得到若干个子CCD；
求取每个子CCD在离散时刻t的入瞳辐亮度；其中，离散时刻t为CCD相机在一个积分时间段中的某一时刻；
根据子CCD在t时刻的入瞳辐亮度求取子CCD在t时刻的焦面辐照度并滤波；
将每个CCD内部的所有子CCD的滤波后的焦面辐照度取平均值，得到每个CCD相机在t时刻的焦面辐照度。
3.根据权利要求2所述一种三线阵立体测绘卫星成像仿真方法，其特征在于，所述求取每个子CCD在离散时刻t的入瞳辐亮度，具体为：
将输入的几何模型剖分得到多个物方三角形，并对所述多个物方三角形编号；所述几何模型为Lidar数据；
将每个所述物方三角形的三个顶点坐标分别投射到CCD相机上，得到多个像方三角形；
统计落在每个像方三角形内的像点坐标，并记录所述像方三角形内每个像点对应的多个物方三角形的编号；
根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与所述物方三角形的交点；
将距离相机光学节点最近的交点作为实际交点；
获取预先输入的三视Pleiades立体卫星影像中与CCD相机视角对应的影像上与实际交点对应位置处的影像灰度值，并反算为实际交点入瞳辐亮度，得到子CCD的入瞳辐亮度。
4.根据权利要求3所述一种三线阵立体测绘卫星成像仿真方法，其特征在于，所述根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与所述物方三角形的交点的步骤中，具体为根据卫星轨道运行参数、卫星姿态参数和相机安装模型计算每个像点对应的主光线向量，并转换到物方坐标系。
5.根据权利要求4所述一种三线阵立体测绘卫星成像仿真方法，其特征在于，所述卫星轨道运行参数按照以下方法获取：
计算卫星轨道周期T；具体为：
式中，G为万有引力常数；M为地球质量；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
计算第1个离散时刻卫星原点对应的实际时间T0；具体为：
T0＝Tdsd-T/4；
式中，Tdsd为降交点时刻；
计算第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)；具体为：
θ＝(90-dθ×i)/180×π；
X＝R×cos(θ)；
Y＝0；
Z＝R×sin(θ)；
式中，θ为地标原点的地心角；(X,Y,Z)为第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)中任一组数据；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
将第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)绕X顺时针旋转I-90度后绕Z轴逆时针旋转L度，得到第一次旋转后的第i个离散时刻的卫星原点位置序列(Xi',Yi',Zi')；I为卫星轨道倾角，L为待拍摄目标区域中心的经度；
对第一次旋转后的第i个离散时刻的卫星原点位置序列(Xi',Yi',Zi')绕Z轴旋转i×dL度，得到第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi”,Yi”,Zi”)；dL为0.1秒内地球的自转角度；dL＝(1/2400)°；
根据第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi”,Yi”,Zi”)中任一个卫星原点坐标(Xi,Yi,Zi)对应的初始大地坐标序列中(Li,Bi,Hi)的初始大地坐标(L,B,H)计算第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)；具体为：
L＝atan(Y/X)
H＝P/cos(B)-N；
式中，(L,B,H)为(Li,Bi,Hi)序列的任一组数据，(X,Y,Z)为对应(L,B,H)的(Xi,Yi,Zi)序列中的数据；a＝6378137；
将第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)的经度Li加上ΔL，得到卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L'i,Bi,Hi)；其中，Li为第i个离散时刻的卫星原点的经度，Bi为第i个离散时刻的卫星原点的纬度，Hi为第i个离散时刻的卫星原点的高度；i为离散时刻的编号，i＝0,1,L N；Lc为卫星原点的星下点通过目标区域中心点的经度，Bc为卫星原点的星下点通过目标区域中心点的纬度；ΔL＝Lc-L+RandL；
RandL＝1.0*rand()/RAND_MAX*(SW-GW)/2/6371000/(2π)*360；
式中，GW为输入的几何模型宽度；SW为Pleiades立体卫星影像的地面幅宽；rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值；
将卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L'i,Bi,Hi)转换成卫星原点地心直角坐标序列(Xi”',Yi”',Zi”')，i＝0,1,L N，具体为：
式中，(L',B,H)为(L'i,Bi,Hi)序列的任一组数据，(X”',Y”',Z”')为对应(Xi”',Yi”',Zi”')序列的任一组数据；
将地标原点的地心角θ加上阈值，求取新的卫星原点地心直角坐标序列(Xi””,Yi””,Zi””)；
根据卫星原点地心直角坐标序列(Xi”',Yi”',Zi”')和新的卫星原点地心直角坐标序列(Xi””,Yi””,Zi””)计算卫星原点速度序列(Vxi,Vyi,Vzi)：
Vxi＝(Xi””-Xi”')/(T*0.00002/360)
Vyi＝(Yi””-Yi”')/(T*0.00002/360)
Vzi＝(Zi””-Zi”')/(T*0.00002/360)
式中，T为轨道周期；
根据第i个离散时刻卫星原点位置、速度和两者对应的实际时间Ti，得到卫星轨道运行参数(Ti,X′i,Yi',Z′i,Vxi,Vyi,Vzi)；式中，Ti＝T0+i*0.1，T0为第1个离散时刻对应的实际时间。
6.根据权利要求4所述一种三线阵立体测绘卫星成像仿真方法，其特征在于，卫星姿态参数的获取包括以下步骤：
将t+Randt时刻的卫星姿态参数作为t时刻的卫星姿态参数，其中t表示第i个离散时刻对应的时间；式中，Randt表示时间随机变量，Randt＝1.0*rand()/RAND_MAX*(SW-GH)/2/6371000/(2π)*T；
式中，T为轨道周期，GH为高度，SW为仿真影像地面幅宽，rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值。
7.一种三线阵立体测绘卫星成像仿真系统，其特征在于，包括：
焦面辐照度获取模块，用于获取CCD相机在积分时间段内各个离散时刻的焦面辐照度；
积分模块，用于对各个离散时刻的焦面辐照度在积分时间段上积分，和在级数上积分，得到CCD相机的焦面辐照度；
转换模块，用于将CCD相机的焦面辐照度转换为CCD相机的入瞳辐亮度；
影像灰度值计算模块，用于根据CCD相机的入瞳辐亮度计算所述积分时间段对应的行影像的影像灰度值；
噪声影像灰度值获取模块，用于根据CCD相机的入瞳辐亮度求取信噪比，并根据所述信噪比计算所述积分时间段对应的行影像的噪声影像灰度值；
最终影像灰度值获取模块，用于将所述噪声影像灰度值与所述影像灰度值相加得到所述积分时间段对应的行影像的最终影像灰度值；
循环及输出模块，用于循环调用焦面辐照度获取模块、积分模块、转换模块、影像灰度值计算模块、噪声影像灰度值获取模块和最终影像灰度值获取模块，计算下一积分时间段对应的行影像的最终影像灰度值，直到计算完所有积分时间段对应的行影像的最终影像灰度值，根据所有行影像的最终影像灰度值输出所有行组成的整幅影像。
8.根据权利要求7所述一种三线阵立体测绘卫星成像仿真系统，其特征在于，所述焦面辐照度获取模块包括：
划分子模块，用于将CCD相机中每个CCD划分为多个子CCD，得到若干个子CCD；
入瞳辐亮度获取子模块，用于求取每个子CCD在离散时刻t的入瞳辐亮度；其中，离散时刻t为CCD相机在一个积分时间段中的某一时刻；
子CCD焦面辐照度获取子模块，用于根据子CCD在t时刻的入瞳辐亮度求取子CCD在t时刻的焦面辐照度并滤波；
焦面辐照度获取子模块，将每个CCD内部的所有子CCD的滤波后的焦面辐照度取平均值，得到每个CCD相机在t时刻的焦面辐照度。
9.根据权利要求8所述一种三线阵立体测绘卫星成像仿真系统，其特征在于，所述入瞳辐亮度获取子模块包括：
剖分单元，用于将输入的几何模型剖分得到多个物方三角形，并对所述多个物方三角形编号；
投射单元，用于将每个所述物方三角形的三个顶点坐标分别投射到CCD相机上，得到多个像方三角形；
统计及记录单元，用于统计落在每个像方三角形内的像点坐标，并记录所述像方三角形内每个像点对应的多个物方三角形的编号；
转换单元，用于根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与所述物方三角形的交点；
实际交点获取单元，用于将距离相机光学节点最近的交点作为实际交点；
子CCD入瞳辐亮度获取单元，用于获取预先输入的三视Pleiades立体卫星影像中与CCD相机视角对应的影像上与实际交点对应位置处的影像灰度值，并反算为实际交点入瞳辐亮度，得到子CCD的入瞳辐亮度。
10.根据权利要求9所述一种三线阵立体测绘卫星成像仿真系统，其特征在于，所述转换单元在根据每个像点对应的多个物方三角形的编号，将每个像点对应的主光线向量转换到物方坐标系，得到多个与所述物方三角形的交点的步骤中，具体为根据卫星轨道运行参数、卫星姿态参数和相机安装模型计算每个像点对应的主光线向量，并转换到物方坐标系。
11.根据权利要求10所述一种三线阵立体测绘卫星成像仿真系统，其特征在于，所述卫星轨道运行参数按照以下方法获取：
计算卫星轨道周期T；具体为：
式中，G为万有引力常数；M为地球质量；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
计算第1个离散时刻卫星原点对应的实际时间T0；具体为：
T0＝Tdsd-T/4；
式中，Tdsd为降交点时刻；
计算第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)；具体为：
θ＝(90-dθ×i)/180×π；
X＝R×cos(θ)；
Y＝0；
Z＝R×sin(θ)；
式中，θ为地标原点的地心角；(X,Y,Z)为第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)中任一组数据；R为卫星测量坐标系各原点坐标到惯性参考系原点的距离；
将第i个离散时刻的卫星原点位置序列(Xi,Yi,Zi)绕X顺时针旋转I-90度后绕Z轴逆时针旋转L度，得到第一次旋转后的第i个离散时刻的卫星原点位置序列(Xi',Yi',Zi')；I为卫星轨道倾角，L为待拍摄目标区域中心的经度；
对第一次旋转后的第i个离散时刻的卫星原点位置序列(Xi',Yi',Zi')绕Z轴旋转i×dL度，得到第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi”,Yi”,Zi”)；dL为0.1秒内地球的自转角度；dL＝(1/2400)°；
根据第二次旋转后的第i个离散时刻的卫星原点坐标序列(Xi”,Yi”,Zi”)中任一个卫星原点坐标(Xi,Yi,Zi)对应的初始大地坐标序列中(Li,Bi,Hi)的初始大地坐标(L,B,H)计算第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)；具体为：
L＝atan(Y/X)
H＝P/cos(B)-N；
式中，(L,B,H)为(Li,Bi,Hi)序列的任一组数据，(X,Y,Z)为对应(L,B,H)的(Xi,Yi,Zi)序列中的数据；a＝6378137；
将第i个离散时刻的卫星原点初始大地坐标序列(Li,Bi,Hi)的经度Li加上ΔL，得到卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L'i,Bi,Hi)；其中，Li为第i个离散时刻的卫星原点的经度，Bi为第i个离散时刻的卫星原点的纬度，Hi为第i个离散时刻的卫星原点的高度；i为离散时刻的编号，i＝0,1,L N；Lc为卫星原点的星下点通过目标区域中心点的经度，Bc为卫星原点的星下点通过目标区域中心点的纬度；ΔL＝Lc-L+RandL；
RandL＝1.0*rand()/RAND_MAX*(SW-GW)/2/6371000/(2π)*360；
式中，GW为输入的几何模型宽度；SW为Pleiades立体卫星影像的地面幅宽；rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值；
将卫星原点的星下点通过目标区域中心点坐标(Lc,Bc)的卫星原点大地坐标序列(L'i,Bi,Hi)转换成卫星原点地心直角坐标序列(Xi”',Yi”',Zi”')，i＝0,1,L N，具体为：
式中，(L',B,H)为(L'i,Bi,Hi)序列的任一组数据，(X”',Y”',Z”')为对应(Xi”',Yi”',Zi”')序列的任一组数据；
将地标原点的地心角θ加上阈值，求取新的卫星原点地心直角坐标序列(Xi””,Yi””,Zi””)；
根据卫星原点地心直角坐标序列(Xi”',Yi”',Zi”')和新的卫星原点地心直角坐标序列(Xi””,Yi””,Zi””)计算卫星原点速度序列(Vxi,Vyi,Vzi)：
Vxi＝(Xi””-Xi”')/(T*0.00002/360)
Vyi＝(Yi””-Yi”')/(T*0.00002/360)
Vzi＝(Zi””-Zi”')/(T*0.00002/360)
式中，T为轨道周期；
根据第i个离散时刻卫星原点位置、速度和两者对应的实际时间Ti，得到卫星轨道运行参数(Ti,X′i,Yi',Z′i,Vxi,Vyi,Vzi)；式中，Ti＝T0+i*0.1，T0为第1个离散时刻对应的实际时间。
12.根据权利要求10所述一种三线阵立体测绘卫星成像仿真系统，其特征在于，卫星姿态参数的获取包括以下步骤：
将t+Randt时刻的卫星姿态参数作为t时刻的卫星姿态参数，其中t表示第i个离散时刻对应的时间；式中，Randt表示时间随机变量，Randt＝1.0*rand()/RAND_MAX*(SW-GH)/2/6371000/(2π)*T；
式中，T为轨道周期，GH为高度，SW为仿真影像地面幅宽，rand()为C++的均匀分布随机整数生成器生成的随机数，RAND_MAX为随机整数生成器生成的随机整数的上限值。
CN201610500333.XA 2016-06-29 2016-06-29 一种三线阵立体测绘卫星成像仿真方法和系统 Active CN106126839B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610500333.XA CN106126839B (zh) 2016-06-29 2016-06-29 一种三线阵立体测绘卫星成像仿真方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610500333.XA CN106126839B (zh) 2016-06-29 2016-06-29 一种三线阵立体测绘卫星成像仿真方法和系统

Publications (2)

Publication Number Publication Date
CN106126839A CN106126839A (zh) 2016-11-16
CN106126839B true CN106126839B (zh) 2019-07-09

Family

ID=57285552

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610500333.XA Active CN106126839B (zh) 2016-06-29 2016-06-29 一种三线阵立体测绘卫星成像仿真方法和系统

Country Status (1)

Country Link
CN (1) CN106126839B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108470099B (zh) * 2018-03-15 2021-11-02 长沙天玖卫星科技有限公司 光学成像类小卫星成像能力分析与姿态控制指标分析方法
CN108961319B (zh) * 2018-07-10 2021-11-19 中国科学院长春光学精密机械与物理研究所 双线阵tdi空间相机对动态飞机运动特性的分析方法
CN114820581B (zh) * 2022-05-26 2023-03-24 清华大学 轴对称光学成像并行仿真方法及装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103438900A (zh) * 2013-07-25 2013-12-11 航天恒星科技有限公司 三线阵相机影像协同绝对辐射定标和校正方法
CN103914808A (zh) * 2014-03-14 2014-07-09 国家测绘地理信息局卫星测绘应用中心 一种资源三号卫星三线阵影像和多光谱影像的拼接方法
CN105528500A (zh) * 2016-01-19 2016-04-27 国家测绘地理信息局卫星测绘应用中心 一种分米级星载tdi ccd立体测绘相机成像仿真方法和系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2455359C (en) * 2004-01-16 2013-01-08 Geotango International Corp. System, computer program and method for 3d object measurement, modeling and mapping from single imagery

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103438900A (zh) * 2013-07-25 2013-12-11 航天恒星科技有限公司 三线阵相机影像协同绝对辐射定标和校正方法
CN103914808A (zh) * 2014-03-14 2014-07-09 国家测绘地理信息局卫星测绘应用中心 一种资源三号卫星三线阵影像和多光谱影像的拼接方法
CN105528500A (zh) * 2016-01-19 2016-04-27 国家测绘地理信息局卫星测绘应用中心 一种分米级星载tdi ccd立体测绘相机成像仿真方法和系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
《亚m级卫星TDI CCD立体测绘相机成像仿真》;岳庆兴，唐新明，高小明;《武汉大学学报.信息科学版》;20150331;第40卷(第3期);第327-332页
《卫星三线阵TDI CCD相机成像仿真研究》;岳庆兴;《中国博士学位论文全文数据库》;20120415;第I140-27页

Also Published As

Publication number Publication date
CN106126839A (zh) 2016-11-16

Similar Documents

Publication Publication Date Title
CN105528500B (zh) 一种分米级星载tdi ccd立体测绘相机成像仿真方法和系统
Henriksen et al. Extracting accurate and precise topography from LROC narrow angle camera stereo observations
Bosch et al. A multiple view stereo benchmark for satellite imagery
US7944547B2 (en) Method and system of generating 3D images with airborne oblique/vertical imagery, GPS/IMU data, and LIDAR elevation data
CN103913148B (zh) 航天tdi ccd相机全链路数值仿真方法
CN106780712B (zh) 联合激光扫描和影像匹配的三维点云生成方法
CN104215261B (zh) 大视场反射式自由曲面空间相机畸变标定方法
CN101852623B (zh) 一种卫星光学遥感相机内方元素在轨检校方法
CN106803892A (zh) 一种基于光场测量的光场高清晰成像方法
CN105913435B (zh) 一种适用于大区域的多尺度遥感影像匹配方法及系统
CN107705329A (zh) 基于几何约束的高分辨率光学卫星凝视影像配准方法
CN109100719B (zh) 基于星载sar影像与光学影像的地形图联合测图方法
CN106126839B (zh) 一种三线阵立体测绘卫星成像仿真方法和系统
CN108982901A (zh) 一种匀速旋转体的转速测量方法
CN105444778B (zh) 一种基于成像几何反演的星敏感器在轨定姿误差获取方法
CN107451957B (zh) 一种星载tdi cmos相机成像仿真方法及设备
CN103129752A (zh) 一种基于地面导航的光学遥感卫星姿态角误差动态补偿方法
CN107067437A (zh) 一种基于多视几何和光束法平差的无人机定位系统及方法
CN103617649A (zh) 一种基于相机自标定技术的河工模型地形测量方法
CN109166143A (zh) 一种大区域网立体测绘卫星影像匹配方法
CN108195736A (zh) 一种三维激光点云提取植被冠层间隙率的方法
CN108562900B (zh) 一种基于高程校正的sar图像几何配准方法
CN105023281A (zh) 基于点扩散函数波前修正的星点像质心计算方法
JP6899915B2 (ja) 標高データグリッドのためのシャドウキャスティング
CN116740288B (zh) 一种融合激光雷达、倾斜摄影的三维重建方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant