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

    全火星模型中震波傳播特征的數(shù)值模擬研究

    2020-07-30 10:35:12鄧迪肖萬博王彥賓
    北京大學學報(自然科學版) 2020年4期
    關(guān)鍵詞:面波波場震源

    鄧迪 肖萬博 王彥賓

    全火星模型中震波傳播特征的數(shù)值模擬研究

    鄧迪 肖萬博 王彥賓?

    北京大學地球與空間科學學院地球物理學系, 北京 100871; ?通信作者, E-mail: ybwang@pku.edu.cn

    采用基于交錯網(wǎng)格的傅里葉偽譜與有限差分混合方法, 求解彈性波動方程, 根據(jù)地球化學分析得到的兩個火星理論結(jié)構(gòu)模型, 模擬二維全火星模型中 P-SV 波和 SH 波的傳播過程。根據(jù)理論地震圖和波場快照, 討論全火星模型中震波的傳播過程以及各種震相的產(chǎn)生和演變, 分析模型內(nèi)部火星殼厚度以及火星核幔邊界深度對震波傳播的影響。結(jié)果表明, 在低速火星殼內(nèi)部多重反射波及轉(zhuǎn)換波的相干疊加會形成很強的波列, 其特征受火星殼厚度的影響較大, 在切向分量上可以更清晰地觀測到核幔邊界的反射震相。

    全火星模型; 震波傳播; 數(shù)值模擬; 偽譜與有限差分混合方法

    火星是太陽系中與地球有最多相似性的類地行星。研究火星的內(nèi)部結(jié)構(gòu)和物質(zhì)組成, 對認識火星的形成和演化, 進而研究地球乃至太陽系的形成和演化歷程具有非常重要的科學意義[1], 因此這也成為空間探測的重要研究課題。

    地震波穿透性強, 對速度界面敏感, 具有較高的分辨能力, 普遍認為地震學是確定行星內(nèi)部結(jié)構(gòu)最好的地球物理學工具, 在地球內(nèi)部結(jié)構(gòu)研究中發(fā)揮著非常重要的作用。在美國阿波羅登月計劃期間, 航天員在月球表面布設(shè) 4 臺月震儀, 記錄了大量月震數(shù)據(jù)。通過對這些數(shù)據(jù)的分析, 揭示了月球的內(nèi)部結(jié)構(gòu)特征[2]。

    火星地震學的研究始于 1976 年美國的海盜計劃, 但遠不如月震探測成功。海盜 1 號著陸器于1976 年 7 月 20 日降落在火星上, 但是著陸器上的地震儀無法解鎖, 沒有任何數(shù)據(jù)返回。海盜 2 號著陸器于 1976 年 9 月 3 日登陸火星上的烏托邦平原, 著陸器上的地震儀連續(xù)工作 19 個月, 但除在第 80 個火星日記錄到一個可能的火星震事件外, 沒有其他發(fā)現(xiàn)[3]。該事件發(fā)生時, 沒有記錄到風數(shù)據(jù), 因此不能排除風噪聲的干擾。沒有記錄到其他事件的原因可能是地震儀在遠震體波的頻帶寬度中靈敏性不足[4],也可能是由于儀器的位置導(dǎo)致其對風噪聲的高靈敏度[5?6]。

    隕石撞擊是引起行星內(nèi)部震波傳播的另一個潛在來源。在阿波羅計劃 1969—1977 年對月球的 8年觀測中, 4 個臺站觀測到 13000 多個月震事件, 其中包含約 1800 多個隕石撞擊事件[7]。在月球上, 隕石直接落在月球表面并產(chǎn)生月震信號, 由于缺乏大氣層, 因此僅通過表面位移就可以檢測到所有的撞擊。但是, 火星的大氣層能非常有效地防止隕石撞擊, 到達火星表面的 10kg 以下隕石的數(shù)量僅為入射到大氣層隕石數(shù)量的 10%, 因此在火星上檢測到隕石撞擊的數(shù)量很少[3]。

    對火星核結(jié)構(gòu)的研究目前主要通過重力觀測進行。Yoder 等[8]通過火星全球探勘者號(Mars Global Surveyor)的無線電跟蹤分析, 測得火星太陽引力下的潮汐變形勒夫數(shù)2。觀察得到的2值為 0.153± 0.017, 排除了火星核為固體鐵芯的可能性, 至少外核是液體。之后, 研究人員對 2001 年火星奧德賽號(Mars Odyssey)、2003 年機遇號火星漫游車(Mars Exploration Rover Opportunity)和 2005 年火星勘測軌道飛行器計劃(Mars Reconnaissance Orbiter)發(fā)回的跟蹤數(shù)據(jù)進行一系列的分析, 確定2值為 0.169± 0.006[9?10]。

    Sohl 等[11]提出火星內(nèi)部結(jié)構(gòu)的兩個模型, 他們在流體靜力平衡和封閉傳熱的條件下進行求解, 得到質(zhì)量、靜水壓力、重力、溫度和熱流密度的徑向分布, 還獲得火星模型的分層密度和火星震波速結(jié)構(gòu)。兩個模型分別是滿足地球物理約束(即極慣性矩達到最大可能值)的模型 A 以及與來自火星的SNC 隕石的地球化學約束一致的模型 B。

    2018 年 11 月 26 日, 洞察號火星地球物理探測器降落在埃律西昂平原, 并在火星表面部署地震實驗儀器(Seismic Experiment for Internal Structure, SEIS), 用于研究火星內(nèi)部結(jié)構(gòu)。地震儀由長周期三軸寬頻帶儀和三軸短周期儀組成, 覆蓋 0.01~50Hz 的頻率范圍。與海盜號地震儀相比, SEIS 檢測火星震的分辨率有所提升。另一項與海盜號不同的重大改進是, 地震儀通過機械臂直接部署在火星表面, 并通過高效的隔熱和隔風保護, 使其免受溫度和風變化的影響。基于現(xiàn)有的火星知識, 預(yù)計 SEIS 每年可檢測到幾十次震中距 40°以內(nèi)、矩震級大于 3級的火星震[12]。2019 年 4 月 23日, 法國國家空間研究中心(Centre National d’Etudes Spatiales, CNES)宣布洞察號首次在火星上檢測到發(fā)生于 2019 年 4月 6 日的火星震信號。目前, 科學家仍在研究該火星震的數(shù)據(jù), 以期確認其震源信息[13]。在洞察號發(fā)射之前, 針對可能觀測到的火星震數(shù)據(jù), 人們利用地震學方法開展了相關(guān)研究, 如火星震的活動性、火星震定位、震波傳播、火星淺層結(jié)構(gòu)和殼幔結(jié)構(gòu)等, 為洞察號的設(shè)計與發(fā)射以及火星震觀測數(shù)據(jù)的分析奠定了基礎(chǔ)[14?20]。

    地震波數(shù)值模擬研究在地震學中起到非常重要的作用, 有助于理解復(fù)雜地震波的傳播過程和解釋地震觀測數(shù)據(jù)。Furumura 等[21]將傅里葉偽譜法與有限差分法相結(jié)合, 對 1999 年臺灣集集地震進行三維強震地面運動的數(shù)值模擬。馬德堂等[22]運用傅里葉偽譜法與有限元法的混合方法, 對彈性波動方程進行求解。魏星等[23]運用傅里葉偽譜法與四階精度有限差分方法的混合方法, 基于傅里葉偽譜法精度高、內(nèi)存消耗少的前提, 發(fā)揮四階精度有限差分方法對人工邊界和自由界面的處理優(yōu)勢, 對二維模型進行彈性波場的數(shù)值模擬計算。秦艷芳等[24]將魏星等的方法推廣到三維非均勻介質(zhì)地震波傳播模擬計算中, 在確保精度和效率的情況下, 成功地實現(xiàn)對三維沉積盆地模型的并行模擬計算。Wang等[25]和 Jiang 等[26]應(yīng)用交錯網(wǎng)格傅里葉偽譜法與有限差分法的混合方法求解彈性動力學方程, 分別模擬計算全月球模型中 P-SV 波和 SH 波的傳播。

    本文對二維全火星模型中火星震波傳播進行數(shù)值模擬, 針對已有的理論火星結(jié)構(gòu)模型, 模擬火星震波傳播的理論地震圖和波場快照, 討論火星殼厚度和核幔邊界深度對火星震波傳播的影響, 可為解釋未來可能獲取的火星震數(shù)據(jù)以及火星結(jié)構(gòu)和火星震發(fā)生機制的研究提供方法參考。

    1 基于交錯網(wǎng)格的傅里葉偽譜與有限差分混合方法和全火星模型

    對于各向同性的彈性介質(zhì), 在柱坐標系(,,)下, 假設(shè)方向上所有變量保持不變, 在(,)平面內(nèi), 以速度和應(yīng)力形式表示的 P-SV 波和 SH 波的二維波動方程為

    對于各向同性的彈性介質(zhì), 其本構(gòu)關(guān)系為

    其中,v,vv分別為垂向、徑向和切向的速度分量,為介質(zhì)密度,f,ff分別為垂向、徑向和切向的體力分量,,,,為各應(yīng)力張量的分量,和為拉梅常數(shù)。

    本文中網(wǎng)格劃分采用如圖 1 所示的交錯網(wǎng)格, 應(yīng)力和速度分量在不同位置離散化, 彼此之間有半個網(wǎng)格間距。半徑方向采用等間距網(wǎng)格離散化, 對于每個固定的半徑, 橫向網(wǎng)格點的數(shù)量是相同的, 橫向網(wǎng)格的間距隨著深度的增加而逐漸減小。采用基于交錯網(wǎng)格的有限差分與傅里葉偽譜的混合方法, 對波動方程組進行求解。

    在橫向上, 通過傅里葉偽譜法計算空間變量對的偏導(dǎo)數(shù)[27]:

    其中,=1, 2, …,–1,()表示在方向上離散的空間變量,()表示其相應(yīng)的波數(shù)域的傅里葉變換,表示空間離散的網(wǎng)格節(jié)點數(shù),和分別表示在方向上和波數(shù)域的離散間隔。

    在半徑方向上, 空間變量()對的偏導(dǎo)數(shù)由交錯網(wǎng)格的四階精度有限差分格式計算[28]:

    在模型區(qū)域內(nèi)共需要處理兩個邊界條件, 通過滿足表面處牽引力為零的自由邊界條件, 將自由表面引入模型中。對于模型的中心區(qū)域, 采用 Cerjan等[29]提出的帶有 20 個網(wǎng)格點的“緩沖區(qū)”的吸收邊界條件。式(1)中的震源體力采用 Wang 等[30?31]給出的二維柱坐標的線源形式:

    其中, (0,0)為震源中心位置;M,,,M分別為震源矩張量的分量;()和()表示震源體力的空間分布方式, 采用 Herrmann[32]提出的偽函數(shù)進行近似計算。

    本文選取以往研究中廣泛采用的 Sohl 等[11]提出的火星內(nèi)部結(jié)構(gòu)參考模型 A 和模型 B, 兩個模型的火星殼厚度分別為 110 和 250km, 火星幔上層與火星幔下層由橄欖石?尖晶石過渡區(qū)分隔開, 火星幔下層又分為尖晶石層和尖晶石層兩層。因此, 火星內(nèi)部分為 5 層, 由淺至深分別是火星殼、火星幔上層(橄欖石層)、尖晶石層、尖晶石層和火星核, 各層的 P 波速度、S 波速度和密度分布如圖 2所示。本文的值模型采用 Toyokuni 等[33]給出的值作為參考, 即液體核心內(nèi)的值為 1000000, 其他區(qū)域的值均為600。

    實線代表模型A, 虛線代表模型B

    2 全火星模型的數(shù)值模擬

    由于火星震的震源機制沒有任何觀測資料的約束, 本文在模擬過程中將典型的剪切地震震源作為參考。本研究的震源機制為二維雙力偶線源, 取M= –1.0,M=1.0, M=M=0, 相當于一個 45°傾角的傾滑斷層。模型的方向離散為 2048 個網(wǎng)格點, 范圍為 0~2π;方向離散為 512 個網(wǎng)格點, 范圍為0~3390km。因此,方向的網(wǎng)格間距為 6.62km,方向的網(wǎng)格間距最小為 0.02km, 最大為 10.40km, 位于火星模型表面。根據(jù) Knapmeyer 等[34]構(gòu)建的包含 13681 個火星震的目錄, 火星震事件的震源深度最大可達 100km, 大多數(shù)事件的震源深度小于 60km, 因此本文模擬的震源深度為 60km。震源時間函數(shù)的寬度為 10s, 計算中使用的時間間隔?由穩(wěn)定性條件確定, 即受模型中最小網(wǎng)格間距?min與最大波速max之比的約束。時間步長?滿足

    其中,為常系數(shù)0.26。中心附近的?min為0.02km,max為10.5km/s, 因此得出的?為0.0004s。該值在實際計算中過小, 因此需要設(shè)法加大?。本文參考Wang等[30]提出的波數(shù)域濾波算法, 解決鄰近球心=0處橫向網(wǎng)格間距過小的問題。在計算橫向?qū)?shù)時應(yīng)用一個平滑因子, 即應(yīng)用低通濾波器濾除高頻波數(shù)成分, 由此得到的?0.049, 比濾波前的時間步長提高約100倍。本研究取?值為0.04s, 計算的總步數(shù)為75000, 能得到3000s的理論地震圖。

    在研究切向分量波傳播時, 由于該分量上只有S 波, 且 S 波在火星液核內(nèi)不傳播, 因此對模擬參數(shù)做了調(diào)整, 將火星核內(nèi)的網(wǎng)格點剔除。方向離散為 2048 個網(wǎng)格點, 范圍為 0~2π;方向離散為 550個網(wǎng)格點, 范圍為 1190~3390km, 因此方向的網(wǎng)格間距為 5.0km, 同時方向的網(wǎng)格能完整覆蓋 SH波的傳播路徑。

    2.1 模型 A 的理論地震圖與波場快照

    圖 3 展示模型 A 在火星表面的合成理論地震圖以及主要震相的理論射線到達時間, 可以看到, 由于火星核的存在, P 波傳播到核幔邊界時速度會明顯減小, 并改變路徑方向, 因此存在著直達 P 波的影區(qū), 影區(qū)范圍的大小與火星核半徑有關(guān)。3 個分量上均能觀測到互相對應(yīng)的直達 S 波, 但只有徑向和垂向分量能觀測到直達 P 波, 同時能觀測到 P 波的多重反射波 PP 和 PPP 等。徑向和垂向分量上均能看到傳播經(jīng)過火星核的 PKP 震相, 只能看到很弱的入射到核幔邊界發(fā)生反射的震相, 如 PcP, PcS 和ScP 等。切向分量上波形單一, 只能觀測到 S 波, 多重反射波 SS 和 SSS 最為清晰明顯, 同時能觀測到 ScS, sScS, ScS2 和 ScS3 等核幔邊界反射震相。由于模型 A 中火星殼厚度較小, 因此 ScS 和 sScS 震相的到時相差不大。在直達 S 波和 SS 波之間, 可以看到很強的低速火星殼內(nèi)部產(chǎn)生的多重反射震相。同時可以看到, 3 個分量的理論地震圖上面波均比較清晰明顯, 徑向和垂向分量上的振幅較大, 切向分量上的勒夫面波沒有另外兩個分量上的瑞利面波強。隨著面波的傳播, 波列越來越長, 頻散現(xiàn)象越來越明顯。

    圖 4 展示震源深度為 60km 的火星震激發(fā)的 P波和 SV 波在全火星模型 A 中傳播的波場快照, 可以看到, 120s 時直達 P 波和 SV 波在火星幔上層向下傳播, 具有清晰的波前, 同時能看到來自地表的直達 P 波和 SV 波的反射波和轉(zhuǎn)換波, pP 和 sP 震相的波前依次在 P 波之后向下傳播。在火星殼中, 速度比較慢的面波開始產(chǎn)生。在低速的火星殼中, 由于P 波和 SV 波的多重反射和轉(zhuǎn)換作用, 形成直達波之后的較長波列。360s 時, 直達 SV 波到達核幔邊界, 可以看到比較清晰的 P 和 pP 震相的核幔邊界反射波, 面波和直達波之后多重反射和轉(zhuǎn)換波的波列變得更強。420s 時, 來自核幔邊界的直達 SV 波的反射波和轉(zhuǎn)換 P 波向上傳播, 由于 P 波影區(qū)的存在, 核幔邊界附近可以看到彎曲的繞射 P 波。720s 時, 可以看到整個波場存在各種“Y”字形 P 波和 SV 波的多重反射波, P 波已經(jīng)傳播至火星核的另一面, 產(chǎn)生轉(zhuǎn)換的 SV 波。

    圖 5 展示 SH 波在全火星模型 A 中傳播的位移波場快照, 能夠很清晰地看到直達 SH 波、多重反射 SS 和 SSS 震相以及傳播到核幔邊界反射的 ScS震相。180s 時, 直達 SH 波和經(jīng)過表面反射的 SH波向下傳播, 波場非常清晰。420s 時, 直達 SH 波到達核幔邊界, 并產(chǎn)生反射SH波向上傳播, ScS 和sScS 震相非常清晰。630s 時, 由于影區(qū)的存在, 核幔邊界附近產(chǎn)生明顯的繞射波。750s 時, 反射 SH波已傳回地表, 能看到火星殼內(nèi)傳播的勒夫面波和低速火星殼內(nèi)由于多重反射產(chǎn)生的很長的波列。

    2.2 模型 B 的理論地震圖與波形快照

    圖 6 展示模型 B 在火星表面的合成理論地震圖以及各個主要震相的理論射線到達時間。與模型 A相比, 模型 B 的火星核半徑更大, 因此直達 P 波的影區(qū)范圍也更大。切向上, 由于火星核半徑更大, 直達 SH 波的傳播距離更短, ScS 震相的到時也比模型 A 更早。同時, 由于模型 B 的火星殼更厚, 低速區(qū)內(nèi) SH 波的傳播距離更長, 因此 ScS 與 sScS 震相的到時相差更大。二次反射的 ScS2 震相和三次反射的 ScS3 震相等也是如此。由于切向上只有 SH波, 徑向和垂向上 P 波和 SV 波會發(fā)生震相轉(zhuǎn)換, 因此在切向上更容易估算波的傳播距離和傳播時間。

    圖 7 展示震源深度為 60km 的火星震激發(fā)的 P波和 SV 波在全火星模型 B 中傳播的波場快照。180 s 時, 直達 P 波傳播進入火星核, 比模型 A 更早。330s 時, 直達 SV 波到達核幔邊界。360s 時, 由于模型 B 的火星核半徑更大, 因此更早地產(chǎn)生彎曲的繞射波。900 s時, 波已傳至另一側(cè)火星表面, 可以看到很明顯的多重反射震相。

    圖 8 展示 SH 波在全火星模型 B 中傳播的位移波場快照。270s 時, 直達 SH 波到達火星幔上層與下層的分界面。由于模型 B 的火星殼更厚, 因此向下傳播先經(jīng)殼幔邊界反射, 再經(jīng)火星表面二次反射向下傳播的 SmSS 和 sSmSS 震相與直達 SH 波和 sS波的距離更遠。630s 時, 開始產(chǎn)生核幔邊界繞射波。690s 時, 直達 SH 波經(jīng)過核幔邊界的反射波傳回地表。

    2.3 模型 A 與模型 B 理論地震圖對比

    為了討論不同界面位置對火星震波傳播的影響, 將模型 A 與模型 B 的理論地震圖進行對比分析, 如圖 9 所示??梢钥闯? 由于模型 A 的火星殼厚度更小, 低速區(qū)更薄, 因此直達 P 波的到時相對更早。震中距離較近時, 波形差異很小; 隨著距離增加, 直達波的傳播深度增加, 受火星殼幔速度結(jié)構(gòu)的影響, 其波形與到時差稍有增加。多重反射震相受到的影響較大, 這是因為多重反射震相 PP, PPP, SS 和 SSS 等在傳播過程中多次穿過低速的火星殼, 受火星殼厚度的影響很大。隨著反射次數(shù)增加, 其到時差別越來越大。PcP 和 PKP 震相的到時均與火星核半徑有關(guān), P 波傳播到火星核內(nèi)之后速度變慢, 因此對于 PcP 震相, 火星核半徑更大的模型 B 到時更早; 對于 PKP 震相, 火星核內(nèi)傳播路徑更長的模型 B 到時更晚。同時, 模型 A 的面波頻散要強一些, 面波的幅度更大, 持續(xù)時間更長, 面波更發(fā)育。切向上, 火星核半徑更大, 直達 SH 波的傳播距離更短, 模型 B 中 ScS 震相的到時也相對更早, 在圖 9中能看到明顯的差異。同時, 由于模型 B 的火星殼更厚, 因此 ScS 與 sScS 震相的到時相差也更大。與瑞利面波相比, 切向上勒夫面波的差異較小。

    3 結(jié)論和討論

    本文采用基于交錯網(wǎng)格的有限差分與傅里葉偽譜法的混合方法, 對目前廣泛采用的兩個全火星模型開展震波傳播特征的數(shù)值模擬研究, 得到以下結(jié)論。

    1)火星殼厚度對火星震波的傳播起著重要的作用, 在低速火星殼內(nèi)部, 多重反射波與轉(zhuǎn)換波相干疊加會形成很強的波列?;鹦潜砻娴亩嘀胤瓷洳〞啻谓?jīng)過低速的火星殼, 受火星殼厚度的影響很大。

    2)核幔邊界深度影響核幔邊界反射波的到達時間和振幅。由于在切向分量上能夠更清晰地觀測到核幔邊界反射震相, 因此利用切向分量記錄到的核幔邊界反射波可以更方便地進行核幔邊界特征的研究。

    3)通過對比模型 A 與模型 B 的波形圖和波場快照, 可知模型 A 具有持續(xù)時間更長、頻散更強的面波, 這與 Toyokuni 等[33]的結(jié)論相符, 即模型差異嚴重地影響面波的波列, 較薄的火星殼模型中面波的波列更長。

    本文數(shù)值模擬的震源時間函數(shù)的主周期為 10 s, 在以后的研究中, 可以考慮不同的主周期, 進一步研究火星殼對不同成分面波傳播的影響。模擬過程中采用的震源為二維線源, 與實際的三維點源相比, 由于存在波形和幾何擴散的差異, 二者的模擬結(jié)果有差異, 可以利用二維線源到三維點源結(jié)果的轉(zhuǎn)換進行近似的校正[30]。另外, 后續(xù)研究中可以嘗試對火星震波在三維空間中的傳播進行數(shù)值模擬, 并且討論橫向非均勻火星殼和地形變化對火星震波傳播的影響。

    [1] 歐陽自遠, 肖福根. 火星探測的主要科學問題. 航天器環(huán)境工程, 2011, 28(3): 205?217

    [2] Lognonné P, Johnson C. Planetary seismology. Trea-tise on Geophysics, 2015, 10: 65–120

    [3] Anderson D L, Miller W F, Latham G V, et al. Seis-mology on Mars. Journal of Geophysical Research, 1977, 82(28): 4524?4546

    [4] Goins N R, Lazarewicz A R. Martian seismicity. Geo-physical Research Letters, 1979, 6(4): 368?370

    [5] Nakamura Y, Anderson D L. Martian wind activity detected by a seismometer at Viking Lander 2 site. Geophysical Research Letters, 1979, 6(6): 499?502

    [6] Solomon S C, Anderson D L, Banerdt W B. Scientific rationale and requirements for a global seismic net-work on Mars [R]. LPI Tech. Rpt. 91-02. Houston: Lunar and Planetary Institute, 1991

    [7] Nakamura Y. New identification of deep moonquakes in the Apollo lunar seismic data. Physics of the Earth and Planetary Interiors, 2003, 139: 197?205

    [8] Yoder C F, Konopliv A S, Yuan D N, et al. Fluid core size of Mars from detection of the Solar tide. Science, 2003, 300: 299?230

    [9] Konopliv A S, Asmar S W, Folkner W M, et al. Mars high resolution gravity fields from MRO, Mars sea-sonal gravity, and other dynamical parameters. Icarus, 2011, 211: 401?428

    [10] Konopliv A S, Park R S, Folkner W M. An improved JPL Mars gravity field and orientation from Mars orbiter and lander tracking data. Icarus, 2016, 274: 253?260

    [11] Sohl F, Spohn T. The interior structure of Mars: imp-lications from SNC meteorites. Journal of Geophysi-cal Research, 1997, 102: 1613?1635

    [12] Lognonné P, Banerdt W B, Giardini, et al. SEIS: insight’s seismic experiment for internal structure of Mars. Space Science Reviews, 2019, 215: 1?170

    [13] Brown D, Johnson A, Good A. NASA’s InSight Lander captures audio of first likely ‘Quake’ on Mars [EB/OL]. (2019?04?24) [2019?05?01]. https://www. nasa.gov/press-release/nasa-s-insight-lander-captures- audio-of-first-likely-quake-on-mars

    [14] Clinton J F, Giardini D, B?se M, et al. The mars-quake service — building a Martian seismicity catalo-gue for InSight. Space Sci Rev, 2018, 214: 1?33

    [15] Knapmeyer-Endrun B, Ceylan S, van Driel M, Crustal S-wave velocity from apparent incidence angles: a case study in preparation for InSight. Space Sci Rev, 2018, 214: 1?40

    [16] Zheng Y, Nimmo F, Lay T. Seismological implications of a lithospheric low seismic velocity zone in Mars. Phys Earth Planet Inter, 2015, 240: 132–141

    [17] Khan A, van Driel M, B?se M, et al. Single-station and single-event marsquake location and inversion for structure using synthetic Martian waveforms. Phys Earth Planet Inter, 2016, 258: 28–42

    [18] Lognonné P, Karakostas F, Rolland L, et al. Modeling of atmospheric-coupled Rayleigh waves on planets with atmosphere: from Earth observation to Mars and Venus perspectives. J Acoust Soc Am, 2016, 140(2): 1447–1468

    [19] Bissig F, Khan A, van Driel M, et al. On the detectability and use of normal modes for determining interior structure of Mars. Space Sci Rev, 2018, 214: 1?28

    [20] Bozdag E, Ruan Y, Metthez N, et al. Simulations of seismic wave propagation on Mars. Space Sci Rev, 2017, 211(1/2/3/4): 571–594

    [21] Furumura T, Koketsu K, Wen K L. Parallel PSM/FDM hybrid simulation of ground motions from the 1999 Chi-Chi, Taiwan, earthquake. Pure and Applied Geo-physics, 2002, 159(9): 2133?2146

    [22] 馬德堂, 朱光明. 有限元法與偽譜法混合求解彈性波動方程. 地球科學與環(huán)境學報, 2004, 26(1): 61?64

    [23] 魏星, 王彥賓, 陳曉非. 模擬地震波場的偽譜和高階有限差分混合方法. 地震學報, 2010, 32(4): 392?400

    [24] 秦艷芳, 王彥賓. 地震波傳播的三維偽譜和高階有限差分混合方法并行模擬. 地震學報, 2012, 34(2): 147?156

    [25] Wang Y, Takenaka H, Jiang X, et al. Modelling two-dimensional global seismic wave propagation in a laterally heterogeneous whole-Moon model. Geophy-sical Journal International, 2013, 192(3): 1271?1287

    [26] Jiang X, Wang Y, Qin Y, et al. Global SH-wave pro-pagation in a 2D whole Moon model using the parallel hybrid PSM/FDM method. Earthquake Science, 2015, 28(3): 163?174

    [27] Witte D C, Richards P G. The pseudospectral method for simulating wave propagation. Computational Acoustics, 1990, 3: 1?18

    [28] Levander A R. Fourth-order finite-difference P-SV seismograms. Geophysics, 1988, 53(11): 1425?1436

    [29] Cerjan C, Kosloff D, Kosloff R, et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equations. Geophysics, 1985, 50(4): 705?708

    [30] Wang Y, Takenaka H, Furumura T. Modelling seismic wave propagation in a two-dimensional cylindrical whole-earth model using the pseudospectral method. Geophysical Journal International, 2001, 145(3): 689?708

    [31] Wang Y, Takenaka H. SH-wavefield simulation for a laterally heterogeneous whole-Earth model using the pseudospectral method. Science China Earth Scien-ces, 2011, 54(12): 1940?1947

    [32] Herrmann R B. SH-wave generation by dislocation source —a numerical study. Bulletin of the Seis-mological Society of America, 1979, 69: 1?15

    [33] Toyokuni G, Ishihara Y, Takenaka H. Preliminary modeling of global seismic wave propagation in the whole Mars // 42nd Lunar and Planetary Science Conference. Woodlands, 2011: Abstract no. 1631

    [34] Knapmeyer M, Oberst J, Hauber E, et al. Working models for spatial distribution and level of Mars’ seis-micity. Journal of Geophysical Research, 2006, 111: 11006

    Numerical Modeling of Global Seismic Wave Propagation in the Whole Mars Models

    DENG Di, XIAO Wanbo, WANG Yanbin?

    Department of Geophysics, School of Earth and Space Sciences, Peking University, Beijing 100871; ? Corresponding author, E-mail: ybwang@pku.edu.cn

    The pseudospectral and finite difference hybrid method on staggered grid is applied to solve seismic wave equations for two Martian structure models derived from geochemical analysis. The numerical modeling is used to calculate P-SV and SH wave propagation inside 2-D whole Mars models. The generation and propagation of various seismic phases in the whole Mars models are shown by synthetic seismograms and wavefield snapshots. Effect of Martian crustal thickness and the depth of Martian core-mantle boundary on seismic wave propagation is analyzed with synthetic seismograms. Multiple reflections and conversions of seismic waves and their constructive interference inside the low-velocity Martian crust form reverberating wave trains, which are strongly affected by the thickness of Martian crust. Seismic reflections from core-mantle boundary can be clearly identified from the calculated transverse component seismogram.

    whole Mars model; seismic wave propagation; numerical modeling; pseudospectral and finite difference hybrid method

    10.13209/j.0479-8023.2020.029

    國家自然科學基金(41930103)資助

    2019?05?16;

    2019?10?13

    猜你喜歡
    面波波場震源
    gPhone重力儀的面波頻段響應(yīng)實測研究
    地震研究(2021年1期)2021-04-13 01:04:56
    自適應(yīng)相減和Curvelet變換組合壓制面波
    彈性波波場分離方法對比及其在逆時偏移成像中的應(yīng)用
    震源的高返利起步
    交錯網(wǎng)格與旋轉(zhuǎn)交錯網(wǎng)格對VTI介質(zhì)波場分離的影響分析
    地震學報(2016年1期)2016-11-28 05:38:36
    基于Hilbert變換的全波場分離逆時偏移成像
    可控震源地震在張掖盆地南緣逆沖斷裂構(gòu)造勘探中的應(yīng)用
    同步可控震源地震采集技術(shù)新進展
    旋轉(zhuǎn)交錯網(wǎng)格VTI介質(zhì)波場模擬與波場分解
    淺析工程勘探的面波勘探方法
    河南科技(2014年8期)2014-02-27 14:07:40
    最新美女视频免费是黄的| 国产精品香港三级国产av潘金莲| 亚洲电影在线观看av| 国产伦在线观看视频一区| 日韩大尺度精品在线看网址| 成熟少妇高潮喷水视频| 黄色片一级片一级黄色片| 亚洲国产精品成人综合色| 精品免费久久久久久久清纯| 黄色视频,在线免费观看| 欧美黑人欧美精品刺激| 美女免费视频网站| 一级a爱片免费观看的视频| 亚洲欧美日韩无卡精品| 亚洲 国产 在线| 亚洲 欧美 日韩 在线 免费| 老司机在亚洲福利影院| 亚洲国产日韩欧美精品在线观看 | 久9热在线精品视频| 国产精华一区二区三区| 亚洲人成77777在线视频| 男女床上黄色一级片免费看| 国产精品香港三级国产av潘金莲| 欧美在线黄色| 久久国产乱子伦精品免费另类| 神马国产精品三级电影在线观看 | 亚洲国产欧美日韩在线播放| 亚洲精品国产区一区二| 黄色 视频免费看| 亚洲午夜精品一区,二区,三区| 久久久久国产精品人妻aⅴ院| 一级毛片高清免费大全| 成人一区二区视频在线观看| 一级毛片女人18水好多| 黄色a级毛片大全视频| 欧美一级a爱片免费观看看 | 国产在线精品亚洲第一网站| 亚洲国产精品久久男人天堂| 国产成人啪精品午夜网站| 成人18禁高潮啪啪吃奶动态图| 在线观看日韩欧美| 精品久久久久久久末码| cao死你这个sao货| 在线观看午夜福利视频| 十八禁人妻一区二区| 十八禁人妻一区二区| 最近在线观看免费完整版| 在线观看免费午夜福利视频| 欧美日韩福利视频一区二区| 人人妻人人澡人人看| 久久久国产成人免费| 精品熟女少妇八av免费久了| bbb黄色大片| 久久久国产成人精品二区| 午夜福利成人在线免费观看| 99热这里只有精品一区 | 无限看片的www在线观看| 国产亚洲精品av在线| 91在线观看av| 中文字幕精品免费在线观看视频| 欧美绝顶高潮抽搐喷水| 亚洲一区中文字幕在线| 少妇 在线观看| 国产一区二区三区视频了| 好男人电影高清在线观看| 18禁观看日本| 99国产极品粉嫩在线观看| 亚洲第一av免费看| 欧美激情 高清一区二区三区| 亚洲av成人av| 女生性感内裤真人,穿戴方法视频| 久久久国产精品麻豆| 欧美日韩亚洲综合一区二区三区_| 精品国产乱码久久久久久男人| 啪啪无遮挡十八禁网站| 久久亚洲真实| 午夜福利视频1000在线观看| 国产熟女午夜一区二区三区| 国内精品久久久久久久电影| 成人国语在线视频| 欧美乱妇无乱码| 亚洲精品在线观看二区| 亚洲精品一区av在线观看| 日韩三级视频一区二区三区| 亚洲人成电影免费在线| 国产精品二区激情视频| 亚洲精品美女久久久久99蜜臀| АⅤ资源中文在线天堂| 亚洲狠狠婷婷综合久久图片| 18禁美女被吸乳视频| 1024视频免费在线观看| 中文字幕人妻熟女乱码| 国产97色在线日韩免费| 久久伊人香网站| 久久青草综合色| 久久久久免费精品人妻一区二区 | 97人妻精品一区二区三区麻豆 | 夜夜爽天天搞| 久久性视频一级片| 黄色片一级片一级黄色片| 欧美黄色淫秽网站| 国产v大片淫在线免费观看| 日韩欧美在线二视频| 人人澡人人妻人| 国产亚洲精品综合一区在线观看 | 国产精品亚洲美女久久久| 亚洲一区中文字幕在线| 国产午夜福利久久久久久| ponron亚洲| 午夜福利欧美成人| 国产av不卡久久| 少妇裸体淫交视频免费看高清 | 色综合站精品国产| 久久久久久久午夜电影| 免费女性裸体啪啪无遮挡网站| av超薄肉色丝袜交足视频| 女警被强在线播放| 国产视频内射| 国产激情久久老熟女| 精品一区二区三区av网在线观看| 国产精品电影一区二区三区| 一本大道久久a久久精品| 国产精品98久久久久久宅男小说| 亚洲国产高清在线一区二区三 | 国产av不卡久久| 叶爱在线成人免费视频播放| 又大又爽又粗| 一进一出抽搐动态| 91麻豆av在线| 制服人妻中文乱码| 日日摸夜夜添夜夜添小说| 天天一区二区日本电影三级| 婷婷精品国产亚洲av| 少妇粗大呻吟视频| 美女午夜性视频免费| 2021天堂中文幕一二区在线观 | 欧美日韩福利视频一区二区| 99久久精品国产亚洲精品| www.熟女人妻精品国产| 人人澡人人妻人| 在线观看舔阴道视频| 午夜激情av网站| 国产亚洲精品久久久久久毛片| 18禁观看日本| 97人妻精品一区二区三区麻豆 | 欧美乱妇无乱码| 久久久久久亚洲精品国产蜜桃av| 九色国产91popny在线| 亚洲自偷自拍图片 自拍| 欧美日韩精品网址| 成人av一区二区三区在线看| 国产精品久久久av美女十八| 18禁国产床啪视频网站| 久久久精品欧美日韩精品| 亚洲av第一区精品v没综合| 91国产中文字幕| 女人被狂操c到高潮| 怎么达到女性高潮| 亚洲最大成人中文| 久久婷婷成人综合色麻豆| 女人高潮潮喷娇喘18禁视频| avwww免费| 国产欧美日韩精品亚洲av| 日韩视频一区二区在线观看| 91国产中文字幕| 亚洲人成伊人成综合网2020| 日本三级黄在线观看| 91老司机精品| 成人永久免费在线观看视频| 老熟妇乱子伦视频在线观看| 91字幕亚洲| 岛国视频午夜一区免费看| 少妇裸体淫交视频免费看高清 | 88av欧美| 搡老妇女老女人老熟妇| 欧美日韩一级在线毛片| 亚洲国产精品久久男人天堂| 亚洲国产欧美网| 亚洲成人免费电影在线观看| 91成人精品电影| 变态另类丝袜制服| 亚洲国产欧美一区二区综合| 18禁黄网站禁片午夜丰满| 国产高清videossex| 亚洲中文日韩欧美视频| 变态另类成人亚洲欧美熟女| 欧美av亚洲av综合av国产av| 久久欧美精品欧美久久欧美| 18禁国产床啪视频网站| 欧美zozozo另类| 久久热在线av| 99re在线观看精品视频| 精品无人区乱码1区二区| 手机成人av网站| 波多野结衣高清作品| www国产在线视频色| 国产亚洲精品久久久久久毛片| 露出奶头的视频| 美女扒开内裤让男人捅视频| 国产精品综合久久久久久久免费| 久久精品91无色码中文字幕| 真人做人爱边吃奶动态| 1024香蕉在线观看| 色av中文字幕| av天堂在线播放| 在线观看午夜福利视频| 少妇熟女aⅴ在线视频| 亚洲男人天堂网一区| 国产区一区二久久| 999久久久国产精品视频| 久久亚洲真实| 熟妇人妻久久中文字幕3abv| 免费高清视频大片| 狂野欧美激情性xxxx| 亚洲精品在线美女| 日本免费a在线| 制服人妻中文乱码| 午夜福利18| 高清毛片免费观看视频网站| 一本大道久久a久久精品| 中文字幕人妻丝袜一区二区| 欧美中文日本在线观看视频| 精品电影一区二区在线| 国产成人精品久久二区二区91| 高潮久久久久久久久久久不卡| 久久国产乱子伦精品免费另类| 成人手机av| 久久精品国产综合久久久| 久久亚洲精品不卡| 最新在线观看一区二区三区| 黄网站色视频无遮挡免费观看| 国内毛片毛片毛片毛片毛片| 免费在线观看日本一区| 国产精品免费一区二区三区在线| 欧美一级毛片孕妇| av中文乱码字幕在线| 十八禁人妻一区二区| 成人永久免费在线观看视频| 国产精品一区二区三区四区久久 | 999久久久国产精品视频| av超薄肉色丝袜交足视频| 色播在线永久视频| 亚洲精品粉嫩美女一区| 91老司机精品| 国产高清激情床上av| 人人妻人人澡人人看| 一级黄色大片毛片| 亚洲七黄色美女视频| 国产激情偷乱视频一区二区| 国产男靠女视频免费网站| 欧美乱妇无乱码| 亚洲午夜精品一区,二区,三区| 色综合站精品国产| 免费在线观看日本一区| 国产一卡二卡三卡精品| 欧美国产日韩亚洲一区| 美女午夜性视频免费| 无人区码免费观看不卡| 999精品在线视频| 老鸭窝网址在线观看| 色播在线永久视频| 国产成年人精品一区二区| 欧美日韩亚洲国产一区二区在线观看| 人人妻人人澡欧美一区二区| 精品熟女少妇八av免费久了| 不卡一级毛片| 嫁个100分男人电影在线观看| 欧美性长视频在线观看| 久久中文字幕人妻熟女| 热re99久久国产66热| 日本精品一区二区三区蜜桃| 国产麻豆成人av免费视频| 久久国产精品影院| 日本 欧美在线| 国产精品电影一区二区三区| 亚洲av美国av| 日韩欧美国产在线观看| 久久狼人影院| 婷婷丁香在线五月| 国产成人系列免费观看| 国产视频内射| 日本五十路高清| 好男人在线观看高清免费视频 | 国产极品粉嫩免费观看在线| 身体一侧抽搐| 久久九九热精品免费| 亚洲国产精品合色在线| or卡值多少钱| 日日摸夜夜添夜夜添小说| 欧美激情 高清一区二区三区| 一卡2卡三卡四卡精品乱码亚洲| 男人舔奶头视频| 91国产中文字幕| 欧美精品亚洲一区二区| 亚洲国产精品成人综合色| 好看av亚洲va欧美ⅴa在| 国产一区二区在线av高清观看| 观看免费一级毛片| 久久精品aⅴ一区二区三区四区| 久久久国产成人精品二区| 亚洲 国产 在线| 制服诱惑二区| 久久婷婷成人综合色麻豆| 久久久久国内视频| 一边摸一边做爽爽视频免费| 久久精品91无色码中文字幕| 国产精品98久久久久久宅男小说| 欧美黑人巨大hd| 国产精品乱码一区二三区的特点| 欧美激情高清一区二区三区| 成人永久免费在线观看视频| 日韩欧美免费精品| 亚洲自偷自拍图片 自拍| 国产一区二区三区在线臀色熟女| 看黄色毛片网站| 在线天堂中文资源库| 黄片播放在线免费| 波多野结衣高清作品| 亚洲天堂国产精品一区在线| 日本 av在线| 国产v大片淫在线免费观看| 女人高潮潮喷娇喘18禁视频| 桃色一区二区三区在线观看| 人人妻,人人澡人人爽秒播| 国产精品美女特级片免费视频播放器 | 国产精品一区二区三区四区久久 | 一a级毛片在线观看| 欧美黄色淫秽网站| 在线视频色国产色| 91麻豆av在线| 精华霜和精华液先用哪个| 国产三级黄色录像| 成人18禁高潮啪啪吃奶动态图| 黄色女人牲交| 国产人伦9x9x在线观看| 国产欧美日韩一区二区三| 亚洲国产精品久久男人天堂| 男女午夜视频在线观看| 99久久精品国产亚洲精品| 自线自在国产av| 国产精品美女特级片免费视频播放器 | 欧美大码av| 国产精品久久久久久精品电影 | 国产黄a三级三级三级人| 国产成人av激情在线播放| 久久久久久久久久黄片| 黄片播放在线免费| 天堂√8在线中文| 淫妇啪啪啪对白视频| 亚洲,欧美精品.| 制服人妻中文乱码| 成人亚洲精品av一区二区| 色在线成人网| 成人三级做爰电影| 一a级毛片在线观看| 老司机福利观看| 巨乳人妻的诱惑在线观看| 最新在线观看一区二区三区| 国产在线观看jvid| 亚洲av电影不卡..在线观看| 久久99热这里只有精品18| 99热这里只有精品一区 | 国语自产精品视频在线第100页| 欧美成人性av电影在线观看| 欧美不卡视频在线免费观看 | 亚洲男人的天堂狠狠| 美女大奶头视频| 少妇的丰满在线观看| 成人三级黄色视频| 99久久综合精品五月天人人| 日韩 欧美 亚洲 中文字幕| 精品日产1卡2卡| 一级a爱片免费观看的视频| 亚洲第一欧美日韩一区二区三区| av欧美777| 欧洲精品卡2卡3卡4卡5卡区| 国产一卡二卡三卡精品| 高潮久久久久久久久久久不卡| 十八禁网站免费在线| 一本精品99久久精品77| 免费高清视频大片| 男男h啪啪无遮挡| 色综合站精品国产| 亚洲一区二区三区不卡视频| 狠狠狠狠99中文字幕| 少妇熟女aⅴ在线视频| 国产成人精品无人区| 欧美日韩瑟瑟在线播放| 一个人观看的视频www高清免费观看 | 后天国语完整版免费观看| 欧美亚洲日本最大视频资源| 99久久国产精品久久久| 色在线成人网| 国产成人啪精品午夜网站| 国产高清激情床上av| www.熟女人妻精品国产| 国产精品二区激情视频| 一区二区三区国产精品乱码| 制服丝袜大香蕉在线| 色婷婷久久久亚洲欧美| 久久精品成人免费网站| 精品国产超薄肉色丝袜足j| 91国产中文字幕| 丝袜在线中文字幕| 人人妻人人澡人人看| 国产av在哪里看| 亚洲成国产人片在线观看| www日本黄色视频网| 亚洲av电影不卡..在线观看| 国产私拍福利视频在线观看| 99精品在免费线老司机午夜| 欧美性猛交╳xxx乱大交人| 熟妇人妻久久中文字幕3abv| 日本 欧美在线| 国产亚洲精品久久久久久毛片| 在线观看免费日韩欧美大片| 日韩免费av在线播放| 免费在线观看亚洲国产| 美女午夜性视频免费| 久久婷婷人人爽人人干人人爱| 校园春色视频在线观看| 亚洲,欧美精品.| 亚洲 欧美一区二区三区| 午夜福利免费观看在线| 欧美绝顶高潮抽搐喷水| 久久精品91无色码中文字幕| 丝袜在线中文字幕| xxxwww97欧美| 日本 av在线| 最近最新中文字幕大全电影3 | √禁漫天堂资源中文www| 国产私拍福利视频在线观看| 亚洲精品av麻豆狂野| 亚洲性夜色夜夜综合| 亚洲va日本ⅴa欧美va伊人久久| 夜夜爽天天搞| 色哟哟哟哟哟哟| 久久久精品国产亚洲av高清涩受| 在线视频色国产色| 正在播放国产对白刺激| 啦啦啦 在线观看视频| netflix在线观看网站| 中文字幕另类日韩欧美亚洲嫩草| 欧美黑人巨大hd| 国产精品99久久99久久久不卡| 黄片大片在线免费观看| 午夜福利高清视频| 最新在线观看一区二区三区| 久久精品国产亚洲av香蕉五月| 久久午夜综合久久蜜桃| 色综合站精品国产| 女性被躁到高潮视频| 丝袜人妻中文字幕| 大型黄色视频在线免费观看| 性欧美人与动物交配| 久久婷婷人人爽人人干人人爱| 淫妇啪啪啪对白视频| 淫秽高清视频在线观看| 色婷婷久久久亚洲欧美| 91国产中文字幕| 自线自在国产av| 中出人妻视频一区二区| 免费在线观看完整版高清| 一级毛片高清免费大全| 亚洲av熟女| 一区二区三区激情视频| 日韩一卡2卡3卡4卡2021年| 丁香欧美五月| 久久人妻av系列| 欧美性猛交黑人性爽| 在线十欧美十亚洲十日本专区| 久久久久国内视频| 大型黄色视频在线免费观看| 色av中文字幕| 老汉色av国产亚洲站长工具| 一本大道久久a久久精品| 中亚洲国语对白在线视频| 久久久久久免费高清国产稀缺| x7x7x7水蜜桃| av欧美777| 亚洲欧美日韩高清在线视频| 天天添夜夜摸| 一级黄色大片毛片| 国产成人av教育| 精品国产一区二区三区四区第35| 夜夜躁狠狠躁天天躁| 国产精品九九99| 一本大道久久a久久精品| 一边摸一边抽搐一进一小说| www.自偷自拍.com| 国产成人啪精品午夜网站| 长腿黑丝高跟| 好看av亚洲va欧美ⅴa在| 久久天躁狠狠躁夜夜2o2o| 国产高清视频在线播放一区| 国产精品亚洲av一区麻豆| 亚洲激情在线av| 热re99久久国产66热| 亚洲国产精品成人综合色| 色哟哟哟哟哟哟| 亚洲一区高清亚洲精品| 黄片小视频在线播放| 亚洲国产欧美网| 日本 欧美在线| 国内少妇人妻偷人精品xxx网站 | 很黄的视频免费| 日本 av在线| 国产激情欧美一区二区| 日本a在线网址| 久久久久久免费高清国产稀缺| 成年免费大片在线观看| 精品乱码久久久久久99久播| 美国免费a级毛片| a在线观看视频网站| 99国产极品粉嫩在线观看| 亚洲第一电影网av| 欧美日韩中文字幕国产精品一区二区三区| 白带黄色成豆腐渣| 一级作爱视频免费观看| 亚洲国产看品久久| 日日爽夜夜爽网站| 亚洲国产精品成人综合色| 一区二区三区高清视频在线| 欧美日韩亚洲国产一区二区在线观看| 日韩三级视频一区二区三区| 成人国产综合亚洲| 日本一区二区免费在线视频| 久久精品aⅴ一区二区三区四区| 丁香六月欧美| 精品免费久久久久久久清纯| 欧美亚洲日本最大视频资源| 亚洲欧美精品综合久久99| 黄片小视频在线播放| 美国免费a级毛片| 嫁个100分男人电影在线观看| 亚洲精品一区av在线观看| 中文字幕最新亚洲高清| 国产精品av久久久久免费| 一级黄色大片毛片| 国产视频内射| 99热只有精品国产| 亚洲va日本ⅴa欧美va伊人久久| 久久午夜亚洲精品久久| 国产在线观看jvid| 成人一区二区视频在线观看| 欧美黄色片欧美黄色片| 亚洲一区二区三区不卡视频| 久久九九热精品免费| 成人三级黄色视频| 免费观看精品视频网站| 人成视频在线观看免费观看| 成人av一区二区三区在线看| 国产激情久久老熟女| 亚洲无线在线观看| 麻豆成人av在线观看| 国产成人一区二区三区免费视频网站| 国产欧美日韩一区二区三| 好看av亚洲va欧美ⅴa在| 青草久久国产| 免费看a级黄色片| 亚洲人成网站高清观看| 国内毛片毛片毛片毛片毛片| 亚洲中文字幕日韩| 久久久久国产一级毛片高清牌| 国产av不卡久久| 国产乱人伦免费视频| 国产国语露脸激情在线看| 免费av毛片视频| 亚洲午夜理论影院| 青草久久国产| 亚洲国产看品久久| 人人妻人人看人人澡| 高清毛片免费观看视频网站| 国产精品久久久av美女十八| 可以在线观看毛片的网站| 最近最新免费中文字幕在线| 男女午夜视频在线观看| 怎么达到女性高潮| 黄频高清免费视频| 久久久久久九九精品二区国产 | 国产亚洲精品久久久久久毛片| 久久久国产成人精品二区| 国产欧美日韩精品亚洲av| 免费无遮挡裸体视频| 三级毛片av免费| 黄片小视频在线播放| 在线十欧美十亚洲十日本专区| 国产精品美女特级片免费视频播放器 | 亚洲欧美精品综合一区二区三区| 免费在线观看成人毛片| 午夜福利一区二区在线看| 久久精品亚洲精品国产色婷小说| 亚洲一卡2卡3卡4卡5卡精品中文| 久久婷婷人人爽人人干人人爱| 香蕉久久夜色| 99精品久久久久人妻精品| 在线永久观看黄色视频| 亚洲自偷自拍图片 自拍| 在线观看日韩欧美| 国产欧美日韩一区二区三| 欧美绝顶高潮抽搐喷水| 可以在线观看的亚洲视频| 色综合亚洲欧美另类图片| 欧美绝顶高潮抽搐喷水| 国产精品美女特级片免费视频播放器 | 变态另类成人亚洲欧美熟女| 在线av久久热| 国产亚洲欧美在线一区二区| 麻豆成人av在线观看| 久久久久国产一级毛片高清牌| 亚洲av日韩精品久久久久久密| 亚洲电影在线观看av| 国产激情久久老熟女| 免费看美女性在线毛片视频| 99久久久亚洲精品蜜臀av|