多源布格重力異常數(shù)據(jù)貝葉斯融合算法及其在川滇地區(qū)的應(yīng)用

李紅蕾， 陳石*， 莊建倉(cāng)， 張貝， 盧紅艷

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

0 引言

布格重力異常數(shù)據(jù)是研究地殼結(jié)構(gòu)和變形的“化石”，是研究地殼結(jié)構(gòu)、斷裂幾何形態(tài)、構(gòu)造邊界和強(qiáng)震孕育環(huán)境的重要基礎(chǔ)地球物理資料(Fu et al., 2014; Motta et al., 2019).近年來(lái)，隨著重力觀測(cè)技術(shù)的不斷進(jìn)步以及觀測(cè)手段的日趨多樣化，布格重力異常的精度和分辨率不斷提高(章傳銀等，2009).然而，布格重力異常由于處理方法和起算基準(zhǔn)等的不同，不同觀測(cè)方式獲得的異常數(shù)據(jù)之間可能存在明顯的系統(tǒng)偏差和顯著的噪聲水平差異.為充分利用多種來(lái)源的布格異常模型，開發(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ù)，以此實(shí)現(xiàn)對(duì)不同來(lái)源數(shù)據(jù)的分析、提取及融合(Pawlowski, 1994; Lane and Peljo， 2004; Oliveira et al.，2013；孫石達(dá)等，2020；孟波和賈真，2020).與數(shù)據(jù)特征擬合類融合算法相比，等效源類融合算法具有以場(chǎng)求源，場(chǎng)源結(jié)合的優(yōu)勢(shì).傳統(tǒng)等效源類融合算法有效實(shí)施的關(guān)鍵是根據(jù)數(shù)據(jù)和場(chǎng)源模型誤差特性確定等效源反演的正則化參數(shù)(Martinez and Li, 2016)，該類正則化參數(shù)一般可以由廣義交叉驗(yàn)證(GCV)或L曲線準(zhǔn)則(Golub et al.，1979；Hansen，1992)獲得，然而，當(dāng)存在誤差特性未知的多個(gè)數(shù)據(jù)源且不同數(shù)據(jù)源之間存在系統(tǒng)性偏差時(shí)，傳統(tǒng)等效源類融合算法的有效實(shí)施變得十分困難.

本文從貝葉斯原理出發(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)化準(zhǔn)則—貝葉斯信息量準(zhǔn)則(Akaike Bayesian Information Criterion ，ABIC準(zhǔn)則)(Akaike，1977)，有效解決傳統(tǒng)融合算法中單一正則化參數(shù)造成的多個(gè)數(shù)據(jù)難以有效融合的問題.

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

1 方法原理

1.1 融合問題

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

圖1 多源重力觀測(cè)網(wǎng)及等效源融合算法示意圖 其中F1及F2代表斷裂；圓點(diǎn)代表衛(wèi)星重力模型格網(wǎng)點(diǎn)，圓點(diǎn)顏色代表了重力異常值的大??；P1及P2代表兩條重力剖面，P1及P2上的三角形代表地面實(shí)測(cè)點(diǎn)，三角形顏色代表了重力異常值的 大??；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)之一是如何最大程度的利用好重力觀測(cè)網(wǎng)中不同類型不同精度的多源觀測(cè)成果，通過不同觀測(cè)信息的相互約束、取長(zhǎng)補(bǔ)短，融合構(gòu)建出局部高精度的融合異常.為此，本文基于等效源反演和貝葉斯不確定性量化理論，發(fā)展了一種新的基于等效源模型的多源重力數(shù)據(jù)貝葉斯優(yōu)化融合算法，用以有效解決不同類型不同精度的多源數(shù)據(jù)有效融合的問題.

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

1.2 算法原理

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

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

=(d1,…,dN)T，

(1)

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

m=(m1,…,mM)T，

(2)

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

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

Gm=d,

(3)

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

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

(4)

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

(5)

(6)

(7)

(8)

2 模型實(shí)驗(yàn)

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

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

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

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

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

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

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

2.2 多源重力基準(zhǔn)偏差修正

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

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

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

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

3 實(shí)際應(yīng)用

3.1 研究背景

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

圖5 研究區(qū)內(nèi)主要構(gòu)造特征與實(shí)測(cè)重力點(diǎn)位分布 其中，紅色實(shí)線為斷裂，白色點(diǎn)線為重力實(shí)測(cè)剖面，其中白色實(shí)心圓為實(shí)測(cè)重力位置，黃色實(shí)心圓為5級(jí)以上地震位置；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寧會(huì)斷裂，F(xiàn)5-2 元謀斷裂， F6安寧河斷裂，F(xiàn)7普渡河斷裂，F(xiàn)8曲江斷裂，F(xiàn)9石屏—建水?dāng)嗔眩現(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).

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

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

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

圖6顯示了地面實(shí)測(cè)布格異常(藍(lán)色實(shí)線)與衛(wèi)星模型計(jì)算布格異常(綠色實(shí)線)數(shù)據(jù)分布.從圖中可以看出，與地面實(shí)測(cè)異常相比，衛(wèi)星模型異常中的高頻擾動(dòng)信息明顯偏多.為了進(jìn)一步分析地面實(shí)測(cè)和衛(wèi)星模型布格異常差值的屬性，我們對(duì)上述兩種不同手段獲得的相同點(diǎn)的布格異常差值進(jìn)行了統(tǒng)計(jì)分析，如圖7所示.

圖6 WGM2012模型和地面實(shí)測(cè)布格重力異常 (a) 剖面P1； (b) 剖面P2.其中藍(lán)色實(shí)線代表剖面點(diǎn)實(shí)測(cè)布格重力異常，黑色實(shí)線代表剖面點(diǎn)地形高度， 橙色實(shí)線代表剖面點(diǎn)自由空氣異常，綠色實(shí)線代表剖面點(diǎn)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)星模型計(jì)算異常和地面實(shí)測(cè)異常差值統(tǒng)計(jì)特征來(lái)看，差值分布服從均值為-120.88 mGal、方差為13.57 mGal、置信區(qū)間為-0.14～0.14 mGal的正態(tài)分布.進(jìn)一步從圖7b兩種異常的差值分布來(lái)看，異常差值具有良好的高斯分布特性.

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

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

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

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

圖9 川滇地區(qū)布格重力異常圖 (a) WGM2012模型布格重力異常(修正基準(zhǔn)差)； (b) 20 km尺度標(biāo)準(zhǔn)布格異常.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é)方差矩陣計(jì)算公式獲得融合所得布格異常不確定度在12.25～30.67 mGal之間，其中地面觀測(cè)約束較多的區(qū)域不確定度較小，地面觀測(cè)缺乏的區(qū)域不確定度較大.

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

圖10 川滇地區(qū)布格重力異常剖面圖 (a) 剖面P1； (b) 剖面P2.其中藍(lán)色實(shí)線代表剖面點(diǎn)實(shí)測(cè)布格重力異常，綠色實(shí)線代表剖面點(diǎn)WGM2012模型布格重力異常(修正基準(zhǔn)差)， 橙色實(shí)線代表剖面點(diǎn)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è)計(jì)了一種新的多源重力異常貝葉斯優(yōu)化融合算法.通過模型測(cè)試了該算法的可行性、有效性和優(yōu)越性.并應(yīng)用該算法對(duì)川滇地區(qū)WGM2012衛(wèi)星格網(wǎng)數(shù)據(jù)和實(shí)測(cè)剖面重力數(shù)據(jù)進(jìn)行了實(shí)驗(yàn)，獲得了研究區(qū)內(nèi)匹配20km結(jié)構(gòu)異常分辨能力的標(biāo)準(zhǔn)融合布格重力異常數(shù)據(jù).

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

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

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

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

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

致謝感謝中國(guó)科學(xué)臺(tái)陣項(xiàng)目提供的實(shí)測(cè)重力資料，感謝兩位匿名審稿人提供的寶貴建議和幫助.