陳 成，金立新，邊少鋒，李松林

1. 海軍工程大學(xué)導(dǎo)航工程系，湖北 武漢 430033； 2. 中鐵第一勘察設(shè)計(jì)院集團(tuán)有限公司，陜西 西安 710043； 3. 甘肅鐵道綜合工程勘察院有限公司，甘肅 蘭州 730000

在大地測(cè)量和制圖學(xué)中，除了大地緯度之外，還有6種常用的輔助緯度，它們都有著廣泛的應(yīng)用[1-3]。如等量緯度常用于墨卡托投影，地心緯度和歸化緯度常用于橢球的幾何學(xué)，等角緯度、等距離緯度和等面積緯度常用于雙重投影[4-6]。在進(jìn)行投影轉(zhuǎn)換等分析計(jì)算時(shí)，經(jīng)常用到它們與大地緯度間的解析展開式，眾多國內(nèi)外學(xué)者對(duì)此展開了研究[7]。地心緯度、歸化緯度和大地緯度間正反解函數(shù)關(guān)系式較簡(jiǎn)潔，文獻(xiàn)[8]利用拉格朗日共軛級(jí)數(shù)法嚴(yán)格得到了相應(yīng)的無窮級(jí)數(shù)展開。而其余4種輔助緯度與大地緯度間關(guān)系式較復(fù)雜，除等量緯度、等角緯度與大地緯度間正解具有明顯的函數(shù)關(guān)系式外，一般表現(xiàn)為第一偏心率的冪級(jí)數(shù)形式，或者大地緯度的三角級(jí)數(shù)形式。由于直接推導(dǎo)的計(jì)算量較大，文獻(xiàn)[9—10]也只分別得到了展開至e′6、e8階等量緯度關(guān)于大地緯度的反解展開式。文獻(xiàn)[11]首次系統(tǒng)研究了輔助緯度與大地緯度間的展開式，并給出了至e8階的正反展開式。由于計(jì)算機(jī)代數(shù)系統(tǒng)可以進(jìn)行符號(hào)運(yùn)算，一些較復(fù)雜的人工推導(dǎo)可以交由計(jì)算機(jī)完成。因此，借助Mathematica、Maxima等計(jì)算機(jī)代數(shù)軟件，文獻(xiàn)[12—13]及其他文獻(xiàn)(http:∥geographiclib.sourceforge.net/html/auxlat.html;https:∥www.academia.edu/7580468/Funciones_de_Latitud)得到了更高階的等角緯度、等距離緯度和等面積緯度關(guān)于大地緯度正反解展開式。然而，除了地心緯度和歸化緯度，其他幾種輔助緯度與大地緯度間的無窮展開仍然沒有給出，上述文獻(xiàn)致力于直接推導(dǎo)求解有限階的輔助緯度展開式，也沒有給出展開式系數(shù)的一般公式。雖然這類無窮展開表現(xiàn)為冪級(jí)數(shù)和三角函數(shù)組合的形式，但利用傅里葉級(jí)數(shù)方法仍然很難得到具體的無窮展開式。另外，由于等面積緯度與大地緯度間的關(guān)系式存在多層函數(shù)嵌套，也很難得到等面積緯度與大地緯度間的無窮展開式。由此，本文擬通過泰勒展開定理和拉格朗日反演定理，推導(dǎo)等量緯度、等角緯度和等距離緯度與大地緯度間的無窮展開式，并取CGCS2000參考橢球，對(duì)輔助緯度正反解展開式進(jìn)行算例分析。

1 等量緯度展開式

由文獻(xiàn)[4,14]，等量緯度q與大地緯度B的關(guān)系式為

q=gd-1(B)-earctanh(esinB)

(1)

