多源布格重力異常數據貝葉斯融合算法及其在川滇地區的應用

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

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

0 引言

布格重力異常數據是研究地殼結構和變形的“化石”，是研究地殼結構、斷裂幾何形態、構造邊界和強震孕育環境的重要基礎地球物理資料(Fu et al., 2014; Motta et al., 2019).近年來，隨著重力觀測技術的不斷進步以及觀測手段的日趨多樣化，布格重力異常的精度和分辨率不斷提高(章傳銀等，2009).然而，布格重力異常由于處理方法和起算基準等的不同，不同觀測方式獲得的異常數據之間可能存在明顯的系統偏差和顯著的噪聲水平差異.為充分利用多種來源的布格異常模型，開發有效地多源重力數據融合算法變得十分重要(Bai et al.,2016; Zhao et al.,2019).

常見的重力數據融合算法可以分為兩大類，一種為直接從數據特征出發擬合給定基函數參數的方法，如：插值融合法、最小二乘配置法、球諧函數融合法、回歸分析法等(Moritz,1978；Hwang and Parsons, 1995；王灼華等，2017；郭飛霄等，2017；謝曦霖等，2018；Ke et al., 2019；劉金釗等，2020；陳勐韜等，2020；劉福香等，2021)；另一種為以等效源為基礎的數據融合方法，其基本原理是通過構建特定的等效源幾何模型，利用反演方法確定等效單元體參數，以此實現對不同來源數據的分析、提取及融合(Pawlowski, 1994; Lane and Peljo， 2004; Oliveira et al.，2013；孫石達等，2020；孟波和賈真，2020).與數據特征擬合類融合算法相比，等效源類融合算法具有以場求源，場源結合的優勢.傳統等效源類融合算法有效實施的關鍵是根據數據和場源模型誤差特性確定等效源反演的正則化參數(Martinez and Li, 2016)，該類正則化參數一般可以由廣義交叉驗證(GCV)或L曲線準則(Golub et al.，1979；Hansen，1992)獲得，然而，當存在誤差特性未知的多個數據源且不同數據源之間存在系統性偏差時，傳統等效源類融合算法的有效實施變得十分困難.

本文從貝葉斯原理出發，借鑒在氣象、水文等領域的數據融合思想(Jaynes, 2003；Lahoz and Schneide, 2014; Khaki et al., 2017)，發展了一種新的多源重力數據貝葉斯優化融合算法.該算法在傳統等效源融合算法的基礎上(Martinez and Li, 2016)，通過引入貝葉斯參數優化準則—貝葉斯信息量準則(Akaike Bayesian Information Criterion ，ABIC準則)(Akaike，1977)，有效解決傳統融合算法中單一正則化參數造成的多個數據難以有效融合的問題.

本文首先給出了多源重力貝葉斯融合算法的基本原理和求解方程.其次，通過模型試驗對融合算法的可行性及有效性進行了驗證.再次，以川滇地區兩條實測重力剖面布格重力異常數據為約束，基于新算法對該地區已有高分辨率(2′×2′)的WGM2012模型布格重力異常數據進行有效修正，獲得了川滇地區20 km格網分辨率的高精度融合布格異常.新算法有效實現了高精度陸地實測重力剖面數據與高分辨率全球重力場模型數據之間的融合，解決了已有地球重力場模型在川滇地區由于地面先驗觀測缺乏、全球解算和地形校正過度等原因而導致的數據精度過低的問題(Balmino et al., 2012).并以所得融合布格重力異常為依據對川滇地區布格重力異常形態與線性構造之間的對應關系進行了簡要分析.最后，對本文提出的數據融合算法和認識進行了總結和討論.

1 方法原理

1.1 融合問題

