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

層狀大地表面中心回線瞬變電磁響應(yīng)特征

2015-05-03 07:45:26李永博
物探化探計(jì)算技術(shù) 2015年5期

吳 瓊, 李永博, 李 貅, 晉 達(dá)

(1.中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所,廊坊 065000;2.長安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院,西安 710054;3.中石化中原油田分公司 物探研究院,濮陽 457001)

?

層狀大地表面中心回線瞬變電磁響應(yīng)特征

吳 瓊1, 李永博1, 李 貅2, 晉 達(dá)3

(1.中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所,廊坊 065000;2.長安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院,西安 710054;3.中石化中原油田分公司 物探研究院,濮陽 457001)

中心回線裝置下水平地層的時(shí)間域電磁場表達(dá)式為雙重積分,內(nèi)層為漢克爾型積分,外層為正弦或余弦積分,線性數(shù)字濾波法是求解此類特殊積分的有效方法。這里采用線性數(shù)字濾波法分別計(jì)算了均勻半空間、兩層和三層介質(zhì)共七種地電斷面的電磁響應(yīng),分析了不同類型響應(yīng)曲線的形態(tài)特征,研究了中間層厚度變化和底層電阻率變化對(duì)時(shí)間域響應(yīng)曲線形態(tài)的影響。通過對(duì)比與分析,得出了一些規(guī)律性的認(rèn)識(shí),對(duì)于瞬變電磁數(shù)據(jù)的處理解釋具有一定的指導(dǎo)意義。

漢克爾積分變換; 正弦變換; 余弦變換; 線性數(shù)字濾波; 瞬變電磁響應(yīng)

0 引 言

瞬變電磁法具有工作效率高、電性分辨力強(qiáng)、探測(cè)深度大、抗干擾能力強(qiáng)等顯著優(yōu)勢(shì),因而被廣泛應(yīng)用于礦產(chǎn)勘探、工程勘察、環(huán)境災(zāi)害地質(zhì)調(diào)查等多個(gè)領(lǐng)域[1-4]。時(shí)間域瞬變電磁響應(yīng)的求取通常有兩種方案:①直接法,根據(jù)一定的初始條件和邊界條件在時(shí)間域中直接求解;②間接法,即頻-時(shí)轉(zhuǎn)換,根據(jù)頻譜分析理論,首先在頻率域求解給定場源的電磁響應(yīng),再通過傅里葉反變換求得相應(yīng)的時(shí)間域電磁響應(yīng)[5]。因頻率域電磁場表達(dá)式較為簡單,故其應(yīng)用更為普遍。

實(shí)現(xiàn)頻-時(shí)轉(zhuǎn)換的主要方法有離散傅里葉變換法、延遲譜法、逆拉普拉斯變換法、數(shù)字濾波法等[5],離散傅里葉變換法因其所需的頻率和核函數(shù)抽樣次數(shù)多、計(jì)算量大、速度慢而極少被采用,延遲譜法的晚期響應(yīng)不太穩(wěn)定,逆拉普拉斯變換法的計(jì)算量大、計(jì)算速度較慢,數(shù)字濾波法相比前三種方法,在層狀大地情況下具有計(jì)算精度高、速度快且晚期響應(yīng)穩(wěn)定的優(yōu)點(diǎn),是當(dāng)前較常用的一種計(jì)算方法。

在線性數(shù)字濾波算法方面,諸多學(xué)者開展了系列的研究工作,D.Guptasarma[6]提出了一種用于計(jì)算時(shí)間域激電響應(yīng)的線性數(shù)字濾波算法,并給出了21點(diǎn)濾波系數(shù);Aderson[7]提出了滯后褶積快速濾波算法;王華軍[8]通過快速漢克爾變換理論導(dǎo)出了正弦變換和余弦變換數(shù)字濾波算法,給出了兩組250點(diǎn)濾波系數(shù),并驗(yàn)證了這組濾波算法具有計(jì)算速度快、精度高的優(yōu)點(diǎn)[9]。

這里采用250點(diǎn)正弦變換濾波算法,分別計(jì)算了一層(均勻半空間)、兩層和三層水平介質(zhì)的時(shí)間域電磁響應(yīng),總結(jié)了典型地電斷面的電磁響應(yīng)曲線特征,并研究了中間層厚度變化和底層電阻率變化對(duì)時(shí)間域響應(yīng)曲線形態(tài)的影響,得到了一些規(guī)律性的認(rèn)識(shí)。

