CN117725354A - 一种大数据量重力与重力梯度联合的快速正反演方法及系统 - Google Patents

Abstract

本发明公开了一种大数据量重力与重力梯度联合的快速正反演方法及系统，属于重力梯度测量技术领域，其包括正演过程和反演过程，正演过程包括基于BTTB结构矩阵的灵敏度矩阵压缩，及二维快速傅里叶变换的快速正演；反演是在高斯‑牛顿法基础上，采用预条件共轭梯度法求解增量方程，从而实现快速且高精度的重力与全张量重力梯度联合反演。本发明不需要牺牲精度来实现重力与重力梯度的快速正演，能节省更多的存储资源和计算时间；在实现高精度反演的同时兼顾了反演方法与快速正演的结合，实现了高精度大数据量重力与全张量重力梯度的快速联合反演。

Description

一种大数据量重力与重力梯度联合的快速正反演方法及系统

技术领域

本发明涉及重力梯度测量技术领域，尤其涉及一种大数据量重力与重力梯度联合的快速正反演方法及系统。

背景技术

近年来，随着重力梯度测量技术的发展，重力与重力梯度数据的联合反演方法得到了广泛的关注。由于重力和重力梯度数据所包含的地下密度异常体的信息不同，重力数据具有更多的低频信息，重力梯度数据则具有更多的高频信息，对浅部异常体边界的刻画更清晰。因此，重力与重力梯度联合反演相比重力或重力梯度单分量的反演具有更多的信息源，可以减少解模型的不确定性，这有助于提高模型的准确性和可靠性，对于异常体边界的刻画也更清晰。

对重力和重力梯度联合反演问题，前人已经进行了大量研究，证明了重力数据和重力梯度分量进行联合反演，可以得到更为收敛的反演结果。不同的重力和重力梯度分量组合方式，也会产生不同的反演效果，其中重力数据与全张量重力梯度进行联合反演，反演效果往往优于其他组合方式。虽然联合反演可以获得更好的反演结果，但同时也带来了计算资源需求上的挑战，数据量的增加意味着需要消耗更多内存和计算时间。现有技术目前具体存在以下缺点：

1、目前反演往往会利用完整的灵敏度矩阵来进行反演，极大地限制了反演的剖分数量；

2、一些快速正演方法如GPU并行正演和合并单元格等技术无法有效地应用到反演计算当中；

3、一些基于小波压缩等频率域方法需要牺牲精度为代价来提升计算速度；

4、一些数据降维和模型降维等技术都会牺牲计算精度。

发明内容

为解决海量数据联合反演计算内存大、时间长的问题，本发明提出一种快速算法。相较于采用小波压缩等频域方法以牺牲精度为代价提升计算速度，本发明提出的算法在保证精度的前提下实现了加速。

为了实现上述目的，本发明所采用的技术方案如下：一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，包括正演过程和反演过程，其中，正演过程包括以下步骤：

获取网格剖分参数和测点坐标数据，基于BTTB结构矩阵计算得到压缩灵敏度矩阵，并存储；

对模型密度矩阵进行零填充处理，并实现与压缩的灵敏度矩阵相同维度，通过二维快速傅里叶变换得到有效的正演数据；

所述反演过程包括以下步骤：

获取网格剖分参数和测点数据，计算压缩灵敏度矩阵，并将各分量的压缩灵敏度矩阵拼合成组合灵敏度矩阵；

将组合灵敏度矩阵用于预条件共轭梯度法中的增量方程求解，通过迭代计算输出反演结果。

进一步的，所述的获取网格剖分参数和测点坐标数据，基于BTTB结构矩阵计算得到灵敏度矩阵，压缩并存储的步骤，具体包括：

S1、获取网格剖分参数、模型单元密度向量、观测点坐标；

S2、计算BTTB结构矩阵所对应的相对坐标矩阵(Δx，Δy，△z)；

S3、根据相对坐标矩阵同时计算重力和全张量重力梯度分量的正演灵敏度矩阵；

S4、将压缩后的灵敏度矩阵进行二维快速傅里叶变换得到频率域的灵敏度矩阵，并存储。

进一步的，所述的对模型密度矩阵进行零填充处理，并实现与压缩的灵敏度矩阵相同维度，通过二维快速傅里叶变换得到有效的正演数据的步骤，具体包括：

S5、将模型单元密度向量按层排列为三维矩阵，并用零填充到与压缩灵敏度矩阵相同维度；

S6、对填充后的三维模型密度矩阵按层做二维快速傅里叶变换，得到频率域的模型密度矩阵；

S7、将频率域的压缩灵敏度矩阵和模型密度矩阵进行点乘操作，得到包含正演结果的频率域数据矩阵；

S8、正演结果的频率域数据矩阵经过傅里叶逆变换得到空间域数据矩阵，从空间域数据矩阵中提取有效正演数据。

进一步的，在步骤S7中，所述频率域数据矩阵的计算公式为：

其中，为每个深度层对应的压缩灵敏度矩阵，/>为每个深度层对应的模型密度矩阵，/>为正演观测数据，/>代表正向FFT运算，/>代表逆FFT运算，'·'表示Hadamard积，即同维度矩阵中元素对应相乘的操作。

进一步的，所述获取网格剖分参数和测点数据，计算压缩灵敏度矩阵，并将各分量的压缩灵敏度矩阵拼合成组合灵敏度矩阵，及将组合灵敏度矩阵用于预条件共轭梯度法中的增量方程求解，通过迭代计算输出反演结果的步骤，具体包括：

S9、获取网格剖分参数、重力和全张量重力梯度数据，计算相应的压缩灵敏度矩阵；

S10、将重力和全张量重力梯度的压缩灵敏度矩阵与所述正演数据拼合成联合反演的组合灵敏度矩阵；

S11、根据吉洪诺夫正则化方法，加入对数界限法约束物性上下限，得到反演迭代求解的目标函数；

S12、将总目标函数按高斯-牛顿法的原理转换为增量方程的形式，高斯-牛顿法包含一次正演矩阵乘法用快速正演替代；

S13、利用预条件共轭梯度法求解每一步高斯-牛顿法迭代的增量，预条件共轭梯度法每次迭代包含一次正演矩阵乘法和一次矩阵转置乘法；

S14、预条件共轭梯度法求得模型增量后，在高斯-牛顿法中用有限步长防止迭代不收敛；

S15、设定阈值和收敛条件，判断迭代收敛后输出反演结果。

进一步的，在步骤S11中，总的目标函数为：

式中，为数据目标函数；/>是模型目标函数，/>为正则化参数，最后一项是界限函数，/>是密度约束上下界限，/>是界限参数，表示界限函数在总目标函数所占的比重。

进一步的，在联合反演时数据目标函数为公式（3）：

其中，不同分量以如下形式组合：

，/>，/>

模型目标函数如公式（4）：

其中，表示重力和重力梯度分量的拼合灵敏度矩阵，/>为数据加权项，/>为拼合的观测数据向量，/>为模型加权矩阵，/>为参考模型。

进一步的，在步骤S12中，增量方程的转换过程为：在第次迭代时，把解模型表示为/>，利用当前解模型/>以及解模型变化方向/>构成新的解模型/>，表示为：

为了最小化总目标函数（2），将其改写为：

利用泰勒级数展开，在时，得到增量方程：

式中，，/>，/>，。

进一步的，在步骤S14中，采用有限步长法确定合适的步长：

其中，是最大搜索步长，由公式（14）计算：

式中，参数γ表征步长的减少量，参数γ规定在（0，1）范围内；每次迭代后的界限参数进行更新，如公式（11）所示：

使用新的界限参数继续进行反演，直至总目标函数收敛。

一种大数据量重力与重力梯度联合的快速正反演系统，其特征在于，包括正演模块和反演模块，其中，正演模块包括：

灵敏度矩阵单元，用于获取网格剖分参数和测点坐标数据，并基于BTTB结构矩阵计算得到压缩灵敏度矩阵，并存储；

正演数据单元，对模型密度矩阵进行零填充处理，并实现与压缩的灵敏度矩阵相同维度，通过二维快速傅里叶变换得到有效的正演数据；

所述反演模块包括：

灵敏度矩阵组合单元，用于获取网格剖分参数和测点数据，计算压缩灵敏度矩阵，并将各分量的压缩灵敏度矩阵拼合成组合灵敏度矩阵；

反演结果计算单元，用于将组合灵敏度矩阵用于预条件共轭梯度法中的增量方程求解，通过迭代计算输出反演结果。

本发明的有益效果是：

1、考虑到海量数据和精细划分的需求，采用规则矩形棱柱体剖分和对应规整测网布设，利用重力和重力梯度灵敏度矩阵的BTTB结构进行压缩，实现各分量的快速正演；

2、反演过程采用高斯-牛顿法外循环迭代求解模型，内部则利用共轭梯度法迭代求解增量方程，过程中存在的大量矩阵运算，可利用结构化矩阵快速正演对其进行加速；

3、本发明不需要牺牲精度来实现重力与重力梯度的快速正演；

4、本发明的正反演方法能节省更多的存储资源和计算时间；

5、本发明在实现高精度反演的同时兼顾了反演方法与快速正演的结合，实现了高精度大数据量重力与全张量重力梯度的快速联合反演。

附图说明

图1为本发明快速正演流程图。

图2为本发明实施例灵敏度矩阵压缩示意图。

图3为本发明实施例频率域卷积示意图。

图4为本发明反演流程图。

图5为本发明对比例2各分量误差图。

图6为本发明对比例3精细剖分组合模型。

图7为本发明对比例3各分量正演数据等值线图。

图8为本发明对比例3在不同深度层面的水平截面图。

具体实施方式

为了使本领域的普通技术人员能更好的理解本发明的技术方案，下面结合附图和实施例对本发明的技术方案做进一步的描述。

参照附图1~8所示的一种大数据量重力与重力梯度联合的快速正反演方法及系统，包括正演过程和反演过程，其中，正演部分包括基于BTTB结构矩阵的灵敏度矩阵压缩和二维快速傅里叶变换的快速正演。

实施例

正演过程包括以下步骤：

1）基于BTTB结构矩阵的灵敏度矩阵压缩

如图2（A）所示，为单层3×3的示例模型，观测点坐标与单元中心坐标一一对应，进行快速正演的流程如图1所示。

S2、根据模型单元中心坐标和观测点坐标计算BTTB结构矩阵所对应的相对坐标矩阵(，/>，/>，在该实施例中，相对坐标/>是不变的，因此只关注图2（B、C）中的相对坐标(，图2（C）为压缩后的灵敏度矩阵，其与完整的灵敏度矩阵（图2（B））的关系即为块Toeplitz矩阵压缩前后的关系。其中图2（B）所示的完整灵敏度矩阵相同颜色的块矩阵是完全相同的，且都是Toeplitz矩阵。

需要强调的是，根据第二步中BTTB结构矩阵所对应的相对坐标矩阵是重力和全张量重力梯度分量计算压缩灵敏度矩阵都需要的，且是完全相同的，只是各分量正演公式不同。

S3、由于重力与全张量重力梯度分量根据正演公式计算灵敏度矩阵时所需的相对坐标矩阵相同，因此这里只需要计算一次相对坐标矩阵即可，最后代入相应分量的正演灵敏度矩阵计算公式既能得到相应的压缩灵敏度矩阵（图2（C））。根据图2（B、C）的关系可知理论的灵敏度矩阵压缩比为。

S4、将压缩后的灵敏度矩阵进行二维快速傅里叶变换得到频率域的灵敏度矩阵，并存储下来。

至此完成了图3（A）部分的计算，将进行二维快速傅里叶变换的快速正演。根据Toeplitz矩阵的性质，只需在频率域将压缩后的灵敏度矩阵与模型密度矩阵进行卷积运算，即可获得与原始正演结果完全一致的结果。

2）二维快速傅里叶变换的快速正演

S5、将模型密度向量排列为矩阵，并用零填充到与压缩灵敏度矩阵相同维度，得到的模型密度矩阵。

S6、对填充后的三维模型密度矩阵按层做二维快速傅里叶变换，得到频率域的模型密度矩阵（图3（B））。

S7、将频率域的压缩灵敏度矩阵和模型密度矩阵进行点乘操作，得到包含正演结果的频率域数据矩阵，计算公式为：

其中，为每个深度层对应的压缩灵敏度矩阵，/>为每个深度层对应的模型密度矩阵，/>为正演观测数据，/>代表正向FFT运算，/>代表逆FFT运算。

S8、公式（1）中正演结果的频率域数据矩阵经过傅里叶逆变换得到空间域数据矩阵，从空间域数据矩阵中提取有效正演数据，图3（C）中d所在的区域。

反演部分是在高斯-牛顿法基础上，采用预条件共轭梯度法求解增量方程，从而实现快速且高精度的重力与全张量重力梯度联合反演。预条件共轭梯度反演过程如下：

S9、获取网格剖分参数、重力和全张量重力梯度数据，计算相应的压缩灵敏度矩阵；

S10、将重力和全张量重力梯度的压缩灵敏度矩阵与所述正演数据拼合成联合反演的组合灵敏度矩阵；

S11、算法利用吉洪诺夫正则化方法来平衡数据目标函数与模型目标函数，以确保反演的有效性。为避免在反演结果中出现不合理的密度值，通常需要采用物性上下限约束来控制反演结果的取值范围。

本发明采用了对数界限法，通过引入界限函数对超出约束范围的结果进行惩罚，以确保反演结果保持在合理的范围内。加入对数项后，总的目标函数如下所示：

式中，为数据目标函数；/>是模型目标函数，/>为正则化参数，最后一项是界限函数，/>是密度约束上下界限，/>是界限参数，表示界限函数在总目标函数所占的比重。在联合反演时数据目标函数为公式（3）：

其中，不同分量以如下形式组合：

，/>，/>

模型目标函数如公式（4）：

其中，表示重力和重力梯度分量的拼合灵敏度矩阵，/>为数据加权项，/>为拼合的观测数据向量，/>为模型加权矩阵，/>为参考模型。

S12、为求解该非线性方程组，本发明采用高斯-牛顿法进行迭代求解。在第次迭代时，把解模型表示为/>，利用当前解模型/>以及解模型变化方向/>构成新的解模型，表示为

为了最小化总目标函数（2），将其改写为：

利用泰勒级数展开，在时，得到增量方程：

式中，，/>，/>，。

S13、利用预条件共轭梯度法求解公式（7）得到解模型的变化方向。在确定初始模型/>后，就可以计算/>的初始值，以平衡对数界限函数与目标函数中的其余两项。随后，反演算法的主要任务是解决公式（7）。由于公式（7）的求解计算成本较高，因此采用共轭梯度法。

共轭梯度算法是一种仅利用一阶导数信息的方法，它克服了最速下降法收敛速度较慢的问题，同时避免了牛顿法需要存储和计算Hesse矩阵以及矩阵求逆的不便之处。该算法具有存储要求低、快速收敛、高稳定性等优点，在求解大规模线性方程时是最有用的方法之一。应用共轭梯度法来解公式（7）带来两大优势：

首先，公式（7）左侧的系数矩阵中包括，/>，/>和/>，这些都是稀疏矩阵，共轭梯度法在处理稀疏矩阵时表现出显著的计算效率。

其次，在不要求精确解的情况下，共轭梯度法能够在较短的时间内提供一个在允许误差范围内的解，这与寻求获得一个近似解的目标相符。

为了便于表述，将公式（7）改写成一个简单形式:

其中，

,/>,/>。

可以证明对任意，矩阵/>都是对称且正定的。因此，采用共轭梯度法来求解公式（7）是满足共轭梯度法求解条件的。

S14、高斯-牛顿迭代法是在极小点的估计值附近做泰勒级数展开式，进而找到极小点的下一个估计值，通过反复迭代直到求得合适的解模型。

然而，高斯-牛顿法在某种程度上只是确定了解模型的搜索方向，其搜索步长还需另外确定。如果步长太大，会使得迭代不收敛，从而发散，步长太小，会使得收敛速度很慢。因此，本发明采用有限步长法来确定一个合适的步长。

其中，是最大搜索步长，由公式（14）计算：

式中，参数γ表征步长的减少量，以确保解模型密度值不靠近上限和下限。参数γ规定在（0，1）范围内，每次迭代后的界限参数进行更新，如公式（11）所示：

使用新的界限参数继续进行反演，直至总目标函数收敛。

S15、反演算法中的阈值越小，迭代次数越多，但并不是迭代次数越多精度越高。因此，在实际应用中，需要合理设置迭代终止条件。初始值的选择对最终的反演结果不产生影响，仅会影响迭代次数。在满足收敛条件后，迭代停止，输出反演结果。

本发明还提供一种大数据量重力与重力梯度联合的快速正反演系统，包括正演模块和反演模块，其中，正演模块包括：

基于BTTB结构矩阵的灵敏度矩阵压缩单元，用于计算压缩灵敏度矩阵以降低数据量消耗的内存，后将压缩灵敏度矩阵进行存储；

二维快速傅里叶变换的快速正演数据单元，用于将模型密度向量排列成三维矩阵，并零填充至与压缩的灵敏度矩阵相同维度，后经二维快速傅里叶变换在频率域计算得到有效的正演数据；

所述反演模块包括：

灵敏度矩阵组合单元，用于获取网格剖分参数和测点数据，计算压缩灵敏度矩阵，并将各分量的压缩灵敏度矩阵拼合成组合灵敏度矩阵；

反演结果计算单元，用于将组合灵敏度矩阵用于预条件共轭梯度法中的增量方程求解，通过迭代计算输出反演结果。

对比例1：系列模型传统正演与本发明快速正演算法对比

为了对比本发明快速正演算法与传统正演算法的加速效果与处理数据的能力，进行如下正演对比试验(计算机配置：AMD Ryzen 9 5900HX CPU @ 3.30GHz, 内存 32GB)，需要强调的是，传统正演算法只有通过循环计算每个测点的累加值才能得到正演数据，并未存储灵敏度矩阵，因为内存所限无法存储完整的灵敏度矩阵。而按照本发明的快速正演方法，灵敏度矩阵是被压缩存储起来的，并不需要重复计算，因此，能极大地减少正演计算耗时，同时还能处理精细剖分的大数据量重力与重力梯度分量正演，对比结果如表1所示。

表1 系列模型正演耗时对比

对比例2：球体模型正演精度对比

将该快速算法的正演结果与球体解析公式的正演结果进行了比较，以验证基于矩形棱柱体剖分的数值正演与解析公式正演的理论值之间的误差。接下来，将这一结果与常规矩形棱柱体剖分数值正演算法（未使用灵敏度矩阵压缩和频率域卷积）的结果进行了对比。

球体模型的矩形棱柱体剖分模型源网的大小为2 km×2 km×1 km，网格剖分数量为200×200×100，网格剖分间距均为10 m。球体模型的中心坐标为(1 km,1 km,0.5 km)，模型密度为。通过计算分割后的体积，转换为理论球体模型的半径为/>。通过计算七个分量的误差，各分量误差如图5所示，从图5各分量误差图可以看出，/>分量的误差量级为/>，而其余六个梯度分量的误差量级均在/>左右。这些误差主要是由于矩形棱柱体剖分无法完全拟合球体模型所导致的。

为了更直观地比较误差与各分量数据的相对关系，计算各分量的相对均方根误差RRMS（公式（12））和最大相对误差MRE（公式（13）），计算结果如表2所示。其中，分量的最大相对误差最大，为0.1349%，/>和/>分量的相对均方根误差最大，为0.0089%，可以看出误差相对解析正演数据来说是非常小的。

表 2 球体模型各分量正演相对均方根误差和最大相对误差

对比例3：精细剖分大数据量组合模型反演

为了验证本发明提出的一种大数据量重力与重力梯度联合的快速正反演算法的可行性以及应对精细剖分大数据量组合模型的能力，本实施例采用一个多块体的精细剖分组合模型（图6）反演，反演流程如图4所示：

获取网格剖分参数、重力和全张量重力梯度数据，计算相应的压缩灵敏度矩阵，该模型总共包含7个不同形状，不同大小和不同深度的块体。该模型的区域为（0~2560m,0~2560m,0~1280m），模型剖分数量为256×256×64，单个矩形棱柱体大小为10m×10m×20m，观测点数量为256×256个，剩余密度均为1。由此可计算出单个重力或重力梯度分量的灵敏度矩阵维度为（65536,4194304），按双精度浮点型存储的内存需求超过2000GB。而压缩后的灵敏度矩阵维度为（511,511,64），同样按双精度浮点型存储的内存需求不到134MB，完整全部压缩灵敏度矩阵计算的耗时为3.23秒。该模型的重力与重力梯度分量正演数据如图7所示。

将重力和全张量重力梯度的压缩灵敏度矩阵拼合成联合反演的组合灵敏度矩阵，七个分量组合也仅需不到1GB的内存。

反演流程中，收敛阈值均设定为0.01，上述双重循环的反演算法，最终经过47次高斯-牛顿法迭代，反演时间为8438.1492秒，迭代收敛时的数据拟合均方根误差为。最后通过与理论模型对比计算模型拟合差均方根误差为/>，根据几个目标体底部深度不同，在不同深度层面上做了水平截面图，如图8所示，可以看出，无论是深度还是水平上都能很好地刻画出异常体。从而证明本发明的重力和重力梯度联合快速算法是准确有效的，对于大数据量的反演速度也是令人满意的。

本发明的原理是：该方法包括正演部分和反演部分，正演部分包括基于BTTB结构矩阵的灵敏度矩阵压缩和二维快速傅里叶变换的快速正演；反演部分是在高斯-牛顿法基础上，采用预条件共轭梯度法求解增量方程，从而实现快速且高精度的重力与全张量重力梯度联合反演，从而实现在保证精度的前提下实现了加速。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。

Claims (10)

1.一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，包括正演过程和反演过程，其中，正演过程包括以下步骤：

获取网格剖分参数和测点坐标数据，基于BTTB结构矩阵计算得到灵敏度矩阵，压缩并存储；

对模型密度矩阵进行零填充处理，并实现与压缩的灵敏度矩阵相同维度，通过二维快速傅里叶变换得到有效的正演数据；

所述反演过程包括以下步骤：

获取网格剖分参数和测点数据，计算压缩灵敏度矩阵，并将各分量的压缩灵敏度矩阵拼合成联合反演的组合灵敏度矩阵；

将组合灵敏度矩阵用于预条件共轭梯度法中的增量方程求解，通过迭代计算输出反演结果。

2.根据权利要求1所述的一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，所述的网格剖分参数和测点坐标数据，基于BTTB结构矩阵计算得到灵敏度矩阵，压缩并存储的步骤，具体包括：

S1、获取网格剖分参数、模型单元密度向量、观测点坐标；

S2、计算BTTB结构矩阵所对应的相对坐标矩阵(Δx，Δy，△z)；

S3、根据相对坐标矩阵同时计算重力和全张量重力梯度分量的正演灵敏度矩阵；

S4、将压缩后的灵敏度矩阵进行二维快速傅里叶变换得到频率域的灵敏度矩阵，并存储。

3.根据权利要求2所述的一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，所述的对模型密度矩阵进行零填充处理，并实现与压缩的灵敏度矩阵相同维度，通过二维快速傅里叶变换得到有效的正演数据的步骤，具体包括：

S5、将模型单元密度向量按层排列为三维矩阵，并用零填充到与压缩灵敏度矩阵相同维度；

S6、对填充后的三维模型密度矩阵按层做二维快速傅里叶变换，得到频率域的模型密度矩阵；

S7、将频率域的压缩灵敏度矩阵和模型密度矩阵进行点乘操作，得到包含正演结果的频率域数据矩阵；

S8、正演结果的频率域数据矩阵经过傅里叶逆变换得到空间域数据矩阵，从空间域数据矩阵中提取有效正演数据。

4.根据权利要求3所述的一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，在步骤S7中，所述频率域数据矩阵的计算公式为：

其中，为每个深度层对应的压缩灵敏度矩阵，/>为每个深度层对应的模型密度矩阵，/>为正演观测数据，/>代表正向FFT运算，/>代表逆FFT运算。

5.根据权利要求1所述的一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，所述获取网格剖分参数和测点数据，计算压缩灵敏度矩阵，并将各分量的压缩灵敏度矩阵拼合成联合反演的组合灵敏度矩阵，及将组合灵敏度矩阵用于预条件共轭梯度法中的增量方程求解，通过迭代计算输出反演结果的步骤，具体包括：

S9、获取网格剖分参数、重力和全张量重力梯度数据，计算相应的压缩灵敏度矩阵；

S10、将重力和全张量重力梯度的压缩灵敏度矩阵拼合成联合反演的组合灵敏度矩阵；

S11、根据吉洪诺夫正则化方法，加入对数界限法约束物性上下限，得到反演迭代求解的目标函数；

S12、将总目标函数按高斯-牛顿法的原理转换为增量方程的形式，高斯-牛顿法包含一次正演矩阵乘法用快速正演替代；

S13、利用预条件共轭梯度法求解每一步高斯-牛顿法迭代的增量，预条件共轭梯度法每次迭代包含一次正演矩阵乘法和一次矩阵转置乘法；

S14、预条件共轭梯度法求得模型增量后，在高斯-牛顿法中用有限步长防止迭代不收敛；

S15、设定阈值和收敛条件，判断迭代收敛后输出反演结果。

6.根据权利要求5所述的一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，在步骤S11中，总的目标函数为：

式中，为数据目标函数；/>是模型目标函数，/>为正则化参数，最后一项是界限函数，是密度约束上下界限，/>是界限参数，表示界限函数在总目标函数所占的比重。

7.根据权利要求6所述的一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于：在联合反演时数据目标函数为公式（3）：

其中，不同分量以如下形式组合：

，/>，/>

模型目标函数如公式（4）：

其中，表示重力和重力梯度分量的拼合灵敏度矩阵，/>为数据加权项，/>为拼合的观测数据向量，/>为模型加权矩阵，/>为参考模型。

8.根据权利要求5所述的一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，在步骤S12中，增量方程的转换过程为：在第次迭代时，把解模型表示为/>，利用当前解模型/>以及解模型变化方向/>构成新的解模型/>，表示为：

为了最小化总目标函数（2），将其改写为：

利用泰勒级数展开，在时，得到增量方程：

式中，，/>，/>，。

9.根据权利要求5所述的一种大数据量重力与重力梯度联合的快速正反演方法，其特征在于，在步骤S14中，采用有限步长法确定合适的步长：

其中，是最大搜索步长，由公式（10）计算：

式中，参数γ表征步长的减少量，参数γ规定在（0，1）范围内；每次迭代后的界限参数进行更新，如公式（11）所示：

使用新的界限参数继续进行反演，直至总目标函数收敛。

10.一种大数据量重力与重力梯度联合的快速正反演系统，其特征在于，包括正演模块和反演模块，其中，正演模块包括：

灵敏度矩阵单元，用于获取网格剖分参数和测点坐标数据，并基于BTTB结构矩阵计算得到压缩灵敏度矩阵，并存储；

正演数据单元，对模型密度矩阵进行零填充处理，并实现与压缩的灵敏度矩阵相同维度，通过二维快速傅里叶变换得到有效的正演数据；

所述反演模块包括：

灵敏度矩阵组合单元，用于获取网格剖分参数和测点数据，计算压缩灵敏度矩阵，并将各分量的压缩灵敏度矩阵拼合成组合灵敏度矩阵；

反演结果计算单元，用于将组合灵敏度矩阵用于预条件共轭梯度法中的增量方程求解，通过迭代计算输出反演结果。