遥感 · 2022年 9月 30日 0

空间坐标与投影系统

目录

空间坐标和地图投影

GIS处理的是空间信息，而所有对空间信息的量算都是基于某个坐标系统的，因此GIS中坐标系统的定义是GIS系统的基础，正确理解GIS中的坐标系统就变得尤为重要。坐标系统又可分为两大类：地理坐标系统、投影坐标系统。GIS中的坐标系定义由基准面和地图投影两组参数确定，而基准面的定义则由特定椭球体及其对应的转换参数确定，因此欲正确定义GIS系统坐标系，首先必须弄清地球椭球体(Ellipsoid)、大地基准面(Datum)及地图投影(Projection)三者的基本概念及它们之间的关系。

地球椭球体(Ellipsoid)

地球表面是一个凸凹不平的表面，而对于地球测量而言，地表是一个无法用数学公式表达的曲面，这样的曲面不能作为测量和制图的基准面。假想一个扁率极小的椭圆，绕大地球体短轴旋转所形成的规则椭球体称之为地球椭球体。地球椭球体表面是一个规则的数学表面，可以用数学公式表达，所以在测量和制图中就用它替代地球的自然表面。因此就有了地球椭球体的概念。地球椭球体有长半径和短半径之分，长半径(a)即赤道半径，短半径(b)即极半径。f=（a-b）/a为椭球体的扁率，表示椭球体的扁平程度，a、b、f被称为地球椭球体的三要素。

img

对地球椭球体而言，其围绕旋转的轴叫地轴。地轴的北端称为地球的北极，南端称为南极；过地心与地轴垂直的平面与椭球面的交线是一个圆，这就是地球的赤道；过英国格林威治天文台旧址和地轴的平面与椭球面的交线称为本初子午线。以地球的北极、南极、赤道和本初子午线等作为基本要素，即可构成地球椭球面的地理坐标系统。

img
img

大地基准面（Geodetic datum）

大地基准面（Geodetic datum），设计用为最密合部份或全部大地水准面的数学模式。它由椭球体本身及椭球体和地表上一点视为原点间之关系来定义。此关系能以 6个量来定义，通常（但非必然）是大地纬度、大地经度、原点高度、原点垂线偏差之两分量及原点至某点的大地方位角。

我们把地球椭球体和基准面结合起来看，在此我们把地球比做是“马铃薯”，表面凸凹不平，而地球椭球体就好比一个“鸭蛋”，那么按照我们前面的定义，基准面就定义了怎样拿这个“鸭蛋”去逼近“马铃薯”某一个区域的表面，X、Y、Z轴进行一定的偏移，并各自旋转一定的角度，大小不适当的时候就缩放一下“鸭蛋”，那么通过如上的处理必定可以达到很好的逼近地球某一区域的表面。

因此，从这一点上也可以很好的理解，每个国家或地区均有各自的基准面，我们通常称谓的北京54坐标系、西安80坐标系实际上指的是我国的两个大地基准面。我国参照前苏联从1953年起采用克拉索夫斯基(Krassovsky)椭球体建立了我国的北京54坐标系，1978年采用国际大地测量协会推荐的1975地球椭球体（IAG75）建立了我国新的大地坐标系–西安80坐标系，目前大地测量基本上仍以北京54坐标系作为参照，北京54与西安80坐标之间的转换可查阅国家测绘局公布的对照表。 WGS1984基准面采用WGS84椭球体，它是一地心坐标系，即以地心作为椭球体中心，目前GPS测量数据多以WGS1984为基准。

椭球体与基准面之间的关系是一对多的关系，也就是基准面是在椭球体基础上建立的，但椭球体不能代表基准面，同样的椭球体能定义不同的基准面。地球椭球体和基准面之间的关系以及基准面是如何结合地球椭球体从而实现来逼近地球表面的可以通过下图一目了然。

img

投影坐标系统（Projected Coordinate Systems ）

