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

    基于隨機采樣的頻率域多路徑波場模擬與偏移成像

    2017-06-29 02:17:40
    石油物探 2017年3期
    關(guān)鍵詞:多路徑波場震源

    李 博

    (1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.中國科學(xué)院地質(zhì)與地球物理研究所油氣資源研究重點實驗室,北京100029)

    基于隨機采樣的頻率域多路徑波場模擬與偏移成像

    李 博1,2

    (1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.中國科學(xué)院地質(zhì)與地球物理研究所油氣資源研究重點實驗室,北京100029)

    采用射線追蹤方法的常規(guī)Kirchhoff深度偏移不能全面描述焦散區(qū)的波現(xiàn)象,有可能導(dǎo)致偏移成像產(chǎn)生畸變,利用此類方法進行復(fù)雜介質(zhì)條件下的偏移成像,不可避免地存在焦散問題。詳細推導(dǎo)了二次震源的波場模擬理論方法與區(qū)域分解的實施方案,在確保局部空間不存在射線交叉的前提下進行區(qū)域分解;基于惠更斯二次震源理論進行波場外推,逐步完成整個區(qū)域的波場延拓。提出了基于隨機稀疏采樣與低秩分解的波場模擬高效算法,借助于GPU的計算能力實現(xiàn)了頻率域多路徑波場模擬與成像,大幅提高了計算效率。數(shù)模實驗結(jié)果表明,該方法能正確處理焦散問題,實現(xiàn)多路徑情況下反射界面的正確成像。

    偏移成像;波場外推;焦散;格林函數(shù);低秩分解;隨機采樣

    Kirchhoff積分法深度偏移通常采用的旅行時是單值走時,即一個炮點針對地下的成像點只能有一條傳播路徑,復(fù)雜介質(zhì)條件下,這種方法很難滿足實際需求。若考慮多路徑問題,積分法偏移的理論會遇到焦散問題。需要引入相空間的概念才能區(qū)分多條路徑的波場。很多學(xué)者對此問題進行了深入研究,發(fā)現(xiàn)在計算高頻近似下的體波傳播算法中,關(guān)鍵是要產(chǎn)生空間局部的方向波場才能從根本上解決焦散問題[13]。HILLN[4]和HALE[5]提出高斯束偏移方法在炮集中對不同檢波點進行加窗局部傾斜疊加,利用炮集的τ-p譜進行偏移成像,可以完全解決焦散問題。后來GAO等[6]和SUN等[7]通過選擇τ-p譜的部分能量進行偏移,實現(xiàn)有選擇性的波場傳播的快速束和控制束偏移。李輝等[8]提出了地震數(shù)據(jù)的特征高斯波包(CGP)分解方法,在Radon域選取出地震數(shù)據(jù)的時空和方向特征,大幅減少了高斯波包框架的數(shù)量,實現(xiàn)了方向波場在Gabor域的稀疏化表達。王雄文等[9]將常規(guī)的線性Radon變換(LRT)轉(zhuǎn)化為壓縮感知問題,在壓縮感知的理論框架下實現(xiàn)了高分辨率的方向波分解。這些方法的核心都是如何正確描述空間局部方向波場疊加效應(yīng)的近似表達,以提高波場方向分辨率、能量分解更準(zhǔn)確、特征更突出為研究目標(biāo)。

    本文提出頻率域多路徑疊加的波場模擬及成像方法,無需方向波分解傳播,也能解決焦散問題。首先對震源局部范圍內(nèi)無焦散區(qū)域進行波場模擬,然后將已知波場作為二次震源逐步進行波場外推。實際上,在空間同一個位置完成多次射線的疊加,既避免了射線穿過焦散區(qū)域,又實現(xiàn)了多路徑的波場模擬。實際測試結(jié)果表明本文方法可以有效解決焦散問題,獲得正確的成像結(jié)果。

    1 頻率域多路徑波場模擬基本原理

    基于惠更斯二次震源原理,波場能夠由空間中的一個閉合曲面上的二次震源進行重構(gòu),其重構(gòu)的表達式為Kirchhoff積分公式:

    式中:S表示二次震源的閉合曲面(如圖1所示);r′表示曲面S上的點;U(r′;r0)表示曲面S上的波場;n(r′)表示二次震源r′處曲面S 的單位法向量;G(r′;r)表示射線從r到r′的Green函數(shù)。

    公式(1)經(jīng)常被應(yīng)用在散射波場描述和成像中。在均勻介質(zhì)中,Green函數(shù)G(r′;r)是柱形或球形貝塞爾函數(shù)組成的顯示解析表達形式[10-12],并且滿足平移和旋轉(zhuǎn)不變性。這種類型的波場模擬快速算法已經(jīng)成功應(yīng)用于各個領(lǐng)域,如快速多極算法[13-17]等,實際上,這些算法的本質(zhì)是如何提高Kirchhoff積分核的矩陣運算。對于非均勻介質(zhì)而言,Green函數(shù)G(r′;r)目前沒有解析表達式,且不滿足平移和旋轉(zhuǎn)不變性,因此均勻介質(zhì)中的快速算法不能簡單地推廣到非均勻介質(zhì)的情況。

    本文提出一種非均勻介質(zhì)中Green函數(shù)與Kirchhoff積分數(shù)值計算方法,能夠快速實現(xiàn)公式(1)的波場計算。

    在公式(1)中,U(r′;r0)和是二次震源曲面上的波場函數(shù),首先我們必須計算G(r′;r)和在曲面S 上的值。根據(jù)幾何光學(xué)的高頻漸進假設(shè),G(r′;r)和有如下形式:

    其中,τ表示地震波旅行時。若僅僅保留頻率ω的高階項,忽略低階項可得:

    將公式(4)代入公式(1)得到波場的積分近似表達:

    由于τ(r′;r)滿足程函方程,代入公式(5)可得:

    式中:θ(r′;r)是從r到r′的射線入射方向與曲面S上r′處的單位法向量的夾角(見圖1)。

    圖1 惠更斯二次震源波場構(gòu)建原理

    從公式(5)和公式(6)可見,需要首先計算G(r′;r)和cos[θ(r′;r)]才能最終完成波場的計算。分析公式(5)需要r到r′的走時τ、振幅A和入射角θ(r′;r),射線的起點應(yīng)該從二次震源的r′開始,即射線從r′發(fā)出到r。令θ(r′;r)=π-θ*(r;r′),那么射線從r′到r的出射角為θ*(r;r′),根據(jù)震源與接收點的互易原理,有G(r′;r)=G(r;r′),可以得出:

    代入公式(5)也可獲得波場的積分表達式:

    可見,公式(8)中,利用二次震源曲面上的點r′處的波場U(r′;r0),Green函數(shù)G(r;r′),法向梯度Δr′U(r′;r0)·n(r′)和出射角cos[θ*(r;r′)]就可以實現(xiàn)波場積分核函數(shù)。公式(8)給出了通過已知波場進行外推的方法,只要保證局部范圍內(nèi)無“焦散”現(xiàn)象,就可以通過射線追蹤、高頻近似等常規(guī)方法進行計算。

    2 基于低秩特征分解的惠更斯掃描積分法

    理論上,我們需要在無限大的平面上進行二次震源的積分,而實際上只能在有限范圍內(nèi)的數(shù)據(jù)上實現(xiàn)。根據(jù)局部射線無“焦散”問題,在沒有射線交叉或多次到達的情況下,該物理過程滿足Sommerfeld輻射邊界條件,因此截斷二次震源的無限大平面會對截斷位置的波場產(chǎn)生邊界效應(yīng),但對中間部分的波場影響很小,因此本文方法仍然適用于有限范圍二次震源的情況。下面介紹如何設(shè)定網(wǎng)格參數(shù)并進行波場計算。

    給定速度v(r)和頻率ω,我們可以估計最小波長λmin=2πvmin/ω,這里vmin表示區(qū)域內(nèi)的最小速度,根據(jù)區(qū)域的范圍和最小波長λmin確定采樣點數(shù),保證每個波長范圍內(nèi)有8~10個采樣點即可。本文方法基于幾何光學(xué)理論,相關(guān)的部分內(nèi)容都可以在稀疏網(wǎng)格上進行,當(dāng)涉及到與頻率相關(guān)的計算時才需要在全部的網(wǎng)格上進行。

    離散化的波場積分公式(8)可以表達為如下形式:

    式中:Ms表示二次震源的個數(shù),用表示二次震源集合,用表示需要計算波場的區(qū)域內(nèi)網(wǎng)格點;h表示網(wǎng)格間距 。如圖2所示的密網(wǎng)格部分是需要輸出波場值的區(qū)域,若直接計算波場值將面臨一個O(N3)或O(N5)計算的難題,計算效率難以滿足實際應(yīng)用。

    我們提出一種基于低秩特征分解的Kirchhoff積分算法,可以大幅度提高計算效率并借助GPU的計算能力實現(xiàn)波場的快速計算。算法的實現(xiàn)過程如下。

    波場計算的公式(9)實際上可以簡化為大型矩陣的乘法和加法,為此,首先做如下變量定義:

    其中,f,g,U 為列向量,ρ和ρ1為矩陣。這樣公式(9)可以簡化為:

    從公式(10)可以看出,波場計算轉(zhuǎn)換為兩個矩陣向量積之和。每個矩陣向量積都是大規(guī)模的計算,需要耗費大量的計算機時。因此,我們將問題轉(zhuǎn)化為矩陣向量積的快速算法問題。

    在光滑背景速度假設(shè)條件下,局部區(qū)域內(nèi)的Green函數(shù)具有一定的冗余性,即ρ=[G(ri;sj)]應(yīng)該是一個低秩的稠密矩陣。由于cos[θ*(ri;sj)]也是一個光滑的可壓縮函數(shù),因此ρ1應(yīng)該也是低秩稠密矩陣。下面僅討論公式(10)的第1項,第2項用同樣的方法可以實現(xiàn)快速計算。

    公式(10)的第1項可表達為:

    其中,Green函數(shù)矩陣[G(ri;sj)](i=1,2,…,n;j=1,2,…,m)是由旅行時表、振幅表和出射角余弦表計算得到;[fj](j=1,2,…,m)由二次震源位置的已知波場計算得到。我們的主要任務(wù)是尋找一種近似關(guān)系滿足如下的條件:

    式中:J和I分別表示矩陣G的行和列。如果公式(12)成立,那么不難得出:

    式中:W是一個|J|×|I|的矩陣,由于G是一個低秩矩陣,我們可以找到部分行和列參與計算,縮?。麶|和|I|的長度,從而能夠減少計算量。如果能夠找到矩陣W,那么可以利用公式(13)在精度范圍內(nèi)近似計算出結(jié)果。

    為了得到矩陣W,我們采用隨機采樣的方法,對計算區(qū)域進行劃分,如圖2所示。然后通過隨機采樣的局部空間可以重新組成一個新的矩陣Gs,它是G的一個子矩陣。由于局部區(qū)域內(nèi)的旅行時、振幅、出射角都是低頻光滑函數(shù),因此局部具有信息冗余,導(dǎo)致矩陣Gs是一個低秩的矩陣,在隨機采樣點數(shù)達到一定數(shù)量后將大于矩陣的秩,此時認為Gs的信息已經(jīng)能夠反映這個局部區(qū)域內(nèi)的特征,利用矩陣Gs來構(gòu)造我們需要的矩陣W是可行的。下面詳細介紹矩陣W的獲取方法。

    圖2 二次震源分布與局部區(qū)域的隨機采樣

    假設(shè)Gs是一個隨機采樣后ns×ms的矩陣,包含ns行和ms列數(shù)據(jù),首先對Gs進行可截斷的正交三角(QR)分解:

    其中,Q1為正交矩陣,R1是上三角矩陣,對角線元素分別對應(yīng)矩陣Gs的特征值,并且按照大小降序排列,這樣可以剪裁矩陣R1對角元素小于門檻值e的行和列,R1成為|J|×|J|的矩陣Rc,Q1變?yōu)閚s×|J|的矩陣Qc,令Gc=QcRc,則Gs被裁減為ns×|J|的矩陣,可見通過QR分解可以選擇能夠表述矩陣Gs的一個子矩陣Gc=G(∶,J)。

    采用相似的方法處理Gs的共軛矩陣G*s,同樣利用QR分解處理G*s:

    選取大于門檻值e的行和列得到一個子矩陣G*r:

    那么有:

    式中:Gr和Qr為|I|×ms的矩陣;Rr為|I|×|I|的方陣。裁減后的子矩陣可以幫助我們找到一個Gs的近似矩陣,類似于公式(12):

    那么此時W的矩陣維度已經(jīng)變?yōu)椋麶|×|I|。不難得出利用廣義逆矩陣的奇異值分解方法可以獲得W 的表達形式:

    上述方法我們稱為低秩特征分解方法,要求Gs矩陣具有低秩的特征,才能通過裁減矩陣降低計算量。回到公式(11),光滑速度模型背景條件下,對于單一頻率ω而言,G(ri;sj)是一個大型的低秩稠密矩陣,可以通過上述低秩特征分解的方法實現(xiàn)矩陣向量積的運算:

    實際計算時首先計算公式(20)最右側(cè)的3項矩陣向量乘積,即:

    其中,m為二次震源的個數(shù),其計算量級為|J|×|I|×m。然后計算:

    其中,n為局部區(qū)域內(nèi)的空間離散網(wǎng)格點個數(shù),其計算量為|J|×n。

    這樣總體計算量為:

    若采用公式(11)直接進行矩陣向量積的計算,則計算量為m×n。當(dāng)矩陣G為低秩矩陣時有|J|《m,|I|《n,因此可以大幅度降低整體計算量。此外,公式(21)和(22)的矩陣乘法和加法適合于GPU的多核并行計算模式,從而能夠進一步提高計算效率。

    3 頻率域多路徑疊前偏移成像

    每一種地震波模擬方法都可以開發(fā)研究出一種相應(yīng)的地震偏移成像算法,基本原理就是采用CLEARBOUT提出的波場成像條件[18]:

    式中:I(ω,r)表示頻率域成像結(jié)果;PS(ω,r)和PR(ω,r)分別表示震源波場與接收點波場??紤]炮點震源位置在r0,成像空間位置用r表示的單炮偏移成像問題,將成像空間分成若干水平層狀區(qū)域,如圖3所示。其中,L表示震源與二次震源的層位置。

    圖3 成像空間分解方法

    依據(jù)前文描述的基于低秩特征分解的二次震源波場構(gòu)建方法,并在第1層內(nèi)采用鏡像邊界條件方法實現(xiàn)震源波場表達形式[19-20]。

    當(dāng)在地表炮點位置處,有層位L=0,第1層內(nèi)的震源波場可以表示為:

    在成像空間中二次震源以下的區(qū)域?qū)游籐>0時:

    式中:r′是二次震源位置;nr′表示二次震源個數(shù);h表示成像網(wǎng)格間距;PS(r0)表示在震源r0處的子波函數(shù)。公式(25)和公式(26)形成逐層遞推的波場延拓公式,第1層和其它層的波場構(gòu)建方法不同。

    下面利用同樣的方式構(gòu)建接收點波場表達,不同之處只是把每個接收點當(dāng)作一個二次震源進行波場構(gòu)建。

    當(dāng)在地表接收點位置L=0時存在nR個檢波器接收信號,那么此時的波場有如下表達:

    在成像空間中二次震源以下的區(qū)域L>0時:

    將(25)式、(26)式、(27)式和(28)式代入成像條件公式(24)可以得到大步長分層延拓的成像公式:

    式中:nS表示單炮數(shù)量:NL表示區(qū)域邊界數(shù)量:nω表示離散采樣的頻率數(shù)量。這樣我們可以通過逐層延拓波場,逐層成像的方式完成整個區(qū)域的偏移成像算法,偽代碼流程如圖4所示。

    圖4 多路徑偏移算法偽代碼流程

    4 數(shù)模實驗

    為了驗證多路徑波場模擬方法,我們設(shè)計了按照正弦函數(shù)變化的速度模型(圖5),其變化規(guī)律描述為:

    模型的長度為2.5km,深度為2.5km,成像網(wǎng)格間距為10m,震源位置為(1.0km,0.1km)分別對應(yīng)(x,z)坐標(biāo)。全局射線追蹤可以發(fā)現(xiàn),速度的非均勻性導(dǎo)致射線存在交叉現(xiàn)象,即“焦散”現(xiàn)象。本文選用20Hz頻率參數(shù)進行波場構(gòu)建,對比分析單路徑和多路徑條件下的波場差異,并統(tǒng)計計算時間。依據(jù)本文的多路徑模擬方法,需要進行區(qū)域劃分,我們將兩個由二次震源組成的平面放在如圖5所示的層位位置,保證在每個區(qū)域內(nèi)不發(fā)生焦散現(xiàn)象。圖5中的紅色線條表示從震源發(fā)出的射線路徑,在通過焦散區(qū)的時候射線發(fā)生了交叉現(xiàn)象。圖5中的藍色線條表示二次震源的位置。

    計算結(jié)果如圖6所示,在焦散以外的區(qū)域,單路徑模擬可以得到正確結(jié)果。但在焦散區(qū),多路徑波場十分復(fù)雜,多條射線干涉現(xiàn)象明顯,常規(guī)的單路徑方法已不能獲得正確波場。多路徑波場模擬的正確性在后文的偏移成像中可以得到驗證。

    圖5 數(shù)值速度模型、射線路徑與二次震源的區(qū)域分解方案

    圖6 圖5所示模型的單路徑波場(a)與多路徑波場(b)

    計算效率對比結(jié)果如表1所示。

    表1 計算效率對比結(jié)果

    可見,基于低秩分解的多路徑模擬算法采用CPU計算時效率比單路徑方法有所降低,但采用GPU計算后,計算效率已經(jīng)與單路徑方法相當(dāng),如果直接采用矩陣向量積進行計算,即使采用GPU加速仍然不能獲得理想的效果。因此,本文的低秩分解優(yōu)化算法成功提高了計算效率,達到能夠滿足實際應(yīng)用需求的目的。

    為了測試本文成像方法的優(yōu)缺點,我們設(shè)計了2個層次的實驗,分別是脈沖響應(yīng)實驗和簡單模型實驗。

    首先,利用脈沖響應(yīng)測試算子的優(yōu)劣,分別與同類偏移算法進行對比。我們?nèi)匀徊捎脠D5所示的速度模型做脈沖響應(yīng)計算,在震源位置設(shè)置一個接收點數(shù)據(jù),并給予雙程旅行時延遲時間,然后作為地震接收數(shù)據(jù)輸入,完成成像過程。并將之與單程波算法結(jié)果進行比較,效果如圖7所示。

    從圖7a,圖7c,圖7d,圖7e和圖7f可見,本文方法在焦散區(qū)的成像有“蝴蝶狀”的脈沖響應(yīng),與單程波方法相似,而單路徑成像方法并沒有這種現(xiàn)象。由于單程波方程求解能夠自然解決波場的多路徑問題,因此我們認為“蝴蝶狀”脈沖響應(yīng)是正確結(jié)果,同時也驗證了圖6b中多路徑波場模擬的正確性。

    第二步,我們設(shè)計的速度模型包含一個水平反射層和一個速度異常體,如圖8a所示。模型寬度為1.5km,深度為2.5km,在深度1.7km處有一個速度反射界面,模型中心有一個速度漸變異常區(qū),導(dǎo)致射線彎曲和交叉,出現(xiàn)焦散現(xiàn)象。我們利用有限差分時間域正演模擬單炮記錄,檢波器與震源位于同一深度0.02km,震源的水平坐標(biāo)為0.75km,檢波器分布在整個橫向網(wǎng)格點。其中一個單炮記錄如圖8b所示。

    圖7 脈沖響應(yīng)測試

    圖8 速度模型(a)與單炮記錄(b)

    可見,單炮記錄中有明顯的分層信息,并不是一個簡單的反射同相軸,這正是焦散引起的波場畸變,用常規(guī)的單路徑模擬的成像算法無法獲得真實反射界面的形態(tài)特征。我們分別用單路徑和多路徑模擬方法及單路徑和多路徑偏移成像方法進行驗證,各自的波場與成像結(jié)果如圖9和圖10所示。

    圖9 圖8所示模型的單路徑波場(a)與多路徑波場(b)

    從圖9a可以看出,焦散區(qū)的單路徑波場已經(jīng)發(fā)生畸變,不能正確反映實際波場現(xiàn)象,而本文方法(圖9b)則能夠正確反映焦散區(qū)的波場。從圖10a所示的成像結(jié)果看,單路徑偏移成像方法在異常體下方的成像結(jié)果與實際模型不符。多路徑偏移成像結(jié)果與實際模型位置吻合更好(圖10b),因此本文的多路徑偏移方法對于非均勻介質(zhì)的適應(yīng)性更強。

    圖10 單路徑偏移(a)與多路徑偏移(b)成像結(jié)果

    5 結(jié)論

    本文介紹了一種復(fù)雜介質(zhì)中的基于二次震源多路徑波場模擬方法,該方法采用隨機采樣和低秩分解的Green函數(shù)并利用GPU實現(xiàn)焦散區(qū)波場的快速正確計算,提高了頻率域射線類深度偏移方法的成像精度。數(shù)模實驗結(jié)果顯示:常規(guī)單路徑的Kirchhoff疊前深度偏移成像位置發(fā)生畸變,本文方法能夠獲得正確的成像結(jié)果。該方法在解決復(fù)雜地質(zhì)體的成像問題(如我國西部、南方的復(fù)雜深度成像問題)中有廣泛的應(yīng)用前景。目前,本文方法在推向大規(guī)模實際生產(chǎn)應(yīng)用的過程中還面臨著如地震波的走時、振幅與傳播角度的計算精度問題,以及如何最優(yōu)化設(shè)計二次震源的位置與數(shù)量等理論問題,需要進一步研究。

    致謝:感謝美國密歇根州立大學(xué)錢建良教授在本文研究過程中給予的理論支持和寶貴建議,感謝中國科學(xué)院地質(zhì)與地球物理研究所劉洪研究員提供的大力支持與幫助。

    [1] 蔡杰雄,方伍寶,楊勤勇.高斯束深度偏移的實現(xiàn)與應(yīng)用研究[J].石油物探,2012,51(5):469-475

    CAI J X,F(xiàn)ANG W B,YANG Q Y.Realization and application of Gaussian beam depth migration[J].Geophysical Prospecting for Petroleum,2012,51(5):469-475

    [2] 崔興福,劉衛(wèi)東,劉桂寶,等.平面波偏移、分角度成像與AVA道集生成[J].石油物探,2007,46(6):615-620

    CUI X F,LIU W D,LIU G B,et al.Plan wave migra-tion,angle Imaging and AVA gather production[J].Geophysical Prospecting for Petroleum,2007,46(6):615-620

    [3] 袁茂林,黃建平,李振春,等.局部角度域高斯束偏移參數(shù)優(yōu)選研究[J].石油物探,2015,54(5):602-612

    YUAN M L,HUANG J P,LI Z C,et al.Parameters optimization in local angle-domain Gaussian beam migration[J].Geophysical Prospecting for Petroleum,2015,54(5):602-612

    [4] HILLN R.Gaussian beam migration[J].Geophysics,1990,55(11):1416-1428

    [5] HALE D.Migration by the Kirchhoff slant stack and Gaussian beam methods[R].Center for Wave Phenomena,Colorado School of Mines,1992,CWP-126

    [6] GAO F,ZHANG P,WANG B,et al.Fast beam migration a step toward interactive imaging[J].Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006,2470-2473

    [7] SUN H,SCHUSTER G T.2-D wavepath migration[J].Geophysics,2001,66(5):1528-1537

    [8] 李輝,王華忠,馮波,等.特征高斯波包疊前深度偏移方法[J].地球物理學(xué)報,2014,57(7):2258-2268

    LI H,WANG H Z,F(xiàn)ENG B,et al.Efficient pre-stack depth migration method using characteristic Gaussian packets[J].Chinese Journal of Geophysics,2014,57(7):2258-2268

    [9] 王雄文,王華忠.基于壓縮感知的高分辨率平面波分解方法研究[J].地球物理學(xué)報,2014,57(9):2946-2960

    WANG X W,WANG H Z.A research of high-resolution plane-wave decomposition based on compressed sensing[J].Chinese Journal of Geophysics,2014,57(9):2946-2960[10] CHENG H,CRUTCHFIELD W Y,GIMBUTAS Z,et al.A wideband fast multipole method for the Hemlholtz equation in three dimensions[J].Journal of Computational Physics,2006,216(1):300-325

    [11] CHEW W C,LU C C.The use of Huygens’equivalence principle for solving the volume integral equation of scattering[J].IEEE Transactions on Antennas &Propagation,1993,41(7):897-904

    [12] ROKHLIN V.Diagonal forms of translation operators for the helmholtz equation in three dimensions[J].Applied &Computational Harmonic Analysis,1993,1(1):82-93

    [13] ENGQUISTB,YING L.Fast directional multilevel algorithms for oscillatory kernels[J].SIAM Journal on Scientific Computing,2007,29(4):1710-1737

    [14] GREENGARD L,ROKHLIN V.A fast algorithm for particle simulations[J].Journal of Computational Physics,1987,73(2):325-348

    [15] MESSNER M,SCHANZ M,DARVE E.Fast direc-tional multilevel summation for oscillatory kernels based on Chebyshev interpolation[J].Journal of Computational Physics,2012,231(4):1175-1196

    [16] MICHIELSSEN E,BOAG A.A multilevel matrix decomposition algorithm for analyzing scattering from large structures[J].Antennas &Propagation IEEE Transactions on,1996,44(8):1086-1093

    [17] ROKHLIN V.Rapid solution of integral equations of scattering theory in two dimensions[J].Journal of Computational Physics,1990,86(2):414-439

    [18] CLAERBOUT J F.Toward a unified theory of reflector mapping[J].Geophysics,1971,36(3):467-481

    [19] BERKHOUT A J,WAPENAAR C P A.One-way versions of the Kirchhoff integral[J].Geophysics,1989,54(4):460-467

    [20] BLEISTEIN N.On the imaging of reflectors in the earth[J].Geophysics,1987,52(7):931-942

    (編輯:顧石慶)

    Multipath seismic simulation and imaging in frequency domain based on random sampling method

    LI Bo1,2
    (1.Sinopec Geophysical Research Institute,Nanjing211103,China;2.Key Laboratory of Petroleum Resources Research,Institute of Geology and Geophysics,Chinese Academy of Sciences,Beijing100029,China)

    Conventional Kirchhoff prestack depth migration cannot handle the wave fields with caustic phenomenon;consequently,the subsurface reflector imaging could be distorted.In the complex media,Kirchhoff migration based on ray tracing must be inevitably affected by the caustic problem.In this paper,a multipath wave field simulation method based on the secondary source Huygens theory in complex media is proposed.Firstly,we ensure that there is no ray-cross in the local scope for domain decomposition.Then we use the secondary source Kirchhoff-Huygens integral method to gradually complete wave field extrapolation for the entire imaging area.In this way,we can correctly handle the caustic phenomena and get the correct imaging reflector.The wave field simulation method and domain decomposition strategy of secondary sources definition we described here have included large-scale matrix calculation that is the main bottleneck of efficiency.We propose a wavefield simulation algorithm based on random sampling and low rank decomposition method for large-scale matrix multiplication.On the other hand,we use GPU as the computing devices to realize the multipath seismic simulation and imaging in frequency domain,which greatly improve the computational efficiency.The numerical simulation test results show caustic could be removed and the correct reflector imaging is realized with this method.

    seismic imaging,wave field extrapolation,caustic,Green function,low-rank decomposition,random sampling

    P631

    A

    1000-1441(2017)03-0382-08

    10.3969/j.issn.1000-1441.2017.03.008

    李博.基于隨機采樣的頻率域多路徑波場模擬與偏移成像[J].石油物探,2017,56(3):382-389

    LI Bo.Multipath seismic simulation and imaging in frequency domain based on random sampling method[J].Geophysical Prospecting for Petroleum,2017,56(3):382-389

    2016-10-09;改回日期:2017-03-16。

    李博(1981—),男,博士,高級工程師,主要從事復(fù)雜地震波傳播與成像算法、地球物理方法高性能計算方法研究。

    國家重點研發(fā)計劃項目(2016YFC0601101)資助。

    This research is financially supported by the Nation Key R &D Program of China(Grant No.2016YFC0601101).

    猜你喜歡
    多路徑波場震源
    多路徑效應(yīng)對GPS多普勒測速的影響
    基于5.8G射頻的多路徑識別技術(shù)應(yīng)用探討
    彈性波波場分離方法對比及其在逆時偏移成像中的應(yīng)用
    震源的高返利起步
    交錯網(wǎng)格與旋轉(zhuǎn)交錯網(wǎng)格對VTI介質(zhì)波場分離的影響分析
    基于Hilbert變換的全波場分離逆時偏移成像
    可控震源地震在張掖盆地南緣逆沖斷裂構(gòu)造勘探中的應(yīng)用
    基于5.8GHz多路徑精確識別方案研究
    同步可控震源地震采集技術(shù)新進展
    旋轉(zhuǎn)交錯網(wǎng)格VTI介質(zhì)波場模擬與波場分解
    五月天丁香电影| 纵有疾风起免费观看全集完整版| 精品一区二区三卡| 亚洲精品成人av观看孕妇| 好男人在线观看高清免费视频| 一级二级三级毛片免费看| 国国产精品蜜臀av免费| 各种免费的搞黄视频| 亚州av有码| 国产真实伦视频高清在线观看| 国产成人freesex在线| 老女人水多毛片| 国产人妻一区二区三区在| 观看美女的网站| 菩萨蛮人人尽说江南好唐韦庄| 2022亚洲国产成人精品| 人妻制服诱惑在线中文字幕| 久久久久久九九精品二区国产| 成人高潮视频无遮挡免费网站| 国产免费又黄又爽又色| 国产91av在线免费观看| 国产精品熟女久久久久浪| 亚洲精品国产成人久久av| 久久久久久久久大av| 蜜臀久久99精品久久宅男| 在线天堂最新版资源| 国产欧美亚洲国产| 国产毛片a区久久久久| 在线观看一区二区三区激情| 久久久久久久久久久丰满| 亚洲国产精品成人久久小说| 国产在视频线精品| av播播在线观看一区| 国内精品美女久久久久久| 免费观看在线日韩| 特级一级黄色大片| 少妇的逼水好多| 欧美高清性xxxxhd video| 高清欧美精品videossex| 精品少妇久久久久久888优播| 一本色道久久久久久精品综合| 欧美日韩精品成人综合77777| 精品人妻偷拍中文字幕| 全区人妻精品视频| 免费观看在线日韩| 国产成人freesex在线| 午夜免费观看性视频| 伊人久久国产一区二区| 成人鲁丝片一二三区免费| 水蜜桃什么品种好| 美女内射精品一级片tv| 听说在线观看完整版免费高清| 亚洲av电影在线观看一区二区三区 | 又爽又黄无遮挡网站| 欧美zozozo另类| 国产成人91sexporn| 日韩av不卡免费在线播放| 九九在线视频观看精品| 亚洲精华国产精华液的使用体验| 国产69精品久久久久777片| 亚洲国产精品成人综合色| 国产精品久久久久久精品电影| 日韩人妻高清精品专区| 免费看光身美女| 亚洲久久久久久中文字幕| 国精品久久久久久国模美| 99九九线精品视频在线观看视频| 欧美日韩在线观看h| 国产精品国产三级专区第一集| 男的添女的下面高潮视频| 麻豆成人av视频| a级一级毛片免费在线观看| 熟妇人妻不卡中文字幕| 下体分泌物呈黄色| 日日摸夜夜添夜夜爱| 亚洲av福利一区| 国产一区二区在线观看日韩| 免费观看av网站的网址| 97超碰精品成人国产| 又黄又爽又刺激的免费视频.| 亚洲国产最新在线播放| 久久久久久久午夜电影| 精品国产乱码久久久久久小说| 麻豆久久精品国产亚洲av| 久久精品国产亚洲av涩爱| 国产免费又黄又爽又色| 女人被狂操c到高潮| 热99国产精品久久久久久7| 亚洲欧美日韩卡通动漫| 超碰av人人做人人爽久久| 一级毛片电影观看| 99久久中文字幕三级久久日本| 亚洲精品,欧美精品| 汤姆久久久久久久影院中文字幕| 天天躁夜夜躁狠狠久久av| 夜夜看夜夜爽夜夜摸| 成年版毛片免费区| 伦精品一区二区三区| 亚洲人成网站在线观看播放| 成年女人在线观看亚洲视频 | 黄色配什么色好看| a级一级毛片免费在线观看| 欧美97在线视频| 国产精品麻豆人妻色哟哟久久| 内地一区二区视频在线| 国产精品蜜桃在线观看| 少妇 在线观看| 99久国产av精品国产电影| 禁无遮挡网站| 性色avwww在线观看| 六月丁香七月| 搡老乐熟女国产| 男女边吃奶边做爰视频| 男人舔奶头视频| 国产乱来视频区| 少妇熟女欧美另类| 亚洲欧美日韩另类电影网站 | 午夜福利高清视频| 麻豆久久精品国产亚洲av| 22中文网久久字幕| 日日啪夜夜撸| 街头女战士在线观看网站| 国产精品国产三级专区第一集| 2022亚洲国产成人精品| av在线蜜桃| 在线看a的网站| 久久热精品热| 欧美精品人与动牲交sv欧美| 亚洲精品自拍成人| 在线观看av片永久免费下载| 久久久亚洲精品成人影院| 在线观看av片永久免费下载| 天堂中文最新版在线下载 | 亚洲精品视频女| 麻豆久久精品国产亚洲av| 久久99蜜桃精品久久| 少妇人妻 视频| 少妇的逼好多水| 久久99热这里只频精品6学生| 一区二区三区四区激情视频| 亚洲四区av| 亚洲欧美日韩另类电影网站 | av播播在线观看一区| 亚洲熟女精品中文字幕| 久久久久久久久大av| 国产精品三级大全| 免费观看性生交大片5| 青春草亚洲视频在线观看| 听说在线观看完整版免费高清| 男女国产视频网站| 国产黄a三级三级三级人| 日韩强制内射视频| 日本av手机在线免费观看| 久久人人爽人人爽人人片va| 听说在线观看完整版免费高清| 国产综合精华液| 午夜老司机福利剧场| 丝袜脚勾引网站| 99热这里只有是精品在线观看| 熟妇人妻不卡中文字幕| 99热这里只有是精品50| 日本熟妇午夜| 成人毛片60女人毛片免费| 又大又黄又爽视频免费| 日韩三级伦理在线观看| 性插视频无遮挡在线免费观看| 国产精品久久久久久精品电影| 亚洲内射少妇av| 69人妻影院| 在线 av 中文字幕| 韩国高清视频一区二区三区| 亚洲精品一区蜜桃| 中国美白少妇内射xxxbb| 老司机影院毛片| 亚洲av国产av综合av卡| 欧美bdsm另类| 欧美一区二区亚洲| av在线播放精品| 熟女人妻精品中文字幕| 美女高潮的动态| 又爽又黄a免费视频| 免费黄频网站在线观看国产| 看黄色毛片网站| 三级国产精品欧美在线观看| 亚洲欧美日韩无卡精品| 免费电影在线观看免费观看| 国产国拍精品亚洲av在线观看| 自拍欧美九色日韩亚洲蝌蚪91 | 美女高潮的动态| 99久久人妻综合| 午夜免费男女啪啪视频观看| 中国美白少妇内射xxxbb| 国产成人精品福利久久| 精品国产露脸久久av麻豆| 国产男女超爽视频在线观看| 亚洲国产高清在线一区二区三| 日韩一本色道免费dvd| 久久久久精品久久久久真实原创| 国产成年人精品一区二区| 高清av免费在线| 青春草亚洲视频在线观看| 免费观看性生交大片5| 欧美日韩综合久久久久久| 欧美日韩综合久久久久久| 全区人妻精品视频| 纵有疾风起免费观看全集完整版| 91aial.com中文字幕在线观看| 美女主播在线视频| 大码成人一级视频| 卡戴珊不雅视频在线播放| 最近2019中文字幕mv第一页| 国产一级毛片在线| 国产精品麻豆人妻色哟哟久久| 国产伦理片在线播放av一区| 人体艺术视频欧美日本| 自拍偷自拍亚洲精品老妇| 国产伦精品一区二区三区视频9| 99久久精品热视频| 一级二级三级毛片免费看| 国产在视频线精品| 人妻系列 视频| 久久人人爽人人片av| 狂野欧美白嫩少妇大欣赏| 女人被狂操c到高潮| 男人添女人高潮全过程视频| 国产高清不卡午夜福利| 人妻夜夜爽99麻豆av| 97人妻精品一区二区三区麻豆| 美女国产视频在线观看| 亚洲精品国产色婷婷电影| 大片免费播放器 马上看| 日韩 亚洲 欧美在线| 男女无遮挡免费网站观看| 亚洲经典国产精华液单| 久久精品国产自在天天线| 偷拍熟女少妇极品色| 在线 av 中文字幕| 超碰97精品在线观看| 国产乱人偷精品视频| 岛国毛片在线播放| 亚洲精品日韩在线中文字幕| 亚洲av成人精品一二三区| 啦啦啦中文免费视频观看日本| 一级黄片播放器| 极品少妇高潮喷水抽搐| 在线播放无遮挡| 国产免费一级a男人的天堂| 日韩一本色道免费dvd| 亚洲伊人久久精品综合| 久久久久久久精品精品| 国产成人免费无遮挡视频| av在线播放精品| 亚洲成人一二三区av| 可以在线观看毛片的网站| 国产高清不卡午夜福利| 精品人妻视频免费看| 精品午夜福利在线看| 一二三四中文在线观看免费高清| 国产精品一区www在线观看| 精品亚洲乱码少妇综合久久| 免费看不卡的av| 午夜福利高清视频| 国产高清国产精品国产三级 | 国产 精品1| av在线亚洲专区| 一区二区av电影网| 国模一区二区三区四区视频| 国产 一区 欧美 日韩| 亚洲成色77777| 久久久午夜欧美精品| 国产精品三级大全| 高清欧美精品videossex| 日韩一区二区三区影片| 老师上课跳d突然被开到最大视频| 久久久久久久久久久免费av| 欧美亚洲 丝袜 人妻 在线| 国产亚洲91精品色在线| 精品久久久精品久久久| 亚洲国产成人一精品久久久| 精品视频人人做人人爽| 三级男女做爰猛烈吃奶摸视频| 亚洲真实伦在线观看| 国产男女内射视频| 精品久久久精品久久久| 国产精品久久久久久精品古装| 国产色婷婷99| 少妇高潮的动态图| 久久99精品国语久久久| 一级二级三级毛片免费看| 建设人人有责人人尽责人人享有的 | 亚洲,欧美,日韩| 91久久精品电影网| av又黄又爽大尺度在线免费看| 亚洲va在线va天堂va国产| 在线亚洲精品国产二区图片欧美 | 国产91av在线免费观看| 欧美日本视频| 天堂中文最新版在线下载 | 色哟哟·www| 精品久久久精品久久久| 极品少妇高潮喷水抽搐| 夫妻性生交免费视频一级片| 日韩伦理黄色片| 国产男人的电影天堂91| 18禁在线无遮挡免费观看视频| 久热久热在线精品观看| 不卡视频在线观看欧美| 一区二区三区精品91| 嫩草影院精品99| 亚洲欧美日韩另类电影网站 | av国产精品久久久久影院| 午夜激情福利司机影院| 最近手机中文字幕大全| 女的被弄到高潮叫床怎么办| 欧美人与善性xxx| 欧美亚洲 丝袜 人妻 在线| 日本猛色少妇xxxxx猛交久久| 最近手机中文字幕大全| 大又大粗又爽又黄少妇毛片口| 夫妻性生交免费视频一级片| 国产精品无大码| 久久综合国产亚洲精品| 国产成人aa在线观看| 80岁老熟妇乱子伦牲交| 国产精品一区二区在线观看99| 国产午夜精品久久久久久一区二区三区| 成人特级av手机在线观看| 人妻少妇偷人精品九色| 精品久久久久久久人妻蜜臀av| 欧美精品一区二区大全| 欧美国产精品一级二级三级 | 人妻 亚洲 视频| 插逼视频在线观看| 国产成人免费观看mmmm| 国产成人a∨麻豆精品| 欧美成人精品欧美一级黄| 尾随美女入室| 精华霜和精华液先用哪个| 噜噜噜噜噜久久久久久91| 韩国av在线不卡| 欧美日本视频| 国产乱人视频| 九色成人免费人妻av| 在线 av 中文字幕| 午夜精品一区二区三区免费看| 在线看a的网站| 又粗又硬又长又爽又黄的视频| 五月天丁香电影| 波多野结衣巨乳人妻| 只有这里有精品99| 久久久久久久精品精品| 免费观看的影片在线观看| 国产综合懂色| 日本免费在线观看一区| 亚洲综合色惰| 插阴视频在线观看视频| 精品亚洲乱码少妇综合久久| 网址你懂的国产日韩在线| 少妇被粗大猛烈的视频| 午夜亚洲福利在线播放| 国产69精品久久久久777片| 国产伦理片在线播放av一区| 久久久午夜欧美精品| 激情 狠狠 欧美| 久久人人爽人人片av| 国产精品人妻久久久久久| 免费高清在线观看视频在线观看| 欧美国产精品一级二级三级 | 亚洲欧美成人精品一区二区| 男女边摸边吃奶| 欧美国产精品一级二级三级 | 成年女人看的毛片在线观看| 一级毛片 在线播放| 亚洲综合精品二区| 日本一本二区三区精品| 午夜福利网站1000一区二区三区| 在线观看三级黄色| 久久久久精品性色| 可以在线观看毛片的网站| 国产爱豆传媒在线观看| 亚洲国产欧美人成| 午夜激情福利司机影院| 欧美高清成人免费视频www| 婷婷色麻豆天堂久久| 噜噜噜噜噜久久久久久91| 久久久久性生活片| 日日摸夜夜添夜夜添av毛片| 亚洲精品日韩在线中文字幕| 99九九线精品视频在线观看视频| 成年免费大片在线观看| 麻豆成人av视频| eeuss影院久久| 免费观看无遮挡的男女| 久久精品国产鲁丝片午夜精品| 国产男女超爽视频在线观看| 99久久中文字幕三级久久日本| 高清欧美精品videossex| 99热这里只有是精品50| 亚洲欧美一区二区三区黑人 | 国产伦精品一区二区三区视频9| 在线 av 中文字幕| 国产一区二区三区综合在线观看 | 成年免费大片在线观看| 伦精品一区二区三区| 亚洲美女视频黄频| 久久国内精品自在自线图片| 精品久久久噜噜| 男的添女的下面高潮视频| 中国国产av一级| 国产精品久久久久久久久免| 人人妻人人澡人人爽人人夜夜| 精品午夜福利在线看| 午夜激情福利司机影院| 十八禁网站网址无遮挡 | 成年av动漫网址| 国产有黄有色有爽视频| 一级片'在线观看视频| 亚洲精品第二区| 久久精品国产亚洲av涩爱| 少妇人妻 视频| 亚洲国产色片| 免费看光身美女| xxx大片免费视频| 久久99热6这里只有精品| 最近中文字幕2019免费版| 国产亚洲午夜精品一区二区久久 | 91午夜精品亚洲一区二区三区| 午夜福利视频精品| 美女被艹到高潮喷水动态| 26uuu在线亚洲综合色| 国产亚洲精品久久久com| 欧美变态另类bdsm刘玥| 欧美日韩综合久久久久久| 免费不卡的大黄色大毛片视频在线观看| 亚洲成人久久爱视频| 成年版毛片免费区| 成人特级av手机在线观看| 极品教师在线视频| 亚洲aⅴ乱码一区二区在线播放| 丰满少妇做爰视频| 成人欧美大片| 国产伦精品一区二区三区视频9| 国产免费一区二区三区四区乱码| 99久久精品一区二区三区| 男女下面进入的视频免费午夜| 水蜜桃什么品种好| 一个人看的www免费观看视频| 午夜激情久久久久久久| 热re99久久精品国产66热6| 日本-黄色视频高清免费观看| 美女xxoo啪啪120秒动态图| videos熟女内射| 日本午夜av视频| 国产成人aa在线观看| 久久久久久久久久人人人人人人| 色哟哟·www| 国产在线一区二区三区精| 欧美亚洲 丝袜 人妻 在线| 精品久久久精品久久久| 在线精品无人区一区二区三 | 婷婷色综合大香蕉| 色婷婷久久久亚洲欧美| 成人免费观看视频高清| 有码 亚洲区| 又爽又黄a免费视频| 国产一区二区三区av在线| 国产精品久久久久久精品电影小说 | 久久久久九九精品影院| 伦精品一区二区三区| 久久久国产一区二区| 国产综合懂色| 少妇人妻一区二区三区视频| 亚洲欧美清纯卡通| 人人妻人人爽人人添夜夜欢视频 | 亚洲精品影视一区二区三区av| 青春草视频在线免费观看| av免费观看日本| 久久精品夜色国产| 国产有黄有色有爽视频| av在线天堂中文字幕| 大又大粗又爽又黄少妇毛片口| 免费观看a级毛片全部| 免费观看的影片在线观看| 人妻制服诱惑在线中文字幕| 综合色av麻豆| 99热这里只有是精品50| 亚洲成人中文字幕在线播放| 成年人午夜在线观看视频| 亚洲国产日韩一区二区| 精品久久国产蜜桃| 欧美日韩亚洲高清精品| 午夜日本视频在线| 日韩亚洲欧美综合| 欧美最新免费一区二区三区| 在线观看国产h片| 成人美女网站在线观看视频| 国内揄拍国产精品人妻在线| 99热全是精品| 人妻夜夜爽99麻豆av| 午夜老司机福利剧场| 久久精品国产亚洲av涩爱| 有码 亚洲区| 别揉我奶头 嗯啊视频| 日韩欧美精品v在线| 精品一区二区三区视频在线| 国产黄片美女视频| 一级二级三级毛片免费看| 日韩视频在线欧美| 美女国产视频在线观看| 欧美丝袜亚洲另类| 在线天堂最新版资源| 色播亚洲综合网| 亚洲无线观看免费| 我要看日韩黄色一级片| 80岁老熟妇乱子伦牲交| 亚洲av免费在线观看| 中文字幕制服av| 国产 一区精品| 亚洲国产精品成人久久小说| 全区人妻精品视频| 日日摸夜夜添夜夜爱| 十八禁网站网址无遮挡 | 久久久久久久亚洲中文字幕| 日韩精品有码人妻一区| 国产91av在线免费观看| eeuss影院久久| 久久人人爽人人片av| 我的老师免费观看完整版| 欧美成人一区二区免费高清观看| 国产成人freesex在线| 青春草亚洲视频在线观看| 赤兔流量卡办理| 中文精品一卡2卡3卡4更新| 尤物成人国产欧美一区二区三区| 黄片wwwwww| 成人鲁丝片一二三区免费| 成人毛片a级毛片在线播放| 亚洲国产最新在线播放| 亚洲经典国产精华液单| 午夜免费鲁丝| 日日摸夜夜添夜夜添av毛片| a级毛片免费高清观看在线播放| 亚洲av成人精品一区久久| 国产精品一及| 久久精品国产亚洲av天美| eeuss影院久久| 中文字幕制服av| 日韩欧美精品v在线| 麻豆精品久久久久久蜜桃| 我的老师免费观看完整版| 男男h啪啪无遮挡| 国产伦理片在线播放av一区| 一级毛片电影观看| 成人一区二区视频在线观看| 国产一区亚洲一区在线观看| 尤物成人国产欧美一区二区三区| 26uuu在线亚洲综合色| 男人狂女人下面高潮的视频| 一区二区三区精品91| 成年免费大片在线观看| freevideosex欧美| 亚洲精品第二区| 亚洲精品亚洲一区二区| 狂野欧美激情性bbbbbb| av网站免费在线观看视频| 日韩亚洲欧美综合| 亚洲国产高清在线一区二区三| 热99国产精品久久久久久7| 麻豆久久精品国产亚洲av| 大香蕉97超碰在线| 午夜免费男女啪啪视频观看| 精品午夜福利在线看| videos熟女内射| 国产精品一及| 观看免费一级毛片| 97超视频在线观看视频| 国产综合懂色| 97超视频在线观看视频| 中文在线观看免费www的网站| 美女国产视频在线观看| 国产日韩欧美亚洲二区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 亚洲av不卡在线观看| 日本色播在线视频| 国产精品国产三级国产专区5o| 国产永久视频网站| 丝袜脚勾引网站| 在线观看av片永久免费下载| 国内精品宾馆在线| 制服丝袜香蕉在线| 在线 av 中文字幕| 国产女主播在线喷水免费视频网站| 国产老妇伦熟女老妇高清| 亚洲精品aⅴ在线观看| 国产片特级美女逼逼视频| 久久热精品热| 下体分泌物呈黄色| 男女边摸边吃奶| 天堂俺去俺来也www色官网| 深夜a级毛片| 国产一级毛片在线| 欧美bdsm另类| 直男gayav资源| 亚洲av.av天堂| 国产一区有黄有色的免费视频| 爱豆传媒免费全集在线观看| 成人美女网站在线观看视频| 日本wwww免费看| 97在线人人人人妻| 99热网站在线观看| 涩涩av久久男人的天堂| 亚洲精品乱久久久久久| 中文字幕久久专区| 欧美少妇被猛烈插入视频| 亚洲精华国产精华液的使用体验|