内容摘要：在地球物理学中，多通过观测场上的信号，来反演源的信息。但因为每次观测都包含误差，通常需要多观测一些数据作为冗余观测，来压制噪声干扰。最小二乘准则用来作为估计准则确定模型。Geoist工具包提供了Misfit类来解决此类问题。无论线性问题还是非线性问题均适用。
1、三点定震中问题
今天我们以地震震中估计这个经典问题为例。我们知道，地震会触发地震波，其在地球介质中传播的波类型不同，波速也有差异。从地震仪检测到这些波的到达时刻（图1），就可以估计出地震波是从什么位置发射出来的。
图1 地震波到时示意图
因此，通过P、S波的到时差可以计算震中位置，每个台站的P和S波到时差大小乘波速差，就可以画一个圆圈，边界位置就是可能的震中位置。理论上有三个台站接收到地震波信号，量出到达时间，通过三个圆的交汇点就可以计算出震中的位置。
图2 震中估计原理示意图
2、程序实现
在Geoist的inversion模块里面，我们可以从Misfit类继承，本文问题的派生类为Homogeneous，根据具体问题的不同，实现init初始化类和predicted、jacobian两个类后即可完成一个实际问题类的定义。
下面的代码中，即是Homogeneous的完整定义，其中p是带估计的值，这里为震中，用x，y坐标表示。也可以是一个向量。
class Homogeneous(Misfit):
def __init__(self, ttres, recs, vp, vs):
super().__init__(data=ttres, nparams=2, islinear=False)
self.recs = np.array(recs)
self.vp = vp
self.vs = vs
def predicted(self, p):
"Calculate the predicted travel time data given a parameter vector."
x, y = p
alpha = 1/self.vs - 1/self.vp
pred = alpha*np.sqrt((self.recs[:, 0] - x)**2 +
(self.recs[:, 1] - y)**2)
return pred
def jacobian(self, p):
"Calculate the Jacobian matrix for the inversion."
x, y = p
alpha = 1/self.vs - 1/self.vp
sqrt = np.sqrt((self.recs[:, 0] - x)**2 +
(self.recs[:, 1] - y)**2)
jac = np.empty((self.ndata, self.nparams))
jac[:, 0] = -alpha*(self.recs[:, 0] - x)/sqrt
jac[:, 1] = -alpha*(self.recs[:, 1] - y)/sqrt
return jac
有了解决实际问题的具体派生类后，之后就可以套用这套反演框架来实现求解，通过创建类的实例对象：
solver = Homogeneous(traveltime, recs, vp, vs)
得到solver，然后用config来选择求解算法，fit开始计算，estimate_为最优解。因为该问题是非线性的，可以采用马奎特方法这里的标签为lwvmarq，任意给一个求解算法的初值即可求解。方法如下：
initial = (1, 1)
estimate = solver.config('levmarq', initial=initial).fit().estimate_
我们做了一个完整的例子，详见图3所示。
图3 示例程序的震中估计结果
一句话总结：今天的例子很简单，重点期望大家了解Misfit这个类的使用方法，通过这个入口大家可以了解Geoist背后的一些知识，在此基础上来实现更复杂问题的建模和求解。如果这些内容能加速你解决实际问题，那么Geoist的存在和发展就有意义了。
附完整代码如下：
import numpy
from geoist.pfm import giutils
from geoist.inversion.geometry import Square
import matplotlib.pyplot as plt
from geoist.vis import giplt
from geoist.inversion import ttime2d
from geoist.inversion.epic2d import Homogeneous
# Make a velocity model to calculate traveltimes
area = (0, 10, 0, 10)
vp, vs = 2, 1
model = [Square(area, props={'vp': vp, 'vs': vs})]
src = (5, 5)
srcs = [src, src, src, src]
recs = [(1, 2),(3,6),(4,7),(2,8)]
#giutils.connect_points(src, rec_points)
ptime = ttime2d.straight(model, 'vp', srcs, recs)
stime = ttime2d.straight(model, 'vs', srcs, recs)
# Calculate the residual time (S - P) with added noise
traveltime, error = giutils.contaminate(stime - ptime, 0.05, percent=True,
return_stddev=True)
solver = Homogeneous(traveltime, recs, vp, vs)
initial = (1, 1)
estimate = solver.config('levmarq', initial=initial).fit().estimate_
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.title('Epicenter + %d recording stations' % (len(recs)))
plt.axis('scaled')
giplt.points(src, '*y', label="True")
giplt.points(recs, '^r', label="Stations")
giplt.points(initial, '*b', label="Initial")
giplt.points([estimate], '*g', label="Estimate")
giplt.set_area(area)
plt.legend(loc='lower right', shadow=True, numpoints=1, prop={'size': 12})
plt.xlabel("X")
plt.ylabel("Y")
ax = plt.subplot(1, 2, 2)
plt.title('Travel-time residuals + error bars')
s = numpy.arange(len(traveltime)) + 1
width = 0.3
plt.bar(s - width, traveltime, width, color='g', label="Observed",
yerr=error)
plt.bar(s, solver.predicted(), width, color='r', label="Predicted")
ax.set_xticks(s)
plt.legend(loc='upper right', shadow=True, prop={'size': 12})
plt.xlabel("Station number")
plt.ylabel("Travel-time residual")
plt.show()
网友评论