A Homogenous Pixel-weighted Interferometric Phase Filtering Method for Time-series INSAR

HUANG Junsong1, ZENG Qiming1,†, JIAO Jian1, GAO Sheng1,2

1. Institute of Remote Sensing and Geographical Information System, Peking University, Beijing 100871; 2. Institute of Electronics, Chinese Academy of Sciences, Beijing 100190; † Correspondence author, E-mail: qmzeng@pku.edu.cn

Abstract In order to effectively use the temporal information of multi-temporal INSAR and ensure the homogeneity of pixels in the window, a homogenous pixels-weighted interferometric phase filtering method is proposed. The method firstly uses the goodness-of-fit test to identify the statistically homogenous pixels (SHP) in the multi-temporal SAR images, and then the weighted interferometric phase filtering is performed just within those SHP, which ensures that the central pixels is not affected by the surrounding heterogeneous pixels. The experimental results based on real SAR interferograms show that compared with classic Goldstein filter and Lee filter, the proposed filter has advantages in visual effect, phase residues reduction, and phase derivative standard deviation, and achieves the original intention of the design of this filter. Key words time-series INSAR; weighted interferometric phase filtering; homogenous pixels; Goldstein filter; Lee filter

相位滤波是从时序 INSAR (Interferometric Synthetic Aperture Radar)干涉图中获取地表形变的一项关键技术, 旨在去除干涉图中的相位噪声, 提高相[1–2]位解缠精度, 改善形变估计结果 。传统的滤波[2–4]方法(如均值滤波、中值滤波、自适应滤波 、频率域滤波[5–6]、小波域滤波[7–10]等及其改进版本)对

SAR 图像空间信息的挖掘比较彻底(例如利用空间中像元的邻近性关系、信噪比大小、相干性大小等构造滤波权重), 在干涉图滤波中应用广泛。不足之处在于, 这些滤波方法未有效地利用时序 INSAR提供的时间维度信息。此外, 已有研究表明, 由于不同地物失相干程度不同, 在多时相 INSAR 干涉图中

携带的噪声水平也不同, 就矩形窗口类滤波(如均值滤波或干涉图多视)而言, 随着窗口扩大, 噪声相位的零均值随机分布不再成立, 导致滤波精度降低, 也就是说, 由像元的异质性导致的噪声水平差异对中心像元的相位估计是不利的[11], 因此有必要保证参与滤波的像元是同质像元。同质像元(statistically homogeneous pixels)的概念由 Ferretti 等[12]提出, 是基于像元间的统计特征而定义的。同质像元在时间序列上具有相同或相似的散射行为, 通常由同一类地物构成的像元间可以互称同质像元。

为了减少相位噪声, 提高相位解缠精度, 有效地利用多时相 INSAR 的时间维度信息, 并保证滤波窗口内像元的同质性, 在对比分析国内外各滤波方法的基础上, 本文提出一种适用于时序INSAR的同质像元加权干涉相位滤波方法。在本方法中, 同质像元的识别是利用多时相 SAR 图像的时间序列振幅信息, 通过拟合优度检验来进行的。此外, 在构造滤波器时, 考虑到各像元与中心像元的同质程度不一致, 本方法利用拟合优度检验计算出来的概率作为权重。为了避免解缠带来的复杂性和低效率,选择在相位的实部及虚部开展滤波。

虽然同质像元识别方法已广泛应用于时序 INSAR分布式散射体处理中, 对基于相位的实部及虚部滤波也有研究, 但针对已有滤波方法中存在的问题或缺陷, 将同质像元的思想应用到加权干涉相位滤波中尚属首次, 滤波的有效性需要实验验证。为此, 本文利用2007 年 1月到 2010 年 10月覆盖同一地区的27景SAR数据, 开展相位滤波实验, 并选用经典 Goldstein 滤波和Lee滤波作为参考方法, 从目视效果、残余点数及相位导数标准偏差方面评价各滤波器的滤波效果。

1 改进的滤波方法

本文提出的基于同质像元识别的加权相位滤波方法包含同质像元识别及加权相位滤波两个核心步骤。同质像元识别是滤波的前提, 多时相 SAR 数据时间信息的利用在此体现。同质像元识别之后,加权相位滤波则只在窗口内的同质像元间进行, 可以纯化参与滤波的像元, 保证像元的同质性。图1是本文改进的滤波方法的算法流程, 其中虚线框内的部分为同质像元识别流程, 检测概率 PK-S 是 K-S (Kolmogorov-smirnov)检验的输出概率, α是由实验确定的显著性水平。下面介绍图1中的关键步骤。

