大地水準面

大地水準面（德語：Geoid）是指地球重力場^{[註 1]}中，與處於自由靜止狀態的平均海水面相重合^[3]^:49或最為接近^[4]^:42的重力等位面。這一概念最早由德國大地測量學家卡爾·弗里德里希·高斯在1828年提出。當時，高斯以「地球的數學表面」^[3]^:73來指稱與重力方向相垂直、且與靜止的平均海水面相重合的幾何表面^[5]^:4，並提出將其作為高程系統的基準面^[6]^:74。其後，高斯的學生利斯廷於1873年創造出了「Geoid」一詞，用以描述高斯所提出的數學表面。^[7]

在大地測量學中，大地水準面被視作是地球的物理形狀和數學形狀。^[4]^:2由於自然的地形表面形態過於複雜，大地測量學通常是將重力場中整體形狀與地形表面最為接近的等位面作為地球的形狀進行研究。^[8]^:226處於靜力平衡狀態下的平均海水面被視作是符合這一標準的重力等位面，這一假想的海水面不受潮汐、風浪及大氣壓的變化影響，僅在地球引力和因地球自轉產生的離心慣性力的作用下保持平衡。^[9]^[4]^:41將該平均海水面所處的重力等位面延伸到陸地內部，形成的閉合曲面即為大地水準面，其所包圍的形體又被稱為大地體。^[10]^:29

1849年，英國物理學家斯托克斯提出了計算大地水準面的斯托克斯方法，使大地水準面的形狀能夠通過其整個表面上的重力觀測值確定。^[10]^:8但在傳統的大地測量中，通過大範圍的重力觀測，以及對這些重力觀測值進行積分而求得大地水準面的方式遇到了諸多困難。^[8]^:294-295隨着20世紀中後期衛星重力測量等技術逐漸發展成熟，採用球諧級數表達的地球重力場模型逐漸成為了描述大地水準面的主要方式^[11]，全球尺度下的大地水準面也得以以分米級的精度測定^[12]。莫洛堅斯基、布耶哈馬等大地測量學家，亦在斯托克斯方法的基礎上提出了通過地面測量技術確定大地水準面的新方法。^[5]^:34

大地水準面是測量外業所依據的基準面，在測量學中具有重要地位。^[13]在各類高程系統中，日常中稱為海拔高程的正高系統是基於大地水準面建立的。^[10]^:42如何確定大地水準面的形狀，是物理大地測量學研究的關鍵問題之一。^[14]^[15]^[16]

定義

經典定義

1873年，利斯廷提出了大地水準面的經典定義，即大地水準面是與平均海水面相重合的重力等位面。^[7]在1975年，美國大地測量學家理查德·拉普（英語：Richard H. Rapp）進一步完善了這一定義，他將重力位 $W$ 視為除大氣層外的地球質量所引起的重力位 $W_{g}$ 、大氣質量所引起的重力位 $W_{a}$ ，以及受日、月的引潮力作用而引起的重力位 $W_{t}$ 的和：^[17]

W=W_{g}+W_{a}+W_{t}

不同於 $W_{g}$ ，後兩項與地球內部的質量分佈無關，因此可以被相對精確地計算。拉普認為，定義大地水準面的重力位應當僅包含 $W_{g}$ 和 $W_{a}$ ，即大地水準面應當是使 $W_{g}+W_{a}$ 等於常數的重力等位面，而因引潮力、風和洋流等其他因素影響產生的變化則不在考慮範圍之內。這種定義方式認為，通過多年的驗潮數據可以得到一個平均海水面，其上的任意一點的位置與海洋面的長期平均位置相重合。^[18]

然而，20世紀60年代以來對平均海水面的研究發現，其與大地水準面間存在最大可達2米的系統性的偏移^[19]，這就使大地水準面的經典定義暴露出了缺陷。因此，馬瑟（英語：Mather R. S.）等人在1978年提出，大地水準面應重新定義為與平均海水面最為密合的重力等位面。^[20]

相對論定義

