# Robust 3D Joint Inversion of Gravity and Magnetic Data: A High-Performance Computing Approach
## 2 Methods
### Regularized Least Squares 3D Joint Inversion with Gramian-Based Constraint  
The following parametric functional expressed in its matrix form was used to make joint inversions of two types of geophysical data:  
\[ \begin{split} P(\
icefrac{{\alpha^{(1)}}}{{\alpha^{(2)}}}, \
icefrac{{\beta^{(1)}}}{{\beta^{(2)}}},\
icefrac{{\gamma^{(1)}}}{{\gamma^{(2) }}},\
icefrac{{\eta^{(1)}}}{{\eta^{(2)}}},\
icefrac{{\mu^{(1)}}}{{\mu^{(2)}}} )(m^{(1)},\ m^{(2)})\\ =\|W_{d}[A(m)-d_{0}]\|_{2}^{2}+\alpha\|L(m)\|_{2}^{2}+\beta\|W_{m} \big{(}m-m_{apr}\big{)}\|_{2}^{2}\\ +\gamma\|M(m)\|_{2}^{2}+\eta\|N(m)\|_{2}^{2}+\mu G(\
abla m, \
abla m)\end{split}  \tag{2}\]  
This is composed of the mismatch functional \(\|W_{d}[A(m)-d_{0}]\|_{2}^{2}\), a smoothing regularizer \(\|L(m)\|_{2}^{2}\), an _a priori_ model approximation regularizer \(\|W_{m}\big{(}m-m_{apr}\big{)}\|_{2}^{2}\), a structural direction regularizer \(\|M(m)\|_{2}^{2}\), a verticality regularizer \(\|N(m)\|_{2}^{2}\) and the structural coupling regularizer using Gramian \(G(\
abla m,\
abla m)\). Each of the elements represents the concatenated vectors and matrices of each input data type; in this case, gravimetric (1) and magnetic (2) data. \[m= \begin{bmatrix}m^{(1)}\\ m^{(2)}\end{bmatrix} ;m_{apr}= \begin{bmatrix}m^{(1)}_{apr}\\ m^{(2)}_{apr}\end{bmatrix} ;d= \begin{bmatrix}d^{(1)}_{0}\\ d^{(2)}_{0}\end{bmatrix} ;A(m)= \begin{bmatrix}A^{(1)}(m^{(1)})\\ A^{(2)}(m^{(2)})\end{bmatrix} ; \tag{3}\] \[G(\
abla m,\
abla\mathbf{m})= \begin{bmatrix}G\
abla m^{(1)}, \
abla m^{(2)}\\ G\
abla m^{(2)},\
abla m^{(1)}\end{bmatrix} ;\] \[A= \begin{bmatrix}A^{(1)}&\vdots\\ \vdots&A^{(2)}\end{bmatrix} ;\;W_{d}= \begin{bmatrix}W^{(1)}_{d}&\vdots\\ \vdots&W^{(2)}_{d}\end{bmatrix} ;W_{m}= \begin{bmatrix}W^{(1)}_{m}&\vdots\\ \vdots&W^{(2)}_{m}\end{bmatrix} ;\] \[\alpha= \begin{bmatrix}\alpha^{(1)}I&\vdots\\ \vdots&\alpha^{(2)}I\end{bmatrix} ;\;\beta= \begin{bmatrix}\beta^{(1)}I&\vdots \\ \vdots&\beta^{(2)}I\end{bmatrix} ;\gamma= \begin{bmatrix}\gamma^{(1)}I&\vdots\\ \vdots&\gamma^{(2)}I\end{bmatrix} ;\] \[\eta= \begin{bmatrix}\eta^{(1)}I&\vdots\\ \vdots&\eta^{(2)}I\end{bmatrix} ;\;\;\mu= \begin{bmatrix}\mu^{(1)}I&\vdots\\ \vdots&\mu^{(2)}I\end{bmatrix} \]  
where \(I\) is an identity matrix, \(A\) is the linear operator for the direct modeling of gravimetry and magnetometry respectively, \(m\) is the vector encoding the model parameters or physical property distribution of a 3D domain (subsurface), \(d_{0}\) is the vector containing the measured data in a 2D domain (surface), \(m_{apr}\) is the vector containing _a priori_ information about the distribution of properties in the subsurface, \(W_{d}\) is the weights diagonal matrix of inverse standard deviation assigned to each data and \(W_{m}\) is the weights diagonal matrix of inverse standard deviation assigned to each cell in the model. The regularization parameters \(\alpha^{(1)}\), \(\alpha^{(2)}\), \(\beta^{(1)}\), \(\beta^{(2)}\), \(\gamma^{(1)}\), \(\gamma^{(2)}\), \(\eta^{(1)}\), \(\eta^{(2)}\), \(\mu^{(1)}\), \(\mu^{(2)}\) determine the weight of the regularization functionals, which allows varying importance to be given to the regularizers in the inversion solutions for each specific case, thus having more control over the geometries of the contrasts obtained in the inversion models, thereby leaving the freedom to choose what kind of results are expected to the user's objectivity.