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

基于協整分析的供熱管道泄漏定位方法

2025-06-26 00:00:00石光輝陳鵬申鵬飛鄭治中王鑄李建剛劉帥張啟亮田雯雯
中國測試 2025年6期
關鍵詞:模型

中圖分類號:TB9;TU995.3 文獻標志碼:A文章編號:1674-5124(2025)06-0067-08

A leakage location method for heating pipeline based on cointegration analysis

SHI Guanghui1,CHEN Peng1,SHEN Pengfei1,ZHENG Zhizhong1,WANG Zhu, LI Jiangang1,LIU Shuai2, ZHANG Qiliang2, TIAN Wenwen2 (1.Taiyuan Heating Group Co.,Ltd.,Taiyuan O30oo0,China; 2.SchoolofControland ComputerEngineering,North China Electric Power University,Baoding O710o0, China)

Abstract: To solve the problem of traditional leak localization methods being unable to cope with the nonstationary characteristics of signals and eliminate the impact of changes in operating conditions on algorithms, a heating pipeline leak localization method based on cointegration and vector error corrction models is proposed.Firstly,perform stationarity tests on the pressure signals of each station to obtain the single integer order of the variables; Then,by conducting cointegration analysis on the signal sequence,calculating the cointegration coeficient,constructing a vector error correction model,and calculating the predicted residual of the model; Propose a threshold setting method suitable for variable operating conditions during leak detection, construct dynamic thresholds to extract leak time points,and accurately locate leak points based on the delay time of each station. Through actual pipeline network data testing experiments,the results show that the proposed algorithm can simultaneously analyze the pressure at multiple sites, effectively suppress the response caused by changes in working conditions,accurately extract the corresponding time of pressure sudden changes,and control the error between theactual delay time and the theoretical delaytime within O.77 seconds. The leakage point positioning error can be controlled within 862 meters, which has practical application value. Keywords: heating pipe network; non stationary; cointegration analysis; vector error correction model; leakage location

0 引言

