Learning-based distribution system state estimation (DSSE) methods typically depend on sufficient fully labeled data to construct mapping functions. However, collecting historical labels (state variables) can be challenging and costly in practice, resulting in performance degradation for these methods. To fully leverage low-cost unlabeled historical measurement data, this article proposes a Gaussian process (GP)-regularized semi-supervised learning method for DSSE models, aiming at achieving feasible estimation precision using minimal state variable labels while also providing valuable interval estimation of state variables. Firstly, a structure-informed graphical encoder is established to generate appropriate node embeddings. A tailored GP-regularized learning method is then developed to model the intermediate latent space using these embeddings. It constructs the unlabeled embeddings by a weighted combination of labeled space vectors, with the weights determined by a kernel function, thereby forming additional supervision and regularizing the learning process of the network in the latent space. The regularized embeddings are then fed into a decoder to yield estimation outcomes. This procedure enables the proposed method to model intrinsic correlations across measurement data, as well as capture essential patterns related to DSSE even using extremely limited state labels. The trained DSSE models can thus adapt to the domain of new measurements. Lastly, the decoder outcomes and related latent embeddings are processed through a composite GP kernel to further derive the interval estimation of state variables, enabling uncertainty quantification. Experimental results demonstrate the effectiveness of the proposed method in handling extremely limited historical state labels and accurately quantifying the uncertainty of state variables.