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

基于多中點弦測法的鋼軌波磨測量不確定度分析

2021-01-09 02:41:08殷華萬靈
鐵道科學與工程學報 2020年12期
關鍵詞:測量

殷華,萬靈

基于多中點弦測法的鋼軌波磨測量不確定度分析

殷華1,萬靈2

(1. 江西農業大學 軟件學院,江西 南昌 330045;2. 江西農業大學 工學院,江西 南昌 330045)

鋼軌波磨的準確測量一直是鐵路工務中的難點,雖然從理論上多中點弦測法已被證實能夠對鋼軌波磨進行測量,但目前缺乏計量方法。為了準確評估多中點弦測法的誤差分布,在無法取得鋼軌波磨真值的前提下,以雙中點弦建立了誤差傳遞模型并采用蒙特卡羅方法計算測量不確定度。研究結果表明:多中點弦測量結果誤差大小與原始波磨幅值呈正相關、與超限波長呈負相關,從整體來看,雙中點弦測量相對誤差約為7%,能夠滿足實際鐵路工務應用。

鋼軌波磨;弦測法;不確定度;誤差分析

波磨病害廣泛存在于各既有線路上,會導致噪聲與振動,嚴重時還會造成機車脫軌,帶來嚴重后果。雖然國內外諸多學者們對鋼軌波磨進行了長時間的研究,從不同方面解釋了其成因,但由于輪軌之間作用較為復雜,目前尚無哪一種理論能夠完整解釋波磨產生的機理,波磨病害沒有辦法避免。因此,對鋼軌表面狀態進行日常巡檢,在保養維護過程中對波磨病害加以關注,以期盡可能早的發現問題、抑制其進一步惡化成為當前鐵路工務養護的重點內容之一[1?3]。由于鐵路自建成起就固定在路基上不可移動,現場環境復雜多變且缺乏恒定的基準,而鋼軌波磨的幅值又通常在1 mm以下,這些因素都給鋼軌波磨的快速測量帶來了困難。因此,長期以來業界主要采用動態的慣性法和靜態的弦測法對鋼軌波磨進行檢測[4?5]。慣性法由加速度傳感器來完成,在測量設備高速行進的過程中對探測到的加速度信號進行二次積分得到波磨幅值;而弦測法采用平直尺(通常為1 m)來完成,以手工測量的方式得到數據后逐段拼接得到測量結果。這2種方法各有優劣:動態的慣性法測量成本高且重復性稍低,靜態的弦測法重復性較好但測量效率低下。為了能夠兼顧檢測效率與精度,殷華等[6?7]對傳統的弦測方法進行改進,提出了基于多中點弦測量原理的鋼軌波磨動態檢測小車,經過理論分析與現場試驗已被證實能夠用于實現對鋼軌波磨進行快速、有效的測量。但由于目前尚無鋼軌波磨弦測的計量標準,鋼軌波磨的真值無法準確得到,并且鋼軌表面狀態的復雜性、小車制造過程中弦中點位置的偏差、弦長的偏差以及傳感器本身的不確定性都會對測量結果造成一定影響。因此,深入分析多中點弦測量結果的不確定度來源,明晰在多個因素的影響下弦測結果的誤差分布,對鋼軌波磨多中點弦測量理論的發展、實際工程應用過程中減少測量結果誤差,得到盡可能準確的波磨測量結果具有重要的理論及現實意義。

1 鋼軌波磨多中點弦測原理

將一根剛性的弦靜止于鋼軌表面,弦的首尾與鋼軌表面接觸,在弦的中點位置安裝一個測量傳感器,采用接觸或非接觸的方法來得到其與鋼軌表面的垂直距離,此即為中點弦測法,其測量值可由式(1)表達:

式中:表示測量弦的長度,1和X分別表示該弦在鋼軌表面起點、終點的位置坐標,2表示安裝在弦中點處的測量傳感器的位置坐標,()軌面不平順函數。

根據文獻[8?9],鋼軌波磨波長通常在1 m以下且具有一定的準周期性,若設其周期為,即()(),則由式(1)可知,中點弦測法所得結果同樣存在周期性。對其進行頻域變換,可得傳遞函數為

為了克服這種缺陷,文獻[6?7]提出了多中點弦測量的思想,根據單一的中點弦測法位置不存在偏差這一特點,設計了一種基于雙中點弦測量原理的波磨檢測小車(如圖1),將多個中點測量弦的結果進行疊加后再進行逆濾波計算得到波磨幅值。

