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

基于能量流的風電機組異常運行狀態(tài)預警軟件設計

2019-05-16 08:22:08李林鑫劉德順覃事剛
計算機應用與軟件 2019年5期

李林鑫 劉德順 覃事剛 張 帆

1(湖南科技大學機電工程學院 湖南 湘潭 411201)2(中南大學機電工程學院 湖南 長沙 410083)

0 引 言

近年來,通過對風電SCADA系統(tǒng)提供的大量數(shù)據(jù)進行分析,來對風電機組運行狀態(tài)識別、故障診斷與預警已經成為了一個研究熱點[1]。一方面,采用歷史SCADA數(shù)據(jù)建立研究對象的正常運行狀態(tài)模型,來對比此后某一段時間的數(shù)據(jù)變化,實現(xiàn)異常檢測。針對風電機組葉片[2]、齒輪箱[3-4]、發(fā)電機軸承[5]等關鍵部件或子系統(tǒng)開展運行狀態(tài)識別與異常預警。Kusiak等[6]采用霍特林T2統(tǒng)計方法,對風電機組的有功功率運行數(shù)據(jù)進行分析,以此來識別整機的運行狀態(tài)。另一方面,基于歷史SCADA數(shù)據(jù)建立研究對象的正常運行狀態(tài)預測模型,來對比實際運行狀態(tài)數(shù)據(jù)變化,實現(xiàn)異常檢測。Schlechtingen等[7]采用自適應神經模糊推理系統(tǒng)對風電SCADA數(shù)據(jù)進行處理,從而實現(xiàn)對風力發(fā)電機運行狀態(tài)的監(jiān)測。孟玲霞等[8]基于流行學習的預警方法、張少敏等[9]基于Storm實時流數(shù)據(jù)處理和Spark內存批處理技術,對風電機組齒輪箱故障預警進行研究。郭慧東等[10]提出了基于SCADA數(shù)據(jù)利用LSSVM與GMM相結合的風電機組變槳系統(tǒng)劣化狀態(tài)在線識別方案。梁穎等[11]提出了基于回歸預測模型和SCADA系統(tǒng)相配合的在線評估方案。

本文根據(jù)能量流概念對風電機組進行系統(tǒng)劃分[12],通過分析風電機組各個子系統(tǒng)輸入輸出參數(shù)之間的關系,建立風電機組運行狀態(tài)健康指標,進而利用跨平臺型良好的Qt開發(fā)軟件來實現(xiàn)對風電機組異常運行狀態(tài)預警的軟件設計。

1 基于能量流參數(shù)關系在線預警原理

1.1 能量流分析

風電機組作為一種復雜機電系統(tǒng),根據(jù)其內部能量傳遞關系,可以將其分為風能吸收環(huán)節(jié)、葉片機械能傳遞環(huán)節(jié)、發(fā)電機電能轉換環(huán)節(jié)以及變頻變流環(huán)節(jié)。根據(jù)這四個環(huán)節(jié)將風電機組劃分四個前后關聯(lián)的子系統(tǒng)。如圖1所示,氣流帶動風力發(fā)電機葉輪轉動,風能轉換成葉輪系統(tǒng)轉動的機械能,然后通過輪轂與發(fā)電機的連接機構將機械能傳遞給發(fā)電機轉子,切割磁力線,產生電能,最后再經由變頻器將電能傳輸?shù)诫娋W。圖中虛線區(qū)域表示風電機組內部能量流變化,每個子系統(tǒng)的輸入輸出參數(shù)與前后子系統(tǒng)的輸入輸出參數(shù)相關聯(lián),對每個子系統(tǒng)的相關信息進行分析,即可獲得風電機組當前的運行狀態(tài)。在本軟件設計中,葉輪系統(tǒng)、傳動系統(tǒng)、發(fā)電機系統(tǒng)以及變頻器系統(tǒng)為四個基本子系統(tǒng),在這四個基本子系統(tǒng)之下還可以根據(jù)實際采集到的SCADA數(shù)據(jù)進行細分,如:在變槳風力發(fā)電機中,葉輪系統(tǒng)中還包括變槳系統(tǒng)等輔助子系統(tǒng)。

