Three methods of GNSS velocity calculation

Position difference

Position difference is the simplest and crudest method of velocity estimation , We only need to know the position estimation of the two epochs $\overrightarrow{r_u}$ And epoch interval $t_k-t_{k-1}$ The speed of the receiver can be estimated by the following formula .
$$\overrightarrow{v_u}=\frac{\overrightarrow{r_u}(t_k)-\overrightarrow{r_u}(t_{k-1})}{t_k-t_{k-1}}$$
Here's the picture , We go through $P_2$ Point and $P_1$ The position of the point can estimate the velocity , But the disadvantages of this method are also obvious ,

• We need to know Two Location information of epochs
• The result of the calculation is not the speed of the facts , But from $P_1$ To $P_2$ Average speed of , So this speed has a certain lag
• Therefore, this method should consider the dynamics of the carrier and the frequency of position output , Whether it meets the application scenario
• There's another problem , The velocity estimation and position estimation obtained by this method are not independent and uncorrelated
  

Velocity estimation based on pseudo range rate

Is there a better way ？ The Doppler shift received by the receiver can be expressed as the following formula .
$$Doppler=f_r-f_t=-\frac{f_t}{c}[(v_s-v_u)\bullet \frac{r_s-r_u}{||r_s-r_u||}]$$
among $\bullet$ For point by point , We multiply both sides of the above formula by the wavelength of the carrier $\lambda$, So the formula above becomes ,
$$\dot{\rho }=(v_s-v_u)\bullet E$$
among $E$ Is the line of sight vector $\frac{r_s-r_u}{||r_s-r_u||}$, Consider the clock drift of the receiver $\delta \dot{t_u}$ And measurement noise $\xi
$, The formula above can be completely written as ,
$$\dot{\rho }=(v_s-v_u)\bullet E+c\delta \dot{t_u}+\xi$$

Here we can estimate the receiver speed , Let's look at it briefly rtklib Source code , A key behavior vs[j]=rs[j+3+i*6]-x[j];, From this we can roughly see rtklib This is the method used .

//  The following code is located in pntpos.c  Medium resdop  function 
 /* LOS (line-of-sight) vector in ECEF */
        cosel=cos(azel[1+i*2]);
        a[0]=sin(azel[i*2])*cosel;
        a[1]=cos(azel[i*2])*cosel;
        a[2]=sin(azel[1+i*2]);
        matmul("TN",3,1,3,1.0,E,a,0.0,e);
        
        /* satellite velocity relative to receiver in ECEF */
        for (j=0;j<3;j++) {
    
            vs[j]=rs[j+3+i*6]-x[j];
        }

There is also a problem with this method , Use this line of code vs[j]=rs[j+3+i*6]-x[j]; When calculating the pseudo range rate , Among them x For the unknown , So this method must be the same as solving the position , It needs to be iterated . We can look at the function of solving speed and verify .

//  The following code comes from pntpos.c estvel  function 
    for (i=0;i<MAXITR;i++) {
    
        
        /* range rate residuals (m/s) */
        if ((nv=resdop(obs,n,rs,dts,nav,sol->rr,x,azel,vsat,err,v,H))<4) {
    
            break;
        }

This for Loop is used to iteratively solve the speed ,MAXITR Is the maximum number of iterations .

Velocity estimation based on Doppler frequency shift

Can we not iterate ？ The answer is yes .
First , The Doppler shift of a satellite can be written in the following form ,
$$d_s=\dot{{\rho }_s}-v_s\bullet E$$
We bring it into the pseudo range rate formula ,
$$d_s=-E\bullet v_u+c\delta \dot{t_u}+\xi$$
We have obtained the position of the satellite and receiver through position calculation , therefore $E$ Known quantity , and $d_s$ To measure , So let's put the observation equations of each satellite together , Then we get a linear constant equation , It can be solved directly by least square or Kalman filter .
$$\begin{gathered} d_{s1}=-E_1\bullet v_u+c\delta \dot{t_u}+\xi \\ {d}_{s2}=-E_2\bullet v_u+c\delta \dot{t_u}+\xi \\ ... \\ {d}_{sm}=-E_m\bullet v_u+c\delta \dot{t_u}+\xi \end{gathered}$$

Use rtklib Solve the velocity

First , We download rtklib The executable of ,

function rtkpost And input ephemeris file and observation file , Click on Execute

Okay , Solution complete ！