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

    貝葉斯同化重力反演方法構(gòu)建龍門(mén)山地殼密度模型

    2021-04-13 17:44:39李紅蕾陳石莊建倉(cāng)張貝石磊
    地球物理學(xué)報(bào) 2021年4期
    關(guān)鍵詞:參考模型重力反演

    李紅蕾, 陳石*, 莊建倉(cāng), 張貝, 石磊

    1 中國(guó)地震局地球物理研究所, 北京 100081 2 北京白家疃地球科學(xué)國(guó)家野外科學(xué)觀測(cè)研究站, 北京 100095 3 日本數(shù)理統(tǒng)計(jì)研究所, 東京 190-8562

    0 引言

    地球物理反演技術(shù)是研究地殼內(nèi)部結(jié)構(gòu)和性質(zhì)的重要工具.由于地球物理反演問(wèn)題先天具有多解性,通常反演結(jié)果在數(shù)學(xué)上適定,但可能與實(shí)際地質(zhì)情況并不相符(Li and Oldenburg, 1998; Williams et al., 2004; Howe, 2009; Dirian et al.,2016).現(xiàn)今在很多研究熱點(diǎn)地區(qū),利用多種手段獲得的地殼結(jié)構(gòu)模型和認(rèn)識(shí)越來(lái)越多.如何充分融合已有的模型,發(fā)揮各種地球物理手段的優(yōu)勢(shì),減小反演結(jié)果的不確定性,無(wú)疑是地球物理聯(lián)合反演研究要解決的核心問(wèn)題(Welford et al., 2018; Yang et al., 2012).

    如果將所有已有模型解釋結(jié)果看作先驗(yàn)認(rèn)識(shí),當(dāng)今地學(xué)家面臨的首要挑戰(zhàn)是如何從已有模型數(shù)據(jù)中提取盡可能多的有用信息以獲取新的見(jiàn)解.與常規(guī)技術(shù)相比,數(shù)據(jù)同化提供了一套新模式,用于發(fā)現(xiàn)常規(guī)技術(shù)不容易揭示的數(shù)據(jù)和模型結(jié)構(gòu)之間的關(guān)系(Bergen et al., 2019).數(shù)據(jù)同化將觀測(cè)數(shù)據(jù)、計(jì)算模型和先驗(yàn)推斷集成到一個(gè)系統(tǒng)中,通過(guò)對(duì)系統(tǒng)中不確定性的估計(jì)和控制,來(lái)提高后驗(yàn)估計(jì)模型的精度.在當(dāng)今這個(gè)大數(shù)據(jù)時(shí)代中,數(shù)據(jù)同化方法已成功地應(yīng)用到了多個(gè)科學(xué)領(lǐng)域的建模、數(shù)據(jù)分析和預(yù)測(cè)中(Riggelsen et al.,2007; de Wit et al.,2013).

    重力反演是探測(cè)地球深部構(gòu)造的有效手段之一,具有對(duì)地下物質(zhì)密度變化敏感、水平分辨能力強(qiáng)、成本低等優(yōu)勢(shì)(Kearey et al., 2002; Hinze et al., 2013).然而,由于觀測(cè)數(shù)據(jù)誤差和核函數(shù)算子本身的性質(zhì),重力反演具有多解性(Green, 1975; Fedi and Rapolla, 1999).除了通過(guò)增加觀測(cè)數(shù)據(jù)和減小觀測(cè)誤差,在一定程度上降低反演的多解性.之外,在反演中加入先驗(yàn)信息進(jìn)行約束是更有效的方法(Li and Oldenburg,1998;陳石等,2010;李紅蕾等,2016).近年來(lái),重力約束條件越來(lái)越受到重視,許多地球物理學(xué)者都提出了各種各樣的約束條件以及相應(yīng)的反演方法.約束總體上可以分為兩類(lèi),即數(shù)學(xué)結(jié)構(gòu)約束和參考模型約束.前者包括最小構(gòu)造約束(Last and Kubik, 1983)、深度加權(quán)約束(Li and Oldenburg, 1998)、平坦度和光滑度約束(Boulanger and Chouteau, 2001)等,后者主要包括地震波速轉(zhuǎn)化密度約束(王新勝等,2013)、地質(zhì)約束(Nabighian et al., 2005)、巖性約束(Geng et al., 2019)等.在傳統(tǒng)的三維重力反演算法中,上述約束和參考都是通過(guò)阻尼因子被引入到反演方程中,并通過(guò)廣義交叉驗(yàn)證(GCV)(Golub et al.,1979)或L曲線準(zhǔn)則(Hansen,2001)獲取的單個(gè)阻尼因子實(shí)現(xiàn)對(duì)觀測(cè)數(shù)據(jù)和先驗(yàn)信息權(quán)重的平衡.然而,當(dāng)在先驗(yàn)假設(shè)(光滑先驗(yàn)、深度加權(quán))和參考權(quán)重設(shè)置增多的情況下,傳統(tǒng)算法很難依據(jù)數(shù)據(jù)特征實(shí)現(xiàn)對(duì)多個(gè)約束信息權(quán)重參數(shù)的優(yōu)選,這很容易造成反演結(jié)果偏離實(shí)際(Li and Oldenburg, 1998, 2000; Williams, 2008).

    Akaike(1977)最早利用貝葉斯方法來(lái)解決氣象領(lǐng)域的數(shù)據(jù)同化問(wèn)題.這種方法已成為數(shù)據(jù)同化的重要手段,其是以數(shù)據(jù)驅(qū)動(dòng)模式量化反演參數(shù),通過(guò)最小化Akaike貝葉斯信息量(ABIC)實(shí)現(xiàn)對(duì)觀測(cè)數(shù)據(jù)及先驗(yàn)約束不確定性的估計(jì)與控制,從而顯著提高模型估計(jì)精度,減少模型不確定性.同樣方法隨后成功應(yīng)用到時(shí)-空傳染型余震序列模型(ETAS模型)的參數(shù)估計(jì)(Ogata and Katsura, 1993)、三維地震層析成像、應(yīng)力圖像反演(Terakawa and Matsu′ura, 2008)、GPS斷層數(shù)據(jù)反演(Fukahata et al., 2004; Fukahata and Wright, 2008)、重力平差參數(shù)優(yōu)化(Chen at al., 2019)等多個(gè)地球物理學(xué)領(lǐng)域.

    本文探索通過(guò)數(shù)據(jù)同化手段,將貝葉斯同化策略引入到傳統(tǒng)的重力約束反演中,設(shè)計(jì)一種可以融合多種先驗(yàn)推斷的重力異常貝葉斯同化反演算法,以期得到更加準(zhǔn)確的最大后驗(yàn)概率意義下的模型參數(shù)估計(jì).文中第2部分給出了三維重力貝葉斯同化反演的求解方程和ABIC目標(biāo)函數(shù)的構(gòu)建;第3部分設(shè)計(jì)了兩個(gè)綜合實(shí)驗(yàn)?zāi)P?,測(cè)試了該方法在解決多種已知先驗(yàn)參考約束和深度加權(quán)約束中超參數(shù)和后驗(yàn)?zāi)P凸烙?jì)問(wèn)題的優(yōu)化效果,驗(yàn)證了該方法的可行性和有效性;第4部分和第5部分是應(yīng)用該方法和龍門(mén)山及周邊地區(qū)2′×2′的高精度WGM2012布格重力異常數(shù)據(jù)反演了該區(qū)的地殼密度模型, 并對(duì)龍門(mén)山地區(qū)地殼內(nèi)低密度分布、隆升變形機(jī)制、強(qiáng)震震源區(qū)環(huán)境等進(jìn)行了深入分析;最后,第6部分對(duì)本文研究方法和認(rèn)識(shí)進(jìn)行了總結(jié)和討論.

    1 方法原理

    一般空間域三維重力位場(chǎng)正演問(wèn)題,可以離散化并線性化為:

    d=Gm,

    (1)

    其中,G為N×M的二維核函數(shù)矩陣或格林函數(shù)矩陣,矩陣中每個(gè)元素值與模型中待計(jì)算密度單元和觀測(cè)點(diǎn)的位置相關(guān),m為M×1的一維密度單元矩陣,d為N×1的一維觀測(cè)矩陣.方程(1)中已知異常d求解m則稱(chēng)為重力反演問(wèn)題,一般觀測(cè)N?M,直接反演在數(shù)學(xué)上是沒(méi)有唯一解的.Tikhonov and Arsenin(1977)引入最小模約束后,求解方程(1)可以變?yōu)椋?/p>

    minφ=‖Gm-d‖2+λ‖m‖2,

    (2)

    其中λ是阻尼因子.

    公式(2)具有數(shù)學(xué)上的唯一解.當(dāng)有最小構(gòu)造、光滑度、深度加權(quán)、參考模型等多個(gè)約束存在時(shí),傳統(tǒng)反演方法通過(guò)式(3)將其加入到反演目標(biāo)函數(shù)中:

    (3)

    為解決上述問(wèn)題,我們提出了重力貝葉斯同化反演算法,其是在綜合考慮重力觀測(cè)數(shù)據(jù)誤差及模型約束不確定性的基礎(chǔ)上,通過(guò)引入ABIC準(zhǔn)則實(shí)現(xiàn)最大化邊際概率分布的數(shù)據(jù)及模型權(quán)重超參數(shù)自動(dòng)優(yōu)選的反演方法.算法的基本原理如下:給定特定模型下重力觀測(cè)數(shù)據(jù)的條件概率分布和模型的先驗(yàn)概率分布,根據(jù)貝葉斯公式,結(jié)合了觀測(cè)數(shù)據(jù),模型的后驗(yàn)概率密度函數(shù)(pdf)表示為:

    (4)

    其中,p(*|*)為條件概率密度函數(shù).如果數(shù)據(jù)和模型誤差都以正態(tài)分布,則重力觀測(cè)數(shù)據(jù)的條件概率為:

    (5)

    考慮最小模型和多個(gè)先驗(yàn)約束的模型先驗(yàn)概率分布可表示為:

    (6)

    其中,Cd為向量(d-Gm)的方差矩陣,Cm為模型約束算子方差矩陣.Cmr為(m-mrefi)量的方差矩陣.在λd、λms、λri已知的情況下,最大化后驗(yàn)分布概率求解m等價(jià)于最小化:

    (7)

    式中λ0、λms、λri為數(shù)據(jù)和模型權(quán)重系數(shù),又稱(chēng)為反演超參數(shù),其可以通過(guò)引入ABIC準(zhǔn)則來(lái)求解,即最小化統(tǒng)計(jì)量:

    (8)

    2 實(shí)驗(yàn)測(cè)試

    為測(cè)試上述方法的可行性和有效性,我們分別設(shè)計(jì)了兩組仿真測(cè)試模型,來(lái)檢驗(yàn)貝葉斯同化反演策略的實(shí)際效果.

    2.1 參考模型約束

    引入適當(dāng)?shù)膮⒖寄P图s束能在垂向分層上改善重力異常反演結(jié)果.通常認(rèn)為參考模型越準(zhǔn)確,反演結(jié)果越接近真實(shí)情況.但在實(shí)際問(wèn)題中,可參考的模型不止一個(gè),其誤差、不確定性差異都不盡相同,同時(shí)不同參考模型之間也往往不一致.對(duì)于多參考模型約束問(wèn)題,特別是在先驗(yàn)參考模型誤差未知情況下,如何分配每個(gè)參考模型的權(quán)重問(wèn)題一直是研究熱點(diǎn)(Hansen, 1994; Vogel, 1996).下面就本文提出的方法如何分配權(quán)重和改善反演結(jié)果進(jìn)行測(cè)試.

    這個(gè)數(shù)值實(shí)驗(yàn)中,真實(shí)模型在y方向沒(méi)有變化,在y=0處的垂直切片,模型起伏界面上方藍(lán)色部分密度為0.1 g·cm-3,下方紅色部分密度為0.2 g·cm-3,起伏界面上下密度差為0.1 g·cm-3(如圖1a).圖1b是仿真的重力異常數(shù)據(jù),在理論真實(shí)模型正演異常的基礎(chǔ)上,加入了5%的高斯誤差.

    在測(cè)試中,設(shè)計(jì)了兩個(gè)簡(jiǎn)單的參考模型,每個(gè)模型與真實(shí)模型都存在一定的誤差.其中,參考模型M1密度起伏界分界面形態(tài)與真實(shí)模型一致,但界面兩側(cè)的密度差比真實(shí)模型小10%,如圖1c所示;參考模型M2中的密度起伏界面位置與真實(shí)模型不一致,界面兩側(cè)的密度差真實(shí)模型大10%,如圖1d所示.此外,參考模型1和2中還分別加入了2%和1%高斯誤差.

    在模型測(cè)試過(guò)程中,我們選擇了四種不同的反演方案,方案1是僅有參考模型M1優(yōu)化約束反演;方案2是僅有參考約束模型M2優(yōu)化約束反演;方案3是參考模型M1和M2聯(lián)合非優(yōu)化約束反演;方案4是參考模型M1和M2聯(lián)合優(yōu)化約束反演.每種反演都實(shí)現(xiàn)了模型異常與數(shù)據(jù)異常的擬合,但反演的模型結(jié)果差異性顯著,詳見(jiàn)圖2所示.

    圖1 不同參考模型約束下反演算法測(cè)試(a) 真實(shí)模型; (b) 正演重力異常; (c) 參考模型M1; (d) 參考模型M2.Fig.1 Tests of inversion algorithm on different reference models with constraints(a) True model; (b) Forward modeling gravity anomaly; (c) Reference model M1; (d) Reference model M2.

    圖2 不同參考模型約束下反演結(jié)果(a) 參考模型M1約束反演結(jié)果; (b) 參考模型M2約束反演結(jié)果; (c) 參考模型M1和M2聯(lián)合非優(yōu)化約束反演結(jié)果; (d) 參考模型M1和M2聯(lián)合優(yōu)化約束反演結(jié)果,白色實(shí)線為真實(shí)模型密度分界面.Fig.2 Inversion results of different reference models with constraints(a) Reference model M1; (b) Reference model M2; (c) Reference models M1 and M2 combined with non-optimization constraint; (d) Reference models M1 and M2 combined with optimization constraint. The white solid line is the real model density interface.

    從圖2中的結(jié)果可以看出,圖2a的反演模型相比圖2b結(jié)果界面更清晰,但上下密度差偏差較大,圖2b的界面凸起部分效果比圖2a差.圖2c的聯(lián)合反演結(jié)果,相比前兩個(gè)無(wú)論在界面起伏和上下密度差方面都有一定程度改善,而與圖2d的貝葉斯同化聯(lián)合反演結(jié)果相比,在上下界面位置和密度差方面圖2d明顯優(yōu)于圖2c.

    在該模型測(cè)試中,參考模型1和2的歸一化權(quán)重參數(shù)分別為λ1和λ2.圖2a—d的反演結(jié)果,對(duì)應(yīng)的參數(shù)詳見(jiàn)表1.從表1中的結(jié)果可以看出,在圖2d的最優(yōu)化過(guò)程中,參考模型1的權(quán)重更大,對(duì)應(yīng)的ABIC值也最小.

    表1 不同參考約束下反演統(tǒng)計(jì)參數(shù)及ABIC參數(shù)特征Table 1 Inversion statistical parameters and ABIC parameter characteristics under different reference constraints

    結(jié)合表1中的反演統(tǒng)計(jì)參數(shù)結(jié)果來(lái)看,參考模型1約束反演優(yōu)化得到的超參數(shù)λ1為129.055,以此計(jì)算得到的ABIC值為-1355.069;參考模型2約束反演優(yōu)化得到的超參數(shù)λ2為160.062,計(jì)算得到的ABIC值為-1406.770;將上述兩個(gè)模型單獨(dú)優(yōu)化得到的超參數(shù)直接代入多模型約束,計(jì)算得到的ABIC值為-3217.861;將兩個(gè)模型同時(shí)優(yōu)化反演得到的最優(yōu)λ1和λ2分別為1798.5和1094.8,最小化ABIC值為-3700.195.

    將表1反演統(tǒng)計(jì)參數(shù)與圖2中的反演結(jié)果結(jié)合來(lái)看,與單個(gè)參考模型和非優(yōu)化權(quán)重的多個(gè)參考模型約束反演相比,ABIC最小化(ABIC最小值為-3700.195)給出的最優(yōu)超參數(shù)可有效降低反演的卡方值,提高反演的準(zhǔn)確性和精度.與手動(dòng)指定的模型約束權(quán)重進(jìn)行反演計(jì)算得到的反演結(jié)果(圖2c)相比,基于優(yōu)化ABIC的反演算法可依據(jù)觀測(cè)數(shù)據(jù)來(lái)選擇提取不同參考之間的有用信息,并實(shí)現(xiàn)不同參考模型與反演結(jié)果之間的數(shù)據(jù)同化.

    2.2 深度加權(quán)約束

    在式(1)的核函數(shù)矩陣G中,每個(gè)元素?cái)?shù)值大小都與觀測(cè)點(diǎn)和模型單元之間的距離平方成正比,在不同深度位置上,模型單元的核函數(shù)大小數(shù)值差異明顯,淺層單元的數(shù)值遠(yuǎn)大于深層.體現(xiàn)在反演結(jié)果中,常常會(huì)出現(xiàn)異常的變化僅集中在淺部模型單元上,通常稱(chēng)之為重力位場(chǎng)反演的趨膚效應(yīng)(Li and Oldenburg,1998).

    (9)

    其中,dz為垂向間隔,α和r的大小直接決定了近地表頂層壓制作用的大小,然而在實(shí)際應(yīng)用中α和r的值也主要根據(jù)經(jīng)驗(yàn)指定(Commer, 2011; 秦朋波和黃大年,2016).

    本節(jié)的測(cè)試模型如圖3所示,與圖2模型在y方向設(shè)置相同,即在y方向沒(méi)有變化,抽取模型在y=0處的垂直切片如圖3a所示.真實(shí)模型在深度方向由兩個(gè)規(guī)則長(zhǎng)方體組成,周邊藍(lán)色區(qū)域密度為零,每個(gè)長(zhǎng)方體的密度大小同為0.2 g·cm-3,長(zhǎng)方體的頂面埋深分別為2 km和11 km,厚度均為4 km.

    我們將在真實(shí)模型正演重力異常的基礎(chǔ)上添加了5%白噪聲的模擬數(shù)據(jù)作為觀測(cè)重力異常.由于重力反演很難分辨深度方向上的異常體,因此,需要給定一定的先驗(yàn)信息作為參考模型.在本次測(cè)試中,我們選擇了局部參考模型作為約束,即假設(shè)在x=0位置處的單元體密度為已知.一般在實(shí)際地殼結(jié)構(gòu)反演中,測(cè)井或接收函數(shù)方法可以提供此類(lèi)的局部模型信息.

    該模型測(cè)試中,我們給出了三種不同深度加權(quán)參數(shù)的反演結(jié)果,分別對(duì)應(yīng)了三種不同的深度加權(quán)參數(shù),如表2中α、β參數(shù)值所示,表2中不同的深度加權(quán)參數(shù)值分別對(duì)應(yīng)了圖3b、c、d的反演結(jié)果.當(dāng)表2中加權(quán)參數(shù)α值過(guò)大時(shí),得到的反演結(jié)果圖3b中高密度部分主要集中在底部,異常體的深度與真實(shí)模型埋深存在一定的差異;而當(dāng)加權(quán)參數(shù)β過(guò)大時(shí),得到的圖3c的反演結(jié)果其高密度主要集中在淺部,異常體埋深同樣存在較大差異.通過(guò)本文算法得到的最優(yōu)化深度加權(quán)參數(shù)反演結(jié)果如圖3d所示,異常體的深度信息較為準(zhǔn)確,得到的反演結(jié)果埋深和形態(tài)與真實(shí)模型基本一致.由此可得,與主觀指定的深度加權(quán)參數(shù)反演結(jié)果相比,本文提出的貝葉斯同化反演方法,在已知局部參考信息下,給出的深度加權(quán)更合理,反演的異常體深度與真實(shí)模型更加一致.

    圖3 深度加權(quán)參數(shù)優(yōu)化反演算法測(cè)試(a) 真實(shí)模型; (b) 過(guò)衰減加權(quán)系數(shù)反演結(jié)果; (c) 欠衰減加權(quán)系數(shù)反演結(jié)果; (d) 優(yōu)化加權(quán)系數(shù)反演結(jié)果.Fig.3 Tests of inversion algorithm with depth weighted parameter optimization(a) Real model; (b) Inversion result of over-attenuation weighted coefficient; (c) Inversion result of under-attenuation weighted coefficient; (d) Inversion results of optimized weighted coefficient.

    表2 不同深度超參數(shù)約束下反演統(tǒng)計(jì)及ABIC參數(shù)特征Table 2 Inversion statistics and ABIC parameter characteristics under different depth hyperparameter constraints

    綜合以上兩個(gè)模型的測(cè)試結(jié)果,可以認(rèn)為本文算法具備同化多個(gè)參考模型和優(yōu)化深度加權(quán)參數(shù)的能力.下面我們進(jìn)一步應(yīng)用該方法,選擇地球物理觀測(cè)資料豐富、地殼結(jié)構(gòu)參考模型較多的龍門(mén)山地區(qū)進(jìn)行實(shí)際數(shù)據(jù)測(cè)試.

    3 龍門(mén)山地殼模型

    3.1 龍門(mén)山構(gòu)造背景

    龍門(mén)山地區(qū)位于南北地震帶中南部位,擁有復(fù)雜的構(gòu)造邊界條件、活動(dòng)斷層系統(tǒng).其西部是活躍的青藏高原邊界,東部是穩(wěn)定的華南地臺(tái),是青藏高原主體向東擠出構(gòu)造邊界,地殼變形強(qiáng)烈,結(jié)構(gòu)復(fù)雜.地塊內(nèi)部褶皺斷裂廣泛發(fā)育,其中包括了多個(gè)大地構(gòu)造區(qū)邊界斷裂和控制強(qiáng)震活動(dòng)的活斷層:鮮水河斷裂帶、龍門(mén)山斷裂帶、東昆侖斷裂帶、龍泉山斷裂帶、龍日壩斷裂帶、毛爾蓋分支斷層等(徐錫偉等,2008).此外,沿著龍門(mén)山斷裂帶,還存在著約5 km的強(qiáng)烈高程梯度帶(如圖4所示).

    圖4 龍門(mén)山地區(qū)主要構(gòu)造特征與地震活動(dòng)紅色實(shí)線為斷裂,四條黑色實(shí)線為跨震中研究剖面,黃色圈為5級(jí)以上地震位置,白色實(shí)心圓為中國(guó)地震觀測(cè)網(wǎng)臺(tái)站位置,F(xiàn)1:鮮水河斷裂帶;F2:龍門(mén)山斷裂帶;F3:龍泉山斷裂帶;F4-1:龍日壩斷裂帶;F4-2:毛爾蓋分支斷層;F4:東昆侖斷裂帶(徐錫偉等,2008).Fig.4 Main tectonic features and seismic activity in and around the Longmenshan areaThe red solid lines represent faults, black solid lines are profiles through epicenters, yellow circles represent earthquakes of M5 or greater, and white solid circles are stations of China earthquake observation network. F1: Xianshuihe fault zone; F2: Longmenshan fault zone; F3: Longquanshan fault zone; F4-1: Longriba fault zone; F4-2: Maoergai branch fault; F4: East kunlun fault zone (Xu et al., 2008).

    同時(shí),該地區(qū)也屬于中國(guó)地震科學(xué)實(shí)驗(yàn)場(chǎng)區(qū),過(guò)去幾十年中以中國(guó)地震科學(xué)臺(tái)陣項(xiàng)目為代表的大規(guī)模地球物理觀測(cè)相繼開(kāi)展,已積累了大量的地球物理探測(cè)工作成果(如Yao et al., 2008; Li et al., 2011; Zhang et al., 2011; 王緒本等,2018).這些深部成果對(duì)該區(qū)的殼幔結(jié)構(gòu)及其相關(guān)動(dòng)力學(xué)問(wèn)題達(dá)成了部分共識(shí),但在一些基本問(wèn)題上依舊存在爭(zhēng)議,如在該區(qū)的地殼上地幔變形機(jī)制的解釋上,就有殼幔連續(xù)變形機(jī)制(England and Molnar, 1997)、塊體擠出機(jī)制 (Tapponnier et al., 1982,2001)、下地殼流機(jī)制(Clark and Royden, 2000; Royden et al., 1997)等多種模式.

    此外,龍門(mén)山地區(qū)地震多發(fā),如2008年汶川8.0級(jí)和2013年蘆山7.0級(jí)強(qiáng)震.雖然國(guó)內(nèi)外研究學(xué)者對(duì)該區(qū)大震震源區(qū)的深部孕震結(jié)構(gòu)等進(jìn)行了大量的研究,但對(duì)于地震孕育深部地殼結(jié)構(gòu)特征及相關(guān)動(dòng)力學(xué)機(jī)制還存在分歧.如:Zhang等(2011)通過(guò)地震層析結(jié)果,認(rèn)為龍門(mén)山斷裂帶兩側(cè)高低速異常的邊界是該區(qū)強(qiáng)震的凹凸區(qū);而房立華等(2013)和王夫運(yùn)等(2015)的地震測(cè)深剖面顯示斷裂帶下方中滑脫層是產(chǎn)生研究區(qū)地震的主要原因;張培震等(2008)、楊文采等(2015)和王緒本等(2018)的多種研究結(jié)果表明研究區(qū)強(qiáng)震的發(fā)生與地殼內(nèi)部存在的低速高導(dǎo)層有關(guān).

    密度作為地球內(nèi)部構(gòu)造最重要的參數(shù)之一,能夠很好地反應(yīng)地下物質(zhì)的運(yùn)動(dòng)和變化,高精度的地殼三維密度結(jié)構(gòu)對(duì)該區(qū)構(gòu)造形成及演化、地震孕育等深部動(dòng)力學(xué)過(guò)程的深入認(rèn)識(shí)具有重要作用.雖然前人已經(jīng)在該研究區(qū)內(nèi)進(jìn)行了一定的重力密度探測(cè)研究工作,這些探測(cè)成果為我們理解和認(rèn)識(shí)研究區(qū)地幔變形及強(qiáng)震風(fēng)險(xiǎn)源區(qū)的地殼結(jié)構(gòu)和物性特征研究提供了重要的深部資料(姜文亮和張景發(fā),2011;申重陽(yáng)等, 2015),但這些成果主要集中在二維.已有三維地殼密度數(shù)據(jù)分辨率和精度較低(方劍和許厚澤,1997; 楊文采等; 2015;李紅蕾等,2017;Li et al., 2017).不足以識(shí)別地殼和上地殼深度的細(xì)節(jié)特征,也不能為解決該區(qū)殼幔變形、地震孕育及相關(guān)動(dòng)力學(xué)過(guò)程爭(zhēng)議提供很好的論據(jù)(王椿鏞等, 2016).

    3.2 研究區(qū)已有地殼結(jié)構(gòu)

    眾所周知,重力數(shù)據(jù)水平分辨率高,但垂向分辨率低,在實(shí)際反演過(guò)程中,要想得到有地質(zhì)意義的結(jié)果,還需要添加深部參考約束.相對(duì)重力方法來(lái)講,地震成像方法具有較好的垂向分辨率,重力反演中將地震成像結(jié)果作為參考約束,可以集兩種數(shù)據(jù)之長(zhǎng),提高重力深部結(jié)構(gòu)的探測(cè)能力.本文選擇的研究區(qū),在地震學(xué)研究方面,已有體波成像(如Zhang et al., 2011; Wang et al., 2017)、面波成像(如Yao et al., 2008; Li et al., 2011)、接收函數(shù)反演(如Bao et al., 2015; Liu et al., 2014)、地震測(cè)深成像(如張先康等, 2007; 王夫運(yùn)等,2008)等多個(gè)結(jié)果.然而,由于不同學(xué)者使用的研究數(shù)據(jù)及方法不同,獲得的參考地震模型結(jié)果在數(shù)據(jù)分布、分辨率、誤差及物性等方面常存在較大差異,如表3所示.

    表3 研究區(qū)內(nèi)已有部分地震波速成果Table 3 Partial seismic wave velocities in the study area available

    如果將已有的模型解釋結(jié)果看作先驗(yàn)認(rèn)識(shí),那么如何根據(jù)這些先驗(yàn)去改進(jìn)反演模型的估計(jì)?這是本文提出的重力異常貝葉斯同化反演新算法應(yīng)該回答的問(wèn)題.我們將以多種不同類(lèi)型的地震速度轉(zhuǎn)換密度作為已知先驗(yàn)參考,利用高精度的衛(wèi)星重力場(chǎng)模型數(shù)據(jù),通過(guò)新算法實(shí)現(xiàn)對(duì)不同先驗(yàn)參考約束的取長(zhǎng)補(bǔ)短,剔除偏差數(shù)據(jù)在反演計(jì)算中的作用,以期得到同時(shí)符合計(jì)算系統(tǒng)不同先驗(yàn)假設(shè)的最優(yōu)高精度地殼結(jié)構(gòu)模型.同時(shí)測(cè)試算法在實(shí)際重力資料反演中的效果,為研究該區(qū)的地殼變形機(jī)制、強(qiáng)震孕育環(huán)境及相關(guān)動(dòng)力學(xué)提供有益參考.

    3.3 參考模型建立

    轉(zhuǎn)換為密度的兩個(gè)參考模型信息在30 km深度切片密度結(jié)構(gòu)如圖5所示.圖中給出了接收函數(shù)轉(zhuǎn)換密度(圓點(diǎn))和地震層析結(jié)果模型(底圖)之間的橫向密度信息差異.在圖5中,地震層析成像和接收函數(shù)資料有明顯差異,具體表現(xiàn)在馬爾康以西、成都以東、康定以南層析成像轉(zhuǎn)換密度結(jié)果與接收函數(shù)轉(zhuǎn)換密度結(jié)果異常大小及分布相差較大.

    圖5 地震層析成像和接收函數(shù)轉(zhuǎn)換密度結(jié)果(30 km水平切片)實(shí)心圓填充的顏色代表接收函數(shù)轉(zhuǎn)換密度,底圖是層析成像轉(zhuǎn)換密度.Fig.5 Transformed density model from tomography and receive function (horizontal section at 30 km depth)The colors of solid circles represent the transformed density of the receiving function result, and the base map is the transformed density of the tomographic result.

    3.4 重力異常特征

    本研究的三維密度結(jié)構(gòu)反演,布格重力異常選擇最新的WGM2012模型數(shù)據(jù),該模型空間分辨率高達(dá)2′.BGI(Bureau Gravimétrique Internationa)官方給出的WGM2012模型資料顯示,該模型融合了多種重力數(shù)據(jù),同時(shí)使用了高分辨率的ETOPO1高程數(shù)據(jù)進(jìn)行地形改正,考慮了異常計(jì)算過(guò)程中的多種不確定性,評(píng)估顯示在高原地區(qū)的平均精度優(yōu)于5 mGal(Balmino et al.,2012).

    圖6a為我們基于WGM2012模型采用50 km高斯低通濾波得到的該區(qū)布格重力異常.下面將使用該布格重力異常進(jìn)行反演,同時(shí)在研究區(qū)下方0~60 km深度,以水平間隔5′×5′(約為8 km)和垂向間隔4 km的分辨率構(gòu)建三維密度模型空間,該模型初始局部參考約束分別來(lái)自于3.3節(jié)地震層析成像和接收函數(shù)轉(zhuǎn)換密度結(jié)構(gòu)(圖5).圖6b是地震層析成像速度結(jié)果轉(zhuǎn)換得到的密度模型正演得到重力異常.

    圖6a結(jié)果顯示了與深部殼幔界面過(guò)渡的相一致的東西向特征,總的來(lái)看青藏高原東南部川西高原地區(qū)都為負(fù)異常區(qū),四川盆地整體呈現(xiàn)結(jié)構(gòu)清晰的正異常帶,與構(gòu)造邊界線符合良好,即沿龍門(mén)山斷裂帶附近有一條橫貫的重力梯級(jí)帶,其走向與地表斷裂位置符合較好.相比圖6b的速度轉(zhuǎn)換密度結(jié)構(gòu)正演重力異常,布格重力異常(圖6a)在水平向上反映出了更多的短波長(zhǎng)特征,這可能與地殼內(nèi)部介質(zhì)密度橫向結(jié)構(gòu)的變化相關(guān).

    4 反演結(jié)果與分析

    根據(jù)上一節(jié)構(gòu)建的兩個(gè)參考密度模型和布格重力異常,我們完成了貝葉斯同化重力反演.下面我們分別對(duì)從結(jié)果的可靠性和剖面特征兩方面進(jìn)行論述.

    4.1 可靠性分析

    (1)模型殘差

    通過(guò)前文所述的反演方法,我們通過(guò)貝葉斯優(yōu)化,獲得了ABIC最小值對(duì)應(yīng)的反演結(jié)果,模型中四個(gè)超參數(shù)分別為:λ1=139.254、λ2=122.965、α=31.337、β=0.45.圖7中給出了反演模型的正演異常值和其與觀測(cè)值的差異特征,圖7b的異常殘差結(jié)果標(biāo)準(zhǔn)差為±2.5 mGal.

    從圖7a中可以看出,最終密度異常正演所得異常理論值與觀測(cè)剩余值形態(tài)具有較好的一致性,與地震模型正演結(jié)果相比有更多短波長(zhǎng)特征.圖7b中觀測(cè)異常和反演密度模型正演理論異常得到數(shù)據(jù)擬合剩余殘差主要是高頻誤差,標(biāo)準(zhǔn)差略?xún)?yōu)于WGM20121布格異常的5 mGal精度.

    (2)水平結(jié)構(gòu)

    為了進(jìn)一步檢驗(yàn)反演結(jié)果的可靠性,我們仿照?qǐng)D5的結(jié)果,取30 km深度切片,對(duì)比接收函數(shù)轉(zhuǎn)換密度結(jié)果與反演結(jié)果之間的差異,如圖8所示.可以看出:以龍門(mén)山斷裂帶為界(F2),川西高原與四川盆地高低密度分界清晰.與地震層析結(jié)果相比,聯(lián)合反演所得密度分布形態(tài)與接收函數(shù)具有較好的一致性;此外,川西高原顯著低密度異常的背景下出現(xiàn)了星狀分布的小尺度橫向密度差異特征,四川盆地內(nèi)部密度相對(duì)川西呈高密度結(jié)構(gòu)特點(diǎn),也伴隨小尺度相間分布的結(jié)構(gòu)特征.

    圖6 (a) 布格重力異常; (b) 層析成像轉(zhuǎn)化密度模型正演重力異常Fig.6 (a) Bouguer gravity anomalies; (b) Forward modeling gravity anomalies from the topography transformed density model

    圖7 (a) 反演模型結(jié)果正演重力異常; (b) 觀測(cè)異常和反演模型正演異常的差異特征Fig.7 (a) Forward gravity anomalies from the inversion model; (b) The difference between the observed anomalies and forward modeling anomalies from the inversion model

    圖8 基于ABIC重力異常貝葉斯同化反演密度結(jié)果(30 km水平切片)Fig.8 Density inversion result derived from the Bayesian assimilation of gravity anomalies based on ABIC (horizontal slice at 30 km depth)

    圖9 反演模型與參考模型點(diǎn)垂向?qū)Ρ?a) P1點(diǎn);(b) P2點(diǎn);(c) P3點(diǎn);(d) P4點(diǎn);(e) P5點(diǎn);(f) P6點(diǎn).藍(lán)色折線為接收函數(shù)轉(zhuǎn)化密度結(jié)果;紅色折線為層析成像轉(zhuǎn)化密度結(jié)果;黑色空心圓為反演密度結(jié)果,黑色短線為密度估計(jì)誤差.Fig.9 Vertical comparison between the inversion model and the reference model(a) Point 1; (b) Point 2; (c) Point 3; (d) Point 4; (e) Point 5; (f) Point 6. The blue step-lines are the transformed density result from receiving function. The red broken lines are the transformed density result from tomography. The black hollow circles are the inversion density result. The black short lines are the density estimation errors.

    (3)垂直結(jié)構(gòu)

    在垂向分層結(jié)構(gòu)的檢測(cè)方面,我們?cè)诓煌瑯?gòu)造區(qū),分別挑選了6個(gè)接收函數(shù)臺(tái)站位置下方的一維垂向密度結(jié)構(gòu)進(jìn)行比較.圖9a—f給出了接收函數(shù)參考模型(藍(lán)色實(shí)線)、地震層析參考模型(紅色實(shí)線)和同化反演模型(黑點(diǎn))的一維結(jié)果.總體來(lái)說(shuō),三者結(jié)果差異不大,最終反演結(jié)果在不同深度位置綜合了兩種參考模型信息,可說(shuō)明重力同化反演結(jié)果的垂向分辨能力可以從參考模型中獲得.

    具體分析圖9,在上地殼淺部位置,反演結(jié)果與各參考分層差異不大,但在深部中下地殼深度,P3和P6點(diǎn)反演結(jié)構(gòu)更多地參考了接收函數(shù)參考模型結(jié)果,由此可見(jiàn)具體同化重力反演的結(jié)果,并非為簡(jiǎn)單的參考模型平均,而是在實(shí)際重力異常、參考模型特點(diǎn)等先驗(yàn)信息基礎(chǔ)上,給出的最大后驗(yàn)概率估計(jì)優(yōu)選的結(jié)果.

    4.2 剖面特征

    在反演模型結(jié)果綜合對(duì)比分析基礎(chǔ)上,我們進(jìn)一步通過(guò)圖4過(guò)震中位置的四條垂直剖面特征,來(lái)分析本文結(jié)果給出的地殼密度垂向結(jié)構(gòu)特征.圖10分別是AA′、BB′ 、CC′ 和DD′剖面位置(圖4)的密度結(jié)構(gòu),采用相同比例色標(biāo),紅色表示高密度、藍(lán)色表示低密度.從圖10給出的結(jié)果來(lái)看,地殼密度結(jié)構(gòu)縱向分層界面清晰,不同活動(dòng)地塊下密度分成界面構(gòu)造形態(tài)存在顯著差異.龍門(mén)山斷裂帶(F2)下方的密度等值線變化劇烈,可能反映了該區(qū)地殼內(nèi)部密度結(jié)構(gòu)復(fù)雜、變形強(qiáng)烈,該地區(qū)的密度變化特征較好的刻畫(huà)了盆山耦合環(huán)境下的殼幔接觸前緣.對(duì)比圖10a、b和圖10c、d,可發(fā)現(xiàn)地殼橫向密度不均勻特征明顯,剖面下方東西向密度變化劇烈(AA′和BB′剖面),構(gòu)造變形多,結(jié)構(gòu)復(fù)雜;南北向密度變化則相對(duì)較緩(CC′和DD′剖面),殼幔構(gòu)造變形及結(jié)構(gòu)相對(duì)簡(jiǎn)單,這與GPS觀測(cè)得到的青藏高原主要以向東運(yùn)動(dòng)和地表體現(xiàn)的隆升變形特征相一致(張培震等,2008).

    從圖10的局部特征看,殼內(nèi)低密度層刻畫(huà)的也很明顯.例如,鮮水河斷裂帶(F1)、東昆侖斷裂帶(F4)、龍日壩斷裂帶(F4-1)和毛爾蓋斷裂帶(F4-2)下方中地殼內(nèi)存在帶狀低密度層分布(AA′、CC′和DD′剖面).綜合地震層析(雷建設(shè)等, 2009)、接收函數(shù)(鄭晨等,2016)、大地電磁測(cè)深(王緒本等,2018)等結(jié)果在該區(qū)低波速、高泊松比、低電阻率等的認(rèn)識(shí),可將這些低密度層解釋為與中地殼的黏滯性地殼流從高原腹地自北西向南東運(yùn)移過(guò)程有關(guān).

    圖10a、b所示的AA′和BB′剖面結(jié)果顯示,龍門(mén)山斷裂帶(F2)下方地殼密度變化強(qiáng)烈,地層發(fā)生強(qiáng)烈的揉皺變形.從AA′剖面可以看出,汶川地震斷層為高角度斷層,該地震所在的斷層從地面切割向下延伸至約30 km,斷裂形態(tài)與密度分層界面突變跳躍的位置一致.同時(shí),斷層下方約35 km深度處存在中地殼高低密度接觸前緣.從BB′剖面可以看出,蘆山地震同樣發(fā)生在密度分層界面陡變帶向下的延伸面上,蘆山地震斷裂向下延伸面與10~30 km深度的密度差異界面具有較好的吻合特征.

    圖10 沿AA′、BB′ 、CC′ 和DD′密度結(jié)果剖面圖(剖面位置如圖4所示)Fig.10 Density results along profiles AA′, BB′, CC′and DD′ sections(Profile positions are shown in Fig.4)

    5 結(jié)論與討論

    本文基于貝葉斯原理,通過(guò)引入ABIC準(zhǔn)則替代傳統(tǒng)的目標(biāo)函數(shù)最小化方法,給出了一種可以有效同化多個(gè)參考模型和優(yōu)選深度加權(quán)參數(shù)的重力反演新方法.隨后經(jīng)過(guò)多組模型測(cè)試了該方法的可行性、收斂性和有效性.測(cè)試結(jié)果表明,通過(guò)新方法得到的反演模型與真實(shí)模型之間差異性最小,反演結(jié)果合理.最后應(yīng)用該方法對(duì)龍門(mén)山地區(qū)的布格重力異常進(jìn)行反演,其結(jié)果對(duì)于地殼垂向分層、局部密度分布特征及深大斷裂孕震特征都有較好地反映.本文得到的主要研究結(jié)論如下:

    (1)基于綜合模型測(cè)試結(jié)果表明,無(wú)論是在全局參考模型還是在局部參考模型情況下,貝葉斯同化反演算法均可實(shí)現(xiàn)ABIC最小化,并獲得最優(yōu)的參數(shù)估計(jì)結(jié)果.在參考模型較多、模型參數(shù)增多的情況下,不但可以減小人為調(diào)參的難度,還可以更多地吸取多個(gè)參考模型的有效信息,獲得多個(gè)先驗(yàn)約束下的最優(yōu)反演結(jié)果.這對(duì)于在已有先驗(yàn)約束較豐富的地區(qū)開(kāi)展綜合反演研究無(wú)疑是十分必要的.

    (2)從實(shí)際資料的反演測(cè)試結(jié)果看,龍門(mén)山地區(qū)整體地殼密度變化顯著,與南北向密度變化相比,東西向密度變化更加劇烈,構(gòu)造變形更加復(fù)雜.此外下地殼低密度層呈現(xiàn)局部分布特征,結(jié)合地震、大地電磁等研究成果認(rèn)為這些局部分布的下地殼低密度層反映了黏滯性地殼流從高原腹地自北西向南東流入東緣的軌跡,支持了該地區(qū)下地殼流隆升變形假設(shè).

    (3)通過(guò)實(shí)際跨震中剖面特征的分析上來(lái)看,在龍門(mén)山斷裂帶下方,蘆山地震和汶川地震斷層形態(tài)與中上地殼密度分層界面陡變位置具有較好的吻合性.據(jù)此推斷,該區(qū)大地震容易發(fā)生在中上地殼強(qiáng)烈的密度界面陡變位置.與蘆山地震相比,汶川地震斷裂帶下方分布有低密度松潘甘孜地塊中地殼俯沖前緣,該俯沖前緣中地殼低密度層的存在可能是造成不同地震發(fā)震機(jī)制的主要原因.

    綜上所述,本文提出的貝葉斯反演算法,反演過(guò)程不依賴(lài)人的主觀認(rèn)識(shí),完全依靠數(shù)據(jù)說(shuō)話,通過(guò)非線性?xún)?yōu)化算法可以實(shí)現(xiàn)全自動(dòng)化地反演參數(shù)調(diào)節(jié)和模型同化,提取不同先驗(yàn)約束之間的有效信息,可以充分發(fā)揮地震學(xué)方法的垂向分辨能力優(yōu)勢(shì)和重力異常蘊(yùn)含的水平密度結(jié)構(gòu)變化特征.

    在對(duì)強(qiáng)震震源區(qū)的結(jié)構(gòu)研究方面,本文得到的密度剖面結(jié)果顯示,蘆山和汶川震中區(qū)位置都存在顯著的中地殼密度界面陡變特征,其斷層位置與密度分層界面陡變區(qū)都具有較好的一致性.同時(shí)龍門(mén)山斷裂帶下方滑脫層的存在與該斷裂帶下方密度分層界面也有較好的對(duì)應(yīng)(房立華等, 2013; 王夫運(yùn)等, 2015),因此,此類(lèi)高精度的地殼密度結(jié)構(gòu)無(wú)疑更有益于研究大地震孕育環(huán)境的結(jié)構(gòu)與物性特征.

    本文提出的貝葉斯重力反演新算法,是解決多種不同分辨率和時(shí)空尺度先驗(yàn)參考模型同化反演的有效手段,但新算法中仍存在一些難題尚需解決,比如:新算法中多個(gè)超參數(shù)優(yōu)選過(guò)程的計(jì)算量遠(yuǎn)遠(yuǎn)大于傳統(tǒng)單個(gè)超參數(shù)優(yōu)化的計(jì)算量;當(dāng)新算法中優(yōu)化超參數(shù)增多、參考模型之間矛盾較大時(shí),算法可能存在一定的不收斂風(fēng)險(xiǎn).未來(lái),我們將繼續(xù)針對(duì)優(yōu)化反演算法,提升模型計(jì)算能力等方面做進(jìn)一步的改進(jìn)和深入研究.

    致謝感謝中國(guó)科學(xué)臺(tái)陣項(xiàng)目提供的地震數(shù)據(jù),感謝中國(guó)地震局地球物理研究所王興臣博士提供的川滇地殼接收函數(shù)結(jié)果,感謝中國(guó)地震局地球物理研究所李永華研究員和兩位匿名審稿人提供的寶貴建議和幫助.

    猜你喜歡
    參考模型重力反演
    瘋狂過(guò)山車(chē)——重力是什么
    反演對(duì)稱(chēng)變換在解決平面幾何問(wèn)題中的應(yīng)用
    基于低頻軟約束的疊前AVA稀疏層反演
    基于自適應(yīng)遺傳算法的CSAMT一維反演
    仰斜式重力擋土墻穩(wěn)定計(jì)算復(fù)核
    基于環(huán)境的軍事信息系統(tǒng)需求參考模型
    語(yǔ)義網(wǎng)絡(luò)P2P參考模型的查詢(xún)過(guò)程構(gòu)建
    一張紙的承重力有多大?
    疊前同步反演在港中油田的應(yīng)用
    重力異常向上延拓中Poisson積分離散化方法比較
    av黄色大香蕉| 蜜桃久久精品国产亚洲av| 狂野欧美白嫩少妇大欣赏| 国产美女午夜福利| 纯流量卡能插随身wifi吗| 日韩一区二区视频免费看| 岛国毛片在线播放| 欧美日韩在线观看h| 欧美国产精品一级二级三级 | 女人十人毛片免费观看3o分钟| 亚洲精品成人av观看孕妇| 日韩视频在线欧美| 边亲边吃奶的免费视频| 国产精品欧美亚洲77777| 婷婷色综合www| 久久精品国产亚洲av天美| 久久久久国产网址| 日本与韩国留学比较| 一级a做视频免费观看| 亚洲av免费高清在线观看| 三级国产精品片| 人妻 亚洲 视频| 国产av精品麻豆| 久久久亚洲精品成人影院| 老师上课跳d突然被开到最大视频| 午夜免费男女啪啪视频观看| 中文字幕制服av| 嘟嘟电影网在线观看| 国产精品一二三区在线看| 亚洲精品成人av观看孕妇| 亚洲国产色片| 18禁在线无遮挡免费观看视频| 男女国产视频网站| www.色视频.com| 高清黄色对白视频在线免费看 | 美女cb高潮喷水在线观看| 国产淫片久久久久久久久| 香蕉精品网在线| 中文字幕久久专区| 国产av精品麻豆| 国产精品国产三级专区第一集| 精品一区二区免费观看| 久久久久久久精品精品| 久久精品夜色国产| 精品久久久噜噜| 免费看光身美女| 国产免费福利视频在线观看| 久久99热这里只有精品18| 欧美另类一区| 毛片女人毛片| 亚洲国产精品国产精品| 日韩,欧美,国产一区二区三区| 国产免费福利视频在线观看| 99久久中文字幕三级久久日本| 久久久久人妻精品一区果冻| 亚洲国产精品成人久久小说| 久久久a久久爽久久v久久| 人妻一区二区av| 免费观看在线日韩| 亚洲av成人精品一二三区| 一级毛片电影观看| 国产精品久久久久久精品电影小说 | 亚洲丝袜综合中文字幕| 欧美三级亚洲精品| 色婷婷av一区二区三区视频| 91久久精品国产一区二区成人| 国产 精品1| 精品亚洲乱码少妇综合久久| 黄色一级大片看看| av国产免费在线观看| 日本wwww免费看| 91久久精品国产一区二区三区| 男女啪啪激烈高潮av片| 男女下面进入的视频免费午夜| 欧美另类一区| 国产美女午夜福利| 赤兔流量卡办理| 黄色一级大片看看| 国产精品成人在线| 精品久久国产蜜桃| 国产精品嫩草影院av在线观看| 国产一区有黄有色的免费视频| 在线精品无人区一区二区三 | 高清欧美精品videossex| 香蕉精品网在线| 在线观看免费视频网站a站| 国产精品一区二区在线观看99| 亚洲真实伦在线观看| 亚洲自偷自拍三级| 日本一二三区视频观看| 乱码一卡2卡4卡精品| 在线精品无人区一区二区三 | 国模一区二区三区四区视频| 亚洲美女搞黄在线观看| 国产av一区二区精品久久 | 成人国产av品久久久| 丝袜喷水一区| h日本视频在线播放| 七月丁香在线播放| 午夜福利视频精品| 成人美女网站在线观看视频| 日韩在线高清观看一区二区三区| 国产成人freesex在线| 在线天堂最新版资源| 国产熟女欧美一区二区| 日韩中字成人| 亚洲aⅴ乱码一区二区在线播放| 欧美xxxx性猛交bbbb| 国产色爽女视频免费观看| 九草在线视频观看| 久久精品国产自在天天线| 狂野欧美激情性xxxx在线观看| 2021少妇久久久久久久久久久| 欧美极品一区二区三区四区| 日日摸夜夜添夜夜添av毛片| 日产精品乱码卡一卡2卡三| 欧美精品国产亚洲| 一本久久精品| 最近最新中文字幕大全电影3| 色综合色国产| 又大又黄又爽视频免费| 美女福利国产在线 | 99国产精品免费福利视频| 亚洲激情五月婷婷啪啪| av免费在线看不卡| 毛片女人毛片| 国产色爽女视频免费观看| 精品酒店卫生间| av在线观看视频网站免费| 国产精品欧美亚洲77777| 国内精品宾馆在线| 在线亚洲精品国产二区图片欧美 | 免费少妇av软件| 亚洲国产精品国产精品| 日韩在线高清观看一区二区三区| 老司机影院毛片| 自拍偷自拍亚洲精品老妇| 亚洲伊人久久精品综合| 一区二区av电影网| 夫妻午夜视频| 免费观看性生交大片5| 欧美性感艳星| 一级片'在线观看视频| 啦啦啦在线观看免费高清www| 18禁裸乳无遮挡动漫免费视频| 777米奇影视久久| 一本色道久久久久久精品综合| 国产成人精品婷婷| 中国三级夫妇交换| 亚洲精品国产色婷婷电影| 日韩av免费高清视频| 最近中文字幕高清免费大全6| 男人爽女人下面视频在线观看| 99热全是精品| 丝瓜视频免费看黄片| 在线观看一区二区三区激情| 永久网站在线| 亚洲内射少妇av| 51国产日韩欧美| 联通29元200g的流量卡| 1000部很黄的大片| 少妇人妻 视频| 中文乱码字字幕精品一区二区三区| 国产色爽女视频免费观看| 精品国产露脸久久av麻豆| 久久精品久久久久久噜噜老黄| 久久99热这里只频精品6学生| 亚洲成人一二三区av| 久久久精品免费免费高清| 国产一区二区三区av在线| 一级二级三级毛片免费看| 成年免费大片在线观看| av天堂中文字幕网| 一级毛片 在线播放| 国产精品久久久久久av不卡| 一级爰片在线观看| 久久热精品热| 美女xxoo啪啪120秒动态图| 婷婷色综合www| 美女中出高潮动态图| 高清黄色对白视频在线免费看 | 有码 亚洲区| 一本一本综合久久| 日韩一区二区视频免费看| 交换朋友夫妻互换小说| 久久久久人妻精品一区果冻| 亚洲欧美日韩另类电影网站 | 国产女主播在线喷水免费视频网站| 日韩三级伦理在线观看| 亚洲欧美中文字幕日韩二区| 欧美日本视频| 搡女人真爽免费视频火全软件| 卡戴珊不雅视频在线播放| 亚洲人成网站在线播| 国产精品.久久久| 精品国产三级普通话版| 精品久久久精品久久久| 久久婷婷青草| 亚洲av国产av综合av卡| 久久久久久久精品精品| 99久久人妻综合| 成人免费观看视频高清| 18禁裸乳无遮挡动漫免费视频| 大香蕉97超碰在线| 小蜜桃在线观看免费完整版高清| 成年美女黄网站色视频大全免费 | 91久久精品国产一区二区成人| 精品久久久久久久末码| 最后的刺客免费高清国语| 日韩 亚洲 欧美在线| 51国产日韩欧美| av线在线观看网站| 一级毛片aaaaaa免费看小| 黄色配什么色好看| 日韩欧美 国产精品| 街头女战士在线观看网站| 久久久久久伊人网av| 一级a做视频免费观看| 亚洲婷婷狠狠爱综合网| 51国产日韩欧美| a级毛色黄片| av黄色大香蕉| 99视频精品全部免费 在线| av免费在线看不卡| 亚洲欧美一区二区三区国产| av.在线天堂| 免费久久久久久久精品成人欧美视频 | 日日撸夜夜添| 亚洲欧美一区二区三区国产| 一个人看视频在线观看www免费| 在线观看av片永久免费下载| 亚洲精品国产色婷婷电影| 亚洲精品第二区| 亚洲三级黄色毛片| 国产真实伦视频高清在线观看| 丰满迷人的少妇在线观看| 亚洲av综合色区一区| av黄色大香蕉| 久久婷婷青草| 人妻一区二区av| 亚洲精品视频女| 亚洲美女视频黄频| 激情 狠狠 欧美| 一级毛片我不卡| 亚洲精品第二区| av专区在线播放| 高清午夜精品一区二区三区| 国产爱豆传媒在线观看| 亚洲欧洲国产日韩| 国产成人免费观看mmmm| 日韩中文字幕视频在线看片 | 亚洲不卡免费看| av专区在线播放| av天堂中文字幕网| 亚洲国产高清在线一区二区三| 欧美日韩视频精品一区| 国产亚洲午夜精品一区二区久久| 亚洲综合精品二区| 韩国av在线不卡| 成年女人在线观看亚洲视频| 国产精品一区二区在线观看99| 国产精品无大码| 国产黄片视频在线免费观看| 最近最新中文字幕大全电影3| 日本色播在线视频| 婷婷色综合www| 成人午夜精彩视频在线观看| 26uuu在线亚洲综合色| 午夜老司机福利剧场| 男人添女人高潮全过程视频| 深夜a级毛片| 欧美国产精品一级二级三级 | 日韩视频在线欧美| 午夜福利高清视频| 国产精品爽爽va在线观看网站| 亚洲国产成人一精品久久久| 国产精品久久久久久久电影| 美女xxoo啪啪120秒动态图| 国产91av在线免费观看| 精品人妻偷拍中文字幕| 日本av手机在线免费观看| 新久久久久国产一级毛片| 国产永久视频网站| 观看美女的网站| av在线app专区| 日韩一本色道免费dvd| 黄色怎么调成土黄色| av黄色大香蕉| 91aial.com中文字幕在线观看| 国产成人a∨麻豆精品| 中文字幕久久专区| 免费播放大片免费观看视频在线观看| 丝袜脚勾引网站| 国产伦在线观看视频一区| 亚洲伊人久久精品综合| 熟女人妻精品中文字幕| 黄色配什么色好看| 亚洲国产欧美人成| 午夜精品国产一区二区电影| 午夜福利在线观看免费完整高清在| 日韩强制内射视频| 18禁动态无遮挡网站| av.在线天堂| 在线观看国产h片| 亚洲图色成人| 欧美精品一区二区免费开放| 插逼视频在线观看| a级毛片免费高清观看在线播放| 纯流量卡能插随身wifi吗| 亚洲欧美一区二区三区黑人 | 黑人猛操日本美女一级片| 亚洲av综合色区一区| 国产精品久久久久久久电影| 久久99热这里只频精品6学生| 99热这里只有是精品50| 少妇丰满av| av女优亚洲男人天堂| 欧美人与善性xxx| 亚洲精品乱久久久久久| 干丝袜人妻中文字幕| 欧美日韩综合久久久久久| 多毛熟女@视频| 久久久久久久国产电影| av网站免费在线观看视频| 日本av手机在线免费观看| 国产日韩欧美在线精品| 国产午夜精品久久久久久一区二区三区| 亚洲婷婷狠狠爱综合网| 久久久欧美国产精品| 成人国产麻豆网| 亚洲欧美一区二区三区国产| 网址你懂的国产日韩在线| 日本猛色少妇xxxxx猛交久久| 嘟嘟电影网在线观看| 亚洲av日韩在线播放| 国产综合精华液| 在线观看一区二区三区激情| av免费观看日本| 久久国内精品自在自线图片| 99re6热这里在线精品视频| 日韩大片免费观看网站| 三级国产精品片| 在线天堂最新版资源| 丰满迷人的少妇在线观看| 老女人水多毛片| 日本av免费视频播放| 久久久色成人| 不卡视频在线观看欧美| 国产亚洲欧美精品永久| 91精品国产九色| 国产熟女欧美一区二区| 成人影院久久| 老女人水多毛片| 亚洲熟女精品中文字幕| 在线免费十八禁| 五月开心婷婷网| 伦理电影免费视频| 免费少妇av软件| 又大又黄又爽视频免费| 热re99久久精品国产66热6| 国产伦理片在线播放av一区| 啦啦啦在线观看免费高清www| 中文字幕制服av| 网址你懂的国产日韩在线| 我要看日韩黄色一级片| 妹子高潮喷水视频| 成人毛片a级毛片在线播放| 亚洲高清免费不卡视频| 免费看不卡的av| 欧美亚洲 丝袜 人妻 在线| 免费看av在线观看网站| 在线观看一区二区三区激情| av在线蜜桃| 最近中文字幕2019免费版| 亚洲天堂av无毛| 色哟哟·www| 久久国内精品自在自线图片| 国精品久久久久久国模美| 涩涩av久久男人的天堂| 伦精品一区二区三区| 亚洲国产av新网站| 国产精品久久久久久av不卡| 97在线人人人人妻| 熟女电影av网| 国产亚洲最大av| 国产中年淑女户外野战色| 久久av网站| 波野结衣二区三区在线| 街头女战士在线观看网站| 久久99热这里只频精品6学生| 国产精品99久久99久久久不卡 | 只有这里有精品99| 国产视频首页在线观看| 韩国av在线不卡| 亚洲在久久综合| 亚洲欧美日韩无卡精品| 日本av手机在线免费观看| 嘟嘟电影网在线观看| h日本视频在线播放| 国产成人免费观看mmmm| 又粗又硬又长又爽又黄的视频| 五月玫瑰六月丁香| 中文字幕精品免费在线观看视频 | 国产日韩欧美在线精品| a级一级毛片免费在线观看| 午夜福利高清视频| 亚洲精品国产色婷婷电影| 亚洲四区av| 夜夜爽夜夜爽视频| 亚洲av日韩在线播放| 精品一区二区免费观看| 久久99热这里只频精品6学生| 一级av片app| 国产探花极品一区二区| 免费人妻精品一区二区三区视频| 亚洲精品国产色婷婷电影| 老司机影院成人| 久久国产精品男人的天堂亚洲 | 1000部很黄的大片| 国内少妇人妻偷人精品xxx网站| 成人国产av品久久久| 欧美+日韩+精品| 大香蕉久久网| 国产精品无大码| 国产精品三级大全| 99久国产av精品国产电影| 一级av片app| 日日摸夜夜添夜夜添av毛片| 免费高清在线观看视频在线观看| 18禁在线播放成人免费| 九九爱精品视频在线观看| 国产精品久久久久久久电影| 又粗又硬又长又爽又黄的视频| 国产精品99久久99久久久不卡 | 特大巨黑吊av在线直播| 波野结衣二区三区在线| 精品久久久久久久久亚洲| 国产亚洲欧美精品永久| 亚洲成人手机| 99久久人妻综合| 伊人久久精品亚洲午夜| 天天躁夜夜躁狠狠久久av| 成人国产麻豆网| 青春草亚洲视频在线观看| 性色avwww在线观看| 亚洲av不卡在线观看| 国产精品国产av在线观看| 制服丝袜香蕉在线| 成人国产麻豆网| 国产黄片视频在线免费观看| 国产淫片久久久久久久久| 亚洲三级黄色毛片| 精品国产一区二区三区久久久樱花 | 精品人妻偷拍中文字幕| 亚洲av男天堂| 亚洲欧美清纯卡通| 亚洲av在线观看美女高潮| 国产精品一区二区三区四区免费观看| 成人黄色视频免费在线看| 91精品伊人久久大香线蕉| 插阴视频在线观看视频| 久久久久久久国产电影| 久久鲁丝午夜福利片| 国产深夜福利视频在线观看| 激情五月婷婷亚洲| av卡一久久| 91精品伊人久久大香线蕉| 97精品久久久久久久久久精品| 老师上课跳d突然被开到最大视频| 精品一区二区三卡| 最近中文字幕2019免费版| av播播在线观看一区| 亚洲精品日韩av片在线观看| 天堂8中文在线网| 国产精品一二三区在线看| 老女人水多毛片| 日本-黄色视频高清免费观看| 免费播放大片免费观看视频在线观看| 久久青草综合色| 日韩免费高清中文字幕av| 中文乱码字字幕精品一区二区三区| 国产精品伦人一区二区| 午夜视频国产福利| 亚洲激情五月婷婷啪啪| 日本猛色少妇xxxxx猛交久久| av在线老鸭窝| av天堂中文字幕网| 婷婷色综合大香蕉| 美女福利国产在线 | 久久久色成人| av视频免费观看在线观看| 日韩国内少妇激情av| 自拍偷自拍亚洲精品老妇| 亚洲精品国产成人久久av| 国产精品一二三区在线看| 人人妻人人爽人人添夜夜欢视频 | 身体一侧抽搐| 高清av免费在线| 亚洲欧美日韩卡通动漫| 久久久久久久国产电影| 欧美最新免费一区二区三区| 少妇人妻久久综合中文| 一本久久精品| 久久久久久人妻| 夜夜看夜夜爽夜夜摸| 亚洲一级一片aⅴ在线观看| 久久人人爽人人片av| 亚洲欧美日韩东京热| 99热国产这里只有精品6| 久久韩国三级中文字幕| 久久久久久久大尺度免费视频| 国产成人精品福利久久| 欧美日韩精品成人综合77777| 一边亲一边摸免费视频| 乱系列少妇在线播放| av视频免费观看在线观看| 91aial.com中文字幕在线观看| 女人十人毛片免费观看3o分钟| 国产高清三级在线| 日本爱情动作片www.在线观看| 免费看光身美女| 国产女主播在线喷水免费视频网站| 黄色日韩在线| 久久精品国产自在天天线| 在线观看国产h片| 国产精品久久久久成人av| 妹子高潮喷水视频| 国产在线免费精品| 久久久久视频综合| 国产精品伦人一区二区| 寂寞人妻少妇视频99o| 免费观看无遮挡的男女| 久久亚洲国产成人精品v| 国产欧美另类精品又又久久亚洲欧美| 国产精品久久久久久久久免| 亚洲精品乱码久久久久久按摩| a 毛片基地| 国产有黄有色有爽视频| 人人妻人人添人人爽欧美一区卜 | 久久久亚洲精品成人影院| 中文乱码字字幕精品一区二区三区| 精品熟女少妇av免费看| 啦啦啦视频在线资源免费观看| 不卡视频在线观看欧美| 国产黄片视频在线免费观看| 麻豆成人av视频| 91久久精品国产一区二区三区| 乱码一卡2卡4卡精品| av免费观看日本| 国产黄色视频一区二区在线观看| 日韩欧美 国产精品| 在线 av 中文字幕| 国产有黄有色有爽视频| 啦啦啦中文免费视频观看日本| h视频一区二区三区| 大片电影免费在线观看免费| 亚洲欧美日韩卡通动漫| 久久久久久久精品精品| 少妇的逼水好多| 黄片无遮挡物在线观看| 成人国产麻豆网| av在线app专区| 各种免费的搞黄视频| 黄色欧美视频在线观看| 亚洲欧美日韩无卡精品| 国产片特级美女逼逼视频| 夜夜看夜夜爽夜夜摸| 王馨瑶露胸无遮挡在线观看| 国产有黄有色有爽视频| 国产亚洲精品久久久com| 我的老师免费观看完整版| 国产精品一区www在线观看| 汤姆久久久久久久影院中文字幕| 99热全是精品| 亚洲天堂av无毛| 国产91av在线免费观看| 久久99热这里只有精品18| 中文字幕人妻熟人妻熟丝袜美| 成人国产av品久久久| 又大又黄又爽视频免费| 国产黄色视频一区二区在线观看| 三级经典国产精品| 久久久久久久精品精品| 精品久久久噜噜| 国产精品免费大片| 久久久久久久亚洲中文字幕| 久久人人爽人人爽人人片va| 久热这里只有精品99| 色哟哟·www| 国产精品人妻久久久久久| 舔av片在线| 欧美激情国产日韩精品一区| 五月玫瑰六月丁香| 一级毛片aaaaaa免费看小| 亚洲国产欧美在线一区| 午夜免费观看性视频| 一级毛片电影观看| 久久毛片免费看一区二区三区| 狂野欧美激情性xxxx在线观看| 日韩在线高清观看一区二区三区| 国产精品女同一区二区软件| 亚洲欧美日韩东京热| 免费观看无遮挡的男女| h日本视频在线播放| 最后的刺客免费高清国语| 免费观看在线日韩| 亚洲国产精品成人久久小说| 亚洲av中文av极速乱| 国产精品久久久久久精品古装| 少妇被粗大猛烈的视频| 日日啪夜夜爽| 色吧在线观看| 在线观看免费视频网站a站| 噜噜噜噜噜久久久久久91| 99国产精品免费福利视频|