基于图像分割技术产生滤波器的自适应分解方法

摘要

本发明提出一种基于图像分割技术产生滤波器的自适应分解方法，包括以下步骤：第一步，对待分析信号进行时频变换，得到对应的时频系数和时频平面；第二步，对时频平面进行阈值分割，以得到二值化图；第三步，标记二值化图中的连通域；第四步，根据连通域生成一组时频滤波器，并根据时频滤波器对时频系统进行滤波，输出滤波结果；第五步，对滤波结果进行时频逆变换，以得到分解结果。本发明的方法在待分解信号各频率成分采样时间足够长时，可以具备任意频率分辨率，与常用的经验模态分解算法相比，大大提高了频率分辨率，可以对频率密集分布的信号进行自适应分解，并以此计算各个成分的瞬时频率和瞬时幅值，从而扩宽自适应分解算法的应用范围。

说明书

技术领域

本发明涉及信号处理技术领域，特别涉及一种基于图像分割技术产生滤波器的自适应分解方法。

背景技术

目前，自适应分解算法在信号处理领域得到非常广泛的应用。其典型代表经验模态分解(Empirical Mode Decomposition，EMD)的发布论文迄今已经在谷歌学术中有了超过12000次的引用，由此可见一斑。究其原因，主要是因为EMD方法与Hilbert变换结合，提供了一种便利的计算信号各频率成分瞬时频率和瞬时幅值的途径。作为定量评价信号特征常用的参数，瞬时频率和瞬时幅值被广泛应用。瞬时频率在诸如雷达、声纳、机械系统故障诊断的工程领域广泛应用。同样瞬时幅值应用于模态参数确定和无损检测检测领域。

当前对瞬时频率和瞬时幅值进行计算的方法主要可分为两类：一类是基于内积计算的时频变换方法，例如，短时傅里叶变换、小波变换、框架分解等；另一类是结合自适应分解算法(例如，EMD，局部均值分解(Local Mean Decompositio，LMD)和经验小波变换分解(Empirical Wavelet transform,EWT))和瞬时频率和瞬时幅值的瞬时算法(例如，Hilbert变换和能量算子法)。

由于分解波形的多样性，前者的不同分量更容易聚积在时频平面中的对应区域，易于分开。但由于内积计算是在某段待分析信号与分解波形之间进行的，所以这不是一种真正的瞬时计算方法，计算结果不能精确反映信号频率和幅值的瞬时变化。对于后者，其优点在于计算过程包括了两个阶段：第一阶段将信号分解为不同的单独分量，例如，在EMD方法中的本征模函数(Intrinsic Mode Function，IMF)；第二阶段采用瞬态方法计算每个分量的瞬时频率和瞬时幅值，例如Hilbert变换和能量算子法。所以采用此方法能更准确计算信号的瞬时频率和瞬时幅值。但此方法的频率分辨率不够理想。迄今，对于LMD和EWT方法尚未有关于频率分辨的报道。对于经典的自适应分解算法EMD方法，已经有文献报道，但当两个分量的频率比值大于0.75时(f_低/f_高)，该方法便不能将两者分开，从而限制自适应分解算法的应用范围。

发明内容

本发明旨在至少解决上述技术问题之一。

为此，本发明的目的在于提出一种基于图像分割技术产生滤波器的自适应分解方法，该方法在待分解信号各频率成分采样时间足够长时，可以具备任意频率分辨率，与常用的经验模态分解算法相比，大大提高了频率分辨率，可以对频率密集分布的信号进行自适应分解，并以此计算各个成分的瞬时频率和瞬时幅值，从而扩宽自适应分解算法的应用范围。

为了实现上述目的，本发明的实施例提出了一种基于图像分割技术产生滤波器的自适应分解方法，包括以下步骤：S1：对待分析信号进行时频变换，得到对应的时频系数和时频平面；S2：对所述时频平面进行阈值分割，以得到二值化图；S3：标记所述二值化图中的连通域；S4：根据所述连通域生成一组时频滤波器，并根据所述时频滤波器对所述时频系统进行滤波，输出滤波结果；以及S5：对所述滤波结果进行时频逆变换，以得到分解结果。