集中供熱管網是現代城市的基礎設施之一,在其使用過程中,由于自然環境或人為損壞等原因導致管路焊縫缺陷、腐蝕和老化[1,使得泄漏事故頻繁發生,造成大量經濟和資源損失。因此開展供熱管道泄漏檢測與定位研究具有重要意義。

目前常用的管道泄漏檢測方法主要通過對管內壓力、流量等數據進行處理[2,提取泄漏特征實現泄漏檢測與定位,結果具有較高準確性。其中,基于壓力信號的檢測方法可靠性較高,且檢測距離較遠,更適用于實際工程應用,是目前應用最為廣泛的檢測方法[3]。文獻[4]基于測試站點所測壓力數據,對其進行小波變換數據預處理,以探索負壓波方法應用于供熱管網中的可行性。文獻[5]利用一維卷積神經網絡算法建立供熱管道泄漏診斷模型判斷管道是否發生泄漏。文獻[6]綜合考慮管道流量和壓力數據特點,建立了多個管道泄漏系數預測模型,實現管道泄漏系數預測、泄漏量估計。上述方法雖能在一定程度上實現管道泄漏判斷,但大多假定變量為平穩變量,無法準確描述非平穩變量之間的關系,而在實際工業生產過程中壓力信號由于工況變化、設備劣化等情況[]極易出現非平穩特性,因此有必要尋找一種新的泄漏檢測方法,消除數據中的非平穩特性,更準確地實現檢測與定位。

協整理論[8-9]是由Engle和Granger于20世紀80年代提出的,其在經濟學領域作為一種非平穩信號處理方法被廣泛使用。如果多項非平穩信號間存在協整關系,對其進行協整計算可去除信號中所含非線性趨勢和環境等因素的影響,從而得到一組協整殘差,該殘差代表了各項信號間的長期動態均衡關系,通過判斷殘差是否平穩來說明被監測對象的狀態[10]。近幾年協整理論被擴展到工程問題中,文獻[11]通過采用協整分析方法得到攜帶有風力發電機運行狀態信息的協整殘差,有效降低了環境和運行因素對SCADA數據的影響,提高了風力發電機運行狀態識別精度。

綜上,本文提出一種基于協整分析和向量誤差修正模型的供熱管道泄漏定位方法,該方法可有效應對管道壓力信號非平穩特性,同時考慮了多變量間可能存在的長期相互依賴關系,利用協整分析消除變量間的共同趨勢,有效抑制外部因素的影響。

1基礎理論

1.1 協整理論

協整理論能夠挖掘兩組或多組非平穩數據間的長期均衡關系。通過協整計算可以得到這些非平穩時間序列的某種線性組合,使原始非平穩序列轉化為一個平穩時間序列。常見的協整關系分析方法有Engle-Granger和Johasen:Engle-Granger方法主要用于分析兩變量間的單一協整關系;而對多變量、多協整關系進行檢驗通常選用基于向量自回歸模型的Johasen方法。

對于 k 維向量 yt=[y1t,y2t,…,ykt] ,若 yt 中每個變量的單整階數均為 d ,且存在非零向量 α 使αTyt~I(d-b) 0t 的各變量間為 (d,b) 階協整,即 yt~CI(d,b) ,其中 α 是協整向量。

1.2 向量誤差修正模型

向量誤差修正模型作為一種時間序列預測模型,同時考慮了長期協整關系與短期誤差修正[12]。通過對變量進行短期誤差修正可以使其保持長期均衡狀態,協整約束的引入提高了預測模型的穩定性,能夠有效應對非平穩變量建模時可能存在的偽回歸問題。

向量自回歸(vectorautoregression,VAR)模型是向量誤差修正模型的構造基礎。設 k 維時間序列向量 yt=[y1t,y2t,…,ykt] , t=1,2,…,T ,對 yt 各變量有yit~I(1)( (一階單整序列), i=1,2,…,k ,設 ?yt 不存在外生變量影響,構建 p 階VAR模型如下:

yt=A1yt-1+A2yt-2+…+Apyt-p+ut

式中: yt-p (2 模型變量的 p 階滯后變量;p 滯后階數;Ap -系數矩陣;ut 隨機擾動項。

對式(1)進行差分處理得:

式中, 為單位矩陣。

檢驗 yt 各分量間的協整關系,若 yt 各分量間存在協整關系,則有 πyt-1~I(0) (平穩序列),即有向量誤差修正模型:

et-1yt-1

式中: Φl,β ——系數矩陣;et-1 一 -誤差修正項。

在式(5的基礎上,可得模型的預測結果為:

yt=yt-1+Δy

則模型預測殘差為:

1.3泄漏檢測與定位原理

1.3.1 泄漏檢測原理

當管道破損發生泄漏時,泄漏點處的流體存在突然損失的現象,從而引起管內泄漏點附近局部流體密度瞬時減小,造成該點所在位置壓力突降,在泄漏點處產生了一個能夠以某種速度向外傳播的負壓波[13]

因此對管道正常使用時各站點采集到的壓力信號進行協整分析,隨后構建向量誤差修正預測模型,當管道出現泄漏時,模型預測值與實際值之間的殘差幅值及波動程度均明顯增大,可通過設定的閾值檢測管道是否泄漏。

但現有協整分析模型在閥值設定時主要基于統計過程控制(statistical process control, SPC)方法[14],難以有效應對多工況問題,易造成故障誤報,降低監測模型的準確性。因此本文在泄漏檢測時提出一種適用于變工況的閾值設定方法,構建動態閾值并實現泄漏時間點提取。

假設第i個窗口殘差數據的均值記為 μi ,標準差記為 si 。則該窗口所對應的閾值上下限值為:

UCLii+ksi

LCLii-ksi

式中, k 為閾值系數,在本文中取 k=2.5 。

由此可計算每個窗口的閥值上下限值,得到動態閥值曲線,進而根據判定準則來對壓力波信號實時監控。

1.3.2 泄漏定位原理

通過在管道多個地方設置壓力變送器,實時監測管道內流體壓力,通過所建立的協整分析模型對接收到的壓力信號進行分析,提取水擊波產生時間,水擊波在到達各個壓力變送器的時間點不同,根據其到達時間節點差以及水擊波在管內的傳播速度就能夠推算得到泄漏點的具體位置。泄漏點定位原理簡圖如圖1所示。

圖1泄漏點定位原理圖

圖中 L 為管道長度,距離首站X處發生泄漏, ν 為負壓波在介質中的傳播速度,首末兩站點檢測到負壓波信號的時間差為 ΔT ,即可由下式求得泄漏點距首站的距離 X

X=(L+νΔT)/2

其中,負壓波在管道內的傳播速度 u 可用下式計算:

式中: u(t) —流體溫度為 T 時壓力波在管道中的傳播速度, m/s K(t) ——流體溫度為 T 時流體體積彈性模量,Pa ρ(t) 1 -流體溫度為 T 時的密度, kg/m3

E 輸送管道管壁的彈性模量,Pa;

D 管徑,m;

e 管道壁厚, m 。

u(t),K(t),ρ(t) 均受管道內流體溫度的影響,因此管道內流體的溫度在較大程度上直接影響了管道內負壓波的傳播速度,較為精確的流體溫度能夠大幅度提高泄漏點定位精度。

1.4基于協整和向量誤差修正模型的管道泄漏定位方法

基于協整和向量誤差修正模型的供熱管道泄漏定位步驟如下。

1)構建訓練集。在管道正常使用情況下,以各站點壓力信號采集器收集的壓力數據構建訓練集,得到正常運行狀態下變量的時間序列 {yj},j= 1,2,…J 。

2)檢驗變量平穩性。使用增廣迪基-富勒(argumentDickey-Fuller,ADF)檢驗[I5]對各時間序列做平穩性檢驗,得到相應的單整階數,判斷序列是否為一階單整序列。

3)變量的Johansen協整檢驗。確定變量的協整關系數目r。

4)構建向量誤差修正模型。計算協整系數矩陣 β 、調整系數矩陣i等模型參數。

5)構建測試集。采集管道處于正常使用狀態及泄漏使用狀態時的數據,構建測試集。

