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

溫度作用下斜拉橋撓度的時間多尺度分析

2020-03-06 05:34:02趙丹陽
哈爾濱工業大學學報 2020年3期
關鍵詞:大氣效應信號

黃 僑, 趙丹陽, 任 遠, 許 翔

(東南大學 交通學院, 南京 211189)

撓度是橋梁健康監測系統實時監測的關鍵指標之一[1],也是評定橋梁結構安全與否的關鍵參數[2-3]. 長期監測得到的撓度數據通常是由車輛及環境因素(溫度、風)等綜合作用下的總效應[4]. 為了更好地掌握橋梁的運營狀況,深入分析撓度檢測數據中的各項響應分量是十分必要的[5]. 研究表明,環境是影響主梁撓度變化的主要因素[6],結構的損傷信號往往會被其“淹沒”,以致很難直接根據實測信號對結構的安全狀態作出合理的評價[7].

國內外已有學者針對溫度作用對橋梁撓度的影響進行過相關研究. 文獻[2]利用長期監測數據分析了鋼箱梁斜拉橋撓度隨溫度的變化特征,并得到了恒載撓度評估指標的基準值. 文獻[8]用中心移動平均法分離撓度的溫度效應. 此外,也有部分學者對日周期上溫度效應進行過研究. 文獻[9]建立有限元模型對青馬大橋的溫度場及撓度等響應進行模擬,發現太陽輻射對主梁豎向撓度影響明顯. 文獻[4]基于日溫差和日溫差效應的線性關系,提出了分離監測信號溫度效應的方法. 文獻[10-12]分析了長期監測數據,認為撓度與溫度并非一般的線性關系.

以上研究認為溫度作用對撓度的影響較大,但大多針對大氣溫度與撓度的關系進行分析,然而日周期中的主梁溫度梯度、太陽輻射等因素對撓度的影響不可忽略. 此外,基于實測數據的不同時間尺度下撓度溫度效應的研究尚欠缺,因此撓度的溫度效應仍需更深入的分析,從而達到為橋梁評估結果的可靠性判斷提供依據的目的. 本文基于有效的長期監測數據,通過多尺度分析,對撓度響應信號在溫度作用的頻率處進行重構,在不同時間尺度上實現撓度的溫度效應分離. 結合日溫度作用頻率處撓度響應的特點,同時考慮太陽輻射的影響,考察日溫度效應與大氣溫度的關系,實現對日溫度效應的預測.

1 撓度響應中溫度效應分離方法

1.1 撓度響應的多時間尺度特征

橋梁在正常運營過程中由長期監測系統采集到的主梁撓度數據是在溫度、車輛、風等共同作用下的綜合響應,其中影響最大的是溫度和汽車作用[2],而這兩種作用的作用周期和影響因素明顯不同. 日溫度作用的周期是天,受到大氣溫度和太陽輻射的共同影響;年溫度作用的周期是年,主要受到大氣溫度的影響;汽車作用的時間較短,且與地區經濟發展水平、人口密度等相關.

橋梁上的不同作用具有多時間尺度的特點,其產生的響應也具有不同的周期. 通常視各類效應在時間尺度上不存在耦合關系,即某一作用產生的撓度不同時受橋梁上其他作用的影響. 因此,可選擇加法模型[13]作為撓度響應增量ΔD和溫度T的表達,即

(1)

式中:DT為溫度效應;DP為汽車、風等其他作用效應;DR為系統的測試誤差;Td為日溫度作用;Ty為年溫度作用;Tf為驟然降溫. 與溫度增量信號相對應,結構的溫度效應DT可以認為是由日溫度效應、年溫度效應和驟然降溫效應組成.

由于撓度響應在時間尺度上具有多尺度特征[14-15],且溫度與其他作用的作用周期差別較大,因此以溫度作用周期對應的頻率為基礎,對撓度信號進行多尺度分析是將溫度效應分離的有效手段.

1.2 多尺度分析理論

多尺度分析(multi-resolution analysis, MRA)是將信號分解為一個概略部分和一系列細節部分[16]. 其中概略部分對應于信號的低頻部分,而細節部分對應于信號的高頻部分. 現實中采集到信號一般是能量有限的,所有能量有限的函數集合所形成的平方可積函數空間是L2(R). 對于信號S,多尺度分析即將其分解為概略部分A1及細節部分D1,針對A1部分再次分解,得到概略部分A2和細節部分D2,如此類推.

(2)

cj,k=〈f(t),φj,k(t)〉.

(3)

式中cj,k稱為尺度展開系數.

(4)

dj,k=〈f(t),ψj,k(t)〉.

(5)

式中dj,k稱為小波展開系數.

將L2(R)按照空間組合展開,有

(6)

其中J為任意設定的尺度.

對于f(t)∈L2(R),可以展開為

(7)

當J→時,式(7)則變為

(8)

對于多尺度分解問題亦可從頻帶劃分的角度來理解. 如果把原始信號X(ω)的頻譜總帶寬0~π定義為空間V0,經過第1級分解,V0被劃分低頻的V1子空間(頻帶為0~π/2)和高頻的W1子空間(頻帶為π/2~π). 經過第2級分解,V1空間被劃分為低頻的V2空間(頻帶為0~π/4)和高頻的W2空間(頻帶為π/4~π/2),如此類推.

小波分析的數學理論基礎是多尺度分析理論. 由多尺度子空間上尺度基所產生的小波基是構造正交小波基的框架[16]. 實際上,式(8)即為離散正交小波變換. 離散正交小波變換的多尺度分析實質是對頻率域的劃分,且按2的整數次冪變化逐次降低分辨率[17-18]. 小波多尺度分析的系數能夠反映信號在時域和頻域上的局部化特征.

1.3 溫度效應的多尺度分離方法

下面將采用小波多尺度分析方法對斜拉橋主梁撓度響應中的溫度效應進行分離,分離步驟包括小波函數的選取、分解層數的選擇以及溫度效應的重構等[19],本文對實測信號的分離步驟如圖1所示.

選取合適的小波函數對于信號成分的準確分離具有重要意義,小波函數選取應對小波函數的正交性、對稱性、消失矩和支撐長度進行考量,表1中歸納了常見的小波函數及其特征.

由表1可知Symlet(symN)小波系列具有較好的正交性、對稱性、足夠的消失矩和支撐長度,能夠精細除噪、快速變換、減少分析重構時的相位失真并能反映信號的奇異性.

圖1 實測信號的分離步驟

表1 常見小波函數的特征

此外,小波多尺度分析實質上是對頻率域的劃分,分析層數需要根據信號源及目標信號的頻率范圍共同確定. 由于撓度響應中車輛的作用效應具有奇異性,應選擇消失矩階數稍高的小波函數. 綜上,本文選擇分解層數為5層的Symlet(sym5)小波函數進行小波多尺度分析.

在溫度效應的多尺度分離中,確定小波函數以及小波分解層數之后,再使用快速傅里葉變換對實測溫度和撓度信號進行頻譜分析,找出其特征頻率. 利用已定參數的小波多尺度分析方法對原始撓度信號進行分解,最后,利用特征頻率所在層的小波系數對撓度信號進行重構,實現溫度效應的分離.

2 溫度效應分離示例

2.1 數據來源

本文選取南京長江第三大橋為背景橋進行研究. 此橋為半漂浮結構體系的雙塔雙索面鋼塔鋼箱梁斜拉橋,跨徑布置為:63 m+257 m+648 m+257 m+63 m=1 288 m. 鋼箱梁采用正交異性板流線型扁平鋼箱梁,梁高3.2 m,寬(含風嘴)37.2 m;斜拉索采用高強度平行鋼絲外擠包高密度聚乙烯形式,全橋共8×21對,上下游對稱分布;橋墩采用混凝土;基礎采用鉆孔灌注樁;索塔為鋼箱形式. 該橋健康監測系統于2006年投入使用,包括溫濕度、風速、應力、索力、撓度、整體位移和結構動力特性監測子系統.

該橋主梁上的溫度傳感器測點位于南塔岸測0號塊主梁,用于監測主梁處大氣溫度. 溫度傳感器型號為Young41372,量程為-50~50 ℃,測量精度為0.3 ℃. 撓度監測子系統由沿梁體布設的封閉式連通管傳感器組成,用于監測各點相對于基點的相對變形量. 撓度傳感器型號為RoseMount 3051S,精度為0.025%. 全橋共72個主梁撓度測點,位于各斜拉索與主梁錨固處,采樣頻率為10 Hz. 上游側撓度測點如圖2所示. 作為示例,圖中標注了南塔江側的撓度測點名稱分別為SJ3、SJ5、SJ7、SJ9、SJ11、SJ13、SJ15、SJ17、SJ19、HL;溫度測點名稱為TEM.

