80、5四、84座標系七參數轉換算法

時間  2019-11-10
標籤 座標系 參數 轉換 算法 简体版
原文   原文鏈接

1、爲何要進行座標轉換

   咱們所在地球是一個不規則的橢球，地表凹凸不平，地底密度不均，所以很難用一個簡單模型來歸納。國際上根據建模座標系的原點不一樣分爲參心座標系和地心座標系，其中參心座標系是指參考橢球的幾何中心座標基準，地心，一般分爲：參心空間直角座標系（以x，y，z爲其座標元素）和參心大地座標系（以B，L，H爲其座標元素），好比我國的北京54座標系和西安80座標系；座標系是用地球質心做爲原點，一般分爲地心空間直角座標系(以x，y，z爲其座標元素)和地心大地座標系(以B，L，H爲其座標元素)，好比GPS採用的WGS84座標系。git

   除了原點之外，任何一個橢球的模型還會有長軸、短軸和扁率三個參數，我整理了下54,80,84三個座標的橢球模型以下表:spa

座標系code

長軸（單位：米）htm

短軸（單位：米）blog

扁率get

北京54座標系數學

6378245it

6356863io

1/298.3table

西安80座標系

6378140

6356755

1/298.25722101

WGS84座標系

6378137

6356752.314

1/298.25722356

 

        因此，咱們可見這三個座標系是原點不一，橢球模型不一的座標系，就如同三個達芬奇的雞蛋，雖然都是橢球狀，可是形狀仍是不同的。因此若是直接作兩個座標系之間的轉換是明顯不行的。

   那麼有沒有辦法進行轉換呢？答案是確定的。在數學上，任何兩個空間直角座標系均可以經過七個參數，即XYZ三個軸的平移，繞XYZ三個軸的旋轉角以及兩個座標系的大小比例因子的轉換後重疊在一塊兒。所以，咱們在作着三個座標系的轉換基本思路就是在空間直角座標系裏，經過七參數的偏移校訂，實現座標系之間參數的轉換。

 2、怎樣座標轉換

        好了，那咱們如何轉換呢。首先，我再區分一下地理上的空間直角座標系、大地座標系和平面座標系的區別。

   空間直角座標系和數學的空間直角座標系基本相同，是以一個點爲原點，三個互相垂直的直線做爲三個軸延伸出來一個座標系，一般表示爲（x,y,z）；

   大地座標系是其地面上一點的大地經度L爲大地起始子午面與該點所在的子午面所構成的二面角，由起始子午面起算，向東爲正，稱爲東經（0~180），向西爲負，稱爲西經（0~180）；大地緯度B是通過該點做橢球面的法線與赤道面的夾角，由赤道面起算，向北爲正，稱爲北緯（0~90），向南爲負，稱爲南緯（0~90）；大地高H是地面點沿橢球的法線到橢球面的距離，一般表示爲（B,L,H）；

   平面座標系主要是我國測繪領域裏使用的座標系，屬於參心大地座標系，有一個座標原點（54位原蘇聯的普爾科沃，80位咱們陝西省涇陽縣永樂鎮），Z軸平行於地球質心指向地極原點方向，大地起始子午面平行於格林尼治平均天文臺子午面，X軸在大地起始子午面內與Z軸垂直指向經度0方向，Y軸與Z、X軸成右手座標系。由於大地高程以1956年青島驗潮站求出的黃海平均水面爲基準，通常狀況默認爲0，因此一般表示爲（X,Y）。

 

