基于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.