空间直角坐标系与空间大地坐标系地相互转换及其C++源程序

1、空间直角坐标系与空间大地坐标系的相互转换1. 空间直角坐标系/笛卡尔坐标系坐标轴相互正交的坐标系被称作笛卡尔坐标系。三维笛卡尔坐标系也被称为 空间直角坐标系。在空间直角坐标系下，点的坐标可以用该点所对应的矢径在三 个坐标轴上的投影长度来表示，只有确定了原地、三个坐标轴的指向和尺度，就 定义了一个在三维空间描述点的位置的空间直角坐标系。以椭球体中心0为原点，起始子午面与赤道面交线为 X轴，在赤道面上与 X轴正交的方向为丫轴，椭球体的旋转轴为Z轴构成右手坐标系O.XYZ，在该 坐标系中，P点的位置用X,Y,Z表示。在测量应用中，常将地球空间直角坐标系的坐标原点选在地球质心 （地心坐 标系）或参考椭

2、球中心（参心坐标系），z轴指向地球北极，x轴指向起始子午面 与地球赤道的交点，y轴垂直于XOZ面并构成右手坐标系。空间直角坐标系2. 空间大地坐标系由于空间直角坐标无法明确反映出点与地球之间的空间关系， 为了解决这一 问题，在测量中引入了大地基准，并据此定义了大地坐标系。大地基准指的是用 于定义地球参考椭球的一系列参数，包括如下常量：2.1椭球的大小和形状2.2椭球的短半轴的指向：通常与地球的平自转轴平息。2.3椭球中心的位置：根据需要确定。若为地心椭球，贝U其中心位于地球质 心。2.4本初子午线：通过固定平极和经度原点的天文子午线，通常为格林尼治 子午线。以大地基准为基础建立的坐标系被称为大

3、地坐标系。由于大地基准又以参考 椭球为基准，因此，大地坐标系又被称为椭球坐标系。大地坐标系是参心坐标系， 其坐标原点位于参考椭球中心，以参考椭球面为基准面，用大地经度L、纬度B和大地高H表示地面点位置。过地面点 P的子午面与起始子午面间的夹角叫 P 点的大地经度。由起始子午面起算，向东为正，叫东经(0 180 )向西为负， 叫西经(0 -180。过P点的椭球法线与赤道面的夹角叫 P点的大地纬度。由 赤道面起算，向北为正，叫北纬(0 90 )向南为负，叫南纬(0 -90 。从 地面点P沿椭球法线到椭球面的距离叫大地高。大地坐标坐标系中，P点的位置用L,B表示。如果点不在椭球面上，表示点的位置除

4、L,B外，还要附加另一参数 大地咼 H。空间大地坐标系3. 空间直角坐标与大地坐标间的转换3.1大地坐标转换为空间直角坐标换公式为：X =疝 COMy=(+jEf)cosBsin2o将同一坐标系下的大地坐标（B、L、H）转换成空间直角坐标（X、丫、Z）的转b为参考椭球的短半轴;式中N为卯酉圈的半径，a为参考椭球的长半轴; e为参考椭球的第一偏心率；并且有若点在椭球面上，则大地高 H=0 ,上式可简化为:X 二y = cosfisinl3.2空间直角坐标转换为空间大地坐标将同一坐标系下的空间直角坐标（X、丫、Z）转换为空间大地坐标（B、L、H）的公式为:YZ(N + H)Zarctao:firc

5、tanf fJ(xW)n(i7)+h-N(l-X) sia J3在使用上式进行空间直角坐标到大地坐标的转换过程中，由于计算大地纬度口时用到大地高W,而计算大地高时又需要用到大地纬度口. 因此不能直接由空 间直角坐标计算出大地坐标，而需要采用迭代计算的方法。具体计算时，可先根 据下式求出大地纬度口的初值：ZTx1+r2然后利用该初值来求出H、N的初值，再利用所求出的H和N初值再次求出B 值.如此反复，直至求出的及日、W收敛为止。4算例本文根据以上公式在 Microsoft VC+6.0环境下编写了一段程序(见附 录)o算例中的坐标采用的是武汉大学信息学部友谊广场上的某点的大地坐标作 为已知值，然

6、后经过转换函数CRDGEODETICtoCRDCARTESEAN(pcg, pcc,dSemiMajorAxis,dFlatning)把大地坐标转换为空间直角坐标得到坐标 X、丫、乙由得到的空间直角坐标X、Y、Z，经过转换函数CRDCARTESIANtoCRDGEODETIC(pcc,pcg,dSemiMajorAxis,dFlat ning)把空间直角坐标还原成空间大地坐标，计算结果如下图所示：计算结果从上图可以看出结果比较满意，高程和精度基本能完全还原，而纬度还原后 有较大的误差，在测量中这种误差不允许的，需要修改算法，完善结果。可能引 起的原因有可能是由于纬度计算公式并不完善， 还有可能

7、是由于计算机的截断误 差引起的，还要找时间继续修改、完善。5.心得体会这次编程自认为很简单，但真动手自己亲自编写，还是或多或少遇到了一些 问题，并分析问题，最终解决问题。虽然这次作业很简单，但经过自己这样一步 一步的编写出来，还是有很多收获，加强了运用VC+编写程序的能力，也充分 认识到了学习VC+的重要性，更找到了自己的一些缺点与不足。6附录(程序源代码)#in clude#in cludeusing n amespace std;#define M_PI 3.1415926typedef struct tagCRDCARTESIANdouble x;double y;double z;CR

8、DCARTESIAN;typedef CRDCARTESIAN *PCRDCARTESIAN;typedef struct tagCRDGEODETICdouble Ion gitude;double latitude;double height;CRDGEODETIC;typedef CRDGEODETIC *PCRDGEODETIC;void DMS_RAD(double DMS,double *Rad)int Deg,Mi n;double Sec;Deg=(i nt)DMS;Min=(i nt)(DMS-Deg)*100);Sec=(DMS-Deg)*100-Mi n)*100;*Ra

