The mathematical description of many physical and engineering problems requires the solution of partial differential equations with different operators or discontinuous coefficients, leading to computational domains partitioned into several, possibly time-dependent, regions. In this context, the fictitious domain approach coupled with distributed Lagrange multipliers (FD-DLM) represents a powerful technique. Applications of FD-DLM include particulate flow problems and, more recently, elliptic interface and fluid-structure interaction (FSI) problems. While particularly effective for complex geometries, FD-DLM formulations give rise to large and ill-conditioned saddle point linear systems, making efficient preconditioning essential. Suitable preconditioners have been proposed in the literature. For instance, in the context of FSI, block preconditioners obtained by neglecting selected blocks of the system matrix have been investigated, using both sparse direct solvers and multigrid-based iterative techniques. However, ensuring robustness in the presence of possibly large coefficient jumps or singular blocks remains a key challenge. We address these issues by proposing an augmented Lagrangian-based preconditioner. We also perform a spectral analysis that derives bounds for the eigenvalues of the preconditioned matrix and shows that these bounds are independent of discretization parameters when the ideal preconditioner is employed. Cheaper and modified variants are also considered to reduce computational costs. The robustness and efficiency of the preconditioner, when used in combination with flexible GMRES, are validated through extensive two- and three-dimensional numerical experiments using the C++ finite element library DEAL.II.