1. 引言

政通晶石孪生平台实现了大至地球级别、小到城市部件的多层级数字仿真能力，该能力背后所依赖的基础理论知识之一便是《大地测量学》和《地理信息系统》中均有讲解的“空间坐标参考”相关知识。这里，我们就简要介绍一下晶石孪生平台其底层所涉及到的有关地理科学的知识。更深层次的说，这些技术的基础来源于计算机科学的《计算机图形学》或者数学的《立体几何》。

2. 空间坐标参考2.1. 大地坐标系首先，必须得理解地球椭球这个概念，这里作者直接用武汉大学《大地测量学基础》（孔详元、郭际明、刘宗全）的解释：

大地经纬度坐标系是地理坐标系的一种，也就是我们常说的经纬度坐标+高度。

经纬度坐标用的虽然多，但是很多人并没有理解经纬度的几何意义：纬度是一种线面角度，是坐标点P的法线与赤道面的夹角（注意这个法线不一定经过球心）；经度是面面角，是坐标点P所在的的子午面与本初子午面的夹角。这也是为什么经度范围是-180 ~ +180，纬度范围却是-90 ~ +90：

2.2. 地心地固坐标系

地心地固坐标系就是我们常用的笛卡尔空间直角坐标系了。

这个坐标系以椭球球心为原点，本初子午面与赤道交线为X轴，赤道面上与X轴正交方向为Y轴，椭球的旋转轴(南北极直线)为Z轴。显然，这是个右手坐标系：

大地坐标系和地心地固坐标系两者都是表达的都是空间中某点P，只不过一个是经纬度坐标（BLH），一个是笛卡尔坐标（XYZ）；

两个坐标系之间是可以相互转换的。

2.3. Web墨卡托投影坐标系

无论是大地坐标还是地心坐标，本质上来说都是一种曲面坐标，或者立体坐标。然而在很多情况下我们更愿意使用平面坐标，并且单位最好与常用的长度单位（米）一致。

所以就产生了从曲面到平面的转换，这个过程也叫做投影，转换的结果也就是投影平面坐标系。

也许很难理解投影具体的算法，但是必须了解投影的一点特性是：投影后发生了形变，只能有限的保证几点相同的特性。

以我们国内情况来说，最常用的三种投影平面坐标系是横轴墨卡托投影，高斯-克吕格投影和UTM投影。

本质上来说，高斯-克吕格投影和UTM投影其实都是横轴墨卡托投影，横轴墨卡托投影也是用的最为广泛的地图投影方式。但是在GIS，尤其是WebGIS领域中，横轴墨卡托投影的使用远没有Web墨卡托投影方式用的多。最重要的原因是Web墨卡托投影的转换算法比横轴墨卡托投影要简单很多，符合Web的轻量化的特点。

2.4. 站心坐标系

无论是大地经纬度坐标系，还是地心地固坐标系，甚至于Web墨卡托投影坐标系，它们都是一种基于全局的坐标系，他们的坐标值都是很大的值，出于数值精度的考虑来说，这样的值是不方便进行空间计算的。

所以很多时候需要一个以观察者为中心的局部坐标系，这种坐标系最初多用于GPS中，为测站为中心点，所以叫做站心坐标系，这个中心点也叫做站心点。

当以选取的站心点为坐标原点，三个坐标轴定义为X轴指东、Y轴指北，Z轴指天，就是ENU（东北天）站心坐标系。这样，从地心地固坐标系转换成的站心坐标系，就会成为一个符合常人对地理位置认知的局部坐标系。同时，只要站心点位置选的合理（通常可选取地理表达区域的中心点），其站心附近的空间坐标将会从一个很大的值变换成一个较小的值，非常适合于站心附近空间范围内的空间计算。

注意站心天向(法向量)与赤道面相交不一定会经过球心。

3. 空间坐标转换3.1. 大地经纬度坐标与地心地固坐标的的转换3.1.1. 大地坐标BLH->地心坐标XYZ

将P点所在的子午椭圆放在平面上，以圆心为坐标原点，建立平面直接坐标系：

对照地心地固坐标系，很容易得出：

