CN110929427A - 一种遥感卫星视频成像快速仿真方法 - Google Patents

一种遥感卫星视频成像快速仿真方法 Download PDF

Info

Publication number
CN110929427A
CN110929427A CN201911338559.4A CN201911338559A CN110929427A CN 110929427 A CN110929427 A CN 110929427A CN 201911338559 A CN201911338559 A CN 201911338559A CN 110929427 A CN110929427 A CN 110929427A
Authority
CN
China
Prior art keywords
point
earth
ground
latitude
longitude
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.)
Pending
Application number
CN201911338559.4A
Other languages
English (en)
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.)
Jilin University
Original Assignee
Jilin University
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 Jilin University filed Critical Jilin University
Priority to CN201911338559.4A priority Critical patent/CN110929427A/zh
Publication of CN110929427A publication Critical patent/CN110929427A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4007—Scaling of whole images or parts thereof, e.g. expanding or contracting based on interpolation, e.g. bilinear interpolation

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种遥感卫星视频成像快速仿真方法，涉及遥感卫星视频成像技术，主要包括：输入仿真目标点经纬度
Figure DEST_PATH_IMAGE002
；
Figure DEST_PATH_IMAGE004
表示地面凝视目标点经纬度；根据输入经纬度
Figure DEST_PATH_IMAGE002A
，利用STK辅助生成卫星轨道位置；根据轨道位置，进行对目标点视频成像过程姿态规划；根据已有的位置、姿态信息，进行像物映射；根据求得的经纬度
Figure DEST_PATH_IMAGE007
，对应高程图，得到此点对应的高程gao_h；采用二维拉格朗日插值方法进行地球曲面重构，对像平面上各个像素点分别定位得到各个点的灰度值，完成整个影像仿真过程。本发明大大减少了计算量，缩短了仿真时间，提高了仿真速度。在整体系统中对计算机性能要求可降低，为高帧频大幅宽视频实时仿真提供了可能性。

Description

