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

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

?

VK-Stokes核函數(shù)的性能分析

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

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

為了獲取高精度大地水準(zhǔn)面、克服重力數(shù)據(jù)不能實(shí)現(xiàn)全球覆蓋所帶來(lái)的問(wèn)題，常應(yīng)用修正Stokes核函數(shù)。它能改善標(biāo)準(zhǔn)Stokes核函數(shù)特性，使其在較小積分范圍內(nèi)達(dá)到較高計(jì)算精度。本文基于地球重力場(chǎng)位系數(shù)模型EGM2008，分析了修正核函數(shù)(VK-Stokes)特性及其遠(yuǎn)區(qū)截?cái)嗾`差。計(jì)算分析表明：在近區(qū)0.2°積分半徑內(nèi)，VK-Stokes與標(biāo)準(zhǔn)Stokes和WG-Stokes核函數(shù)值接近；但隨著積分半徑增加，VK-Stokes較標(biāo)準(zhǔn)Stokes核函數(shù)收斂更快，且其遠(yuǎn)區(qū)截?cái)嗾`差數(shù)值也相對(duì)較小。由此可見(jiàn)，應(yīng)用VK-Stokes核函數(shù)，既可實(shí)現(xiàn)在較小積分范圍內(nèi)提高計(jì)算能力，又能抑制其遠(yuǎn)區(qū)截?cái)嗾`差影響。

大地水準(zhǔn)面；Stokes核函數(shù)；VK-Stokes核函數(shù)；截?cái)嗾`差影響

1 引 言

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

修正Stokes核函數(shù)已在多個(gè)大地水準(zhǔn)面計(jì)算實(shí)例中得以應(yīng)用。2005年，A.Ellmann應(yīng)用WG-Stokes函數(shù)和VK-Stokes核函數(shù)，在波羅的海沿岸地區(qū)進(jìn)行了分析[13]。2011年，W.E. Featherstone等人在建立最新澳大利亞重力水準(zhǔn)面模型AUSGeoid09時(shí)，對(duì)球面Stokes核函數(shù)(S-Stokes)、WG-Stokes和FEO-Stokes進(jìn)行了應(yīng)用分析[14]。2012年，Y.M. Wang等在美國(guó)重力大地水準(zhǔn)面模型(USGG2009)建立中，采用了WG-Stokes核函數(shù)[15]。2012年，Huang Jianliang 綜合了WG-Stokes和限階次Stokes函數(shù)，提出新的修正核函數(shù)(MDBK)，并將其應(yīng)用于加拿大最新重力大地水準(zhǔn)面模型建立[16]。2013年，傅露等結(jié)合DNSC08- CRA模型中的美國(guó)近海測(cè)高重力數(shù)據(jù)，分析比較了5種修正Stokes核函數(shù)的計(jì)算精度，認(rèn)為修正Stokes核函數(shù)可有效改善計(jì)算精度[17]。為了進(jìn)一步提高局部區(qū)域(似)大地水準(zhǔn)面計(jì)算精度，更恰當(dāng)?shù)厥褂肰K-Stoke核函數(shù)，本文分析了VK-Stokes函數(shù)特性及其遠(yuǎn)區(qū)的截?cái)嗾`差，并與標(biāo)準(zhǔn)Stokes和WG-Stokes核函數(shù)進(jìn)行比較。

2 大地水準(zhǔn)面差距計(jì)算

所謂大地水準(zhǔn)面差距是指大地水準(zhǔn)面與參考橢球面間的距離。通過(guò)大地水準(zhǔn)面差距，可實(shí)現(xiàn)大地高與正高間的轉(zhuǎn)化。大地水準(zhǔn)面差距N的計(jì)算表達(dá)式為[1-2]：

(1)

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

將積分區(qū)劃分為近區(qū)(σ0)和遠(yuǎn)區(qū)σ-(σ0),則有：

(2)

基于移去-恢復(fù)法，則有：

(3)

式中，NM為低階位系數(shù)模型計(jì)算大地水準(zhǔn)面差距；Δgres為用地球重力場(chǎng)模型移去長(zhǎng)波項(xiàng)后，剩余重力異常值。

低階位系數(shù)模型計(jì)算重力異常公式如下：

(4)

盡管采用移去-恢復(fù)法，基于地球重力場(chǎng)位系數(shù)模型移去長(zhǎng)波項(xiàng)，但是剩余重力異常中還會(huì)時(shí)常存在殘余長(zhǎng)波信號(hào)。另外，剩余重力異常在有限積分范圍內(nèi)進(jìn)行解算，其遠(yuǎn)區(qū)影響通常被直接忽略，引入了一定誤差。考慮應(yīng)用VK-Stokes核函數(shù)，對(duì)于解決上述問(wèn)題有一定程度改善。

3 修正的Stokes核函數(shù)

3.1WG-Stokes核函數(shù)

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

(5)

式中，Ls為截?cái)嚯A次；Pn(cosΨ)為L(zhǎng)egendre多項(xiàng)式。

或?qū)懗桑?/p>

(6)

(7)

(8)

3.2VK-Stokes核函數(shù)

1987年，Vanicek和Kleusberg基于Molodensky思想修正了WG-Stokes核函數(shù)，使截?cái)嗾`差上限最小，給出了VK-Stokes核函數(shù)，其具體形式為：

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

