Classically, the gauge field strength is a curvature of a connection, the same way that the Riemann tensor is. Since $F^a_{\mu\nu}$ lives in the adjoint representation of the gauge group, you can define a 4-index object very analogous to the Riemann tensor:
$$\mathcal{F}^a{}_{b\mu\nu} \equiv F^c_{\mu\nu} f_c{}^a{}_b,$$
where $f_c{}^a{}_b$ are the structure constants defined via
$$[t_c, t_b] = f_c{}^a{}_b \, t_a,$$
modulo factors of $i$ if you care to insert them. In any case, the object $\mathcal{F}^a{}_{b\mu\nu}$ contains information about parallel transport around infinitesimal loops, in the same sense that the Riemann tensor does. But the vector being translated is not a tangent vector; instead it is a vector in gauge space (or "internal" space), whose components are defined via expansion in the generators $t_a$, as in $V = V^a \, t_a$.
So, $\mathcal{F}^a{}_{b\mu\nu}$ describes the change in $V = V^a \, t_a$ as it is parallel-transported (via the covariant derivative $D_\mu \equiv \partial_\mu + A_\mu$, again modulo factors of $i$, etc.) around a small parallelogram in the directions $\partial_\mu, \partial_\nu$.
The posed problem is formidable. If we had closed analytical expressions of the gauge invariant components of the reduced space of a gauge system, it would be a step towards the solution of the Yang-Mills Millennium problem. (We would then need to quantize this space, which is a formidable problem by itself).
Huebschmann, Rudolph and Schmidt have performed the above exercise on a (very) simplified case of a lattice gauge theory consisting of a single plaquette; and as you can see, even this case is extremely difficult. The reduced phase space in this case is
$$M_{red} = T^*T/W$$
Where $T$ is the maximal torus of the gauge group $G$, $ T^*T$ is its cotangent bundle and $W$ is the Weyl group.
Moreover, the reduced space is not a manifold. It has the structure of what is called a stratified symplectic space, which is a disjoint union of strata of the action of the gauge groups on the unreduced space. The approximate computation of the energy levels in this problem entails the computation of tunneling probabilities between the strata.
Returning to the original question in the continuum
Some of the invariants of Yang-Mills fields are can be given by characteristic classes, please see the review by Eguchi, Gilkey and Hansen (section 6). These characteristic classes are gauge differential forms consisting of traces of polynomials of the Field strengths, which in addition, when integrated over appropriate cycles of the space-time manifold, give topological invariants. For example the second Chern class
$$c_2 = \frac{1}{8 \pi^2} \left (\mathrm{tr} F \wedge F - \mathrm{tr}F \wedge \mathrm{tr}F \right )$$
However, these invariants do not separate the reduced gauge invariant space (i.e., do not exhaust the gauge invariant coordinates of the reduced space).
Another type of invariants based on holonomies involves the Yang-Mills connection (vector potential) and not the field strengths. These invariants have the form Wilson loops.
$$\mathrm{tr}\mathrm{P}e^{\int_{\Gamma} A}$$
where the integration is over a closed path $\Gamma$ and $\mathrm{P}$ denotes path ordering.
The collection of the Wilson loops over all possible paths in space-time does separate the reduced phase space. Actually, this representation is the basis of the loop formulation of the Yang-Mills theory, please see the following review by Loll. However, this approach poses great difficulties in quantization.
Best Answer
Note that, for example, \begin{align} [A_\mu,\partial_\nu]f&=A_\mu\partial_\nu f-\partial_\nu(A_\mu f)\\ &=A_\mu\partial_\nu f-\partial_\nu(A_\mu)f-A_\mu\partial_\nu f\\ &=-f\partial_\nu A_\mu\,. \end{align} So you don't get terms like $A_\mu\partial_\nu$.