徐紫阳，朱自强

（中南大学地球科学与信息物理学院，湖南长沙 410083）

重磁勘探是一种重要勘探地球物理方法，在矿产资源勘察和地质调查中有着重要的作用[1-3]，三维物性反演是重磁资料处理和解释中至关重要的一个环节。但是，由于观测数据量的局限性以及不可避免的噪声干扰，使得反演结果常常是病态的。为了克服这一问题，众多学者通过引入各种约束条件来降低反演的多解性，使得病态问题趋于良态。Tikhonov 等[4]提出在反演过程中通过加入正则化项来获得稳定的解；Li 等[5]针对位场反演过程中在深部方向上的分辨率低的问题，通过引入深度加权函数来消除这种“趋肤效应”；兰学毅等[6]通过收集不同地质先验信息进行基于先验信息约束的重磁三维反演，梁生贤[7]利用反演迭代过程中的拟合残差计算互相关系数，在此基础上，采用互相关系数与深度加权对重力反演进行约束求解。在进行地球物理数据处理和解释的过程中，不同的地球物理方法所对应的物性参数不同，误差和噪声来源也不同；相较单一物性反演而言，多种地球物理方法的联合反演能够更好的综合应用不同类型的地球物理数据，增强反演稳定性和精确度[8-11]。因此，多种勘探地球物理方法的联合反演也越来越受到重视。

交叉梯度约束的结构耦合方法自Gallardo 和Meju[12]提出以来，联合反演应用效果良好，迅速得到了广泛的应用。Emilia 和Fregoso 等[13]将交叉梯度函数引入到重力观测数据和磁法观测数据联合反演工作中，数值计算结果表明其在横向和深度分辨率方面得到改进。Wang 等[14]于球坐标下构造交叉梯度算子，实现了大尺度重磁联合反演。修春晓等[15]研究地质体的结构特性，实现重磁联合反演，保证了反演结果与先验构造信息具有较强的一致性。侯宇健等[16]利用交叉梯度约束实现了三维极化率/电阻率的联合反演；闫政文等[17]对多种物性参数进行了交叉梯度联合反演工作。

本文通过引入交叉梯度算子，基于重磁多约束单一反演开展三维联合反演工作，反演过程使用预条件共轭梯度法进行求解。通过建立不同三维地质理论模型，分别开展单一反演和联合反演研究工作，并且将单一反演结果和联合反演结果进行了对比分析。

1 重磁正反演理论

1.1 正则化约束

在重力和磁法勘探过程中，通过由地下密度或者磁性异常体在地表所引起的重力场和磁场的变化来计算获得地下地质异常密度和磁化率分布情况，从而确定地下地质异常体的赋存情况。

重磁正反演问题可以用下式来表示：

其中：d 为长度为Nd=nx×ny的观测数据向量，可以是磁法或重力数据，nx×ny分别表示测区范围内的测线与测点数；m 为长度为Nm=nx×ny×nz的模型参数向量，nz是地下网格划分层数；G 是反演中灵敏度矩阵或核矩阵。正演问题是已知G 和m 的值，求取观测数据d的过程；反演问题则是通过地表观测数据d，推导计算地下模型参数m。

地表观测数据是有限的离散的数据，且采集过程存在观测误差和噪声干扰，同时在反演过程中，地下场源参数个数远远大于观测数据个数，这导致反演问题的结果具有多解性和非唯一性。为了提高反演精确度和稳定性，本文将Tikhonov 正则化思想引入到反演计算中进化求解：

其中：φ（m）为数据拟合泛函；S（m）为模型稳定泛函；α 代表正则化参数，也称阻尼参数，它通过调节数据拟合泛函与模型目标函数之间的平衡，以防止数据过度拟合或模型过度平滑。

其中数据拟合函数可表示为：

1.2 预条件共轭梯度法

共轭梯度法的基本思想是把共轭性与最速下降法相结合，利用已知点处的梯度构造一组共轭方向，沿着这个方向去搜索目标函数极小值点。根据共轭方向的性质，共轭梯度法具有二次终止性。但是，由于方程组本身是欠定的，模型参数规模较大，导致雅可比矩阵的条件数较大、求解过程不稳定。本文通过引入预条件因子，对迭代法进行改性，以改善反演收敛次数多，反演不稳定等问题。

假设观测数据为b，系数矩阵为A，模型参数为m，初始模型参数为m0。则对于共轭梯度法中的方程ATAm=ATb 可改进为

其中：S 为预条件因子，在理想条件下，S 是ATA逆的近似值，因此SATA 约等于单位矩阵。与原始对称正定矩阵ATA 相比，SATA 的特征值更加聚类。因此，减少了条件数，并加快了CG 方法的收敛速度。但是，以这种方式求解S 计算量较大，求解较困难。在实际应用中，用式（8）中的深度加权函数作为预条件因子进行计算。

2 基于交叉梯度法重磁联合反演

2.1 交叉梯度函数基本形式

交叉梯度函数最早由Gallardo 和Meju[12]提出，因其应用限制较小，且应用效果良好，迅速在地球物理学中得到广泛应用。两种物性参数的三维形式交叉梯度函数如下：

其中：mg和mm为联合反演中的物性参数。在重磁联合反演中，mg代表密度参数、mm代表磁化率参数。式（6）为上式的具体展开形式：

实测数据是离散的且地下网格剖分也是离散的。因此，在应用交叉梯度函数时，还需要对其进行离散化处理。本文采用泰勒展开处理，只保留一次项，简化可得式（7）：

式（7）在零点处展开，其形式可进一步简化为：

以中心差分代替微分，Bg可表示为如下形式：

Bm的表示式也可由式（8）形式得到。

2.2 联合反演目标函数构建

多约束条件的单一反演算法以及交叉梯度函数，进行重磁联合反演研究工作，构建其联合反演目标函数。

交叉梯度项可以表示为：

综合可得联合反演目标函数如下：

其中：Φg为重力反演项；Φm为磁法反演项；β 为重力反演项的加权因子；γ 为磁法反演项的加权因子；μ表示交叉梯度项的加权因子。

3 理论模型试算

本文采用MATLAB 语言开发了重力、磁法联合反演算法，为了检验算法的有效性，通过建立长方体模型进行数值模拟反演研究。

设置一个理论的地下三维物性模型（见图1），地下场源空间大小为200 m×200 m×100 m，将地下场源空间划分为独立的nx×ny×nz=20×20×20=8 000 个模型单元体，每个单元格规模为10 m×10 m×5 m。设置一个棱柱体异常模型，其中心坐标为（0，0，30）。异常体x 轴方向，即东向规模为80 m，y 轴方向，即北向长为80 m，棱柱体顶部埋深为20 m，底部埋深为40 m。地下物性异常体与围岩间的剩余密度为ρ=1 g/cm3，地磁场倾角50°，磁偏角0°，无剩磁存在，物性异常体与围岩间的相对磁化率设置为κ=1 SI。

图1 单一棱柱体模型示意图Fig.1 The single-prism model

根据所设置单一模型体，分别进行重力和磁法正演计算，得到异常体在地表所引起的重力异常响应和磁异常响应，用所得重磁异常数据分别开展重磁单一反演工作和重磁联合反演工作，从而得到地下介质密度和磁化率分布情况。

通过重磁两种方法的单一反演结果与联合反演结果进行对比（见图2），可以看出无论是对于密度反演结果还是磁化率反演结果，联合反演相比于单一反演，在边界的分辨能力上都具有较大的提升，对异常体形态大小的描述更加接近真实情况。

图2 单一模型单独反演与联合反演结果对比图Fig.2 The comparison between separate and joint inversion of the single prism model

4 结论

（1）本文对交叉梯度函数进行讨论，基于重磁单一反演，开展了基于预条件共轭梯度法的重磁交叉梯度三维联合反演工作，该算法收敛速度快，稳定性好，可较好的应用于三维反演计算。

（2）模型反演结果表明，交叉梯度联合反演算法较单一反演而言，可以挖掘不同物性信息之间的联系，更为准确的刻画出地下异常体位置信息和规模，降低反演问题的多解性，使反演结果更加收敛。