衛(wèi)星導(dǎo)航定位算法與程序設(shè)計(jì)-實(shí)驗(yàn)報(bào)告

級(jí)測(cè)繪工程專業(yè)之吉白夕凡創(chuàng)作創(chuàng)作時(shí)間：二零二一年六月三十天衛(wèi)星導(dǎo)航定位算法與法式設(shè)計(jì)實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱：衛(wèi)星導(dǎo)航基本法式設(shè)計(jì)班級(jí)：學(xué)號(hào)：姓名：實(shí)驗(yàn)時(shí)間：6月28日~6月30中國(guó)礦業(yè)大年夜學(xué)目錄實(shí)驗(yàn)一時(shí)空基準(zhǔn)變換2一、實(shí)驗(yàn)?zāi)康?二、實(shí)驗(yàn)內(nèi)容2三、實(shí)驗(yàn)過(guò)程創(chuàng)作時(shí)間：二零二一年六月三十天四、實(shí)驗(yàn)感想??實(shí)驗(yàn)二????RINEX文件讀寫??一、實(shí)驗(yàn)?zāi)康?二、實(shí)驗(yàn)內(nèi)容7三、實(shí)驗(yàn)過(guò)程7實(shí)驗(yàn)三衛(wèi)星軌道計(jì)算12一、實(shí)驗(yàn)?zāi)康?2二、實(shí)驗(yàn)內(nèi)容12三、實(shí)驗(yàn)過(guò)程12四、實(shí)驗(yàn)感想15創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天實(shí)驗(yàn)一時(shí)空基準(zhǔn)變換一、實(shí)驗(yàn)?zāi)康?、加深對(duì)時(shí)空系統(tǒng)及其之間變換關(guān)系的理解2、掌握常常使用時(shí)空基準(zhǔn)之間的變換模型與軟件實(shí)現(xiàn)3、每人自力達(dá)成實(shí)驗(yàn)規(guī)定的內(nèi)容二、實(shí)驗(yàn)內(nèi)容本實(shí)驗(yàn)內(nèi)容包含：內(nèi)容一：編程實(shí)現(xiàn)GPS起點(diǎn)1980年1月6日0時(shí)對(duì)應(yīng)的儒略日內(nèi)容二：編程實(shí)現(xiàn)11月27日對(duì)應(yīng)的GPS周數(shù)與一周內(nèi)的秒數(shù)內(nèi)容三：在WGS84橢球的條件下,編程實(shí)現(xiàn)傍邊央子午線為117度時(shí),計(jì)算高斯坐標(biāo)x=3548910.811290287,y=179854.6172135982對(duì)應(yīng)的經(jīng)緯度坐標(biāo)？?jī)?nèi)容四：WGS84橢球下,概略x=2408000;y=4698000;z=3566000處的地平坐標(biāo)系坐標(biāo)為:e=704.8615;n=114.8683;u=751.9771的點(diǎn)對(duì)應(yīng)的直角坐標(biāo)為幾多？三、實(shí)驗(yàn)過(guò)程1.針對(duì)第一、二部份內(nèi)容：1.1解決思路：先成立”TimeStruct.h的頭”文件,將格里高利歷、GPS時(shí)間結(jié)構(gòu)、儒略日時(shí)間構(gòu)造共構(gòu)造體的方式放在里面；在成立“TimeTr的”頭文件,成立類“CTimeTr”,創(chuàng)立變量“GPSTime”、“Time”、”JulDay”而,且申明函數(shù)“TIME2JUL”、“TIME2GTIME”等,用這些函數(shù)分別實(shí)現(xiàn)所需要的變換.1.2詳細(xì)的實(shí)現(xiàn)函數(shù)：“TIME2JUL”函數(shù)：doubleCTimeTr::TIME2JUL( )//TIMETime,JULIANDAY&JulDay{doublem,y;doubleD;//h=Time.byHour+Time.byMinute/60.0+Time.dSecond/3600.00;if(Time.byMonth<=2){創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天y=Time.wYear1;m=Time.byMonth+12;}else{y=Time.wYear;m=Time.byMonth;}D=floor(365.25*(y+4716))+floor(30.6001*(m+1))+Time.byDay+Time.byHour/24.01537.5;JulDay.lDay=int(D);JulDay.lSecond=Dint(JulDay.lDay);return0;}“TIME2GTIME”：voidCTimeTr::TIME2GTIME( ){doubleJD;longm,y;intWN;doubleWsecend;//UT=Time.byHour+Time.byMinute/60.0+Time.dSecond/3600.00;if(Time.byMonth<=2){y=Time.wYear1;m=Time.byMonth+12;}else{y=Time.wYear;m=Time.byMonth;}JD=int(365.25*y)+int(30.6001*(m+1))+Time.byDay+Time.byHour/24.0+1720981.5;WN=floor((JD2444244.5)/7.0);GpsTime.lWeek=WN;Wsecend=(JD2444244.57*WN)*604800;GpsTime.lSecond=Wsecend;}1.3實(shí)驗(yàn)結(jié)果：創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天針對(duì)第三部份內(nèi)容：2.1解決思路：運(yùn)用實(shí)驗(yàn)指導(dǎo)書(shū)中供給的matlab高斯反算的代碼,進(jìn)行解算；將高斯反算的公式直接輸成matlab代碼,繞后在函數(shù)“function[B,L]=gauss_fansuanx,y,L0）”中,將坐標(biāo)x=3548910.811290287,y=179854.6172135982,L0=117,帶入函數(shù)的坐邊,即可獲取所需要的經(jīng)緯度.2.2主要函數(shù)的代碼：function[B,L]=gauss_fansuan(x,y,L0)a=6378137;f=1/298.257223563;b=aa*f;c=a^2/b;e=sqrt(a^2b^2)/a;e1=sqrt(a^2b^2)/b;Beta0=1(3/4)*e1^2+(45/64)*e1^4(175/256)*e1^6+(11025/16384)*e1^8;Beta2=Beta01;Beta4=(15/32)*e1^4(175/384)*e1^6+(3675/8192)*e1^8;Beta6=(35/96)*e1^6+(735/2048)*e1^8;Beta8=(315/1024)*e1^8;B0=x/(c*Beta0);aa0=(a*cos(B0))/sqrt(1e^2*sin(B0)^2);l0=y/aa0;N=a*sqrt(1e^2*sin(B0)^2);創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天t=tan(B0);in=e1*cos(B0);a1=N*cos(B0);a2=(1/2)*N*sin(B0)*cos(B0);a3=(1/6)*N*cos(B0)^3*(1t^2+in^2);a4=(1/24)*N*sin(B0)*cos(B0)^3*(5t^2+9*in^2+4*in^4);a5=(1/120)*N*cos(B0)^5*(518*t^2+t^4+14*in^258*in^2*t^2);a6=(1/720)*N*sin(B0)*cos(B0)^5*(6158*t^2+t^4);F_xB=(c*Beta2+(c*Beta4+(c*Beta6+c*Beta8*cos(B0)^2)*cos(B0)^2)*cos(B0)^2)*sin(B0)*cos(B0);F_xBl=a2*l0^2+a4*l0^4+a6*l0^6;F_yBl=a3*l0^3+a5*l0^5;B1=(xF_xBF_xBl)/(c*Beta0);aa1=(a*cos(B1))/sqrt(1e^2*(sin(B1))^2);l1=(yF_yBl)/aa1;B0=B1;aa0=aa1;l0=l1;F_xB=(c*Beta2+(c*Beta4+(c*Beta6+c*Beta8*cos(B0)^2)*cos(B0)^2)*cos(B0)^2)*sin(B0)*cos(B0);F_xBl=a2*l0^2+a4*l0^4+a6*l0^6;F_yBl=a3*l0^3+a5*l0^5;B1=(xF_xBF_xBl)/(c*Beta0);aa1=(a*cos(B1))/sqrt(1e^2*(sin(B1)^2));l1=(yF_yBl)/aa1;endL=rad2deg(l1)+L0;B=rad2deg(B1);創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天四、實(shí)驗(yàn)感想本次試驗(yàn)是花時(shí)間許多的一次實(shí)驗(yàn),對(duì)于時(shí)間變換的部份所有都是自己下手將matlab代碼寫成“C++”的類,進(jìn)行實(shí)現(xiàn)的.此中碰到的較大年夜的困難是儒略日向UTC變換的部份,這部份的函數(shù)步伐許多,重點(diǎn)是在一開(kāi)始的時(shí)間構(gòu)造里面,各時(shí)間各部份的數(shù)據(jù)種類大年夜多界說(shuō)的是“int型”的,但是在進(jìn)行計(jì)算的時(shí)候有許多的小數(shù),需要用到浮點(diǎn)型的函數(shù),這部份用了許多的時(shí)間.在做這個(gè)實(shí)驗(yàn)的時(shí)候,第一天花了時(shí)間主假如變換代碼,使法式?jīng)]有缺點(diǎn),可以正常的運(yùn)轉(zhuǎn)出來(lái),體現(xiàn)黑框框,但是還只有個(gè)別功能可以用,可以運(yùn)轉(zhuǎn)出正確的結(jié)果；次日時(shí)間主假如花在改正函數(shù)上頭,可以使所寫的功能都能運(yùn)轉(zhuǎn)出正確的結(jié)果.經(jīng)過(guò)做時(shí)間變換的實(shí)驗(yàn),使自己發(fā)生了第一次親身編寫“C++”代碼的經(jīng)歷,并且所有缺點(diǎn)的解決所有都是自己解決,收獲好多.實(shí)驗(yàn)二RINEX文件讀寫一、實(shí)驗(yàn)?zāi)康?、深入認(rèn)識(shí)RINEX文件格式2、進(jìn)一步提升MATLAB法式設(shè)計(jì)能力3、掌握N文件、O文件、SP3文件的基本讀寫技巧二、實(shí)驗(yàn)內(nèi)容本實(shí)驗(yàn)內(nèi)容包含：1、任選IGS站,下載N文件、O文件與SP3文件；2、編程實(shí)現(xiàn)N文件讀入,并采用中文標(biāo)明出主要參數(shù)的名稱及作用；4、編程實(shí)現(xiàn)O文件讀入,并采用中文標(biāo)明出主要參數(shù)的名稱及作用；5、編程實(shí)現(xiàn)SP3文件讀入,并采用中文標(biāo)明出主要參數(shù)的名稱及作用；三、實(shí)驗(yàn)過(guò)程1、針對(duì)第一部份內(nèi)容：編程實(shí)現(xiàn)N文件讀入,并采用中文標(biāo)明出主要參數(shù)的名稱及作用創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天1.1、解決思路：依據(jù)“GPSeasy開(kāi)”源代碼供給的函數(shù),依據(jù)實(shí)驗(yàn)要求讀取了N文件的內(nèi)容,先用“rinexe函”數(shù),將N文件讀取成“eph.dat文”件,而后再用“get_eph函”數(shù)將“eph.dat文”件讀取成“Eph矩”陣,此矩陣中包含了N文件中的數(shù)據(jù),在最后用“fprintf函數(shù)”將所需要的數(shù)據(jù)輸出成”.TXT”文件即可.1.2、主要函數(shù)代碼:“get_eph函”數(shù)：functioneph=get_eph(ephemeridesfile)fide=fopen(ephemeridesfile);[eph,count]=fread(fide,Inf,'double');noeph=count/22;eph=reshape(eph,22,noeph);“rinexe函”數(shù)（部份）：functionrinexe(ephemerisfile,outputfile)fide=fopen(ephemerisfile);head_lines=0;while1head_lines=head_lines+1;line=fgetl(fide);answer=findstr(line,'ENDOFHEADER');if~isempty(answer),break;end;end;head_lines主函數(shù)中輸出結(jié)果得函數(shù)部份：af0=data(19);%衛(wèi)星中差M0=data(3);roota=data(4);deltan=data(5);ecc=data(6);omega=data(7);cuc=data(8);cus=data(9);crc=data(10);crs=data(11);i0=data(12);創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天idot=data(13);toe=data(18);af1=data(20);%對(duì)所要輸出的參數(shù)賦值fprintf(fid,'\n衛(wèi)星編號(hào)：%d\n衛(wèi)星鐘差：%d\n平近點(diǎn)角距：%d\n軌道長(zhǎng)半軸的平方根：%d\n均勻運(yùn)動(dòng)修正量：%d\n軌道偏愛(ài)率：%d\n近地點(diǎn)角距：%d\n緯度幅角的余弦調(diào)解項(xiàng)更正的振幅',prn,af0,M0,roota,deltan,ecc,omega,cuc);fprintf(fid,'緯度幅角的正弦調(diào)解項(xiàng)更正的振幅：%d\n軌道半徑的余弦調(diào)解項(xiàng)更正的振幅：%d\n軌道半徑的正弦調(diào)解項(xiàng)更正的振幅：%d\n軌道傾角：%d\n軌道傾角改動(dòng)率：%d\n星歷參照時(shí)刻：%d\n',cus,crc,crs,i0,idot,toe)fclose(fid);1.3、輸出結(jié)果2、針對(duì)第二部份內(nèi)容：編程實(shí)現(xiàn)O文件讀入,并采用中文標(biāo)明出主要參數(shù)的名稱及作用；2.1、實(shí)現(xiàn)思路：經(jīng)過(guò)matlab的函數(shù)“fopen讀”取O文件,獲取O文件的指針,經(jīng)過(guò)“anheader函”數(shù)將文件中的接收機(jī)大年夜致地點(diǎn)”approx_XYZ1”天,線的偏移值”ant_delta1觀”測(cè),值種類“Obs_types1等”讀入成為matlab的矩陣,而后經(jīng)過(guò)循環(huán),利用“grabdata函”數(shù)將所需要的歷元的觀察文件挨次輸出來(lái),最后經(jīng)過(guò)“fprintf函數(shù)”,將所需要的數(shù)據(jù)挨次打印出來(lái).2.2、主要函數(shù)：“anheader函”數(shù)：function[Obs_types,ant_delta,ifound_types,approx_XYZ]=anheader(file)fid=fopen(file,'rt');eof=0;ifound_types=0;Obs_types=[];ant_delta=[];approx_XYZ=[];while1%Gobblingtheheaderline=fgetl(fid);answer=findstr(line,'ENDOFHEADER');if~isempty(answer),break;end;if(line==1),eof=1;break;end;answer=findstr(line,'ANTENNA:DELTAH/E/N');創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天if~isempty(answer)fork=1:3[delta,line]=strtok(line);del=str2num(delta);ant_delta=[ant_deltadel];end;endanswer=findstr(line,'APPROXPOSITIONXYZ');if~isempty(answer)fork=1:3[app_XYZ,line]=strtok(line);del=str2num(app_XYZ);approx_XYZ=[approx_XYZdel];end;endanswer=findstr(line,'#/TYPESOFOBSERV');if~isempty(answer)[NObs,line]=strtok(line);NoObs=str2num(NObs);fork=1:NoObs[ot,line]=strtok(line);Obs_types=[Obs_typesot];end;ifound_types=1;end;end;%fclose(fid);“grabdata函”數(shù)：functionObs=grabdata(fid,NoSv,NoObs)%GRABDATAPositionedinaRINEXfileataselectedepochreadsobservationsofNoSvsatellitesgloballinObs=zeros(NoSv,NoObs);ifNoObs<=5%ThiswilltypicalbeTurboSIIdataforu=1:NoSvlin=fgetl(fid);fork=1:NoObsObs(u,k)=str2num(lin(2+16*(k1):16*k2));endendelse%ThiswilltypicalbeZ12dataObs=Obs(:,[12345]);%Wecancelthelasttwocolumns6and7NoObs=5;foru=1:NoSv創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天lin=fgetl(fid);lin_doppler=fgetl(fid);fork=1:NoObs%%1ifisempty(str2num(lin(1+16*(k1):16*k2)))==1,Obs(u,k)=nan;else%Obs(u,k)=str2num(lin(1+16*(k1):16*k2));endObs(u,NoObs)=str2num(lin(65:78));endendend四、實(shí)驗(yàn)感想這部份實(shí)驗(yàn)是我在以前做的,以前自己有看過(guò)“gps_easy有”關(guān)的代碼,看過(guò)相關(guān)的“N文件”“O文件”讀寫函數(shù),并且學(xué)會(huì)了怎樣調(diào)用這些函數(shù),對(duì)里面的輸出量有了一點(diǎn)的認(rèn)識(shí),所以我自己的主要工作就是運(yùn)用了“fprintf函數(shù)”,將讀取到matlab中的矩陣寫入TXT文檔中,這部份工作量不是很大年夜,但較存心義.實(shí)驗(yàn)三衛(wèi)星軌道計(jì)算一、實(shí)驗(yàn)?zāi)康?、進(jìn)一步熟習(xí)N文件的讀入2、掌握開(kāi)普勒參數(shù)計(jì)算衛(wèi)星軌道的過(guò)程3、編程實(shí)現(xiàn)采用廣播星歷計(jì)算衛(wèi)星軌道4、掌握MATLAB函數(shù)調(diào)用步伐二、實(shí)驗(yàn)內(nèi)容本實(shí)驗(yàn)內(nèi)容包含：1、調(diào)試時(shí)間變換函數(shù),熟習(xí)內(nèi)容,備主函數(shù)調(diào)用2、調(diào)試廣播星歷導(dǎo)航文件的讀入法式,備主函數(shù)調(diào)用3、依據(jù)衛(wèi)星地點(diǎn)計(jì)算公式編寫主函數(shù),同時(shí)調(diào)用時(shí)間變換、星歷讀取等的子函數(shù)來(lái)共同達(dá)成衛(wèi)星地點(diǎn)的計(jì)算,最后輸出計(jì)算結(jié)果4、理清法式各模塊的功能構(gòu)造創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天三、實(shí)驗(yàn)過(guò)程1、實(shí)驗(yàn)思路：在老師供給的“SPP”文件中,直接利用衛(wèi)星地點(diǎn)計(jì)算函數(shù),進(jìn)行衛(wèi)星地點(diǎn)的計(jì)算,將利用“Gps.cpp文”件中的”GetGpsPosition函”數(shù),利用此中的迭代求解衛(wèi)星地點(diǎn)部份,用“cout直”接將衛(wèi)星迭代后的地點(diǎn)直接輸出,由于星歷文件中有許多的星歷,所以利用循環(huán)語(yǔ)句,將求解出來(lái)的衛(wèi)星地點(diǎn)挨次輸出出來(lái).2、主要函數(shù)boolCGps::GetGpsPosition( ){GPSTIMEts;GPSTIMEtr;GPSTIMEts0;GPSTIMEoGTime;GPSTIMEnGTOC;vector<GpsSendPosition>SendSignPosition;GpsPosGpsPTemp;GpsSendPositionSdSignPoTemp;intnTheFitPoint=0;PositionpTemp;GpsPosition.clear( );if(oData.size( )==0){strErr=("PRN=%d沒(méi)有對(duì)應(yīng)星歷");returnfalse;}if(nData.size( )==0){returnfalse;}cout<<fixed;for(inti=0;i<oData.size( );i++){double*VX=newdouble[4];////for(intj=0;j<oData[i].ObserveInfo.ObserveSatSum;j++){//TimeToGpsTime(oData[i].ObserveInfo.ObserveTime,oGTime);tr=oGTime;//創(chuàng)作時(shí)間：二零二一年六月三十天創(chuàng)作時(shí)間：二零二一年六月三十天//FindTheBestFitTime(oData[i].ObserveInfo.ObserveTime,oData[i].oObserveData[j].PRN,nTheFitPoint);if(nTheFitPoint==1){//stringjuge;//juge.("PRN=%d沒(méi)?有?D對(duì)?應(yīng)?|星?歷¤¨2",oData[i].oObserveData[j].PRN);//AfxMessageBox(juge);continue;}ts.lWeek=tr.lWeek;ts0.lWeek=tr.lWeek;//ts.lSecond=tr.lSecond0.075;do{ts0.lSecond=ts.lSecond;SatellitePosition(nData[nTheFitPoint],ts0,pTemp);ts.lSecond=tr.lSecondsqrt(pow((pTemp.XXx0),2)+pow((pTemp.YYy0),2)+pow((pTemp.ZZz0),2))/c;}while(fabs(ts0.lSecondts.lSecond)>1e007);////earthrot(tr.lSecondts.lSecond,