国产日韩欧美一区二区三区三州_亚洲少妇熟女av_久久久久亚洲av国产精品_波多野结衣网站一区二区_亚洲欧美色片在线91_国产亚洲精品精品国产优播av_日本一区二区三区波多野结衣 _久久国产av不卡

?

基于VMD和IBA-LSSVM的短期風(fēng)電功率預(yù)測(cè)

2021-11-24 10:35:56陳澤坤
關(guān)鍵詞:電功率蝙蝠分量

王 瑞,陳澤坤,逯 靜

(1.河南理工大學(xué)電氣工程與自動(dòng)化學(xué)院, 河南 焦作 454000; 2.河南理工大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院, 河南 焦作 454000)

近年來為緩解能源緊張,減少環(huán)境污染,新能源的開發(fā)和利用持續(xù)增加,其中風(fēng)能是具有巨大利用價(jià)值的新型清潔環(huán)保能源,被廣泛應(yīng)用于風(fēng)力發(fā)電。未來的電力系統(tǒng)必然將為可持續(xù)的全球經(jīng)濟(jì)增長(zhǎng)提供更高滲透率的清潔能源,然而大量清潔能源的不斷接入給電力系統(tǒng)提出了前所未有的挑戰(zhàn)[1]。對(duì)風(fēng)電場(chǎng)的功率進(jìn)行短期預(yù)測(cè),可以使電力調(diào)度部門能夠提前根據(jù)風(fēng)電功率變化,及時(shí)調(diào)整調(diào)度計(jì)劃,保證電能質(zhì)量,降低電力系統(tǒng)運(yùn)行成本,這是減輕風(fēng)電對(duì)電網(wǎng)造成不利影響、提高電網(wǎng)中風(fēng)電裝機(jī)比例的一種有效途徑[2]。

目前,國(guó)內(nèi)外對(duì)于風(fēng)力發(fā)電功率的預(yù)測(cè)已有了深入研究,按照預(yù)測(cè)模型的不同,可分為物理方法、統(tǒng)計(jì)方法和學(xué)習(xí)方法。其中,物理方法需要對(duì)所在風(fēng)電場(chǎng)進(jìn)行建模[3],由于受氣象預(yù)報(bào)更新頻率的影響,該方法更適合中期風(fēng)電功率預(yù)測(cè)。統(tǒng)計(jì)方法包括回歸分析法[4]、指數(shù)平滑法、時(shí)間序列法[5-6]和灰色預(yù)測(cè)法。這種預(yù)測(cè)模型計(jì)算簡(jiǎn)單,但隨著預(yù)測(cè)時(shí)間的增加預(yù)測(cè)精度會(huì)快速下降,且不能很好地適應(yīng)非線性影響因素。

學(xué)習(xí)方法包括人工神經(jīng)網(wǎng)絡(luò)法[7-9]、決策樹[10-11]和支持向量機(jī)[12-14]等。其中,最小二乘支持向量機(jī)(least squares support vector machine,LSSVM)具有預(yù)測(cè)精度高、計(jì)算簡(jiǎn)單等優(yōu)點(diǎn)。但是在實(shí)際的風(fēng)電功率預(yù)測(cè)中,單一的預(yù)測(cè)模型存在局限性,無法取得最佳的預(yù)測(cè)效果,因此目前多采用組合預(yù)測(cè)模型,如周松林[15]引入粒子群算法(particle swarm optimization,PSO)優(yōu)化LSSVM參數(shù)尋優(yōu),可有效縮短搜索時(shí)間,但存在過早收斂的問題。趙鳳展等[16]采用蝙蝠算法(bat algorithm, BA)優(yōu)化LSSVM模型,與PSO-LSSVM相比搜索過程具有更好的收斂性,但存在不能保持優(yōu)化能力等問題。此外,組合預(yù)測(cè)模型中經(jīng)常會(huì)用信號(hào)分解的方法,目的是將原始序列分解成一系列子模態(tài)以降低非平穩(wěn)性,對(duì)分解的序列分別建立預(yù)測(cè)模型并重組來實(shí)現(xiàn)最終預(yù)測(cè),如姜貴敏等[17]通過集成經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition,EEMD)將功率歷史數(shù)據(jù)分解為一系列相對(duì)平穩(wěn)的子序列,解決了經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)方法所產(chǎn)生的模態(tài)混疊現(xiàn)象。變分模態(tài)分解(variational mode decomposition,VMD)是一種非遞歸、變模式的分解方法,克服了EEMD遞歸求解的缺點(diǎn),具有更好的諧波分離效果[18]。

