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

李紅蕾， 陳石*， 莊建倉(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ì)的重要工具.由于地球物理反演問題先天具有多解性，通常反演結(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í)越來越多.如何充分融合已有的模型，發(fā)揮各種地球物理手段的優(yōu)勢(shì)，減小反演結(jié)果的不確定性，無疑是地球物理聯(lián)合反演研究要解決的核心問題(Welford et al., 2018; Yang et al., 2012).

如果將所有已有模型解釋結(jié)果看作先驗(yàn)認(rèn)識(shí)，當(dāng)今地學(xué)家面臨的首要挑戰(zhàn)是如何從已有模型數(shù)據(jù)中提取盡可能多的有用信息以獲取新的見解.與常規(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)中，通過對(duì)系統(tǒng)中不確定性的估計(jì)和控制，來提高后驗(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).除了通過增加觀測(cè)數(shù)據(jù)和減小觀測(cè)誤差，在一定程度上降低反演的多解性.之外，在反演中加入先驗(yàn)信息進(jìn)行約束是更有效的方法(Li and Oldenburg，1998；陳石等，2010；李紅蕾等，2016).近年來，重力約束條件越來越受到重視，許多地球物理學(xué)者都提出了各種各樣的約束條件以及相應(yīng)的反演方法.約束總體上可以分為兩類，即數(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)的三維重力反演算法中，上述約束和參考都是通過阻尼因子被引入到反演方程中，并通過廣義交叉驗(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ǐng)域的數(shù)據(jù)同化問題.這種方法已成為數(shù)據(jù)同化的重要手段，其是以數(shù)據(jù)驅(qū)動(dòng)模式量化反演參數(shù)，通過最小化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)域.

本文探索通過數(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)Ｐ?，測(cè)試了該方法在解決多種已知先驗(yàn)參考約束和深度加權(quán)約束中超參數(shù)和后驗(yàn)?zāi)Ｐ凸烙?jì)問題的優(yōu)化效果，驗(yàn)證了該方法的可行性和有效性；第4部分和第5部分是應(yīng)用該方法和龍門山及周邊地區(qū)2′×2′的高精度WGM2012布格重力異常數(shù)據(jù)反演了該區(qū)的地殼密度模型, 并對(duì)龍門山地區(qū)地殼內(nèi)低密度分布、隆升變形機(jī)制、強(qiáng)震震源區(qū)環(huán)境等進(jìn)行了深入分析；最后，第6部分對(duì)本文研究方法和認(rèn)識(shí)進(jìn)行了總結(jié)和討論.

1 方法原理

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

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則稱為重力反演問題，一般觀測(cè)N?M，直接反演在數(shù)學(xué)上是沒有唯一解的.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)反演方法通過式(3)將其加入到反演目標(biāo)函數(shù)中：

(3)

為解決上述問題，我們提出了重力貝葉斯同化反演算法，其是在綜合考慮重力觀測(cè)數(shù)據(jù)誤差及模型約束不確定性的基礎(chǔ)上，通過引入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ù)，又稱為反演超參數(shù)，其可以通過引入ABIC準(zhǔn)則來求解，即最小化統(tǒng)計(jì)量：

(8)

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

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

2.1 參考模型約束

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

這個(gè)數(shù)值實(shí)驗(yàn)中，真實(shí)模型在y方向沒有變化，在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è)試過程中，我們選擇了四種不同的反演方案，方案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é)果差異性顯著，詳見圖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è)無論在界面起伏和上下密度差方面都有一定程度改善，而與圖2d的貝葉斯同化聯(lián)合反演結(jié)果相比，在上下界面位置和密度差方面圖2d明顯優(yōu)于圖2c.

在該模型測(cè)試中，參考模型1和2的歸一化權(quán)重參數(shù)分別為λ1和λ2.圖2a—d的反演結(jié)果，對(duì)應(yīng)的參數(shù)詳見表1.從表1中的結(jié)果可以看出，在圖2d的最優(yōu)化過程中，參考模型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é)果來看，參考模型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é)合來看，與單個(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ù)來選擇提取不同參考之間的有用信息，并實(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)反演的趨膚效應(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方向沒有變化，抽取模型在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ù)方法可以提供此類的局部模型信息.

該模型測(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ù)α值過大時(shí)，得到的反演結(jié)果圖3b中高密度部分主要集中在底部，異常體的深度與真實(shí)模型埋深存在一定的差異；而當(dāng)加權(quán)參數(shù)β過大時(shí)，得到的圖3c的反演結(jié)果其高密度主要集中在淺部，異常體埋深同樣存在較大差異.通過本文算法得到的最優(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) 過衰減加權(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)參考模型較多的龍門山地區(qū)進(jìn)行實(shí)際數(shù)據(jù)測(cè)試.

3 龍門山地殼模型

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

龍門山地區(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)的活斷層：鮮水河斷裂帶、龍門山斷裂帶、東昆侖斷裂帶、龍泉山斷裂帶、龍日壩斷裂帶、毛爾蓋分支斷層等(徐錫偉等，2008).此外，沿著龍門山斷裂帶，還存在著約5 km的強(qiáng)烈高程梯度帶(如圖4所示).

圖4 龍門山地區(qū)主要構(gòu)造特征與地震活動(dòng)紅色實(shí)線為斷裂，四條黑色實(shí)線為跨震中研究剖面，黃色圈為5級(jí)以上地震位置，白色實(shí)心圓為中國(guó)地震觀測(cè)網(wǎng)臺(tái)站位置，F(xiàn)1：鮮水河斷裂帶；F2：龍門山斷裂帶；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ó)地震科學(xué)臺(tái)陣項(xiàng)目為代表的大規(guī)模地球物理觀測(cè)相繼開展，已積累了大量的地球物理探測(cè)工作成果(如Yao et al., 2008; Li et al., 2011; Zhang et al., 2011; 王緒本等，2018).這些深部成果對(duì)該區(qū)的殼幔結(jié)構(gòu)及其相關(guān)動(dòng)力學(xué)問題達(dá)成了部分共識(shí)，但在一些基本問題上依舊存在爭(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)等多種模式.

此外，龍門山地區(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)通過地震層析結(jié)果，認(rè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é)過程的深入認(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é)過程爭(zhēng)議提供很好的論據(jù)(王椿鏞等, 2016).

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

眾所周知，重力數(shù)據(jù)水平分辨率高，但垂向分辨率低，在實(shí)際反演過程中，要想得到有地質(zhì)意義的結(jié)果，還需要添加深部參考約束.相對(duì)重力方法來講，地震成像方法具有較好的垂向分辨率，重力反演中將地震成像結(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)該回答的問題.我們將以多種不同類型的地震速度轉(zhuǎn)換密度作為已知先驗(yàn)參考，利用高精度的衛(wèi)星重力場(chǎng)模型數(shù)據(jù)，通過新算法實(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ì)算過程中的多種不確定性，評(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)建三維密度模型空間，該模型初始局部參考約束分別來自于3.3節(jié)地震層析成像和接收函數(shù)轉(zhuǎn)換密度結(jié)構(gòu)(圖5).圖6b是地震層析成像速度結(jié)果轉(zhuǎn)換得到的密度模型正演得到重力異常.

圖6a結(jié)果顯示了與深部殼幔界面過渡的相一致的東西向特征，總的來看青藏高原東南部川西高原地區(qū)都為負(fù)異常區(qū)，四川盆地整體呈現(xiàn)結(jié)構(gòu)清晰的正異常帶，與構(gòu)造邊界線符合良好，即沿龍門山斷裂帶附近有一條橫貫的重力梯級(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)模型殘差

通過前文所述的反演方法，我們通過貝葉斯優(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)差略優(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所示.可以看出：以龍門山斷裂帶為界(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é)果.總體來說，三者結(jié)果差異不大，最終反演結(jié)果在不同深度位置綜合了兩種參考模型信息，可說明重力同化反演結(jié)果的垂向分辨能力可以從參考模型中獲得.

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

4.2 剖面特征

在反演模型結(jié)果綜合對(duì)比分析基礎(chǔ)上，我們進(jìn)一步通過圖4過震中位置的四條垂直剖面特征，來分析本文結(jié)果給出的地殼密度垂向結(jié)構(gòu)特征.圖10分別是AA′、BB′ 、CC′ 和DD′剖面位置(圖4)的密度結(jié)構(gòu)，采用相同比例色標(biāo)，紅色表示高密度、藍(lán)色表示低密度.從圖10給出的結(jié)果來看，地殼密度結(jié)構(gòu)縱向分層界面清晰，不同活動(dòng)地塊下密度分成界面構(gòu)造形態(tài)存在顯著差異.龍門山斷裂帶(F2)下方的密度等值線變化劇烈，可能反映了該區(qū)地殼內(nèi)部密度結(jié)構(gòu)復(fù)雜、變形強(qiáng)烈，該地區(qū)的密度變化特征較好的刻畫了盆山耦合環(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)低密度層刻畫的也很明顯.例如，鮮水河斷裂帶(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ān).

圖10a、b所示的AA′和BB′剖面結(jié)果顯示，龍門山斷裂帶(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é)論與討論

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

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

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

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

綜上所述，本文提出的貝葉斯反演算法，反演過程不依賴人的主觀認(rèn)識(shí)，完全依靠數(shù)據(jù)說話，通過非線性優(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í)龍門山斷裂帶下方滑脫層的存在與該斷裂帶下方密度分層界面也有較好的對(duì)應(yīng)(房立華等, 2013; 王夫運(yùn)等, 2015)，因此，此類高精度的地殼密度結(jié)構(gòu)無疑更有益于研究大地震孕育環(huán)境的結(jié)構(gòu)與物性特征.

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

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