圖1 風電機組能量流示意圖

根據(jù)風電機組能量流動方向確定每個子系統(tǒng)對應的SCADA輸入輸出參數(shù),表1展示了基本子系統(tǒng)與其對應的主要輸入輸出參數(shù)。對這些輸入輸出參數(shù)之間關系建模,以此來描述風電機組及其子系統(tǒng)運行狀態(tài)。在表1中還可以看到每個中間子系統(tǒng)的輸入輸出變量都具有承前啟后的特性,即上一個子系統(tǒng)的輸出參數(shù)作為下一個子系統(tǒng)的輸入?yún)?shù)。

表1 子系統(tǒng)名稱及其主要輸入輸出參數(shù)

本軟件按照能量流理論將風電機組劃分為四個子系統(tǒng)后,通過計算子系統(tǒng)的運行狀態(tài)健康指標,并分析其變化閾值,從而實現(xiàn)風電機組運行狀態(tài)的在線預警。

1.2 風電機組運行狀態(tài)預警模型

風電機組運行狀態(tài)預警大致分為狀態(tài)識別與預警兩個步驟。首先需要對風電機組運行狀態(tài)進行識別,提出風電機組運行狀態(tài)健康指標,然后根據(jù)健康指標進行統(tǒng)計分析,計算閾值,實現(xiàn)風電機組運行狀態(tài)預警。

對于實際風電機組及其各個子系統(tǒng)來說,其包含的輸入輸出參數(shù)之間的關系是復雜的,但都可以采用多項式逼近的方法來表示每個子系統(tǒng)中輸入輸出參數(shù)之間的關系。這樣,某輸入x與輸出y之間的關系可以表示為:

y=a0+a1x+a2x2+…+anxn

(1)

式中:n為多項式階數(shù),a0、a1、a2、…、an為常系數(shù)。

對于風電機組及其子系統(tǒng),基于SCADA參數(shù)歷史數(shù)據(jù)應用最小二乘法,可以計算出式(1)中的系數(shù)a0、a1、a2、…、an,根據(jù)SCADA參數(shù)當前數(shù)據(jù)可以計算出模型系數(shù)記作b0、b1、b2、…、bn,于是,可以得到風電機組及其子系統(tǒng)運行狀態(tài)健康指標[13]:

(2)

式中:xmax、xmin為當前數(shù)據(jù)集中,輸入變量x的最大值與最小值。這里,健康指標C表示風電機組某個時刻的數(shù)據(jù)集計算出來的運行模型相對于正常運行時刻的數(shù)據(jù)集計算出來的標準運行模型的差異。健康指標C的值越大,說明該時刻運行狀態(tài)相對于正常運行狀態(tài)差異越大,發(fā)生異常的可能性越大;C的值越小,說明該時刻運行狀態(tài)相對于正常運行狀態(tài)差異越小,發(fā)生異常的可能性越小,即運行正常的可能性越大。

由于風速具有很強的隨機性,風電機組運行狀態(tài)也表現(xiàn)出隨機性,其健康指標C實際上是一個隨機變量。各個時刻的健康指標組成新的數(shù)據(jù)集Dc={C1,C2,…,Cn},并具有一定的隨機分布規(guī)律。采用核密度估計法對其進行處理,構造其密度函數(shù),一維核密度估計函數(shù)的一般式為:

(3)

式中:K(x)為核函數(shù);n為數(shù)據(jù)集Dc中數(shù)據(jù)量;h為核密度函數(shù)的帶寬。帶寬h的取值影響著核密度估計的精確度,如果h太小或太大都會對計算結果產生劇烈影響。

選擇高斯核函數(shù),則可以對數(shù)據(jù)集Dc統(tǒng)計分析得到核密度估計函數(shù)為:

