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

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

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

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

4、地面點P沿橢球法線到橢球面的距離叫大地高。 大地坐標(biāo)坐標(biāo)系中，P點的位 置用L,B表示。如果點不在橢球面上，表示點的位置除L,B外，還要附加另一參 數(shù)一一大地高Ho空間大地坐標(biāo)系3 .空間直角坐標(biāo)與大地坐標(biāo)間的轉(zhuǎn)換3.1 大地坐標(biāo)轉(zhuǎn)換為空間直角坐標(biāo)將同一坐標(biāo)系下的大地坐標(biāo)（B、L、H）轉(zhuǎn)換成空間直角坐標(biāo)（X、Y、Z）的轉(zhuǎn)換 公式為：X = : NZ=N(l-/)+HsinB =式中N為卯酉圈的半徑，a為參考橢球的長半軸；b為參考橢球的短半軸; e為參考橢球的第一偏心率；并且有若點在橢球面上，則大地高 H=0,上式可簡化為:X = N8sBecnLY-NcosBsinLZ = N/l-e%inB

5、3.2 空間直角坐標(biāo)轉(zhuǎn)換為空間大地坐標(biāo)將同一坐標(biāo)系下的空間直角坐標(biāo)（X、Y、Z）轉(zhuǎn)換為空間大地坐標(biāo)（B、L、H）的公式為:Z = arctanR)展由g - Z(N+H)W+r2)畫(>/)+siaB在使用上式進(jìn)行空間直角坐標(biāo)到大地坐標(biāo)的轉(zhuǎn)換過程中，由于計算大地緯度 口時用到大地高IV ,而計算大地高時又需要用到大地緯度口. 因此不能直接由空 間直角坐標(biāo)計算出大地坐標(biāo)，而需要采用迭代計算的方法。具體計算時，可先根 據(jù)下式求出大地緯度口的初值：2B = arctan(一)卡產(chǎn)然后利用該初值來求出H、N的初值，再利用所求出的 M口也值再次求出B 值.如此反復(fù)，直至求出的及日、IV收斂為止。4

6、 .算例本文根據(jù)以上公式在Microsoft VC+6.0環(huán)境下編寫了一段程序(見附錄)。 算例中的坐標(biāo)采用的是武漢大學(xué)信息學(xué)部友誼廣場上的某點的大地坐標(biāo)作為已 知值，然后經(jīng)過轉(zhuǎn)換函數(shù) CRDGEODETICtoCRDCARTESEAN(pcg, pcc, dSemiMajorAxis,dFlatning)把大地坐標(biāo)轉(zhuǎn)換為空間直角坐標(biāo)得到坐標(biāo)X、Y、Z。由得到的空間直角坐標(biāo)X、Y、Z,經(jīng)過轉(zhuǎn)換函數(shù)CRDCARTESIANtoCRDGEODETIC(pcc,pcg,dSemiMajorAxis,dFlatning)把空間直 角坐標(biāo)還原成空間大地坐標(biāo)，計算結(jié)果如下圖所示：計算結(jié)果從上圖可以看出結(jié)果

7、比較滿意，高程和精度基本能完全還原，而緯度還原后 有較大的誤差，在測量中這種誤差不允許的，需要修改算法，完善結(jié)果?？赡芤?起的原因有可能是由于緯度計算公式并不完善， 還有可能是由于計算機的截斷誤 差引起的，還要找時間繼續(xù)修改、完善。5 .心得體會這次編程自認(rèn)為很簡單，但真動手自己親自編寫，還是或多或少遇到了一些 問題，并分析問題，最終解決問題。雖然這次作業(yè)很簡單，但經(jīng)過自己這樣一步 一步的編寫出來，還是有很多收獲，加強了運用 VC+編寫程序的能力，也充分 認(rèn)識到了學(xué)習(xí)VC+的重要性，更找到了自己的一些缺點與不足。6 .附錄(程序源代碼)#include<iostream>#incl

