实值球谐函数从原理到代码：摆脱复数，玩转球谐光照

做球谐光照的时候，我第一次在代码里撞见“球谐函数”这四个字，满屏的复数差点把我劝退。当时教程里全是 (Y_l^m)，公式里带着 (e^{im\varphi})，展开就是一堆 (a+bi)，看着就头疼。后来把“球谐函数变实值”这件事彻底捋了一遍，才发现原理并不复杂——就是把一组复指数基底重新线性组合成三角函数基底，让每个基函数在实数域里都能直接用。这篇文章把前置知识、构造方法、代码实现和踩坑经验一次讲清楚，适合正在做球谐光照、量子化学基组、球面信号分析，或者被球谐函数困扰的同学参考。

1. 球谐函数为什么天生带复数

1.1 从定义里揪出虚数单位

球谐函数的经典定义长这样：

[
Y_l^m(\theta,\varphi)=N_{lm}P_l^m(\cos\theta)e^{im\varphi}
]

其中 (P_l^m(\cos\theta)) 是连带勒让德函数，(N_{lm}) 是归一化因子，(l) 是阶数（0, 1, 2...），(m) 是磁量子数（(-l \le m \le l)）。这个定义在量子力学里太常用了，因为它是球面拉普拉斯算子的本征函数，也是角动量平方算符和角动量 z 分量算符的共同本征函数。

问题就出在最后那个 (e^{im\varphi}) 上。用欧拉公式展开：

[
e^{im\varphi}=\cos(m\varphi)+i\sin(m\varphi)
]

只要 (m \ne 0)，这个函数在 (\varphi) 方向上就是复数的。这就是球谐函数“天生带复数”的根源。

实际计算中，(m=0) 的项比较幸运，(e^{i0}=1)，函数退化为实数。但 (m \ne 0) 的项真没法绕开虚数。如果你直接拿着 (Y_l^m) 去做数值计算，矩阵是复的，存储要翻倍，求解也慢，很多场景下还得手动取实部或虚部，非常难受。

1.2 复指数项的物理意义

物理学家当初没有刻意回避复数，是因为 (e^{im\varphi}) 在量子力学里有明确的含义：它描述粒子绕 z 轴转动的相位。(m) 的正负对应两种相反的旋转方向，代表角动量 z 分量的大小和方向。使用复值基函数，角动量守恒律在数学上表示得非常简洁。

但工程上，复数不是免费的午餐。我们在做图形学光照、球面信号拟合、形状分析时，往往不关心粒子内部角动量的方向，只需要一组正交函数基去展开定义在球面上的信号。用复数基底意味着每个系数都有实部和虚部，写缓存、写序列化、写渲染 shader 都要多处理一层，容易出错，性能也差。

于是“实值化”（real spherical harmonics）就成了刚需：找到一组实系数正交基，张成和复球谐函数完全相同的函数空间，但每个基函数都是纯实数。

1.3 实值化之后实际上改变了什么

需要说明的是，实值化不是简单地取 (Y_l^m) 的实部或虚部。如果直接取实部，你得到的那组函数虽然也是实函数，但并不同时具备正交归一性，不能直接当基函数用。正确做法是把 (m) 和 (-m) 两个方向的复函数做线性组合，让虚部恰好抵消，同时保证归一化。

换句话说，实值球谐不是对单个复数函数做截断，而是对复函数空间重新选基。这个理解非常重要，后面所有推导、代码、验证都建立在它之上。

2. 实值化原理：把复指数拆成三角函数

2.1 靠欧拉公式构造正余弦组合

核心思路一句话：利用欧拉公式，把两个互为共轭的复指数组合成余弦或正弦。

欧拉公式告诉我们：

[
\cos(m\varphi)=\frac{e^{im\varphi}+e^{-im\varphi}}{2}
]

[
\sin(m\varphi)=\frac{e^{im\varphi}-e^{-im\varphi}}{2i}
]

