Magnetics¶
The geomagnetic field can be ranked as the longest studied of all the geophysical properties of the earth. In addition, magnetic survey, has been used broadly in diverse realm e.g., mining, oil and gas industry and environmental engineering. Although, this geophysical application is quite common in geoscience; however, we do not have modular, welldocumented and welltested opensource codes, which perform forward and inverse problems of magnetic survey. Therefore, here we are going to build up magnetic forward and inverse modeling code based on two common methodologies for forward problem  differential equation and integral equation approaches.
First, we start with some backgrounds of magnetics, e.g., Maxwell’s equations. Based on that secondly, we use differential equation approach to solve forward problem with secondary field formulation. In order to discretzie our system here, we use finite volume approach with weak formulation. Third, we solve inverse problem through GaussNewton method.
Backgrounds¶
Maxwell’s equations for static case with out current source can be written as
where \(\vec{B}\) is magnetic flux (\(T\)) and \(U\) is magnetic potential and \(\mu\) is permeability. Since we do not have any source term in above equations, boundary condition is going to be the driving force of our system as given below
where \(\vec{n}\) means the unit normal vector on the boundary surface (\(\partial \Omega\)). By using seocondary field formulation we can rewrite above equations as
where \(\vec{B}_s\) is the secondary magnetic flux and \(\vec{B}_0\) is the background or primary magnetic flux. In practice, we consider our earth field, which we can get from International Geomagnetic Reference Field (IGRF) by specifying the time and location, as \(\vec{B}_0\). And based on this background fields, we compute secondary fields (\(\vec{B}_s\)). Now we introduce the susceptibility as
Since most materials in the earth have lower permeability than \(\mu_0\), usually \(\chi\) is greater than 0.
Note
Actually, this is an assumption, which means we are not sure exactly this is true, although we are sure, it is very rare that we can encounter those materials. Anyway, practical range of the susceptibility is \(0 < \chi < 1 \).
Since we compute secondary field based on the earth field, which can be different from different locations in the world, we can expect different anomalous responses in different locations in the earth. For instance, assume we have two susceptible spheres, which are exactly same. However, anomalous responses in Seoul and Vancouver are going to be different.
Since we can measure total fields ( \(\vec{B}\)), and usually have reasonably accurate earth field (\(\vec{B}_0\)), we can compute anomalous fields, \(\vec{B}_s\) from our observed data. If you want to download earth magnetic fields at specific location see this website (noaa).
What is our data?¶
In applied geophysics, which means in practice, it is common to refer to measurements as “the magnetic anomaly” and we can consider this as our observed data. For further descriptions in GPG materials for magnetic survey. Now we have the simple relation ship between “the magnetic anomaly” and the total field as
where \(\theta\) is the angle between total and anomalous fields, \(\hat{B}_o\) is the unit vector for \(\vec{B}_o\). Equivalently, we can use the vector dot product to show that the anomalous field is approximately equal to the projection of that field onto the direction of the inducing field. Using this approach we would write
This is important because, in practice we usually use a total field magnetometer (like a proton precession or optically pumped sensor), which can measure only that part of the anomalous field which is in the direction of the earth’s main field.
Sphere in a whole space¶
Forward problem¶
Differential equation approach¶
\[ \begin{align}\begin{aligned}\mathbf{A}\mathbf{u} = \mathbf{rhs}\\\mathbf{A} = \Div(\MfMui)^{1}\Div^{T}\\\mathbf{rhs} = \Div(\MfMui)^{1}\mathbf{M}^f_{\mu_0^{1}}\mathbf{B}_0  \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC}\\\mathbf{B}_s = (\MfMui)^{1}\mathbf{M}^f_{\mu_0^{1}}\mathbf{B}_0\mathbf{B}_0 (\MfMui)^{1}\Div^T \mathbf{u}\end{aligned}\end{align} \]
Mag Differential eq. approach¶

class
SimPEG.PF.Magnetics.
Problem3D_DiffSecondary
(mesh, **kwargs)[source]¶ Bases:
SimPEG.Problem.BaseProblem
Secondary field approach using differential equations!Optional Properties:
 model (
Model
): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)  mu (
PhysicalProperty
): Magnetic Permeability (H/m), a physical property, Default: 1.25663706144e06  muMap (
Mapping
): Mapping of Magnetic Permeability (H/m) to the inversion model., a SimPEG Map  mui (
PhysicalProperty
): Inverse Magnetic Permeability (m/H), a physical property  muiMap (
Mapping
): Mapping of Inverse Magnetic Permeability (m/H) to the inversion model., a SimPEG Map
Other Properties:
 muDeriv (
Derivative
): Derivative of Magnetic Permeability (H/m) wrt the model.  muiDeriv (
Derivative
): Derivative of Inverse Magnetic Permeability (m/H) wrt the model.

