VK-Stokes核函數的性能分析

榮 敏,周 巍,吳富梅,任紅飛

1.信息工程大學地理空間信息學院，河南 鄭州，450052；2.西安測繪研究所，陜西 西安，710054；3. 測繪信息技術總站，陜西 西安，710054

?

VK-Stokes核函數的性能分析

榮 敏1，2，3,周 巍3,吳富梅2,任紅飛2

1.信息工程大學地理空間信息學院，河南 鄭州，450052；2.西安測繪研究所，陜西 西安，710054；3. 測繪信息技術總站，陜西 西安，710054

為了獲取高精度大地水準面、克服重力數據不能實現全球覆蓋所帶來的問題，常應用修正Stokes核函數。它能改善標準Stokes核函數特性，使其在較小積分范圍內達到較高計算精度。本文基于地球重力場位系數模型EGM2008，分析了修正核函數(VK-Stokes)特性及其遠區截斷誤差。計算分析表明：在近區0.2°積分半徑內，VK-Stokes與標準Stokes和WG-Stokes核函數值接近；但隨著積分半徑增加，VK-Stokes較標準Stokes核函數收斂更快，且其遠區截斷誤差數值也相對較小。由此可見，應用VK-Stokes核函數，既可實現在較小積分范圍內提高計算能力，又能抑制其遠區截斷誤差影響。

大地水準面；Stokes核函數；VK-Stokes核函數；截斷誤差影響

1 引 言

應用Stokes公式計算大地水準面時，需要在全球范圍內進行求解；但由于數據、成本和計算效率等原因，實現全球積分不太現實[1-15]。為了克服困難、提高計算精度和效率，學者們提出了多種改善方法。其中，將積分域劃分為球冠(近區)和剩余部分(遠區)兩部分[5]；再將重力異常值按照頻譜劃分為高頻和低頻兩部分，用高精度地球重力場模型移去低頻影響，對剩余重力異常采用Stokes公式進行解算。若直接采用標準Stokes核函數，可能會產生波長扭曲現象，因此需要修正標準Stokes核函數[5]。修正Stokes核函數，一則可以改變標準Stokes函數特性，二則還可使其遠區截斷誤差迅速減小，最為理想的狀態是其值為0。

修正Stokes核函數已在多個大地水準面計算實例中得以應用。2005年，A.Ellmann應用WG-Stokes函數和VK-Stokes核函數，在波羅的海沿岸地區進行了分析[13]。2011年，W.E. Featherstone等人在建立最新澳大利亞重力水準面模型AUSGeoid09時，對球面Stokes核函數(S-Stokes)、WG-Stokes和FEO-Stokes進行了應用分析[14]。2012年，Y.M. Wang等在美國重力大地水準面模型(USGG2009)建立中，采用了WG-Stokes核函數[15]。2012年，Huang Jianliang 綜合了WG-Stokes和限階次Stokes函數，提出新的修正核函數(MDBK)，并將其應用于加拿大最新重力大地水準面模型建立[16]。2013年，傅露等結合DNSC08- CRA模型中的美國近海測高重力數據，分析比較了5種修正Stokes核函數的計算精度，認為修正Stokes核函數可有效改善計算精度[17]。為了進一步提高局部區域(似)大地水準面計算精度，更恰當地使用VK-Stoke核函數，本文分析了VK-Stokes函數特性及其遠區的截斷誤差，并與標準Stokes和WG-Stokes核函數進行比較。

2 大地水準面差距計算

所謂大地水準面差距是指大地水準面與參考橢球面間的距離。通過大地水準面差距，可實現大地高與正高間的轉化。大地水準面差距N的計算表達式為[1-2]：

(1)

式中，R為地球平均半徑；γ0為正常重力；Δg為大地水準面上重力異常值；Ψ為計算點與積分點之間的球面角距；積分單元dσ=sinΨdΨdα；S(Ψ)表示Stokes函數。

將積分區劃分為近區(σ0)和遠區σ-(σ0),則有：

(2)

基于移去-恢復法，則有：

(3)

式中，NM為低階位系數模型計算大地水準面差距；Δgres為用地球重力場模型移去長波項后，剩余重力異常值。

低階位系數模型計算重力異常公式如下：

(4)

盡管采用移去-恢復法，基于地球重力場位系數模型移去長波項，但是剩余重力異常中還會時常存在殘余長波信號。另外，剩余重力異常在有限積分范圍內進行解算，其遠區影響通常被直接忽略，引入了一定誤差。考慮應用VK-Stokes核函數，對于解決上述問題有一定程度改善。

3 修正的Stokes核函數

3.1WG-Stokes核函數

WG-Stokes核函數是1969年由Wong和Gore給出，它是直接將球面Stokes核函數值剔除其低階部分。用SWG(Ψ)表示WG-Stokes核函數，則有：

(5)

式中，Ls為截斷階次；Pn(cosΨ)為Legendre多項式。

或寫成：

(6)

