CN109191501B - 一种多源高空间分辨率影像的自动几何配准方法 - Google Patents

Abstract

一种多源高空间分辨率影像的自动几何配准方法，步骤(1)、使用参考影像对待配准影像进行重采样处理；步骤(2)、对处理后的待配准影像划分网格；步骤(3)、计算每个网格内的待配准影像和零矩阵影像之间的互信息量，互信息量最小值所对应的区域作为待配准影像中的特征匹配区域；步骤(4)、使用特征匹配区域遍历参考影像，计算每个网格内的特征匹配区域和参考影像之间的归一化互信息量，将归一化互信息量最大值所对应的区域作为参考影像匹配区域；步骤(5)、提取特征匹配区域和参考影像匹配区域的中心点坐标，计算出最终的转换参数；步骤(6)待配准影像使用最终的转换参数，得到影像配准后的结果。同名点分布均匀、配准精度高、速度快。

Description

一种多源高空间分辨率影像的自动几何配准方法

技术领域

本发明涉及测绘地理信息技术领域，特别涉及一种多源高空间分辨率影像的自动几何配准方法。

背景技术

随着遥感技术的发展，不同时间、空间和光谱分辨率的遥感数据被广泛应用于各行各业，为我国的城市规划、城市建设、城市管理和生态文明等提供了重要的数据源。不同数据的几何配准是实现多源遥感影像综合应用的前提和关键步骤。当前，遥感影像配准方法主要包括人工配准和自动配准两种。

人工配准是在专业软件中导入待配准影像和参考影像，基于人工识别方式寻找影像之间的同名点，利用四参数模型计算待配准影像坐标信息来实现影像的配准，该方法工作量大、费时费力，且配准的精度会因人工选点的准确度而产生一定影响。

自动配准方法是利用算法程序，自动识别待配准影像和参考影像的特征区域来确定两景影像的同名点，并通过同名点计算四参数模型来实现影像的配准。目前常用的自动配准方法包括基于特征和基于灰度的自动配准方法，基于特征的自动配准方法是通过提取边缘、纹理和区域等特征信息进行图像之间的匹配，常见的有Moravec、Forstner、Harris、SIFT等算法，该类方法适用于匹配区域较小的图像，当匹配区域较大时，则这类方法的匹配效率明显降低，且对于特征不明显的区域，误匹配会明显增多，配准精度受到影响；基于灰度的自动配准方法主要依赖图像上的灰度统计信息，现有的方法会因不同传感器影像的灰度特征差异而导致配准精度不高，且基于全局的灰度统计也会让大区域影像的匹配效率不高。因此，对于高空间分辨率影像，匹配区域通常较大，需要一种既能保证配准精度，又能快速完成的自动几何配准方法。

发明内容

本发明所要解决的技术问题是针对现有技术的现状，提供一种配准精度高、同名点分布均匀和配准速度快的多源高空间分辨率影像的自动几何配准方法。

本发明解决上述技术问题所采用的技术方案为：一种多源高空间分辨率影像的自动几何配准方法，其特征在于：包括以下步骤：

步骤(1)、输入待配准影像和参考影像，使用参考影像对待配准影像进行重采样处理，得到N×M待配准影像；

步骤(2)、对处理后得到的N×M待配准影像划分网格，划分成为k×k的网格，其中，单个网格矩阵为INT(N/k)×INT(M/k)，网格数为k²；其中，k为正整数，k≥2，INT()函数为取整函数；

步骤(3)、使用j×j的零矩阵影像遍历步骤(2)中的每个网格内对应的待配准影像，分别计算每个网格内j×j区域的待配准影像和零矩阵影像之间的互信息量，将每个网格内j×j区域的互信息量最小值所对应的区域作为该网格内的待配准影像中的特征匹配区域；其中j为正整数，j×j区域远小于单个网格的INT(N/k)×INT(M/k)区域；

步骤(4)、分别使用每个网格内待配准影像的特征匹配区域遍历参考影像，并分别计算每个网格内待配准影像的特征匹配区域和参考影像之间的归一化互信息量，将归一化互信息量的最大值所对应的区域作为该网格内待配准影像的特征匹配区域相对应的参考影像匹配区域；

步骤(5)、提取每个网格内的待配准影像中的特征匹配区域的中心点坐标和相对应的参考影像匹配区域的中心点坐标，形成k²对坐标点对，计算出满足残差要求的转换参数，作为最终的转换参数；

步骤(6)、在步骤(1)中输入的待配准影像中使用步骤(5)中最终的转换参数，得到影像配准后的结果。

具体的，所述步骤(1)的重采样处理的计算公式为：

其中，(N_原始,M_原始)为输入的待配准影像的分辨率，(N′,M′)为参考影像的分辨率，(N,M)为重采样处理后的待配准影像的分辨率，INT()函数为取整函数。

作为改进，所述步骤(3)中的每个网格内j×j区域的待配准影像和零矩阵影像之间的互信息量的计算公式为：

其中，U为待配准影像的矩阵，[0]为零矩阵影像的矩阵，I(U；[0])为待配准影像和零矩阵影像之间的互信息值，H(U)为待配准影像的信息熵，H([0])为零矩阵影像的信息熵，H(U,[0])为待配准影像和零矩阵影像之间的联合信息熵，P(u)为待配准影像的概率密度，P([0])为零矩阵影像的概率密度，P(u,[0])为待配准影像和零矩阵影像之间的联合概率密度；

设

其中，u_x,y为待配准影像的矩阵U中的像素点，x为像素点u_x,y对应在待配准影像的矩阵U中的横坐标，y为像素点u_x,y对应在待配准影像的矩阵U中的纵坐标，h(u_x,y)为待配准影像的像素直方图统计值；h(u_x,y,[0])为待配准影像的像素点和零矩阵影像之间的像素联合直方图统计值；z为待配准影像像素的分类数目；其中z为预设值；

根据公式(2)，将公式(1)的互信息量的计算公式转换为：

进一步的，所述步骤(4)中计算每个网格内待配准影像的特征匹配区域和参考影像之间的互信息量的计算公式如下：

其中，W为特征匹配区域的矩阵，V为参考影像的矩阵，I(W；V)为特征匹配区域和参考影像之间的互信息值，H(W)为特征匹配区域的信息熵，H(V)为参考影像的信息熵，H(W,V)为特征匹配区域和参考影像之间的联合信息熵，P(w)是特征匹配区域的概率密度，P(v)为参考影像的概率密度，P(w,v)为特征匹配区域和参考影像之间的联合概率密度；

设

其中，

为特征匹配区域的矩阵W中的像素点，x₀为像素点

对应在特征匹配区域的矩阵W中的横坐标，y₀为像素点

对应在特征匹配区域的矩阵W中的纵坐标，

为参考影像的矩阵V中的像素点，x₁为像素点

对应在参考影像矩阵V中的横坐标，y₁为像素点

对应在参考影像矩阵V中的纵坐标，

为特征匹配区域影像的像素直方图统计值，

为参考影像的像素直方图统计值，

为特征匹配区域影像的像素和参考影像之间的像素联合直方图统计值；z₁为特征匹配区域影像和参考影像像素的分类数目；其中z₁为预设值；

因此根据公式(4)和公式(5)，待配准影像的特征匹配区域与参考影像之间的互信息量计算公式为：

对公式(6)归一化处理，得到待配准影像的特征匹配区域和参考影像之间的归一化互信息量计算公式为：

其中，

为特征匹配区域的像素点

的最大值，

为特征匹配区域的像素点

的最小值，

为参考影像的像素点

的最大值，

为参考影像的像素点

的最小值。

具体的，所述步骤(5)中的具体步骤为：

步骤(5-1)、采用平面四参数转换模型对k²对坐标点对进行计算，计算出转换参数，并计算出每对坐标点对的转换残差；

步骤(5-2)、判断是否有坐标点对的转换残差大于阈值，如有，则剔除该坐标点对，得到剩余的坐标点对，转入步骤(5-3)，反之，则输出该转换参数；

步骤(5-3)、使用步骤(5-1)中的方法，使用剩余的坐标点对计算转换参数并计算出每对坐标点对的转换残差，直至每对坐标点对的转换残差均小于阈值时，则输出满足残差要求时所对应的转换参数，作为最终的转换参数。

与现有技术相比，本发明的优点在于：通过将待配准影像划分成k×k的网格，使待配准影像和参考影像之间的匹配同名点分布均匀，并且通过计算各网格内的待配准影像与零矩阵影像之间的互信息量，提取出待配准影像的特征匹配区域，有效降低了计算量，提高了计算速度，并通过计算特征匹配区域和参考影像之间的归一化互信息量，有效提高特征匹配区域和参考影像的配准精度，因此配准效率高。

附图说明

图1为本发明实施例的配准方法的流程图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

如图1所示，一种多源高空间分辨率影像的自动几何配准方法包括以下步骤：

步骤(1)、输入待配准影像和参考影像，使用参考影像对待配准影像进行重采样处理，得到N×M待配准影像；本实施例中，N＝18210，M＝12760；

其中，重采样处理的计算公式为：

其中，(N_原始,M_原始)为输入的待配准影像的分辨率，(N′,M′)为参考影像的分辨率，(N,M)为重采样处理后的待配准影像的分辨率，INT()函数为取整函数；

步骤(2)、对处理后得到的N×M待配准影像划分网格，划分成为k×k的网格，其中，单个网格矩阵为INT(N/k)×INT(M/k)，网格数为k²；其中，k为正整数，k≥2，INT()函数为取整函数；

对N×M的待配准影像分别从上至下、从左至右进行划分，对划分后待配准影像右侧和下侧多余不能划分的部分舍弃，本实施例中，待配准影像的矩阵为18210×12760，k＝6，因此单个网格矩阵为3035×2126，网格数为36，待配准影像的12757～12760列的影像不进行划分；该步骤将待配准影像划分成网格，主要目的是使待配准影像和参考影像之间的匹配坐标点分布均匀。

步骤(3)、使用j×j的零矩阵影像遍历步骤(2)中的每个网格内对应的待配准影像，分别计算每个网格内j×j区域的待配准影像和零矩阵影像之间的互信息量，将每个网格内j×j区域的互信息量最小值所对应的区域作为该网格内的待配准影像中的特征匹配区域；其中j为正整数，j×j区域远小于单个网格的INT(N/k)×INT(M/k)区域；本实施例中，j＝16；

其中，每个网格内j×j区域的待配准影像和零矩阵影像之间的互信息量的计算公式为：

其中，U为待配准影像的矩阵，[0]为零矩阵影像的矩阵，I(U；[0])为待配准影像和零矩阵影像之间的互信息值，H(U)为待配准影像的信息熵，H([0])为零矩阵影像的信息熵，H(U,[0])为待配准影像和零矩阵影像之间的联合信息熵，P(u)为待配准影像的概率密度，P([0])为零矩阵影像的概率密度，P(u,[0])为待配准影像和零矩阵影像之间的联合概率密度；

设

其中，u_x,y为待配准影像的矩阵U中的像素点，x为像素点u_x,y对应在待配准影像的矩阵U中的横坐标，y为像素点u_x,y对应在待配准影像的矩阵U中的纵坐标，h(u_x,y)为待配准影像的像素直方图统计值；h(u_x,y,[0])为待配准影像的像素点和零矩阵影像之间的像素联合直方图统计值；z为待配准影像像素的分类数目；其中z为预设值，本实施例中z＝8；

根据公式(2)，将公式(1)的互信息量的计算公式转换为：

该步骤获取了待配准影像各网格内待匹配区域的信息量大小，若区域的互信息值较大，说明该区域与零矩阵影像较为相似，信息量较少，反之则信息量较多。该步骤的主要目的是为了提取待配准影像的关键信息，用少量关键信息作为配准的依据，能有效降低计算量，提高算法的计算速度。

步骤(4)、分别使用每个网格内待配准影像的特征匹配区域遍历参考影像，并分别计算每个网格内待配准影像的特征匹配区域和参考影像之间的归一化互信息量，将归一化互信息量的最大值所对应的区域作为该网格内待配准影像的特征匹配区域相对应的参考影像匹配区域；

计算每个网格内待配准影像的特征匹配区域和参考影像之间的互信息量的计算公式如下：

其中，W为特征匹配区域的矩阵，V为参考影像的矩阵，I(W；V)为特征匹配区域和参考影像之间的互信息值，H(W)为特征匹配区域的信息熵，H(V)为参考影像的信息熵，H(W,V)为特征匹配区域和参考影像之间的联合信息熵，P(w)是特征匹配区域的概率密度，P(v)为参考影像的概率密度，P(w,v)为特征匹配区域和参考影像之间的联合概率密度；

设

其中，

为特征匹配区域的矩阵W中的像素点，x₀为像素点

对应在特征匹配区域的矩阵W中的横坐标，y₀为像素点

对应在特征匹配区域的矩阵W中的纵坐标，

为参考影像的矩阵V中的像素点，x₁为像素点

对应在参考影像矩阵V中的横坐标，y₁为像素点

对应在参考影像矩阵V中的纵坐标，

为特征匹配区域影像的像素直方图统计值，

为参考影像的像素直方图统计值，

为特征匹配区域影像的像素和参考影像之间的像素联合直方图统计值；z₁为特征匹配区域影像和参考影像像素的分类数目；其中z₁为预设值，本实施例中，z₁＝8；

因此根据公式(4)和公式(5)，待配准影像的特征匹配区域与参考影像之间的互信息量计算公式为：

对公式(6)归一化处理，得到待配准影像的特征匹配区域和参考影像之间的归一化互信息量计算公式为：

其中，

为特征匹配区域的像素点

的最大值，

为特征匹配区域的像素点

的最小值，

为参考影像的像素点

的最大值，

为参考影像的像素点

的最小值。

步骤(5)、提取每个网格内的待配准影像中的特征匹配区域的中心点坐标和相对应的参考影像匹配区域的中心点坐标，形成k²对坐标点对，计算出满足残差要求的转换参数，作为最终的转换参数；本实施例中，k²＝36；

具体步骤为：

步骤(5-1)、采用平面四参数转换模型对k²对坐标点对进行计算，计算出转换参数，并计算出每对坐标点对的转换残差；

步骤(5-2)、判断是否有坐标点对的转换残差大于阈值，如有，则剔除该坐标点对，得到剩余的坐标点对，转入步骤(5-3)，反之，则输出该转换参数；本实施例中，阈值为2个像元；

步骤(5-3)、使用步骤(5-1)中的方法，使用剩余的坐标点对计算转换参数并计算出每对坐标点对的转换残差，直至每对坐标点对的转换残差均小于阈值时，则输出满足残差要求时所对应的转换参数，作为最终的转换参数。

步骤(6)、在步骤(1)中输入的待配准影像中使用步骤(5)中最终的转换参数，得到影像配准后的结果。

该方法适用于多源高空间分辨率影像的配准，如：GF数据，ZY数据及无人机影像，通过将待配准影像划分网格，使待配准影像和参考影像之间的匹配同名点分布均匀，提高了配准的稳定性和精度；并且通过计算每个网格内j×j区域的待配准影像和零矩阵影像之间的互信息量，并根据互信息量提取出待配准影像中信息量较大的特征匹配区域，解决了现有配准算法中的直接采用整景的待配准影像与参考影像进行匹配，计算量大的问题，有效降低了计算量，提高了算法的效率；且通过计算特征匹配区域和参考影像之间的归一化互信息量，寻找同名区域，减少了数据自身引起的不确定性及误差，有效提高了配准精度。

Claims (5)

1.一种多源高空间分辨率影像的自动几何配准方法，其特征在于：包括以下步骤：

步骤(1)、输入待配准影像和参考影像，使用参考影像对待配准影像进行重采样处理，得到N×M待配准影像；

步骤(2)、对处理后得到的N×M待配准影像划分网格，划分成为k×k的网格，其中，单个网格矩阵为INT(N/k)×INT(M/k)，网格数为k²；其中，k为正整数，k32，INT()函数为取整函数；

步骤(3)、使用j×j的零矩阵影像遍历步骤(2)中的每个网格内对应的待配准影像，分别计算每个网格内j×j区域的待配准影像和零矩阵影像之间的互信息量，将每个网格内j×j区域的互信息量最小值所对应的区域作为该网格内的待配准影像中的特征匹配区域；其中j为正整数，j×j区域远小于单个网格的INT(N/k)×INT(M/k)区域；

步骤(4)、分别使用每个网格内待配准影像的特征匹配区域遍历参考影像，并分别计算每个网格内待配准影像的特征匹配区域和参考影像之间的归一化互信息量，将归一化互信息量的最大值所对应的区域作为该网格内待配准影像的特征匹配区域相对应的参考影像匹配区域；

步骤(5)、提取每个网格内的待配准影像中的特征匹配区域的中心点坐标和相对应的参考影像匹配区域的中心点坐标，形成k²对坐标点对，计算出满足残差要求的转换参数，作为最终的转换参数；

步骤(6)、在步骤(1)中输入的待配准影像中使用步骤(5)中最终的转换参数，得到影像配准后的结果。

2.根据权利要求1所述的多源高空间分辨率影像的自动几何配准方法，其特征在于：所述步骤(1)的重采样处理的计算公式为：

其中，(N_原始,M_原始)为输入的待配准影像的分辨率，(N′,M′)为参考影像的分辨率，(N,M)为重采样处理后的待配准影像的分辨率，INT()函数为取整函数。

3.根据权利要求1所述的多源高空间分辨率影像的自动几何配准方法，其特征在于：所述步骤(3)中的每个网格内j×j区域的待配准影像和零矩阵影像之间的互信息量的计算公式为：

其中，U为待配准影像的矩阵，[0]为零矩阵影像的矩阵，I(U；[0])为待配准影像和零矩阵影像之间的互信息值，H(U)为待配准影像的信息熵，H([0])为零矩阵影像的信息熵，H(U,[0])为待配准影像和零矩阵影像之间的联合信息熵，P(u)为待配准影像的概率密度，P([0])为零矩阵影像的概率密度，P(u,[0])为待配准影像和零矩阵影像之间的联合概率密度；

设

P([0])＝1，

其中，u_x,y为待配准影像的矩阵U中的像素点，x为像素点u_x,y对应在待配准影像的矩阵U中的横坐标，y为像素点u_x,y对应在待配准影像的矩阵U中的纵坐标，h(u_x,y)为待配准影像的像素直方图统计值；h(u_x,y,[0])为待配准影像的像素点和零矩阵影像之间的像素联合直方图统计值；z为待配准影像像素的分类数目；其中z为预设值；

根据公式(2)，将公式(1)的互信息量的计算公式转换为：

。

4.根据权利要求1所述的多源高空间分辨率影像的自动几何配准方法，其特征在于：所述步骤(4)中计算每个网格内待配准影像的特征匹配区域和参考影像之间的互信息量的计算公式如下：

其中，W为特征匹配区域的矩阵，V为参考影像的矩阵，I(W；V)为特征匹配区域和参考影像之间的互信息值，H(W)为特征匹配区域的信息熵，H(V)为参考影像的信 息熵，H(W,V)为特征匹配区域和参考影像之间的联合信息熵，P(w)是特征匹配区域的概率密度，P(v)为参考影像的概率密度，P(w,v)为特征匹配区域和参考影像之间的联合概率密度；

设

其中，

为特征匹配区域的矩阵W中的像素点，x₀为像素点

对应在特征匹配区域的矩阵W中的横坐标，y₀为像素点

对应在特征匹配区域的矩阵W中的纵坐标，

为参考影像的矩阵V中的像素点，x₁为像素点

对应在参考影像矩阵V中的横坐标，y₁为像素点

对应在参考影像矩阵V中的纵坐标，

为特征匹配区域影像的像素直方图统计值，

为参考影像的像素直方图统计值，

为特征匹配区域影像的像素和参考影像之间的像素联合直方图统计值；z₁为特征匹配区域影像和参考影像像素的分类数目；其中z₁为预设值；

因此根据公式(4)和公式(5)，待配准影像的特征匹配区域与参考影像之间的互信息量计算公式为：

对公式(6)归一化处理，得到待配准影像的特征匹配区域和参考影像之间的归一化互信息量计算公式为：

其中，

为特征匹配区域的像素点

的最大值，

为特征匹配区域的 像素点

的最小值，

为参考影像的像素点

的最大值，

为参考影像的像素点

的最小值。

5.根据权利要求1所述的多源高空间分辨率影像的自动几何配准方法，其特征在于：所述步骤(5)中的具体步骤为：

步骤(5-1)、采用平面四参数转换模型对k²对坐标点对进行计算，计算出转换参数，并计算出每对坐标点对的转换残差；

步骤(5-2)、判断是否有坐标点对的转换残差大于阈值，如有，则剔除该坐标点对，得到剩余的坐标点对，转入步骤(5-3)，反之，则输出该转换参数；

步骤(5-3)、使用步骤(5-1)中的方法，使用剩余的坐标点对计算转换参数并计算出每对坐标点对的转换残差，直至每对坐标点对的转换残差均小于阈值时，则输出满足残差要求时所对应的转换参数，作为最终的转换参数。

Images