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

基于改進雷達圖綜合指標的渠道系統PI控制參數優化

2022-11-07 07:35:58李梟華管光華
節水灌溉 2022年10期
關鍵詞:優化

李梟華,管光華

(武漢大學水資源與水電工程科學國家重點實驗室,武漢 430072)

0 引言

隨著水資源稀缺與供需矛盾的日益明顯,節約水資源、提高水資源的利用率已成為實現水資源可持續利用的關鍵。中國是農業大國,2020 年全國農業用水占用水總量的62.1%,灌溉水有效利用系數為0.565[1],湖北省某灌區渠系水利用系數為0.713[2],均有較大提升空間。因此,灌區自動化建設亟待發展,以期提高灌區的用水效率、降低運行成本。

明渠系統作為灌溉輸配水的主要手段,其自動化是灌區管理信息化和自動化的重要標志。渠系自動控制目標是為灌區用戶提供更加足量及時、安全可靠的供水服務[3],其核心是控制算法,即從輸入信息(水位、流量)到輸出的控制動作(閘門開度),以維持相應的控制目的(如下游常水位)的計算邏輯和方法。國內外渠道控制算法可主要分為3類:以傳遞函數為基礎的經典控制算法,以PID類反饋控制算法為代表成果,如P、PI、PID、PIF 等算法;以狀態空間為基礎的現代控制算法,以LQR、MPC 為代表;以人的思維方式為基礎研究復雜系統控制的智能控制算法,如魯棒控制、自適應算法和神經網絡等[4]。由于現代控制算法和智能控制算法尚處于試驗和研究階段,在實際工程中,原理簡單、魯棒性強的PID類經典控制算法應用較為常見,尤其是分布式比例積分(PI)類反饋控制算法研究最為廣泛[5-9]。近年來的研究趨勢是用智能方法結合經典的PID類經典算法,以獲得更良好、更全面的控制性能。韓延成等[10]推導了RBF人工神經網絡整定PID輸水控制算法,張雨萌等[11]采用粒子群算法進行高效PI 控制參數尋優,葉雯雯[12]采用自適應PI 控制,針對多種優化目標進行分析、比選,比常規PI算法的閘門超調量降低了9%~21%。

渠道系統PI 控制的首要難點是確定控制參數,需要給定一個評價控制性能優劣的量化指標,按照該指標最大化或最小化進行優化求解。目前在渠道控制性能指標領域相關研究較少。美國土木工程師協會(ASCE)渠道控制算法工作組于1998 年提出了一系列基本控制性能指標[13],對渠系控制性能進行優劣衡量。但該系列指標均只針對水位、流量或閘門的其中一項。而實際工程中僅僅追求單項性能最優,往往會導致其他方面性能的大幅惡化。管光華等[14]對現有性能指標進行去量綱處理,提出了綜合指標GI,并以GI 為優化準則整定了PI 控制參數,并與單一指標作為優化準則的結果進行了對比,結果顯示追求單一指標的最優難以取得滿意結果。

為兼顧控制系統的多方面性能,實現安全可靠、性能均衡的控制目的,本文從輸配水渠道系統綜合性能評價的問題出發,選取了水位、流量、閘門開度和穩定過渡時間等多個優化目標,提出了基于改進雷達圖的綜合指標實現多目標綜合評價,并以漳河灌區第四干渠為例,驗證了該綜合指標作為PI控制參數優化目標的有效性。

1 材料與方法

1.1 PI控制算法

PID 算法,即比例積分微分算法,是工業控制中常見的反饋控制算法。在渠道控制算法中,PID 控制算法因簡單有效而應用最廣。PID 控制算法包括3 個反饋環節,分別是比例環節,積分環節和微分環節,分別對應比例參數Kp,積分參數Ki和微分參數Kd。其中,微分環節對干擾較為敏感,容易造成超調過大、系統失穩的問題,因而在許多工業控制算法中,常使用沒有引入微分環節的PI控制器[9]。在本文中的反饋算法邏輯為:將傳感器檢測得到的水位y(t)反饋到控制中心,該水位和預先設定的目標水位ytarget做比較,其差值為水位誤差e(t),其作為PI 算法的輸入量,進行比例、積分兩個環節的線性累加計算,得到輸出量為流量u(t),原理見下式:

式中:t為運行時間或仿真時間,min。

將上式進行離散處理可得位置式PI控制的離散表達式:

考慮到渠道系統的水動力過程非線性強,復雜性高的特點,對式(3)在時間步上進行作差,得到增量式PI 控制表達式:

即:

式中:DT為實際采樣步長或仿真步長,min。

增量式PI 控制的優點為運算工作量小,手動∕自動切換沖擊小,誤動作影響小[15]。這些優點更加契合輸配水渠道系統的控制要求。

1.2 渠道系統控制結構

由于渠道水波傳遞的滯后性,不能只依靠反饋控制來進行水位的調節,否則容易出現水位變化超限甚至失穩。在取水流量已知的情況下,通常利用取水流量計劃的信息,在取水動作發生的同時,更新目標流量,計算實時的流量校正量,從而使過閘流量向目標流量靠攏,稱此步驟為前饋校正,其前饋校正流量Δu'計算式為:

綜合前饋和反饋后的控制過閘流量為:

計算得到控制過閘流量后,需要根據當地閘門過流特性反算出符合該過閘流量的閘門開度,然后讓執行機構(人工∕自動)執行閘門動作。

式中:G為閘門開度,m;F為當地過閘流量-閘門開度計算公式。

整體控制邏輯見圖1:位于渠池下游控制點的水位傳感器將監測的下游水位實時反饋給PI 反饋控制器,計算得到相應的反饋流量增量,加上前饋校正流量,兩者之和即為該時刻的流量變化量,可得下一時刻的流量,按此流量反算出相應的閘門開度,指導當地人工或者電機執行,以達到控制下游水位在目標值附近的控制目的。其中,PI 控制器計算邏輯見圖2。

圖1 渠道控制結構示意圖Fig.1 Schematic diagram of channel control structure

圖2 PI控制器計算示意圖Fig.2 PI controller calculation diagram

1.3 多目標選取

下游常水位運行模式中,控制對象是每個渠池的上游閘門,控制目標是使該渠池下游水位保持在目標值附近(目標值通常為設計水位,即恒定常數),見圖3。

圖3中,通過控制閘門維持下游水位在目標值附近,保持為常數。Qmax為渠道加大流量,對應的折線為加大流量下的水面線;Q0為最小流量,一般是生態流量,其對應著最小流量下的水面線。

圖3 下游常水位控制模式示意圖Fig.3 Schematic diagram of downstream normal water level control mode

灌溉渠道系統中的分水建筑物大多為分水閘,由于經濟性和實用性,分水閘閘門要保持過流穩定,就要求閘前水位盡量保持在恒定值,也有利于減少分水閘閘門動作調整量。因閘前壅水效應的存在,渠池下游段通常能保持較為水平的水位,所以大多數分水閘都分布在下游段,這就要求閘前水位要盡量保持在一個恒定值,才能保證供水的精確性和穩定性。

在衡量供水渠道系統的控制代價方面,渠道系統通過控制閘門開度來調整過閘流量,進而控制渠池下游水位,而閘門動作需要人力啟閉,半自動電力輔助啟閉,或者全自動電動閘門啟閉,需要消耗的人力和電力就是主要的運行代價。另外,閘門動作量越大,動作次數越頻繁,相關的部件磨損就越快。與此類似,在下游常水位的主要目標之外,從安全性和穩定性來考慮,通常希望流量變化量不要過大。

因為控制參數變化會導致系統的控制結果在不同方面的性能上發生變化,并且這種變化通常是互相矛盾的,這就需要提出一個兼顧各方面控制性能的綜合性能指標,來衡量、評價其控制結果,進而優化控制參數。針對供水渠道系統,從水位、流量、閘門開度和穩定時間4 個方面,選取以下5 個性能指標作為優化目標,用于供水渠道系統PI 控制參數優化[13]。

(1)最大絕對誤差MAE(Maximum absolute error)描述了水位波動極值與目標值的差距,用于衡量水位控制的安全性:

式中:yt為t時刻下游測量水位,m;ytarget為控制點目標水位,m;Tmax為總仿真時長或總運行時長,min;DT為測量水位采樣步長或仿真設置步長,min。

