常國賓,歐陽永忠,李勝全,金際航,李 科
(海軍海洋測繪研究所,天津300061)
等量緯度是地圖投影理論中的一種重要輔助緯度,它是大地緯度的函數,在地圖制圖、大地測量和地球物理等領域應用十分廣泛[1-6]。一般由大地緯度求解等量緯度的過程稱為等量緯度的正解問題,反之,則為反解問題。等量緯度的反解函數是復雜的隱函數,求解方法主要有迭代法和直接法,直接法便于采用和理論分析。在公開發表的文獻中,等量緯度反解問題的直接求解方法大致分為4種,分別為等量緯度的麥克勞林級數展開[6]、偏心率的麥克勞林級數展開[7-8]、拉格朗日級數展開[5,9]、埃爾米特插值法[10-13],其中,用埃爾米特插值法得到的公式其精度為最高。泰勒級數是麥克勞林級數的廣義形式,麥克勞林級數是一種在初值為0處的特殊泰勒級數,當要處理的值距離0較遠時,麥克勞林級數存在收斂較慢的問題,而泰勒級數可以通過合理地選擇展開初值避免此問題。鑒于此,本文嘗試利用泰勒級數展開法求解等量緯度的反解問題,其中展開初值的選取方法為:在假設偏心率為0的情況下,由等量緯度解析得到球心緯度作為大地緯度的初值,將球心緯度代入等量緯度的正解公式得到等量緯度初值。將等量緯度反解隱函數在等量緯度和大地緯度的初值處進行泰勒級數展開就可以得到等量緯度的反解公式,該公式為等量緯度和偏心率的函數。整個推導過程由強大的Mathematica計算機代數系統[10]完成,推導結果由計算機存儲,無需人工推導和記憶,且可以根據精度和計算成本靈活選擇泰勒級數展開的階數。試算結果表明,該方法具有明顯的精度優勢。
等量緯度與大地緯度之間的微分關系為

式中,q為等量緯度;B為大地緯度;e為第一偏心率。以B為積分變量,對上式進行積分可以得到

特別的,當e為0時,大地緯度B變為球面緯度φ,則有

因此有

對上式進行處理,即可以得到φ與B的關系式,一般為如下形式

式(3)和式(5)即為等量緯度正解展開公式,文獻[10]利用Mathematica計算機代數系統對上述系數進行了推導,糾正之前人工推導的系數中存在的錯誤。
假設等量緯度反解函數為如下形式

則有

華棠對式(6)在q=0(對應的,B=0)處進行泰勒級數展開(即麥克勞林級數展開),并利用式(7)所示的關系,得到了實用的等量緯度反解公式。
邊少鋒等從新的角度對等量緯度反解公式進行了研究,在文獻[10]中,首先根據三角級數回求公式確定式(6)的具體形式為

然后利用埃爾米特插值法確定式(9)中的各系數。
試算結果表明邊少鋒等推求的公式其精度優于前人得到的公式。
麥克勞林級數是泰勒級數在初值為0處的特殊形式。如果所處理的值距離0較遠,則麥克勞林級數收斂速度較慢,甚至可能發散。而如果合理地選擇初值,泰勒級數則具有較快的收斂速度,在相同的階數下,具有更高的函數近似精度。因此,相比傳統方法中對式(6)進行麥克勞林級數展開,在合適的初值處對式(6)進行泰勒級數展開,理論上可以得到更優的結果。
對一函數進行泰勒展開,首先需要確定展開初值,在等量緯度反解問題中,也就是確定B和q的初值B0和q0。
令e=0,則根據式(1)可以得到

或者

將大地緯度初值代入式(5),然后再代入式(3)則可以得到B0對應的等量緯度初值q0,即

求式(6)在B0和q0處的各階導數值

對式(6)在B=B0、q=q0處進行N階泰勒級數展開得

