999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

川滇塊體東邊界斷裂帶應(yīng)力演化特征及地震危險性數(shù)值

2024-04-02 06:05:20李彥欣董培育
大地測量與地球動力學(xué) 2024年4期
關(guān)鍵詞:模型

李彥欣 董培育 趙 斌

1 中國地震局地震研究所,武漢市洪山側(cè)路40號,430071

2 中國地震局地震大地測量重點實驗室,武漢市洪山側(cè)路40號,430071

川滇塊體東邊界主要由鮮水河斷裂帶、安寧河-則木河斷裂帶、大涼山斷裂帶等多條斷裂帶組成,是我國大陸內(nèi)部構(gòu)造活動性最強的大型左旋走滑斷裂帶系。1725年以來,該地區(qū)先后發(fā)生6級以上強震17次,其中7級以上破壞性強震8次[1]。與鮮水河斷裂帶交叉相鄰的龍門山斷裂帶是一條巨型逆沖推覆構(gòu)造帶,在2008年和2013年相繼發(fā)生MS8.0汶川地震和MS7.0蘆山地震。這2條大型斷裂帶系構(gòu)成一個“Y”字形,將青藏高原東南緣地區(qū)分割為川滇塊體、巴顏喀拉塊體和華南塊體3大活動塊體[2],如圖1所示。

LMSF:龍門山斷裂帶;XSHF:鮮水河斷裂帶;LTF:理塘斷裂帶;ANHF:安寧河斷裂帶;ZMHF:則木河斷裂帶;DLSF:大涼山斷裂帶;LJ-XJHF:麗江-小金河斷裂帶;YLXF:玉龍希斷裂帶;藍(lán)色震源球表示1893年以來發(fā)生在該區(qū)域的MS>6.5強震;紅色震源球表示2022-09-05瀘定MS6.8地震;黑線表示主要斷裂;紅色虛線橢圓區(qū)域表示安寧河-則木河斷裂帶上地震空區(qū)

2022-09-05研究區(qū)再次發(fā)生MS6.8地震,震中位于四川省瀘定縣(29.59°N,102.08°E),發(fā)震斷層為鮮水河斷裂帶南段(磨西段)。此次地震為一次典型的左旋走滑型地震,走向、傾角和滑動角分別約為254°、73°和178°,破裂面總長度達(dá)40 km,最大滑動量約為1.83 m[3-4]。該地震破裂區(qū)內(nèi)的前次地震發(fā)生于1786年,距今已超過200 a[1]。定量計算川滇塊體東邊界主要斷裂帶應(yīng)力積累特征,尤其是瀘定地震破裂區(qū)的應(yīng)力特征,不僅有助于理解該區(qū)域地震孕育和發(fā)生的機理,同時可以為分析區(qū)域未來地震危險性提供定量參考依據(jù)。

應(yīng)力場的演化主要受控于震間長期構(gòu)造活動的加載和地震引起的同震及震后粘彈性松弛效應(yīng)的擾動。對震間長期構(gòu)造應(yīng)力場積累過程的認(rèn)識目前主要依賴于大地測量觀測數(shù)據(jù),通過解算斷層滑動速率及閉鎖程度等信息間接估算應(yīng)力積累速率[5]。也可基于區(qū)域?qū)嶋H地質(zhì)構(gòu)造特征建立數(shù)值計算模型,在合理的邊界條件約束下,直接進(jìn)行定量模擬計算,獲取研究區(qū)域各斷層上應(yīng)力積累速率值[6]。對于地震引起的同震和震后庫侖應(yīng)力變化,從物理上考慮,震源區(qū)斷層破裂會對周邊區(qū)域產(chǎn)生不同程度的應(yīng)力加載或者卸載效應(yīng),從而導(dǎo)致后續(xù)地震提前或者延遲發(fā)生,即庫侖破裂應(yīng)力變化理論[7]。綜合考慮震間長期構(gòu)造應(yīng)力積累和歷史地震應(yīng)力擾動效應(yīng),能夠更全面地定量描述不同區(qū)域應(yīng)力積累速率特征或斷層各區(qū)段應(yīng)力積累特征,從而進(jìn)一步解釋地震的孕育和發(fā)生[8-9]。

本文建立三維粘彈性有限元數(shù)值計算模型,以區(qū)域長期平均 GNSS速度場為約束,采用Pylith開源軟件包進(jìn)行計算。根據(jù)各斷裂帶不同的幾何特征,模擬其震間長期平均構(gòu)造應(yīng)力積累速率,同時考慮周邊歷史強震對其造成的同震和震后庫侖應(yīng)力變化,定量描述構(gòu)造應(yīng)力場的演化過程,探究此次瀘定地震的發(fā)震成因,并綜合分析區(qū)域構(gòu)造環(huán)境特征及未來強震危險性。

1 三維有限元數(shù)值計算模型

1.1 模型構(gòu)建

結(jié)合研究區(qū)域活動斷裂的空間分布和瀘定地震的震中位置,確定研究區(qū)域范圍為98°~106°E、26°~33°N,主要包含川滇塊體、巴顏喀拉塊體和華南塊體3個活動塊體。據(jù)此,建立區(qū)域三維粘彈性有限元模型(圖2),縱向上分為4層,分別為上地殼、中地殼、下地殼和上地幔,其中斷層深度約為20 km,均切穿上地殼。模型剖分為四面體單元,在斷裂帶附近單元尺寸約為3 km,在其他區(qū)域約為20~30 km,共劃分413 401個單元。根據(jù)地震地質(zhì)調(diào)查結(jié)果和前人研究成果,給出各斷裂帶的幾何特征[10],如表1所示。

表1 主要活動斷裂的幾何參數(shù)

圖2 三維有限元模型

1.2 本構(gòu)關(guān)系與介質(zhì)參數(shù)

三維有限元數(shù)值計算模型中,上地殼設(shè)置為彈性材料。采用USTClith2.0模型給出的上地殼S波、P波速度及密度數(shù)據(jù)[11],為充分考慮其橫向不均勻性,將其插值到模型上地殼不同網(wǎng)格中,計算彈性模量和泊松比等參數(shù)。彈性模量的變化范圍為(7.5~9.8)×1010Pa,泊松比為0.20~0.29,其中華南塊體內(nèi)彈性模量最高。模型中下地殼及上地幔設(shè)置為粘彈性介質(zhì),采用線性Maxwell流變學(xué)模型來模擬巖石圈短期內(nèi)的彈性效應(yīng)和長期的粘滯性效應(yīng)。參考前人研究成果,給出各塊體不同圈層的介質(zhì)參數(shù)信息(表2)。

表2 模型中介質(zhì)參數(shù)

1.3 邊界條件

使用Wang等[12]給出的中國大陸約30 a的平均GNSS速度場結(jié)果,將其插值到模型邊界上,假設(shè)速度不隨深度變化,作為模型邊界條件(圖3),設(shè)置模型表面自由、底部切向自由但垂向約束。已有研究證明,川滇地區(qū)軟弱的中下地殼運動速度快于脆性上地殼[13],如果假定上地殼與中下地殼的運動速度無差異,則模擬速度值與觀測值吻合度較低。為解決此問題,在二維平面模型中通常采用等效拖曳力的方法[14];三維模型中可以在下地殼側(cè)邊界上增大速度邊界值[6]。參考前人做法,本文在川滇塊體下地殼側(cè)邊界施加一個速度增量,通過試錯法確定該增量為4 mm/a。

圖3 模型邊界條件設(shè)置

1.4 歷史地震信息

本文主要選擇MS>6.5的歷史地震進(jìn)行計算,但考慮到2014年鮮水河斷裂帶中段康定地區(qū)發(fā)生MS6.3地震,其震源區(qū)緊鄰瀘定地震震源區(qū),可能對瀘定地區(qū)產(chǎn)生一定影響,因此計算中包含該次地震。2008年汶川地震、2013年蘆山地震和2022年瀘定地震都有豐富的觀測數(shù)據(jù)支撐,前人已通過各種反演方法給出這些地震的滑動模型[4,15-17]。而對于早期的地震,由于缺乏數(shù)據(jù)資料,多采用經(jīng)驗公式[18]估算其破裂信息(長度、寬度和滑動量)。具體歷史地震信息如表3所示。

表3 參與計算的地震

2 計算結(jié)果

計算GNSS速度場約束下區(qū)域形變特征、各主要斷層的震間應(yīng)力積累速率和過去130 a研究區(qū)域9次歷史強震同震和震后粘彈性松弛效應(yīng)產(chǎn)生的應(yīng)力變化,以定量描述研究區(qū)域主要斷裂帶的應(yīng)力演化過程。

2.1 震間構(gòu)造應(yīng)力積累速率

將GNSS速度場邊界條件加載下的模型運行至平衡狀態(tài),對比地表GNSS觀測值與模擬值。設(shè)未經(jīng)改正的速度邊界條件的計算模型為M0,增加下地殼速度邊界條件的模型為M1,二者的模擬值與觀測值之間的RMS如表4所示(單位mm/a)。可以看出,M1模型精度優(yōu)于M0模型。圖4(a)顯示,模擬速度值在方向和大小上與觀測值擬合較好,殘差基本小于2 mm/a,個別點位殘差達(dá)到4~5 mm/a(圖4(b)),主要集中在鮮水河斷裂帶和玉龍希斷裂帶圍限區(qū)域,該區(qū)域是上下地殼速度差異最大的集中區(qū),也是前人研究中加載拖曳力最大的區(qū)域[14]。本文通過試錯法在該區(qū)域西側(cè)的下地殼區(qū)域加大邊界條件,結(jié)果得到有效改善,基本反映了現(xiàn)今構(gòu)造變形特征,驗證了模型的合理性。

表4 RMS結(jié)果對比

在M1模型的基礎(chǔ)上計算構(gòu)造加載作用下各斷裂帶的庫侖應(yīng)力積累速率。由于上地殼介質(zhì)和邊界條件的縱向均一性,計算結(jié)果在上地殼范圍內(nèi)是均勻分布的,故僅給出10 km深度處的數(shù)值,如圖5所示。可以看出,鮮水河斷裂帶、理塘斷裂帶和大涼山斷裂帶南段的應(yīng)力積累速率較高,達(dá)到1.0~1.5 kPa/a;安寧河斷裂帶、則木河斷裂帶和大涼山斷裂帶北段次之,約為0.5~1.0 kPa/a;龍門山斷裂帶、玉龍希斷裂帶和麗江-小金河斷裂帶震間應(yīng)力積累速率低于0.5 kPa/a。

圖5 主要斷層庫侖應(yīng)力年積累速率

2.2 瀘定地震同震庫侖應(yīng)力變化

基于2022年瀘定地震同震破裂模型[4],計算其產(chǎn)生的同震應(yīng)力變化張量,并投影到各斷層上(幾何參數(shù)見表1),得到此次地震產(chǎn)生的同震ΔCFS(圖6)。結(jié)果顯示,ΔCFS增加區(qū)域主要集中在震源破裂區(qū)四周,最大值可達(dá) MPa量級。麗江-小金河斷裂帶、安寧河斷裂帶以及大涼山斷裂帶北段都有0.01~0.1 MPa量級的加載效應(yīng),超過應(yīng)力觸發(fā)效應(yīng)的閾值(0.01 MPa),可能會增加這些地區(qū)的地震風(fēng)險。此外,與鮮水河斷裂帶相交的龍門山斷裂帶西南段靠近破裂地區(qū)的ΔCFS有所減小,約為-0.5 MPa。

圖6 瀘定地震引起的周邊各斷裂帶庫侖應(yīng)力變化

2.3 歷史大地震發(fā)展演化過程

2.3.1 歷史地震引起的庫侖應(yīng)力變化

斷層上庫侖應(yīng)力演化的另一個重要因素是歷史地震破裂產(chǎn)生的同震和震后應(yīng)力擾動效應(yīng)。計算過去130 a鮮水河斷裂帶、龍門山斷裂帶發(fā)生的9次大地震導(dǎo)致的各斷層累積ΔCFS,如圖7所示。結(jié)果顯示,鮮水河斷裂帶上1893年乾寧、1904年道孚、1923年仁達(dá)、1973年爐霍以及1981年道孚地震皆導(dǎo)致破裂區(qū)兩端產(chǎn)生MPa量級的ΔCFS增量。其中,1955年康定地震在龍門山斷裂帶南部產(chǎn)生了顯著的應(yīng)力卸載效應(yīng),最大可達(dá)到-0.8 MPa(圖7(d))。此后,鮮水河斷裂帶上1973年和1981年2次強震持續(xù)卸載了龍門山斷裂帶南端的庫侖應(yīng)力(圖7(e)、7(f))。2008年汶川地震導(dǎo)致破裂區(qū)南北端ΔCFS加載約0.03~0.07 MPa,在龍門山斷裂帶南部,即2013年蘆山地震附近也存在應(yīng)力加載效應(yīng)(圖7(g)),表明其對蘆山地震的觸發(fā)作用。在2022年瀘定地震發(fā)生前,9次歷史地震在破裂區(qū)產(chǎn)生的總ΔCFS約為0.07 MPa(圖7(j))。安寧河-則木河斷裂帶、理塘斷裂帶、麗江-小金河斷裂帶、玉龍希斷裂帶南段、大涼山斷裂帶等ΔCFS變化不大(圖7(j)),表明其受9次歷史地震的影響相對較小。

