我國陸海統一似大地水準面構建的三維重力矢量法

邢志斌，李姍姍

信息工程大學，河南 鄭州 450000

目前確定(似)大地水準面的方法主要有幾何法、重力法以及組合法[1]。早期幾何法主要以天文重力水準的方式確定似大地水準面[2]，現在多以擬合GPS/水準點為主，雖然可以獲得很高的精度，但僅限于小范圍、GPS/水準點豐富的地區[3]，其計算結果為區域(似)大地水準面。幾何法的缺陷在于需要進行大量的水準測量，耗費巨大的人力、物力，海洋測高雖然能以很高的精度確定大地水準面，但不適用于陸地及陸海交界處[4]。重力法主要是解算物理大地測量邊值問題，再由Bruns公式轉換為大地水準面高或高程異常[5]，計算結果為絕對大地水準面，可統一全球高程基準[6]，適用范圍廣，但計算復雜、計算量大，需要與本國高程基準擬合。組合法則是將幾何法與重力法相結合確定(似)大地水準面，其本質是將重力似大地水準面數值擬合到GPS/水準點上，與本國高程基準相統一[7-8]，結合了GPS/水準精度高、重力(似)大地水準面分辨率高的優點，但其擬合過程并不嚴密。世界主要國家和地區如美國、歐洲、日本、加拿大及我國所建立的高精度(似)大地水準面模型多以組合法建立[4,9-17]。由于重力場水平分量-垂線偏差對地形高頻信息敏感，因此首先利用三維重力矢量-格網垂線偏差與格網重力異常，聯合格網高程數據求得格網中點間相對高程異常(高程異常差)，然后通過GPS/水準點的控制，構成緊密的幾何條件，進行嚴密的平差，從而可以獲得高分辨率、高精度似大地水準面的數值模型[18-19]。該方法的優點是充分吸收了高頻的地形信息，彌補了我國重力場數據的不足；同時以確定相對高程異常代替以往的直接確定高程異常的方法，而且這種格網相對高程異常構成了嚴密的幾何條件，可以通過平差消除矛盾、提高精度，比以往采用數值擬合的簡單方法提高了嚴密度。論文基于此，展開深入研究，旨在提高我國陸海統一的似大地水準面模型的建模精度。

1 格網高程異常差及其平差模型的建立

1.1 格網垂線偏差計算高程異常差的基本原理

在計算高程異常差時所用的垂線偏差數據是只包含子午方向和卯酉方向分量的格網數據，根據天文重力水準方法計算似大地水準面的原理可得高程異常差與垂線偏差的關系為[20]

(1)

以任一格網為例說明利用垂線偏差計算相鄰格網中點高程異常差的方法，如圖1所示。

圖1中，黑色三角形表示格網中點?？紤]重力異常改正項[21]，子午方向和卯酉方向相鄰格網中點之間的高程異常差分別為

(2)

(3)

圖1 垂線偏差與高程異常差的關系Fig.1 Using deflection of the vertical to calculate height anomaly difference

1.2 計算高程異常差的平差模型

(4)

式(4)可寫成矩陣形式為

(5)

相應的法方程為

(6)

法方程的詳細構建與解算過程見文獻[20]。經解算后得

(7)

2 基于三維重力矢量的我國似大地水準面模型的建立

2.1 數據準備

為建立我國1′×1′陸海統一的似大地水準面模型，本文采用以下數據：分辨率為1′×1′的格網垂線偏差、格網重力異常、高程數據以及全國6600多個GPS/水準數據。

垂線偏差數據主要是采用“移去-恢復”技術[22]，根據邊值理論由重力、地形以及地球重力場模型確定[23-25]。垂線偏差精度優于1.5″，格網點間距以1.8 km計算，由此產生的誤差低于2 cm，加入GPS/水準控制以后，可滿足建立厘米級似大地水準面的要求。我國陸地子午與卯酉垂線偏差分量如圖2、圖3所示。重力異常數據主要利用我國及周邊81萬多個離散重力異常數據經格網化所得，如圖4所示；地形數據由分辨率3″×3″的SRTM數據取平均后處理成分辨率1′×1′的地形數據，如圖5所示。共收集包括國家GPS A、B級點、工程點，以及部分地區測量的C級網點在內的陸地GPS/水準點6700多個，如圖6所示，剔除粗差點后實際上采用了6600多個計算并檢核了似大地水準面模型。

