amflab.m

这是国外关于卫星导航方面一书的源代码

   % widerspruch = double-difference(obs) - double-difference(calc) bzgl. L1 in [cycles]
   wid1=(lss1(isat)-lss2(isat)-refs)/lambda1+(s2l1(t,rovindx)-s2l1(t,refs2indx)-s1l1(t,ksat+versatz)+s1l1(t,refs1indx));
   
   % amf-Wert f黵 den aktuellen Satelliten zur aktuellen Epoche f黵 L1 
   
   
   ddamfobs=ddamfobs+cos(2.*pi.*(wid1-(as2(isat-1,1).*dxcube1+as2(isat-1,2).*dycube1+as2(isat-1,3).*dzcube1))); 
   used_obs1=used_obs1+1;

end

% Dasselbe noch f黵 L2 
if mit_l2==1
if s2l2(t,rovindx)~=0 & s2l2(t,refs2indx)~=0 & s1l2(t,ksat+versatz)~=0 & s1l2(t,refs1indx)~=0
   
   % widerspruch = double-difference(obs) - double-difference(calc) bzgl. L2 in [cycles]
   wid2=(lss1(isat)-lss2(isat)-refs)/lambda2+(s2l2(t,rovindx)-s2l2(t,refs2indx)-s1l2(t,ksat+versatz)+s1l2(t,refs1indx));

   % amf-Wert f黵 den aktuellen Satelliten zur aktuellen Epoche f黵 L2
   ddamfobs=ddamfobs+cos(2.*pi.*(wid2-(as2(isat-1,1).*dxcube2+as2(isat-1,2).*dycube2+as2(isat-1,3).*dzcube2))); 
   used_obs2=used_obs2+1;
   
end      % Ende des IF-Blocks [s2l2(t,rovindx)~=0 & ...]
end      % Ende des IF-Blocks [mit_l2==1]



end      % Ende des IF-Blocks [ if isat>1 ]

else     % else zu [ if abs(s1ko_z)>cos((90-el_mask)/180*pi)& ..., also Elev.maske ]: 
   isat=isat-1;       % Wenn durch Elevationsmaske gefallen, zur點ksetzen von isat
end      % Ende des IF-Blocks [ if abs(s1ko_z)>cos((90-el_mask)/180*pi)& ..., also Elev.maske ]

end		% Ende des IF-Blocks [ if sum(rovmask)>0 ], also ob Sat auf beiden Stat. gemessen.

end      % Ende der Schleife ueber alle Satelliten
%------------------------------------------------------------------------------------------------

% Ergebnisse und Statistik nach der Epochenabarbeitung
% ----------------------------------------------------

if used_obs1>2        % wenn mindestens 3 double differences, dann:
   
   ddamfall=ddamfall+ddamfobs;                  % Aufakkumulation aller Epochen-Ergebnisse
   used_all=used_all+used_obs1+used_obs2;       % Anzahl aller aufakkumulierten Epochenergebnisse
   
