李永波， 陳石*， 李紅蕾， 張貝

1 中國地震局地球物理研究所， 北京 100081 2 北京白家疃地球科學(xué)國家野外觀測科學(xué)研究站， 北京 100095

0 引言

重力異常能很好地反映地殼內(nèi)部介質(zhì)的密度變化，因此常被用于礦產(chǎn)勘探(Jordens et al., 2014)、地形起伏(Silva et al., 2014;Pallero et al., 2015)以及地殼結(jié)構(gòu)(van der Meijde et al., 2013;Steffen et al., 2017;Xu et al., 2017)等多個與地球內(nèi)部密度分布相關(guān)的領(lǐng)域.利用重力場數(shù)據(jù)提取密度分布信息或特殊密度體結(jié)構(gòu)的過程稱為重力反演.通過重力反演手段獲得的地殼三維密度結(jié)構(gòu)，可用于研究地殼深部孕震區(qū)物性結(jié)構(gòu)特征(陳石等, 2014;李紅蕾等, 2021b).由于重力反演具有較好的橫向分辨能力，適合與接收函數(shù)、地震層析成像等具備垂向分辨能力的手段開展聯(lián)合反演(Lelièvre et al., 2012;張佩等, 2019)，以提高研究目標(biāo)物性異常體整體的分辨能力.

本文研究區(qū)域為四川省南部和云南省北部交界地區(qū)(103°E—106°E,27°N—30°N)，在大地構(gòu)造上位于青藏高原東南緣.如圖1所示，研究區(qū)主要包括了四川盆地南部褶皺帶邊緣、大涼山褶皺帶邊緣以及四川盆地西南區(qū)域，研究區(qū)內(nèi)有華鎣山斷裂、五蓮峰斷裂以及龍泉山斷裂等活動構(gòu)造.根據(jù)1970年以來的地震目錄記載，研究區(qū)域內(nèi)5級以上中強(qiáng)地震主要集中于四川盆地西南緣的盆山交界地帶.在馬邊地區(qū)，1936年4月26日和5月16日相隔20多天分別發(fā)生MW6.7和MW6.8兩次地震，1971年8月17日前后連續(xù)發(fā)生一系列5級左右地震，最大達(dá)5.8級；在大關(guān)地區(qū)，最大地震是1974年5月11日發(fā)生的MS7.1地震，此次地震造成了巨大的人員傷亡和財產(chǎn)損失.但近年來榮縣、自貢、涪陵等地由于頁巖氣開發(fā)、深井采鹽等工業(yè)化活動，導(dǎo)致盆地內(nèi)部地震活動增加.特別是2019年6月17日長寧MS6.0地震、2021年9月16日瀘縣MS6.0地震都發(fā)生在盆地內(nèi)部，雖然這兩次地震震中位置周邊都存在斷裂構(gòu)造，但近幾十年來并未發(fā)生超過6.0級的地震.因此，這兩次6.0級地震到底是人類活動誘發(fā)的地震，還是構(gòu)造地震？逐漸成為了地球物理工作者非常關(guān)注的熱點研究問題(Sun et al., 2017;He et al., 2019a,b;Lei et al., 2019a,b;Li et al., 2021).針對該地區(qū)已有大量的研究成果，讓我們對這一區(qū)域有了一定程度的認(rèn)識，但是這些研究主要集中于地震學(xué)方法(Wang et al., 2016;Long et al., 2020;Yang et al., 2020;Zuo et al., 2020;李大虎等, 2021;孫權(quán)等, 2021;張杰等, 2020).而本研究嘗試?yán)弥亓?shù)據(jù)反演該區(qū)域的地殼密度結(jié)構(gòu)，進(jìn)一步分析該地區(qū)多個地震活動區(qū)與物性結(jié)構(gòu)之間的關(guān)系.

2019年6月17日長寧MS6.0地震后，中國地震局組織了實地科考，并在川南地區(qū)完成了多條重力探測剖面，獲得了陸地實測重力數(shù)據(jù)，經(jīng)過異常改正之后得到該區(qū)域的布格重力異常.但是由于陸地重力測網(wǎng)形態(tài)不規(guī)則，單純利用實測重力數(shù)據(jù)難以填補(bǔ)測量空區(qū)導(dǎo)致的數(shù)據(jù)缺失問題.為此，本文在反演過程中結(jié)合陸地實測重力異常和WGM2012布格重力異常模型，通過將地表實測的布格重力異常與WGM2012布格重力異常模型采用等效源方法進(jìn)行融合，有效剔除了WGM2012布格重力異常模型的系統(tǒng)誤差，壓制了高頻干擾，同時還提高了該區(qū)域布格重力異常的精度和分辨率(李紅蕾等, 2021a)，獲得了可靠性更高的布格重力異常數(shù)據(jù).

