发布时间 : 星期日 文章高斯投影正反算公式更新完毕开始阅读5776611cd4d8d15abf234e14
高斯投影坐标正反算
一、基本思想：
高斯投影正算公式就是由大地坐标（L，B）求解高斯平面坐标（x，y），而高斯投影反算公式则是由高斯平面坐标（x，y）求解大地坐标（L，B）。
二、计算模型：
基本椭球参数： 椭球长半轴a 椭球扁率f
椭球短半轴：b?a(1?f)
a2?b2椭球第一偏心率 ：e?
aa2?b2椭球第二偏心率 ：e?? b高斯投影正算公式：此公式换算的精度为0.001m
x?X??NN2??sinBcosB?l?simBcos3B(5?t2?9?2?4?4)l??4242???24???
N5246??sinBcosB(61?58t?t)l720???6y?NN3223??cosB?l???cosB(1?t??)l???6???3N5242225???cosB(5?18t?t?14??58?t)l120???5
其中：角度都为弧度
B为点的纬度，l???L?L0，L为点的经度，L0为中央子午线经度；
N为子午圈曲率半径，N?a(1?esinB)； t?tanB；
22?12?2?e?2cos2B
????180??3600
其中X为子午线弧长：
1616??X?a0B?sinBcosB?(a2?a4?a6)?(2a4?a6)sin2B?a6sin4B?
33??a0,a2,a4,a6,a8为基本常量，按如下公式计算： m23535?a?m??m?m?m8046?02816128??a?m2?m4?15m?7m68?2223216?m437? ?m6?m8?a4?81632?m6m8?a??632?16??a?m88?128?m0,m2,m4,m6,m8为基本常量，按如下公式计算：
m0?a(1?e2);m2?3279em0;m4?5e2m2;m6?e2m4;m8?e2m6； 268
高斯投影反算公式：此公式换算的精度为0.0001’’.
B?Bf??tftf2MfNfy?2tf24MfN3f?5?3t2f224??2f?9?ftf?y720MfN5f46y?61?90t2f?45tf?yyy322l??1?2t???ff?NfcosBf6N3cosBffy5422?5?28t2?2?f?24tf?6?f?8ftf?5120NfcosBfL?l?L0
其中：
L0为中央子午线经度。
Bf为底点纬度，也就是当x?X时的子午线弧长所对应的纬度。按照子午线弧长公式：
X?a0B?aaa2asin2B?4sin4B?6sin6B?8sin8B，迭代进行计算； 24681初始开始时设：Bf?Xa0 以后每次迭代按下式计算：
Bif?1?X(X?F(Bif))a0 a6a8a2a4iiiiF(B)??sin2Bf?sin4Bf?sin6Bf?sin8Bf2468if重复迭代至Bif?1?Bif??为止。
Nf?a(1?esinBf)； Mf?a(1?e)(1?esinBf)
222?3222?12tf?tanBf；
22??2?ecosBf f
海福特椭球(1910) 我国52年以前基准椭球 a=6378388m b=6356911.9461279m α=0.33670033670 球 球
克拉索夫斯基椭球(1940 Krassovsky) 北京54坐标系基准椭球 a=6378245m b=6356863.018773m α=0.33523298692
1975年I.U.G.G推荐椭球(国际大地测量协会1975) 西安80坐标系基准椭a=6378140m b=6356755.2881575m α=0.0033528131778
WGS-84椭球(GPS全球定位系统椭球、17届国际大地测量协会) WGS-84 GPS 基准椭a=6378137m b=6356752.3142451m α=0.00335281006247
三、程序代码函数：
/************高斯投影正算函数***************
输入 : double a ,f 椭球参数，B,L为大地坐标,L0为中央子午线的经度，单位为弧度，x,y为高斯平面坐标,y加上了500000常量 返回：none
******************************************/
void gaosiforward(double a,double f,double B,double L,double L0,double &x,double &y) {
double b, c,e1, e2; //短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率 double l, W,N, M, daihao;//W为常用辅助函数，N为子午圈曲率半径，M为卯酉圈曲率半径 double X;//子午线弧长，高斯投影的坐标 double ruo, ita, sb, cb,t; double m[5],n[5]; //计算一些基本常量 {
b=a*(1-f);
e1=sqrt(a*a-b*b)/a; e2=sqrt(a*a-b*b)/b; c=a*a/b;
m[0]=a*(1-e1*e1); m[1]=3*(e1*e1*m[0])/2.0;
}
}
m[2]=5*(e1*e1*m[1])/4.0; m[3]=7*(e1*e1*m[2])/6.0; m[4]=9*(e1*e1*m[3])/8.0;
n[0]=m[0]+m[1]/2+3*m[2]/8+5*m[3]/16+35*m[4]/128; n[1]=m[1]/2+m[2]/2+15*m[3]/32+7*m[4]/16; n[2]=m[2]/8+3*m[3]/16+7*m[4]/32; n[3]=m[3]/32+m[4]/16;
n[4]=m[4]/128; /////by kjh 2014.5.22 把改成了
//由纬度计算子午线弧长 {
X=n[0]*B-sin(B)*cos(B)*((n[1]-n[2]+n[3])+(2*n[2]-(16*n[3]/3.0))*sin(B)*sin(B)+16*n[3]*p}
l=L-L0;//弧度 ita=e2*cos(B); sb=sin(B); cb=cos(B);
W=sqrt(1-e1*e1*sb*sb); N=a/W; t=tan(B);
ruo=(180/Pi)*3600;
x=(X+N*sb*cb*l*l/2+N*sb*cb*cb*cb*(5-t*t+9*ita*ita+4*ita*ita*ita*ita)*l*l*l*l/24+N*sb*cby=(N*cb*l+N*cb*cb*cb*(1-t*t+ita*ita)*l*l*l/6+N*cb*cb*cb*cb*cb*(5-18*t*t+t*t*t*t+14*ita*y=y+500000;
ow(sin(B),4)/3.0);
*cb*cb*cb*cb*(61-58*t*t+t*t*t*t)*l*l*l*l*l*l/720); ita-58*ita*ita*t*t)*l*l*l*l*l/120);
/**************高斯反算函数***************
输入 : double a ,f 椭球参数， x,y为高斯平面坐标，L0为中央子午线的经度； B,L为大地坐标，单位为弧度 *返回：none
*****************************/
void gaosibackward(double a,double f,double x,double y,double L0,double &B,double &L) {
double b, c,e1, e2; //短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率 double Bf,itaf,tf,Nf,Mf,Wf; double l; double m[5],n[5]; y=y-500000; //计算一些基本常量 {
b=a*(1-f);