高斯投影正反算公式

高斯投影坐标正反算一、基本思想：高斯投影正算公式就是由大地坐标（L，B）求解高斯平面坐标（x，y），而高斯投影反算公式则是由高斯平面坐标（x，y）求解大地坐标（L，B）。二、计算模型：基本椭球参数：椭球长半轴a椭球扁率f椭球短半轴：(1)baf椭球第一偏心率：22abea椭球第二偏心率：22abeb高斯投影正算公式：此公式换算的精度为0.001m6425644223422)5861(cossin720)495(cos24cossin2lttBBNltBsimBNlBBNXx5222425532233)5814185(cos120)1(cos6cosltttBNltBNlBNy其中：角度都为弧度B为点的纬度，0lLL，L为点的经度，0L为中央子午线经度；N为子午圈曲率半径，1222(1sin)NaeB；tantB；222coseB1803600其中X为子午线弧长：2402464661616sincos()(2)sinsin33XaBBBaaaaaBaB02468,,,,aaaaa为基本常量，按如下公式计算：200468242684468686883535281612815722321637816323216128mammmmmmammmammmmama02468,,,,mmmmm为基本常量，按如下公式计算：22222020426486379(1);;5;;268maememmemmemmem；高斯投影反算公式：此公式换算的精度为0.0001’’.2222243246532235242225053922461904572012cos6cos5282468120cosfffffffffffffffffffffffffffffttBByttyMNMNtyttyMNyyltNBNBytttNBLlL其中：0L为中央子午线经度。fB为底点纬度，也就是当xX时的子午线弧长所对应的纬度。按照子午线弧长公式：68240sin2sin4sin6sin82468aaaaXaBBBBB，迭代进行计算；初始开始时设：10fBXa以后每次迭代按下式计算：106824(())()sin2sin4sin6sin82468iiffiiiiifffffBXXFBaaaaaFBBBBB重复迭代至1iiffBB为止。1222(1sin)ffNaeB；32222(1)(1sin)ffMaeeBtanfftB；222cosffeB海福特椭球(1910)我国52年以前基准椭球a=6378388mb=6356911.9461279mα=0.33670033670克拉索夫斯基椭球(1940Krassovsky)北京54坐标系基准椭球a=6378245mb=6356863.018773mα=0.335232986921975年I.U.G.G推荐椭球(国际大地测量协会1975)西安80坐标系基准椭球a=6378140mb=6356755.2881575mα=0.0033528131778WGS-84椭球(GPS全球定位系统椭球、17届国际大地测量协会)WGS-84GPS基准椭球a=6378137mb=6356752.3142451mα=0.00335281006247三、程序代码函数：/************高斯投影正算函数***************输入:doublea,f椭球参数，B,L为大地坐标,L0为中央子午线的经度，单位为弧度，x,y为高斯平面坐标,y加上了500000常量返回：none******************************************/voidgaosiforward(doublea,doublef,doubleB,doubleL,doubleL0,double&x,double&y){doubleb,c,e1,e2;//短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率doublel,W,N,M,daihao;//W为常用辅助函数，N为子午圈曲率半径，M为卯酉圈曲率半径doubleX;//子午线弧长，高斯投影的坐标doubleruo,ita,sb,cb,t;doublem[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;/////bykjh2014.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]*pow(sin(B),4)/3.0);}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*cb*cb*cb*cb*cb*(61-58*t*t+t*t*t*t)*l*l*l*l*l*l/720);y=(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*ita-58*ita*ita*t*t)*l*l*l*l*l/120);y=y+500000;}/**************高斯反算函数***************输入:doublea,f椭球参数，x,y为高斯平面坐标，L0为中央子午线的经度；B,L为大地坐标，单位为弧度*返回：none*****************************/voidgaosibackward(doublea,doublef,doublex,doubley,doubleL0,double&B,double&L){doubleb,c,e1,e2;//短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率doubleBf,itaf,tf,Nf,Mf,Wf;doublel;doublem[5],n[5];y=y-500000;//计算一些基本常量{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;}//计算Bf{doubleBf1,Bfi0,Bfi1,FBfi;Bf1=x/n[0];Bfi0=Bf1;Bfi1=0;FBfi=0;intnum=0;do{num=0;FBfi=0.0-n[1]*sin(2*Bfi0)/2.0+n[2]*sin(4*Bfi0)/4.0-n[3]*sin(6*Bfi0)/6.0;Bfi1=(x-FBfi)/n[0];if(fabs(Bfi1-Bfi0)(Pi*pow(10.0,-8)/(36*18))){num=1;Bfi0=Bfi1;}}while(num==1);Bf=Bfi1;}tf=tan(Bf);Wf=sqrt(1-e1*e1*sin(Bf)*sin(Bf));Nf=a/Wf;Mf=a*(1-e1*e1)/(Wf*Wf*Wf);itaf=e2*cos(Bf);B=Bf-tf*y*y/(2*Mf*Nf)+tf*(5+3*tf*tf+itaf*itaf-9*itaf*itaf*tf*tf)*pow(y,4)/(24*Mf*pow(Nf,3))-tf*(61+90*tf*tf+45*pow(tf,4))*pow(y,6)/(720*Mf*pow(Nf,5));l=y/(Nf*cos(Bf))-(1+2*tf*tf+itaf*itaf)*pow(y,3)/(6*pow(Nf,3)*cos(Bf))+(5+28*tf*tf+24*pow(tf,4)+6*itaf*itaf+8*itaf*itaf*tf*tf)*pow(y,5)/(120*pow(Nf,5)*cos(Bf));L=l+L0;}2014-5-22'输入:doublea,f椭球参数，B,L为大地坐标,L0为中央子午线的经度，单位为弧度，x,y为高斯平面坐标,y加上了常量PrivateFunctiongaosiforward(ByValaAsDouble,ByValfAsDouble,ByValBAsDouble,ByValLAsDouble,ByValL0AsDouble)AsDouble()Dimx,y,xy(2)AsDoubleDimbb,c,e1,e2AsDouble'短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率Dimll,W,N,M,daihaoAsDouble'W为常用辅助函数，N为子午圈曲率半径，M为卯酉圈曲率半径DimxxAsDouble'子午线弧长，高斯投影的坐标Dimruo,ita,sb,cb,tAsDoubleDimmm(5),nn(5)AsDoublebb=a*(1-f)e1=Math.Sqrt(a*a-bb*bb)/ae2=Math.Sqrt(a*a-bb*bb)/bbc=a*a/bbmm(0)=a*(1-e1*e1)mm(1)=3*(e1*e1*mm(0))/2.0mm(2)=5*(e1*e1*mm(1))/4.0mm(3)=7*(e1*e1*mm(2))/6.0mm(4)=9*(e1*e1*mm(3))/8.0nn(0)=mm(0)+mm(1)/2+3*mm(2)/8+5*mm(3)/16+35*mm(4)/128nn(1)=mm(1)/2+mm(2)/2+15*mm(3)/32+7*mm(4)/16nn(2)=mm(2)/8+3*mm(3)/16+7*mm(4)/32nn(3)=mm(3)/32+mm(4)/16nn(4)=mm(4)/128xx=nn(0)*B-Sin(B)*Cos(B)*((nn(1)-nn(2)+nn(3))+(2*nn(2)-(16*nn(3)/3.0))*Sin(B)*Sin(B)+16*nn(3)*Pow(Sin(B),4)/3.0)ll=L-L0'弧度ita=e2*Cos(B)sb=Sin(B)cb=Cos(B)W=Sqrt(1-e1*e1*sb*s