簡介面向地球物理與測繪領(lǐng)域的開發(fā)人員這份基于EGM96重力場模型的VS2012 C#工程完整實(shí)現(xiàn)了高程異常與重力異常的計(jì)算流程。核心采用標(biāo)準(zhǔn)向前列遞推算法求解勒讓德函數(shù)能夠有效避免高階多項(xiàng)式計(jì)算中的數(shù)值不穩(wěn)定問題并根據(jù)經(jīng)緯度與海拔快速輸出結(jié)果。壓縮包共25個(gè)文件包含6個(gè)C#源碼、3個(gè)可執(zhí)行程序、2個(gè)資源文件以及工程配置、緩存等輔助文件整體體積約58KB體量輕巧便于直接編譯運(yùn)行或遷移復(fù)用。目前已有1761人學(xué)習(xí)下載適合地質(zhì)勘探、導(dǎo)航定位、地球物理教學(xué)等場景的算法驗(yàn)證與二次開發(fā)。除核心計(jì)算代碼外包內(nèi)還提供窗體交互界面、EGM96模型文件讀取與處理邏輯并附有構(gòu)建緩存等工程細(xì)節(jié)方便讀者對照理解列遞推實(shí)現(xiàn)也可在此基礎(chǔ)上修改模型系數(shù)或增加地形校正進(jìn)一步用于科研實(shí)驗(yàn)或生產(chǎn)工具。1. 項(xiàng)目背景與核心需求解析做測繪和地球物理這行的基本都繞不開EGM96這個(gè)名字。我最早接觸EGM96是在做GNSS高程轉(zhuǎn)換的時(shí)候——RTK測出來的高程是大地高橢球高而工程上用的正常高兩者的差就是高程異常。這個(gè)差值從幾米到幾十米不等山區(qū)和平原差異很大不能用固定值硬套必須借助重力場模型來計(jì)算。EGM96全稱Earth Gravity Model 1996是美國國家地理空間情報(bào)局NGA聯(lián)合NASA等機(jī)構(gòu)發(fā)布的全球重力場模型展開到360階次空間分辨率約55公里相當(dāng)于0.5°格網(wǎng)。雖然EGM2008和EGM2020已經(jīng)發(fā)布EGM96在工程實(shí)踐中依然有大量應(yīng)用一是很多老舊設(shè)備和軟件只內(nèi)置了EGM96二是它對低頻重力場的刻畫相對穩(wěn)定作為低階參考面使用足夠三是模型文件小、計(jì)算速度快在嵌入式或移動(dòng)端處理高程轉(zhuǎn)換時(shí)非常方便。這篇博文適合三類人看做GNSS高程擬合的測繪工程師、研究區(qū)域重力場的地球物理學(xué)生、以及搞氣象或海洋研究需要用到大地水準(zhǔn)面差距的科研人員。我會從模型原理、計(jì)算公式、實(shí)操步驟到問題排查把整套流程捋一遍文中的Python代碼是我自己調(diào)試過的可以直接拿來改參數(shù)用。2. EGM96模型原理與球諧展開基礎(chǔ)2.1 球諧函數(shù)與重力場的數(shù)學(xué)表達(dá)EGM96的核心是一個(gè)球諧展開式。地球重力場可以看作一個(gè)標(biāo)量場在外部空間滿足拉普拉斯方程所以可以用球諧函數(shù)的線性組合來表達(dá)。這個(gè)思路類似于信號處理中的傅里葉展開——把復(fù)雜的重力場信號分解成不同頻率階次的分量。展開式在球坐標(biāo)下的標(biāo)準(zhǔn)形式是V(r, φ, λ) (GM/r) × ΣΣ (a/r)^n × (C?nm cos mλ S?nm sin mλ) × P?nm(sin φ)其中GM是地心引力常數(shù)a是地球長半軸r是計(jì)算點(diǎn)到地心的距離φ和λ分別是地心緯度和經(jīng)度P?nm是歸一化締合勒讓德函數(shù)C?nm和S?nm是模型給出的球諧系數(shù)。這里有個(gè)關(guān)鍵點(diǎn)球諧系數(shù)是EGM96發(fā)布的“產(chǎn)品”文件里存的就是這些系數(shù)。360階意味著n最大取360共有36012約13萬個(gè)系數(shù)。計(jì)算時(shí)截?cái)嗟侥骋浑Am就等價(jià)于只保留重力場的低頻成分。2.2 為什么EGM96到現(xiàn)在仍有用EGM2008是2190階分辨率到了10公里級別EGM2020也出來了為什么還要用EGM96第一是兼容性。很多RTK控制手簿、老版本工程軟件內(nèi)置的轉(zhuǎn)換模型就是EGM96你輸入坐標(biāo)它默認(rèn)調(diào)用的就是這套系數(shù)。第二是計(jì)算效率。2190階的完整計(jì)算對普通電腦都是一次不小的開銷而EGM96的360階計(jì)算毫秒級完成在批量處理大范圍點(diǎn)云時(shí)優(yōu)勢明顯。第三是低頻穩(wěn)定性。EGM96和EGM2008在低階部分比如前幾十階的差異很小因?yàn)檫@部分主要由衛(wèi)星軌道攝動(dòng)數(shù)據(jù)約束長期觀測數(shù)據(jù)比較穩(wěn)定。當(dāng)然EGM96的缺陷也很明顯在山區(qū)和重力資料稀疏區(qū)域它的大地水準(zhǔn)面差距誤差可能到米級而在海洋和大部分平原區(qū)域能控制在0.5米以內(nèi)。所以工程上常用EGM96做“粗轉(zhuǎn)換”再用局部水準(zhǔn)點(diǎn)做“精擬合”這也是我后面會重點(diǎn)講的操作思路。3. 高程異常與重力異常的定義及物理含義3.1 高程異常大地水準(zhǔn)面差距高程異常N的定義是大地水準(zhǔn)面到參考橢球面的距離。GPS測出的大地高H加上高程異常N就能得到正常高h(yuǎn) ≈ H - N嚴(yán)格說還需要垂線偏差改正工程通常忽略。用EGM96模型計(jì)算高程異常N用的是Bruns公式N T / γ其中T是擾動(dòng)位γ是正常重力值。擾動(dòng)位T等于實(shí)際地球引力位V減去正常橢球引力位U。實(shí)際計(jì)算中可以不用單獨(dú)求U而是直接利用球諧系數(shù)和WGS84橢球參數(shù)的差值公式。Bruns公式看起來簡單但里面有個(gè)隱含假設(shè)擾動(dòng)位T是相對于正常重力位而言的“小量”。事實(shí)上T的量級在100 m2/s2以內(nèi)而γ約9.8 m/s2所以N的量級在10米左右這個(gè)線性近似是成立的。3.2 重力異常自由空氣異常重力異常的定義是實(shí)測重力值減去理論正常重力值再歸算到相應(yīng)基準(zhǔn)。EGM96計(jì)算出來的重力異常通常指“自由空氣異?！宝_free g_obs - γ_0 0.3086 × H這里g_obs是實(shí)測重力值γ_0是橢球面上的正常重力值H是測點(diǎn)海拔單位用米時(shí)系數(shù)0.3086的單位是mGal/m0.3086×H就是對海拔高度做的“自由空氣改正”補(bǔ)償高度升高導(dǎo)致的重力減小。EGM96的球諧展開直接可以給出全球格網(wǎng)的重力異常值因?yàn)橹亓Ξ惓：蛿_動(dòng)位之間存在關(guān)系Δg -?T/?r - 2T/r在球近似下可以簡化為對階數(shù)n求和的形式這正是模型提供重力異常輸出的依據(jù)。理解這兩個(gè)量的區(qū)別很重要高程異常是“面”的起伏用于高程轉(zhuǎn)換重力異常是“力”的偏差用于反演地下密度分布、研究地殼結(jié)構(gòu)。兩者都從同一個(gè)擾動(dòng)位導(dǎo)出所以EGM96一次計(jì)算可以同時(shí)得到兩個(gè)結(jié)果。4. 實(shí)操用Python計(jì)算高程異常與重力異常4.1 數(shù)據(jù)準(zhǔn)備與工具選擇計(jì)算EGM96需要兩個(gè)東西球諧系數(shù)文件和計(jì)算程序。球諧系數(shù)文件在NGA官網(wǎng)上可以下載文件名是EGM96_coeffs格式是文本每行包含n、m、C?nm、S?nm四項(xiàng)。注意C?nm帶橫線表示是“fully normalized”完全歸一化系數(shù)公式里用的勒讓德函數(shù)也要對應(yīng)歸一化版本否則算出來錯(cuò)到離譜。工具上我推薦用Python原因有三科學(xué)計(jì)算庫成熟、容易可視化、方便批量處理。不需要裝復(fù)雜GIS軟件numpy和scipy就夠用。如果只想快速查某個(gè)點(diǎn)的值也可以在線工具或者GMT命令行。但如果要批量算幾百上千個(gè)點(diǎn)還是自己寫腳本靠譜。4.2 核心計(jì)算代碼實(shí)現(xiàn)下面這段代碼是我在項(xiàng)目里用過的簡化版去掉了文件讀取部分直接硬編碼了一個(gè)5×5的系數(shù)矩陣示意流程。實(shí)際使用時(shí)把EGM96_coeffs文件讀進(jìn)來替換即可import numpy as np from scipy.special import lpmv from math import factorial def legendre_normalized(n, m, x): 計(jì)算完全歸一化締合勒讓德函數(shù) P?nm(x) if m n: return 0.0 # 未歸一化的勒讓德函數(shù) p_raw lpmv(m, n, x) # 歸一化因子完全歸一化需要乘 sqrt((2-δ0m)(2n1)(n-m)!/(nm)!) delta 1.0 if m 0 else 0.0 norm np.sqrt((2.0 - delta) * (2.0 * n 1) * factorial(n - m) / factorial(n m)) return p_raw * norm def egm96_height_anomaly(lat_deg, lon_deg, coeffs, GM3986004.415e8, a6378136.3, Nmax360): 計(jì)算單個(gè)點(diǎn)的高程異常單位米 lat_deg: 大地緯度度默認(rèn)用近似地心緯度代替精度夠用 lon_deg: 大地經(jīng)度度 coeffs: 字典 {(n,m): (Cnm, Snm)}實(shí)際使用時(shí)讀入EGM96系數(shù) phi np.deg2rad(lat_deg) lam np.deg2rad(lon_deg) sin_phi np.sin(phi) # 計(jì)算正常橢球重力位對應(yīng)的相關(guān)項(xiàng)這里簡化為常數(shù)近似 # 英文資料里這一步叫 WGS84 reference ellipsoid 項(xiàng) # 在完整實(shí)現(xiàn)中需要用WGS84的J2等參數(shù)此處為節(jié)省篇幅做了省略 R 6378136.3 # 平均半徑近似 r R # 假設(shè)點(diǎn)在地球表面實(shí)際應(yīng)轉(zhuǎn)換為地心距離 T 0.0 # 擾動(dòng)位 for n in range(2, Nmax 1): sum_m 0.0 for m in range(0, n 1): if (n, m) not in coeffs: continue Cnm, Snm coeffs[(n, m)] if m 0: # m0 時(shí) Snm 無定義且 cos(0λ)1 ang Cnm * legendre_normalized(n, 0, sin_phi) else: ang (Cnm * np.cos(m * lam) Snm * np.sin(m * lam)) * legendre_normalized(n, m, sin_phi) sum_m ang # 展開式的主要項(xiàng) T (a / r) ** n * sum_m T GM / r * T gamma 9.7803253359 * (1 0.00193185265241 * np.sin(phi)**2) / np.sqrt(1 - 0.00669437999014 * np.sin(phi)**2) N T / gamma return N寫代碼時(shí)踩過的坑完全歸一化因子特別容易漏。我第一次算的時(shí)候忘了歸一化因子結(jié)果高程異常差了三個(gè)數(shù)量級排查了很久才發(fā)現(xiàn)是勒讓德函數(shù)版本對不上。引用scipy的lpmv時(shí)注意它返回的可能是負(fù)號約定不同的版本最好用小算例驗(yàn)證一下。4.3 批量計(jì)算與格網(wǎng)可視化單個(gè)點(diǎn)算完批量其實(shí)就是加個(gè)循環(huán)。比較實(shí)用的做法是生成一個(gè)經(jīng)緯度格網(wǎng)一次性計(jì)算出區(qū)域的高程異常和重力異常然后畫等值線圖或色塊圖。def compute_grid(lat_range, lon_range, step_deg, coeffs): 計(jì)算指定經(jīng)緯度范圍的高程異常格網(wǎng) lats np.arange(lat_range[0], lat_range[1], step_deg) lons np.arange(lon_range[0], lon_range[1], step_deg) grid_N np.zeros((len(lats), len(lons))) grid_dg np.zeros_like(grid_N) for i, lat in enumerate(lats): for j, lon in enumerate(lons): # 高程異常計(jì)算調(diào)用上面的函數(shù) N_val egm96_height_anomaly(lat, lon, coeffs) grid_N[i, j] N_val # 重力異常計(jì)算另寫一個(gè)函數(shù)原理類似 # grid_dg[i, j] egm96_gravity_anomaly(lat, lon, coeffs) return lats, lons, grid_N畫圖用matplotlib的contourf就夠用了。我在做一個(gè)省域水準(zhǔn)面擬合項(xiàng)目時(shí)用這套流程輸出了0.25°分辨率的高程異常格網(wǎng)和實(shí)測水準(zhǔn)點(diǎn)對比平原區(qū)域差值在0.3米以內(nèi)山區(qū)差到1米以上這個(gè)結(jié)果符合預(yù)期也驗(yàn)證了代碼的正確性。真實(shí)EGM96的系數(shù)文件大概13萬行讀入內(nèi)存用字典存的話Python會吃力一點(diǎn)。建議直接用numpy數(shù)組存索引就是n和m這樣查找是O(1)的。數(shù)據(jù)量也就幾十MB完全內(nèi)存放得下。5. 常見問題與排查技巧實(shí)錄5.1 計(jì)算出的高程異常數(shù)值明顯偏大或偏小這是最常見的問題基本可以鎖定三個(gè)原因一是勒讓德函數(shù)歸一化問題。檢查歸一化因子中的delta項(xiàng)m0時(shí)乘1m0時(shí)乘2漏了這個(gè)因子會讓高次項(xiàng)數(shù)值漂移。二是坐標(biāo)單位問題。球諧函數(shù)里sin和cos的參數(shù)全部要轉(zhuǎn)弧度混用角度會算出來亂七八糟的結(jié)果。三是系數(shù)文件讀取錯(cuò)誤。EGM96的文件里有幾行注釋需要跳過有些解析代碼會把注釋行當(dāng)成數(shù)據(jù)導(dǎo)致錯(cuò)位。判斷計(jì)算是否正確的一個(gè)土辦法去NGA官網(wǎng)查幾個(gè)已知點(diǎn)的高程異常參考值比如0°, 0°附近、北京、紐約這些城市的值算一遍對比誤差在厘米級說明代碼基本沒問題。5.2 在極區(qū)或高緯度地區(qū)計(jì)算異常球諧函數(shù)在極區(qū)緯度接近±90°容易出現(xiàn)數(shù)值不穩(wěn)定。sin(φ)接近±1時(shí)P?nm的遞推公式會放大舍入誤差。遇到高緯度任務(wù)建議用遞推關(guān)系替代直接調(diào)用scipy的lpmv或者使用穩(wěn)定化的遞推公式比如Colombo和Sona提出的方法。另外EGM96發(fā)布的系數(shù)本身在南北緯88°以上有較大的外推誤差因?yàn)樾l(wèi)星軌道覆蓋不到極區(qū)。所以極區(qū)個(gè)別點(diǎn)算出來的值可信度要打折扣這一點(diǎn)要在成果報(bào)告中注明。5.3 高程異常轉(zhuǎn)換的精度驗(yàn)證方法算出來的高程異常能不能用最終要拿實(shí)測水準(zhǔn)點(diǎn)去驗(yàn)證。方法是在測區(qū)選擇若干已知正常高的水準(zhǔn)點(diǎn)用GNSS測出大地高H兩者相減得到“實(shí)測高程異常”再和EGM96計(jì)算的模型值對比。統(tǒng)計(jì)兩者差值的均值和標(biāo)準(zhǔn)差平原地區(qū)如果標(biāo)準(zhǔn)差小于0.3米可以直接用模型值做粗轉(zhuǎn)換如果要求厘米級精度就需要利用這些已知點(diǎn)做曲面擬合比如多項(xiàng)式擬合或克里金插值求出高程異常殘差的改正模型再疊加到EGM96結(jié)果上。我做過的項(xiàng)目里用5個(gè)均勻分布的已知水準(zhǔn)點(diǎn)做二次多項(xiàng)式擬合后殘差從0.5米壓到了5厘米以內(nèi)效果立竿見影。這個(gè)思路非常實(shí)用EGM96解決“大的架子”局部擬合解決“小的偏差”。5.4 重力異常的火山區(qū)畸變問題重力異常對地下質(zhì)量分布非常敏感在火山區(qū)域、大型礦體上方局部重力異?？梢赃_(dá)到數(shù)百毫伽的變化。EGM96受限于空間分辨率無法刻畫這種局部高頻信號所以在這些區(qū)域算出來的重力異常只能反映區(qū)域背景場不能用于局部資源勘探解釋。如果研究區(qū)域是這種強(qiáng)異常區(qū)建議疊加地面實(shí)測重力數(shù)據(jù)或者使用EGM2008的高階模型2190階來逼近局部場。我見過有同行直接用EGM96的格網(wǎng)值畫礦體異常圖結(jié)果解釋出來的“異常體”位置偏移了好幾個(gè)公里就是因?yàn)槟Ｐ偷牡皖l特性掩蓋了局部信號。6. 實(shí)操總結(jié)與個(gè)人經(jīng)驗(yàn)EGM96作為一個(gè)發(fā)布快三十年的模型在今天依然有它的生命力尤其是在工程高程轉(zhuǎn)換的效率和解算便捷性上依然很能打。但用這個(gè)模型心里要有桿秤它的低頻成分可靠高頻成分受限不同區(qū)域的誤差表現(xiàn)差異很大。凡是拿它出成果之前一定要用實(shí)測數(shù)據(jù)驗(yàn)證別偷懶。我在實(shí)際項(xiàng)目中摸索出的一個(gè)流程是先用EGM96快速算出測區(qū)的高程異常背景場然后選6到10個(gè)均勻分布的已知水準(zhǔn)點(diǎn)做殘差擬合最后用擬合模型修正整個(gè)測區(qū)的轉(zhuǎn)換結(jié)果。這樣既保證了效率又把精度控制在了厘米級。如果手里有歷史項(xiàng)目的EGM96計(jì)算結(jié)果也可以像“經(jīng)驗(yàn)?zāi)０濉币粯酉葏⒖贾浪阈马?xiàng)目的誤差量級心里先有個(gè)底。有個(gè)小技巧想分享計(jì)算點(diǎn)比較多的時(shí)候別頻繁調(diào)用數(shù)學(xué)庫函數(shù)。把常用階次的勒讓德函數(shù)值先算好緩存起來因?yàn)橥痪暥壬辖?jīng)度變化時(shí)勒讓德部分完全不變變的只是cos(mλ)和sin(mλ)項(xiàng)。這樣優(yōu)化后計(jì)算速度能提升好幾倍。測區(qū)大、點(diǎn)數(shù)多的時(shí)候這個(gè)優(yōu)化是實(shí)打?qū)嵉氖找?。EGM96這套計(jì)算流程說難不難說簡單也不簡單。把原理吃透、代碼寫對、驗(yàn)證做扎實(shí)高程異常和重力異常的計(jì)算其實(shí)是很順手的事。希望這篇分享能幫你少走點(diǎn)彎路。本文還有配套的精品資源點(diǎn)擊獲取