劉纘武，劉世晗，黃 歐

1.信息工程大學(xué)理學(xué)院，河南鄭州450001；2.鄭州市市政工程總公司，河南鄭州450007；3.美國俄亥俄州立大學(xué)地球科學(xué)學(xué)院，俄亥俄哥倫布43210

1 引 言

由于當(dāng)今衛(wèi)星跟蹤技術(shù)、衛(wèi)星測高以及重力場求解技術(shù)的發(fā)展，國內(nèi)外正在把高階次重力場模型向更高階次擴(kuò)展［1－5］。顯然，超高階次球諧位系數(shù)模型的構(gòu)制與應(yīng)用，與超高階次締合勒讓德函數(shù)及其導(dǎo)數(shù)的計(jì)算密切相關(guān)。并且，隨著航天技術(shù)的不斷發(fā)展，需要構(gòu)建的位系數(shù)模型的階次越來越高。因此，研究計(jì)算超高階次締合勒讓德函數(shù)及其導(dǎo)數(shù)值的方法是非常必要的［4，6－7］。

由締合勒讓德函數(shù)及其關(guān)于θ的導(dǎo)數(shù)構(gòu)成的球諧級(jí)數(shù)表達(dá)式為

和

式中，（r，θ，λ）是空間點(diǎn)的球面坐標(biāo)，r是點(diǎn)至地心的距離，θ是地心余維，λ是經(jīng)度；M是球諧展開式的最大階數(shù)；Enm（r，λ）是與階n和次m有關(guān)的（r，λ）的函數(shù)；和是n階m次第一類完全正?；喓侠兆尩潞瘮?shù)和它們的導(dǎo)數(shù)。

例如，在點(diǎn)（r，θ，λ）處，地球外部擾動(dòng)位T（r，θ，λ）的M階球諧級(jí)數(shù)表達(dá)式為［8］

這里

式中，

式中

按式（4）～式（8）遞推計(jì)算，當(dāng)M＝2 500在接近兩極時(shí)函數(shù)值達(dá)到10－4000量級(jí)［4，6］。而一般個(gè)人計(jì)算機(jī)的允許數(shù)值范圍約在10－310和10＋310之間［10］。為了避免下溢，文獻(xiàn)［4］修改遞推式（4）和式（7）為

計(jì)算式（1）采用如下Horner求和技術(shù)［4，12］

式中，

式（2）可以類似計(jì)算。

的值當(dāng)M＝2 700，θ→0°時(shí)將超過10500。為了避免上溢，遞推過程中乘以一個(gè)下調(diào)因子［4］10－280。

2 遞推計(jì)算中的壓縮因子和復(fù)合Horner求和技術(shù)

分析式（9）中的系數(shù)因子t、αnm和βnm不難發(fā)現(xiàn)，對確定的θ和m，當(dāng)n很大時(shí)，βnm接近于1；當(dāng)θ→0°時(shí)t→1，且它不隨m和n的變化而變化；而隨著n的增大αnm趨于2（對于固定的m），且當(dāng)m＜n＜1.666 7 m時(shí)αnm＞2.5。因此，每一列遞推，n由m增大到n＝1.666 7 m，每一步遞推值將增大超過αnm／βnm＝1.5倍。例如，當(dāng)n由m＝ 2 700增大到n＝1.666 7 m（約4 500）時(shí)，如果不顧及t的微小影響，遞推值將遠(yuǎn)遠(yuǎn)超過初始值的9×10316（實(shí)際超過1×10500）倍以上。

鑒于以上分析，在遞推中插入壓縮因子q，式（9）變?yōu)?/p>

式中，q是與θ有關(guān)的介于1和2之間的一個(gè)數(shù)。顯然，與式（9）相比，式（12）降低了遞推值增大的速度。

當(dāng)θ逐漸增大到90°時(shí)，t＝cosθ逐漸減小到0，式（9）右端第一項(xiàng)作用越來越小，遞推值主要依賴于第二項(xiàng)。為避免下溢，這時(shí)取q＝1。由此，定義q為

式中，ρ是一個(gè)介于0.5和1之間的常數(shù)；k是一個(gè)正整數(shù)。令

式（12）改寫為

式中，

式（1）的值如下計(jì)算

式中，

式（18）和式（19）合稱為復(fù)合Horner求和技術(shù)。

式（2）中f（1）可以類似計(jì)算。

式（20）和式（21）也可合稱為復(fù)合Horner求和技術(shù)。

3 數(shù)值檢驗(yàn)

圖1 ｜（θ）｜（?n，m≤3 600）的最大值的對數(shù)值曲線Fig.1 Logarithm plot of maximum values of｜（θ）｜（?n，m≤3 600）

圖2 ｜（θ）｜（?n，m≤3 600）的最小值的對數(shù)值曲線Fig.2 Logarithm plot of minimum values of｜（θ）｜（?n，m≤3 600）

式中，

式（23）、式（24）類似于式（20）、式（21）計(jì)算。

從圖3看到，新算法求得的完全正則化締合勒讓德函數(shù)值具有較高的精度，同時(shí)圖1～圖3也說明新算法有較好的穩(wěn)定性。

Belikov遞推法的計(jì)算公式為

式中，

跨階次遞推法的計(jì)算公式為

式中，

由于式（25）和式（26）的遞推系數(shù)都小于1，因此比式（4）的穩(wěn)定性要強(qiáng)。但在接近兩極時(shí)的值仍然達(dá)到10－4000量級(jí)，超出一般計(jì)算機(jī)的允許數(shù)值范圍而下溢。根據(jù)式（25）和式（26）的結(jié)構(gòu)，不能像式（9）那樣處理。不過，式（25）可改寫為

