一种基于Non-Local Means的地基差分干涉雷达相位滤波改进方法

冯丽源 邓云开 聂祥飞 杨 鸿

（1.重庆三峡学院电子与信息工程学院，三峡库区地质环境监测与灾害预警重庆市重点实验室，重庆 404000；2.北京理工大学信息与电子学院，北京 100081）

1 引言

滑坡是我国常见的地质灾害之一，对其开展预防和预测是保障人员安全和国家经济财产的重要手段之一。地基差分干涉雷达通过微波遥感技术对隐患区域进行全天时、全天候的实时形变监测，已经在工程应用中得到了较好的应用［1］。但地基差分干涉雷达系统受各种扰动因素的影响，其原始干涉相位图中通常含有较多的相位噪声，严重影响后续数据处理，而干涉相位滤波技术可有效减少噪声干扰，改善相位图质量。

干涉相位图中，噪声相位是空变的，在空间上非平稳变化。受空变噪声影响，干涉相位图中不同区域的相位质量是不同的。研究人员对相位噪声滤波算法进行了大量研究，2016 年，岳顺等人［2］探究相位导数对地基合成孔径雷达（GB-SAR）干涉相位图的影响，提出了一种基于相位导数的GB-SAR影像自适应滤波方法，较好的抑制了噪声干扰；2018 年，胡晋山等人［3］对非局部算法进行改进，有效保留了干涉相位图的纹理和细节信息；2019 年，王念秦等人［4］证明了相位标准差模型在保持干涉图的细节和边缘信息具有明显优势；2020 年，樊伟等人［5］比较分析了几种高效的去噪算法，为形变数据处理指引了方向。

多视滤波器［6］被誉为最大似然下的最优滤波器，但该算法滤波窗口相对固定，从而使条纹密集区域的细节信息损失严重。Goldstein 滤波［7］作为工程中常用的滤波处理方法，可较好地保持相位条纹连续性，但该算法的幂指数相对固定，一定程度上减弱了算法对噪声的抑制效果。相对而言，Lee 滤波算法［8］的滤波精度较高，该算法采用了16 个可选择的方向窗对中心像素进行相位估计，但在实际处理中依然不能满足条纹方向的随机性，容易造成条纹断裂。维纳滤波［9］作为具有最优性能的最小均方误差准则，可有效滤除地基差分干涉雷达相位图中的噪声。Non-Local Means 滤波算法［10］通过比较整幅图像与待恢复区域窗口内各像素值来确定图像的相似块，然后将图像块聚类后进行统一滤波，该算法通过描述相位图结构特征之间的相似性，能更好的保护地基差分干涉雷达相位图的结构信息，最大程度地保留图像细节信息，为相位解缠提供有利条件，但传统Non-Local Means 滤波算法的平滑参数相对固定，容易使干涉相位图出现过滤波和欠滤波现象［11］。

针对以上算法的缺点与不足，结合干涉相位中噪声空变的特点，本文提出一种改进的自适应Non-Local Means 组合滤波算法，有效抑制了干涉相位图噪声，保持了干涉相位连续性，最终提升了相位图边缘模糊度。本算法利用相干系数构建自适应平滑参数模型，实现了窗口的自适应，通过与维纳滤波进行有效结合，改善了算法在高相干区域所引起的干涉相位图边缘模糊问题。最后通过实测数据验证了该算法的可行性。

2 干涉相位模型

2.1 相位域模型

干涉相位图是由两幅雷达图像复共轭相乘后取相位获得，两幅SAR 图像的相干性将影响高程测量的精度。两幅相干SAR 图像u1和u2的相干性可通过相干系数进行定义［12］：

上式（1）中，φi为获取的地形相位；φ0为相干系数相位；|ρ|为实相关系数，该值越大，图像相位质量越好。

若随机散射点满足复高斯分布模型，多视情况下分布式目标干涉相位的概率密度函可以表示为［13］：

上式中，β=|ρ|cos(φ-φ0)；N为视数；Γ(·) 表 示Gamma函数；F为高斯超几何函数。

由下式（3）计算出干涉相位的方差：

由干涉相位概率密度推导出加性噪声模型为：

式中，φb为相位测量值；φa为缠绕相位值；n表示方差为的0均值噪声。

2.2 复数域模型

为了保证相位解缠的顺利进行，需保留干涉相位的跳变点。干涉相位在-π 到π 内呈周期分布，相位域滤波会在一定程度上减少相位跳变点个数，为了有效保留相位跳变位置（信号高频位置），可将相位变换到复数域进行处理［13］，复数域表达式为：

实部和虚部的可表示为［14］：

式中，Nc=和αi表示与φa相互独立的加性噪声。

3 算法原理

3.1 Non-Local Means算法

