The multi-dimensional joint probability distribution of uncertain mechanical strength parameters is crucial for reliability analysis in geotechnical engineering. However, due to the limited number of test data samples, the analysis faces significant challenges. This paper presents a modeling method for joint probability distribution based on multidimensional Gaussian copula, which is suitable for small sample conditions. Firstly, the copula theory is introduced, and a construction method for multidimensional Gaussian copula is proposed. Secondly, measurement data for 192 sets of coal seam roof strength parameters were collected from 24 coal mines in China. For the four-dimensional Gaussian copula, analyses of goodness of fit and simulation error were performed using the Pearson method, Kendall method, and Spearman method, respectively. Finally, the impact of different copula functions on the fitting capability of the original data is discussed. The results indicate that the established multi-dimensional Gaussian copula joint distribution model possesses an advantage in characterizing the uncertainty of relevant parameters. This model also offers an effective approach for studying the uncertainty of coal mine strength parameters. Different construction methods for multidimensional Gaussian copulas yield varying simulation outcomes. The Pearson method exhibits the most favorable fitting effect for constructing correlation parameters. The two rank correlation methods have the advantages of broad applicability, excellent fitting effect, and minimal variation in simulation error. For the bivariate copula fitting ability of uncertain parameters, the Gaussian Copula exhibits the best fitting ability under both positive and negative correlation scenarios. When the correlation is weak, Plackett copula and Frank copula can also serve as alternative candidate copula functions. These findings offer a significant theoretical foundation for conducting reliability analysis of a coal seam roof under limited test data conditions.