一种GF‑4卫星序列图像自动相对配准方法与流程

本发明涉及遥感图像几何处理技术，具体的说，涉及一种针对高分辨率静止卫星序列图像的自动相对配准技术。

背景技术：

遥感图像的几何预处理一直是遥感数据处理的主要课题之一。新的遥感卫星投入使用后，对新数据的预处理便是紧要的问题。遥感图像的预处理通常包含几何预处理与辐射预处理，其中的几何预处理由于每个卫星的成像设计与轨道的不同，需要针对具体的卫星图像开发专门的几何预处理流程。几何预处理通常包含系统几何校正与几何精纠正两步主要处理。系统几何校正通过卫星成像的严格几何模型，建立卫星图像像素点与对应的地面点之间的几何转换关系，并将卫星图像映射到某一地图投影下。经过系统几何校正的卫星图像具有统一的地图投影，可依据地图投影进行多幅图像的镶嵌，但同一区域多幅图像按照地图投影叠加后，由于系统几何校正的精度有限，通常无法达到逐像素的严格配准，并且卫星成像过程还会发生设备老化、器械抖动及姿态偏差等造成的各种非系统误差。要实现图像之间的逐像素配准，通常的方法是在系统几何校正的结果上进行几何精纠正。几何精纠正将匹配获得的控制点数据用于拟合各种数学纠正模型，再利用得到的数学模型纠正图像，达到修正图像内部的局部几何变形的目的，常用的数学纠正模型有多项式模型、样条函数模型及Delaunay三角网模型等。

高分四号卫星(以下简称GF-4)是中国于2015年12月发射的一颗地球同步轨道卫星，搭载空间分辨率为50米的全色、多光谱相机和400米分辨率的中波红外相机，采用面阵凝视方式成像，成像间隔快至20秒，具备高时间、高空间分辨率的优势。自2016年2月3号国防科技工业局公布首批图像以来，GF-4已经获取了中国及周边区域大量数据，在检测森林火灾、洪涝灾害等方面发挥着重要作用。GF-4卫星图像的预处理同样包括几何与辐射两个部分，几何预处理包括系统成像模型的构建、控制点匹配与几何精纠正等，目标是实现同一地图投影下成像数据的逐像素配准。辐射处理包括辐射定标、云及云下阴影的检测等，目标是使得图像像素值能精确地描述地表辐射状况。

GF-4卫星数据预处理技术研发在充分利用之前卫星数据处理成果的基础上，还要考虑GF-4卫星图像本身的特性，研发专门的处理算法。面临的关键问题主要有两个：一是高分辨率静止卫星图像内部几何变形的校正，该问题对于公里级低分辨率静止卫星图像并不突出，但对于50米的高分辨率静止卫星图像而言，对能否实现逐像素配准关系密切，因为更高分辨率下，图像内部的几何变形情况更加复杂；二是如何利用基于序列图像的处理技术提高GF-4图像预处理结果精度，GF-4卫星的面阵凝视成像方式，成像间隔快至20秒，短时间内可以获取同一地理区域下大量的图像，考虑到云的运动特性，利用序列图像对于提高云及云下阴影的检测精度非常必要。

对于上述两个关键问题，如果能在不经过系统几何校正的基础上实现图像之间的快速配准，配准后的数据将会成为GF-4图像预处理过程中的关键中间数据产品，尤其是对于工程化生产云及云下阴影的检测数据产品而言，需要一个针对GF-4卫星图像的序列图像的自动、快速配准算法。

技术实现要素：

本发明的目的是针对GF-4静止卫星图像预处理应用，提供一种自动的序列图像自动配准技术，特别是对于GF-4在凝视模式下对同一区域进行观测获取的序列图像的云及云下阴影检测数据产品的生产，提供序列图像逐像素配准数据准备。本技术基于成熟的SIFT的自动匹配算法，根据GF-4卫星图像的预处理需求与成像几何特性定制的自动、快速配准处理流程。

本发明的基本思路为：GF-4静止卫星对同一地理区域获取的序列图像数据，在没有系统几何模型的情况下，仅利用卫星记录的图像中心点及四个角点的近似经纬度坐标确定图像的大致相对位置关系，然后根据图像定位误差范围确定图像分块自动匹配检索范围，从分块的图像区域中获得至少3组不同区域的控制点，利用这些控制点拟合出一个线性函数。对于同一地理区域，不同时间获取的n(n＞2)幅图像相邻两幅图像之间，都能得到一个线性函数，利用这n-1个线性函数，实现n幅图像之间逐像素的配准。需要注意的是，这里的序列图像配准并不对图像本身进行变换与重采样，配准的结果只有n-1组线性函数系数，在后续的GF-4图像云与云下阴影检测等预处理中利用该组系数实现n幅图像的逐像素配准。