2.2 原始數據預處理

長期監測系統受到不可抗力因素影響時,如電力系統中斷、噪聲、通信故障,會導致采集數據的質量下降[20]. 為得到有效的數據,應對原始數據信號進行預處理.

圖2 傳感器的布設

長期監測系統的數據預處理中,對于原始數據單點缺失現象,可利用樣條插值法填補缺失值,該方法利用當前數據信息,保持數據曲線的性質. 對于原始數據的非連續野值,應利用統計學方法進行粗差剔除,常用方法為“3σ法則”. 一般的常數測量值符合正態分布,其中σ代表標準差,μ代表均值. 該法則認為測量值的取值基本集中在(μ-3σ,μ+3σ)區間. 對于噪聲含量較高的數據信號,應進行去噪處理. 小波分析理論在時間和頻率上具有很好的局部性,對信號的去噪效果良好.

因此,對于長期監測系統中的原始數據,本文依次采用三次樣條差值法對缺失的數據點進行補齊、采用3σ法則進行野值剔除、采用小波濾波方法進行去噪. 本文以2007年8月24日跨中撓度數據為例說明典型數據預處理過程,并將預處理結果繪于圖3.

(a)跨中撓度的原始信號

(c)跨中撓度的野值剔除、數據解算結果

(b)跨中撓度缺失數據處理結果

(d)經消噪后的跨中撓度

圖3 原始數據預處理效果

Fig.3 Preprocessing effect of initial data

由圖3(a)、3(b)可知,原始數據中密集的遺漏值掩蓋了真實撓度變化區域,通過樣條插值法處理缺失數據,數據質量提高. 值得說明的是,原始信號向撓度數據的轉換遵循一定的規則,轉換公式為yi=[(vi-vi0)-(vk-vk0)]×A,其中yi為測點i處的撓度,vi為測點i測量值,vi0為測點i零值,vk為基準點測量值,vk0為基準點零值,A為解算系數,這里A=1 500.613. 將圖3(b)中的原始信號經野值剔除、數據解算后得到實際撓度數據,如圖3(c)所示. 可知實際撓度數據合理,可大致看出汽車荷載和溫度對撓度造成的影響. 由圖3(d)可知,利用小波濾波的方法可以有效去除高頻信號成分,使撓度信號更為合理.

2.3 撓度響應的多尺度分離

如2.2節所示,本文所用數據均按照上述步驟經過預處理. 以主跨跨中測點為例,選取1 a的實測信號對大氣溫度和撓度信號進行分析. 首先快速傅里葉變換,并以1 h的信號平均值作為該時段的特征值,得到其頻譜分別如圖4、5所示.

圖4 大氣溫度實測信號頻譜

由圖4、5可知,撓度和溫度信號的幅值主要分布在0.05 Hz以下的低頻部分(此處視小時為秒進行頻率單位刻畫). 在1.14×10-4Hz(周期365×24 h=8 760 h)附近,撓度信號主要受年溫度等緩慢變化的作用影響. 由于該頻率處溫度和撓度信號均有明顯的集中幅值分布且干擾較小,可認為該頻率是年溫度作用頻率. 在4.17×10-2Hz(周期24 h)附近,溫度信號和撓度信號均有明顯的集中幅值分布,可認為該頻率是日溫度作用頻率.

圖5 跨中撓度響應實測信號頻譜

選取2007年該橋跨中處的實測撓度數據,使用小波多尺度分析方法分別在日溫度作用和年溫度作用頻率上重構原始撓度信號,得到日溫度作用產生的撓度響應(稱為日溫度效應),和年溫度作用產生的撓度響應(稱為年溫度效應),如圖6(a)、6(b)所示. 選取2007年8月份連續幾日內分離出的日溫度效應和年溫度效應進行對比,如圖6(c)所示.

(a)日溫度效應

(b)年溫度效應

(c)溫度效應對比

