CN108090878B - 基于差异图和补偿滤波的干涉相位滤波方法 - Google Patents

基于差异图和补偿滤波的干涉相位滤波方法 Download PDF

Info

Publication number
CN108090878B
CN108090878B CN201711304057.0A CN201711304057A CN108090878B CN 108090878 B CN108090878 B CN 108090878B CN 201711304057 A CN201711304057 A CN 201711304057A CN 108090878 B CN108090878 B CN 108090878B
Authority
CN
China
Prior art keywords
interferometric phase
image
phase image
remaining
pixel
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.)
Active
Application number
CN201711304057.0A
Other languages
English (en)
Other versions
CN108090878A (zh
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.)
Hunan Tripod Quantum Technology Co Ltd
Original Assignee
Hunan Tripod Quantum Technology Co Ltd
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 Hunan Tripod Quantum Technology Co Ltd filed Critical Hunan Tripod Quantum Technology Co Ltd
Priority to CN201711304057.0A priority Critical patent/CN108090878B/zh
Publication of CN108090878A publication Critical patent/CN108090878A/zh
Application granted granted Critical
Publication of CN108090878B publication Critical patent/CN108090878B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G06T5/70—
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/10—Image enhancement or restoration by non-spatial domain filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20021—Dividing image into blocks, subimages or windows

Abstract

本发明提出一种基于差异图和补偿滤波的干涉相位滤波方法，先对输入的干涉相位图进行基于差异图的频域滤波提取干涉相位的主要分量，从而降低了残余干涉相位图中干涉条纹的密度，便于后续处理。之后，对残余干涉相位图进行补偿滤波。补偿滤波采用以残余干涉相位图的相位导数方差图作为自适应滤波参数的BM3D滤波，能更好的提取干涉相位的细节分量。本发明包含频域滤波和空域滤波，适合处理含有密集干涉条纹的干涉相位，能在有效滤除噪声的同时不破坏密集干涉条纹的连续性。

Description