非局部定义为空间中存在但彼此不相邻的相似图像块。将获得的相似图像块首先进行聚类，然后进行统一滤波是非局部的主要思想［15］。Non-Local Means 滤波算法的提出是图像去噪领域又一新的里程碑［16］。该算法可通过图像的相似性确定权值大小，将相似图像块进行加权平均可得到一个新的像素［17］。Non-Local Means 算法较好的运用了图像的冗余性，提升了图像质量，同时，使图像纹理信息更加清晰。Non-Local Means算法表达式为：

式中，f(x)表示去噪后的图像；w(x，y)表示权值；v(y)为原含噪图像；S为像素x的邻域。上式中，权值的确定是算法的关键，可由v(x)和v(y)的欧氏距离的平方确定，可表示为：

上式中，U(x)为归一化系数；h表示平滑参数。算法的权值和平滑参数的选择由欧氏距离决定。在实际算法的处理过程中，需设定搜索窗和邻域窗来确定像素权值，通过邻域窗在搜索窗内进行滑动，并与固定邻域窗口进行比较，确定邻域间相似度，计算权值，并以此方法遍历整幅图像。如图1所示，大窗口是以目标像素x为中心的搜索窗口，两个灰色小窗口分别是以x、y为中心的邻域窗口。其中以y为中心的邻域窗口在搜索窗口滑动，通过计算两个邻域窗口间的相似度，为y赋以权值w(x，y)。

3.2 维纳滤波算法

维纳滤波被称为最小均方误差下的最优估计器，对干涉相位图中的加性噪声具有明显的抑制作用［18］。维纳滤波中噪声的统计特性不随时间的变化而变化，可较好的去除高斯白噪声［19］。一个系统中，若原信号为x（n），含噪信号表示为y（n），则表达式如下：

式中，v(n)表示噪声信号。维纳滤波取得的x的估计值可定义为，其中，H表示维纳滤波器，可定义为：

上式中，Qyy(w1，w2)为噪声信号的自相关函数的傅里叶变换；Qyx(w1，w2)为噪声信号与真实信号的互相关函数的傅里叶变换。

3.3 本文算法

地基差分干涉雷达获取的影像数据被称为单视复数据，即SLC（single look complex）数据，SLC 数据包括雷达波振幅和相位两部分。针对地基差分干涉雷达相位噪声的空变特性，以及传统Non-Local Means 算法在平滑参数选择上的固定性，提出一种改进的自适应Non-Local Means 组合滤波算法。改进算法利用相干系数构建数据模型，实现滤波参数的自适应，在有效滤除相位噪声的前提下保持了相位连续性。另一方面，干涉图中不仅含有乘性噪声，而且含有加性噪声，而维纳滤波可在最优化准则下对干涉图的加性白噪声进行有效滤除。因此，本文将改进的自适应Non-Local Means 算法与维纳滤波进行有效结合，提高了噪声的滤除效率，改善了图像边缘模糊度。算法具体步骤为：

1）首先，将同一监测区域的两幅地基差分干涉雷达影像SLC1 和SLC2 进行复共轭相乘，获取相位后，生成干涉相位图，两幅影像的复数值用C1和C2分别表示：

两幅影像的相位差可表示为：

式中，Ai（i=1，2）表示复图像的灰度值；φi表示复图像的相位值。

2）然后，计算两幅影像的相干系数为：

3）通过相干值确定的参数模型为：

4）确定搜索窗口，选取邻域窗口。本文选择17×17的搜索窗口，根据式（17）得到的平滑参数的值自适应的选取邻域窗口大小，当平滑参数较小时，选择5×5的邻域窗口，可保留相位图更多的细节信息；当平滑参数较大时，选择7×7 的邻域窗口，可去除更多的相位噪声。最后，利用该点高斯邻域的相位均值替代去噪后的相位值。

5）确定参数之后对干涉图进行自适应滤波。

6）完成自适应滤波之后，结合维纳滤波进行处理。首先，获取复图像C1的互相关矩阵A1进行傅里叶变换得到Qyx(w1，w2)，然后获取C2的自相关矩阵A2进行傅里叶变换得到Qyy(w1，w2)，通过式（12）计算维纳滤波器H(w1，w2)，选取5×5的滤波窗口，对干涉图的实部和虚部分别进行维纳滤波，最后，将实部和虚部结合得到干涉相位图。综上所述，本文提出的改进的自适应Non-Local Means 组合滤波算法流程图如图2所示。

4 数据处理结果评估与分析

4.1 滤波效果评估指标

图像质量评估主要从客观和主观两个方面进行。主观评价主要是利用人主观上的认知意识和人眼视觉系统，通过对图像空间几何特征的识别进行评判；在客观评价方法中，一方面可通过残差点来判断是否具有较好的噪声滤除效果，另一方面，可通过计算边缘保持指数来评估图像是否较好的保留了边缘信息［20］；在客观评价方法中，主要看噪声的抑制效果，可通过剖线法和定量评价的方法进行评估。

1）边缘保持指数

边缘保持指数（edge preservation index，EPI）有效衡量了图像对边缘信息的保持能力，表达式如下［21］：