1985年，瑞典大地測量學家布耶哈馬基於相對論重力提出了相對論大地水準面的概念，即一個使其上所有精密時鐘的運行速率相同且與平均海水面最為密合的封閉曲面。^[21]在2019年，德國不來梅大學的丹尼斯·菲利普（德語：Dennis Philipp）利用新的測量結果，以及無含時紅移勢能（英語：time-independent redshift potential）等新的理論，對相對論大地水準面的概念給出了數學上的表述方法。這一方法給出的相對論水準面與經典大地水準面之間存在着毫米級的差異。^[22]

形狀

從數學上看，大地水準面是一個連續但不規則的閉合曲面，它與經過這一曲面的鉛垂線處處正交。^[23]^:49地球內部的質量分佈，特別是外層的質量分佈的不均勻性，使得大地水準面的形狀變得特別複雜。^[10]^:29在發生密度突變的區域，大地水準面的曲率會表現出不連續性。因此，大地水準面並非是解析曲面。^[4]^:41

大地水準面的形狀可以用橢球面近似，該橢球面也被稱為參考橢球面。以WGS84坐標系使用的橢球為例，其長半軸位於赤道面上，長度約6,378,137米；短半軸指向極點，長度約6,356,752米，比長半軸短約21千米。^[24]相對於WGS84橢球，大地水準面的起伏從最高點到最低點的範圍僅不到200米。在這一差距被放大的情況下，大地水準面的形狀則與馬鈴薯或梨相似。^[25]儘管大地水準面是一個不規則的曲面，但它仍比真實的地形表面更加平滑，後者的高程差跨越了近20千米。

從赤道觀察大地水準面（紅線）相對於參考橢球面（黑線）的形狀，兩者間的差值為同一緯線上大地水準面差距的平均值

從地極處觀察大地水準面（紅線）相對於參考橢球面（黑線）的形狀，兩者間的差值為同一經線上大地水準面差距的平均值

起伏

大地水準面起伏又稱大地水準面高或大地水準面差距，是指大地水準面上的某一點，與其沿法線投影到參考橢球面上的相應位置之間的距離。^[23]^:83大地水準面高的正極值出現在新幾內亞，為+85.4米；負極值出現在印度南部，為−107.0米。^[26]

大地水準面起伏的幅度和空間分佈受地球內部密度分佈的直接影響。在頻譜分析中，大地水準面起伏的分佈可分為長波（尺度為上千公里，幅度為±100米左右）、中波（尺度為上百公里，幅度為±10米）和短波（尺度為數十公里，幅度為±1米）部分。^[26]其中，長波部分反映地幔及地球內部的密度異常分佈^[27]，而中短波部分則反映岩石圈和地形的影響。在地球物理學研究中，大地水準面的起伏又可分為地形的不規則部分、結晶基底起伏、莫霍面起伏、岩石圈底面起伏，以及因地幔和地核部分的橫向密度不均勻造成的起伏等。^[28]

正異常

大地水準面高在北大西洋、歐洲、非洲的西部和南部、印度洋西南部、安地斯山脈、西太平洋及赤道太平洋部分等地表現出正值^[29]^[30]，世界上所有的大型地盾、顯生宙地台和地質熱點幾乎也都位於這些地區^[31]。其中，分佈較廣而幅度較高的有兩個區域，橫跨了大西洋和非洲的正異常區域被稱為大西洋-非洲大地水準面高值區（英語：Atlantic-African geoid high），太平洋中西部的正異常區域則被稱為赤道太平洋大地水準面高值區（英語：Equatorial Pacific geoid high）。西太平洋的正異常區域在中生代早期曾是俯衝雜岩帶所處的位置。其他與俯衝帶不相關的正異常區域則可能曾是大陸板塊的聚集處，或是古板塊的中心等會發生地幔絕熱現象的地區。地幔的絕熱會使得其本身的溫度升高，其產生的抬升和岩漿作用能夠產生大地水準面的正異常。^[31]

負異常