基于差异图和补偿滤波的干涉相位滤波方法
技术领域
本发明属于遥感和信号处理的交叉技术领域，特别涉及一种适用于含有密集干涉条纹的干涉相位滤波方法。
背景技术
对两幅合成孔径雷达图像共轭相乘，并对相乘结果取相位便得到干涉相位。干涉相位是合成孔径雷达干涉测量中重要的物理量，其质量的好坏将决定最终产品——数字高程模型或地形形变量的精度。然而，受去相关因素的影响，干涉相位中总是存在严重的空变噪声。空变噪声不仅会引入残差点还会破坏干涉相位的分布，从而增加了后续相位解缠的难度，最终导致产品精度的降低，因此必须通过干涉相位滤波过程予以滤除。
BM3D滤波算法是一个基于图像非局部相似性的高性能滤波算法，它将结构相似的二维图像块组成三维数组，然后对数组进行联合滤波，从而达到了经典图像滤波算法无法企及的效果。
但是将BM3D滤波算法应用于干涉相位滤波会存在2个问题：(1)对于密集干涉条纹，由于BM3D滤波算法的滤波窗口可能会跨越多个干涉条纹，从而破坏干涉条纹的连续性；(2)BM3D滤波算法需要能表征图像噪声方差的自适应参数作为滤波输入参数，不同的滤波输入参数将导致不同的滤波结果，但现有的BM3D滤波算法并未给出该滤波输入参数的计算方法。
发明内容
针对BM3D滤波算法应用于干涉相位滤波时存在的2个问题，本发明提出一种基于差异图和补偿滤波的干涉相位滤波方法。
本发明先对输入的干涉相位图进行基于差异图的频域滤波提取干涉相位的主要分量(密集干涉条纹)，从而降低了残余干涉相位图中干涉条纹的密度，便于后续处理。之后，对残余干涉相位图进行补偿滤波。补偿滤波采用以残余干涉相位图的相位导数方差图作为自适应滤波参数的BM3D滤波，能更好的提取干涉相位的细节分量。因此，本发明即能有效去噪又能保证密集干涉条纹的连续性。
基于差异图和补偿滤波的干涉相位滤波方法，包括以下步骤：
S1.对于待滤波的干涉相位图，计算干涉相位图的差异图；
设待滤波的干涉相位图即输入的干涉相位图的大小为Nx×Ny像素，将干涉相位图像素的像素值([-π,π])量化成13个值，量化采用为间隔的均匀量化，得到量化后的像素值集合G。最后得到一个新的映射关系：f:X×Y→G，其中，X＝{1,2,3,…,Nx}，Y＝{1,2,3,…,Ny}，
干涉相位图的差异图DM(Gi,Gj)的计算公式为：
其中，i是干涉相位图中任一像素点，j是以像素点i为中心且与像素点i相邻的8个像素点中的任一像素点，Gi是像素点i量化后的像素值；Gj是像素点j量化后的像素值，p(Gi,Gj)是干涉相位图的归一化灰度共生矩阵GLCM(Gi,Gj)；
其中GLCM(Gi,Gj)的计算方法如下：
P(Gi,Gj,0°和180°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,k-m＝0,l-n|＝1}
P(Gi,Gj,45°和225°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,(k-m＝1,l-n＝-1)or(k-m＝-1,l-n＝1)}
P(Gi,Gj,90°和270°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,k-m|＝1,l-n＝0}
P(Gi,Gj,135°和315°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,(k-m＝1,l-n＝1)or(k-m＝-1,l-n＝-1)}
其中，#{·}是计数算子，(k,l)和(m,n)是干涉相位图中像素值分别为Gi和Gj(Gi,Gj∈G)对应的两个像素点的坐标。
基于计算得到的频数可以得到下面的矩阵：
其中，θ∈{0°和180° 45°和225° 90°和270° 135°和315°}。
因此，可以得到4个矩阵Pθ，将4个矩阵求平均就得到灰度共生矩阵GLCM(Gi,Gj)。
S2.将干涉相位图转换成复干涉相位图，并对复干涉相位图进行分块；
其中，是干涉相位图，ψ是复干涉相位图。
对复干涉相位图进行分块，将其分成尺寸相同的多个正方形图像块，且列方向的相邻正方形图像块之间有重叠区域，行方向的相邻正方形图像块之间也有重叠区域。各正方形图像块在列方向与其相邻的正方形图像块之间重叠的区域大小与在行方向与其相邻的正方形图像块之间重叠的区域大小相同。在对复干涉相位图进行分块时，对于复干涉相位图的边缘区域，如出现待分块的行或者列不足分块的正方形图像块大小时，通过扩展将待分块的行或者列补足像素以满足正方形图像块大小(可以采用镜像、周期延拓等扩展方式，不同扩展方式造成的滤波结果差异可以忽略不计)，待滤波结束后截取原有效区域(去除扩展区域)即可。
正方形图像块尺寸的选择会影响到滤波的速度和精度：尺寸大效率高但精度低；尺寸小效率低精度高。本发明中取正方形图像块尺寸为32×32像素，该尺寸能在效率和精度之间进行较好的折中。
S3.对S2分块得到的各正方形图像块进行图像块滤波；
S3.1对各正方形图像块进行二维傅里叶正变换得到对应的正方形图像块的二维频谱：
B＝FFT2(P)
其中，P是正方形图像块，B是正方形图像块的二维频谱，FFT2(·)二维傅里叶正变换。
S3.2对各正方形图像块的二维频谱进行取绝对值、平滑、归一化和以正方形图像块对应的干涉相位图的差异图的均值为幂指数的指数操作得到与其对应的权重图像块。
计算正方形图像块二维频谱的绝对值A：
A＝abs(B)
其中，abs(·)表示取绝对值。然后对正方形图像块二维频谱的绝对值A进行平滑：
S＝smooth(A)
其中，S是平滑后的正方形图像块二维频谱的绝对值，smooth(·)是二维高斯滤波器且其窗口大小根据实际情况确定，本发明实施方式中为7×7像素。接着对S进行归一化：
N＝normalization(S)
其中，normalization(·)是归一化算子。最后，对N取幂指数得到权重图像块W：
W＝Nα
其中：
其中，是正方形图像块对应的干涉相位图的差异图的均值，计算公式为：
其中，mean(·)是取均值算子，DMw是正方形图像块对应的干涉相位图的差异图。
S3.3计算各正方形图像块其对应的权重图像块与正方形图像块二维频谱的Hadamard积(矩阵对应元素相乘)，从而完成干涉相位主要分量提取：
BW＝B·W
BW是权重图像块与正方形图像块二维频谱的Hadamard积。
S4.计算频域滤波后的复干涉相位图；
得到所有正方形图像块二维频谱与其对应权重图像块的Hadamard积后，对所有
Hadamard积进行二维傅里叶逆变换，得到频域滤波后的复干涉相位图像块PW：
PW＝IFFT2(BW)
其中，IFFT2(·)是二维傅里叶逆变换。
正方形图像块之间重叠区域的频域滤波结果来自不同的正方形图像块，因此采用对所有重叠区域的频域滤波结果进行取均值的方法，从而得到频域滤波后的复干涉相位图ψfrequency。
S5.计算空域滤波后的残余复干涉相位图；
S5.1计算残余复干涉相位图ψresidual，残余干涉相位图残余干涉相位图的相位导数方差图VM：
ψresidual＝ψ-ψfrequency
其中，angle(·)是取相位算子；
残余干涉相位图的相位导数方差图VM的计算方法如下：
残余干涉相位图的相位导数方差图VM的大小和残余干涉相位图的大小一样，即为Nx×Ny像素。对于残余干涉相位图中的每一个像素点(s,t)((s,t)∈Nx×Ny)，可以计算该像素点的相位导数方差值VM(s,t)：
其中，和的计算公式为：
在计算残余干涉相位图边缘像素点即残余干涉相位图第一行、残余干涉相位图最后一行、残余干涉相位图第一列和残余干涉相位图最后一列上各像素点的相位导数方差值时，需要先将残余干涉相位图由Nx×Ny像素扩展成(Nx+2)×(Ny+2)像素，即在残余干涉相位图第一列的左边和最后一列的右边各增加一列、在残余干涉相位图第一行的上部和最后一行的下部各增加一行，扩展可以采用镜像、周期延拓等扩展方式，不同扩展方式造成的滤波结果差异可以忽略不计。
当得到残余干涉相位图中每一个像素点的相位导数方差值后，可以得到残余干涉相位图的相位导数方差图VM：
S5.2计算残余复干涉相位图的实部图像和虚部图像：
ψresidual_real＝real(ψresidual)
ψresidual_imag＝imag(ψresidual)
其中，ψresidual_real是残余复干涉相位图的实部图像，ψresidual_imag是残余复干涉相位图的虚部图像，real(·)是取实部算子，imag(·)是取虚部算子。
S5.3计算空域滤波后的残余复干涉相位图ψspatial的实部图像ψspatial_real和虚部图像ψspatial_imag：
其中，BM3D(·)是BM3D滤波器(Kostadin Dabov,Alessandro Foi,VladimirKatkovnik,and Karen Egiazarian,Image Denoising by Sparse 3-D Transform-DomainCollaborative Filtering,IEEE Transactions on Image Processing,vol.16,no.8,2007)，是残余干涉相位图的相位导数方差图的均值，即
S5.4计算空域滤波后的残余复干涉相位图ψspatial：
S6.计算滤波后的干涉相位图。
S6.1计算滤波后的复干涉相位图ψfiltered：
ψfiltered＝ψspatial+ψfrequency
S6.2计算滤波后的干涉相位图：
本发明具有以下技术效果：
本发明包含频域滤波和空域滤波，适合处理含有密集干涉条纹的干涉相位，能在有效滤除噪声的同时不破坏密集干涉条纹的连续性。
本发明的创新点在于：
1、频域滤波方面，采用干涉相位的差异图来确定频域滤波中的幂指数。“一种干涉相位滤波方法(专利申请号为201510853164.3)”中频域滤波的幂指数是由复干涉相位图的噪声标准差和干涉相位图的归一化倒置相位导数方差图复合而成。比较而言，噪声标准差和归一化倒置相位导数方差图只包含了图像噪声水平的大小，而差异图是由灰度共生矩阵衍生而来，它不仅包含了图像噪声水平的大小，还包含了图像的方向、幅度变化和纹理等综合信息。干涉相位与地形对应，因此干涉相位具有很强的方向性和纹理特征，所以差异图在提取干涉相位的主要分量上更加合理有效。
2、空域滤波方面，采用了以残余干涉相位图的相位导数方差图为自适应滤波参数的BM3D滤波。“一种干涉相位滤波方法(专利申请号为201510853164.3)”中的空域滤波采用的是固定滤波参数的均值滤波，去噪能力不能随残余干涉相位图中的噪声水平自适应变化，在噪声较强时不能有效去除噪声，而在噪声较弱时易损失掉相位细节信息；而残余干涉相位图的相位导数方差图是随残余干涉相位图中的噪声水平自适应变化的，因此本发明的空域滤波是自适应的，能有效去噪同时减少相位细节信息的损失。
附图说明
图1是本发明的流程图；
图2是复干涉相位图分块示意图；
图3是待滤波的干涉相位图；
图4为BM3D方法得到的滤波后的干涉相位图；
图5为本发明得到的滤波后的干涉相位图；
图6为图3，图4和图5的性能对比图。
具体实施方式
参照图1，为本发明基于差异图和补偿滤波的干涉相位滤波方法的流程图，包括以下步骤：
S1.计算干涉相位图的差异图；
设待滤波的干涉相位图(即输入的干涉相位图)的大小为Nx×Ny像素，将干涉相位图像素的像素值([-π,π])量化成13个值，量化采用为间隔的均匀量化，得到量化后的像素值集合G。最后得到一个新的映射关系：f:X×Y→G，其中，X＝{1,2,3,…,Nx}，Y＝{1,2,3,…,Ny}，干涉相位图的差异图DM(Gi,Gj)的计算公式为：
其中，i是干涉相位图中任一像素点，j是以像素点i为中心且与像素点i相邻的8个像素点中的任一像素点，Gi是像素点i量化后的像素值；Gj是像素点j量化后的像素值，p(Gi,Gj)是干涉相位图的归一化灰度共生矩阵GLCM(Gi,Gj)。
其中GLCM(Gi,Gj)的计算方法如下：
分别计算干涉相位图中每一个像素点与其相邻8个像素点在0°和180°方向上的频数P(Gi,Gj,0°和180°)、在45°和225°方向上的频数P(Gi,Gj,45°和225°)、在90°和270°方向上的频数P(Gi,Gj,90°和270°)以及在135°和315°方向上的频数P(Gi,Gj,135°和315°)：
P(Gi,Gj,0°和180°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,k-m＝0,l-n|＝1}
P(Gi,Gj,45°和225°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,(k-m＝1,l-n＝-1)or(k-m＝-1,l-n＝1)}
P(Gi,Gj,90°和270°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,k-m|＝1,l-n＝0}
P(Gi,Gj,135°和315°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,(k-m＝1,l-n＝1)or(k-m＝-1,l-n＝-1)}
其中，#{·}是计数算子，(k,l)和(m,n)是干涉相位图中像素值分别为Gi和Gj(Gi,Gj∈G)对应的两个像素点的坐标。
基于计算得到的频数可以得到下面的矩阵：
其中，θ∈{0°和180° 45°和225° 90°和270° 135°和315°}。
因此，可以得到4个矩阵Pθ，将4个矩阵求平均就得到灰度共生矩阵GLCM(Gi,Gj)。
S2.将干涉相位图转换成复干涉相位图，并对复干涉相位图进行分块；
其中，是干涉相位图，ψ是复干涉相位图。
对复干涉相位图进行分块，将其分成尺寸相同的多个正方形图像块，且列方向的相邻正方形图像块之间有重叠区域，行方向的相邻正方形图像块之间也有重叠区域。各正方形图像块在列方向与其相邻的正方形图像块之间重叠的区域大小与在行方向与其相邻的正方形图像块之间重叠的区域大小相同。在对复干涉相位图进行分块时，对于复干涉相位图的边缘区域，如出现待分块的行或者列不足分块的正方形图像块大小时，通过扩展将待分块的行或者列补足像素以满足正方形图像块大小(可以采用镜像、周期延拓等扩展方式，不同扩展方式造成的滤波结果差异可以忽略不计)，待滤波结束后截取原有效区域(去除扩展区域)即可。
正方形图像块尺寸的选择会影响到滤波的速度和精度：尺寸大效率高但精度低；尺寸小效率低精度高。本发明中取正方形图像块尺寸为32×32像素，该尺寸能在效率和精度之间进行较好的折中。
本实施例复干涉相位图分块和重叠区域的示意图如图2所示，图像左上角的O(1,1)点为坐标原点，图中简单显示了4个32×32像素大小的正方形图像块。复干涉相位图分块的重叠在行方向和列方向进行：每一行的正方形图像块在列方向重叠32×14像素大小的区域，每一列的正方形图像块在行方向重叠14×32像素大小的区域。当分块到复干涉相位图最后一行和最后一列时，可能会出现分块的行或列不足32像素的情况，可以通过扩展将行或列补足32像素(可以采用镜像、周期延拓等扩展方式，不同扩展方式造成的滤波结果差异可以忽略不计)，待滤波结束后截取原有效区域。分块结束后，对于每个正方形图像块完成步骤S3。
S3.对S2分块得到的各图像块进行图像块滤波；
S3.1对各正方形图像块进行二维傅里叶正变换得到对应的正方形图像块的二维频谱：
B＝FFT2(P)
其中，P是正方形图像块，B是正方形图像块的二维频谱，FFT2(·)二维傅里叶正变换。
S3.2对各正方形图像块的二维频谱进行取绝对值、平滑、归一化和以正方形图像块对应的干涉相位图的差异图的均值为幂指数的指数操作得到与其对应的权重图像块。
计算正方形图像块二维频谱的绝对值A：
A＝abs(B)
其中，abs(·)表示取绝对值。然后对正方形图像块二维频谱的绝对值A进行平滑：
S＝smooth(A)
其中，S是平滑后的正方形图像块二维频谱的绝对值，smooth(·)是二维高斯滤波器且其窗口大小根据实际情况确定，本发明实施方式中为7×7像素。接着对S进行归一化：
N＝normalization(S)
其中，normalization(·)是归一化算子。最后，对N取幂指数得到权重图像块W：
W＝Nα
其中：
其中，是正方形图像块对应的干涉相位图的差异图的均值，计算公式为：
其中，mean(·)是取均值算子，DMw是正方形图像块对应的干涉相位图的差异图。
S3.3计算各正方形图像块其对应的权重图像块与正方形图像块二维频谱的Hadamard积(矩阵对应元素相乘)，从而完成干涉相位主要分量提取：
BW＝B·W
BW是权重图像块与正方形图像块二维频谱的Hadamard积。
S4.计算频域滤波后的复干涉相位图；
得到所有正方形图像块二维频谱与其对应权重图像块的Hadamard积后，对所有
Hadamard积进行二维傅里叶逆变换，得到频域滤波后的复干涉相位图像块PW：
PW＝IFFT2(BW)
其中，IFFT2(·)是二维傅里叶逆变换。
正方形图像块之间重叠区域的频域滤波结果来自不同的正方形图像块，因此采用对所有重叠区域的频域滤波结果进行取均值的方法，从而得到频域滤波后的复干涉相位图ψfrequency。
S5.计算空域滤波后的残余复干涉相位图；
S5.1计算残余复干涉相位图ψresidual，残余干涉相位图残余干涉相位图的相位导数方差图VM：
ψresidual＝ψ-ψfrequency
其中，angle(·)是取相位算子；
残余干涉相位图的相位导数方差图VM的计算方法如下：
残余干涉相位图的相位导数方差图VM的大小和残余干涉相位图的大小一样，即为Nx×Ny像素。对于残余干涉相位图中的每一个像素点(s,t)((s,t)∈Nx×Ny)，可以计算该像素点的相位导数方差值VM(s,t)：
其中，和的计算公式为：
在计算残余干涉相位图边缘像素点(第一行、最后一行、第一列和最后一列)的相位导数方差值时，需要先将残余干涉相位图由Nx×Ny像素扩展成(Nx+2)×(Ny+2)像素，即在残余干涉相位图第一列的左边和最后一列的右边各增加一列、在残余干涉相位图第一行的上部和最后一行的下部各增加一行(可以采用镜像、周期延拓等扩展方式，不同扩展方式造成的滤波结果差异可以忽略不计)。
当得到残余干涉相位图中每一个像素点的相位导数方差值后，可以得到残余干涉相位图的相位导数方差图VM：
S5.2计算残余复干涉相位图的实部图像和虚部图像：
ψresidual_real＝real(ψresidual)
ψresidual_imag＝imag(ψresidual)
其中，ψresidual_real是残余复干涉相位图的实部图像，ψresidual_imag是残余复干涉相位图的虚部图像，real(·)是取实部算子，imag(·)是取虚部算子。
S5.3计算空域滤波后的残余复干涉相位图ψspatial的实部图像ψspatial_real和虚部图像ψspatial_imag：
其中，BM3D(·)是BM3D滤波器(Kostadin Dabov,Alessandro Foi,VladimirKatkovnik,and Karen Egiazarian,Image Denoising by Sparse 3-D Transform-DomainCollaborative Filtering,IEEE Transactions on Image Processing,vol.16,no.8,2007)，是残余干涉相位图的相位导数方差图的均值，即
S5.4计算空域滤波后的残余复干涉相位图ψspatial：
S6.计算滤波后的干涉相位图。
S6.1计算滤波后的复干涉相位图ψfiltered：
ψfiltered＝ψspatial+ψfrequency
S6.2计算滤波后的干涉相位图：
图3示出一幅待滤波的干涉相位图。该干涉相位截自美国夏威夷岛上的“大岛”的实测数据，大小为250×3000像素，其中横坐标是列索引，单位是像素；纵坐标是行索引，单位是像素，图3中存在大量的噪声和密集的干涉条纹。图4为采用传统的BM3D方法对图3进行处理后得到的滤波后的干涉相位图。可以看到，图4中图像下半部很多密集干涉条纹的连续性遭到了破坏。
图5为采用本发明提供的基于差异图和补偿滤波的干涉相位滤波方法对图3进行处理后得到的滤波后的干涉相位图。通过图4和图5对比可以看到，本发明得到的滤波后的干涉相位图中，噪声不但得到了有效的滤除，密集干涉条纹也得到了很好的保持。
图6为图3，图4和图5的性能对比图。该性能对比图包含了一个公开通用的性能指标：残差点(D.Ghiglia and M.Pritt.Two-Dimensional Phase Unwrapping:Theory,Algorithms and Software[M].New York,USA:Wiley,1998.)。残差点的计算采用下述方法：对于滤波后的干涉相位图中任意的4个相邻的像素点：我们分别计算两两像素点的缠绕梯度值：
其中，C{·}是缠绕算子，其取值范围为[-π,π)，Δ1，Δ2，Δ3，Δ4是缠绕梯度值。然后对4个缠绕梯度值求和：
Δtotal＝Δ1+Δ2+Δ3+Δ4
如果Δtotal≠0，称像素点(i,j)为残差点，残差点数量越少，滤波后的干涉相位图质量越好。可以看出，相比传统的方法，本发明得到了更少的残差点同时保持了干涉条纹的连续性，更有利于后续的相位解缠。

