BNU-LQP
0%

球谐函数(Spherical Harmonics)-近代数学中最著名的函数之一

发表于 更新于

前言

球谐函数是近代数学中最著名的函数之一，在光学、天体运动、量子力学中的应用十分广泛。而在地球重力场模型构建过程中，球谐函数及其导数与定积分的计算是必不可少的环节，其精度和稳定性直接影响着重力场模型的构建。

球谐函数的数学理论

球坐标系下的Laplace方程

球谐函数的数学概念和性质在每一本数学物理方法书中都有非常详细的推到和论述。球谐函数是傅里叶级数的高维类比，由一组表示球体表面的基函数构成。同时，球谐函数是Laplace算子角动量在三个纬度上的特征函数。三维直角坐标系下，Laplace方程如下：

将其转化为球坐标系，带入,于是得到球坐标系下的Laplace方程：

球谐函数

方程的理论推导

对于任意一个球坐标系下的函数,我们将其分离变量，可以得到

  • 表示距离部分；
  • ，而这一部分就是球谐函数。

我们将分离后的函数带入上一节的方程，化简后可以得到：

方程的左边是关于R的函数，方程的右边是关于的函数，显然两侧的等式要成立，左侧和右侧的函数需要同时等于一个常数，可以得到两个新方程：

其中第二个方程就称为球函数方程，令,进一步分离变量，得到了一个极为特殊的方程：阶连带勒让德方程：

勒让德方程是物理（引力场、电磁辐射、量子力学）和工程领域里面常常遇到的一类常微分方程。这个方程的解，也就是的表达式为连带勒让德函数，记为。

其中关于连带勒让德多项式，有以下几点需要注意： - 是对求导数。 - 是连带勒让德函数，其定义域是。 - 分别是连带勒让德函数的阶(degree)和次(order), 。 - 广义的勒让德多项式是复数，而连带勒让德则是指实数部分。

以下连带勒让德函数的前2阶表达式：
0 0 1
1 -1
1 0
0 0
2 -2
2 -1
2 0
2 1
2 2

将其带入，即可得到球坐标系下的解，这个解就是球谐函数:

利用上面连带勒让德函数的结论，我们可以推导出前二阶的球谐表达式：
0 0
1 -1
1 0
1 1
2 -2
2 -1
2 0
2 1
2 2

Note:我们一般都不会直接求解，计算复杂度和重复度都很高，在计算高阶球谐函数的时候，更不能Hardcoding,所以一般情况下，都是通过递推迭代的方法得到。

球谐函数的基本性质

球谐函数在球面具有正交性，完备性，以及旋转不变性

  • 正交性：说明基函数两两相互正交。
  • 完备性：说明基函数的线性叠加可以逼近任意一个球体（函数）。
  • 旋转不变性：球谐函数只与有关，与半径无关。故球谐函数是旋转对称的。

球谐（基）函是在复数域的基础上定义的，但是在应用过程中，一般只关心实数域内定义的球函数。所以这篇文章只关心实球面调和（Real Spherical Harmonics）。下面给出了前10阶实球面函数展开后的可视化结果（红色表示正数；蓝色表示负数）。

不同阶次的实球面展开图

球谐函数的逐阶可视化

手动实现（Hardcoding）

根据上述推导出来的公式，我们可以写出对应的勒让德多项式、球谐函数等特殊函数。具体代码如下所示，也可参考这一篇文章氢原子波函数可视化（Python）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#连带勒让德多项式
def LegendrePolynomials(N, x):
    """
    N 为勒让德多项式的阶数
    x 为自变量 sym.symbol对象
    """
    if N == 0:
        return 1
    if N == 1:
        return x
    p0 = LegendrePolynomials(0, x)
    p1 = LegendrePolynomials(1, x)
    assert N >= 2
    for i in range(1, N):
        p = (2 * i + 1) / (i + 1) * x * p1 - i / (i + 1) * p0
        p0 = p1
        p1 = p
    return sym.simplify(p1)

# Factorial function
# 分数阶乘函数
def Factorial(n):
    if n == 0:
        return 1
    elif n < 0:
        print('there is a number < 0 in the Factorial function!')
        return 1
    else:
        s = 1
        while n > 0:
            s *= n
            n -= 1
        return s
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# Associated Laguerre polynomials 
# 输出对应k阶参数p的关联拉盖尔多项式系数列表
def Associate_Laguerre(k, p):
    c_list = []
    for j in range(k+1):
        c_u = Factorial(k+p)
        c_l = Factorial(k-j)*Factorial(p+j)*Factorial(j)
        c = math.pow(-1, j)*c_u/c_l
        c_list.append(c)            
    return c_list

