Ill‐posed seismic inverse problems are often solved using Tikhonov‐type regularization, that is, incorporation of damping and smoothing to obtain stable results. This typically results in overly smooth models, poor amplitude resolution, and a difficult choice between plausible models. Recognizing that the average of parameters can be better constrained than individual parameters, we propose a seismic tomography method that stabilizes the inverse problem by projecting the original high‐dimension model space onto random low‐dimension subspaces and then infers the high‐dimensional solution from combinations of such subspaces. The subspaces are formed by functions constant in Poisson Voronoi cells, which can be viewed as the mean of parameters near a certain location. The low‐dimensional problems are better constrained, and image reconstruction of the subspaces does not require explicit regularization. Moreover, the low‐dimension subspaces can be recovered by subsets of the whole dataset, which increases efficiency and offers opportunities to mitigate uneven sampling of the model space. The final (high‐dimension) model is then obtained from the low‐dimension images in different subspaces either by solving another normal equation or simply by averaging the low‐dimension images. Importantly, model uncertainty can be obtained directly from images in different subspaces. Synthetic tests show that our method outperforms conventional methods both in terms of geometry and amplitude recovery. The application to southern California plate boundary region also validates the robustness of our method by imaging geologically consistent features as well as strong along‐strike variations of San Jacinto fault that are not clearly seen using conventional methods.