In computerized ionospheric tomography (CIT) with ground-based GNSS, the voxels without satellite-receiver ray traversing cannot be reconstructed directly. We present a CIT algorithm based on virtual reference stations (VRSs), called VRS–CIT, to decrease the number of unilluminated voxels and improve the precision of the estimated ionospheric electron density (IED). The VRSs are set at the nodes of grids with a 0.5° × 0.5° resolution in longitude and latitude. We generate the virtual observations with the observations from nearby six or three stations selected according to azimuths and distances. The generation utilizes multi-quadric surface fitting with six stations and triangular linear interpolation with three stations. With the virtual observations added, the IED distribution is reconstructed by the multiplicative algebraic reconstruction technique with the initial values obtained from IRI-2016. The performance of VRS–CIT is examined by using the data from 127 GNSS stations located in 24–46° N and 122–146° E to derive the IED every 30 min. The study focuses on April 29, 2014, with the adaptability of VRS–CIT analyzed by 12 days, evenly distributed around equinoxes and solstices of 2014. The accuracy of the virtual observation is about 1 TECU. Comparing to that derived from CIT with only real observations, the unsolvability of VRS–CIT declined by 4–12% for the whole region, and for the main area, the improvement can be up to 70%. Taking two IED profiles from radio occultation as reference measurements, the mean absolute error (MAE) of IED by VRS–CIT decreases by 6.88% and 8.43%, respectively. Comparing with slant total electron content (STEC) extracted from five additional GNSS stations, the MAE and the root mean square error of the estimated STEC can be reduced up to 17.24% and 33.81%, respectively.