CN103208101A - 一种基于局部信噪比的干涉图滤波方法 - Google Patents

Abstract

本发明涉及一种基于局部信噪比的干涉图滤波装置，所述方法通过窗口在所述干涉图中的滑动，依次对滑动窗口内的干涉图进行滤波，其中，在所述滑动窗口内对干涉图进行滤波包括：对所述滑动窗口内的干涉图的干涉相位进行分解，得到所述干涉图在所述滑动窗口内的估计频率分量；根据所述干涉相位和所述估计频率分量计算出所述干涉图在所述滑动窗口内的局部信噪比；根据所述局部信噪比计算出滤波参数；以及根据所述滤波参数对所述滑动窗口内的干涉图进行滤波。本发明能准确估计干涉图中局部噪声的程度，在提供传统滤波方法降噪能力的同时，也降低了滤波对干涉图相位信息造成的损失，体现了一种整体优越性。

Description

一种基于局部信噪比的干涉图滤波方法

技术领域

本发明涉及图像处理领域，尤其涉及一种基于局部信噪比的InSAR干涉图的滤波方法。

背景技术

合成孔径雷达干涉测量技术（Interferometric Synthetic ApertureRadar,INSAR）是以同一地区的两张SAR图像为基本处理数据，通过求取两幅SAR图像的相位差，获取干涉图像，然后经相位解缠，从干涉条纹中获取地形高程数据的空间对地观测新技术。

在合成孔径雷达干涉测量（InSAR）的数据处理流程中，对干涉相位图进行滤波是一个非常重要的步骤。如上所述，干涉相位图用于后续的相位解缠，干涉相位图的质量会影响后续相位解缠的复杂程度，决定了利用InSAR技术生成数字高程模型以及地表形变监测的精度，对数字高程产品以及地表形变监测的应用都有影响。由于干涉相位图通常是条纹状的，通常也将其称之为干涉条纹图。请注意，下文中的干涉图、干涉相位图、干涉条纹图都意指同一事物。综上所述，必须采用有效的滤波方法对干涉图进行滤波降噪处理。

目前，较为常用的干涉相位图的滤波方法是Goldstein频域滤波算法，它是由Goldstein和Werner在1998年提出的。该方法针对干涉图中的噪声往往是高频分量，而干涉条纹的空间频率则相对稳定这一特点，基于噪声点和一般信息点的频率特性的不同，将滤波处理从常用的空间域转换到频率域，分析处理噪声显著的频谱特性，从而实现对干涉相位图的频域滤波。

该方法的缺陷是滤波因子需要人为设定，而滤波因子的大小直接影响了滤波程度的大小；进一步，由于干涉图在不同的区域具有不同的噪声,因此通常情况下很难确定其滤波因子。具体而言，该方法中的滤波参数α是通过经验人为设定的,取值范围例如在0到1之间。实验结果表明，滤波效果随α的增大而增强，但是带来的负面影响是使图像的分辨率降低，图像的一些区域将变得模糊。所以，α的取值应该适当，在达到较好的降噪性能的同时要保证图像仍然具有较高的分辨率，使两者之间达到平衡。然而，通常情况下干涉图局部窗口内的噪声程度是未知的，因此如何选取准确的滤波参数α提高频域滤波的自适应性和准确度是一个值得研究的问题。这也是目前Goldstein频域滤波算法最大的缺陷。

所以，如何准确估计InSAR干涉相位图中的噪音程度非常重要，它是确定滤波尺度的重要指标，直接影响了干涉相位的滤波效果。

发明内容

本发明的目的是提供一种能够克服上述缺陷的基于局部信噪比的干涉图滤波方法。

本发明提供了一种基于局部信噪比的干涉图滤波装置，所述方法通过窗口在所述干涉图中的滑动，依次对滑动窗口内的干涉图进行滤波，其中，在所述滑动窗口内对干涉图进行滤波包括：对所述滑动窗口内的干涉图的干涉相位进行分解，得到所述干涉图在所述滑动窗口内的估计频率分量；根据所述干涉相位和所述估计频率分量计算出所述干涉图在所述滑动窗口内的局部信噪比；根据所述局部信噪比计算出滤波参数；以及根据所述滤波参数对所述滑动窗口内的干涉图进行滤波。

优选地，所述对所述滑动窗口内的干涉图的干涉相位进行分解，得到所述干涉图在所述滑动窗口内的估计频率分量的步骤包括：识别所述干涉相位的极大值点和极小值点；根据所述极大值点和所述极小值点，计算极大值包络面和极小值包络面；计算所述极大值包络面和极小值包络面的平均值；以及根据所述滑动窗口内的干涉相位和所述平均值计算所述干涉图在所述滑动窗口内的估计频率分量。

优选地，所述根据所述滑动窗口内的干涉相位和所述平均值计算所述干涉图在所述滑动窗口内的估计频率分量的步骤包括：从所述干涉相位中减去所述平均值，得到所述干涉图在所述滑动窗口内的估计频率分量；或者判断所述频率分量是否满足预定条件；如果不满足，则用所述频率分量代替所述干涉相位，从所述识别所述滑动窗口内的干涉相位的极大值点和极小值点的步骤开始迭代，得出满足所述预定条件的第一频率分量，将所述第一频率分量作为所述估计频率分量。

优选地，所述对所述滑动窗口内的干涉图的干涉相位进行分解，得到所述干涉图在所述滑动窗口内的估计频率分量的步骤还包括：根据所述干涉相位和所述估计频率分量得到所述干涉相位的余量；用所述余量代替所述干涉相位，从所述识别所述滑动窗口内的干涉相位的极大值点和极小值点的步骤开始预定次数的迭代，得出满足所述预定条件的一个或多个第二频率分量，将所述一个或多个第二频率分量作为所述估计频率分量。

优选地，所述局部信噪比是根据所述干涉相位和所述估计频率分量的关系计算得出的，所述局部信噪比是伪信噪比。

优选地，所述滤波参数是对所述局部信噪比进行归一化处理得到的。

优选地，所述根据所述滤波参数对所述滑动窗口内的干涉图进行滤波的步骤包括：将所述滑动窗口内的干涉图从空间域变换到频域，得到所述干涉图的频谱；用所述滤波参数对所述频谱进行频域平滑滤波；以及将所述频域平滑滤波后的频谱变换到空间域。

优选地，所述极大值点和极小值点是通过像素点邻域搜索来识别的，所述极大值点和极小值点分别有一个或多个。

优选地，所述极大值包络面和所述极小值包络面是通过Delaunay三角剖分插值得到的。

优选地，所述干涉图是通过合成孔径雷达干涉测量InSAR得到的。

本发明通过对滑动窗口中的干涉图进行相位分解，利用分解结果得到滑动窗口中的局部信噪比，根据局部信噪比构造针对每个滑动窗口的滤波参数，自适应地对滑动窗口内的干涉图进行滤波。本发明能准确估计干涉图中局部噪声的程度，在提供传统滤波方法降噪能力的同时，也降低了滤波对干涉图相位信息造成的损失，体现了一种整体优越性。

附图说明

图1是根据本发明实施例的基于局部信噪比的干涉图滤波方法的流程图。

具体实施方式

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

鉴于此问题，本发明的干涉图滤波方法利用经验模态分解（EMD）技术分解干涉相位得到干涉图的局部信噪比，通过该局部信噪比确定噪声程度，进而对干涉图进行对噪声程度有针对性的滤波。

举例说明，根据本发明的方法可以采用二维EMD技术对干涉相位进行分解，得到相位图在不同频率成分的分量，其中高频分量主要包含了原始相位的噪声信息，低频分量包含了原始相位大尺度的变化趋势。通过获得的不同频率成分的相位分量，计算原始干涉图的局部信噪比。用局部信噪比作为衡量噪声程度的参数，把它设为滤波因子。利用上一步计算得到的滤波因子对原始干涉相位图进行频域加权滤波。

在此，为了更好地理解本发明，首先对经验模态分解进行介绍。