Fig.6 Separation results of thermal deflection at different time scales

由圖6(a)可知,分離的日溫度效應較平穩,相鄰時刻的幅值較為接近,且夏季的波動幅值明顯大于冬季. 由圖6(b)可知,分離出的年溫度效應的變化較小,且與原始信號的總體變化趨勢一致. 由圖6(c)可知,日溫度效應保留了日溫度作用下撓度響應的細節特征,頻率明顯高于年溫度效應. 日溫度效應以0 mm為零點,以24 h為周期上下波動,波動幅值受太陽輻射強度影響;年溫度效應在短期內變化較小,反應撓度效應的長期變化趨勢. 這符合日溫度作用和年溫度作用的周期特性和季節特性,證實了該分離方法的有效性. 此外,年溫度效應的變化是[27.04 mm,66.28 mm],日溫度效應的變化是[-88.67 mm,99.70 mm]. 與年溫度相比,日溫度效應具有幅值較大,季節性較強的特點. 考慮到日溫度作用中包含太陽輻射作用,因此在日溫度效應的分析中要綜合考慮.

3 日溫度效應分析與預測

3.1 日溫度效應特點

日溫度效應由日溫度作用引起,同時受大氣溫度和太陽輻射的影響. 太陽輻射能夠在結構內部產生溫度梯度,與長期緩慢變化的年溫度作用相比,這種溫度梯度的變化更為迅速,對結構受力的影響更大. 大地表面接收的太陽輻射強度在一年中隨季節不斷改變,在一天中隨時間不斷改變[21]. 因此為了更準確地分析溫度效應,應該對太陽輻射影響的特點進行分析.

3.1.1 年周期中太陽輻射影響的特點

為研究不同季節間太陽輻射對日溫度效應影響的差異,選取2007年冬季和夏季各55 d的日溫度效應數據進行分析,如圖7所示. 日溫度效應變化范圍為冬季[-52.72 mm,44.56 mm],夏季[-88.67 mm,98.93 mm]. 日溫度效應的變化幅值在同一季節中相對平穩,而在不同季節間差別較大.

圖7 冬季和夏季日溫度效應對比

Fig.7 Comparison of diurnal thermal deflection in winter and summer

3.1.2 日周期中太陽輻射影響的特點

文獻[22]中規定了不同緯度處、不同空氣透明度等級條件下,夏季的太陽總輻射強度. 背景橋處于北緯32°地區,大氣透明度等級為5級,夏季水平面單位面積上的太陽輻射強度日變化曲線如圖7所示. 夏季和冬季的太陽輻射日變化趨勢相同,而太陽輻射強度不同[23-24]. 為研究日周期內太陽輻射對日溫度效應的影響,選取2007年冬季、夏季中天氣晴朗的兩天(1月25日和8月25日)日溫度效應同繪于圖8.

圖8 太陽輻射強度與典型日溫度效應

Fig.8 Solar radiation intensity and typical diurnal thermal deflection

由圖8可知,日周期內太陽輻射強度的變化曲線類似余弦曲線,并于中午12:00時達到最大. 日周期中撓度響應在0:00時至8:00時間沒有明顯的波動,8:00時開始日溫度效應明顯增長,主梁出現上拱. 冬季主梁上拱至15:00時達到最大,夏季則上拱至17:00時達到最大,隨后主梁上拱逐漸減小. 撓度的日溫度效應變化與太陽輻射強度的變化規律相似,但相對太陽輻射強度的變化存在一定的滯后.

取2007年冬季、夏季天氣晴朗的各55 d的數據為樣本,經統計可知,撓度日溫度效應相對于太陽輻射強度的滯后性表現為冬季滯后3 h,夏季滯后5 h. 若不考慮滯后影響,將上述典型日的日溫度效應與太陽輻射強度變化曲線繪于圖9.

圖9 不考慮滯后日溫度效應與太陽輻射的相關性

Fig.9 Correlation between thermal deflection and solar radiation without lag

由圖9可知,撓度的日溫度效應與太陽輻射強度呈明顯的正相關關系,其中冬季的相關系數均值為R2=0.795 8, 夏季的相關系數均值為R2=0.880 4.