將上式改寫(xiě)為：

(10)

(11)

(12)

(13)

求解待定系數(shù)tn(Ψ0)(n=0,1……,Lm)時(shí)，若使截?cái)嗾`差上限最小，則必須令下式最小。

(14)

(15)

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

4 遠(yuǎn)區(qū)截?cái)嗾`差

遠(yuǎn)區(qū)截?cái)嗾`差是指剩余重力異常的遠(yuǎn)區(qū)影響，其計(jì)算表達(dá)式為：

(16)

WG-Stokes核函數(shù)相應(yīng)截?cái)嗾`差影響計(jì)算公式為：

(17)

VK-Stokes核函數(shù)相應(yīng)截?cái)嗾`差影響計(jì)算公式為：

(18)

基于移去-恢復(fù)法，遠(yuǎn)區(qū)影響一般可由地球重力場(chǎng)位系數(shù)模型近似求得，其主要誤差源于位系數(shù)誤差和有限截?cái)嚯A次引起的誤差。位系數(shù)誤差隨模型建立而引入，為已知量。

5 實(shí)驗(yàn)分析

5.1 核函數(shù)特性分析

為了更好地應(yīng)用VK-Stokes函數(shù)，基于EGM2008位系數(shù)模型，比較分析其與標(biāo)準(zhǔn)Stokes核函數(shù)以及WG-Stokes核函數(shù)，與截?cái)嚯A次、修正階次和積分半徑間的關(guān)系，見(jiàn)圖1(a)～(f)所示。圖1(a)和(b)分別為近區(qū)、標(biāo)準(zhǔn)Stokes函數(shù)與截?cái)嗖⑿拚?0和120階的WG-Stokes函數(shù)及VK-Stokes函數(shù)，隨著積分半徑變化的情況；圖1(c)和(d)分別為遠(yuǎn)區(qū)、標(biāo)準(zhǔn)Stokes函數(shù)與截?cái)嗖⑿拚?0和120階的WG-Stokes函數(shù)及VK-Stokes函數(shù)，隨著積分半徑變化的情況；圖1(e)為近區(qū)，相同修正階次、不同截?cái)嚯A次下，VK-Stokes函數(shù)隨積分半徑變化的情況；圖1(f)為近區(qū)，相同截?cái)嚯A次、不同修正階次下，VK-Stokes核函數(shù)隨積分半徑變化的情況。

圖1 Stokes核函數(shù)特性圖

由圖1(a)和(b)可見(jiàn)，積分半徑約在0.2°以?xún)?nèi)，三函數(shù)值最為接近。隨著積分半徑增加，Stokes核函數(shù)與標(biāo)準(zhǔn)Stokes核函數(shù)值差異增大。截?cái)嚯A次越高，VK-Stokes和WG-Stokes收斂速度則越快。VK-Stokes核函數(shù)是在WG-Stokes核函數(shù)基礎(chǔ)上，使得截?cái)嗌舷逎M(mǎn)足最小，因此與WG-Stokes核函數(shù)在近區(qū)較為接近。

由圖1(c)和(d)可見(jiàn)，遠(yuǎn)區(qū)、VK-Stokes較WG-Stokes收斂快，且比標(biāo)準(zhǔn)Stokes核函數(shù)計(jì)算值小。由此看來(lái)，VK-Stokes核函數(shù)能夠減弱遠(yuǎn)區(qū)影響。

由圖1(e)和(f)可見(jiàn)，截?cái)嚯A次對(duì)VK-Stokes核函數(shù)的影響比修正階次對(duì)其影響大。選用較低修正階次，可提高計(jì)算速度。

5.2 截?cái)嗾`差系數(shù)分析

基于EGM2008位系數(shù)模型，分析截?cái)嗾`差系數(shù)與修正階次、截?cái)嚯A次以及積分半徑的關(guān)系。圖2(a)和(b)分別為L(zhǎng)s=Lm=20和120階，不同積分半徑(1°、3°和6°)下，VK-Stokes核函數(shù)相應(yīng)的截?cái)嗾`差系數(shù)隨模型階次變化圖。圖2(c)為積分半徑6°，Ls=Lm=20，WG-Stokes和VK-Stokes函數(shù)相應(yīng)的截?cái)嗾`差系數(shù)隨模型階次變化圖。

圖2 截?cái)嗾`差系數(shù)圖

