It sounds like what you would like is a rigorous version of Witten's original Feynman integral & Wilson loop approach. This is not a totally unreasonable thing to ask for, since QFTs have been rigorously constructed along these lines, by demonstrating the existence of moments of the Euclidean-signature Feynman integrals. It would indeed be wonderful to just run some Monte Carlo simulations and see knot invariants emerge.
But there are major obstacles to this program for the particular QFT Witten was thinking about. It is not known even at a physical level of rigor how to define the SU(2) Chern-Simons theory in this fashion. It's an open research topic. In fact, this question is still under active investigation by physicists, even for the far simpler abelian Chern-Simons theory. (See Bietenholz, Nishimura, Sodano '02 and DeMarco, Wen '19.)
Roughly, what one would like is a lattice gauge theory which flows to Chern-Simons of the desired level at large distances. There's a number of obstacles here. I'm only going to mention some of them.
The first sign of the trouble is that, in the continuum, the Euclidean Chern-Simons action is purely imaginary and unbounded, rather than real and pleasantly bounded like the Yang-Mills action. Euclidean path integral constructions usually depend on having Gaussian-ish behavior, and there's no reason to expect that here.
To make matters worse, there isn't any exact lattice version of the Chern-Simons class. You can approximate it in various ways, but the errors in the approximation will make the lattice theory non-topological.
Even the abelian Chern-Simons theory (which has a quadratic classical action) has problems on the lattice. There's a natural discretization of the classical action, due to Frohlich & Marchetti, but it has a fatal flaw: It has zero modes, which means that the exponential $e^{-S(A)}$ is not integrable, even on a finite lattice, where the integration measure $dA$ is well-defined.
There's a theorem due to Berruto, Diamanti, & Sidano which says, in effect, that any lattice action for a gauge field which is gauge-invariant, linear in derivatives, and odd under parity reversal (as the cubic Chern-Simons action should be) will share these problems. You can't fix this by screwing around with the discretization, and there's not much reason to expect better behavior in the case of non-abelian gauge groups.
There is a reasonable way out of this problem: One can add to the lattice Chern-Simons action a parity-even term, which lifts the degeneracies. The leading order parity-even term you can add is the Maxwell action. The abelian theory is basically OK at this point. The Maxwell term's effect on long-distance physics is negligible, so you basically have a definition of abelian Chern-Simons. I think one could even construct it along the lines of Driver and Gross. (Frohlich & King sketch the idea.) But the analogue of the Jones polynomials here are just linking numbers.
The nonabelian case is probably worse: The Maxwell term that showed up in regularization gets replaced by the Yang-Mills Lagrangian. The basic story isn't expected to change: the long distance scaling limit should be dominated by Chern-Simons. So that's good. But now you've got to show that 3d Yang-Mills theory makes sense. This is expected to be easier than the Millenium Prize problem, but it shares a lot of the same features.
So as not to end on a down-note, let me mention that Yang-Mills theory is not the only possible lattice QFT that leads to Chern-Simons. There's a quite large line of research (starting, I think, with
Fradkin & Shaposnik) in which fairly accessible QFTs, like the 3d Thirring model, are shown to give rise to Chern-Simons at long distances. I don't know to what degree the constructive QFT community has taken up these challenges; I've seen it mainly from lattice and condensed matter people.
