python - 将两层模型拟合到 python 中的风剖面数据

Question

我正在尝试将模型拟合到我的风剖面数据集，即不同高度 z 的风速值 u(z)。

该模型由两部分组成，我现在将其简化为:

u(z) = ust/k * ln(z/z0)     for z < zsl

u(z) = a*z + b              for z > zsl

在对数模型中，ust和z0是自由参数，k是固定的。 zsl 是表层的高度，也是先验未知的。

我想使这个模型适合我的数据，并且我已经尝试过不同的方法。到目前为止我得到的最好结果是:

def two_layer(z,hsl,ust,z0,a,b):
    return ust/0.4*(np.log(z/z0)) if z<hsl else a*z+b

two_layer = np.vectorize(two_layer)

def two_layer_err(p,z,u):
    return two_layer(z,*p)-u

popt, pcov ,infodict, mesg, ier = optimize.leastsq(two_layer_err,[150.,0.3,0.002,0.1,10.],(wspd_hgt,prof),full_output=1)

# wspd_hgt are my measurements heights and 
# prof are the corresponding wind speed values

这为我提供了所有参数的合理估计，zsl 除外，它在拟合过程中未更改。我想这与用作阈值而不是函数参数的事实有关。有什么方法可以让 zsl 在优化过程中发生变化？

我用 numpy.piecewise 尝试了一些东西，但效果不是很好，也许是因为我不太了解它，或者我可能完全离开这里因为它不适合我的事业。

对于这个想法，如果轴反转(z 相对于 u 绘制)，风剖面看起来像这样:

Answer 1

我想我终于找到了解决此类问题的方法，这是我在回答 similar question 时遇到的。 .

解决方案似乎是在两个模型之间的切换处实现一个约束，说明 u1 == u2。由于没有发布数据，我无法在您的模型上进行尝试，因此我将展示它如何适用于其他模型，您可以根据自己的情况进行调整。我使用一个 scipy 包装器解决了这个问题，我编写了一个名为 symfit 的包装器，使此类问题的拟合更像 pythonic。但如果您愿意，也可以使用 scipy 中的 SLSQP 算法来执行相同的操作。

from symfit import parameters, variables, Fit, Piecewise, exp, Eq
import numpy as np
import matplotlib.pyplot as plt

t, y = variables('t, y')
m, c, d, k, t0 = parameters('m, c, d, k, t0')

# Help the fit by bounding the switchpoint between the models
t0.min = 0.6
t0.max = 0.9

# Make a piecewise model
y1 = m * t + c
y2 = d * exp(- k * t)
model = {y: Piecewise((y1, t <= t0), (y2, t > t0))}

# As a constraint, we demand equality between the two models at the point t0
# Substitutes t in the components by t0
constraints = [Eq(y1.subs({t: t0}), y2.subs({t: t0}))]

# Read the data
tdata, ydata = np.genfromtxt('Experimental Data.csv', delimiter=',', skip_header=1).T

fit = Fit(model, t=tdata, y=ydata, constraints=constraints)
fit_result = fit.execute()
print(fit_result)

plt.scatter(tdata, ydata)
plt.plot(tdata, fit.model(t=tdata, **fit_result.params).y)
plt.show()

我认为您应该能够根据您的情况调整此示例!

编辑:根据评论中的要求，也可以要求匹配的衍生品。为此，上面的示例需要以下附加代码:

from symfit import Derivative

dy1dt = Derivative(y1, t)
dy2dt = Derivative(y2, t)
constraints = [
    Eq(y1.subs({t: t0}), y2.subs({t: t0})),
    Eq(dy1dt.subs({t: t0}), dy2dt.subs({t: t0}))
]

这应该可以解决问题!因此，从编程的角度来看，这是非常可行的，尽管根据模型的不同，这实际上可能没有解决方案。