Bodenlänge/-breite unter einem Satelliten (kartesische Koordinaten) zu einer bestimmten Epoche

Von dem Punkt an, an dem Sie sind, müssen Sie:

1. Berechnen Sie die XYZ-Koordinaten Ihrer Bodenstation
2. Finden Sie die Linie im Raum zwischen dem Satelliten und der Bodenstation
3. Finden Sie die XYZ-Koordinaten des Punktes am Schnittpunkt dieser Linie mit dem Ellipsoid, das durch alle Punkte mit Höhen von 300 km definiert ist
4. Finden Sie den Breiten- und Längengrad dieses Punktes.

Für (1) können Sie diese Matlab-Funktion verwenden:

function [X,Y,Z] = ll2xyz(lat,lon,alt)
    b=6356752.3141;
    a=6378137.0;
    lat=lat*pi/180;
    %transformation between geografic and geocentric latitude
    gclat=atan(b^2*tan(lat)/a^2);
    lon=lon*pi/180;
    R=sqrt(1./(tan(gclat).^2/b^2 + 1/a^2));
    Z=R.*tan(gclat);
    Z=Z+alt.*sin(lat);
    R=R+alt.*cos(lat);
    X=R.*cos(lon);
    Y=R.*sin(lon);

Für (2) und (3) reicht diese Funktion aus

function [X Y Z] = rectaxelip(x1,x2,y1,y2,z1,z2,alt)
%rectaxelip calculates the point XYZ were a straight line that pass by points P1
%and P2, intersects a ellipsoid alt meters bigger than WGS84 
%ellispsoid (on each axis)

%The director cosines of the line are 
    d=sqrt((x1-x2).^2.*(y1-y2).^2.*(z1-z2).^2);
    L=(x1-x2)./d;
    M=(y1-y2)./d;
    N=(z1-z2)./d;
    %its parametric  for is
    %x = x1 + L * t
    %y = y1 + M * t
    %z = z1 + N * t
    %Now we look for the intersection
    %the ellipsoid would have axes bp and ap

    b=6356752.3141; 
    a=6378137.0;
    bp=b+alt; 
    ap=a+alt;

    %The equation of such ellipsod would be
    %  x^2/ap^2+y^2/ap^2+z^2/bp^2=1
    %subtituting and reorganizng you get a quadratic equation like
    %   A x^2 + B x + C = 0
    %with
    A=(L./ap).^2+(M./ap).^2+(N./bp).^2;
    B=2*((x1.*L./(ap.^2))+(y1.*M./(ap.^2))+(z1.*N./(bp.^2)));
    C=((x1./ap).^2+(y1./ap).^2+(z1./bp).^2)-1;
    ta=(-B+sqrt(B.^2 - (4*A.*C)))./(2*A);
    tb=(-B-sqrt(B.^2 - (4*A.*C)))./(2*A);

    %then
    x1a = x1 + L .* ta;
    y1a = y1 + M .* ta;
    z1a = z1 + N .* ta;
    x1b = x1 + L .* tb;
    y1b = y1 + M .* tb;
    z1b = z1 + N .* tb;
    %Now the distance from each point to the satellite is
    da=sqrt((x1a-x2).^2+(y1a-y2).^2+(z1a-z2).^2);
    db=sqrt((x1b-x2).^2+(y1b-y2).^2+(z1b-z2).^2);
    ESta=da<db;
    t=ta.*ESta+tb.*(~ESta);

    %then
    X = x1 + L .* t;
    Y = y1 + M .* t;
    Z = z1 + N .* t;

Und für (4) brauchen Sie nur eine Funktion, die Lat-Lon-Koordinaten aus XYZ-Positionen berechnet, und diese würde es tun:

function [lat,lon,alt] = xyz2ll(x,y,z1)
%For some reason the error spikes at +-45
a=6356752.3141; 
b=6378137.0;

signos=sign(z1);
z1=-abs(z1);

lon=atan2(y,x)*180/pi; %Longitude
d1=sqrt(x.^2+y.^2);
z2=-sqrt(1./((1/(a^2))+(1./(b^2*(z1./d1).^2))));
d2=z2./(z1./d1);
p2=-a*d2./(b^2*sqrt(1-(d2.^2/(b^2)))); %slope of the tangent to the ellipse in the point
for i=1:5
    %points in the palne
    dp=(z1-z2-(d1./p2)+(p2.*d2))./(p2-(1./p2));
    zp=(p2.*dp)+z2-(p2.*d2);
    %points in the ellipse
    z2=-sqrt(1./((1/(a^2))+(1./(b^2*(zp./dp).^2))));
    d2=z2./(zp./dp);
    p2=-a*d2./(b^2*sqrt(1-(d2.^2/(b^2))));
end
lat=-90-(atan(p2)*180/pi);   %geographic latitude
lat=abs(lat).*signos;
ecuador=find(abs(z1)<0.001);
lat(ecuador)=0;
alt=sqrt(d1.^2+z1.^2)-sqrt(d2.^2+z2.^2);
alt(ecuador)=d1(ecuador)-a;

Und dann sind Sie bereit und haben die Bodenposition des Punktes, an dem das GPS-Signal, das die Basisstation erreicht, die 300 km Höhe überquert.

Ich habe diese Funktionen selbst geschrieben und ausgiebig getestet, damit ich weiß, dass sie funktionieren, aber seien Sie sich bewusst, dass die letzte Funktion trotz ihrer sehr guten Genauigkeit aus einem Grund, den ich nicht viel untersucht habe, große Fehler bei Breitengraden zu erzeugen scheint von +-45°.

Das Letzte, was zu beachten ist. Wenn Sie mehr Genauigkeit wünschen und die Berechnung der XYZ-Koordinaten der GPS-Satelliten vermeiden möchten, können Sie die sp3-Ephemeridendateien verwenden. Sie können hier einfach heruntergeladen werden . Es gibt eine Datei pro Tag, und jede Datei enthält die XYZ-Positionen jedes Satelliten, die alle 15 Minuten tabelliert werden. Um jederzeit Positionen zu erhalten, bietet die Spline-Interpolation eine sehr gute Lösung, ziemlich genau die gleiche wie die empfohlene Polynom-Interpolation nach Neville.