GPS系統的前身為美軍研制的一種子午儀衛星定位系統(Transit），1958年研制，64年正式投入使用。該系統用5到6顆衛星組成的星網工作，每天最多繞過地球13次，并且無法給出高度信息，在定位精度方面也不盡如人意。然而，子午儀系統使得研發部門對衛星定位取得了初步的經驗，并驗證了由衛星系統進行定位的可行性，為GPS系統的研制埋下了鋪墊。由于衛星定位顯示出在導航方面的巨大優越性及子午儀系統存在對潛艇和艦船導航方面的巨大缺陷。美國海陸空三軍及民用部門都感到迫切需要一種新的衛星導航系統[13]。 1973年12月 ,美國國防部批準它的陸海空三軍聯合研制新的衛星導航系統: NAVSTAR/GPS。它是英文“Navigation Satellite Timing and Ranging/Global Positioning System” 的縮寫詞。其意為 “衛星測時測距導航/全球定位系統”,簡稱 GPS。這個系統向有適當接受設備的全球范圍用戶提供精確、 連續的三維位置和速度信息 ,并且還廣播一種形式的世界協調時(U TC) 。通過遍布全球的（21+3）GPS導航衛星，向全球范圍內的用戶全天候提供高精度的導航、跟蹤定位和授時服務。目前，GPS已在地形測量，交通管理，導航，野外勘探，空間宇宙學等諸多領域得到了廣泛的應用[11]。 目前全球共有4大GPS系統，分別是： 美國 GPS ，由美國國防部于 20 世紀 70 年代初開始設計、 研制 ，于1993 年全部建成。1994 年 ，美國宣布在 10 年內向全世界免費提供 GPS使用權 ，但美國只向外國提供低精度的衛星信號。 歐盟 “伽利略”，1999 年 歐洲提出計劃 ,準備發射 30 顆衛星 ，組成 “伽利略” 衛星定位系統。 俄羅斯 “格洛納斯”，尚未部署完畢。始于上世紀 70年代，需要至少 18 顆衛星才能確保覆蓋俄羅斯全境;如要提供全球定位服務 ，則需要 24 顆衛星。 中國“北斗”2003 年我國北斗一號建成并開通運行 ，不同于 GPS， “北斗” 的指揮機和終端之間可以雙向交流。四川大地震發生后 ,北京武警指揮中心和四川武警部隊運用 “北斗” 進行了上百次交流。北斗二號系列衛星今年起將進入組網高峰期 ，預計在 2015 年形成由三十幾顆衛星組成的覆蓋全球的系統。

1.2 本課題研究的意義和方法

GPS系統是一個很龐大的系統，包含了天文，地理，計算機，電磁學，通信學，信息學等等。通過本文對GPS的學習研究，最重要的還是要學習其原理：衛星運動原理；衛星定位原理；衛星跟蹤原理等等。通過基礎原理的學習，一方面，可以使我們更進一步的理解衛星運動，定位的實現方法；通過仿真，進一步了解簡單定位的方法及其在仿真平臺上的實現途徑；另一方面，也可以培養我們自學的能力，訓練仿真模擬的技巧和方法。

至今，基本上完成了課題的要求，通過不斷的注入既定參數，可以更加詳細，直觀的理解基本的定位原理和實現方法！

1.3 GPS前景

GPS導航定位以其定位精度高、觀測時間短、測站間無需通視、可提供三維坐標、操作簡便、全天候作業、功能多、應用廣泛等特點著稱。

用GPS信號可以進行海、空和陸地的導航、導彈的制導、大地測量和工程測量的精密定位、時間的傳遞和速度的測量等。對于測繪領域，GPS衛星定位技術已經用于建立高精度的全國性的大地測量控制網，測定全球性的地球動態參數；用于建立陸地海洋大地測量基準，進行高精度的海島陸地聯測以及海洋測繪；用于檢測地球板塊運動狀態和地殼形變；用于工程測量，成為建立城市與工程控制網的主要手段。用于測定航空航天攝影瞬間相機位置，實現僅有少量的地面控制或無地面控制的航測快速成圖，導致地理信息系統、全球環境遙感監測的技術革命[4]。

目前，GPS、GLONASS、INMARSAT等系統都具備了導航定位功能，形成了多元化的空間資源環境。這一多元化的空間資源環境，促使國際民間形成了一個共同的策略，即一方面對現有系統充分利用，一方面積極籌建民間GNSS系統，待2011年左右，GNSS純民間系統建成，全球?將形成GPS/GLONASS/GNSS三足鼎立之勢，才能從根本上擺脫對單一系統的依賴，形成國際共有、國際共享的安全資源環境。世界才可以將衛星導航作為單一導航手段的最高應用境界。國際民間的這一策略，反過來又影響和迫使美國對其GPS使用政策作出更開放的調整。多元化的空間資源環境的確立，給GPS的發發展應用創造了一個前所未有的良好的國際環境。

?

?

第二章 GPS測量原理

GPS導航系統的基本原理是測量出已知位置的衛星到用戶接收機之間的距離 ,然后綜合多顆衛星的數據就可知道接收機的具體位置。要達到這一目的 ,衛星的位置可以根據星載時鐘所記錄的時間在衛星星歷中查出。而用戶到衛星的距離則通過紀錄衛星信號傳播到用戶所經歷的時間 ,再將其乘以光速得到(由于大氣層電離層的干擾 ,這一距離并不是用戶與衛星之間的真實距離 ,而是偽距( PR) :當GPS衛星正常工作時 ,會不斷地用 1 和 0 二進制碼元組成的偽隨機碼(簡稱偽碼)發射導航電文。GPS系統使用的偽碼一共有兩種 ,分別是民用的 C/ A 碼和軍用的 P( Y)碼。C/ A 碼頻率 1. 023MHz ,重復周期一毫秒 ,碼間距 1 微秒 ,相當于 300m; P 碼頻率10. 23MHz ,重復周期266. 4 天 ,碼間距0. 1 微秒 ,相當于 30m。而 Y碼是在 P碼的基礎上形成的 ,保密性能更佳。

GPS導航系統衛星部分的作用就是不斷地發射導航電文。然而 ,由于用戶接受機使用的時鐘與衛星星載時鐘不可能總是同步 ,所以除了用戶的三維坐標 x、 y、 z外 ,還要引進一個Δt 即衛星與接收機之間的時間差作為未知數 ,然后用 4 個方程將這 4個未知數解出來。所以如果想知道接收機所處的位置 ,至少要能接收到 4 個衛星的信號。

2.1 偽距測量的原理

GPS定位采用的是被動式單程測距。它的信號發射書機由衛星鐘確定，收到時刻是由接收機鐘確定，這就在測定的衛星至接收機的距離中，不可避免地包含著兩臺鐘不同步的誤差和電離層、對流層延遲誤差影響，它并不是衛星與接受機之間的實際距離，所以稱之為偽距。

偽距定位法是利用全球衛星定位系統進行導航定位的最基本的方法，其基本原理是：在某一瞬間利用GPS接收機同時測定至少四顆衛星的偽距，根據已知的衛星位置和偽距觀測值，采用距離交會法求出接收機的三維坐標和時鐘改正數。它的優點是速度快、無多值性問題，利用增加觀測時間可以提高定位精度；缺點是測量定位精度低，但足以滿足部分用戶的需要。

?

?

用戶位置

地心位置

?

圖5-5 用戶可見的衛星分布

SatellitePosition =?????????????????????

??1.0e+004 *

??? 1.7746??? 1.7572??? 0.7365??? 0.0001

?? -1.2161??? 0.9732??? 2.1091??? 0.0001

??? 2.2883??? 0.2975?? -1.4132??? 0.0001

?? -2.3882?? -0.8836??? 0.7869???????? 0

?? -0.3681?? -2.5255??? 0.7012???????? 0

?? -0.1323??? 2.7059?? -0.0000??? 0.0001

?? -1.5424??? 0.4273??? 2.1704??? 0.0001

??? 1.4582??? 0.4926?? -2.1144???????? 0

?? -2.4090??? 0.6583??? 0.7365???????? 0

??? 1.0759??? 1.1599?? -2.1757?????? ??0

??? 1.3324?? -1.8178??? 1.4392??? 0.0001

?? -1.3225??? 1.7688?? -1.3835???????? 0

??? 1.2127?? -0.9774??? 2.1091??? 0.0001

?? -1.8766??? 0.3838?? -1.9058???????? 0

?? -2.4302?? -0.9295?? -0.7516???????? 0

??? 0.4443??? 1.8187?? -1.8463???????? 0

?? -0.6718??? 2.4562?? -0.7869???????? 0

??? 1.4000?? -1.3073??? 1.9058??? 0.0001

?? -1.4992??? 0.4223?? -2.1091???????? 0

?? -0.8706?? -2.0297?? -1.3835???????? 0

?? -1.7722?? -0.6021?? -1.9217???????? 0

?? -0.2481?? -1.5812?? -2.1704???????? 0

?? -0.5421??? 1.7915??? 1.9217??? 0.0001

??? 1.9376?? -1.5756?? -0.7365??? 0.0001

仿真程序四：用可見衛星計算用戶位置

程序見附錄

多于四顆的衛星用最小二乘法逼近計算。對于前面的假設用戶位置（6400，3352，5410）進行計算，結果如下：

?

圖5-6 用戶位置的計算

由圖13可見，根據最小二乘法的原理，計算出來了用戶在一個時刻的位置坐標，并把它表示在了這個天球坐標系當中。

calculaterecord =

? 1.0e+003 *

?????? 0? ???????0???????? 0

??? 7.3399??? 3.8857??? 6.1388

??? 6.3034??? 3.2976??? 5.3346

??? 6.3914??? 3.3470??? 5.4034

??? 6.4009??? 3.3525??? 5.4107

??? 6.4001??? 3.3520??? 5.4101

先后經過了六次的迭代算法，吧用戶的計算位置一步一步的逼近了用戶的實際位置，根據部同的精度要求，我們運算的量的大小也有不同。這點，可以根據程序當中的參數Error的設定而有所出入。

?

圖5-7 用戶位置的計算

通過對用戶仿真計算的進一步放大，可以很清楚地看到：在每一次迭代的過程當中，我們都實時的把每一步迭代的值也表示在了同一個坐標軸里面：目的就是很直觀的反應出我們基于最小二乘法的偽距算法的主要思路。如下圖所示：

?

圖5-8 用戶位置的計算

?? 上圖15當中的黑色箭頭從上到下一次表示的是第一次到第六次計算出來的用戶的位置，第七箭頭表示的用戶的實際位置（用白色的柱體表示）。

?

?

程序一（繪制衛星的軌道平面）代碼

function drawsatelliteorbit

a=26560;% 衛星軌道的長半軸.

e=0.02;%e是橢圓的偏心率.

E=[0:0.1:2*pi];

x=a*(cos(E)-e);

y=a*sqrt((1-e^2))*sin(E);

z=0*E;

DtoR=2*pi/360;

A1=[32.8? 92.8? 152.8? 212.6? 272.8? 332.8];%衛星星座數據.

for k=1:6

??? A=A1(k)*DtoR;%升交點的經度

??? B=55*DtoR;%軌道的傾角

??? C=pi/100;%近地點幅角

??? %總共有6個衛星軌道平面

??? R1=[cos(A)? -sin(A)? 0;

??????? sin(A)? cos(A)?? 0;

??????? 0??????? 0?????? 1;];

??? R2=[1???????? 0?????? 0;

??????? 0??????? cos(B)? -sin(B);

??????? 0??????? sin(B)?? cos(B);];

??? R3=[cos(C)??? -sin(C)??? 0;

??????? sin(C)??? cos(C)???? 0;

??????? 0????????? 0??????? 1;];

??? L1=length(E);

??? R312=R1*R2*R3;

??? Ans=R312*[x;y;z;];

??? x1=Ans(1,:);

??? y1=Ans(2,:);

??? z1=Ans(3,:);

??? plot3c(x1,y1,z1,k);

??? boxplot3(0,0,0,200,200,200,7);

??? hold on;

? axis equal;

? axis off;

end

程序二（單顆衛星不同時刻的動態仿真）代碼

clear; clc;close all;

DtoR=2*pi/360;

jiaostep=360/24*DtoR;

?j=0;

??? a=26560;% 衛星軌道的長半軸.

e=0.02;%e是橢圓的偏心率.

E=[0:0.1:2*pi];

x=a*(cos(E)-e);

y=a*sqrt((1-e^2))*sin(E);

z=0*E;

??? drawearth(0);

hold on;

for time=1:12

??? axis on;

A1=32.8 ;%衛星星座數據.

??? A=A1*DtoR;%升交點的經度

??? B=55*DtoR;%軌道的傾角

??? C=pi/100;%近地點幅角

??? R1=[cos(A)? -sin(A)? 0;

??????? sin(A)? cos(A)?? 0;

??????? 0??????? 0?????? 1;];

??? R2=[1???????? 0?????? 0;

??????? 0??????? cos(B)? -sin(B);

??????? 0??????? sin(B)?? cos(B);];

??? R3=[cos(C)??? -sin(C)??? 0;

??????? sin(C)??? cos(C)???? 0;

??????? 0????????? 0??????? 1;];

??? L1=length(E);

??? R312=R1*R2*R3;

??? Ans=R312*[x;y;z;];

??? x1=Ans(1,:);

??? y1=Ans(2,:);

??? z1=Ans(3,:);

??? plot3c(x1,y1,z1,2);

??? hold on;

? axis equal;

? axis on;

? grid on;

ctable=10;

??????? A=A1*DtoR;

??????? B=55*DtoR;

??????????? C=ctable*DtoR+time*2*pi/24;???????? %近地點幅角

????????????? x=a*(cos(C)-e);

y=a*sqrt((1-e^2))*sin(C);

z=0*C;

?R1=[cos(A)??? -sin(A)?? ????0;

????? sin(A)??? cos(A)?????? 0;

??????? 0?????????? 0?????? 1;];

??? R2=[1???????? 0??????? 0;

??????? 0??????? cos(B)? -sin(B);

??????? 0??????? sin(B)?? cos(B);];

??? R3=[cos(C)??? -sin(C)???? 0;

??????? sin(C)??? cos(C)????? 0;

??????? 0???????? ?0??????? 1;];

??? L1=length(E);

??? R312=R1*R2*R3;

??? Ans=R312*[x;y;z;];

??? x1=Ans(1,:);

??? y1=Ans(2,:);

??? z1=Ans(3,:);

??? drawsatellite(x1,y1,z1,6);??

??? M(time)=getframe;

??? M(time+1) = getframe;

?? end

axis on;???????????

程序三（衛星在某個時刻的全軌道平面的分布和可見衛星）代碼

clear;clc;close all;

a=26560;????????????????????????????????????????? % 衛星軌道的長半軸.

e=0.02;???

temp=0;%便于把衛星的指標點順序放到矩陣指標當中。

%e是橢圓的偏心率.

E=[0:0.1:2*pi];

x=a*(cos(E)-e);

y=a*sqrt((1-e^2))*sin(E);

z=0*E;

timenow=0;%單位是小時

global SatellitePosition

global ttum

SatellitePosition=ones(24,4);

figure(1);

drawearth(0);??

ttum=ones(24,4);

hold on;

DtoR=2*pi/360;

A1=[32.8? 92.8? 152.8? 212.6? 272.8? 332.8];

drawsatelliteorbit;????? ??????????????????????????????????????????????????????

ctable=[10 50 160 260;?? ???????????????????????????????%平均近地角

??????? 80 180 220 320;

??????? 10 130 250 340;

??????? 50 150 170 300;

??????? 100 210 310 340;

??????? 120 140 240 350;];

??? simple=1;

???? for k=1:6

??????? A=A1(k)*DtoR;??????????????????????????????? %升交點經度

???? ???B=55*DtoR;?????????????????????????????????? %軌道傾角

??????? for m=1:4

??????????? C=ctable(k,m)*DtoR+timenow*2*pi/24;???????? %近地點幅角

????????????? x=a*(cos(C)-e);

y=a*sqrt((1-e^2))*sin(C);

z=0*C;

?R1=[cos(A)??? -sin(A)?????? 0;

????? sin(A)??? cos(A)??? ???0;

??????? 0?????????? 0?????? 1;];

R2=[1???????? 0??????? 0;

??????? 0??????? cos(B)? -sin(B);

??????? 0??????? sin(B)?? cos(B);];

R3=[cos(C)??? -sin(C)???? 0;

??????? sin(C)??? cos(C)????? 0;

??????? 0????????? 0??????? 1;];

??? L1=length(E);

??? R312=R1*R2*R3;

??? Ans=R312*[x;y;z;];

??? x1=Ans(1,:);

??? y1=Ans(2,:);

??? z1=Ans(3,:);

?? drawsatellite(x1,y1,z1,k);?????????????????????????????

??? temp=temp+1;

??? SatellitePosition(temp,:)=[x1 y1 z1 1];

??? ttum(temp,:)=[x1 y1 z1 1];

??? hold on;

????? ??end

???? end

?%用余弦定理計算地心到用戶和用戶到衛星證件的夾角要是兩者證件的夾角小于90°,就認為改衛星是不可見的.

???? earthcenterpos=[0 0 0];

???? userposition=[6400? 3352? 5410];?????????????? %假設一個用戶的位置

???? for k=1:24

???????? temp=SatellitePosition(k,1:3)-userposition;

???????? Dist1=temp*temp';

???????? temp=userposition-earthcenterpos;

???????? Dist2=temp*temp';

???????? temp=SatellitePosition(k,1:3)-earthcenterpos;

???????? Dist3=temp*temp';

???????? jiajiao=acos((Dist2+Dist3-Dist1)/2/sqrt(Dist3)/sqrt(Dist2));

???????? if(jiajiao>=pi/2)

?????? ??????SatellitePosition(k,4)=0;

???????????? %ttum(k,:)=SatellitePosition(k,1:3);

???????? end

???? end

?? %? ttum????

???? figure(2)

??? drawearth(0);????????????????????????????????????

???? hold on;

???? drawsatelliteorbit;?

???????? j=1;

????? for k=1:6

????????? for m=1:4

???????? if(SatellitePosition(j,4)==1)

???????????? tempx=SatellitePosition(j,1);

???????????? tempy=SatellitePosition(j,2);

???????????? tempz=SatellitePosition(j,3);

???????????? drawsatellite(tempx,tempy,tempz,k);???????????

? ???????end

????????? j=j+1;

????????? end

????? end

???? tempx=userposition(1);

???? tempy=userposition(2);

???? tempz=userposition(3);

???? cube=100;

???? boxplot3(tempx,tempy,tempz,cube,cube,cube,0); %標記出用戶的位置

程序四（用可見衛星計算用戶的位置）代碼（mian3）

%用余弦定理計算地心到用戶和用戶到衛星證件的夾角要是兩者證件的夾角小于90°,就認為該衛星是不可見的.

???? earthcenterpos=[0 0 0];

???? userposition=[6400? 3352? 5410];?????????????? %假設一個用戶的位置

???? for k=1:24

???? SatellitePosition(k,:)=ttum(k,:);

???? end

???? for k=1:24

???????? temp=SatellitePosition(k,1:3)-userposition;

???????? Dist1=temp*temp';

???????? temp=userposition-earthcenterpos;

???????? Dist2=temp*temp';

???????? temp=SatellitePosition(k,1:3)-earthcenterpos;

???????? Dist3=temp*temp';

???????? jiajiao=acos((Dist3+Dist2-Dist1)/2/sqrt(Dist3)/sqrt(Dist2));

? ???????if(jiajiao>=pi/2)

???????????? SatellitePosition(k,4)=0;

???????? end

???? end

???? figure(1)

??? drawearth(0);??????????????????????????????????

???? hold on;

???? drawsatelliteorbit;?

????? tempx=userposition(1);

???? tempy=userposition(2);

??? ?tempz=userposition(3);

???????? cube=55;

???????? boxplot3(tempx,tempy,tempz,cube,cube,cube,7);????????????

???????? hold on; %先畫出用戶的位置，以便和后面的做比較

????? j=1;

????? for k=1:6

????????? for m=1:4

???????? if(SatellitePosition(j,4)==1)

???????????? tempx=SatellitePosition(j,1);

???????????? tempy=SatellitePosition(j,2);

???????????? tempz=SatellitePosition(j,3);

???????????? drawsatellite(tempx,tempy,tempz,k);????????????

???????? end

????????? j=j+1;

????????? end

????? end

????? tempx=userposition(1);

??? ?tempy=userposition(2);

???? tempz=userposition(3);

???? hold on;

???? SatellitePosition=[SatellitePosition;userposition 0];????

?? %用變量caluserpos來記錄迭代計算的結果.

???? caluserpos=calculateuserposition(SatellitePosition);

???? [m,n]=size(caluserpos);

???? for(k=3:m);

???????? tempx=caluserpos(k,1);

???????? tempy=caluserpos(k,2);

???????? tempz=caluserpos(k,3);

???????? cube=50;

???????? boxplot3(tempx,tempy,tempz,cube,cube,cube,0);????????????

???????? hold on;

???? end

以下是一個和上面程序配套使用的功能函數，把它命名為calculateuserposition.m。并且要和上面的程序放在一個文件夾當中。

function caluserposition=calculateuserposition(SatellitePosition)

global SatellitePosition?????????????

?%單獨運行此程序時候，為調用前面程序而得的satelliteposotion結果，應該聲明satelliteposotion為全局變量，即應該在satelliteposition前面加上global

r1=6400;?????????? ?????????????????????????????????????????%地球半徑

c=300000;

deltat=0.00001;

satelliteposnew=ones(1,3);

vissatnum=0;

calculateok=1;

for k=1:24

??? if(SatellitePosition(k,4)==1)

??????? vissatnum=vissatnum+1;

??????? satelliteposnew=[satelliteposnew;SatellitePosition(k,1:3)];

??? end%if

end%for

satelliteposnew(1,:)=[];

if(vissatnum<4)

??? calculateok=0;

??? caluserposition=[0 0 0 0];

??? return

end

prange=ones(1,vissatnum);

userpos=SatellitePosition(25,1:3);

for n=1:vissatnum??? prange(1,n)=sqrt((satelliteposnew(n,:)-userpos)*(satelliteposnew(n,:)-userpos)')+c*deltat;

end

calculaterecord=[0??? 0?? 0];

xyz0=[0 0 0];

deltat0=0;

wxyz=satelliteposnew;

Error=1000;

computertime=0;

while((Error>0.00001)&(computertime<1000))

??? computertime=computertime+1;

??? r=ones(1,vissatnum);

??? for n=1:vissatnum

??????? r(1,n)=sqrt((wxyz(n,:)-xyz0)*(wxyz(n,:)-xyz0)')+deltat0*c;

??? end

??? deltap=r-prange;

??? A=ones(vissatnum,3);

??? for n=1:vissatnum

??? A(n,:)=(wxyz(n,:)-xyz0)./r(1,n);

??? end%for

??? H=[A ones(vissatnum,1)];

??? deltax=inv(H'*H)*H'*deltap';

??? tempdeltax=deltax(1:3)

??? Error=max(abs(tempdeltax));

??? xyz0=xyz0+deltax(1:3,:)';

??? if(computertime<1000)

??????? calculaterecord=[calculaterecord;xyz0];

??? end

??? deltat0=deltax(4,1)/(-c);

end%while

calculaterecord;

if(computertime==1000)

??? caluserposition=[99999999999 99999999999 999999999999];

else

??? caluserposition=[xyz0;calculaterecord];

end

總結

以上是生活随笔為你收集整理的GPS卫星运动及定位matlab仿真的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。