• <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
    99久久综合免费| 免费av中文字幕在线| 日韩大码丰满熟妇| 肉色欧美久久久久久久蜜桃| 女性生殖器流出的白浆| 美女高潮到喷水免费观看| 国产精品.久久久| 1024香蕉在线观看| 国产欧美日韩一区二区三区在线| 青青草视频在线视频观看| 十分钟在线观看高清视频www| 欧美人与性动交α欧美精品济南到| 欧美精品一区二区大全| 少妇被粗大的猛进出69影院| 精品亚洲成a人片在线观看| 亚洲欧美激情在线| 99热全是精品| 国产精品一区二区精品视频观看| 99久久99久久久精品蜜桃| 国产成人免费无遮挡视频| 午夜免费鲁丝| av国产精品久久久久影院| 50天的宝宝边吃奶边哭怎么回事| 永久免费av网站大全| 午夜福利乱码中文字幕| 国产精品一区二区精品视频观看| 18禁裸乳无遮挡动漫免费视频| 日韩三级视频一区二区三区| 欧美日韩黄片免| 啦啦啦中文免费视频观看日本| 久久久久国产一级毛片高清牌| 视频区图区小说| 国产一区二区三区av在线| 亚洲人成电影观看| 嫁个100分男人电影在线观看| 久久九九热精品免费| 亚洲精品一二三| 精品人妻在线不人妻| 日韩电影二区| avwww免费| 啦啦啦免费观看视频1| 欧美精品亚洲一区二区| 两性夫妻黄色片| 亚洲欧美一区二区三区黑人| 国产成人系列免费观看| 男人爽女人下面视频在线观看| 91大片在线观看| 色婷婷av一区二区三区视频| 亚洲国产欧美一区二区综合| 久久精品人人爽人人爽视色| av视频免费观看在线观看| 自线自在国产av| 少妇 在线观看| 手机成人av网站| 91精品伊人久久大香线蕉| √禁漫天堂资源中文www| 自拍欧美九色日韩亚洲蝌蚪91| 夜夜夜夜夜久久久久| 日本wwww免费看| 亚洲av欧美aⅴ国产| 免费高清在线观看日韩| 久久久水蜜桃国产精品网| 国产欧美亚洲国产| 国产男女超爽视频在线观看| 亚洲九九香蕉| 日本a在线网址| 天天操日日干夜夜撸| 精品国产乱码久久久久久男人| 黄片小视频在线播放| 国产野战对白在线观看| 纯流量卡能插随身wifi吗| 亚洲欧美精品综合一区二区三区| av免费在线观看网站| 欧美性长视频在线观看| 高清欧美精品videossex| 超碰成人久久| 日韩 亚洲 欧美在线| 欧美精品亚洲一区二区| 91九色精品人成在线观看| 久久久久精品国产欧美久久久 | 老鸭窝网址在线观看| 国产黄色免费在线视频| 久久精品久久久久久噜噜老黄| 视频区欧美日本亚洲| 午夜福利视频在线观看免费| 我要看黄色一级片免费的| 成人三级做爰电影| 欧美亚洲 丝袜 人妻 在线| 男女免费视频国产| 精品一品国产午夜福利视频| 新久久久久国产一级毛片| 亚洲全国av大片| 婷婷色av中文字幕| 亚洲av成人不卡在线观看播放网 | 久久精品成人免费网站| 久久精品亚洲av国产电影网| 免费观看a级毛片全部| 精品国产一区二区久久| 欧美日韩亚洲国产一区二区在线观看 | av在线老鸭窝| 国产视频一区二区在线看| 天天躁日日躁夜夜躁夜夜| 中国国产av一级| 老司机午夜福利在线观看视频 | 国产91精品成人一区二区三区 | 亚洲少妇的诱惑av| 亚洲精品国产区一区二| 亚洲精品第二区| 建设人人有责人人尽责人人享有的| 免费在线观看视频国产中文字幕亚洲 | 老司机在亚洲福利影院| 色老头精品视频在线观看| 久久久久久久精品精品| 多毛熟女@视频| svipshipincom国产片| 黑人猛操日本美女一级片| 国产亚洲欧美精品永久| 搡老乐熟女国产| 亚洲av成人一区二区三| a级毛片黄视频| 啦啦啦啦在线视频资源| 精品国产超薄肉色丝袜足j| 伦理电影免费视频| 精品一区二区三区av网在线观看 | 欧美午夜高清在线| 亚洲精品在线美女| 人妻久久中文字幕网| www.熟女人妻精品国产| 国产精品一区二区在线不卡| 欧美国产精品va在线观看不卡| 天天躁夜夜躁狠狠躁躁| 久久亚洲国产成人精品v| 一级片免费观看大全| 国产男人的电影天堂91| 国产1区2区3区精品| av免费在线观看网站| tube8黄色片| 国产一区二区激情短视频 | 啦啦啦视频在线资源免费观看| 午夜激情av网站| 国产一区有黄有色的免费视频| 色老头精品视频在线观看| 亚洲全国av大片| 1024香蕉在线观看| 老鸭窝网址在线观看| 国产97色在线日韩免费| 精品卡一卡二卡四卡免费| 天堂俺去俺来也www色官网| 国产真人三级小视频在线观看| 91成年电影在线观看| 亚洲激情五月婷婷啪啪| 国产成人免费观看mmmm| 黄片播放在线免费| 亚洲精品成人av观看孕妇| 精品国产超薄肉色丝袜足j| 午夜福利视频在线观看免费| 久久久水蜜桃国产精品网| 亚洲av欧美aⅴ国产| 欧美日韩福利视频一区二区| 亚洲综合色网址| 黑人猛操日本美女一级片| 精品久久久精品久久久| 久久久久久久国产电影| 成年动漫av网址| 国产在线免费精品| 免费看十八禁软件| 人妻久久中文字幕网| 色播在线永久视频| 亚洲成人免费电影在线观看| 成年av动漫网址| 久久久精品免费免费高清| 老司机影院毛片| 18在线观看网站| 一区二区av电影网| 成人影院久久| 老司机在亚洲福利影院| 亚洲国产毛片av蜜桃av| 捣出白浆h1v1| 国产精品成人在线| 狂野欧美激情性xxxx| 美女福利国产在线| 亚洲va日本ⅴa欧美va伊人久久 | 国产精品久久久av美女十八| 宅男免费午夜| 精品国产一区二区三区久久久樱花| 亚洲黑人精品在线| 99精品久久久久人妻精品| 2018国产大陆天天弄谢| 久久国产亚洲av麻豆专区| 欧美精品亚洲一区二区| 久久久国产欧美日韩av| 在线天堂中文资源库| 久久天躁狠狠躁夜夜2o2o| 亚洲黑人精品在线| 久久人妻熟女aⅴ| 国产av精品麻豆| 亚洲av日韩在线播放| 亚洲中文av在线| 成年人黄色毛片网站| 国产精品 国内视频| 一本色道久久久久久精品综合| www.精华液| a 毛片基地| 亚洲少妇的诱惑av| 岛国毛片在线播放| 老鸭窝网址在线观看| 亚洲伊人久久精品综合| 啦啦啦中文免费视频观看日本| a级毛片黄视频| 99久久人妻综合| 国产欧美日韩精品亚洲av| 国产精品久久久久久人妻精品电影 | 91成人精品电影| 午夜福利,免费看| 国产成人精品无人区| 91麻豆av在线| 精品国产国语对白av| 97在线人人人人妻| 视频区欧美日本亚洲| 老司机午夜十八禁免费视频| netflix在线观看网站| 99国产极品粉嫩在线观看| 中文字幕另类日韩欧美亚洲嫩草| av网站在线播放免费| 手机成人av网站| a级片在线免费高清观看视频| 免费高清在线观看视频在线观看| 91老司机精品| 精品少妇一区二区三区视频日本电影| 97在线人人人人妻| 国产欧美亚洲国产| 51午夜福利影视在线观看| 日韩一卡2卡3卡4卡2021年| 国产精品一区二区在线不卡| 人妻人人澡人人爽人人| 男人舔女人的私密视频| 国产免费av片在线观看野外av| 视频区欧美日本亚洲| 99国产精品一区二区蜜桃av | 我的亚洲天堂| 人妻人人澡人人爽人人| 一级黄色大片毛片| 国产又爽黄色视频| 一个人免费在线观看的高清视频 | 国产成人av激情在线播放| 亚洲欧美日韩另类电影网站| 精品亚洲成a人片在线观看| 啦啦啦 在线观看视频| 男女下面插进去视频免费观看| 99香蕉大伊视频| 久久影院123| 亚洲一区中文字幕在线| 中文字幕制服av| 国产一区二区激情短视频 | 亚洲成人手机| 丝袜在线中文字幕| 人妻一区二区av| 午夜福利在线观看吧| av网站在线播放免费| 一本综合久久免费| 91成人精品电影| 亚洲av美国av| 亚洲欧美一区二区三区黑人| 热99re8久久精品国产| 两个人看的免费小视频| 欧美日韩av久久| 丝袜喷水一区| 日韩视频一区二区在线观看| 久久性视频一级片| 国产精品亚洲av一区麻豆| 99久久综合免费| 国产成人欧美| 中文字幕人妻丝袜一区二区| 这个男人来自地球电影免费观看| 欧美中文综合在线视频| 国产亚洲精品一区二区www | 一级片免费观看大全| 中文字幕色久视频| 欧美日韩国产mv在线观看视频| 91字幕亚洲| 久久av网站| 国产在线免费精品| 最新的欧美精品一区二区| 国产精品成人在线| 日韩制服骚丝袜av| 人妻久久中文字幕网| 久久国产精品人妻蜜桃| 一本综合久久免费| 手机成人av网站| 老熟女久久久| 亚洲精品粉嫩美女一区| 国产av又大| 一级毛片女人18水好多| 国产无遮挡羞羞视频在线观看| 精品乱码久久久久久99久播| 免费高清在线观看视频在线观看| 黑人猛操日本美女一级片| 老熟妇乱子伦视频在线观看 | 最近中文字幕2019免费版| 别揉我奶头~嗯~啊~动态视频 | 精品亚洲乱码少妇综合久久| 久久亚洲精品不卡| 久久人妻熟女aⅴ| 久久久久视频综合| 亚洲av片天天在线观看| 女性被躁到高潮视频| 欧美+亚洲+日韩+国产| 18禁国产床啪视频网站| 性色av乱码一区二区三区2| 亚洲人成电影免费在线| 成人影院久久| www日本在线高清视频| www.自偷自拍.com| 亚洲av美国av| 视频在线观看一区二区三区| 少妇的丰满在线观看| 午夜日韩欧美国产| 极品少妇高潮喷水抽搐| 色94色欧美一区二区| 欧美日韩一级在线毛片| 91老司机精品| av免费在线观看网站| 搡老岳熟女国产| 国产av精品麻豆| 一本—道久久a久久精品蜜桃钙片| 91九色精品人成在线观看| 精品一区二区三卡| 嫁个100分男人电影在线观看| 女人爽到高潮嗷嗷叫在线视频| 欧美97在线视频| 如日韩欧美国产精品一区二区三区| 精品国内亚洲2022精品成人 | 热99久久久久精品小说推荐| 亚洲av国产av综合av卡| 丰满人妻熟妇乱又伦精品不卡| 老司机在亚洲福利影院| 丁香六月天网| 欧美xxⅹ黑人| 狂野欧美激情性xxxx| 免费日韩欧美在线观看| 少妇精品久久久久久久| 啪啪无遮挡十八禁网站| av天堂久久9| 建设人人有责人人尽责人人享有的| 宅男免费午夜| 香蕉国产在线看| 久久青草综合色| 亚洲少妇的诱惑av| 久久国产精品大桥未久av| 狂野欧美激情性xxxx| 97人妻天天添夜夜摸| 少妇被粗大的猛进出69影院| 国产一区二区三区av在线| 一级毛片精品| 国产在线一区二区三区精| 日本a在线网址| 多毛熟女@视频| av电影中文网址| 国产真人三级小视频在线观看| 一级片免费观看大全| 国产日韩一区二区三区精品不卡| 免费在线观看影片大全网站| 青春草视频在线免费观看| 免费观看人在逋| 别揉我奶头~嗯~啊~动态视频 | 老汉色∧v一级毛片| av欧美777| 精品国产乱码久久久久久男人| 波多野结衣一区麻豆| 男女高潮啪啪啪动态图| 久久久精品区二区三区| 男女床上黄色一级片免费看| 久久精品亚洲av国产电影网| 窝窝影院91人妻| 十八禁高潮呻吟视频| 狠狠狠狠99中文字幕| 正在播放国产对白刺激| 欧美日韩亚洲高清精品| 国产一区二区三区在线臀色熟女 | 国产成人免费观看mmmm| 搡老熟女国产l中国老女人| 国产av一区二区精品久久| 成年人黄色毛片网站| 久久综合国产亚洲精品| 韩国高清视频一区二区三区| 精品人妻熟女毛片av久久网站| 性少妇av在线| 亚洲自偷自拍图片 自拍| 波多野结衣av一区二区av| 免费观看人在逋| 婷婷丁香在线五月| 美女国产高潮福利片在线看| 男人添女人高潮全过程视频| 久久免费观看电影| 制服人妻中文乱码| 国产精品久久久久久精品古装| 美女高潮到喷水免费观看| 性少妇av在线| 国产精品熟女久久久久浪| 久久av网站| 可以免费在线观看a视频的电影网站| 亚洲国产精品成人久久小说| 一本一本久久a久久精品综合妖精| 亚洲av男天堂| 欧美成人午夜精品| 久久女婷五月综合色啪小说| 国产成人啪精品午夜网站| 久久久久久人人人人人| 99国产精品一区二区蜜桃av | 国产精品久久久av美女十八| 中文字幕另类日韩欧美亚洲嫩草| 亚洲国产中文字幕在线视频| 老汉色∧v一级毛片| 日韩视频一区二区在线观看| 黑人猛操日本美女一级片| 久久久久久久国产电影| 人成视频在线观看免费观看| 国产99久久九九免费精品| 好男人电影高清在线观看| 天天操日日干夜夜撸| 我的亚洲天堂| 午夜影院在线不卡| 国产av又大| av不卡在线播放| 老汉色av国产亚洲站长工具| 国产成人一区二区三区免费视频网站| 91精品伊人久久大香线蕉| 一区二区日韩欧美中文字幕| 女警被强在线播放| 国产精品一区二区精品视频观看| 亚洲精品粉嫩美女一区| 免费久久久久久久精品成人欧美视频| 飞空精品影院首页| 永久免费av网站大全| 国产精品欧美亚洲77777| 啦啦啦啦在线视频资源| 国产三级黄色录像| www.精华液| 亚洲 国产 在线| 大陆偷拍与自拍| 搡老乐熟女国产| 女人高潮潮喷娇喘18禁视频| 国产精品久久久久久人妻精品电影 | 国产真人三级小视频在线观看| tocl精华| www日本在线高清视频| 日韩,欧美,国产一区二区三区| 久久久国产一区二区| 日本五十路高清| 最近最新免费中文字幕在线| 色播在线永久视频| 国产精品久久久av美女十八| 日本av免费视频播放| 天天躁日日躁夜夜躁夜夜| 久久毛片免费看一区二区三区| 久久精品久久久久久噜噜老黄| 国产精品99久久99久久久不卡| 久久久久精品人妻al黑| 可以免费在线观看a视频的电影网站| 天堂俺去俺来也www色官网| 国产亚洲精品久久久久5区| 高潮久久久久久久久久久不卡| 午夜福利视频在线观看免费| 男人操女人黄网站| 欧美激情久久久久久爽电影 | 中文字幕最新亚洲高清| 一级黄色大片毛片| 日本精品一区二区三区蜜桃| 中文精品一卡2卡3卡4更新| 一二三四在线观看免费中文在| 国产伦理片在线播放av一区| 99久久国产精品久久久| 国产视频一区二区在线看| 满18在线观看网站| 人人澡人人妻人| 国产精品国产av在线观看| 最新在线观看一区二区三区| 手机成人av网站| 99精国产麻豆久久婷婷| 精品国产乱子伦一区二区三区 | 精品国产一区二区久久| 麻豆乱淫一区二区| 在线看a的网站| 一本—道久久a久久精品蜜桃钙片| 热re99久久精品国产66热6| 成年美女黄网站色视频大全免费| 亚洲色图综合在线观看| 国产精品久久久人人做人人爽| 少妇人妻久久综合中文| 12—13女人毛片做爰片一| 啪啪无遮挡十八禁网站| 久久精品国产综合久久久| 久久女婷五月综合色啪小说| 两个人看的免费小视频| 国产不卡av网站在线观看| 老司机午夜福利在线观看视频 | 91字幕亚洲| 亚洲av欧美aⅴ国产| 男男h啪啪无遮挡| 在线永久观看黄色视频| 性少妇av在线| 欧美日韩黄片免| 亚洲av欧美aⅴ国产| 丝袜美足系列| 国产伦理片在线播放av一区| 多毛熟女@视频| 欧美av亚洲av综合av国产av| 中文字幕人妻熟女乱码| 精品少妇一区二区三区视频日本电影| 最近中文字幕2019免费版| 亚洲精品久久午夜乱码| 男人添女人高潮全过程视频| 精品久久久久久电影网| 女人被躁到高潮嗷嗷叫费观| 色精品久久人妻99蜜桃| 99国产精品一区二区三区| 日韩熟女老妇一区二区性免费视频| 午夜成年电影在线免费观看| 在线观看人妻少妇| 在线观看免费午夜福利视频| 国产无遮挡羞羞视频在线观看| 亚洲欧美一区二区三区黑人| 欧美日韩国产mv在线观看视频| 美女高潮到喷水免费观看| 高潮久久久久久久久久久不卡| 国产不卡av网站在线观看| 国产精品欧美亚洲77777| 精品一区二区三卡| 中文字幕色久视频| 国产亚洲av片在线观看秒播厂| 日本撒尿小便嘘嘘汇集6| 成年av动漫网址| 99久久精品国产亚洲精品| 女人高潮潮喷娇喘18禁视频| 国产成人免费观看mmmm| 黑人巨大精品欧美一区二区蜜桃| 久久香蕉激情| 又紧又爽又黄一区二区| 交换朋友夫妻互换小说| 99九九在线精品视频| 日韩三级视频一区二区三区| 正在播放国产对白刺激| 美女中出高潮动态图| 首页视频小说图片口味搜索| 99久久人妻综合| 精品熟女少妇八av免费久了| 999精品在线视频| 亚洲人成电影观看| 狂野欧美激情性bbbbbb| 在线看a的网站| 老司机深夜福利视频在线观看 | 一级毛片电影观看| 成年动漫av网址| 老司机福利观看| 精品国产一区二区三区四区第35| 少妇的丰满在线观看| 久久久国产欧美日韩av| 国产成人免费观看mmmm| 十八禁高潮呻吟视频| 日韩制服骚丝袜av| 午夜免费成人在线视频| 18禁黄网站禁片午夜丰满| 18禁观看日本| 日韩有码中文字幕| 国产欧美日韩一区二区三区在线| av福利片在线| 狠狠精品人妻久久久久久综合| 三级毛片av免费| 老熟妇仑乱视频hdxx| 黄片大片在线免费观看| 欧美日韩精品网址| 人成视频在线观看免费观看| 国产成+人综合+亚洲专区| 性高湖久久久久久久久免费观看| xxxhd国产人妻xxx| 国产不卡av网站在线观看| 欧美另类亚洲清纯唯美| 在线观看一区二区三区激情| 十分钟在线观看高清视频www| 成人国产一区最新在线观看| tocl精华| 一级黄色大片毛片| 最近最新免费中文字幕在线| 黄色视频在线播放观看不卡| 免费观看av网站的网址| 男男h啪啪无遮挡| 国产精品一区二区在线观看99| 王馨瑶露胸无遮挡在线观看| 最近中文字幕2019免费版| 欧美日韩国产mv在线观看视频| 中文字幕人妻丝袜制服| 热99国产精品久久久久久7| 99国产精品一区二区蜜桃av | 亚洲欧美一区二区三区黑人| 色婷婷久久久亚洲欧美| 久久久久久久久久久久大奶| 日韩 亚洲 欧美在线| 大香蕉久久成人网| 欧美日韩亚洲高清精品| 日韩有码中文字幕| 国产野战对白在线观看| 午夜激情久久久久久久| 国产精品国产av在线观看| 国产日韩一区二区三区精品不卡| 亚洲精品美女久久久久99蜜臀| 窝窝影院91人妻| 亚洲性夜色夜夜综合| www日本在线高清视频| 国产成人精品久久二区二区91| 色老头精品视频在线观看| 亚洲av欧美aⅴ国产| 亚洲欧美精品综合一区二区三区| 丝袜人妻中文字幕| 久久久久久免费高清国产稀缺| 亚洲一卡2卡3卡4卡5卡精品中文|