本文在現(xiàn)有研究基礎(chǔ)上,提出了一種基于VMD和改進(jìn)蝙蝠算法(improved bat algorithm,IBA)優(yōu)化的LSSVM預(yù)測(cè)模型(VMD-IBA-LSSVM模型),并利用該模型對(duì)寧夏某風(fēng)電場(chǎng)的發(fā)電功率進(jìn)行預(yù)測(cè),通過與其他幾種典型模型進(jìn)行比較來驗(yàn)證模型的有效性。

1 VMD方法

VMD是基于EMD提出的一種自適應(yīng)、完全非遞歸的模態(tài)變分和信號(hào)處理的方法。與EMD相比,VMD對(duì)噪聲和采樣誤差具有更強(qiáng)的魯棒性。VMD的分解過程是變分問題的求解過程,可分為變分問題的構(gòu)造和求解過程[19]。

1.1 變分問題的構(gòu)造

步驟1對(duì)于每個(gè)模態(tài)函數(shù),采用Hilbert變換計(jì)算相關(guān)的解析信號(hào),以獲得單側(cè)頻譜。

步驟2通過混合一個(gè)調(diào)諧到各自中心頻率的指數(shù)項(xiàng),將各個(gè)模態(tài)的頻譜調(diào)制到相應(yīng)的基頻帶。

步驟3由解調(diào)信號(hào)的高斯平滑度,得到的約束變分問題:

(1)

式中:{uk}、{ωk}——子信號(hào)及其相應(yīng)的中心頻率集合;k——子信號(hào)總數(shù);t——采樣時(shí)刻;δ(t)——狄拉克分布;f(t)——一個(gè)序列。

1.2 變分問題的求解

步驟1將約束性變分問題重構(gòu)為非約束性變分問題,增廣的拉格朗日表達(dá)式為

(2)

式中:λ(t)——拉格朗日乘法算子;α——二次懲罰因子。

步驟2通過交替更新uk,n+1、ωk,n+1和λn+1求解增廣后拉格朗日表達(dá)式中的“鞍點(diǎn)”。VMD更新過程如下:

(3)

(4)

(5)

步驟3對(duì)于給定判別精度e>0,若滿足式(6),則VMD收斂,停止更新。

(6)

2 基于IBA優(yōu)化的LSSVM預(yù)測(cè)模型

2.1 LSSVM

LSSVM是一種新型支持向量機(jī)方法,LSSVM采用最小二乘線性系統(tǒng)作為損失函數(shù),代替?zhèn)鹘y(tǒng)的支持向量機(jī)采用的二次規(guī)劃方法。利用等式約束取代SVM中的不等式約束,將原問題轉(zhuǎn)化為一個(gè)解線性方程組的問題。LSSVM的優(yōu)化問題[20]可以轉(zhuǎn)化為

(7)

式中:w——權(quán)向量;γ——正則化參數(shù);ek——誤差變量;φ(xk)——xk在特征空間的映射;b——偏置。

式(7)可采用拉格朗日乘數(shù)法把原問題優(yōu)化,同時(shí)根據(jù)KKT最優(yōu)條件求解。本文選擇RBF函數(shù)作為L(zhǎng)SSVM的核函數(shù)。此外,超參數(shù)(C,σ)的選取對(duì)模型的預(yù)測(cè)結(jié)果有顯著影響,本文采用IBA對(duì)LSSVM模型參數(shù)尋優(yōu),建立最優(yōu)預(yù)測(cè)模型。

2.2 IBA

BA是一種仿生尋優(yōu)算法,模擬自然界中微蝙蝠的回聲定位行為。BA將回聲定位理想化,將蝙蝠種群初始化為一組隨機(jī)解,然后通過調(diào)節(jié)蝙蝠發(fā)出的聲波頻率更新個(gè)體的脈沖速率和脈沖響度迭代搜尋最優(yōu)解,且在最優(yōu)解周圍通過隨機(jī)飛行產(chǎn)生局部新解,加強(qiáng)了局部搜索。

蝙蝠個(gè)體更新其聲波頻率、速度與位置,公式表述為

fi=fmin+(fmax-fmin)β

