球谐展开与谱分析


SHTOOLS(Tools for working with spherical harmonics) 4.2 是一个专门用来处理球谐函数的工具包,该工具包有 Fortran 版本和 Python 版本,本文主要针对 Python 版本的使用来进行说明。由于使用 python 2.x 版本会出现一些错误,例如使用 read_icgem_gfc 函数来读取 ICGEM_GFC 发布的静态地球重力场位系数文件(icgem2.0 版本)会出现错误,因此推荐使用 python 3.x 版本的 SHTOOLS。

表面球谐函数

通常情况下我们用到的表面球谐(Surface Spherical Harmonics)函数基本为实值函数,因此在没有特殊说明的情况下,默认表面球谐函数为实函数。在介绍表面球谐函数之前,我们先引入平方可积函数的基本概念。

平方可积函数

若一个球面函数 $f(\theta,\phi)$ 为平方可积(square-integrable)函数,那么该函数满足 $$\int_\Omega f^2(\theta,\phi) \, d\Omega < \infty$$,其中 $\Omega$ 代表整个球面域,即 $$d \Omega = \sin \theta d \theta d \phi $$。

表面球谐函数的基本概念

球面上任意一个平方可积函数 $f(\theta,\phi)$ 都可以用一系列正交的表面球谐函数来表示,即 $$f(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^l F_{lm}Y_{lm}(\theta,\phi)$$,其中 $F_{lm}$ 为球谐系数,$Y_{lm}(\theta,\phi)$ 为表面球谐函数,$\theta$ 为余纬,$\phi$ 为经度,$l$ 和 $m$ 分别为球谐阶数和次数。为了方便计算,大多数情况下我们会用到正规化的球谐系数 $f_{lm}$ 和正规化(正规化为 $4\pi$,又称为 geodesy 正规化或者完全正规化)的表面球谐函数 $\mathcal{Y_{lm}}(\theta,\phi)$,使其满足 $$f(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^l f_{lm} \mathcal{Y_{lm}}(\theta,\phi)$$。正规化的表面球谐函数的表达式为 $$\mathcal{Y_{lm}}(\theta,\phi) = \left \lbrace \begin{array}{ll} \mathcal{P_{lm}}(\cos \theta) \cos m \phi & \mbox{if $m \ge 0$} \\ \mathcal{P_{l|m|}}(\cos \theta) \sin |m| \phi & \mbox{if $m < 0$}, \end{array} \right. $$,其中 $\mathcal{P_{lm}}(\mu)$ 为正规化缔合勒让德(normalized associated Legendre)函数,其表达式为 $$\mathcal{P_{lm}}(\mu) = \sqrt{\left(2-\delta_{m0}\right) \left(2l+1\right)\frac{(l-m)!}{(l+m)!}}\, P_{lm}(\mu)$$,其中 $\delta_{m0}$ 为 Kronecker Delta 函数,$P_{lm}(\mu)$ 为缔合勒让德(associated Legendre)函数,其表达式为 $$P_{lm}(\mu) = (-1)^m(1-\mu^2)^{m/2} \frac{d^m}{d\mu^m} P_l(\mu)$$,其中 $P_l(\mu)$ 为勒让德函数,其表达式为 $$P_l(\mu) = \frac{1}{2^l l!}\frac{d^l}{d\mu^l}(\mu^2-1)^l$$。

表面球谐函数的正交性

任意阶次的表面球谐函数满足 $$\int_\Omega \mathcal{Y_{lm}}(\theta,\phi) \mathcal{Y_{l’m’}}(\theta,\phi) d\Omega = 4 \pi \delta_{ll’} \delta_{mm’}$$,这说明表面球谐函数之间满足正交性。利用表面球谐函数的正交性,可得球谐系数 $$f_{lm} = \frac{1}{4\pi} \int_\Omega f(\theta,\phi) \, \mathcal{Y_{lm}}(\theta,\phi) \, d\Omega$$。

功率谱

Parseval 定理:球面函数 $f(\theta,\phi)$ 的平方在整个球面域的积分等于球面的面积乘以球谐系数的平方和。根据 Parseval 定理,我们得到 $$\frac{1}{4\pi} \int_\Omega f^2(\theta,\phi) \, d\Omega
= \sum_{l=0}^{\infty} S_{ff}(l)$$,其中 $\sum_{l=0}^{\infty} S_{ff}(l)$ 为总功率,$S_{ff}(l)$ 为 $l$ 阶的功率,其表达式为 $$S_{ff}(l) = \sum\limits_{m=-l}^l f^2_{lm}$$。由于 $S_{ff}(l)$ 为阶数 $l$ 的函数,因此$S_{ff}(l)$ 即为功率谱(Power spectrum)。与上面类似,根据 Parseval 定理,球面函数 $f(\theta,\phi)$ 和 $g(\theta,\phi)$ 满足 $$\frac{1}{4\pi} \int_\Omega f(\theta,\phi) \, g(\theta,\phi)\, d\Omega
= \sum_{l=0}^{\infty} S_{fg}(l)$$,其中 $\sum_{l=0}^{\infty} S_{fg}(l)$ 为总交叉功率,$S_{fg}(l)$ 为 $l$ 阶的交叉功率,其表达式为 $$S_{fg}(l) = \sum\limits_{m=-l}^l f_{lm}g_{lm}$$。由于 $S_{fg}(l)$ 为阶数 $l$ 的函数,因此 $S_{fg}(l)$ 即为交叉功率谱( Cross Power spectrum)。功率谱密度(Power Spectrum Density)的定义为 $$\mbox{power per $lm$} = \frac{S(l)}{(2l+1)}$$,其中 $S(l)$ 代表 $S_{ff}(l)$ 和 $S_{fg}(l)$。

Admittance 和 Correlation

Admittance 和 Correlation 的定义和物理意义分别如下:

  • Admittance:$Z(l) = \dfrac{S_{fg}(l)}{S_{gg}(l)}$,其中 $S_{fg}(l)= \sum \limits_{m = -\ell}^{\ell}f_{lm}g_{lm}$,$S_{gg}(l) = \sum \limits_{m = -\ell}^{\ell}g_{lm}^2$。Admittance 其实表示向量 $\mathbf A$ 在另一个向量 $\mathbf B$ 的投影与向量 $\mathbf B$ 的比值。可以把 $f_{lm}$ 和 $g_{lm}$ 看成是 $2l+1$ 维的向量,那么如果 $f_{lm}=\alpha g_{lm}$ 则 $Z(l) = \alpha$。
  • Correlation:$\gamma(l) = \frac{S_{fg}(l)}{\sqrt{S_{ff}(l)S_{gg}(l)}}$,其含义很明确,指的的是两个函数在 $l$ 阶的相关系数。

窗口球谐函数

有时我们并不知道全球的数据,而是只知道某一个区域中的数据,例如只知道全球的陆地起伏数据(假设海底的起伏数据未知)。如果对该区域中的数据进行球谐展开,那么需要用窗口球谐函数(本质上是 Slepian 基函数)将谱功率集中在该区域。根据该区域的大小和形状,若窗口的最大阶数为 $L$,那么我们可得到 $(L+1)^2$ 个正交 Slepian 基函数(总功率从大到小排列),记为 $h_{\alpha}(\theta,\phi)$,$\alpha$ 为 $1…(L+1)^2$ 的整数。若该区域中的数据可用 $s(\theta,\phi)$ 来表示,那么将数据投影到各个基函数可得 $$s(\theta,\phi) = \sum_{\alpha=1}^{(L+1)^2} s_{\alpha}h_{\alpha}(\theta,\phi)$$,而 $s(\theta,\phi)$ 又可以表示为 $$s(\theta,\phi) = \sum_{l=0}^{L}\sum_{m=-l}^{l}s_{lm}\mathcal{Y_{lm}}(\theta,\phi)$$。由于 Slepian 基函数的正交性,我们有 $$\int_{\Omega}h_{\alpha}(\theta,\phi) h_{\beta}(\theta,\phi)d\Omega = 4\pi\delta_{\alpha\beta}$$,其中 $$h_{\alpha}(\theta,\phi) = \sum_{l=0}^{L}\sum_{m=-l}^{l}h_{\alpha,lm}\mathcal{Y_{lm}}(\theta,\phi)$$,$$\delta_{\alpha\beta} = \left \lbrace \begin{array}{ll} \sum_{l=0}^{L}\sum_{m=-l}^{l} h_{\alpha,lm} h_{\beta,lm} & \mbox{if $\alpha= \beta$}\\ 0 & \mbox{if $\alpha \neq \beta$} \end{array} \right.$$。根据 Slepian 基函数的正交性,我们最终可以求得 $$s_{\alpha} = \dfrac{\sum_{l=0}^{L}\sum_{m=-l}^{l}h_{\alpha,lm}s_{lm}}{\sum_{l=0}^{L}\sum_{m=-l}^{l}h_{\alpha,lm}^2}$$,考虑到 $$\sum_{l=0}^{L}\sum_{m=-l}^{l}h_{\alpha,lm}s_{lm} = \sum_{l=0}^{L}\sum_{m=-l}^{l}h_{\alpha,lm}f_{lm}\;,$$ 因此 $$s_{\alpha} = \dfrac{\sum_{l=0}^{L}\sum_{m=-l}^{l}h_{\alpha,lm}f_{lm}}{\sum_{l=0}^{L}\sum_{m=-l}^{l}h_{\alpha,lm}^2}\;.$$ 在完成窗口球谐展开后,没有必要用所有 $(L+1)^2$ 个正交 Slepian 基函数来进行球谐恢复,那么用多少个比较合适?Shannon 数给出了一个参考,其定义为 $$N = (L+1)^2\dfrac{A}{4\pi}$$,其中 A 为区域的面积。所以根据投影系数和 Shannon 数,球谐恢复近似为 $$s(\theta,\phi) \approx \sum_{\alpha=1}^{N} s_{\alpha}h_{\alpha}(\theta,\phi)$$。

用球谐函数分析全球地形起伏数据

可视化全球地形起伏数据

下面的代码将全球地形起伏数据 Etopo5 ($5^{’}$x$5^{’}$)可视化。该数据共有 2161×4320 个网格点,其中纬度范围为 $90^{\circ}S$ 到 $90^{\circ}N$;经度范围为 $0^{\circ}E$ 到 $360^{\circ}E-5^{’}$。

输出结果为

全球地形起伏

使用 imshow 作图时要注意原点的设置,默认原点设置在左上角(upper)。因为 etopo5.dat 的第一行数据为南极点,最后一行数据为北极点,因此需要将原点设置在左下角(lower)。 为了方便,在作图之前,需要将 $360^{\circ}$ 对应的列数据(与 $0^{\circ}$ 度对应的列数据相等)补充到网格数据的后面。extent 范围参数是根据二维数组的 index 决定的,例如本例中纬度的 index 范围为 $0\sim2160$,对应 $-90^{\circ}\sim 90^{\circ}$;经度的 index 范围为 $0\sim4320$,对应 $0^{\circ}\sim360^{\circ}$。

将全球地形起伏数据展开到球谐函数

将地形起伏数据展开成球谐函数(Spherical Harmonic)之前,需要将地形起伏数据按行进行反转(使数据的第一行对应 $90^{\circ}N$ 并格式化为 nxn 或 nx2n 的形式,其中球谐函数的最高阶数默认为 $l_{max} = n/2-1$,n 为网格沿纬度方向的维度(必须为偶数)。在此处 n 为 2161,因此需要将最后一行数据去掉,使 n 为 2160,由此 $l_{max}$ 为 1079。

输出结果为

从球谐系数来看 $C_{0,0} = -2.39001111e+03$,这说明地球的平均地形起伏为 $-2390m$。

功率谱(Power Spectrum)计算

输出结果为

power spectrum per l and per lm

对全球网格数据滤波

输出结果为

滤波后的全球地形起伏

用球谐函数分析月球地形模型

月球地形球谐模型的最高阶次为 2600,根据 $l_{max} = n/2-1$ 可知沿经度方向网格点的维度为 5202。该球谐模型是通过 LOLA(Lunar Orbiter Laser Altimeter) 获得测高数据并将其转换为月球主轴坐标系下网格化的地形数据(分辨率为 $1/64^{\circ}$),然后对地形数据球谐展开而得到。

输出结果为

从球谐系数来看 $C_{0,0} = 1.73715120e+06$,这正好对应月球的平均半径。

月球地形的功率谱
999 阶的月球地形

球谐函数在不同坐标系下的转换

以坐标系绕 z 轴(地球自转轴)的旋转为例,简单介绍球谐系数转换的基本原理。最初球面上的网格数据可以用球谐函数 $$f(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^l F_{lm}Y_{lm}(\theta,\phi)$$ 表示。若 x 轴绕 z 轴逆时针旋转 $\alpha$ 得到新坐标系,那么新坐标系下球面上的球谐函数为 $$g(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^l F_{lm}Y_{lm}(\theta,\phi+\alpha) = \sum_{l=0}^{\infty} \sum_{m=-l}^l F_{lm}^{\prime}Y_{lm}(\theta,\phi)$$。通过展开表面球谐函数我们得到以下关系
$$\left [
\begin{matrix}
\mathcal{C_{lm}^{\prime}}\\
\mathcal{S_{lm}^{\prime}}
\end{matrix}
\right ]=\left [
\begin{matrix}
\cos m \alpha & \sin m \alpha \\
-\sin m \alpha & \cos m \alpha
\end{matrix}
\right ]
\left [
\begin{matrix}
\mathcal{C_{lm}}\\
\mathcal{S_{lm}}
\end{matrix}
\right ]
$$。旋转一个坐标系的方式有多种,例如用欧拉角来实现(见下图).

欧拉角

用欧拉角来实现坐标系的旋转可根据旋转轴的顺序分为六种方式,每种方式包含 intrinsic rotations 和 extrinsic rotations。

  1. z-x’-z″ (intrinsic rotations) or z-x-z (extrinsic rotations)
  2. x-y’-x″ (intrinsic rotations) or x-y-x (extrinsic rotations)
  3. y-z’-y″ (intrinsic rotations) or y-z-y (extrinsic rotations)
  4. z-y’-z″ (intrinsic rotations) or z-y-z (extrinsic rotations)
  5. x-z’-x″ (intrinsic rotations) or x-z-x (extrinsic rotations)
  6. y-x’-y″ (intrinsic rotations) or y-x-y (extrinsic rotations)

SHTOOLs 默认采用第四种 intrinsic 旋转方式。以月球的地形模型为例,我们保持月球不动,首先将原坐标系 CS 沿 z 轴逆时针旋转 $30^{\circ}$ 得到一个新的坐标系 CS1,然后将坐标系 CS1 沿 y 轴逆时针旋转 $40^{\circ}$ 得到坐标系 CS2,最后将坐标系 CS2 沿 z 轴逆时针旋转 $50^{\circ}$ 得到最终的坐标系 CS3。

输出结果为

球谐函数在不同坐标系下的转换

功率谱在坐标变换下的不变性

输出结果为

功率谱在坐标变换下的不变性

对网格类的一些操作

输出结果为

网格类的操作
网格类的操作的功率谱

用球谐函数分析局部数据

多窗口谱分析(Multitaper Spectral Analysis)主要包括两类,第一类是有限频谱(bandlimited)分析,第二类是有限空间(spacelimited)分析。

  • 有限频谱分析:bandlimited 指的是在频谱域上找到一系列球谐系数(假设 L 阶,则共有 $(L+1)^2$ 个系数)使其构成一个窗口球谐函数 $h(\theta,\phi)$,使得集中于窗口的能量占全部球域能量的比值 $\lambda$ 最大,即 $$\lambda = \dfrac{\int_R h^2(\theta,\phi)d\Omega}{\int_{\Omega} h^2(\theta,\phi)d\Omega}$$,其中 $R$ 代表窗口区域。  
  • 有限空间分析:spacelimited 指的是在频谱域上找到一系列球谐系数(假设 L 阶,则共有 $(L+1)^2$ 个系数)使其构成一个截断的窗口阶梯函数(真实的窗口阶梯函数在窗口外的值域为 0,并且可球谐展开为无穷阶;但实际应用中无穷阶可用大于 $L$ 的 $L_h$ 近似替代),使得不超过阶数 L 的功率谱占总功率谱的比值 $\lambda$ 最大,即 $$\lambda = \dfrac{\sum_{l=0}^{L}\sum_{m=-l}^{l}h^2_{lm}}{\sum_{l=0}^{\infty}\sum_{m=-l}^{l}h^2_{lm}}$$。

下面分别用有限频谱(bandlimited)分析法来对球盖内的数据进行球谐展开并分析和用有限空间(spacelimited)分析法对陆地起伏数据进行球谐展开并分析。注意:下面将每一个窗口球谐函数都记为一个 taper。

球盖窗口

球盖(Spherical Cap)窗口可通过内建的方法来构建。由于球盖内的数据和球盖外的数据都是已知的,因此我们采用有限频谱分析法来对球盖内的数据进行球谐展开并分析。对于角半径大小为 $\theta_0$ 的球盖,选择 L 大小的方法有两种:

  • 如果只需要一个 taper,那么可通过给定 taper 的集中因子(concentration factor) $\lambda$ 来确定 L 的大小,例如选择集中因子为 0.99;
  • 如果要选取 K 个集中度较好的 taper,那么可通过公式 $L = \sqrt{\frac{2K}{1-\cos\theta_0}}-1$ 来确定 L 的大小。

构建一个角半径为 $30^{\circ}$,L 为 60 的球盖。

输出结果为

结果显示该球盖窗口共有 3721 个 taper。通过上面的公式可知至少需要 249 个集中度较好的 taper 才能较好地进行球谐恢复。对前 20 个 taper 作图

输出结果为

球盖的前 20 个窗口

计算集中度大于 0.99 的 taper 个数和每个 taper 的集中度

输出结果为 179。对前 20 个 taper 的功率谱作图

输出结果为

球盖前 20 个 taper 的功率谱

旋转球盖到指定位置

将球盖旋转到 $45^{\circ}N, 100^{\circ}E$,并且仅旋转前 249 个 taper

输出结果为

旋转后的球盖前 20 个窗口

打印旋转后窗口的基本信息

输出结果为

球谐展开位于球盖中的全球地形起伏数据

将窗口内的起伏数据投影到多个 taper 中,得到每一个 taper 的投影系数。

输出结果为

球谐展开球盖中的数据后又将球谐系数展开到网格

从上图中可以看到球盖内部的球谐恢复很好,而球盖边缘的球谐恢复发生了畸变,因此在使用球盖局部球谐分析时,要使用一个较大的球盖来覆盖感兴趣的区域,以避免区域的边缘发生畸变。

非球盖窗口

通常情况下我们感兴趣的区域并非是一个球盖,而是一个形状不规则的球面多边形。对于这种情况,我们需要用非球盖窗口来进行球谐展开和分析。非球盖窗口可通过遮盖位于某一特定区域的数据来实现,例如通过遮盖海洋区域将其作为一个非球盖窗口。假如我们只知道窗口内的数据,想要对窗口内的数据进行球谐展开,那么使用有限空间分析更加合适。对于窗口面积为 A 的非球盖窗口,选择 L 大小的方法有两种:

  • 如果只需要一个 taper,那么可通过给定 taper 的集中因子(concentration factor) $\lambda$ 来确定 L 的大小,例如选择集中因子为 0.99;
  • 如果要选取 K 个集中度较好的 taper,那么可通过公式 $L = \sqrt{\frac{4\pi K}{A}}-1$ 来确定 L 的大小。

全球陆地面积约占全球面积的 $29\%$,若取 L 为 60,那么通过上面的公式可知至少需要 1079 个集中度较好的 taper。下面根据全球地形起伏来划分窗口:若一个点的地形起伏为正,则认为该点位于陆地,反之则位于海洋。

输出结果为

设置海洋窗口和陆地窗口的阶数,打印窗口的基本信息并分别对海洋窗口和陆地窗口作图

输出结果为

海洋起伏的前 20 个窗口
陆地起伏的前 20 个窗口

对海洋窗口和陆地窗口的功率谱作图

输出结果为

海洋起伏前 20 个窗口的功率谱
陆地起伏前 20 个窗口的功率谱

通过遮盖海洋起伏来构建陆地起伏的窗口

输出结果为

为了数值上的方便,我们将陆地窗口以外的地形起伏设置为 0

输出结果为

将窗口内的陆地起伏数据进行球谐展开并截断到 60 阶,然后再恢复到网格

输出结果为

遮盖海洋起伏后的陆地起伏

将窗口内的起伏数据投影到多个陆地 taper 中,得到每一个陆地 taper 的投影系数

输出结果为

球谐展开陆地数据后又将球谐系数展开到网格

从上图中可以看到陆地窗口内部的球谐恢复很好。与球盖窗口一样,陆地窗口边缘的球谐恢复也会发生畸变,因此在使用非球盖局部球谐分析时,最好使用一个较大的窗口来覆盖感兴趣的区域,以避免区域的边缘发生畸变。一般沿区域边缘向外延拓一段距离来建立窗口,距离的大小要根据感兴趣区域的大小来定,例如对于亚洲区域,可向外延拓 $5^{\circ}$。


参考

关于 “球谐展开与谱分析” 的 1 个意见

发表评论

电子邮件地址不会被公开。 必填项已用*标注