式中，e是橢球第一偏心率，為一個(gè)小參數(shù)，gd-1(x)=arctanh(sinx)，為古德曼函數(shù)的反函數(shù)[15-16]。當(dāng)大地緯度趨近于極點(diǎn)時(shí)，等量緯度趨向于無窮大。在進(jìn)行數(shù)值計(jì)算時(shí)，為了避免舍入誤差，gd-1(x)通常修正為gd-1(x)=arcsinh(tanx)。

根據(jù)反雙曲函數(shù)的級(jí)數(shù)展開，等量緯度可展開成無窮級(jí)數(shù)

(2)

等量緯度與大地緯度的反解可通過等角緯度與大地緯度的反解展開式和等量緯度與等角緯度的閉合關(guān)系式而得到，但是等量緯度與大地緯度的反解也有直接關(guān)系式。為求得等量緯度與大地緯度的反解展開，應(yīng)用拉格朗日反演定理[17]，有

(3)

式中，gd(x)=arcsin(tanhx)，為古德曼函數(shù)[15-16]。同樣，為了避免舍入誤差，取計(jì)算式gd(x)=arctan(sinhx)。

現(xiàn)在來求得一個(gè)有用的級(jí)數(shù)恒等式，利用組合函數(shù)和冪級(jí)數(shù)的性質(zhì)[16]，可得

(4)

式中，系數(shù)ζmk用遞推形式表示為

(5)

部分系數(shù)為

(6)

令

(7)

將式(7)和式(4)代入式(3)得

(8)

通過逐次微分運(yùn)算，可以建立遞推式

(9)

進(jìn)一步可得系數(shù)遞推式

(10)

(11)

(12)

式(10)中的系數(shù)遞推關(guān)系可用圖解表示，如圖1所示。

圖遞推關(guān)系Fig.1 Recursion relations for

(13)

因此，有

(14)

或者

(15)

式中，系數(shù)

(16)

[ηmk]=

(17)

(18)

式中，e′為橢球第二偏心率。

(19)

文獻(xiàn)[9—10]分別給出了到e′6、e8階等量緯度與大地緯度反解展開式的精度，文獻(xiàn)[10]還討論了反解展開式的精度。文獻(xiàn)[18—20]也分別用數(shù)值方法、解析方法對(duì)大地緯度反解等量緯度作了一定的研究。利用第二偏心率與第一偏心率的關(guān)系式和恒等式sech2q=1-tanh2q，可驗(yàn)證文獻(xiàn)[9—10]的結(jié)果分別與本文展開至e6、e8的結(jié)果一致。

2 等角緯度展開式

2.1 古德曼函數(shù)的泰勒展開

根據(jù)泰勒展開定理，古德曼函數(shù)在有一增量Δ時(shí)，可表達(dá)為

(20)

式中，gd(m)(x)為古德曼函數(shù)對(duì)自變量的m階導(dǎo)數(shù)。

利用導(dǎo)數(shù)遞推公式和數(shù)學(xué)歸納法，容易得到

(21)

式中，系數(shù)G11=-H11=1，當(dāng)k=0或k>m時(shí)Gmk=Hmk=0，其他情況可由下列遞推公式計(jì)算

(22)

或者

(23)

特別地，有

(24)

以及一些低階系數(shù)

(25)

(26)

反正弦函數(shù)也有類似的泰勒展開，可用于求解等面積緯度的展開式。但由于等面積緯度與大地緯度的關(guān)系式存在多層函數(shù)嵌套，很難求得一個(gè)簡(jiǎn)潔的無窮展開式，需借助Mathematica或者M(jìn)axima軟件來求解一定階的級(jí)數(shù)展開式。

2.2 等角緯度的正解展開式

根據(jù)文獻(xiàn)[4,14]，等角緯度φ與大地緯度的關(guān)系式為

φ=gdq=gdgd-1B-earctanh(esinB)

(27)

將式(2)、式(4)代入式(27)，得

(sinB)2m-kgd(k)(gd-1B)

(28)

不難得到

(29)

(30)

sin 2(p-r)B+sin 2(p+r+1)B

(31)

