Time-frequency Analysis of BGabor-NSPWVD Algorithm With Strong Robustness and High Sharpening Concentration
-
摘要: 针对短时傅里叶变换（Short-time Fourier transform，STFT）、Gabor变换和魏格纳-维尔分布（Wigner-Ville distribution，WVD）出现的时频分辨率模糊和交叉项干扰，以及目前一些主流改进算法如STFT-WVD和Gabor-WVD存在的频率分量三维幅度失真，且抗噪性能及鲁棒性能不理想等问题，提出基于局部二值化、归一化处理再结合的二值化Gabor-归一化WVD（Binarized Gabor-normalized WVD，BGabor-NWVD）和二值化Gabor-归一化伪平滑WVD（Binarized Gabor-normalized smoothed pseudo WVD，BGabor-NSPWVD）算法.数值仿真实验结果表明，BGabor-NWVD和BGabor-NSPWVD算法较好地抑制了交叉项干扰，具有较高的时频锐化聚集度，且两种算法的抗噪性能和鲁棒性也较为理想.基于本文方法对硬质合金顶锤工作时产生的疑似破裂信号进行时频分析，在抑制噪声和交叉项的同时能够较为准确地寻找传感器的频率判别窗口，为金属破裂监测设备数据采集卡提供有效的阈值参考.
-
关键词:
- 时频分析 /
- 交叉项抑制 /
- 鲁棒性 /
- BGabor-NWVD /
- BGabor-NSPWVD
Abstract: For the time-frequency resolution fuzziness and cross-term interference caused by short-time Fourier transform (STFT), Gabor transform and Wigner-Ville distribution (WVD), this paper proposes the binarized Gabor-normalized WVD (BGabor-NWVD) and binarized Gabor-normalized smoothed pseudo WVD (BGabor-NSPWVD) algorithms, which combine partial binaryzation and normalization. Experiment results show that the BGabor-NWVD and BGabor-NSPWVD algorithms can effectively reduce the cross-term interference and have high sharpening resolution of high TF concentration degree. The algorithms can suppress the cross-term, and the anti-noise performance and robustness of the two algorithms are also ideal. Based on the methods, a TF analysis of the suspected rupture signal generated during the work of the carbide torch is carried out to suppress the noise and the cross-term. At the same time, the frequency discrimination window of the sensor can be found accurately, which can be used as threshold reference for the data acquisition board of the metal rupture monitoring equipment.-
Key words:
- Time-frequency analysis /
- cross term inhibition /
- robustness /
- BGabor-NWVD /
- BGabor-NSPWVD
-
时频分析(Time-frequency analysis, TFA)是从时间-频率的角度来处理非平稳信号的方法, 其基本思想是建立时间和频率的联合函数, 描述信号在时间-频率平面的能量密度或强度[1-2].它将信号从时间域变换到时间-频率联合域, 不同频率分量的时间关联特性能够在时-频平面上有效地表示出来.信号在任意时刻的能量都聚集在此瞬时频率附近, 基于线性方法的逆变换则可以重构其等效的时间域信号[3].目前时频分析方法已经广泛应用于自动化控制、信号处理、数据通信、地质勘探、基础物理、工业生产等各个领域[4-5].在合成人造金刚石加工环节, 自动检测顶压机顶锤是否破裂是该项生产过程的重要步骤之一, 可有效降低人工凭经验听音辨别的误差.该金属破裂信号属于典型的非平稳信号, 利用合适的时频分析方法能够对此类信号进行有效的时频表示, 针对信号的破裂局部信息给出对应的高幅值频率分布, 为数据采集卡的频率判别窗口提供阈值依据.本文提出基于二值化Gabor的归一化Wigner-Ville分布(Wigner-Ville distrivution, WVD)和归一化伪平滑WVD (BGabor-NWVD, BGabor-NSPWVD)时频分析算法, 具有良好的高锐化时频聚集度和鲁棒性, 对金属破裂样本信号有较好的去噪效果, 可以得出有参考意义的时间和频率联合分布.
1. 时频分析
时频分析通过设计时间和频率的幅度或能量密度关联函数, 将一维的时间序列信号以二维的时间-频率密度函数形式表示, 旨在揭示信号中包含多少频率分量, 以及每一分量随时间的变化规律.时频分析方法可以分为线性型时频表示、非线性型时频表示和二次型时频表示. 1)典型的线性时频表示有短时傅里叶变换(Short-time Fourier transform, STFT)[6-7]、S变换(S transform, ST)[8]、和小波变换(Wavelet transform, WT)[9]等, 线性时频表示会受到不确定性原理的制约, 存在时频模糊等问题. 2)非线性型时频表示包括基于经验模型分解(Empirical mode decomposition, EMD)的希尔伯特-黄变换(Hilbert-Huang transformation, HHT)[10]、集合经验模型分解(Ensemble empirical mode decomposition, EEMD)及其改进算法[11-12]、同步压缩变换(Synchrosqueezing transform, SST)[13]等.基于EMD和EEMD的HHT具有良好的窄带自适应性, 但缺乏严谨的数学支撑, 且信号频率带宽较大时, 其时频聚集度有待提高. SST方法是Daubechies于2011年提出的一种时频重排技术, 其变换核可以基于小波变换或者基于STFT变换, SST支持信号重构, 其时频输出具有高锐化聚集度的特点[13], 但直接SST存在复杂频率成分交叉点附近的描述模糊和鲁棒性弱的缺点.改善SST的鲁棒性, 可采用与压缩感知稀疏方法[14]相结合的方法, 此部分本文作者另有相关文章阐述. 3)二次型时频分布则是一种更为严格的时频表示, 能够描述信号的瞬时功率谱密度, 可从时间-频率-能量三者联合分布的维度来了解信号的特点.二次型时频分布一般包括Cohen类时频分布[6]、WVD[15], 以及对WVD方法的改进, 如STFT-WVD[16]、STFT-SPWVD[17]和NSTFT-WVD[18]等.
常用的STFT方法由Gabor于1946年提出[7], 即加窗的傅里叶变换.窗函数可以选择矩形窗、三角窗、汉宁(Hanning)窗、海明(Hamming)窗、布拉克曼(Blackman)窗、恺撒(Kaiser)窗等, 当STFT的窗函数选择高斯窗时, 加窗短时傅里叶变换称为Gabor变换.根据Heisenberg测不准原理, 高斯函数窗口面积已达到测不准原理下界, 是时域窗口面积达到最小的函数, Gabor变换是最优的STFT.文献[19]采用稀疏分析窗的方法来求解离散Gabor变换, 一定程度提高了时频的聚集度. STFT和Gabor变换使用的是大小和形状固定的滑动窗口, 变换基函数为非正交系, 对于突变和非平稳信号, 不能精确分解周期比时间窗大的低频信号, 且高频的时频分辨率比较差, 需要辅以其他高时频聚集度方法加以改进, WVD就是较好的选择, 该方法具有理论上最高时频分辨率和许多优良的数学性质.
WVD是一种基础并且重要的二次型时频分布, 具有较高的时频聚集度, 应用广泛. Wigner分布在1932年被提出, 并应用于量子力学, 但没有引起重视, 直到Ville在1948年将它应用于信号分析领域, 并称之为Wigner-Ville分布, 开始了它的快速发展时期.尤其自上世纪八十年代以来, 陆续发表了相当数量的WVD论文, 对WVD的定义、性质、时频表示、各种改进进行了大量的讨论, 以其良好的时变特性, 特别是可直接由其精确定义"瞬时频率"而广泛用于非平稳信号的分析.同时需要注意的是, 这种方法虽然对边缘特性、瞬时频率和局域刻画等都有很好的描述, 其时频图上的频率成分有高锐化的聚焦, 但WVD不满足叠加原理, 其变换过程中产生新的频率交叉项分量成为该方法的瓶颈[18].交叉项的存在使得时频图上的分布混乱, 额外产生大量的噪声成分, 充斥在真实频率之间, 干扰对信号的识别分析, 极大地限制了WVD分布的实际应用.
综上所述, 采用Gabor与WVD结合对非平稳信号进行时频分析, 可发挥各自的极限优点, 去除WVD交叉项的同时, 得到较高的局部时频分辨率, 满足对非平稳信号进行时频分析的高锐化时频聚集度和强鲁棒性要求.本文对Gabor和WVD分别进行改进后再结合, 提出BGabor-NWVD和BGabor-NSPWVD算法, 同时兼顾线性时频表示与二次型时频表示的优点, 能够获得较好的时频聚集度、消除WVD交叉项干扰、抑制噪声, 并具有强鲁棒性, 同时能在金属破裂样本信号的时频分析中有较好的应用.
2. BGabor-NWVD和BGabor-NSPWVD
2.1 Gabor与WVD的优缺点
STFT与Gabor变换是常用的时频分析方法之一.其优点是, 对信号具有一定程度的自适应性, 能够如实还原输入信号的时频分布, 对信号的低频和高频两端都有较为一致的局部描述, 得到的非平稳信号的时频分布较为平滑, 分析多分量信号的时候, 没有交叉项的虚假频率成分干扰. STFT变换的定义如式(1), 把STFT的窗函数变为高斯函数, 就得到式(2)的Gabor变换[7].
$ \begin{align} &{\rm STFT}_x(t, f)=\int x(\tau)g^{*}(\tau-t){\rm e}^{-j2\pi\tau f}{\rm d}\tau \end{align} $
(1) $ \begin{align} &G_x(t, f)=\int x(\tau){\rm e}^{-\pi (\tau -t)^2}{\rm e}^{-j2\pi \tau f}{\rm d}\tau \end{align} $
(2) Gabor的时频分辨率由窗函数的时域大小直接决定, 一旦窗口函数选定, 其时频分辨率就已确定, 不随时间和频率的变化而变化.为了提高时间分辨率, 窗函数的时间宽度应该尽量短, 但为了提高频率分辨率, 窗函数的时间宽度则应尽可能长, 受不确定性原理的约束, 时间分辨率和频率分辨率是一对矛盾体.在处理非平稳信号的过程中, 对于高频信息, 需要用较窄的窗函数进行分析, 而对于低频信息, 则需要用较宽的窗函数进行分析.而Gabor的窗函数确定以后, 只能以一种固定分辨率进行时频分析, 无法兼顾高频信息和低频信息[20-21].这种方法的缺点很明显, 其时频聚集性不佳, 局部的频带粗糙, 难以精确显示时频特点, 如图 1所示, 需要通过其他方法提高其时频聚集度.
图 1 信号$f_1$的理想时频、Gabor和WVD对比图Fig. 1 Ideal time-frequency spectrum, Gabor, WVD of $f_1$WVD是最基本且严格的时间和频率二维联合函数, 可以看作信号在时间和频率平面上的能量密度解.虽然时频表示的线性特性是我们所希望具备的重要特性, 但因为能量本身就是一种二次型表示, 因此, 当反映信号的时间-能量分布时, 采用WVD这样的二次型时频分布则更加合适. WVD属于Cohen类双线性时频分布的一种, 可在时域和频域同时揭示信号的能量分布, 并且其物理意义明确.近年来, 凭借着其优越的时频聚集性, WVD被广泛应用于信号分析和处理领域, 尤其在非平稳信号的瞬时频率估计、信号的相干检测和时变滤波等诸多领域, 是最常用的双线性时频分布之一. WVD分布定义为:
$ \begin{align} &W_z(t, f)=\int z(t+\frac{\tau}{2})z^{*}(t-\frac{\tau}{2}){\rm e}^{-j2\pi \tau f}{\rm d}\tau \end{align} $
(3) 其中, $z(t)$为目标信号的解析信号, $\tau$是积分变量, $t$是时移, $f$是频率.式(3)没有使用窗函数, 避免了线性时频表示时间和频率分辨率相互制约的矛盾. WVD的时间带宽积达到了不确定原理给出的下界, 可以呈现较其他时频分布更好的时频聚集特性.但WVD在处理复杂信号, 尤其是频率相近的加性信号的时候, 由自身变换过程所引入的交叉项问题同样很严重.设$z(t)=z_1(t)+z_2(t)$, 则$z(t)$的WVD如式(4).
$ \begin{align} W_z(t, f)=\, &W_{z1}(t, f)+W_{z2}(t, f)+\nonumber\\&2{\rm Re}\{W_{z1, z2}(t, f)\} \end{align} $
(4) 其中, 等号右边的$W_{z1}(t, f)$与$W_{z2}(t, f)$分别为$z_1(t)$和$z_2(t)$的WVD, 即变换所要求解的时频信息项, 第3项$2{\rm Re}\{ W_{z1, z2}(t, f)\}$为交叉项.由式(4)可知, 两个信号和的WVD并不等于它们各自WVD之和, 交叉项的存在给WVD的应用带来不小的困扰.
$ \begin{align} f_1(t)=\, &\cos(t^2+t+\cos(t))+\nonumber\\&\cos(4\pi t-\sin(0.8\pi t)) \end{align} $
(5) 式(5)是信号$f_1(t)$的时域表达式, 图 1为信号$f_1(t)$理想时频表示图、Gabor和WVD的时频分布图.从图中可以看出, Gabor变换的时频分辨率不高, 而WVD方法的分辨率尽管略高于Gabor变换, 但存在严重的交叉项干扰问题.因此, 单独的Gabor和WVD都不理想, 需要进行一定的改进, 以克服时频分辨率和交叉项的影响.
2.2 Gabor与WVD的结合思想
Gabor-WVD变换的基本思想是利用Gabor和WVD各自的优点, 通过两者重叠运算来增强Gabor和WVD信息项, 抑制由WVD变换产生的交叉项部分, 以达到Gabor-WVD变换在保持良好的时频聚集特性的同时, 具有消除交叉项的效果. Gabor-WVD变换定义了$G(t, f)$与$W(t, f)$这两个过程变量, 其任意函数表达式如式(6)所示[18].
$ \begin{align} &GW(t, f)=p(G(t, f), W(t, f)) \end{align} $
(6) 其中, $p(x, y)$为任意函数, 例如当$p(x, y)=x^{a}y^{b}$时, $GW(t, f)=G^{a}(t, f)W^{b}(t, f)$; 当$p(x, y)=x+y$时, $GW(t, f)=G(t, f)+W(t, f)$. Gabor-WVD得出的结果可以很好地描述非平稳信号的二维时频分布, 并且具有良好的聚集性, 对于WVD产生的交叉项也有较好的抑制.但本文作者在文献[18]中已经明确, STFT-WVD不能真实反映信号的三维幅值, 如本文第三部分的数值实验仿真例子所示. Gabor-WVD的三维输出时频表示存在较大的误差和干扰, 具体体现在信息项上的交叉项无法消除.针对交叉项和三维幅值不能正确显示的问题, 需要对Gabor-WVD进一步改进.由于交叉项成因与变换过程中的两个乘项表现为时频轴上的频率相邻与中间频段位置的量[18], 因此, 消除交叉项的思路, 可采用设定阈值来清除指定频段上的虚假分量, 或是利用其他合适的数值进行替代, 或是对该项增添幂指数进行消除.实现Gabor-WVD的若干种方法如式(7)~(9).
$ \begin{align} &GW(t, f)={\rm min}\{ G(t, f), \vert W(t, f)\vert\} \end{align} $
(7) $ \begin{align} &GW(t, f)=W(t, f)\{\vert G(t, f)\vert >c\} \end{align} $
(8) $ \begin{align} &GW(t, f)={G^a(t, f)}{W^b(t, f)} \end{align} $
(9) 根据函数$p(x, y)$形式的不同, $GW(t, f)$的运算方式也有所不同.由于式(7)只取Gabor变换与WVD后的数值中的较小值, 故称为最小值法.式(8)将$c$设置为交叉项消除阈值, 对Gabor数组中的部分特定数据取0或1, 故称为二值化法.式(9)设置$a$、$b$为幂指数, 通过幂指数来消除交叉项, 故称为幂调节系数法.
2.3 BGabor-NWVD和BGabor-NSPWVD
尽管Gabor-WVD方法能够有效消除WVD产生的虚假分量, 但并不能消除叠加在信息项上的分量.针对此问题, 本文对Gabor和WVD同时进行改进, 提出BGabor-NWVD和BGabor-NSPWVD算法. BGabor-NWVD算法的步骤为:
步骤1. 首先各自对输入信号进行Gabor和WVD变换, 分别得到数组$A$和数组$B$.
步骤2. 根据时域波形的幅度, 对Gabor数组进行二值化处理(Binarization), 得到数组$Gabor\_ {1}$, 同时, 根据WVD绝对值数组中的最大值对WVD数组进行归一化处理(Normalization), 得到数组$WVD\_{1}$.
步骤3. 两数组进行点除, 得到新的临时数组$Y$, 根据其本身的大小, 设置阈值$K_1$、$K_2$.
步骤4. 修正点除后的结果数组$B$, 记录数组$Y$中大于$K_1$的位置, 在数组$B$中将相同位置的元素置0, 得到数组$B{'}$.
步骤5. 对数组$Y$进行修正, 将数组$Y$中大于$K_1$的元素置1, 将数组$Y$中小于$K_2$的元素也置1, 得到新数组$Y{'}$.
步骤6. 数组$B{'}$点除$Y{'}$, 输出BGabor-NWVD数组, 算法结束.
数值实验结果表明, 这一改进算法对信息项上交叉项的叠加有较好的抑制效果. SPWVD是对WVD的一种加窗平滑改进, SPWVD本身即具有消除WVD交叉项的作用, 将BGabor与SPWVD相结合, 也取得了较好的效果.需要特别说明的是, 与SPWVD结合后, 尽管牺牲了一定的二维时频聚集度, 但大幅度提高了非平稳信号三维时频表示的准确度. BGabor-NSPWVD算法流程与BGabor-NWVD算法类似, 将信号的WVD换成SPWVD即可, BGabor-NWVD变换算法流程框图如图 2所示.
3. 数值仿真实验
3.1 时频分布对比
本节采用2个构造函数来进行多个时频分析方法的测试实验, 重点讨论每种方法的时间和频率聚集度, 以及每个频率分量的幅值分布是否正确.构造的2个函数属于易产生交叉项的多频分量信号, 如式(10)和式(11), 分别为较为复杂的四分量线性调频信号和三分量正弦频率信号, 数值实验的方法采用STFT、WVD、SPWVD、Gabor-WVD、Gabor-SPWVD、BGabor-NWVD和BGabor-NSPWVD.图 3为三分量和四分量构造函数的二维时频图, 采用本文改进后的BGabor-NWVD和BGabor-NSPWVD, 可以看出, 这两种方法都能准确地显示函数的时频分布, 没有交叉项的干扰, 且时间-频率轴都有较好的聚集度.在进行BGabor-NWVD和BGabor-NSPWVD处理时, 同时比较使用最小值法、二值化法和幂指数调节法实现Gabor-WVD和Gabor-SPWVD的结果, 讨论各方法的三维幅值分布是否准确.线性调频四频率分量信号表示如式(10)所示, $m$为一般参数, 可以根据需要调整, 这里选取$m=0.2.$
图 3 四分量$f_2$和三分量$f_3$的二维时频图Fig. 3 Two-dimensional time-frequency diagram of four components signal $f_2$ and three components signal $f_3$$ \begin{align} f_2(t)=\sum\limits_{j=0}^3 {\rm e}^{(0.06+0.12j)t+j2\pi (\frac{1}{2} mt^2)} \end{align} $
(10) 正弦三频率分量信号如式(11)所示.
$ \begin{align} f_3(t)=\sum\limits_{j=0}^2 {\rm cos}((0.3+0.2j)\pi\cdot (t-10{\rm sin}(0.3\pi t))) \end{align} $
(11) 图 3中, BGabor-NWVD和BGabor-NSPWVD方法都能较好地实现复杂信号的二维时频分布, 消除了交叉项的干扰, BGabor-NWVD较BGabor-NSPWVD的二维时频聚集度略好, 采用NSPWVD方法牺牲了少许的平面聚集性, 但这并不影响时频表现和分布的判断.下面讨论这些方法的时频三维分布情况, 首先看基本的Gabor、WVD及WVD的改进SPWVD, 如图 4所示.
图 4 四分量$f_2$的三维时频图Fig. 4 Three-dimensional time-frequency diagram of four components signal $f_2$图 4为函数$f_2$的Gabor、WVD和SPWVD方法得到三维时频分布图, 其中, 图 4(a)中的Gabor方法没有交叉项, 但是底端呈扩散状分布, 频率聚集性差; 图 4(b)中WVD的交叉项干扰严重; 图 4(c)中SPWVD没有交叉项, 但是频率项上的幅度有尖端失真, 需要予以消除.
3.2 BGabor-NWVD
针对图 4出现的问题, 采用二值化改进的Gabor与WVD、SPWVD相结合的算法, 时频分析三维效果如图 5所示.
图 5 基于Gabor和WVD的四分量$f_2$(上)和三分量$f_3$ (下)的三维时频比较图Fig. 5 Three-dimensional time-frequency diagram of four-components $f_2$ (upper) and three-components $f_3$ (bottom) based on Gabor and WVD将Gabor与WVD直接相结合, 按照式(7)~(9)分别进行Gabor-WVD最小值法、二值化法和幂系数调节法的三维时频分布画图, 如图 5(a)、5(b)和5(c)所示, 3个图中函数$f_2$和$f_3$的频率分量幅度出现了较大失真.由此需要对Gabor-WVD算法进一步改进, 先将Gabor二值化(BGabor), 再结合归一化的WVD (NWVD), 得到BGabor-NWVD算法.实验仿真函数$f_2$和$f_3$的三维时频分布如图 5(d), 较好地克服了Gabor-WVD方法的幅值失真问题, 其幅值顶端依旧存在模糊现象, 还需要进一步加以改进.
3.3 BGabor-NSPWVD
SPWVD具有较好地平滑效果, 考虑将Gabor与SPWVD相结合, 可以进一步优化幅值模糊现象, 如图 6所示.
图 6 基于Gabor和SPWVD的四分量$f_2$(上)和三分量$f_3$ (下)的三维时频比较图Fig. 6 Three-dimensional time-frequency diagram of four-components $f_2$ (upper) and three-components $f_3$ (bottom) based on Gabor and SPWVDGabor与SPWVD相结合, 可以很大程度上平滑信号的幅度数值, 使其时频表现更接近真实的各频率分量幅度.但图 6中, 图 6(a)、6(b)和6(c)的幅度仍然存在一些干扰和噪声分量, 其中, 图 6(a)幅度顶端较好, 但是幅度的底端范围较宽, 出现能量泄漏现象. 图 6(d)的BGabor-NSPWVD算法效果最佳, 幅度刻画良好, 仿真函数$f_2$和$f_3$各分量的顶端和低端没有出现明显失真, 时间-频率-幅度分布清晰准确. BGabor-NSPWVD和Gabor-NWVD算法都能去除信息项上的叠加残余, 对交叉项有较好地抑制效果.两种算法相比较, BGabor-NSPWVD对信号分量在时频平面的幅度还原效果更优于Gabor-NWVD, 其时间-频率-幅度的三维分布具有良好的表现能力.
4. 鲁棒性实验
4.1 二维时频分布对比
BGabor-NWVD和BGabor-NSPWVD算法不仅对复杂信号的交叉项有较理想的消除作用, 而且对附着在信号上的高斯白噪声同样有较强的抑制效果, 两种算法在时频分析处理效果上具有强鲁棒性的特点, 利于处理实际的含噪信号.构造具有多频分量的实验仿真信号, 加入-10dB至20dB的高斯白噪声, 对比Gabor、WVD、SPWVD、Gabor-WVD、Gabor-SPWVD、BGabor-NWVD和BGabor-NSPWVD方法的二维时频分析效果.由于篇幅限制, 本文仅列出加入2dB噪声时函数$f_4$的二维时频分布, 如图 7所示.
图 7 含噪信号$f_4$的二维时频分布比较(SNR = 2dB)Fig. 7 The two-dimensional time-frequency distribution of the noisy signal $f_4$ (SNR = 2dB)图 7(b)为改进前的Gabor方法, 即高斯变换核STFT, 时频聚集度较差且噪声分布干扰严重.图 7(c)的WVD方法存在严重的交叉项干扰, 需要去除存在于真实频率之间的虚假频率分量.图 7(d)的SPWVD较好地去除了交叉项和噪声干扰, 但是由于采用伪平滑的手段, 使得该方法的时频聚集度较差.图 7(e)、7(f)和7(g)是结合后的Gabor-WVD方法, 其中幂系数调节法去除交叉项的效果略好, 但是这三种方法对噪声的抑制和去除交叉项都不满足要求.图 7(h)采用进一步结合后的BGabor-NWVD方法, 较好地去除了交叉项和噪声干扰, 具有最佳的二维时频分布表现.
图 7(i)、7(j)和7(k)是Gabor与SPWVD相结合的三种方法, 对交叉项的抑制尚可, 但对噪声干扰的去除效果还需提高. 图 7(l)为BGabor-NSPWVD方法, 对交叉项的抑制和噪声的去除都较为理想, 时频聚集度略逊于图 7(h)的BGabor-NWVD方法, 但是对交叉项的抑制和鲁棒性而言, BGabor-NSPWVD和BGabor-NWVD都是可以采用的方法.
4.2 时频聚集度评价
文献[20]给出了评价时频聚集度的量化公式, 如式(12), $n$为时间窗长度, $\omega$为频率, 式(12)通过调整窗口参数, 求得最大的时频分布的第四幂范数与第二幂范数之商来评价时频聚集度.向实验函数$f_4$添加-10dB至20dB的高斯白噪声, 求解函数$f_4$的$E_{JP}$数值, 用以衡量各个方法的聚集度分布, 如表 1所示.将表 1的数值用图 8的折线趋势来表示, 能够较为直观地比较每种方法的聚集度效果.
图 8 时频聚集度参数$E_{JP}$评价比较Fig. 8 Comparison of time-frequency aggregation degree evaluation on $E_{JP}$表 1 仿真函数$f_4$在不同噪声条件下各方法的聚集度$E_{JP}$数值比较Table 1 The $E_{JP}$ numerical comparison of experimental function $f_4$ in different noise conditions不同方法 SNR = -10 SNR = -5 SNR = 0 SNR = 5 SNR = 10 SNR = 15 SNR = 20 Gabor $9.35\, \times\, 10^{-7}$ $2.03\, \times\, 10^{-6}$ $5.71\, \times\, 10^{-6}$ $1.22\, \times\, 10^{-5}$ $1.65\, \times\, 10^{-5}$ $1.87\, \times\, 10^{-5}$ $1.96\, \times\, 10^{-5}$ WVD $3.10\, \times\, 10^{-6}$ $3.50\, \times\, 10^{-6}$ $8.21\, \times\, 10^{-6}$ $2.62\, \times\, 10^{-5}$ $4.74\, \times\, 10^{-5}$ $5.95\, \times\, 10^{-5}$ $6.59\, \times\, 10^{-5}$ SPWVD $8.35\, \times\, 10^{-7}$ $5.26\, \times\, 10^{-6}$ $7.55\, \times\, 10^{-6}$ $1.18\, \times\, 10^{-5}$ $1.21\, \times\, 10^{-5}$ $1.27\, \times\, 10^{-5}$ $1.27\, \times\, 10^{-5}$ Gabor-WVD (二值化) $4.61\, \times\, 10^{-6}$ $1.21\, \times\, 10^{-5}$ $1.00\, \times\, 10^{-4}$ $1.22\, \times\, 10^{-4}$ $1.19\, \times\, 10^{-4}$ $1.19\, \times\, 10^{-4}$ $1.24\, \times\, 10^{-4}$ Gabor-WVD (幂系数调节) $1.83\, \times\, 10^{-6}$ $2.85\, \times\, 10^{-6}$ $1.12\, \times\, 10^{-5}$ $3.06\, \times\, 10^{-5}$ $5.43\, \times\, 10^{-5}$ $6.33\, \times\, 10^{-5}$ $6.68\, \times\, 10^{-5}$ Gabor-WVD (最小值) $2.76\, \times\, 10^{-6}$ $3.35\, \times\, 10^{-6}$ $7.43\, \times\, 10^{-6}$ $2.43\, \times\, 10^{-5}$ $5.92\, \times\, 10^{-5}$ $7.92\, \times\, 10^{-5}$ $8.85\, \times\, 10^{-5}$ Gabor-SPWVD (二值化) $2.24\, \times\, 10^{-6}$ $1.23\, \times\, 10^{-5}$ $3.43\, \times\, 10^{-5}$ $4.78\, \times\, 10^{-5}$ $4.77\, \times\, 10^{-5}$ $4.85\, \times\, 10^{-5}$ $4.86\, \times\, 10^{-5}$ Gabor-SPWVD (幂系数调节) $1.76\, \times\, 10^{-6}$ $3.75\, \times\, 10^{-6}$ $1.14\, \times\, 10^{-5}$ $2.04\, \times\, 10^{-5}$ $2.70\, \times\, 10^{-5}$ $3.07\, \times\, 10^{-5}$ $3.22\, \times\, 10^{-5}$ Gabor-SPWVD (最小值) $2.02\, \times\, 10^{-6}$ $4.04\, \times\, 10^{-6}$ $1.49\, \times\, 10^{-5}$ $3.13\, \times\, 10^{-5}$ $3.81\, \times\, 10^{-5}$ $3.94\, \times\, 10^{-5}$ $3.99\, \times\, 10^{-5}$ BGabor-NWVD $6.73\, \times\, 10^{-6}$ $4.37\, \times\, 10^{-5}$ $1.47\, \times\, 10^{-4}$ $1.53\, \times\, 10^{-4}$ $1.42\, \times\, 10^{-4}$ $1.39\, \times\, 10^{-4}$ $1.44\, \times\, 10^{-4}$ BGabor-NSPWVD $3.48\, \times\, 10^{-6}$ $3.26\, \times\, 10^{-5}$ $4.66\, \times\, 10^{-5}$ $6.16\, \times\, 10^{-5}$ $6.11\, \times\, 10^{-5}$ $6.08\, \times\, 10^{-5}$ $6.03\, \times\, 10^{-5}$ $ \begin{align} E_{JP}=\frac{{\sum\limits_{n=0}^{N-1}\sum\limits_{\omega =0}^{W-1}\vert Q(n, \omega)\vert ^4}}{\left(\sum\limits_{n=0}^{N-1}\sum\limits_{\omega =0}^{W-1}\vert Q(n, \omega)\vert ^2\right)^2} \end{align} $
(12) 时频聚集度评价图 8中, 信噪比(SNR)大于0dB时, BGabor-NWVD和BGabor-NSPWVD方法的聚集度$E_{JP}$开始平稳, 且其值高于其他方法, 随着SNR的变化, 这两种改进算法皆具有较高的时频聚集度和强鲁棒性, BGabor-NWVD的聚集度优于BGabor-NSPWVD算法, 与图 7的二维时频分布效果图一致.
5. 硬质合金顶锤破裂信号监测
在人造金刚石合成加工过程中, 硬质合金顶锤在交变的外载应力及热应力作用下, 易发生弹性形变和塑性形变.如由于疲劳损坏而产生微观裂纹, 这些裂纹不断发展贯穿许多晶粒成为宏观裂纹, 使顶锤断面进而发生横向或纵向的破裂或压溃[22].作为高脆性材料, 受材料自身特点的正常性破坏和结构、工艺、人为操作等非正常性破坏等因素影响, 硬质合金顶锤对于微观裂纹极为敏感, 当产生微观裂纹时, 由裂纹成核、裂纹扩展延伸并迅速发展成宏观裂纹, 产生强烈的声裂发射频率段信号[23].大量实验研究表明, 声发射源主要有塑性形变(滑移和孪生)、断裂(裂纹的形成和扩展、第二相质点或夹杂物)、相变(马氏体相变、共晶反映等)、磁效应和表面效应等[22].图 9(a)为硬质合金顶锤的工作现场图, 图 9(b)为硬质合金顶锤的三维模拟效果图.
图 9 硬质合金顶锤的现场实物图与三维模拟图Fig. 9 Carbide anvil physical site map and 3D simulation figure通过检测硬质合金顶锤工作时的疑似破裂信号, 对比正常状态下信号频率特性及裂纹产生时的信号频率特征, 判别顶锤破裂的发生与否.目前存在的难点是微小破裂过程中释放的应力能太小, 声发射信号相当微弱, 在顶锤工作恶劣环境下, 受复杂噪声信号干扰, 难以有效地设置采集频率窗口和区分破裂频率成分.为有效地接收到破裂信号, 需要传感器设置合适的频率检测阈值, 可采用时频分析的手段来分析疑似破裂信号, 为判断是否破裂提供参考依据.本文利用BGabor-NWVD和BGabor-NSPWVD算法对金属破裂疑似样本信号进行分析, 获取其时频联合分布特点, 找到疑似金属破裂信号段的时频分布表示.
图 10中, 图 10(a)、10(d)、10(g)为疑似破裂样本信号的时域波形, 图 10(b)、10(e)、10(h)为采用BGabor-NWVD算法得到的时频分析结果, 图 10(c)、10(f)、10(i)为采用BGabor-NSPWVD算法得到的时频分析结果.图 10(c)、10(f)、10(i)的时频集中区域较图 10(b)、10(e)、10(h)更为明显, 在疑似破裂发生时段, 3个样本信号的频率分布集中范围, 大致分布在100kHz这个较为明显的频率区域.根据BGabor-NSPWVD算法得到的金属破裂频率窗口, 将其设置为传感器的破裂判断阈值, 用以自动在线监测硬质合金顶锤工作时是否发生疑似破裂, 提高了判别几率, 取得了较好的效果.
6. 结论
针对STFT、Gabor和WVD出现的时频分辨率模糊和存在交叉项等缺点, 以及一些结合算法如STFT-WVD和Gabor-WVD出现的三维幅度失真, 抗噪性能及鲁棒性能还不理想的问题, 本文提出BGabor-NWVD和BGabor-NSPWVD算法.通过对复杂线性调频信号和多分量的正弦频率分量信号进行数值仿真实验, BGabor-NWVD和BGabor-NSPWVD算法在抑制了交叉项的同时, 具有较高锐化时频分辨率, 两种算法的抗噪性能和鲁棒性也较为理想.由四分量线性调频信号和三分量正弦信号的数值仿真实验可知, BGabor-NWVD的二维时频表示优于BGabor-NSPWVD, BGabor-NSPWVD的三维时频表示优于BGabor-NWVD. BGabor-NWVD和BGabor-NSPWVD算法综合了STFT、Gabor、WVD和SPWVD各自的频率自适应性和良好的时频表示, 其时频聚集度评价参数$E_{JP}$高于其他方法, 具有高锐化频率聚集度优点的同时, 能够真实还原信号频率分量的幅度.通过硬质合金顶锤工作时产生的疑似破裂样本信号进行时频分析, 本文方法可以较为准确地寻找传感器的频率判别窗口, 为金属破裂监测设备数据采集卡提供有效的阈值参考.
-
图 1 信号$f_1$的理想时频、Gabor和WVD对比图
Fig. 1 Ideal time-frequency spectrum, Gabor, WVD of $f_1$
图 3 四分量$f_2$和三分量$f_3$的二维时频图
Fig. 3 Two-dimensional time-frequency diagram of four components signal $f_2$ and three components signal $f_3$
图 4 四分量$f_2$的三维时频图
Fig. 4 Three-dimensional time-frequency diagram of four components signal $f_2$
图 5 基于Gabor和WVD的四分量$f_2$(上)和三分量$f_3$ (下)的三维时频比较图
Fig. 5 Three-dimensional time-frequency diagram of four-components $f_2$ (upper) and three-components $f_3$ (bottom) based on Gabor and WVD
图 6 基于Gabor和SPWVD的四分量$f_2$(上)和三分量$f_3$ (下)的三维时频比较图
Fig. 6 Three-dimensional time-frequency diagram of four-components $f_2$ (upper) and three-components $f_3$ (bottom) based on Gabor and SPWVD
图 7 含噪信号$f_4$的二维时频分布比较(SNR = 2dB)
Fig. 7 The two-dimensional time-frequency distribution of the noisy signal $f_4$ (SNR = 2dB)
图 8 时频聚集度参数$E_{JP}$评价比较
Fig. 8 Comparison of time-frequency aggregation degree evaluation on $E_{JP}$
图 9 硬质合金顶锤的现场实物图与三维模拟图
Fig. 9 Carbide anvil physical site map and 3D simulation figure
表 1 仿真函数$f_4$在不同噪声条件下各方法的聚集度$E_{JP}$数值比较
Table 1 The $E_{JP}$ numerical comparison of experimental function $f_4$ in different noise conditions
不同方法 SNR = -10 SNR = -5 SNR = 0 SNR = 5 SNR = 10 SNR = 15 SNR = 20 Gabor $9.35\, \times\, 10^{-7}$ $2.03\, \times\, 10^{-6}$ $5.71\, \times\, 10^{-6}$ $1.22\, \times\, 10^{-5}$ $1.65\, \times\, 10^{-5}$ $1.87\, \times\, 10^{-5}$ $1.96\, \times\, 10^{-5}$ WVD $3.10\, \times\, 10^{-6}$ $3.50\, \times\, 10^{-6}$ $8.21\, \times\, 10^{-6}$ $2.62\, \times\, 10^{-5}$ $4.74\, \times\, 10^{-5}$ $5.95\, \times\, 10^{-5}$ $6.59\, \times\, 10^{-5}$ SPWVD $8.35\, \times\, 10^{-7}$ $5.26\, \times\, 10^{-6}$ $7.55\, \times\, 10^{-6}$ $1.18\, \times\, 10^{-5}$ $1.21\, \times\, 10^{-5}$ $1.27\, \times\, 10^{-5}$ $1.27\, \times\, 10^{-5}$ Gabor-WVD (二值化) $4.61\, \times\, 10^{-6}$ $1.21\, \times\, 10^{-5}$ $1.00\, \times\, 10^{-4}$ $1.22\, \times\, 10^{-4}$ $1.19\, \times\, 10^{-4}$ $1.19\, \times\, 10^{-4}$ $1.24\, \times\, 10^{-4}$ Gabor-WVD (幂系数调节) $1.83\, \times\, 10^{-6}$ $2.85\, \times\, 10^{-6}$ $1.12\, \times\, 10^{-5}$ $3.06\, \times\, 10^{-5}$ $5.43\, \times\, 10^{-5}$ $6.33\, \times\, 10^{-5}$ $6.68\, \times\, 10^{-5}$ Gabor-WVD (最小值) $2.76\, \times\, 10^{-6}$ $3.35\, \times\, 10^{-6}$ $7.43\, \times\, 10^{-6}$ $2.43\, \times\, 10^{-5}$ $5.92\, \times\, 10^{-5}$ $7.92\, \times\, 10^{-5}$ $8.85\, \times\, 10^{-5}$ Gabor-SPWVD (二值化) $2.24\, \times\, 10^{-6}$ $1.23\, \times\, 10^{-5}$ $3.43\, \times\, 10^{-5}$ $4.78\, \times\, 10^{-5}$ $4.77\, \times\, 10^{-5}$ $4.85\, \times\, 10^{-5}$ $4.86\, \times\, 10^{-5}$ Gabor-SPWVD (幂系数调节) $1.76\, \times\, 10^{-6}$ $3.75\, \times\, 10^{-6}$ $1.14\, \times\, 10^{-5}$ $2.04\, \times\, 10^{-5}$ $2.70\, \times\, 10^{-5}$ $3.07\, \times\, 10^{-5}$ $3.22\, \times\, 10^{-5}$ Gabor-SPWVD (最小值) $2.02\, \times\, 10^{-6}$ $4.04\, \times\, 10^{-6}$ $1.49\, \times\, 10^{-5}$ $3.13\, \times\, 10^{-5}$ $3.81\, \times\, 10^{-5}$ $3.94\, \times\, 10^{-5}$ $3.99\, \times\, 10^{-5}$ BGabor-NWVD $6.73\, \times\, 10^{-6}$ $4.37\, \times\, 10^{-5}$ $1.47\, \times\, 10^{-4}$ $1.53\, \times\, 10^{-4}$ $1.42\, \times\, 10^{-4}$ $1.39\, \times\, 10^{-4}$ $1.44\, \times\, 10^{-4}$ BGabor-NSPWVD $3.48\, \times\, 10^{-6}$ $3.26\, \times\, 10^{-5}$ $4.66\, \times\, 10^{-5}$ $6.16\, \times\, 10^{-5}$ $6.11\, \times\, 10^{-5}$ $6.08\, \times\, 10^{-5}$ $6.03\, \times\, 10^{-5}$ 
下载: 导出CSV
-
[1] Volaric I, Sucic V, Stankovic S. A data driven compressive sensing approach for time-frequency signal enhancement. Signal Processing, 2017, 141:229-239 doi: 10.1016/j.sigpro.2017.06.013 [2] Daubechies I. The wavelet transform, time-frequency localization and signal analysis. IEEE Transactions on Information Theory, 1990, 36(5):961-1005 doi: 10.1109/18.57199 [3] Kayhan A S. Difference equation representation of chirp signals and instantaneous frequency/amplitude estimation. IEEE Transactions on Signal Processing, 1996, 44(12):2948-2958 doi: 10.1109/78.553470 [4] 储昭碧, 张崇巍, 冯小英.基于自适应滤波器的时频分析.自动化学报, 2009, 35(11):1420-1428 http://www.aas.net.cn/CN/abstract/abstract13594.shtmlChu Zhao-Bi, Zhang Chong-Wei, Feng Xiao-Ying. Adaptive filter-based time-frequency analysis. Acta Automatica Sinica, 2009, 35(11):1420-1428 http://www.aas.net.cn/CN/abstract/abstract13594.shtml [5] 杨默涵, 陈万忠, 李明阳.基于总体经验模态分解的多类特征的运动想象脑电识别方法研究.自动化学报, 2017, 43(5):743-752 http://www.aas.net.cn/CN/abstract/abstract19052.shtmlYang Mo-Han, Chen Wan-Zhong, Li Ming-Yang. Multiple feature extraction based on ensemble empirical mode decomposition for motor imagery EEG recognition tasks. Acta Automatica Sinica, 2017, 43(5):743-752 http://www.aas.net.cn/CN/abstract/abstract19052.shtml [6] Cohen L. Time-frequency distributions-a review. Proceedings of the IEEE, 1989, 77(7):941-981 doi: 10.1109/5.30749 [7] Gabor D. Theory of communication. Journal of the Institution of Electrical Engineers, 1946, 93:429-549 http://d.old.wanfangdata.com.cn/OAPaper/oai_doaj-articles_20930838401474143ef29e11d60c01a5 [8] Mansinha L, Stockwell R G, Lowe R P. Pattern analysis with two-dimensional spectral localisation:applications of two-dimensional S transforms. Physica A:Statistical Mechanics and Its Applications, 1997, 239(1-3):286-295 doi: 10.1016/S0378-4371(96)00487-6 [9] Daubechies I. Ten Lectures on Wavelets. Philadelphia, Pa, USA:SIAM, 1992. 39-72 [10] Huang N E, Shen Z, Long S R, Wu M C, Shih H H, Zheng Q, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 1998, 454(1971):903-995 doi: 10.1098/rspa.1998.0193 [11] Wu Z H, Huang N E. Ensemble empirical mode decomposition:a noise-assisted data analysis method. Advances in Adaptive Data Analysis, 2011, 1(1):1-41 http://d.old.wanfangdata.com.cn/Periodical/dianzixb201305033 [12] Zhang J, Yan R Q, Gao R X, Feng Z H. Performance enhancement of ensemble empirical mode decomposition. Mechanical Systems and Signal Processing, 2010, 24(7):2104-2123 doi: 10.1016/j.ymssp.2010.03.003 [13] Daubechies I, Lu J F, Wu H T. Synchrosqueezed wavelet transforms:an empirical mode decomposition-like tool. Applied and Computational Harmonic Analysis, 2011, 30(2):243-261 doi: 10.1016/j.acha.2010.08.002 [14] Zhou W X, Shao Z F, Diao C Y, Cheng Q M. High-resolution remote-sensing imagery retrieval using sparse features by auto-encoder. Remote Sensing Letters, 2015, 6(10):775-783 doi: 10.1080/2150704X.2015.1074756 [15] Boashash B, Black P. An efficient real-time implementation of the Wigner-Ville distribution. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1987, 35(11):1611-1618 doi: 10.1109/TASSP.1987.1165070 [16] 王见, 李金同, 卢华玲, 尹爱军.采用STFT-Wigner变换抑制Wigner-Ville分布交叉项.重庆大学学报, 2013, 36(8) 15-18, 25 http://d.old.wanfangdata.com.cn/Periodical/cqdxxb201308003Wang Jian, Li Jin-Tong, Lu Hua-Ling, Yin Ai-Jun. Using STFT-Wigner transform to suppress the cross terms in Wigner-Ville distribution. Journal of Chongqing University, 2013, 36(8) 15-18, 25 http://d.old.wanfangdata.com.cn/Periodical/cqdxxb201308003 [17] 付卫红, 王璐, 贾坤, 路贵朝.基于STFT与SPWVD的跳频参数盲估计算法.华中科技大学学报(自然科学版), 2014, 42(9):59-63 http://d.old.wanfangdata.com.cn/Periodical/hzlgdxxb201409013Fu Wei-Hong, Wang Lu, Jia Kun, Lu Gui-Chao. Blind parameter estimation algorithm for frequency hopping signals based on STFT and SPWVD. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2014, 42(9):59-63 http://d.old.wanfangdata.com.cn/Periodical/hzlgdxxb201409013 [18] 郝国成, 陈忠昌, 赵娟, 曾佐勋, 刘辉, Sibgatulin V G, 等.基于NSTFT-WVD变换的芦山slMs7.0级地震前后地球天然脉冲电磁场信号时频分析.地学前缘, 2016, 23(1) 276-286Hao Guo-Cheng, Chen Zhong-Chang, Zhao Juan, Zeng Zuo-Xun, Liu Hui, Sibgatulin V G, et al. Time-frequency analysis of the Earth's natural pulse electromagnetic field signal before and after the Lushan slMs7.0 earthquake based on NSTFT-WVD transform. Earth Science Frontiers, 2016, 23(1) 276-286 [19] Zhou J, Fang X Y, Tao L. A sparse analysis window for discrete Gabor transform. Circuits, Systems, and Signal Processing, 2017, 36(10):4161-4180 doi: 10.1007/s00034-017-0510-0 [20] Jones D L, Parks T W. A high resolution data-adaptive time-frequency representation. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1990, 38(12):2127-2135 doi: 10.1109/29.61539 [21] Allen J. Short term spectral analysis, synthesis, and modification by discrete Fourier transform. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1977, 25(3):235-238 doi: 10.1109/TASSP.1977.1162950 [22] 王民, 王松顺.硬质合金顶锤和压缸的破裂特征、破坏因素的分析和提高使用寿命的途径.硬质合金, 1994, 11(4):231-240 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK199400560973Wang Min, Wang Song-Shun. The breakage characteristic and breakdown factors analysis of cemented carbide anvil and pressure die and the way of raise their sevice life. Gemented Carbide, 1994, 11(4):231-240 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK199400560973 [23] Ibrahim M E, Zhuang W Z L. Nondestructive inspection of fatigue crack propagation beneath supersonic particle deposition coatings during fatigue testing. International Journal of Fatigue, 2017, 102:149-157 doi: 10.1016/j.ijfatigue.2017.05.005 期刊类型引用(9)
1. 滕召胜，梁成斌，唐求，张雷鹏，成达. K-S变换及其电网超谐波时频分析应用. 湖南大学学报(自然科学版). 2024(06): 128-136 . 
百度学术2. 杨露露，孙文，孙薇，孙维，张译航，李艺. 输油泵故障智能诊断研究综述. 石油化工设备. 2024(05): 55-61 . 百度学术3. Liangliang Dong，Qian Xiao，Yanjie Jia，Tianhai Fang. Review of research on intelligent diagnosis of oil transfer pump malfunction. Petroleum. 2023(02): 135-142 . 必应学术4. 郝国成，张必超，锅娟，张雅冰，石光耀，王盼盼，张薇. 高质量LMSCT时频分析算法及其在雷达信号目标检测中的应用. 上海交通大学学报. 2022(02): 231-241 . 百度学术5. 郝国成，冯思权，王巍，凌斯奇，谭淞元. 基于NGWarblet-WVD的高质量时频分析方法. 自动化学报. 2022(10): 2526-2536 . 本站查看6. 张帆，杨晓忠，吴立飞，韩晓明，王树波. 基于短时傅里叶变换和卷积神经网络的地震事件分类. 地震学报. 2021(04): 463-473+533 . 百度学术7. 江莉，尚文擎，周军妮，卫铭斐，王燕妮. 一种用于地震信号分析的二阶挤压小波变换算法. 电子与信息学报. 2021(12): 3710-3717 . 百度学术8. 邬蒙蒙. 改进的完备经验模态分解与WVD变换相结合在油气检测中的应用. 油气地球物理. 2020(01): 26+73-76 . 百度学术9. 张海强，王楠，刘松，李柄奎. 基于参数优化VMD和SPWVD的轨道波磨辨识方法. 铁路计算机应用. 2020(06): 18-24 . 百度学术其他类型引用(2)
-
下载:
计量
- 文章访问数: 3267
- HTML全文浏览量: 347
- PDF下载量: 387
- 被引次数: 11
