秦佳峰，周 超，林 穎，白德盟，鄭文杰

（國網(wǎng)山東省電力公司電力科學(xué)研究院，山東 濟南 250003）

0 引言

隨著電力系統(tǒng)信息化程度的提高，電網(wǎng)的運行和監(jiān)測數(shù)據(jù)均呈指數(shù)級增長［1-3］。巨量的數(shù)據(jù)為智能調(diào)度和決策提供了信息基礎(chǔ)，但也對數(shù)據(jù)處理手段提出了更高的要求［4-6］。變壓器作為電網(wǎng)中的重要一環(huán)，其運行數(shù)據(jù)對后續(xù)的設(shè)備狀態(tài)評估至關(guān)重要，但是由于設(shè)備故障、通信故障、操作失誤等原因，實際取得的運行數(shù)據(jù)往往存在缺失的情況，從而降低后續(xù)數(shù)據(jù)挖掘算法的性能。因此，修復(fù)缺失的變壓器運行數(shù)據(jù)對提升智能電網(wǎng)的運行效率，提高電網(wǎng)系統(tǒng)的可靠性和自愈能力至關(guān)重要。

關(guān)于電力系統(tǒng)缺失數(shù)據(jù)的修復(fù)，目前已有多種方法［7-9］。張曉星等［7］基于聚類的思想計算電力負荷曲線的特征函數(shù)，利用該特征函數(shù)修復(fù)缺失數(shù)據(jù)，但這種方法對數(shù)據(jù)的日曲線相似度要求很高。嚴(yán)英杰等［8］將輸變電設(shè)備狀態(tài)數(shù)據(jù)看成時間序列，用差分整合移動平均自回歸模型（Autoregressive Integrated Moving Average，ARIMA）擬合并迭代檢驗的方法修復(fù)缺失數(shù)據(jù)，但這種方法利用的信息較少，且不適合連續(xù)缺失的情況。劉沅昆等［9］提出了一種基于Pearson 相關(guān)系數(shù)的歷史數(shù)據(jù)挖掘恢復(fù)方法，但這種方法本質(zhì)上是在做線性回歸，很難刻畫出數(shù)據(jù)的實際波動情況。

通過分析真實的變壓器運行數(shù)據(jù)，希望找到一種缺失數(shù)據(jù)修復(fù)方法，不僅能利用不同日曲線之間的相似性，并且適用不同的缺失點分布情況。函數(shù)型主成分分析（Functional Principal Component Analysis，F(xiàn)PCA）作為一種高維特征提取方法，能夠通過對同一觀測對象的重復(fù)測量，在低維函數(shù)空間上表示出數(shù)據(jù)最主要的波動情況［10-11］，能夠很好地擬合出數(shù)據(jù)波動的整體趨勢。小波變換作為信號處理的常用方法，能夠通過對函數(shù)的多尺度細分，聚焦函數(shù)的局部細節(jié)，提取出其中的有效信息［12］。

因此，結(jié)合FPCA和小波變換處理函數(shù)型數(shù)據(jù)的優(yōu)勢，提出了一種變壓器缺失數(shù)據(jù)修復(fù)方法。該方法能夠自動提取數(shù)據(jù)波動的特征，并使恢復(fù)所得的數(shù)據(jù)在局部上更符合實際情形［13］。

1 FPCA原理

1.1 根據(jù)FPCA構(gòu)建模型

具體到變壓器的運行數(shù)據(jù)修復(fù)，考慮其中一類運行數(shù)據(jù)，可以將第i天內(nèi)變壓器的該類運行數(shù)據(jù)視為一個觀測對象Yi，傳感器在第j個觀測點對應(yīng)的時間點tj測得的數(shù)值Yij為該時刻的觀測值，記第i天的觀測點數(shù)量為ni。由于同一個變壓器在連續(xù)的幾天內(nèi)運行數(shù)據(jù)的波動趨勢有一定的相似性，因此可以將N天的數(shù)據(jù)Yij，1 ≤i≤N，1 ≤j≤ni視為一組縱向數(shù)據(jù)進行分析。首先，利用FPCA 方法估計函數(shù)型數(shù)據(jù)的低頻部分Xi(t)。根據(jù)Karhunen-Loève 表示［14］，Xi(t)可表示為

式中：μ(t)為均值函數(shù)；φk(t)為協(xié)方差曲面的特征函數(shù)，即主成分；αik為Xi(t)在主成分上的得分。

利用FPCA方法分別估計均值函數(shù)μ(t)、特征函數(shù)φk(t)、系數(shù)αik。由于通過FPCA方法估計的Xi(t)忽略了數(shù)據(jù)高頻部分的信息，利用小波刻畫函數(shù)數(shù)據(jù)細節(jié)的能力，通過對殘差函數(shù)ξi(t)做小波變換，提取出其中高頻部分的信息。因此整體模型可以改寫為

