最小二乘代数法圆拟合大坐标下结果异常问题求助

问题现象 #

当前采用基于最小二乘法（LSF）的代数圆拟合方案，核心逻辑为构建Ax=B矩阵直接求解，因对该方法推导逻辑掌握清晰，将其作为主用圆拟合方案。
但方案存在明确异常：当输入点坐标数值升高至约500及以上时，拟合结果出现畸变，求解得到的圆心坐标甚至可达$10^8$量级，完全偏离预期范围。
核心实现代码（简化版，仅保留求解逻辑，未附绘图部分）如下：

import numpy as np
import math
import matplotlib as mpl
import matplotlib.pyplot as plt

for num in range(N):
    point = input('Input coordinate of point {}: '.format(num+1))

    # 拆分x、y坐标
    point_x = float(point.split()[0])
    points_x.append(point_x)

    # 计算x相关项求和
    sum_x += point_x
    sum_x_squared += point_x**2
    sum_x_cubed += point_x**3

    # 计算y相关项求和
    point_y = float(point.split()[1])
    points_y.append(point_y)
    sum_y += point_y
    sum_y_squared += point_y**2
    sum_y_cubed += point_y**3

    # 计算交叉项求和
    sum_xy += point_x * point_y
    sum_xsq_y += point_x**2 * point_y
    sum_x_ysq += point_x * point_y**2

    points.append((point_x, point_y))

# 求导后构建矩阵形式方程
# matrix1为系数矩阵，对应A、B、C三个参数的系数
matrix1 = np.array([[sum_x_squared, sum_xy, sum_x],
                   [sum_xy, sum_y_squared, sum_y],
                   [sum_x, sum_y, len(points)]])

# matrix2为常数项矩阵
matrix2 = np.array([[sum_x*(sum_x_squared + sum_y_squared)],
                   [sum_y*(sum_x_squared + sum_y_squared)],
                   [sum_x_squared + sum_y_squared]])

# 求解矩阵
center_x = np.linalg.solve(matrix1, matrix2)[0]/2
center_y = np.linalg.solve(matrix1, matrix2)[1]/2
radius = math.sqrt(np.linalg.solve(matrix1, matrix2)[2] + (center_x)**2 + (center_y)**2)/2

print('Equation of best fit circle: (x-{})^2 + (y-{})^2 = {}^2'.format(center_x, center_y, radius))

本次测试输入为8个坐标值在970~1030区间的点：

Number of points: 8
Input coordinate of point 1: 987.091 1028.443
Input coordinate of point 2: 992.504 1021.98
Input coordinate of point 3: 991.764 1013.589
Input coordinate of point 4: 985.328 1008.178
Input coordinate of point 5: 976.959 1008.871
Input coordinate of point 6: 971.523 1015.297
Input coordinate of point 7: 972.204 1023.683
Input coordinate of point 8: 978.459 1029.105

实际求解得到的拟合圆方程完全偏离预期：

Equation of best fit circle: (x-[1.13510241e+08])^2 + (y-[1.18346499e+08])^2 = 163981720.40047127^2

拟合得到的圆完全落在输入点围成的约束范围外，不符合使用预期。

问题根因 #

• 核心公式实现错误：代码中提前定义了sum_x_cubed、sum_y_cubed、sum_xsq_y、sum_x_ysq四个高次项、交叉项求和变量，但构建右端常数矩阵matrix2时完全没有使用这些变量，错误地用sum_x*(sum_x_squared + sum_y_squared)这类低次项乘积代替了正规方程要求的高次项求和，直接导致求解逻辑错误。
• 数值稳定性缺陷：即便修正公式错误，当输入点坐标绝对值较大时，构建的正规方程矩阵条件数会急剧升高，浮点运算的截断误差会被大幅放大，同样可能出现量级级别的求解偏差。

修复方案 #

1. 修正矩阵构建逻辑
   最小二乘代数圆拟合基于展开后的圆一般方程：$x^2+y2 = Ax + By + C$，其中圆心为$(A/2, B/2)$，半径为$\frac{1}{2}\sqrt{A^2+B2+4C}$。对残差求导得到的正规方程中，左端系数矩阵matrix1的写法是正确的，但右端matrix2的三个元素需要修正为：
   • 第一行：sum_x_cubed + sum_x_ysq（即所有点的$x^3+xy2$求和）
   • 第二行：sum_xsq_y + sum_y_cubed（即所有点的$x^2y+y3$求和）
   • 第三行：sum_x_squared + sum_y_squared（原写法正确，无需修改）
2. 增加坐标中心化预处理提升数值稳定性
   在计算所有求和项前，先对输入点做坐标平移，将点集质心移动到坐标原点，从根源上降低系数矩阵的条件数：
   • 先计算所有输入点的x坐标均值mean_x、y坐标均值mean_y
   • 对每个点计算平移后的局部坐标：x_local = x - mean_x，y_local = y - mean_y
   • 用局部坐标计算所有求和项，求解得到局部坐标系下的圆心(cx_local, cy_local)和半径r（半径不受平移影响）
   • 将圆心坐标转换回原坐标系：cx = cx_local + mean_x，cy = cy_local + mean_y

修正后用本次测试用例计算，得到的拟合圆圆心约为(982.0, 1018.7)，半径约为10.8，完全落在输入点的约束范围内，结果符合预期。

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