地球椭球体表面也是个曲面，而我们日常生活中的地图及量测空间通常是二维平面，因此在地图制图和线性量测时首先要考虑把曲面转化成平面。由于球面上任何一点的位置是用地理坐标（λ，φ）表示的，而平面上的点的位置是用直角坐标（χ，у）或极坐标（r）表示的，所以要想将地球表面上的点转移到平面上，必须采用一定的方法来确定地理坐标与平面直角坐标或极坐标之间的关系。这种在球面和平面之间建立点与点之间函数关系的数学方法，就是地图投影方法。

每一个投影坐标系统都必定会有Geographic Coordinate System（地理坐标系统）。那么我们从这一角度上解释一下投影和投影所需要的必要条件：将球面坐标转化为平面坐标的过程便是投影过程；投影所需要的必要条件是：第一、任何一种投影都必须基于一个椭球（地球椭球体），第二、将球面坐标转换为平面坐标的过程（投影算法）。简单的说投影坐标系是地理坐标系+投影过程。

我们从透视法（地图投影方法的一种）角度来直观的理解投影, 几何透视法是利用透视的关系，将地球体面上的点投影到投影面（借助的几何面）上的一种投影方法。如假设地球按比例缩小成一个透明的地球仪般的球体，在其球心或球面、球外安置一个光源，将球面上的经纬线投影到球外的一个投影平面上。
img

可以按照变形性质将投影方法如下分类：等角投影（Conformal Projection） 、 等积投影（Equal Area Projection）、等距投影（Equidistant Projection）、等方位投影（True-direction Projection）四种。每种投影根据其名称就可以知道其方法保证了数据的那些几何属性，在实际应用过程中应根据需求来选取某种投影。

接下来我们来看看我们国家通常采用的投影——高斯—克吕格（Gauss-Kruger)投影，是一种“等角横切圆柱投影”。 设想用一个圆柱横切于球面上投影带的中央经线，按照投影带中央经线投影为直线且长度不变和赤道投影为直线的条件，将中央经线两侧一定经差范围内的球面正形投影于圆柱面。然后将圆柱面沿过南北极的母线剪开展平，即获高斯-克吕格投影平面。高斯—克吕格投影后，除中央经线和赤道为直线外，其他经线均为对称于中央经线的曲线。高斯—克吕格投影没有角度变形，在长度和面积上变形也很小，中央经线无变形，自中央经线向投影带边缘，变形逐渐增加，变形最大处在投影带内赤道的两端。按一定经差将地球椭球面划分成若干投影带,这是高斯投影中限制长度变形的最有效方法。

img
img

高斯－克吕格投影以6度或3度分带，每一个分带构成一个独立的平面直角坐标网，投影带中央经线投影后的直线为X轴（纵轴，纬度方向），赤道投影后为Y轴（横轴，经度方向），为了防止经度方向的坐标出现负值，规定每带的中央经线西移500公里，即东伪偏移值为500公里，由于高斯-克吕格投影每一个投影带的坐标都是对本带坐标原点的相对值，所以各带的坐标完全相同，因此规定在横轴坐标前加上带号，如(4231898,21655933)其中21即为带号，同样所定义的东伪偏移值也需要加上带号，如21带的东伪偏移值为21500000米。

我国目前常用的坐标系有北京54、西安80、WGS-84和CGCS2000四种坐标系，接下来详细介绍这几种坐标系及其相互之间的转换方式。

国内常用投影坐标系

我国大中比例尺地图均采用高斯-克吕格投影，其通常是按6度和3度分带投影，1:2.5万－1:50万比例尺地形图采用经差6度分带，1:1万比例尺的地形图采用经差3度分带。具体分带法是：6度分带从本初子午线开始，按经差6度为一个投影带自西向东划分，全球共分60个投影带，带号分别为1－60；3度投影带是从东经1度30秒经线开始，按经差3度为一个投影带自西向东划分，全球共分120个投影带。

为了便于地形图的测量作业，在高斯-克吕格投影带内布置了平面直角坐标系统，具体方法是，规定中央经线为X轴，赤道为Y轴，中央经线与赤道交点为坐标原点，x值在北半球为正，南半球为负，y值在中央经线以东为正，中央经线以西为负。由于我国疆域均在北半球，x值均为正值，为了避免y值出现负值，规定各投影带的坐标纵轴均西移500km，中央经线上原横坐标值由0变为500km。为了方便带间点位的区分，可以在每个点位横坐标y值的百千米位数前加上所在带号，如20带内A点的坐标可以表示为YA=20745921.8m。