surveyPair
¶ alias of
BaseMagSurvey

modelPair
¶ alias of
BaseMagMap

MfMuI
¶

MfMui
¶

MfMu0
¶

getRHS
(m)[source]¶  \[\mathbf{rhs} = \Div(\MfMui)^{1}\mathbf{M}^f_{\mu_0^{1}}\mathbf{B}_0  \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC}\]

getA
(m)[source]¶ GetA creates and returns the A matrix for the Magnetics problem
The A matrix has the form:
\[\mathbf{A} = \Div(\MfMui)^{1}\Div^{T}\]

fields
(m)[source]¶ Return magnetic potential (u) and flux (B) u: defined on the cell center [nC x 1] B: defined on the cell center [nF x 1]
After we compute u, then we update B.
\[\mathbf{B}_s = (\MfMui)^{1}\mathbf{M}^f_{\mu_0^{1}}\mathbf{B}_0\mathbf{B}_0 (\MfMui)^{1}\Div^T \mathbf{u}\]

Jvec
(*args, **kwargs)[source]¶ Computing Jacobian multiplied by vector
By setting our problem as
\[\mathbf{C}(\mathbf{m}, \mathbf{u}) = \mathbf{A}\mathbf{u}  \mathbf{rhs} = 0\]And taking derivative w.r.t m
\[ \begin{align}\begin{aligned}\nabla \mathbf{C}(\mathbf{m}, \mathbf{u}) = \nabla_m \mathbf{C}(\mathbf{m}) \delta \mathbf{m} + \nabla_u \mathbf{C}(\mathbf{u}) \delta \mathbf{u} = 0\\\frac{\delta \mathbf{u}}{\delta \mathbf{m}} =  [\nabla_u \mathbf{C}(\mathbf{u})]^{1}\nabla_m \mathbf{C}(\mathbf{m})\end{aligned}\end{align} \]With some linear algebra we can have
\[ \begin{align}\begin{aligned}\nabla_u \mathbf{C}(\mathbf{u}) = \mathbf{A}\\\nabla_m \mathbf{C}(\mathbf{m}) = \frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u}  \frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}}\end{aligned}\end{align} \]\[ \begin{align}\begin{aligned}\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} = \frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[\Div \diag (\Div^T \mathbf{u}) \dMfMuI \right]\\\dMfMuI = \diag(\MfMui)^{1}_{vec} \mathbf{Av}_{F2CC}^T\diag(\mathbf{v})\diag(\frac{1}{\mu^2})\\\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} = \frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[ \Div \diag(\M^f_{\mu_{0}^{1}}\mathbf{B}_0) \dMfMuI \right]  \diag(\mathbf{v})\mathbf{D} \mathbf{P}_{out}^T\frac{\partial B_{sBC}}{\partial \mathbf{m}}\end{aligned}\end{align} \]In the end,
\[\frac{\delta \mathbf{u}}{\delta \mathbf{m}} =  [ \mathbf{A} ]^{1}\left[ \frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u}  \frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} \right]\]A little tricky point here is we are not interested in potential (u), but interested in magnetic flux (B). Thus, we need sensitivity for B. Now we take derivative of B w.r.t m and have
\[ \begin{align}\begin{aligned}\frac{\delta \mathbf{B}} {\delta \mathbf{m}} = \frac{\partial \mathbf{\mu} } {\partial \mathbf{m} } \left[ \diag(\M^f_{\mu_{0}^{1} } \mathbf{B}_0) \dMfMuI \  \diag (\Div^T\mathbf{u})\dMfMuI \right ]\\  (\MfMui)^{1}\Div^T\frac{\delta\mathbf{u}}{\delta \mathbf{m}}\end{aligned}\end{align} \]Finally we evaluate the above, but we should remember that
Note
We only want to evalute
\[\mathbf{J}\mathbf{v} = \frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}}\mathbf{v}\]Since forming sensitivity matrix is very expensive in that this monster is “big” and “dense” matrix!!

Jtvec
(*args, **kwargs)[source]¶  Computing Jacobian^T multiplied by vector.\[ \begin{align}\begin{aligned}(\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} = \left[ \mathbf{P}_{deriv}\frac{\partial \mathbf{\mu} } {\partial \mathbf{m} } \left[ \diag(\M^f_{\mu_{0}^{1} } \mathbf{B}_0) \dMfMuI \  \diag (\Div^T\mathbf{u})\dMfMuI \right ]\right]^{T}\\  \left[\mathbf{P}_{deriv}(\MfMui)^{1}\Div^T\frac{\delta\mathbf{u}}{\delta \mathbf{m}} \right]^{T}\end{aligned}\end{align} \]
where
\[\mathbf{P}_{derv} = \frac{\partial \mathbf{P}}{\partial\mathbf{B}}\]Note
Here we only want to compute
\[\mathbf{J}^{T}\mathbf{v} = (\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} \mathbf{v}\]

