基于LS-TSVM與降維ESPRIT諧波檢測的研究

柏 林，唐 智，劉小峰，劉子軍

(1.重慶大學(xué) 機械傳動國家重點實驗室，重慶 400044；2.重慶聚星儀器有限公司，重慶 400044)

高速列車供電系統(tǒng)的電流與電壓波形產(chǎn)生畸變后，對輸電、供電系統(tǒng)、用電負(fù)載產(chǎn)生一定的影響，其電力諧波[1]可能會造成自控系統(tǒng)的二次保護誤動作，影響列車正常運行，造成晚點，影響公眾出行。因此，對電網(wǎng)中諧波進(jìn)行快速、準(zhǔn)確地檢測對于維護電力系統(tǒng)安全、經(jīng)濟運行都具有重要意義。

目前，針對電網(wǎng)諧波的檢測方法有傅里葉變換、小波變換以及神經(jīng)網(wǎng)絡(luò)等，但各種方法均存在一定的局限性。如基于小波變換的諧波檢測方法[2]，由于不同尺度的小波函數(shù)在頻域中存在互相干擾，對頻率相近的諧波分離困難；基于神經(jīng)網(wǎng)絡(luò)的諧波檢測方法[3]需要大量的訓(xùn)練樣本，且在電網(wǎng)頻率波動情況下的檢測結(jié)果存在較大誤差；參數(shù)化諧波檢測方法，如拓展Prony方法[4]具有較高的頻率分辨率，但模型參數(shù)計算較復(fù)雜，計算量大，且抗噪能力較差；基于FFT的諧波檢測方法[5]由于電網(wǎng)基波頻率變化，導(dǎo)致在非同步采樣或非整周期截斷的情況下產(chǎn)生較嚴(yán)重的頻譜泄漏，使測量精度降低。

快速發(fā)展的電氣化鐵路，相對于其他諧波源，其諧波產(chǎn)生量和諧波影響范圍都更加突出，且以高次諧波居多，再加上其快速移動的用電特性，其諧波檢測有其特殊性，在保證檢測效率的前提下，應(yīng)具有較高的諧波分辨能力。本文針對傳統(tǒng)諧波檢測方法的缺陷和牽引供電系統(tǒng)的特點，提出基于雙子支持向量機結(jié)合多級維納濾波器ESPRIT的牽引供電系統(tǒng)諧波檢測方法。仿真實驗表明：該方法在確保分辨率、精度前提下，具有較好的魯棒性和較少的計算量。采用該方法對某型試驗列車各主要環(huán)節(jié)進(jìn)行諧波檢測，取得了較為理想的測試結(jié)果。該研究方法有利于未來電氣化鐵路諧波檢測標(biāo)準(zhǔn)的順利實施。

1 最小二乘雙子支持向量機

雙支持向量機(TSVM)把傳統(tǒng)SVM[6]中運算量較大的規(guī)劃問題一分為二，將原來的一個求解問題轉(zhuǎn)化為一對問題進(jìn)而求最優(yōu)。采用了超松弛下降法進(jìn)行求解，該方法相較于最小二乘算法來說，具有運算量大、耗時長等缺點，且最小二乘算法已能運用在SVM以及TSVM中，故提出結(jié)合最小二乘算法的TSVM[7]。

對線性樣本集(xi，yi)進(jìn)行回歸預(yù)測問題，TSVM的思路是尋求一對線性函數(shù)。

( 1 )

式中：w1，w2∈Rn；b1，b2∈R。通過對f1(x)與f2(x)進(jìn)行平移，可得控制區(qū)域的上界限與下界限，在盡量壓縮上下界限的同時，保證輸入的樣本(xi，yi)能在上下界之間。利用該樣本集訓(xùn)練回歸函數(shù)。

( 2 )

使f(x)預(yù)測值與實際觀測值之間的誤差ε允許在誤差范圍內(nèi)較好地擬合樣本集。最小二乘的核心是把SVM中以不等式形式存在的約束條件轉(zhuǎn)化為用等式進(jìn)行約束，并且懲罰項用一次項替換為二次項，修改了不等式的約束條件，用一個等式約束來替代，最小二乘TSVM求解兩個優(yōu)化問題。

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

( 3 )

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

( 4 )

