三种坐标间转换的雅可比矩阵数值导数计算方法_王解先

第24卷第4期2004年11月

大地测量与地球动力学

JOU RNAL OF GEODESY AND GEODYNAM ICS

Vol. 24, No. 4

Nov. , 2004

文章编号:1671-5942(2004) 04-0019-05

三种坐标间转换的雅可比矩阵数值导数计算方法

王解先 徐志京

(同济大学测量系, 上海 200092)

X

大地坐标以及平面坐标与大地高之间相互转换的雅可比矩阵的数摘 要 提出了一种用于空间大地直角坐标、

值导数的计算方法。该方法可以用较容易且精确的替代直接方法来实现不同类型间的坐标转换。在空间大地直角坐标与平面坐标与大地高之间的转换实验测试过程中, 只要选取合适的增量, 数值导数的精度都会满足要求, 并且数值导数的精度的一致性也得到了验证。另外, 该方法比采用直接方法计算速度显著加快, 数值导数的精度、一致性和减少计算耗时在实验测试中效果非常好。

关键词 坐标转换 雅可比矩阵 数值导数 精度 计算方法

中图分类号:P226+. 3 文献标识码:A

METHOD OF CALCULATING NUMERICAL DERIVATIVE OF

JACOBIN MATRIX FOR TRANSFORM AMONG SPACE CARTESIAN, GEODETIC AND PLAN COORDINATES

Wang Jiex ian and Xu Zhijing

(Department of Surveying o f Tongj i University , Shanghai 200092)

Abstract

In this paper, w e present a method of calculating the numerical derirative of Jacobin matrix for

transform among space Cartesian, geodetic and plan coordinates. It is accurate and easy to use in parameter subs-i tution and transformation from one type of coordinates to another. T he computation by this method will be faster dramatically than in ordinary w ay to transform space Cartesian to plan coordinates. The accuracy of numerical partial is tested and is high enough.

Key words:coordinates transformation, Jacobin matrix , numerical derivative, accuracy, calculating method

或者用带有大地高的平面坐标系中的(x y h ) 表示, 由于理论研究和工程计算的需要, 常常涉及不同坐标方式的转换问题。不同坐标间的转换都有具体的转换公式, 例如, 点坐标(X Y Z ) 可以采用下面的公式转换成(B L h) T 的形式:

T

T

1 引言

在大地测量学中, 一个点的坐标可以用不同的坐标方式表示, 例如, 在空间大地直角坐标系中用(X Y Z) T 表示, 在大地坐标系中用(B L h ) T 表示,

X 收稿日期:2004-06-10

基金项目:国家自然科学基金(40234039)

男, 年生, 教授,

20

大地测量与地球动力学24卷

X Y =Z

(N +h) cos (B ) cos (L ) (N +h) cos (B ) sin (L )

(1)

不能轻易写出, 因为转换的公式既长又复杂, 从(X Y Z) T 到(x y h) T 转换雅可比矩阵也就更加难以推导。在本文中, 我们将介绍一种非常容易的数值导数方法来计算这些转换偏导数, 该方法的精度和计算时间也将进行检验。

(N (1-e 2) +h) sin (B )

式中N =(2)

1-e sin (B)

a 是大地坐标系对应椭球之长半轴, e 2是椭球的第一偏心率的平方。如果我们选定WGS84椭球, 那么a =6378137m , e 2=0. [1**********]29。然而点坐标从(B L h) 转换到(x y h ) 的公式或者从(X Y Z ) T 转换到(X g Y g h) T 的公式则更加冗长和复杂。

各种形式的坐标之间的雅可比矩阵也是非常有用的。例如, 如果一种形式的坐标的统计信息已知, 我们想要计算另外一种坐标形式的统计信息时, 或者如果一种坐标的增量已知, 我们想获得另外一种坐标形式的增量等情况就会用到。根据公式(1) 我们可以把从(X Y Z) T 到(B L h ) T 转换雅可比矩阵写成:

=J =这里=5B ==cos (B) cos (L ) -(N +h) sin (B) cos (L ) 5B

-(N +h) cos (B ) sin (L ) cos (B ) cos (L )

(3)

T

T

2 数值导数的概念

若y =f (x 1, x 2, , , x i , , ) , y 到x i 的导数的定义为:=i lim D y 0

f (x 1, x 2, , , x i +D , , ) -f (x 1, x 2, , , x i , , )

(4)

当D 足够小时, 这一导数可以近似写成数值导数的形式:U i

f (x 1, x 2, , , x i +D , , ) -f (x 1, x 2, , , x i , , )

