地球物理学报  2021, Vol. 64 Issue (9): 3232-3245   PDF    
多源布格重力异常数据贝叶斯融合算法及其在川滇地区的应用
李红蕾1,2, 陈石1,2, 庄建仓3, 张贝1,2, 卢红艳1,2     
1. 中国地震局地球物理研究所, 北京 100081;
2. 北京白家疃国家地球科学野外科学观测研究站, 北京 100095;
3. 日本数理统计研究所, 东京 190-8562
摘要:布格重力异常数据是研究地壳结构和变形的重要依据,由于观测手段、处理方法等存在差异,不同手段获得的布格重力异常可能存在数据噪声水平及起算基准不一致,从而引起多源布格重力异常数据难以有效融合的难题.为此,本文在传统等效源重力数据融合方法的基础上,通过引入贝叶斯参数优化准则,提出了一种新的多源布格重力异常融合算法.随后对川滇地区两条实测重力剖面布格异常和WGM2012模型计算布格异常开展了融合实验,结果表明:新算法不仅能显著降低不同异常中的非相干噪声,而且能将多源异常的系统偏差有效纠正.与已有布格重力异常相比,融合所得布格重力异常特征与川滇地区主要线性构造分布形态具有更好的一致性.本文给出的川滇地区融合布格重力异常同时兼顾了地表实测重力的高精度和卫星重力场全球模型的高分辨率,可为研究川滇地区地壳结构和构造区划提供更加有价值的地球物理基础资料.
关键词: 等效源反演      贝叶斯优化      数据融合      布格重力异常      川滇地区     
A Bayesian data fusion algorithm for multi-source Bouguer gravity anomalies and its application in the Sichuan-Yunnan region
LI HongLei1,2, CHEN Shi1,2, ZHUANG JianCang3, ZHANG Bei1,2, LU HongYan1,2     
1. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
2. Beijing Baijiatuan Earth Sciences National Observation and Research Station, Beijing 100095, China;
3. The Institute of Statistical Mathematics, Research Organization of Information and Systems, Tokyo 190-8562, Japan
Abstract: Bouguer gravity anomaly data are important sources for the study of crustal structure and deformation. Due to different observation means, processing methods and so on, there may be obvious inconsistent in noise level and datum among Bouguer anomalies which obtained by different approaches. This leads to the difficulty of effective integration of multi-source Bouguer gravity anomaly data. Therefore, based on the traditional gravity data fusion method that is based on the equivalent-source principle, this paper proposes a new data fusion algorithm for multi-source Bouguer gravity anomaly by introducing the Bayesian optimization criteria. At the same time, a numerical experiment was carried out on two actual gravity profile Bouguer anomalies and the Bouguer anomaly calculated by the WGM2012 model in the Sichuan-Yunnan area. The results show that the new algorithm can not only significantly reduce the incoherent noise in different anomalies, but also effectively correct the system deviation of multi-source anomalies. Compared with the existing results of Bouguer gravity anomaly, the characteristics of the Bouguer gravity anomaly obtained by this data fusion algorithm are more consistent with the distribution of main linear structures in the Sichuan-Yunnan area. The Bouguer gravity anomaly model presented in this paper incorporates to both the high accuracy of ground measured gravity and the high resolution of global satellite gravity model, and provides valuable basic geophysical data for the study of crustal structure and tectonic zoning in the Sichuan-Yunnan region.
Keywords: Inversion for equivalent-source    Bayesian optimization    Data fusion    Bouguer gravity anomaly    The Sichuan-Yunnan region    
0 引言

布格重力异常数据是研究地壳结构和变形的“化石”,是研究地壳结构、断裂几何形态、构造边界和强震孕育环境的重要基础地球物理资料(Fu et al., 2014; Motta et al., 2019).近年来,随着重力观测技术的不断进步以及观测手段的日趋多样化,布格重力异常的精度和分辨率不断提高(章传银等,2009).然而,布格重力异常由于处理方法和起算基准等的不同,不同观测方式获得的异常数据之间可能存在明显的系统偏差和显著的噪声水平差异.为充分利用多种来源的布格异常模型,开发有效地多源重力数据融合算法变得十分重要(Bai et al., 2016; Zhao et al., 2019).

