高斯投影正反算公式新

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