基于云控制的国产陆地观测卫星姿轨精化方法与流程

本发明涉及一种基于广义“云控制”信息的国产陆地观测卫星姿轨精化方法，其中基于参考底图的全自动匹配控制点和基于姿轨的卫星影像条带平差是该方法的关键技术。

背景技术：

随着高分系列(如高分一号、高分二号、高分六号等)、资源系列(资源一号、资源三号)等国产卫星的不断发射，对地观测的卫星影像数据快速积累，为相关科学研究和实际工程应用提供了大量卫星影像。大比例尺测图需求对于卫星影像的几何定位精度提出越来越高的要求，因此，利用广义多源的地面控制信息作为几何参考，通过几何检校和几何定向，从而提高卫星影像的直接对地定位精度，是国产卫星影像地面处理的重要任务。

几何检校是提高卫星影像直接对地定位精度的有效手段。传统的几何检校主要依赖专门设计的地面检校场。美国、法国等卫星影像处理技术较为先进的国家，已在世界范围内建立了地面检校场，我国也建立了嵩山检校场、中卫检校场等。这些地面检校场通常是专门设计的，检校场范围内一般具备多种地形(如平地、山地)且具有一定的高差，而且均匀的布设了大量高精度的地面控制点。利用卫星过境时拍摄的影像便可以进行相应的几何检校。然而，地面检校场建设成本大、维护费用高，且几何检校易受天气影像，具体而言，当天气条件较差时(如雨雾天影像不清晰、有云覆盖时不能拍摄有效影像)难以获取满足检校要求的影像。利用控制点的几何定向也是比较常用的方法，但是控制点通常获取周期长，并不是适用于常态的卫星影像几何处理。

技术实现要素：

本发明主要解决现有的卫星影像地面处理系统中几何定位精度不够高的问题。

基于广义云控制的几何检校和几何定向可有效解决以上问题。云控制包括数字正射影像(digitalorthophotomap,dom)和数字表面模型(digitalsurfacemodel,dsm)或数字高程模型(digitalelevationmodel)，通过卫星影像与它们的匹配便可获取大量的控制点，在此基础上进行区域网平差便可解算卫星影像的精确几何定向参数，从而实现高精度的几何定位。

本方法基于已有地理信息数据，实现了gf1、gf1b/c/d、gf2、gf6、zy3-02等不同分辨率的国产光学卫星1a级标准景全色影像与参考底图之间的控制点自动量测，获取控制点数据，通过标准景的分景信息对1a级影像进行逻辑拼接，恢复原始条带影像信息，结合条带影像的行时、gps轨道、姿态等数据，进行基于原始姿轨的条带平差，解算精确的姿态参数。同时，分析像点残差的分布规律，并根据残差规律对卫星的ccd线阵进行内方位几何检校。最后根据检校后的内方位元素和解算的精确轨道、姿态(即外方位元素)，重新拟合标准景1a影像的有理多项式参数(rationalpolynomialcoefficients,rpc)。本发明的算法优化精度高，适用范围广，计算速度快，具有较好的通用性。

本发明提出了一种基于已有dom和dem/dsm构成的广义云控制信息，实现国产陆地观测卫星，包括高分一号(gf-1)、高分一号b/c/d(gf-1b/c/d)、高分二号(gf-2)、高分六号(gf-6)和资源三号02(zy3-02)等卫星姿轨精化的方法，以提高卫星影像直接对地定位精度。所提方法为全自动，处理效率高，支持并行处理，适合各种传感器类型的卫星影像。

本发明的上述技术问题主要是通过下述技术方案得以解决：

本发明的方法以已有参考数字正射影像和数字表面模型或数字高程模型作为几何参考，通过全自动匹配控制点和基于姿轨数据的卫星影像条带平差，对国产卫星影像进行姿轨精化，提高卫星影像的几何精度，可用于各种型号的卫星影像，其核心过程包括以下步骤：

步骤1，数据解压：读取整轨分景1a级影像内所有卫星影像文件，卫星影像原始文件格式为tar.gz，将其解压缩到临时目录下；

步骤2，分景控制点匹配：对步骤1中解压缩得到的文件目录，判断目录中是否存在全色影像和多光谱影像tiff文件，若不存在，则返回错误信息；若存在，读取该影像对应的属性信息，包括影像行列数、拍摄时间、轨道号和条带号，然后分别对全色影像和多光谱影像进行控制点匹配，生成该景影像的匹配结果文件，即像点与控制点数据，匹配的控制点数据与卫星传感器的观测数据构成区域网平差的输入数据，以下简称平差数据；

