CN103942775B - 基于最大核密度估计的相位相关亚像素匹配方法 - Google Patents

Abstract

本发明涉及一种基于最大核密度估计的相位相关亚像素匹配方法，具体为：①影像间互功率谱Q(u，v)的计算；②分离Q(u，v)实部与虚部，应用相位边缘滤波去噪，求取相位角矩阵ψ(u，v)；③对相位角矩阵解缠，得到对应的平面矩阵ψ_unwrap(u，v)；④使用MKDE方法对平面矩阵进行拟合，得到平面方程参数x₀和y₀，并根据相位角与空域像素的比例关系，计算获得亚像素匹配结果。与现有技术相比，本发明获得的匹配结果具有良好的稳健性，能够有效避免相位角矩阵中粗差对相位角平面拟合的影响，提高匹配精度，同时显著减小了亚像素匹配结果的像素锁现象。

Description

基于最大核密度估计的相位相关亚像素匹配方法

技术领域

本发明涉及遥感影像领域，尤其是涉及一种基于最大核密度估计的相位相关亚像素匹配方法。

背景技术

遥感影像的亚像素精确匹配是摄影测量与遥感领域的关键技术之一。通过影像间的亚像素匹配，能够精确的获取立体视差、地物位移、地表形变场等数据，在DEM生成、影像融合、地表形变监测等领域有着极其重要的应用。

基于相位相关的亚像素匹配方法，具有抗噪声性能好、受影像灰度影响小、匹配精度高等特点，近年来得到了广泛的关注。其中，Hoge等使用基于奇异值分解(singular value decomposition，SVD)的方法，将相位矩阵分解为两个一维向量，通过直线拟合的方式对亚像素位移值进行了计算；Nagashima等使用二维sinc函数对互功率谱峰值进行插值拟合(peak evaluation formula，PEF)，获得了亚像素精度；Liu等采用最小二乘和快速最大能量密度估计(2D-LS)的方法，对相位角平面进行拟合，通过相位角平面斜率计算了亚像素位移。

然而，现有的相位相关亚像素匹配方法依然存在三个问题：①通过相位相关获得的相位矩阵通常含有较大的粗差，传统的最小二乘拟合方法难以对相位角精确求解；②影像的过采样、相位矩阵的奇异值分解等方法，均有计算复杂度高，抗粗差性差等缺点；③亚像素匹配结果会产生像素锁现象(pixel locking effect)，直接影响匹配精度。上述的问题直接影响到基于相位相关的亚像素匹配精度和匹配结果的稳健性，并对运算量提出了较高的要求。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于最大核密度估计的相位相关亚像素匹配方法，对相位解缠后的相位角平面进行精确的拟合，通过拟合平面方程，得到了影像间的亚像素匹配结果，匹配结果具有良好的稳健性。

本发明的目的可以通过以下技术方案来实现：

一种基于最大核密度估计的相位相关亚像素匹配方法，包括以下步骤：

1)获取参考影像m(i，j)与待匹配影像s(i，j)，计算两个影像间在频域中的互功率谱矩阵Q(u，v)：

$Q(u,v)=\frac{M(u,v)S^{*}(u,v)}{\left|M(u,v)S^{*}(u,v)\right|}=\exp \{-i({\mathrm{ux}}_0+{\mathrm{vy}}_0)\}$

其中，M(u，v)、S(u，v)分别为m(i，j)、s(i，j)的傅里叶变换，且满足M(u，v)＝S(u，v)exp{-i(ux₀+vy₀)}，S(u，v)与S^*(u，v)互共轭，x₀、y₀为参考影像m(i，j)与待匹配影像s(i，j)在空域中行列方向存在的位移参数；

2)分离Q(u，v)的实部和虚部，获得实部矩阵Q_real(u，v)和虚部矩阵Q_imag(u，v)，并分别采用相位边缘滤波方法对实部和虚部进行滤波去噪，得到的Q_{real-filtered}(u，v)和Q_{imag-filtered}(u，v)：

