内容摘要：书接上文，在实际的反演问题中，观测的重力异常都是包含噪声污染的。这种含噪声的数据会不会影响反演结果呢？还有，重力异常都是相对量，起算基准不同可能存在系统偏差，这种情况下应该怎么开展反演？别急，这节课我们一起来说。
1、噪声污染数据
实际的观测重力异常数据都包含噪声，会不会影响反演结果呢？我们首先生成一个含噪声的观测数据，还记得前面讲过的如何在理论数据中添加噪声吧？采用contaminate函数即可，代码如下：
field = giutils.contaminate(np.array(field0).ravel(), 0.05, percent = True)
图1
在我们上一节课的模型中，有五个参数必须事先给定，分别是
am = 1000.0， ax = 100.0， ay = 100.0， az = 100.0；深度加权参数z0 = 10000，beta = 1.0 。
在该模型的超参数搜索范围设定如下：
regul_params = [10**i for i in range(-7, 10, 1)]
程序运行后L-Curve优化结果如图2所示。
图2 L-Curve曲线（log-log）
曲线的拐点还是很漂亮的，这里注意data norm的部分很大，因为，包含噪声，反演模型并没有完全拟合观测的包含噪声的异常数据结果，模型认为是噪声部分就不再拟合了。
当然，如果感觉L-Curve找到的值不合适，也可以用objectives访问指定
参数的反演结果。
values = np.fromiter(density3d.objectives[8].p_, dtype=np.float)
那么，这个模型参数下的异常拟合结果如何呢？看图3左是不是还是与真实的无噪声结果不一样？是啊，这就是噪声的影响。
图3 模型重力异常和剩余异常
最后，看看反演得到的密度结构吧。注意图4中的浅部异常特征，由于观测数据中包含噪声，还是反映到反演密度结构模型中了，但是仅在浅部。
图4 反演密度结构特征
2、重力异常中的系统误差
前面在重力异常数据处理中，就讲过去趋势的方法。如果你反演使用的是剩余布格重力异常，那这个异常与真实模型可能存在一个系统误差。那么，这时候可以采用零均值反演，方法是设置movemean参数。
Density3D(np.array(field.T).ravel(), [xp, yp, zp], mesh, movemean = True)
如果开启这个参数，在异常输入模型和每次统计dnorm的时候，都可以去掉一个均值。我们在原来的异常里面，加上300mGal的偏差，来试试效果。
图5 自动去均值反演效果
反演过程中的去均值信息：
starting inversion
remove the mean value of anomaly...
(1000, 1200)
The mean 1.5485088421307581e-07 has been removed.
The mean -6.324391081798571e-08 has been removed.
The mean -1.1204735714905507e-06 has been removed.
The mean -2.3719417504490064e-06 has been removed.
The mean -2.94958342768048e-06 has been removed.
The mean -4.714087561386026e-06 has been removed.
The mean 2.127315945655184e-05 has been removed.
The mean 0.0001171615023231798 has been removed.
The mean 0.0002833518745389145 has been removed.
The mean 0.0002645561867915358 has been removed.
The mean -0.0022071256106698587 has been removed.
The mean -0.011567825931014084 has been removed.
The mean -0.04628737534019467 has been removed.
The mean -0.18042321669168682 has been removed.
The mean -0.7621151534379003 has been removed.
The mean -2.105420470427456 has been removed.
The mean -2.2641299498413847 has been removed.
The mean -0.5292459700177187 has been removed.
Hyperparameter Lambda value is 10000
res mean=-0.0463; std=3.5358
同样参数，如果不开启去均值设置后，反演结果如图6所示。由于开启了深度加权，剩余密度为正的异常体被压到底部，上面还出现了一些伴生的负异常。
图6 同样参数不开启去均值反演的效果
而对比两者的L-Curve曲线，差别不大。
图5结果的L-Curve曲线
图6结果的L-Curve曲线
一句话总结：今天讨论完后，是不是觉得距离解决实际问题又进了一步，如果每天你觉得懂得多了一点，这就是进步啦。我们这个例子只有1000个单元，1200个观测点，如果10000个单元的时候，每次反演可要一会才能完成。重力反演如果你做多了，最麻烦的其实就是调整参数了，这里面是预设6个参数，加上一个自动搜索，说实在的如果给不好一时半会也不会出来一个相对靠谱的结果。所以，调参是最折磨人的。小哥所在的课题组，近年来正在发展一种贝叶斯优化方法，可以对该反演问题中的全部参数自动调节到最优。如果你关注本教程，我们会以番外篇的形式来讨论。
附完整代码如下：
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 26 17:45:42 2020
@author: chens
"""
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
field0=np.mat(kk)*np.transpose(np.mat(density.ravel()))
field = giutils.contaminate(np.array(field0).ravel(), 0.05, percent = True)
#画图
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
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) #np.eye(nz*ny*nx)
datamisfit = inv3d.Density3D(np.array(field.T).ravel(), [xp, yp, zp], mesh)
regul_params = [10**i for i in range(-7, 10, 1)]
density3d = LCurve(datamisfit, regul, regul_params, loglog=True)
_ = density3d.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:\deninvbxyz.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')
