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

    基于新的差分結構的時-空域高階有限差分波動方程數值模擬方法

    2016-06-30 07:39:18張保慶周輝陳漢明盛善波
    地球物理學報 2016年5期
    關鍵詞:數值模擬

    張保慶, 周輝, 陳漢明, 盛善波

    1 油氣資源與探測國家重點實驗室, CNPC物探重點實驗室, 中國石油大學(北京), 北京 102249 2 東方地球物理公司研究院大港分院, 天津大港 300280 3 中國石油天然氣勘探開發(fā)公司, 北京 100034

    基于新的差分結構的時-空域高階有限差分波動方程數值模擬方法

    張保慶1, 2, 周輝1*, 陳漢明1, 盛善波3

    1 油氣資源與探測國家重點實驗室, CNPC物探重點實驗室, 中國石油大學(北京), 北京102249 2 東方地球物理公司研究院大港分院, 天津大港300280 3 中國石油天然氣勘探開發(fā)公司, 北京100034

    摘要傳統(tǒng)的高階有限差分波動方程數值模擬方法采用高階差分算子近似空間偏導數,能有效抑制空間頻散.然而,傳統(tǒng)的有限差分法僅采用二階差分算子近似時間偏導數,這使得地震波場沿時間外推的精度較低.當采用較大的時間采樣間隔,傳統(tǒng)的有限差分法模擬波場會出現明顯的時間頻散,甚至不穩(wěn)定.本文基于新的差分結構和中心網格剖分,發(fā)展了一種空間任意偶數階精度、時間四階和六階精度的時空域有限差分方法.基于對離散后的頻散關系進行泰勒展開,本文推導了時空域高階有限差分算子的差分系數.相速度分析表明時間四階、六階精度的差分方法能顯著地減小傳統(tǒng)時間二階精度差分方法的時間頻散.在相同的精度下與傳統(tǒng)差分法比較,本文發(fā)展的時間四階、六階有限差分方法的計算效率比傳統(tǒng)方法高.均勻和非勻均介質中的波場數值模擬實驗進一步證實本文研究的時空高階有限差分方法的優(yōu)越性.

    關鍵詞波動方程; 時空域; 高階有限差分; 數值模擬

    1引言

    波動方程數值模擬是當前熱門的逆時偏移、最小二乘偏移和全波形反演技術的基礎.波動方程數值模擬的精度和效率在很大程度上決定了波動方程偏移和反演的精度和效率.當前,波動方程數值模擬的方法主要包括基于傅里葉變換的偽譜法,基于非規(guī)則網格剖分的有限元類方法以及基于差分近似的有限差分方法.偽譜法采用傅里葉變換計算空間偏導數,能完全壓制由于空間離散造成的數值頻散,但其缺點是計算效率較低(謝桂生等,2005),并且由于傅里葉變換的全局性,吸收邊界的處理較為麻煩.有限元類方法剖分網格更靈活,能準確地模擬各種不規(guī)則邊界,從而避免矩形網格剖分造成的階梯繞射.然而,受制于巨大的計算量,目前有限元類方法在逆時偏移和全波形反演中很少使用.相比之下,有限差分法無疑是一種更靈活簡便的數值算法,已被廣泛用于各類偏微分方程的數值求解,在波動方程偏移和反演中也大受歡迎.

    Alterman and Karal(1968)最早采用有限差分法模擬二維層狀介質中的彈性波場,Kelly et al.(1976)則討論了彈性波動方程數值模擬的兩種有限差分方法,分別對應中心網格有限差分法和交錯網格有限差分法.Virieux(1984, 1986)正式將Yee(1966)的交錯網格離散方法引入到波動方程數值模擬中,并求解一階速度-應力形式的聲波和彈性波動方程.早期的交錯網格有限差分法僅采用二階差分算子分別近似空間和時間偏導數項,當采用較大的網格間距時,數值頻散嚴重,需要校正(吳國忱和王華忠,2005).為了提高有限差分數值模擬的精度,Levander(1988)采用四階精度的差分算子近似空間偏導數,顯著地改善了二階差分算子的計算精度,形成了空間四階精度、時間二階精度的差分方法,本文將其簡記為(4,2)差分法.與二階差分方法相比,(4,2)差分方法更好地平衡了計算精度和計算效率,因此被廣泛用于地震波場的數值模擬(Graves, 1996;Moczo et al., 2000;Shi et al., 2012).為了進一步提高空間差分算子的精度,許多學者進一步討論了高階有限差分數值模擬方法(Holberg, 1987;Fornberg, 1987;Kristek et al. 2010;Dong et al., 2013;杜啟振等,2015).通過增加差分算子的長度,空間精度可以達到任意偶數(2M)階(劉洋等,1988),本文將傳統(tǒng)的高階有限差分方法記為(2M,2).

    盡管傳統(tǒng)的高階有限差分方法能有效地抑制空間頻散,但長期以來很少有學者關注時間頻散.事實上,地震波在時間和空間里傳播,因網格離散造成的數值頻散同時包括空間頻散和時間頻散,傳統(tǒng)的高階差分(2M,2)方法時間外推的精度只有二階,當采用較大的時間步長時會出現明顯的時間頻散,甚至不穩(wěn)定.近些年,許多旨在提高地震波場時間外推精度的數值模擬方法得以發(fā)展.這些方法主要包括基于傅立葉變換的矩陣低秩分解方法(Fomel et al., 2013;Song et al., 2013;Fang et al., 2014)和基于頻散關系的時-空域有限差分方法(Finkelstein and Kastner, 2007;Tan and Huang, 2014).基于矩陣低秩分解的數值模擬方法在波場外推的每個時間步至少需要兩次反傅里葉變換,計算量較大.因此,Song等(2013)和Fang等(2014)在矩陣低秩分解的基礎上分別發(fā)展了基于中心網格和基于交錯網格的有限差分方法,進一步提高了計算效率.Tan和Huang(2014)基于新的差分結構發(fā)展了時間四階和六階精度、空間任意偶數階精度的交錯網格有限差分方法,分別記為(2M,4)和(2M,6).頻散分析表明(2M,4)和(2M,6)交錯網格有限差分方法能顯著地提高傳統(tǒng)(2M,2)交錯網格有限差分方法的精度,并且比傳統(tǒng)方法具有更好的穩(wěn)定性.國內學者董良國等(2000)在求解彈性波動方程的一階速度-應力形式時將速度分量對時間的奇數階偏導數轉化成應力分量對空間的偏導數,發(fā)展了時間四階精度、空間任意偶數階精度的交錯網格有限差分法.楊慶節(jié)等(2014)推導了可導函數任意次導數的任意階精度的差分近似及相應的差分系數,進一步完善了董良國等(2000)的高階差分法.

    高階差分算子的差分系數可利用泰勒多項式展開進行確定,也可采用一些優(yōu)化的方法確定差分系數,如基于最小二乘的有限差分方法(Liu et al., 2014; 張杰和吳國忱,2015),基于正則化的有限差分方法(Wang et al., 2014)以及基于非線性優(yōu)化的有限差分方法(梁文全等,2015).時-空域有限差分方法的差分系數隨速度而變換,當速度劇烈變化時,基于優(yōu)化的方法需要對每個不同的速度值進行一次優(yōu)化,引入了額外的計算量.相比之下,基于泰勒展開的方法給出了有限差分系數的解析表達式,使得計算差分系數的計算量基本可以忽略.

    本文基于新的差分結構,發(fā)展了基于中心網格剖分的空間任意2M階,時間四階和六階精度的有限差分方法,與Tan和Huang(2014)討論的(2M,4)和(2M,6)交錯網格有限差分法相比,本文的差分方法適用于離散二階偏導數,因此適用于模擬二階波動方程.當不考慮介質的密度變化時,本文研究的(2M,4)和(2M,6)中心網格有限差分方法具有更高的計算效率.

    2波動方程的差分離散格式

    2.1時間四階精度的差分方法

    地震波在二維聲學介質中傳播的控制方程為

    (1)

    其中t表示時間,xoz組成二維(2D)笛卡爾坐標系,v為地震波的傳播速度,p為壓力.采用如下差分格式離散方程(1),公式為

    (2)

    其中上標表示時間網格點的標號,下標表示空間網格點的標號,a0,0,a1,1和am,0(m=1,2,…,M)表示差分系數,h為網格大小,Δt為時間步長,M為空間差分算子長度.方程(2)為一種新的差分結構,與傳統(tǒng)的差分結構相比,額外包含了非軸向的四個空間網格節(jié)點.Tan和Huang(2014)的研究表明,額外包含非軸向的四個網格節(jié)點剛好使得差分近似的截斷誤差為o(Δt4)+o(h2M)(符號o表示高階無窮小),因此時間差分能達到四階精度,空間差分能達到2M階精度.本文將基于新的差分結構(2)并采用泰勒展開的方法確定方程(2)中的差分系數.對方程(2)兩端做傅里葉變換,并利用傅里葉變換的時移性質可得到如下頻散關系為

    (3)

    其中γ=vΔt/h為網格比,kx,kz為2D離散波數分量,ω為角頻率.利用ω=vk替換ω并對方程(3)中的余弦函數做泰勒展開得到:

    (4)

    (5)

    將方程(5)代入(4)并按h2j項整理得到

    (6)

    對比方程(6)兩端h2j項的系數可得到:

    (7)

    (8)

    (9)

    (10)

    將方程(10)代入方程(7)和(9)求得差分系數為

    (11)

    利用方程(11)中的差分系數,差分格式(2)能達到空間2M階精度,時間四階精度.

    2.2時間六階精度的差分方法

    Tan and Huang (2014)發(fā)展了時間六階精度的交錯網格有限差分法,與其時間四階精度的差分方法相比,其差分結構里又額外包含了八個網格節(jié)點.采用泰勒展開確定差分系數時,額外的網格節(jié)點使得求取差分系數的方程適定(即方程的個數剛好等于差分系數的個數),并且差分截斷誤差為o(Δt6)+o(h2M),因此時間差分能達到六階精度,空間差分能達到2M階精度.類似地,為發(fā)展時間六階精度的中心網格差分方法,本文采用如下的差分格式為

    (12)

    其中d0,0,d1,2和dm,0(m=1,2,…,M)為差分系數.類似地,對方程(12)兩端做傅里葉變換,并利用其時移性質得到如下頻散關系為

    (13)

    將ω=vk代入方程(13),并對方程(13)兩端的余弦函數做泰勒展開得到

    (14)

    對比(14)式兩端h2j項的系數可得到:

    (15)

    (16)

    γ2j-2, (j=1,2,…,M).

    (17)

    為獲得時間外推的六階精度,方程(16)兩端混

    (18)

    聯立方程(15)、(17)和(18)得到空間任意偶數階,時間六階精度的差分系數為

    (19)

    3頻散分析

    根據頻散關系(3)和(13)求得角頻率ω,進而得到數值相速度為vnum=ω/k,定義數值相速度和真實地震波傳播速度之間的比值δ=vnum/v,用該比值來描述頻散程度.該比值越接近于1,頻散越小,越遠離1頻散越嚴重.另外,空間頻散造成地震波的相位滯后,即δ<1;而時間頻散則造成地震波的相位超前,即δ>1(Levander,1988;Tan and Huang,2014).據此可大致區(qū)分時間頻散和空間頻散.根據以上思路,時間四階差分算子的δ的表達式為

    (20)

    而時間六階差分算子的δ的表達式為

    (21)

    圖1顯示了當網格比取不同值時,傳統(tǒng)時間二階差分方法和本文推導的時間四階、六階差分方法的頻散誤差示意圖,空間差分算子長度統(tǒng)一取為2M=16.當網格比取較小值γ=0.2時,傳統(tǒng)的時間二階差分方法仍然存在較明顯的時間頻散,如圖1a所示.當網格比增加到γ=0.4(實際上是時間采樣步長增加一倍),傳統(tǒng)的方法頻散誤差顯著增加,如圖1b所示.相比之下,時間四階精度的差分方法在γ=0.2時頻散誤差很小,當時間步長增加一倍時,頻散誤差也沒有顯著增加,如圖1c和圖1d所示.對比圖1d和圖1a也可得出結論,當傳統(tǒng)方法的時間步長取為時間四階精度差分方法的一半時,其頻散誤差仍然比時間四階精度的差分方法大.圖1e和圖1f則分別顯示了γ=0.2和γ=0.4時時間六階精度差分方法的頻散誤差,很明顯,時間六階精度的差分方法壓制時間頻散的效果最好.圖1g和圖1h則分別顯示了γ=0.6時本文研究的時間四階和六階精度差分算子的頻散誤差,在高波數部分四階差分算子出現了較明顯的頻散,而六階差分算子的精度更高,誤差更小.對比圖1h和圖1a可以得出結論,即使傳統(tǒng)時間二階差分方法所取的時間步長為時間六階差分方法的1/3,其頻散誤差仍然比后者嚴重.

    4穩(wěn)定性分析

    (22)

    取尼奎斯特波數(kx,kz)=(π/h,π/h)代入方程(22)得到:

    (23)

    將方程(11)中a0,0的表達式代入(23)式,經過整理得到(2M, 4)差分方法的穩(wěn)定性條件為

    (24)

    其中int表示直接截取浮點數的整數部分.同理可得到(2M,6)差分方法的穩(wěn)定性條件為

    圖1 頻散示意圖(a) 傳統(tǒng)(16,2)方法,γ=0.2; (b) 傳統(tǒng)(16,2)方法,γ=0.4; (c) 新的(16,4)方法,γ=0.2; (d) 新的(16,4)方法γ=0.4; (e) 新的(16,6)方法,γ=0.2; (f) 新的(16,6)方法γ=0.4; (g) 新的(16,4)方法γ=0.6; (h) 新的(16,6)方法γ=0.6.Fig.1 Dispersion relations(a) The traditional (16,2) scheme with γ=0.2 and (b) with γ=0.4; (c) The new (16,4) scheme with γ=0.2 and (d) with γ=0.4; (e) The new (16,6) scheme with γ=0.2 and (f) with γ=0.4, (g) The new (16,4) scheme with γ=0.6, and the new (16,6) scheme with γ=0.6.

    (25)

    圖2 傳統(tǒng)時間二階差分方法,新的時間四階、六階差分方法取不同差分算子長度時對應的最大網格比Fig.2 The allowed maximum Courant-Friedrichs-Lewy (CFL) number with ranging stencil length from 4 to 32, for the traditional temporal 2nd-order, and our temporal 4th- and sixth-order FD schemes.

    圖2進一步比較了取不同差分算子長度時,傳統(tǒng)的時間二階差分,本文所研究的時間四階和六階差分方法能取得的最大網格比.從圖中可以看出,時間四階和六階差分方法顯著地改善了傳統(tǒng)二階差分方法的穩(wěn)定性,時間六階差分算子的穩(wěn)定性比時間四階差分算子的穩(wěn)定性略好.

    5計算效率分析

    時間四階和六階精度差分算子基于新的差分結構,該差分結構除了包含軸向2M+1個離散點外,還包含少數的非軸向離散點,與傳統(tǒng)的差分格式相比,引入了額外的計算量.然而,時間高階差分方法改善了傳統(tǒng)二階方法的精度,在相同的精度下比較,傳統(tǒng)方法需要使用更小的時間步長.整體而言,時間高階差分方法能提高計算效率.通過一系列的頻散分析可以得到,本文討論的時間四階差分算子采用約2.7倍于傳統(tǒng)方法的時間步長時,仍然達到相似的精度;而時間六階差分算子則可采用約4倍于傳統(tǒng)方法的時間步長時,仍然能達到可比擬的精度.表1列出了傳統(tǒng)的時間二階差分算子、本文討論的時間四階、六階差分算子包含的離散點數,每個時間步內計算空間偏導數的浮點運算次數以及達到相似精度時所采用的時間步長.當差分算子長度2M=16時,與傳統(tǒng)的時間二階差分方法相比,本文發(fā)展的時間四階差分算子和時間六階差分算子的理論加速比分別為2.4和3.0.

    表1 傳統(tǒng)時間二階差分算子與時間四階、

    6數值例子

    本文的第一個數值模擬例子模擬均勻介質中的聲波波場,介質的地震波傳播速度為2000 m·s-1,網格大小為20 m,采用主頻為12.5 Hz的雷克子波作為震源,置于模型的中間,差分算子的長度2M=16.圖3顯示了采用不同差分方法模擬得到的3.0 s的波場切片.圖3a采用傳統(tǒng)的時間二階差分算子計算,時間步長為5.0 ms,模擬結果出現了顯著的相位扭曲,說明傳統(tǒng)的時間二階差分方法的精度低;圖3b是時間步長減半為2.5 ms時傳統(tǒng)方法模擬的結果,其數值頻散仍然比較明顯;進一步減小時間步長為2 ms,傳統(tǒng)的時間二階精度差分法的精度有所提高,但仍然可見微弱的時間頻散,如圖3c所示;圖3d顯示的波場切片為傳統(tǒng)的二階方法模擬,采用的時間步長已經縮小4倍達到1.25 ms,未見明顯的數值頻散.圖3e和3f是分別采用本文發(fā)展的時間四階、六階差分方法模擬得到的波場切片,采用時間步長為5.0 ms.對比圖3b和3e可以發(fā)現,盡管傳統(tǒng)的二階差分方法采用的時間步長為時間四階差分方法的一半,但其數值頻散仍然比四階方法嚴重;對比圖3c和3e可以發(fā)現,當四階差分方法采用的時間步長為傳統(tǒng)方法的2.5倍時,其壓制數值頻散的效果仍然更好,這與前文的頻散分析結論是一致的.進一步對比圖3d和3f,兩圖中的波場均未見明顯的數值頻散,而此時時間六階差分方法采用的時間步長為傳統(tǒng)二階方法的4倍,這與前文的頻散分析結論也是一致的.

    本文的第二個數值算例模擬在二維BP鹽丘中傳播的聲波波場,圖4為相應的速度模型,該模型的最大速度差異比約為3.4.離散采用的網格大小為15 m,采用主頻為12.5 Hz的雷克子波作為震源,位于(x,z)=(7500,300) m處,數值模擬采用的差分長度為2M=16.分別采用傳統(tǒng)的時間二階差分方法、本文發(fā)展的時間四階、六階差分方法模擬聲波的傳播.為對比各種差分方法模擬的結果,將傳統(tǒng)時間二階差分方法采用微小的時間步長0.15 ms模擬的地震記錄作為參考解(見圖5).

    圖3 采用不同差分方法模擬均勻介質中3 s時的波場切片,(a)傳統(tǒng)時間二階差分方法的模擬結果,時間步長5.0 ms,(b)傳統(tǒng)時間二階差分方法的模擬結果,時間步長2.5 ms,(c)傳統(tǒng)時間二階差分方法的模擬結果,時間步長 2.0 ms,(d)傳統(tǒng)時間二階差分方法的模擬結果,時間步長1.25 ms,(e)時間四階差分方法模擬的結果,時間步長5.0 ms,(f)時間六階差分方法模擬的結果,時間步長5.0 msFig.3 Acoustic snapshots at 3.0 s computed by (a) the traditional temporal 2nd-order FD scheme with the time step of 5.0 ms, (b) with the time step of 2.5 ms, (c) with the time step of 2.0 ms, and (d) with the time step of 1.25 ms, (e) by our temporal 4th-order and (f) 6th-order FD schemes with the time step of 5.0 ms

    圖4 BP鹽丘模型Fig.4 BP salt model

    圖5 傳統(tǒng)時間二階差分方法采用微小時間步長0.15 ms生成的參考地震記錄Fig.5 Reference seismic recording generated by traditional temporal 2nd-order FD scheme with a tiny time step of 0.15 ms

    圖6 (x,z)=(12,0) km處的單道地震記錄對比,時間段4.8~5.2 s,(a)時間二階、四階和六階差分方法的模擬結果,采用相同的時間步長1.5 ms,(b)時間二階差分方法采用0.75 ms,0.5 ms時間采樣間隔的模擬結果,以及時間四階、六階差分方法采用1.5 ms采樣間隔的模擬結果Fig.6 Recorded seismograms at (x,z)=(12,0) km during 4.8~5.2 s, (a) computed by temporal 2nd-order, 4th-order and 6th-order FD schemes with the same time step of 1.5 ms, (b) computed by temporal 2nd-order FD scheme with the time step of 0.75 ms, 0.5 ms, and computed by the temporal 4th-order and 6th-order FD schemes with the time step of 1.5 ms

    圖6中對比了遠偏移距接收點(x,z)=(12,0) km處的振動圖,時間區(qū)段為4.8~5.2 s.圖6a上、中、下三個子圖分別對比了傳統(tǒng)時間二階差分方法、本文發(fā)展的時間四階、六階差分方法模擬的地震道與參考道的波形差異,三種方法均采用1.5 ms的時間步長.由圖可知,傳統(tǒng)的時間二階差分方法模擬的結果與參考解存在明顯的差異,而本文發(fā)展的時間四階、六階差分方法與參考解吻合較好.圖6b上、中兩子圖顯示了傳統(tǒng)時間二階差分方法模擬的地震道,時間步長分別采用0.75 ms和0.5 ms,分別為圖6a中時間步長的1/2倍和1/3倍.圖6b中的第三個子圖則顯示了時間四階、六階差分方法模擬的地震道.仔細對比圖6b中的三個子圖,可發(fā)現采用1.5 ms時間步長的四階、六階差分方法模擬的結果幾乎一致,并且與參考解的吻合程度略好于采用0.75 ms時間步長的二階差分方法.然而,與理論的頻散分析結論有些差異的是,當傳統(tǒng)的時間二階差分方法采用的時間步長為0.5 ms時(見圖6b中第二個子圖),其吻合參考解的效果略微勝過采用1.5 ms時間步長的六階差分方法.我們也對比了其他地震道,發(fā)現了一致的規(guī)律.這說明模擬非均勻介質中的地震波場,時間四階精度搭配空間高階精度的差分方法已經足夠,而時間六階差分算子的優(yōu)勢并未體現.

    7結論

    本文基于中心網格剖分發(fā)展了適用于二階聲波方程的時間四階、六階,空間任意偶數階的有限差分數值模擬方法,并基于泰勒多項式展開推導了差分系數的解析表達式.與傳統(tǒng)的時間二階差分方法相比,本文發(fā)展的時間高階差分方法具有更好的穩(wěn)定性,允許更大的時間采樣步長.頻散分析表明本文發(fā)展的時間四階、六階精度的差分方法能顯著地改善傳統(tǒng)的時間二階精度差分方法的精度,在相同的頻散誤差限下,本文發(fā)展的時間四階、六階精度差分方法允許更大的時間步長,因此提高了波場數值模擬的計算效率.在相似的頻散條件下,理論的計算量分析表明本文發(fā)展的(16,4)和(16,6)差分算子的計算效率分別是傳統(tǒng)(16,2)方法的2.4倍和3.0倍.然而,非均勻介質中波場模擬的數值實驗證實,時間六階精度差分方法模擬的結果與時間四階精度差分方法模擬的結果幾乎完全一致,考慮到時間四階差分方法包含的浮點運算更少,因此,時間四階精度搭配空間高階精度更可取.此外,本文發(fā)展的時-空高階差分方法可以很容易地推廣至三維情況.

    References

    Alterman Z, Karal F. 1968. Propagation of elastic waves in layered media by finite difference methods.Bull.Seism.Soc.Am., 58(1): 367-398.Dong C H, Wang S X, Zhao J G, et al. 2013. Numerical experiment and analysis of the differential acoustic resonance spectroscopy for elastic property measurements.J.Geophys.Eng., 10(5): 1-10.

    Dong L G, Ma Z T, Cao J Z, et al. 2000. The staggered-grid high-order difference method of one-order elastic equation.ChineseJ.Geophys. (in Chinese), 43(6): 856-864.

    Du Q Z, Guo C F, Gong X F. 2015, Hybrid PS/FD numerical simulation and stability analysis of pure P-wave propagation in VTI media.ChineseJ.Geophys. (in Chinese), 58(4): 1290-1304.Fang G, Fomel S, Du Q Z, et al. 2014. Lowrank seismic-wave extrapolation on a staggered grid.Geophysics, 79(3): T157-T168.

    Finkelstein B, Kastner R. 2007. Finite difference time domain dispersion reduction schemes.J.Comput.Phys., 221(1): 422-438.Fomel S. Ying L X, Song X L. 2013. Seismic wave extrapolation using lowrank symbol approximation.Geophys.Prospect., 61(3): 526-536.

    Fornberg B. 1987. The pseudospectral method: comparisons with finite differences for the elastic wave equation.Geophysics, 52(4): 483-501.

    Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences.Bull.Seism.Soc.Am., 86(4): 1091-1106.Holberg O. 1987. Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena.Geophys.Prospect., 35(6): 629-655.

    Kelly K, Ward R, Treitel S, et al. 1976. Synthetic seismograms: A finite-difference approach.Geophysics, 41(1): 2-27.

    Kristek J, Moczo P, Galis M. 2010. Stable discontinuous staggered grid in the finite-difference modeling of seismic motion.Geophy.J.Int., 183(3): 1401-1407.Levander A R. 1988. Fourth-order finite-difference P-SV seismograms.Geophysics, 53(11): 1425-1436.Liang W Q, Wang Y F, Yang C C. 2015. Determination on the implicit finite-difference operator based on optimization method in time-space domain.Geophy.Prosp.Petroleum, 54(3): 254-259.

    Liu H F, Dai N X, Niu F, et al. 2014. An explicit time evolution method for acoustic wave propagation.Geophysics, 79(3): T117-T124.

    Liu Y, Li C C, Mu Y G. 1998. Finite-difference numerical modeling of any even-order accuracy.OilGeophys.Prosp., 33(1): 1-10.Moczo P, Kristek J, Halada L.2000. 3D fourth-order staggered-grid finite-difference schemes: Stability and grid dispersion.Bull.Seism.Soc.Am., 90(3): 587-603.

    Shi R Q, Wang S X, Zhao J G. 2012. An unsplit complex-frequency-shifted PML based on matched Z-transform for FDTD modelling of seismic wave equations.J.Geophys.Eng. 9(2): 218-229.

    Song X L, Fomel S, Ying L X. 2013. Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation.Geophy.J.Int., 193(2): 960-969.

    Tan S R, Huang L J. 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation.Geophy.J.Int., 197(2): 1250-1267.

    Virieux J. 1984. SH-wave propagation in heterogeneous media: Velocity-stress finite-difference method.Geophysics, 49(11): 1933-1957.

    Virieux J. 1986. P-SV wave propagation in heterogeneous media: Velocity stress finite-difference method.Geophysics, 51(4): 889-901.

    Wang Y F, Liang W, Nashed Z, et al. 2014. Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method.Geophysics, 79(5): T277-T285.

    Wu G C, Wang H Z. 2005. Analysis of numerical dispersion in wave-field simulation.ProgressinGeophysics. (in Chinese), 20(1), 58-65.Xie G S, Liu H, Zhao L G. 2005. Parallel Algorithm based on the multithread technique for pseudospectal modeling of seismic wave.ProgressinGeophysics. (in Chinese), 20(1), 17-23.

    Yang Q J, Liu C, Geng M X, et al. 2014. Staggered-grid finite-difference scheme and coefficients deduction of any number of derivatives.J.JilinUniv. (Earthscienceedition), 44(1): 375-385.

    Yee K S. 1966. Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media.IEEETrans.AntennasPropagat, 14(3): 302-307.

    Zhang J, Wu G C. 2015, Seismic modeling method of two phase medium utilizing staggered grid with optimal difference coefficients.ProgressinGeophysics. (in Chinese), 30(5), 2301-2311.

    附中文參考文獻

    董良國,馬在田,曹景中等. 2000. 一階彈性波方程交錯網格高階差分解法. 地球物理學報. 43(6): 856-864.

    杜啟振,郭成鋒,公緒飛. 2015. VTI介質純P波混合法正演模擬及穩(wěn)定性分析. 地球物理學報. 58(4): 1290-1304.

    梁文全,王彥飛,楊長春. 2015. 基于優(yōu)化方法的時間空間域隱格式有限差分算子確定方法. 石油物探. 54(3): 254-259.

    謝桂生,劉洪,趙連功. 2005. 偽譜法地震波正演模擬的多線程并行計算. 地球物理學進展. 20(1): 17-23.

    劉洋,李承楚,牟永光. 1998. 任意偶數階精度有限差分法數值模擬. 石油地球物理勘探. 33(1): 1-10.

    吳國忱,王華忠. 2005. 波場模擬中的數值頻散分析與校正策略. 地球物理學進展. 20(1): 58-65.

    楊慶節(jié),劉財,耿美霞等. 2014. 交錯網格任意階導數有限差分格式及差分系數推導. 吉林大學學報(地球科學版). 44(1): 375-385.

    張杰,吳國忱. 2015. 基于優(yōu)化差分系數的雙相介質交錯網格正演模擬. 地球物理學進展. 30(5): 2301-2311.

    (本文編輯劉少華)

    Time-space domain high-order finite-difference methods for seismic wave numerical simulation based on new stencils

    ZHANG Bao-Qing1,2,ZHOU Hui1*, CHEN Han-Ming1, SHENG Shan-Bo3

    1StateKeyLaboratoryofPetroleumResourcesandProspecting,CNPCKeyLabofGeophysicalExploration,ChinaUniversityofPetroleum,Beijing102249,China2DagangBranchofGRI,BGPInc.,CNPC,TianjinDagang300280,China3ChinaNationalOil&GasExplorationandDevelopmentCorporation,Beijing100034,China

    AbstractThe traditional finite-difference (FD) seismic wave simulation scheme adopts high-order FD operators to discretize the spatial derivatives, and the second-order FD operator to discretize the temporal derivative. Therefore, the traditional high-order FD method only achieves high-order accuracy in space, but exhibits low-order accuracy in time. When a relatively large time step is applied, the traditional FD method suffers from visible temporal dispersion and even instability. This paper develops new time-space domain high-order FD methods that attain arbitrary even-order accuracy in space, fourth- and sixth-order accuracy in time. The new FD methods are developed based on new FD stencils and a centered-grid. The FD coefficients are determined from the discrete dispersion relation using the Taylor-series expansion (TE) approach. Dispersion analysis indicates that our temporal fourth- and sixth-order FD methods improve the accuracy of the traditional temporal second-order FD method significantly. Computational cost analysis demonstrates that our temporal high-order FD methods are more efficient than the traditional temporal second-order method. Numerical simulation of seismic waves in homogeneous and heterogeneous media further validates the effectiveness of our high-order FD methods.KeywordsSeismic wave equation; Time-space domain; High-order finite-difference; Numerical simulation

    基金項目國家重點基礎研究發(fā)展計劃(973計劃)(2013CB228603),國家自然科學基金(41174119)和中石油物探新方法新技術研究(2014A-3609)聯合資助.

    作者簡介張保慶,男,1973年生,高級工程師,主要從事地震資料處理等工作. E-mail:zhangbaoqing1640@sina.com *通訊作者周輝,教授,博士生導師,主要從事地震勘探和地質雷達探測等地球物理資料的處理和正反演研究. E-mail: huizhou@cup.edu.cn

    doi:10.6038/cjg20160523 中圖分類號P631

    收稿日期2015-08-04,2016-04-22收修定稿

    張保慶, 周輝, 陳漢明等. 2016. 基于新的差分結構的時-空域高階有限差分波動方程數值模擬方法. 地球物理學報,59(5):1804-1814,doi:10.6038/cjg20160523.

    Zhang B Q, Zhou H, Chen H M, et al. 2016. Time-space domain high-order finite-difference methods for seismic wave numerical simulation based on new stencils.ChineseJ.Geophys. (in Chinese),59(5):1804-1814,doi:10.6038/cjg20160523.

    猜你喜歡
    數值模擬
    基于AMI的雙色注射成型模擬分析
    錐齒輪精密冷擺輾成形在“材料成型數值模擬”課程教學中的應用
    科教導刊(2016年28期)2016-12-12 06:22:00
    基于氣象信息及風場信息的風機輪轂處風速預測
    鉆孔灌注樁樁底沉渣對樁體承載特性影響的模擬分析
    西南地區(qū)氣象資料測試、預處理和加工研究報告
    科技資訊(2016年18期)2016-11-15 08:01:18
    張家灣煤礦巷道無支護條件下位移的數值模擬
    科技視界(2016年18期)2016-11-03 23:14:27
    張家灣煤礦開切眼錨桿支護參數確定的數值模擬
    科技視界(2016年18期)2016-11-03 22:57:21
    跨音速飛行中機翼水汽凝結的數值模擬研究
    科技視界(2016年18期)2016-11-03 20:38:17
    姚橋煤礦采空區(qū)CO2防滅火的數值模擬分析
    雙螺桿膨脹機的流場數值模擬研究
    科技視界(2016年22期)2016-10-18 14:53:19
    久久九九热精品免费| ponron亚洲| 99精品在免费线老司机午夜| 国产成+人综合+亚洲专区| 白带黄色成豆腐渣| 网址你懂的国产日韩在线| 国内精品一区二区在线观看| 免费无遮挡裸体视频| 亚洲国产高清在线一区二区三| 夜夜看夜夜爽夜夜摸| 91字幕亚洲| 中文字幕人成人乱码亚洲影| 九九久久精品国产亚洲av麻豆| avwww免费| 久久99热这里只有精品18| 日韩欧美三级三区| 精品国产美女av久久久久小说| 2021天堂中文幕一二区在线观| 精品一区二区三区视频在线观看免费| 又黄又爽又免费观看的视频| 18+在线观看网站| 免费观看精品视频网站| 最近最新中文字幕大全电影3| 特大巨黑吊av在线直播| 国产 一区 欧美 日韩| 精品久久久久久久久久免费视频| 观看免费一级毛片| 亚洲成a人片在线一区二区| 欧美丝袜亚洲另类 | 级片在线观看| av国产免费在线观看| 欧美最黄视频在线播放免费| 久久久国产成人免费| 久久久久久国产a免费观看| а√天堂www在线а√下载| 亚洲欧美日韩无卡精品| 天堂网av新在线| 亚洲精华国产精华精| 90打野战视频偷拍视频| 国产精品av视频在线免费观看| 在线观看免费视频日本深夜| 中文在线观看免费www的网站| 欧美成狂野欧美在线观看| 噜噜噜噜噜久久久久久91| 丁香欧美五月| 伊人久久精品亚洲午夜| 高清在线国产一区| 69人妻影院| 岛国在线免费视频观看| 免费大片18禁| 国产黄色小视频在线观看| 国产精品av视频在线免费观看| 日韩中文字幕欧美一区二区| 国产激情偷乱视频一区二区| 成人鲁丝片一二三区免费| 国产精品女同一区二区软件 | 九九热线精品视视频播放| 亚洲成人精品中文字幕电影| 老熟妇乱子伦视频在线观看| 国产又黄又爽又无遮挡在线| 毛片女人毛片| 综合色av麻豆| 国产真人三级小视频在线观看| 久久精品91蜜桃| 又黄又爽又免费观看的视频| 成人欧美大片| 日本三级黄在线观看| 桃红色精品国产亚洲av| 最后的刺客免费高清国语| 成人性生交大片免费视频hd| e午夜精品久久久久久久| 国产亚洲精品久久久com| 午夜福利视频1000在线观看| 男女做爰动态图高潮gif福利片| 欧美不卡视频在线免费观看| 99热这里只有是精品50| 少妇的丰满在线观看| 欧美成人a在线观看| 99精品欧美一区二区三区四区| 亚洲成人久久性| 香蕉久久夜色| 国产精品电影一区二区三区| 亚洲国产中文字幕在线视频| 少妇人妻精品综合一区二区 | 久久6这里有精品| 欧美bdsm另类| 人人妻人人澡欧美一区二区| 亚洲 欧美 日韩 在线 免费| 9191精品国产免费久久| 真实男女啪啪啪动态图| 可以在线观看的亚洲视频| 欧美xxxx黑人xx丫x性爽| 99久久无色码亚洲精品果冻| 免费电影在线观看免费观看| 在线看三级毛片| 九九热线精品视视频播放| 俺也久久电影网| 99久久综合精品五月天人人| 老鸭窝网址在线观看| 亚洲av美国av| 午夜福利欧美成人| 亚洲一区二区三区色噜噜| 99久久九九国产精品国产免费| 国产成人a区在线观看| 精品无人区乱码1区二区| 给我免费播放毛片高清在线观看| 久久久成人免费电影| 成人精品一区二区免费| 日韩av在线大香蕉| 国产爱豆传媒在线观看| aaaaa片日本免费| 亚洲欧美日韩无卡精品| 午夜免费观看网址| 精品久久久久久久久久免费视频| 精品乱码久久久久久99久播| 欧美最黄视频在线播放免费| 国产一区二区三区在线臀色熟女| av天堂中文字幕网| 给我免费播放毛片高清在线观看| 亚洲精品一区av在线观看| 午夜视频国产福利| 中文在线观看免费www的网站| 色尼玛亚洲综合影院| 性欧美人与动物交配| 99riav亚洲国产免费| 国产精品女同一区二区软件 | 国产在视频线在精品| av女优亚洲男人天堂| av视频在线观看入口| 亚洲精品美女久久久久99蜜臀| 欧美bdsm另类| 精品人妻偷拍中文字幕| 国产精品久久久久久久电影 | 宅男免费午夜| 午夜影院日韩av| 不卡一级毛片| 亚洲av五月六月丁香网| 看片在线看免费视频| 99久久成人亚洲精品观看| 两人在一起打扑克的视频| 国产淫片久久久久久久久 | 免费无遮挡裸体视频| 怎么达到女性高潮| 亚洲avbb在线观看| 国内精品久久久久久久电影| 免费看十八禁软件| 制服丝袜大香蕉在线| 亚洲国产欧洲综合997久久,| 亚洲精品美女久久久久99蜜臀| 国产精品av视频在线免费观看| 男人舔奶头视频| av福利片在线观看| 日韩欧美在线乱码| 久久国产精品人妻蜜桃| 老司机福利观看| 欧美日韩福利视频一区二区| 中国美女看黄片| or卡值多少钱| 少妇的逼好多水| 波多野结衣高清作品| 亚洲成人中文字幕在线播放| 成人国产一区最新在线观看| 又黄又粗又硬又大视频| 精品久久久久久久久久免费视频| 亚洲avbb在线观看| 国产日本99.免费观看| 欧美xxxx黑人xx丫x性爽| 亚洲欧美一区二区三区黑人| 成人特级av手机在线观看| 国产成人福利小说| 亚洲av成人av| 长腿黑丝高跟| 亚洲18禁久久av| 一进一出好大好爽视频| 国产精品,欧美在线| 黄色女人牲交| 久久精品国产亚洲av香蕉五月| 色噜噜av男人的天堂激情| 亚洲精品成人久久久久久| 国产一区二区在线观看日韩 | 亚洲激情在线av| 久久精品国产清高在天天线| 草草在线视频免费看| 18禁黄网站禁片免费观看直播| 日本与韩国留学比较| 伊人久久大香线蕉亚洲五| 国产高潮美女av| 好男人电影高清在线观看| 日韩欧美国产一区二区入口| 欧美又色又爽又黄视频| 午夜福利在线在线| 欧美又色又爽又黄视频| 最好的美女福利视频网| 91av网一区二区| 午夜福利欧美成人| 国产欧美日韩一区二区三| 国产私拍福利视频在线观看| 五月玫瑰六月丁香| 亚洲av日韩精品久久久久久密| 国产高清激情床上av| 亚洲成人精品中文字幕电影| www.999成人在线观看| 国产精品三级大全| 老汉色av国产亚洲站长工具| 国产一区二区激情短视频| 噜噜噜噜噜久久久久久91| 最近视频中文字幕2019在线8| 亚洲性夜色夜夜综合| 国产高清有码在线观看视频| 久久亚洲真实| 久久天躁狠狠躁夜夜2o2o| 在线观看午夜福利视频| 亚洲国产色片| 99热6这里只有精品| 免费看日本二区| 校园春色视频在线观看| 99久久成人亚洲精品观看| 欧美最黄视频在线播放免费| 午夜免费男女啪啪视频观看 | 国产在视频线在精品| 99久久成人亚洲精品观看| 日本一本二区三区精品| 国产伦精品一区二区三区四那| 久久精品国产综合久久久| 日本熟妇午夜| 国产亚洲精品久久久久久毛片| av在线蜜桃| 国产野战对白在线观看| 黄色女人牲交| 99精品欧美一区二区三区四区| 有码 亚洲区| 深爱激情五月婷婷| 亚洲最大成人手机在线| 成人性生交大片免费视频hd| e午夜精品久久久久久久| 日本三级黄在线观看| 国产成人a区在线观看| 免费在线观看日本一区| 国产熟女xx| 欧美激情久久久久久爽电影| 伊人久久精品亚洲午夜| xxx96com| 日本撒尿小便嘘嘘汇集6| 91在线精品国自产拍蜜月 | 亚洲欧美日韩卡通动漫| 在线免费观看的www视频| 伊人久久大香线蕉亚洲五| 国产亚洲精品一区二区www| 内地一区二区视频在线| 久久久色成人| 日韩欧美免费精品| 国产蜜桃级精品一区二区三区| ponron亚洲| 国产毛片a区久久久久| 黑人欧美特级aaaaaa片| 可以在线观看的亚洲视频| 成年免费大片在线观看| 久久伊人香网站| 丁香六月欧美| 色老头精品视频在线观看| 成年女人看的毛片在线观看| 一级a爱片免费观看的视频| 久久久成人免费电影| 一个人免费在线观看的高清视频| 一个人观看的视频www高清免费观看| 久久久色成人| 色av中文字幕| 少妇的逼水好多| 网址你懂的国产日韩在线| 国产高清视频在线播放一区| eeuss影院久久| 国产成人福利小说| 免费看十八禁软件| 性欧美人与动物交配| 高潮久久久久久久久久久不卡| 成人高潮视频无遮挡免费网站| 国产极品精品免费视频能看的| 日韩成人在线观看一区二区三区| 黄片小视频在线播放| 国产精品女同一区二区软件 | 国产精品一区二区三区四区免费观看 | 丰满人妻一区二区三区视频av | 日本一二三区视频观看| 久久久国产成人精品二区| 中文字幕久久专区| av天堂中文字幕网| 97碰自拍视频| 美女黄网站色视频| 免费人成视频x8x8入口观看| 免费在线观看成人毛片| 可以在线观看毛片的网站| 少妇高潮的动态图| 国产三级中文精品| 深夜精品福利| 在线观看美女被高潮喷水网站 | 真人做人爱边吃奶动态| 日本熟妇午夜| 国产高清视频在线播放一区| 国产一区在线观看成人免费| 十八禁网站免费在线| 亚洲欧美一区二区三区黑人| 国产欧美日韩一区二区精品| 欧美一级毛片孕妇| 国产伦精品一区二区三区四那| 日韩精品中文字幕看吧| 一个人看的www免费观看视频| 757午夜福利合集在线观看| 每晚都被弄得嗷嗷叫到高潮| 午夜福利18| 久久欧美精品欧美久久欧美| 亚洲美女视频黄频| 日韩 欧美 亚洲 中文字幕| 天堂√8在线中文| 99国产精品一区二区蜜桃av| 99热6这里只有精品| 欧美大码av| 久久精品91蜜桃| 久久精品夜夜夜夜夜久久蜜豆| 国产精品,欧美在线| 国产伦人伦偷精品视频| 可以在线观看的亚洲视频| 黄色丝袜av网址大全| 内射极品少妇av片p| 成年女人看的毛片在线观看| 国产真实伦视频高清在线观看 | 午夜影院日韩av| 99精品久久久久人妻精品| 中文资源天堂在线| 99热只有精品国产| 九色成人免费人妻av| 亚洲国产精品成人综合色| 久久精品国产99精品国产亚洲性色| 看免费av毛片| 夜夜看夜夜爽夜夜摸| 男人舔女人下体高潮全视频| 亚洲性夜色夜夜综合| tocl精华| 亚洲不卡免费看| 内射极品少妇av片p| 嫩草影视91久久| 国产极品精品免费视频能看的| 成人永久免费在线观看视频| 色哟哟哟哟哟哟| 高清日韩中文字幕在线| 岛国在线免费视频观看| 欧美黄色淫秽网站| 99久久99久久久精品蜜桃| 亚洲精品亚洲一区二区| 麻豆国产97在线/欧美| 国产亚洲精品久久久com| 日韩欧美国产在线观看| 久久久久国产精品人妻aⅴ院| av天堂在线播放| 琪琪午夜伦伦电影理论片6080| 午夜福利视频1000在线观看| 高清在线国产一区| 国产精品1区2区在线观看.| 村上凉子中文字幕在线| 午夜亚洲福利在线播放| 在线看三级毛片| 99久久精品一区二区三区| av在线蜜桃| 黄片小视频在线播放| 高清毛片免费观看视频网站| 欧美另类亚洲清纯唯美| 激情在线观看视频在线高清| 精华霜和精华液先用哪个| 成人午夜高清在线视频| 欧美黄色淫秽网站| 国产伦在线观看视频一区| 女人高潮潮喷娇喘18禁视频| 亚洲狠狠婷婷综合久久图片| 色视频www国产| 757午夜福利合集在线观看| 18禁黄网站禁片免费观看直播| 日日干狠狠操夜夜爽| 校园春色视频在线观看| 色老头精品视频在线观看| 国产一区二区在线观看日韩 | 一级作爱视频免费观看| 丁香六月欧美| 每晚都被弄得嗷嗷叫到高潮| 少妇的逼好多水| 国产麻豆成人av免费视频| 少妇的逼水好多| 欧美又色又爽又黄视频| 搡女人真爽免费视频火全软件 | 女生性感内裤真人,穿戴方法视频| 亚洲无线在线观看| 午夜福利在线观看吧| 国产淫片久久久久久久久 | 国产亚洲精品综合一区在线观看| 成人精品一区二区免费| 日本熟妇午夜| 亚洲 欧美 日韩 在线 免费| 欧美又色又爽又黄视频| 男女下面进入的视频免费午夜| 日韩av在线大香蕉| 久久伊人香网站| 岛国在线观看网站| 欧美在线黄色| 母亲3免费完整高清在线观看| 国产精品 欧美亚洲| 免费av不卡在线播放| 免费av毛片视频| 亚洲精品久久国产高清桃花| 天堂影院成人在线观看| 亚洲人成网站在线播放欧美日韩| 久久久久久久久大av| 一本久久中文字幕| 99热这里只有是精品50| 久久久久久人人人人人| 琪琪午夜伦伦电影理论片6080| aaaaa片日本免费| 欧美乱码精品一区二区三区| 美女大奶头视频| 免费人成视频x8x8入口观看| 高清在线国产一区| 精品一区二区三区视频在线观看免费| 久久久久久久久久黄片| 午夜精品久久久久久毛片777| av专区在线播放| 国产三级中文精品| 天美传媒精品一区二区| 色尼玛亚洲综合影院| 精品电影一区二区在线| 日本免费a在线| 国产成人啪精品午夜网站| av天堂中文字幕网| 哪里可以看免费的av片| 免费在线观看成人毛片| 757午夜福利合集在线观看| 精品国产美女av久久久久小说| 亚洲最大成人手机在线| 在线国产一区二区在线| 亚洲精品成人久久久久久| 日韩欧美精品免费久久 | 久久伊人香网站| 每晚都被弄得嗷嗷叫到高潮| 欧美日韩乱码在线| 级片在线观看| 男人舔奶头视频| 免费人成视频x8x8入口观看| 久久久精品大字幕| 嫩草影视91久久| 欧美在线一区亚洲| 国产综合懂色| 亚洲 国产 在线| 岛国在线免费视频观看| 精品人妻一区二区三区麻豆 | 午夜久久久久精精品| 我的老师免费观看完整版| 叶爱在线成人免费视频播放| 在线观看美女被高潮喷水网站 | or卡值多少钱| 欧美3d第一页| 日本黄色视频三级网站网址| 偷拍熟女少妇极品色| svipshipincom国产片| 国产一区二区三区视频了| 亚洲精品粉嫩美女一区| 首页视频小说图片口味搜索| 亚洲精品一卡2卡三卡4卡5卡| 深夜精品福利| 国产亚洲欧美98| 欧美最新免费一区二区三区 | 久久精品国产亚洲av涩爱 | 中出人妻视频一区二区| 欧美精品啪啪一区二区三区| 18禁国产床啪视频网站| 可以在线观看的亚洲视频| 日韩精品中文字幕看吧| 波多野结衣高清作品| 国产日本99.免费观看| av福利片在线观看| 丰满的人妻完整版| 亚洲人成电影免费在线| 又黄又爽又免费观看的视频| 国产成年人精品一区二区| xxxwww97欧美| 亚洲中文字幕日韩| 欧美性猛交╳xxx乱大交人| av黄色大香蕉| 午夜福利在线在线| 免费人成视频x8x8入口观看| 99热精品在线国产| 色噜噜av男人的天堂激情| 麻豆成人午夜福利视频| 波多野结衣巨乳人妻| 欧美bdsm另类| 国产淫片久久久久久久久 | 午夜精品久久久久久毛片777| 一本综合久久免费| 他把我摸到了高潮在线观看| 一级黄片播放器| 国产免费男女视频| 伊人久久大香线蕉亚洲五| 舔av片在线| 欧美3d第一页| 国产伦一二天堂av在线观看| 国产蜜桃级精品一区二区三区| 黄色片一级片一级黄色片| 国产精品国产高清国产av| 中文字幕人妻丝袜一区二区| 久久草成人影院| 1000部很黄的大片| 亚洲aⅴ乱码一区二区在线播放| 99精品在免费线老司机午夜| eeuss影院久久| 午夜福利视频1000在线观看| 久久精品国产综合久久久| 美女高潮喷水抽搐中文字幕| 女人被狂操c到高潮| 精品国产亚洲在线| 亚洲欧美一区二区三区黑人| 亚洲av成人精品一区久久| 99久久精品国产亚洲精品| 国产亚洲精品一区二区www| 国产成年人精品一区二区| 国产v大片淫在线免费观看| 国产精品亚洲一级av第二区| 叶爱在线成人免费视频播放| 国产探花在线观看一区二区| 最新中文字幕久久久久| 亚洲人成伊人成综合网2020| 精品不卡国产一区二区三区| 亚洲av一区综合| 免费观看精品视频网站| 男女做爰动态图高潮gif福利片| 女警被强在线播放| 午夜两性在线视频| 亚洲av第一区精品v没综合| 国产激情欧美一区二区| 国产免费av片在线观看野外av| 美女高潮的动态| 老熟妇仑乱视频hdxx| 桃红色精品国产亚洲av| 亚洲欧美日韩东京热| 色综合站精品国产| 国产欧美日韩精品一区二区| www.色视频.com| 色综合欧美亚洲国产小说| 日韩欧美国产一区二区入口| 青草久久国产| 99久久九九国产精品国产免费| 黄色视频,在线免费观看| 亚洲欧美日韩高清在线视频| www.色视频.com| 精品人妻一区二区三区麻豆 | 老司机福利观看| 精品人妻1区二区| 女警被强在线播放| 1000部很黄的大片| 99riav亚洲国产免费| 手机成人av网站| 午夜亚洲福利在线播放| 看黄色毛片网站| 国产黄a三级三级三级人| 真实男女啪啪啪动态图| 国产高清激情床上av| 精品日产1卡2卡| 成人av一区二区三区在线看| 欧美成狂野欧美在线观看| 婷婷亚洲欧美| 91在线精品国自产拍蜜月 | 亚洲熟妇中文字幕五十中出| 亚洲天堂国产精品一区在线| 搡女人真爽免费视频火全软件 | 69人妻影院| 欧美日本视频| 国产精品日韩av在线免费观看| 人妻丰满熟妇av一区二区三区| 国产精品,欧美在线| 日韩欧美在线二视频| 亚洲精品成人久久久久久| 女人被狂操c到高潮| 日本与韩国留学比较| 精品无人区乱码1区二区| 国产乱人视频| 欧美日韩一级在线毛片| 两个人的视频大全免费| tocl精华| 国产极品精品免费视频能看的| 在线观看午夜福利视频| 一个人看的www免费观看视频| 亚洲欧美一区二区三区黑人| 成人永久免费在线观看视频| 亚洲av电影在线进入| 9191精品国产免费久久| 亚洲狠狠婷婷综合久久图片| 国产欧美日韩精品一区二区| 一二三四社区在线视频社区8| 国产成人a区在线观看| 美女高潮喷水抽搐中文字幕| 观看免费一级毛片| 宅男免费午夜| 亚洲天堂国产精品一区在线| 免费av毛片视频| av黄色大香蕉| 免费电影在线观看免费观看| 中文在线观看免费www的网站| 国产精品久久久人人做人人爽| 亚洲精品乱码久久久v下载方式 | 午夜激情福利司机影院| 男人舔奶头视频| 午夜激情福利司机影院| 婷婷精品国产亚洲av在线| 久久久久性生活片| 欧美日本亚洲视频在线播放| 久久久国产成人精品二区| 亚洲精品日韩av片在线观看 | 亚洲成av人片在线播放无| 嫩草影院精品99| 一本综合久久免费| 国产伦人伦偷精品视频|