999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于正交變換與置信域的量測方差估計與權重設置算法

2016-12-27 06:04:56盧志剛田莎莎
電工技術學報 2016年21期

盧志剛 田莎莎,2 邵 奇,3 吳 杰

(1.燕山大學電力電子節能及傳動控制河北省重點實驗室 秦皇島 066004 2.國網滄州供電公司 滄州 061000 3.國網天津市電力公司物資公司 天津 300304)

?

基于正交變換與置信域的量測方差估計與權重設置算法

盧志剛1田莎莎1,2邵 奇1,3吳 杰1

(1.燕山大學電力電子節能及傳動控制河北省重點實驗室 秦皇島 066004 2.國網滄州供電公司 滄州 061000 3.國網天津市電力公司物資公司 天津 300304)

在狀態估計實際應用中,量測方差獲取和權重設置存在一定的困難。伴隨狀態估計運算量越來越繁重,現有量測方差估計算法的收斂性無法得到保證。為此提出了一種基于正交變換與置信域的量測方差估計和權重設置算法。利用正交變換降低迭代系數矩陣條件數,提高量測方差估計的數值穩定性;根據正交變換的化簡結果,結合拉格朗日乘子法,建立新的量測方差估計模型;對于迭代過程的中間數據,采用置信域作為目標函數的約束,確保每次估計結果落在置信區間內;最后,通過IEEE 14/30/118標準系統算例證明了該算法的有效性。

狀態估計 正交變換 置信域 量測方差 權重 收斂性

0 引言

狀態估計是能量管理系統的重要組成部分并得到廣泛應用。但是在實際中仍然存在值得研究的問題,例如狀態估計的合格率與精度[1]。

電力系統利用最小二乘法對冗余量測量擬合得到狀態量的最優估計值。由監測設備采集并經由通信系統到達狀態估計器的量測量存在不同程度的誤差。對于量測量精度不相等的問題,可以通過加權來提高估計的準確性[2-5]。根據最小二乘原理,權重取量測方差的倒數。權重與量測方差不僅是保證估計精度的重要前提,還是不良數據辨識與抗差狀態估計[6]的基礎參數,應重視并合理設置這兩個參數。

文獻[7]探討了確定權重的原則,指出數值穩定性與估計準確性需要同時兼顧。具體的設置方法目前主要有組合參數法[4,8]、量測分類設置法[9]和樣本估計法[10-12]。第一種方法應用較早,主要利用儀表精度、轉換誤差等設備參數組合設置權重。該方法沒有考慮通信噪聲的影響,也沒有考慮裝置逐年老化產生的誤差偏移。即便如此,作為遠端的狀態估計維護人員也很難獲取在現場運行的裝置精度參數。第二種方法將量測分為注入、支路等類型,綜合層次分析等模糊評價方法確定權重。第三種方法通過量測殘差與誤差的數學模型對樣本數據計算得到量測方差和權重。文獻[10,11]根據殘差靈敏度矩陣搭建殘差與誤差模型。但殘差靈敏度矩陣不滿秩限制了該法的應用。在此基礎上,文獻[12]簡化了殘差與誤差的數學模型,建立了自適應權重更新算法,即迭代重加權量測方差估計算法(Iterative Re-weighted Measurement Variance Estimation,IRMVE)。狀態估計求解運用牛頓法及正則方程組(Normal Equations)法。牛頓法容易發生數值病態問題,現場運行有時會發生不收斂現象[13]。IRMVE法迭代計算量大,不收斂現象較為嚴重,限制了應用范圍。

針對狀態估計數值病態問題,等式約束法[14,15]、正交變換法[16-19]、置信域法[20-22]等給出了解決方案。等式約束法針對電力系統零注入節點較多的情況,將虛擬量測處理成等式約束,主要克服了虛擬量測權重系數設置過大導致的數值問題。正交變換法直接降低迭代過程中系數矩陣的條件數,顯著提高數值穩定性。而置信域法將狀態空間處理成不等式約束,在置信域內,可使二次模型更好地適應非線性目標函數,是提高狀態估計收斂性的重要手段。

本文首先研究現有量測方差算法,針對收斂性以及數學模型不滿秩等問題做出了改進。在采用正交變換提高迭代數值穩定性的同時,推導出新的誤差方差計算模型,并利用狀態置信域保障單次迭代計算結果的可信性。最后通過IEEE 14/30/118節點標準系統算例進行驗證。