取當(dāng)且僅當(dāng)m為奇數(shù)且l=(m+1)/2時(shí)δml=0，其余情況為δml=1，以及

(32)

利用恒等式sinh(gd-1B)=tanB和cosh(gd-1B)=secB，連同式(20)、式(28)—式(31)，可得

(33)

式中，系數(shù)

δmlζ2l,m-2lHls-ζ2l-1,m-2l+1Gls

(34)

式中，相對(duì)于式(31)系數(shù)Cpr中的整數(shù)m、k，此處應(yīng)用m+s-l-1、l-s替代；展開到e8階的系數(shù)參考文獻(xiàn)[11]，展開到e10階的系數(shù)和精度分析參考文獻(xiàn)[13]，本文計(jì)算結(jié)果均與其一致，并補(bǔ)充e12階的系數(shù)

(35)

2.3 等角緯度的反解展開式

將式(27)代入式(14)，并利用恒等式tanhq=sinφ和sechq=cosφ，有

(36)

因此

(37)

式中

(38)

或者

(39)

展開到e10階的系數(shù)和精度分析仍然可以參考文獻(xiàn)[13]，本文計(jì)算結(jié)果與其一致，并補(bǔ)充e12階的系數(shù)

(40)

3 等距離緯度展開式

根據(jù)文獻(xiàn)[4,21]，等距離緯度ψ與大地緯度的關(guān)系式為

(41)

式中，X為子午線弧長，R為等距離半徑，分別為(取橢球常半軸為單位長度)

(42)

(43)

根據(jù)冪級(jí)數(shù)展開或者傅里葉級(jí)數(shù)展開方法[21-22]，可得

(44)

(45)

式中，系數(shù)x0=B,r0=1，

(46)

(47)

從式(44)、式(46)可以看出，子午線弧長X展開式關(guān)于大地緯度線性項(xiàng)的系數(shù)就是等距離半徑R。因此，根據(jù)級(jí)數(shù)除法公式[16]，有

(48)

式中，系數(shù)

b0=B

(49)

(50)

或者

(51)

由于rk=xk0(k≥1)，式(51)中行列式第一列rkx0-r0xk關(guān)于B的線性項(xiàng)rkx0-r0xk0B為0，對(duì)bm的計(jì)算沒有影響，這也是顯然的。因此，在計(jì)算bmk時(shí)，可以去除所有關(guān)于B線性項(xiàng)和正弦項(xiàng)sin 2lB(l≠k)，再通過行列式的逐次運(yùn)算，可進(jìn)一步得到

(52)

特別地

(53)

展開到e10階的系數(shù)和精度分析參考文獻(xiàn)[13]，本文計(jì)算結(jié)果與其一致，并補(bǔ)充e12階的系數(shù)

(54)

對(duì)于等距離緯度關(guān)于大地緯度的反解展開式，可以根據(jù)正解公式，運(yùn)用拉格朗日反演方法，或者直接建立常微分方程

(55)

(56)

等距離緯度反解展開式的計(jì)算量很大，也很難得到比較簡(jiǎn)潔的遞推公式，必須借助Mathematica或者M(jìn)axima軟件計(jì)算得到一定階的展開式。展開到e10階的系數(shù)和精度分析可參考文獻(xiàn)[13]，本文補(bǔ)充e12階的系數(shù)

(57)

4 算例分析

為了進(jìn)一步驗(yàn)證本文輔助緯度與大地緯度間無窮展開式的正確性，可將本文方法(遞推法)得到的高階展開式與計(jì)算機(jī)代數(shù)系統(tǒng)直接推導(dǎo)(直接法)得到的結(jié)果進(jìn)行比較。其中，在利用計(jì)算機(jī)代數(shù)系統(tǒng)直接推導(dǎo)正解展開式時(shí)是對(duì)原函數(shù)進(jìn)行泰勒展開，而在反解時(shí)直接應(yīng)用拉格朗日反演定理。顯然，本文求解展開式系數(shù)的遞推方法也可以利用計(jì)算機(jī)代數(shù)系統(tǒng)編程實(shí)現(xiàn)。為了節(jié)省程序運(yùn)行時(shí)間，應(yīng)用拉格朗日反演定理時(shí)應(yīng)先把三角級(jí)數(shù)乘積化簡(jiǎn)成倍角形式再進(jìn)行求導(dǎo)，而在遞推求解等距離緯度正解展開式時(shí)，應(yīng)用如下遞推公式替代式(46)、式(47)

