K-S变换及其电网超谐波时频分析应用

K-S Transformation and Its Application to Superharmonic Time-Frequency Analysis in Power System

摘要

依据FFT→优化窗→IFFT思路，突破线性时频变换的窗函数积分性能桎梏，实现高性能优化窗函数的线性时频变换应用，建立新型时频变换算法——K-S变换.对信号x(t)的FFT频谱向量进行频移处理后，与该频移点下Kaiser优化窗的频谱向量进行Hadamard乘积，再将乘积结果进行FFT逆变换（IFFT），构造出K-S变换复时频矩阵，由此获得x(t)的时间-频率-幅值、时间-频率-相位三维信息；给出逆变换的数学推导与局部性质、线性性质和变分辨率特性；0～150 kHz电网的稳态与时变超谐波信号仿真实验表明，K-S变换的时域、频域分辨能力均优于流行的短时傅里叶变换、S变换，具有优良的变分辨率性能；0～40 kHz超谐波信号的实测证明，基于K-S变换的超谐波电压幅值测量绝对误差均小于0.032 3 V.

Abstract

According to the idea of FFT（Fast Fourier Transform）→ optimization window → IFFT（Inverse Fast Fourier Transform）， the performance constraints of window function integration in linear time-frequency transformation are broken through， achieve the application of high-performance optimization window function in linear time-frequency transformation is achieved， and establish a new time-frequency transformation algorithm， K-S transformation is establisged. After frequency shifting for the FFT spectrum vector of signal x（t） is carried out， Hadamard multiplication with the spectrum vector of Kaiser optimization window at the frequency shift point， and then inverse FFT transformation （IFFT） on the product result is performed， the K-S transform complex time-frequency matrix is constructed to obtain three-dimensional time-frequency amplitude and time-frequency phase information of x（t）. The mathematical derivation， local properties， linear properties， and variable resolution properties of inverse transformation are given； simulation experiments on steady state and time-varying superharmonic signals from 0~150 kHz power grids show that the K-S transform has better resolution in both time and frequency domains than popular short-time Fourier transform and S-transform， and has excellent variable resolution performance. The actual measurement of 0~40 kHz superharmonic signals proves that the absolute error of superharmonic voltage amplitude measurement based on K-S transformation is less than 0.032 3 V.

Graphical abstract

关键词

K-S变换 / 时频分析 / Kaiser优化窗 / 变分辨率特性 / 电网超谐波

时频分析即时域、频域联合域分析，利用时间和频率联合函数，将一维时域信号映射到二维时频空间，以实现信号局部分析、瞬态特征刻画，是时变非平稳信号分析研究热点，越来越受到国内外重视^［1-2］.时频分析可获得复杂信号时间域和频率域联合分布特征信息，刻画信号幅值随频率变化和信号频率随时间变化的特性.如电网基波和谐波的频率、幅值、相位、起止时刻和持续时间等参数.快速且准确的时频分析方法是测试信号分析、故障诊断、模式识别、生物医学工程等领域工程应用迫切需要解决的科学问题^［3］.以电网监测为例，时频分析为电力系统谐波潮流计算、谐波电能计量、绝缘设备监测提供科学依据，对定位谐波源、治理电能质量扰动、确保电网安全经济运行具有重要意义^［4］.

根据时频分析联合函数表达形式，现有时频分析方法可分为线性、双线性、非线性三大类^［5］：1） 以短时傅里叶变换（Short-Time Fourier Transform，STFT）、小波变换（Wavelet Transform，WT）、S变换（S-Transform，ST）为代表的线性时频分析方法^［6-8］，计算量相对较小，是当前时频分析的主流方法，但时频聚集性能受限，时、频分辨率不能兼顾.2） Wigner-Ville及其改进分布等双线性时频分析方法^［9-10］，从能量角度刻画信号的时频分布，比取模值平方的线性时频分析方法更能刻画信号能量随时间、频率的变化，具有较高的时频聚集性能，但受交叉项干扰，动态信号分析的分辨率低，抑制交叉项则计算量很大，未能在高分辨率与抑制交叉项两方面兼顾. 3） Hilbert-Huang变换、信号稀疏分解等非线性时频分析方法^［11-12］，具有高时频聚集性能，可以获得较高分辨率，但计算量大、效率低，工程应用受限.

