地面上两点间方位角和距离计算的实用公式_王瑞

　第21卷　第5期V o l . 21　N o . 5海军工程大学学报　2009年10月　　　　JO U RN A L O F NA V A L U NI VERSI T Y O F ENG IN EERIN G 　　　　　O ct . 2009　

地面上两点间方位角和距离计算的实用公式

王　瑞1, 2, 纪　兵, 边少锋11

(1. 海军工程大学电气与信息工程学院, 武汉430033; 2. 海军司令部航海保证部, 天津300042)

摘　要:在涉及距离计算的很多工程应用时, 经常需要计算地面上位置已知的两点间的方位角和距离, 解算这一问题时一般是将地球近似成旋转椭球来进行, 而实际情况是计算点并不在参考椭球上, 而是距参考椭球有大地高H , 因此使用理论距离公式计算误差太大。为解决该问题, 提出了利用相似椭球法构造新的参考椭球, 并在前人研究的基础上借助计算机代数系统M a thematica 进行了推导, 得到了简单、实用的距离和方位角计算公式。最后, 进行了数据验证, 结果表明了该方法的正确性。

关键词:大地主题问题; 相似椭球; 大地高; 大地线; M athematica

中图分类号:T N 911　　　　文献标志码:A 　　　　文章编号:1009-3486(2009) 05-0022-05

Practical formulae for calculating azimuth and distance between two

points on the ground

WANG Rui 1, 2, JI Bing , BIAN Shao -feng 11

(1. College of Electrical and Information Engineering , Naval Univ . of Engineering , Wuhan 430033, China ;

2. The Navigation Guarantee Dept . , Navy Headquarters , Tianjin 300042, China )

A bstract :As fo r calculating the azim uth and distance betw een the transmitter and the receiver , the general method is to take the earth as ro ta ry ellipsoid approxim ately , and then calculate the azim uth and distance o n its surface , but the real point po sitio n is no t on the surface . As the re is a distance of geodetic height H betw een point position and the ellipsoid , the theo retic distance fo rmula can not be used directly . In o rder to solve the problem , a new similar ellipsoid w as constructed . Then , based on the former research and w ith the help of the co mputer alg ebra sy stem M athematica , a simple and prac -tical formula for calculating the azimuth and distance w as g iven . The data com putation result prove s the cor rectness of the method .

Key words :geode tic problem ; similar ellipsoid ; geodetic heig ht ; geodesic line ; M athematica

在无线电电波传播距离计算、利用无线电电波测向定位[1]、利用时域散射场进行目标方位探测[2]、大地测量以及航海、航空航程计算等领域内, 经常会遇到求解地球上两坐标已知的点位间的大地线长度及方位角的问题, 这是经典大地主题的反解问题。当两点间距离超过一定范围(10km 以上) 时, 对地面上同一坐标系统中两点间的距离计算, 就不能按照直线来处理, 而要考虑地球的曲率, 同时为了计算结果的精确, 也不能视地球为圆球, 一般的做法是将地球近似为有一定扁率的椭圆绕短半轴旋转得到的旋转椭球, 则其上地面两点间最短距离是大地线。许多学者对椭球面上距离的计算进行了深入的探索和研究, 提出了许多行之有效的解决方法[3-6]。笔者也对该问题进行了探索[7-10], 并且借助于Mathem a -

收稿日期:2008-12-20; 修回日期:2009-04-26。

基金项目:国家自然科学基金资助项目(40774002) ; 国家杰出青年科学基金资助项目(40125013) 。

作者简介:王　瑞(1965-) , 男, 博士生, 主要研究方向为航海导航和卫星测高数据的处理与应用, E -mail :ruiw1965@

126. com 。

　第5期　王　瑞等:地面上两点间方位角和距离计算的实用公式[11, 12]　　　·23·tica 计算机代数系统, 利用He rm ite 插值法推导

出了较以往方法更简洁的非迭代的正反计算公式。

但以往在研究这一问题时无一例外都是在参考椭球

上进行的, 而实际应用时所求点并不在参考椭球上,

所求点与参考椭球面有一距离, 即大地高H 。同时,

由于地形起伏变化的原因, 所求两点间的大地高往往

不一致, 其实际情况如图1所示, 这导致用理论计算图1　大地线相对参考椭球的关系Fig . 1　Relation of g eodesic line to reference ellipsoid 公式直接计算误差太大, 不能用于工程实践。为了解决这一问题, 文中巧妙地利用相似椭球法构造了一个新的参考椭球, 在此基础上利用法截线方位角加入改正项计算大地线方位角, 再借助于Mathe -matica 计算机系统解决了椭球积分问题, 从而推导出了可以方便工程应用的实用公式。

1　相似椭球的构造与椭球参数的解算

设地面上两已知点坐标分别为P 1(B 1, L 1, H 1) 和P 2(B 2, L 2, H 2) , 坐标系统所在参考椭球的长半

H 1+H 2、与原参2

考椭球相似的椭球, 由于P 1和P 2的大地高一般是不相等的, 而新构造的椭球夹在这两点之间, 所以可轴为a 原、短半轴为b 原、第一偏心率为e 。为了解算, 笔者构造了过两点间高度为H 近似地将两已知点在相似椭球上的投影点间距离当作两点间的实际距离, 在新的参考椭球上进行相应的解算。所谓相似椭球是指中心重合、三轴重合、偏心率相等的椭球, 则由相似性(第二等式是近似相等) 有:

原原=≈, a 原b 原N

则可以得到:

Δa =原H =H N 22-e sin B =(1e sin B ) H ; 2(1) (2)

原原Δb =H =H =H -e 1-e sin B =(1-e 2-e 2sin 2B ) H 。(3) N N 22

式中:Δa , Δb 为由于大地高变化而导致的相较于原椭球元素的变化量, 则构造后的参考椭球元素为

12a =a 原+Δa , b =b 原+Δb ; N 为卯酉圈曲率半径, 计算该值时, 用两已知点间纬度的平均值B (该2

12点理论上应与取H 的点位于同一点, 实际上很难做到, 但由于地球椭球的偏心率很小, 接近2

原, 从而可以计算出相似1-e sin B

椭球的相关参数。由相似椭球的特性知, 构造后的相似椭球与原椭球的偏心率相等, 也为e , 至此完成于圆球, 故这两点不重合导致的误差可以忽略) 代入公式N =了相似椭球的构造与椭球参数的解算。

2　两点间方位角和距离的求解

由于参考椭球的特殊性, 对椭球面元素的解算比较复杂, 因此在解算大地主题问题时, 为了计算的方便, 一般是将椭球面元素按一定的原则投影到球面[2, 5]上, 在球面上解算大地问题, 再将求得的球面元素按投影关系换算为椭球面上相应的元素。贝塞尔在解这问题时提出了如下的投影条件:

1) 椭球面大地纬度投影到球面为归化纬度, 且tan u =-e tan B 。其中:B , u 分别为大地纬度和归化纬度; e 为参考椭圆第一偏心率。

2) 椭球面大地线投影到球面为大圆弧。

·24·海　军　工　程　大　学　学　报　　　　　　　　　　　第21卷　

3) 大地方位角投影后数值不变。

由以上的投影条件, 略去推导, 给出椭球面上与球面上对应元素的微分关系:

d S =a 基本方程, 其中S 、σ分别为大地线长和大圆弧。

2. 1　方位角的求取

假设地球上两已知点的大地坐标分别为P 1(B 1, L 1, H 1) 和P 2(B 2, L

2, H 2) , 则本文的目的是求取这两点间的大地线长S 以及方位角A 1和A 2。推导的思路是:

先求从P 1到P 2点的方位角, 再利用大地线长公式求出两点间

的大地线长; 由于大地线与法截线非常接近, 故可以先利用两点

大地坐标先求出法截线方位角, 再加前人已给出由球面上大圆

弧表示的法截线与大地线之间的截面差改正即可求出大地方位

角(见图2) 。通常情况下, 大地线位于相对法截线之间, 并且靠

近正法截线, 它分相对法截线的夹角约为2∶1, 即γ=2μ。具体

求取过程如下。

首先, 直接给出椭球面上两点处的卯酉圈曲率半径计算公

式(推导略去) :

N i =图2　法截线与大地线的关系Fig . 2　Relationship betw een g eode sic line and no rmal sec tion line 1-e cos u d σ。(4) 　　上式为椭球面大地线长与球面大圆弧的微分关系式, 亦称为贝塞尔微分方程, 是解大地主题问题的, 　i =1, 2-e sin B i

以及换算成直角坐标系后的P 2点的坐标:X 2=N 2co s B 2co s (L 2-L 1) ,

Y 2=N 2cos B 2sin (L 2-L 1) ,

Z 2=N 2(1-e 2) sin B 2。(5) (6)

式中:e =a -b ) /a ; a 、b 分别为参考椭球的长短半轴。由此可计算出P 1点处的法截线方位角A ′1:

A ′1=arctan (Y 2/((Z 2+N 1·e 2·sin B 1) cos B 1-X 2sin B 1) ) 。

在此, 笔者用A ′1近似代替大地方位角后可求出大地线后向延长至赤道处的近似方位角A ′0:

A ′0=arcsin (co s u 1sin A ′1) 。

　　再依据投影条件(1) , 略去推导直接给出已知点从椭球投影到球面上对应的归化纬度:

u i =arctan (-e tan B i ) , 　i =1, 2(9)

以及投影到球面后, 在球面上由直角三角形正弦定理可得的由归化纬度和赤道处的方位角表示的球面弧长:

i =arcsin 0) σ(sin u i /cos A ′, 　i =1, 2。(10)

[4]　　由于可以取近似ΔS /N 1≈Δσ, 同时Δσ=σ2-σ1, 而熊介、陈健和晁定波[6]已经给出了大地线方位(7) (8)

角A 1与法截线方位角A ′1之差, 截面差的具体表达式经变形后有:22ΔA 1=A 1-A ′1′cos 2B 1sin A 1co s A ′1-tan B 1e 4N 16N 1

2e ′co s 2B 1sin A ′1(cos A ′1-tan B 1) 。642(11)

　　则由式(11) 可以计算出法截线方位角的改正量, 有了这个改正量就可以计算赤道处更准确的法截线方位角来代替大地线方位角A 0, 即

A 0=arcsin (cos u 1sin (A ′1+ΔA 1) ) 。

则由式(12) 求出的A 0可以计算P 1点处精确的大地线方位角:(12)

　第5期　王　瑞等:地面上两点间方位角和距离计算的实用公式A 1=a rcsin (sin A 0/co s u 1) 。　　　·25·(13)

　　同理, 有A 2=a rcsin (sin A 0/cos u 2) , 则可以计算出P 2点处的大地线方位角, 其中A ′2与ΔA 2的计算与A ′1与ΔA 1的计算方法一样(应用时需判断其所在的象限角) 。

2. 2　大地线长度的计算

利用以上计算的加入改正项大地线方位角后, 可以得到更精确的球面弧长的表达式:

sin σ=sin u /co s A 0。

代入式(4) 并经整理和变形可得:

d S =a =a

式中:b =a 22(14) 　　有了球面弧长后, 利用式(4) 的微分关系就可以计算投影换算回椭球面上大地线的长度, 将式(14) 1-e (1-sin u ) d σ=a -e 2-e +e co s A 0·sin σd σ1+εsin σd σ。2221+e ′cos A 0sin σd σ=b (15) -e ; ε=e ′co s A 0, e ′为椭球第二偏心率, 且有e ′=e /(1-e ) 。

从式(15) 可以看出, 大地线长是由球面弧长表示的微分公式。因为它是一椭圆积分, 该式无法直接积分解出, 为此笔者利用计算机代数系统Mathem atica 强大的符号计算功能, 通过级数展开命令将该式

6按二项式定理展开至ε项并积分, 则有:

S =b

式中:

1

β

=0

004-800-6432-2560256-10241024-3∫21σSeries [σ21+εsin σ,{0, ε, 6}]d σ=b [ασ+βsin 2σ+γsin 4σ+δsin 6σ]|σ。1σ(16) 1642。

　　式(16) 即为两点间大地线长的实用计算公式, 式中球面大圆弧σ1, σ2可由加入改正项后的弧长公式σi =arcsin (sin u i /cos A 0) (i =1, 2) 方便地计算出。至此完成了借助于球面弧长解算椭球面上大地线长和方位角的过程。

3　数据验证

为了验证以上推导公式的正确性, 文中进行了数据试算。由于大地主题问题的特殊性, 无法直接验证, 为此仅验证本文方法与传统方法不计大地高直接计算的结果之差。参照文献[12]取已知点坐标为P 1(40. 043244°, 115. 166667°) , P 2(23. [**************]°, 122. [1**********]°) , 计算结果如表1所示。

表1　本文方法与传统方法计算结果比较

Tab . 1　Comparison of calculation by traditional m ethod and this paper

(H 1, H 2)

(-100, 300)

(0, 350)

(100, 200)

(200, 500) 直接解算/m 2000000. 002000000. 002000000. 002000000. 00本文方法/m 2000031. 392000054. 892000047. 062000109. 72差值/m 31. 3954. 8947. 06109. 72

　　数据验证结果表明, 如果不计大地高的影响, 直接用传统公式对大地线长进行计算的话, 其误差是很大的, 不能满足精密工程应用的需要, 而使用相似椭球推导后的公式则可以弥补这一缺陷, 从而提高

·26·海　军　工　程　大　学　学　报　　　　　　　　　　　第21卷　工程应用的精度水平。

4　结束语

文中对椭球面上两点间距离和方位角的计算问题进行了研究, 针对以往公式中没有考虑大地高的情形, 利用相似椭球法构造一新的参考椭球进行了解算, 并依据椭球面与球面的投影关系, 在球面上解算相关参数, 简化了计算过程, 使得该公式更适合工程实践应用。同时, 在推导大地线方位角时利用投影条件式(3) , 巧妙地借助于法截线方位角加入截面差改正项来计算大地线方位角, 得到了较精确的大地线方位角计算公式。数据验证结果显示, 文中大地线长计算公式比不考虑大地高的传统计算方法, 精度要高许多, 可以提高工程应用的计算精度。本文的推导过程是借助于计算机代数系统M athematica 完成的, 笔者从中发现这种先进的数学工具可以解决以往数学推演中人工难以实现的问题。该工具具有较好的推广应用价值。

参考文献(References ) :

[1]　刘庆彬, 刘金霞. 基于无线电测向的台站定位[J ]. 电波科学学报, 2003, 18(2) :211-215.

LI U Qing -bin , L IU Jin -xia . Station lo ca tion based o n r adio directio n -finding [J ]. Chinese Journal of Radio Scie nce , 2003, 18(2) :211-215. (in Chine se )

[2]　刘　波, 高本庆, 薛　峰. 利用时域散射场进行目标方位探测研究[J ]. 电波科学学报, 2005, 20(6) :741-744.

LI U Bo , G AO Ben -qing , XU E Feng . Study of target o rientatio n identification using scattering field in time domain [J ]. Chinese Jo urna l of Radio Science , 2005, 20(6) :741-744. (in Chinese )

[3]　胡江强, 杨盐生, 李铁山. 恒向线航向和航程的精确计算[J ]. 大连海事大学学报, 2005, 31(2) :11-14.

H U Jia ng -qiang , YA NG Yan -sheng , LI T ie -shan . A ccurate calcula tion o f the co urse and distance in r humb -line [J ]. Journal o f Dalian M aritime U niver sity , 2005, 31(2) :11-14. (in Chinese )

[4]　熊　介. 椭球大地测量学[M ]. 北京:解放军出版社, 1988.

[5]　朱华统, 黄继文. 椭球大地计算[M ]. 北京:八一出版社, 1993.

[6]　陈　健, 晁定波. 椭球大地测量学[M ]. 北京:测绘出版社, 1991.

[7]　纪兵, 边少锋. 大地主题问题的非迭代新解[J ]. 测绘学报, 2007, 36(3) :269-273.

JI Bing , BIA N Shao -feng . T he new no n -itera tive solution to the g eodetic pro blem [J ]. A cta G eode tica et Car to -g raphica Sinica , 2007, 36(3) :269-273. (in Chine se )

[8]　BIA N Shao -feng , CH EN Yong -bing . Solving and inve rse problem o f a meridian arc in ter ms of co mputer alg ebra

sy stem [J ]. Jo urnal of Surv eying Eng ineering , 2006, 132(1) :7-10.

[9]　边少锋, 柴洪洲, 金际航. 大地坐标系与大地基准[M ]. 北京:国防工业出版社, 2005.

[10]　纪　兵, 边少锋. 三轴椭球表面积的计算[J ]. 海军工程大学学报, 2008, 20(4) :21-24.

JI Bing , BIA N Shao -fe ng . Ca lculatio n of surface ar ea of t riaxial ellipsoid [J ]. Jour nal o f N av al U niver sity of Eng i -neering , 2008, 20(4) :21-24. (in Chine se )

[11]　边少锋, 许江宁. 计算机代数系统与大地测量数学分析[M ]. 北京:国防工业出版社, 2004.

[12]　洪维恩. 数学运算大师M athema tica 4[M ]. 北京:人民邮电出版社, 2002.