1.1 同质像元识别

同质像元通常是由同一类地物构成的像元, 代表时间序列上相同或相似的散射行为, 因此对同质像元的识别可以抽象为以下数学问题: 用拟合优度检验判定两组时间序列振幅观测值 Mp1 与 Mp2 (如式(1)所示)是否符合相同统计分布[12–13]。如果判定两个像元的时间序列振幅来自同一统计分布, 则接受二者是同质像元的假设, 否则拒绝。式(1)中, N为 SAR 图像的数量, mp1,1 和 mp2,1分别为像元 p1 和p2 在第一景 SAR 图像上的振幅, 依此类推。 N   M  m | m m , m ,  , m , p1 p1,1 p1,2 p1,  (1) 。 M  m | m  m p2,1, m , , mp2, p2 p2,2 N拟合优度检验分为参数化检验和非参数化检验。非参数化检验不需要对分布做先验假设, 具有鲁棒性相对较好的检测效果。通常认为单视 SAR振幅符合瑞利分布, 但在以往的研究中也发现有非瑞利分布[14]。为避免假设错误, 采用一种鲁棒性较好的非参数化拟合优度检验办法——K-S 检验[15]。K-S 检验的原理可以简要概括如下: 当 Mp1 和 M p2的累积分布函数(式(2)和(3))之差的绝对值取最大值DN (式(4))时, 判断DN是否落在某区间之内, 如果是, 则两组样本符合相同统计分布, 否则不符合。0, m m ， p1,1 k S p1( M p1) , m  m  mp1, , (2) p1,k k 1 N 1, m  mp1, ; N 0, m  mp2,1 , k S p2( M p2) , m  m  mp2,k , (3) p2,k 1 N 1,  m mp2,n ; DN  max S ( M )  S (M ) 。 (4)  M  p1 p1 p2 p2 DN的累积分布概率函数近似表达如下: P (D  D )  1  P ( D ) N K-S p1,p2

 0.11   1 P  D( Ne  0.12  ) , (5) K-S    Ne   Ne =N/2, N为 SAR 图像的景数, P ( t )  2  (  1) j 1  K-S j 1 e 2 j 2 t2 , 满足 PK-S(0)=1, PK-S()=0。当DN≤DTHR时, 即

P(DN ≤D)≤P(DN ≤Dthr)时, 两个样本服从相同统计分布, 而阈值 Dthr 与显著性水平α 相关, 二者关系有α =1P(DN≤DTHR)。式(5)化简后得到式(6)。

当 K-S 检验的输出概率 PK-S(DP1,P2) ≥ α 时, 接受零 假设 H0, 表明两个像元p1和 p2来自同一统计分布,是同质像元; 反之则接受非零假设 H1, 表明两个像元 p1 和 p2 来自不同统计分布, 不是同质像元。当N>8 时, K-S 检验结果非常可靠[12]。

1.2 加权干涉相位滤波

由于 INSAR 干涉相位丢失了相位整周期而缠绕在[−,+]之间, 是不连续的, 因此在利用诸如均

值滤波、中值滤波或其他低通滤波开展相位滤波时, 为了恢复真实相位, 要先进行相位解缠。但是,在滤波之前开展相位解缠, 大量相位噪声的存在会给解缠带来更多误差。因此, 本文选择在相位解缠之前进行相位滤波, 以便去除相位噪声, 提高相位解缠精度。事实上, 许多 INSAR 和时序 INSAR 工具包采取的也是这种策略, 例如时序 INSAR 工具包STAMPS 中的 Goldstein 滤波、INSAR 工具包 DORIS中的 spatialcov 滤波和 spectral 滤波都在相位解缠之前进行[16]。在相位解缠之前进行滤波, 要求构造的滤波器能够合理地处理缠绕的相位。Lee 等[3]针对缠绕和解缠的相位, 分别推导了对应的自适应滤波器。本文避免直接对缠绕的相位进行滤波, 主要考虑到缠绕的相位是不连续的, 并且能处理缠绕相位的滤波器的设计比较复杂, 存在较大难度。最终, 选择对相位的实部及虚部进行平均加权滤波(式(7)), 既能保证参与滤波的元素(即相位实部及虚部)是连续的,滤波器操作简单, 也可避免解缠。滤波后的相位实部及虚部需要重新组合成相位。

式中, Iˆp 和 Qˆ 分别为中心像元 p 的相位虚部及实

