The high-precision construction of geomagnetic reference map is the basis of geomagnetic matching navigation. In order to address the problem of modelling the geomagnetic reference maps with different line spacing and different geomagnetic anomaly features, this paper combines the geomagnetic anomaly features with the principal component analysis method to integrate the geomagnetic standard deviation, mean geomagnetic field value, mean cumulative gradient value, kurtosis coefficient, and standard deviation of gradient, and the five geomagnetic feature parameters, as a measure of the change of geomagnetic anomaly features, and then thins out the measured 500 m-high-resolution data. The thinned line data with different line spacings are modelled using five modelling methods, namely, multifaceted function, minimum curvature, Kriging interpolation, local polynomials and radial basis function, and the areas of gentle geomagnetic variations and complex areas are selected for accuracy analysis. The experimental results show that the algorithms should be reasonably selected in the light of the line spacing, geomagnetic complexity and the modelling region, and the overall modelling accuracy of the multifaceted function is the highest; the multifaceted function has the highest modelling accuracy at different line spacings in the gentle region; the multifaceted function has the highest modelling accuracy in the complex region at line spacings from 2 to 5 km and at spacings greater than 8 km lines, and the minimum curvature modelling accuracy is the highest at line spacings from 5 to 7 km lines.