空間直角坐標(biāo)系與空間大地坐標(biāo)系的相互轉(zhuǎn)換及其C++源程序

1、空間直角坐標(biāo)系與空間大地坐標(biāo)系的相互轉(zhuǎn)換1. 空間直角坐標(biāo)系/笛卡爾坐標(biāo)系坐標(biāo)軸相互正交的坐標(biāo)系被稱(chēng)作笛卡爾坐標(biāo)系。三維笛卡爾坐標(biāo)系也被稱(chēng)為空間直角坐標(biāo)系。在空間直角坐標(biāo)系下，點(diǎn)的坐標(biāo)可以用該點(diǎn)所對(duì)應(yīng)的矢徑在三個(gè)坐標(biāo)軸上的投影長(zhǎng)度來(lái)表示，只有確定了原地、三個(gè)坐標(biāo)軸的指向和尺度，就定義了一個(gè)在三維空間描述點(diǎn)的位置的空間直角坐標(biāo)系。以橢球體中心0為原點(diǎn)，起始子午面與赤道面交線(xiàn)為 X軸，在赤道面上與X軸正交的方向?yàn)檠?軸，橢球體的旋轉(zhuǎn)軸為Z軸構(gòu)成右手坐標(biāo)系O.XYZ在該坐標(biāo)系中，P點(diǎn)的位置用X,Y,Z表示。在測(cè)量應(yīng)用中，常將地球空間直角坐標(biāo)系的坐標(biāo)原點(diǎn)選在地球質(zhì)心（地心坐標(biāo)系）或參考橢球中心（參心坐

2、標(biāo)系），z軸指向地球北極，x軸指向起始子午面與地球赤道的交點(diǎn)，y軸垂直于XOZ 面并構(gòu)成右手坐標(biāo)系?？臻g直角坐標(biāo)系2. 空間大地坐標(biāo)系由于空間直角坐標(biāo)無(wú)法明確反映出點(diǎn)與地球之間的空間關(guān)系，為了解決這一問(wèn)題，在測(cè)量中引入了大地基準(zhǔn)，并據(jù)此定義了大地坐標(biāo)系。大地基準(zhǔn)指的是用于定義地球參考橢球的一系列參數(shù)，'I X 1 I y '包括如下常量：2.1橢球的大小和形狀 :2.2橢球的短半軸的指向：通常與地球的平自轉(zhuǎn)軸平息。2.3橢球中心的位置：根據(jù)需要確定。若為地心橢球，則其中心位于地球質(zhì)心。2.4本初子午線(xiàn)：通過(guò)固定平極和經(jīng)度原點(diǎn)的天文子午線(xiàn)，通常為格林尼治子午線(xiàn)。以大地基準(zhǔn)為基礎(chǔ)建

3、立的坐標(biāo)系被稱(chēng)為大地坐標(biāo)系。 由于大地基準(zhǔn)又以參考橢球?yàn)榛鶞?zhǔn)，因此， 大地坐標(biāo)系又被稱(chēng)為橢球坐標(biāo)系。 大地坐標(biāo)系是參心坐標(biāo)系，其坐標(biāo)原點(diǎn)位于參考橢球中心，以參 考橢球面為基準(zhǔn)面，用大地經(jīng)度 L、緯度B和大地高H表示地面點(diǎn)位置。過(guò)地面點(diǎn) P的子午面與起 始子午面間的夾角叫P點(diǎn)的大地經(jīng)度。由起始子午面起算，向東為正，叫東經(jīng)（0° 180°）,向西 為負(fù)，叫西經(jīng)（0° -180 °）。過(guò)P點(diǎn)的橢球法線(xiàn)與赤道面的夾角叫 P點(diǎn)的大地緯度。由赤道面起 算，向北為正，叫北緯（0° 90°）,向南為負(fù)，叫南緯（0° -90 °）。

4、從地面點(diǎn)P沿橢球法線(xiàn)到 橢球面的距離叫大地高。大地坐標(biāo)坐標(biāo)系中，P點(diǎn)的位置用L,B表示。如果點(diǎn)不在橢球面上，表示點(diǎn)的位置除L,B夕卜，還要附加另一參數(shù)一一大地高 Ho空間大地坐標(biāo)系3. 空間直角坐標(biāo)與大地坐標(biāo)間的轉(zhuǎn)換的轉(zhuǎn)換公式為：式中N為卯酉圈的半徑，a為參考橢球的長(zhǎng)半軸；b為參考橢球的短半軸；e為參考橢球的第 一偏心率；并且有若點(diǎn)在橢球面上，則大地高 H=0,上式可簡(jiǎn)化為：3.2空間直角坐標(biāo)轉(zhuǎn)換為空間大地坐標(biāo)將同一坐標(biāo)系下的空間直角坐標(biāo)（X、Y、Z）轉(zhuǎn)換為空間大地坐標(biāo)（B、L、H）的公式為：在使用上式進(jìn)行空間直角坐標(biāo)到大地坐標(biāo)的轉(zhuǎn)換過(guò)程中，由于計(jì)算大地緯度口時(shí)用到大地高IV,而計(jì)算大地高時(shí)又

