有两种方法可以做到这一点。
- 使用非线性求解器
- 线性化问题并以最小二乘法解决
设定
因此，据我所知，您知道在4个不同点处的F，a，b和c，并且想要对模型参数X，Y和Z求逆。我们有3个未知数和4个观测数据点，因此这个问题太确定了。因此，我们将在最小二乘意义上进行求解。
在这种情况下，使用相反的术语更为常见，因此让我们翻转方程式。代替：
F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ
让我们写：
F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)
我们知道
F，
X，
Y，并
Z在4个不同点（例如
F_0, F_1, ... F_i）。
我们只是在更改变量的名称，而不是方程式本身。（这比我想的要容易得多。）
线性解
实际上有可能使该方程线性化。您可以轻松地解决
a^2，
b^2，
a b cos(c)，和
a bsin(c)。为了使操作更简单，让我们再次重新标记：
d = a^2e = b^2f = a b cos(c)g = a b sin(c)
现在，方程式变得简单得多：
F_i = d + e X_i + f Y_i + gZ_i。这很容易做到最小二乘线性反演
d，
e，
f，和
g。然后
a，我们可以从中获取
b，和
c：
a = sqrt(d)b = sqrt(e)c = arctan(g/f)
好吧，让我们用矩阵形式写出来。我们将翻译4个观测值（我们将编写的代码将采用任意数量的观测值，但现在让我们保持具体）：
F_i = d + e X_i + f Y_i + g Z_i
进入：
|F_0| |1, X_0, Y_0, Z_0| |d||F_1| = |1, X_1, Y_1, Z_1| * |e||F_2| |1, X_2, Y_2, Z_2| |f||F_3| |1, X_3, Y_3, Z_3| |g|
或：（
F = G * m我是地球物理学家，所以我们将其
G用于“格林函数”和
m“模型参数”。通常，我们也将
d“数据”而不是
F）。
在python中，这将转换为：
def invert(f, x, y, z): G = np.vstack([np.ones_like(x), x, y, z]).T m, _, _, _ = np.linalg.lstsq(G, f) d, e, f, g = m a = np.sqrt(d) b = np.sqrt(e) c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees return a, b, c
非线性解
您也可以使用
scipy.optimize@Joe建议使用来解决此问题。其中最易访问的函数
scipy.optimize是
scipy.optimize.curve_fit默认情况下使用Levenberg-
Marquardt方法。
Levenberg-
Marquardt是一种“爬坡”算法（在这种情况下，它会走下坡路，但是无论如何都使用该术语）。从某种意义上说，您可以初步估计模型参数（所有参数，默认情况下为
scipy.optimize），然后沿
observed- predicted参数空间的坡度向下倾斜至底部。
警告：
选择正确的非线性反演方法，初步猜测并调整方法的参数在很大程度上是一门“黑手艺”。您只能通过做来学习它，并且在很多情况下，事情将无法正常进行。如果您的参数空间相当平滑（应该是这样），那么Levenberg-
Marquardt是一个很好的通用方法。在其他情况下，还有很多其他方法（包括遗传算法，神经网络等，除了更常见的方法（如模拟退火）之外）。我在这里不打算深入探讨这一部分。
有一个常见的陷阱，一些优化工具包试图纠正这一错误，而
scipy.optimize没有解决。如果模型参数的大小不同（例如
a=1, b=1000,c=1e-8），则需要重新缩放比例，以使它们的大小相似。否则
scipy.optimize，“爬山”算法（如LM）将无法准确计算局部梯度的估算值，并且会给出非常不准确的结果。现在，我假设
a，
b以及
c具有相对类似的量值。另外，请注意，基本上所有非线性方法都需要您进行初始猜测，并且对此猜测很敏感。我将其保留在下面（将其作为
p0kwarg传递给它
curve_fit），因为默认值
a,b, c = 1, 1, 1是对的相当准确的猜测
a, b, c = 3, 2, 1。
避免了警告，
curve_fit希望传递一个函数，进行观测的一组点（作为单个
ndim x npoints数组）和观测值。
因此，如果我们这样编写函数：
def func(x, y, z, a, b, c): f = (a**2 + x * b**2 + y * a * b * np.cos(c) + z * a * b * np.sin(c)) return f
在将其传递给之前，我们需要将其包装为接受稍有不同的参数
curve_fit。
简而言之：
def nonlinear_invert(f, x, y, z): def wrapped_func(observation_points, a, b, c): x, y, z = observation_points return func(x, y, z, a, b, c) xdata = np.vstack([x, y, z]) model, cov = opt.curve_fit(wrapped_func, xdata, f) return model
这两种方法的独立示例：
为了给您一个完整的实现，下面是一个示例
- 生成随机分布的点以评估函数，
- 在这些点上评估功能（使用设置的模型参数），
- 给结果增加噪音，
- 然后使用上述线性方法和非线性方法对模型参数进行求逆。
import numpy as npimport scipy.optimize as optdef main(): nobservations = 4 a, b, c = 3.0, 2.0, 1.0 f, x, y, z = generate_data(nobservations, a, b, c) print 'Linear results (should be {}, {}, {}):'.format(a, b, c) print linear_invert(f, x, y, z) print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c) print nonlinear_invert(f, x, y, z)def generate_data(nobservations, a, b, c, noise_level=0.01): x, y, z = np.random.random((3, nobservations)) noise = noise_level * np.random.normal(0, noise_level, nobservations) f = func(x, y, z, a, b, c) + noise return f, x, y, zdef func(x, y, z, a, b, c): f = (a**2 + x * b**2 + y * a * b * np.cos(c) + z * a * b * np.sin(c)) return fdef linear_invert(f, x, y, z): G = np.vstack([np.ones_like(x), x, y, z]).T m, _, _, _ = np.linalg.lstsq(G, f) d, e, f, g = m a = np.sqrt(d) b = np.sqrt(e) c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees return a, b, cdef nonlinear_invert(f, x, y, z): # "curve_fit" expects the function to take a slightly different form... def wrapped_func(observation_points, a, b, c): x, y, z = observation_points return func(x, y, z, a, b, c) xdata = np.vstack([x, y, z]) model, cov = opt.curve_fit(wrapped_func, xdata, f) return modelmain()
