Measuring the 3D displacement fields provide essential information regarding the Earth crust interaction and the mantle rheology. The interferometric synthetic aperture radar (InSAR) has an appropriate capability in revealing the displacements of the Earth’s crust. Although, it measures the real 3D displacements in the line of sight (LOS) direction. The 3D displacement vectors can be retrieved through multiple InSAR measurements acquired from at least three independent imaging geometries in a theoretical manner. However, this is a physically ill-posed inverse problem and consequently, the retrieving process of the components in 3D displacements become sensitive to observation errors, especially in the northern component due to the near-polar orbiting of SAR missions. Combining different datasets regarding this issue requires proper treatment of the weight of observations, which otherwise will have a negative effect on both the precision and accuracy of the estimated 3D displacement field. In retrieving the 3D displacement fields through InSAR technique, we deal with two major issues, integration of inhomogeneous precision of observations and instability of the estimation problem. These facts constitute the motivations to address the Tikhonov regularization (TR) and least squares variance component estimation (LS-VCE). In this article, to overcome these drawbacks, the regularized least squares variance component estimation (RLS-VCE) is proposed for retrieving the 3D displacement vectors. Usually, the number of InSAR observations in relation to the three unknowns of 3D displacements for each pixel is not enough to apply VCE. Therefore, observations of some neighborhood cells are taken into account to increase the redundancy of stochastic model. In this context, a moving frame including a window of 3 × 3 pixels is considered to increase the number of observations and consequently, the degree of freedom of stochastic model. To assess the efficiency of the proposed method, the RADAR dataset of the Envisat and ALOS missions for the 17 June 2007 eruption of Kilauea volcano on Hawaiian island are applied. To validate the results of the proposed method, co-event displacement vectors of 19 GNSS stations around the Kilauea volcano are used. Furthermore, the 3D displacements of GNSS stations are applied for detrending the displacements of InSAR from systematic or random disturbing effects (e.g. orbit errors, curvature and topography of the Earth, atmosphere, etc.) through fitting a two variates linear or quadratic polynomial. Comparing the co-event retrieved 3D displacement vectors through RLS-VCE method and GNSS measurements indicates that the componential RMSE of northern displacements decreases drastically to 2.2 cm from 11.7 cm (for range displacements and primary weights). This is approximately equivalent to 80% improvement in the accuracy of estimating the northern component of displacement. The overall RMSE of retrieving 3D displacement vectors decrease from 7.8 cm to 2.6 cm, which is equal to 66% improvement. Achieving to this overall accuracy and for northern component is of major interest for all disciplines of geoscience dealing with 3D surface deformation analysis. Results indicate that retrieving the 3D displacement vectors through applying the RLS-VCE method has a meaningful improvement on the precision and accuracy of the results, the northern-southern component in special.