3D Focusing Inversion of Full Tensor Magnetic Gradiometry Data with Gramian Regularization

Abstract

Full tensor magnetic gradiometry (FTMG) is becoming a practical method for exploration due to recent advancements in superconducting quantum interference device (SQUID) technology. This paper introduces an efficient method of 3D modeling and inversion of FTMG data. The forward modeling uses single-point Gaussian integration with pulse basis functions to compute the volume integrals representing the second spatial derivatives of the magnetic potential. The inversion is aimed at recovering both the magnetic susceptibility and magnetization vectors. We have introduced a 3D regularized focusing inversion technique that utilizes Gramian regularization and a moving sensitivity domain approach. We have also developed a new method of magnetization vector decomposition into induced and remanent parts. The case study includes applying the developed inversion method and computer code to interpret a helicopter-borne FTMG survey carried out over the Thompson Nickel Belt. We have analyzed and separately inverted the observed FTMG and total magnetic intensity (TMI) data using the developed 3D inversion methods to obtain the subsurface susceptibility and magnetization vector models. Furthermore, we present a comparison of the inversions utilizing the FTMG data and the TMI data.

1. Introduction

The measurement of magnetic vector data using orthogonal fluxgate magnetometers is heavily influenced by the Earth’s magnetic field and can be greatly affected by the orientation of the instruments. Therefore, due to the instability of airborne platforms, cesium vapor magnetometers have been preferred as they directly measure the total magnetic intensity (TMI) and are not affected by instrument orientation.

Despite this, direct measurements of magnetic tensors have advantages, as they contain directional sensitivity, which is important for determining rock magnetization vector and studying remanent magnetization. Recently, full tensor magnetic gradiometers based on superconducting quantum interference devices (SQUIDs) have been developed and are now commercially used for geophysical surveying (e.g., [1,2,3,4]). The quantum magnetic tensor (QMAGT) system, developed by Dias Airborne of Saskatoon, Saskatchewan, Canada in collaboration with Supracon AG of Jena, Thuringia, Germany, is a helicopter-borne magnetic survey system that measures the FTMG data using the SQUID, operating within a liquid helium bath (https://www.diasgeo.com/technology-innovation/full-tensor-magnetic-gradiometry-ftmg/, accessed on 21 June 2023).

Inverting the data obtained from full tensor magnetic gradiometry into a magnetization vector through 3D inversion remains a challenging problem [4,5]. In this paper, we present recent advances in FTMG data modeling and inversion and illustrate the developed methods using the results of a helicopter-borne QMAGT FTMG survey conducted over the Thompson Nickel Belt near Thompson, Manitoba, Canada.

We also demonstrate that the inversion of FTMG data has several advantages over traditional TMI data for mineral exploration. TMI data usually require significant post-processing before inversion, such as arbitrary filtering, but FTMG data are naturally high-frequency and require only minimal post-processing, such as leveling and outlier removal, before inversion. Furthermore, the FTMG data are less sensitive to instrument orientation than TMI data because magnetic gradients are mainly caused by localized sources rather than the Earth’s background field or regional trends. However, the depth sensitivity of FTMG data is limited, typically to about 600 m below the surface, whereas TMI data can penetrate several kilometers below the surface.

Several publications discussed the methods of the FTMG data modeling and inversion (e.g., [4,6,7,8,9]). For completeness, we present below the basic principles forming the foundations of these methods.

2. Modeling of FTMG Data for the Models without Remanent Magnetization

In what follows, we assume first that there is no remanent magnetization, the impact of self-demagnetization is insignificant, and the magnetic susceptibility is uniform in all directions. Given these assumptions, the magnetization intensity I(r) can be directly linked to the inducing magnetic field, H⁰(r), through the magnetic susceptibility, χ(r):

$$I\left(r\right)=\chi \left(r\right)H^{\mathbf{0}}\left(r\right).$$

(1)

We discretize the 3D earth model into a grid of N_m cells, each of constant magnetic susceptibility. Following [10], the magnetic potential can be expressed in discrete form as follows:

$$U\left({\mathbf{r}}^{\prime }\right)={\sum }_{k=1}^{N_m}{\chi }_k{{\displaystyle \iiint }}_{D_k}^{}\psi \left({\mathbf{r}}^{\prime },r\right)dv,$$

(2)

where $\psi \left({\mathbf{r}}^{\prime },r\right)$ is the whole-space Green’s function for the magnetic potential:

$$\psi \left({\mathbf{r}}^{\prime },r\right)={\nabla }^{\prime }\frac{1}{\left|r-{\mathbf{r}}^{\prime }\right|}.$$

As we will show, all components of the magnetic field can be computed as the spatial derivatives of Equation (2). For example, the magnetic field, H(r), is the first derivative of the magnetic potential:

$$H\left(r\right)=\nabla U\left(r\right).$$

(3)

Substituting Equation (2) into (3), we arrive at the following discrete form:

$$H\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}{\chi }_k{{\displaystyle \iiint }}_{D_k}^{}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [l-\frac{3\left(l·\left(r-{\mathbf{r}}^{\prime }\right)\right)\left(r-{\mathbf{r}}^{\prime }\right)}{{\left|r-{\mathbf{r}}^{\prime }\right|}^2}\right]dv,$$

(4)

where ${\mathbf{r}}^{\prime }=\left(x^{\prime },y^{\prime },z^{\prime }\right)$ denotes the point of observation; $r=\left(x,y,z\right)$ denotes the point of source location; and $l=\left(l_x,l_y,l_z\right)$ and $H_0$ are the direction of magnetization and the absolute value of the inducing magnetic field, H⁰, respectively.

