球谐函数的基础数学理论

作者: 离原春草 | 来源:发表于2021-08-15 23:48 被阅读0次

球谐函数在图形学光照计算等领域有着重要应用，因为目前在实际工作中接触较少，所以对其的理解仅仅停留在表面，本着越是基础的东西，其重要性越高的想法，特此开篇文章对其背后的数学理论进行拆解，拆解过程参考了大量其他同学的工作，相应链接在文末的参考文献中有列出，引用过程中如有表述不清晰的内容，可以通过原文辅助阅读。

1. 调和函数（harmonic function）

调和函数指的是一种特殊的二阶连续可导函数（简称C2，在某个定义域存在二阶导数，且二阶导数连续），数学符号用 $f:U\rightarrow R$ 表达，其中 $U$ 是 $R^n$ （表示n维实数域）的一个开子集（相当于一维数据空间中的开区间），其特殊在于需要满足拉普拉斯方程（下面有介绍），用（笛卡尔坐标系下）数学表达式来描述的话，就是对于任意 $(x_1, x_2, ..., x_n) \in U$ ，需要满足如下的二阶偏微分方程：
$\frac{\partial^2f}{\partial x_1^2} + \frac{\partial^2f}{\partial x_2^2} + ... + \frac{\partial^2f}{\partial x_n^2} = 0$

