李德強

(臨沂市國土資源局測繪院，山東 臨沂 276000)

利用改進的Stokes-Helmert邊值理論解算高精度似大地水準面

李德強

(臨沂市國土資源局測繪院，山東 臨沂 276000)

利用改進的Stokes-Helmert邊值問題實現(xiàn)了臨沂市厘米級似大地水準面精化。首先，基于Stokes理論和Molodensky理論，聯(lián)合精密確定地表及其外部擾動引力位的嚴密解算理論，給出Stokes-Helmert邊值問題的數(shù)學(xué)描述，以及直接地形影響和間接地形影響的嚴密理論表達式；然后，利用多源觀測資料，根據(jù)“移去-恢復(fù)”技術(shù)構(gòu)建臨沂高精度的重力似大地水準面模型；最后，利用GPS/水準高程異常對重力似大地水準面模型進行控制擬合，求得最終的大地水準面模型，其外符合精度達到1.6 cm。

似大地水準面精化；Stokes-Helmert法；多面函數(shù)擬合；地形影響

大地水準面的確定目前仍基于求解經(jīng)典大地測量邊值問題，最具有代表性的是Stokes理論和Molodesnky理論。為保證大地水準面外部擾動位為諧函數(shù)，Stokes理論要求大地水準面外部無質(zhì)量。將大地水準面外部的地形質(zhì)量移去，再以一定的方式對移去的地形質(zhì)量進行補償，是解算Stokes邊值問題處理地形影響的基本方法。然而，無論用何種方式移去地形質(zhì)量都將使大地水準面發(fā)生變化，產(chǎn)生間接影響。

Molodensky級數(shù)不必考慮地球質(zhì)量的影響，無需作與地球密度相關(guān)的各種改正，它擺脫了Stokes方法需要已知地形密度分布的困難。但事實上地面及其外部擾動位又必然因地形起伏而變，因此其導(dǎo)出的級數(shù)解包含了地形影響項Gn(n=1，2，…)。計算Gn過于復(fù)雜耗時，其級數(shù)各項正負相間，收斂性目前也尚無定論，這就大大降低了Molodensky級數(shù)的實用效果。

利用Stokes-Helmert邊值問題精化似大地水準面對地形影響處理技術(shù)提出了更高要求。目前大地水準面計算中地形影響的計算通常將參考面近似為平面，所有數(shù)據(jù)處理問題都可以在直角坐標系中進行理論推導(dǎo)和計算。然而，對于大區(qū)域乃至全球性的尺度問題，需要利用球坐標進行表達與處理，此時建立在直角坐標系中的數(shù)據(jù)處理方法已不再適用。因此，本文在已有工作的基礎(chǔ)上，推導(dǎo)球坐標的Stokes-Helmert公式和嚴密地形影響計算公式，實現(xiàn)臨沂厘米級大地水準面精化。首先基于確定地球外部擾動重力場的經(jīng)典理論，借鑒國內(nèi)外精密確定地表及其外部擾動引力位的嚴密解算理論，給出Stokes-Helmert邊值問題的定義及Helmert重力場概念；然后根據(jù)實際算例，精確計算第2類Helmert凝聚法中由地形和凝聚層質(zhì)量所產(chǎn)生的直接影響，以及地形及相應(yīng)引力變化的間接影響，構(gòu)建臨沂高精度的似大地水準面模型。

1 Stokes-Helmert邊值問題基本原理

Helmert重力場中的Stokes邊值問題可表示為

(1)

(2)

(3)

式中，g為向量的純量積，將式(3)中的第3項略去，作泰勒級數(shù)展開，保留線性項，則有

(4)

(5)

式(6)在Sg外部成立，用于St上的P點有

(6)

式中，r為地心距離，δA為對重力的直接地形影響。設(shè)點Q為點P沿正常重力線在正常橢球面Se上的對應(yīng)點，則γP可表示為

(7)

式中，N表示大地水準面高；H表示地形高，略去H3等高次項。引進空間改正，令

(8)

則式(7)可寫為

(9)

(10)

式中

(11)

(12)

(13)

2 地形影響

2.1 直接地形影響

