2019, Vol. 62
2. 中国地震局地球物理研究所, 北京 100081;
3. 日本国立统计数理研究所, 东京 1908562
2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
3. The Institute of Statistical Mathematics, Tokyo 1908562, Japan
对重力观测值进行固体潮、气压、地形、中间层、高度和正常场等一系列校正后获得的布格重力异常是进行地球物理反演与解释的重要资料.重力观测值经过初步整理后，再进行高度校正和正常场校正，得到自由空气重力异常.其中，高度校正和正常场校正仅与测点高度和纬度有关，可通过相关公式直接计算得出，进一步对自由空气重力异常进行地形校正和中间层校正，最终得到布格重力异常.地形校正和中间层校正不仅与测点高程有关，还与近地表岩石密度分布情况密切相关.岩石密度取值不当会导致布格重力异常存在误差，影响后续反演与解释工作.
实际区域重力资料整理中，一般采用地壳平均密度2.67g·cm-3作为布格校正密度，但平均地壳密度常与实际近地表岩石密度存在差异，使得校正不完全得到的布格重力异常并不能很好反映实际地下密度结构.野外露头采样与钻井测量可直观反映测区岩石密度分布(Hammer, 1965; Gibb and Thomas, 1980; LaFehr, 1983; Sastry, 2005; Singh et al., 2006).但当研究区面积较大时，有限的露头采样和钻井测量资料并不能全面反映测区地表岩石密度的变化情况.因此，许多学者提出了众多估算近地表密度的方法和技术.Nettleton(1939)提出了确定近地表平均密度的剖面法，通过采用假定平均密度校正得到的布格重力异常与地形起伏的相关性情况估计近地表平均密度.此后，部分学者对最小相关性理论进行了不同程度的改进(Parasnis, 1979; Rikitake et al., 1965; Saha et al., 2006; Jilinski et al., 2013).Tontini等(2007)则根据观测数据径向功率谱斜率估计近地表平均密度.Eshaghzadeh等(2015)结合分型理论和统计学方法估计近地表平均密度.然而由于实际地表密度是变化的，因此平均密度估计方法仍然很难得到校正完全的布格重力异常.
Vajk(1956)解释了地形校正与中间层校正中采用变密度的重要性.Grant等(1962)将近地表密度视为光滑函数，通过对地形与剩余布格重力异常的相关性取最小进行横向变密度估计.Sissons(1981)提出了已知地下重力测量资料时的纵向变密度估计法.还有部分学者提出了基于滑动窗口的横向变密度估计方法(Marc and Jacques, 1983; Rimbert et al., 1987; Moribayashi, 1990).对于大面积研究区，地表地形起伏与深部地质构造存在一定的相关性，采用最小相关理论得到的布格重力异常在一定程度上损失了深部异常信息.针对这一问题，有学者提出了布格异常相对于地表地形起伏光滑分布的假设条件，并采用光滑函数近似拟合布格重力异常.Fukao和Murata分别通过分段常函数和二维三次B样条函数近似拟合布格重力异常(Fukao et al., 1981; Murata, 1993)，Murata同时采用了基于ABIC准则的客观贝叶斯方法(Akaike, 1977)确定数据拟合和先验约束之间的权重系数，同时得到每个滑动窗口内的平均密度和布格异常，但未考虑不同窗口之间布格异常的光滑程度.Nawa在Murata方法的基础上作出改进，得到了考虑了不同窗口之间布格重力异常光滑程度的横向变密度分布和布格重力异常(Nawa et al., 1997).
本文基于Nawa等(1997)方法同样假设布格重力异常相对于近地表地形光滑分布，从剖面自由空气重力异常数据出发，研究了自动优化估算剖面近地表岩石横向密度的技术.该方法通过地表重力测量，可以给出近地表岩层参考密度信息，这对于构建浅层地壳密度模型，精化地壳模型的属性信息等都具有一定意义.本文通过均匀剖分模型和非均匀剖分模型数据试验，验证方法的有效性和准确性，同时讨论样条系数个数对估计结果准确性的影响，以及异常噪声水平不同时，ABIC准则对数据拟合和先验光滑约束之间权重系数自动调整的效果.最后通过云南地区的两条剖面实际数据试验验证本文方法的有效性，并分析其近地表密度分布特征.
1 方法原理 1.1 近地表密度估计的反演方程已知测点自由空气重力异常时，需作地形校正和中间层校正来获取布格重力异常.地形校正是为了去除测点周围地形起伏对测点产生的引力效应，无论测点周围地形起伏情况如何，起伏地形都将导致测点重力值减小.经过地形校正后，各测点周围一定范围内的起伏地形的影响被消除，但各个不同高程的测点距离大地水准面之间物质层的影响依然存在，因此需作中间层校正去除这一影响(Wang, 2003; Zeng, 2005).本文所述近地表的范围是指从地表面到大地水准面之间的物质层.
假设剖面各测点的自由空气重力异常值[F1, F2, …, FN]T已知，N为测点总数.根据高精度地形网格数据，通过常规公式可计算得到各个测点的单位密度地形校正值[T1, T2, …, TN]T和中间层校正值[M1, M2, …, MN]T(Kane, 1962; Nagy, 1966; 王谦身, 2003; 曾华霖, 2005).对剖面作等长度的分区，每个分区具有相应的密度值ρi，则剖面的近地表密度结构可表示为向量形式[ρ1, ρ2, …, ρK]T，K为密度分区总数.根据各个测点与所在密度分区的对应关系，将单位密度地形校正值和单位密度中间层校正值变形为矩阵形式TN×K、MN×K后，可得如下关系式：
|
(1) |
式(1)就是近地表密度估计的反演方程.其中，左侧FN为已知的实际自由空气异常数据，右侧HN×K为综合考虑了地形校正和中间层校正的单位密度校正系数.右侧BN和ρK是反演的模型参数，分别为未知的布格重力异常和近地表密度结构.公式(1)的反演方程是非线性方程，需要采用非线性算法优化求解.
1.2 布格重力异常的三次B样条光滑约束假设布格重力异常具有光滑性且与起伏地形并不相关，采用具有良好光滑性质的三次B样条函数(Inoue, 2012)拟合测点布格重力异常，进而对公式(1)的估计问题起到光滑约束作用.
将剖面划分为(Mx-1)个等长度的区间，共Mx个样条节点.测点的布格重力异常可用三次B样条函数表示为
|
(2) |
Ri(x)为M个分段样条基函数，si为M个样条系数，其中M=Mx+2(Inoue, 2012).将(2)式整理为矩阵形式：
|
(3) |
EN×M中的每一个元素为Ep, ξ=R4-k(rxp)，p=1, 2, …, N，ξ=i+k(k=0, 1, 2, 3)，i为测点所在的样条区间，rx为测点所在样条区间上的归一化坐标(Inoue, 2012).样条基函数矩阵EN×M仅与样条剖分方式以及测点位置有关，而样条系数向量sM×1未知.此时，(1)式中反演模型参数转变为sM和ρK，反演方程(1)可整理为：
|
(4) |
假设布格重力异常拟合曲线具有一阶、二阶光滑的性质.公式(5)、(6)是衡量连续布格重力异常曲线一阶、二阶光滑程度的变量，分别称作平坦度和光滑度：
|
(5) |
|
(6) |
对其进行离散化处理，可分别转化为矩阵形式：‖D1(Mx-1)×MsM×1‖2和‖D2(Mx-2)×MsM×1‖2，其中：
|
(7) |
|
(8) |
EMx×M为样条节点坐标对应的样条基函数矩阵.
1.3 模型参数的最大后验估计基于客观贝叶斯分析的反演理论(Tarantola, 1987)，通过贝叶斯公式从概率统计的角度将地球物理反演方程、模型参数的先验约束信息联系起来，以获得模型参数的后验估计，同时得到模型参数的统计特征.
实测自由空气异常中的噪声服从正态分布，且各测点的布格重力异常值与拟合曲线上函数值之差也服从正态分布，反演方程(4)可看成均值为0，协方差为Σd2=σd2IN×N的多元正态分布.给定模型参数时，数据拟合项[F-(Hρ+Es)]的条件概率密度函数可表示为
|
(9) |
假设布格重力异常具有一阶、二阶光滑的性质，一阶导数分布[D1(Mx-1)×MsM×1]服从均值为0，协方差为Σ12=σ12I(Mx-1)(Mx-1)的(Mx-1)维正态分布，二阶导数[D2(Mx-2)×MsM×1]服从均值为0，协方差为Σ22=σ22I(Mx-2)(Mx-2)的(Mx-2)维正态分布，将两个分布综合表示为
|
(10) |
为后续计算方便，引入变量
|
(11) |
则光滑约束项的概率密度函数为
|
(12) |
其中[]+代表矩阵的正定部分，其中D=
根据式(9)、(12)及贝叶斯公式(Tarantola, 1987)，模型参数的后验概率密度函数可表示为
|
(13) |
(13) 式取最大时对应的模型参数sM和ρK即为最佳估计值，这一过程称为模型参数的最大后验估计.(13)式取最大即为指数部分取最大，也即对下式取最小：
|
(14) |
其中，a=[F1, F2, …, FN, 0, 0, …, 0]T，v是模型参数向量v=[ρ1, ρ2, …, ρK, s1, s2, …, sM]T，Z为(N+2Mx-3)×(K+M)维矩阵：
|
(15) |
当Z矩阵中的超参数w1、w2确定时，根据广义反演理论，模型参数向量v的最大后验估计值为
|
(16) |
根据模型参数的后验概率密度函数：
|
(17) |
可得模型参数估计误差的协方差矩阵Σm2=
超参数w1、w2代表了估计问题中数据拟合和光滑约束之间的权重，当超参数w1、w2的值确定时，可得到确定模型参数的最大后验估计解v*.本文采用Akaike提出的基于最大熵原理的贝叶斯准则确定超参数的取值(Akaike, 1980)：
|
(18) |
(18) 式可看作是观测数据的对数似然函数.通过(18)式取最小确定超参数的值.将估计密度值的个数也视为超参数.将(9)和(12)式代入(18)式，可得
|
(19) |
(19)式是关于数据拟合项标准差σd2及超参数w1、w2的函数.为求其最小值，首先通过对σd2求偏导，得到后验估计噪声方差
|
(20) |
此时(19)式可表示为
|
(21) |
其次，(21)式取最小可看作是关于超参数w1、w2的非线性最优化问题，对于非线性最优化问题有多种求解方法，例如单纯形函数法、拟牛顿法和模拟退火法等(Kowalik and Osborne, 1968; Davidon, 1975; Kirkpatrick et al., 1983).本文选择单纯形函数算法(Kowalik and Osborne, 1968)求解该优化问题.
2 模型数据试验及参数讨论根据近地表密度分区均匀和分区不均匀情况，将模型分为均匀剖分模型和非均匀剖分模型两种.不均匀分布的测点位置、剖面地形和单位密度校正系数、光滑连续的布格重力异常曲线以及近地表密度结构直接给出.剖面地形及单位密度校正系数分别是从某区域实际地形数据中提取、计算得到.正演计算得到原始自由空气重力异常，并添加指定标准差的高斯噪声，确定样条系数个数后，进行近地表密度结构估计.本文还针对样条系数个数的取值问题、噪声水平不同时ABIC准则对数据拟合和光滑约束之间权重自动调整的有效性问题，进行了参数讨论.
2.1 均匀剖分模型试验本文设计了如图 1所示的均匀剖分模型，测点数为100，剖面地形起伏情况如图 1a所示，单位密度校正系数如图 1b所示.根据给定的光滑连续曲线，提取相应测点的理论布格异常值，如图 1c所示，同时给出了如图 1d所示的理论近地表密度结构，共划分为10个等长度的小区块，密度分布范围在2.2~2.7 g·cm-3之间.根据理论模型信息，可正演计算原始自由空气异常，添加标准差为0.2 mGal的高斯噪声，得到如图 1e所示的含噪声原始自由空气重力异常.
|
图 1 均匀剖分模型 (a)剖面地形; (b)剖面单位密度校正系数; (c)剖面各测点理论布格重力异常; (d)理论近地表密度结构; (e)含噪声的原始自由空气异常(噪声标准差为0.2 mGal). Fig. 1 Uniform partition model (a) Topography of the profile; (b) The correction coefficient per unit density of the profile; (c) The theoretical Bouguer anomaly of the profile; (d) The theoretical density distribution; (e) Original Free-air anomaly with Gaussian noise (The standard deviation of Gaussian noise is 0.2 mGal). |
样条系数个数的选取接近测点数时，可获得稳定的估计结果(后续章节有详细描述)，因此，设计样条系数个数为120.根据给定模型数据，首先采用本文前述方法进行超参数搜索，当ABIC值最小时，对应的超参数为w1=8.78×10-5，w2=4.11×10-1，超参数搜索路径如图 2所示.当超参数取该值时，得到如图 3a所示的估计近地表密度结构，图中同时给出了误差上下限并同时得到如图 3b中的估计光滑布格重力异常，估计的光滑连续布格重力异常曲线与各测点理论布格重力值整体上拟合较好，图 3b中同时给出了采用常密度2.67 g·cm-3校正后的布格重力异常，显然，采用估计近地表密度结构校正后的布格重力异常与采用常密度校正后的布格重力异常存在明显偏差，最大可以达到15 mGal左右，如图 3c所示，可见考虑近地表密度结构横向变化计算得到的布格重力异常包含了更丰富的重力信息.根据估计密度值与估计光滑布格异常，可正演计算得到自由空气重力异常，如图 3d所示，其与含噪声的原始自由空气重力异常之间的最大偏差为0.39 mGal，估计的高斯噪声标准差为0.197 mGal，接近理论噪声标准差.观察图 3a中的估计误差上下限可见，180 km以西地区近地表密度的估计误差有轻微增大的现象，我们认为这与该区域地形起伏不大有关，地形起伏不够剧烈导致了单位密度校正系数较小，因而本文中的估计方法对密度差异不够敏感.
|
图 2 均匀剖分模型试验的超参数搜索图 Fig. 2 The search path of hyperparameters in uniform partition model test (The hyperparameters are 8.78×10-5 and 4.11×10-1, and the minimum ABIC value was 244.31) |
|
图 3 均匀剖分模型试验的估计结果 (a)估计近地表密度结构(红色实线)，估计近地表密度结构加减估计误差(灰色区域)，理论近地表密度结构(黑色实线)；(b)估计光滑连续布格重力异常(红色实线)，理论布格重力异常(黑色实心圆点)，采用常密度校正结果(蓝色实心圆点)；(c)估计光滑连续布格异常曲线上相应测点处的布格异常值与常密度2.67 g·cm-3校正得到的布格异常值之差；(d)估计结果正演自由空气异常(红色菱形线)，含有噪声的原始自由空气异常(黑色方形线). Fig. 3 The estimated results of uniform partition model (a) Estimated density value (red solid line), estimated density value with estimated error (gray area), theoretical density value (black solid line); (b) Estimated smooth continuous Bouguer gravity anomaly (red solid line), theoretical Bouguer gravity anomaly (black solid dot), correction result by constant density value of 2.67 g·cm-3 (blue line with solid dot); (c) The difference between the Bouguer anomaly at each observation station read from smooth continuous Bouguer gravity anomaly curve and corrected by constant value of 2.67 g·cm-3; (d) Free-air anomaly forward calculated by estimated result (red line with rhombus), original Free-air anomaly with Gaussian noise (black line with square). |
考虑到实际近地表岩石密度分布具有不均匀性，本文设置了密度分区不均匀的模型.剖面地形数据(同图 1a)、单位密度校正系数(同图 1b)与理论布格异常(同图 1c)与均匀剖分模型试验保持一致，理论近地表密度结构如图 4a所示，同样添加标准差为0.2 mGal的高斯噪声，得到含噪声的原始自由空气重力异常如图 4b所示.
|
图 4 非均匀剖分模型 (a)理论近地表密度结构；(b)含有高斯噪声的原始自由空气异常(噪声标准差为0.2 mGal). Fig. 4 Un-uniform partition model (a) Theoretical density distribution; (b) Original Free-air anomaly with Gaussian noise (the standard deviation of Gaussian noise is 0.2 mGal). |
确定超参数的值为w1=2.57×10-4，w2=5.3×10-1，得到如图 5a所示的估计近地表密度结构及误差上下限；同时得到如图 5b中的估计光滑连续布格重力异常，与测点理论布格重力异常值拟合程度较好，图 5b中同时给出了采用常密度2.67 g·cm-3校正后的布格重力异常，采用估计近地表密度结构校正得到的布格重力异常与采用常密度校正后的布格重力异常之间的偏差，最大可以达到13 mGal左右，如图 5c所示，可见估计布格重力异常可揭示更精细的重力信息.根据估计密度值与估计光滑布格异常，可正演计算得到自由空气重力异常, 如图 5d所示，其与含噪声的原始自由空气重力异常之间的最大偏差为0.94 mGal，估计的高斯噪声标准差为0.3 mGal.
|
图 5 非均匀剖分模型试验的估计结果 (a)估计近地表密度结构(红色实线)，近地表密度结构加减估计误差(灰色区域)，理论近地表密度结构(黑色实线)；(b)估计光滑连续布格重力异常(红色实线)，理论布格重力异常(黑色实心圆点)，采用常密度校正结果(蓝色实心圆点)；(c)估计光滑连续布格异常曲线上相应测点处的布格异常值与常密度2.67 g·cm-3校正得到的布格异常值之差；(d)估计结果正演自由空气异常(红色菱形线)，含有噪声的原始自由空气异常(黑色方形线). Fig. 5 The estimated result of non-uniform partition model (a) Estimated density value (red solid line), estimated density value with estimated error (gray area), theoretical density value (black solid line); (b) Estimated smooth continuous Bouguer gravity anomaly (red solid line), theoretical Bouguer gravity anomaly (black solid dot), correction result by constant density value of 2.67 g·cm-3 (blue line with solid dot); (c) The difference between the Bouguer anomaly at each observation station read from smooth continuous Bouguer gravity anomaly curve and corrected by constant value of 2.67 g·cm-3; (d) Free-air anomaly forward calculated by estimated result (red line with rhombus), original Free-air anomaly with Gaussian noise (black line with square). |
为了进一步检验样条系数个数不同对估计结果的影响以及ABIC准则控制超参数自动搜索的有效性，本文分别改变了均匀剖分模型试验中的样条系数个数和异常噪声水平，分析估算结果.
2.3.1 样条系数个数影响为了检验样条系数个数不同对估计结果的影响，使样条系数的个数从50变化到200，最小ABIC值的变化情况，如图 6a所示；估计噪声标准差的变化情况，如图 6b所示；各个密度子区间的估计值以及估计误差的变化情况，如图 7所示.可以看出，随着样条系数个数的增加，最小ABIC值逐渐减小，当样条系数个数小于100时，估计标准差、估计各个子区间密度值与估计密度误差有小范围的波动，而当样条系数的个数大于100时，这些变量的值趋于稳定.可见，足够多的样条系数可以保证估计结果的准确性，但是样条系数个数设置过大，会增加计算量.因此，样条系数个数接近测点数为最佳.
|
图 6 (a) 最小ABIC值随M的变化情况；(b)估计标准差随M的变化情况 Fig. 6 (a) The changes of the minimum ABIC value with the increase of M value; (b) The changes of the estimated standard deviation of data residual with the increase of M value |
|
图 7 各个子区块的估计密度值(A)及估计密度误差随M的变化情况(B)(黑色实心方形：估计值；黑色实线：理论值) Fig. 7 With the increase of M value, the changes of estimated density in each subdomain (A; black line with square: estimated value; black solid line: theoretical value) and estimated density error (B) |
为验证噪声水平不同时，ABIC准则对数据拟合和先验光滑约束之间权重的自动调整，本文改变了均匀剖分模型试验中的噪声水平，设置高斯噪声标准差分别为σd=0.4 mGal和0.6 mGal，相应的超参数搜索结果为(w1=2.48×10-4，w2=8.56×10-1)，(w1=8.11×10-4，w2=2.52)，相应的估计噪声标准差为0.32 mGal和0.56 mGal.从超参数及估计噪声标准差值的变化可见，随着噪声水平的增加，超参数的值增加，估计噪声标准差也逐渐增大，且都接近理论值，说明模型先验光滑约束在反演估计中所占权重增大.可见在近地表密度估计过程中，ABIC准则可根据噪声水平的不同，自动调整数据拟合和光滑约束之间的权重.
图 8a展示了根据不同噪声水平的原始自由空气异常所估计的近地表密度结构；图 8b展示了估计误差，可看出随着噪声水平的增加，估计误差逐渐增大，当噪声水平达到0.6 mGal时，估计密度误差可达到0.1 g·cm-3；图 8c展示了不同的估计光滑布格重力异常，均与理论模型吻合程度较好.
|
图 8 不同噪声水平的原始自由空气异常对应的估计结果 (a)估计密度值(图中红色、紫色和蓝色实线)，理论密度值；(b)估计密度误差；(c)估计光滑连续布格重力异常，理论布格重力异常. Fig. 8 Different estimated results of original Free-air anomaly with different noise level (a) Estimated density value; (b) Estimated density error; (c) Estimated smooth continuous Bouguer anomaly, theoretical Bouguer anomaly. |
噪声水平、单位密度校正系数的大小均会影响估计结果的稳定性.观察图 1b中180 km以西地区单位密度校正系数值可知，近地表密度差异为0.1 g·cm-3时，由近地表地形起伏及中间层影响产生的引力效应最小达0.8 mGal左右，因此当噪声标准差接近这一值时，噪声的影响会大过地形起伏和中间层的影响，降低算法对密度差异的敏感性，导致算法不稳定，因此在理论模型试验中，不易添加过高水平的噪声.
3 实际数据试验对中国地震局在云南地区测量完成的两条绝对重力控制下的相对重力剖面进行实际数据试验，剖面位置如图 9所示，两条剖面均呈西南-北东走向.研究区属于南北地震带南段，位于由扬子地块、川滇地块、思茅地块组成的青藏高原东南缘构造活动区.受多期次、反复的地质构造运动影响，研究区形成了极其复杂的地质构造环境，广泛存在大面积出露的前寒武纪基底以及晚古生代以来巨大的近南北向断裂(缪以琨等, 1986；崔作舟等, 1987；葛碧如和杨科佑, 1990).在数次构造运动中，丰富的岩浆活动使得该区域的多数断裂被岩浆岩填充，显著的岩浆特征以及丰富的矿产资源使得攀西地区得到了国内外地质和地球物理学家的关注.
|
图 9 (a) 研究区地形图及(b)研究区地质构造简图 Fig. 9 (a) The topography map of research area; (b) The geological structure sketch of research area |
剖面AA′自云县附近至会东附近，在南涧—弥渡一带穿越红河断裂带.投影后的一维剖面全长约44.6 km，共198个测点.为了估计该剖面的近地表密度结构，我们设置了200个样条节点，将剖面剖分成25个密度区块，得到估计密度结果以及估计布格重力异常，如图 10中d1—f1所示.剖面BB′自思茅附近至七甸附近，在墨江—元江一带穿越红河断裂带.投影后的一维剖面全长约35.3 km，共195个测点.为了估计该剖面的近地表密度结构，我们设置了200个样条节点，将剖面剖分成20个密度区块.得到估计密度结果以及估计布格重力异常，如图 10中d2—f2所示.
|
图 10 两条重力剖面的已知信息及估计结果(左列：AA′剖面，右列：BB′剖面) (a1)—(f1)和(a2)—(f2)为AA′和BB′剖面的地形图、单位密度校正系数、原始自由空气重力异常、估计密度值、估计密度误差、估计布格重力异常与Geosoft软件校正得到的布格重力异常(图d1与d2中，灰色区域为估计密度值加减估计误差结果，黑色箭头为主要断裂带在剖面上的投影；图f1与f2中，黑色实线为估计布格重力异常，黑色点划线为采用Geosoft软件校正得到的布格重力异常). Fig. 10 The known information and estimated results of two gravity profiles (left column: profile AA′, right column: profile BB′) (a1)—(f1) in left column and (a2)—(f2) in right column respectively represented the topography, correction parameter per unit, original Free-air gravity anomaly, estimated density, estimated density error and estimated Bouguer gravity anomaly of profile AA′ and profile BB′. (in fig. d1 and d2: gray areas are estimated density results adding or subtracting estimated density error, black arrows are projection of faults; in fig. f1 and f2: the black full lines are estimated Bouguer gravity anomaly, the black dotted lines are Bouguer gravity anomaly corrected by const method in Geosoft Oasis software). |
根据图 9所示的云南地区地质图上反映的岩性特征，可见前寒武纪、古生代、中生代在两条剖面所经过的地区均有覆盖，新生代地层覆盖范围较少.AA′剖面上前寒武纪地层主要分布在剖面西端云县以西地区以及剖面东端永仁以东—会东以西地区，估计密度值在2.53~2.74 g·cm-3之间；BB′剖面上前寒武纪地层在元江以东—晋宁一带广泛出露，除去元江—玉溪之间出现高达2.79~2.83 g·cm-3的密度值，其余部分的估计结果在2.74 g·cm-3左右.古生代地层在AA′剖面弥渡—祥云附近以及会东地区有少量分布，估计密度值在2.58~2.6 g·cm-3，在BB′剖面晋宁以东出露，估计密度值在2.7~2.74 g·cm-3.中生代地层在AA′剖面南涧—永仁附近出露，估计密度为2.55~2.73 g·cm-3；在BB′剖面西侧—墨江之间大面积出露，估计密度值在2.45~2.65 g·cm-3之间.对比图 9中两条剖面所经过的地层年代一致的区域，可见估计的近地表密度值相差不大，两条剖面中前寒武纪近地表地层密度大多在2.53~2.75 g·cm-3之间，部分地区出现高达2.8 g·cm-3左右的极大值；古生代近地表地层密度基本在2.58~2.74 g·cm-3；中生代近地表地层密度基本都在2.45~2.73 g·cm-3之间.两条剖面的估计近地表密度结构整体表现出东侧略高于西侧的特征，与图 9显示的东侧地层年代老于西侧地层年代大致吻合.估计的近地表地层密度结果也可体现出图 9中显示的岩浆岩，AA′剖面西端云县及其以东附近有前寒武纪岩浆岩以及古生代岩浆岩出露，估计近地表密度结果在2.59~2.7 g·cm-3.
根据文献查阅及项目收集的云南省地层密度资料(Shellnutt and Zhou, 2007; Jian et al., 2009; Hou et al., 2013; Marian et al., 2013)，云南省前寒武纪地层密度范围在2.5~2.74 g·cm-3，平均密度值为2.71 g·cm-3；古生代寒武纪-泥盆纪地层密度变化范围是2.53~2.72 g·cm-3，石炭纪-二叠纪地层密度变化范围为2.7~2.87 g·cm-3；中生代的三叠纪-侏罗纪地层密度分布在2.5~2.68 g·cm-3，白垩纪地层密度变化范围较大，在2.02~2.78 g·cm-3之间，盐岩平均密度为2.25 g·cm-3，砂岩、泥岩平均密度为2.6 g·cm-3，中生代近地表地层密度整体分布在2.25~2.78 g·cm-3.不同年代地层的估计密度范围与实测资料中大致吻合.在AA′剖面永仁以东至会东以西地区，以及BB′剖面元江—杨武地区及峨山地区存在高达2.8 g·cm-3以上的估计密度值，推测该高值与该区域附近出露的华力西期辉绿岩、闪长岩及二叠系玄武岩有关(Shellnutt and Zhou, 2007；Hou eta al., 2013；石磊等，2015).
红河断裂是攀西地区同时跨越两条剖面的主要断裂带(分别在AA′剖面南涧地区、BB′剖面墨江至元江之间地区穿越)，两条剖面上均显示为断裂带左侧密度值相对低，右侧密度值相对高，BB′剖面红河断裂左右两侧的密度差异大于AA′剖面.AA′剖面上红河断裂带西侧估计密度在2.6 g·cm-3，东侧估计密度在2.84 g·cm-3，BB′剖面上红河断裂西侧估计密度在2.4 g·cm-3，东侧估计密度在2.8 g·cm-3左右.除红河断裂以外，在剖面上几处重要断裂所在区域均可见估计密度结果有一定的跳变，这一特征可作为地表断裂的辅助识别标志.
4 结论本文基于前人的研究成果，提出了一种基于重力贝叶斯分析的剖面近地表岩石横向变密度估计方法，采用基于客观贝叶斯分析的ABIC方法自动搜索数据拟合和先验约束之间的权重系数.均匀剖分模型和非均匀剖分模型试验验证了方法的可靠性，使用该方法可在得到较可靠的近地表岩石密度分布的同时，得到光滑连续分布的剖面布格重力异常.经过进一步的均匀剖分模型试验，我们得到了如下认识：(1)增大样条系数的个数，可得到更稳定、更准确的估计密度值，样条系数个数的确定接近测点数时可得到最佳估计结果；(2)随着异常噪声水平的增大，先验光滑约束相对于数据拟合的权重在ABIC准则的控制下会自动增大，证实了ABIC准则对数据拟合和先验约束之间权重系数自动调整的有效性.云南地区的两条重力剖面估计结果显示，前寒武纪地层近地表密度值在2.53~2.75 g·cm-3之间，古生代地层近地表密度值在2.58~2.74 g·cm-3，中生代地层近地表密度在2.45~2.73 g·cm-3，这一特征与区域物性资料吻合；对于部分估计密度值高达2.8 g·cm-3的结果，我们认为与研究区广泛出露的辉绿岩、闪长岩及二叠系玄武岩有关；估计密度结果中出现跳变的地方可指示地表断裂构造；估计布格重力异常具有良好的光滑性.该方法取得的结果不仅可为云南地区重力反演与解释工作提供包含更精细重力信息的布格异常，而且可进一步为该区域近地表地质构造分析提供可靠的密度约束.在平原地区，由于地形起伏不够剧烈，因此较小的单位密度校正系数可能会降低该方法对密度差异的敏感性，估计密度结果可能不稳定，关于平坦地区下一维近地表密度估计方法还有待改善.
致谢 衷心感谢两位匿名评审专家提出的宝贵意见.
Akaike H. 1980. Likelihood and the Bayes Procedure (with Discussion), Bayesian Statistics. Valencia, Spain: University Press.
|
Cui Z Z, Lu D Y, Chen J P, et al. 1987. The deep structural and tectonic features of the crust in panxi area. Chinese Journal of Geophysics (in Chinese), 30(6): 566-580. |
Davidon W C. 1975. Optimally conditioned optimization algorithms without line searches. Mathematical Programming, 9(1): 1-30. DOI:10.1007/BF01681328 |
Eshaghzadeh A, Kalantari R S, Hekmat Z M. 2015. Optimum density determination for bouguer correction using statistical methods:a case study from north of Iran. International Journal of Advanced Geosciences, 3(2): 25-29. DOI:10.14419/ijag.v3i2.4988 |
Fukao Y, Yamamoto A, Nozaki K. 1981. A method of density determination for gravity correction. Journal of Physics of the Earth, 29(2): 163-166. DOI:10.4294/jpe1952.29.163 |
Ge B R, Yang K Y. 1990. Mesozoic-cenozoic tectonic features in panzhihua-xichang area. Acta Geoghysics Sinica (in Chinese), 33(1): 64-69, 125-126. |
Gibb R A, Thomas M D. 1980. Density determinations of basic volcanic rocks of the Yellowknife supergroup by gravity measurements in mine shafts-Yellowknife, Northwest Territories. Geophysics, 45(1): 18-31. DOI:10.1190/1.1441036 |
Grant F S, Elsaharty A F. 1962. Bouguer gravity corrections using a variable density. Geophysics, 27(5): 616-626. DOI:10.1190/1.1439071 |
Hammer S. 1965. Density determinations by underground gravity measurements-Sequel. Geophysics, 30(6): 1133-1134. DOI:10.1190/1.1439698 |
Hou T, Zhang Z C, Encarnacion J, et al. 2013. The role of recycled oceanic crust in magmatism and metallogeny:Os-Sr-Nd isotopes, U-Pb geochronology and geochemistry of picritic dykes in the Panzhihua giant Fe-Ti oxide deposit, central Emeishan large igneous province, SW China. Contributions to Mineralogy and Petrology, 165(4): 805-822. DOI:10.1007/s00410-012-0836-3 |
Inoue H. 2012. A least-squares smooth fitting for irregularly spaced data:finite-element approach using the cubic B-spline basis. Geophysics, 51(11): 2051-2066. DOI:10.1190/1.1442060 |
Jian P, Liu D Y, Kröner A, et al. 2009. Devonian to Permian plate tectonic cycle of the Paleo-Tethys Orogen in southwest China (Ⅱ):insights from zircon ages of ophiolites, arc/back-arc assemblages and within-plate igneous rocks and generation of the Emeishan CFB province. Lithos, 113(3-4): 767-784. DOI:10.1016/j.lithos.2009.04.006 |
Jilinski P, Fontes S L, Meju M A. 2013. Estimating optimum density for regional Bouguer reduction by morphological correlation of gravity and bathymetric maps:examples from offshore south-eastern Brazil. Geo-Marine Letters, 33(1): 67-73. DOI:10.1007/s00367-012-0312-0 |
Kane M F. 1962. A comprehensive system of terrain corrections using a digital computer. Geophysics, 27(4): 455-462. DOI:10.1190/1.1439044 |
Kirkpatrick S, Gelatt C D Jr, Vecchi M P. 1983. Optimization by simulated annealing. Science, 220(4598): 671-680. DOI:10.1126/science.220.4598.671 |
Kowalik J S, Osborne M R. 1968. Methods of unconstrained optimization problems. Mathematics of Computation, 8(24): 425-427. DOI:10.2307/2004893 |
LaFehr T R. 1983. Rock density from borehole gravity surveys. Geophysics, 48(3): 341-356. DOI:10.1190/1.1441472 |
Marc B, Jacques L. 1983. Determination directe des densites du sol et de remblais a partir de mesures gravimetriques. Bulletin of the International Association of Engineering Geology, 26(1): 171-173. DOI:10.1007/BF02594213 |
Miao Y K, Ren Z D, Chen X Y, et al. 1986. Tectonic environment and tectonic development of Panxi Rift in Sichuan. Earth Science-Journal of Wuhan College of Geology (in Chinese), 11(6): 631-637. |
Moribayashi S. 1990. A new method for variable density correction of gravity data. BUTSURI-TANSA (Geophys. Explor.), 43: 97-106. |
Munteanu M, Yao Y, Wilson A H, et al. 2013. Panxi region (South-West China):Tectonics, magmatism and metallogenesis. A review. Tectonophysics, 608: 51-71. DOI:10.1016/j.tecto.2013.09.008 |
Murata Y. 1993. Estimation of optimum average surficial density from gravity data-an objective bayesian approach. Journal of Geophysical Research:Solid Earth, 98(B7): 12097-12109. DOI:10.1029/93jb00192 |
Nagy D. 1966. The gravitational attraction of a right rectangular prism. Geophysics, 31(2): 362-371. DOI:10.1190/1.1439779 |
Nawa K, Fukao Y, Shichi R, et al. 1997. Inversion of gravity data to determine the terrain density distribution in southwest Japan. Journal of Geophysical Research:Solid Earth, 102(B12): 27703-27719. DOI:10.1029/97jb02543 |
Nettleton L L. 1939. Determination of density for reduction of gravimeter observations. Geophysics, 4(3): 176-183. DOI:10.1190/1.0403176 |
Parasnis D S. 1979. Principles of Applied Geophysics. 3rd ed. London: Chapman and Hall.
|
Rikitake T, Tajima H, Izutuya S, et al. 1965. Gravimetric and geomagnetic studies of Onikobe area. Bul Earthq Res Inst, 43: 241-267. |
Rimbert F, Erling J C, Lakshmanan J. 1987. Variable density Bouguer processing of gravity data from Herault, France. First Break, 5(1): 9-13. DOI:10.3997/1365-2397.1987001 |
Saha D, Chatterjee S M, Baruah B N. 2006. Determination of optimum average density for bouguer correction-a case study in Tirasujanpur Area, Himachal Pradesh, India.//6th International Conference & Exposition on Petroleum Geophysics. Kolkata: Society of Petroleum Geophysicists India, 89-93.
|
Sastry R G, Chaudhary N, Israil M. 2005. Is there a need for variable density option in making combined mass correction to gravity data acquired over high relief?. Current Science, 88(2): 269-273. |
Shellnutt J G, Zhou M F. 2007. Permian peralkaline, peraluminous and metaluminous A-type granites in the Panxi district, SW China:Their relationship to the Emeishan mantle plume. Chemical Geology, 243(3-4): 286-316. DOI:10.1016/j.chemgeo.2007.05.022 |
Shi L, Lou H, Wang Q S, et al. 2015. Gravity field characteristics and crust density structure in the Panxi region, China. Chinese Journal of Geophysics (in Chinese), 58(7): 2402-2412. DOI:10.6038/cjg20150717 |
Singh P, Pandey R S, Gupta S K, et al. 2006. Bouguer reduction with lateral variable surface densities-a model based case study in geologically complex frontier area.//6th International Conference & Exposition on Petroleum Geophysics. Kolkata: Society of Petroleum Geophysicists India, 72-77.
|
Sissons B A. 1981. Densities determined from surface and subsurface gravity measurements. Geophysics, 46(11): 1568-1571. DOI:10.1190/1.1441163 |
Tarantola A. 1987. Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation. Amsterdam: Elsevier Science Publishing Co.
|
Tontini F C, Graziano F, Cocchi L, et al. 2007. Determining the optimal Bouguer density for a gravity data set:implications for the isostatic setting of the Mediterranean Sea. Geophysical Journal International, 169(2): 380-388. DOI:10.1111/j.1365-246X.2007.03340.x |
Vajk R. 1956. Bouguer corrections with varying surface density. Geophysics, 21(4): 1004-1020. DOI:10.1190/1.1438292 |
Wang Q S. 2003. Gravitology (in Chinese). Beijing: Seismological Press.
|
Zeng H L. 2005. Gravity Field and Gravity Exploration (in Chinese). Beijing: Geological Publishing House.
|
崔作舟, 卢德源, 陈纪平, 等. 1987. 攀西地区的深部地壳结构与构造. 地球物理学报, 30(6): 566-580. DOI:10.3321/j.issn:0001-5733.1987.06.003 |
葛碧如, 杨科佑. 1990. 攀西地区中、新生代构造特征. 地球物理学报, 33(1): 64-69, 125-126. DOI:10.3321/j.issn:0001-5733.1990.01.007 |
缪以琨, 任祖德, 陈小源, 等. 1986. 四川攀西裂谷的构造环境及构造发展. 地球科学-武汉地质学院学报, 11(6): 631-637. |
石磊, 楼海, 王谦身, 等. 2015. 攀西地区重力场特征及地壳密度结构. 地球物理学报, 58(7): 2402-2412. DOI:10.6038/cjg20150717 |
王谦身. 2003. 重力学. 北京: 地震出版社.
|
曾华霖. 2005. 重力场与重力勘探. 北京: 地质出版社.
|