Claims (6)

1.基于差异图和补偿滤波的干涉相位滤波方法，其特征在于，包括以下步骤：
S1.对于待滤波的干涉相位图，计算干涉相位图的差异图；
其中，i是干涉相位图中任一像素点，j是以像素点i为中心且与像素点i相邻的8个像素点中的任一像素点，Gi是像素点i量化后的像素值；Gj是像素点j量化后的像素值，p(Gi,Gj)是干涉相位图的归一化灰度共生矩阵GLCM(Gi,Gj)；
S2.将干涉相位图转换成复干涉相位图，并对复干涉相位图进行分块；
其中，是干涉相位图，ψ是复干涉相位图；
对复干涉相位图进行分块，将其分成尺寸相同的多个正方形图像块，且列方向的相邻正方形图像块之间有重叠区域，行方向的相邻正方形图像块之间也有重叠区域；各正方形图像块在列方向与其相邻的正方形图像块之间重叠的区域大小与在行方向与其相邻的正方形图像块之间重叠的区域大小相同；
S3.对S2分块得到的各正方形图像块进行频域滤波；
S3.1对各正方形图像块进行二维傅里叶正变换得到对应的正方形图像块的二维频谱；
S3.2对各正方形图像块的二维频谱进行取绝对值、平滑、归一化和以正方形图像块对应的干涉相位图的差异图的均值为幂指数的指数操作得到与其对应的权重图像块；
S3.3计算各正方形图像块对应的权重图像块与正方形图像块二维频谱的Hadamard积，从而完成干涉相位主要分量提取；
S4.计算频域滤波后的复干涉相位图；
得到所有正方形图像块二维频谱与其对应权重图像块的Hadamard积后，对所有Hadamard积进行二维傅里叶逆变换，得到频域滤波后的复干涉相位图像块PW：
PW＝IFFT2(BW)
其中，IFFT2(·)是二维傅里叶逆变换；BW是权重图像块与正方形图像块二维频谱的Hadamard积；
正方形图像块之间重叠区域的频域滤波结果来自不同的正方形图像块，采用对所有重叠区域的频域滤波结果进行取均值的方法，从而得到频域滤波后的复干涉相位图ψfrequency；
S5.计算空域滤波后的残余复干涉相位图；
S5.1计算残余复干涉相位图ψresidual，残余干涉相位图残余干涉相位图的相位导数方差图VM：
S5.2计算残余复干涉相位图的实部图像ψresidual_real和虚部图像ψresidual_imag：
S5.3计算空域滤波后的残余复干涉相位图ψspatial的实部图像ψspatial_real和虚部图像ψspatial_imag：
其中，BM3D(·)是BM3D滤波器，是残余干涉相位图的相位导数方差图的均值，即
S5.4计算空域滤波后的残余复干涉相位图ψspatial：
S6.计算滤波后的干涉相位图；
S6.1计算滤波后的复干涉相位图ψfiltered：
ψfiltered＝ψspatial+ψfrequency
S6.2计算滤波后的干涉相位图：
其中，angle(·)是取相位算子。
2.根据权利要求1所述的基于差异图和补偿滤波的干涉相位滤波方法，其特征在于，S1中，GLCM(Gi,Gj)的计算方法如下：
设待滤波的干涉相位图即输入的干涉相位图的大小为Nx×Ny像素，将干涉相位图像素的像素值[-π,π]量化成13个值，量化采用为间隔的均匀量化，得到量化后的像素值集合G；最后得到一个新的映射关系：f:X×Y→G，其中，X＝{1,2,3,…,Nx}，Y＝{1,2,3,…,Ny}，
分别计算干涉相位图中每一个像素点与其相邻8个像素点在0°和180°方向上的频数P(Gi,Gj,0°和180°)、在45°和225°方向上的频数P(Gi,Gj,45°和225°)、在90°和270°方向上的频数P(Gi,Gj,90°和270°)以及在135°和315°方向上的频数P(Gi,Gj,135°和315°)：
P(Gi,Gj,0°和180°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,k-m＝0,|l-n|＝1}
P(Gi,Gj,45°和225°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,(k-m＝1,l-n＝-1)or(k-m＝-1,l-n＝1)}
P(Gi,Gj,90°和270°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,|k-m|＝1,l-n＝0}
P(Gi,Gj,135°和315°)＝#{[(k,l),(m,n)]∈(Nx,Ny)×(Nx,Ny)|
f(k,l)＝Gi,f(m,n)＝Gj,(k-m＝1,l-n＝1)or(k-m＝-1,l-n＝-1)}
其中，#{·}是计数算子，(k,l)和(m,n)是干涉相位图中像素值分别为Gi和Gj对应的两个像素点的坐标，其中，Gi,Gj∈G；
基于计算得到的频数可以得到下面的矩阵：
其中，θ∈{0°和180°，45°和225°，90°和270°，135°和315°}；
依此能够得到4个矩阵Pθ，将4个矩阵求平均就得到灰度共生矩阵GLCM(Gi,Gj)。
3.根据权利要求1所述的基于差异图和补偿滤波的干涉相位滤波方法，其特征在于，S2中，正方形图像块尺寸为32×32像素。
4.根据权利要求1所述的基于差异图和补偿滤波的干涉相位滤波方法，其特征在于，S3.2中，计算正方形图像块二维频谱的绝对值A：
A＝abs(B)
其中，abs(·)表示取绝对值；B＝FFT2(P)，B是正方形图像块的二维频谱；P是正方形图像块，FFT2(·)二维傅里叶正变换；
对正方形图像块二维频谱的绝对值A进行平滑：
S＝smooth(A)
其中，S是平滑后的正方形图像块二维频谱的绝对值，smooth(·)是二维高斯滤波器；
接着，对S进行归一化：
N＝normalization(S)
其中，normalization(·)是归一化算子；
最后，对N取幂指数得到权重图像块W：
W＝Nα
其中：
其中，是正方形图像块对应的干涉相位图的差异图的均值，计算公式为：
其中，mean(·)是取均值算子，DMw是正方形图像块对应的干涉相位图的差异图。
5.根据权利要求1所述的基于差异图和补偿滤波的干涉相位滤波方法，其特征在于，S5.1中，残余复干涉相位图ψresidual，残余干涉相位图的计算方法如下：
ψresidual＝ψ-ψfrequency
其中，angle(·)是取相位算子。
6.根据权利要求1所述的基于差异图和补偿滤波的干涉相位滤波方法，其特征在于，S5.1中，残余干涉相位图的相位导数方差图VM的计算方法如下：
残余干涉相位图的相位导数方差图VM的大小和残余干涉相位图的大小一样，即为Nx×Ny像素；对于残余干涉相位图中的每一个像素点(s,t),其中(s,t)∈Nx×Ny，可以计算该像素点的相位导数方差值VM(s,t)：
其中，和的计算公式为：
在计算残余干涉相位图边缘像素点即残余干涉相位图第一行、残余干涉相位图最后一行、残余干涉相位图第一列和残余干涉相位图最后一列上各像素点的相位导数方差值时，需要先将残余干涉相位图由Nx×Ny像素扩展成(Nx+2)×(Ny+2)像素，即在残余干涉相位图第一列的左边和最后一列的右边各增加一列、在残余干涉相位图第一行的上部和最后一行的下部各增加一行；
当得到残余干涉相位图中每一个像素点的相位导数方差值后，可以得到残余干涉相位图的相位导数方差图VM：
CN201711304057.0A 2017-12-11 2017-12-11 基于差异图和补偿滤波的干涉相位滤波方法 Active CN108090878B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711304057.0A CN108090878B (zh) 2017-12-11 2017-12-11 基于差异图和补偿滤波的干涉相位滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711304057.0A CN108090878B (zh) 2017-12-11 2017-12-11 基于差异图和补偿滤波的干涉相位滤波方法

