17_Geoist三维重力反演_4

内容摘要：书接上文，说了半天还没有开始写一个完整的反演例子。这次我们以一个最简单的例子，来讨论一下如何通过反演得到一个三维密度结构。

1、问题定义

我们从一个最简单的最小构造模型入手。下式中右端第二项是正则部分(mnorm)，第一项为观测数据部分(dnorm)。

上式中为正则化因子，可以通过L-Curve算法找到一个最佳值，平衡mnorm和dnorm之间的关系。G为核函数，在反演之前可以得到，d是重力异常，m是要反演获得的密度结构。

问题简单来说：已知G，d，求m，优化，最小化。

2、程序实现

这样一个最简单的反演如何实现呢？核心代码如下所示：

print('starting inversion')
from geoist.inversion.regularization import Damping
from geoist.inversion.hyper_param import LCurve
from geoist.pfm import inv3d
datamisfit = inv3d.Density3D(np.array(field.T).ravel(), [xp, yp, zp], mesh)
regul = Damping(datamisfit.nparams)
regul_params = [10**i for i in range(-10, 5, 1)]
density3d = LCurve(datamisfit, regul, regul_params)
_ = density3d.fit()
print(density3d.regul_param_)
density3d.plot_lcurve()

上面代码重点看Density3D这个类，创建后实例对象为datamisfit，这部分可以理解为dnorm，Damping是一个生成最小构造正则项的类，实例化后的对象为regul。然后，指定一个要搜索的正则参数列表regul_params，最后，调用LCurve类，生成一个density3d对象。

运行density3d对象的fit接口，开始反演。LCurve支持多线程，当模型较大的时候可以多线程运行，njobs参数对应开的线程数量。

density3d = LCurve(datamisfit, regul, regul_params, njobs=4)

本例子最佳的值为0.01，曲线如图1所示。

图1 LCurve曲线

以反演得到的结果，导出后，用Meshtools3D画图如下：

图2 反演模型
真实模型如图3所示。
图3 真实模型

如果要显示一下残差和反演异常，代码如下：

predicted = density3d[0].predicted()
residuals = density3d[0].residuals()

如图4所示，这个理想的无噪声异常，残差很小了（图4右）。

图4 模型正演结果和残差

一句话总结：今天简单介绍了一下最简单情况下的重力正则化反演，没有深度加权，光滑等更复杂的正则项。别急，后面我们一点点深入讨论。

附完整代码：

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"
area = (-20000, 20000, -40000, 40000, 4000, 32000) #ny nx nz
shape = (10, 20, 5) #nz nx ny
mesh = PrismMesh(area, shape)
density=np.zeros(shape)
density[3:8,9:12,1:4]=1.0 # z x y
mesh.addprop('density', density.ravel())
mesh.dump(meshfile, densfile, 'density') #输出网格到磁盘，MeshTools3D可视化
#生成核矩阵
kernel=[] 
narea = (-28750, 28750,-48750, 48750)
nshape = (30, 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.figure()
plt.title('gz')
plt.axis('scaled')
levels = giplt.contourf(yp * 0.001, xp * 0.001, field, nshape, 15)
cb = plt.colorbar()
giplt.contour(yp * 0.001, xp * 0.001, field, nshape,
                levels, clabel=False, linewidth=0.1)
plt.show()
###反演
print('starting inversion')
from geoist.inversion.regularization import Damping
from geoist.inversion.hyper_param import LCurve
from geoist.pfm import inv3d
datamisfit = inv3d.Density3D(np.array(field.T).ravel(), [xp, yp, zp], mesh)
regul = Damping(datamisfit.nparams)
regul_params = [10**i for i in range(-10, 5, 1)]
density3d = LCurve(datamisfit, regul, regul_params)
_ = density3d.fit()
print(density3d.regul_param_)
density3d.plot_lcurve()

predicted = density3d[0].predicted()
residuals = density3d[0].residuals()
plt.figure(figsize=(16, 8))
plt.subplot(1, 2, 1)
plt.axis('scaled')
plt.title('inversed gravity anomlay (mGal)')
levels = giplt.contourf(yp * 0.001, xp * 0.001, predicted, nshape, 15)
cb = plt.colorbar(orientation='horizontal')
giplt.contour(yp * 0.001, xp * 0.001, predicted, nshape,
                levels, clabel=False, linewidth=0.1)

plt.subplot(1, 2, 2)
plt.axis('scaled')
plt.title('residual (mGal)')
levels = giplt.contourf(yp * 0.001, xp * 0.001, residuals, nshape, 15)
cb = plt.colorbar(orientation='horizontal')
giplt.contour(yp * 0.001, xp * 0.001, residuals, nshape,
                levels, clabel=False, linewidth=0.1)

print('res mean={:.4f}; std={:.4f}'.format(residuals.mean(), residuals.std()))

densinv = r"d:\deninv.txt"
values = np.fromiter(density3d.estimate_, dtype=np.float)
reordered = np.ravel(np.reshape(values, mesh.shape), order='F')
np.savetxt(densinv, 1000.*reordered, fmt='%.8f')    

normfile = r"d:\norms.txt"
with open(normfile, 'w') as f:
    for i in range(len(density3d.dnorm)):
        f.write('{} {} {}\n'.format(density3d.mnorm[i], density3d.dnorm[i]
                                    ,regul_params[i]))