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

2015年尼泊爾MW7.9地震對強余震的觸發作用和對周邊斷層的影響冰

2016-08-27 11:06:22黃祿淵楊樹新陳連旺李玉江米琦
地震研究 2016年2期
關鍵詞:余震模型

黃祿淵 楊樹新 陳連旺 李玉江 米琦

摘要:基于EDGRN/EDCMP程序,計算了尼泊爾MW7.9地震引起的庫侖應力變化對強余震的觸發作用和對周邊斷層的加/卸載作用。結果顯示:4月25日MW6.7和4月26日MW6.8強余震均受到明顯觸發作用,地震產生的應力場擾動有利于2015年5月12日MW7.2強余震的發生。單獨反演破裂模型和聯合反演破裂模型的計算結果均顯示:余震的分布和庫侖應力變化的正值區有很好的對應關系,分別有85.5%和75.8%的余震落在庫侖應力變化正值區。雅魯藏布江縫合帶西段的西側受到超過10kPa的庫侖應力加載。喀喇昆侖斷裂、班公錯怒江縫合帶、格林錯斷裂的庫侖應力加載約為1kPa。龍木錯-邦達錯斷裂、瑪尼玉樹斷裂以及嘉黎斷裂受本次地震的影響較小。在震中的西北和東北方向,藏南典型南北向伸展構造受到了庫侖應力加載作用,量值約為1kPa。

關鍵詞:尼泊爾MW7.9地震;庫侖應力變化;應力觸發;周邊斷層

中圖分類號:P315.72 文獻標識碼:A 文章編號:1000-0666(2016)02-0187-09

0 引言

2015年4月25日14時11分26.3秒,在尼泊爾境內發生Mw7.9強烈地震,震中位置(28.2°N,84.7°E),震源深度為15km。根據中國地震局地球物理研究所(2015)發布的快速反演結果,斷層面上最大滑動量約2.1m,而根據Hayes(2015)的結果,斷層面上最大滑動量約為3.1m。張勇等(2015)后續又增加了GSN全球臺網中震中距更近的資料和2個GPS臺站的同震位移數據,進行了地震波資料單獨反演與聯合反演,得到的斷層面上的最大滑動量分別為3.2m和5.2m。4月26日我國西藏自治區日喀則地區的吉隆縣、定日縣以及聶拉木縣發生了多次地震,其中最大震級為MS5.9,給日喀則地區造成了一定的經濟損失。4月25日14時45分和4月26日15時9分尼泊爾陸續發生了MW6.7和MW6.8強余震,5月12日在尼泊爾主震的西南部又發生了MW7.2強余震。

2015年尼泊爾MW7.9地震發生于青藏高原南部喜馬拉雅構造帶上,該帶在印度板塊與歐亞板塊的陸陸俯沖碰撞和推擠下,發生過多次大地震,如1934年BiharMS8.2地震(Sapkota et al,2013),1950年察隅Ms8.7地震(Gupta,Gahalaut,2014)和2005年KashmirMS7.6地震(Avouac et al,2006)。印度板塊相對于歐亞板塊以約46mm/a的速度向北運動,使得青藏高原內部產生南北向的縮短和東西向的拉張,在喜馬拉雅構造帶上應力高度集中,尼泊爾MS7.9地震是由歐亞板塊與印度板塊持續的南北向俯沖擠壓造成的,深入研究該地震鄰近區域及周圍斷層的庫侖應力加,卸載作用,對于理解青藏高原的動力學模型以及青藏高原內部的地震危險性分析具有重要意義。

目前地震應力觸發與斷層相互作用的研究多以庫侖破裂應力變化為基礎。Stein和Lisowski(1983)首次提出庫侖應力概念,發現1979年加州Homestead Valley大地震的余震集中在庫侖應力大于3bar的區域。這類研究的廣泛開展始于1992年美國Landers MS7.3地震,King等(1994)計算Landers地震對后續Big Bear地震產生的庫侖應力增量為2.2~2.9bar,認為Landers地震加速了BigBear地震發生。Jaumé和sykes(1996)計算發現舊金山灣地區的許多斷層位于1906年舊金山大地震產生的應力影區,使下一次大地震的發生延遲數十年。萬永革等(2000,2002)討論我國幾次復雜地震中的應力觸發問題,并對地震靜態應力觸發模型做了全球檢驗。沈正康等(2003)討論了東昆侖活動斷裂帶大地震之間的粘彈性應力觸發。汶川地震后,許多學者(Parsons et al,2008;萬永革等,2009;單斌等,2009;石耀霖,曹建玲,2010;李玉江等,2013)給出了汶川地震對周圍斷層活動的影響。

