CN103823219B - 自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法 - Google Patents

Abstract

本发明公开了一种自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，该方法包括将干涉相位变换到复数域并对复数域的实部和虚部做非局部滤波，对得到的复干涉相位的实部和虚部计算，得到干涉相位，选出残差点最少的滤波后的干涉相位并记录迭代次数参数m＝1；增长搜索窗口、邻域窗口尺寸，m＝m+1，将残差点最少的第m-1次滤波后干涉相位作为待滤波图像进行非局部滤波，选出残差点最少的第m次滤波后干涉相位，获取最佳邻域窗口值；判断是否继续迭代，若继续迭代，则增长搜索窗口、邻域窗口尺寸；若停止迭代，则输出滤波后的干涉相位。

Description

自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法

技术领域

本发明涉及干涉合成孔径雷达信号处理领域，尤其涉及一种自适应迭代的非局部干涉相位滤波方法。

背景技术

干涉合成孔径雷达(InterferometricSyntheticApertureRadar，InSAR)利用两个通道雷达回波的干涉相位信息提取地形高程、运动目标速度、地表形变等各种信息，将合成孔径雷达(SAR)的测量拓展到三维空间。其具有全天时、全天候、高精度的特点，因此在地形测绘、冰川研究、海洋测绘以及地面沉降监测等多个领域都有广泛的应用。

干涉SAR在获取相位信息的过程中，干涉相位会受到热噪声去相干、基线或几何去相干、时间去相干等随机误差的影响，具体表现为相位噪声。低质量的干涉条纹将会影响后续的干涉相位展开及DEM反演的准确性，

最终给高程测量带来误差，影响测绘精度。因此，为获取高质量的干涉相位，必须对干涉条纹进行滤波处理。

针对该问题，前人提出了多种相位滤波方法，这些算法可以分为两类：空域滤波算法，如圆周期均值滤波(EichelPH，GhigliaDC，JakowatzCVJr，etal..SpotlightSARinterferometryforterrainelevationmappingandinterferometricchangedetection.SandiaNationalLabsTechReport，WashingdonDC，1993，US-DOE.)、Lee滤波(LeeJS，PapathanassiouP，AinsworthTL，etal..AnewtechniquefornoisefilteringofSARinterferometricphaseimages.IEEETrans.onGeoscienceandRemoteSensing，1998，36(5)：1456-1465.)等，以及变换域滤波算法，如Goldstein滤波(GoldsteinRMandWernerCL.Radarinterferogramfilteringforgeophysicalapplications.GeophysicalResearchLetters，1998，25(21)：4035-4038.)、小波变换滤波(MartinezCLandFabregasX.ModelingandreductionofSARmterferometricphasenoiseinthewaveletdomain.IEEETrans.onGeoscienceandRemoteSensing，2002，40(12)：2553-2566.)等。其中，变换域滤波算法在实际中应用更为广泛，其基于信号与噪声在变换域中所处的位置不同这一假设来进行噪声抑制，在很多情况下可以得到满意结果。然而，在相干性较差或地形变化剧烈的情况下，由于细节信息与噪声信息在变换域中所处位置无法完全区分，因此上述算法不能同时满足去除噪声和保持相位细节两方面的要求。为得到高质量的干涉相位，有必要研究一种适用于低相干及地形变化剧烈区域的相位滤波方法。

发明内容

(一)要解决的技术问题

为解决上述问题，本发明提供了一种自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，以克服传统的基于像素的处理方法不能兼顾去噪和细节保持的缺点，从而提高干涉相位的精度。

(二)技术方案

根据本发明的一个方面，提供一种自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，包括如下步骤：

步骤S1：将干涉相位φ变换到复数域；

步骤S2：设定初始搜索窗口尺寸参数S，初始邻域窗口参数N，对复数域的实部进行非局部滤波处理，得到多个滤波后复干涉相位的实部

步骤S3：对复数域的虚部进行非局部滤波处理，得到多个滤波后的复干涉相位的虚部

步骤S4：对多个滤波后的复干涉相位的实部和虚部进行计算，得到多个滤波后的干涉相位为从多个滤波的干涉相位的结果中选出残差点最少的滤波后的干涉相位，同时得到该残差点最少的第1次去除噪声的滤波后的干涉相位对应的邻域窗口尺寸参数值n，并记录迭代次数参数m＝1；

步骤S5：设定搜索窗口尺寸参数增长步长为r，邻域窗口尺寸增长步长为t，令S＝S+r，N＝(n，n+t，n+2t)，m＝m+1，将得到的残差点最少的第m-1次去除噪声的滤波后的干涉相位作为待滤波的图像，计算待滤波的图像的残差点数目RES1，并对待滤波的图像进行非局部滤波处理，并选出残差点最少的第m次去除噪声的滤波后的干涉相位，确定残差点最少的第m次去除噪声的滤波后的干涉相位对应的邻域窗口参数值N0；

步骤S6：考察邻域窗口参数值N＝n+k大于残差点最少的第m次去除噪声的滤波后的干涉相位所对应的邻域窗口参数值N0＝n+2t，k＞2t且为偶数时的第m次去除噪声的滤波后的干涉相位，获取最佳邻域窗口值；如果残差点最少的第m次去除噪声的滤波后的干涉相位所对应的邻域窗口参数值N0≠n+2t，则记录该残差点最少的第m次去除噪声的滤波后的干涉相位的残差点数目RES2，本次迭代结束，转到步骤S7；

步骤S7：判断是否继续迭代，若继续迭代，则转到步骤S5；若停止迭代，则输出滤波后的干涉相位。

(三)有益效果

本发明基于非局部去噪的思想，利用图像的冗余信息进行滤波，根据干涉相位的自身特性，利用干涉相位噪声标准差对权值计算中的衰减系数进行估计，实现了衰减系数的自适应，之后通过调节窗口尺寸及自动迭代，在滤除噪声的同时能够更好的利用相位的纹理信息，使得干涉相位细节结构不被破坏，从而克服了传统的基于像素的滤波方法不能兼顾去噪和细节保持的缺点，实现高精度的干涉相位滤波，为高精度的干涉测量提供了条件，为干涉合成孔径雷达干涉相位滤波提供了一种有效的方法。

附图说明

图1为本发明实施例干涉合成孔径雷达自适应迭代的非局部相位滤波的流程图。

图2为本发明干涉合成孔径雷达自适应迭代的非局部相位滤波的实施例流程图。

图3为意大利Etna火山的实测干涉合成孔径雷达干涉相位图。

图4a-图4d为几种滤波方法对图3的干涉相位进行滤波后的结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。且在附图中，以简化或是方便标示。再者，附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。

如图1示出本发明提供的一种自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法包括如下步骤：

步骤S1：将干涉相位φ变换到复数域；

步骤S2：设定初始搜索窗口尺寸参数S，初始邻域窗口参数N，对复数域的实部进行非局部滤波处理，得到多个滤波后复干涉相位的实部

步骤S3：对复数域的虚部进行非局部滤波处理，得到多个滤波后的复干涉相位的虚部

步骤S4：对多个滤波后的复干涉相位的实部和虚部进行计算，得到多个滤波后的干涉相位为从多个滤波的干涉相位的结果中选出残差点最少的滤波后的干涉相位，同时得到该残差点最少的第1次去除噪声的滤波后的干涉相位对应的邻域窗口尺寸参数值n，并记录迭代次数参数m＝1；

步骤S5：设定搜索窗口尺寸参数增长步长为r，邻域窗口尺寸增长步长为t，令S＝S+r，N＝(n，n+t，n+2t)，m＝m+1，将得到的残差点最少的第m-1次去除噪声的滤波后的干涉相位作为待滤波的图像，计算待滤波的图像的残差点数目RES1，并对待滤波的图像进行非局部滤波处理，并选出残差点最少的第m次去除噪声的滤波后的干涉相位，确定残差点最少的第m次去除噪声的滤波后的干涉相位对应的邻域窗口参数值N0；

步骤S6：考察邻域窗口参数值N＝n+k大于残差点最少的第m次去除噪声的滤波后的干涉相位所对应的邻域窗口参数值N0＝n+2t，k＞2t且为偶数时的第m次去除噪声的滤波后的干涉相位，获取最佳邻域窗口值；如果残差点最少的第m次去除噪声的滤波后的干涉相位所对应的邻域窗口参数值N0≠n+2t，则记录该残差点最少的第m次去除噪声的滤波后的干涉相位的残差点数目RES2，本次迭代结束，转到步骤S7；

步骤S7：判断是否继续迭代，若继续迭代，则转到步骤S5；若停止迭代，则输出滤波后的干涉相位。

图2为本发明自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法的实施例流程图。本发明中的参数选取可选择但不限于该实施例中的参数选取情况。如图2所示，本实施例包括以下步骤：

步骤S1：将干涉相位φ变换到复数域e^jφ＝cosφ+jsinφ中，分别取复数域e^jφ的实部Re{e^jφ}＝cosφ、虚部Im{e^jφ}＝sinφ；j为虚数，φ为滤波前的干涉相位；

步骤S2：设定窗口尺寸参数(S，N)＝(1，1)，(1，2)，(1，3)，其中，所述搜索窗口尺寸参数S对应搜索窗口大小(2S+1)×(2S+1)，且参数S的初始值为小于4的自然数；所述邻域窗口尺寸参数N对应邻域窗口大小(2N+1)×(2N+1)，且参数N取值为多个逐渐增大的自然数。

对复数域的实部进行非局部滤波处理，得到三个滤波后复干涉相位的实部

上述对复数域的实部进行非局部滤波处理，得到滤波后的坐标i处的像素值NL(υ(i))表示为：

$\begin{array}{c}\mathrm{NL}\left(\mathrm{\υ}\left(i\right)\right)={\mathrm{\Σ}}_j\mathrm{\ω}(i,j)\mathrm{\υ}\left(j\right)={\mathrm{\Σ}}_j\mathrm{\υ}\left(j\right)\·\exp (-{\left|\right|N_i-N_j\left|\right|}^2/h^2)/Z_i,j\∈S_i\\ Z_i={\mathrm{\Σ}}_j\exp (-{\left|\right|N_i-N_j\left|\right|}^2/h^2)\end{array}$

其中，NL表示非局部(Non-local)的缩写，υ(i)表示图像坐标i处的像素值，S_i表示以像素υ(i)为中心的搜索窗口，j表示搜索窗口S_i内的各个像素的坐标，ω(i，j)表示为恢复像素值υ(i)而对像素υ(j)赋予的权值，N_i与N_j分别表示以像素υ(i)和υ(j)为中心的邻域像素块，zi为归一化常数，用以保证h表示非局部滤波处理中的衰减参数，非局部滤波处理中的衰减系数h的值是根据待滤波的干涉相位的相位噪声及邻域窗口大小设置：在条纹较稀疏即干涉相位每变化2π所经历的像素在30个像素以上的情况下，衰减系数h值在5到几百的范围内得到的滤波结果变化不大，衰减系数h值取值为经典值10σ_n；而当条纹较密集即干涉相位每变化2π所经历的像素在30个像素以下时，最优衰减系数h值与噪声标准差及邻域窗口大小近似呈线性，且衰减系数h值在最优值附近一个±1到±5的范围内得到的滤波结果变化不大，即与理想干涉相位的均方误差值变化在0.1以内，即具有一定的鲁棒性，为采用拟合近似提供了条件，通过多项式拟合，最优衰减系数h值表示为：h＝1.3σ_n+0.5·(2N+1)+c₃，其中，σ_n表示噪声标准差，c₃为常数项，c₃取值在-0.5到-1之间；

步骤S3：对复数域的虚部进行步骤S2的非局部滤波处理，得到三个滤波后的复干涉相位的虚部

步骤S4：对滤波后的复干涉相位的实部和虚部进行计算，得到三个滤波后的干涉相位为 为实部，为虚部，从三个滤波的干涉相位的结果中选出残差点最少的滤波后的干涉相位；

步骤S5：设定窗口尺寸参数(S，N)＝(2，1)，(2，3)，(2，5)，将得到的残差点最少的滤波后的干涉相位作为待滤波的图像，对待滤波的图像进行与步骤S1-步骤S3相同的非局部滤波处理，并选出残差点最少的滤波结果，同时得到该残差点最少的滤波结果对应的窗口尺寸参数(S，N)＝(2，n)，n为邻域窗口尺寸参数值；

步骤S6：本次迭代中采用的搜索窗口尺寸参数值S＝S0≥3时，将上一次迭代(S0＝3时对应步骤S5，S0＞3时对应步骤S7)得到的残差点最少的滤波结果作为待滤波的图像，计算待滤波的图像的残差点数目RES1，并分别计算(S，N)＝(S()，n)，(S0，n+2)，(S0，n+4)时的滤波结果，确定本次残差点最少的滤波结果邻域窗口尺寸对应的参数N设为N0；

步骤S7：如果本次迭代中残差点最少的滤波结果所对应的邻域窗口参数值N0≠n+4，记录该滤波结果的残差点数目RES2，本次迭代结束，转到步骤S8，若N0＝n+4，则继续考察N＞N0时的滤波结果，用以获取最佳邻域窗口值；

当N0＝n+4时，继续考察N＞N0时的滤波结果以获取最佳邻域窗口值的具体操作为：计算S0，n+k，k＞4且为偶数时的滤波结果，其中，k表示邻域窗口边长增加了2k个像素，更新N0值，直至满足N0≠n+k(即残差点最少的滤波结果对应的邻域窗口大小不是目前的最大值)或者N＝n+k，n+(k-2)，n+(k-4)时的残差点数目不变(即随邻域窗口增大，滤波效果基本不变)这两者中的任一条件时，本次迭代结束，转到步骤S8；

步骤S8：判断是否继续迭代，若4≤S＜10且本次迭代的待滤波图像与滤波结果的残差点数目RES1，RES2满足：(RES1-RES2)／RES1＜C，其中，C为常数，且0＜C＜0.5，停止迭代，得到滤波后的干涉相位；或当S≥10时，停止迭代，得到滤波后的干涉相位；否则，S0＝S0+1，n＝N0，转到步骤S6；若4≤S＜10且本次迭代的待滤波图像与滤波结果的残差点数目RES1，RES2满足：(RES1-RES2)／RES1＞C，则S0＝S0+1，n＝N0，转到步骤S6；若S≥10，停止迭代，得到滤波后的干涉相位。

下面通过实测数据验证了本发明方法的有效性。在本实例中实施本发明方法所采用的硬件设施为32-bit3.10-GHzCore^TMi5-2400处理器，软件设施为WindowsXPProfessional操作系统及MatlabR2010a。图3为意大利Etna火山的干涉合成孔径雷达干涉相位图。图4a-图4d为几种滤波方法对图3的干涉相位进行滤波后的结果，分别利用圆周期均值滤波的结果，Goldstein滤波的结果，DWT滤波的结果以及本发明方法滤波后的干涉相位图。其中图4a为圆周期均值滤波结果，滤波窗口为7×7，图4b为Goldstein滤波结果，分块大小取8×8，图4c为DWT滤波结果，小波变换层数为2，图4d为本发明方法滤波结果。可以看出，圆周期均值滤波后在条纹密集处有一定的交叠或间断，破坏了相位细节，Goldstein滤波去噪效果很差，在大部分区域改善效果很小，且在干涉相位左侧部分有较明显的结构破坏。DWT滤波后仍存在较明显的噪声，且在条纹密集处也存在交叠或间断现象。本发明方法滤波后的干涉相位最为平滑，且基本恢复了该区域的干涉条纹，且条纹密集处结构保持良好。滤波后干涉相位图中的残差点数目可用于评价去噪效果的好坏，计算出滤波前后的残差点数目如表1。从表中同样可以看出本发明方法滤波后剩余的残差点数目最少，去噪效果最好。

表1实测干涉相位滤波前后残差点数目比较

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，包括如下步骤：

步骤S1：将干涉相位φ变换到复数域；

步骤S2：设定初始搜索窗口尺寸参数S，初始邻域窗口参数N，对复数域的实部进行非局部滤波处理，得到多个滤波后复干涉相位的实部

步骤S3：对复数域的虚部进行非局部滤波处理，得到多个滤波后的复干涉相位的虚部

步骤S4：对多个滤波后的复干涉相位的实部和虚部进行计算，得到多个滤波后的干涉相位为从多个滤波的干涉相位的结果中选出残差点最少的滤波后的干涉相位，同时得到该残差点最少的第1次去除噪声的滤波后的干涉相位对应的邻域窗口尺寸参数值n，并记录迭代次数参数m＝1；

步骤S5：设定搜索窗口尺寸参数增长步长为r，邻域窗口尺寸增长步长为t，令S＝S+r，N＝(n，n+t，n+2t)，m＝m+1，将得到的残差点最少的第m-1次去除噪声的滤波后的干涉相位作为待滤波的图像，计算待滤波的图像的残差点数目RES1，并对待滤波的图像进行非局部滤波处理，并选出残差点最少的第m次去除噪声的滤波后的干涉相位，确定残差点最少的第m次去除噪声的滤波后的干涉相位对应的邻域窗口参数值N0；

步骤S6：考察邻域窗口参数值N＝n+k大于步骤S5得到的残差点最少的第m次去除噪声的滤波后的干涉相位所对应的邻域窗口参数值N0＝n+2t，k＞2t且为偶数时的第m次去除噪声的滤波后的干涉相位，获取最佳邻域窗口值；如果步骤S5得到的残差点最少的第m次去除噪声的滤波后的干涉相位所对应的邻域窗口参数值N0≠n+2t，则记录该残差点最少的第m次去除噪声的滤波后的干涉相位的残差点数目RES2，本次迭代结束，转到步骤S7；

步骤S7：判断是否继续迭代，若继续迭代，则转到步骤S5；若停止迭代，则输出滤波后的干涉相位。

2.根据权利要求1所述的自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，其特征在于，所述复数域e^jφ＝cosφ+jsinφ中，分别取复数域e^jφ的实部Re{e^jφ}＝cosφ、虚部Im{e^jφ}＝sinφ；j为虚数，φ为滤波前的干涉相位。

3.根据权利要求1所述的自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，其特征在于，所述搜索窗口尺寸参数S，其对应的搜索窗口大小为(2S+1)×(2S+1)，且参数S的初始值为小于4的自然数；邻域窗口尺寸参数N，其对应邻域窗口大小(2N+1)×(2N+1)，且参数N取值为多个逐渐增大的自然数。

4.根据权利要求1所述的自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，其特征在于，所述对复数域的实部进行非局部滤波处理，得到滤波后的坐标i处的像素值NL(v(i))表示为：

NL(v(i))＝∑_jw(i，j)v(j)＝∑_jv(j)·exp(-||N_i-N_j||²/h²)/Z_i，j∈S_i

Z_i＝∑_jexp(-||N_i-N_j||²/h²)

其中，NL表示非局部，v(i)表示图像坐标i处的像素值，S_i表示以像素v(i)为中心的搜索窗口，j表示搜索窗口S_i内的各个像素的坐标，w(i，j)表示为恢复像素值v(i)而对像素v(j)赋予的权值，N_i与N_j分别表示以像素v(i)和v(j)为中心的邻域像素块，Z_i为归一化常数，用以保证h为非局部滤波处理中的衰减系数。

5.根据权利要求4所述的自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，其特征在于，根据待滤波的干涉相位的相位噪声及邻域窗口大小设置非局部滤波处理中的衰减系数h的值的步骤为：

在条纹较稀疏即干涉相位每变化2π所经历的像素在30个像素以上的情况下，衰减系数h值在5到几百的范围内得到的滤波结果变化不大，衰减系数h值取值为经典值10σ_n；而当条纹较密集即干涉相位每变化2π所经历的像素在30个像素以下时，衰减系数h值的最优值与噪声标准差及邻域窗口大小近似呈线性，且衰减系数h值在最优值附近一个±1到±5的范围内得到的滤波结果变化不大，即与理想干涉相位的均方误差值变化在0.1以内，即具有一定的鲁棒性，为采用拟合近似提供了条件，通过多项式拟合，最优衰减系数h值表示为h＝1.3σ_n+0.5·(2N+1)+c₃，其中，σ_n表示噪声标准差，c₃为常数项，c₃取值在-0.5到-1之间。

6.根据权利要求1所述的自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，其特征在于，所述当N0＝n+2t时，继续考察N＞N0时的第m次去除噪声的滤波后的干涉相位以获取最佳邻域窗口值的具体操作为：计算N＝n+k，k＞2t且为偶数时的第m次去除噪声的滤波后的干涉相位，其中，k表示邻域窗口边长增加了2k个像素，更新N0值及残差点最少的第m次去除噪声的滤波后的干涉相位，直至残差点最少的第m次去除噪声的滤波后的干涉相位对应的邻域窗口大小不是目前的最大值即满足N0≠n+k时，或者随邻域窗口增大，滤波效果基本不变即满足N＝n+k，n+(k-2)，n+(k-4)时的残差点数目不变时，记录该N0值对应的第m次去除噪声的滤波后的干涉相位的残差点数目RES2，本次迭代结束，转到步骤S7。

7.根据权利要求1所述的自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法，其特征在于，所述判断是否继续迭代的具体操作为：若4≤S＜10且待滤波图像与残差点最少的第m次去除噪声的滤波后的干涉相位的残差点数目RES1，RES2满足：(RES1-RES2)/RES1＜C，其中，C为常数，且0＜C＜0.5，停止迭代，输出残差点最少的第m次去除噪声的滤波后的干涉相位作为最终结果；若4≤S＜10且待滤波图像与残差点最少的第m次去除噪声的滤波后的干涉相位的残差点数目RES1，RES2满足：(RES1-RES2)/RES1＞C，则继续迭代，n＝N0，转到步骤S5；若S≥10，则停止迭代，输出残差点最少的第m次去除噪声的滤波后的干涉相位作为最终结果。