Given a 2+1D Abelian K matrix Chern-Simon theory (with multiplet of internal gauge field $a_I$) partition function:
$$ Z=\exp\left[i\int\big( \frac{1}{4\pi} K_{IJ} a_I \wedge d a_J + a \wedge * j(\ell_m)+ a \wedge * j(\ell_n)\big)\right] $$
with anyons (Wilson lines) of $j(\ell_m)$ and $j(\ell_n)$.
One can integrate out internal gauge field $a$ to get a Hopf term, which we interpret as the braiding statistics angle, i.e. the phase gained of the full wave function of the system when we do the full braiding between two anyons:
$$ \exp\left[i\theta_{ab}\right]\equiv\exp\left[i 2 \pi\ell_{a,I}^{} K^{-1}_{IJ} \ell_{b,J}^{}\right] $$ see also this paper and this paper.
I would like to know the way(s) to obtain braiding statistics of anyons from a Non-Abelian Chern-Simon theory? (generically, it should be a matrix.) How to obtain this braiding matrix from Non-Abelian Chern-Simon theory?