地球重力場模型計算點向徑取值對模型重力場元影響

邢志斌，李姍姍，范 雕，張 馳

（信息工程大學 地理空間信息學院，鄭州 450001）

地球重力場信息在遠程武器的精確打擊，航天器的精密定軌、慣性導航系統的精確導航、地球物理現象的研究、資源勘探以及地球環境問題研究等方面起著重要作用[1-5]。地球重力場模型可為以上科學研究、工程實踐等提供必要的數據支持。鑒于地球重力場模型的重要意義，世界主要國家無不投入巨大精力研究如何快速精確獲取地球重力場信息用于構建地球重力場模型。時至今日，重力場模型已由最初的低分辨率、低階次發展為目前的高分辨率、超高階次，重力場信息的獲取手段也由原來的陸地重力測量發展為陸地、海上、空中和太空[6]，數據的類型也由最初的單一數據發展為多源數據的融合。在利用超高階次地球重力場模型計算重力場元時，模型重力場元計算精度受計算點高程的影響，在精度要求較低時，為提高計算效率，通常不考慮高程信息影響，即在球近似下計算，此時，位系數模型中的R r=1[7]。然而，更多的時候對精度要求較高，此時必須考慮高程項的影響，許多文獻在計算模型重力場元時為提高計算效率通常會將分母中的計算點向徑近似為r=R+H[8-10（]事實上美國國家地理空間情報局（NGA）提供的程序中采用的是計算點真實向徑[12]），顧及計算精度與效率，可將計算點高程在球面做泰勒級數展開至一次項[8-9]、二次項或更高次項[10]。由于比例因子R取值地球長半軸[13（]也有文獻采用地球平均半徑[14]），單純的將計算點向徑近似為r=R+H，將會帶來千米級距離誤差，同時由于計算區域格網點坐標一般是以大地坐標 (B,L,H)表示，而位系數模型是以球坐標 (r,φ,λ)表示[13]，計算時需要將大地坐標轉化為球坐標，在轉換時同一緯圈B不同的高度H對應不同的φ，快速計算的前提是同一緯圈對應的φ相等，這些近似誤差將不可避免的引入到地球重力場元的計算中。本文針對計算點向徑取值展開研究，利用 EGM2008地球重力場模型以模型垂線偏差計算為例分析了r的不同取值造成的距離誤差及對地心緯度φ的影響，并與實測垂線偏差進行比對，為顧及計算精度與效率，提出用等效高程代替實際高程進行泰勒級數展開快速計算模型垂線偏差，并分析了其適用性。

1 模型垂線偏差計算

1.1 模型垂線偏差計算公式

計算子午方向和卯酉方向模型垂線偏差的嚴密公式為[10-11]：

式中，R通常為地球長半軸，表示第 ,( )i j個計算點的地心球坐標，為完全正常化的球諧系數，為完全正常化的締合勒讓德函數，n表示階，m表示次。

通常我們會使rij=R+Hij，Hij表示計算點的高程，在僅考慮計算效率的情況下，忽略高程項的影響，可直接采用球近似下垂線偏差的計算公式[7]：

顯然式(2)忽略地形影響會降低計算精度，為此一般將式(2)在球面對計算點高程做泰勒級數展開至一階項得[10,15]垂線偏差子午分量：

和卯酉分量：

為進一步提高精度，可將式(2)按泰勒級數展開至二階項、三階項、四階項，具體展開式及計算過程見文獻[10]。

1.2 大地坐標與地心球坐標轉換

計算重力場元時，通常給出大地坐標(B,L,H)，而式(1)是在地心球坐標系(r,φ,λ)下表達，因此需要將 (B,L,H)轉換到(r,φ,λ)下，再利用式(1)計算模型垂線偏差。首先由式(5)將大地坐標系 (B,L,H)轉換到地心空間直角坐標系 (X,Y,Z)下：

三個坐標系之間的關系如圖1所示。