Okada(1992)的彈性半空間位錯解析解,具有數學上的解析性和簡潔性,但只能計算半無限空間均勻介質的同震形變。Wang等(2003)發展的EDGRN/EDCMP程序,可以計算考慮自重的層狀均勻介質中點源位錯產生的位移的格林函數,把破裂面離散成許多離散的點位錯,通過線性疊加的方法計算同震形變。本文采用EDGRN/EDCMP程序,計算了尼泊爾MS7.9地震在后續3次強余震破裂面上的靜態庫侖應力變化,并定量分析本次地震對周邊活動斷裂帶的應力加/卸載作用,從而判斷尼泊爾MW7.9地震對周圍斷層地震活動性的影響。

1 研究區域及數據

1.1 研究區域與構造背景

青藏高原南緣的喜馬拉雅構造帶是印度板塊與歐亞板塊的俯沖、碰撞邊界,這一過程開始于距今約55~65Ma,至今仍處于強烈推擠過程中。喜馬拉雅構造帶由MCT(主中央斷裂)、MBT(主邊界斷裂)和MFT(主前緣斷裂)3個巨大逆斷裂系組成,其南部最新活動邊界的MFT活動性最強(鄧起東等,2014),此次地震發生在MFT或者鄰近的位置(USGS,2015)。青藏高原次級塊體拉薩塊體、羌塘塊體距離此次地震震中位置較近,塊體之間被活動斷裂系所分割,歷史記錄顯示在這些活動斷裂上發生了多次7級以上地震。由于Okada(1992)彈性位錯表達式基于平面假設,未能考慮地球曲率影響,因此選擇的研究區域不宜過大,本文的研究區域和區域內主要斷層見圖1。

1.2 震源模型及計算參數

尼泊爾地震發生后的2小時9分,張勇等(2015)利用遠震地震波數據快速確定和發布了此次地震的破裂過程模型,此后又根據更豐富的資料,增加了GSN全球臺網中震中距更近的數據以及同震位移明顯的兩個臺站的GPS數據對破裂模型進行了更新(圖2)。聯合反演得到的靜態滑動量分布比地震波單獨反演結果的滑動分布范圍有所縮小,但滑動量幅度增加,最大滑動量由地震波單獨反演的3.2m提高到聯合反演的5.2m。此外,Hayes(2015)利用波形反演也給出有限斷層模型,但其使用的是遠震的數據,考慮到近場信息對于確定震源破裂過程意義重大,本文的計算采用張勇等(2015)用單獨地震波反演和地震波、GPS聯合反演得到的兩個破裂模型。

本文的分層地殼模型參數參考自Crust2.0(Bassin et al,2000)給出的研究區沉積層、上地殼、中地殼和下地殼分層的地震波速和密度數據,參數如表1所示。

2 庫侖應力變化計算原理

根據庫侖破裂準則,巖石趨于破裂時庫侖破裂應力增量為(Harris,1998)

△σf=△τ+μ'△n. (1)式中,μ'為視摩擦系數,△τ為剪應力變化量,△σn為法向正應力變化量。

筆者采用(1)式計算靜態庫侖應力變化,對應力符號規定如下:△τ與破裂滑動方向一致為正,△σn以拉為正。μ'的選取,一般認為滑動速率比較大的走滑斷層、正斷層等具有較小的摩擦系數,常取0.2~0.4,滑動速率較小的逆沖斷層具有較大摩擦系數,常取0.6~0.8(Parsons et al,1999)。地震引起的擾動應力張量的6個獨立分量(△σxx,△σyy,△σzz,△σxy,△σyz,△σxz),可以通過二階張量坐標變換(陳連旺等,2008)計算任意斷層面上△σn,和該斷層面上任意滑動方向的△τ,代入(1)式得到庫侖應力變化△σf。如果△σf為正,認為有利于斷層的滑動或后續地震的發生,△σf為負則認為不利于斷層的滑動或后續地震的發生。

3 尼泊爾MW7.9地震與后續余震的關系