(4)

式中:最優(yōu)帶寬hopt的計算公式為:

(5)

(6)

2 在線預警軟件的設計與實現(xiàn)

2.1 應用場景

現(xiàn)代大型風電機組大多裝配了相應的SCADA系統(tǒng),對風電機組進行狀態(tài)監(jiān)控。信息收集與儲存是其核心業(yè)務,而在風電機組運行狀態(tài)報警方面,對風電機組各子系統(tǒng)的物理聯(lián)系缺乏有效分析與解釋,導致SCADA系統(tǒng)只能在故障發(fā)生時給出報警信息,而無法提前預警以避免事故發(fā)生。利用SCADA系統(tǒng)提供的大量歷史數(shù)據(jù),對風電機組運行狀態(tài)進行分析是本軟件的主要功能。本軟件應用架構如圖2所示,通過訪問SCADA系統(tǒng)服務器的數(shù)據(jù)庫獲取實時數(shù)據(jù)進行分析。

圖2 預警軟件在SCADA系統(tǒng)中的應用構架

在核心業(yè)務模塊中,給出了一個風電機組實時閾值計算方法,讓原本固定不變的閾值隨SCADA數(shù)據(jù)的變化而變化。除了訪問SCADA系統(tǒng)服務器數(shù)據(jù)庫以外,本軟件完全獨立于SCADA系統(tǒng)本身進行獨立的風電機組運行狀態(tài)預警。即對于不同的SCADA系統(tǒng),只要正確設置其數(shù)據(jù)存儲格式,都可以通過本軟件進行訪問分析。

成功連接數(shù)據(jù)庫后,軟件會自動讀取SCADA數(shù)據(jù),并將數(shù)據(jù)傳遞給軟件的其余部件,通過人機交互接口顯示當前數(shù)據(jù)變化。作為一個基于SCADA數(shù)據(jù)的在線預警軟件,本軟件不需要在運行電腦上額外配置環(huán)境變量,不需要增加任何物理設備,對于風電機組維護成本的影響幾乎為零。

2.2 軟件結構

軟件結構如圖3所示,可以看出,本軟件分為數(shù)據(jù)準備層、數(shù)據(jù)操作層、繪圖操作層與輔助功能層等四個層次。

圖3 軟件結構圖

數(shù)據(jù)準備層:本層為數(shù)據(jù)訪問模塊,具備軟件外部數(shù)據(jù)訪問以及軟件內部數(shù)據(jù)分流功能。首先,訪問SCADA服務器中的數(shù)據(jù)庫作為本軟件的數(shù)據(jù)來源,以驅動軟件內置各子系統(tǒng)的數(shù)學模型;然后將讀取到的SCADA數(shù)據(jù)按照軟件不同功能模塊傳輸給繪圖操作層、數(shù)據(jù)操作層。

數(shù)據(jù)操作層:本層包含兩個數(shù)據(jù)處理模塊:一個是基于能量流參數(shù)關系的風電機組運行狀態(tài)健康指標計算模塊;另一個是基于統(tǒng)計學的非參數(shù)核密度估計的閾值計算模塊。通過接收來自數(shù)據(jù)準備層的數(shù)據(jù)后進行計算分析。首先計算實時風電機組運行狀態(tài)指標值,傳輸給繪圖操作層以及運行狀態(tài)指標實時閾值分析模塊,在運行狀態(tài)指標實時閾值分析模塊中,進行實時閾值分析,新讀取的數(shù)據(jù)與實時更新的歷史數(shù)據(jù)集閾值進行對比,其對比結果同樣傳遞給繪圖操作層。