In our chosen coordinate system, we consider positive y as the north direction, positive x as the east direction, and positive z as the downward direction. When conducting a magnetic survey, the inclination (I), declination (D), and azimuth (A) of the inducing magnetic field (measured in degrees) can be determined by referencing the International Geomagnetic Reference Field (IGRF) model. Under these assumptions, the components of the unit vector in the direction of the inducing magnetic field are as follows:

$$l_x=cos\left(I\right)sin\left(D-A\right),$$

(5)

$$l_y=cos\left(I\right)cos\left(D-A\right),$$

(6)

$$l_z=sin\left(I\right).$$

(7)

Previous studies (e.g., [11]) have provided analytical solutions for calculating the volume integral described in Equation (4) for right rectangular prisms with magnetic susceptibility. However, we choose to perform the volume integral computation numerically using single-point Gaussian integration with pulse basis functions, as proposed by [12]. This method dramatically speeds up the processing time compared to the conventional approach based on prismatic cell approximation while yielding a very accurate result. In this case, $r=\left(x,y,z\right)$ denotes the cell center. We assume constant discretization of ∆x, ∆y, and ∆z in the x, y, and z directions, respectively.

Thus, Equation (4) can be simplified as follows:

$$H\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}{\chi }_k\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [l-\frac{3\left(l·\left(r-{\mathbf{r}}^{\prime }\right)\right)\left(r-{\mathbf{r}}^{\prime }\right)}{{\left|r-{\mathbf{r}}^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z,$$

(8)

From Equation (8), we can derive discrete expressions for the scalar components of the magnetic field:

$$H_x\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}{\chi }_k\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [l_x-\frac{3t\left(x_k-x^{\prime }\right)}{{\left|r-{\mathbf{r}}^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z,$$

(9)

$$H_y\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}{\chi }_k\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [l_y-\frac{3t\left(y_k-y^{\prime }\right)}{{\left|r-{\mathbf{r}}^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z,$$

(10)

$$H_z\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}{\chi }_k\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [l_z-\frac{3t\left(z_k-z^{\prime }\right)}{{\left|r-{\mathbf{r}}^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z,$$

(11)

where

$$t=l_x\left(x_k-x^{\prime }\right)+l_y\left(x_k-x^{\prime }\right)+l_z\left(x_k-x^{\prime }\right).$$

(12)

In the case of studying the full magnetic tensor, one has to calculate the second derivatives of the magnetic potential as follows:

$$H_{\alpha \beta }\left(r\right)=\frac{{\partial }^2}{\partial \alpha \partial \beta }U\left(r\right), \alpha ,\beta =x,y,z.$$

(13)

The second spatial derivatives of the magnetic potential form a symmetric tensor with zero trace:

$$\widehat{\mathbf{H}}=\left. [\begin{array}{ccc}{H}_{xx}& H_{xy}& H_{xz}\\ H_{yx}& H_{yy}& H_{yz}\\ H_{zx}& H_{zy}& H_{zz}\end{array}\right], H_{xx}+H_{yy}+H_{zz}=0,$$

(14)

where

$$H_{\alpha \beta }\left(r\right)=\frac{\partial H_{\alpha }}{\partial \beta }, \alpha ,\beta =x,y,z.$$

(15)

This implies that of the nine tensor components, only five are independent.

After some algebra, we find discrete forms for each component of the magnetic tensor [13]:

$$H_{xx}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{\left. [-l_x\left(x_k-x^{\prime }\right)-t\right]r^2+2{\left(x^{\prime }-x_k\right)}^{2}t}{r^7}\right.+\left. [l_x-\frac{3t\left(x_k-x^{\prime }\right)}{r^2}\right]\left.\frac{\left(x^{\prime }-x_k\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(16)

$$H_{xy}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{-l_y\left(x_k-x^{\prime }\right)r^2+2\left(y^{\prime }-y_k\right)\left(x^{\prime }-x_k\right)t}{r^7}\right.+\left. [l_x-\frac{3t\left(x_k-x^{\prime }\right)}{r^2}\right]\left.\frac{\left(y^{\prime }-y\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(17)

$$H_{xz}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{-l_z\left(x_k-x^{\prime }\right)r^2+2\left(z^{\prime }-z_k\right)\left(x^{\prime }-x_k\right)t}{r^7}\right.+\left. [l_x-\frac{3t\left(x_k-x^{\prime }\right)}{r^2}\right]\left.\frac{\left(z^{\prime }-z_k\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(18)

$$H_{yx}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{\left. [-l_x\left(y_k-y^{\prime }\right)-t\right]r^2+2\left(x^{\prime }-x_k\right)\left(y^{\prime }-y\right)t}{r^7}\right.+\left. [l_y-\frac{3t\left(y_k-y^{\prime }\right)}{r^2}\right]\left.\frac{\left(x^{\prime }-x_k\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(19)

$$H_{yy}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{\left. [-l_y\left(y_k-y^{\prime }\right)-t\right]r^2+2{\left(y^{\prime }-y_k\right)}^{2}t}{r^7}\right.+\left. [l_y-\frac{3t\left(y_k-y^{\prime }\right)}{r^2}\right]\left.\frac{\left(y^{\prime }-y_k\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(20)

$$H_{yz}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{-l_z\left(y_k-y^{\prime }\right)r^2+2\left(z^{\prime }-z_k\right)\left(y^{\prime }-y_k\right)t}{r^7}\right.+\left. [l_y-\frac{3t\left(y_k-y^{\prime }\right)}{r^2}\right]\left.\frac{\left(z^{\prime }-z_k\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(21)

$$H_{zx}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{-l_x\left(x_k-x^{\prime }\right)r^2+2\left(x^{\prime }-x_k\right)\left(z^{\prime }-z_k\right)t}{r^7}\right.+\left. [l_z-\frac{3t\left(z_k-z^{\prime }\right)}{r^2}\right]\left.\frac{\left(x^{\prime }-x_k\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(22)

$$H_{zy}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{-l_y\left(x_k-x^{\prime }\right)r^2+2\left(y^{\prime }-y_k\right)\left(z^{\prime }-z_k\right)t}{r^7}\right.+\left. [l_z-\frac{3t\left(z_k-z^{\prime }\right)}{r^2}\right]\left.\frac{\left(y^{\prime }-y_k\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(23)

$$H_{zz}\left({\mathbf{r}}^{\prime }\right)=3H_0{\sum }_{k=1}^{N_m}{\chi }_k\left\{\frac{\left. [-l_z\left(z_k-z^{\prime }\right)-t\right]r^2+2{\left(z^{\prime }-z_k\right)}^{2}t}{r^7}\right.+\left. [l_z-\frac{3t\left(z_k-z^{\prime }\right)}{r^2}\right]\left.\frac{\left(z^{\prime }-z_k\right)}{r^5}\right\}\Delta x\Delta y\Delta z,$$

(24)

where t is defined by Equation (12), and

$$r=\sqrt{{\left(x_k-x^{\prime }\right)}^2+{\left(y_k-y^{\prime }\right)}^2+{\left(z_k-z^{\prime }\right)}^2}.$$

Equations (16) through (24) provide the basis for computing the FTMG data for the models without the remanent magnetization.

3. Magnetization Vector Decomposition into Induced and Remanent Parts

Conventional inversions for magnetic susceptibility only consider induced magnetization, which limits their ability to provide a comprehensive geological understanding [4,14,15]. In order to obtain an adequate understanding of geology and mineralization and recover information about remanent magnetization, magnetic data must be inverted towards a magnetization vector model. This paper presents a new method for determining remanent magnetization using the magnetization vector decomposition technique. This method decomposes the full magnetization vector into the induced (inline relative to the inducing) and remanent components (Figure 1). The amplitude of the inline magnetization can provide models similar to susceptibility since paramagnetic and ferromagnetic materials tend to align parallel to the inducing field, while diamagnetic materials tend to align in the opposite direction. The remanent magnetization can manifest itself as a vector pointing away from the inducing field (Figure 1). We illustrate this novel method for direct imaging of remanent magnetization, which is useful for determining optimal exploration strategies.

Remanent magnetization (or remanence) is a permanent magnetization of a rock obtained in the past when the Earth’s magnetic field had a different magnitude and direction than it has today. It follows that the total intensity of magnetization, I(r), is linearly related to both the induced, M_ind, and remanent, M_rem, magnetizations:

$$I\left(r\right)=H_0\left. [M_{ind}\left(r\right)+M_{rem}\left(r\right)\right],$$

(25)

where the induced magnetization vector, M_ind, is linearly proportional to the inducing magnetic field, H⁰(r), through the magnetic susceptibility, χ(r):

$$M_{ind}\left(r\right)=\chi \left(r\right)H^0\left(r\right)/H_0=\chi \left(r\right)l,$$

(26)

and l the unit vector in the direction of the inducing magnetic field, defined by Equations (5)–(7). We should note that we have defined the magnetization vectors (both the induced, M_ind, and remanent, $M_{rem}$) as unitless for convenience of derivations.

The Koenigsberger ratio, Q, is the ratio of the absolute values of the remanent magnetization to the induced magnetization [16]:

$$Q=\frac{\left|M_{rem}\right|}{\left|M_{ind}\right|}.$$

(27)

For Koenigsberger ratios greater than 1, the remanent magnetization vector is the predominant contribution to the total intensity of magnetization [17].

We can rewrite Equation (25) as follows:

$$I\left(r\right)=H_0\mathit{M},$$

(28)

where M is the dimensionless magnetization vector:

$$M\left(r\right)=M_{ind}\left(r\right)+M_{rem}\left(r\right).$$

(29)

At the same time, one can also consider a simplified alternative method to decompose the full magnetization vector into the induced and remanent components. This involves finding the magnitude of the projection of the vector onto the inducing field direction. This is termed the inline component, which is a proxy for induced magnetization. We then calculate the magnitude of vector elements that are not parallel, which we term the perpendicular component. These elements are not necessarily perpendicular, but they are not parallel to the inducing field, and this “perpendicular” term is the proxy for remanent magnetization.

4. Modeling of FTMG Data for the Models with Remanent Magnetization

For the modeling of magnetic data, we discretize the 3D earth model into a grid of N_m cells, each of a constant magnetization vector. Following [10], the magnetic field can be expressed in discrete form as follows:

$$H\left({\mathbf{r}}^{\prime }\right)=-H_{\mathbf{0}}{\sum }_{k=1}^{N_m}{{\displaystyle \iiint }}_{V_k}^{}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [M_k-\frac{3\left(M_k·\left(r-r^{\prime }\right)\right)\left(r-{\mathbf{r}}^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]dv,$$

(30)

where $M_k=\left(M_{xk},M_{yk},M_{zk}\right)$ is the magnetization vector of the k-th cell.

Closed-form solutions for the volume integral in Equation (30) over right rectangular magnetic prisms have been previously presented (e.g., [11]). As discussed in [12], we prefer to evaluate the volume integral numerically using single-point Gaussian integration with pulse basis functions because this numerical solution is as accurate as the analytic solution, provided that the depth to the center of the cell exceeds twice the dimension of the cell. In this case, $r=\left(x,y,z\right)$ now denotes the center of the cell. We assume constant discretization of ∆x, ∆y, and ∆z in the x, y, and z directions, respectively. It follows that Equation (30) can be simplified as follows:

$$H\left({\mathbf{r}}^{\prime }\right)=-H_{\mathbf{0}}{\sum }_{k=1}^{N_m}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [M_k-\frac{3\left(M_k·\left(r-r^{\prime }\right)\right)\left(r-{\mathbf{r}}^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z.$$

(31)

From Equation (31), we obtain discrete expressions for the vector components of the magnetic field:

$$H_x\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [M_{xk}-\frac{3t_k\left(x_{\mathit{k}}-x^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z,$$

(32)

$$H_y\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [M_{yk}-\frac{3t_k\left(y_{\mathit{k}}-y^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z,$$

(33)

$$H_z\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [M_{zk}-\frac{3t_k\left(z_{\mathit{k}}-z^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z,$$

(34)

where

$$t_k=M_{xk}\left(x_k-x^{\prime }\right)+M_{yk}\left(y_k-y^{\prime }\right)+M_{zk}\left(z_k-z^{\prime }\right).$$

(35)

We can write Equations (32)–(34) in a compact form as follows:

$$H_{\alpha }\left({\mathbf{r}}^{\prime }\right)=-H_0{\sum }_{k=1}^{N_m}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [M_{\alpha k}-\frac{3t_k\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z, \alpha =x,y,z.$$

(36)

The second spatial derivatives of the magnetic potential,

$$H_{\alpha \beta }\left(r\right)=\frac{{\partial }^2}{\partial \alpha \partial \beta }U\left(r\right), \alpha ,\beta =x,y,z,$$

(37)

form a symmetric magnetic tensor

$$\widehat{\mathbf{H}}=\left. [\begin{array}{ccc}{H}_{xx}& H_{xy}& H_{xz}\\ H_{yx}& H_{yy}& H_{yz}\\ H_{zx}& H_{zy}& H_{zz}\end{array}\right]$$

(38)

with zero trace

$$H_{xx}+H_{yy}+H_{zz}=0,$$

(39)

where

$$H_{\alpha \beta }=\frac{\partial H_{\alpha }}{\partial \beta }, \alpha ,\beta =x,y,z,$$

(40)

As for the case without the remanent magnetization considered above, this implies that of the nine tensor components, only five are independent.

By introducing $r_k=r=\sqrt{{\left(x_k-x^{\prime }\right)}^2+{\left(y_k-y^{\prime }\right)}^2+{\left(z_k-z^{\prime }\right)}^2}$ and after some algebra, we find discrete forms for each component of the magnetic tensor:

$$H_{\mathsf{\alpha }\mathsf{\beta }}\left({\mathbf{r}}^{\prime }\right)=\frac{\partial H_{\alpha }\left(r^{\prime }\right)}{\partial {\beta }^{\prime }}=-H_0\frac{\partial }{\partial {\beta }^{\prime }}\sum _{k=1}^{N_m}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^3}\left. [M_{\alpha k}-\frac{3t_k\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]\Delta x\Delta y\Delta z.$$

(41)

We now take the derivatives:

$$\frac{\partial }{\partial {\beta }^{\prime }}\left(\frac{1}{{\left|r-r^{\prime }\right|}^3}\left. [M_{\alpha k}-\frac{3t_k\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]\right)=\frac{3}{{\left|r-{\mathbf{r}}^{\prime }\right|}^5}\left. [\left({\beta }_{\mathit{k}}-{\beta }^{\prime }\right)M_{\alpha k}-\frac{5t_k\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right)\left({\beta }_{\mathit{k}}-{\beta }^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}+M_{\beta k}\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right)+t_k{\delta }_{\alpha \beta }\right],$$

(42)

where

$${\delta }_{\alpha \beta }=\left\{\begin{array}{c}1, \alpha =\beta \\ 0, \alpha \ne \beta \end{array}\right.$$

Substituting Equation (42) into (41), we obtain

$$H_{\mathsf{\alpha }\mathsf{\beta }}\left({\mathbf{r}}^{\prime }\right)=-3H_0{\sum }_{k=1}^{N_m}\frac{1}{{\left|r-{\mathbf{r}}^{\prime }\right|}^5}\left. [\left({\beta }_{\mathit{k}}-{\beta }^{\prime }\right)M_{\alpha k}-\frac{5t_k\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right)\left({\beta }_{\mathit{k}}-{\beta }^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}+M_{\beta k}\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right)+t_k{\delta }_{\alpha \beta }\right]\Delta x\Delta y\Delta z.$$

(43)

According to (35), we can write

$$t_k={{\displaystyle \sum }}_{\gamma =x,y,z}^{}\left({\gamma }_{\mathit{k}}-{\gamma }^{\prime }\right)M_{\gamma k}.$$

(44)

and

$$M_{\alpha k}={{\displaystyle \sum }}_{\gamma =x,y,z}^{}{\delta }_{\alpha \gamma }{M}_{\gamma k},M_{\beta k}={{\displaystyle \sum }}_{\gamma =x,y,z}^{}{\delta }_{\beta \gamma }{M}_{\gamma k}.$$

(45)

Substituting (44) and (45) into (43), we have

$$H_{\mathsf{\alpha }\mathsf{\beta }}\left({\mathbf{r}}^{\prime }\right)=-3H_0{\sum }_{k=1}^{N_m}\frac{\Delta x\Delta y\Delta z}{{\left|r-{\mathbf{r}}^{\prime }\right|}^5}{{\displaystyle \sum }}_{\gamma =x,y,z}^{}\left. [\left({\beta }_{\mathit{k}}-{\beta }^{\prime }\right){\delta }_{\alpha k}+\left. [{\delta }_{\alpha \beta }-\frac{5\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right)\left({\beta }_{\mathit{k}}-{\beta }^{\prime }\right)}{{\left|r-r^{\prime }\right|}^2}\right]\left({\gamma }_{\mathit{k}}-{\gamma }^{\prime }\right)+\left({\alpha }_{\mathit{k}}-{\alpha }^{\prime }\right){\delta }_{\beta \gamma }\right]M_{\gamma k}.$$

(46)

We introduce a sensitivity kernel for magnetic tensor as follows:

$$G_{\alpha \beta k}^{\gamma }=\left({\beta }_k-{\beta }^{\prime }\right){\delta }_{\alpha \gamma }+\left. [{\delta }_{\alpha \beta }-\frac{5\left({\alpha }_k-{\alpha }^{\prime }\right)\left({\beta }_k-{\beta }^{\prime }\right)}{r-{r^{\prime }}^2}\right]\left({\gamma }_k-{\gamma }^{\prime }\right)+\left({\alpha }_k-{\alpha }^{\prime }\right){\delta }_{\beta \gamma }.$$

(47)

Using this notation, we can rewrite Equation (46) in the following compact form:

$$H_{\alpha \beta }\left({\mathbf{r}}^{\prime }\right)=-3H_0{\sum }_{k=1}^{N_m}\frac{\Delta x\Delta y\Delta z}{{\left|r-{\mathbf{r}}^{\prime }\right|}^5}{{\displaystyle \sum }}_{\gamma =x,y,z}^{}{G}_{\alpha \beta k}^{\gamma }{M}_{\gamma k}, \alpha ,\beta =x,y,z.$$

(48)

where $M_{\gamma k}$ are the components of the magnetization vector:

$$M_k={\left. [M_{xk},M_{yk},M_{zk}\right]}^T.$$

(49)

Equations (36) and (48) are the key equations that we need for solving both modeling and inversion for remanent magnetization. Note also that any implementation of these equations includes both topography and variable flight height, obviating the need for upward continuation of the data.

5. Principles of Inversion

Regardless of the iterative scheme used, most regularized inversions seek to minimize the Tikhonov parametric functional, $P^{\alpha }\left(m\right)$:

$$P^{\alpha }\left(m\right)=\varphi \left(m\right)+\alpha s\left(m\right)\to min,$$

(50)

where $\varphi \left(m\right)$ is a misfit functional of the observed and predicted potential field data, $s\left(m\right)$ is a stabilizing functional, and $\alpha$ is the regularization parameter that balances the misfit and stabilizing functionals [10,18]. Data and model weights can be introduced to Equation (50) through data and model weighting matrices. The purpose of these weighting matrices is to adjust the inverse problem in logarithmic space, thereby decreasing the range of variation in both the data and model parameters.

In our specific implementation, the weights assigned to the data and model are determined by considering their integrated sensitivity, as originally introduced in [19]. The concept of sensitivity weighting in geophysical inverse problems was also discussed in a number of publications, e.g., [20,21] and many others. The interested reader can find a detailed description of the theory of sensitivity weighting applied to general linear and nonlinear inverse problems in the textbooks by Zhdanov [12,18]. The sensitivity-based weighting functions ensure that the various components of the observed data exhibit equal sensitivity to cells situated at different depths and horizontal positions. Consequently, the sensitivity-based weighting functions inherently incorporate necessary adjustments for the vertical and horizontal distribution of susceptibility. This distinction represents a key difference between our approach and the geometric depth weighting functions devised by [22]. Another critical issue of the potential field inversion is the nonuniqueness of the inverse model. This question was addressed in hundreds of published papers and monographs, including those cited above. The most general approach to overcoming inherent nonuniqueness is based on the regularization theory. The nonuniqueness of the inversion is significantly reduced by bringing the a priori information via a regularization. For example, in our paper, we use Gramian and focusing regularization, which forces a level of structural similarity between the magnetization vector components.

Every geological constraint is expressed as a form of regularization, which can be measured by adjusting data weights, establishing upper and lower bounds for the model, determining model weights, incorporating prior knowledge, and selecting an appropriate stabilizing functional. The stabilizing functional incorporates information about the category of models utilized in the inversion process. The selection of a suitable stabilizing functional should rely on the user’s geological understanding. In the following section, we will provide a concise overview of various smooth and focusing stabilizers, showcasing the outcomes obtained from the 3D inversion of magnetic vector and tensor data for each approach. A minimum norm (MN) stabilizer will seek to minimize the norm of the difference between the current model and an a priori model:

$$s_{MN}\left(m\right)={{\displaystyle \iiint }}_V^{}{\left(m-m_{apr}\right)}^{2}dv,$$

(51)

and usually produces a relatively smooth model (where m = m(r) is some function describing the model parameters).

The first derivative (FD) stabilizer implicitly introduces smoothness as the first spatial derivatives of the model parameters:

$$s_{FD}\left(m\right)={{\displaystyle \iiint }}_V^{}{\left(\nabla m-\nabla m_{apr}\right)}^{2}dv,$$

(52)

and can result in spurious oscillations and artifacts when the model parameters are discontinuous. A combination of stabilizers (51) and (52) is often used (e.g., [22]).

In reality, the majority of geological formations do not possess smooth magnetization distributions. Geology often displays distinct boundaries with sharp contrasts in susceptibility, such as those observed between an ore deposit and its surrounding rock, or at a discontinuity. Consequently, the application of stabilizers (51) and (52), or their combinations, leads to results that lack physical relevance to the actual geological conditions. To address this limitation, ref. [23] introduced focusing stabilizers that enable the recovery of models with sharper boundaries and greater contrast. In this paper, we consider the minimum support (MS) stabilizer:

$$s_{MS}\left(m\right)={{\displaystyle \iiint }}_V^{}\frac{{\left(m-m_{apr}\right)}^2}{{\left(m-m_{apr}\right)}^2+e^2}dv,$$

(53)

where e is a focusing parameter introduced to avoid singularity when $m=m_{apr}$.

The minimum support stabilizer aims to minimize the volume that deviates from the a priori model, thereby promoting the recovery of compact bodies. Consequently, any smooth distribution of model parameters with only minor deviations from the a priori model is penalized.

While variations of Equation (53) were derived in [10,12], we base our solution on the re-weighted regularized conjugate gradient method (RRCG) [10,18], which is easier to implement numerically.

In numerical dressing, the corresponding parametric functional can be written as follows:

$$P^{\alpha }\left(m,d\right)={\left(W_{d}Am-W_{d}d\right)}^T\left(W_{d}Am-W_{d}d\right)+\alpha {\left(W_e{W}_{m}m-W_e{W}_m{m}_{apr}\right)}^T\left(W_e{W}_{m}m-W_e{W}_m{m}_{apr}\right),$$

(54)

where W_d is the weighting matrix for the data, and W_m is a diagonal matrix for weighting the model parameters based on integrated sensitivity, ($F^{T}F$), where $F$ is the Frechet derivative of the forward operator:

$$W_m=diag{\left(F^{T}F\right)}^{1/4}.$$

(55)

Matrix W_e is a diagonal matrix, determined by the discrete values of the model parameters, m_i, representing the action of the minimum support stabilizer [10]:

$$W_e=diag\left. [w_e\right]=diag\left. [\frac{1}{{\left(m_i^2+e^2\right)}^{1/2}}\right],$$

(56)

where e is a small number (focusing parameter).

The minimization problem (50) can be reformulated using a space of weighted parameters:

$$m^w=WM=W_e{W}_{m}m,$$

(57)

where $W=W_e{W}_m$.

Equation (54) can be rewritten as follows:

$$P^{\alpha }\left(m,d\right)={\left(A^w{m}^w-d\right)}^T\left(A^w{m}^w-d\right)+\alpha {\left(m^w-m_{apr}^w\right)}^T\left(m^w-m_{apr}^w\right),$$

(58)

where $A^w$ is a new forward operator in the space of weighted parameters, which can be related to the forward operator A in the original space as

$$A^w=A{W}^{-1}=A{W}_m^{-1}{W}_e^{-1}.$$

(59)

The algorithm of the re-weighted regularized conjugate gradient method to solve the minimization of parametric functional (58) is given in Appendix A.

The inversion process continues through iterations until one of the following conditions is met: the residual error reaches a predetermined threshold, the decrease in error between consecutive iterations is smaller than the preset threshold, or the maximum number of iterations is reached. Once the inversion is completed, the quality of the results is assessed based on the data misfit and through visual examination of the model. As we discussed above, in practical applications of the inversion method, some boundary conditions must be imposed on the variations of the model parameters:

$$m_i^{-}\le m_i\le m_i^{+},$$

(60)

where $m_i^{-}$ and $m_i^{+}$ are the lower and upper limits of the model parameter $m_i$. However, during the process of minimization of the Tikhonov parametric functional, we can obtain the values of the model parameters outside the above boundaries. In order to limit the interval of possible values of the inverse problem solution, one can introduce a new model parameter $\tilde{m}$ with the property that the corresponding original model parameter m will always remain within the imposed above boundaries $m_i^{-}$ and $m_i^{+}$. One way to solve this problem is by using the following logarithmic transformation to arrive at new model parameters:

$${\tilde{m}}_i=\ln \frac{m_i-m_i^{-}}{m_i^{+}-m_i}.$$

(61)

The corresponding inverse transformation is given by the following formula:

$$m_i=\ln \frac{m_i^{-}-m_i^{+}\exp \left({\tilde{m}}_i\right)}{1+\exp \left({\tilde{m}}_i\right)}.$$

(62)

It is evident that conducting the inversion process in the logarithmic space of model parameters, denoted as $\tilde{m}$, guarantees that the transformed value of each ${\tilde{m}}_i$, regardless of its magnitude, will always remain within the predefined intervals defined by $m_i^{-}$ and $m_i^{+}$. The details of the RRCG algorithm used to solve the inverse problem in the logarithmic space of model parameters are provided in the Appendix A.

6. Gramian Stabilizer

Inverting for the magnetization vector is a more challenging problem than inverting for scalar magnetic susceptibility, because we have three unknown scalar components of the magnetization vector for every cell. At the same time, there is inherent correlation between the different components of the magnetization vector. The different scalar components have similar spatial variations and represent the same zones of anomalous magnetization. Therefore, it is possible to expect that the different components of the magnetization vector should be mutually correlated [15,24]. It was demonstrated in [25] that one can enforce the correlation between the different model parameters by using the Gramian constraints. Following the cited paper, we have included the Gramian constraint in Equation (50) as follows:

$$P^{\alpha }\left(m\right)=\varphi \left(m\right)+\alpha s\left(m\right)+\beta {{\displaystyle \sum }}_{\gamma =x,y,z}^{}{S}_G\left(m_{\gamma },{\chi }_{eff}\right),$$

(63)

where m is the 3N_m length vector of magnetization vector components; $m_{\gamma }$ is the N_m length vector of the $\gamma$ component of magnetization vector, $\gamma =x,y,z$; ${\chi }_{eff}$ is the N_m length vector of the effective magnetic susceptibility, defined as the magnitude of the magnetization vector,

$${\chi }_{eff}=\sqrt{M_x^2+M_y^2+M_z^2}.$$

(64)

Functional S_G is the Gramian constraint,

$$S_G\left(m_{\gamma },{\chi }_{eff}\right)=\left|\begin{array}{cc}\left(m_{\gamma },m_{\gamma }\right)& \left(m_{\gamma },{\chi }_{eff}\right)\\ \left({\chi }_{eff},m_{\gamma }\right)& \left({\chi }_{eff},{\chi }_{eff}\right)\end{array}\right|,$$

(65)

where ${\left(\ast ,\ast \right)}_{}$ denotes the L₂ inner product operation [18].

Using the Gramian constraint (65), we enhance a direct correlation between the scalar components of the magnetization vector with ${\chi }_{eff}$, which is computed at the previous iteration of an inversion and is updated on every iteration. The advantage of using the Gramian constraint in the form of Equation (65) is that it does not require any a priori information about the magnetization vector (e.g., direction, the relationship between different components, etc.).

The minimization problem in Equation (63) is solved using the re-weighted regularized conjugate gradient (RRCG) method described above.

Finally, the regularized inversion outlined above can be applied to very-large-scale problems by incorporating the concept of the moving sensitivity domain [9,26]. According to this concept, the sensitivity matrix for the entire large-scale model could be constructed via the superposition of sensitivity domains for all observation points. The sensitivities are also computed “on the flight” during the computations without storing them in the computer memory. This approach makes it possible to invert the data collected by surveys of several thousand line km on a modest PC cluster as one inversion run.

For example, in our case study presented below, we employ a 10 km size of the sensitivity domain, which results in the sensitivity matrix taking up less than 100 GB of memory even on a fine grid discretization (25 m × 25 m laterally with a reasonable number of vertical layers—less than 50). Calculation of the sensitivity for this sensitivity domain takes only a few minutes. Calculation of the sensitivity matrix for an entire area—say 100 km to a side with the same discretization—can take hours or days on a single CPU or GPU with ~200 GB memory.

7. Determining the Remanent Magnetization

The separation of the magnetization vector into the parts responsible for remanent and induced polarization is a challenging problem if no a priori information is available. However, this problem can be solved approximately, assuming that the susceptibility only inversion produced the correct magnetic susceptibility distribution in the rock formations.

Taking into account Equations (25) and (28), we can find the vector of remanent magnetization M_rem as follows:

$$M_{rem}\left(r\right)=M\left(r\right)-M_{ind}\left(r\right),$$

(66)

where according to Equation (26),

$$M_{ind}\left(r\right)=\chi \left(r\right)l,$$

(67)

where l is the unit vector in the direction of the inducing magnetic field defined by Equations (5)–(7).

The unknown susceptibility, $\chi \left(r\right)$, can be found by susceptibility only inversion of the observed data. On the next step, we invert the data for the magnetization vector, $M\left(r\right)$. Finally, we can find the remanent magnetization by substituting Equation (67) into (66):

$$M_{rem}\left(r\right)=M\left(r\right)-\chi \left(r\right)l.$$

(68)

We should emphasize that the value of remanent magnetization determined by Formula (68) is correct if we use the true values of magnetic susceptibility $\chi \left(r\right)$. In the case when $\chi \left(r\right)$ is found by susceptibility only inversion, Formula (68) can be used as some approximation for the remanent magnetization.

Based on Equation (68), one can determine the Koenigsberger ratio, Q:

$$Q\approx \frac{\left|M_{rem}\right|}{\left|M_{ind}\right|}=\frac{\left|M\left(r\right)-\chi \left(r\right)l\right|}{\left|\chi \left(r\right)l\right|}=\left|\frac{M\left(r\right)}{\chi \left(r\right)}-l\right|.$$

(69)

If we assume that the earth model has the susceptibility of free space, ${\chi }_0$, the vector of induced magnetization is constant,

$$M_{ind}\left(r\right)={\chi }_{0}l,$$

(70)

Therefore, we arrive at the following equation for the remanent magnetization:

$$M_{rem}\left(r\right)=M\left(r\right)-{\chi }_{0}l,$$

(71)

We can approximately calculate the Koenigsberger ratio, Q, as follows:

$$Q\approx \left|\frac{M\left(r\right)}{{\chi }_0}-l\right|.$$

(72)

8. Model Study: Inversion of the TMI and FTMG Data for Synthetic Model

We first conducted a synthetic model study using the developed method and computer code. Figure 2 represents the vertical sections of the magnetized target and the inversion results. The anomaly in red in panel (A) has a susceptibility of 0.1 SI and amplitude of dimensionless magnetization vector 0.1 in a homogeneous background of zero elsewhere. Magnetization vector anomaly has a magnitude of 0.1 in the direction shown by the white arrows. The inclination is 75 degrees, and the declination is 35 degrees. The inducing magnetic field (indicated by the black arrow at the top of Figure 2) has an inclination of 45 degrees and a declination of 5 degrees. This is 30 degrees offset from the direction of the magnetization vector in both inclination and declination, which indicates the presence of “remanent” magnetization in this model. We have calculated the synthetic observed TMI and FTMG data using the vector magnetization model described above.

We ran inversions of the synthetic data for the magnetic susceptibility model first. Panel (B) shows the susceptibility model obtained from the inversion of the synthetic TMI data. Panel (C) shows the susceptibility model obtained from the inversion of synthetic FTMG data. One can see that TMI data inversion has difficulties recovering the correct susceptibility model; the image is diffused and slightly shifted down. At the same time, the FTMG inversion solves this problem more accurately (Panel (C)).

In the next step of this numerical experiment, we run inversion for the magnetization vector model. Panel (D) shows the magnetization vector model obtained from the inversion of synthetic FTMG data. The white arrows show the full magnetization vector, while the color map represents the amplitude of this vector. We can see the reasonable correlation between the inverse magnetization vector and the actual magnetization vector model shown in Panel (A).

Finally, Panels (E) and (F) show the color maps of the inline and “remanent” magnetization amplitudes, respectively. The stronger values in Panel (F) indicate that the anomalous body possesses “remanent” magnetization, responsible for a deviation of the magnetization vector from the inducing magnetic field indicated by the black arrow at the top of Figure 2.

9. Case Study: Inversion of the TMI and FTMG Data in the Thompson Nickel Belt in Northern Manitoba, Canada

Dias Airborne collected approximately 10,655 line-km of QMAGT FTMG data in the Thompson Nickel Belt in northern Manitoba, Canada. The survey aimed to determine the gross geomorphology of the Thompson Nickel Belt which is part of the fifth-largest nickel camp in the world and the associated Ni sulfide mineralization [27,28]. This paper focuses on a small area of interest (AOI) in the northern extent of Thompson Nickel Belt near Mystery Lake, covered by approximately 100 line km. The line spacing was 50 m, with sampling every ~2.5 m, and the flight height was around 40 m. The digital terrain model (DTM) and flight lines over the AOI are shown in Figure 3.

9.1. Susceptibility and Magnetization Vector Models Produced by Inversion of FTMG Data

We conducted separate inversions of the FTMG data to obtain models for both susceptibility and magnetization vector. The inversions used a 20 m × 40 m grid with 48 vertical layers ranging from 10–95 m, but the models were ultimately limited to a depth of 600 m based on the sensitivity study. The Tikhonov regularization method with minimum support stabilizer was employed to produce focused images of magnetized structures [10,18,23,29]. The focusing regularization in Equation (53) was used in the inversions with a focusing epsilon of 0.001 which was not varied during the inversion.

Figure 4 displays observed and predicted data for selected FTMG components from the magnetization vector inversion, with the total L2 norm misfit between observed and predicted data for all six magnetic tensor components being less than 10%.

In this study, the inversion parameters for the magnetization vector model were the X, Y, and Z components of the triaxial vector. The inverse models were represented by three characteristics of the full magnetization vector, namely amplitude, inline amplitude (representing induced magnetization), and approximate amplitude of the remanent magnetization. Additionally, the projection of the full magnetization vector was shown as black arrows over the amplitude in the vertical section.

Figure 5, Figure 6, Figure 7 and Figure 8 give a 3D overview of the inverted models. Figure 5 is a 3D isobody view of the inverted susceptibility model with values greater than 0.1 SI shown in green. This isobody is the mystery ultramafic intrusion. The ultramafic body itself is not an economic resource although it does contain some nickel. Figure 6 is a 3D isobody view of the inverted inline magnetization model, which is a proxy for induced magnetization. Values above 0.02 A/m are shown in yellow and represent the same ultramafic intrusion. We see the enhanced resolution of the body versus the susceptibility model shown in Figure 5. Koenigsberger ratios ($Q$) corresponding to the ultramafic ranged from 1 to 2, indicating a degree of remanence, which is consistent with the difference we see in the susceptibility and inline magnetization models. Figure 7 is a 3D isobody view of the remanent magnetization model, with values over 0.04 A/m shown in red. These red isobodies are potentially breccia-associated mineralization within the Mystery Lake ultramafic coupled with minor sedimentary sulfides which occurs within the package of Pipe Formation rocks. $Q$ values were much higher in the breccia-associated mineralized zone, ranging from 5 to 5.5. These values were a factor in our choice of cutoff values for the isobodies shown. We may also see some circulation of the field in this component due to the approximate decomposition to inline and perpendicular magnetization. This is a powerful result of the MVI versus the susceptibility inversion in that we can image the target mineralization directly with the various products of the MVI. Finally, Figure 8 shows a composite view of the isobodies presented in Figure 5, Figure 6 and Figure 7.

Figure 9 depicts the S–N profiles of the inverted susceptibility and magnetization vector models. The ultramafic intrusion is clearly visible in the susceptibility and inline magnetization models as high values indicated by warm colors. However, in the remanent magnetization model, the intrusion appears as low values, indicated by cool colors, while the country rocks are shown as high values. This technique allows for better identification of potential mineralization targets by imaging the lithologies near the intrusion.

Figure 10 illustrates the W–E profiles of the magnetic models. While the susceptibility model identifies a high susceptibility anomaly near the surface, the magnetization vector models provide a more detailed structure. Specifically, there is a subvertical feature directly west of the susceptibility anomaly at a depth of several hundred meters.

9.2. Comparison of the Inversion Results Using TMI and FTMG Data

Figure 11, Figure 12 and Figure 13 compare the S–N profiles of the susceptibility, inline, and remanent magnetization models, respectively, inverted using the FTMG data and the TMI data which were computed using the tensor data. One can see that the images produced from FTMG data are more focused and compact and outline the prospective target with better clarity than those generated by inverting the TMI data. We can see the TMI inversions adequately resolve the base of the anomalies; however, the FTMG inversions provide significantly better resolution in the near surface. The inversion of FTMG data also results in stronger magnetic property contrast between the magnetized formations and the host rocks. Thus, the advantages of using the FTMG data are clear. The improvement in resolution is comparable to that of airborne gravity gradiometry (AGG) over the vertical gravity total field [30].

10. Conclusions

We have developed and applied a regularized inversion scheme employing a focusing stabilizer and Gramian constraints to invert FTMG data. We have also introduced a new approach to determining the remanent magnetization from the magnetic field and magnetic gradiometry data. These new methods were applied to the airborne magnetic and magnetic gradiometry survey data acquired in the Thompson Nickel Belt. The significant ultramafic formation in the region is distinctly visible in the data, and the surrounding lithologies are also directly imaged in the remanent magnetization transformation of the full magnetization vector. This technique provides a more comprehensive view of the geology and mineralization of the surveyed region. The study has shown the superiority of the inversion of FTMG data over traditional TMI data, as well as the benefits of imaging various components of the magnetization vector.