Remanent magnetization is a common rock physical property, and constructing a strong remanent magnetic field source model is crucial for deep geological structure and mineral resource exploration. Sparse inversion of magnetic anomaly modulus is an effective method for generating remanent magnetization models with high resolution, which essentially minimizes the L₀ norm of the model. However, the computation of the L₀ norm is challenging, and traditional methods approximate it using convex functions. While this allows for stable gradient calculation of the objective function, it is difficult to control the volume of the inverted model. This paper proposes a magnetic anomaly modulus inversion optimization model that includes data fitting, model constraints, and L₀ norm constraints. The method utilizes an iterative hard thresholding approach to achieve stable iterative calculation of the L₀ norm, enabling clear delineation of the inversion model boundary under strong remanent magnetization conditions while controlling its scale. The study shows that this method ensures accurate data fitting and achieves asymptotic linear convergence of model scale, leading to a strong remanent magnetization inversion model of the target volume. The method’s effectiveness and stability were validated through theoretical model experiments. It was then applied to the Yeshan area in Jiangsu, where a magnetic model of the mineralization body was constructed, providing scientific evidence for deep mineral exploration.