这里来回顾一下微分方程的相关知识，单个变量下，也就是一元变量情况下，函数与函数各阶导数组成的微分方程叫做常微分方程：
$F(x, y, y^{'}, y^{''}, ... ,y^{(n)}) = 0$

多元函数而言，函数以及函数对各个自变量的各阶偏导数组成的微分方程叫做偏微分方程：
$F(x, y, u_x^{'}, u_y^{'}, u_x^{''}, u_y^{''}, u_{xy}^{''}, ... ) = 0$

这个公式也经常以如下的形式出现（ $\Delta$ 称为拉普拉斯算子， $\nabla$ 称为向量微分算子，也就是nabla算子）：
$\Delta f = \nabla^2f =div \cdot grad ~ u = 0$

其中 $\Delta$ 叫做拉普拉斯算子，光看定义太抽象，我们来举个例子吧，下面两个函数都是二元的调和函数：
$f(x,y) = ln(x^2+y^2), ~ f(x, y) = e^xsin(y)$

2. 拉普拉斯方程（Laplace's equation）

拉普拉斯方程也被称为调和方程、位势方程，这是一种偏微分方程，因为其可以用势函数的形式来描述电磁场、引力场、流场（统称为保守场或者有势场）的性质而被广泛应用。

笛卡尔坐标系下的表述形式前面已经写过了，下面给出球面坐标系下的拉普拉斯方程表述形式：
$\Delta f = \frac{1}{r^2}\frac{\partial (r^2\frac{\partial f}{\partial r})}{\partial r}+\frac{1}{r^2 sin\theta}\frac{\partial (sin\theta\frac{\partial f}{\partial \theta})}{\partial \theta}+\frac{1}{r^2 sin^2\theta}\frac{\partial ^2f}{\partial\phi^2} = 0$
这个方程也常用如下的简化形式来代替：
$\triangledown^2 \phi = 0$
或者
$div ~ grad ~ \phi = 0$
其中div指的是向量场（指的是空间中的每一点都有一个对应的带长度的向量）的散度（divergence），grad表示的是标量场的梯度（gradient）。

散度是向量分析中常用的向量算子，用于实现向量场到标量场的转换映射，也就是说，经过散度算子处理后，得到的是一个标量场（每一点有一个不带方向的数值）。以静电场为例，空间中的电场强度是一个向量场，电场线正出负归，在正电荷附近，对应的散度为正值，且电荷带电量越大，散度越大，负电荷附近则反之，其散度为负值，且电荷带电量越大，散度绝对值越大。更为通用的概括是，散度可以看成是向量场在某一点的通量密度，当散度大于0的时候，就表示该点有流量留出，此时这一点可以被称为源点，当散度小于0的时候，表示此点有流量流入，此时此点被称为汇点，散度为0，表示该点无流入也无流出，如果整个向量场的散度都是0，那么这个向量场可以称为无源场。

对于某个向量场 $F = P(x,y,z) \cdot i + Q(x,y,z) \cdot j + R(x,y,z) \cdot k$ 而言，其散度可以通过如下公式求得：
$div ~ F = \frac{\partial P}{\partial x} + \frac{\partial Q}{\partial y} + \frac{\partial R}{\partial z}$

梯度是对多元函数的导数的一种描述，单元函数（只有一个自变量）的导数是标量值函数，而多元（多个自变量）函数的导数则是一个向量值函数，这里多元函数的导数，我们也称为多元函数的梯度，多元函数f在点P处的梯度指的是以f在P处的偏微分作为分量的向量，如一个三维空间函数 $f(x,y,z)$ ，其梯度函数可以用如下的形式来表述：
$\triangledown f = (\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z})$

单元函数的导数对应的是函数在某一点切线的斜率，对应到梯度上，如果多元函数在某点P的梯度不为0的话，那么计算出来的梯度方向指的是这个函数在P点处增长最快的方向（超平面的切线），而梯度的长度则是函数在此点处的增长率（超平面的斜率）。

举个例子，如果某个房间内的温度用一个函数来表示，那么这个函数在三维空间中的梯度就对应于房间中某点处温度上升最快的方向，而其长度则对应于温度增长率。

可以看到，一个多元函数的标量场，经过梯度转化后，得到的是一个向量场。

3. 球谐函数(Spherical Harmonics Function)

从调和函数的定义我们可以看到，所谓的调和函数，实际上就是拉普拉斯方程的解，而我们日常所说的球谐函数（Spherical Harmonics Function）实际上就是拉普拉斯方程在球坐标系空间下的解。

拉普拉斯方程是一个偏微分方程，而解偏微分方程常用的策略是分离变量法，即将偏微分方程分解成几个常微分方程进行求解，下面我们通过将半径跟角度进行分离来进行求解。

设 $f(r, \theta, \phi) = R(r)Y(\theta, \phi)$ ，将之代入前面的拉普拉斯方程，可以得到：
$\frac{Y}{r^2} \frac{d}{\text{d}r}(r^2 \frac{\text{d}R}{\text{d}r}) + \frac{R}{r^2 sin \theta}\frac{\partial}{\partial \theta}(sin \theta \frac{\partial Y}{\partial \theta})+\frac{R}{r^2 sin^2 \theta}\frac{\partial ^2 Y}{\partial ^2 \phi} = 0$

上面公式乘上 $\frac{r^2}{YR}$ 之后可以得到：
$\frac{1}{R} \frac{d}{\text{d}r}(r^2 \frac{\text{d}R}{\text{d}r}) = - \frac{1}{Y sin \theta}\frac{\partial}{\partial \theta}(sin \theta \frac{\partial Y}{\partial \theta}) - \frac{1}{Y sin^2 \theta}\frac{\partial ^2 Y}{\partial ^2 \phi} = \lambda$

对于上面公式中后面的等式，我们继续使用分离变量法，令（这里是假设Y具有可以分离的形式，当然这个假设不是必然成立的，只是为了简化计算而给出的，只有一些特殊的函数才具有这种假设的可分离的形式） $Y(\theta, \phi) = \Theta(\theta) \cdot \Phi(\phi)$ ，代入前面公式可以得到：
$-\frac{\Theta}{\Phi \Theta}\frac{\partial ^2 \Phi}{\partial ^2 \phi} = \lambda \cdot sin^2 \theta + \frac{sin \theta \cdot \Phi}{\Phi \Theta}\frac{\partial}{\partial \theta}(sin \theta \frac{\partial \Theta}{\partial \theta})$

简化后，令左右两边均等于 $m^2$ ，可以得到：
$\frac{1}{\Phi}\frac{\partial ^2 \Phi}{\partial ^2 \phi} = -m^2$
$\lambda \cdot sin^2 \theta + \frac{sin \theta}{\Theta}\frac{\partial}{\partial \theta}(sin \theta \frac{\partial \Theta}{\partial \theta}) = m^2$

一个先验知识是m是一个复数常量（怎么得到的？），且由于 $\Phi$ 是一个周期函数，其周期可以整除 $2\pi$ ，因此m就会是一个整数，而 $\Phi$ 则是复数指数 $e^{\pm im \phi}$ 的线性组合，Y的常规解出现在极点，也就是 $\theta = 0, \pi$ 的时候，而在上面的第二个方程中求解 $\Theta$ 时的常规状态出现在Sturm-Liouville problem的边界点上，在这个边界点中会将 $\lambda = l(l+1)$ ，其中l是非负整数，且 $l \geq |m|$ ，此外，将上面公式中的 $cos \theta$ 用t来替代，就能够得到勒让德公式（Legendre equation），而勒让德公式的解就是伴随勒让德多项式 $P_l^m(cos \theta )$ 的倍数。

对于满足前面假设的Y，对于给定的 $l$ ，我们总共有 $2l+1$ 个独立解，这些角度上的解可以表示为三角函数的乘积，这里可以用复数指数与伴随勒让德多项式来表示：
$Y_l^m(\theta, \phi) = N e^{im \phi}P_l^m(cos \theta)$
其中这个解需要满足：
$r^2 \nabla^2Y_l^m(\theta, \phi) = -l(l+1)Y_l^m(\theta, \phi)$
上述公式中的 $Y_l^m$ 就被称为一个m阶（order）l度（degree）的球谐函数， $P_l^m$ 就是一个伴随勒让德多项式，N是一个归一化的常量， $\theta, \phi$ 则代表着球上的经纬度