《電子技術應用》
您所在的位置:首頁 > 模擬設計 > 設計應用 > 電力系統低頻振蕩穩定監測分析方法綜述
電力系統低頻振蕩穩定監測分析方法綜述
2015《電子技術應用》智能電網增刊
高 潔1,汪 佳2,高曙光1
(1.西南交通大學 電氣工程學院,四川 成都 610031; 2. 四川省電力公司計量中心,四川 成都 610045)
摘要: 對電力系統低頻振蕩穩定監測分析方法進行了回顧綜述。首先針對系統在大擾動激勵以及環境激勵下不同的測量響應,將低頻振蕩分析方法分為兩大部分;然后根據算法流程將各部分方法進一步分為時域、頻域法以及線性、非線性方法。在分類的基礎上,對部分重點方法的適用條件、基本原理以及應用進行了詳細闡述,探討了進一步研究的方向,并為之提供了有益的思路。
Abstract:
Key words :

  高  潔1,汪  佳2,高曙光1

  (1.西南交通大學 電氣工程學院,四川 成都 610031;2. 四川省電力公司計量中心,四川 成都 610045)

  摘  要: 對電力系統低頻振蕩穩定監測分析方法進行了回顧綜述。首先針對系統在大擾動激勵以及環境激勵下不同的測量響應,將低頻振蕩分析方法分為兩大部分;然后根據算法流程將各部分方法進一步分為時域、頻域法以及線性、非線性方法。在分類的基礎上,對部分重點方法的適用條件、基本原理以及應用進行了詳細闡述,探討了進一步研究的方向,并為之提供了有益的思路。

  關鍵詞: 低頻振蕩;穩定監測;大擾動激勵;環境激勵

0 引言

  低頻振蕩是電力系統中存在的固有現象。隨著電網規模的不斷擴大,運行極限的不斷逼近,弱阻尼振蕩不僅限制了區域電網間的功率傳輸,而且可能影響系統穩定、造成大規模停電事故。因此,低頻振蕩在線監測是保障系統動態穩定性的重要環節。近幾年,隨著廣域測量技術的發展,國內外電力系統中已經有大量投入同步相量測量裝置(PUM),對全網實時功角數據及電網中發生的所有異常工況進行準確的記錄,是電力系統安全穩定運行中有力的技術支持[1]。因此,利用PUM的實測軌跡獲取低頻振蕩的模態信息,不僅避免了人為對系統動態仿真來求取受擾軌跡的不便性,還能很好地反映實際系統的動態信息,并進一步發展為在線監測,具有廣泛的實用前景。由于導致系統隨機波動的輸入激勵無法測量,基于實測輸出數據進行低頻振蕩在線分析成為了研究重點,也是對相關科研工作者的一個重要挑戰。

001.jpg

  低頻振蕩的分析研究類似于振動力學問題的研究。振動力學問題的研究主要集中在結構(工程結構、機械結構等)于外界激勵的作用下產生的受迫振蕩上。由于動靜荷載的長期作用、環境侵蝕等原因,降低了結構抵抗正常荷載的能力,更無法抵抗自然災害,因此,利用模態分析實現結構健康監測顯得尤為重要。但是由于大型工程結構組成材料蕓雜、規模龐大、邊界條件復雜,傳統實驗模態識別技術的局限性開始突顯,首先正常工作環境下難以對結構施加有效的人為激勵;其次局部施加的人為激勵容易對結構造成損壞,因此,大型結構模態分析主要通過測量結構的振動信號提取有關特征參數,并得到其實際工作狀態[2]。而電力系統在正常運行條件下也很難量測到系統的輸入激勵,因此,該領域的研究也主要集中在利用量測輸出信號進行系統識別。由此可以發現,兩個領域的研究具有極大的相似性:結構相當于電力傳輸網絡;結構激勵相當于電力系統勵磁或加載應用;結構振動量測相當于電力系統的功率、電壓、電流以及頻率量測等。因此,有理由設想:對于電力系統低頻振蕩問題可以借鑒振動分析領域的類似方法進行在線分析,而且成功交叉應用的可能性很大。目前,兩個領域的交叉應用研究已經得到相關學者的關注,并取得了一定的進展。

  電力系統的測量響應可分為環境激勵響應以及大擾動激勵響應(瞬態響應),根據測量響應的不同可以將低頻振蕩在線分析方法進行有效的分類[3]。環境激勵響應常常基于正常運行狀態下的系統在運行點附近是線性的假設,負荷變化等產生的自然激勵可以近似認為是隨機高斯白噪聲,而系統的輸出在本質上是隨機的。大擾動激勵響應描述測量系統在遭受大的擾動或故障后的瞬態反映,該響應通常的特點是系統頻率或其他變量產生較大的偏差,例如在輸電線路功率流等。基于環境激勵下的低頻振蕩分析方法對系統模態頻率估計相對于系統阻尼要更容易些,而基于瞬態響應的低頻振蕩分析方法則更側重于描述系統的阻尼振蕩行為[4]。圖1[3]描述了某一電力系統在環境激勵下以及瞬態操作下的不同輸出響應曲線。

  參考文獻[4]針對電力系統區間振蕩信號的有效分析方法進行了總結,根據信號以及方法的不同特性進行了合理分類,并將所給示例的方法歸到不同類別之中。本文采用相同分類的一個樹形圖將基于量測數據的低頻振蕩分析方法進行有效分類,其中包括一些最新應用的方法,如圖2所示。