關(guān)于殘差函數(shù)ξi(t)的處理將在下文敘述。同時，將第i天的第j個測量時間點tj記為Tij，將觀測值Yi(tj)記為Yij。

1.2 均值函數(shù)及協(xié)方差曲面估計

采用加權(quán)最小二乘的方法估計均值函數(shù)μ(t)。對某一固定的時間點t，假設(shè)β0(t)和β1(t)為兩個參數(shù)，其中β0(t)為均值函數(shù)μ(t)的估計，β1(t)為線性修正函數(shù)的斜率，因此Yij-β1(t-Tij)即為t時刻Yij對應(yīng)的修正值。取高斯核函數(shù)為權(quán)重值，因此，可通過求解如下極小化問題得到對均值函數(shù)μ(t)的估計為

式中：β0、β1為關(guān)于時間t的參數(shù)函數(shù)；hμ是帶寬，為了使算法有更強的自適應(yīng)性，采用廣義交叉檢驗的方法選擇帶寬；k(·)是高斯核函數(shù)，其表達式為

通過求解式（3）中的極小化問題，得到對均值函數(shù)μ(t)的估計為

對協(xié)方差曲面的估計同樣可以采用加權(quán)最小二乘的方法。與估計均值函數(shù)μ(t)一樣，對某一組固定的時間點(s，t)，β0(s，t)，β11(s，t)和β12(s，t)都是待估計的參數(shù)，其中β0(s，t)為對協(xié)方差曲面G(s，t)的估計，β11(s，t)和β12(s，t)分別為兩個方向上線性修正函數(shù)的斜率。

首先估計在s≠t時的協(xié)方差曲面，可以通過求解如下極小化問題得到，即為

式中：hG為帶寬，可通過廣義交叉檢驗的方法選擇帶寬；k2(·，·)為一個二元高斯核函數(shù)。k2(·，·)和Gi(Tij，Til)的表達式分別為：

通過求解式（6）中的極小化問題，得到在s≠t時協(xié)方差曲面的估計為

當(dāng)s=t時，由于協(xié)方差曲面在垂直對角線方向上的形狀更接近于二次曲線［15］，在垂直對角線方向上用二次函數(shù)修正，而在對角線方向上仍采用線性函數(shù)修正，將坐標(biāo)軸順時針旋轉(zhuǎn)45°，即為

結(jié)合式（9）和式（11），可得對協(xié)方差曲面的估計(s，t)。

1.3 主成分的選取及主成分得分的估計

主成分是協(xié)方差曲面的特征函數(shù)，因此可以通過求解式（13）的積分方程得到對特征函數(shù)的估計。

具體的，通過將協(xié)方差曲面離散化，得到K組特征值和特征函數(shù)1 ≤k≤K，需要選取其中最重要的k個主成分。為此，定義方差解釋比（Fraction of Variance Explained，F(xiàn)VE）為

給定一個閾值θ，選擇前k個最大的特征值和對應(yīng)的特征函數(shù)，使得AFVE(k)≥θ，完成對主成分的選取。

最后用條件期望來估計主成分得分αik，記

其中δjl在j=l時為1，j≠l時為0，Γ為t的定義域為

該函數(shù)能夠較好地近似函數(shù)型數(shù)據(jù)的整體波動趨勢，但是在局部細節(jié)方面存在過于平滑的現(xiàn)象，因此需要分析真實值與估計函數(shù)做差得到的殘差函數(shù)ξi(t)，通過用小波變換提取其中的高頻有效信息，得到更為精確的估計。

2 小波變換原理

2.1 根據(jù)小波變換構(gòu)建模型

假設(shè)殘差函數(shù)ξi(t)是一個包含噪聲的一維信號，即ξi(t)可以表示為

式中：?i(t)為真實信號；ei(t)為噪聲。將ξi(t)看成是一個閉區(qū)間上的有界連續(xù)函數(shù)，因此可進一步假設(shè)ξi(t)是平方可積的，根據(jù)離散小波變換的理論［16］，ξi(t)展開為

式中：Φ(t)和ψ(t)分別為尺度函數(shù)和小波基函數(shù)；cj0(k)和dj(k)分別為對應(yīng)尺度下的系數(shù)。

由于實際應(yīng)用中，得到的ξi(t)是一列離散的采樣點，根據(jù)采樣定理，尺度無法無限細分下去，因此具體計算中將ξi(t)展開為

式中：j0=0；j1為分解層數(shù)。

2.2 去噪方案的選取

由于小波基、分解層數(shù)和閾值方案的選擇會對去噪效果產(chǎn)生影響，對不同類型的數(shù)據(jù)可能需要選擇不同的參數(shù)進行去噪，隨機選取真實數(shù)據(jù)點中的10%作為測試集，根據(jù)不同參數(shù)組合在測試集上的表現(xiàn)選擇最優(yōu)的去噪方案。