圖1 雙中點弦測小車

此時,弦測結果及其頻域特性可以分別表示為式(4)和式(5):

式中:表示所使用中點弦的數量。按照這種思想,總能找到至少一個合理的弦的數量和長度的組合使得對目標區間內的所有待測波長均不存在無響應的零點。以上述雙中點弦測量小車所用的174 mm和292 mm結構為例,圖2為其幅值增益波長響應曲線:若采用單中點弦結構,不論是174 mm還是292 mm的弦長都存在不止一個無響應的零點;而當采用174+292 mm雙中點弦進行組合測量后,幅值增益曲線不存在零點,且更加平緩。隨后,再利用頻域逆濾波的方法[10],根據該幅值增益曲線計算出波磨的真實幅值。在理想情況下,無論組成測量小車的弦長為多少,只要其幅值增益波長響應曲線不存在無響應的零點,且不考慮振動與系統裝配精度的影響,雙中點弦測量結果與原始波磨幅值對比在理論上誤差最大不會超過0.02 mm。因此,雙中點弦測法理論上能夠得到準確的鋼軌波磨測量 結果。

圖2 雙中點弦測量幅頻特性曲線

2 測量不確定度評定方法

根據前述可知,采用多中點弦測量的方法改善了幅值增益?波長響應曲線,消除了測量時無響應的零點,使弦測法準確測量鋼軌波磨變為了可能。但上述分析均在理想情況下進行的,而在實際工程中,中點弦結構的制造偏差、鋼軌表面狀態的復雜、傳感器本身的測量不確定性及其他相關因素都會對測量結果帶來一定的誤差,在無法得到鋼軌波磨真值的前提下,上述因素究竟會對多中點弦測量造成多大的影響,多中點弦法測量結果的不確定性有多少,這些都是從理論上必須考慮的問題。

通常,要得到一個系統的不確定度主要有2種方法:基于《測量不確定度表示指南》提出的GUM方法和基于蒙特卡洛傳遞分布(MCM)方法[11]。采用GUM方法得到待測對象不確定度需要經過建立準確的模型、利用偏導求靈敏系數等關鍵步驟,但隨著測量系統的越來越復雜,很多時候無法建立統一的數學模型,更無法求得靈敏系數,特別在輸出量概率分布明顯不對稱時,GUM方法還有可能得到不正確的結果。而蒙特卡洛方法是根據待測系統的不確定度來源及分布特點,選用合適的概率分布模型隨機產生模擬輸入的數據代入計算,從而得到分布傳播規律的一種數值方法,適用于具有一個或以上輸入量、單一輸出量的測量模型。其基本步驟為:

1) 根據待測系統的特點,建立測量模型=(1,2,…,X),其中表示系統的輸出,X(=1, 2, …,)表示對輸出有影響的各個輸入量。

2) 根據X的特點,選定合適分布類型

3) 設定參與蒙特卡洛實驗的樣本的大小,選擇每個X的個樣本,代入到前面的測量模型中,可以獲得的多個輸出。

4) 將得到的多個進行排序,得到的分布函數及約定概率下的包含區間[y,y]。

由于多中點弦測法模型復雜,在計算的過程中會涉及到多個影響參數的輸入、頻域變換、逆濾波等過程,無法建立精確統一的數學模型。因此,對其不確定度的評定采用蒙特卡洛方法得到。

3 多中點弦測量不確定度來源

由于已完成了雙中點弦測量樣機的制作且能夠正常獲取數據,故下面以具有5個測頭的雙中點弦測量小車為基礎進行分析。測量小車采用5個基恩士公司推出的CMOS點激光位移傳感器IL-100構成雙中點測量弦,其發射波長為655 nm的紅色半導體激光,線性度最高可達±0.025%,輸出1~5 V的模擬電壓,對應距離測量范圍為?5 mm~+5 mm。模數轉換器選用ADI公司的AD7606芯片,其具有雙極性16位分辨率和8路信號同步采樣,避免了異步采樣延時誤差的產生。

3.1 數據采集系統誤差分析

雖然傳感器在出廠時已經標注了精度,但是測量信號經過后端調理電路、采樣、量化等過程后都會引入各種誤差,因此需要對整個采集系統的誤差進行重新評定。對安裝在測量小車上的5個傳感器器進行校準,由于波磨幅值較小,校準通過測量5 mm標準量塊進行,得到如表1所示數據。