繪圖操作層:從數(shù)據(jù)準備層與數(shù)據(jù)操作層獲取數(shù)據(jù)信號,根據(jù)實時的數(shù)據(jù)信號繪制實時圖像,包含基本數(shù)據(jù)圖表與狀態(tài)數(shù)據(jù)圖表兩個繪圖模塊。基本數(shù)據(jù)圖表模塊由主界面基本參數(shù)繪圖與鼠標點擊參數(shù)繪圖模塊組成;狀態(tài)數(shù)據(jù)圖表框架由實時運行狀態(tài)健康指標繪制模塊、實時健康指標閾值繪制模塊與預警狀態(tài)繪圖模塊三部分組成。其中主界面基本繪圖模塊為軟件主界面基本參數(shù)實時顯示,而鼠標點擊參數(shù)繪圖則是將本軟件所分析的各個參數(shù)單獨繪制成圖,便于運維人員在閱讀數(shù)據(jù)的同時,更直觀掌握當前風電機組參數(shù)變化情況。狀態(tài)數(shù)據(jù)圖表框架中包含的三個模塊是本軟件設計的核心內容,它們反映了當前風電機組實時運行狀態(tài)、實時閾值變化以及風電機組是否處于預警狀態(tài)。當風電機組運行狀態(tài)指標突破健康指標閾值時,預警狀態(tài)繪圖模塊將根據(jù)繪圖策略更改當前主界面顯示狀態(tài),提醒工作人員注意風電機組運行狀態(tài)。

輔助功能層:作為人機交互層,本層的設計以便于工作人員隨時查看風電機組運行狀態(tài)預警記錄為目的,設計了預警日志輸出模塊、參數(shù)設置模塊以及運行狀態(tài)提醒模塊。當本軟件監(jiān)測到風電機組可能存在異常運行狀態(tài)時,預警日志模塊就會在軟件主界面打印一條預警事件的記錄,工作人員隨時可以通過主界面了解風電機組的預警情況。該模塊還提供了日志文檔功能,可以將預警日志輸出到文本文檔,保存至工作人員指定目錄下。對于不同風電機組,其初始模型參數(shù)必然有所差異,因此,利用參數(shù)設置模塊對不同風電機組的初始運行狀態(tài)模型參數(shù)進行設置。運行狀態(tài)提醒模塊是本軟件層的最后一個模塊,當風電機組運行狀態(tài)出現(xiàn)異常預警時,從無異常狀態(tài)切換至異常,與異常預警消失后,從異常提醒狀態(tài)恢復到無異常狀態(tài)。

2.3 基于Qt的軟件模型算法實現(xiàn)

Qt作為一個集成開發(fā)環(huán)境,為開發(fā)人員提供了豐富的類庫以及其特有的信號與槽機制[14],易實現(xiàn)風電機組運行狀態(tài)數(shù)據(jù)的在線分析,并輸出預警結果。本軟件在數(shù)據(jù)分析過程中包含了一個時間窗口模型與兩個數(shù)據(jù)處理方法。

2.3.1時間窗口模型

為了實現(xiàn)數(shù)據(jù)的實時更新,引入一個滑動時間窗口模型來實現(xiàn)對獲取到的數(shù)據(jù)建立一個動態(tài)數(shù)據(jù)集D,如圖4所示。

圖4 滑動時間窗口模型

圖4中,n表示數(shù)據(jù)集D的時間寬度,不同的SCADA系統(tǒng)提供的數(shù)據(jù)采樣頻次不同。本文采用的SCADA系統(tǒng)數(shù)據(jù)采樣頻率為每秒一次,因此數(shù)據(jù)集寬度n在數(shù)值上與時間計數(shù)相同,即數(shù)據(jù)集D包含n個元素,就表示數(shù)據(jù)集D包含n秒鐘的數(shù)據(jù),D(ti-1)表示ti-1時刻的數(shù)據(jù)集,Δd表示n個時間單位后新增加的數(shù)據(jù)子集。此滑動時間窗口模型的具體算法設計詳見算法1。

算法1Sliding Time Window Model(d)

輸入:當前時刻待分析SCADA數(shù)據(jù);

輸出:動態(tài)數(shù)據(jù)集D;

