It is possible to introduce the gauge field in a QFT purely on geometric arguments. For simplicity, consider QED, only starting with fermions, and seeing how the gauge field naturally emerges. The observation is that the derivative of the Dirac field doesn't have a well-defined transformation, because: $$n^\mu \partial_\mu \,\psi = \lim\limits_{\epsilon\rightarrow 0}\big[\psi(x+\epsilon n)-\psi(x)\big],$$ i.e. the derivative combines two fields at different spacetime points (having different transformation rules). We need to introduce a parallel transporter ##U(y,x)## that transforms as $$U(y,x) \rightarrow e^{ig\alpha(y)}U(y,x) e^{-ig \alpha(x)},$$ such that we can adapt the defintion of the derivative into a covariant derivative, that transforms in a well-defined way: $$n^\mu D_\mu \,\psi = \lim\limits_{\epsilon\rightarrow 0}\big[\psi(x+\epsilon n)- U(x+\epsilon n,x)\psi(x)\big].$$ From geometric arguments, it is straightforward to show that the parallel transporter is a Wilson line: $$U(y,x) = \mathcal{P}\,e^{ig\int\limits_x^y dz^\mu\, A_\mu(z)},$$ which introduces a new field, namely the gauge field ##A_\mu##. See e.g. Peskin & Schroeder Chapter 15 for more detail. However.. Where the interaction term ##\bar{\psi} A_\mu \psi## emerged in a natural way, I totally don't see how the kinetic terms emerge. The standard way to proceed, is to consider a Wilson loop (a Wilson line on a closed path), and use Stokes' theorem: $$\text{exp}\left\{ig\oint_\mathcal{C}dx^\mu\, A_\mu \right\} = \text{exp}\left\{ig\int_\Sigma dx^\mu \wedge dx^\nu\,\left(\partial_\mu A_\nu - \partial_\nu A_\mu \right)\right\},$$ where of course ##\partial_\mu A_\nu - \partial_\nu A_\mu \equiv F_{\mu\nu}##. In Peskin & Schroeder, they then consider a small rectangular loop, and see that in the limit ##\epsilon\rightarrow 0##, ##F_{\mu\nu}## is invariant. But what's the point? I mean, the transformation law for ##A_\mu## is easily calculated from the definition of the Wilson loop: $$A_\mu\rightarrow A_\mu+\partial_\mu \alpha,$$ making ##F_{\mu\nu}## invariant by definition: $$F_{\mu\nu}\rightarrow \partial_\mu A_\nu - \partial_\nu A_\mu +\square \alpha-\square \alpha. $$ I would have liked to see a calculation, starting from a particular loop parameterisation, that naturally leads to the correct kinetic terms in the Lagrangian, as was the case for the interaction term. In other words $$\text{exp}\left\{ig\oint_\mathcal{C}dx^\mu\, A_\mu \right\} \leadsto -\frac{1}{4}\left(F_{\mu\nu}\right)^2,$$ but I have no idea how to do it. Or is the idea simply 'look I have found some quadratic derivative terms that are invariant, now let me fiddle a bit and put its square in ##\mathcal{L}##'? If yes, then why do Peskin & Schroeder bother calculating a loop parameterisation (p484), if using Stokes' theorem would have been enough to find ##F_{\mu\nu}## somewhere?
The point is just to show you the geometrical meaning of the field strength tensor. This construction demonstrates that the field strength tensor is a measure of curvature, in the differential geometry sense. The justification for the kinetic term is this: if you want the new field ##A_\mu## to be dynamical, then you have to add to your Lagrangian all dimension-4 operators that involve ##A_\mu## and are consistent with gauge symmetry and all the other symmetries. There are two terms that satisfy these requirements: ##\bar \psi \gamma^\mu D_\mu \psi## and ##F_{\mu \nu} F^{\mu \nu}##. Both of these terms must appear or else your theory will not be renormalizable.