10560 怎样在球面上「均匀」排列许多点？（上）

2017-03-25 13:38:59

　　在这个春不暖花未开的时节，我又跟一道数学题杠上了：怎样在球面上「均匀」地排列许多点呢？

　　这是一个很有实际意义的问题。比如我要测量地球上陆地的总面积。如果能在地球表面「均匀」地取 n 个点，那么我只要简单地数一下其中落在陆地上的点的个数 m，就可以知道陆地面积约占地球总表面积的 m/n 了。注意，按照一定的经纬度间隔，取经纬线的交点是不行的，因为这样取出的点不均匀：两极附近的点比赤道处更密集，如下图所示。

　　要严格解决这个问题，首先要把直观感觉上的「均匀」用数学语言定义出来。一种定义方法是：让各点之间距离的小值大。这样定义出来的问题叫 Tammes problem，是密铺问题（Packing problems）的一个特例。很不幸的是，密铺问题往往没有很优雅的解。另一种定义方法是：把各个点看成同种电荷，让整个系统的电势能小。这个问题可以通过模拟电荷的运动来解决，但计算复杂度非常高，而且只能得到数值解。

　　经过一番搜索，我在 StackOverflow 上找到了一种令我惊艳的近似解法：Evenly distributing n points on a sphere。这种解法只用几个简单的公式，就给出了每个点的坐标。设球面半径为 1，一共要取 N 个点，则第 n 个点的坐标 (x_n, y_n, z_n) 由下列公式给出：

　　★　 $\begin{align} z_n &= (2n-1)/N - 1 & (1) \\ x_n &= \sqrt{1 - z_n^2} \cdot \cos(2 \pi n \phi) & (2) \\ y_n &= \sqrt{1 - z_n^2} \cdot \sin(2 \pi n \phi) & (3) \end{align}$

其中常数 $\phi = (\sqrt{5} - 1) / 2 \approx 0.618$ 正是黄金分割比。

　　用这套公式生成 1000 个点的效果如下，它惊人地符合我们对「均匀」的预期：

　　公式 (1~3) 生成的点阵，称为「菲波那契网格」（Fibonacci lattice 或 Fibonacci grid）。至于为什么叫菲波那契网格，目前可以简单地用「黄金分割比也出现在菲波那契数列中」来解释，在下文中你还会见到菲波那契数的出现。这篇文章（Measurement of areas on a sphere using Fibonacci and latitude–longitude lattices）说明，使用菲波那契网格测量球面上不规则图形的面积，与用经纬网格（并加权）相比，误差可以减小 40%。

　　仔细观察上图中点的分布，可以发现它其实有好几条特点：

均匀：宏观上看，各处点的密度都差不多；
密集：各点之间没有较大的缝隙，以容纳更多的点而不减小点的间距；
混乱：点的排列似乎有规律，但具体是什么规律又说不出来。

我们在说「均匀」的时候，其实是暗含了这三条性质的。

　　公式 (1~3) 对于 $\phi$ 的值非常敏感。哪怕 $\phi$ 稍微偏离黄金分割比一丁点儿，作出的图效果就不好。例如，下面三个图分别是 $\phi$ 取 0.616、0.617 和 0.619 时得到的点阵。可以看到， $\phi = 0.617$ 时点阵在两极处形成螺旋，在赤道附近形成条纹，不满足「混乱」；而左、右图两中的点不仅形成螺旋或条纹，而且间距还很大，连「密集」都不满足。

　　一个很自然的问题就是：为什么公式 (1~3) 有这么神奇的效果呢？我们先来看看它到底做了什么。

　　条公式 z_n = (2n-1)/N-1 ，说明各个点的竖坐标成等差数列。这相当于把球面切成相同厚度的 N 层，并在每一层的厚度中点处的表面上取一个点。注意，各层的厚度相同，但纬度的跨度是不同的：两极处纬度的跨度更大。这样切出来的各层有一个性质：侧面积都相等。这是因为各层的侧面可以近似看成环面，在纬度为 $\theta$ 处，环面的半径为 $\cos\theta$ ，而环面的宽度为 $2/(N\cos\theta)$ （思考：为什么要除以 $\cos\theta$ ？），故各环面的面积均为 $4\pi/N$ 。这个性质保证了点阵分布在宏观上的均匀性：不管在什么纬度，都是每 $4\pi/N$ 面积上有一个点。

　　第二、三条公式 $\begin{align} x_n &= \sqrt{1 - z_n^2} \cdot \cos(2 \pi n \phi) \\ y_n &= \sqrt{1 - z_n^2} \cdot \sin(2 \pi n \phi) \end{align}$ ，实际上就是指明了各个点的经度成等差数列。也可以这样形象地理解：从每一个点到下一个点，首先沿着经线向上爬，使得竖坐标增加 2/N ；然后沿着纬线转 $\phi$ 圈。当然， $\phi \approx 0.618$ 比半圈大，只绕 $1 - \phi \approx 0.382$ 圈也是可以的，如下图所示。这两个值对应的角度分别为 222.5 度和 137.5 度，它们是把 360 度黄金分割后的两部分，称为「黄金角」（Golden angle）。

（图片取自：Distributing Many Points on a Sphere）

　　正是 $\phi = (\sqrt{5} - 1) / 2 \approx 0.618$ 这个常数让产生的点阵具有「密集」和「混乱」两条性质。这是为什么呢？我们观察前文中 $\phi = 0.619$ 时的点阵，它明显形成几乎沿经线方向的条纹，仔细数一下的话，这样的条纹共有 21 条。原来，0.619 十分接近于分数 13/21（注意啦，13 和 21 都是菲波那契数哦），后者的值为 $0.\dot{6}1904\dot{7}$ 。若每产生一个点转 0.619 圈，那么产生 21 个点一共就转过了差不多正好 13 圈。第 n 个点和第 n + 21 个点的经度十分接近，这就是条纹的来源。再看 $\phi = 0.616$ 时的点阵，它形成 13 条螺旋。这是因为 $0.616 \approx 8/13 = 0.\dot{6}1538\dot{4}$ （注意，8 和 13 也是菲波那契数），但误差比 0.619 和 13/21 大一些，所以条纹旋转了起来。从这两个例子我们可以发现，当 $\phi$ 接近分数 p/q 时，就会形成 q 条条纹或螺旋。要满足「密集」和「混乱」，就要选取不接近任何有理数的 $\phi$ 值。

　　网上有许多科普资料，说明「黄金分割比是难用有理数逼近的无理数」，或称「无理的无理数」，所以黄金角是能使得点阵密集、混乱的转角。说明方法一般是把黄金分割比 $\phi$ 写成连分式（Continued fraction）：

要用有理数逼近 $\phi$ ，可以在连分式的任意一层截断，比如舍弃上式中省略号的部分。如果某一层分母的整数部分比较大，那么在这一层舍弃小数部分，误差就比较小。而黄金分割比的连分式中，每一层的分母都是 1，截断误差相对就比较大，故难用有理数逼近。

　　但这种解释，我总觉得不够直接。比如，我还想弄清下列问题：

当连分式中出现较大的项时，点阵会呈现怎样的分布？
为什么当 $\phi = 0.617$ 时，点阵在两极和赤道处呈现出不同的有规律分布，而 $\phi$ 取黄金分割比时则呈现一片混乱？
$\phi = (\sqrt{5} - 1) / 2 \approx 0.618$ 真的是优值吗？我在试探的过程中发现 $\phi = \sqrt{2} - 1 \approx 0.414$ 也不错（见下图对比），它比黄金分割比差在哪里？