In order to solve the problem that it is difficult to obtain high-precision elevation outliers in complex areas, this paper proposes an elevation anomaly fitting method based on IHHO-LSSVM. Firstly, the Harris Hawk Optimization algorithm is improved using nonlinear convergence factors, jump distances, and adaptive weights; Then, the improved HHO algorithm is used to provide more accurate regularization parameters and kernel functions for the Least Squares Support Vector Machine elevation anomaly fitting model; Finally, to verify the adaptability of the elevation anomaly combination model in complex terrain, the root mean square error of the elevation anomaly values was used as the evaluation basis, and experiments were conducted using engineering case data from two different terrains. The results show that in the bridge strip area and karst surface area, compared with the HHO-LSSVM method and LSSVM method, the IHHO-LSSVM method has higher external conformity accuracy, stronger stability, and wider adaptability. The accuracy of the bridge strip area reaches 0.0101m, and the karst surface area reaches 0.0125m, which can provide certain reference value for the establishment of GNSS elevation anomaly fitting models.