移去大地水準面外部地形質(zhì)量的影響稱為地形的直接影響。計算直接影響的經(jīng)典公式為

(14)

式中，G為重力常數(shù)；ρ為地殼密度；l0=2rsinψ/2；Hp和H分別為地面點P和流動點的正常高，R=r。從嚴格意義上來講，式(15)只適用于遠區(qū)積分，即l0≥H。Sj?berg于1995年將直接影響球諧展開至H2，Nahavandchi等1998年將其擴展到H3，Sj?berg于2000年將直接影響寫成曲面積分模式。

對于大區(qū)域的乃至全球尺度問題，發(fā)展橢球坐標系內(nèi)高精度的地形改正方法尤其重要。本文基于國際通用的GRS80橢球采用Tesseroid單元體引力效應(yīng)計算地形改正的方法。對于以地心為圓心，高度H=0、H=h1的兩個球面，其半徑可近似表示為：r1=Rradii+0，r2=Rradii+h1，Rradii表示所選等質(zhì)圓的半徑(這里選地球的平均半徑)。則球面扇形柱體的重力位可表示為

(15)

(16)

由于橢球積分困難，扇形柱體的重力位不能解析計算，但是可通過數(shù)值計算得到近似解，式(15)可以進行級數(shù)展開。

2.2 間接地形影響

恢復(fù)地形質(zhì)量帶來的影響稱為間接影響，計算間接影響的經(jīng)典計算公式為

(17)

式(17)右端第一項稱為主項，第二項稱為次項。Sj?berg于1995年將間接影響球諧展開至H2，Nahavandchi等于1998年將其擴展到H3。然而由于公式需要對整個地球進行積分，這在計算中是不實際的。這里將給出間接改正的嚴密球近似公式。假設(shè)任意一點P的地形引力可寫為

(18)

(19)

式(19)核函數(shù)f(H,t)可表示為

(20)

3 數(shù)據(jù)準備及離散重力數(shù)據(jù)格網(wǎng)化

3.1 數(shù)據(jù)準備

本文研究區(qū)域范圍為34.3°N—36.3°N,117.4°E—119.2°E，地形數(shù)據(jù)采用SRTM3模型。圖1為所用的SRTM3地形數(shù)據(jù)等值線圖，高程最大值為1125 m，平均高程為151.7 m。

圖1 SRTM3地形數(shù)據(jù)

臨沂市境內(nèi)相對重力數(shù)據(jù)為995個，其點位分布如圖2圓點所示。為了檢驗地面重力數(shù)據(jù)的精度，采用該區(qū)18個長期觀測的絕對重力觀測數(shù)據(jù)(其點位分布如圖2中星號點所示)對相對重力測量數(shù)據(jù)進行檢核。經(jīng)分析，測區(qū)相對重力測量數(shù)據(jù)標準差為0.02 mGal。

圖2 地面重力測量點位分布

目前應(yīng)用最為廣泛的全球地球重力場模型是EGM2008地球重力位系數(shù)模型。利用位系數(shù)計算參考重力場的公式為

(21)

本區(qū)共有165個GPS/水準點，分布如圖3所示。其坐標系統(tǒng)采用2000國家大地坐標系，高程基準采用1985國家高程基準。任取其中142個點(圓點標記)作為高一級的控制點計算大地水準面，其余23個點(星號點標記)作為外部檢驗點，以檢驗求得的大地水準面的精度。

3.2 離散重力數(shù)據(jù)格網(wǎng)化

重力測量的結(jié)果是一些分布不規(guī)則的離散點重力值。為了便于使用快速傅里葉變換(FFT)等科學(xué)計算方法，離散重力數(shù)據(jù)需要進行格網(wǎng)化。然而由于空間重力異常的變化較大，直接對空間重力異常進行格網(wǎng)化將產(chǎn)生較大的誤差，因此，通常在進行格網(wǎng)化之前對重力異常進行歸算。本文重力歸算采用布格重力異常進行格網(wǎng)化，具體計算步驟如圖4所示。

圖3 GPS/水準點分布

圖4 重力異常數(shù)據(jù)格網(wǎng)化計算流程

