10560 怎樣在球面上「均勻」排列許多點？（上）

01-27

　　在這個春不暖花未開的時節，我又跟一道數學題杠上了：怎樣在球面上「均勻」地排列許多點呢？

　　這是一個很有實際意義的問題。比如我要測量地球上陸地的總面積。如果能在地球表面「均勻」地取 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}nz_n &= (2n-1)/N - 1 & (1) nx_n &= sqrt{1 - z_n^2} cdot cos(2 pi n phi) & (2) ny_n &= sqrt{1 - z_n^2} cdot sin(2 pi n phi) & (3)nend{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$ 處，環面的半徑為 $costheta$ ，而環面的寬度為 $2/(Ncostheta)$ （思考：為什麼要除以 $costheta$ ？），故各環面的面積均為 $4pi/N$ 。這個性質保證了點陣分布在宏觀上的均勻性：不管在什麼緯度，都是每 $4pi/N$ 面積上有一個點。

　　第二、三條公式 $begin{align}nx_n &= sqrt{1 - z_n^2} cdot cos(2 pi n phi) ny_n &= sqrt{1 - z_n^2} cdot sin(2 pi n phi)nend{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}1904dot{7}$ 。若每產生一個點轉 0.619 圈，那麼產生 21 個點一共就轉過了差不多正好 13 圈。第 n 個點和第 n + 21 個點的經度十分接近，這就是條紋的來源。再看 $phi = 0.616$ 時的點陣，它形成 13 條螺旋。這是因為 $0.616 approx 8/13 = 0.dot{6}1538dot{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$ 也不錯（見下圖對比），它比黃金分割比差在哪裡？

這幾個問題，我將在下一篇專欄中詳細分析。

　　順便說一下，在平面上也可以討論怎樣生成「均勻、密集、混亂」的點陣。一個與公式 (1~3) 類似的解答如下：

　　　 $begin{align}nx_n &= sqrt{n} cos(2pi n phi) & (4)ny_n &= sqrt{n} sin(2pi n phi) & (5)nend{align}$

即第 n 個點落在半徑為 $sqrt{n}$ 的圓上，連續兩點之間也是轉過一個黃金角。半徑 $sqrt{n}$ 中的根號保證了均勻（證明留給讀者），黃金角保證了密集和混亂。(4,5) 兩個公式描述的排列規律，在自然界中很常見，一個例子就是題圖中向日葵的花序。

　　與公式 (1~3) 不同，(4,5) 兩式可以生成無窮多個點，它們排成一個越來越大的圓盤。這種無限使得 (4,5) 兩式容易揭示更多規律，所以我下一篇中的分析將從平面點陣的分布入手。