高时频聚集是时频分析的性能需求，卓越的变分辨率特性才能获得优良的时频聚集性能.但现有线性时频分析方法建立在窗函数无穷积分须为常数基础上，这一限制严重制约了变换算法的变分辨率性能.突破时频变换对窗函数面积的桎梏，是窗函数应用理论与时频分析方法亟待解决的问题.

本文在已有优化窗函数^［13］研究基础上，深化和发展快速傅里叶变换（Fast Fourier Transform，FFT）、STFT、WT、ST等线性时频分析思想，突破逆变换，形成高性能时频分析新方法——K-S变换，揭示变换算法的局部性质、线性性质和变分辨率特性，复杂电网动态超谐波仿真与测试证明变换算法的准确性.

1 K-S变换思路

FFT→STFT→WT→ST的递进思路，展示了线性时频分析方法不断优化的进展.在STFT与WT结合、改进基础上发展的ST^［14］，是目前国内外信号时频分析领域最热门的方法之一.ST凭借其在STFT、WT基础上构造的Gaussian窗可随被测信号频率自适应调节的优势，获得了大量应用^［8］.但STFT、ST推导逆变换的充要条件是时频变换算法中的窗函数无穷积分为常数，基于此才能通过FFT逆变换（Inverse Fast Fourier Transform，IFFT）对信号进行还原^［15］.而这种“窗函数面积限制”恰恰也是STFT、ST等时频分析方法未能获得卓越变分辨率特性的桎梏.

目前，各种广义ST均对其窗函数进行限制以获得逆变换，窗函数因此不能达到更佳的自适应时频调节状态^［15-16］.为实现窗函数更好的时频调节性能，部分学者探索了分段设计多组ST参数的改进方法，通过参数分段，提高特征频率下的时频变换性能.这种相互拼接的方式，在变换算法的自适应性和计算量方面显然有一定的局限性^［17］.

卓越的变分辨率特性需要卓越的窗函数性能作为保障，突破现有时频分析方法对窗函数无穷积分的桎梏，才能将性能更优越的窗函数应用于时频变换、取得更好的变分辨率特性.

2 K-S变换原理

Kaiser优化窗^［13］基于Kaiser窗主瓣与旁瓣的能量比例近乎最大、主瓣宽度和旁瓣高度之间的比重可自由选择、对信号加权灵活等优越性，宜于动态信号时频分析应用.

2.1 Kaiser窗及其优化

Kaiser窗函数如下^［13］：

w_k(t,f)=\frac{I_{\mathrm{0}}[\alpha (f)\sqrt[]{\mathrm{1}-{(t/T)}^{\mathrm{2}}}]}{I_{\mathrm{0}}[\alpha (f)]}, |t|\le T（1）

式中：\alpha (f)为频率f调节因子；I_{\mathrm{0}}(\cdot )为第一类零阶数Bessel函数，函数的幂级数展开为：

I_{\mathrm{0}}(\cdot )=\mathrm{1}+{\displaystyle \sum _{m=\mathrm{1}}^{\mathrm{\infty }}}\left[\frac{{(\cdot /\mathrm{2})}^m}{m!}\right]（2）

其N点采样的FFT频谱表示为：

W_k(\xi ,f)=\frac{N}{I_{\mathrm{0}}[\alpha (f)]}\frac{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\sqrt[]{{\alpha }^{\mathrm{2}}(f)-{(N\xi /\mathrm{2})}^{\mathrm{2}}})}{\sqrt[]{{\alpha }^{\mathrm{2}}(f)-{(N\xi /\mathrm{2})}^{\mathrm{2}}}}（3）

设c_{\mathrm{1}}、c_{\mathrm{2}}为比例常数，定义调节因子：

\alpha (f)=\frac{c_{\mathrm{1}}}{f+c_{\mathrm{2}}}（4）

即定义窗宽随被测信号频率f的调节机制.这一定义下的频率调节因子\alpha (f)随频率f的调节具有自适应能力，形成依据被测信号频率特征改变参数特性的Kaiser优化窗.在实际计算中，Kaiser优化窗调节因子式（4）的离散表达式为：

\alpha \left[f\right]=\frac{c_{\mathrm{1}}}{\left(Nf/f_{\mathrm{s}}\right)+c_{\mathrm{2}}}（5）

式中：(Nf/f_{\mathrm{s}})为实际计算中表示频率f的频率点；f_s为采样率；N为采样点数.

图1为Kaiser优化窗与目前线性时频分析应用最广泛的Gaussian窗^［9］的幅频特性仿真实验比较.其中选取控制参数c_{\mathrm{1}}=\mathrm{40}\mathrm{000}、c_{\mathrm{2}}=\mathrm{400}，采样率f_{\mathrm{s}}=\mathrm{409.6} kHz，采样长度N=\mathrm{4}\mathrm{096}.由图1可见.

1） Kaiser优化窗、Gaussian窗均可通过调节因子实现随频率调节变分辨率.

2） Kaiser优化窗的时频聚集性能优于Gaussian窗，主瓣随被测信号频率的增加而变窄，时频聚集性能更优，对高频信号有更佳的频率分辨率；反之，Gaussian窗频谱的主瓣随被测信号频率增加而变宽，时频聚集性能、频率分辨率逐步降低.

3） 高时频聚集性能，需要频率分辨率高的窗函数，即窗函数的频谱要有尽量窄的主瓣，Gaussian窗极端情况下是矩形窗，完全丧失了时间检测能力；而Kaiser优化窗在f=2 kHz处旁瓣峰值电平为-295.44 dB，能够保证一定时间分辨率，同时其频谱有较窄的主瓣，在宽频范围具有卓越的自适应变分辨率特性.

2.2 K-S变换模型

K-S变换定义如下：

K(\tau ,f)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(t)w_k(\tau -t,f){\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }ft}\mathrm{d}t

={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(t)\frac{I_{\mathrm{0}}\left[\alpha (f)\sqrt[]{\mathrm{1}-[{(\tau -t)/T]}^{\mathrm{2}}}\right]}{I_{\mathrm{0}}[\alpha (f)]}{\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }ft}\mathrm{d}t（6）

式中：x(t)\in L^{\mathrm{2}}(R)为被测信号；t为时间；f为频率；\tau为时移因子；w_k(t,f)为Kaiser优化窗.信号x(t)\in L^{\mathrm{2}}(R)的连续小波变换W(\tau ,g)为：

W(\tau ,g)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(t)w_{\mathrm{w}}(t-\tau ,g)\mathrm{d}t（7）

式中：w_{\mathrm{w}}(t,g)为母小波的伸缩时移变换；g为控制小波w_{\mathrm{w}}(t,g)宽度的尺度因子.

从算法原理上，K-S变换可以看作连续小波变换W(\tau ,g)的一种相位修正，即信号x(t)的K-S变换可表达为一相位因子与其连续小波变换的乘积：

K(\tau ,f)={\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }f\tau }W(\tau ,f)（8）

式中：小波母函数w_{\mathrm{w}}(t,f)为改进Kaiser优化窗与一个复向量的乘积.

w_w(t,f)=w_k(t,f){\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }ft}=

\frac{I_{\mathrm{0}}\left[\alpha (f)\sqrt[]{\mathrm{1}-{(t/\mathrm{T})}^{\mathrm{2}}}\right]}{I_{\mathrm{0}}[\alpha (f)]}{\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }ft}（9）

将式（7）和式（9）代入式（8）可得x(t)的K-S变换表达式为：

K(\tau ,f)={\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }f\tau }{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(t)w_{\mathrm{w}}(t-\tau ,f)\mathrm{d}t=

{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(t)\frac{I_{\mathrm{0}}\left[\alpha (f)\sqrt[]{\mathrm{1}-[{(t-\tau )/\mathrm{T}]}^{\mathrm{2}}}\right]}{I_{\mathrm{0}}[\alpha (f)]}{\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }ft}\mathrm{d}t（10）

式中：时移因子\tau为控制Kaiser优化窗在时间轴上的位置参数.

式（9）中的小波母函数无须满足小波变换中零均值容许性条件.据此，K-S变换可看作小波变换的一种自适应、智能化、简洁化扩展.

进一步，由卷积定理可获得K-S变换的FFT表达.对于信号x(t)\in L^{\mathrm{2}}(R)，K-S变换K(\tau ,f)由其FFT运算X(f)表示为：

K(\tau ,f)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}X(\xi +f)W_{\mathrm{k}}(\xi ,f){\mathrm{e}}^{\mathrm{j}\mathrm{2}\mathrm{\pi }\xi \tau }\mathrm{d}\xi（11）

式中：X(\xi +f)为信号x(t)的FFT频谱进行频移后的频谱，W_{\mathrm{k}}(\xi +f)为Kaiser优化窗的频谱.记：

p(\tau ,f)=x(\tau ){\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }f\tau }（12）

则变换式（6）的卷积形式可表达为：

K(\tau ,f)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}p(t,f)w_{\mathrm{k}}(\tau -t,f)\mathrm{d}t=

p(\tau ,f)\cdot w_k(\tau ,f)（13）

式中：*表示卷积.

设P(\xi ,f)为p(\tau ,f)关于时间\tau的FFT，则对p(\tau ,f)和w_{\mathrm{k}}(\tau ,f)作关于时间\tau的FFT，有：

P(\xi ,f)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}p(\tau ,f){\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }\xi \tau }\mathrm{d}\tau =

{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(\tau ){\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }(\xi +f)\tau }\mathrm{d}\tau =X(\xi +f)（14）

Kaiser优化窗频谱为：

W_k(\xi ,f)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{w}_k(\tau ,f)e^{-\mathrm{j}\mathrm{2}\mathrm{\pi }\xi \tau }\mathrm{d}\tau（15）

对式（13）进行FFT，由卷积定理可得：

G(\xi ,f)=X(\xi +f)W_k(\xi ,f)（16）

式中：G(\xi ,f)为K(\tau ,f)关于时间\tau的FFT.

对式（16）进行IFFT（FFT逆变换）可得到K-S变换的FFT表达即式（10）.

可见，K-S变换可由其信号FFT与Kaiser优化窗的频率谱实现.具体如下：1）将信号x（t）的FFT频谱向量进行频移，再与该频移点下Kaiser优化窗的频谱向量进行Hadamard乘积；2）将频谱乘积结果进行IFFT即实现该频率点下K-S变换，构造出K-S变换复时频矩阵；3）提取K-S变换对信号x（t）的时间-频率-幅值、时间-频率-相位三维信息.具体地，对 式（10）进行离散化处理，由FFT和IFFT对算法进行快速实现，变换结果的复时频矩阵\mathit{K}[n,m]可由幅值矩阵\mathit{A}[n,m]和相位矩阵\mathit{\varphi }[n,m]表示为：

\mathit{K}[n,m]=\mathit{A}[n,m]{\mathrm{e}}^{\mathrm{j}\varphi [n,m]}（17）

其中，

\mathit{A}[n,m]=\sqrt[]{\mathrm{R}\mathrm{e}{\left[\mathit{K}[n,m]\right]}^{\mathrm{2}}+\mathrm{I}\mathrm{m}{\left[\mathit{K}[n,m]\right]}^{\mathrm{2}}}（18）

\mathit{\varphi }[n,m]=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{\mathrm{R}\mathrm{e}\left[\mathit{K}[n,m]\right]}{\mathrm{I}\mathrm{m}\left[\mathit{K}[n,m]\right]}（19）

幅值矩阵\mathit{A}[n,m]和相位矩阵\mathit{\varphi }[n,m]的列向量分别表示信号在某一频率处的幅值和相位随时间变化的分布信息，行向量则分别表示信号某一采样时刻的幅值和相位随频率变化的分布.因此，通过K-S变换不仅可以获得某一时刻的频率信息，还可获得在某一频率上的信号幅值信息.

2.3 逆变换

逆变换才能保证信号从一维时域到二维时频域变换不存在能量损失，确保变换算法的科学性.STFT、ST能够推导逆变换的根本条件是时频变换算法中窗函数的时域积分值为常数.而这种窗函数面积限制正是制约STFT、ST变分辨率特性的关键.突破时频分析领域非限定面积窗函数的应用难关，可为具有卓越性能窗函数的更广泛应用提供依据.

将式（6）在时域对Kaiser优化窗时移因子进行积分，即将时频域变换结果在时域整合，可得到信号频域信息.通过时移因子积分变换，确定相关的频率点参数系数，继而通过IFFT可得到原始信号.

信号x(t)的K-S逆变换表达式为：

x(t)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left\{\frac{\mathrm{1}}{{\beta }_f}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}K(\tau ,f)\mathrm{d}\tau \right\}{\mathrm{e}}^{\mathrm{j}\mathrm{2}\mathrm{\pi }ft}\mathrm{d}f（20）

式中：{\beta }_f为Kaiser优化窗函数在频率点f的系数参数，即当Kaiser优化窗调节因子的值为α（f）时，Kaiser优化窗的面积，具体表达如下：

{\beta }_f={\int }_{-T}^T\frac{I_{\mathrm{0}}\left[\alpha \left(f\right)\sqrt[]{\mathrm{1}-{\left(v/\mathrm{T}\right)}^{\mathrm{2}}}\right]}{I_{\mathrm{0}}\left[\alpha \left(f\right)\right]}\mathrm{d}v,\begin{array}{cc}& \left|v\right|\le \mathrm{T}\end{array}（21）

事实上，

{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}K(\tau ,f)\mathrm{d}\tau ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left\{{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(t)w_k(\tau -t,f){\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }ft}\mathrm{d}t\right\}\mathrm{d}\tau =

{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(t){\mathrm{e}}^{-i\mathrm{2}\mathrm{\pi }ft}{\int }_{t-T}^{t+T}\frac{I_{\mathrm{0}}\left[\alpha (f)\sqrt[]{\mathrm{1}-[{(\tau -t)/T]}^{\mathrm{2}}}\right]}{I_{\mathrm{0}}[\alpha (f)]}\mathrm{d}\tau \mathrm{d}t（22）

令v=\tau -t，\left|v\right|\le T，则式（22）可改写成：

\begin{array}{l}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}K(\tau ,f)\mathrm{d}\tau =\\ {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x(t)e^{-\mathrm{j}\mathrm{2}\mathrm{\pi }ft}{\int }_{-T}^T\frac{I_{\mathrm{0}}\left[\alpha (f)\sqrt[]{\mathrm{1}-{(v/\mathrm{T})}^{\mathrm{2}}}\right]}{I_{\mathrm{0}}[\alpha (f)]}\mathrm{d}v\mathrm{d}t=\end{array}

{\beta }_{f}X(f)（23）

利用IFFT与式（23）可得：

{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left\{\frac{\mathrm{1}}{\beta }{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}K\left(\tau ,f\right)\mathrm{d}\tau \right\}{\mathrm{e}}^{\mathrm{j}\mathrm{2}\mathrm{\pi }ft}\mathrm{d}f={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}X\left(f\right){\mathrm{e}}^{\mathrm{j}\mathrm{2}\mathrm{\pi }ft}\mathrm{d}f=

x\left(t\right)（24）

即式（20）成立，K-S逆变换可利用IFFT得出.

此外，由于信号K-S变换可由其FFT，结合Kaiser优化窗频率谱而快速实现.因此，亦可利用信号K-S变换的FFT表达来对信号进行快速重构.具体地，1）对信号K-S变换的FFT表达式（10）进行FFT运算，其结果为原始信号x(t)的FFT频谱进行频移后在该频率点的频谱X(\xi +f)与Kaiser优化窗频谱W_k(\xi ,f)的向量Hadamard积，即式（16）；2）将其与窗频谱向量元素的倒数向量作Hadamard积，得到原始信号FFT频谱频移向量\mathit{X}(\xi +f)；3）通过频移与IFFT运算可对原始信号进行重构.

2.4 K-S变换算法流程

K-S变换流程如图2所示.算法实现步骤如下：

1）对被测信号x（t）采样并进行自适应抽样，计算抽样数据 x ［n］的FFT，得到 X ［k］；

2）阈值滤波、非特征信息剔除，得到特征信息点 X ［m］；

3）通过FFT获得频率点m处Kaiser优化窗的频率谱W［m］；

4）对于频率点m，计算步骤2）中X［m］与步骤3）中 W ［m］的向量Hadamard积；

5）对步骤4）中获得的乘积向量进行IFFT运算，得到频率点m的K-S时频变换；

6）重复步骤3）~5），得到信号所有频率点对应的复时频矩阵.最终，通过式（18）和式（19）提取信号的时频幅值信息和时频相位信息.

3 K-S变换的基本性质

K-S变换为复数变换，因此K-S矩阵为复数矩阵，K-S变换的结果可表示为式（17）的形式.式（18）展示的幅值矩阵称为K-S模矩阵.

3.1 局部性质

参见式（23），信号x（t）的K-S变换K(\tau ,f)沿时间轴积分，等效于x（t）的傅里叶变换X(f)与Kaiser优化窗频率系数的乘积，此即信号x(t)的时间局部频谱.另一方面，K(\tau ,f)可以通过式（20）恢复到信号x（t）.可见，对于非平稳信号，K-S变换可看作一种广义傅里叶变换.K-S变换从时域变换到时频域，再变换到频域，然后回到时域，整个变换过程是无损、可逆的，而傅里叶变换是K-S变换的运算核心.

信号x（t）经K-S变换，得到的K-S二维矩阵中的行为时间、列为频率，且矩阵各行所对应的为该时间点的局部频谱.因此K-S变换结果在时频平面上有直观、易于理解的展示.

3.2 线性性质

若有x(t)\in L^{\mathrm{2}}(R)、y(t)\in L^{\mathrm{2}}(R)两信号，k_{\mathrm{1}}、k_{\mathrm{2}}为任意常数，设信号x（t）、y(t)的K-S变换分别为K_x （τ，f）、K_y （τ，f），并设x（t）、y（t）的线性运算为：

z(t)=k_{\mathrm{1}}x(t)+k_{\mathrm{2}}y(t)（25）

则z(t)的K-S变换为：

K_z(\tau ,f)=k_{\mathrm{1}}{K}_x(\tau ,f)+k_{\mathrm{2}}{K}_y(\tau ,f)（26）

可见，K-S变换对于时间信号是线性运算.

若x（t）由原始信号n（t）与加性噪声h（f）构成，即：

x(t)=h(t)+n(t)（27）

则x(t)的K-S变换为：

K_x(\tau ,f)=K_h(\tau ,f)+K_n(\tau ,f)（28）

式中：K_h （τ，f）、K_n （τ，f）分别为h（f）、n（t）的K-S变换结果.可见，变换结果中不存在交叉项，不受交叉项影响，K-S变换不仅有效提高了信号联合时频分析的分辨率，而且能有效抑制信号的加性噪声影响，这一性质大力增强了K-S变换的应用价值.

3.3 变分辨率特性

引入调节因子α（f）后，可根据信号频率自适应调整Kaiser优化窗宽度，参见式（4），频率f通过调节因子α（f）作用于Kaiser优化窗，为获得优良的调控效果，检测频率f较高时，c₁取值需较大，而c₂则用于对窗函数的性能进行微调.

频率调节因子α（f）随检测频率f增大而逐渐减小.当f较小时，Kaiser优化窗的窗形相应较窄，旁瓣变低，主瓣较宽，有较好的时间分辨率.随着检测频率f增大，优化窗的窗形变宽，主瓣变窄，逐渐往频率分辨率更优的方向变化.当α（f）=0时，Kaiser优化窗退化成矩形窗，主瓣变窄，可提供更优的频率分辨率.特别地，K-S变换中时频分辨率变化速率可由窗控制参数c₁决定：较大c₁对应着较快的转变速率，c₂则可针对某一特定频率下检测需求有针对性地对时频分辨率进行相应调整.对0~150 kHz频率范围复杂信号的仿真实验表明，随着信号频率增加，Kaiser优化窗的宽度变宽，如图3所示.由于窗函数灵活的变分辨率特性，K-S变换既可获得信号某一时刻准确的频率信息，又可获得信号在某一频率上的准确幅值信息，具优良变分辨率特性.

4 基于K-S变换的电网超谐波检测

为证实K-S变换算法在较宽频率范围具有可靠的变分辨率特性和时频分辨能力，利用电网超谐波信号^［18］进行仿真与测试实验.依据IEEE 159： 2019^［19］，合成的仿真实验信号包括稳态超谐波信号和时变超谐波信号，采样频率为409.6 kHz，信号分解方法包括经典STFT、ST和本文建立的K-S变换算法.STFT窗函数采用标准差为35的Gaussian窗.Kaiser优化窗函数控制参数与图1（a）中参数相同.

4.1 超谐波检测仿真实验研究

4.1.1 稳态超谐波信号仿真实验

对于稳态超谐波信号：

\begin{array}{l}{x}_{\mathrm{1}}(t)=\mathrm{0.8}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{10}\times {\mathrm{10}}^{\mathrm{3}}\times \mathrm{2}\mathrm{\pi }t\right)+\\ \mathrm{0.6}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{80}\times {\mathrm{10}}^{\mathrm{3}}\times \mathrm{2}\mathrm{\pi }t\right)+\end{array}

\mathrm{0.3}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{150}\times {\mathrm{10}}^{\mathrm{3}}\times \mathrm{2}\mathrm{\pi }t\right)（29）

三种算法得到的三维时频谱、能量聚集情况分别如图4~图5所示.可见，K-S变换相对STFT、ST有一下特征：1）对稳态宽频范围超谐波有更突出的频率检测能力，能准确得到各次超谐波的幅值信息；2）K-S变换的频率检测能力更佳，能精确检测出信号中存在的10 kHz、80 kHz和150 kHz超谐波.

4.1.2 时变超谐波信号仿真实验

对于时变超谐波信号：

x_{\mathrm{2}}(t)=\mathrm{0.8}\mathrm{s}\mathrm{i}\mathrm{n}\left(b_{\mathrm{1}}\left(t\right)\times {\mathrm{10}}^{\mathrm{3}}\times \mathrm{2}\mathrm{\pi }t\right)+

\mathrm{0.8}\mathrm{s}\mathrm{i}\mathrm{n}\left(b_{\mathrm{2}}\left(t\right)\times {\mathrm{10}}^{\mathrm{3}}\times \mathrm{2}\mathrm{\pi }t\right)（30）

式中：

b_{\mathrm{1}}\left(t\right)=\left\{\begin{array}{c}\mathrm{15}\begin{array}{cc}& t\in \left[\mathrm{0,0.5}T\right]\end{array}\\ \mathrm{10}\begin{array}{cc}& t\in \left[\mathrm{0.5}T,T\right]\end{array}\end{array}\right.（31）

b_{\mathrm{2}}\left(t\right)=\left\{\begin{array}{c}\mathrm{100}\begin{array}{cc}& t\in \left.\left[\mathrm{0},\right.\mathrm{0.25}T\right)\bigcup \left[\mathrm{0.75}T,T\right]\end{array}\\ \mathrm{80}\begin{array}{cc}\begin{array}{ccc}& t\in \left[\mathrm{0.25}T,\mathrm{0.75}T\right]& \begin{array}{cc}& \end{array}\end{array}& \end{array}\end{array}\right.（32）

三种算法得到的时变超谐波三维时频谱、能量聚集情况分别如图6和图7所示.可见：K-S变换对时变超谐波有最突出的频率检测能力和最少的频谱泄漏，能准确获得时变超谐波的时间和幅值信息.

4.2 超谐波信号测试实验

为验证K-S变换的准确性，构建了图8所示的测试实验平台.其中信号源采用Agilent（安捷伦）的 33 500B波形发生器，信号经调理电路后送16位A/D转换器ADS8568，A/D转换器采样频率设置为400 kHz.采集得到的离散信号送TMS320C6748中进行数据预处理和K-S变换快速算法处理，LCD显示器显示时频参数分析结果.也可将采样数据直接交送上位机进行K-S变换处理和显示.

为检验算法的实际性能与工程应用可行性，大力简化计算量，形成了高效且准确的K-S变换快速算法：1）依据信号频率与频段分析需求，建立了信号抽样机制，有效减少低频信号和信号低频段分析时的数据长度，减少了信号运算量；2）探究变换过程中得到的二维时频复数矩阵所具有的关联特性，通过阈值滤波和感兴趣频率确定等，建立了特征频率保留机制及非特征信息剔除机制.具体简化过程不赘述.受信号源频率限制，仅利用实验平台分别对 1 kHz、5 kHz、10 kHz、20 kHz、40 kHz五个频率分量的谐波和超谐波信号进行5次测量，得到的各频率点幅值测量数据如表1所示.由表1可见，测量最低绝对误差为0.000 1 V，所有测量结果中绝对误差值均优于0.032 3V.

5 结论

本文突破了线性时频变换对窗函数面积要求为1的桎梏，建立了窗函数无穷积分面积无约束要求的K-S变换，给出K-S变换及其逆变换的数学推导，介绍了变换算法的局部性质、线性性质和变分辨率特性，从理论上确保了变换算法的科学性，获得了宽频复杂信号时频分析的高性能能量聚集性和频率自适应调节能力.0~150 kHz电网超谐波的仿真实验与0~40 kHz电网谐波、超谐波的测试实验证明了K-S变换的可行性与准确性.K-S变换为宽频复杂动态信号的时频联合分析提供了一种新的高效且准确的方法.

参考文献

[1] 张玮， 王平 ．多跳频信号参数估计时频分析算法［J］．探测与控制学报， 2023， 45（1）：113-118．

[2] ZHANG W， WANG P ．Multi-frequency hopping signal parameter estimation method of time-frequency analysis［J］．Journal of Detection & Control，2023，45（1）：113-118．（in Chinese）

[3] 郝国成，谈帆，程卓，等 ．强鲁棒性和高锐化聚集度的BGabor-NSPWVD时频分析算法［J］．自动化学报，2019，45（3）：566-576．

[4] HAO G C， TAN F， CHENG Z，et al ．Time-frequency analysis of BGabor-NSPWVD algorithm with strong robustness and high sharpening concentration［J］．Acta Automatica Sinica，2019， 45（3）： 566-576．（in Chinese）

[5] 刘会杰， 高新海， 郭汝江 ．一种低副瓣无混叠的线性调频信号时频分析方法［J］．电子与信息学报， 2019， 41（11）： 2614-2622．

[6] LIU H J， GAO X H， GUO R J ．A time-frequency analysis method for linear frequency modulation signal with low sidelobe and nonaliasing property［J］．Journal of Electronics & Information Technology，2019，41（11）：2614-2622．（in Chinese）

[7] KORN P ．Some uncertainty principles for time-frequency transforms of the Cohen class［J］．IEEE Transactions on Signal Processing， 2005， 53（2）：523-527．

[8] 黄昱丞， 郑晓东， 栾奕， 等 ．地震信号线性与非线性时频分析方法对比［J］．石油地球物理勘探，2018， 53（5）： 975-989．

[9] HUANG Y C， ZHENG X D， LUAN Y，et al ．Comparison of linear and nonlinear time-frequency analysis on seismic signals［J］．Oil Geophysical Prospecting，2018，53（5）：975-989．（in Chinese）

[10] GU Y H， BOLLEN M H J ．Time-frequency and time-scale domain analysis of voltage disturbances［J］．IEEE Transactions on Power Delivery，2000，15（4）：1279-1284．

[11] YOON W K， DEVANEY M J ．Reactive power measurement using the wavelet transform［J］．IEEE Transactions on Instrumentation and Measurement，2000，49（2）：246-252．

[12] 张淑清，李盼，张立国，等 ．广义S变换的参数优化及在电能质量分析中的应用［J］．中国科学：技术科学，2016，46（6）：593-601．

[13] ZHANG S Q， LI P， ZHANG L G，et al ．Parameter optimized method for generalized S-transform algorithm and its application in the analysis of power quality disturbances［J］．Scientia Sinica （Technologica），2016，46（6）：593-601．（in Chinese）

[14] 李思源，徐天吉 ．基于Wigner-Ville分布与Chrip-Z变换的高分辨时频分析方法［J］．石油地球物理勘探，2022，57（1）：168-175．

[15] LI S Y， XU T J ．A new high-resolution time-frequency analysis method based on Wigner-Ville distribution and Chrip-Z transform［J］．Oil Geophysical Prospecting，2022，57（1）：168-175．（in Chinese）

[16] DJUKANOVIĆ S ．An accurate method for frequency estimation of a real sinusoid［J］．IEEE Signal Processing Letters，2016， 23（7）： 915-918．

[17] HUANG N E， SHEN Z， LONG S R，et al ．The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］．Proceedings of the Royal Society of London Series A，1998，454（1971）：903-998．

[18] 蒋海峰，张曼，赵斌炎，等 ．基于改进Hilbert-Huang变换的电网故障诊断［J］．电工技术学报，2019，34（S1）：336-342．

[19] JIANG H F， ZHANG M， ZHAO B Y，et al ．Power grid fault diagnosis based on improved Hilbert-Huang transform［J］．Transactions of China Electrotechnical Society，2019，34（S1）：336-342．（in Chinese）

[20] YAO W X， TENG Z S， TANG Q，et al ．Measurement of power system harmonic based on adaptive Kaiser self-convolution window［J］．IET Generation，Transmission & Distribution，2016，10（2）：390-398．

[21] STOCKWELL R G， MANSINHA L， LOWE R P ．Localization of the complex spectrum：the S transform［J］．IEEE Transactions on Signal Processing，1996，44（4）：998-1001．

[22] MOUKADEM A， OULD ABDESLAM D， DIETERLEN A ．Time–Frequency Domain for Segmentation and Classification of Non-Stationary Signals［M］.2014．

[23] HE S F， LI K C， ZHANG M. A real-time power quality disturbances classification using hybrid method based on S-transform and dynamics［J］. IEEE Transactions on Instrumentation and Measurement，2013， 62（9）： 2465-2475．

[24] LI J M， TENG Z S， TANG Q，et al ．Detection and classification of power quality disturbances using double resolution S-transform and DAG-SVMs［J］. IEEE Transactions on Instrumentation and Measurement， 2016， 65（10）： 2302-2312．

[25] RONNBERG S， BOLLEN M， AMARIS H，et al ．On waveform distortion in the frequency range of 2 kHz–150 kHz—review and research challenges［J］．Electric Power Systems Research，2017，150：1-10．

[26] IEEE Std 1159.2—2009.IEEE recommended practice for monitoring electric power quality［C］. 2009.

基金资助

国家自然科学基金资助项目（52077067,51377049）