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

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

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

1 等量緯度展開式

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

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

(1)

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

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

(2)

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

(3)

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

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

(4)

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

(5)

部分系數(shù)為

(6)

令

(7)

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

(8)

通過逐次微分運算，可以建立遞推式

(9)

進一步可得系數(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)

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

2 等角緯度展開式

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

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

(20)

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

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

(21)

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

(22)

或者

(23)

特別地，有

(24)

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

(25)

(26)

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

2.2 等角緯度的正解展開式

根據(jù)文獻[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)

取當且僅當m為奇數(shù)且l=(m+1)/2時δ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)

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

(35)

2.3 等角緯度的反解展開式

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

(36)

因此

(37)

式中

(38)

或者

(39)

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

(40)

3 等距離緯度展開式

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

(41)

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

(42)

(43)

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

(44)

(45)

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

(46)

(47)

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

(48)

式中，系數(shù)

b0=B

(49)

(50)

或者

(51)

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

(52)

特別地

(53)

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

(54)

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

(55)

(56)

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

(57)

4 算例分析

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

(58)

(59)

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

表2 改進法求解等量緯度反解、等角緯度反解和等距離緯度正解展開Mathematica程序運行時間

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

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

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

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

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

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