• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    近似解析離散化方法的粘彈聲波方程數(shù)值模擬及波場(chǎng)特征分析

    2017-06-29 02:17:40段焱文桂志先
    石油物探 2017年3期
    關(guān)鍵詞:波場(chǎng)聲波振幅

    汪 勇,段焱文,王 婷,桂志先,高 剛

    (油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(長(zhǎng)江大學(xué)),湖北武漢430100)

    近似解析離散化方法的粘彈聲波方程數(shù)值模擬及波場(chǎng)特征分析

    汪 勇,段焱文,王 婷,桂志先,高 剛

    (油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(長(zhǎng)江大學(xué)),湖北武漢430100)

    地層的粘彈性對(duì)地震波產(chǎn)生的吸收和衰減規(guī)律非常復(fù)雜,所以研究地震波在粘彈性介質(zhì)中的傳播規(guī)律對(duì)地震勘探有著重要的意義。近似解析離散化方法是近年來(lái)出現(xiàn)的一種有限差分?jǐn)?shù)值模擬方法,籍此探討了幾種近似解析離散化方法的原理和過(guò)程。通過(guò)與一維粘彈聲波方程解析解的對(duì)比,認(rèn)為改進(jìn)的近似解析離散化方法(INAD)的精度、穩(wěn)定性和計(jì)算效率更具優(yōu)勢(shì),并將其首次應(yīng)用于二維粘彈聲波方程數(shù)值模擬。通過(guò)數(shù)值模擬,重點(diǎn)分析了聲波振幅、頻率、頻帶寬度隨品質(zhì)因子及傳播距離的衰減規(guī)律,粘彈介質(zhì)中的地震波不僅因波前擴(kuò)散導(dǎo)致振幅衰減,還因吸收作用導(dǎo)致振幅和頻率同時(shí)衰減,利用其衰減特征可以指導(dǎo)含油氣儲(chǔ)層的預(yù)測(cè)。

    粘彈介質(zhì);近似解析離散方法;數(shù)值模擬;品質(zhì)因子;吸收衰減

    Keywords:viscoelastic media,nearly-analytic discrete method,numerical simulation,Qvalue,absorption and attenuation

    地震波場(chǎng)數(shù)值模擬是由已知的巖石結(jié)構(gòu)和物理參數(shù)建立地震地質(zhì)模型,依據(jù)不同的算法模擬地震波在地下介質(zhì)中的傳播過(guò)程,計(jì)算在地面或地下各觀測(cè)點(diǎn)的地震記錄的一種方法,是勘探地震學(xué)和地震波理論的重要基礎(chǔ)。目前,常用的數(shù)值模擬方法主要有射線追蹤法[1]和波動(dòng)方程法,其中,波動(dòng)方程法有偽譜法[23]、有限元法[4]、邊界元法[5]、譜元法[6]和有限差分法。自20世紀(jì)60年代ALTERMAN等[7]首先將有限差分方法用于地震波場(chǎng)的數(shù)值模擬以來(lái),該方法已經(jīng)成為地震波場(chǎng)數(shù)值模擬中應(yīng)用最為廣泛的一種方法,差分格式也由早期的低階差分發(fā)展到高階差分,由常規(guī)網(wǎng)格發(fā)展到交錯(cuò)網(wǎng)格、旋轉(zhuǎn)交錯(cuò)網(wǎng)格、可變網(wǎng)格等[8-16]。近似解析離散化(Nearly-Analytic Discrete,以下簡(jiǎn)稱(chēng)NAD)[17-18]方法是近些年發(fā)展起來(lái)的一種有限差分?jǐn)?shù)值模擬方法,基本思想是在時(shí)間上采用泰勒公式展開(kāi),在空間上利用截?cái)嗟奶├展綐?gòu)造高階插值函數(shù)來(lái)逼近空間偏導(dǎo)數(shù),與其它有限差分法的根本區(qū)別是該方法利用空間節(jié)點(diǎn)的位移和梯度來(lái)共同逼近變量的空間偏導(dǎo)數(shù),從而提高計(jì)算精度,在粗網(wǎng)格下可有效地壓制數(shù)值頻散,提高計(jì)算效率。YANG等[19-22]、盧明[23]、宋國(guó)杰[24-25]和TONG等[26]對(duì)該方法進(jìn)行了研究和改進(jìn),使得該方法逐步完善。

    粘彈性介質(zhì)是介于流相粘滯性介質(zhì)和固相完全彈性介質(zhì)之間的一種介質(zhì)。開(kāi)爾文粘彈性介質(zhì)模型是地震學(xué)中一種重要的模型,側(cè)重于固相的特征,對(duì)它的研究和數(shù)值模擬受到廣大學(xué)者的重視。CARCIONE等[27-28]對(duì)粘彈地層中傳播的聲波進(jìn)行了數(shù)值模擬。ROBERTSSON等[29]利用速度-應(yīng)力方程組,采用交錯(cuò)網(wǎng)格有限差分法對(duì)二維粘彈性介質(zhì)中的彈性波進(jìn)行了數(shù)值模擬研究,并且將該算法推廣到三維介質(zhì)中。杜啟振[30]利用偽譜法對(duì)各向異性粘彈介質(zhì)中的波場(chǎng)進(jìn)行了數(shù)值模擬。奚先等[31]利用波動(dòng)方程交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬方法,對(duì)地震波在二維粘彈性隨機(jī)介質(zhì)中的傳播進(jìn)行了模擬并研究了其波場(chǎng)特征。宋常瑜等[32]采用交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬方法進(jìn)行了井間地震粘彈性波場(chǎng)數(shù)值模擬。唐啟軍等[33]將Von Karman型的隨機(jī)各向同性背景引入粘彈性單斜各向異性波動(dòng)方程,并應(yīng)用交錯(cuò)網(wǎng)格技術(shù)進(jìn)行模擬。

    本文以一維開(kāi)爾文粘彈模型為例,首先討論4種近似解析離散化方法的基本思想,推導(dǎo)了差分方程的建立過(guò)程,介紹了空間導(dǎo)數(shù)的求取和方程穩(wěn)定性,然后進(jìn)行數(shù)值模擬,依據(jù)數(shù)值模擬精度和計(jì)算效率,確定了適合粘彈性聲波方程數(shù)值模擬的方法。并首次將近似解析離散化方法應(yīng)用到二維粘彈聲波方程的數(shù)值模擬中,通過(guò)對(duì)波場(chǎng)快照和地面地震記錄的分析,討論了粘彈介質(zhì)中聲波的衰減規(guī)律。

    1 方法原理

    地震勘探的非完全彈性介質(zhì)的研究中,經(jīng)常使用一種稱(chēng)為粘彈性體或開(kāi)爾文佛各特(Kelvin-Voigt)體的模型。這種模型假定,介質(zhì)的應(yīng)力包括兩個(gè)部分,一部分是彈性應(yīng)力,另一部分是粘滯應(yīng)力,其中彈性應(yīng)力和應(yīng)變成正比,粘滯應(yīng)力和應(yīng)變的時(shí)間變化率成正比。也就是說(shuō),應(yīng)力與應(yīng)變不再是簡(jiǎn)單的瞬時(shí)關(guān)系,而與應(yīng)變的歷史有關(guān),其位移滿(mǎn)足的矢量波動(dòng)方程表示為:

    式中:ρ為介質(zhì)密度;U=(u,v,w)T為介質(zhì)質(zhì)點(diǎn)x,y和z方向的位移;grad表示梯度算子;div表示散度算子2為拉普拉斯算子;λ和μ為介質(zhì)的拉梅系數(shù);λ′和μ′為描述粘滯性質(zhì)的兩個(gè)參數(shù),稱(chēng)為粘滯系數(shù)。

    1.1 4種NAD方法

    近似解析離散化是一種新的數(shù)值模擬方法,下面以一維粘彈介質(zhì)聲波方程為例,討論該類(lèi)方法的基本思想及4種NAD差分方程建立過(guò)程。

    由公式(1)可以得到一維粘彈介質(zhì)聲波波動(dòng)方程:

    其中,u表示質(zhì)點(diǎn)x方向的位移,η=λ+2μ,η′=λ′+2μ′。

    設(shè)滿(mǎn)足聲波方程的解為U,有如下定義:式中:u表示質(zhì)點(diǎn)位移大??;ux表示位移u在x方向的梯度;V表示U對(duì)時(shí)間的一階偏導(dǎo)數(shù);P表示U對(duì)時(shí)間的二階偏導(dǎo)數(shù)。U,V和P之間滿(mǎn)足:

    式中:m,k分別表示U對(duì)x和t的偏導(dǎo)階數(shù)。

    時(shí)間上利用截?cái)嗟奶├展綄?duì)U和V展開(kāi)可以得到:

    式中:i表示U和V的空間網(wǎng)格點(diǎn)坐標(biāo);n和n+1表示U和V的時(shí)間網(wǎng)格坐標(biāo)。P和P對(duì)時(shí)間的一階偏導(dǎo)數(shù)可以由波動(dòng)方程(2)轉(zhuǎn)化為U和V對(duì)空間的高階偏導(dǎo),即:

    由公式(5)至公式(8)可知,只需要知道U和V及它們的空間偏導(dǎo)數(shù)在第n時(shí)間層的值,即可計(jì)算第n+1時(shí)間層的值,其中,i表示空間網(wǎng)格坐標(biāo)。在公式(7)和公式(8)中,有位移U和V的空間高階偏導(dǎo)數(shù),如根據(jù)NAD類(lèi)方法的思想,這些空間導(dǎo)數(shù)是利用待求空間網(wǎng)格點(diǎn)周?chē)?7個(gè)點(diǎn)的位移及梯度共同逼近,而不是僅僅采用位移近似,這里只列舉其中兩個(gè)u對(duì)x的高階導(dǎo)數(shù)((9)式,(10)式)進(jìn)行說(shuō)明,其余高階偏導(dǎo)數(shù)的近似函數(shù)見(jiàn)文獻(xiàn)[19]。

    從公式(9)和公式(10)可以看出,要計(jì)算u對(duì)空間x的高階導(dǎo)數(shù),不僅需要x方向左右兩個(gè)點(diǎn)的位移,還需要這兩個(gè)點(diǎn)的位移梯度ux,這正是NAD類(lèi)方法與其它方法的不同之處。由于加上了梯度值,使得近似精度更高,更適合大尺度模型下粗網(wǎng)格的數(shù)值模擬[19]。

    在公式(7)和公式(8)中,有V對(duì)空間的高階導(dǎo)數(shù),如果采用公式(11)表示的一階向后差分近似,則稱(chēng)為近似解析離散化方法(NADM)[19]。式中:m表示空間偏導(dǎo)階數(shù)。如果質(zhì)點(diǎn)速度V也采用類(lèi)似(9)式和(10)式的展開(kāi)式,則稱(chēng)為改進(jìn)的近似解析離散化(Improving Nearly-Analytic Discrete,簡(jiǎn)稱(chēng)INAD)方法[22]。如果將公式(5)改為公式(12),則稱(chēng)為精化的近似解析離散化方法(Refined Nearly-Analytic Discrete,簡(jiǎn)稱(chēng)RNAD)[24],僅僅是將公式(5)中上一時(shí)刻的速度場(chǎng)替換為當(dāng)前速度場(chǎng),用來(lái)計(jì)算新時(shí)間層的位移。

    如果計(jì)算位移場(chǎng)U時(shí)的速度項(xiàng)V由n和n+1時(shí)刻速度場(chǎng)的加權(quán)組合來(lái)代替,則稱(chēng)為加權(quán)近似解析離散化方法(Weighted Nearly Analytic Discrete,簡(jiǎn)稱(chēng)WNAD)[25],即采用公式(13)對(duì)位移場(chǎng)進(jìn)行時(shí)間層推進(jìn)。

    式中:φ為加權(quán)系數(shù),取值范圍為0~1。

    1.2 算法穩(wěn)定性

    穩(wěn)定性條件是有限差分?jǐn)?shù)值模擬中一個(gè)非常重要的問(wèn)題,如果不滿(mǎn)足穩(wěn)定性條件,則前面各層的舍入誤差將會(huì)影響到后面各層的值,如果誤差的影響越來(lái)越大,會(huì)使得有限差分的結(jié)果完全錯(cuò)誤。文獻(xiàn)[19,22,24,25]利用Fourier方法研究了以上4種方法的穩(wěn)定性,并給出了穩(wěn)定性條件。假設(shè)Δt為時(shí)間步長(zhǎng),Δx為空間網(wǎng)格間距,最大的速度值為vmax,穩(wěn)定性條件為:

    穩(wěn)定性條件方面,NADM方法最為嚴(yán)格,空間步長(zhǎng)相同時(shí),所需的時(shí)間網(wǎng)格最小。INAD方法最為寬松,所需要滿(mǎn)足穩(wěn)定性條件的時(shí)間網(wǎng)格最大,更有利于進(jìn)行長(zhǎng)采樣時(shí)間的數(shù)值模擬。

    1.3 算法精度和效率分析

    為了比較NADM,RNAD,INAD和WNAD 4種近似解析離散化方法的數(shù)值模擬精度,設(shè)置了一維粘彈介質(zhì)模型,模型長(zhǎng)度16km,聲波速度4 000m/ s。模型品質(zhì)因子,其中,ω表示縱波源中心圓頻率。品質(zhì)因子與吸收系數(shù)成反比,其值越大,表示能量耗散越小。當(dāng)Q趨近于無(wú)窮大時(shí),介質(zhì)變?yōu)橥耆珡椥越橘|(zhì),粘彈聲波方程則變?yōu)榱送耆珡椥月暡ǚ匠獭?/p>

    公式(2)表示的一維粘彈介質(zhì)聲波方程,其解析解可以表示為:式中:A為振幅;波數(shù)k和吸收系數(shù)α都與頻率ω有關(guān)。

    數(shù)值模擬中,取dt=0.8ms,dx=20m。在模型的8km處加載公式(15)表示的震源,可得圖1所示的2s時(shí)刻的波場(chǎng)快照和相對(duì)誤差曲線,相對(duì)誤差定義如下:

    式中:u(tn,xi)表示精確解析解表示數(shù)值解。采用WNAD方法時(shí),公式(13)中的加權(quán)系數(shù)φ取0.5[25]。

    圖1 一維粘彈聲波方程數(shù)值模擬波場(chǎng)快照(a)和相對(duì)誤差(b)

    從圖1可以看出:①I(mǎi)NAD和WNAD兩種方法數(shù)值模擬的結(jié)果與解析解的波剖面基本重合,模擬精度較高,而NADM和RNAD的相對(duì)誤差較大;②數(shù)值模擬結(jié)果顯示,隨著傳播距離的增加,振幅隨之發(fā)生衰減,而波長(zhǎng)隨之增加。采用同樣的方法,分析不同模型參數(shù)下的相對(duì)誤差,結(jié)果如表1所示。

    表1 采用不同計(jì)算方法和模型時(shí)位移的平均相對(duì)誤差

    從表1可以看出,4種方法的相對(duì)誤差都隨時(shí)間步長(zhǎng)或空間步長(zhǎng)的增加而增大。整體而言,INAD結(jié)果的相對(duì)誤差最小,在空間步長(zhǎng)為50m時(shí),平均相對(duì)誤差為10.31%,說(shuō)明了該方法的可靠性和實(shí)用性。

    下面從所需存儲(chǔ)空間和計(jì)算時(shí)間兩方面分析以上4種方法的計(jì)算效率。由于RNAD,WNAD和INAD 3種方法在速度場(chǎng)V的時(shí)間層推進(jìn)上沒(méi)有區(qū)別,只是在計(jì)算位移場(chǎng)U 時(shí)采用的速度場(chǎng)不同,RNAD用的是n+1時(shí)刻的速度場(chǎng),INAD用的是n時(shí)刻的速度場(chǎng),WNAD用的是加權(quán)的速度場(chǎng),所以這3種方法存儲(chǔ)數(shù)組的數(shù)量沒(méi)有差別。而NADM方法由于不需要按公式(6)進(jìn)行速度場(chǎng)時(shí)間層推進(jìn),而是用位移場(chǎng)的向后差分進(jìn)行計(jì)算,所以存儲(chǔ)數(shù)組個(gè)數(shù)減少了2個(gè)。從計(jì)算時(shí)間來(lái)看,由于INAD具有最大的庫(kù)郎數(shù),在取同樣空間步長(zhǎng)時(shí),它可采用較大的時(shí)間步長(zhǎng),所以在同樣的采樣時(shí)間內(nèi),其計(jì)算次數(shù)和時(shí)間最少。例如計(jì)算上述模型1時(shí),NADM,RNAD,WNAD和INAD可取的時(shí)間步長(zhǎng)分別為0.8,1.2,1.7,2.2ms,它們的時(shí)間推進(jìn)層數(shù)分別約為2 500,1 667,1 176和909次,INAD計(jì)算時(shí)間約是NADM,RNAD和WNAD的36%,54%和77%。

    綜合來(lái)看,對(duì)粘彈聲波方程進(jìn)行數(shù)值模擬時(shí),INAD方法具有較高的精度和計(jì)算效率,所以本文利用該方法進(jìn)行數(shù)值模擬,并進(jìn)一步分析一維粘彈介質(zhì)中聲波的傳播規(guī)律。模型長(zhǎng)度設(shè)置為4km,聲波速度4 000m/s,空間步長(zhǎng)dx=10m,時(shí)間步長(zhǎng)dt=0.5ms。數(shù)值模擬時(shí),改變品質(zhì)因子和信號(hào)頻率:①Q(mào)→∞,f=50Hz;②Q=100,f=50Hz;③Q=100,f=20Hz;④Q=10,f=50Hz。在x=2km處加載雷克子波震源信號(hào),得到2 000,2 500,3 000,3 500m處質(zhì)點(diǎn)的振動(dòng)曲線和振幅譜如圖2所示,振動(dòng)曲線的主頻及振幅見(jiàn)表2。

    圖2 一維粘彈聲波方程數(shù)值模擬振動(dòng)曲線(a,b,c,d)及其對(duì)應(yīng)的頻譜(e,f,g,h)

    一維模型中傳播的聲波為平面波,不存在波前擴(kuò)散,振幅和頻率的衰減都是由介質(zhì)的粘彈性導(dǎo)致的,分析圖2和表2,結(jié)論如下。

    1)當(dāng)Q→∞時(shí),隨傳播距離和時(shí)間的增加,振幅和頻率沒(méi)有衰減,粘彈聲波方程等價(jià)于完全彈性聲波方程,同時(shí)也證明了INAD方法模擬結(jié)果的可靠性。

    2)比較②和④說(shuō)明,當(dāng)震源主頻不變時(shí),Q的大小決定了振幅和頻率的衰減速度,當(dāng)Q=10時(shí),振幅發(fā)生最大的衰減,衰減至0.004;頻率也降低得最快,降至8Hz。

    3)比較②和③可以看出,當(dāng)Q一定時(shí),衰減系數(shù)近似地與頻率的平方成正比,高頻信號(hào)衰減得快,低頻衰減得慢。②中的頻率降低了26Hz,振幅降低了90%,而③中的頻率降低了8Hz,振幅降低了74%。所以在地震勘探中,地震波隨著傳播距離的增大,高頻成分很快被衰減,深部地震信號(hào)只保留較低的頻率成分。

    4)同公式(14)表示的含義相同,振幅隨傳播距離呈指數(shù)衰減,開(kāi)始衰減速度快,隨傳播距離的增大,衰減速率降低。頻率的衰減也是如此。如模型④中的最初500m,頻率降低了36Hz,振幅降低了98%,而最后500m,頻率只降低了2Hz,振幅降低50%。

    表2 不同Q及f時(shí)的信號(hào)主頻及振幅

    2 二維模型試算

    二維粘彈介質(zhì)聲波波動(dòng)方程為:

    與一維情況相同,二維粘彈介質(zhì)聲波方程INAD方法的差分公式依然為公式(5)和公式(6),只是公式中P和有所變化,如下式所示:

    其中,

    且品質(zhì)因子Q與λ,μ,λ′和μ′有關(guān),其大小由1.3節(jié)中公式給出。

    2.1 均勻粘彈聲波方程模擬

    設(shè)置模型長(zhǎng)度和深度均為2 000m,vP=4 000m/s。為了比較,計(jì)算Q=20,50,100和趨近于∞(即完全彈性聲波模擬)4種情況。時(shí)間步長(zhǎng)Δt=0.5ms,空間網(wǎng)格Δx=Δz=10m,模型中間(x=1 000m,z=1 000m)位置激發(fā)主頻50Hz雷克子波,獲得的200ms波場(chǎng)快照如圖3所示,圖4為地表接收到的地震記錄,記錄長(zhǎng)度500ms。

    圖3 200ms時(shí)刻波場(chǎng)快照

    圖4 z=0處地震記錄

    從圖3和圖4所示的波場(chǎng)快照和地震記錄可以看出:①模擬波場(chǎng)非常清晰,沒(méi)有數(shù)值頻散;②同一時(shí)刻,粘彈聲波模擬和完全彈性(Q→∞)聲波模擬的波前位置相同,證實(shí)了方法的可靠性;③隨著Q的減小,地震波波長(zhǎng)增加,波數(shù)降低;頻率降低,同相軸變波的振幅歸一化到0~1,得到炮檢距 振幅曲線(圖5)。從圖5可以看出:①當(dāng)采用完全彈性聲波模擬時(shí)(Q→∞),與一維情況不同,振幅隨地震波的傳播發(fā)生波前擴(kuò)散,振幅與傳播距離呈反比關(guān)系;②當(dāng)采用粘彈聲波模擬時(shí),除存在波前擴(kuò)散外,還存在粘彈介質(zhì)對(duì)地震波的吸收衰減,使得衰減加??;③衰減速度由Q決定,Q越小,衰減越快。寬,子波延續(xù)時(shí)間增加;④粘彈性介質(zhì)中,地震波隨傳播距離的增加,振幅呈指數(shù)衰減,地震記錄中遠(yuǎn)道的振幅要明顯小于近道。

    為了分析地震波衰減規(guī)律,將圖4中的4個(gè)反射

    圖5 炮檢距 振幅曲線

    對(duì)圖4所示的地震記錄,提取x=1 000m位置處質(zhì)點(diǎn)的單道記錄并進(jìn)行頻譜分析,得到不同Q情況下的振動(dòng)曲線的振幅譜(圖6)。圖6顯示了不同Q時(shí)頻率衰減大小,完全彈性聲波模擬時(shí)記錄信號(hào)的主頻依然為50Hz,頻帶范圍為32~70Hz。隨著Q的減小,粘彈聲波信號(hào)的高頻成分減少得越多,Q=20時(shí)的頻率降低至14Hz,頻帶范圍減小為8~19Hz,這與實(shí)際地震記錄中的頻譜特征一致。

    圖6 振幅譜曲線(x=1 000m,z=0)

    2.2 多層水平層狀模型模擬

    利用INAD方法對(duì)4層層狀模型分別進(jìn)行粘彈和完全彈性聲波方程模擬試算,模型長(zhǎng)度和深度均為2 000m,各層巖石物理參數(shù)如表3所示。

    表3 4層水平層狀介質(zhì)模型參數(shù)

    為了比較粘彈性和完全彈性介質(zhì)中地震波場(chǎng)的差異,第2層和第3層的Q分別取50,100。時(shí)間步長(zhǎng)Δt=0.5ms,空間網(wǎng)格Δx=Δz=10m,模型地表中間(x=1 000m,z=0)位置激發(fā)主頻為40Hz的雷克子波,獲得不同時(shí)刻的波場(chǎng)快照(圖7,圖8,圖9)和地表接收的單炮記錄(圖10),記錄長(zhǎng)度600ms。在數(shù)值模擬中,采用PML吸收邊界處理,算法見(jiàn)文獻(xiàn)[33]和文獻(xiàn)[34],這里不再贅述。

    圖7至圖10中,P1,P2和P3代表來(lái)自3個(gè)反射界面的反射P波。圖7中,因?yàn)榈?層設(shè)置為完全彈性介質(zhì),所以粘彈聲波和彈性聲波模擬快照中的P1相同。而透射進(jìn)入第2層、第3層地層的P波則出現(xiàn)了差異,正是由于粘彈介質(zhì)對(duì)地震波的吸收衰減,導(dǎo)致了粘彈聲波模擬結(jié)果中的P2和P3反射波振幅明顯衰減,同相軸變寬,子波延續(xù)時(shí)間增加,波長(zhǎng)增大。

    圖7 t=300ms波場(chǎng)快照

    圖8 t=400ms波場(chǎng)快照

    圖9 t=500ms波場(chǎng)快照

    圖10 地震記錄

    圖11 x=1 500m處單道地震記錄

    圖11是從上述兩次模擬地面記錄中抽取的單道地震記錄(x=1 500m),由于直達(dá)波能量太強(qiáng),所以顯示時(shí)間從300ms開(kāi)始。波前擴(kuò)散導(dǎo)致彈性聲波記錄的振幅隨傳播時(shí)間和距離逐步減弱,而粘彈聲波記錄的振幅衰減速度增大,這與波場(chǎng)快照中的規(guī)律一致。圖11顯示了第1層反射波(350ms左右)振幅,可以看出,彈性聲波的振幅要略小于粘彈聲波的振幅,這是由于第2層變?yōu)檎硰椊橘|(zhì)后,導(dǎo)致P波的垂直反射系數(shù)增大。為了進(jìn)一步分析兩次模擬記錄在頻譜上的差別,利用廣義S變換,對(duì)這兩道地震記錄進(jìn)行時(shí)頻分析,結(jié)果如圖12和圖13所示。

    從時(shí)頻分析結(jié)果中提取3層反射波主頻和振幅,結(jié)果如表4所示。表4中數(shù)據(jù)表明:①當(dāng)采用彈性聲波模擬時(shí),子波信號(hào)的主頻和頻帶范圍基本不變,實(shí)際上這并不能準(zhǔn)確地描述地震波在地下的傳播規(guī)律。在傳播過(guò)程中發(fā)生的振幅衰減是由波前擴(kuò)散導(dǎo)致的,在2 000m厚度模型中,振幅衰減了約50%。②采用粘彈聲波模擬能更加精確地描述地震波在粘彈介質(zhì)中的傳播特征,子波主頻和頻帶寬度均隨傳播深度的增加而減小,衰減速度與介質(zhì)的品質(zhì)因子和子波頻率有關(guān)。振幅除了波前擴(kuò)散外,衰減會(huì)變得更劇烈,本模型中,振幅衰減了約98%。③雖然子波頻率降低和頻帶變窄會(huì)導(dǎo)致地震資料縱向分辨率降低,但這些頻率衰減特征正好可以指導(dǎo)含油氣儲(chǔ)層的預(yù)測(cè)工作(因?yàn)閹r石含流體后,會(huì)使其粘彈性質(zhì)更加明顯)。

    圖12 x=1 500m處單道地震記錄時(shí)頻分析結(jié)果(完全彈性聲波模擬)

    圖13 x=1 500m處單道地震記錄時(shí)頻分析結(jié)果(粘彈聲波模擬)

    表4 彈性和粘彈性聲波模擬的3層反射波主頻及振幅

    3 結(jié)束語(yǔ)

    粘彈介質(zhì)地震波場(chǎng)的數(shù)值模擬研究對(duì)地震勘探具有重要的意義,有助于人們更清楚地認(rèn)識(shí)地震波在巖石和油氣儲(chǔ)集區(qū)的傳播規(guī)律,具有較強(qiáng)的理論和現(xiàn)實(shí)意義。本文將一種新的數(shù)值模擬方法應(yīng)用到粘彈介質(zhì)聲波方程的數(shù)值模擬中,研究結(jié)果表明:①改進(jìn)的近似解析離散化(INAD)有限差分算法在粗網(wǎng)格下具有較高的精度,更有利于大尺度模型的數(shù)值模擬,為粘彈介質(zhì)的地震波場(chǎng)數(shù)值模擬提供了新的思路和方法;②采用完全彈性聲波模擬時(shí),振幅只隨波前擴(kuò)散發(fā)生衰減,而子波信號(hào)的頻率和頻帶范圍基本不變;③采用粘彈聲波模擬時(shí),振幅在波前擴(kuò)散的基礎(chǔ)上,衰減幅度增加;同時(shí),子波主頻和頻帶寬度隨傳播深度的增加而減小,衰減速度與介質(zhì)的品質(zhì)因子和子波頻率有關(guān)。

    由于近似解析離散化方法在空間上可求得位移的最高5階空間導(dǎo)數(shù),且粘彈介質(zhì)波動(dòng)方程中存在速度項(xiàng),使得本文中采用的INAD方法的差分精度為時(shí)間二階和空間四階,在下一步研究中,需要進(jìn)一步研究時(shí)間差分精度更高的近似解析離散化方法,提高粘彈介質(zhì)的模擬精度。其次,可將該類(lèi)方法在其它粘彈介質(zhì)模型的全波場(chǎng)數(shù)值模擬中進(jìn)一步推廣,擴(kuò)大該類(lèi)方法的適用范圍。

    [1] CERVENY V.Seismic rays and ray intensities in inhomogeneous anisotropic media[J].Geophysical Journal International,1972,29(1):1-33

    [2] KOSLOFF D D,BAYSAL E.Forward modeling by a Fourier method[J].Geophysics,1982,47(10):1402-1412

    [3] RESHEF M,KOSLOFF D D,EDWARDS M,et al.Three dimensional elastic modeling by the Fourier method[J].Geophysics,1988,53(9):1184-1193

    [4] DRAKE L A.Rayleigh waves at a continental boundary by the finite element method[J].Bulletin of the Seismological Society of America,1972,62(5):1259-1268

    [5] BOUCHON M,COUTANT O.Calculation of synthetic seismograms in a laterally varying medium by the boundary element-discrete wavenumber method[J].Bulletin of the Seismological Society of America,1994,84(6):1869-1881

    [6] KOMATISSCH D,BARNES C,TROMP J.Simulation of anisotropic wave propagation based upon a spectral element method[J].Geophysics,2000,65(4):1251-1260

    [7] ALTERMAN Z,KARAL F C.Propagation of seismic wave in layered media by finite difference methods[J].Bulletin of the Seismological Society of America,1968,58(1):367-398

    [8] KELLY K R,WARD R W,TREITEL S,et al.Synthetic seismograms:a finite-difference approach[J].Geophysics,1976,41(1):2-27

    [9] VIRIEUX J.P-SV wave propagation in heterogeneous media:velocity-stress finite-difference method[J].Geophysics,1986,51(4):889-901

    [10] GRAVES R W.Simulation seismic wave propagation in 3Delastic media using staggered-grid finite differences[J].Bulletin of the Seismological Society of America,1996,86(4):1091-1106

    [11] 董良國(guó),馬在田,曹景忠,等.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法[J].地球物理學(xué)報(bào),2000,43(3):411-419 DONG L G,MA Z T,CAO J Z,et al.The staggeredgrid high-order difference method of first-order elastic equation[J].Chinese Journal of Geophysics,2000,43(3):411-419

    [12] 董良國(guó),馬在田,曹景忠.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法穩(wěn)定性研究[J].地球物理學(xué)報(bào),2000,43(6):856-864 DONG L G,MA Z T,CAO J Z.The stability study of the staggered-grid high-order difference method of first-order elastic equation,[J].Chinese Journal of Geophysics,2000,43(6):856-864

    [13] 董良國(guó).彈性波數(shù)值模擬中的吸收邊界條件[J].石油地球物理勘探,1999,34(1):46-56

    DONG L G.Absorptive boundary condition in elasticwave numerical modeling[J].Oil Geophysical Prospecting,1999,34(1):46-56

    [14] 王書(shū)強(qiáng),楊頂輝,楊寬德.彈性波方程的緊致差分方法[J].清華大學(xué)學(xué)報(bào)(自然科學(xué)版),2002,42(8):1128-1131

    WANG S Q,YANG D H,YANG K D.Compact finite difference scheme for elastic equations[J].Journal of Tsinghua University(Science and Technology),2002,42(8):1128-1131

    [15] SAENGER E H,GOLD N,SHAPIRO S A.Modeling the propagation of elastic waves using a modified finitedifference grid[J].Wave Motion,2000,31(1):77-92

    [16] SAENGER E H,BOHLEN T.Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[J].Geophysics,2004,69(2):583-591

    [17] KONDOH Y,HOSAKA Y,ISHII K.Kernel optimum nearly-analytical discrimination algorithm applied to parabolic and hyperbolic equations[J].Computers &Mathematics with Applications,1994,27(3):59-90

    [18] 楊頂輝,滕吉文,張中杰.三分量地震波場(chǎng)的近似解析離散模擬技術(shù)[J].地球物理學(xué)報(bào),1996,39(增刊1):283-291

    YANG D H,TENG G W,ZHANG Z J.Nearly-analytical discretization modeling technique of 3-component seismic wave-fields[J].Chinese Journal of Geophysics,1996,39(S1):283-291

    [19] YANG D H,TENG J W,ZHANG Z J,et al.A nearlyanalytic discrete method for acoustic and elastic wave equation[J].Bulletin of the Seismological Society of A-merica,2003,93(2):882-890

    [20] YANG D H,LU M,WU R S,et al.An optimal nearlyanalytic discrete method for 2Dacoustic and elastic wave equations[J].Bulletin of the Seismological Society of America,2004,94(5):1982-1991

    [21] YANG D H,SONG G J,LU M.Optimally accurate nearly-analytic discrete scheme for wave-filed simulation in 3Danisotropic media[J].Bulletin of the Seismological Society of America,2007,97(5):1557-1569

    [22] YANG D H,SONG G J,CHEN S,et al.An improved nearly analytical discrete method:an efficient tool to simulate the seismic response of 2-D porous structures[J].Journal of Geophysics &Engineering,2007,4(1):40-52

    [23] 盧明.改進(jìn)的近似解析離散化方法及彈性波波場(chǎng)模擬[D].北京:清華大學(xué),2004

    LU M.The optimum nearly analytic discrete methodand elastic wave numerical simulations[D].Beijing:Tsinghua University,2004

    [24] 宋國(guó)杰.三維彈性波方程的改進(jìn)近似解析離散化方法及波場(chǎng)模擬[D].北京:清華大學(xué),2008

    SONG G J.The improved nearly analytic discrete method for 3Delastic equations and its numerical simulation[D].Beijing:Tsinghua University,2008

    [25] 宋國(guó)杰.基于WNAD方法的非一致網(wǎng)格算法及其彈性波場(chǎng)模擬[J].地球物理學(xué)報(bào),2010,53(8):1985-1992

    SONG G J.Non-uniform grid algorithm based on the WNAD method and elastic wave-field simulations[J].Chinese Journal of Geophysics,2010,53(8):1985-1992

    [26] TONG P,YANG D H,HUA B L.High accuracy wave simulation-Revised derivation,numerical analysis and testing of a nearly analytic integration discrete method for solving acoustic equation[J].International Journal of Solids &Structures,2011,48(1):56-70

    [27] CARCIONE J M,KOSLOFF D,KOSLOFF R.Wave propagation simulation in a linear viscoacoustic medium[J].Geophysical Journal International,1988,93(2):393-401

    [28] CARCIONE J M,KOSLOFF D,KOSLOFF R.Viscoacoustic wave propagation simulation in the Earth[J].Geophysics,1988,53(6):769-777

    [29] ROBERTSSON J O A,BLANCH J O.Viscoelastic finite-difference modeling[J].Geophysics,1994,59(9):1444-1456

    [30] 杜啟振.各向異性黏彈性介質(zhì)偽譜法波場(chǎng)模擬[J].物理學(xué)報(bào),2004,53(12):4428-4434

    DU Q Z.Wavefield forward modeling with the pseudo-spectral method in viscoelastic and azimuthally anisotropic media[J].Acta Physica Sinica,2004,53(12):4428-4434

    [31] 奚先,姚姚.二維粘彈性隨機(jī)介質(zhì)中的波場(chǎng)特征分析[J].地球物理學(xué)進(jìn)展,2004,19(3):608-615

    XI X,YAO Y.The analysis of the wave field characteristics in 2-D viscoelastic random medium[J].Progress in Geophysics,2004,19(3):608-615

    [32] 宋常瑜,裴正林.井間地震粘彈性波場(chǎng)特征的數(shù)值模擬研究[J].石油物探,2006,45(5):508-513

    SONG C Y,PEI Z L.Numerical simulation of viscoelastic wevefield exploration[J].Geophysical Prospecting for Petroleum,2006,45(5):508-513

    [33] 唐啟軍,韓立國(guó),王恩利,等.基于隨機(jī)各向異性背景的粘彈性單斜介質(zhì)二維三分量正演模擬[J].西北地震學(xué)報(bào),2009,30(1):35-39

    TANG Q J,HAN L G,WANG E L.2-D three-component seismic modeling for viscoelastic monoclinic media based on background of random isotropic media[J].Northwestern Seismological Journal,2009,30(1):35-39

    [34] 單啟銅,樂(lè)友喜.PML邊界條件下二維粘彈性介質(zhì)波場(chǎng)模擬[J].石油物探,2007,46(2):126-130

    SHAN Q T,YUE Y X.Wavefield simulation of 2-D viscoelastic medium in perfectly matched layer boundary[J].Geophysical Prospecting for Petroleum,2007,46(2):126-130

    (編輯:顧石慶)

    新書(shū)介紹

    《裂縫各向異性地震等效介質(zhì)理論》(精裝),王赟,劉媛媛,張美根著,2017年3月1日由科學(xué)出版社出版。本書(shū)系統(tǒng)介紹了地震各向異性領(lǐng)域近30年來(lái)發(fā)展相對(duì)成熟的4種裂縫誘導(dǎo)地震各向異性的等效介質(zhì)理論Hudson擾動(dòng)理論、Thomsen等徑孔隙理論、Schoenberg線性滑移理論和Chapman多尺度裂縫等效介質(zhì)理論;并通過(guò)典型裂縫介質(zhì)建模與地震波場(chǎng)的有限元數(shù)值模擬、分析,對(duì)比了4種理論的異同及裂縫介質(zhì)的地震波場(chǎng)特征,重點(diǎn)討論了地震波在多組多尺度裂縫介質(zhì)中的傳播。該書(shū)可作為固體地球物理地震專(zhuān)業(yè)本科生、研究生的教材,也可供工程技術(shù)人員參考。

    Numerical simulation and the wave field characteristics analysis of viscoelastic acoustic wave equation based on the nearly-analytic discrete method

    WANG Yong,DUAN Yanwen,WANG Ting,GUI Zhixian,GAO Gang
    (Key Laboratory of Exploration Technologies for Oil and Gas Resources,Ministry of Education(Yangtze University),Wuhan 430100,China)

    Viscoelasticity of the media has a very complicated relationship with the absorption and attenuation of seismic wave energy.Therefore,the research on the propagation law of seismic waves in the stratigraphic equivalent media carries great significance to the seismic exploration.Nearly-analytic discrete method is a finite-difference numerical simulation method developed in recent years.In this paper,the author firstly explored the principles and process of four different nearly-analytic discrete methods.In comparison with the one-dimensional analytical solutions,the improved nearly-analytic discrete method(INAD)was proved to excel in precision,stability and computational efficiency.Then INAD was applied to the numerical simulation of 2Dviscoelastic acoustic wave equation for the first time.With the simulated viscoelastic wave field,the author focused on analyzing the propagation law of the amplitude,frequency and bandwidth of the acoustic wave with the change of Qvalues and propagation distance The amplitude of the seismic wave in the viscoelastic medium reduces not only due to the wavefront diffusion,but also the absorption which attenuates the frequency at the same time.These attenuation characteristics can guide the prediction of the oil and gas reservoir.The amplitude of the seismic wave in the viscoelastic medium reduces not only due to the wavefront diffusion,but also the absorption which attenuates the frequency at the same time.These attenuation characteristics can guide the prediction of the oil and gas reservoir.

    P631

    A

    1000-1441(2017)03-0362-11

    10.3969/j.issn.1000-1441.2017.03.006

    汪勇,段焱文,王婷,等.近似解析離散化方法的粘彈聲波方程數(shù)值模擬及波場(chǎng)特征分析[J].石油物探,2017,56(3):362-372

    WANG Yong,DUAN Yanwen,WANG Ting,et al.Numerical simulation and the wave field characteristics analysis of viscoelastic acoustic wave equation based on the nearly-analytic discrete method[J].Geophysical Prospecting for Petroleum,2017,56(3):362-372

    2016-03-24;改回日期:2016-10-18。

    汪勇(1979—),男,副教授,博士,現(xiàn)主要從事地震波場(chǎng)特征研究和儲(chǔ)層預(yù)測(cè)研究。

    國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)(2013CB228605)和中國(guó)石油科技創(chuàng)新基金項(xiàng)目(2015D-5006-0301)聯(lián)合資助。

    This research is financially supported by the National Key Basic Research Project of China(973Program)(Grant No.2013CB228605)and the Science and Technology Innovation Foundation of CNPC(Grant No.2015D-5006-0301).

    猜你喜歡
    波場(chǎng)聲波振幅
    彈性波波場(chǎng)分離方法對(duì)比及其在逆時(shí)偏移成像中的應(yīng)用
    愛(ài)的聲波 將愛(ài)留在她身邊
    聲波殺手
    十大漲跌幅、換手、振幅、資金流向
    十大漲跌幅、換手、振幅、資金流向
    自適應(yīng)BPSK在井下鉆柱聲波傳輸中的應(yīng)用
    交錯(cuò)網(wǎng)格與旋轉(zhuǎn)交錯(cuò)網(wǎng)格對(duì)VTI介質(zhì)波場(chǎng)分離的影響分析
    基于Hilbert變換的全波場(chǎng)分離逆時(shí)偏移成像
    十大漲跌幅、換手、振幅、資金流向
    “聲波驅(qū)蚊”靠譜嗎
    亚洲不卡免费看| 中国国产av一级| a级毛片免费高清观看在线播放| 欧美人与性动交α欧美精品济南到 | 成人午夜精彩视频在线观看| 80岁老熟妇乱子伦牲交| 91精品一卡2卡3卡4卡| 久久久久久久久久成人| 欧美激情极品国产一区二区三区 | 久久久久国产网址| 国产色爽女视频免费观看| 国产免费又黄又爽又色| 国产精品麻豆人妻色哟哟久久| 国产欧美日韩一区二区三区在线 | 日本vs欧美在线观看视频| 国产黄色免费在线视频| 热99久久久久精品小说推荐| 国产精品99久久久久久久久| 在线播放无遮挡| 久久精品国产亚洲av涩爱| 在线观看人妻少妇| 午夜视频国产福利| 欧美激情极品国产一区二区三区 | 国产免费一级a男人的天堂| 国产黄色免费在线视频| 国产欧美日韩综合在线一区二区| 日韩制服骚丝袜av| 国国产精品蜜臀av免费| 男女国产视频网站| 91精品国产国语对白视频| av国产久精品久网站免费入址| 男女国产视频网站| 黑人猛操日本美女一级片| 黑人猛操日本美女一级片| 一级二级三级毛片免费看| 国产精品一区二区在线不卡| 日本av手机在线免费观看| 毛片一级片免费看久久久久| 久久国内精品自在自线图片| 狂野欧美白嫩少妇大欣赏| 香蕉精品网在线| 亚洲国产欧美日韩在线播放| 亚洲av不卡在线观看| 丁香六月天网| 内地一区二区视频在线| 香蕉精品网在线| 乱码一卡2卡4卡精品| 久热这里只有精品99| 国产欧美亚洲国产| 精品国产一区二区三区久久久樱花| 成年美女黄网站色视频大全免费 | videossex国产| 精品一区二区三区视频在线| 韩国av在线不卡| 18禁观看日本| 久久久a久久爽久久v久久| 极品少妇高潮喷水抽搐| 99久久精品一区二区三区| 午夜影院在线不卡| 欧美激情 高清一区二区三区| 精品人妻在线不人妻| 亚洲,欧美,日韩| 一级毛片aaaaaa免费看小| 伦精品一区二区三区| 午夜免费观看性视频| 午夜免费观看性视频| 看十八女毛片水多多多| 夫妻性生交免费视频一级片| 欧美日韩综合久久久久久| 日本黄大片高清| 欧美精品人与动牲交sv欧美| 欧美人与性动交α欧美精品济南到 | 中文字幕最新亚洲高清| 久久97久久精品| 亚洲精品日本国产第一区| 一本一本综合久久| 亚洲欧美清纯卡通| 亚洲精品456在线播放app| 色哟哟·www| 精品久久国产蜜桃| 国产免费视频播放在线视频| 一本一本综合久久| 18禁动态无遮挡网站| 日韩电影二区| 一区二区日韩欧美中文字幕 | 国产高清三级在线| kizo精华| 欧美日韩一区二区视频在线观看视频在线| 久久国产精品大桥未久av| 亚洲av欧美aⅴ国产| 国产黄色免费在线视频| 日日撸夜夜添| 国产欧美另类精品又又久久亚洲欧美| 色哟哟·www| 美女内射精品一级片tv| 最近2019中文字幕mv第一页| 熟妇人妻不卡中文字幕| 久热这里只有精品99| 伦精品一区二区三区| 国产免费福利视频在线观看| 黄片无遮挡物在线观看| av免费在线看不卡| 黄片播放在线免费| av国产精品久久久久影院| 日韩亚洲欧美综合| 中文字幕制服av| 亚洲第一av免费看| www.av在线官网国产| 国产成人精品在线电影| 婷婷色麻豆天堂久久| 亚洲伊人久久精品综合| 97在线视频观看| 美女xxoo啪啪120秒动态图| 午夜福利影视在线免费观看| 在线观看www视频免费| 国产成人精品一,二区| 99re6热这里在线精品视频| 91久久精品电影网| 亚洲怡红院男人天堂| 多毛熟女@视频| 91精品伊人久久大香线蕉| 女性被躁到高潮视频| 香蕉精品网在线| 蜜桃在线观看..| 熟女电影av网| av在线观看视频网站免费| 国产欧美另类精品又又久久亚洲欧美| 欧美少妇被猛烈插入视频| 有码 亚洲区| 欧美日韩视频高清一区二区三区二| av电影中文网址| 嘟嘟电影网在线观看| 亚洲少妇的诱惑av| 十八禁网站网址无遮挡| 午夜福利网站1000一区二区三区| 婷婷色综合大香蕉| 国产精品免费大片| 男男h啪啪无遮挡| 久久久久网色| 国产视频内射| 一个人免费看片子| 91精品三级在线观看| 成人毛片a级毛片在线播放| 桃花免费在线播放| 免费久久久久久久精品成人欧美视频 | 久久毛片免费看一区二区三区| 一级毛片黄色毛片免费观看视频| 国产成人精品无人区| av播播在线观看一区| 最近最新中文字幕免费大全7| 精品少妇久久久久久888优播| 99热6这里只有精品| 日本91视频免费播放| 久久鲁丝午夜福利片| 国产av码专区亚洲av| 国产精品 国内视频| 我的女老师完整版在线观看| 激情五月婷婷亚洲| 精品卡一卡二卡四卡免费| 国产69精品久久久久777片| 久久精品熟女亚洲av麻豆精品| 精品酒店卫生间| 2022亚洲国产成人精品| 91aial.com中文字幕在线观看| 夫妻午夜视频| 午夜福利视频精品| a级毛色黄片| 精品亚洲成a人片在线观看| 国产精品国产av在线观看| 纵有疾风起免费观看全集完整版| 午夜av观看不卡| 亚洲欧洲日产国产| 最近中文字幕高清免费大全6| 一本一本综合久久| 国产精品久久久久久久久免| 国产日韩欧美在线精品| 欧美97在线视频| 亚洲,欧美,日韩| 亚洲在久久综合| av播播在线观看一区| 国产伦精品一区二区三区视频9| 99九九线精品视频在线观看视频| 欧美xxxx性猛交bbbb| 少妇精品久久久久久久| 大香蕉久久网| 日本-黄色视频高清免费观看| 菩萨蛮人人尽说江南好唐韦庄| 黄色毛片三级朝国网站| 婷婷色av中文字幕| 精品人妻在线不人妻| 久久鲁丝午夜福利片| 91精品国产国语对白视频| 欧美成人精品欧美一级黄| 麻豆成人av视频| 免费高清在线观看视频在线观看| 精品午夜福利在线看| 亚洲精品色激情综合| 一级a做视频免费观看| 99久国产av精品国产电影| 国产亚洲精品第一综合不卡 | 国产色爽女视频免费观看| 考比视频在线观看| 22中文网久久字幕| 欧美精品国产亚洲| 国产精品嫩草影院av在线观看| 国产亚洲精品第一综合不卡 | 亚洲第一区二区三区不卡| 欧美日韩综合久久久久久| 亚洲国产av新网站| 国产成人免费无遮挡视频| 黑人高潮一二区| 少妇猛男粗大的猛烈进出视频| 一级黄片播放器| 国产乱人偷精品视频| 国产精品99久久久久久久久| 久久99精品国语久久久| 亚洲欧洲精品一区二区精品久久久 | 久久久久久久久久成人| 亚洲激情五月婷婷啪啪| 一级a做视频免费观看| 高清不卡的av网站| 91精品国产九色| 欧美人与善性xxx| 精品人妻偷拍中文字幕| 国产一区二区在线观看日韩| 18禁在线播放成人免费| 妹子高潮喷水视频| 亚洲国产成人一精品久久久| 久热这里只有精品99| 亚洲性久久影院| 久久国产精品大桥未久av| 亚洲怡红院男人天堂| 丰满少妇做爰视频| 男人爽女人下面视频在线观看| 午夜影院在线不卡| 边亲边吃奶的免费视频| 少妇丰满av| 99久久综合免费| 亚洲精品456在线播放app| 日韩欧美一区视频在线观看| 大陆偷拍与自拍| 国产成人免费观看mmmm| 精品久久久精品久久久| 一边摸一边做爽爽视频免费| 人人妻人人爽人人添夜夜欢视频| 国产免费现黄频在线看| 纵有疾风起免费观看全集完整版| 少妇丰满av| 成人无遮挡网站| 曰老女人黄片| 秋霞伦理黄片| 中国美白少妇内射xxxbb| 黄片无遮挡物在线观看| 国产黄片视频在线免费观看| 国产 精品1| 久久毛片免费看一区二区三区| 国产成人freesex在线| 亚洲精品久久久久久婷婷小说| 丰满迷人的少妇在线观看| 中文乱码字字幕精品一区二区三区| 青春草视频在线免费观看| a级毛片在线看网站| 最近最新中文字幕免费大全7| 亚洲国产色片| 亚洲av成人精品一区久久| 久久久久国产网址| 一区二区av电影网| 久久 成人 亚洲| 一级黄片播放器| 一个人看视频在线观看www免费| 午夜av观看不卡| av国产精品久久久久影院| 亚洲人成网站在线播| 看十八女毛片水多多多| 最新中文字幕久久久久| 免费观看无遮挡的男女| 乱码一卡2卡4卡精品| 国产白丝娇喘喷水9色精品| www.色视频.com| 毛片一级片免费看久久久久| 亚洲少妇的诱惑av| 高清黄色对白视频在线免费看| a级毛片免费高清观看在线播放| 曰老女人黄片| 人妻一区二区av| 中文欧美无线码| 国产免费视频播放在线视频| 新久久久久国产一级毛片| 寂寞人妻少妇视频99o| 少妇的逼好多水| 国产视频内射| 99re6热这里在线精品视频| 极品人妻少妇av视频| 精品卡一卡二卡四卡免费| 丁香六月天网| 成人毛片a级毛片在线播放| 国产亚洲午夜精品一区二区久久| 亚洲精品成人av观看孕妇| 午夜福利,免费看| 一区二区三区精品91| 中文字幕av电影在线播放| 久久久国产欧美日韩av| 精品熟女少妇av免费看| 国产国拍精品亚洲av在线观看| 丝袜脚勾引网站| 欧美老熟妇乱子伦牲交| 国产极品粉嫩免费观看在线 | 乱码一卡2卡4卡精品| 精品人妻偷拍中文字幕| 伦理电影免费视频| 亚洲av综合色区一区| 成人亚洲欧美一区二区av| 国产精品人妻久久久影院| 国产精品一区www在线观看| 亚洲精品亚洲一区二区| 久久综合国产亚洲精品| 九九爱精品视频在线观看| 99国产综合亚洲精品| 国产极品粉嫩免费观看在线 | 中文字幕免费在线视频6| 少妇高潮的动态图| 99热网站在线观看| 老司机亚洲免费影院| 中国美白少妇内射xxxbb| 精品久久久久久久久亚洲| 亚洲欧美日韩卡通动漫| 一个人免费看片子| 两个人的视频大全免费| 久久午夜综合久久蜜桃| 成年av动漫网址| 国产女主播在线喷水免费视频网站| 观看美女的网站| 国产亚洲精品第一综合不卡 | 男女边摸边吃奶| 国产成人免费无遮挡视频| 免费人成在线观看视频色| a级毛片在线看网站| 一二三四中文在线观看免费高清| 国产免费又黄又爽又色| 黄片无遮挡物在线观看| av国产久精品久网站免费入址| 日韩一区二区视频免费看| 日本与韩国留学比较| 亚洲精品国产av成人精品| 一区二区日韩欧美中文字幕 | 亚洲中文av在线| av又黄又爽大尺度在线免费看| 国产精品99久久99久久久不卡 | 亚洲激情五月婷婷啪啪| 午夜免费鲁丝| 新久久久久国产一级毛片| freevideosex欧美| 国产高清国产精品国产三级| 九九爱精品视频在线观看| 亚洲av欧美aⅴ国产| 啦啦啦视频在线资源免费观看| 日本欧美视频一区| 中文欧美无线码| 一本大道久久a久久精品| 久久亚洲国产成人精品v| 一区二区av电影网| 精品人妻熟女av久视频| 人妻制服诱惑在线中文字幕| 中文精品一卡2卡3卡4更新| 婷婷成人精品国产| 久久亚洲国产成人精品v| 日韩大片免费观看网站| 日韩亚洲欧美综合| 久久精品国产亚洲av涩爱| 菩萨蛮人人尽说江南好唐韦庄| 国产片特级美女逼逼视频| 精品酒店卫生间| 男的添女的下面高潮视频| 亚洲,一卡二卡三卡| av在线app专区| 蜜桃在线观看..| 日本午夜av视频| 亚洲精品日本国产第一区| 欧美日本中文国产一区发布| 一级,二级,三级黄色视频| 高清av免费在线| 丝袜在线中文字幕| 99re6热这里在线精品视频| 街头女战士在线观看网站| videossex国产| 国产不卡av网站在线观看| 亚洲一级一片aⅴ在线观看| 下体分泌物呈黄色| 久久久久久久久久久丰满| 婷婷色麻豆天堂久久| 国产黄片视频在线免费观看| 涩涩av久久男人的天堂| 91aial.com中文字幕在线观看| 久久久国产一区二区| 欧美 日韩 精品 国产| 91久久精品电影网| 美女内射精品一级片tv| 亚洲欧美成人精品一区二区| 免费不卡的大黄色大毛片视频在线观看| 国产亚洲一区二区精品| 一区二区日韩欧美中文字幕 | 最近中文字幕2019免费版| 青春草亚洲视频在线观看| 麻豆精品久久久久久蜜桃| 美女中出高潮动态图| 精品国产乱码久久久久久小说| 欧美xxxx性猛交bbbb| 99re6热这里在线精品视频| 精品久久国产蜜桃| 香蕉精品网在线| 一级,二级,三级黄色视频| 男男h啪啪无遮挡| 午夜福利,免费看| 国产男女内射视频| 欧美成人午夜免费资源| 成人手机av| 2021少妇久久久久久久久久久| 岛国毛片在线播放| 黑人猛操日本美女一级片| 国产高清三级在线| 搡女人真爽免费视频火全软件| 99热全是精品| 你懂的网址亚洲精品在线观看| 国产成人精品婷婷| 欧美日韩亚洲高清精品| 一级,二级,三级黄色视频| 一本久久精品| 十八禁网站网址无遮挡| 欧美国产精品一级二级三级| av女优亚洲男人天堂| 九草在线视频观看| 成人影院久久| 亚洲三级黄色毛片| 国产精品秋霞免费鲁丝片| 日韩一区二区三区影片| 三级国产精品欧美在线观看| 午夜久久久在线观看| videos熟女内射| 亚洲精品,欧美精品| 午夜激情av网站| 人体艺术视频欧美日本| 久久青草综合色| 免费看不卡的av| 国产国拍精品亚洲av在线观看| av在线app专区| 久久精品国产自在天天线| 久久久久国产精品人妻一区二区| 国产白丝娇喘喷水9色精品| 精品少妇久久久久久888优播| 一级爰片在线观看| 日本wwww免费看| 丰满饥渴人妻一区二区三| 免费高清在线观看日韩| 黄色怎么调成土黄色| tube8黄色片| 九九久久精品国产亚洲av麻豆| 18禁动态无遮挡网站| 国产av码专区亚洲av| 国产免费视频播放在线视频| 18禁观看日本| 成人综合一区亚洲| 久久精品久久久久久噜噜老黄| 性色avwww在线观看| 伦理电影大哥的女人| 亚洲综合色惰| 中文字幕免费在线视频6| 日韩av不卡免费在线播放| a级毛片黄视频| 亚洲国产精品一区二区三区在线| 亚洲欧洲日产国产| 日本欧美视频一区| 99久久中文字幕三级久久日本| 国产深夜福利视频在线观看| 一区在线观看完整版| 国产视频内射| 婷婷色麻豆天堂久久| 久久久久国产精品人妻一区二区| 精品人妻熟女毛片av久久网站| 精品一区在线观看国产| av专区在线播放| 少妇的逼好多水| 亚洲av综合色区一区| 只有这里有精品99| 亚洲四区av| 国产高清不卡午夜福利| 国产免费视频播放在线视频| 看免费成人av毛片| 黄片播放在线免费| 80岁老熟妇乱子伦牲交| 日本猛色少妇xxxxx猛交久久| 日韩视频在线欧美| 日本vs欧美在线观看视频| 国产爽快片一区二区三区| 成人黄色视频免费在线看| 另类精品久久| 国产白丝娇喘喷水9色精品| 黑人高潮一二区| 国产日韩欧美亚洲二区| 美女国产视频在线观看| 亚洲欧美一区二区三区国产| 高清黄色对白视频在线免费看| av不卡在线播放| 熟女电影av网| 夜夜骑夜夜射夜夜干| 欧美日韩综合久久久久久| 毛片一级片免费看久久久久| 九草在线视频观看| 啦啦啦视频在线资源免费观看| 岛国毛片在线播放| 制服人妻中文乱码| 中国美白少妇内射xxxbb| 精品久久国产蜜桃| 黄色欧美视频在线观看| 日韩大片免费观看网站| 99久久人妻综合| 免费高清在线观看视频在线观看| 熟女人妻精品中文字幕| 一级,二级,三级黄色视频| 欧美国产精品一级二级三级| 一区二区日韩欧美中文字幕 | 三级国产精品欧美在线观看| 女性被躁到高潮视频| 亚洲经典国产精华液单| 一级a做视频免费观看| 午夜免费观看性视频| 免费久久久久久久精品成人欧美视频 | 你懂的网址亚洲精品在线观看| 18禁在线播放成人免费| 国产av一区二区精品久久| 狂野欧美激情性xxxx在线观看| 日日啪夜夜爽| 最近最新中文字幕免费大全7| 一个人免费看片子| 国产精品人妻久久久影院| 在线亚洲精品国产二区图片欧美 | 精品99又大又爽又粗少妇毛片| 久久久久精品久久久久真实原创| 欧美xxxx性猛交bbbb| 欧美日韩成人在线一区二区| 999精品在线视频| 日产精品乱码卡一卡2卡三| 一级毛片我不卡| 久久精品人人爽人人爽视色| 中文字幕av电影在线播放| 夫妻性生交免费视频一级片| 亚洲情色 制服丝袜| 十八禁网站网址无遮挡| 久久精品久久精品一区二区三区| 黄色一级大片看看| 国产免费一区二区三区四区乱码| 插逼视频在线观看| 2021少妇久久久久久久久久久| 成人毛片60女人毛片免费| 国产精品久久久久成人av| 人体艺术视频欧美日本| 日韩成人av中文字幕在线观看| 亚洲国产日韩一区二区| 男女高潮啪啪啪动态图| 18+在线观看网站| 精品一区在线观看国产| 人人妻人人爽人人添夜夜欢视频| 久久久久久久精品精品| 国产精品.久久久| 插阴视频在线观看视频| 免费观看在线日韩| 水蜜桃什么品种好| 人成视频在线观看免费观看| 国产在视频线精品| 亚洲成人av在线免费| 免费观看的影片在线观看| 在线天堂最新版资源| 2021少妇久久久久久久久久久| 边亲边吃奶的免费视频| 丝袜脚勾引网站| freevideosex欧美| 国产精品国产三级国产专区5o| 汤姆久久久久久久影院中文字幕| 成人二区视频| 国产伦精品一区二区三区视频9| 亚洲av成人精品一区久久| 亚洲国产精品专区欧美| 日本欧美视频一区| 国产片内射在线| 久久av网站| 91精品伊人久久大香线蕉| 中国美白少妇内射xxxbb| 亚洲欧美日韩卡通动漫| 美女脱内裤让男人舔精品视频| 国产老妇伦熟女老妇高清| av在线观看视频网站免费| 三级国产精品片| av福利片在线| 亚洲图色成人| 五月玫瑰六月丁香| kizo精华| 日日摸夜夜添夜夜爱| 亚洲成人av在线免费| 在线 av 中文字幕| 亚洲欧美精品自产自拍| 欧美xxxx性猛交bbbb| 欧美日韩成人在线一区二区| 在线亚洲精品国产二区图片欧美 | 亚洲av不卡在线观看| 在线观看免费日韩欧美大片 | 大香蕉久久网| 日韩免费高清中文字幕av| 亚洲av二区三区四区| 老司机影院成人| 汤姆久久久久久久影院中文字幕| 久久久久久久大尺度免费视频| 午夜福利在线观看免费完整高清在| 水蜜桃什么品种好| 最黄视频免费看| 午夜福利视频在线观看免费| 国产亚洲欧美精品永久|