CN104103048B - 一种合成孔径雷达干涉图频域滤波的方法 - Google Patents

Abstract

本发明实施例公开了一种合成孔径雷达干涉图频域滤波的方法。所述方法包括：在对干涉图进行滤波前，利用滤波窗口中的SAR强度图和相干系数图生成联合搜索矩阵，在搜索矩阵中搜索与待滤波像素点统计特性相似的临近像素点，形成一个由统计特性决定的自适应邻域，再对这个自适应邻域内的像素进行滤波计算。本发明解决了现有滤波方法中因为滤波窗口固定造成的滤波结果计算不准确的情况，有效的提高了率滤波效果，尤其是干涉图中不同地物边缘区域的滤波效果。

Description

一种合成孔径雷达干涉图频域滤波的方法

技术领域

本发明涉及图像处理技术领域，尤其涉及一种合成孔径雷达干涉图频域滤波的方法。

背景技术

合成孔径雷达(Synthetic Aperture Radar,SAR)是一种主动式高分辨率空间微波遥感成像传感器，它具有全天时、全天候、高分辨率成像的能力，而且作用距离远，成像范围大，能够穿透植被和云层，解决了光学遥感图像容易受天气干扰的难题。合成孔径雷达干涉测量(Synthetic Aperture Radar Interferometry,InSAR)技术成功地综合了合成孔径雷达的成像原理和干涉测量技术，利用传感器的系统参数、姿态参数和轨道参数之间的几何关系等精确测量地表某一点的三维空间位置及其微小变化，是一种高精度对地观测新技术，它除了在地形测绘和地表形变监测方面具有重要的应用价值外，在极地冰川运动、海洋洋流运动、陆地分类、农业和资源调查等方面也具有重要的应用价值和潜力。

InSAR数据处理主要环节包括：基线估算、图像配准、干涉图生成、平地相位去除、相位解缠以及目标高程确定等。其中，相位解缠是InSAR技术获取地面高程、地表形变等信息的关键处理步骤之一。但目前相位解缠仍是InSAR应用的一个瓶颈问题，主要原因是干涉图质量对相位解缠的精度和可靠性影响较大。由于影响干涉图质量的因素有很多，如系统的热噪声、由于叠掩、阴影引起的时间去相关、空间去相关、多普勒中心去相关以及大气等因素，干涉相位图中常常存在严重的噪声。如果相位噪声不能得到有效的去除，一方面会直接传播到InSAR生成的地形图或形变结果中，另一方面会导致大量的相位残差点从而阻碍相位解缠的进程。因此必须要对相位干涉图进行滤波以降低相位噪声，提高其质量。

目前较为常用的滤波方法是Goldstein于1998年提出的频域滤波算法。方法对待滤波干涉图在滑动窗口内进行傅立叶变化，通过进一步分析噪声的高频响应，利用噪声不同于一般信号的频谱特性，对频率域的系数进行处理实现降噪的目的。但是该方法存在滤波窗口固定、缺乏自适应性的缺陷。在真实的SAR图像中，不同地物的后向散射系数不同，干涉形成的条纹图也会产生一定的差异。如果滤波过程中的固定矩形窗口内含有不同的地物类型，不加区分地采用矩形窗内的所有像素进行滤波计算，将造成滤波结果的不准确。因此，Goldstein滤波过程中滤波窗口的选取非常重要，对干涉相位的滤波效果有直接影响。

目前，Goldstein算法的核心步骤——傅立叶变换是在一个固定的矩形窗口中进行的。但是在真实的SAR图像中，不同地物的后向散射系数不同，干涉形成的条纹图也会产生一定的差异。如果固定矩形窗口内含有不同的地物类型，不加区分地采用矩形窗内的所有像素进行滤波参数估计，将造成参数估计的不准确，得到的滤波结果也会产生误差。因此，在对每一个像素进行滤波的过程中，如何选取合适的滤波窗口是一个值得研究的问题。

发明内容

有鉴于此，本发明实施例提供了一种合成孔径雷达干涉图频域滤波的方法，能够解决现有滤波方法中因为滤波窗口固定造成的滤波结果计算不准确的情况，有效的提高了率滤波效果，尤其是干涉图中不同地物边缘区域的滤波效果。

在第一方面，本发明实施例提供了一种合成孔径雷达干涉图频域滤波的方法，包括：

根据滤波窗口中的合成孔径雷达SAR强度图和相干系数图生成联合搜索矩阵；

根据预设初始点的邻域区域内像素点的统计结果确定所述联合搜索矩阵内的第一种子点的估计值；

