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

考慮發射電流波形的地–空瞬變電磁三分量響應研究

2021-11-03 07:03:38侯彥威郭建磊司銀女
煤田地質與勘探 2021年5期
關鍵詞:影響研究

侯彥威,郭建磊,司銀女,姜 濤,寧 輝

考慮發射電流波形的地–空瞬變電磁三分量響應研究

侯彥威,郭建磊,司銀女,姜 濤,寧 輝

(中煤科工集團西安研究院有限公司,陜西 西安 710077)

地–空瞬變電磁法在煤炭采空區勘探等領域受到了越來越多的關注,有必要研究發射電流波形參數對地–空瞬變電磁三分量響應特征影響,為地–空瞬變電磁數據處理與解釋提供理論依據。以梯形波為例,首先,研究不同發射電流波形的頻譜分布情況;然后,基于三維時域有限差分正演研究發射波形的上升沿時間、脈寬和關斷時間對地–空瞬變電磁三分量磁場響應的影響。結果表明:上升沿時間對三分量二次場響應基本不產生影響;關斷時間對三分量二次場響應的影響主要集中在0.2 ms之前,且關斷時間越長對純異常響應影響越大;脈寬對三分量二次場響應的影響主要集中在0.1 ms之后,且脈寬越短對純異常響應影響越大。三維采空區模型結果表明:關斷時間、脈寬對三分量異常場和背景場響應的影響特征基本一致;通過三分量純異常場響應多測道圖和時間道圖可以判斷異常體的分布范圍和深度。研究成果可為地–空瞬變電磁激勵源波形的參數選取提供有價值的理論借鑒。

地–空瞬變電磁;發射電流波形;三維時域有限差分;三分量

地-空瞬變電磁法[1]將接收探頭搭載于無人機上進行測量,該方法有效結合了航空瞬變電磁法[2]和地面瞬變電磁法[3]的優勢,在環境地質調查[4]、煤田采空區勘查[5]等領域具有較大優勢。該方法最早起源于1988年,M. N. Nabighian[6]基于水平電偶源提出該方法;P. Elliot[7]研制了FLAIRTEM系統,該系統鋪設發射線圈邊長較大,在地形特別復雜的地區施工較為困難;R. S. Smith等[8]對航–空、地–空、地面瞬變電磁法進行了對比研究,證明地–空瞬變電磁法具有勘探深度大、信噪比高的優點;T. Mogi等[9]研制了在地表使用電性源作激勵源、使用直升機攜帶磁通門磁力儀采集磁場信號方式的GREATEM系統;嵇艷鞠等[10]研制了用無人機飛艇作為載體搭載電磁接收系統、采集磁場時間導數的無人飛艇長導線源時域地–空電磁勘探系統,解決了我國飛行管制的問題;李肅義等[11]研究了小波去噪方法,有效壓制了噪聲,提高了電磁數據的電阻率成像質量;李貅團隊[12-14]從發射系統、解釋等方面提出了多源發射系統、全域視電阻率定義、逆合成孔徑成像、擬地震成像及地-井與地-空聯合解釋等技術,進一步提升了地-空瞬變電磁法的理論水平;李賀和孫懷鳳[15-16]從時間域出發實現了三維有限元和三維時域有限差分正演,滿足了復雜模型計算需求;覃慶炎[17]、趙越等[18]分析了發射波形對航空瞬變電磁和淺海瞬變電磁全波形響應特征的影響。

上述均沒有涉及發射電流波形對地–空瞬變電磁三分量響應特征的影響。理想的瞬變電磁激勵源為周期性雙極方波電流,受到儀器裝置的限制發射電流上升和關斷都需要一定的時間[19],關斷時間也會造成探測盲區[20]。基于此,研究發射電流波形對地–空瞬變電磁三分量響應特征影響,為地-空瞬變電磁數據處理與解釋提供理論依據。

1 正演理論基礎

在均勻、各向同性、有耗、非磁性、無源媒質中,Maxwell方程組為:

式中:為電場強度;為磁感應強度;為磁場強度;為介質電導率;表示介電常數;為時間。

地球物理勘探一般忽略位移電流,為構成顯式時間迭代格式,在式(1b)中加入虛擬介電常數項,將激勵源電流密度加入Maxwell方程組的安培環路定理實現激勵源的加載,在有源區域式(1b)變為:

地-空瞬變電磁法一般采用電性源進行一次場激發[1],為便于正演計算假定電性源加載在直角坐標系的軸方向(圖1)。

圖1 電性源與相鄰晶胞網格的位置

必須考慮低頻條件下正演結果的準確性,對磁場各分量采用低頻近似進行處理,將式(1a)、式(2)、式(1d)在直角坐標系下展開,得到各分項形式:

對式(3)和式(4)用差分代替微分,由于Euler前向差分對離散時間步的要求比較嚴格,故空間離散采用后向差分,均勻網格剖分中電(磁)場的時間采樣恰好在兩相鄰磁(電)場采樣時刻的中心,故時間離散采用中心差分,可以得到6個分量的迭代公式。其中,無源區域的電磁場迭代方程見參考文獻[21],有源區域E分量迭代格式為:

式中:,,分別為,,方向的網格數。

采用第一類邊界條件進行模型剖分(即在邊界處將電場和磁場強制性賦零),需要將整個計算模型剖分得盡可能大。計算過程中需要滿足時間域和空間域的穩定性條件,如下式。

在激勵源中加入階躍電流進行電流源加載,電流源考慮上升沿、持續時間和下降沿。本次采用梯形波作為激勵源,激勵源電流發射示意圖如圖2所示。

激勵源函數為:

式中:1為上升沿時間;2為上升沿時間加脈寬;3為上升沿時間加脈寬加關斷時間。通過修改1、2、3可有效模擬不同上升沿時間、脈寬及關斷時間等參數對瞬變電磁三分量響應的影響。

建立如圖3所示的地-空瞬變電磁均勻半空間模型。參數如下:地層電導率0.01 S/m,激勵源長度500 m,電流方向沿軸正向、大小1 A,采用梯形波(上升沿0.1 μs、持續時間100 ms、關斷時間0.1 μs)進行發射,接收測點位于(300 m,-50 m)處,飛行高度50 m,二次場采樣時間10 ms,剖分均勻網格數為221×221×200,網格尺寸為10 m。將三維計算結果與一維解析解進行對比,得到響應對比圖和誤差對比圖(圖4)。由圖4發現,三分量響應一維解析解與數值解的誤差均在6%以下,滿足計算精度要求。

圖3 測點與激勵源相對位置俯視圖

2 考慮發射電流波形的響應特征分析

2.1 電流波形分辨能力

一定寬度的脈沖信號的頻率響應占據一定的頻帶寬且二者之間存在定量對偶關系,對激勵源函數進行傅里葉變換可以得到激勵源函數的頻譜方程。分析發射波形的上升沿時間、脈寬和關斷時間對頻率分布影響,相關時間參數見表1。對表1所述的不同激勵源進行傅里葉變換得到頻率分布圖(圖5)。由圖5a、圖5c發現,縮短上升沿時間和關斷時間后其振幅幅值稍微增大,但差別很小;改變上升沿時間和關斷時間基本不影響激勵源的頻率分布情況。由圖5b發現,改變脈寬時間明顯改變過零點帶寬和初始振幅,脈寬越小其初始振幅越小、過零點帶寬越大、高頻成分越多,高頻成分越多其探測分辨率越高。經計算發現過零點帶寬與脈寬時間呈反相關關系(過零點帶寬等于脈寬時間倒數)。

圖4 均勻半空間模型FDTD解與解析解計算結果對比曲線

表1 激勵源時間參數

2.2 改變梯形波發射波形參數的響應特征

基于均勻半空間模型(圖3)研究激勵源的不同參數對地–空瞬變電磁三分量全波形響應和二次場響應特征的影響,二次場響應起始時刻以發射電流完全關斷時刻開始計算。

研究上升沿時間變化對三分量響應的影響,設置如下發射源參數:上升沿時間分別為1、20、50、100 μs,脈寬10 ms,關斷時間1 μs。三維計算后得到全波形響應(圖6a)和二次場響應(圖6b)。

圖6a顯示不同上升沿時間的分量響應幅值為負,關斷瞬間響應值趨于零,全波形響應整體與發射波形特征相似;分量響應的一次場響應幅值為正,關斷后響應幅值由正變負后續逐漸趨于零;分量與分量響應特征一致,但響應幅值極值大于分量。不同上升沿時間的三分量全波形響應主要在上升沿時間和脈寬的早期存在差別,上升沿時間越短其一次場響應幅值極值越早達到。對于二次場響應,分量響應早期數值為負值,分量響應數值為負值,為方便成圖對計算的三分量二次場響應數據取絕對值后進行雙對數繪圖得到圖6b。圖6b所示不同上升沿時間的三分量二次場響應曲線完全重合,說明上升沿時間對二次場響應基本不產生影響。

為研究關斷時間變化對三分量響應的影響,設置如下發射源參數:關斷時間分別為1、20、50、100 μs,上升沿時間1 μs,脈寬10 ms。三維計算后得到全波形響應(圖7a)和二次場響應(圖7b)。

圖7a顯示不同關斷時間的全波形響應特征與發射波形相似,關斷后響應曲線出現細微分離,但不明顯。圖7b顯示不同關斷時間對三分量二次響應特征產生較大影響,關斷時間越長其三分量二次場早期響應幅值越小,0.2 ms以后響應曲線基本重合,說明關斷時間對三分量二次場響應的影響主要集中在0.2 ms之前,隨著時間的延遲其影響逐漸減弱。關斷時間為1 μs的分量響應“跳變”時間道最晚,對比不同關斷時間的分量響應的“跳變”時間道所示關斷時間越短分量響應的“跳變點”越靠晚期時間道。

圖5 不同參數激勵源的頻率分布

為研究脈寬變化對三分量響應的影響,設置如下發射源參數:脈寬分別為1、5、10、15、20 ms,上升沿時間1 μs,關斷時間1 μs。三維計算后得到全波形響應(圖8a)和二次場響應(圖8b)。

圖8a所示全波形響應特征與發射波形相似,脈寬越長其一次場響應越長,不同脈寬全波形響應特征一致。圖8b發現脈寬變化對三分量二次場響應產生較大影響,脈寬越短,三分量二次場晚期響應幅值越小,0.1 ms以前早期響應曲線基本重合,說明關斷時間對三分量二次場響應的影響主要集中在0.1 ms之后;通過分量的“跳變”時間道可見脈寬為1 ms的分量響應“跳變”比其他脈寬的分量響應“跳變”靠晚期時間道;脈寬達到10 ms及以上時三分量二次場響應曲線幾乎重合,說明脈寬達到一定時間后其對二次場響應的影響較為微弱。

圖8 脈寬變化對三分量響應的影響

2.3 三維采空區充水模型響應特征

為研究采空區勘探的激勵源設計及參數選取問題,建立圖9所示的地–空瞬變電磁三維充水采空區模型,參數如下:激勵源長度500 m,電流方向沿軸正向,電流大小1 A,接收點位于(300 m, –50 m)處,飛行高度50 m,地層電導率為0.01 S/m,充水采空區電導率1 S/m,頂界面距地面100 m,大小為150 m×150 m×100 m。地-空瞬變電磁主要依據二次場數據進行解釋,影響二次場響應的主要因素為關斷時間和脈寬,因此,接下來主要分析關斷時間和脈寬對充水采空區模型的二次場響應的影響。將帶異常體的響應定義為異常場響應,不帶異常體的響應定義為背景場響應,異常場響應減去背景場響應得到純異常響應。

圖9 三維采空區模型俯視圖

為研究關斷時間的影響,設置如下發射源參數:關斷時間分別為1、20、100 μs,上升沿時間1 μs,脈寬10 ms。三維計算后得到三分量二次場響應(圖10a)和三分量純異常響應(圖10b)。

圖10 關斷時間變化對三分量響應的影響

圖10a所示關斷時間對三分量異常場響應和背景場響應的影響特征基本一致,影響主要集中在早期,關斷時間越長其早期響應幅值越大,且分量響應的“跳變”靠晚期時間道。圖10b所示關斷時間越長對三分量純異常響應影響越大,不同關斷時間響應數據比對表現為關斷時間越長純異常響應開始時間和達到純異常極值的時間越早。

為研究脈寬的影響,設置如下發射源參數:脈寬時間分別為1、10、20 ms,上升沿時間1 μs,關斷時間1 μs。三維計算后得到三分量二次場響應(圖11a)和三分量純異常響應(圖11b)。

由圖11a發現,脈寬對三分量異常場響應和背景場響應的影響特征基本一致,影響主要集中在晚期,脈寬越短其晚期響應幅值越大,且分量響應“跳變”越靠晚期時間道。由圖11b發現,脈寬為1 ms對三分量純異常場響應影響較大,表現為三分量純異常場響應整體幅值偏高;脈寬為10 ms和20 ms時三分量純異常場響應曲線基本重合,說明脈寬達到一定時間后其對二次場響應的影響較為微弱甚至不產生影響。

以脈寬10 ms、上升沿時間1 μs、關斷時間1 μs的激勵源進行三維計算得到三維采空區模型的純異常響應多測道圖(圖12)和時間道圖(圖13),通過多測道圖可知和分量呈現下凹形狀,分量呈“S”型,多測道圖和時間道圖呈現良好的對應關系,分量和分量呈現負值的圈閉,分量呈現左側負值、右側正值,并且異常主要反映在50~60道之間。上述模型說明通過純異常多測道圖和時間道圖判斷異常體的分布范圍和深度。

3 結論

a. 基于傅里葉變換研究不同發射電流波形的頻譜分布情況,發現上升沿時間和關斷時間基本不影響激勵源的頻率分布情況,脈寬與零點帶寬呈反相關關系,脈寬越短其高頻成分越多。

b. 基于三維時域有限差分正演研究上升沿時間、脈寬和關斷時間對三分量響應的影響,發現上升沿時間對三分量響應基本不產生影響;關斷時間的影響主要集中在0.2 ms之前;脈寬的影響主要集中在0.1 ms之后。

圖11 脈寬變化對三分量響應的影響

圖12 三分量純異常響應多測道圖

圖13 三分量純異常響應時間道圖

c. 通過三維采空區模型發現:關斷時間、脈寬對三分量異常場和背景場響應的影響特征基本一致;其中關斷時間越長其對三分量純異常場影響越大;脈寬越短對二次場影響越大,當脈寬達到一定時間后其影響較為微弱甚至不產生影響;通過三分量純異常場響應多測道圖和時間道圖可以判斷異常體的分布范圍和深度。

d. 不同地層電導率、接收高度、偏移距及復雜地質模型等條件下不同發射電流波形參數對地–空瞬變電磁三分量響應特征的影響,有待進一步研究。

[1] 張瑩瑩,李貅. 地空瞬變電磁法研究進展[J]. 地球物理學進展,2017,32(4):1735–1741.

ZHANG Yingying,LI Xiu. Research progress on ground-airborne transient electromagnetic method[J]. Progress in Geophysics,2017,32(4):1735–1741.

[2] 許洋. 半航空電磁一維正反演研究[D]. 成都:成都理工大學,2014.

XU Yang. Study about 1D forward and inversion of SAEM[D]. Chengdu:Chengdu University of Technology,2014.

[3] 侯彥威. TEM探測深部煤層上覆多電性層的OCCAM反演[J]. 煤田地質與勘探,2018,46(6):169–173.

HOU Yanwei. OCCAM inversion for detecting overlying multiple electrical layers above deep coal seams by TEM[J]. Coal Geology & Exploration,2018,46(6):169–173.

[4] 吳啟龍. 半航空瞬變電磁視電阻率成像及在復雜地形區域隧道勘察中的應用[D]. 濟南:山東大學,2019.

WU Qilong. Semi-airborne transient electromagnetic apparent resistivity imaging and its application in tunnel survey in complex terrain areas[D]. Jinan:Shandong University,2019.

[5] 張慶輝,田忠斌,林君,等. 時域電性源地空電磁系統在煤炭采空積水區勘查中的應用[J]. 煤炭學報,2019,44(8):2509–2515.

ZHANG Qinghui,TIAN Zhongbin,LIN Jun,et al. Application of time domain electrical source ground airborne electromagnetic system in goaf water exploration[J]. Journal of China Coal Society,2019,44(8):2509–2515.

[6] NABIGHIAN M N. Electromagnetic methods in applied geophysics:Volume 1,theory[M]. Society of Exploration Geophysicists,1988.

[7] ELLIOTT P. The principles and practice of FLAIRTEM[J]. Exploration Geophysics,1998,29(1/2):58–59.

[8] SMITH R S,ANNAN P,MC GOWAN P D. A Comparison of data from airborne,semi-airborne,and ground electromagnetic systems[J]. Geophysics,2001,66(5):1379–1385.

[9] MOGI T,TANAKA Y,KUSUNOKI K,et al. Development of grounded electrical source airborne transient EM (GREATEM)[J]. Exploration Geophysics,1998,29(1/2):61–64.

[10] 嵇艷鞠,王遠,徐江,等. 無人飛艇長導線源時域地空電磁勘探系統及其應用[J]. 地球物理學報,2013,56(11):3640–3650.

JI Yanju,WANG Yuan,XU Jiang,et al. Development and application of the grounded long wire source airborne electromagnetic exploration system based on an unmanned airship[J]. Chinese Journal of Geophysics,2013,56(11):3640–3650.

[11] 李肅義,林君,陽貴紅,等. 電性源時域地空電磁數據小波去噪方法研究[J]. 地球物理學報,2013,56(9):3145–3152.

LI Suyi,LIN Jun,YANG Guihong,et al. Ground-airborne electromagnetic signals de-noising using a combined wavelet transform algorithm[J]. Chinese Journal of Geophysics,2013,56(9):3145–3152.

[12] 張瑩瑩,李貅,姚偉華,等. 多輻射場源地空瞬變電磁法多分量全域視電阻率定義[J]. 地球物理學報,2015,58(8):2745–2758.

ZHANG Yingying,LI Xiu,YAO Weihua,et al. Multi-component full field apparent resistivity definition of multi-source ground-airborne transient electromagnetic method with galvanic sources[J]. Chinese Journal of Geophysics,2015,58(8):2745–2758.

[13] 李貅,張瑩瑩,盧緒山,等. 電性源瞬變電磁地空逆合成孔徑成像[J]. 地球物理學報,2015,58(1):277–288.

LI Xiu,ZHANG Yingying,LU Xushan,et al. Inverse synthetic aperture imaging of ground-airborne transient electromagnetic method with a galvanic source[J]. Chinese Journal of Geophysics,2015,58(1):277–288.

[14] 武軍杰. 地–井與地–空瞬變電磁聯合解釋方法研究[D]. 西安:長安大學,2018.

WU Junjie. Study on joint interpretation of borehole and ground-airborne TEM data[D]. Xi’an:Chang’an University,2018.

[15] 孫懷鳳,李貅,李術才,等. 考慮關斷時間的回線源激發TEM三維時域有限差分正演[J]. 地球物理學報,2013,56(3):1049–1064.

SUN Huaifeng,LI Xiu,LI Shucai,et al. Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time[J]. Chinese Journal of Geophysics,2013,56(3):1049–1064.

[16] 李賀. 直接時間域矢量有限元瞬變電磁三維正演模擬[D]. 西安:長安大學,2016.

LI He. Three-dimensional transient electromagnetic forward modeling in the direct time-domain by vector finite element[D]. Xi’an:Chang’an University,2016.

[17] 覃慶炎. 發射電流波形變化對固定翼航空瞬變電磁響應的影響[J]. 地球物理學進展,2013,28(5):2673–2679.

QIN Qingyan. Effect of transmitter current waveform variation on the fixed-wing airborne transient electromagnetic response[J]. Progress in Geophysics,2013,28(5):2673–2679.

[18] 趙越,許楓,李貅,等. 淺海瞬變電磁全波形響應特征及探測能力分析[J]. 地球物理學報,2019,62(4):1526–1540.

ZHAO Yue,XU Feng,LI Xiu,et al. Exploration capability of transmitter current waveform on shallow water TEM response[J]. Chinese Journal of Geophysics,2019,62(4):1526–1540.

[19] 姜升. 瞬變電磁發射機電流波形改善技術研究[D]. 重慶:重慶大學,2017.

JIANG Sheng. Research on improvement for quality of current waveform of time-domain electromagnetic transmitter[D]. Chongqing:Chongqing University,2017.

[20] 張軍. 礦井超淺層高分辨率瞬變電磁探測技術[J]. 煤田地質與勘探,2020,48(4):219–225.

ZHANG Jun. The high-resolution transient electromagnetic detection technology for ultra-shallow layer in coal mine[J]. Coal Geology & Exploration,2020,48(4):219–225.

[21] 郭建磊. 陣列源瞬變電磁法隧道高分辨超前探測研究[D]. 西安:長安大學,2017.

GUO Jianlei. Transient electromagnetic method using source array for tunneling ahead geological prospecting of high resolution[D]. Xi’an:Chang’an University,2017.

Research on three-component response of ground-airborne TEM considering emission current waveform

HOU Yanwei, GUO Jianlei, SI Yinnyu, JIANG Tao, NING Hui

(Xi’an Research Institute Co. Ltd., China Coal Technology and Engineering Group Corp., Xi’an 710077, China)

The ground-airborne TEM has received more and more attention in coal goaf exploration and other fields. In order to provide a theoretical basis for the processing and interpretation of ground-airborne transient electromagnetic data, it is necessary to study the influence of the transmitted current waveform parameter on the three-component response characteristics of ground-airborne TEM. Taking trapezoidal waves as an example, this article firstly studies the frequency distribution of different emission current waveform, and then investigates the effects of the rising edge time, pulse width and turn-off time of the transmitted waveform on the three-component magnetic field response of ground-air transient electromagnetic based on the three-dimensional finite difference time domain method(3D-FDTD). The results show that the rising edge time basically has no effect on the three-component secondary field response; the turn off time effect on the three-component secondary field response is mainly concentrated before 0.2 ms, and the longer the turn off time, the greater the impact on the pure abnormal response; the pulse width has the main effect on the three-component secondary field response after 0.1 ms, and the shorter the pulse width, the greater the impact on the pure abnormal response. The results of the three-dimensional goaf model show that the characteristics of the off-time and pulse width on the response of the three-component abnormal field and background field are basically the same; the distribution range and depth of anomalous objects can be judged by the three-component pure anomaly field response multi-track map and time track map. The research results will provide some valuable theoretical references for the selection of the parameters of the ground-airborne TEM excitation source waveform.

ground-airborne TEM; emission current waveform; 3D-FDTD; three-component

移動閱讀

語音講解

P631

A

1001-1986(2021)05-0238-09

2021-01-06;

2021-06-22

陜西省自然科學基礎研究計劃項目青年基金項目(2020JQ-994);中煤科工集團西安研究院有限公司科技創新基金項目(2019XA YMS30)

侯彥威,1983年生,男,河南商丘人,碩士,副研究員,研究方向為煤田電磁法勘探. E-mail:houyanwei@cctegxian.com

