  The characterization and modeling of polarized foregrounds has become a
critical issue in the quest for primordial $B$-modes. A typical method to
proceed is to factorize and parametrize the spectral properties of foregrounds
and their scale dependence (i.e. assuming that foreground spectra are well
described everywhere by their sky average). Since in reality foreground
properties vary across the Galaxy, this assumption leads to inaccuracies in the
model that manifest themselves as biases in the final cosmological parameters
(in this case the tensor-to-scalar ratio $r$). This is particularly relevant
for surveys over large fractions of the sky, such as the Simons Observatory
(SO), where the spectra should be modeled over a distribution of parameter
values. Here we propose a method based on the existing ``moment expansion''
approach to address this issue in a power-spectrum-based analysis that is
directly applicable in ground-based multi-frequency data. Additionally, the
method uses only a small set of parameters with simple physical interpretation,
minimizing the impact of foreground uncertainties on the final $B$-mode
constraints. We validate the method using SO-like simulated observations,
recovering an unbiased estimate of the tensor-to-scalar ratio $r$ with standard
deviation $\sigma(r)\simeq0.003$, compatible with official forecasts. When
applying the method to the public BICEP2/Keck data, we find an upper bound
$r<0.06$ ($95\%\,{\rm C.L.}$), compatible with the result found by BICEP2/Keck
when parametrizing spectral index variations through a scale-independent
frequency decorrelation parameter. We also discuss the formal similarities
between the power spectrum-based moment expansion and methods used in the
analysis of CMB lensing.
