2. 山东科技大学 测绘学院，山东 青岛 266590
2. Geomatics College,Shandong University of Science and Technology,Qingdao 266590,China
1 引 言
InSAR图像由于成像方式的原因，不可避免地会产生噪声，而大量的颗粒状的相干斑降低了图像质量，影响了目标的检测、分类和识别，因此噪声去除对图像后期处理非常重要[1, 2, 3, 4, 5, 6, 7]。目前常用的降噪方法有多视滤波、自适应均值滤波、自适应中值滤波、自适应小波滤波和自适应圆周期中值滤波等[8, 9]。如多视滤波法可以抑制噪声，但是这种方法同时平滑了影像数据，是以损失图像空间分辨率为代价的，自适应均值滤波、自适应中值滤波和自适应圆周期中值滤波是传统滤波方法的改进，在去除噪声的同时可以保留更多的图像细节。还有一些滤波方法，比如粒子滤波[10, 11, 12]，或基于某种量化指标，如文献[13]是基于信噪比的滤波方法，文献[14]是基于特征保持的滤波方法。
本文提出利用泛函分析中的光滑子函数进行滤波，这是因为光滑子函数具有保凸性和光滑性等很好的性质[15, 16]，可以用来对函数进行磨光处理。在其与曲面相结合时，可将带有角点的函数磨光而又与原来的函数相切，在图像发生突变处使图像变平滑，产生“峰削尖、谷填底”的效果。
将InSAR图像的相位考虑为空间中的曲面，而斑点噪声使曲面发生突变，从而利用光滑子函数的性质对其进行滤波。对每一个点求出其4个方向的梯度，利用这4个梯度来判断这个点是否为噪声点。如果是噪声点，将噪声点所在的窗口进行曲面拟合，然后用光滑子函数对拟合曲面进行平滑处理；如果不是噪声点，则不进行平滑处理。通过这种自适应滤波方法，可以较好地抑制斑点噪声同时，又较好地保留了边缘信息并使相位偏移量较小。本文采用真实数据进行了滤波试验，并将自适应光滑子函数滤波和自适应圆周期中值滤波、自适应均值滤波、自适应中值滤波、自适应小波滤波的效果进行了分析和比较。
2 光滑子函数滤波原理 2.1 曲面构造在三维坐标系中，x、y分别表示图像的横坐标和纵坐标，Φ(x,y)表示图像的相位，给定图像中m×n个点。选定M个关于x的连续函数gr(x)，r=1,2,…,M，M＜m。同时，选定N个关于y的连续函数fs(y)，s=1,2,…,N，N＜n。且两组函数间线性无关，组成基函数。这些基函数可以构成以{ars}为系数的曲面族
根据最小二乘法则，式中系数[ars]M×N 可以表示为如下矩阵形式[17]
式中，G=[gr(xi)]m×M；F=[fs(yj)]n×N。将ars代入式(1)即得满足最小二乘原则的拟合曲面。通常情况下r和s在整数值2以内足以满足实际需要的精度。 2.2 自适应光滑子函数滤波
采用泛函分析中的二维光滑子函数对拟合曲面进行磨光处理，二维光滑子函数形式如下
式中，h为图像中心到边缘的距离。
如果没有噪声，干涉图的相位值具有很好的连续性。但是斑点噪声的存在破坏了这种连续性，滤波的目的就是恢复这种连续性。噪声的存在，使噪声点与周围的相位之间产生较大的梯度，因此可以利用梯度来判断干涉图中的突变点。
基于梯度的自适应法就是由干涉图的梯度来区分噪声点和非噪声点，进而决定是否要进行滤波处理。图像当中每一点的梯度为向量，本文中求出每一点梯度的模，并将梯度的模称之为图像在这一点的梯度。将带有噪声的干涉图进行分割，对于中心点，其梯度按照方向来分有4种，分别是：水平方向、左上右下方向、垂直方向、右上左下方向。公式如下
水平方向梯度 Gh[f(i,j)]=|f(i,j-1)+f(i,j+1)-2f(i,j)|
左上右下方向梯度 Gl[f(i,j)]=|f(i-1,j-1)+f(i+1,j+1)-2f(i,j)|
垂直方向梯度 Gv[f(i,j)]=|f(i-1,j)+f(i+1,j)-2f(i,j)|
左下右上方向梯度 Gr[f(i,j)]=|f(i+1,j-1)+f(i-1,j+1)-2f(i,j)|
对于中心点是否为噪声，采用这个点的上述4个方向梯度进行判断[18]。具体过程是，设置阈值k，求出4个梯度的最小值，如果最小值大于k，说明4个方向梯度都比较大，则认为该点为噪声点，将这个小的曲面进行拟合，并与二维光滑子函数相结合，进行平滑；如果最小值不大于k，可以认为不是噪声点，不参与平滑，从而得到噪声去除后的干涉图。数据处理基本流程如图 1所示。
|
| 图 1 数据处理流程 Fig. 1 Data flow processing |
对InSAR干涉图用4种滤波方法进行滤波去噪处理，通过边缘保持指数、峰值信噪比、均方根误差、相位标准偏差和不连续点等量化指标来对图像质量进行比较和评价[19, 21]。
3.1 边缘保持指数边缘保持指数(edge preserve index,EPI)反映了对图像进行滤波后的边缘信息保持能力，其值越大,表明滤波器的边缘保持能力越强。其表达式为
式中，Φs(i,j)是滤波后的相位值；Φo(i,j)是原始干涉图相位值。 3.2 峰值信噪比
峰值信噪比(peak signal-to-noise ratio,PSNR)是反映图像质量好坏的较为通用的指标。其值越大，去噪后的图像质量越高。其表达式为
式中，
。
3.3 均方根误差均方根误差(root mean square error,RMS)用来衡量滤波后相位图和参考相位图产生的偏差，其值越小代表滤波器的保真性越好。其表达式为
式中，Φs(i,j)是滤波后的相位值；Φo(i,j)是参考干涉图相位值；N是滑动窗口像元个数。 3.4 相位标准偏差
相位标准差(phase standard deviation,PSD)反映了相位的平滑程度，其值越小，说明相位越平滑。其计算公式为
式中，φ(i,j)是干涉图相位值；φ(i,j)是所选滑动窗口干涉图相位值的平均值。 3.5 不连续点个数
噪声去除后，干涉图像上的不连续点个数也是噪声去除能力的重要指标之一。通过滤波使图像不连续点减少到可接受的水平，为下一步相位解缠准备。
4 试验结果与分析 4.1 试验结果为了对光滑子函数滤波和现有滤波方法进行比较分析，本文采用了时间为2009-01-10和2009-02-25的北京ALOS数据，从中选取了400像素×400像素进行了试验。图 2给出了原始干涉图及滤波后的图像，其中原始干涉图中，梯度的最小值为2.050×10-8，滤波时阈值k取为5.6×10-5，表 1给出了各个降噪算法滤波结果的量化指标。
|
| 图 2 原始干涉图的降噪结果 Fig. 2 Filter results of original interferogram |
| 量化指标 | RMS | PSNR | EPI | PSD |
| 原始图像 | — | — | — | 1.814 4 |
| 自适应圆周期中值滤波 | 1.546 1 | 102.111 0 | 0.435 6 | 1.814 1 |
| 自适应中值滤波 | 1.430 1 | 103.671 0 | 0.357 7 | 1.696 2 |
| 自适应均值滤波 | 1.224 1 | 106.781 8 | 0.317 6 | 1.428 5 |
| 自适应小波滤波 | 1.319 8 | 110.910 9 | 0.397 5 | 1.627 5 |
| 自适应光滑子函数滤波 | 1.090 4 | 109.095 0 | 0.413 3 | 1.366 0 |
在原始干涉图和滤波后的干涉图中,选取第200行的相位值进行比较研究,其剖面图如图 3所示。
|
| 图 3 模拟的和滤波后的干涉图剖面( 第200 行) Fig. 3 Cross section over the simulated and filtered interferogram (row 200) |
从图 2(a)中可以看到，图像条纹比较清晰，但条纹的边缘仍有很多斑点噪声。经过光滑子函数滤波和其他4种方法滤波后，图像干涉条纹清晰，斑点噪声明显减少，说明这几种方法都可以有效地去除斑点噪声，同时又保留了图像的细节特征。
从图 3的剖面图可以看出，与其他滤波方法相比，经本文提出的滤波方法滤波后，干涉图的条纹曲线变化平稳且毛刺较少。这说明光滑子函数滤波方法不但能去除大部分噪声，而且对于较大颗粒的斑点噪声效果也非常好。
通过表 1中的量化指标可以对这4种滤波方法进行定量比较。在这5种滤波方法中，光滑子函数滤波方法的RMS最小，这说明用这种方法滤波后的相位图和参考相位图产生的偏差最小，即滤波器的保真性是这几种方法中最好的。从PSNR来看，经过自适应光滑子函数滤波后，PSNR要略低于小波滤波但是高于其他3种方法，说明这种方法在去噪后图像质量方面比小波滤波稍弱但是要好于其他3种滤波方法。EPI是指滤波器保持边缘信息的能力，从数据来看本文提出的方法在滤波后低于自适应圆周期中值滤波，高于其他3种滤波方法，这说明光滑子函数滤波在保持边缘信息方面的效果不如自适应圆周期中值滤波，但是要好于其他3种滤波方法。经过5种滤波方法去噪后，PSD均有下降，又以自适应光滑子函数滤波下降幅度最大，比原始图像减少了24.7％，自适应圆周期中值滤波降低很少，只减少了1.7％。这说明滤波后相位平滑程度均有所增加，但本文方法的相位平滑能力较强。
对干涉图去除噪声后，干涉图上的不连续点个数也是衡量滤波的重要指标值之一。不连续点个数在一定程度上反映了去噪后干涉图上的相位值高的颗粒噪声的数量多少。滤波后不连续点个数越少说明滤波方法抗相位畸变能力越强。表 2给出了原始及滤波去噪后，真实数据干涉图上断点的个数和通过滤波不连续点减少比例。
| 滤波方法 | 不连续点个数 | 减少百分比/(%) |
| 原始干涉图 | 13 062 | |
| 自适应圆周期中值滤波 | 5429 | 58.44 |
| 自适应中值滤波 | 4663 | 64.30 |
| 自适应均值滤波 | 1 | 99.99 |
| 自适应小波滤波 | 7685 | 41.17 |
| 自适应光滑子函数滤波 | 2 | 99.98 |
从表 2可以看出通过5种方法滤波后，不连续点都明显减少，自适应光滑子函数滤波和自适应均值滤波后不连续点最少，要明显少于小波滤波、自适应圆周期中值滤波和自适应中值滤波后的不连续点。这说明自适应光滑子函数滤波方法抗相位畸变的能力和自适应均值滤波相当，又要明显强于其他3种滤波方法。图 4给出原始干涉图及滤波后干涉图不连续点的分布图。
|
| 图 4 原始和滤波后干涉图不连续点分布图(白色为不连续点) Fig. 4 The original and denoised interferogram breakpoint distribution |
本文提出的基于梯度的光滑子函数进行滤波的方法是一种新的滤波方法。首先通过每个点的梯度来判断是否为噪声，如果是噪声，通过将包含噪声点的窗口进行曲面拟合，然后将拟合曲面与光滑子函数进行结合来平滑滤波，如果不是噪声则不进行平滑。在一定程度上消除了斑点噪声对图像造成的突变，从而达到了抑制斑点噪声的目的。利用真实数据进行试验，并与几种常用的滤波方法比较分析来看，本文的滤波方法在图像的保真性和相位平滑能力这几方面都是有优势的，在滤波后的图像质量方面好于自适应均值、中值及圆周中值滤波，在边缘保持能力方面要好于自适应中值、均值及小波滤波，抗相位畸变方面效果也较好。另外，光滑子函数还具有对空间目标进行平滑的作用，对以后的目标检测，分类和识别也具有一定作用。
| [1] | WANG A L, ZHANG Y, GU Y F. Simultaneous Speckle Reduction and SAR Image Compression Using Multiwavelet Transform[J]. Journal of Electronic Science and Technology of China, 2007, 5(2): 163-166. |
| [2] | FERNANDEZ S A, LOPEZ C A. On the Estimation of the Coefficient of Variation for Anisotropic Diffusion Speckle Filtering[J]. IEEE Transactions on Image Processing, 2006, 15(9): 2694-2701. |
| [3] | GOLDSTEIN T, BRESSON X, OSHER S. Geometric Applications of the Split Bregman Method: Segmentation and Surface Reconstruction [J]. Journal of Scientific Computing, 2010, 10(1): 272-293. |
| [4] | PI Yiming, YANG Jianyu, FU Yusheng, et al. The Principle of Sythetic Aperture Radar Imaging[M]. Chengdu: Press of University of Electronic Science and Technology, 2007. (皮亦鸣,杨建宇,付毓生,等. 合成孔径雷达成像原理[M]. 成都:电子科技大学出版社,2007). |
| [5] | MAITRE H. Traitement des Images de Radar à Synthèse dOuyerture[M]. SUN Hong, translate. Beijing: Publishing House of Electronics Industry, 2005. (MAITRE H. 合成孔径雷达图像处理[M]. 孙洪, 译. 北京: 电子工业出版社, 2005.) |
| [6] | LIAO Mingsheng, LIN Hui, ZHANG Zuxun, et al. Adaptive A1gorithm for Filtering Interferometric Phase Noise[J]. Journal of Remote Sensing, 2003, 7(2): 98-105. (廖明生,林晖,张祖勋, 等. InSAR干涉条纹图的复数空间自适应滤波[J]. 遥感学报, 2003, 7(2): 98-105.) |
| [7] | HAN Chunming, GUO Huadong, WANG Changlin. The Essence of SAR Image Speckle Suppression[J]. Journal of Remote Sensing, 2002,6(6) :470-474. (韩春明, 郭华东, 王长林. SAR图像斑点噪声抑制的本质[J]. 遥感学报, 2002, 6(6): 470-474.) |
| [8] | LIAO M S, LIN H, ZHANG Z X, et al. Adaptive Algorithm for Filtering Interferometric Phase Noise[J]. Journal of Remote Sensing, 2003, 7(2): 98-105. |
| [9] | WANG Xingwang. The Filter Method of InSAR Phase[D]. Changsha: Central South University, 2008. (王兴旺. 合成孔径雷达干涉图相位滤波方法的研究[D]. 长沙: 中南大学, 2008.) |
| [10] | TIAN Shijun, CHEN Jun, PI Yiming. Particle Filtering and Its Application for Kinematic GPS Positioning with High Speed Movement[J]. Acta Geodaetica et Cartographica Sinica, 2008, 36(3): 274-278. (田世君,陈俊,皮亦鸣.粒子滤波在高动态GPS定位中的应用[J]. 测绘学报, 2008, 36(3): 274-278.) |
| [11] | NIE Jianliang, YANG Yuanxi, WU Fumei. An Algorithm of Dynamic Precise Point Positioning Based on Modified Particle Filtering[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(4): 338-343. (聂建亮,杨元喜, 吴富梅. 一种基于改进粒子滤波的动态精密单点定位算法[J]. 测绘学报, 2010,39(4): 338-343.) |
| [12] | GONG Yisong, GUI Qingming, LI Baoli, et al. Design of Particle Filtering Algorithm Based on Mean Shift and Its Application in Navigation Data Processing[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(supplement):120-126. (宫轶松, 归庆明, 李保利,等. 基于均值漂移的粒子滤波算法设计及其在导航数据处理中的应用[J]. 测绘学报, 2011, 40(增刊): 120-126.) |
| [13] | SUN Qian, ZHU Jianjun, LI Zhiwei, et al. A New Adaptive InSAR Interferogram Filter Based on SNR[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(5): 437-443. (孙倩,朱建军,李志伟,等. 基于信噪比的InSAR 干涉图自适应滤波[J]. 测绘学报,2009,38(5): 437-443.) |
| [14] | YANG Shenbin, LI Bingbai, SHEN Shuanghe, et al. Structure Retaining Linear Multi Channel SAR Image Speckle Filter[J]. Acta Geodaetica et Cartographica Sinica, 2006, 35(4): 364-370. (杨沈斌,李秉柏,申双和,等. 基于特征保持的线性多通道最优求和SAR 图像滤波算法[J]. 测绘学报, 2006, 35(4): 364-370.) |
| [15] | YAO Zongyuan, QIU Chunhui, ZHONG Chunping. A New Leray Formula for Smooth Functions on Bounded Domains in Cn[J]. Science in China: Series A, 2002, 45(5): 610-619. |
| [16] | MAZROUI A, SBIBIH D, TIJINI A. A Simple Method for Smoothing Functions and Compressing Hermite Data[J]. Advances in Computational Mathematics, 2005,10(3): 279-297. |
| [17] | ZHANG Mengjun, SHU Hong, LIU Yan, et al. An Adaptive Thresholding Approach Based on Spatial Curved Surface Fitting[J]. Geomatics and Information Science of Wuhan University, 2006,31(5): 395-398. (张孟军,舒红,刘艳, 等. 基于空间曲面拟合的自适应阈值选取方法[J]. 武汉大学学报: 信息科学版, 2006, 31(5): 395-398.) |
| [18] | BOVIK A C. On Detecting Edges in Speckle Imagers[J]. IEEE Transactions on Signal Processing, 1988, 36(10):1618-1627. |
| [19] | GHIGLIA D C, MASTIN G A, ROMERO L A. Cellular-autom Atemethod for Phase Unwrapping[J]. Journal of the Optical Society of America,1987, 4(1): 267-280. |
| [20] | LIU Guoxiang. Monitoring of Ground Deformations with Radar Interferometry[M]. Beijing: Surveying and Mapping Press, 2006. |
| [21] | GOLDSTEIN R M, WERNER C L. Radar Interferogram Filtering for Geophysical Applications[J]. Geophysical Research Letters, 1998, 25(21):4035-4038. |