(a)~(i)為自1893年以來9次大地震依次破裂造成的累積ΔCFS;(j)為2022年瀘定地震前各斷層上ΔCFS;黑色線條為模擬地震的破裂位置

2.3.2 斷層上總庫侖應(yīng)力演化

將震間構(gòu)造加載與歷史地震產(chǎn)生的應(yīng)力變化相結(jié)合,得到川滇塊體東邊界區(qū)域主要斷裂帶的總ΔCFS演化結(jié)果(圖8)。與單獨考慮歷史地震產(chǎn)生的應(yīng)力擾動效應(yīng)相比(圖7(f)),2008年汶川地震震源區(qū)位于歷史地震產(chǎn)生的應(yīng)力影區(qū),而考慮構(gòu)造應(yīng)力加載效應(yīng),震源區(qū)的應(yīng)力積累狀態(tài)在地震前是正值(圖8(f))。在鮮水河斷裂等應(yīng)力積累較高的地區(qū),震間構(gòu)造應(yīng)力積累有很大的貢獻(xiàn)。從圖8(b)~8(e)中可以看出,與強震的影響相比,1923年M7.3仁達(dá)、1955年M7.5康定、1973年M7.6爐霍地震的破裂段庫侖應(yīng)力有顯著增長,這3次強震的破裂中心的庫侖應(yīng)力增加主要依賴于構(gòu)造加載作用。但前序強震造成破裂段兩側(cè)庫侖應(yīng)力增加,從而導(dǎo)致的觸發(fā)效應(yīng)也十分顯著,比如1904年M7.0道孚地震和1981年M6.9道孚地震皆發(fā)生在前序破裂導(dǎo)致的庫侖應(yīng)力增長區(qū)域。此外,在圖7(j)中,安寧河-則木河斷裂帶、麗江-小金河斷裂帶北段、理塘斷裂帶以及大涼山斷裂帶應(yīng)力變化不明顯,但它們在圖8(j)中表現(xiàn)出較高的應(yīng)力積累量,分別約為0.04~0.08 MPa、0.05 MPa、0.12~0.19 MPa、0.10~0.15 MPa,說明這些區(qū)域的應(yīng)力積累主要受構(gòu)造應(yīng)力加載作用。其中,鮮水河斷裂帶、理塘斷裂帶北段、大涼山斷裂帶南段在百年時間尺度內(nèi)應(yīng)力積累尤為顯著(>0.1 MPa),可能有較大的地震活動潛力。

(a)~(i)為自1893年以來9次大地震依次破裂產(chǎn)生的總ΔCFS;(j)為2022年瀘定地震前各斷層上總ΔCFS

3 討 論

3.1 瀘定地震的發(fā)震成因

基于以上結(jié)果,進(jìn)一步討論前述9次歷史地震和構(gòu)造應(yīng)力加載作用與2022年瀘定地震之間的關(guān)系。震間構(gòu)造應(yīng)力積累速率的計算結(jié)果表明,瀘定地震震源區(qū)庫侖應(yīng)力加載速率約為0.75 kPa/a(圖5),該結(jié)果與Luo等[8]的計算結(jié)果(0.77~0.91 kPa/a)相近。基于本文計算結(jié)果,進(jìn)一步計算百年時間范圍內(nèi)的構(gòu)造應(yīng)力積累量為0.097 MPa。歷史地震記錄顯示,該地區(qū)上次強震發(fā)生于1786年,距離2022年瀘定地震已有236 a,假設(shè)震間構(gòu)造應(yīng)力積累速率保持穩(wěn)定,在瀘定地震復(fù)發(fā)周期內(nèi)應(yīng)力積累量約為0.177 MPa。前序9次強震累積產(chǎn)生庫侖應(yīng)力增加量約為0.07 MPa(圖7(j)和圖9),其中1955年康定地震與破裂區(qū)距離較近,對其造成約0.047 MPa的庫侖應(yīng)力加載效應(yīng)。歷史地震在瀘定地區(qū)產(chǎn)生的應(yīng)力擾動量達(dá)到了理論地震觸發(fā)閾值(0.01 MPa),說明附近歷史地震對此次瀘定地震起到了促進(jìn)作用。基于本文計算結(jié)果認(rèn)為,此次瀘定地震的發(fā)生主要受構(gòu)造應(yīng)力加載作用驅(qū)動,但周圍歷史強震也對其有不可忽視的影響。

