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

        ?

        裂縫儲(chǔ)層地震物理模擬研究

        2017-06-28 14:46:24王玲玲魏建新黃平狄?guī)妥?/span>秦菽苑
        石油科學(xué)通報(bào) 2017年2期
        關(guān)鍵詞:入射角方位角方位

        王玲玲,魏建新*,黃平,狄?guī)妥專剌脑?/p>

        1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室, CNPC物探重點(diǎn)實(shí)驗(yàn)室,北京 102249 2 中國(guó)石油西南油氣田公司勘探開(kāi)發(fā)研究院,成都 610041

        裂縫儲(chǔ)層地震物理模擬研究

        王玲玲1,魏建新1*,黃平2,狄?guī)妥?,秦菽苑1

        1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室, CNPC物探重點(diǎn)實(shí)驗(yàn)室,北京 102249 2 中國(guó)石油西南油氣田公司勘探開(kāi)發(fā)研究院,成都 610041

        根據(jù)中國(guó)西南部四川盆地中部龍崗地區(qū)的地質(zhì)構(gòu)造,制作出一個(gè)小裂縫尺度、多種裂縫參數(shù)、成比例的裂縫儲(chǔ)層地震物理模型。采用薄片模擬裂縫的方法制作可變參數(shù)的小尺寸裂縫帶,模擬一定區(qū)域內(nèi)所有裂縫的等效特征。該模型是對(duì)裂縫物理模型技術(shù)的一個(gè)新嘗試。通過(guò)進(jìn)行數(shù)據(jù)采集和處理,得到分方位數(shù)據(jù)。然后分析4個(gè)裂縫密度變化的裂縫帶其中心位置處振幅隨入射角的變化情況,最后采用振幅和頻率衰減梯度兩種屬性的方位各向異性預(yù)測(cè)裂縫發(fā)育區(qū)。結(jié)果顯示,裂縫帶中心位置處振幅具有方位各向異性的特征,方位振幅的擬合結(jié)果表明預(yù)測(cè)的裂縫方位是90°,預(yù)測(cè)的4個(gè)裂縫密度值是增加的趨勢(shì),都與模型設(shè)計(jì)相符。和振幅屬性相比,頻率衰減梯度屬性的方位各向異性對(duì)小尺度裂縫更為敏感。

        裂縫;物理模型;各向異性;裂縫預(yù)測(cè)

        0 引言

        中國(guó)油氣資源豐富,但地質(zhì)條件十分復(fù)雜,資源開(kāi)發(fā)難度大且成本高,長(zhǎng)期處于油氣資源短缺的狀態(tài),石油和天然氣消費(fèi)量對(duì)外依存度逐年增加,因此,在引進(jìn)國(guó)外油氣資源的同時(shí)要加強(qiáng)國(guó)內(nèi)油氣資源研究[1]。國(guó)內(nèi)油氣資源中非常規(guī)油氣資源占的比率高且種類多,其中,致密油儲(chǔ)層是一種主要類型,具有特低滲透率、特低孔隙度、儲(chǔ)層物性差、非均質(zhì)性強(qiáng)等特征。裂縫是這類油氣儲(chǔ)層的主要存儲(chǔ)空間和滲流通道。研究裂縫發(fā)育程度和裂縫參數(shù)特征對(duì)非常規(guī)油氣的甜點(diǎn)預(yù)測(cè)起著十分重要的作用。

        地震物理模擬作為一種地震勘探的正演手段,在裂縫儲(chǔ)層識(shí)別和定量預(yù)測(cè)中起到重要作用。很多學(xué)者利用地震物理模擬技術(shù)進(jìn)行裂縫研究。模型從單一層的HTI(Traverse Isotropy with Horizontal axis)介質(zhì)模型[2-3]和包含一個(gè)HTI介質(zhì)層的多層模型[4-8]、逐步發(fā)展到包含HTI裂縫塊的模型[9],可以預(yù)見(jiàn)模擬的裂縫區(qū)將會(huì)越來(lái)越小。由于裂縫形態(tài)和尺度差異較大、裂縫填充物屬性多變、裂縫空間結(jié)構(gòu)和所處地質(zhì)構(gòu)造復(fù)雜多樣等因素,隨著勘探要求的提高,傳統(tǒng)的包含裂縫塊的地震物理模型和層狀垂直裂縫模型(HTI介質(zhì))已不能滿足要求,探索和提高地震物理模型制作和模擬技術(shù)勢(shì)在必行。

        現(xiàn)今,各種裂縫預(yù)測(cè)方法層出不窮,方法之間存在著差異和優(yōu)劣。喻岳鈺等利用瞬時(shí)頻域衰減屬性的方位各向異性特征,較為成功的預(yù)測(cè)了碳酸鹽巖裂縫[10]。姜傳金等利用縱波阻抗、頻率衰減的方位各向異性特征,準(zhǔn)確預(yù)測(cè)了火山巖的裂縫發(fā)育情況[11]。王洪求等利用振幅屬性、旅行時(shí)、旅行時(shí)差、AVO梯度等屬性的方位各向異性預(yù)測(cè)裂縫[12]。劉依謀等介紹了旅行時(shí)、速度、振幅、衰減、頻率和相位等屬性的方位各向異性原理和進(jìn)展,給出了未來(lái)發(fā)展趨勢(shì),即綜合利用地震波能量衰減、頻率、相位方位各向異性提高裂縫預(yù)測(cè)精度[13]。由于裂縫塊和裂縫層的尺度都較大,各種方法預(yù)測(cè)結(jié)果的差異并不明顯,使得采用簡(jiǎn)單的物理模型優(yōu)選預(yù)測(cè)方法存在很大的局限性。

        在三維地震物理模型制作技術(shù)[14]的基礎(chǔ)上,本文探討裂縫儲(chǔ)層模型制作技術(shù),根據(jù)實(shí)際儲(chǔ)層特征,構(gòu)建了一個(gè)小裂縫帶尺度的裂縫儲(chǔ)層地震物理模型,并對(duì)振幅隨偏移距和入射角變化進(jìn)行了分析,然后采用振幅和頻率衰減梯度兩個(gè)屬性的方位各向異性進(jìn)行裂縫預(yù)測(cè)。

        1 裂縫儲(chǔ)層地震物理模型構(gòu)建

        1.1 模型制作難題探討

        裂縫儲(chǔ)層地震物理模型是地震物理模型技術(shù)中一個(gè)新的發(fā)展方向,將儲(chǔ)層構(gòu)造特征與裂縫參數(shù)特征相結(jié)合,在制作過(guò)程中遇到了一些問(wèn)題,此處探討解決方法。

        ① 空間幾何相似比問(wèn)題。一方面,單條裂縫尺度非常小,裂縫延伸長(zhǎng)度和張開(kāi)度最小可達(dá)到毫米級(jí)別,按照幾何相似比原理[15],實(shí)驗(yàn)室中無(wú)法制作出如此小的裂縫,使得模型制作過(guò)程中空間尺度的相似性受到挑戰(zhàn)。此外,裂縫性油氣儲(chǔ)層預(yù)測(cè)中,研究單條裂縫的意義并不大,并且儲(chǔ)層裂縫常常是成群、成帶分布,因此常采用等效的方式制作裂縫帶,模擬某一區(qū)域所有裂縫的總體特征。另一方面,本文中裂縫目的層的埋深一般在三千米左右,而目的層的裂縫大小分布范圍從幾毫米到幾百米,構(gòu)建模型時(shí),較難找到合適的尺度比例因子兼顧這兩種尺寸,對(duì)模型的尺度比例有很大的限制性。綜上所述,根據(jù)現(xiàn)有物理模型實(shí)驗(yàn)技術(shù)和條件,經(jīng)分析本模型選用1:10 000,實(shí)驗(yàn)室1 mm的裂縫和地層尺度相當(dāng)野外實(shí)際10 m,這樣能模擬米級(jí)以上的小尺度裂縫帶和單條裂縫。

        ② 多參數(shù)裂縫帶的制作和邊界影響問(wèn)題。將環(huán)氧樹脂和硅橡膠混合得到的基質(zhì)材料層和滲透型紙片相互疊合,制作出不同裂縫參數(shù)的裂縫帶,也稱之為片狀疊合等效裂縫巖石物理模塊。這種裂縫帶的邊界能量過(guò)強(qiáng),裂縫識(shí)別和預(yù)測(cè)過(guò)程中容易掩蓋裂縫帶的內(nèi)部信息。為了減弱裂縫帶邊界的影響,一種理論上可行的方法是混合基質(zhì)材料的速度盡量與地層速度一致,減小混合基質(zhì)材料層和地層的波阻抗之間的差異。

        1.2 地質(zhì)背景

        四川盆地是中國(guó)重要油氣資源盆地之一,中部地區(qū)構(gòu)造平緩、東高西低,自東向西逐漸傾伏。龍崗地區(qū)大安寨段含油氣構(gòu)造是本文主要研究區(qū)域。大安寨段厚80~110 m,埋深2 300 m以上,為一套黑色頁(yè)巖(生油巖)與介殼灰?guī)r(儲(chǔ)集巖)互層。介殼灰?guī)r層薄、致密、非均質(zhì)性強(qiáng),屬于特低孔、特低滲孔隙-裂縫性儲(chǔ)層[16-19]。特低孔滲儲(chǔ)層很不均勻,“低中有高”[20]。儲(chǔ)集空間類型主要為次生溶蝕孔、洞、裂縫及微孔隙和微裂縫。高角度(60°~90°)構(gòu)造裂縫是該區(qū)儲(chǔ)層裂縫的主要類型,裂縫分布比較規(guī)則,產(chǎn)狀穩(wěn)定,常成組出現(xiàn)。

        1.3 裂縫地震模型設(shè)計(jì)

        模型設(shè)計(jì)思路有3個(gè):

        ① 依照龍崗9井區(qū)和高淺1井區(qū)大安寨段儲(chǔ)層的地質(zhì)模式、地震特征以及裂縫發(fā)育特征設(shè)計(jì)裂縫儲(chǔ)層地震物理模型(圖1)。共有6個(gè)地層,分別為第一層、沙一層、過(guò)渡層1、大安寨層、過(guò)渡層2以及第六層。表1給出各地層參數(shù)。

        ② 實(shí)際地質(zhì)模型和簡(jiǎn)化裂縫模型組合。兩個(gè)模型除了大安寨目的層有所不同,其余地質(zhì)構(gòu)造相同。如圖2所示,實(shí)際地質(zhì)模型在整個(gè)大模型的南半段,目的層模擬野外斷層、薄互層以及尖滅等,有9組裂縫帶,7個(gè)交叉裂縫,7組群縫,還有一些隨機(jī)分布的裂縫帶;簡(jiǎn)化裂縫模型在整個(gè)大模型的北半段,目的層是一個(gè)厚度約為190 m的厚層。

        ③ 簡(jiǎn)化裂縫模型。以模擬裂縫系統(tǒng)為主,放置單一裂縫參數(shù)變化的多組裂縫帶。將8組不同參數(shù)的裂縫帶規(guī)則分布在簡(jiǎn)化裂縫區(qū)。裂縫帶參數(shù)包括裂縫帶長(zhǎng)度、裂縫帶寬度、裂縫密度、兩個(gè)裂縫帶之間的水平距離、裂縫方位角、裂縫面傾角、階梯狀裂縫的寬度以及一些裂縫群的規(guī)模。

        圖3給出圖2中4條線的垂向剖面的示意圖,此處只有地層結(jié)構(gòu),沒(méi)有裂縫。設(shè)計(jì)的垂向示意圖經(jīng)過(guò)三維坐標(biāo)采集的坐標(biāo)數(shù)據(jù)修正。每幅圖的下半部分是目的層的垂向4倍放大圖,以便于清楚地看到薄互層,并且能夠看到導(dǎo)致薄互層層數(shù)變化的尖滅點(diǎn)(圖3(c))。

        空間上模型以1:10 000的比例尺度進(jìn)行設(shè)計(jì)制作,速度相似比為1:2,采樣間隔相似比為1:5 000,頻率相似比與采樣間隔相反(5 000:1)。表1顯示模型地層的密度和速度參數(shù),是相似比后的數(shù)據(jù)。裂縫物理模型尺寸為100 cm長(zhǎng),100 cm寬,20 cm高。模擬10 000 m × 10 000 m × 2 000 m的區(qū)域和深度范圍,模擬的實(shí)際深度從2 000 m到4 000 m。分析時(shí)均采用轉(zhuǎn)換成野外參數(shù)的數(shù)據(jù)。

        圖1 三維裂縫儲(chǔ)層地震物理模型Fig. 1 3D fractured reservoir seismic physical model

        表1 模型各地層速度和密度參數(shù)Table 1 Parameters of velocity and density in strata

        1.4 裂縫帶制作和測(cè)試

        在物理模型制作前需要制作裂縫帶。圖4是在幾何相似原理和動(dòng)力學(xué)相似原理的基礎(chǔ)上,基于等效介質(zhì)理論[21-23],制作具有不同裂縫參數(shù)的的裂縫帶。采用a、b、c三種不同配比的環(huán)氧樹脂和硅橡膠混合基質(zhì)材料制作不同參數(shù)的裂縫帶,其中a的基質(zhì)速度更接近圍巖速度。除了裂縫密度變化的裂縫帶,其余幾組的裂縫密度均為0.95 條/m。采用超聲波透射測(cè)試法[24]多次測(cè)量得到縱橫波速度數(shù)據(jù),參見(jiàn)表2(經(jīng)過(guò)相似比轉(zhuǎn)換為野外實(shí)際參數(shù))。除了簡(jiǎn)化裂縫區(qū)第6組縫面傾角變化的裂縫帶,其余都是具有水平對(duì)稱軸的橫向各向同性介質(zhì)。

        圖4為簡(jiǎn)化裂縫區(qū)前七組裂縫帶。第三組是裂縫密度變化的裂縫帶(圖4(c)),裂縫帶的尺度為180 × 180 × 50 m3,裂縫密度、縱波速度及縱波各向異性顯示在表3中,隨著裂縫密度增加,速度各向異性增大。第五組裂縫方位角變化的裂縫帶(圖4(e))和第六組縫面傾角變化的裂縫帶(圖4(f))尺度都是130 × 130 × 60 m3。定義裂縫走向與圖2中X方向的夾角為裂縫方位,Z方向是裂縫高度方向。第五組方位角從0°增加到90°,間隔為15°。第六組縫面傾角分別為5°、15°、30°、45°、60°、75°和 90°。

        圖5給出復(fù)雜裂縫區(qū)的裂縫帶的分布情況。大多數(shù)裂縫是高角度裂縫,少部分是中角度裂縫,幾乎沒(méi)有低角度裂縫。9組裂縫帶中最邊上靠近模型邊緣處的裂縫是縫面傾角最小的裂縫帶,約43°。9組裂縫帶中裂縫帶的寬度從6 m到14 m不等,裂縫帶長(zhǎng)度都是200 m,高度都是50 m。大多數(shù)裂縫帶與斷層之間存在斜交的關(guān)系;交叉裂縫的形態(tài)很多,如十字結(jié)構(gòu)、井字結(jié)構(gòu)等,裂縫的寬度約為1.5 m,長(zhǎng)度從200 m到400 m不等;群縫的最大寬度是378.2 m,最小是178.8 m。

        圖2 大安寨層頂界面水平裂縫分布圖Fig. 2 Fracture distribution on the top of Da’anzhai layer

        1.4 裂縫物理模型的制作地震物理模型采用逐層澆筑的方式進(jìn)行的。為了便于放置

        裂縫,澆筑流程與地層沉積次序相反。由第一層開(kāi)始澆筑,然后逐漸向下,順序是第一層、沙一層、過(guò)渡層1、大安寨層(簡(jiǎn)化裂縫區(qū)和復(fù)雜裂縫區(qū))、過(guò)渡層2以及第六層。圖6顯示物理模型的澆筑流程。和6個(gè)地層相對(duì)應(yīng),總共有6大步驟,照片中是每一層的底界面。以步驟一中圖6a為例,它不僅是第一層的底界面還是沙一層的頂界面。當(dāng)澆筑到目的層時(shí),先澆筑復(fù)雜裂縫區(qū),再澆筑簡(jiǎn)化裂縫區(qū)。在復(fù)雜裂縫區(qū),裂縫垂直粘合到大安寨頂界面上,然后一層層的澆筑薄互層。第三個(gè)薄互層(圖 6f)是一個(gè)尖滅層,嚴(yán)格按照設(shè)計(jì)要求確定尖滅線的位置。為了獲得更為準(zhǔn)確的地層參數(shù),每一層的澆筑過(guò)程中預(yù)留下一個(gè)地層樣品,對(duì)各地層樣品進(jìn)行測(cè)試得到地層速度信息,然后對(duì)設(shè)計(jì)參數(shù)進(jìn)行修正。

        圖3 圖2 中四條線的垂向剖面圖Fig. 3 Vertical profiles of four lines in figure 2

        表2 裂縫密度為0.95條/m時(shí),縱橫波速度和各向異性Table 2 Velocities and anisotropy of the fractured zones when the fracture density is 0.95 cracks per meter

        2 地震數(shù)據(jù)數(shù)據(jù)采集和處理

        采用超聲脈沖反射波法進(jìn)行寬方位地震數(shù)據(jù)采集。為了減弱面波、提高采集速度,模型數(shù)據(jù)的采集是在水箱中進(jìn)行。根據(jù)目的層的深度設(shè)計(jì)水層的厚度為200 mm,用來(lái)模擬一個(gè)2 000 m的低速層環(huán)境。采集時(shí)使用的柱塞狀壓電換能器的發(fā)射端直徑和接收端直徑分別為9 mm和14 mm(轉(zhuǎn)換成野外參數(shù)為90 m和140 m),頻率為0.23 MHz,相似比后為46 Hz。目的層波長(zhǎng)為120.7 m,換能器輻射的近場(chǎng)范圍大致由換能器輻射半徑的平方除以波長(zhǎng)來(lái)計(jì)算[25],遠(yuǎn)場(chǎng)區(qū)域波場(chǎng)穩(wěn)定,振幅變化規(guī)律明顯(逐漸減小),利于檢測(cè)裂縫。本實(shí)驗(yàn)在縱波情況下,按照目的層波長(zhǎng)計(jì)算激發(fā)換能器的近場(chǎng)距離為45 m×45 m /120.7 m=16.78 m,遠(yuǎn)遠(yuǎn)小于目的層的深度,故本實(shí)驗(yàn)測(cè)試到的是遠(yuǎn)場(chǎng)數(shù)據(jù)。

        設(shè)計(jì)的寬方位(近全方位)采集觀測(cè)系統(tǒng)參數(shù)如表4,觀測(cè)系統(tǒng)橫縱比是0.84,采樣間隔為1 ms,滿覆蓋次數(shù)為224次(圖7a),記錄長(zhǎng)度為4 s。圖7b給出采集數(shù)據(jù)的偏移距和方位角的交匯圖。偏移距大于1 500 m并小于3 000 m時(shí),數(shù)據(jù)的方位角信息比較全。數(shù)據(jù)偏移距小于1 500 m時(shí),方位角80°~100°的數(shù)據(jù)缺失,偏移距大于3 000 m數(shù)據(jù)方位角缺失更嚴(yán)重。

        模型的地震數(shù)據(jù)有比較豐富的地震波場(chǎng),從淺層到深層都有較高的信噪比,各層反射界面明顯,但是波場(chǎng)記錄中也伴隨著不同種類的干擾波。對(duì)原始數(shù)據(jù)加載觀測(cè)系統(tǒng)和道編輯后進(jìn)行一系列的疊前處理。處理過(guò)程包括觀測(cè)系統(tǒng)加載、真振幅恢復(fù)、高通濾波、衰減諧振、預(yù)測(cè)反褶積、速度分析、拉東濾波、道集分選、分方位角度疊加、疊后和疊前時(shí)間偏移。振幅補(bǔ)償處理在利用振幅信息研究裂縫儲(chǔ)層時(shí)非常重要,真振幅恢復(fù)的目的是盡量對(duì)地震波能量的衰減和畸變進(jìn)行補(bǔ)償和校正。圖8是某一單炮記錄真振幅恢復(fù)前后的對(duì)比。單炮記錄的信噪比較高,各地層反射清晰可見(jiàn)。通過(guò)對(duì)比發(fā)現(xiàn),深層(如第六層頂界面)和遠(yuǎn)偏移距(橢圓框中所示)的振幅得到補(bǔ)償和恢復(fù)。在模型內(nèi)部固-固界面上存在著多次轉(zhuǎn)換方式的干擾波,通過(guò)Radon 濾波將這種干擾波從波場(chǎng)中去除。

        3 振幅隨入射角和方位角變化

        3.1 換能器指向性能量補(bǔ)償

        每一個(gè)換能器都有特定的指向性特征,物理模型數(shù)據(jù)分析會(huì)涉及到振幅變化,尤其是在振幅隨著角度的變化(AVO)分析時(shí)需要進(jìn)行換能器指向性能量補(bǔ)償。換能器在應(yīng)用于地震物理模型數(shù)據(jù)采集之前,都需要根據(jù)實(shí)驗(yàn)確定其指向性特征,并根據(jù)實(shí)驗(yàn)得到的結(jié)果在數(shù)據(jù)處理前進(jìn)行能量補(bǔ)償。圖9a為換能器指向性示意圖。換能器的振幅隨入射角增大急劇減小,在入射角度大于60°外幾乎沒(méi)有能量。這種不同角度上振幅的差異,是指向性能量補(bǔ)償?shù)幕A(chǔ)。根據(jù)換能器指向性得到的與入射角有關(guān)的校正系數(shù):

        cc為校正系數(shù),θ是入射角,相關(guān)系數(shù)R2=0.995。

        圖9b和圖9c顯示的裂縫密度變化裂縫帶a1(CDP位置Inline448,Crossline964)的振幅隨著入射角的變化情況,b圖為換能器指向性能量補(bǔ)償前后道集,可見(jiàn)不同角度上的振幅得到了明顯的補(bǔ)償,尤其是大入射角的位置能量補(bǔ)償更為明顯。抽取裂縫帶第一個(gè)波形(圖9b中2線)負(fù)振幅的絕對(duì)值(圖9c),并將數(shù)據(jù)進(jìn)行四次多項(xiàng)式最小二乘擬合,可以看到開(kāi)角補(bǔ)償后振幅總體能量增強(qiáng),隨入射角的增大補(bǔ)償越多。

        圖4 不同裂縫參數(shù)照片和示意圖,裂縫參數(shù)分別是裂縫帶長(zhǎng)度(a)、裂縫帶寬度(b)、裂縫密度(c)、裂縫帶間距(d)、裂縫方位(e)、縫面傾角(f) 和階梯狀裂縫帶Fig. 4 Photos and sketchs of fractured zones with different fracture parameters. The fracture parameters are lengths of fractured zones (a), widths of fractured zones (b), fracture densities (c), spacings between two fractured zones (d), azimuth angles (e), fracture surface dip angles (f) and widths of terraced fractured zones (g)

        表3 不同裂縫密度裂縫帶的參數(shù)Table 3 Parameters of fractured zones with varied fracture densities

        圖5 復(fù)雜裂縫區(qū)裂縫帶(a)、交叉裂縫(b)和群縫(c)Fig. 5 Fracture zones (a), intersected fractures (b) and fracture swarms (c) in the complex area

        3.2 六個(gè)方位上振幅隨入射角和方位變化

        方位P波裂縫檢測(cè)方法是對(duì)某一給定位置點(diǎn)的一系列不同炮檢方位角的道集進(jìn)行分析,查找不同方位角上的振幅變化。以第三組裂縫密度變化的裂縫帶為例,裂縫帶編號(hào)增加,裂縫密度是增大的,分別是0.319、0.631、1.027和 1.331。圖10a和10b給出4個(gè)裂縫帶中心位置處CDPs(Common Depth Points)面元的偏移距和方位角分布,是分析中必不可少的環(huán)節(jié),為接下來(lái)的整體數(shù)據(jù)的裂縫預(yù)測(cè)提供參數(shù)設(shè)定參考,如提供數(shù)據(jù)選擇范圍、最佳面元設(shè)置。4個(gè)CDP位置處有相同的偏移距分布和相似的方位角分布,保證相同的分方位條件下、4個(gè)裂縫帶在同一方位的能量均衡。

        正北方向(X)為方位角0°方向,地震資料按方位角不同進(jìn)行分方位處理, 0°~30°和180°~210°歸為一個(gè)方位,依次類推30°為一個(gè)方位角范圍,共得到6個(gè)方位。中心角度分別為15°、45°、75°、105°、135°和165°。圖11是裂縫帶中心位置CDPs處的15°方位的振幅隨入射角變化道集,入射角從0°增加到42°,間隔為3°。道集上裂縫帶顯示出多個(gè)“波峰-波谷”的特征(疊后剖面上顯示為“串珠”狀特征),裂縫帶a1位于1 950~2 000 ms處。

        圖6 裂縫物理模型澆筑過(guò)程Fig. 6 Construction procedure of the fractured physical model

        圖7 (a)工區(qū)覆蓋次數(shù)和(b)偏移距-方位角交匯圖Fig. 7 Fold (a) and Offset-azimuth cross plot (b)

        圖8 某一單炮記錄真振幅恢復(fù)前后對(duì)比Fig. 8 A single shot records before and after true amplitude recovery

        表4 三維寬方位采集觀測(cè)系統(tǒng)參數(shù)Table 4 3D wide-azimuth geometry parameters

        圖9 換能器指向性示意圖(a)、換能器開(kāi)角校正前后入射角道集(b)和換能器開(kāi)角校正前后裂縫帶頂界面(圖b中2線)振幅隨入射角變化曲線(c)Fig. 9 The schematic diagram of the transducer directionality (a), the gathers of different incident angles before and after the opening angle correction of transducers (b), the top boundary of fractured zone (line 2 in Figure b) before and after the opening angle correction of transducers (c)

        提取紅色橢圓位置的波谷的最大振幅隨入射角變化曲線(圖12),由于不同方位上數(shù)據(jù)的不均勻性,有效的入射角道集范圍不同。總的來(lái)說(shuō),各個(gè)方位數(shù)據(jù)的振幅隨著入射角的增大而減小。同一入射角,不同方位上振幅也存在差異,中心角度15°和165°振幅值相對(duì)其他幾個(gè)方位角更小,中心角度75°和105°的振幅值最大,中心角度45°和135°振幅處于中間。說(shuō)明振幅具有方位各向異性特征。進(jìn)一步將入射角18°的振幅以極坐標(biāo)的方式顯示出來(lái)并按照對(duì)稱性將180°~360°數(shù)據(jù)補(bǔ)全,然后對(duì)數(shù)據(jù)進(jìn)行擬合。如圖13所示。擬合曲線類似“花生”形狀,a1~a4四個(gè)裂縫帶擬合曲線的長(zhǎng)軸長(zhǎng)度(同一方位兩個(gè)振幅之和)分別為0.161 8、0.164 0、0.129 8、0.111 9,其方向?qū)?yīng)的方位分別為94°、88°、86°和82°方位,總的來(lái)說(shuō)在90°方位附近,與模型設(shè)計(jì)的裂縫走向一致。

        圖10 裂縫密度變化的裂縫帶中心CDP位置處偏移距(a)和方位角(b)分布Fig. 10 Offset (a) and azimuth (b) of bins at CDPs in the central positions of fractured zones with varied fracture densities

        圖11 方位15°,裂縫帶位置疊前AVA道集Fig. 11 The pre-stack AVA gathers of four fractured zones at the central CDPs with azimuth 15°

        圖12 四個(gè)裂縫密度不同的裂縫帶在不同方位上振幅絕對(duì)值隨入射角變化曲線(線條為擬合值)Fig. 12 The variation curve of the absolute values of amplitude with the incident angles in different azimuths (Lines of different colors are the fi tting curve of the amplitude)

        圖13 入射角為18°時(shí),振幅隨方位角變化。擬合函數(shù)格式為amp=a1sin(b1θ+c1)+a2sin (b2θ+c2)。θ為入射角,amp為振幅值。Fig. 13 The variation of amplitude with the azimuth when the incident angle is 18°. Curves of different colors are the fi tting curve of the amplitude

        3.3 所有方位上振幅隨入射角和方位變化

        在6個(gè)方位道集分析的基礎(chǔ)上對(duì)所有方位的數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析,用余弦擬合方法[26]預(yù)測(cè)裂縫密度和方位,余弦擬合計(jì)算量小,適用于擬合多種地震屬性。本實(shí)驗(yàn)中,提取裂縫帶頂界面的振幅并將其進(jìn)行余弦擬合以獲得裂縫參數(shù)。裂縫主方位是超面元內(nèi)所有CDP的擬合結(jié)果,超面元大小是3 × 3 (37.5 m × 37.5 m),目標(biāo)CDP在9個(gè)面元的中心位置。圖14顯示裂縫帶中心位置CDP處裂縫預(yù)測(cè)方位的玫瑰圖。反射振幅余弦擬合出的長(zhǎng)軸和正北方位角0°的方向之間的夾角定義為裂縫方位(與設(shè)計(jì)的方位一致)。裂縫帶3-a1預(yù)測(cè)得到的裂縫方位是75°~90°;裂縫帶3-a2預(yù)測(cè)得到的裂縫主方位是90°~105°,此外還有少量的75°~90°方位角的裂縫;裂縫帶3-a3和3-a4的裂縫主方位是75°~90°,此外還有少量的90°~105°方位角的裂縫。因此,玫瑰圖顯示的裂縫帶的預(yù)測(cè)的方位角在90°附近,與設(shè)計(jì)的裂縫走向相符合。裂縫密度是另一個(gè)需要關(guān)注的參數(shù),此處裂縫密度是無(wú)量綱的,由于玫瑰圖中的色標(biāo)只能粗略地表示裂縫密度,并不詳細(xì),所以需要詳細(xì)分析密度值。4個(gè)裂縫帶預(yù)測(cè)得到的裂縫密度分別為0.621 247、 0.729 767、0.742 916 和0.813 922,可以看出4個(gè)值的增加趨勢(shì)與設(shè)計(jì)值的趨勢(shì)相吻合。值得注意的是,設(shè)計(jì)的裂縫密度是基于等效介質(zhì)理論,垂直裂縫面方向上的線密度,而預(yù)測(cè)值是面密度。

        4 地震屬性方位各向異性裂縫預(yù)測(cè)

        振幅屬性和衰減屬性的方位變化可以檢測(cè)裂縫方向和預(yù)測(cè)裂縫密度[27-29]。將采集的數(shù)據(jù)以等方位覆蓋次數(shù)的方式劃分方位角,選取的方位角范圍分別為0°~23°,22°~46°,45°~90°,90°~135°,134°~158°,157°~180°,得到中心方位角分別為11.5°,34°,67.5°,112.5°,146°,168.5°。圖15和16是由振幅屬性和頻率衰減梯度屬性的方位各向異性預(yù)測(cè)的裂縫密度分布結(jié)果。對(duì)于簡(jiǎn)化裂縫區(qū)的裂縫帶,振幅屬性的方位各向異性追蹤的是裂縫帶的邊界,內(nèi)部信息不明顯,預(yù)測(cè)到的裂縫位置與設(shè)計(jì)(圖2)相吻合,a基質(zhì)制作的裂縫帶的邊界響應(yīng)更弱,說(shuō)明保證基質(zhì)速度和圍巖速度一致來(lái)降低邊界影響是可行的。第一組小尺寸1-a1、1-b1和1-c1(50 ×42 ×100 m3)也有所指示,1-c1由于受到了干擾,指示不明顯。第二組2-a1、2-b1和2-c1(200 ×10 ×100 m3)因裂縫帶寬度10 m過(guò)小,所以基本上看不到密度指示。針對(duì)實(shí)際裂縫區(qū),振幅屬性預(yù)測(cè)的密度結(jié)果能量團(tuán)尺寸較大,對(duì)斷層和多個(gè)裂縫的整體有指示作用。頻率衰減梯度屬性的方位各向異性對(duì)對(duì)簡(jiǎn)化裂縫區(qū)的裂縫帶給出的響應(yīng)比較雜亂。對(duì)實(shí)際裂縫區(qū)中裂縫帶、群縫、交叉裂縫都有響應(yīng),橢圓選區(qū)是群縫的響應(yīng),矩形是交叉裂縫的指示,與已知分布能夠一一對(duì)應(yīng),但也出現(xiàn)了個(gè)別無(wú)法對(duì)應(yīng)的能量團(tuán),這主要是由于薄互層和地層尖滅線等復(fù)雜的地質(zhì)構(gòu)造、這些地質(zhì)構(gòu)造與裂縫相互作用等因素產(chǎn)生的干擾能量引起的。總的來(lái)說(shuō),頻率衰減梯度屬性對(duì)小尺度裂縫比較敏感,效果優(yōu)于振幅屬性。

        圖14 不同裂縫密度變化裂縫帶玫瑰圖Fig. 14 Rose maps of fractured zones with varied fracture densities

        5 討論

        圖15 振幅屬性的方位各向異性預(yù)測(cè)的裂縫密度Fig. 15 Fracture density detected by azimuthal anisotropy of amplitude attribute

        野外裂縫常常以組系形式出現(xiàn),微裂縫、中尺度裂縫的尺度是從毫米到米級(jí)。這在基于相似比的地震物理模擬中是很難模擬和制作的。本實(shí)驗(yàn)制作的裂縫帶是地震尺度級(jí)別,是某一區(qū)域內(nèi)的等效模型。模型研究裂縫帶尺度、裂縫密度、縫面傾角、裂縫方位角、多組裂縫帶之間的距離等參數(shù),不能研究裂縫張開(kāi)度參數(shù)和裂縫流體填充的情況。

        經(jīng)分析認(rèn)為進(jìn)一步研究方向有兩個(gè),一是天然裂縫的詳盡描述和模型制作,包括張開(kāi)度參數(shù)變化的裂縫和含流體的裂縫。另一個(gè)是不明顯邊界裂縫帶的建模,野外不明顯邊界裂縫帶居多,盡量減小邊界效應(yīng)對(duì)地震響應(yīng)和裂縫預(yù)測(cè)結(jié)果的影響。一種可行的方法是保證裂縫帶混合基質(zhì)材料的速度與圍巖一致。另一種可能的方法在不明顯邊界的裂縫帶制作方面,盡管很難實(shí)現(xiàn)。

        圖16 頻率衰減梯度屬性的方位各向異性預(yù)測(cè)的裂縫密度Fig. 16 Fracture density detected by azimuthal anisotropy of frequency attention gradient attribute

        6 結(jié)論

        這個(gè)模型是將裂縫帶和模擬實(shí)際地層的地震物理模型結(jié)合的一個(gè)初步試探。以四川中部龍崗地區(qū)大安寨段儲(chǔ)層為地質(zhì)背景,構(gòu)建裂縫地震物理模型。此模型能夠模擬多種裂縫參數(shù)的裂縫帶。

        對(duì)裂縫密度變化的4個(gè)裂縫帶中心位置的疊前道集進(jìn)行振幅隨入射角和方位角分析,得到裂縫帶處疊前道集數(shù)據(jù)具有振幅方位各向異性特征,滿足了裂縫的各向異性特征。擬合得到的裂縫方位與設(shè)計(jì)相吻合,同時(shí)裂縫密度變化情況也與設(shè)計(jì)的相符。在此基礎(chǔ)上,分析振幅屬性和頻率衰減梯度屬性的方位各向異性得到的裂縫密度分布,認(rèn)為頻率衰減梯度屬性在預(yù)測(cè)小尺度的裂縫方面優(yōu)勢(shì)明顯。

        致謝

        感謝中國(guó)石油西南油氣田分公司勘探開(kāi)發(fā)研究院在模型構(gòu)建過(guò)程中提供大量地質(zhì)構(gòu)造和地球物理方面珍貴資料和協(xié)助。感謝宋瑾芝、袁三一、張福宏等人的處理工作和有益的討論。

        [1] 賈承造, 龐雄奇, 姜福杰. 中國(guó)油氣資源研究現(xiàn)狀與發(fā)展方向[J]. 石油科學(xué)通報(bào), 2016, 01: 2-23. [JIA C Z, PANG X Q, JIANG F J. Research status and development directions of hydrocarbon resources in China[J]. Petroleum Science Bulletin, 2016, 01: 2-23.]

        [2] DEWANGAN P, TSVANKIN I, BATZLE M, et al. PS-wave moveout inversion for tilted TI media: A physical-modeling study[J]. Geophysics, 2006, 71: D135-D143.

        [3] ZHU Y P, TSVANKIN I, DEWANGAN P, et al. Physical modeling and analysis of P-wave attenuation anisotropy in transversely isotropic media[J]. Geophysics, 2007, 72: D1-D7.

        [4] WANG S X, LI X Y, QIAN Z P, et al. Physical modelling studies of 3-D P-wave seismic for fracture detection[J]. Geophysical Journal International, 2007, 168: 745-756.

        [5] VARELA I, MAULTZSCH S, CHAPMAN M, et al. Fracture density inversion from a physical geological model using azimuthal AVO with optimal basis functions[C] // 79th Annual International Meeting, SEG, ExpandedAbstracts. Houston, 2009: 2075-2079.

        [6] YIN Z H, LI X Y, DI B R, et al. Effects of offset-depth ratio on fracture detection - a physical modeling study[C] // 81st Annual International Meeting, SEG, ExpandedAbstracts. San Antonio, 2011: 2014-2018.

        [7] MAHMOUDIAN F, MARGRAVE F, WONG J, et al. Fracture orientation and intensity from AVAz inversion: A physical modeling study[C]// 83rd Annual International Meeting, SEG, ExpandedAbstracts. Houston, 2013: 483-487.

        [8] MAHMOUDIAN F, MARGRAVE G F, WONG J, et al. Azimuthal amplitude variation with offset analysis of physical modeling data acquired over an azimuthally anisotropic medium[J]. Geophysics, 2015, 80: C21-C35.

        [9] EKANEM A M, WEI J X, LI X Y, et al. P-wave attenuation anisotropy in fractured media: A seismic physical modelling study[J]. Geophysical Prospecting, 2013, 61(Supp. 1): 420-433.

        [10] 喻岳鈺, 楊長(zhǎng)春, 王彥飛, 等. 瞬時(shí)頻域衰減屬性及其在碳酸鹽巖裂縫檢測(cè)中的應(yīng)用[J]. 地球物理學(xué)進(jìn)展, 2009, 24(5): 1717-1722. [YU Y Y, YANG C C, WANG Y F, et al. P-wave azimuthal attenuation attributes in wavelet-scale domain and its application to fracture detection in carbonate[J]. Progress in geophysics (in Chinese) , 2009, 24(5): 1717-1722.]

        [11] 姜傳金, 鞠林波, 張廣穎, 等. 利用地震疊前數(shù)據(jù)預(yù)測(cè)火山巖裂縫的方法和效果分析-以松遼盆地北部徐家圍子斷陷營(yíng)城組火山巖為例[J]. 地球物理學(xué)報(bào), 2011, 54(2): 515-523. [JIANG C J, JU L B, ZHANG G Y. The method and effect analysis of volcanic fracture prediction with prestack seismic data-An example from the volcanic rocks of Yingcheng formation in Xujiaweizi fault depression, north of Songliao Basin[J]. Chinese Journal of Geophysics (in Chinese), 2011, 54(2): 515-523.]

        [12] 王洪求, 楊午陽(yáng), 謝春輝, 等. 不同地震屬性的方位各向異性分析及裂縫預(yù)測(cè)[J]. 石油地球物理勘探, 2014, 49(5): 925-931. [WANG H Q, YANG W Y, XIE C H, et al. Azimuthal anisotropy analysis of different seismic attributes and fracture prediction[J]. Oil Geophysical Prospecting, 2014, 49(5): 925-931.]

        [13] 劉依謀, 印興耀, 張三元, 等. 寬方位地震勘探技術(shù)新進(jìn)展[J]. 石油地球物理勘探, 2014, 49(3): 596-610. [LIU Y M, YIN X Y, ZHANG S Y, et al. Recent advances in wide-azimuth seismic exploration[J]. Oil Geophysical Prospecting, 2014, 49(3): 596-610.]

        [14] 魏建新, 牟永光, 狄?guī)妥? 三維地震物理模型的研究. 石油地球物理勘探, 2002, 37( 6): 556-561. [ WEI J X, MOU Y G, DI B R. Study of 3-D seismic physical model[J]. Oil Geophysical Prospecting, 2002, 37(6): 556-561.]

        [15] 牟永光. 三維復(fù)雜介質(zhì)地震物理模擬[M]. 北京: 石油工業(yè)出版社, 2003: 4-5. [MOU Y G. Seismic physical modeling for 3-D complex media[M]. Beijing: Petroleum Industry Press, 2003: 4-5.]

        [16] 蔣裕強(qiáng), 漆麟, 鄧海波, 等. 四川盆地侏羅系油氣成藏條件及勘探潛力[J]. 天然氣工業(yè), 2010, 30(3): 22-26. [JIANG Y Q, QI L, DENG H B, et al. Hydrocarbon accumulation conditions and exploration potentials of the Jurassic reservoirs in the Sichuan Basin[J]. Natural Gas Industry, 2010, 30(3): 22-26.]

        [17] 張福宏, 鄒定永, 陳玲, 等. 川中地區(qū)大安寨致密油儲(chǔ)層預(yù)測(cè)技術(shù)研究[J]. 天然氣勘探與開(kāi)發(fā), 2013, 36(4): 44-48. [ZHANG H F, ZOU Y D, CHEN L, et al. Prediction technology of Daanzhai tight-oil reservoir, central Sichuan Basin[J]. Natural Gas Exploration & Development, 2013, 36(4): 44-48.]

        [18] 廖群山, 胡華, 林建平, 等. 四川盆地川中侏羅系致密儲(chǔ)層石油勘探前景[J]. 石油與天然氣地質(zhì), 2011, 12(6): 815-838. [LIAO Q S, HU H, LIN J P, et al. Petroleum exploration prospect of the Jurassic tight reservoirs in central Sichuan Basin[J]. Oil & Gas Geology, 2011, 12(6): 815-838.]

        [19] 楊躍明, 楊家靜, 楊光, 等. 四川盆地中部地區(qū)侏羅系致密油研究新進(jìn)展[J]. 石油勘探與開(kāi)發(fā), 2016, 43(6): 873-882. [YANGY M, YANG J J, YANG G, et al. New research progress of Jurassic tight oil in central Sichuan Basin[J]. Petroleum Exploration and Development, 2016, 43(6): 873-882.]

        [20] 梁狄剛, 冉隆輝, 戴彈申, 等. 四川盆地中北部侏羅系大面積非常規(guī)石油勘探潛力的再認(rèn)識(shí)[J]. 石油學(xué)報(bào), 2011, 32(1): 8-17. [LIANG D G, RAN L H, DAI T S, et al. A re-recognition of the prospecting potential of Jurassic large-area and non-conventional oils in the central-northern Sichuan Basin[J]. Acta Petrolei Sinica, 2011, 32(1): 8-17.]

        [21] HUDSON J A. Wave speeds and attenuation of elastic waves in material containing cracks[J]. Geophysical Journal of the Royal Astronomical Society, 1981, 64: 133-150.

        [22] THOMSEN L. Weak elastic anisotropy[J]. Geophysics, 1986, 51: 1954-1966.

        [23] RUGER A. Reflection coefficient and azimuthal AVO analysis in anisotropic media [Ph. D. thesis]. Colorado: Colorado School of Mines, 1996.

        [24] 魏建新, 狄?guī)妥? 地震物理模型超聲測(cè)試技術(shù)中的幾個(gè)問(wèn)題[J]. 勘探地球物理進(jìn)展, 2003, 26(4): 294-300. [WEI J X, DI B R. Several problems of the ultrasonic measurement techniques of seismic physical model[J]. Progress in Exploration Geophysics, 2003, 26(4): 294-300.]

        [25] NDT Resource Center. Radiated Fields of Ultrasonic Transducers [EB/OL]. https: //www. nde-ed. org/EducationResources/Community-College/Ultrasonics/EquipmentTrans/radiatedf i elds. htm.

        [26] XIE C H, YONG X S, YANG W Y, et al. The comparison between the Cosine fi tting and AVAZ in fracture detection[C]// Beijing 2014 International Geophysical Conference & Exposition. Beijing, 2014: 614-616.

        [27] CARCIONE J M. A model for seismic velocity and attenuation in petroleum source rocks[J]. Geophysics, 2000, 65: 1080-1092.

        [28] CHICHININA T I, OBOLENTSEVA I R, RONQUILLO-JARILLO G, et al. Attenuation anisotropy of P- and S-waves: Theory and laboratory experiment[J]. Journal of Seismic Exploration, 2007, 16: 235-264.

        [29] LI F, LYU B, QI J, et al. Seismic attenuation anisotropy and corresponding frequency versus azimuth (FVAz) attribute[C]// 86th Annual International Meeting, SEG, ExpandedAbstracts. Dallas, 2016, 352-356.

        Study of seismic physical modelling of fractured reservoirs

        WANG Lingling1, WEI Jianxin1, HUANG Ping2, DI Bangrang1, QIN Shuyuan1
        1 State Key Laboratory of Petroleum Resources and Prospecting, CNPC Key Laboratory of Geophysical Exploration, China University of Petroleum-Beijing, Beijing 102249, China
        2 Exploration and Development Research Institute, PetroChina Southwest Oil & Gas fi eld Company, Chengdu 610041, China

        A small fracture-scale, multi-fractured parameter proportional model of a fractured reservoir has been constructed in the laboratory based on the geological background of Longgang in the central Szechwan Basin, southeast China. Small-scale fractured zones with variable fracture parameters are constructed to simulate the equivalent characteristics of all fractures in certain regions and paper sheets are used to simulate fractures. The model is a trial of the fractured physical model technique. Then the data acquisition and processing are carried out to obtain azimuthal data. Finally, we analyse amplitude at central CDPs of four fractured zones with different fracture densities versus incident angles and the azimuthal anisotropy of amplitude and frequency attenuation gradient attributes are used to detect the development of fractures. The results show that the amplitude of the fractured zones has characteristics of azimuthal anisotropy and the fi tting results of the azimuthal amplitude show the predicted orientations are around 90o and four fracture densities have an increasing trend, in accordance with the model design. The azimuthal anisotropy of frequency attenuation gradient attributes is more sensitive to small fractures, comparing with amplitude.

        fractures; physical model; anisotropy; fracture detection

        10.3969/j.issn.2096-1693.2017.02.020

        (編輯 付娟娟)

        王玲玲, 魏建新, 黃平, 狄?guī)妥? 秦菽苑. 裂縫儲(chǔ)層地震物理模擬研究. 石油科學(xué)通報(bào), 2017, 02: 210-227

        WANG Lingling, WEI Jianxin, HUANG Ping, DI Bangrang, QIN Shuyuan. Study of seismic physical modelling of fractured reservoirs. Petroleum Science Bulletin, 2017, 02: 210-227.doi:10.3969/j.issn.2096-1693.2017.02.020

        *通信作者, weijx@cup.edu.cn

        2017-5-28

        中國(guó)國(guó)家自然科學(xué)基金(41474112);國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展規(guī)劃(973)計(jì)劃項(xiàng)目(2013CB228601);中石油重大科技專項(xiàng)(XNS14JS2013-042)資助

        猜你喜歡
        入射角方位角方位
        一般三棱鏡偏向角與入射角的關(guān)系
        認(rèn)方位
        幼兒園(2021年12期)2021-11-06 05:10:20
        探究無(wú)線電方位在無(wú)線電領(lǐng)航教學(xué)中的作用和意義
        卷宗(2021年2期)2021-03-09 07:57:24
        近地磁尾方位角流期間的場(chǎng)向電流增強(qiáng)
        預(yù)制圓柱形鎢破片斜穿甲鋼靶的破孔能力分析*
        用經(jīng)典定理證明各向異性巖石界面異常入射角的存在
        借助方位法的拆字
        說(shuō)方位
        幼兒100(2016年28期)2016-02-28 21:26:17
        基于TMS320C6678的SAR方位向預(yù)濾波器的并行實(shí)現(xiàn)
        向量?jī)?nèi)外積在直線坐標(biāo)方位角反算中的應(yīng)用研究
        河南科技(2015年18期)2015-11-25 08:50:14
        欧美色图50p| 中文字幕日本人妻久久久免费| 国产午夜精品一区二区三区软件| 精品视频入口| 手机av在线观看视频| 亚洲一区二区三区高清在线| 免费看黄a级毛片| 91福利视频免费| 黑人一区二区三区高清视频| av素人中文字幕在线观看| 激情第一区仑乱| 在线观看亚洲AV日韩A∨| 国产av熟女一区二区三区蜜臀| 久久亚洲道色综合久久| 欧美怡红院免费全部视频| 欧美伊人亚洲伊人色综| 最全精品自拍视频在线| 2021国产精品视频网站| 一本一道av无码中文字幕 | 亚洲中文字幕久在线| 久热这里只有精品99国产| 在线免费观看亚洲毛片| 日本av在线一区二区| 少妇无码一区二区三区免费| 亚洲国产欧美另类va在线观看| 亚洲av资源网站手机在线| 台湾佬中文网站| 国产主播一区二区三区在线观看 | 国产免费av片在线观看播放| 一区二区在线观看视频亚洲| 免费成人电影在线观看 | 亚洲熟妇网| 亚洲一区二区三区在线高清中文| 强开少妇嫩苞又嫩又紧九色| 国产精品亚洲五月天高清| 在线视频一区二区观看| 亚洲av无码专区国产乱码4se| 国产精品毛片无码| 宅男久久精品国产亚洲av麻豆 | 免费无码肉片在线观看| 中文字幕成人精品久久不卡|