式中，φs(m，n)表示滤波后图像相位；φ0(m，n)表示原始图像相位值。EPI 的范围介于0～1 之间，该指标越接近1，滤波处理图像和原图像在保持边缘信息上越接近，进一步说明边缘保持效果越好，反之，则边缘效果保持越差。

2）残差点

在一幅干涉相位图中，残差点的数量决定了相位解缠能否顺利进行，相位图的残差点数量越多，越不利于相位解缠，通过干涉图的优化滤波可以有效提高图像质量，减少残差点个数，提升算法运算效率。残差点是对相邻4个点形成的闭环进行积分运算得到，如图3所示。

如上图所示，若想正确判定左上角的像元是否为残差点，需求得4 个相邻像元点的相位差的模和q，计算公式如下：

式中，Δi（i=1，2，3，4）表示相邻点之间的相位差；w(x)=x+2kp（-p≤w(x) ≤p）。若模和q不为0，则定义为残差点（q大于零为正残差点，反之，为负残差点）；若模和q等于0，则不是残差点。

4.2 实验及滤波结果分析

本文的实验数据来自重庆市彭水县马岩1号监测区，实验设备来源于北京理工大学雷达技术研究所开发的边坡监测雷达系统。图4所示为监测区域（红框区域）的雷达光学图像，获取的两幅SLC 影像数据大小为301×3702 像素，如图5 所示，其中，B 线表示剖面横断面位置，图5 中的相位呈现周期性变化，且噪声在纹理区域较为密集，不仅覆盖了大量细节信息，而且给后续处理造成困难。

实验选取三种典型的滤波算法与本文算法进行对比，实验结果如图6 所示。为进一步保证实验的客观性，在前期实验基础上，验证选取11×11窗口大小作为各算法最优处理窗口。从实验结果来看，自适应中值滤波算法对相位图的噪声滤除效果较差，原因在于该区域的相位噪声密度较大，而自适应中值滤波算法在处理噪声较为密集的图像区域时容易造成图像碎片化，从而影响了噪声滤除效果；Goldstein 滤波的噪声滤除效果优于自适应中值滤波，且图像清晰度和细节信息的保持能力也比传统Non-local Means 算法更好；传统Non-local Means算法的滤波效果相对较好，主要原因在于该算法可较好的利用图像的冗余信息，有效提升噪声滤除效率。结果显示，将其余三种算法与本文算法在宏观上进行对比，可以明显观察到，本文算法对噪声的抑制效果最好，且较好保留了噪声密集区域的有效细节信息。

经四种滤波方法处理后取相位图第230列数据（原始干涉相位图的B 线位置）的剖面横断面散点图，如图7所示，所选横断面穿过了噪声密度最大的区域，能代表性的衡量各算法的去噪能力。从图7可以看出，自适应中值滤波算法和传统Non-local Means 算法含有较多的毛刺，而Goldstein 算法所含毛刺较少，这是由Goldstein 算法的滤波特性所决定的，Goldstein 算法在处理噪声密度较大的形变连续区域时，使相位图中信号和噪声的频带相互分离，可有效提升算法处理效果，同时，从实验结果可以看出，Goldstein 算法的剖断面图像随相位的周期性变化较为明显。虽然Goldstein 算法的处理效果具有一定优势，但从实验结果来看，本文算法的剖断面毛刺更少，且平滑程度更高，相比于其余三种算法，优势更加明显。

四种算法滤波前后评估数值如表1 所示，为了能够准确客观的对本文算法的去噪能力进行评估，选取残差点和边缘保持指数作为衡量噪声抑制效果的有效评价指标。在残差点方面，四种算法均能取得不同程度减少噪声干扰，残差点滤除率（表中至上而下）分别为9.17%，28.9%，21.6%和53.4%。在边缘保持方面，Goldstein 算法最差，本文算法最优，原因在于维纳滤波可对随机噪声进行有效滤除，并且，通过自适应的设定边缘滤波参数，可有效改善图像模糊度，保留图像边缘信息。

表1 滤波算法评估结果Tab.1 The filtering algorithm evaluates the results

5 结论

本文提出了一种抑制干涉相位图噪声的有效方法，该方法充分利用Non-local Means 算法在图像去噪方面的优势，针对地基差分干涉雷达相位图噪声的空变特性，将Non-local Means 算法与维纳滤波算法进行有效结合，利用相干系数实现了Non-local Means 算法对滤波窗口的自适应选取，改进后的自适应Non-Local Means 算法更好的保留了干涉相位图中有用的细节信息，使相位连续性得到保持，随后，将组合滤波思想融入到改进的自适应Non-Local Means 算法中，不仅有效弥补了单一算法的不足，而且较好的改善了图像边缘模糊度。最后，将本文算法与其余三种滤波算法在实测干涉数据下进行综合分析，通过比较残差点个数、边缘保持指数以及滤波效果图，论证了本文算法的有效性。