考慮系統(tǒng)諧波阻抗改變的諧波責(zé)任定量劃分方法

陳 靜，符 玲，臧天磊，何正友

（西南交通大學(xué) 電氣工程學(xué)院，四川 成都 610031）

0 引言

隨著智能電網(wǎng)的發(fā)展，越來(lái)越多的現(xiàn)代電力電子裝置接入電網(wǎng)，給電網(wǎng)注入了大量諧波［1-2］。而且現(xiàn)代電子產(chǎn)品對(duì)電網(wǎng)諧波的敏感性更強(qiáng)，對(duì)電能質(zhì)量的要求更高，因此對(duì)電網(wǎng)進(jìn)行諧波治理顯得尤為重要。針對(duì)諧波污染問(wèn)題，國(guó)際上提出了一種“獎(jiǎng)懲性方案”，但沒有給出具體的定量規(guī)范。因此，為了保證該方案實(shí)施得有據(jù)可依，必須對(duì)各諧波源的諧波責(zé)任進(jìn)行定量劃分［3-5］。

諧波責(zé)任劃分的關(guān)鍵是系統(tǒng)諧波阻抗的辨識(shí)?，F(xiàn)有的諧波阻抗辨識(shí)的研究可以分為兩大類：波動(dòng)量法［6-9］和線性回歸法［10-14］。 波動(dòng)量法是用諧波電壓波動(dòng)量與電流波動(dòng)量的比值的符號(hào)特征進(jìn)行估計(jì)。回歸法是根據(jù)等效電路中母線處諧波電壓和諧波電流的關(guān)系列出對(duì)應(yīng)方程，求解回歸系數(shù)，從而得到系統(tǒng)諧波阻抗的值。但上述研究方法大部分都基于系統(tǒng)諧波不變的假設(shè)，而考慮到系統(tǒng)阻抗改變的研究甚少。在實(shí)際運(yùn)行的系統(tǒng)中，設(shè)備的投切、無(wú)功補(bǔ)償?shù)母淖兊榷紩?huì)導(dǎo)致系統(tǒng)諧波阻抗改變，而系統(tǒng)諧波阻抗的改變會(huì)引起用戶注入諧波的改變，使用現(xiàn)有方法無(wú)法準(zhǔn)確地劃分諧波源的諧波責(zé)任。文獻(xiàn)［15］討論了系統(tǒng)諧波阻抗的變化對(duì)公共耦合點(diǎn)PCC（Point of Common Coupling）處諧波含量的影響，分X／R恒定和X／R變動(dòng)2種情況，研究諧波阻抗波動(dòng)對(duì)PCC處諧波電流、諧波電壓及諧波功率的影響，定性地分析了系統(tǒng)諧波阻抗變化對(duì)PCC處系統(tǒng)側(cè)和用戶側(cè)諧波貢獻(xiàn)量的影響，但并沒有定量地計(jì)算諧波責(zé)任的大小。因此有必要進(jìn)行系統(tǒng)諧波阻抗改變情況下的諧波責(zé)任定量劃分研究。

在上述2類現(xiàn)有的系統(tǒng)諧波阻抗辨識(shí)的研究中，回歸法的應(yīng)用更廣泛。在線性回歸法的基礎(chǔ)上，研究學(xué)者提出了二元線性回歸［10］、復(fù)線性最小二乘回歸［11］、偏最小二乘回歸［14］等方法。 但這些方法要求自變量是精確無(wú)誤差的，或者是其測(cè)量誤差與因變量的測(cè)量誤差相比可以忽略不計(jì)。但在實(shí)際情況中諧波電流和諧波電壓的測(cè)量量都存在誤差。而整體最小二乘回歸［16］能夠很好地解決這個(gè)問(wèn)題。同時(shí)考慮到數(shù)據(jù)中存在粗差或異常值，本文采用穩(wěn)健整體最小二乘回歸［17］辨識(shí)系統(tǒng)諧波阻抗。

基于上述分析，本文在考慮系統(tǒng)諧波阻抗改變的條件下，首先采用小波變換模極大值法檢測(cè)系統(tǒng)諧波阻抗改變點(diǎn)，以此對(duì)測(cè)量數(shù)據(jù)進(jìn)行分段，在每一段內(nèi)采用穩(wěn)健整體最小二乘回歸法辨識(shí)系統(tǒng)諧波阻抗，并計(jì)算每一段的諧波責(zé)任，最后根據(jù)定義的總諧波責(zé)任指標(biāo)，求得諧波源的總諧波責(zé)任。

1 總諧波責(zé)任指標(biāo)

諧波責(zé)任劃分的前提是建立用于評(píng)價(jià)諧波責(zé)任的指標(biāo)。在系統(tǒng)諧波阻抗改變的條件下需要在傳統(tǒng)諧波責(zé)任算式的基礎(chǔ)上定義新的諧波責(zé)任劃分指標(biāo)。

假設(shè)PCC處接有m個(gè)諧波源負(fù)荷，則PCC處的h次諧波電壓是由這m個(gè)諧波源在PCC處產(chǎn)生的諧波電壓和系統(tǒng)諧波電壓疊加形成的。其電壓相量圖如圖1所示。