表1 對5 mm量塊的測量結果

通過計算可以得到5個測頭10個樣本的均值與方差??紤]到傳感器誤差通常滿足正態分布且樣本數量較少,故可以使用式(6)中的2個公式來分別估計每個測頭均值及方差區間,并以此為依據來作為實際測量時的誤差。

3.2 多中點弦結構誤差影響

從前述多中點弦測量原理可知,要利用多中點弦結構對鋼軌波磨進行測量必須滿足以下2個條件:弦的長度須與理論設計值一致,每根測量弦必須嚴格為中點弦,否則傳遞函數會受到影響,繼而在頻域逆濾波的過程中產生誤差。但在雙中點弦測量設備制造和安裝過程中,由于機械結構加工精度和傳感器手動安裝限制,雙測量弦的長度、中點測頭的位置都會有所偏差,那么這些誤差需要控制在什么樣的一個范圍內才會對最終結果影響最小,也是必須考慮得問題。

3.2.1 弦長偏差的影響

式(8)中:表示待測波磨的目標波長。從上式可以看出,誤差大小受測量弦長、待測目標波長及弦長誤差1,2影響。以樣機中采用的174 mm和292 mm弦長為例,給出1,2在不同目標波長下的誤差Δ變化情況(圖3),考慮到在實際機械加工與安裝過程中,不可能出現5 mm以上的誤差,因此1,2取值限定在5 mm之內。

(a) 174 mm弦長;(b) 292 mm弦長

從圖3可以看出,不論是174 mm的弦還是292 mm的弦只要在制造和裝配的過程中出現了少許誤差,其對幅值增益?波長響應曲線的影響是較大的。特別是對于200 mm波長以下的波磨,其最大的幅值增益誤差可超過0.4,但對于400 mm以上的波磨,幅值增益誤差較小。

3.2.2 弦中點偏差的影響

除了弦長誤差外,弦中點誤差同樣會對最終測量結果造成影響,因此,在制造裝配的過程中同樣要對弦中點的位置進行嚴格限定。根據式(4)~(5),一旦中點測頭的實際位置偏離了弦的中點,那么中點弦就變成了偏弦,其不再滿足中點弦測法中相位不發生偏移的特性。此時,雙弦測法的傳遞函數將變為式(9):

式中:原始的設計弦長分別是1,2,但中點位置發生了偏移,形成了1,2的偏弦結構,且≠,≠。同樣以174 mm+292 mm的雙弦組合為例,討論中點偏差在5 mm之內時的測量結果誤差如圖4。

從圖(4)中分析可知,當中點偏差達到2 mm時,幅值增益偏差就超過0.1而相位偏差接近于0.5弧度,測量結果與實際結果偏差嚴重,因此,弦中點誤差不應該超過2 mm。另外,在理想情況下只要5個測頭均垂直于鋼軌表面可以得到較為理想測量弦長和雙弦結構。但由于安裝過程中會出現誤差且小車在鋼軌表面推行時,軌道不平順會對測量小車的姿態造成一定影響,故此時測頭實際上是垂直于測量弦的。由于鋼軌波磨的幅值通常不會大于5 mm,為了模擬在鋼軌表面走行時測量小車姿態改變帶來的測頭的變化范圍,采用在實驗室中墊入5 mm標準量塊方法,近似可以估算小車各測頭之間距離變化范圍如表2。

(a) 最大幅值增益誤差;(b) 最大相位偏差

表2 5 mm量塊下測量弦變化范圍

考慮到鋼軌表面波磨滿足準正弦特征,其幅值呈均勻分布,因此在小車姿態影響下雙測量弦的長度同樣會符合均勻分布,雙弦的概率密度函數如式(10),并以此作為蒙特卡羅分析的輸入。

3.3 外界溫度及其他影響

考慮到測量小車為鋁合金材質,而我國鐵路軌道分布較為廣泛,當外界環境發生變化時,小車本身的熱脹冷縮同樣會使得弦長發生一定的變換。設小車的原始長度為,則根據熱膨脹計算公式(11)可以計算機出環境溫度每增加1度,小車長度變 化為:

′=×(1+23.8×0.000 001) (11)

另外,還有編碼器引起的采樣誤差、受溫度變化待測鋼軌同樣會發生膨脹等等,但通過理論計算機同樣發現其帶來的變化較小,故此處暫時不予考慮。

