Abstract: In order to effectively address the challenge of obtaining high-precision elevation outliers in complex geographical areas, this paper proposes an innovative elevation anomaly fitting method based on IHHO-LSSVM. The study begins with an improved Harris Hawk Optimization (HHO) algorithm through the implementation of nonlinear convergence factors, optimized jump distances, and adaptive weights. These improvements significantly enhance the algorithm’s ability to escape local optima and improve convergence efficiency, thereby providing a more robust optimization framework for subsequent model parameter tuning. Subsequently, the improved HHO algorithm is employed to determine more accurate regularization parameters and kernel functions for the Least Squares Support Vector Machine (LSSVM) elevation anomaly fitting model. This optimization process ensures that the LSSVM model achieves higher precision and better generalization capabilities in elevation anomaly fitting tasks. To thoroughly validate the adaptability and robustness of the proposed elevation anomaly combination model in complex terrains, extensive experiments were conducted using engineering case data from two distinct geographical regions: a bridge strip area and a karst surface area. The evaluation was based on the Root Mean Square Error (RMSE) of the elevation anomaly values as the primary metric, with additional consideration given to computational efficiency and model stability. The experimental results demonstrate that in both the bridge strip area and karst surface area, the IHHO-LSSVM method outperforms the conventional HHO-LSSVM and standard LSSVM methods in terms of external conformity accuracy, stability, and adaptability. Specifically, the IHHO-LSSVM method achieves remarkable accuracy levels of 0.0101 meters in the bridge strip area and 0.0125 meters in the karst surface area, representing significant improvements over traditional methods. Furthermore, the proposed method exhibits superior stability across different terrain types, with reduced variance in prediction errors. These findings not only highlight the superior performance of the proposed method but also provide valuable insights and a reliable reference for the establishment of GNSS elevation anomaly fitting models in various complex terrains. The study contributes to the field of geodetic surveying by offering a more precise and robust solution for elevation anomaly fitting, particularly in challenging geographical conditions.