New algorithms for spherical harmonic analysis of area mean values over blocks delineated by equiangular and Gaussian grids

Abstract

Spherical harmonic analysis is widely used in all aspects of geoscience. Exact quadrature methods are available for the spherical harmonic analysis of band-limited point values at the grid points of equiangular and Gaussian grids. However, no similarly exact quadrature methods are available for the spherical harmonic analysis of area mean values over the blocks delineated by these grids. In this study, new algorithms appropriate for the exact spherical harmonic analysis of the band-limited area mean values over the blocks delineated by equiangular and Gaussian grids are proposed. For band-limited data, precision that is between that of the least-squares estimation method and of the approximate quadrature methods can be achieved by using the new algorithms. Regarding the computational complexity, fewer operations are needed by the new methods as compared to those needed by the least-squares estimation method and the approximate quadrature methods in the preparation stage when the maximum degree of the spherical harmonic analysis is very large. Simulation experiments are performed to compare the ability to recover the spherical harmonic coefficients by using the least-squares estimation method, the approximate quadrature methods and these new algorithms from aliased data with aliasing components of realistic magnitudes. The results suggest that these new algorithms, with time complexity one order less than that of the least-squares estimation method in the solving stage, perform roughly the same as the least-squares estimation method in recovering spherical harmonic coefficients from the aliased data.

Appendices

Appendix A: Transformation from ${\mathbf{c}}^{\mathbf{^{\prime}}}$ to $\mathbf{c}$

The relationship between ${\mathbf{c}}^{{^{\prime}}}$ and $\mathbf{c}$ can be denoted as:

$${\mathbf{M}}_{\mathbf{m}+1}{\mathbf{M}}_{\mathbf{m}+2}{\mathbf{M}}_{\mathbf{m}+3}\cdots {\mathbf{M}}_{{{\varvec{l}}}_{\mathbf{m}\mathbf{a}\mathbf{x}}}\mathbf{c}={\mathbf{c}}^{\mathbf{^{\prime}}}$$

(A1)

or written as:

$$\mathbf{c}={\mathbf{M}}_{{{\varvec{l}}}_{\mathbf{m}\mathbf{a}\mathbf{x}}}^{-1}\cdots {\mathbf{M}}_{\mathbf{m}+3}^{-1}{\mathbf{M}}_{\mathbf{m}+2}^{-1}{\mathbf{M}}_{\mathbf{m}+1}^{-1}{\mathbf{c}}^{\mathbf{^{\prime}}}$$

(A2)

${\mathbf{M}}_{\mathbf{m}+\mathbf{k}}^{-1}$ can be written as:

$${\mathbf{M}}_{\mathbf{m}+\mathbf{k}}^{-1}=\left[\begin{array}{llll}{{M}^{^{\prime}}}_{\left(m+k\right)}^{11}& {{M}^{^{\prime}}}_{\left(m+k\right)}^{12}& \cdots & {{M}^{^{\prime}}}_{\left(m+k\right)}^{1,{l}_{{\rm max}}-m+1}\\ {{M}^{^{\prime}}}_{\left(m+k\right)}^{21}& \ddots & & \vdots \\ \vdots & & \ddots & \vdots \\ {{M}^{^{\prime}}}_{\left(m+k\right)}^{{l}_{{\rm max}}-m+\mathrm{1,1}}& \cdots & \cdots & {{M}^{^{\prime}}}_{\left(m+k\right)}^{{l}_{{\rm max}}-m+1,{l}_{{\rm max}}-m+1}\end{array}\right]$$

(A3)

with

