A seismic attenuation estimation approach is proposed based on multivariate variational mode decomposition (MVMD). MVMD, as a multivariable or multichannel signal processing tool, can extract several predefined multivariable modulation oscillations from the seismic data. Intrinsic mode functions (IMFs) obtained by MVMD may highlight some seismic attenuation information in certain frequency bands with lateral continuity. For improving the accuracy of absorption gradient estimation, we define an absorption gradient by using the low-frequency and high-frequency segments of the amplitude spectrum extracted from the time-frequency map, which is generated by applying the Teager energy separation algorithm to each IMF. A weighted factor is used for highlighting the different attenuation contributions of each IMF. Model test and field data examples show that the proposed method can give more accurate absorption gradient estimation results and yield better hydrocarbon-prone interpretations than the traditional methods.