Interval quantification for multibody systems can provide an accurate dynamic prediction and a robust reliability design. In order to achieve a robust numerical model, multiple interval uncertain parameters should be considered in the uncertainty propagation of multibody systems. The response bounds obtained by the bivariate Chebyshev method (BCM) present an intensive deterioration with the increase of time history in the interval dynamic analysis. To circumvent this problem, a novel method that combines the bivariate Chebyshev polynomial and local mean decomposition (BC-LMD) is proposed in this paper. First, the multicomponent response of the system was decomposed into the sum of several mono-component responses and a residual response, and the corresponding amplitude and phase of the mono-component were obtained. Then, the bivariate function decomposition was performed on the multi-dimensional amplitude, phase, and residual to transform a high-dimensional problem into several one-dimensional and two-dimensional problems. Subsequently, a low order Chebyshev polynomial can be used to construct surrogate models for the multi-dimensional amplitude, phase, and residual responses. Then, the entire coupling surrogate model of the system can be established, and the response bounds of the system can be enveloped. Illustrative examples of a slider-crank mechanism and a double pendulum are presented to demonstrate the effectiveness of the proposed method. The numerical results indicate that, compared to the BCM, BC-LMD can present a tight envelope in the long time-dependent dynamic analysis under multiple interval parameters.