CN114066749A - 相位相关抗噪位移估计方法、设备及存储介质 - Google Patents

相位相关抗噪位移估计方法、设备及存储介质 Download PDF

Info

Publication number
CN114066749A
CN114066749A CN202111228184.3A CN202111228184A CN114066749A CN 114066749 A CN114066749 A CN 114066749A CN 202111228184 A CN202111228184 A CN 202111228184A CN 114066749 A CN114066749 A CN 114066749A
Authority
CN
China
Prior art keywords
phase
singular vector
displacement
image
interference
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202111228184.3A
Other languages
English (en)
Inventor
余志武
张云生
杨振
戴吾蛟
邢磊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
National Engineering Laboratory for High Speed Railway Construction Technology
Original Assignee
Central South University
National Engineering Laboratory for High Speed Railway Construction Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University, National Engineering Laboratory for High Speed Railway Construction Technology filed Critical Central South University
Priority to CN202111228184.3A priority Critical patent/CN114066749A/zh
Publication of CN114066749A publication Critical patent/CN114066749A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种相位相关抗噪位移估计方法、设备及存储介质，引入相位标准偏差，在进行真实干涉相位估计前，首先将信噪比非常低的高频区域自动地剔除，以免在进行干涉相位估计时对真实信号产生过多干扰，大大提高了抗噪性能；同时，在进行解缠后相位直线拟合时，引入高斯函数定权，给予干涉相位低频部分更高的权值，高频部分低的权值，既保证了观测值的数量，同时也保证了观测值的质量，在此基础上，使用残差进行二次定权拟合直线，进一步减少了异常值的干扰，从而大大提高了抗噪性能。相较于现有方法，本发明具有很强的抗噪性能，可以在成像条件很差的情况下准确地估计出位移量。

Description

