多尺度形态梯度和标记分水岭的时频谱图分割

郭海涛，徐 雷，赵红叶，田原嫄，焦圣喜

(1.内蒙古大学电子信息工程学院，内蒙古 呼和浩特 010021；2.东北电力大学电气工程学院，吉林 吉林 132012；3.东北电力大学自动化工程学院，吉林 吉林 132012；4.东北电力大学机械工程学院，吉林 吉林 132012)



多尺度形态梯度和标记分水岭的时频谱图分割

郭海涛1，2，徐 雷3，赵红叶3，田原嫄4，焦圣喜3

(1.内蒙古大学电子信息工程学院，内蒙古 呼和浩特 010021；2.东北电力大学电气工程学院，吉林 吉林 132012；3.东北电力大学自动化工程学院，吉林 吉林 132012；4.东北电力大学机械工程学院，吉林 吉林 132012)

时频谱图干扰强，目标之间、目标与干扰之间有重叠，其分割是重要而困难的问题.提出一种基于图像熵定义的时频谱图多尺度形态梯度图像融合方法，将该方法与标记分水岭分割结合形成一种基于多尺度形态梯度和标记分水岭的时频谱图分割方法.实验结果表明，与基于单尺度形态梯度和标记分水岭的分割方法相比，新方法实用性更强;与Otsu法相比，新方法分割更准确.

时频谱图;多尺度梯度图像;图像熵;水岭算法;图像分割

众所周知，利用时频分析方法可以获得时间信号(包括噪声)的时频分布.将该时频分布作为图像，即为信号的时频谱图.时频谱图的处理有助于对信号进行时频分析，而信号的时频分析有广阔的应用前景.信号的种类(应用的领域)不同，时频分析的方法也可以不同，因此时频谱图的种类多种多样，但其处理本质上属于图像处理范畴.时频谱图分割的目的是分离出时频曲线.显然，时频谱图的分割是时频谱图低层处理的核心与关键.

时频谱图分割是重要而困难的问题，有关其分割研究的报道甚少.不同领域的学者以不同的应用为背景研究这个问题.文献[1]利用Otsu阈值化方法分割柴油机振动信号时频谱图，文献[2]利用Otsu阈值化方法结合区域分裂与增长实现心音信号时频谱图的分割.文献[1-2]中的方法适用于直方图为理想双峰的时频谱图的分割，而对于直方图为非理想双峰的时频谱图的分割效果较差.文献[3]利用腐蚀、膨胀、开闭等数学形态学运算对二值化后的心音信号时频谱图进行分割，其方法不适用于灰度时频谱图.文献[4]利用BP神经网络分割时频谱图，该方法存在以下原理上缺欠:(1)需要大量的学习样本;(2)样本的选择无理论指导;(3)方法的性能与多种因素有关，存在不确定性.文献[5]利用分水岭算法分割海豚音等自然信号的Capon时频谱图，文献[6]利用分水岭算法分割语音信号时频谱图，文献[7]应用文献[6]的方法分割步态周期时频谱图.文献[5-7]中的方法容易得到连续的单像素宽边界，但其梯度算子容易受到干扰和量化误差的影响而产生局部极小值，从而导致过分割.为了解决过分割，文献[8-9]将分水岭算法与区域合并、区域增长结合使用，取得较好的效果.这些方法均是基于分水岭算法的，笔者因此受到启示，考虑利用分水岭算法分割时频谱图，但使用的分水岭方法是基于多尺度形态梯度的，这一点与前述方法不同.

考虑到噪声和过多伪极小值点的存在可导致分水岭算法时频谱图严重的过分割，因此先对时频谱图进行滤波.这里利用文献[10]中的形态学开闭重构滤波器.分水岭算法图像分割的性能在很大程度上依赖于待分割图像的梯度图像[11]，因此梯度图像的形成是分水岭算法图像分割中的重要环节.考虑到形态梯度图像使图像灰度级阶跃更为急剧、图像边缘更加突出，多尺度梯度图像能够更好地保留图像边缘细节，因此在时频谱图分割中利用多尺度形态梯度图像.为了形成时频谱图的多尺度形态梯度图像，笔者提出一种基于图像熵定义的时频谱图多尺度形态梯度图像融合方法.上述时频谱图的滤波和梯度图像的形成统称为时频谱图的预处理.

