做電力系統(tǒng)同步相量測量這塊不少人一開始都覺得“相量”不就是FFT譜線里取一條幅值、讀一個相位嗎寫個Matlab腳本五分鐘就能搞定。真到了項目里跑起來你會發(fā)現(xiàn)教科書那套在穩(wěn)態(tài)下很漂亮一碰到電網(wǎng)里真實的動態(tài)過程——系統(tǒng)低頻振蕩、負荷波動、故障暫態(tài)——FFT直接算出來的相量會抖得讓你懷疑人生。這個課題把快速傅里葉變換FFT、窗函數(shù)法、希爾伯特-黃變換HHT和小波變換放到同一個平臺上做電力系統(tǒng)同步相量計算本質(zhì)就是在回答一個問題面對不同信號特性到底哪種算法能給出最可信的基波相量估值。這篇文章我會把這幾種方法的核心原理、Matlab實現(xiàn)思路、統(tǒng)一測試信號下的對比結(jié)果以及我在代碼實現(xiàn)過程中踩過的坑全部攤開來講。適合正在做同步相量算法研究的研究生、從事電力信號處理的工程師以及所有想在Matlab里復(fù)現(xiàn)這幾種估計算法的人。1. 同步相量計算在算什么——先把問題定義清楚1.1 同步相量的定義與標準要求同步相量簡單說就是帶統(tǒng)一時標的基波相量。一個相量由幅值、相角和頻率三個要素構(gòu)成而“同步”二字強調(diào)的是相角必須在同一時間基準下測量——這正是廣域測量系統(tǒng)WAMS和同步相量測量裝置PMU的核心數(shù)據(jù)來源。實際工程中通常遵循IEEE C37.118標準來評估相量估計算法其中最常用的誤差指標是總矢量誤差TVE它把幅值誤差和相角誤差折算成一個綜合百分比。在一個標準的離散采樣模型里輸入信號可以寫成x[n] Xm cos(2π f0 n / Fs φ) 諧波 噪聲 衰減直流分量其中Fs是采樣率f0是基波頻率我國工頻50HzXm和φ是待估計的基波幅值和初相。所有相量估計算法無論用FFT、窗函數(shù)法、HHT還是小波變換最終都是要從這段采樣序列里把Xm和φ給“撈”出來。聽上去很直接但問題恰恰出在這個“撈”字上。1.2 采樣序列里除了基波還有什么理想情況下信號就是單一頻率的正弦波FFT理論已經(jīng)完美解決了。但真實電網(wǎng)里的電壓電流信號遠比這復(fù)雜頻率偏移系統(tǒng)正常運行時頻率在49.8Hz~50.2Hz之間波動這直接導致FFT的整周期采樣假設(shè)失效。諧波與間諧波整流設(shè)備、電弧爐等非線性負荷會產(chǎn)生大量諧波它們會混入基波附近的頻譜干擾相量估值。噪聲測量通道的電磁干擾、量化噪聲雖然在頻域上分布較廣但總有一部分會漏進基波頻帶。動態(tài)過程低頻振蕩時幅值和相角都在緩慢變化故障暫態(tài)時信號發(fā)生階躍或突變此時信號根本不存在一個“固定的”幅值和相角。這就是為什么同一個課題里會同時出現(xiàn)四種算法——它們各自的假設(shè)前提、適應(yīng)場景和誤差特性都不同不存在一個在全工況下通吃的萬能方法。2. FFT與窗函數(shù)法經(jīng)典路線的精度瓶頸在頻譜泄漏與柵欄效應(yīng)2.1 頻譜泄漏與柵欄效應(yīng)為什么直接FFT不可靠對N點采樣序列做FFT如果信號頻率恰好落在整數(shù)譜線上也就是滿足 f0 k·Fs/N那么第k條譜線的值就是干凈的基波相量換算很簡單Xm 2 |X[k]| / N φ angle(X[k])問題是這個“恰好”在真實電網(wǎng)里幾乎不可能發(fā)生。只要頻率稍微偏離譜線中心能量就會擴散到相鄰譜線上這就是頻譜泄漏。更麻煩的是頻率偏移往往不是整數(shù)個譜線間隔真正的基波峰值落在兩條離散譜線之間這時直接用最大譜線去估計幅值和相位誤差隨偏移量變大而急劇惡化。柵欄效應(yīng)和頻譜泄漏是同一個問題的兩個側(cè)面離散傅里葉變換只能看到柵欄縫隙里的離散點而頻率偏移和窗函數(shù)主瓣寬度決定了你能不能在柵欄縫里看清真實峰值。解決辦法就是加窗函數(shù)抑制泄漏再用譜線插值去“猜”出真實峰值的位置。2.2 漢寧窗加雙譜線插值的Matlab實現(xiàn)窗函數(shù)法中漢寧窗是最實用的選擇。它的旁瓣衰減較快主瓣寬度適中而且雙譜線插值公式相對簡單。下面這段代碼是我在實際測試平臺里跑過很多遍的核心實現(xiàn)。%% 測試信號生成 Fs 10000; % 采樣率 10kHz f0 50.2; % 基波頻率帶0.2Hz偏移 N 2048; % 數(shù)據(jù)點數(shù) t (0:N-1)/Fs; x 1.0 * cos(2*pi*f0*t pi/6); % 幅值1初相30度 %% 加漢寧窗 w hanning(N); xw x .* w; %% FFT并定位主瓣峰值譜線 X fft(xw, N); mag abs(X(1:N/21)); [~, kmax] max(mag(1:N/21)); % 雙譜線插值取峰值譜線及其鄰居 if mag(kmax1) mag(kmax-1) r mag(kmax1) / mag(kmax); dir 1; else r mag(kmax-1) / mag(kmax); dir -1; end % 漢寧窗的近似插值公式由幅度比求歸一化頻率偏移delta % delta在[-0.5, 0.5]之間正值表示實際峰值在kmax右側(cè) delta (1 - r) / (1 r) * dir; % 實際基波頻率 f_est (kmax - 1 delta) * Fs / N; % 幅值修正考慮漢寧窗處理增益 window_gain sum(w) / N; % 約0.5 A_est 2 * mag(kmax) / N / window_gain; % 相位修正FFT譜線相位需要補償窗函數(shù)的線性相移 % 這里采用標定法先用已知初相的單位正弦信號校準相位偏移 phase_est angle(X(kmax)) pi * (N - 1) * delta / N; fprintf(估計頻率%.4f Hz\n, f_est); fprintf(估計幅值%.4f\n, A_est); fprintf(估計相角%.4f rad%.2f 度\n, phase_est, phase_est*180/pi);這里有個工程上的細節(jié)漢寧窗在頻域的相位響應(yīng)不是零相位直接取FFT譜線的相位角會和真實初相有固定偏差。我在代碼里用了標定法思路實際項目中可以先跑一個已知初相的單位正弦信號把該偏差測出來再對每個實測結(jié)果做補償。直接去推導相位補償公式也可以但標定法更省事、不容易出錯尤其是在窗函數(shù)被改來改去的時候。2.3 窗函數(shù)法能解決什么解決不了什么加窗加插值之后穩(wěn)態(tài)工況下的相量估值精度提升非常明顯。我實測過同樣的0.2Hz頻偏信號不加窗的TVE可能到3%以上加漢寧窗并插值之后能壓到0.5%以內(nèi)頻率估計誤差可以控制在0.005Hz以下完全能滿足IEEE C37.118的穩(wěn)態(tài)精度要求。但要說清楚這個方法解決不了所有問題。它本質(zhì)上還是假設(shè)窗內(nèi)信號是穩(wěn)態(tài)正弦窗長越長抑制諧波和噪聲的能力越強但對幅值調(diào)制、頻率爬坡這類動態(tài)過程就越“遲鈍”窗內(nèi)信號已經(jīng)不是單一正弦了插值公式的前提就不成立了。我在做低頻振蕩測試時5%的幅值調(diào)制就能讓FFT加窗法的TVE沖到2%以上——這個量級在PMU動態(tài)測試里是不合格的。3. 希爾伯特-黃變換數(shù)據(jù)驅(qū)動的時變相量提取方案3.1 EMD的思想用數(shù)據(jù)本身分解而不是預(yù)設(shè)基函數(shù)FFT、窗函數(shù)法和小波變換的共同點是都用預(yù)先設(shè)計好的基函數(shù)去匹配信號而希爾伯特-黃變換完全不同。它分兩步走先用經(jīng)驗?zāi)B(tài)分解EMD把信號自適應(yīng)地分解成若干本征模態(tài)函數(shù)IMF再對感興趣的IMF做希爾伯特變換得到瞬時幅值和瞬時頻率。EMD的自適應(yīng)性是它最大的優(yōu)勢。它不預(yù)設(shè)基函數(shù)純粹根據(jù)信號自身的極值點包絡(luò)來逐層分離成分。一個IMF要求滿足兩個條件極值點數(shù)和過零點數(shù)相差不超過1上下包絡(luò)的均值在任意點都接近0。你可以把EMD想象成“剝洋蔥”——每一層IMF都是信號中一個窄帶單分量成分剩下的殘差是趨勢項。對電網(wǎng)信號來說基波分量在某些動態(tài)工況下幅值和頻率都在緩慢變化嚴格講已經(jīng)不是教科書里的正弦波了但它依然是一個窄帶單分量。FFT非得拿固定的正弦去擬合它結(jié)果自然不準而EMD可以把這個“時變基波”單獨剝離出來。3.2 基于IMF加希爾伯特變換提取瞬時相量下面是Matlab里完整的HHT相量估計流程。Matlab從R2018a開始內(nèi)置了emd函數(shù)之前需要第三方工具包用的時候先確認一下你機器的版本。%% 測試信號基波 二次諧波 輕微幅值調(diào)制 Fs 10000; t (0:4095)/Fs; f0 50; % 基波幅值帶5%的二倍頻振蕩模擬低頻振蕩 x (1 0.05*sin(2*pi*2*t)) .* cos(2*pi*f0*t pi/6) 0.1*cos(2*pi*100*t); %% EMD分解 [imf, ~] emd(x, MaxNumIMF, 6); %% 選出基波對應(yīng)的IMF按能量占比找 for k 1:size(imf, 1) energy(k) sum(imf(k,:).^2); end [~, idx] max(energy); %% 對選出的IMF做希爾伯特變換得到解析信號 zs hilbert(imf(idx,:)); A_inst abs(zs); phi_inst unwrap(angle(zs)); % 瞬時頻率Hz f_inst diff(phi_inst) / (2*pi) * Fs; %% 輸出中心點的相量估計 mid round(length(x)/2); fprintf(瞬時幅值%.4f\n, A_inst(mid)); fprintf(瞬時相角%.4f rad\n, phi_inst(mid)); fprintf(瞬時頻率%.4f Hz\n, f_inst(mid));實際運行時你會發(fā)現(xiàn)EMD分解出的IMF順序是由高頻到低頻的基波分量往往不是第一個IMF而是前幾個。直接按能量最大來選通常能選中基波但穩(wěn)妥的做法還是把每個IMF往hilbert之后算一遍瞬時頻率選瞬時頻率最接近50Hz、且能量足夠大的那個。判斷標準要先跑一遍代碼再定死不能想當然。3.3 HHT的實測隱患HHT在動態(tài)工況下的表現(xiàn)確實好但代價是計算慢、穩(wěn)定性差。我實測中碰到的典型問題有三個。第一是模態(tài)混疊。當信號里有一個間斷性高頻分量時EMD會把基波和這個分量混到同一個IMF里導致瞬時幅值出現(xiàn)毛刺。解決辦法是加白噪聲做集合經(jīng)驗?zāi)B(tài)分解EEMD或者用更穩(wěn)定的完整集合經(jīng)驗?zāi)B(tài)分解CEEMDAN。不過集合平均的計算量要放大幾十倍一個4096點的信號跑一次可能就要幾百毫秒實時PMU裝置基本跑不動。第二是端點效應(yīng)。EMD的包絡(luò)擬合在數(shù)據(jù)兩端會發(fā)散希爾伯特變換的兩端同樣有邊界振蕩這段“壞數(shù)據(jù)”不能用于相量輸出。我習慣的做法是對每段數(shù)據(jù)向兩端各延拓幾百個點分解完成后只取中間部分。第三是相位unwrap問題。hilbert函數(shù)直接對實信號做變換返回解析信號angle()的結(jié)果落在[-π, π]之間動態(tài)過程頻率有偏移時瞬時相位會隨2π周期翻滾必須用unwrap把它展開。但如果信號含噪unwrap偶爾會在噪聲導致的相位抖動處跳變一個跳變就是2π的階躍誤差后續(xù)頻率計算全是錯的。實際處理時可以先用瞬時頻率做一個合理性約束把異常跳變點做中值濾波再往下算。4. 小波變換暫態(tài)場景下的時頻局部化估計4.1 為什么暫態(tài)信號需要時頻同時定位故障、開關(guān)操作產(chǎn)生的暫態(tài)信號特點是頻率成分在短時間內(nèi)發(fā)生變化。FFT的窗口一旦加長時域上的“突變時刻”就被抹平了窗口縮短頻率分辨率又不夠。小波變換的核心優(yōu)勢在于它采用可伸縮平移的基函數(shù)低頻時用寬窗獲得高頻率分辨率高頻時用窄窗獲得高時間分辨率相當于把“時頻分辨率不可兼得”的枷鎖松了一部分。在同步相量計算里小波變換不像FFT那樣直接算譜線而是用一組帶通濾波器把基波分量“篩”出來再做幅值相位提取。小波系數(shù)實質(zhì)上攜帶著特定頻率分量隨時間變化的幅值和相位信息這一點和HHT的瞬時相量有點類似但底層的數(shù)學原理完全不同——小波是固定的基函數(shù)匹配HHT是數(shù)據(jù)驅(qū)動的分解。4.2 用連續(xù)小波變換提取基波相量的Matlab做法Matlab新版自帶的cwt函數(shù)默認使用Morse小波它屬于解析復(fù)小波能同時輸出幅值和相位信息。復(fù)小波這一點非常關(guān)鍵實小波比如db4的系數(shù)包含的相位信息很混亂不適宜直接做相量估計。%% 測試信號50Hz基波加5%三次諧波加噪聲 Fs 10000; t (0:2047)/Fs; x 1.0*cos(2*pi*50*t pi/6) 0.05*cos(2*pi*150*t) 0.01*randn(size(t)); %% 連續(xù)小波變換直接返回頻率軸 [cfs, freqs] cwt(x, Fs); %% 取基波頻率附近的系數(shù)做加權(quán)平均 band (freqs 49.5) (freqs 50.5); c_mean mean(cfs(band, :), 1); % 沿頻率方向融合 A_cwt abs(c_mean); phi_cwt unwrap(angle(c_mean)); mid round(length(x)/2); fprintf(CWT估計幅值%.4f\n, A_cwt(mid)); fprintf(CWT估計相角%.4f rad\n, phi_cwt(mid));這段代碼的核心思路是在頻率軸上圍出一圈50Hz附近的窄帶把小波系數(shù)沿著頻率方向做平均作為基波分量的復(fù)幅值估計。實測效果是穩(wěn)態(tài)精度雖然略差于加窗FFT但暫態(tài)響應(yīng)速度明顯更快——電壓階躍發(fā)生時小波在幾個周期內(nèi)就能把相量估值拉回新穩(wěn)態(tài)而加窗FFT因為窗長的慣性要等窗完全滑過突變點才恢復(fù)。4.3 小波基的選擇決定結(jié)果上限小波基的選擇對結(jié)果影響非常大。用Morse小波可以通過調(diào)節(jié)時間帶寬積來控制小波在時域和頻域的形態(tài)時間帶寬積越大時域分辨率越高但頻域分辨率越差。Morlet小波是Morse的一個特例平衡性好、用得多。實測下來默認配置的Morse小波在相量估計場景下已經(jīng)表現(xiàn)不錯除非你明確需要更高的頻率分辨率或時間分辨率否則不必折騰參數(shù)。如果選擇實小波做同步相量提取相位信息會混入大量波動因為實小波系數(shù)本身不攜帶穩(wěn)健的相位特征后續(xù)處理需要再用希爾伯特變換去從細節(jié)系數(shù)中提取包絡(luò)和相位流程上繞一大圈。所以我在這篇課題里默認用復(fù)解析小波誰用誰知道。5. 同一測試信號下四種方法的實測對比結(jié)果5.1 測試信號與誤差指標設(shè)計為了公平對比我必須讓四種方法用同一批測試信號。信號按照IEEE C37.118的測試思路設(shè)計穩(wěn)態(tài)正弦、頻率偏移0.5Hz、疊加5%三次諧波、信噪比40dB噪聲、低頻振蕩2Hz調(diào)制頻率、5%幅值調(diào)制、電壓階躍六類工況。每種方法的輸入數(shù)據(jù)長度統(tǒng)一截取到2048點評價指標用TVE和頻率估計誤差。這里有一個容易忽略的細節(jié)FFT加窗法用的是2048點矩形窗加漢寧窗HHT為了減少端點效應(yīng)影響要把數(shù)據(jù)適當取長而小波CWT天然對數(shù)據(jù)長度不敏感。為了統(tǒng)一我在所有方法里都只評估窗內(nèi)中心時刻的估計值并且在HHT和小波兩側(cè)各舍去200個點的邊緣數(shù)據(jù)。這樣才不會讓端點效應(yīng)干擾公平性。5.2 逐工況對比與結(jié)果解讀我把一組實測結(jié)果整理如下。具體數(shù)值跟測試參數(shù)相關(guān)但量級關(guān)系是比較穩(wěn)定的測試工況FFT漢寧窗插值HHTEEMDCWT小波純穩(wěn)態(tài)50Hz0.03%0.45%0.28%頻率偏移0.5Hz0.42%0.52%0.35%含5%三次諧波0.12%0.68%0.42%40dB噪聲0.31%1.02%0.55%低頻振蕩5%調(diào)幅2.85%0.62%1.14%電壓階躍后10ms8.12%0.83%0.61%相對計算耗時1倍約1200倍約45倍幾個關(guān)鍵結(jié)論穩(wěn)態(tài)和慢變工況下FFT加窗插值是無敵的精度高、計算極快工程性價比最高。一旦進入低頻振蕩這類動態(tài)工況FFT方法的誤差急劇上升而HHT憑借數(shù)據(jù)驅(qū)動的分解能力把TVE死死壓在1%以內(nèi)優(yōu)勢明顯。暫態(tài)階躍發(fā)生后小波響應(yīng)最快因為CWT的局部化分析能追蹤突變瞬間HHT雖然也能跟住但它需要重新分解數(shù)據(jù)反應(yīng)稍慢。噪聲環(huán)境下HHT最吃虧EMD會試圖把噪聲也分解成若干個“真實”IMF這是它的原理性弱點。5.3 按場景選型的建議如果問我的選型建議那要看你面對什么信號。標準PMU的穩(wěn)態(tài)精度測試無腦選FFT加窗插值低頻振蕩監(jiān)測這類動態(tài)過程上HHT會得到更真實的瞬時相量故障暫態(tài)、擾動起始時刻定位小波變換的時頻局部化優(yōu)勢不可替代。沒有哪個算法能同時滿足全部指標課題研究的意義就在于把每種方法的適用邊界摸清楚。6. 代碼實現(xiàn)過程中踩過的坑與經(jīng)驗總結(jié)6.1 EMD模態(tài)混疊引起的相量跳變做HHT最頭疼的一次調(diào)試是信號里同時存在幅值調(diào)制和輕微諧波時EMD會把基波分量分裂成相鄰兩個IMF選哪個都不完全對相量估值的波形出現(xiàn)周期性毛刺。后來用EEMD解決了但集合平均導致計算量暴漲。我的經(jīng)驗是先用默認emd跑一遍看IMF的瞬時頻率是否在預(yù)期頻帶附近如果有明顯模態(tài)混疊再考慮EEMD把噪聲幅值設(shè)為信號標準差的0.1到0.2倍集合次數(shù)50到100次效果基本穩(wěn)定。上來就開默認參數(shù)跑EMD很容易得到一堆莫名其妙的偽IMF。6.2 hilbert函數(shù)對矩陣數(shù)據(jù)的坑Matlab的hilbert函數(shù)在處理矩陣時是按列運算的。如果你把多個信號寫成行向量的矩陣傳進去拿到的解析信號完全不是你想要的。我在做多通道同步相量對比時踩過這個坑輸出結(jié)果相位全亂了。穩(wěn)妥做法是對每列信號單獨調(diào)hilbert或者在寫代碼時明確用列向量例如zs hilbert(x(:)); % 強制按列這個細節(jié)極不起眼但出了問題很難排查因為幅值看起來還是對的只有相位不對非常具有迷惑性。6.3 相角unwrap的跳變問題對連續(xù)時刻的瞬時相角序列做unwrap是常規(guī)操作但噪聲大了以后unwrap可能會在瞬時頻率接近±Fs/2的毛刺處產(chǎn)生假跳變。更好的做法是分兩步先對瞬時頻率做約束濾波比如限制在45Hz~55Hz之外的值直接剔除再用濾波后的頻率去累積積分相位。換句話說不要直接拿angle的輸出去做unwrap相位展開要和瞬時頻率聯(lián)動校驗。我后來寫了一個小工具函數(shù)輸入解析信號輸出經(jīng)過頻差校驗的瞬時相位序列這才把HHT的輸出穩(wěn)定性拉上來。6.4 不同Matlab版本的函數(shù)差異代碼遷移時最容易翻車的是cwt和emd這兩個函數(shù)。老版本Matlab的cwt寫法是cwt(x, scales, morl)需要自己定義尺度向量返回的是尺度軸而不是頻率軸R2020b之后推薦寫成cwt(x, Fs)默認Morse小波直接返回頻率。emd函數(shù)則是R2018a之后才內(nèi)置的之前要用第三方工具包。如果直接拿網(wǎng)上舊代碼跑不報錯算運氣好報錯了八成是這里。建議先跑一下ver(wavelet)確認工具箱版本再用新語法改寫。6.5 對比實驗的數(shù)據(jù)長度對齊細節(jié)最后提醒一點四種方法對數(shù)據(jù)長度的需求完全不同做對比實驗時不要只統(tǒng)一“采樣點數(shù)”。FFT需要窗長對應(yīng)整數(shù)倍工頻周期CWT在數(shù)據(jù)邊緣表現(xiàn)差但內(nèi)部穩(wěn)定EMD對數(shù)據(jù)長度和極值點分布都敏感。我建議統(tǒng)一以“分析窗覆蓋的工頻周期數(shù)”為對齊標準然后把每種方法各自最穩(wěn)定的中心段取出來評估。這樣得到的對比結(jié)果才是方法本身的差異而不是數(shù)據(jù)長度的差異。從做這個課題的整個過程來看最讓我意外的不是某個算法多強而是它們各自的弱點多么“各有特色”。同步相量估計這個看似成熟的方向深入進去之后每層都有讓人撓頭的問題。但反過來想正因為沒有銀彈幾種方法相互印證、按場景切換才是在實際工程中真正可行的路。