(8)

vi,t=vi,t-1+(xi,t-x*)fi

(9)

xi,t=xi,t-1+vi,t

(10)

式中:fmax、fmin——最大與最小頻率值;v——個(gè)體速度;x——個(gè)體位置;β——[0,1]之間的隨機(jī)數(shù)。

在局部尋優(yōu)過程中,每只蝙蝠更新公式如下:

Xnew=Xold+εAt(ε∈[-1,1])

(11)

式中:Xold——當(dāng)前最優(yōu)解;ε——隨機(jī)數(shù);At=〈Ai,t〉——全部個(gè)體在第t次迭代時(shí)的平均脈沖響度。

當(dāng)蝙蝠發(fā)現(xiàn)目標(biāo)逼近時(shí),會(huì)更新發(fā)射脈沖的發(fā)射速率和響度:

Ai,t+1=αAi,t

(12)

ri,t-1=ri,0(1-e-γt)

(13)

式中:Ai——脈沖響度;α——脈沖音響衰減系數(shù);ri——脈沖發(fā)射速率;γ——搜索頻度增強(qiáng)系數(shù)。

根據(jù)上述BA的優(yōu)化原理可知,其參數(shù)更新方式相對(duì)固定,個(gè)體本身缺乏變異機(jī)制,致使存在后期收斂慢、收斂精度低以及容易陷入局部極小值等問題。為了克服BA的缺點(diǎn),本文采用改進(jìn)的慣性權(quán)重、自適應(yīng)頻率與變異機(jī)制來改善尋優(yōu)過程[21],形成IBA。

公式(9)中,速度更新時(shí)添加慣性權(quán)重以改進(jìn)速度更新的方向,使得種群中個(gè)體可以有效地跳出局部最優(yōu)點(diǎn):

(14)

式中:g——慣性權(quán)重因子;fit——適應(yīng)度函數(shù);N——種群數(shù)量。

當(dāng)蝙蝠i的適應(yīng)度低于平均適應(yīng)度時(shí),該蝙蝠將賦予較低的權(quán)重,增強(qiáng)其尋優(yōu)全局更優(yōu)解的能力;而蝙蝠i的適應(yīng)度高于平均適應(yīng)度時(shí),該蝙蝠將賦予較高的權(quán)重,增加其跳出局部最優(yōu)解的機(jī)會(huì),間接增大該蝙蝠尋優(yōu)全局最優(yōu)解的能力,主要模仿粒子群算法的慣性權(quán)重策略。

在BA中,存在一些蝙蝠已經(jīng)處在最優(yōu)解邊緣的,仍采用與平時(shí)一致的頻率來尋優(yōu),最終會(huì)影響其尋找最優(yōu)解的機(jī)會(huì),因此采用自適應(yīng)頻率調(diào)整:

(15)

range=di,max-di,min

(16)

(17)

式中:di——第i個(gè)解到最優(yōu)解距離;range——最大距離與最小距離的差值。

由此,速度更新公式變?yōu)?/p>

vij,t=vij,t-1+(xij,t-xj,*)fj

(18)

當(dāng)蝙蝠都趨于收斂時(shí),部分種群陷入了局部最優(yōu)解,此時(shí)增加變異機(jī)制,跳出該局部最優(yōu)解。本文采用的方法是產(chǎn)生一個(gè)隨機(jī)數(shù),當(dāng)這個(gè)隨機(jī)數(shù)大于變異概率的時(shí)候,對(duì)蝙蝠重新初始化。

2.3 改進(jìn)的IBA-LSSVM預(yù)測(cè)模型的建立

為了進(jìn)一步提高LSSVM模型的預(yù)測(cè)精度和速度,需要對(duì)模型的參數(shù)進(jìn)行優(yōu)化,優(yōu)化步驟如下:

步驟1設(shè)定LSSVM模型中懲罰參數(shù)C、核參數(shù)σ的取值范圍。

步驟2種群基本參數(shù)化,設(shè)定種群數(shù)個(gè)體xi(i=1,2,…,N)、脈沖頻率最大值fmax和最小值fmin、脈沖響度Ai、脈沖發(fā)射率ri、空間維度d和最大迭代次數(shù)MI。

步驟3初始化種群中每只蝙蝠個(gè)體的位置xi和速度vi,其中蝙蝠i的位置代表著參數(shù)C和σ。

