Start with the linear problem: Let $K \subset \Gamma(E)$ denote the kernel of $P$ and $\hat{K} \subset \Gamma(F)$ the cokernel. Assume that there are inner products on the vector bundles $E$ and $F$ and let $K^\perp$ and $\hat{K}^\perp$ the orthogonal complements with respect to the induced $L^2$ inner product on $\Gamma(E)$ and $\Gamma(F)$.
Then the operator $\hat{P} = \hat\pi\circ P$ restricted to $K^\perp$ is 1-1 and onto. So there is an inverse operator $\hat{Q}: \hat{K}^\perp \rightarrow K^\perp$. It's easily confirmed that $\hat{Q}$ is linear. Extend it to a linear map $Q: \Gamma(F) \rightarrow \Gamma(E).
It is now straightforward to use a priori elliptic estimates (for example, using pseudodifferential operators to construct a parametrix as described by Igor Khavkine) to establish that $Q$ is a bounded operator between the appropriate Sobolev spaces. So $Q$ is the Green's operator.
Next, the smooth tame estimates, where $P$ and the inner products on $E$ and $F$ are defined in terms of a Riemannian metric. Here, things gets trickier. You can get these estimates from the pseudodifferential approach by estimating the constant in the a priori estimates in terms of the symbol of the operator. This is described in a different setting (strictly hyperbolic systems) in a Duke Math Journal paper by Bryant-Griffiths-Yang. It is however better, if possible, to find a proof of the a priori estimates without using pseudodifferential operators. My favorite approach (I'm sure there are others) is to use Moser iteration (not the same as Nash-Moser iteration). I'm sure there's an explanation of how to get smooth tame estimates out of this, but I'm not sure where. Maybe in Hamilton's original Ricci curvature paper. Or his even older harmonic maps paper.