如圖2(a)和(b)可見(jiàn)，VK-Stokes核函數(shù)截?cái)嗾`差系數(shù)受截?cái)嚯A次影響較大。隨著積分半徑Ψ0增加，收斂速度加快；但當(dāng)積分半徑超過(guò)3°時(shí)，其收斂速度減緩。

由圖2(c)可知，VK-Stokes核函數(shù)截?cái)嗾`差系數(shù)比WG-Stokes核函數(shù)截?cái)嗾`差系數(shù)收斂快。

5.3 截?cái)嗾`差分析

基于EGM2008地球重力場(chǎng)位系數(shù)模型，以B=30.72211，L=110.40191,H=1045.3m為例，采用VK-Stokes核函數(shù)，選擇Ls=Lm=20，積分半徑1°和3°，計(jì)算遠(yuǎn)區(qū)截?cái)嗾`差,見(jiàn)圖3(a)所示；選擇Ls=120；Lm=20和Ls=Lm=120，積分半徑1°和3°，分別計(jì)算遠(yuǎn)區(qū)截?cái)嗾`差，見(jiàn)圖3(b)和(c)。

圖3 遠(yuǎn)區(qū)截?cái)嗾`差

由圖3(a)可見(jiàn)，截?cái)嗖⑿拚?0階，其截?cái)嗾`差在幾個(gè)厘米量級(jí)，這對(duì)于建立1cm精度大地水準(zhǔn)面來(lái)說(shuō)不恰當(dāng)，還需增大截?cái)嚯A次或者積分半徑來(lái)提高計(jì)算精度。

由圖3(b)和(c)可知，截?cái)嘀?20階，其截?cái)嗾`差在1cm以?xún)?nèi)。積分半徑對(duì)VK-Stokes的截?cái)嗾`差影響較小，因此可選較小積分半徑。隨著積分半徑增加，VK-Stokes核函數(shù)作用減弱。在實(shí)際截?cái)嗾`差當(dāng)中，還包括重力測(cè)量誤差以及地球重力場(chǎng)位系數(shù)等誤差，這些誤差為固定誤差。

6 小 結(jié)

區(qū)域或局部(似)大地水準(zhǔn)面精化工作中，常受到數(shù)據(jù)量稀少的限制，長(zhǎng)波以及系統(tǒng)性誤差污染，給計(jì)算帶來(lái)不便，影響計(jì)算精度。通過(guò)修正Stokes核函數(shù)，可增強(qiáng)近區(qū)計(jì)算能力，削弱遠(yuǎn)區(qū)影響，抑制其他污染源影響。

VK-Stokes核函數(shù)既可有效利用有限的地面重力測(cè)量數(shù)據(jù)，還能減弱遠(yuǎn)區(qū)影響，但其作用會(huì)隨著積分范圍增大而減緩。從計(jì)算效率角度考慮，選用較高截?cái)嚯A次和較小積分半徑更為合適。修正階次對(duì)計(jì)算精度影響不大，但其階次選取影響計(jì)算速度，因此可選用較低階次。

相對(duì)而言，VK-Stokes核函數(shù)截?cái)嗾`差系數(shù)比WG-Stokes和標(biāo)準(zhǔn)Stokes核函數(shù)截?cái)嗾`差系數(shù)收斂快，且截?cái)嗾`差影響數(shù)值也較小，因此選用VK-Stokes核函數(shù)相對(duì)合理。值得注意的是，本實(shí)驗(yàn)只反映了單項(xiàng)截?cái)嗾`差，尚未顧及重力異常觀測(cè)誤差，以及參考模型位系數(shù)誤差等因素，在實(shí)際應(yīng)用當(dāng)中還需進(jìn)一步分析研究。

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

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

[3]李建成．我國(guó)現(xiàn)代高程測(cè)定關(guān)鍵技術(shù)若干問(wèn)題的研究及進(jìn)展[J]．武漢大學(xué)學(xué)報(bào)· 信息科學(xué)版，2007，32(11)：980-987．

[4]李建成．最新中國(guó)陸地?cái)?shù)字高程基準(zhǔn)模型：重力似大地水準(zhǔn)面CNGG2011[J]．測(cè)繪學(xué)報(bào)，2012，41(5)：651-660．

[5]魏子卿，王剛．用地球位模型和GPS /水準(zhǔn)數(shù)據(jù)確定我國(guó)大陸似大地水準(zhǔn)面[J]．測(cè)繪學(xué)報(bào)，2003，32(1)：1-5．[6]魏子卿．大地水準(zhǔn)面短議[J]．地理空間信息，2009，7(1)：1-3．

[7]黃謨濤，翟國(guó)君，管錚等．海洋重力場(chǎng)測(cè)定及其應(yīng)用[M]．北京：測(cè)繪出版社，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]傅露，褚永海．區(qū)域大地水準(zhǔn)面確定中Stokes核函數(shù)的應(yīng)用[J]．大地測(cè)量與地球動(dòng)力學(xué)，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。

國(guó)家自然科學(xué)基金資助項(xiàng)目(41174018; 41304022;41474015)。

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

P223

A