GPS time series analysis (2) matlab language analysis
1. Simple GPS time series
load E:\RMTE.txt yuanshi_dataWeidu=RMTE(:,4)';%Assigning values to original coordinates yuanshi_dataJingdu=RMTE(:,3)'; yuanshi_dataGaochen=RMTE(:,5)'; yuanshi_data=[yuanshi_dataWeidu;yuanshi_dataJingdu;yuanshi_dataGaochen]; [col,raw]=size(yuanshi_data); llh_Z_xyz=llh2xyz(yuanshi_data);%Longitude-Latitude Elevation Conversion xyz Planar coordinates, output format: 3 lines n column %After the first conversion xyz Coordinates minus first observed xyz Coordinates, composition xyz Variable Matrix llh_Z_xyz_X=llh_Z_xyz(1,:)-mean(llh_Z_xyz(1,:)); llh_Z_xyz_Y=llh_Z_xyz(2,:)-mean(llh_Z_xyz(2,:)); llh_Z_xyz_Z=llh_Z_xyz(3,:)-mean(llh_Z_xyz(3,:)); llh_Z_xyz_XYZ=[llh_Z_xyz_X;llh_Z_xyz_Y;llh_Z_xyz_Z]; %This is the second step, converting the first time to a good one xyz Convert to NEU_xyz,Arrange the matrix of changes into 3*n Column Form,Arrange in dot order llh_Z_xyz3=reshape(llh_Z_xyz_XYZ,col*raw,1); %Take latitude and longitude of the first observation as llh Parameter, will xyz Coordinate conversion to NEU_xyz Variable Matrix latitude=yuanshi_data(1,1); longitude=yuanshi_data(2,1); altitude=yuanshi_data(3,1); neu=xyz2enu([latitude longitude altitude],llh_Z_xyz3);%take xyz Coordinate conversion NEU_xyz coordinate NEU_xyz=reshape(neu,col,raw);%NEU_xyz After conversion, arrange in 3 lines n column
RMTE_NEU_xyz_VariableQuantity; load E:\RMTE.txt E=RMTE(:,1);%vector a All values are 1 Time=RMTE(:,2);%station GPS Time Series Observation Time NEU_xyz_WE=NEU_xyz(1,:)*1000;%East-West Direction NEU_xyz_NS=NEU_xyz(2,:)*1000;%North-South Direction NEU_xyz_GC=NEU_xyz(3,:)*1000;%Elevation Direction %Finding Origin Coordinate Residual in East-West Direction subplot(3,1,1); p=polyfit(Time,NEU_xyz_WE',1); z=polyval(p,Time); plot(Time,NEU_xyz_WE,'.black',Time,z,'r') ylabel('Latitude(mm)'); xlabel('Time'); hTit(1)=title('RMTE Original sequence diagram of east-west direction of station');
2 Noise Item Residual Distribution Item Graph
RMTE_NEU_xyz_VariableQuantity; load E:\RMTE.txt E=RMTE(:,1);%vector a All values are 1 Time=RMTE(:,2);%Wuhan Station GPS Time Series Observation Time NEU_xyz_WE=NEU_xyz(1,:);%East-West Direction NEU_xyz_NS=NEU_xyz(2,:);%North-South Direction NEU_xyz_GC=NEU_xyz(3,:);%Elevation Direction for i=1:1:length(Time) sin_2piTime(i)=sin(2*pi*Time(i)); end for j=1:1:length(Time) cos_2piTime(j)=cos(2*pi*Time(j)); end for m=1:1:length(Time) sin_4piTime(m)=sin(4*pi*Time(m)); end for n=1:1:length(Time) cos_4piTime(n)=cos(4*pi*Time(n)); end B=[E Time sin_2piTime' cos_2piTime' sin_4piTime' cos_4piTime'];%B matrix abcdef_NEU_xyz_NS=(inv(B'*B))*(B'*NEU_xyz_NS');%Calculating North-South Direction Parameters abcdef Value of abcdef_NEU_xyz_WE=(inv(B'*B))*(B'*NEU_xyz_WE');%Finding East-West Direction Parameters abcdef Value of abcdef_NEU_xyz_GC=(inv(B'*B))*(B'*NEU_xyz_GC');%Calculating Elevation Direction Parameters abcdef Value of %Correcting noise coordinate residuals in east-west direction subplot(3,1,1); V_NEU_xyz_WE=(NEU_xyz_WE'- B*abcdef_NEU_xyz_WE)*1000; plot(Time,V_NEU_xyz_WE,'.black') ylabel('Latitude(mm)'); xlabel('Time'); hTit(1)=title;
3 GPS Time Series Analysis Model
GPS time series contains a lot of useful information, which needs to be based on the principles and methods of time series analysis. Currently, the main models introduced are power spectrum analysis, maximum likelihood estimation, wavelet analysis and space-time filtering.
3.1 Power Spectrum Analysis
Power spectrum analysis is an effective method to study time series from frequency domain by decomposing time series into amplitudes in several different frequency ranges using Fourier transform method.By studying power spectrum analysis, the noise type of time series can be obtained, and the period term and amplitude of time series can also be obtained. It is a better method to analyze the time series of GPS station position.
The noise time series after deducting the fitted values from the original GPS time series can be used to qualitatively analyze the noise type using the spectrum.Spectrum analysis is a method of analyzing signals in frequency domain.The power spectrum of noise time series in the spectral domain can be expressed as power law
function[a,Pf,f]=Power_spectrum_analysis(t,st) %This is a function using the FFT function to caculate a signal's Fourier Translation %Input is the time and the signal vectors ,the length of time must greater than 2 %Output is the frequency and the signal spectrum t1=t(1);t2=t(end); N=length(st); Fs=N/(t2-t1); window=boxcar(N); %Rectangular window [Pf,f]=periodogram(st,window,N,Fs); %direct method %InP(f) = InP0 - aInf %Solving spectral index k and Po n=length(Pf); m=length(f); if n==m for i=1:n-1 B(i,1)=1; B(i,2)=-log(f(i+1)); L(i)=log(Pf(i+1)); end a=inv(B'*B)*B'*L'; end
%Finding the Power Spectrum of Noise Residual Terms N=length(Time); T=0; for i=1;N T=T+Time(i+1)+Time(i); %Find Period T end anv_T=T/N; FS=1/anv_T; V_NEU_xyz_WE=NEU_xyz_WE'- B*abcdef_NEU_xyz_WE;%Finding the power spectrum of noise residuals in east-west direction [a1,Pf1,f1]=Power_spectrum_analysis(Time,V_NEU_xyz_WE)
figure PX_GC=100000*Pf3; F_GC=0.28*f3; % subplot(3,1,3); loglog(F_GC,PX_GC,'black') axis([0.08,100,0.00009,10]); ylabel('power/(mm^2/cpy)'); xlabel('frequency/cpy)'); legend('Vertical');
3.2 Maximum Likelihood Estimation (MLE)
The maximum likelihood method is a non-linear least squares method to find the model parameters that are closest to the time series.By adjusting the covariance matrix to maximize the likelihood function, the noise model closest to the time series can be obtained, so that the noise amplitude in the time series can be estimated by the covariance matrix.Compared with spectral analysis, maximum likelihood estimation can quantitatively calculate the noise size.
Assume that one-way position time series of GPS stations is:
3.3 Wavelet Analysis
Wavelet analysis is a time-frequency analysis method of signal. It has the characteristics of multi-resolution analysis, and has the ability to characterize signal characteristics in both time and frequency domains. The size of time window can change with the change of frequency window. With the help of wavelet method, noise with different characteristics can be separated and estimated to some extent.Wavelet multiresolution analysis can decompose the period term of GPS site coordinate time series to extract anniversary and half anniversary items. Therefore, this method is called anniversary mathematical filter and half anniversary mathematical filter.
3.4 Elevation Direction AR, ARMA Model
ARMA series models are an effective means to study time series. ARMA models study system characteristics when only system output is known from the perspective of system identification.The ARMA series is divided into AR, MA and ARMA models.AR model and ARMA model are mainly used to model elevation time series.The ARMA model is:
4 GPS coordinate time series space-time filtering method
There are some common errors in GPS coordinate time series, which are often referred to as common mode error (CME). In order to improve the accuracy of GPS time series, some research needs to be done.There are many methods to extract common mode errors in GPS site coordinate time series. At present, principal component analysis (PCA) and (KLE) are mainly used for analysis.
4.1 PCA Method
Principal Component Analysis is an effective tool for modern data analysis.PCA recombines the original correlated observation data into a set of unrelated vectors, eliminates the correlated dimensions in the original independent variable and transforms them into low-dimensional vector space by linear transformation while minimizing the loss of data information. The transformed low-dimensional space has orthogonal principal components, which synthesize the maximum amount of information of the original data, so the data is dividedThe results of the analysis will not have much impact and are more spatially valid.It avoids the ill-conditioned problem of the normal equation and can obtain accurate estimates of the parameters, and reveals some rules and structural features hidden behind the data.However, if the local noise is high, the common mode error can not be accurately distinguished by the PCA.
When the cumulative contribution rate reaches this threshold, the first P values correspond to the main mode components.The principal component is actually the projection of the residual vector of a spatial site in the direction of the corresponding eigenvector. A set of principal components reflects the coordinates of the residual vector at a corresponding time on the corresponding spatial axis expressed in the direction of the eigenvector.The first P principal component corresponds to a large eigenvalue, that is, it contains most of the contribution values to the time series variance of the residuals. It has the most information in the whole region and often reflects the common trend of the entire network.The common mode error is determined by the first P principal component.
%Solving principal components, common mode errors between stations CME %list New_Time_VGC That is, the elevation direction correlation coefficient X_VGC X_VGC=[BMCL_New_Time_VGC BRN2_New_Time_VGC CHLM_New_Time_VGC DLPA_New_Time_VGC DNGD_New_Time_VGC GRHI_New_Time_VGC JMLA_New_Time_VGC KKN4_New_Time_VGC LHAZ_New_Time_VGC NPGJ_New_Time_VGC ODRE_New_Time_VGC RMTE_New_Time_VGC SMKT_New_Time_VGC TPLJ_New_Time_VGC XZDX_New_Time_VGC XZGE_New_Time_VGC ]; %Data related to elevation direction X_VGC Standardized to X_VGC_BZ X_VGC_COV = cov(X_VGC); %Finding covariance matrix [X_VGC_Xv , X_VGC_Xd] = eig(X_VGC_COV); %Eigenvalues and eigenvectors corresponding to covariance matrices X_VGC_GXL = X_VGC_Xd / sum(X_VGC_Xd); X_VGC_GXL =( fliplr(X_VGC_GXL'))'; %Common mode error in elevation direction CME Extraction, to select a variable with a contribution rate and an eigenvalue greater than 1 %according to X_VGC_tent A conclusion greater than 1 indicates that the first four variables should be selected x1,x2,x3,x4. %The method of extracting principal components is to multiply the eigenvector matrix by the normalized variable,All variables here are standardized. % for j = 1:length(X_VGC_Xv) % for m = 1:length(X_VGC) % Y1_VGC(m,1)=0; % Y1_VGC(m,1) = Y1_VGC(m,1)+ X_VGC_Xv(j,1)* X_VGC(m,j); % end % end Y1_VGC1(:,1) = X_VGC_Xv(1,1)* X_VGC(:,1); Y1_VGC2(:,1) = X_VGC_Xv(2,1)* X_VGC(:,2); Y1_VGC3(:,1) = X_VGC_Xv(3,1)* X_VGC(:,3); Y1_VGC4(:,1) = X_VGC_Xv(4,1)* X_VGC(:,4); Y1_VGC5(:,1) = X_VGC_Xv(5,1)* X_VGC(:,5); Y1_VGC6(:,1) = X_VGC_Xv(6,1)* X_VGC(:,6); Y1_VGC7(:,1) = X_VGC_Xv(7,1)* X_VGC(:,7); Y1_VGC8(:,1) = X_VGC_Xv(8,1)* X_VGC(:,8); Y1_VGC9(:,1) = X_VGC_Xv(9,1)* X_VGC(:,9); Y1_VGC10(:,1) = X_VGC_Xv(10,1)* X_VGC(:,10); Y1_VGC11(:,1) = X_VGC_Xv(11,1)* X_VGC(:,11); Y1_VGC12(:,1) = X_VGC_Xv(12,1)* X_VGC(:,12); Y1_VGC13(:,1) = X_VGC_Xv(13,1)* X_VGC(:,13); Y1_VGC14(:,1) = X_VGC_Xv(14,1)* X_VGC(:,14); Y1_VGC15(:,1) = X_VGC_Xv(15,1)* X_VGC(:,15); Y1_VGC16(:,1) = X_VGC_Xv(16,1)* X_VGC(:,16); Y1_VGC(:,1) = Y1_VGC1(:,1)+Y1_VGC2(:,1)+Y1_VGC3(:,1)+Y1_VGC4(:,1)+Y1_VGC5(:,1)+Y1_VGC6(:,1)+Y1_VGC7(:,1)+Y1_VGC8(:,1)+Y1_VGC9(:,1)+Y1_VGC10(:,1)+Y1_VGC11(:,1)+Y1_VGC12(:,1)+Y1_VGC13(:,1)+Y1_VGC14(:,1)+Y1_VGC15(:,1)+Y1_VGC16(:,1);
4.2 KLE Method
KLE is similar to PCA in that it is also an orthogonal decomposition-based method. Unlike PCA, KLE standardizes the covariance matrices of residual time series matrices by using correlation matrices instead of covariance matrices to compute orthogonal vector bases.KLE is insensitive to local noise and can extract the required signal from time series with a lot of noise, which can better make up for the shortcomings of PCA.
PCA can extract the spatial characteristics of time series well. KLE method can extract continuous signals from time series with strong local noise. PCA/KLE method combined with PCA and KLE is used to filter GPS coordinate time series to remove common mode error in GPS time series.Assuming there are n stations in the area, the time series of each station are arranged in columns, which can be expressed as a matrix of m x n:
%KLE Method for finding common mode error %Solving principal components, common mode errors between stations CME %list New_Time_VGC That is, the elevation direction correlation coefficient X_VGC X_VGC=[BMCL_New_Time_VGC BRN2_New_Time_VGC CHLM_New_Time_VGC DLPA_New_Time_VGC DNGD_New_Time_VGC GRHI_New_Time_VGC JMLA_New_Time_VGC KKN4_New_Time_VGC LHAZ_New_Time_VGC NPGJ_New_Time_VGC ODRE_New_Time_VGC RMTE_New_Time_VGC SMKT_New_Time_VGC TPLJ_New_Time_VGC XZDX_New_Time_VGC XZGE_New_Time_VGC ]; %Do principal component analysis and solve X_VGC_pc Orthogonal unitary eigenvector matrix, X_VGC_la Score matrix, X_VGC_tent characteristic value X_VGC_BZ = zscore(X_VGC); R_X_VGC = corrcoef(X_VGC); %Correlation coefficient matrix for elevation direction [R_X_VGC_LG_v,R_X_VGC_LG_d] = eig (R_X_VGC); %Eigenvalues and eigenvectors corresponding to correlation coefficient matrices R_X_VGC_LG_w = R_X_VGC_LG_d / sum(R_X_VGC_LG_d);%Contribution ratio of each principal component R_X_VGC_LG_w = (fliplr(R_X_VGC_LG_w'))'; Y1_VGC1(:,1) = R_X_VGC_LG_v(1,1)* X_VGC(:,1); Y1_VGC2(:,1) = R_X_VGC_LG_v(2,1)* X_VGC(:,2); Y1_VGC3(:,1) = R_X_VGC_LG_v(3,1)* X_VGC(:,3); Y1_VGC4(:,1) = R_X_VGC_LG_v(4,1)* X_VGC(:,4); Y1_VGC5(:,1) = R_X_VGC_LG_v(5,1)* X_VGC(:,5); Y1_VGC6(:,1) = R_X_VGC_LG_v(6,1)* X_VGC(:,6); Y1_VGC7(:,1) = R_X_VGC_LG_v(7,1)* X_VGC(:,7); Y1_VGC8(:,1) = R_X_VGC_LG_v(8,1)* X_VGC(:,8); Y1_VGC9(:,1) = R_X_VGC_LG_v(9,1)* X_VGC(:,9); Y1_VGC10(:,1) = R_X_VGC_LG_v(10,1)* X_VGC(:,10); Y1_VGC11(:,1) = R_X_VGC_LG_v(11,1)* X_VGC(:,11); Y1_VGC12(:,1) = R_X_VGC_LG_v(12,1)* X_VGC(:,12); Y1_VGC13(:,1) = R_X_VGC_LG_v(13,1)* X_VGC(:,13); Y1_VGC14(:,1) = R_X_VGC_LG_v(14,1)* X_VGC(:,14); Y1_VGC15(:,1) = R_X_VGC_LG_v(15,1)* X_VGC(:,15); Y1_VGC16(:,1) = R_X_VGC_LG_v(16,1)* X_VGC(:,16); Y1_VGC(:,1) = Y1_VGC1(:,1)+Y1_VGC2(:,1)+Y1_VGC3(:,1)+Y1_VGC4(:,1)+Y1_VGC5(:,1)+Y1_VGC6(:,1)+Y1_VGC7(:,1)+Y1_VGC8(:,1)+Y1_VGC9(:,1)+Y1_VGC10(:,1)+Y1_VGC11(:,1)+Y1_VGC12(:,1)+Y1_VGC13(:,1)+Y1_VGC14(:,1)+Y1_VGC15(:,1)+Y1_VGC16(:,1);
Totem Xu Wei 23 original articles published, 5 praised, 945 visits Private letter follow