(2)絕對值誤差積分IAE(Integral of absolute error)描述了系統過渡過程中水位波動的累積量,用于衡量水位控制的平穩性:

(3)絕對流量變化積分IAQ(Integrated absolute discharge change)描述了系統過渡過程中流量變化的累積量,用于衡量流量控制的平穩性:

式中:t1為取水口開始取水時刻,即系統開始發生變化的時刻,min;t2為系統進入穩定狀態的時刻,穩定狀態一般指水位變化和閘門動作都在t2時刻后一直小于某個設定閾值的狀態,min;Q(t)為t時刻渠池上游閘門過閘流量,m3∕s。

(4) 絕對閘門開度積分IAW(Integrated absolute gate movement)描述了系統過渡過程中閘門開度變化的累積量,用于衡量閘門控制的平穩性:

式中:W(t)為t時刻渠池上游閘門開度,m。

(5)穩定狀態過渡時間Tstable(Time of transition to stable state)指系統從受外部擾動導致狀態開始變化時,受控制算法進行控制應對擾動,直到回復到穩定狀態時所用的時間:

可以看出,上述5 個優化目標都采用最小化的優化選項,即優化目標數值越小,則表示系統對應的性能越優。

1.4 基于改進雷達圖的綜合指標計算

在確定供水渠道系統控制性能評價的多個目標后,需提出一個綜合指標的計算式來對上述5 個性能指標進行綜合計算。本文借鑒用雷達圖來描述系統控制綜合性能的思路[16],改進了雷達圖繪制中的數據處理方法,提出采用各項性能指標與初篩樣本平均值之比作為雷達圖中心到各頂點的距離,并以該改進后的雷達圖面積作為綜合指標。

上述各優化目標的指標量綱不同且數值數量級差別較大,直接使用會對優化結果的有效性和均衡性造成影響[13]。例如,某一指標的值遠大于其他指標,則優化結果就會明顯傾向于使得該指標最小的優化方向,而忽視了其他指標,導致優化后控制參數過度追求該方面的控制性能,而在其他方面的性能明顯低于期望。為消除上述差異對優化結果的影響,需要對優化目標進行標準化處理。

在數據處理之前,需要先得到符合基本管理要求的初篩樣本。初篩樣本是指符合初步篩選標準的控制參數樣本及與之相應的過程數據及其計算得到的性能指標。初步篩選的標準一般為,水位不超過設計超高,流量不超過設計流量且不低于生態流量,閘門開度變化在允許值上下限內等。具體數值視渠道規模、工況類型和管理要求而定,后文將針對具體算例進行說明。

數據標準化最常用的方法是min-max 標準化法和Z-score標準化法[17,18]。但前者要求確定各項性能指標的最小值和最大值,而控制對象為渠道系統這一復雜非線性系統時,難以給出其在水位、流量等方面的指標極值,故不適用;后者利用樣本均值和方差進行標準化處理后的數據弱化了原數據的直觀意義,且計算結果有正有負,而負數無法在雷達圖上表示,也不適用于本文對象。因此,在標準化處理時,本文提出使用初篩樣本平均數作為分母進行標準化處理,即:

式中:X為性能指標原值;Xˉ為初篩樣本該項性能指標的平均值;X'為處理后的性能指標計算值,作為雷達圖的繪制依據。

將標準化處理后的各項優化指標用于繪制雷達圖,見圖4。定義雷達圖中心到各頂點的距離為l,其長度為該指標與初篩樣本平均值之比(l>0),中心與各頂點連線的夾角相等,均為θ=72°。

圖4 綜合性能雷達圖Fig.4 Comprehensive performance radar chart