式(11)~式(15)即為由q求解B的過程,整個推導過程由Mathematica計算機代數系統自動完成,推導結果由計算機存儲,無需人工推導和記憶。公式為偏心率和等量緯度的函數,可以根據不同的參考橢球模型代入不同的偏心率,此外,還可以根據對精度和計算成本的不同要求選擇不同的階數N。
為了驗證本文方法的有效性,對等量緯度正、反解試算,其中反解過程中分別采用文獻[10,12]的方法和本文的方法,以比較兩種方法精度。驗證過程為:給定一大地緯度B,將B代入式(5),然后再代入式(3),得到q,然后分別將q代入式(8)、式(9)以及式(11)~式(15)得到 B1和B2,分別比較B1和B的差、B2和B的差。對幾種不同參考橢球的橢球參數進行了試算,得到的結論是相同的,簡便起見,只列出 WGS-84橢球基本常數 e=0.006 694 379 990 14的試算結果,如表1所示,其他參考橢球的試算結果不再贅述。

表1 WGS-84橢球等量緯度反解數據驗證
由表1的試算結果可以發現,本文方法在階數取到4時,其精度就已經全面超越傳統方法的精度,當階數取到5時,本文方法的精度幾乎在所有取值處都超過傳統方法一個數量級以上。由階數取6和階數取7的試算結果對比可以發現,當階數取到6以上時,其精度達到穩定,計算機的舍入誤差超過泰勒級數的截斷誤差,成為主要部分,精度不再隨階數的增加而提高。
本文采用泰勒級數展開法對等量緯度反解問題進行了新的研究,通過合理選擇展開初值,本文方法可以得到更為優越的反解精度。本文方法具有如下幾個方面的優勢:一是推導過程由Mathematica計算機代數系統完成,推導結果由計算機存儲,無需人工推導和記憶;二是推導的公式是偏心率的函數,可以靈活地應用于不同的參考橢球模型;三是在選定偏心率時,公式是等量緯度的函數,可以直接將等量緯度代入,無需迭代;四是精度明顯優于傳統方法;五是可以根據精度和計算成本的要求,靈活選擇泰勒級數展開的階數。本文方法還可以用于等面積緯度、等距離緯度等其他輔助緯度的反解問題,并可用于地圖投影等領域。
[1]SNYDER J.Map Projections-A Working Manual[M].Washington DC:Government Printing Office,1987.
[2]YANG Q,SNYDER J,TOBLER W.Map Projection Transformation:Principles and Applications[M].Florida:CRC Press,1999.
[3]BOWRING B.The Direct and Inverse Problems for Short Geodesic Lines on the Ellipsoid[J].Surveying and Mapping,1981,41(2):135-141.
[4]熊介.橢球大地測量學[M].北京:解放軍出版社,1988.
[5]楊啟和.地圖投影變換原理與方法[M].北京:解放軍出版社,1989.
[6]華棠.海圖數學基礎[M].北京:海潮出版社,1985.
[7]邊少鋒,紀兵.等距離緯度等量緯度和等面積緯度展開式[J].測繪學報,2007,36(2):218-223.
[8]BIAN Shao Feng,CHEN Yong Bing.Solving an Inverse Problem of a Meridian Arc in terms of Computer Algebra System[J].Journal of Surveying Engineering,2006,132(1):153-155.
[9]王瑞,李厚樸.輔助緯度反解公式的拉格朗日級數法推演[J].海洋測繪,2008,28(3):18-23.
[10]邊少鋒,許江寧.計算機代數系統與大地測量數學分析[M].北京:國防工業出版社,2004.
[11]紀兵,邊少鋒.大地主題問題的非迭代新解[J].測繪學報,2007,36(3):269-273.
[12]李厚樸,邊少鋒.等量緯度展開式的新解法[J].海洋測繪,2007,27(4):6-10.
[13]李厚樸,邊少鋒.輔助緯度反解公式的埃爾米特插值法新解[J].武漢大學學報:信息科學版,2008,33(6):623-626.