步驟4計(jì)算每只蝙蝠的適應(yīng)度值,尋找當(dāng)前時(shí)刻最優(yōu)解。以蝙蝠位置對(duì)應(yīng)的參數(shù)訓(xùn)練LSSVM模型,然后選取訓(xùn)練集進(jìn)行訓(xùn)練和測(cè)試,按照降序排列,找到當(dāng)前最優(yōu)解。

步驟5分別按照式(8)(9)和(18)更新種群中各個(gè)蝙蝠個(gè)體的脈沖發(fā)射率、所在位置和飛行速度,采用改進(jìn)的慣性權(quán)重、自適應(yīng)頻率與變異機(jī)制來改善尋優(yōu)過程。

步驟6出現(xiàn)第一個(gè)隨機(jī)數(shù)rand1,當(dāng)rand1>ri時(shí),根據(jù)式(11)更新出局部最優(yōu)解Xnew。

步驟7出現(xiàn)第二個(gè)隨機(jī)數(shù)rand2,若rand2

步驟8對(duì)蝙蝠的適應(yīng)度重新排序,確認(rèn)當(dāng)前最優(yōu)值,重復(fù)迭代過程,直到滿足設(shè)定的終止條件,停止循環(huán)并輸出全局最優(yōu)解。

3 VMD-IBA-LSSVM模型

VMD-IBA-LSSVM模型的建模流程為:①利用VMD將具有非線性、隨機(jī)性的原始風(fēng)電功率序列分解為一系列平穩(wěn)的模態(tài)分量,根據(jù)各個(gè)子模態(tài)的近似熵的分析結(jié)果進(jìn)行子序列重組,對(duì)于每個(gè)重組后模態(tài)分量,結(jié)合歷史氣象數(shù)據(jù),分別建立第2.3節(jié)中的IBA-LSSVM模型; ②疊加各子模態(tài)模型預(yù)測(cè)值,得到最終的風(fēng)電功率預(yù)測(cè)值。建模過程如圖1所示。

圖1 風(fēng)電功率預(yù)測(cè)流程Fig.1 Flow chart of wind power prediction

輸入變量的選取對(duì)預(yù)測(cè)精度有直接影響,本文將輸入原始數(shù)據(jù)集分為歷史氣象數(shù)據(jù)和歷史功率數(shù)據(jù)兩類。其中歷史氣象數(shù)據(jù)包含多種風(fēng)速等,歷史功率數(shù)據(jù)作為輸入變量的主體,需要進(jìn)行VMD分解,而風(fēng)速數(shù)據(jù)只需要加入分解后的子序列中。原始數(shù)據(jù)集中各個(gè)數(shù)據(jù)量綱不同,為了提高風(fēng)電功率預(yù)測(cè)精度,需要進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化處理,將原始數(shù)據(jù)線性化轉(zhuǎn)換到[0,1]的范圍,歸一化公式如下:

(19)

4 算例分析

4.1 算例概況

a.數(shù)據(jù)來源。以寧夏某風(fēng)電場(chǎng)2017年1月1—31日的樣本作為研究數(shù)據(jù),數(shù)據(jù)包括24 h風(fēng)電功率數(shù)據(jù)和日環(huán)境數(shù)據(jù)。數(shù)據(jù)采樣時(shí)間間隔為15 min,總共96個(gè)點(diǎn)。在實(shí)際建模過程中,將2017年1月1—28日數(shù)據(jù)作為L(zhǎng)SSVM模型的訓(xùn)練集,其余數(shù)據(jù)作為模型的測(cè)試集。

b.預(yù)測(cè)模型。本文VMD-IBA-LSSVM模型和BP、SVM、LSSVM、VMD-BP、VMD-SVM、VMD-LSSVM、IBA-LSSVM 7種對(duì)比模型。8種模型在輸入量不變的情況下進(jìn)行對(duì)比預(yù)測(cè)。