采用上述格網(wǎng)化步驟對地面離散數(shù)據(jù)進行格網(wǎng)化，格網(wǎng)化范圍為34.3°N—36.3°N,117.4°E—119.2°E。臨沂市境內(nèi)采用地面觀測得到的空間異常，臨沂市外空白區(qū)域采用EGM2008模型重力異常進行填充。格網(wǎng)化過程中的地形數(shù)據(jù)采用SRTM3模型。按圖4所述的方法，將離散重力數(shù)據(jù)格網(wǎng)化成2.5′×2.5′分辨率的網(wǎng)格點。圖5給出了格網(wǎng)化后的空間異常等值線圖，其最大值為50.74mGal，最小值為-6.61mGal，平均值為9.36mGal。

圖5 格網(wǎng)化后的地面空間異常

4 臨沂大地水準面模型精化

4.1 臨沂重力似大地水準面計算

采用2.5′×2.5′格網(wǎng)空間重力異常作為輸入數(shù)據(jù)，以EGM2008作為參考重力場，再利用第2類Helmert凝聚法計算大地水準面中的各類地形位及地形引力的影響，即牛頓地形質(zhì)量引力位和凝聚層位間的殘差地形位的間接影響，以及Helmert重力異常由地形質(zhì)量引力位和凝聚層位所產(chǎn)生的引力影響。圖6給出了由Stokes-Helmert法計算的重力似大地水準面，單位為m。

圖6 臨沂市重力似大地水準面等值線

為了評定重力似大地水準面的精度，采用171個GPS/水準資料與其進行比較。結(jié)果表明，重力似大地水準面與GPS/水準高程異常比較，差值最大值為11.43cm，標準差為4.52cm。同樣，采用171個GPS/水準資料與EGM2008模型計算的大地水準面高進行比較，差值最大為21.04cm，標準差為6.67cm。由結(jié)果可見，由Stokes-Helmert方法計算的重力似大地水準面比模型結(jié)果相比，精度有顯著提高。

4.2GPS水準數(shù)據(jù)與重力似大地水準面高擬合

利用高精度GPS/水準數(shù)據(jù)對上述建立的重力似大地水準面進行控制擬合，求得最終的大地水準面模型。由GPS水準測量得到的高程異常及重力似大地水準面起伏之間的擬合模型通?？杀磉_為

H-h-N=Ax+v

(21)

式中，H為GPS測量的大地測高；h為水準測量得到的正常高；N為重力似大地水準面起伏；x為基準不一致的擬合參數(shù)；A為系數(shù)矩陣，取決于擬合模型。本文計算中擬合模型Ax選用多面函數(shù)擬合模型，其高程異常函數(shù)可表示為

(22)

式中，αi為待定系數(shù)；Q(x,y；xi,yi)為核函數(shù)；(x,y)為未測點坐標；(xi,yi)為已測點坐標。核函數(shù)一般可取

(23)

式中，δ稱為平滑因子，為了得到較好的逼近效果，應(yīng)作一定的試算后確定，一般可取作0；b為一個可供選擇的非零實數(shù)，一般取1/2，即核函數(shù)為正雙曲形。

本文采用的GPS/水準點共165個。任意選取142個均勻分布的GPS水準點作為控制點，其他23個點作為檢驗點?？刂泣c和檢核點的點位分布如圖3所示。采用曲面函數(shù)擬合法擬合大地水準面的外符合精度，最大值為3.18 cm，最小值為-2.79 cm，標準差為1.61 cm。如圖7所示。

圖7 臨沂市似大地水準面模型

5 結(jié) 語

本文從討論地形數(shù)據(jù)在大地水準面確定中的作用出發(fā)，分析比較了Stokes和Molodensky兩種經(jīng)典方法處理地形影響存在的問題，給出了Stokes-Helmert邊值問題的定義及數(shù)學(xué)描述，以及地形直接影響、間接影響的理論表達式，介紹了解算Stokes-Helmert邊值問題的詳細步驟和流程。

以格網(wǎng)化后的重力場模型為基礎(chǔ)，采用Stokes-Helmert邊值問題，建立臨沂地區(qū)2.5′×2.5′分辨率的重力似大地水準面模型。利用GPS水準高程異常對重力似大地水準面模型進行控制擬合，求得最終的大地水準面模型。經(jīng)檢核點檢驗可知，擬合外符合精度達到1.6 cm。

