Expanding on my previous comment. There are efficient algorithms to solve linear systems of the form
$$
(K \otimes M+I_T\otimes \Sigma)\operatorname{vec}(X) = \operatorname{vec}(B),
$$
so in practice you can just use those instead of an explicit inverse. Explicit inverses are overrated, in some fields. :)
The trick is that these linear systems are equivalent to the matrix equation
$$
MXK^\top + \Sigma X I_T^\top = B,
$$
which is known as generalized Sylvester equation.
There are around efficient algorithms to solve these equations in $O(\text{larger_dimension_of_B}^3)$ For instance, see https://dl.acm.org/citation.cfm?id=146929 or https://people.cs.umu.se/isak/recsy/ . Basically the idea is reducing all coefficients to triangular form simultaneously using two QZ decompositions, and then the resulting system can be solved by back-substitution.
In your case you can get a slight simplification, by exploiting the fact that one of the coefficients is an identity. Set $\Sigma X = Y$, to get
$$
(M\Sigma^{-1})YK^\top + Y = B,
$$
which is another classical matrix equation (a discrete-time Sylvester equation), and can be solved with a slightly faster algorithm, which is also implemented in a Matlab one-liner as Y = dlyap(-M/Sigma, K', -B)
. Unfortunately, a solver for this variant is not available in Numpy/Scipy, up to my knowledge; it has only a solver for the classical Sylvester equation, which you can use with the approach described in @Reid.Atcheson's answer.
In the case in which $M\Sigma^{-1}$ and $K^\top$ have all eigenvalues inside the unit circle, you can also write a closed-form solution for this equation as the infinite series
$$
Y = \sum_{k=0}^\infty (-M\Sigma^{-1})^kB(K^\top)^k.
$$
This follows from applying the Neumann series $(I-M)^{-1}=\sum_{k=0}^\infty M^k$ to the Kroneckerized matrix $(I_{T^2} - K\otimes (-M\Sigma^{-1}))$.