一种遥感卫星视频成像快速仿真方法
技术领域
本发明涉及遥感卫星视频成像技术，具体涉及一种遥感卫星视频成像快速仿真方法。
背景技术
目前，随着航空航天技术的发展，使用遥感卫星对地观测越来越普遍，凝视成像地面快速仿真技术是实现在轨应用的重要前提和必要环节。为了直观的分析和评测凝视成像图像质量、便于后续数据分析，成像过程理论建模和仿真分析在遥感器的全生命周期内，发挥着十分重要的作用。近年来，在国防需求的牵引下，国内越来越多的高等院校和科研院所都参与到遥感成像仿真领域的研究中，并取得了不小的进展。现有技术根据应用需求设计载荷物理参数指标，根据卫星运动轨迹及成像目标地理位置计算视频凝视成像时刻姿态机动角度，通过分析能量在大气、场景和遥感器之间传递时发生的各种相互作用，基于全物理或半物理平台，生成理想视频凝视仿真图像。
现有仿真技术根据应用需求设计载荷物理参数指标，根据卫星平台设计指标，仿真理想状态下卫星运动轨迹并进行摄动因素改正，根据卫星运动轨迹及成像目标地理位置自动计算卫星凝视成像时刻的姿态激动角度，并进行偏流角修正。根据卫星运动轨迹、姿态数据及相机设计参数等,分析能量在大气、场景和遥感器之间传递时发生的各种相互作用，生成理想成像仿真图像。同时现有技术部分借助全物理或半物理仿真平台，如气浮台等模拟卫星本体作为控制对象，验证仿真过程中姿态控制系统方案设计的正确性与可靠性。
现有仿真技术主要针对于理论仿真，不具备大范围面阵、视频仿真模式。而视频成像仿真具有大幅宽、多帧频特点，不可避免的出现数据计算量大、所需内存大等特点，对计算机速度、内存要求较高。按照现有技术方案进行视频仿真，仿真时间过长，难以满足工程需求。
发明内容
本发明针对上述问题，提供了一种遥感卫星视频成像快速仿真方法。
本发明采用的技术方案是：一种遥感卫星视频成像快速仿真方法，包括以下步骤：
一种遥感卫星视频成像快速仿真方法，其特征在于，包括以下步骤：
S1，输入仿真目标点经纬度
Figure 316813DEST_PATH_IMAGE001
；
Figure 776745DEST_PATH_IMAGE002
表示地面凝视目标点经纬度；
S2，根据输入的经纬度
Figure 147683DEST_PATH_IMAGE003
，利用STK辅助生成卫星轨道位置；输入
轨道参数及起止时间，根据视频帧频设置时间间隔，使用可视模块分析选取目标点对卫星可观测时间段内的轨道位置，导出并保存文件；
S3，根据轨道位置，对目标点视频成像过程进行姿态规划；
S4，根据已有的位置、姿态信息，进行像物映射；
S5，根据步骤S4求得的经纬度
Figure 34868DEST_PATH_IMAGE004
，对应高程图，得到此点对应的
高程gao_h；
Figure 858467DEST_PATH_IMAGE005
表示未对应高程信息时得到的相机视线与地球交点的经纬度；gao_h表示步骤S4求得的点
Figure 122089DEST_PATH_IMAGE006
在高程图中对应的此点的高程；
S6，根据步骤S5得到的gao_h，再次进行步骤4，求得
Figure 81955DEST_PATH_IMAGE007
，再次进
行步骤S5，重复此过程直到满足求得的高程不再明显变化，此时的经纬度为最终视线与地球交点
Figure 264675DEST_PATH_IMAGE008
；
Figure 215051DEST_PATH_IMAGE009
为对应高程信息后得到的相机视线与地球交点的经纬度；
Figure 141419DEST_PATH_IMAGE010
为对应高程信息并经过多次迭代后满足迭代终止要求的相机视线与地球交点的经纬度；
S7，采用二维拉格朗日插值方法进行地球曲面重构，像素行列号
Figure 831157DEST_PATH_IMAGE011
作
为已知，地面经纬度
Figure 450357DEST_PATH_IMAGE012
作为节点坐标；
Figure 123915DEST_PATH_IMAGE013
表示像素点的位置信息，x、y均为整数；
S8，按照步骤S7的方法计算每个像素点对应的地面点经纬度，对应原始图像，定位地面物点；
S9，根据物点在原始图像中对应点的行列号进行灰度提取，作为像面像素点
Figure 588395DEST_PATH_IMAGE014
的灰度值；
Figure 522853DEST_PATH_IMAGE014
表示位置
Figure 188320DEST_PATH_IMAGE015
处的像素p；
S10，重复步骤S8至S9，对像平面上各个像素点分别定位得到各个点的灰度值，完成整个影像仿真过程。
进一步地，所述步骤S3具体包括：
定义
Figure 473808DEST_PATH_IMAGE016
为地心惯性系到航天器轨道系的旋转矩阵，其中轨道升交点赤
经为
Figure 853230DEST_PATH_IMAGE017
，轨道倾角为
Figure 907774DEST_PATH_IMAGE018
，升交点幅角为
Figure 744143DEST_PATH_IMAGE019
；从惯性系转换到轨道系按照如下顺序进行四次转换：绕z轴旋转
Figure 251348DEST_PATH_IMAGE020
、绕x轴旋转
Figure 198575DEST_PATH_IMAGE021
、绕z轴旋转
Figure 107625DEST_PATH_IMAGE022
，最后左乘坐标轴反向矩阵，计算方法如下：
Figure 973950DEST_PATH_IMAGE023
。
定义
Figure 109396DEST_PATH_IMAGE024
为地心地固系到地心惯性系的旋转矩阵，按照如下公式计算
Figure 453790DEST_PATH_IMAGE024
：
Figure 856827DEST_PATH_IMAGE025
(2)。
其中
Figure 894053DEST_PATH_IMAGE026
分别为岁差矩阵、章动矩阵、真恒星时旋转矩阵、
极移矩阵；
真恒星时旋转矩阵计算方法如(4)，其中
Figure 516796DEST_PATH_IMAGE027
为格林威治真时间时：
Figure 930459DEST_PATH_IMAGE028
。
定义相机光轴指向与本体坐标系Z轴重合，影像中心像方单位矢量为
Figure 548523DEST_PATH_IMAGE029
；对地凝视期间，相机光轴一直稳定对准地面目标，即
Figure 897596DEST_PATH_IMAGE030
始终与航天器位置指向地面目标的矢量
Figure 601109DEST_PATH_IMAGE031
重合：
Figure 428251DEST_PATH_IMAGE032
(4)。
Figure 431979DEST_PATH_IMAGE033
为航天器在地固系下位置，
Figure 545429DEST_PATH_IMAGE034
为地面目标点D在地固系下的位置；
Figure 378649DEST_PATH_IMAGE035
Figure 134115DEST_PATH_IMAGE036
(5)。
其中
Figure 602137DEST_PATH_IMAGE037
为轨道坐标系到地心地固系的旋转矩阵，
Figure 417646DEST_PATH_IMAGE038
Figure 236697DEST_PATH_IMAGE039
分别为绕x、y、z轴的姿态角，
Figure 264696DEST_PATH_IMAGE040
从本体系转到轨道系下的向量，定义
Figure 977437DEST_PATH_IMAGE041
整理式(5)可得：
Figure 839214DEST_PATH_IMAGE042
。
由此可计算横滚角和俯仰角为：
Figure 4616DEST_PATH_IMAGE043
(7) 。
更进一步地，所述步骤S4具体包括：
求得地面点后，由已有高分辨率图像反演地面点反射率，进行相应像
点灰度求解；对于面阵CCD，本体坐标系下的像元观测向量写为如下形式：
Figure 475787DEST_PATH_IMAGE044
。
Figure 43034DEST_PATH_IMAGE045
为像面上像素行列值，
Figure 934767DEST_PATH_IMAGE046
为CCD中心像元行列值，
Figure 728411DEST_PATH_IMAGE047
为CCD
像元尺寸，
Figure 363791DEST_PATH_IMAGE048
为探测器焦距；
将
Figure 660912DEST_PATH_IMAGE049
由卫星本体坐标系转换到地心地固系，需要进行本体-轨道、轨道-惯性、惯性-地固3次转换，公式如下：
Figure 989125DEST_PATH_IMAGE050
（9）。
Figure 863540DEST_PATH_IMAGE051
为地心惯性系转换到地心地固系的旋转矩阵，
Figure 443557DEST_PATH_IMAGE052
为轨道坐标系转换到地心惯性系的旋转矩阵，
Figure 985397DEST_PATH_IMAGE053
为按照xyz转序的欧拉旋转矩阵的逆矩阵；
求得像点
Figure 861342DEST_PATH_IMAGE054
在地固系下的观测向量
Figure 223053DEST_PATH_IMAGE055
后，计算
Figure 75603DEST_PATH_IMAGE055
与地球椭球面交点；
假设地面点坐标为M(X,Y,Z)，得到如下方程：
Figure 737528DEST_PATH_IMAGE056
（10）。
其中P(XP,YP,ZP)为卫星在地固系下的位置坐标，由卫星辅助数据给
出，
Figure 141965DEST_PATH_IMAGE057
为比例系数；另外，由于M点在地球椭球上，M点坐标同时满足地球椭球方程：
Figure 131918DEST_PATH_IMAGE058
（11）。
其中
Figure 647213DEST_PATH_IMAGE059
a、b分别为地球赤道半径和极半径，h为地面高程，初始时刻可设为0；
联立式(10)、(11)，求解未知数
Figure 39011DEST_PATH_IMAGE060
，当有两个解时保留较小值；将
Figure 879928DEST_PATH_IMAGE060
代
回式(10)，得到地面点M的地固系坐标(X,Y,Z)，进而转换为地理坐标
Figure 216231DEST_PATH_IMAGE061
。
更进一步地，所述步骤S7具体包括：
设已知
Figure 174698DEST_PATH_IMAGE062
，其中
Figure 280057DEST_PATH_IMAGE063
其中m,n为总行
列数，即选取m*n个像素点，按照步骤S6中方法计算得到对应的地面点经
Figure 167241DEST_PATH_IMAGE064
，则二维拉格朗日多项式为：
Figure 990841DEST_PATH_IMAGE065
(14)。
点
Figure 113518DEST_PATH_IMAGE066
为已知地面经纬度的像素的行列号，
Figure 214329DEST_PATH_IMAGE067
对应地物的
经纬度，
Figure 397049DEST_PATH_IMAGE068
为待求像素的行列号；
选取不同个数插值点，以式（14）进行二维插值，得到每点对应的地
面经纬度。
更进一步地，所述步骤S8具体包括：
按照步骤S7的方法计算每个像素点对应的地面点经纬度，对应原始图像，定位地面物点；
Figure 583310DEST_PATH_IMAGE069
Figure 509678DEST_PATH_IMAGE070
分别为物点
Figure 324050DEST_PATH_IMAGE071
在原始图像中对应点行列号，
Figure 331800DEST_PATH_IMAGE072
分别为原始图像初始位置经纬度，
Figure 129992DEST_PATH_IMAGE073
分别为原始图像像素点尺度，即每个像素点覆盖的经纬度范围。
本发明的优点：
本发明采用二维拉格朗日插值方法结合传统物像映射模型，采用部分插值方法代替逐点计算，大大减少了计算量，缩短了仿真时间，提高了仿真速度。在整体系统中对计算机性能要求可降低，为高帧频大幅宽视频实时仿真提供了可能性。
除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
图1是本发明的观测方向与地球椭球交点图；
图2是本发明的仿真过程框图；
图3是本发明的仿真实验90s内三轴姿态角变化图；
图4是本发明的仿真实验光轴指向点与目标点距离图；
图5是本发明的仿真实验凝视成像仿真图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
参考图1和图2，如图1和图2所示，一种遥感卫星视频成像快速仿真方法，包括以下步骤：
一种遥感卫星视频成像快速仿真方法，其特征在于，包括以下步骤：
S1，输入仿真目标点经纬度
Figure 204258DEST_PATH_IMAGE001
；
Figure 404296DEST_PATH_IMAGE002
表示地面凝视目标点经纬度；
S2，根据输入的经纬度
Figure 69763DEST_PATH_IMAGE003
，利用STK辅助生成卫星轨道位置；输入
轨道参数及起止时间，根据视频帧频设置时间间隔，使用可视模块分析选取目标点对卫星可观测时间段内的轨道位置，导出并保存文件；
S3，根据轨道位置，对目标点视频成像过程进行姿态规划；
S4，根据已有的位置、姿态信息，进行像物映射；
S5，根据步骤S4求得的经纬度
Figure 824093DEST_PATH_IMAGE004
，对应高程图，得到此点对应的
高程gao_h；
Figure 357842DEST_PATH_IMAGE005
表示未对应高程信息时得到的相机视线与地球交点的经纬度；gao_h表示步骤S4求得的点
Figure 287752DEST_PATH_IMAGE006
在高程图中对应的此点的高程；
S6，根据步骤S5得到的gao_h，再次进行步骤4，求得
Figure 983176DEST_PATH_IMAGE007
，再次进
行步骤S5，重复此过程直到满足求得的高程不再明显变化，此时的经纬度为最终视线与地球交点
Figure 129861DEST_PATH_IMAGE008
；
Figure 936143DEST_PATH_IMAGE009
为对应高程信息后得到的相机视线与地球交点的经纬度；
Figure 720559DEST_PATH_IMAGE010
为对应高程信息并经过多次迭代后满足迭代终止要求的相机视线与地球交点的经纬度；
S7，采用二维拉格朗日插值方法进行地球曲面重构，像素行列号
Figure 852463DEST_PATH_IMAGE011
作
为已知，地面经纬度
Figure 581385DEST_PATH_IMAGE012
作为节点坐标；
Figure 332303DEST_PATH_IMAGE013
表示像素点的位置信息，x、y均为整数；
S8，按照步骤S7的方法计算每个像素点对应的地面点经纬度，对应原
始图像，定位地面物点；
S9，根据物点在原始图像中对应点的行列号进行灰度提取，作为像面
像素点
Figure 361439DEST_PATH_IMAGE014
的灰度值；
Figure 8452DEST_PATH_IMAGE014
表示位置
Figure 755828DEST_PATH_IMAGE015
处的像素p；
S10，重复步骤S8至S9，对像平面上各个像素点分别定位得到各个点
的灰度值，完成整个影像仿真过程。
所述步骤S3具体包括：
定义
Figure 546323DEST_PATH_IMAGE016
为地心惯性系到航天器轨道系的旋转矩阵，其中轨道升交点赤
经为
Figure 164386DEST_PATH_IMAGE017
，轨道倾角为
Figure 372514DEST_PATH_IMAGE018
，升交点幅角为
Figure 216973DEST_PATH_IMAGE019
；从惯性系转换到轨道系按照如下顺序进行四次转换：绕z轴旋转
Figure 168749DEST_PATH_IMAGE020
、绕x轴旋转
Figure 782264DEST_PATH_IMAGE021
、绕z轴旋转
Figure 426872DEST_PATH_IMAGE022
，最后左乘坐标轴反向矩阵，计算方法如下：
Figure 617682DEST_PATH_IMAGE023
。
定义
Figure 248514DEST_PATH_IMAGE024
为地心地固系到地心惯性系的旋转矩阵，按照如下公式计算
Figure 106749DEST_PATH_IMAGE024
：
Figure 296159DEST_PATH_IMAGE025
(2)。
其中
Figure 708686DEST_PATH_IMAGE026
分别为岁差矩阵、章动矩阵、真恒星时旋转矩阵、
极移矩阵；
真恒星时旋转矩阵计算方法如(4)，其中
Figure 143210DEST_PATH_IMAGE027
为格林威治真时间时：
Figure 855951DEST_PATH_IMAGE028
。
定义相机光轴指向与本体坐标系Z轴重合，影像中心像方单位矢量为
Figure 717728DEST_PATH_IMAGE029
；对地凝视期间，相机光轴一直稳定对准地面目标，即
Figure 883130DEST_PATH_IMAGE030
始终与航天器位置指向地面目标的矢量
Figure 980399DEST_PATH_IMAGE031
重合：
Figure 423012DEST_PATH_IMAGE032
(4)。
Figure 314745DEST_PATH_IMAGE033
为航天器在地固系下位置，
Figure 344274DEST_PATH_IMAGE034
为地面目标点D在地固系下的位置；
Figure 245234DEST_PATH_IMAGE035
Figure 807934DEST_PATH_IMAGE036
(5)。
其中
Figure 136147DEST_PATH_IMAGE037
为轨道坐标系到地心地固系的旋转矩阵，
Figure 10562DEST_PATH_IMAGE038
Figure 325000DEST_PATH_IMAGE039
分别为绕x、y、z轴的姿态角，
Figure 866840DEST_PATH_IMAGE040
从本体系转到轨道系下的向量，定义
Figure 975741DEST_PATH_IMAGE041
整理式(5)可得：
Figure 603031DEST_PATH_IMAGE042
。
由此可计算横滚角和俯仰角为：
Figure 219695DEST_PATH_IMAGE043
(7) 。
所述步骤S4具体包括：
传感器的成像几何模型反映的是地面点的三维空间坐标与相应像点的像方坐标之间的数学关系。进行图像仿真必须建立可靠的成像几何模型获取每个像素点对应的地面点地固系坐标，即求解任意像点
Figure 881621DEST_PATH_IMAGE074
与投影中心连线与地球椭球的交点。求得地面点后，由已有高分辨率图像反演地面点反射率，进行相应像点灰度求解；对于面阵CCD，本体坐标系下的像元观测向量写为如下形式：
Figure 286057DEST_PATH_IMAGE044
。
Figure 10431DEST_PATH_IMAGE045
为像面上像素行列值，
Figure 56884DEST_PATH_IMAGE046
为CCD中心像元行列值，
Figure 448682DEST_PATH_IMAGE047
为CCD
像元尺寸，
Figure 24020DEST_PATH_IMAGE048
为探测器焦距；
将
Figure 235690DEST_PATH_IMAGE049
由卫星本体坐标系转换到地心地固系，需要进行本体-轨道、轨道-惯性、惯性-地固3次转换，公式如下：
Figure 820255DEST_PATH_IMAGE050
（9）。
Figure 925614DEST_PATH_IMAGE051
为地心惯性系转换到地心地固系的旋转矩阵，
Figure 314264DEST_PATH_IMAGE052
为轨道坐标系转
换到地心惯性系的旋转矩阵，
Figure 403442DEST_PATH_IMAGE053
为按照xyz转序的欧拉旋转矩阵的逆矩阵；
求得像点
Figure 401485DEST_PATH_IMAGE054
在地固系下的观测向量
Figure 626930DEST_PATH_IMAGE055
后，计算
Figure 685016DEST_PATH_IMAGE055
与地球椭球面交点；
假设地面点坐标为M(X,Y,Z)，得到如下方程：
Figure 995912DEST_PATH_IMAGE056
（10）。
其中P(XP,YP,ZP)为卫星在地固系下的位置坐标，由卫星辅助数据给
出，
Figure 797646DEST_PATH_IMAGE057
为比例系数；另外，由于M点在地球椭球上，M点坐标同时满足地球椭球方程：
Figure 877597DEST_PATH_IMAGE058
（11）。
其中
Figure 231218DEST_PATH_IMAGE059
a、b分别为地球赤道半径和极半径，h为地面
高程，初始时刻可设为0；
联立式(10)、(11)，求解未知数
Figure 403311DEST_PATH_IMAGE060
，当有两个解时保留较小值；将
Figure 867791DEST_PATH_IMAGE060
代
回式(10)，得到地面点M的地固系坐标(X,Y,Z)，进而转换为地理坐标
Figure 943194DEST_PATH_IMAGE061
。
所述步骤S7具体包括：
设已知
Figure 467716DEST_PATH_IMAGE062
，其中
Figure 628570DEST_PATH_IMAGE063
其中m,n为总行
列数，即选取m*n个像素点，按照步骤S6中方法计算得到对应的地面点经
Figure 896741DEST_PATH_IMAGE064
，则二维拉格朗日多项式为：
Figure 951284DEST_PATH_IMAGE065
(14)。
点
Figure 522074DEST_PATH_IMAGE066
为已知地面经纬度的像素的行列号，
Figure 29279DEST_PATH_IMAGE067
对应地物的
经纬度，
Figure 466252DEST_PATH_IMAGE068
为待求像素的行列号；
选取不同个数插值点，以式（14）进行二维插值，得到每点对应的地
面经纬度。
所述步骤S8具体包括：
按照步骤S7的方法计算每个像素点对应的地面点经纬度，对应原始图像，定位地面物点；
Figure 375302DEST_PATH_IMAGE069
Figure 382573DEST_PATH_IMAGE070
分别为物点
Figure 111494DEST_PATH_IMAGE071
在原始图像中对应点行列号，
Figure 987046DEST_PATH_IMAGE072
分别为原始图像初始位置经纬度，
Figure 625969DEST_PATH_IMAGE073
分别为原始图像像素点尺度，即每个像素点覆盖的经纬度范围。
仿真实验：
本发明选取一地形平坦区域作为目标区域，选取区域内东经120°，北纬36°作为目标点进行视频成像时卫星三轴姿态计算及成像仿真，仿真参数如表1。仿真系统整体框图如图2，首先确定地面目标点，利用STK（Satellite Tool Kit，卫星工具包）辅助建立卫星轨道，使目标点对卫星可观测。然后按照本发明方法计算凝视过程所需的横滚角、俯仰角，并利用STK 辅助验证姿态计算的准确性。计算偏航角，修正由于轨道运动、地球自转等引起的图像旋转。然后建立成像几何模型作为直接法，引入二维拉格朗日插值算法作为快速算法，分别计算两种方法的地面点经纬度，对比像素误差和用时，最后生成仿真图像。
表1 设计实例仿真参数
Figure 663195DEST_PATH_IMAGE075
图3所示为仿真实验的90s内姿态角变化。由图3可知，90s凝视过程（0s是凝视开始时刻）中卫星主要依靠俯仰角的大角度姿态机动能力实现凝视观测，横滚角变化很小。为修正目标区域图像慢旋，偏航角也有较大变化。
为验证规划的姿态角准确性，将姿态角转化为四元数后倒入STK，得到光轴指向的地面点。计算地面点与目标点距离，结果如图4所示，在初始和结束位置距目标点较远，但均小于4e-003m，在中间位置距目标点最近。这是由于中间位置星地距离也最小，俯仰角约为0。图4结果表明所求姿态能够满足凝视要求，能以较小的误差指向目标点。
采用本文的姿态数据、几何模型与快速算法进行成像仿真实验。电脑采用win7 64位系统，2.60GHz处理器，4.00G安装内存。快速算法分别选取不同个数的插值点，进行二维拉格朗日插值，以典型成像几何模型直接（表2中简称直接法）求解结果为标准，分析单帧仿真时间及像素误差，结果如表2所示：
表2 快速仿真像素误差及时间
Figure 20359DEST_PATH_IMAGE076
根据表2可得，随着插值点的增加，像素误差会随之减小，但同时运算时间也会增加。以不超过0.3个像素误差为衡量标准，考虑到时间成本，选取16个插值点的快速算法为最优方法。由于无需逐点迭代高程，此时的运算时间为14.879s，仅占直接法时间的12.549%，减少了约87%的仿真时间。
通过上述方法得到的凝视仿真图像如图5所示。图5中的（a）（b）（c）分别为第0s、45s、90s、时的计算机仿真图像。分析姿态角可知第45s时俯仰角最接近0，此时卫星位于距离目标最近的轨道位置，因此图（b）看到的是正视图像，视场最小；而（a）（c）俯仰角不为0，光轴与地面夹角小于90°，视场比正视图像大，除目标中心点外能看到更多地物，但同时地面分辨率也会减小。由图可以看出，除视场大小有差别，三幅图的观测场景基本重合，表明通过本文规划的三轴姿态能够实现目标区域无旋控制，采用二维插值结合经典成像几何模型能够快速生成目标区域的计算机仿真图像。
本发明针对现有技术方案计算量大，仿真时间长，对计算机性能要求高等缺点，在现有技术基础上，依据工程项目实际所需，针对地形平坦区域进行大幅宽视频成像仿真，提出传统仿真模型结合拉格朗日插值算法，进行遥感影像快速成像仿真，降低计算机性能要求。
本发明主要技术点在于将传统仿真算法与二维插值算法相结合，达到快速仿真的目的。
本发明采用二维拉格朗日插值方法结合传统物像映射模型，采用部分插值方法代替逐点计算，大大减少了计算量，缩短了仿真时间，提高了仿真速度。在整体系统中对计算机性能要求可降低，为高帧频大幅宽视频实时仿真提供了可能性。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种遥感卫星视频成像快速仿真方法，其特征在于，包括以下步骤：
S1，输入仿真目标点经纬度
Figure 282986DEST_PATH_IMAGE001
；
Figure 322486DEST_PATH_IMAGE002
表示地面凝视目标点经纬度；
S2，根据输入的经纬度
Figure 656516DEST_PATH_IMAGE003
，利用STK辅助生成卫星轨道位置；输入
轨道参数及起止时间，根据视频帧频设置时间间隔，使用可视模块分析选取目标点对卫星可观测时间段内的轨道位置，导出并保存文件；
S3，根据轨道位置，对目标点视频成像过程进行姿态规划；
S4，根据已有的位置、姿态信息，进行像物映射；
S5，根据步骤S4求得的经纬度
Figure 641045DEST_PATH_IMAGE004
，对应高程图，得到此点对应的
高程gao_h；
Figure 693315DEST_PATH_IMAGE005
表示未对应高程信息时得到的相机视线与地球交点的经纬度；gao_h表示步骤S4求得的点
Figure 208610DEST_PATH_IMAGE006
在高程图中对应的此点的高程；
S6，根据步骤S5得到的gao_h，再次进行步骤4，求得
Figure 787359DEST_PATH_IMAGE007
，再次进
行步骤S5，重复此过程直到满足求得的高程不再明显变化，此时的经纬度为最终视线与地球交点
Figure 565959DEST_PATH_IMAGE008
；
Figure 715312DEST_PATH_IMAGE009
为对应高程信息后得到的相机视线与地球交点的经纬度；
Figure 237560DEST_PATH_IMAGE010
为对应高程信息并经过多次迭代后满足迭代终止要求的相机视线与地球交点的经纬度；
S7，采用二维拉格朗日插值方法进行地球曲面重构，像素行列号
Figure 670815DEST_PATH_IMAGE011
作
为已知，地面经纬度
Figure 620317DEST_PATH_IMAGE012
作为节点坐标；
Figure 178337DEST_PATH_IMAGE013
表示像素点的位置信息，x、y均为整数；
S8，按照步骤S7的方法计算每个像素点对应的地面点经纬度，对应原
始图像，定位地面物点；
S9，根据物点在原始图像中对应点的行列号进行灰度提取，作为像面
像素点
Figure 878178DEST_PATH_IMAGE014
的灰度值；
Figure 41306DEST_PATH_IMAGE014
表示位置
Figure 20763DEST_PATH_IMAGE015
处的像素p；
S10，重复步骤S8至S9，对像平面上各个像素点分别定位得到各个点
的灰度值，完成整个影像仿真过程。
2.根据权利要求1所述的遥感卫星视频成像快速仿真方法，其特征在
于，所述步骤S3具体包括：
定义
Figure 534921DEST_PATH_IMAGE016
为地心惯性系到航天器轨道系的旋转矩阵，其中轨道升交点赤
经为
Figure 274338DEST_PATH_IMAGE017
，轨道倾角为
Figure 291973DEST_PATH_IMAGE018
，升交点幅角为
Figure 645594DEST_PATH_IMAGE019
；从惯性系转换到轨道系按照如下顺序进行四次转换：绕z轴旋转
Figure 506102DEST_PATH_IMAGE020
、绕x轴旋转
Figure 173844DEST_PATH_IMAGE021
、绕z轴旋转
Figure 157237DEST_PATH_IMAGE022
，最后左乘坐标轴反向矩阵，计算方法如下：
Figure 885021DEST_PATH_IMAGE023
定义
Figure 232826DEST_PATH_IMAGE024
为地心地固系到地心惯性系的旋转矩阵，按照如下公式计算
Figure 235417DEST_PATH_IMAGE024
：
Figure 227644DEST_PATH_IMAGE025
(2)
其中
Figure 1696DEST_PATH_IMAGE026
分别为岁差矩阵、章动矩阵、真恒星时旋转矩阵、
极移矩阵；
真恒星时旋转矩阵计算方法如(4)，其中
Figure 181005DEST_PATH_IMAGE027
为格林威治真时间时：
Figure 315183DEST_PATH_IMAGE028
定义相机光轴指向与本体坐标系Z轴重合，影像中心像方单位矢量为
Figure 427495DEST_PATH_IMAGE029
；对地凝视期间，相机光轴一直稳定对准地面目标，即
Figure 605404DEST_PATH_IMAGE030
始终与航天器位置指向地面目标的矢量
Figure 599905DEST_PATH_IMAGE031
重合：
Figure 147561DEST_PATH_IMAGE032
(4)
Figure 239014DEST_PATH_IMAGE033
为航天器在地固系下位置，
Figure 213923DEST_PATH_IMAGE034
为地面目标点D在地固系下的位置；
Figure 774349DEST_PATH_IMAGE035
Figure 860116DEST_PATH_IMAGE036
(5)
其中
Figure 806076DEST_PATH_IMAGE037
为轨道坐标系到地心地固系的旋转矩阵，
Figure 217466DEST_PATH_IMAGE038
Figure 501073DEST_PATH_IMAGE039
分别为绕x、y、z轴的姿态角，
Figure 187269DEST_PATH_IMAGE040
从本体系转到轨道系下的向量，定义
Figure 128680DEST_PATH_IMAGE041
整理式(5)可得：
Figure 570026DEST_PATH_IMAGE042
由此可计算横滚角和俯仰角为：
Figure 964098DEST_PATH_IMAGE043
(7) 。
3.根据权利要求1所述的遥感卫星视频成像快速仿真方法，其特征在
于，所述步骤S4具体包括：
求得地面点后，由已有高分辨率图像反演地面点反射率，进行相应像
点灰度求解；对于面阵CCD，本体坐标系下的像元观测向量写为如下形式：
Figure 267035DEST_PATH_IMAGE044
Figure 62952DEST_PATH_IMAGE045
为像面上像素行列值，
Figure 940778DEST_PATH_IMAGE046
为CCD中心像元行列值，
Figure 618884DEST_PATH_IMAGE047
为CCD
像元尺寸，
Figure 850146DEST_PATH_IMAGE048
为探测器焦距；
将
Figure 874471DEST_PATH_IMAGE049
由卫星本体坐标系转换到地心地固系，需要进行本体-轨道、轨道-惯性、惯性-地固3次转换，公式如下：
Figure 532986DEST_PATH_IMAGE050
（9）
Figure 26284DEST_PATH_IMAGE051
为地心惯性系转换到地心地固系的旋转矩阵，
Figure 61236DEST_PATH_IMAGE052
为轨道坐标系转
换到地心惯性系的旋转矩阵，
Figure 362904DEST_PATH_IMAGE053
为按照xyz转序的欧拉旋转矩阵的逆矩阵；
求得像点
Figure 333265DEST_PATH_IMAGE054
在地固系下的观测向量
Figure 189226DEST_PATH_IMAGE055
后，计算
Figure 621344DEST_PATH_IMAGE055
与地球椭球面交点；
假设地面点坐标为M(X,Y,Z)，得到如下方程：
Figure 246361DEST_PATH_IMAGE056
（10）
其中P(XP,YP,ZP)为卫星在地固系下的位置坐标，由卫星辅助数据给
出，
Figure 877369DEST_PATH_IMAGE057
为比例系数；另外，由于M点在地球椭球上，M点坐标同时满足地球椭球方程：
Figure 751784DEST_PATH_IMAGE058
（11）
其中
Figure 394118DEST_PATH_IMAGE059
a、b分别为地球赤道半径和极半径，h为地面
高程，初始时刻可设为0；
联立式(10)、(11)，求解未知数
Figure 998275DEST_PATH_IMAGE060
，当有两个解时保留较小值；将
Figure 169493DEST_PATH_IMAGE060
代
回式(10)，得到地面点M的地固系坐标(X,Y,Z)，进而转换为地理坐标
Figure 609833DEST_PATH_IMAGE061
。
4.根据权利要求1所述的遥感卫星视频成像快速仿真方法，其特征在
于，所述步骤S7具体包括：
设已知
Figure 524699DEST_PATH_IMAGE062
，其中
Figure 248942DEST_PATH_IMAGE063
其中m,n为总行
列数，即选取m*n个像素点，按照步骤S6中方法计算得到对应的地面点经
Figure 653378DEST_PATH_IMAGE064
，则二维拉格朗日多项式为：
Figure 705648DEST_PATH_IMAGE065
(14)
点
Figure 798107DEST_PATH_IMAGE066
为已知地面经纬度的像素的行列号，
Figure 252222DEST_PATH_IMAGE067
对应地物的
经纬度，
Figure 889876DEST_PATH_IMAGE068
为待求像素的行列号；
选取不同个数插值点，以式（14）进行二维插值，得到每点对应的地
面经纬度。
5.根据权利要求1所述的遥感卫星视频成像快速仿真方法，其特征在
于，所述步骤S8具体包括：
按照步骤S7的方法计算每个像素点对应的地面点经纬度，对应原始图像，定位地面物点；
Figure 429442DEST_PATH_IMAGE069
Figure 827057DEST_PATH_IMAGE070
分别为物点
Figure 135678DEST_PATH_IMAGE071
在原始图像中对应点行列号，
Figure 209813DEST_PATH_IMAGE072
分别为原始图像初始位置经纬度，
Figure 767834DEST_PATH_IMAGE073
分别为原始图像像素点尺度，即每个像素点覆盖的经纬度范围。
CN201911338559.4A 2019-12-23 2019-12-23 一种遥感卫星视频成像快速仿真方法 Pending CN110929427A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911338559.4A CN110929427A (zh) 2019-12-23 2019-12-23 一种遥感卫星视频成像快速仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911338559.4A CN110929427A (zh) 2019-12-23 2019-12-23 一种遥感卫星视频成像快速仿真方法