北京54坐标系(BJZ54)

北京54坐标系为参心坐标系1，大地上的一点可用经度L54、纬度M54和大地高H54定位，它是以克拉索夫斯基椭球为基础，经局部平差后产生的坐标系。北京54坐标系采用了前苏联的克拉索夫斯基椭球参数，并与前苏联1942年坐标系进行联测，通过计算建立了我国大地坐标系， 1954年北京坐标系可以认为是前苏联1942年坐标系的延伸。它的原点不在北京而是在前苏联的普尔科沃。

北京54坐标系，属参心坐标系，长轴6378245m，短轴6356863，扁率1/298.3。

西安80坐标系

1978年4月在西安召开全国天文大地网平差会议，确定重新定位，建立我国新的坐标系。为此有了1980年国家大地坐标系。1980年国家大地坐标系采用地球椭球基本参数为1975年国际大地测量与地球物理联合会第十六届大会推荐的数据，即IAG75地球椭球体。该坐标系的大地原点设在我国中部的陕西省泾阳县永乐镇，位于西安市西北方向约60公里，故称1980年西安坐标系，又简称西安大地原点。基准面采用青岛大港验潮站1952－1979年确定的黄海平均海水面（即1985国家高程基准）。
西安80坐标系，属参心坐标系1，长轴6378140m，短轴6356755，扁率1/298.25722101。

2000国家大地坐标系

英文缩写为CGCS2000。2000国家大地坐标系是全球地心坐标系1 在我国的具体体现，其原点为包括海洋和大气的整个地球的质量中心。

2000国家大地坐标系采用的地球椭球参数如下：长半轴a=6378137m，扁率f=1/298.257222101。

WGS84坐标系

WGS－84坐标系（World Geodetic System）是一种国际上采用的地心坐标系1。坐标原点为地球质心，其地心空间直角坐标系的Z轴指向国际时间局（BIH）1984.0定义的协议地极（CTP）方向，X轴指向BIH1984.0的协议子午面和CTP赤道的交点，Y轴与Z轴、X轴垂直构成右手坐标系，称为1984年世界大地坐标系。这是一个国际协议地球参考系统（ITRS），是目前国际上统一采用的大地坐标系。GPS广播星历是以WGS-84坐标系为根据的。
WGS84坐标系，长轴6378137.000m，短轴6356752.314，扁率1/298.257223563。 由于采用的椭球基准不一样，并且由于投影的局限性，使的全国各地并不存在一至的转换参数。对于这种转换由于量较大，有条件的话，一般都采用GPS联测已知点，应用GPS软件自动完成坐标的转换。

通用GIS软件中的常用坐标系投影方式

在ArcGIS中，关于坐标系的区别：
在ArcGIS的Coordinate Systems目录下，分为Geographic Coordinate Systems 和 Projected Coordinate Systems，如下：

img

分别表示地理坐标系统和投影坐标系统。

以西安80为例，坐标系统文件命名方式有以下几种：

名称EPSG代码2解释
GCS_Xian_19804610西安80地理坐标系
Xian 1980 3 Degree GK CM 75E2370三度分带法的西安80坐标系，投影方式为高斯-克吕格投影，中央经线在东经75度，横坐标前加带号。
Xian 1980 3 Degree GK CM 75E2349三度分带法的西安80坐标系，投影方式为高斯-克吕格投影，横坐标前加带号。
Xian 1980 GK CM 75E2338六度分带法的西安80坐标系，投影方式为高斯-克吕格投影，中央经线为东经75度，横坐标前不加带号。
Xian 1980 GK Zone 132327六度分带法的西安80坐标系，投影方式为高斯-克吕格投影，横坐标前加带号。

同样的，我们可以以这种方法来区分北京54坐标系和CGCS2000在软件中的表示。
接下来我们以实际的地理位置（以武汉大学珞珈山为例），查看不同投影坐标系的显示方式：

img
坐标系东坐标西坐标
Xian 1980 3 Degree GK CM 114E3379772.618534945.1555
Xian 1980 3 Degree GK Zone 383379772.61838534945.16
Xian 1980 GK CM 117E3382664.038247012.8279
Xian 1980 GK Zone 203382664.03820247012.83

上表可以看出同一个地理位置，表示方法是不一样的，大家平时需要注意这个问题。

通过坐标推断投影带类型

在工作中，我们可能会遇到将一些地图文件(如CAD文件)导入到ArcGIS中，现在我们可以根据上面的规则来推断它的类型。以上文中的位置为例，地图上某个点的坐标为：

3382664.038，20247012.83

这时我们可以检查它的横坐标，为20247012.83，显而易见横坐标前面有带号且带号为20。西安80的三度带范围为25至45，六度带范围为13至23。因此我们可以判断此地图文件的投影带为六度带，且带号为20，并可进一步推出它的中央经线为东经117度。

我们再来看一个坐标点：

3379772.618，534945.1555

它的横坐标前是不加带号的，此时我们可以通过比例尺来判断，我国规定，一般对于1:2.5万～1:10万的地图采用六度带，对于1:1万或更大比例尺的地形图采用三度带。因此这类地图文件可以检查比例尺大小来判断地图是六度带还是三度带。

1 : 参心坐标系、地心坐标系:物体均有其质心，处处密度相等的物体的质心在其几何中心。所以，地球只有一个质心，只是测不测的精确的问题而已。由地球的唯一性和客观存在，以地球质心为旋转椭球面的中心的坐标系，叫地心坐标系，且唯一。当然，由于a、b两个值的不同，就有多种表达方式，例如，CGCS2000系，WGS84系等。（注：地心坐标系又名协议地球坐标系，与GPS中的瞬时地球坐标系要对应起来。）但是又有一个问题——政治问题，地图是给一个国家服务的，那么这地图就要尽可能描述准确这个国家的地形地貌，尽量减小误差，至于别国就无所谓。所以，就可以人为的把地球的质心“移走”，将局部的表面“贴到”该国的国土，使之高程误差尽量减小到最小。这个时候，就出现了所谓的“参心坐标系”。即椭球中心不在地球质心的坐标系。

我国常用的参心系及对应椭球：

  • 北京54坐标系：克拉索夫斯基椭球体
  • 西安80坐标系：IAG75椭球体

我国常用的地心系及对应椭球：

  • WGS84坐标系：WGS84椭球体
  • CGCS2000坐标系：CGCS2000椭球体（事实上，CGCS2000椭球和WGS84椭球极为相似，偏差仅有0.11mm，完全可以兼容使用）

2 : European Petroleum Survey Group (EPSG)，负责维护并发布坐标参照系统的数据集参数，以及坐标转换描述，该数据集被广泛接受并使用。目前已有的椭球体，投影坐标系等不同组合都对应着不同的ID号，这个号在EPSG中被称为EPSG code，它代表特定的椭球体、单位、地理坐标系或投影坐标系等信息。

坐标系的选用在测绘工作中十分重要，当然熟悉坐标系统知识也是我们测绘工作者的基本要求。近几年，随着测绘科学的发展，3S技术相继引入，作业难度降低，测绘精度逐步提高。国内测绘部门先后采用了几种坐标系统，为了充分利用现有测绘数据，节约成本，我们需要掌握坐标转换方法，接下来，我们看一下我国常用坐标系的转换方法。

坐标转换

坐标系统是测量工作中定位的基础，坐标系统有很多形式和基准，不同历史时期所建立和使用的坐标系是不同的。随着科学技术的进步，测量方法和观测技术不断改进，采用的参考椭球及定位方式也逐步完善和精化。为更加精确的确定点位信息并综合利用各种信息，常常用到坐标转换。

坐标转换通常包含两层含义：基准变换和坐标系变换。

坐标系变换

坐标系变换就是在同一地球椭球下，空间点的不同坐标表示形式进行变换。包括大地坐标系与空间直角坐标系的相互转换、空间直角坐标系与质心坐标系的转换以及大地坐标系与高斯平面坐标系的转换。

基准变换

基准变换是指空间点在不同椭球的坐标转换。可用空间的三参数或者七参数实现不同椭球间空间直角坐标系或者不同椭球间大地坐标系的转换。

坐标转换的基本方法

坐标转换的基本方法分为大地坐标与空间直角坐标的转换、大地坐标与高斯平面直角坐标的转换，以及不同大地坐标系之间的转换。其中大地坐标与高斯平面直角坐标的转换是大地坐标(B，L)向高斯平面直角坐标(X，Y)的转换，称为高斯正算。由高斯平面直角坐标(X， Y)向大地坐标(B，L)的转换，称为高斯反算。高斯正反算只能在同一参考椭球下进行。

对于不同大地坐标系之间的转换则需要注意：在不同的坐标系统之间, 由于椭球参数不同, 两个椭球之间没有一种统一的方法实现坐标转换。但是, 在两个椭球所指的同一区域内, 由于椭球面弯曲度较小, 该区域同名点在不同的椭球系上存在一定的曲面数学关系, 因此可以通过区域转换模型进行坐标转换。一般常用的转换方法是四参数转换法和七参数转换法。

全国及省级范围的坐标转换选择二维七参数转换模型;省级以下的坐标转换可选择三维四参数模型或平面四参数模型。对于相对独立的平面坐标系统与2000国家大地坐标系的转换可采用平面四参数模型或多项式回归模型。但是最通用的方法是布尔莎七参数转换法, 也称综合转换, 所谓综合法即就是在相似变换(布尔莎七参数转换)的基础上, 再对空间直角坐标残差进行多项式拟合, 系统误差通过多项式系数得到消弱, 使统一后的坐标系框架点坐标具有较好的一致性, 从而提高坐标转换精度。

布尔莎七参数转换法涉及到的七个参数为：X平移，Y平移，Z平移，X旋转，Y旋转，Z旋转，尺度变化m。此七个参数可以通过在需要转化的区域里选取3个以上的转换控制点对而获取。

如果区域范围不大，最远点间的距离不大于30Km，可以用三参数（莫洛登斯基模型），即X平移，Y平移，Z平移，而此时将X旋转，Y旋转，Z旋转，尺度变化 m 视为0。所以三参数只是七参数的一种特例。三参数只需通过1个控制点对就能获取。

七参数公式如下:

img
\begin{bmatrix}
 X\\
 Y\\
Z
\end{bmatrix}
=
\begin{bmatrix}
 \Delta X\\
 \Delta Y\\
\Delta Z
\end{bmatrix}
+(1+m)\begin{bmatrix}
  1&  \theta_z& -\theta_y\\
  -\theta_z & 1 &\theta_x \\
  \theta_y&  -\theta_x& 1
\end{bmatrix}
\begin{bmatrix}
 X_0\\
 Y_0\\
Z_0
\end{bmatrix}

坐标系转换的实施步骤

  1. 收集、整理转换区域内重合点成果。
  2. 分析、选取用于计算坐标转换参数的重合点。
  3. 确定坐标转换参数计算方法与坐标转换模型。
  4. 两坐标系下重合点坐标形式的转换。若采用平面四参数转换模型，则要将重合点的两坐标系坐标转换成同一投影带的高斯平面坐标（至少需要两对重合点坐标）；若采用布尔莎七参数转换模型（至少需要三对重合点坐标），则要将重合点的两坐标系坐标转换成个坐标系下的空间直角坐标。
  5. 根据确定的转换方法与转换模型利用最小二乘法初步计算坐标转换参数。
  6. 分析重合点坐标转换残差，根据转换残差剔除粗差点。
    坐标转换残差满足精度要求时，计算最终的坐标转换参数并估计坐标转换精度。
  7. 根据计算得到的转换参数，转换得到转换点的目标坐标系坐标。

坐标转换常见误区

坐标转换就是旋转平移

坐标转换不仅仅是旋转平移，因为椭球参数不一致，以及地球表面不是平面，坐标转换往往涉及到旋转、平移、缩放、比例因子等问题。我们所说的“旋转平移”有点接近平面四参数转换法，它只是一种近似的转换方法，只适用于小范围、精度要求不高的平面坐标转换。在使用该法转前，要先将转换点的原坐标系坐标换算为计算转换参数时所在投影带(中央子午线)的高斯平面坐标(X，Y)。再根据两个以上重合点计算出转换参数(2个平移参数、1个旋转参数、1个尺度参数)。工作中，一般根据工作区范围和精度要求选择转换方法，通常采用布尔莎七参数法。

地面上同一个点只有一个经纬度坐标

这种认识是错误的。我们所说的经纬度坐标就是大地坐标，地面上同一个点，对于不同的参考椭球，它的大地坐标是不一样的。值得注意的是，2000坐标系的大地坐标是地心经度和地心纬度。地心经度和大地经度定义一样，地心纬度是指地面点的向径(地面点与质心的连线)与赤道面的夹角。其它参考椭球的经度和纬度定义是一样的。由于，椭球参数不一致，因此地面上同一点，在不同的参考椭球下面，大地坐标是不一样的。

坐标转换就是经纬度转换为平面直角坐标系

经纬度坐标转换为平面直角坐标只是坐标转换的一种，坐标转换还包括很多种类型，经纬度坐标与高斯平面直角坐标之间的坐标转换有两种情况：一种是在同一参考椭球下转换，另一种是在不同参考椭球下转换。在同一参考椭球下转换可直接利用高斯正反算进行；在不同参考椭球下转换，首先 要将两套坐标转换为同一参考椭球下坐标，再利用高斯正反算进行转换计算。

坐标转换就是换带计算

坐标转换是一个广泛的概念，它包括很多种类，涉及到不同的椭球参数。而换带计算只是坐标换算的一种，它是指利用高斯正反算，将地面点在特定中央子午线经度和带宽下的高斯平面直角坐标，换算为相邻投影带或任意中央子午线投影带下的高斯平面直角坐标。换带计算只涉及中央子午线经度和带宽，不涉及其它参考椭球。因此，不能简单的将坐标转换当做换带计算。

常用的坐标系及其 EPSG 编码

  1. 地理坐标系 EPSG编码 坐标系 说明 4326 WGS 1984 GPS采用的坐标系 4214 Beijing 1954 北京1954坐标系 4610 Xian 1980 西安1980坐标系 4490 CGCS2000/Gauss-Kruger 国家2000
  2. 投影坐标系 EPSG编码 坐标系 说明 3857 WGS 1984/Pseudo-Mercator(web墨卡托投影) WGS 1984 Web墨卡托投影坐标系(曾经代码EPSG:3785、EPSG:900913) 102025 Asia_North_Albers_Equal_Area_Conic 亚洲北部阿尔伯斯等积圆锥投影坐标系 4510 CGCS2000/Gauss-Kruger CM 123E(国家2000) CGCS2000地理坐标系下中央经线在123度E的高斯克吕格投影坐标系(between 120°E and 126°E) 4528 CGCS2000/Gauss-Kruger (between 118°30’E and 121°30’E) 32650 WGS 84/笛卡尔二维坐标系 Between 114°E and 120°E

七参数模型介绍 与 python 实现

七参数转换法

七参数法（包括布尔莎模型，一步法模型，海尔曼特等）是解决此问题的比较严密和通用的方法1。其涉及到的七个参数为：X平移，Y平移，Z平移，X旋转，Y旋转，Z旋转，尺度变化m。此七个参数可以通过在需要转化的区域里选取3个以上的转换控制点对而获取。

如果区域范围不大，最远点间的距离不大于30Km，可以用三参数（莫洛登斯基模型），即X平移，Y平移，Z平移，而此时将X旋转，Y旋转，Z旋转，尺度变化m视为0。所以三参数只是七参数的一种特例。三参数只需通过1个控制点对就能获取。

七参数公式如下:

img
\begin{bmatrix}
 X\\
 Y\\
