亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種基于有限動態(tài)源的烈度估計方法

        2016-11-08 03:03:33許力生張旭魏強李春來
        地球物理學報 2016年10期
        關鍵詞:模型

        許力生, 張旭, 魏強, 李春來

        中國地震局地球物理研究所, 北京 100081

        ?

        一種基于有限動態(tài)源的烈度估計方法

        許力生, 張旭, 魏強, 李春來

        中國地震局地球物理研究所, 北京100081

        首先,基于地震烈度與震級和震中距或震源距的經(jīng)驗關系,考慮實際地震斷層的有限性以及幾何學和運動學特性,提出了一種基于有限動態(tài)源模型的烈度預估方法.然后,利用數(shù)值實驗討論了體現(xiàn)幾何學和運動學特性的參數(shù)的作用以及斷層的傾角和破裂速度的影響;最后,將這種方法應用于近年來發(fā)生的兩次破壞性地震,討論了這種方法的實用性與不足之處.在實時獲取破裂過程圖像成為可能的今天,用這種方法預估烈度對于震后救援十分重要.關鍵詞地震烈度; 有限動態(tài)源; 估計方法; 數(shù)值試驗; 實際應用

        1 引言

        傳統(tǒng)的地震烈度已廣泛應用于地震動空間特征和震害空間特征的定性描述. 地震烈度的概念誕生于現(xiàn)代地震儀出現(xiàn)之前,但卻提供了一種描述復雜地面運動的簡單方法.更有意義的是,在許多情況下烈度還成為定量描述災害性地震地震動水平的唯一參數(shù)(Wald et al., 1999a, 1999b).

        雖然地震烈度起初屬于定性參數(shù),但試圖建立其與地面加速度和速度這樣的定量參數(shù)的關系的工作早已開始(Gutenberg and Richter, 1942, 1956; Murphy and O′Brien, 1977).然而,由于地震烈度的原始定義,實際上很難建立其與地面加速度或者速度這些定量參數(shù)的簡單關系,所以迄今為止類似的工作還在繼續(xù)(李大華和左慧強,1991;Wald et al., 1999a; 陳鯤等,2014).況且,地震烈度還具有明顯的區(qū)域特征(Wald et al., 1999a; 汪素云等,2000;Karim and Yamazaki, 2002; Wu et al., 2003; Atkinson and Kaka, 2007; 俞言祥等,2013).

        建立烈度和地震動參數(shù)之間的定量關系的目的大體有兩個,一是通過烈度來計算地震動參數(shù),這對那些沒有地震動參數(shù)的地區(qū)尤其必要;二是通過地震動參數(shù)計算烈度,這對那些沒有標志性震害的地區(qū)尤其重要.前者可以為工程設計服務,后者至少可以為應急救災服務.需要強調(diào)的是,就因為后者,近十多年來出現(xiàn)了一個新的重要的研究方向,即破壞性地震烈度的快速預測或計算(Teng et al., 1997; Wald et al., 1999b; Wu et al., 2001; 李山有等,2002;王玉石等,2008;Moratto et al., 2009).

        到目前為止,基于實測加速度或者速度峰值的預測技術似乎成為快速預測烈度的主流技術(Teng et al., 1997; Wald et al., 1999b; Wu et al., 2001; Moratto et al., 2009).但這種技術需要足夠的強震記錄或地震記錄作前提.顯然,面對沒有足夠記錄設備的地區(qū),這種技術便無能為力.于是,有人開始探索利用震源模型計算地面運動,進而用合成地面運動預測烈度的方法(徐劍俠等,2015),但這種技術似乎離不開高性能的計算設備,且高頻地震動模擬仍很困難.

        面對無記錄設備的地區(qū),一個較好的選擇是利用烈度衰減經(jīng)驗關系.事實上,世界上關于烈度衰減規(guī)律的研究工作已經(jīng)有相當?shù)姆e累(Bormann, 2012; 俞言祥等,2013).但這些工作都是針對地震點源或線源,至少沒有考慮有限震源的動態(tài)破裂過程.我們知道,一個有限震源總是可以考慮成多個點源的集合,因此基于點源的烈度衰減規(guī)律拓廣到有限震源的烈度估計是完全可能的.所以,在這里我們將嘗試一種新方法.這種方法源于點源的烈度衰減經(jīng)驗關系,但充分考慮震源的有限性以及它的幾何學和運動學特性.

        實時確定震源破裂過程的技術已逐漸成熟(Zhang et al., 2014b),因此近實時確定烈度分布也必將成為現(xiàn)實.

        2 烈度公式的推廣

        (1)

        (2)

        (3)

        (4)

        (5)

        IE=3.6588+1.3626MS-3.5406log(R+13),

        (6)

        IX=3.6113+1.4347MS-3.8477log(R+13)

        (7)

        和IQ=3.3682+1.2746MS-3.3119log(R+9)

        (8)分別描述中國東部、新疆地區(qū)和青藏地區(qū)的烈度衰減規(guī)律,而用 IM=3.9440+1.0710MS-2.8450log(R+7)

        (9)描述中強地震區(qū)的烈度衰減規(guī)律.(6)—(9)式中R為震中距.而澳大利亞則采用

        I1=1.64MS-1.70lnR+4.0

        (10)

        I2=1.41MS-1.18lnR-0.0044R+2.18

        (11)

        分別描述板內(nèi)地震和板間地震的烈度衰減規(guī)律(Borman, 2012).注意,(10)和(11)式中R為震源距.換言之,(3)式中的Iij可以根據(jù)具體地域選擇使用.

        我們還知道,面波震級和標量地震矩具有如下經(jīng)驗關系(Lay and Wallace, 1995):

        lg(M0)=1.5MS+9.1,

        (12)

        而矩震級則可以利用標量地震矩M0(N·m)由(13)式確定(Shearer,1999),

        (13)

        所以,點源的烈度衰減關系可以用面波震級來計算,也可以用標量地震矩或矩震級來計算.由此可見,有限動態(tài)源的烈度可以由通過波形反演確定的標量地震矩的空間分布即有限動態(tài)源破裂模型加以確定.

        3 數(shù)值檢驗

        3.1經(jīng)驗常數(shù)的影響

        第2節(jié)給出的烈度估計公式涉及到兩個經(jīng)驗性參數(shù)μ和n,這兩個參數(shù)主要依賴于震源的幾何學和運動學特性,因此其取值應依賴于大量的實例.不過,在還沒有大量應用實例的情況下,通過數(shù)值實驗為其選擇一個合適值或取值范圍是權宜之計.這里我們選用(8)式描述的衰減關系,借助于簡單的線源模型,通過數(shù)值試驗考察這兩個參數(shù)對烈度分布的影響,并為μ和n選擇合適的值或取值范圍.

        選用的線源模型分別如圖1所示的單側破裂和雙側破裂兩種情形,均由15個點源構成,點源間距為1 km,震源深度均為5 km,單側破裂從左向右,雙側破裂從中間位置分別向兩側.為了避免破裂速度的影響,這里設置破裂速度等于S波速度,均為3 km·s-1.為了簡單起見,這里令Mw=7且地震矩平均分配于各點源.

        為了測試n的影響,令μ=0.1,而n=1,2,3,4和5,利用基于(8)式的經(jīng)驗關系得到如圖2所示的結果.可以看出,無論是單側破裂還是雙側破裂,隨著n的增大,烈度在破裂傳播方向的斜前方隨距離的衰減逐漸增大,即n越大衰減越快,可見n主要控制烈度依賴于方位的衰減特性.相互比較,注意到n在2和4之間似乎更符合通常情況下等震線的幾何特征.因此,在沒有足夠?qū)嵗_定n的情況下,不妨取n=3.

        圖1 線源破裂模型黑色箭頭示意破裂方向,紅色五角星表示起始破裂點.破裂速度均為3 km·s-1.(a)單側破裂;(b)雙側破裂.Fig.1 Line-type rupture modelBlack arrows indicate rupture directions, and red stars refer to initiation point. Rupture velocities for both are 3 km·s-1. (a) Unilateral rupture; (b) Bilateral rupture.

        圖2 μ固定n變化時預估烈度的空間分布(a1—a5)對應圖1a所示的單側破裂;(b1—b5)對應圖1b所示的雙側破裂.Fig.2 Spatial distributions of the estimated intensity in case where n is changing while μ is fixed(a1—a5) correspond to the unilateral rupture shown in Fig.1a, and (b1—b5) correspond to the bilateral rupture shown in Fig.1b.

        為了測試μ的影響,令n=3,但μ=0.05,0.1,0.15,0.20和0.25,計算得到如圖3所示的結果.可以看出,無論是單側破裂還是雙側破裂,隨著μ的增大,烈度在破裂的方向上逐漸增強,且烈度區(qū)域逐漸擴展,可見μ主要控制破裂方向?qū)α叶鹊姆糯笞饔?相互比較,注意到μ在0.1與0.2之間的取值似乎更符合通常情況下的等震線幾何特征.因此,在沒有足夠?qū)嵗那闆r下,不妨取μ=0.1.

        由此看來,μ決定破裂方向上烈度的放大作用,而n則調(diào)制烈度在不同方向的衰減特性.

        3.2傾角和破裂速度的影響

        真實的地震斷層不是一個線源,但大多數(shù)情況下可以看成是一個具有一定傾角的面源.傾角不同,斷層面上點源的深度不同,而不同深度的點源在同一場點引起的烈度也不同.同時,斷層的破裂也具有一定的速度,破裂速度的變化必然引起烈度的變化.因此借助于簡單的面源模型、通過數(shù)值試驗認識斷層傾角和破裂速度對烈度的影響是非常必要的.

        為了測試斷層傾角的影響,這里建立如圖4所示的單側破裂和雙側破裂模型,并令μ=0.1和n=3.斷層平面由15個點源構成,分為3層,每層5個點源,最淺層震源深度為2 km, 并讓斷層平面沿最淺層發(fā)生改變,傾角分別為15°、30°、45°、60°和75°.其他參數(shù)同上,計算結果如圖5所示.可以看出,傾角引起的變化主要在極震區(qū),傾角增大使極震區(qū)縮??;隨著距離的增大,傾角的影響逐漸減小.

        圖3 μ變化n固定時預估烈度的空間分布(a1—a5)對應圖1a中的單側破裂;(b1—b5)對應圖1b中的雙側破裂.Fig.3 Spatial distributions of the estimated intensity in case where μ is changing while n is fixed(a1—a5) correspond to the unilateral rupture shown in Fig.1a, and (b1—b5) correspond to the bilateral rupture shown in Fig.1b.

        圖4 面源破裂模型(a) 單側破裂; (b) 雙側破裂(參看圖1).Fig.4 Plane-type rupture model(a) Unilateral rupture; (b) Bilateral rupture (also see Fig.1).

        為了測試破裂速度的影響,固定斷層傾角為45°,同時令μ=0.1和n=3,但破裂速度由小到大發(fā)生變化,分別為Vr=1 km·s-1,2 km·s-1,3 km·s-1,4 km·s-1和5 km·s-1,計算結果如圖6所示.可以看出,破裂速度的影響主要表現(xiàn)在破裂方向上,隨著破裂速度的增加,在破裂方向上烈度增大,烈度區(qū)域相應擴展,這種特征在極震區(qū)尤其明顯.

        由此看來,斷層的傾角和破裂速度對最終烈度的分布都有影響,尤其在極震區(qū).

        圖5 面源傾角對預估烈度的影響(a1—a5)對應傾角分別為15°、30°、45°、60°和75°的單側破裂;(b1—b5)對應傾角分別為15°、30°、45°、60°和75°的雙側破裂.Fig.5 Influence of the dip of source-plane on the estimated intensities(a1—a5) correspond to the unilateral rupture with the dips 15°, 30°, 45°, 60° and 75°, respectively; (b1—b5)correspond to the bilateral rupture with the dips 15°, 30°, 45°, 60° and 75°, respectively.

        圖6 面源破裂速度對預估烈度的影響(a1—a5)對應破裂速度分別為1 km·s-1、2 km·s-1、3 km·s-1、4 km·s-1和5 km·s-1的單側破裂;(b1—b5)對應破裂速度分別為1 km·s-1、2 km·s-1、3 km·s-1、4 km·s-1和5 km·s-1的雙側破裂.Fig.6 Influence of the rupture velocity of plane source on the estimated intensities(a1—a5)correspond to the unilateral rupture with the rupture velocities of 1 km·s-1, 2 km·s-1, 3 km·s-1, 4 km·s-1 and 5 km·s-1, respectively; (b1—b5)correspond to the bilateral rupture with the rupture velocities of 1 km·s-1, 2 km·s-1, 3 km·s-1, 4 km·s-1 and 5 km·s-1, respectively.

        4 對實際震例的應用

        4.1蘆山MS7.0地震

        2013年4月20日在我國四川蘆山發(fā)生了MS7.0地震,這是繼2008年汶川MS8.0地震之后在龍門山斷裂帶發(fā)生的又一次破壞性地震,造成約200人死亡.地震發(fā)生數(shù)日后,中國地震局會同其他機構發(fā)布了這次地震的烈度圖(孟令媛等,2014).為了檢驗本文提出的有限動態(tài)源的烈度預估方法,我們將這種方法應用于這次地震并進行比較討論.

        一個符合實際情況的烈度預估必須基于一個可靠的有限動態(tài)源模型.關于這次地震的震源破裂模型已經(jīng)有多個研究結果(張勇等,2013;王衛(wèi)民等,2013;劉成利等,2013;趙翠萍等,2013;Hao et al., 2013; Zhang et al.,2014a),但由于使用資料和方法的差異,結果之間也有明顯不同.所以,仍有必要利用其他資料和方法進一步確認這次地震的有限動態(tài)源模型.與已有研究不同的是,這里采用從區(qū)域地震記錄(鄭秀芬等,2009)中提取的120條基階Rayleigh波視震源時間函數(shù)以及26個站點的三分向靜態(tài)GPS同震位移作為觀測資料.面波視震源時間函數(shù)比體波視震源時間函數(shù)具有更高的時空分辨能力(Lay and Wallace, 1995),近場GPS資料更有利于約束靜態(tài)位錯分布.另外,這里假設所有的子斷層具有相同的震源機制,大大減少了未知數(shù),有利于解的穩(wěn)定性.

        反演結果如圖7所示,起始破裂點(許力生等,2013b)以南滑動量較大,以西沿斷層面較深的區(qū)域也具有較大的滑動量,表明破裂有向南和向西擴展的優(yōu)勢,最大滑動量達~1.8 m.總體上,各個方向均有破裂擴展,但東北和西北方向相對較弱.需要說明的是,在主要破裂區(qū)之外,位于西北方向的較深區(qū)域還有一些破裂,但這部分破裂在時間上較晚,在空間上較深,在強度上較弱,因此不能十分確認.

        跟早期我們利用體波反演的結果(張勇等,2013)相比,本研究的反演結果中向深部和西南方向破裂的特征更加明顯.早期的結果中最大位錯點基本位于起始破裂位置,而本研究的結果中最大位錯點在起始破裂點的西南且較深的位置.跟我們后來利用遠場和近場加速度資料聯(lián)合反演的結果(Zhang et al., 2014a)相比,本研究的反演結果展示了一個向西南和深部擴展的破裂過程,而聯(lián)合反演的結果則展示了一個向東北和淺部擴展的破裂過程.同樣,王衛(wèi)民等(2013)利用遠場P波和SH波反演結果也展示了一個向東北和淺部擴展的破裂過程.趙翠萍等(2013)利用遠場P波反演的結果卻展示了一個從深部向淺部的破裂過程.劉成利等(2013)利用遠場P波的反演結果展示了一個近乎以起始點為中心的圓盤形破裂.Hao等(2013)利用遠場P波和面波以及近場強震波形反演的結果揭示了和本研究最近似的破裂過程,即沿斷層走向向西南和深部擴展的破裂過程.很顯然,面波資料發(fā)揮了積極的作用.

        圖7 反演得到的蘆山MS7.0地震震源模型與余震分布(a) 蘆山MS7.0地震震源模型的三維展示.五角星表示起始破裂點,粗黑線為上邊界,顏色表示位錯量; (b) 蘆山MS7.0地震等效震源與余震震中分布.五角星表示主震震中位置,顏色表示等效震源的深度.灰色點表示余震震中.Fig.7 The inverted rupture model of the Lushan MS7.0 earthquake and the distribution of the aftershocks(a) The rupture model in 3-D frame. Star refers to the initiation point, thick black line shows the upper bounder of the fault plane, the color indicates the slip amount; (b) The epicenters of the equivalent sources of the Lushan MS7.0 earthquake and its aftershocks. Star refers to the epicenter of the Lushan MS7.0 earthquake, and the color indicates the source depths. Grey dots show the aftershocks.

        關于蘆山MS7.0地震的震源過程,不同的作者利用不同的方法和資料組合得到了總體特征大體相同但細節(jié)上仍存在差異的結果(張勇等,2013;王衛(wèi)民等,2013;劉成利等,2013;趙翠萍等,2013;Hao et al., 2013; Zhang et al., 2014a),但根據(jù)我們多次利用多種資料的反演結果以及其他研究結果,圓盤形破裂為主但在西南方向和沿斷層面深度方向的破裂略占優(yōu)勢的破裂模型應當能夠反映這次地震的主體特征,因此我們基于這樣一個有限動態(tài)源模型預估這次地震的烈度分布.

        需要說明的是,我們還沒有烈度隨震源距的衰減關系,但為了體現(xiàn)震源深度以及斷層傾角對烈度的影響,這里人為地把(8)式中的震中距改為震源距,補充計算了在震源距代替震中距情況下的烈度分布.

        如圖8a所示,如果不考慮震源深度和斷層傾角,極震區(qū)烈度較大,極震區(qū)面積也較大,最大烈度達到IX,很明顯這是由于較深的點源(圖7b中藍色部分)人為抬升到地表所致.如圖8b所示,如果考慮震源深度和斷層傾角,較深的點源的影響明顯減弱,極震區(qū)烈度適度減小,極震區(qū)面積也相應縮小,最大烈度降為VIII+.似乎烈度隨震源距的衰減關系更恰當.

        從圖8b可以看出,預估烈度能夠反映出斷層的有限性和破裂方向性的影響.斷層以西的烈度高于斷層以東且衰減較慢,震中以南的烈度略高于震中以北且衰減較慢,這符合斷層近乎西傾,破裂向南且向下傳播的震源特征.然而,與圖8c所示的實際調(diào)查烈度相比差別仍然明顯.從調(diào)查烈度看,斷層以東的烈度高于斷層以西且衰減較慢,這與預估烈度恰好相反.如果不考慮其他因素,實際調(diào)查烈度特征應該反映一個向東南傾斜且朝西南破裂的有限動態(tài)源模型,而這樣一個模型的傾向和走向恰好與實際的震源模型相反.不過,或許我們不能否認,至少余震的分布更支持本研究確定的烈度分布(圖8).

        4.2魯?shù)镸S6.5地震

        2014年8月3日16時30分(北京時間)在云南省昭通市魯?shù)榭h發(fā)生了MS6.5地震,導致600余人死亡.已有的研究發(fā)現(xiàn),發(fā)震斷層并不在當?shù)氐闹饕獦嬙烨覜]有明顯的地表破裂(徐錫偉等,2014),近場地面運動加速度峰值空間分布不規(guī)則(徐錫偉等,2014;陳鯤等,2015),部分遠場波形的復雜性難以解釋(張勇等,2014),余震非線型分布(房立華等,2014;王未來等,2014;徐甫坤等,2014;張廣偉等,2014).這似乎意味著,這次地震的震源過程比較復雜.

        關于這次地震的破裂過程或震源模型已經(jīng)有多個研究結果(張勇等,2014, 2015;許力生等,2014;劉成利等,2014),但這些結果還不能很好地解釋已有觀測.本研究展示一種新的不同于已有的反演結果,為這次地震提供一種新的震源模型.

        不同于已有的反演研究,我們采用從區(qū)域Love波資料(鄭秀芬等,2009)中提取的視震源時間函數(shù)作為觀測資料,根據(jù)我們最新確定的余震空間分布以及最新確定的震源機制解(許力生等,2014)構建如圖9所示的彎曲斷層模型,并以聯(lián)合地方臺記錄和巧家臺陣記錄借助于非線性方法——逆時成像技術——確定的震源位置(27.089°/103.352°/9.5 km)(許力生等,2014)為起始破裂點.根據(jù)余震定位的震中分布,以間隔約1 km的尺度構建斷層模型與地表的彎曲交線,然后從地表向深部取13層,相鄰兩層間的傾角從90°開始,按1°間隔逐漸變成78°,相鄰兩層沿斷層面向下間距均為1 km.

        圖8 蘆山MS7.0地震預估烈度與調(diào)查烈度的比較(a) 基于震中距的預估烈度; (b) 基于震源距的預估烈度; (c) 調(diào)查烈度.黑色空心圓為余震.Fig.8 Comparison of the investigated intensities with the estimated ones for the Lushan MS7.0 earthquake(a) The forecasted intensities while the epicenter-distance is adopted; (b) The forecasted intensities while the hypocenter-distance is adopted; (c) The investigated intensities. The black open circles are aftershocks.

        需要說明的是,這里使用的主震和余震位置(圖9b)均為我們最新利用逆時成像技術(許力生等,2013a, 2013c),并使用當?shù)氐呐_站與中國地震局地球物理研究所在此布設的巧家流動臺陣的地震記錄聯(lián)合確定的.余震時間跨度為2014年8月1日至2015年3月1日.在這個時間段得以絕對定位的事件2193次,得以雙差定位的事件1284次.我們注意到,無論是絕對定位結果還是雙差定位結果,與以往發(fā)表的結果(房立華等,2014;王未來等,2014)或多或少均有差異.首先,從我們聯(lián)合定位的結果可以直接看出余震的“L”型分布,而以往的結果并非如此;其次,雙差定位后余震均呈“L”型分布,但本文的結果中兩翼的夾角要鈍于以往的結果.

        圖9a展示了反演得到的靜態(tài)位錯分布.可以看出,主要破裂區(qū)位于斷層彎曲及其偏東南的部位和起始破裂點的上方,近乎單側破裂,最大位錯達1.1 m.早期我們利用遠場P波和SH波以及區(qū)域臺網(wǎng)的長周期體波,并且在假設地震發(fā)生在兩個相互垂直的斷層面的情況下反演得到的震源過程表明,地震起始破裂點位于近東西向的斷層,但很快遷移到近南北向的斷層,以南北向斷層的破裂為主(張勇等,2014).劉成利等(2014)利用長周期區(qū)域地震波分別反演了近南北向和近東西向的兩個斷層面的位錯分布,最終認定此次地震發(fā)生在近南北向的斷層上,主要破裂區(qū)位于起始破裂點以南較淺的區(qū)域.雖然我們早期使用了相互垂直的共軛斷層模型(張勇等,2014),而現(xiàn)在使用了彎曲斷層模型;雖然早期使用的資料為遠場體波和區(qū)域體波,而現(xiàn)在使用的資料是勒夫波視震源時間函數(shù),但反演結果一致表明,破裂起始于近東西向的斷層,但很快遷移到近南北向的斷層,而且南北向的斷層為主要的發(fā)震斷層.

        圖9 反演得到的魯?shù)镸S6.5地震震源模型與余震分布(a) 魯?shù)镸S6.5地震震源模型的三維展示; (b) 魯?shù)镸S6.5地震等效震源與余震震中分布(參看圖7).Fig.9 The inverted rupture model of the Ludian MS6.5 earthquake and the distribution of the aftershocks(a) The rupture model in 3-D frame; (b) The epicenters of the equivalent sources of the Ludian MS6.5 earthquake and its aftershocks (also see Fig.7).

        圖10 魯?shù)镸S6.5地震預估烈度與調(diào)查烈度的比較(a) 基于震中距的預估烈度; (b) 基于震源距的預估烈度; (c) 調(diào)查烈度.黑色空心圓為余震.Fig.10 Comparison of the investigated intensities with the estimated ones for the Ludian MS6.5 earthquake(a) The forecasted intensities while the epicenter-distance is adopted; (b) The forecasted intensities while the hypocenter-distance is adopted; (c) The investigated intensities. The black open circles are aftershocks.

        雖然,魯?shù)镸S6.5地震的發(fā)震斷層比較復雜,還有待更深入的研究,但我們認為基于共軛斷層模型和彎曲斷層模型的反演結果都一致地揭示了這次地震的主要特征.因此,不妨用本研究得到的震源模型計算這次地震的烈度.

        圖10展示了利用有限動態(tài)源模型計算的烈度分布.如果不考慮震源深度和斷層傾角,如圖10a所示,極震區(qū)烈度略大,極震區(qū)面積也略大,最大烈度達到近IX;如果考慮震源深度和斷層傾角,如圖10b所示,較深的點源的影響明顯減弱,極震區(qū)明顯縮小,最大烈度也相應減小.

        從圖10b所示的預估烈度可以看出,在震中東北和東南方向烈度衰減較慢, 而在西南和西北方向則較快,這種特征可以被近東西向斷層的東向破裂和近南北向斷層的南向破裂解釋.然而,如圖10c所示的調(diào)查烈度則顯示,震中西南和西北烈度衰減較慢,而震中東南和東北則較快,似乎恰好與預估烈度相反.另外, 預估烈度橢圓的長軸在南東—北西方向,而調(diào)查烈度橢圓的長軸在南南東—北北西方向,二者相差約15°左右.我們注意到,實際調(diào)查烈度似乎與早期定位的“線型”余震分布符合得更好,而最新的余震定位結果似乎更支持本研究預估的烈度分布.更為重要的是,來自實際加速度記錄的強地面運動特征可為本研究確定的烈度分布提供更有利的支持(Hu et al., 2016).

        對上述兩個實際震例的應用表明,預估烈度特征能夠反映有限動態(tài)源的震源特性,但與實際調(diào)查烈度相比仍有不可忽視的差別.我們認為造成這種差別的原因是兩方面的,一方面,預估方法本身沒有考慮場地效應;另一方面,實際調(diào)查烈度的信息不充分或不獨立或不客觀.

        5 討論與結論

        烈度是對地震復雜地面運動的簡單描述,更是對震害的直接反映(Wald et al., 1999a, 1999b).當一次破壞性地震發(fā)生后,地震造成的震害最受關注.在第一時間獲得烈度的空間特征相當于間接地獲得了震害的空間分布.因此,快速獲取烈度信息對于應急救援十分重要.

        烈度的概念誕生于現(xiàn)代地震儀器出現(xiàn)之前,烈度的經(jīng)驗關系似乎缺乏直接的物理意義,然而它能夠反映地震的強弱或震害的大小,因此非常實用且十分有效.多少年來,地震學家一直在致力于地震烈度與地震大小以及震中距或震源距關系的研究,但由于數(shù)字地震觀測和數(shù)字地震學的發(fā)展時間十分有限,所以上述經(jīng)驗關系大多基于點源假設(汪素云等,2000; Bormann, 2012; 俞言祥等, 2013).因此,探究基于近乎符合實際的有限動態(tài)源的烈度快速預估方法十分必要.

        經(jīng)典的烈度衰減關系基于地震點源,而實際的地震震源不但具有一定尺度且往往具有一定的幾何學屬性和運動學屬性.然而,實際地震的震源總可以看成是眾多點源的集合,這為我們基于點源的經(jīng)驗關系預估有限動態(tài)源的烈度分布提供了可能.

        本研究提出的預估方法通過場點的最終烈度取決于該點的最大烈度的假設考慮了震源的有限性,通過引入經(jīng)驗性常數(shù)μ考慮了由于破裂傳播的多普勒效應引起的烈度放大作用,通過引入經(jīng)驗性常數(shù)n考慮相對于破裂傳播方向的不同方位的烈度衰減特性,還通過引入破裂速度考慮了破裂速度相對于S波速度的大小對烈度的影響,因此本方法充分考慮了有限動態(tài)源對烈度的幾乎所有影響.這種影響通過數(shù)值試驗可以得到清楚的闡釋.

        本研究將提出的方法應用于近年來發(fā)生在我國的兩次中強地震,2013年4月20日的四川蘆山MS7.0地震和2014年8月3日的云南魯?shù)镸S6.5地震.首先,關于這兩次地震震源過程模型有較多研究,而且我們?nèi)栽陉P注這兩次地震的震源動態(tài)模型的研究;其次,這兩次地震的震源模型的有限性以及幾何學特性和運動學特性具有代表性,尤其是魯?shù)镸S6.5地震的共軛斷層或者彎曲斷層,便于體現(xiàn)這些屬性對烈度的影響;最后,這兩個地震都具有完整的調(diào)查烈度,便于對比討論.

        我們國家還沒有廣泛使用的烈度隨震源距的衰減關系,但事實上深度不同的震源在同一場點的烈度必然不同.為了討論問題的方便,我們在預估烈度時不但計算了烈度隨震中距的變化,也人為地將經(jīng)驗關系中的震中距當作震源距計算了烈度隨震源距的變化.很顯然,只考慮震中距必然會放大較深震源對烈度的作用.由此看來,一個烈度與震源距的經(jīng)驗關系更恰當.不過,考慮到把震中距人為地改為震源距不會影響震源的有限性以及幾何學和運動學特性對烈度空間分布的影響,所以本研究姑且把基于震源距的預估烈度和實際調(diào)查烈度進行了對比分析.

        對比蘆山MS7.0地震的預估烈度和調(diào)查烈度,我們注意到,預估烈度能夠充分反映斷層的有限性以及幾何學和運動學特性對烈度的作用,但不同于實際的調(diào)查烈度.對比魯?shù)镸S6.5地震的預估烈度和調(diào)查烈度,同樣發(fā)現(xiàn),預估烈度能夠充分反映斷層的有限性以及幾何學和運動學特性對烈度的作用,但仍有別于實際的調(diào)查烈度.我們認為造成這種差異的原因是兩方面的,一方面,預估方法本身沒有考慮場地效應;另一方面,實際調(diào)查烈度的信息不獨立和/或不充分和/或不客觀.

        綜上所述,本研究提出的烈度預估方法是對已有烈度衰減關系的拓展,繼承了經(jīng)典關系的實用和高效的優(yōu)點,也充分考慮了實際地震震源的有限性以及其幾何學屬性和運動學屬性對烈度的影響.這一點不但從本文的數(shù)值試驗可以看出,而且從對實際震例的應用也可得以體會.

        需要特別說明的是,預估結果的準確性或可靠性歸根到底依賴于已有的實際調(diào)查烈度.作為本方法核心的基于點源的經(jīng)驗衰減關系歸根到底來源于以前的實際調(diào)查烈度,實際調(diào)查烈度的錯誤必然導致經(jīng)驗關系的錯誤;本方法中經(jīng)驗性常數(shù)μ和n歸根到底也依賴于已有的實際調(diào)查烈度,從錯誤的實際調(diào)查烈度無法得到正確的經(jīng)驗性常數(shù).

        最后需要強調(diào)的是,本文提出的方法只考慮了震源效應,沒能考慮場地的作用.

        致謝俞言祥研究員為烈度計算提供了指導;中國地震局地球物理研究所國家測震臺網(wǎng)數(shù)據(jù)備份中心(doi:10.11998/SeisDmc/SN)、云南、四川、重慶、青海、貴州、甘肅、陜西、西藏、湖北、寧夏、湖南、廣西、河南、山西、內(nèi)蒙古地震臺網(wǎng)為本研究提供地震波形數(shù)據(jù);地震科學探測臺陣項目為巧家臺陣提供了6臺甚寬帶地震儀.

        Atkinson G M, Kaka S L. 2007. Relationships between Felt intensity and instrumental ground motion in the central United States and California.Bull.Seism.Soc.Amer., 97(2): 497-510. Bormann P. 2012. New Manual of Seismological Observatory Practice (NMSOP-2). IASPEI, GFZ German Research Centre for Geosciences, Potsdam.

        Chen K, Yu Y X, Gao M T, et al. 2014. Research on method of estimating ground motion parameters in western China using intensity data.SeismologyandGeology(in Chinese), 36(4): 1043-1052.

        Chen K, Yu Y X, Gao M T, et al. 2015. ShakeMap of peak ground acceleration for 2014 Ludian, Yunnan,MS6.5 earthquake.ActaSeismologicaSinica(in Chinese), 37(3): 429-436.

        Fang L H, Wu J P, Wang W L, et al. 2014. Relocation of the aftershock sequence of theMS6.5 Ludian earthquake and its seismogenic structure.SeismologyandGeology(in Chinese), 36(4): 1173-1185.

        Gutenberg B, Richter C F. 1942. Earthquake magnitude, intensity, energy, and acceleration.Bull.Seism.Soc.Amer., 32(3): 163-191. Gutenberg B, Richter C F. 1956. Earthquake magnitude, intensity, energy, and acceleration (second paper).Bull.Seism.Soc.Amer., 46(2): 105-145.

        Hao J L, Ji C, Wang W M, et al. 2013. Rupture history of the 2013MW6.6 Lushan earthquake constrained with local strong motion and teleseismic body and surface waves.Geophys.Res.Lett., 40(20): 5371-5376.

        Hu J J, Zhang Q, Jiang Z J, et al. 2016. Characteristics of strong ground motions in the 2014MS6.5 Ludian earthquake, Yunnan, China.J.Seismol., 20(1): 361-373, doi: 10.1007/s10950-015-9532-x.

        Karim K R, Yamazaki F. 2002. Correlation of JMA instrumental seismic intensity with strong motion parameters.EarthquakeEngineering&StructuralDynamics, 31(5): 1191-1212.

        Lay T, Wallace T C. 1995. Modern Global Seismology. San Diego: Academic Press.

        Li D H, Zuo H Q. 1991. Conversion between seismic intensity and peak ground shaking.ActaSeismologicaSinica(in Chinese), 13(1): 32-40.

        Li S Y, Jin X, Chen X, et al. 2002. Rapid reporting of peak strong motion and seismic intensity.EarthquakeEngineeringandEngineeringVibration(in Chinese), 22(6): 1-7.

        Liu C L, Zheng Y, Ge C, et al. 2013. Rupture process of theM7.0 Lushan earthquake.ScienceChina:EarthSciences(in Chinese), 43(6): 1020-1026. Liu C L, Zheng Y, Xiong X, et al. 2014. Rupture process ofMS6.5 Ludian earthquake constrained by regional broadband seismograms.ChineseJ.Geophys. (in Chinese), 57(9): 3028-3037, doi: 10.6038/cjg20140927.

        Meng L Y, Zhou L Q, Liu J. 2014. Estimation of the near-fault strong ground motion and intensity distribution of the 2013 Lushan, Sichuan,MS7.0 earthquake.ChineseJ.Geophys. (in Chinese), 57(2): 441-448, doi: 10.6038/cjg20140210.

        Moratto L, Gosta G, Suhadolc P. 2009. Real-time generation of ShakeMaps in the southeastern Alps.Bull.Seism.Soc.Amer., 99(4): 2489-2501. Murphy J R, O′Brien L J. 1977. The correlation of peak ground acceleration amplitude with seismic intensity and other physical parameters.Bull.Seism.Soc.Amer., 67(3): 877-915. Shearer P M. 1999. Introduction to Seismology. Cambridge: Cambridge University Press. Teng T L, Wu L D, Shin T C, et al. 1997. One minute after: Strong-motion map, effective epicenter and effective magnitude.Bull.Seism.Soc.Amer., 87(5): 1209-1219.

        Wald D J, Quitoriano V, Heaton T H, et al. 1999a. Relationships between peak ground acceleration, peak ground velocity, and modified Mercalli intensity in California.EarthquakeSpectra, 15(3): 557-564.

        Wald D J, Quitoriano V, Dengler L A, et al. 1999b. Utilization of the Internet for rapid community intensity maps.Seism.Res.Lett., 70(6): 680-597.

        Wang S Y, Yu Y X, Gao A J, et al. 2000. Development of Attenuation Relations for Ground Motion in China.EarthquakeResearchinChina(in Chinese), 16(2): 99-106.

        Wang W L, Wu J P, Fang L H, et al. 2014. Double difference location of the LudianMS6.5 earthquake sequences in Yunnan province in 2014.ChineseJ.Geophys. (in Chinese), 57(9): 3042-3051, doi: 10.6038/cjg20140929.

        Wang W M, Hao J L, Yao Z X. 2013. Preliminary result for rupture process of Apr. 20, 2013, Lushan Earthquake, Sichuan, China.ChineseJ.Geophys. (in Chinese), 56(4): 1412-1417, doi: 10.6038/cjg20130436.

        Wang Y S, Zhou Z H, Wang W. 2008. A hypothesis testing based method for seismic intensity rapid assessment by strong ground motion parameters.JournalofEarthquakeEngineeringandEngineeringVibration(in Chinese), 28(5): 49-54.

        Wu Y M, Shin T C, Chang C H. 2001. Near real-time mapping of peak ground acceleration and peak ground velocity following a strong earthquake.Bull.Seism.Soc.Amer., 91(5): 1218-1228. Wu Y M, Teng T L, Shin T C, et al. 2003. Relationship between peak ground acceleration, peak ground velocity, and intensity in Taiwan.Bull.Seism.Soc.Amer., 93(1): 386-396.

        Xu F K, Li J, Su Y J. 2014. Relocations of Yunnan LudianMS6.5 earthquake sequences in 2014.JournalofSeismologicalResearch(in Chinese), 37(4): 515-522. Xu J X, Zhang Z G, Dai W J, et al. 2015. Preliminary simulation of seismic wave propagation and the intensity map for the 25 April 2015 Nepal earthquake.ChineseJ.Geophys. (in Chinese), 58(5): 1812-1817, doi: 10.6038/cjg20150531.

        Xu L S, Du H L, Yan C, et al. 2013a. A method for determination of earthquake hypocentroid: time-reversal imaging technique (I)—Principle and numerical tests.ChineseJ.Geophys. (in Chinese), 56(4): 1190-1206, doi: 10.6038/cjg20130414. Xu L S, Yan C, Zhang X, et al. 2013b. Where did the LushanMS7.0 earthquake occur in the world?ChineseJ.Geophys. (in Chinese), 56(9): 2982-2993, doi: 10.6038/cjg20130912.

        Xu L S, Yan C, Zhang X, et al. 2013c. A method for determination of earthquake hypocentroid: Time-reversal imaging technique (II)—An examination based on people-made earthquakes.ChineseJ.Geophys. (in Chinese), 56(12): 4009-4027, doi: 10.6038/cjg20131207.

        Xu L S, Zhang X, Yan C, et al. 2014. Analysis of the Love waves for the source complexity of the LudianMS6.5 earthquake.ChineseJ.Geophys. (in Chinese), 57(9): 3006-3017, doi: 10.6038/cjg20140925.

        Xu X W, Jiang G Y, Yu G H, et al. 2014. Discussion on seismogenic fault of the LudianMS6.5 earthquake and its tectonic attribution.ChineseJ.Geophys. (in Chinese), 57(9): 3060-3068, doi: 10.6038/cjg20140931.

        Yu Y X, Li S Y, Xiao L. 2013. Development of ground motion attenuation relations for the new seismic hazard map of China.TechnologyforEarthquakeDisasterPrevention(in Chinese), 8(1): 24-33.

        Zhang G W, Lei J S, Liang S S, et al. 2014. Relocations and focal mechanism solutions of the 3 August 2014 Ludian, YunnanMS6.5 earthquake sequence.ChineseJ.Geophys. (in Chinese), 57(9): 3018-3027, doi: 10.6038/cjg20140926.

        Zhang Y, Xu L S, Chen Y T. 2013. Rupture process of the Lushan 4.20 earthquake and preliminary analysis on the disaster-causing mechanism.ChineseJ.Geophys. (in Chinese), 56(4): 1408-1411, doi: 10.6038/cjg20130435.

        Zhang Y, Wang R J, Chen Y T, et al. 2014a. Kinematic rupture model and hypocenter relocation of the 2013MW6.6 Lushan earthquake constrained by strong-motion and teleseismic data.Seism.Res.Lett., 85(1): 15-22.

        Zhang Y, Wang R J, Zschau J, et al. 2014b. Automatic imaging of earthquake rupture processes by iterative deconvolution and stacking of high-rate GPS and strong motion seismograms.J.Geophys.Res.:SolidEarth, 119(7): 5633-5650, doi: 10.1002/2013JB010469.

        Zhang Y, Xu L S, Chen Y T, et al. 2014. Rupture process of the 3 August 2014 Ludian, Yunnan,MW6.1 (MS6.5) earthquake.ChineseJ.Geophys. (in Chinese), 57(9): 3052-3059, doi: 10.6038/cjg20140930.

        Zhang Y, Chen Y T, Xu L S, et al. 2015. The 2014MW6.1 Ludian, Yunnan, earthquake: A complex conjugated ruptured earthquake.ChineseJ.Geophys. (in Chinese), 58(1): 153-162, doi: 10.6038/cjg20150113.

        Zhao C P, Zhou L Q, Chen Z L. 2013. Source rupture process of LushanMS7.0 earthquake, Sichuan, China and its tectonic implications.Chin.Sci.Bull., 58(28-29): 3444-3450.

        Zheng X F, Ouyang B, Zhang D N, et al. 2009. Technical system construction of Data Backup Centre for China Seismograph Network and the data support to researches on the Wenchuan earthquake.ChineseJ.Geophys. (in Chinese), 52(5): 1412-1417, doi: 10.3969/j.issn.0001-5733.2009.05.031.

        附中文參考文獻

        陳鯤, 俞言祥, 高孟潭等. 2014. 中國西部地區(qū)利用烈度數(shù)據(jù)估計地震動參數(shù)的方法. 地震地質(zhì), 36(4): 1043-1052.

        陳鯤, 俞言祥, 高孟潭等. 2015. 2014年云南魯?shù)镸S6.5地震峰值加速度震動圖. 地震學報, 37(3): 429-436.

        房立華, 吳建平, 王未來等. 2014. 云南魯?shù)镸S6.5地震余震重定位及其發(fā)震構造. 地震地質(zhì), 36(4): 1173-1185.

        李大華, 左慧強. 1991. 地震烈度與地震動峰值的轉(zhuǎn)換. 地震學報, 13(1): 32-40.

        李山有, 金星, 陳先等. 2002. 地震動強度與地震烈度速報研究. 地震工程與工程振動, 22(6): 1-7.

        劉成利, 鄭勇, 葛粲等. 2013. 2013年蘆山7.0級地震的動態(tài)破裂過程. 中國科學: 地球科學, 43(6): 1020-1026. 劉成利, 鄭勇, 熊熊等. 2014. 利用區(qū)域?qū)掝l帶數(shù)據(jù)反演魯?shù)镸S6.5級地震震源破裂過程. 地球物理學報, 57(9): 3028-3037, doi: 10.6038/cjg20140927.

        孟令媛, 周龍泉, 劉杰. 2014. 2013年四川蘆山MS7.0地震近斷層強地面運動模擬及烈度分布估計. 地球物理學報, 57(2): 441-448, doi: 10.6038/cjg20140210.

        汪素云, 俞言祥, 高阿甲等. 2000. 中國分區(qū)地震動衰減關系的確定. 中國地震, 16(2): 99-106.

        王未來, 吳建平, 房立華等. 2014. 2014年云南魯?shù)镸S6.5地震序列的雙差定位. 地球物理學報, 57(9): 3042-3051, doi: 10.6038/cjg20140929.

        王衛(wèi)民, 郝金來, 姚振興. 2013. 2013年4月20日四川蘆山地震震源破裂過程反演初步結果. 地球物理學報, 56(4): 1412-1417, doi: 10.6038/cjg20130436.

        王玉石, 周正華, 王偉. 2008. 基于假設檢驗的地震動強度(烈度)速報方法. 地震工程與工程振動, 28(5): 49-54.

        徐甫坤, 李靜, 蘇有錦. 2014. 2014年云南魯?shù)?.5級地震序列重定位研究. 地震研究, 37(4): 515-522.

        徐劍俠, 張振國, 戴文杰等. 2015. 2015年4月25日尼泊爾地震波場傳播及烈度初步模擬分析. 地球物理學報, 58(5): 1812-1817, doi: 10.6038/cjg20150531.

        許力生, 杜海林, 嚴川等. 2013a. 一種確定震源中心的方法: 逆時成像技術(一)——原理與數(shù)值實驗. 地球物理學報, 56(4): 1190-1206, doi: 10.6038/cjg20130414.

        許力生, 嚴川, 張旭等. 2013b. 蘆山MS7.0地震究竟發(fā)生在哪里?地球物理學報, 56(9): 2982-2993, doi: 10.6038/cjg20130912.

        許力生, 嚴川, 張旭等. 2013c. 一種確定震源中心的方法: 逆時成像技術(二)——基于人工地震的檢驗. 地球物理學報, 56(12): 4009-4027, doi: 10.6038/cjg20131207.

        許力生, 張旭, 嚴川等. 2014. 基于勒夫波的魯?shù)镸S6.5地震震源復雜性分析. 地球物理學報, 57(9): 3006-3017, doi: 10.6038/cjg20140925.

        徐錫偉, 江國焰, 于貴華等. 2014. 魯?shù)?.5級地震發(fā)震斷層判定及其構造屬性討論. 地球物理學報, 57(9): 3060-3068, doi: 10.6038/cjg20140931.

        俞言祥, 李山有, 肖亮. 2013. 為新區(qū)劃圖編制所建立的地震動衰減關系. 震災防御技術, 8(1): 24-33. 張廣偉, 雷建設, 梁姍姍等. 2014. 2014年8月3日云南魯?shù)镸S6.5級地震序列重定位與震源機制研究. 地球物理學報, 57(9): 3018-3027, doi: 10.6038/cjg20140926.

        張勇, 許力生, 陳運泰. 2013. 蘆山4.20地震破裂過程及其致災特征初步分析. 地球物理學報, 56(4): 1408-1411, doi: 10.6038/cjg20130435.

        張勇, 許力生, 陳運泰等. 2014. 2014年8月3日云南魯?shù)镸W6.1(MS6.5)地震破裂過程. 地球物理學報, 57(9): 3052-3059, doi: 10.6038/cjg20140930.

        張勇, 陳運泰, 許力生等. 2015. 2014年云南魯?shù)镸W6.1地震: 一次共軛破裂地震. 地球物理學報, 58(1): 153-162, doi: 10.6038/cjg20150113.

        趙翠萍, 周連慶, 陳章立. 2013. 2013年四川蘆山MS7.0 級地震震源破裂過程及其構造意義. 科學通報, 58(20): 1894-1900.

        鄭秀芬, 歐陽飚, 張東寧等. 2009. “國家數(shù)字測震臺網(wǎng)數(shù)據(jù)備份中心”技術系統(tǒng)建設及其對汶川大地震研究的數(shù)據(jù)支撐. 地球物理學報, 52(5): 1412-1417, doi: 10.3969/j.issn.0001-5733.2009.05.031.

        (本文編輯何燕)

        A method for estimating the earthquake intensity caused by a finite-dynamic source

        XU Li-Sheng, ZHANG Xu, WEI Qiang, LI Chun-Lai

        InstituteofGeophysics,ChinaEarthquakeAdministration,Beijing100081,China

        At first, a method to be used to quickly estimate the intensity caused by a destructive earthquake is proposed based on empirical relationship between intensity and magnitude as well as epicentral or hypocentral distance, considering the finiteness and the kinematic and geometrical feature of real faults. Then, the numerical tests are employed to illustrate the effects of the parameters reflecting the kinematic and geometrical feature and the influence of fault dip and rupture velocity. Finally, the method is applied to two destructive earthquakes which have occurred in our country recent years in order to illustrate its practical value and shortage. At the present time when real-time obtaining of the rupture process image is coming true, the method by which the intensity is able to be quickly obtained appears especially significant.

        Earthquake intensity; Finite-dynamic source; Estimating method; Numerical test; Practical application

        10.6038/cjg20161015.

        地震行業(yè)專項(201408014)和中國地震局地球物理研究所基本業(yè)務費(DQJB14B01)聯(lián)合資助.

        許力生,男,研究員,主要從事地震學研究.E-mail:xuls@cea-igp.ac.cn

        10.6038/cjg20161015

        P315

        2016-02-01,2016-03-22收修定稿

        許力生, 張旭, 魏強等. 2016. 一種基于有限動態(tài)源的烈度估計方法. 地球物理學報,59(10):3684-3695,

        Xu L S, Zhang X, Wei Q, et al. 2016. A method for estimating the earthquake intensity caused by a finite-dynamic source.ChineseJ.Geophys. (in Chinese),59(10):3684-3695,doi:10.6038/cjg20161015.

        猜你喜歡
        模型
        一半模型
        一種去中心化的域名服務本地化模型
        適用于BDS-3 PPP的隨機模型
        提煉模型 突破難點
        函數(shù)模型及應用
        p150Glued在帕金森病模型中的表達及分布
        函數(shù)模型及應用
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權M-估計的漸近分布
        3D打印中的模型分割與打包
        69天堂国产在线精品观看| 精品伊人久久大香线蕉综合| 中文字幕无码日韩专区免费| 婷婷亚洲国产成人精品性色 | 久久久久人妻精品一区二区三区| 日韩毛片免费无码无毒视频观看| 欧美另类在线视频| 成年女人18毛片毛片免费| 国产一区二区黄色网页| 少妇人妻精品一区二区三区| 中文在线天堂网www| 国产精品自在在线午夜出白浆| 美艳善良的丝袜高跟美腿| 99久久精品费精品国产一区二| 国产精品黄网站免费观看| 国产一区二区三区资源在线观看 | 男女做羞羞事的视频网站| 日本熟妇hdsex视频| 亚洲区在线| 国产91熟女高潮一曲区| 日韩在线观看入口一二三四| 女人让男人桶爽30分钟| 成人欧美在线视频| 一区二区三区在线观看视频| 久久国语露脸国产精品电影| 好男人日本社区www| 亚洲一区二区三区在线观看播放| 91久久精品一区二区三区大全| 国产三级在线观看完整版| 国产欧美日韩a片免费软件| 日本精品久久久久中文字幕1| 亚洲国产精品婷婷久久| 欧美人与动牲交a精品| 亚洲熟妇网| 亚洲一区二区三区免费的视频| 人人做人人爽人人爱| 亚洲综合色丁香婷婷六月图片| 国产一区二区三区亚洲精品| 亚洲中文字幕精品乱码2021| 亚洲一线二线三线写真| 国内视频一区|