所以，如果我把 (Y_l^m) 和 (Y_l^{-m}) 按特定比例加起来，虚部就会互相抵消，剩下干净的 (\cos(m\varphi)) 或 (\sin(m\varphi))。

不过这里有一个细节：连带勒让德函数在 (m) 和 (-m) 之间不是简单地相等，而是满足：

[
P_l^{-m}(x)=(-1)^m\frac{(l-m)!}{(l+m)!}P_l^m(x)
]

所以做线性组合时，不仅要考虑 (e^{im\varphi}) 和 (e^{-im\varphi}) 的系数，还要把连带勒让德函数的比例关系放进归一化因子。我在代码实现时习惯直接用绝对值的 (m) 来构造，避免符号推导绕晕。

常用的一种实值化定义是：

• 当 (m>0) 时，实球谐正比于 (P_l^m(\cos\theta)\cos(m\varphi))
• 当 (m<0) 时，实球谐正比于 (P_l^{|m|}(\cos\theta)\sin(|m|\varphi))
• 当 (m=0) 时，直接取 (P_l^0(\cos\theta))

比例系数由正交归一条件确定。经过这样处理后的实球谐，才是真正的一组实正交基。

2.2 前几阶展开：从球函数看坐标

实值化的好处，从低阶项就能直观感受到。以常用的无 Condon-Shortley 相位约定为例，前三阶实球谐可以写成：

[
Y_{00}=\frac{1}{2\sqrt{\pi}}
]

[
Y_{1,-1}=\sqrt{\frac{3}{4\pi}}\frac{y}{r},\quad
Y_{10}=\sqrt{\frac{3}{4\pi}}\frac{z}{r},\quad
Y_{11}=\sqrt{\frac{3}{4\pi}}\frac{x}{r}
]

也就是说，一阶实球谐的三个分量恰好对应于笛卡尔坐标的 (x/r)、(y/r)、(z/r)。这个性质在实际中太有用了——它意味着方向光、法线方向、漫反射响应等几何量可以直接和低阶球谐系数建立对应关系。

二阶的情况更复杂一点，但也同样有直观的几何意义：

[
Y_{2,-2} \propto \frac{xy}{r^2},\quad
Y_{2,-1} \propto \frac{yz}{r^2},\quad
Y_{20} \propto \frac{3z^2-r^2}{r^2},\quad
Y_{21} \propto \frac{xz}{r^2},\quad
Y_{22} \propto \frac{x^2-y^2}{r^2}
]

搞化学的同学看到这一组应该很熟悉，它们就是 d 轨道的五个分量：(d_{xy})、(d_{yz})、(d_{z^2})、(d_{xz})、(d_{x^2-y^2})。这也是为什么量子化学计算里，实球谐比复球谐更常用——分子轨道的空间形状可以直接用实函数画出来，不需要在复数空间里看图。

2.3 正交归一性和完备性

实值化之后，这组实函数依然满足正交归一性：

[
\int_0^{2\pi}\int_0^\pi Y_{lm}(\theta,\varphi)Y_{l'm'}(\theta,\varphi)\sin\theta,d\theta,d\varphi=\delta_{ll'}\delta_{mm'}
]

这一点是整个实球谐应用的基石。有了正交归一性，球面上的任意平方可积函数 (f(\theta,\varphi)) 都可以用实球谐展开：

[
f(\theta,\varphi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}c_{lm}Y_{lm}(\theta,\varphi)
]

系数通过投影得到：

[
c_{lm}=\int_0^{2\pi}\int_0^\pi f(\theta,\varphi)Y_{lm}(\theta,\varphi)\sin\theta,d\theta,d\varphi
]

实际工程中取一个有限的 (l_{\max})，比如 (l_{\max}=2) 或 (3)，就能得到一个低维近似。图形学里的漫反射球谐光照，用的就是 (l_{\max}=2) 的 9 个系数。

3. Python 实操：写一套可用的实值球谐函数

3.1 准备工作与约定

代码层面，最省事的方式是使用 SciPy 的连带勒让德函数 scipy.special.lpmv。它支持非整数阶、整数阶参数，能直接算 (P_l^m(\cos\theta))。

但这里有个大坑，必须提前说明：SciPy 的 lpmv 默认包含 Condon-Shortley 相位，也就是连带勒让德函数定义里带着 ((-1)^m)。如果你直接拿它去构造“图形学约定”的实球谐，一阶项会差一个负号。我在写代码时比较喜欢先去掉这个相位：

python复制import numpy as np
from scipy.special import lpmv, gammaln

def real_sh(l, m, theta, phi, remove_cs=True):
    """
    计算实值球谐函数。

    参数:
        l: 阶数 (0, 1, 2, ...)
        m: 磁量子数 (-l <= m <= l)
        theta: 极角，范围 [0, pi]，从 z 轴正方向算起
        phi:   方位角，范围 [0, 2*pi]，从 x 轴正方向算起
        remove_cs: 是否去掉 Condon-Shortley 相位，默认 True

    返回:
        与 theta, phi 同形状的实数值数组
    """
    m_abs = abs(m)
    Plm = lpmv(m_abs, l, np.cos(theta))

    if remove_cs:
        Plm = Plm * ((-1) ** m_abs)

    # 用 gammaln 计算归一化因子，避免大 l 时阶乘溢出
    ln_norm = 0.5 * (
        np.log(2 * l + 1)
        - np.log(4 * np.pi)
        + gammaln(l - m_abs + 1)
        - gammaln(l + m_abs + 1)
    )
    norm = np.exp(ln_norm)

    if m_abs == 0:
        return norm * Plm
    elif m > 0:
        return np.sqrt(2) * norm * Plm * np.cos(m_abs * phi)
    else:
        return np.sqrt(2) * norm * Plm * np.sin(m_abs * phi)

用这个函数算一阶项，很快能看到结果：

python复制theta = np.array([0.3, 0.8, 1.2])
phi = np.array([0.5, 1.0, 2.0])

print(real_sh(1, 1, theta, phi))   # 应该正比于 sin(theta) * cos(phi)，即 x/r
print(real_sh(1, -1, theta, phi))  # 应该正比于 sin(theta) * sin(phi)，即 y/r
print(real_sh(1, 0, theta, phi))   # 应该正比于 cos(theta)，即 z/r

如果装了 pyshtools，它提供的 real_sph_harm 和这版代码应该高度一致，属于同一个符号约定。

3.2 从复值球谐线性组合的实现

如果你想验证“实值化就是复指数的线性组合”，也可以直接写一个基于复球谐的版本。SciPy 新版本里 scipy.special.sph_harm_y 可以算复球谐，老版本叫 scipy.special.sph_harm，参数顺序是 sph_harm(m, l, theta, phi)。我习惯自己用 lpmv 拼一个，因为依赖更少：

python复制def complex_sh(l, m, theta, phi):
    """复数球谐函数，采用带 Condon-Shortley 相位的定义"""
    m_abs = abs(m)
    Plm = lpmv(m_abs, l, np.cos(theta))

    if m < 0:
        # P_l^{-m} = (-1)^m * (l-m)! / (l+m)! * P_l^m
        from scipy.special import factorial
        Plm = ((-1) ** m_abs) * factorial(l - m_abs) / factorial(l + m_abs) * Plm

    ln_norm = 0.5 * (
        np.log(2 * l + 1)
        - np.log(4 * np.pi)
        + gammaln(l - m + 1)
        - gammaln(l + m + 1)
    )
    return np.exp(ln_norm) * Plm * np.exp(1j * m * phi)

有了复球谐，就能按照线性组合公式构造实球谐。但我必须提醒你：不同文献对“实值化线性组合”的符号定义并不统一，有的地方 (Y_{1,1}) 算出来是正 (x/r)，换个公式就变成负 (x/r)。这不是算法错误，是约定不同。跨库、跨论文对照时，一定要先做符号测试，不要盲目相信公式。

