基于LS-TSVM与降维ESPRIT谐波检测的研究

柏 林，唐 智，刘小峰，刘子军

(1.重庆大学 机械传动国家重点实验室，重庆 400044；2.重庆聚星仪器有限公司，重庆 400044)

高速列车供电系统的电流与电压波形产生畸变后，对输电、供电系统、用电负载产生一定的影响，其电力谐波[1]可能会造成自控系统的二次保护误动作，影响列车正常运行，造成晚点，影响公众出行。因此，对电网中谐波进行快速、准确地检测对于维护电力系统安全、经济运行都具有重要意义。

目前，针对电网谐波的检测方法有傅里叶变换、小波变换以及神经网络等，但各种方法均存在一定的局限性。如基于小波变换的谐波检测方法[2]，由于不同尺度的小波函数在频域中存在互相干扰，对频率相近的谐波分离困难；基于神经网络的谐波检测方法[3]需要大量的训练样本，且在电网频率波动情况下的检测结果存在较大误差；参数化谐波检测方法，如拓展Prony方法[4]具有较高的频率分辨率，但模型参数计算较复杂，计算量大，且抗噪能力较差；基于FFT的谐波检测方法[5]由于电网基波频率变化，导致在非同步采样或非整周期截断的情况下产生较严重的频谱泄漏，使测量精度降低。

快速发展的电气化铁路，相对于其他谐波源，其谐波产生量和谐波影响范围都更加突出，且以高次谐波居多，再加上其快速移动的用电特性，其谐波检测有其特殊性，在保证检测效率的前提下，应具有较高的谐波分辨能力。本文针对传统谐波检测方法的缺陷和牵引供电系统的特点，提出基于双子支持向量机结合多级维纳滤波器ESPRIT的牵引供电系统谐波检测方法。仿真实验表明：该方法在确保分辨率、精度前提下，具有较好的鲁棒性和较少的计算量。采用该方法对某型试验列车各主要环节进行谐波检测，取得了较为理想的测试结果。该研究方法有利于未来电气化铁路谐波检测标准的顺利实施。

1 最小二乘双子支持向量机

双支持向量机(TSVM)把传统SVM[6]中运算量较大的规划问题一分为二，将原来的一个求解问题转化为一对问题进而求最优。采用了超松弛下降法进行求解，该方法相较于最小二乘算法来说，具有运算量大、耗时长等缺点，且最小二乘算法已能运用在SVM以及TSVM中，故提出结合最小二乘算法的TSVM[7]。

对线性样本集(xi，yi)进行回归预测问题，TSVM的思路是寻求一对线性函数。

( 1 )

式中：w1，w2∈Rn；b1，b2∈R。通过对f1(x)与f2(x)进行平移，可得控制区域的上界限与下界限，在尽量压缩上下界限的同时，保证输入的样本(xi，yi)能在上下界之间。利用该样本集训练回归函数。

( 2 )

使f(x)预测值与实际观测值之间的误差ε允许在误差范围内较好地拟合样本集。最小二乘的核心是把SVM中以不等式形式存在的约束条件转化为用等式进行约束，并且惩罚项用一次项替换为二次项，修改了不等式的约束条件，用一个等式约束来替代，最小二乘TSVM求解两个优化问题。

s.t.Y-(Xw1+eb1)=-ε1e-ξ1

( 3 )

s.t. (Xw2+eb2)-Y=-ε2e-ξ2

( 4 )

式中：c1，c2，c3，c4，ε1，ε2是正的参数。对于式( 3 )的优化，目标函数中除去正则项为‖Y-(Xw1+eb1)‖2，极小化中间项意味着要求函数f1(x)尽可能拟合训练样本，后面的松弛项结构限制条件要求训练样本的输出尽可能大于控制下线f1(x)-ε1，同样对于优化式( 4 )，目标函数除正则项外要求f2(x)尽可能小于控制上限f2(x)+ε2。为使参数ε1、ε2达到最小值，从而使训练样本被“夹”在这对控制上下界函数之间。分别求对w1与b1的偏导，令其为零，即

( 5 )

( 6 )

获得正确的权值和偏置后，便求得回归预测函数f(x)。从上述推导可见，LS-TSVM计算过程中并未出现参数迭代求解，故对应的运算量较小。

2 LS-ESPRIT频率提取方法

ESPRIT原理是正弦信号的频率以及幅值不会因为时移而发生变化，故可用该方法估计其频率成分，进而估计出待测信号幅值及初始相位。

假设s(n)是由K个正弦分量加上观测噪声组成

( 7 )

式中：K为谐波成分的数目；ωk、ak、φk对应谐波k的频率、幅值以及相位；u(n)表示信号中的噪声成分。忽视直流部分，可将式( 7 )改写成指数形式