步骤3，逻辑条带拼接与条带像点构建：整轨分景1a影像中所有影像匹配完成后，根据分景信息中提供的每景影像的起始行与终止行，对整条轨道的逻辑条带拼接，并根据逻辑条带信息对步骤2获取的分景影像的像点进行统一坐标转换处理，恢复整轨的像点；

步骤4，条带平差解算精确定向参数：利用步骤2生成的整条带的平差数据，首先根据卫星的参数进行内方位元素初始化，并设置原始观测数据、自检校参数、定向片、漂移改正以及观测值权重，然后基于卫星影像姿轨参数进行线阵影像定向片法的条带平差，得到整条轨道影像的精确外方位元素和像点残差；

步骤5，基于步骤4计算得到的整轨影像的外方位元素，内插得到分景1a影像中的全色影像和多光谱影像每一行的精确外方位元素，基于地形无关算法重新拟合得到新的rpc系数，生成新的rpc文件。

进一步的，步骤2中控制点匹配采用初始地形约束的多层级金字塔影像匹配方法，该算法的具体实现方式如下，

(2.1)基于初始地形约束的影像几何变形改正：①在参考正射影像上确定匹配的基准窗口，根据正射影像具有的地理信息，得到基准窗口4个角点的平面坐标，然后在全球srtm上内插其高程z；②利用影像的rpc参数将基准窗口的4个角点投影至影像上，构成匹配窗口γ′；③窗口纠正，利用仿射变换建立匹配窗口γ′与基准窗口γ之间的变换关系，并将γ′纠正至γ，从而改正几何变形；

(2.2)多层级金字塔影像匹配：以一定的缩放系数对卫星影像和参考正射影像分别建立金字塔影像；首先在参考正射影像上提取harris特征点；以每个特征点为中心确定基准窗口，按上述基于初始地形约束的影像几何变形改正方法，将基准窗口投影至卫星影像并进行几何改正；利用互相关系数作为相似性测度匹配同名点，匹配从影像金字塔的顶层开始，在每一级金字塔影像匹配后，利用匹配的控制点对该影像进行基于rpc参数的单片定向，根据定向的像点残差剔除粗差点度；然后将解算的定向参数约束下一级影像的匹配。

进一步的，步骤4中基于卫星影像姿轨参数进行的线阵影像定向片法条带平差的具体实现方式如下，

(4.1)按一定行间隔选定定向片：沿着卫星飞行方向，即影像从上到下方向，按一定的时间间隔δt选定若干影像行，这些选定的影像行即定向片orientationimage，它们的外方位元素作为未知数在平差过程中求解，其它成像时刻的外方位元素通过内插算法计算得到；

(4.2)利用像点、控制点、初始内外方位元素，根据线阵影像共线条件方程列立误差方程式，并进而构建法方程；

(4.3)解算法方程得到未知数的解，并统计像点残差。

进一步的，步骤(4.1)中所述内插算法包括拉格朗日内插、分段多项式内插。

进一步的，步骤5中基于地形无关算法重新拟合得到新的rpc系数的具体实现方式如下，

(5.1)影像高程范围统计：计算影像四个角点的经纬度坐标，统计全球srtm在该范围内最大高程与最小高程，作为该影像范围的高程范围，等间隔设定若干个高程值；

(5.2)虚拟三维格网的建立：在影像的像方以若干个像素为间距，进行二维格网划分，根据每个格网点的像点坐标和不同高程分层的高程值，利用平差得到的定向参数计算它们的物方坐标，并转换为经纬度坐标，得到虚拟的三维空间格网；

(5.3)rpc参数拟合：利用上述虚拟空间格网点，拟合rpc参数。

本发明具有如下优点：适用范围广，包括高分一号(gf-1)、高分一号b/c/d(gf-1b/c/d)、高分二号(gf-2)、高分六号(gf-6)和资源三号02(zy3-02)等卫星姿轨精化均适用。本发明方法可以有效提高卫星影像几何精度，大幅降低影像残余误差，另外本方法支持全自动运行，支持并行处理，可以保证整体运算的高效性和稳定性。

附图说明

图1为基于初始地形约束的影像几何变形改正原理示意图；

图2为影像金字塔示意图；

图3整轨影像逻辑拼接示意图；

图4为本发明的总体流程图；

