This paper concerns a spectral estimation problem for multivariate (i.e., vector-valued) signals defined on a multidimensional domain, abbreviated as M-2. The problem is posed as solving a finite number of trigonometric moment equations for a nonnegative matricial measure, which is well known as the covariance extension problem in the literature of systems and control. This inverse problem and its various generalizations have been extensively studied in the past three decades, and they find applications in diverse fields such as modeling and system identification, signal and image processing, robust control, circuit theory, etc. In this paper, we address the challenging M-2 version of the problem, and elaborate on a solution technique via convex optimization with the tau-divergence family. As a major contribution of this work, we show that by properly choosing the parameter of the divergence index, the optimal spectrum is a rational function, that is, the solution is a spectral density which can be represented by a finite-dimensional system, as desired in many practical applications.