logo资料库

坐标系转换公式:各种地理坐标相互转换.pdf

第1页 / 共25页
第2页 / 共25页
第3页 / 共25页
第4页 / 共25页
第5页 / 共25页
第6页 / 共25页
第7页 / 共25页
第8页 / 共25页
资料共25页,剩余部分请下载后查看
坐标系转换公式 青岛海洋地质研究所 戴勤奋 译 (Email: qddqinfen@cgs.gov.cn) 本 文 译 自 国 际 石 油 技 术 软 件 开 放 公 司 ( http://www.posc.org ) 文 献 “ Coordinate Conversions and Transformations including Formulas ”part 2 ––“Formulas for Coordinate Operations other than Map Projections”。原文献由 EPSG(欧洲石油勘探组织)编写,初稿时 间 1995,最新更新时间 2004 年 1 月。世界各国参考椭球体参数及基准面转换参数,可参考 (http://www.globalsecurity.org/military/library/policy/army/fm/6-2/appe2.htm#tabe_22)。 1.前言 坐标系转换问题在工作中经常会遇到,例如,在陆地和海洋的地震勘探中,当今最便捷 的定位方法是 GPS 卫星定位,可是 GPS 定位数据是 WGS84 坐标系上的数据,而各国采用 的往往是早先建立的国家坐标系,为了避免出现矛盾,就需要将 WGS84 定位数据转换到国 家坐标系上。 如果将实测定位坐标或说原大地坐标用 XYZ 直角坐标(地心坐标系)来表示,而不是 用经纬度,那么坐标系转换就会容易得多。然而,欲获取地球表面任一点的 XYZ 值,必须 已知该点相对大地坐标系椭球面的高度,当然你可以假设该点位于椭球面上,高度为零,但 这种情况在实际中很少见。因此,必须考虑坐标点的高度,且记:是相对椭球面的高度。坐 标点相对 WGS84 椭球面的高度可以通过 GPS 测量得到,但相对其它大地坐标系椭球面的 高度却不容易直接得到。大地测量获取的通常是与重力相关的高度值,也就是相对国家高程 基准的高度,高程基准是特定海区的平均海平面,由长期观测数据的综合推算而得,采用传 统方法联测建立的水准点都是以该大地水准面为基准的。因此,如果大地水准面相对椭球面 的高度已知,坐标点相对椭球面的高度就不难得到了。可是在全球大部分地区,大地水准面 相对椭球面的高度数据往往精度不足。不过全球、及全球部分地区与国家已经建立了大地水 准面的数学模型,随着卫星和地面重力数据的不断积累,大地水准面模型的精度正在不断提 高。通过高精度大地水准面数据,就可将实测高程转换成相对椭球面的高度,实现坐标系的 转换。 不过,如无法得到坐标点相对椭球面的高度,可以假设高度值为零,一般不会导致水平 坐标值产生太大的偏差。 在卫星探测早期阶段,由于大地坐标系之间的关系还未能明确定义,并且数据本身的精 度也不高,通常采用 dX,dY,dZ(两椭球参心差值)三参数法进行坐标系转换。该方法假 定两个大地坐标系的直角坐标轴相互平行,当然这种假设通常是不成立的。对一个国家或地 域的局部地区来说,该假设引起的误差可以忽略,一般小于数据的观测精度。然而,随着数 据的不断积累、认知的不断深入及探测方法精度的不断提高,人们逐渐发现,当精度要求较 高时,三参数转换法既不适合在全球范围应用,也不适合地域广大的国家与地区应用。 对石油勘探而言,在特定的勘探许可区内,三参数转换法也许完全能满足精度要求,但 并不能由此假定此区块的转换参数就一定适合邻近地区。 最简单的坐标系转换方法就是上述三参数法,通过两坐标系的原点位移实现,莫洛金斯
基(Molodenski)提出了相应三参数的直接转换方法。方法假定原坐标系与新坐标系的坐标 轴相互平行,正如前面已提到的,该假设不一定成立,由此得到的转换结果只能达到中等精 度,对范围大的区域尤其如此。 赫尔默特(Helmert)7 参数转换法提高了转换精度。由于三个旋转参数有两种相反的符 号协定,EPSG(欧洲石油勘探组织)将其分为两种不同的转换方法,其一称位置矢量法, 另一称坐标框架法,其中位置矢量法也称布尔莎-沃尔夫(Bursa-Wolf)转换法。赫尔默特方 法的关键在于转换参数的符号要与遵循的约定一致。鉴于赫尔默特方法的平移和旋转参数之 间具有很强的相关性,有碍于实际应用,莫洛金斯基-巴德克斯(Molodenski-Badekas)提出 了改进的赫尔默特 7 参数转换法,避免了上述相关性问题。 根据研究区内一系列已知点的大地坐标或网格坐标改正量进行插值,也是一种坐标系转 换方法。北美 1927 基准面(基于 Clarke 1966 椭球体)与北美 1983 基准面(基于 GRS 1980 椭球体)之间的坐标系转换就是其中一例。北美测量控制网是用传统大地测量方法建立的, 由于早期的仪器精度不足、网平差不完善等因素,基于 Clarke 1966 椭球体、并且只有一个 基准点(位于堪萨斯州 Meades Ranch)的老坐标网精度低、且误差分布不均匀;新坐标网 采用了卫星技术、现代先进的测量仪器和电子计算机技术,其精确度与可靠性完全能得到保 障。由此造成北美大陆网内,不同地区、甚至不同位置点的转换参数都有可能不一致,所以 如仅采用莫洛金斯基(Molodenski)和赫尔默特(Helmert)方法对付上述新、旧坐标系的 转换显然不合适,为此需要用到 EPSG(欧洲石油勘探组织)所谓的“双线性插值”转换技 术。到北美 NAD83 的坐标转换就是通过格网双线性插值实现的,其中采用了美国海岸带与 大地测绘局(US Coast & Geodetic Survey)的 NADCON 控制点网。注:美国以西经为正, 而 EPSG 文献中 NAD27 与 NAD83 坐标系的设定均以东经为正;加拿大的网格文件格式也 被澳大利亚与新西兰采用;英国采用北向与东向的双线性网格插值。 此外,经纬度多项式也可以用于坐标系转换,挪威在海岸带调查中,就采用这种方法进 行新(ED87— 欧洲 1987 基准面)、旧(ED50— 欧洲 1950 基准面)坐标系之间的转换。挪 威地调局 Statens Kartwerk 发表的文献中列出了包含 15 个系数的经纬度 4 次多项式展开公 式。 上述格网插值及多项式拟合方法更适合于早期基准面与新建基准面之间的坐标值转换。 坐标系转换时,选择正确的转换参数符号非常重要,转换前应该明确转换的“从 (From)”……“到(To)” ……,避免符号混淆。 2.坐标变换方法 2.1 大地坐标与地心坐标的变换 大地经纬度坐标(纬度j ,经度 l )可以用地心直角坐标 X、Y、Z 表示,其中,直角 坐标系原点位于地心;Z 轴为极轴,向北为正;X 轴穿过本初子午线与赤道的交点;Y 轴穿 过赤道与东经 90°的交点。 本文设定坐标系的零经线为格林威治子午线,如果定义不一致,在使用各公式前首先将
零经线转换到格林威治子午线。 设椭球长半轴为 a ,短半轴为 b ,扁率倒数为 1/f ,那么 X Y Z = nj l ( +h)coscos = nj l +h)cos (sin = n ) (( 1-e+h)sin 2 j 式中:n 为纬度j 处的卯酉圈曲率半径, n = a 2 /(1sin e ) 20.5 j j 和 l 分别为坐标点的纬度和经度, h 为相对椭球面的高度, e 为椭球第一偏心率, 22 ea =-= 22 ()/ baf 2 f 2 (注意:h 为相对椭球面的高度,也就是通过 GPS 卫星定位就可观测到的高度值,而不是通 常的与重力相关的大地测量高程值。重力相关的高程(H)通常是相对海平面,或某一水准面 的高度。如果重力高程 H 已知,那么在使用以上公式时必须将其转换成椭球高程 h。h = H + N,其中 N 为大地水准面相对椭球面的高度,N 有时为负值。大地水准面是近视于海平面的 重力面。WGS84 椭球的 N 值在 -100 米(斯里兰卡)到 +60(北大西洋)之间。不过国家坐 标系的椭球面与大地水准面的相对高度一般不容易得到。) 反之,X、Y、Z 地心直角坐标也可以转换成经纬度大地坐标: =+ jn = l = h X 12220.5 j [(sin)/( Y tgZeX 1 ( ) tg Y X / lj n secsec + ) ] 通过迭代计算 式中 l 从格林威治本初子午线起算。 实例请见 3.2。 3.坐标系转换方法 3.1 偏移量 数学上,如果一维坐标系的原点沿轴的正向移动距离 A(A > 0),那么转换公式为: X = newold X A 然而,在坐标系转换中,通常将偏移量视作校正值,符号与数学上相反,EPSG 自 1999 年 起采用该协定至今。这样,原坐标值加上一个校正参数就得到新坐标值: X t = + A X s 其中,Xs 和 Xt 分别为原坐标系和新坐标系的坐标值;A 为转换参数,通过该参数将坐标值 从原坐标系转换到新坐标系。 3.1.1 垂向偏移量(EPSG 坐标运算方法 9616) - - - - - -
如果两坐标系的坐标轴向定义不一致,例如将原坐标系中的高度(向上为正),转换为 新坐标系的深度(向下为正),或转换成不同的坐标单位,公式 Xt = Xs + A 需要调整为: = tss XXUAUm U *) A [( + (*)]*( / t ) 式中:Xt = 新坐标系的坐标值, Xs = 原坐标系的坐标值, A 为新坐标系的原点在原坐标系中的坐标值, m 为单位方向因子(如果两坐标系轴向一致 m=1;如果两坐标系轴向不一致 m = -1), Us、Ut 和 UA 分别为原坐标、新坐标和偏移参数的坐标单位与米的比值。 3.2 三参数转换(EPSG 坐标运算方法 9603) 假设两椭球体的长、短轴相互平行,零经线为格林威治本初子午线,从原坐标系转换到 新坐标系的三平移参数为 dX,dY,dZ,那么转换公式为: t X Y t Z t = + XdX s + = YdY s = + ZdZ s 实例: 本例结合了 2.1 中所述的大地坐标与地心坐标的转换方法。 已知位于北海的一个 WGS84 大地坐标点,由 GPS 卫星定位: s j 纬度 l 经度 s 椭球高 = 5348'33.82''N = 207'46.38''E = sh 73.0 m 需要将其转换为 ED50 大地坐标,相应的椭球体为 International 1924。北海区从 WGS84 转 换到 ED50 的转换三参数为:dX = +84.87m ,dY = +96.49m,dZ = +116.95m。 首先将 WGS84 大地坐标转换为地心直角坐标: s X Y s Z s = 3771793.97 = 140253.34 m = 5124304.35 m m 根据上述地心坐标转换方法,得到 ED50 的地心直角坐标: t X Y t Z t = = =+ 3771793.9784.873771878.84 =+ 140253.3496.49140349.83 =+ 5124304.35116.955124421.30 = m m m 利用 2 中定义的反变换方法,得到 ED50 的大地坐标:
t j 纬度 l 经度 t 椭球高 = 5348'36.565''N = 207'51.477''E = th 28.02 m 其中椭球高从 International 1924 椭球面起算,如果换算到海平面高程需要作大地水准面高度 校正。 3.3 简化莫洛金斯基(Molodenski)转换 莫洛金斯基(Molodenski)推出的转换公式,可将上述三参数方法的计算步骤合而为一, 公式的简化形式非常适合三参数坐标系转换: = jj t = ll t = h t s + j + l + hdh s s d d 式中: =- jjljljj r ddXdYdZadffda ''(.sin.cos.sin.sin.cos[..].sin2)/(.sin1 + =- j llln ''(.sin.cos)/(.cos.sin1'') ddXdY =++++ jljlj dhdXdYdZadffdada .cos.cos.cos.sin.sin(..).sin -++ + '') j 2 其中:dX,dY,dZ 为两椭球参心差值,也就是椭球体原点平移参数 r 为原椭球体纬度j 处的子午圈曲率半径 =- r (1)/(1sin ae e 22 23/2 ) j n 为原椭球体纬度j 处的卯酉圈曲率半径 n = a 221/2 ) /(1sin e j da 为新椭球体与原椭球体的长半轴之差 = da a t a s df 为新椭球体与原椭球体的扁率之差 df = -= ff tst 1/(1/ s f )1/(1/ f ) djl j l 和d 为 和 的偏差值,以弧度为单位。 实例: 已知位于北海的一个 WGS84 大地坐标点,由 GPS 卫星定位: s j 纬度 l 经度 s 椭球高 = 5348'33.82''N = 207'46.38''E = sh 73.0 m 需要将其转换为 ED50 大地坐标,相应的椭球体为 International 1924。 北海区从 WGS84 转换到 ED50 的转换三参数为: dX = +84.87m dY = +96.49m dZ = +116.95m 椭球参数为: - - - - -
WGS 1984 a = 6378137.0 米 1/f = 298.2572236 International 1924 a = 6378388.0 米 1/f = 297.0 计算得到: da = 6378137 - 6378388 = -251 df = 0.003352811 – 0.003367003 = -1.41927E-05 计算得到: = = = - j d l d dh 2.545'' 5.097'' 44.98 m 计算得到 ED50(基于 International 1924 椭球体)的大地坐标值为: j l t 纬度 经度 椭球高h = 5348'36.565'' = 207'51.477'' = 28.02 m t t N E 3.4 赫尔默特(Helmert)转换 从一个大地坐标系转换到另一个大地坐标系(俗称为基准面转换)一般需要经过三个环 节:大地坐标到地心坐标 → 地心坐标到地心坐标 → 地心坐标到大地坐标。 其中的中间环节,地心坐标到地心坐标,通常又称为 7 参数赫尔默特(Helmert)转换, 将转换公式用 7 参数矩阵表示,即得到著名的布尔莎-沃尔夫(Bursa-Wolf)公式: T X Y TZX Z T 1 =+ MRRYdY * S + RRdZ Y Z + RRdX 1 1 X Y + * X S Z S 式中,(Xs,Ys,Zs)为原坐标系中的点坐标,(Xt,Yt,Zt)为新坐标系中的点坐标。 转换参数的定义不唯一, 可引伸出不同的转换方法,其中“位置矢量转换”方法(EPSG 坐标运算方法 9606)在欧洲石油勘探行业中广泛采用,国际大地测量协会(IAG)也通过 ISO 19111 建议采用该方法的参数定义: (dX,dY,dZ):两坐标系的原点平移矢量(平移参数),原坐标系中的点位置矢量加上原点 平移矢量即得到该点在新坐标系中的位置矢量。平移参数也就是原坐标系的原点在新坐标系 中的坐标值。 (RX, RY, RZ):位置矢量的旋转角(旋转参数)。参数符号约定如下:从直角坐标系原点, 沿轴正向看,位置矢量绕轴顺时针旋转为正。从原坐标系转换到新坐标系,如果绕 Z 轴的旋 转角度为正,那么转换后坐标点的经度将增大。公式中的角度单位,本文要求是弧度。 M:位置矢量的比例因子(尺度比参数),位置矢量从原坐标系转换到新坐标系的尺度伸缩 量。M = (1 + dS*10-6),其中 dS 为尺度校正量,以百万分之一计(ppm)。 实例: - Ø ø Ø ø Ø ø Œ œ Œ œ Œ œ - Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ - Ł ł º ß º ß º ß
本例结合了 2.1 中所述的大地坐标与地心坐标的转换方法及位置矢量转换方法。 要求从 WGS72 转换到 WGS84,转换参数如下: dX = 0.000 m dY = 0.000 m dZ = +4.5 m RX = 0.000″ RY = 0.000″ RZ = +0.554″ dS = +0.219 ppm 已知 WGS72 大地坐标: S j 纬度 l 经度 S 椭球高 = 5500'00''N = 400'00''E = Sh m 0 利用大地坐标到地心坐标的转换公式得到: S m X Y S Z = 3657660.66 = 255768.55 m = 5201382.11 m 利用 7 参数位置矢量转换方法得到: m S T X Y T Z T = 3657660.78 = 255778.43 m = 5201387.75 m 利用地心坐标到大地坐标的转换公式得到点在 WGS84 上的三维坐标值: j 纬度 T l 经度 T 椭球高 = = Th 5500'00.090''N 400'00.554''E = + 3.22 m 位置矢量转换法虽然在欧洲石油勘探行业中普遍应用,但其符号约定并未得到全球的 认同。美国石油勘探行业则采用另一种符号约定,其区别在于旋转参数的约定上,称作“坐 标框架旋转”约定(EPSG 坐标运算方法 9607),转换公式为: T X Y TZX Z T =- 1 MRRYdY * S + + RRdZ Y + Z RRdX 1 Y + * 1 X X S Z S 转换参数定义为: (dX,dY,dZ):两坐标系的原点平移矢量(平移参数),原坐标系中的点位置矢量加上原点 平移矢量即得到该点在新坐标系中的位置矢量。平移参数也就是原坐标系的原点在新坐标系 中的坐标值。 (RX, RY, RZ):坐标参考框架的旋转角(旋转参数)。参数符号约定如下:从直角坐标系原 - Ø ø Ø ø Ø ø Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ - Ł ł º ß º ß º ß
点,沿轴正向看,坐标参考框架绕轴顺时针旋转为正。从原坐标系转换到新坐标系,如果绕 Z 轴的旋转角度为正,那么转换后坐标点的经度将变小。公式中的角度单位,本文要求是弧 度。 M:位置矢量的比例因子(尺度比参数),位置矢量从原坐标系转换到新坐标系的尺度伸缩 量。M = (1 + dS*10-6),其中 dS 为尺度校正量,以百万分之一计(ppm)。 如果不考虑旋转因素,位置矢量法与坐标框架法是一致的。需要注意的是,同样的旋 转变换,在第一种方法中转换参数为正,而在第二种方法中却是负的。所以,在基准面转换 前,必须明确旋转参数的符号约定,可以说,参数的符号是与坐标系转换方法相关联的。 前面位置矢量法的例子同样可通过坐标框架法计算,只不过输入的转换参数略有差异: dX = 0.000 m dY = 0.000 m dZ = +4.5 m RX = 0.000″ RY = 0.000″ RZ = -0.554″ dS = +0.219 ppm 注意:与位置矢量法相比,只有旋转参数的符号变号;此外,国际大地测量协会(IAG) 在 ISO 19111 标准中建议使用位置矢量法的符号约定。 3.5 莫洛金斯基-巴德卡斯(Molodenski-Badekas)转换 为了消除赫尔默特(Helmert)方法中平移与旋转参数之间的强相关性,引入了另一旋 转中心点,也就是旋转中心由原来的地心坐标系原点,改为一个特定的位置,转换公式变为: + T X Y TZXSP Z T 1 + =-+- MRRYYYdY RRdZ Y P + * RRdX Z + 1 X Y 1 * X S X P Z S Z P X P Z P 参数定义如下: (dX,dY,dZ):两坐标系的原点平移矢量(平移参数),原坐标系中的点位置矢量加上原点 平移矢量即得到该点在新坐标系中的位置矢量。平移参数也就是原坐标系的原点在新坐标系 中的坐标值。 (RX, RY, RZ):坐标参考框架的旋转角(旋转参数)。参数符号约定如下:从直角坐标系原 点,沿轴正向看,坐标参考框架绕轴顺时针旋转为正。从原坐标系转换到新坐标系,如果绕 Z 轴的旋转角度为正,那么转换后坐标点的经度将变小。公式中的角度单位,本文要求是弧 度。 (XP, YP, ZP):坐标参考框架的旋转中心点,在原直角坐标系中定义。 M:位置矢量的比例因子(尺度比参数),位置矢量从原坐标系转换到新坐标系的尺度伸缩 - - Ø ø Ø ø Ø ø Ø ø Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ Œ œ - - Ł ł º ß º ß º ß º ß
分享到:
收藏