{//start algorithm

步驟1初始化數(shù)據(jù)集基本參數(shù)

int n;

//設置窗寬;

new QVectorD(n);

//初始化

QVector容器D用于存儲數(shù)據(jù);

D<

//將數(shù)據(jù)d傳遞給D;

步驟2控制數(shù)據(jù)集大小

while(D.size()>=n){

//當D的大小達到窗寬n時,開始循環(huán);

D.removefirst();

//移除D中第一個元素。

…;

}

//end while;

}

//end algorithm;

在算法1所示的時間窗口模型實現(xiàn)過程中,利用while循環(huán)可以保證動態(tài)數(shù)據(jù)集的大小固定,當數(shù)據(jù)集大小達到時間窗口設定值時,刪除第一個數(shù)據(jù),同時新進的數(shù)據(jù)存入數(shù)據(jù)集的末尾。基于時間窗口模型,可以實現(xiàn)風電機組運行狀態(tài)健康指標與預警閾值實時動態(tài)分析。

2.3.2風電機組運行狀態(tài)健康指標算法

第1節(jié)已經討論了關于風電機組運行狀態(tài)健康指標的計算方法,在計算健康指標的過程中,包含了風電機組實時運行狀態(tài)建模,以及實時健康指標計算兩部分。首先利用最小二乘多項式擬合來求解風電機組實時運行狀態(tài)模型的參數(shù),然后根據(jù)式(2)計算健康指標,具體算法設計詳見算法2。

算法2Healthy Criterion algorithm(x,y)

輸入:目標子系統(tǒng)原始輸入輸出參數(shù)SCADA數(shù)據(jù);

輸出:目標子系統(tǒng)運行狀態(tài)健康指標C;

{

//start algorithm;

步驟1初始化數(shù)據(jù)集基本參數(shù)

new QVectorX(length);

new QVectorY(length);

CriterionCalculation(int n, B[n], x, y)

{

//n為擬合階次;length為時間窗寬;

//B[n]為標準多項式擬合參數(shù)b0, b1, b2, …, bn;

X<

Y<

//輸入數(shù)據(jù)傳遞給X、Y;

while(X.size()>=length;){

X.removefirst();

Y.removefirst();

步驟2利用最小二乘法計算式中參數(shù);

LeastSquares(X,Y,n){

…;

return A[n];

};

//A=(a0, a1, …, an)

步驟3根據(jù)式計算風電機組運行狀態(tài)健康指標C;

calculationC(A[n],B[n]){

…;

return C;

}

}

//end while;

}

//end CriterionCalculation();

}

//end algorithm;

在算法2的自定義函數(shù)calculationC(A[n],B[n])中,根據(jù)式(2)計算風電機組運行狀態(tài)健康指標時,采用Romberg積分算法[15]來解決其中的積分問題。Romberg算法是建立在Simpson公式、Newton-Cotes公式基礎上的一種加速積分算法。利用Richardson外推加速收斂技術,使該算法可以在不增加計算量的前提下提高誤差精度。

2.3.3風電機組運行狀態(tài)預警閾值算法

利用統(tǒng)計學中非參數(shù)估計原理,可以對算法2計算出來的風電機組健康指標C進行閾值計算。在這一部分中,要實現(xiàn)風電機組運行狀態(tài)預警閾值的實時分析。首先,收集實時健康指標C,建立一個動態(tài)數(shù)據(jù)集;其次,根據(jù)式(4)與式(5)計算健康指標C的密度函數(shù);最后,根據(jù)式(6)計算閾值。這部分算法詳見算法3。基于算法3,軟件可以完成實時預警閾值的計算,返回預警閾值所需要的參數(shù),再根據(jù)式(6)便可以計算出隨SCADA數(shù)據(jù)變化的風電機組運行狀態(tài)的預警閾值。

算法3Warning Threshold algorithm(c)

輸入:待分析風電機組運行狀態(tài)健康指標C;

輸出:實時風電機組預警閾值Cth;