1 迭代重加權算法

該算法主要由以下3個步驟構成:

1)確定量測樣本總數,量測樣本通過相同時段不同時刻的量測數據構成。對每個時刻樣本數據進行狀態估計獲得單個時刻樣本的量測殘差,量測模型為

z=h(x)+v

(1)

式中,z、v分別為m階的量測向量與量測誤差向量;狀態向量x含有n=2N-1個元素;N為系統的節點數。利用加權最小二乘原理,目標函數為

J(Δx)=[Δz-HΔx]TW[Δz-HΔx]

(2)

建立正則方程組求解

Δx=(HTWH)-1HTWΔz

(3)

當全部樣本都經過狀態估計的計算獲得了量測殘差,則稱之為實現了一次“全計算”。

2)計算量測方差。運用加權殘差與加權誤差方程

rw=Swvw

(4)

Sw=I-W1/2H(HTWH)-1HTW1/2

(5)

式中,rw、 vw、 Sw分別為m維加權殘差向量、m維加權誤差向量和m階的加權殘差靈敏度矩陣。式(4)建立起殘差與誤差的關系。然而,加權殘差靈敏度矩陣Sw的秩為m-n, 不能夠用求逆的方法求解式(4)。這給求量測方差帶來一定的困難。對此有簡化的誤差協方差的模型[10]如式(6)所示。

Cov(r)=SRvST=SRv

(6)

式中,Rv=E[vvT]; Cov(r)的對角元素Rii與Rv的對角元素Rvii可通過S對角元Sii建立如下關系:

Rii=SiiRvii

(7)

殘差靈敏度矩陣對角元Sii在理論上是非負數,但是在數值計算過程中,該對角元如果越來越小,會發生“變號”形成負值,導致計算結果迅速變壞。

3)根據新的權系數矩陣,判斷是否收斂。若收斂,則停止輸出最終的量測方差和權系數。若不收斂,則回到步驟1),重復上述過程,直到收斂為止。這就是迭代重加權量測方差估計算法(IRMVE)。

IRMVE運算并不穩定,不收斂情況會在前幾次全計算中出現。直接原因在于權系數的差異會隨全計算次數增多而擴大,同時殘差靈敏度對角元發生“變號”。根本原因在于正則方程組法本身容易發生數值病態,這限制了IRMVE向更大規模節點系統應用。

2 基于正交變換與置信域的量測方差估計與權重設置方法

2.1 目標函數的處理

針對正則方程組法給狀態估計帶來的數值問題,正交變換法能夠使迭代過程中的系數矩陣的條件數降低,從根本上提高數值穩定性。簡要分析正交變換原理,其形式為

(8)

Givens法是一種重要的正交變換方法。Givens法通過挖掘雅克比矩陣的稀疏條件,可以顯著提高正交變換的速度。IEEE 118節點系統單次狀態估計能夠控制到幾十ms范圍內[16]。

IRMVE需要數以千次甚至萬次的迭代運算,中間結果不可能一一輸出用以辨識其準確性。對此,可以給出狀態的某個值得信賴的區域作為約束,從而建立全局收斂的狀態估計模型,保證計算結果能夠分布在可信區域內。并且,可以將正交變換與置信域相結合。

首先給出某一狀態初值,在迭代過程中,如果狀態增量滿足二次準則

(9)

式中,q(·)為二次準則目標函數;J(·)為原目標函數;δ為狀態增量的步長;ρ為置信域半徑。通過將步長因子參數化,使式(9)變為等式約束問題。

(10)

忽略式(10)的二次項,并應用拉格朗日乘子法,得到[21]

(11)

考慮正則方程組式(3),有

(HTR-1H+μkI)δ(μk)=HTR-1[z-h(xk)]

(12)

將式(12)進一步分解為

(13)

(14)

2.2 量測方差計算模型

本文引入正交變換,一方面是為了提高數值穩定性;另一方面,借正交變換對殘差靈敏度矩陣進行化簡,給出更明確的量測方差計算模型。

將式(11)代入式(4),推導過程如下:

得到

(15)

式中,On為n階全零方陣;Im-n為m-n階單位陣。用QT左乘式(15)兩邊得

(16)

(17)

整理得

(18)

(19)

(20)

基礎量測誤差可以看做是一組確定的值。根據假設:量測已經通過剔除顯著粗差,粗差水平并不高。從而建立如下求解量測方差模型。

(21)

