In this study, we propose an interpretable deep network for variational pansharpening (VP), named VP-Net. Different from traditional priors using linear operators, such as the gradient, we construct a prior based on the similarity between panchromatic (PAN) and high-resolution multispectral (HRMS) images by a nonlinear operator that can be learned through a deep network. Considering the spectral difference of various satellite multispectral (MS) imaging platforms, we specifically seek the aforementioned similarity from the PAN image and the intensity of the HRMS image to reduce the spectral distortion. Based on this prior, we propose a novel VP model by further incorporating a data fidelity term from the low-resolution MS image. Specifically, we build the VP-Net by unrolling the variable splitting method for an optimal solution to this model. Consequently, all modules in VP-Net have clear physical meanings and strong generalization capabilities. Meanwhile, all parameters and the aforementioned nonlinear operator are learned in VP-Net, avoiding the difficulty of selecting optimal handcrafted parameters in traditional methods. Therefore, VP-Net not only achieves an optimal balance between spatial and spectral qualities but also has a strong generalization capability across different types of training and testing data. In the experiment, we first demonstrate the superiority of the proposed method over the current state of the arts in terms of both visual effect and quantitative analysis on different satellite datasets. Moreover, we carry out an normalized difference vegetation index (NDVI) experiment to demonstrate its potential in remote sensing.