0 引言

地球物理反演是数据处理与解释的重要手段，如何解决反演结果多解性问题一直是研究热点，联合多种数据共同反演是降低多解性的一种重要手段(Gallardo and Meju, 2003; 于鹏等, 2007; Zhdanov et al., 2012; 高级和张海江, 2016; 殷长春等, 2018).处理重磁及其张量数据的反演方法主要是基于正则化理论的反演方法，包括光滑反演(Li and Oldenburg, 1996, 1998)、紧致反演(Last and Kubik, 1983)、聚焦反演(Zhdanov, 2002, 2009; 秦朋波和黄大年, 2016; 高秀鹤和黄大年, 2017)、聚类C均值反演(Sun and Li, 2015)、多元反演(Lin and Zhdanov, 2018)等.Myers (1982)从数学角度推导了协克里金理论，Asli等(2000)首次直接利用协克里金方法反演重力数据，并在缺少地质信息的情况下，提出用Ⅴ-Ⅴ图的方法来拟合变差函数，利用变差函数求出密度协方差函数，进而实现协克里金反演.从此，协克里金反演方法被广泛用于反演地球物理数据，包括井中雷达速度数据(Gloaguen et al., 2005)、重力数据(Shamsipour et al., 2010)、磁数据(Shamsipour et al., 2011)、以及重力梯度数据(Geng et al., 2014; 耿美霞等, 2016).Chasseriau和Chouteau (2003)给出了基于初步解释模型计算空间中任意方向变差函数的方法，并联合处理井中重力数据和地面重力数据，提高了反演精度和分辨率.前面提到的方法都是基于平稳、发散的协方差矩阵实现的协克里金反演，这导致了反演结果发散，Shamsipour (2013)将非平稳协方差函数(Higdon et al., 1999)引入到协克里金反演中处理重力数据，与基于常规协方差函数的协克里金方法对比，反演效果明显提高，但结果仍然发散.协克里金反演得到的图像发散问题限制了该方法的适用性，比如具有密度突变的地质构造，常规协克里金法无法得到清晰的接触面.本文致力于改善协克里金法反演结果发散问题，利用密度阈值对常规协方差函数进行修剪，得到非平稳、相对聚焦的协方差矩阵，迭代执行协克里金反演程序，最终改善了协克里金方法的反演效果，得到了相对聚焦的反演图像.

1 反演原理及方法 1.1 协克里金反演原理

将地下3D模型空间划分成具有均一密度值的M个规则立方体，地面观测点数为N，基于牛顿定律计算由地下密度异常体引起的重力或重力梯度异常：

(1)

其中，列向量d是地面观测数据，即重力或重力梯度异常观测数据，其维数是观测点数N；列向量ρ是密度，其维数为地下空间剖分的立方体数M；矩阵A表示与观测数据对应的正演算子，是一个N×M维矩阵.由于观测数据与密度是线性相关的，根据地质统计学知识，其自协方差矩阵也是线性相关的，其交叉协方差矩阵与密度自协方差矩阵也是线性相关的(Myers, 1982; Asli et al., 2000)：

(2)

(3)

式中，N×N维矩阵C_dd为观测数据协方差；M×M维矩阵C_ρρ为密度协方差；N×N维矩阵C₀为观测数据误差协方差；N×M维矩阵C_dρ为观测数据与密度之间的交叉协方差.根据地质统计学知识，假设均值E(d)=E(ρ)=0，计算密度向量ρ的估计方差如下(Shamsipour et al., 2010)：

(4)

每个块体的密度估计方差位于(4)式主对角线上. ρ和ρ^*分别为真实密度向量和估计密度向量；Λ为加权系数矩阵.为了使估计方差最小，令(4)式偏导等于0：

(5)

则有

(6)

由(6)式求解加权系数矩阵Λ，将问题转化成求解目标函数最优解问题，如下：

(7)

α为阻尼因子，其阻尼最小二乘解为

(8)

由(8)式得到最佳加权系数矩阵Λ.从而，可以计算估计密度

(9)

按照上述步骤，只要能够获得密度协方差矩阵C_ρρ，即可实现协克里金反演计算.因此，密度协方差矩阵的准确程度决定了最终反演结果的优劣.

另外，重磁位场数据反演得到的结果具有严重的趋肤效应(Fedi et al., 2005)，为了解决这个问题，深度加权函数广泛应用于正则化反演中.基于地质统计学的协克里金法的反演结果同样存在趋肤效应，为了克服趋肤效应，本文引入基于灵敏度矩阵的深度加权函数W (Zhdanov, 2002)，其人为干扰因素少，是一个对角矩阵，对角元素为

(10)

在协克里金反演中加入深度加权函数的方法比较简单，形式上只要改变C_ρρ即可

(11)

1.2 阈值修剪获得密度协方差矩阵及迭代反演方法

概率论中，期望值分别为E(X)与E(Y)的两个随机变量X和Y之间的协方差定义为

(12)

由于存在假设条件均值E(d)=E(ρ)=0，计算理论密度协方差矩阵：

(13)

由公式(13)可知，若真实模型密度ρ_r已知，则理论密度协方差为

(14)

由(14)式，我们能够获得协方差矩阵的性质：若某一块体的密度ρ_i=0，则与其相关的所有协方差都为0，即

(15)

通过以上性质可以看出，理论模型密度集中分布时，其理论密度协方差矩阵也应是集中分布的.实际工作中，我们的目的是反演密度分布，不可能通过(14)式获得理论密度协方差矩阵，因此，根据式(15)表示的理论协方差的性质，我们设置密度阈值修剪常规协方差矩阵，以改善常规协方差矩阵的性质，并配合迭代反演方法进一步发展协克里金反演方法，具体反演流程如图 1所示.由反演流程图可以看到，经过阈值ρ_thr多次修剪的协方差矩阵打破了常规协方差矩阵平稳、发散的性质；每次调整阈值后，上一步得到的C_ρρ^k会被进一步修剪并应用于协克里金反演，这可以充分利用前面的计算结果并得到稳定的密度分布；对于密度阈值的选取，首先，选择一个远远小于理论异常值的初始密度阈值，为了保证反演结果稳定，阈值增量Δρ_thr也远远小于理论异常值，在迭代反演过程中，逐渐增大阈值，直到所得的密度模型的最大值达到设置的密度上限ρ_max，输出反演结果.

1.3 联合反演方法

本文联合重力数据g与重力梯度数据T_zz共同反演地下密度分布，由于两种数据对应同一种物性参数，即密度ρ，且两种数据与密度都是线性相关，这种情况联合反演比较容易实现，只要将两种数据同时作为观测数据d带入反演程序即可

(16)

系数λ₁、λ₂的大小取决于相应观测数据的幅值、观测数据中噪声水平等因素. λ₁、λ₂的大小决定两种数据在反演中所起作用的大小：系数越大，相应数据在反演中所起的作用越大，反之亦然.当某一系数为0，说明相应的观测数据在反演中不起作用，此时变成另一种数据的单独反演.与上述观测数据对应，联合反演的正演算子A写成如下形式：

(17)

其中，系数λ₁、λ₂取值与公式(16)相同，A_g表示与重力数据g对应的正演算子，A_Tzz表示与重力梯度数据T_zz对应的正演算子.

2 模型测试

我们将基于常规协方差函数、理论协方差函数(公式(14))以及阈值约束改进的协方差函数反演两种模型(1)双棱柱体模型；(2)倾斜岩脉模型，并分析反演结果.两个模型的地下空间网格剖分和观测网格相同，都是将地下三维密度空间划分为19×21×15=5985个大小为100 m×100 m×100 m的立方体单元，19×21个地面观测数据位于每个立方体单元中心.通过正演计算获得理论的重力异常g和重力梯度异常T_zz，并在理论异常中加入标准差等于幅值的5%的高斯噪声.

首先设置双棱柱体模型，如图 2a，两个棱柱体大小均为500 m×400 m×500 m，密度差均为1 g·cm^-3，顶部埋深均为300 m，此模型地面观测重力异常和重力梯度数据异常T_zz中均包含5%高斯噪声，分别如图 2b和2c.首先，为了估计密度协方差矩阵，我们建立初步解释模型，基于初步解释模型计算沿坐标轴x、y、z方向的变差函数

(18)

此处，h为滞后距，N为数据对的个数，ρ(x_i)和ρ(x_i+h)为一个数据对，ρ(y_i)和ρ(y_i+h)为一个数据对，ρ(z_i)和ρ(z_i+h)为一个数据对，V_x(h)为x方向滞后距为h的变差函数值，V_y(h)为y方向滞后距为h的变差函数值，V_z(h)为z方向滞后距为h的变差函数值，从而可以画出x、y、z方向变差函数与滞后距h的关系图，即三个方向的变差函数图.选择高斯模型分别拟合上述三个方向的变差函数图，高斯模型的形式为

(19)

γ(h)为高斯模型函数，C₀为块金常数，C为基台值，h为滞后距，a为变程.拟合x、y、z方向的变差函数图的三个高斯模型的块金常数平均值C₀=2000 (kg·m^-3)²，平均基台值C=25000 (kg·m^-3)²，x、y、z方向的变程分别为a_x=400 m，a_y=500 m，a_z=500 m，则通过(20)式可以获得任意方向的变程(Chasseriau and Chouteau, 2003)

(20)

θ表示任意方向的方向向量在x-y平面的投影与x轴的夹角；α表示任意方向的方向向量与z轴的夹角.用高斯模型拟合任意方向的变差函数

(21)

通过协方差函数与变差函数、方差之间的关系，我们能够求出任意方向两个密度块体的协方差，即

(22)

其中，方差σ²=C.在整个模型空间通过(22)式扫描，能够求出常规密度协方差矩阵C_ρρ，协克里金反演方法能够顺利实现.

双棱柱体模型重力数据g和重力梯度数据T_zz中所含噪声水平相同，即两种数据的可信度相同，但重力数据的幅值远小于重力梯度数据，为了联合反演，我们令公式(13)、(14)中的系数λ₁、λ₂分别取100、1.通过常规协克里金方法反演得到的密度分布如图 3a，能够识别出密度异常体浅部的大致位置，但密度值明显偏小，反演出的最大密度值约为实际密度的一半，且异常体底部发散严重，两个异常体之间的部分甚至连成一片，无法区分它们的边界.为了验证理论协方差函数的正确性，我们假设能够通过真实的密度值计算出理论协方差矩阵，并用于协克里金反演，结果如图 3b所示，能够反演出密度异常的准确位置和理论异常值，由于观测数据中噪声的影响使得反演的密度值略小，证明了理论协方差的可靠性.在实际反演工作中，我们很难获得理论协方差函数，通过调整密度阈值修剪密度协方差矩阵，逐渐改进密度协方差的性质，迭代完成协克里金反演.设置阈值初始值ρ_thr=0.01 g·cm^-3，阈值改变量Δρ_thr=+0.01 g·cm^-3，直到阈值ρ_thr达到0.3 g·cm^-3时，反演的密度最大值达到密度上限1 g·cm^-3，停止迭代反演得到的结果如图 3c，异常体中心位置识别准确，边界清晰，反演结果收敛集中，尤其在异常体底部，与常规方法相比分辨率明显改善，异常体的密度值也更接近真实值.

本文从模型拟合程度和数据拟合程度两个方面来定量分析反演效果.如图 3中的反演模型(a)(b)(c)与理论模型的均方根误差表明模型拟合程度；模型预测数据g和T_zz与观测数据g和T_zz残差的标准差代表数据拟合程度，定量分析结果如表 1所示.相对于常规算法，基于阈值修剪的迭代算法反演得到的均方根误差更小(0.0992 g·cm^-3 < 0.178 g·cm^-3)，定量说明改进方法反演得到的模型更准确.通过数据拟合差可以看到，原始方法反演得到的模型产生的预测数据与实际观测数据残差的标准差比实际噪声标准差小很多(1.1523 mGal < 1.3641 mGal; 3.1128 E < 3.9103 E)，证明常规方法存在过拟合现象，而改进方法对应的标准差更接近实际噪声的标准差(1.4800 mGal≈1.3641 mGal; 3.8719 E≈3.9103 E)，这说明改进方法更合理.

第二个测试模型为倾斜岩脉模型，模型密度差为1 g·cm^-3，顶面埋深为300 m，如图 4a所示，模型正演的重力异常g和重力梯度异常T_zz如图 4b、4c，其中均包含5%高斯噪声.用与模型一相同的方法计算密度协方差矩阵，平均块金常数C₀=1000 (kg·m^-3)²，平均基台值C=15000 (kg·m^-3)²，沿岩脉倾斜方向的变程为1400 m，沿y轴方向的变程为200 m，沿垂直于前两个轴的方向的变程为300 m.

联合重力与重力梯度T_zz数据共同进行常规协克里金反演，系数λ₁、λ₂分别取100、1，结果如图 5所示，常规方法反演结果能够大致识别出倾斜岩脉的位置和倾向，最大密度值远小于真实密度，深部的密度分布发散.由理论协方差同样能够反演出完全正确的位置再一次说明理论密度协方差的可靠性.基于阈值修剪的方法迭代执行协克里金反演，设置初始阈值ρ_thr=0.01 g·cm^-3，阈值改变量Δρ_thr=+0.01 g·cm^-3，逐步迭代反演，直到阈值ρ_thr达到0.29 g·cm^-3，反演结果的最大密度值达到密度上限1 g·cm^-3，反演的密度分布能够准确识别出异常体位置和倾角，图像聚焦收敛，这体现了迭代反演的优势以及反演相对复杂模型体的能力.定量分析反演结果如表 2所示，从模型拟合项可以看出，基于阈值修剪的改进方法获得的模型与理论模型匹配更好，从数据拟合差项可以看出，改进方法的预测数据与观测数据残差标准差与噪声标准差匹配更好，这与模型一的结论是一致的.

3 实际数据应用

文顿盐丘位于美国路易斯安那州西南部，2008年美国Bell Geospace公司在此处进行了航空重力梯度测量，飞行高度约为80 m，测区面积为196.2 km²，共53条测线，两侧测线间距250m，中间位置加密测线，所以中间位置测线间距为125 m.另外还有17条联络测线，联络测线间距为1000 m.我们截取与文顿盐丘岩盖的位置对应的4 km×4 km区域内的数据进行网格化，网格大小为100 m×100 m，共41×41个观测数据.同时，将地下空间划分为41×41×20个长方体，每个长方体大小为100 m×100 m×50 m.与文顿盐丘有关的地质资料(Coker, 2007; Ennen, 2012; Gao et al., 2017)显示，盐丘与围岩密度基本相等，大概为2.2 g·cm^-3，因此，我们用2.2 g·cm^-3给重力梯度数据做地形改正，公开可用的重力数据为布格异常，对布格重力异常做去背景场处理，得到的重力梯度数据和重力数据作为待反演数据，如图 6所示.引起重力及重力梯度T_zz异常的地下密度异常体是盐丘顶部沉积的岩盖，密度约2.75 g·cm^-3，故岩盖异常体的密度差为0.55 g·cm^-3.

首先建立初步解释模型，高斯模型拟合三个方向的变差函数，三个方向高斯模型平均块金常数C₀=3000 (kg·m^-3)²，平均基台值C=10000(kg·m^-3)²，三个方向变程a_x=1400 m，a_y=1400 m，a_z=400 m.利用基于阈值约束的迭代协克里金法联合反演重力和重力梯度数据，设置初始密度阈值ρ_thr=0.01 g·cm^-3，阈值改变量Δρ_thr=+0.005 g·cm^-3，当密度阈值增长到0.210 g·cm^-3时，反演密度分布的最大值达到密度上限0.55 g·cm^-3，停止迭代计算得到反演结果切片如图 7所示，岩盖异常体的密度聚焦分布，边界清晰，中心埋深250 m，岩盖厚度约为200 m，东西方向延伸约为1500 m，南北方向埋深约为1300 m，反演结果与已知的地质信息相符(Ennen, 2012).为了对比效果，利用常规协克里金方法的反演结果与本文改进方法的反演结果做对比，常规协克里金方法反演结果切片如图 8所示，能够识别出岩盖异常体的水平位置，但垂直方向分辨率差，恢复的密度异常值远远小于岩盖的已知密度异常.通过对比图 7与图 8的反演结果，能够看出，改进的迭代协克里金反演方法的优越性.

4 结论

常规方法估计的密度协方差矩阵平稳、发散，本文基于一致性假设，推导的理论密度协方差函数表明，当真实模型密度分布聚焦时，密度协方差矩阵是紧致、非平稳的.实际工作中，几乎无法获得理论密度协方差矩阵，因此，本文提出利用密度阈值对常规方法估计的协方差矩阵进行处理，改善了协方差矩阵的性质.将处理后的协方差矩阵带入到协克里金反演中迭代计算，提高了常规协克里金方法的反演效果.理论模型试验表明，迭代协克里金反演方法得到的图像聚焦分布，与常规协克里金方法相比，分辨率大大提高，这弥补了常规协克里金方法只能得到光滑、发散图像的问题.定量分析说明，与常规方法相比，基于阈值的迭代协克里金方法得到的密度模型与理论模型匹配更好；另外，迭代协克里金方法反演得到的模型预测数据与观测数据残差的标准差与噪声标准差更相近，而常规方法反演得到模型的预测数据与观测数据残差的标准差小于噪声标准差，说明常规方法存在过拟合现象.用本文方法联合反演文顿盐丘重力和重力梯度数据，得到的密度异常分布与已知地质信息匹配较好.

致谢  感谢Bell Geospace公司提供文顿盐丘地区的实测数据.