This is explained nicely in arXiv:1012.0461 [cond-mat.stat-mech] (section 2.1) in the case where there are no thermal effects. I leave out the thermal effects in the following. Since you know your paper better than me you should be able to generalise this.
The projection comes from the incompressibility of the fluid. The Navier-Stokes equations read
$$ \partial_t v^\alpha + v^\beta \partial_\beta v^\alpha = \nu \partial_\beta \partial_\beta v^\alpha - \partial_\alpha P + f^\alpha \, , \qquad \partial_\beta v^\beta = 0 \, .$$
$\nu$ is the kinematic viscosity, $\partial_t = \partial/\partial t$, $\partial_\beta = \partial/\partial x_\beta$ and upper indices indicate the components of the corresponding vector. Repeated indices are summed over.
The pressure can be eliminated from these equations by taking the divergence of Navier-Stokes equations,
$$\partial_\alpha \left(\partial_t v^\alpha + v^\beta \partial_\beta v^\alpha = \nu \partial_\beta \partial_\beta v^\alpha - \partial_\alpha P + f^\alpha\right) \, .$$
Taking incompressibility into account, this leads to
$$\partial_\alpha \partial_\alpha P = -\partial_\alpha v^\beta \partial_\beta v^\alpha + \partial_\alpha f^\alpha \, , \text{ or } \, P = - \frac{1}{\nabla^2} \left(\partial_\gamma v^\beta \partial_\beta v^\gamma - \partial_\gamma f^\gamma \right)\, .$$
Then Navier-Stokes equations are written as
$$ \partial_t v^\alpha + \left(\delta_{\alpha,\gamma}-\partial_\alpha \frac{1}{\nabla^2} \partial_\gamma \right) v_\beta \partial_\beta v^\gamma = \nu \partial_\beta \partial_\beta v^\alpha +\left(\delta_{\alpha,\gamma}-\partial_\alpha \frac{1}{\nabla^2} \partial_\gamma \right) f^\gamma \, , $$
after inserting the second expression for the pressure. Finally, since the velocity field is incompressible,
$$ \left(\delta_{\alpha,\gamma}-\partial_\alpha \frac{1}{\nabla^2} \partial_\gamma \right) v^\gamma = v^\alpha \, ,$$
we can define new fields as
$$ \hat{v}^\alpha = \left(\delta_{\alpha,\gamma}-\partial_\alpha \frac{1}{\nabla^2} \partial_\gamma \right) v^\gamma \, , \text{ and } \hat{f}^\alpha = \left(\delta_{\alpha,\gamma}-\partial_\alpha \frac{1}{\nabla^2} \partial_\gamma \right) f^\gamma \, ,$$
in terms of which Navier-Stokes equations become
$$ \partial_t \hat{v}^\alpha + \hat{v}^\beta \partial_\beta \hat{v}^\alpha = \nu \partial_\beta \partial_\beta \hat{v}^\alpha + \hat{f}^\alpha \, , \qquad \partial_\beta v^\beta = 0 \, .$$
The pressure term has been completely eliminated. The price to pay is that the velocity field is implicitly constrained to be incompressible and that the force is replaced by $\hat{f}$.
If you compute the correlator of $\hat{f}$ in Fourier space you get,
$$ \langle \hat{f}^\alpha(p) \hat{f}^\beta(q) \rangle = \int_{x,y} \text{e}^{i px + q y} \left(\delta_{\alpha,\gamma_1}-\partial_\alpha \frac{1}{\nabla^2} \partial_{\gamma_1} \right) \left(\delta_{\beta,\gamma_2}-\partial_\beta \frac{1}{\nabla^2} \partial_{\gamma_2} \right) \langle f^{\gamma_1}(x) f^{\gamma_2}(y) \rangle$$
$$= \left(\delta_{\alpha,\gamma_1}- \frac{p^\alpha p^{\gamma_1}}{p^2} \right) \left(\delta_{\beta,\gamma_2}-\frac{q^\beta q^{\gamma_2}}{q^2} \right) \int_{x,y} \text{e}^{i px + q y} \langle f^{\gamma_1}(x) f^{\gamma_2}(y) \rangle$$
Finally, inserting
$$ \int_{x,y} \text{e}^{i px + q y} \langle f^{\gamma_1}(x) f^{\gamma_2}(y) \rangle = \delta_{\gamma_1,\gamma_2} \delta(p+q) D_0(p) \, ,$$
and contracting the indices leads to your expression for the correlator of $\hat{f}$.