另外，根据本发明上述实施例的基于图像分割技术产生滤波器的自适应分解方法还可以具有如下附加的技术特征：

在一些示例中，在所述S1中，所述时频变换的方式包括：短时傅里叶变换或小波变换。

在一些示例中，所述短时傅里叶变换的变换公式为：

$>{\mathrm{STFT}}_w\left(s\right)[m,n]=\underset{k=0}{\overset{K-1}{\Sigma }}s\left[k\right]w[k-l]e^{-j2\pi n(k-m)/K},$>

其中，STFT_w(s)[m,n]为点(m,n)处的短时傅里叶变换系数，m为时频平面的时间变量，n为频率变量，s[k]为待分析的离散信号，k为抽样点的时间变量，K为采样时间长度，w[k-l]为窗函数，且||w||＝1。

在一些示例中，在所述S2中，所述得到二值化图的变换公式为：

$>\hat{c}(m,n)=\left(\begin{array}{cc}1& c(m,n)>T_0\\ 0& c(m,n)\le T_0\end{array}\right),$>

其中，为所述二值化图，c(m,n)为所述时频系数，T₀为分割阈值。

在一些示例中，所述S3进一步包括：S31：逐行扫描所述二值化图，将所述二值化图的每一行中连续的值为1的像素点组成的一个序列作为一个团，并记下所述团的起点和终点以及所述团所在的行号；S32：对除第一行外的所有行中的团进行扫描，如果当前团与前一行中的所有团都没有重合区域，则给所述当前团一个新的标号，如果当前团仅与上一行中一个团有重合区域，则用所述当前团的标号进行标记，如果当前团与上一行的2个以上的团有重叠区域，则给所述当前团赋予联通团中的最小标号，并将上一行中的几个团的标记写入等价对；S33：将所述等价对转换为等价序列，并对每一个等价序列赋予一个相同标号；S34：遍历开始团的标记，查找等价序列，并给予所述团新的标记；S35：将每个团的标号填入所述二值化图像中。

在一些示例中，在所述S4中，所述滤波结果为：

$>{\hat{C}}_i-F_i\left(C\right),$>

其中，(i＝1,2,3,…,P)为所述滤波结果，P为时频滤波器个数，F_i为第i个时频滤波器，C为所述时频系数的矩阵。

在一些示例中，在所述S5中，所述时频逆变换的变换公式为：

$>s\left[k\right]=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\mathrm{STFT}}_w\left(s\right)[m,n]w[k-M]e^{j2\pi n(k-m)/K},$>

其中，STFT_w(s)[m,n]为点(m,n)处的短时傅里叶变换系数，m为时频平面的时间变量，n为频率变量，s[k]为待分析的离散信号，k为抽样点的时间变量，K为采样时间长度，w[k-M]为窗函数，且||w||＝1。

根据本发明实施例的基于图像分割技术产生滤波器的自适应分解方法，结合时频变换和逆变换方法，将图像分割方法引入，用于构造时频滤波器，来实现信号成分的分解，且在待分解信号各频率成分采样时间足够长时，可以具备任意频率分辨率，与常用的经验模态分解算法相比，大大提高了频率分辨率，可以对频率密集分布的信号进行自适应分解，并以此计算各个成分的瞬时频率和瞬时幅值，从而扩宽自适应分解算法的应用范围。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是根据本发明实施例的基于图像分割技术产生滤波器的自适应分解方法的流程图；

图2是根据本发明一个具体实施例的示例信号的傅里叶变换时频图；

图3是根据本发明一个具体实施例的阈值分割后得到的二值化图的图像；以及

图4是根据本发明一个具体实施例的基于图像分割技术产生滤波器的自适应分解方法分解后基于Hilbert变换得到的瞬时频率和瞬时幅值的示意图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

以下结合附图描述根据本发明实施例的基于图像分割技术产生滤波器的自适应分解方法。

图1是根据本发明一个实施例的基于图像分割技术产生滤波器的自适应分解方法的流程图。如图1所示，该方法包括以下步骤：

步骤S1：对待分析信号进行时频变换，得到对应的时频系数和时频平面。

其中，时频变换的方式例如为短时傅里叶变换或小波变换。本实施例中所用的时频变换方法为短时傅里叶变换，其变换公式为：

$>{\mathrm{STFT}}_w\left(s\right)[m,n]=\underset{k=0}{\overset{K-1}{\Sigma }}s\left[k\right]w[k-l]e^{-j2\pi n(k-m)/K},$>

其中，STFT_w(s)[m,n]为点(m,n)处的短时傅里叶变换系数，m为时频平面的时间变量，n为频率变量，s[k]为待分析的离散信号，k为抽样点的时间变量，K为采样时间长度，w[k-l]为窗函数，且||w||＝1。

当然，本发明的时频变换方法不仅限于此，可以为任意变换结果为二维时频平面的时频变换算法，此处仅是以示例性为目的的描述。

步骤S2：对时频平面进行阈值分割，以得到二值化图。

在本发明的一个实施例中，二值化图的表达式例如为：

$>\hat{c}(m,n)=\left(\begin{array}{cc}1& c(m,n)>T_0\\ 0& c(m,n)\le T_0\end{array}\right),$>

其中，为二值化图，c(m,n)为时频系数，T₀为分割阈值。

步骤S3：标记二值化图中的连通域。

具体地，步骤S3进一步包括：

S31：逐行扫描二值化图，将二值化图的每一行中连续的值为1的像素点组成的一个序列作为一个团，并记下团的起点和终点以及团所在的行号。

S32：对除第一行外的所有行中的团进行扫描，如果当前团与前一行中的所有团都没有重合区域，则给当前团一个新的标号，如果当前团仅与上一行中一个团有重合区域，则用当前团的标号进行标记，如果当前团与上一行的2个以上的团有重叠区域，则给当前团赋予联通团中的最小标号，并将上一行中的几个团的标记写入等价对。

S33：将等价对转换为等价序列，并对每一个等价序列赋予一个相同标号。

S34：遍历开始团的标记，查找等价序列，并给予团新的标记。

S35：将每个团的标号填入二值化图像中。

步骤S4：根据连通域生成一组时频滤波器，并根据时频滤波器对时频系统进行滤波，输出滤波结果。

在本发明的一个实施例中，滤波结果例如为：

$>{\hat{C}}_i=F_i\left(C\right),$>

其中，(i＝1,2,3,…,P)为滤波结果，也即滤波后时频系数矩阵，P为时频滤波器个数，F_i为第i个时频滤波器，C为时频系数的矩阵。

步骤S5：对滤波结果进行时频逆变换，以得到分解结果。

在本发明的一个实施例中，对滤波结果进行时频逆变换的变换公式为：

$>s\left[k\right]=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\mathrm{STFT}}_w\left(s\right)[m,n]w[k-M]e^{j2\pi n(k-m)/K},$>

其中，STFT_w(s)[m,n]为点(m,n)处的短时傅里叶变换系数，m为时频平面的时间变量，n为频率变量，s[k]为待分析的离散信号，k为抽样点的时间变量，K为采样时间长度，w[k-M]为窗函数，且||w||＝1。

为了便于更好地理解本发明实施例的基于图像分割技术产生滤波器的自适应分解方法，以下结合附图和具体的实施例对本发明的方法进行进一步具体描述。

在本实施例中，该方法在matlab软件平台下进行实现。以一个具体的示例信号为例对该方法进行说明。该示例信号的表达式为：

s(t)＝s₁(t)+s₂(t)+s₃(t),0≤t≤1

s₁＝(t-t²)sin(1000·2π·t)

s₂＝(t-t²)sin(750·2π·t)

s₃＝(t-t²)sin(400·2π·t+200sin(2πt))，

即该示例信号中包括3个分量，其采样频率f_s为100kHz。

基于此，在本实施例中，该方法包括以下步骤：

步骤1：对示例信号进行时频变换。本示例中采用最为普通的短时傅里叶变换对示例信号进行时频变换，具体的变换公式为：

$>{\mathrm{STFT}}_w\left(s\right)[m,n]=\underset{k=0}{\overset{K-1}{\Sigma }}s\left[k\right]w[k-l]e^{-j2\pi n(k-m)/K},$>

其中，k代表时间，m(m＝0,1,2,…,M-1)代表时频变换中的时间移动，n(n＝0,1,2,…,N-1)代表时频变换中的频率变化。进一步地，得到的变换结果(即时频系数)c(m,n)例如图2所示。由图2可知，3个不同频率成分对应的时频系数在时频平面中分别聚积在不同区域。

步骤2：阈值分割步骤1中得到的时频平面产生二值化图。具体计算公式为：

$>\hat{c}(m,n)=\left(\begin{array}{cc}1& c(m,n)>T_0\\ 0& c(m,n)\le T_0\end{array}\right),$>

其中，为计算后的二值化图像的像素数值，T₀为分割阈值。此处，分割阈值设定为0.01，则变换结果如图3所示，由图3可知，阈值分割后，形成了3个独立的连通域，在其中包含了3个频率成分对应的时频系数。

步骤3：标记步骤2中得到二值化图中的连通域，3个连通区域被分别标记，便于程序自动完成分解过程。具体的标记连通域算法过程如下：

(1)逐行扫描步骤2中得到的二值化图，把每一行中连续的值为1的像素点组成的一个序列称为一个团，并记下它的起点和终点以及其所在的行号。

(2)对于除了第一行外的所有行里的团进行扫描，如果它与前一行中的所有团都没有重合区域，则给它一个新的标号；如果仅与上一行中一个团有重合区域，则将用该团的标号进行标记；如果其与上一行的2个以上的团有重叠区域，则给当前团赋予联通团中的最小标号，并将上一行的这几个团的标记写入等价对。

(3)将等价对转换为等价序列，每一个序列赋予一个相同标号。赋予每个等价序列一个标号，从1开始。

(4)遍历开始团的标记，查找等价序列，给予它们新的标记。

(5)将每个团的标号填入标记图像中。

步骤4：根据连通域产生一组时频滤波器，并对步骤1中得到的时频系数进行滤波，具体滤波公式为：

$>{\hat{C}}_i=F_i\left(C\right),$>

其中，(i＝1,2,3,…,P,P为时频滤波器个数)为滤波结果，也即滤波后时频系数矩阵，F_i为第i个时频滤波器，C为步骤1中得到的时频系数矩阵。

在本实施例中，产生了3组滤波器，滤波后，相应的也就得到了3组不同的时频系数。

步骤5：对步骤4中得到的滤波结果进行与步骤1中对应的时频逆变换，以得到分解结果，具体的计算公式为：

$>s\left[k\right]=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\mathrm{STFT}}_w\left(s\right)[m,n]w[k-M]e^{j2\pi n(k-m)/K},$>

其中，STFT_w(s)[m,n]为点(m,n)处的短时傅里叶变换系数，m为时频平面的时间变量，n为频率变量，s[k]为待分析的离散信号，k为抽样点的时间变量，K为采样时间长度，w[k-M]为窗函数，且||w||＝1。

进一步地，在得到分解结果之后，基于Hilbert变换计算分解结果的瞬时频率和瞬时频率，其结果如图4所示。由图4(a)可知，对于瞬时频率，除了端点位置，计算值与设定值吻合；由图4(b)可知，对于瞬时幅值，计算值与设定值曲线基本重合在一起。瞬时频率和瞬时幅值的计算结果表明了该算法有效地将3个成分分离。

综上，根据本发明实施例的基于图像分割技术产生滤波器的自适应分解方法，结合时频变换和逆变换方法，将图像分割方法引入，用于构造时频滤波器，来实现信号成分的分解，且在待分解信号各频率成分采样时间足够长时，可以具备任意频率分辨率，与常用的经验模态分解算法相比，大大提高了频率分辨率，可以对频率密集分布的信号进行自适应分解，并以此计算各个成分的瞬时频率和瞬时幅值，从而扩宽自适应分解算法的应用范围。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同限定。