循环中值平均滤除遥测地形粗差

摘要  对于地形遥测数据中引入的脉冲噪声，传统的滤波算法往往需要选择滤波门限，自适应能力不强，或者不能在滤除噪声的同时有效地保护信号数据，或者运算复杂。针对以上问题，文中提出了一种快速的脉冲噪声自适应滤除算法，该算法基于数理统计思想，对局部数据（滤波窗口）进行均值和方差估计，并根据估计结果自动选择检噪门限，进而实现噪声检测和平滑。实验结果显示，在脉冲噪声密度小于5%时，该算法的滤波信噪比增益远高于常用滤波算法和其他同类算法；对平稳变化信号的处理效果较好。关键词  地形遥测 粗差探测 自适应滤波 数理统计 1       引言地形遥测在地图制作、国土资源调查和地形辅助导航等民用和国防领域都有广泛的应用。随着科技的进步，地形测绘已由传统的人工测绘为目前利用星载和机载传感器自动测绘数字地图。在传感器遥测地形数据的过程中，不可避免地会引入某种程度的脉冲噪声，对于高精度传感器而言，脉冲噪声是测图过程的主要误差源，必须加以探测和处理。以平面扫描型机载侧视激光雷达为例，其扫描过程是在载机飞行过程中，对下方地形从左至右逐行扫描并获取对地斜距，根据测距数据、飞行位置和取向可生成航带区域的数字高程地图。因此，机载激光雷达往往用于数字地图测绘和地形辅助导航。如果把雷达测距数据按扫描行顺序排列，则可视作距离数据图像。因此，图像处理的各种方法都可应用于这类地形区域覆盖型传感器。目前常见的图像脉冲噪声去噪工具有中值滤波器（Median Filter）、粗差探测算子（如Robert算子和Laplacian算子）^[1]、低通滤波器^[2]，以及各种自适应滤波器。^[3]利用滤波窗口中间数据的平均来估计窗口均值，并采用软门限探测脉冲噪声，实质是中值和均值混合型滤波器。文献^[4]采用了与文献^[3]类似的脉冲噪声探测方法。文献^[5]则利用模糊隶属度函数对窗口数据加权平均，通过迭代得到收敛的加权均值，进而实现一种不等权的均值滤波。文献^[6] 提出了一种自适应模糊检噪平滑算法（AFND），该算法借鉴文献^[3][4] 方法首先对粗差进行探测，然后对确定为噪声的像点使用改进的文献^[5] 方法进行模糊加权均值平均以实现噪声点平滑；文献^[6] 对文献^[5]方法的改进主要体现在两个方面：一是基于局部参数估计在滤波前进行检噪预处理，二是将正态信号的分布性质与模糊隶属度函数的参数联系起来，从而解决检噪和滤波中的参数自动设定问题。中值滤波、均值滤波等平滑类算法不能提供对信号数据的有效保护；而粗差探测算子（硬门限）和采用软门限的检测类算法的检噪效果取决于门限或门限参数的选取是否适当，低通滤波同样需要对信号频带的先验估计，实质上其自适应能力都十分有限；AFND算法尽管在前两个方面做得较好，但运算复杂度高，实时性得不到保证。因此，对于需要较高滤波输出信噪比和自动化程度的实时地形测绘工作（如地形辅助导航），上述方法均不能得到令人满意的滤波性能。本文将地形数据图像的局部参数估计与传统的均值/中值滤波的简单算法有机地融合起来，构造一种自适应检测和平滑脉冲噪声的快速算法，该算法可以自动估计数据图像中的脉冲噪声级别是否下降到了允许范围，在条件不满足时自动循环迭代。下文从讨论如何利用图像的局部统计特性进行粗差检测着手，提出了循环中值平均检噪平滑算法，并与文中提到的几种其他滤波器进行了性能对比实验。2    脉冲噪声的统计法检测不失一般性，假定无污染遥测图像数据服从正态分布。对大量高程数字地图和激光测距图像信号的直方图统计表明，其数据近似服从正态分布，因此这一假设是合理的。对于局部图像（滤波窗口）中的信号数据，可看作是在该局部范围内无限密集采样总体的一组取样，因而窗口中的无污染信号数据是局部正态总体Ｔ的一组样本，且Ｔ～Ｎ（m,σ²）。类似地，窗口内的脉冲噪声数据可视为局部噪声总体Ｔ’的一组取样。为讨论方便，假设全局图像被密度为P_n且服从（a, b）上均匀分布的脉冲噪声所污染，对于均匀分布，局部噪声总体Ｔ’具有与全局噪声总体同样的分布。这样，检噪就成为随机分布总体的两类划分问题，此类问题多应用概率分类法解决。滤波窗口内Ｔ和Ｔ’的概率密度函数和概率分类的判决门限如图1所示：图 1 窗口数据的概率密度和判决门限 从图1可见，总体Ｔ和Ｔ’的分布往往是交迭的，无污染数据的判决门限界定为m±L，超出此范围的数据被判为脉冲噪声。对于两总体的交迭分布，无论判决门限如何选取，总存在两种错误判决的可能：第一类错误：x∈T，却被误判为x∈T’，设其概率为Ｐ₁；第二类错误：x∈T’，却被误判为x∈T，设其概率为Ｐ₂。Ｐ₁实质上表征了概率分类对信号数据的保护能力，而Ｐ₂则是对概率分类检出脉冲噪声能力的度量。其定量描述为：     （1） 公式（1）中的φ是标准正态概率分布函数，由公式可得到Ｐ_１、Ｐ_２和Ｌ的定性关系：（１）在σ、Ｐ_n、a、b不变的前提下，Ｌ增大时，Ｐ_１减小、Ｐ_２增大；Ｌ减小时，Ｐ_１增大、Ｐ_２减小；无论Ｌ如何选取，都不可能同时减小Ｐ_１和Ｐ_２。（２）判决门限Ｌ的选取与总体Ｔ的方差有很大关系，在保持P₁不变的前提下，T的方差σ越小，L值就越小。相对于整幅图像而言，局部总体Ｔ由于离散程度σ较小，因而在一定的P₁下使得Ｌ较小，可获得较小的P₂。由于图像的总错分概率是所有局部错分概率的平均，因而采用局部概率分类能获得比全局概率分类更优的检噪效果。 为方便，固定各窗口的门限比常数Ｒ＝Ｌ／σ，并根据统计推断原则将R取为2σ~3σ，这样可使得仅约0.26%~5%的数据点将被作为脉冲噪声进行平滑处理。这样，在脉冲噪声密度不高的情形下（高精度数据采集传感器工作在正常情况下），检噪算法对于未污染数据的保护性能较好。对于各滤波窗口中心点，用局部概率分类法进行噪声判断，需要估计局部均值m和局部方差σ，本文采用循环中值平均方法实现局部参数估计。3 基于循环中值平均进行检噪平滑对用于地形辅助导航的激光图像进行粗差滤波的要求是既能去除图像以外的脉冲噪声，同时又要保持图像的细节不被改变；此外还要求算法实时性强。AFND算法虽有较好的滤波效果，但涉及指数和迭代运算，运算速度低。中值滤波是最常用和简单的粗差滤波方法，但与以上激光像滤波要求相比也有差距。第一，中值滤波法及其改进型^[7][8]虽能减少图像中的脉冲噪声，但同时也改变了有效的图像数据。对此，^[5]采用了先进行噪声判断，再作针对性处理的方法。第二，中值滤波需要大量的排序操作，实时性还不够好。对此，文献^[3]采用的方法是取窗口中除去最大和最小值之后的中间数据的平均来代替中值。从这两种方法得到启发，提出利用中值平均进行局部参数估计，并进行统计法检噪滤波的方法，算法如下：（1）以当前输入检噪点I（i , j）为中心选取n×n的窗口（一般可取3×3窗口）；（2）找出窗内的最大值max、最小值min，并出窗内所有像点值之和S以及像点值平方和S2；（3）用窗内中段数据估计窗口非噪声数据的均值：      m=(S-max-min)/(n×n-2)；（4）用窗内中段数据估计窗口非噪声数据的方差：      σ=(S2-max²-min²)/(n×n-2)-m²；（5）粗差检测和平滑输出：      如果 |I(i,j)-m|>2×σ：则输出O (i , j)=m      否则：输出O (i , j)=I(i , j)该算法假设在一个滤波窗口内的高值脉冲和低值脉冲数量各不超过一个，当窗内脉冲数量不符合这一假设时，均值和方差估计误差增大，造成漏检。但由于每次估计都排除了最强脉冲的干扰，使估计更接近理想并检测平滑一定数量的噪声，如果反复执行这一过程，参数估计误差将趋于下降，滤波效果得到提高。因此设置了循环滤波方案，即当第一轮滤波检测出的噪声点总数与图像总像点数之比大于某一数值（循环阈值）时，进行第二轮滤波。以下为讨论方便，称这种方法为CMAND算法（Circular Median Average for Noise Detecting）。理论上，3×3的滤波窗口适合于粗差点的密度小于10%的情形，考虑到噪声分布位置的不均匀并通过实验，将循环阈值设置为0.05。４ 实验结果及讨论为了检验CMAND算法的性能，利用两组受到脉冲噪声不同程度污染的激光地形测距图像（见图2、图3所示，以灰度高低表示距离远近）将CMAND算法与其他算法进行滤波性能比较。由于从视觉直观上看不出各种滤波算法在处理效果上的明显差异，为了综合衡量滤波算法的数据保护和去噪效果，采用滤波信噪比增益（GSNR）作为滤波效果的客观评价性能指标，GSNR定义为滤波器输出信噪比与输入信噪比的比值，当噪声总能量远小于信号能量时，GSNR简化为输入噪声功率与输出噪声功率之比，即：公式中x、x_I、x_o分别代表无污染图像、脉冲噪声污染图像和滤波输出图像的像点数据。

图2   第一组滤波输入图像

图3  第二组滤波输入图像采用各种滤波算法，分别对上述两组输入图像进行处理，输出图像的GSNR值对比见表1、表2所示。表中Median表示中值滤波，Laplacian表示拉普拉斯检测算子，FWA表示模糊加权均值滤波器^[5]，AFND表示自适应模糊检噪平滑算法^[6]，CMAND表示本文提出的循环中值平均检噪平滑算法。 表1　各种滤波器在不同粗差密度下的GSNR值（对图2滤波）

DensityFilter	1%	2%	4%	6%	8%	10%	12%	16%	20%
Median	3.91	7.28	12.64	15.72	18.46	22.91	24.54	24.38	23.98
Laplacian	1.56	3.03	6.58	8.96	12.50	14.98	18.23	24.27	29.63
FWA	3.57	6.41	10.96	13.14	15.54	18.95	15.34	12.68	12.34
AFND	5.71	12.53	21.16	23.83	27.14	42.09	29.34	25.89	19.60
CMAND	12.19	26.00	31.14	39.17	47.92	54.59	26.42	30.16	18.19

    表2   各种滤波器在不同粗差密度下的GSNR值（对图3滤波）

DensityFilter	1%	2%	4%	6%	8%	10%	12%	16%	20%
Median	9.17	17.04	27.54	34.70	44.33	39.35	29.90	39.31	40.74
Laplacian	2.28	4.45	8.82	13.74	17.10	22.58	26.65	34.16	40.74
FWA	13.19	24.21	43.01	42.68	60.34	33.05	30.39	24.83	18.29
AFND	19.75	37.62	58.23	71.94	87.01	72.83	59.14	46.96	28.18
CMAND	44.68	70.59	99.38	114.1	105.3	128.8	59.22	40.23	43.30

 由表可见，AFND和CMAND算法的GSNR值在脉冲噪声密度低于10%时单调上升且远高于其余滤波算法，CMAND算法的表现尤为优异；而在噪声密度高于10%时滤波性能逐渐下降并接近常用的中值滤波和拉普拉斯算法。同时说明，AFND和CMAND算法对于较平稳过渡的图像具有较高的滤波增益。另外，CMAND算法是上述算法中速度最快的一种，在CELERON 333处理器和VC++ 6.0平台上，对180×180点的激光像进行滤波，几乎在瞬间完成，处理时间远低于0.1秒。因此CMAND是一种实时性和滤波效果都不错的算法。 ５ 结束语本文针对激光雷达对地形的遥测数据图像中所存在的脉冲噪声，提出利用中值平均粗略估计信号的局部统计参数，并根据概率统计理论设置局部检噪门限以检测和平滑噪声，通过循环迭代进一步降低噪声密度。在此基础上，提出了一种具有优良去噪性能和自适应能力的快速粗差检测型滤波算法（CMAND）。与脉冲去噪传统算法和其他粗差检测型滤波算法的对比实验表明，该算法在对低密度脉冲噪声的处理效果和处理速度上表现优异。为了对算法的适用场合作进一步的了解，将该算法应用于对混杂脉冲噪声的实验数字地图和一维平稳变化信号（对数字地图的一维采样）进行处理，获得了与以上实验类似的结果；算法作用于具有剧烈突变的人像图像，则性能下降较多并在低噪情况下接近中值滤波。由于篇幅所限，这里不再一一列举实验结果。综上所述，CMAND算法是适用于对混杂低密度脉冲噪声的平稳过渡信号进行高精度自动滤波处理的一种快速算法。由于激光测距雷达具有良好的虚警控制和测距精度，利用CMAND算法对其遥测的地形数据进行预处理，能实时获得远优于一般算法的滤波输出。 文献1   张永生. 遥感图像信息系统. 北京，出版社, 2000.2   Dan E Dudgern, Russell M Mersereau. MULTIDIMENSIONAL DIGITAL SIGNAL PROCESSING.    Prentice-Hall, 1984.李树涛, 王耀南. 图像椒盐噪声的非线性自适应滤除. 图形图像学报, 2000, 5(12): 999~1001.张政, 张宇等. 基于均值操作的快速自适应滤波器. 中国图形图像学报, 2000, 5(6): 530~534.蔡靖, 杨晋生, 丁润涛. 模糊加权均值滤波器. 中国图形图像学报, 2000, 5(1): 52~56.殷蔚明, 柳健, 田金文. 遥测地形数据中粗差的自动滤除. 遥感技术与应用,2001, 16(3): 158~162.Lin H M, Willson. Median filters with adaptive length.   IEEE CASI, 1988, 35(6): 675~690.8   Florencio D A F, Schafer R W. Decision-based median filter using local signal statistics.     Proc SPIE Symp Visual Comm Image Processing, Chicago, 1994, 2308: 268~275.