圖1 ( B , L, H)、( X , Y, Z)和 ( r ,φ, λ)之間關系Fig.1 Relationship among(B , L, H),(X , Y, Z)and (r ,φ, λ)

1.3 模型垂線偏差計算公式分析

1.3.1 向徑誤差分析

采用式(1)計算模型垂線偏差時，由式(6)確定每個計算點的向徑，將會耗費大量的時間，通常會令ri=R+Hi，由于R的值為地球長半軸，因此除赤道區域以外，在其他區域采用這種近似將會帶來很大的距離誤差。如果令ri=R+Hi，R=6371000.0m，表示地球平均半徑，這樣會在中緯度地區改善距離精度，但在赤道和兩極地區則又會失真。此時，為提高計算效率和距離精度，可計算出同一緯圈橢球面向徑r0再加高程Hi，即r=r0+Hi。

為比較以上三種近似帶來的誤差影響，分別令：

采用5′×5′全球陸地和海底地形格網數據（由ETOPO1數據平滑后得到，如圖2所示）計算每個點的向徑，并與由式(6)計算的向徑值進行比較，統計結果如表1、圖2～5所示。

表1 不同近似值與真值之差統計量（單位：m）Tab.1 Differences between different approximations and true values (unit: m)

圖2 全球陸地與海底地形Fig.2 Global land and seabed topography

圖3 參考橢球面向徑和地面高程之和與真值之差 ( vr0 +H- r)Fig.3 Differences between the summation of reference ellipsoid and elevation and true values (vr0 +H- r)

圖4 地球長半軸和地面高程之和與真值之差(vR+H-r)Fig.4 Differences between the summation of semi-major axis and elevation and true values (vR+H-r)

圖5 地球平均半徑和地面高程與真值之差(v )Fig.5 Differences between the summation of Earth' mean radius and elevation and true values (v )

由以表1和圖2～5可以看出，采用不同的近似值，距離誤差相差很大。在三個近似值中，近似值（1）的精度最高，與真值相比只有厘米級的誤差，近似值（2）、（3）與真值相比，誤差達到了km級，并且這些誤差具有系統性，即同一緯圈附近的誤差基本一致，近似值（1）在赤道地區誤差最小，向兩極逐漸增大，而近似值（3）造成的誤差在南北緯30°附近最小，兩極地區最大。

1.3.2 不同向徑對球心緯度影響分析

由式(5)(6)可知，地心緯度φ與高程H有關，因此，同一緯圈B如果格網點的高程不同，則φ不同。但是在模型重力場元計算中，我們通常認為同一緯圈對應的φ相同，為分析這種近似誤差，分別計算格網點真實高度（海深）處的地心緯度φ以及對應的參考橢球面上的地心緯度φ0，統計結果如表2、圖6所示。

表2 橢球面地心緯度與計算點地心緯度之差統計量（單位：″）Tab.2 Differences between ellipsoid and calculation points' geocentric latitude (unit: arc second)

圖6 參考橢球面地心緯度與實際點地心緯度之差 ( vφ0-φ)Fig.6 Differences between references ellipsoid surfaces'geocentric latitude and real geocentric latitude (vφ0-φ)

由表2、圖6可知，由地球上實際高程（海深）造成的地心緯度誤差很小，絕對值最大值也就 1″左右，這對于最高分辨率5′×5′的地球重力場模型來說是可以忽略不計的。由圖6可以看出，地心緯度誤差具有一定的系統性，在北半球海洋區域多為正值，而南半球海洋區域多為負值，并且相比于參考橢球面地心緯度φ0，高程為正會使得地心緯度有增大的趨勢。

1.4 模型垂線偏差快速計算

向徑不同取值之間存在 km級的差距，采用式(3)(4)快速計算模型垂線偏差將會產生誤差，但是逐點計算又會耗費大量時間。此時可采用等效高程Hi′j代替式(3)(4)中Hij，即垂線偏差子午分量：和卯酉分量：