在反演方法方面，本文結(jié)合該區(qū)已有的地震學(xué)觀測和模型結(jié)果，采用貝葉斯原理進(jìn)行重力約束反演，利用馬爾可夫鏈蒙特卡羅(MCMC)算法進(jìn)行最優(yōu)化求解.在研究過程中，首先通過模型測試驗證了該方法的有效性，然后將其運用到實際資料處理中，獲得了川南地區(qū)地殼的三維密度結(jié)構(gòu)，同時還對結(jié)果的不確定性給出了估計.

1 反演方法和模型測試

傳統(tǒng)的重力反演為了解決多解性問題，一般采用的是正則化方法，使用最為廣泛的是Tikhonov正則化(Calvetti and Somersalo, 2018)，其優(yōu)點是可以在不改變似然函數(shù)的條件下通過在目標(biāo)函數(shù)中額外增加一個約束項來對參數(shù)進(jìn)行懲罰，從而降低反演多解性.常用的約束項包括：最小模型約束，空間光滑約束以及深度加權(quán)約束等(Li and Oldenburg, 1998;Boulanger and Chouteau, 2001).通過引入正則化約束可以得到最優(yōu)解，但并不意味著這就是唯一可能的解(Muoz and Rath, 2006)，同時也很難判斷這個結(jié)果是否與真實的密度結(jié)構(gòu)相符.如何在獲得最優(yōu)解的同時考慮到解的其他可能性，得到對模型參數(shù)更加準(zhǔn)確而全面的估計是我們在反演過程中重點關(guān)心的問題.

圖1 研究區(qū)地形與MS5.0以上地震分布 黑色線為反演結(jié)果重點討論的八條剖面所在位置；紅色線為斷層分布.Fig.1 Topography and distribution of earthquakes greater than MS5.0 in the study area The black lines are the positions at which cross-sections of inverted density model will be discussed later; The red lines indicate the faults.

近年來，貝葉斯統(tǒng)計學(xué)模型在多個地球物理領(lǐng)域得到了廣泛應(yīng)用，可用于評估震源模型(Minson et al., 2013)、地球物理反演(Ray et al., 2018; Blatter et al., 2019;Figueiredo et al., 2019;Gebraad et al., 2020)、重力數(shù)據(jù)平差(Chen et al., 2019)等.在重力反演方面，Hightower等(2020)利用貝葉斯模型反演復(fù)雜地區(qū)密度分布，李紅蕾等(2021b)利用經(jīng)驗貝葉斯方法反演龍門山地區(qū)地殼三維密度結(jié)構(gòu).相比于傳統(tǒng)的確定性反演方法，貝葉斯反演方法可以更好地將先驗認(rèn)識引入模型，同時還可以將觀測數(shù)據(jù)的不確定性以及先驗參數(shù)的不確定性傳遞給后驗分布，獲得求解參數(shù)可能的分布范圍及其在每一區(qū)間分布的概率，能更加完整的反映所求參數(shù)的性質(zhì).獲得后驗分布樣本后，我們可以利用后驗均值、后驗中位數(shù)或者最大后驗來提取最適合的參數(shù)估計值.

在貝葉斯反演框架中，反演問題如果想要獲得具有實際意義的解，需要設(shè)置合理的先驗分布.貝葉斯反演方法可以靈活的引入先驗性息，同時在以觀測數(shù)據(jù)控制的似然函數(shù)作用下獲取新的認(rèn)識(Blatter et al., 2019).若先驗信息已經(jīng)較好地反映了真實的參數(shù)信息，則似然函數(shù)不會對后驗分布產(chǎn)生太大影響.若先驗信息對于參數(shù)的認(rèn)識是不準(zhǔn)確的，那么在新的觀測數(shù)據(jù)的作用下會使得后驗分布在先驗的基礎(chǔ)上向正確的方向偏移(Gallagher et al., 2009).

1.1 基本原理

在三維密度反演中，通常的做法是將反演區(qū)域劃分為正六面體單元，并給每一個單元體一個密度值，在反演過程中，單元體的位置和大小不變，反演單元體的密度值(Li and Oldenburg, 1998;Boulanger and Chouteau, 2001).用m表示場源體的密度模型，如果將反演目標(biāo)區(qū)域劃分成N個六面體單元，那么模型參數(shù)可以表述為m=[m1,m2,m3,…,mN]T.同時用d表示布格重力異常矢量，假設(shè)有P個觀測數(shù)據(jù)，那么重力異?？梢员硎緸閐=[d1,d2,d3,…,dP]T.

布格重力異常通常與密度異常分布有關(guān)，在密度體空間位置確定的條件下，重力異常與密度值之間的關(guān)系可用線性泛函表達(dá)，其數(shù)學(xué)關(guān)系可以表示為

d=Gm,

(1)

其中，G為核函數(shù)矩陣；m為場源密度模型.

依據(jù)貝葉斯理論，可以將三維密度反演問題表示為

(2)

式中，m為所求參數(shù)；d為觀測數(shù)據(jù)；P(d|m)為似然函數(shù)，表示數(shù)據(jù)擬合程度；π(m)表示先驗認(rèn)識，獨立于觀測數(shù)據(jù)；P(d)為歸一化因子；P(m|d)表示后驗分布.在(2)式中，分母項表示的歸一化因子為常數(shù)，后驗條件概率分布可以簡化為：

P(m|d)∝P(d|m)π(m).

(3)

似然函數(shù)反映了反演模型參數(shù)與觀測數(shù)據(jù)的擬合程度，在貝葉斯理論框架下，似然函數(shù)可以表示為如下形式：

式中，Cd為協(xié)方差矩陣.

在貝葉斯反演框架中，先驗認(rèn)識的選取會影響最終的反演結(jié)果.本文將研究區(qū)內(nèi)的多個已有觀測的地震接收函數(shù)一維速度結(jié)構(gòu)轉(zhuǎn)化為密度結(jié)構(gòu)作為參考模型，依此參考模型設(shè)定參數(shù)的先驗分布.其基本思想是：參考模型能反映參數(shù)分布的基本特征，參數(shù)先驗圍繞參考模型在一定范圍內(nèi)進(jìn)行調(diào)整，且其差值服從正態(tài)分布.因此，參數(shù)的先驗分布可以表達(dá)如下：

(5)

式中，mref為參考密度，由速度結(jié)構(gòu)轉(zhuǎn)換得到；Cm為模型協(xié)方差矩陣，反映我們對已有認(rèn)識的信賴程度.在此基礎(chǔ)上，參數(shù)的后驗分布可以表示為

(6)

我們可以看到，在(6)式右端項括號中的部分與傳統(tǒng)反演方法的目標(biāo)函數(shù)是類似的，因此，在貝葉斯框架下也可以靈活的引入正則約束.

1.2 MCMC最優(yōu)化

本文在反演中，與李紅蕾等(2021b)采用的經(jīng)驗貝葉斯方法求解模型參數(shù)所不同，而是選擇了完全貝葉斯反演方法，并基于MCMC求解方法，獲得參數(shù)的后驗概率分布.MCMC方法通過直接對后驗參數(shù)進(jìn)行采樣獲取后驗樣本，當(dāng)樣本數(shù)量足夠大時便可以用樣本統(tǒng)計特征值近似代替參數(shù)的實際分布特征，如后驗均值、后驗眾數(shù)以及后驗標(biāo)準(zhǔn)差等.在本文研究中，我們將后驗分布的均值作為參數(shù)的點估計值，同時利用參數(shù)的標(biāo)準(zhǔn)差來衡量參數(shù)自身的不確定性大小.

常見的MCMC采樣方法包括Metropolis-Hastings抽樣(Metropolis et al., 1953;Hitchcock, 2003; Ardilly and Tillé, 2006)、Gibbs抽樣(Metropolis et al., 1953; Ardilly and Tillé, 2006)、Hamiltonian抽樣(Hanson, 2001)以及No-U-Turn抽樣(Neal, 2011; Homan and Gelman, 2014)等.本文采用NUTS(No-U-Turn Sampling)算法進(jìn)行采樣(Salvatier et al., 2016)，NUTS是對Hamiltonian Monte Carlo方法的進(jìn)一步改進(jìn)，相對于Metropolis-Hastings與Gibbs等隨機(jī)游走型抽樣方法而言，NUTS方法在處理多維問題時更容易達(dá)到收斂狀態(tài)(Homan and Gelman, 2014).

1.3 模型測試

重力反演由于其多解性特征，在獲取反演結(jié)果的時通常需要參考模型作為先驗認(rèn)識，反演結(jié)果的準(zhǔn)確性很大程度上取決于先驗信息的可靠程度，當(dāng)先驗信息不足以給出有用的信息時，反演結(jié)果往往也能非常好的擬合觀測數(shù)據(jù)，但很難反映真實的參數(shù)特征.貝葉斯反演方法基于參考模型產(chǎn)生參數(shù)的先驗分布，通過調(diào)整先驗分布的特征以控制參考模型對反演結(jié)果的影響程度.

首先，構(gòu)建了一個維度為16×4×8(nx×ny×nz)的離散化網(wǎng)格模型，該模型中有兩個不同剩余密度和幾何形態(tài)的異常體，如圖2a所示.為了便于分析，模型在y方向上密度不變，僅在x方向和z方向上設(shè)置了密度異常體.在模型正演結(jié)果的基礎(chǔ)上加入5%的高斯噪聲獲得的重力異常如圖2b所示.然后，在貝葉斯方法的反演效果測試中，通過設(shè)置參數(shù)先驗分布的標(biāo)準(zhǔn)差來控制反演結(jié)果對于參考模型的依賴程度，并分析參數(shù)不確定性與參數(shù)空間分布特征之間的關(guān)系.

在構(gòu)建參考模型時，將真實模型作為基準(zhǔn)模型，并針對不同位置給與不同的偏差，其背景密度與真實模型相同，梯形部分密度比真實模型小0.2 g·cm-3，淺部高密度層比真實模型大0.2 g·cm-3.假設(shè)先驗分布是以參考模型數(shù)值大小作為均值的正態(tài)分布，通過設(shè)置該正態(tài)分布的標(biāo)準(zhǔn)差來反映參考模型的可信賴程度，標(biāo)準(zhǔn)差越小，說明參考模型的可信賴程度越高.

如圖3所示，當(dāng)先驗分布標(biāo)準(zhǔn)差設(shè)置為0.02 g·cm-3時，淺部參數(shù)恢復(fù)效果較好，隨著深度增加，后驗參數(shù)受參考模型的影響減小(圖3a)，殘差隨著深度增加而變大(圖3b)，參數(shù)恢復(fù)效果減弱.當(dāng)先驗分布標(biāo)準(zhǔn)差設(shè)置為0.01 g·cm-3時，參數(shù)的恢復(fù)效果整體上具有明顯的提升(圖3c), 密度殘差減小(圖3d).

上述實驗結(jié)果顯示，可通過設(shè)置先驗分布標(biāo)準(zhǔn)差的大小來反映對參考模型的信賴程度，進(jìn)而控制參考模型對反演結(jié)果的影響.

基于MCMC方法的完全貝葉斯反演，因為要通過大量的采樣來生成參數(shù)的后驗分布樣本，因此可以給出反演結(jié)果的不確定性.圖4為上述實驗結(jié)果的不確定性特征，圖4a為先驗分布標(biāo)準(zhǔn)差為0.02 g·cm-3時的不確定性特征，圖4b為先驗分布標(biāo)準(zhǔn)差為0.01 g·cm-3時的不確定性特征.從這兩幅圖中，我們可以總結(jié)出關(guān)于參數(shù)不確定性特征分布的三個特點：(1)反演結(jié)果的不確定度隨著深度的增加而增加，這與重力響應(yīng)的敏感性隨著密度體埋深增加而降低這一特征是相符的；(2)先驗分布特征同樣會影響后驗結(jié)果的不確定性，先驗信息分布范圍越小，反演結(jié)果的不確定性也就越小;(3)參考模型本身的準(zhǔn)確度會影響反演結(jié)果的不確定性，我們可以看到，在同一深度處，由于參考模型梯形部分本身的密度與真實密度存在差異，導(dǎo)致該處反演結(jié)果的不確定性相比于背景處明顯增大.

基于貝葉斯原理的重力約束反演在先驗分布合理時可以提高反演結(jié)果的準(zhǔn)確性，同時對參數(shù)的不確定性做出評價，得到對所求參數(shù)更加全面的估計.因此，此方法對于先驗認(rèn)識較為充足的研究區(qū)具有很好的適用性.

圖2 密度模型及重力異常 (a)真實模型;(b)模型正演重力異常，其中每一個數(shù)據(jù)點添加了其自身幅值5%大小的高斯噪聲.Fig.2 Density model and gravity anomaly(a) True model; (b) The gravity anomaly produced by density model, each datum has been contaminated by uncorrelated Gaussian noise of 5% of the datum magnitude.

圖3 密度模型反演結(jié)果及殘差分布 (a)先驗分布標(biāo)準(zhǔn)差為0.02 g·cm-3時的反演結(jié)果；(b) 先驗分布標(biāo)準(zhǔn)差為0.02 g·cm-3時的模型殘差分布； (c) 先驗分布標(biāo)準(zhǔn)差為0.01 g·cm-3時的反演結(jié)果；(d) 先驗分布標(biāo)準(zhǔn)差為0.01 g·cm-3時的模型殘差分布.Fig.3 Inversion results and residual distribution(a) Inversion results when the prior model standard deviation is 0.02 g·cm-3; (b) The residual distribution when the prior model standard deviation is 0.02 g·cm-3; (c) Inversion results when the prior model standard deviation is 0.01 g·cm-3; (d) The residual distribution when the prior model standard deviation is 0.01 g·cm-3.

圖4 密度模型反演結(jié)果的不確定性特征分布(a) 先驗分布標(biāo)準(zhǔn)差為0.02 g·cm-3時反演結(jié)果的不確定性特征； (b) 先驗分布標(biāo)準(zhǔn)差為0.01 g·cm-3時反演結(jié)果的不確定性特征.Fig.4 Uncertainty of inversion results(a) Uncertainty characteristics of inversion results when the standard deviation of the prior model is 0.02 g·cm-3; (b) Uncertainty characteristics of inversion results when the standard deviation of the prior model is 0.01 g·cm-3.

2 重力數(shù)據(jù)和參考模型

2.1 布格重力異常

WGM2012布格重力異常模型是基于EGM2008重力場模型和DTU10重力場模型提供的自由空氣異常數(shù)據(jù)為基礎(chǔ)，采用ETOPO1高分辨率地形模型進(jìn)行重力校正得到的全球布格重力異常模型，該布格重力異常模型具有分辨率高、覆蓋范圍廣的特征(Balmino et al., 2012;Pavlis et al., 2012)，但是全球重力場模型數(shù)據(jù)源以衛(wèi)星資料為主，在中國陸地區(qū)域缺少地表觀測資料的約束，這導(dǎo)致WGM2012布格重力異常模型與陸地實測布格重力異常之間存在較大的系統(tǒng)誤差(李紅蕾等, 2021a).

本文在研究過程中，將地表實測的9條重力剖面數(shù)據(jù)(位置如圖5所示)與WGM2012布格重力異常模型采用等效源方法進(jìn)行融合(李紅蕾等, 2021a)，WGM2012布格重力異常分辨率為2′×2′，而陸地相鄰實測重力點之間的距離約2～5 km.在數(shù)據(jù)融合過程中，充分結(jié)合兩種數(shù)據(jù)的優(yōu)勢，得到的重力異常在陸地測網(wǎng)分布區(qū)分辨率可達(dá)7 km.從圖5中可以看出，布格異常從東北到西南方向呈逐漸降低之勢，與地形起伏具有明顯的負(fù)相關(guān)特征，較好地反映出該區(qū)的盆山結(jié)構(gòu)具有重力均衡特點.同時，該融合布格重力異常模型包含了較豐富的深部地殼構(gòu)造信息.

2.2 參考模型

重力數(shù)據(jù)對橫向密度變化有很好的敏感性，但是缺乏垂向分辨能力，在實際反演中通常需要加入垂向密度變化信息作為先驗認(rèn)識以提高反演結(jié)果的垂向分辨能力.而地震接收函數(shù)則可以提供觀測臺站下方的一維垂向速度結(jié)構(gòu)模型，且對速度變化界面十分敏感.因此，在本文中我們首先將收集到的區(qū)域內(nèi)地震一維接收函數(shù)結(jié)果(Wang et al., 2017;鄭晨等, 2016)，采用地震波速度-密度經(jīng)驗公式(Brocher, 2005)將速度模型轉(zhuǎn)換為密度模型作為此次反演的參考模型.

圖6展示了參考模型在10 km深度處的密度分布，其中藍(lán)色倒三角為研究區(qū)內(nèi)具有地震接受函數(shù)結(jié)果的地震臺站分布，共計182個觀測臺站.為最大程度的保留接收函數(shù)結(jié)果攜帶的垂向變化信息，在構(gòu)建參考模型時垂向分層直接依據(jù)接收函數(shù)結(jié)果劃分，本文選取了從地表向下沿伸的前30層數(shù)據(jù)，其最大深度為49.2 km，在水平向上按照0.2°×0.2°的網(wǎng)格大小進(jìn)行插值.從圖6中可以看出在本文研究的幾次地震震中區(qū)及周邊都存在地震臺站分布，因此以該參考模型提供的垂向密度分層信息為約束獲得的反演結(jié)果，對研究孕震區(qū)構(gòu)造特征具有較好的可信度.

圖5 研究區(qū)融合布格重力異常 黑色圓圈為陸地重力測量點(共9條測線).Fig.5 Fusion Bouguer gravity anomaly in the study area Dark circles indicate the location of terrestrial gravity measurement points (9 survey lines in total).

圖6 參考模型10 km深度處水平切片及地震觀測臺站位置(藍(lán)色倒三角)Fig.6 Horizontal slice at 10km depth of reference model and the location of seismic observation stations (Blue inverted triangles)

圖7 研究區(qū)剖面位置對應(yīng)的重力異常殘差特征Fig.7 Residual characteristics of gravity anomaly corresponding to profile positions in the study area

圖8 重力異常殘差分布 圖中右下方為殘差統(tǒng)計結(jié)果.Fig.8 Residual gravity between the observed and predicted gravity grids The statistical results of residuals are shown at the bottom right of theFigure.

圖9 研究區(qū)三維密度反演結(jié)果與參考模型差異統(tǒng)計直方圖(a) 所有參數(shù)差異統(tǒng)計直方圖； (b) 10 km以內(nèi)參數(shù)差異統(tǒng)計直方圖； (c) 10 km到29.5 km深度內(nèi)參數(shù)差異統(tǒng)計直方圖； (d) 29.5 km到49.2 km深度內(nèi)參數(shù)差異統(tǒng)計直方圖.Fig.9 Statistical Characteristics of the difference between the inversion results and the reference model(a) Statistical histogram of all parameter differences; (b) Statistical histogram of parameter difference within 10 km; (c) Statistical histogram of parameter differences within the depth of 10 km to 29.5 km; (d) Statistical histogram of parameter differences at depths from 29.5 km to 49.2 km.

圖10 研究區(qū)不同深度處密度結(jié)構(gòu)水平切片F(xiàn)ig.10 Horizontal slices of density structure at different depths in the study area

圖11 AA′-HH′八個垂直剖面密度分布特征(剖面位置如圖1所示)Fig.11 Density distribution characteristics of eightcross-sections AA′-HH′ (cross-section positions are shown in Fig.1)

圖12 反演結(jié)果不確定性估計Fig.12 Uncertainty estimation of inversion results

3 反演結(jié)果

參數(shù)的先驗分布是以參考模型為均值的正態(tài)分布.本研究在實際資料反演過程中發(fā)現(xiàn)，當(dāng)先驗分布標(biāo)準(zhǔn)差設(shè)置過大時，反演結(jié)果預(yù)測重力異?？梢苑浅：玫臄M合觀測數(shù)據(jù)，但是反演結(jié)果在空間上的變化特征相比于參考模型容易產(chǎn)生較大的偏離，同時導(dǎo)致反演結(jié)果的分辨率降低.當(dāng)先驗分布標(biāo)準(zhǔn)差設(shè)置過小時，反演得到的結(jié)果與參考模型高度一致，但是存在數(shù)據(jù)擬合效果較差的問題.

通過對實際數(shù)據(jù)多次計算發(fā)現(xiàn)，當(dāng)參數(shù)先驗分布標(biāo)準(zhǔn)差設(shè)置為0.2 g·cm-3時，反演結(jié)果在其預(yù)測重力異常擬合觀測數(shù)據(jù)的前提下，也能基本保留參考模型攜帶的密度變化特征.因此，本文最終采用的是先驗分布標(biāo)準(zhǔn)差為0.2 g·cm-3時的反演結(jié)果.

3.1 對比分析

(1)重力異常對比

圖7給出了在圖4中所示的過四次強(qiáng)震震中位置的8條剖面對應(yīng)的重力異常數(shù)據(jù)殘差特征.其中，綠色線為參考模型正演獲得的重力異常與融合重力異常之間的差異，可以看出參考模型正演的重力異常與融合布格異常之間存在較明顯的差異，其最大值可達(dá)60 mGal；而紅色曲線為反演結(jié)果正演重力異常與融合重力異常之間的殘差.在這8條剖面的對比結(jié)果中，反演結(jié)果正演重力異常整體上能較好地擬合該區(qū)的布格重力異常場，但在剖面A-A′與E-E′兩條曲線的大關(guān)位置附近出現(xiàn)了較大的殘差信號，通過與該區(qū)位置的布格重力異常數(shù)據(jù)對比，我們發(fā)現(xiàn)這些點誤差較大的原因是該區(qū)域存在變化劇烈的高頻異常信號.

為了進(jìn)一步分析反演結(jié)果的殘差特征，圖8給出了整個研究區(qū)范圍內(nèi)的重力反演模型擬合殘差空間分布及殘差統(tǒng)計直方圖，結(jié)果顯示重力異常擬合殘差服從正態(tài)分布，且該正態(tài)分布的均值為0.0026 mGal，標(biāo)準(zhǔn)差為2.5495 mGal.對比圖5中的異常數(shù)值變化范圍，該反演結(jié)果給出的三維密度結(jié)構(gòu)模型可以提供與重力觀測異常較為一致的解釋.

(2)反演結(jié)果與參考模型差異對比

圖9展示了密度參數(shù)反演結(jié)果與參考模型之間差異的統(tǒng)計特征，從圖9a可以看出反演結(jié)果給出的密度結(jié)構(gòu)與參考模型之間的差異主要集中于-0.4 g·cm-3到0.3 g·cm-3之間，參數(shù)差異的均值為-0.0150 g·cm-3.說明參考模型相對于反演結(jié)果的密度差均值為負(fù)，即參考模型整體上略偏大.圖9(b—d)分別展現(xiàn)了參考模型差異在不同深度上的分布特征，通過對比我們發(fā)現(xiàn)淺層參數(shù)差異與深層參數(shù)差異具有不同的特性，對于10 km以上的淺層，反演得到的密度參數(shù)整體上略大于參考模型，而在中深層反演結(jié)果整體上略小于參考模型，特別是圖9c中在-0.3 g·cm-3左右出現(xiàn)了一個小高峰，說明在10～29.5 km深度范圍內(nèi)，參考模型在一些位置上的密度值與本文反演給出的密度結(jié)果之間存在較為顯著的差異.

圖6和圖10a分別為參考模型和反演結(jié)果在10 km深度的密度結(jié)構(gòu)水平切片，通過對比可以看出，本文反演結(jié)果雖然在整體結(jié)構(gòu)上與參考模型基本一致，但是仍然存在很多細(xì)節(jié)差異.由于布格重力異常場的橫向分辨率優(yōu)勢，通過基于貝葉斯理論的重力約束反演可以獲得更多重力數(shù)據(jù)所攜帶的物性橫向變化信息，進(jìn)一步提高物性模型的分辨率.

(3)水平密度結(jié)構(gòu)特征

本文分別給出了反演結(jié)果在10 km、20 km、30 km和40 km四個不同深度處的密度結(jié)構(gòu)水平切片(圖10).結(jié)果顯示在10 km深度處的淺部地殼，密度從西南到東北方向大體上呈減小的趨勢，山地地區(qū)密度整體上高于盆地內(nèi)部，尤其在四川盆地西南周緣的盆山結(jié)合帶部位密度變化高梯度帶明顯，說明這些地方的淺部地殼構(gòu)造變形復(fù)雜.隨著深度增加，密度變化趨勢發(fā)生反轉(zhuǎn)，如40 km深度剖面所展示的這樣，密度從西南方向的最小值一直增加到東北方向的最大值，與地形變化趨勢正好相反，這與Airy重力均衡假說是基本符合的.

(4)垂直密度結(jié)構(gòu)特征

為了更多好地對比幾次中強(qiáng)震震源區(qū)密度結(jié)構(gòu)特征，本文在大關(guān)、馬邊、長寧以及瀘縣四個震源區(qū)分別繪制了沿EW向和SN向的密度剖面，如圖11所示.圖11a和e所示的A-A′剖面與E-E′剖面展示了大關(guān)地震震源區(qū)下方的密度分布特征，通過A-A′剖面我們看到大關(guān)地震震源下方密度沿東西向具有較為明顯的橫向變化，推測可能與1974年的大關(guān)MS7.1地震的構(gòu)造背景有關(guān).圖11b和f所示的B-B′與F-F′剖面展示了長寧震源區(qū)下方的密度結(jié)構(gòu)，B-B′剖面下方密度分布平滑，沒有明顯的橫向不均勻性，而在F-F′剖面上，我們可以看到震源區(qū)下方正好有一個條形的高密度體，說明長寧地震也有可能受構(gòu)造活動影響.

圖11c和g所示的C-C′剖面與G-G′剖面反映了馬邊地震下方密度分布特征，從圖中可以看到馬邊地震序列正好位于高密度突出體的尖部，即馬邊地震發(fā)生位置與構(gòu)造特征具有很強(qiáng)的一致性.圖11d和h所示的D-D′與H-H′兩條剖面可以發(fā)現(xiàn)瀘縣MS6.0地震下方在東西方向與南北方向密度分布均未體現(xiàn)出明顯的橫向不均勻性，這說明2021年9月16日的瀘縣MS6.0地震的深部地殼構(gòu)造變形特征不明顯.

3.2 不確定性分析

利用貝葉斯理論進(jìn)行反演的優(yōu)勢在于可以對參數(shù)進(jìn)行不確定性分析，圖12展示了AA′剖面的不確定性特征，從圖中可以看出參數(shù)不確定性最大可達(dá)0.011 g·cm-3.參數(shù)不確定性整體上呈現(xiàn)出隨深度增加的特性，這與模型測試結(jié)果是符合的，說明基于貝葉斯原理的重力反演方法得到的密度參數(shù)后驗點估計值的精度隨著深度增加有所降低.同時，我們發(fā)現(xiàn)剖面兩端參數(shù)不確定性相比中間部分略大，這一特性可能與邊界部分相關(guān)聯(lián)的重力數(shù)據(jù)觀測量少于中間區(qū)域有關(guān)，即有效觀測數(shù)據(jù)的豐富程度也能影響反演結(jié)果的可靠性.

我們還注意到頂部出現(xiàn)不確定性較大的異常特征，這主要是受層厚影響，在反演模型構(gòu)建時，淺層厚度小，導(dǎo)致相同密度變化引起的重力異常小，因此其不確定性較大.隨著深度增加，盡管層厚也增加，但由于此時深度起主導(dǎo)作用，參數(shù)不確定性隨著深度增加而變大.

4 結(jié)論與討論

本文以陸地實測布格重力異常與WGM2012布格重力異常模型融合得到的布格重力異常作為研究依據(jù)，收集了該區(qū)已有的地震接收函數(shù)結(jié)果作為反演先驗?zāi)Ｐ徒⒌幕A(chǔ)，利用貝葉斯原理構(gòu)建密度參數(shù)后驗?zāi)Ｐ?，并通過MCMC采樣算法提取密度參數(shù)的后驗樣本，獲得了川南地區(qū)的地殼三維密度結(jié)構(gòu)模型，并評估了反演結(jié)果的不確定性特征.通過對研究區(qū)內(nèi)4次地震震中位置的8條剖面進(jìn)行分析，本文給出了研究區(qū)內(nèi)4次地震震源區(qū)的密度結(jié)構(gòu)特征.

本文主要研究結(jié)論如下：

(1)通過對地面實測剖面數(shù)據(jù)和WGM2012布格重力異常模型采用等效源方法進(jìn)行融合，得到的融合布格重力異常在測點分布區(qū)的空間分辨率可達(dá)7 km，異常變化范圍-300 mGal至-80 mGal，馬邊、大關(guān)、長寧地震震中都位于重力異常梯度帶區(qū)域.

(2)基于貝葉斯統(tǒng)計學(xué)模型，以地震接收函數(shù)結(jié)果為參考，反演獲得了川南地區(qū)地殼三維密度結(jié)構(gòu)，布格異常擬合殘差均值小于3 mGal，密度參數(shù)不確定性最高達(dá)0.011 g·cm-3.

(3)對比研究區(qū)內(nèi)馬邊、大關(guān)、長寧及瀘縣四個6級以上地震集中分布區(qū)的震源物性結(jié)構(gòu)，發(fā)現(xiàn)馬邊和大關(guān)兩次地震震源區(qū)密度結(jié)構(gòu)不均勻性較強(qiáng)，長寧地震震源區(qū)具備一定的不均勻性特征，但在瀘縣地震震源區(qū)及深部沒有發(fā)現(xiàn)明顯的密度不均體.

(4)2019年的長寧MS6.0地震和2021年瀘縣MS6.0地震震級較大，除了斷層活化或誘發(fā)因素外，在震源區(qū)底部及周邊可能存在深部孕震構(gòu)造背景.

(5)密度參數(shù)的不確定性大小與多種因素有關(guān)，主要為密度體埋深、觀測數(shù)據(jù)豐富程度以及密度體體積大小.隨著密度體埋深增加，重力響應(yīng)的敏感性降低，密度參數(shù)的不確定性增加.在相同深度時，研究區(qū)邊緣的密度體所對應(yīng)的有效觀測數(shù)據(jù)量相對研究區(qū)中部的密度體要少一些，進(jìn)而導(dǎo)致參數(shù)不確定性呈現(xiàn)出邊緣位置略大于研究區(qū)中心位置的特征.密度體體積越大，重力響應(yīng)越大，不確定性也就越小.

致謝感謝三位匿名審稿人對本文提出的建設(shè)性意見和建議.感謝GEOIST開源Python軟件包(https:∥cea2020.gitee.io/geoistdoc)為本文模型測試和反演提供的程序支持.本文部分圖件由pyGMT軟件繪制(Wessel et al., 2019).本文研究中使用的部分陸地實測重力數(shù)據(jù)由湖北省地震局申重陽研究員課題組提供，在此一并感謝.