Koordinatentransformation: ECI-ECEF
Verfasst: 18.08.2006 - 11:41
So ich sitzt jetzt schon den zweiten Tag dran und versuche, ECI-Koordinaten in ECEF-Koordinaten umzurechen.
Ich habe die Daten von GPS-Satteliten mit GPS-Sekunden (z.B. 838388720) und Quartanionen. Nun habe ich schon die Quartanionen in das ECI-System umgerechnet, das dürfte passen. Anschließend rechne ich die ECI Koordinaten in ECEF Koordinaten um und zum Schluß die ECEF-Koordinaten in Spherische Koordinaten(long, lat).
Der Fehler müsste sich bei der Umrechnung ECI->ECEF eingeschlichen haben, ich denke mal in der Zeitumrechnung (Earth Rotation), hat mir da vielleicht jemand ne Idee? Oder andere Formeln. Denn die Latidude stimmt ja, nur der Punkt bewegt sich auf der Longitudinalen, dies sollte er aber nicht, da dieser Punkt eine Stadt ist.
hier mein MatLab-Program:
load P.dat %Datei mit GPS-Positionsdaten
P
for k=1:20
sec=P(k,1)
q0=P(k,2)
q1=P(k,3)
q2=P(k,4)
q3=P(k,5)
%E -> unit vector along axis of rotation
sigma=2*acos(q0)
e1=q1/(sin (sigma/2))
e2=q2/(sin (sigma/2))
e3=q3/(sin (sigma/2))
E=[e1;e2;e3]
% Antenna in inertial system:
%Ai=[0;0;1]
% Antenna in body fixed system
Ab=[0;0;1]
% Transformation Matrix
R=[(2*q0^2-1+2*q1^2) 2*(q1*q2+q0*q3) 2*(q1*q3-q0*q2); 2*(q1*q2-q0*q3) (2*q0^2-1+2*q2^2) 2*(q3*q2+q0*q1); 2*(q1*q3+q0*q2) 2*(q3*q2-q0*q1) (2*q0^2-1+2*q3^2)]
% Rotation from body fixed system to inertial system,
% E defined in body fixed system
%Ai=R*Ab
%Ab=R'*Ai
% E defined in inertial system
%Ab=R*Ai
%Ai*R'*Ab
Ai= R*Ab
x1=Ai(1,1)
x2=Ai(2,1)
x3=Ai(3,1)
x=[x1;x2;x3]
% ECI -> ECEF
% Julian Date: the number of days since noon January 1, 4713 BC including the fraction of day. It thus provides a continuous time scale.
JD=2444244.5+(7*1024)+(1/86400)*sec
% Greenwich mean Sidereal Time:
d=JD-2451545.0
T=(JD-2451545.0)/36525
GMST=4.894961215+(6.300388101*d)+(6.770708127e-006)*T^3-(0.4508185458e-009)*T^3
%1) Precession -> D= transformation matrix for precession
%zetaA=0.01118086087*T+(0.1463555541e-005)*T^2+(0.8725676634e-007)*T^3
%zA=0.01118086087*T+(0.5307158406e-005)*T^2+(0.8825063438e-007)*T^3
%ThetaA=(0.9717173458e-002)*T-(0.2068457571e-005)*T^2-(0.2028121073e-006)*T^3
% Transformation Matrix
%Rz1=[cos(-zA) sin(-zA) 0;-sin(-zA) cos(-zA) 0;0 0 1]
%Ry=[cos(ThetaA) 0 -sin(ThetaA);0 1 0;sin(ThetaA) 0 cos(ThetaA)]
%Rz2=[cos(-zetaA) sin(-zetaA) 0;-sin(-zetaA) cos(-zetaA) 0;0 0 1]
%D=Rz1*Ry*Rz2
%2) Nutation -> C= transformation matrix for nutation
%my=2.181661565-33.75464593*T
%eta=3.506366467+1256.660924*T
%epsilon0=0.4090928041-(0.2269655249e-003)*T-(0.2860400719e-008)*T^2+(0.8789672039e-008)*T^3
%deltaepsilon=(4.537856054e-005)*cos(my)+(0.3490658504e-005)*cos(eta)
epsilon=epsilon0+deltaepsilon
deltapsi=(-8.377580412e-005)*sin(my)-(0.6981317008e-005)*sin(eta)
% Transformation Matrix
%Rx1=[1 0 0;0 cos(-epsilon) sin(-epsilon);0 -sin(-epsilon) cos(-epsilon)]
%Rz=[cos(-deltapsi) sin(-deltapsi) 0;-sin(-deltapsi) cos(-deltapsi) 0;0 0 1]
%Rx2=[1 0 0;0 cos(epsilon0) sin(epsilon0);0 -sin(epsilon0) cos(epsilon0)]
%C=Rx1*Rz*Rx2
%3) Earth Rotation -> B= transformation matrix for Earth`s rotation
% Greenwich Mean Sidereal Time
GAST=GMST+deltapsi*cos(epsilon)
B=[cos(GAST) sin(GAST) 0;-sin(GAST) cos(GAST) 0;0 0 1]
%4) Polar Motion -> A= transformation matrix for polar motion
% The polar motion is neglected since it`s effect is insignificant
% Transformation Matrix T: ECEF - ECI
%T=B*C*D
T=B
a=T'*x
% Umrechnung von ECEF->Sphere (x,y,z gegeben)
x=a(1,1)
y=a(2,1)
z=a(3,1)
r=(x^2+y^2+z^2)^0.5
p=(x^2+y^2)^0.5
% sphärische Koordinaten
lat1=(acos(z/r))*180/pi
lat=90-lat1
long=(acos(x/p))*180/pi
kord(:, k)=[lat long]
Ich habe die Daten von GPS-Satteliten mit GPS-Sekunden (z.B. 838388720) und Quartanionen. Nun habe ich schon die Quartanionen in das ECI-System umgerechnet, das dürfte passen. Anschließend rechne ich die ECI Koordinaten in ECEF Koordinaten um und zum Schluß die ECEF-Koordinaten in Spherische Koordinaten(long, lat).
Der Fehler müsste sich bei der Umrechnung ECI->ECEF eingeschlichen haben, ich denke mal in der Zeitumrechnung (Earth Rotation), hat mir da vielleicht jemand ne Idee? Oder andere Formeln. Denn die Latidude stimmt ja, nur der Punkt bewegt sich auf der Longitudinalen, dies sollte er aber nicht, da dieser Punkt eine Stadt ist.
hier mein MatLab-Program:
load P.dat %Datei mit GPS-Positionsdaten
P
for k=1:20
sec=P(k,1)
q0=P(k,2)
q1=P(k,3)
q2=P(k,4)
q3=P(k,5)
%E -> unit vector along axis of rotation
sigma=2*acos(q0)
e1=q1/(sin (sigma/2))
e2=q2/(sin (sigma/2))
e3=q3/(sin (sigma/2))
E=[e1;e2;e3]
% Antenna in inertial system:
%Ai=[0;0;1]
% Antenna in body fixed system
Ab=[0;0;1]
% Transformation Matrix
R=[(2*q0^2-1+2*q1^2) 2*(q1*q2+q0*q3) 2*(q1*q3-q0*q2); 2*(q1*q2-q0*q3) (2*q0^2-1+2*q2^2) 2*(q3*q2+q0*q1); 2*(q1*q3+q0*q2) 2*(q3*q2-q0*q1) (2*q0^2-1+2*q3^2)]
% Rotation from body fixed system to inertial system,
% E defined in body fixed system
%Ai=R*Ab
%Ab=R'*Ai
% E defined in inertial system
%Ab=R*Ai
%Ai*R'*Ab
Ai= R*Ab
x1=Ai(1,1)
x2=Ai(2,1)
x3=Ai(3,1)
x=[x1;x2;x3]
% ECI -> ECEF
% Julian Date: the number of days since noon January 1, 4713 BC including the fraction of day. It thus provides a continuous time scale.
JD=2444244.5+(7*1024)+(1/86400)*sec
% Greenwich mean Sidereal Time:
d=JD-2451545.0
T=(JD-2451545.0)/36525
GMST=4.894961215+(6.300388101*d)+(6.770708127e-006)*T^3-(0.4508185458e-009)*T^3
%1) Precession -> D= transformation matrix for precession
%zetaA=0.01118086087*T+(0.1463555541e-005)*T^2+(0.8725676634e-007)*T^3
%zA=0.01118086087*T+(0.5307158406e-005)*T^2+(0.8825063438e-007)*T^3
%ThetaA=(0.9717173458e-002)*T-(0.2068457571e-005)*T^2-(0.2028121073e-006)*T^3
% Transformation Matrix
%Rz1=[cos(-zA) sin(-zA) 0;-sin(-zA) cos(-zA) 0;0 0 1]
%Ry=[cos(ThetaA) 0 -sin(ThetaA);0 1 0;sin(ThetaA) 0 cos(ThetaA)]
%Rz2=[cos(-zetaA) sin(-zetaA) 0;-sin(-zetaA) cos(-zetaA) 0;0 0 1]
%D=Rz1*Ry*Rz2
%2) Nutation -> C= transformation matrix for nutation
%my=2.181661565-33.75464593*T
%eta=3.506366467+1256.660924*T
%epsilon0=0.4090928041-(0.2269655249e-003)*T-(0.2860400719e-008)*T^2+(0.8789672039e-008)*T^3
%deltaepsilon=(4.537856054e-005)*cos(my)+(0.3490658504e-005)*cos(eta)
epsilon=epsilon0+deltaepsilon
deltapsi=(-8.377580412e-005)*sin(my)-(0.6981317008e-005)*sin(eta)
% Transformation Matrix
%Rx1=[1 0 0;0 cos(-epsilon) sin(-epsilon);0 -sin(-epsilon) cos(-epsilon)]
%Rz=[cos(-deltapsi) sin(-deltapsi) 0;-sin(-deltapsi) cos(-deltapsi) 0;0 0 1]
%Rx2=[1 0 0;0 cos(epsilon0) sin(epsilon0);0 -sin(epsilon0) cos(epsilon0)]
%C=Rx1*Rz*Rx2
%3) Earth Rotation -> B= transformation matrix for Earth`s rotation
% Greenwich Mean Sidereal Time
GAST=GMST+deltapsi*cos(epsilon)
B=[cos(GAST) sin(GAST) 0;-sin(GAST) cos(GAST) 0;0 0 1]
%4) Polar Motion -> A= transformation matrix for polar motion
% The polar motion is neglected since it`s effect is insignificant
% Transformation Matrix T: ECEF - ECI
%T=B*C*D
T=B
a=T'*x
% Umrechnung von ECEF->Sphere (x,y,z gegeben)
x=a(1,1)
y=a(2,1)
z=a(3,1)
r=(x^2+y^2+z^2)^0.5
p=(x^2+y^2)^0.5
% sphärische Koordinaten
lat1=(acos(z/r))*180/pi
lat=90-lat1
long=(acos(x/p))*180/pi
kord(:, k)=[lat long]