留言板

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

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

以地心纬度为变量的常用纬度正反解符号表达式

焦晨晨 李厚朴 边少锋 钟业勋 潘雄

焦晨晨, 李厚朴, 边少锋, 钟业勋, 潘雄. 以地心纬度为变量的常用纬度正反解符号表达式[J]. 武汉大学学报 ( 信息科学版), 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454
引用本文: 焦晨晨, 李厚朴, 边少锋, 钟业勋, 潘雄. 以地心纬度为变量的常用纬度正反解符号表达式[J]. 武汉大学学报 ( 信息科学版), 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454
JIAO Chenchen, LI Houpu, BIAN Shaofeng, ZHONG Yexun, PAN Xiong. Symbolic Expressions for Direct and Inverse Solutions to Common Latitudes with the Variable of Geocentric Latitude[J]. Geomatics and Information Science of Wuhan University, 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454
Citation: JIAO Chenchen, LI Houpu, BIAN Shaofeng, ZHONG Yexun, PAN Xiong. Symbolic Expressions for Direct and Inverse Solutions to Common Latitudes with the Variable of Geocentric Latitude[J]. Geomatics and Information Science of Wuhan University, 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454

以地心纬度为变量的常用纬度正反解符号表达式

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

国家自然科学基金 41871376

国家自然科学基金 41671459

国家自然科学基金 41771487

湖北省杰出青年科学基金 2019CFA086

详细信息
    作者简介:

    焦晨晨,硕士生,主要研究方向为大地测量和地图投影。jiaocchen@163.com

    通讯作者: 李厚朴,博士,副教授。lihoupu1985@126.com
  • 中图分类号: P282

Symbolic Expressions for Direct and Inverse Solutions to Common Latitudes with the Variable of Geocentric Latitude

Funds: 

The National Natural Science Foundation of China 41871376

The National Natural Science Foundation of China 41671459

The National Natural Science Foundation of China 41771487

the Natural Science Foundation for Distinguished Young Scholars of Hubei Province 2019CFA086

More Information
    Author Bio:

    JIAO Chenchen, postgraduate, specializes in geodesy and map projection. E-mail: jiaocchen@163.com

    Corresponding author: LI Houpu, PhD, associate professor. E-mail: lihoupu1985@126.com
图(2)
计量
  • 文章访问数:  217
  • HTML全文浏览量:  61
  • PDF下载量:  23
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-02-23
  • 刊出日期:  2022-05-05

以地心纬度为变量的常用纬度正反解符号表达式

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

    国家自然科学基金 41871376

    国家自然科学基金 41671459

    国家自然科学基金 41771487

    湖北省杰出青年科学基金 2019CFA086

    作者简介:

    焦晨晨,硕士生,主要研究方向为大地测量和地图投影。jiaocchen@163.com

    通讯作者: 李厚朴,博士,副教授。lihoupu1985@126.com
  • 中图分类号: P282

摘要: 借助具有强大符号运算功能的计算机代数系统Mathematica,推导了地图投影学中等距离纬度、等角纬度、等面积纬度与地心纬度之间的正反解直接展开式,并将式中的系数统一表示成关于椭圆偏心率e和椭球第三扁率n的幂级数形式。与以往反解方法不同的是,采用符号迭代法进行以地心纬度为变量的等距离纬度、等角纬度、等面积纬度的反解,并使用最大差异值作为衡量精度的标准。算例分析表明,以地心纬度为变量的常用纬度展开式在结构和形式上与以大地纬度为变量的辅助纬度保持一致,基于第三扁率n的幂级数表达式具有更紧凑的形式和更好的收敛性,其直接展开式的精度分别优于(1×10-8)″和(1×10-10)″,可以满足大地测量和地图投影精密计算的需要。

English Abstract

焦晨晨, 李厚朴, 边少锋, 钟业勋, 潘雄. 以地心纬度为变量的常用纬度正反解符号表达式[J]. 武汉大学学报 ( 信息科学版), 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454
引用本文: 焦晨晨, 李厚朴, 边少锋, 钟业勋, 潘雄. 以地心纬度为变量的常用纬度正反解符号表达式[J]. 武汉大学学报 ( 信息科学版), 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454
JIAO Chenchen, LI Houpu, BIAN Shaofeng, ZHONG Yexun, PAN Xiong. Symbolic Expressions for Direct and Inverse Solutions to Common Latitudes with the Variable of Geocentric Latitude[J]. Geomatics and Information Science of Wuhan University, 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454
Citation: JIAO Chenchen, LI Houpu, BIAN Shaofeng, ZHONG Yexun, PAN Xiong. Symbolic Expressions for Direct and Inverse Solutions to Common Latitudes with the Variable of Geocentric Latitude[J]. Geomatics and Information Science of Wuhan University, 2022, 47(5): 707-714. doi: 10.13203/j.whugis20190454
  • 在大地测量及地图投影的相关理论中,等距离纬度ψ、等角纬度φ和等面积纬度ϑ一般以大地纬度B为自变量进行表达和运算[1-4]。然而,在地图投影、空间大地测量和地球物理等领域的相关理论中,地心纬度也是常用辅助变量之一。例如,在椭球面日晷投影[5]中,从椭球面上投影到球面上是基于地心纬度进行分析;在空间大地测量的理论问题尤其是几何问题中,以地心纬度为自变量可以使问题得到简化;在卫星轨道、大地测高等问题中,地心纬度也有着重要的作用。通常,在辅助纬度或辅助纬度函数相互变换时,用到的方法主要为反解变换法和正解变换法[6]

    国内外地图投影研究者在正解变换法上进行了深入的研究并取得丰硕的成果[7-17],不仅解决了投影变换在不同椭球下的问题,而且推导了大地纬度表示的常用纬度之间以及常用纬度函数之间的展开式[11-1216],给出了常用纬度之间差异极值的符号表达式[13],并且得到了以归化纬度、地心纬度为变量的子午线弧长的正反解展开式[17]。但在地图投影变换中涉及到的以地心纬度为变量的辅助纬度变换方面的问题,相关研究十分匮乏。事实上,在地图投影中常常会直接用到地心纬度进行计算。

    鉴于此,本文以具有强大的符号运算功能的计算机代数系统Mathematica作为公式推导的重要辅助手段,推导出以地心纬度为变量的常用纬度的正解展开式,并采用符号迭代法推导其反解展开式,以基于偏心率e和第三扁率n的幂级数形式表示,给出了理论与精度分析。

    • 以椭球长半轴ax轴,短半轴by轴,建立平面直角坐标系。则椭圆标准方程为:

      x2a2+y2b2=1

      P为椭球面上任意一点,连接OP,则OP=ρP点的向径,POX=ϕP点的地心纬度。过P点作y轴的垂线PQ,则PQ=r为纬线圈半径。将极坐标x=ρcosϕy=ρsinϕ代入式(1)可得:

      ρ=aba2sin2ϕ+b2cos2ϕ

      设椭球第一偏心率为e,并记Mϕ=ρ2(ϕ)+(ρ')2(ϕ),结合极坐标中弧长微分公式得dX=ρ2(ϕ)+(ρ')2(ϕ)dϕ=Mϕdϕ,则以地心纬度为变量的Mϕ和纬线圈半径r为:

      Mϕ=a1-e21-(2-e2)e2cos2ϕ(1-e2cos2ϕ)3r=acosϕ1-e21-e2cos2ϕ

      根据地图投影理论,椭球面上由赤道至地心纬度ϕ处的子午线弧长X为:

      X=0ϕMϕdϕ=a1-e20ϕ1-(2-e2)e2cos2ϕ(1-e2cos2ϕ)3dϕ

      设一幅角为ψ、半径为R的圆弧与子午线弧长在数值上等同,则有等距离纬度ψ为:

      ψ=XR=Xac0c0=1-14e2-364e4-5256e6- 17516 384e8-44165 536e10

      借助Mathematica,略去推导过程,得到以地心纬度为变量的等距离纬度ψ正解表达式,并在e=0处将其展开为e的幂级数形式,取到e10项,整理得到:

      ψ(ϕ)=ϕ+a2sin(2ϕ)+a4sin(4ϕ)+ a6sin(6ϕ)+a8sin(8ϕ)+a10sin(10ϕ)

      其中,

      a2= 18e2+ 116e4+ 531 024e6+ 952 048e8+ 1 35932 768e10a4=- 1256e4- 1256e6+ 58 192e8+ 214 096e10a6=- 51 024e6- 152 048e8- 1 811262 144e10a8=- 261131 072e8- 26165 536e10,a10=- 9211 310 720e10

      已知椭球第三扁率n通常作为一种小展开参数应用于大地测量计算中,n的表达式及其与偏心率e的关系如下:

      n=a-ba+b=1-1-e21+1-e2

      将系数转化为第三扁率n的幂级数形式,取至n5,地心纬度计算等距离纬度的直接表达式为:

      ψ(ϕ)=ϕ+a2'sin(2ϕ)+a4'sin(4ϕ)+ a6'sin(6ϕ)+a8'sin(8ϕ)+a10'sin(10ϕ )

      其中,

      a2'=12n+1316n3-1532n5a4'=-116n2+3332n4,a6'=-516n3+349256n5a8'=-261512n4,a10'=-9211 280n5

      显然,基于第三扁率n展开的系数比基于第一偏心率e的更为简洁,且前者只需展开至n5即可满足地图投影的要求,收敛性更好。

    • 根据地图投影理论,结合式(3),可得等量纬度q关于地心纬度ϕ的函数关系式为:

      q=0ϕMϕrdϕ=0ϕ1-(2-e2)e2cos2ϕ(1-e2cos2ϕ)cosϕdϕ

      若将地球视为圆,则e=0,有:

      φ=2arctan(exp(q))-π2

      借助Mathematica,略去推导过程,得到以地心纬度为变量的等角纬度φ正解表达式,并在e=0处将其展开为e的幂级数形式,取到e10项,整理得到:

      φ(ϕ)=ϕ+b2sin(2ϕ)+b4sin(4ϕ)+ b6sin(6ϕ)+b8sin(8ϕ)+b10sin(10ϕ )

      其中,

      b2=124e4+596e6+591 152e8+53911 520e10b4=-148e4-160e6-192 304e8-75 760e10b6=-1160e6-252 688e8-717 680e10b8=-5532 256e8-4111 520e10b10=-1123 040e10

      结合式(7),将系数转化为第三扁率n的幂级数形式,并取至n5,由地心纬度计算等角纬度的直接表达式为:

      φ(ϕ)=ϕ+b2'sin(2ϕ)+b4'sin(4ϕ)+ b6'sin(6ϕ)+b8'sin(8ϕ)+b10'sin(10ϕ )

      其中,

      b2'=23n2+23n3-29n4-1445n5b4'=-13n2+415n3+4345n4-445n5b6'=-25n3+2105n4+124105n5b8'=-55126n4-16105n5,b10'=-2245n5
    • 根据地图投影理论,在椭球面上,由赤道到地心纬度ϕ及单位经差所围成的曲边梯形的面积即为等面积纬度函数[12],一般用F表示:

      F=0ϕMϕr dϕ= a2(1-e2)0ϕcosϕ1-(2-e2)e2cos2ϕ(1-e2cos2ϕ)2dϕ

      由赤道至纬度ϑ及单位经差所围曲边梯形面积与F相等,由球面积分公式得:

      sinϑ=FR2=Fa2k0

      式中,k0=1-13e2-115e4-135e6-163e8-199e10ϑ为等面积纬度;R为半径。

      借助Mathematica,省略推导步骤,得到以地心纬度为变量的等面积纬度正解展开式,在e=0处进行展开,取至e10项,整理为:

      ϑ(ϕ)=ϕ+c2sin(2ϕ)+c4sin(4ϕ)+ c6sin(6ϕ)+c8sin(8ϕ)+c10sin(10ϕ)

      其中,

      c2=16e2+790e4+2815 040e6+27 869604 800e8+ 593 20714 968 800e10c4=1180e4+1252e6+10 6691 814 400e8+ 507 84159 875 200e10c6=-13145 360e6-8 6691 814 400e8- 537 259119 750 400e10c8=-5 9333 628 800e8-81 22923 950 080e10c10=-80 011119 750 400e10

      结合式(7),将系数转化为第三扁率n的幂级数形式,并取至n5,由地心纬度计算等角纬度的直接表达式为:

      ϑ(ϕ)=ϕ+c2'sin(2ϕ)+c4'sin(4ϕ)+ c6'sin(6ϕ)+c8'sin(8ϕ)+c10'sin(10ϕ )

      其中,

      c2'=23n-445n2+62105n3+7784 725n4- 193 082467 775n5c4'=445n2-32315n3+12 33814 175n4+92 696467 775n5c6'=-5242 835n3-1 61814 175n4+612 536467 775n5c8'=-5 93314 175n4-8 32466 825n5c10'=-320 044467 775n5
    • 在地图投影的理论推导中,经常会遇到过程略有些繁琐且不便于理解和计算的函数的反解问题,如超越函数、反三角函数等。国内外许多学者采用如埃尔米特插值法、拉格朗日级数法等解决方法,并取得相当好的效果。然而迭代法因为过程较为简单,容易理解,且可以控制精度,被广泛应用于解算复杂问题的数值解,而符号迭代法不同之处在于解算的结果为符号表达式,更加便于进行理论研究[18]。本文借助Mathematica完成采用符号迭代法进行反解问题的分析。

    • 对式(6)变形,可得迭代关系式:

      ϕi+1=ψ-a2sin(2ϕi)-a4sin(4ϕi)- a6sin(6ϕi)-a8sin(8ϕi)-a10sin(10ϕi)

      由于地心纬度ϕ和等距离纬度ψ相差很小,所以取初值ϕ0=ψ代入式(17),进行第一次迭代,得:

      ϕ1=ψ-a2sin(2ψ)-a4sin(4ψ)- a6sin(6ψ)-a8sin(8ψ)-a10sin(10ψ)

      ϕ1代入式(17)再次迭代,并将迭代结果进行偏心率e的幂级数展开,展开至e10,可得第二次迭代结果ϕ2,再次将ϕ2代入式(17)进行迭代,依此类推,直至反解展开式中的各项系数完全相同,终止迭代。以等距离纬度ψ计算地心纬度ϕ的直接展开式为例,第1次迭代结果:

      ϕ1=ψ-18sin(2ψ)e2+1256(-16sin(2ψ)+sin(4ψ))e4+(-53sin(2ψ)+4sin(4ψ)+5sin(6ψ))e61 024+[-80(76sin(2ψ)+sin(4ψ)-12sin(6ψ))+261sin(8ψ)]e8131 072+ (-54 360sin(2ψ)-6 720sin(4ψ)+9 055sin(6ψ))e101 310 720+(5 220sin(8ψ)+921sin(10ψ))e101 310 720

      略去第2~4次迭代结果,第5次、第6次迭代结果为:

      ϕ5=ϕ6=ψ-18sin(2ψ)e2+1256(-16sin(2ψ)+5sin(4ψ))e4+ sin(2ψ)[-51+20csc(2ψ)sin(4ψ)]e61 024+sin(2ψ)[12csc(2ψ)(-sin(2ψ)+sin(6ψ))]e61 024+ (-17 376sin(2ψ)+6 640sin(4ψ)+288sin(6ψ)+283sin(8ψ))e8393 216+ (-307 825sin(2ψ)+112 000sin(4ψ))e107 864 320+[20(93sin(6ψ)+566sin(8ψ))+1 301sin(10ψ)]e107 864 320

      可以发现,第6次的迭代结果与第5次的迭代结果相同。因此,以地心纬度ϕ为变量的等距离纬度ψ的反解展开式为:

      ϕ(ψ)=ψ+α2sin(2ψ)+α4sin(4ψ)+ α6sin(6ψ)+α8sin(8ψ)+α10sin(10ψ)

      其中,

      α2=-18e2-116e4-1032 048e6-1814 096e8- 61 5651 572 864e10α4=5256e4+5256e6+41524 576e8+17512 288e10α6=12 048e6+34 096e8+31131 072e10α8=283393 216e8+283196 608e10α10=1 3017 864 320e10

      将系数转化为第三扁率n的幂级数形式,并取至n5,由地心纬度计算等角纬度的直接表达式为:

      ϕ(ψ)=ψ+α2'sin(2ψ)+α4'sin(4ψ)+ α6'sin(6ψ)+α8'sin(8ψ)+α10'sin(10ψ )

      其中,

      α2'=-12n-2332n3+4991 536n5α4'=516n2-596n4α6'=132n3-77128n5α8'=2831 536n4, α10'=1 3017 680n5
    • 采用反解思路,迭代3次后即可得到以地心纬度ϕ为变量的等角纬度φ的反解展开式为:

      ϕ(φ)=φ+β2sin(2φ)+β4sin(4φ)+ β6sin(6φ)+β8sin(8φ)+β10sin(10φ)

      其中,

      β2=-124e4-596e6-29576e8-13288e10β4=148e4+160e6+232 304e8+71 152e10β6=1160e6+3448e8+1256e10β8=8332 256e8+1256e10, β10=1311 520e10

      将系数转化为第三扁率n的幂级数形式,并取至n5,由地心纬度计算等角纬度的直接表达式为:

      ϕ(φ)=φ+β2'sin(2φ)+β4'sin(4φ)+ β6'sin(6φ)+β8'sin(8φ)+β10'sin(10φ)

      其中,

      β2'=-23n2-23n3+49n4+29n5β4'=13n2-415n3-2345n4+6845n5β6'=25n3-2435n4-4635n5β8'=83126n4-8063n5β10'=5245n5
    • 采用反解思路,迭代6次后即可得到以地心纬度ϕ为变量的等面积纬度ϑ的反解展开式为:

      ϕ(ϑ)=ϑ+γ2sin(2ϑ)+γ4sin(4ϑ)+ γ6sin(6ϑ)+γ8sin(8ϑ)+γ10sin(10ϑ)

      其中,

      γ2=-16e2-790e4-1372 520e6-79 7111 814 400e8- 4 470 689119 750 400e10γ4=145e4+833 780e6+35 0271 814 400e8+ 333 09719 958 400e10γ6=-2922 680e6-3 0191 814 400e8-483 577239 500 800e10γ8=7191 209 600e8+144 037119 750 400e10γ10=7 177239 500 800e10

      将系数转化为第三扁率n的幂级数形式,取至n5,由地心纬度计算等角纬度的直接表达式为:

      ϕ(ϑ)=ϑ+γ2'sin(2ϑ)+γ4'sin(4ϑ)+γ6'sin(6ϑ)+γ8'sin(8ϑ)+γ10'sin(10ϑ)

      其中,

      γ2'=-23n+445n2-158315n3-2 10214 175n4+ 109 042467 775n5γ4'=1645n2-16945n3+93414 175n4-7 256155 925n5γ6'=-2322 835n3+92214 175n4-25 28666 825n5γ8'=7194 725n4+26818 711n5γ10'=14 354467 775n5
    • 本文借助Mathematica,选用CGCS2000坐标系[19]参考椭球常数对以地心纬度为变量的等距离纬度ψ、等角纬度φ和等面积纬度ϑ进行正反解精度分析,用以验证推导的公式,其中,a=6 378 137 mf=1/298.257 222 101

      地心纬度与大地纬度的表达式为:

      ϕ=arctan((1-e2)tanB)

      基本思路如下:首先,取定大地纬度B,将B代入式(28)计算得到地心纬度ϕ0的理论值,将B分别代入文献[20]中大地纬度表示的等距离纬度、等角纬度和等面积纬度正解表达式中,得到其理论值分别为ψ0φ0ϑ0;然后,将ϕ0代入式‍(6)、式‍(11)和式(15),即可得到基于偏心率e的以地心纬度为变量的等距离纬度ψ1、等角纬度φ1和等面积纬度ϑ1,将常用纬度理论值与计算值相减,便可得到不同正解表达式的计算误差,记为ΔψΔφΔϑ;最后,将得到的ψ0φ0ϑ0分别代入反解式(22)、式(24)和式(26),便可得到地心纬度计算值ϕ(ψ)ϕ(φ)ϕ(ϑ),将地心纬度理论值与计算值相减,便可得到不同反解表达式的计算误差,记为Δϕ(ψ)Δϕ(φ)Δϕ(ϑ)。从理论上讲,ϕ0ψ0φ0ϑ0ϕ(ψ)ϕ(φ)ϕ(ϑ)ψ1φ1ϑ1应该对应相等,由于推导的公式存在一定的误差,所以本文选用最大差异值作为精度评定标准。同理,可得到基于第三扁率n的计算误差。

      以大地纬度B为横坐标,得到以地心纬度为变量的辅助纬度正反解公式的绝对值计算误差ΔΔϕ曲线,如图 1图 2所示。

      图  1  基于e的常用纬度正、反解绝对值计算误差曲线

      Figure 1.  Absolute Error Curves of Common Latitude Direct and Inverse Solutions Based on e

      图  2  基于n的常用纬度正、反解绝对值计算误差曲线

      Figure 2.  Absolute Error Curves of Common Latitude Direct and Inverse Solutions Based on n

      图 1图 2可知:(1)利用本文公式计算的结果进行辅助纬度正反解误差很小,基于e的常用纬度正反解最大计算误差分别约为1.64×10-92.03×10-9,均小于1×10-8,基于n的常用纬度正反解绝对值最大计算误差分别约为5.26×10-114.4×10-11,均小于1×10-10,完全满足大地测量和地图投影精密计算的需要;(2)基于n的幂级数展开式的计算误差与基于e的幂级数展开式相比,误差更小,精度更高;(3)随着纬度升高,基于e的常用纬度正反解表达式绝对值计算误差总体呈先上升后下降的趋势,在70°左右达到最大误差。

    • 本文建立了适合计算机代数系统下进行分析的常用纬度变换的数学模型,借助具有强大符号运算功能的计算机代数系统Mathematica,推导了以地心纬度为变量的常用纬度正反解直接展开式,并将其分别进行基于椭球偏心率e和第三扁率n的幂级数展开,丰富了地图投影理论。

      1)本文推导的公式与以大地纬度为变量的辅助纬度正反解公式在结构和形式上保持一致,并给出了基于e和基于n的幂级数符号表达式,各纬度间正反解公式系数简洁,结构一致,便于计算机程序实现。

      2)相较于基于e的幂级数展开式的计算误差,基于n的幂级数展开至n5比基于e的幂级数展开至e10的计算精度更高。由此可说明,基于第三扁率n的幂级数表达式具有更好的收敛性和更高的精度;导出的公式计算误差分别小于1×10-81×10-10,完全满足大地测量和地图投影精密计算的需要。

参考文献 (20)

目录

    /

    返回文章
    返回