p部估计值, M为中心像元的同质像元数量, PK-S(DP, j)为中心像元 p 与像元 j 的 K-S 检测概率。由于窗口内参与滤波的像元都与中心像元互为同质像元, 因此像元的同质性有保障。此外, 考虑到 K-S 检测概率 PK-S(DP, j)刻画了像元 j 与中心像元 p 的同质程度, 因此取其概率作为滤波的权重。

2实验与结果分析2.1实验方案

本研究构建的实验方案如图2 所示, 阴影框部分示意本文改进的滤波方法的核心步骤(完整的算法流程见图 1)。将 2007 年 1 月至 2010 年 10 月覆盖北京西北部地区同一视角27景配准后的ENVISAT ASAR数据(由中国地震局地质研究所提供)作为输入数据(即N=27)。在干涉图生成过程中, 这 27 景ASAR 数据按照短时空基线原则生成 54 个干涉对。

分别用本方法、Goldstein滤波和 Lee 滤波, 对这54对干涉图进行相位滤波, 并根据目视效果、相位残差点及相位导数标准偏差, 对比和评价各滤波器的滤波效果。

2.2 评价指标

相位滤波效果的评价指标包括 3 个方面: 目视效果、相位残差点(又称相位奇点)数量及相位导数标准偏差。在目视效果方面, 要兼顾噪声去除(即平滑度)效果和相位细节的保真性。相位残差点的数量影响相位解缠精度, 点数越少, 对相位解缠越有利[4,17]。相位残差点其定义如下: 如图 3 所示,顺时针(或逆时针)计算 4 个像元的相位差 Δ1, Δ2, Δ3 和 Δ4 (式(8)), 如果 Δ1, Δ2, Δ3 和 Δ4 的值超过[−, +], 则对其值加/减 2, 使之重新缠绕到[−, +]之间, 记缠绕后的值分别为 Δ1, Δ2, Δ3和 Δ4。如

其中, xi j和 yi 分别为干涉图沿 x, y 方向的梯度值; , ,j xm 和 ym 分别为 k×k 窗口内 x, y 方向的梯度均值; ,n ,n Zm 为窗口中心像元的相位导数标准偏差, 是直接,n评价干涉图质量的有效测度, 其值越大, 说明干涉图质量越差[19]。

2.3 结果及对比分析

同质像元的识别采用 K-S 检验, 涉及显著性水平 α 及确定窗口的大小。一方面, 参与滤波的同质像元数量需要达到一定的规模, 因为数量太少会导致滤波不够准确, 因此滤波窗口不宜过小。另一方面, 实验区域包含地形起伏的山区, 如果滤波窗口太大, 相位细节特征会被平滑。兼顾像元数量和地形起伏两方面, 并考虑 SAR 图像分辨率(4.5 m ×7.8 m)和处理效率, 最终将窗口大小确定为 25×9, 对应的地面尺寸约为 100 m × 100 m, 大体上相当于研究区地表覆盖基本单元的尺度。对于显著性水平 ,其值越大, 表示两个像元互为同质像元的条件越严格, 反之则越宽松。实验中测试了0.05, 0.15, 0.25, 0.45, 0.65, 0.70, 0.80 和 0.90 等覆盖低、中、高区间的取值, 并统计同质像元的数量。从表1看出,当 α 取值为 0.65 时, 平均同质像元数只有36.00 个;当 α取值为 0.70 及以上时, 平均同质像元数更少,只有 7.52 个。因此, 当 α 取值为 0.65 及以上时, 滤波窗口内的像元数都太少, 不利于相位滤波; 而当α取值为 0.25 及以下时, 检测条件偏宽松, 容易将很多非同质像元检测为同质像元。最终, 将α 值确定为 0.45, 保证了比较严格的同质像元判别条件,此时的平均同质像元数为 74.31 个, 相对适中, 地物间的区分度也较好, 能较好地与地表相对应。

图 4 显示 α 取值为 0.45 时 K-S 检验识别出来的同质像元数量分布, 左下角灰度最深处为北京市中心城区, 建筑物密集, 是良好的点状散射体源,

1246 同质像元数一般较少, 大部分像元的同质像元为50个左右; 中间灰度较浅处为北京市郊区, 建筑物稀疏, 包含大量未开发的土地、草地等, 是良好的面状散射体源, 同质像元数一般较多, 大部分像元的同质像元在 150~200 个之间。

完成同质像元识别后, 利用构造的加权相位滤波器开展滤波实验。实验中不仅分析本文算法对干涉图的滤波效果, 还与 Goldstein 滤波和 Lee 滤波进行比较, 以便验证本文滤波算法的有效性和优越性。

