亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        最小二乘疊前時(shí)間偏移在地震數(shù)據(jù)規(guī)則化中的應(yīng)用

        2019-01-25 10:31:38龔仁彬王從鑌胥小馬吳海莉
        石油地球物理勘探 2019年1期
        關(guān)鍵詞:規(guī)則化算子反演

        吳 丹 龔仁彬 王從鑌 胥小馬 吳海莉

        (中國(guó)石油勘探開發(fā)研究院西北分院計(jì)算機(jī)所,北京 100191)

        0 引言

        在三維地震勘探中,受地表障礙物、觀測(cè)設(shè)備、施工條件等因素的限制,采集到的數(shù)據(jù)往往出現(xiàn)不規(guī)則、缺道、壞道等現(xiàn)象。地震數(shù)據(jù)的這種不規(guī)則性會(huì)嚴(yán)重影響后續(xù)的多次波壓制[1]、偏移成像[2]、四維地震監(jiān)測(cè)[3]等環(huán)節(jié)。為了消除地震數(shù)據(jù)不規(guī)則的影響,一般采用地震數(shù)據(jù)規(guī)則化方法重建缺失數(shù)據(jù)。但是,由于數(shù)據(jù)規(guī)則化本質(zhì)上是一個(gè)欠定的反演問(wèn)題,因此一直以來(lái)都是地震數(shù)據(jù)處理領(lǐng)域十分具有挑戰(zhàn)性的研究課題。

        目前,地震數(shù)據(jù)規(guī)則化方法主要分為三類:第一類是在信號(hào)變換域引入稀疏約束恢復(fù)地震數(shù)據(jù),如Radon變換[4]、傅里葉變換[5-8]、小波變換[9]、曲波變換[10-11]等;第二類是基于褶積算子的預(yù)測(cè)濾波器方法,根據(jù)已知數(shù)據(jù)提取局部地震屬性(如局部?jī)A角),再擴(kuò)展該信息以恢復(fù)缺失數(shù)據(jù)[12-18];第三類是模型驅(qū)動(dòng)方法,利用已知的地下模型信息,采用積分延拓算子重建缺失數(shù)據(jù),常見的延拓算子包括炮檢距延拓[19]、炮點(diǎn)延拓[20]、方位角延拓[21]等。其中前兩類方法屬于數(shù)據(jù)驅(qū)動(dòng)的規(guī)則化方法,

        以往基于模型驅(qū)動(dòng)的數(shù)據(jù)規(guī)則化方法雖然利用了已知的地下模型信息,但是處理精度往往過(guò)度依賴于模型的準(zhǔn)確性。當(dāng)?shù)叵履P痛嬖谝欢ㄕ`差時(shí),規(guī)則化結(jié)果不理想。本文采用疊前時(shí)間偏移和反偏移算子作為積分延拓算子,在成像空間保留了一個(gè)維度的冗余信息(如地表炮檢距),使映射到成像空間的過(guò)程中不丟失模型的誤差信息; 再通過(guò)引入更加穩(wěn)定的預(yù)條件算子,可以在成像空間更好地提取和應(yīng)用反射波同相軸特征,從而改善數(shù)據(jù)規(guī)則化質(zhì)量。最后通過(guò)理論模型和實(shí)際資料處理驗(yàn)證本文方法能否夠有效改善數(shù)據(jù)規(guī)則化效果,改進(jìn)最終成像質(zhì)量。

        1 方法原理

        1.1 最小二乘偏移理論

        基于一階Born近似理論,反射地震數(shù)據(jù)d可以通過(guò)地下反射系數(shù)模型m線性映射得到

        d=Lm

        (1)

        式中L為線性映射算子,也稱為Born正演算子,可以是疊前/疊后反時(shí)間偏移算子,也可以是疊前/疊后反深度偏移算子。該算子與速度模型有關(guān),由于規(guī)則化過(guò)程中速度模型固定,因此式(1)忽略了速度這個(gè)參數(shù)。最小二乘偏移最小化目標(biāo)函數(shù)為

        (2)

        式中:dobs為觀測(cè)數(shù)據(jù);D為正則化算子;μ為阻尼因子。式(2)的顯式解為

        m=(LTL+μDTD)-1LTdobs

        (3)

        定義一個(gè)新的向量p,滿足

        m=Pp

        (4)

        則式(2)變?yōu)?/p>

        (5)

        式中P=D-1為預(yù)條件算子。上式的解為

        p=(PTLTLP+μI)-1PTLTdobs

        (6)

        式中I為單位矩陣。將式(6)代入式(4),可得

        m=P(PTLTLP+μI)-1PTLTdobs

        (7)

        Clearbout[22]指出,合理的預(yù)條件算子比對(duì)應(yīng)的正則化算子可以更快地加速迭代反演的收斂速度。

        由于Born正演(或反偏移)算子L一般無(wú)法用矩陣顯式表達(dá),而且模型和數(shù)據(jù)的規(guī)模很大,因此一般采用共軛梯度法近似求解式(7),得到分辨率更高、振幅保真度更好的反射系數(shù)模型。在相同的理論框架下,選擇不同的反偏移算子,可以構(gòu)建不同的最小二乘偏移算法,如最小二乘疊前時(shí)間偏移[23-24]、最小二乘逆時(shí)偏移[25-26]等。

        1.2 地震數(shù)據(jù)規(guī)則化

        當(dāng)觀測(cè)數(shù)據(jù)不規(guī)則時(shí),式(5)可以修改為

        (8)

        式中K為Mask算子,該算子的作用是將缺失地震道置0,非缺失地震道保持不變。通過(guò)共軛梯度法最小化目標(biāo)函數(shù)式(8)得到最優(yōu)解popt,再應(yīng)用預(yù)條件算子和反偏移算子,則可以重建缺失地震數(shù)據(jù)

        dreg=LPpopt

        (9)

        重建后的地震數(shù)據(jù)可以表示為

        dreg=LP(PTLTKTKLP+μI)-1PTLTKTdobs

        (10)

        分三種情況分析式(10)中各項(xiàng)對(duì)數(shù)據(jù)規(guī)則化的影響。

        第一種情況,忽略Hessian矩陣求逆過(guò)程(式(10)中的求逆項(xiàng)),即直接采用偏移+反偏移進(jìn)行數(shù)據(jù)重建,式(10)變?yōu)?/p>

        dreg=LLTKTdobs

        (11)

        這種做法雖然也能夠起到數(shù)據(jù)插值的作用,但是振幅信息不可靠。通過(guò)設(shè)計(jì)一個(gè)匹配濾波器可以在一定程度上解決這個(gè)問(wèn)題,但是匹配濾波器的估計(jì)本身也是一個(gè)欠定的反問(wèn)題,存在多解性。

        第二種情況,忽略預(yù)條件項(xiàng),即直接采用最小二乘偏移+反偏移進(jìn)行數(shù)據(jù)重建,式(10)變?yōu)?/p>

        dreg=L(LTKTKL)-1LTKTdobs

        (12)

        這種做法可以解決偏移+反偏移中存在的振幅匹配問(wèn)題,但是無(wú)法消除數(shù)據(jù)缺失帶來(lái)的成像噪聲,該噪聲反投影到數(shù)據(jù)域,會(huì)降低插值數(shù)據(jù)的信噪比。

        第三種情況,完整地求解式(10),即采用預(yù)條件最小二乘偏移+反偏移進(jìn)行數(shù)據(jù)重建,這種做法既考慮了振幅不匹配的問(wèn)題,又可以壓制數(shù)據(jù)不規(guī)則帶來(lái)的噪聲。

        實(shí)際上,預(yù)條件算子的選取對(duì)數(shù)據(jù)規(guī)則化效果的影響很大。本文采用因果積分算子作為預(yù)條件算子,該算子的作用是使振幅隨炮檢距平滑變化。當(dāng)速度比較準(zhǔn)確時(shí),共成像點(diǎn)道集拉平,滿足振幅平滑變化的條件;當(dāng)速度誤差較大時(shí),相鄰道之間的振幅差異也不會(huì)太大,因此也基本能滿足條件;當(dāng)速度誤差很大時(shí),需要進(jìn)一步采用更加合適的預(yù)條件算子,以保留原始數(shù)據(jù)中的剩余旅行時(shí)信息。

        1.3 疊前時(shí)間偏移/反偏移

        綜合考慮算子的靈活性、計(jì)算效率以及預(yù)處理階段對(duì)規(guī)則化的需求,本文采用Kirchhoff疊前時(shí)間反偏移作為線性正演算子L,該算子可以表示為

        (13)

        式中:m(x,h,τ)為時(shí)間偏移共炮檢距道集;d(y,H,t)為疊前共炮檢距數(shù)據(jù);x和h分別為成像域中心點(diǎn)坐標(biāo)和炮檢距;τ為成像域雙程旅行時(shí);y和H分別為數(shù)據(jù)域中心點(diǎn)坐標(biāo)和炮檢距,且H∈[h-Δh,h+Δh],Δh為炮檢距間隔;w為加權(quán)函數(shù),合適的加權(quán)函數(shù)可以加快反演的收斂速度。在各向同性介質(zhì)和直射線假設(shè)下,雙程旅行時(shí)滿足雙平方根方程

        (14)

        式中v為均方根速度。

        正演算子式(13)的共軛轉(zhuǎn)置LT為Kirchhoff疊前時(shí)間偏移,其表達(dá)式為

        (15)

        1.4 速度的影響

        疊前時(shí)間偏移共成像點(diǎn)道集屬于擴(kuò)展模型的一種,該模型除了和物理位置有關(guān),還允許模型隨炮檢距變化。因此,不論背景速度是否正確,使用該模型均可以模擬出和觀測(cè)數(shù)據(jù)相同的反射波旅行時(shí)信息,都可以在不發(fā)生周波跳躍的前提下正確地匹配波形信息[27],即不論速度正確與否,隨著迭代次數(shù)的增加,目標(biāo)函數(shù)式(2)中的第一項(xiàng)(數(shù)據(jù)擬合項(xiàng))總是可以逐漸趨近于零。但是,如果速度模型不準(zhǔn)確,式(2)中的正則化算子需要適應(yīng)速度不準(zhǔn)確帶來(lái)的旅行匹配誤差。

        2 數(shù)值試驗(yàn)

        2.1 Marmousi模型測(cè)試

        圖1a為深度域Marmousi層速度模型,圖1b為Dix公式轉(zhuǎn)換后的時(shí)間域均方根速度模型。圖1c是公開的波動(dòng)方程有限差分正演共炮點(diǎn)道集[28]經(jīng)過(guò)炮檢距分選后的(時(shí)間,中心點(diǎn),炮檢距)域的疊前數(shù)據(jù)體。圖1d為圖1c數(shù)據(jù)隨機(jī)剔除50%地震道后的疊前數(shù)據(jù)體。

        圖1 Marmousi模型及其正演數(shù)據(jù)體

        使用均方根速度場(chǎng)(圖1b)和不同的疊前數(shù)據(jù)體進(jìn)行疊前時(shí)間偏移和最小二乘疊前時(shí)間偏移試驗(yàn)。首先對(duì)比炮檢距為-0.325km的共炮檢距成像剖面。圖2a和圖2b為使用完整的地震數(shù)據(jù)(圖1c)分別進(jìn)行常規(guī)疊前時(shí)間偏移和最小二乘疊前時(shí)間偏移得到的共炮檢距成像剖面,可以看出:盡管Marmousi模型構(gòu)造復(fù)雜,速度橫向變化劇烈,但是時(shí)間偏移仍然可以大致準(zhǔn)確地對(duì)地下構(gòu)造形態(tài)進(jìn)行成像;經(jīng)過(guò)最小二乘反演之后的成像結(jié)果分辨率得到顯著提高。圖2c、圖2d和圖2e為使用不完整數(shù)據(jù)(圖1d)分別進(jìn)行常規(guī)疊前時(shí)間偏移、最小二乘疊前時(shí)間偏移和預(yù)條件最小二乘疊前時(shí)間偏移得到的共炮檢距成像剖面,可以看出,數(shù)據(jù)缺失給成像結(jié)果帶來(lái)許多噪聲,降低了信噪比,經(jīng)過(guò)最小二乘反演,雖然可以提高分辨率,但是并不能消除數(shù)據(jù)缺失的影響。通過(guò)加入合理的先驗(yàn)信息(比如本文所采用的共成像點(diǎn)道集同相軸沿炮檢距平滑變化),既可以提高成像分辨率,又可以壓制成像噪聲。

        圖3a和圖3b為采用完整數(shù)據(jù)分別進(jìn)行常規(guī)疊前時(shí)間偏移和最小二乘疊前時(shí)間偏移的共成像點(diǎn)道集,可以看出:盡管模型復(fù)雜,速度橫向變化劇烈,但是采用Dix公式轉(zhuǎn)化而來(lái)的均方根速度場(chǎng)(圖1b),共成像點(diǎn)道集基本可以被拉平,該信息可以作為先驗(yàn)信息約束反演成像;另外,經(jīng)過(guò)最小二乘反演,成像分辨率有一定的提升,而且淺層大炮檢距的拉伸效應(yīng)也在一定程度上得到壓制。圖3c、圖3d、圖3e為采用不完整地震數(shù)據(jù)分別進(jìn)行常規(guī)疊前時(shí)間偏移、最小二乘疊前時(shí)間偏移和預(yù)條件最小二乘疊前時(shí)間偏移的共成像點(diǎn)道集,可以看出,和共炮檢距成像剖面類似,數(shù)據(jù)缺失在成像道集中也引入噪聲,通過(guò)最小二乘反演能夠提高分辨率,但是無(wú)法壓制這類噪聲,利用共成像點(diǎn)道集同相軸沿炮檢距平滑變化這個(gè)先驗(yàn)信息,再進(jìn)行預(yù)條件最小二乘疊前時(shí)間偏移,不僅能夠提高分辨率,而且能夠有效壓制數(shù)據(jù)缺失帶來(lái)的噪聲。

        (a)完整數(shù)據(jù)常規(guī)疊前時(shí)間偏移;

        (b)完整數(shù)據(jù)最小二乘疊前時(shí)間偏移;

        (c)不完整數(shù)據(jù)常規(guī)疊前時(shí)間偏移;

        (d)不完整數(shù)據(jù)最小二乘疊前時(shí)間偏移;

        (e)不完整數(shù)據(jù)預(yù)條件最小二乘疊前時(shí)間偏移

        (a)完整數(shù)據(jù)常規(guī)疊前時(shí)間偏移;

        (b)完整數(shù)據(jù)最小二乘疊前時(shí)間偏移;

        (c)不完整數(shù)據(jù)常規(guī)疊前時(shí)間偏移;

        (d)不完整數(shù)據(jù)最小二乘疊前時(shí)間偏移;

        (e)不完整數(shù)據(jù)預(yù)條件最小二乘疊前時(shí)間偏移

        對(duì)上述成像結(jié)果應(yīng)用疊前時(shí)間反偏移即可以進(jìn)行數(shù)據(jù)規(guī)則化。圖4為5.300km位置處的CMP道集規(guī)則化結(jié)果,圖5為規(guī)則化結(jié)果與原始道集誤差。

        圖4 Marmousi模型不同方法數(shù)據(jù)規(guī)則化的5.300km處CMP道集對(duì)比

        圖5 Marmousi模型不同方法數(shù)據(jù)規(guī)則化的5.300km處CMP道集誤差對(duì)比

        由圖4和圖5可以看出,直接采用常規(guī)偏移方法也可以對(duì)數(shù)據(jù)進(jìn)行插值,但是,由于沒有考慮Hessian矩陣的影響,插值結(jié)果的振幅信息不準(zhǔn)確,需要考慮振幅匹配問(wèn)題,而這個(gè)問(wèn)題也是一個(gè)欠定問(wèn)題,存在多解性。通過(guò)最小二乘反演可以解決振幅匹配問(wèn)題,但是在數(shù)據(jù)缺失嚴(yán)重的區(qū)域,沒有先驗(yàn)信息做約束,容易將成像結(jié)果中的噪聲反投影到數(shù)據(jù)空間,使插值結(jié)果信噪比較低。預(yù)條件最小二乘偏移規(guī)則化可以同時(shí)解決振幅不匹配和信噪比的問(wèn)題,使得插值結(jié)果振幅屬性更加可靠、信噪比更高。

        需要注意的是,對(duì)于本文方法,如果均方根速度比較準(zhǔn)確,可以采用沿著炮檢距方向的平滑算子作為預(yù)條件算子,但是如果均方根速度不準(zhǔn)確,則需要采用更加復(fù)雜的預(yù)條件算子,比如雙曲Radon變換、非穩(wěn)態(tài)平面波構(gòu)建濾波器[8]等。盡管如此,通過(guò)偏移算子的作用,成像空間依然會(huì)比數(shù)據(jù)空間更簡(jiǎn)單,因此有助于簡(jiǎn)化預(yù)條件算子的設(shè)計(jì),這也是模型驅(qū)動(dòng)類數(shù)據(jù)規(guī)則化方法的意義所在。

        2.2 實(shí)際資料處理

        為了進(jìn)一步驗(yàn)證本文方法的有效性,對(duì)海上實(shí)際不規(guī)則地震數(shù)據(jù)進(jìn)行規(guī)則化處理。圖6a為海上實(shí)際觀測(cè)的疊前數(shù)據(jù)體,受拖纜漂移等因素的影響,觀測(cè)數(shù)據(jù)中存在明顯的地震道缺失;圖6b為疊加速度分析得到的速度場(chǎng),本文使用該速度場(chǎng)進(jìn)行時(shí)間域成像,并將其作為模型先驗(yàn)信息進(jìn)行數(shù)據(jù)規(guī)則化處理。

        使用圖6b所示的速度場(chǎng),采用常規(guī)疊前時(shí)間偏移、最小二乘疊前時(shí)間偏移、預(yù)條件最小二乘疊前時(shí)間偏移進(jìn)行成像。圖7是炮檢距為-0.425km的三種方法共炮檢距成像剖面,可以看出,受數(shù)據(jù)不規(guī)則的影響,常規(guī)偏移結(jié)果中存在明顯的成像噪聲。經(jīng)過(guò)最小二乘反演,成像分辨率得到顯著提高,但是信噪比依然較低。通過(guò)在共成像點(diǎn)道集引入平滑約束,進(jìn)行預(yù)條件最小二乘反演之后,不僅可以提高成像分辨率,而且可以有效壓制成像噪聲。

        圖8為常規(guī)疊前時(shí)間偏移、最小二乘疊前時(shí)間偏移、預(yù)條件最小二乘疊前時(shí)間偏移后得到的位于2.360km處的炮檢距域共成像點(diǎn)道集,可以看出,共成像點(diǎn)道集基本拉平,說(shuō)明圖6b所示的速度場(chǎng)比較準(zhǔn)確。另外,也能夠明顯的看出,受數(shù)據(jù)缺失的影響,常規(guī)疊前時(shí)間偏移成像道集的信噪比較低(圖8a)。

        圖6 海上實(shí)際觀測(cè)疊前數(shù)據(jù)體(a)及其疊加速度場(chǎng)(b)

        (a)常規(guī)疊前時(shí)間偏移;

        (b)最小二乘疊前時(shí)間偏移;

        (c)預(yù)條件最小二乘疊前時(shí)間偏移

        圖7 海上實(shí)際數(shù)據(jù)炮檢距為-0.425m的不同方法共炮檢距偏移剖面

        通過(guò)最小二乘反演進(jìn)行成像,雖然可以提高分辨率,但是無(wú)法壓制數(shù)據(jù)缺失帶來(lái)的噪聲(圖8b)。進(jìn)一步在成像道集引入同相軸沿炮檢距方向平滑變化這一先驗(yàn)信息,可以顯著提高道集的信噪比,同時(shí)也能夠保持反演帶來(lái)的分辨率(圖8c)。

        最后,應(yīng)用疊前時(shí)間反偏移對(duì)缺失數(shù)據(jù)進(jìn)行重建。圖9a為2.360km處的CMP道集,該道集存在明顯的數(shù)據(jù)缺失。圖9b~圖9d分別為采用常規(guī)疊前時(shí)間偏移、最小二乘疊前時(shí)間偏移和預(yù)條件最小二乘疊前時(shí)間偏移的規(guī)則化結(jié)果。從結(jié)果可以看出,對(duì)常規(guī)疊前時(shí)間偏移結(jié)果進(jìn)行反偏移,雖然也可以起到數(shù)據(jù)插值的作用,但是由于沒有考慮Hessian矩陣的影響,插值出來(lái)的振幅信息不可靠。經(jīng)過(guò)最小二乘反演,雖然可以解決振幅信息不匹配的問(wèn)題,但是受不規(guī)則數(shù)據(jù)引入噪聲的影響,插值結(jié)果的信噪比較低。通過(guò)預(yù)條件化技術(shù)既可以恢復(fù)振幅信息,又能夠提高規(guī)則化后數(shù)據(jù)的信噪比。

        (a)常規(guī)疊前時(shí)間偏移;

        (b)最小二乘疊前時(shí)間偏移;

        (c)預(yù)條件最小二乘疊前時(shí)間偏移

        圖9 海上實(shí)際數(shù)據(jù)不同方法數(shù)據(jù)規(guī)則化的CMP道集

        3 結(jié)論與討論

        本文在最小二乘反演的理論框架下,采用最小二乘疊前時(shí)間偏移進(jìn)行數(shù)據(jù)規(guī)則化,利用合理的預(yù)條件算子改善數(shù)據(jù)規(guī)則化的處理效果。從理論模型和實(shí)際資料處理結(jié)果可以得到如下幾點(diǎn)結(jié)論及認(rèn)識(shí):

        (1)疊前時(shí)間偏移共成像點(diǎn)道集是一種擴(kuò)展模型,可以緩解模型驅(qū)動(dòng)類數(shù)據(jù)規(guī)則化方法對(duì)模型的過(guò)度依賴;

        (2)最小二乘疊前時(shí)間偏移可以提高成像分辨率,改善振幅保真度,是本文規(guī)則化方法的中間結(jié)果,可以用作質(zhì)量監(jiān)控和后續(xù)的成像結(jié)果分析;

        (3)成像空間比數(shù)據(jù)空間的信息更加簡(jiǎn)單,更有利于設(shè)計(jì)預(yù)條件算子;

        (4)合理的預(yù)條件算子可以改善規(guī)則化的效果。當(dāng)然,在復(fù)雜地表、復(fù)雜地質(zhì)條件下,也可以采用更加復(fù)雜的成像算子,以簡(jiǎn)化成像空間預(yù)條件算子的設(shè)計(jì),從而進(jìn)一步改善規(guī)則化效果,但是相應(yīng)地計(jì)算成本會(huì)增加,對(duì)速度模型的依賴性也會(huì)變強(qiáng)。

        猜你喜歡
        規(guī)則化算子反演
        反演對(duì)稱變換在解決平面幾何問(wèn)題中的應(yīng)用
        擬微分算子在Hp(ω)上的有界性
        論黨的民主集中制的規(guī)則化
        理論探索(2021年6期)2021-05-26 17:06:22
        各向異性次Laplace算子和擬p-次Laplace算子的Picone恒等式及其應(yīng)用
        數(shù)據(jù)規(guī)則化技術(shù)的研究及應(yīng)用
        數(shù)據(jù)規(guī)則化技術(shù)在焉耆盆地老資料處理中的應(yīng)用
        一類Markov模算子半群與相應(yīng)的算子值Dirichlet型刻畫
        采用移動(dòng)掃描成像法實(shí)現(xiàn)高壓氣瓶?jī)?nèi)壁凹坑缺陷檢測(cè)
        基于低頻軟約束的疊前AVA稀疏層反演
        基于自適應(yīng)遺傳算法的CSAMT一維反演
        国产黑色丝袜在线观看视频| 亚洲精品久久久久avwww潮水| 亚洲色欲久久久综合网| 日韩在线看片免费人成视频| 一本大道久久精品一本大道久久| 丰满老熟女性生活视频| 国产日韩精品欧美一区喷水| 日本高清aⅴ毛片免费| 97超级碰碰碰久久久观看| 天堂久久一区二区三区| 无套中出丰满人妻无码| 国产露脸精品产三级国产av| 高潮喷水无遮挡毛片视频| av高潮一区二区三区| 99久久久无码国产精品秋霞网| 又硬又粗又大一区二区三区视频| 99热这里只有精品久久6| 久久人妻中文字幕精品一区二区| 美女高潮黄又色高清视频免费 | 久久精品片| 久久亚洲精品一区二区| 亚洲一区二区国产激情| 国产办公室沙发系列高清| 国产精品久久1024| 偷拍熟女露出喷水在线91| 在线无码中文字幕一区| 国自产偷精品不卡在线| 国色天香精品亚洲精品| 亚洲不卡在线免费视频| 亚洲中文字幕无码一久久区| 一本色道久久综合亚洲精品小说| 一区二区三区精品偷拍| 美女国产毛片a区内射| 亚洲乱亚洲乱少妇无码99p| 人妻无码中文专区久久AV| 亚洲av高清一区二区在线观看| 小蜜被两老头吸奶头在线观看| 乱子伦av无码中文字幕| 国产免费人成视频在线观看播放播| 成人中文乱幕日产无线码| 国产欧美日韩视频一区二区三区|