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

青島地區鉆孔傾斜儀地鐵干擾定量分析

2025-03-01 00:00:00李瀟岳龍徐清風臧藝博韓幫杰王哲李煒
地震工程學報 2025年2期
關鍵詞:模態信號方法

摘要: 受青島地鐵運行影響,青島地區鉆孔傾斜儀存在高頻干擾,每日干擾出現時間與地鐵運行時間一致,通過定性分析青島地區五套鉆孔傾斜儀數據,發現距離地鐵2 km內的傾斜儀會受到地鐵干擾,當距離超過15 km時,不再受地鐵干擾。為定量分析地鐵的干擾特征,利用變分模態分解(VMD)方法對青島地區受地鐵干擾的傾斜儀數據進行分解,同時結合廣義S變換,以驗證VMD分解層數K值的準確性及去噪的有效性,分析得出青島地區受不同地鐵線路干擾的鉆孔傾斜儀均具有5個干擾模態分量。去除地鐵干擾分量并重構數據后,數據信噪比提升,有效提高了數據可用性。

關鍵詞: 鉆孔傾斜儀; 地鐵; VMD; 去噪

中圖分類號: P315"" """文獻標志碼:A"" 文章編號: 1000-0844(2025)02-0383-10

DOI:10.20000/j.1000-0844.20231106001

Quantitative analysis of subway interference on

borehole tiltmeter in the Qingdao area

LI Xiao1, YUE Long1, XU Qingfeng1, ZANG Yibo1, HAN Bangjie1, WANG Zhe1, LI Wei2

(1. Qingdao Earthquake Monitoring Center, Shandong Earthquake Agency, Qingdao 266071, Shandong, China;

2. Qingdao Emergency Management Affairs Service Center, Qingdao 266034, Shandong, China)

Abstract:

Due to the subway's operations, high-frequency interference occurs on the tiltmeters in the Qingdao area. Furthermore, the daily interference occurs simultaneously with the subway's operation time. After conducting a qualitative analysis of data obtained from five sets of borehole tiltmeters in the Qingdao area, the results reveal that the tiltmeters within 2 km of the subway are disturbed by the subway's operations. However, when the distance exceeds 15 km, the tiltmeters are no longer affected by subway interference. At the same time, to quantitatively analyze the interference characteristics of the subway, the variational mode decomposition (VMD) method was used in this study to decompose the data of tiltmeters affected by the Qingdao subway interference. The accuracy of the K value of the VMD decomposition layer and the effectiveness of denoising were also verified using the generalized S transform. The findings indicate that the borehole tiltmeters undergoing interference from different subway lines in Qingdao have five interference mode components. Upon removing these components and reconstructing the data, the signal-to-noise ratio of the data is improved, effectively improving data availability.

Keywords:

borehole tiltmeter; subway; VMD; denoising

0 引言

地球物理場前兆異常判別對地震預測預報有重要意義。一些學者利用不同方法研究了不同地震的震前異常變化,張維辰等1對姑咱臺的鉆孔應變數據進行時頻分析,提取到蘆山MS7.0地震前高頻畸變信號。唐磊等2利用分量應變數據進行主張應變方向動態分析,初步得到門源6.9級地震前異常變化。隨著經濟社會發展,地球物理場觀測受到的干擾越來越多,如何排除環境干擾,提取真實有效的前兆異常具有重大社會需求和深遠的科學意義3-4。因此,需要研究更有效的噪音壓制方法,池成全5系統地研究了固體潮、氣溫、氣壓、水位對應變觀測的影響并加以去除。

隨著地鐵等現代交通工具的快速發展,地鐵運行導致流入大地的雜散電流等干擾因素對于地電阻率、地電場、地磁等觀測手段的影響也越來越普遍。趙俊香等6以上海崇明地震臺地電場觀測數據為研究對象,利用 Welch 功率譜估計方法,對 ZD9A-Ⅱ地電場分采樣數據進行分析,得出地鐵干擾的頻率范圍,并選擇適合的小波基進行數據分解與重構,較好地還原了原始信號。樊曉春等7研究了地鐵對地電阻率觀測的影響;謝凡等8研究了天津地區城市軌道直流牽引交通系統對地磁觀測的干擾;王同利等9研究了北京城市軌道交通對地電場觀測干擾的影響,目前針對地鐵對地球物理觀測數據干擾的研究主要集中在對電磁儀器的干擾研究方面,而針對地鐵運行對鉆孔傾斜儀的影響還未開展系統研究。

青島地震監測中心站(以下稱青島臺)安裝有體應變儀、四分量鉆孔應變儀、豎直擺鉆孔傾斜儀和氣象三要素儀等四套地球物理觀測儀器,臺站觀測點距離青島地鐵二號線燕兒島地鐵站0.7 km左右。根據青島臺日常數據分析和異常核實結果,確定青島地鐵2號線對青島臺豎直擺鉆孔傾斜儀造成影響,導致固體潮曲線加粗。而同場地觀測的其他兩套儀器未受到地鐵干擾。

針對非平穩態隨機數據處理,常用分析方法包括小波變換法10、基于經驗模態分解(Empirical Mode Decomposition,EMD)的希爾伯特黃變換法等11,而小波分析的基函數缺乏自適應性,傳統的EMD方法又存在模態混疊現象和端點效應。因此,Dragomiretskiy等12提出了變分模態分解方法(Variational Mode Decomposition,VMD)方法,相較于傳統的小波分解、EMD分解等信號分解方法,VMD方法是一種自適應、完全非遞歸的方法,其可以根據信號的實際情況來確定分解個數,同時適合復雜度高且非線性比較強的非平穩時間序列信號。該方法能夠避免計算迭代過程中的端點效應和虛假分量問題。VMD方法大量應用于醫療、電力、自動控制等不同領域的一維信號處理,大部分相關領域的研究均表明該方法比小波分解、EMD分解有著更好的處理效果13

本文以青島地區為例,利用VMD方法研究青島地區鉆孔傾斜儀受地鐵干擾情況與干擾特征,并去除地鐵干擾。

1 觀測儀器概況

青島地區共有五套CZB-1c型豎直擺鉆孔傾斜儀,分別位于青島臺、黃島臺、膠南臺、即墨臺、萊西臺。其中青島臺儀器安裝深度57 m,黃島臺儀器安裝深度30 m,膠南臺儀器安裝深度83 m,即墨臺儀器安裝深度67 m,萊西臺儀器安裝深度88 m[14。五個臺站與地鐵的位置關系如圖1所示,臺站與地鐵線路的距離以及周圍地鐵線路的運行時間如表1所列。

2 干擾分析

選取青島臺鉆孔傾斜儀2022年5月18—21日東西分量和北南分量的分鐘值數據[圖2(a)、圖2(c)]和一階差分數據[圖2(b)、圖2(d)]進行分析可知,東西分量脈動加粗,且相對于原始數據,這種干擾影響在差分數據上更加明顯。將鉆孔傾斜儀每日脈動加粗的時間與地鐵2號線運行時間對比,二者基本一致。

為進一步確定地鐵運行的影響范圍,本文選取了其他四套鉆孔傾斜儀2022年1月24—30日的數據,其中黃島臺距離地鐵1號線1 km,膠南臺距離地鐵13號線約2 km,即墨臺距離最近的地鐵11號線約15 km,萊西臺距離最近的地鐵11號線約40 km。四套傾斜儀的原始數據和一階差分曲線如圖3~6所示。

如圖3和圖4所示,距離地鐵2 km左右的膠南臺和距離地鐵1 km左右的黃島臺鉆孔傾斜儀在地鐵運行時間內均受明顯的地鐵干擾。

如圖5所示,即墨臺鉆孔傾斜儀原始曲線存在高頻干擾,一階差分數據每日00:00—24:00幅度一致。為判斷該干擾是否為地鐵干擾,選取2014年9月10—15日青島地區無地鐵運行時即墨臺傾斜儀數據做對比(圖6)。如圖6所示,2014年即墨臺原始數據的高頻干擾,每日00:00—24:00一階差分數據幅值沒有明顯變化,且幅值與2022年基本一致,證明地鐵開通運行對即墨臺鉆孔傾斜儀無影響。如圖7所示,萊西臺的北南分量和東西分量一階差分數據每日00:00—24:00幅值無明顯變化,表明地鐵運行對萊西臺鉆孔傾斜儀無影響。由此綜合推斷當鉆孔傾斜儀觀測井距離地鐵15 km以上,不受地鐵運行的影響。另外三個臺站則受地鐵活動影響,表2是受干擾臺站鉆孔傾斜儀的干擾幅度統計表。由表2可知,隨著臺站與地鐵的距離增加,地鐵干擾幅度呈下降趨勢。

3 VMD方法原理

本文將采用變分模態分解方法將青島、黃島、膠南臺鉆孔傾斜儀數據分解成各模態分量,并對青島臺分離出來的模態分量做頻譜分析和重構信號去噪。

VMD是一種信號分解的估計方法,實現過程包括:(1)迭代搜索變分模型最優解,將原始信號分解為若干個模態分量,而VMD假定分離出來的各模態信號均為集中在各自中心頻率附近的窄帶信號;(2)根據分量窄帶條件建立約束優化問題,進而計算出模態分量的中心頻率;(3)根據每一個模態分量各自的中心頻率和帶寬,實現各分量信號的有效剝離。

傳統的經驗模態分解缺乏嚴密的數學基礎做支撐,即在數學上沒有嚴格的證明;而變分模態分解是一種變分問題,數學基礎嚴密,適用于處理地震數據等非平穩態數據。對于有噪音的信號應具有良好的穩定性;而對于將要研究的原始信號,假定它由多個有限帶寬的模態分量(或稱IMF分量) νk(t)構成,對其分解得出的任意一個模態分量,它的中心頻率為ω(t)。問題求解的約束條件則是所有求解出來的模態分量之和等于原始信號。VMD算法是將輸入的非穩定非線性信號分解成多個本征模態IMF分量,再將信號分解構造成變分約束問題,其過程能夠分成變分問題的約束和變分問題的求解兩個過程。以輸入的原始信號等于分解分量的和,以及每個分解的分量預估帶寬最小為目標,構造變分問題的約束,公式為:

min{uk},{ωk}∑Nk=1t(δ(t)+jπt)·uk(t)e-jωkt22

s.t. ∑Nk=1uk=f(t)(1)

式中:f(t)為輸入的原始信號;N為分解層數;uk\,ωk分別為第k個分解得到的模態分量以及其中心頻率;δ(t)為單位脈沖函數;t為t時刻的偏導數。為了求解式(1)的最優解,引入二次懲罰函數因子和拉格朗日乘數,能夠得到以下公式:

L[(uk),(ωk),λ]=α∑Nk=1tδ(t)+jπt·uk(t)e-jωkt22+

f(t)-∑Nk=1uk(t)22+λ(t),f(t)-∑Nk=1uk(t)(2)

式中:λ為拉格朗日乘法算子;α為懲罰因子。計算過程采用乘法算子交替方向法迭代計算式(2)的最優解,直到達到停止迭代的條件。停止迭代的條件為:

∑Nk=1un+1k-unk22/unk22lt;ε (3)

式中:ε為收斂容差。達到式(3)的條件即可停止迭代,輸出N個窄帶本征模態函數。

雖然VMD方法具有諸多優勢,但也存在一定的局限性,例如VMD方法并沒有完全解決模態分解存在的端點效應,而且在計算過程中,需要提前設定兩個重要參數:分解層數K和懲罰因子。一旦參數設置不當,甚至具有物理意義的復雜原始信號,也會影響分解效果。本文在設置參數過程中,需要通過大量試驗,分析信號分解效果,同時結合廣義S變換時頻譜特征,確定最優參數。

4 信號分解與去噪

4.1 信號分解

選取2022年5月18—19日青島臺鉆孔傾斜儀東西分量數據、2022年8月8—9日黃島臺鉆孔傾斜儀北南分量數據、2022年1月30—31日膠南臺鉆孔傾斜儀北南分量數據為研究對象,三個臺站原始波形和一階差分結果,如圖8所示。

如圖8所示,3個臺站鉆孔傾斜儀每日存在高頻干擾。利用VMD方法對青島臺鉆孔傾斜儀東西分量數據進行模態分解,經過參數測試,選擇分解層數K為7,懲罰因子為2 000時,分解出7個模態分量。然后對7個模態分量進行傅里葉變換,得到各分量的頻譜圖,結果如圖9所示。圖中u1~u7分量分別是分解得到的7個模態分量,u1分量為固體潮部分;u2分量為00:00一直持續到24:00,證明u2分量與地鐵運行無關;u3~u7分量為06:00之后高頻干擾開始出現,并一直持續到23:23。圖9(b)是各個分量的頻譜圖,從頻譜圖中可以看出各分量均有優勢頻率和一定的頻帶范圍,且頻譜無混疊,5個高頻干擾優勢頻率分別是0.007 8 Hz、0.006 3 Hz、0.004 4 Hz、0.003 7 Hz、0.002 1 Hz。

同樣,利用VMD方法對黃島臺和膠南臺鉆孔傾斜儀北南分量數據進行模態分解,經過參數測試,選擇分解層數K為7,懲罰因子為2 000,分解出7個模態分量,然后對這7個模態分量進行傅里葉變換,得到各分量的頻譜圖,結果如圖10、11所示。由圖10可知,黃島臺鉆孔傾斜儀數據中地鐵的5個高頻干擾優勢頻率分別是0.007 8 Hz、0.006 3 Hz、0.005 0 Hz、0.003 8 Hz、0.002 0 Hz。膠南臺鉆孔傾斜儀數據中地鐵的5個高頻干擾優勢頻率分別是0.007 8Hz、0.006 4 Hz、0.005 8 Hz、0.004 3 Hz、0.003 8 Hz(圖11)。通過對比圖9~11可知,u1分量為原始信號分解得到的固體潮部分,其余6個分量為細節項。6個細節分量中,u2在地鐵開始運行時刻前后沒有明顯變化,其余5個分量在地鐵運行前沒有干擾;而地鐵運行后,出現明顯干擾,干擾時間段和地鐵運行時間高度相關,因此地鐵運行導致的模態分量主要集中在u3~u7。由分解結果可知,干擾青島臺、黃島臺和膠南臺的地鐵線路不同,但地鐵干擾的模態分量均為5個,分解結果具有一定相似性。

4.2 去噪分析

由圖8、9可知,2022年5月18日06:00地鐵干擾開始出現,干擾信號的模態分量06:00之后開始出現并一直持續到23:23。如前所述,u1、u2分量為非地鐵干擾部分,u3~u7分量為地鐵干擾部分。因此在信號重構時,將u3~u7分量部分舍棄,其他部分重新相加,即可得到去除干擾的數據。

為了分析VMD去噪效果,對去噪前后的數據分別做廣義S變換時頻分析,得到時頻如圖12(b)和圖12(d)所示。圖12(a)是原始數據,圖12(c)是利用VMD方法去噪后的數據。從圖12(b)中可以看出,06:00地鐵開始運行后,時頻譜圖顯示地鐵高頻干擾的優勢頻率段與VMD分解得到的干擾模態分量的數量、優勢頻率范圍、頻率幅值一致,驗證了VMD方法分解層數K設置準確,表明該方法可以將地鐵干擾的5個模態分量分離出來并去噪,提高了信號的信噪比。

5 結論與討論

青島臺鉆孔傾斜儀受距離0.7 km的青島地鐵2號線的影響,在地鐵運行時段內,東西分量脈動加粗,北南分量原始數據曲線干擾不明顯,一階差分處理后,也存在明顯的脈動加粗現象。通過綜合分析地鐵對青島地區五套鉆孔傾斜儀的干擾影響,得到以下結論:

(1) 通過分析青島臺、黃島臺、膠南臺、即墨臺、萊西臺的鉆孔傾斜儀數據的原始數據和一階差分數據,結合不同臺站離地鐵線的距離和地鐵運行時間等要素,可知位于地鐵2 km范圍內的鉆孔傾斜儀會受地鐵運行的干擾;在2~15 km之間,現有觀測條件無法確定地鐵的運行是否會對鉆孔傾斜儀造成干擾;但可以確定的是,當距離地鐵線超過15 km時,鉆孔傾斜儀基本不受地鐵運行的干擾。

(2) 本文利用VMD方法對受地鐵干擾的青島臺、黃島臺和膠南臺鉆孔傾斜儀數據進行分解,當分解層數設置為7,懲罰因子為2 000時,能夠將地鐵干擾的5個模態分量與有效信號模態分量分離。將干擾模態分量壓制并重構,得到去除地鐵干擾的有效數據。為了驗證分解和去噪效果,對于原始數據進行廣義S變換得到時頻譜,分析得出地鐵高頻干擾的優勢頻率段與VMD方法分解得到的地鐵干擾模態分量的數量、頻率范圍和幅度均一致,也證明VMD方法參數設置的準確性和去噪的有效性,該去噪方法可以為地震異常識別提供更加真實可信的基礎數據。

參考文獻(References)

[1] 張維辰,朱凱光,池成全,等.基于小波變換的2013年蘆山MS7.0地震前姑咱臺鉆孔應變異常時頻分析[J].地震學報,2019,41(2):230-238.

ZHANG Weichen,ZHU Kaiguang,CHI Chengquan,et al.Time-frequency analyses for borehole strain anomaly at Guzan station before 2013 Lushan MS7.0 earthquake based on wavelet transform[J].Acta Seismologica Sinica,2019,41(2):230-238.

[2] 唐磊,邱澤華,李玉江,等.四分量鉆孔應變觀測近地表應力應變狀態的判定與分析[J/OL].武漢大學學報(信息科學版),1-15[2023-10-12].https://doi.org/10.13203/j.whugis20220397.

TANG Lei,QIU Zehua,LI Yujiang,et al.Determination and analysis of near-surface stress-strain state of 4-component borehole strainmeter[J/OL].Geomatics and Information Science of Wuhan University(Information Science Edition),1-5[2023-10-12].https://doi.org/10.13203/j.whugis20220397.

[3] 蔡佩蕊,陳偉,林立峰,等.基于STFT方法對沿海寬頻帶傾斜儀的噪聲分析[J].地震工程學報,2020,42(2):396-402.

CAI Peirui,CHEN Wei,LIN Lifeng,et al.Noise analysis of coastal broadband tiltmeter based on STFT method[J].China Earthquake Engineering Journal,2020,42(2):396-402.

[4] 于紫凝.鉆孔應變觀測數據的震前異常提取與評價方法研究[D].長春:吉林大學,2021.

YU Zining.Study on anomaly extraction and evaluation method of borehole strain observation data before earthquake[D].Changchun:Jilin University,2021.

[5] 池成全.鉆孔應變前兆觀測數據分析與異常提取研究[D].長春:吉林大學,2020.

CHI Chengquan.Study on data analysis and anomaly extraction of borehole strain precursor observation[D].Changchun:Jilin University,2020.

[6] 趙俊香,朱培育,畢波,等.崇明地電場地鐵干擾分析[J].地震地磁觀測與研究,2016,37(2):124-129.

ZHAO Junxiang,ZHU Peiyu,BI Bo,et al.The analysis of subway disturbance for Chongming geoelectric field[J].Seismological and Geomagnetic Observation and Research,2016,37(2):124-129.

[7] 樊曉春,李偉,袁慎杰,等.地鐵對井下地電阻率觀測的影響分析:以江寧臺為例[J].地震工程學報,2020,42(3):680-687.

FAN Xiaochun,LI Wei,YUAN Shenjie,et al.Influence of the subway on deep-well resistivity observation:a case study of Jiangning station[J].China Earthquake Engineering Journal,2020,42(3):680-687.

[8] 謝凡,滕云田,徐學恭,等.直流牽引城市軌道交通系統對地磁觀測干擾研究:以天津軌道交通為例[J].地球物理學進展,2011,26(2):732-738.

XIE Fan,TENG Yuntian,XU Xuegong,et al.Magnetic perturbation caused by D. C. traction power supply system for urban rail transit:a case study from Tianjin railway transit in China[J].Progress in Geophysics,2011,26(2):732-738.

[9] 王同利,胡樂銀,崔博聞,等.北京城市軌道交通對地電場觀測的干擾影響[J].地震地質,2013,35(4):887-893.

WANG Tongli,HU Leyin,CUI Bowen,et al.A preliminary analysis on the interference generated by urban railway transit in Beijing to the geoelectric field observations[J].Seismology and Geology,2013,35(4):887-893.

[10] 楊志鵬,陳秀清,余洋洋,等.基于小波分解與同步擠壓變換的西昌小廟臺形變典型干擾的時頻響應特征研究[J].地震工程學報,2022,44(5):1192-1206.

YANG Zhipeng,CHEN Xiuqing,YU Yangyang,et al.Time-frequency response characteristics for the typical deformation interference in Xiaomiao station based on wavelet decomposition and the synchrosqueezing transform[J].China Earthquake Engineering Journal,2022,44(5):1192-1206.

[11] 高涵,袁希平,甘淑,等.利用HHT-EEMD方法分析云南區域GNSS應變時序孕震信息[J].測繪學報,2022,51(9):1899-1910.

GAO Han,YUAN Xiping,GAN Shu,et al.Analysis of seismogenic information of GNSS strain time series based on HHT-EEMD method in Yunnan region[J].Acta Geodaetica et Cartographica Sinica,2022,51(9):1899-1910.

[12] DRAGOMIRETSKIY K,ZOSSO D.Variational mode decomposition[J].IEEE Transactions on Signal Processing,2014,62(3):531-544.

[13] 喬云,李瓊,錢浩東,等.基于VMD與改進小波閾值的地震信號去噪方法研究[J].物探化探計算技術,2021,43(6):690-696.

QIAO Yun,LI Qiong,QIAN Haodong,et al.Seismic signal denoising method based on VMD and improved wavelet threshold[J].Computing Techniques for Geophysical and Geochemical Exploration,2021,43(6):690-696.

[14] 岳龍,徐清風,劉云,等.短周期氣壓波對青島臺體應變的影響分析[J].大地測量與地球動力學,2019,39(9):977-981.

YUE Long,XU Qingfeng,LIU Yun,et al.Analysis on the influence of short-period atmospheric pressure wave on volumetric strain of Qingdao station[J].Journal of Geodesy and Geodynamics,2019,39(9):977-981.

(本文編輯:任 棟)

猜你喜歡
模態信號方法
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
基于LabVIEW的力加載信號采集與PID控制
國內多模態教學研究回顧與展望
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 婷婷色丁香综合激情| 欧美日韩国产精品综合| 波多野一区| 最新国产精品鲁鲁免费视频| 黄色污网站在线观看| 久久久黄色片| 欧美综合成人| 亚洲天堂网在线播放| 99精品视频九九精品| 永久成人无码激情视频免费| 日韩欧美一区在线观看| 国产成人综合网| www精品久久| 99re在线免费视频| 国产视频a| 91在线播放免费不卡无毒| 欧美日在线观看| 伊人久久大香线蕉综合影视| 久久国产精品麻豆系列| 中文字幕在线不卡视频| 成年人国产网站| 一级毛片基地| 99免费视频观看| 国产区网址| 国产精品成人一区二区不卡| 亚洲一区二区三区麻豆| 久久男人视频| 波多野结衣一区二区三区AV| 18禁黄无遮挡网站| 亚洲三级片在线看| 色男人的天堂久久综合| 国产午夜小视频| 国产精品99久久久| av在线无码浏览| 无码精油按摩潮喷在线播放 | 亚洲中文字幕国产av| 国产黄色片在线看| 在线免费看片a| 丁香婷婷综合激情| 日韩福利视频导航| av无码久久精品| 免费一级毛片不卡在线播放 | 9999在线视频| 71pao成人国产永久免费视频| 2018日日摸夜夜添狠狠躁| 色视频国产| аⅴ资源中文在线天堂| 伊人中文网| 国产制服丝袜无码视频| 亚洲成人播放| 在线观看亚洲精品福利片| 四虎成人在线视频| 无码免费视频| 国产一级小视频| 久久婷婷色综合老司机| 亚洲精品制服丝袜二区| 四虎精品黑人视频| 免费看久久精品99| 成人毛片免费在线观看| 亚洲第一极品精品无码| 国产男人天堂| 91成人在线免费视频| 看国产毛片| 中文纯内无码H| 国产精品欧美激情| 成人a免费α片在线视频网站| 欧美色99| 久久精品人人做人人综合试看| 日本国产在线| 狠狠色香婷婷久久亚洲精品| a毛片在线| 日本国产在线| 国产91av在线| 国产小视频网站| 久久a级片| 九九九九热精品视频| www.91中文字幕| 色综合中文字幕| 国产精品手机视频一区二区| 午夜成人在线视频| 99青青青精品视频在线| 欧美色图第一页|