图5为控制点自动匹配结果的示意图。图中，每个小三角形均为匹配点；

图6为整轨卫星影像，按景为单元进行并行匹配的示意图。不同的计算节点，同时并行匹配不同景的卫星影像；

图7为姿轨精化效果示意图。图中，左边为原始影像，接边差为180米，右边为精化纠正后的影像，接边差大幅减小为4米。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

本发明提供的技术方案是，一种以已有参考数字正射影像和数字表面模型或数字高程模型作为基础，通过全自动匹配控制点和基于姿轨数据的卫星影像条带平差，对国产卫星影像进行姿轨精化，提高影像的精度，可用于各种型号的卫星影像的方法。如图4所示，包括以下步骤：

步骤1，数据解压。读取整轨分景1a级影像内所有卫星影像文件，卫星影像原始文件格式为tar.gz，将其解压缩到临时目录下。

步骤2，分景控制点匹配。对步骤1中解压缩得到的文件目录，判断目录中是否存在全色影像和多光谱影像tiff文件，若不存在，则返回错误信息；若存在，读取该影像对应的属性信息，包括影像行列数、拍摄时间、轨道号和条带号等，然后分别对全色影像和多光谱影像进行控制点自动匹配。

控制点匹配采用初始地形约束的多层级金字塔影像匹配方法，该算法特点如下：

(2.1)基于初始地形约束的影像几何变形改正。卫星影像由于拍摄时成像角度和地形起伏的影响，地面物体在不同影像上的形状不同，即呈现一定的几何畸变，比如地面垂直的方角在影像呈现锐角或者钝角，这将大大降低影像匹配的成功率和精度。本发明按图1所示的方法进行影像几何畸变改正：①在参考正射影像上确定匹配的基准窗口(例如23×23像素，如图1中的γ)，根据正射影像具有地理信息可得到基准窗口4个角点的平面坐标(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4)然后在全球srtm上内插其高程z；②利用影像的rpc参数将基准窗口的4个角点投影至影像上，构成匹配窗口γ′；③窗口纠正，利用仿射变换建立匹配窗口γ′与基准窗口γ之间的变换关系，图1中p(x,y)是原始卫星影像匹配窗口中的点，p(x,y)是参考影像匹配窗口对应的点，图1公式中(x,y)和(x’,y’)分别是卫星影像上的点仿射变换前和仿射变换后的坐标，a0、b0是仿射变换中的平移分量，a1、b1、a2、b2是仿射变换中的旋转和缩放分量，利用匹配窗口γ′与基准窗口γ之间对应的4组角点，即(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4)和(x′1,y′1)、(x′2,y′2)、(x′3,y′3)、(x′4,y′4)，根据图1公式建立4个方程组，解算得到平移分量a0、b0以及旋转和缩放分量a1、b1、a2、b2，并将γ′纠正至γ，从而改正几何变形；

(2.2)多层级金字塔影像匹配。以2为缩放系数对卫星影像和参考正射影像分别建立金字塔影像，如附图2所示。首先在参考正射影像上提取harris特征点；以每个特征点为中心确定基准窗口，按上述基于初始地形约束的影像几何变形改正方法，将基准窗口投影至卫星影像并进行几何改正；利用互相关系数作为相似性测度匹配同名点。匹配从影像金字塔的顶层开始，在每一级金字塔影像匹配后，利用匹配的控制点对该影像进行基于rpc参数的单片定向，根据定向的像点残差剔除粗差点(像点残差超过4倍的像点中误差)，避免其影响定向的精确度；然后将解算的定向参数约束下一级影像的匹配。

匹配过程中，算法首先根据全色影像和多光谱影像的地理范围在参考底图数据库检索参考影像，参考底图数据库由数字正射影像和数字高程模型组成，通常而言，可使用公众地理信息数据，例如全球数字高程模型shuttleradartopographymission(srtm)或aw3d30(alosworld3d–30m)、landsatetm+数字正射影像等。利用高精度的参考底图数据，可以解决传统基于地面定标场的几何检校方法的以下问题：(1)地面定标场建设成本高；(2)几何检校依赖于天气条件，当有云遮挡或空气质量较差导致不能有效获取清晰高的卫星影像时，则不能进行几何检校；(3)检校周期长，不利于常态化检校。卫星的具有一定的重复周期，一般为数月，而基于地面定标场的几何检校需要卫星过境拍摄的影像数据，故检校周期也较长；(4)全局性差。卫星几何检校参数有可能与摄影的地理位置(如纬度)相关，地面检校场相对于全球范围而言仍为局部范围，故利用其得到的检校参数不能满足全球范围，从而导致精度下降。参考底图由已有的全球范围地理信息大数据组成，成本低，分布全球，可随时进行任意位置的几何检校且可得到符合全球范围的精确几何检校参数，故具有技术优势。