Goldstein 滤波是先将窗口内的像元通过快速傅里叶变换转换到频域空间, 然后用选定的核函数对频谱振幅做平滑处理, 将平滑处理后的结果作为频谱的权重, 并进行傅里叶逆变换, 其本质是一种

频域内的加权滤波[5,20]。Lee 滤波则是一种带方向检测的自适应滤波, 利用相干性刻画的噪声水平,在高噪声水平处加大滤波强度, 以便获取足够平滑的结果; 在低噪声水平处降低滤波强度, 以便保留细节信息。此外, 为保留条纹信息, Lee 滤波通过搜寻条纹方向, 并沿条纹方向进行滤波[3]。作为两种经典滤波器, Goldstein 滤波和 Lee 滤波在 INSAR 测量中广泛应用, 因此, 与二者的对比具有重要参考价值。

由于本文面向的是时序 INSAR 中相位滤波方法的改进, 因此 Goldstein滤波参数α的设定参照著名时序 INSAR 工具包 STAMPS 给 Goldstein 滤波设定的值, 即 α=0.8 (α=0时无滤波, α=1时为最强滤波)。本实验中, 除保证三者滤波窗口一致(25×9)外, 不讨论参数的最优化问题。

图 5 显示与图 4相同区域的原始差分干涉图及滤波后的差分干涉图(13800×2300 像元)。从目视效果看, 本文方法与 Lee 滤波接近, 相对于原始干涉图有较大程度的平滑(图 5(a)中红色框所示), 条纹边界也更加清晰; Goldstein 滤波在目视效果上有所改善, 但不明显。

相位残差点数量越少, 对解缠越有利。相位导数标准偏差刻画干涉图的质量, 其值也是越小越好。从表 2 可以看到, 本文方法在相位残差点去除及相位导数标准偏差(Zm,n)两项指标上均为最优, 残差点减少率达到 65.8%, Zm,n平均值为 0.53, 远优于Goldstein 滤波的 30.6%和 0.65, 明显优于 Lee 滤波的 56.18%和 0.57。

图 6 显示与图 4相同区域的相位导数标准偏差Zm,n在原始差分干涉图及滤波后干涉图上的分布(与表 2 对应): 原始干涉图大部分在 1.0 以下; 本文方法滤波后干涉图大部分在 0.6 以下; Goldstein 滤波大部分在 0.8 以下; Lee 滤波大部分也在 0.6 以下, 但 0.4以下的蓝色区域像元数比本文方法少。

综合上述结果可知, 本文方法在目视效果、相位残差点数量及相位导数标准偏差方面均优于 Lee滤波和 Goldstein 滤波, 比原始干涉图有较大的改善。在这 3 个指标上, 本文方法与 Lee 滤波方法比较接近, 可能原因是 Lee 滤波是带方向检测的自适应滤波, 在图 5 中, 由于条纹信息比较清楚, Lee 滤波能够很好地找到条纹方向, 因此效果较好; 而本文方法也可在某种程度上视为一种自适应滤波, 滤波像元随窗口内同质像元数量而变, 且只在同质像元间进行加权平均, 因此既能保证滤波效果平滑,又能兼顾条纹信息。Goldstein滤波在本实验中效果不如前二者的原因, 可能在于频域内噪声与信号重叠较多, Goldstein 滤波器未能很好地加以区分。需要指出的是, 本文滤波方法计算速率相对较慢(劣于 Goldstein 滤波和 Lee 滤波), 将大量时间消耗在同质像元识别过程中, 存在较大的改进空间。

3 结论

本文提出适用于时序 INSAR 的同质像元加权干涉相位的滤波方法, 利用 27 景覆盖北京地区的时间序列 SAR 图像识别同质像元, 并在同质像元内对相位实部及虚部进行加权相位滤波。该方法纯化了窗口内参与滤波的像元特性, 保证了像元的同质性, 且滤波方法简单, 不需要解缠。通过 54 对干涉图的滤波实验, 得到以下结论。

1) 与经典的 Goldstein 滤波和 Lee 滤波相比,本文方法滤波后的干涉图在目视效果、残差点去除及相位导数标准偏差 3 个指标上均较优, 说明加入时间信息的相位滤波方法具有效性和优越性。

2) 同质像元的识别与滤波窗口大小及显著性水平相关, 二者的取值应适中, 应结合研究区地表覆盖类型、地表覆盖基本单元的尺寸以及图像分辨

率等综合考虑。

3) 地表像元的同质像元数量与其所属的地表覆盖类型相关。城区建筑物密集, 散射强烈, 同质像元数较少; 郊区的地表像元是良好的分布式散射体, 同质像元数较多。