$\left\{\begin{array}{c}{Q}_{\mathrm{real}-\mathrm{filtered}}(u,v)=\frac{1}{s^2}\underset{i,j\∈s}{\mathrm{\Σ}}{Q}_{\mathrm{real}}(i,j)\\ Q_{\mathrm{imag}-\mathrm{filtered}}(u,v)=\frac{1}{s^2}\underset{i,j\∈s}{\mathrm{\Σ}}{Q}_{\mathrm{imag}}(i,j)\end{array}\right.$

式中，s为滤波模板直径；

3)根据Q_{real-filtered}(u，v)和Q_{imag-filtered}(u，v)按反正切得到缠绕的相位角矩阵，对相位角矩阵进行解缠获取相位角平面ψ_unwrap(u，v)：

${\mathrm{\ψ}}_{\mathrm{unwrap}}(u,v)=\mathrm{unwrap}\{{\tan }^{-1}\left(\frac{Q_{\mathrm{real}-\mathrm{filtered}}(u,v)}{Q_{\mathrm{imag}-\mathrm{filtered}}(u,v)}\right)\};$

4)采用最大核密度估计方法对相位角平面进行拟合，估计参数x₀、y₀，并根据相位角与空域像素的比例关系，计算获得亚像素匹配结果。

所述的步骤3)中，对相位角矩阵进行解缠采用的方法为基于最小网络费用流的二维相位解缠算法。

所述的步骤4)中，采用最大核密度估计方法对相位角平面进行参数估计的具体流程为：

401)从候选数据点中随机选取一个数据子集；

402)由数据子集计算相位角平面的初始估计参数；

403)用相位角平面的初始估计参数计算全部候选点的残差，并计算残差的核密度

${\hat{f}}^{*}=\frac{1}{\mathrm{nh}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}K\left(\frac{r_i}{h}\right)$

式中，h为带宽，K为核函数，{r_i)_{i＝1，...，n}为候选点的残差，n为残差个数；

404)对η个随机数据子集重复步骤402)、403)后执行步骤405)；

405)选取核密度最大的数据子集对应的平面拟合参数作为最优候选参数，子集对应的候选点作为最优候选点；

406)对获取的最优候选参数和最优候选点，应用两步法尺度估计获取最优平面拟合参数。

所述的核函数的具体表达式为：

$K_e\left(u\right)=\left\{\begin{array}{ccc}\frac{3}{4}(1-u^2)& \mathrm{if}& u^2<1\\ 0& & \mathrm{otherwise}\end{array}\right..$

所述的步骤4)中，相位角与空域像素的比例系数为N/2π，N为离散傅里叶变换参数。

与现有技术相比，本发明通过相位边缘滤波方法对相位角矩阵进行滤波去噪，对得到的相位角矩阵进行二维相位解缠，使用基于最大核密度估计(maximumkernel density estimator，MKDE)的无参数模型参数估计方法对相位角平面进行精确拟合，获取相位角平面参数，得到影像间的亚像素匹配结果。与快速最大能量密度估计(quick maximum density power estimator，QMPDE)等基于无参数核密度估计的方法相比，MKDE在避免了子集的均值漂移(mean shift)遍历操作，在减小了运算量和计算复杂度的同时，达到了与之相当的精度。相对于LMedS、M估计和RANSAC等稳健估计方法，MKDE具有更强的抗粗差特性和稳健性。

实验结果表明，本发明获得的匹配结果具有良好的稳健性，能够有效避免相位角矩阵中粗差对相位角平面拟合的影响，同时显著减小了匹配结果的像素锁现象。在使用实验影像时，其影像间的匹配精度可以优于1/20像素。另外，本发明在具有较高匹配精度和抗粗差特性的同时，对Pixel-Locking现象的影响也有较好的抑制作用。

附图说明

图1为本发明实施例中的实验影像数据；

其中，(1a)为ASTER影像，(1b)为TH-1影像；

图2为相位角矩阵的滤波和解缠过程示意图；

其中，(2a)为使用相位边缘滤波前的相位角矩阵，(2b)为相位边缘滤波后的相位角矩阵，(2c)为相位缠绕的相位角矩阵的三维视图，(2d)为二位解缠后的相位角矩阵的三维视图；

图3为相位角平面的拟合过程示意图；