Jtvec_approx
(m, v, f=None)¶ Approximate effect of transpose of J(m) on a vector v.
Parameters:  m (numpy.array) – model
 v (numpy.array) – vector to multiply
 f (Fields) – fields
Return type: Returns: JTv

Jvec_approx
(m, v, f=None)¶ Approximate effect of J(m) on a vector v
Parameters:  m (numpy.array) – model
 v (numpy.array) – vector to multiply
 f (Fields) – fields
Return type: Returns: approxJv

counter
= None¶

curModel
¶ Setting the curModel is depreciated.
Use SimPEG.Problem.model instead.

deleteTheseOnModelUpdate
= []¶

deserialize
(value, trusted=False, verbose=True, **kwargs)¶ Creates HasProperties instance from serialized dictionary
This uses the Property deserializers to deserialize all JSONcompatible dictionary values into their corresponding Property values on a new instance of a HasProperties class. Extra keys in the dictionary that do not correspond to Properties will be ignored.
Parameters:
 value  Dictionary to deserialize new instance from.
 trusted  If True (and if the input dictionary has
'__class__'
keyword and this class is in the registry), the new HasProperties class will come from the dictionary. If False (the default), only the HasProperties class this method is called on will be constructed.  verbose  Raise warnings if
'__class__'
is not found in the registry or of there are unused Property values in the input dictionary. Default is True.  Any other keyword arguments will be passed through to the Property deserializers.

equal
(*args, **kwargs)¶ Determine if two HasProperties instances are equivalent
Equivalence is determined by checking if all Property values on two instances are equal, using
Property.equal
.

ispaired
¶ True if the problem is paired to a survey.

mapPair
¶ alias of
IdentityMap

mapping
¶ Setting an unnamed mapping has been depreciated in v0.4.0. Please see the release notes for more details.

mesh
= None¶

model
¶ Inversion model.

mu
¶ Magnetic Permeability (H/m)

muDeriv
¶ Derivative of Magnetic Permeability (H/m) wrt the model.

muMap
¶ Mapping of Magnetic Permeability (H/m) to the inversion model.

mui
¶ Inverse Magnetic Permeability (m/H)

muiDeriv
¶ Derivative of Inverse Magnetic Permeability (m/H) wrt the model.

muiMap
¶ Mapping of Inverse Magnetic Permeability (m/H) to the inversion model.

needs_model
¶ True if a model is necessary

pair
(d)¶ Bind a survey to this problem instance using pointers.

serialize
(*args, **kwargs)¶ Serializes a HasProperties instance to dictionary
This uses the Property serializers to serialize all Property values to a JSONcompatible dictionary. Properties that are undefined are not included. If the HasProperties instance contains a reference to itself, a
properties.SelfReferenceError
will be raised.Parameters:
 include_class  If True (the default), the name of the class
will also be saved to the serialized dictionary under key
'__class__'
 Any other keyword arguments will be passed through to the Property serializers.
 include_class  If True (the default), the name of the class
will also be saved to the serialized dictionary under key

solverOpts
= {}¶

summary
()¶

survey
¶ The survey object for this problem.

unpair
()¶ Unbind a survey from this problem instance.