002.jpg

1 大擾動激勵下模態參數識別方法

  大擾動激勵響應是指系統受到短路、斷線故障等大擾動激勵后的反應。該響應的特點是系統參數將發生較大偏差,例如頻率。假設快速故障或干擾產生后,系統瞬態代表真實系統的脈沖響應,則瞬態分析的目的是通過瞬態振蕩頻率和阻尼估計對系統穩定性進行評估。相關文獻對大擾動激勵下模態參數識別方法的分類常常基于系統瞬態反應過程呈線性的假設。而圖2對該類分析方法則分為線性及非線性。

  1.1 線性分析方法

  線性分析方法假設故障或干擾后系統保持線性,并通過一些衰減正弦曲線的線性組合來描述瞬態響應數據的數學模型。假設系統的量測輸出y(t)是由n個加權的衰減曲線組成(其中幅值為Bi,復頻率為λi),λi可以分解得到角頻率ωi和衰減因子αi,其數學描述如下:

  )P_BG$1XH`0`JP~58T1@D$8.png

  線性分析方法的不同在于采用不同的途徑去獲得該數學模型,同時由圖2可知該類方法又可分為時域分析方法及頻域分析方法。

  1.1.1 時域分析方法

  時域分析方法通過建立時間序列的線性模型,進行數據分析。圖2中列舉了4種不同的時域分析法:Prony方法、矩陣束算法(Matrtx Peneil,MP)、Hankel Total Least squares(HTLS)方法、特征系統實現算法(Eigensystem Realization Algorithm,ERA)。這些方法大多應用于單信號分析,部分算法需要構造Hankel矩陣。

  (1)Prony方法

  Prony方法是利用Prony分析來確定先前所描述的線性模型,其過程分為兩步,首先建立描述等間距采樣數據的線性預測模型,獲得了模型的回歸系數;然后對模型多項式的特征值進行求解,進而得到不同的模態頻率、幅值、衰減因子等信息。參考文獻[5]首次將Prony方法引入電力系統低頻振蕩分析中,并仿真驗證了其有效性。由于Prony方法不需要建立系統數學模型,相較于特征值分析法,對大電網振蕩模式分析具有獨特的優勢[6],是目前應用最為廣泛的一種低頻振蕩分析方法,并且基于Prony辨識結果設計出阻尼控制器[7-8]也得到了廣泛的應用。

  但是Prony方法辨識精度受模型階數等因素的影響很大且抗噪性能較差。針對噪聲干擾,參考文獻[9-12]分別提出結合低通濾波器、模糊濾波、小波變換、經驗模態分解數據預處理方法,對測量信號中的噪聲及高頻成分進行有效抑制與濾除,從而改善了Prony方法的辨識效果。針對模型階數選取問題,參考文獻[9]提出利用自回歸模型構造采樣數據的Hankel矩陣,采用不同階次下矩陣的行列式比估計系統模型的最佳降維階數。參考文獻[13]利用二階樣本矩陣的奇異值分布特征來估計系統的實際階數。參考文獻[14]對Prony方法分析低頻振蕩的有效性進行了驗證,討論了信號噪聲及非平穩特性對Prony方法的影響,提出計算均方差確定算法階數以及縮小數據分析窗口避免信號非平穩性的影響等方法,對Prony算法進行了改進。

  (2)矩陣束算法

  MP算法[15-16]需要構造兩個束函數Y1、Y2,根據其總體特征值等價于函數一般特征值的原理,對系統特征根進行求取。該方法首先對測量輸出矩陣進行奇異值分解(SVD),矩陣行與列的長度取決于數據長度以及矩陣束參數。SVD分解后,篩選出n個奇異值,并由此數據截斷得到新的奇異矩陣,通過反向計算構造束矩陣Y1、Y2,因此計算出系統特征值。MP方法將矩陣維數進行了簡化,同時具有一定的抗噪能力。

  (3)特征系統實現算法

  特征系統實現算法[17]是1984年由JUANG J N等人提出的一種多輸入、多輸出的時域整體模態參數識別方法,該方法以多點激勵得到的脈沖響應函數為基礎,構造 Hankel 矩陣,構成最小階的系統實現,并將該實現變換為特征值規范型。最初以脈沖響應函數為輸入數據在航空航天復雜結構中得到廣泛應用,后來聯合隨機減量技術 (random decrement,RD)[18]或自然激勵技術(Natural Excitation Technique,NExT)[19]后,可進行環境激勵下的模態識別,并把研究成果應用于土木、橋梁等復雜結構中[20-21]。參考文獻[22]最早提出將ERA引入電力系統領域,通過對低頻振蕩瞬態穩定響應的分析,驗證其有效性,同時具有一定的抗噪能力,而且該方法對延時數據以及基本模式的辨識也同樣有效。參考文獻[23]利用標準的電力系統暫態穩定模型對Steiglitz-McBride、ERA以及Prony 三種算法的模態識能力進行了比較,發現Prony算法與ERA算法對單輸入單輸出系統的識別結果相似,且優于Steiglitz-McBride算法。

  現有研究表明,特征系統實現算法理論推導嚴密,只需很短的數據就可以精確識別系統的模態參數,減小了解算量,同時具有很好的抗噪效果,是目前公認的最完善、先進的模態參數識別方法之一。然而ERA在電力系統領域的應用研究不夠深入,例如其適用性和穩定性,還有待進一步探討;同時與NExT技術相結合,可應用于低頻振蕩環境激勵響應下的模態識別。

  1.1.2 頻域分析法

  大擾動激勵響應中的線性分析方法的第二種分類是頻域分析法。這類方法主要通過對響應信號的頻譜分析得到數據模型參數。參考文獻[24]提出通過滑動窗口的FFT變換得到模態的相對幅值,進而求出這些幅值所對應的模態阻尼。參考文獻[25]通過分析滑窗前后相應譜分量的變化,識別出模式的衰減特性進而求解低頻振蕩的模態參數。該方法原理簡單,識別迅速,對噪聲的魯棒性很好,但無法識別密集模態,對阻尼比識別誤差較大。

  1.2 非線性分析方法

  大擾動激勵下的分析方法的第二類是非線性分析方法,如圖2所示。該類方法認為電力系統其本質是非線性的,并假設小擾動或故障后的激發系統響應也是非線性的。在這種情況下,非線性主要是體現為瞬態響應中頻率成分的交互干擾。其中研究最多的是小波變換算法以及Hilbert-Huang transform (HHT)算法。

  (1)小波變換

  小波分析是一種分析非平穩信號的有效工具,與短時傅里葉變換的恒定分辨率不同,在時域與頻域都具有良好的局部分辨能力,在低頻處達到頻率細分,高頻處達到時間細分,能廣泛適用于時頻信號分析的要求。已經在信號處理、圖像處理、故障診斷、損傷識別等非線性領域得到了廣泛的應用,參考文獻[27]研究基于Morlet小波提取時頻平面上的小波脊點,計算該點的小波變換系數得到低頻振蕩模式。小波分析能夠很好地反映振蕩頻率的時變性,也具有一定抗噪性,但實際應用中小波脊特性較差,當噪聲較大或模態密集時難以提取。參考文獻[28]提出將小波變換和SVD相結合,用提升小波系數SVD的頻率向量來識別各階振蕩模式的頻率,提高算法的抗噪性能。

  (2)HHT

  Hilbert-Huang變換(HHT)[29]是1998 年由HUANG N E提出的一種能自適應處理非線性、非平穩信號的新方法。該方法的基本過程是先將時間信號進行 EMD 分解,產生一組具有不同特征時間尺度的 IMF,實現各模態分量的有效分離,然后運用Hilbert 變換得到各分量的瞬時振幅以及瞬時頻率等。HHT在電力系統低頻振蕩分析領域已經展開了研究應用[30-31],但HHT方法存在篩選次數難以確定、模式混疊以及Hilbert譜分析存在局限等諸多問題,對低頻振蕩模式參數的辨識精度將造成一定的影響。