3.1 對3次強余震的觸發作用

2015年4月25日和4月26日在主震震中附近陸續發生了MW6.7和MW6.8強余震,5月12日在尼泊爾主震的西南部又發生了MW7.2強余震。表2給出了這3次強余震的破裂信息,筆者據此計算了余震破裂面上的庫侖應力變化。需要說明的是,大地震發生后的大量余震會產生可以觀測到的地表位移場(萬永革等,2005),因此余震會對周圍區域的應力場產生擾動,所以筆者在計算中考慮了4月25日MW6.7強余震對后續兩余震的影響以及4月26日MW6.8強余震對5月12日MW7.2強余震的影響。計算中,MW6.7和MW6.8強余震的錯動方向按照表2給出,位錯大小、破裂長度、破裂寬度按照經驗公式(Wells,Coppersmith,1994)給出。

表2顯示,2015年4月25日MW6.7和4月26日MW6.8地震破裂面上庫侖應力變化均超過0.01MPa;5月12日MW7.2地震破裂面上的庫侖應力變化均為正值,但獨立反演模型計算的庫侖應力變化大于一般認為的庫侖應力觸發閾值0.01MPa,而聯合反演模型計算的庫侖應力變化低于0.01MPa。表2僅給出了取視摩擦系數為0.6的庫侖應力值,筆者還測試了視摩擦系數分別為0.4和0.8的情況,結果類似。兩種震源模型的計算結果都支持尼泊爾MW7.9地震觸發了4月25日MW6.7和4月26日MW6.8余震,單獨反演模型的計算結果支持MW7.9地震觸發了5月12日MW7.2余震,聯合反演模型的計算結果顯示MW7.9地震產生的應力場擾動有利于2015年5月12日MW7.2余震的產生。

3.2 庫侖應力變化與余震分布的關系

尼泊爾MW7.9地震發生后,截至6月3日,共記錄到其周邊發生124次MW>4.0余震,余震的平均震源深度約15km,幾次震級較大的余震震源深度也接近15km,因此本文采用15km作為庫侖應力變化的計算深度。有研究采用最優破裂面來計算庫侖應力變化與余震的對應關系,但最優破裂面的取向取決于區域背景應力場和地震破裂擾動應力場的迭加(King et al,1994),而往往只能通過震源機制解或者實測應力獲得區域背景應力場的方向,對背景應力場大小的獲取也局限于淺表的地應力測量(周龍壽等,2013)。因此,對區域背景應力場的估計往往很粗略,故本文僅采用主震破裂面作為庫侖應力變化的投影面,計算結果見圖3。

圖3a顯示,有85.5%的余震落在庫侖應力變化正值區,其中78%的余震落在變化大于0.01MPa的區域。圖3b顯示,有75.8%的余震落在庫侖應力變化正值區,其中32.3%的余震落在變化大于0.01MPa的區域。兩震源破裂模型均顯示,余震的分布和庫侖應力變化的正值區有很好的對應關系,圖3b所示的庫侖應力變化正值區面積比圖3a小,這可能與聯合反演破裂模型的滑動分布范圍小于獨立反演破裂模型有關。

圖3僅給出了取視摩擦系數為0.6的庫侖應力值,筆者還測試了視摩擦系數分別為0、0.4、0.8的情況,庫侖應力變化與余震分布的對應關系沒有發生大的改變,例如,采用圖3a的震源破裂模型取視摩擦系數為0和0.8進行計算,分別有81.5%和80.7%的余震落在變化正值區,其中分別有75.8%和76.6%的余震落在變化大于0.01MPa的區域。

4 尼泊爾MW7.9地震對周邊斷層的影響

為研究尼泊爾MW7.9地震的發生對附近斷裂構造的影響,選取圖1中的8條主要活動斷裂,計算尼泊爾MW7.9地震震源破裂在這些斷層面上產生的庫侖應力變化。斷層面的走向、傾向、滑動角規定了庫侖應力變化計算中的“接收斷層”,通過對震源機制解取平均,并根據參考文獻(萬永革等,2010;馬杏垣,1989;彭小龍,王道永,2013;鄧起東,張培震,2002;肖根如等,2010;張培震等,2003)得到具體的斷層參數,見表3。雅魯藏布江斷裂沿印度與歐亞板塊碰撞的縫合帶發育,該區海拔高、生存環境惡劣,并且該區可用的地震資料稀少、研究程度很低,尹安(2001)在雅魯藏布江縫合帶發現了大量逆沖推覆構造及其相伴產生的韌性剪切變形帶,但張培震等(2003)根據橫跨青藏高原南部的GPS觀測發現了該斷裂活動的跡象,其運動方式為右旋,滑動速率為(5±3)mm/a,因此本文中其滑動角范圍為90°~180°。

按照表3各條斷層的走向、傾向和滑動角以及摩擦系數范圍計算了各條斷層0~20km深度范圍的庫侖應力變化,以深度范圍里的極值(深度的采樣間隔為2km)來表示斷層應力狀態的改變,圖4僅給出各斷層取各自摩擦系數范圍中間值的結果,其中雅魯藏布江縫合帶滑動角取中間值135°。為了考察不同滑動角、摩擦系數帶來的不確定性,圖5給出了單獨反演模型計算的雅魯藏布江縫合帶西段斷層面在不同滑動角和摩擦系數條件下的庫侖應力變化,可以看出不同滑動角、摩擦系數會影響斷層上的庫侖應力正值區的面積大小。

由圖4可以看出,2種震源模型計算的斷層面上的庫侖應力變化整體特征較為一致。雅魯藏布江縫合帶西段的庫侖應力變化正負交替,其西側的最大庫侖應力加載超過了10kPa,其東側的最大庫侖應力卸載約為0.1kPa。龍木錯-邦達錯斷裂西段的西側受到輕微庫侖應力卸載作用,龍木錯-邦達錯斷裂西段的東側、龍木錯-邦達錯斷裂東段、瑪尼玉樹斷裂的庫侖應力加載約為0.1kPa。喀喇昆侖斷裂、班公錯怒江縫合帶、格林錯斷裂的庫侖應力加載約為1kPa。嘉黎斷裂的庫侖應力變化正負交替,加/卸載量值均小于0.1kPa,受本次地震的影響較小。

青藏高原一個最顯著的地質特征是在藏南及高原腹地廣泛發育的南北向伸展構造如藏南拆離系(STDS)。在先前的討論中,表3并未給出藏南內部南北向的正斷層,根據張進江和丁林(2003)獲取了藏南南北向伸展構造的大致產狀,作為簡化計算了走向北東,傾角、滑動角的典型正斷層的庫侖應力變化,以此來代表本次地震對藏南的南北向伸展構造的影響,結果見圖6。可以看出,在尼泊爾Mw7.9地震的西北和東北方向,藏南的南北向伸展構造受到了庫侖應力加載作用,量值約為1kPa。

5 不確定性討論

現有的觀測資料還不足以對靜態庫侖應力變化計算提供足夠的約束,如介質的縱向分層和橫向的不均勻性、初始應力場的未知、震源破裂過程反演的不確定性、孔隙水壓力的差異以及斷層摩擦性質的差異。單斌等(2009)討論了斷層參數不確定性較高的岷江斷裂取不同傾角和滑動角對計算的影響,石耀霖和曹建玲(2010)采用修正的庫侖應力計算方法,比較了不同地震破裂模型和視摩擦系數條件下的庫侖應力分布差異。王力偉和陳棋福(2010)分析了不同計算模型、斷層位置、走向、傾角、摩擦系數等因素對庫侖應力變化計算的影響。

地震產生的庫侖應力變化強烈依賴于震源模型,如圖3a、b所示,在近場庫侖應力變化正值區與余震的對應程度存在一定差異。此外,Hayes(2015)在震后1小時發布了震源破裂模型,該模型與張勇等(2015)模型的方法和使用的數據都不一樣,筆者在相同的程序和相同的摩擦系數和分層介質參數條件下計算了Hayes模型的庫侖應力變化,發現僅49.2%的余震落在庫侖應力變化正值區。由于沒有足夠的同震形變資料,也沒有直接的標準判斷不同破裂模型的優劣,如果僅以應力變化與后續地震直接的對應性作為標準,張勇等(2015)利用地震波獨立單獨反演的破裂模型更優。

石耀霖和曹建玲(2010)發現不同的視摩擦系數雖然有一定影響,但不會對庫侖應力整體分布圖像有太大影響。視摩擦系數的變化改變了正應力變化在庫侖應力變化計算中的權重,對于斷層面上的庫侖應力變化,按照表3中各斷層的摩擦系數變化范圍做了測試,發現視摩擦系數的影響主要體現在庫力應力量值的變化,也可能會導致庫侖應力正值區面積的變化,如圖5所示,隨著視摩擦系數的增加,斷層上庫侖應力正值區面積增加。此外,滑動角的不確定性對庫侖應力的變化影響也較大,對于雅魯藏布江縫合帶西段,當滑動角為135。時庫侖應力正值區面積最大。因此,想要獲得對地震活動性有指示意義的結果,需要結合更多資料縮小斷層參數和視摩擦系數的不確定范圍。

研究區域包含了拉薩、羌塘和板塊邊界等構造塊體,EDGRN/EDCMP軟件只能計算橫向均勻、縱向分層模型的格林函數,因此比較不同地殼模型對庫侖應力變化的影響,根據Crust2.0給出的拉薩塊體、羌塘塊體和板塊邊界的地殼模型參數,計算不同的分層地殼模型的同震形變。采用不同地殼模型,結果的差異主要體現在庫侖應力變化的量值,例如對于雅魯藏布江縫合帶,不同地殼模型的計算結果的平均相對變化量不超過12.8%。此外,由圖5可以看出,地殼模型差異對計算結果的影響較斷層面上滑動角的不確定性產生的影響要小,因此作為一級近似,本文采用橫向均勻、縱向分層的模型進行計算是可以接受的。需要指出的是,真實地殼具有很強的橫向不均勻性,地表和莫霍面等圈層也具有較大的起伏,需要在將來采用數值方法建立更為接近真實的模型,更精確地計算庫侖應力變化。

6 結論與討論

本文根據地震靜態觸發原理,利用已發布的有限斷層模型,計算了尼泊爾MW7.9地震引起的同震庫侖應力變化對余震的觸發作用,并分析了對鄰近主要活動斷層的影響。主要結論如下:

(1)尼泊爾MW7.9地震對4月25日MW6.7和4月26日的MW6.8強余震有顯著觸發作用,MW7.9地震產生的應力場擾動有利于5月12日MW7.2余震的產生。

(2)單獨反演和聯合反演兩個震源破裂模型均顯示,余震的分布和庫侖應力變化的正值區有很好的對應關系,分別有85.5%和75.8%的余震落在庫侖應力變化正值區。

(3)雅魯藏布江縫合帶西段的西側受到超過10kPa的庫侖應力加載,地震活動性增強。喀喇昆侖斷裂、班公錯怒江縫合帶、格林錯斷裂的庫侖應力加載約為1kPa。龍木錯-邦達錯斷裂、瑪尼玉樹斷裂以及嘉黎斷裂的庫侖應力加/卸載小于0.1kPa,受本次地震的影響較小。

(4)在尼泊爾MW7.9地震的西北和東北方向,藏南典型南北向伸展構造受到了庫侖應力加載作用,量值約為1kPa。

震源滑動模型對同震應力變化的計算影響較大,尤其在離震源較近的近場區域。本文采用的兩個震源破裂模型在近場的庫侖應力變化有較大差別,但遠場斷層面上庫侖應力變化的整體特征比較一致。此外,Hayes(2015)利用遠場地震波數據反演得到了一個震源破裂模型,由于缺乏近場的地震波數據約束,使用該模型在實際計算中得到的近場庫侖應力變化與余震對應的并不好,僅有不超過半數的余震落在庫侖應力變化正值區。特別地,本文使用的兩個震源模型來源于同一作者,它們的差別僅在于反演過程中是否考慮GPS數據,因此只是在滑動量的幅值和分布面積上有一些差別。但每次大地震發生后,不同作者給出的破裂模型往往具有很大差異,如2010年Men-tawai MW7.8地震反演工作中,不同研究得到的最大滑動量可介于3.5~20m之間(Lay et al,2011;Yue et al,2014)。因此,為了提高計算結果的可靠性,將來需要使用基于更多近場地震波數據和大地測量數據反演得到的斷層破裂模型。

斷層的視摩擦系數、滑動角具有較大的不確定性,在地震活動性分析中,應該重視參數的不確定性引起的庫侖應力變化正值區的改變。本文采用了橫向均勻、縱向分層模型,真實地球介質的橫向不均勻性對庫侖破裂應力變化的計算結果有一定影響,需要在未來利用數值手段建立更接近真實的模型來計算庫侖應力變化。此外,在研究庫侖應力變化與余震分布的對應關系時,將來建立在余震精定位和更完整的地震序列基礎上的分析工作將更具可靠性。

受限于目前還沒有可以利用的大量余震的震源破裂信息,本文僅考慮了二次強余震對庫侖應力變化的影響,大量的余震活動對該區域的應力變化產生多大影響?根據萬永革等(2005)的研究結果,Landers地震后發生大量余震活動,產生的位移場方向與主震大體一致達到厘米量級,不足主震產生位移場的10%。據此可推論,大量余震對庫侖應力變化的空間模式改變不會很大,但對絕對量值會有一定改變。下地殼、上地幔的黏彈性松弛效應,也對庫侖破裂應力變化結果產生一些影響,但本文討論的是地震發生后短期內的彈性響應,時間尺度遠低于流變介質松弛時間。另外,本文未考慮背景應力場,動態應力觸發,流體作用等因素,相關內容需繼續深入研究,本文的研究僅從一個角度提供一個有意義的參考。

猜你喜歡
余震模型
一半模型
基于指數函數的川滇地區余震序列衰減規律研究
“超長待機”的余震
哈哈畫報(2022年5期)2022-07-11 05:57:48
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
生死之間的靈魂救贖——《余震》和《云中記》的倫理問題
阿來研究(2019年2期)2019-03-03 13:35:00
3D打印中的模型分割與打包
本土化改編與再創——從小說《余震》到電影《唐山大地震》
三次8級以上大地震的余震活動特征分析*
地震研究(2015年4期)2015-12-25 05:33:44
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 无码国内精品人妻少妇蜜桃视频 | 亚洲成人高清在线观看| 国产视频a| 国产成人一区免费观看| 亚洲精品无码高潮喷水A| 五月天综合网亚洲综合天堂网| 国产99视频在线| 国产波多野结衣中文在线播放| 久久综合伊人77777| 色哟哟国产精品一区二区| 国产欧美日韩在线一区| 99无码中文字幕视频| 久久久久无码精品| 亚洲午夜片| 看av免费毛片手机播放| 亚洲成人网在线观看| 亚洲国产欧美中日韩成人综合视频| 欧美精品在线免费| 久久精品中文无码资源站| 亚洲欧美精品一中文字幕| 伊人久久综在合线亚洲2019| 成人福利在线看| 成年女人a毛片免费视频| 精品国产一区91在线| 中字无码av在线电影| 在线色国产| 国产av剧情无码精品色午夜| 亚洲人在线| 中文字幕佐山爱一区二区免费| 国产网站免费| 中文字幕啪啪| аⅴ资源中文在线天堂| 国产一区二区三区夜色| 国产99精品视频| 中文字幕日韩视频欧美一区| 亚洲国产综合精品一区| 欧美一区二区福利视频| 97综合久久| 日韩经典精品无码一区二区| 999国内精品久久免费视频| 精品福利视频导航| 四虎永久免费地址在线网站| 亚洲男人的天堂在线观看| 欧美精品一二三区| 青草视频久久| 蜜臀AVWWW国产天堂| 色综合天天娱乐综合网| 午夜日韩久久影院| 中文字幕第1页在线播| 全部免费特黄特色大片视频| 久久大香香蕉国产免费网站| av在线无码浏览| 91精品国产丝袜| 久久夜色精品| 色偷偷一区二区三区| 911亚洲精品| 国产经典免费播放视频| 操美女免费网站| 伊人久久久久久久| 亚洲欧美成人影院| 国产一二视频| 草草影院国产第一页| 日韩在线第三页| 在线观看国产网址你懂的| 欧美黄网站免费观看| 精品国产电影久久九九| 日韩av在线直播| 日韩中文无码av超清| 91免费国产高清观看| 国产极品美女在线| 久久96热在精品国产高清| 91在线精品麻豆欧美在线| 国产综合无码一区二区色蜜蜜| AV在线天堂进入| 国产性生交xxxxx免费| 欧美精品伊人久久| 九色免费视频| 无码中文字幕精品推荐| a天堂视频| 久久精品波多野结衣| 色婷婷成人网| 色首页AV在线|