根據拉格朗日乘子法得到目標函數

(22)

式中,λ為m-n階系數向量。

可得到極值條件

(23)

每個樣本按照正交變換與置信域法進行狀態估計,利用上述模型計算誤差,通過多樣本誤差獲得量測方差,進而得到各個量測的權系數,計算方法如下:

(24)

(25)

式中,i為量測量序號;K為樣本總數。考慮到K個樣本誤差自由度為K, 所以式(24)除數為K而非K-1。 樣本個數并非越多越好,當樣本個數達到一定程度時,計算精度就不再提高。為提高計算效率,應該合理設置樣本容量。最終得到了基于正交變換與置信域的量測方差估計算法(Orthogonal Transformation and Trust Region based Measurement Variance Estimation,OTTRMVE)。

2.3 OTTRMVE算法流程

準備工作:確定最小樣本容量值K。 利用同一格式存儲相應個數的電力系統運行數據與網絡結構參數。

1)利用已知信息初始化方差矩陣和權系數矩陣W。 如果方差及權系數未知,則將權系數統一為1。輸入其他初始化參數。

2)時刻計數器k=1, 全計算次數計數器t=1。

3)根據式(14)的正交變換與置信域法狀態估計,對當前樣本數據進行計算。判斷狀態估計是否收斂,收斂則跳轉至步驟4);不收斂繼續步驟3)。如果狀態估計迭代次數達到上限,則結束。

4)記錄k時刻的量測殘差。判斷時刻計數器k是否小于樣本容量值。若小于,則計數器k++, 返回至步驟3);如果等于,跳轉至步驟5)。

5)根據式(24)、式(25)計算全計算次數t水平下的權系數。如果t=1, 則t++, 并返回至步驟3);如果t>1, 判斷權系數是否收斂。如果收斂則結束,并輸出權系數與誤差方差向量;如果不收斂,則跳轉至步驟6)。

6)判斷t是否到達上限。如果小于上限,則t++, 返回步驟3);如果達到上限,則結束。流程如圖1所示。

圖1 OTTRMVE算法流程Fig.1 The flowchart of the OTTRMVE algorithm

3 算例仿真

3.1 仿真環境與算例設置

仿真軟件為Java(Eclipse),虛擬機(JVM)最大內存512M。計算機為雙核2.0GHz便攜式計算機。本文采用IEEE 14/30/118節點三個標準系統算例,各個系統基本的信息見表1。

表1 測試系統與量測配置信息Tab.1 Test systems and measurement configurations

采用“全量測”的配置方式,即已知所有電壓幅值、注入有功無功、支路一端的有功無功。這主要是避免對杠桿量測與關鍵量測問題的討論。樣本量測數據模擬負荷變化緩慢情況下時間上足夠接近的系統運行數據。量測編號為IEEE默認編號。

3.2 IEEE 14節點系統

量測標準差設置如下:電壓幅值:0.004,除母線3、8(0.080)與母線5、9(0.001);注入功率:0.01,除母線4(0.05)與母線7、10(0.001);支路功率:0.008,除支路3、5、9(0.1)與支路10(0.001)。OTTRMVE與IRMVE法收斂對比如圖2、圖3所示。

圖2 IEEE 14節點OTTRMVE收斂判據變化Fig.2 The convergence criterion variation of OTTRMVE for IEEE 14-bus system

圖3 IEEE 14節點IRMVE收斂判據變化Fig.3 The convergence criterion variation of IRMVE for IEEE 14-bus system

對于相同的收斂精度,當全計算次數達到第3次時,OTTRMVE法即收斂,收斂判據不再發生變化,而IRMVE法則呈現出發散趨勢。值得指出的是,這樣的收斂精度是較為寬松的,特別是對電壓幅值標幺值而言。

IEEE 14節點OTTRMVE法數值估計精度如圖4所示。量測誤差標準差真值設置的部分突變點,可以反映OTTRMVE法對于突變值的計算精度。從圖4可以看出,OTTRMVE法能夠較為準確地反映誤差方差的突變。OTTRMVE法保證每一步狀態估計的結果收斂到狀態置信域當中。在突變幅度較大的情況下仍能較好地反映誤差方差的變化。

圖4 IEEE 14節點OTTRMVE法估計精度Fig.4 The estimation accuracy of OTTRMVE for IEEE 14 -bus system

