中国大地坐标系地球旋转椭球子午线弧长的严密精确解

中国大地坐标系地球旋转椭球子午线弧长的严密精确解
田红亮;朱大林;秦红玲
【摘要】由于椭球子午线弧长经典解存在3个缺陷,在大地测量学上会引入无法接受的误差.根据定积分的换元法,推导以大地纬度为自变量的地球旋转椭球子午线弧长精确解析解.按照2000中国大地坐标系、GRS 80和WGS 84所定义的2个基本椭圆常数(长半轴、短半轴),给出相应3种椭球子午线弧长.以椭圆第一偏心率为自变量,将解析解与文献[2]近似级数解做比较.研究结果表明,当椭圆第一偏心率较大时,文献[2]近似级数解的误差较大.在Matlab R2009b语言中,开发了2套命令文件.解析解为子午线弧长计算的实用化和误差控制提供了理论依据.%The classical solution for the ellipsoid meridian arc length, which results in the unacceptable warp in geodesy, has 3 limitations. According to transforming element method of fixed integral, the precise analytic expression for earth rotating ellipsoid meridian arc length was deduced by considering the geodetic latitude an independent variable. By virtue of 2 basic elliptic constants (larger semi axis and smaller semi axis) defined by China Geodetic Coordinate System 2000 (CGCS 2000), GRS 80 and WGS 84, three relevant ellipsoid meridian arc length were provided. Taking account of the elliptic first eccentricity as an independent variable, the analytic answer was compared with the approximate series solution of reference[2], The research results make clear that the windage of approximate series solution of reference [2] becomes larger as the elliptic first eccentricity increases. Two sets of command files were compiled in Matlab R2009b
language. The analytic solution can be the theoretical basis of its utility and error control in calculating the meridian arc length.
【期刊名称】《三峡大学学报（自然科学版）》
【年(卷),期】2012(034)004
【总页数】6页(P81-85,112)
【关键词】地球旋转椭球;子午线弧长;解析解;定积分的换元法
【作者】田红亮;朱大林;秦红玲
【作者单位】三峡大学机械与材料学院,湖北宜昌 443002;三峡大学机械与材料学院,湖北宜昌 443002;三峡大学机械与材料学院,湖北宜昌 443002
【正文语种】中文
【中图分类】P22
1 椭球子午线弧长经典解的3个缺陷
从赤道开始到任意大地纬度B的椭球子午线弧长为
式中，Φ为大地纬度.
经典解［1]，首先将子午线曲率半径 M＝a（1－e2）×（1－e2sin2Φ）－1.5按二项式定理展开级数，取至8次项；然后将展开的8次项代入式（1），得
式中，5个系数为
式中，5个递归系数满足
但经典解式（2）存在3个缺陷：①在给定精度的条件下，不能事先确定子午线曲率半径M应展开至少几项级数，且并非取级数项越多越好（第4节将说明）；②没有从高等数学的严格意义上证明该级数解是否收敛、一致收敛、绝对收敛或发散，没有给出误差估计；③计算精度低，在表达全球尺度数据时，数据质量难以保证，进而影响其实用化.
程鹏飞和文汉江等［2]根据文献［1]的级数展开方法，直接给出由赤道到北极子午线弧长（在下文式（4）中，令B＝π／2）的近似级数解式（21），但没有给出该近似公式的具体推导过程，且该近似公式只适用于大地纬度B＝π／2这1个点，
不能求解B∈（0，π／2）更大区间的子午线弧长.为在动态地球的客观环境中创建数字中国，中国政府决定从2008年7月1日起正式启用中国大地坐标系2000（CGCS 2000）.中国若由西安80坐标系更换为CGCS 2000，则中国境内地面奌大地纬度的变动区间为－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]并纠正其错误的基础上，根据定积分的换元法，推导地球旋转椭球子午线弧长的精确解析解，且此解析解无中间变量.
2 椭球子午线弧长的绝对无偏差解析显式解
由式（9）得且根据式（7）、式（5），可将式（1）改写为
式中，c为极点处的子午线曲率半径；e′为椭圆的第二偏心率，e′≥0，且
式中，a、b分别为椭圆的长半轴和短半轴；e为椭圆的第一偏心率，0≤e≤1；f 为椭圆的扁率.
将式（7）代入式（6），得
将式（5）、式（9）代入式（4），得
根据定积分的换元法，选择以下变量代换
式（11）的微分［6]为
由式（11）可得归化纬度为
通过变量代换式（11），式（10）可改写为
将式（7）代入式（14），得
根据定积分的换元法，选择以下变量代换
可将式（15）改写为
根据定积分对于积分区间具有可加性，式（17）可改写为
式中，E（e）为第二类完全椭圆积分［7]，E（si nφ，e）为勒让德第二类椭圆积分，且
3 3种椭球子午线弧长的实例计算及比较
在式（4）中，令B＝π／2可得由赤道到北极子午线弧长的近似解［2]为
表1给出了CGCS 2000、GRS 80及 WGS 84［8]所采用的椭球常数［2，9].根
据表1数据、式（18）和式（21），可计算椭球的子午线弧长，见表2，其中文
献［2]解指文献［2]提供的数据；文献［2]校正解指本文按照式（21）计算的数据；校正解误差指文献［2]校正解与本文解的绝对误差.根据校正解误差，可见文献［2]解是合理的.值得指出，表2只验证了椭圆第一偏心率取6 378 137＝0.081 819 191 042 815这一个值时较合理.图1为椭球子午线弧长随大地纬度的变化情况，
图中近似3条直线几乎重合，说明根据式（18）计算的数字解有意义，分析表2
的数据可看出3种椭球参数（表1的长、短半轴）微小变化带来的子午线弧长变化.
表1 CGCS 2000椭球、GRS 80及WGS 84椭球基本常数比较参数CGCS 2000 GRS 80 WGS 84长半轴a／m 6 378 137 6 378 137 6 378 137短半轴b／m 6 356 752.314 140 355 8 6 356 752.314 140 347 4 6 3 56 752.314 245 179 5 表2 CGCS 2000椭球与GRS 80和WGS 84椭球推导的子午线弧长比较大地纬度B／° CGCS 2000的子午线弧长X／m GRS 80的子午线弧长X／m WGS 84的子午线弧长X／m 0 0 0 0 5 552 885.451 040 19 552 885.451 040 20 552 885.451 058 36 10 1 105 854.833 198 45 1 105 854.833 198 45 1 105 854.833 234 37 15 1 658 989.589 347 67 1 658 989.589 347 67 1 658 989.589 400 55 20 2 212 366.254 102 98 2 212 366.254 102 97 2 212 366.254 171 63 25 2 766 054.169 063 11 2 766 054.169 063 10 2 766 054.169 146 03 30 3 320 113.397 845 02 3 320 113.397 845 02 3 320 113.397 940 38 35 3 874 592.901 589 18 3 874 592.901 589 17 3 874 592.901 694 94 40 4 429 529.030 236 59 4 429 529.030 236 58 4 429 529.030 35 051 45 4 984 944.377 858 00 4 984 944.377 857 99 4 984 944.377 977 74 50 5 540 847.041 560 97 5 540 847.041 560 96 5 540 847.041 684 15 55 6 097 230.312 999 93 6 097 230.312 999 92 6 097 230.313 124 18 60 6 654 072.819 367 44 6 654 072.819 367 43 6 654 072.819 490 51 65 7 211 339.117 188 21 7 211 339.117 188 20 7 211 339.117 308 01 70 7 768 980.727 655 52 7 768 980.727 655 51 7 768 980.727 770 19 75 8 326 937.587 172 34 8 326 937.587 172 34 8 326 937.587 280 35 80 8 885 139.871 836 76 8 885 139.871 836 75 8 885 139.871 936 87 85 9 443 510.140 574 85 9 443 510.140 574 84 94 43
510.140 666 25 90 10 001 965.729 230 46（本文解） 10 001 965.729 230 46（本文解） 10 001 965.729 312 72（本文解）10 001 965.729 140 73（文献［1]解） 10 001 965.729 140 73（文献［1]解） 10 001 965.729 222 99（文
献［1]解）10 001 965.729 323（文献［2]解） 10 001 965.729 323（文献［2]解） 10 001 965.729 405（文献［2]解）10 001 965.729 322 96（文献［2]校
正解） 10 001 965.729 322 95（文献［2]校正解） 10 001 965.729 405 22
（文献［2]校正解）文献［1]误差／m －0.000 089 73 －0.000 089 73 －0.000 089 73校正解误差／m 0.000 092 50 0.000 092 49 0.000 092 50
图1 椭球子午线弧长随大地纬度的变化
为了在e∈［0，1]更大区间内验证式（21）是否完全合理，下面将讨论此问题.
4 文献［2]近似级数解的缺陷及可能的补救方法
取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]近似级数解与本文严密精确解都吻
合得特别好，但当e∈［0.64，0.8]较大时，文献［1－2]近似级数解的误差都较大，当e∈［0.8，1]更大时，文献［2]近似级数解的误差特别大（未给出）.而且，根据图2（a），本文解与文献［1]解随椭圆第一偏心率增大而减小；但根据图2（b），文献［2]解随椭圆第一偏心率增大而增大.
图2 椭球子午线弧长随椭圆第一偏心率的变化
导致文献［2]近似级数解误差大的原因是：根据式（9）知，因变量e′是自变量e
的单调增加的函数，因此当自变量e较小时，e′较小，则文献［2]近似级数解式（21）误差就小，一种极限情况是当e→0时，e′→0，根据式（18）知，根据式（21）知，两者结果一样；但当e较大时，e′较大（一种极限情况是当e→1时，
e′→＋∞），则文献［2]近似级数解式（21）误差就大，且取级数项越多，误差越大，这证明了第1节所说的“并非取级数项越多越好”的观点.
使文献［2]近似级数解式（21）误差较小的一个可能补救方法是：当e较大时，
以e为自变量展开级数，不以e′为自变量展开级数.
5 结论
根据定积分的换元法，推导了地球旋转椭球子午线弧长的解析解，为子午线弧长计算的实用化和误差控制提供了理论依据.本文的理论公式严密精确解，计算工作量小；而级数展开［10－12]、数值积分［13]、Hermite插值［14]、Horner求和
技术［15]（该算法至少可以计算到3600完全阶次的球谐级数式，可把高阶次重
力场模型向更高阶次扩展）都是近似解，计算工作量大，且本文第4节表明“并
非取级数项越多越好”.本文的求解思路有可能较好解决与椭球表面积相等的球半
径R2［2]、椭球面梯形图幅面积等［1]的理论计算问题.
在本文工作的基础上，再做一些理论工作（例如第一类完全椭圆积分、勒让德第一类椭圆积分等），这些理论公式可供实际工程直接使用［16－21].文献［16－20]研究两球体单峰接触问题，这种将复杂工程问题过于简单化导致的结果是：机械结构的理论振型与实验振型有时不一致，机械结构的固有频率与实验固有频率之间的相对误差较大，例如文献［18]预测的相对误差在－19.2%～16.8%之间.事实上，
在实际工程中，接触压力的椭球—Hertz分布将在两个物体中产生与所提出的椭圆接触面相协调的弹性位移［21]，相应地大量存在以第二类完全椭圆积分式（19）、勒让德第二类椭圆积分式（20），以及本文尚未涉及的第一类完全椭圆积分勒让
德第一类椭圆积分这4类积分为中间自变量，衍生出许多显函数，可能会使螺栓
结合部［22]的理论预测结果与实验数据更接近.
必须指出，文献［23]中的式（50）应为
参考文献：
［1]孔祥元，郭际明，刘宗泉.大地测量学基础［M].武昌：武汉大学出版社，2002：67－73.
［2]程鹏飞，文汉江，成英燕，等.2000国家大地坐标系椭球参数与GRS 80和WGS 84的比较［J].测绘学报，2009，38（3）：189－194.
［3]廖瑞祥，王刚，邹良超.基于GIS的三峡库区滑坡空间数据库系统设计［J].三峡大学学报：自然科学版，2011，33（1）：24－27.
［4]王海栋，柴洪洲，王敏.多波束测深数据的抗差Kriging拟合［J].测绘学报，2011，40（2）：238－242，248.
［5]过家春，赵秀侠，徐丽，等.基于第二类椭圆积分的子午线弧长公式变换及解
算［J].大地测量与地球动力学，2011，31（4）：94－98.
［6]同济大学数学系.高等数学（上册）［M].6版.北京：高等教育出版社，2011：116－117.
［7]Mindlin Raymond pliance 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]陈俊勇.中国现代大地基准——中国大地坐标系统2000（CGCS 2000）及其
框架［J].测绘学报，2008，37（3）：269－271.
［10]牛卓立.以空间直角坐标为参数的子午线弧长计算公式［J].测绘通报，2001（11）：14－15.
［11]边少锋，纪兵.等距离纬度等量纬度和等面积纬度展开式［J].测绘学报，2007，36（2）：218－223.
［12]李厚朴，边少锋，陈良友.等面积纬度函数和等量纬度变换的直接解算公式
［J].武汉大学学报·信息科学版，2011，36（7）：843－846.
［13]刘修善.计算子午线弧长的数值积分法［J].测绘通报，2006，（5）：4－6. ［14]李厚朴，边少锋.辅助纬度反解公式的Hermite插值法新解［J].武汉大学学报·信息科学版，2008，33（6）：623－626.
［15]刘缵武，刘世晗，黄欧.超高阶次勒让德函数递推计算中的压缩因子和Horner求和技术［J].测绘学报，2011，40（4）：454－458.
［16]田红亮，朱大林，秦红玲，等.结合部法向载荷解析解修正与定量实验验证［J].农业机械学报，2011，42（9）：213－218.
［17]田红亮，朱大林，秦红玲.固定接触界面法向静弹性刚度［J].应用力学学报，2011，28（3）：318－322.
［18]田红亮，方子帆，朱大林，等.固定接触界面切向静弹性刚度问题研究［J].应用力学学报，2011，28（5）：458－464.
［19]田红亮，朱大林，秦红玲.结合面静摩擦因数分形模型的建立与仿真［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]田红亮，赵春华，朱大林，等.金属材料结合部法切向刚度修正与实验验证［J].农业机械学报，2012，43（6）：207－214.
［23]田红亮，朱大林，秦红玲.地球旋转椭球等面积纬度函数与等量纬度的相互近
似变换［J].三峡大学学报：自然科学版，2012，34（2）：71－75.。