本文創新地提出用初篩樣本平均值來進行標準化處理,其主要優勢有:其一,平均值更能代表初篩樣本參數組合的平均表現,而管理者通常希望得到的優化結果能都在各方面都優于平均表現。用指標原值與平均值之比繪制雷達圖,當中心點到該方向頂點距離為l=1 時,代表著該項指標等于樣本平均值,而l<1 時,表示該項指標優于平均表現,符合管理者對于“各項指標均較優”的期望。若采用原值與最大值之比,則所有情況均滿足0<l<1,不能反映出管理者的期望;其二,采用樣本平均值比用最大值更可靠,因為最大值遠離樣本,其分布存在較大的不確定性,而平均值受樣本單值分布的影響小,所以將平均值用于性能評價的參考具有更好的可靠性;另外,如果采用最大值,部分指標的最大值均遠離樣本中心,會導致繪制所得的雷達圖各段l長度過小(<0.5),不利于結果的呈現。綜上,用初篩樣本平均值來進行標準化處理在代表性、可靠性和結果呈現3個方面都具有優勢。

基于改進后的雷達圖,本文提出用該圖形面積ARD(Area of radar diagram)作為渠道控制性能綜合評價的指標,其計算式如下:

雷達圖面積受l1,l2,l3,l4,l5的共同影響,能夠全面、綜合地反映控制性能的優劣。從運算法則來看,直接用各項指標之和作為綜合指標,容易出現“大數吃小數”的現象,弱化了數值較小的指標的影響力;用各項指標之積作為綜合指標,容易過度強化了變異幅度大的指標的影響力;用雷達圖面積ARD作為綜合指標,其計算式包含加法和乘法運算,能夠改進簡單直接的加和、乘積作為綜合指標的缺點,如避免純加法運算導致“大數吃小數”[19]。因此,所提出的綜合指標ARD能夠兼顧不同數值大小和變異范圍的各指標的影響,更加合理地評價綜合性能。

2 結果與分析

2.1 測試算例——漳河灌區四干渠

測試算例所用的渠段取自漳河灌區四干渠,起始四干進水閘,終至田家沖節制閘,總長13.62 km,由5 個節制閘劃分為4 個渠池。其渠道系統示意圖見圖5,仿真渠段幾何參數及水力參數見表1。

表1 漳河灌區四干渠參數表Tab.1 Parameter table of the fourth main channel in Zhanghe Irrigation Area

圖5 漳河灌區四干渠示意圖Fig.5 Schematic diagram of the fourth trunk canal in Zhanghe irrigation area

根據四干渠進口閘2019 年流量監測記錄,在5-8 月供水時其特征流量為6 m3∕s。又根據位于渠池3 下游的高店分水口設計流量為2 m3∕s,設置典型供水工況為T=2 h 時,高店分水口開始階躍分水。總仿真時長為6 h,仿真步長為DT=3 min。

仿真平臺在渠道內采用一維圣維南方程組來描述水力特性。因缺少各閘門過流實測率定數據,渠池間的節制閘過閘流量計算式用美國中亞利桑那調水工程CAP 公式描述其過流特性。CAP公式如下:

參數優化方法采用基本的網格法,以便看出各指標和綜合指標的分布規律。網格尋優法主要是通過選擇參數空間中的代表點來代表一定區間的參數空間,從而避免了對參數空間的全局搜索,進而減少了計算成本,以在一定程度上尋優得到合適的參數。本算例中初次優化網格大小為1,二次優化網格大小為0.1。

根據四干渠的規模大小和設計工況的流量變幅,本算例中對控制參數初步篩選的標準為:下游水位不超過目標水位3 cm,末尾渠道(渠池4)流量不少于3.5 m3∕s。本算例僅對水位超高和下游最小流量要求作為篩選標準,對于其他工程,還可以考慮閘門運動限速,最低水位要求等等。另外,因四干渠流量設計值較大(33 m3∕s),其運行流量遠小于該設計值,故不需限制。

因評價對象為多渠池系統,而每個渠池都有對應的性能指標,故在計算系統的整體性能指標時,最大水位偏差和穩定過渡時間取4個渠池的最大值,水位偏差累積量、流量變化累積量和閘門動作累積量取4個渠池的平均數。

2.2 基于綜合指標的多目標優化結果

初次優化網格范圍設置為:Kp=1~20,Ki=1~20,共400 個參數組合。采用Matlab 進行水動力過程仿真計算,以及對控制過程數據處理后進行三維曲面圖繪制,得到各單項性能指標分布情況見圖6(a)~圖6(e)。經初步篩選后,將滿足標準的初篩樣本各單項性能指標的平均值用于標準化處理,計算得到綜合指標ARD分布情況見圖6(f)。

