求基于三个球面点计算球心的Python/PowerShell实现方案

需求说明 #

输入三个三维空间中的球面点（示例格式如下），编写脚本计算并输出球心坐标：

$p1 = 121021.7391,27709572.4368,-2369575.3890
$p2 = 746442.4255,27095391.7051,-2388765.7650
$p3 = 766543.7492,27435797.7212,-2571188.0087

现有基础与疑问 #

• 擅长PowerShell，但不确定所需的向量运算、线性方程组求解等数学函数是否存在；
• 对Python了解不多；
• 已实现三角形重心计算，但不知道如何求平面垂线及等距点，询问是否有现成库或函数可用。

解决方案 #

PowerShell实现（手动实现核心运算） #

PowerShell没有内置线性代数库，但可以基于.NET基础数学方法，手动实现向量运算和线性方程组求解。核心思路是利用球心到三个点距离相等的条件，推导线性方程组后用克莱姆法则求解。

# 定义三个点（转为浮点数组）
$p1 = [double[]]@(121021.7391,27709572.4368,-2369575.3890)
$p2 = [double[]]@(746442.4255,27095391.7051,-2388765.7650)
$p3 = [double[]]@(766543.7492,27435797.7212,-2571188.0087)

# 推导前两个线性方程（距离相等展开）
$a1 = 2*($p2[0]-$p1[0])
$b1 = 2*($p2[1]-$p1[1])
$c1 = 2*($p2[2]-$p1[2])
$d1 = ($p2[0]**2 + $p2[1]**2 + $p2[2]**2) - ($p1[0]**2 + $p1[1]**2 + $p1[2]**2)

$a2 = 2*($p3[0]-$p1[0])
$b2 = 2*($p3[1]-$p1[1])
$c2 = 2*($p3[2]-$p1[2])
$d2 = ($p3[0]**2 + $p3[1]**2 + $p3[2]**2) - ($p1[0]**2 + $p1[1]**2 + $p1[2]**2)

# 计算平面法向量（p2-p1叉乘p3-p1），用于第三个方程（球心在平面垂线上）
$cross_x = ($p2[1]-$p1[1])*($p3[2]-$p1[2]) - ($p2[2]-$p1[2])*($p3[1]-$p1[1])
$cross_y = ($p2[2]-$p1[2])*($p3[0]-$p1[0]) - ($p2[0]-$p1[0])*($p3[2]-$p1[2])
$cross_z = ($p2[0]-$p1[0])*($p3[1]-$p1[1]) - ($p2[1]-$p1[1])*($p3[0]-$p1[0])
$d3 = $cross_x*$p1[0] + $cross_y*$p1[1] + $cross_z*$p1[2]

# 用克莱姆法则求解3阶线性方程组
# 计算系数矩阵行列式
$det = $a1*($b2*$cross_z - $c2*$cross_y) - $b1*($a2*$cross_z - $c2*$cross_x) + $c1*($a2*$cross_y - $b2*$cross_x)

# 替换第一列，计算x的行列式
$mat_x = @(
    @($d1, $b1, $c1),
    @($d2, $b2, $c2),
    @($d3, $cross_y, $cross_z)
)
$det_x = $mat_x[0][0]*($mat_x[1][1]*$mat_x[2][2] - $mat_x[1][2]*$mat_x[2][1]) - `
         $mat_x[0][1]*($mat_x[1][0]*$mat_x[2][2] - $mat_x[1][2]*$mat_x[2][0]) + `
         $mat_x[0][2]*($mat_x[1][0]*$mat_x[2][1] - $mat_x[1][1]*$mat_x[2][0])

# 替换第二列，计算y的行列式
$mat_y = @(
    @($a1, $d1, $c1),
    @($a2, $d2, $c2),
    @($cross_x, $d3, $cross_z)
)
$det_y = $mat_y[0][0]*($mat_y[1][1]*$mat_y[2][2] - $mat_y[1][2]*$mat_y[2][1]) - `
         $mat_y[0][1]*($mat_y[1][0]*$mat_y[2][2] - $mat_y[1][2]*$mat_y[2][0]) + `
         $mat_y[0][2]*($mat_y[1][0]*$mat_y[2][1] - $mat_y[1][1]*$mat_y[2][0])

# 替换第三列，计算z的行列式
$mat_z = @(
    @($a1, $b1, $d1),
    @($a2, $b2, $d2),
    @($cross_x, $cross_y, $d3)
)
$det_z = $mat_z[0][0]*($mat_z[1][1]*$mat_z[2][2] - $mat_z[1][2]*$mat_z[2][1]) - `
         $mat_z[0][1]*($mat_z[1][0]*$mat_z[2][2] - $mat_z[1][2]*$mat_z[2][0]) + `
         $mat_z[0][2]*($mat_z[1][0]*$mat_z[2][1] - $mat_z[1][1]*$mat_z[2][0])

# 计算球心坐标
$ox = $det_x / $det
$oy = $det_y / $det
$oz = $det_z / $det

Write-Host "球心坐标：($($ox.ToString("F4")), $($oy.ToString("F4")), $($oz.ToString("F4")))"

Python实现（利用numpy库简化运算） #

Python的numpy库提供了成熟的线性代数工具，无需手动实现矩阵运算，直接构建方程组求解即可。

import numpy as np

# 定义三个点
p1 = np.array([121021.7391, 27709572.4368, -2369575.3890])
p2 = np.array([746442.4255, 27095391.7051, -2388765.7650])
p3 = np.array([766543.7492, 27435797.7212, -2571188.0087])

# 构建线性方程组的系数矩阵A和常数项向量b
A = 2 * np.array([p2 - p1, p3 - p1, np.cross(p2 - p1, p3 - p1)])
b = np.array([
    np.dot(p2, p2) - np.dot(p1, p1),
    np.dot(p3, p3) - np.dot(p1, p1),
    np.dot(np.cross(p2 - p1, p3 - p1), p1)
])

# 求解方程组
center = np.linalg.solve(A, b)

print(f"球心坐标：({center[0]:.4f}, {center[1]:.4f}, {center[2]:.4f})")

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