c.功率預(yù)測(cè)模型參數(shù)。在建立IBA-LSSVM模型過程中,懲罰參數(shù)C、核參數(shù)σ的尋優(yōu)范圍均為[0.001,1 000],IBA參數(shù)設(shè)置為:種群數(shù)大小N=10,脈沖頻率最大值fmax=2,最小值fmin=0,解的維度d=2,脈沖響度Ai=0.9,初始脈沖發(fā)射率r0=0.5,最大迭代次數(shù)MI=1 000,慣性權(quán)重因子最大值wmax=0.5、最小值wmin=0.2。對(duì)比模型BP隱含層神經(jīng)元數(shù)量為20,學(xué)習(xí)率為0.01,學(xué)習(xí)目標(biāo)為0.001,迭代次數(shù)為5 000。對(duì)比模型SVM中的模型學(xué)習(xí)參數(shù)C與ε與LSSVM中參數(shù)選取均通過網(wǎng)格搜索法優(yōu)化得出,參數(shù)范圍為[2-8,28],迭代步長(zhǎng)為0.1。

4.2 評(píng)價(jià)標(biāo)準(zhǔn)

采用均方根誤差(RMSE)和平均絕對(duì)誤差(MAE)作為模型預(yù)測(cè)精度評(píng)價(jià)指標(biāo)。

4.3 原始風(fēng)電功率序列的分解結(jié)果

VMD的參數(shù)設(shè)置如下:懲罰參數(shù)C=2 000;模態(tài)函數(shù)個(gè)數(shù)K由于當(dāng)其大于5后子序列趨于相似,因此取K=5;初始中心頻率ω=0;收斂判據(jù)t=10-7。VMD分解效果如圖2和圖3所示。

圖2 原始風(fēng)電功率序列及VMD分解結(jié)果Fig.2 Original wind power sequence and VMD decomposition results

圖3 原始風(fēng)電功率序列及模態(tài)函數(shù)重構(gòu)結(jié)果Fig.3 Original wind power sequence and modal function reconstruction results

從圖2可以看出,原始風(fēng)電功率被分解為多個(gè)存在不同波動(dòng)性的子序列,若分別對(duì)各個(gè)子序列分別建立模型進(jìn)行預(yù)測(cè),不僅增加了任務(wù)量,而且忽略了各個(gè)子序列的相關(guān)性。本文采用近似熵度量各個(gè)子序列的復(fù)雜度,將具有相關(guān)性的序列進(jìn)行重組,合并成新的序列,形成趨勢(shì)分量、細(xì)節(jié)分量、隨機(jī)分量,這樣不僅可以有效縮短運(yùn)算時(shí)間,而且可以更好地突顯同類序列的特性[22]。

根據(jù)各分量的近似熵,圖2(b)中第一分量為0.01數(shù)量級(jí),圖2(c)中第2分量和圖2(d)中第3分量為0.1數(shù)量級(jí),圖2(e)中第4分量和圖2(f)中第5分量為1數(shù)量級(jí),按照數(shù)量級(jí)不同可以將子序列重新分組。即第1分量作為趨勢(shì)分量,第2和第3分量組成細(xì)節(jié)分量,第4和第5分量組成隨機(jī)分量,重構(gòu)的新序列如圖3所示。

分析圖3重構(gòu)后的分量序列,3種分量都具有各自明顯的特點(diǎn),其中趨勢(shì)分量波動(dòng)平緩,可以將原始序列的總體波動(dòng)趨勢(shì)表現(xiàn)出來;細(xì)節(jié)分量規(guī)律性較強(qiáng),可以很好地表征出原始序列的細(xì)節(jié)波動(dòng);隨機(jī)分量以很好地表征出原始序列的細(xì)節(jié)波動(dòng),隨機(jī)分量具有一定的隨機(jī)性和波動(dòng)性表明了一些不確定因素造成的波動(dòng)。

4.4 預(yù)測(cè)結(jié)果對(duì)比分析

采用上述8種預(yù)測(cè)模型對(duì)2017年1月29—31日進(jìn)行提前24 h風(fēng)電功率預(yù)測(cè),最后一天的預(yù)測(cè)結(jié)果如圖4所示,可以看出VMD和IBA-LSSVM模型相較于其他對(duì)比模型可以更好地對(duì)真實(shí)曲線的波動(dòng)規(guī)律進(jìn)行預(yù)測(cè),整體趨勢(shì)與真實(shí)曲線更為貼合。而單一的預(yù)測(cè)模型在風(fēng)電功率預(yù)測(cè)中其曲線明顯滯后于實(shí)際曲線,采用VMD對(duì)原始序列進(jìn)行分解,分別對(duì)各個(gè)分量建立預(yù)測(cè)模型可以改善這一現(xiàn)象,預(yù)測(cè)精度有不同程度的提高。