圖9 瀘定地震震源區(qū)庫侖應(yīng)力變化

3.2 區(qū)域未來地震危險性

鮮水河斷裂帶上地震活躍,而安寧河-則木河斷裂帶相對平靜,其中安寧河段M≥7歷史地震的離逝時間已有487 a,則木河段為170 a[1]。前人根據(jù)地震活動圖像、小震精定位信息及其他形變測量數(shù)據(jù)等綜合分析認(rèn)為,安寧河-則木河斷裂帶存在2個典型的地震空區(qū)(圖1中紅色虛線橢圓圍限區(qū)域),位于27.5°~28.9°N(冕寧-西昌段),長度約為140 km,估計潛在地震震級為7.4級[1]。汶川地震后,GNSS觀測數(shù)據(jù)揭示安寧河-則木河斷裂帶和大涼山斷裂帶中南段完全閉鎖,在安寧河-則木河地震空區(qū)段存在MW6.8~7.5的強震風(fēng)險[19-20]。基于本文計算的構(gòu)造應(yīng)力增長速率(約0.5~0.7 kPa/a),在空區(qū)1(安寧河斷裂)487 a構(gòu)造應(yīng)力積累量約為300 kPa,空區(qū)2(則木河斷裂)170 a構(gòu)造應(yīng)力積累量約為87 kPa。歷史地震在這2個空區(qū)產(chǎn)生的應(yīng)力積累量較小,約為數(shù)個kPa。在地震復(fù)發(fā)周期內(nèi)的總應(yīng)力積累量分別為310 kPa和83 kPa。未來應(yīng)進(jìn)一步關(guān)注該區(qū)域發(fā)生地震的可能性,尤其在安寧河斷裂空區(qū)段可能具有較大的地震危險性。

盡管大涼山斷裂帶受歷史地震應(yīng)力擾動效應(yīng)不明顯,且其受到2008年汶川地震和2013年蘆山地震的微弱卸載效應(yīng)(圖7),然而據(jù)本文計算,大涼山南北段震間平均構(gòu)造應(yīng)力增長速率分別約為0.6 kPa/a和1.0 kPa/a,其中南段最近一次歷史強震離逝時間已有千年[21],構(gòu)造應(yīng)力積累量約為1.0 MPa。汶川地震后的GNSS觀測數(shù)據(jù)反演結(jié)果表明,該段完全閉鎖,有MW7.5強震風(fēng)險[20],其較高的構(gòu)造應(yīng)力積累速率表明該斷裂帶的地震危險性值得關(guān)注。

4 結(jié) 語

綜合考慮地質(zhì)、地球物理等基礎(chǔ)資料,建立川滇塊體東邊界三維粘彈性有限元模型,以GNSS速度場為模型約束邊界條件,計算不同斷裂帶構(gòu)造應(yīng)力積累速率,并同時考慮周邊1893年以來9次歷史強震造成的同震和震后庫侖應(yīng)力變化,綜合分析研究區(qū)域各主要斷裂帶應(yīng)力演化特征。結(jié)果表明:

1)鮮水河斷裂帶、理塘斷裂帶和大涼山斷裂帶南段的應(yīng)力積累速率較高,可達(dá)到1.0~1.5 kPa/a;安寧河斷裂帶、則木河斷裂帶和大涼山斷裂帶北段應(yīng)力積累速率次之,約為0.5~1.0 kPa/a;龍門山斷裂帶、玉龍希斷裂帶和麗江-小金河斷裂帶震間庫侖應(yīng)力積累速率低于0.5 kPa/a。

2)瀘定地震除造成鮮水河斷裂帶震源破裂區(qū)附近庫侖應(yīng)力顯著增加外,同時造成麗江-小金河斷裂帶以及大涼山斷裂帶北端0.01~0.1 MPa量級的應(yīng)力加載效應(yīng);與鮮水河斷裂帶交界處的龍門山斷裂帶西南段庫侖應(yīng)力減小約0.5 MPa。瀘定地震震源區(qū)在地震復(fù)發(fā)周期內(nèi)積累的構(gòu)造應(yīng)力加載約為0.177 MPa,9次歷史地震產(chǎn)生的應(yīng)力擾動效應(yīng)約為0.07 MPa。據(jù)此認(rèn)為,瀘定地震的發(fā)生主要受控于構(gòu)造應(yīng)力的加載驅(qū)動,同時也受到周邊歷史強震的觸發(fā)效應(yīng)影響。

3)青藏高原東南緣主要斷裂帶的總庫侖應(yīng)力演化結(jié)果表明,鮮水河斷裂帶應(yīng)力積累速率高,震間應(yīng)力積累效應(yīng)顯著,但其中一些強震破裂區(qū)的庫侖應(yīng)力變化也很大程度上受同震或震后的影響。安寧河-則木河斷裂帶、麗江-小金河斷裂帶北段、理塘斷裂帶、大涼山斷裂帶以及玉龍希斷裂帶南段主要受構(gòu)造應(yīng)力加載作用,其中安寧河斷裂帶地震空區(qū)在前次地震離逝時間(約487 a)內(nèi)的總應(yīng)力積累量約為0.3 MPa,大涼山斷裂帶南段和理塘斷裂帶歷史地震離逝時間均為千年左右,總應(yīng)力積累量分別約為1.0 MPa和1.0~1.5 MPa,可能有較大的地震危險性。

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 91精品综合| 九九这里只有精品视频| 成人一区专区在线观看| 成人av手机在线观看| 国产综合亚洲欧洲区精品无码| 欧美综合区自拍亚洲综合天堂| 在线看免费无码av天堂的| 97视频免费看| 免费大黄网站在线观看| 无码福利视频| 成人精品午夜福利在线播放| 无码专区在线观看| 亚洲中文字幕在线观看| 亚洲国产欧美自拍| 天天做天天爱夜夜爽毛片毛片| 99久久国产综合精品女同| 亚洲第一黄色网| 伊人福利视频| 在线人成精品免费视频| 欧美亚洲国产视频| 国产成人精品一区二区不卡| 天天综合网亚洲网站| 中文字幕亚洲无线码一区女同| 午夜成人在线视频| 黄色网在线| 青青热久免费精品视频6| 婷婷亚洲天堂| 韩日午夜在线资源一区二区| 久久www视频| 日本午夜精品一本在线观看| 在线日韩一区二区| 91高清在线视频| 久久国产精品波多野结衣| 色天天综合| 亚洲成人高清无码| 国产拍揄自揄精品视频网站| 福利在线免费视频| 精品国产电影久久九九| 91青青草视频| 婷婷色一二三区波多野衣| 91免费国产在线观看尤物| 又爽又大又光又色的午夜视频| 免费国产好深啊好涨好硬视频| 四虎影视库国产精品一区| 婷婷午夜天| 国产美女久久久久不卡| 亚洲天堂啪啪| 任我操在线视频| 亚洲欧美日韩动漫| 91精品视频网站| 1769国产精品免费视频| 99在线视频免费| 欧美日韩综合网| 中国美女**毛片录像在线| 国产一区二区三区精品欧美日韩| 日本亚洲国产一区二区三区| 丁香五月婷婷激情基地| 97se亚洲综合在线| 日韩不卡免费视频| 狠狠色香婷婷久久亚洲精品| 老司机午夜精品视频你懂的| 国产a v无码专区亚洲av| 波多野结衣一区二区三区四区视频| 99re免费视频| 亚洲中文字幕国产av| 午夜日b视频| 国内精品一区二区在线观看| 亚洲狠狠婷婷综合久久久久| 免费视频在线2021入口| 亚洲国产精品成人久久综合影院 | 亚洲色图欧美视频| 午夜视频www| 就去色综合| 免费a在线观看播放| 一级做a爰片久久免费| 91青青在线视频| 四虎国产永久在线观看| 精品夜恋影院亚洲欧洲| a毛片免费在线观看| 嫩草影院在线观看精品视频| 国产真实自在自线免费精品| 亚洲欧美在线综合图区|