那么，关键问题在于求子午面直角坐标系的x,y。过P点作原椭球的法线Pn，他与子午面直角坐标系X轴的夹角为B；过P点作子午椭圆的切线，它与X轴的夹角为(90°+B):

根据椭圆的方程，位于椭圆的P点满足：

对x求导，有：

又根据解析几何可知，函数曲线（椭圆）某一点（就是P点）的倒数为该点切线的斜率，也就是正切值：

联立公式（2）(3)，可得：

其中，e为椭圆第一偏心率：

令Pn的距离为N，那么显然有：

根据（4）式可得

将其带入（1）式，可得到椭球上P点的坐标为

那么唯一的未知量就是Pn的长度N了，将（4）式带入到椭圆方程式(1.2)：

化简，得：

联立式（5）式（6），得：

通过式（5）式（6），可以计算椭球上某一点的坐标。但这个点并不是我们真正要求的点，我们要求的点P(B,L,H)是椭球面沿法向量向上H高度的点：

P点在椭球面上的点为

,那么根据矢量相加的性质，有：

其中，

也就是式(5)，而n是

在椭球面的法线单位矢量。

矢量在任意位置的方向都是一样的，那么我们可以假设存在一个单位球(球的半径为单位1)，将法线单位矢量移动到球心位置，可得法线单位矢量为：

因此有：

其中：

3.1.2. 地心坐标XYZ->大地坐标BLH

根据式（8），可知：

因此有：

不过纬度Ｂ就不是那么好算了，首先需要计算法线Pn在赤道两侧的长度。根据图（1），有：

与式（4-3）比较可得：

显然，由于：

有：

接下来如下图所示，对图1做辅助线：

有：

因而可得：

这个式子两边都有待定量B，需要用迭代法进行求值。具体可参看代码实现，初始的待定值可取

。

大地纬度B已知，那么求高度H就非常简单了，直接根据式（8）中的第三式逆推可得：

汇总三式，可得：

3.1.3. 具体实现

根据前面的推导过程，可以用个人擅长的编程语言进行编码。

其最关键的还是计算大地纬度B时的迭代过程，其余的计算都只是套公式。数值计算中的很多算法都是采用迭代趋近的方法来趋近一个最佳解。

在作者的个人测试中，最后的运行结果如下：

3.2. 大地经纬度坐标系与Web墨卡托坐标系的转换

Web墨卡托投影是横轴墨卡托投影的特化版，要完全搞清楚Web墨卡托投影就必须得先搞清楚横轴墨卡托投影，不过横轴墨卡托投影实在太复杂了，但是我们可以定性地去理解。

它的计算过程大概可以这样理解：

在X方向上，为了保证投影到平面后经线和纬线仍然垂直，那么每条纬线都会按照赤道周长展开，也就是由于原点位于平面中心，那么可以算得X轴的取值范围：[-20037508.3427892,20037508.3427892]。经度与投影后X长是简单的线性关系；在Y方向上，则需要借助于墨卡托投影公式。为了保证投影的结果是正方形，那么就把Y轴的取值范围也取值成[-20037508.3427892,20037508.3427892]之间。这样做没什么道理，纯粹是为了希望投影的结果是正方形，便于切片。最后，通过墨卡托投影公式进行反算，得到的经纬度范围就是[-85.05112877980659，85.05112877980659]。也就是这种投影方式，大于这个范围是失效的。

参考Cesium的具体实现如下：