2.2 分區似大地水準面的確定

依據上述數據的質量以及分布情況，將全國分為6個區域，分區確定區域似大地水準面模型。區域1包括華東、華中、華南地區；區域2為云貴地區；區域3是藏南地區；區域4包括新疆南部及西藏北部地區；新疆北部地區作為區域5；區域6為東北地區。不同區域范圍以及控制點、檢核點分布如圖7—圖12所示，分區似大地水準面模型如圖13所示，檢核點的殘差統計如表1、圖14所示。

表1 殘差統計特征

圖7反映出區域1 GPS/水準點較多，尤其是東部沿海地區，但整體上分布不均，我國臺灣地區沒有任何數據，而中南部地區GPS/水準點相對匱乏，但由GPS/水準點檢核后發現該地區似大地水準面的精度并沒有受此影響。該區域似大地水準面的精度達到了3.9 cm。

由圖8可知，區域2 GPS/水準點較少，分布極其不均，圖5也顯示該區域地勢復雜，因此選擇大部分的點作為控制點。邊界區域都布設了控制點，個別地區控制點分布集中、距離較近，這是因為在不斷選點計算時發現，這幾個點單獨做檢核點時殘差較大，但做控制點時，其對其他檢核點的殘差并沒有大的影響，因此排除了這些點是粗差點的可能。值得注意的是沿24°N緯線分布的3個GPS/水準點附近沒有任何控制點，但并沒有影響似大地水準面的精度，由圖4也可以看出該地區重力場的變化并不明顯。

圖9顯示區域3是我國地勢最高的區域，GPS/水準點的分布并不均勻，為保證精度與可靠性著重在邊界地區、橫斷山脈附近布設較多控制點。但仍有部分地區缺少GPS/水準點，這些地區的似大地水準面模型既沒有得到合理控制也沒有得到有效檢核，幸運的是這些地區重力異常的變化并不像喜馬拉雅山脈地區變化的那么劇烈，為保證可靠性對北部山脈地區施加大量的控制點。該區域似大地水準面模型經過GPS/水準點的檢核表明在GPS/水準點有效控制的區域精度達到了5 cm。

由圖10可知，區域4 GPS/水準點多集中分布在東部地區，在該區域東南角處，控制點選取較為密集，主要是因為此處是階梯的交界處，也正是由于如此密集的選點，使得該地區的似大地水準面的精度得到了保證。經GPS/水準點檢核區域4的似大地水準面模型精度達到了4 cm。

由圖11可知，區域5的GPS/水準點較少，分布多集中在中部地區，由圖4可知中部地區的重力異常變化較為劇烈，因此適度增加控制點有利于增加該區域似大地水準面的可靠性。經GPS/水準點檢核區域5的似大地水準面模型精度優于6 cm。

圖12表明區域6 GPS/水準點主要集中在中西部地區，東部地區較為稀少，且重力場變化較為劇烈，因此該區域需要較多的控制點，以保證似大地水準面的可靠性，在中西部地區，控制點多為三角網的角點處。該區域似大地水準面模型的整體精度達到了4 cm，并且大部分地區選用了較少的控制點，但是在東部地區受重力場變化、地形等因素的影響，絕大部分的GPS/水準點都選作了控制點。

2.3 全國陸海統一似大地水準面確定

綜合上述陸地控制點，并聯合海島礁聯測GPS/水準點進行整體平差計算，建立了我國陸海統一的似大地水準面模型，如圖15所示。其中為消除不同區域控制點對似大地水準面整體精度的影響，適當調整了部分的檢核點作為控制點，最終采用共計2444個GPS/水準點做控制，經4241個GPS/水準點檢核，結果表明我國似大地水準面的整體精度達到4.01 cm。檢核點殘差統計如圖16、圖17、表2所示。

圖16 殘差分布統計Fig.16 The distribution of residual

最大值殘差特征值/m殘差范圍/個最小值平均值RMS

2.4 相對精度評定