6預測。將測試集數據輸入已建立的向量誤差修正模型,得到各壓力信號的預測值及模型預測殘差。

7)設定動態預警閾值。利用所提閾值設定方法,設定多工況預警閾值。

8)泄漏檢測。觀測模型預測殘差是否超過設定的閾值,判定管道是否泄漏

9)泄漏點定位。依據各信號監測熱力站與實驗熱力站間的實際距離及管道內壓力信號傳播速度,計算壓力波從實驗熱力站到監測熱力站的理論延遲時間,將實際延遲時間與理論延遲時間作差,進而結合速度求得實際距離誤差。

2 管道泄漏定位

2.1供熱管道結構與實驗內容

本文以太古一級網217#、221#、222#、234#、243#共五座熱力站為測試對象進行相應實驗,其位置分布情況如圖2所示1。為模擬管道泄漏的突發情況,選取217#熱力站作為實驗對象,通過調整閥門開度模擬管道泄漏時壓力波的變化。為采集實驗過程中的壓力信號、實時監測壓力變化情況,為每座熱力站配備了一部頻率為 20Hz 的壓力信號采集器。由于實驗過程中全網平衡及附近補水會對壓力信號產生較大程度的干擾,所以停止周邊熱力站補水并退出全網平衡。

圖2實驗熱力站分布圖

根據圖2所示站點分布圖,測量可以得到 217# 熱力站與其余各監測站點之間的距離,然后結合式(12)計算出的管道內壓力波傳播速度即可得到壓力波從217#站點傳播到其余各熱力站的理論延遲時間,結果如表1所示。

表1監測熱力站與實驗熱力站之間的實際距離與延遲時間

具體實驗步驟如下:將五臺壓力采集器分別與相應試驗熱力站一次網就地壓力表接口連接,并進行上電、聯網校時操作。對217#熱力站站內進行泄水實驗,將供水進站閥門調至零,出站閥門開度關至較小位置,通過調節除污器DN125排污球閥開度來模擬管道泄漏狀態。用壓力采集器進行波形采集,共計進行4次泄放實驗,具體操作如下所示:

1)泄放球閥開度調至 100% ,持續 3min 后關閉。

2)泄放球閥開度調至 25% ,持續 2min 0

3)泄放球閥開度由 25% 調至 50% ,持續 2min 。

4)泄放球閥開度由 50% 調至 75% ,持續 2min 后關閉。

實驗監測到的壓力波動曲線如圖3和圖4所示。

圖3217#熱力站壓力波動
圖4信號監測熱力站壓力波動

由圖3和圖4可知,217#熱力站到各站的距離遠近不一,所以各個站點采集到壓力波動的時間節點不同,并且由于壓力波經過的管道長度不同、分支數不同,導致壓力衰減程度也有所不同。

