【Objective】Improving the spatial prediction accuracy of soil attributes is of great significance for achieving accurate fertilization of farmland and protecting the ecological environment. 【Method】Soil organic matter (SOM) data was collected from 1773 samples from soil surface layer (0-20cm) of cultivated land in Luanping County, Hebei Province. The optimal environmental variables were screened through a stepwise regression analysis method. Multiple linear regression (MLR), ordinary kriging (OK), random forest (RF), Bayesian regularized neural network (BRNNBP), and the corresponding three integrated models combined with a geostatistical model (MLR-OK, RF-OK and BRNNBP-OK) were utilized to predict SOM content via the training set including 1426 sampling points. Also, the prediction accuracy of each method was compared with 347 sampling points of the testing set. Autocorrelation analysis was carried out based on the residual of the integrated model to evaluate the fitting effect of the model. 【Result】Results showed that the range of SOM content in the study area was 8.62～35.64 g·kg^-1, and the coefficient of variation was 20.26%, which showed a moderate spatial variation. High concentrations of SOM were mainly distributed in the northeast and southeast areas with higher altitudes, while relative low concentrations of SOM were mostly observed in the southwest and central valleys of the study area. Elevation, slope and temperature selected by stepwise regression were closely related to SOM content (P<0.001). The lowest average absolute error and the root mean square error of the BRNNBP-OK model were 2.162 g·kg^-1and 2.801 g·kg^-1, respectively. Compared with the OK, MLR, RF, BRNNBP, MLR-OK and RF-OK models, the goodness of fit of the BRNNBP-OK model increased by 1.84%～43.72%, making it an optimal model for SOM spatial prediction. Compared with the single model, the nugget coefficient of the integrated model residual was greater than 0.75, and the Moran's I was less than 0 and numerically closer to 0, indicating that the spatial autocorrelation of the integrated model residual was weakened and the residual presented a more discrete spatial distribution. At the same time, the accuracy of all models was significantly correlated with Moran's index of model residuals. 【Conclusion】In this study, the integrated model fitted more trends and the spatial aggregation of model residuals decreased and even tended to be discrete. Thus, the overall prediction accuracy of the integrated models was improved.