With the development of industrialization and urbanization, heavy metal contamination in agricultural soils tends to accumulate rapidly and harm human health. Visible and near-infrared (Vis-NIR) spectroscopy provides the feasibility of fast monitoring of the variation of heavy metals. This study explored the potential of fractional-order derivative (FOD), the optimal band combination algorithm and different mathematical models in estimating soil heavy metals with Vis-NIR spectroscopy. A total of 80 soil samples were collected from an agriculture area in Suzi river basin, Liaoning Province, China. The spectra for mercury (Hg), chromium (Cr), and copper (Cu) of the samples were obtained in the laboratory. For spectral preprocessing, FODs were allowed to vary from 0 to 2 with an increment of 0.2 at each step, and the optimal band combination algorithm was applied to the spectra after FOD. Then, four mathematical models, namely, partial least squares regression (PLSR), adaptive neural fuzzy inference system (ANFIS), random forest (RF) and generalized regression neural network (GRNN), were used to estimate the concentration of Hg, Cr and Cu. Results showed that high-order FOD had an excellent effect in highlighting hidden information and separating minor absorbing peaks, and the optimal band combination algorithm could remove the influence of spectral noise caused by high-order FOD. The incorporation of the optimal band combination algorithm and FOD is able to further mine spectral information. Furthermore, GRNN made an obvious improvement to the estimation accuracy of all studied heavy metals compared to ANFIS, PLSR, and RF. In summary, our results provided more feasibility for the rapid estimation of Hg, Cr, Cu and other heavy metal pollution areas in agricultural soils.