• <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ì)波場模擬與波場分解
    or卡值多少钱| 精品久久久久久久末码| 99热6这里只有精品| 99热这里只有精品一区| 在线十欧美十亚洲十日本专区| 岛国在线观看网站| 亚洲欧美日韩无卡精品| 久久香蕉国产精品| 精品熟女少妇八av免费久了| 级片在线观看| 网址你懂的国产日韩在线| 免费在线观看影片大全网站| 99精品欧美一区二区三区四区| 国产精品自产拍在线观看55亚洲| 亚洲国产色片| 淫妇啪啪啪对白视频| 国产黄片美女视频| 欧美在线黄色| 偷拍熟女少妇极品色| 久久久国产成人免费| 国产av一区在线观看免费| 日日夜夜操网爽| 免费无遮挡裸体视频| 最好的美女福利视频网| 中文字幕熟女人妻在线| 国产精品嫩草影院av在线观看 | 欧美日本亚洲视频在线播放| 国产一区二区亚洲精品在线观看| 国产真实乱freesex| 欧美大码av| 成人高潮视频无遮挡免费网站| 亚洲专区中文字幕在线| xxxwww97欧美| 搡老熟女国产l中国老女人| 老汉色∧v一级毛片| 免费电影在线观看免费观看| 日本一二三区视频观看| 九色国产91popny在线| 搞女人的毛片| 无遮挡黄片免费观看| 国产一区二区亚洲精品在线观看| 黄色片一级片一级黄色片| 国产高清videossex| 国产精品乱码一区二三区的特点| 亚洲国产精品999在线| 少妇的逼好多水| a在线观看视频网站| 一卡2卡三卡四卡精品乱码亚洲| 波多野结衣巨乳人妻| 亚洲中文字幕一区二区三区有码在线看| 男人舔女人下体高潮全视频| 精品久久久久久久毛片微露脸| 毛片女人毛片| 脱女人内裤的视频| 亚洲中文字幕一区二区三区有码在线看| 99久久精品热视频| 一区二区三区免费毛片| 日韩欧美在线二视频| 国产伦一二天堂av在线观看| 亚洲真实伦在线观看| 色综合站精品国产| 国内精品久久久久久久电影| 久久久久免费精品人妻一区二区| 国产爱豆传媒在线观看| 综合色av麻豆| 国产日本99.免费观看| 午夜视频国产福利| 国产精品影院久久| 噜噜噜噜噜久久久久久91| 久久久国产成人免费| 国产成人av激情在线播放| 国产亚洲精品久久久久久毛片| 亚洲美女视频黄频| 一区二区三区免费毛片| 一本精品99久久精品77| 欧美最新免费一区二区三区 | 国产欧美日韩一区二区三| 两个人视频免费观看高清| e午夜精品久久久久久久| 亚洲成av人片在线播放无| 91字幕亚洲| 国产一区二区激情短视频| 成人av在线播放网站| 精品一区二区三区av网在线观看| 精品一区二区三区视频在线观看免费| 国产精品av视频在线免费观看| а√天堂www在线а√下载| 日本黄色片子视频| 国产精品国产高清国产av| 国产一区在线观看成人免费| 欧美黑人巨大hd| 欧美极品一区二区三区四区| 在线a可以看的网站| 久久99热这里只有精品18| 免费搜索国产男女视频| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 亚洲最大成人手机在线| 国产 一区 欧美 日韩| 亚洲专区国产一区二区| 亚洲片人在线观看| 免费观看的影片在线观看| 真实男女啪啪啪动态图| 国产高清视频在线播放一区| 极品教师在线免费播放| 亚洲av电影不卡..在线观看| 麻豆成人av在线观看| 色综合亚洲欧美另类图片| 亚洲五月天丁香| 亚洲avbb在线观看| 亚洲中文字幕日韩| 19禁男女啪啪无遮挡网站| 51国产日韩欧美| 99热精品在线国产| 午夜免费成人在线视频| 91在线精品国自产拍蜜月 | 免费电影在线观看免费观看| 在线观看免费午夜福利视频| 搞女人的毛片| 国产精品免费一区二区三区在线| 成年人黄色毛片网站| 脱女人内裤的视频| 亚洲国产精品sss在线观看| 国产伦人伦偷精品视频| 中文在线观看免费www的网站| 小说图片视频综合网站| 有码 亚洲区| 很黄的视频免费| 无限看片的www在线观看| 人人妻人人澡欧美一区二区| 夜夜爽天天搞| 高潮久久久久久久久久久不卡| 国产精品亚洲一级av第二区| 国产欧美日韩精品亚洲av| 欧美日韩一级在线毛片| 国内精品一区二区在线观看| 精品国产美女av久久久久小说| 男人的好看免费观看在线视频| 国产高清三级在线| 老熟妇仑乱视频hdxx| 亚洲无线在线观看| 又爽又黄无遮挡网站| 激情在线观看视频在线高清| 99热只有精品国产| 在线看三级毛片| 熟妇人妻久久中文字幕3abv| 国产精品久久久人人做人人爽| 乱人视频在线观看| 蜜桃亚洲精品一区二区三区| 少妇高潮的动态图| 宅男免费午夜| 日本黄色片子视频| 亚洲电影在线观看av| 久久精品影院6| 精品一区二区三区视频在线 | 91av网一区二区| 亚洲熟妇熟女久久| 成年人黄色毛片网站| 日日夜夜操网爽| 人人妻,人人澡人人爽秒播| 两个人视频免费观看高清| 亚洲国产中文字幕在线视频| 国产高清视频在线播放一区| 欧美zozozo另类| 级片在线观看| 看片在线看免费视频| 日韩免费av在线播放| 天天添夜夜摸| 色尼玛亚洲综合影院| 在线天堂最新版资源| 老汉色av国产亚洲站长工具| 亚洲第一欧美日韩一区二区三区| 变态另类成人亚洲欧美熟女| 亚洲无线在线观看| 欧美国产日韩亚洲一区| 国产爱豆传媒在线观看| 日韩av在线大香蕉| 亚洲人成网站高清观看| 村上凉子中文字幕在线| 18禁裸乳无遮挡免费网站照片| 日本一本二区三区精品| av国产免费在线观看| 成人无遮挡网站| 怎么达到女性高潮| 久久久久亚洲av毛片大全| 香蕉丝袜av| 国产精品久久电影中文字幕| 又粗又爽又猛毛片免费看| 观看美女的网站| 日本三级黄在线观看| 日本成人三级电影网站| 亚洲国产欧美人成| 国产精品 欧美亚洲| 国产99白浆流出| 床上黄色一级片| e午夜精品久久久久久久| 手机成人av网站| 岛国视频午夜一区免费看| 国产精品美女特级片免费视频播放器| 特级一级黄色大片| 国产精品一区二区三区四区久久| e午夜精品久久久久久久| 亚洲专区国产一区二区| 日韩欧美免费精品| 国产探花在线观看一区二区| 啦啦啦免费观看视频1| tocl精华| 手机成人av网站| 亚洲内射少妇av| 久久草成人影院| 亚洲人成网站在线播放欧美日韩| 欧美性猛交╳xxx乱大交人| 色av中文字幕| 我的老师免费观看完整版| 国内毛片毛片毛片毛片毛片| 99精品在免费线老司机午夜| 国产主播在线观看一区二区| 一进一出抽搐动态| 久久国产精品影院| 嫩草影院精品99| 一卡2卡三卡四卡精品乱码亚洲| 国产精品亚洲一级av第二区| 在线观看美女被高潮喷水网站 | 九色国产91popny在线| 精品电影一区二区在线| 嫩草影院入口| 亚洲国产日韩欧美精品在线观看 | 亚洲人成网站高清观看| 国产真人三级小视频在线观看| 高清毛片免费观看视频网站| 成人高潮视频无遮挡免费网站| 精品乱码久久久久久99久播| 少妇的丰满在线观看| 中文亚洲av片在线观看爽| 久久久国产精品麻豆| 99久国产av精品| 午夜两性在线视频| 免费大片18禁| 一级毛片女人18水好多| 一区福利在线观看| bbb黄色大片| 99久久九九国产精品国产免费| 波多野结衣高清作品| 女人被狂操c到高潮| 成人特级黄色片久久久久久久| 精品久久久久久久毛片微露脸| 国产久久久一区二区三区| 麻豆久久精品国产亚洲av| 18禁在线播放成人免费| 天堂网av新在线| 日韩成人在线观看一区二区三区| 久久久国产成人免费| 手机成人av网站| 国产成人福利小说| 中文在线观看免费www的网站| 波多野结衣高清无吗| 精品久久久久久久末码| 成人高潮视频无遮挡免费网站| 天堂动漫精品| 99热这里只有精品一区| 亚洲精品亚洲一区二区| 亚洲最大成人手机在线| h日本视频在线播放| 欧美日韩中文字幕国产精品一区二区三区| 色尼玛亚洲综合影院| 国产乱人视频| 国产亚洲精品久久久com| 香蕉久久夜色| 国产黄a三级三级三级人| 精品99又大又爽又粗少妇毛片 | 女同久久另类99精品国产91| 日本免费a在线| 久久人妻av系列| 国产午夜精品久久久久久一区二区三区 | 深爱激情五月婷婷| 日韩欧美国产在线观看| 国产真实伦视频高清在线观看 | 欧美成人性av电影在线观看| 色综合婷婷激情| 久久精品国产亚洲av香蕉五月| 美女 人体艺术 gogo| 变态另类成人亚洲欧美熟女| 香蕉丝袜av| 啦啦啦观看免费观看视频高清| 亚洲美女黄片视频| 在线观看免费午夜福利视频| 亚洲国产高清在线一区二区三| 天天添夜夜摸| 窝窝影院91人妻| 精品国产美女av久久久久小说| 欧美丝袜亚洲另类 | 白带黄色成豆腐渣| 51午夜福利影视在线观看| 午夜日韩欧美国产| 90打野战视频偷拍视频| 国产av麻豆久久久久久久| 在线免费观看的www视频| 欧洲精品卡2卡3卡4卡5卡区| 亚洲自拍偷在线| 国产成人aa在线观看| 欧美三级亚洲精品| 十八禁人妻一区二区| 国产精品一区二区三区四区免费观看 | 欧美一级毛片孕妇| 97碰自拍视频| 欧美色视频一区免费| 国产视频一区二区在线看| 岛国在线免费视频观看| 国内精品一区二区在线观看| 岛国在线免费视频观看| 久久久久国产精品人妻aⅴ院| 中文字幕熟女人妻在线| 99久久99久久久精品蜜桃| 少妇的逼好多水| 国产精品一区二区免费欧美| 男女午夜视频在线观看| 51国产日韩欧美| 人妻夜夜爽99麻豆av| 亚洲精品久久国产高清桃花| 黄色日韩在线| 色精品久久人妻99蜜桃| 变态另类成人亚洲欧美熟女| 中国美女看黄片| 国产精品国产高清国产av| 欧美日韩乱码在线| 亚洲成av人片在线播放无| 成年女人永久免费观看视频| 国产伦精品一区二区三区四那| 久久中文看片网| 国产av一区在线观看免费| 日本 av在线| 亚洲人成伊人成综合网2020| 婷婷精品国产亚洲av| 搡女人真爽免费视频火全软件 | 国产一区二区在线观看日韩 | 久久精品91蜜桃| 欧美日本亚洲视频在线播放| 亚洲无线在线观看| 男人和女人高潮做爰伦理| 国产精品久久久久久久久免 | 国产一区二区亚洲精品在线观看| 蜜桃亚洲精品一区二区三区| 韩国av一区二区三区四区| 全区人妻精品视频| 97人妻精品一区二区三区麻豆| 老鸭窝网址在线观看| 啦啦啦观看免费观看视频高清| 国产av在哪里看| 久久九九热精品免费| 中文字幕av在线有码专区| 少妇丰满av| 狂野欧美激情性xxxx| 成人特级黄色片久久久久久久| 欧美一区二区国产精品久久精品| 怎么达到女性高潮| 国内揄拍国产精品人妻在线| 久久九九热精品免费| 在线天堂最新版资源| 丰满乱子伦码专区| 色尼玛亚洲综合影院| 男女做爰动态图高潮gif福利片| 欧美日韩福利视频一区二区| 国产精品1区2区在线观看.| 搡老妇女老女人老熟妇| 成人欧美大片| www.色视频.com| 校园春色视频在线观看| 亚洲精品亚洲一区二区| 日韩欧美精品免费久久 | 91九色精品人成在线观看| 国产欧美日韩一区二区三| 精品无人区乱码1区二区| 在线观看av片永久免费下载| 成人无遮挡网站| 国产伦精品一区二区三区四那| 有码 亚洲区| svipshipincom国产片| 99久久久亚洲精品蜜臀av| 一本一本综合久久| 白带黄色成豆腐渣| 欧美黄色片欧美黄色片| 国产激情欧美一区二区| 搡女人真爽免费视频火全软件 | 最新美女视频免费是黄的| 午夜福利高清视频| 亚洲国产精品成人综合色| 观看免费一级毛片| 十八禁网站免费在线| 亚洲第一欧美日韩一区二区三区| 午夜a级毛片| 国产中年淑女户外野战色| 夜夜看夜夜爽夜夜摸| 亚洲成av人片免费观看| 白带黄色成豆腐渣| 久久久精品欧美日韩精品| 国产亚洲精品一区二区www| 欧美午夜高清在线| 国产麻豆成人av免费视频| 好男人电影高清在线观看| 久久精品91无色码中文字幕| 午夜免费男女啪啪视频观看 | 深夜精品福利| 国产精品亚洲av一区麻豆| 欧美又色又爽又黄视频| 青草久久国产| 国内精品一区二区在线观看| 久久亚洲真实| 亚洲国产日韩欧美精品在线观看 | 亚洲美女视频黄频| 久久久久性生活片| 国产伦在线观看视频一区| 男女床上黄色一级片免费看| 99精品在免费线老司机午夜| 黄片小视频在线播放| 欧美色欧美亚洲另类二区| 亚洲人成网站在线播放欧美日韩| 久久久久久久精品吃奶| 久久国产精品影院| 亚洲aⅴ乱码一区二区在线播放| 一级毛片女人18水好多| 舔av片在线| 脱女人内裤的视频| 琪琪午夜伦伦电影理论片6080| 久久久久精品国产欧美久久久| av天堂在线播放| 国产精品久久久久久久久免 | 精品一区二区三区av网在线观看| 亚洲欧美日韩无卡精品| 亚洲国产色片| 国产三级在线视频| 麻豆成人午夜福利视频| 99国产极品粉嫩在线观看| 亚洲成av人片免费观看| 一个人看的www免费观看视频| 久久精品91无色码中文字幕| 日本熟妇午夜| 国产精品久久久久久人妻精品电影| 小蜜桃在线观看免费完整版高清| 91麻豆精品激情在线观看国产| 日本一本二区三区精品| 香蕉av资源在线| 日本黄色片子视频| 久久精品国产亚洲av香蕉五月| 久久久国产成人精品二区| 中文字幕人成人乱码亚洲影| 在线观看免费午夜福利视频| 一区二区三区激情视频| 在线免费观看不下载黄p国产 | 欧美大码av| 成人三级黄色视频| 精品熟女少妇八av免费久了| 亚洲乱码一区二区免费版| tocl精华| 亚洲美女黄片视频| 欧美丝袜亚洲另类 | 欧美中文综合在线视频| 国产三级在线视频| 国产欧美日韩精品一区二区| 午夜福利在线在线| 亚洲国产色片| 欧美三级亚洲精品| 在线看三级毛片| 亚洲专区中文字幕在线| 俺也久久电影网| 精品国产亚洲在线| 嫩草影院入口| ponron亚洲| 男人的好看免费观看在线视频| 精品福利观看| 久9热在线精品视频| 国产视频一区二区在线看| 欧美日韩一级在线毛片| 精品福利观看| 亚洲七黄色美女视频| 成年版毛片免费区| а√天堂www在线а√下载| 在线观看免费视频日本深夜| 欧美黑人巨大hd| netflix在线观看网站| 国产精品一区二区三区四区久久| 尤物成人国产欧美一区二区三区| 99在线人妻在线中文字幕| 久久久精品大字幕| 亚洲 国产 在线| 婷婷丁香在线五月| 亚洲精品成人久久久久久| 久久久久精品国产欧美久久久| 3wmmmm亚洲av在线观看| 久久久久性生活片| 日本与韩国留学比较| 啪啪无遮挡十八禁网站| 一区二区三区国产精品乱码| 老司机深夜福利视频在线观看| 久久精品综合一区二区三区| 久久精品国产亚洲av涩爱 | 可以在线观看毛片的网站| 欧美性感艳星| 亚洲成av人片在线播放无| www.999成人在线观看| 最新美女视频免费是黄的| 久久精品国产亚洲av香蕉五月| 免费在线观看日本一区| 国产精品一区二区三区四区久久| 国产精品99久久久久久久久| 欧美性感艳星| 国产高清videossex| 麻豆国产97在线/欧美| 高潮久久久久久久久久久不卡| 国产精品爽爽va在线观看网站| 国产av麻豆久久久久久久| 日韩欧美一区二区三区在线观看| 麻豆一二三区av精品| 国产亚洲欧美在线一区二区| 成人国产综合亚洲| 18禁美女被吸乳视频| 国产99白浆流出| 久久久色成人| 日本黄色视频三级网站网址| 俺也久久电影网| 日日摸夜夜添夜夜添小说| 免费观看人在逋| 日韩欧美在线乱码| 一区二区三区免费毛片| 三级男女做爰猛烈吃奶摸视频| 一进一出抽搐动态| 亚洲 国产 在线| 在线观看日韩欧美| 最新中文字幕久久久久| 欧美极品一区二区三区四区| 国产成人福利小说| 一区二区三区国产精品乱码| 最近最新免费中文字幕在线| 波野结衣二区三区在线 | 欧美日韩亚洲国产一区二区在线观看| 女人高潮潮喷娇喘18禁视频| 免费在线观看成人毛片| 18禁国产床啪视频网站| 午夜久久久久精精品| 90打野战视频偷拍视频| 韩国av一区二区三区四区| 十八禁人妻一区二区| 亚洲 欧美 日韩 在线 免费| 国产视频内射| 国产三级在线视频| 亚洲在线自拍视频| 少妇的逼好多水| 老熟妇乱子伦视频在线观看| 精品免费久久久久久久清纯| 成人av一区二区三区在线看| tocl精华| 男人舔女人下体高潮全视频| 国产精品久久电影中文字幕| a级毛片a级免费在线| 久久久久免费精品人妻一区二区| 男人的好看免费观看在线视频| 国产极品精品免费视频能看的| 成年免费大片在线观看| 国产午夜精品久久久久久一区二区三区 | 乱人视频在线观看| 亚洲无线观看免费| 中文字幕熟女人妻在线| 亚洲国产高清在线一区二区三| www.www免费av| 村上凉子中文字幕在线| 成熟少妇高潮喷水视频| 3wmmmm亚洲av在线观看| 国内久久婷婷六月综合欲色啪| 亚洲精品粉嫩美女一区| 国产亚洲精品综合一区在线观看| 久久精品综合一区二区三区| 国产精品亚洲一级av第二区| 桃红色精品国产亚洲av| 美女免费视频网站| 久久伊人香网站| 国产精品爽爽va在线观看网站| a级一级毛片免费在线观看| 母亲3免费完整高清在线观看| 又紧又爽又黄一区二区| 日韩av在线大香蕉| 国产精品影院久久| 国产69精品久久久久777片| 波多野结衣高清作品| 精品久久久久久久人妻蜜臀av| 国模一区二区三区四区视频| 夜夜看夜夜爽夜夜摸| 欧美成人免费av一区二区三区| 两个人的视频大全免费| 麻豆国产av国片精品| 三级男女做爰猛烈吃奶摸视频| 国产精品女同一区二区软件 | 日韩欧美三级三区| 国产探花极品一区二区| 最近最新中文字幕大全免费视频| 1000部很黄的大片| 在线观看免费午夜福利视频| 午夜福利高清视频| 国产97色在线日韩免费| 黄片小视频在线播放| 少妇的逼水好多| 国产亚洲av嫩草精品影院| 在线观看av片永久免费下载| 制服人妻中文乱码| 久久精品影院6| 又紧又爽又黄一区二区| 1024手机看黄色片| 色综合亚洲欧美另类图片| 最新美女视频免费是黄的| 国产高清有码在线观看视频| 国产精品1区2区在线观看.| bbb黄色大片| 国产午夜精品久久久久久一区二区三区 | 97超视频在线观看视频| 天堂网av新在线| 国产国拍精品亚洲av在线观看 | 日韩欧美 国产精品| 18禁黄网站禁片午夜丰满|