Potential field methods play a crucial role in mineral exploration, aiding in mapping sub-surface physical properties like density and magnetic susceptibility. However, modeling and interpreting potential field data can be challenging due to non-uniqueness, leading to uncertainty about the true subsurface distribution. Additionally, gravity and magnetic data lack inherent depth information, often resulting in recovered models toward the surface. To address these limitations, this study combines statistical analysis with geological constraints. Geological knowledge is incorporated during the inversion process as a constraint to recover better depth, shape, and size consistent with geological knowledge. This counteracts the decaying kernel influence on observations with depth and distance. Also, it reduces the number of iterations for better computational efficiency and geologically plausible solutions. However, gravity and magnetic data often contain noise and outliers that can lead to inaccurate models. To mitigate this, we introduce a statistical approach, i.e., the probabilistic bound method, that minimizes the impact of outliers on the traditional weighted least squares inversion. The effectiveness of this combined approach is evaluated using synthetic datasets simulating simple to complex scenarios at varying depths and noise levels. Our results demonstrate significant improvements in delineating the shape and depth of subsurface structures compared to traditional methods. Applying the technique to real field data confirms its potential to produce more compact and improved subsurface models. This approach holds significant promise for enhancing the characterization of ore bodies and other geological features based on potential field data, particularly in recovering their shapes, geometries, and depths.