絕對精度檢核并未考慮GPS/水準點的誤差累積，因此，整體上4 cm的精度只能說明所建立的似大地水準面與GPS/水準點的“吻合”情況較好。在計算過程中，地形復雜、重力場變化劇烈的地區，多是增加了GPS/水準點，增強對該區域似大地水準面的控制，考慮到絕對精度表示的是似大地水準面與本國GPS/水準點擬合以后的精度，在這種強制擬合中忽視了GPS/水準點中水準測量誤差累積的影響，為此評定了似大地水準面模型的相對精度，通過做差消除了部分水準測量累積誤差影響。

2.4.1 相對精度評定方法

相對精度評定的具體步驟如下：

(1) 任選一檢核點A，以A點為中心選擇一定半徑范圍的檢核點Bi(i=1,2,…,n，n為以A點為中心一定半徑范圍內的檢核點的個數)。

(3)A點的高程異常真值ζA與Bi點的高程異常真值ζi做差，得

vi=ζA-ζi

(4) 步驟(2)與步驟(3)計算所得的殘差之間做差，稱為相對高程異常差，得

dvi=v′-vi

對dvi做誤差統計，得到似大地水準面模型的相對精度。

2.4.2 似大地水準面相對精度分析

選擇檢核點0.5°半徑范圍的其他檢核點，按照上述方法共統計了49 477組相對高程異常差，統計結果如表3、圖18、圖19所示。不同統計半徑下的分區統計如表4所示，分區精度統計結果如圖20所示。

圖18 相對高程異常差分布統計Fig.18 Distribution of differences of relative height anomaly

最大值殘差特征值/m殘差范圍/個最小值平均值RMS

表4 不同統計半徑下分區相對精度統計

通過相對精度的評定，可以發現東部地區似大地水準面精度優于西部地區精度，這主要與西部地區缺少必要的觀測數據有關。值得注意的是統計半徑0.5°時，區域2的相對精度最高，達到了2 cm，但由于該地區GPS/水準點的數量最少，并且多數都作為控制點，造成了參與檢核的相對高程異常較少，因此進一步擴大檢核時選點的半徑至1°、2°。由表4可知隨著選點半徑的增大，區域2的精度不斷降低，其他4個區域的精度基本穩定，因此可認為這4個區域先前的精度評定是可信的?？紤]到區域2 GPS/水準點嚴重不足，可以認為統計半徑為0.5°情況下得出的2 cm精度并不可信，當選點半徑擴大至2°已包含了區域2絕大部分的檢核點，因此認為在該選點半徑情況下得出的6.7 cm的精度較為可靠。

圖2 我國陸地1′×1′分辨率子午垂線偏差Fig.2 China land grid vertical deflections of meridian direction with a resolution of 1′×1′

圖3 我國陸地1′×1′分辨率卯酉垂線偏差Fig.3 China land grid vertical deflections of prime direction with a resolution of 1′×1′

圖4 我國陸地1′×1′分辨率重力異常Fig.4 China land grid gravity anomaly with a resolution of 1′×1′

圖5 我國陸地1′×1′SRTM地形數據Fig.5 China land SRTM terrain data with a resolution of 1′×1′

圖6 我國陸地GPS/水準點分布Fig.6 Distribution of China land GPS/leveling points

圖7 區域1的控制點、檢核點分布Fig.7 Distribution of control and check points of area 1

圖8 區域2的控制點、檢核點分布Fig.8 Distribution of control and check points of area 2

圖9 區域3的控制點、檢核點分布Fig.9 Distribution of control and check points of area 3

圖10 區域4的控制點、檢核點分布Fig.10 Distribution of control and check points of area 4

圖11 區域5的控制點、檢核點分布Fig.11 Distribution of control and check points of area 5

圖12 區域6的控制點、檢核點分布Fig.12 Distribution of control and check points of area 6

圖13 陸地分區似大地水準面模型及精度Fig.13 The quasi-geoid model precision of different areas

圖15 陸海統一的似大地水準面模型Fig.15 China land and ocean quasi-geoid model

圖17 殘差大小統計Fig.17 Statistic of residual

圖19 全國陸地相對精度分區統計Fig.19 Statistic of precision of different areas

圖20 相對高程異常差大小統計Fig.20 The statistic of differences of relative height anomaly

3 結 論

