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

基于高階統計量的砂泥巖薄互層油氣檢測方法

2022-02-25 08:00:08陳偉劉威
長江大學學報(自科版) 2022年1期
關鍵詞:信號模型

陳偉,劉威

1.油氣資源與勘探技術教育部重點實驗室(長江大學),湖北 武漢 430100

2.長江大學地球物理與石油資源學院,湖北 武漢 430100

隨著油氣勘探程度的提高,新勘探目標更加隱蔽,而隱蔽性油氣藏也蘊含著潛在的、豐富的油氣資源。面對當前我國石油和天然氣對外依存度較高的國際環境,為了實現油氣自給的目標,加大對隱蔽性油氣藏的勘探開發勢在必行。砂泥巖薄互層儲層厚度薄,在埋藏較深的外部環境下,其反射信號非常微弱,一般很難直接在噪聲背景下識別出有效反射信號,因此砂泥巖薄互層油氣藏是一種典型的隱蔽性油氣藏。目前,砂泥巖薄互層油氣藏的勘探開發是世界級難題。從地震勘探的角度,識別薄互層層間多次波相互干涉形成的弱反射信號對砂泥巖薄互層油氣藏的勘探開發具有重要意義。

WIDESS[1]針對薄層的調諧效應開展了研究,具體給出了薄層的定義和識別薄層的極限厚度,這些研究為后續的薄層研究提供了理論認識。但是WIDESS只對單個薄層的厚度開展了研究,沒有考慮不同巖性的薄層相互疊加形成薄互層之后的振幅響應特征。KOEFOED等[2]基于褶積模型開展了薄互層的數值模擬研究,認為薄互層層間疊加后的振幅綜合響應與單個薄層的厚度之間存在一定的近似關系。KALLWEIT等[3]基于時域拾取進一步研究了薄層的厚度與地震波長之間的關系。BROWN[4]在流體均勻分布的假設情況下,進一步開展了薄層的厚度估計,但所提出的方法不適用于無井地區的資料。PARTYKA等[5]以及MARFURT等[6]基于離散傅里葉變換的譜分解方法對薄層的厚度進行了估計,并結合低頻陰影效應對含流體區域的定量刻畫提供了理論支持。由于離散傅里葉變換的時頻分析分辨率遠低于小波變換等更高級別的時頻分析工具的分辨率,因此后續又有眾多學者基于新的時頻分析算法開展了薄互層地震資料譜分解的研究。PORTNIAGUINE[7,8]和CASTAGNA[8],PURYEAR[9],CHOPRA等[10]和LI等[11]將反演理論和譜分解相結合,進一步提高了薄層的識別能力。PURYEAR等[12]定量分析了薄層厚度估計與譜反演之間的聯系,形成了薄互層定量研究的基礎理論。NOWAK等[13]結合AVO技術和互譜分解技術,首次證實可以識別厚度小于1/8波長的薄層,在一定的假設條件下,該方法可以識別任意厚度的薄層,但是假設條件過于苛刻,在實際資料處理中很難滿足。

上述學者對薄互層的研究側重于對薄層厚度的計算,很少對薄互層的地震響應特征進行研究,且大部分集中于數值模擬方法。薄互層地震數值模擬存在過于理想化、參數單一、入射條件特殊等不足,影響了研究結果的普遍性。相反,地震物理模型具有靈活、模擬條件真實以及約束條件少等特點,有利于對薄互層的地震響應特征展開研究[14]。地震波在地球介質中傳播時,經歷不同的地層會產生不同的地震響應[15,16];另外,地震信號是廣義高斯分布的時間序列,特別是薄互層儲層介質所產生的地震信號的高階譜特征具有非高斯性[17]。上述結論為薄互層地震信號的檢測和識別提供了堅實的基礎。從理論上說,地震信號既可以被認為是一種確定性信號,又可以被認為是一種隨機性信號。在實際工作中,由于地下的地質構造、介質速度和密度等參數的分布是未知的,地震波在穿過隨機分布的介質到達地表或地下的檢波器時,被記錄后的地震信號更多地呈現出或然性特征,而且地震信號是隨時間變化的,因此可以把地震信號看作隨機信號,把地震層位數據看成是隨機過程。根據隨機過程的理論,隨機過程的統計特征可以用隨機過程(時間序列)所有樣點的統計特性代替各個樣點的統計特性[18,19],在地球物理中常用的統計量是三階累積量和四階累積量[20-25]。

