中國大地坐標系地球旋轉(zhuǎn)橢球子午線弧長的嚴密精確解

田紅亮 朱大林 秦紅玲

（三峽大學 機械與材料學院，湖北 宜昌 443002）

1 橢球子午線弧長經(jīng)典解的3個缺陷

從赤道開始到任意大地緯度B的橢球子午線弧長為

式中，Φ為大地緯度.

經(jīng)典解［1]，首先將子午線曲率半徑M＝a（1－e2）×（1－e2sin2Φ）－1.5按二項式定理展開級數(shù)，取至8次項；然后將展開的8次項代入式（1），得

式中，5個系數(shù)為

但經(jīng)典解式（2）存在3個缺陷：①在給定精度的條件下，不能事先確定子午線曲率半徑M應展開至少幾項級數(shù)，且并非取級數(shù)項越多越好（第4節(jié)將說明）；②沒有從高等數(shù)學的嚴格意義上證明該級數(shù)解是否收斂、一致收斂、絕對收斂或發(fā)散，沒有給出誤差估計；③計算精度低，在表達全球尺度數(shù)據(jù)時，數(shù)據(jù)質(zhì)量難以保證，進而影響其實用化.

程鵬飛和文漢江等［2]根據(jù)文獻［1]的級數(shù)展開方法，直接給出由赤道到北極子午線弧長（在下文式（4）中，令B＝π／2）的近似級數(shù)解式（21），但沒有給出該近似公式的具體推導過程，且該近似公式只適用于大地緯度B＝π／2這1個點，不能求解B∈（0，π／2）更大區(qū)間的子午線弧長.為在動態(tài)地球的客觀環(huán)境中創(chuàng)建數(shù)字中國，中國政府決定從2008年7月1日起正式啟用中國大地坐標系2000（CGCS 2000）.中國若由西安80坐標系更換為CGCS 2000，則中國境內(nèi)地面奌大地緯度的變動區(qū)間為－1.6″～＋0.7″，因此對理論計算精度的要求越來越高［3－4].

過家春和趙秀俠等［5]提出基于第二類橢圓積分的子午線弧長公式，但該文獻［5]存在7個公式錯誤：第1個錯誤是原文式（1）

應為第2個錯誤是原文第95頁左邊第11～12行a（1－應為選用變量代換t＝sinφ，原文式（5）可化為E（φ，k）＝；第3個錯誤是原文式（6）E（x，k）應為E（x，k）第4個錯誤是原文式（9）應為第5個錯誤是原文式（16）

應為

第6個錯誤是原文式（17）

應為

第7個錯誤是原文圖2的標注E（e，φ）應為E（φ，e）.

本文在分析文獻［5]并糾正其錯誤的基礎(chǔ)上，根據(jù)定積分的換元法，推導地球旋轉(zhuǎn)橢球子午線弧長的精確解析解，且此解析解無中間變量.

2 橢球子午線弧長的絕對無偏差解析顯式解

式中，c為極點處的子午線曲率半徑；e′為橢圓的第二偏心率，e′≥0，且

式中，a、b分別為橢圓的長半軸和短半軸；e為橢圓的第一偏心率，0≤e≤1；f為橢圓的扁率.

將式（7）代入式（6），得

將式（5）、式（9）代入式（4），得

根據(jù)定積分的換元法，選擇以下變量代換

式（11）的微分［6]為

由式（11）可得歸化緯度為

通過變量代換式（11），式（10）可改寫為

將式（7）代入式（14），得

根據(jù)定積分的換元法，選擇以下變量代換

可將式（15）改寫為

根據(jù)定積分對于積分區(qū)間具有可加性，式（17）可改寫為

式中，E（e）為第二類完全橢圓積分［7]，E（sinφ，e）為勒讓德第二類橢圓積分，且

3 3種橢球子午線弧長的實例計算及比較

在式（4）中，令B＝π／2可得由赤道到北極子午線弧長的近似解［2]為

表1給出了CGCS 2000、GRS 80及 WGS 84［8]所采用的橢球常數(shù)［2，9].根據(jù)表1數(shù)據(jù)、式（18）和式（21），可計算橢球的子午線弧長，見表2，其中文獻［2]解指文獻［2]提供的數(shù)據(jù)；文獻［2]校正解指本文按照式（21）計算的數(shù)據(jù)；校正解誤差指文獻［2]校正解與本文解的絕對誤差.根據(jù)校正解誤差，可見文獻［2]解是合理的.值得指出，表2只驗證了橢圓第一偏心率取6 378 137＝0.081 819 191 042 815這一個值時較合理.圖1為橢球子午線弧長隨大地緯度的變化情況，圖中近似3條直線幾乎重合，說明根據(jù)式（18）計算的數(shù)字解有意義，分析表2的數(shù)據(jù)可看出3種橢球參數(shù)（表1的長、短半軸）微小變化帶來的子午線弧長變化.

