Texture characterization of natural images using the mathematical framework of multifractal (MF) analysis, enables the study of the fluctuations in the regularity of image intensity. Although successfully applied in various contexts, the use of MF analysis has so far been limited to the independent analysis of a single image, while the data available in applications are increasingly multivariate. This paper addresses this limitation and proposes a joint Bayesian model and associated estimation procedure for MF parameters of multivariate images. It builds on a recently introduced generic statistical model that enabled the Bayesian estimation of MF parameters for a single image and relies on the following original key contributions: First, we develop a novel Fourier domain statistical model for a single image that permits the use of a likelihood that is separable in the MF parameters via data augmentation. Second, a joint Bayesian model for multivariate images is formulated in which prior models based on gamma Markov random elds encode the assumption of the smooth evolution of MF parameters between the image components. The design of the likelihood and of conjugate prior models is such that exploitation of the conjugacy between the likelihood and prior models enables an ecient estimation procedure that can handle a large number of data components. Numerical simulations conducted using sequences of multifractal images demonstrate that the proposed procedure significantly outperforms previous univariate benchmark formulations at a competitive computational cost.