内容摘要：书接上文，了解了反演的目的意义后，我们开始重力反演的第一个例子，并介绍Geoist中inv3d接口的使用方法。以及通过实例了解使用Geoist反演功能需要知道的相关类和函数使用方法。

1、基本概念

（1）观测方程

其中，表示场源，表示观测到的异常场，为核函数。

通常，用矩阵表示，而和是一维向量。反演问题是为已知，求的过程。正演反之。

一般在实际计算过程中，如果观测点数量为，场源的剖分数量为，用矩阵表示的维数为行×列。

对于已知场源求异常的正演过程来说，上述方程计算是稳定的，没有问题，关键在于反演。

如果，方程为典型的欠定类型。这反演问题是没有唯一解的，因此，必须通过正则化等处理后才能求解。

当然，如果大家问，是否能通过增加观测点的数量来让方程数量增加，而达到有唯一解的目的呢？理论上来讲是可以的，但要求观测系统的设计满足求解重力位场泊松方程的边值问题条件。

或者简单的说一下，如果你想反演地球内部的密度分布，必须让观测系统遍布整个地球表面才行。但是这在现实中通常无法做到，因为我们的观测一般仅可以在有限的地表范围内实施，这种情况下，即使观测点数大于场源剖分的单元数量，矩阵的条件数也会非常大，求解还是很困难的，而这时候还需要通过正则化的方法改善核矩阵的性质。

因此，下面我们说一下正则化问题。

（2）正则化

正则化（Regularization），最早由Tikhonov提出，也叫洪吉诺夫正则化。正则化后的反演目标函数数学形式可表示为：

上式中右端第二项可以通常称为正则项，其核心思想是引入一个罚函数，在原来的norm项中加入对模型m的约束。

其矩阵实现，等价于以下形式：

其中，

对比上式中的和上一部分的，扩展后的核矩阵性质发生了变化。给定一个合适的后，反演问题则可以获得唯一解。

（3）求解

通常，反演是在最小二乘解意义下求解，表示为：

方程左端的矩阵部分，一般是对称正定的。这时候可以用共轭梯度、加入预条件子等多种方法求解。论线性方程组的解放当然有很多种，最直接的解法就是矩阵求逆，可以表示为：

重磁位场反演问题，是很难的。就像你把很多个照片给ps到一起，加起来容易，但是要再分开可就难啦。

不止是重磁位场反演，对于大多数的地球物理反演问题本来在数学上是没有唯一解的，但是，工程实践中需要得到一个合理的结果。

因此，可以研究的方向有很多，我们总结几点：

（1）约束条件，也可以说是正则项的选择方法研究。要想改善反演结果，让其更接近于真实的地质情况，控制反演过程的有效方法就是正则项部分，对于很多已知的先验比如模型光滑、深度加权、参数有效范围都有很多研究成果，派生出各种针对特点问题的研究方法。特别是联合反演也可以看成此类，常见的重震联合反演，从几何相似性提出的交叉梯度联合反演等等。

（2）大规模计算，上面我们看到的核函数G，对于重力反演问题是非稀疏的（地震波层析成像的核函数G是稀疏的），因为场源对于每个观测都是有贡献的。稠密矩阵很麻烦，模型规模一大内存就扛不住了，100×100×100这么大规模的场源反演问题，对于PC或普通工作站基本都是常规方法无解的。小波压缩、分块计算等一些手段常常用于解决此类问题。

（3）反演求解，这个方向通常侧重于将线性方程解法应用于改善反演问题求解效率方面，因为对直接矩阵求逆来解方程，一般而言不太现实，常用迭代法更加节省内存开销。比如：代数重构法、预条件子法等解方程组的技巧。

2、设计模型

（1）设计场源模型

要进行实现重力反演、测试算法，首先是要设计一个场源模型，并实现其正演。PrismMesh函数可以生成一个指定场源参数的网格。

from geoist.inversion.mesh import PrismMesh
area = (-47500, 47500, -47500, 47500, 4000, 32000)
shape = (20, 20, 5)
mesh = PrismMesh(area, shape)

生成的mesh对象，可以通过addprop方法添加指定位置的密度或磁化率参数，然后，dump方法支持模型导出。
下图是采用MeshTools3D工具，对导出模型的可视化结果。

图1 场源模型

导出方法代码如下：

meshfile = r"d:\msh.txt"
densfile = r"d:\den.txt"
density=np.zeros(shape)
density[9:12,9:12,2:3]=1.00
mesh.addprop('density', density.ravel())
mesh.dump(meshfile, densfile, 'density')

（2）核函数生成

反演之前，最重要的一步要准备好核函数，也就是上面公式中G这个矩阵，也可以叫敏感度矩阵。每一个场源单元对每一个观测点都可以计算得到一个单位密度的理论异常值，全部组合起来就是一个M单元×N个测点数的大矩阵。可以先用一个空的list对象kernel来准备装该矩阵，代码如下：

kernel=[] 
narea = (-48750, 48750,-48750, 48750)
nshape = (40, 40)
xp, yp, zp = gridder.regular(narea, nshape, z=-1)
for i, layer in enumerate(mesh.layers()):
    for j, p in enumerate(layer):
        x1 = mesh.get_layer(i)[j].x1
        x2 = mesh.get_layer(i)[j].x2
        y1 = mesh.get_layer(i)[j].y1
        y2 = mesh.get_layer(i)[j].y2
        z1 = mesh.get_layer(i)[j].z1
        z2 = mesh.get_layer(i)[j].z2
        den = mesh.get_layer(i)[j].props
        model=[geometry.Prism(x1, x2, y1, y2, z1, z2, {'density': 1000})]
        field = prism.gz(xp, yp, zp, model)
        kernel.append(field)       
kk=np.array(kernel) 

最后，通过两个循环可以得到kk这个numpy的array类型数据。

（3）理论重力异常

有了模型m，有了核函数矩阵G，就可以计算理论重力异常了。下面代码计算得到的field变量，就是正演的理论观测重力异常值。

kk=np.transpose(kernel) 
field=np.mat(kk)*np.transpose(np.mat(density.ravel()))

3、程序实现

重力位场正反演的接口分别封装到pfm模块和inversion模块下，今天谈到的内容，完成程序如下：

import matplotlib.pyplot as plt
import numpy as np
from geoist import gridder
from geoist.inversion import geometry
from geoist.pfm import prism
from geoist.inversion.mesh import PrismMesh
from geoist.vis import giplt

meshfile = r"d:\msh.txt"
densfile = r"d:\den.txt"
#生成场源网格 10km*20km*5km, 500个单元，z方向5层
area = (-47500, 47500, -47500, 47500, 4000, 32000)
shape = (20, 20, 5)
mesh = PrismMesh(area, shape)
density=np.zeros(shape)
density[9:12,9:12,2:3]=1.00
mesh.addprop('density', density.ravel())
mesh.dump(meshfile, densfile, 'density') #输出网格到磁盘，MeshTools3D可视化
#生成核矩阵
kernel=[] 
narea = (-48750, 48750,-48750, 48750)
nshape = (40, 40)
xp, yp, zp = gridder.regular(narea, nshape, z=-1)
for i, layer in enumerate(mesh.layers()):
    for j, p in enumerate(layer):
        x1 = mesh.get_layer(i)[j].x1
        x2 = mesh.get_layer(i)[j].x2
        y1 = mesh.get_layer(i)[j].y1
        y2 = mesh.get_layer(i)[j].y2
        z1 = mesh.get_layer(i)[j].z1
        z2 = mesh.get_layer(i)[j].z2
        den = mesh.get_layer(i)[j].props
        model=[geometry.Prism(x1, x2, y1, y2, z1, z2, {'density': 1000})]
        field = prism.gz(xp, yp, zp, model)
        kernel.append(field)       
kk=np.array(kernel)        
kk=np.transpose(kernel)  #kernel matrix for inversion, 500 cells * 400 points
field=np.mat(kk)*np.transpose(np.mat(density.ravel()))
#画图
plt.title('gz')
levels = giplt.contourf(yp * 0.001, xp * 0.001, field, nshape, 15)
cb = plt.colorbar()
giplt.contour(yp * 0.001, xp * 0.001, field, shape,
                levels, clabel=False, linewidth=0.1)
plt.show()

图2 由模型正演得到的重力异常

一句话总结：反演之前当然要知道如何获得理论重力正演方法，上面仅举了一个重力核函数计算的例子，怎么对一个磁异常进行正演，如何生成重力梯度的核函数等等，别急，后面会逐一讲解。