2 環境激勵下模態參數識別方法

  電力系統環境激勵是指負荷投切、風載、沖擊波等隨機性質的小擾動。在日常運行過程中,即使是正常運行狀態,電力系統也時刻存在環境激勵的影響,其響應振幅小,易于采集,涵蓋的頻率豐富,可以及時準確地反映當前電力系統的運行特性。由于系統發生短路、斷線等大擾動的情況不是時刻發生,因此,利用環境激勵下響應數據辨識系統的振蕩信息,可以對大擾動下振蕩特性辨識進行有益補充,同時也可以用于電力系統的在線監測與狀態評估,具有重要的實際意義和工程應用價值。

  由于環境激勵是一種未知且不可控的激勵源,更難以運用合適的數學公式來表達。結合力學公式可知,“輸入—系統—輸出”中的前兩項均為未知量,這給理論與實際應用中的選用標準問題帶來了新的挑戰,為此國內外相關學者做出了許多突破性的研究[32],尤其是在土木建筑、航空航天、汽車、船舶制造等領域 [33-35]。近年來,環境激勵模態識別方法在電力系統低頻振蕩分析中得到越來越廣泛的重視。基于前面的設想,電力系統低頻振蕩可以借鑒振動分析領域的類似方法進行環境激勵下的模態分析研究,按照識別信號域的不同主要可分為頻域方法與時域方法兩種。

  2.1 頻域分析方法

  環境激勵的模態分析中,頻域法大多利用功率譜密度函數進行估計,較傳統的計算方法有快速傅里葉變換(FFT)[36]和韋爾奇周期圖方法(Welch periodogram methods)[37],它們均屬于非參數化的方法,利用傅氏變換將時域信號轉入頻域,從而得到一系列不同的振動分量信息。

  高階譜(Higher Order Spectral, HOS)方法是處理非最小相位系統和非高斯信號的主要分析工具,可以保存不同頻率間的幅值和相位信息,描述二次相位耦合,主要用于狀態監測以及故障診斷。該方法主要包括雙頻譜、雙相干譜以及三頻譜分析。參考文獻[38-39]提出基于HOS的相位分析方法對感應電機的故障進行診斷,并取得很好的效果;參考文獻[40]利用故障信號相耦合的特性,通過雙頻譜渦輪葉片進行狀態監測;參考文獻[41]將高階譜理論中的非參數直接法應用于電力系統低頻振蕩分析,通過雙譜分析及雙相干系數辨識模式間的二次相位耦合信息,揭示模式間的非線性相關作用以及系統的動態行為。

  最小二乘(LS)算法是最常用的遞歸分析方法,通過最小化誤差的平方和從數學模型尋找數據的最佳函數匹配。在此算法基礎上進行優化后得到的最小均方算法(Least-Mean Squares, LMS)[42]、魯棒遞歸最小二乘法(Robust Recursive Least Square, RRLS)算法[43]及卡爾曼濾波技術[44]等,均可作為自適應濾波算法應用于電力系統機電振蕩模態估計中。參考文獻[43]將RRLS算法應用于環境激勵及大擾動激勵下的系統模態識別,通過不同的仿真數據驗證了該方法的適用性,針對17機電力系統模型結合蒙泰卡羅方法,對RRLS方法、傳統的RLS方法以及LMS方法進行比較研究,發現RRLS對模態頻率估計的能力要更勝一籌。

  2.2 時域分析方法

  (1)時間序列法

  時間序列法是一種利用線性輸入輸出模型對固有響應數據進行描述的一種方法。參考文獻[4]提出最早的非遞歸方法是自回歸(AR)模型估計中的Yule-Walker(YW)方程法,該方法經過進一步修正被用于自回歸滑動平均(ARMA)模型參數估計中。早在1997年,PIERRE J W等人將AR模型用于電力系統環境激勵下的類噪聲響應信號辨識,并與大擾動激勵下響應信號的Prony辨識結果相比較,驗證其有效性[45]。參考文獻[46]提出將ARMA模型用于環境激勵下響應數據的模態辨識,并與基于AR模型的辨識結果進行比較,發現兩種方法的頻率估計受模型行列數的影響,而阻尼估計受模型階數的影響,而且隨著階數越高,阻尼辨識結果與Prony的辨識結果越靠近。基于AR模型和ARMA模型的識別方法適用于白噪聲激勵,分辨率較高,可用于在線模態分析,但存在實際應用中模型定階困難等缺陷。為了辨識出所有的模態與噪聲,模型階數的選取要足夠大,但又不能因太大而導致計算效率低下。

  (2)隨機子空間

  隨機子空間法(SSI)[47]由PEETERS B等人于1995 年提出,并首次應用于大型結構參數識別。該算法首先由輸出響應數據構造Hankel 矩陣或者Toeplitz矩陣,通過矩陣截斷構造線性子空間,降低矩陣維數,并利用截斷頻率對系統的模態參數進行估計。根據辨識方法性質的不同,隨機子空間識別法可進一步分為基于協方差驅動隨機子空間識別法(Cov-SSI)和基于數據驅動隨機子空間識別法(Data-SSI)兩種[26],其中Data-SSI法將Hankel 矩陣進行 QR 分解,求得投影矩陣后進行SVD 分解,獲得卡爾曼濾波狀態向量,在狀態確定的情況下將識別問題轉變為系統矩陣的線性最小二乘問題得到系統失穩模態參數;與Cov-SSI法相比,避免了多次協方差矩陣的計算,計算效率更高,更適用于自動識別。

  隨機子空間法基于系統處于平穩隨機信號激勵下的假設,適用于環境激勵條件下低頻振蕩模態參數的識別[48-50]。該方法與其他頻域識別方法相比,識別過程可直接作用于時域數據,無需傅氏變換等處理,沒有頻率分辨率誤差的問題,非常適合識別接近空間模態的系統;該方法不但能準確識別系統的頻率,而且能識別出系統的模態振型和阻尼,在土木結構、航空航天、大型機械結構等領域都有成功應用。國內相關學者將該方法引入低頻振蕩大擾動激勵下模態參數識別中[51-52],但在環境激勵條件下低頻振蕩模態參數識別的應用并未深入。

  以上兩種模態識別的時域方法均可直接對環境激勵下的低頻振蕩類噪聲信號進行模態參數求取,然而這種一步求取的時域方法種類畢竟有限,而且許多傳統的時域模態識別方法無法直接對類噪聲信號進行有效辨識。因此可以通過信號預處理,得到與脈沖響應函數近似的互相關函數、自由振蕩響應等,之后再將時域經典模態識別算法進行擴展運用。下面將介紹兩種預處理算法:隨機減量法和自然激勵技術。

  2.3 預處理算法

  (1)隨機減量法

  隨機減量法是利用樣本平均的方法,去掉響應中的隨機成分,從而獲得一定初始激勵下的自由響應[53]。由于很多傳統的實驗模態識別方法要求系統的響應數據為自由振動響應數據,故可以先用隨機減量法將脈動時域響應處理為自由振動響應,然后根據自由響應的數學模型建立特征方程,求出特征根后再估算各階模態參數。

  該方法已成功用于多個工程結構的模態參數識別工作[54-55],參考文獻[56]首先采用隨機減量法從環境激勵下低頻振蕩類噪聲信號中提取自由衰減響應數據,然后采用Prony算法對該數據進行模態參數辨識,其結果與實測擾動后系統響應的分析結果一致。但該方法僅在理論上適合白噪聲激勵,且所用的是單通道信號,存在模態丟失的現象。

  (2)自然激勵技術

  自然激勵技術(Natural Excitation Technique, NExT)[57]是一種環境激勵條件下利用互相關函數近似獲得脈沖響應的有效方法,擴展了傳統以脈沖響應函數進行模態識別的方法應用。該方法的基本思想是:在白噪聲激勵條件下,線性系統兩個響應點之間的互相關函數與任意一點的脈沖響應函數具有近似的解析表達形式。NExT法利用相關函數作為識別計算輸入,具有一定的抗噪能力,但在辨識參數時沒有自己的計算公式,需要借助傳統的模態分析方法。NExT 技術幾乎在所有和結構動態分析有關領域得到廣泛應用[58-60],它避免了對橋梁、船舶、高層建筑、運載火箭等大型結構進行人為激勵的不便,僅利用大地脈動、車輛、風等自然激勵進行大型結構工程的模態參數識別,簡單快捷,并且使識別結果更加符合實際情況。在電力系統領域,參考文獻[61]中采用NExT與HHT相結合作用于電力系統低頻振蕩類噪聲信號的辨識,除此,NExT技術在該領域的應用研究并未深入。

3 結論

  電力系統相關研究人員所面臨的一個重要挑戰是輸入激勵未知條件下的系統識別,這更是電網正常運行監測的一個重要挑戰。隨著電力系統規模擴大并接近其運行極限,區域間的弱阻尼將產生低頻振蕩甚至導致系統失穩,因此,電網的實時監測顯得尤為必要。在正常運行期間,負載變化將引起的功率轉移和其他系統變量的隨機波動,而全網負載的波動要全部進行測量是有難度的,同時使用標準的輸入輸出模型進行分析是有局限性的,因此通過實際振動測量信號進行分析更具有實際意義以及工程應用價值。

  本文根據測量響應的不同,將低頻振蕩穩定監測分析方法進行了分類與簡單介紹,這些方法均基于系統輸入激勵未知的假設,且多為隨機高斯白噪聲。目前,基于環境激勵下的低頻振蕩分析存在激勵不平穩、模態識別精度較差、虛假模態較多等問題,其研究也有待進一步深入。由于結構工程振蕩分析研究已經數十年,尤其是環境激勵下的模態參數識別方法的研究已經比較成熟,因此可被應用到電力系統低頻振蕩研究中,令其穩定監測分析更加深入。

參考文獻

  [1]陳實,許勇,王正風,等. 電網實時動態監測技術及應用[M]. 北京:中國電力出版社,2010.

  [2]祁泉泉. 基于振動信號的結構參數識別系統方法研究[D]. 北京:清華大學,2011.

  [3]THAMBIRAJAH J,BAROCIO E,Thornhill N F. Comparative review of methods for stability monitoring in electrical power systems and vibrating structures[J]. IET Gener. Transm. Distrib., 2010, 4. (10): 1086-1103.

  [4]MESSINA A R,TRUDNOWSKI D,PIERRE J.Inter-area oscillations in power systems - a nonlinear and nonstationary perspective[C].Springer,2009:  1-36.

  [5]HAUER J F, DEMEURE C J, SCHARF L L. Initial Results in Prony Analysis of Power System Response Signals[J]. IEEE Trans.  Power Systems, 1990. (5):80-89.

  [6]GRUND C E, PASERBA J J, HAUER J F, et al. Comparison of Prony and eigen analysis for power system control design[J]. IEEE Transactions on Power Systems, 1993, 8(3): 964-971.

  [7]TRUDNOWSKI D J, SMITH J R, SHORT T A, et al.An application of prony methods in PSS design for multimachine systems[J]. IEEE Trans. Power Syst., 1991, 6(1): 118-126.

  [8]蘆晶晶,郭劍,田芳,等. 基于Prony方法的電力系統振蕩模式分析及PSS參數設計[J].電網技術,2004,25(15):31-34.

  [9]鞠平,謝歡,孟遠景,等.基于廣域測量信息在線辨識低頻振蕩[J].中國電機工程學報,2005,25(22):56-60.

  [10]李大虎,曹一家.基于模糊濾波和Prony算法的低頻振蕩模式在線辨識方法[J].電力系統自動化,2007,31(1):14-1.


