技術背景
做過位置數據處理的小伙伴基本上都會遇到坐標轉換，而基於高斯投影原理的大地坐標轉平面坐標就是其中一種坐標轉換，坐標轉換的目的就是方便后面數據的處理工作，大地坐標轉高斯平面坐標常用的有兩種，即3°帶和6°帶，具體采用哪種根據實際情況而定。
計算原理
6°帶帶號n與相應的中央子午線L0經度的關系為：
3°帶帶號n’與相應的中央子午線L0’經度的關系為：
設參考橢球的長半軸為 a，第一偏心率為 e，並令：
設中央子午線的經度為 L0，再記：
則高斯投影正算公式為：
其中：
沒學過測繪的同學對以上原理不是很理解，也很正常，大家可以查閱相關《大地測量學》書籍，具體武大版還是礦大版的，沒什么區別。
具體實現
具體實現平台依然是IBM的Eclipse軟件，編程語言為Java，下面是以3°帶為例，進行高斯投影坐標正算，具體內容請看代碼
1 package package1; 2 3 public class BLH_xyh { 4 5 public static double a = 6378137; 6 public static double e = Math.sqrt(0.0066943799013); 7 public static double scale_wide = 3; 8 public static double π = 3.14159265358979323846; 9 10 public static void main(String[] args) { 11 // TODO 自動生成的方法存根 12 Point3d xyh = the_coordinates_are_counting(36.0307523111,120.184664478,9.7065); 13 System.out.println(xyh.getX()+","+xyh.getY()+","+xyh.getZ()); 14 } 15 public static Point3d the_coordinates_are_counting(double B,double L,double H){ 16 double A_ = 1 17 +3*e*e/4 18 +45*e*e*e*e/64 19 +175*e*e*e*e*e*e/256 20 +11025*e*e*e*e*e*e*e*e/16384 21 +43659*e*e*e*e*e*e*e*e*e*e/65536; 22 double B_ = 3*e*e/4+15*e*e*e*e/16 23 +525*e*e*e*e*e*e/512 24 +2205*e*e*e*e*e*e*e*e/2048 25 +72765*e*e*e*e*e*e*e*e*e*e/65536; 26 double C_ = 15*e*e*e*e/64 27 +105*e*e*e*e*e*e/256 28 +2205*e*e*e*e*e*e*e*e/4096 29 +10395*e*e*e*e*e*e*e*e*e*e/16384; 30 double D_ = 35*e*e*e*e*e*e/512 31 +315*e*e*e*e*e*e*e*e/2048 32 +31185*e*e*e*e*e*e*e*e*e*e/13072; 33 /* 34 double E_ = 315*e*e*e*e*e*e*e*e/16384 35 +3465*e*e*e*e*e*e*e*e*e*e/65536; 36 double F_ = 693*e*e*e*e*e*e*e*e*e*e/13072; 37 */ 38 39 double α = A_*a*(1-e*e); 40 double β = -B_*a*(1-e*e)/2; 41 double γ = C_*a*(1-e*e)/4; 42 double δ = -D_*a*(1-e*e)/6; 43 /* 44 double ε = E_*a*(1-e*e)/8; 45 double ζ = -F_*a*(1-e*e)/10; 46 */ 47 48 double C0 = α; 49 double C1 = 2*β+4*γ+6*δ; 50 double C2 = -8*γ-32*δ; 51 double C3 = 32*δ; 52 53 double x,y,sign; 54 double scale_number = Math.floor(L/scale_wide); 55 if(L > (scale_number * scale_wide + scale_wide/2)){ 56 scale_number =scale_number + 1; 57 sign = -1; 58 }else{ 59 sign = 1; 60 } 61 62 double L0 = scale_wide*scale_number; 63 double l = Math.abs(L-L0); 64 B = B*π/180; 65 l = l*π/180; 66 double t = Math.tan(B); 67 double m0 = Math.cos(B)*l; 68 double η = Math.sqrt(e*e*Math.pow(Math.cos(B),2)/(1-e*e)); 69 double N = a/Math.sqrt(1-e*e*Math.pow(Math.sin(B), 2)); 70 71 double X0 = C0*B+Math.cos(B)*(C1*Math.sin(B)+C2*Math.pow(Math.sin(B),3)+C3*Math.pow(Math.sin(B), 5)); 72 73 x = X0 74 +N*t*m0*m0/2 75 +N*t*m0*m0*m0*m0*(5-t*t+9*η*η+4*η*η*η*η)/24 76 +N*t*m0*m0*m0*m0*m0*m0*(61-58*t*t+t*t*t*t)/720; 77 y = N*m0+ 78 N*m0*m0*m0*(1-t*t+η*η)/6+ 79 N*m0*m0*m0*m0*m0*m0*(5-18*t*t+t*t*t*t+14*η*η-58*η*η*t*t)/120; 80 81 y = y*sign+500000; 82 83 double h = H; 84 85 Point3d xyh = new Point3d(x,y,h); 86 return xyh; 87 } 88 }
其中Point3d的定義如下：
1 package package1; 2 3 public class Point3d { 4 private double x; 5 private double y; 6 private double z; 7 public Point3d(double x,double y,double z){ 8 this.x=x; 9 this.y=y; 10 this.z=z; 11 } 12 public double getX() { 13 return x; 14 } 15 public void setX(double x) { 16 this.x = x; 17 } 18 public double getY() { 19 return y; 20 } 21 public void setY(double y) { 22 this.y = y; 23 } 24 public double getZ() { 25 return z; 26 } 27 public void setZ(double z) { 28 this.z = z; 29 } 30 }
測試數據為36.0307523111,120.184664478,9.7065，運行結果如下：
至此結束
致謝
感謝山東科技大學北斗星光創客興趣學習小組的王老師對於原理文檔的整理以及鄭** C++代碼的技術分享！
參考文檔
1、山東科技大學”北斗星光創客”興趣學習小組GNSS技術文檔