(58)

(59)

表1 本文方法(遞推法)與直接法求解輔助緯度展開式Mathematica程序運(yùn)行時(shí)間

表2 改進(jìn)法求解等量緯度反解、等角緯度反解和等距離緯度正解展開Mathematica程序運(yùn)行時(shí)間

由表1、表2可以看出，求解等量緯度反解時(shí)，本文方法一般比直接法計(jì)算用時(shí)短，但超過e40階時(shí)直接法可能更快一些，這是因?yàn)榈攘烤暥冉馕鍪?1)并不復(fù)雜，但在改進(jìn)程序后，本文方法計(jì)算速度遠(yuǎn)遠(yuǎn)優(yōu)于直接法；求解等角緯度正解時(shí)，直接法計(jì)算用時(shí)遠(yuǎn)遠(yuǎn)大于本文方法，說明等角緯度進(jìn)行直接泰勒展開并不是一種高效的運(yùn)算，耗費(fèi)了大量時(shí)間；求解等角緯度反解時(shí)，取e0～e20階本文方法計(jì)算用時(shí)短，取更高階時(shí)直接法用時(shí)短，這是因?yàn)橹苯臃ㄖ苯永昧苏庹归_式系數(shù)，經(jīng)過改進(jìn)后，本文方法計(jì)算用時(shí)大大縮短，小于直接法；求解等距離緯度正解時(shí)，直接法計(jì)算用時(shí)大于本文方法，但直接法在改進(jìn)后用時(shí)縮短，也比本文方法更快，說明Mathematica內(nèi)部對(duì)級(jí)數(shù)除法運(yùn)算優(yōu)化地較好。綜合來說，進(jìn)行具體系數(shù)運(yùn)算時(shí)，本文方法不僅提供了一種遞推計(jì)算方法，也加快了運(yùn)算。但更重要的是，本文分析的是展開式及其系數(shù)的一般形式，這是直接法所無法比擬的。

圖2 等量緯度隨大地緯度B∈(0,90°)變化的曲線Fig.2 Curve of the isometric latitude varying with the geodetic latitude B∈(0,90°)

由圖2可看出，在B∈(0,90°)時(shí)，等量緯度隨大地緯度單調(diào)遞增，在接近極點(diǎn)時(shí)，曲線斜率非常大。理論上在極點(diǎn)處等量緯度為無窮大，但進(jìn)行數(shù)值計(jì)算時(shí)不可能處理無窮大量；16位和34位數(shù)字精度下所能計(jì)算的等量緯度最大值分別為qmax16=38.018 293 995 274 90，qmax34=78.115 872 564 187 653 490 898 757 927 628 82，可得其比值qmax34/qmax16≈2.05；因此，提高數(shù)字精度可有效增加等量緯度計(jì)算范圍，顯然，這也增加了數(shù)值計(jì)算精度。由圖3可以看出，由于舍入誤差影響，在古德曼反函數(shù)未修正時(shí)(gd-1B=arctanh(sinB))，等量緯度隨大地緯度變化曲線在趨近極點(diǎn)時(shí)成折線，在一定區(qū)間內(nèi)不再變化；等量緯度計(jì)算范圍減小，最大值僅為18.708 264 496 564 56，在B>89.999 999 396 289 99°時(shí)根本無法計(jì)算；而在古德曼反函數(shù)修正后，等量緯度函數(shù)曲線是正常、連續(xù)的；這些是基于16位數(shù)字精度運(yùn)算情況下的，對(duì)34位數(shù)字精度情況也有類似結(jié)果。