这里, (x y h ) T 到(X Y Z) T 转换导数关系就满足这一条件。

D x D y =D h

D Y

D Z

D X

(5)

式中, (D x D y D h ) T 是(x y h ) T 的小增量, (D X D Y D Z ) T 是(X Y Z) T 的小增量, 我们定义:

=

a 11a 21a 31

a 12a 22a 32

T

a 13a 23a T

(6)

=cos (B ) sin (L ) -(N +h) sin (B ) sin (L ) =(N +h) cos (B ) cos (L ) =cos (B ) sin (L ) ===(1-e 2) sin (B) +(N (1-e 2) +h ) cos (B ) 0sin (B)

2

对测量学专家来说, 从(X Y Z) 到(x y h) 的转换过程非常熟悉。如果在X 分量上添加一个小的改变D , 被转换的平面直角坐标和大地高将出现一个小的增量($x $y $h ) T , 参照公式(5) , 我们得到:

a 11=, a 21=, a 31=(7)

如果在Y 和Z 分量上分别添加一个小的改变D , 那么我们同样可以得到a 12, a 22, a 32, a 13, a 23, a 33。这就是数值方式计算雅可比矩阵的主要思路。采用这一方法我们可以很容易地计算出用于3种形式的坐标之间转换的雅可比矩阵。下面我们讨论该方法的精度。

=(1-e 2sin 2(B ) ) 3/2

从(X Y Z) T 到(B L h) T 转换是一个反复的过程,

但雅可比矩阵的逆矩阵=J -1可以轻易获

得(参见公式(3) ) 。

L ) T ) T 3 数值导数的精度

根据数值导数的定义, D 应该趋向于零。但是, D

第4期王解先等:三种坐标间转换的雅可比矩阵数值导数计算方法

21

所采用的计算机的性能。公式(3) 可以用来计算, 同样可以计算如表1中的数值导数。这

里我们选取了地球上的3个点来研究数值导数的精度。

表1 用于测试的3个点坐标Tab. 1 Coordinates of three points for test

Point A B C

B (b ) 14585

L (b ) 291291291

H (m) 101010

5/5

表2和表3中所取的结果是在B 和L 选择D =10-9rad 、h 选择D =0. 01m 的基础上取得的。如果我们选取更小的D , 所得到的结果会更好, 但是, 我们认为这样的结果已经满足要求了。

下面我们检验经常用于GPS 导航定位中的从(X Y Z) T 到(x y h) T 之间坐标转换中的数值导数的精度。

表3 B 点和C 点的数值导数误差(%)

T ab. 3 Erro rs o f the numerical derivative at point B and C (%)

B 0. 06910672

0. 069118970. 0700461B 0. 00126980. 00126960. 0012314

Po int B L 0. 069572610. 0695535 -Po int C L

0. 00126850. 0012706 -h 0. 000004-0. 0000080. 000003h

0. 0000040. 0000010. 000001

对于点A , 根据公式(3) 和WGS84椭球体, 我们得到:=

-39626. [1**********]4. 114220. 358103231. [1**********]8. 85070-0. [1**********]7. 39622

0. 017(8)

采用数值导数的方法, 对于B 和L 如果我们选择D =10-9(弧度, 相当于6mm ) ; 对h , 选择D =0. 01m , 可得:Nu

X Y Z 5/5X Y Z

4 数值导数的精度和一致性检验

为了检验数值导数的精度和一致性, 我们把从(X Y Z) T 到(x y h) T 转换的直接方式和数值方式进行比较。U TM 选用的中央子午线为290b 。如果我们选择一个点, 它的经度为291b , 高度为10m , 它的纬度在0b 到90b 之间变动。

如果X 、Y 和Z 3个分量的增量均为100m 。对平面坐标和大地高的计算采用直接方式且认为是标准的, 图1到图3表明了直接方式与数值导数方式的误差。

从图1到图3我们可以看出, 直接方式和数值方式求导的误差值不超过1m m 。我们曾经测试了增量小于1m 时的情况, 两者的误差值几乎为零。同样的, 图4到图6表明了3个分量的增量均为1000m 的情况, 从中我们可以看出, 该情况下直接方式和数值方式求导的误差值也不超过0. 1m 。

下面我们来检验混合增量时直接方式和数值方式的求导误差情况。

选定表1中的A 、B 和C 3个点, 假设所有的增量之和为d , 由计算机产生3个随机值r 1, r 2, r 3, 那么X 、Y 和Z 3个分量的增量为:d

r 1

r 1+r 2+r 3