依次将所述预设初始点的邻域区域内的每一像素点作为搜索点，确定所述搜索点与所述第一种子点之间的欧式距离；

如果所述搜索点与所述第一种子点之间的欧式距离小于或等于第一阈值，则将所述搜索点归入自适应邻域，并将所述搜索点作为自适应邻域搜索点进行自适应邻域的像素点搜索；

当所述自适应邻域搜索点的邻域区域内的像素点全部搜索完毕，或者所述自适应邻域内的像素点达到最大预定阈值时，计算所述自适应邻域内的像素点的联合矩阵均值，并将所述联合矩阵均值作为所述第一种子点的重估计值；

将所述预设初始点的邻域区域和所述自适应邻域搜索点的邻域区域内的像素点中除所述自适应邻域内的像素点之外的其他像素点进行重搜索，如果所述其他像素点中第一像素点的值与所述第一种子点的重估计值之间的差值小于或等于第二阈值，则将所述第一像素点加入所述自适应邻域；如果所述其他像素点中第一像素点的值与所述第一种子点的重估计值之间的差值大于第二阈值，则将所述第一像素点作为背景点；

根据所述自适应邻域和所述背景点重构滤波窗口；

根据所述重构的滤波窗口内自适应邻域的像素点的相干系数确定滤波参数；

根据所述滤波参数，对重构的滤波窗口内的像素点进行滤波处理。

优选的，所述将所述搜索点作为自适应邻域搜索点进行自适应邻域的像素点搜索具体为：

将所述搜索点作为自适应邻域搜索点，根据所述自适应邻域搜索点的邻域区域内像素点的统计结果确定所述联合搜索矩阵内的第二种子点；

依次确定所述自适应邻域搜索点的邻域区域内的每一像素点与所述第二种子点之间的欧式距离；

将与所述第二种子点之间的欧式距离小于或等于第一阈值的所述初始搜索点的邻域区域内的搜索点归入自适应邻域，并作为自适应邻域搜索点进行自适应邻域的像素点搜索。

优选的，在所述根据所述自适应邻域和所述背景点重构滤波窗口之前，所述方法还包括：

计算所述自适应邻域内的像素点的干涉相位均值；

获取全部背景点，将所述干涉相位均值作为所述背景点的相位值。

优选的，所述第一阈值的确定方法具体为：

其中，T₁为第一阈值，g(m,n)为搜索矩阵中第一种子点的估计值，σ_u为所述搜索矩阵的像素点的标准差，为所述搜索矩阵的像素点的均值。

优选的，所述根据所述重构的滤波窗口内自适应邻域的像素点的相干系数确定滤波参数具体为：

其中，α为滤波参数，为自适应邻域的像素点的相干系数均值。

优选的，所述根据所述滤波参数，对重构的滤波窗口内的全部像素进行滤波处理具体为：

H(u,v)＝S{|Z_s(u,v)|}^α·Z_s(u,v)；

其中，Z_s(u,v)为干涉图频谱，S{}为平滑算子，α为滤波参数，H(u,v)为平滑处理后的频谱。

优选的，在所述根据滤波窗口中的合成孔径雷达SAR强度图和相干系数图生成联合搜索矩阵之前，所述方法还包括：

确定滤波窗口的大小和初始位置。

进一步优选的，所述方法还包括：

变更所述滤波窗口的初始位置；

基于变更后的滤波窗口进行所述合成孔径雷达干涉图频域滤波。

通过应用本发明实施例提供的合成孔径雷达干涉图频域滤波的方法，在对干涉图进行滤波前，利用滤波窗口中的SAR强度图和相干系数图生成联合搜索矩阵，在搜索矩阵中搜索与待滤波像素点统计特性相似的临近像素点，形成一个由统计特性决定的自适应邻域，再对这个自适应邻域内的像素进行滤波计算。本发明解决了现有滤波方法中因为滤波窗口固定造成的滤波结果计算不准确的情况，有效的提高了率滤波效果，尤其是干涉图中不同地物边缘区域的滤波效果。

附图说明

图1为本发明实施例一提供的一种合成孔径雷达干涉图频域滤波的方法流程图；

图2为本发明用于冰川区域干涉图滤波结果对比分析图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为便于对本发明实时的理解，下面将结合附图以具体实施例做进一步的解释说明，实施例并不构成对本发明实施例的限定。

下面以图1为例详细说明本发明实施例一提供的合成孔径雷达干涉图频域滤波的方法，图1为本发明实施例一提供的一种合成孔径雷达干涉图频域滤波的方法流程图。其步骤包括：

步骤100，确定滤波窗口的大小和初始位置；

具体的，根据干涉图大小选取合适的滤波窗口。滤波窗口通常设定为32×32像素，也可以根据需要设定为64×64或128×128。同时确定并记录滤波窗口起始点在干涉图中的像素位置。

步骤101，根据滤波窗口中的合成孔径雷达(Synthetic Aperture Radar，SAR)强度图和相干系数图生成联合搜索矩阵；

具体的，获取滤波窗口中的SAR强度图A和相干系数图C，根据SAR强度图A和相干系数图C生成联合搜索矩阵g可以如下式表示：

g＝A*C＝[A_1,1C_1,1,A_1,2C_1,2,A_1,3C_1,3...,A_n,nC_n,n] (式1)

步骤102，根据预设初始点的邻域区域内像素点的统计结果确定所述联合搜索矩阵内的第一种子点的估计值；

具体的，在滤波窗口中预设有初始点，对预设初始点的邻域区域内的像素点进行统计可以确定联合搜索矩阵的种子值。

例如，以初始点为中心，将其8邻域区域内的像素点的联合矩阵中值作为种子点的初始值。

在其他例子中，还可以采用4邻域或其他领域区域内的像素点的其他统计结果来确定联合搜索矩阵的种子点。

步骤103，依次将所述预设初始点的邻域区域内的每一像素点作为搜索点，确定所述搜索点与所述第一种子点之间的欧式距离；

具体的，将上述邻域区域内的每一个像素点都作为一个搜索点，依次计算每个搜索点与种子点之间的欧氏距离。

步骤104，如果所述搜索点与所述第一种子点之间的欧式距离小于或等于第一阈值，则将所述该搜索点归入自适应邻域；

具体的，在上述邻区区域内的搜索点与第一阈值的关系满足如下条件时，将该搜索点归入自适应邻域。

||g(k,l)-g(m,n)||≤T₁ (式2)

其中，g(m,n)为第一种子点的值，g(k,l)为搜索点的值，T₁为第一阈值。

按照上述方法遍历邻区区域内的全部搜索点，将所有满足条件的搜索点都归入自适应邻域。比如在上述以初始点为中心的8邻域区域内，存在第一搜索点满足上述条件，则将第一搜索点归入自适应邻域。

步骤105，将所述搜索点作为自适应邻域搜索点进行自适应邻域的像素点搜索；

具体的，将上述满足条件的搜索点作为自适应邻域搜索点，根据所述自适应邻域搜索点的8邻域区域内像素点的统计结果确定联合搜索矩阵内的第二种子点；具体方法同上述步骤102；

随后，依次将自适应邻域搜索点的邻域区域内的每一像素点再作为搜索点，确定所述搜索点与所述第二种子点之间的欧式距离；

如果自适应邻域搜索点的邻域区域内存在第二搜索点与第二种子点之间的欧式距离小于等于第一阈值，将第二搜索点也归入自适应邻域。

按照上述方法遍历前述步骤104中所有归入自适应邻域内的搜索点。

在本方案中，第一阈值的选取非常重要。

在一个优选的例子中，第一阈值T₁可以通过如下方式来确定。

其中，T₁为第一阈值，g(m,n)为搜索矩阵中第一种子点的估计值，σ_u为所述搜索矩阵的标准差，为所述搜索矩阵的均值。

通过本步骤105，可以增大满足自适应邻域条件的像素的搜索范围。对初始点的8邻域内，满足加入自适应邻域的像素点如第一像素点的8邻域进行搜索。当然也可再对第一像素点的8邻域内，满足加入自适应邻域的像素点的8邻域进行搜索。。。直到再没有新的像素点满足一定的阈值条件，或者自适应邻域内包含的像素达到了预先设定的最大限制N_max。

步骤106，当所述自适应邻域搜索点的邻域区域内的像素点全部搜索完毕，或者所述自适应邻域内的像素点达到最大预定阈值时，计算所述自适应邻域内的像素点的联合矩阵均值，并将所述均值作为所述第一种子点的重估计值；

具体的，在预设初始点的邻域区域和所述自适应邻域搜索点的邻域区域内的像素点全部搜索完毕，或者当自适应邻域内的像素点达到最大预定阈值时，确定的所有自适应邻域内的像素点的联合矩阵均值更新为第一种子点的重估计值这样能够提高种子值的估计精度。

步骤107，将所述预设初始点的邻域区域和所述自适应邻域搜索点的邻域区域内的像素点中除所述自适应邻域内的像素点之外的其他像素点进行重估计；

步骤108，判断所述其他像素点中每一像素点的值与所述第一种子点的重估计值之间的差值是否小于或等于第二阈值；

具体的，根据所述预设初始点的邻域区域和所述自适应邻域搜索点的邻域区域内的像素点中除所述自适应邻域内的像素点之外的其他像素点判断，每一其他像素点的值g(u,v)与第一种子点的重估计值之间的差值是否满足：

(式3）

其中，T₂为第二阈值，可以根据预设来确定。

如果满足差值小于等于第二阈值，执行步骤109；否则执行步骤110。

步骤109，如果所述其他像素点中第一像素点的值与所述种子点的重估计值之间的差值小于等于第二阈值，则将所述第一像素点加入所述自适应邻域；

步骤110，如果所述其他像素点中第一像素点的值与所述种子点的重估计值之间的差值大于第二阈值，则将所述第一像素点作为背景点；

步骤111，计算所述自适应邻域内的像素点的干涉相位均值；

步骤112，获取全部背景点，将所述干涉相位均值作为所述背景点的相位值。

具体的，上述步骤111和112的目的在于，在经过上述各步骤后，被包含在自适应邻域内的像素点被认为具有相似的统计特性。滤波窗口内的像素点可以分为两类，一部分被包含在待滤波点的自适应邻域内，一部分则是背景点。将自适应邻域内的像素点保持原干涉相位值不变，将背景点的相位值则全部更新自适应邻域内像素点的干涉相位均值。经过这样的处理，使得在后续重构矩形窗口之后的接下来滤波时，实际参与滤波计算的只是自适应邻域内的像素点，而背景点在经过替换后没有实际计算意义，但是满足了窗口傅立叶的要求。

步骤113，根据所述自适应邻域和所述背景点重构滤波窗口；

步骤114，根据所述重构的滤波窗口内自适应邻域的像素点的相干系数确定滤波参数；

具体的，在重构滤波窗口之后，滤波参数α由重构的滤波窗口内所有自适应邻域内像素点的相干系数确定：

(式4)

其中为自适应邻域内像素点的相干系数均值。

步骤115，根据所述滤波参数，对重构的滤波窗口内的像素点进行滤波处理。

具体的，在重构的滤波窗口内进行Goldstein滤波。

H(u,v)＝S{|Z_s(u,v)|}^α·Z_s(u,v)

其中，Z_s(u,v)为干涉图频谱，S{}为平滑算子，α为滤波参数，H(u,v)为平滑处理后的频谱。

由此完成重建的滤波窗口中每一个像素点的滤波。

在步骤115之后，可以变更所述滤波窗口的初始位置，再基于变更后的滤波窗口进行所述合成孔径雷达干涉图频域滤波。

具体的，可以将滤波窗口向下一列移动一个像素；如果达到最大列，则移至下一行第一个像素，得到新的滤波窗口。

在新的滤波窗口内重复上述步骤102-115，直至完成整幅干涉图的滤波。

本发明实施例提供的合成孔径雷达干涉图频域滤波的方法，能够保证参与滤波计算的像素点具有相似的统计特性，从而提高了滤波效果。相比传统Goldstein滤波算法中，当固定滤波窗口内包含有不同地物类型时，不加区分地采用所有这些像素进行滤波计算，本发明实施例提供的方法更为合理，尤其会避免了在不同地物的边缘区域造成的滤波结果误差。

下述表1中给出了冰川区域干涉图经过本发明提供的滤波方案进行滤波后的结果，与传统Goldstein算法的定量对比分析，选取的指标为相位残差点

和相位标准差。

表1

图2中示出了冰川区域干涉图滤波结果对比分析图。

通过应用本发明实施例提供的合成孔径雷达干涉图频域滤波的方法，在对干涉图进行滤波前，利用滤波窗口中的SAR强度图和相干系数图生成联合搜索矩阵，在搜索矩阵中搜索与待滤波像素点统计特性相似的临近像素点，形成一个由统计特性决定的自适应邻域，再对这个自适应邻域内的像素进行滤波计算。本发明解决了现有滤波方法中因为滤波窗口固定造成的滤波结果计算不准确的情况，有效的提高了率滤波效果，尤其是干涉图中不同地物边缘区域的滤波效果。

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

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

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

Claims (6)

1.一种合成孔径雷达干涉图频域滤波的方法，其特征在于，所述方法包括：

根据滤波窗口中的合成孔径雷达SAR强度图和相干系数图生成联合搜索矩阵；

根据预设初始点的邻域区域内像素点的统计结果确定所述联合搜索矩阵内的第一种子点的估计值；

依次将所述预设初始点的邻域区域内的每一像素点作为搜索点，确定所述搜索点与所述第一种子点之间的欧式距离；

如果所述搜索点与所述第一种子点之间的欧式距离小于或等于第一阈值，则将所述搜索点归入自适应邻域，并将所述搜索点作为自适应邻域搜索点进行自适应邻域的像素点搜索；

当所述自适应邻域搜索点的邻域区域内的像素点全部搜索完毕，或者所述自适应邻域内的像素点达到最大预定阈值时，计算所述自适应邻域内的像素点的联合矩阵均值，并将所述联合矩阵均值作为所述第一种子点的重估计值；

将所述预设初始点的邻域区域和所述自适应邻域搜索点的邻域区域内的像素点中除所述自适应邻域内的像素点之外的其他像素点进行重搜索，如果所述其他像素点中第一像素点的值与所述第一种子点的重估计值之间的差值小于或等于第二阈值，则将所述第一像素点加入所述自适应邻域；如果所述其他像素点中第一像素点的值与所述第一种子点的重估计值之间的差值大于第二阈值，则将所述第一像素点作为背景点；

根据所述自适应邻域和所述背景点重构滤波窗口；

根据所述重构的滤波窗口内自适应邻域的像素点的相干系数确定滤波参数；

根据所述滤波参数，对重构的滤波窗口内的像素点进行滤波处理；

所述将所述搜索点作为自适应邻域搜索点进行自适应邻域的像素点搜索具体为：

将所述搜索点作为自适应邻域搜索点，根据所述自适应邻域搜索点的邻域区域内像素点的统计结果确定所述联合搜索矩阵内的第二种子点；

依次确定所述自适应邻域搜索点的邻域区域内的每一像素点与所述第二种子点之间的欧式距离；

将与所述第二种子点之间的欧式距离小于或等于第一阈值的初始搜索点的邻域区域内的搜索点归入自适应邻域，并作为自适应邻域搜索点进行自适应邻域的像素点搜索；

在所述根据所述自适应邻域和所述背景点重构滤波窗口之前，所述方法还包括：

计算所述自适应邻域内的像素点的干涉相位均值；

获取全部背景点，将所述干涉相位均值作为所述背景点的相位值；

在根据所述滤波参数，对重构的滤波窗口内的像素点进行滤波处理之后，变更所述滤波窗口的初始位置，再基于变更后的滤波窗口进行所述合成孔径雷达干涉图频域滤波。

2.根据权利要求1所述的方法，其特征在于，所述第一阈值的确定方法具体为：

<mrow> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>=</mo> <msub> <mi>&amp;sigma;</mi> <mi>u</mi> </msub> <mo>/</mo> <mover> <mi>u</mi> <mo>&amp;OverBar;</mo> </mover> <mo>&amp;CenterDot;</mo> <mo>|</mo> <mo>|</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> <mo>;</mo> </mrow>

其中，T₁为第一阈值，g(m,n)为搜索矩阵中第一种子点的估计值，σ_u为所述搜索矩阵的像素点的标准差，为所述搜索矩阵的像素点的均值。

3.根据权利要求1所述的方法，其特征在于，所述根据所述重构的滤波窗口内自适应邻域的像素点的相干系数确定滤波参数具体为：

<mrow> <mi>&amp;alpha;</mi> <mo>=</mo> <mn>1</mn> <mo>-</mo> <msub> <mover> <mi>r</mi> <mo>&amp;OverBar;</mo> </mover> <mi>s</mi> </msub> <mo>;</mo> </mrow>

其中，α为滤波参数，为自适应邻域的像素点的相干系数均值。

4.根据权利要求1所述的方法，其特征在于，所述根据所述滤波参数，对重构的滤波窗口内的全部像素进行滤波处理具体为：

H(u，v)＝S{|Z_s(u，v)|}^a·Z_s(u，v)；

其中，Z_s(u,v)为干涉图频谱，S{}为平滑算子，α为滤波参数，H(u,v)为平滑处理后的频谱。

5.根据权利要求1所述的方法，其特征在于，在所述根据滤波窗口中的合成孔径雷达SAR强度图和相干系数图生成联合搜索矩阵之前，所述方法还包括：

确定滤波窗口的大小和初始位置。

6.根据权利要求5所述的方法，其特征在于，所述方法还包括：

变更所述滤波窗口的初始位置；

基于变更后的滤波窗口进行所述合成孔径雷达干涉图频域滤波。