由表2可見,隨著樣本數量的增加,OTTRMVE法估計值的相對誤差方差呈現出先減小而后增大的趨勢。根據前文分析,并非樣本數量越多估計精度越高。對于一定的估計精度要求,樣本數量存在一個最小值。由表2可知,該值為100左右。合理規劃樣本個數,可以有效縮短估計時間,提高效率。

表2 IEEE 14節點系統OTTRMVE試驗結果Tab.2 Test results of OTTRMVE for IEEE 14-bus system

3.3 IEEE 30節點系統

IEEE 30節點系統量測規模明顯增大,系數矩陣規模達到214維。用以考察中等維度水平下OTTRMVE法的計算精度,如圖5所示。

圖5 IEEE 30節點OTTRMVE法估計精度Fig.5 The estimation accuracy of OTTRMVE for IEEE 30-bus system

量測標準差真值設置按照注入有功、注入無功、支路有功、支路無功、電壓幅值排序。每一類量測的標準差真值設置彼此不同。事實上,量測序號可以任意排序,這對估計結果沒有影響。從圖5可以看出,絕大部分誤差標準差估計值與真值相吻合,結果精度令人滿意。

3.4 IEEE 118節點系統

IEEE 118節點系統雅克比矩陣階數為712,狀態估計迭代計算量顯著增加,用以考察較高維度水平下IRMVE與OTTRMVE法的收斂性能。IRMVE和OTTRMVE全計算結果見表3。

表3 IEEE 118節點系統方差估計結果對比Tab.3 The comparison of variance estimate results for IEEE 118-bus system

全計算第3次OTTRMVE法收斂,而IRMVE法彈出錯誤并結束。從表3看出,截止到前3次全計算,兩種方法耗時大致相同,OTTRMVE稍快一些。這是因為系統節點數目越多,雅克比矩陣稀疏程度越高。與小系統相比,大系統更能發揮Givens正交變換法的效率優勢。

盡管最大全計算次數是10次,但是IRMVE法第3次權系數就已不收斂,而OTTRMVE法在第3次收斂。

OTTRMVE法的最大最小方差在合理區間,其最大最小方差之比要遠小于IRMVE法。IRMVE的“最小方差”是一個負值,這是殘差靈敏度矩陣對角元出現負值造成的,與表3所示負值相對應的殘差靈敏度矩陣對角元為-3.29×10-7。負值是IRMVE數值不穩定的一個顯著現象。與之形成對比的是,OTTRMVE法始終未出現負值,顯示良好的穩定性。

4 結論

本文針對現有方差估計算法在收斂性上存在的問題,提出了基于正交變換與置信域的量測方差估計與權重設置算法。該方法采用正交變換法提高單步迭代的數值穩定性。通過化簡殘差靈敏度矩陣,得到一種新的量測誤差計算模型。同時,將置信域作為目標函數的約束條件可確保計算結果在合理區間。仿真結果證明了該算法對于大規模量測系統仍然具有較快的處理速度,可以為狀態估計初始化權重以及確定量測方差提供一種解決思路。下一步將對杠桿量測與關鍵量測對算法的影響進行研究。

[1] 顧錦汶.對有關標準中關于狀態估計的一些不同意見[J].電力系統自動化,2014,38(2):134-135. Gu Jinwen.Disagreements discussion on state estimation in related standards[J].Automation of Electric Power Systems,2014,38(2):134-135.

[2] Schweppe F C,Wildes J.Power system static-state estimation,part I:exact model[J].IEEE Transactions on Power Apparatus and Systems,1970,89(1):120-125.

[3] Robert E L,William F T.State estimation in power systems,part I:theory and feasibility[J].IEEE Transactions on Power Apparatus and Systems,1970,89(3):345-352.

[4] 于爾鏗.電力系統狀態估計[M].北京:水利電力出版社,1985:54-60.

[5] Abur A,Expósito A G.Power system state estimation,theory and implementation[M].New York:Marcel Dekker,2004:20-58.

[6] 李靜,羅雅迪,趙昆,等.考慮大規模風電接入的快速抗差狀態估計研究[J].電力系統保護與控制.2014,42(22):113-118. Li Jing,Luo Yadi,Zhao Kun,et al.Research of fast and robust state estimation considering large-scale wind power intergration[J].Power System Protection and Control,2014,42(22):113-118.

[7] 李碧君,薛禹勝,顧錦汶,等.狀態估計中選取量測權值的新原則[J].電力系統自動化,2000,24(8):10-14. Li Bijun,Xue Yusheng,Gu Jinwen,et al.A new criterion of determining measurement weights in power system state estimation[J].Automation of Electric Power Systems,2000,24(8):10-14.

[8] Vanslyck L S,Allemong J J.Operating experience with the AEP state estimator[J].IEEE Transactions on Power Systems,1988,3(2):521-526.

[9] 黃知超,謝霞,王斌.結合模糊綜合評判與決策的電力系統狀態估計[J].電力系統保護與控制,2015,43(7):65-69. Huang Zhichao,Xie Xia,Wang Bin.Power system state estimation combined with fuzzy comprehensive evaluation and decision-making[J].Power System Protection and Control,2015,43(7):65-69.

[10]劉廣一,于爾鏗,夏祖治.量測系統誤差方差的估計與修正[J].中國電機工程學報,1990,10(6):31-39. Liu Guangyi,Yu Erkeng,Xia Zuzhi.Estimation and update of measurement system error variance[J].Proceedings of the CSEE,1990,10(6):31-39.

[11]Liu Guangyi,Yu Erkeng,Song Y H.Noval algorithms to estimate and adaptively update measurement error variance using power system state estimation results[J].Electric Power Systems Research,1998,47:57-64.

[12]Shan Zhong,Ali A.Auto tuning of measurement weights in WLS state estimation[J].IEEE Transactions on Power System,2004,19(4):2006-2013.

[13]鄭偉業,吳文傳,張伯明,等.基于內點法的交直流混聯系統抗差狀態估計[J].電力系統保護與控制,2014,42(21):1-7. Zheng Weiye,Wu Wenchuan,Zhang Boming,et al.Robust state estimatior for AC/DC hybrid power system based on an interior point method[J].Power System Protection and Control,2014,42(21):1-7.

[14]郭燁,張伯明,吳文傳,等.直角坐標下含零注入約束的電力系統狀態估計修正牛頓法[J].中國電機工程學報,2012,32(19):96-100. Guo Ye,Zhang Boming,Wu Wenchuan,et al.Power system state estimation solution with zero injection constraints using modified Newton method[J].Proceedings of the CSEE,2012,32(19):96-100.

[15]郭燁,吳文傳,張伯明,等.極坐標下含零注入約束的電力系統狀態估計的修正牛頓法與快速解耦估計[J].中國電機工程學報,2012,32(22):113-117. Guo Ye,Wu Wenchuan,Zhang Boming,et al.Power system state estimation solution with zero injection constraints using modified Newton method and fast decoupled method in polar coordinate[J].Proceedings of the CSEE,2012,32(22):113-117.

[16]Pandian A,Parthasarathy K,Soman S A.Towards faster givens rotations based power system state estimator[J].IEEE Transactions on Power Systems,1999,14(3):837-843.

[17]Pires R C,Costa A S,Mili L.Iteratively reweighted least-squares state estimation trough givens rotations[J].IEEE Transactions on Power Systems,1999,14(4):1499-1506.

[18]杜正春,牛振勇,方萬良.基于分塊QR分解的一種狀態估計算法[J].中國電機工程學報,2003,23(8):51-55. Du Zhengchun,Niu Zhenyong,Fang Wanliang.A block QR based power system state estimation algorithm[J].Proceedings of the CSEE,2003,23(8):51-55.

[19]郭瑞鵬,邵學儉,韓禎祥.基于分塊吉文斯旋轉的電力系統狀態估計[J].中國電機工程學報,2006,26(12):26-31. Guo Ruipeng,Shao Xuejian,Han Zhenxiang.A blocked Givens rotations based algorithm for power system state estimation[J].Proceedings of the CSEE,2006,26(12):26-31.

[20]Costa A S,Salgado R,Haas P.Globally convergent state estimation based on Givens rotations[J].IEEE Transactions on Power Systems,2014,29(5):2381-2390.

[21]Pajic S,Clements K A.Power system state estimation via globally convergent methods[J].IEEE Transactions on Power Systems,2005,20(4):1683-1689.

[22]Sousa A A,Torres G L,Caizares C A.Robust optimal power flow solution using trust region and interior-point methods[J].IEEE Transactions on Power Systems,2011,26(2):487-499.

Measurement Variance Estimation and Weights Configuration Algorithm Based on Orthogonal Transformation and Trust Region

Lu Zhigang1Tian Shasha1,2Shao Qi1,3Wu Jie1