太陽輻射對日溫度效應的影響在年周期內具有季節性差異,在日周期內具有時刻間差異,因此在對日溫度效應的分析中應充分考慮.

3.2 日溫度效應的預測

預測未來某周期內橋梁結構響應的信息,對于橋梁運營狀態的預警有積極意義. 為實現日溫度效應的預測,應先定量研究日周期內大氣溫度和太陽輻射對日溫度效應的影響. 圖10為典型日2007年1月25日、8月25日的大氣溫度和日溫度效應的關系散點圖(橫坐標為大氣溫度,縱坐標為日溫度效應,折線中各點對應各時刻的平均大氣溫度與日溫度效應).

根據圖10可將日溫度效應分為升溫段:冬季9:00—15:00,夏季9:00—17:00;降溫段:冬季16:00—22:00,夏季18:00—22:00;平穩段:22:00—次日8:00. 平穩段基本處于夜間,太陽輻射的影響很小,大氣溫度和溫度效應的變化都較小;升溫段和降溫段分別由兩條曲線構成,主梁撓度變化較大. 由圖10中折線可知,在日周期內各時刻的日溫度效應與大氣溫度表現為線性相關,日溫度效應可由多折線模型描述.

(a)1月25日

(b)8月25日

Fig.10 Correlation between diurnal thermal deflection and ambient temperature

為考察上述3個區段內大氣溫度與日溫度效應之間的定量關系,分別對2007年冬季、夏季天氣晴朗的55 d中各個時刻的大氣溫度與日溫度效應進行線性擬合. 其單位溫度效應(即線性擬合系數)的平均絕對值及線性擬合決定系數見表2.

由表2可知各時刻間大氣溫度與日溫度效應的線性相關性良好,可使用線性擬合方法對各時刻間的日溫度效應進行擬合[25].

日溫度效應實質上是日周期內各時刻溫度效應的線性累加,滿足多折線模型特征,可表示為

(9)

其中:Dk為1 d中k時刻的日溫度效應,mm;D0為1 d中0:00時的日溫度效應,mm; ΔDi為i時刻的單位溫度效應絕對值, mm/℃;ΔTi為i時刻大氣溫度的變化,℃.

隨機抽取2008年冬季和夏季中天氣晴朗的幾天,同樣采用時間多尺度分析方法將日溫度效應從總撓度響應中分離出來,再利用上述多折線模型對其日溫度效應在日周期內進行預測,預測結果如圖11所示.

表2 各時刻間大氣溫度與日溫度效應的相關性

Tab.2 Correlation between thermal deflection and ambient temperature at different times

時刻ΔDi/(mm·℃-1)冬季夏季R2冬季夏季1:0012.21316.0060.9270.9282:0011.39312.7380.8520.8713:009.6259.1280.7940.7584:005.2000.6490.8840.8945:002.2665.7180.8160.9216:0010.65719.9710.8890.8737:0016.23824.8110.9100.8778:0016.71921.0900.8840.7859:0016.23320.0870.8190.86910:0017.84923.7800.9040.97411:0016.09922.6160.9450.94412:0016.10821.0470.9630.93713:0015.43118.5890.9670.96414:0010.47216.9640.9530.94615:006.40210.6620.9190.92616:003.3409.5320.8870.92217:008.5714.5890.8760.85818:0011.53823.6660.9090.94619:0022.05228.1620.9240.92820:0023.92329.0370.9280.93421:0023.95228.2110.9710.94622:0019.75325.1850.9450.95123:0014.97622.1440.8660.91324:0013.73818.9900.8980.888

圖11 日溫度效應多折線模型預測效果

Fig.11 Prediction of diurnal thermal deflection by multi-fold line model

由圖11可知平穩段內日溫度效應的預測絕對誤差在2 cm以內,升溫段和降溫段內日溫度效應的預測絕對誤差在1 cm以內,日周期內預測平均相對誤差為8.86%,決定系數R2=0.973 4. 預測日溫度效應能反應實際日溫度效應的變化趨勢,且預測值與實際值能較好吻合,因此應用該方法在日周期內對日溫度效應進行預測具有一定的準確性.

4 結 論