常见的重力数据融合算法可以分为两大类,一种为直接从数据特征出发拟合给定基函数参数的方法,如:插值融合法、最小二乘配置法、球谐函数融合法、回归分析法等(Moritz, 1978Hwang and Parsons, 1995王灼华等,2017郭飞霄等,2017谢曦霖等,2018Ke 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., 1979Hansen,1992)获得,然而,当存在误差特性未知的多个数据源且不同数据源之间存在系统性偏差时,传统等效源类融合算法的有效实施变得十分困难.

本文从贝叶斯原理出发,借鉴在气象、水文等领域的数据融合思想(Jaynes, 2003Lahoz 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沿XY方向观测间隔的最小水平尺寸分别约为XY方向总长度的1/8和1/4,根据采样定理,给出等效源模型最佳水平尺度分别应为XY方向总长度的1/4和1/2,据此我们建立起如图 1中B1—B8所示8个均匀分布的单元体组成的等效源模型.下一节我们将结合该模型,给出多源重力贝叶斯优化融合算法的基本原理和求解方程.

1.2 算法原理

根据重力异常叠加原理,相同观测点重力异常的场源积分信号应该具有一致性.这使得通过构建公共的等效场源体(如图 1,B1—B8)实现对多种来源重力异常信号的分析、提取及融合成为可能.相关融合原理简单描述如下:

假设同一组观测区域,不同手段获得了n种类型的观测数据,观测向量d可以表示为:

(1)

其中,S1S2…,Sn表示异常数据的n种不同来源.在我们的算法种,等效源模型通过直立棱柱单元体进行离散化,假设单元体个数为M,模型参数m可表示为如下一维列向量:

(2)

其中,每个直立棱柱体深度、厚度及水平尺度的确定取决于观测数据分布及等效源模型对多种观测数据的拟合情况.

观测点处的重力异常可以表示为(3)式的线性方程组形式:

(3)

式中,mM×1的一维密度模型向量,dN×1的一维观测矩阵,GN×M的二维核函数矩阵,矩阵中每个元素可以由直立棱柱体重力正演公式获得.

通过多组观测构建出统一的等效源模型,以统一的等效源模型为依据去除观测中不遵循一致关系的噪声和偏差,提取出不同观测中的真实信号,是基于等效源模型实现多源数据融合的基本思想.其中,等效源模型的构建是关键,其可归结为对方程(3)中模型参数m的求解.当考虑有n种不同的观测,假设不同观测之间的基准偏差为ci,根据不适定非唯一性反演理论(Tikhonov and Arsenin, 1977),考虑模型拟合差和模型估计误差之和最小建立反演目标函数:

(4)

其中,di为多源异常数据,Gi为其对应的核函数矩阵,Br为最小约束矩阵,λdiλmci称为反演超参数.λdiλm分别为对数据方差σdi2和模型方差σm2的估计值,其大小代表了方程的正则化水平.其中,λdi与观测数据拟合程度相关,λdi越大代表了数据拟合程度越高;λm与先验假设的准确程度相关,λm越大代表了先验认识的准确程度越高.根据差异性原理,正则化超参数λdiλm最优估计值为预期数据失配时的复制值,即其理想最优值需要与数据和先验噪声水平一致(Hansen,1998).然而,当有多种观测或模型约束参与计算时,在对其噪声水平没有准确认识的前提下,正则化超参数最优值的确定十分困难.

为解决上述问题,Akaike(1977)根据数据同化的基本思想,提出了基于贝叶斯信息准则的参数优化算法,可以实现不同观测及模型约束噪声水平的准确估计.本文我们以此为基础,提出了新的基于等效源的多源重力贝叶斯优化数据融合算法,该算法是在贝叶斯反演方法的基础上发展起来的,算法的基本实现过程如下:

假设观测数据向量(di-Gim)的方差σdi2和模型向量Brm的方差σm2满足高斯分布,根据贝叶斯原理,模型向量m的后验概率密度函数可表示为:

(5)

其中P(d|m)为似然函数,可以表示为P(d|m)=P(m)为先验概率密度函数,可以表示为P(m)=.其中Cdi为方差σdi2的方阵,Cm为方差σm2的方阵,根据等价原理,求使得后验概率最大的m就变成了求使得E(m)最小的mE(m)为:

(6)

根据贝叶斯信息准则,需要建立与数据边缘概率密度函数∫P(d|m)P(m)dm有关的反演目标函数,将多种观测数据、模型及超参数统一到一个计算系统中.新的反演目标函数如(7)式所示:

(7)

其中Nh为超参数的个数.求取新目标函数最小的过程就是求使得模型最大后验概率估计和反演系统不确定性估计同时达到最优的过程.新目标函数最优化求解可以通过求导获得,具体求解方法可以采用非线性最优化迭代方法,如单纯形法或拟牛顿方法.将非线性最优化获得的最优超参数λdi*λm*ci*代入E(m*)即可获得模型最大后验概率估计.

将最优超参数λdi*λm*ci*和最优解估计代入到公式(3),获得融合重力异常表达式为:

(8)

融合异常误差协方差矩阵可表示为∑d2=,其对角线元素值即为融合异常的不确定度.

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的等间隔格网将观测网下方区域进行离散化,所得结果分别如图 3abc所示.从图 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 不同算法所得融合重力偏差统计参数及ABIC参数特征 Table 1 Statistical parameters and ABIC parameter characteristics of fusion gravity by different algorithms

此外,从表 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统计数据可知,优化权重等效源融合算法在多源重力数据异常恢复及高频噪声剔除方面具有更好的优势.

2.2 多源重力基准偏差修正

本节的测试模型与2.1节模型及观测网设置相同,即在同一个观测网中依旧存在两种不同的观测,两种观测依旧加入了与前面测试相同噪声水平的高斯白噪声,不同的是两种观测数据除了在噪声水平上有差异之外,还存在100 mGal的基准偏差.

图 4中我们同样测试了三种不同融合方案,分别为传统三次样条插值融合、基准差未作为超参数的优化权重等效源融合和基准差作为超参数的优化权重等效源融合,所得结果分别如图 4abc所示.传统插值算法得到的融合异常结果(图 4a),虽然在形态上与真实异常(图 2a)具有一定的一致性,但融合异常与真实异常偏差统计参数(=101.56 mGal)远超过了容许误差(2σ1=30 mGal).与传统插值算法相比,基准差未作为超参数的优化权重等效源算法所得融合结果(图 4b),在无基准差测点处的数据可以得到较好的恢复,但在有基准差的其他点数据存在较大偏差.基准差作为超参数的优化等效源算法所得融合异常结果(图 4c),异常的形态、幅度范围与图 2a中真实异常基本一致,且表 2 Test3偏差统计参数显示所得融合异常具有较好的精度.上述测试结果表明,同时将基准差和数据噪声估计作为超参数的贝叶斯优化融合算法可以同时避免数据噪声和基准不一致对融合结果的影响.

图 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 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 (WGM 2012).

图 6显示了地面实测布格异常(蓝色实线)与卫星模型计算布格异常(绿色实线)数据分布.从图中可以看出,与地面实测异常相比,卫星模型异常中的高频扰动信息明显偏多.为了进一步分析地面实测和卫星模型布格异常差值的属性,我们对上述两种不同手段获得的相同点的布格异常差值进行了统计分析,如图 7所示.

图 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.

图 7a卫星模型计算异常和地面实测异常差值统计特征来看,差值分布服从均值为-120.88 mGal、方差为13.57 mGal、置信区间为-0.14~0.14 mGal的正态分布.进一步从图 7b两种异常的差值分布来看,异常差值具有良好的高斯分布特性.

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.292c=-111.75.该优化计算结果表明,在20 km尺度深部结构异常研究中,WGM2012卫星模型计算布格重力异常的噪声水平为15.64 mGal,地面布格重力异常噪声水平为2.78 mGal,两者之间的基准偏差为-111.75 mGal.在此基础上,我们进一步对优化超参数的特征进行了分析,并给出了优化权重参数特征图,如图 8所示.其中,图 8acfj分别给出了当有任意三个超参数固定时,ABIC值随超参数λd1λd2λm1c不同取值的变化.图 8bdeghi则分别对应了当有任意两个超参数固定时,ABIC值随另外两个超参数组合不同取值的变化.图 8中红色实心圆给出的最优参数取值与拟牛顿迭代法优化获得的最优超参数估计值具有很好的一致性,且对应的ABIC最小值均为2420.0138,进一步验证了本文非线性优化所得最优超参数的准确性.将上述最优超参数代入公式(6),通过求解获取最优等效源模型估计,将最优模型估计代入公式(8)即可获得优化融合布格异常(图 9b).

图 8 优化权重贝叶斯融合算法参数特征图 (a) λd2λm1c固定时,ABIC取值超参数λd1取值变化示意图;(b) λm1c固定时,ABIC值随超参数λd1和λd2组合取值变化示意图;(c) λd1λm1c固定时,ABIC值随超参数λd2取值变化示意图;(d) λd2c固定时,ABIC值随超参数λd1λm1组合取值变化示意图;(e) λd1c固定时,ABIC值随超参数λd2λm1组合取值变化示意图;(f) λd1λd2c固定时,ABIC值随超参数λm1取值变化示意图;(g) λd2λm1固定时,ABIC值随超参数λd1c组合取值变化示意图;(h) λd1λm1固定时,ABIC值随超参数λd2c组合取值变化示意图;(i)λd2λd1固定时,ABIC值随超参数λm1c组合取值变化示意图; (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.
图 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.
3.4 布格重力异常特征

图 9a卫星模型计算布格异常数据显示出了潜在的长波特征以及明显的高频成分.其中的高频数据除了观测误差影响之外,很多是极浅的地表质量变化产生的,这两者对于地球深部结构异常研究来说都应该被看作噪声,当噪声处理不当时,深部的长波异常特征很容易被高频噪声污染和掩盖.与图 9a卫星模型异常相比,通过贝叶斯优化计算获得的20 km尺度的融合布格异常(图 9b)总体呈现出清晰的南高北低的趋势性变化,反映出了更多的长波长特征,可更好的反映出实际深部结构的横向变化特征.与图 9a中卫星计算布格异常相比,图 9b中融合异常等值线与板块边界和线性断裂分布具有更好的一致性.

根据融合异常误差协方差矩阵计算公式获得融合所得布格异常不确定度在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) (WGM 2012C) 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下载使用.

致谢  感谢中国科学台阵项目提供的实测重力资料,感谢两位匿名审稿人提供的宝贵建议和帮助.
References
Akaike H. 1977. On entropy maximization principle. Application of Statistics, 543: 27-41.
Bai Y L, Dong D D, Wu S G, et al. 2016. A wavelet transformation approach for multi-source gravity fusion: applications and uncertainty tests. Journal of Applied Geophysics, 128: 18-30. DOI:10.1016/j.jappgeo.2016.03.003
Balmino G, Vales N, Bonvalot S, et al. 2012. Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies. Journal of Geodesy, 86(7): 499-520. DOI:10.1007/s00190-011-0533-4
Braitenberg C, Shum C K. 2017. Geodynamic implications of temporal gravity changes over Tibetan Plateau. Italian Journal of Geosciences, 136(1): 39-49. DOI:10.3301/IJG.2015.38
Chan M T, Yang W C, Yan P, et al. 2020. Two-dimensional spatial data fusion method based on SVD decomposition. Progress in Geophysics (in Chinese), 35(3): 940-949. DOI:10.6038/pg2020DD0351
Chen S, Zheng Q Y, Xu W M. 2015. Joint optimal inversion of gravity and seismic data to estimate crustal thickness of the southern section of the north-south seismic belt. Chinese Journal of Geophysics (in Chinese), 58(11): 3941-3951. DOI:10.6038/cjg20151105
Deng Q D, Zhang P Z, Ran Y K, et al. 2002. Basic characteristics of active tectonics of China. Science in China Series D: Earth Sciences, 46(4): 356-372.
Fu G Y, Gao S H, Freymueller J T, et al. 2014. Bouguer gravity anomaly and isostasy at western Sichuan Basin revealed by new gravity surveys. Journal of Geophysical Research: Solid Earth, 119(4): 3925-3938. DOI:10.1002/2014JB011033
Fu G Y, Zhu Y Q, Gao S H, et al. 2013. Discrepancies between free air gravity anomalies from EGM2008 and the ones from dense gravity/GPS observations at west Sichuan Basin. Chinese Journal of Geophysics (in Chinese), 56(11): 3761-3769. DOI:10.6038/cjg20131117
Golub G H, Heath M, Wahba G. 1979. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21(2): 215-223. DOI:10.1080/00401706.1979.10489751
Guo F X, Sun Z M, Zhai Z H, et al. 2017. Tectonic stress field in Bohai bay area inversed by multi-source gravity observation data. Geomatics and Information Science of Wuhan University (in Chinese), 42(1): 136-142.
Hansen P C. 1992. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Review, 34(4): 561-580. DOI:10.1137/1034115
Hansen P C. 1998. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. Vancouver, BC, Canada: Society for Industrial and Applied Mathematics.
Hwang C, Parsons B. 1995. Gravity anomalies derived from Seasat, Geosat, ERS-1 and TOPEX/POSEIDON altimetry and ship gravity: a case study over the Reykjanes Ridge. Geophysical Journal International, 122(2): 551-568. DOI:10.1111/j.1365-246X.1995.tb07013.x
Jaynes E T. 2003. Probability Theory: The Logic of Science. Cambridge: Cambridge University Press.
Ke B G, Zhang L M, Xu J, et al. 2019. Determination of the mean dynamic ocean topography model through combining multi-source gravity data and DTU15 MSS around China's coast. Advances in Space Research, 63(1): 203-212. DOI:10.1016/j.asr.2018.10.040
Khaki M, Hoteit I, Kuhn M, et al. 2017. Assessing sequential data assimilation techniques for integrating GRACE data into a hydrological model. Advances in Water Resources, 107: 301-316. DOI:10.1016/j.advwatres.2017.07.001
Lahoz W A, Schneider P. 2014. Data assimilation: making sense of Earth Observation. Frontiers in Environmental Science, 2: 16.
Lane R, Peljo M. 2004. Estimating the pre-mining gravity and gravity gradient response of the Broken Hill Ag-Pb-Zn Deposit. //74th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1-4.
Li H L, Fang J, Braitenberg C. 2017. Lithosphere density structure beneath the eastern margin of the Tibetan Plateau and its surrounding areas derived from GOCE gradients data. Geodesy and Geodynamics, 8(3): 147-154. DOI:10.1016/j.geog.2017.02.007
Liu F X, Wang W Y, Ji X L, et al. 2021. The fusion with gravity and magnetic potential field data at multi-dimension and multi-scale. Chinese Journal of Geophysics (in Chinese), 64(4): 1453-1470. DOI:10.6038/cjg2021O0235
Liu J Z, Liang X H, Ye Z R, et al. 2020. Combining multi-source data to construct full tensor of regional airborne gravity gradient disturbance. Chinese Journal of Geophysics (in Chinese), 63(8): 3131-3143. DOI:10.6038/cjg2020O0044
Martinez C, Li Y G. 2016. Denoising of gravity gradient data using an equivalent source technique. Geophysics, 81(4): G67-G79. DOI:10.1190/geo2015-0379.1
Meng B, Jia Z. 2020. Equivalent source optimization strategy for gravity anomaly profile transformations. Progress in Geophysics (in Chinese), 35(2): 623-633. DOI:10.6038/pg2020DD0287
Moritz H. 1978. Least-squares collocation. Reviews of Geophysics, 16(3): 421-430. DOI:10.1029/RG016i003p00421
Motta J G, De Souza Filho C R, Carranza E J M, et al. 2019. Archean crust and metallogenic zones in the Amazonian Craton sensed by satellite gravity data. Scientific Reports, 9: 2565. DOI:10.1038/s41598-019-39171-9
Oliveira V C Jr, Barbosa V C F, Uieda L. 2013. Polynomial equivalent layer. Geophysics, 78(1): G1-G13. DOI:10.1190/geo2012-0196.1
Pawlowski R S. 1994. Green's equivalent-layer concept in gravity band-pass filter design. Geophysics, 59(1): 69-76. DOI:10.1190/1.1443535
Sun S D, Du J S, Chen C, et al. 2020. Nonlinear equivalent source method for transformation and inversion of total-field magnetic anomaly. Chinese Journal of Geophysics (in Chinese), 63(1): 351-361. DOI:10.6038/cjg2020M0325
Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. New York: Winston: 1-30.
Wang Q S, Teng J W, Wang G J, et al. 2007. The correction for special pattern of Bouguer gravity anomaly in Heng Duan Mts area by using satellite gravity. Progress in Geophysics (in Chinese), 22(2): 345-352.
Wang Z H, Jin H L, Fu G Y, et al. 2017. EGM2008 model and the measured ground gravity data fusion and its applicaiton in the western Sichuan region. Progress in Geophysics (in Chinese), 2017, 32(3): 1051-1058, doi: 10.6038/pg20170314.
Xie X L, Xu C J, Wen Y M, et al. 2018. A refined least squares collocation method based on multiquadric function. Geomatics and Information Science of Wuhan University (in Chinese), 43(4): 592-598.
Zhang C Y, Guo C X, Chen J Y, et al. 2009. EGM 2008 and its application analysis in Chinese Mainland. Acta Geodaetica et Cartographica Sinica (in Chinese), 38(4): 283-289.
Zhao J X, Wan J H, Sun Q T, et al. 2019. Multi-source ocean gravity anomaly data fusion processing method. //IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium. Yokohama, Japan: IEEE, 8275-8278.
陈勐韬, 杨文采, 颜萍, 等. 2020. 基于SVD分解的二维空间数据融合方法. 地球物理学进展, 35(3): 940-949. DOI:10.6038/pg2020DD0351
陈石, 郑秋月, 徐伟民. 2015. 南北地震带南段地壳厚度重震联合最优化反演. 地球物理学报, 58(11): 3941-3951. DOI:10.6038/cjg20151105
邓起东, 张培震, 冉勇康, 等. 2002. 中国活动构造基本特征. 中国科学(D辑), 32(12): 1020-1030.
付广裕, 祝意青, 高尚华, 等. 2013. 川西地区实测自由空气重力异常与EGM2008模型结果的差异. 地球物理学报, 56(11): 3761-3769. DOI:10.6038/cjg20131117
郭飞霄, 孙中苗, 翟振和, 等. 2017. 利用多源重力观测数据反演渤海湾地区构造应力场. 武汉大学学报·信息科学版, 42(1): 136-142.
刘福香, 王万银, 纪晓琳, 等. 2021. "多维"多尺度重磁位场数据融合方法. 地球物理学报, 64(4): 1453-1470. DOI:10.6038/cjg2021O0235
刘金钊, 梁星辉, 叶周润, 等. 2020. 融合多源数据构建区域航空重力梯度扰动全张量. 地球物理学报, 63(8): 3131-3143. DOI:10.6038/cjg2020O0044
孟波, 贾真. 2020. 面向剖面重力异常转换的等效源优化策略. 地球物理学进展, 35(2): 623-633. DOI:10.6038/pg2020DD0287
孙石达, 杜劲松, 陈超, 等. 2020. 基于等效源的总强度磁异常非线性处理方法. 地球物理学报, 63(1): 351-361. DOI:10.6038/cjg2020M0325
王谦身, 滕吉文, 王光杰, 等. 2007. 应用卫星重力信息对横断山系地区布格重力异常特异分布的纠正. 地球物理学进展, 22(2): 345-352. DOI:10.3969/j.issn.1004-2903.2007.02.003
王灼华, 金红林, 付广裕, 等. 2017. EGM2008模型与地表观测重力数据的融合及其在川西地区的应用. 地球物理学进展, 32(3): 1051-1058. DOI:10.6038/pg20170314
谢曦霖, 许才军, 温扬茂, 等. 2018. 一种基于多面函数的改进最小二乘配置方法. 武汉大学学报·信息科学版, 43(4): 592-598.
章传银, 郭春喜, 陈俊勇, 等. 2009. EGM 2008地球重力场模型在中国大陆适用性分析. 测绘学报, 38(4): 283-289. DOI:10.3321/j.issn:1001-1595.2009.04.001