斜率公式k(完整word版)高斯投影正反算公式_新

高斯投影坐标正反算

一、相关概念
大地坐标系由大地基准面和地图投影确 定，由地图投影到特定椭
圆柱面后在南北两极剪开展开而成，是对地球表面的逼近，各国或地
区 有各自的大地基准面，我国目前主要采用的基准面为：
84基准面，为GPS基准面，17届国际大地 测量协会上
推荐，椭圆柱长半轴a=6378137m，短半轴
b=6356752.3142 451m；
2.西安80坐标系，1975年国际大地测量协会上推荐，椭圆柱长
半轴a=6 378140m，短半轴b=6356755.2881575m;
3.北京54坐标系，参照前苏联 克拉索夫斯基椭球体建立,椭圆柱
长半轴a=6378245m, 短半轴b=6356863.018773m;
通常所说的高斯投影有三种，即投影后：
a) 角度不变（正角投影），投影后经线和纬线仍然垂直；
b) 长度不变；
c) 面积不变；
大地坐标一般采用高斯正角投影，即在地球球心放一点光源，地
图 投影到过与中央经线相切的椭圆柱面上而成；可分带投影，按中央
经线经度值分带，有每6度一带或每3 度一带两种(起始带中央经线
经度为均为3度，即：6度带1带位置0-6度，3度带1带位置1.5- 4.5
度)，即所谓的高斯-克吕格投影。
图表 11高斯投影和分带


地球某点经度（L）为过该点和地球自转轴的半圆与子午线所在
半圆夹角，东半球为东经，西半 球为西经；地球某点纬度（B）为所
在水平面法线与赤道圆面的线面角。
正算是已知大地坐标 （L，B），求解高斯平面坐标（X，Y）,为
确保Y值为正，Y增加500公里；反算则是由高斯平面 坐标（X，Y）
求解大地坐标（L，B）。
二、计算模型：
地球椭球面由椭圆绕地球自转轴旋转180度而成。

图表 1 椭圆
椭圆长半轴
a
，椭圆短半轴b, 椭圆方程：
(1)

图表 2椭球面
椭球面方程：

***************************************
与网 上充斥的将函数关系先展开为泰勒级数，再依据投影规则确
定各参数不同，本文直接依据空间立体三角函 数关系得出结果。
*****
（一）正算
由图表1，
由方程式(1),

令，可得
在图表2中，
，则
由椭圆方程，令

可知：
正算依据公式(4)、(5)、(6)、(7)得到结果，其中
a:地球椭球长半轴；
b:地球椭球短半轴；
B:该点纬度；
L:该点经度减去中央经线L
0
后的值；
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.;换算成弧度
double L=(m_L-6.0*L0+3.0)*3.; 换算成弧
度
double xita=atan(b*b*tan(B)aacos(L));
double dxita=0.000001;
double xi=dxita;
x=0.0;
double c=a*abb;
while(xi {


}
x*=a;
y=a*cos(xita)*t an(L)sqrt(c*sin(xita)*sin(xita)+cos(xita)*c
os( xita));
y+=500000.0;
}

**************高斯反算函数***************
输入 : double a ,double b，(B,L)为大地坐标, L0为带号（6度
带），(x,y)为高斯平面坐标,y加上了500000常量
*返回：none
x+=dxitasqrt(c*sin(xi)*sin(xi)+co s(xi)*cos(xi));
xi+=dxita;
*****************************
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 bb;
while(X< x a)
{


}
double r= asqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi));
double Y=y-500000.0;

L=atan(Yrcos(xi));
L=L*1803.14159265-3.0+6.0*L0;

B=atan(m_a*m_a*tan(xi)*cos(L)m_bm_b);
B=B*1803.14159265;
}
X+=dxisqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi));
xi+=dxi;