為提高計算效率，r0或r可提前計算出，可采用向量運算的方式一次性計算出計算區域的等效高程，向量運算的計算方法可參考文獻[10][17]。

為分析模型垂線偏差快速計算方法的適用性，選擇我國青藏高原一 6°×6°地勢變化劇烈區域，采用式(7)(8)，由EGM2008地球重力場模型計算分辨率5′×5′的子午方向和卯酉方向格網垂線偏差，與(1)式計算的模型垂線偏差作對比。計算區域高程如表3、圖7所示。

表3 研究區域高程統計（單位：m）Tab.3 Statistics of study area elevations (unit: m)

圖7 計算區域高程Fig.7 Elevations of study areas

根據表3、圖7可知，該區域地勢變化相當復雜，平均高度接近4000 m，但最低高度只有27 m，而最高高度接近7000 m。將前50行格網分成一部分，后70行格網分成一部分，兩部分高程變化差距很大，前50行格網的高程標準差只有百米量級，而后70行格網則高達近2 km。展開至不同階次泰勒級數計算模型垂線偏差與嚴密公式計算結果殘差如圖8～11所示（左圖為子午方向、右圖為卯酉方向），整體統計如表4所示，分區統計如表5所示。

由表4可知，就整個區域而言，展開階次對精度沒有實質性的影響，并且子午方向的精度一直低于卯酉方向的精度。對模型垂線偏差計算精度進行分區統計時（如表5所示），前50行的精度遠遠優于后70行的精度，并且展開至二階項以后其精度趨于穩定。

后 70行計算結果與文獻[10]給出的結果有一定的出入，分析原因主要有：

1）文獻[10]中的計算點向徑取值r=R+H要大于計算點實際向徑km級，這相當于計算的是高空的重力場元在一定程度上平滑了重力場信號，重力場元與高程可近似為線性關系，因此泰勒級數展開至不同階次，對計算精度提升較大；

表4 展開至不同階次垂線偏差整體精度統計（單位：″）Tab.4 Accuracy statistic of extended to different orders for all study areas (unit: arcsec)

表5 展開至不同階次垂線偏差分區精度統計（單位：″）Tab.5 Accuracy statistic of extended to different orders for different areas (unit: arcsec)

圖8 垂線偏差展開至一階項與模型真值之差Fig.8 Differences between true values and extended to first-order terms

圖9 垂線偏差展開至二階項與模型真值之差Fig.9 Differences between true values and extended to second-order terms

圖10 垂線偏差展開至三階項與模型真值之差Fig.10 Differences between true values and extended to third-order terms

圖11 垂線偏差展開至四階項與模型真值之差Fig.11 Differences between true values and extended to forth-order terms

2）文獻[10]中計算區域的地勢復雜程度遠低于圖7區域，重力場變化的劇烈程度也低于圖7區域；

3）式(7)(8)泰勒級數展開至一階項，本質上是利用徑向梯度延拓距離H′，一方面重力場模型中很難含有地表的梯度信息，另一方面H′過大會造成延拓發散。

由此可知，利用式(7)(8)快速計算模型垂線偏差并不完全適用于任意地區，在地勢變化復雜的區域需要綜合考慮各方面的因素。

2 實驗對比分析

2.1 近似值與模型垂線偏差真值比較

為分析計算點向徑誤差對模型垂線偏差的影響，選擇圖7所示區域，分別將向徑設為：1）r=r0+Hi；2）

計算分辨率 5′×5′的子午方向和卯酉方向格網垂線偏差，與采用真實向徑計算的模型垂線偏差作對比。統計結果如表6、圖12～19所示。

根據圖12、圖13可知，該地區垂線偏差變化比較劇烈，最大值可達1′多。由表6可知，當r=r0+Hi時計算的模型垂線偏差與真值之差非常小，RMS不會超過0.01″（圖14、圖15所示），遠低于當前一等天文垂線偏差的測量誤差（0.3″），這說明在計算模型垂線偏差時這種近似是可行的。