侯彥威,郭建磊,司銀女,等.考慮發射電流波形的地–空瞬變電磁三分量響應研究[J]. 煤田地質與勘探,2021,49(5):238–246. doi: 10. 3969/j. issn. 1001-1986. 2021. 05. 026

HOU Yanwei,GUO Jianlei,SI Yinnyu,et al. Research on three-component response of ground-airborne TEM considering emission current waveform[J]. Coal Geology & Exploration,2021,49(5):238–246. doi: 10. 3969/j. issn. 1001-1986. 2021. 05. 026

(責任編輯 聶愛蘭)

猜你喜歡
影響研究
FMS與YBT相關性的實證研究
是什么影響了滑動摩擦力的大小
2020年國內翻譯研究述評
遼代千人邑研究述論
哪些顧慮影響擔當?
當代陜西(2021年2期)2021-03-29 07:41:24
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
新版C-NCAP側面碰撞假人損傷研究
沒錯,痛經有時也會影響懷孕
媽媽寶寶(2017年3期)2017-02-21 01:22:28
擴鏈劑聯用對PETG擴鏈反應與流變性能的影響
中國塑料(2016年3期)2016-06-15 20:30:00
主站蜘蛛池模板: 99爱在线| 亚洲va在线∨a天堂va欧美va| 国产尤物jk自慰制服喷水| 欧美在线视频不卡第一页| 狠狠色狠狠色综合久久第一次| 91原创视频在线| 亚洲AV无码乱码在线观看裸奔| 奇米精品一区二区三区在线观看| 国产成人精品视频一区视频二区| 不卡无码h在线观看| 亚洲中字无码AV电影在线观看| 秋霞一区二区三区| 国内精品久久九九国产精品| 成人免费一区二区三区| 亚洲视频a| 九九久久精品免费观看| 午夜啪啪网| 中文字幕调教一区二区视频| 国产精品欧美激情| 粗大猛烈进出高潮视频无码| 亚洲欧美精品日韩欧美| 91久久偷偷做嫩草影院电| 97精品伊人久久大香线蕉| 国产成人精品一区二区三区| 免费无码网站| 亚洲日韩在线满18点击进入| 美女一级毛片无遮挡内谢| a级毛片免费播放| 亚洲中文字幕97久久精品少妇| 特级精品毛片免费观看| 天天色天天操综合网| 国产玖玖视频| 国产xx在线观看| 永久免费AⅤ无码网站在线观看| 成人国产一区二区三区| 日韩国产无码一区| 国产成人高清精品免费软件| 在线观看国产精美视频| 国产无码网站在线观看| 亚洲精品你懂的| 久久国产V一级毛多内射| 亚洲高清资源| 精品国产电影久久九九| 国产精品刺激对白在线| 久久人搡人人玩人妻精品| 高清欧美性猛交XXXX黑人猛交 | 激情爆乳一区二区| 97人人做人人爽香蕉精品| 熟女日韩精品2区| 四虎成人免费毛片| 香蕉久人久人青草青草| 欧美日韩中文字幕在线| a毛片在线| 欧美翘臀一区二区三区| 国产高清免费午夜在线视频| 国产手机在线小视频免费观看| 国产swag在线观看| 日韩欧美在线观看| 萌白酱国产一区二区| 免费无码网站| 91精品国产丝袜| 99久久精品国产麻豆婷婷| 免费一级毛片不卡在线播放| 在线观看视频99| 亚洲欧美国产视频| 真实国产精品vr专区| 国产在线观看第二页| 国产丝袜啪啪| 亚洲成人黄色在线| 黄色一级视频欧美| 国产亚洲精品97AA片在线播放| 久久亚洲国产一区二区| 青青草原国产免费av观看| 日韩在线2020专区| 四虎在线观看视频高清无码| 国产麻豆精品久久一二三| 国产麻豆va精品视频| 毛片在线播放a| 九色视频在线免费观看| 亚洲精品动漫在线观看| 国产91小视频在线观看| 色噜噜中文网|