The acquisition of a 3D displacement field can help to understand the crustal deformation pattern of seismogenic faults and deepen the understanding of the earthquake nucleation. The data for 3D displacement field extraction are usually from GPS/interferometric synthetic aperture radar (InSAR) observations, and the direct solution method is usually adopted. We proposed an iterative least squares for virtual observation (VOILS) based on the maximum a posteriori estimation criterion of Bayesian theorem to correct the errors caused by the GPS displacement interpolation process. Firstly, in the simulation examples, both uniform and non-uniform sampling schemes for GPS observation were used to extract 3D displacement. On the basis of the experimental results of the reverse fault, the normal fault with a strike-slip component, and the strike-slip fault with a reverse component, we found that the VOILS method is better than the direct solution method in both horizontal and vertical directions. When a uniform sampling scheme was adopted, the percentages of improvement for the reverse fault ranged from 3% to 9% and up to 70%, for the normal fault with a strike-slip component ranging from 4% to 8% and up to 68%, and for the strike-slip fault with a reverse component ranging from 1% to 8% and up to 22%. After this, the VOILS method was applied to extract the 3D displacement field of the 2008 Mw 7.9 Wenchuan earthquake. In the East–West (E) direction, the maximum displacement of the hanging wall was 1.69 m and 2.15 m in the footwall. As for the North–South (N) direction, the maximum displacement of the hanging wall was 0.82 m for the southwestern, 0.95 m for the northeastern, while that of the footwall was 0.77 m. In the vertical (U) direction, the maximum uplift was 1.19 m and 0.95 m for the subsidence, which was significantly different from the direct solution method. Finally, the derived vertical displacements were also compared with the ruptures from field investigations, indicating that the VOILS method can reduce the impact of the interpolated errors on parameter estimations to some extent. The simulation experiments and the case study of the 3D displacement field for the 2008 Wenchuan earthquake suggest that the VOILS method proposed in this study is feasible and effective, and the degree of improvement in the vertical direction is particularly significant.