假設有這樣一組由衛星模型格網點和實測重力剖面點組成的重力觀測網，如圖1所示.測網中圓點代表了衛星模型格網點位置，從圓點的位置分布可以看出，衛星模型格網點的點位分布均勻且全區覆蓋但點距間隔較大，圖中圓點顏色代表了重力異常值的大小，從其顏色分布可以看出，該測網中衛星重力異常從NW到SE顏色逐漸由藍變紅，異常值呈現趨勢性平緩變化.此外，圖1中三角形給出了實測重力剖面點的點位分布，從圖中可以看出，該測網中的實測重力剖面點沿著與P1、P2兩條測線垂直F1、F2斷裂密集分布；其顏色同樣代表了異常值的大小，從衛星重力異常(圓點)和地面實測重力異常(三角形)顏色分布可以看出，地面實測重力剖面異常可以給出更多重力變化的細節特征，但其容易受到局部高頻噪聲的干擾(如P1測線中F1斷裂附近的地面三角點異常受局部干擾較大).此外，從兩種數據相同區域的異常值對比來看，兩種異常值之間還存在一個明顯的系統偏差.圖1中的重力測網示例圖從實際問題的角度出發，很好地給出了同一觀測系統中，不同觀測手段在數據分布、分辨率和噪聲水平上的復雜性.為了更好的利用同一觀測系統中的多源數據，發展針對多源復雜觀測的科學數據融合算法成為當前多源重力數據應用中亟待解決的重要問題之一.

圖1 多源重力觀測網及等效源融合算法示意圖 其中F1及F2代表斷裂；圓點代表衛星重力模型格網點，圓點顏色代表了重力異常值的大小；P1及P2代表兩條重力剖面，P1及P2上的三角形代表地面實測點，三角形顏色代表了重力異常值的 大小；B1—B8代表等效源.Fig.1 Schematic diagram of multi-source gravity observation network and equivalent source fusion algorithm F1 and F2 are faults. Dots are grid nodes of the satellite gravity model, and the color of dots represents the magnitude of gravity anomalies. P1 and P2 represent two gravity profiles. Triangles on P1 and P2 represent the measurement points on the ground, and the triangle color represents the magnitude of gravity anomalies. B1—B8 represent the equivalent sources.

多源重力數據融合算法面臨的重大挑戰之一是如何最大程度的利用好重力觀測網中不同類型不同精度的多源觀測成果，通過不同觀測信息的相互約束、取長補短，融合構建出局部高精度的融合異常.為此，本文基于等效源反演和貝葉斯不確定性量化理論，發展了一種新的基于等效源模型的多源重力數據貝葉斯優化融合算法，用以有效解決不同類型不同精度的多源數據有效融合的問題.

根據等效源類數據融合算法的基本原理，等效源劃分尺寸的大小可以通過觀測間隔的最小尺寸給出，圖1沿X和Y方向觀測間隔的最小水平尺寸分別約為X和Y方向總長度的1/8和1/4，根據采樣定理，給出等效源模型最佳水平尺度分別應為X和Y方向總長度的1/4和1/2，據此我們建立起如圖1中B1—B8所示8個均勻分布的單元體組成的等效源模型.下一節我們將結合該模型，給出多源重力貝葉斯優化融合算法的基本原理和求解方程.

1.2 算法原理

根據重力異常疊加原理，相同觀測點重力異常的場源積分信號應該具有一致性.這使得通過構建公共的等效場源體(如圖1，B1—B8)實現對多種來源重力異常信號的分析、提取及融合成為可能.相關融合原理簡單描述如下：

假設同一組觀測區域，不同手段獲得了n種類型的觀測數據，觀測向量d可以表示為：

=(d1,…,dN)T，

(1)

其中，S1，S2…，Sn表示異常數據的n種不同來源.在我們的算法種，等效源模型通過直立棱柱單元體進行離散化，假設單元體個數為M，模型參數m可表示為如下一維列向量：

m=(m1,…,mM)T，

(2)

其中，每個直立棱柱體深度、厚度及水平尺度的確定取決于觀測數據分布及等效源模型對多種觀測數據的擬合情況.

觀測點處的重力異常可以表示為(3)式的線性方程組形式：

Gm=d,

(3)

