*OLDER EDIT. (Elementary derivation)* I realized that my original answer was actually *overkill* for the question. The said integral in question follows from the definition of the *multivariate Gamma function*

\begin{equation*}
\Gamma_p(a) := \int_{A > 0} \exp(-\mbox{tr}(A))\det(A)^{a-(p+1)/2}(dA),
\end{equation*}
where $\Re(a)>(p-1)/2$.

From this it follows (by a change of variables) that for a positive definite matrix $S$,
\begin{equation*}
\int_{A > 0} \exp(-\mbox{tr}(S^{-1}A))\det(A)^{a-(p+1)/2}(dA) = \Gamma_p(a)\det(S)^a,
\end{equation*}
so that with $S=(I-U)^{-1}$ we obtain the integral in question.

Of course, to complete the picture it may be helpful to express $\Gamma_p(a)$ in more elementary terms. Chapter 2 of Muirhead's book provides these details. I cite the result that provides this expression.

**Theorem (Muirhead (1982), Thm 2.1.2)** Let $\Re(a) > (p-1)/2$. Then,
$$
\Gamma_p(a) = \pi^{p(p-1)/4}\prod_{j=1}^p \Gamma(a - (j-1)/2)
$$

(Hint: To prove the above, write the *Cholesky decomposition* $A=T'T$, with these change of variables, the original Gamma integral factorizes into a product of Gaussian and Gamma integrals.)

The part that I recalled below provides yet another representation that expresses the multiplicative determinantal lhs in terms of an infinite sum.

**OLDER STUFF**

This is actually somewhat classical knowledge. Here are two related pointers.

A *partition* $\tau=(t_1,\ldots,t_m)$ is a vector of nonnegative integers listed in increasing order, and $|\tau|$ denotes $t_1+\cdots+t_m$. The *generalized Pochhammer* symbol $(a)_\tau$ is defined as
\begin{equation*}
\newcommand{\risingf}[2]{{{#1}}^{\overline{{#2}}}}
(a)_\tau := \frac{\Gamma_d(a+\tau)}{\Gamma_d(a)} = \prod_{l=1}^m \risingf{\bigl(a - \tfrac{1}{2}(l-1)\bigr)}{t_l}
\end{equation*}

Let $C_\tau(X)$ be the *Zonal Polynomial* with signature partition $\tau$. Then, the following representation exists

For a matrix $U$ satisfying $\|U\| < 1$, we have the following "binomial-theorem"

\begin{equation}
\frac{1}{|I-U|^a} = \sum_{k\ge 0}\sum_{|\tau| = k} \frac{(a)_\tau C_\tau(U)}{k!}.
\end{equation}

Using representations for these Zonal polynomials, one can obtain the integral representation mentioned in the original post.

More directly, you can look at Chapter 7 of *R. Muirhead, "Aspects of Multivariate Statistical Theory"*, where you'll see that actually, $|I-U|^{-a}={}_1F_0(a;U)$, a matrix argument hypergeometric function. I've to run now, if I get a chance I'll clean up my answer and fill in the details.