其中，(3a)为ASTER多光谱实验影像，(3b)为平面拟合前的相位角矩阵，(3c)为使用最小二乘方法拟合的相位角平面结果，(3d)为使用MKDE拟合的相位角平面结果；

图4为ASTER实验影像匹配结果；

图5为TH-1实验影像匹配结果；

图6为TH-1实验影像匹配结果的Pixel-Locking现象。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

一种基于最大核密度估计的相位相关亚像素匹配方法，是基于二维平面拟合的相位相关亚像素匹配的发展。其理论基础是二维傅立叶变换的平移性质，即参考影像与待匹配影像在空域内的平移，可以被表达为频域中二维傅立叶变换式的线性相位差。影像间相位相关得到的相位角矩阵，在解缠后，对应相位角空间中的一个二维平面。通过对相位角平面的精确拟合，可以得到了相位角平面在两个正交方向的倾角，进而对亚像素匹配结果求解。本方法的具体流程为：①影像间互功率谱Q(u，v)的计算；②分离Q(u，v)实部与虚部，应用相位边缘滤波去噪，求取相位角矩阵ψ(u，v)；③对相位角矩阵解缠，得到对应的平面矩阵ψ_unwrap(u，v)；④使用MKDE方法对平面矩阵进行拟合，得到平面方程参数x₀和y₀，并根据相位角与空域像素的比例关系(乘以比例系数N/2π，N为离散傅里叶变换参数)，计算获得亚像素匹配结果。

1、基于平面拟合的相位相关原理

获取参考影像m(i，j)与待匹配影像s(i，j)，参考影像m(i，j)与待匹配影像s(i，j)在空域中行列方向存在的水平位移和垂直位移x₀、y₀，即：

m(i，j)＝s(i+x₀，j+y₀) (1)

那么其在频域中其对应关系及互功率谱矩阵Q(u，v)可以表达为：

M(u，v)＝S(u，v)exp{-i(ux₀+vy₀)} (2)

$Q(u,v)=\frac{M(u,v)S^{*}(u,v)}{\left|M(u,v)S^{*}(u,v)\right|}=\exp \{-i({\mathrm{ux}}_0+{\mathrm{vy}}_0)\}---\left(3\right)$

其中，M(u，v)、S(u，v)分别为m(i，j)、s(i，j)的傅里叶变换，S(u，v)与S^*(u，v)互共轭。

将式(3)变换到极坐标下可得：

ψ(u，v)＝∠Q(u，v)＝ux₀+vy₀ (4)

根据文献“Estimating sub-pixel shifts directly from the phase difference[C]//ImageProcessing，2005.ICIP 2005.IEEE International Conference on.Genova：IEEE，2005.I-1057-60”的研究，式(4)对应一个平面方程，其中为ψ(u，v)矩阵Q(u，v)实部和虚部的夹角矩阵，即相位角矩阵。根据ψ(u，v)与u和v的关系，可以对式(4)的平面进行拟合，进而求解方程的参数x₀和y₀，从而在避免插值和傅里叶逆变换的情况下，得到亚像素位移值。

2、相位角矩阵的二维解缠

相位角矩阵ψ(u，v)自身为二维相位缠绕矩阵，在x和y方向有2π相位的缠绕。在进行平面拟合之前，需要对相位角矩阵进行解缠。但是受到互功率谱的噪声影响，直接对Q(u，v)进行解缠会导致相位角平面含有较多噪声，引起平面拟合误差。在进行相位解缠前，对互功率谱Q(u，v)的实部和虚部进行分离，使用相位边缘滤波方法(phase fringe filtering)对Q_real(u，v)和Q_real(u，v)分别滤波。对滤波后的得到的Q_{real-filtered}(u，v)和Q_{imag-filtered}(u，v)矩阵按反正切得到缠绕的相位角矩阵，最后对相位角矩阵进行解缠获取相位角平面ψ_unwrap(u，v)。

