Galactic foregrounds are the main obstacle to observations of the cosmic microwave background (CMB) $B$-mode polarization. In addition to obscuring the inflationary $B$-mode signal by several orders of magnitude, Galactic foregrounds have non-trivial spectral signatures that are partially unknown and distorted by averaging effects along the line-of-sight, within the pixel/beam window, and by various analysis choices (e.g., spherical harmonic transforms and filters). Statistical moment expansion methods provide a powerful tool for modeling the effective Galactic foreground emission resulting from these averaging effects in CMB observations, while blind component separation treatments can handle unknown foregrounds. In this work, we combine these two approaches to develop a new semi-blind component separation method at the intersection of parametric and blind methods, called constrained moment ILC (cMILC). This method adds several constraints to the standard ILC method to de-project the main statistical moments of the Galactic foreground emission. Applications to maps are performed in needlet space and when compared to the NILC method, this helps significantly reducing residual foreground contamination (bias, variance, and skewness) in the reconstructed CMB $B$-mode map, power spectrum, and tensor-to-scalar ratio. We consider sky-simulations for experimental settings similar to those of LiteBIRD and PICO, illustrating which trade-offs between residual foreground biases and degradation of the constraint on $r$ can be expected within the new cMILC framework. We also outline several directions that require more work in preparation for the coming analysis challenges.