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

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

0 引言

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

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

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

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

1 方法原理

1.1 融合問題

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

圖1 多源重力觀測網(wǎng)及等效源融合算法示意圖 其中F1及F2代表斷裂；圓點代表衛(wèi)星重力模型格網(wǎng)點，圓點顏色代表了重力異常值的大??；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.

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

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

1.2 算法原理

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

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

=(d1,…,dN)T，

(1)

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

m=(m1,…,mM)T，

(2)

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

觀測點處的重力異?？梢员硎緸?3)式的線性方程組形式：

Gm=d,

(3)

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

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

(4)

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

(5)

(6)

(7)

(8)

2 模型實驗

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

2.1 優(yōu)化權(quán)重等效源融合降噪

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

圖2 不同模型異常及模型噪聲特征 (a) 區(qū)域模型正演理論異常，其中黑色實線為剖面測點位置； (b) 混入噪聲的區(qū)域模型異常； (c) 混入噪聲的剖面測點觀測模型異常和區(qū)域模型異常； (d) 區(qū)域模型異常與剖面測點異常噪聲特征.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.

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

圖3 不同融合算法測試結(jié)果 (a) 三次樣條插值算法融合重力結(jié)果； (b) 等權(quán)等效源融合算法重力結(jié)果； (c) 優(yōu)化權(quán)重等效源融合算法重力結(jié)果.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中的統(tǒng)計結(jié)果可以看出，與傳統(tǒng)的插值算法(Test1)和等權(quán)等效源融合算法(Test2)相比，優(yōu)化權(quán)重等效源融合算法(Test3)具有更好的數(shù)據(jù)融合能力.Test3中融合異常與真實異常的差值大小范圍小于25 mGal，標(biāo)準差小于6.5 mGal，遠遠優(yōu)于Test1中三次樣條插值融合和Test2中等權(quán)等效源融合的結(jié)果.此外，優(yōu)化權(quán)重等效源融合算法給出的優(yōu)化權(quán)重(λd1=14.62、λd2=2.10)與數(shù)據(jù)噪聲水平(σd1=15 ms-2、σd2=2 ms-2)具有較好的一致性，由此可以看出，我們的優(yōu)化算法可以實現(xiàn)對不同重力觀測噪聲水平的有效估計，以此給出的權(quán)重可有效避免高頻噪聲對融合數(shù)據(jù)的污染和對單一數(shù)據(jù)的過度依賴，有效提高融合異常的精度與可信度.綜合圖3融合結(jié)果和表1統(tǒng)計數(shù)據(jù)可知，優(yōu)化權(quán)重等效源融合算法在多源重力數(shù)據(jù)異?；謴?fù)及高頻噪聲剔除方面具有更好的優(yōu)勢.

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

2.2 多源重力基準偏差修正

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

圖4中我們同樣測試了三種不同融合方案，分別為傳統(tǒng)三次樣條插值融合、基準差未作為超參數(shù)的優(yōu)化權(quán)重等效源融合和基準差作為超參數(shù)的優(yōu)化權(quán)重等效源融合，所得結(jié)果分別如圖4a、b、c所示.

圖4 不同融合算法測試結(jié)果 (a) 三次樣條插值融合算法有基準差重力結(jié)果； (b) 有基準差的優(yōu)化權(quán)等效源融合算法重力結(jié)果； (c) 無基準差的優(yōu)化權(quán)重等效源融合算法重力結(jié)果.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 不同算法所得融合重力偏差統(tǒng)計參數(shù)及ABIC參數(shù)特征Table 2 Statistical parameters of gravity deviation and ABIC parameter characteristics obtained by different algorithms

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

3 實際應(yīng)用

3.1 研究背景

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

圖5 研究區(qū)內(nèi)主要構(gòu)造特征與實測重力點位分布 其中，紅色實線為斷裂，白色點線為重力實測剖面，其中白色實心圓為實測重力位置，黃色實心圓為5級以上地震位置；F1-1南汀河斷裂東支，F(xiàn)1-2南汀河斷裂西支，F(xiàn)2畹町—安寧斷裂，F(xiàn)3瀾滄江斷裂，F(xiàn)4-1紅河斷裂北段，F(xiàn)4-2紅河斷裂中段，F(xiàn)5-1寧會斷裂，F(xiàn)5-2 元謀斷裂， F6安寧河斷裂，F(xiàn)7普渡河斷裂，F(xiàn)8曲江斷裂，F(xiàn)9石屏—建水?dāng)嗔?，F(xiàn)10白云—山花斷裂(鄧起東等，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).

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

3.2 研究區(qū)衛(wèi)星及地面重力數(shù)據(jù)處理

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

圖6顯示了地面實測布格異常(藍色實線)與衛(wèi)星模型計算布格異常(綠色實線)數(shù)據(jù)分布.從圖中可以看出，與地面實測異常相比，衛(wèi)星模型異常中的高頻擾動信息明顯偏多.為了進一步分析地面實測和衛(wèi)星模型布格異常差值的屬性，我們對上述兩種不同手段獲得的相同點的布格異常差值進行了統(tǒng)計分析，如圖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衛(wèi)星模型計算異常和地面實測異常差值統(tǒng)計特征來看，差值分布服從均值為-120.88 mGal、方差為13.57 mGal、置信區(qū)間為-0.14～0.14 mGal的正態(tài)分布.進一步從圖7b兩種異常的差值分布來看，異常差值具有良好的高斯分布特性.

圖7 WGM2012模型和地面實測布格重力異常之差 (a) 差值統(tǒng)計特征； (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 融合參數(shù)估計

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

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

圖8 優(yōu)化權(quán)重貝葉斯融合算法參數(shù)特征圖 (a) λd2、λm1、c固定時，ABIC取值超參數(shù)λd1取值變化示意圖；(b) λm1、c固定時，ABIC值隨超參數(shù)λd1和λd2組合取值變化示意圖；(c) λd1、λm1、c固定時，ABIC值隨超參數(shù)λd2取值變化示意圖；(d) λd2、c固定時，ABIC值隨超參數(shù)λd1和λm1組合取值變化示意圖；(e) λd1、c固定時，ABIC值隨超參數(shù)λd2和λm1組合取值變化示意圖；(f) λd1、λd2、c固定時，ABIC值隨超參數(shù)λm1取值變化示意圖；(g) λd2、λm1固定時，ABIC值隨超參數(shù)λd1和c組合取值變化示意圖；(h) λd1、λm1固定時，ABIC值隨超參數(shù)λd2和c組合取值變化示意圖；(i)λd2、λd1固定時， ABIC值隨超參數(shù)λm1和c組合取值變化示意圖; (j) λd1、λd2、λm1時，ABIC值隨超參數(shù)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衛(wèi)星模型計算布格異常數(shù)據(jù)顯示出了潛在的長波特征以及明顯的高頻成分.其中的高頻數(shù)據(jù)除了觀測誤差影響之外，很多是極淺的地表質(zhì)量變化產(chǎn)生的，這兩者對于地球深部結(jié)構(gòu)異常研究來說都應(yīng)該被看作噪聲，當(dāng)噪聲處理不當(dāng)時，深部的長波異常特征很容易被高頻噪聲污染和掩蓋.與圖9a衛(wèi)星模型異常相比，通過貝葉斯優(yōu)化計算獲得的20 km尺度的融合布格異常(圖9b)總體呈現(xiàn)出清晰的南高北低的趨勢性變化，反映出了更多的長波長特征，可更好的反映出實際深部結(jié)構(gòu)的橫向變化特征.與圖9a中衛(wèi)星計算布格異常相比，圖9b中融合異常等值線與板塊邊界和線性斷裂分布具有更好的一致性.

圖9 川滇地區(qū)布格重力異常圖 (a) WGM2012模型布格重力異常(修正基準差)； (b) 20 km尺度標(biāo)準布格異常.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.

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

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

圖10 川滇地區(qū)布格重力異常剖面圖 (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 結(jié)論

本文基于等效源反演原理，通過引入貝葉斯參數(shù)優(yōu)化方法，設(shè)計了一種新的多源重力異常貝葉斯優(yōu)化融合算法.通過模型測試了該算法的可行性、有效性和優(yōu)越性.并應(yīng)用該算法對川滇地區(qū)WGM2012衛(wèi)星格網(wǎng)數(shù)據(jù)和實測剖面重力數(shù)據(jù)進行了實驗，獲得了研究區(qū)內(nèi)匹配20km結(jié)構(gòu)異常分辨能力的標(biāo)準融合布格重力異常數(shù)據(jù).

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

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

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

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

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

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