通过轨道数据拟合椭圆

Question

我已经为围绕太阳运行的行星的 (x,y,z) 坐标生成了一堆数据。现在我想通过这些数据拟合一个椭圆。

我尝试做的事情：

我根据五个参数创建了一个虚拟椭圆：定义大小和形状的半长轴和偏心率以及围绕椭圆旋转的三个欧拉角。由于我的数据并不总是以原点为中心，我还需要平移需要额外三个变量 (dx、dy、dz) 的椭圆。一旦我用这八个变量初始化这个函数，我就会得到 N 个位于这个椭圆上的点。（N = 我绘制椭圆的数据点数）我计算这些虚拟点与实际数据的偏差，然后使用 some 最小化方法最小化此偏差，以找到这八个变量的最佳拟合值。

我的问题是最后一部分：最小化偏差并找到变量的值。

为了尽量减少偏差，我使用 scipy.optimize.minimize 来尝试近似最佳拟合变量，但它做得不够好：

Here is an image 我最适合的一个看起来像什么，而且初始猜测非常准确。（蓝色 = 数据，红色 = 适合）

Here is the entire code.（不需要数据，它会生成自己的虚假数据）

简而言之，我用这个scipy函数：

initial_guess = [0.3,0.2,0.1,0.7,3,0.0,-0.1,0.0]
bnds = ((0.2, 0.5), (0.1, 0.3), (0, 2*np.pi), (0, 2*np.pi), (0, 2*np.pi), (-0.5,0.5), (-0.5,0.5), (-0.3,0.3)) #reasonable bounds for the variables
result = optimize.minimize(deviation, initial_guess, args=(data,), method='L-BFGS-B', bounds=bnds, tol=1e-8) #perform minimalisation
semi_major,eccentricity,inclination,periapsis,longitude,dx,dy,dz = result["x"]

要最小化这个误差（或偏差）函数：

    def deviation(variables, data):
    """
    This function calculates the cumulative seperation between the ellipse fit points and data points and returns it
    """
    num_pts = len(data[:,0])
    semi_major,eccentricity,inclination,periapsis,longitude,dx,dy,dz = variables
    dummy_ellipse = generate_ellipse(num_pts,semi_major,eccentricity,inclination,periapsis,longitude,dz,dy,dz)
    deviations = np.zeros(len(data[:,0]))
    pair_deviations = np.zeros(len(data[:,0]))
    # Calculate separation between each pair of points
    for j in range(len(data[:,0])):
        for i in range(len(data[:,0])):
            pair_deviations[i] = np.sqrt((data[j,0]-dummy_ellipse[i,0])**2 + (data[j,1]-dummy_ellipse[i,1])**2 + (data[j,2]-dummy_ellipse[i,2])**2)
            deviations[j] = min(pair_deviations) # only pick the closest point to the data point j.
    total_deviation = sum(deviations)
    return total_deviation

（我的代码可能有点混乱和低效，我是新手）

我可能在编码时犯了一些逻辑错误，但我认为这归结为 scipy.minimize.optimize 函数。我不太了解它是如何工作的以及对它有什么期望。我还被推荐在处理这么多变量时尝试马尔可夫链 Monte Carlo。我确实看过 emcee，但现在有点超出我的理解范围了。

Answer 1

首先，您的 objective 函数中有一个拼写错误，导致无法优化其中一个变量：

dummy_ellipse = generate_ellipse(...,dz,dy,dz)

应该是

dummy_ellipse = generate_ellipse(...,dx,dy,dz)

另外，去掉 sqrt 并最小化欧几里德距离的平方和使得优化器在数值上更容易一些。

您的 objective 函数也并非处处可微，因为 min()，正如 BFGS 求解器所假设的那样，因此其性能将不是最优的。

此外，从解析几何的角度处理问题可能会有所帮助：3d 中的椭圆定义为两个方程的解

f1(x,y,z,p) = 0
f2(x,y,z,p) = 0

其中p是椭圆的参数。现在，为了使参数适合数据集，您可以尝试最小化

F(p) = sum_{j=1}^N [f1(x_j,y_j,z_j,p)**2 + f2(x_j,y_j,z_j,p)**2]

总和超过数据点。

更好的是，在这个问题公式中你可以使用 optimize.leastsq，这在最小二乘问题中可能更有效。

通过轨道数据拟合椭圆

Fitting an ellipse through orbital data

python

ellipse

minimization

astronomy

scipy