大地水準面高在印度及印度洋北部，以及南極洲附近表現出負值。其中又以印度洋北部的區域表現得尤其顯著，該區域也被稱為印度洋大地水準面低值區（英語：Indian Ocean geoid low）。^[32]對於這些負異常形成的原因，存在有眾多推測，如軟流層中的淺層物質流^[32]、板塊俯衝在中地幔引起的地幔對流^[33]，以及地幔和地核的邊界處發生的凹陷^[34]等。同時，這些推測也都同意，印度洋大地水準面低值區出現的部分原因是特提斯洋在中生代的俯衝帶來的密度異常。^[35]

測定

測定大地水準面進行的技術手段有多種，包括在地面、船舶和飛機中直接進行重力測量，或是通過天文方法和衛星導航系統進行水準測量以獲取垂線偏差或大地水準面高，又或是利用雷達和人造衛星進行衛星測高來描繪地形表面的高度，以及通過對衛星軌道攝動情況的分析、衛星跟蹤衛星以及衛星重力梯度測量等衛星重力探測技術直接反演地球重力場的中長波分量等。^[36]^[37]^[38]^[6]^:318-319其所依據的原理是對地面的重力測量數據進行積分（如斯托克斯方法），或將地球重力場以球諧級數的形式表達（如拉普拉斯方法）。^[39]

在人造衛星出現前，斯托克斯方法是測定大地水準面的主要原理，但這種方式在應用過程中存在諸多困難：斯托克斯公式中的積分難以收斂、全球性的觀測數據難以獲得等。^[11]1957年斯普特尼克衛星的升空，揭示了通過人造衛星反演地球重力場的可能性。^[40]從此，使用拉普拉斯方法構建全球的重力場模型，成為了大尺度上表達大地水準面的主要方式。^[11]21世紀後，諸如CHAMP、GOCE和GRACE等新一代重力測量衛星的出現，得以使對地球重力場的研究以空前的精度和空間解像度進行，還使得對地球重力場的時變監測成為可能。^[41]^[42]

斯托克斯方法

喬治·斯托克斯在1849年提出的斯托克斯定理指出：^[43]在大地水準面的外部沒有質量分佈的情況下，其形狀和其外部的重力場可以由地球的質量、自轉角速度以及大地水準面上的重力觀測值唯一確定，而與其內部的質量分佈無關。^[23]^:17這一定理說明，大地水準面的形狀可以通過地面上的重力觀測數據確定。而利用斯托克斯方法確定大地水準面的形狀，需要確定大地水準面高相對於參考橢球面的高程差（即大地水準面高）以及這一高程差在參考橢球面上的位置。^[8]^:284這分別可以用斯托克斯公式和維寧·曼尼茲公式進行計算。

值得注意的是，實際上的重力觀測數據是在地形表面，而非大地水準面上得到的。因此，在得到地表的觀測數據後，還需要除去大地水準面外的質量對重力觀測值的影響，並將地面上重力觀測值的位置改正到大地水準面上，這一部分的內容又被稱為重力歸算。^[8]^:264通過斯托克斯方法，可以求得較高精度的高程異常並隨之確定相應的大地水準面高，常用於局部重力場模型的建立。^[23]^:329

斯托克斯公式

斯托克斯公式給出了大地水準面高與全球大地水準面上的重力異常及其分佈的數學關係：^[43]

N={R \over 4\pi {\bar {\gamma }}}\iint \limits _{\sigma }\Delta g\,S(\psi )\operatorname {d} \!\sigma

上式中， $N$ 是大地水準面高， $R$ 是地球的平均半徑^[23]^:87， ${\bar {\gamma }}$ 是全球正常重力的平均值。而對於積分部分， ${\begin{smallmatrix}\iint \limits _{\sigma }\end{smallmatrix}}$ 表示在整個單位球面上的積分， $\operatorname {d} \!\sigma$ 是立體角元素。^[23]^:89 $\Delta g$ 是單位角元素上的重力異常，由重力測量數據給出的重力位減去計算得到的正常重力位後得到。 $S(\psi )$ 則為斯托克斯函數^[44]，該項是單位球面上的被計算點與重力異常觀測值所在的角元素之間的夾角 $\psi$ 的函數。^[23]^:94

