我认为以下可能是一个有用的近似值，它随点数线性而不是二次缩放，并且很容易实现：
- 计算点的质心 M
- 找到距离M最大的点P0
- 找到到P0的最大距离的点P1
- 用 P0 和 P1 之间的距离近似最大直径
这可以通过重复步骤 3 N 次来概括，
并取 PN-1 和 PN
之间的距离
步骤 1 可以有效地将 M 近似为经度和纬度的平均值，当距离“小”并且两极距离足够远时，这是可以的。其他步骤可以使用精确的距离公式执行，但如果点的坐标可以近似为位于平面上，它们会更快。一旦找到“远对”（希望是距离最大的对），就可以用精确的公式重新计算它的距离。
一个近似的例子如下：如果 φ(M) 和 λ(M) 是质心的纬度和经度，计算为 Σφ(P)/n 和 Σλ(P)/n，
- x(P) = (λ(P) - λ(M) + C) cos(φ(P))
- y(P) = φ(P) - φ(M) [这只是为了清楚起见，也可以简单地为 y(P) = φ(P)]
其中 C 通常为 0，但如果点集穿过 λ=±180° 线，则可以为 ± 360°。要找到最大距离，您只需找到
- max((x(PN) - x(PN-1))2 + (y(P N) - y(PN-1))2)
（你不需要平方根，因为它是单调的）
相同的坐标变换可用于重复步骤 1（在新坐标系中）以获得更好的起点。我怀疑如果满足某些条件，上述步骤（不重复步骤 3）总是会导致“真正的远距离对”（我的术语）。如果我只知道哪些条件...
编辑：
我讨厌建立在别人的解决方案上，但总得有人这样做。
仍然保持上述 4 个步骤，可选（但可能有益，取决于点的典型分布）重复第 3 步，
并关注solution of Spacedman，
在 3D 中进行计算克服了与两极的接近和距离的限制：
- x(P) = sin(φ(P))
- y(P) = cos(φ(P)) sin(λ(P))
- z(P) = cos(φ(P)) cos(λ(P))
（唯一的近似是这仅适用于完美球体）
质心由x(M) = Σx(P)/n等给出，
并且要寻找的最大值是
- max((x(PN) - x(PN-1))2 + (y(P N) - y(PN-1))2 + (z(PN) - z(PN-1))2)
所以：您首先将球面坐标转换为笛卡尔坐标，然后从质心开始，至少分两步（第 2 步和第 3 步）找到距前一点最远的点。只要距离增加，您就可以重复第 3 步，也许重复次数最多，但这不会使您远离局部最大值。如果这些点遍布整个地球，那么从质心开始也没有多大帮助。
编辑 2：
我学了足够多的 R 来写下算法的核心（数据分析的好语言！）
对于平面近似，忽略λ=±180°线周围的问题：
# input: lng, lat (vectors)
rad = pi / 180;
x = (lng - mean(lng)) * cos(lat * rad)
y = (lat - mean(lat))
i = which.max((x - mean(x))^2 + (y )^2)
j = which.max((x - x[i] )^2 + (y - y[i])^2)
# output: i, j (indices)
在我的 PC 上，查找索引 i 和 j 所需的时间不到 1000000 个点。
下面的 3D 版本有点慢，但适用于任何点分布（并且不适用于穿越λ=±180°线时需要修正）：
# input: lng, lat
rad = pi / 180
x = sin(lat * rad)
f = cos(lat * rad)
y = sin(lng * rad) * f
z = cos(lng * rad) * f
i = which.max((x - mean(x))^2 + (y - mean(y))^2 + (z - mean(z))^2)
j = which.max((x - x[i] )^2 + (y - y[i] )^2 + (z - z[i] )^2)
k = which.max((x - x[j] )^2 + (y - y[j] )^2 + (z - z[j] )^2) # optional
# output: j, k (or i, j)
k 的计算可以省略（即，结果可以由i 和j 给出），具体取决于数据和要求。另一方面，我的实验表明，再计算一个索引是没有用的。
应该记住，在任何情况下，结果点之间的距离都是估计值，它是集合“直径”的下限，尽管它通常是直径本身（如何 通常取决于数据。）
编辑 3：
不幸的是，平面近似的相对误差在极端情况下可能高达 1-1/√3 ≅ 42.3%，即使非常罕见，这也可能是不可接受的。可以修改该算法以获得大约 20% 的上限，这是我通过罗盘和直尺得出的（解析解很麻烦）。修改后的算法找到一对具有局部最大距离的点，然后重复相同的步骤，但这次从第一对的中点开始，可能会找到不同的对：
# input: lng, lat
rad = pi / 180
x = (lng - mean(lng)) * cos(lat * rad)
y = (lat - mean(lat))
i.n_1 = 1 # n_1: n-1
x.n_1 = mean(x)
y.n_1 = 0 # = mean(y)
s.n_1 = 0 # s: square of distance
repeat {
s = (x - x.n_1)^2 + (y - y.n_1)^2
i.n = which.max(s)
x.n = x[i.n]
y.n = y[i.n]
s.n = s[i.n]
if (s.n <= s.n_1) break
i.n_1 = i.n
x.n_1 = x.n
y.n_1 = y.n
s.n_1 = s.n
}
i.m_1 = 1
x.m_1 = (x.n + x.n_1) / 2
y.m_1 = (y.n + y.n_1) / 2
s.m_1 = 0
m_ok = TRUE
repeat {
s = (x - x.m_1)^2 + (y - y.m_1)^2
i.m = which.max(s)
if (i.m == i.n || i.m == i.n_1) { m_ok = FALSE; break }
x.m = x[i.m]
y.m = y[i.m]
s.m = s[i.m]
if (s.m <= s.m_1) break
i.m_1 = i.m
x.m_1 = x.m
y.m_1 = y.m
s.m_1 = s.m
}
if (m_ok && s.m > s.n) {
i = i.m
j = i.m_1
} else {
i = i.n
j = i.n_1
}
# output: i, j
可以用类似的方式修改 3D 算法。可以（在 2D 和 3D 情况下）从第二对点的中点（如果找到）重新开始。在这种情况下，上限是“留给读者练习”:-)。
修改后的算法与（过于）简单的算法的比较表明，对于正态分布和方形均匀分布，处理时间几乎翻倍，平均误差从 0.6% 降低到 0.03%（顺序量级）。从中点进一步重新开始会导致平均误差稍好一些，但几乎等于最大误差。
编辑 4：
我还得研究this article，但看起来我用指南针和直尺找到的 20% 实际上是 1-1/√(5-2√3) ≅ 19.3%