参考底图搜索完成后，对全色影像进行多层级金字塔影像匹配，然后对多光谱影像进行多层级金字塔影像匹配。匹配完成后，匹配效果如图5所示，每个三角形即为一个匹配点，然后生成该景影像的匹配结果文件，即像点与控制点数据(从参考数据中获取的地面坐标)。匹配是自动并行处理的，不同的计算节点分别对不同景的卫星影像进行自动匹配处理，如图6所示。匹配的控制点数据与卫星传感器的观测数据(包括由星载gps观测的轨道位置、陀螺仪和星敏感器感测的姿态信息、行时信息)构成区域网平差的输入数据，以下简称平差数据。

步骤3，逻辑条带拼接与条带像点构建。整轨分景1a影像中所有影像匹配完成后，根据分景信息中提供的每景影像的起始行与终止行，对整条轨道的逻辑条带拼接，如附图3所示，对于第i景影像中的点(x,y)，点的坐标都是相对于该景影像，但我们需要对整轨影像一起处理，那么需要将点(x,y)的坐标由相对于i景影像变换为整轨影像统一的坐标系下。图中第i景影像有起始行si的信息，那么该景影像中的点(x,y)坐标就转换为(x,y+si)，对所有点的坐标根据起始行和终止行信息统一做转换，便把所有景影像都变换到统一的坐标系下，即进行了条带的拼接。并根据逻辑条带信息对步骤2获取的分景影像的像点进行统一坐标转换处理，恢复整轨的像点。

步骤4，条带平差解算精确定向参数。利用步骤2与步骤3生成的整条带的平差数据，首先根据卫星的参数进行内方位元素(卫星相机的焦距、像主点偏移值，这些值由卫星运营单位提供)初始化，并设置原始观测数据、自检校参数、定向片、漂移改正以及观测值权重，然后基于卫星影像姿轨参数进行的线阵影像定向片法条带平差。具体过程如下：

(4.1)按一定行间隔选定定向片。沿着卫星飞行方向，即影像从上到下方向(对应不同的影像行)，按一定的时间间隔δt选定若干影像行(根据起始时间t0和时间间隔δt可计算选定影像的成像时间，由于线阵影像每一行都已知成像时间，故可对应为影像行)，这些选定的影像行即定向片(orientationimage)，它们的外方位元素作为未知数在平差过程中求解，其它成像时刻的外方位元素通过内插算法(如拉格朗日内插、分段多项式内插等)计算得到；

(4.2)利用像点、控制点、初始内外方位元素，根据线阵影像共线条件方程列立误差方程式，并进而构建法方程；

(4.3)解算法方程得到未知数的解，并统计像点残差。

平差完成后，得到整条轨道影像的精确外方位元素和像点残差。

步骤5，基于步骤4计算得到的整轨影像的外方位元素，内插得到分景1a影像中的全色影像(前视、后视和下视)和多光谱影像每一行的精确外方位元素，基于“地形无关算法”重新拟合得到新的rpc系数，生成新的rpc文件。具体拟合过程为：

(1)影像高程范围统计。计算影像四个角点的经纬度坐标，统计全球srtm在该范围内最大高程与最小高程，作为该影像范围的高程范围，等间隔设定若干个高程值(一般为9个)；

(2)虚拟三维格网的建立。在影像的像方以200个像素为间距，进行二维格网划分，根据每个格网点的像点坐标和不同高程分层的高程值，利用平差得到的定向参数计算它们的物方坐标，并转换为经纬度坐标，得到虚拟的三维空间格网；

(3)rpc参数拟合。利用上述虚拟空间格网点，拟合rpc参数。

如果需要处理的数据不是整轨影像，而是单景影像。同步骤1和2，首先进行数据的解压缩和判断全色影像以及多光谱影像是否存在。然后在参考底图文件夹中搜索参考影像，分别对全色影像和多光谱影像进行多层级金字塔影像匹配。匹配完成后，进行纠正和rpc系数拟合。

本发明实施例的精化结果如图7所示，左边为原始影像，接边差为180米，右边为精化纠正后的影像，接边差大幅减小为4米。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。