• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    實(shí)測(cè)鐵磁諧振時(shí)間序列的非線性動(dòng)力學(xué)分析

    2016-04-07 10:36:00黃艷玲司馬文霞
    電工技術(shù)學(xué)報(bào) 2016年5期
    關(guān)鍵詞:龐加萊相空間鐵磁

    黃艷玲 司馬文霞 楊 鳴 楊 慶 袁 濤

    (輸配電裝備及系統(tǒng)安全與新技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué)) 重慶 400044)

    ?

    實(shí)測(cè)鐵磁諧振時(shí)間序列的非線性動(dòng)力學(xué)分析

    黃艷玲司馬文霞楊鳴楊慶袁濤

    (輸配電裝備及系統(tǒng)安全與新技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué))重慶400044)

    摘要在僅有實(shí)測(cè)單一電壓時(shí)間序列的情況下,為準(zhǔn)確識(shí)別出電力系統(tǒng)鐵磁諧振類型,運(yùn)用坐標(biāo)延遲方法重構(gòu)與原系統(tǒng)拓?fù)湟饬x上等價(jià)的相空間。采用相平面、龐加萊截面和關(guān)聯(lián)維數(shù)3種基于重構(gòu)相空間的非線性動(dòng)力學(xué)分析工具表達(dá)電壓時(shí)間序列的運(yùn)動(dòng)特征,由此識(shí)別其所屬的鐵磁諧振類型。對(duì)3例典型的中性點(diǎn)不接地系統(tǒng)電磁式電壓互感器鐵磁諧振實(shí)測(cè)時(shí)間序列進(jìn)行了動(dòng)力學(xué)特征分析。其中兩例時(shí)間序列的相平面軌跡、龐加萊截面均呈現(xiàn)周期運(yùn)動(dòng)特性,關(guān)聯(lián)維數(shù)估計(jì)值分別為1.015 0±0.003 9、1.006 1±0.000 6,由此判別它們作周期運(yùn)動(dòng);另一例時(shí)間序列的相平面軌跡、龐加萊截面均未呈現(xiàn)周期或準(zhǔn)周期運(yùn)動(dòng)特性,關(guān)聯(lián)維數(shù)估計(jì)值為2.300 2±0.061 2,由此判別其作混沌運(yùn)動(dòng)。實(shí)測(cè)電壓時(shí)間序列受到多種因素的影響,綜合以上3種分析方法所表達(dá)的特征,對(duì)所發(fā)生的鐵磁諧振類型作出的判斷更準(zhǔn)確和更有說(shuō)服力。

    關(guān)鍵詞:鐵磁諧振時(shí)間序列非線動(dòng)力學(xué)相空間重構(gòu)關(guān)聯(lián)維數(shù)電磁式電壓互感器

    Nonlinear Dynamic Analysis of the Measured Ferroresonance Time Series

    HuangYanlingSimaWenxiaYangMingYangQingYuanTao

    (State Key Laboratory of Power Transmission Equipment and System Safety and New Technology Chongqing UniversityChongqing400044China)

    AbstractWith only a single voltage time series recorded in fact,for accurately identifying the ferroresonance types,the delay coordinates method was applied to obtain the reconstructed phase space that is topologically equivalent to the phase space of the original system.Nonlinear dynamics methods based on the reconstructed phase space,e.g.phase plane,Poincaré section and correlation dimension,were used to express the dynamic characteristics of the time series and identify the type of ferroresonance based on them.The nonlinear dynamic characteristics analysis was conducted to 3 typical cases of measured time series of ferroresonance occurred in inductive voltage transformer in isolated neutral system.The phase plane trajectories and Poincaré sections of two time series cases of them show periodic motion characteristics,their correlation dimension estimation values respectively are 1.015 0±0.003 9,1.006 1±0.000 6;therefore their motion modes were identified as fundamental mode and subharmonic mode.The phase plane trajectory and Poincaré section of another time series case do not show periodic or quasiperiodic motion characteristics,its correlation dimension estimation value is 2.300 2±0.061 2;therefore its motion mode was identified as chaotic mode.The measured voltage time series is affected by many factors,the identification result of ferroresonance mode occurred is more accurate and convictive through synthesizing the features characterized by the above 3 analysis methods.

    Keywords:Ferroresonance,time series,nonlinear dynamics,phase space construction,correlation dimension,inductive voltage transformer

    0引言

    鐵磁諧振是一種復(fù)雜的電現(xiàn)象,它可能發(fā)生在由電力變壓器或電磁式電壓互感器(PT)及帶鐵心的電抗器等非線性電感元件與電容元件構(gòu)成的回路中[1]。非線性電感元件的飽和勵(lì)磁特性使鐵磁諧振行為呈現(xiàn)非線性動(dòng)力學(xué)特性,其穩(wěn)態(tài)響應(yīng)模式高度依賴于系統(tǒng)參數(shù)和初始條件,初始條件稍有變化就可能使響應(yīng)模式發(fā)生改變[2]。穩(wěn)態(tài)響應(yīng)模式可能是周期、準(zhǔn)周期或混沌模式。在合適的條件組合下,很多設(shè)備操作行為及單相接地故障消除瞬間[3]均可能激發(fā)鐵磁諧振,由此產(chǎn)生的過(guò)電壓或過(guò)電流嚴(yán)重危害著變電站設(shè)備的絕緣,甚至造成變壓器或PT的燒毀,繼而引發(fā)大面積的停電事故。因此,當(dāng)鐵磁諧振發(fā)生時(shí),很有必要識(shí)別出諧振類型,為采取有效的過(guò)電壓抑制措施和電網(wǎng)事故分析提供指導(dǎo)。這對(duì)保障電網(wǎng)的安全運(yùn)行具有重要意義。

    常規(guī)的線性數(shù)學(xué)方法不適于分析非線性的鐵磁諧振行為,Z.Emin和其他學(xué)者將基于相空間的非線性動(dòng)力學(xué)分析工具相平面[1,2,4,5,10]、龐加萊截面[1,2,4-10]和Lyapunov指數(shù)[4,8,10]等運(yùn)用到鐵磁諧振行為的特征分析和識(shí)別中,再結(jié)合功率譜密度分析[1,6,9]對(duì)識(shí)別結(jié)果進(jìn)行驗(yàn)證。在此基礎(chǔ)上,為進(jìn)一步獲取系統(tǒng)發(fā)生諧振的參數(shù)取值范圍以及諧振類型,基于分岔理論的分岔圖[2,5,8,9]被用于預(yù)測(cè)系統(tǒng)隨某一參量變化可能出現(xiàn)的諧振類型。

    盡管基于相空間的非線性動(dòng)力學(xué)分析工具對(duì)鐵磁諧振行為的特征分析和識(shí)別取得了很好的效果,但當(dāng)前所研究的鐵磁諧振數(shù)據(jù)主要是通過(guò)建立諧振回路的微分方程組進(jìn)行數(shù)值求解[1,4,6-10]或是應(yīng)用ATP-EMTP軟件搭建回路仿真[5]所得,與實(shí)際情況有一定的差距。因建立的數(shù)學(xué)仿真模型是系統(tǒng)結(jié)構(gòu)的簡(jiǎn)化表達(dá),鐵心損耗及系統(tǒng)對(duì)地電容等系統(tǒng)參數(shù)難以準(zhǔn)確獲取,尤其對(duì)電感元件的非線性特性模擬的精度有待提高。只有對(duì)正在發(fā)生的鐵磁諧振類型進(jìn)行識(shí)別,才能為制定合適的過(guò)電壓抑制策略提供幫助,但實(shí)際觀測(cè)和獲取的常常只是鐵磁諧振系統(tǒng)單一狀態(tài)變量(如電壓的時(shí)間序列),而以上基于相空間的研究方法[1,2,4-10]至少需要兩個(gè)狀態(tài)變量的時(shí)間序列數(shù)據(jù),因此需要借助其他方法先構(gòu)造相空間。

    因一個(gè)變量隨時(shí)間的變化隱含了整個(gè)系統(tǒng)的運(yùn)動(dòng)規(guī)律,根據(jù)F.Takens提出的嵌入定理[11],從一維時(shí)間序列可重構(gòu)一個(gè)與原動(dòng)力系統(tǒng)在拓?fù)湟饬x上等價(jià)的相空間。因此,可對(duì)單一變量的時(shí)間序列運(yùn)用相空間重構(gòu)方法構(gòu)造相空間,通過(guò)分析重構(gòu)相空間軌線的分布或結(jié)構(gòu)得到原系統(tǒng)的動(dòng)力學(xué)特征。

    本文在僅有某一中性點(diǎn)不接地系統(tǒng)PT鐵磁諧振實(shí)測(cè)電壓時(shí)間序列的情況下,運(yùn)用相空間重構(gòu)方法獲得系統(tǒng)相空間。在此基礎(chǔ)上,運(yùn)用相平面、龐加萊截面和關(guān)聯(lián)維數(shù)3種基于相空間重構(gòu)的非線性動(dòng)力學(xué)分析工具表達(dá)系統(tǒng)的動(dòng)力學(xué)特征,以識(shí)別鐵磁諧振類型。

    1非線性動(dòng)力學(xué)分析方法

    1.1相空間重構(gòu)

    系統(tǒng)相空間是由狀態(tài)變量xi(i=1,2,…,n)張成的空間Rn, 相空間的一個(gè)點(diǎn)對(duì)應(yīng)系統(tǒng)的一個(gè)狀態(tài)。N.H.Parkard等[12]于1980年提出了基于單一變量時(shí)間序列的兩種相空間重構(gòu)方法:導(dǎo)數(shù)重構(gòu)法和坐標(biāo)延遲重構(gòu)法。因?qū)崪y(cè)數(shù)據(jù)混有噪聲,對(duì)其求取的微分誤差較大,因此實(shí)際應(yīng)用中主要使用坐標(biāo)延遲重構(gòu)法進(jìn)行一維時(shí)間序列的相空間重構(gòu)。

    1.1.1坐標(biāo)延遲重構(gòu)法

    設(shè)有一單變量的時(shí)間序列u={ui|i=1,2,…,N}, 對(duì)于某一給定的嵌入維數(shù)m和延遲時(shí)間τ,其相空間向量序列為

    V={Vk|uk,uk+τ,…,uk+(m-1)τ}

    (1)

    式中,k=1,2,…,n;n=N-(m-1)τ。

    由此,重構(gòu)相空間的軌線分布或結(jié)構(gòu)便反映了系統(tǒng)的運(yùn)動(dòng)特征。將重構(gòu)相空間的某兩個(gè)向量的取值序列分別作為橫、縱坐標(biāo)繪制在平面中就構(gòu)成重構(gòu)相平面。系統(tǒng)作周期運(yùn)動(dòng)時(shí),其相平面軌跡是一條不斷重復(fù)的閉合曲線[2];而作混沌運(yùn)動(dòng)時(shí),相平面軌跡永不重復(fù)[1]。

    1.1.2坐標(biāo)延遲重構(gòu)參數(shù)選取

    坐標(biāo)延遲重構(gòu)法參數(shù)即延遲時(shí)間τ和嵌入維數(shù)m的選取很關(guān)鍵。當(dāng)前對(duì)它們的選取有兩種觀點(diǎn)[13],一種觀點(diǎn)認(rèn)為二者互不相關(guān),對(duì)τ和m的選取可獨(dú)立進(jìn)行,分別有不同的確定方法。最佳τ值的求取方法主要有自相關(guān)函數(shù)法[14]和互信息法[15]。最佳m值的估計(jì)方法有計(jì)算吸引子不變量法(或G-P法)[16]、奇異值分解法[17]、偽最近鄰法[18]和偽最近鄰Cao氏改進(jìn)法[19]。另一種觀點(diǎn)認(rèn)為τ和m相互關(guān)聯(lián),即τ和m的關(guān)系與重構(gòu)相空間的嵌入窗τw(τw=(m-1)τ)密切相關(guān)。D.Kugiurmtzis[20]于1996年首次提出嵌入窗思想,基于嵌入窗思想的典型算法有C-C法[21]和自動(dòng)嵌入式相空間重構(gòu)法[22]。

    本文在確定相空間重構(gòu)參數(shù)時(shí)采用第一種觀點(diǎn),即τ和m的值采用不同的方法確定。采用互信息法獲取最佳τ值,因互信息刻畫的是兩個(gè)變量間的普遍(或總體)相關(guān)性,適用于非線性系統(tǒng)。本文采用的計(jì)算過(guò)程無(wú)需事先確定最優(yōu)m值,但通過(guò)計(jì)算一組m值下的關(guān)聯(lián)維數(shù)可間接獲取最優(yōu)m值(即G-P法確定最優(yōu)m值,實(shí)現(xiàn)方法見(jiàn)1.3.1節(jié))。

    1.2龐加萊截面

    龐加萊截面是相平面軌跡的一種簡(jiǎn)化表達(dá),是對(duì)相平面軌跡以一定的頻率采樣獲取的相點(diǎn)圖。

    若對(duì)相平面軌跡采樣的頻率為系統(tǒng)外加激勵(lì)函數(shù)的頻率(如電源頻率),則對(duì)作周期運(yùn)動(dòng)的系統(tǒng),在龐加萊截面上表現(xiàn)為一個(gè)不動(dòng)點(diǎn)或是多個(gè)孤立點(diǎn)[1];對(duì)作準(zhǔn)周期運(yùn)動(dòng)的系統(tǒng),龐加萊截面是一條分布著無(wú)限個(gè)不同點(diǎn)的封閉曲線[6,7];而當(dāng)系統(tǒng)作混沌運(yùn)動(dòng)時(shí),龐加萊截面呈現(xiàn)為局限于平面某一特定區(qū)域的一些隨機(jī)分布的相點(diǎn)集合[1]。

    1.3關(guān)聯(lián)維數(shù)

    判斷一個(gè)時(shí)間序列是否具有混沌特性最為常用的動(dòng)態(tài)不變量有Lyapunov指數(shù)[4,8,23]和吸引子維數(shù)[7,23]。由于實(shí)驗(yàn)獲取的時(shí)間序列數(shù)據(jù)ui受到多種因素(如采樣精度、采樣長(zhǎng)度、噪聲和誤差等)的影響,僅憑其最大Lyapunov指數(shù)λmax的正負(fù)就對(duì)系統(tǒng)的運(yùn)動(dòng)特性作出判斷,很可能不準(zhǔn)確[24]。自1983年P(guān).Grassberger等[25]提出針對(duì)一維時(shí)間序列數(shù)據(jù)的關(guān)聯(lián)維數(shù)計(jì)算方法(簡(jiǎn)稱G-P算法)后,關(guān)聯(lián)維數(shù)作為吸引子維數(shù)的一種度量得到了廣泛應(yīng)用。本文采用關(guān)聯(lián)維數(shù)作為判別系統(tǒng)混沌特性的動(dòng)態(tài)不變量。

    1.3.1關(guān)聯(lián)維數(shù)G-P算法

    給定一尺度r和嵌入維數(shù)m,定義相空間距離小于r的點(diǎn)對(duì)數(shù)與所有點(diǎn)對(duì)數(shù)的比例為關(guān)聯(lián)積分Cm(r),即

    (2)

    式中,n為相空間的點(diǎn)數(shù);‖·‖表示求兩個(gè)重構(gòu)向量的歐氏范數(shù);θ(x)為Heaviside階躍函數(shù)。

    (3)

    對(duì)于無(wú)限長(zhǎng)序列(n→∞)和足夠小的r值,若Cm(r)與r呈冪律關(guān)系,即Cm(r)∝rD2, 則定義關(guān)聯(lián)維數(shù)D2為

    (4)

    實(shí)際所獲取的數(shù)據(jù)量有限,無(wú)法滿足極限條件n→∞, 而且尺度r的取值須滿足式(5)。

    min{‖Vi-Vj‖}

    (5)

    式中,i≠j,i,j=1,2,…,n。

    因此,計(jì)算關(guān)聯(lián)維數(shù)采取的常規(guī)步驟為:

    1)在式(5)的范圍內(nèi)取一尺度序列ri(i=1,2,…,l),計(jì)算與ri對(duì)應(yīng)的關(guān)聯(lián)積分序列Cm(ri)。

    2)求出序列Cm(ri)和ri的對(duì)數(shù)序列l(wèi)nCm(r)i和lnri,以序列l(wèi)nri為橫坐標(biāo)、序列l(wèi)nCm(r)i為縱坐標(biāo)繪制雙對(duì)數(shù)曲線。

    3)將雙對(duì)數(shù)曲線線性區(qū)(或無(wú)標(biāo)度區(qū))的斜率作為D2的估計(jì)值。

    因雙對(duì)數(shù)曲線或多或少呈現(xiàn)S形,很難通過(guò)目測(cè)或經(jīng)驗(yàn)確定其線性區(qū)域。文獻(xiàn)[26]指出一個(gè)可信的無(wú)標(biāo)度區(qū)必須是清晰可辨的,表現(xiàn)在lnCm(r)對(duì)lnr曲線的局部斜率dsm對(duì)lnr曲線存在隨m增大而趨于飽和的平坦區(qū)域。

    為使估算所得D2盡量準(zhǔn)確,本文計(jì)算關(guān)聯(lián)維數(shù)采取的步驟與常規(guī)的不同之處為:

    1)計(jì)算一組m(m=2~10)的對(duì)數(shù)序列l(wèi)nCm(r)i。

    2)分別求出m=2~10下lnCm(r)對(duì)lnr曲線的局部斜率dsm。

    3)若dsm對(duì)lnr曲線存在隨m值增大而趨于飽和的平坦區(qū)域,則將各m值下dsm對(duì)lnr曲線的平坦區(qū)域作為無(wú)標(biāo)度區(qū),并將此區(qū)域?qū)?yīng)的數(shù)據(jù)點(diǎn)對(duì)(lnCm(r),lnr)進(jìn)行直線擬合得到序列D2_m,由此獲取D2估計(jì)值。同時(shí),D2_m趨于飽和時(shí)對(duì)應(yīng)的m值可作為該時(shí)間序列重構(gòu)的最優(yōu)m值,即G-P法確定最優(yōu)m值。

    局部斜率dsm的計(jì)算式為

    (6)

    式中,i=1,2,…,l-1。

    對(duì)于隨機(jī)過(guò)程的時(shí)間序列,其D2隨m的增大而增大,不會(huì)趨于飽和;而確定性系統(tǒng)的D2隨m的增大而逐漸趨于飽和。不同的D2飽和值對(duì)應(yīng)于不同的系統(tǒng)運(yùn)動(dòng)狀態(tài)[27]:D2=1,系統(tǒng)作周期運(yùn)動(dòng);D2=2,系統(tǒng)作準(zhǔn)周期或擬周期運(yùn)動(dòng);D2>2或D2不為整數(shù),系統(tǒng)作混沌運(yùn)動(dòng)。

    1.3.2關(guān)聯(lián)維數(shù)G-P優(yōu)化算法

    對(duì)于長(zhǎng)度為n的數(shù)據(jù)集,在所有的點(diǎn)對(duì)中直接搜尋距離小于r的點(diǎn)對(duì)所需的計(jì)算量為O(n2/2)。如果需要在多個(gè)嵌入維下估算D2,則計(jì)算耗時(shí)還要增加。由于在估算D2時(shí),只有相空間的鄰近點(diǎn)對(duì)對(duì)Cm(r)值有貢獻(xiàn)。因此,如果事先大致確定每個(gè)相點(diǎn)的近鄰分布位置,可大大縮短搜索鄰近點(diǎn)對(duì)的耗時(shí)?;谝陨纤枷?,文獻(xiàn)[28-30]先后提出采用盒子網(wǎng)格的方法幫助尋找距離小于r的鄰近點(diǎn)對(duì)。文獻(xiàn)[30]提出的基于直方圖的盒子輔助(box-assisted)近鄰搜索算法,解決了將相空間的點(diǎn)放入盒子而不浪費(fèi)存儲(chǔ)空間的技術(shù)難題,而且易于實(shí)現(xiàn)。本文應(yīng)用這一算法計(jì)算Cm(r)。

    2典型鐵磁諧振時(shí)間序列的非線性動(dòng)力學(xué)分析

    圖1 重慶市先鋒變電站一次主接線簡(jiǎn)圖Fig.1 Primary wiring simplified diagram of Xianfeng substation in Chongqing city

    本文所研究的鐵磁諧振時(shí)間序列是重慶市先鋒變電站的實(shí)測(cè)數(shù)據(jù)。該變電站的一次主接線簡(jiǎn)圖如圖1所示,圖中‘/’表示隔離開關(guān)或接地刀開關(guān)。該站有兩臺(tái)3繞組變壓器,繞組接線方式為YN/yn0/d11,電壓等級(jí)為110 kV/35 kV/10 kV,35 kV側(cè)中性點(diǎn)不接地。為監(jiān)視母線絕緣,各電壓等級(jí)的Ⅰ、Ⅱ段母線均安裝有PT。為實(shí)現(xiàn)對(duì)該站3個(gè)電壓等級(jí)的過(guò)電壓監(jiān)測(cè),分別在10 kV開關(guān)柜母線,35 kV戶外開關(guān)場(chǎng)母線,2#主變壓器套管110 kV側(cè)安裝了電壓傳感器[31]。自2006年11月16日開始實(shí)施對(duì)10 kV系統(tǒng)Ⅰ、Ⅱ段母線和2#主變壓器110 kV系統(tǒng)的過(guò)電壓監(jiān)測(cè),自2008年6月23日開始實(shí)施對(duì)35 kV系統(tǒng)Ⅰ、Ⅱ段母線的過(guò)電壓監(jiān)測(cè)。10 kV Ⅱ段母線和35 kV Ⅰ段母線上均監(jiān)測(cè)到鐵磁諧振過(guò)電壓,本文選取具有代表性的鐵磁諧振電壓序列數(shù)據(jù)進(jìn)行特征分析,為進(jìn)一步研究其分類識(shí)別方法奠定基礎(chǔ)。

    本文對(duì)原始采集數(shù)據(jù)進(jìn)行了下采樣處理和濾波去噪。①下采樣。原始數(shù)據(jù)的采樣頻率為200 kHz,對(duì)于鐵磁諧振的穩(wěn)態(tài)過(guò)程而言過(guò)高,因其所含成分的最高頻率不大于1 kHz。為使被分析的數(shù)據(jù)展現(xiàn)電壓變量更多的變化規(guī)律,又不使計(jì)算量過(guò)大,對(duì)原始數(shù)據(jù)采取整數(shù)倍抽取的方式獲得待分析的時(shí)間序列。抽取倍數(shù)的確定原則為:最低抽取頻率設(shè)為10 kHz,抽取頻率必須為50 Hz的整數(shù)倍,抽取后的相平面軌跡與未抽取前的一致。因此,本文采用實(shí)驗(yàn)法確定合適的抽取倍數(shù):以抽取倍數(shù)的取值范圍2、4、5、8、10、16、20依次進(jìn)行抽取,并對(duì)比抽取前、后的相平面軌跡,若發(fā)現(xiàn)有明顯差異,則將前一次采用的抽取倍數(shù)作為最終的信號(hào)抽取倍數(shù)。②濾波去噪。因?qū)崪y(cè)數(shù)據(jù)不可避免地含有噪音,對(duì)下采樣處理所得的時(shí)間序列采用sym8小波5層啟發(fā)式閾值估計(jì)法去噪。

    1)時(shí)間序列1

    時(shí)間序列1是2007年2月14日02∶05∶41于10 kV Ⅱ段母線處采集的A相原始數(shù)據(jù)進(jìn)行5倍抽取和濾波處理所得的數(shù)據(jù),數(shù)據(jù)長(zhǎng)度為10 000。

    圖2和圖3分別為時(shí)間序列1的波形圖和重構(gòu)相平面圖(N=10 000,τ=26)。圖2呈現(xiàn)明顯的周期性,過(guò)電壓幅值達(dá)2.802 6 (pu)。圖3所示的相平面軌跡是一條近似橢圓形的封閉曲線,表明時(shí)間序列1作周期運(yùn)動(dòng)。圖4為對(duì)應(yīng)于圖3的龐加萊截面,采樣頻率為電源頻率。本文中如不作特別說(shuō)明,對(duì)相平面軌跡采樣的頻率均與電源頻率相同。圖4中的相點(diǎn)幾乎重合,可近似為一個(gè)孤立點(diǎn),表明時(shí)間序列1作基頻周期運(yùn)動(dòng),頻率與電源頻率相同。

    圖2 時(shí)間序列1的波形圖Fig.2 Waveform plot of time series 1

    圖3 時(shí)間序列1的重構(gòu)相平面圖Fig.3 Reconstructed phase plane plot of time series 1

    圖4 時(shí)間序列1的龐加萊截面Fig.4 Poincaré section of time series 1

    圖5為時(shí)間序列1的lnCm(r)對(duì)lnr雙對(duì)數(shù)曲線圖,圖6為對(duì)應(yīng)于圖5的dsm對(duì)lnr曲線圖。圖5和圖6中顏色由淺到深的線條分別表示m=2、3、…、10下的相應(yīng)曲線,后文類似圖形中各線條的含義與此相同。尺度r的取值序列為:ri=0.06×1.311 6i-1,i=1,2,…,16。 由圖5可知,各條lnCm(r)對(duì)lnr曲線的下半部分線性度較好,且在lnr取值區(qū)間[-2.270 9,-1.185 9],m=3~10的lnCm(r)對(duì)lnr曲線基本平行。在圖6中可觀察到各條dsm對(duì)lnr曲線均存在清晰的平坦區(qū)域,且在lnr取值區(qū)間[-2.270 9,-1.475 2],m=3~10的dsm對(duì)lnr曲線基本重合。

    圖5 時(shí)間序列1的lnCm(r)對(duì)lnr曲線Fig.5 Curves lnCm(r) versus lnr of time series 1

    圖6 時(shí)間序列1的dsm對(duì)lnr曲線Fig.6 Curves dsm versus lnr of time series 1

    m=2~10時(shí),D2_m分別為1.024 3、1.013 3、1.017 5、1.018 9、1.018 8、1.017 0、1.017 4、1.013 6、1.012 2。m≥3時(shí),D2_m趨于穩(wěn)定,由此估算所得D2=1.015 0±0.003 9。該序列重構(gòu)所需的最小嵌入維數(shù)為3。

    時(shí)間序列1的D2相對(duì)于整數(shù)1的誤差低于2%,非常小,因此其D2理論值應(yīng)為整數(shù)1,作周期運(yùn)動(dòng)。這與從相平面軌跡和龐加萊截面得出的結(jié)論一致。

    對(duì)與時(shí)間序列1同一時(shí)間同一位置獲取的B相、C相數(shù)據(jù)也進(jìn)行了動(dòng)力學(xué)特征分析,其運(yùn)動(dòng)模式與時(shí)間序列1的相同,D2分別為1.007 6±0.008 7,1.028 8±0.011 5。

    根據(jù)識(shí)別結(jié)果可知,2007年2月14日 02∶05∶41于10 kV Ⅱ段母線監(jiān)測(cè)的鐵磁諧振類型是基頻周期諧振,其過(guò)電壓幅值為2.802 6(pu),一般是不危險(xiǎn)的,但若持續(xù)時(shí)間長(zhǎng),可能導(dǎo)致PT絕緣破壞,甚至燒毀或避雷器爆炸。

    2)時(shí)間序列2

    時(shí)間序列2是2009年7月16 日21∶52∶26于35 kV Ⅰ段母線處采集的A相原始數(shù)據(jù)進(jìn)行10倍抽取和濾波處理所得的數(shù)據(jù),數(shù)據(jù)長(zhǎng)度為10 000。

    圖7和圖8分別為時(shí)間序列2的波形圖和重構(gòu)相平面圖(N=10 000,τ=44)。從圖7可知,t2內(nèi)的電壓取值與t1內(nèi)的基本相同,t2之后的電壓值也有重復(fù)的跡象,可初步斷定時(shí)間序列2具有周期性。t1和t2內(nèi)有4 400個(gè)采樣點(diǎn),t1=t2=220 ms(采樣頻率為20 kHz)。過(guò)電壓幅值為-1.876 2 (pu)。圖8所示的相平面軌跡是一條由多條軌道組成的封閉曲線,不同的軌道之間有交點(diǎn)。這些特征表明時(shí)間序列2作周期運(yùn)動(dòng),其周期比較長(zhǎng)。

    圖7 時(shí)間序列2的波形圖Fig.7 Waveform plot of time series 2

    圖8 時(shí)間序列2的重構(gòu)相平面圖Fig.8 Reconstructed phase plane plot of time series 2

    圖9 時(shí)間序列2的龐加萊截面Fig.9 Poincaré section of time series 2

    與圖8對(duì)應(yīng)的龐加萊截面如圖9所示。經(jīng)仔細(xì)分析發(fā)現(xiàn),最早出現(xiàn)在龐加萊截面上的11個(gè)相點(diǎn)呈離散分布,但分布在龐加萊截面次對(duì)角線兩端的相點(diǎn)相距較近。此時(shí)僅憑圖9難以判斷時(shí)間序列2的運(yùn)動(dòng)性質(zhì)。此后獲取的相點(diǎn)依次分布在原11個(gè)相點(diǎn)的附近,甚至與原有相點(diǎn)重合。如果時(shí)間序列長(zhǎng)度為4 400(11個(gè)工頻周期)的整數(shù)倍,則這11個(gè)相點(diǎn)分布區(qū)域包含的相點(diǎn)數(shù)相同。結(jié)合相平面軌跡的分析結(jié)論,可斷定這11個(gè)相點(diǎn)區(qū)域分別表示11個(gè)孤立的離散點(diǎn)(用不同的標(biāo)記表示)。這表明時(shí)間序列2作分頻周期運(yùn)動(dòng),其頻率為工頻的1/11。受測(cè)量誤差的影響,11個(gè)孤立離散區(qū)域處相點(diǎn)的分散性不同。

    對(duì)時(shí)間序列1和2的分析表明,測(cè)量誤差的存在導(dǎo)致龐加萊截面上應(yīng)重疊的相點(diǎn)未完全重疊,給分析帶來(lái)困難。因此,對(duì)周期比較長(zhǎng)的時(shí)間序列,運(yùn)用龐加萊截面時(shí)需非常謹(jǐn)慎,時(shí)間序列數(shù)據(jù)要足夠長(zhǎng),至少要包含2個(gè)運(yùn)動(dòng)周期。

    圖10為時(shí)間序列2的lnCm(r)對(duì)lnr雙對(duì)數(shù)曲線圖。尺度r的取值序列為:ri=0.06×1.311 6i-5,i=1,2,…,16。 從圖10可知,m=2~5時(shí)lnCm(r)對(duì)lnr曲線的形狀差異較大,m≥6時(shí)lnCm(r)對(duì)lnr曲線的形狀逐漸一致,且m=9與m=10的基本重合。與圖10相對(duì)應(yīng)的dsm對(duì)lnr曲線圖如圖11所示。從圖11可觀察到,m≥6時(shí)dsm對(duì)lnr曲線均存在明顯的平坦區(qū)域,且在lnr取值區(qū)間[-3.084 7,-2.270 9],dsm對(duì)lnr曲線平坦且基本重合。

    圖10 時(shí)間序列2的lnCm(r)對(duì)lnr曲線Fig.10 Curves lnCm(r) versus lnr of time series 2

    圖11 時(shí)間序列2的dsm對(duì)lnr曲線Fig.11 Curves dsm versus lnr of time series 2

    m=6~10時(shí),D2_m分別為1.006 3、1.006 7、1.006 7、1.005 7、1.005 6,波動(dòng)極小,由此估算可得D2=1.006 1±0.000 6。該序列重構(gòu)所需的最小嵌入維數(shù)為6。

    時(shí)間序列2的D2相對(duì)于整數(shù)1的誤差低于0.7%,非常小,因此其D2理論值應(yīng)為整數(shù)1,作周期運(yùn)動(dòng)。這與從相平面軌跡和龐加萊截面得出的結(jié)論一致。

    對(duì)與時(shí)間序列2同一時(shí)間同一位置獲取的B相、C相數(shù)據(jù)也進(jìn)行了動(dòng)力學(xué)特征分析,其運(yùn)動(dòng)模式與時(shí)間序列2的相同,D2分別為1.001 3±0.001 0,1.018 4±0.005 3。

    根據(jù)識(shí)別結(jié)果可知,2009年7月16日21∶52∶26于35 kV Ⅰ段母線監(jiān)測(cè)的鐵磁諧振類型是分頻周期諧振,過(guò)電壓幅值為-1.876 2 (pu),電壓本身是不危險(xiǎn)的,但分頻引起的PT過(guò)飽和電流,若長(zhǎng)時(shí)間持續(xù),易導(dǎo)致PT高壓熔絲熔斷,或引起PT因嚴(yán)重過(guò)熱而燒毀,甚至爆炸。

    3)時(shí)間序列3

    時(shí)間序列3是2009年4月18日01∶29∶58于35 kV Ⅰ段母線處采集的A相原始數(shù)據(jù)進(jìn)行10倍抽取和濾波處理所得的數(shù)據(jù),數(shù)據(jù)長(zhǎng)度為10 000。

    圖12和圖13分別為時(shí)間序列3的波形圖和重構(gòu)相平面圖(N=10 000,τ=57)。時(shí)間序列3的各波峰之間呈現(xiàn)此消彼長(zhǎng)的現(xiàn)象,有一定的重復(fù)性,但又不完全重復(fù)。過(guò)電壓幅值為-2.127 5 (pu)。圖13呈現(xiàn)的相平面軌跡不是一條封閉曲線,但相點(diǎn)的運(yùn)動(dòng)軌道被限制在一個(gè)固定的區(qū)域,即其運(yùn)動(dòng)軌跡并不是無(wú)限制擴(kuò)散和完全雜亂無(wú)章。這些特征表明時(shí)間序列3不是作周期運(yùn)動(dòng),也不是作隨機(jī)運(yùn)動(dòng)。

    圖12 時(shí)間序列3的波形圖Fig.12 Waveform plot of time series 3

    圖13 時(shí)間序列3的重構(gòu)相平面圖Fig.13 Reconstructed phase plane plot of time series 3

    與圖13對(duì)應(yīng)的龐加萊截面如圖14所示。圖14的相點(diǎn)呈不均勻的離散分布,未形成一條近似閉曲線,這表明時(shí)間序列3不是作準(zhǔn)周期運(yùn)動(dòng)。因圖14中的相點(diǎn)數(shù)較少,未呈現(xiàn)自相似的分層結(jié)構(gòu),因此不能判斷時(shí)間序列3是否作混沌運(yùn)動(dòng)。

    圖14 時(shí)間序列3的龐加萊截面Fig.14 Poincaré section of time series 3

    時(shí)間序列3的lnCm(r)對(duì)lnr雙對(duì)數(shù)曲線圖如圖15所示,對(duì)應(yīng)于圖15的dsm對(duì)lnr曲線圖如圖16所示。尺度r的取值序列為:ri=0.06×1.3116i-1,i=1,2,…,16。 從圖15可知,在lnr取值區(qū)間[-0.643 4,0.170 3],雙對(duì)數(shù)曲線隨m的增大而逐漸平行,即斜率值逐漸達(dá)到飽和。相應(yīng)地,在圖16中可清晰的觀察到dsm對(duì)lnr曲線存在隨m的增大而逐漸飽和的平坦區(qū)域,且m=8~10時(shí)dsm對(duì)lnr曲線的平坦區(qū)域基本重合。

    圖15 時(shí)間序列3的lnCm(r)對(duì)lnr曲線Fig.15 Curves lnCm(r) versus lnr of time series 3

    圖16 時(shí)間序列3的dsm對(duì)lnr曲線Fig.16 Curves dsm versus lnr of time series 3

    m=2~10時(shí),D2_m分別為1.678 4、1.837 0、2.118 5、2.341 2、2.361 4、2.264 4、2.238 9、2.302 4、2.312 1。m≥5時(shí),D2_m趨于穩(wěn)定,由此估算可得D2=2.300 2±0.061 2。該序列重構(gòu)所需的最小嵌入維數(shù)為5。

    時(shí)間序列3的D2>2,不為整數(shù),且小數(shù)部分較大(>0.230 0),表明其具有混沌運(yùn)動(dòng)特性。從龐加萊截面得出的其不是作準(zhǔn)周期運(yùn)動(dòng)的判斷,可驗(yàn)證這一結(jié)論。

    對(duì)與時(shí)間序列3同一時(shí)間同一位置獲取的B相、C相數(shù)據(jù)也進(jìn)行了動(dòng)力學(xué)特征分析,其運(yùn)動(dòng)模式與時(shí)間序列3的相同,D2分別為2.223 6±0.043 6,2.408 3±0.082 4。

    根據(jù)識(shí)別結(jié)果可知,2009年4月18日01∶29∶58于35 kV Ⅰ段母線監(jiān)測(cè)的鐵磁諧振類型是混沌諧振,其過(guò)電壓幅值僅稍微超過(guò)2 (pu),電壓本身也是不危險(xiǎn)的。對(duì)時(shí)間序列3的頻譜分析得知,其含24 Hz頻率分量的幅值達(dá)1.161 4 (pu),即含較大幅值的分次諧波。因此,若長(zhǎng)時(shí)間持續(xù),有與分頻諧振類似的危害。

    3結(jié)論

    1)為準(zhǔn)確識(shí)別實(shí)測(cè)鐵磁諧振一維電壓時(shí)間序列所屬的諧振類型,本文運(yùn)用相空間重構(gòu)方法構(gòu)造系統(tǒng)相空間。在此基礎(chǔ)上,采用相平面、龐加萊截面和關(guān)聯(lián)維數(shù)3種非線性動(dòng)力學(xué)分析方法,以不同的形式表達(dá)其動(dòng)力學(xué)特征,由此對(duì)其所屬的鐵磁諧振類型作出了準(zhǔn)確判斷。

    2)對(duì)3例典型的中性點(diǎn)不接地系統(tǒng)PT鐵磁諧振實(shí)測(cè)時(shí)間序列進(jìn)行了動(dòng)力學(xué)特征分析。其中兩例時(shí)間序列的相平面軌跡、龐加萊截面均呈現(xiàn)周期運(yùn)動(dòng)特性,關(guān)聯(lián)維數(shù)估計(jì)值分別為1.015 0±0.003 9、1.006 1±0.000 6,由此判斷它們作周期運(yùn)動(dòng)。另一例時(shí)間序列的相平面軌跡、龐加萊截面均未呈現(xiàn)周期或準(zhǔn)周期運(yùn)動(dòng)特性,關(guān)聯(lián)維數(shù)估計(jì)值為2.300 2±0.061 2,由此判斷其作混沌運(yùn)動(dòng)。

    3)鑒于實(shí)測(cè)電壓時(shí)間序列數(shù)據(jù)受到多種因素(如采樣精度、噪音和誤差等)的影響,龐加萊截面和關(guān)聯(lián)維數(shù)估計(jì)值與理想情況有一定差距,要對(duì)其運(yùn)動(dòng)特性作出準(zhǔn)確判斷,有必要綜合多種方法所表達(dá)的特征。

    4)中性點(diǎn)不接地系統(tǒng)PT鐵磁諧振屬于零序性質(zhì),發(fā)生鐵磁諧振時(shí)對(duì)任一相電壓數(shù)據(jù)進(jìn)行非線性動(dòng)力學(xué)分析,均能對(duì)其類型作出準(zhǔn)確判斷。

    5)本文的研究為鐵磁諧振類型的在線自動(dòng)識(shí)別和諧振過(guò)電壓的抑制策略研究奠定了一定的基礎(chǔ),為實(shí)測(cè)鐵磁諧振數(shù)據(jù)的識(shí)別研究提供了新思路,豐富了鐵磁諧振的研究實(shí)例。

    參考文獻(xiàn)

    [1]Emin Z,Al Zahawi B A T,Auckland D W,et al.Ferroresonance in electromagnetic voltage transformers:a study based on nonlinear dynamics[J].IEE Proceedings-Generation,Transmission and Distribution,1997,144(4):383-387.

    [2]Mork B A,Stuehm D L.Application of nonlinear dynamics and chaos to ferroresonance in distribution systems[J].IEEE Transactions on Power Delivery,1994,9(2):1009-1017.

    [3]周麗霞,尹忠東,鄭立.配網(wǎng)PT鐵磁諧振機(jī)理與抑制的試驗(yàn)研究[J].電工技術(shù)學(xué)報(bào),2007,22(5):153-158.

    Zhou Lixia,Yin Zhongdong,Zheng Li.Research on principle of PT resonance in distribution power system and its suppression[J].Transactions of China Electrotechnical Society,2007,22(5):153-158.

    [4]張博,魯鐵成,杜曉磊.中性點(diǎn)接地系統(tǒng)鐵磁諧振非線性動(dòng)力學(xué)分析[J].高電壓技術(shù),2007,33(1):31-35.

    Zhang Bo,Lu Tiecheng,Du Xiaolei.Nonlinear dynamic analysis of ferroresonance in neutral-grounded power system[J].High Voltage Engineering,2007,33(1):31-35.

    [5]Val E M,Dudurych I,Redfern M A.Characterization of ferroresonant modes in HV substation with CB grading capacitors[J].Electric Power Systems Research,2007,77(11):1506-1513.

    [6]Araujo A E A,Soudack A C,Marti J R.Ferroresonance in power systems:chaotic behaviour[J].Generation,Transmission and Distribution,IEE Proceedings C,1993,140(3):237-240.

    [7]Emin Z,Al Zahawi B A T,Yu K T,et al.Quantification of the chaotic behavior of ferroresonant voltage transformer circuits[J].Circuits and Systems I:Fundamental IEEE Transactions on Theory and Applications,2001,48(6):757-760.

    [8]Mozaffari S,Henschel S,Soudack A C.Chaotic ferroresonance in power transformers[J].IEE Proceedings-Generation,Transmission and Distribution,1995,142(3):247-250.

    [9]Kieny C.Application of the bifurcation theory in studying and understanding the global behavior of a ferroresonant electric power circuit[J].IEEE Transactions on Power Delivery,1991,6(2):866-872.

    [10]司馬文霞,鄭哲人,楊慶,等.用參數(shù)不匹配混沌系統(tǒng)的脈沖同步方法抑制鐵磁諧振過(guò)電壓[J].電工技術(shù)學(xué)報(bào),2012,27(6):218-225.

    Sima Wenxia,Zheng Zheren,Yang Qing,et al.Suppressing the ferroresonance overvoltage by impulsive synchronization of chaotic systems with parameter mismatched[J].Transactions of China Electrotechnical Society,2012,27(6):218-225.

    [11]Takens F.Detecting strange attractors in turbulence[M]//Lecture Notes in Mathematica.Berlin:Springer Press,1981:366-381.

    [12]Packard N H,Crutchfield J P,F(xiàn)armer J D,et al.Geometry from a time series[J].Physical Review Letters,1980,45:712-715.

    [13]行鴻彥,龔平,徐偉.嵌入窗方法確定混沌系統(tǒng)重構(gòu)參數(shù)的仿真研究[J].系統(tǒng)仿真學(xué)報(bào),2013,25(6):1219-1225.

    Xing Hongyan,Gong Ping,Xu Wei.Simulation research of chaos system reconstruction parameters based on embedded window[J].Journal of System Simulation,2013,25(6):1219-1225.

    [14]Albano A M,Muench J,Schwartz C,et al.Singular-value decomposition and the Grassberger-Procaccia algorithm[J].Physical Review A,1988,38(6):3017-3026.

    [15]Fraser A M,Swinney H L.Independent coordinates for strange attractors from mutual information[J].Physical Review A,1986,33(2):1134-1440.

    [16]Grassberger P,Procaccia I.Measuring the strangeness of strange attractors[J].Physica D Nonlinear Phenomena,1983,9(1/2):189-208.

    [17]Broomhead D S,King G P.Extracting qualitative dynamics from experimental data[J].Physica D,1986,20(2/3):217-236.

    [18]Kennel M B,Brown R,Abarbanel H D I.Determining embedding dimension for phase-space reconstruction using a geometrical construction[J].Physical Review A,1992,45(6):3403-3411.

    [19]Cao L Y.Practical method for determining the minimum embedding dimension of a scalar time series[J].Physica D,1997,110(1/2):43-50.

    [20]Kugiumtzis D.State space reconstruction parameters in the analysis of chaotic time series—the role of the time window length[J].Physica D:Nonlinear Phenomena,1996,95(1):13-28.

    [21]Kim H S,Eykholt R,Salas J D.Nonlinear dynamics,delay times,and embedding windows[J].Physica D:Nonlinear Phenomena,1999,127(1/2):48-60.

    [22]Otani M,Jones A.Automated embedding and creep phenomenon in chaotic time series[EB/OL].http://users.cs.cf.ac.uk/Antonia.J.Jones/UnpublishedPapers/Creep.pdf/,2000-10-14/2011-03-15.

    [23]Franca L,Savi M A.Distinguishing periodic and chaotic time series obtained from an experimental nonlinear pendulum[J].Nonlinear Dynamics,2001,26(3):253-271.

    [24]劉秉正,彭建華.非線性動(dòng)力學(xué)[M].北京:高等教育出版社,2005.

    [25]Grassberge P,Procaccia I.Characterization of strange attractors[J].Physical Review Letters,1983,50(5):346-349.

    [26]Kantz H,Schreiber T.Nonlinear time series analysis[M].Cambridge:Cambridge University Press,1997.

    [27]張雨.時(shí)間序列的混沌和符號(hào)分析及實(shí)踐[M].長(zhǎng)沙:國(guó)防科技大學(xué)出版社,2007.

    [28]Theiler J.Efficient algorithm for estimating the correlation dimension from a set of discrete points[J].Physical Review A,1987,36(9):4456-4462.

    [29]Grassberger P.An Optimized box-assisted algorithm for fractal dimensions[J].Physics Letters A,1990,148(1/2):63-68.

    [30]Schreiber T.Efficient neighbor searching in nonlinear time-series analysis[J].International Journal of Bifurcation and Chaos,1995,5(2):349-358.

    [31]杜林,李欣,司馬文霞,等.110 kV變電站過(guò)電壓在線監(jiān)測(cè)系統(tǒng)及其波形分析[J].高電壓技術(shù),2012,38(3):535-543.

    Du Lin,Li Xin,Sima Wenxia,et al.Overvoltage on-line monitoring system for 110 kV substation and its waveforms analysis[J].High Voltage Engineering,2012,38(3):535-543.

    黃艷玲女,1978年生,博士,講師,研究方向?yàn)殡娏ο到y(tǒng)過(guò)電壓。

    E-mail:huangyanling@cqu.edu.cn(通信作者)

    司馬文霞女,1965年生,教授,博士生導(dǎo)師,研究方向?yàn)殡娏ο到y(tǒng)的防雷與過(guò)電壓防護(hù)、特殊環(huán)境中外絕緣放電特性及機(jī)理。

    E-mail:cqsmwx@cqu.edu.cn

    作者簡(jiǎn)介

    中圖分類號(hào):TM864

    收稿日期2015-03-04改稿日期2015-06-21

    國(guó)家自然科學(xué)基金(51177182)、中央高?;究蒲袠I(yè)務(wù)費(fèi)(0209005206004)和國(guó)家創(chuàng)新研究群體科學(xué)基金(51321063)資助。

    猜你喜歡
    龐加萊相空間鐵磁
    關(guān)于兩類多分量海森堡鐵磁鏈模型的研究
    龐加萊偵察術(shù)
    中外文摘(2022年7期)2022-05-17 09:36:42
    龐加萊偵察術(shù)
    青年文摘(2021年20期)2021-12-11 18:45:12
    龐加萊偵查術(shù)
    束團(tuán)相空間分布重建技術(shù)在西安200 MeV質(zhì)子應(yīng)用裝置的應(yīng)用
    你好,鐵磁
    非對(duì)易空間中的三維諧振子Wigner函數(shù)
    你好,鐵磁
    基于相空間重構(gòu)的電磁繼電器電性能參數(shù)預(yù)測(cè)研究
    推挽式直流變換器的龐加萊映射圖分析
    热99re8久久精品国产| 成人二区视频| 青春草亚洲视频在线观看| 亚洲成人中文字幕在线播放| 国产精品福利在线免费观看| 亚洲成人精品中文字幕电影| 精华霜和精华液先用哪个| 精品人妻视频免费看| 日韩欧美在线乱码| 免费av观看视频| 免费看光身美女| 伦精品一区二区三区| 97超碰精品成人国产| 最近最新中文字幕大全电影3| 蜜桃久久精品国产亚洲av| 中文字幕制服av| 国产91av在线免费观看| 国产在线一区二区三区精 | 国产69精品久久久久777片| 日本免费在线观看一区| 18禁裸乳无遮挡免费网站照片| 国产精品福利在线免费观看| 久久99蜜桃精品久久| av播播在线观看一区| 国内少妇人妻偷人精品xxx网站| 国产一级毛片七仙女欲春2| 两个人的视频大全免费| 亚洲av中文av极速乱| 五月玫瑰六月丁香| 日本欧美国产在线视频| 尤物成人国产欧美一区二区三区| av福利片在线观看| 欧美精品国产亚洲| 久久久久性生活片| 午夜视频国产福利| 欧美日本亚洲视频在线播放| 日日干狠狠操夜夜爽| 免费看av在线观看网站| 大香蕉久久网| 免费电影在线观看免费观看| 精品欧美国产一区二区三| 直男gayav资源| 欧美激情久久久久久爽电影| 亚洲av男天堂| 精品一区二区三区视频在线| 波多野结衣巨乳人妻| 欧美一区二区精品小视频在线| 中文精品一卡2卡3卡4更新| 国产精华一区二区三区| 小蜜桃在线观看免费完整版高清| 亚洲真实伦在线观看| 久久久精品欧美日韩精品| 午夜激情福利司机影院| 久久精品国产亚洲av涩爱| 黄片无遮挡物在线观看| 国产免费视频播放在线视频 | 秋霞在线观看毛片| 久久久久久久久久久丰满| 亚洲在线自拍视频| 天天一区二区日本电影三级| 亚洲成色77777| www.av在线官网国产| 大香蕉久久网| 男的添女的下面高潮视频| 国产精品一区二区三区四区久久| 水蜜桃什么品种好| 麻豆成人av视频| 日本与韩国留学比较| 超碰av人人做人人爽久久| 女人被狂操c到高潮| 国产免费男女视频| 国产不卡一卡二| 久久韩国三级中文字幕| 国产免费又黄又爽又色| 波多野结衣高清无吗| 日韩强制内射视频| 99久国产av精品| 亚洲在久久综合| 精华霜和精华液先用哪个| 日日摸夜夜添夜夜爱| 九九爱精品视频在线观看| 亚洲经典国产精华液单| 亚洲av成人av| 免费黄色在线免费观看| 人妻制服诱惑在线中文字幕| 日韩欧美三级三区| 国产中年淑女户外野战色| 18+在线观看网站| 日本五十路高清| 偷拍熟女少妇极品色| 婷婷色av中文字幕| 小说图片视频综合网站| 精品久久国产蜜桃| 嘟嘟电影网在线观看| 国产精品国产高清国产av| 亚洲美女视频黄频| 国产淫片久久久久久久久| 赤兔流量卡办理| 日产精品乱码卡一卡2卡三| 久久综合国产亚洲精品| 亚洲不卡免费看| 夫妻性生交免费视频一级片| 亚洲18禁久久av| 精品人妻熟女av久视频| www.av在线官网国产| 久久99热这里只有精品18| 最后的刺客免费高清国语| 亚洲av电影不卡..在线观看| 中文字幕人妻熟人妻熟丝袜美| 成年女人永久免费观看视频| 久久久精品94久久精品| 在现免费观看毛片| 韩国高清视频一区二区三区| 亚洲精品乱码久久久v下载方式| 九九爱精品视频在线观看| 国产毛片a区久久久久| 在线天堂最新版资源| 国产真实伦视频高清在线观看| 欧美又色又爽又黄视频| 草草在线视频免费看| 精品午夜福利在线看| 九色成人免费人妻av| 国产亚洲最大av| 超碰av人人做人人爽久久| 少妇猛男粗大的猛烈进出视频 | 久久久久国产网址| 免费av不卡在线播放| 汤姆久久久久久久影院中文字幕 | eeuss影院久久| 一区二区三区免费毛片| 国产成人午夜福利电影在线观看| 色噜噜av男人的天堂激情| 一级毛片我不卡| 青春草视频在线免费观看| 久久欧美精品欧美久久欧美| 综合色丁香网| 亚洲成人久久爱视频| 日本黄大片高清| 在线观看av片永久免费下载| 婷婷六月久久综合丁香| 亚洲自偷自拍三级| 日日干狠狠操夜夜爽| 韩国av在线不卡| 中文字幕av在线有码专区| 欧美日韩一区二区视频在线观看视频在线 | 国产又黄又爽又无遮挡在线| 可以在线观看毛片的网站| 精品国产一区二区三区久久久樱花 | 日韩欧美精品免费久久| 少妇的逼水好多| www日本黄色视频网| 成人无遮挡网站| 夜夜爽夜夜爽视频| 国产亚洲5aaaaa淫片| 一二三四中文在线观看免费高清| 国产色婷婷99| 自拍偷自拍亚洲精品老妇| 国产成人a∨麻豆精品| 欧美一级a爱片免费观看看| 男人舔奶头视频| 91在线精品国自产拍蜜月| av黄色大香蕉| 久久久亚洲精品成人影院| 久久99热这里只频精品6学生 | 日韩中字成人| 在线免费十八禁| 国产伦一二天堂av在线观看| 亚洲av福利一区| 国产视频内射| 成年女人看的毛片在线观看| 两个人视频免费观看高清| 国产大屁股一区二区在线视频| 中文天堂在线官网| 国产成人a∨麻豆精品| 精品国内亚洲2022精品成人| 亚洲欧美成人综合另类久久久 | 国产一区二区三区av在线| 精品久久久久久久末码| 国产精品野战在线观看| 国产不卡一卡二| 久久久久网色| 免费av不卡在线播放| 亚洲欧洲日产国产| 中文字幕人妻熟人妻熟丝袜美| 又粗又爽又猛毛片免费看| 欧美zozozo另类| 国产成人精品一,二区| 欧美又色又爽又黄视频| 一本一本综合久久| 色综合色国产| 噜噜噜噜噜久久久久久91| 久久精品国产亚洲av天美| 两个人的视频大全免费| av国产久精品久网站免费入址| 欧美潮喷喷水| 久久久久久伊人网av| 韩国av在线不卡| 夫妻性生交免费视频一级片| 深爱激情五月婷婷| 啦啦啦观看免费观看视频高清| 亚洲欧美日韩无卡精品| 国产老妇伦熟女老妇高清| 国产免费男女视频| 久久久久久大精品| 久久久久久久亚洲中文字幕| 男女下面进入的视频免费午夜| 日韩亚洲欧美综合| 一级黄色大片毛片| 一个人免费在线观看电影| 亚洲人成网站高清观看| 国内揄拍国产精品人妻在线| 99热这里只有是精品在线观看| 亚洲成人av在线免费| a级毛片免费高清观看在线播放| 亚洲人成网站在线播| 欧美一区二区精品小视频在线| 亚洲av熟女| 亚洲丝袜综合中文字幕| 我的女老师完整版在线观看| 99久国产av精品| 亚洲精品aⅴ在线观看| 久久6这里有精品| 成人一区二区视频在线观看| 久久久久久伊人网av| 韩国高清视频一区二区三区| 久久久久久九九精品二区国产| av视频在线观看入口| 成人欧美大片| 免费不卡的大黄色大毛片视频在线观看 | 久久久久久久国产电影| 边亲边吃奶的免费视频| 男插女下体视频免费在线播放| 91精品一卡2卡3卡4卡| 久久久久久久久大av| 建设人人有责人人尽责人人享有的 | 丰满人妻一区二区三区视频av| 桃色一区二区三区在线观看| 国产乱人视频| 亚洲va在线va天堂va国产| 久久久欧美国产精品| 久久久久久久久大av| 欧美3d第一页| 亚洲av中文av极速乱| 最近手机中文字幕大全| 亚洲成色77777| 成人漫画全彩无遮挡| 国产高清三级在线| 亚洲伊人久久精品综合 | 亚洲中文字幕一区二区三区有码在线看| 观看美女的网站| 99热精品在线国产| 黑人高潮一二区| 亚洲av不卡在线观看| 国产精品综合久久久久久久免费| 午夜免费激情av| 久久久久国产网址| 日韩成人av中文字幕在线观看| 夜夜爽夜夜爽视频| 美女国产视频在线观看| 亚洲不卡免费看| 亚洲精品,欧美精品| 亚洲最大成人手机在线| 老司机影院毛片| 一级毛片久久久久久久久女| 久久久久久久久中文| 色吧在线观看| 高清在线视频一区二区三区 | 午夜激情欧美在线| 中文天堂在线官网| 久久久色成人| 人体艺术视频欧美日本| 色播亚洲综合网| 欧美高清成人免费视频www| 国产高清视频在线观看网站| 99热这里只有精品一区| 精华霜和精华液先用哪个| 午夜激情欧美在线| 日韩人妻高清精品专区| 久久久国产成人免费| av天堂中文字幕网| 99久久无色码亚洲精品果冻| 看免费成人av毛片| 一夜夜www| 亚洲人成网站在线观看播放| 三级国产精品片| 我的老师免费观看完整版| 亚洲五月天丁香| 免费看光身美女| 噜噜噜噜噜久久久久久91| 夜夜爽夜夜爽视频| 国产精品一区二区三区四区久久| 十八禁国产超污无遮挡网站| 美女国产视频在线观看| 国产乱人偷精品视频| www日本黄色视频网| 青春草国产在线视频| 欧美精品一区二区大全| 国产精品久久久久久精品电影小说 | 中国国产av一级| 日韩视频在线欧美| 亚洲国产欧洲综合997久久,| 久久久久久久亚洲中文字幕| 亚洲av免费在线观看| 亚洲国产精品久久男人天堂| 少妇被粗大猛烈的视频| 18禁动态无遮挡网站| 国产精品福利在线免费观看| 国产麻豆成人av免费视频| 国产极品精品免费视频能看的| 性插视频无遮挡在线免费观看| 免费看光身美女| 亚洲av电影在线观看一区二区三区 | 在线免费十八禁| 亚洲国产高清在线一区二区三| 亚洲五月天丁香| 欧美激情久久久久久爽电影| 一二三四中文在线观看免费高清| 99久久九九国产精品国产免费| 一个人免费在线观看电影| 在线播放无遮挡| 成年女人永久免费观看视频| 午夜视频国产福利| 国产视频首页在线观看| 欧美高清成人免费视频www| 国产亚洲午夜精品一区二区久久 | 国产久久久一区二区三区| 91aial.com中文字幕在线观看| 在线免费观看不下载黄p国产| 中文天堂在线官网| 亚洲成人久久爱视频| 人妻夜夜爽99麻豆av| 日本一二三区视频观看| 寂寞人妻少妇视频99o| 日本一二三区视频观看| 寂寞人妻少妇视频99o| 午夜免费激情av| 亚洲成人av在线免费| 嫩草影院入口| 国产高清不卡午夜福利| 国语对白做爰xxxⅹ性视频网站| 欧美+日韩+精品| 久久久久久久久久黄片| 中文字幕精品亚洲无线码一区| 久久久久久久久久成人| 精品久久久久久电影网 | 国产av码专区亚洲av| 久久久久久久午夜电影| 中文字幕制服av| 中文天堂在线官网| 卡戴珊不雅视频在线播放| 小蜜桃在线观看免费完整版高清| www.色视频.com| 水蜜桃什么品种好| 不卡视频在线观看欧美| 亚洲av.av天堂| 免费在线观看成人毛片| 看十八女毛片水多多多| 国产精品三级大全| 国产91av在线免费观看| 日韩av不卡免费在线播放| 一二三四中文在线观看免费高清| 亚洲丝袜综合中文字幕| 亚洲美女视频黄频| 少妇熟女欧美另类| 一级黄色大片毛片| 99九九线精品视频在线观看视频| 级片在线观看| 日本与韩国留学比较| 男人和女人高潮做爰伦理| 国产免费男女视频| 精品久久久久久久久av| 欧美高清性xxxxhd video| 欧美一区二区精品小视频在线| 免费搜索国产男女视频| 国产精品麻豆人妻色哟哟久久 | av国产免费在线观看| 日本一二三区视频观看| 高清视频免费观看一区二区 | 我要看日韩黄色一级片| 老师上课跳d突然被开到最大视频| 国产极品天堂在线| 99久久精品国产国产毛片| 不卡视频在线观看欧美| 亚洲内射少妇av| 亚洲国产精品久久男人天堂| 国产亚洲5aaaaa淫片| 亚洲精品影视一区二区三区av| 精品久久久久久久人妻蜜臀av| 男女国产视频网站| 欧美激情国产日韩精品一区| 国产免费视频播放在线视频 | av免费在线看不卡| 免费观看精品视频网站| or卡值多少钱| 国产成人freesex在线| 在线观看美女被高潮喷水网站| 一个人看视频在线观看www免费| 日日干狠狠操夜夜爽| 国产精品.久久久| 嫩草影院入口| 在线天堂最新版资源| 99久国产av精品国产电影| 亚洲av日韩在线播放| 午夜福利视频1000在线观看| 成人毛片a级毛片在线播放| 最近2019中文字幕mv第一页| 欧美日本亚洲视频在线播放| 欧美另类亚洲清纯唯美| 日韩欧美三级三区| 精品久久久久久久人妻蜜臀av| 舔av片在线| 国产探花极品一区二区| 天天躁日日操中文字幕| 亚洲欧美精品自产自拍| 精品无人区乱码1区二区| 日产精品乱码卡一卡2卡三| 久久久国产成人精品二区| 国产不卡一卡二| 国产av一区在线观看免费| 最近最新中文字幕免费大全7| 亚洲在线观看片| 97在线视频观看| 色5月婷婷丁香| 亚洲av熟女| 在线观看66精品国产| 大香蕉97超碰在线| 国产乱人视频| 极品教师在线视频| 七月丁香在线播放| a级一级毛片免费在线观看| 欧美性感艳星| 亚洲成色77777| 国产亚洲精品久久久com| ponron亚洲| 99热这里只有是精品在线观看| 97人妻精品一区二区三区麻豆| 欧美成人午夜免费资源| 亚洲av.av天堂| 大又大粗又爽又黄少妇毛片口| 丰满乱子伦码专区| 国产成人freesex在线| 精品久久久久久久末码| 1000部很黄的大片| 建设人人有责人人尽责人人享有的 | av在线播放精品| 黄色日韩在线| 亚洲国产色片| 国产精品福利在线免费观看| 精品久久久久久成人av| 国产精品一区二区三区四区久久| 国产精品不卡视频一区二区| 男女视频在线观看网站免费| 三级国产精品欧美在线观看| 三级国产精品片| 久久精品久久久久久噜噜老黄 | 久久久精品欧美日韩精品| 国内精品一区二区在线观看| 亚洲欧美清纯卡通| 亚洲av中文av极速乱| 国产精品蜜桃在线观看| 真实男女啪啪啪动态图| 有码 亚洲区| 欧美人与善性xxx| 亚洲婷婷狠狠爱综合网| 国产精品综合久久久久久久免费| 又爽又黄无遮挡网站| 国产精品av视频在线免费观看| 亚洲熟妇中文字幕五十中出| 免费人成在线观看视频色| av在线天堂中文字幕| 看十八女毛片水多多多| 麻豆成人午夜福利视频| 深夜a级毛片| 午夜久久久久精精品| 国产黄a三级三级三级人| 联通29元200g的流量卡| 亚洲国产欧洲综合997久久,| 午夜精品一区二区三区免费看| 国产伦一二天堂av在线观看| 亚洲在久久综合| 两个人视频免费观看高清| 免费无遮挡裸体视频| 欧美性猛交╳xxx乱大交人| 成人一区二区视频在线观看| 亚州av有码| 夜夜爽夜夜爽视频| 亚洲熟妇中文字幕五十中出| 亚洲人成网站在线播| 男女下面进入的视频免费午夜| 亚洲美女搞黄在线观看| 久久久成人免费电影| 乱码一卡2卡4卡精品| 色吧在线观看| 又黄又爽又刺激的免费视频.| 国产精品人妻久久久影院| 亚洲天堂国产精品一区在线| 99热网站在线观看| 国产大屁股一区二区在线视频| 亚洲av男天堂| 国产一区有黄有色的免费视频 | 免费看光身美女| 日本av手机在线免费观看| 欧美日韩国产亚洲二区| 青春草亚洲视频在线观看| 人人妻人人看人人澡| 九九热线精品视视频播放| 免费一级毛片在线播放高清视频| 国产成人福利小说| 日韩一区二区视频免费看| 亚洲欧洲日产国产| 国产欧美日韩精品一区二区| 亚洲内射少妇av| 波多野结衣巨乳人妻| 综合色av麻豆| 久久亚洲国产成人精品v| 狂野欧美激情性xxxx在线观看| 国产成人午夜福利电影在线观看| 久久久精品欧美日韩精品| 日韩成人伦理影院| 免费观看人在逋| 日韩,欧美,国产一区二区三区 | 亚洲国产精品合色在线| 日日摸夜夜添夜夜爱| 欧美不卡视频在线免费观看| 久久久久网色| 26uuu在线亚洲综合色| 久久精品熟女亚洲av麻豆精品 | 日韩亚洲欧美综合| 亚洲精品成人久久久久久| 啦啦啦观看免费观看视频高清| 国产黄片视频在线免费观看| 午夜福利网站1000一区二区三区| 欧美极品一区二区三区四区| 九色成人免费人妻av| 岛国毛片在线播放| 国产成人aa在线观看| 亚洲美女视频黄频| kizo精华| 午夜福利视频1000在线观看| 亚洲国产精品sss在线观看| 国产日韩欧美在线精品| 国产高清有码在线观看视频| 永久免费av网站大全| 中文字幕免费在线视频6| 男女边吃奶边做爰视频| 夫妻性生交免费视频一级片| 三级国产精品片| 在线免费观看的www视频| 亚洲性久久影院| 毛片女人毛片| 少妇被粗大猛烈的视频| 成人午夜精彩视频在线观看| 国产成人a∨麻豆精品| 亚洲精品自拍成人| 国产亚洲av片在线观看秒播厂 | 国语对白做爰xxxⅹ性视频网站| 自拍偷自拍亚洲精品老妇| 又黄又爽又刺激的免费视频.| 一夜夜www| 一区二区三区高清视频在线| 日本三级黄在线观看| 免费观看精品视频网站| 少妇猛男粗大的猛烈进出视频 | 午夜福利网站1000一区二区三区| 久久久久久久国产电影| 91久久精品国产一区二区三区| 精品国产一区二区三区久久久樱花 | 在线播放无遮挡| 午夜福利视频1000在线观看| www日本黄色视频网| 麻豆一二三区av精品| 18禁在线无遮挡免费观看视频| 亚洲精品久久久久久婷婷小说 | 菩萨蛮人人尽说江南好唐韦庄 | 国产精品日韩av在线免费观看| 国产精品1区2区在线观看.| 国产视频内射| 亚洲av成人av| 国产精品1区2区在线观看.| 成年免费大片在线观看| 国产成人精品婷婷| 国产精品1区2区在线观看.| 久久99热6这里只有精品| 亚洲av成人av| 免费av不卡在线播放| 欧美丝袜亚洲另类| 99久久中文字幕三级久久日本| 两性午夜刺激爽爽歪歪视频在线观看| 乱系列少妇在线播放| 亚洲综合色惰| 免费看光身美女| 最近的中文字幕免费完整| 我的女老师完整版在线观看| 伊人久久精品亚洲午夜| 1024手机看黄色片| 精品久久久久久电影网 | 亚洲欧美清纯卡通| 看非洲黑人一级黄片| 免费av观看视频| 精品欧美国产一区二区三| 亚洲怡红院男人天堂| 欧美激情国产日韩精品一区| 精品人妻熟女av久视频| 色吧在线观看| 精品酒店卫生间| 桃色一区二区三区在线观看| 99在线人妻在线中文字幕| 欧美日韩综合久久久久久| 欧美成人精品欧美一级黄| 九九热线精品视视频播放| 亚洲欧美精品自产自拍| 亚洲精品乱码久久久v下载方式| 免费大片18禁| 国产免费福利视频在线观看| 欧美丝袜亚洲另类|