為了分析等量緯度反解精度，設(shè)有一大地緯度B，可由解析式(1)計(jì)算得到等量緯度，再利用反解式(14)得到另一大地緯度B′，因此B′-B就是反解理論值與實(shí)際值的誤差。等量緯度反解相對(duì)誤差(取對(duì)數(shù)log10(B′-B)/B)隨大地緯度變化的曲線如圖4所示。

由圖4可以看出，等量緯度反解相對(duì)誤差隨展開式階數(shù)逐漸減小，在16位數(shù)字精度下，取到e8階時(shí)反解式相對(duì)誤差最大絕對(duì)值約為1.34×10-11，取到e10階時(shí)約為9.04×10-14，取到e12階時(shí)約為1.14×10-15，已經(jīng)接近數(shù)字精度，取更高階時(shí)精度幾乎不再增加(Fortran的16位機(jī)器精度約為2.22×10-16)；而在34位數(shù)字精度下，取到e8、e10和e12階時(shí)反解式最大相對(duì)誤差約為1.34×10-11、9.00×10-14和6.03×10-16，比16位精度時(shí)略好，取到e28階時(shí)接近數(shù)字精度(Fortran的34位機(jī)器精度約為1.93×10-34)。另外，等量緯度正解展開式取到e8、e10和e12時(shí)最大相對(duì)誤差分別為5.49×10-13、3.34×10-15和5.44×10-16(16位精度時(shí))，比反解精度高，不過一般直接用定義的解析式直接計(jì)算等量緯度；其余輔助緯度與大地緯度之間展開式的相對(duì)誤差情況類似，本文不再贅述，文獻(xiàn)[4,13]也有一定論述。在大地測(cè)量和制圖學(xué)的實(shí)際應(yīng)用中，取到e8或e10階一般已經(jīng)滿足精度要求，而要求精度較高時(shí)，為了避免舍入誤差，一般取到e12階。

圖3 等量緯度隨大地緯度B∈(89.999 99°,90°)變化的曲線Fig.3 Curve of the isometric latitude varying with the geodetic latitude B∈(89.999 99°,90°)

圖4 等量緯度反解誤差隨大地緯度變化的曲線Fig.4 Relative error of the inverse solution of the isometric latitude

5 結(jié) 論

本文從無窮級(jí)數(shù)理論出發(fā)，詳細(xì)推導(dǎo)了旋轉(zhuǎn)橢球情況下等量緯度、等角緯度和等距離緯度關(guān)于大地緯度和參考橢球第一偏心率的無窮級(jí)數(shù)公式，主要表現(xiàn)為遞推形式。計(jì)算結(jié)果表明，展開至e6、e8階的等量緯度反解式與文獻(xiàn)[9—10]等結(jié)果是一致的，展開至e10階輔助緯度展開式也與文獻(xiàn)[13]結(jié)果一致；借助Mathametica進(jìn)一步檢驗(yàn)了e0～e40階展開式，并比較了本文方法與利用計(jì)算機(jī)代數(shù)系統(tǒng)直接推導(dǎo)展開式的程序運(yùn)行時(shí)間，不僅說明本文方法是正確的，也可以加快展開式的求解運(yùn)算；利用Fortran對(duì)輔助緯度正反解進(jìn)行了數(shù)值分析，說明增加數(shù)字精度可以增加等量緯度計(jì)算范圍，也略增了數(shù)值計(jì)算精度，若要避免16位、32位數(shù)字精度運(yùn)算的舍入誤差，展開式分別需要展開到e12階、e28階。本文嚴(yán)格推導(dǎo)的輔助緯度關(guān)于大地緯度的無窮展開公式，是一種一般形式，也是一種有效、快速的算法，對(duì)于完善制圖學(xué)的數(shù)學(xué)體系具有重要參考意義。