Image segmentation using similarity or dissimilarity measures between probability distributions has been of great research interest in recent years. It is shown that the cross-bin metrics such as EMD is superior to the bin-wise metrics. However, existing segmentation approaches involving EMD are limited to univariate distributions, or one-dimensional marginal distributions of multidimensional features. This paper presents a novel segmentation method based on the variational EMD (VEMD) model, which can exploit joint distributions of multidimensional features. This method formulates the segmentation problem as the minimization of the EMD-based functional, which measures the distance between the foreground (resp. background) distribution and the reference foreground (resp. background) distribution. Using the simplex method and theory of shape derivative, we minimize the functional and obtain the gradient descent flow. We use a Gaussian filtering level-set method to obtain the numerical solution, in which the level-set re-initialization and smoothness constraint commonly imposed by the contour length are not necessary. Experiments show that the proposed method outperforms the state-of-the-art segmentation methods in the presence of illumination changes and noise.