將5個站點采集數據納入一個19 282×5 的矩陣 y=[y1,y2,…,y5] 中,其中 y1,y2,…,y5 每列數據分別為 217#,221#,222#,234#,243#* 站點所采集的壓力數據,使用0-1歸一化方法對原始數據進行處理。令 1~1900 組數據為訓練集用于建立預測模型, 1901~ 19282為測試集用于測試模型性能。

2.2向量誤差修正模型的建立

基于訓練集建立向量誤差修正模型,首先分別檢驗5列壓力數據是否存在單位根,繼而對它們之間的協整關系進行分析,最后建立向量誤差修正模型。

1)ADF單位根檢驗。根據協整理論,要建立向量誤差修正模型,用于建模的各時間序列均應為一階單整序列,否則模型可能存在偽回歸現象。因此選用ADF方法來對5組壓力數據進行單位根檢驗。若原序列無單位根,則該序列是平穩時間序列;否則原序列不平穩,此時須進一步對其差分做單位根檢驗,當差分序列是平穩序列時原序列即為一階單整序列。檢驗結果見表2。

表2中 y1,…,y5 分別為217#、221#、222#、234#、243#站點所采集的壓力數據, Δ 表示原變量的一階差分。由表2可知,5組數據原始序列均為非平穩序列,而相應一階差分序列是平穩的,即5組壓力數據均滿足一階單整要求,可進行協整分析。

表2各變量ADF檢驗結果

2)變量的Johansen協整檢驗。通過該檢驗確定5組變量間是否存在協整關系以及協整關系的數目。Johansen協整檢驗使用假設檢驗的方式,最多存在 r 個協整關系為其原假設 H0 。 r 由0開始逐漸增加,依次檢驗直到接受原假設為止。檢驗結果見表3。

表3Johansen檢驗結果

由表3可知當r為4時檢驗統計量小于臨界值,此時原假設成立,可得5組變量之間存在4個協整關系。

3)構建向量誤差修正模型。模型數學表達如式(5)所示。計算得到由4列協整向量組成的協整系數矩陣為:

2.3 泄漏檢測

基于測試集數據,利用所建模型對各站點壓力進行預測,閾值求取時滑動窗口長度設為2000,預測殘差如圖5所示。

圖5模型預測殘差曲線

由圖5可知,基于動態閾值各站點均能有效檢測到四次管道壓力變化情況,其中在約 12000~ 16000的數據區間中,各站點殘差均較大,處于閾值限邊緣,是因為該段時間內泄放球閥開度由 25% 調至 50% ,持續 2min 又由 50% 調至 75% ,管道工況變化過快也過于頻繁,導致數據波動較大,因此該區間內一次閥門變化過程算法并未檢測出來。

依據管道泄漏評判準則,超過閾值時所對應的時間即為潛在的時間指標,根據圖5所示結果提取壓力變化時間節點如表4所示。

將表4中除217#熱力站以外站點提取到的泄漏時間與相應的217#熱力站提取時間作差,得到每組實驗對應的實際延遲時間,并將其與表1中的理論延遲時間進行對比,得到相應的時間誤差如表5所示。

表4泄漏信號提取時間

由表5可知,監測熱力站與泄漏熱力站之間的實際延遲時間與理論延遲時間誤差最大僅為 0.762s 因此該方法可以有效提取泄漏時間點,具有良好的泄漏檢測效果。

表5監測熱力站與實驗熱力站之間的實際延遲時間與理論延遲時間誤差

為驗證動態預警閾值的有效性,將所提方法與基于統計過程控制的固定閾值進行比較。若隨機變量 X~N(μ,σ2) ,由正態分布相關理論可知隨機變量X 位于 (μ-2.5σ,μ+2.5σ) 區間內的概率為:

P(μ-2.5σ

當 X 持續超出式(14)所給區間,認為此時過程出現故障。因此基于訓練數據預測殘差,求取殘差均值 ?μ 及殘差標準差 σ ,以 μ+2.5σ 為上限, μ-2.5σ 為下限。帶有固定閾值的殘差預警圖如圖6所示。

圖6帶有固定閾值的殘差預警圖

對比圖5和圖6可以看到,隨著管道狀態的變化,動態閾值能夠實時更新,識別泄漏突變點及后續管道處于穩定運行時的狀態,從而有效提取泄漏時間點;但固定閾值由管道處于正常狀態下的數據確定,因此當管道發生泄漏后,即使后續管道壓力逐漸趨于平穩,但由于工況發生變化,預測殘差依舊處于閾值限以外,無法提取泄漏時間點,實際應用價值較低。

2.4 泄漏點定位

分別計算所提算法測量距離與實際距離之間的誤差,得到結果如表6所示。

表6監測熱力站與實驗熱力站之間所測距離與實際距離誤差

由表6可知,算法所提取時間與管道內壓力波傳播速度計算得到的站點距離和實際距離之間的最大誤差為 861.51m ,能起到有效的定位作用。

3結束語

本文提出一種基于協整分析和向量誤差修正模型的供熱管網泄漏定位方法。將各站點所采集壓力信號作為模型輸入變量,分析了各變量間協整關系,在協整的基礎上建立了向量誤差修正模型。通過實際管網泄漏實驗數據,驗證了所提方法在泄漏檢測時的有效性,并依據所提取的泄漏發生時刻完成了管道泄漏定位任務。該方法優勢如下:

1)所提方法能夠有效應對壓力信號非平穩特性,通過同時對多組數據進行分析,去除數據中的共同趨勢,降低外界環境與工況變化的影響,準確實現管道狀態監測;2)為實現多工況狀態監測,提高算法準確性,在協整分析之后提出一種動態閾值設定方法,能夠有效實現泄漏時間點提取。3)所提算法能夠有效提取壓力突變點對應的時間,實際延遲時間與理論延遲時間誤差均能控制在0.77s內,將泄漏點定位誤差控制在 862m 以內,具有實際應用價值。

參考文獻

[1]張曼,張立申,王隨林,等.供熱管道泄漏流場/聲源特性及 其變化規律[J].消防科學與技術,2021,40(3):312-318. ZHANGM,ZHANGLS,WANGSL,etal.Characteristics and variation of flowing/acoustic source of leaking heating pipeline[J]. Fire Science and Technology, 2021, 40(3): 312- 318.

[2]HENRIE M, CARPENTERP,NICHOLASRE.Pipeline leak detection handbook[J].Boston: Gulf Professional Publishing, 2016:1-17.

[3]聶維,江竹,劉伯相,等.一維卷積長短期記憶神經網絡的管 道泄漏檢測方法[J].中國農村水利水電,2022(1):147-152. NIE W,JIANG Z,LIU B X,et al. One-dimensional convolutional and long short-term memory neural network method for pipeline leak detection[J]. China Rural Water and Hydropower, 2022(1): 147-152.

[4]石光輝,齊衛雪,陳鵬,等.負壓波與小波分析定位供熱管道 泄漏[J].振動與沖擊,2021,40(14):212-218. SHI G H, QI W X, CHEN P, et al. Experimental study of leakage location in a heating pipeline based on the negative pressurewave and wavelet analysis[J]. Journal of Vibration and Shock,2021,40(14): 212-218.

[5]馬廣興,曲波,常琛,等.基于CNN的供熱管道泄漏識別方 法研究[J].電子測量技術,2022,45(16):34-41. MA G X,QU B, CHANG C,et al.Research on leakage identification method of heating pipeline based on CNN[J]. Electronic Measurement Technology,2022,45(16): 34-41.

[6]馬云路,鄭堅欽,梁永圖.基于特征提取的輸油管道泄漏系 數預測[J].中國安全生產科學技術,2022,18(10):130-135. MAYL,ZHENGJQ,LIANGYT.Prediction onleakage coefficient of oil pipeline based on feature extraction[J]. Journalof Safety Scienceand Technology,2022,18(1O):130- 135.

[7]王印松,劉勇,孫天舒.基于協整系數矩陣的電動調節閥劣 化分析[J].控制工程,2021,28(12):2293-2298. WANG Y S,LIU Y,SUNT S.Degradation analysisof electric control valve based on cointegration coefficient matrix[J].ControlEngineeringofChina,2021,28(12):2293- 2298.

[8]ENGLE R F,GRANGER C W J. Cointegration and errorcorrection:representation,estimationandtesting,[J]. Econometrica, 1987,55(2): 251-276.

[9]JOHANSEN S. Statistical analysis of cointegration vectors[J]. Journal ofEconomicDynamicsamp;Control,1988,12(2-3):231- 254.

[10]張超,武越,鄭思遙.協整分析下的風力發電機狀態監測方 法[J].機械設計與制造,2023,(12):85-87 ZHANG C,WUY,ZHENG SY.Monitoring method ofwind turbine state under cointegration analysis[J].Machinery Designamp; Manufacture,2023,(12):85-87

[11] 汪千程,蘇春,文澤軍.基于協整分析的風力機多工況監測 與故障診斷[J].中國機械工程,2022,33(13):1596-1603. WANG QC, SU C, WEN Z J. Multi-condition monitoring and faultdiagnosis of wind turbines based on cointegration analysis[J]. China Mechanical Engineering, 2022,33(13): 1596-1603.

[12]林愛華,沈利生.向量誤差修正模型中誤差修正系數的辨識 [J].統計與決策,2020,36(8):27-31. LINA H,SHENL S.Identification of error correction coefficient invector error correction model[J].Statisticsand Decision,2020,36(8): 27-31.

[13]程家銘,張漢國.輸油管道負壓波法測漏原理及實現[J].石 油機械,2002(9):28-30. CHENG JM, ZHANG H G. Leak hunting of oil pipeline by negative pressure wave[J]. China Petroleum Machinery, 2002(9):28-30.

[14]劉長良,張書瑤,王梓齊.基于改進KNN回歸算法的風電機 組齒輪箱狀態監測[J].中國測試,2021,47(1):153-159. LIU CL, ZHANG S Y, WANG Z Q. Condition monitoring of windturbine gearbox based on improved KNN regression algorithm[J].ChinaMeasurementamp; Test,2021,47(1):153- 159.

[15] DIAO Y S, CAO Y D,SUN Y T. Structural damage identificationbasedonARmodelcoefficientsand cointegration for Offshore platform under environmental variations[J].EngineeringMechanics,2017,34(2): 179-188.

(編輯:譚玉龍)

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产成熟女人性满足视频| 国产精品一区二区国产主播| 女高中生自慰污污网站| 国产在线精品99一区不卡| 一本一道波多野结衣一区二区| 久久国产成人精品国产成人亚洲 | 亚洲浓毛av| 久久国产精品77777| 久久网综合| 亚洲国产欧美国产综合久久 | 日韩无码黄色| 免费啪啪网址| 亚洲AV无码不卡无码| 国产91视频免费观看| 国产欧美日韩专区发布| 色欲不卡无码一区二区| 农村乱人伦一区二区| 999国产精品| 在线va视频| 国产在线第二页| 国产导航在线| 成年片色大黄全免费网站久久| 麻豆国产在线观看一区二区 | 亚洲一级毛片| 免费一级毛片| 精品少妇人妻av无码久久| 中文字幕一区二区人妻电影| 久久午夜影院| 国产精品 欧美激情 在线播放 | 国产情侣一区| 免费看的一级毛片| 国产99视频在线| 九色视频线上播放| 欧美日韩午夜| 久操中文在线| 日韩精品毛片| 日本欧美一二三区色视频| 国产九九精品视频| 亚洲一区二区在线无码| 99在线视频免费| 伊人91在线| 性色生活片在线观看| 在线观看av永久| 国产福利在线免费| 精品国产香蕉在线播出| 欧美在线视频不卡| 国产亚洲精久久久久久无码AV| 四虎成人精品在永久免费| 日韩精品资源| 欧洲成人在线观看| 欧美另类精品一区二区三区| 欧美中文字幕在线二区| 久久www视频| 美女内射视频WWW网站午夜| 夜夜高潮夜夜爽国产伦精品| 999精品在线视频| 久视频免费精品6| 国产精品久线在线观看| 秘书高跟黑色丝袜国产91在线| 国产情侣一区二区三区| 久久成人18免费| 人妻无码一区二区视频| 久久精品娱乐亚洲领先| 久操线在视频在线观看| 青青青国产视频手机| 久久国产毛片| 精品免费在线视频| 亚洲无码高清视频在线观看| 狠狠做深爱婷婷久久一区| 永久免费无码成人网站| 999福利激情视频| 欧美精品亚洲精品日韩专| 国产激情第一页| 国产精品久久自在自线观看| 99久久精彩视频| 影音先锋亚洲无码| 国产欧美视频综合二区| 国产精品亚洲一区二区三区z| 国产va免费精品| 国产乱子伦视频三区| 亚洲国产无码有码| 九色综合伊人久久富二代|