表1 CGCS 2000橢球、GRS 80及WGS 84橢球基本常數(shù)比較

表2 CGCS 2000橢球與GRS 80和WGS 84橢球推導的子午線弧長比較

圖1 橢球子午線弧長隨大地緯度的變化

為了在e∈［0，1]更大區(qū)間內(nèi)驗證式（21）是否完全合理，下面將討論此問題.

4 文獻［2]近似級數(shù)解的缺陷及可能的補救方法

取CGCS 2000、GRS 80及WGS 84橢圓公共長半軸a＝6 378 137為例，橢球的子午線弧長如圖2（a）（e∈［0，0.96]）、圖2（b）（e∈［0，0.8]）所示.可見，當e∈［0，0.64]較小時，文獻［1－2]近似級數(shù)解與本文嚴密精確解都吻合得特別好，但當e∈［0.64，0.8]較大時，文獻［1－2]近似級數(shù)解的誤差都較大，當e∈［0.8，1]更大時，文獻［2]近似級數(shù)解的誤差特別大（未給出）.而且，根據(jù)圖2（a），本文解與文獻［1]解隨橢圓第一偏心率增大而減小；但根據(jù)圖2（b），文獻［2]解隨橢圓第一偏心率增大而增大.

圖2 橢球子午線弧長隨橢圓第一偏心率的變化

導致文獻［2]近似級數(shù)解誤差大的原因是：根據(jù)式（9）知，因變量e′是自變量e的單調(diào)增加的函數(shù)，因此當自變量e較小時，e′較小，則文獻［2]近似級數(shù)解式（21）誤差就小，一種極限情況是當e→0時，e′→0，根據(jù)式（18）知，根據(jù)式（21）知，兩者結(jié)果一樣；但當e較大時，e′較大（一種極限情況是當e→1時，e′→＋∞），則文獻［2]近似級數(shù)解式（21）誤差就大，且取級數(shù)項越多，誤差越大，這證明了第1節(jié)所說的“并非取級數(shù)項越多越好”的觀點.

使文獻［2]近似級數(shù)解式（21）誤差較小的一個可能補救方法是：當e較大時，以e為自變量展開級數(shù)，不以e′為自變量展開級數(shù).

5 結(jié) 論

根據(jù)定積分的換元法，推導了地球旋轉(zhuǎn)橢球子午線弧長的解析解，為子午線弧長計算的實用化和誤差控制提供了理論依據(jù).本文的理論公式嚴密精確解，計算工作量小；而級數(shù)展開［10－12]、數(shù)值積分［13]、Hermite插值［14]、Horner求和技術(shù)［15]（該算法至少可以計算到3600完全階次的球諧級數(shù)式，可把高階次重力場模型向更高階次擴展）都是近似解，計算工作量大，且本文第4節(jié)表明“并非取級數(shù)項越多越好”.本文的求解思路有可能較好解決與橢球表面積相等的球半徑R2［2]、橢球面梯形圖幅面積等［1]的理論計算問題.

在本文工作的基礎(chǔ)上，再做一些理論工作（例如第一類完全橢圓積分、勒讓德第一類橢圓積分等），這些理論公式可供實際工程直接使用［16－21].文獻［16－20]研究兩球體單峰接觸問題，這種將復雜工程問題過于簡單化導致的結(jié)果是：機械結(jié)構(gòu)的理論振型與實驗振型有時不一致，機械結(jié)構(gòu)的固有頻率與實驗固有頻率之間的相對誤差較大，例如文獻［18]預測的相對誤差在－19.2%～16.8%之間.事實上，在實際工程中，接觸壓力的橢球—Hertz分布將在兩個物體中產(chǎn)生與所提出的橢圓接觸面相協(xié)調(diào)的彈性位移［21]，相應地大量存在以第二類完全橢圓積分式（19）、勒讓德第二類橢圓積分式（20），以及本文尚未涉及的第一類完全橢圓積分勒讓德第一類橢圓積分這4類積分為中間自變量，衍生出許多顯函數(shù)，可能會使螺栓結(jié)合部［22]的理論預測結(jié)果與實驗數(shù)據(jù)更接近.