Publications (1)

Publication Number Publication Date
CN110929427A true CN110929427A (zh) 2020-03-27

Family

ID=69860677

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911338559.4A Pending CN110929427A (zh) 2019-12-23 2019-12-23 一种遥感卫星视频成像快速仿真方法

Country Status (1)

Country Link
CN (1) CN110929427A (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111781854A (zh) * 2020-07-02 2020-10-16 北京世冠金洋科技发展有限公司 一种天体仿真的运行控制方法、装置及电子设备
CN112498745A (zh) * 2020-11-16 2021-03-16 长光卫星技术有限公司 一种低轨光学卫星立体成像自主规划方法
CN113094371A (zh) * 2021-04-14 2021-07-09 嘉兴毕格智能科技有限公司 一种用户自定义坐标系的实现方法
CN113612927A (zh) * 2021-10-08 2021-11-05 广东海纬地恒空间信息技术有限公司 一种基于摄影测量与遥感技术的数码成像系统
CN114152267A (zh) * 2021-02-26 2022-03-08 武汉大学 一种火星轨道相机影像仿真方法及系统
CN114440885A (zh) * 2021-12-24 2022-05-06 中国人民解放军战略支援部队信息工程大学 一种静止轨道遥感卫星定位方法及装置
CN114578385A (zh) * 2022-01-18 2022-06-03 自然资源部第二海洋研究所 一种基于cocts四元逐点摆扫方式的几何定位方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102147249A (zh) * 2010-02-08 2011-08-10 同济大学 基于直线特征的星载光学线阵影像精确纠正处理方法
US20120293669A1 (en) * 2011-04-25 2012-11-22 Skybox Imaging, Inc. Systems and methods for overhead imaging and video
CN103679711A (zh) * 2013-11-29 2014-03-26 航天恒星科技有限公司 一种遥感卫星线阵推扫光学相机在轨外方位参数标定方法
CN105136164A (zh) * 2015-08-13 2015-12-09 航天恒星科技有限公司 考虑卫星综合运动的凝视成像仿真及质量评估方法及装置
CN108268735A (zh) * 2018-01-29 2018-07-10 浙江大学 基于多源遥感卫星融合数据的地表土壤水分降尺度方法
US20190333643A1 (en) * 2018-04-26 2019-10-31 Vektor Medical, Inc. Machine learning using clinical and simulated data
CN110598834A (zh) * 2019-09-19 2019-12-20 吉林大学 一种双目视觉探测系统结构优化方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102147249A (zh) * 2010-02-08 2011-08-10 同济大学 基于直线特征的星载光学线阵影像精确纠正处理方法
US20120293669A1 (en) * 2011-04-25 2012-11-22 Skybox Imaging, Inc. Systems and methods for overhead imaging and video
CN103679711A (zh) * 2013-11-29 2014-03-26 航天恒星科技有限公司 一种遥感卫星线阵推扫光学相机在轨外方位参数标定方法
CN105136164A (zh) * 2015-08-13 2015-12-09 航天恒星科技有限公司 考虑卫星综合运动的凝视成像仿真及质量评估方法及装置
CN108268735A (zh) * 2018-01-29 2018-07-10 浙江大学 基于多源遥感卫星融合数据的地表土壤水分降尺度方法
US20190333643A1 (en) * 2018-04-26 2019-10-31 Vektor Medical, Inc. Machine learning using clinical and simulated data
CN110598834A (zh) * 2019-09-19 2019-12-20 吉林大学 一种双目视觉探测系统结构优化方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张刘等: ""光学遥感卫星凝视成像姿态规划及快速仿真方法"", 《吉林大学学报（工学版）》 *
杨蕊等: ""遥感面阵凝视图像并行超分辨重建方法"", 《红外与激光工程》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111781854A (zh) * 2020-07-02 2020-10-16 北京世冠金洋科技发展有限公司 一种天体仿真的运行控制方法、装置及电子设备
CN112498745A (zh) * 2020-11-16 2021-03-16 长光卫星技术有限公司 一种低轨光学卫星立体成像自主规划方法
CN112498745B (zh) * 2020-11-16 2022-06-21 长光卫星技术股份有限公司 一种低轨光学卫星立体成像自主规划方法
CN114152267A (zh) * 2021-02-26 2022-03-08 武汉大学 一种火星轨道相机影像仿真方法及系统
CN113094371A (zh) * 2021-04-14 2021-07-09 嘉兴毕格智能科技有限公司 一种用户自定义坐标系的实现方法
CN113612927A (zh) * 2021-10-08 2021-11-05 广东海纬地恒空间信息技术有限公司 一种基于摄影测量与遥感技术的数码成像系统
CN113612927B (zh) * 2021-10-08 2022-02-11 广东海纬地恒空间信息技术有限公司 一种基于摄影测量与遥感技术的数码成像系统
CN114440885A (zh) * 2021-12-24 2022-05-06 中国人民解放军战略支援部队信息工程大学 一种静止轨道遥感卫星定位方法及装置
CN114440885B (zh) * 2021-12-24 2023-10-13 中国人民解放军战略支援部队信息工程大学 一种静止轨道遥感卫星定位方法及装置
CN114578385A (zh) * 2022-01-18 2022-06-03 自然资源部第二海洋研究所 一种基于cocts四元逐点摆扫方式的几何定位方法
CN114578385B (zh) * 2022-01-18 2024-05-14 自然资源部第二海洋研究所 一种基于cocts四元逐点摆扫方式的几何定位方法

Similar Documents

Publication Publication Date Title
CN110929427A (zh) 一种遥感卫星视频成像快速仿真方法
CN104848860B (zh) 一种敏捷卫星成像过程姿态机动规划方法
CN102565797B (zh) 一种针对聚束模式星载sar图像的几何校正方法
CN102654576B (zh) 基于sar图像和dem数据的图像配准方法
CN106124170A (zh) 一种基于高精度姿态信息的相机光轴指向计算方法
CN109948852B (zh) 一种敏捷卫星的同轨多点目标成像任务规划方法
CN101114022A (zh) 无姿态信息条件下的航空多光谱扫描仪几何粗校正方法
CN101750067B (zh) 一种成像式地球敏感器地球扁率修正方法
CN106909161B (zh) 一种敏捷卫星零偏流角成像的最优姿态机动规划方法
CN102426025B (zh) 遥感卫星姿态机动时偏流修正角的仿真分析方法
CN105444778B (zh) 一种基于成像几何反演的星敏感器在轨定姿误差获取方法
CN105136164B (zh) 考虑卫星综合运动的凝视成像仿真及质量评估方法及装置
CN102519433B (zh) 一种利用rpc反演星载线阵传感器几何定标参数方法
CN105928525A (zh) 一种卫星对月定标的姿态确定方法
CN108344396A (zh) 一种敏捷卫星斜条带成像模式姿态计算方法
CN102346922A (zh) 航天遥感载荷成像几何畸变三维可视化仿真方法
CN107450582A (zh) 一种基于星上实时规划的相控阵数传引导控制方法
CN107132542A (zh) 一种基于光学和多普勒雷达的小天体软着陆自主导航方法
CN103644918A (zh) 卫星对月探测数据定位处理方法
CN106251282A (zh) 一种机械臂采样环境仿真图的生成方法及装置
CN105444781A (zh) 星载自主引导成像地面验证方法
CN113264201A (zh) 一种敏捷卫星主动推扫成像实现方法
CN106897962B (zh) 大视场空间对地观测图像的圆盘投影和拼接方法
CN104976991B (zh) 一种考虑卫星姿态变化的三线阵影像像方偏差的获取方法
CN108489483B (zh) 一种船载星光定向仪单星次优修正算法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20200327

RJ01 Rejection of invention patent application after publication