由于大多数所有要分析的数据都不是本征模函数，在任意时间点上，数据可能包含多个波动模式，这就是简单的希尔伯特变换不能完全表征一般数据的频率特性的原因。于是需要对原数据进行EMD分解来获得本征模函数。

EMD分解方法是基于以下假设条件：⑴数据至少有两个极值，一个最大值和一个最小值；⑵数据的局部时域特性是由极值点间的时间尺度唯一确定；⑶如果数据没有极值点但有拐点，则可以通过对数据微分一次或多次求得极值，然后再通过积分来获得分解结果。这种方法的本质是通过数据的特征时间尺度来获得本征波动模式，然后分解数据。这种分解过程可以形象地称之为“筛选（sifting）”过程。

分解过程是：找出原数据序列X(t）所有的极大值点并用三次样条插值函数拟合形成原数据的上包络线；同样，找出所有的极小值点，并将所有的极小值点通过三次样条插值函数拟合形成数据的下包络线，上包络线和下包络线的均值记作ml，将原数据序列X(t）减去该平均包络ml，得到一个新的数据序列h，X(t)-ml=hl

由原数据减去包络平均后的新数据，若还存在负的局部极大值和正的局部极小值，说明这还不是一个本征模函数，需要继续进行“筛选”。

图1是根据本发明实施例的基于局部信噪比的干涉图滤波方法的流程图。

在步骤101，识别干涉图在滑动窗口内的干涉相位的极大值点和极小值点。

滤波方式可以采用滑动窗口的方式进行，在选定滑动窗口的大小后，滑动窗口在干涉图中沿行方向自左向右滑动，到达边缘后跳至下一行，从而依次对滑动窗口内的干涉相位进行滤波。滑动窗口的大小可以设定为32×32，窗口滑动的重叠覆盖为16。另外，本征模函数估计值F_estimation被初始化，令F_estimation=I。

此步骤旨在识别滑动窗口内的干涉相位I的局部极大值点和局部极小值点，一个滑动窗口内通常存在多个局部极大值点和多个局部极小值点。局部极值的选取可以采用像素点邻域搜索。例如，在像素点的8邻域内进行搜索，比较当前像素点和最邻近的8个像素点的大小。如果当前像素点的相位值大于周围8个邻近像素点的相位值，则判定当前像素点为极大值点；同理，如果当前像素点的相位值小于周围8个邻近像素点的相位值，则判定当前像素点为极小值点。这样，可以找出多个极大值点和多个极小值点。

在步骤102，计算极大值包络面（UE）和极小值包络面（LE）。

找出极值点之后，对全部的极大值点和全部的极小值点分别进行曲面拟合，计算出极大值包络面（UE）和极小值包络面（LE）。例如，可以采用Delaunay三角剖分插值的方法得到极大值点包络曲面UE和极小值点包络曲面LE。由于Delaunay三角剖分插值的运算效率较高，在插值点较多时的运算速度相对较快，出于综合考虑不失为一种优选的插值方法。然而，也可以采用其他具有相同效果的插值方法。

在步骤103，计算均值包络面。

将极大值包络面UE和极小值包络面LE的曲面数据进行平均得到均值包络面ME，如公式（4）所示。

ME=(UE+LE)/2                       (4)

在步骤104，从滑动窗口内的干涉图的干涉相位中减去均值包络面。

具体地，从原始数据中减去均值包络面ME，得到第一个本征模函数(Intrinsic Mode Function,IMF)的估计值F_estimation，如公式(5)所示：

F_estimate=I-ME                          (5)

在步骤105，判断步骤104中计算得出的F_estimation是否满足预定条件。

例如，如果F_estimation的均值包络面的干涉相位值在每一个点处小于预先设定的阈值ε，则该F_estimation即是想要得到的IMF分量F_i；如果不是，则用F_estimation代替原始数据I，然后从步骤101开始进行迭代计算，直至得到满足预定条件的IMF分量F_i为止。

在步骤106，通过余量R_i进行多次分解，得到多个频率分量。

具体地，从原始干涉相位中减去满足预定条件的IMF分量，得到原始数据的余量信息R_i，如下面的公式6所示。

R_i=I-F_i                         (6)

余量包含了原始干涉图的大尺度信息，将R_i作为原始数据重复步骤二开始的过程，得到第二个满足条件的IMF分量。如此循环ｎ次，可以得到原始数据的ｎ个满足条件的IMF分量以及最终的余量。循环次数n可以根据要求预先人为设定，达到分解次数n后迭代运算终止。最终，原始数据可以表示为不同尺度的IMF分量和第n次分解后的余量R_n之和：

$I=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F}_i+R_n---\left(7\right)$

这样，经过二维EMD处理之后，窗口内的原始干涉相位被分解成了不同频率成分的分量，其中，第一个分离出来的分量F₁是原始数据的高频分量。而干涉相位中的噪声是非连续、非平稳的，同样属于高频成分，因此F₁主要包含了原始数据的噪声成分。

在步骤107，计算滑动窗口内的伪信噪比（Pseudo SNR）。

通过这些分量信息可以计算出原始数据的局部信噪比SNR，从而更加精确地估计统计窗口内的噪声程度。SNR的定义为信号与噪声的功率谱之比，而通常情况图像的信噪比可以近似为信号与噪声的方差之比，但是噪声的难以估计导致这一步很难进行。

这里，采用二维EMD分离出的F₁分量和原始数据I的关系，构造出一个新的参数——伪信噪比(Pseudo_SNR)，如公式(8)所示，这个新定义的参数能够更加准确地计算干涉相位中噪声所占的比重。

$\mathrm{Pseudo}\_\mathrm{SNR}=10\mathrm{lo}{g}_{10}({\mathrm{\σ}}_{\mathrm{\φ}\left(I\right)}^2/{\mathrm{\σ}}_{\mathrm{\φ}\left(F_1\right)}^2)---\left(8\right)$

为了提高计算效率和简化计算过程，优选地，可以只采用步骤108得到的IMF分量F₁来计算局部信噪比，即在公式（7）中取n=1即可。然而，本领域技术人员应当理解，为了更精确的计算结果，当然也可以采用步骤109中经过多次分解得到的多个IMF分量来构造伪信噪比。

在步骤108，根据伪信噪比计算出滤波参数。

众所周知，大信噪比的区域噪声程度低，小信噪比区域噪声程度高。在实际的SAR干涉图中，信噪比高的区域噪声程度低，因此需要的滤波强度相对较小；而信噪比低的区域噪声程度高，对应需要相对较大的滤波强度。根据上述特点，新的滤波参数α'可以根据Pseudo_SNR而设定：

α'=1-[exp^Pseudo_SNR/max(exp^Pseudo_SNR)]        (9)

在步骤112，用新的滤波参数对滑动窗口内的干涉相位图进行滤波。

首先，对滑动窗口中的原始干涉图进行离散傅立叶变换，得到原始干涉图的频谱：

Z(u,v)=F(z(r,a))             (10)

接着，采用新的滤波参数α'对取出的干涉图的频谱Z(u,v)进行平滑处理，得到平滑处理后的频谱：

H(u,v)=S{|Z(u,v)|}^α'·Z(u,v)          (11)

其中S{}为平滑算子。

最后，进行傅立叶逆变换把平滑后的频谱H(u,v)变化到空间域，从而完成了在滑动窗口内的局部滤波。

z(r,a)=F^-1(H(u,v))                     (12)

在步骤109，转移到下一滑动窗口重复整个滤波流程，直到完成整幅干涉相位图的滤波为止。

具体地，根据第一步设定的窗口大小和重叠转移到下一窗口重新进行从第二步开始的滤波流程，直到在原始干涉图的所有窗口完成上述处理，从而完成整幅干涉相位图的滤波。

在滤波结束，选取相位残差点和相位标准差作为评价参数，定量分析滤波结果。应当理解，根据本发明的干涉图滤波方法不限于应用在利用InSAR技术中，而是可以应用在各种场合下对干涉图的滤波中。

表1为唐古拉山冬克玛底冰川干涉图采用本发明的滤波方法与现有技术的定量分析对比，选取的指标为相位残差点和相位标准差。

表1

从表1可以看出，根据本发明的基于信噪比的InSAR干涉相位图滤波，能准确估计干涉图中局部噪声的程度，带来更佳的滤波效果。采用本发明的滤波结果，在边缘细节保持方面有很大的提高，在低噪声区域的降噪效果也更佳。总体而言，提升了传统Goldstein滤波算法降噪能力的同时，也降低了滤波对干涉图相位信息造成的损失，体现了本发明的滤波算法的整体优越性。

专业人员应该还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

结合本文中所公开的实施例描述的方法或算法的步骤可以用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器（RAM）、内存、只读存储器（ROM）、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。

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

Claims (10)

1.一种基于局部信噪比的干涉图滤波方法，所述方法通过窗口在所述干涉图中的滑动，依次对滑动窗口内的干涉图进行滤波，其中，在所述滑动窗口内对干涉图进行滤波包括：

对所述滑动窗口内的干涉图的干涉相位进行分解，得到所述干涉图在所述滑动窗口内的估计频率分量；

根据所述干涉相位和所述估计频率分量计算出所述干涉图在所述滑动窗口内的局部信噪比；

根据所述局部信噪比计算出滤波参数；以及

根据所述滤波参数对所述滑动窗口内的干涉图进行滤波。

2.根据权利要求1所述的方法，其中，所述对所述滑动窗口内的干涉图的干涉相位进行分解，得到所述干涉图在所述滑动窗口内的估计频率分量的步骤包括：

识别所述干涉相位的极大值点和极小值点；

根据所述极大值点和所述极小值点，计算极大值包络面和极小值包络面；

计算所述极大值包络面和极小值包络面的平均值；以及

根据所述滑动窗口内的干涉相位和所述平均值计算所述干涉图在所述滑动窗口内的估计频率分量。

3.根据权利要求2所述的方法，其中，所述根据所述滑动窗口内的干涉相位和所述平均值计算所述干涉图在所述滑动窗口内的估计频率分量的步骤包括：

从所述干涉相位中减去所述平均值，得到所述干涉图在所述滑动窗口内的估计频率分量；或者

判断所述频率分量是否满足预定条件；

如果不满足，则用所述频率分量代替所述干涉相位，从所述识别所述滑动窗口内的干涉相位的极大值点和极小值点的步骤开始迭代，得出满足所述预定条件的第一频率分量，将所述第一频率分量作为所述估计频率分量。

4.根据权利要求2所述的方法，其中，所述对所述滑动窗口内的干涉图的干涉相位进行分解，得到所述干涉图在所述滑动窗口内的估计频率分量的步骤还包括：

根据所述干涉相位和所述估计频率分量得到所述干涉相位的余量；

用所述余量代替所述干涉相位，从所述识别所述滑动窗口内的干涉相位的极大值点和极小值点的步骤开始预定次数的迭代，得出满足所述预定条件的一个或多个第二频率分量，将所述一个或多个第二频率分量作为所述估计频率分量。

5.根据权利要求1所述的方法，其中，所述局部信噪比是根据所述干涉相位和所述估计频率分量的关系计算得出的，所述局部信噪比是伪信噪比。

6.根据权利要求1所述的方法，其中，所述滤波参数是对所述局部信噪比进行归一化处理得到的。

7.根据权利要求1所述的方法，其中，所述根据所述滤波参数对所述滑动窗口内的干涉图进行滤波的步骤包括：

将所述滑动窗口内的干涉图从空间域变换到频域，得到所述干涉图的频谱；

用所述滤波参数对所述频谱进行频域平滑滤波；以及

将所述频域平滑滤波后的频谱变换到空间域。

8.根据权利要求2所述的方法，其中，所述极大值点和极小值点是通过像素点邻域搜索来识别的，所述极大值点和极小值点分别有一个或多个。

9.根据权利要求2所述的方法，其中，所述极大值包络面和所述极小值包络面是通过Delaunay三角剖分插值得到的。

10.根据权利要求1-9之一所述的方法，其中，所述干涉图是通过合成孔径雷达干涉测量InSAR得到的。

Images