美文网首页
番外篇10_Geoist之几种正则方法对比

番外篇10_Geoist之几种正则方法对比

作者: 地学小哥 | 来源:发表于2020-05-08 10:38 被阅读0次

内容摘要:重磁位场反演因为要面对求解病态方程问题,因此需要通过正则化方法来求解。如果你看的文献多了,会发现正则化方法类别很多,如何对比之间的效果,在合适的时候挑选合适的正则化工具呢?今天我们就Geoist实现的几种正则化方法,来看看效果有什么不一样。本节也会顺带介绍一下一个Geoist支持的L1范数意义下的TV正则约束。

1、正则化模型

我们前面已经讲过,反演的目标函数通常可以分为两部分,如下表示:

\phi = {||W_d(Gm-d)||}^2+\lambda {||W_m(m-m_{ref})||^2}

公式右端第二部分称为正则项,这个部分的定义方法可以变化很多种形式。今天讨论也主要针对这个右端第二项来展开。

首先,我们对正则化约束的数学公式写法,先简单定义一下,因为这样讲起来更清楚。

(1)最小模型约束

\phi = {||W_d(Gm-d)||}^2+\lambda {||m||^2}

(2)最小模型和深度加权约束

\phi = {||W_d(Gm-d)||}^2+\lambda {||W_m m||^2}

W_m = \alpha_sW_rI,其中,W_r 为深度加权。

(3)光滑和深度加权约束

\phi = {||W_d(Gm-d)||}^2+\lambda {||W_m m||^2}

W_m = \alpha_xB_xW_r+\alpha_yB_yW_r+\alpha_zB_zW_r

(4)最小模型、光滑和深度加权约束

\phi = {||W_d(Gm-d)||}^2+\lambda {||W_m m||^2}

W_m = \alpha_sW_rI+\alpha_xB_xW_r+\alpha_yB_yW_r+\alpha_zB_zW_r

(5)TV和深度加权约束

\phi = {||W_d(Gm-d)||}^2+\lambda {||W_m m||}

W_m = \alpha_xB_xW_r+\alpha_yB_yW_r+\alpha_zB_zW_r

注意,这里的第二项是L1范数的正则项。TV约束在求解过程中要用到非线性方法,优点是有利于保存图像的尖锐变化部分或边缘效果。

2、效果对比

下面我们分别看一下以上五种不同正则项的反演效果,下面五个结果都是通过LCurve方法得到的,拟合残差都很小,这里就不放残差图了。

图1 最小模型约束反演效果 图2 最小模型和深度加权约束 图3 光滑和深度加权约束 图4 最小模型、光滑和深度加权约束 图5 TV和深度加权约束

上面的效果差异,大家可以对照公式来理解,让子弹飞一会,详细的评论,我们回头再写吧。

3、实现代码

# -*- 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,TotalVariation
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 = field0 #giutils.contaminate(np.array(field0).ravel(), 0.05, percent = True)

#零均值
print(field.mean())
#field = field + 300. # field.mean()
#画图

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 = 300.0
ay = 300.0
az = 300.0
z0 = 10000 
beta = 1.0 #1.0

regul0 = Damping(nz*ny*nx)

wdepth = np.diag(1./(np.array(depthz)+z0)**beta)
sm0 = am*np.eye(nz*ny*nx)*wdepth
regul1 = Smoothness(sm0)

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)))
regul2 = Smoothness(sm) #np.eye(nz*ny*nx)

sm1 = np.vstack((az*np.dot(sz.T,wdepth),
                ay*np.dot(sy.T,wdepth),
                ax*np.dot(sx.T,wdepth)))
regul3 = Smoothness(sm1)

# L1 norm non-linear
regul4 = TotalVariation(0.0001, sm1)

#regul = Damping(nz*ny*nx)
datamisfit = inv3d.Density3D(np.array(field.T).ravel(), [xp, yp, zp], mesh
                             , movemean = False)

regul_params = [10**i for i in range(-5, 2, 1)]
# 目标函数
density3d = LCurve(datamisfit, regul4, regul_params, loglog=False)

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('levmarq', initial=np.ones(datamisfit.nparams), maxit=100, maxsteps=100, tol=10**-4).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:\deninv_reg4.txt"
#values = np.fromiter(density3d.objectives[8].p_, dtype=np.float)
values = np.fromiter(density3d.p_, dtype=np.float)
reordered = np.ravel(np.reshape(values, mesh.shape), order='F')
np.savetxt(densinv, reordered, fmt='%.8f')    

相关文章

  • 番外篇10_Geoist之几种正则方法对比

    内容摘要:重磁位场反演因为要面对求解病态方程问题,因此需要通过正则化方法来求解。如果你看的文献多了,会发现正则化方...

  • 第七天正则表达式

    正则表达式的理念: 正则的使用方法: 正则里string两个方法: 使用正则找字母的几种写法: 正则的转义: 正则...

  • Regex 正则表达式

    1. 优秀 正则方法对比,及优缺点http://blog.jobbole.com/109403/ 正则规则http...

  • 第六周第四天笔记

    1. 验证正则表达式中的全局g对几种方法的影响 test校验方法:正则的方法总结:1)不添加全局g时,lastIn...

  • 千分符

    把连续数字转换成逗号分隔数值的几种方法 方法一:使用正则表达式 语法: String(Number)...

  • 3-3 正则扩展

    和ES5的对比,以下方面 构造函数的变化 正则方法的扩展(ES6中调用了正则对象的方法) u修饰符 y修饰符 s修...

  • 吴恩达深度学习笔记(34)-你不知道的其他正则化方法

    除了L2正则化和随机失活(dropout)正则化,还有几种方法可以减少神经网络中的过拟合: 一.数据扩增 假设你正...

  • 几种人脸检测方法的对比分析

    简单的对比了下几种人脸检测(物体检测方法) [1]FCOSCenterFaceKPNetLFFDFSAF确定正负样...

  • PHP爬取网页

    主要流程就是获取整个网页,然后正则匹配(关键的)。 PHP抓取页面的主要方法,有几种方法是网上前辈的经验,现在还没...

  • Python数据科学:正则化方法!

    本文主要介绍,Python数据科学:正则化方法。正则化方法的出现,通过收缩方法(正则化方法)进行回归。 正则化方法...

网友评论

      本文标题:番外篇10_Geoist之几种正则方法对比

      本文链接:https://www.haomeiwen.com/subject/rwvgwhtx.html