維寧·曼尼茲公式

在由斯托克斯公式求得大地水準面相對於參考橢球面的高程差後，還需通過垂線偏差改正，將大地水準面上的天文坐標轉換成參考橢球面上的地理坐標。這項改正能夠通過通過維寧·曼尼茲公式計算：^[23]^:114

{\begin{Bmatrix}\xi \\\eta \end{Bmatrix}}={1 \over 4\pi {\bar {\gamma }}}\iint \limits _{\sigma }\Delta g\,{\begin{Bmatrix}\cos {\alpha }\\\sin {\alpha }\end{Bmatrix}}\,S(\psi )\operatorname {d} \!\sigma

上式中的 $\xi$ 表示天文坐標相對於地理坐標的緯度差， $\eta$ 表示天文坐標相對於地理坐標的經度差， $\alpha$ 則是天文坐標相對於地理坐標的方位角。

拉普拉斯方法

利用引力位在質體外部滿足拉普拉斯方程的特性^[23]^:5，可以將其展開成球諧級數進行逼近，球諧級數中的各項係數由地球的各項物理性質確定。這些物理性質，如地球的扁率和重力的不規則性等，會對人造衛星在空間中的運動造成影響，即衛星軌道的攝動。因此通過對衛星軌道進行動力學分析，可以反演出地球的重力場模型。^[45]^:149又由於人造衛星對重力場模型中的高階項並不敏感^[6]^:319，人造衛星主要用於獲得重力場模型中的低階項，而高階項還需結合地面的重力觀測數據等其他方式獲得。^[16]

拉普拉斯方程的球諧級數解

通過求解球面坐標中的拉普拉斯方程 $\Delta V=0$ ，單位球外部（ $r>1$ ）的引力位可以展開為：^[23]^:59

V={\sum _{n=0}^{\infty }}{\sum _{m=0}^{n}}\left[A_{nm}{R_{nm}(\theta ,\lambda ) \over r^{n+1}}+B_{nm}{S_{nm}(\theta ,\lambda ) \over r^{n+1}}\right]

上式中各個量的含義如下：

$(r,\theta ,\lambda )$ 是空間中某特定點的球坐標， $r\$ 是點的向徑， $\theta$ 和 $\lambda \$ 分別是點的極角和方位角

$A_{nm}$ 和 $B_{nm}$ 是通過對單位球上的觀測數據積分確定的待定係數

$R_{nm}(\theta ,\lambda )=P_{nm}(\cos \theta )\cos {m\lambda }$ ， $S_{nm}(\theta ,\lambda )=P_{nm}(\cos \theta )\sin {m\lambda }$ ，是與締合勒讓德多項式有關的面諧函數^[23]^:29

上式中的向徑 $r$ 是被計算點相對於單位球的相對半徑，而非其相對於地心的真實距離。

衛星動力學中的球諧級數

在衛星動力學中，球諧級數常寫成如下形式：^[23]^:59

V={GM \over r}\{{1-\sum _{n=1}^{\infty }}{\sum _{m=0}^{n}}{({a \over r})}^{n}\left[J_{nm}R_{nm}(\theta ,\lambda )+K_{nm}S_{nm}(\theta ,\lambda )\right]\}

上式中各個量的含義如下：

$(r,\theta ,\lambda )$ 是空間中某特定點的球坐標， $r\$ 是空間中點的地心距離， $\theta$ 和 $\lambda \$ 分別是點的地心緯度的餘角和經度

$GM$ 為重力場模型的地心引力常數
$J_{nm}$ 和 $K_{nm}$ 是與地球動力學有關的參數，其與 $A_{nm}$ 和 $B_{nm}$ 的關係為

{\begin{cases}A_{nm}=-GMa^{n}J_{nm}\\B_{nm}=-GMa^{n}K_{nm}\end{cases}}