式中：c1，c2，c3，c4，ε1，ε2是正的參數(shù)。對于式( 3 )的優(yōu)化，目標(biāo)函數(shù)中除去正則項為‖Y-(Xw1+eb1)‖2，極小化中間項意味著要求函數(shù)f1(x)盡可能擬合訓(xùn)練樣本，后面的松弛項結(jié)構(gòu)限制條件要求訓(xùn)練樣本的輸出盡可能大于控制下線f1(x)-ε1，同樣對于優(yōu)化式( 4 )，目標(biāo)函數(shù)除正則項外要求f2(x)盡可能小于控制上限f2(x)+ε2。為使參數(shù)ε1、ε2達(dá)到最小值，從而使訓(xùn)練樣本被“夾”在這對控制上下界函數(shù)之間。分別求對w1與b1的偏導(dǎo)，令其為零，即

( 5 )

( 6 )

獲得正確的權(quán)值和偏置后，便求得回歸預(yù)測函數(shù)f(x)。從上述推導(dǎo)可見，LS-TSVM計算過程中并未出現(xiàn)參數(shù)迭代求解，故對應(yīng)的運算量較小。

2 LS-ESPRIT頻率提取方法

ESPRIT原理是正弦信號的頻率以及幅值不會因為時移而發(fā)生變化，故可用該方法估計其頻率成分，進(jìn)而估計出待測信號幅值及初始相位。

假設(shè)s(n)是由K個正弦分量加上觀測噪聲組成

( 7 )

式中：K為諧波成分的數(shù)目；ωk、ak、φk對應(yīng)諧波k的頻率、幅值以及相位；u(n)表示信號中的噪聲成分。忽視直流部分，可將式( 7 )改寫成指數(shù)形式

( 8 )

復(fù)正弦部分可表示為

( 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)

由自相關(guān)、互相關(guān)矩陣定理可得