9、d=(Deg+Mi n/60.0+Sec/3600.0)/180.0*M_PI;return;void RAD_DMS(double Rad,double *DMS)int Deg,Mi n;double Sec;double AR,AM;AR=Rad;if (Rad0)AR=-Rad;AR=AR+1.0e-10;AR=AR*180.0/M_PI;Deg=(i nt)AR;AM=(AR-Deg)*60.0;Min=(i nt)AM;Sec=(AM-Mi n)*60;*DMS=Deg+Mi n/100.0+Sec/10000.0;if(Radx;double Y=pcc-y;double Z=p

10、cc-z;R=sqrt(X*X+Y*Y);B0=ata n2(Z,R);while (1)N=dSemiMajorAxis/sqrt(1.0-dFlatte nin g*(2-dFlatte nin g)*si n(B0)*si n(BO)；B_=ata n2(Z+N*dFlatte nin g*(2-dFlatte nin g)*si n(B0),R); if(fabs(B_-B0)height=R/cos(B_)-N;RAD_DMS(B _，&pcg-latitude);RAD_DMS(L_,&pcg-lo ngitude);return true;boolCRDGEODETICtoCRD

11、CARTESEAN(PCRDGEODETICpcg,PCRDCARTESIAN pcc,double dSemiMajorAxis,double dFlatte ning)double N;double B_,L_;double B=pcg-latitude;double L=pcg- Ion gitude;double H=pcg-height;DMS_RAD(B,&B_);DMS_RAD(L，&L_);N=dSemiMajorAxis/sqrt(1.0-dFlatte nin g*(2-dFlatte nin g)*si n(B_)* si n( B_);pcc-x=(N+H)*cos(B

12、_)*cos(L_);pcc-y=(N+H)*cos(B_)*si n(L_);pcc-z=(N*(1.0-dFlatte nin g*(2-dFlatte nin g)+H)*si n(B_); return true;void mai n()PCRDCARTESIAN pcc=new CRDCARTESIAN;PCRDGEODETIC pcg=new CRDGEODETIC;double rad;rad=(30*3600+31*60+40.23)/3600; pcg-latitude=rad;rad=(114*3600+21*60+20.51)/3600;pcg-height=41;pcg- Ion gitude=rad;double dSemiMajorAxis=6378137;double dFlat nin g=1/298.257223563;coutstd:fixed;coutvv转换前已知的大地坐标：endl;coutH( 高度)=heightL(度)=longitudevvB(纬度)=latitudeendl;CRDGEODETICtoCRDCARTESEAN（pcg,dSemiMajorAxis,dFlat nin g）;