在這一表達式中，參數 $A_{nm}$ 、 $B_{nm}$ 、 $J_{nm}$ 和 $K_{nm}$ 均由地球本身的特性決定，包含了地球的轉動慣量和動力學形狀等信息，是引力場模型中的待解參數。^[23]^:61-63

模型

重力場模型

大地水準面是通過重力場模型進行表達的，對全球的重力場進行描述的模型即為全球重力場模型。^[46]在建立重力場模型的過程中，重力位 $W$ 被分成引力位 $V$ 和離心力位 $\Phi$ 兩部分（即 $W=V+\Phi$ ），其中離心力位 $\Phi$ 能夠根據地球的自轉角速度 $\omega$ 和點在空間中的坐標 $(x,y,z)$ 直接計算得到：^[23]^:47

\Phi ={1 \over 2}{\omega }^{2}(x^{2}+y^{2})

引力場位的確定比離心力位的確定要複雜得多，引力場模型的建立是重力場模型建立過程中的主要部分。因此，在不嚴謹的表述中，有時也將引力場模型直接稱為重力場模型。^[46]

球諧表達式

在引力場模型中，引力位 $V$ 在地球外部（ $r>a$ ）表達為如下形式：^[47]

V={\frac {GM}{r}}{\sum _{n=0}^{n_{\text{max}}}}\left({\frac {a}{r}}\right)^{n}{\sum _{m=0}^{n}}{\overline {P}}_{nm}(\sin \phi )\left[{\overline {C}}_{nm}\cos m\lambda +{\overline {S}}_{nm}\sin m\lambda \right]

上式中的各個量的含義如下：

$(r,\phi ,\lambda )$ 是空間中某特定點的球坐標， $r\$ 是點的地心距離， $\phi \$ 和 $\lambda \$ 分別是點的地心緯度和經度

$GM$ 為重力場模型的地心引力常數
$n_{\text{max}}$ 為該重力場模型的最大階數

$a$ 為地球（參考橢球體）的長半軸
${\overline {P}}_{nm}$ 是 $n\$ 階 $m\$ 次的完全正規化締合勒讓德多項式
${\overline {C}}_{nm}$ 和 ${\overline {S}}_{nm}$ 是由測量數據確定的該重力場模型的完全正規化係數

其中，當 $n=0$ 時， ${\overline {P}}_{00}(\sin \phi )=1$ ， ${\overline {C}}_{00}=1$ ，則求和符號右側所有項的和為 $1$ ；當 $n=1$ 時， ${\overline {C}}_{1,m}={\overline {S}}_{1,m}=0$ ，則求和符號右側所有項的和為 $0$ . ^[47]因此，上式可進一步表達為如下形式：

V={\frac {GM}{r}}\left(1+{\sum _{n=2}^{n_{\text{max}}}}\left({\frac {a}{r}}\right)^{n}{\sum _{m=0}^{n}}{\overline {P}}_{nm}(\sin \phi )\left[{\overline {C}}_{nm}\cos m\lambda +{\overline {S}}_{nm}\sin m\lambda \right]\right)

在得到引力位 $V$ 和橢球面坐標為 $(r,\phi ,\lambda )$ 的關係後，再將其與離心力位模型 $\Phi$ 相加，即可得到重力位模型 $W$ ，最後得到表達為 $W(r,\phi ,\lambda )=W_{0}$ 的大地水準面的形狀。^[23]^:48 在實際的模型建立的過程中，也常取引力位中的數階最大項作為正常重力場中的引力位，再將餘下的引力位以擾動位的形式進行表達，而正常重力場中的等位面即為參考橢球面。^[8]^:207-212

空間解像度

根據奈奎斯特採樣定理，球諧函數表達式中的階次越高，其表達的範圍也就越精細。對於以球諧函數來表達的重力場模型，其空間解像度通常以最短半波長表示。在地球表面，模型的空間解像度 $\psi _{\text{min}}$ 與最大階數 $n_{\text{max}}$ 的關係為：^[46]