(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進(jìn)行頻率計算，但該過程需計算樣本協(xié)方差矩陣的逆以及特征值分解，同樣需要較大的運算量。高次諧波廣泛存在于牽引電力系統(tǒng)中，為準(zhǔn)確檢測出諧波成分并且提升算法時效，有必要在得到諧波成分的頻率信息之后再進(jìn)行TSVM分析。

本文為降低ESPRIT計算量，采用基于多級維納濾波器(MSWF)[8-9]對信號進(jìn)行前置降維處理。圖1為MSWF結(jié)構(gòu)。

圖1 多級維納濾波器結(jié)構(gòu)

(20)

逐級反推得到標(biāo)量維納權(quán)序列為

(21)

總的最佳權(quán)為

WMWF=TDWD

(22)

MSWF的前向部分與后向部分均含有迭代，若用D表示迭代數(shù)，降維則是將MN-1維的信號映射為D維度矢量，即對濾波的計算量進(jìn)行了降低。

采用TSVM結(jié)合ESPRIT的諧波檢測算法整體步驟如下：

步驟1用采集得到的時域信號構(gòu)成陣列接收陣

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

式中：M為陣元數(shù)；M+n為總采樣數(shù)目。

步驟2當(dāng)連續(xù)3次分解以后，期望數(shù)據(jù)di與觀測數(shù)據(jù)Xi的互相關(guān)系數(shù)小于閾值時，則降維結(jié)束。根據(jù)式(20)、式(21)計算降維矩陣TD和綜合矩陣WD，利用MSWF降維技術(shù)對空域信號接收矩陣進(jìn)行降維，觀測輸入信號經(jīng)降維矩陣TD處理之后得到的D維向量XD=TDWDS(n)。

步驟3利用前降維后的信號子空間XD取代式(10)中的被分析時空域信號S(n)，根據(jù)式(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)。根據(jù)參數(shù)a1k與a2k可得諧波的幅值A(chǔ)k與相位φk。

(24)

求解a1k與a2k可變換為求式(23)的線性回歸問題，通過TSVM求解回歸系數(shù)w，即可得到每個頻率成分對應(yīng)的Ak與φk。

步驟5采用閾值去除方式，把幅值小于閾值的成分篩掉，從而減少虛假頻率的干擾。

4 綜合性能仿真法分析

這里從算法的魯棒性、精度、時效3個維度對其進(jìn)行測試分析，并與其他算法進(jìn)行對比，用于算法測試的計算機硬件參數(shù)：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點。設(shè)c1=c2=0.000 1，ε1=ε2=0.01，陣元數(shù)200，快拍數(shù)280，信號加高斯白噪聲。

圖2 仿真信號時域波形

本文方法的主要特點是在ESPRIT進(jìn)行頻率估計過程中加入多級維納濾波器進(jìn)行降維，減少運算量。將采用維納濾波器降維的ESPRIT頻率估計方法與常規(guī)ESPRIT頻率估計進(jìn)行對比，用MATLAB編寫程序，計算結(jié)果和計算時長見表1。

表1 降維前后ESPRIT及FFT比較

表1中計算時長采用20次實驗的平均值作為最終結(jié)果。從表1數(shù)據(jù)可以看出，采用維納濾波器的ESPRIT頻率估計方法和常規(guī)ESPRIT頻率估計方法在精度上比較接近，但在計算時長上，采用維納濾波器降維的ESPRIT頻率估計大約需要0.007 219 s，而不采用維納濾波器的ESPRIT頻率估計方法大約需要0.103 042 s，采用維納濾波器降維的ESPRIT能有效減少頻率估計等待時長，提高計算效率。

除了運算效率，表1中我們還將該方法與傳統(tǒng)的FFT頻率計算結(jié)果進(jìn)行了比較，從表1可以看出FFT在計算時長上與降維的ESPRIT方法相差不大，但是在相同的采樣數(shù)據(jù)長度(480點)的情況下，F(xiàn)FT的計算精度遠(yuǎn)不及ESRIT。

本文還研究不同噪聲環(huán)境對降維ESPRIT頻率估計精度的影響，對采集到的信號分別加入30、50、70 dB的高斯白噪聲進(jìn)行計算，其結(jié)果如圖3所示。

圖3 噪聲對計算精度的影響

根據(jù)計算結(jié)果可以看出，在添加30、50、70 dB的高斯白噪聲后，運算結(jié)果差距并不大，故本文的頻率估計算法除了運算效率高以外，同時也具有較好的魯棒性。

用表1得到的頻率信息構(gòu)造最小二乘雙子支持向量機訓(xùn)練的輸入向量，繼而求出諧波的幅值和相位。本文還將TSVM算法和SVM算法進(jìn)行了對比，對比結(jié)果見表2。

表2 TSVM與標(biāo)準(zhǔn)SVM對比

由表2數(shù)據(jù)可知，TSVM與SVM在計算諧波幅值和相位時計算精度差距不大，但TSVM算法在計算時長上占有絕對的優(yōu)勢。TSVM計算時長大約為SVM計算時長的四分之一，這也和理論研究相吻合。TSVM把規(guī)劃問題一分為二后，規(guī)劃問題的復(fù)雜度僅為原來的1/2。假設(shè)樣本數(shù)量為n，SVM的運算量為8n3，則TSVM的運算量為2n3。

運用本文方法完整檢測出諧波的頻率、幅值、相位后，還與矩陣束法[12]進(jìn)行了比較。采用矩陣束算法計算所得結(jié)果見表3。

表3 矩陣束檢測結(jié)果

采用上述兩方法對20組數(shù)據(jù)統(tǒng)計絕對誤差的均方根值，以降低算法對比的隨機性，其結(jié)果見表4。

表4 絕對誤差的均方根值

從表4的結(jié)果分析可知：當(dāng)信噪比為30 dB時，兩種方法的精度較好，但本文算法略高于矩陣束算法。文中方法在降維過程中若MSWF分解的級數(shù)多于待測信號中的諧波數(shù)目，可能引入虛假頻率，在20次實驗中，虛假頻率并不會持續(xù)存在，其幅值常低于0.01。20次連續(xù)仿真中，本文算法最高運行時長為0.009 358 s，最少運行時長為0.005 597 s，運行時長的均方根值為0.007 6 s；而矩陣束法對應(yīng)時長為0.028 0、0.014 5、0.018 7 s。相比而言，矩陣束法明顯比本文算法的時效性差，當(dāng)待檢信號諧波頻率較高時，根據(jù)采樣定理需要提高采樣頻率，此時數(shù)據(jù)點數(shù)增多，由于本文算法進(jìn)行了降維，應(yīng)用價值更高。

5 應(yīng)用

研究列車供電系統(tǒng)諧波監(jiān)測方法不但能時刻監(jiān)測列車電力系統(tǒng)的狀態(tài)，還能用于設(shè)計者檢驗列車諧波特性是否滿足設(shè)計要求。圖4為試驗列車主電路圖。該方案采用的硬件設(shè)備有德國CRONOS-PL以及美國NI公司的數(shù)據(jù)采集卡，下位機由imc devices對電力數(shù)據(jù)采集通道進(jìn)行配置，上位機采用具有數(shù)據(jù)顯示等相關(guān)功能的LabView程序，并對記錄的數(shù)據(jù)采用MATLAB離線分析。

圖5是列車運行速度為300 km/h時采集的牽引電機電流信號，從圖5時域波形可以看出信號中有較為嚴(yán)重的畸變，從采集的信號中截取200 ms進(jìn)行分析，加漢寧窗，得其頻譜圖，如圖6(a)所示，截取40 ms數(shù)據(jù)應(yīng)用本文算法進(jìn)行分析處理，計算諧波參數(shù)，得到頻譜圖，如圖6(b)所示。

為了進(jìn)一步驗證ESPRIT-TSVM方法的實用性，采用該方法對動車電力系統(tǒng)進(jìn)行諧波分析，同時列出了FFT算法的分析結(jié)果。

從圖6可以看出，采用FFT得到的諧波成分與采用本文方法提取的頻譜圖數(shù)值較為接近，在一定程度上也佐證了本文方法的準(zhǔn)確性，實測發(fā)現(xiàn)采用本文方法計算得到的諧波成分予以剔除，優(yōu)于直接采用FFT法。

圖4 試驗列車主電路圖

圖5 列車速度300 km/h電流信號

(a)傅里葉頻譜

(b)本文提取譜圖6 列車頻譜圖

圖7(a)為本方案上位機測試軟件的主界面，主要信息包含各個測點的諧波占比、諧波譜圖等，采用圖表加文字的方式更加清晰；圖7(b)為列車時速380 km時的網(wǎng)側(cè)電流分析結(jié)果，主要包含33次以下的奇次諧波，諧波成分中23、25、27、29、33次諧波所占比重較大，若不加以控制，可能會導(dǎo)致較大干擾電流的產(chǎn)生。

(a)測試軟件主界面

(b)380 km/h網(wǎng)側(cè)電流諧波分布圖7 測試軟件主界面及網(wǎng)側(cè)電流諧波分布

從圖8可以看出，在列車運行速度小于210 km/h時，隨著車速的增加，THD值迅速降低，但仍比要求的10%要高；在列車速度超過210 km/h后THD值變化比較平穩(wěn)，表明此時THD受車速的影響開始變小，但普遍大于8%；在車速大于450 km/h后，THD又呈現(xiàn)出增加的趨勢。

圖8 不同速度下的網(wǎng)側(cè)電流諧波畸變率

6 結(jié)論

本文在支持向量機結(jié)合ESPRIT諧波檢測的方法上進(jìn)行改進(jìn)，提出LS-TSVM結(jié)合多級維納濾波器降維ESPRIT的檢測方法。利用多級維納濾波器降維的ESPRIT快速估計出諧波的頻率，用得到的頻率信息構(gòu)造最小二乘雙子支持向量機的輸入向量得到諧波的幅值與相位。經(jīng)過仿真試驗獲得如下結(jié)論：

(1)在相同的電網(wǎng)諧波數(shù)據(jù)量之下，該算法與傳統(tǒng)的FFT等方法相比具有更好的計算精度。

(2)通過采用多級維納濾波器降維后，降低了程序的運行時長，提高了諧波檢測的效率。

(3)本文方法在不同的噪聲環(huán)境均保持可觀的精度，具有較好的魯棒性。

該方法能準(zhǔn)確檢測出電網(wǎng)諧波的頻率、幅值以及相位，具有較好的實用價值，并且已在某型列車電網(wǎng)中進(jìn)行應(yīng)用，反映良好。

參考文獻(xiàn)：

[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].儀器儀表學(xué)報,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]占勇,丁屹峰,程浩忠,等.電力系統(tǒng)諧波分析的穩(wěn)健支持向量機方法研究[J].中國電機工程學(xué)報,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]丁前軍,王永良,張永順.一種多級維納濾波器的快速實現(xiàn)——迭代相關(guān)相減法[J].通信學(xué)報,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]黃磊,吳順君,張林讓,等.快速子空間分解方法及其維數(shù)的快速估計[J].電子學(xué)報,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]黃磊.快速子空間估計方法研究及其在陣列信號處理中的應(yīng)用[D].西安:西安電子科技大學(xué),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]朱瑞可,李興源,趙睿,等.矩陣束算法在同步電機參數(shù)辨識中的應(yīng)用[J].電力系統(tǒng)自動化,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.