你完全可以先跑通一个最简单的验证：在球面上均匀采样足够多点，分别用 real_sh 算 (Y_{11})，然后乘上 (\sin\theta\cos\varphi) 做积分，看看符号是正还是负，和你的应用场景是否需要对齐。

3.3 归一化和正交性自检

写完代码别急着用，先做两个自检：归一化检查、正交性检查。

归一化检查的思路是对整个球面积分 (Y_{lm}^2)，结果应该是 1：

python复制def check_normalization(l, m, n_theta=400, n_phi=400):
    theta = np.linspace(0, np.pi, n_theta)
    phi = np.linspace(0, 2 * np.pi, n_phi)
    Theta, Phi = np.meshgrid(theta, phi, indexing='ij')
    f = real_sh(l, m, Theta, Phi)
    # 球面积分：f^2 * sin(theta) 在 theta-phi 网格上的二重积分
    integral = np.sum(f**2 * np.sin(Theta))
    integral *= (np.pi / (n_theta - 1)) * (2 * np.pi / (n_phi - 1))
    return integral

print(check_normalization(2, 1))  # 应该接近 1.0

正交性检查类似，换成两个不同 ((l,m)) 组合的乘积积分，结果应该是 0：

python复制def check_orthogonality(l1, m1, l2, m2, n_theta=400, n_phi=400):
    theta = np.linspace(0, np.pi, n_theta)
    phi = np.linspace(0, 2 * np.pi, n_phi)
    Theta, Phi = np.meshgrid(theta, phi, indexing='ij')
    f1 = real_sh(l1, m1, Theta, Phi)
    f2 = real_sh(l2, m2, Theta, Phi)
    integral = np.sum(f1 * f2 * np.sin(Theta))
    integral *= (np.pi / (n_theta - 1)) * (2 * np.pi / (n_phi - 1))
    return integral

print(check_orthogonality(2, 1, 1, -1))  # 应该接近 0.0

这几行代码是我每换一个球谐库必跑的测试。现实情况是：很多库的符号约定、归一化方式存在细微差异，数据移植时肉眼不容易发现，一跑积分立刻现形。

4. 落地应用：实值球谐能解决哪些实际问题

4.1 图形学里的球谐光照

实值球谐在图形学里最出名的应用是球谐光照。把环境光 (L(\theta,\varphi)) 编码成一组 SH 系数，再和余弦加权 BRDF 做卷积，漫反射光照只需要前几阶系数就能算得很准。(l_{\max}=2) 时是 9 个系数，(l_{\max}=3) 时是 16 个系数，采样和存储成本都很低。

为什么图形学非要用实值球谐而不是复值球谐？因为渲染管线里最终要输出的是 RGB 颜色，是实数值。用实基函数，系数直接存储为 float，不用处理复数；在 shader 里重建光照时，直接乘实系数再求和，不需要对虚部做额外处理。

还有一个很实际的原因：可视化。实球谐低阶项的几何意义清楚，你可以直接调试某个方向上的值是否符合预期。复球谐则很难直观检查，出了问题不好定位。

4.2 量子化学与分子轨道

在量子化学计算里，原子轨道基函数常用高斯型函数乘实球谐函数来构造。学过结构化学的同学应该记得 d 轨道的五种形状：(d_{xy})、(d_{yz})、(d_{xz})、(d_{z^2})、(d_{x^2-y^2})，这五个轨道本质上是 (l=2) 实球谐函数的五种方位角分布。

选择实球谐而不是复球谐，最大原因是计算能标量分子轨道时，实数矩阵对角化比复数矩阵快得多，存储也少。另外，实函数可以直接可视化，直接画等值面，不需要做复数到实数的映射。化学软件里几乎清一色使用实球谐作为角向基函数。

4.3 球面信号分析与形状描述

除了图形学和化学，实球谐在生物医学影像、地球物理、天文学里也常见。比如脑皮层表面形状分析，把每个顶点位置或某种标量属性定义为球面函数，用实球谐展开得到一组形状描述子。这类描述子天然具有旋转不变性质（前提是取模或做对齐），便于做分类和检索。