1 理論基礎(chǔ)

一個(gè)半徑為a的圓形線圈水平放置在水平地層表面,通入電流的大小為I0,角頻率為ω,磁導(dǎo)率為μ。水平地層的電導(dǎo)率由上到下依次為σ1、σ2、…、σn-1、σn,各地層厚度依次為H1、H2、…、Hn-1、∞。采用柱坐標(biāo)系,設(shè)原點(diǎn)在圓回線中心點(diǎn),取z軸向下為正。當(dāng)z=0時(shí),接收線圈位于地表,此時(shí)各電磁場分量在頻率域的積分表達(dá)式為[1]:

(1)

(2)

(3)

當(dāng)忽略位移電流,并在中心點(diǎn)接收時(shí):

Eφ=Hr=0

(4)

(5)

公式(5)為各向同性水平地層表面圓回線源中心點(diǎn)的頻率域電磁場表達(dá)式,在回線中心點(diǎn)只能觀測(cè)到磁場的垂直分量。在均勻半空間這一特殊情況下,回線源中心點(diǎn)頻率域垂直磁場的積分表達(dá)式可化簡為式(6)。

(6)

諧變電流條件下的電磁場量Hz(ω)與階躍電流條件下的電磁場量Hz(t)滿足如下關(guān)系:

(7)

利用歐拉公式、Hz(ω)的復(fù)數(shù)性質(zhì)、δ(ω)函數(shù)的積分性質(zhì)、積分區(qū)間的對(duì)稱性、ReHz(ω)和lmHz(ω)的奇偶性及其頻率特性等[1],可得:

(8)

其相應(yīng)的時(shí)間導(dǎo)數(shù)為:

(9)

文中水平層狀介質(zhì)的時(shí)間域響應(yīng),采用250點(diǎn)正弦變換濾波算法,計(jì)算公式為[9-10]式(10)。

(10)

2 模型計(jì)算

2.1 均勻半空間

為研究時(shí)間域響應(yīng)曲線隨回線半徑的變化規(guī)律,令I(lǐng)0=1A、ρ1=1 Ω·m,回線半徑a分別取五個(gè)不同的值,求得的時(shí)間域響應(yīng)曲線如圖1所示。

圖1 均勻半空間時(shí)間域響應(yīng)(不同回線半徑)

為研究時(shí)間域響應(yīng)曲線隨電阻率的變化規(guī)律,令I(lǐng)0=1 A、a=100 m,均勻半空間電阻率ρ1分別取五個(gè)不同的值,求得的時(shí)間域響應(yīng)曲線如圖2所示。

圖2 均勻半空間時(shí)間域響應(yīng)(不同電阻率)

均勻半空間時(shí)間域響應(yīng)曲線具有如下規(guī)律性特征:在雙對(duì)數(shù)坐標(biāo)下,早期響應(yīng)曲線近乎水平,晚期響應(yīng)曲線為向下傾斜的直線,與時(shí)間軸的夾角約為68.2°。回線源半徑越小,在中心點(diǎn)的早期響應(yīng)幅值越大,越早進(jìn)入晚期,均勻半空間電阻率越小,在中心點(diǎn)的早期響應(yīng)幅值越小,越晚進(jìn)入晚期。

2.2 兩層模型

兩層模型分為G型和D型兩種地電斷面,其參數(shù)取值如圖3所示,線性數(shù)字濾波法求得的時(shí)間域響應(yīng)見圖4。

圖3 兩層模型示意圖

圖4 兩層模型時(shí)間域響應(yīng)

