通過matlab計算衛星位置


衛星星歷是描述衛星運動軌道的信息。也可以說衛星星歷就是一組對應某一時刻的軌道參數及其變率。有了衛星星歷就可以計算出任意時刻的衛星位置及其速度。GPS衛星星歷分為預報星歷和后處理星歷。預報星歷又稱廣播星歷。

GPS廣播星歷參數共有16個,其中包括1個參考時刻,6個對應參考時刻的開普勒軌道參數和9個反映攝動力影響的參數。這些參數通過GPS衛星發射的含有軌道信息的導航電文傳遞給用戶。

 

1.星歷參考時刻 : 

2.長半軸平方根:

3.偏心率:

4.參考歷元下平近點角:

5.近地點角距:

6.軌道傾角:

7.本周初始歷元的升交點赤經:

8.平運動差(由精密星歷計算得到的衛星平均角速度與按給定參數計算所得的平均角速度之差):n

9.軌道傾角變化率(弧度/秒):

10.升交點赤經變化率(弧度/秒):

11.緯度幅角的余弦調和項改正的振幅(弧度):

12.緯度幅角的正弦調和項改正的振幅(弧度):

13.軌道半徑的余弦調和項改正的振幅(m):

14.軌道半徑的正弦調和項改正的振幅(m):

15.軌道傾角的余弦調和項改正的振幅(弧度):

16.軌道傾角的正弦調和項改正的振幅(弧度):

以下為matlab程序:

%the homework 1 of chapater 1
%student :Taylen
%time :2018/9/25
%程序功能:根據所提供的星歷參數,計算此衛星在信號發射時刻t( GPS時間)239050.7223s
%時的時空位置

%GPS廣播星歷參數共有16個,其中包括1個參考時刻,6個對應參考時刻的開普勒軌道參數
%和9個反映攝動力影響的參數。這些參數通過GPS衛星發射的含有軌道信息的導航電文傳遞給用戶。
%時間參數
pra1 = 244800; %te---星歷參考時刻,即星歷表參考歷元(s)

%開普勒六參數
pra2 = -1.064739758; %M0---按參考歷元te計算的平近點角(弧度)
pra3 = 0.005912038265; %e---軌道偏心率
pra4 = 5153.65531; %sqrt(a)長半軸平方根
pra5 = -1.717457876; %w0---近地點角距
pra6 = 0.9848407943; %i0---按參考歷元計算的軌道傾角(弧度)
pra7 = 1.038062244; %Ω0---本周初始歷元的升交點赤經(弧度)

%軌道攝動九參數
pra8 = 4.249105564e-9; % Δ𝑛---平運動差(弧度)
pra9 = 7.422851197e-51; % dot_i---軌道傾角變化率
pra10 = -8.151768125e-9; % dot_Omega---升交點赤經變化率
pra11 = 3.054738045e-7; % Cuc---升交點角距的調和改正項振幅
pra12 = 2.237036824e-6; % Cus---升交點角距的調和改正項振幅
pra13 = 350.53125; % Crc---衛星地心距的調和改正項振幅
pra14 = 2.53125; % Crs---衛星地心距的調和改正項振幅
pra15 = -8.381903172e-8; % Cic---軌道傾角的調和改正項振幅
pra16 = 8.940696716e-8; % Cis---軌道傾角的調和改正項振幅

t=239050.7223; %此衛星在信號發射時刻t ( GPS時間)

miu = 3.986005e14; % 地心引力常數
we = 7.2921151467e-5; % 地球自轉角速度
F = -4.442807633e-10; % Constant, [sec/(meter)^(1/2)]

% 計算歸化時間
te = pra1;
tk = t - te;
% 計算衛星的平均角速度
a = pra4.^2;
n0 = sqrt(miu/a.^3);
n = n0 + pra8;
% 計算平近點角
Mk = pra2 + n * tk;
% 計算偏近點角Ek
%迭代計算:相鄰兩次計算差之絕對值值<1e-15時結束迭代計算,Ek的迭代初始值為0
Ek0=0;Ek=Mk;
while abs(Ek0-Ek)>1e-15
Ek0 = Ek;
Ek=Mk + pra3 * sin(Ek0);
end
% 計算衛星鍾差相對論校正值 Compute relativistic correction term

% 計算真近點角:取值在(-pi,pi]
v1=sqrt(1-pra3^2)*sin(Ek);
v2=cos(Ek)-pra3;
vk=atan(v1/v2);
if v1>0 & v2<0
Vk = pi + vk;%第二象限
elseif v1<0 & v2<0
Vk = vk - pi;%第三象限
end
% 計算升交點角距
Faik = Vk + pra7;
% 計算升交點角距改正值
Sigmauk=pra11*cos(2*Faik)+pra12*sin(2*Faik);
% 計算衛星地心向徑改正值
Sigmark=pra13*cos(2*Faik)+pra14*sin(2*Faik);
% 計算衛星軌道傾角改正值
Sigmaik=pra15*cos(2*Faik)+pra16*sin(2*Faik);
% 修正升交點角距
uk=Faik+Sigmauk;
% 修正衛星地心相徑
rk=a*(1-pra3*cos(Ek))+Sigmark;
% 修正衛星軌道傾角
ik=pra6+Sigmaik+tk*pra9;
% 計算衛星軌道平面直角坐標
xk=rk*cos(uk);
yk=rk*sin(uk);
% 計算升交點赤經
Omegak=pra5+(pra10-we)*tk-we*te;
% 將衛星坐標由軌道平面轉換到ECEF坐標系
Xk=xk*cos(Omegak)-yk*cos(ik)*sin(Omegak);
Yk=xk*sin(Omegak)+yk*cos(ik)*cos(Omegak);
Zk=yk*sin(ik);




免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM