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