5、需要用到大地緯度口.因此不能直接由空間直角坐標(biāo)計(jì)算出大地坐標(biāo)，而需要采用迭代計(jì)算的方法。具體計(jì)算時(shí)，可先根據(jù)下式求出大地緯度口的初值：然后利用該初值來(lái)求出H N的初值，再利用所求出的H和N初值再次求出B值如此反復(fù)，直至 求出的及日、V收斂為止。4. 算例本文根據(jù)以上公式在 Microsoft VC+6.0 環(huán)境下編寫(xiě)了一段程序（見(jiàn)附錄）。算例中的坐標(biāo)采 用的是武漢大 學(xué)信息學(xué)部友誼廣場(chǎng) 上的某點(diǎn)的大地坐標(biāo) 作為已知值，然后經(jīng)過(guò)轉(zhuǎn)換函數(shù) CRDGEODETICtoCRDCARTESEAN（pcg, pcc, dSemiMajorAxis,dFlat ning 把大地坐標(biāo)轉(zhuǎn)換為空間直角 坐標(biāo)得到坐

6、標(biāo)X、Y、乙由得到的空間直角坐標(biāo) X、Y、Z,經(jīng)過(guò)轉(zhuǎn)換函數(shù) vpx -.%. CRDCARTESIANtoCRDGEODETIC（pcc,pcg,dSemiMajorAxis,dFlat ning把空間直角坐標(biāo)還原成空間大廠(chǎng)一一_.v'X地坐標(biāo)，計(jì)算結(jié)果如下圖所示：計(jì)算結(jié)果從上圖可以看出結(jié)果比較滿(mǎn)意，高程和精度基本能完全還原，而緯度還原后有較大的誤差，在 '7I 測(cè)量中這種誤差不允許的，需要修改算法，完善結(jié)果?？赡芤鸬脑蛴锌赡苁怯捎诰暥扔?jì)算公式 并不完善，還有可能是由于計(jì)算機(jī)的截?cái)嗾`差引起的，還要找時(shí)間繼續(xù)修改、完善。I'5. 心得體會(huì)這次編程自認(rèn)為很簡(jiǎn)單，但真動(dòng)手

7、自己親自編寫(xiě)，還是或多或少遇到了一些問(wèn)題，并分析問(wèn)題， 最終解決問(wèn)題。雖然這次作業(yè)很簡(jiǎn)單， 但經(jīng)過(guò)自己這樣一步一步的編寫(xiě)出來(lái)，還是有很多收獲，加強(qiáng)了運(yùn)用VC+編寫(xiě)程序的能力，也充分認(rèn)識(shí)到了學(xué)習(xí) VC+勺重要性，更找到了自己的一些缺點(diǎn)與 不足。6. 附錄（程序源代碼）#in clude<iostream>#in clude<cmath>using n amespace std;#define M_PI 3.1415926typedef struct tagCRDCARTESIANdouble x;double y;double z;CRDCARTESIAN;typedef

8、 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;廣：j、Mi n=(i nt)(DMS-Deg)*100);Sec=(DMS-Deg)*100-Mi n)*100;, I I 'I2i-r

9、、. | N"X'*Rad=(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 (Rad<0)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(Rad&l

10、t;0)*DMS=-*DMS;return;bool CRDCARTESIANtoCRDGEODETIC(PCRDCARTESPANPCRDGEODETIC pcg,double dSemiMajorAxis,double dFlattening)double B0,R,N;double B_,L_;-double X=pcc->x;double Y=pcc->y;7廠(chǎng)丁 丨double Z=pcc->z;R=sqrt(X*X+Y*Y); , I I '2f z'B0=ata n2(Z,R);while (1)N=dSemiMajorAxis/sqrt(1.0

11、-dFlatte nin g*(2-dFlatte nin g)*si n(B0)*si n(B0);B_=ata n2(Z+N*dFlatte nin g*(2-dFlatte nin g)*si n(B0),R);if(fabs(B_-B0)<1.0e-10)break;B0=B_;L_=ata n2(Y,X);pcg->height=R/cos(B_)-N;RAD_DMS(B _,&pcg->latitude);RAD_DMS(L_,&pcg->lo ngitude);return true;pcc,doublebool CRDGEODETICto

12、CRDCARTESEAN(PCRDGEODEcg(PCRDCARTESIAN dSemiMajorAxis,double dFlattening)double N;double B_,L_;double B=pcg->latitude;double L=pcg-> Ion gitude;double H=pcg->height;| /A二二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-&g

13、t;x=(N+H)*cos(B_)*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_);7 ci '、/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=

14、(114*3600+21*60+20.51)/3600;pcg->height=41;pcg-> Ion gitude=rad;double dSemiMajorAxis=6378137;cout<<std:fixed;coutvv"轉(zhuǎn)換前已知的大地坐標(biāo)："<<endl;cout«"H（高度）="<vpcg->height<v" "<v"L（經(jīng)度）="<<pcg->longitude<<" "&l

15、t;<"B（ 緯度）="<<pcg->latitude<<e ndl;CRDGEODETICtoCRDCARTESEAN（pcg, pcc, dSemiMajorAxis,dFlat nin g）;cout«e ndl;coutvv"轉(zhuǎn)換后的空間直角坐標(biāo)："cout<<"x="<<pcc->x<<" "vv"y="vvpcc->yv<" "vv"z="vvpcc->z<v" "<<e ndl;CRDCARTESIANtoCRDGEODETIC（pcc,pcg,dSemiMajorAxis,dFlat nin g）;coutvve ndl;coutvv"由轉(zhuǎn)換后的空間直角坐標(biāo) x,y,z還原成大地坐標(biāo)："<<endl;co