1、高斯投影坐標正反算一、相關概念大地坐標系由大地基準面和地圖投影確定，由地圖投影到特定橢圓柱面后在南北兩極 剪開展開而成，是對地球表面的逼近，各國或地區有各自的大地基準面，我國目前主要采 用的基準面為：1 .WGS8裾準面，為GPSS準面，17屆國際大地測量協會上推薦，橢圓柱長半軸a=6378137m 短半軸 b=6356752.3142451m2 .西安80坐標系，1975年國際大地測量協會上推薦，橢圓柱長半軸 a=6378140m短 半軸 b=6356755.2881575m;3 .北京54坐標系，參照前蘇聯克拉索夫斯基橢球體建立，橢圓柱長半軸a=6378245m, 短半軸 b=635686

2、3.018773m;通常所說的高斯投影有三種，即投影后：a）角度不變（正角投影），投影后經線和緯線仍然垂直；b）長度不變；c）面積不變；大地坐標一般采用高斯正角投影，即在地球球心放一點光源，地圖投影到過與中央經 線相切的橢圓柱面上而成；可分帶投影，按中央經線經度值分帶，有每6度一帶或每3度一帶兩種（起始帶中央經線經度為均為3度，即：6度帶1帶位置0-6度，3度帶1帶位置1.5-4.5 度），即所謂的高斯-克呂格投影。圖表11高斯投影和分帶地球某點經度（L）為過該點和地球自轉軸的半圓與子午線所在半圓夾角，東半球為東經，西半球為西經；地球某點緯度（B）為所在水平面法線與赤道圓面的線面角。正算是已知

3、大地坐標（L, B）,求解高斯平面坐標（X, Y），為確保Y值為正，Y增加 500公里；反算則是由高斯平面坐標（X, Y）求解大地坐標（L, B）。、計算模型:地球橢球面由橢圓繞地球自轉軸旋轉 180度而成圖表1橢圓橢圓長半軸a,橢圓短半軸b,橢圓方程:圖表2橢球面橢球面方程:2 y a2+2 X b2 +2 z a2=1*與網上充斥的將函數關系先展開為泰勒級數，再依據投影規則確定各參數不同，本文直接依據空間立體三角函數關系得出結果。*/（一）正算由圖表 1, tan（B）=由方程式（1）, -dya2 Xo=-2 ?dxb2 y0令Xo = b sin。， y0 = a - cosa-可得卜

4、 n _ xo _ b-sincc _ b tana _ b2 tan (B) 在圖表 2 中，y0 a-cosa-cos (L) a cos (L) a2 cos (1)，則由橢圓方程，令k ='y = r，cos。(r?cos 0)2(r?sin 0)22+22a2b可知:正算依據公式（4）、（5）、（6）、（7）得到結果，其中a:地球橢球長半軸；b:地球橢球短半軸；B：該點緯度；L:該點經度減去中央經線L0后的值;X:大地x坐標值;Y:大地y坐標值。(6)式積分按積分原理由計算機求積分(二)反算由式(4)可得，三、程序代碼函數:/* 高斯投影正算函數 *輸入：double a ,

5、double b , (m_B,m_L)為大地坐標，L0為帶號(6度帶)，(x,y)為高斯平面坐標,y加上了 50000師量返回：none*/void gaosiforward (double a, double b, double m_B, double m_L, double L0, double&x, double &y)double B=m_B*3.14159265/180; / 換算成弧度double L=(m_L-6.0*L0+3.0)*3.14159265/180;/ 換算成弧度double xita=atan(b*b*tan(B)/a/a/cos(L);doub

6、le dxita=0.000001;double xi=dxita;x=0.0;double c=a*a/b/b;while(xi<xita)x+=dxita/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi);xi+=dxita;x*=a;y=a*cos(xita)*tan(L)/sqrt(c*sin(xita)*sin(xita)+cos(xita)*cos(xita);y+=500000.0;/* 高斯反算函數*輸入 : double a ,double b ， (B,L) 為大地坐標, L0 為帶號(6度帶)，(x,y) 為高斯 平面坐標，y加上了 50

7、000師量* 返回： none*void gaosibackward ( double a, double b, double x, double y, double L0, double&B, double &L)double dxi=0.000001;double xi=dxi;double X=0.0;double c= a* a/ b/b;while(X< x/ a)X+=dxi/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi);xi+=dxi;double r= a/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi);double Y=y-50000