留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

兰勃脱等角圆锥投影反解不同算法的解析

丁士俊 李鹏鹏 邹进贵 金银龙

丁士俊, 李鹏鹏, 邹进贵, 金银龙. 兰勃脱等角圆锥投影反解不同算法的解析[J]. 武汉大学学报 ● 信息科学版, 2022, 47(9): 1452-1459. doi: 10.13203/j.whugis20200383
引用本文: 丁士俊, 李鹏鹏, 邹进贵, 金银龙. 兰勃脱等角圆锥投影反解不同算法的解析[J]. 武汉大学学报 ● 信息科学版, 2022, 47(9): 1452-1459. doi: 10.13203/j.whugis20200383
DING Shijun, LI Pengpeng, ZOU Jingui, JIN Yinlong. Analysis of Different Algorithms for Reverse Solution of Lambert Conic Conformal Projection[J]. Geomatics and Information Science of Wuhan University, 2022, 47(9): 1452-1459. doi: 10.13203/j.whugis20200383
Citation: DING Shijun, LI Pengpeng, ZOU Jingui, JIN Yinlong. Analysis of Different Algorithms for Reverse Solution of Lambert Conic Conformal Projection[J]. Geomatics and Information Science of Wuhan University, 2022, 47(9): 1452-1459. doi: 10.13203/j.whugis20200383

兰勃脱等角圆锥投影反解不同算法的解析

doi: 10.13203/j.whugis20200383
基金项目: 

国家自然科学基金 41871373

武汉大学“教育教学改革”建设引导专项 2022

详细信息
    作者简介:

    丁士俊,博士,教授,研究方向为现代测量数据处理理论、测量与地图投影理论与方法。shjding@sgg.whu.edu.cn

    通讯作者: 李鹏鹏,博士。276149876@qq.com
  • 中图分类号: P282

Analysis of Different Algorithms for Reverse Solution of Lambert Conic Conformal Projection

Funds: 

The National Natural Science Foundation of China 41871373

the Wuhan University Construction of Education and Teaching Reform 2022