相位相关抗噪位移估计方法、设备及存储介质
技术领域
本发明属于图像区域匹配技术领域，尤其涉及一种基于PSD+SVD分解及带权直线拟合的相位相关抗噪位移估计方法、设备及存储介质。
背景技术
1975年，Kuglin和Hines首次提出数字图像相位相关(phase-correlation)算法，并将其应用到影像配准中，其相较于数字图像灰度相关，具有一定的抗噪性能、对光线变化不敏感等优势，因此被很多学者用来计算图像间位移。最原始的相位相关方法只能得到整像素的位移，随着相位相关研究的深入，通过相位直线拟合、相位平面拟合、齿状信号条纹数目估计、频率估计、局部上采样等方式的使用，使得相位相关方法可以得到亚像素精度的位移。
虽然相位相关方法相对于灰度相关方法具有更好的抗噪性能，但在某些场景中，如在高频拍摄、环境光线不足等场景中，所获得的图像存在信噪比极低的问题。为了进一步提高相位相关的抗噪性能，有学者提出使用SVD(Singular Value Decomposition，奇异值分解)分解、RANSAC(Random Sample Consensus，随机抽样一致)算法以及最小二乘二维相位拟合方法(例如文献“A subspace identification extension to the phasecorrelation method[MRI application][J].Medical I maging，IEEE Transactions on，2003，22(2)：277-280)等进一步解决相位相关位移估计中的噪声问题，提高相位相关位移估计抗噪性能。虽然这些方法可以进一步提高相位相关的抗噪性能，但仍存在一定的问题，比如使用SVD分解恢复的干涉相位和原始不受噪声干扰的相位具有一定差异；使用RANSAC算法估计位移对于某些受噪声干扰特别严重的情况可能不适用。
总的来说，目前所有方法在成像条件较好的情况下都能达到较高的精度，但随着成像条件的逐渐变差，精度均会有一定的损失，有些方法甚至会出现估计不到位移情况，即使使用了抗噪方法，如SVD分解、RANSAC算法等，也不能保证完全消除噪声干扰。
发明内容
本发明的目的在于提供一种相位相关抗噪位移估计方法、设备及存储介质，以解决低质量成像时噪声干扰问题，本发明在经典相位相关基础上，引入相位标准偏差(Phasestandard deviation，简称PSD)对两张影像求得的互功率谱干涉条纹进行评价，剔除干涉效果很差的区域，以增加相关矩阵的置信度，减少影像噪声的干扰。
本发明是通过如下的技术方案来解决上述技术问题的：一种相位相关抗噪位移估计方法，包括以下步骤：
步骤1：读取发生位移前、后的影像，分别对两幅影像进行二维傅里叶变换得到对应的频谱图；对所述频谱图进行频谱中心化，使频谱图中低频部分分布于频谱图中心区域，高频部分分布于频谱图边缘区域；
步骤2：对经所述步骤1处理后的两幅频谱图分别进行加窗，再对加窗后的两幅频谱图进行相位干涉，得到干涉相位图；
步骤3：计算所述干涉相位图中每个元素对应的相位标准偏差，并由所有相位标准偏差生成相位标准偏差图；
步骤4：以所述相位标准偏差图的中心元素作为种子点，设置种子生长阈值，进行种子生长，得到所述干涉相位图中干涉条纹清晰区域；
步骤5：获取所述干涉条纹清晰区域的最小外接矩阵，对所述最小外接矩阵所包含的区域进行裁剪，得到相应区域块；
步骤6：采用SVD分解对所述区域块进行秩1估计，提取SVD分解后最大特征值对应的左奇异向量和右奇异向量，并分别计算所述左奇异向量和右奇异向量的相位；
步骤7：对所述左奇异向量和右奇异向量的相位分别进行相位解缠，对左奇异向量和右奇异向量相位解缠后的相位分别进行第一次定权，对第一次定权后的相位及其对应的位置序号进行直线拟合，得到第一次拟合直线；
根据所述第一次拟合直线模型的残差，对左奇异向量和右奇异向量相位解缠后的相位进行第二次定权，对第二次定权后的相位及其对应的位置序号进行直线拟合，得到第二次拟合直线；
步骤8：根据左奇异向量和右奇异向量对应的第二次拟合直线的斜率，再结合所述影像的长和宽，计算发生位移后影像相对于发生位移前影像在x方向和y方向的平移量。
进一步地，所述步骤2中，采用汉宁窗对每幅所述频谱图进行加窗，所述汉宁窗的计算公式为：
Figure BDA0003315070490000021
其中，Whn(i,j)表示频谱图第i行第j列位置处的汉宁窗函数值，(i,j)分别表示影像中某像素点的行列号，rows表示影像中像素点的总行数，cols表示影像中像素点的总列数。
进一步地，所述步骤2中，干涉相位图的计算公式为：
Figure BDA0003315070490000022
其中，Q(u,v)表示干涉相位图，即相关矩阵，F1(u,v)表示发生位移前影像经二维傅里叶变换后对应的频谱图，F2(u,v)表示发生位移后影像经二维傅里叶变换后对应的频谱图，
Figure BDA0003315070490000031
表示取F2(u,v)的共轭，x0表示发生位移后影像相对于发生位移前影像在x方向的位移，y0表示发生位移后影像相对于发生位移前影像在y方向的位移，u表示像素坐标系的横坐标，v表示像素坐标系的纵坐标。
进一步地，所述步骤3中，计算相位标准偏差的具体实现过程为：
步骤3.1：定义滑动窗口大小w*w，其中w表示边长且为奇数；
步骤3.2：对所述干涉相位图中每个元素求取相位，得到对应的相位矩阵M；
步骤3.3：构造与所述相位矩阵M维度一致的相位标准偏差矩阵P，对所述相位标准偏差矩阵P中的所有元素进行初始化；
步骤3.4：对于所述相位矩阵M中的每个元素，以该元素为中心，取w*w滑动窗口区域内的所有元素构成一个新矩阵M'，按照以下公式计算对应位置的相位标准偏差：
Figure BDA0003315070490000032
其中，
Figure BDA0003315070490000033
表示对应位置的相位标准偏差，对应位置是指滑动窗口的中心位置，
Figure BDA0003315070490000034
表示新矩阵M'中第i行第j列元素对应的相位的梯度，
Figure BDA0003315070490000035
表示新矩阵M'中第i行第j列元素对应的相位。
进一步地，所述步骤4中，种子生长的具体实施过程为：
步骤4.1：提取出所述种子点的八领域内像素点，计算八领域内每个像素点位置对应的值与该种子点位置对应的值之间的距离；
步骤4.2：判断八领域内每个像素点对应的所述距离是否小于所述种子生长阈值，如果是，则将该像素点归入到种子点集合中；
步骤4.3：以所述种子点集合内的每个像素点作为新的种子点，重复步骤4.1～4.2，直到没有新的像素点归入到种子点集合中为止，所述种子点集合内所有像素点所构成的区域即为所述干涉条纹清晰区域。
优选地，所述种子生长阈值为所述相位标准偏差图的方差的一半，所述相位标准偏差图是一个矩阵，相位标准偏差图的方差是指对矩阵中所有元素计算方差。
进一步地，所述步骤6中，采用SVD分解对所述区域块进行秩1估计的具体公式为：
Figure BDA0003315070490000041
其中，A表示区域块，即待SVD分解的矩阵，
Figure BDA0003315070490000042
表示SVD分解得到的奇异矩阵，u1表示第一特征向量对应的左奇异向量，v1表示第一特征向量对应的右奇异向量，σ1＞σ2＞…＞σk＞…，σk(k＝2,3,…)均接近0。
优选地，所述步骤7中，相位解缠为一维相位解缠，具体计算公式为：
Figure BDA0003315070490000043
其中，
Figure BDA0003315070490000044
表示左奇异向量或右奇异向量中第i个相位解缠后的相位，
Figure BDA0003315070490000045
表示左奇异向量或右奇异向量的第一个相位，mod{}表示取余运算，
Figure BDA0003315070490000046
表示左奇异向量或右奇异向量中第i个待解缠相位的下一个相位，
Figure BDA0003315070490000047
表示左奇异向量或右奇异向量中第i个待解缠相位，s表示左奇异向量或右奇异向量中元素数量。
进一步地，所述步骤7中，采用高斯函数进行第一次定权，所述高斯函数的计算公式为：
Figure BDA0003315070490000048
其中，s表示左奇异向量或右奇异向量中的元素数量，由左奇异向量的相位或右奇异向量的相位进行相位解缠后的所有相位构成左奇异向量或右奇异向量的解缠相位向量，i表示解缠相位向量中的第i个位置，i＝0,1,2,…,s-1，G(i)表示第一次定权时解缠相位向量中第i个位置对应的权值。
优选地，所述步骤7中，采用残差进行第二次定权，残差计算公式为：
Figure BDA0003315070490000049
其中：
Figure BDA0003315070490000051
表示左奇异向量或右奇异向量中第i个相位解缠后的相位，
Figure BDA0003315070490000052
表示将i代入第一次拟合直线的直线方程中所求得的值，由左奇异向量的相位或右奇异向量的相位进行相位解缠后的所有相位构成左奇异向量或右奇异向量的解缠相位向量，i表示解缠相位向量中的第i个位置，i＝0,1,2,…,s-1，di表示对应残差。
进一步地，所述步骤8中，发生位移后影像相对于发生位移前影像在x方向和y方向的平移量的计算公式分别为：
Figure BDA0003315070490000053
其中，x0表示发生位移后影像相对于发生位移前影像在x方向的平移量，y0表示发生位移后影像相对于发生位移前影像在y方向的平移量，rows表示影像中像素点的总行数，cols表示影像中像素点的总列数，ku表示右奇异向量对应的第二次拟合直线的斜率，kv表示左奇异向量对应的第二次拟合直线的斜率。
本发明还提供一种相位相关抗噪位移估计设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如上所述相位相关抗噪位移估计方法。
本发明还提供一种存储介质，其上存储有计算机程序，该程序被处理器执行时实现如上所述相位相关抗噪位移估计方法。
有益效果
与现有技术相比，本发明的优点在于：
本发明所提供的一种相位相关抗噪位移估计方法，引入相位标准偏差，在进行真实干涉相位估计前，首先将信噪比非常低的高频区域自动地剔除，以免在进行干涉相位估计时对真实信号产生过多干扰，大大提高了抗噪性能；同时，在进行解缠后相位直线拟合时，引入高斯函数定权，给予干涉相位低频部分更高的权值，高频部分低的权值，既保证了观测值的数量，同时也保证了观测值的质量，在此基础上，使用残差进行二次定权拟合直线，进一步减少了异常值的干扰，从而大大提高了抗噪性能。
相较于现有方法，本发明具有很强的抗噪性能，可以在成像条件很差的情况下准确地估计出位移量。
附图说明
为了更清楚地说明本发明的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一个实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明实施例中相位相关抗噪位移估计方法的流程图；
图2是本发明实施例中发生位移前、后的两幅影像，其中(a)为发生位移前影像，(b)为发生位移后影像；
图3是本发明实施例中滑动窗口示意图；
图4是本发明实施例中干涉相位图以及干涉相位图对应的PSD图，其中(a)为干涉相位图，(b)为对应的PSD图；
图5是本发明实施例中干涉条纹清晰区域的相位图，其中(a)为干涉条纹清晰的矩形区域，(b)为采用SVD分解进一步去噪后的相位图；
图6是本发明实施例中相位解缠后对应的相位，其中(a)为左奇异向量对应相位，(b)为右奇异向量对应相位；
图7是本发明实施例中第二次拟合直线，其中(a)为左奇异向量对应相位的第二次拟合直线，(b)为右奇异向量对应相位的第二次拟合直线。
具体实施方式
下面结合本发明实施例中的附图，对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面以具体地实施例对本申请的技术方案进行详细说明。下面这几个具体的实施例可以相互结合，对于相同或相似的概念或过程可能在某些实施例不再赘述。
如图1所示，本实施例所提供的一种相位相关抗噪位移估计方法，包括以下步骤：
步骤1：读取发生位移前、后的影像，分别对两幅影像进行二维傅里叶变换得到对应的频谱图；对频谱图进行频谱中心化，使频谱图中低频部分分布于频谱图中心区域，高频部分分布于频谱图边缘区域。
对于一张完整影像，裁剪其左上角部分用来模拟发生位移前影像，裁剪其右下角部分用来模拟发生位移后影像，两次裁剪的位置不同实现了位移模拟，再给裁剪后的两幅影像添加高斯噪声、斑点噪声等以模拟恶劣成像条件，得到本实施例所要采用的两幅影像，如图2所示。本实施例中，拍摄距离为10m左右，x方向的位移为-5.5pixel，y方向的位移为-5.5pixel。本实施例中，像素坐标系是以发生位移前影像的左上角为坐标原点，以左上角指向右上角为x方向，以左上角指向左下角为y方向。
通过二维傅里叶变换将影像在空间域上的像素坐标变换到频率域的频谱坐标，但是二维傅里叶变换得到的频谱图是高频部分在中间，低频部分在四周边缘区域，为了把能量集中起来，利用二维傅里叶变换的平移特性对频谱图进行频谱中心化，使频谱图中低频部分分布于频谱图中心区域，高频部分分布于频谱图边缘区域。
步骤2：对经步骤1处理后的两幅频谱图分别进行加窗，再对加窗后的两幅频谱图进行相位干涉，得到干涉相位图，如图4(a)所示。
为了减少频谱泄漏，对信号进行加窗操作。傅里叶变换的频率分辨率主要受窗函数的主瓣宽度影响，而泄漏程度则依赖于主瓣和旁瓣的相对幅值大小。矩形窗有最小的主瓣宽度，但是在最常见的窗中，矩形窗的旁瓣最大。因此，矩形窗的频率分辨率最高，而频谱泄漏则最大。不同的窗函数就是在频率分辨率和频谱泄漏中作一个折中选择。同时考虑频谱分辨率和频谱泄漏，本实施例采用汉宁窗对每幅所述频谱图进行加窗，汉宁窗的计算公式为：
Figure BDA0003315070490000071
其中，Whn(i,j)表示频谱图第i行第j列位置处的汉宁窗函数值，(i,j)分别表示影像中某像素点的行列号，rows表示影像中像素点的总行数，cols表示影像中像素点的总列数。
本实施例中，干涉相位图的计算公式为：
Figure BDA0003315070490000072
其中，Q(u,v)表示干涉相位图(即两幅影像的相关矩阵)，F1(u,v)表示发生位移前影像经二维傅里叶变换后对应的频谱图，F2(u,v)表示发生位移后影像经二维傅里叶变换后对应的频谱图，
Figure BDA0003315070490000073
表示取F2(u,v)的共轭，x0表示发生位移后影像相对于发生位移前影像在x方向的位移，y0表示发生位移后影像相对于发生位移前影像在y方向的位移，u表示像素坐标系的横坐标，v表示像素坐标系的纵坐标。
步骤3：计算干涉相位图中每个元素对应的相位标准偏差，并由所有相位标准偏差生成相位标准偏差图。
干涉相位图中的每个元素对应一个相位标准偏差，本实施例利用滑动窗口计算相位标准偏差，具体实现过程为：
步骤3.1：定义滑动窗口大小w*w，其中w表示边长且为奇数；
步骤3.2：对干涉相位图中每个元素求取相位，得到与干涉相位图对应的相位矩阵M；干涉相位图实质上为一个复数矩阵；
步骤3.3：构造与相位矩阵M维度一致(或大小一致)的相位标准偏差矩阵P，对相位标准偏差矩阵P中的所有元素进行初始化，相位标准偏差矩阵P中所有元素可以初始化为0或其他值；
步骤3.4：对于相位矩阵M中的每个元素(边界w/2位置的元素除外，从相位矩阵M第w/2行、第w/2列开始，如图3所示)，以该元素为中心，取w*w滑动窗口区域内的所有元素构成一个新矩阵M'，按照以下公式计算对应位置的相位标准偏差：
Figure BDA0003315070490000081
其中，
Figure BDA0003315070490000082
表示对应位置的相位标准偏差，对应位置是指滑动窗口的中心位置(即窗口中心位置对应到干涉相位图中相应位置的相位标准偏差)，
Figure BDA0003315070490000083
表示新矩阵M'中第i行第j列元素对应的相位的梯度(该梯度是指第i行第j列对应的行方向、列方向梯度均值)，
Figure BDA0003315070490000084
表示新矩阵M'中第i行第j列元素对应的相位。
当滑动窗口大小为3*3时，如图3所示，左边的滑动窗口的中心位置为A，所求得的相位标准偏差即为位置A的相位标准偏差；右边的滑动窗口的中心位置为B，所求得的相位标准偏差即为位置B的相位标准偏差；下边的滑动窗口的中心位置为C，所求得的相位标准偏差即为位置C的相位标准偏差。
引入相位标准偏差对干涉相位图进行评价，从而剔除干涉效果很差的区域，以增加干涉相位图的置信度，减少影像噪声的干扰，计算得到的相位标准偏差图如图4(b)所示。
步骤4：以相位标准偏差图的中心元素作为种子点，设置种子生长阈值，进行种子生长，得到干涉相位图中干涉条纹清晰区域。
本实施例中，种子生长的具体实施过程为：
步骤4.1：提取出种子点的八领域内像素点(即在相位标准偏差图中，与种子点相邻的八个元素)，计算八领域内每个像素点位置对应的值与该种子点位置对应的值之间的距离；
步骤4.2：判断八领域内每个像素点对应的距离是否小于种子生长阈值，如果是，则将该像素点归入到种子点集合中；
步骤4.3：以种子点集合内的每个像素点作为新的种子点，重复步骤4.1～4.2，直到没有新的像素点归入到种子点集合中为止，种子点集合内所有像素点所构成的区域即为干涉条纹清晰区域，如图5(a)所示。
本实施例中，种子生长阈值为相位标准偏差图的方差的一半，相位标准偏差图是一个矩阵，相位标准偏差图的方差是指对矩阵中所有元素计算方差。
步骤5：获取干涉条纹清晰区域的最小外接矩阵，对该最小外接矩阵所包含的区域进行裁剪，得到相应区域块。
步骤6：采用SVD分解对步骤5的区域块进行秩1估计，提取SVD分解后最大特征值对应的左奇异向量和右奇异向量，并分别计算左奇异向量和右奇异向量的相位。
采用SVD分解对区域块进行秩1估计的具体公式为：
Figure BDA0003315070490000091
其中，A表示区域块，即待SVD分解的矩阵，
Figure BDA0003315070490000092
表示SVD分解得到的奇异矩阵，u1表示第一特征向量对应的右奇异向量，v1表示第一特征向量对应的右奇异向量，即σ1对应的两个奇异向量，σ1＞σ2＞…＞σk＞…，σk(k＝2,3,…)均接近0。u1,u2,u3,…均是SVD分解得到的左奇异向量，v1,v2,v3…均是SVD分解得到的右奇异向量，本实施例仅取σ1对应的u1、v1进行后续相位解缠步骤。
利用SVD分解对区域块(即相关矩阵)进行秩1估计，将二维相位转换成一维，降低了相位解缠的难度。采用SVD分解进一步去噪后的相位图如图5(b)所示。
步骤7：对左奇异向量和右奇异向量的相位分别进行相位解缠，对左奇异向量和右奇异向量相位解缠后的相位分别进行第一次定权，对第一次定权后的相位及其对应的位置序号进行直线拟合，得到第一次拟合直线；根据所述第一次拟合直线模型的残差，对左奇异向量和右奇异向量相位解缠后的相位进行第二次定权，对第二次定权后的相位及其对应的位置序号进行直线拟合，得到第二次拟合直线。
本实施例中，相位解缠为一维相位解缠，具体计算公式为：
Figure BDA0003315070490000101
其中，
Figure BDA0003315070490000102
表示左奇异向量或右奇异向量中第i个相位解缠后的相位，
Figure BDA0003315070490000103
表示左奇异向量或右奇异向量的第一个相位，mod{}表示取余运算，
Figure BDA0003315070490000104
表示左奇异向量或右奇异向量中第i个待解缠相位的下一个相位，
Figure BDA0003315070490000105
表示左奇异向量或右奇异向量中第i个待解缠相位，s表示左奇异向量或右奇异向量中元素数量，根据式(4)，左奇异向量的s对应m，右奇异向量的s对应n。
如图6所示相位解缠后对应的相位，相位解缠后的相位是一个个离散点，每个相位离散点均对应一个残差，对相位离散点以及相位离散点在左奇异向量或右奇异向量中的位置进行直线拟合，在直线拟合前，采用残差进行定权，利用带权的解缠相位直线拟合进行位移估计。本实施例中，采用高斯函数进行第一次定权，高斯函数的计算公式为：
Figure BDA0003315070490000106
其中，s表示左奇异向量或右奇异向量中的元素数量，由左奇异向量的相位或右奇异向量的相位进行相位解缠后的所有相位构成左奇异向量或右奇异向量的解缠相位向量，i表示解缠相位向量中的第i个位置，i＝0,1,2,…,s-1，G(i)表示第一次定权时解缠相位向量中第i个位置对应的权值，第一次定权过程就是G(i)*第i个位置对应的相位的过程。
对第一次定权后的相位和该相位在左奇异向量或右奇异向量中对应的位置序号进行直线拟合，得到第一次拟合直线。
根据第一次拟合直线，采用残差对左奇异向量和右奇异向量相位解缠后的相位进行第二次定权，残差计算公式为：
Figure BDA0003315070490000107
其中：
Figure BDA0003315070490000108
表示左奇异向量或右奇异向量中第i个相位解缠后的相位，
Figure BDA0003315070490000109
表示将i代入第一次拟合直线的直线方程中所求得的值，由左奇异向量的相位或右奇异向量的相位进行相位解缠后的所有相位构成左奇异向量或右奇异向量的解缠相位向量，i表示解缠相位向量中的第i个位置，i＝0,1,2,…,s-1，di表示对应残差。
对第二次定权后的相位和该相位在左奇异向量或右奇异向量中对应的位置序号进行直线拟合，得到第二次拟合直线，即得到左奇异向量对应的第二次拟合直线，右奇异向量对应的第二次拟合直线，如图7所示。
步骤8：根据左奇异向量和右奇异向量对应的第二次拟合直线的斜率，再结合影像的长和宽，计算发生位移后影像相对于发生位移前影像在x方向和y方向的平移量。
发生位移后影像相对于发生位移前影像在x方向和y方向的平移量的计算公式分别为：
Figure BDA0003315070490000111
其中，x0表示发生位移后影像相对于发生位移前影像在x方向的平移量，y0表示发生位移后影像相对于发生位移前影像在y方向的平移量，rows表示影像中像素点的总行数，cols表示影像中像素点的总列数，ku表示右奇异向量对应的第二次拟合直线的斜率，kv表示左奇异向量对应的第二次拟合直线的斜率。
具体计算结果如表1所示，表1还给出几种经典抗噪相位相关位移估计方法，由表1可知，在拍摄影像噪声干扰严重的情况下，本发明可以高精度第估计出场景中的真实位移。
表1不同估计方法的对比
Figure BDA0003315070490000112
本实施例还提出了一种相位相关抗噪位移估计设备，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上所述相位相关抗噪位移估计方法的步骤。
示例性的，所述计算机程序可以被分割成一个或多个模块/单元，所述一个或者多个模块/单元被存储在所述存储器中，并由所述处理器执行，以完成本发明。所述一个或多个模块/单元可以是能够完成特定功能的一系列计算机程序指令段，该指令段用于描述所述计算机程序在所述计算机设备中的执行过程。
所述计算机设备可以是桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备。所述计算机设备可包括，但不仅限于，处理器、存储器。本领域技术人员可以理解，相位相关抗噪位移估计设备仅仅是设备的示例，并不构成对设备的限定，可以包括比设备更多或更少的部件，或者组合某些部件，或者不同的部件，例如所述设备还可以包括输入输出设备、网络接入设备、总线等。
所述处理器可以是中央处理单元(Central Processing Unit，CPU)，还可以是其他通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现成可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
所述存储器可用于存储所述计算机程序和/或模块，所述处理器通过运行或执行存储在所述存储器内的计算机程序和/或模块，以及调用存储在存储器内的数据，实现所述指针式计量器具表盘图像转换系统的各种功能。所述存储器可主要包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序(比如声音播放功能、图像播放功能等)等；存储数据区可存储根据手机的使用所创建的数据(比如音频数据、电话本等)等。此外，存储器可以包括高速随机存取存储器，还可以包括非易失性存储器，例如硬盘、内存、插接式硬盘，智能存储卡(Smart Media Card，SMC)，安全数字(SecureDigital，SD)卡，闪存卡(Flash Card)、至少一个磁盘存储器件、闪存器件、或其他易失性固态存储器件。
所述计算机程序被处理器执行时实现所述相位相关抗噪位移估计方法的步骤。
所述相位相关抗噪位移估计方法设备集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实现上述实施例方法中的全部或部分流程，也可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质可以包括：能够携带所述计算机程序代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，Random Access Memory)、电载波信号、电信信号以及软件分发介质等。需要说明的是，所述计算机可读介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括电载波信号和电信信号。
以上所揭露的仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或变型，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种相位相关抗噪位移估计方法，其特征在于，包括以下步骤：
步骤1：读取发生位移前、后的影像，分别对两幅影像进行二维傅里叶变换得到对应的频谱图；对所述频谱图进行频谱中心化，使频谱图中低频部分分布于频谱图中心区域，高频部分分布于频谱图边缘区域；
步骤2：对经所述步骤1处理后的两幅频谱图分别进行加窗，再对加窗后的两幅频谱图进行相位干涉，得到干涉相位图；
步骤3：计算所述干涉相位图中每个元素对应的相位标准偏差，并由所有相位标准偏差生成相位标准偏差图；
步骤4：以所述相位标准偏差图的中心元素作为种子点，设置种子生长阈值，进行种子生长，得到所述干涉相位图中干涉条纹清晰区域；
步骤5：获取所述干涉条纹清晰区域的最小外接矩阵，对所述最小外接矩阵所包含的区域进行裁剪，得到相应区域块；
步骤6：采用SVD分解对所述区域块进行秩1估计，提取SVD分解后最大特征值对应的左奇异向量和右奇异向量，并分别计算所述左奇异向量和右奇异向量的相位；
步骤7：对所述左奇异向量和右奇异向量的相位分别进行相位解缠，对左奇异向量和右奇异向量相位解缠后的相位分别进行第一次定权，对第一次定权后的相位及其对应的位置序号进行直线拟合，得到第一次拟合直线；
根据所述第一次拟合直线模型的残差，对左奇异向量和右奇异向量相位解缠后的相位进行第二次定权，对第二次定权后的相位及其对应的位置序号进行直线拟合，得到第二次拟合直线；
步骤8：根据左奇异向量和右奇异向量对应的第二次拟合直线的斜率，再结合所述影像的长和宽，计算发生位移后影像相对于发生位移前影像在x方向和y方向的平移量。
2.如权利要求1所述的相位相关抗噪位移估计方法，其特征在于，所述步骤2中，采用汉宁窗对每幅所述频谱图进行加窗，所述汉宁窗的计算公式为：
Figure FDA0003315070480000011
其中，Whn(i,j)表示频谱图第i行第j列位置处的汉宁窗函数值，(i,j)分别表示影像中某像素点的行列号，rows表示影像中像素点的总行数，cols表示影像中像素点的总列数。
3.如权利要求1所述的相位相关抗噪位移估计方法，其特征在于，所述步骤2中，干涉相位图的计算公式为：
Figure FDA0003315070480000021
其中，Q(u,v)表示干涉相位图，F1(u,v)表示发生位移前影像经二维傅里叶变换后对应的频谱图，F2(u,v)表示发生位移后影像经二维傅里叶变换后对应的频谱图，
Figure FDA0003315070480000022
表示取F2(u,v)的共轭，x0表示发生位移后影像相对于发生位移前影像在x方向的位移，y0表示发生位移后影像相对于发生位移前影像在y方向的位移，u表示像素坐标系的横坐标，v表示像素坐标系的纵坐标。
4.如权利要求1所述的相位相关抗噪位移估计方法，其特征在于，所述步骤3中，计算相位标准偏差的具体实现过程为：
步骤3.1：定义滑动窗口大小w*w，其中w表示边长且为奇数；
步骤3.2：对所述干涉相位图中每个元素求取相位，得到对应的相位矩阵M；
步骤3.3：构造与所述相位矩阵M维度一致的相位标准偏差矩阵P，对所述相位标准偏差矩阵P中的所有元素进行初始化；
步骤3.4：对于所述相位矩阵M中的每个元素，以该元素为中心，取w*w滑动窗口大小区域内的所有元素构成一个新矩阵M'，按照以下公式计算对应位置的相位标准偏差：
Figure FDA0003315070480000023
其中，
Figure FDA0003315070480000024
表示对应位置的相位标准偏差，对应位置是指滑动窗口的中心位置，
Figure FDA0003315070480000025
表示新矩阵M'中第i行第j列元素对应的相位的梯度，
Figure FDA0003315070480000026
表示新矩阵M'中第i行第j列元素对应的相位。
5.如权利要求1所述的相位相关抗噪位移估计方法，其特征在于，所述步骤4中，种子生长的具体实施过程为：
步骤4.1：提取出所述种子点的八领域内像素点，计算八领域内每个像素点位置对应的值与该种子点位置对应的值之间的距离；
步骤4.2：判断八领域内每个像素点对应的所述距离是否小于所述种子生长阈值，如果是，则将该像素点归入到种子点集合中；
步骤4.3：以所述种子点集合内的每个像素点作为新的种子点，重复步骤4.1～4.2，直到没有新的像素点归入到种子点集合中为止，所述种子点集合内所有像素点所构成的区域即为所述干涉条纹清晰区域；
优选地，所述种子生长阈值为所述相位标准偏差图的方差的一半，所述相位标准偏差图是一个矩阵，相位标准偏差图的方差是指对矩阵中所有元素计算方差。
6.如权利要求1所述的相位相关抗噪位移估计方法，其特征在于，所述步骤6中，采用SVD分解对所述区域块进行秩1估计的具体公式为：
Figure FDA0003315070480000031
其中，A表示区域块，即待SVD分解的矩阵，
Figure FDA0003315070480000032
表示SVD分解得到的奇异矩阵，u1表示第一特征向量对应的左奇异向量，v1表示第一特征向量对应的右奇异向量，σ1＞σ2＞…＞σk＞…，σk(k＝2,3,…)均接近0；
优选地，所述步骤7中，相位解缠为一维相位解缠，具体计算公式为：
Figure FDA0003315070480000033
其中，
Figure FDA0003315070480000034
表示左奇异向量或右奇异向量中第i个相位解缠后的相位，
Figure FDA0003315070480000035
表示左奇异向量或右奇异向量的第一个相位，mod{}表示取余运算，
Figure FDA0003315070480000036
表示左奇异向量或右奇异向量中第i个待解缠相位的下一个相位，
Figure FDA0003315070480000037
表示左奇异向量或右奇异向量中第i个待解缠相位，s表示左奇异向量或右奇异向量中元素数量。
7.如权利要求1所述的相位相关抗噪位移估计方法，其特征在于，所述步骤7中，采用高斯函数进行第一次定权，所述高斯函数的计算公式为：
Figure FDA0003315070480000038
其中，s表示左奇异向量或右奇异向量中的元素数量，由左奇异向量的相位或右奇异向量的相位进行相位解缠后的所有相位构成左奇异向量或右奇异向量的解缠相位向量，i表示解缠相位向量中的第i个位置，i＝0,1,2,…,s-1，G(i)表示第一次定权时解缠相位向量中第i个位置对应的权值；
优选地，所述步骤7中，采用残差进行第二次定权，残差计算公式为：
Figure FDA0003315070480000041
其中：
Figure FDA0003315070480000042
表示左奇异向量或右奇异向量中第i个相位解缠后的相位，
Figure FDA0003315070480000043
表示将i代入第一次拟合直线的直线方程中所求得的值，由左奇异向量的相位或右奇异向量的相位进行相位解缠后的所有相位构成左奇异向量或右奇异向量的解缠相位向量，i表示解缠相位向量中的第i个位置，i＝0,1,2,…,s-1，di表示对应残差。
8.如权利要求1～7中任一项所述的相位相关抗噪位移估计方法，其特征在于，所述步骤8中，发生位移后影像相对于发生位移前影像在x方向和y方向的平移量的计算公式分别为：
Figure FDA0003315070480000044
其中，x0表示发生位移后影像相对于发生位移前影像在x方向的平移量，y0表示发生位移后影像相对于发生位移前影像在y方向的平移量，rows表示影像中像素点的总行数，cols表示影像中像素点的总列数，ku表示右奇异向量对应的第二次拟合直线的斜率，kv表示左奇异向量对应的第二次拟合直线的斜率。
9.一种相位相关抗噪位移估计设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其特征在于：所述处理器执行所述程序时实现如权利要求1～8中任一项所述相位相关抗噪位移估计方法。
10.一种存储介质，其上存储有计算机程序，其特征在于：该程序被处理器执行时实现如权利要求1～8中任一项所述相位相关抗噪位移估计方法。
CN202111228184.3A 2021-10-21 2021-10-21 相位相关抗噪位移估计方法、设备及存储介质 Pending CN114066749A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111228184.3A CN114066749A (zh) 2021-10-21 2021-10-21 相位相关抗噪位移估计方法、设备及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111228184.3A CN114066749A (zh) 2021-10-21 2021-10-21 相位相关抗噪位移估计方法、设备及存储介质

