正常重力位

正常重力位（Normal gravity potential）是大地測量學中用於對地球的真實重力位進行近似的數學工具，是一個規則的、較為簡單的重力位函數。正常重力是正常重力位的梯度。

地球的真實重力位在研究地球形狀及其外部重力場的過程中是待解的未知量，且對其直接計算需要了解地球內部的質量分布，理論上被無法精確求得，且反演過程計算複雜。因此，研究中常選用一個規則的正常橢球體對地球的形狀進行近似，將其產生的重力場中的重力位稱為正常重力位，而把真實重力位與正常重力位的差異稱作擾動位。選擇適當的正常重力位，可以使擾動位成為微小量，便於以線性近似的方式對其進行求解。在地球重力場中，正常重力位可以占到到真實重力位的 $99.999 5\%$。

## 確定要求

正常重力位應當具有以下特性：

- 旋轉對稱性，即與經度無關
- 相對於赤道面對稱
- 在橢球面上是常數

這些特性保證了正常重力位是規則分布的，利用這些性質能夠簡化複雜的計算過程。最後一個條件保證了還正常橢球體的表面是一個重力等位面。

## 數學表達

正常重力位 $U$ 包含兩部分，一部分是因正常橢球體的質量而產生的引力位 $V$，另一部分是因正常橢球體的繞軸自轉而產生的離心力位 $\Phi$：

$$U = V + \Phi = V + {1\over2}\omega^2(x^2+y^2)$$

上式中，$\omega$ 表示地球自轉的角速度，$(x,y,z)$ 是正常橢球體外部空間中某點的笛卡爾坐標，該坐標系的 $Z$ 軸與橢球體的自轉軸平行或重合。其中的離心力位 $\Phi$ 可以由點的坐標直接求得，不必展開為級數，只需對引力位 $V$ 進行展開即可。

### 斯托克斯方法

採用球面作為正常橢球體的近似，將橢球體外部（$r>a$）的引力位展開為邊界面（$r = a$）上的球諧級數，則其表達式為：

$$V = {\sum_{n=0}^{\infin}}\frac{1}{r^{n+1}} {\sum_{m=0}^n}P_{nm}(\cos\theta)\left[A_{nm}\cos m\lambda+B_{nm}\sin m\lambda\right]$$

或

$$V = \frac{GM}{r} {\sum_{n=0}^{\infin}}\left(\frac{a}{r}\right)^n {\sum_{m=0}^n}P_{nm}(\cos\theta)\left[C_{nm}\cos m\lambda+S_{nm}\sin m\lambda\right]$$

上式中各項符號的物理或數學意義如下：

- $(r, \theta, \lambda)$是空間中某特定點的球坐標，$r$ 是空間中某點的地心距離， $\theta$ 和 $\lambda$ 分別是該點的極距和經度
- $GM$ 為正常橢球體的地心引力常數
- $a$ 為橢球的長半徑
- $P_n$ 是 $n\ $階勒讓德多項式，$P_{nm}$ 是 $n\ $階 $m\ $次締合勒讓德多項式
- $A_{nm}$ 和 $B_{nm}$、$C_{nm}$ 和 $S_{nm}$ 是由積分求得的球諧係數，且積分表達式為：

$$\begin{Bmatrix} A_{nm} \\ B_{nm} \end{Bmatrix} = (2-\delta_{m,0}){(n-m)!\over(n+m)!}G \iiint\limits_{\text{earth}} r'^n P_{nm}(\cos\theta') \begin{Bmatrix} \cos m\lambda' \\ \sin m\lambda' \end{Bmatrix}\operatorname{d}\!M' \quad$$

或

