## 2 Methods
### Multivariate Variational Mode Decomposition  
The goal of multivariate variational mode decomposition (MVMD) is to analyze a multivariate input signal containing \(M\) number of data channels into a predefined number \(K\) of multivariate components. As an extension of the VMD method, \(K\) multivariate components that will fully construct the input signal \(\mathbf{x}(t)\) are searched:  
\[\mathbf{x}(t)=\sum\limits_{k=1}^{K}\mathbf{u}_{k}(t). \tag{8}\]  
where the \(k\)th multivariate component. \(\mathbf{u}_{k}(t)=[u_{k,1}(t),u_{k,2}(t),\ldots,u_{k,M}(t)]\) is a vector with \(M\) components.  
MVMD takes the multivariate input signal as a whole and tries to seek \(K\) number of multivariate components from the input signal, with a minimum sum of band widths [24]. The optimization problem for MVMD can be naturally expanded from Equation (2) as follows:  
\[\begin{split}\min_{u_{k},\omega_{k}}&\left\{\sum_{k=1} ^{K}\sum_{m=1}^{M}\|\partial_{t}\Big{[}\Big{(}\delta(t)+\frac{i}{\pi t}\Big{)} \times u_{k,m}(t)\Big{]}e^{-i\omega_{k}t}\|_{2}^{2}\right\},\\ & s.t.\quad\sum_{k=1}^{K}u_{k,m}(t)=x_{i}(t),i=1,2,\ldots,M.\end{split} \tag{9}\]  
In Equation (9), the target function becomes the sum of the bandwidth of all of the modes over all of the channels, and the central frequencies \(\omega_{k}\) in each mode are the same among different channels. The target function in Equation (8) can be reformulated as an augmented Lagrangian:  
\[\begin{split}& L(u_{k},\omega_{k},\lambda)=\alpha\sum_{k=1}^{K} \sum_{m=1}^{M}\|\partial_{t}\Big{[}\Big{(}\delta(t)+\frac{i}{\pi t}\Big{)} \times u_{k,m}(t)\Big{]}e^{-i\omega_{k}t}\|_{2}^{2}\\ &+\sum_{m=1}^{M}\|x_{m}(t)-\sum_{k=1}^{K}u_{k,m}(t)\|_{2}^{2}+ \sum_{m=1}^{M}\left\langle\lambda_{m}(t),x_{m}(t)-\sum_{k=1}^{K}u_{k,m}(t) \right\rangle\.\end{split} \tag{10}\]  
where the Lagrangian multiplier and the penalty weighting factor play the same role as those in VMD. The ADMM procedure is also used for the solution of Equation (10), as formulated in Equations (10)-(12), which is slightly different from that in VMD:  
\[\begin{split}\hat{a}_{k,m}^{n+1}&=\frac{\hat{x}_{m} -\sum_{j\
eq k}\hat{u}_{j,m}+\big{(}\hat{\lambda}_{l,m}/2\big{)}}{1+2\alpha( \omega-\omega_{k})^{2}},\end{split} \tag{11}\]  
\[\begin{split}\omega_{k}^{n+1}&=\frac{\sum_{m=1}^{M }\int_{0}^{\infty}\omega\big{|}\hat{a}_{k,m}(\omega)\big{|}^{2}d\omega}{\sum_{ m=1}^{M}\int_{0}^{\infty}\big{|}\hat{a}_{k,m}(\omega)\big{|}^{2}d\omega}, \end{split} \tag{12}\]  
\[\begin{split}\lambda^{n+1}&=\lambda^{n}+\tau\Bigg{(} x_{m}(\omega)-\sum_{k}u_{k,m}^{n+1}\Bigg{)},m=1:M\end{split} \tag{13}\]  
Until convergence:  
\[\begin{split}\sum_{m}\sum_{k}\|\hat{a}_{k,m}^{n+1}-\hat{a}_{k,m}^ {n}\|_{2}^{2}/\|\hat{a}_{k,m}^{n}\|_{2}^{2}<\varepsilon.\end{split} \tag{14}\]  
The associated explanation can be located in the original work of [PERSON] and [PERSON] [24] and is not repeated here.  
Figure 1 gives a schematic illustration of the algorithm of MVMD for a 2D seismic profile. In the algorithm, MVMD unifies the quantity of the decomposed component and their center frequency. It can not only increase the seismic signal similarity, but also the calculation efficiency by exchanging the time through memory space. The decomposition results reflect the seismic geological structures on various scales.