8、ude<cmath>using namespace std;#define M_PI 3.1415926typedef struct tagCRDCARTESIANdouble x;double y;double z;CRDCARTESIAN;typedef CRDCARTESIAN *PCRDCARTESIAN;typedef struct tagCRDGEODETICdouble longitude;double latitude;double height;CRDGEODETIC;typedef CRDGEODETIC *PCRDGEODETIC;void DMS_RAD(d

9、ouble DMS,double *Rad)int Deg,Min;double Sec;Deg=(int)DMS;Min=(int)(DMS-Deg)*100);Sec=(DMS-Deg)*100-Min)*100;*Rad=(Deg+Min/60.0+Sec/3600.0)/180.0*M_PI;return;void RAD_DMS(double Rad,double *DMS)(int Deg,Min;double Sec;double AR,AM;AR=Rad;if (Rad<0)AR=-Rad;AR=AR+1.0e-10;AR=AR*180.0/M_PI;Deg=(int)A

10、R;AM=(AR-Deg)*60.0;Min=(int)AM;Sec=(AM-Min)*60;*DMS=Deg+Min/100.0+Sec/10000.0;if(Rad<0)*DMS=-*DMS;return;boolCRDCARTESIANtoCRDGEODETIC(PCRDCARTESIANpcc,PCRDGEODETIC pcg,double dSemiMajorAxis,double dFlattening)(double B0,R,N;double B_,L_;double X=pcc->x;double Y=pcc->y;double Z=pcc->z;R=

11、sqrt(X*X+Y*Y);B0=atan2(Z,R);while (1) (N=dSemiMajorAxis/sqrt(1.0-dFlattening*(2-dFlattening)*sin(B0)*sin(B0);B_=atan2(Z+N*dFlattening*(2-dFlattening)*sin(B0),R);if(fabs(B_-B0)<1.0e-10)break;B0=B_;L_=atan2(Y,X);pcg->height=R/cos(B_)-N;RAD_DMS(B_,&pcg->latitude);RAD_DMS(L_,&pcg->lo

12、ngitude);return true;boolCRDGEODETICtoCRDCARTESEAN(PCRDGEODETICpcg,PCRDCARTESIAN pcc,double dSemiMajorAxis,double dFlattening)(double N;double B_,L_;double B=pcg->latitude;double L=pcg->longitude;double H=pcg->height;DMS_RAD(B,&B_);DMS_RAD(L,&L_);N=dSemiMajorAxis/sqrt(1.0-dFlattenin

13、g*(2-dFlattening)*sin(B_)*sin(B_);pcc->x=(N+H)*cos(B_)*cos(L_);pcc->y=(N+H)*cos(B_)*sin(L_);pcc->z=(N*(1.0-dFlattening*(2-dFlattening)+H)*sin(B_);return true;void main()PCRDCARTESIAN pcc=new CRDCARTESIAN;PCRDGEODETIC pcg=new CRDGEODETIC;/B=30.31.40.23L=114.21.20.51 h=41double rad;rad=(30*36

14、00+31*60+40.23)/3600;pcg->latitude=rad;rad=(114*3600+21*60+20.51)/3600;pcg->height=41;pcg->longitude=rad;double dSemiMajorAxis=6378137;double dFlatning=1/298.257223563;cout<<std:fixed;cout<<”轉(zhuǎn)換前已知的大地坐標(biāo)："<<endl;cout<<"H(高 度尸"<<pcg->height<&l

15、t;""<<"L( 經(jīng)pcc,度尸"<<pcg->longitude<<""<<"B(緯度尸"<<pcg->latitude<<endl;CRDGEODETICtoCRDCARTESEAN(pcg,dSemiMajorAxis,dFlatning);cout<<endl;cout<<"轉(zhuǎn)換后的空間直角坐標(biāo)："cout<<"x="<<pcc->x<<""<<"y="<<pcc->y<<""<<"z="<<pcc->z<<""<<endl;CRDCARTESIANtoCRDGEODETIC(pcc,pcg,dSemiMajorAxis,dFlatning);cout<<endl;cout<<"由轉(zhuǎn)換后的空間直