Your second-order differential operator appears when one takes the variation of the scalar curvature of the sphere by a symmetric tensor $h\mapsto d/dt|_{t=0}\mathrm{Sc}_{g+th}$. In constant sectional curvature, such operators were studied for instance by Calabi (60'), and then in constant scalar curvature by Ebin (69'). Here is what I suggest to solve your question in the special case of $S^{2}$.

If you have a Riemmanian manifold with constant sectional curvature, $(M,g)$, consider the second order differential operator $H_{g}=\frac{1}{2}{(d^{\nabla^{g}}d^{\nabla^{g}}_{V}+d^{\nabla^{g}}_{V}d^{\nabla^{g}}})-\frac{1}{2}\kappa g\wedge$. 
$\kappa$ here is the sectional curvature of the space, so in your case it can be taken to be $\kappa=1$. The wedge product here is known as the Kulkarni-Numizu product of symmetric tensor fields (there is a more generlized version of this wedge product, operating on "double forms", also adressed by Kulkarni, but this notion will do). $d^{\nabla^{g}}$ and $d^{\nabla^{g}}_{V}$ here are the exterior covariant derviatvie of the first and second index of a symmetric tensor, respectively, where we think of symmetric tensors as examples of vector-valued differential forms $\Omega^{1}(M ; T^{*}M)$. The image of $H_{g}$ lies in $\Omega^{2}(M;\Lambda^{2}T^{*}M)$, and in the case where $M$ is two-dimensional every element in this space is fully determined by its "scalar curvature", namely if $\sigma\in\Omega^{2}(M;\Lambda^{2}T^{*}M)$ then $\sigma=\frac{1}{4}(tr_{g}tr_{g}\sigma) g\wedge g$.

In the case where $\kappa=1$, a direct calcultion shows that taking the trace twice from $H_{g}$ yields the operator $tr_{g}tr_{g}H_{g}=-\frac{1}{2}\Delta_{g}tr_{g}+\mathrm{div}\mathrm{div}-tr_{g}$, so solving your equation is equivelent of solving $H_{g}h=\frac{1}{4} f g\wedge g$, where $g$ is the metric of the sphere. Replacing $h=H^{*}_{g}\psi$ for $\psi\in\Omega^2(M;\Lambda^{2}T^{*}M)$ where $H_{g}^{*}$ is the formal $L^{2}$ dual of $H_{g}$ yields the equation $H_{g}H^{*}_{g}\psi=\frac{1}{4} f g\wedge g$. Note how since $H_{g}$ operates on symemtric tensors, the image of $H^{*}_{g}$ is a symmetric tensor in $\Omega^{1}(M;T^{*}M)$. 

Another calculation then shows that the principle symbol of $H_{g}H^{*}_{g}$ in the case where $M$ is two dimensional is $|\xi|^{4}$. Thus this is an elliptic fourth order differential opertor (a "bilaplcian"), and so the equation $H_{g}H^{*}_{g}\psi=\frac{1}{4} f g\wedge g$ is solvable for $\psi$ if and only if $ \frac{1}{4} f g\wedge g$ is orthogonal to the kernel of $H_{g}H^{*}_{g}$. By duality, this kernel is equal to $\mathrm{ker} H^{*}_{g}$. I am not sure if this kernel is trivial when $M$ is simply connected. If not, then this orthogonality yields a compatability condition which $f$ must satisfy in order for the equation to have a solution.