一種低頻振蕩快速衰減信號(hào)分量參數(shù)的提取方法
【技術(shù)領(lǐng)域】
[0001 ]本發(fā)明涉及一種電力系統(tǒng)的Prony法低頻振蕩的計(jì)算方法，特別是一種Prony法低 頻振蕩快速衰減信號(hào)分量參數(shù)的提取方法。
【背景技術(shù)】
[0002] 電力系統(tǒng)運(yùn)行時(shí)，由于擾動(dòng)會(huì)發(fā)生頻率為0.1Hz~2.5Hz的低頻振蕩，嚴(yán)重影響電 力系統(tǒng)的穩(wěn)定性，危及電網(wǎng)及其相關(guān)設(shè)備的安全運(yùn)行。因此必須嚴(yán)密監(jiān)測(cè)可能出現(xiàn)的低頻 振蕩現(xiàn)象，及時(shí)處理，防止造成破壞系統(tǒng)穩(wěn)定性的后果。Prony分析是目前分析低頻振蕩最 常用最有效的方法之一，通過對(duì)信號(hào)進(jìn)行Prony分析，可以直接得到反映低頻振蕩的幅值、 相位、頻率及衰減因子等參數(shù)。
[0003] Prony算法用一系列的任意頻率、衰減因子、幅值和初相的指數(shù)函數(shù)的線性組合來 擬合一個(gè)函數(shù)，不用再通過頻域響應(yīng)來求解，其計(jì)算量大為減少。也就是說該數(shù)學(xué)模型可以 由一組衰減的正弦分量組成，是一種使用線性方程組來求解非線性問題的分析方法。
[0004] 假設(shè)輸入信號(hào)x(n)有N個(gè)采樣點(diǎn)，即輸入信號(hào)為χ(0)，···，χ(Ν-1)。對(duì)該輸入信號(hào)建 立信號(hào)模型如下：
[0005]
(1)
[0006] 式中，是x(n)的近似值，ρ為信號(hào)模型的階數(shù)，ρ值與信號(hào)分量數(shù)有關(guān)，bdPZl 的表達(dá)式為：
[0007] (2)
[0008] ( ⑶
[0009]式中，Ai為幅值，為初相，cii為衰減因子，fi為頻率，Δ t為采樣時(shí)間間隔。
[0010]為了求式（1)中21的值，由21構(gòu)造如下的特征多項(xiàng)式
[0011]
(4)
[0012] 則21是特征多項(xiàng)式(4)構(gòu)成的如下特征方程的根
[0013]
(5)
[0014] 式中，a()=l。
[0015] 根據(jù)式（1)和式(5)，可以推導(dǎo)出i(?)滿足遞推的常系數(shù)線性差分方程，該差分方 程為：
[0016]
(6)
[0017] 式(6)中是實(shí)際測(cè)量數(shù)據(jù)x(n)的近似值，它們之間存在誤差e(n)，即
[0018]
(7)
[0019]式⑻代入式(7)，得
[0020]
(8)
[0021] 對(duì)式(8)進(jìn)行最小二乘估計(jì)，使得誤差平方彳
i小，則導(dǎo)致一組非線性 方程，求解困難。為了實(shí)現(xiàn)線性估計(jì)，令
[0022]
(9)
[0023] 則式(8)可寫成
[0024] (10) 1=1
[0025] 式(10)寫成矩陣形式為
[0026]
[0027] 式（11)為特征方程系數(shù)求解方程組，它的方程個(gè)數(shù)為N-p，未知數(shù)個(gè)數(shù)為p，一般情 況下方程個(gè)數(shù)多于未知數(shù)個(gè)數(shù)，可以用最小二乘法計(jì)算估計(jì)值。求解式（11)的線性最小二 乘法稱為擴(kuò)展Prony法。 N-1
[0028] 采用最小二乘法，使[I我《) |2最小，得到最小二乘法的法方程 ?.Ρ
[0029]
(12)
[0030]式中，樣本函數(shù)r(i, j)為(下式中上標(biāo)(*)表示共輒）
[0031]
(13) n=p
[0032] 最小誤差能量為
[0033]
(14)
[0034]用高斯消去法解式（12)，就可以求出a=[ai a2…apil'a的精確求解是Prony法 的關(guān)鍵，求出的a代入特征方程(5)，求出z = [Z1 Z2…zP]T后，就可以計(jì)算出信號(hào)分量的頻 率和衰減因子。式(5)是一元高次方程，一般采用QR分解法求解。
[0035]信號(hào)分量的頻率和衰減因子計(jì)算公式如下：
[0036]
(15)
[0037] 式中，arc tan、In、Im、Re分別為反正切函數(shù)、自然對(duì)數(shù)函數(shù)、取復(fù)數(shù)虛部函數(shù)、取復(fù) 數(shù)實(shí)部函數(shù)。
[0038] 計(jì)算出信號(hào)模型的Zl代入式(1)，得到關(guān)于未知數(shù)b的線性方程如下：
[0039] Zb = x (16)
[0040] 式中
[0041]
[0042]
[0043]
[0044]式（17)的Z是NXp維的范德蒙矩陣。由于21各不相同，范德蒙矩陣Z的各列線性獨(dú) 立，為滿列秩的。求解式（16)計(jì)算出b后，就可以計(jì)算信號(hào)分量的幅值和初相了。式(16)也是 方程個(gè)數(shù)多于未知數(shù)個(gè)數(shù)的方程組，仍然用最小二乘法計(jì)算估計(jì)值。
[0045]幅值和初相計(jì)算公式為
[0046]
(20)
[0047] 由于事先不知道低頻振蕩中信號(hào)分量的個(gè)數(shù)，即式（11)中的p不確定，因此需要把 信號(hào)模型的階數(shù)取得大些（即取為Pe>P)，再通過奇異值分解法計(jì)算出方程組系數(shù)矩陣的有 效秩r，即可得到信號(hào)模型的實(shí)際階數(shù)p = r。這樣把式(11)改寫為式(21)
[0048]
[0049] 式(21)的系數(shù)矩陣為擴(kuò)展階矩陣
[0050]
[0051] 擴(kuò)展階矩陣XeSW-pdXpJ介矩陣，用奇異值分解法可以求出該擴(kuò)展階矩陣的有 效秩r，即可得到信號(hào)模型的實(shí)際階數(shù)p = r。
[0052] 如圖1和圖2所示，現(xiàn)有電力系統(tǒng)的Prony法低頻振蕩的分析方法，主要包括以下步 驟：
[0053] A、讀入濾波后的測(cè)量數(shù)據(jù)
[0054]測(cè)量數(shù)據(jù)中含有高頻信號(hào)及噪聲，需先對(duì)測(cè)量數(shù)據(jù)濾波去掉高頻信號(hào)及噪聲，才 能進(jìn)行Prony法低頻振蕩分析。
[0055] B、奇異值分解法信號(hào)模型實(shí)際階數(shù)計(jì)算
[0056] 根據(jù)測(cè)量數(shù)據(jù)形成擴(kuò)展階矩陣，所取信號(hào)模型階數(shù)Pe3大于信號(hào)模型實(shí)際階數(shù)p (pe取滿足< LM2」的整數(shù)，符號(hào)" U "為向下取整符號(hào)），然后采用奇異值分解 法對(duì)擴(kuò)展階矩陣)yi行計(jì)算，求出擴(kuò)展階矩陣的所有奇異值，此奇異值為非負(fù)的，并按下列 順序排列：
[0057] on > 〇22 > ··· > 〇hh > 0 (23)
[0058] 式中
[0059] h=min(N-pe ,pe) (24)
[0060] 求出所有奇異值并排順后，按下式對(duì)奇異值歸一化
[0061] 〇kk〇 = okk/〇n,l <k<h (25)
[0062] 選擇一個(gè)較小的正數(shù)作為閾值，使〇kkQ大于此閾值的最大整數(shù)k取為擴(kuò)展階矩陣Xe 的有效秩r，則信號(hào)模型的實(shí)際階數(shù)p = r。
[0063] C、特征方程系數(shù)計(jì)算
[0064] 計(jì)算時(shí)，一般需要把信號(hào)模型的階數(shù)m取得比信號(hào)模型的實(shí)際階數(shù)p大一些，否則 計(jì)算結(jié)果的精度很差。取方程未知數(shù)個(gè)數(shù)m大于信號(hào)模型實(shí)際階數(shù)p(m取滿足 UW2j的整數(shù)，符號(hào)"U"為向下取整符號(hào)），重新列方程（η)得到新的特征方 程系數(shù)求解方程組如下：
[0065]
[0066] 此方程組的方程數(shù)(N-m)大于未知數(shù)個(gè)數(shù)m，可以采用最小二乘法對(duì)未知數(shù)進(jìn)行估 計(jì)，計(jì)算出對(duì)應(yīng)特征方程的系數(shù)。由于方程(26)是線性方程，形成相應(yīng)法方程后采用高斯消 去法求解。
[0067] D、頻率和衰減因子計(jì)算
[0068] 信號(hào)模型的階數(shù)為m時(shí)的特征方程為：
[0069]
(27)
[0070] 式中，加二！。
[0071] 式(27)是一元高次方程，其解為實(shí)數(shù)或共輒復(fù)數(shù)，一般采用QR分解法求解，得到z 后，就可以計(jì)算出信號(hào)分量的頻率和衰減因子。
[0072] E、幅值和初相計(jì)算
[0073] 對(duì)應(yīng)式(26 )，求解信號(hào)模型的系數(shù)bi的公式(17)和（18)相應(yīng)變?yōu)椋?br>[0074]
(28)
[0075] b = [bi b2 ··· bm]T (29)
[0076] 把求出的z代入式(28)，并根據(jù)式（16)計(jì)算出b后，就可以計(jì)算信號(hào)分量的幅值和 初相。式（16)的方程數(shù)也大于未知數(shù)個(gè)數(shù)，仍然可以用最小二乘法對(duì)未知數(shù)進(jìn)行估計(jì)。計(jì)算 出信號(hào)分量按幅值由大到小排序，取前P個(gè)信號(hào)分量作為最終的信號(hào)分量結(jié)果。
[0077] 計(jì)算出的信號(hào)分量為如下復(fù)指數(shù)函數(shù)，
[0078]
(3:0)
[0079] -般情況下，復(fù)指數(shù)函數(shù)成對(duì)出現(xiàn)，即復(fù)指數(shù)函數(shù)和它的共輒值成對(duì)出現(xiàn)。函數(shù) VKt)的共輒值為
[0080]
(3 1 )
[0081 ]信號(hào)分量Y i(t)和信號(hào)分量x〃i(t)疊加為正弦分量如下：
[0082]
(32)
[0083] 由于計(jì)算特征方程系數(shù)矩陣時(shí)，未知數(shù)個(gè)數(shù)m大于信號(hào)模型實(shí)際階數(shù)p，計(jì)算出信 號(hào)分量比實(shí)際的信號(hào)分量多。這些多出的分量幅值很小，可以根據(jù)幅值剔除這些多余的分 量。
[0084] 以上傳統(tǒng)分析方法求頻率和衰減因子時(shí)，如果低頻振蕩中含有衰減較快的分量， 分析就非常困難，分析得到的衰減較快信號(hào)分量的參數(shù)誤差很大，甚至分析結(jié)果中不包含 某些衰減較快的信號(hào)分量，丟失這些衰減較快信號(hào)分量的信息。

【發(fā)明內(nèi)容】

[0085]為解決現(xiàn)有技術(shù)存在的上述問題，本發(fā)明要提出一種Prony法低頻振蕩快速衰減 信號(hào)分量參數(shù)的提取方法，以便解決準(zhǔn)確提取低頻振蕩中快速衰減分量的問題。
[0086]為了實(shí)現(xiàn)上述目的，本發(fā)明提出了采用Prony法分段分析方法來解決準(zhǔn)確提取低 頻振蕩中快速衰減分量的問題。
[0087] 本發(fā)明的技術(shù)方案如下：一種Prony法低頻振蕩快速衰減信號(hào)分量參數(shù)的提取方 法，包括以下步驟：
[0088] A、讀入濾波后的測(cè)量數(shù)據(jù)
[0089]按照采樣間隔Δ t和采樣時(shí)長(zhǎng)T讀入濾波后的測(cè)量數(shù)據(jù)，所述的采樣時(shí)長(zhǎng)T>10s，采 樣間隔 At = 0.