目前似大地水準面建模的通用方法是按照Molodensky理論直接計算一點的高程異常，再用GPS/水準擬合系統差。由于大地水準面不能有效地吸收地形等高頻信息并且難以利用GPS/水準數據對大地水準面做真正地實際控制，因此很難再大幅度提高似大地水準面精度。本文基于重力場水平分量-垂線偏差對地形信息敏感的特點，在根據邊值理論通過重力與地形數據確定格網垂線偏差模型的基礎上，利用三維重力矢量-格網垂線偏差與格網重力異常，并聯合格網高程數據求得格網中點間相對高程異常，通過GPS/水準點的控制，構成緊密的幾何條件，進行嚴密平差，從而計算得到了我國分辨率為1′×1′的陸海統一的似大地水準面模型，并進行了整體與分區的絕對精度與相對精度的評定。檢核結果表明采用本文方法構建的我國似大地水準面模型的整體絕對精度達到了4 cm，相對精度優于7 cm。本文提出的基于三維重力矢量構建我國大地水準面的方法，充分利用了垂線偏差對地形信息敏感的特點，以高分辨率地形數據彌補重力場數據的不足，以確定相對高程異常代替以往的直接確定高程異常的方法，類似于GPS從直接定位到差分定位的轉變，而且這種格網相對高程異常構成了嚴密的幾何條件，可以通過平差消除矛盾提高精度，比以往采用數值擬合的簡單方法提高了嚴密度，為我國似大地水準面模型的建立提供了一種新的思路與方法。

參考文獻：

[1] 蔣光偉, 郭春喜, 田曉靜, 等. 基于多種方法構建似大地水準面模型推估特性分析[J]. 地球物理學進展, 2014, 29(1): 51-56. DOI: 10.6038/pg20140107.

JIANG Guangwei, GUO Chunxi, TIAN Xiaojing, et al. The Analysis of Extrapolated Capability of Quasi-geoid Based on Muti-methods[J]. Progress in Geophysics, 2014, 29(1): 51-56. DOI: 10.6038/pg20140107.

[2] HOFMANN-WELLENHOF B, MORITZ H. Physical Geodesy[M]. New York: SpringerWienNewYork, 2005.

[3] KAVZOGLU T, SALA M H. Modelling Local GPS/Levelling Geoid Undulations Using Artificial Neural Networks[J]. Journal of Geodesy, 2005, 78(9): 520-527.

[4] 李建成, 陳俊勇, 寧津生, 等. 地球重力場逼近理論與中國2000似大地水準面的確定[M]. 武漢: 武漢大學出版社, 2003.

LI Jiancheng, CHEN Junyong, NING Jinsheng, et al. The Theory of Earth’s Gravity Field Approximation and Determination of the China Quasi-geoid 2000[M]. Wuhan: Wuhan University Press, 2003.

[6] GERLACH C, RUMMEL R. Global Height System Unification with GOCE: a Simulation Study on the Indirect Bias Term in the GBVP Approach[J]. Journal of Geodesy, 2013, 87(1): 57-67.

[7] JIANG Ziheng, DUQUENNE H. On the Combined Adjustment of a Gravimetrically Determined Geoid and GPS Levelling Stations[J]. Journal of Geodesy, 1996, 70(8): 505-514.

[8] KOTSAKIS C, SIDERIS M G. On the Adjustment of Combined GPS/Levelling/Geoid Networks[J]. Journal of Geodesy, 1999, 73(8): 412-421.

[9] SMITH D A, ROMAN D R. GEOID99 and G99SSS: 1-Arc-Minute Geoid Models for the United States[J]. Journal of Geodesy, 2001, 75(9-10): 469-490.

[10] ROMAN D R, WANG Yanming, HENNING W, et al. Assessment of the New National Geoid Height Model,GEOID03[J]. Surveying and Land Information Science: Journal of the American Congress on Surveying and Mapping, 2004, 64(3): 153-162.

[11] WANG Y M, SALEH J, LI X, et al. The US Gravimetric Geoid of 2009 (USGG2009): Model Development and Evaluation[J]. Journal of Geodesy, 2012, 86(3): 165-180.

[12] DENKER H. The European Gravimetric Quasigeoid EGG2008[C]∥American Geophysical Union, Spring Meeting 2009. New York: AGU, 2009.

[13] ODERA P A, FUKUDA Y, KUROISHI Y. A High-resolution Gravimetric Geoid Model for Japan from EGM2008 and Local Gravity Data[J]. Earth, Planets and Space, 2012, 64(5): 361-368.