圖1 PCC處諧波電壓相量示意圖Fig.1 Schematic diagram of harmonic voltage phasor at PCC

圖1中，為PCC處的h次諧波電壓；為第i個(gè)諧波源單獨(dú)作用時(shí)在PCC處產(chǎn)生的 h 次諧波電壓；γi（i=1，2，…，m）為第 i個(gè)諧波源產(chǎn)生的諧波電壓與PCC處諧波電壓相量的夾角；為系統(tǒng)h次諧波電壓；γ0為系統(tǒng)諧波電壓與PCC處諧波電壓相量的夾角。

圖1中的相量關(guān)系可以表示為：

設(shè)為第i個(gè)諧波源的h次系統(tǒng)諧波阻抗，為第i個(gè)等效諧波電流源，則有，因此由式（1）可得：

根據(jù)諧波責(zé)任的定義式［18］，第i個(gè)諧波源在PCC處產(chǎn)生的h次諧波責(zé)任定義為該諧波源產(chǎn)生的諧波電壓在PCC處諧波電壓上的投影與PCC處諧波電壓模值的比值，即：

在系統(tǒng)諧波阻抗變化的情況下，假設(shè)在關(guān)注時(shí)間段內(nèi)系統(tǒng)諧波阻抗變化了N-1次，即在該時(shí)間段內(nèi)系統(tǒng)諧波阻抗有N個(gè)值。按照時(shí)間順序?qū)?shù)據(jù)點(diǎn)分為N段，每段對(duì)應(yīng)不同的系統(tǒng)諧波阻抗值。則第j（j=1，2，…，N）段數(shù)據(jù)的諧波責(zé)任為，設(shè)該段的時(shí)間長(zhǎng)度為h（j），則定義關(guān)注時(shí)間內(nèi)總的諧波責(zé)任為：

2 基于小波變換模極大值法的數(shù)據(jù)分段

諧波責(zé)任劃分的前提是求取系統(tǒng)諧波阻抗。在實(shí)際系統(tǒng)中，系統(tǒng)諧波阻抗隨著系統(tǒng)運(yùn)行方式、負(fù)荷或無(wú)功補(bǔ)償?shù)母淖兌兓?，即系統(tǒng)諧波阻抗會(huì)隨時(shí)間而發(fā)生階段性的變化。因此，在這種情況下，需要尋求新的諧波阻抗的辨識(shí)方法。

在實(shí)際系統(tǒng)中，諧波電壓或諧波電流的波形并不能完全反映系統(tǒng)諧波阻抗的改變情況，因此本文采用系統(tǒng)諧波阻抗的粗估值來(lái)對(duì)數(shù)據(jù)進(jìn)行分段。實(shí)際情況下系統(tǒng)諧波阻抗的基本改變趨勢(shì)可以用圖2進(jìn)行模擬。 在（tj，tj+1）（j=1，2，…，N-1）時(shí)間段內(nèi)系統(tǒng)諧波阻抗基本不變，現(xiàn)有的系統(tǒng)諧波阻抗求取方法也大多基于這一假設(shè)。但在t1、t2、…、tN-1時(shí)間點(diǎn)上系統(tǒng)諧波阻抗發(fā)生了突變，用現(xiàn)有的方法計(jì)算系統(tǒng)諧波阻抗會(huì)產(chǎn)生較大的偏差。因此，本文首先對(duì)數(shù)據(jù)進(jìn)行分段處理，再在不同段內(nèi)估計(jì)系統(tǒng)諧波阻抗，而分段的關(guān)鍵就是分辨出系統(tǒng)諧波阻抗的改變時(shí)間點(diǎn) t1、t2、…、tN-1。

圖2 系統(tǒng)諧波阻抗改變示意圖Fig.2 Schematic diagram of system harmonic impedance change

在突變點(diǎn)tj時(shí)刻，由于實(shí)際系統(tǒng)的復(fù)雜性以及系統(tǒng)改變暫態(tài)過(guò)程的存在，實(shí)際的系統(tǒng)諧波阻抗的改變不會(huì)呈現(xiàn)出如圖2所示的理想狀態(tài)，而是一個(gè)漸變的過(guò)程，因此系統(tǒng)諧波阻抗改變的界限可能并不明晰，因此需要找到一種有效的方法對(duì)數(shù)據(jù)的改變點(diǎn)進(jìn)行自適應(yīng)檢測(cè)。由于小波變換模極大值在檢測(cè)信號(hào)奇異性方面具有良好性能，本文采用小波變換模極大值法檢測(cè)系統(tǒng)諧波阻抗改變點(diǎn)。

2.1 小波變換模極大值法檢測(cè)原理

小波變換模極大值可以用來(lái)表征信號(hào)的局部奇異性。小波變換模極大值法是根據(jù)不同變換尺度下信號(hào)奇異點(diǎn)與噪聲的模極大值不同的特性來(lái)檢測(cè)奇異值點(diǎn)的。該方法對(duì)輸入信號(hào)的要求低，而且無(wú)需信號(hào)的數(shù)學(xué)模型，是進(jìn)行奇異點(diǎn)檢測(cè)和定位的有力工具［19］。