#include <iostream>//#include <eigen3/Eigen/Eigen>//#include <osgEarth/GeoData>using namespace std;const double epsilon = 0.000000000000001;const double pi = 3.14159265358979323846;const double d2r = pi / 180;const double r2d = 180 / pi;const double a = 6378137.0;  //椭球长半轴const double f_inverse = 298.257223563;   //扁率倒数const double b = a - a / f_inverse;//const double b = 6356752.314245;   //椭球短半轴const double e = sqrt(a * a - b * b) / a;//墨卡托范围[-PI, PI]->大地纬度范围[-PI/2, PI/2]static double mercatorAngleToGeodeticLatitude(double mercatorAngle){ return pi / 2.0 - (2.0 * atan(exp(-mercatorAngle))); //return 2.0 * atan(exp(mercatorAngle)) - pi / 2.0;}//Web墨卡托投影所支持的最大纬度（北和南）static double maximumLatitude = mercatorAngleToGeodeticLatitude(pi);//大地纬度范围[-PI/2, PI/2]->墨卡托范围[-PI, PI]static double geodeticLatitudeToMercatorAngle(double latitude){ // Clamp the latitude coordinate to the valid Mercator bounds. if (latitude > maximumLatitude) {  latitude = maximumLatitude; } else if (latitude < -maximumLatitude) {  latitude = -maximumLatitude; } double sinLatitude = sin(latitude); return 0.5 * log((1.0 + sinLatitude) / (1.0 - sinLatitude));} void Blh2Wmc(double &x, double &y, double &z){ x = x * d2r * a; y = geodeticLatitudeToMercatorAngle(y * d2r) * a;}void Wmc2Blh(double &x, double &y, double &z){ //var oneOverEarthSemimajorAxis = this._oneOverSemimajorAxis; x = x / a * r2d; y = mercatorAngleToGeodeticLatitude(y / a) * r2d; }int main(){ double x = 113.6; double y = 38.8; double z = 100;     printf("%.10lf\n", maximumLatitude * r2d); printf("原大地经纬度坐标：%.10lf\t%.10lf\t%.10lf\n", x, y, z); Blh2Wmc(x, y, z); printf("Web墨卡托坐标：%.10lf\t%.10lf\t%.10lf\n", x, y, z); Wmc2Blh(x, y, z); printf("转回大地经纬度坐标：%.10lf\t%.10lf\t%.10lf\n", x, y, z);}

最终运行的结果：

3.3. 地心地固坐标系(ECEF)与站心坐标系(ENU)的转换3.3.1. 原理

令选取的站心点为P，其大地经纬度坐标为

，对应的地心地固坐标系为

。地心地固坐标系简称为ECEF，站心坐标系简称为ENU。

3.3.1.1. 平移

通过2.4节的图可以看出，ENU要转换到ECEF，一个很明显的图形操作是平移变换，将站心移动到地心。根据站心点P在地心坐标系下的坐标

，根据《计算机图形学》的知识，可以很容易推出ENU转到ECEF的平移矩阵：

反推之，ECEF转换到ENU的平移矩阵就是T的逆矩阵：

3.3.1.2. 旋转

另外一个需要进行的图形变换是旋转变换，其旋转变换矩阵根据P点所在的经度L和纬度B确定。这个旋转变换有点难以理解，需要一定的空间想象能力，但是可以直接给出如下结论:

当从ENU转换到ECEF时，需要先旋转再平移，旋转是先绕X轴旋转再绕Z轴旋转；当从ECEF转换到ENU时，需要先平移再旋转，旋转是先绕Z轴旋转再绕X轴旋转。

根据《计算机图形学》旋转变换的相关知识，绕X轴旋转矩阵为：

绕Z轴旋转矩阵为：

从ENU转换到ECEF的旋转矩阵为：

根据三角函数公式：

有：

将(14)、(15)带入(13)中，则有：

而从ECEF转换到ENU的旋转矩阵为：

旋转矩阵是正交矩阵，根据正交矩阵的性质:正交矩阵的逆矩阵等于其转置矩阵，那么可直接得：

3.3.1.3. 总结

将上述公式展开，可得从ENU转换到ECEF的图形变换矩阵为：

而从ECEF转换到ENU的图形变换矩阵为：

3.3.2. 实现

接下来用代码实现这个坐标转换，选取一个站心点，以这个站心点为原点，获取某个点在这个站心坐标系下的坐标。作者个人测试，运行结果如下：

4. 扩展

以上基础知识不仅是晶石平台在地理信息空间方面的坐标转换所要用到，也是一个3DGIS开发人员必须了解的，对该知识的理解更有助于我们在三维图形中的具体编码实现。

但是，晶石平台的三维图形渲染有自己的图形空间坐标转换算法，它们与以上讲解的地理信息空间坐标转换有类似之处，也有不少区别的地方。这些区别的核心作用是让真实空间地物在晶石平台中得到更好的渲染仿真。

作者:charlee

来源:微信公众号:政通技术团队

出处: