正常重力位

正常重力位（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$。