Publications (2)

Publication Number Publication Date
CN108090878A CN108090878A (zh) 2018-05-29
CN108090878B true CN108090878B (zh) 2018-11-09

Family

ID=62174556

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711304057.0A Active CN108090878B (zh) 2017-12-11 2017-12-11 基于差异图和补偿滤波的干涉相位滤波方法

Country Status (1)

Country Link
CN (1) CN108090878B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110060220A (zh) * 2019-04-26 2019-07-26 中国科学院长春光学精密机械与物理研究所 基于改进bm3d算法的图像去噪方法及系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103208101A (zh) * 2013-03-28 2013-07-17 中国科学院对地观测与数字地球科学中心 一种基于局部信噪比的干涉图滤波方法
CN103871030A (zh) * 2014-02-17 2014-06-18 中国科学院电子学研究所 一种干涉图像的滤波方法及设备
CN105469368A (zh) * 2015-11-30 2016-04-06 中国人民解放军国防科学技术大学 一种干涉相位滤波方法
CN106485677A (zh) * 2016-09-30 2017-03-08 湖南鼎方电子科技有限公司 一种快速高效干涉相位滤波方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN100493444C (zh) * 2007-12-12 2009-06-03 中国科学院上海光学精密机械研究所 高分辨率光学相干层析成像方法
CN101881831B (zh) * 2010-06-24 2012-07-18 中国人民解放军信息工程大学 基于差分滤波的多波段InSAR相位解缠方法
WO2014003598A1 (en) * 2012-06-28 2014-01-03 Intel Corporation Inter-carrier interference phase noise compensation

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103208101A (zh) * 2013-03-28 2013-07-17 中国科学院对地观测与数字地球科学中心 一种基于局部信噪比的干涉图滤波方法
CN103871030A (zh) * 2014-02-17 2014-06-18 中国科学院电子学研究所 一种干涉图像的滤波方法及设备
CN105469368A (zh) * 2015-11-30 2016-04-06 中国人民解放军国防科学技术大学 一种干涉相位滤波方法
CN106485677A (zh) * 2016-09-30 2017-03-08 湖南鼎方电子科技有限公司 一种快速高效干涉相位滤波方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Interference imaging with a spatial spiral phase filter;Steverin furhapter et al.;《Proceedings of SPIE-The International Society for Optical Engineering》;20050930;第1-5页 *
基于Gabor滤波器的干涉相位图滤波算法;贾凌春 等;《四川兵工学报》;20131231;第34卷(第12期);第123-127页 *