\psi _{\text{min}}={\pi R \over n_{\text{max}}}\approx {20\,000\,{\text{km}} \over n_{\text{max}}}

其中的 $R$ 為地球的平均半徑。通過這一公式，可以求出現有重力場模型常採用的空間解像度：


最大階數	最短半波長 / km
100	200.00
200	100.00
360	55.55
2190	9.13

在美國國家地理空間情報局（NGA）發佈的EGM96模型中， $n_{\text{max}}=360$ ，即該模型給出了最高完全階次至360的重力場模型係數，其表達的大地水準面的解像度可達55千米或110千米（取決於對解像度定義的不同）。而在2009年，NGA又發佈了最大階數 $n_{\text{max}}=2190$ ，最高完全階次為2159次的EGM08模型，其大地水準面的解像度可達12千米，精度在全球範圍內可達±15厘米。^[48]^[12]

現有的地球重力場模型

現有的地球重力場模型由國際地球重力場模型中心（ICGEM）負責收集和存檔。^[49]截至2020年，ICGEM存儲的階數超過2000階的超高階靜態全球重力場模型有如下幾個：


模型	發佈時間	最高階數	數據來源	參考文獻
EGM2008	2009年	2190	A、G、S(Grace)	^[48]
EIGEN-6C4	2014年	2190	A、G、S(Goce、Grace、Lageos)	^[50]
GECO	2015年	2190	EGM2008、S(Goce)	^[51]
SGG-UGM-1	2018年	2159	EGM2008、S(Goce)	^[52]
XGM2019e_2159	2019年	760/2190/5540	A、G、S(GOCO06s)、T	^[53]
註：其中A表示測高數據，S表示為衛星測量數據，G表示地表測量數據（包括陸地、船舶和航空測量），T表示測圖數據

高程系統

大地水準面是正高系統的參考面，正高被定義為地面上任意一點沿其所在的鉛垂線到大地水準面的距離。^[23]^:166在重力水準測量中，觀測到大地水準面上一點 $O$ 與任意一點 $B$ 之間的重力位差 ${\begin{smallmatrix}\int \limits _{OB}g\operatorname {d} \!h\end{smallmatrix}}$ ，再求得點 $B$ 沿鉛垂線到大地水準面的路徑上的重力平均值 $g_{\text{m}}^{B}$ ，可計算得點 $B$ 的正高為：^[10]^:42

H^{B}={1 \over g_{\text{m}}^{B}}\int \limits _{OB}g\operatorname {d} \!h

一個正高系統必須與一個重力位已知的大地水準面相聯繫，這個重力位即為大地水準面重力位。

與大地高的關係

在大地測量學以外的領域，正高也被稱為海拔高度。^[23]^:172而在進行空間位置的計算（如計算GPS等全球衛星導航系統中的衛星位置）時，使用的高度則是相對於參考橢球面的距離，稱為大地高或橢球高。兩者間的轉換需要通過大地水準面高進行：^[54]

h=H+N

其中， $h$ 表示大地高（橢球高）， $H$ 表示正高， $N$ 表示大地水準面高。

與正常高的關係

正高的計算需要用到重力平均值 $g_{\text{m}}^{B}$ ，但該值無法在地球內部的質量和重力分佈未知的情況下得到，所以正高無法被精確求得。^[8]^:3511945年，蘇聯大地測量學家莫洛堅斯基提出了利用正常高來代替正高的計算方式。^[55]在這一高程系統中，由地球表面向下量取正常高而獲得的表面被稱為似大地水準面，由正常橢球面向上量取正常高而獲得的表面則被稱為似地形表面，而似大地水準面和正常橢球面之間的差距又被稱為高程異常。^[56]似大地水準面與水準面不同，它沒有具體的物理意義，更不是重力等位面，只是用於計算的輔助面。^[23]^:294大地水準面和似大地水準面在山區的差異可達4米，但兩者在平均海水面上是重合的。^[8]^:352