此內容為AET網站原創,未經授權禁止轉載。
主站蜘蛛池模板: 狠狠色噜噜狠狠狠狠色综合网 | 亚洲爆乳无码一区二区三区 | 国产96在线 | 黄色国产大片 | haose08永久免费视频 | 久久成人免费观看全部免费 | 中文天堂网 | 国内午夜免费鲁丝片 | 天天澡夜夜澡狠狠澡 | 久久久精品免费视频 | 制服女子校生在线调教 | 久久久噜噜噜久久网 | 久久久久久88色愉愉 | 欧美中文字幕无线码视频 | 黄网站在线观看视频 | 操美女在线 | 欧美日在线 | 久久综合丝袜长腿丝袜 | 天天色综合3 | 欧美日韩一区二区三区视频 | 国产无遮挡床戏视频免费 | 青春草在线观看精品免费视频 | 免费二级c片观看 | 天天操人人射 | 天堂中文在线免费观看 | 国产高清一级毛片在线人 | 国产免费播放一区二区 | 看免费毛片天天看 | 综合五月激情 | 手机在线精品视频 | 亚洲羞羞视频 | 国产噜噜噜精品免费 | 欧美xxx性 | 两性午夜欧美高清做性 | 国产丝袜大长腿精品丝袜美女 | 天天插天天狠天天透 | 欧美精品在线观看视频 | 最新国产精品精品视频 | 香蕉官网 | 国产高清视频网站 | 欧美一区二区另类有声小说 |