The term "Yang-Mills theory" in the mass gap problem refers to a particular QFT. It is believed that this QFT (meaning its Hilbert space of states and its observable operators) should be defined in terms of a measure on the space of connections on $\mathbb{R}^4$; roughly speaking, the moments of this measure are the matrix elements of the operator-valued distributions. (People also use the term "gauge theory" to refer to any QFT, like QCD, which has a Yang-Mills sub-theory.)
The mass gap problem really has two aspects: First, one has to construct an appropriate measure $d\mu$ on some space of connections. Then, one has to work out which functions on the space of connections are integrable with respect to this measure, and show that the corresponding collection of operators includes an energy operator (i.e. a generator of time translations) which has a gap in its spectrum.
You'll have to read the literature to really learn anything about this stuff, but I can make a few points to help you on your way. Be warned that what follows is a caricature. (Hopefully, a helpful one for people trying to learn this stuff.)
About the Measure:
First, the measure isn't really defined on the space of connections. Rather, it should be defined on space $\mathcal{F}$ of continuous linear functionals on the space $\mathcal{S}$ of smooth rapidly-vanishing $\mathfrak{g}$-valued vector fields on $\mathbb{R}^4$, where $\mathfrak{g}$ is the Lie algebra of the gauge group $G$. The space ;$\mathcal{F}$ contains the space of connections, since any connection on $\mathbb{R}^4$ can be written as $d$ plus a $\mathfrak{g}$-valued $1$-form and paired with a vector field via the Killing form, but it also has lots of more "distributional" elements.
We're supposed to get $d\mu$ as the "infinite-volume/continuum limit" of a collection of regularized measures. This means that we are going to write $\mathcal{S}$ as an increasing union of chosen finite-dimensional vector spaces $\mathcal{S}(V,\epsilon)$; these spaces are spanned by chosen vector fields which have support in some finite-volume $V \subset \mathbb{R}^4$ and which are slowly varying on distance scales smaller than $\epsilon$. (You should imagine we're choosing a wavelet basis for $\mathcal{S}$.) Then we're going to construct a measure $d\mu_\hbar(V,\epsilon)$ on the finite dimensional dual vector space $\mathcal{F}(V,\epsilon) = \mathcal{S}(V,\epsilon)' \subset \mathcal{F}$; these measures will have the form $\frac{1}{\mathcal{Z}(V,\epsilon,\hbar)} e^{-\frac{1}{\hbar}S_{V,\epsilon,\hbar}(A)}dA$. Here, $dA$ is Lesbesgue-measure on $\mathcal{F}(V,\epsilon)$, and $S_{V,\epsilon,\hbar}$ is some discretization of the Yang-Mills action, which is adapted to the subspace $\mathcal{F}(V,\epsilon)$. (Usually, people use some version of Wilson's lattice action.)
Existence of the Yang-Mills measure means that one can choose $S_{V,\epsilon}$ as a function of $V$,$\epsilon$, and $\hbar$ so that the limit $d\mu_\hbar$ exists as a measure on $\mathcal{F}$ as $vol(V)/\epsilon \to \infty$. We also demand that the $\hbar \to 0$ limit of $d\mu_\hbar$ is supported on the space of critical points of the classical Yang-Mills equations. (We want to tune the discretized actions to fix the classical limit.)
About the Integrable Functions:
Generally speaking, the functions we'd like to integrate should be expressed in terms of the "coordinate functions" which map the $A$ to $A(f)$, where $f$ is one of the basis elements we used to define the subspaces $\mathcal{S}(V,\epsilon)$. You should imagine that $f$ is a bump vector field, supported near $x \in \mathbb{R}^4$ so that these functions approximate the map sending a $\mathfrak{g}$-valued $1$-form to the value $A_{i,a}(x)$ of its $(i,a)$-th component.
There are three warnings to keep in mind:
First, we only want to look at functions on $\mathcal{F}$ which are invariant under the group of gauge transformations. So the coordinate functions themselves are not OK. But gauge invariant combinations, like the trace of the curvature at a point, or the holonomy of a connection around a loop are OK.
Second, when expressing observables in terms of the coordinate functions, you have to be careful, because the naive classical expressions don't always carry over. The expectation value of the function $A \mapsto A_{i,a}(x)A_{j,b}(y)$ with respect to $d\mu_\hbar$ (for $\hbar \neq 0$) is going to be singular as $x \to y$. This is OK, because we were expecting these moments to define the matrix elements of operator-valued distributions. But it means we have to be careful when considering the expectation values of functions like $A \mapsto A(x)^2$. Some modifications may be required to obtain well-defined quantities. (The simplest example is normal-ordering, which you can see in many two-dimensional QFTs.)
Finally, the real problem. Yang-Mills theory should confine. This means, very very roughly, that there are some observables which make sense in the classical theory but which are not well-defined in the quantum theory; quantum mechanical effects prevent the phenomena that these observables describe. In the measure theoretic formulation, you see this by watching the expection values of these suspect observables diverge (or otherwise fail to remain integrable) as you approach the infinite-volume limit.
About the Operators:
In classical Yang-Mills theory, the coordinate observables $A \mapsto A_{i,a}(x)$ satisfy equations of motion, the Yang-Mills equations. Moreover, in classical field theory, for pure states, the expectation value of a product of observables $\mathcal{O}_1\mathcal{O}_2...$ is the product of the individual expectation values. In quantum YM, the situation is more complicated: coordinate observables may not have be well-defined, thanks to confinement, and in any case, observables only satisfy equations of motion in a fairly weak sense: If $\mathcal{O}_1$ is an expression which would vanish in the classical theory thanks to the equations of motion, then then the expectation value of $\mathcal{O}_1\mathcal{O}_2...$ is a distribution supported on the supports of $\mathcal{O}_2...$.