圖4 不同模型1月31日功率預(yù)測(cè)曲線Fig.4 Power prediction curves of different models on 31st January

以上8種預(yù)測(cè)模型2017年1月29—31日預(yù)測(cè)誤差及3日預(yù)測(cè)誤差平均值如表1所示。

表1 不同模型預(yù)測(cè)誤差對(duì)比

從表1預(yù)測(cè)誤差值可以看出,VMD-IBA-LSSVM模型的RMSE和MAE均低于其他模型,證明該模型預(yù)測(cè)效果相比于其他模型更好。同時(shí)相較于單一的預(yù)測(cè)模型BP、SVM和LSSVM,經(jīng)過VMD進(jìn)行數(shù)據(jù)分解,建立組合預(yù)測(cè)模型的VMD-BP、VMD-SVM和VMD-LSSVM其預(yù)測(cè)性能都有不同程度的提高,其中RMSE精度分別提高了4.95%、10.99%和11.34%;MAE精度分別提高了8.96%、9.68%和8.78%。在運(yùn)行時(shí)間方面,基礎(chǔ)模型BP、SVM和LSSVM平均訓(xùn)練用時(shí)分別為9 s、384 s和127 s,雖然LSSVM耗時(shí)較長(zhǎng),但是預(yù)測(cè)效果最好。此外,采用IBA對(duì)LSSVM模型優(yōu)化后,其訓(xùn)練時(shí)間縮短為48 s。分析對(duì)比LSSVM與IBA-LSSVM預(yù)測(cè)模型,在采用了IBA后,其預(yù)測(cè)精度具有明顯的提高,在此基礎(chǔ)上利用IBA-LSSVM預(yù)測(cè)模型對(duì)VMD技術(shù)分解的原始序列產(chǎn)生的各個(gè)變量分別進(jìn)行建模,從而進(jìn)一步提高了預(yù)測(cè)精度,RMSE和MAE分別僅為13.29%和8.80%。

5 結(jié) 論

a.針對(duì)具有非線性的原始風(fēng)電功率序列,采用VMD分解為一系列平穩(wěn)的模態(tài)分量,利用近似熵度量各個(gè)子序列的復(fù)雜度,將具有相關(guān)性的序列進(jìn)行重組,降低了預(yù)測(cè)模型的計(jì)算規(guī)模。

b.針對(duì)LSSVM模型參數(shù)難以選取,影響模型預(yù)測(cè)性能,提出了采用IBA優(yōu)化LSSVM的核參數(shù)選取,利用其全局尋優(yōu)能力強(qiáng)、收斂精度高等優(yōu)點(diǎn),為重組后的子序列分別構(gòu)建IBA-LSSVM預(yù)測(cè)模型,提高了模型預(yù)測(cè)精度。

c.VMD-IBA-LSSVM模型能有效對(duì)風(fēng)電場(chǎng)的功率進(jìn)行短期預(yù)測(cè),從而可以及時(shí)調(diào)整調(diào)度計(jì)劃,降低電力系統(tǒng)運(yùn)行成本,具有一定的應(yīng)用價(jià)值。下一步應(yīng)對(duì)風(fēng)電機(jī)組特性進(jìn)行深入分析并考慮更多的環(huán)境信息作為輸入樣本,從而進(jìn)一步提高其預(yù)測(cè)精度。

猜你喜歡
電功率蝙蝠分量
基于PCC-CNN-GRU的短期風(fēng)電功率預(yù)測(cè)
帽子的分量
輕松上手電功率
你會(huì)計(jì)算電功率嗎
一物千斤
智族GQ(2019年9期)2019-10-28 08:16:21
解讀電功率
論《哈姆雷特》中良心的分量
分量
蝙蝠
蝙蝠女
安乡县| 蕲春县| 读书| 赣州市| 郓城县| 临清市| 广元市| 海丰县| 克什克腾旗| 本溪市| 大港区| 治多县| 嘉黎县| 南平市| 青神县| 太仆寺旗| 巫溪县| 东乌珠穆沁旗| 天柱县| 方山县| 闵行区| 杨浦区| 珲春市| 张家界市| 银川市| 会东县| 沅陵县| 龙州县| 玛曲县| 凤冈县| 淮阳县| 敖汉旗| 梁河县| 舟山市| 彭水| 凤冈县| 灌南县| 兴仁县| 深水埗区| 新安县| 桂东县|