下面,筆者在改進和提高前人試驗技術的基礎上,制作了人工多孔砂巖薄層,觀測了含氣、含水和含油3種情況下薄互層的地震響應特征,并利用高階統計量(三階矩和三階累積量及Wigner雙譜)對相應特征進行了分析。

1 方法原理

1.1 高階統計量與高階譜

設x是隨機變量,則E(xk)(k=1,2,3,…)表示x的k階矩,其中E表示對隨機變量求取期望。

設隨機變量x的概率密度函數為f(x),則x的特征函數φ(v)為:

(1)

式中:v是任意的具有與隨機變量x相同分布的隨機變量。

將式(1)兩邊取自然對數并按泰勒級數展開,得到:

(2)

其中:

ψ(v)=lnφ(v)

(3)

(4)

式中:ψ(v)為累積量生成函數;Ck為k階累積量。

該次研究需要用到的是三階累積量,其物理意義是描述概率分布的非對稱性,即偏度。

設有長度為N的零均值非平穩時間序列{x(n)}(n=1,2,…,N),其時頻雙譜和三譜為:

B(n,ω1,ω2)=DFT{cum3[W(N)x(n)]}

(5)

T(n,ω1,ω2,ω3)=DFT{cum4[W(N)x(n)]}

(6)

式中:B(n,ω1,ω2)和T(n,ω1,ω2,ω3)分別表示在時刻n的雙譜和三譜;DFT表示離散傅里葉變換;cum3和cum4分別表示對序列求三階和四階累積量,即式(4)中的C3和C4;W(n)表示在時刻n的時窗。

該次研究在具體計算雙譜時,為了得到隨時間變化的光滑的高階譜,根據信號的實際情況設計了相應的時窗。對地震信號來說,應以覆蓋住地震子波為時窗的設計原則。在實際計算時,筆者采用的是Parzen時窗。

1.2 地震信號的高階譜時頻分析方法

地震信號是一種典型的非平穩和時變信號,其統計特性隨時間變化。常規的傅里葉變換只能單獨獲取地震信號時間域或者頻率域的振幅變化規律,而不能反映某個時刻某個頻率分量的振幅特征。而對于地震信號而言,往往要同時考慮時間和頻率信息,時頻分析工具正好可以完成這樣的任務。時頻分析能夠在二維平面上同時表現地震信號的時間、頻率和振幅信息。目前,在地震勘探中采用的時頻分析方法主要有:短時傅里葉變換法、帶通濾波法、表波變換法和復數道技術等。在上述方法中,除了復數道技術常用于地震剖面的特殊處理和解釋外,其他方法的應用并不廣泛。為了解決碳酸鹽巖儲層中廣泛發育的孔洞、裂縫、裂隙等引起的非均勻儲集體識別問題,很多人將高階譜分析方法引入到地震信號的時頻分析中,結果證明,高階譜時頻分析方法能夠較好地識別碳酸鹽巖非均勻儲集體。

設x(t)是一個確定的解析信號,k階的Wigner高階矩譜定義為[26]:

(7)

上述定義可以被解釋為一個k維局部矩函數mkx(t,τ1,τ2,…,τk)的k維傅里葉變換Wkx(t,ω1,ω2,…,ωk),即:

(8)

其中:

單個地震道的時頻雙譜是包含時間、兩個頻率共3個參數的四維數據體。該次研究在實現時頻雙譜時的流程為:①構造滑動時窗;②對每一窗口內的時間序列,先求取均值,用原序列減去均值得到近似的零均值序列;③求雙譜,對雙譜做切片處理,提取特征譜;④形成時頻雙譜剖面。

2 地震物理模型設計及數據采集和處理

2.1 模型設計

該次研究的物理模型是在借鑒A油田的油區資料后建立的,主要模擬埋深在2000m左右且層厚小于5m的砂泥巖薄互層。A油田的薄互層儲層段傾角較小,地層平緩,斷裂和背斜不發育,近似為水平層,在許多鉆井中儲層厚度小于5m。該物理模型由1個蓋層、7個砂巖層和8個泥巖層組成(見圖1)。模型中薄互層和蓋層上面有一個水層,模擬陸上勘探地下低速帶地層或海上勘探時的海水層,這樣可使得該次研究的研究結果更具有普遍性。同時,從室內物理實驗角度來看,這一低速水層有很好的耦合條件,可減少面波干擾及改進地震資料的可重復性,有利于地震監測的觀測。

圖1 地震物理模型示意圖Fig.1 Schematic diagram of seismic physical model

該模型中涉及到的薄互層走向沿著水平面,即薄互層中各層不傾斜,也沒有背斜和斷層等構造。模型中采用細小的SiO2砂粒(顆粒數為150~300目)與環氧樹脂攪拌混合的方法合成1~2mm厚人工砂巖。截至目前,在國內外已發表的物理模型實驗文獻中還沒見過使用薄至1~2mm厚的多孔材料。該次研究物理模型的主要參數設計為:①蓋層用環氧樹脂制作,層厚50.2mm,模擬地層厚度約為150m,速度2440m/s,密度1.2g/cm3;②薄互層由有機玻璃和人工砂層組成,有機玻璃片平均厚1.22mm,模擬地層厚度約為3.66m,速度約2500m/s,密度1.4g/cm3;人工砂層平均厚1.434mm,模擬地層厚度約為4.3m,砂層孔隙度約20%,砂層速度約為2600m/s,密度1.8g/cm3,7層薄互層總厚度約19.8mm,模擬地層厚度約為60m;③水層厚185mm,模擬地層厚度約為555m,速度約為1470m/s,密度1.0g/cm3。物理模型與野外地質體的參數對應關系如表1所示。

表1 模型與野外地質體的各變量對應關系

2.2 地震數據采集

整個試驗過程在水箱內進行,即把物理模型放入水中(見圖2),震源和接收器放在水面上按野外觀測方式移動,水面到模型頂面之間的水層看作是低速層。地震觀測類似于海上地震觀測,但震源不是使用空氣槍,而是使用美國PR5077可調脈沖發生器,震源頻率為105kHz(相當于野外實際地震觀測頻率35Hz)。薄互層含氣、含水和含油物理模型3次地震數據觀測時跨3個月,每次試驗的儀器條件完全一致,觀測系統保持一致,并確保每次試驗時模型的測線位置盡量一致。圖2給出了二維地震單邊觀測示意圖,炮點在右側,接收點從右向左移動。二維地震單邊觀測的觀測系統為0-60-1008m,參數如表2所示。

圖2 地震資料采集裝置示意圖Fig.2 Schematic diagram of seismic data acquisition device

表2 地震觀測系統參數表

具體觀測過程為:

1)第1次于某年11月5日進行了含氣模型地震觀測。往砂層中注入空氣后把模型放入真空恒溫箱內,用50℃溫度烘烤24h,然后密封四周進行地震觀測。

2)第2次于同年12月5日進行了含水模型地震觀測。往砂層中注入水,去掉四周密封,把模型放入盛水容器內后一起放在真空罐內,用抽真空法使薄互層全飽含水,反復抽真空和稱重以確定全飽和為止(7層薄互層總共注入186g水)。最后將模型放入水中進行地震觀測。

3)第3次于次年1月5日進行了含油模型地震觀測。含油狀態為砂層全飽和油。首先把砂層中的水盡量排出,方法與含氣狀態一樣,把模型放入真空恒溫箱內,用50℃溫度烘烤,每隔2h抽一次真空(歷時0.5h),24h后稱一次重量,以確定排出水的量。7層模型10d共排出約148g水,在后3d中幾乎很少有水排出(將薄互層內殘余的38g水全部排出的可能性很小,加上時間也不允許,所以只能讓這部分水留在模型內),最后浸入柴油中,用抽真空法飽和油,飽和油量174g。最后將模型四周密封后放入水中進行地震觀測。

2.3 地震物理模型數據處理及結果初步分析

對采集到的物理模型數據處理流程為:真振幅恢復→脈沖反褶積→速度分析→疊加→疊后時間偏移→帶通濾波。

該次研究設計的物理模型中的薄互層單層砂巖的厚度小于或近似等于1/8波長(約4m),處理后的結果(見圖3)表明,用直接的時差或振幅分析方法無法分辨薄互層中的單層,其原因是每個單層頂底反射系數是相反的,產生的反射信號時差很小,近似于反向疊加,信號疊加干涉結果是相互削弱。因此,圖3中的地震反射波是薄互層的整體地震響應,不論砂層的波阻抗差與圍巖相比是大、是小或接近,都不能分辨出單層界面的反射波,而薄互層段每個反射波中都包含了各單層反射能量的貢獻量。

圖3 疊后地震剖面Fig. 3 Post-stack seismic profile

從疊后地震剖面能夠看出,含不同流體的薄互層的地震響應不同(圖3中紅色方框所示):含氣地震剖面中薄互層頂層同向軸較穩定,但是頂層以下同向軸分布散亂;含水和含油地震剖面中薄互層的頂層和底層同向軸均較穩定,但是中間層有缺失的跡象。上述特征可以初步作為含氣薄互層和含水、含油薄互層的地震響應差別,但是僅僅根據該響應差別無法分辨出含水和含油薄互層。

3 薄互層油氣檢測

無論是高階矩還是高階累積量都是針對一維數據計算得到。該次研究將地震剖面的橫向層位(該次研究稱為橫向同向軸)信息看成是一維的隨機過程,利用高階統計量對其進行分析。圖4是原始地震剖面的放大圖,可以看出由于物理模型實驗的誤差造成含氣、含水和含油薄互層地震剖面中目標層位(薄互層頂層)的不一致。根據對應關系,選取目標層位的3組不同的橫向層位數據,分別是含氣地震剖面上的第40、80和120個層位以及含水、含油地震剖面上的第25、65和105個層位(圖4中黑色箭頭所示)。根據高階矩的定義分別計算各自的三階矩(見圖5),其中圖5(a)~(c)分別是含氣地震剖面上3個橫向層位數據對應的三階矩結果,圖5(d)~(f)對應含水地震剖面的結果,圖5(g)~(i)對應含油地震剖面的結果。分析得出,含氣薄互層的橫向同向軸的三階矩能量比較分散,盡管在中心處有能量集中,但是大部分能量分散在水平對稱軸、垂直對稱軸以及對角線位置;而含水和含油薄互層的橫向同向軸的三階矩能量比較集中,大部分能量集中在中心位置,能量沿水平對稱軸、垂直對稱軸以及對角線有發散的跡象,但是很弱。另外需要指出的是,含氣、含水和含油薄互層對應的三階矩剖面都具有穩定的形態,這一點可以作為識別含氣、含水與含油薄互層的標志。

將原始地震剖面的放大圖(圖4)與相應的三階矩剖面(圖5)進行對比可以看出,含氣地震剖面的第40個、80個和120個層位比較散亂,在視覺上的體現是含噪較強,偽波較多,其相應的三階矩剖面(圖5(a)~(c))中的能量團較多且沿水平對稱軸、垂直對稱軸以及對角線3個方向進行分布;而含水和含油地震剖面的第25和125個層位均較穩定,其對應的三階矩剖面(圖5(d)和(f)及圖5(g)和(i))中能量團少且非常集中(最明顯的是居中的大能量團),且沿3個方向的能量發散不明顯;第65個層位盡管比較散亂且不連續,但是其對應的三階矩剖面(圖5(e)和(h))依然呈現出能量大部分集中在中心的特點。

圖4 地震記錄中目標層(薄互層)放大圖Fig. 4 Enlarged view of the target layer (thin interbedded layer) in seismic records

圖5 橫向同向軸的三階矩Fig. 5 Third-order moment of transverse isotropic axis

三階矩剖面能夠很好地將含氣和含水、含油薄互層的特征區分開,試著計算3種疊后剖面的三階累積量。圖6是根據式(4)計算得到的原始疊后地震記錄的三階累積量圖,可以看出,3種不同流體對應的三階累積量分布形態比較對稱,整個能量背景比較平緩。含氣和含水、含油地震剖面對應的三階累積量剖面均很有規律且基本類似,即都是以對角線為界兩邊成波浪形片狀起伏,每個波浪平面上的數值又各有起伏,認為是同一形態的三階累積量。盡管如此,依然能夠看到含氣(圖6(a))與含水(圖6(b))、含油剖面(圖6(c))的區別。從圖6各分圖對比可以看到,第17道至第22道高值區的范圍依次增大,差別最大的區域在坐標點[25,25]附近,表現為從第17道的接近零值到第22道整體呈低頻狀峰值凸起。另外,從圖6還可以發現不同流體的三階累積量在時移量最大的地方均存在明顯的低值區,即在坐標點[-25,-25]和[25,25]兩處出現最大時移量。從整個三階累積量圖上(圖6)難以觀察具體的細節,對合成地震信號的三階累積量沿著時移平面做對角線的切片,得到如圖7所示的曲線。圖7清晰地顯示了含油剖面的三階累積量剖面的坐標點[25,25]處存在一個明顯的低值(如圖7(c)黑色箭頭所示),該特征在圖6中被掩蓋不易發現。另外,從圖7還能很明顯地看出含氣地震記錄的同向軸的三階累積量沿對角線的切片呈余弦曲線分布,而含水和含油地震記錄呈正弦曲線分布。

圖6 原始疊后地震記錄的三階累積量Fig. 6 Third-order cumulants of the original post-stack seismic records

圖7 地震記錄三階累積量沿對角線的切片Fig. 7 Diagonal slice of third-order cumulants of seismic records

圖8是含氣、含水和含油地震記錄的雙譜分布圖,雙譜的幅值主要分布在頻率點坐標的中間區域,近似呈柱狀形態,這些數據分布也呈尖脈沖狀,相互之間沒有平滑的過渡,連續性差。圖8中含氣、含水和含油數據的雙譜都呈束狀分布,形態比較類似,在束狀雙譜的根部包裹較多數值偏低的峰值,從下往上逐漸減少,頂部的峰值大小存在差別。對比含氣、含水和含油雙譜圖可以發現,含氣地震記錄的雙譜中的峰值更參差不齊,而含水和含油地震記錄的雙譜峰值比較均一。因此可以得到如下認識:計算的含水、含油及含氣地震記錄的雙譜基本類似,但是在某些方面也存在一定區別,主要表現在雙譜峰值的分布特點上。對比三階累積量和雙譜可以發現,在分析不同類型含流體薄互層的單道數據的高階統計量特征時,三階累積量相比雙譜能夠更好地體現出數據相互之間的差異。

圖8 地震記錄雙譜圖Fig. 8 Bispectrum of seismic records

從地震記錄(圖4)上可以清晰地看出,薄互層附近存在大量的尾波噪聲(該噪聲屬于高斯型噪聲),嚴重影響了實際結果的分析,而高階統計量時頻分析方法具有分析非高斯過程的能力,可以壓制高斯型噪聲,對于較低信噪比的地震記錄有很強的適用性,其中Wigner方法最具代表性。由于多次散射波一致性很差,而高階Wigner分布能自動抑制一部分噪聲的影響并增強時頻分布的集聚性,所以選擇高階Wigner 分布來進行分析是非常合理的。為了保證數據精確度,該次研究采用了Wigner三譜方法。

圖9是圖4中目標層位對應的Wigner三譜時頻分布圖。從圖中可以看出,計算出的Wigner三譜時頻分布圖與原始疊后地震剖面相比,能夠精確地確定薄互層的頂層和底層的縱向位置,即時頻分布圖中紅色所示區域,說明高階譜時頻分布方法能夠找出地震剖面中蘊含的由非均質介質引起的地震信息。地震信號是一種廣義高斯分布的信號,其Wigner三譜值較小,而由砂巖薄互層產生的地震波主要表現為低頻的、相位畸變的非線性特性,是一種相對于正常的廣義高斯分布的地震信號的異常信號,其Wigner三譜值較大。因此,利用Wigner三譜將能夠較容易地識別出由砂巖薄互層等非均質產生的低頻和相位畸變的地震波,從而獲知薄互層在地下的分布區域。在Wigner三譜時頻圖上也能較容易地分辨出含氣和含水含油薄互層的差異,即含氣薄互層的響應在縱向上不連續,在橫向上分辨率更低,且頻率范圍在負值區域內(圖中紅色虛線以左);含水和含油薄互層的響應在縱向上非常連續,在橫向上分辨率更高,且頻率范圍在正值區域內(圖中紅色虛線以右)。需要說明的是,通過Wigner三譜時頻分布圖依然不能區分含水和含油薄互層。

圖9 圖4中目標層位對應的Wigner三譜時頻分布圖Fig. 9 The time-frequency distribution of Wigner trispectrum corresponding to the target layer in FIG.4

4 結論

