《(完整word版)高斯投影正反算公式-新.doc》由会员分享，可在线阅读，更多相关《(完整word版)高斯投影正反算公式-新.doc（7页珍藏版）》请在金锄头文库上搜索。
1、高斯投影坐标正反算一、相关概念大地坐标系由大地基准面和地图投影确定，由地图投影到特定椭圆柱面后在南北两极剪开展开而成，是对地球表面的逼近，各国或地区有各自的大地基准面，我国目前主要采用的基准面为：1.WGS84基准面，为GPS基准面，17届国际大地测量协会上推荐，椭圆柱长半轴a=6378137m，短半轴b=6356752.3142451m；2.西安80坐标系，1975年国际大地测量协会上推荐，椭圆柱长半轴a=6378140m，短半轴b=6356755.2881575m;3.北京54坐标系，参照前苏联克拉索夫斯基椭球体建立,椭圆柱长半轴a=6378245m, 短半轴b=6356863.01877
2、3m;通常所说的高斯投影有三种，即投影后：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 椭圆椭圆长半轴，椭圆短半轴b, 椭圆方程： (1)图表 2椭球面椭球面方程：/*与网上充斥的将函数关系先展开为泰勒级数，再依据投影规则确定各参数不同，本文直接依据空间立体三角函数关系得出结果。*/（一）正算由图表1， 由方程式(1),令，可得在图表2中，则由椭圆方程，令可知：正算依据公式(4)、(5)、(6)、(7)得到结果，其中a:地球椭球长半轴；b:地球椭球短半轴；B:该点纬度；L:该点经度减
4、去中央经线L0后的值；X:大地x坐标值；Y: 大地y坐标值。（6）式积分按积分原理由计算机求积分。（二）反算由式（4）可得，三、程序代码函数：/*高斯投影正算函数*输入 : double a ,double b，(m_B,m_L)为大地坐标, L0为带号（6度带），(x,y)为高斯平面坐标,y加上了500000常量返回：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
5、_L-6.0*L0+3.0)*3.14159265/180; /换算成弧度double xita=atan(b*b*tan(B)/a/a/cos(L);double dxita=0.000001;double xi=dxita;x=0.0;double c=a*a/b/b;while(xixita)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;/*高斯反算函
6、数*输入 : double a ,double b，(B,L)为大地坐标, L0为带号（6度带），(x,y)为高斯平面坐标,y加上了500000常量*返回：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-500000.0;L=atan(Y/r/cos(xi);L=L*180/3.14159265-3.0+6.0*L0;B=atan(m_a*m_a*tan(xi)*cos(L)/m_b/m_b);B=B*180/3.14159265;