表6 垂線偏差不同近似值與模型真值之差統計（單位：″）Tab.6 Differences between different approximations and model true values (unit: arcsec)

圖12 子午方向模型垂線偏差真值Fig.12 Model vertical deflection true values in meridian direction

圖13 卯酉方向模型垂線偏差真值Fig.13 Model vertical deflection true values in prime direction

圖14 子午方向模型垂線偏差殘差1Fig.14 Residual 1 of model vertical deflection in meridian direction

圖15 卯酉方向模型垂線偏差殘差1Fig.15 Residual 1 of model vertical deflection in prime direction

圖16 子午方向模型垂線偏差殘差2Fig.16 Residual 2 of model vertical deflection in meridian direction

圖17 卯酉方向模型垂線偏差殘差2Fig.17 Residual 2 of model vertical deflection in prime direction

圖18 子午方向模型垂線偏差殘差3Fig.18 Residual 3 of model vertical deflection in meridian direction

圖19 卯酉方向模型垂線偏差殘差3Fig.19 Residual 3 of model vertical deflection in prime direction

2.2 模型值與實測值比

為比較計算點向徑r的不同取值對模型垂線偏差計算結果的影響，分別計算了25個實測天文大地垂線偏差點在四種取值情況下的模型垂線偏差，與其真值比較，殘差統計結果如表7、圖20、圖21所示。

由表7、圖20、圖21可知，在子午方向上，采用計算點向徑真值與近似值 1（r=r0+Hi）計算的模型垂線偏差與實測數據相比，精度是一致的，相比于其他兩種近似（ri=R+Hi和垂線偏差計算精度有了很大的提高，并且這樣的精度與EGM2008模型在美國本土的精度是基本一致的。但是在卯酉方向，采用向徑真值和近似值 1計算垂線偏差精度卻有了略微的降低，分析原因可能是受我國西高東低地勢的影響，由于僅采用了25個實測點進行驗證，因此造成這種現象的原因還需做進一步研究。

圖20 子午方向垂線偏差殘差（單位：″）Fig.20 Residual of vertical deflection in meridian direction (unit：arc second)

表7 垂線偏差模型值與實測值之差統計（單位：″）Tab.7 Differences between model values and measured values (unit: arc second)

圖21 卯酉方向垂線偏差殘差（單位：″）Fig.21 Residual of vertical deflection in prime direction (unit：arcsec)

3 結 論

本文以模型垂線偏差計算為例，分析了計算點向徑不同取值造成的距離誤差及對地心緯度的影響，進而通過實驗分析了距離誤差對模型垂線偏差計算精度的影響。實驗表明：1）r=r0+Hi距離誤差造成的影響可忽略不計；2）r=R+Hi在赤道地區的影響很小，向兩極逐漸增大；3）r=+Hi在南北緯30°地區影響最小；4）實際計算點的地心緯度與計算點對應的參考橢球面上點的地心緯度之差最大值在 1″左右，對垂線偏差造成的影響絕對值最大不超過0.07″，這些影響對模型重力場元計算可忽略不計，在格網模型重力場元計算時同一緯圈可只計算一次締合勒讓德函數，為快速計算帶來了極大的方便。與25個實測點垂線偏差進行比較，r取真值與r=r0+Hi計算的子午方向垂線偏差結果一致，且精度最高，達到了EGM2008重力場模型在美國本土的精度，但是卯酉方向垂線偏差計算的精度卻略低于r=R+Hi計算的精度，由于用來比較的實測數據太少，因此還很難分析出精度變差的原因。用計算點等效高程Hi′j代替實際高程Hij，可在原有程序的基礎上利用式(7)(8)快速計算格網模型垂線偏差，但是其適用范圍是有限的，并不適用于地勢變化劇烈地區，如若進一步提高計算精度可采用在橢球面展開的方法[15]。