1. 概述
要解決這個問題首先得理解地球橢球這個概念，這里直接用武漢大學《大地測量學基礎》（孔詳元、郭際明、劉宗全）的解釋吧：
大地經緯度坐標系是地理坐標系的一種，也就是我們常說的經緯度坐標+高度。經緯度坐標用的雖然多，但是很多人並沒有理解經緯度的幾何意義：緯度是一種線面角度，是坐標點P的法線與赤道面的夾角（注意這個法線不一定經過球心）；經度是面面角，是坐標點P所在的的子午面與本初子午面的夾角。這也是為什么經度范圍是-180 ~ +180，緯度范圍卻是-90 ~ +90：
地心地固坐標系就是我們常用的笛卡爾空間直角坐標系了。這個坐標系以橢球球心為原點，本初子午面與赤道交線為X軸，赤道面上與X軸正交方向為Y軸，橢球的旋轉軸(南北極直線)為Z軸。顯然，這是個右手坐標系：
顯然，兩者都是表達的都是空間中某點P，只不過一個是經緯度坐標（BLH），一個是笛卡爾坐標（XYZ）；兩者是可以相互轉換的。
2. 推導
2.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點在橢球面上的點為\(P_0\)，那么根據矢量相加的性質，有：
其中，\(P_0\)也就是式(5)，而n是\(P_0\)在橢球面的法線單位矢量。
矢量在任意位置的方向都是一樣的，那么我們可以假設存在一個單位球(球的半徑為單位1)，將法線單位矢量移動到球心位置，可得法線單位矢量為：
因此有：
其中：
2.2. XYZ->BLH
根據式（8），可知：
因此有：
不過緯度Ｂ就不是那么好算了，首先需要計算法線Pn在赤道兩側的長度。根據圖1，有：
與式（4-3）比較可得：
顯然，由於：
有：
接下來如下圖所示，對圖1做輔助線：
有：
因而可得：
這個式子兩邊都有待定量B，需要用迭代法進行求值。具體可參看代碼實現，初始的待定值可取\(tanB = \frac{z}{\sqrt{x^2+y^2}}\)。
大地緯度B已知，那么求高度H就非常簡單了，直接根據式（8）中的第三式逆推可得：
匯總三式，可得：
3. 實現
根據前面的推導過程，具體的C/C++代碼實現如下：
#include <iostream>
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;
void Blh2Xyz(double &x, double &y, double &z)
{
double L = x * d2r;
double B = y * d2r;
double H = z;
double N = a / sqrt(1 - e * e * sin(B) * sin(B));
x = (N + H) * cos(B) * cos(L);
y = (N + H) * cos(B) * sin(L);
z = (N * (1 - e * e) + H) * sin(B);
}
void Xyz2Blh(double &x, double &y, double &z)
{
double tmpX = x;
double temY = y ;
double temZ = z;
double curB = 0;
double N = 0;
double calB = atan2(temZ, sqrt(tmpX * tmpX + temY * temY));
int counter = 0;
while (abs(curB - calB) * r2d > epsilon && counter < 25)
{
curB = calB;
N = a / sqrt(1 - e * e * sin(curB) * sin(curB));
calB = atan2(temZ + N * e * e * sin(curB), sqrt(tmpX * tmpX + temY * temY));
counter++;
}
x = atan2(temY, tmpX) * r2d;
y = curB * r2d;
z = temZ / sin(curB) - N * (1 - e * e);
}
int main()
{
double x = 113.6;
double y = 38.8;
double z = 100;
printf("原大地經緯度坐標：%.10lf\t%.10lf\t%.10lf\n", x, y, z);
Blh2Xyz(x, y, z);
printf("地心地固直角坐標：%.10lf\t%.10lf\t%.10lf\n", x, y, z);
Xyz2Blh(x, y, z);
printf("轉回大地經緯度坐標：%.10lf\t%.10lf\t%.10lf\n", x, y, z);
}
其最關鍵的還是計算大地緯度B時的迭代過程，其余的計算都只是套公式。數值計算中的很多算法都是采用迭代趨近的方法來趨近一個最佳解。最后的運行結果如下：