從時(shí)間域響應(yīng)曲線形態(tài)來看,三者在早期階段基本重合,在雙對(duì)數(shù)坐標(biāo)下,為近水平的直線,隨著時(shí)間的推移,G型逐漸向下偏離均勻半空間響應(yīng)曲線,D型先向下后向上偏離均勻半空間響應(yīng)曲線。由于瞬變電磁場服從熱傳導(dǎo)方程的規(guī)律,隨著時(shí)間的增加,場向深處傳播過程中以“煙圈”形式逐漸向外擴(kuò)散,時(shí)間域響應(yīng)逐漸衰減,而衰減曲線的斜率對(duì)應(yīng)著地層的信息,G型的第二層電阻率大于第一層電阻率,故時(shí)間域響應(yīng)衰減速率增大,表現(xiàn)為時(shí)間域響應(yīng)曲線在均勻半空間響應(yīng)曲線的下方,D型的第二層電阻率小于第一層電阻率,故時(shí)間域響應(yīng)衰減速率減小,表現(xiàn)為時(shí)間域響應(yīng)曲線在均勻半空間響應(yīng)曲線的上方。響應(yīng)曲線的兩種不同形態(tài),是由于層狀模型高電阻率目標(biāo)或低電阻率目標(biāo)的電磁場衰減特性引起的。由此可見,瞬變電磁響應(yīng)曲線反映了地下介質(zhì)的電性變化特征,即當(dāng)電磁場由低阻層進(jìn)入高阻層時(shí)衰減加快,當(dāng)電磁場由高阻層進(jìn)入低阻層時(shí)衰減放慢。

2.3 三層模型

三層模型分為A型、K型、H型和Q型四種地電斷面,其參數(shù)取值如圖5所示。線性數(shù)字濾波法求得的時(shí)間域響應(yīng)見圖6。

圖5 三層模型示意圖

圖6 三層模型時(shí)間域響應(yīng)

從時(shí)間域響應(yīng)曲線形態(tài)來看,五者形態(tài)各異。在雙對(duì)數(shù)坐標(biāo)下,隨著時(shí)間的推移,A型逐漸向下偏離均勻半空間響應(yīng)曲線,在中晚期衰減變快;K型和均勻半空間響應(yīng)曲線首尾重合,在中期衰減稍快;H型圍繞均勻半空間響應(yīng)曲線先“下凹”后“上凸”,在晚期逐漸向均勻半空間響應(yīng)曲線靠近;Q型呈“下凹”、“上凸”、“下凹”、“上凸”規(guī)律性變化,尤其在晚期,曲線衰減極為緩慢。總體來看,A型是GG型組合,H型是DG型組合,K型是GD型組合,Q型是DD型組合,因此在中晚期階段,曲線的形態(tài)變化仍是與電磁場是由低阻進(jìn)入高阻還是由高阻進(jìn)入低阻兩種不同電性界面性質(zhì)有關(guān)。

2.4 改變中間層厚度模型

令G型、D型地電斷面的第一層厚度與A型、H型、K型、Q型地電斷面的第二層厚度分別取25m、50m、100m,同類型模型的其他參數(shù)取值保持固定,計(jì)算出的時(shí)間域響應(yīng)曲線如圖7所示。

由圖7(a)、7(b)可見,G型、D型地電斷面的時(shí)間域響應(yīng)曲線在早期與均勻半空間響應(yīng)曲線基本重合,隨著時(shí)間的推移,G型曲線逐漸向下偏離均勻半空間響應(yīng)曲線,D型曲線先向下后向上偏離均勻半空間響應(yīng)曲線,G型和D型的第一層厚度越小,曲線偏離越明顯,在晚期階段,三種不同第一層厚度的G型和D型時(shí)間域響應(yīng)曲線均具有重合趨勢(shì)。由圖7(c)、7(e)可見,A型、K型地電斷面的時(shí)間域響應(yīng)曲線在早期與具有相同地電參數(shù)的G型時(shí)間域響應(yīng)曲線基本重合,隨著時(shí)間的推移,A型曲線逐漸向下偏離G型響應(yīng)曲線,K型曲線逐漸向上偏離G型響應(yīng)曲線,不同第二層厚度的A型和K型時(shí)間域響應(yīng)曲線均比較接近,很難分辨。由圖7(d)、7(f)可見,H型、Q型地電斷面的時(shí)間域響應(yīng)曲線在早期與具有相同地電參數(shù)的D型時(shí)間域響應(yīng)曲線基本重合,隨著時(shí)間的推移,H型曲線逐漸向下偏離D型響應(yīng)曲線,Q型曲線先向下后向上偏離D型響應(yīng)曲線,H型和Q型的第二層厚度越小,曲線偏離越明顯,在晚期階段,三種不同第二層厚度的H型和Q型時(shí)間域響應(yīng)曲線均具有重合趨勢(shì)。

2.5 改變底層電阻率模型

令G型、D型、A型、K型、H型、Q型六種地電斷面的底層電阻率分別取3個(gè)不同的值,同類型模型的其他參數(shù)取值保持固定,計(jì)算出的時(shí)間域響應(yīng)曲線如圖8所示。

由圖8(a)可見,第二層電阻率越大,G型曲線向下偏離均勻半空間響應(yīng)曲線越明顯,但當(dāng)?shù)诙与娮杪试龃蟮揭欢ǔ潭群螅淙≈档牟煌跁r(shí)間域響應(yīng)曲線中難以辨別。由圖8(b)可見,第二層電阻率越小,D型曲線先向下后向上偏離均勻半空間響應(yīng)曲線幅度越大,在晚期階段,第二層電阻率不同的三種D型時(shí)間域響應(yīng)呈直線形態(tài),且平行于均勻半空間情況下的晚期響應(yīng)。由圖8(c)可見,第三層電阻率不同的三種A型時(shí)間域響應(yīng)曲線基本重合,無明顯差異。由圖8(d)可見,第三層電阻率越大,H型曲線向下偏離D型響應(yīng)曲線越明顯,但當(dāng)?shù)谌龑与娮杪试龃蟮揭欢ǔ潭群螅淙≈档牟煌跁r(shí)間域響應(yīng)曲線上反映得不再明顯。由圖8(e)可見,第三層電阻率越小,K型曲線偏離G型響應(yīng)曲線越明顯,在晚期階段,第三層電阻率不同的三種K型時(shí)間域響應(yīng)曲線幾乎平行。由圖8(f)可見,第三層電阻率越小,Q型曲線先向下后向上偏離D型響應(yīng)曲線幅度越大,第三層電阻率不同的三種Q型時(shí)間域響應(yīng)曲線差異明顯,并在晚期階段,呈相互平行的直線。

圖7 改變中間層厚度模型時(shí)間域響應(yīng)

圖8 改變底層電阻率模型時(shí)間域響應(yīng)

3 結(jié) 論

通過對(duì)均勻半空間、兩層和三層介質(zhì)共七種地電斷面的電磁響應(yīng)曲線進(jìn)行對(duì)比分析,得到以下幾點(diǎn)認(rèn)識(shí)和結(jié)論:

1)在雙對(duì)數(shù)坐標(biāo)下,早期的時(shí)間域響應(yīng)曲線為平行于時(shí)間軸的直線;晚期的時(shí)間域響應(yīng)曲線為向下傾斜的直線,它與時(shí)間軸的夾角約為68.2°。

2)回線源半徑越小,在中心點(diǎn)的早期響應(yīng)幅值越大,越早進(jìn)入晚期;地層電阻率越小,在中心點(diǎn)的早期響應(yīng)幅值越小,越晚進(jìn)入晚期。

3)當(dāng)由低阻層進(jìn)入高阻層(G型)時(shí),電磁響應(yīng)衰減加快;當(dāng)由高阻層進(jìn)入低阻層(D型)時(shí),電磁響應(yīng)衰減放慢。

4)A型與G型地電斷面的電磁響應(yīng)非常接近,不易分辨;K型和均勻半空間響應(yīng)曲線首尾重合,在中期衰減稍快,反映了中間層的高阻特征;H型圍繞均勻半空間響應(yīng)曲線先“下凹”后“上凸”,在中期衰減較慢,反映了中間層的低阻特征;Q型呈“下凹”、“上凸”、“下凹”、“上凸”規(guī)律性變化,尤其在晚期,曲線衰減極為緩慢。時(shí)間域電磁法對(duì)高阻層的分辨能力遠(yuǎn)沒有對(duì)低阻層的分辨能力強(qiáng)。

5)G型和D型地電斷面的第一層厚度越小,電磁響應(yīng)曲線偏離均勻半空間響應(yīng)曲線越明顯,在晚期階段,不同第一層厚度的G型和D型時(shí)間域響應(yīng)曲線均具有重合趨勢(shì);不同第二層厚度的A型和K型時(shí)間域響應(yīng)曲線均比較接近,很難分辨;H型和Q型的第二層厚度越小,曲線偏離具有相同地電參數(shù)的D型時(shí)間域響應(yīng)曲線越明顯,在晚期階段,不同第二層厚度的H型和Q型時(shí)間域響應(yīng)曲線均具有重合趨勢(shì)。

