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
來源:微信公衆号:政通技術團隊
出處:https://mp.weixin.qq.com/s/lUkBxdIt_oBC2H7W6bmBRA