分形插值在風(fēng)速時間序列中的應(yīng)用

郭秀婷，朱昶勝，張生財，趙奎鵬

（1.蘭州理工大學(xué)計算機與通信學(xué)院，蘭州 730050；2.蘭州理工大學(xué)理學(xué)院，蘭州 730050）

0 引言

可再生能源發(fā)電具有不依賴化石燃料、減少溫室氣體排放、保護環(huán)境等一系列優(yōu)點，因此得到了快速發(fā)展；其中風(fēng)能作為一種低成本、分布域廣的清潔能源得到了各國的高度重視［1］。但風(fēng)的間歇性和隨機性給風(fēng)力發(fā)電的規(guī)劃及管理帶來了很多困擾［2-3］，而風(fēng)速預(yù)測和風(fēng)電功率預(yù)測是應(yīng)對風(fēng)電大規(guī)模并網(wǎng)運行問題的重要手段，因此風(fēng)速預(yù)測和風(fēng)電功率預(yù)測成為全球研究的熱點［4-6］。受天氣因素、地理環(huán)境因素、傳感器精度、電磁干擾以及通信故障等因素影響，采集到的風(fēng)速數(shù)據(jù)會出現(xiàn)一定的異常值甚至一段時間的缺失值，但在研究過程中這些異常值和缺失值卻至關(guān)重要，因此研究適用于風(fēng)速的異常值檢測及插值算法具有重要的意義［7-8］。

陳偉等［8］在處理風(fēng)電異常數(shù)據(jù)時采用三次樣條插值進行數(shù)據(jù)重構(gòu)，但自然風(fēng)速時間序列是一個波動強烈的非線性時間序列，Lagrange、三次樣條和線性插值等傳統(tǒng)插值算法中，任意兩個信息點之間用直線或者光滑的曲線連接，掩蓋了相鄰插值點間的局部變化特征，因此用于波動強烈的風(fēng)速時間序列效果并不理想；武艷強等［9］提出的多點三次樣條插值算法對于頻率相對穩(wěn)定的低頻和中頻曲線插值效果比較理想，但對于包含多種頻率混合的非線性風(fēng)速曲線卻不適合；楊茂等［10］提出用支持向量機重構(gòu)風(fēng)電數(shù)據(jù)，但由于參數(shù)的敏感性因素導(dǎo)致誤差較大；陳偉等［11］提出一種粒子群優(yōu)化的最小二乘支持向量機對風(fēng)速數(shù)據(jù)進行重構(gòu)，但將缺失風(fēng)速置為零后的時間序列參與模型訓(xùn)練，如果連續(xù)缺失數(shù)據(jù)較多，則參與模型訓(xùn)練的數(shù)據(jù)本身誤差較大，無法保證插值效果，且計算量大；自然風(fēng)速時間序列具有典型的分形特征［12-13］；文獻［14-15］用分形插值算法及分段分形插值算法處理TurbSim 仿真軟件在特定條件下生成的模擬風(fēng)速時間序列，這種理想風(fēng)速時間序列不受瞬時天氣條件的影響，與實測風(fēng)速有一定的差距。

本文將自適應(yīng)粒子群算法優(yōu)化的分形插值算法應(yīng)用到甘肅省某風(fēng)電場實測風(fēng)速時間序列的連續(xù)缺失值的重構(gòu)問題，實驗證明針對波動強烈、連續(xù)缺失數(shù)據(jù)量大的風(fēng)速時間序列，分形插值算法可以取得比傳統(tǒng)插值算法更高的精度。

1 分形插值原理

分形插值的概念是由美國數(shù)學(xué)家Barnsley于1986年提出來的。分形插值［16］與多項式插值和曲線插值不同，它是由一類特殊的迭代函數(shù)系統(tǒng)（Iterated Function System，IFS）產(chǎn)生的，它是針對集合而非針對點的，能夠反映出相鄰已知相關(guān)點之間的局部特征，得到比傳統(tǒng)插值方法更高的精度。

1.1 數(shù)據(jù)集

假設(shè)點集 {(xi，yi)∈R2|i=0，1，…，N}，其中{(xi，yi)∈R2|i=0，1，…，N}，對應(yīng)于這個數(shù)據(jù)集合的插值函數(shù)是一個連續(xù)函數(shù)f：[x0，xN]→R，使得：

點(xi，yi)∈R2稱為插值點，函數(shù)f內(nèi)插該數(shù)據(jù)且f的圖形穿過該插值點。吸引子G={(x，f(x))；x∈[x0，xN]}是內(nèi)插數(shù)據(jù)的連續(xù)函數(shù)的圖像。

1.2 構(gòu)造IFS

考慮IFS{R2；Wn，n=1，2，…，N}，其中Wn是具有如下形式的仿射變換：

將已知點代入方程（4），不斷迭代，使系統(tǒng)收斂到其構(gòu)造的吸引子G附近，即完成了插值運算。

1.3 分形維數(shù)

分形維數(shù)是刻畫分形特征的定量指標(biāo)，它可以衡量兩個分形的相似程度，是衡量分形物體不規(guī)則程度的參數(shù)。對于風(fēng)速時間序列，分形維數(shù)表示研究時間段內(nèi)風(fēng)速波動的強烈程度，分形維數(shù)越大表明風(fēng)速波動越劇烈［17］。分形維數(shù)的計算遵循如下定理：假定數(shù)據(jù)集{(xi，yi)}不共線且則滿足G的分形維數(shù)D是如下方程的唯一常數(shù)解：

否則，滿足G的分形維數(shù)等于1。

式（5）表明分形維數(shù)D不依賴于因變量（數(shù)據(jù)集），而自由變量垂直比例因子dn較為關(guān)鍵，它能影響分形維數(shù)D的確定。目前常用盒維數(shù)法求解分形維數(shù)，即采用邊長為δi(i=1，2，…)的正方形的網(wǎng)格覆蓋分形體，數(shù)出覆蓋分形體所需盒子的數(shù)目N(δi)，隨著 limδi→0，得到一組數(shù)據(jù)(lnδi-1，lnN(δi))(i=1，2，…，M)，用最小二乘法擬合，得到的擬合直線的斜率即為分形體的分形維數(shù)。

2 自適應(yīng)變異粒子群優(yōu)化算法

2.1 粒子群優(yōu)化算法

粒子群優(yōu)化（Particle Swarm Optimization，PSO）算法是一種群體智能的優(yōu)化算法，源于對鳥類捕食行為的研究。這種算法最早由Kennedy 和Eberhart 在1995 年提出，是一種非常成熟的求解優(yōu)化問題的方法。該算法及其改進算法已在各種行業(yè)得到廣泛應(yīng)用［18-20］，并取得了很好的效果。粒子群算法收斂速度快，具有很強的通用性。算法中每個粒子都代表問題的一個潛在最優(yōu)解，用位置、速度和適應(yīng)度值三項指標(biāo)表示粒子的特征。位置用來存儲用來優(yōu)化的變量，在本文中表示分形插值的垂直比例因子參數(shù)值，適應(yīng)度值表示粒子的優(yōu)劣，粒子的速度決定粒子移動的方向和距離，速度隨著自身及其他粒子的移動經(jīng)驗進行動態(tài)調(diào)整。最初隨機初始化粒子，計算每個粒子的適應(yīng)度值，根據(jù)適應(yīng)度值更新個體極值和群體極值，然后粒子在解空間中運動，通過跟蹤個體極值和群體極值更新個體的位置和速度，粒子每更新一次位置，就重新計算一次適應(yīng)度值，通過比較新粒子的適應(yīng)度值和個體極值、群體極值的適應(yīng)度值更新個體極值和群體極值，最終的群體極值即最優(yōu)解。

每次迭代過程中，粒子通過個體極值和群體極值更新自身的位置和速度，公式如下：

其中：Vi為粒子的速度，Xi為粒子的位置，Pi為個體極值，Pg為群體極值，ω為慣性權(quán)重，c1和c2為加速度因子，是非負的常數(shù)，k為當(dāng)前迭代次數(shù)，r1和r2為分布在［0，1］區(qū)間的隨機數(shù)，γ為約束因子。

其中ω的計算公式如下：

其中：ωstart和ωend分別為初始慣性權(quán)重和迭代至最大次數(shù)時的慣性權(quán)重，分別取值ωstart=0.9、ωend=0.4，Tmax為最大迭代次數(shù)。

本文中參數(shù)取值為c1=c2=1.49，γ=0.5，最大迭代次數(shù)Tmax=100，種群大小p=15。

2.2 自適應(yīng)變異粒子群算法

粒子群算法雖然具有收斂快、通用性強的優(yōu)點，但同時存在易早熟收斂、搜索精度低等缺點，因此，在本文中提出將遺傳算法中的變異思想引入PSO算法，即在粒子每次更新之后，引入變異算子，以一定的概率重新初始化粒子。變異操作拓展了在迭代中不斷縮小的搜索空間，使粒子能夠跳出當(dāng)前搜索到的最佳位置，在更大的空間中開展搜索，保持了粒子的多樣性，同時提高算法搜索到更優(yōu)值的可能性。本文中采用的變異算子公式如下：

其中r為在［0，1］區(qū)間的隨機數(shù)。

3 衡量分形插值結(jié)果的指標(biāo)

3.1 插值誤差

插值曲線與實際曲線之間的誤差稱為插值誤差。本次實驗采用均方根誤差（Root Mean Square Error，RMSE）對插值結(jié)果進行評價：

3.2 納什效率系數(shù)

納什效率系數(shù)（Nash-Sutcliffe Efficiency coefficient，NSE），用于表示模擬結(jié)果與實際結(jié)果的吻合程度，本文用于驗證插值曲線與實際曲線的吻合程度，吻合程度越高表示插值效果越好，公式如下：

4 實驗結(jié)果與討論

所有實驗運行環(huán)境均為：1.60 GHz CPU（8 CPUs），8.00 GB RAM，Windows 10，Matlab R2019a。

4.1 實驗數(shù)據(jù)統(tǒng)計分析

采用甘肅省酒泉風(fēng)電場2017 年1 月9 日09：55 至22：20的間隔5 min 的150條數(shù)據(jù)集和2017 年1 月10 日00：00 至14：30 間隔5 min 的175 條數(shù)據(jù)集進行實驗，簡稱Dataset A 和Dataset B，這兩組數(shù)據(jù)統(tǒng)計指標(biāo)即：均值（Mean），最小值（Min），最大值（Max）和標(biāo)準(zhǔn)方差（Standard deviation，Std.），如表1和圖1所示，由此可知該組數(shù)據(jù)表現(xiàn)出非線性和較強的波動性。在兩組數(shù)據(jù)集中分別抽取30 個數(shù)據(jù)點和40 個數(shù)據(jù)點作為插值初始點。

表1 實驗數(shù)據(jù)集的統(tǒng)計指標(biāo)Tab.1 Statistical indexes of experimental datasets

圖1 實驗數(shù)據(jù)集的曲線Fig.1 Curves of experimental datasets

4.2 參數(shù)的設(shè)定

4.2.1 垂直比例因子

由式（4）可知，有一個自由變量dn(|dn|＜1)，即垂直比例因子，它是影響分形插值精度的關(guān)鍵因素，當(dāng)dn=0 時，分形插值函數(shù)就變成分段線性插值函數(shù)。圖2 為選取不同的垂直比例因子時圖形呈現(xiàn)的不同形態(tài)，由圖2可以看出隨著|dn|的增大，分形插值擬合曲線越來越粗糙，波動幅度越來越大，曲線的局部特征也更加明顯。分形維數(shù)是刻畫風(fēng)速序列波動的強烈程度的指標(biāo)，圖3 驗證了|dn|越大，分形維數(shù)D越大，即風(fēng)速曲線波動越強烈。

圖2 垂直比例因子對曲線的影響Fig.2 Influence of vertical scaling factor on curve

圖3 分形維數(shù)隨垂直比例因子的變化Fig.3 Change of fractal dimension with vertical scaling factor

但|dn|并不是越大越好，即不代表插值曲線局部特征越復(fù)雜越好，圖4 為RMSE 隨dn取值的變化，由此可知，dn嚴重影響分形插值的誤差。

圖4 RMSE隨垂直比例因子的變化Fig.4 Change of RMSE with vertical scaling factor

本文采用自適應(yīng)變異粒子群算法對垂直比例因子dn進行尋優(yōu)，粒子群中每個粒子包含垂直比例因子，且以式（9）RMSE為適應(yīng)度值，即最優(yōu)垂直比例因子取值滿足以下公式：

4.2.2 迭代次數(shù)

表2 給出了迭代次數(shù)對分形維數(shù)、插值誤差、運行時間及插值點數(shù)的影響。由表2 可知，隨著迭代次數(shù)的增多，分形維數(shù)基本不變，即曲線波動的強烈程度不隨迭代次數(shù)的改變而改變，插值誤差也基本不受迭代次數(shù)的影響，但迭代次數(shù)越多，插值得到的數(shù)據(jù)點數(shù)越大，曲線細節(jié)會更豐富，但不會變光滑。圖5 給出了迭代次數(shù)分別為2 和4 時的圖形，由局部放大圖可以明顯看到迭代2 次和迭代4 次插值得到的數(shù)據(jù)點的密集程度的差別。迭代次數(shù)n與生成的插值點數(shù)m之間的關(guān)系［17］為：