Publications (1)

Publication Number Publication Date
CN114066749A true CN114066749A (zh) 2022-02-18

Family

ID=80235209

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111228184.3A Pending CN114066749A (zh) 2021-10-21 2021-10-21 相位相关抗噪位移估计方法、设备及存储介质

Country Status (1)

Country Link
CN (1) CN114066749A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115937282A (zh) * 2023-01-10 2023-04-07 郑州思昆生物工程有限公司 一种荧光图像的配准方法、装置、设备和存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103824286A (zh) * 2014-02-14 2014-05-28 同济大学 一种svd-ransac亚像素相位相关匹配方法
CN110657742A (zh) * 2019-09-30 2020-01-07 广东工业大学 含水层形变信号分离方法、装置、设备及可读存储介质
US10740341B1 (en) * 2017-04-04 2020-08-11 Cyber Atomics, Inc. Single- and multi-variate tensor spectral analysis

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103824286A (zh) * 2014-02-14 2014-05-28 同济大学 一种svd-ransac亚像素相位相关匹配方法
US10740341B1 (en) * 2017-04-04 2020-08-11 Cyber Atomics, Inc. Single- and multi-variate tensor spectral analysis
CN110657742A (zh) * 2019-09-30 2020-01-07 广东工业大学 含水层形变信号分离方法、装置、设备及可读存储介质

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115937282A (zh) * 2023-01-10 2023-04-07 郑州思昆生物工程有限公司 一种荧光图像的配准方法、装置、设备和存储介质