$\left\{\begin{array}{c}{Q}_{\mathrm{real}-\mathrm{filtered}}(u,v)=\frac{1}{s^2}\underset{i,j\&Element;s}{\mathrm{\&Sigma;}}{Q}_{\mathrm{real}}(i,j)\\ Q_{\mathrm{imag}-\mathrm{filtered}}(u,v)=\frac{1}{s^2}\underset{i,j\&Element;s}{\mathrm{\&Sigma;}}{Q}_{\mathrm{imag}}(i,j)\end{array}\right.---\left(5\right)$

${\mathrm{\&psi;}}_{\mathrm{unwrap}}(u,v)=\mathrm{unwrap}\left\{{\tan }^{-1}\left(\frac{Q_{\mathrm{real}-\mathrm{filtered}}(u,v)}{Q_{\mathrm{imag}-\mathrm{filtered}}(u,v)}\right)\right\}---\left(6\right)$

其中，s为滤波模板直径。使用的二维相位解缠算法为文献“A Novel PhaseUnwrapping Method Based on Network Programming[J].IEEE Transactions onGeoscience and Remote Sensing，1998，36(3)：813-821”中提出的基于最小网络费用流(minimum cost network flow，MCNF)的相位解缠算法。

3、最大核密度估计

最大核密度估计(MKDE)是一种基于无参数核密度估计技术的稳健估计方法。MKDE方法假设内点在全部候选点中仅占据相对多数，并且呈现高斯分布，当模型完美拟合时，在残差空间中0残差的点的概率密度应该尽可能高。因此，MKDE可以表达为：

$\widehat{\mathrm{\&theta;}}=\underset{\mathrm{\&theta;}}{\arg \max }{\hat{f}}_{\mathrm{\&theta;}}^{*}---\left(7\right)$

为一组n个候选点的残差{r_i}_{i＝1，…，n}的核密度，其表达式为：

${\hat{f}}^{*}=\frac{1}{\mathrm{nh}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}K\left(\frac{r_i}{h}\right)---\left(8\right)$

其中，h为带宽，K为核函数。针对平面模型的拟合，其候选点的残差空间为一维，在此情况下选用文献“MDPE：A very robust estimator for model fitting andrange image segmentation[J].International Journal of Computer Vision，2004，59(2)：139-166”中的Epanechnikov kernel作为核函数的表达为：

$K_e\left(u\right)=\left\{\begin{array}{ccc}\frac{3}{4}(1-u^2)& \mathrm{if}& u^2<1\\ 0& & \mathrm{otherwise}\end{array}\right.---\left(9\right)$

在带宽参数为h，随机数据子集选取个数为η的情况下，使用MKDE进行平面模型参数估计的流程如下：

(1)从候选数据点中随机选取一个数据子集；

(2)由数据子集计算平面模型的初始估计参数；

(3)用平面模型的初始估计参数计算全部候选点的残差，并计算残差的核密度

(4)对η个随机数据子集重复(2)(3)步操作；

(5)选取核密度最大的子集对应的平面拟合参数作为最优候选参数，子集对应的候选点作为最优候选点；

(6)对获取的最优候选参数和最优候选点，应用文献“Maximum kernel densityestimator for robust fitting[C]//Acoustics，Speech and Signal Processing，2008.ICASSP2008.IEEE International Conference on.Las Vegas：IEEE，2008.3385-3388”中的两步法尺度估计(two-step scale estimator，TSSE)获取最优平面拟合参数。

以下通过具体实验说明本发明。本实施例中所使用的实验影像数据采用文献“A fast direct Fourier-based algorithm for subpixel registration of images[J].IEEETransactions on Geoscience and Remote Sensing，2001，39(10)：2235-2243”中的方法生成已知真实亚像素视差的影像对。遥感影像数据来源于ASTER卫星和TH-1卫星的多光谱和全色影像，如图1所示，实验所用影像参数如表1所示。

表1

影像参数	ASTER	TH-1
			分辨率	15m	2m
尺寸	256×256Pixel	128×128Pixel
			波段	3N	PAN

实验分析

图2给出了是使用基于MKDE的相位相关匹配方法进行的影像匹配实验的计算过程及其结果，待匹配影像相对于参考影像在X方向和Y方向的真实视差值分别为6.25像素和4.75像素。

表2显示了基于MKDE、SVD、PEF和2D-LS的四种相位相关方法对图2对应的影像进行匹配的结果对比。由表2的结果可知，基于相位角平面拟合的MKDE和2D-LS方法整体优于使用PEF和SVD匹配方法，其中，MKDE方法精度优于2D-LS方法，在X方向和Y方向的匹配误差均小于1/20像素。

表2

真实视差值	PEF	SVD	2D-LS	MKDE
					X＝6.25	6.1639	6.1636	6.2074	6.2263
Y＝4.75	4.8043	4.8193	4.7474	4.7896

(1)相位角平面拟合

图(3a)为ASTER多光谱实验影像，待匹配影像相对于参考影像在X方向和Y方向的真实视差值分别为1.3像素和2.7像素。受多光谱不同波段影像低相关性的影响，(3b)中的相位角矩阵包含大量粗差。由(3c)和(3d)的结果可知，使用最小二乘方法在平面拟合时都会受到大量粗差的影响，造成拟合平面与真实平面存在误差倾角，直接影响匹配精度。而MKDE具有良好的抗粗差性，平面拟合过程中有效的避免了相位角矩阵粗差的影响，有效的提高了匹配精度。

表3给出了上述两种方法对平面拟合RMS及对应的匹配结果精度。由表3可知，MKDE方法在低相关影像其相位角矩阵含有大量粗差的情况下，平面拟合的RMS小于最小二乘(LS)的结果，其匹配精度显著优于LS匹配方法，其在X方向和Y方向的匹配误差小于1/15像素。

表3

平面拟合方法	LS	MKDE
			平面拟合RMS	1.7030	1.1726
X＝1.3	1.0773	1.3184
			Y＝2.7	2.4214	2.6367

(2)匹配精度

图4、图5和图6分别给出了使用ASTER和TH-1两组模拟实验影像，应用MKDE、SVD、PEF和2D-LS四种相位相关方法的匹配精度对比。ASTER实验影像组包含10个影像对，其视差真实值以0.25像素递增，TH-1实验影像组各包含10个影像对，其视差真实值以0.1像素递增。

由图4、图5可知，基于MKDE的相位相关匹配精度优于基于PEF、SVD和2D-LS的相位相关匹配方法，在两组实验中匹配误差均优于1/20像素。值得注意的是，实验中使用基于PEF、SVD和LS的相位相关匹配的误差，出现了明显的接近正弦函数的周期性趋势。相对而言，基于MKDE的相位相关匹配误差的周期性趋势远小于上述方法。其匹配误差的周期性趋势体现为：在真实视差值为整像素和半像素时，影像匹配误差达到最小，而在真实视差值为0.2～0.4像素和0.6～0.8像素时，匹配误差会显著增大，从而使匹配误差随真实视差值的递增而产生周期性的规律。这种在亚像素匹配过程中，获得的亚像素位移值趋向于整像素和半像素而引起的***性误差现象被称为Pixel-Locking现象。

(3)Pixel Locking现象的抑制

为了对Pixel Locking现象引起的误差进行分析，本文使用基于MKDE、SVD的PEF的三种相位相关匹配方法，对模拟实验影像进行匹配实验。使用TH-1实验影像，包含50个影像对，在X方向其视差真实值以0.1像素递增。图6给出了上述三种方法匹配误差随真实视差值变化的对比结果。

图6表明了，在使用TH-1实验影像时，使用PEF和SVD方法的匹配误差随真实视差值变化出现了明显的大幅度的周期性变化趋势，而使用MKDE方法的匹配误差亦随真实视差值变化出现了周期性变化趋势，但受影响幅度远小于PEF和SVD方法，匹配误差总体趋于稳定。

在亚像素匹配过程中，当真实视差值非整像素或半像素时，Pixel-Locking现象会对匹配精度产生显著影响，其匹配误差的变化趋势具有***周期性。基于MKDE的相位相关匹配方法的匹配误差明显减小，但只在一定程度上抑制了Pixel-Locking现象的影响，说明除了匹配方法以外，影像自身的特性也对Pixel-Locking现象产生影响。实验结果表明基于MKDE的相位相关匹配方法在具有较高匹配精度和抗粗差特性的同时，对Pixel-Locking现象的影响也有较好的抑制作用。

Claims (4)

1.一种基于最大核密度估计的相位相关亚像素匹配方法，其特征在于，包括以下步骤：

1)获取参考影像m(i,j)与待匹配影像s(i,j)，计算两个影像在频域中的互功率谱矩阵Q(u,v)：

$Q(u,v)=\frac{M(u,v)S^{*}(u,v)}{\left|M(u,v)S^{*}(u,v)\right|}=\exp \{-i({\mathrm{ux}}_0+{\mathrm{vy}}_0)\}$

其中，M(u,v)、S(u,v)分别为m(i,j)、s(i,j)的傅里叶变换，且满足M(u,v)＝S(u,v)exp{-i(ux₀+vy₀)}，S(u,v)与S^*(u,v)互共轭，x₀、y₀为参考影像m(i,j)与待匹配影像s(i,j)在空域中行列方向存在的位移参数；

2)分离Q(u,v)的实部和虚部，获得实部矩阵Q_real(u,v)和虚部矩阵Q_imag(u,v)，并分别采用相位边缘滤波方法对实部和虚部进行滤波去噪，得到的Q_{real-filtered}(u,v)和Q_{imag-filtered}(u,v)：