$${{M}^{^{\prime}}}_{m+k}^{p,q}=\left\{\begin{array}{l} -\frac{\left(m+k-2\right)}{{a}_{m+k-1,m}} \quad \left(p=k-1;q=k+1;k>1\right) \\ \frac{\left(m+k+1\right)}{{a}_{m+k,m}} \quad \, \left(p=q=k+1\right) \\ 1 \quad \quad \quad \left(p=q\ne k+1\right)\\ 0 \quad \quad \quad \mathrm{otherwise}\end{array}\right.$$

(A4)

We denote the vectors ${\mathbf{d}}^{(\mathrm{k})}$ as:

$${\mathbf{d}}^{\left(\mathrm{k}\right)}={\left[\begin{array}{llll}{d}_{mm}^{\left(k\right)}& {d}_{m+1,m}^{\left(k\right)}& \cdots & {d}_{{l}_{{\rm max}},m}^{\left(k\right)}\end{array}\right]}^{{\varvec{T}}} \left(0\le k\le {l}_{{\rm max}}-m\right)$$

(A5)

$${\mathbf{d}}^{\left(0\right)}={\mathbf{c}}^{\mathbf{^{\prime}}}$$

(A6)

$${\mathbf{d}}^{\left(\mathrm{k}+1\right)}={\mathbf{M}}_{\mathbf{m}+\mathbf{k}}^{-1}{\mathbf{d}}^{\left(\mathrm{k}\right)}$$

(A7)

$${\mathbf{d}}^{\left({l}_{{\rm max}}-m\right)}=\mathbf{c}$$

(A8)

Equation (A7) can be expressed as:

$$\left\{\begin{array}{l}\begin{array}{l}{d}_{m+k-2,m}^{\left(k+1\right)}={d}_{m+k-2,m}^{\left(k\right)}-\frac{\left(m+k-2\right)}{{a}_{m+k-1,m}}{d}_{m+k,m}^{\left(k\right)} \left(k>1\right) \\ {d}_{m+k,m}^{\left(k+1\right)}=\frac{\left(m+k+1\right)}{{a}_{m+k,m}}{d}_{m+k,m}^{\left(k\right)} \left(k>0\right) \end{array}\\ {d}_{l,m}^{\left(k+1\right)}={d}_{l,m}^{\left(k\right)} \left(l\ne m+k;l\ne m+k-2;m\le l<{l}_{{\rm max}}\right)\end{array}\right.$$

(A9)

Note that in the first two equations of Equation (A9), only two elements of ${\mathbf{d}}^{(\mathrm{k}+1)}$ are involved. Hence, the multiplication of ${\mathbf{M}}_{\mathbf{m}+\mathbf{k}}^{-1}$ with ${\mathbf{d}}^{(\mathrm{k})}$ can be replaced with four multiplication operations and one subtraction operation. By repeated use of Equation (A9), one can obtain:

$${\bar{C}}_{m+k,m}=\left\{\begin{array}{l}{\bar{C}^{\prime}}_{m+k,m}-\frac{m}{{a}_{m+1,m}}{{\bar{C}}^{^{\prime}}}_{m+2,m} \left(k=0\right)\\ \frac{\left(m+k+1\right)}{{a}_{m+k,m}}{\bar{C}^{\prime}}_{m+k,m}-\frac{\left(m+k\right)}{{a}_{m+k+1,m}}{{\bar{C}}^{^{\prime}}}_{m+k+2,m} \left(0<k<{l}_{{\rm max}}-m\right)\\ \frac{\left(m+k+1\right)}{{a}_{m+k,m}}{\bar{C}^{\prime}}_{m+k,m} \left(k={l}_{{\rm max}}-m\right)\end{array}\right.$$

(A10)

Appendix B: Estimation of ${u}_{{l}_{1}m}^{i}$

According to Eq. (75), the definition of ${u}_{{l}_{1}m}^{i}$ can be given as:

$$\sum_{i=0}^{N-1}{u}_{{l}_{1}m}^{i}\left({\bar{P}}_{{l}_{2}m}^{i+1}{\mathrm{sin}}^{2}\left({\theta }_{i+1}\right)-{\bar{P}}_{{l}_{2}m}^{i}{\mathrm{sin}}^{2}\left({\theta }_{i}\right)\right)=\sum_{i=0}^{N-1}{v}_{{l}_{1}m}^{i}{\bar{P}}_{{l}_{1}m}^{i}{\bar{P}}_{{l}_{2}m}^{i}{\mathrm{sin}}^{2}\left({\theta }_{i}\right)$$

(B1)

Note that the left side of equation (B1) can be expanded as:

$$\sum_{i=0}^{N-1}{u}_{{l}_{1}m}^{i}\left({\bar{P}}_{{l}_{2}m}^{i+1}{\mathrm{sin}}^{2}\left({\theta }_{i+1}\right)-{\bar{P}}_{{l}_{2}m}^{i}{\mathrm{sin}}^{2}\left({\theta }_{i}\right)\right)={u}_{{l}_{1}m}^{N-1}{\bar{P}}_{{l}_{2}m}^{N}{\mathrm{sin}}^{2}\left({\theta }_{N}\right)+\sum_{i=1}^{N-1}\left({u}_{{l}_{1}m}^{i-1}-{u}_{{l}_{1}m}^{i}\right){\bar{P}}_{{l}_{2}m}^{i}{\mathrm{sin}}^{2}\left({\theta }_{i}\right)-{u}_{{l}_{1}m}^{0}{\bar{P}}_{{l}_{2}m}^{0}{\mathrm{sin}}^{2}\left({\theta }_{0}\right)=\sum_{i=1}^{N-1}\left({u}_{{l}_{1}m}^{i-1}-{u}_{{l}_{1}m}^{i}\right){\bar{P}}_{{l}_{2}m}^{i}{\mathrm{sin}}^{2}\left({\theta }_{i}\right)$$

(B2)

Note that the second equation of (B2) takes advantage of the following equation:

$$\mathrm{sin}{\theta }_{N}=\mathrm{sin}{\theta }_{0}=0$$

(B3)

On the right-hand side of equation (B1) ${v}_{{l}_{1}m}^{0}=0$. We find that the following equations satisfy Equation (B1):

$${u}_{{l}_{1}m}^{i-1}-{u}_{{l}_{1}m}^{i}={v}_{{l}_{1}m}^{i} \left(i=\mathrm{1,2},3\cdots ,N-1\right)$$

(B4)

In Equation (B4), the number of unknowns ${u}_{{l}_{1}m}^{i}\left(i=\mathrm{0,1},\mathrm{2,3}\cdots ,N-1\right)$ is $N$ while there are $\left(N-1\right)$ equations. Consequently, there are an infinite number of solutions for ${u}_{{l}_{1}m}^{i}\left(i=\mathrm{0,1},\mathrm{2,3}\cdots ,N-1\right)$ in Equation (B4). One solution could be:

$$\left\{\begin{array}{l}{u}_{{l}_{1}m}^{0}=0 \\ {u}_{{l}_{1}m}^{i}={u}_{{l}_{1}m}^{i-1}-{v}_{{l}_{1}m}^{i} \left(i=\mathrm{1,2},3\cdots ,N-1\right)\end{array}\right.$$

(B5)

There are other useful solutions as given below. Note that from Equation (B4), we have:

$${u}_{{l}_{1}m}^{i}={u}_{{l}_{1}m}^{0}-\sum_{k=1}^{i}{v}_{{l}_{1}m}^{k}$$

(B6)

$${u}_{{l}_{1}m}^{N-1-i}={u}_{{l}_{1}m}^{N-1}+\sum_{k=1}^{i}{v}_{{l}_{1}m}^{N-k}$$

(B7)

If $\left({l}_{1}-m\right)$ is odd, the following equations hold:

$${\bar{P}}_{{l}_{1}m}^{i}=-{\bar{P}}_{{l}_{1}m}^{N-1-i}$$

(B8)

$${v}_{{l}_{1}m}^{i}=-{v}_{{l}_{1}m}^{N-1-i}$$

(B9)

From Equation (B4), we have:

$${u}_{{l}_{1}m}^{0}-{u}_{{l}_{1}m}^{N-1}=\sum_{k=1}^{N-1}{v}_{{l}_{1}m}^{k}=0 \quad \quad \left({l}_{1}-m\right) \mathrm{is odd}$$

(B10)

If we set ${u}_{{l}_{1}m}^{0}$ to 0, then ${u}_{{l}_{1}m}^{0}={u}_{{l}_{1}m}^{N-1}=0$. Thus, from Equation (B6) and (B7), we have:

$${u}_{{l}_{1}m}^{i}=-\sum_{k=1}^{i}{v}_{{l}_{1}m}^{k}=\sum_{k=1}^{i}{v}_{{l}_{1}m}^{N-k}={u}_{{l}_{1}m}^{N-1-i}$$

(B11)

If $\left({l}_{1}-m\right)$ is even, the following equations hold:

$${\bar{P}}_{{l}_{1}m}^{i}={\bar{P}}_{{l}_{1}m}^{N-1-i}$$

(B12)

$${v}_{{l}_{1}m}^{i}={v}_{{l}_{1}m}^{N-1-i}$$

(B13)

From Equation (B4), we have:

$${u}_{{l}_{1}m}^{0}-{u}_{{l}_{1}m}^{N-1}=\sum_{i=1}^{N-1}{v}_{{l}_{1}m}^{k}\ne 0 \quad \quad \left({l}_{1}-m\right) \mathrm{is even}$$

(B14)

If we set ${u}_{{l}_{1}m}^{0}$ to $1/2\sum_{i=1}^{N-1}{v}_{{l}_{1}m}^{k}$, then ${u}_{{l}_{1}m}^{0}=-{u}_{{l}_{1}m}^{N-1}=1/2\sum_{i=1}^{N-1}{v}_{{l}_{1}m}^{k}$. Thus, from Equations (B6) and (B7), we have:

$${u}_{{l}_{1}m}^{i}=\frac{1}{2}\sum_{i=1}^{N-1}{v}_{{l}_{1}m}^{k}-\sum_{k=1}^{i}{v}_{{l}_{1}m}^{k}=\frac{1}{2}\sum_{i=1}^{N-1}{v}_{{l}_{1}m}^{k}+\sum_{k=1}^{i}{v}_{{l}_{1}m}^{N-k}=-{u}_{{l}_{1}m}^{N-1-i}$$

(B15)

In summary, if $\left({l}_{1}-m\right)$ is odd and ${u}_{{l}_{1}m}^{0}=0$, we have:

$${u}_{{l}_{1}m}^{i}={u}_{{l}_{1}m}^{N-1-i}$$

(B16)

If $\left({l}_{1}-m\right)$ is even and ${u}_{{l}_{1}m}^{0}=1/2\sum_{i=1}^{N-1}{v}_{{l}_{1}m}^{k}$, we have:

$${u}_{{l}_{1}m}^{i}=-{u}_{{l}_{1}m}^{N-1-i}$$

(B17)

Suppose $\mathbf{e}\mathbf{^{\prime}}$ in Eq. (71) can be denoted as:

$${\mathbf{e}}^{\mathbf{^{\prime}}}={\left[\begin{array}{llll}{e}^{0}& {e}^{1}& \cdots & {e}^{N-1}\end{array}\right]}^{{\varvec{T}}}$$

(B18)

Then, the following equation can be used to reduce the number of operations if Equations (B16) and (B17) are used to derive the corresponding ${u}_{{l}_{1}m}^{i}$:

$$\sum_{i=0}^{N-1}{u}_{{l}_{1}m}^{i}{e}^{i}=\left\{\begin{array}{l}\sum_{i=0}^{\frac{N}{2}-1}{u}_{{l}_{1}m}^{i}\left({e}^{i}+{e}^{N-1-i}\right) \quad \quad \left({l}_{1}-m\right)\mathrm{ is odd}\\ \sum_{i=0}^{\frac{N}{2}-1}{u}_{{l}_{1}m}^{i}\left({e}^{i}-{e}^{N-1-i}\right) \quad \quad \left({l}_{1}-m\right)\mathrm{ is even}\end{array}\right.$$

(B19)

Appendix C: Orthogonality and the aliasing-free area of GLQ

The discrete orthogonal equation used in the GLQ method can be written as (Rexer and Hirt 2015; Sneeuw 1994):

$$\sum_{i=0}^{N}{w}_{{l}_{1}m}^{GLQ,i}{\bar{P}}_{{l}_{1}m}^{i}{\bar{P}}_{{l}_{2}m}^{i}={\delta }_{{l}_{1}{l}_{2}}$$

(C1)

where ${\theta }_{i}(i=\mathrm{0,1},\cdots N)$ is the zero value of the Legendre polynomials ${\bar{P}}_{N+\mathrm{1,0}}(\mathrm{cos}\theta )$. Equation (C1) always holds when both ${l}_{1}$ and ${l}_{2}$ are smaller than $N$. However, when one of ${l}_{1}$ and ${l}_{2}$ is equal to or larger than $N$, Equation (C1) sometimes holds. This can be proved as shown below.

The product of two associated Legendre functions can be expanded to (Sneeuw 1994):

$${\bar{P}}_{{l}_{1}m}\left(x\right){\bar{P}}_{{l}_{2}m}\left(x\right)=\sum_{n=0}^{{l}_{1}+{l}_{2}}{a}_{n}{x}^{n}$$

(C2)

The associated Legendre functions are orthogonal as:

$${\int }_{0}^{\pi }{\bar{P}}_{{l}_{1}m}\left(\mathrm{cos}\theta \right){\bar{P}}_{{l}_{2}m}\left(\mathrm{cos}\theta \right)\mathrm{sin}\theta d\theta =2\left(2-{\delta }_{m0}\right){\delta }_{{l}_{1}{l}_{2}}$$

(C3)

According to the characteristics of Gauss–Legendre quadrature (Sneeuw 1994; Ye and Shen 2005), as long as $({l}_{1}+{l}_{2})$ is no larger than $2N+1$, Gauss–Legendre quadrature can be used to calculate the integral on the left-hand side of equation (C2) exactly:

$${\int }_{0}^{\pi }\sum_{n=0}^{{l}_{1}+{l}_{2}}{a}_{n}{x}^{n}dx=\sum_{i=0}^{N}{w}_{{l}_{1}m}^{Q,i}\left(\sum_{n=0}^{{l}_{1}+{l}_{2}}{a}_{n}{x}_{i}^{n}\right)$$

(C4)

where ${w}_{{l}_{1}m}^{Q,i}$ are the known quadrature weights. By combining Equations (C1), (C2) and (C3), we have:

$$2\left(2-{\delta }_{m0}\right){\delta }_{{l}_{1}{l}_{2}}=\sum_{i=0}^{N}{w}_{{l}_{1}m}^{Q,i}{\bar{P}}_{{l}_{1}m}^{i}{\bar{P}}_{{l}_{2}m}^{i}$$

(C5)

with the prerequisite that:

$${l}_{1}+{l}_{2}\le 2N+1$$

(C6)

One can see that by replacing ${w}_{{l}_{1}m}^{Q,i}$ in Equation (C4) with $2\left(2-{\delta }_{m0}\right){w}_{{l}_{1}m}^{GLQ,i}$, we obtain Equation (C1) with the prerequisite given in Equation (C5). Equations (C4) and (C5) are derived here to prove the existence of an ‘aliasing-free area’, as mentioned in Sect. 3.3.