我自己做过的球面散射测量数据拟合也踩过这个坑：数据在球面上分布不均匀，直接做傅里叶变换不行，用实球谐展开后，低阶系数就能反映整体趋势，高阶系数对应局部细节。实球谐的 (Y_{1,-1},Y_{10},Y_{11}) 直接对应三个笛卡尔分量，拟合得到的系数还能解释为平均方向，省了不少事。

4.4 密度估计与预计算辐射传输

PRT（Precomputed Radiance Transfer）这类技术在游戏引擎里出现过很多轮，核心思路是把光照传输函数也投影到球谐基上。由于光照传输矩阵在球谐基下往往比较稀疏，用实球谐可以把矩阵元素直接存成浮点数，引擎运行时做向量-矩阵乘法就行了。

这类应用的共同特征是：信号是实值、目标是加速计算、需要低阶近似。只要看到满足这三个条件，实球谐基本就是首选。

5. 踩坑记录：符号约定、数值稳定性和旋转问题

5.1 符号约定差异对比

我见过太多人在球谐系数上栽跟头，最后查出来只是符号约定不同。这里把常见的约定差异列一下：

这不是谁对谁错的问题。你只需要记住：用别人的代码前，先拿低阶项和坐标做对比，确认符号是否和你的目标一致。

5.2 高 l 与极点的数值坑

实球谐在高阶时（(l>50)）会遇到数值稳定性问题。连带勒让德函数 (P_l^m) 在阶数高时数值范围很大，归一化因子如果直接算阶乘很容易溢出。我的建议是全程用对数形式计算归一化因子，也就是前面代码里 gammaln 的写法，避免 factorial(l - m) 和 factorial(l + m) 直接相除。

另一个容易忽略的坑是极点和球面上的“经度奇异”。当 (\theta=0) 或 (\theta=\pi) 时，(\sin\theta=0)，此时 (\varphi) 没有意义，而球谐函数里包含 (\cos(m\varphi)) 和 (\sin(m\varphi))，在极点上所有不同 (m) 的项必须相等或为零。数值上，如果你在极点直接采样，通常会算出 NaN 或者因浮点误差产生随机数。处理办法是：在球面积分时避开极点，或者把极点上的函数值单独设置为 (m=0) 项的值。

5.3 旋转操作时复值基反而更顺手

实球谐在静态基底下很好用，但一旦涉及旋转，问题就来了。实球谐的旋转矩阵既不是稀疏的，也不是容易构造的。如果你需要在旋转后的球面上做函数展开，或者旋转一组 SH 系数，直接对实基做旋转非常麻烦。

我的习惯是：先把实系数转成复球谐系数，在复基底下用 Wigner D 矩阵完成旋转，再转回实基。虽然多了一次转换，但在数学上更干净，也不容易出错。如果只做小角度旋转或者只旋转少数几个方向，也可以直接用球谐函数的重采样方法，在旋转后的方向上重新求值再投影。

5.4 采样密度不足导致伪系数

球谐展开的系数是用球面积分算出来的。如果采样点太少，尤其是高纬地区采样不够，就会出现混叠，低阶系数被高阶成分污染。这个现象和信号处理里的采样定理是同一回事，只不过这里的“频域”变成了球谐阶数 (l)。

实操建议：做环境贴图编码时，至少要用 32x32 以上的经纬均匀采样；如果信号本身有尖锐的高亮区域，先做低通滤波或者用更高 (l_{\max}) 展开，再截断到目标阶数，否则低阶系数质量会很差。

我个人在实际操作中的体会是，实值球谐最关键的并不是那几个公式，而是一套可以自检的工程流程。拿到一个新库，先用低阶项对照坐标方向，再跑归一化积分和正交性积分，最后才放心使用。符号约定不同不用慌，统一约定后所有结果是一致的。如果你正在做光照、分子轨道、球面数据分析这类事情，建议把上面这组代码和验证函数直接留在项目里，后面换环境、换库、换论文公式时都用得上。