内容摘要：书接上文，今天是最后一讲，小哥准备和大家再讨论一下上下界约束的反演方法。因为，在反演开始之前，我们对目标地质体的物性特征会有一定的先验认识。期望反演得到的密度结构在一定的有效范围内，但是前面介绍的反演方法，没有让我们加入约束的地方。今天我们介绍通过障碍函数的方法，来实现反演模型的上下界约束问题。
1、上下界约束
上下界约束的控制，可以通过障碍函数法实现。基本思想就是在反演过程中，如果结果超出上下界，就给目标函数一个很大的值，总是不让这个结果超出给定的范围。这个障碍函数法属于非线性最优化方法范畴，具体实现直接应用了Scipy里面的minimize函数，选择了trust-constr这个方法来引用障碍函数，详细算法可以参考文献（Conn，2000）。
在具体实现上，在config方法中，使用tcbound参数配置即可，如下所示：
density3d.config('tcbound', bounds = bounds, nparams = datamisfit.nparams, x0 = x0).fit()
对比上节课的代码发现哪里不一样了吗？原来的反演计算直接density3d.fit()即可，但是这个中间多了config这部分，默认是linear，还有levmarq和newton法等，根据不同需求可以灵活选择，但是参数配置可能不一样噢。tcbound这个需设置三个参数，第一个bounds是每个单元的上下界，第二个参数nparams是模型单元个数，第三个x0是初始值。
详细的写法，请参考后面完整代码。
2、效果测试
接下来我们看看上下界约束的反演效果，小哥实验了一下，速度真心不快，18个值，1000个单元的反演问题，算一遍半个小时吧。都快等得快睡着了，结果出来了，如图1所示。看到这个结果眼前一亮，真实模型图2这样子，中心有一个剩余密度为
的异常体，对全部模型加了0-1范围的约束，图1看上去还不错，我们也进行了多个维度的切片，知识下底界有点明显差别，其它地方都挺好。
图1 上下界约束反演效果
图2 真实模型
从残差效果看（图3），效果也不错，这个结果是两天来，最接近真实模型的结果了，看来我们的上下界约束还是起作用啦！
图3 反演结果的拟合残差
本节测试的例子很简单，没有什么异常体走向上的变化，还没有噪声干扰，也许这太理想化了，后续我们将进一步测试程序，期望广大科研工作者，不用再去重复发明轮子，在这些代码的基础上，任意扩展、节省时间，做更多有意义的事情。
一句话总结：说到这里，我们把密度结构反演问题基础讲完了。由于篇幅问题（实际是程序还没写）其实还有很多问题没讨论，比如：小波压缩实现大规模反演问题；多参数的贝叶斯优化问题；参考模型的引入问题等等。万事开头难，留点不完美吧，如果你感兴趣，关注小哥，后面番外篇再见！
参考文献：Conn, A. R., Gould, N. I., & Toint, P. L. Trust region methods. 2000. Siam. pp. 19.
完整代码如下：
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.pfm import giutils
from geoist.inversion.mesh import PrismMesh
from geoist.vis import giplt
from geoist.inversion.regularization import Smoothness,Damping
from geoist.inversion.pfmodel import SmoothOperator
from geoist.inversion.hyper_param import LCurve
from geoist.pfm import inv3d
meshfile = r"d:\msh.txt"
densfile = r"d:\den.txt"
#生成场源网格 NS-40km, EW-80km, 500个单元，z方向10层
area = (-20000, 20000, -40000, 40000, 2000, 32000) #NS EW Down
shape = (10, 20, 5) #nz ny nx
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) #NS, EW
nshape = (30, 40) #NS, EW
depthz = []
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)
depthz.append((z1+z2)/2.0)
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(figsize=(16, 8))
plt.subplot(1, 2, 1)
plt.title('gravity anomlay')
plt.axis('scaled')
levels = giplt.contourf(yp * 0.001, xp * 0.001, field0, nshape, 15)
cb = plt.colorbar(orientation='horizontal')
giplt.contour(yp * 0.001, xp * 0.001, field0, nshape,
levels, clabel=False, linewidth=0.1)
plt.subplot(1, 2, 2)
plt.title('gravity anomlay with 5% noise')
plt.axis('scaled')
levels = giplt.contourf(yp * 0.001, xp * 0.001, field, nshape, 15)
cb = plt.colorbar(orientation='horizontal')
giplt.contour(yp * 0.001, xp * 0.001, field, nshape,
levels, clabel=False, linewidth=0.1)
plt.show()
print('starting inversion')
smop = SmoothOperator()
nz = shape[0]
ny = shape[1]
nx = shape[2]
p = np.eye(nz*ny*nx).reshape(-1,nz,ny,nx)
sx = smop.derivation(p, component = 'dx').reshape(nz*ny*nx,-1)
sy = smop.derivation(p, component = 'dy').reshape(nz*ny*nx,-1)
sz = smop.derivation(p, component = 'dz').reshape(nz*ny*nx,-1)
am = 1000.0
ax = 100.0
ay = 100.0
az = 100.0
z0 = 10000
beta = 1.0 #1.0
wdepth = np.diag(1./(np.array(depthz)+z0)**beta)
sm = np.vstack((am*np.eye(nz*ny*nx)*wdepth,
az*np.dot(sz.T,wdepth),
ay*np.dot(sy.T,wdepth),
ax*np.dot(sx.T,wdepth)))
regul = Smoothness(sm)
datamisfit = inv3d.Density3D(np.array(field.T).ravel(), [xp, yp, zp], mesh
, movemean = False)
regul_params = [10**i for i in range(-8, 10, 1)]
density3d = LCurve(datamisfit, regul, regul_params, loglog=True)
initial = np.zeros(nz*ny*nx) #np.ones(datamisfit.nparams)
minval = initial
maxval = initial + 1.0
bounds = list(zip(minval, maxval))
x0 = initial + 1.0
_ = density3d.config('tcbound', bounds = bounds, nparams = datamisfit.nparams, x0 = x0).fit()
print('Hyperparameter Lambda value is {}'.format(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:\deninvmean.txt"
values = np.fromiter(density3d.p_, dtype=np.float)
reordered = np.ravel(np.reshape(values, mesh.shape), order='F')
np.savetxt(densinv, reordered, fmt='%.8f')
网友评论