$$\begin{Bmatrix} C_{nm} \\ S_{nm} \end{Bmatrix} = (2-\delta_{m,0}){(n-m)!\over(n+m)!}{1 \over M} \iiint\limits_{\text{earth}} {\left({r' \over a}\right)}^n P_{nm}(\cos\theta') \begin{Bmatrix} \cos m\lambda' \\ \sin m\lambda' \end{Bmatrix}\operatorname{d}\!M' \quad$$

其中 $\delta_{m,0}$ 為克羅內克δ函數，當 $m=0$ 時為零，其他情況下為一。

展開後的球諧級數具有以下性質：

- 正常重力場具有旋轉對稱性，即其產生的正常重力位與經度無關，因此上式中面諧函數的部分可被簡化作 $A_{n,0}P_n(\cos\theta)$ 或 $C_{n,0}P_n(\cos\theta)$
- 當坐標系的原點與橢球體的質心重合時，可證明該級數中的一階係數為零，即 $A_{1,m} = B_{1,m} = 0$ 或 $C_{1,m} = S_{1,m} = 0$
- 保留的階數根據觀測資料的精度以及對正常重力位要求的精度確定，一般只需保留至8階項

根據上述性質，該表達式可進一步簡化為：

$$V = {\sum_{n=0}^{8}}\frac{1}{r^{n+1}}A_{n,0}P_n(\cos\theta)$$

或

$$V = \frac{GM}{r}{\sum_{n=0}^{8}}\left(\frac{a}{r}\right)^n C_{n,0}P_n(\cos\theta)$$

### 拉普拉斯方法

利用引力位 $V$ 在邊界面外部滿足拉普拉斯方程的性質，將地球的真實引力位可展開為球諧級數，保留其中的頭幾項作為正常重力位的引力位部分從而確定正常重力位的方法被稱為拉普拉斯方法。通過選取不同的大地水準面重力位值，可以得到不同的正常重力位等位面，從中選取一個最接近於大地水準面的，這一曲面即為產生正常重力位的質體的表面。

#### 重力位的橢球諧級數

將正常重力位直接展開成橢球面 $S_0$ 上的級數，稱為橢球諧級數，形式較球諧函數更為複雜：

$$U = V + \Phi = \sum_{n}^{\infin}{Q_n\left(i{u \over E}\right) \over Q_n\left(i{b \over E}\right)}A_nP_n(\sin\beta) + {1\over2}\omega^2(u^2+E^2)\cos^2\beta$$

上式中各項符號的物理或數學意義與球諧函數有所不同：

- $(u, \beta, \lambda)$ 是空間中某特定點的橢球坐標，$u$ 是該坐標所在橢球的短軸半長， $\beta$ 和 $\lambda\ $分別是該點的歸化緯度和經度
- $a$、$b$ 分別是正常橢球面 $S_0$ 的半長軸和半短軸，$E = \sqrt{a^2-b^2}$ 是橢球的線性偏心率（即半焦距）
- $Q_n$是 $n$ 階第二類勒讓德多項式

當 $u = b$ 時，由重力位 $U_0$ 所決定的等位面應當與正常橢球面 $S_0$ 相重合，此時有：

$$U_0 = V_0 + \Phi_0 = \sum_{n}^{\infin}A_nP_n(\sin\beta) + {1\over2}\omega^2a^2\cos^2\beta$$

#### 引力位的橢球諧級數

若且唯若所有含 $P_n(\sin\beta)$ 的項均為零時，對於任意 $\beta$ 值該公式都成立。對其頭三項進行展開：

$$\left(A_0+{1\over3}\omega^2a^2-U_0\right)P_0(\sin\beta) + A_1P_1(\sin\beta) + \left(A_2-{1\over3}\omega^2a^2\right)P_2(\sin\beta) + \sum_{n=3}^{\infin}A_nP_n(\sin\beta) = 0$$

因此，有

- $$A_0=U_0-{1\over3}\omega^2a^2$$
- $$A_1=0$$
- $$A_2={1\over3}\omega^2a^2$$
- $$A_3=A_4=\cdots=0$$

得到引力位部分的橢球諧級數表達式為：

$$V = \left(U_0-{1\over3}\omega^2a^2\right){Q_0\left(i{u \over E}\right) \over Q_0\left(i{b \over E}\right)} + {1\over3}\omega^2a^2{Q_2\left(i{u \over E}\right) \over Q_2\left(i{b \over E}\right)}P_2(\sin\beta) = \left(U_0-{1\over3}\omega^2a^2\right){\tan^{-1}{E \over u} \over \tan^{-1}{E \over b}} + {1\over3}\omega^2a^2{q \over q_0}P_2(\sin\beta)$$

其中

- $$q={1\over2}\left[\left(1+3\frac{u^2}{E^2}\right)\tan^{-1}\frac{E}{u}-3\frac{u}{E}\right]$$
- $$q_0={1\over2}\left[\left(1+3\frac{b^2}{E^2}\right)\tan^{-1}\frac{E}{b}-3\frac{b}{E}\right]$$

#### 大地水準面重力位

考慮到橢圓坐標 $(u,\beta)$ 與向徑 $r$ 存在如下轉換關係：

$$r^2=u^2+E^2\cos^2\beta$$

利用該關係對下式進行線性化，得：

$$\frac{1}{u}=\frac{1}{r}+\Omicron(\frac{1}{r^3})$$

$$\tan^{-1}\frac{E}{u}=\frac{E}{u}+\Omicron(\frac{1}{u^3})=\frac{E}{r}+\Omicron(\frac{1}{r^3})$$

利用這兩項關係式，可以得到線性化後的引力位函數：

$$V = \left(U_0-{1\over3}\omega^2a^2\right){E \over \tan^{-1}(E/b)}\frac{1}{r} + \Omicron(\frac{1}{r^3})$$

比較與

$$V = \frac{GM}{r} + \Omicron(\frac{1}{r^3})$$

得出大地水準面重力位的表達式：

$$U_0=\frac{GM}{E}\tan^{-1}\frac{E}{b}+{1\over3}\omega^2a^2$$

#### 正常重力位計算公式

將大地水準面重力位的表達式代入原正常重力位的計算公式中，得：

$$U(u,\beta)=\frac{GM}{E}\tan^{-1}\frac{E}{u} +{1\over2}\omega^2a^2\frac{q}{q_0}\left(\sin^2\beta-{1\over3}\right) +{1\over2}\omega^2(u^2+E^2)\cos^2\beta$$

## 物理性質

### 球諧係數

球諧係數的積分公式中包含了正常橢球體內質量的分布關係，且積分範圍是整個正常橢球體，因此球諧係數與正常橢球體的某些物理性質相關。

#### 零階項

當 $n=0$ 時，球諧係數只有一項：

$$A_{0,0} = G \iiint\limits_{\text{earth}} \operatorname{d}\!M' = GM$$

或

$$C_{0,0} = {1 \over M} \iiint\limits_{\text{earth}} \operatorname{d}\!M' = 1$$

即球諧係數的零階項反映了正常橢球體的地心引力常數或總質量。

#### 一階項

當 $n=1$ 時，球諧係數有三項：

$$A_{1,0} = G \iiint\limits_{\text{earth}} z' \operatorname{d}\!M' = Gz_M$$

$$A_{1,1} = G \iiint\limits_{\text{earth}} x' \operatorname{d}\!M' = Gx_M$$

$$B_{1,1} = G \iiint\limits_{\text{earth}} y' \operatorname{d}\!M' = Gy_M$$

或

$$C_{1,0} = {1 \over Ma} \iiint\limits_{\text{earth}} z' \operatorname{d}\!M' = {1 \over a}z_M$$

$$C_{1,1} = {1 \over Ma} \iiint\limits_{\text{earth}} x' \operatorname{d}\!M' = {1 \over a}x_M$$

$$S_{1,1} = {1 \over Ma} \iiint\limits_{\text{earth}} y' \operatorname{d}\!M' = {1 \over a}y_M$$

其中，$\vec{r}_M = \left(x_M,y_M,z_M\right)$ 表示橢球的質心坐標。當坐標系的原點與橢球質心重合時，$x_M=y_M=z_M=0$，所以球諧係數的一階項一般都為零。

#### 二階項

當 $n=2$ 時，球諧係數有五項（僅以 $A_{m,n}$ 和 $B_{m,n}$ 為例）：

$$A_{2,0} = G \iiint\limits_{\text{earth}} {1\over2}(2z'^2-x'^2-y'^2) \operatorname{d}\!M'$$

$$A_{2,1} = G \iiint\limits_{\text{earth}} x'z' \operatorname{d}\!M'$$

$$B_{2,1} = G \iiint\limits_{\text{earth}} y'z' \operatorname{d}\!M'$$

$$A_{2,2} = G \iiint\limits_{\text{earth}} {1\over4}(x'^2-y'^2) \operatorname{d}\!M'$$

$$B_{2,2} = G \iiint\limits_{\text{earth}} {1\over2}x'y' \operatorname{d}\!M'$$

其中的五項積分，既可以四極矩張量 $\mathbf{Q}$ 或慣性張量 $\mathbf{I}$ 表達：

$$\iiint\limits_{\text{earth}} {1\over2}(2z'^2-x'^2-y'^2)\operatorname{d}\!M' = {1\over2}Q_{zz} = {1\over2}(I_{xx}+I_{yy})-I_{zz}$$

$$\iiint\limits_{\text{earth}} x'z'\operatorname{d}\!M' = {1\over3}Q_{xz} = -I_{xz}$$

$$\iiint\limits_{\text{earth}} y'z'\operatorname{d}\!M' = {1\over3}Q_{yz} = -I_{yz}$$

$$\iiint\limits_{\text{earth}} {1\over4}(x'^2-y'^2)\operatorname{d}\!M' = {1\over12}(Q_{xx}-Q_{yy}) = {1\over4}(I_{yy}-I_{xx})$$

$$\iiint\limits_{\text{earth}} {1\over2}x'y'\operatorname{d}\!M' = {1\over6}Q_{xy} = -{1\over2}I_{xy}$$

若且唯若坐標系的各坐標軸與地球的主慣性軸重合時，$I_{xz}=I_{xy}=I_{yz}=0$，因此亦有 $A_{2,1}=B_{2,1}=B_{2,2}=0$。反過來，坐標軸的選擇又決定了這三個二階項係數的值：當坐標系的 $Z$ 軸指向協議地球極時，受極移等因素的影響，這一指向與地球的瞬時主慣性軸並不重合，因此 $A_{2,1}$ 與 $B_{2,1}$ 的並不為零；而 $X$ 軸的指向（通常是本初子午線）則決定了 $B_{2,2}$ 的數值。另外，$A_{2,2}$ 由赤道的形狀決定：當正常橢球為對稱的旋轉體時，赤道是圓形，此時 $I_{xx}=I_{yy}$，即有$A_{2,2}=0$。