4 基于蒙特卡洛的不確定度分析

要對中點弦測結果的不確定度進行準確的評定,模擬生成的軌道與實際軌道的相符程度至關重要,在這方面國內外學者進行了很多研究,產生了諸如二次濾波法、白噪聲濾波法等方法,但這些方法或是針對某種特殊類型不平順進行模擬、又或是計算過程過于復雜不具有通用性,導致不能在實際中快速運用??紤]到鋼軌波磨的嚴重程度采用粗糙度譜來進行評價,Nielsen等學者綜合上述方法,提出利用多個正、余弦波疊加的方法來模擬軌道不平順,其核心是通過采用ISO3095標準中對鋼軌踏面粗糙度的評價方法,構造符合粗糙度譜標準的鋼軌波磨[12]。

設當前在頻帶內隨機生成復合波由下式表示,并據此求得均方根為R

RR之間的關系為,為每個頻率成份的長度,y為第個頻帶內符合ISO標準的幅值,則可由式(13)表示。

將上式簡化可yy之間的關系

由此,只需要求得就可以得到滿足ISO3095中鋼軌踏面粗糙度要求的軌道模擬波形。該算法受到軌道踏面粗糙度評價標準的制約,可以得到滿足同一粗糙度水平下的多個隨機波形,且當需要在某一頻段產生特定的波磨病害時,只需要改變t值即可。利用上述方法構造鋼軌波磨模擬波形,同時按照鐵路工務的實際工況,設定溫度在范圍?20~+60 ℃之間均勻分布,雙中點弦測小車的信號采集系統為正態分布,分布參數在式(7)中隨機選擇,受鋼軌波磨影響小車姿態改變帶來的各傳感器之間的間隔變化為均勻分布。為了更可能的接近真實的鋼軌現場工況,在生成的波形中隨機加入中心波長超限波形并疊加信噪比為20 dB的白噪聲。共生成模擬軌道500 000段,每段取10 000個點,共產生了5×109個測試數據點進行蒙特卡羅分析,由此取得誤差的均值和概率分布。圖5為某段生成的模擬軌道波形及幅值的95%概率區間,從圖中可以看出幅值約為±0.5 mm。

(a) 模擬波形;(b) 模擬波形的95%概率區間

圖5 生成的模擬鋼軌波磨波形

Fig. 5 Waveform of Simulation of track irregularity

由于已經從傳感器誤差、數據采集系統誤差、弦長及弦中點誤差方面對174+292 mm的雙中點弦測量小車進行了不確定度來源分析,確定了誤差的概率分布,故在此同樣利用174+292 mm的雙中點弦結構在模擬軌道上進行測量。圖6為測量后,復原波形與原始波形之差,從圖6(a)中可以看出除首尾部分外,最大相差不超過0.02 mm(這是因為首尾部分在頻域逆濾波時被截斷導致,實際應用中可以忽略),而從圖6(b)中亦可以看出誤差的95%概率區間在±0.02 mm內,相對誤差為4%。

(a) 復原誤差;(b) 復原誤差的95%概率區間

考慮到鋼軌波磨的波長恒定機理,分別設定40~400 mm中單一中心波長超限,超限幅值在信噪比0~150 dB之間隨機選擇,在不確定度來源、概率分布、樣本數量與之前一致的前提下,得到的鋼軌波磨原始幅值與雙中點弦逆濾波結果的標準差及95%概率分布如表3所示。從表中分析可知,受到各種不確定的因素影響,不論是復合波磨波長超限還是單一波磨波長超限,弦測結果會與模擬軌道的波磨真值有所偏差:超限波長越短,相對波長誤差越大,超限波長越長,相對誤差越??;模擬軌道波磨幅值越大,弦測結果偏差越大;而模擬軌道波磨幅值越小,弦測結果偏差越??;但從整體來看,多中點弦測結果相對誤差均值約為7%,該誤差能夠滿足鐵路工務中對鋼軌波磨測量需求。

表3 波長超限原始波磨幅值與弦測小車測量誤差對比

5 結論

1) 測量弦的弦長誤差、弦中點位置偏差都會對鋼軌波磨的測量結果產生影響,同時受到各種不確定的因素影響,不論是復合波磨波長超限還是單一波磨波長超限,弦測結果都會與模擬軌道的波磨真值有所偏差, 偏差越大影響越多。