validate
(*args, **kwargs)¶ Call all registered class validator methods
These are all methods decorated with
@properties.validator
. Validator methods are expected to raise an error if they fail.
 model (

class
SimPEG.PF.BaseMag.
BaseMagSurvey
(**kwargs)[source]¶ Bases:
SimPEG.Survey.BaseSurvey
Base Magnetics Survey

rxLoc
= None¶ receiver locations

rxType
= None¶ receiver type

Qfx
¶

Qfy
¶

Qfz
¶

projectFields
(u)[source]¶ This function projects the fields onto the data space.
Especially, here for we use total magnetic intensity (TMI) data, which is common in practice.
First we project our B on to data location
\[\mathbf{B}_{rec} = \mathbf{P} \mathbf{B}\]then we take the dot product between B and b_0
\[\text{TMI} = \vec{B}_s \cdot \hat{B}_0\]

projectFieldsDeriv
(*args, **kwargs)[source]¶ This function projects the fields onto the data space.
\[\frac{\partial d_\text{pred}}{\partial \mathbf{B}} = \mathbf{P}\]Especially, this function is for TMI data type

counter
= None¶

dobs
= None¶

dpred
(m, f=None)¶ Create the projected data from a model. The fields, f, (if provided) will be used for the predicted data instead of recalculating the fields (which may be expensive!).
\[d_\text{pred} = P(f(m))\]Where P is a projection of the fields onto the data space.
Note
To use survey.dpred(), SimPEG requires that a problem be bound to the survey. If a problem has not been bound, an Exception will be raised. To bind a problem to the Data object:
survey.pair(myProblem)

dtrue
= None¶

eps
= None¶

eval
(f)¶ This function projects the fields onto the data space.
\[d_\text{pred} = \mathbf{P} f(m)\]

evalDeriv
(f)¶ This function s the derivative of projects the fields onto the data space.
\[\frac{\partial d_\text{pred}}{\partial u} = \mathbf{P}\]

getSourceIndex
(sources)¶

isSynthetic
¶ Check if the data is synthetic.

ispaired
¶

makeSyntheticData
(m, std=0.05, f=None, force=False)¶ Make synthetic data given a model, and a standard deviation.
Parameters:  m (numpy.array) – geophysical model
 std (numpy.array) – standard deviation
 u (numpy.array) – fields for the given model (if precalculated)
 force (bool) – force overwriting of dobs

mesh
¶ Mesh of the paired problem.

mtrue
= None¶

nD
¶ Number of data

nSrc
¶ Number of Sources

pair
(p)¶ Bind a problem to this survey instance using pointers

prob
¶ The geophysical problem that explains this survey, use:
survey.pair(prob)

residual
(m, f=None)¶ Parameters:  m (numpy.array) – geophysical model
 f (numpy.array) – fields
Return type: Returns: data residual
The data residual:
\[\mu_\text{data} = \mathbf{d}_\text{pred}  \mathbf{d}_\text{obs}\]

srcList
¶ Source List

srcPair
¶ alias of
BaseSrc

std
= None¶

unpair
()¶ Unbind a problem from this survey instance

vnD
¶ Vector number of data

Mag Integral eq. approach¶
Mag analytic solutions¶

SimPEG.PF.MagAnalytics.
spheremodel
(mesh, x0, y0, z0, r)[source]¶ Generate model indicies for sphere  (x0, y0, z0 ): is the center location of sphere  r: is the radius of the sphere  it returns logical indicies of cellcenter model

SimPEG.PF.MagAnalytics.
MagSphereAnaFun
(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag='total')[source]¶ test Analytic function for Magnetics problem. The set up here is magnetic sphere in wholespace assuming that the inducing field is oriented in the xdirection.
 (x0, y0, z0)
 (x0, y0, z0 ): is the center location of sphere
 r: is the radius of the sphere
\[\mathbf{H}_0 = H_0\hat{x}\]

SimPEG.PF.MagAnalytics.
CongruousMagBC
(mesh, Bo, chi)[source]¶ Computing boundary condition using Congrous sphere method. This is designed for secondary field formulation.
>> Input
 mesh: Mesh class
 Bo: np.array([Box, Boy, Boz]): Primary magnetic flux
 chi: susceptibility at cell volume
\[\vec{B}(r) = \frac{\mu_0}{4\pi} \frac{m}{ \ \vec{r}  \vec{r}_0\^3}[3\hat{m}\cdot\hat{r}\hat{m}]\]

SimPEG.PF.MagAnalytics.
MagSphereAnaFunA
(x, y, z, R, xc, yc, zc, chi, Bo, flag)[source]¶ Computing boundary condition using Congrous sphere method. This is designed for secondary field formulation. >> Input mesh: Mesh class Bo: np.array([Box, Boy, Boz]): Primary magnetic flux Chi: susceptibility at cell volume
\[\vec{B}(r) = \frac{\mu_0}{4\pi}\frac{m}{\ \vec{r}\vec{r}_0\^3}[3\hat{m}\cdot\hat{r}\hat{m}]\]

SimPEG.PF.MagAnalytics.
IDTtoxyz
(Inc, Dec, Btot)[source]¶ Convert from Inclination, Declination, Total intensity of earth field to x, y, z

SimPEG.PF.MagAnalytics.
MagSphereFreeSpace
(x, y, z, R, xc, yc, zc, chi, Bo)[source]¶ Computing the induced response of magnetic sphere in freespace.
>> Input x, y, z: Observation locations R: radius of the sphere xc, yc, zc: Location of the sphere chi: Susceptibility of sphere Bo: Inducing field components [bx, by, bz]*H0