式中，m為M×1的一維密度模型向量，d為N×1的一維觀測矩陣，G為N×M的二維核函數矩陣，矩陣中每個元素可以由直立棱柱體重力正演公式獲得.

通過多組觀測構建出統一的等效源模型，以統一的等效源模型為依據去除觀測中不遵循一致關系的噪聲和偏差，提取出不同觀測中的真實信號，是基于等效源模型實現多源數據融合的基本思想.其中，等效源模型的構建是關鍵，其可歸結為對方程(3)中模型參數m的求解.當考慮有n種不同的觀測，假設不同觀測之間的基準偏差為ci，根據不適定非唯一性反演理論(Tikhonov and Arsenin, 1977)，考慮模型擬合差和模型估計誤差之和最小建立反演目標函數：

(4)

為解決上述問題，Akaike(1977)根據數據同化的基本思想，提出了基于貝葉斯信息準則的參數優化算法，可以實現不同觀測及模型約束噪聲水平的準確估計.本文我們以此為基礎，提出了新的基于等效源的多源重力貝葉斯優化數據融合算法，該算法是在貝葉斯反演方法的基礎上發展起來的，算法的基本實現過程如下：

(5)

(6)

(7)

(8)

2 模型實驗

本節分別設計了兩組仿真測試模型，以此來檢驗多源重力異常貝葉斯優化融合算法對不同觀測噪聲水平準確量化，基準偏差精確求解，提高融合數據精度的能力.

2.1 優化權重等效源融合降噪

本節在1000 km×1000 km的區域范圍內，設計了一個由均勻格網和非均勻剖面測點組成的重力觀測網.在該區域觀測網下方給定了4個埋深和分布各不相同的異常體，異常體密度分別為0.4 g·cm-3、0.425 g·cm-3、0.45 g·cm-3和0.5 g·cm-3，區域背景密度為0 g·cm-3，區域模型和剖面觀測測點分辨率分別為25 km和5 km.其中，區域模型正演理論異常及實測點位分布如圖2a所示，在此理論異常的基礎上，加入均值為0 mGal、標準差為15 mGal高斯噪聲，混入高斯噪聲的區域模型異常如圖2b所示.圖2c分別給出了混入噪聲的剖面測點正演異常和同樣測點的區域模型異常，其中藍色和紅色實線代表了較高分辨率(約5 km)的測點觀測異常，其噪聲水平為2 mGal；黑色和綠色實線代表了較低分辨率(約20 km)的區域模型異常，其噪聲水平為15 mGal.藍色和紅色實線代表的實測異常與黑色和綠色實線代表的區域模型異常相比，不存在虛假的大的異常波動，并且可以給出異常源更多細節變化信息.區域模型異常和實測剖面異常噪聲水平及分布特征如圖2d所示.

圖2 不同模型異常及模型噪聲特征 (a) 區域模型正演理論異常，其中黑色實線為剖面測點位置； (b) 混入噪聲的區域模型異常； (c) 混入噪聲的剖面測點觀測模型異常和區域模型異常； (d) 區域模型異常與剖面測點異常噪聲特征.Fig.2 Different model anomalies and model noise characteristics (a) Theoretical anomaly of regional model. The solid black lines are measurement profiles； (b) Regional model anomaly mixed with noise； (c) Observation model anomaly and regional model anomaly mixed with noise； (d) Characteristics of regional model anomaly and profile measurement point anomaly.

在算法測試過程中，我們選擇了三種不同的融合方案，分別為三次樣條插值融合、等權等效源融合和優化權重等效源融合，等效源融合算法是以20 km的等間隔格網將觀測網下方區域進行離散化，所得結果分別如圖3a、b、c所示.從圖3中的融合結果可以看出，不同融合算法所得融合結果差異性顯著.圖3a融合異常與圖2a正演理論異常相比，高頻噪聲依舊存在，且從表1的統計結果可以看出融合異常值與真實值相比偏差較大，高達55 mGal；與圖3a相比，圖3b的高頻噪聲明顯減少，但局部融合異常會存在較大偏差，如右上角的融合異常與真實異常相比存在較大偏差；圖3c優化權重等效源方法獲得的融合異常，相比前兩個結果，在異常形態、偏差及高頻噪聲剔除方面都有較好的改善.