2) 波磨幅值與弦測結果偏差密切相關;以174+292 mm的雙中點弦結構為例,在實際檢測設備達到文中所列的裝配精度后,通過蒙特卡羅分析得到平均相對誤差約為7%,此時,雙中點弦測法適合鐵路工務中對軌道波磨檢測。

3) 若采用三中點弦及以上的多中點弦結構,測頭數量的增加、測量弦數量的增加使得測量結果受更多因素影響,因此必須按照文中方法對其測量不確定度進行重新評定,確定滿足實際需求的工程裝配精度。

[1] Sato Y, Matsumoto A, Knothe K. Review on rail corrugation studies[J]. Wear, 2002, 253(1): 130?139.

[2] Grassie S L. Rail corrugation: Characteristics, causes, and treatments[J]. Proceedings of the Institution of Mechanical Engineers, Part F: Journal of Rail & Rapid Transit, 2009, 223(6): 581?596.

[3] Vuong T T, Meehan P A, Eadle D T, et al. Investigation of a transitional wear model for wear and wear-type rail corrugation prediction[J]. Wear, 2011, 271(1/2): 287?298.

[4] 姜子清, 司道林, 李偉, 等. 高速鐵路鋼軌波磨研究[J].中國鐵道科學, 2014, 35(4): 9?14. JIANG Ziqing, SI Daolin, LI Wei, et al. On rail corrugation of high speed railway[J]. China Railway Science, 2014, 53(4): 9?14.

[5] 劉伶萍, 杜鶴亭, 楊愛紅. 鋼軌波浪磨耗檢測系統的研究開發[J]. 中國鐵道科學, 2002, 23(6): 65?69. LIU Lingping, DU Heting, YANG Aihong.Development of rail corrugation inspection system[J]. China Railway Science, 2002, 23(6): 65?69.

[6] 殷華, 朱洪濤, 魏暉, 等. 基于中點弦測模型的鋼軌波磨量值估計[J]. 振動、測試與診斷, 2016, 36(5): 954? 959. YIN Hua, ZHU Hongtao, WEI Hui, et al. Discussion on estimate rail corrugation amplitude based upon midpoint chord model[J]. Journal of Vibration, Measurement & Diagnosis, 2016, 36(5): 954?959.

[7] 殷華, 朱洪濤, 王志勇, 等. 基于多弦模型的軌道短波不平順測量研究[J]. 振動與沖擊, 2017, 36(14): 178? 182, 193. YIN Hua, ZHU Hongtao, WANG Zhiyong, et al. Rail short-wave irregularity measurement based upon a multi- midpoint chord model[J]. Journal of Vibration and Shock, 2017, 36(14): 178?182, 193.

[8] 曹亮, 許玉德, 周宇. 城市軌道交通鋼軌波浪形磨耗特征分析[J]. 城市軌道交通研究, 2010(2): 46?48, 52. CAO Liang, XU Yude, ZHOU Yu. Characteristics of rail corrugation in urban mass transit[J].Urban Mass Transit, 2010(2): 46?48, 52.

[9] 沈鋼, 張學華, 郭滿鴻. 地鐵曲線鋼軌波浪型磨耗的測量分析[J]. 城市軌道交通研究, 2011(4): 53?54, 58. SHEN Gang, ZHANG Xuehua, GUO Manhong. Measurement and analysis of rail corrugation on curved track of metro system[J].Urban Mass Transit, 2011(4): 53?54, 58.

[10] 程櫻, 許玉德, 周宇. 三點偏弦法復原軌面不平順波形的理論及研究[J]. 華東交通大學學報, 2011(1): 42?46. CHENG Ying, XU Yude, ZHOU Yu. Theory and research of asymmetrical chord offset method of restoring a wave form of track irregularity[J]. Journal of East China Jiaotong University, 2011(1): 42?46.

[11] 劉園園, 楊健, 趙希勇, 等. GUM法和MCM法評定測量不確定度對比分析[J]. 計量學報, 2018, 39(1): 135? 139. LIU Yuanyuan, YANG Jian, ZHAO Xiyong, et al. Comparative analysis of uncertainty measment evaluation with GUM and MCM[J]. Acat Metrologica Sinica, 2018, 39(1): 135?139.

[12] Hiensch M, Nielsen J C O, Verheijen E. Rail corrugation in the Netherlands—measurements and simulations[J]. Wear, 2002, 253(12): 140?149.