Z
\end{bmatrix}
=
\begin{bmatrix}
 \Delta X\\
 \Delta Y\\
\Delta Z
\end{bmatrix}
+(1+m)\begin{bmatrix}
  1&  \theta_z& -\theta_y\\
  -\theta_z & 1 &\theta_x \\
  \theta_y&  -\theta_x& 1
\end{bmatrix}
\begin{bmatrix}
 X_0\\
 Y_0\\
Z_0
\end{bmatrix}

式中, \Delta X, \Delta Y, \Delta Z 为三个平移参数; \theta_x, \theta_y, \theta_z 为3个旋转参数，m 为尺度参数。将七参数作为未知数，并且按参数线性化，则可以转换为：

img

“`latex
\begin{bmatrix}
X-X(X_0)\
Y-Y(Y_0)\
Z-Z(Z_0)

\end{bmatrix}

\begin{bmatrix}
1 &0 &0 &0 &-(1-m)Z_0 &(1+m)Y_0 &X_0 + \theta_z Y_0-theta_y Z_0\
0 &1 &0 &(1+m)Z_0 &0 &-(1+m)X_0 &-\theta_z X_0 + Y_0+theta_x Z_0\
0 &0 &1 &-(1+m)Y_0 &(1+m)X_0 &0 &\theta_y X_0 – \theta_x Y_0 + Z_0
\end{bmatrix}
\begin{bmatrix}
\Delta X\
\Delta Y\
\Delta Z\
\theta_x\
\theta_y\
\theta_z\
m
\end{bmatrix}

其中，X(X_0), Y(Y_0), Z(Z_0) 由七参数公式给出，利用最小二乘法法即可对该线性化的方程进行迭代求解。

## 七参数转换法实例
下面以python代码为例<sup>2</sup> ，利用已有四个控制点坐标，根据最小二乘来求解七参数模型的参数，并利用一个检查点来验证其转换精度。
> 导入数据

使用4个同名点的原坐标(x, y, z)和目标坐标（x’, y’, z’)来求解七参数

python
import numpy as np

source_points = [
(3381400.980, 395422.030, 32.956),
(3381404.344, 395844.239, 32.207),
(3382149.810, 396003.592, 33.290),
(3382537.793, 395985.359, 33.412)
]
target_points = [
(3380968.194, 539468.888, 13.875),
(3380977.154, 539890.934, 13.179),
(3381724.612, 540040.47, 14.273),
(3381724.636, 540040.485, 14.282)
]

归一化处理便于模型的快速迭代

ave_src = np.mean(np.array(source_points), axis=0)
ave_tar = np.mean(np.array(target_points), axis=0)
source_points -= ave_src
target_points -= ave_tar

> 正向计算

根据坐标转换公式来列立系数方程定义函数返回系数矩阵。使用4个同名点的坐标转换公式作为误差方程进行最小二乘平差，根据坐标转换公式来列立系数方程 **W_x=b**

python

定义函数返回系数矩阵 B， l

定义函数： point2matrix，

通过给定的同名点坐标列立误差方程B系数阵的部分

x, y, z： 原坐标值

args: 七参数误差值[Delta_X, Delta_Y, Delta_Z, theta_x, theta_y, theta_z, m]

返回值： W系数阵的部分

def point2matrix(x, y, z, args):
array = [
[1, 0, 0, 0, -(1+args[6])z, (1+args[6])y, x + args[5]y – args[4]z],
[0, 1, 0, (1+args[6])z, 0, -(1+args[6])x, -args[5]x + y + args[3]z],
[0, 0, 1, -(1+args[6])y, (1+args[6])x, 0, args[4]x – args[3]y + z]
]
return np.array(array)

定义函数： points2W

通过同名点序列列立误差方程B系数阵的整体

x, y, z： 同名点序列

args: 七参数误差值[Delta_X, Delta_Y, Delta_Z, theta_x, theta_y, theta_z, m]

返回值： W系数阵

def points2W(points, args):
big_mat = None
for (x, y, z) in points:
mat = point2matrix(x, y, z, args)
if big_mat is None:
big_mat = mat
else:
big_mat = np.concatenate((big_mat, mat))

return big_mat

定义函数： points2b

通过同名点坐标转换关系列立误差方程B系数阵的整体

x, y, z： 同名点的原坐标和目标坐标对组成的序列

args: 七参数误差值[Delta_X, Delta_Y, Delta_Z, theta_x, theta_y, theta_z, m]

返回值： b系数阵

def points2b(source, target, args):
big_mat = [0] * len(source) * 3

for i, ((x1, y1, z1), (x2, y2, z2)) in enumerate(zip(source, target)):
    (x0, y0, z0) = ordinationConvert(x1, y1, z1, args)
    big_mat[3*i + 0] = x2 - x0
    big_mat[3*i + 1] = y2 - y0
    big_mat[3*i + 2] = z2 - z0

return np.array(big_mat).transpose()

def ordinationConvert(x1, y1, z1, args):
x2 = args[0] + (1+args[6])( x1 + args[5]y1 – args[4]z1) y2 = args[1] + (1+args[6])(-args[5]x1 + y1 + args[3]z1)
z2 = args[2] + (1+args[6])( args[4]x1 – args[3]*y1 + z1)
return (x2, y2, z2)

> 系数求解

组成法方程 (W'W) x = (W'b) 并利用最小二乘求解 x

python
Args = np.array([0, 0, 0, 0, 0, 0, 0], dtype=’float64′)
parameters = np.array([1, 1, 1, 1, 1, 1, 1])

当七参数的误差值之和大于1e-10时，迭代运算得到更精确的结果

while np.fabs(np.array(parameters)).sum() > 1e-10 :
W = points2W(source_points, Args)
b = points2b(source_points, target_points, Args)
qxx = np.linalg.inv(np.dot(W.transpose(), W))
parameters = np.dot(np.dot(qxx, W.transpose()), b)
Args += parameters

打印七参数

print(“Args:”)
print(np.round(Args, 3))

七参数输出结果为

[−0.−0.−0.−0.0.−0.072−0.215]

> 模型精度评定和结果验证

评定最小二乘模型的精度结果，并使用预留的已知坐标同名点来验证七参数模型的效果

python

检查点坐标

source_test_points = [
(3381402.058, 395657.940, 32.728)
]

target_test_points = [
(3380972.424, 539704.811, 13.665)
]

归一化处理

source_test_points = np.array(source_test_points – ave_src)

单位权标准差，即所得模型的坐标精度

sigma0 = np.sqrt((b*b).sum() / 2)

参数标准差，即所得模型的七参数精度

sigmax = sigma0 * np.sqrt(np.diag(qxx))
print(‘单位权中误差: %.3f’ % (sigma0))
print(‘参数中误差:’)
print(np.round(sigmax,3))
(x2, y2, z2) = ordinationConvert(source_test_points[0, 0], source_test_points[0, 1], source_test_points[0, 2], Args) + ave_tar
print(‘模型预测结果: ‘)
print(‘[(%.3f, %.3f, %.3f)]’%(x2, y2, z2))
print(‘真实结果: ‘)
print(target_test_points)

> 精度检查

目标点与预测值对比

python
单位权中误差: 162.711
参数中误差:
[ 81.356 81.356 81.356 0.65 0.311 0.191 0.15 ]
模型预测结果:
[(3380987.508, 539711.311, 13.654)]
真实结果:
[(3380972.424, 539704.811, 13.665)]
“`

1

七参数转换是针对将地理坐标系所对应的空间直角坐标转换为另一坐标系的空间直接坐标。空间直角坐标系的原点位于地球参考椭球的中心，Z轴与地球自转轴平行并指向参考椭球的北极，X轴指向参考椭球的本初子午线，Y轴与X轴和Z轴相互垂直最终构成一个右手系。大地坐标系是以大地基准为基础建立起来的，大地基准又以参考椭球为基础，由此大地坐标系又被称为椭球坐标系。

2

代码存在中文注释，在Python3中可正常运行，若使用Python2，则需要去掉中文注释来正常运行。

参考文献:

Post Views: 172

标签：

您可能还喜欢...

发表回复

目录

← TOC