可供 選擇的 小波基［17-19］包 括db1-8、sym2-6、coif1-4，分解層數(shù)為1～8，閾值方案包括閾值選取規(guī)則和小于閾值小波系數(shù)的處理方法兩個部分。

2.2.1 閾值選取規(guī)則

閾值選取規(guī)則主要包括四種方式：

1）無偏風(fēng)險估計閾值（rigrsure），算法為：

b）若Aeta＜Acrit，則λheur=λsqt；若Aeta≥Acrit，則λheur=min{λrigr，λsqt}。

4）極大極小閾值（minimax）為

2.2.2 小于閾值小波系數(shù)的處理方法

小于閾值小波系數(shù)的處理方法主要包括兩種：

1）硬閾值處理為

2）軟閾值處理為

3 變壓器運行數(shù)據(jù)缺失點修復(fù)方案

變壓器運行數(shù)據(jù)缺失點修復(fù)方案主要包括三個部分，首先用FPCA方法得到觀測對象的初次估計函數(shù)；然后估計殘差函數(shù)，并用交叉驗證的方法找到最優(yōu)的小波去噪方案，得到對殘差函數(shù)的估計；最后將殘差函數(shù)的估計作為修正函數(shù)，與FPCA初次估計函數(shù)結(jié)合，得到觀測對象在整個定義域中的估計值，用該估計值完成缺失點的修復(fù)。

以變壓器第i天內(nèi)某一類運行數(shù)據(jù)作為一個觀測對象Yi，考慮連續(xù)n天內(nèi)同一變電站的同一類數(shù)據(jù)，其中1 ≤d≤n。將第d天的觀測對象Yd的真實觀測值集合記為Sd={Yd1，Yd2，···，Ydnd}，nd為連續(xù)觀測n天中的第d天，計算第d天的估計函數(shù)，包括9個步驟。

步驟1）隨機選取其中10%的點作為測試集，記Ud={Ydt1，Ydt2，···，Ydtp}，p=[0.1×nd]，其 余90% 的點作為訓(xùn)練集，記為Vd={Ydtp+1，Ydtp+2，···，Ydtnd}。以第d天的訓(xùn)練集Vd，以及除第d天外所有的真實測量值作為訓(xùn)練集，用FPCA 方法建模，其中t的單位為分鐘，考慮到實際情況，所有可能的觀測點都是整數(shù)分鐘。

步驟3）由式（6），得到s≠t時的協(xié)方差曲面的估計函數(shù)，再由式（11），得到s=t時協(xié)方差曲面的估計函數(shù)，從而得到整個協(xié)方差曲面的估計函數(shù)(s，t)。

步驟4）根據(jù)式（13），通過將協(xié)方差曲面離散化的方法，得到K組特征函數(shù)和對應(yīng)的特征值的估計選取前k個最大的特征值及其對應(yīng)的特征函數(shù)，使得式（14）定義的AFVE(k)大于給定的閾值θ。得到k個主成分的估計

步驟5）根據(jù)式（16），估計出觀測對象Yd在主成分(t)上的得分。綜合步驟1）—步驟5），得到對觀測對象Yd的初次估計函數(shù)，即低頻部分的估計為

步驟6）估計殘差函數(shù)(t)。為方便敘述，將前述Sd中的觀測時間點記為，同理，記，tdnd}。當(dāng)t∈時，則殘差函數(shù)為當(dāng)t?時，補充定義td0=1，td(nd+1)=1440，則必存在p，使得t∈(tdp，td(p+1))。定義則當(dāng)Agap較小時，可用移動平均來估計該點的殘差函數(shù)的值，當(dāng)Agap較大時，對殘差函數(shù)的估計意義不大，反而可能造成最終估計的錯誤，因此直接置為零。綜上，對殘差函數(shù)(t)的估計為

步驟7）利用MATLAB軟件中的小波工具箱對殘差函數(shù)進行小波變換，并選擇最優(yōu)的去噪方案。

首先選定一個去噪方案，用wden 函數(shù)進行小波去噪。該函數(shù)的四個參數(shù)“wname”、“n”、“tptr”、“sorh”分別對應(yīng)于小波基、分解層數(shù)、閾值選取規(guī)則以及軟硬閾值。最后一個參數(shù)“scal”表示定義的閾值是否需要重新調(diào)整，默認為“one”，即不用重新調(diào)整。輸入去噪方案對應(yīng)的參數(shù)以及信號(t)之后，就可得到去噪后的殘差函數(shù)(t)。

選取不同的去噪方案，重復(fù)步驟7），記錄最小的RRMSEi及對應(yīng)的去噪方案。

