GNSS and InSAR data are integrated to extract the 3-D surface deformations, which are of great significance for studying geological hazards. In this study, two major problems are focused on integration. For one thing, we propose an Iterated Almost Unbiased Estimation (IAUE) method to estimate the variance components of GNSS and InSAR for the case where the estimation of variance components of multi-source data by traditional variance component estimation (VCE) methods may be negative and inaccurate. For another, considering that heterogeneous data errors may lead to unstable three-dimensional solutions, we propose adding the Laplacian smoothness constraint (LSC) to the function model, which can smooth the solutions by minimizing the second derivative of the displacements. These two methods are abbreviated as IAUE-LSC. In the simulation experiment, the performance of traditional Helmert variance component estimation (HVCE) is first compared with IAUE. IAUE can not only converge more quickly but also avoid negative variances. Furthermore, we find that the excessively large relative error ratio between GNSS and InSAR is an essential factor leading to the instability of the three-dimensional solutions. IAUE-LSC method is immune to this instability and can obtain more stable results. In addition, the 2018 Hawaii case demonstrates that IAUE achieves improvements of 2.58 cm, 2.77 cm, and 7.69 cm in the east, north, and up directions relative to the traditional weighted least squares (WLS) method. While the combined IAUE-LSC achieve improvements of 2.29 cm, 0.32 cm, and 1.68 cm compared to the IAUE alone.