所述的n幅GF-4卫星序列图像，限定为50米空间分辨率的全色、多光谱图像，图像四个角点的经纬度差异不超过±0.3度，n幅图像可以是GF-4卫星凝视模式下获取的序列图像，也可以是不同天获取的，按照获取时间先后排列的序列图像。

本发明的技术方案提供的GF-4卫星序列图像自动配准方法，其特征在于包括以下实施步骤：

A数据预处理，生成不精确的云与阴影掩膜；

B将序列图像中的所有图像都划分成4行4列，共16个网格区域，结合每个图像的掩膜数据与四个角点的近似地理坐标，对两幅图像的每个网格区域的数据进行基于二级网格划分的快速自动匹配，获取一组控制点数据；

C对两幅图像的每个网格区域匹配获取的控制点数据中选取一个代表控制点；

D如果序列中两幅图像的16个网格区域中获得了2个或者更多的代表控制点，利用这些代表控制点拟合出一个线性函数，线性函数的参数作为结果输出。

上述实施步骤的特征在于：

步骤A中所述数据预处理，指算法程序对输入的序列图像进行数据完整性检查、地理覆盖范围检查、按获取时间排序、以及程序运行的一些准备初始化处理。生成不精确的云与阴影掩膜通过标记过亮与过暗像素完成，掩膜排除的区域不参与后续的自动匹配，可提高算法的速度与鲁棒性。

步骤B中所述的自动匹配，是在序列图像中的两幅图像之间，采用二级网格划分的策略，对于一组包含n幅图像的序列图像，第一级网格划分先将每幅图像划分为4行4列，共16个网格块，自动匹配过程在(n-1)×16个块中独立并行运算，每个一级网格块再进行二级网格划分，划分为20行20列，共400个网格，在两幅图像之间，遍历其中一个图像的二级网格，根据图像四个角点近似坐标，对应到另一幅图像对应的16个二级网格区域进行自动匹配，当有3个二级网格都成功匹配得到控制点后，停止该块二级网格遍历并输出所有控制点数据。通过二级网格划分后，自动匹配仅在一个128×128像素区域与512×512像素区域展开，运算量大幅减少，利用基于SIFT特征的自动匹配算法，结合RANSAC算法剔除误匹配点，可快速获得控制点数据。

步骤C中所述代表控制点，是从通过对两幅图像的每个网格区域匹配获取的控制点数据对进行误差分析与数据拟合，利用线性最小二乘法，选取所有控制点中误差最小的点。

步骤D中，如果获得了2个控制点，则严格求解出线性函数的两个参数，如果获得了多余2个的控制点，则采用线性最小二乘拟合获得线性函数的两个参数。若n个序列图像全部处理成功，则算法最终输出n-1组线性函数参数。将线性函数的参数作为最终输出结果，序列图像通过线性函数实现逐像素的相对配准，处理过程中不对图像本身进行插值与重采样。

本发明与现有技术相比有如下特点：本发明提供了一种简单的GF-4序列图像自动配准解决方案，直接通过线性函数实现GF-4序列图像逐像素位置关系的配准，算法结构适合大数据量的快速并行处理，具有处理流程简单、无需人工交互、运算速度快的特点，而且在算法的稳定性以及适用性上都具有优势。对于GF-4卫星数据预处理中的序列图像云及云下阴影检测，以及GF-4卫星数据系统几何模型在轨修正，提供了关键的数据预处理支撑。

附图说明：

图1是GF-4卫星序列图像自动匹配流程图

图2是二级网格划分示意图

图3是二级网格单个块的自动匹配示意图

具体实施方式：

本技术的思想是使用线性函数描述GF-4卫星序列图像位置偏移关系，达到逐像素配准，其必要条件为：GF-4卫星在不同时间获取的同一区域的多幅图像之间仅存在旋转与平移的位置关系，其内部复杂的几何畸变近似一致。该必要条件与GF-4卫星成像特性有关，具体分析如下：

首先，GF-4卫星采用静止卫星轨道。GF-4卫星是中国第一颗地球同步轨道遥感卫星，运行在距地36000公里的地球同步轨道，定点位置105.6°E。静止卫星相对于地球的位置是固定的，其成像几何是不变的，对地球可观测范围内的任意一点到卫星传感器成像点的几何关系都是固定的。静止卫星位置的固定性，保证了对同一区域的多个观测图像，在图像中心点与四个角点坐标相同的情况下，所有图像的系统成像几何模型是相同的，图像的系统畸变与空间分辨率都是一致的。

