Consider a supercritical d-type branching process $Z_n^{i} = $ $ (Z_n^i(1), \cdots, Z_n^i(d)), \,n\geq 0,$ in an i.i.d. environment $\xi =(\xi_0, \xi_1, \ldots)$, starting with one particle of type $i$, whose offspring distributions of generation $n$ depend on the environment $\xi_n$ at time $n$. In the deterministic environment case, the famous Kesten-Stigum (1966) theorem states essentially that, if the mean matrix of the offspring distribution has spectral radius $\rho >1$, then for all $i,j =1, \ldots, d$, almost surely $ \lim_{n\rightarrow \infty} \frac{Z_{n}^{i}(j)}{\rho^n }$ exists and is finite; moreover, the limit variables are non-degenerate if and only if $\mathbb{E} \left(Z_1^{i}(j) \log^{+} Z_1^{i}(j)\right)<+\infty$ for all $i$ and $j$. The extension to the random environment case with $d=1$ has been done by Athreya and Karlin (1971) and Tanny (1988). Extending the Kesten-Stigum theorem to the random environment case with $d>1$ is a long-standing problem. The main objective of this paper is to resolve this delicate problem. In particular, under simple conditions, we prove that for any $1\le i,j\le d$, as $n \rightarrow +\infty$, $Z_n^{i}(j)/\mathbb{E}_\xi Z_n^{i}(j) \rightarrow W^{i}$ in probability, where $W^i$ is a non-negative random variable, $ \mathbb{E}_\xi Z_n^{i}(j)$ is the conditional expectation of $Z_n^{i}(j)$ given the environment $\xi$, which diverges to $\infty$ with geometric rate in the sense that $ \frac{1}{n} \log \mathbb{E}_\xi Z_n^{i}(j) \rightarrow \gamma >0$ almost surely, $\gamma$ being the Lyapunov exponent of the mean matrices of the offspring distributions; moreover $W^i$ are non-degenerate for all $i $ if and only if $\mathbb{E}\left(\frac{Z_1^{i}(j)}{M_0(i,j)}\log^{+}\frac{Z_1^{i}(j)}{M_0(i,j)}\right)<+\infty$ for all $i$ and $j$, where $M_0(i,j)$ is the conditioned mean of the number of children of type $j$ produced by a particle of type $i$ at time $0$, given the environment $\xi$. The key idea of the proof is the introduction of a non-negative martingale $(W_n^{i})$ which converges a.s. to $W^i$, and which reduces to the well-known fundamental martingale in the deterministic environment case. In addition, we prove that the direction $Z_n^{i}/ \| Z_n^{i} \|$ converges in law conditioned on the explosion event $\{\|Z_n^i\| \rightarrow +\infty\}$. The case of stationary and ergodic environment is also considered. Our results open ways to prove important properties such as central limit theorems with convergence rate and large deviation asymptotics.