$\left\{\begin{array}{c}{Q}_{real-filtered}\left(u,v\right)=\frac{1}{s^2}\underset{i,j\&Element;s}{\mathrm{\&Sigma;}}{Q}_{real}\left(i,j\right)\\ Q_{imag-filtered}\left(u,v\right)=\frac{1}{s^2}\underset{i,j\&Element;s}{\mathrm{\&Sigma;}}{Q}_{imag}\left(i,j\right)\end{array}\right.$

式中，s为滤波模板直径；

3)根据Q_{real-filtered}(u,v)和Q_{imag-filtered}(u,v)按反正切得到缠绕的相位角矩阵，对相位角矩阵进行解缠获取相位角平面ψ_unwrap(u,v)：

${\mathrm{\&psi;}}_{unwrap}\left(u,v\right)=unwrap\left\{{\tan }^{-1}\left(\frac{Q_{real-filtered}\left(u,v\right)}{Q_{imag-filtered}\left(u,v\right)}\right)\right\};$

4)采用最大核密度估计方法对相位角平面进行拟合，估计参数x₀、y₀，并根据相位角与空域像素的比例关系，计算获得亚像素匹配结果；

所述的步骤4)中，采用最大核密度估计方法对相位角平面进行参数估计的具体流程为：

401)从候选数据点中随机选取一个数据子集；