其次，GF-4卫星采用面阵凝视方式成像。目前卫星搭载的高分辨率多光谱传感器多采用线阵CCD传感器，线阵CCD传感器推扫或摆扫式的成像方式，单幅成像需要花费一定的时间，期间卫星姿态的抖动与轨道的偏移都会引入几何畸变，导致不同时间获取的同一区域数据内部几何畸变不一致。不同于线阵CCD传感器，面阵方式成像时间是瞬时的，图像的成像过程不会引入新的几何畸变，卫星的抖动与传感器的抖动都很难对成像造成影响。

最后，GF-4卫星成像定位精度高。通过指向控制，GF-4可实现对中国及周边地区进行自由观测，也可采用凝视模式对幅宽400千米的固定区域进行持续观测，其定位精度达±0.1度，并且GF-4卫星的观测控制一般限定了成像中心点确定的情况下，固定成像的四个角点，使得GF-4卫星具有获得同一固定区域序列图像的能力。

GF-4卫星实际获取的数据具有很好的序列特性。自2016年2月3号国防科技工业局公布首批图像以来，GF-4卫星已经获取了中国及周边区域大量数据，其中包含大量的可构成序列的数据，通过数据抽样验证，在GF-4卫星凝视模式下获取的同一天相近时刻的序列图像，通过线性函数关系得到的配准精度最大误差不超过2个像素，对于不同天获取的序列图像，通过线性函数关系得到的配准精度平均误差也在2个像素以内，个别区域误差超过2个像素，估计是由大气差异导致的光学几何畸变，这些误差对于后续的云及云检测应用、以及系统几何模型修正研究的影响可以忽略。

采用本发明实现GF-4卫星序列图像自动匹配的实施流程如图1所示，现结合附图对其进行描述。

处理单元111数据预处理与掩膜处理，数据预处理针对GF-4卫星序列图像，对于GF-4数据发布的50米空间分辨率的可见光与近红外图像，有两种方式可以获得序列的图像：一是GF-4相机在凝视工作模式下获取的同一天、相近时间内获取的连续多幅图像组成的序列；二是非凝视模式下，不同天获取的、同一地区的多幅图像。在GF-4卫星的实际运营中，获取第二种非凝视模式构成的序列图像数据很容易。因为灾害监控作为GF-4的重要任务，当某一区域发生灾害时，GF-4可以在经纬度误差0.1度的范围内对固定区域进行多天、多次观测，从而构成序列图像。数据预处理算法程序对输入的序列图像进行数据完整性检查、地理覆盖范围检查、按获取时间排序、以及程序运行的一些准备初始化处理。在检查数据完整性的过程中统计各图像各波段的直方图，根据蓝光波段的直方图统计，利用自动阈值法标记过亮与过暗像素，作为不精确的云与阴影掩膜，掩膜数据按照8位Byte型保存为外部临时文件。

处理单元112自动匹配，该处理是算法流程的关键技术，该处理过程通过两个关键步骤实现：一是SIFT(Scale-invariant feature transform)特征提取；二是特征匹配。

第一步：SIFT特征提取。

SIFT算法是在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量。SIFT特征是指在图像尺度空间中明显不同于周围和相临尺度上同区域点的局部极值，是差分高斯函数在图像局部的最大/最小值(Lowe，2004)。SIFT特征点在灰度值上表现为比其周围点和相临尺度上同区域点的灰度值都大或者都小的点。目前David Lowe的SIFT特征提取算法程序(Lowe，2005)较为成熟。

第二步：特征匹配。

特征匹配是从两个图像提取的SIFT特征点中寻找相同特征点的过程。对于每一个特征点，根据图像间的位置关系生成每一个候选匹配点的粗匹配点，形成粗匹配点对。然后以每一对粗匹配点为中心，在匹配的两幅图像上分别提取模板窗口(M×M)和搜索窗口(N×N，N＞M)，这里，M取值128，N取值512，让模板窗口在搜索窗口内移动，计算每一个位置与搜索窗口的相似度，如果相似度最大值大于某一阈值T，这里T取值0.65，则将相应匹配位置作为该对控制点的精确匹配结果，并记录相似度的最大值。

相似性判别准则采用常用的归一化相关系数(NCC)：

由于直接利用NCC进行模板匹配计算量较大，采用NCC的加速算法Fast NCC(J.P Lewis，1995)大幅度减少匹配时间。

