维直角坐标系直接转换算法的研究与实现.doc

测绘与空间地理信息GEOMATICS SPATIAL INFORMATION TECHNOLOGY第 35 卷 第 5 期2012 年 5 月Vol 35，No 5May ，20123维直角坐标系直接转换算法的研究与实现 达1 ，韩晓冰2 ，房龙1刘(1 黑龙江第二测绘工程院，黑龙江 哈尔滨 150081; 2 黑龙江测绘计量仪器检定站，黑龙江 哈尔滨 150081)摘 要: 测绘工作中经常碰到 3 维直角坐标系的 3 维转换问题，本文通过简单的几何平移旋转缩放算法，对 3 维十六参数的直角坐标系变换，进行了理论描述和几何证明以及数据验证。关键词: 3 维直角坐标系; 坐标转换; 变换矩阵中图分类号: P226 + 3文献标识码: B文章编号: 1672 5867( 2012) 05 0175 04The Research and Implementation of 3 Dimensional CartesianCoordinate System Direct Conversion AlgorithmLIU Da1 ，HAN Xiao bing2 ，FANG Long1( 1 The 2nd Heilongjiang Surveying and Mapping Engineering Institute，Harbin 150081，China;2 The Heilongjiang Surveying Instruments Calibration Station，Harbin 150081，China)Abstract: 3 dimensional Cartesian coordinate system conversion is often used in surveying and mapping This paper provides theo-retical description，geometric proof and data verification to 3 dimensional 16 parameter Cartesian coordinate system conversion through simple geometric translation，rotating and scaling algorithmKey words: 3 dimensional Cartesian coordinate system; coordinate conversion0引言在测绘生产中，经 常 会遇到坐标系间的坐标转换问 题，同一坐 标 系 统下进行投影换带计算，有 相 应 计 算 公 式，其计算也是 严 密 的，对 于 不 同坐标系间的投影换带， 因为其椭球参数不同，而涉及到 3 维直角坐标系 ( 空间直 角坐标系) 的 3 维转换。在 3 维直角坐标系转换中，多采 用 7 参 数 Bursa Wolf 模 型，Molodensky 模 型 和 武 测 模 型1，随 着 移 动 测 图 系 统 ( Mobile Mapping System，简 称 MMS) 技术的成熟 和 应 用，对 于 运 动 载 体 ( 飞 机、轮 船、汽 车等) 姿态的测量也越来越多。这就需要准确计算转换 参数模型，适应任意旋转角的坐标转换2。图 1 任意点在不同 3 维直角坐标系例图Fig 1 The position of arbitrary point in different3 dimensional Cartesian coordinate system1直接转换算法的理论基础假设任意 3 维 直 角 坐 标 系 下，不在同一直线上三点 A，B，C 在两个空间坐标系: 坐标系 1 ( o1x1y1z1) 和坐标系2( o2x2y2z2 ) 下的坐标分别为: A ( ax1，ay1，az1 ) ，B ( bx1， by1，bz1) ，C( cx1，cy1，cz1) 和 A( ax2，ay2，az2) ，B( bx2，by2， bz2) ，C( cx2，cy2，cz2) ( 如图 1 所示) 。1) 以坐标系 2 为基准对坐标系 1 进行缩放，设其缩放 系数为 S( S 求法见后) ;2) 将坐标系 1 与坐标系 2 变换为以 A 点为坐标原点，以 AB 为 Y 轴，并以 ABC 平面内过 A 点且垂直于 AB 的直线为 X 轴，以过 A 点且垂直于 ABC 平面的直线为 Z 轴的 新坐标系( 后简称中间坐标系) ( 如图 2 所示) ;3) 进行变换后此时坐标系 1 与坐标系 2 两个坐标系收稿日期: 2011 11 15作者简介: 刘 达( 1979 ) ，男，河南长垣人，工程师，武汉大学测绘工程专业硕士研究生，主要从事测绘生产软件开发与测绘生产辅助研究工作。BC2 = SQRT( ( BX2 CX2) 2 + ( BY2 CY2) 2 + ( BZ2 CZ2) 2) ;CA2 = SQRT( ( CX2 AX2) 2 + ( CY2 AY2) 2 + ( CZ2 AZ2) 2) 。并求出对应缩放比例，取边长缩放平均值即:S = ( AB1 / AB2 + BC1 / BC2 + CA1 / CA2) /3 0如果进行坐标 系 1 到 坐 标 系 2 的 缩 放 变 换，则 需 对ABC 在做标系 1 下的坐标进行如下变换:X1 = X1 / S; Y1 = Y1 / S; Z1 = Z1 / S如果 进 行 坐 标 系 2 到 坐 标 系 1 的 缩 放 变 换 则 需 对ABC 在坐标系 2 下的坐标进行如下变换:X2 = X1S; Y2 = Y2S; Z2 = Z2S2 2 XYZ 平移将坐标系原点 O 移至点 A，此时坐标系 1 内任意点坐 标发生的相应变化为:X = X AX1; Y = Y AY1; Z = Z AZ1坐标系 2 内任意点发生的相应变化为:X = X AX2; Y = Y AY2; Z = Z AZ2移动后如图 4 所示，即坐标系 1 下的平移参数:MX1 = AX1; MY1 = AY1; MZ1 = AZ1坐标系 2 下的平移参数:MX2 = AX2; MY2 = AY2; MZ2 = AZ2。平移后:AX1 = AX1 + MX1; AY1 = AY1 + MY1; AZ1 = AZ1 + MZ1; BX1 = BX1 + MX1; BY1 = BY1 + MY1; BZ1 = BZ1 + MZ1; CX1 = CX1 + MX1; CY1 = CY1 + MY1; CZ1 = CZ1 + MZ1; AX2 = AX2 + MX2; AY2 = AY2 + MY2; AZ2 = AZ2 + MZ2; BX2 = BX2 + MX2; BY2 = BY2 + MY2; BZ2 = BZ2 + MZ2; CX2 = CX2 + MX2; CY2 = CY2 + MY2; CZ2 = CZ2+ MZ2。图 2 等价坐标系例图Fig 2 Equivalent coordinate system完全重合，即为等价坐标系( 如图 2 所示) 。若再对坐标系 1 进行坐标系 2 的逆变换，或对坐标系2 进行坐标系 1 的逆变换，即可实现坐标系 1 与 坐 标 系 2 的相互变换( 以坐标系 2 到坐标系 1 为例，如图 3 所示) 。 整个变换以三角形 ABC 为静态参考。图 3 坐标系 1 转换为坐标系 2 例图Fig 3 The conversion of coordinate system1 to coordinate system 22直接转换的实现与相关参数的求解根据变换理论基础，该变换过程所需参数 为: 坐 标 系缩放系数 S; 坐标系 1 下三个平移参数 MX1，MY1，MZ1; 坐标系 2 下三个平移参数 MX2，MY2，MZ2; 坐标系 1 三个旋 转参数 RX1，RY1，RZ1; 坐标系 2 三个旋转参数 RX2，RY2，RZ2; 下面分别给出坐标系的变换实现与参数求法。2 1 坐标系尺度的统一求出 AB，BC，CA 在坐标系 1 和坐标系 2 的长度，求相 应边在不同坐标系 下 的 长 度 AB1，BC1，CA1，AB2，BC2，CA2。其中，AB1，BC1，CA1 为 AB，BC，CA 在坐标系 1 下的 长度，其求法为:AB1 = SQRT( ( AX1 BX1) 2 + ( AY1 BY1) 2 + ( AZ1 BZ1) 2) ;BC1 = SQRT( ( BX1 CX1) 2 + ( BY1 CY1) 2 + ( BZ1 CZ1) 2) ;CA1 = SQRT( ( CX1 AX1) 2 + ( CY1 AY1) 2 + ( CZ1 AZ1) 2) 。其中，AB2，BC2，CA2 为 AB，BC，CA 在 坐 标 系 2 下 的 长度，其求法为:AB2 = SQRT( ( AX2 BX2) 2 + ( AY2 BY2) 2 + ( AZ2图 4 坐标系平移例图Fig 4 The translation of coordinate system2 3XYZ 旋转1) 绕 Z 轴旋转绕 Z 轴旋转。其中，设绕 Z 轴旋转角度 RZ 等于在平面 XOY 上 AB 与 OY 的夹角，则坐标系 1 内任意点绕 Z 轴旋转角度 RZ1 = QFWJ( BX1，BY1，AX1，AY1) pi0 5( 其 中函数 QFWJ 在后面的算法描述中有具体算法) ; 坐标系2 内任意点旋转角度为 RZ2 = QFWJ( BX2，BY2，AX2，AY2) pi0 5( 旋转示意图如图 5) 。坐标系 1 下任意点旋转 公式为: BZ2) 2) ;第 5 期刘 达等: 3 维直角坐标系直接转换算法的研究与实现177CY1 = cos( rx1) CY1 + sin( rx1) CZ1;CZ1 = sin( rx1) CY1 + cos( rx1) CZ1;CX1 = CX1。坐标系 2 下任意点旋转公式为:y = cos( rx2) y + sin( rx2) zz = sin( rx2) y + cos( rx2) z旋转后:AY2 = cos( rx2) AY2 + sin( rx2) AZ2;AZ2 = sin( rx2) AY2 + cos( rx2) AZ2;AX2 = AX2;BY2 = cos( rx2) BY2 + sin( rx2) BZ2;BZ2 = sin( rx2) BY2 + cos( rx2) BZ2;BX2 = BX2;CY2 = cos( rx2) CY2 + sin( rx2) CZ2;CZ2 = sin( rx2) CY2 + cos( rx2) CZ2;CX2 = CX2。3) 绕 Y 轴旋转( 7)( 8)图 5 坐标系旋转例图Fig 5 The rotation of coordinate systemx = cos( rz1) x + sin( rz1) yy = sin( rz1) x + cos( rz1) y( 1)( 2)则旋转后:AX1 = cos( rz1) AX1 + sin( rz1) AY1; AY1 = sin( rz1) AX1 + cos( rz1) AY1; AZ1 = AZ1;BX1 = cos( rz1) BX1 + sin( rz1) BY1;BY1 = sin( rz1) BX1 + cos( rz1) BY1;BZ1 = BZ1;CX1 = cos( rz1) CX1 + sin( rz1) CY1; CY1 = sin( rz1) CX1 + cos( rz1) CY1; CZ1 = CZ1。坐标系 2 下任意点旋转公式为:x = cos( rz2) x + sin( rz2) yy = sin( rz2) x + cos( rz2) y旋转后:AX2 = cos( rz1) AX1 + sin( rz1) AY1; AY2 = sin( rz1) AX1 + cos( rz1) AY1; AZ2 = AZ1;BX2 = cos( rz1) BX2 + sin( rz1) BY2; BY2 = sin( rz1) BX2 + cos( rz1) BY2; BZ2 = BZ2;CX2 = cos( rz1) CX2 + sin( rz1) CY2;CY2 = sin( rz1) CX2 + cos( rz1) CY2;CZ2 = CZ2。2) 绕 X 轴旋转绕 Y 轴旋转角度 RY 等于在平面 XOZ 上 AC 与 OX 的夹角，坐 标 系 1 内 任 意 点 绕 Y 轴 旋 转 角 度 RY1 = QFWJ( CZ1，CX1，AZ1，AX1) pi0 5; 坐标系 2 内任意点旋转角度为 RY2 = QFWJ( CZ2，CX2，AZ2，AX2) pi0 5，旋转后三角形 ABC 在 XOY 平面上，此时坐标 系 1 与 坐 标 系 2等价( 如图 2 所示) ，坐标系 1 下任意点旋转公式为:z = cos( ry1) z + sin( ry1) xx = sin( ry1) z + cos( ry1) x( 9)( 10)( 3)( 4)则旋转后:AZ1 = cos( ry1) AZ1 + Sin( ry1) AZ1;AX1 = sin( ry1) AZ1 + cos( ry1) AX1;AY1 = AY1;BZ1 = cos( ry1) BZ1 + sin( ry1) BX1;BX1 = sin( ry1) BZ1 + cos( ry1) BX1;BY1 = BY1;CZ1 = cos( ry1) CZ1 + sin( ry1) CX1;CX1 = sin( ry1) CZ1 + cos( ry1) CX1;CY1 = CY1。坐标系 2 下任意点旋转公式为:z = cos( ry2) z + sin( ry2) xx = sin( ry2) z + cos( ry2) x则旋转后:AZ2 = cos( ry2) AZ2 + sin( ry2) AZ2;AX2 = sin( ry2) AZ2 + cos( ry2) AX2;AY2 = AY2;BZ2 = cos( ry2) BZ2 + sin( ry2) BX2;BX2 = sin( ry2) BZ2 + cos( ry2) BX2;BY2 = BY2;CZ2 = cos( ry2) CZ2 + sin( ry2) CX2; CX2 = sin( ry2) CZ2 + cos( ry2) CX2; CY2 = CY2。( 11)( 12)与绕 Z 轴旋转类似，将坐标系统绕 X 轴旋转。其中，设绕 X 轴旋转角度 RX 等于在平面 YOZ 上 AB 与 OY 的夹 角，坐 标 系 1 内 任 意 点 绕 X 轴 旋 转 角 度 RX1 = QFWJ( BY1，BZ1，AY1，AZ1) pi0 5; 坐 标 系 2 内 任 意 点 旋 转角度为 RX2 = QFWJ( BY2，BZ2，AY2，AZ2) pi0 5，旋转后 AB 与 OY 轴重合，坐标系 1 下任意点旋转公式为:y = Math cos( rx1) y + Math sin( rx1) zz = Math sin( rx1) y + Math cos( rx1) z则旋转后:AY1 = cos( rx1) AY1 + sin( rx1) AZ1;AZ1 = sin( rx1) AY1 + cos( rx1) AZ1;AX1 = AX1;BY1 = cos( rx1) BY1 + sin( rx1) BZ1;BZ1 = sin( rx1) BY1 + cos( rx1) BZ1;( 5)( 6)3变换的矩阵表示和相关算法描述在经过 1，2，3 变换后，坐 标 系 1 与 坐 标 系 2 等 价，如BX1 = BX1;果需要进行坐标系 1 到坐标系 2 的变换，只需对坐标系 1在 1，2，3 的变换基础上，进行坐标系 2 的 3，2，1 的逆变换即可。具体变换如下: 将坐标系 1 围绕 Y 轴旋转 ry2;将坐标系 1 围绕 X 轴旋转 rx2; 将 坐 标 系 1 围 绕 Z轴旋转 rz2; 将 坐 标 系 1 在 XYZ 上 分 别 平 移 MX， MY， MZ。坐标系 2 变换为 坐 标 系 1 与 坐 标 系 1 到 坐 标系 2 的变换原理相同，不再赘述。旋转参数。x2，y2，z2 为点 ( x1，y1，z1) 在坐标系 2 内 转 换后的坐标。图 6 3 维直角坐标系转换矩阵Fig 6 The conversion matrix of 3 dimensionalCartesian coordinate system3 1算法的矩阵表示如图 6 所示。其中，x1，y1，z1 为任意点在坐标系 1 内的坐标值; s 为坐标系 1 到坐标系 2 的缩放系数; R1，R2 为图 7 正旋转矩阵Fig 7 Positive rotation matrix图 8 逆旋转矩阵Fig 8 Negative rotation matrix其中，x1，ry1，rz1，rx2，ry2，rz2 为坐标系 1 与坐标系 2 变换为中间坐标系的旋转参数。3 2 求方位角算法描述public double qfwj ( double x1，double y1，double x2，doubley2)double jd = 0 0;if( x1 x2) ( y1y2) )jd = Math Atn( ( y1 y2) / ( x1 x2) ) ;if( ( x1 x2) ( y1 y2) )jd = Math PI* 2 0 + Math Atan( ( y1 y2) / ( x1 x2) ) ;if( ( x1 = = x2) ( y1 y2) )jd = Math PI /2 0;if( x1 = = x2) ( y1 y2) )jd = Math PI* 1 5; |if( x1 x2)jd = Math PI + Math Atan ( ( y1 y2 ) / ( x1 x2) ) ;return jd;4算法验证在 CAD 下构造十个点( 其中点 1，2，3 不 在 一 条 直 线上，用来反算直接转换参数; 点 4 10 用来验 证 转 换 的 正确性) 坐标，其旋转前后及转换后坐标见表 1。表 1 3 维十六参数坐标转换检测表( 单位: m)Tab 1 The verification table of 3 dimensional 16 parameter Cartesian coordinate system conversion( unit: m)原始坐标旋转结果用直接参数法推算结果点号XYZXYZXYZ123745 058 71364 806 01222 683 7376 442 9240 660 0681 517 3201 500 045 700 0324 900 0958 053 81 590 215 01 391 821 8579 812 3469 400 8963 781 1 27 2141 149 602 1 56 142 3958 053 81 590 215 01 391 821 8579 812 3469 400 8963 781 1 27 214 1 149 602 1 56 142 3( 下转第 182 页)测绘与空间地理信息2012 年182视线上下最近 的 两 个 刻 划，垂直角的变化可 以 通 过 以 下公式算出:水准测量的误差来源很多，很 广 泛。各种环境因素和客 观条件复杂多 样，外 界 环 境 千 变 万 化，有 的 影 响 大，有 的 影响小。但无论如何，为了提高观测精度，在实际 操 作 中 应灵活掌握，考虑周全，合理调度，正确运用仪器、克 服 不 利条件。在国家 一、二 等 水 准 测 量 规 范 中 严 格 规 定 跨 河水准的多次读数，多组数，多测回观测。其目的 就 是 要 应用概率论的 原 理，通过大量的观测值取平 均 数 的 方 法 消除各种误差带来的影响，使 观测成果接近真值。在 外 业的实际作 业 中，应 严 格 按 照 规 范 执 行。在 呈 像 条 件 较 差的情况下，更应该增加观测光段数、测回数及组数来提 高观测精度。 = h D( 11)式中: = 206 265; D 为仪器距离标尺的距离( 单位: m) ; h为标尺相邻刻划间的距离( 单位: m) 。由公式可以计算出，当距离为 10 m 时，采用刻划间隔 距离为 2 mm 的 光 学 水 准 标 尺，相邻刻划间的垂直角为 41 25，因此近标尺读数的垂直角最大值即为 41 25。实 验证明，对水平视线附近的两 个刻划线分别测量、计 算， 得到近尺高差最大相差 0 8 mm。由此可 见，水 平 视 线 量 取的结果精度 比 较 高，所以无论选取哪个刻 划 都 是 在 精 度允许范围内的。2 7 测距误差测距的精度在整体短距离跨河水准计算中是可以 忽 略不计的。测距误差的大小是由于测量方法的不同而不 同的，本次测量测距是使用全 站仪的测距功能。测 出 两 点间的斜距及垂直角，在内业进行改正计算。参考文献:1国 家 测 绘 总 局 国家 三角测量与精密 导 线测量 规 范S 北京: 测绘出版社，1975孟庆云，侯和轩 光电三角高程在铁路既有线上的试验 和精度分析J 兰州铁道学院学报，1995徐晖 论三 角高程测量中的大气折光影 响J 工 程 勘 察，1993( 6) : 34 37编辑: 吴 迪233结束语通过以上的论述可知，一、二等测距三角高 程 法 跨 河( 上接第 178 页)续表 1Tab 1( Continued)原始坐标旋转结果用直接参数法推算结果点号XYZXYZXYZ4946 399 7614 246 1134 500 01 126 668 1794 8