402)由数据子集计算相位角平面的初始估计参数；

403)用相位角平面的初始估计参数计算全部候选点的残差，并计算残差的核密度

${\hat{f}}^{*}=\frac{1}{nh}\underset{i=1}{\overset{n}{\&Sigma;}}K\left(\frac{r_i}{h}\right)$

式中，h为带宽，K为核函数，{r_i}_i＝1,...,n为候选点的残差，n为残差个数；

404)对η个随机数据子集重复步骤402)、403)后执行步骤405)；

405)选取核密度最大的数据子集对应的平面拟合参数作为最优候选参数，子集对应的候选点作为最优候选点；

406)对获取的最优候选参数和最优候选点，应用两步法尺度估计获取最优平面拟合参数。

2.根据权利要求1所述的一种基于最大核密度估计的相位相关亚像素匹配方法，其特征在于，所述的步骤3)中，对相位角矩阵进行解缠采用的方法为基于最小网络费用流的二维相位解缠算法。

3.根据权利要求1所述的一种基于最大核密度估计的相位相关亚像素匹配方法，其特征在于，所述的核函数的具体表达式为：

$K_e\left(u\right)=\left\{\begin{array}{ccc}\frac{3}{4}(1-u^2)& if& u^2<1\\ 0& & otherwise\end{array}\right..$

4.根据权利要求1所述的一种基于最大核密度估计的相位相关亚像素匹配方法，其特征在于，所述的步骤4)中，相位角与空域像素的比例系数为N/2π，N为离散傅里叶变换参数。