圖3 不同融合算法測試結果 (a) 三次樣條插值算法融合重力結果； (b) 等權等效源融合算法重力結果； (c) 優化權重等效源融合算法重力結果.Fig.3 Test results of different fusion algorithms (a) Fusion gravity by the cubic spline interpolation algorithm； (b) Fusion gravity by the equivalent-source algorithm with equal weights； (c) Fusion gravity by the equivalent source algorithm with optimized weights.

此外，從表1中的統計結果可以看出，與傳統的插值算法(Test1)和等權等效源融合算法(Test2)相比，優化權重等效源融合算法(Test3)具有更好的數據融合能力.Test3中融合異常與真實異常的差值大小范圍小于25 mGal，標準差小于6.5 mGal，遠遠優于Test1中三次樣條插值融合和Test2中等權等效源融合的結果.此外，優化權重等效源融合算法給出的優化權重(λd1=14.62、λd2=2.10)與數據噪聲水平(σd1=15 ms-2、σd2=2 ms-2)具有較好的一致性，由此可以看出，我們的優化算法可以實現對不同重力觀測噪聲水平的有效估計，以此給出的權重可有效避免高頻噪聲對融合數據的污染和對單一數據的過度依賴，有效提高融合異常的精度與可信度.綜合圖3融合結果和表1統計數據可知，優化權重等效源融合算法在多源重力數據異常恢復及高頻噪聲剔除方面具有更好的優勢.

表1 不同算法所得融合重力偏差統計參數及ABIC參數特征Table 1 Statistical parameters and ABIC parameter characteristics of fusion gravity by different algorithms

2.2 多源重力基準偏差修正

本節的測試模型與2.1節模型及觀測網設置相同，即在同一個觀測網中依舊存在兩種不同的觀測，兩種觀測依舊加入了與前面測試相同噪聲水平的高斯白噪聲，不同的是兩種觀測數據除了在噪聲水平上有差異之外，還存在100 mGal的基準偏差.

圖4中我們同樣測試了三種不同融合方案，分別為傳統三次樣條插值融合、基準差未作為超參數的優化權重等效源融合和基準差作為超參數的優化權重等效源融合，所得結果分別如圖4a、b、c所示.

圖4 不同融合算法測試結果 (a) 三次樣條插值融合算法有基準差重力結果； (b) 有基準差的優化權等效源融合算法重力結果； (c) 無基準差的優化權重等效源融合算法重力結果.Fig.4 Test results of different fusion algorithms (a) Cubic spline interpolation fusion algorithm with datum drift； (b) Optimized weight equivalent source fusion algorithm with datum drift； (c) Optimized weight equivalent-source fusion algorithm without datum drift.

表2 不同算法所得融合重力偏差統計參數及ABIC參數特征Table 2 Statistical parameters of gravity deviation and ABIC parameter characteristics obtained by different algorithms

以上兩個模型的測試結果表明了新算法在重力數據噪聲水平估計、基準偏差測定、多源數據融合方面的可行性、有效性和優越性.下面我們將選擇地球物理探測需求和重力觀測資料較多的川滇地區，測試算法在實際重力資料融合中的效果，并為相關應用研究提供有益參考.

3 實際應用

3.1 研究背景

本文研究區域位于川滇地區，該地區為青藏高原東南向擠出的構造邊界，地殼結構復雜，內部褶皺斷裂廣泛發育，包括了豐富的大地構造區邊界斷裂和控制強震活動的活斷層，如紅河斷裂、安寧河斷裂、曲江斷裂等(如圖5所示).布格重力對該區地殼結構、斷裂信息等研究具有重要作用，目前該區布格重力異常及相關構造結構研究成果主要來源于衛星重力場模型(陳石等，2015; Li et al., 2017)，但由于缺少地面約束導致已有衛星模型數據精度有限，相關構造及變形的研究精度受到較多限制(付廣裕等，2013).

圖5 研究區內主要構造特征與實測重力點位分布 其中，紅色實線為斷裂，白色點線為重力實測剖面，其中白色實心圓為實測重力位置，黃色實心圓為5級以上地震位置；F1-1南汀河斷裂東支，F1-2南汀河斷裂西支，F2畹町—安寧斷裂，F3瀾滄江斷裂，F4-1紅河斷裂北段，F4-2紅河斷裂中段，F5-1寧會斷裂，F5-2 元謀斷裂， F6安寧河斷裂，F7普渡河斷裂，F8曲江斷裂，F9石屏—建水斷裂，F10白云—山花斷裂(鄧起東等，2002).Fig.5 Major structural features and gravity measurement points in the study area Red solid lines are faults. White dotted lines represent gravity measured profiles. White solid circles represent measured gravity positions. Yellow solid circles represent earthquakes with magnitude greater than 5. F1-1:east branch of Nantinghe fault, F1-2:west branch of Nantinghe fault, F2:Wanding-Anning fault, F3: Lancangjiang fault, F4-1: northern segment of Red River fault, F4-2: middle segment of Red River fault, F5-1: Huining fault, F5-2: Yuanmou fault, F6:Anning river fault, F7: Pudu river fault, F8: Qujiang fault, F9: Shiping-Jianshui fault, F10:Baiyun-Shanhua fault (Deng et al., 2002).

近年來，國家重力臺網中心、中國地震科學臺陣探測等項目在川滇地區布設了多個重力觀測網.地面實測重力具有高精度的特點，可以實現對衛星模型的有效約束(Fu et al., 2014).中國地震科學臺陣計劃I期觀測項目在研究區內實測重力點位如圖5中白色實心圓所示.從圖5中的測點分布來看，由于地形復雜，在研究區內依舊存在很多重力測量空區，對于其中的測量空區，不能采用直接外推或擬合方法進行填補.現今通用的多種高階地球重力場模型數據(http:∥icgem.gfz-potsdam.de/home)具有全球覆蓋的優勢，可以作為重力測量空區一種很好的補充(王謙身等，2007).通過球諧計算、濾波、校正等處理手段可以獲得一定精度的重力異常用于大尺度的重力異常應用研究(Braitenberg and Shum, 2017)，但隨著高精度地球內部結構和變形研究的需要，此類處理手段獲得的異常數據不可避免的會代入一些系統誤差、全球平均誤差以及濾波算法本身帶來的誤差等，這些都會對高精度重力異常的應用研究產生影響.因此，綜合利用研究區內多源重力數據成果，在考慮數據時空分布和觀測場及背景場誤差的基礎上，采用以場求源，場源結合的等效源算法，依據貝葉斯參數優化理論，發揮不同重力觀測手段的優勢，融合構建出研究區內高精度的布格重力數據，對于該地區地殼結構和構造區劃等諸多前沿科學問題的研究具有重要的意義.

3.2 研究區衛星及地面重力數據處理

首先，我們收集了中國地震科學臺陣計劃在川滇地區沿NE-SW走向、平均采樣間隔約為2 km的兩條實測重力剖面，沿剖面共獲得了394個局部高精度流動重力和GPS觀測值.通過對各觀測點正常場改正、高度改正、布格板改正及地形改正計算得到各測點布格重力異常，結果如圖6中藍色實線所示.同時，收集BGI(Bureau Gravimétrique Internationa)官方網站發布的2′× 2′的WGM2012模型計算布格重力異常作為補充進行融合計算.WGM2012模型計算布格重力異常是通過全球重力模型EGM2008和DTU10計算獲得，模型同時引入了高精度的全球高程模型 ETOP01模型來計算地表起伏的貢獻(Balmino et al.,2012).通過比較圖6中地面實測布格重力異常(藍色實線)和相同點位的WGM2012模型布格異常(綠色實線)可以看出，兩種不同來源的布格重力異常數據除了存在覆蓋范圍、分辨率及精度的差異之外，還存在著高達百mGal的系統偏差，也可以稱之為基準偏差.造成其基準偏差的原因中，地形校正參數的差異是其中很重要的一方面，包括了地形校正參數分布范圍、分辨率、精度等的差異.這種基準偏差在傳統的融合和濾波處理中很難被修正，而我們的算法可以通過將基準偏差作為求解超參數，通過優化算法對其進行準確求解，以此大大提高融合結果精度.

圖6顯示了地面實測布格異常(藍色實線)與衛星模型計算布格異常(綠色實線)數據分布.從圖中可以看出，與地面實測異常相比，衛星模型異常中的高頻擾動信息明顯偏多.為了進一步分析地面實測和衛星模型布格異常差值的屬性，我們對上述兩種不同手段獲得的相同點的布格異常差值進行了統計分析，如圖7所示.

圖6 WGM2012模型和地面實測布格重力異常 (a) 剖面P1； (b) 剖面P2.其中藍色實線代表剖面點實測布格重力異常，黑色實線代表剖面點地形高度， 橙色實線代表剖面點自由空氣異常，綠色實線代表剖面點WGM2012模型布格異常.Fig.6 Bouguer gravity anomalies derived from WGM2012 model and ground measurements (a) Profile P1； (b) Profile P2. The blue solid line represents the observed Bouguer gravity anomaly at the profile point (OBS), the black solid line represents the topography height at the profile point (TOPO), the orange solid line represents the free air anomaly at the profile point (FGA), and the green solid line represents the WGM2012 model Bouguer anomaly at the profile point (WGM2012).

從圖7a衛星模型計算異常和地面實測異常差值統計特征來看，差值分布服從均值為-120.88 mGal、方差為13.57 mGal、置信區間為-0.14～0.14 mGal的正態分布.進一步從圖7b兩種異常的差值分布來看，異常差值具有良好的高斯分布特性.

圖7 WGM2012模型和地面實測布格重力異常之差 (a) 差值統計特征； (b) 差值分布范圍.Fig.7 The difference between Bouguer gravity anomalies derived from WGM2012 model and ground measurement (a) Statistical characteristics of differences； (b) Range of difference distribution.

3.3 融合參數估計

基于上述分析可得，研究區域內地面實測和WGM2012模型計算布格重力噪聲服從高斯分布的情況下，我們同時以兩種數據噪聲水平估計值λd1、λd2，先驗約束權重估計值λm1以及系統偏差c作為超參數.利用研究區內高分辨率的WGM2012模型計算格網布格重力和394個局部高精度實測剖面重力異常進行融合計算.由于實測剖面最大觀測點距可達10 km，根據觀測數據分布特征和采樣定理，我們在99°E—103°E、22°N—27°N水平范圍和10 km以內的深度范圍內，以水平尺寸為20 km和垂向高度為100 m的等效源單元對研究區域進行離散化構建融合模型.

在融合模型求解過程中，通過擬牛頓最優化迭代算法求解獲得的四個最優超參數分別為：λd1=1/15.642、λd2=1/2.782、λm1=1/0.292、c=-111.75.該優化計算結果表明，在20 km尺度深部結構異常研究中，WGM2012衛星模型計算布格重力異常的噪聲水平為15.64 mGal，地面布格重力異常噪聲水平為2.78 mGal，兩者之間的基準偏差為-111.75 mGal.在此基礎上，我們進一步對優化超參數的特征進行了分析，并給出了優化權重參數特征圖，如圖8所示.其中，圖8a、c、f、j分別給出了當有任意三個超參數固定時，ABIC值隨超參數λd1、λd2、λm1、c不同取值的變化.圖8b、d、e、g、h、i則分別對應了當有任意兩個超參數固定時，ABIC值隨另外兩個超參數組合不同取值的變化.圖8中紅色實心圓給出的最優參數取值與擬牛頓迭代法優化獲得的最優超參數估計值具有很好的一致性，且對應的ABIC最小值均為2420.0138，進一步驗證了本文非線性優化所得最優超參數的準確性.將上述最優超參數代入公式(6)，通過求解獲取最優等效源模型估計，將最優模型估計代入公式(8)即可獲得優化融合布格異常(圖9b).

圖8 優化權重貝葉斯融合算法參數特征圖 (a) λd2、λm1、c固定時，ABIC取值超參數λd1取值變化示意圖；(b) λm1、c固定時，ABIC值隨超參數λd1和λd2組合取值變化示意圖；(c) λd1、λm1、c固定時，ABIC值隨超參數λd2取值變化示意圖；(d) λd2、c固定時，ABIC值隨超參數λd1和λm1組合取值變化示意圖；(e) λd1、c固定時，ABIC值隨超參數λd2和λm1組合取值變化示意圖；(f) λd1、λd2、c固定時，ABIC值隨超參數λm1取值變化示意圖；(g) λd2、λm1固定時，ABIC值隨超參數λd1和c組合取值變化示意圖；(h) λd1、λm1固定時，ABIC值隨超參數λd2和c組合取值變化示意圖；(i)λd2、λd1固定時， ABIC值隨超參數λm1和c組合取值變化示意圖; (j) λd1、λd2、λm1時，ABIC值隨超參數c取值變化示意圖.Fig.8 Schematic diagrams of Bayesian parameter features of Bayesian fusion algorithm with optimal weights (a) ABIC value changed by hyperparameter λd1 with fixed λd2,λm1 and c；(b) ABIC value changed by combination of hyperparameter λd1 and λd2 with fixed λm1 and c； (c) ABIC value changed by hyperparameter λd2 with fixed λd1,λm1 and c; (d) ABIC value changed by combination of hyperparameter λd1 and λm1 with fixed λd2 and c; (e) ABIC value changed by combination of hyperparameter λd2 and λm1 with fixed λd1 and c; (f) ABIC value changed by hyperparameter λm1 with fixed λd1,λd2 and c； (g) ABIC value changed by combination of hyperparameter λd1 and c with fixed λd2 and λm1; (h) ABIC value changed by combination of hyperparameter λd2 and c with fixed λd1 and λm1； (i) ABIC value changed by combination of hyperparameter λm1 and c with fixed λd1 and λd2； (j) ABIC value changed by hyperparameter c with fixed λd1, λd2 and λm1.

3.4 布格重力異常特征

圖9a衛星模型計算布格異常數據顯示出了潛在的長波特征以及明顯的高頻成分.其中的高頻數據除了觀測誤差影響之外，很多是極淺的地表質量變化產生的，這兩者對于地球深部結構異常研究來說都應該被看作噪聲，當噪聲處理不當時，深部的長波異常特征很容易被高頻噪聲污染和掩蓋.與圖9a衛星模型異常相比，通過貝葉斯優化計算獲得的20 km尺度的融合布格異常(圖9b)總體呈現出清晰的南高北低的趨勢性變化，反映出了更多的長波長特征，可更好的反映出實際深部結構的橫向變化特征.與圖9a中衛星計算布格異常相比，圖9b中融合異常等值線與板塊邊界和線性斷裂分布具有更好的一致性.

圖9 川滇地區布格重力異常圖 (a) WGM2012模型布格重力異常(修正基準差)； (b) 20 km尺度標準布格異常.Fig.9 Bouguer gravity anomaly map in the Sichuan-Yunnan region (a) Bouguer gravity anomaly of WGM2012 model (with corrected datum drift)； (b) Standard Bouguer anomaly of 20 km resolution.

根據融合異常誤差協方差矩陣計算公式獲得融合所得布格異常不確定度在12.25～30.67 mGal之間，其中地面觀測約束較多的區域不確定度較小，地面觀測缺乏的區域不確定度較大.

圖10給出了實測剖面P1和P2觀測點位置的優化融合異常(橙色實線)、衛星模型異常(綠色實線)和地面觀測異常(藍色實線)對比曲線.與地面觀測異常(藍色實線)相比，優化融合異常(橙色實線)不僅可以有效消除地面實測異常中虛假的大的異常波動，而且可以較好地體現出實測異常的趨勢性變化.與衛星模型異常(綠色實線)相比，優化融合異常(橙色實線)不僅可以修正衛星模型異常的基準偏差，而且可以有效去除衛星模型異常中的高頻擾動，使得更符合地質實際的深部異常長波特征得以有效保留.

圖10 川滇地區布格重力異常剖面圖 (a) 剖面P1； (b) 剖面P2.其中藍色實線代表剖面點實測布格重力異常，綠色實線代表剖面點WGM2012模型布格重力異常(修正基準差)， 橙色實線代表剖面點20 km尺度融合后的布格異常.Fig.10 Profile of Bouguer gravity anomaly in Sichuan-Yunnan region (a) Profile P1； (b) Profile P2. The blue solid line represents the observed Bouguer gravity anomaly (OBS) at the measurement point, the green solid line represents the WGM2012 model Bouguer gravity anomaly (corrected reference difference) (WGM2012C) at the measurement point, and the orange solid line represents fused Bouguer gravity anomaly(FBGA) of 20 km resolution at the measurement point.

4 結論

本文基于等效源反演原理，通過引入貝葉斯參數優化方法，設計了一種新的多源重力異常貝葉斯優化融合算法.通過模型測試了該算法的可行性、有效性和優越性.并應用該算法對川滇地區WGM2012衛星格網數據和實測剖面重力數據進行了實驗，獲得了研究區內匹配20km結構異常分辨能力的標準融合布格重力異常數據.

本文得到的主要研究結論如下：

(1)綜合測試結果表明，貝葉斯優化融合算法給出的最優權重與數據噪聲水平具有較好的一致性，以此給出的最優等效源估計可有效避免高頻噪聲對融合數據的污染，廣泛提高融合數據的可信度與準確度.將基準差作為超參數的貝葉斯優化融合算法可有效避免基準不統一對融合結果的影響.與傳統的插值算法相比，貝葉斯優化融合算法在多源重力融合、高頻噪聲估計與剔除、基準歸算中具有更好的優勢.

(2)通過新的多源重力貝葉斯優化融合算法，以20 km分辨率的等間距格網對川滇地區WGM2012衛星模型數據和實測重力剖面數據進行融合，獲得了該區匹配20 km結構異常分辨能力的標準融合布格異常，標準異常與地面實測異常和WGM2012模型異常相比，非相干高頻誤差大大減少，不同觀測基準偏差得以有效消除，所得標準布格異常更有利于地球深部結構異常的研究.

(3)在20 km尺度構造異常的研究中，WGM2012衛星模型布格重力異常噪聲水平為15.64 mGal，地面布格重力異常噪聲水平為2.78 mGal.衛星模型布格重力與地面布格重力基準偏差為-111.75 mGal.融合所得標準布格異常不確定度在12.25～30.67 mGal之間，其中在地面觀測約束較多的附近區域不確定度較小，在地面觀測缺乏的區域附近不確定度較大.

本文提出的多源重力貝葉斯優化融合算法，除了能有效解決由于數據噪聲水平和起算基準差異造成的多源重力異常難以有效融合的難題之外，也存在一些問題，比如：優化超參數選取與設置問題、參數優選中計算量的問題等.我們將繼續針對融合參數優化，融合計算能力提升等方面做進一步的改進和深入研究.本文融合的川滇地區布格重力異常模型可以通過https:∥gitee.com/cea2020/geodataset/blob/master/bug/BGAmodel.dat下載使用.

致謝感謝中國科學臺陣項目提供的實測重力資料，感謝兩位匿名審稿人提供的寶貴建議和幫助.