[14] HUANG Jianliang, VéRONNEAU M. Canadian Gravimetric Geoid Model 2010[J]. Journal of Geodesy, 2013, 87(8): 771-790.

[15] http:∥www.nrcan.gc.ca/earth-sciences/geomatics/geodetic-reference-systems/9054.

[16] 李建成. 最新中國陸地數字高程基準模型: 重力似大地水準面CNGG2011[J]. 測繪學報, 2012, 41(5): 651-660, 669.

LI Jiancheng. The Recent Chinese Terrestrial Digital Height Datum Model: Gravimetric Quasi-geoid CNGG2011[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(5): 651-660, 669.

[17] 束蟬方, 李斐, 李明峰, 等. 應用Bjerhammar方法確定GPS重力似大地水準面[J]. 地球物理學報, 2011, 54(10): 2503-2509.

SHU Chanfang, LI Fei, LI Mingfeng, et al. Determination of GPS/gravity Quasi-geoid Using the Bjerhammar Method[J]. Chinese Journal of Geophysics, 2011, 54(10): 2503-2509.

[18] 朱雷鳴. 基于相對高程異常差平差法的區域似大地水準面精化[D]. 鄭州:信息工程大學, 2011.

ZHU Leiming. Refining of Quasi-geoid Based on the Relative Height Anomaly Difference Adjustment[D]. Zhengzhou:Information Engineering University, 2011.

[19] 孫文. 我國高分辨率三維重力場與似大地水準面研究[D]. 鄭州:信息工程大學, 2014.

SUN Wen. Determination of High Resolution of 3D Gravity Field and Quasi-Geoid of China[D]. Zhengzhou:Information Engineering University, 2014.

[20] 邢志斌, 李姍姍, 王偉, 等. 利用垂線偏差計算高程異常差法方程的快速構建方法[J]. 武漢大學學報(信息科學版), 2016, 41(6): 778-783.

XING Zhibin, LI Shanshan, WANG Wei, et al. Fast Approach to Constructing Normal Equation during the Time of Calculating Height Anomaly Difference by Using Vertical Deflections[J]. Geomatics and Information Science of Wuhan University, 2016, 41(6): 778-783.

[21] 吳曉平. 似大地水準面的定義及在空中測量中涉及的問題[J]. 測繪科學, 2006, 31(6): 24-25.

WU Xiaoping. Title Definition of Quasi-geoid and Some Questions Encountered in Airborne Gravimetry[J]. Science of Surveying and Mapping, 2006, 31(6): 24-25.

[22] 寧津生, 郭春喜, 王斌, 等. 我國陸地垂線偏差精化計算[J]. 武漢大學學報(信息科學版), 2006, 31(12): 1035-1038.

NING Jinsheng, GUO Chunxi, WANG Bin, et al. Refined Determination of Vertical Deflection in China Mainland Area[J]. Geomatics and Information Science of Wuhan University, 2006, 31(12): 1035-1038.

[23] 范昊鵬, 李姍姍. 局部區域模型重力異?？焖偎惴ㄑ芯縖J]. 大地測量與地球動力學, 2013, 33(6): 28-30, 35.

FAN Haopeng, LI Shanshan. Study on a Fast Algorithm for Model Gravity Anomalies in Local Areas[J]. Journal of Geodesy and Geodynamics, 2013, 33(6): 28-30, 35.

[24] 李姍姍, 吳曉平, 張傳定, 等. 顧及地形與完全球面布格異常梯度項改正的區域似大地水準面精化[J]. 測繪學報, 2012, 41(4): 510-516.

LI Shanshan, WU Xiaoping, ZHANG Chuanding, et al. Regional Quasi-geoid Refining Considering Corrections of Terrain and Complete Spherical Bouguer Anomaly’s Gradient Term[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(4): 510-516.

[25] 邊少鋒, 薛芳俠. 論地形垂線偏差中央區貢獻的計算[J]. 測繪學報, 1997, 26(1): 33-36, 57.

BIAN Shaofeng, XUE Fangxia. Discussion on the Calculation of Plumb Line Deviation[J]. Acta Geodaetica et Cartographica Sinica, 1997, 26(1): 33-36, 57.