必須指出，文獻［23]中的式（50）應為

［1]孔祥元，郭際明，劉宗泉.大地測量學基礎(chǔ)［M].武昌：武漢大學出版社，2002：67－73.

［2]程鵬飛，文漢江，成英燕，等.2000國家大地坐標系橢球參數(shù)與GRS 80和WGS 84的比較［J].測繪學報，2009，38（3）：189－194.

［3]廖瑞祥，王 剛，鄒良超.基于GIS的三峽庫區(qū)滑坡空間數(shù)據(jù)庫系統(tǒng)設(shè)計［J].三峽大學學報：自然科學版，2011，33（1）：24－27.

［4]王海棟 ，柴洪洲，王 敏.多波束測深數(shù)據(jù)的抗差Kriging擬合［J].測繪學報，2011，40（2）：238－242，248.

［5]過家春，趙秀俠，徐 麗，等.基于第二類橢圓積分的子午線弧長公式變換及解算［J].大地測量與地球動力學，2011，31（4）：94－98.

［6]同濟大學數(shù)學系.高等數(shù)學（上冊）［M].6版.北京：高等教育出版社，2011：116－117.

［7]Mindlin Raymond David.Compliance of Elastic Bodies in Contact［J].ASME Journal of Applied Mechanics，1949，16（3）：259－268.

［8]白建軍，孫文彬，趙學勝.基于QTM的 WGS－84橢球面層次剖分及其特點分析［J].測繪學報，2011，40（2）：243－248.

［9]陳俊勇.中國現(xiàn)代大地基準——中國大地坐標系統(tǒng)2000（CGCS 2000）及其框架［J].測繪學報，2008，37（3）：269－271.

［10]牛卓立.以空間直角坐標為參數(shù)的子午線弧長計算公式［J].測繪通報，2001（11）：14－15.

［11]邊少鋒，紀 兵.等距離緯度等量緯度和等面積緯度展開式［J].測繪學報，2007，36（2）：218－223.

［12]李厚樸，邊少鋒，陳良友.等面積緯度函數(shù)和等量緯度變換的直接解算公式［J].武漢大學學報·信息科學版，2011，36（7）：843－846.

［13]劉修善.計算子午線弧長的數(shù)值積分法［J].測繪通報，2006，（5）：4－6.

［14]李厚樸，邊少鋒.輔助緯度反解公式的Hermite插值法新解［J].武漢大學學報·信息科學版，2008，33（6）：623－626.

［15]劉纘武，劉世晗，黃 歐.超高階次勒讓德函數(shù)遞推計算中的壓縮因子和Horner求和技術(shù)［J].測繪學報，2011，40（4）：454－458.

［16]田紅亮，朱大林，秦紅玲，等.結(jié)合部法向載荷解析解修正與定量實驗驗證［J].農(nóng)業(yè)機械學報，2011，42（9）：213－218.

［17]田紅亮，朱大林，秦紅玲.固定接觸界面法向靜彈性剛度［J].應用力學學報，2011，28（3）：318－322.

［18]田紅亮，方子帆，朱大林，等.固定接觸界面切向靜彈性剛度問題研究［J].應用力學學報，2011，28（5）：458－464.

［19]田紅亮，朱大林，秦紅玲.結(jié)合面靜摩擦因數(shù)分形模型的建立與仿真［J].應用力學學報，2011，28（2）：158－162.

［20]Tian Hongliang，Li Bin，Liu Hongqi，et al.A New Method of Virtual Material Hypothesis－Based Dynamic Modeling on Fixed Joint Interface in Machine Tools［J].ELSEVIER International Journal of Machine Tools ＆Manufacture，2011，51（3）：239－249.

［21]田紅亮，朱大林，方子帆，等.赫茲接觸129年［J].三峽大學學報：自然科學版，2011，33（6）：61－71.

［22]田紅亮，趙春華，朱大林，等.金屬材料結(jié)合部法切向剛度修正與實驗驗證［J].農(nóng)業(yè)機械學報，2012，43（6）：207－214.

［23]田紅亮，朱大林，秦紅玲.地球旋轉(zhuǎn)橢球等面積緯度函數(shù)與等量緯度的相互近似變換［J].三峽大學學報：自然科學版，2012，34（2）：71－75.