80、54、84座標系七引數轉換演算法及Java程式碼
一、為什麼要進行座標轉換
我們所在地球是一個不規則的橢球，地表凹凸不平，地底密度不均，因此很難用一個簡單模型來概括。國際上根據建模座標系的原點不同分為參心座標系和地心座標系，其中參心座標系是指參考橢球的幾何中心座標基準，地心，通常分為：參心空間直角座標系（以x，y，z為其座標元素）和參心大地座標系（以B，L，H為其座標元素），比如我國的北京54座標系和西安80座標系；座標系是用地球質心作為原點，通常分為地心空間直角座標系(以x，y，z為其座標元素)和地心大地座標系(以B，L，H為其座標元素)，比如GPS採用的WGS84座標系。
除了原點以外，任何一個橢球的模型還會有長軸、短軸和扁率三個引數，我整理了下
座標系 | 長軸（單位：米） | 短軸（單位：米） | 扁率 |
北京54座標系 | 6378245 | 6356863 | 1/298.3 |
西安80座標系 | 6378140 | 6356755 | 1/298.25722101 |
WGS84座標系 | 6378137 | 6356752.314 | 1/298.25722356 |
所以，我們可見這三個座標系是原點不一，橢球模型不一的座標系，就如同三個達芬奇的雞蛋，雖然都是橢球狀，但是形狀還是不一樣的。所以如果直接做兩個座標系之間的轉換是明顯不行的。
那麼有沒有辦法進行轉換呢？答案是肯定的。在數學上，任何兩個空間直角座標系都可以通過七個引數，即XYZ三個軸的平移，繞XYZ三個軸的旋轉角以及兩個座標系的大小比例因子的轉換後重疊在一起。因此，我們在做著三個座標系的轉換基本思路就是在空間直角座標系裡，通過七引數的偏移校正，實現座標系之間引數的轉換。
二、怎樣座標轉換
好了，那我們如何轉換呢。首先，我再區分一下地理上的空間直角座標系、大地座標系和平面座標系的區別。
空間直角座標系和數學的空間直角座標系基本相同，是以一個點為原點，三個互相垂直的直線作為三個軸延伸出來一個座標系，通常表示為（x,y,z）；
大地座標系是其地面上一點的大地經度L為大地起始子午面與該點所在的子午面所構成的二面角，由起始子午面起算，向東為正，稱為東經（0~180），向西為負，稱為西經（0~180）；大地緯度B是經過該點作橢球面的法線與赤道面的夾角，由赤道面起算，向北為正，稱為北緯（0~90），向南為負，稱為南緯（0~90）；大地高H是地面點沿橢球的法線到橢球面的距離
平面座標系主要是我國測繪領域裡使用的座標系，屬於參心大地座標系，有一個座標原點（54位原蘇聯的普爾科沃，80位我們陝西省涇陽縣永樂鎮），Z軸平行於地球質心指向地極原點方向，大地起始子午面平行於格林尼治平均天文臺子午面，X軸在大地起始子午面內與Z軸垂直指向經度0方向，Y軸與Z、X軸成右手座標系。因為大地高程以1956年青島驗潮站求出的黃海平均水面為基準，一般情況預設為0，所以通常表示為（X,Y）。
下面我以一個例子說明下如何進行WGS84座標系轉西安80座標系
WGS84座標系是大地座標系，首先要轉換到空間直角座標系，公式如圖所示：
其中a,b是WGS84橢球長短軸。程式碼為：
//第一步轉換，大地座標系換換成空間直角座標系
private void BLHtoXYZ(double B, double L, double H, double aAxis, double bAxis) {
double dblD2R = Math.PI / 180;
double e1 = Math.sqrt(Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / aAxis;
double N = aAxis / Math.sqrt(1.0 - Math.pow(e1, 2) * Math.pow(Math.sin(B * dblD2R), 2));
targetX = (N + H) * Math.cos(B * dblD2R) * Math.cos(L * dblD2R);
targetY = (N + H) * Math.cos(B * dblD2R) * Math.sin(L * dblD2R);
targetZ = (N * (1.0 - Math.pow(e1, 2)) + H) * Math.sin(B * dblD2R);
}第二步是進行空間直角座標系裡七引數轉換，公式如圖：
其中△X，△Y，△Z是座標平移量，R（ω）是旋轉矩陣，（1+m）是比例因子，程式碼如下：
//第二步轉換，空間直角座標系裡七引數轉換
Ex = transParaSeven.m_dbXRotate / 180 * Math.PI;
Ey = transParaSeven.m_dbYRotate / 180 * Math.PI;
Ez = transParaSeven.m_dbZRotate / 180 * Math.PI;
double newX = transParaSeven.m_dbXOffset + transParaSeven.m_dbScale * targetX + targetY * Ez - targetZ * Ey + targetX;
double newY = transParaSeven.m_dbYOffset + transParaSeven.m_dbScale * targetY - targetX * Ez + targetZ * Ex + targetY;
double newZ = transParaSeven.m_dbZOffset + transParaSeven.m_dbScale * targetZ + targetX * Ey - targetY * Ex + targetZ;第三步是空間直角座標系轉大地座標系，轉後的大地座標系為80的大地座標系，如圖：
其中a,b是西安80橢球長短軸，程式碼如下：
//第三步轉換，空間直接座標系轉換為大地座標系
private void XYZtoBLH(double X, double Y, double Z, double aAxis, double bAxis) {
double e1 = (Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / Math.pow(aAxis, 2);
double e2 = (Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / Math.pow(bAxis, 2);
double S = Math.sqrt(Math.pow(X, 2) + Math.pow(Y, 2));
double cosL = X / S;
double B = 0;
double L = 0;
L = Math.acos(cosL);
L = Math.abs(L);
double tanB = Z / S;
B = Math.atan(tanB);
double c = aAxis * aAxis / bAxis;
double preB0 = 0.0;
double ll = 0.0;
double N = 0.0;
//迭代計算緯度
do {
preB0 = B;
ll = Math.pow(Math.cos(B), 2) * e2;
N = c / Math.sqrt(1 + ll);
tanB = (Z + N * e1 * Math.sin(B)) / S;
B = Math.atan(tanB);
}
while (Math.abs(preB0 - B) >= 0.0000000001);
ll = Math.pow(Math.cos(B), 2) * e2;
N = c / Math.sqrt(1 + ll);
targetZ = Z / Math.sin(B) - N * (1 - e1);
targetB = B * 180 / Math.PI;
targetL = L * 180 / Math.PI;
}第四步是進行高斯變換，將大地座標轉為平面座標，公式如圖：
有點複雜，也一時說不清，具體可以查查高斯變換，這個還是很多的。程式碼如下：
//第四部轉換，高斯變換，大地座標系轉平面座標系，84轉80
private void gaussBLtoXY(double mX,double mY){
double m_aAxis = Axis_WGS84_a; //參考橢球長半軸
double m_bAxis = Axis_WGS84_b; //參考橢球短半軸
double m_dbMidLongitude = transParaSeven.daihao*3;//中央子午線經度 濟南117 威海123 巴州 87
double m_xOffset = 500000;
double m_yOffset = 0.0;
try{
//角度到弧度的係數
double dblD2R = Math.PI / 180;
//代表e的平方
double e1 = (Math.pow(m_aAxis, 2) - Math.pow(m_bAxis, 2)) / Math.pow(m_aAxis, 2);
//代表e'的平方
double e2 = (Math.pow(m_aAxis, 2) - Math.pow(m_bAxis, 2)) / Math.pow(m_bAxis, 2);
//a0
double a0 = m_aAxis * (1 - e1) * (1.0 + (3.0 / 4.0) * e1 + (45.0 / 64.0) * Math.pow(e1, 2) + (175.0 / 256.0) * Math.pow(e1, 3) + (11025.0 / 16384.0) * Math.pow(e1, 4));
//a2
double a2 = -0.5 * m_aAxis * (1 - e1) * (3.0 / 4 * e1 + 60.0 / 64 * Math.pow(e1, 2) + 525.0 / 512.0 * Math.pow(e1, 3) + 17640.0 / 16384.0 * Math.pow(e1, 4));
//a4
double a4 = 0.25 * m_aAxis * (1 - e1) * (15.0 / 64 * Math.pow(e1, 2) + 210.0 / 512.0 * Math.pow(e1, 3) + 8820.0 / 16384.0 * Math.pow(e1, 4));
//a6
double a6 = (-1.0 / 6.0) * m_aAxis * (1 - e1) * (35.0 / 512.0 * Math.pow(e1, 3) + 2520.0 / 16384.0 * Math.pow(e1, 4));
//a8
double a8 = 0.125 * m_aAxis * (1 - e1) * (315.0 / 16384.0 * Math.pow(e1, 4));
////緯度轉換為弧度表示
//B
double B = mX * dblD2R;
//l
double l = (mY - m_dbMidLongitude) * dblD2R;
////X
double X = a0 * B + a2 * Math.sin(2.0 * B) + a4 * Math.sin(4.0 * B) + a6 * Math.sin(6.0 * B) + a8 * Math.sin(8.0 * B);
//
double ll = Math.pow(Math.cos(B), 2) * e2;
double c = m_aAxis * m_aAxis / m_bAxis;
//N
double N = c / Math.sqrt(1 + ll);
//t
double t = Math.tan(B);
double p = Math.cos(B) * l;
double dby = X + N * t * (1 + ((5.0 - t * t + (9.0 + 4.0 * ll) * ll) + ((61.0 + (t * t - 58.0) * t * t + (9.0 - 11.0 * t * t) * 30.0 * ll) + (1385.0 + (-31111.0 + (543 - t * t) * t * t) * t * t) * p * p / 56.0) * p * p / 30.0) * p * p / 12.0) * p * p / 2.0;
double dbx;
dbx = N * (1.0 + ((1.0 - t * t + ll) + ((5.0 + t * t * (t * t - 18.0 - 58.0 * ll) + 14 * ll) + (61.0 + (-479.0 + (179.0 - t * t) * t * t) * t * t) * p * p / 42.0) * p * p / 20.0) * p * p / 6.0) * p;
mTargetX = dbx + m_xOffset;
mTargetY = dby + m_yOffset;
}
catch (Exception ex) {
}
}OK，也就搞定了。其他的轉換也就和這個差不多了，只是步驟相反而已