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

李紅蕾， 陳石*， 莊建倉， 張貝， 石磊

1 中國地震局地球物理研究所， 北京 100081 2 北京白家疃地球科學國家野外科學觀測研究站， 北京 100095 3 日本數理統計研究所， 東京 190-8562

0 引言

地球物理反演技術是研究地殼內部結構和性質的重要工具.由于地球物理反演問題先天具有多解性，通常反演結果在數學上適定，但可能與實際地質情況并不相符(Li and Oldenburg, 1998; Williams et al., 2004; Howe, 2009; Dirian et al.,2016).現今在很多研究熱點地區，利用多種手段獲得的地殼結構模型和認識越來越多.如何充分融合已有的模型，發揮各種地球物理手段的優勢，減小反演結果的不確定性，無疑是地球物理聯合反演研究要解決的核心問題(Welford et al., 2018; Yang et al., 2012).

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

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

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

本文探索通過數據同化手段，將貝葉斯同化策略引入到傳統的重力約束反演中，設計一種可以融合多種先驗推斷的重力異常貝葉斯同化反演算法，以期得到更加準確的最大后驗概率意義下的模型參數估計.文中第2部分給出了三維重力貝葉斯同化反演的求解方程和ABIC目標函數的構建；第3部分設計了兩個綜合實驗模型，測試了該方法在解決多種已知先驗參考約束和深度加權約束中超參數和后驗模型估計問題的優化效果，驗證了該方法的可行性和有效性；第4部分和第5部分是應用該方法和龍門山及周邊地區2′×2′的高精度WGM2012布格重力異常數據反演了該區的地殼密度模型, 并對龍門山地區地殼內低密度分布、隆升變形機制、強震震源區環境等進行了深入分析；最后，第6部分對本文研究方法和認識進行了總結和討論.

1 方法原理

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

d=Gm，

(1)

其中，G為N×M的二維核函數矩陣或格林函數矩陣，矩陣中每個元素值與模型中待計算密度單元和觀測點的位置相關，m為M×1的一維密度單元矩陣，d為N×1的一維觀測矩陣.方程(1)中已知異常d求解m則稱為重力反演問題，一般觀測N?M，直接反演在數學上是沒有唯一解的.Tikhonov and Arsenin(1977)引入最小模約束后，求解方程(1)可以變為：

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

(2)

其中λ是阻尼因子.

公式(2)具有數學上的唯一解.當有最小構造、光滑度、深度加權、參考模型等多個約束存在時，傳統反演方法通過式(3)將其加入到反演目標函數中：

(3)

為解決上述問題，我們提出了重力貝葉斯同化反演算法，其是在綜合考慮重力觀測數據誤差及模型約束不確定性的基礎上，通過引入ABIC準則實現最大化邊際概率分布的數據及模型權重超參數自動優選的反演方法.算法的基本原理如下：給定特定模型下重力觀測數據的條件概率分布和模型的先驗概率分布，根據貝葉斯公式，結合了觀測數據，模型的后驗概率密度函數(pdf)表示為：

(4)

其中，p(*|*)為條件概率密度函數.如果數據和模型誤差都以正態分布，則重力觀測數據的條件概率為：

(5)

考慮最小模型和多個先驗約束的模型先驗概率分布可表示為：

(6)

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

(7)

式中λ0、λms、λri為數據和模型權重系數，又稱為反演超參數，其可以通過引入ABIC準則來求解，即最小化統計量：

(8)

2 實驗測試

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

2.1 參考模型約束

引入適當的參考模型約束能在垂向分層上改善重力異常反演結果.通常認為參考模型越準確，反演結果越接近真實情況.但在實際問題中，可參考的模型不止一個，其誤差、不確定性差異都不盡相同，同時不同參考模型之間也往往不一致.對于多參考模型約束問題，特別是在先驗參考模型誤差未知情況下，如何分配每個參考模型的權重問題一直是研究熱點(Hansen, 1994; Vogel, 1996).下面就本文提出的方法如何分配權重和改善反演結果進行測試.

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

在測試中，設計了兩個簡單的參考模型，每個模型與真實模型都存在一定的誤差.其中，參考模型M1密度起伏界分界面形態與真實模型一致，但界面兩側的密度差比真實模型小10%，如圖1c所示；參考模型M2中的密度起伏界面位置與真實模型不一致，界面兩側的密度差真實模型大10%，如圖1d所示.此外，參考模型1和2中還分別加入了2%和1%高斯誤差.

在模型測試過程中，我們選擇了四種不同的反演方案，方案1是僅有參考模型M1優化約束反演；方案2是僅有參考約束模型M2優化約束反演；方案3是參考模型M1和M2聯合非優化約束反演；方案4是參考模型M1和M2聯合優化約束反演.每種反演都實現了模型異常與數據異常的擬合，但反演的模型結果差異性顯著，詳見圖2所示.

圖1 不同參考模型約束下反演算法測試(a) 真實模型； (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 不同參考模型約束下反演結果(a) 參考模型M1約束反演結果; (b) 參考模型M2約束反演結果; (c) 參考模型M1和M2聯合非優化約束反演結果; (d) 參考模型M1和M2聯合優化約束反演結果，白色實線為真實模型密度分界面.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中的結果可以看出，圖2a的反演模型相比圖2b結果界面更清晰，但上下密度差偏差較大，圖2b的界面凸起部分效果比圖2a差.圖2c的聯合反演結果，相比前兩個無論在界面起伏和上下密度差方面都有一定程度改善，而與圖2d的貝葉斯同化聯合反演結果相比，在上下界面位置和密度差方面圖2d明顯優于圖2c.

在該模型測試中，參考模型1和2的歸一化權重參數分別為λ1和λ2.圖2a—d的反演結果，對應的參數詳見表1.從表1中的結果可以看出，在圖2d的最優化過程中，參考模型1的權重更大，對應的ABIC值也最小.

表1 不同參考約束下反演統計參數及ABIC參數特征Table 1 Inversion statistical parameters and ABIC parameter characteristics under different reference constraints

結合表1中的反演統計參數結果來看，參考模型1約束反演優化得到的超參數λ1為129.055，以此計算得到的ABIC值為-1355.069；參考模型2約束反演優化得到的超參數λ2為160.062，計算得到的ABIC值為-1406.770；將上述兩個模型單獨優化得到的超參數直接代入多模型約束，計算得到的ABIC值為-3217.861；將兩個模型同時優化反演得到的最優λ1和λ2分別為1798.5和1094.8，最小化ABIC值為-3700.195.

將表1反演統計參數與圖2中的反演結果結合來看，與單個參考模型和非優化權重的多個參考模型約束反演相比，ABIC最小化(ABIC最小值為-3700.195)給出的最優超參數可有效降低反演的卡方值，提高反演的準確性和精度.與手動指定的模型約束權重進行反演計算得到的反演結果(圖2c)相比，基于優化ABIC的反演算法可依據觀測數據來選擇提取不同參考之間的有用信息，并實現不同參考模型與反演結果之間的數據同化.

2.2 深度加權約束

在式(1)的核函數矩陣G中，每個元素數值大小都與觀測點和模型單元之間的距離平方成正比，在不同深度位置上，模型單元的核函數大小數值差異明顯，淺層單元的數值遠大于深層.體現在反演結果中，常常會出現異常的變化僅集中在淺部模型單元上，通常稱之為重力位場反演的趨膚效應(Li and Oldenburg,1998).

(9)

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

本節的測試模型如圖3所示，與圖2模型在y方向設置相同，即在y方向沒有變化，抽取模型在y=0處的垂直切片如圖3a所示.真實模型在深度方向由兩個規則長方體組成，周邊藍色區域密度為零，每個長方體的密度大小同為0.2 g·cm-3，長方體的頂面埋深分別為2 km和11 km，厚度均為4 km.

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

該模型測試中，我們給出了三種不同深度加權參數的反演結果，分別對應了三種不同的深度加權參數，如表2中α、β參數值所示，表2中不同的深度加權參數值分別對應了圖3b、c、d的反演結果.當表2中加權參數α值過大時，得到的反演結果圖3b中高密度部分主要集中在底部，異常體的深度與真實模型埋深存在一定的差異；而當加權參數β過大時，得到的圖3c的反演結果其高密度主要集中在淺部，異常體埋深同樣存在較大差異.通過本文算法得到的最優化深度加權參數反演結果如圖3d所示，異常體的深度信息較為準確，得到的反演結果埋深和形態與真實模型基本一致.由此可得，與主觀指定的深度加權參數反演結果相比，本文提出的貝葉斯同化反演方法，在已知局部參考信息下，給出的深度加權更合理，反演的異常體深度與真實模型更加一致.

圖3 深度加權參數優化反演算法測試(a) 真實模型； (b) 過衰減加權系數反演結果； (c) 欠衰減加權系數反演結果； (d) 優化加權系數反演結果.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 不同深度超參數約束下反演統計及ABIC參數特征Table 2 Inversion statistics and ABIC parameter characteristics under different depth hyperparameter constraints

綜合以上兩個模型的測試結果，可以認為本文算法具備同化多個參考模型和優化深度加權參數的能力.下面我們進一步應用該方法，選擇地球物理觀測資料豐富、地殼結構參考模型較多的龍門山地區進行實際數據測試.

3 龍門山地殼模型

3.1 龍門山構造背景

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

圖4 龍門山地區主要構造特征與地震活動紅色實線為斷裂，四條黑色實線為跨震中研究剖面，黃色圈為5級以上地震位置，白色實心圓為中國地震觀測網臺站位置，F1：鮮水河斷裂帶；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).

同時，該地區也屬于中國地震科學實驗場區，過去幾十年中以中國地震科學臺陣項目為代表的大規模地球物理觀測相繼開展，已積累了大量的地球物理探測工作成果(如Yao et al., 2008; Li et al., 2011; Zhang et al., 2011; 王緒本等，2018).這些深部成果對該區的殼幔結構及其相關動力學問題達成了部分共識，但在一些基本問題上依舊存在爭議，如在該區的地殼上地幔變形機制的解釋上，就有殼幔連續變形機制(England and Molnar, 1997)、塊體擠出機制 (Tapponnier et al., 1982,2001)、下地殼流機制(Clark and Royden, 2000； Royden et al., 1997)等多種模式.

此外，龍門山地區地震多發，如2008年汶川8.0級和2013年蘆山7.0級強震.雖然國內外研究學者對該區大震震源區的深部孕震結構等進行了大量的研究，但對于地震孕育深部地殼結構特征及相關動力學機制還存在分歧.如：Zhang等(2011)通過地震層析結果，認為龍門山斷裂帶兩側高低速異常的邊界是該區強震的凹凸區；而房立華等(2013)和王夫運等(2015)的地震測深剖面顯示斷裂帶下方中滑脫層是產生研究區地震的主要原因；張培震等(2008)、楊文采等(2015)和王緒本等(2018)的多種研究結果表明研究區強震的發生與地殼內部存在的低速高導層有關.

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

3.2 研究區已有地殼結構

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

表3 研究區內已有部分地震波速成果Table 3 Partial seismic wave velocities in the study area available

如果將已有的模型解釋結果看作先驗認識，那么如何根據這些先驗去改進反演模型的估計？這是本文提出的重力異常貝葉斯同化反演新算法應該回答的問題.我們將以多種不同類型的地震速度轉換密度作為已知先驗參考，利用高精度的衛星重力場模型數據，通過新算法實現對不同先驗參考約束的取長補短，剔除偏差數據在反演計算中的作用，以期得到同時符合計算系統不同先驗假設的最優高精度地殼結構模型.同時測試算法在實際重力資料反演中的效果，為研究該區的地殼變形機制、強震孕育環境及相關動力學提供有益參考.

3.3 參考模型建立

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

圖5 地震層析成像和接收函數轉換密度結果(30 km水平切片)實心圓填充的顏色代表接收函數轉換密度，底圖是層析成像轉換密度.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 重力異常特征

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

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

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

4 反演結果與分析

根據上一節構建的兩個參考密度模型和布格重力異常，我們完成了貝葉斯同化重力反演.下面我們分別對從結果的可靠性和剖面特征兩方面進行論述.

4.1 可靠性分析

(1)模型殘差

通過前文所述的反演方法，我們通過貝葉斯優化，獲得了ABIC最小值對應的反演結果，模型中四個超參數分別為：λ1=139.254、λ2=122.965、α=31.337、β=0.45.圖7中給出了反演模型的正演異常值和其與觀測值的差異特征，圖7b的異常殘差結果標準差為±2.5 mGal.

從圖7a中可以看出，最終密度異常正演所得異常理論值與觀測剩余值形態具有較好的一致性，與地震模型正演結果相比有更多短波長特征.圖7b中觀測異常和反演密度模型正演理論異常得到數據擬合剩余殘差主要是高頻誤差，標準差略優于WGM20121布格異常的5 mGal精度.

(2)水平結構

為了進一步檢驗反演結果的可靠性，我們仿照圖5的結果，取30 km深度切片，對比接收函數轉換密度結果與反演結果之間的差異，如圖8所示.可以看出：以龍門山斷裂帶為界(F2)，川西高原與四川盆地高低密度分界清晰.與地震層析結果相比，聯合反演所得密度分布形態與接收函數具有較好的一致性；此外，川西高原顯著低密度異常的背景下出現了星狀分布的小尺度橫向密度差異特征，四川盆地內部密度相對川西呈高密度結構特點，也伴隨小尺度相間分布的結構特征.

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

圖7 (a) 反演模型結果正演重力異常； (b) 觀測異常和反演模型正演異常的差異特征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重力異常貝葉斯同化反演密度結果(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 反演模型與參考模型點垂向對比(a) P1點；(b) P2點；(c) P3點；(d) P4點；(e) P5點；(f) P6點.藍色折線為接收函數轉化密度結果；紅色折線為層析成像轉化密度結果；黑色空心圓為反演密度結果，黑色短線為密度估計誤差.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)垂直結構

在垂向分層結構的檢測方面，我們在不同構造區，分別挑選了6個接收函數臺站位置下方的一維垂向密度結構進行比較.圖9a—f給出了接收函數參考模型(藍色實線)、地震層析參考模型(紅色實線)和同化反演模型(黑點)的一維結果.總體來說，三者結果差異不大，最終反演結果在不同深度位置綜合了兩種參考模型信息，可說明重力同化反演結果的垂向分辨能力可以從參考模型中獲得.

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

4.2 剖面特征

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

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

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

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

5 結論與討論

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

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

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

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

綜上所述，本文提出的貝葉斯反演算法，反演過程不依賴人的主觀認識，完全依靠數據說話，通過非線性優化算法可以實現全自動化地反演參數調節和模型同化，提取不同先驗約束之間的有效信息，可以充分發揮地震學方法的垂向分辨能力優勢和重力異常蘊含的水平密度結構變化特征.

在對強震震源區的結構研究方面，本文得到的密度剖面結果顯示，蘆山和汶川震中區位置都存在顯著的中地殼密度界面陡變特征，其斷層位置與密度分層界面陡變區都具有較好的一致性.同時龍門山斷裂帶下方滑脫層的存在與該斷裂帶下方密度分層界面也有較好的對應(房立華等, 2013; 王夫運等, 2015)，因此，此類高精度的地殼密度結構無疑更有益于研究大地震孕育環境的結構與物性特征.

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

致謝感謝中國科學臺陣項目提供的地震數據，感謝中國地震局地球物理研究所王興臣博士提供的川滇地殼接收函數結果，感謝中國地震局地球物理研究所李永華研究員和兩位匿名審稿人提供的寶貴建議和幫助.