給定平方可積的信號(hào) x（t），則 x（t）的小波變換定義為：

其中，ψ（t）為一個(gè)基本小波或母小波；a為尺度因子；b 為時(shí)移因子；ψa，b（t）為基本小波通過(guò)平移和伸縮產(chǎn)生的一族函數(shù)，稱為小波基函數(shù)。

當(dāng)小波函數(shù)為平滑函數(shù)的一階導(dǎo)數(shù)時(shí)，信號(hào)小波變換模的局部極值點(diǎn)即對(duì)應(yīng)于信號(hào)的突變點(diǎn)，即小波變換模的局部極大值對(duì)應(yīng)于信號(hào)x（t）的急劇變化點(diǎn)；的局部極小值對(duì)應(yīng)于信號(hào)x（t）的緩慢變化點(diǎn)。因此，求奇異點(diǎn)的位置就可以轉(zhuǎn)化為求WTx的模的局部極大值。

用小波變換對(duì)電力信號(hào)進(jìn)行時(shí)頻分析時(shí)，必須能夠提取信號(hào)的突變成分。最常用的小波基是Daubechies（db）小波，該小波具有正交、時(shí)頻緊支撐、高正則性和具有Mallat快速算法等特點(diǎn)［20］，對(duì)于檢測(cè)信號(hào)的奇異性具有很好的效果。但是db小波具有近似的對(duì)稱性，在信號(hào)重構(gòu)時(shí)會(huì)引起一定的相移。而Symlets（sym）小波是在db小波的基礎(chǔ)上提出的，在保持了db小波基本性能的基礎(chǔ)上提高了小波的對(duì)稱性，減少了重構(gòu)時(shí)的相移［20］。

Symlets函數(shù)系通常表示為 sym N（N=2，3，…，8）的形式。 本文選用 sym4 小波，其尺度函數(shù) φ（t）［20］和小波圖形ψ（t）如圖3所示。

圖3 sym4小波的尺度函數(shù)和小波函數(shù)圖形Fig.3 Scale function waveform and wavelet function waveform of sym4 wavelet

2.2 數(shù)據(jù)分段原理

在實(shí)際系統(tǒng)中，系統(tǒng)諧波阻抗的突變點(diǎn)界限并不明顯，因此本文采用加窗的方法，將對(duì)突變點(diǎn)的辨識(shí)轉(zhuǎn)化為對(duì)突變窗的辨識(shí)，以減少辨識(shí)的誤差。選取窗長(zhǎng)L，根據(jù)式（2），在該窗內(nèi)若系統(tǒng)諧波阻抗值不發(fā)生突變，則對(duì)諧波電壓和諧波電流數(shù)據(jù)進(jìn)行一元線性回歸擬合，擬合曲線的斜率作為系統(tǒng)諧波阻抗的粗略估計(jì)值，而且在時(shí)間段（tj，tj+1）內(nèi)斜率值近似相等；若系統(tǒng)諧波阻抗發(fā)生了突變，則擬合得到的斜率值相較于臨近窗的斜率值也會(huì)發(fā)生突變。由此可以得到系統(tǒng)諧波阻抗的大致變化曲線。

為了確定系統(tǒng)諧波阻抗突變的時(shí)間窗，對(duì)上述系統(tǒng)諧波阻抗的變化曲線進(jìn)行小波變換，其中小波基為sym4小波。由于高頻段可以反映出信號(hào)的突變點(diǎn)，因此取小波變換后得到的高頻段d1和d2進(jìn)行分析。設(shè)定閾值，其中σ為高頻段信號(hào)的標(biāo)準(zhǔn)差，N為采樣點(diǎn)數(shù)據(jù)。在閾值T以下的認(rèn)為是噪聲，在T以上的即認(rèn)為是信號(hào)本身具有的突變。因此，選取在閾值T以上的數(shù)據(jù)點(diǎn)，并確定其對(duì)應(yīng)的時(shí)間，即為系統(tǒng)諧波阻抗發(fā)生突變的時(shí)間窗。

將系統(tǒng)諧波阻抗發(fā)生突變對(duì)應(yīng)的窗內(nèi)的數(shù)據(jù)剔除，對(duì)剔除后的數(shù)據(jù)進(jìn)行分段，并將2個(gè)突變點(diǎn)之間的數(shù)據(jù)歸為同一段。

3 系統(tǒng)諧波阻抗估計(jì)及諧波責(zé)任劃分

用上述原理對(duì)初始測(cè)量數(shù)據(jù)進(jìn)行分段后，每一段內(nèi)的系統(tǒng)諧波阻抗是相同的，因此在每一段內(nèi)分別求取系統(tǒng)諧波阻抗的值?？紤]到自變量和因變量都存在測(cè)量誤差，以及測(cè)量數(shù)據(jù)存在異常值，本文選用穩(wěn)健整體最小二乘回歸法求取系統(tǒng)諧波阻抗。

3.1 穩(wěn)健整體最小二乘回歸原理

設(shè)多元線性回歸方程為：