1)為更好地分析溫度對斜拉橋主梁撓度的影響,結合撓度響應的多時間尺度特征,對其溫度效應進行多尺度分離和分析. 采用小波多尺度分析方法對撓度響應信號進行不同時間尺度內的分解,能得到不同時間尺度上的溫度效應,實現撓度的溫度效應分離,可以對橋梁狀態評估的數據選取提供參考.

2)對背景橋撓度實測數據進行多尺度分離及分析,結果表明日溫度效應的季節性明顯,不考慮滯后時與太陽輻射強度呈明顯相關關系. 在定量研究橋梁日溫度效應時,需充分考慮太陽輻射引起的季節性差異和時刻間差異.

3)日周期內各時刻間的日溫度效應與大氣溫度間具有良好的線性相關性,本文對其關系提出了多折線擬合模型,綜合考慮大氣溫度和太陽輻射對日溫度效應的影響,獲得不同季節日周期內任意時刻的日溫度效應,實現對溫度效應的可靠預測.

4)所提出的多折線預測模型可為斜拉橋主梁撓度的溫度響應預測提供新思路,對基于監測數據的病害診斷和安全評估等相關研究有積極意義.

猜你喜歡
大氣效應信號
大氣的呵護
軍事文摘(2023年10期)2023-06-09 09:15:06
鈾對大型溞的急性毒性效應
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
懶馬效應
今日農業(2020年19期)2020-12-14 14:16:52
完形填空二則
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
應變效應及其應用
大氣古樸揮灑自如
大氣、水之后,土十條來了
新農業(2016年18期)2016-08-16 03:28:27
基于LabVIEW的力加載信號采集與PID控制
主站蜘蛛池模板: 亚洲国产精品无码久久一线| 亚洲日韩国产精品无码专区| 亚洲丝袜第一页| 国产性爱网站| 91精品人妻互换| 麻豆AV网站免费进入| 欧美一级高清片欧美国产欧美| 伊人91在线| 无码精品国产dvd在线观看9久| 农村乱人伦一区二区| 国产91透明丝袜美腿在线| 日本高清有码人妻| 国产永久在线观看| 国产人免费人成免费视频| 亚洲中文字幕在线观看| 经典三级久久| 日本午夜影院| 国产精品爆乳99久久| 麻豆精品在线| 日本人真淫视频一区二区三区| 国产人在线成免费视频| 欧美亚洲日韩不卡在线在线观看| 五月婷婷导航| 欧美精品亚洲精品日韩专区| 小说 亚洲 无码 精品| a级毛片在线免费| 亚洲精品视频免费看| 最新无码专区超级碰碰碰| 色婷婷色丁香| 亚洲AV色香蕉一区二区| 国产亚洲精品97AA片在线播放| 亚洲二三区| 成人午夜福利视频| 亚洲性网站| 日韩a在线观看免费观看| AV天堂资源福利在线观看| a级毛片毛片免费观看久潮| AV天堂资源福利在线观看| 一区二区日韩国产精久久| 欧美日韩另类在线| 国产网友愉拍精品视频| 色欲不卡无码一区二区| 亚洲欧洲日本在线| 久久五月天国产自| 久久精品人人做人人综合试看| 欧美精品在线观看视频| 亚洲人成网址| 国产激情第一页| 国产日本一区二区三区| 日本不卡在线视频| 亚洲av无码久久无遮挡| 国产无码制服丝袜| 国产一级视频久久| 精品伊人久久久大香线蕉欧美| 一级毛片免费播放视频| 久久窝窝国产精品午夜看片| 2018日日摸夜夜添狠狠躁| 欧美国产在线看| 97精品国产高清久久久久蜜芽| 超清人妻系列无码专区| 日韩专区欧美| 青青久视频| 国产精品爽爽va在线无码观看| 日韩精品中文字幕一区三区| 视频二区亚洲精品| 在线免费看黄的网站| 欧美有码在线| 婷婷六月在线| 中美日韩在线网免费毛片视频| 2022国产91精品久久久久久| 手机成人午夜在线视频| 波多野结衣视频一区二区| 日本精品一在线观看视频| 国产精品成人一区二区不卡| 免费又爽又刺激高潮网址 | 国产精品熟女亚洲AV麻豆| 国禁国产you女视频网站| 亚洲精品免费网站| 亚洲bt欧美bt精品| 无码在线激情片| 中文字幕首页系列人妻| 国产精品原创不卡在线|