To solve the problems of high calculation amounts and low classification accuracy of Gaussian process regression in hyperspectral image classification, a spectral-spatial joint classification algorithm for hyperspectral images based on improved Gaussian process regression was proposed. A subset of samples was selected using maximum variance as the index to narrow the calculation range of the Gaussian process regression parameter solution, and a square root matrix decomposition method was introduced to predict the model results for incoming added samples, all of which effectively improve the efficiency of calculation. A spatial-spectral correlation distance of neighborhood pixels was redefined in the pixel neighbor space based on spatial-spectral feature information. In addition, a space-spectrum correlation distance integrated with spatial neighbor information was used as the weight to measure the similarity of neighborhood pixels. These increase the probability that similar features would be classified as neighbors, thus improving the accuracy of feature classification. Simulation experiments were conducted on two sets of hyperspectral datasets from Indian Pines and Pavia University. Experimental results show that, compared with other similar algorithms, the proposed algorithm improves overall classification accuracy, average classification accuracy, and the Kappa coefficient by at least 2.3%, 1.4%, and 1.07%, respectively. Compared with the model algorithm prior to enhancements, the improved algorithm not only achieves higher overall classification accuracy but also considerably reduces the running time.