Uncertainty evaluation for rail corrugation measurement based upon multi-chord method

YIN Hua1, WAN Ling2

(1. School of Software, Jiangxi Agricultural University, Nanchang 330045, China;2. College of Engineering, Jiangxi Agricultural University, Nanchang 330045, China)

An effective method to measure the rail corrugation accurately on railway engineering has been quite a task in railway maintenance management for dozens of years. Although multi-chord based method has been proved can be measured rail corrugation easily, it lacks metrological method. Due to the true value of rail corrugation cannot be obtained, the error transfer model was established as well as uncertainty was calculated by Monte Carlo Method(MCM)for evaluating the error distribution of the multi-midpoint chord method. The data simulation and test results show that the error of multi-chord based method is positively correlated with the amplitude and negatively correlated with the wavelength, the relative error of double chord method is approximately 7%.Thus, it is appropriate for engineering application.

rail corrugation; chord measurement method; uncertainty; error analysis

U216.3;TH17

A

1672 ? 7029(2020)12 ? 3036 ? 09

10.19713/j.cnki.43?1423/u.T20200110

2020?02?14

江西省自然科學基金資助項目(20192BAB206032,20202BABL214041);國家自然科學基金地區科學基金資助項目(51468042)

萬靈(1986?),女,江西吉安人,講師,博士,從事結構智能監測、災害預警等方向的研究;E?mail:wanlingstar@126.com

(編輯 涂鵬)

猜你喜歡
測量
測量重量,測量長度……
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
測量的樂趣
二十四節氣簡易測量
日出日落的觀察與測量
滑動摩擦力的測量與計算
測量
測量水的多少……
主站蜘蛛池模板: 中文字幕免费在线视频| 亚洲男人在线天堂| 美女被操黄色视频网站| 国产新AV天堂| 欧美精品一区二区三区中文字幕| 国产在线自乱拍播放| 东京热av无码电影一区二区| 二级特黄绝大片免费视频大片| 波多野结衣在线一区二区| 成年人国产网站| 成人一区专区在线观看| 国产精品亚洲片在线va| 日本手机在线视频| 九九热在线视频| 亚洲性一区| 国产不卡在线看| 麻豆精品在线播放| 99精品久久精品| 国产精品私拍99pans大尺度| 性色生活片在线观看| 四虎影视国产精品| 国产凹凸一区在线观看视频| 国产a v无码专区亚洲av| 日韩国产一区二区三区无码| 亚洲无码视频喷水| 18禁黄无遮挡免费动漫网站| 亚洲成人www| 国产在线啪| 日韩精品视频久久| 中文字幕自拍偷拍| 影音先锋丝袜制服| 在线欧美a| 91久久夜色精品| 婷婷开心中文字幕| 国产精品亚洲一区二区在线观看| 72种姿势欧美久久久大黄蕉| 久久婷婷国产综合尤物精品| 日韩精品久久无码中文字幕色欲| 国产婬乱a一级毛片多女| 99在线视频免费观看| 久久婷婷五月综合色一区二区| 40岁成熟女人牲交片免费| 国产人免费人成免费视频| 亚洲人成网址| 精品人妻AV区| 2021国产精品自产拍在线| 99在线视频网站| 精品综合久久久久久97超人| 黄色三级网站免费| 男人的天堂久久精品激情| 国产成人三级| 99热这里只有成人精品国产| 热99精品视频| 国产免费久久精品99re丫丫一| 91福利免费| 国产成人精品高清在线| 国产精品太粉嫩高中在线观看| 久久黄色免费电影| 毛片在线看网站| 国产欧美精品专区一区二区| 综合色88| 亚洲91在线精品| 欧美精品色视频| 日韩中文无码av超清| 亚洲天堂在线免费| 老司国产精品视频| 呦女亚洲一区精品| 亚洲黄色视频在线观看一区| 97久久超碰极品视觉盛宴| 亚洲综合在线网| 亚洲精品午夜无码电影网| 国产91九色在线播放| 国产黄在线观看| 一级毛片免费不卡在线| 精品久久国产综合精麻豆| 国产精品亚洲日韩AⅤ在线观看| 久久久久久久久亚洲精品| 亚洲天堂视频网站| 四虎精品国产AV二区| 欧美www在线观看| 怡春院欧美一区二区三区免费| 精品欧美一区二区三区在线|