1)具有廣義高斯分布的地震信號,當記錄中混有由非均勻介質(例如薄互層等)導致的非高斯的或非線性的低頻信號時,能夠利用廣義高斯分布地震時間序列的基本性質,用隨時間變化的時頻雙譜分析、Wigner時頻雙譜分析等方法檢測出這些異常信號。

2)該次研究旨在證明高階統計量方法是一種高效的檢測儲層非均質性的方法。該方法不需要井資料的約束,方法簡單,可以快速實現且效果明顯。尤其是對于薄互層這種具有橫向非均質性的介質,利用高階統計量方法來進行油氣的識別具有很好的效果。

3)當應用高階統計量分析常見異常數據時,要從不同的角度即各階累積量、雙譜以及Wigner三譜等進行全方位研究,因為不同的異常現象(含氣、含水和含油)會具有類似的某一種高階統計量特征,但是不同現象之間又會存在另外一種高階統計量的差異。以該原則為基礎,該次研究通過分析單道地震映像數據高階統計量特征,對比分析并歸納總結了在勘探地球物理中薄互層這種非均質介質的各類高階統計量分布特征。高階統計量作為一種輔助手段,在資料解釋時可減少人為的感性判斷,提高了解釋的準確性。

4)含氣薄互層與含水、含油薄互層存在明顯的地震響應差異;三階矩是區分含水和含油薄互層的關鍵統計量,而用其他高階統計量無法區分這兩者的差異。

致謝:文中物理模型實驗在中國石油大學(北京)CNPC物探重點實驗室完成。

猜你喜歡
信號模型
一半模型
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
重要模型『一線三等角』
完形填空二則
重尾非線性自回歸模型自加權M-估計的漸近分布
孩子停止長個的信號
3D打印中的模型分割與打包
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 精品久久久久久成人AV| 72种姿势欧美久久久大黄蕉| 国产精品亚洲一区二区在线观看| 亚洲综合激情另类专区| 中文国产成人精品久久| 亚洲欧洲日产无码AV| 高h视频在线| 亚洲视频无码| 国产一区二区三区日韩精品| 四虎国产精品永久在线网址| 亚洲第一极品精品无码| 日韩精品无码不卡无码| 四虎成人免费毛片| 99re热精品视频中文字幕不卡| 色爽网免费视频| 国产99视频精品免费视频7| 亚洲二区视频| 高清大学生毛片一级| 亚洲精品无码久久毛片波多野吉| 在线欧美国产| 91视频99| 99久久精彩视频| 国产黄色片在线看| 五月婷婷导航| 国内精品视频| 国产第一页免费浮力影院| 亚洲一区无码在线| 国产黄色视频综合| 亚洲男人在线天堂| 国产一区三区二区中文在线| 国产无码高清视频不卡| 在线观看91精品国产剧情免费| 草草影院国产第一页| 色综合成人| 亚洲,国产,日韩,综合一区| 久久精品只有这里有| 少妇精品在线| 2020精品极品国产色在线观看 | 一本综合久久| 欧美三级日韩三级| 日韩欧美高清视频| a天堂视频| 亚洲精品无码抽插日韩| 成人午夜网址| 9久久伊人精品综合| 青草视频久久| 999福利激情视频 | 美女无遮挡免费视频网站| 国产微拍一区| 午夜国产理论| 国产欧美视频在线| 9cao视频精品| 国产成人亚洲精品蜜芽影院| 污网站免费在线观看| 毛片基地美国正在播放亚洲 | 青草国产在线视频| 人妻少妇乱子伦精品无码专区毛片| 亚洲,国产,日韩,综合一区| 成人韩免费网站| 91视频国产高清| 97超碰精品成人国产| 高清无码手机在线观看| 欧美日本在线播放| 国产99热| 日韩无码黄色| 亚洲婷婷丁香| 99精品一区二区免费视频| 99久久免费精品特色大片| 噜噜噜久久| 亚洲一区二区三区国产精品 | 免费观看亚洲人成网站| 国产精品永久免费嫩草研究院| 91原创视频在线| 成人蜜桃网| 91人妻在线视频| 国产自无码视频在线观看| 成人蜜桃网| 欧美精品伊人久久| 全部免费毛片免费播放 | 九色视频一区| 九九九久久国产精品| 在线视频亚洲色图|