Also Published As

Publication number Publication date
CN108090878A (zh) 2018-05-29

Similar Documents

Publication Publication Date Title
Mäkinen et al. Exact transform-domain noise variance for collaborative filtering of stationary correlated noise
WO2016086699A1 (zh) 一种结合局部频率估计的小波域InSAR干涉相位滤波方法
JP2014154141A5 (zh)
Singh et al. Comparative performance analysis of various wavelet and nonlocal means based approaches for image denoising
CN103873855B (zh) 一种基于人类视觉基本特性的立体图像客观质量评价方法
CN102289670A (zh) 一种具有光照鲁棒性的图像特征提取方法
CN108090878B (zh) 基于差异图和补偿滤波的干涉相位滤波方法
CN103903228A (zh) 一种基于hwd变换的非局部图像去噪方法
CN102196155B (zh) 基于Surfacelet变换的系数自适应收缩视频去噪方法
Kadiri et al. Magnitude-phase of the dual-tree quaternionic wavelet transform for multispectral satellite image denoising
Wu et al. Attenuating seismic noise via incoherent dictionary learning
CN104240201A (zh) 基于群-轮廓小波变换断口图像去噪和增强方法
CN107146206B (zh) 基于四维块匹配滤波的高光谱遥感图像去噪方法
CN106296591A (zh) 基于马氏距离的非局部均匀数字图像去噪方法
Bala et al. Image denoising method using curvelet transform and wiener filter
Thai et al. Performance evaluation of high dynamic range image tone mapping operators based on separable non-linear multiresolution families
Chang et al. New interpolation algorithm for image inpainting
CN106485677B (zh) 一种快速高效干涉相位滤波方法
CN103747268B (zh) 一种分层自适应阈值视频去噪方法
CN108198140A (zh) 基于ncsr模型的三维协同滤波去噪方法
CN107371013B (zh) 基于色相角和差通道的立体图像质量评价方法
CN105894565A (zh) 一种基于cuda的无镜立体图像重建方法
Fang et al. A novel InSAR phase denoising method via nonlocal wavelet shrinkage
Oleiwi et al. High frequency coefficient effect on image based on contourlet transformation
Nayak et al. Object tracking in curvelet domain with dominant curvelet subbands

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
GR01 Patent grant
GR01 Patent grant