其中，y 為因變量；x1、x2、…、xp為 p 個(gè)自變量；b0、b1、…、bp為p+1個(gè)待定參數(shù)；自變量和因變量有n次測(cè)量值（n＞p+1）。在傳統(tǒng)的最小二乘法中，假設(shè)x1、x2、…、xp的測(cè)量值沒有誤差，以沿y軸方向的殘差平方和最小來(lái)求取參數(shù) b0、b1、…、bp的值，即目標(biāo)函數(shù)為：

而整體最小二乘回歸是在考慮自變量和因變量均存在誤差的條件下回歸得到模型的參數(shù)。即尋找參數(shù) b0、b1、…、bp使得數(shù)據(jù)點(diǎn)（xi1，xi2，… ，xip，yi）到方程所在超平面的距離的殘差平方和最小，其目標(biāo)函數(shù)為：

記 Xi=［1，xi1，xi2，…，xip］，β0=［b1，b2，…，bp］T，β=［b0，b1，b2，…，bp］T，則目標(biāo)函數(shù)轉(zhuǎn)化為：

上式為多元函數(shù)的最小值求解問(wèn)題，本文用奇異值分解（SVD）法處理，求解原理和過(guò)程見文獻(xiàn)［16］。

由于測(cè)量數(shù)據(jù)中粗差或異常值的存在，在整體最小二乘回歸的基礎(chǔ)上，本文采用穩(wěn)健整體最小二乘回歸［17］來(lái)獲取更穩(wěn)健的參數(shù)估計(jì)值。其具體計(jì)算過(guò)程如下。

（2）計(jì)算觀測(cè)點(diǎn)到超平面的距離：

（3）計(jì)算所有距離的標(biāo)準(zhǔn)差，根據(jù)拉依達(dá)準(zhǔn)則，當(dāng)時(shí)認(rèn)為是異常點(diǎn)，并將之剔除，否則予以保留。

（4）用保留下來(lái)的點(diǎn)重新計(jì)算待求參數(shù)的整體最小二乘回歸解。

（5）重復(fù)步驟（2）和（3）進(jìn)行數(shù)據(jù)的篩選，直到所有觀測(cè)點(diǎn)到超平面的距離都在規(guī)定的閾值內(nèi)。

（6）用整體最小二乘回歸法計(jì)算最終的參數(shù)估計(jì)值，即為穩(wěn)健整體最小二乘回歸的解。

3.2 系統(tǒng)諧波阻抗及諧波責(zé)任劃分原理

根據(jù)文獻(xiàn)［14］，式（1）可以簡(jiǎn)化為：

在實(shí)際系統(tǒng)中測(cè)量量為饋線諧波電流和PCC處諧波電壓，即。所有饋線諧波電流和諧波電壓都有n次測(cè)量值，即有n組樣本點(diǎn)：。 由式（11）可知，以為自變量，以為因變量進(jìn)行回歸，得到的回歸系數(shù)即為系統(tǒng)諧波阻抗的模值。對(duì)這n組樣本點(diǎn)進(jìn)行加窗處理后，用上述數(shù)據(jù)分段原理對(duì)數(shù)據(jù)進(jìn)行分段處理。在每一段內(nèi)用穩(wěn)健整體最小二乘回歸法求取系統(tǒng)諧波阻抗。在求得每段的系統(tǒng)諧波阻抗后，根據(jù)式（3）求取每段諧波責(zé)任，最后用式（4）的定義求得總諧波責(zé)任的值。本文諧波責(zé)任劃分的流程如圖4所示。

圖4 諧波責(zé)任劃分流程圖Fig.4 Flowchart of harmonic contribution determination

4 仿真驗(yàn)證

4.1 三饋線系統(tǒng)算例仿真

在三饋線系統(tǒng)中進(jìn)行諧波責(zé)任劃分仿真分析，PCC處接有3個(gè)諧波源，其諾頓等效電路如圖5所示，PCC左側(cè)為系統(tǒng)諧波源。各電氣元件的初始參數(shù)如表1所示，表中所有電氣量均為標(biāo)幺值。

圖5 三饋線系統(tǒng)Fig.5 Three-feeder system

表1 電氣元件初始參數(shù)Table 1 Initial parameters of electrical components

在三饋線系統(tǒng)中，以諧波源1為研究對(duì)象對(duì)諧波責(zé)任劃分的過(guò)程作具體說(shuō)明。仿真產(chǎn)生n=14400組樣本點(diǎn)數(shù)據(jù)。設(shè)定諧波電流源I1的幅值和相角都以初始值為中心值，作標(biāo)準(zhǔn)差為0.1的正態(tài)波動(dòng)。為了模擬系統(tǒng)諧波阻抗的變化，設(shè)定Z3的值在［3001，5001，9001，11001］處發(fā)生改變，改變后的值分別為［2+j2π，4+j3π，2+j5π，1+j8π］。 根據(jù)電路原理在MATLAB中仿真產(chǎn)生n組PCC處的諧波電壓和諧波電流數(shù)據(jù)。其中諧波電壓和諧波電流的幅值分別如圖 6（a）和圖 6（b）所示，以諧波電流幅值為橫坐標(biāo)，以諧波電壓幅值為縱坐標(biāo)作出散點(diǎn)圖如圖6（c）所示（圖6中各電壓、電流均為標(biāo)幺值）。由圖6（c）可以看出，PCC處的諧波電壓和諧波電流存在某種線性關(guān)系，而且這種關(guān)系可以被分為幾種不同的狀態(tài)，即對(duì)應(yīng)于不同的系統(tǒng)諧波阻抗的狀態(tài)。但是圖中不同狀態(tài)的界限并不明晰，因此本文用加窗結(jié)合小波變換模極大值法檢測(cè)數(shù)據(jù)的改變點(diǎn)。

圖6 PCC處樣本數(shù)據(jù)圖Fig.6 Diagrams of sampled data at PCC

取窗長(zhǎng)L=10，對(duì)上述n個(gè)數(shù)據(jù)點(diǎn)進(jìn)行加窗處理。在每個(gè)窗內(nèi)以諧波電流數(shù)據(jù)為自變量，以諧波電壓數(shù)據(jù)為因變量進(jìn)行最小二乘回歸?；貧w得到的斜率值作為系統(tǒng)諧波阻抗的粗略估計(jì)值，以此可以得到系統(tǒng)諧波阻抗的粗略變化趨勢(shì)圖，如圖7所示（圖中諧波阻抗為標(biāo)幺值，后同）。

圖7 系統(tǒng)諧波阻抗變化趨勢(shì)圖Fig.7 Trend of system harmonic impedance change

用小波變換模極大值法對(duì)上述數(shù)據(jù)點(diǎn)的改變點(diǎn)進(jìn)行檢測(cè)。得到的高頻部分的d1和d2的波形和原波形比較如圖8所示。

圖8 小波變換模極大值法檢測(cè)圖Fig.8 Test results by wavelet transform modulus maximum method

結(jié)合d1、d2的波形，根據(jù)閾值檢測(cè)波形突變點(diǎn)，檢測(cè)超出閾值的點(diǎn)為［301，501，901，1101］，即系統(tǒng)諧波阻抗發(fā)生改變的數(shù)據(jù)窗為（3001，3010）、（5001，5010）、（9001，9010）、（11001，11010）。 與實(shí)驗(yàn)的初始設(shè)定比較，發(fā)現(xiàn)該方法能夠準(zhǔn)確地檢測(cè)出系統(tǒng)諧波阻抗變化所在的窗。

檢測(cè)出系統(tǒng)諧波阻抗改變點(diǎn)后，對(duì)原始測(cè)量數(shù)據(jù)進(jìn)行分段。將系統(tǒng)諧波阻抗發(fā)生改變的數(shù)據(jù)窗內(nèi)的數(shù)據(jù)剔除，將PCC處的諧波電壓和諧波電流分為5 段，分別是（1，3000）、（3011，5000）、（5011，9000）、（9011，11000）、（11 011，14400）。 在這 5 段內(nèi)分別采用穩(wěn)健整體最小二乘回歸估計(jì)系統(tǒng)諧波阻抗的值，計(jì)算結(jié)果如表2所示（表中諧波阻抗為標(biāo)幺值）。由表2看出，穩(wěn)健整體最小二乘回歸的計(jì)算結(jié)果與理論值較為接近。因此可以用穩(wěn)健整體最小二乘回歸計(jì)算系統(tǒng)諧波阻抗的值。

表2 三饋線系統(tǒng)諧波阻抗計(jì)算結(jié)果Table 2 Calculated harmonic impedances of three-feeder system

在計(jì)算得到每一段的系統(tǒng)諧波阻抗的值后，根據(jù)式（3）和式（4）可以分別求出每一段的諧波責(zé)任和總的諧波責(zé)任（THI），結(jié)果如表3所示。

表3 三饋線系統(tǒng)諧波責(zé)任劃分值Table3 Results of harmonic contribution determination for three-feeder system

由表3的結(jié)果可以看出，計(jì)算得到的諧波責(zé)任值和理論值相差不大。在三饋線系統(tǒng)中的仿真是對(duì)本文方法實(shí)施的一個(gè)具體說(shuō)明，從仿真結(jié)果可以證明本文方法的有效性?？紤]到該系統(tǒng)比較簡(jiǎn)單，無(wú)法完全體現(xiàn)實(shí)際系統(tǒng)的復(fù)雜性，本文還選用了較為復(fù)雜的IEEE 13節(jié)點(diǎn)系統(tǒng)進(jìn)行仿真分析，驗(yàn)證本文方法的準(zhǔn)確性。

4.2 IEEE 13節(jié)點(diǎn)系統(tǒng)算例仿真

在如圖9所示的IEEE 13節(jié)點(diǎn)系統(tǒng)中進(jìn)行算例仿真。該系統(tǒng)有2臺(tái)發(fā)電機(jī)、12條支路（包含7臺(tái)變壓器和5條短線路）、7個(gè)PQ節(jié)點(diǎn)。系統(tǒng)各元件和負(fù)荷的參數(shù)參照文獻(xiàn)［21］。以節(jié)點(diǎn)3為關(guān)注母線，設(shè)定節(jié)點(diǎn)8、10、11處負(fù)荷為諧波源負(fù)荷，即諧波源8、10、11，為了模擬系統(tǒng)側(cè)的諧波，在節(jié)點(diǎn)3處也注入了諧波電流。系統(tǒng)諧波阻抗改變可能是由于系統(tǒng)運(yùn)行方式改變、負(fù)荷改變或無(wú)功補(bǔ)償改變引起的。本算例選用無(wú)功補(bǔ)償改變的情況進(jìn)行分析。設(shè)定系統(tǒng)的無(wú)功補(bǔ)償分別為 4 500 kvar、5 500 kvar、6 500 kvar，利用MATLAB進(jìn)行諧波責(zé)任劃分的仿真測(cè)試。

圖9 IEEE 13節(jié)點(diǎn)系統(tǒng)Fig.9 IEEE 13-bus system

仿真時(shí)，將諧波源視為PQ已知的恒定負(fù)荷，用Newton-Raphson法進(jìn)行潮流計(jì)算，生成基波潮流，再根據(jù)文獻(xiàn)［21］的典型諧波電流頻譜計(jì)算出3個(gè)諧波源的注入諧波電流，以5次諧波為例進(jìn)行計(jì)算。設(shè)定注入的諧波電流作標(biāo)準(zhǔn)差為0.2的正態(tài)波動(dòng)，仿真產(chǎn)生n=14 400組數(shù)據(jù)，每4800次仿真改變一次系統(tǒng)無(wú)功補(bǔ)償量。計(jì)算諧波潮流，得到測(cè)量量的模擬數(shù)據(jù)，即節(jié)點(diǎn)3的諧波電壓和3條饋線的諧波電流數(shù)據(jù)。為了使測(cè)量數(shù)據(jù)更加符合實(shí)際情況，在仿真得到的測(cè)量數(shù)據(jù)上加信噪比SNR=30的高斯白噪聲。用諧波源8所接饋線上的諧波電流數(shù)據(jù)和節(jié)點(diǎn)3上的諧波電壓數(shù)據(jù)檢測(cè)系統(tǒng)諧波阻抗的突變點(diǎn)。由于該算例中系統(tǒng)波動(dòng)較大，因此選取窗長(zhǎng)L=30，采用本文方法檢測(cè)系統(tǒng)諧波阻抗改變窗并對(duì)數(shù)據(jù)進(jìn)行分段處理，對(duì)每一段精確求取系統(tǒng)諧波阻抗，進(jìn)而求得諧波源的總諧波責(zé)任。

加窗處理后得到的系統(tǒng)諧波阻抗，即為諧波源8的等效諧波阻抗，其粗略估計(jì)值見圖10。采用小波變換模極大值法檢測(cè)數(shù)據(jù)，得到突變點(diǎn)為［481，961］，則系統(tǒng)諧波阻抗發(fā)生改變的窗為（4 801，4 830）和（9601，9630），與原始設(shè)定相同，因此本文的檢測(cè)方法是有效的。

圖10 諧波源8等效諧波阻抗變化趨勢(shì)圖Fig.10 Trend of equivalent harmonic impedance change of harmonic source 8

由以上結(jié)果可以將數(shù)據(jù)分為（1，4800）、（4831，9600）、（9631，14400）3 段。 對(duì)每段用穩(wěn)健整體最小二乘回歸計(jì)算系統(tǒng)諧波阻抗，并將其計(jì)算結(jié)果與最小二乘LS（Least Squares）回歸［12］、穩(wěn)健最小二乘RLS（Robust Least Squares）回歸［13］和整體最小二乘TLS（Total Least Squares）回歸［16］的計(jì)算結(jié)果比較，3個(gè)諧波源負(fù)荷的計(jì)算結(jié)果如圖11所示。

圖11 系統(tǒng)諧波阻抗計(jì)算結(jié)果比較圖Fig.11 Comparison of calculated systemharmonic impedances

由圖11可以看出，對(duì)諧波源負(fù)荷8、10、11而言，本文方法在每個(gè)時(shí)間段內(nèi)的等效諧波阻抗估計(jì)值都比其他方法更接近于理論值，精確度更高，證明了穩(wěn)健整體最小二乘回歸在諧波阻抗辨識(shí)中更具優(yōu)勢(shì)。

由于諧波責(zé)任劃分的前提是系統(tǒng)諧波阻抗的準(zhǔn)確估計(jì)，在系統(tǒng)諧波阻抗計(jì)算更接近于理論值的情況下，諧波源的諧波責(zé)任劃分結(jié)果也就更準(zhǔn)確。因此，諧波責(zé)任劃分的結(jié)果與系統(tǒng)諧波阻抗的計(jì)算結(jié)果類似，都驗(yàn)證了本文的穩(wěn)健整體最小二乘回歸的計(jì)算結(jié)果優(yōu)于其他方法，在此不再贅述，根據(jù)式（12）計(jì)算得到3個(gè)諧波源的諧波責(zé)任劃分誤差如表4所示。

表4 諧波責(zé)任劃分誤差Table 4 Errors of harmonic contribution determination

其中，ui（j）和 ui_th（j）分別為第 i個(gè)諧波源在第 j段內(nèi)諧波責(zé)任計(jì)算值和理論值。

由表4可以看出，相較于其他方法，本文方法的計(jì)算誤差最小，因此本文方法的準(zhǔn)確性更高。

根據(jù)分段求得的諧波責(zé)任可以求取總的諧波責(zé)任，用本文方法分別計(jì)算3個(gè)諧波源的諧波責(zé)任。另外，將計(jì)算結(jié)果與不分段情況下用穩(wěn)健整體最小二乘回歸計(jì)算得到的總諧波責(zé)任進(jìn)行比較，結(jié)果如表5所示。

表5 總諧波責(zé)任計(jì)算值Table 5 Calculated general harmonic contributions

從表5的結(jié)果可以看出，諧波源負(fù)荷8的總諧波責(zé)任最大，諧波源負(fù)荷11的總諧波責(zé)任最小。不分段情況下的計(jì)算結(jié)果與理論值相差較大，而將分段計(jì)算的結(jié)果與理論值比較，發(fā)現(xiàn)本文方法的計(jì)算結(jié)果與理論值很接近，能準(zhǔn)確地進(jìn)行系統(tǒng)阻抗改變情況下的諧波責(zé)任劃分。另外，不分段情況下主要諧波責(zé)任源被誤判為諧波源10，這主要是因?yàn)椴煌沃g系統(tǒng)諧波阻抗的水平改變比較劇烈，從而導(dǎo)致不分段情況下的計(jì)算值與理論值相差較大。因此，不分段情況下計(jì)算會(huì)導(dǎo)致諧波責(zé)任的劃分結(jié)果誤差大，嚴(yán)重情況下還可能導(dǎo)致主要諧波源的誤判，由此證明了分段計(jì)算的必要性。

5 結(jié)論

在系統(tǒng)諧波阻抗發(fā)生改變的情況下，本文采用小波變換模極大值檢測(cè)法和穩(wěn)健整體最小二乘法結(jié)合的方法劃分諧波源的諧波責(zé)任。通過(guò)算例測(cè)試得到如下結(jié)論：在系統(tǒng)諧波阻抗改變的情況下，采用小波變換模極大值法能夠準(zhǔn)確地檢測(cè)到系統(tǒng)諧波阻抗改變點(diǎn)，以此對(duì)測(cè)量數(shù)據(jù)進(jìn)行準(zhǔn)確地分段；在此基礎(chǔ)上，對(duì)每段的數(shù)據(jù)采用穩(wěn)健整體最小二乘回歸估計(jì)系統(tǒng)諧波阻抗，相較于其他方法能夠獲得更高的精確度。同時(shí)，根據(jù)定義求解總諧波責(zé)任，與不分段情況比較，分段計(jì)算能夠更加準(zhǔn)確地劃分各諧波源的諧波責(zé)任。

［1］趙偉，姜飛，涂春鳴，等.電動(dòng)汽車充電站入網(wǎng)諧波分析［J］.電力自動(dòng)化設(shè)備，2014，34（11）：61-66.ZHAOWei，JIANG Fei，TU Chunming，et al.Harmonic currents of grid-connected EV charging station［J］.Electric Power Automation Equipment，2014，34（11）：61-66.

［2］楊源，林圣，臧天磊，等.基于改進(jìn)稀疏表示法的諧波源定位［J］.電網(wǎng)技術(shù)，2013，37（5）：1279-1284.YANG Yuan，LIN Sheng，ZANG Tianlei，et al.Improved sparse representation-based location of harmonic sources in power system［J］.Power System Technology，2013，37（5）：1279-1284.

［3］許加柱，龐麗忠，張志文，等.多諧波源用戶的諧波責(zé)任分?jǐn)偟亩糠治觯跩］.電力自動(dòng)化設(shè)備，2012，32（12）：38-42.XU Jiazhu，PANG Lizhong，ZHANG Zhiwen，et al.Quantitative analysis for harmonic responsibility proration among multiple harmonic sources［J］.Electric Power Automation Equipment，2012，32（12）：38-42.

［4］王瑜，臧天磊，符玲，等.考慮背景諧波電壓變化的多諧波源諧波責(zé)任劃分［J］.電力系統(tǒng)自動(dòng)化，2015，39（18）：55-61.WANG Yu，ZANG Tianlei，F(xiàn)U Ling，et al.Harmonic contribution partition of multiple harmonic sources considering background harmonic voltage fluctuation［J］.Automation of Electric Power Systems，2015，39（18）：55-61.

［5］馬宏忠，徐剛，宋樹平，等.配電網(wǎng)中諧波電流責(zé)任劃分定量分析［J］.電力自動(dòng)化設(shè)備，2014，34（6）：44-49.MA Hongzhong，XU Gang，SONG Shuping，et al.Quantitative analysis of harmonic current responsibility in distribution network［J］.Electric Power Automation Equipment，2014，34（6）：44-49.

［6］YANG H，PIROTTE P，ROBERT A.Harmonic emission levels of industrial loads statistical assessment［C］∥CIGRE Conference on Large High Voltage Electric Systems.Paris，F(xiàn)rance：CIGRE，1996：36-306.

［7］DAN A M.Identification of individual harmonic sources and evaluation their contribution in the harmonic distortion level［C］∥IEEE Power&Energy Society General Meeting.Calgary，Canada：IEEE，2009：1-6.

［8］龔華麟，肖先勇，劉亞梅，等.基于主導(dǎo)波動(dòng)量篩選原理的用戶諧波發(fā)射水平估計(jì)方法［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2010，30（4）：22-27.GONG Hualin，XIAO Xianyong，LIU Yamei，et al.A method for assessing customer harmonic emission level based on the dominant fluctuation filtering principle［J］.Proceedings of the CSEE，2010，30（4）：22-27.

［9］惠錦，楊洪耕，葉茂清.基于阻抗歸一化趨勢(shì)判別的諧波發(fā)射水平估計(jì)［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2011，31（10）：73-80.HUI Jin，YANG Honggeng，YE Maoqing.Assessing harmonic emission level based on the impedance gathering trend discrimination［J］.Proceedings of the CSEE，2011，31（10）：73-80.

［10］張巍，楊洪耕.基于二元線性回歸的諧波發(fā)射水平估計(jì)方法［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2004，24（6）：54-57.ZHANGWei，YANG Honggeng.A method for assessing harmonic emission level based on binary linear regression［J］.Proceedings of the CSEE，2004，24（6）：54-57.

［11］賈秀芳，華回春，曹東升，等.基于復(fù)線性最小二乘法的諧波責(zé)任定量劃分［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2013，33（4）：149-155.JIA Xiufang，HUA Huichun，CAO Dongsheng，et al.Determining harmonic contributions based on complex least squaresmethod［J］.Proceedings of the CSEE，2013，33（4）：149-155.

［12］MAZIN H E，XUW，HUANG B.Determining the harmonic impacts of multiple harmonic-producing loads［J］.IEEE Transactions on Power Delivery，2011，26（2）：1187-1195.

［13］車權(quán)，楊洪耕.基于穩(wěn)健回歸的諧波發(fā)射水平估計(jì)方法［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2004，24（4）：39-42，49.CHE Quan，YANG Honggeng.Assessing the harmonic emission level based on robust regression method ［J］.Proceedings of the CSEE，2004，24（4）：39-42，49.

［14］李麗，馬宏忠，姜寧，等.基于改進(jìn)偏最小二乘法的系統(tǒng)諧波阻抗及諧波發(fā)射水平估算［J］.電力系統(tǒng)保護(hù)與控制，2011，39（1）：92-95.LI Li，MA Hongzhong，JIANG Ning，et al.Assessing harmonic impedance and the harmonic emission level based on improved partial least-squares regressionmethod［J］.Power System Protection and Control，2011，39（1）：92-95.

［15］任立志，陳勇.諧波阻抗變化對(duì)公共耦合點(diǎn)處諧波貢獻(xiàn)量的影響［J］.電力建設(shè)，2014，35（11）：7-12.REN Lizhi，CHEN Yong.Influence ofharmonic impedance changes on PCC harmonic contribution［J］.Electric Power Construction，2014，35（11）：7-12.

［16］RAHMAN M D，KAI-BOR Y.Total least squares approach for frequency estimation using linear prediction［J］.IEEE Transactions on Acoustics，Speech，and Signal Processing，1987，35（10）：1440-1454.

［17］官云蘭，周世健，張立亭，等.穩(wěn)健整體最小二乘直線擬合［J］.工程勘察，2012，40（2）：60-62.GUAN Yunlan，ZHOU Shijian，ZHANG Liting，et al.A robust method for fitting a line to point clouds based on TLS ［J］.Geotechnical Investigation&Surveying，2012，40（2）：60-62.

［18］華回春，賈秀芳，安海清.背景諧波波動(dòng)情況下的諧波責(zé)任定量計(jì)算［J］.電測(cè)與儀表，2014，51（7）：1-8.HUA Huichun，JIA Xiufang，AN Haiqing.Harmonic contributions determination under the background harmonic fluctuations［J］.Electrical Measurement&Instrumentation，2014，51（7）：1-8.

［19］胡銘，陳珩.基于小波變換模極大值的電能質(zhì)量擾動(dòng)檢測(cè)與定位［J］.電網(wǎng)技術(shù)，2001，25（3）：12-16.HU Ming，CHEN Heng.Detection and location of power quality disturbances using wavelet transform modulusmaxima［J］.Power System Technology，2001，25（3）：12-16.

［20］韋力強(qiáng).基于小波變換的信號(hào)去噪研究［D］.長(zhǎng)沙：湖南大學(xué)，2007.WEI Liqiang.Signal denoising research based on wavalet transformation［D］.Changsha：Hunan University，2007.

［21］BURCH R，CHANG G，GRADY M，et al.Test systems for harmonics modeling and simulation ［J］.IEEE Transactions on Power Delivery，1999，14（2）：579-587.