步驟8）將訓(xùn)練集擴充為整個觀測集，重復(fù)步驟1）—步驟5），更新初次估計函數(shù)(t)，再更新對殘差函數(shù)(t)的估計。即用重復(fù)步驟6），得到更新后的殘差函數(shù)為

這樣就完成了對觀測對象Yd的修復(fù)。

4 算例分析

以某220 kV 主變壓器某月的高壓側(cè)-I 數(shù)據(jù)為例，對本文所提缺失點修復(fù)方法進行測試?？紤]到實際應(yīng)用中既存在零散的隨機缺失點，又存在連續(xù)缺失的情況，因此分別對這兩種情況做了測試。

首先考慮缺失點零散分布的情況。以該月前十天的數(shù)據(jù)為例進行分析，1 日實際獲取了997 個觀測點，隨機選取其中100 個點作為測試集，假設(shè)這些觀測點缺失，以此為參照，評價修復(fù)效果。

經(jīng)過步驟1）—步驟5）后，得到的初次估計函數(shù)如圖1 所示。真實測量值及其變化趨勢用黑點及黑色虛線連接，假設(shè)缺失的100 個點用藍色星狀點表示，初次估計函數(shù)用紅色實線表示。

圖1 初次估計函數(shù)

可以看到，初次估計函數(shù)已經(jīng)能夠較好地擬合真實觀測值，但是在局部細節(jié)上還有欠缺。在步驟1）—步驟5）的基礎(chǔ)上，繼續(xù)實施步驟6）—9）后，得到的最終估計函數(shù)如圖2 所示。黑點和藍色星狀點的含義同圖1，紅色實線表示最終估計函數(shù)。

圖2 最終估計函數(shù)

結(jié)果表明，在對初次估計函數(shù)加上修正函數(shù)之后，最終估計函數(shù)能夠很好地擬合真實測量值。以測試集上的RRMSEd作為指標(biāo)，將這種方法同三次樣條插值方法進行對比。重復(fù)進行了10 次實驗，每次的測試集都是隨機選取，對比兩種方法的均方根誤差（Root Mean Squard Error，RMSE）如表1所示。

表1 高壓側(cè)-I上的RMSE對比

結(jié)果顯示，經(jīng)本文方法得出的估計函數(shù)在測試集上的RMSE 要明顯小于三次樣條插值法的RMSE，10次實驗平均能夠使RMSE下降21.3%，本文的缺失點修復(fù)方法更加穩(wěn)定可靠。

當(dāng)缺失點為大塊的連續(xù)缺失時，普通的插值方法已經(jīng)沒法給出一個合理的估計，反而會添加錯誤數(shù)據(jù)。因此仍然以該月1 日至10 日的數(shù)據(jù)作為一組進行分析，假設(shè)第700～799 個觀測值缺失，利用所述的方法，對該段的缺失值進行恢復(fù)?；謴?fù)的結(jié)果如圖3所示。

圖3 連續(xù)缺失時的恢復(fù)曲線

為了檢驗在所考慮數(shù)據(jù)集上的整體效果，測試該方法應(yīng)用于其他11 類數(shù)據(jù)上的效果，同高壓側(cè)-I的做法一樣，隨機選取100 個真實觀測值作為測試集，重復(fù)實驗10 次，分別計算本文方法和三次樣條方法在測試集上的平均RMSE，結(jié)果如表2 所示。測試結(jié)果顯示，本文方法在不同類型的變壓器運行數(shù)據(jù)恢復(fù)上均有較好的表現(xiàn)。對三次樣條差值方法有很大程度的提高。除高壓側(cè)-Q與中壓側(cè)-Q外，RMSE下降的百分比均超過15%。

表2 12類數(shù)據(jù)上的RMSE對比

5 結(jié)語

利用FPCA 對數(shù)字信號整體特征的提取能力和小波變換對局部信息的提取能力，對變壓器運行數(shù)據(jù)的缺失數(shù)據(jù)修復(fù)這一實際問題，提出了基于FPCA和小波變換的變壓器運行數(shù)據(jù)缺失數(shù)據(jù)修復(fù)方法。

使用FPCA 方法能夠?qū)崿F(xiàn)對需修復(fù)數(shù)據(jù)集的初次估計函數(shù)；并利用交叉驗證方法找到最優(yōu)的小波去噪方案，得到對殘差函數(shù)的估計；通過將殘差函數(shù)的估計作為修正函數(shù)，與FPCA 初次估計函數(shù)結(jié)合，完成缺失點的修復(fù)。

該方法能夠根據(jù)少量歷史數(shù)據(jù)自動尋找數(shù)據(jù)的波動特征，并在此基礎(chǔ)上對局部細節(jié)進行優(yōu)化，從而完成對缺失數(shù)據(jù)的恢復(fù)。該方法對缺失點是連續(xù)分布還是離散分布沒有要求，具有適用性廣、恢復(fù)精度高的特點。