Aiming at the problems that the traditional single ellipsoid fitting fails to eliminate the errors of the geomagnetic sensor and does not consider the geomagnetic anomaly information, etc., this paper proposes a Robust Least Squares combined with the Levenberg-Marquardt (LM-RLS) algorithm to fit an accurate ellipsoid to achieve accurate data calibration. According to the geomagnetic sensor error model, the total amount of geomagnetism is used as the constraint, the fitted ellipsoid coefficients of Robust Least Squares (RLS) are used as the initial values, and the LM algorithm is used to calculate the magnetic field dispersion points to the nearest point of the fitted ellipsoid surface. Then, the initial calibration value is solved by the criterion of the minimum weighted sum of squares of the closest point distance. Finally, the accurate ellipsoid is obtained to realize data calibration. The results of the dynamic experiment show that after the ellipsoid is precisely fitted by the LM-RLS algorithm, the sum of the absolute values of the distance from the magnetic field discrete data to the ellipsoid surface is reduced by 48.33% (RLS) and 43.10% (RWTLS). The relative errors of the magnetic field in the X, Y, and Z axes are reduced from 0.48%, 0.78%, 1.48% (RLS) and 0.30%, 0.58%, 1.14% (RWTLS) to 0.10%, 0.04%, 0.16 %, respectively. The accuracy of geomagnetic calibration is improved by 68.06% (RLS) and 61.02% (RWTLS) in the dynamic experiment. The algorithm improves the accuracy of geomagnetic measurement, which provides a new idea for the calibration of three-axis geomagnetic sensors.