More Information
    Author Bio:

    DING Shijun, PhD, professor, specializes in survey data processing, theories and methods of surveying and map projection.E-mail: shjding@sgg.whu.edu.cn

    Corresponding author: LI Pengpeng, PhD. E-mail: 276149876@qq.com
  • 摘要: 在测量与地图制图中,等量纬度求解大地纬度是一种常见的投影反解计算,就该反解问题的几种不同算法进行研究,包括迭代法、等量纬差求解大地纬度的级数展开式及等量纬度求解大地纬度的直接算法。利用Mathematica对后两种算法的计算公式进行了详细推导,给出了其高阶系数展开式,同时对现有算法中存在的问题进行了解析。兰勃脱等角投影算例表明,所推导的公式其计算精度可达(1×10-7)″~(1×10-8)″,完全满足测量与地图投影高精度的要求。
  • 表  1  兰勃脱等角割圆锥投影正解

    Table  1.   Forward Calculation for Lambert Conic Conformal Projection with Two Standard Parallels

    点号 B L X/m Y/m
    1 23°30'25.369 43″ 46°50'47.284 55″ 2 946 710.860 07 1 188 342.791 39
    2 28°00'45″ 45°30'25″ 3 444 391.831 74 1 049 914.561 05
    下载: 导出CSV

    表  2  兰勃脱等角割圆锥投影反解

    Table  2.   Reverse Calculation for Lambert Conic Conformal Projection with Two Standard Parallels

    点号 迭代法(ε=5×10-14 rad) 等量纬差法 直接法
    1 B 23°30'25.369 429 99″ 23°30'25.369 429 99″ 23°30'25.369 429 76″
    L 46°50'47.284 550 09″ 46°50'47.284 550 09″ 46°50'47.284 550 09″
    |dB| (1×10-8)″ (1×10-8)″ (2.4×10-7)″
    2 B 28°00'45.000 000 034″ 28°00'45.000 000 035″ 28°00'44.999 999 945″
    L 45°30'25.000 000 121″ 45°30'25.000 000 122″ 45°30'25.000 000 122″
    |dB| (3.4×10-8)″ (3.5×10-8)″ (4.5×10-8)″
    下载: 导出CSV

    表  3  等量纬差反解大地纬度的比较

    Table  3.   Comparison of Reverse Solutions to Geodetic Latitude byConformal Latitude Difference

    计算方法 1号点 2号点 系数项
    B |dB| B |dB|
    文献[12]算法(式(37)) 23°30'25.369 433 56″ (3.6×10-6)″ 28°00'45.368 167 85″ 0.681 68″ *t1~t5
    23°30'25.369 330 48″ (9.9×10-5)″ 28°00'44.505 910 06″ 0.494 08″ *t1~t3
    本文算法 23°30'25.369 429 99″ (1.0×10-8)″ 28°00'45.000 000 035″ (3.5×10-8)″ t1~t5
    23°30'25.369 429 32″ (6.8×10-7)″ 28°00'45.999 998 692″ (1.3×10-6)″ t1~t3
    注:带*的为文献[12]给出的系数
    下载: 导出CSV
  • [1] 杨启和. 地图投影变换原理与方法[M]. 北京: 解放军出版社, 1989

    Yang Qihe. The Principle and Method of Map Projection Transformation[M]. Beijing: The People's Liberation Army Press, 1989
    [2] 杨启和. 测量和地图学中应用的六种纬度及其变换关系式[J]. 测绘科技通讯, 1995(3): 14-19 https://www.cnki.com.cn/Article/CJFDTOTAL-CHKJ199503004.htm

    Yang Qihe. Six Latitudes and Their Transformation Relations Used in Surveying and Cartography[J]. Geomatics Technology and Equipment, 1995(3): 14-19 https://www.cnki.com.cn/Article/CJFDTOTAL-CHKJ199503004.htm
    [3] 李厚朴, 边少锋, 陈良友. 等面积纬度函数和等量纬度变换的直接解算公式[J]. 武汉大学学报·信息科学版, 2011, 36(7): 843-846 http://ch.whu.edu.cn/article/id/588

    Li Houpu, Bian Shaofeng, Chen Liangyou. The Direct Calculating Formulae for Transformations Between Authalic Latitude Function and Isometric Latitude[J]. Geomatics and Information Science of Wuhan University, 2011, 36(7): 843-846 http://ch.whu.edu.cn/article/id/588
    [4] 李厚朴, 边少锋, 刘敏. 地图投影中三种纬度间变换直接展开式[J]. 武汉大学学报·信息科学版, 2013, 38(2): 217-220 http://ch.whu.edu.cn/article/id/6108

    Li Houpu, Bian Shaofeng, Liu Min. Direct Expansions of Transformations Between Three Kinds of Latitudes Used in Map Projection[J]. Geomatics and Information Science of Wuhan University, 2013, 38(2): 217-220 http://ch.whu.edu.cn/article/id/6108
    [5] Snyder J P. Map Projections Used by the U S Geological Survey[M]. Washington, USA: United States Government Printing Office, 1982
    [6] Snyder J P. Map Projections―A Working Manual[M]. Washington, USA: United States Government Printing Office, 1987
    [7] Grafarend E W, You R J, Syffus R. Map Projections[M]. Berlin, Germany: Springer, 2014
    [8] 杨启和, 杨晓梅. 测量和地图学中应用的三种纬度函数及其反解变换的线性插值方法[J]. 测绘学报, 1997, 26(1): 92-94 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB701.016.htm

    Yang Qihe, Yang Xiaomei. Linear Interpolation Method of Three Latitude Functions and Their Inverse Solutions Used in Surveying Cartography[J]. Acta Geodaetica et Cartographic Sinica, 1997, 26(1): 92-94 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB701.016.htm
    [9] 郑彤, 边少锋. 子午线弧长反问题新解[J]. 武汉大学学报·信息科学版, 2007, 32(3): 255-258 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200703016.htm

    Zheng Tong, Bian Shaofeng. New Solutions on Inverse Problem of Meridian Arc[J]. Geomatics and Information Science of Wuhan University, 2007, 32(3): 255-258 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200703016.htm
    [10] Craig R. Auxiliary Latitude Formulas: Finding the Coefficients Numerically and Symbolically[C]//Wolfram Technology Conference, Champaign, USA, 2006
    [11] 李国藻, 杨启和, 胡定荃. 地图投影[M]. 北京: 解放军出版社, 1993

    Li Guozao, Yang Qihe, Hu Dingquan. Map Projections [M]. Beijing: The People's Liberation Army Press, 1993
    [12] 孔祥元, 郭际明, 刘宗泉. 大地测量学基础[M]. 武汉: 武汉大学出版社, 2010

    Kong Xiangyuan, Guo Jiming, Liu Zongquan. Foundation of Geodesy[M]. Wuhan: Wuhan University Press, 2010
    [13] International Association of Oil & Gas Producers. Coordinate Conversions and Transformations Including Formulas[EB/OL]. [2018-01-25]. http//www. epsg. org/guides/index. html
    [14] 边少锋, 纪兵. 等距离纬度等量纬度和等面积纬度展开式[J]. 测绘学报, 2007, 36(2): 218-223 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB200702017.htm

    Bian Shaofeng, Ji Bing. The Expansions of Rectifying Latitude, Conformal Latitude and Authalic Latitude[J]. Acta Geodaetica et Cartographica Sinica, 2007, 36(2): 218-223 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB200702017.htm
    [15] 王瑞, 李厚朴. 辅助纬度反解公式的Lagrange级数法推演[J]. 海洋测绘, 2008, 28(3): 18-23 https://www.cnki.com.cn/Article/CJFDTOTAL-HYCH200803007.htm

    Wang Rui, Li Houpu. The Derivation of the Inverse Expansions for Auxiliary Latitudes by Lagrange Series Method[J]. Hydrographic Surveying and Charting, 2008, 28(3): 18-23 https://www.cnki.com.cn/Article/CJFDTOTAL-HYCH200803007.htm
    [16] 李厚朴, 边少锋. 辅助纬度反解公式的Hermite插值法新解[J]. 武汉大学学报·信息科学版, 2008, 33(6): 623-626 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200806020.htm

    Li Houpu, Bian Shaofeng. Derivation of Inverse Expansions for Auxiliary Latitudes by Hermite Interpolation Method[J]. Geomatics and Information Science of Wuhan University, 2008, 33(6): 623-626 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200806020.htm
    [17] 过家春, 李厚朴, 庄云玲, 等. 依不同纬度变量的子午线弧长正反解公式的级数展开[J]. 测绘学报, 2016, 45(5): 560-565

    Guo Jiachun, Li Houpu, Zhuang Yunling, et al. Series Expansion for Direct and Inverse Solutions of Meridian in Terms of Different Latitude Variables[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(5): 560-565
    [18] 陈成, 边少锋, 刘强. 纬度正反解问题的傅里叶级数解算方法[J]. 海洋测绘, 2017, 37(4): 19-23

    Chen Cheng, Bian Shaofeng, Liu Qiang. A Fourier Series Method for Forward and Inverse Solutions of Latitudes[J]. Hydrographic Surveying and Charting, 2017, 37(4): 19-23
    [19] 边少锋, 李厚朴. 大地测量计算机代数分析[M]. 北京: 科学出版社, 2018

    Bian Shaofeng, Li Houpu. Computer Algebra Analysis on Geodesy[M]. Beijing: Science Press, 2018
    [20] Bian S F, Chen Y B. Solving an Inverse Problem of a Meridian Arc in Terms of Computer Algebra System[J]. Journal of Surveying Engineering, 2006, 132(1): 7-10 doi:  10.1061/(ASCE)0733-9453(2006)132:1(7)
    [21] 丁大正. Mathematica基础与应用[M]. 北京: 电子工业出版社, 2013

    Ding Dazheng. Mathematica Foundation and Application[M]. Beijing: Publishing House of Electronics Industry, 2013
  • [1] 叶彤, 李厚朴, 钟业勋, 金立新.  常用纬度与归化纬度差异极值分析表达式 . 武汉大学学报 ● 信息科学版, 2022, 47(3): 473-480. doi: 10.13203/j.whugis20200004
    [2] 焦晨晨, 李厚朴, 边少锋, 钟业勋, 潘雄.  以地心纬度为变量的常用纬度正反解符号表达式 . 武汉大学学报 ● 信息科学版, 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454
    [3] 宋伟伟, 何成鹏, 辜声峰.  不同纬度区域电离层增强PPP-RTK性能分析 . 武汉大学学报 ● 信息科学版, 2021, 46(12): 1832-1842. doi: 10.13203/j.whugis20210243
    [4] 王谦, 赵学胜, 王政, 刘青平.  一种改进的QTM地址码与经纬度坐标转换算法 . 武汉大学学报 ● 信息科学版, 2020, 45(2): 303-308, 316. doi: 10.13203/j.whugis20170390
    [5] 刘晓刚, 徐婧林, 张素琴, 管斌.  地磁日变数据确定中顾及纬度和经度方向影响的双因子定权方法 . 武汉大学学报 ● 信息科学版, 2020, 45(10): 1547-1554. doi: 10.13203/j.whugis20180454
    [6] 朱永兴, 谭述森, 明锋, 崔先强.  顾及经纬度方向异性的电离层TEC IDW插值及精度分析 . 武汉大学学报 ● 信息科学版, 2019, 44(11): 1605-1612. doi: 10.13203/j.whugis20180233
    [7] 景一帆, 杨元喜, 曾安敏, 明锋.  北斗区域卫星导航系统定位性能的纬度效应 . 武汉大学学报 ● 信息科学版, 2017, 42(9): 1243-1248. doi: 10.13203/j.whugis20150011
    [8] 曾小牛, 李夕海, 刘志刚, 杨晓君, 刘代志.  低纬度磁异常化极及分量换算的正则化方法 . 武汉大学学报 ● 信息科学版, 2016, 41(3): 388-394. doi: 10.13203/j.whugis20140342
    [9] 张晓平, 边少锋, 李忠美.  极区高斯投影与日晷投影的比较 . 武汉大学学报 ● 信息科学版, 2015, 40(5): 667-672. doi: 10.13203/j.whugis20140128
    [10] 姚朝龙, 罗志才, 刘立龙, 周波阳.  顾及地形起伏的中国低纬度地区湿延迟与可降水量转换关系研究 . 武汉大学学报 ● 信息科学版, 2015, 40(7): 907-912. doi: 10.13203/j.whugis20130409
    [11] 李厚朴, 边少锋, 刘敏.  地图投影中三种纬度间变换直接展开式 . 武汉大学学报 ● 信息科学版, 2013, 38(2): 217-220.
    [12] 李胜全, 李厚朴, 边少锋.  拉格朗日投影与常用等角投影间解析变换的复变函数表示 . 武汉大学学报 ● 信息科学版, 2012, 37(11): 1382-1385.
    [13] 李厚朴, 边少锋, 陈良友.  等面积纬度函数和等量纬度变换的直接解算公式 . 武汉大学学报 ● 信息科学版, 2011, 36(7): 843-846.
    [14] 李厚朴, 边少锋.  辅助纬度反解公式的Hermite插值法新解 . 武汉大学学报 ● 信息科学版, 2008, 33(6): 623-626.
    [15] 童晓冲, 张永生, 贲进.  经纬度坐标与QTM编码的三向互化算法及其精度评价标准 . 武汉大学学报 ● 信息科学版, 2006, 31(1): 27-30.
    [16] 程芦颖.  虚拟地形投影与地球同步卫星定位解全集 . 武汉大学学报 ● 信息科学版, 2005, 30(7): 629-631.
    [17] 丁佳波.  采用双重投影法的椭球面日晷投影 . 武汉大学学报 ● 信息科学版, 2000, 25(2): 183-185.
    [18] 刘彩璋, 黄杰安, 戴四军.  新天文常数、系统对天文经纬度和方位角的影响 . 武汉大学学报 ● 信息科学版, 1988, 13(3): 47-55.
    [19] 刘彩璋.  固体潮对天文经纬度的影响 . 武汉大学学报 ● 信息科学版, 1986, 11(1): 104-107.
    [20] 许厚泽.  应用高斯-克吕格投影解大地主题的初步意见 . 武汉大学学报 ● 信息科学版, 1957, 1(0): 69-85.
  • 加载中
计量
  • 文章访问数:  78
  • HTML全文浏览量:  25
  • PDF下载量:  21
  • 被引次数: 0
出版历程
  • 收稿日期:  2022-09-16
  • 刊出日期:  2022-09-05

兰勃脱等角圆锥投影反解不同算法的解析

doi: 10.13203/j.whugis20200383
    基金项目:

    国家自然科学基金 41871373

    武汉大学“教育教学改革”建设引导专项 2022

    作者简介:

    丁士俊,博士,教授,研究方向为现代测量数据处理理论、测量与地图投影理论与方法。shjding@sgg.whu.edu.cn

    通讯作者: 李鹏鹏,博士。276149876@qq.com
  • 中图分类号: P282

摘要: 在测量与地图制图中,等量纬度求解大地纬度是一种常见的投影反解计算,就该反解问题的几种不同算法进行研究,包括迭代法、等量纬差求解大地纬度的级数展开式及等量纬度求解大地纬度的直接算法。利用Mathematica对后两种算法的计算公式进行了详细推导,给出了其高阶系数展开式,同时对现有算法中存在的问题进行了解析。兰勃脱等角投影算例表明,所推导的公式其计算精度可达(1×10-7)″~(1×10-8)″,完全满足测量与地图投影高精度的要求。

English Abstract

丁士俊, 李鹏鹏, 邹进贵, 金银龙. 兰勃脱等角圆锥投影反解不同算法的解析[J]. 武汉大学学报 ● 信息科学版, 2022, 47(9): 1452-1459. doi: 10.13203/j.whugis20200383
引用本文: 丁士俊, 李鹏鹏, 邹进贵, 金银龙. 兰勃脱等角圆锥投影反解不同算法的解析[J]. 武汉大学学报 ● 信息科学版, 2022, 47(9): 1452-1459. doi: 10.13203/j.whugis20200383
DING Shijun, LI Pengpeng, ZOU Jingui, JIN Yinlong. Analysis of Different Algorithms for Reverse Solution of Lambert Conic Conformal Projection[J]. Geomatics and Information Science of Wuhan University, 2022, 47(9): 1452-1459. doi: 10.13203/j.whugis20200383
Citation: DING Shijun, LI Pengpeng, ZOU Jingui, JIN Yinlong. Analysis of Different Algorithms for Reverse Solution of Lambert Conic Conformal Projection[J]. Geomatics and Information Science of Wuhan University, 2022, 47(9): 1452-1459. doi: 10.13203/j.whugis20200383
  • 在测量与地图投影中,常常涉及到不同辅助纬度正反解计算及其相互变换[1-10]。等量纬度是一种重要的辅助纬度,在兰勃脱投影与正轴墨卡托投影中,已知等量纬度求解大地纬度属于投影反解计算。由地图投影可知[11-12],等量纬度是大地纬度的函数,其反解函数关系式是一个较为复杂的隐函数,其解算方法一般采用基于正解公式的迭代算法[513],缺点是计算效率较低,不便于理论与数值分析。杨启和[1]采用泰勒级数展开式将等量纬差展开为大地纬差的幂级数,导出了等量纬差求解大地纬差的级数展开式,该方法的关键是如何计算等量纬度与大地纬度的近似值,但对此没有加以阐述;孔祥元等[12]在兰勃脱投影反解中运用了该计算方法,给出了在投影坐标原点纬度处展开的计算方法,但该算法有时会导致计算不收敛,存在较大的计算误差。另外一类算法即直接法,早期杨启和[1-2]对此作了研究,通过复杂的级数展开及求导计算,运用拉格朗日级数方法,给出了等量纬度求解大地纬度的直接算法的计算过程,但由于手工推导过程复杂,未能给出其反解直接算法的通用计算公式,仅列出基于正解公式和特定参考椭球参数下反解系数的数值表达式。十多年来,学者们借助于Mathematica计算代数系统,基于泰勒级数、拉格朗日级数、Hermite插值及傅里叶级数等方法[14-18],对辅助纬度正反解算法进行了大量的研究,其中,王瑞等[15]利用拉格朗日级数原理将导出的反解公式与杨启和[2]给出的计算公式从形式上作了比较,两者存在一定的偏差。王瑞等[15]给出了基于拉格朗日级数法等量纬度正反解的系数关系式,但缺少详细的推导与论证过程。

    基于此,本文借助于Mathematica计算机代数系统[19-21]对上述后两种算法从理论上进行了重新推导,导出了更为精确的计算公式,给出了其高阶系数展开式。通过算例将本文算法与上述后两种算法进行计算分析与比较,其计算精度可达(1×10-7)″~(1×10-8)″,完全满足测量与地图投影高精度的要求。

    • 由地图投影可知,大地纬度B是等量纬度q的函数,因此可以表示为:

      B=f(q)

      B=B0+ΔB=f(q0+Δq)B0=f(q0),采用泰勒级数将式(1)展开为:

      B=B0+ΔB=f(q0)+dBdq0Δq+12d2Bdq02Δq2+ 13!d3Bdq03Δq3++1n!dnBdq0nΔqn (2)

      式中,n为序号。

      tn=1n!dnBdq0n,则式(2)可改写为:

      B=f(q0)+t1Δq+t2Δq2+t3Δq3++tnΔqn

      由等量纬度与大地纬度的微分关系可知:

      dq=MNcosBdB

      式中,M=a(1-e2)/(1-e2sin2B)3为子午线曲率半径,a为椭球长半径,e为椭球第一偏心率;N=a/1-e2sin2B为卯酉线曲率半径。由式(4)可得:

      dBdq=(1+η2)cosB

      式中,η2=(e')2cos2Be'为椭球第二偏心率,且满足(e')2=e2/(1-e2)

      由式(5)可知,一阶导数dB/dq是大地纬度B的函数,可利用复合函数依次求各阶导数,即:

      dnBdqn=d(dn-1Bdqn-1)dBdBdq

      由于高阶导数计算比较复杂,借助Mathematica计算各阶导数,令t=tanB,舍弃六阶导数以上高阶导数项,忽略tη及其乘积6次以上高次项,则有:

      B=B0+t1Δq+t2Δq2+t3Δq3++t6Δq6

      式中,各系数分别为:

      t1=(1+η02)cosB0t2=12t0(-1-4η02-3η04) cos2B0t3=16(-1+t02-5η02+13t02η02-7η04) cos3B0t4=124t0(5-t02+56η02-40t02η02+154η04) cos4B0t5=1120(5-18t02+t04+61η02-418t02η02+ 210η04) cos5B0t6=1720t0(-61+58t02-t04-1 324η02+ 2 632t02η02-7 153η04) cos6B0

      文献[1]推算的结果为:

      t1=(1+η02) cosB0t2=12t0(-1-4η02-3η04) cos2B0t3=16(-1+t02-5η02+13t02η02-7η04) cos3B0t4=124t0(5-t02+56η02-40t02η02) cos4B0t5=1120(5-18t02+t04) cos5B0t6=1720t0(-61+58t02-t04) cos6B0

      将式(8)与式(9)进行比较可知,上述各系数的主项基本一致,式(8)保留了tη及其乘积的4次项系数,其系数项的表达式更精确。

    • 由地图投影理论可知,等量纬度与大地纬度的闭合式为:

      q=lntan(π4+B2)(1-esinB1+esinB)e2

      如果令e=0,则地球椭球体变为球体,椭球面大地纬度B变为球面纬度φ,则有:

      q=lntan(π4+φ2)

      由式(10)和式(11)可得:

      tan(π4+φ2)=tan(π4+B2)(1-esinB1+esinB)e2

      由式(12)可得:

      φ=2arctantan(π4+B2)(1-esinB1+esinB)e2-π2

      式中,φ称为球面等角纬度。

      由于椭球第一偏心率e是一个较小的数值,等角纬度φ与大地纬度B相差也很小,因此可将式(13)展开为e的幂级数形式:

      φ(B,e)=φ(B,0)+φe|e=0 e+122φe2|e=0 e2+ 13!3φe3|e=0 e3++1n!nφen|e=0 en

      借助Mathematica计算各导数,令e=0,得到各阶导数为:

      φe|e=0=3φe3|e=0=5φe5|e=0= 7φe7|e=0=9φe9|e=0=02φe2|e=0=-sin2B4φe4|e=0=-52(2sin2B-sin4B)6φe6|e=0=-32(45 sin2B-42 sin4B+13 sin6B) 8φe8|e=0=14(-7 868 sin2B+9 758 sin4B- 5 532 sin6B+1 237 sin8B)10φe10|e=0=-452(4 704 sin2B-6 696 sin4B+ 5 079 sin6B-2 096 sin8B+367 sin10B)

      将式(15)导数代入式(14),合并同类项得到表达式为:

      φ=B+P2sin2B+P4sin4B++P2nsin(2nB) (16)

      式中,各系数分别为:

      P2=-e22-5e424-3e632-281e85 760-7e10240-P4=5e448+7e680+697e811 520+93e102 240+P6=-13e6480-461e813 440-1 693e1053 760-P8=1 237e8161 280+131e1010 080+P10=-367e10161 280-

      就CGCS2000椭球而言,则有:

      φ=B-692.339 74 sin2B+ 0.968 33 sin4B-(1.69×10-3)sin6B+(3.21×10-6)sin8B-(6.31×10-9)sin10B (18)

      式中,系数P10系数项的最大值为(6.31×10-9),因此计算时取至P8项,其计算精度可达到(1×10-7)~(1×10-8)″。

      综上,由等量纬度的闭合式(10)可直接求解等量纬度,也可利用式(11)与式(16)求解等量纬度。首先,由大地纬度计算等角纬度;然后,利用等角纬度计算大地纬度,从而实现等量纬度的正解。

    • 等量纬度求解大地纬度可由闭合式(10)得到如下迭代式:

      B=2 arctaneq(1+esinB1-esinB)e2-π2

      式中,e为自然对数的底,e = 2.718 281 828 459…;e‍为椭球第一偏心率。

      在大地测量与地图投影中,辅助纬度(如归化纬度、地心纬度、等角纬度、等距纬度与等积纬度)通常可表达为如同式(16)三角函数表达式[1-25]。由式(16)可得:

      B=φ+f(B)

      式中,

      f(B)=-(P2sin2B+P4sin4B++ P2nsin(2nB) (21)

      由文献[2]可知,为了求得式(16)的反解公式,由式(20)利用拉格朗日级数公式可得:

      B=φ+f(φ)+12!d(f2(φ))dφ+ 12!d2(f3(φ))dφ2+1n!dn-1(fn(φ))dφn-1 (22)

      进一步可得到等角纬度解大地纬度的计算式:

      B=φ+a2sin2φ+a4sin4φ++a2nsin(2nφ)

      为了得到式(23),运用Mathematica计算式(22)中各阶导数,将其导数代入式(22),通过三角函数约化将三角函数的幂级数转换为三角函数倍角函数,合并同类项,整理后得到式(23)。

      本文研究表明,如果计算的系数精确至e10项,则取n5;若系数项取至e8,则取n4即可。本文取n=5,经推导得到反解式(23)与正解式(16)的系数关系为:

      a2=-P2-P2P4-P4P6+12P23+P2P42- 12P22P6+13P23P4-112P25a4=-P4+P22-2P2P6+4P22P4-43P24a6=-P6+3P2P4-3P2P8-32P23+92P2P42+ 9P22P6-272P23P4+278P25a8=-P8+2P42+4P2P6-8P22P4+83P24a10=-P10+5P2P8+5P4P6-252P2P42- 252P22P6+1256P23P4-12524P25

      式(24)可将系数精确展开至e10项。如果仅将系数展开到e8项,式(24)可简化为:

      a2=-P2-P2P4+12P23a4=-P4+P22-2P2P6+4P22P4-43P24α6=-P6+3P2P4-32P23a8=-P8+2P42+4P2P6-8P22P4+83P24

      文献[1-2]利用拉格朗日级数公式,当n=4时,取:

      f(B)=- (P2sin2B+P4sin4B+ P6sin6B+P8sin8B) (26)

      通过复杂的人工演算得到反解式:

      B=φ+a2sin2φ+a4sin4φ+a6sin6φ+a8sin8φ

      其中,

      a2=-P2-P2P4-P4P6+12P23+P2P42- 12P22P6-18.3P23P4a4=-P4+P22-2P2P6+4P22P4-43P24a6=-P6+3P2P4-3P2P8-32P23+ 92P2P42+9P22P6-12.5P23P4a8=-P8+2P42+4P2P6-8P22P4+83P24

      将式(17)各系数代入式(24)可得:

      a2=e22+5e424+e612+13e8360+3e10160+a4=7e448+29e6240+811e811 520+81e102 240+a6=7e6120+81e81 120+3 029e1053 760+a8=4 279e8161 280+883e1020 160+a10=2 087e10161 280+

      比较式(24)与式(28)可知,文献[2]推演的正反解系数关系与本文推导的结果低次项系数相同,高次项系数存在一定差异,在文献[15]中也给出相同的结论。本文通过分析比较,文献[215]导出的系数关系只能精确到e8系数项,两者高阶系数与式(24)中的高阶系数(e10系数项)都存在一定的偏差,其原因在于利用拉格朗日级公式求解时级数项数取值不同,本文推导的系数关系式(24)与文献[19]给出的公式一致,该正反解系数关系式可精确计算到e10项,因此该系数关系更为严密,在大地测量文献资料中,常称之为三角级数回求公式。本文推导的等角纬度反解系数关系式(29)与文献[18-19]导出的表达式相同。由此可以验证辅助纬度正反解系数关系式(24)的正确性。

      采用椭球第三偏心率可简化式(29)的系数表达式:

      a2=2e3-2e323-2e33+116e3445+26e3545a4=7e323-8e335-227e3445+2 704e35315a6=56e3315-136e3435-1 262e35105a8=4 279e34630-332e3535a10=4 174e35315

      其中,

      e3=1-1-e21+1-e2

      e3称为椭球第三偏心率,e34达到10-12量级,e35可达10-14量级。式(29)与式(30)比较而言,显然式‍(30)计算收敛速度更快。就CGCS2000椭球而言有:

      B=φ+692.339 09 sin2φ+1.355 55 sin4φ+ (3.64×10-3) sin6φ+(1.11×10-5) sin8φ+ (3.59×10-8) sin10φ+ (32)

      式中,系数a10项的最大值为(3.59×10-8),因此计算取至a8项即可。

    • 兰勃脱等角割圆锥投影,其投影参数为:椭球长半径a、偏率f、第一标准纬度B1、第二标准纬度B2、原点经度L0、原点纬度B0、原点东偏移量FE和北偏移量FN

      大地坐标(B,L)计算投影直角坐标(X,Y)正解算法[12-13]

      X=FN+ρ0-ρ cosγY=FE+ρ sinγ
      β=1q2-q1ln(N1cosB1N2cosB2)K=N1cosB1βe-βq1K=N2cosB2βe-βq2ρ=Ke-βq,B0Bρ0ρλ=L-L0γ=βλ

      式中,q=ln[tan(π4+B2)(1-esinB1+esinB)e2]N=a1-e2sin2Be‍为椭球第一偏心率,e=2f-f2;由B1B2可以计算出N1N2;由纬度B0B1B2B分别计算等量纬度q0q1q2q

      由投影直角坐标(X,Y)反解大地坐标(B,L)算法,计算大地纬度采用如下3种算法:

      1)迭代法:

      B(i+1)=2atan[eq(1+esinB(i)1-esinB(i))e/2]-π2

      式中,e为自然对数的底;e‍为椭球第一偏心率。迭代时,取初始值B(0)=2atan(eq)-π/2,直到B(i+1)-B(i)ε迭代结束,迭代阈值取ε=5×10-14

      2)等量纬差求解大地纬度:

      B=B0+t1Δq+t2Δq2+t3Δq3+t4Δq4+ t5Δq5+t6Δq6+ (36)

      式中,各系数参见式(8)。文献[12]将式(7)大地纬度在投影坐标原点纬度B0处展开,由兰勃脱投影则有等量纬差:

      Δq=-1βlnρρ0

      本文算法取值为:

      B0=2atan(eq)-π/2q0=ln[tan(π4+B02)(1-esinB01+esinB0)e2] Δq=q-q0

      3)等量纬度反解大地纬度直接算法:

      φ=2 atan(eq)-π/2B=φ+a2sin2φ+a4sin4φ++a2nsin(2nφ)

      计算大地经度:

      L=L0+γ/β
      x=ρ0-(X-FN)y=Y-FEρ=±x2+y2,β q=-ln(ρ/K)/βγ=atan(y/x)

      式中,βKρ0的计算与正解相同。

    • 在测量与制图中,兰勃脱投影在世界范围内被许多国家尤其是“一带一路”沿线及周边国家广泛采用,本文以沙特Ainel Abd 1970大地坐标系为例,该大地坐标系的投影参数如下:参考椭球为International 1924;a=6 378 388.0 m,1/f= 297.0;原点纬度为B0=24°N,经度为L0=45°E;第一标准纬度为B1=21°N;第二标准纬度为B2=27°N;坐标格网东偏移量为FE=1 000 000.0 m;北偏移量为FN=3 000 000.0 m。

      为了便于计算与分析,假定1号点纬度与原点纬度相差约30',2号点纬度与原点纬度约相差4°。兰勃脱等角割圆锥投影正解计算结果如表 1所示。

      表 1  兰勃脱等角割圆锥投影正解

      Table 1.  Forward Calculation for Lambert Conic Conformal Projection with Two Standard Parallels

      点号 B L X/m Y/m
      1 23°30'25.369 43″ 46°50'47.284 55″ 2 946 710.860 07 1 188 342.791 39
      2 28°00'45″ 45°30'25″ 3 444 391.831 74 1 049 914.561 05

      按照本文推导的计算式进行投影反解,采用迭代法、等量纬差法以及直接法3种方法进行大地纬度的反解,计算结果见表 2。从表 2可以看出,当迭代法的阈值取ε=(1×10-8)″=5×10-14rad时(迭代5~6次即可满足),3种算法的计算精度基本相当,纬度的计算精度可达(1×10-7)″~(1×10-8)″。因此文献[12]等量纬差法与本文算法比较结果见表 3,由于坐标原点B0不是纬度B的近似值,故q0也不是q的近似值,文献[12]算法计算收敛速度慢,在有限项级数展开的情况下,计算点纬度与原点纬度相差越大,其误差越大(如2号点),若要提高计算精度,必须增加级数展开的项数,这样使得计算复杂化,难以把握。

      表 2  兰勃脱等角割圆锥投影反解

      Table 2.  Reverse Calculation for Lambert Conic Conformal Projection with Two Standard Parallels

      点号 迭代法(ε=5×10-14 rad) 等量纬差法 直接法
      1 B 23°30'25.369 429 99″ 23°30'25.369 429 99″ 23°30'25.369 429 76″
      L 46°50'47.284 550 09″ 46°50'47.284 550 09″ 46°50'47.284 550 09″
      |dB| (1×10-8)″ (1×10-8)″ (2.4×10-7)″
      2 B 28°00'45.000 000 034″ 28°00'45.000 000 035″ 28°00'44.999 999 945″
      L 45°30'25.000 000 121″ 45°30'25.000 000 122″ 45°30'25.000 000 122″
      |dB| (3.4×10-8)″ (3.5×10-8)″ (4.5×10-8)″

      表 3  等量纬差反解大地纬度的比较

      Table 3.  Comparison of Reverse Solutions to Geodetic Latitude byConformal Latitude Difference

      计算方法 1号点 2号点 系数项
      B |dB| B |dB|
      文献[12]算法(式(37)) 23°30'25.369 433 56″ (3.6×10-6)″ 28°00'45.368 167 85″ 0.681 68″ *t1~t5
      23°30'25.369 330 48″ (9.9×10-5)″ 28°00'44.505 910 06″ 0.494 08″ *t1~t3
      本文算法 23°30'25.369 429 99″ (1.0×10-8)″ 28°00'45.000 000 035″ (3.5×10-8)″ t1~t5
      23°30'25.369 429 32″ (6.8×10-7)″ 28°00'45.999 998 692″ (1.3×10-6)″ t1~t3
      注:带*的为文献[12]给出的系数

      本文算法实质上是将大地纬度在球面等角纬度处展开,等角纬度与大地纬度相差较小,可作为大地纬度的近似值,故收敛速度快,运用t1~t5系数项计算,其精度可达(1×10-8)″,若仅采用t1~t3系数项,其精度可达(1×10-6)″。

    • 本文借助于Mathematica计算机代数系统强大的计算功能,导出了等量纬差求解大地纬度级数展开式以及等量纬度反解大地纬度直接算法计算式。得出以下结论:(1)推导等量纬差反解大地纬度的级数展开式,并与文献[1-2]推导的公式比较,其系数关系的主项基本相同,本文推导的结果保留tη及其乘积4次项,计算精度更高。(2)指出了文献[12]等量纬差反解大地纬度计算方法的缺陷,给出了正确算法模型,本文算法的计算精度高,收敛速度快。(3)‍利用拉格朗日级数公式,从理论上详细地推导了辅助纬度正反解三角函数展开式的系数关系式,该系数关系式可将三角函数的系数项精确展开至e10项,同时指出了本文推导公式与文献‍[1-215]给出公式之间的差异性,并给出了合理的解释。运用本文导出的系数关系式,得到了等量纬度解大地纬度的三角函数关系式(即直接法公式)。(4)通过兰勃脱割圆锥投影算例计算,迭代法、本文给出的等量纬差法与直接法三者计算精度基本相当,迭代法计算效率相对较低,等量纬差法算法模型相对复杂,而直接法模型关系式相对简单,尤其是采用椭球第三偏心率表达的系数关系式更为简洁,计算效率更高。

参考文献 (21)

目录

    /

    返回文章
    返回