The Chinese remainder theorem (CRT) provides an efficient way to reconstruct an integer from its remainders modulo several integer moduli, and has been widely applied in signal processing and information theory. Its multidimensional extension (MD-CRT) generalizes this principle to integer vectors and integer matrix moduli, enabling reconstruction in multidimensional signal processing scenarios. However, since matrices are generally non-commutative, the multidimensional extension introduces new theoretical and algorithmic challenges. When all matrix moduli are diagonal, the system is equivalent to applying the one-dimensional CRT independently along each dimension. This work first investigates whether non-diagonal (non-separable) moduli offer fundamental advantages over traditional diagonal ones. We show that under the same determinant constraint, non-diagonal matrices do not increase the dynamic range but yield more balanced and better-conditioned sampling patterns. More importantly, they generate lattices with longer shortest vectors, leading to higher robustness to vector remainder errors, compared to diagonal ones. To further improve the robustness, we develop a multi-stage robust MD-CRT framework that improves the robustness level without reducing the dynamic range. Due to the multidimensional nature and modulo matrix forms, it is challenging and not straightforward to extend the existing one-dimensional multi-stage robust CRT. In this paper, we obtain a new condition for matrix moduli, which can be easily checked, such that a multi-stage robust MD-CRT can be implemented. Both theoretical analysis and simulation results demonstrate that the proposed multi-stage robust MD-CRT achieves stronger error tolerance and more reliable reconstruction under erroneous vector remainders than that of single-stage robust MD-CRT.