本发明涉及一种面向侧视中分辨率卫星的几何成像构建方法，包括以下步骤：步骤1，导入卫星影像元数据以及卫星星历数据，确定卫星轨道及姿态拟合方程以及数据链节点的参数初值以及其他参数初值；步骤2，根据基准已包含投影信息的影像元数据自动初步获得图像控制点初值；步骤3，手工设定至少三个控制点，通过灰度匹配的方法自动匹配图像控制点；然后自动提出控制点缺失区域，对缺失区域进行控制点分布改化；步骤4，将完成步骤3的控制点分布改化,并进行坐标变换，并构建卫星的几何成像模型；步骤5，根据步骤1中已经获取的参数初值和步骤4中的卫星的几何成像模型进行模型参数改化迭代解算；步骤6，输出经纠正后的模型结果及精度指标。
1.一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，包括以下步骤：步骤1，导入卫星影像元数据以及卫星星历数据，确定卫星轨道及姿态拟合方程以及数据链节点的参数初值以及其他参数初值；
步骤2，根据基准已包含投影信息的影像元数据自动初步获得图像控制点初值；
步骤3，手工设定至少三个控制点，通过灰度匹配的方法自动匹配图像控制点；然后根据设定的步骤1已经设定的参数初值所需的条件，依据卫星轨道及姿态拟合数据节点位置自动检测图像控制点分布并自动提出控制点缺失区域，对缺失区域进行控制点分布改化；
步骤4，将完成步骤3的控制点分布改化,并进行坐标变换，即由投影坐标转到地固地心三维直角坐标系，再结合卫星轨道拟合数据规划到近似卫星本体坐标系下Porb=Morb(t)·PECR，并构建卫星的几何成像模型；
步骤5，根据步骤1中已经获取的参数初值和步骤4中的卫星的几何成像模型进行模型参数改化迭代解算；
2.根据权利要求1所述的一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，所述的步骤1中，方程与数据链以及其他参数初值的确定包括以下步骤：步骤1.1，三则二次多项式分别拟合卫星在地固地心坐标系下的位置矢量变化，其中参数X0,X1,X2,Y0,Y1,Y2,Z0,Z1,Z2通过元数据中的数据节点依最小二乘计算；
步骤1.2，三则二次多项式分别拟合卫星在地固地心坐标系下的速度矢量变化，其中参数VX0,VX1,VX2,VY0,VY1,VY2,VZ0,VZ1,VZ2通过元数据中的数据节点依最小二乘计算；
步骤1.3，三条数据链，每个结点包括卫星t时刻在本体坐标系下的姿态角以及该姿态角变化速率，每条数据链包括7～9个数据结点，三个方向的旋角数据节点可以分别表述为：对t时刻某方向姿态角的插值按照两点三次Hermite插值计算，对于λ角，在t时刻的插值计算公式为：步骤1.4，将卫星标称的传感器安放倾角作为偏移常量的初值，并且按照标称的星下分辨率r与标称的卫星飞行高度h计算数字传感器的等效焦距初值:
3.根据权利要求1所述的一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，所述步骤2中，手工设定部分控制点即手工选取3-6个控制点，分别位于待纠正影像中央以及四个角附近。
4.根据权利要求1所述的一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，所述的步骤3中，进行控制点分布改化即进行有目的的对步骤3中的控制点缺失区域进行控制点添加，具体操作步骤如下：步骤3.1，建立二分图B={P,C,E}，C为控制点节点集合，每个控制点占两个节点，P为待求参数集合，每个参数占一个节点，若控制点对应节点ci与参数pj有关联，则ci与pj间有一条边ek，否则无边；其中共有控制点个数n，参数个数m，0≤i≤2n,0≤j≤m；则二部图B中其P中的点只与相应时间段内成像范围内的C中的点之间有边；
步骤3.2，依据匈牙利算法搜索二部图B的最大匹配；
步骤3.3，标记步骤3.2中生成最大匹配的结果，其中P集中没有找到匹配的参数点需要对应的增加控制点，标记该类点的控制范围，优先在被多次重复标记的区域内添加一个控制点，刷新二分图的结构重复执行标记步骤3.2中生成最大匹配的结果，直到找到P到C的完全相容匹配。
5.根据权利要求1所述的一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，所述步骤4中，进行控制点坐标变换，构建卫星的几何成像模型原理可被描述为：卫星本体坐标系下指向已知地面点的向量经过空间相似变换之后转至像方坐标系，其与像方坐标系下指向已知点对应图像点的向量共线，其方程如下：其 中 ：
式中角度值 由步骤1.3中插值得到，g(sample)为畸变模型，本步骤采用五次多项式用于拟合线阵传感器变形以及折光等因素造成的畸变：g(x)=g0+g1x+g2x2+g3x3+g4x4+g5x5。
6.根据权利要求1所述的一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，所述的步骤5中，进行模型参数改化迭代解算包括：步骤6.1，建立误差方程，构建线性系统用于迭代求解；
步骤6.2，规定卫星几何模型参数的解算顺序，循环迭代求解；
7.根据权利要求6所述的一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，所述的步骤6.1的具体操作步骤如下：步骤6.11，根据步骤4中所述成像模型，建立修正模型参数的原始误差方程式为：定义G(sample)=sample+g(sample)，则G′(sample)为G(sample)的一阶导数，即G′
(x)=1+g1+2g2x+3g3x+4g4x+5g5x方程的系数分别为：
q7=fb1+fb3；q8=-fcosκt；q9=-G(sample)；
则误差方程可表示为：V=AΔ-L,则直接使用最小二乘解算改正数其法方程为：式中 为Δ的估值；
步骤6.12，参数解算中当单次求解多个参数时使用普修正解法增强解算的可靠性，以避免法方程病态造成的无法收敛或收敛到错误值，即在改正数迭代中单次迭代又引入一层内层迭代用以解算法方程，内层迭代改写法方程为： 式中为第i次计算的改正数，I为单位阵，该解法为最小二乘的无偏估计。
8.根据权利要求6所述的一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，所述的步骤6.2的具体操作步骤如下：步骤6.21，解算等效焦距f粗略值；
步骤6.22，解算角度偏移常量 ω0,κ0粗略值；
步骤6.23，解算畸变模型系数g0,g1,g2,g3,g4,g5粗略值，同时改化角度偏移常量ω0,κ0；
步骤6.24，重复执行步骤6.21至步骤6.23至控制点残差减小速度小于10-2；
步骤6.25，修正所有角度数据链节点数据 ωi,Vωi,κi,Vκi；
步骤6.26，精化计算畸变模型系数g0,g1,g2,g3,g4,g5至点位残差不再变小。
9.根据权利要求6所述的一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，所述的步骤6.3的具体操作步骤如下：步骤6.31，使用二次多项式通用模型计算点位正算的概略图像位置，解算方程如下式：步骤6.32，用概略位置的扫描线坐标对应时间，t=rt·line(rt为扫描线的扫描速率)，根据模型参数中的数据链内插瞬时外方位元素 ωt,κt,Xst,YstXst,VXt,VYt,VZt；
步骤6.33，用初始的瞬时外方位元素代入成像模型，根据line坐标的残差值迭代改化扫描线位置，计算公式如下：linei=linei-1+Vline， 分别为转入卫星本体坐标系的地面点坐标中的两个坐标分量；
步骤6.34，用精化扫描时刻的瞬时外方位元素代入成像畸变模型，迭代求出采样点列坐标，从而完成模型机精确正算，畸变差模型改正数计算公式如下：samplei=samplei-1+Vsample。
一种面向侧视中分辨率卫星的几何成像构建方法
技术领域
[0001] 本发明属于影像处理领域，涉及一种面向侧视中分辨率卫星的几何成像构建方法。
背景技术
[0002] 自然灾害发生的时间和地点总是有极大的随机性，随着现代人类社会的进步，快速及时的对自然环境灾害进行响应并提供相应的援助是减少灾害损失的最有效手段，因而对区域性自然灾害的监控变得尤为重要。建立监测体系的核心方式是建立一套反应时间迅速的遥感卫星系统，在成本较低的小卫星星座成为主流发展方向的今天，使用少数几颗大视场角并以侧视观测为成像方式的中分辨率小卫星星座作为环境减灾监控体系的信息获取方式成为既保证响应速度又节省成本的最佳选择。而遥感数据在供使用之前需要赋予一定的地理参考信息，即通过特定的方式对其进行几何处理，几何处理的核心几何成像模型的建立。
[0003] 架设于小卫星上的中分辨率宽视场角侧视观测的线阵传感器采集到的影像有以下几点主要特征：扫描范围大，受地球曲率影像严重，由于采用侧视观测这一影响更加严重；视场角大，在单条扫描线上完全以中心投影方式呈现，而与平行投影相差较大；卫星平台稳定性一般，对于成本较低的小卫星，由于各种因素卫星姿态抖动较严重，无法用平行投影进行拟合。而现阶段对中低分辨率影像的几何处理手段主要依赖经典的多项式拟合方式，即进行像点对像点的单一方式进行配准。但是采用线阵传感器的大视场角的中分辨率卫星，由于覆盖面积大，地形以及地球曲率对变形畸变的影响，用单一的多项式拟合精度的非常差，不能满足大多数情况的制图需求。同时其也不适用仿射变换模拟其成像过程，该模型的应用条件限制在长焦距，小视场角的高分辨率卫星；又由于卫星平台本身因素，姿态抖动模式较为复杂，并且卫星起算数据精度有限，无法用直接线性变换模型（DLT）或有理多项式模型（RFM）作为几何模型描述成像过程。描述宽视场角成像一般只能靠传统的共线方程模型，但是描述线阵传感器的共线方程对传感器外定向参数的变化没有明确定义，并且对参数变化的描述和答解又涉及到了对控制点分布等一系列数据要求。
发明内容
[0004] 本发明主要是解决现有技术所存在的技术问题；提供了一种充分利用卫星元数据同时考虑到地球曲率，地形起伏以及传感器形变，卫星平台抖动等综合因素造成的畸变，提高了几何处理的精度，并提出了相应算法保证解算数据的完整性的一种面向侧视中分辨率卫星的几何成像构建方法。
[0005] 本发明的上述技术问题主要是通过下述技术方案得以解决的：
[0006] 一种面向侧视中分辨率卫星的几何成像构建方法，其特征在于，包括以下步骤：
[0007] 步骤1，导入卫星影像元数据以及卫星星历数据，确定卫星轨道及姿态拟合方程以及数据链节点的参数初值以及其他参数初值；
[0008] 步骤2，根据基准已包含投影信息的影像元数据自动初步获得图像控制点初值；
[0009] 步骤3，手工设定至少三个控制点，通过灰度匹配的方法自动匹配图像控制点；然后根据设定的步骤1已经设定的参数初值所需的条件，依据卫星轨道及姿态拟合数据节点位置自动检测图像控制点分布并自动提出控制点缺失区域，对缺失区域进行控制点分布改化；
[0010] 步骤4，将完成步骤3的控制点分布改化,并进行坐标变换，即由投影坐标转到地固地心三维直角坐标系，再结合卫星轨道拟合数据规划到近似卫星本体坐标系下Porb=Morb(t)·PECR，并构建卫星的几何成像模型；
[0011] 步骤5，根据步骤1中已经获取的参数初值和步骤4中的卫星的几何成像模型进行模型参数改化迭代解算；
[0012] 步骤6，输出经纠正后的模型结果及精度指标。
[0013] 在上述的一种面向侧视中分辨率卫星的几何成像构建方法，所述的步骤1中，方程与数据链以及其他参数初值的确定包括以下步骤：
[0014] 步骤1.1，三则二次多项式分别拟合卫星在地固地心坐标系下的位置矢量变化，其中参数X0,X1,X2,Y0,Y1,Y2,Z0,Z1,Z2通过元数据中的数据节点依最小二乘计算；
[0015] 步骤1.2，三则二次多项式分别拟合卫星在地固地心坐标系下的速度矢量变化，其中参数VX0,VX1,VX2,VY0,VY1,VY2,VZ0,VZ1,VZ2通过元数据中的数据节点依最小二乘计算；
[0016] 步骤1.3，三条数据链，每个结点包括卫星该时刻在本体坐标系下的姿态角以及该方位姿态角变化速率，每条数据链包括7～9个数据结点，三个方向的旋角数据节点可以分别表述为：
[0017] 对t时刻某方向姿态角的插值按照两点三次Hermite插值计算，对于λ角，在t时刻的插值计算公式为：
[0018]
[0019]
[0020] 步骤1.4，将卫星标称的传感器安放倾角作为偏移常量的初值，并且按照标称的星下分辨率r与标称的卫星飞行高度h计算数字传感器的等效焦距初值:
[0021] 在上述的一种面向侧视中分辨率卫星的几何成像构建方法，所述步骤2中，手工设定部分控制点即手工选取3-6个控制点，分别位于待纠正影像中央以及四个角附近。
[0022] 在上述的一种面向侧视中分辨率卫星的几何成像构建方法，所述的步骤3中，进行控制点分布改化即进行有目的的对步骤3中的控制点缺失区域进行控制点添加，具体操作步骤如下：
[0023] 步骤3.1，建立二分图B={P,C,E}，C为控制点节点集合，每个控制点占两个节点，P为待求参数集合，每个参数占一个节点，若控制点对应节点ci与参数pj有关联，则ci与pj间有一条边ek，否则无边；其中共有控制点个数n，参数个数m，0≤i≤2n,0≤j≤m；则二部图B中其P中的点只与相应时间段内成像范围内的C中的点之间有边；
[0024] 步骤3.2，依据匈牙利算法搜索二部图B的最大匹配；
[0025] 步骤3.3，标记步骤3.2中生成最大匹配的结果，其中P集中没有找到匹配的参数点需要对应的增加控制点，标记将该类点的控制范围，优先在被多次重复标记的区域内添加一个控制点，刷新二分图的结构重复执行标记步骤3.2中生成最大匹配的结果，直到找到P到C的完全相容匹配。
[0026] 在上述的一种面向侧视中分辨率卫星的几何成像构建方法，所述步骤4中，进行控制点坐标变换，构建卫星的几何成像模型原理可被描述为：卫星本体坐标系下指向已知地面点的向量经过空间相似变换之后转至像方坐标系，其与像方坐标系下指向已知点对应图像点的向量共线，其方程如下：
[0027]
[0028] MRotation
[0029] 其 中 ：
[0030]
[0031] 式中角度值 由步骤1.3中所述数据链插值得到，g(sample)为畸
[0032] 变模型，本步骤采用五次多项式用于拟合线阵传感器变形以及折光等因素造成的畸变：g(x)=g0+g1x+g2x2+g3x3+g4x4+g5x5。
[0033] 在上述的一种面向侧视中分辨率卫星的几何成像构建方法，所述的步骤5中，进行模型参数改化迭代解算包括：
[0034] 步骤6.1，建立误差方程，构建线性系统用于迭代求解；
[0035] 步骤6.2，规定卫星几何模型参数的解算顺序，循环迭代求解；
[0036] 步骤6.3，根据模型以及参数进行控制点正算。
[0037] 在上述的一种面向侧视中分辨率卫星的几何成像构建方法，所述的步骤6.1的具体操作步骤如下：
[0038] 步骤6.11，根据步骤4中所述成像模型，建立修正模型参数的原始误差方程式为：
[0039]
[0040] 定义G(sample)=sample+g(sample)，则G′(sample)为G(sample)的一阶导数，2 3 4
即G′(x)=1+g1+2g2x+3g3x+4g4x+5g5x
[0041] 方程的系数分别为：
[0042]
[0043]
[0044]
[0045]
[0046]
[0047]
[0048]
[0049]
[0050]
[0051]
[0052] p12=0
[0053]
[0054]
[0055] p15=0
[0056] q0=q1=q2=q3=q4=q5=q6=0；
[0057] q7=fb1+fb3；q8=-fcosκt；q9=-G(sample)；
[0058]
[0059]
[0060]
[0061]
[0062]
[0063]
[0064] 为表达简便，令
[0065]
[0066]
[0067]
[0068] 则误差方程可表示为：V=AΔ-L,则直接使用最小二乘解算改正数其法方程为：式中 为Δ的估值；
[0069] 步骤6.12，参数解算中当单次求解多个参数时使用普修正解法增强解算的可靠性，以避免法方程病态造成的无法收敛或收敛到错误值，即在改正数迭代中单次迭代又引入一层内层迭代用以解算法方程，内层迭代改写法方程为： 式中 为第i次计算的改正数，I为单位阵，该解法为最小二乘的无偏估计。
[0070] 在上述的一种面向侧视中分辨率卫星的几何成像构建方法，所述的步骤6.2的具体操作步骤如下：
[0071] 步骤6.21，解算等效焦距f粗略值；
[0072] 步骤6.22，解算角度偏移常量 ,ω0,κ0粗略值；
[0073] 步骤6.23，解算畸变模型系数g0,g1,g2,g3,g4,g5粗略值，同时改化角度偏移常量,ω0,κ0；
[0074] 步骤6.24，重复执行步骤6.21至步骤6.23至控制点残差减小速度小于10-2；
[0075] 步骤6.25，修正所有角度数据链节点数据 ,ωi,Vωi,κi,Vκi；
[0076] 步骤6.26，精化计算畸变模型系数g0,g1,g2,g3,g4,g5至点位残差不再变小。
[0077] 在上述的一种面向侧视中分辨率卫星的几何成像构建方法，所述的步骤6.3的具体操作步骤如下：
[0078] 步骤6.31，使用二次多项式通用模型计算点位正算的概略图像位置，解算方程如下式，其中sample，line，x，y是控制点数据，为已知数，si，li（i=0,1,2，…，4）为该模型待求参数：
[0079]
[0080] 步骤6.32，用概略位置的扫描线坐标对应时间，t=rt·line(rt为扫描线的扫描速率)，根据模型参数中的数据链内插瞬时外方位元素 ωt，κt，Xst,YstXst,VXt,VYt,VZt；
[0081] 步骤6.33，用初始的瞬时外方位元素代入成像模型，根据line坐标的残差值迭代改化扫描线位置，计算公式如下：
[0082] linei=linei-1+Vline， 分别为转入卫星本体坐标系的地面点坐标中的两个坐标分量；
[0083] 步骤6.34，用精化扫描时刻的瞬时外方位元素代入成像畸变模型，迭代求出采样点列坐标，从而完成模型机精确正算，畸变差模型改正数计算公式如下，其中x0为已获得的控制点地理坐标值横坐标值：
[0084]
[0085] 因此，本发明具有如下优点：充分利用卫星元数据同时考虑到地球曲率，地形起伏以及传感器形变，卫星平台抖动等综合因素造成的畸变，提高了几何处理的精度，并提出了相应算法保证解算数据的完整性。
附图说明
[0086] 图1是本发明流程示意图。
[0087] 图2是宽幅卫星单个成像示意图。
[0088] 图3是影像图像控制点初始分布示意图。
[0089] 图4是本发明所指模型输出的控制点残差图(其中残差放大10倍)。
[0090] 图5是多项式模型输出的控制点残差图(其中残差放大10倍)。
具体实施方式
[0091] 下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。
[0092] 实施例：
[0093] 本实例通过带有卫星局部轨道定姿与定位信息的宽幅侧视卫星影像以及其与已知地理参考影像，在卫星传感器安放矩阵未知的前提下求得本发明所提到的几何模型的精确参数，从而为几何精纠正纠正提供准确的映射关系。其流程见图1。
[0094] 本实例选取30m分辨率(12000×12000)环境减灾小卫星1级数据（仅做过辐射纠正，未作任何几何处理。其成像方式为典型的宽幅侧视对扫成像，如图2所示，其中S为投影中心，I为承影的扫描线，I’为扫描线对应地面线，O为单线上的像主点，O’为像主点对应地面点，T为卫星对到，T’为星下轨迹），并参考带有UTM投影信息的Landsat7影像以及对应区域SRTM格式的DEM数据作为参考，使用本发明提到的模型构建方法及通用多项式模型进行对比试验。
[0095] 步骤1，导入卫星影像元数据以及卫星星历数据，确定卫星轨道及姿态拟合方程以及数据链节点的参数初值以及其他参数初值；
[0096] （a）xml格式元数据文件，其中包含数据获取的起始终止时间，根据行数计算扫描速率；
[0097] （b）读取EPH格式的GPS节点数据，用二次多项式分别拟合卫星位置变化以及速率变化，拟合方程表述为：
[0098]
[0099] 以及:
[0100]
[0101] （c）读取ATT格式的姿态角节点数据，将数据直接增加至姿态拟合的数据链中。
[0102] （d）根据标称分辨率和标称卫星航高计算等效焦距初值。
[0103] 步骤2，根据基准已包含投影信息的影像元数据自动初步获得图像控制点初值；
[0104] 步骤3，手工设定至少三个控制点，通过灰度匹配的方法自动匹配图像控制点；然后根据设定的步骤1已经设定的参数初值所需的条件，依据卫星轨道及姿态拟合数据节点位置自动检测图像控制点分布并自动提出控制点缺失区域，对缺失区域进行控制点分布改化,在本实施例中，通过灰度匹配的方法自动匹配图像控制点，匹配图像控制点结果共358个点，然后根据设定的参数所需的条件，建立表示参数和控制点位置关系的二分图，并根据二分图进行最大匹配的解给出需要增设控制点的位置范围以及个数，并手工增加少量控制点，修正后控制点分布如图3示。该控制点可用于通用的二次多项式模型进行解算，解算结果图像点位残差如图4所示。
[0105] 步骤4，将完成步骤3的控制点分布改化,并进行坐标变换，即由投影坐标转到地固地心三维直角坐标系，再结合卫星轨道拟合数据规划到近似卫星本体坐标系下Porb=Morb(t)·PECR，其中PECR为通过高斯投影逆运算获得的控制点地固地心直角坐标系坐标，Moeb为（d）中所述的转换矩阵，Porb为转换获得的卫星本体坐标系坐标。并构建卫星的几何2 3 4 5
成像模型；应当注意的是：对于畸变模型：g(x)=g0+g1x+g2x+g3x+g4x+g5x，其中，该畸变模型是纯数学角度拟合的，各参数可分别解释为：0次项主要与传感器平移等误差有关，1次项主要与传感器安放倾斜等因素有关，2、4次项主要与镜头畸变等因素有关，3、5次项主要与传感器自身形变以及折光等因素有关。
[0106] 在本实施例中，
[0107] （a）根据Landsat7数据中提供的放射变换参数，将匹配以及手工修正过的图像控制点转入参考数据附带的地理参考系统中，即转为UTM分带的平面直角坐标系中；
[0108] （b）将UTM坐标系下的点位由投影变换的逆变换转入大地经纬度坐标系（本实例采用WGS84椭球的地心坐标系），并从DEM数据中以双线性内插方式读取该点位的高程值；
[0109] （c）认为内插高程值为点位大地高，将大地经纬度坐标转入地固地心三维直角坐标系；
[0110] （d）根据控制点所在的图像扫描线行号位置计算成像时间，内插出瞬时时刻的卫星外定向元素，结合卫星所在的位置地球自转线速度得到卫星本体坐标系到地固地心三维直角坐标系的转换矩阵，将地面点进一步转入卫星本体坐标系。转换矩阵为：其中：
[0111]
[0112] 上式中所指含时间变量的参数如：
[0113] Vt=(VXt,VYt,VZt)，Pt=(Xt,Yt,Zt)，Vt'=(VXt',VYt',VZt')，
[0114]
[0115]
[0116] 其中常量ω为地球自转的角速度：
[0117] 步骤5，根据步骤1中已经获取的参数初值和步骤4中的卫星的几何成像模型进行模型参数改化迭代解算；
[0118] 步骤6，输出经纠正后的模型结果及精度指标：计算结果输出图像点位残差分布图如图5所示，其与通用多项式模型解算的图像点残差中误差量化比对如下表所示。
[0119]
[0120] 本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