( 8 )

复正弦部分可表示为

( 9 )

使用信号子空间表示方法可将式( 8 )写成M维等式

S(n)=EΦnA+U(n)

(10)

其中，

S(n)=[s(n)s(n+1) …s(n+M-1)]T

U(n)=[u(n)u(n+1) …u(n+(M-1))]T

Φ=diag{ejω1，ejω2…，ejωK}

由式( 9 )可得

(11)

同理可得

S(n)=S(n+1)=

[s(n+1)s(n+2) …s(n+M)]T

(12)

S′(n)=EΦn+1A+U(n+1)

(13)

由自相关、互相关矩阵定理可得

(14)

(15)

用求信号子空间特征方程的问题求解谐波频率

Rtαi=λiRttαii=1，…，K

(16)

将式(14)、式(15)带入式(16)得

EAAHEHαi=λiEAAHEαii=1，…，K

(17)

即

EAAH(I-λiΦH)EHαi=0i=1，…，K

(18)

由式(18)求得谐波频率为

λi=ejωi

ωi=angle(λi)i=1，…，K

(19)

3 基于降维ESPRIT-TSVM谐波检测方法

尽管能够用ESPRIT进行频率计算，但该过程需计算样本协方差矩阵的逆以及特征值分解，同样需要较大的运算量。高次谐波广泛存在于牵引电力系统中，为准确检测出谐波成分并且提升算法时效，有必要在得到谐波成分的频率信息之后再进行TSVM分析。

本文为降低ESPRIT计算量，采用基于多级维纳滤波器(MSWF)[8-9]对信号进行前置降维处理。图1为MSWF结构。

图1 多级维纳滤波器结构

(20)

逐级反推得到标量维纳权序列为

(21)

总的最佳权为

WMWF=TDWD

(22)

MSWF的前向部分与后向部分均含有迭代，若用D表示迭代数，降维则是将MN-1维的信号映射为D维度矢量，即对滤波的计算量进行了降低。

采用TSVM结合ESPRIT的谐波检测算法整体步骤如下：

步骤1用采集得到的时域信号构成阵列接收阵

S(n)=[s(n)s(n+1) …x(n+M-1)]T

式中：M为阵元数；M+n为总采样数目。

步骤2当连续3次分解以后，期望数据di与观测数据Xi的互相关系数小于阈值时，则降维结束。根据式(20)、式(21)计算降维矩阵TD和综合矩阵WD，利用MSWF降维技术对空域信号接收矩阵进行降维，观测输入信号经降维矩阵TD处理之后得到的D维向量XD=TDWDS(n)。

步骤3利用前降维后的信号子空间XD取代式(10)中的被分析时空域信号S(n)，根据式(13)～式(19)得谐波频率ωi=angle(λi)(i=1，…，K)。

步骤4将步骤3中的谐波频率ωi组成输入向量

z(n)=[cos(ω1n) … cos(ωMn)

sin(ω1n) … sin(ωMn)]T

令w=[a11…a1Ma21…a2M]T，忽略噪声项，则式( 7 )中的被分析信号可表示为

(23)

式中：a1k=Aksin(φk)；a2k=Akcos(φk)。根据参数a1k与a2k可得谐波的幅值Ak与相位φk。

(24)

求解a1k与a2k可变换为求式(23)的线性回归问题，通过TSVM求解回归系数w，即可得到每个频率成分对应的Ak与φk。

步骤5采用阈值去除方式，把幅值小于阈值的成分筛掉，从而减少虚假频率的干扰。

4 综合性能仿真法分析

这里从算法的鲁棒性、精度、时效3个维度对其进行测试分析，并与其他算法进行对比，用于算法测试的计算机硬件参数：CPU i5-2430，4GB运存，仿真软件为

MATLAB 2015b。仿真信号为

f(x)=100sin(2π·49.8t)+5sin(2π·1 444.2t+

π/2)+sin(2π·1 494t+π/3)+5sin(2π·1 543.8t+

π/4)+sin(2π·2 490t)

信真信号时域波形如图2所示。采样频率为16 kHz，采样时宽30 ms，共采样480点。设c1=c2=0.000 1，ε1=ε2=0.01，阵元数200，快拍数280，信号加高斯白噪声。

图2 仿真信号时域波形

本文方法的主要特点是在ESPRIT进行频率估计过程中加入多级维纳滤波器进行降维，减少运算量。将采用维纳滤波器降维的ESPRIT频率估计方法与常规ESPRIT频率估计进行对比，用MATLAB编写程序，计算结果和计算时长见表1。

表1 降维前后ESPRIT及FFT比较

表1中计算时长采用20次实验的平均值作为最终结果。从表1数据可以看出，采用维纳滤波器的ESPRIT频率估计方法和常规ESPRIT频率估计方法在精度上比较接近，但在计算时长上，采用维纳滤波器降维的ESPRIT频率估计大约需要0.007 219 s，而不采用维纳滤波器的ESPRIT频率估计方法大约需要0.103 042 s，采用维纳滤波器降维的ESPRIT能有效减少频率估计等待时长，提高计算效率。

除了运算效率，表1中我们还将该方法与传统的FFT频率计算结果进行了比较，从表1可以看出FFT在计算时长上与降维的ESPRIT方法相差不大，但是在相同的采样数据长度(480点)的情况下，FFT的计算精度远不及ESRIT。

本文还研究不同噪声环境对降维ESPRIT频率估计精度的影响，对采集到的信号分别加入30、50、70 dB的高斯白噪声进行计算，其结果如图3所示。

图3 噪声对计算精度的影响

根据计算结果可以看出，在添加30、50、70 dB的高斯白噪声后，运算结果差距并不大，故本文的频率估计算法除了运算效率高以外，同时也具有较好的鲁棒性。

用表1得到的频率信息构造最小二乘双子支持向量机训练的输入向量，继而求出谐波的幅值和相位。本文还将TSVM算法和SVM算法进行了对比，对比结果见表2。

表2 TSVM与标准SVM对比

由表2数据可知，TSVM与SVM在计算谐波幅值和相位时计算精度差距不大，但TSVM算法在计算时长上占有绝对的优势。TSVM计算时长大约为SVM计算时长的四分之一，这也和理论研究相吻合。TSVM把规划问题一分为二后，规划问题的复杂度仅为原来的1/2。假设样本数量为n，SVM的运算量为8n3，则TSVM的运算量为2n3。

运用本文方法完整检测出谐波的频率、幅值、相位后，还与矩阵束法[12]进行了比较。采用矩阵束算法计算所得结果见表3。

表3 矩阵束检测结果

采用上述两方法对20组数据统计绝对误差的均方根值，以降低算法对比的随机性，其结果见表4。

表4 绝对误差的均方根值

从表4的结果分析可知：当信噪比为30 dB时，两种方法的精度较好，但本文算法略高于矩阵束算法。文中方法在降维过程中若MSWF分解的级数多于待测信号中的谐波数目，可能引入虚假频率，在20次实验中，虚假频率并不会持续存在，其幅值常低于0.01。20次连续仿真中，本文算法最高运行时长为0.009 358 s，最少运行时长为0.005 597 s，运行时长的均方根值为0.007 6 s；而矩阵束法对应时长为0.028 0、0.014 5、0.018 7 s。相比而言，矩阵束法明显比本文算法的时效性差，当待检信号谐波频率较高时，根据采样定理需要提高采样频率，此时数据点数增多，由于本文算法进行了降维，应用价值更高。

5 应用

研究列车供电系统谐波监测方法不但能时刻监测列车电力系统的状态，还能用于设计者检验列车谐波特性是否满足设计要求。图4为试验列车主电路图。该方案采用的硬件设备有德国CRONOS-PL以及美国NI公司的数据采集卡，下位机由imc devices对电力数据采集通道进行配置，上位机采用具有数据显示等相关功能的LabView程序，并对记录的数据采用MATLAB离线分析。

图5是列车运行速度为300 km/h时采集的牵引电机电流信号，从图5时域波形可以看出信号中有较为严重的畸变，从采集的信号中截取200 ms进行分析，加汉宁窗，得其频谱图，如图6(a)所示，截取40 ms数据应用本文算法进行分析处理，计算谐波参数，得到频谱图，如图6(b)所示。

为了进一步验证ESPRIT-TSVM方法的实用性，采用该方法对动车电力系统进行谐波分析，同时列出了FFT算法的分析结果。

从图6可以看出，采用FFT得到的谐波成分与采用本文方法提取的频谱图数值较为接近，在一定程度上也佐证了本文方法的准确性，实测发现采用本文方法计算得到的谐波成分予以剔除，优于直接采用FFT法。

图4 试验列车主电路图

图5 列车速度300 km/h电流信号

(a)傅里叶频谱

(b)本文提取谱图6 列车频谱图

图7(a)为本方案上位机测试软件的主界面，主要信息包含各个测点的谐波占比、谐波谱图等，采用图表加文字的方式更加清晰；图7(b)为列车时速380 km时的网侧电流分析结果，主要包含33次以下的奇次谐波，谐波成分中23、25、27、29、33次谐波所占比重较大，若不加以控制，可能会导致较大干扰电流的产生。

(a)测试软件主界面

(b)380 km/h网侧电流谐波分布图7 测试软件主界面及网侧电流谐波分布

从图8可以看出，在列车运行速度小于210 km/h时，随着车速的增加，THD值迅速降低，但仍比要求的10%要高；在列车速度超过210 km/h后THD值变化比较平稳，表明此时THD受车速的影响开始变小，但普遍大于8%；在车速大于450 km/h后，THD又呈现出增加的趋势。

图8 不同速度下的网侧电流谐波畸变率

6 结论

本文在支持向量机结合ESPRIT谐波检测的方法上进行改进，提出LS-TSVM结合多级维纳滤波器降维ESPRIT的检测方法。利用多级维纳滤波器降维的ESPRIT快速估计出谐波的频率，用得到的频率信息构造最小二乘双子支持向量机的输入向量得到谐波的幅值与相位。经过仿真试验获得如下结论：

(1)在相同的电网谐波数据量之下，该算法与传统的FFT等方法相比具有更好的计算精度。

(2)通过采用多级维纳滤波器降维后，降低了程序的运行时长，提高了谐波检测的效率。

(3)本文方法在不同的噪声环境均保持可观的精度，具有较好的鲁棒性。

该方法能准确检测出电网谐波的频率、幅值以及相位，具有较好的实用价值，并且已在某型列车电网中进行应用，反映良好。

参考文献：

[1]IEC.IEC 61000-4-7.General Guide on Harmonics and Interharmonics Measurement and Instrumentation,for Power Supply Systems and Equipment Connected Thereto[S].2002.

[2]EREN L,UNAL M,DEVANEY M J.Harmonic Analysis via Wavelet Packet Decomposition Using Special Elliptic Half-band Filters[J].IEEE Transactions on Instrumentation & Measurement,2007,56(6):2289-2293.

[3]LIN H C.Intelligent Neural Network-based Fast Power System Harmonic Detection[J].IEEE Transaction on Industrial Electronics,2007,54(1):43-52.

[4]COSTA F F,CARDOSO A J M.Harmonic and Interharmonic Identification Based on Improved Prony’s Method[C]//32nd Annual Conference on IEEE Industrial Electronics.New York:IEEE,2006:1047-1052.

[5]牛胜锁,梁志瑞,张建华,等.基于四项余弦窗三谱线插值FFT的谐波检测方法[J].仪器仪表学报,2012,33(9):2002-2008.

NIU Shengsuo,LIANG Zhirui,ZHANG Jianhua,et al.Harmonic Detection Approach Based on 4-term Cosine Window Triple-spectrum-line Interpolation[J].Chinese Journal of Scientific Instrument,2012,33(9):2002-2008.

[6]占勇,丁屹峰,程浩忠,等.电力系统谐波分析的稳健支持向量机方法研究[J].中国电机工程学报,2004,24(12):43-47.

ZHAN Yong,DING Yifeng,CHENG Haozhong,et al.Robust Support Vector Machine Method for Harmonic Analysis of Power System[J].Proceedings of the CSEE,2004,24(12):43-47.

[7]JAYADEV A,KHEMCHANDANI R,CHANDRA S.Twin Support Vector Machines for Pattern Classification[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,2007,29(5):905-910.

[8]丁前军,王永良,张永顺.一种多级维纳滤波器的快速实现——迭代相关相减法[J].通信学报,2005,26(12):1-7.

DING Qianjun,WANG Yongliang,ZHANG Yongshun.Fast Implementation of a Multistage Wiener—Iterative Correlation Subtraction[J].Journal of Communication,2005,26(12):1-7.

[9]黄磊,吴顺君,张林让,等.快速子空间分解方法及其维数的快速估计[J].电子学报,2005,33(6):977-981.

HUANG Lei,WU Shunjun,ZHANG Linrang,et al.Fast Subspace Decomposition Method and the Fast Estimation of Its Dimension[J].Acta Electronica Sinica,2005,33(6):977-981.

[10]黄磊.快速子空间估计方法研究及其在阵列信号处理中的应用[D].西安:西安电子科技大学,2005.

[11]WEISS A J,GAVISH M.Direction Finding Using ESPRIT with Interpolated Arrays[J].IEEE Transaction on Signal Processing,1991,39(6):1473-1478.

[12]朱瑞可,李兴源,赵睿,等.矩阵束算法在同步电机参数辨识中的应用[J].电力系统自动化,2012,36(6):52-55.

ZHU Ruike,LI Xingyuan,ZHAO Rui,et al.Application of Matrix Pencil Algorithm in Parameter Identification of Synchronous Machines[J].Automation of Electric Power Systems,2012,36(6):52-55.