由式（27）計(jì)算＾Pnm（θ）／um可避免下溢，且若遞推中乘以10－280，可計(jì)算至2 500階不會(huì)上溢。由此轉(zhuǎn)化為就可類似式（11）計(jì)算式（1）。但當(dāng)M＞2 500時(shí)出現(xiàn)上溢，若要進(jìn)一步提高計(jì)算的階次，因＾Pnm（θ）／um由三項(xiàng)構(gòu)成，則壓縮因子結(jié)構(gòu)非常復(fù)雜。

跨階次遞推式（26）若如法處理，則產(chǎn)生系數(shù)βnm／u2和γnm／u2，當(dāng)θ→0°它們急劇變大，雖然避免了下溢，但很快出現(xiàn)上溢。由此分析，利用本文的處理思想，當(dāng)階次很高時(shí)，無法避免遞推值的溢出。

4 結(jié) 論

在個(gè)人計(jì)算機(jī)上執(zhí)行現(xiàn)有的遞推公式計(jì)算超過M＞2 700階次的締合勒讓德函數(shù)會(huì)產(chǎn)生上溢或下溢。通過在基本遞推算法中插入壓縮因子法，避免了溢出現(xiàn)象的發(fā)生。同時(shí)給出了計(jì)算由完全正則化締合勒讓德函數(shù)構(gòu)成的超高階球諧級(jí)數(shù)式的復(fù)合Horner求和技術(shù)。數(shù)值檢驗(yàn)表明新算法具有較高的精度。

毫無疑問，在大地測量中會(huì)有包含真實(shí)擾動(dòng)位系數(shù)在內(nèi)的計(jì)算締合勒讓德函數(shù)級(jí)數(shù)f的更嚴(yán)密的檢測方法。本文重點(diǎn)在于提出締合勒讓德函數(shù)遞推的壓縮因子算法和復(fù)合Horner求和技術(shù)的思想。依據(jù)這個(gè)思路，Belikov遞推法可以類似處理，且遞推計(jì)算更高階的締合勒讓德函數(shù)級(jí)數(shù)及確定相應(yīng)的級(jí)數(shù)求和技術(shù)，都會(huì)涉及更為復(fù)雜的壓縮因子和球諧級(jí)數(shù)求和手段；任意高階的締合勒讓德函數(shù)遞推算法及其相應(yīng)的級(jí)數(shù)求和問題，必須要有新的思想和方法；跨階次遞推法很難利用本文給出的方法處理下溢問題。對以上問題，將進(jìn)一步研究。

［1］ WENZEL G.Ultra－h(huán)igh Degree Geopotential Models GPM98A，B，and C to Degree 1 800［C］∥Joint Meeting of the International Gravity Commission and International Geoid Commission.Trieste：［s.n.］，1998.

［2］ LIU Zuanwu，LU Zhonglian.A Study of Incomplete Geopotential Coefficient Model［J］.Geophysics，1998，41（1）：89－98.

［3］ HUANG Motao，ZHAI Guojun，OUYANG Yongzhong，et al.Analysis and Computation of Ultra High Degree Geopotential Model［J］.Acta Geotaetica et Cartographica Sinica，2001，30（3）：208－213.（黃漠濤，翟國君，歐陽永忠，等.超高階地球位模型的計(jì)算與分析［J］.測繪學(xué)報(bào)，2001，30（3）：208－213.）

［4］ HOLMES S A，F(xiàn)EATHERSTONE W E.A Unified Approach to the Clenshaw Summation and the Recursive Computation of Very High Degree and Order Normalised Associated Legendre Functions［J］.Journal of Geodesy， 2002，76（5）：279－299.

［5］ PAVLIS N K，HOLMES S A，KENYON S C，et al.Earth Gravitational Model 2008（EGM2008）－WGS－84Version［EB／OL］.［2010－02－20］.http：∥earth－info.nga.mil／GandG／WGS－84／gravitymod／egm2008／index.html.

［6］ JEKELI C，LEE J K，KWON J H.On the Computation and Approximation of Ultra－h(huán)igh Degree Spherical Harmonic Seies［J］.Journal of Geodesy，2007，81（9）：603－615.

［7］ WANG Jianqiang，ZHAO Guoqiang，ZHU Guangbin.Contrastive Analysis of Common Computing Methods of Ultra－h(huán)igh Degree and Order Fully Normalized Associated Legendre Function［J］.Journal of Geodesy and Geodynamics，2009，29（2）：126－130.（王建強(qiáng)，趙國強(qiáng)，朱廣彬.常用超高階次締合勒讓德函數(shù)計(jì)算方法對比分析［J］.大地測量與地球動(dòng)力學(xué)，2009，29（2）：126－130.）

［8］ HEISKANEN W A，MOTITZ H.Physical Geodesy［M］.San Francisco：Freeman，1967.

［9］ COLOMBO C.Numerical Methods for Harmonic Analysis on the Sphere：Rep 310［R］.Columbus：Ohio State University，1981.

［10］ OVERTON M L.Numerical Computing with IEEE Floating Point Arithmetic［M］.Philadelphia：SIAM，2001.

［11］ KOOP R，SPTELPSTRA D.On the Computation of the Gravitational and Its First and Second Order Derivatives［J］.Manuscripta Geodaetica，1989（14）：373－382.

［12］ GLEASON D M.Partial Sums of Legendre Series via Clenshaw Summation［J］.Manuscr Geod，1985（10）：115－130.