内容摘要：书接上文，我们还是回到重力反演这个问题上来。通过正则项获得一个反演结果后，细心的朋友会发现还有一些问题，比如反演结果在空间上比较分散，而且高值都趋近于浅层界面。那么如何改善反演结果呢？今天我们就再介绍两个正则项的加法，即深度加权和光滑先验。

1、位场反演的正则设计

对于三维重力反演问题，最小化目标函数可以采用下面形式设计：

其中，通常可以分为几个部分，上一节我们仅用了一个对角阵，实现了最小构造约束。如果加上x，y，z方向的光滑约束，可以用如下形式表示：

上式中，分别是在最小模型，x，y，z三个方向光滑的系数；是深度加权，是一个与模型深度相关的对角阵；为光滑矩阵，通常可以选一阶差分或二阶差分形式。

ok，看完上面的这些正则设计，是不是有点晕？没关系我们直接来看在Geoist中怎么设计这部分。

2、正则项构建

上一节课，我们用Damping函数构建了一个最简单的正则项，其实就是一个对角单位阵，代码如下：

from geoist.inversion.regularization import Damping
regul = Damping(datamisfit.nparams) #nparams为模型单元个数

而上面说的为光滑项稍微有点复杂了，需要先生成差分矩阵。还需要用到pfmodel里面的SmoothOperator函数来搞定。下面是生成代码：

import numpy as np
from geoist.inversion.regularization import Smoothness
from geoist.inversion.pfmodel import SmoothOperator
smop = SmoothOperator()
nz = 2
ny = 2
nx = 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)
print(sx.T)
print(sy.T)
print(sz.T)

对应的分别为:

Bx =
[[ 1. -1.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  1. -1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  1. -1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1. -1.]]
By =
[[ 1.  0. -1.  0.  0.  0.  0.  0.]
 [ 0.  1.  0. -1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  1.  0. -1.  0.]
 [ 0.  0.  0.  0.  0.  1.  0. -1.]]
Bz =
[[ 1.  0.  0.  0. -1.  0.  0.  0.]
 [ 0.  1.  0.  0.  0. -1.  0.  0.]
 [ 0.  0.  1.  0.  0.  0. -1.  0.]
 [ 0.  0.  0.  1.  0.  0.  0. -1.]]

如果是二阶光滑呢？稍微变一下函数中的参数即可，代码如下：

sxx = smop.derivation(p, component = 'dxx').reshape(nz*ny*nx,-1)
syy = smop.derivation(p, component = 'dyy').reshape(nz*ny*nx,-1)
szz = smop.derivation(p, component = 'dzz').reshape(nz*ny*nx,-1)

好了既然这些矩阵都可以生成，那怎么构建正则部分呢？先把模型正则项部分合并到一起，然后再使用Smoothness函数即可。

ax = 0.1
ay = 0.2
az = 3.0
sm = np.vstack((ax*sx.T,ay*sy.T,az*sz.T))
smoothxyz = Smoothness(sm)
s1 = smoothxyz.hessian(np.array([0,0,0,0,0,0,0,0]))

s1的结果如下：

好了，我们看看反演的模型结果：

图1 加了光滑约束的反演结果

图2 只有最小模型约束的结果

对比图1和图2，是不是感觉加了光滑之后，反演的密度结构更平滑了一些？图2是仅加了最小模型约束的结果，两者的拟合残差没有明显区别。下图是图1对应的L-Curve最优结果。

图3 L-Curve的优化结果

3、深度加权

上面把光滑先验的正则加法讲完了，我们再聊聊的设计，为什么放到最后讲呢？因为，这个设计只有位场反演中有，算是一个特色部分吧。

我们都知道重力大小是随着距离的平方衰减的，也就是距离观测点越远，重力越小。而且这个衰减速度非常快。这个的结果就是当你生成核函数的时候，矩阵里面对应每个模型单元对观测点的贡献差别很大。

换句话说，对于浅部的模型单元，对观测平面贡献大，深部的就很小，这导致什么结果呢？当你的模型剖分在深度方向较多时，反演结果总是集中在浅部，这通常在地球物理上，叫“趋肤效应”。

前人们，因此提出了一个深度加权的概念，方式就是在模型正则上做文章。加一个特定的函数，与深度相关，以抵消这种趋肤影响。

形式当然有很多，但是多采用函数类似下面的形式：

上式中，为每个单元的深度，和是两个反演前给定的常数。

深度加权原理讲完了，我们看看代码怎么写：

am = 100.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)

上面代码为深度加权和光滑在一起的正则写法，包括了四部分，第一个是最小模型约束，后面三个分别为x，y，z三个方向的光滑。

深度加权的反演效果如图4，是不是异常体跑到下边去了。

图4 深度加权后的反演效果

真实的模型结果长什么样子呢？图5所示。

图5 真实模型结果

一句话总结：说到这里总结时间到了。大家一定发现了，图4的结果形态上虽然与真实模型差不多了，但是数值上面还差不少。后面我们会继续讨论怎么加反演结果的范围约束。另外，对于含噪声数据的反演还没测试，也在后面一起说吧。

图1结果的完整代码如下：

# -*- 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.inversion.mesh import PrismMesh
from geoist.vis import giplt
from geoist.inversion.regularization import Smoothness
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"
#生成场源网格 10km*20km*5km, 500个单元，z方向5层
area = (-20000, 20000, -40000, 40000, 4000, 32000) #ny nx nz
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)
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')

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 = 1.0
ax = 1.0
ay = 1.0
az = 1.0

sm = np.vstack((am*np.eye(nz*ny*nx),az*sz.T,ay*sy.T,ax*sx.T))
regul = Smoothness(sm)
datamisfit = inv3d.Density3D(np.array(field.T).ravel(), [xp, yp, zp], mesh)
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:\deninvbxyz.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')