{

//start algorithm;

步驟1初始化數(shù)據(jù)基本參數(shù)

QVectorC(length);

ThresholdCalculation(double c, int length){

//c為實時健康指標值;length為時間窗寬;

步驟2輸入數(shù)據(jù)傳遞給C;

C<

while(C.size()>=length; ){

C.removefirst();

步驟3根據(jù)式中最優(yōu)帶寬hopt;

OptimalBandwidth(C){

…;

return h;

};

步驟4計算數(shù)據(jù)集的數(shù)字特征

NumericalCharacteris(){

NuclearDensity(C,h){

…;

//根據(jù)式進行核密度估計分析;

};

…;

return Cmean、S、Sigma;

}

//Cmean、S分別運行狀態(tài)健康指標

//數(shù)據(jù)集的期望與標準差;

//sigma為核密度估計標準差;

…;

return Cth;

//根據(jù)式計算預警閾值Cth;

}

//end while;

}

//end ThresholdCalculation();

}

//end algorithm;

基于時間窗口模型、風電機組運行狀態(tài)健康指標算法與預警閾值算法,可以實現(xiàn)一個完整數(shù)據(jù)處理周期。本軟件設計的風電機組運行狀態(tài)預警流程如圖5所示。

圖5 預警流程圖

在圖5中,SCADA數(shù)據(jù)從開始出發(fā),首先計算出某個風電機組及其子系統(tǒng)輸入輸出參數(shù)關系的運行狀態(tài)健康指標C,再根據(jù)健康指標C的值進行實時統(tǒng)計分析并計算預警閾值Cth,然后做出判斷。若滿足判定條件,即無異常預警,直接進入下一輪計算,若不滿足判定條件,則給出預警提醒,再繼續(xù)下一個周期的分析。

3 實例分析

本軟件采用Qt creator開發(fā)軟件,版本號為5.10.1,操作系統(tǒng)環(huán)境為Windows 10。以某風電場2 MW直驅式風電機組的SCADA系統(tǒng)數(shù)據(jù)為監(jiān)測對象,提取風電機組SCADA系統(tǒng)4天運行數(shù)據(jù)進行在線監(jiān)測模擬。首先將SCADA數(shù)據(jù)導入服務器數(shù)據(jù)庫,運行本軟件并配置好各子系統(tǒng)標準模型參數(shù),如圖6所示,點擊確定生效。然后連接數(shù)據(jù)庫,模擬在線監(jiān)測過程。其中A0、A1分別對應正常運行狀態(tài)下的多項式回歸模型參數(shù)。

圖6 葉輪子系統(tǒng)參數(shù)設置

圖7展示了當風電機組正常運行時,在軟件主界面中可以看到實時功率曲線穩(wěn)定在額定功率附近,總體運行狀態(tài)進度條也與其余子系統(tǒng)進度條一樣,處于一個常規(guī)的狀態(tài)。而當軟件檢測出風電機組即將偏離正常運行狀態(tài)時風電機組運行狀態(tài)進度條數(shù)值將不斷增長,并且在軟件主界面顯示預警消息,如圖8所示。

圖7 運行狀態(tài)主界面(非預警狀態(tài))

圖8 運行狀態(tài)主界面(預警狀態(tài))

從圖8中可以看出,葉輪系統(tǒng)出現(xiàn)異常運行狀態(tài)預警,其預警進度條預警狀態(tài)值已經達到92%。而風電機組總體運行狀態(tài)進度條仍然處于“非常好”的狀態(tài),說明風電機組總體的運行狀態(tài)與子系統(tǒng)運行狀態(tài)不一定完全相同,即風電機組子系統(tǒng)出現(xiàn)異常時風電機組整體運行狀態(tài)并不一定會同時出現(xiàn)異常。由于一個輸入輸出參數(shù)關系代表著一個子系統(tǒng)或者零部件運行狀態(tài)模型,所以基于參數(shù)關系的預警也就提供了發(fā)生異常的零部件方面的信息,這有利于針對性的維護。

在軟件運行過程中,工作人員可以通過點擊主界面中各個系統(tǒng)對應的狀態(tài)進度條查看當風電機組運行狀態(tài)健康指標與實時預警閾值的變化。圖9顯示了風電機組風速與功率之間健康指標以及其預警閾值的變化,橫坐標表示當前時刻到過去360 s的時間跨度,縱坐標表示風電機組運行狀態(tài)健康指標與預警閾值的變化范圍。

圖9 總體運行狀態(tài)評估指標與閾值

4 結 語

風力發(fā)電機運行狀態(tài)預警的實時性對于風電機組預警十分關鍵,本文所設計實時在線預警軟件重點在于構建一個實時預警的基本框架,讓使用者更加方便快捷地閱讀風電機組的實時運行狀態(tài)。本文注重風電機組運行狀態(tài)的時效性與整體性,將原本相對獨立的各項SCADA參數(shù)按照能量流理論進行系統(tǒng)分類,并整理出對應的輸入輸出參數(shù),有利于用戶通過預警信息迅速查找異常零部件。將軟件的數(shù)據(jù)操作設計成獨立的數(shù)據(jù)操作模塊,有利于引入不同的狀態(tài)預警數(shù)學模型,且有較強的擴展性。

主站蜘蛛池模板: 成人福利在线看| 伊人久久福利中文字幕| 四虎国产成人免费观看| 男女男精品视频| 亚洲人成网站日本片| 人妻少妇久久久久久97人妻| 中文字幕在线看视频一区二区三区| 亚洲福利视频一区二区| 国产xx在线观看| 日韩欧美国产中文| 亚洲国内精品自在自线官| 国产99视频在线| 91小视频版在线观看www| 亚洲国产成人在线| 在线免费观看AV| 四虎精品黑人视频| 黄色网站不卡无码| 91免费国产高清观看| 一区二区日韩国产精久久| 91精品国产自产91精品资源| 伊人蕉久影院| 国产精品综合久久久 | 波多野结衣无码AV在线| 欧洲欧美人成免费全部视频| 精品成人免费自拍视频| 国产精品无码作爱| 国产午夜无码专区喷水| 国产不卡网| 香蕉久久永久视频| 99视频在线观看免费| 亚洲h视频在线| 国产精品视频猛进猛出| 四虎亚洲国产成人久久精品| 中文字幕亚洲无线码一区女同| 亚洲成人免费看| 8090成人午夜精品| 欧美色图第一页| 国产精品视频观看裸模| 日韩午夜福利在线观看| 老司国产精品视频91| 18禁黄无遮挡免费动漫网站| 综合色亚洲| 日韩欧美中文字幕在线精品| 成人在线视频一区| 2020国产免费久久精品99| 免费国产在线精品一区| AV熟女乱| 欧美精品xx| 亚洲国产日韩在线成人蜜芽| 久久免费视频6| 欧美精品高清| 青青草原国产av福利网站| 不卡网亚洲无码| 毛片久久久| 精品视频在线一区| 99中文字幕亚洲一区二区| 欧美日韩国产在线人| 国产精品自拍露脸视频| 91精品国产自产91精品资源| 国产噜噜噜| 少妇极品熟妇人妻专区视频| 青青青国产视频手机| 亚洲天堂网2014| 午夜精品久久久久久久2023| 欧美亚洲国产精品第一页| 国产一级毛片yw| 久久一级电影| 69综合网| 国产69囗曝护士吞精在线视频| a免费毛片在线播放| 国产亚洲视频免费播放| 国产在线精品人成导航| a毛片免费看| 热99精品视频| 日韩无码视频播放| 欧美国产在线看| 五月丁香伊人啪啪手机免费观看| 国产高清在线精品一区二区三区| 综合色区亚洲熟妇在线| 国产第四页| 亚洲成人一区二区三区| 欧美α片免费观看|