(7)

(8)

3.2VK-Stokes核函數

1987年，Vanicek和Kleusberg基于Molodensky思想修正了WG-Stokes核函數，使截斷誤差上限最小，給出了VK-Stokes核函數，其具體形式為：

式中，tn(Ψ0)為修正系數；Lm為修正階次。

將上式改寫為：

(10)

(11)

(12)

(13)

求解待定系數tn(Ψ0)(n=0,1……,Lm)時，若使截斷誤差上限最小，則必須令下式最小。

(14)

(15)

采用數值方法計算求解tn(Ψ0)。

4 遠區截斷誤差

遠區截斷誤差是指剩余重力異常的遠區影響，其計算表達式為：

(16)

WG-Stokes核函數相應截斷誤差影響計算公式為：

(17)

VK-Stokes核函數相應截斷誤差影響計算公式為：

(18)

基于移去-恢復法，遠區影響一般可由地球重力場位系數模型近似求得，其主要誤差源于位系數誤差和有限截斷階次引起的誤差。位系數誤差隨模型建立而引入，為已知量。

5 實驗分析

5.1 核函數特性分析

為了更好地應用VK-Stokes函數，基于EGM2008位系數模型，比較分析其與標準Stokes核函數以及WG-Stokes核函數，與截斷階次、修正階次和積分半徑間的關系，見圖1(a)～(f)所示。圖1(a)和(b)分別為近區、標準Stokes函數與截斷并修正至20和120階的WG-Stokes函數及VK-Stokes函數，隨著積分半徑變化的情況；圖1(c)和(d)分別為遠區、標準Stokes函數與截斷并修正至20和120階的WG-Stokes函數及VK-Stokes函數，隨著積分半徑變化的情況；圖1(e)為近區，相同修正階次、不同截斷階次下，VK-Stokes函數隨積分半徑變化的情況；圖1(f)為近區，相同截斷階次、不同修正階次下，VK-Stokes核函數隨積分半徑變化的情況。

圖1 Stokes核函數特性圖

由圖1(a)和(b)可見，積分半徑約在0.2°以內，三函數值最為接近。隨著積分半徑增加，Stokes核函數與標準Stokes核函數值差異增大。截斷階次越高，VK-Stokes和WG-Stokes收斂速度則越快。VK-Stokes核函數是在WG-Stokes核函數基礎上，使得截斷上限滿足最小，因此與WG-Stokes核函數在近區較為接近。

由圖1(c)和(d)可見，遠區、VK-Stokes較WG-Stokes收斂快，且比標準Stokes核函數計算值小。由此看來，VK-Stokes核函數能夠減弱遠區影響。

由圖1(e)和(f)可見，截斷階次對VK-Stokes核函數的影響比修正階次對其影響大。選用較低修正階次，可提高計算速度。

5.2 截斷誤差系數分析

基于EGM2008位系數模型，分析截斷誤差系數與修正階次、截斷階次以及積分半徑的關系。圖2(a)和(b)分別為Ls=Lm=20和120階，不同積分半徑(1°、3°和6°)下，VK-Stokes核函數相應的截斷誤差系數隨模型階次變化圖。圖2(c)為積分半徑6°，Ls=Lm=20，WG-Stokes和VK-Stokes函數相應的截斷誤差系數隨模型階次變化圖。

圖2 截斷誤差系數圖

如圖2(a)和(b)可見，VK-Stokes核函數截斷誤差系數受截斷階次影響較大。隨著積分半徑Ψ0增加，收斂速度加快；但當積分半徑超過3°時，其收斂速度減緩。

由圖2(c)可知，VK-Stokes核函數截斷誤差系數比WG-Stokes核函數截斷誤差系數收斂快。

5.3 截斷誤差分析

基于EGM2008地球重力場位系數模型，以B=30.72211，L=110.40191,H=1045.3m為例，采用VK-Stokes核函數，選擇Ls=Lm=20，積分半徑1°和3°，計算遠區截斷誤差,見圖3(a)所示；選擇Ls=120；Lm=20和Ls=Lm=120，積分半徑1°和3°，分別計算遠區截斷誤差，見圖3(b)和(c)。

圖3 遠區截斷誤差

由圖3(a)可見，截斷并修正至20階，其截斷誤差在幾個厘米量級，這對于建立1cm精度大地水準面來說不恰當，還需增大截斷階次或者積分半徑來提高計算精度。

由圖3(b)和(c)可知，截斷至120階，其截斷誤差在1cm以內。積分半徑對VK-Stokes的截斷誤差影響較小，因此可選較小積分半徑。隨著積分半徑增加，VK-Stokes核函數作用減弱。在實際截斷誤差當中，還包括重力測量誤差以及地球重力場位系數等誤差，這些誤差為固定誤差。

6 小 結

區域或局部(似)大地水準面精化工作中，常受到數據量稀少的限制，長波以及系統性誤差污染，給計算帶來不便，影響計算精度。通過修正Stokes核函數，可增強近區計算能力，削弱遠區影響，抑制其他污染源影響。

VK-Stokes核函數既可有效利用有限的地面重力測量數據，還能減弱遠區影響，但其作用會隨著積分范圍增大而減緩。從計算效率角度考慮，選用較高截斷階次和較小積分半徑更為合適。修正階次對計算精度影響不大，但其階次選取影響計算速度，因此可選用較低階次。

相對而言，VK-Stokes核函數截斷誤差系數比WG-Stokes和標準Stokes核函數截斷誤差系數收斂快，且截斷誤差影響數值也較小，因此選用VK-Stokes核函數相對合理。值得注意的是，本實驗只反映了單項截斷誤差，尚未顧及重力異常觀測誤差，以及參考模型位系數誤差等因素，在實際應用當中還需進一步分析研究。

[1]陸仲連．地球重力場理論與方法[M]．北京：解放軍出版社，1996．

[2]Hofmann-Wellenhof B. and H. Moritz, Physical geodesy [M].second edition, Springer Wien New York, 2006.

[3]李建成．我國現代高程測定關鍵技術若干問題的研究及進展[J]．武漢大學學報· 信息科學版，2007，32(11)：980-987．

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

[5]魏子卿，王剛．用地球位模型和GPS /水準數據確定我國大陸似大地水準面[J]．測繪學報，2003，32(1)：1-5．[6]魏子卿．大地水準面短議[J]．地理空間信息，2009，7(1)：1-3．

[7]黃謨濤，翟國君，管錚等．海洋重力場測定及其應用[M]．北京：測繪出版社，2005.

[8]L．E．Sj?berg，A．Hunegnaw．Some modifications of Stokes formula the account for truncation and potential coefficient errors[J]．Journal of Geodesy，2000，74(3)：232-238．

[9]L．E．Sj?berg．A computational scheme to model the geoid by the modified Stokes formula without gravity reductions[J]．Journal of Geodesy，2003，77(4)：423-432.

[10]W．E．Featherstone，J．D．Evans，J．G．Olliver．A Meissl-modified Vanicek and Kleusberg kernel toreduce the truncation error in gravimetric geoidcomputations[J]．Journal of Geodesy，1998，72(3)：154-160．

[11]P．Vanicek，W．E．Featherstone. Performance of three types of Stokes’s kernel in the combined solution for the geoid[J]．Journal of Geodesy，1998，72(9)：684-697．

[12]J．D．Evans，W．E． Featherstone．Improved convergence rates for the truncation error in gravimetric geoid determination[J]．Journal of Geodesy，2000，74(2)：239-248．

[13]A. Ellmann. Two deterministic and three stochastic modifications of Stokes’s formula：a case study for the Baltic countries[J]．Journal of Geodesy，2005，79(1)：11-23．

[14]W. E. Featherstone， J. F. Kirby， C. Hirt，et al. The AUSGeoid09 model of the Australian Height Datum[J]．Journal of Geodesy，2011，85(3)：133-150．

[15]Y. M. Wang, J.Saleh, X. Li,et al. The US Gravimetric Geoid of 2009(USGG2009):model development and evaluation[J]．Journal of Geodesy，2012，86(3)：165-180．

[16]J．Huang，M．Veronneau．Canadian gravimetric geoid model 2010[J]．Journal of Geodesy，2013，87(9)：771-790．

[17]傅露，褚永海．區域大地水準面確定中Stokes核函數的應用[J]．大地測量與地球動力學，2013，33(2)：110-113．

Performance Analysis of the VK-Stokes Kernel Function

Rong Min1，2，3,Zhou Wei3,Wu Fumei2,Ren Hongfei2

1. Institute of Geospatial Information，Information Engineering University, Zhengzhou 450052,China 2. Xi’an Research Institute of Survey and Mapping, Xi’an 710054,China 3. Technical Division of Surveying and Mapping，Xi’an 710054,China

In order to get high accuracy of the geoid and to overcome the problem that the gravity data cannot cover the whole world, modified Stokes kernel function is often used to deal with the difficulty. It can improve Stokes kernel function characteristics and make it achieve high accuracy in the small integral range. The characteristics and the far-region truncation error of the VK-Stokes kernel function are analyzed in this paper based on the gravity potential coefficient of the model EGM2008. The results show that in the near-region of the integral radius 0.2°，the VK-Stokes value is close to that of Stokes and WG-Stokes. But with the increasing of integral radius, VK-Stokes converges faster than Stokes and the truncation error of the VK-Stokes kernel function is relatively small in the far-region. Thus it shows the VK-Stokes not only improves the ability of computation in the small integral range, but also controls the truncation error in the far-region.

geoid；Stokes kernel function；VK- Stokes kernel function；truncation error

2015-01-13。

國家自然科學基金資助項目(41174018; 41304022;41474015)。

榮敏(1977—)，女，工程師，主要從事重力大地水準面計算方法的研究。

P223

A