% Statistik: RDOP(=sqrt(trace(inv(a'a)))), aufakkumulierte invertierte NGlMatrix

ns2=(as2'*as2)\eye(3);                  % invertierte Normalgleichungsmatrix
rdops2(t)=sqrt(trace(ns2));         % RDOP
n_accu=n_accu+ns2;                  % aufakkumulierte invert. Normalgleichungsmatrix

% !!!!!!!!!! Bestimmung einer epochenweisen L鰏ung (Kinematic) f鋘de hier statt !!!!!
% inclusive Korrektur wegen Antennenh鰄e (also Red. von Antennenphasenzentrum auf Punktmarke) ...

obs_is_used(t)=1;

else     % else zu [if used_obs1>2], d.h. wenn used_obs1<=2 :
   
obs_is_used(t)=0;
% ... nur diese eine Variable, um andere per Maske reduzieren zu k鰊nen (s.u.)

end      % Ende des IF-Blocks [if used_obs1>2]

%------------------------------------------------------------------------------------------------
waitbar(t/nobs1);   
end      % Ende der Schleife ueber alle Beobachtungen nobs1   ### ### ### ### ### ### ### ###
%------------------------------------------------------------------------------------------------
close(h);


% Maskierung: Wo gibt's > 3 Sat., d.h. > 2 doublediff. simultan ? (used_obs1>2)
% -------------------------------------------------------------------------------------------

dgpsmask=(obs_is_used~=0);

if sum(dgpsmask)==0
       errordlg(['Error in program AMFLAB.M during calculation.',...
             ' It seems that less than 4 satellites were observed or the program failed.' ...
             ' Check S2EPO in comparison to S2C1,S2L1 and S2L2, an check the program.'], ...   
       'GPSLab: Break -  1998 zeb');
       return;   
end

% --- Reduzierung der Felder auf die Elemente mit berechneten L鰏ungen -----------------
s1c1=s1c1(logical(dgpsmask),:);   % ...durch Herausstreichen der anderen Beobachtungen
s1epo=s1epo(logical(dgpsmask),:); % -"-
s1l1=s1l1(logical(dgpsmask),:);   
if mit_l2==1
   s1l2=s1l2(logical(dgpsmask),:);   
end

s2c1=s2c1(logical(dgpsmask),:);   % -"-
s2epo=s2epo(logical(dgpsmask),:); % -"-
s2l1=s2l1(logical(dgpsmask),:);   
if mit_l2==1
   s2l2=s2l2(logical(dgpsmask),:);   
end

s1epoms=s1epoms(logical(dgpsmask),:); % -"-
s2epoms=s2epoms(logical(dgpsmask),:); % -"-

rdops2=rdops2(logical(dgpsmask));


satmask1=satmask1(logical(dgpsmask),:);
satmask2=satmask2(logical(dgpsmask),:);


% Ergebnisse und Statistik nach allen Epochen:
% --------------------------------------------

% Normierung der Ambiguity Function aller 21x21x21 Testpositionen
ddamfnorm=ddamfall./used_all;

% Extrahieren der besten und der schlechtesten L鰏ung
% -------------------------------------------------

% Suchen des gr鲞ten Werts (Maximum)
[m1i,m1j]=find(ddamfnorm==max(max(max((ddamfnorm)))));
[m1k,m1l]=ind2sub([21,21],m1j);
best1st=ddamfnorm(m1i,m1k,m1l);

% L鰏chen des Maximums (nur bei Suche des 2.gr鲞ten Werts n鰐ig gewesen ...
%ddamfnorm(m1i,m1k,m1l)=-1;

% Suchen des kleinsten Werts (Minimum)
[m2i,m2j]=find(ddamfnorm==min(min(min((ddamfnorm)))));
[m2k,m2l]=ind2sub([21,21],m2j);
mist2nd=ddamfnorm(m2i,m2k,m2l);

% Wiederherstellen des Maximums
ddamfnorm(m1i,m1k,m1l)=best1st;

% Berechnen der Zuschl鋑e auf die N鋒erungskoordinaten
% ----------------------------------------------------

dxs2=(m1i-11)*raster;
dys2=(m1k-11)*raster;
dzs2=(m1l-11)*raster;

% Check auf ausreichendes Maximum: 50% erreicht ?
% -----------------------------------------------

if best1st<.50
       warndlg(['Attention: In program AMFLAB.M the data processing results in,',...
             ' a maximum of the ambiguity function with a value of < 50%.' ...
             ' Either the solution is out of search cube or the' ...
             ' Observation data are too noisy respectively insufficient. View the plots and the result file.'], ...   
       'GPSLab: Extremely small AMF-value !');
end

% Check auf Ratio
% ---------------
% wenn zweitbeste L鰏ung sehr nahe an der besten (unmittelbar daneben, 1 rasterweite),
% dann ratio<2 erlaubt, weil dann eindeutiges Max. !!!
% Denn: ratio > 1.4 / 2.0 (mit_l2==0) sinnvoll, wenn vom gleichen "Berg"?
%        ?????????????????????? ?????????????????????        ??????????????
% Evtl. noch Check einbauen, ob 1st und 2nd benachbart liegen und  nachfolgende Abfrage
% noch dementsprechend ab鋘dern ...
%
% NEIN! -> JETZT: Vergleich zwischen Max und Min: (evtl. noch zweite Abfrage anpassen, falls eindeutigere L鰏ungen)

% Nur noch abgewandelt in der Ergebnisdatei:
%if ((best1st-mist2nd)<1 & mit_l2==1) | ((best1st-mist2nd)<1 & mit_l2==0)
%       warndlg(['Achtung: Im Programm AMFLAB.M ergab sich beim Berechnen der Ergebnisse',...
%             ' ein Ratio von R<0.5 bei L1+L2, bzw. R<0.5 bei nur-L1. 躡erpr黤en Sie anhand der' ...
%             ' Plots und der Ergebnisdatei, ob der Suchbereich gro?genug gew鋒lt wurde.' ...
%             ' Ein Ratio von 0.5 bedeutet, da?Min und Max nur um die halbe cos-Amplitude' ...
%             ' auseinanderliegen. Dies weist auf ein flaches, wenig pr鋑nantes Maximum hin.' ...
%             ' Weiter durch Best鋞igung von OK ...'], ...   
%       'GPSLab: Warnung');
%end

% Nur beste L鰏ung:
% -----------------

% Varianz, sigma_0, Std.abw. der Unbekannten, Korrelationskoeffizienten
m_res=acos(best1st)/(2*pi)*lambda1;   % mittleres Residuum, hier behandelt wie ein zweifacher mittlerer Fehler
%                                            einer Beobachtung (DDMR = double diff. mean res.).
sigma_0=m_res/2*mean(rdops2);         % theoret. mittl. Positionsfehler (DDMR/2)*RDOP, entspricht (UERE*RDOP)
%                                            also wie mittlerer Fehler der Beob. multipliziert mit RDOP
sigma_0=sqrt(sigma_0^2+(raster/2)^2); % tats鋍hl. mittl. Positionsfehler evtl. mindestens so gro?wie halbe 
%                                            Rasterweite im Positionsraum; Fehlerfortpflanzungsgesetz
%                                            angewandt auf Rasterweite und theoret. sigma_0


ngl=n_accu/sum(dgpsmask);                 % mittlere inv. NGlMatrix, wobei sum(dgpsmask)=Anz. aller verwendeten Epochen
kxxs2=sigma_0^2*ngl;                      % Kovarianzmatrix der Unbekannten
sigxxs2=sqrt(diag(kxxs2));                % Standardabweichungen der Unbekannten
sxs2=sigxxs2(1);
sys2=sigxxs2(2);
szs2=sigxxs2(3);

korrxy=kxxs2(2,1)/(sxs2*sys2);
korrxz=kxxs2(3,1)/(sxs2*szs2);
korryz=kxxs2(3,2)/(sys2*szs2);


% Geoz.kart. Koord. und B, L, H f黵 S2: 


% Zuschl鋑e auf urspr黱gliche N鋒erungskoord. (Punkt-Punkt)
xxs2=xs2+dxs2;
yys2=ys2+dys2;
zzs2=zs2+dzs2;
[bbs2,lls2,hhs2]=xyz2blh(xxs2,yys2,zzs2);


% Berechnung der Koordinatendifferenzen und der Baseline-L鋘ge

dapp_h=app2_h-app1_h;  % H鰄endifferenz aus N鋒erungswerten
strdapph=num2str(dapp_h); % H鰄endifferenz aus N鋒erungswerten (String)
app_bsl=sqrt((xs2-xs1)^2+(ys2-ys1)^2+(zs2-zs1)^2); % Bsl-L鋘ge aus N鋒erungswerten
strappbsl=num2str(app_bsl); % Bsl-L鋘ge aus N鋒erungswerten (String)

m_delta_h=hhs2-app1_h;          % H鰄endifferenz

m_bsl=sqrt((xxs2-xs1).^2+(yys2-ys1).^2+(zzs2-zs1).^2); % Baseline-L鋘ge


% Beobachtungsstatistik
satnrs=s2epo(:,2:size(s2epo,2)); % ohne 1. Spalte (Epochen), nur Sat.nummern
satnrs=satnrs(:); % Vektor aus Matrix 
satnrs=satnrs(logical(satnrs>0)); % Nullen eliminieren (Verk黵zung des Vektors)

[sbbs2,slls2,shhs2,dums,dumA,dumz]=dxyz2neu(sxs2,sys2,szs2,bbs2,lls2);
sbbs2=abs(sbbs2);
slls2=abs(slls2);
shhs2=abs(shhs2);

rdmin=min(rdops2);
rdmax=max(rdops2);
rdmit=mean(rdops2);

[haufig2,n2]=hist(satnrs,1:32);
n2=n2(logical(haufig2>0));
haufig2=haufig2(logical(haufig2>0));


% Flag, da?Berechnung duchgef黨rt wurde und entsprechende Variable vorliegen

amf_is_calc=1;

% Ende  ==================================================================================
% (c) iapg 1999 Zebhauser