圖6 初次優化各指標分布圖Fig.6 Indicator distribution of initial optimization

控制參數影響規律分析:最大水位偏差MAE隨著積分參數Ki的增大,有先減后增的趨勢,與比例參數Kp呈正相關,且其相關性受Ki的影響,Ki越大,MAE隨Kp增大而增大的變化速率越大;水位偏差累積IAE的分布規律與前者相似,只是在積分參數Ki>15 的區域內,IAE隨Ki增加而上升的趨勢相對不明顯。值得注意的是,當Ki=1 與Ki=2 時,由于積分參數過小,PI 控制的積分環節反饋調整量過小,系統無法回到目標穩定狀態,因此,此區域為無效區域,其值不參與討論。流量變化累積IAQ和閘門動作累積IAW由于受過閘流量計算式的關聯,有較強相關性,表現在兩者的分布極其相似,都與積分參數Ki和比例參數Kp的呈正相關,且相關性也類似。穩定過渡時間Tstable隨積分參數Ki有先減后增的趨勢,但增長部分不如減小部分的趨勢明顯,隨比例參數Kp的增大而增大,在Ki=10~15,Kp=1~5的部分范圍內,Tstable最小,穩定在40 min左右,形成矮平臺。綜合指標雷達圖面積ARD隨積分參數Ki的增大而先減后增,隨比例參數Kp的增大而增大。

初次優化得到的參數最優組合為:Kp=1,Ki=9,在此參數組合下,綜合指標ARD最小,其雷達圖面積為0.235,渠道系統在供水工況過渡過程中的最大水位偏差為0.028 8 m,水位偏差累積量為0.192 m,流量變化累積量為0.923 m3∕s,閘門動作累積量為1.06 m,穩定過渡時間為42 min。

二次優化網格范圍設置為:Kp=0.1~2.0,Ki=8.0~10.0,共400個參數組合。得到各性能指標和綜合指標ARD分布情況見圖7。二次優化得到的參數最優組合為:Kp=0.2,Ki=8.2,在此參數組合下,綜合指標ARD最小,其雷達圖面積為0.224,略小于初次優化所得值。渠道系統在供水工況過渡過程中的最大水位偏差為0.028 9 m,水位偏差累積量為0.196 m,流量變化累積量為0.772 m3∕s,閘門動作累積量為0.952 m,前兩者略大于初次優化值,后兩者略小于初次優化值。穩定過渡時間為42 min,與初次優化結果一致。

圖7 二次優化各指標分布圖Fig.7 Indicator distribution of secend optimization

另外,與初次優化結果相比,由于二次優化的參數范圍縮小了10倍,導致所有指標的變化范圍均有了大幅度的縮小。可以注意到,最大水位偏差分布的極差在0.001 m 以內,水位偏差累積的極差在0.02 m以內,流量變化累積的極差在0.4 m3∕s以內,閘門動作累積變化范圍在0.3 m 以內,尤其是穩定過渡時間處于平臺期,均為42 min。綜合指標的變化范圍為0.224~0.226,極差在0.002 以內。這說明了二次優化網格已經足夠小,其控制性能指標再減小的潛力已經太小,沒有必要進行更加精細的參數優化。

以Kp=0.2,Ki=8.2 為最終優化的控制參數,可以得到渠道系統在供水工況下的下游水位、過閘流量和閘門開度的變化過程如圖8~圖10 所示。該系統在較短的時間內,付出較少的閘門動作和流量變化量,將水位控制回復到目標值附近,且在此過程中的水位偏差也較小,符合管理者的預期控制效果。注意上述尋優過程均基于仿真計算,其控制效果的進一步提升需要結合更詳細的渠道參數,通過多次系統調試運行,才能得到更符合渠系實際情況的控制參數,以達到更優的控制效果[20]。

圖8 下游水位偏差變化過程Fig.8 Variation process of downstream water level deviation

圖10 過閘流量變化過程Fig.10 Gate flow change process

3 結 論

本文針對供水渠道系統,選取了最大水位偏差、水位偏差累積、流量變化累積、閘門動作累積和穩定過渡時間5個優化目標,用初篩樣本平均值來改進綜合性能雷達圖,并提出了以該雷達圖面積作為綜合指標。另外,以漳河灌區第四干渠典型供水工況為例,驗證了該綜合指標作為PI 控制中參數優化目標的有效性。主要結論如下:

(1)用樣本平均值改進了雷達圖的繪制方法,相比于用樣本最大值進行標準化處理,改進后的方法在代表性、可靠性和結果呈現3個方面都具有優勢,能夠全面、綜合地反映控制性能的優劣。

圖9 閘門開度變化過程Fig.9 Gate opening change process

(2)提出的雷達圖面積綜合指標,其計算式包含加法和乘法運算,能夠改進簡單直接的加和、乘積作為綜合指標的缺點,兼顧不同數值大小和變異范圍的各指標的影響,能夠合理地評價綜合性能,可為其他供水渠道的控制和管理提供參考。

(3)以漳河灌區四干渠供水工況為例,應用了提出的綜合性能指標雷達圖面積,并經兩次優化后得到參數組合Kp=0.2,Ki=8.2,能使其綜合性能已接近最優。此時雷達圖面積最小,為0.224,最大水位偏差為0.028 9 m,穩定過渡時間為42 min,其水位、流量和閘門開度在過渡過程中的累積量均較小,符合控制預期效果。

猜你喜歡
優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
PEMFC流道的多目標優化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
圍繞“地、業、人”優化產業扶貧
今日農業(2020年16期)2020-12-14 15:04:59
事業單位中固定資產會計處理的優化
消費導刊(2018年8期)2018-05-25 13:20:08
4K HDR性能大幅度優化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 91精品人妻一区二区| 国产成人AV男人的天堂| jizz在线观看| 国内黄色精品| 专干老肥熟女视频网站| 五月天综合网亚洲综合天堂网| 成人日韩精品| 国产成人无码综合亚洲日韩不卡| 日韩精品亚洲人旧成在线| 东京热高清无码精品| 一级毛片网| 新SSS无码手机在线观看| 国产自无码视频在线观看| 色综合天天娱乐综合网| 亚洲第一黄色网| 国产精品99久久久| 国产精品网拍在线| 五月天久久综合| 欧美午夜网| 九九九久久国产精品| 另类欧美日韩| 国产无码高清视频不卡| 青青青国产免费线在| 成年女人a毛片免费视频| 中国国产A一级毛片| 亚洲一区网站| 欧洲亚洲一区| 亚洲婷婷六月| 国产亚洲视频中文字幕视频| 国产成人一区在线播放| 91丝袜在线观看| 色呦呦手机在线精品| 97狠狠操| 在线日韩日本国产亚洲| 中文字幕在线看视频一区二区三区| 青青草原国产免费av观看| 老司机午夜精品网站在线观看| 制服丝袜 91视频| 免费看的一级毛片| 国产精品亚洲αv天堂无码| 综合色88| AV不卡国产在线观看| 重口调教一区二区视频| 国产高清毛片| 一区二区日韩国产精久久| 国产成人精品日本亚洲| 国模私拍一区二区三区| 久久精品国产精品一区二区| 国产网站一区二区三区| 国产精品网拍在线| AV无码国产在线看岛国岛| 一区二区三区精品视频在线观看| 亚洲日韩AV无码一区二区三区人| 人妻91无码色偷偷色噜噜噜| 欧美成a人片在线观看| 亚洲第一成年人网站| 精品自窥自偷在线看| 国产色婷婷视频在线观看| 亚洲欧美一区在线| 热re99久久精品国99热| 国模沟沟一区二区三区| 露脸真实国语乱在线观看| 免费观看精品视频999| 2020精品极品国产色在线观看| 美女免费黄网站| 91网站国产| 亚洲无码精彩视频在线观看 | 亚洲无码高清免费视频亚洲| 国产女人在线| www亚洲精品| 99re这里只有国产中文精品国产精品 | 亚洲欧美日韩动漫| 久青草免费在线视频| 在线观看亚洲人成网站| 国产精品欧美激情| 亚洲三级a| 极品国产在线| 伊人久久精品无码麻豆精品 | 国产精品毛片在线直播完整版| 欧美国产日产一区二区| 国产精品白浆无码流出在线看| 亚洲区视频在线观看|