特征匹配完成后，经常会存在一些误匹配点，剔除误匹配点采用RANSAC(Random Sample Consensus)算法。该方法是从一组包含异常数据的样本数据集中，采用迭代的方法进行拟合模型，多次迭代后获得最佳模型拟合参数，再根据容许误差将所有的匹配点对分为内点和外点，外点就是需要剔除的误匹配点。

考虑到GF-4卫星L1级数据图像单幅像素尺寸达10240×10240，若整幅图像进行SIFT特征提取与匹配将会消耗大量的运算时间与资源，所以采用二级网格处理的策略，大幅减小运算量的同时便于程序的并行处理。二级网格划分示意图见图2。假定一组序列图像有n幅，每幅图像对应有云与阴影粗掩膜，第一级网格划分先将每幅图像划分为4行4列，共16个网格块，则自动匹配过程可划分为(n-1)×16个独立的块来并行运算。一级网格划分后，单个网格的尺寸为2560×2560像素，由于对于单个一级网格自动匹配的目的是获得一个代表控制点，而不是获得尽量多的控制点，所以没必要对整个网格进行图对图的SIFT自动匹配，这样运算量大且存在大量不必要运算。这里对单个网格再次进行二级网格划分，划分为20行20列，共400个网格。在两幅图像之间，遍历其中一个图像的二级网格，对于不含有掩膜标记区域的网格，根据图像四个角点近似坐标，对应到另一幅图像对应的16个二级网格区域，此时自动匹配仅在一个128×128像素区域与512×512像素区域展开，当有3个二级网格都成功匹配得到控制点后，停止二级网格遍历并输出所有控制点数据，从而大幅减少了运算量。二级网格单个块的自动匹配示意图见图3。其中地理坐标粗对应是通过卫星获取图像时的指向信息所记录的图像四个角点的坐标计算得到的大致经纬度坐标，根据GF-4卫星实际数据四个角点的坐标精度情况，确定在一个图像中的二级网格划分中的小块可以在另一个图像的二级网格划分中的4行4列，共16个小块，的范围内找到匹配控制点。二级网格单个块的自动匹配算法伪代码如下：

处理单元113计算代表控制点，通过对两幅图像的每个网格区域匹配获取的控制点数据进行误差分析与数据拟合，利用线性最小二乘法，选取所有控制点中误差最小的一个控制点作为代表控制点。

处理单元114拟合线性配准参数。两幅图像之间的线性配准参数为旋转与平移，通过一阶多项式描述，即y＝ax+b。如果获得了2个控制点，则严格求解出线性函数的两个参数a与b，如果获得了多于2个的控制点，则采用线性最小二乘拟合获得线性函数的两个参数。若n个序列图像处理成功，则算法输出n-1组线性函数参数。

本发明实现的GF-4卫星序列图像自动配准处理流程比传统的卫星数据实现逐像素配准的处理流程要简化。传统的方法通过系统几何校正与几何精纠正两步处理，在将原始卫星图像处理为包含投影坐标的地图的同时实现逐像素的配准，该过程包含多次图像插值与重采样处理，并且为了校正图像局部区域的几何畸变，有时需要采用人机交互的方式添加部分控制点用于几何精纠正，整个处理过程复杂且耗时。本发明的处理流程简化，仅通过图像局部小区域的少量配准控制点来拟合线性变换参数，实现序列图像之间相对的逐像素配准，处理过程不包含图像插值与重采样处理，处理过程简单、速度快，且算法自动，无需人机交互。同时本发明的处理流程充分考虑网格化并行处理，算法流程结构适用于CPU与GPU加速编程。

本发明在GF-4卫星数据预处理中有两个专门的应用：第一个应用在GF-4云检测数据产品生产流程中，作为序列图像云及云下阴影检测处理前必须的数据预处理步骤。GF-4云检测数据产品使用未经系统几何校正的L1级数据作为输入，序列图像云检测主要利用云的运动变化特性区分云与地表，利用本发明得到的序列图像线性变换参数，在不对图像进行插值与重采样的情况下，实现序列图像相对的逐像素配准。第二个应用是为GF-4卫星数据系统几何模型在轨修正提供数据预处理支撑。

本发明的实例已经在PC平台上实现，目前已经交付用户方进行测试与使用。GF-4凝视卫星获取的数据量巨大，本发明充分利用算法结构的并行特性，通过CPU与GPU加速技术，处理速度接近数据IO速度，满足了工程化数据处理的需求。

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。