[1] HECK B, SEITZ K. A Comparison of the Tesseroid, Prism and Point-mass Approaches for Mass Reductions in Gravity Field Modeling[J]. Journal of Geodesy, 2007, 81(2):121-136.

[2] HEISKANEN W A, MORITZ H. Physical Geodesy[M]. San Francisco: Freeman WH, 1967.

[3] MARTINEC Z, VANEK P. Direct Topographical Effect of Helmert’s Condensation for a Spherical Approximation of the Geoid[J]. Manuscr geod, 1994, 19: 257-268.

[4] NAHAVANDCHI H. On Some Methods of Downward Continua-tion of Mean Free-air Gravity Anomaly[J]. IGES Bull, 1998(8):1-17.

[5] NAHAVANDCHI H, SJ?BERG L E. Terrain Corrections to Power 3 in Gravietric Geoid Determination[J]. Journal of Geodesy 1998, 72(3): 124-135.

[6] NAHAVANDCHI H. The Direct Topographical Correction in Gravimetric Geoid Determination by the Stokes-Helmert Method[J]. Journal of Geodesy, 2000, 74(6): 488-496.

[7] OMANG O C D, FORSBERG R. How to Handle Topography in Practical Geoid Determination: Three Examples [J].Lournal of Geodesy, 2000, 74:458-466.

[8] SJ?BERG L E. On the Quasi-geoid to Geoid Separation[J]. Manuscr Geod, 1995, 20:182-192.

[9] SJ?BERG L E. Topographic Effects by the Stokes-Helmert Method of Geoid and Quasi-geoid Determinations[J]. Journal of Geodesy, 2000, 74(2): 255-268.

[11] WICHIENCHAROEN C. The Indirect Effects on the Computation of Geoid Undulations[M]. Columbus: The Ohio State University, 1982.

[14] 郭東美,鮑李峰,許厚澤.中國大陸厘米級大地水準面的地形影響分析[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版),2016, 41(3):342-348.

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

[16] 羅志才,陳永奇,寧津生. 地形對確定高精度局部大地水準面的影響[J].武漢大學(xué)學(xué)報(信息科學(xué)版), 2003, 28(3): 340-344.

[17] 許厚澤. 我國精化大地水準面工作中若干問題的討論[J]. 地理空間信息,2006,5(4)：1-3.

[18] 章傳銀,晁定波,丁劍,等.厘米級高程異常地形影響的算法及特征分析[J].測繪學(xué)報, 2006, 35(4): 340-344.

High Precision Quasi-geoid Computations Using Improved Stokes-Helmert Boundary Value Problem

LI Deqiang

(Linyi Institute of Land Surveying and Mapping, Linyi 276000, China)

Using the improved Stokes-Helmert boundary value problem, centimeter-order quasi-geoid over Linyi city have been established. Firstly, based on the Stokes theory and Molodensky theory, and making use of the achievements of rigorous solution for disturbing gravitational potential near the surface and exterior of the earth，the mathematical forms of Stokes-Helmert boundary value problem was derived, and the computational formula of direct and indirect topographic effect on the geoid were given. Secondly, the multi-source hetemgeneous data and information were combined to compute the gravimetric quasi-geoid model based on the improved Stokes-Helmert method and the well-known “remove-restore” technology. Finally, the combined adjustment of GPS/leveling height anomaly and gravimetric quasi-geoid data were conducted to obtain the final quasi-geoid model. Numerical results showed that the external accuracy reached 1.6 cm.

quasi-geoid computations; Stokes-Helmert method; multi-surface function fitting; topographic effect

李德強.利用改進的Stokes-Helmert邊值理論解算高精度似大地水準面[J].測繪通報，2017(8)：13-18.

10.13474/j.cnki.11-2246.2017.0246.

2017-04-24

李德強(1975—)，男，高級工程師，研究方向為大地測量、地圖制圖與地理信息工程。E-mail: 280696113@qq.com

P22

A

0494-0911(2017)08-0013-06