番外篇8_Geoist之二维多边形拟合

内容摘要：通过任意多定点的二维多边形来拟合地质体，实现重力异常正演计算，在位场异常解释中十分常见。特别是交互式拟合异常曲线的时候，多采用多边形这种几何单元来近视地质体。多边形可以模拟界面的起伏，矿体形态等等。商业软件中Modelvision还可以提供2.5维的多剖面异常模拟与解释。今天我们讨论一下Geoist中的二维多边形正演和反演功能。

1、二维多边形异常计算

Talwani公式

图1 典型的地质剖面 图2 多边形地质体建模 图3 2.5维多边形截面地质体效果图 图4 二维多边形反演盆地基底

附完整代码：

from geoist.inversion import Smoothness1D
from geoist.pfm.basin2d import PolygonalBasinGravity
from geoist.pfm import talwani
from geoist.inversion.geometry import Polygon
from geoist.vis import giplt 
from geoist.pfm import giutils as utils
import numpy as np
import matplotlib.pyplot as plt

xs = np.linspace(0, 100000, 100)[::-1]
depths = (-1e-15*(xs - 50000)**4 + 8000 -
          3000*np.exp(-(xs - 70000)**2/(10000**2)))
depths -= depths.min()  # Reduce depths to zero
props = {'density': -300}
model = Polygon(np.transpose([xs, depths]), props)
x = np.linspace(0, 100000, 100)
z = -100*np.ones_like(x)
data = utils.contaminate(talwani.gz(x, z, [model]), 0.5, seed=0)

misfit = PolygonalBasinGravity(x, z, data, 50, props, top=0)
regul = Smoothness1D(misfit.nparams)
solver = misfit + 1e-4*regul
initial = 3000*np.ones(misfit.nparams)
solver.config('levmarq', initial=initial).fit()

plt.figure()
plt.subplot(2, 1, 1)
plt.plot(x, data, 'ok', label='observed')
plt.plot(x, solver[0].predicted(), '-r', linewidth=2, label='predicted')
plt.legend()
ax = plt.subplot(2, 1, 2)
giplt.polygon(model, fill='gray', alpha=0.5, label='True')
giplt.polygon(solver.estimate_, style='o-r', label='Estimated')
ax.invert_yaxis()
plt.legend()
plt.show()

参考文献：
Talwani, M, Worzel, JL, and Landisman, M, 1959,
Rapid Gravity Computations for Two-Dimensional Bodies with Application to the Mendocino Submarine Fracture Zone, Journal of Geophysical Research, 64:49-61.