1 开闭重构滤波器

利用形态学开闭重构滤波器[10]对时频谱图进行平滑，抑制图像中的噪声及伪极小值点.形态学开闭重构滤波器的步骤如图1所示.首先对图像进行测地腐蚀得到标记图像，然后以原始图像作为掩膜图像对标记图像进行形态学开重构.在此基础上，对开重构后的图像进行测地膨胀，将膨胀得到的结果作为掩膜图像，对开重构后的图像进行形态学闭重构.

图1 开闭重构滤波器

图1中相关操作的定义如下:若f为掩膜图像，b为结构元素，r为标记图像，则测地膨胀和测地腐蚀运算可分别表示为

dr(f)=min{d(f)，r}，er(f)=max{e(f)，r}，

其中d(f)和e(f)分别表示利用结构元素b对图像f作膨胀和腐蚀运算.相应的膨胀重构和腐蚀重构分别表示为

膨胀重构和腐蚀重构均为迭代运算，当运算收敛时迭代终止.开重构与闭重构分别表示为

2 多尺度形态梯度图像融合方法

基于分水岭变换的图像分割方法，其性能在很大程度上依赖于形成梯度图像的梯度算法[11].传统的梯度算子，如高斯滤波器一阶偏微分和形态学梯度算子存在严重的缺陷，它们会受噪声和量化误差的影响，在均匀一致的区域内部产生过多的局部最小值，从而导致过分割现象.针对时频谱图干扰强的特点，笔者提出一种基于图像熵定义的时频谱图多尺度形态梯度图像融合方法.

文献[12-13]定义单尺度形态梯度算子为

Grad(f)=(f⊕b)-(fΘb)，

(1)

其中⊕和Θ分别为膨胀和腐蚀运算，其性能取决于结构元素b的大小.小尺度的结构元素去噪能力弱，但能有效地检测梯度边缘细节;大尺度的结构元素去噪能力强，但梯度图像边缘细节有所丢失.为了更好地利用小尺度结构元素和大尺度结构元素的各自优点，近年来学者们提出平均加权的多尺度形态梯度算子[11].这种简单的取均值进行多尺度形态梯度图像融合的方法不能体现出不同尺度梯度图像的差别，因此不能很好地反映不同尺度下梯度图像的综合信息.结合图像熵的定义，笔者提出一种新的多尺度形态梯度算子.该算子的实现步骤是:分别求取单尺度形态梯度图像，利用图像熵的定义计算各尺度形态梯度图像的加权值，然后通过加权求和将单尺度形态梯度图像融合为多尺度形态梯度图像.多尺度形态梯度图像表示为

(2)

利用(2)式即可融合得到时频谱图多尺度形态梯度图像.与文献[11]不同，这里选择圆盘形结构元素bj.取3个不同的尺度，即取序号j =1，2，3，相应的圆盘形结构元素bj的半径分别为1，2，3.对应于这3个尺度下的形态梯度图像的熵值和加权值分别为E1，E2，E3和a1，a2，a3.

3 标记分水岭分割方法

3.1 传统分水岭分割

分水岭算法是基于拓扑理论的数学形态学分割方法，其基本思想是将图像看作测地学上的拓扑地貌.分水岭算法有很多种，具有代表性的是基于浸没模拟(Immersion Simulation)[13]的过程.该方法将灰度图像看成地形学上的地貌，是一种自下而上形成区域的方法.文献[13]给出的基于浸没模拟的算法级定义是较为严格的递归定义，具有一定的代表性.该算法包括2个部分:第一部分为排序，第二部分为泛洪.算法可简述如下(其中步骤(ⅰ)为排序，步骤(ⅱ)—(ⅵ)为泛洪):

(ⅰ)首先计算图像中各点的梯度，然后扫描整幅图像得到各梯度的概率密度.各像素点在排序数组中的位置由梯度分布的累积概率与该像素点的梯度值计算得到.计算出所有像素点的排序位置并将其存入排序数组.在排序后的数组中，梯度值越低的点存放的位置越靠前.

(ⅱ)像素点按梯度值从低到高的顺序处理，相同梯度值的点作为一个梯度层级.

(ⅲ)处理一个梯度层级hcur(当前层)，将该层中所有邻域已被标识的点加入到一个先进先出的队列中去.

(ⅳ)若先进先出队列非空，则弹出队列的首元素作为当前处理像素.顺序处理当前像素所有高度为hcur的相邻点.若邻点已被标识，则根据该邻点标识刷新当前像素点的标识.若邻点尚未标识，则将该邻点加入到先进先出队列中去.循环执行本步骤直至队列空为止.

(ⅴ)再一次扫描当前梯度级的像素点，检查是否仍有未标识点.此时的未标识点意味着一个新的极小区.因此，若发现未标识点，则将当前区域标识值加1，并将该值赋为未标识点的标识值.然后从该点出发执行与步骤(ⅳ)相同的泛洪步骤，标识该极小区的所有像素点.

(ⅵ)返回步骤(ⅲ)处理下一梯度层级，直至将所有梯度层级都处理完为止.

3.2 标记分水岭分割

3.2.1 标记提取 为了改善传统分水岭过分割现象，采用标记的方法分别提取图像的内部标记和外部标记，辅助修改多尺度形态梯度图像，抑制无意义的局部极小值.

(1)内部标记提取:对多尺度形态梯度图像用Extend-H-maxima变换[15]计算局部最大值，将局部最大值图像转化为二值图像，从而提取内部标记.

(2)外部标记提取:对多尺度形态梯度图像用欧氏距离变换后提取外部标记[15].

3.2.2 图像分割 采用数学形态学中强制最小技术[16]修改多尺度形态梯度图像.具体步骤是:将内部标记和外部标记强制作为多尺度形态梯度图像的局部最小值，强制修改多尺度形态梯度图像，对修正后的梯度图像进行分水岭分割得到最终分割结果.

4 时频谱图分割步骤

时频谱图分割具体步骤为:

(ⅰ)利用形态重构滤波器去除时频谱图中的噪声.结合时频谱图噪声特点，这里选用半径为1的圆盘形结构元素进行开闭重构滤波，滤波中令r=f-1.

(ⅱ)利用(1)式求圆盘形结构元素半径分别为1，2，3时的单尺度形态梯度图像，并利用(2)式对3个单尺度形态梯度图像融合得到多尺度形态梯度图像.

(ⅲ)对步骤(ⅱ)中得到的多尺度形态梯度图像进行Extend-H-maxima变换[15](参数h= 10).

(ⅳ)对步骤(ⅱ)中得到的多尺度形态梯度图像进行欧氏距离变换，提取外部标记.

(ⅴ)对步骤(ⅲ)，(ⅳ)得到的内部标记及外部标记先进行闭操作，后膨胀操作.

(ⅵ)对步骤(ⅴ)中得到的标记图像利用Matlab中bwareaopen函数去除周边产生的孤立像素.

(ⅶ)采用形态学极小值标定技术修正多尺度形态梯度图像.

(ⅷ)对修正后的梯度图像进行分水岭算法分割.

5 实验结果与讨论

图2a是海豚音的Capon时频谱图[17]，水平方向表示时间，垂直方向表示频率.谱图中大体呈现为水平方向的间断的曲线是时频曲线，即为目标，其他部分均为背景.从图2a可以看出，干扰非常严重，存在大量的呈现出不完全纹理现象的强干扰(大体垂直方向的不规则条纹).干扰的强度与目标的强度相当，且与目标有重叠.对于这样的时频谱图，其分割是一个很大的挑战.

按时频谱图分割的步骤(ⅰ)得到去滤波后的时频谱图(图2b).按步骤(ⅱ)得到时频谱图的多尺度形态梯度图像(图2c).按步骤(ⅲ)—(ⅷ)得到新方法分割结果(图2d).将图2d中分割得到的目标置为白色，其他部分置为黑色，得到图2e.

作为比较，这里给出结构元素半径分别为1，2，3时单尺度形态梯度标记分水岭分割后的结果(图2f—h);另外，还给出最常用的Otsu法分割后的结果(图2i)，阈值为98.在上述的分割结果2e—i中仅仅是阈值化的结果，没有去除孤立区，但这并不影响考察分割的结果.总的来看，各种方法分割出的时频曲线均有一定的失真，没有分割出完整的时频曲线;这是因为时频谱图中干扰过于严重，干扰与时频曲线严重重叠.相比较而言，Otsu法分割出的时频曲线严重变宽，失真较大，这是因为时频谱图的直方图不是理想的双峰型，Otsu法不适用造成的.结构元素半径为2和3时单尺度形态梯度标记分水岭分割出的时频曲线残缺严重，丢失过半.在上述这些方法中，新方法和结构元素半径为1时单尺度形态梯度标记分水岭方法的分割效果相对较好.这两者相比，在分割出的时频曲线宽度失真方面，新方法失真较严重;在分割出的时频曲线细节保持方面，新方法较好;在方法的实用性方面，由于合适的单一尺度的选择相对比较困难，盲目性大，而多尺度的选择相对容易，因此新方法更有优势.综合来看，新方法是一个可取的折中.

a 海豚音的Capon时频谱图

b 滤波后时频谱图

c 多尺度形态学梯度图像

d 新方法分割

e 图2d结果重新赋值

f 结构元素半径为1时的分水岭分割

g 结构元素半径为2时的分水岭分割

h 结构元素半径为3时的分水岭分割

i Otsu法分割

6 结语

提出多尺度形态梯度图像的一种融合方法，并将该方法与标记分水岭分割结合形成一种基于多尺度形态梯度和标记分水岭的时频谱图分割方法.新方法思路是:利用图像熵定义对时频谱图3个不同尺度下的形态梯度图像进行融合，获得时频谱图的多尺度形态梯度图像，然后利用标记分水岭方法分割时频谱图的多尺度形态梯度图像.

新方法较最常用的Otsu法更适合于时频谱图分割，较基于单尺度形态梯度和标记分水岭的时频谱图分割方法更具实用性.由于时频谱图中干扰过于严重，干扰与时频曲线严重重叠，因此新方法分割出的时频曲线同样存在一定程度的失真.综合考虑方法的失真程度和实用性，新方法是一个较好的折中选择.

基于多尺度形态梯度和标记分水岭的时频谱图分割方法的分割效果还有待提高.认为以下2个问题需要进一步深入研究:(1)形成时频谱图的多尺度形态梯度图像时，如何根据时频谱图自身的特点选取不同尺度和不同形状的结构元素;(2)标记分水岭分割时，如何更有效地进行时频谱图多尺度形态梯度图像的标记.这2个问题的解决可望大幅度提高时频谱图的分割效果.

[1] 蔡艳平，李艾华，王 涛，等.基于时频谱图与图像分割的柴油机故障诊断[J].内燃机学报，2011，29(2):181-186.

[2] GAVROVSKA A M，PASKAS M P，DUJKOVIC D，et al.Region-Based Phonocardiogram Event Segmentation in Spectrogram Image[C].IEEE:The 10th Symposium on Neural Network Applications in Electrical Engineering，2010:69-72.

[3] GAVROVSKA A M，PASKAM P，RELJIN I S.Determination of Morphologically Characteristic PCG Segments from Spectrogram Image[J].Telfor Journal，2010，2(2):4-77.

[4] SHAFI I，AHMAD J，SHAH S I，et al.Time Frequency Image Analysis Using Neural Networks[C].IEEE:IMACS Multiconference on Computational Engineering in Systems Applications，2006:315-320.

[5] LEPRETTRE B，MARTIN N.Extraction of Pertinent Subsets from Time-Frequency Representations for Detection and Recognition Purposes[J].Signal Processing，2002，82(2):229-238.

[6] STEINBERG R，O’SHAUGHNESSY D.Segmentation of a Speech Spectrogram Using Mathematical Morphology[C].IEEE:IEEE International Conference on Acoustics，Speech and Signal Processing，2008:1 637 -1 640.

[7] YUWONO M，SU S W，MOULTON B D，et al.Gait Cycle Spectrogram Analysis Using a Torso-Attached Inertial Sensor[C].IEEE:Annual International Conference of the IEEE Engineering in Medicine and Biology Society，2012:6 539-6 542.

[9] ZHANG Zhenggang，XU Kailiang，TA Dean，et al.Joint Spectrogram Segmentation and Ridge-Extraction Method for Separating Multimodal Guided Waves in Long Bones[J].Science China-Physics，Mechanics and Astronomy，2013，56(7):1 317-1 323.

[10] MUKHOPADHYAY S，CHANDA B.An Edge Preserving Noise Smoothing Technique Using Multiscale Morphology[J].Signal Processing，2002，82:527-544.

[11] WANG D.A Multiscale Gradient Algorithm for Image Segmentation Using Watersheds[J].Pattern Recognition，1997，30(12):2 043-2 052.

[12] SALEMBIER P.Morphological Multiscale Segmentation for Image Coding[J].Signal Processing，1994，38(3):359-386.

[13] VINCENT L，SOILLE P.Watersheds in Digital Spaces:An Efficient Algorithm Based on Immersion Simulations[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1991，13(6):583-598.

[14] PUN T，KAPUR J N，WESZKA J S.Entropic Thresholding:A New Approach[J].Computer Vision，Graphics and Image Processing，1981，16(3):210-239.

[15] WU Qiang，MERCHANT F A，CASTLEMAN K R.Microscope Image Processing[M].Waltham，Massachusetts:Academic Press Elsevier，2008:135-142.

[16] SOILLE P.Morphological Image Analysis:Principle and Applications[M].New York:Springer-Verlag，2003.

[17] LEPRETTRE B，MARTIN N.Extraction of Pertinent Subsets from Time-Frequency Representations for Detection and Recognition Purposes[J].Signal Processing，2002，82(2):229-238.

(责任编辑 向阳洁)

Time-Frequency Spectrogram Segmentation Using the Multi-Scale Morphological Gradient and the Marked Watershed Algorithm

GUO Haitao1，2，XU Lei3，ZHAO Hongye3，TIAN Yuanyuan4，JIAO Shengxi3

(1.College of Electronic Information Engineering，Inner Mongolia University，Hohhot 010021，China;2.Electrical Engineering College，Northeast Dianli University，Jilin 132012，Jilin China;3.School of Automation Engineering，Northeast Dianli University，Jilin 132012，Jilin China;4.School of Mechanical Engineering，Northeast Dianli University，Jilin 132012，Jilin China)

It is important and difficult to segment a time-frequency spectrogram due to its strong interference as well as serious overlap among targets and between targets and interference.An image entropy based fusion method for multi-scale morphological gradient image of a time-frequency spectrogram is presented.By combining that method with the marked watershed algorithm，a method for segmenting time-frequency spectrogram based on the multi-scale morphological gradient image and the marked watershed algorithm，is obtained.The experiment results show that it is more practical than the segmentation method which is based on single scale morphological gradient and marker watershed，and more accurate than Otsu method.

time-frequency spectrogram;multi-scale morphological gradient image;image entropy;watershed algorithm;image segmentation

1007-2985(2015)04-0012-06

郭海涛(1965—)，男，黑龙江安达人，内蒙古大学电子信息工程学院教授，博士，主要从事图像处理、模式识别等研究.

TP391.4

A

10.3969/j.issn.1007-2985.2015.04.004