(1.Hebei Key Lab of Power Electronics for Energy Conservation and Motor Drive Yanshan University Qinhuangdao 066004 China 2.State Grid Cangzhou Electric Power Supply Company Cangzhou 061000 China 3.State Grid Tianjin Procurement Company Tianjin 300304 China)

In the practical applications of state estimation,it is difficult to obtain the measurement variances and set the weights of measurements.With the computation of state estimation getting harder,the convergence of the existing measurement variance estimation algorithm can not be guaranteed.To solve this problem,an algorithm based on orthogonal transformation and trust region to estimate measurement variance and set weight is proposed in this paper.The numerical stability of measurement variance estimation is improved by the orthogonal transformation method,which can reduce the condition number of the coefficient matrix at each iteration.According to the reduced results of orthogonal transformation,combining with Lagrange multipliers method,a new measurement variance estimation model is established.The trust region is considered as the constraint of the objective function to ensure that each estimation result data is in a confidence interval.The simulation results of the IEEE 14 bus system,IEEE 30 bus system and IEEE 118 bus system show the effectiveness of the proposed algorithm.

State estimation,orthogonal transformation,trust region,measurement variance,weight,convergence

國家自然科學基金(61374098)、教育部高等學校博士學科點專項科研基金(20131333110017)和河北省研究生創新項目(00302-6370032)資助。

2015-06-15 改稿日期2016-04-27

TM734

盧志剛 男,1963年生,教授,博士生導師,研究方向為電力系統經濟運行分析與控制。

E-mail:Zhglu@ysu.edu.cn(通信作者)

田莎莎 女,1989年生,碩士研究生,研究方向為電力系統狀態估計。

E-mail:tianshasha0317@163.com

主站蜘蛛池模板: 91在线无码精品秘九色APP| 2022国产91精品久久久久久| 国产尤物jk自慰制服喷水| 又爽又大又光又色的午夜视频| 69av免费视频| 乱色熟女综合一区二区| 国产一区在线观看无码| 国产精品区视频中文字幕| 婷婷色中文网| 国产精品19p| 国产97视频在线观看| 国产欧美一区二区三区视频在线观看| 国产99在线观看| 天堂成人在线视频| 国产91在线|日本| 亚洲一区国色天香| 青青网在线国产| 国产中文在线亚洲精品官网| 成人av专区精品无码国产| 欧美黄网在线| 国产一区二区三区在线精品专区| 久久精品嫩草研究院| 亚洲男人的天堂久久精品| 亚洲无码视频图片| 国产高清在线丝袜精品一区| 亚洲天堂视频在线观看免费| 国内a级毛片| 久久青草精品一区二区三区| 国产第一色| 国产熟睡乱子伦视频网站| 日本成人在线不卡视频| 青青草国产免费国产| 中文字幕有乳无码| 国产微拍一区二区三区四区| 亚洲一区二区在线无码| 欧美第二区| 91国内视频在线观看| 福利视频一区| 欧美成a人片在线观看| 在线观看视频一区二区| 欧美日韩在线成人| a级毛片免费网站| 亚洲综合专区| 国产亚洲一区二区三区在线| 毛片国产精品完整版| 欧洲熟妇精品视频| 美女一级毛片无遮挡内谢| 免费人成视网站在线不卡| 国产91在线|日本| 日本三级欧美三级| 亚洲色偷偷偷鲁综合| 丁香亚洲综合五月天婷婷| 亚洲中文字幕久久无码精品A| 99热免费在线| 少妇被粗大的猛烈进出免费视频| 91蝌蚪视频在线观看| 免费毛片网站在线观看| 国产精品xxx| 久久精品国产999大香线焦| 又黄又湿又爽的视频| 99re这里只有国产中文精品国产精品| 亚洲国产成人综合精品2020| 91久久国产成人免费观看| 亚洲综合一区国产精品| 国产女人在线观看| 欧美一级黄色影院| 中国成人在线视频| 日日拍夜夜操| 中文字幕永久视频| AV不卡国产在线观看| 国产麻豆aⅴ精品无码| 毛片久久网站小视频| 久久精品中文无码资源站| 伊人国产无码高清视频| 四虎永久在线精品影院| 国产91av在线| 亚洲男人天堂2020| 亚洲性视频网站| 亚洲天堂免费观看| 国产精品欧美亚洲韩国日本不卡| 亚洲swag精品自拍一区| 国产乱子伦视频三区|