6)G型地電斷面的第二層電阻率越大,其響應(yīng)曲線向下偏離均勻半空間響應(yīng)曲線越明顯,但當(dāng)?shù)诙与娮杪试龃蟮揭欢ǔ潭群螅淙≈档牟煌跁r(shí)間域響應(yīng)曲線中難以辨別,在晚期階段,第二層電阻率不同的G型時(shí)間域響應(yīng)曲線具有重合趨勢(shì);D型地電斷面的第二層電阻率越小,其響應(yīng)曲線先向下后向上偏離均勻半空間響應(yīng)曲線幅度越大,在晚期階段,第二層電阻率不同的D型時(shí)間域響應(yīng)呈直線形態(tài),且平行于均勻半空間情況下的晚期響應(yīng)。在中晚期階段,A型和H型時(shí)間域響應(yīng)曲線隨第三層電阻率的變化規(guī)律與G型時(shí)間域響應(yīng)曲線隨第二層電阻率的變化規(guī)律類似,K型和Q型時(shí)間域響應(yīng)曲線隨第三層電阻率的變化規(guī)律與D型時(shí)間域響應(yīng)曲線隨第二層電阻率的變化規(guī)律類似。

上述規(guī)律性認(rèn)識(shí)對(duì)于瞬變電磁資料的處理和解釋,具有一定的參考價(jià)值。

[1] 方文藻,李予國,李貅. 瞬變電磁測(cè)深法原理[M]. 西安: 西北工業(yè)大學(xué)出版社,1993.FANGWZ,LIYG,LIX.Theprincipleoftransientelectromagneticsoundingmethod[M].Xi’an:NorthwesternPolytechnicUniversityPress, 1993. (InChinese)

[2] 樸化榮. 電磁測(cè)深原理[M]. 北京: 地質(zhì)出版社,1990.PIAOHR.Principleofelectromagneticsoundings[M].Beijing:GeologicalPublishingHouse, 1990. (InChinese)

[3] 李貅. 瞬變電磁測(cè)深的理論與應(yīng)用[M]. 西安: 陜西科學(xué)技術(shù)出版社,2002.LIX.Thetheoryandapplicationoftransientelectromagneticsounding[M].Xi’an:ShaanxiScienceandTechnologyPress, 2002. (InChinese)

[4] 牛之璉. 時(shí)間域電磁法原理[M]. 長沙: 中南大學(xué)出版社,2007.NIUZHL.Timedomainelectromagneticmethods[M].Changsha:CentralSouthUniversityPress, 2007. (InChinese)

[5] 陳向斌, 胡社榮, 張超. 瞬變電磁場響應(yīng)計(jì)算的頻-時(shí)域轉(zhuǎn)換方法綜述[J]. 工程地球物理學(xué)報(bào), 2008, 5(2):242-246.CHENXB,HUSR,ZHANGCH.Reviewoffrequency-domainandtime-domaintransformationmethodintransientelectromagneticfieldresponsecomputation[J].ChineseJournalofEngineeringGeophysics, 2008, 5(2): 242-246. (InChinese)

[6]GUPTASARMAD.Computationofthetime-domainresponseofapolarizableground[J].Geophysics, 1982, 47: 1574-1576.

[7]ANDERSONWL.FastHankeltransformsusingrelatedandlaggedconvolutions[J].Assn.Comp.Math.,Trans.,Math.Software, 1982, 8: 344-368.

[8]JOHANSENHK,SORENSENK.FastHankeltransforms[J].GeophysicalProspecting, 1979, 27: 876-901.

[9] 王華軍. 正余弦變換的數(shù)值濾波算法[J]. 工程地球物理學(xué)報(bào), 2004, 1(4): 329-335.WANGHJ.Digitalfilteralgorithmofthesineandcosinetransform[J].ChineseJournalofEngineeringGeophysics, 2004, 1(4): 329-335. (InChinese)

[10]吳瓊. 大回線源電磁場正演與波場變換理論研究[D]. 西安: 長安大學(xué),2012.WUQ.Forwardmodelingoflargelooptransientelectromagneticfieldandwave-fieldtransformtheory’sstudy[D].Xi’an:Chang’anUniversity, 2012. (InChinese)

Response features of central-loop transient electromagnetic field on the layered earth

WU Qiong1, LI Yong-bo1, LI Xiu2, JIN Da3

(1.Institute of Geophysical and Geochemical Exploration, Langfang 065000,China;2. School of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China;3. Geophysical Research Institute, Zhongyuan Oilfield Company, SINOPEC, Puyang 457001, China)

The expressions of time-domain electromagnetic fields with a central-loop on the layered earth are double integrals. The inner integral is Hankel transform while the outer integral is sine transform or cosine transform. Linear digital filtering is an effective method to solve such special integrals. The electromagnetic responses for half space, two-layer and three-layer models have been calculated with linear digital filtering method, and the response curves for the seven types (half space, G, D, A, K, H and Q) have been analyzed, and time-domain responses of layered models with changing interlayer-thickness and underlayer-resistivity have been studied. By comparing and analyzing, some regular cognitions are drawn, which have certain guiding significance for transient electromagnetic data’s processing and interpretation.

Hankel integral transform; sine transform; cosine transform; linear digital filtering; transient electromagnetic response

2015-05-11改回日期:2015-06-08

中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)項(xiàng)目(AS2013P01)

吳瓊(1986-),女,碩士,主要從事電磁勘探理論與方法技術(shù)研究,E-mail:wuqiong8668@163.com。

1001-1749(2015)05-0560-06

P 631.3

A

10.3969/j.issn.1001-1749.2015.05.03

主站蜘蛛池模板: 国产视频你懂得| 人禽伦免费交视频网页播放| 国产JIZzJIzz视频全部免费| 亚洲无码高清一区二区| 国产免费看久久久| 国产欧美精品一区二区| 精品久久久久久中文字幕女 | 四虎成人精品在永久免费| 精品国产一区91在线| 蜜桃臀无码内射一区二区三区| 亚洲成人免费看| 狠狠色丁婷婷综合久久| 91区国产福利在线观看午夜| 免费观看精品视频999| 青青操视频在线| 亚洲天堂日韩在线| 制服丝袜在线视频香蕉| 国产精品尹人在线观看| 国产超碰一区二区三区| 99草精品视频| 日韩精品一区二区深田咏美| 亚洲人妖在线| 国产成人久视频免费| 天天色天天操综合网| 国产玖玖视频| a级毛片一区二区免费视频| 亚洲区一区| 亚洲精品你懂的| 国产精欧美一区二区三区| 免费不卡视频| 日韩毛片免费| 亚洲人成在线精品| 欧美a级在线| 国产鲁鲁视频在线观看| 91亚洲视频下载| 日韩精品成人网页视频在线| 在线国产你懂的| 国产久草视频| 伊大人香蕉久久网欧美| 全部免费毛片免费播放 | 直接黄91麻豆网站| 国产乱子伦一区二区=| 国产精品久久久久久搜索| 久久精品人人做人人| 亚洲成人免费看| 这里只有精品在线| 久久99久久无码毛片一区二区| 国产jizz| 毛片免费高清免费| 欧美精品v日韩精品v国产精品| 欧美三级视频网站| 国产又大又粗又猛又爽的视频| 国产丰满成熟女性性满足视频| 亚洲一级色| 欧美日韩激情| 亚洲一区二区日韩欧美gif| 久久免费精品琪琪| 久久综合色播五月男人的天堂| 国产精品网址在线观看你懂的| 91高清在线视频| 国产精品尹人在线观看| 国产欧美成人不卡视频| 国产玖玖视频| 亚洲最新地址| 国产精品亚洲综合久久小说| 色九九视频| 日韩在线2020专区| 午夜福利视频一区| 欧美日韩在线亚洲国产人| 一级毛片免费高清视频| 大陆精大陆国产国语精品1024 | 日本五区在线不卡精品| 日本欧美成人免费| 国产毛片基地| 老熟妇喷水一区二区三区| 尤物精品国产福利网站| 美女国内精品自产拍在线播放 | 日韩大片免费观看视频播放| 日a本亚洲中文在线观看| 国产成人欧美| 999在线免费视频| 国产成人无码AV在线播放动漫|