Scipy curve_fit求解投影与Helmert变换参数结果不准求助

我正在为某投影求解球面坐标转笛卡尔坐标的最优参数，该投影的数学描述见附图（Helmert变换中省略了半径）。需要求解的参数包括4个投影参数+平面坐标Helmert变换的4个参数，共8个未知量，现有12组球面坐标与对应平面坐标对。

使用scipy.optimize.curve_fit求解时结果不准确：pcov数值极大，角度参数偏差达7度。以下是脚本代码及输出（最优参数值见代码注释）：

import numpy as np
from scipy.optimize import curve_fit

def geographic2ProjAndHelmert(VRadianCoord, lat_0, lon_0, lat_st_1, lat_st_2, a, b, TranslationX, TranslationY):
    #VRadianCoord - [lon1 (x), lat1 (y), lon2, lat2, ...]; a - scale * cos(φ); b - scale * sin(φ)
    
    VRadianCoord=VRadianCoord.reshape((-1, 2))
    
    #scale variables
    a = a * 10
    b = b * 10
    lat_0 = lat_0 / 1000
    lon_0 = lon_0 / 1000
    lat_st_1 = lat_st_1 / 1000
    lat_st_2 = lat_st_2 / 1000
    
    #Calculate Cartesian
    cos1 = np.cos(lat_st_1)
    cos2 = np.cos(lat_st_2)
    CompositeTanCalc0 = np.tan(np.pi/4.0 + lat_0/2.0)
    CompositeTanCalc1 = np.tan(np.pi/4.0 + lat_st_1/2.0)
    CompositeTanCalc2 = np.tan(np.pi/4.0 + lat_st_2/2.0)
    n = (np.log(cos1/cos2)) / (np.log(CompositeTanCalc2/CompositeTanCalc1))
    F = cos1*CompositeTanCalc1**n/n
    r   = F/np.tan(np.pi/4.0 + VRadianCoord[:,1]/2.0)**n
    r_0 = F/CompositeTanCalc0**n
    ProjCoord = np.array([r*np.sin(n*(VRadianCoord[:,0]-lon_0)), r_0-r*np.cos(n*(VRadianCoord[:,0]-lon_0))])
    
    #Calculate Helmert Transformation
    RotationAndScale = np.array([[a, -b], [b, a]])
    EndCoord = np.copy(np.dot(RotationAndScale, ProjCoord).T)
    EndCoord[:,0] = np.add(EndCoord[:,0], TranslationX)
    EndCoord[:,1] = np.add(EndCoord[:,1], TranslationY)
    
    return EndCoord.ravel()

RadCoord = np.array([0.73303829, 1.15191731, 1.01229097, 1.13446401, 0.9250245, 0.97738438, 0.71558499, 1.01229097, 0.87266463, 1.15191731, 0.83775804, 0.99483767, 0.95993109, 1.06465084, 0.71558499, 1.08210414, 0.87266463, 1.06465084, 0.78539816, 1.02974426, 0.80285146, 1.08210414, 0.89011792, 1.01229097])
#format [x1, y1, x2, y2, ...]

Ydata = np.array([897.66512043, -4893.42010515, 585.79116218, -470.03157748, 6683.05857945, -272.07474213, 6149.61630033, -4694.83472251, 550.95105566, -2763.00971689, 6453.46869002, -2176.47012956, 3359.83097409, -508.82821497, 3546.62200096, -4926.08133397, 3761.31345969, -2047.13507678, 5341.79985183, -3415.91999683, 3355.96753336, -3399.9554403, 5595.98729439, -1281.99094122])

p0          = np.radians(55)*1000, np.radians(85)*1000, np.radians(40)*1000, np.radians(69)*1000, -3705.8,  4627, -2887, 7500
#p0_optimal = np.radians(60)*1000, np.radians(90)*1000, np.radians(48)*1000, np.radians(64)*1000, -2705.8,  2627, -1887, 9500

popt, pcov = curve_fit(geographic2ProjAndHelmert, RadCoord, Ydata, p0, maxfev=100000, bounds=([np.radians(30)*1000,np.radians(60)*1000,np.radians(30)*1000, np.radians(40)*1000, -6000, -6000, -6000, -10000],[np.radians(70)*1000,np.radians(100)*1000,np.radians(70)*1000, np.radians(70)*1000, 6000, 6000, 6000, 10000]))
popt[0:4] = np.degrees(popt[0:4])/1000
print("popt = ", popt, "\n")
print("pcov = ", pcov)

输出结果 #

popt =  [   59.20915123    83.65741635    40.57766132    69.72091072
 -2508.47749214  2936.1387793    -73.51463145  8172.85715372]

pcov =  [[ 9.18230672e+12 -2.08906927e+13  1.99021988e+12 -1.51118796e+12
   5.20958799e+13  4.21175110e+13 -1.30500456e+12 -5.27326558e+14]
 [-2.08906927e+13  4.76711630e+13 -5.13324974e+12  3.89771178e+12
  -1.19230289e+14 -9.56985287e+13  1.19399847e+12  1.20179982e+15]
 [ 1.99021988e+12 -5.13324974e+12  6.92449865e+12 -5.25781389e+12
   1.66182697e+13  5.88100103e+12  7.24698818e+12 -1.23109006e+14]
 [-1.51118796e+12  3.89771178e+12 -5.25781389e+12  3.99229002e+12
  -1.26183630e+13 -4.46548828e+12 -5.50267812e+12  9.34775221e+13]
 [ 5.20958799e+13 -1.19230289e+14  1.66182697e+13 -1.26183630e+13
   3.00448606e+14  2.36727395e+14  1.38861694e+12 -3.00208240e+15]
 [ 4.21175110e+13 -9.56985287e+13  5.88100103e+12 -4.46548828e+12
   2.36727395e+14  1.95183572e+14 -7.51770717e+12 -2.41695451e+15]
 [-1.30500456e+12  1.19399847e+12  7.24698818e+12 -5.50267812e+12
   1.38861694e+12 -7.51770717e+12  2.22666466e+13  4.90988220e+13]
 [-5.27326558e+14  1.20179982e+15 -1.23109006e+14  9.34775221e+13
  -3.00208240e+15 -2.41695451e+15  4.90988220e+13  3.03138515e+16]]

可以看到，拟合结果与预期偏差较大，pcov数值异常高，角度参数偏差达7度。

已尝试的解决方法 #

• 变量缩放、调整10的幂次
• 修改sigma参数
• 增减参与拟合的坐标点数量
• 移除Y方向平移参数
• 预处理坐标排除Helmert变换参数：
  1. 将坐标平移至质心
  2. 缩放使点到坐标中心的最大距离为1
  3. 旋转使最远点位于某轴上
     预处理后求解速度提升，但结果仍不准确。

恳请提供技术帮助及替代工具建议。

内容的提问来源于stack exchange，提问作者lukon