其中q為初始插值點的個數(shù)。雖然迭代次數(shù)越多，插值得到的曲線細節(jié)越豐富，但該參數(shù)不宜取得過大，符合實際插值要求即可，迭代次數(shù)過大會增加計算時間，降低計算效率，本次實驗迭代次數(shù)選擇2次即可滿足要求。

圖5 迭代次數(shù)對插值點數(shù)的影響Fig.5 Influence of iteration number on interpolation point number

表2 迭代次數(shù)對各指標(biāo)的影響Tab.2 Influence of iteration number on different indexes

4.3 分形插值與傳統(tǒng)插值的插值結(jié)果的對比研究

分別采用三次樣條插值、Lagrange 插值等傳統(tǒng)插值方法和分形插值算法進行比較分析，Dataset A 數(shù)據(jù)集的插值結(jié)果如圖6 所示。結(jié)果表明：分形插值較好地保持了圖形的形狀以及粗糙程度，能更好地凸顯圖形的局部細節(jié)特征；而Lagrange 插值和三次樣條插值只是保持了圖形的大致形狀，相鄰點間用光滑曲線連接，忽略了圖形局部豐富的信息。

表3是對Dataset A數(shù)據(jù)集插值結(jié)果的各種量化指標(biāo)，三種插值結(jié)果的均值與原始數(shù)據(jù)相差均不大，對于表征曲線波動強烈程度的分形維數(shù)，Lagrange插值與分形插值相當(dāng)，三次樣條插值最差；分形插值誤差最小為0.501 8，與Lagrange插值算法比較，減小了66.52%，與三次樣條插值算法比較，減小了58.57%；分形插值序列與原始序列的吻合程度同樣最佳，達到0.812 0，遠遠超過Lagrange插值和三次樣條插值兩種傳統(tǒng)插值方法。

圖6 Dataset A插值效果比較Fig.6 Interpolation effect comparison of Dataset A

由以上分析可知，對于波動強烈并且連續(xù)缺失數(shù)據(jù)較多的風(fēng)速時間序列，分形插值不僅可以保持圖形的整體形狀和波動特性，而且插值精度更高，與原始曲線吻合更好，與傳統(tǒng)的插值算法相比，分形插值具有明顯的優(yōu)勢。為了進一步證實上述結(jié)論的可靠性，對Dataset B數(shù)據(jù)集插值計算，得到插值結(jié)果如圖7所示。

表3 Dataset A 插值性能對比Tab.3 Interpolation performance comparison of Dataset A

圖7 Dataset B插值效果比較Fig.7 Interpolation effect comparison of Dataset B

表4 給出了Dataset B 數(shù)據(jù)集的各插值結(jié)果。Dataset B 數(shù)據(jù)集的各插值誤差RMSE 和納什效率系數(shù)NSE 的規(guī)律與Dataset A 的規(guī)律相同，進一步說明了分形插值比傳統(tǒng)的插值方法更適合連續(xù)缺失數(shù)據(jù)較多且波動強烈的風(fēng)速序列。

表4 Dataset B 插值性能對比Tab.4 Interpolation performance comparison of Dataset B

5 結(jié)語

自適應(yīng)變異粒子群優(yōu)化的分形插值算法為風(fēng)速序列的連續(xù)缺失數(shù)據(jù)插值提供了一種新途徑，并對分形插值算法中的兩個重要參數(shù)進行討論分析，利用自適應(yīng)變異粒子群算法對垂直比例因子尋優(yōu)，最后對分形插值算法與傳統(tǒng)插值算法的進行了對比研究。通過計算和驗證得出如下結(jié)論：

1）垂直比例因子是影響分形插值精度的重要因素，采用自適應(yīng)變異粒子群優(yōu)化垂直比例因子，分形插值取得了很好的插值效果，證明了其有效性。

2）針對連續(xù)缺失的波動強烈的風(fēng)速序列，分形插值比傳統(tǒng)插值方法更有優(yōu)勢，分形插值不僅可以保持風(fēng)速曲線的整體變化趨勢，而且可以反映數(shù)據(jù)的局部波動特征，并且具有更高的插值精度和吻合度。