下面我以一個例子說明下如何進行WGS84座標系轉西安80座標系

   WGS84座標系是大地座標系，首先要轉換到空間直角座標系，公式如圖所示：

  其中a,b是WGS84橢球長短軸。代碼爲：

 

  1. //第一步轉換，大地座標系換換成空間直角座標系

  2. private void BLHtoXYZ(double B, double L, double H, double aAxis, double bAxis) {

  3. double dblD2R = Math.PI / 180;

  4. double e1 = Math.sqrt(Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / aAxis;

  5.  

  6. double N = aAxis / Math.sqrt(1.0 - Math.pow(e1, 2) * Math.pow(Math.sin(B * dblD2R), 2));

  7. targetX = (N + H) * Math.cos(B * dblD2R) * Math.cos(L * dblD2R);

  8. targetY = (N + H) * Math.cos(B * dblD2R) * Math.sin(L * dblD2R);

  9. targetZ = (N * (1.0 - Math.pow(e1, 2)) + H) * Math.sin(B * dblD2R);

  10. }

      第二步是進行空間直角座標系裏七參數轉換，公式如圖：

      其中△X，△Y，△Z是座標平移量，R（ω）是旋轉矩陣，（1+m）是比例因子，代碼以下：

 

  1. //第二步轉換，空間直角座標系裏七參數轉換

  2. Ex = transParaSeven.m_dbXRotate / 180 * Math.PI;

  3. Ey = transParaSeven.m_dbYRotate / 180 * Math.PI;

  4. Ez = transParaSeven.m_dbZRotate / 180 * Math.PI;

  5.  

  6. double newX = transParaSeven.m_dbXOffset + transParaSeven.m_dbScale * targetX + targetY * Ez - targetZ * Ey + targetX;

  7. double newY = transParaSeven.m_dbYOffset + transParaSeven.m_dbScale * targetY - targetX * Ez + targetZ * Ex + targetY;

  8. double newZ = transParaSeven.m_dbZOffset + transParaSeven.m_dbScale * targetZ + targetX * Ey - targetY * Ex + targetZ;

 

      注意單位爲弧度。

 

      第三步是空間直角座標系轉大地座標系，轉後的大地座標系爲80的大地座標系，如圖：

      

     其中a,b是西安80橢球長短軸，代碼以下：

 

  1. //第三步轉換，空間直接座標系轉換爲大地座標系

  2. private void XYZtoBLH(double X, double Y, double Z, double aAxis, double bAxis) {

  3. double e1 = (Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / Math.pow(aAxis, 2);

  4. double e2 = (Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / Math.pow(bAxis, 2);

  5.  

  6. double S = Math.sqrt(Math.pow(X, 2) + Math.pow(Y, 2));

  7. double cosL = X / S;

  8. double B = 0;

  9. double L = 0;

  10.  

  11. L = Math.acos(cosL);

  12. L = Math.abs(L);

  13.  

  14. double tanB = Z / S;

  15. B = Math.atan(tanB);

  16. double c = aAxis * aAxis / bAxis;

  17. double preB0 = 0.0;

  18. double ll = 0.0;

  19. double N = 0.0;

  20. //迭代計算緯度

  21. do {

  22. preB0 = B;

  23. ll = Math.pow(Math.cos(B), 2) * e2;

  24. N = c / Math.sqrt(1 + ll);

  25.  

  26. tanB = (Z + N * e1 * Math.sin(B)) / S;

  27. B = Math.atan(tanB);

  28. }

  29. while (Math.abs(preB0 - B) >= 0.0000000001);

  30.  

  31. ll = Math.pow(Math.cos(B), 2) * e2;

  32. N = c / Math.sqrt(1 + ll);

  33.  

  34. targetZ = Z / Math.sin(B) - N * (1 - e1);

  35. targetB = B * 180 / Math.PI;

  36. targetL = L * 180 / Math.PI;

  37. }

    第四步是進行高斯變換，將大地座標轉爲平面座標，公式如圖：

    有點複雜，也一時說不清，具體能夠查查高斯變換，這個仍是不少的。代碼以下：

 

  1. //第四部轉換，高斯變換，大地座標系轉平面座標系，84轉80

  2. private void gaussBLtoXY(double mX,double mY){

  3. double m_aAxis = Axis_WGS84_a; //參考橢球長半軸

  4. double m_bAxis = Axis_WGS84_b; //參考橢球短半軸

  5. double m_dbMidLongitude = transParaSeven.daihao*3;//中央子午線經度 濟南117 威海123 巴州 87

  6. double m_xOffset = 500000;

  7. double m_yOffset = 0.0;

  8. try{

  9. //角度到弧度的係數

  10. double dblD2R = Math.PI / 180;

  11. //表明e的平方

  12. double e1 = (Math.pow(m_aAxis, 2) - Math.pow(m_bAxis, 2)) / Math.pow(m_aAxis, 2);

  13. //表明e'的平方

  14. double e2 = (Math.pow(m_aAxis, 2) - Math.pow(m_bAxis, 2)) / Math.pow(m_bAxis, 2);

  15. //a0

  16. 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));

  17. //a2

  18. 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));

  19. //a4

  20. 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));

  21. //a6

  22. 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));

  23. //a8

  24. double a8 = 0.125 * m_aAxis * (1 - e1) * (315.0 / 16384.0 * Math.pow(e1, 4));

  25. ////緯度轉換爲弧度表示

  26. //B

  27. double B = mX * dblD2R;

  28. //l

  29. double l = (mY - m_dbMidLongitude) * dblD2R;

  30. ////X

  31. 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);

  32. //

  33. double ll = Math.pow(Math.cos(B), 2) * e2;

  34. double c = m_aAxis * m_aAxis / m_bAxis;

  35. //N

  36. double N = c / Math.sqrt(1 + ll);

  37. //t

  38. double t = Math.tan(B);

  39. double p = Math.cos(B) * l;

  40. 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;

  41. double dbx;

  42. 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;

  43. mTargetX = dbx + m_xOffset;

  44. mTargetY = dby + m_yOffset;

  45. }

  46. catch (Exception ex) {

  47. }

  48. }

 

    OK，也就搞定了。其餘的轉換也就和這個差很少了，只是步驟相反而已

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程