数学形态滤波（共8篇）
数学形态滤波 篇1
0 引言
机载激光扫描(Li DAR)系统应用广泛,可快速获取地表的三维点云数据[1],在获取的点云数据中,包含有地面点与地物点,地物点主要包括建筑物、植被、电力线、车辆等信息。为获取地面点得到DEM,需要对点云数据进行处理,滤掉地物信息,该项工作称为滤波。
机载Li DAR数据滤波作为目前Li DAR数据处理的主要问题,国内外已经有许多学者对其进行了研究,如基于多次回波的滤波算法[1~4]、数学形态学的滤波方法[5~8]、迭代线性最小二乘内插法[9]、基于多方向的地面滤波算法[10]、基于TIN的滤波算法[11]等。在多次回波滤波算法中,主要考虑了多次回波特性对地物进行滤除,文献[1,2]主要基于双次回波进行滤波,主要通过双次回波特性首先分离出植被点,再进行建筑物等地物的滤除;文献[3]对森林地区的多次回波进行分析并滤波;文献[4]则通过多次回波特性预先剔除部分建筑物与植被点信息,减少数据量。基于数学形态学的滤波方法主要通过数学形态学的膨胀与腐蚀运算进行滤波,该类算法关键在于结构元素的选取,文献[5]通过改进的形态开运算对未填补空白区域的格网化DSM进行滤波,通过迭代计算实现点云分类;文献[6]考虑实际地形变化,根据形态学知识判断点间的高差与距离来实现地面点与地物点的分类,但是对于高差变化较大或者陡峭的地形具有一定局限性;文献[7]通过改进的形态学梯度计算每个点的梯度,再基于梯度值选取特定的点进行迭代开运算,根据梯度直方图减少迭代的次数,通过判断每次开运算后点的高程与原高程的差值是否小于一定的阈值,逐步滤除非地面点;文献[8]在数学形态学基础上,对格网化DSM进行腐蚀运算获得标记图像,对其反复进行测地膨胀实现开重建,利用白顶帽变换得到n DSM,实现地物与地面点的分类;文献[9]则使用线性最小二乘法内插激光脚点,根据拟合残差不服从正态分布的原则,使用迭代线性最小二乘内插法滤波与内插,并逐渐逼近和修正地形,该算法在地形变化较大的区域滤波时间较长,效率较低;文献[10]提出了基于多方向的地面滤波算法,对Li DAR数据内插出格网点,在扫描线上自动选取最低点作为地面点,通过该点的多个扫描方向的坡度,每一个点与最近邻点的高程差,每一个点与邻域最低高程点的高程差一系列的相关关系来判别数据点的类型,在城市地区与森林地区起到一个比较好的滤波效果;文献[11]采用改进后的三角网迭代滤波算法进行滤波,以待处理格网为中心,将待处理格网最低点与领域格网最低点构成三角网,再对格网内的数据与所在三角形的三个顶点进行比较,实现地面点判别,算法的难点在于不同地形阈值的确定。在这些算法中,部分算法仅考虑了双次回波而忽略了三次及以上的回波,或者适用于森林地区具有一定的局限性,或者需要先进行格网化再进行滤波。
分析Li DAR点云数据的多次回波特性,考虑现有算法思想,提出了考虑多次回波的数学形态学机载Li DAR滤波方法,滤波对象为原始离散Li DAR点云数据,首先提取多次回波中的首次回波与中间次回波作为种子点,对其进行区域增长,滤除掉大部分建筑物与植被,再采用较小结构元素进行形态学滤波,得到地面点,通过内插即可生成DEM。通过实例验证了滤波方法的有效性,具有较高的精确度和较好的实用性。
1 研究方法
1.1 多次回波分析
机载Li DAR系统发射激光脉冲,遇到物体时即返回信号,返回的回波信息包括单次回波与多次回波。同一束激光脉冲只发生一次反射则为单次回波,发生多次反射则为多次回波,其中对于多次回波而言,返回的第一个回波信号称为首次回波,返回的最后一个回波信号称为尾次回波,其余的回波信号皆为中间次回波。在实际中,两次回波之间的距离达到一定间距时,第二次回波才有可能被探测到[3]。目前,现有机载Li DAR系统都具备多次回波信息记录的能力,目前ASPRS激光委员会已制定并不断完善点云数据格式,即LAS格式,LAS数据的主要格式见表1,其中,回波次数为当前脉冲产生的回波总数,回波号表示当前数据为第几次回波。若回波次数为3,回波号为2,表示当前数据为第二次回波,当前脉冲共产生了3次回波。
LAS数据具有多次回波信息记录,分析多回波信息可以判断被测目标的类型:
(1)对于植被密集区域:首次回波大多为树冠部分反射形成,最后一次回波一般是较树冠低一些的枝叶或地面的反射信号[3,12],同时树冠部分或者地面也易形成单次回波信息,但是对于植被区域而言,多次回波数据中的首次与中间次回波一定为植被点。
(2)对于建筑区域,由于建筑物表面较为规则,高差较小,通常产生单次回波,在建筑物边缘区域高差较大,易产生多次回波[4],多次回波中的首次回波通常为建筑物边缘的屋顶点或者墙壁点,中间次与末次回波则可能为地面点或者墙壁点。
通过以上分析,多次回波的首次与中间次回波通常为植被点或者建筑物边界点,对其进行区域增长可以滤除大部分的植被点或者大面积的建筑物。
1.2 区域增长
通过对点云数据的多次回波分析,将多次回波的首次与中间次回波数据作为种子点,执行区域增长算法,即可获取大部分的地物点(主要为较大的建筑或者密集植被区域部分点)。由于滤波方法直接基于原始离散点云数据进行处理,点云密度高、数据量大,在进行区域增长之前,需要对点云进行组织,方式通常为重采样[13],如通过使用一个二维的格网覆盖在原始点云上,判断其落入哪一个格网单元内,每一个格网单元只保留一个最低点,当格网单元内没有点落入时,搜寻邻域最临近点作为该格网单元的点。这种重采样方式的缺点在于:当一个格网存在多个点时,由于只保留一个点,造成了精度损失。
为保持原始点云数据原貌,基于重采样思想[13],设计了一种网格分块链表的方法建立空间索引:
(1)计算原始Li DAR数据的平均间距d,构造一个合适的二维矩阵,矩阵的单元格大小与d相当,对每一单元格设置一空的链表,表示当前格网无点云数据;
(2)遍历Li DAR数据点,判断点云落入的格网单元,并对该格网单元加入链表索引,反复执行该步骤,直至所有点云判断完毕。
建立空间索引后,将空间区域按照平面坐标划分为方形小块,并基于行、列值对小块进行索引,在后期区域增长与形态学滤波时,对某一参考点的搜索可通过设置搜寻半径,限制在该参考点附近的子块当中进行,极大地提高了数据处理效率。建立空间索引后,即可对点云数据执行区域增长算法[14]:
(1)遍历点云数据,获取多次回波中的首次与中间次回波数据,将其作为地物点种子点;
(2)将种子点的高程同邻域点高程进行比较,若邻域点与种子点的高差小于给定的阈值,则标记该邻域点为地物点,并将其作为种子点继续进行增长;
(3)反复执行第(2)步的操作,直至所有点已判别并分类完毕。
1.3 形态学滤波
形态学图象处理是基于集合论对图像进行分析的工具,它着重研究图像的几何结构,基本思想是用具有一定形态的结构元素去量度和提取图像中的对应形状以达到对图像分析和识别的目的[9]。数学形态学的基本运算主要为膨胀、腐蚀、开运算、闭运算,若以A、C表示灰值图,A(m,n)表示图A在(m,n)的灰度,B为结构矩阵,B(i,j)为其在(i,j)的值。基本运算定义如下[10]:
定义1
A被B膨胀,记为AB:
定义2
A被B腐蚀,记为:
式中,表示对所有的B(i,j)取[·]中值的最大值。表示对所有的B(i,j)取[·]中值的最小值,结构元素可以为一维直线或者其他形状。对图像A先腐蚀后膨胀则为开运算,先膨胀后腐蚀则为闭运算。
基于数学形态学的滤波方法关键在于结构元素的大小设置,同时还与Li DAR点云数据内地物分布有关。合适的结构元素窗口执行一次开运算可以滤除掉所有的地物点,但是固定窗口大小的结构元素很难滤除不同类型的地物,结构元素窗口过小,只能滤除掉较小地物点(比如车、单株树木、小建筑物、行人等),大面积的房屋则较难滤除,当结构元素窗口过大,则易造成过度滤除,丢失地形特征信息。由于在多次回波中,较大的建筑物或者植被区域,已经通过区域增长进行滤除,剩余的地物多为较小的建筑物或一些低矮地物,因此,采用较小的结构元素窗口进行形态学滤波,获取地面点,避免了较大结构元素进行形态学滤波的过度滤除问题。
1.4 DEM内插
经过区域增长与数学形态学滤波,建筑物与植被被滤除,得到地面点数据,此时地面点存在大量数据空洞,为得到DEM,还需要进行内插处理,逐点内插法[15]以原始离散点云数据为基础,能较好地保留局部地形,该方法以待内插的数据点为中心,使用一个曲面函数去拟合邻域数据点,逐点内插法拟合函数常采用二次曲面进行拟合:
求解出二次曲面参数后,代入待插点的坐标值,即可求出该点的高程。
2 数据试验与分析
选取城市区域约480m×300m大小的原始LAS点云数据对本文提出的滤波方法进行验证,实验数据采用LAS文件原始格式,LAS版本为1.2,区域总点数为255074个,平均每平方米点数为1.77,点云平均间距为0.75m。该数据最大回波次数为4,产生的多次回波数依次为233948、20603、520、3,主要为首次、二次与三次回波,而四次回波数据只有3个。根据原始Li DAR点云数据生成灰度图像,见图1,从图中可以看出,选取的实验区域被大量建筑物与植被覆盖,且建筑物之间也存在大量树木,建筑物与植被等地物较为复杂,选取的数据较为典型。
在多次回波分析中得出多次回波的首次与中间次回波为植被点或者建筑物边界点,首先提取多次回波数据中的首次与中间次回波数据,提取原则为该数据点的回波次数(Number of Returns)与回波号(Return Number)不相等,且回波次数不为1。根据提取的结果数据生成灰度图像,见图2。对比图1可以看出,多次回波数据中的首次与中间次回波数据主要集中在密集植被区域与建筑物边缘处,与分析一致。
以提取的多次回波数据中的首次与中间次回波数据为种子点,执行区域增长算法,经过测试,高程阈值设置为0.3m,遍历点云数据,根据高程一致性原则进行区域增长,得到地物点,见图3。与图1对比,从图3中可以看出,绝大部分建筑物与植被点已被正确提取出来。
对比图3与图1可以看出,还有部分地物点没有滤除,主要为部分植被点及低矮地物点,对区域增长滤波的结果再采用数学形态学的方法进行滤波,结构元素窗口设置为5m,滤波后,采用逐点内插法内插出DEM,生成的DEM见图4。
从图4中可以看出,建筑物与地面点已经被滤除,得到的DEM平滑,基本保持了地表形状,得到的DEM质量较高,说明提出的考虑多次回波影响的形态学滤波方法较为理想,具有一定的实用性。
3 结论
本文通过分析机载Li DAR数据的多次回波特性,基于原始Li DAR点云,提出了一种考虑多次回波影响的数学形态学滤波方法,并用实例LAS点云数据进行算法实验,验证了该算法的适用性与有效性。根据研究,得出以下结论:考虑到多次回波总是发生在建筑物边缘区域或者植被区域,则其首次与中间次回波数据总为地物点,利用该特性可以对其进行区域增长,滤除掉大部分的地物点,对于剩余植被点或者低矮地物可采用较小结构元素使用形态学滤波的方法进行滤除,对其他滤波算法可提供一定的借鉴。不足之处在于对于地形突变的情况无法进行很好的滤波,这也是目前滤波算法普遍存在的问题,需要在其后的滤波算法中进一步改进。
数学形态滤波 篇2
关键词：形态学融合滤波；农业图像；MSR算法；局部像素最大化原则
中图分类号：TP391.41 文献标志码：A 文章编号：1002—1302（2016）01—0394—02
农业图像增强的根本目的是突出图像中的感兴趣信息，弱化其余信息，尽可能提高图像判读、分析的针对性。农业图像受野外成像环境多样性的影响以及在传输、解码过程中会存在不同程度的噪声并且在此过程中图像对比度也有所降低。对于农业图像的预处理，近年来学者们着重在滤波、增强2个方面进行针对性的研究，但是该类成果要么是针对图像中的噪声进行滤波，要么着重于进行图像增强，因而倾向性较为明显。当处理对比度较低且含有噪声的农业图像时，该类算法的处理效果则不尽如人意。因此，要实现对农业图像的有效处理，将图像滤波和增强算法进行有机融合是比较理想的选择。根据这一思路，本研究將形态学滤波与多尺度Retinex（muti-scale retinex，MSR）增强算法有机结合，提出了1种改进型MSR增强算法，即首先提出了1种形态学融合滤波算法对图像进行噪声滤除，然后对滤波后的图像进行MSR增强。
1算法原理
1.1形态学融合滤波
形态学图像处理的基本思路是采用预先设计的不同形状（圆形、矩形、菱形等）、不同尺寸（结构元素半径）的结构元素通过不同的运算方法来对图像进行处理和分析。形态学运算方法最基本的是腐蚀和膨胀运算，令函数F（i，j）表示任意一幅图像，B（i，j）为结构元素，腐蚀和膨胀运算定义为：
膨胀运算能够将图像中处于结构元素范围内的信息进行合并，对于图像中的空洞或凹陷部分（如裂缝）能够进行适当填补，能够滤除图像中负噪声点（噪声点灰度值明显低于图像中其余像素点灰度值），但是对于图像中的正噪声（噪声点灰度值明显高于图像中其余像素点灰度值）则无能为力。腐蚀运算则能够有效去除图像中的正噪声点，对于图像中小于结构元素尺寸且亮度较大的区域能够进行削弱甚至消除。因此，腐蚀和膨胀运算互补性较强，将二者进行有机组合，形成了开启运算、闭合运算：
开启运算能够有效去除图像中呈孤立分布的正噪声点（如图像中孤立存在的斑点、毛刺），从整体上平滑图像，但如果图像中的负噪声点过于密集且彼此间的距离明显小于结构元素尺寸，开运算处理的结果只能是进一步放大图像中负噪声点的分布区域；如果图像中的正噪声点较多，对图像首先进行膨胀运算然后进行腐蚀运算（即闭运算）也难以有效去除该类噪声点。开动、闭合运算的性能是基于采用同一形状、同一尺寸的结构元素得出的，进一步提高噪声滤波性能，最为有效的思路是采用不同尺寸的结构元素。这是因为采用尺寸较大的结构元素能更为有效地去除噪声点，但是会模糊图像；而采用尺寸较小的结构元素，尽管噪声去除能力下降，但能很好地刻画图像中的边缘轮廓信息。基于上述分析，本研究采用不同尺寸的结构元素（图1）将开启、闭合运算有机结合，并采用图像融合的策略，提出了农业图像形态学融合滤波的思路，具体步骤如下。
（1）采用如图1-a所示的尺寸为1的菱形结构元素对含有噪声的农业图像（尺寸大小为M×N）首先进行开启运算，然后进行闭合运算，得到滤波图像1。
（2）采用如圖1-b所示的尺寸为2的菱形结构元素对含有噪声的农业图像首先进行闭合运算，然后进行开启运算，得到滤波图像2。
数学形态滤波 篇3
心电图(Electrocardiogram,ECG)是临床上常规检查方法之一,它对某些疾病尤其是心血管疾病的诊断具有重要意义。心电信号作为心脏电活动在人体体表的表现,信号比较微弱,极易受环境的影响。其含有不同类型的噪声,主要有工频干扰和基线漂移。由于这些噪声与信号混叠,影响了心电各段波形特征的正确识别。为了消除心电信号中的主要干扰,提高检测准确率,人们提出了许多方法对心电信号进行预处理。如有限长单位冲激响应(Finite Impulse Response,FIR)滤波器、自适应滤波器、小波滤波器、 神经网络滤波方法和数学形态学滤波器等。自适应滤波器算法复杂,且需要附加参考信号[1]。小波滤波方法计算量大,处理时间长,不适于对算法实时性要求较高的场合[2,3]。 神经网络滤波方法计算复杂、速度较慢[4]。FIR滤波器虽然结构简单,易于实现,但由于ECG信号和基线漂移的频带相重叠,仅采用FIR滤波器方法在滤除噪声的同时,往往也会损失ECG信号中许多极有诊断价值的波形信息[5]。随着非线性滤波技术的发展,数学形态学提供了一种非常有效的非线性信号处理方法,其建立在积分几何和随机集合论基础上的,根据信号的局部特征对信号进行分析和识别, 可以用于ECG信号的滤波处理,在滤除噪声的同时可以较好地保持必要的心电几何信息不变。
本文概述了FIR滤波器和数学形态学的应用于心电信号滤波的原理,结合FIR滤波器和数学形态学滤波方法的优点,设计出基于FIR滤波和数学形态学的综合滤波方法。 该方法首先采用形式简单的FIR平滑滤波器滤除心电信号的50 Hz工频及其高频谐波,接着将数学形态滤波器应用于滤除基线漂移。本文采用PTB标准数据库[6]进行研究, 该数据库工频干扰为50 Hz,采样频率为1000 Hz。实验结果表明,本文设计的综合滤波方法能够有效地滤除工频干扰和基线漂移,为心电信号进一步的分析提供良好的基础。
1 FIR平滑滤波和数学形态学在心电信号预处理原理
1.1 FIR平滑滤波器
FIR平滑滤波是数字滤波方法中常被人们采用的方法, 该方法算法简单,处理速度快,滤波效果较好[7,8]。使用FIR平滑滤波器对信号滤波时,实际上是拟合了信号中的低频成分,而把高频成分“平滑”出去。由于人体心电信号的频率较低,主要频率范围是0.05~100 Hz,而大部分能量又集中在0.5~45 Hz。可以采用FIR平滑滤波器滤除工频噪声及其高频谐波。
本论文采用PTB标准数据库进行研究,该数据库工频干扰为50 Hz,采样频率为1000 Hz。将FIR平滑滤波器运用于滤除该数据库心电信号的50 Hz工频及其高频谐波,该滤波器阶数必须为N=1000 Hz/50 Hz=20,即该滤波器传递函数为:
该FIR平滑滤波器对50 Hz工频及其高频谐波截止, 且对100 Hz以后信号基本衰减到原信号的10%,即 -10 d B (图1)。将该滤波器应用于PTB数据库中s0001_re.dat文件的滤除噪声后结果见图2。本文设计的20阶FIR平滑滤波器可以有效的滤除心电信号中的50 Hz工频及其高频谐波, 对100 Hz以上的高频噪声抑制效果也不错。
1.2 数学形态学滤波器
FIR平滑滤波器虽然有效的滤除心电信号中的50 Hz工频及其高频谐波,但信号存在基线漂移,对信号的检测及特征提取影响很大。本文将数学形态学运算应用于一维信号处理中,利用这种非线性滤波方法对于噪声带来的奇异点敏感性来滤除心电信号中的噪声,实现对一维的ECG信号进行数学形态学滤波的目的(图3)。
设一维ECG信号的数字化序列为, 形态学结构元为,且有N>M。则ECG信号可经过下面的转化与二维图像相统一 :令二维欧氏空间图像集A(图3中网格部分)的包络为f(n),则信号A在N处取值为
同理,对结构元素序列k(m) 有
则信号f(n) 关于结构元k(m) 的形态学膨胀运算定义为:
信号f(n) 关于结构元k(m) 的形态学腐蚀运算定义为
对一维信号的开运算记为;闭运算记为。图4(a) 给出对一段ECG信号进行开运算处理后的结果,实线为运算前的ECG信号,虚线为对ECG信号开运算后结果。将ECG信号减去开运算后的波形,突出ECG的峰值。图4(b) 给出了对ECG信号进行闭运算处理后的结果,实线为闭运算前的ECG信号,虚线为对ECG信号闭运算后结果。将ECG减去闭运算后的波形,突出ECG谷值。
从图4可以清楚地看出,从原始ECG信号中分别减去其形态学开运算或闭运算后的结果,就可以得到原信号的峰值或谷值,这些波峰或波谷的宽度取决于所选择的结构元宽度。
对于含噪ECG信号而言,如果选择结构元的宽度小于ECG信号所有特征子波形的宽度,则对ECG信号进行开运算和闭运算后,ECG信号的所有特征子波形都会被保留,而信号中混杂的宽度小于结构元宽度的高频干扰则会被滤除。另一方面,有基线漂移的ECG信号可以认为是缓慢变化的信号上叠加宽度相对狭窄的ECG信号,因此,也可以用一组开、 闭运算从原始信号中有选择地除去ECG特征波形,开运算移去正脉冲,而闭运算移去负脉冲。从而得到从原始信号中分离出来的基线漂移信号,再用原始信号减去基线漂移信号后, 得到除去基线漂移的矫正后的ECG信号。
因此,本文选择了幅值为0的直线形的结构元素,设计了包含两组串联的数学形态学滤波器模块,这两个数学形态学滤波模块分别具有不同的结构序列。第一组结构序列宽度较大,其序列宽度大于心电特征波形P、Q、R、S、T宽度, 进行开闭运算的结果,使得这些特征波形都被滤除,只剩下基线漂移信号,再用原信号减去获得的基线干扰信号,即可获得滤除了基线漂移干扰后的心电信号 ;第二组形态滤波器的结构序列宽度较窄,其序列宽度大于高频噪声信号宽度, 而小于ECG特征波形P、Q、R、S、T的宽度,则开闭运算的结果,使得心电信号中的高频噪声信号被去除。这样在经过这两组不同的形态滤波器后,获得了滤除了高频干扰和基线漂移的信号(图5)。
采用上述数学形态学滤波器,对PTB标准心电数据库中的s0001_re.dat的第I导联的部分心电信号数据进行形态学滤波实验。对该信号采用第一组数学形态学滤波器模块滤波后结果见图6,其中开运算的结构元素的宽度M=190, 闭运算的结构元素的宽度M=70。
对滤除基线漂移后的心电信号采用第二组数学形态学滤波模块,此时结构元宽度M为6(图7)。从图8可以看到,第二组数学形态学滤波模块采用结构元宽度M为6对去除基线漂移后的心电信号进行去噪,虽然基本能够去除高频噪声,但是还是有部分宽度较大的噪声不能去除。改变第二组数学形态学滤波模块结构元宽度,加大结构元宽度令M=15,高频噪声去除效果良好,但是形态学滤波器在处理高频干扰时产生了一种近似矩形或梯形的小波动,使得ECG信号在高频小信号范围内产生了失真。
2 综合滤波算法及实验结果
实验结果表明,FIR平滑滤波器算法在处理高频干扰信号时,有很好的处理效果,但无法滤除基线干扰信号。 而形态滤波算法在滤除基线干扰信号时,有较好的效果, 但在滤除高频干扰信号时,则会产生截断误差。本文将这两种滤波算法相结合,提出了基于FIR平滑滤波器法与数学形态学滤波法相结合的综合滤波算法。该算法首先用FIR平滑滤波器对心电信号进行处理,去除50 Hz工频及其高频谐波,同时抑制频率大于100 Hz的高频噪声,获得了去除工频高频和高频噪声的输出信号,然后采用上一节设计的数学形态学滤波器滤除基线漂移,并进一步滤除高频噪声,该算法流程见图8。
采用上述综合滤波算法,在对PTB标准心电数据库中的s0001_re.dat的第I导联的部分心电信号数据进行滤波实验。运用综合滤波处理算法后的对ECG信号进行滤波,滤波效果比单独采用FIR平滑滤波器或数学形态学滤波器好。 从实验结果比较可以看到,该算法不仅能够去除基线漂移、 50Hz工频及其高频谐波,同时可以把宽度很小的噪声滤除, 且不会产生了近似矩形或梯形的小波动,出现高频小信号失真(图9)。
3 结论
本文从理论上论述了FIR平滑滤波和数学形态学应用于心电信号预处理的优缺点,并综合两者的优点,提出综合滤波算法,该算法能有效地滤波工频及其高频谐波、基线漂移和其他高频噪声。
摘要：本文针对心电信号的工频干扰和基线漂移,提出一种基于有限长单位冲激响应(FIR)滤波器和数学形态学的综合滤波方法,该方法首先采用形式简单的FIR平滑滤波器滤除心电信号的50 Hz工频及其高频谐波,接着将数学形态滤波器应用于滤除基线漂移。实验结果表明,本文设计的综合滤波方法能够有效地滤除工频干扰和基线漂移,为心电信号进一步的分析提供良好的基础。
基于拉普拉斯滤波的形态相关器 篇4
因为能够避免复数滤波器的制作和严格对准,联合变换相关器(Joint transform Correlator,JTC)目前已成为使用最为普遍的光学相关方法。传统的JTC输出通常伴有较宽的旁瓣和强烈的中心自相关峰干扰,造成识别度的下降。此外,JTC的输出性能对背景噪声和输入图像的光照变化也很敏感。许多方案已被提出用来改进传统JTC的性能[1,2,3,4,5]。在这些改进型的JTC中,有些方案[3,4,5]提出在JTC装置的联合功率谱(Joint Power Spectrum,JPS)强度平面上引入一些强度型滤波器,从而达到提高相关峰强度、锐度和识别度的目的。联合小波变换相关器(Joint Wavelet Transform Correlator,JWTC)就是这样一类相关器[6,7,8,9]。在JWTC中,一个经适当缩放的小波母函数的频谱强度谱函数被挑选出来对JPS进行滤波操作。滤波后的谱函数再经变换透镜做傅里叶变换,得到最后的相关结果。通过这样的处理,JWTC能够同时完成图像的边缘特征提取和相关操作。目标图像的边缘特征蕴含着模式识别所需的关键信息,因此利用小波滤波的边界提取功能能够进一步提高相关器的识别性能。边缘信息通常包含在图像的频谱高频部分,因此小波母函数的强度滤波处理等价于一个高通或带通滤波。由于拉普拉斯算子的能谱具有很明显的高通结构,因此具有与小波滤波类似的图像边缘提取功能。一种基于拉普拉斯算子滤波的JTC已被报道用来提取图像的边缘特征[10]。
随着近年来可编程光电设备性能的飞速发展,人们提出了一种叫做形态相关器(Morphological Correlation,MC)的非线性相关器[11,12,13]。形态相关函数的最大值点对应于两幅图像的平均绝对误差(Mean-absolute Error,MAE)的最小位置。众所周知,传统线性相关(Linear Correlation,LC)最大值出现在两幅图像平均平方误差(Mean-square Error,MSE)最小的位置处。随着偏离最小误差位置,MAE比MSE增加得越快。此外,同MSE相比,MAE对某些非高斯噪声,如椒盐噪声和冲激噪声,具有较强的抗噪性。MAE的这些特性使MC具有更高的识别度和相关峰锐度,并对椒盐噪声等非高斯噪声产生更强的抑制作用。而且,当输入图像(其中包含参考目标)的照明亮度强于原有参考图像时,MC的相关峰强度将会保持不变,即可实现照明亮度不变识别[12]。MC可以视为多个普通LC的叠加,因此利用传统的JTC型装置就可光学实现MC[13]。
本文提出了一种基于拉普拉斯算子的形态相关器(Laplace Filtering Based MC,LBMC)。这种相关器具有比LC和MC更为优良的输出特性。与原有MC不同,LBMC在执行最后的傅里叶变换之前,要对JPS预先作拉普拉斯能谱滤波处理。本文通过计算机模拟对其性能作了初步的验证,取得了较好的效果。
1 形态相关器
设t(x,y)和r(x,y)分别是灰度输入图像和参考图像。两者间的MC定义如下:
这里⊗是相关运算。tq(x,y)和rq(x,y)分别是经如下阈值分解而成的二值化图层片:
式(1)第一个定义式表明MC可视为所有来自对应相同灰度级q的二值图tq(x,y)和rq(x,y)的线性相关叠加[13]。根据该定义,一种基于JTC的光学置备系统被提出用来实现MC。第二个定义式是MC的原始定义[12]。从中可以看出MC执行的是参考图像与输入图像间的最小取值操作。与LC相比,MC具有如下优点:1)较强的抗椒盐噪声能力;2)能产生更尖锐得相关信号;3)参考图像的自相关峰强度总是高于参考图像与错误图像间的互相关峰强度;4)当输入图像(包含参考目标)的照明亮度强于原有参考图像时,MC的相关峰强度是照明变化不变的。
尽管如此,现有MC仍然存在着一些问题。首先,MC的最小化操作使其相关强度总是低于相同条件下LC的相关强度;其次,虽然MC具有性质2),但当输入错误目标的亮度远高于参考目标时,错误相关峰和正确相关峰之间的差异将会迅速减小,从而导致识别度下降;其三,MC的相关峰锐度仍不足以使其能够精确地定位小目标的位置;最后,MC的抗椒盐噪声能力仍然局限在较低的噪声强度范围内。
2 基于拉普拉斯滤波的形态相关器
拉普拉斯算子如下所示:
图1展示了它的能谱,其通频带分布在四角方向。该强度滤波器能使四角方向的高频分量顺利通过。拉普拉斯算子滤波可视为二阶求导运算,因此其零节点可用来探测边缘位置。除拉普拉斯算子外,罗伯特算子也可用于边缘特征的提取。罗伯特算子有两种,可分别用来探测45°和135°方向的边缘特征。使用单个罗伯特算子不能完整提取所有四个方向的边角信息。显然,使用单个拉普拉斯算子就能完成两个罗伯特算子才能完成的边缘提取任务。
图2是LBMC的光电系统实时制备装置。CCD1首先记录下输入图像t(x,y),并将其传送到计算机中执行阈值分解。拉普拉斯能谱|L(u,v)|2(u,v是空频坐标)和参考图像的阈值分解图层片rq(x,y)事先已经被制备好,因此不会影响系统处理速度。接着,联合图像
被传输到空间光调制器SLM1上。此时,紧贴SLM1的空间光调制器SLM2保持完全透明。经透镜傅里叶变换,CCD2接收到对应于灰度级q的联合功率谱JPSq。按这种方式可陆续得到对应所有灰度级的JPSq。所有JPSq叠加,最后结果传输到SLM1上。此时|L(u,v)|2被显示到SLM2上。滤波后的JPS可用下式表达:
这里Rq(u,v)和Tq(u,v)分别是rq(x,y)和tq(x,y)的傅里叶频谱。对其再做傅里叶变换,得到结果:
式中:“*”是卷积符号,L是拉普拉斯算子。图像边缘的特征提取与最后的形态相关运算经一次傅里叶变换就可同时完成。出现在坐标(x-2d,0)附近的第三项就是LBMC的输出结果。它可以看成是经拉普拉斯算子滤波后的所有二值图片rq(x,y)和tq(x,y)的相关叠加。输出平面中心的直流自相关项可以使用多种方法有效消除[8,14,15,16],这里不再详加讨论。
3 模拟测试与讨论
本节给出了使用Matlab程序,对LBMC的性能进行模拟测试的结果。同时LC和MC的模拟结果也被给出,以便比较。图3(a)给出了一个T90坦克作为参考目标(128×128,64灰度级)。图3(b)是用于测试的输入图像,除参考目标外,其中还包含了一个M1A2坦克作为错误目标。注意错误目标亮度要高于参考目标。图4展示了三种相关器的输出结果。LC的错误信号明显强于正确信号,因此已经失去了识别能力。MC和LBMC都能识别出参考目标。但MC的正确相关信号和错误相关信号差异很小,这导致其识别度不高。在做相关运算之前,LBMC预先提取了目标的边缘特征信息,这使其能够产生强而尖锐正确相关信号。另一方面,LBMC中的形态相关操作又保证其正确的相关信号总是强于错误参考信号。表1列出了这些相关器的自相关强度(Autocorrelation Peak Intensity,API)、峰-旁瓣比(Peak-to-Sidelobe Ratio,PSR)、自相关-互相关峰强度比(Autocorrelation Peak Intensity to Cross Correlation Peak Intensity Ratio,ACR)、半高宽(Ful Width of the Autocorrelation Peak at Half of Its Maximum,FWHM)。注意表中所有API和PSR都依照LC的相应值做了归一化。所有参数都表明LBMC具有更为优良的输出性能。
(a)参考图像;(b)输入图像;(c)被椒盐噪声破坏的输入图像(强度0.45).Fig.3 Experimental result of T90 tank(a)Reference scene;(b)Input scene;(c)Input scene corrupted by salt-and-pepper noise(intensity is 0.45).
下面给出椒盐噪声影响时的模拟结果。图3(c)是受椒盐噪声(强度0.45)干扰的输入图像。图5是相应的三维相关输出。可以看出,MC和LBMC的抗噪性能明显优于LC,LC的信号已被噪声完全湮没,尤其是LBMC,其相关峰尖锐明显,仍能保证ACR为4.35的识别度。表2第一部分给出了不同噪声强度影响下,三种相关器的峰噪比。在大部分测试范围内,LBMC的结果都要远高于另两种相关器。这主要是因为LBMC中的形态相关运算能够有效的压制椒盐噪声的干扰,而拉普拉斯滤波又保证其能产生强而尖锐的相关信号。表2第二部分给出了LC和MC的ACR随噪声强度的变化。在整个测试范围内,MC的ACR值变化不大,但非常低,表现出较差的识别度。尽管LBMC的ACR随噪声增加有大幅的下降,但即使在高强度噪声影响下仍保持了很高的数值。这进一步表明LBMC综合了MC和拉普拉斯滤波的优点。虽然拉普拉斯滤波因其高通性,而对噪声比较敏感,然而形态相关所蕴含的取最小操作成功的抵消了这种敏感性。
我们最后给出照明测试的结果。图6展示了当参考目标T90坦克分别乘上照明因子α=0.4,0.7,0.9,12,3(由左至右)时,三种相关器所产生的相关峰外观形状。LC的相关强度与α2成正比。当α>1,MC的相关峰强度不随照明因子变化。由于拉普拉斯算子滤波对照明因子比较敏感,在α>1部分,LBMC不具有MC的照明不变识别特性。尽管如此,LBMC仍能维持强而尖锐的相关峰和识别度。我们计算了不同照明因子影响情况下MC和LBMC的ACR值(如图7)。虽然LBMC的ACR随着照明因子偏离1而迅速下降,其值仍比普通MC的结果要高出许多。这保证了在几乎整个测试范围内MC能够准确无误地排斥虚假目标。形态相关的定义决定了越亮的像素对相关输出的贡献越大,因此在α>1区间的LBMC相关峰强度下滑速度要比α<1部分和缓。另一方面拉普拉斯滤波对正确相关峰和错误相关峰间的相对差异起到了展宽扩大的效果。两种因素综合导致LBMC的ACR值在α>1时下降趋近于某个恒定值(该值与参考图像的选择有关)。MC的ACR显然不能保证其能毫无差错地完成模式识别任务。我们使用其它参考目标进行了类似的模拟测试,也获得了相同的结论。
4 结论
数学形态滤波 篇5
“飞思卡尔”杯全国大学生智能车竞赛规则明确指出，智能车在赛道上连续跑两圈，并记录其中最好的单圈成绩，这使路径记忆算法成为可能。如图1所示，赛道记忆算法在第一圈以最安全的速度缓慢驶过一圈，并将赛道信息保存下来，第二圈根据保存下来的信息进行车速和转角决策的相应最优化，从而在第二圈取得好成绩。无论智能车的传感器前瞻距离有多远，在跑圈时它都只能预测在一段有限距离内赛道的情况。而采用赛道记忆算法的智能车，在第二圈时已对整个赛道有了全面的认识，从而在相同条件下，将比不使用赛道记忆的智能车更具优势。
第一圈准确记忆赛道信息是第二圈控制策略的基础，是比赛成败的关键。但是在第一圈中不论控制策略如何优秀，赛车总会或多或少的偏离赛道，舵机的转角信息总会出现一定程度的毛刺和扰动等粗大误差，其幅值足以引起处理器的误判。如果不加处理直接应用于第二圈控制，会对赛车造成严重干扰，不能以极限速度跑完比赛或者冲出赛道。通常的处理方法有两种：一是通过阈值比较丢弃干扰值，但同时赛道信息也会同干扰信息一起被丢弃;二是通过低通滤波将干扰平均到多个位置，但同时破坏了赛道原始信息，而且在干扰幅值很大的时候效果也不是很明显。
数学形态学 (Mathematical Morphology) 是一种新型的数字图像处理方法和理论，其主要内容是设计一整套的变换 (运算) 、概念和算法，用以描述图像的基本特征。提供了非常有效的非线性滤波技术，该技术只取决于信号的局部形状特征。因此，它在诸如形状分析、模式识别、视觉校验、计算机视觉等方面，要比传统的线性滤波更为有效。
算法的选取与实验结果对比
数学形态学的运算以腐蚀和膨胀这两种基本运算为基础，引出了其它几个常用的数学形态运算。数学形态学中最常见的基本运算只有七种，分别为：腐蚀、膨胀、开运算、闭运算、击中、细化和粗化，它们是全部形态学的基础。它们的定义如下：
设X代表一个数字图像，我们假定该图像是二值的，即取值只有1或0，则X表示一个二进制信号集合，B是一个简单的紧集合，称为“结构元素”。X被B膨胀和腐蚀的结果可以分别定义为：
而X被B进行开运算和闭运算的结果分别定义为：
在数字图形处理领域中，数学形态学主要用于非线性变形，它可以局部地修改信号的几何特征，并提供有关信号的几何特征信息。根据不同的信号的形态特征，可以采用不同的数学形态学运算对信号进行处理，这些数学形态与运算都被视为数学形态滤波器。在这种应用方法中，每一个信号都被视为适当的维数的欧几里德空间中的集合。数学形态滤波器被定义为集合的运算，它使信号的图形变形，以提供关于其几何结构的数字化信息。对于被视为集合的二进制信号，腐蚀、膨胀、开运算和闭运算是最简单的形态运算。这些滤波器还可以引申到多维信号中去。此时，形态滤波器利用的是灰值图的数学形态运算的定义。下面将探讨如何将数学形态滤波器应用到舵机转角信号 (一维数字信号) 的处理中，实现去除脉冲噪声和减小扰动，以及在单片机上编程实现和快速运算的方法。
数学形态滤波器通常是用在二维图形的处理，为把数学形态滤波器推广到一维的信号的处理中，下面再介绍一下腐蚀、膨胀、开运算和闭运算这一个基本运算在一维信号处理中的定义：
设H、K分别为h[n]和k[n]的定义域，长度分别为N和M，一般N>M。H和K均为整数集合。
h[n]指包含舵机转角信号的数字化序列，k[n]指结构元素序列。
h被k腐蚀：
h被k膨胀：
h被k开：
h被k闭：
采用数字形态滤波方法，还要选用合理的算法。其中，如何选取模板序列的长度是关键，如果模板序列过长会将有用信号当作噪声滤除，过短则达不到滤除噪声的目的。在采样速率一定的情况下，序列的长度与时间成正比，这要求模板的长度要小于模型车的最小转弯时间，大于舵机扰动的最长时间。第一圈让模型车匀速通过，这样处理有两个优点：
1) 可以固定最小转弯时间，从而确定模板的长度。非匀速通过时速转弯时间不定，要求模板长度可变，从而造成后续处理复杂，稳定性不高。
2) 采样序列的顺序可以直接转化为位移量, 便于后续控制策略处理。相对于非匀速通过速度与时间乘积得到的位移, 直接转化得到的位移更准确 (在标准的韩国赛道上, 实验模型车直接转化得到的赛道长度误差小于5cm, 速度与时间乘积得到赛道长度误差在10cm以上) 。
实验系统在图2所示的赛道上，智能车对赛道信息的采样速率为200Hz，以1.5m/s的速度匀速跑完第一圈的数据如图3所示。可以看到在弯道中，舵机的转角信息存在着严重的毛刺和扰动，不能直接用于第二圈的控制策略。图4为matlab中采用3阶巴特沃兹滤波处理后的结果，干扰的抑制效果仍然不理想，而且运算量偏大，单片机难以承受。图5为采用形态学滤波处理后的数据，赛道信息完整准确，可以较好的应用于后续控制策略。
数学形态滤波的快速算法
由于数学形态滤波器只由加法、减法和比较运算构成，其运算相对简单，因此，它很适合于在计算功能相对较弱的单片机上应用并能取得很好的效果。以往单片机由于受存储容量、计算速度及字长的限制而使大多数的数字滤波器较难实现，而形态滤波器则为单片机应用数字滤波器代替以往的模拟滤波器提供了一条新的途径。
由腐蚀的定义可知，欲计算f (n) 的腐蚀值，需要知道该点前w (w为结构元素的宽度) 点的数据;而要计算f (n) 膨胀后的结果，则需要知道该点后w点的数据。由于运算是一个腐蚀运算接着一个膨胀运算后得到的，在长度为L的数据中只有从第w点到第 (L-w+1) 点，才可以得到开运算的结果。
如图6所示，我们定义一个模板序列，该序列的长度和结构元素的宽度相同。该模板的初始值由前w个点的腐蚀值组成。以第n点为例，沿该点向前的方向对模板序列的值进行膨胀运算，运算的结果即为该点的开运算的结果。同时，沿该点向后的方向继续进行腐蚀运算，得到第 (n+1) 点的腐蚀值。将 (n+1) 点的腐蚀值作为模板序列的最后一个点，并将模板序列前 (w-1) 点顺次向前移动一个位置。更新后的模板值即可用来做 (n+1) 点的膨胀运算，得到在 (n+1) 点的开运算值。如此继续下去，就可完成全部的开运算。在做闭运算也可采用类似的方法来提高计算的速度。
这样，对长度为N的一段数据，用M个零作为其结构元素进行处理时，当采用一般方法进行计算，需要进行2×M×N次的比较运算。而当采用快速算法时，能够在比较最大值的同时得到最小值，减少 (N-M+1) ×M次比较运算，使程序的执行速度提高了近一倍，而且，由于采用这种快速算法，可以实现路径记忆信息的实时处理，很大程度上方便了第二圈控制策略的制定，因此，它使得形态滤波这种方法更加适合应用于路径信息的处理中。
实验及结论
通过对不同传感器方案 (光电管和CCD) 的智能车在不同赛道多次实验发现，对于光电管方案和CCD方案的智能车，赛道记忆算法都能一定程度上提高第二圈速度。智能车采用形态学滤波算法处理赛道记忆数据后，不但行驶的稳定性、准确性有了较大的提升，而且没有大幅增加MCU的资源消耗，同时可以支持复杂的控制策略。上述方案具有很强的通用性，适用于不同传感器方案、不同控制算法的智能车。
参考文献
[1].智能车大赛网站, h t t p://www.smartcar.au.tsinghua.edu.cn/
[2].周斌, 刘旺, 林辛凡等, 智能车赛道记忆算法的研究, 电子产品世界, 2006.08
[3].黄开胜, 金华民, 蒋狄南, 韩国智能模型车技术方案分析, 电子产品世界, 2006.03
[4].Chen Ping, Li Qing-min, Design andanalysis of mathematical morphology-baseddigital filters, Proceedings of the CSEE, v25, n11, June2005, 60-5
数学形态滤波 篇6
在很多情况下, 但大多数情况下的单高斯假设EM算法提出了估计的高斯混合参数的要求, 有限混合是一种灵活而强大的概率建模工具[6,7]。它可以被用来提供一个模型在模式识别领域中的聚类。然而该有限混合到图像分割中的应用遇到一些困难;尤其是高斯噪声[8,9,10]。在本文中, 提出了一种结合形态学滤波的EM图像分割算法。首先利用2G-R-B算法将样本RGB彩色图像转换为灰度图, 再将HSI颜色空间的H分量与其进行融合, 以减少噪声对于分割结果的影响, 然后基于EM算法对融合灰度图像进行图像分割, 最后采用数学形态学滤波及最小面积剔除法进行噪声点的处理。最后通过仿真实验结果表明, 本文提出的算法能够有效的对图像进行分割, 有效增强了算法的自适应性, 改善了分割的效果。
1 分割图像的RGB初始化
本文首先对待分割的图像进行RGB二值化处理, 本文采用了RGB, YUV以及HIS对得到的图像进行颜色空间区分, 对采集到的待分割的图像根据R.G.B分量进行首要分析, 将图像中特征明显的区域与背景区分开。
根据色差分布图, 像素点 (x, y) 的灰度值计算公式定义为:
HSI颜色空间模拟人类的视觉特性, 单个像素的特征被分解为色调 (Hue) 、饱和度 (Saturation) 和亮度 (Intensity) [11], 其中, 色调分量表示物体的颜色信息, 而饱和度分量则可以指示出颜色的深浅。所以本文将G-R-B与HSI颜色空间的色调分量H进行融合, 得到融合后的灰度图像R, 以减少由于噪声的干扰而对分割效果的影响[12]。融合方法如式 (2) 所示。
其中, 为融合参数。
2 EM图像分割算法
EM算法 (Expectation-Maximization algorithm) , 又称为期望最大化算法, 可以在包含了隐含变量的概率模型中寻找参数最大后验概率。在模式识别领域的数据聚类问题中经常用其进行参数的估计。在图像分割中, 基于特征空间聚类的EM算法可以有效规避其在混合成分数及迭代初始值方面的局限性, 可以为分类模型提供有效的极大似然估计结果。
设灰度图像的像素集为{x (1) , ..., x (m) }, 其概率密度函数为:
其中K为分类的数目, p表示概率, θ表示分布的均值μ以及方差 。
本文将图像分为区域部分以及背景两部分, 所以式 (3) 可以写为:
其中, 分别为区域与背景的混合先验概率,
对于要进行分割的灰度图像, 将每个像素样例视为独立, 只要能够找到每个像素隐含的类别z, 即属于图像区域还是属于背景, 能够使得 最大, 的最大似然计算如下:
EM算法的流程是首先初始化分布参数θ, 然后采用迭代的方法重复E步骤及M步骤, 直到收敛。
E步骤:计算隐性变量的后验概率即隐性变量的期望 , 若是第一次计算, 则采用初始化的模型参数θ, 否则采用上一次迭代过程中M步骤更新的θ。
M步骤:求取最大似然函数, 将此时的参数θ作为新的参数传递到E步骤。
EM算法流程图如图1所示。
3 改进的形态滤波EM图像分割
由于采用EM图像分割后, 由于图像中的不同的噪声会对分割的图像造成干扰, 影响了分割的质量, 本文在对分割后的图像采用数学形态学滤波处理, 即采用一种改进的数学形态学滤波及最小面积剔除法进行噪声点的处理。
3.1 开-闭滤波器结构
数学形态学变换一般分为灰度和二值两种变换形式, 通常来说, 二值变换组要用于处理集合, 而灰度的变换主要用于处理函数, 基本的数学形态学变换主要包括了膨胀、腐蚀、闭运算和开运算四种形式。
假设f (x) 和g (x) 分别定义在二维离散空间F和G的两个离散函数, 这其中f (x) 表示的是图像的输入图像, g (x) 表示的是结构元素。那么f (x) 对于g (x) 的膨胀和腐蚀函数分别定义为:
那么f (x) 对于g (x) 的闭运算和开运算分别定义为:
本文即利用了开运算来消除了图像分割中的正脉冲噪声, 利用了闭运算过滤掉了分割过程中的负脉冲噪声。本文为了同时消除分割过程中的正负脉冲噪声, 本文提出了一种采用形态开-闭混合结构形式, 即是构成了一种开闭滤波器。
形态开-闭滤波器设计如下:
3.2 本文算法步骤
根据本文提出的改进的算法思想, 给出如下步骤:
第一步:对采集的RGB图像进行二值化计算, 得到灰度图像;
第二步:将RGB图像转换到HSI图像空间, 提取色调 (H) 分量图, 根据式 (2) 将其与灰度图进行融合, 获得融合后的灰度图像;
第三步:采用EM算法对融合灰度图像进行分割, 取分类数K (28) 3, 根据分割结果得到二值化图像B。
为进一步避免光照的影响, 取阈值0.18与0.26对H分量图进行双阈值分割, 得到二值图像E, 并对图像B的分割结果进行修正, 如式 (13) 。
对于分割出来的图像W, 仍然噪声等问题, 采用公式 (12) 形态学开-闭运算对其进行修正, 去除噪声, 以改善分割的效果
4 实验设计与结果分析
本文的所有实验都是在惠普ENVY 14-U005TX (J6M91PA) 进行的测试, CPU处理器采用的是Intel酷睿i5 4210U, 1.7GHz, 内存大小为4GB, 显卡芯片为NVIDIA Ge Force GTX 850M+Intel GMA HD。仿真软件采用的是mat lab 7.0。
本文实验数据来源分为两种, 一种是合成的图像, 一种是自然图像。通过这两种图像来验证本文算法的有效性能。本文通过这两种类型的实验图像比较了传统的标准的EM算法以及和本文提出的改进EM算法进行了比较, 实验结果如图2和图3所示。
从图2中可以看出, 图像主要包括了3个主要区域, 图2 (b) 显示的是标准的EM图像分割算法分割的图像结果, 图中明显看出, 方法难以抑制噪声对其干扰。图2 (c) 显示的是标准的形态学图像分割算法结果, 相比图2 (b) 对噪声具有了一定的抑制, 但是效果还是不强, 还有一些负脉冲噪声干扰, 图2 (d) 显示的是本文提出的改进算法, 从图中可以直观的看出, 分割效果清晰, 对噪声去除效果明显, 精确的图像轮廓非常清晰。
图3显示的是采集的脑电图MRI分割实验。MRI图像中包括了三层结构, 灰度 (MG) , 白质 (WM) 以及脑脊液 (CSF) , 这三种物质混合在一起, 对分割具有一定的难度。本文对其中增加了5%的高斯噪声进行实验。图3 (b) 显示的是标准的EM算法图像分割结果, 由于5%高斯噪声的影响, 标准的EM无法抑制该噪声。图3 (c) 是形态学图像分割结果, 同样的具有一定的去噪效果, 但是仍然效果不强。图3 (d) 显示的是本文提出的改进算法, 从图中可以直观的看出, 特别是在白质区域和脑脊液区域, 能够有效的对轮廓进行分割。
本文为了定性的分析所提方法的优越性, 本文采用了峰值信噪比 (PSNR) 来对比分割的图像质量。PSNR是衡量图像质量的重要指标, 其实最大的信号量与噪声的强度的比值。
其中, L表示的是图像中像素最大的灰度值;MSE表示的是协方差
所以PSNR值越大, 就代表失真越少, 分割图像的质量越高, 如表1所示, 可以看出本文算法PSNR高。
5 结语
图像分割技术是图像处理和计算机视觉领域一项关键技术, 对其研究一直是数字图像处理技术研究中的热点与难点问题。本文提出了一种结合了数学形态学与EM算法图像分割新方法。方法首先对图像进行RGB二值化处理, 然后通过EM算法对其进行分割, 对分割后的图像采用形态学开-闭组合算法进行去噪处理。最后得到高质量的分割图像, 文章同时也通过峰值信噪比 (PSNR) 定性的分析了算法的性能, 从实验结果可以看出, 不管是直观上还是定性上, 本文提出的改进方法都能够较好地抑制噪声的干扰, 能够有效的保留图像的轮廓, 改善分割的质量。
摘要：由于传统胡EM图像分割算法由于在进行聚类运算时, 耗时太长, 影响了图像分割的精度, 为了解决该问题, 本文提出了一种结合形态学滤波的EM图像分割算法。首先利用G-R-B算法将样本RGB彩色图像转换为灰度图, 再将HSI颜色空间的H分量与其进行融合, 以减少噪声对于分割结果的影响, 然后采用基于EM算法对融合灰度图像进行图像分割, 最后利用一种改进的开-闭混合结构数学形态学滤波对分割后的图像进行噪声点去除, 以减少噪声对图像质量的影响。最后通过仿真实验结果表明, 本文提出的算法能够有效的对图像进行分割, 有效增强了算法的自适应性, 改善了分割的效果。
关键词：图像分割,EM算法,数学形态学,滤波技术
参考文献
[1]Fu J C, Chen C C, Chai J W, et al.Image segmentation by EM-based adaptive pulse coupled neural networks in brain magnetic resonance imaging[J].Computerized Medical Imaging and Graphics, 2010, 34 (4) :308-320.
[2]Lu Q, Huang X, Zhang L.A Novel Clustering-Based Feature Representation for the Classification of Hyperspectral Imagery[J].Remote Sensing, 2014, 6 (6) :5732-5753.
[3]李宗民, 公绪超, 刘玉杰.多特征联合建模的视频对象分割技术研究[J].计算机学报, 2013, 36 (11) :2356-2363.
[4]Cardoso M J, Clarkson M J, Ridgway G R, et al.Lo Ad:a locally adaptive cortical segmentation algorithm[J].Neuro Image, 2011, 56 (3) :1386-1397.
[5]Dalmiya S, Dasgupta A, Datta S K.Application of Wavelet based K-means Algorithm in Mammogram Segmentation[J].International Journal of Computer Applications, 2012, 52 (15) :15-19.
[6]杨名宇, 丁欢, 赵博, 等.结合邻域信息的Chan—Vese模型图像分割[J].计算机辅助设计与图形学报, 2011, 23 (3) :413-418.
[7]Velasco-Forero S, Angulo J.Supervised ordering in:Application to morphological processing of hyperspectral images[J].Image Processing, IEEE Transactions on, 2011, 20 (11) :3301-3308.
[8]许新征, 丁世飞, 史忠植, 等.图像分割的新理论和新方法[J].电子学报, 2010, 38 (2) :76-82.
[9]Yao H, Duan Q, Li D, et al.An improved K-means clustering algorithm for fish image segmentation[J].Mathematical and Computer Modelling, 2013, 58 (3) :790-798.
[10]王芳梅, 范虹.利用改进CV模型连续水平集算法的核磁共振乳腺图像分割[J].西安交通大学学报, 2014, 48 (2) :38-43.
[11]Lucchi A, Smith K, Achanta R, et al.Supervoxel-based segmentation of mitochondria in em image stacks with learned shape features[J].Medical Imaging, IEEE Transactions on, 2012, 31 (2) :474-486.
数学形态滤波 篇7
滚动轴承是旋转机械中应用最广的重要部件之一，它对整个旋转机械的工作状态有着重要影响，一旦发生故障可能导致严重后果，因此，滚动轴承的状态监测和故障诊断有着十分重要的意义[1,2]。国内外学者采用自回归模型、神经网络、支持向量机、Hilbert-Huang变换、边际谱等方法成功地对滚动轴承进行了故障诊断[3,4,5]。
在工程实际中，对滚动轴承进行振动信号采集时会受到严重的噪声干扰，尤其是在故障发生的早期。要获得准确的诊断结果，首先要抑制信号的噪声，使有用的故障特征信息更加突出。形态学滤波是在数学形态学变换的基础上发展起来的一种重要的非线性滤波工具[6,7]，利用该滤波器对信号进行滤波时，能够根据待分析信号的局部形状特征，将其分解为具有物理意义的多个部分，并将其与背景剥离，同时保留信号的主要形状特征[8,9]。目前该方法已成功地应用到对齿轮等故障信号的滤波处理中[10]。
1948年Shannon[11]把热力学中熵的概念引入信息论中，定义了信息熵。信息熵可以定量地描述信息的不确定性和复杂性。目前，信息熵理论已经成功地应用于压缩机气阀、电机和发动机的故障诊断。
本文融合数学形态学滤波的高效消噪能力和信息熵理论在信号表征方面的优越性，提出了一种新的滚动轴承故障诊断方法。
1 数学形态学和差分熵
1.1 数学形态学基本变换
数学形态滤波的核心思想是利用起到滤波器作用的结构元素对待分析信号进行形态学迭代运算，从而实现滤波[12]。四种基本的形态学算子分别为腐蚀、膨胀、开和闭。
若f(n）为一维原始离散信号，其定义域为F={0,1，…，N-1}，g(n）为一维离散信号，称其为结构元素，其定义域为G={0,1，…，M-1}，且MN，则f(n）关于g(n）的腐蚀和膨胀运算分别定义如下：
m∈{0,1，…，M-1}
f(n）关于g(n）的开和闭运算分别定义如下：
1.2 差值滤波器
一维信号f(n）分别经过结构元素g(m）膨胀和腐蚀后的差值称为差值滤波，其表达式为
式（5）可以分解为
其中，-f和f-fg正是数学形态学中Top-Hat变换中的黑白帽变换，而这两种变换分别能提取信号中的负正脉冲，所以差值滤波器将两种变换融合，能够更好地对信号进行滤波处理，提取信号中的正负脉冲信息，克服了开闭算子在处理实际信号中需要得到正负脉冲的先验知识的缺点[13]。
在本文中，将利用该滤波器对轴承故障信号进行滤波处理。
1.3 结构元素
形态学滤波的实质就是通过结构元素与待分析信号进行形态学迭代运算，从而实现滤波目的。所以结构元素的长度对形态学滤波的效果有着重要的影响。
常用的结构元素有扁平型、三角型和半圆型等。三角型和半圆型结构元素具有长度和高度两个参数，分别适用于脉冲噪声和随机噪声的滤除。而只有长度参数的扁平型结构元素具有计算简单、需要优化的参数少等优点，同时因其高度为零，可避免对信号幅值的改变，所以获得了广泛的应用[10]，因此在本文中，选用优化参数少且计算简单的扁平型结构元素来对轴承故障信号进行滤波处理。
1.4 差分熵
差分熵（difference entropy,DE）能够凸显信号突变点的信息，设fDIF(n）为滤波器处理后的故障信号，则差分熵的定义如下[14]:
n=0,1，…，N-2
差分熵为
为减小计算量和提高诊断效率，信号长度L不宜过大。
2 数值仿真算例
2.1 仿真信号
为了验证差值滤波器和差分熵两种方法的结合对故障诊断的效果，设计了以下仿真信号进行分析：
其中，x1(t）是频率为16Hz的周期性指数衰减冲击信号，每周期内衰减函数为8e-100tsin(512πt），周期冲击时间间隔为0.0625s，用来模拟故障冲击信号；x2(t）是周期为20Hz和40Hz的叠加低频谐波干扰信号：cos(40πt)+cos(80πt);x3(t）为标准差为1的高斯白噪声，用于模拟强背景噪声。采样频率为2048Hz，采样时间为1s。
2.2 仿真结果分析
图1为模拟故障信号的时域图和0～100Hz低频段功率谱图。由于高斯白噪声和低频谐波的干扰，从图1a中难以分辨出冲击成分；从图1b中可以清楚地看出20Hz和40Hz的低频谐波干扰，而故障频率16Hz及其倍频处的谱峰已经被背景噪声所淹没。
为消除噪声的干扰，提高故障诊断的正确率，根据文献[15]提出的滤波最优效果原则：特征频率处的能量占选定频率段的总能量最大时，滤波效果最好。基于这一原则，选取长度为9的扁平型结构元素的差值算子对含有噪声的模拟故障信号进行滤波处理。图2所示为滤波后的结果，图2a展示了清晰可见的冲击成分，图2b展示了故障频率16Hz及其倍频处的谱峰，峰值也较图1b有很大提高。20Hz和40Hz处的低频谐波干扰和高斯白噪声得到了有效的抑制，说明差值滤波器具有很好的滤波效果。
对滤波后的模拟故障信号进行差分熵的求取（L选定为2048），结果如图3a所示。为了更加准确、高效地进行故障诊断，对差分熵的结果进行了阈值处理，即将小于最大值0.4倍的熵值置为0，如图3b所示。
理论上，在1s的时间内有16次的冲击，从图3b中可以很明显地看出在1s时间内有16束熵值线，初步可以诊断出频率为16Hz的故障。相邻束熵值线间突变平均时间间隔为62.375ms，和理论周期冲击时间间隔0.0625s相差0.2%。由此可以看出误差值很小，说明形态差值滤波和差分熵理论的结合可以有效地诊断出该模拟故障。
3 轴承故障信号分析
为了验证该方法的有效性，本文采用美国Case Western Reserve University提供的两种故障程度的滚动轴承内圈和外圈坑点损伤的故障数据，进行进一步实验分析验证[16]。
该轴承的型号为SKF6250，电机转速为1750r/min，载荷为1.497kW，转频约为29Hz。实验中内圈和外圈两处坑点损伤是由电火花机人工加工制作的，轻微故障程度的损伤直径为0.0178 mm，重度故障程度的损伤直径为5.334mm。内外圈故障的特征频率分别为156Hz和103.4Hz[13]。用加速度传感器对电机驱动端进行数据采集，采样频率为12kHz，采样长度为0.5s。
3.1 轴承内圈故障分析
当滚动轴承出现内圈故障时，故障频率为156Hz，周期冲击时间间隔为6.41ms。
滚动轴承具有径向间隙，承受单边载荷，当故障发生时，振动幅值会随着损伤部分与滚动体发生冲击接触位置的不同而发生周期性变化，即发生幅值调制。按轴旋转频率fr进行调制的振动频率为nZfi±fr(n=1,2，…），以滚动体的公转频率fc进行调制的振动频率为nZfi±fc，其中，Z为滚动体个数，fi为内圈故障频率。
为了更好地显示和观察差分熵故障诊断方法及其优势，截取0.5s的内圈早期故障信号长度中的0.08s数据进行分析。图4是长度为0.08s含有噪声的轴承内圈早期故障信号的时域图和0～1kHz低频段功率谱图。由于噪声的干扰，从图4a中难以识别冲击成分；图4b展示了在故障频率156Hz处存在比较模糊的谱峰，但其倍频处的谱峰被噪声淹没，同时在600Hz处的干扰相对非常大。
为滤除噪声的干扰，提高故障诊断的准确性，根据文献[15]提出的最优滤波效果原则，选取长度为9的扁平结构元素的差值算子对轴承内圈的早期故障信号进行滤波处理，结果如图5所示。
图5a展示出噪声得到了很好的抑制，冲击成分很明显；图5b展示在故障频率156Hz及其二倍频处存在谱峰，峰值较图4b有所提高，能够清晰地反映轴承内圈的故障特征，说明差值滤波器对轴承内圈故障信号具有很好的滤波效果。
对滤波后的轴承内圈早期故障信号进行差分熵的求取（L选定为6000），结果如图6a所示。为了更加准确、有效地进行故障诊断，对差分熵的结果进行阈值处理，即将小于0.4倍最大值的熵值置为0，结果如图6b所示。
理论上，在0.08s时间段内包含12.48次冲击，从图6b中可以很明显地看出在0.08s时间段内有13束熵值线，所以可判断为轴承内圈故障。突变相邻束熵值线间平均时间间隔为6.2838ms，和理论周期冲击时间间隔6.41s相差2.68%。由此可以看出，相邻冲击平均时间间隔误差值非常小，说明形态差值滤波和差分熵理论的结合可以准确地诊断出该故障为轴承内圈的轻度故障。
背景噪声的增强往往伴随着故障程度的劣化，为了考察该方法对轴承内圈重度故障和强背景噪声加大时的有效性，截取0.5s的轴承内圈重度故障信号长度中的0.08s的数据进行分析。
图7是长度为0.08s、含有噪声的轴承内圈重度故障信号的时域图和0～1kHz低频段功率谱图。图7b展示了在故障频率156Hz处存在比较模糊的谱峰，其峰值相对很小，其倍频处的谱峰也被噪声淹没，同时在444Hz和602.5Hz处的干扰谱峰相对非常大。
为滤除噪声的干扰，提高故障诊断的准确性，根据文献[15]提出的最优滤波效果原则，选取长度为15的扁平结构元素的差值算子对轴承内圈的重度故障信号进行滤波处理，结果如图8所示。
对滤波后的轴承内圈重度损伤信号进行差分熵的求取（L选定为6000），结果如图9a所示。为了更加准确、有效地进行故障诊断，阈值的选取遵循以下原则：相邻束线间的时间间隔能够最大程度地符合轴承内圈故障的理论周期冲击时间间隔，所以将小于0.4倍最大值的熵值置为0，结果如图9b所示。
理论上，在0.08s时间段内包含12.48次冲击，而图9b所展示的有些相邻束线间的时间间隔和理论周期冲击6.41ms的时间间隔误差不是很大，但是有些时间间隔很不规律，误差也和理论周期冲击的时间间隔相差很大，由此可得：该方法对轴承内圈重度故障的诊断效果不是很理想。
3.2 轴承外圈故障分析
当滚动轴承出现外圈故障时，故障频率为103.4Hz，周期冲击时间间隔为9.671ms。
当故障发生，在滚动体通过时，会产生冲击振动，由于外圈损伤点以频率Zfo不断产生，所以引起一系列的高频衰减振动，由于阻尼较小，因此脉冲衰减时间比1/（Zfo）要小得多，衰减振动也基本独立，振动频率为Zfo。其中，fo为外圈通过频率。
为了更好地显示和观察差分熵故障诊断方法及其优势，截取0.5s的外圈早期故障信号长度中的0.08s数据进行分析。
图10是长度为0.08s含有噪声的轴承外圈的早期故障信号的时域图和0～1kHz低频段功率谱图。从图10b可以看出，由于噪声的干扰，在故障频率103.4Hz及其倍频处的故障特征被噪声淹没，而在故障频率720Hz处的干扰非常大。
为消除噪声的干扰，提高故障诊断的准确性，根据文献[15]提出的最优滤波效果原则，选取长度为5的扁平结构元素的差值算子对轴承外圈早期故障信号进行滤波处理，结果如图11所示。
图11a展示出噪声得到了很好的抑制，冲击成分很明显；图11b中展示在故障频率103.4Hz及其倍频处存在谱峰，峰值较图10b有很大提高，能够清晰地反映轴承外圈的故障特征，说明差值滤波器对轴承外圈故障信号具有很好的滤波效果。
对滤波后的轴承外圈早期故障信号进行差分熵的求取（L选定为6000），结果如图12a所示。为了更加准确、有效地进行故障诊断，对差分熵的结果进行了阈值处理，即将小于0.4倍最大值的熵值置为0，结果如图12b所示。
理论上，在0.08s时间段内包含8.27次冲击，从图12b中可以很明显地看出在0.08s时间段内有8束熵值线，所以可判断为轴承外圈故障。相邻束线间突变平均时间间隔为6.299ms，和理论周期冲击时间间隔相差0.96%。由此可以看出相邻冲击平均时间间隔误差值非常小，说明形态差值滤波和差分熵理论的结合可以准确地诊断出该故障为轴承外圈的轻度故障。
背景噪声的增强往往伴随着故障程度的劣化，为了考察该方法对轴承外圈重度故障和强背景噪声加大时的有效性，截取0.5s的轴承外圈重度故障信号长度中的0.08s的数据进行分析。
图13是长度为0.08s含有噪声的轴承外圈重度故障信号的时域图和0～1kHz低频段功率谱图。由图13b可以看出，由于噪声的干扰，在故障频率103.4Hz及其倍频处的故障特征被噪声淹没，并且在故障频率628Hz处的干扰非常大。
为消除噪声的干扰，提高故障诊断的准确性，根据文献[15]提出的最优滤波效果原则，选取长度为16的扁平结构元素的差值算子对轴承外圈重度故障信号进行滤波处理，结果如图14所示。
图14a展示出噪声得到了很好的抑制，冲击成分很明显；图14b展示了在故障频率103.4Hz及其二倍频处存在谱峰，峰值较图13b有很大提高，能够清晰地反映轴承外圈的故障特征，说明差值滤波器对轴承外圈故障信号具有很好的滤波效果。
对滤波后的轴承外圈重度故障信号进行差分熵的求取（L选定为6000），结果如图15a所示。为了更加准确、有效地进行故障诊断，阈值的选取遵循以下原则：相邻束线间的时间间隔能够最大程度地符合外圈故障的理论周期冲击时间间隔，所以将小于0.25倍最大值的熵值置为0，结果如图15b所示。
理论上，在0.08s时间段内包含8.27次冲击，而图15b所展示的相邻束线间的时间间隔很不规律，误差也和理论周期冲击9.671ms的时间间隔相差很大，由此可得：该方法对轴承外圈重度故障的诊断效果不太理想。
综上所述，本方法较神经网络无需用大量样本对模型进行训练，只需一个样本；较贝叶斯网络无需大量且复杂的先验知识，只需可能出现的故障频率；较支持向量机无需用其他智能方法对参数进行优化，只需简单地设置一个阈值；此外，本方法涉及不到知识冗余，从而避免了降维处理。所以，本方法非常适用于对滚动轴承早期故障的诊断，为早期故障诊断提供了一个很好的工具。其不足之处是对晚期故障（即伴随着强噪声背景的重度故障）的诊断效果不是很理想，有待进一步深入研究。
4 结论
(1）形态差值滤波器能够很好地滤除轴承内外圈故障信号的噪声，提取出更多的有用故障特征信息。
数学形态滤波 篇8
红外弱小目标的检测一直是红外监视告警系统中的关键技术之一。探测的红外目标因成像距离远、面积小而缺乏形状和结构特征,且由于复杂的红外图像背景而呈现低信噪比、低对比度等特点。若想要很好地检测出此类红外弱小目标,则必须对图像进行预处理,最有效的办法就是对复杂的红外背景进行大幅抑制,然后进行阈值分割检测出目标位置,所以能否高性能地抑制背景决定了后续目标虚检、漏检的概率。
目前,红外图像的背景抑制技术取得了很大的发展,常用的方法包括空域和频域滤波、小波域滤波、SVD滤波、数学形态滤波等方法[1,2,3,4,5,6,7,8],各种方法都有自己的特点,能实现一定的背景抑制。一般的空域和频域滤波能够实现简单的背景抑制,但对红外图像中背景起伏比较大且目标局部信噪比较低时,目标可能当作背景被抑制掉,或者不能完全平滑掉背景边缘;采用小波域等多尺度多分辨率滤波能够得到很好的结果,但算法相对复杂,不能取得很好的时效性;单一的SVD滤波和数学形态滤波等会在背景抑制后使目标的强度被弱化,使得目标与残留背景的对比度变低,后续目标检测时选的目标点过多,误检的机率变大。
因此本文结合红外图像中背景、目标和噪声的特点,提出了一种将改进的SVD和Tophat相结合的红外图像背景抑制算法。首先通过奇异值分解对红外图像进行非线性增强,以提高目标与目标周围背景的对比度,然后再利用形态学滤波中的Tophat变换进行滤波达到背景抑制的目的。该方法不仅能抑制强起伏背景,而且能够增强目标信号,从而达到较好的检测数值指标和视觉效果。
1 基于奇异值分解的红外图像增强
1.1 奇异值分解[9]
设数值矩阵A是M × N的复矩阵,且A的秩为r≤min(M,N),则存在正交(或酉)矩阵UM×M和矩阵VN × N,使得A的SVD可表示为:
式中:,是r阶对角矩阵,Σ = diag(σ1,σ2,…,σr),σi>0(i=1,2,…,r)为矩阵A的非零奇异值;UM × M=[U1,U2,···,UM] 称为A的左奇异矩阵,由M阶列向Ui=[U1i,U2i,···,UMi]T, i = 1,2,⋯,M构成;VN × N=[V1,V2,…,VN]称为A的右奇异矩阵,由N阶列向Vi=[V1i,V2i,…,VNi]T,i = 1,2,…,N构成;AAT的正交单位特征向量组成U,特征值组成 ΓΤΓ,AΤA的正交单位特征向量组成V,特征值与AAT相同,因此式(1)可以改写为:
由式(2)可以得出:用r个非零奇异值对应的r个分量的线性和可实现对矩阵A的重构。
由奇异值分解的性质可以得知,如果对二维红外图像进行SVD分解,则奇异值矩阵中奇异值的大小反映了图像能量的分布。如图1 所示,图1(a)是一幅大小为128×128的真实云天背景的红外场景图像,其中添加信噪比不大于2 的弱小目标;图1(b)为图像对应的奇异值曲线图。通过图1(b)可以得知,图像矩阵的奇异值在前10~20个序号之间急剧下降,后续几乎为零,所以图像的大量信息主要集中在前面十几个较大的奇异值中。
根据图1(b)的奇异值曲线图分别采用前15 个和前100 个奇异值重构的图像,如图1(c)和(e)所示,图1(d)和(f)分别为对应的灰度三维图。从图中信息可以看出,当采用奇异值矩阵中前15 个奇异值重构图像时,重构后的图像就已经携带了大量的图像信息,主要表现为图像的背景等低频分量;当重构图像采用的奇异值增加到100 个时,重构图像变得清晰,增添了图像的细节信息,也就是对应图像的高频分量。因此,选取不同数量的奇异值对图像进行重构,可以得到不同的结果。
1.2 改进奇异值分解的红外图像增强
将奇异值分解用在红外图像背景抑制方面,主要是通过选取有效的奇异值来重构图像的低频背景信息,然后用原始图像减去重构图像就可以实现背景的抑制。但此种方法存在目标强度变弱,对高频信息抑制有限等情况,特别是在复杂红外背景中更是不利于后续目标的分割。因此,本文提出先将红外图像进行增强,然后再背景抑制。
图像增强的方法有很多,本文根据红外图像中目标、背景和杂波的关系,背景主要是低频成分,目标和杂波及其图像轮廓等是高频部分,目的是要增强目标的强度,就是要对图像的高频分量进行增强。而将原图像进行奇异值分解后,根据奇异值分解的性质,目标和杂波及其图像轮廓等高频部分信息主要集中在后续小的奇异值上,因而采用对奇异值进行对数变换的方式来实现。
对数变换是一种常用的非线性变换,对图像奇异值进行对数变换,对数曲线在奇异值比较低的地方斜率大,奇异值较大的地方斜率比较低,因而通过对数变换提升了红外图像在较暗区域的对比度,因而能增强显示出暗部的细节,同时也增强目标信号的能量。对图像奇异值采用对数变换的形式可表达为:
式中:σi是图像的奇异值;λi为变换后的奇异值;k为常数,可以根据输入图像适当调整。由于k值相当于是对非线性变换后再作线性增强,所以当目标灰度值较低时为进一步增大目标能量可以选取大的k值。由于k值的增大也会增强其他的高频分量,所以会影响后续的背景抑制程度,根据本文多次实验比较,在强起伏背景下建议选取3 左右。图2 给出了图1 采用k= 5 和k = 10时的增强结果。
2 改进的奇异值分解和形态滤波相结合的弱小目标背景抑制
2.1 形态学Tophat滤波的背景抑制
Tophat变换称为高帽变换,它是灰度形态学的一种处理方法。该变换使用上部平坦的柱体或平行六面体(像一顶高帽)作为结构元素。图像的高帽变换可以记为:
它是将原图像与进行形态学开运算后的图像相减得到的残差图像[10]。开运算指的是先对图像进行腐蚀然后再膨胀:
通过开运算可以去掉图像上那些与结构元素的形态相吻合的“高帽”结构。通过SVD非线性增强后的红外弱小目标图像的目标点已经远比周围的背景像素要亮,因此可以认为目标点是红外图像的“高帽”,因而通过开运算实现了红外图像背景的预测,再用原图像减去开运算后的图像就可以得到背景抑制后的图像。
2.2 基于改进的奇异值分解和形态滤波的弱小目标背景抑制算法的实现
由于采用对数非线性增强后的红外图像对背景的轮廓等高频分量同时也进行了增强,所以如果继续选取部分奇异值进行重构来抑制背景,抑制后的背景会保留很多轮廓细节,因而本文采用改进的奇异值分解和形态滤波Tophat相结合的弱小目标背景抑制算法进行背景抑制。
本文算法具体实现步骤如下:
(1)读入包含红外弱小目标图像A;
(2)对图像A的数值矩阵进行奇异值分解,[U,Σ,V]= SVD(A), 得到奇异值矩阵:Σ =diag(σ1,σ2,…,σr),σi>0(i=1,2,…,r);
(3)对奇异值σi采用式(3)进行对数变换,得到优化后的奇异值矩阵;
(4)利用 λi构成的奇异值矩阵重构增强后的红外图像;
(5)对增强后的红外图像作Tophat变换,最终得到背景抑制后的目标图像。
3 实验结果与分析
采用128 pixel×128 pixel的云天背景的红外图像,256 灰度级,信杂比为1.5 左右,对比度约为8%。给出了两组采用本文方法进行背景抑制处理的结果,并将其与SVD和Tophat滤波的抑制结果进行了比较。
3.1 参数指标
图像背景抑制结果的好坏通常采用以下几个参数来衡量。它们分别是:
图像信杂比:SCR = |Gt- Gb|/σbc
图像信杂比增益:ISCR = SCRout/SCRin
图像对比度:CR=|Gt-Gb|/|Gt+Gb|×100%
图像对比度增益:ICR=CRout/CRin
图像背景抑制因子:BSF =σin/σout
式中:Gt,Gb代表目标的灰度均值和目标周围一定区域内的灰度均值;σbc是背景的均方差;σin和 σout分别为输入图像的均方差和输出图像的均方差。
3.2 实验结果
为了比较本文算法的优势,特选取两幅不同信杂比和不同对比度的弱小目标图片进行验证。实验结果如图3,图4 所示。两幅图中的目标仅占几个像元,呈现为点状,没有任何结构与形状特征,且图片背景为天空背景,具有强烈起伏的云层,所以两幅图片的性噪比都不高。第一幅图片的目标点处于白云之上,和周围背景相比,强度略高于背景。第二幅图片的目标恰好处于云层的边缘,目标给强烈起伏的云层干扰了,造成目标及周围背景的对比度很低。
将两幅图片分别采用SVD方法和Tophat变换进行实验,然后与本文的方法进行对比,本文奇异值图像增强选择参数k=3,图3(b1)和图4(b1)为SVD方法进行背景抑制,采用文献[8]中提到的偏差指数来选择重构的奇异值个数。图3(c1)和图4(c1)为直接将图像用Tophat变换进行处理的结果。图3(d1)和图4(d1)为本文所提方法处理的结果。为了更加直观地观察处理结果,给出了每种方法处理的结果对应的三维视图,如图3(b2),图3(c2),图3(d2)和图4(b2),图4(c2),图4(d2)所示。评价指标参数值如表1 所示。
3.3 结果分析
由实验结果可以看出,经SVD方法处理后,图像背景能够得到一定的抑制,但目标的强度明显减弱,如果进一步加大抑制,将会使目标变得更加的微弱,不利于后续目标的分割。采用Tophat变换直接处理虽然也能进行背景抑制,但在目标与背景对比度不高的情况下会有一定的误差。从第二幅图用Tophat变换处理的三维结果来看,处于图片背景的边缘也留有一定的强度,它会严重影响后续目标的检测,并且从大量的仿真实验验证得出,Tophat变换只针对信噪比较高的场景有良好的效果。而采用本文的方法处理都能很好地对复杂背景进行抑制,不但云层被平滑掉,连云层的边界也能很好的平滑掉,使得图像的信杂比和对比度得到很大改善,并且用表1 中的评价指标参数值也可以推断出本文的算法具有一定的实用价值。
4 结论
针对复杂背景中红外弱小目标检测的背景抑制问题,利用奇异值分解的方法实现对红外图像的非线性增强,然后将其与Tophat变换相结合,提出了一种新的红外弱小目标背景抑制算法。将算法利用真实红外场景序列图像进行验证,结果证明了该方法抑制背景的有效性。该算法在对比度增益和信杂比增益以及背景抑制程度上都明显优于单一的SVD和Tophat变换,因而可以将该算法作为红外弱小目标检测的理论依据之一,具有一定的实用性。
摘要：红外弱小目标的复杂背景抑制一直是弱小目标检测与跟踪的一个难点。提出一种改进的奇异值分解和形态滤波Tophat变换相结合的红外弱小目标背景抑制算法。首先通过奇异值分解得到原红外图像的奇异值矩阵和左右奇异矩阵,然后通过对奇异值进行对数非线性变换,利用优化后的奇异值矩阵进行重构得到增强对比度的红外图像,最后利用形态滤波中的Tophat变换进行滤波达到背景抑制的目的。实验结果表明,该算法能够很好地实现红外弱小目标图像的背景抑制,并能使目标信号得到保存和增强。
关键词：红外图像,目标检测,背景抑制,奇异值分解,形态滤波
参考文献
[1]彭嘉雄,周文琳.红外背景抑制与小目标分割检测[J].电子学报,1999,27(12):47-51.
[2]秦翰林,刘上乾,周慧鑫,等.采用Gabor核非局部均值的弱小目标背景抑制[J].红外与激光工程,2009,38(4):737-741.
[3]张媛,辛云宏,张春琴.基于时空联合滤波技术的缓慢运动红外弱小目标检测算法[J].光子学报,2010,30(10):2049-2054.
[4]马文伟,赵永强,张国华,等.基于多结构元素形态滤波与自适应阈值分割相结合的红外弱小目标检测[J].光子学报,2008,37(1):1020-1024.
[5]ERCELEBI E,KOC S.Lifting-based wavelet domain adaptive wiener filter for image enhancement[J].IEE Proceedings of Vision,Image and Signal Processing,2006,153(1):31-36.
[6]SELESNICK I W,BARANIUK R G,KINGSBURY N G.The dual-tree complex wavelet transform[J].IEEE Signal Processing Magazine,2005,22(6):123-151.
[7]ZHANG Biyin,ZHANG Tianxu,CAO Zhiguo,et al.Fast new small target detection algorithm based on a modified partial differential equation in infrared clutter[J].Optical Engineering,2007,46(10):1-6.
[8]秦翰林,周慧鑫,刘上乾,等.基于奇异值分解的红外弱小目标背景抑制[J].半导体光电,2009,30(3):473-476.
[9]张贤达.矩阵分析与应用[M].北京:清华大学出版社,2004:341-400.