=

5953612. 73736

0. 3580. 017(9)

-39624. [1**********]3. 80701

103224. [1**********]8. 05378-0. 93344

两者之间的百分比差值如表2所示。

表2 A 点数值导数误差(%)

Tab. 2 Errors of the numerical derivative at point A(%)5/5X

Y Z

B 0. 0066670. 0067750. 005739

L 0. 0057340. 005723

-h 0. 000000-0. 0000020. 000003

偏导数意味着如果(B L h) T 有一个小的增量, (X Y Z) T 的增量将被扩大多少倍。公式(8) 中的最后一个值是指如果h 有一个小的变化D , Z 的变化将被乘上0. 01745。表2中表明数值导数的精度优于0. 01%。当然, 该精度将取决于所选择的D 以及所选点所处的位置。同样的, 表3给出了B 点和C 点的数值导数误差, 从中可以看出, 对于每一个偏导数分量来说, 其误差都不大于1%。

d

r 2

r 1+r 2+r 3

d

r 3

r 1+r 2+r 图7到图9分别表明了直接方式和数值方式之间的误差值。

22

大地测量与地球动力学24卷

第4期王解先等:三种坐标间转换的雅可比矩阵数值导数计算方法

23

5 数值求导的计算时间

如果我们采用直接的方式把(X Y Z) T 转换成(x y h) , 我们需要首先将空间直角坐标转换成大地坐标系(B L h) T , 然后转换(B L ) T 到(x y ) T , 对于所有的点转换程序都是相同的。如果我们采用数值求导的方法转换, 首先计算转换用的雅可比矩阵, 如果所有点处于同一个域(根据要求精度确定) , 那么所有点都可以采用同一个雅可比矩阵进行转换。图10比较了采用两种不同方式转换一定数量点时所需的计算时间。从图中可以看出, 如果转换点的数量越来越多, 数值求导方式将会节省大量的计算时间。当然计算时间将取决于所用计算机性能的好坏。但我们认为相对时间是相同的。图10表明, 如果我们转换10个点, 直接方式将需要357s, 数值方式仅需37s 。

6

T

阵数值导数计算方法, 该方法同传统的直接计算方法相比, 在选择合适增量的前提下, 计算的精度高, 计算方法的一致性好。特别是在海量数据进行坐标转换时, 计算时间远远少于直接方法所用时间。这一方法也可以用于雅可比矩阵难于推导的其他领域。

References

1 朱华统. 常用大地坐标及其变换[M ].北京:解放军出版社, 1990.

1 Zhu Huatong. Common geodeti c coordinates and transformation [M ].Beijing:People . s Liberation Army Press, 1990. (i n Ch-i nese)

2 沈剑华. 数值计算基础[M ].上海:同济大学出版社, 1999. 2 Shen Jianhua. Basic numerical computation [M ]. Shanghai:

Tongji U niversity Press, 1999. (i n Chinese) 3 周长发. 科学与工程数值算法(Visual C +

清华大学出版社, 2002.

3 Zhou Changfa. Sci ence and engineering numerical computation

method(V ersion of Visual C ++) [M]. Bei jing:Tsighua Univer -sity Press, 2002. (in Chinese)

+

版) [M ].北京:

6 结论

由于理论研究和工程计算的需要, 常常涉及不同

坐标方式的转换, 本文提出了坐标转换中的雅可比矩

新书简介

由赖锡安、黄立人、徐菊生等编著的5中国大陆现今地壳运动6一书, 已由地震出版社出版发行。

全书由四部分组成:第一部分为地质背景介绍(第2~3章) ; 第二部分是本书的重点, 论述了应用GPS 技术观测和研究中国大陆地壳水平运动的方法和结果(第4~9章) ; 第三部分总结了利用30多年来精密水准测量资料和地面跨断层观测资料研究地壳垂直运动和断裂带运动的成果(10~12章); 第四部分简要介绍基于力学模式的大地测量数据反演(第13章) 。在最后的附录中给出了我国大陆近1000个GPS 测站上观测到的现今地壳水平运动速度场和大量的参考文献。

书中的研究成果是作者近10多年来承担和参与多项国家课题的研究总结, 同时汇集了国内外相关领域的新成果和新进展。

参与此书编写的有丁国瑜院士以及朱文耀、任金卫、周硕愚、晁定波、赵少荣、牛之俊、游新兆、王敏等。此书可供从事地壳运动和地球动力学研究的科研人员和工程技术人员以及高等院校地球科学专业的师生阅读参考。如需订阅此书, 可与本刊编辑部联系。