Method for the computation of horizontal and vertical protection levels, HPL VPL, that bound up to a given confidence level 1-α the horizontal and vertical components, respectively, of the position estimation error δ of a least squares-based GNSS navigation solution, wherein the HPL and VPL are computed as follows:HPL=k·∥r∥·aH, VPL=k·∥r∥aV,where:∥r∥ is the Euclidean norm of the least squares estimation residuals vector r;aH is defined as:aH=hEE+hNN2+(hEE-hNN2)2+(hEN+hNE2)2,andaV is defined as: aV=√{square root over (hUU)}, wherehEE, hNN, hEN, hNE, hEU, hUE, hUU, hNU, hUN are the spatial components of the dilution of precision matrix DoP (HT·H)−1 of the least squares estimation expressed in the local horizontal coordinate system of the estimated position:(HT·H)Space-1=[hEEhENhEUhNEhNNhNUhUEhUNhUU],H being the n×m observation matrix of the least squares-based navigation solution;and wherek is an isotropic confidence ratio computed by numerically solving the following expression Eq. 3 which links k with n, m and a, n and m being the number of observations and of estimated parameters, respectively, of the least squares estimation:Γ(n2)·∫z∈ℜm,z2≤k21+k2(1-z2)n-m-22·z=Γ(n-m2)·πm2·(1-α)[Eq.3]where Γ is Euler's Gamma function; and m denotes a standard m-dimensional real vector space.