0 引言
重力梯度数据包含丰富的地质信息，与其相关的数据处理及解释方法发展迅速[1]。将协克里金引入到地球物理反演中处理重力梯度数据，具有很多优势：抗噪性好，当观测数据噪声较大时仍能够得到稳定的反演结果；能够以协方差的形式方便地融合先验地质信息；反演结果对变程、块金常数、基台值的选择不敏感[2]。由于这些优势，协克里金得到了广泛应用[3-5]。协克里金反演方法的局限性也很明显，常用的协方差矩阵具有平稳的特征，使整个研究区域中密度异常分布趋势相同；非平稳协方差的引入解决了这个问题[6]。岩脉倾斜角度的先验信息对反演结果影响非常大，前人在模拟试验过程中默认倾斜角度已知，以先验信息的形式将正确的岩脉倾角引入[2, 7-8]，这意味着当有关岩脉倾角的先验信息不充足或者错误时，很难通过协克里金法得到正确的密度分布。
本文致力于研究岩脉倾角信息未知时协克里金反演方法的使用。先设定4种具有代表性的模型，假设岩脉模型倾角未知，通过扫描的方法估计岩脉倾角；再将估计的岩脉角度以先验信息的形式加入到协克里金反演中；最后将扫描倾角的方法用于文顿盐丘地区的航空重力梯度数据，验证方法的有效性。
岩脉的倾斜角要通过参数协方差加入到协克里金反演中，计算参数协方差常用的方法有三种：一是当先验信息充足时，Chasseriau等[8]通过大量已获得的密度点直接计算协方差；二是在缺少先验信息的情况下，Asli等[3]利用V-V图的方法来拟合估计协方差矩阵；三是根据平稳高斯过程的条件，相关函数与方差相乘得到参数协方差函数[6, 9-10]。我们选择第三种方法获得协方差矩阵，它可以通过改变旋转矩阵引入不同的倾角，方便地完成角度扫描。
1 方法原理 1.1 协克里金法计算由地下密度异常体引起的观测重力或重力梯度异常的常用方法是基于牛顿定律，将3D模型空间划分成具有均一密度值的简单几何体：
(1) 式中：N维向量dN为重力梯度异常观测数据；M维向量ρM为密度；M×N维矩阵AN×M为几何因子。由于观测数据与密度是线性相关的，根据地质统计学知识，它们的自协方差矩阵也是线性相关的，其交叉协方差矩阵与密度自协方差矩阵也是线性相关的：
(2) 
(3) 式中：N×N维矩阵Cdd为观测数据协方差；M×M维矩阵Cρρ为密度协方差；N×N维矩阵C0为观测数据误差协方差矩阵；N×M维矩阵Cdρ为观测数据与密度之间的交叉协方差矩阵。假设观测数据的平均值E(d) 和密度参数的平均值E(ρ) 均为0，则估计方差位于E((ρ-ρ*)(ρ-ρ*)T) 的主对角线上[7]：
(4) 式中：ρ和ρ*分别为真实密度和估计密度；Λ为加权系数矩阵。为了使估计方差最小，令式 (4) 偏导等于0，即
(5) 则
(6) 从式 (6) 中我们能够得到最佳加权系数矩阵。从而，可以得到估计密度：
(7) 在三维情况下，对于任意平稳的高斯过程，两点的相关函数为
(8) 式中：si、sj为位置向量；Σ为3×3的核矩阵。
(9) 式中：R为旋转矩阵；Ψ为由椭圆变程决定的对角线矩阵。在三维空间内，设以岩脉倾向为主轴的变程椭球坐标系仅绕第三轴旋转α角度即能与笛卡尔坐标系重合，则
(10) 
(11) 式中：a1为沿岩脉方向的变程；a3为沿z方向的变程；a2为沿垂直于前两个轴的方向的变程。
将式 (8) 与方差σ2相乘，即得到参数协方差函数Cρρ。
通过改变旋转角度，完成角度扫描，容易获得不同岩脉倾角对应的协方差矩阵，得到反演结果。统计预测数据与测量数据残差的标准差，剔除不合理异常后，旋转角度与实际岩脉倾角相同时，标准差最小。当剩余密度为正值时，将观测数据中的负异常值视为不合理异常；当剩余密度为负值时，将观测数据中的正异常值视为不合理异常。
2 模型试验将地下三维密度空间剖分为21×19×10＝3 990个大小为100 m×100 m×100 m的规则立方体单元，21×19＝399个地面测点位于每个密度单元的中心。在理论重力梯度异常Tzz(重力梯度垂直分量) 中加入5%的高斯噪声。为了验证方法的普遍适应能力，选择的4个岩脉模型的倾角分别为0°，45°，90°，135°，岩脉宽度为300 m，相应的模型主视图如图 1所示。
|
| h.深度；y=1 000 m。a.0°倾角岩脉模型；b. 45°倾角岩脉模型；c. 90°倾角岩脉模型；d. 135°倾角岩脉模型。 图 1 3D理论模型密度分布中心切片 Figure 1 Single section of 3D synthetic model |
|
|
为了对比，令4个岩脉模型用于协克里金反演的决定密度协方差的参数都相同：σ2＝500 (kg/m3)2，a1＝500 m，a2＝200 m，a3＝200 m。设定扫描角度θ分别为0°、15°、30°、45°、60°、75°、90°、105°、120°、135°、150°、165°、180°(扫描角度180°与扫描角度0°相同)，多次反演。
根据协克里金反演的密度分布预测Tzz，首先求Tzz预测值与Tzz测量值拟合残差的标准差σ1；然后根据先验信息，剔除密度分布中的负值，再求一组Tzz预测值与Tzz测量值拟合残差的标准差σ2。统计两种情况下标准差随扫描角度的变化情况，如表 1所示，对应的趋势如图 2所示。
| θ/(°) | 0°倾角岩脉模型 | 45°倾角岩脉模型 | 90°倾角岩脉模型 | 135°倾角岩脉模型 | |||||||
| σ1/E | σ2/E | σ1/E | σ2/E | σ1/E | σ2/E | σ1/E | σ2/E | ||||
| 0 | 4.57 | 5.88 | 4.11 | 4.57 | 5.81 | 6.06 | 4.18 | 4.68 | |||
| 15 | 4.48 | 6.47 | 3.92 | 4.62 | 5.25 | 6.28 | 3.89 | 5.03 | |||
| 30 | 4.26 | 7.04 | 3.46 | 4.66 | 4.47 | 6.11 | 3.43 | 5.51 | |||
| 45 | 4.01 | 7.39 | 3.09 | 4.56 | 4.10 | 5.71 | 3.13 | 5.73 | |||
| 60 | 3.94 | 7.62 | 2.96 | 4.65 | 3.91 | 5.36 | 2.99 | 5.76 | |||
| 75 | 3.92 | 7.75 | 2.91 | 4.84 | 3.81 | 5.09 | 2.93 | 5.61 | |||
| 90 | 3.91 | 7.75 | 2.89 | 5.06 | 3.78 | 4.97 | 2.91 | 5.35 | |||
| 105 | 3.92 | 7.63 | 2.91 | 5.12 | 3.81 | 5.05 | 2.93 | 5.02 | |||
| 120 | 3.95 | 7.44 | 2.96 | 5.33 | 3.91 | 5.26 | 2.99 | 4.78 | |||
| 135 | 4.01 | 7.14 | 3.10 | 5.31 | 4.09 | 5.57 | 3.12 | 4.75 | |||
| 150 | 4.25 | 6.81 | 3.44 | 5.18 | 4.47 | 5.95 | 3.45 | 4.92 | |||
| 165 | 4.47 | 6.29 | 3.89 | 4.78 | 5.25 | 6.19 | 3.91 | 4.83 | |||
| 180(0) | 4.57 | 5.88 | 4.11 | 4.57 | 5.81 | 6.06 | 4.18 | 4.68 | |||
| 注：E (厄缶) 为非法定计量单位，1 E=10-9 s-2，下同。 | |||||||||||
|
| a. 0°倾角岩脉模型；b. 45°倾角岩脉模型；c. 90°倾角岩脉模型；d. 135°倾角岩脉模型。 图 2 标准差随扫描角度的变化 Figure 2 Standard deviations changes with scan angles |
|
|
由图 2可知，4个岩脉模型的σ1整体变化趋势一致：σ1整体较小；当θ为0°(180°) 时，σ1最大；当0°≤θ＜90°时，σ1随θ单调递减；当θ为90°时，σ1最小；当90°≤θ＜180°时，σ1随θ单调递增。这种统计特征，可作为判断统计有效性的标志。
σ2的变化特征各不相同 (图 2)。在0°≤θ＜180°范围内：0°倾角岩脉模型只有当θ为0°时，σ2出现一个极小值；45°倾角模型分别在θ为0°和45°时，σ2出现极小值；90°倾角岩脉模型分别在θ为0°和90°时，σ2出现极小值；135°倾角岩脉模型分别在θ为0°和135°时，σ2出现极小值。因此，无论岩脉真实倾角多大，当θ＝0°时都会出现极小值，另一极小值对应的θ为本方法对真实倾角的估值；若除θ＝0°处的极小值外没有其他极小值点，说明两个极小值点在θ＝0°处重合，则本方法的估计倾角为0°。
将估计的岩脉倾角以先验信息的形式引入到协克里金反演过程中，剔除不合理负值，反演得到的密度分布 (图 3) 与真实模型吻合。
鉴于方法对模型试验的成功，我们将方法运用于实际测量数据，进一步验证方法可靠性。Bell Geospace公司在文顿盐丘地区测量了航空重力梯度全张量数据。根据前人研究结果[11-12]，盐丘密度与围岩密度均为2 200 kg/m3，盖层密度约为2 750 kg/m3，盖层与岩丘和围岩的密度差是引起观测异常的主要原因。对数据进行地形校正时，选择校正密度为2 200 kg/m3，盖层岩石的剩余密度为550 kg/m3。我们截取测区中覆盖文顿盐丘的一部分重力梯度垂直分量Tzz数据进行反演计算，测点网格为21×21，如图 4所示。
|
| 图 4 文顿盐丘地区实测全张量Tzz数据 Figure 4 Tzz component data of airborne FTG data measured over the Vinton salt dome |
|
|
将地下三维空间划分成21×21×15个大小为200 m×200 m×100 m的规则块体。决定密度协方差的参数σ2=20 000 (kg/m3)2，a1=800 m, a2=400 m, a3=800 m。设定θ为0°、15°、30°、45°、60°、75°、90°、105°、120°、135°、150°、165°，统计每个角度反演结果的预测数据与观测数据残差的标准差，结果如表 2所示，对应的趋势如图 5所示。原始结果预测数据与观测数据残差标准差的统计规律与模型试验的统计规律相同，证明此次统计结果有效。根据剔除负值后的统计结果，标准差极小值对应的角度为0°和90°，因此，确定倾斜角度为90°。
| θ/(°) | σ1/E | σ2/E |
| 0 | 4.35 | 9.49 |
| 15 | 4.00 | 9.27 |
| 30 | 3.64 | 9.20 |
| 45 | 3.38 | 9.15 |
| 60 | 3.22 | 9.12 |
| 75 | 3.14 | 9.10 |
| 90 | 3.01 | 9.08 |
| 105 | 3.13 | 9.11 |
| 120 | 3.20 | 9.15 |
| 135 | 3.36 | 9.21 |
| 150 | 3.62 | 9.30 |
| 165 | 4.00 | 9.37 |
| 180(0) | 4.35 | 9.49 |
|
| 图 5 文顿盐丘地区实际数据标准差随扫描角度的变化 Figure 5 Standard deviations changes with scan angles for the real observed data over Vinton salt dome |
|
|
保持决定密度协方差的参数不变，将倾斜角度90°作为先验信息加入协克里金理论中再反演，得到的密度分布如图 6所示。密度异常体的中心位置与已知地质信息一致，反演模型异常体的最大密度值为250 kg/m3，比已知盖层剩余密度550 kg/m3约小一半；这是由协克里金法本身的反演理论决定的，它属于光滑反演的一种，没有加入聚焦信息，故反演结果比较发散。
|
| 图 6 极小值标准差对应的90°倾角作为先验信息反演的文顿盐丘密度分布 Figure 6 Recovered Vinton salt dome model based on the prior information of 90° dip angle corresponding to the minimum standard deviation |
|
|
由密度分布得到的预测数据如图 7所示，预测数据与观测数据 (图 4) 吻合。
|
| 图 7 文顿盐丘地区全张量重力梯度垂直分量反演模型预测数据 Figure 7 Tzz component data of airborne FTG data predicted from the recovered model over the Vinton salt dome |
|
|
为了进一步验证方法的有效性，选择一个非极小标准差对应的角度45°，作为倾斜角度引入到协方差中，协克里金反演所得密度分布如图 8所示，与已知的地质信息不符。
|
| a. y=1 050 m的竖直切片；b. x=1 050 m的竖直切片；c. h=300 m的水平切片。 图 8 非极小值标准差对应的45°倾角作为先验信息反演的盐丘密度分布 Figure 8 Recovered Vinton salt dome model based on the prior information of 45° dip angle corresponding to the optional standard deviation |
|
|
1) 本文通过角度扫描，利用协克里金法反演重力梯度垂直分量Tzz，统计预测数据与观测数据残差的标准差，能够获得有关角度的信息。岩脉真实角度非0°时，当0°≤θ < 180°，会存在两个标准差极小值，分别对应0°和真实倾角角度；岩脉真实角度为0°时，当0°≤θ < 180°，只存在一个极小值，对应0°。
2) 对于4种代表性模型，本文方法都能够得到正确的倾角信息。再将倾角信息加入协克里金反演，扩宽了协克里金方法的使用范围。
3) 将确定倾角的角度扫描方法应用于文顿盐丘实测重力梯度数据，也证实了正确地估计倾角能够得到更好的反演结果。
致谢: Bell Geospace公司允许笔者使用其于2008年在文顿盐丘地区测量的航空重力梯度全张量数据，特此致谢。| [1] | 袁园, 黄大年, 余青露, 等. 全张量重力梯度数据误差分析及补偿[J]. 吉林大学学报 (地球科学版), 2014, 44(3): 1003-1011. Yuan Yuan, Huang Danian, Yu Qinglu, et al. Error Analysis and Compensation of Full Tensor Gravity Gradiometer Measurements[J]. Journal of Jilin University (Earth Science Edition), 2014, 44(3): 1003-1011. |
| [2] | Geng M, Huang D, Yang Q, et al. 3D Inversion of Air-borne Gravity-Gradiometry Data Using Cokriging[J]. Geophysics, 2014, 79(4): G37-G47. DOI:10.1190/geo2013-0393.1 |
| [3] | Asli M, Marcotte D, Chouteau M. Direct Inversion of Gravity Data by Cokriging[C]//Kleingeld W, Krige D. Proceedings of the 6th International Geostatistics Congress. 2000:64-73. |
| [4] | Gloaguen E, Marcotte D, Chouteau M, et al. Borehole Radar Velocity Inversion Using Cokriging and Cosimulation[J]. Journal of Applied Geophysics, 2005, 57(4): 242-259. DOI:10.1016/j.jappgeo.2005.01.001 |
| [5] | Shamsipour P, Chouteau M, Marcotte D. 3D Stochastic Inversion of Magnetic Data[J]. Journal of Applied Geophysics, 2011, 73(4): 336-347. DOI:10.1016/j.jappgeo.2011.02.005 |
| [6] | Shamsipour P, Marcotte D, Chouteau M, Rivest M, et al. 3D Stochastic Gravity Inversion Using Nonstationary Covariances[J]. Geophysics, 2013, 78(2): G15-G24. DOI:10.1190/geo2012-0122.1 |
| [7] | Shamsipour P, Marcotte D, Chouteau M, et al. 3D Stochastic Inversion of Gravity Data Using Cokriging and Cosimulation[J]. Geophysics, 2010, 75(1): I1-I10. DOI:10.1190/1.3295745 |
| [8] | Chasseriau P, Chouteau M. 3D Gravity Inversion Using a Model of Parameter Covariance[J]. Journal of Applied Geophysics, 2003, 52(1): 59-74. DOI:10.1016/S0926-9851(02)00240-9 |
| [9] | Higdon D, Swall J, Kern J. Non-Stationary Spatial Modeling[J]. Bayesian statistics, 1999, 6(1): 761-768. |
| [10] | Paciorek C J, Schervish M J. Spatial Modelling Using a New Class of Nonstationary Covariance Functions[J]. Environmetrics, 2006, 17(5): 483-506. DOI:10.1002/(ISSN)1099-095X |
| [11] | Hou Z, Wei X, Huang D, et al. Full Tensor Gravity Gradiometry Data Inversion:Performance Analysis of Parallel Computing Algorithms[J]. Applied Geophysics, 2015, 12(3): 292-302. DOI:10.1007/s11770-015-0495-z |
| [12] | Ennen C. Mapping Gas-Charged Fault Blocks Around the Vinton Salt Dome, Louisiana Using Gravity Gradiometry Data[D]. Houston:University of Houston, 2012. |