# Testing Laguerre polynomial
for i in range(5):
    print(Laguerre(i))
# Testing Associated Laguerre polynomials 
for n in range(5):
    for l in range(n):
        p = l + 0.5
        k = (n-l)//2
        print(Associate_Laguerre(k, p))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
# 生成勒让德多项式和关联勒让德多项式
import math
import numpy as np

# Comb function
# 排列组合函数
def comb(k,l):
    result_u = 1
    temp = k
    for i in range(l):
        result_u *= temp
        temp -= 1
    result_l = math.factorial(l)
    return result_u/result_l

# Legendre Polynomials  
# 输出对应n阶之下所有阶数勒让德多项式系数列表 
def legendre(n):
    a = np.zeros((n, n))
    for l in range(n):
        for k in range(l+1):
            if l==0 and k==0:
                a[l][k] = 1
            elif (l+k) % 2 == 0:
                a_back = comb(l, int((l+k)/2))
                a_u = np.prod(np.arange(1, l+1)+np.ones(l)*k)
                a[l][k] = np.power(-1, (l-k)/2)*a_u*a_back/(2**l*math.factorial(l))
    return a

# Associated Legendre Polynomials
# 输出对应参数为n，m，x的关联勒让德多项式的值 
def associate_legendre(n, m, x):
    leg = legendre(n+1)[-1]
    result_front = (1-x**2)**(m/2)
    ass_leg = leg[m:]
    temp = 0
    for i in range(len(ass_leg)):
        driv_max_num = i + m
        for j in range(m):
            ass_leg[i] *= driv_max_num
            driv_max_num -= 1
        temp += ass_leg[i]*(x**i)
    return result_front*temp

# Testing Legendre Polynomials
print(legendre(6))
# Testing Associated Legendre Polynomials
print(associate_legendre(2, 1, 0.5))

第三方库

在python的scipy库中也提供了特殊函数如广义勒让德函数、球谐函数等，我们也可以通过调用这些特殊函数来可视化球谐基函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# encoding:utf-8

import matplotlib.pyplot as plt, numpy as np
from scipy import special

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection="3d")
# 
theta, phi = np.linspace(0, np.pi, 100), np.linspace(0, 2*np.pi, 100)
THETA, PHI = np.meshgrid(theta, phi)

## Generate animation
def init():
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.get_cmap('jet'), linewidth=0, antialiased=False, alpha=0.5)
    return fig,

def animate(i):
    ax.view_init(elev=30, azim=0 + 5 * i)
    return fig,

ani = animation.FuncAnimation(fig, animate, init_func=init, frames=36, interval=20, blit=True)

ani.save("test.gif", writer='pillow', fps=30)

通过改变阶数()和次数（）的大小，即可得到不同球谐基函数的可视化结果

球谐基函数实部可视化图

其中，红色代表正值，蓝色代表负值，离中心越远的地方，

球谐函数的线性叠加

理论上，球谐基函数构成了一个函数空间集合（），因此，只要给定了一个三维的球面函数，我们就可以通过无穷多个球谐函数的线性组合，近似的表达出来：

通过上述的方程可知，要将给定球面函数（形状）进行球谐展开，问题的关键在于如何计算球谐函数中的系数。若函数的表达式已知，则可以通过积分计算出精确解，若没有解析式或不可积，则需要进行数值求解。比如Monte-Carlo积分等。我们以大地水准面为例。

大地水准面

地球表面是一个崎岖不平、极不规则的表面，对于地球测量而言，这样的表面无法用数学语言表达出来，不能作为测量和制图的基准面。由于地球表面70%以上几乎都被海洋所覆盖，故假想自由静止的水面在陆地上延伸所形成的连续闭合的曲面称为大地水准面，因此大地水准面同时也是重力等位面，该面上各个点的重力加速度相等。GRACE重力卫星通过GPS和K波段测距仪以1cm的精度精确测量两颗子卫星的间距，从而获得地球重力场的信息，经过滤波处理后可以得到地球静态的重力场模型。再通过球谐函数的线性叠加即可得到大地水准面的形态，其大地水准面表达式如下：

其中：是地球长半径，是归一化连带勒让德多项式；

由于大地水准面接近旋转椭球面，而旋转椭球面又是大地测量的基准面，如果采用实际地球的长半径、扁率，引力场参数及旋转角速度作为椭球参数，就能得到一个与大地水准面几何形状和外部重力场都很符合的水准椭球。大地水准面扣除参考椭球面（一般为WGS84椭球体）就可以推算出大地水准面的起伏和重力异常。

大地水准面3D模型示意图