Similar Documents

Publication Publication Date Title
US9262808B2 (en) Denoising of images with nonstationary noise
Sandić-Stanković et al. DIBR synthesized image quality assessment based on morphological wavelets
Yuan et al. Factorization-based texture segmentation
CN111353947A (zh) 磁共振并行成像方法及相关设备
CN111783583B (zh) 基于非局部均值算法的sar图像相干斑抑制方法
CN114677300A (zh) 一种基于双阶段学习框架的高光谱图像深度降噪的方法及系统
Park 2D discrete Fourier transform on sliding windows
Patel et al. Separated component-based restoration of speckled SAR images
Farouj et al. Hyperbolic Wavelet-Fisz denoising for a model arising in Ultrasound Imaging
CN111079893B (zh) 用于干涉条纹图滤波的生成器网络的获取方法和装置
Gao et al. A novel two-step noise reduction approach for interferometric phase images
CN114066749A (zh) 相位相关抗噪位移估计方法、设备及存储介质
CN114155161B (zh) 图像去噪方法、装置、电子设备与存储介质
CN113935246A (zh) 一种信号鲁棒稀疏时频分析方法、终端设备及存储介质
He et al. A parallel alternating direction method with application to compound l1-regularized imaging inverse problems
CN112614081A (zh) 干涉图去噪的方法
CN112435211A (zh) 一种内窥镜图像序列中轮廓稠密特征点描述和匹配的方法
CN116596801A (zh) 一种图像非局部均值去噪方法及装置
Das et al. A concise review of fast bilateral filtering
KR100633555B1 (ko) 신호 복구 방법
Mastriani Denoising based on wavelets and deblurring via self-organizing map for Synthetic Aperture Radar images
Joseph et al. A Novel Denoising Algorithm Based on Superpixel Clustering and Dictionary Learning Approach.
Shi et al. Image denoising through locally linear embedding
Wagenmaker et al. Robust surface reconstruction from gradients via adaptive dictionary regularization
Wang et al. A fast Lee filter algorithm based on Fourier transform

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination