重力异常的反演

重力异常的反演计算，就是根据观测的重力异常求取地下场源形态和密度参数，给出产生此重力异常的剩余密度体在地球内部的分布。

9.5.1 反演的原则

反演问题是重力勘探，特别是资料处理解释中的主要环节之一，须遵循以下原则。

1.从观测异常中选择提取适当的用于反演的场值

进行反演必须采用从观测异常中分离出的单纯由反演目标引起的那部分异常值，可以不直接用观测场值而是先将它们作某些变换。例如，求出异常的导数和换算至不同高度等，再用变换后的场值作为反演用的原始数据。

2.选择适当的模型

为了逼近形状不规则的地质体，一般采用一组（多个）而不是单个模型，每个模型的几何形状要尽可能简单。在计算时关键问题是给定模型的初值。并非任意给一组初值，计算机都会计算出接近实际地质体的模型。所以，要根据实际的地质、钻井及其他地球物理资料给出接近实际地质体模型参数的初值。

3.应用合适的计算方法

反演计算过程能否快速、收敛并得到精确的结果，除了给出合适的初始模型外，还取决于所用的计算方法。应针对不同的问题选择合适的计算方法，而且在计算过程的不同阶段不断调整、修改所用的计算方法。

4.减少多解性（非唯一性）的影响

内特尔顿（1987年）指出（图9-7），与一定的重力异常宽度对应的场源可能的最大深度，就是引起同样宽度异常的点源（球体）的深度。在这个最大深度和地面之间，存在一个可能源的锥形区。同一个异常可以由埋藏很浅的薄透镜体（图9-7的“3”），或不太狭窄的较厚的物体，或球体引起。而且在球体上方较浅的深度上，却有无限多个可能的场源引起这个相同的异常。引起同一异常的不同的场源，有唯一一个共同的特性，那就是它们的剩余质量必须是相同的。

图9-7 引起相同异常的可能源的锥形区

显然，反演计算结果的可靠性是反演方法得以存在的根据，而重力反演中固有的多解性，严重地影响到计算结果的可靠性。采用适当的模型，特别是应用已知资料施加约束，可在一定程度上减少多解性的影响，并取得比较可靠的结果。

9.5.2 反演方法

在实际重力资料反演过程中，有直接计算密度异常体几何参数和物性参数的半定量、定量方法，或是针对密度界面的反演方法。

1.直接法

直接法是直接利用由反演目标引起的局部异常，通过某种积分运算和函数关系，求得与异常分布有关的密度异常体的某些参量，如三度体剩余质量和重心水平坐标。

考虑地下剩余质量体M，它产生的重力异常就是剩余质量体对地面各点单位质量的引力沿垂向方向的分力。如果将观测面当作无限大的水平面，下半空间用半径无限大的半球代替，则因上半空间无剩余质量，下半空间的立体角为2π。由场论中的高斯（Gauss）定理可知，位于任意封闭曲面S以内的质量M产生的引力场强，沿该曲面外法线方向的通量与质量M成正比，比例系数为-4πG。从而由观测面上重力异常的分布就可求出地质体剩余质量M：

而三度体重心水平坐标（x0，y0）可表示为

在具体计算时，不可能遍及全平面，因而要分为下述两部分：一部分为实测异常在有限测区内的求和，另一部分为有限测区外由测区内这部分地质体剩余质量M产生的异常求和。后者可采用如下办法求取：令R表示已经参与计算的异常区的折合半径，在计算区以外的重力异常可近似地用一剩余质量为M、中心埋深为D的均匀球体引起的异常来代替，当R远大于D时，它可近似为折合半径为R的圆周上的异常平均值g（R）与2πR2的乘积。

于是

仿照前面的方法，可以得到三度体重心水平坐标（x0，y0）：

2.特征点法

特征点法（或任意点法）是根据异常曲线上的一些点如特征点（如极大值点、零值点、拐点）的异常以及相应的坐标求取场源体的几何或物性参数，仅适用于剩余密度为常数的几何形体。

以水平圆柱体模型为例，在垂直走向的剖面图上，根据其产生的Δg异常求圆柱体中心埋藏深度D，线密度λ，半径R及顶部埋藏深度h的公式分别为

这里取Δg异常极值处为x=0，x1/n为Δg异常极值的1/n处的x坐标值。

而对铅垂台阶模型，在垂直走向的剖面图上，根据Δg异常求铅垂台阶厚度H-h的公式为

3.选择法

选择法的原理是根据实测重力异常在剖面或平面的分布和变化的基本特征，结合工区的地质、其他地球物理和物性等资料，给出引起异常的初始地质体模型，然后进行正演计算。将计算的理论异常与实测异常进行对比，两者偏差较大时，根据掌握的场源体资料对模型进行修改，重算其理论异常。再次进行对比，如此反复进行，以两种异常的偏差达到要求的误差范围时的理论模型来表示实际的地质体。目前在计算机上用得较多的是最优化选择法反演。

设实测异常为Δgk，k=1，2，…，M，M是用于计算的异常的点数。密度异常体模型引起的理论异常用Δg′k表示，k=1，2，…，M，它是模型体参数和计算点坐标的函数，即

Δg′k=f（xk，yk，zk，b1，b2，…，bn）（9-24）

式中，xk，yk，zk为计算点坐标；b1，b2，…，bn为模型的N个几何或物性参数。

判定实测异常与理论异常的符合程度的依据是最小二乘原则，即让两者偏差的平方和为最小：

函数Φ称为目标函数。在满足上式的条件下求取模型体的参数。求取模型体参数的过程一般需要下列迭代计算：

（1）给模型体初值bji（0），i=1，2，…，n；b的上标j表示第j次迭代，0即表示初值。

（2）根据模型体初值计算理论重力异常Δ。

（3）根据理论重力异常Δg′k与实测异常Δgk的差别计算模型体参数的修改量δbi，并按下式修改模型：

式中，l为迭代次数；b（l）i为本次迭代的初值；b（l+1）i为迭代结果。

（4）判断计算结果是否满足要求。如果不满足，则重复第（2）和第（3）步计算；如果满足，本次迭代得到的b（l+1）i就是最后得到的模型参数。

目前，应用的最优化方法有阻尼最小二乘法、广义逆矩阵法和共轭梯度法等。

4.密度界面迭代反演

1）线性回归法

这是计算密度分界面深度的近似方法。如果界面起伏平缓，可以认为重力变化与界面的起伏近似呈线性关系，即

h=a+b·Δg

式中，h为界面深度；Δg为界面起伏引起的重力异常；a，b是两个系数，它们与重力异常起算点处的界面深度和界面上下物质层的密度差有关。应用上式求深度，至少要知道界面上两个点的深度，以确定a，b两个系数值。如果已知点有n个，它们的深度为hi（i=1，2，…，n），则根据最小二乘原理，为确定系数a，b，应使各点由线性公式计算出的深度h′i与已知点深度hi的偏差的平方和为最小，即

可导出

系数a，b求出后，就可计算出各测点下方界面的深度。

本方法十分简便，二度与三度界面皆可适用，在界面起伏平缓时，误差是不大的。例如，当界面起伏最大倾角小于3°，起伏幅度不超过界面最大深度的1/10时，用这一方法求得结果的最大相对误差不超过7%。即使界面起伏最大倾角达到11°，起伏幅度达到界面最大深度的1/5，带来的最大相对误差也小于8%。

用这一方法求得的界面起伏比实际情况要平缓，越接近隆起的顶部或凹陷的中心，求得的界面深度的误差越大。

2）压缩质面法反演

用若干个宽度相同、顶深不同的直立板组合逼近两度体密度界面，设密度界面密度差为ρ，平均深度H已知，第j个直立板相对界面平均深度的高度为hj。将第i个直立板在地面上产生的重力异常等效成由在直立板处位于平均深度的质量面ρhj产生的：

式中，2a为板水平宽度，σj=ρhj。上式系数Aij可事先算出，在地面xi处观测到的总异常就是每一个质量面异常的叠加：

这是一个线性方程组，由地面各测点的重力异常Δg（xi）和系数Aij可求解出各个质量面的密度σj=ρhj，于是第j个直立板相对界面平均深度的高度为hj=。