本文为时序 INSAR 相位滤波时间信息的有效利用提供了一种新的思路, 未来可在更高效的同质像元识别方法研究、滤波器的构造、最佳滤波参数的选择方面进一步探索, 更好地应用于高精度地表形变反演的研究。 致谢 中国地震局地质研究所提供 SAR 图像数据, 研究工作得到北京大学高性能计算校级公共平台的支持。在此一并致谢。

参考文献

[1] Ferretti A, Prati C, Rocca F. Multibaseline INSAR DEM reconstruction: the wavelet approach. IEEE Transactions on Geoscience & Remote Sensing, 2002,

37(2): 705–715 [2] Bo G, Dellepiane S, Beneventano G. A locally adaptive approach for interferometric phase noise reduction // IGARSS’99 Proceedings. Hamburg, 1999: 264– 266 [3] Lee J S, Papathanassiou K P, Ainsworth T L, et al. A new technique for noise filtering of SAR interferometric phase images. IEEE Transactions on Geoscience & Remote Sensing, 1998, 36(5): 1456–1465 [4] 廖明生, 林珲, 张祖勋, 等. INSAR干涉条纹图的复数空间自适应滤波. 遥感学报, 2003, 7(2): 98–105 [5] Goldstein R M, Werner C L. Radar interferogram filtering for geophysical applications. Geophysical Research Letters, 1998, 25(21): 4035–4038 [6] Baran I, Stewart M P, Kampes B M, et al. A modification to the Goldstein radar interferogram filter. IEEE Transactions on Geoscience & Remote Sensing, 2003, 41(9): 2114–2118 [7] 蔡国林, 李永树, 刘国祥. 小波–维纳组合滤波算法

及其在 INSAR 干涉图去噪中的应用. 遥感学报, 2009, 13(1): 129–136 [8] 赵艳明, 全子一. 一种有效的小波–Wiener滤波去噪算法. 北京邮电大学学报, 2004, 27(4): 41–45 [9] Wang L C, Wang Y N, Mao L P. An algorithm of interferometric phase filter of INSAR based on wavelet analysis and median filter algorithm. Acta Geodaetica et Cartographic Sinica, 2005, 34(2): 108–112 [10] Braunisch H, Wu B I, Kong J A. Phase unwrapping of SAR interferograms after wavelet denoising // IGARSS 2000 Proceedings. Honolulu, 2000: 752–754 [11] Xia Y. Homogeneous pixel selection for distributed scatterers using multitemporal SAR data stacks [D]. Munich: Technical University of Munich, 2016 [12] Ferretti A, Fumagalli A, Novali F, et al. A new algorithm for processing interferometric data-stacks: squeesar. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(9): 3460–3470 [13] Lin K F, Perissin D. Identification of statistically homogeneous pixels based on one-sample test. Remote Sensing, 2017, 9(1): doi: 10.3390/rs9010037 [14] Moser G, Zerubia J, Serpico S B. SAR amplitude probability density function estimation based on a generalized Gaussian model. IEEE Transactions on Image Processing, 2006, 15(6): 1429–1442 [15] Kvam P H, Vidakovic B. Nonparametric statistics with applications to science and engineering. Hoboken: John Wiley & Sons, 2007 [16] Hooper A, Spaans K, Bekaert D, et al. STAMPS/MTI manual. Delft: Delft Institute of Earth Observation and Space Systems, Delft University of Technology [EB/OL]. (2013–09–12) [2015–10–24]. http://homepages. see.leeds.ac.uk/~earahoo/stamps/ [17] 曾琪明. 合成孔径雷达干涉技术研究与应用——以中国台湾集集921地震为例[D]. 北京: 北京大学, 2001 [18] Ghiglia D C, Pritt M D. Two-dimensional phase unwrapping: theory, algorithms, and software. Hokoben: Wiley, 1998 [19] Brown M, Harris C J. Neurofuzzy adaptive modeling and control. Englewood: Prentice Hall, 1994 [20] 严卫东, 倪维平, 赵亦工, 等. 自适应的改进Goldstein干涉相位图滤波算法. 西安电子科技大学学报(自然科学版), 2010, 37(2): 248–253

图 1改进的滤波算法流程Fig. 1 Flowchart of the improved filtering algorithm

图 3 相位残差点积分路线Fig. 3 Integral path of phase residual points

图 2实验方案流程Fig. 2 Flowchart of experimental scheme

图 4 α=0.45 时 K-S检验识别出来的同质像元数量分布(雷达坐标系) Fig. 4 Quantitative distribution of statistically homogenous pixels identified by K-S when α=0.45 (radar coordinates)