利用改進(jìn)的Stokes-Helmert邊值理論解算高精度似大地水準(zhǔn)面

李德強(qiáng)

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

利用改進(jìn)的Stokes-Helmert邊值理論解算高精度似大地水準(zhǔn)面

李德強(qiáng)

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

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

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

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

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

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

1 Stokes-Helmert邊值問(wèn)題基本原理

Helmert重力場(chǎng)中的Stokes邊值問(wèn)題可表示為

(1)

(2)

(3)

式中，g為向量的純量積，將式(3)中的第3項(xiàng)略去，作泰勒級(jí)數(shù)展開(kāi)，保留線(xiàn)性項(xiàng)，則有

(4)

(5)

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

(6)

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

(7)

式中，N表示大地水準(zhǔn)面高；H表示地形高，略去H3等高次項(xiàng)。引進(jìn)空間改正，令

(8)

則式(7)可寫(xiě)為

(9)

(10)

式中

(11)

(12)

(13)

2 地形影響

2.1 直接地形影響

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

(14)

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

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

(15)

(16)

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

2.2 間接地形影響

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

(17)

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

(18)

(19)

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

(20)

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

3.1 數(shù)據(jù)準(zhǔn)備

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

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

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

圖2 地面重力測(cè)量點(diǎn)位分布

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

(21)

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

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

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

圖3 GPS/水準(zhǔn)點(diǎn)分布

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

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

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

4 臨沂大地水準(zhǔn)面模型精化

4.1 臨沂重力似大地水準(zhǔn)面計(jì)算

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

圖6 臨沂市重力似大地水準(zhǔn)面等值線(xiàn)

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

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

利用高精度GPS/水準(zhǔn)數(shù)據(jù)對(duì)上述建立的重力似大地水準(zhǔn)面進(jìn)行控制擬合，求得最終的大地水準(zhǔn)面模型。由GPS水準(zhǔn)測(cè)量得到的高程異常及重力似大地水準(zhǔn)面起伏之間的擬合模型通常可表達(dá)為

H-h-N=Ax+v

(21)

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

(22)

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

(23)

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

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

圖7 臨沂市似大地水準(zhǔn)面模型

5 結(jié) 語(yǔ)

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

以格網(wǎng)化后的重力場(chǎng)模型為基礎(chǔ)，采用Stokes-Helmert邊值問(wèn)題，建立臨沂地區(qū)2.5′×2.5′分辨率的重力似大地水準(zhǔn)面模型。利用GPS水準(zhǔn)高程異常對(duì)重力似大地水準(zhǔn)面模型進(jìn)行控制擬合，求得最終的大地水準(zhǔn)面模型。經(jīng)檢核點(diǎn)檢驗(yàn)可知，擬合外符合精度達(dá)到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] 郭東美,鮑李峰,許厚澤.中國(guó)大陸厘米級(jí)大地水準(zhǔn)面的地形影響分析[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2016, 41(3):342-348.

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

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

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

[18] 章傳銀,晁定波,丁劍,等.厘米級(jí)高程異常地形影響的算法及特征分析[J].測(cè)繪學(xué)報(bào), 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

李德強(qiáng).利用改進(jìn)的Stokes-Helmert邊值理論解算高精度似大地水準(zhǔn)面[J].測(cè)繪通報(bào)，2017(8)：13-18.

10.13474/j.cnki.11-2246.2017.0246.

2017-04-24

李德強(qiáng)(1975—)，男，高級(jí)工程師，研究方向?yàn)榇蟮販y(cè)量、地圖制圖與地理信息工程。E-mail: 280696113@qq.com

P22

A

0494-0911(2017)08-0013-06