CN104268883A - 一种基于边缘检测的时频谱曲线提取方法 - Google Patents

Abstract

本发明公开了一种基于边缘检测的时频谱曲线提取方法，首先计算信号时频分布的k阶原点矩，然后通过边缘检测的方法获得含有时频谱边缘信息的矩阵，根据对该矩阵进行判断和对边缘信息进行合并处理，进而获得时频谱曲线，在根据已有的曲线模型建立参数空间对曲线的参数进行估计，完成对时频谱曲线的提取。与传统的检测方法相比，本发明充分利用了信号时频分布的结构属性，提高了曲线提取的准确性，尤其在低信噪比和有多条时频谱带的情况下，传统提取方法的性能有所下降甚至失效，但本方法仍可以提取出与时频谱曲线。对于时频谱曲线的提取，本发明的方法性能更优。

Description

一种基于边缘检测的时频谱曲线提取方法

技术领域

本发明属于雷达信号处理技术领域，涉及一种基于边缘检测的时频谱曲线提取方法。

背景技术

时频分析提供了信号时间域和频率域的联合分布信息，描述了信号在不同时间和频率的能量和密度。时频分析在通信、自动化、雷达、医学、地球物理和故障诊断等几乎所有工程技术领域中都得到了广泛应用。对时频谱曲线的提取是优化时频分布，点对点计算各个时刻的瞬时频率，从而进一步研究。

目前，提取时频谱曲线的方法主要采用脊线提取和质心曲线提取等方法。这两种方法是利用时频分布特征，在每个时刻对时频分布矩阵进行最值或质心计算获得该时刻的即时频率信息。在有噪声或有多条时频谱带的条件下进行提取时，时频分布在每一时刻的最值或加权平均值会与真值产生较大偏离，使得提取出的曲线不能准确反映即时频率信息。而边缘检测是标识数字图像中亮度变化明显的点，剔除不相关的信息，提取包含图像重要的结构属性的边界信息。这样，有利于标识出时频分布中的奇异点和突变点，滤除部分噪声对原信号时频分布的影响，也可以同时对多条时频谱曲线进行检测。

发明内容

本发明的目的在于提供一种基于边缘检测的时频谱曲线提取方法，解决了已有的时频谱曲线提取方法在噪声显著情况下的不足以及只能提取单条曲线的问题。

本发明所采用的技术方案是按照以下步骤进行：

步骤1：计算信号的时频分布X×Y维矩阵W_s(n,f_m)中各元素的k次方(即增大矩阵各元素之间的差值，有利于边缘检测剔除噪声影响)，得到一个新的X×Y维矩阵W(x,y)＝W_s(n,f_m)^k，其中n表示时间点序列，f_m表示频率点序列，x，y分别表示新矩阵W_X×Y的行和列；

步骤2，对步骤1中的矩阵W(x，y)使用Roberts算子进行边缘检测：

$G(x,y)=\sqrt{{|W(x,y)-W(x+1,y+1)|}^2-{|W(x+1,y)-W(x,y+1)|}^2}$

得到时频谱分布W_s(n,f_m)的梯度幅值G_X×Y，选取梯度幅值G_X×Y的平均值做为阈值，与G_X×Y进行比较，将G_X×Y化为二值边缘矩阵：

$H(x,y)=\left\{\begin{array}{cc}0& G(x,y)<T_g\\ 1& G(x,y)\&GreaterEqual;T_g\end{array}\right.$

得到只含有0和1的二值边缘矩阵H_X×Y；

步骤3，对步骤2中的边缘矩阵H_X×Y处理，计算边缘矩阵H_X×Y中每列的边缘点个数(X为矩阵H_X×Y的行数，每列为1的点为边缘点，为0的不是边缘点。)，逐个选出r(y)为奇数的列，再记录这些列中第一个边缘点的位置，将二值边缘矩阵H_X×Y这些位置元素置0，即删去矩阵中这个边缘点，将每列都检索完毕之后，可以保证矩阵H_X×Y中每列含有的边缘点均为偶数，再按照列优先顺序，即矩阵中各元素按照其在矩阵H_X×Y中行数和列数按照如下公式排序((列数-1)×X+行数)，(X为矩阵行的维度)，将矩阵H_X×Y中每个边缘点排序后的位置记录在向量d_N中,其中N为边缘点的总个数。最后将边缘矩阵H_X×Y边缘点两两合并，得到矩阵Z_X×Y；

步骤4，根据已有的曲线模型y＝f(x,T,θ,A,...)，如正弦曲线模型其中T为曲线周期，θ为曲线的相位，A为曲线幅度，建立I×J维参数空间(T,θ)，I，J分别为参数空间中周期时间点数和相位角度点数，将参数空间中每个点带入曲线模型中对Z_X×Y中的所有曲线点进行计算，即对每个参数空间点(T_i,θ_j)带入曲线模型中，得到再将Z_X×Y中所有边缘点位置的行和列带入参数点的模型中，若边缘点行和列(x,y)满足带入参数的模型，则累加器∑g(T,θ)累加在参数空间(T,θ)的每个点上加1，如此取遍所有参数空间的点，每次均对Z_X×Y中边缘点进行计算，最后得出各参数点的累加结果，找出参数累加器g中局部最大值点，局部最大值的个数为曲线的条数，各局部最大点在累加器中对应的周期值和相位值即为各条曲线的周期T和相位θ的具体数值，即得到曲线的条数以及各条时频谱曲线的参数。

进一步，所述步骤3中排序过程为：

将X×Y维0矩阵Z_X×Y中按照之前列优先排序，将位置为的元素置1，即标记为边缘点。

其中，N为矩阵H_X×Y中边缘点的总数，矩阵Z_X×Y其它位置元素为0，得到时频谱的曲线矩阵Z_X×Y。

本发明的有益效果是能在低信噪比以及多条时频谱的情况下完成对时频谱曲线的提取，获得即时频率。

附图说明

图1为本发明即基于边缘检测的时频谱曲线提取方法的流程图；

图2是针对本发明具体实施方式给出的仿真例子，给出了在信噪比为-5dB时，提取单条时频谱曲线本发明与传统提取方法的性能比较；

图3是针对提取多条时频谱曲线时，本发明与传统提取方法的性能比较。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

如图1所示为本发明基于边缘检测的时频谱曲线提取方法的流程。本发明首先计算信号时频分布各元素的k次方，然后通过边缘检测的方法获得含有时频谱边缘信息的矩阵，根据对该矩阵进行判断和对边缘信息进行合并处理，进而获得时频谱曲线，在根据已有的曲线模型建立参数空间对曲线的参数进行估计，完成对时频谱曲线的提取。与传统的检测方法相比，本发明充分利用了信号时频分布的结构属性，提高了曲线提取的准确性，尤其在低信噪比和有多条时频谱带的情况下，传统提取方法的性能有所下降甚至失效，但本方法仍可以提取出时频谱曲线。对于时频谱曲线的提取，本发明的方法性能更优。

本发明的目的是克服已有的时频谱曲线提取方法在噪声显著情况下的不足以及只能提取单条曲线的约束性，基本思路是对时频分布矩阵计算各元素的高阶次方，通过边缘检测提取出时频谱边缘矩阵，然后对边缘信息矩阵进行边缘合并处理，进而得到时频谱对应的曲线，即获得点对点信号的及时频率。再根据已有的曲线模型y＝f(x,T,θ,a,b...)设置检测参数空间，获得每条曲线的各个参数信息，完成对曲线的提取。

本发明一种基于边缘检测的时频谱曲线提取方法，包括以下步骤：

步骤1，计算信号的时频分布X×Y维矩阵W_s(n,f_m)中各元素的k次方(即增大矩阵各元素之间的差值，有利于边缘检测剔除噪声影响)，得到一个新的X×Y维矩阵W(x,y)＝W_s(n,f_m)^k，其中n表示时间点序列，f_m表示频率点序列，x，y分别表示新矩阵W_X×Y的行和列。

步骤2，对步骤1中的矩阵W(x，y)使用Roberts算子进行边缘检测：

$G(x,y)=\sqrt{{|W(x,y)-W(x+1,y+1)|}^2-{|W(x+1,y)-W(x,y+1)|}^2}$

得到时频谱分布W_s(n,f_m)的梯度幅值G_X×Y，选取梯度幅值G_X×Y的平均值做为阈值，与G_X×Y进行比较，将G_X×Y化为二值边缘矩阵：

$H(x,y)=\left\{\begin{array}{cc}0& G(x,y)<T_g\\ 1& G(x,y)\&GreaterEqual;T_g\end{array}\right.$

得到只含有0和1的二值边缘矩阵H_X×Y。

步骤3，对步骤2中的边缘矩阵H_X×Y处理，计算边缘矩阵H_X×Y中每列的边缘点个数(X为矩阵H_X×Y的行数，每列为1的点为边缘点，为0的不是边缘点)，逐个选出r(y)为奇数的列，再记录这些列中第一个边缘点的位置，将二值边缘矩阵H_X×Y这些位置元素置0，即删去矩阵中这个边缘点，将每列都检索完毕之后，可以保证矩阵H_X×Y中每列含有的边缘点均为偶数，再按照列优先顺序，即矩阵中各元素按照其在矩阵H_X×Y中行数和列数按照如下公式排序((列数-1)×X+行数)，(这里X为矩阵行的维度)，将矩阵H_X×Y中每个边缘点排序后的位置记录在向量d_N中,其中N为边缘点的总个数。最后将边缘矩阵H_X×Y边缘点两两合并，得到时频谱的曲线矩阵Z_X×Y。具体是将X×Y维0矩阵Z_X×Y中按照之前列优先排序，将位置为的元素置1，即标记为边缘点。

其中，N为矩阵H_X×Y中边缘点的总数，矩阵Z_X×Y其它位置元素为0，得到时频谱的曲线矩阵Z_X×Y。

步骤4，根据已有的曲线模型y＝f(x,T,θ,A,...)，如正弦曲线模型其中T为曲线周期，θ为曲线的相位，A为曲线幅度，建立I×J维参数空间(T,θ)，其中I，J分别为参数空间中周期时间点数和相位角度点数，将参数空间中每个点带入曲线模型中对Z_X×Y中的所有曲线点进行计算，即对每个参数空间点(T_i,θ_j)带入曲线模型中，得到再将Z_X×Y中所有边缘点位置的行和列带入参数点的模型中，若边缘点行和列(x,y)满足带入参数的模型，则累加器∑g(T,θ)累加在参数空间(T,θ)的每个点上加1，如此取遍所有参数空间的点，每次均对Z_X×Y中边缘点进行计算，最后得出各参数点的累加结果，找出参数累加器g中局部最大值点，局部最大值的个数为曲线的条数，各局部最大点在累加器中对应的周期值和相位值即为各条曲线的周期T和相位θ的具体数值，即得到曲线的条数以及各条时频谱曲线的周期和相位参数。

图2给出了对某一信噪比为-5dB信号的时频谱曲线提取时，本方法与脊线提取和质心曲线等传统提取方法的比较。在低信噪比的条件下，脊线提取得到的时频谱曲线在某些时间点上会产生跳跃，误差较大；质心曲线提取的曲线由于噪声影响与时频谱曲线也有很大差异；而本方法所提取的曲线与原时频分布更加近似。从仿真结果可以看出，显然本方法的性能要更好。

图3给出了某一信号的时频分布含有多条时频谱带时，本方法与脊线提取和质心曲线等传统提取方法的比较。由于传统方法的算法限制，只能提取出一条与时频谱不一致的曲线，而本方法可以同时提取出多条与时频谱近似的曲线，显然本方法提取曲线的性能更好。

本发明的有益效果：本发明通过对信号的时频分布进行边缘检测来提取时频谱曲线信息，先检测时频分布各元素的k次方矩阵的边缘再对边缘进行处理合并，最后对提取时频谱曲线进行参数估计达到曲线提取目的，解决了常见的质心曲线和脊线检测方法低信噪比检测效果差，无法同时提取多条曲线的缺点，该方法能在低信噪比以及多条时频谱的情况下完成对时频谱曲线的提取，获得即时频率。

下面结合附图和具体实施例对本发明做进一步的说明，并与传统检测方法做比较。

施例1：设某单目标振动由0.125Hz和0.6Hz两个分量组成，其中频率为0.125Hz分量的振动幅度为0.9mm，0.6Hz分量的振动幅度为5.5mm，雷达载波频率α＝500GHz，波长采样频率为500Hz。

对该目标的回波信号在0～30s内进行采样后的信号进行STFT，获得256×15000维的时频分布矩阵，并对该矩阵的各元素求2次方，得到矩阵W_256×15000。用Roberts算子对矩阵W_256×15000进行边缘提取，并二值化得到只含有0和1的边缘矩阵H_256×15000。

计算边缘矩阵H_256×15000中每列的边缘点个数逐个选出r(y)为奇数的列，再记录这些列中第一个边缘点的位置，将二值边缘矩阵H_256×15000这些位置元素置0，即删去矩阵中这个边缘点，再按照列优先顺序，即矩阵中各元素按照其在矩阵H_256×15000中行数和列数按照如下公式排序((列数-1)×256+行数)，将矩阵H_256×15000中每个边缘点的排序后的位置记录在向量d_N中，其中N为边缘点的总个数。最后将边缘矩阵H_256×15000边缘点两两合并，得到矩阵Z_256×15000。具体是将256×15000维0矩阵Z_256×15000中按照之前列优先排序，将位置为的元素置1，即标记为边缘点。

矩阵Z_256×15000其它位置元素为0，得到时频谱的曲线矩阵Z_256×15000。其图像如图2所示，该单目标时频分布曲线只有一条，所以无需知道其具体曲线模型，直接对矩阵Z_256×15000按列求加权和即可得到曲线与时间的关系，即完成了对曲线的提取。

进一步可以再进行FFT可以求出该曲线的频率分量，分别为0.125Hz和0.6Hz，也可以再对获得的曲线进行其他处理以获得其他信息，这里不进行叙述。

图2给出了对信噪比为-5dB该回波信号的时频谱曲线提取时，本方法与脊线提取和质心曲线等传统提取方法的比较。在低信噪比的条件下，脊线提取得到的时频谱曲线在某些时间点上会产生跳跃，误差较大；质心曲线提取的曲线由于噪声影响与时频谱曲线也有很大差异；而本方法所提取的曲线与原时频分布更加近似。从仿真结果可以看出，显然本方法的性能要更好。

施例2：设某目标在距离雷达R₀＝50m处以初速度v＝1.5m/s，加速度a＝0.14m/s²向雷达运动，目标上有两部分A,B均以f_p＝0.1Hz的频率做正弦振动，其中A部分振幅为r_A＝8.2m，B部分振幅为r_B＝5m，雷达载波频率α＝1.5GHz，波长采样频率为1000Hz。则目标的回波为：

$S\left(t\right)=e^{-2\mathrm{\&pi;j}\left({2\mathrm{\&lambda;}}^{-1}{R}_0\left(t\right)\right)}+e^{-2\mathrm{\&pi;j}\left(2{\mathrm{\&lambda;}}^{-1}{R}_A\left(t\right)\right)}+e^{-2\mathrm{\&pi;j}\left({2\mathrm{\&lambda;}}^{-1}{R}_B\left(t\right)\right)}$

其中j为虚数单位，t为时间，为目标主体距离雷达的瞬时距离，为A部分与雷达的瞬时位置关系，R_B(t)＝R₀(t)+r_Bsin(2πf_pt+π)为B部分与雷达的瞬时位置关系，各式中参数如前所述。

对该目标的回波信号在0～30s内进行采样后的信号进行STFT，获得256×30000维的时频分布矩阵，并对该矩阵的各元素求2次方，得到矩阵W_256×30000。用Roberts算子对矩阵W_256×30000进行边缘提取，并二值化得到只含有0和1的边缘矩阵H_256×30000。

计算边缘矩阵H_256×30000中每列的边缘点个数逐个选出r(y)为奇数的列，再记录这些列中第一个边缘点的位置，将二值边缘矩阵H_256×30000这些位置元素置0，即删去矩阵中这个边缘点，再按照列优先顺序，即矩阵中各元素按照其在矩阵H_256×30000中行数和列数按照如下公式排序((列数-1)×256+行数)，将矩阵H_256×30000中每个边缘点的排序后的位置记录在向量d_N中，其中N为边缘点的总个数。最后将边缘矩阵H_256×30000边缘点两两合并，得到矩阵Z_256×30000。具体是将256×30000维0矩阵Z_256×30000中按照之前列优先排序，将位置为的元素置1，即标记为边缘点。

矩阵Z_256×30000其它位置元素为0，得到时频谱的曲线矩阵Z_256×30000。其图像如图3所示，该图像就反映了时频谱线的特征。由于该目标时频分布曲线有多条，所以需要进一步对矩阵Z_256×30000进行直线和曲线检测得到具体时频谱线的特征，如本例中就是对矩阵Z_256×30000使用Hough变换，检测其中直线的参数，再使用正弦曲线检测，获得其中曲线部分的参数，完成对时频谱曲线的提取。

图3给出了该目标回波信号的时频分布，本方法与脊线提取和质心曲线等传统提取方法的比较。其时频分布含有多条时频谱带，由于传统方法的算法限制，只能提取出一条与时频谱不一致的曲线，而本方法可以同时提取出多条与时频谱近似的曲线，显然本方法提取曲线的性能更好。

以上详细描述了本发明的较佳具体实施例，本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (2)

1.一种基于边缘检测的时频谱曲线提取方法,其特征在于按照以下步骤进行：

步骤1：计算信号的时频分布X×Y维矩阵W_s(n,f_m)中各元素的k次方(即增大矩阵各元素之间的差值，有利于边缘检测剔除噪声影响)，得到一个新的X×Y维矩阵W(x,y)＝W_s(n,f_m)^k，其中n表示时间点序列，f_m表示频率点序列，x，y分别表示新矩阵W_X×Y的行和列；

步骤2：对步骤1中的矩阵W(x，y)使用Roberts算子进行边缘检测：

$G(x,y)=\sqrt{{|W(x,y)-W(x+1y+1)|}^2-{|W(x+1,y)-W(x,y+1)|}^2}$

得到时频谱分布W_s(n,f_m)的梯度幅值G_X×Y，选取梯度幅值G_X×Y的平均值做为阈值，与G_X×Y进行比较，将G_X×Y化为二值边缘矩阵：

$H(x,y)=\left\{\begin{array}{cc}0& G(x,y)<T_g\\ 1& G(x,y)\&GreaterEqual;T_g\end{array}\right.$

得到只含有0和1的二值边缘矩阵H_X×Y；

步骤3：对步骤2中的边缘矩阵H_X×Y处理，计算边缘矩阵H_X×Y中每列的边缘点个数(X为矩阵H_X×Y的行数，每列为1的点为边缘点，为0的不是边缘点)，逐个选出r(y)为奇数的列，再记录这些列中第一个边缘点的位置，将二值边缘矩阵H_X×Y这些位置元素置0，即删去矩阵中这个边缘点，将每列都检索完毕之后，可以保证矩阵H_X×Y中每列含有的边缘点均为偶数，再按照列优先顺序，即矩阵中各元素按照其在矩阵H_X×Y中行数和列数按照如下公式排序((列数-1)×X+行数)，(这里X为矩阵行的维度)，将矩阵H_X×Y中每个边缘点排序后的位置记录在向量d_N中,其中N为边缘点的总个数，最后将边缘矩阵H_X×Y边缘点两两合并，得到矩阵Z_X×Y；

步骤4，根据已有的曲线模型y＝f(x,T,θ,A,…)，如正弦曲线模型其中T为曲线周期，θ为曲线的相位，A为曲线幅度，建立I×J维参数空间(T,θ)，I，J分别为参数空间中周期时间点数和相位角度点数，将参数空间中每个点带入曲线模型中对Z_X×Y中的所有曲线点进行计算，即对每个参数空间点(T_i,θ_j)带入曲线模型中，得到再将Z_X×Y中所有边缘点位置的行和列带入参数点的模型中，若边缘点行和列(x,y)满足带入参数的模型，则累加器Σg(T,θ)累加在参数空间(T,θ)的每个点上加1，如此取遍所有参数空间的点，每次均对Z_X×Y中边缘点进行计算，最后得出各参数点的累加结果，找出参数累加器_g中局部最大值点，得到曲线的条数以及各条时频谱曲线的参数。

2.按照权利要求1所述一种基于边缘检测的时频谱曲线提取方法，其特征在于：所述步骤3中排序过程为将X×Y维0矩阵Z_X×Y中按照之前列优先排序，将位置为的元素置1，即标记为边缘点，

其中，N为矩阵H_X×Y中边缘点的总数，矩阵Z_X×Y其它位置元素为0，得到时频谱的曲线矩阵Z_X×Y。