# Dipole recovery from sparse measurements of its field on a cylindrical geometry

Sylvain Chevillard, Juliette Leblond, Konstantinos Mavreas

► To cite this version:

Sylvain Chevillard, Juliette Leblond, Konstantinos Mavreas. Dipole recovery from sparse measurements of its field on a cylindrical geometry. International Journal of Applied Electromagnetics and Mechanics, In press, Special issue, proceedings of the conference ISEM 2017, 59 (1), pp.7. 10.3233/JAE-171033. hal-01618885v2

HAL Id: hal-01618885

https://inria.hal.science/hal-01618885v2

Submitted on 10 Dec 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Dipole recovery from sparse measurements of its field on a cylindrical geometry

Sylvain Chevillard Juliette Leblond

Konstantinos Mavreas

Inria, 2004 route des Lucioles, BP 93, 06902, Sophia Antipolis Cedex, France

## Abstract

We consider the inverse problem of recovering the position and moment of a magnetic dipole from sparse measurements of the field it generates, known on sections of three orthogonal cylinders enclosing it. This problem is motivated by recent measurements performed on Moon rocks, in view of determining their magnetic properties. The key ingredient of the presented method is the use of rational approximation techniques, together with properties of the poles of the approximants, in order to estimate the position of the dipole.

**Keywords:** Inverse problems, magnetic dipole, source estimation, rational approximation, Moon rocks.

## 1 Introduction

Although today the Moon has no global magnetic field, paleomagnetic studies of Moon rocks (available at NASA, from APOLLO missions) [6] reveal magnetic anomalies on its surface, which indicate that it used to have one for billions of years. In order to clarify the evolution of the Lunar magnetic field, geologists at CEREGE† want to select a number of rock samples that possess paleomagnetic recording potentials, namely those that admit a magnetic moment of sufficiently high amplitude, in order to study them more carefully. However, Moon rocks must be protected from any kind of damage whence the inspection protocol must be nondestructive. Moreover, the preliminary selection process is also constrained to be fast, since there are numerous samples to examine and in order to reduce the cost of measurement time. For this purpose, a portable magnetometer (“lunometer”) has been constructed [5]. It furnishes sparse measurements of the components of the magnetic field generated by the sample in a specific cylindrical geometry, from which the order of magnitude of the remanent magnetization is to be estimated. This is the inverse problem that we address.

We first describe in Section 2 the data acquisition process and the physical model. We then decompose the initial problem into two sub-problems: the estimation of the source location and the recovery of the magnetic moment. The moment estimation is a rather easy task when the location is well estimated, as explained in Section 3. In order to recover the location, we propose an approach based on best quadratic rational approximation techniques in Section 4. Finally, some simulations with synthetic data are provided in Section 5, including comparisons with a more naive approach.

<footnote>†Centre Européen de Recherche et d'Enseignement des Géosciences de l'Environnement, CNRS. Aix en Provence, France.</footnote>
## 2 Data geometry and physical model

The lunometer is made of a horizontal circular deposition plate that can be rotated along the vertical axis corresponding to its center. The rock sample is enclosed in a nonmagnetic cubic box, sat on the deposition plate. At some distance $R$ from it, three flux-gate probes are vertically aligned. Each probe measures a given component of the magnetic field where it lies: either the vertical component, or the horizontal component along the direction from the center of the deposition plate to the probe, or the horizontal component orthogonal to the latter. The experimental protocol is explained in details in [5]. To sum up, it is as follows. We represent points of the space in the reference frame of the cube, with the convention that (0, 0, 0) corresponds to the center of the cube. In the beginning of the experiment, the cube lies in the center of the deposition plate, so that the revolution axis of the plate corresponds to the line (0, 0, $\mathbb{R}$) and the three probes lie at coordinates $(R, 0, h_k)$, $k = 1, 2, 3$. The cube is isolated from external electromagnetic fields with a mu-metal shield, so that one can assume that the flux-gate probes actually measure only the field produced by the sample. The deposition plate is rotated around its axis by angle steps of 1° until a full turn has been made. In the reference frame of the cube, this corresponds therefore to measurements of the field performed on three sections of the cylinder centered in 0 with radius $R$ and revolution axis (0, 0, $\mathbb{R}$). On each of these circles, a component of the magnetic field is measured at $N = 360$ points. The measured component is either in the direction of the axis of the cylinder, or tangent or radial to the circle, depending on its height $h_k$. Then the cube is rotated, so that another of its face lies on the deposition plate and the same operation is performed. In the reference frame of the cube, this corresponds to measurements performed on sections of a cylinder with revolution axis, e.g., ($\mathbb{R}, 0, 0$). Finally, the cube is rotated again in order to get measurements corresponding to the cylinder with axis (0, $\mathbb{R}, 0$). Overall, $9N$ measurements are available ($N$ measurements on 3 × 3 circles). This is summed up in Figure 1.

<figure><img src="image_1.png" /><figcaption>Figure 1: Data geometry: measurements are performed on three orthogonal cylinders; for each cylinder, data are indeed available at discrete points of three circles. The places where the field is measured are marked with small dots. Each color indicates which component of the field is indeed measured (black: radial, red: tangential, blue: along the axis).</figcaption></figure>
We make the assumption that the magnetic field produced by the sample can be assimilated to the field produced by a unique pointwise dipolar magnetic source located at $X_d \in \mathbb{R}^3$, with moment $M_d \in \mathbb{R}^3$. Our goal is to recover $X_d$ and $M_d$ whenever the dipolar assumption is verified, and to test whether or not the assumption is too much a simplification to explain the true field.

The underlying magnetic phenomenon is modeled by Maxwell equations in the magneto-static and macroscopic framework and in particular by Biot-Savart law (see, e.g., [3]). Denoting the permeability of the free space with $\mu_0 = 4\pi \cdot 10^{-7} \text{Tm/A}$, the magnetic field at $X \neq X_d$ then has the form

$$ \mathbf{B}_{[X_d, M_d]}(\mathbf{X}) = -\frac{\mu_0 |\mathbf{X} - \mathbf{X}_d|^2 M_d - 3[\mathbf{M}_d \cdot (\mathbf{X} - \mathbf{X}_d)] (\mathbf{X} - \mathbf{X}_d)}{4\pi |\mathbf{X} - \mathbf{X}_d|^5}, \quad (1) $$

where, if $U = (u_1, u_2, u_3)$ and $V = (v_1, v_2, v_3)$, the quantity $U \cdot V$ denotes the scalar product $u_1v_1 + u_2v_2 + u_3v_3$ and $|U|$ denotes the Euclidean norm $(u_1^2 + u_2^2 + u_3^2)^{1/2}$.

### 3 Magnetic moment recovery from estimated dipole position

A very rough estimate of the dipole position consists in considering that $X_d$ is approximately $(0, 0, 0)$. Indeed, if the rock sample is small enough compared to the radius $R$ and uniformly magnetized, it appears in first approximation as a dipole located at its center of mass, which roughly corresponds to the center of the cube $(0, 0, 0)$. This is the method used by the scientists from CEREGE so far [5]. Our contribution, in Section 4, will be to propose a method to better estimate $X_d$ but, for the moment being, let us simply assume that $X_d$ is specified (as the center of the sample or as an estimated point). We see from Equation (1) that the recovery of the magnetic moment $M_d$ from measurements of $B_{[X_d, M_d]}$ at sensors locations is a linear problem. Indeed, we dispose of $9N$ measurements performed at $9N$ different points $X$; each measurement is a component of $B_{[X_d, M_d]}$ along some known direction $v$, meaning that, on the one hand, we know the value $b = B_{[X_d, M_d]} \cdot v$. On the other hand, $X$ and $X_d$ being known, Equation (1) leads to $b = a_1m_1 + a_2m_2 + a_3m_3$ where $(m_1, m_2, m_3)$ is the unknown moment $M_d$ and $a_1, a_2$ and $a_3$ are constants explicitly computed from $v$, $X$ and $X_d$. These relations, for all available measurements, can be arranged into an over-determined system $A M_d = \mathfrak{B}$, where $\mathfrak{B}$ denotes the vector with $9N$ elements formed with the different values of $b$ and $A$ is the $9N \times 3$ matrix formed with the different values of $a_1, a_2$ and $a_3$. We numerically solve it, as is classical, by finding the solution $M_d$ of the least squares problem $M_d = \text{Argmin}_{M \in \mathbb{R}^3} |\mathfrak{B} - AM|^2$, where $| \cdot |$ is the Euclidean norm. This can be done, e.g., using the Moore-Penrose pseudo-inverse [2, Sec. 5.5.4]. In this situation, the Moore-Penrose pseudo-inverse matrix $A^+$ of $A$ is given by the expression $A^+ = (A^t A)^{-1} A^t$ where $A^t$ is the transposed matrix of $A$. The solution to the least squares problem is then given by $M_d = A^+ \mathfrak{B}$.

### 4 Magnetic dipole localization

In view of recovering the location $X_d$ of the dipole, and since the field is known on circles, we adapt and apply the ideas of [1] to our context (see also the references there and [4] for other localization procedures). The key observation is that the denominator of Equation (1) is not polluted by the moment of the dipole, together with the fact that it can be linked to the pole of some rational function whose definition only relies on the data available on each of the 9 circles.

Let us consider the measurements provided on one circle, e.g., the circle where the sensor positions are of the form $X = (x_1, x_2, h)$ with the couple $(x_1, x_2)$ actually varying on the circle
centered at (0,0) of radius R. We introduce the complex number $\xi = x_1 + jx_2$, where j ∈ C is the imaginary number $j^2 = -1$, so that the measurements are indeed known for values of $\xi$ that all satisfy $|\xi| = R$. Accordingly, we note $\mathbf{X}_d = (x_{d,1}, x_{d,2}, x_{d,3})$ and we introduce the complex number $\xi_d = x_{d,1} + jx_{d,2}$. We assume that we are in the generic case of $\xi_d \neq 0$. Moreover the magnetic source that we are looking for is necessary inside the sample, which is itself inside the cubic box surrounded by the circle of radius R where the sensors lie. This implies that $|x_{d,3}| < R$ and $|\xi_d| < R$.

Now, for any $\xi$ satisfying $|\xi| = R$ we have $\bar{\xi} = R^2/\xi$, whence

$$|\mathbf{X} - \mathbf{X}_d|^2 = |\xi - \xi_d|^2 + (h - x_{d,3})^2 = (\xi - \xi_d) \left( \frac{R^2}{\xi} - \bar{\xi}_d \right) + (h - x_{d,3})^2 = -\frac{\bar{\xi}_d}{\xi} p_h(\xi), \quad (2)$$

where $p_h$ is the polynomial of degree 2 in variable $\xi$ given by

$$p_h(\xi) = \xi^2 - \frac{|\xi_d|^2 + (h - x_{d,3})^2 + R^2}{\bar{\xi}_d} \xi + \frac{\xi_d}{\bar{\xi}_d} R^2. \quad (3)$$

Accordingly, if $\mathbf{c} = (c_1, c_2, c_3)$ is an arbitrary vector of $\mathbb{R}^3$ with $\gamma = c_1 + j c_2$, and observing that $x_1 = (\xi + \bar{\xi})/2$ and $x_2 = (\xi - \bar{\xi})/(2j)$, we have, for any $\xi$ satisfying $|\xi| = R$,

$$\mathbf{X} \cdot \mathbf{c} = \frac{1}{2\xi} (\bar{\gamma} \xi^2 + 2 c_3 h \xi + \gamma R^2). \quad (4)$$

Now, the measurement available at point $\mathbf{X}$ is $b = \mathbf{B}[\mathbf{X}_d, \mathbf{M}_d] \cdot \mathbf{v}$ with $\mathbf{v}$ being either $\mathbf{v}_\perp = (0,0,1)$ or $\mathbf{v}_r = (x_1/R, x_2/R, 0)$ or $\mathbf{v}_\tau = (-x_2/R, x_1/R, 0)$. Observing that $\mathbf{X} \cdot \mathbf{v}_\perp = h$, $\mathbf{X} \cdot \mathbf{v}_r = R$, $\mathbf{X} \cdot \mathbf{v}_\tau = 0$, $\mathbf{v}_\perp \cdot \mathbf{c} = c_3$, $\mathbf{v}_r \cdot \mathbf{c} = \mathbf{X} \cdot (c_1/R, c_2/R, 0)$ and $\mathbf{v}_\tau \cdot \mathbf{c} = \mathbf{X} \cdot (-c_2/R, c_1/R, 0)$, we see from Equation (1) together with (2) and (4) that $b^2$ is of the form $\xi q(\xi)/p_h(\xi)^5$ where $q$ is a polynomial of degree 8 that depends on the measured component and the height $h$.

We claim that the two complex roots $\xi_-$ and $\xi_+$ of $p_h$ satisfy $|\xi_-| \le |\xi_d| \le R \le |\xi_+$|, that $\xi_\pm$ and $\xi_d$ share the same complex argument and that, when $h$ varies, $|\xi_-|$ reaches its maximal value $|\xi_d|$ once, at $h = x_{d,3}$. Indeed, doing the change of variable $\xi = t\xi_d/|\xi_d|$, the polynomial $p_h$ rewrites as

$$p_h(\xi) = \frac{\xi_d}{\bar{\xi}_d} (t^2 - \lambda(h)t + R^2) \quad \text{where} \quad \lambda(h) = \frac{|\xi_d|^2 + (h - x_{d,3})^2 + R^2}{|\xi_d|} > 2R.$$

The fact that $\lambda(h) > 2R$ follows from $|\xi_d|^2 + R^2 = (|\xi_d| - R)^2 + 2R|\xi_d|$ and implies that the polynomial $t^2 - \lambda(h)t + R^2$ with real-valued coefficients has positive discriminant. Hence, it admits two real roots $t_-$ and $t_+$ that satisfy $t_+ t_- = R^2$ and $t_+ + t_- = \lambda(h)$. Choosing by convention that $t_- \le t_+$, we hence have $0 < t_- \le R \le t_+$. Since $|\xi| = |t|$, this establishes half of our claim. Now, the explicit formula for a second order equation gives $2t_- = \lambda - \sqrt{\lambda^2 - 4R^2}$, which is a decreasing function of $\lambda$, whence $t_-$ is maximal when $\lambda$ is minimal, which obviously happens only once, for $h = x_{d,3}$. Finally, when $h = x_{d,3}$ we have $\lambda(h)|\xi_d| = |\xi_d|^2 + R^2$ whence $|\xi_d|$ is a root of $t^2 - \lambda(h)t + R^2$, and since $|\xi_d| < R$, it must hold that $|\xi_d| = t_-$, thereby establishing the overall claim.

We can now describe our recovery algorithm for the location $\mathbf{X}_d$ of the dipole: we have seen that the square $b^2$ of the measurements on a given circle at height $h$ is the trace on the circle centered at the origin and of radius R of a rational function with one pole $\xi_-$ of order 5 inside the disc and one pole $\xi_+$ of order 5 outside the disc. Following the same strategy as in [1, Sec. C.1], we use an algorithm of best quadratic rational approximation with constrained poles, in order to recover a numerical estimate of $\xi_-(h)$. For this purpose, we use the RARL2‡ software.

<footnote>‡https://project.inria.fr/rarl2/</footnote>
Using the properties of the root $\xi_-$ with respect to $h$ as discussed above, we select among the estimated $\xi_-(h_k)$ for $k = 1, 2, 3$ the one with maximum modulus, which is expected to be the closest to $\xi_d$, thereby obtaining a rough approximation of $x_{d,1}$ and $x_{d,2}$.

Repeating the same procedure with the circles corresponding to the second position of the cubic box leads a rough approximation of, e.g., $x_{d,2}$ and $x_{d,3}$. Finally, the circles corresponding to the third position provide us with approximations of $x_{d,1}$ and $x_{d,3}$. Overall, we have two estimates for each of the three coordinates of $X_d$. For each coordinate, we retain the estimate of maximal absolute value.

# 5 Numerical simulations with synthetic data

We conducted experiments using Matlab, with several locations $X_d$ and moments $M_d$. For each configuration, we built our data using Equation (1). The parameters values are realistic: $R = 111$ mm and the sections heights are $h_1 = 0$ mm, $h_2 = 15$ mm, $h_3 = 30$ mm. For each example, we estimated the position $X_d$ using the algorithm described in Section 4, leading to an estimate $X_{est}$. Then we estimated the moment using the approach of Section 3: first using the rough estimate $X_d \simeq (0,0,0)$ as is done in [5], and second using the estimate $X_d \simeq X_{est}$ (this gives the estimates $M_{rough}$ and $M_{est}$ respectively).

<table><caption>Table 1: Four examples: the field of a perfect dipole at position $X_d$ (expressed in mm) and moment $M_d$ (expressed in A·m²) is generated on the geometry described in Section 2, so as to simulate real measurements. The quantities $X_{est}$ and $M_{est}$ are our new estimates, whereas $M_{rough}$ is the estimate of the moment obtained with the method currently in use at CEREGE.</caption><thead><tr><th>k</th><th>X<sub>d</sub><sup>(k)</sup></th><th>X<sub>est</sub><sup>(k)</sup></th><th>|X<sub>est</sub><sup>(k)</sup> - X<sub>d</sub><sup>(k)</sup>| / |X<sub>d</sub><sup>(k)</sup>|</th><th>M<sub>d</sub><sup>(k)</sup></th><th>M<sub>rough</sub><sup>(k)</sup></th><th>M<sub>est</sub><sup>(k)</sup></th><th>|M<sub>rough</sub><sup>(k)</sup> - M<sub>d</sub><sup>(k)</sup>| / |M<sub>d</sub><sup>(k)</sup>|</th><th>|M<sub>est</sub><sup>(k)</sup> - M<sub>d</sub><sup>(k)</sup>| / |M<sub>d</sub><sup>(k)</sup>|</th></tr></thead><tbody><tr><td rowspan="3">1</td><td>(28)</td><td>(28)</td><td rowspan="3">3.8%</td><td>(0.05)</td><td>(0.034)</td><td>(0.05)</td><td rowspan="3">29%</td><td rowspan="3">2.4%</td></tr><tr><td>(28)</td><td>(28)</td><td>(0.05)</td><td>(0.034)</td><td>(0.05)</td></tr><tr><td>(28)</td><td>(26.1)</td><td>(0.05)</td><td>(0.039)</td><td>(0.048)</td></tr><tr><td rowspan="3">2</td><td>(46)</td><td>(44.9)</td><td rowspan="3">5.1%</td><td>(0.017)</td><td>(-0.023)</td><td>(0.018)</td><td rowspan="3">58%</td><td rowspan="3">5.2%</td></tr><tr><td>(26)</td><td>(25.4)</td><td>(0.047)</td><td>(0.019)</td><td>(0.046)</td></tr><tr><td>(45)</td><td>(41.7)</td><td>(0.086)</td><td>(0.054)</td><td>(0.081)</td></tr><tr><td rowspan="3">3</td><td>(46)</td><td>(44.9)</td><td rowspan="3">5.1%</td><td>(0.067)</td><td>(0.017)</td><td>(0.065)</td><td rowspan="3">73%</td><td rowspan="3">4.1%</td></tr><tr><td>(26)</td><td>(25.4)</td><td>(0.069)</td><td>(0.032)</td><td>(0.067)</td></tr><tr><td>(45)</td><td>(41.7)</td><td>(0.017)</td><td>(-0.018)</td><td>(0.014)</td></tr><tr><td rowspan="3">4</td><td>(26.5)</td><td>(26.1)</td><td rowspan="3">4.8%</td><td>(-0.026)</td><td>(-0.017)</td><td>(-0.024)</td><td rowspan="3">126%</td><td rowspan="3">6.1%</td></tr><tr><td>(-37.1)</td><td>(-32.4)</td><td>(-0.079)</td><td>(-0.031)</td><td>(-0.074)</td></tr><tr><td>(-89)</td><td>(-87.9)</td><td>(-0.034)</td><td>(0.069)</td><td>(-0.035)</td></tr></tbody></table>

To illustrate what happens when the dipolar assumption is obviously not fulfilled, we built a fifth example, corresponding to the field produced by the dipoles of examples 1 and 4 simultaneously. Our method recovers a position $X_{est}^{(5)} = (34, -32, -88)$ and a moment $M_{est}^{(5)} = (0.259, 0.06, 0.046)$. Figure 2(a) shows that the field produced by the recovered dipole fits usually well the data when the dipolar assumption is satisfied; in contrast, Figure 2(b) shows that it is not the case anymore for the example 5. This gives a practical way of testing afterward whether the dipolar assumption was realistic or not.

One can notice that, while the position $X_d$ is the same in Examples 3 and 4, the obtained
<figure><img src="image_2.png" /><figcaption>Figure 2: For two examples, plot of the data compared to the field that would be produced by the recovered dipoles. Data correspond to the second position of the cubic box and the sensor at height $h_1 = 0$mm (radial component) and is plotted with respect to the angle of the deposition plate (in degrees).</figcaption></figure>

estimate $\mathbf{X}_{\text{est}}$ is not exactly the same: this is expected, since the estimate relies on recovering a rational function that depends not only on $\mathbf{X}_d$ but also on $\mathbf{M}_d$ (even though the denominator remains unchanged). The numerical behavior is hence slightly different and leads to two different results.

The obtained results clearly show that the recovery of the moment highly depends on the quality of the estimate of the position: when $\mathbf{X}_d$ is not close to zero, $\mathbf{M}_{\text{est}}$ is usually much more accurate than $\mathbf{M}_{\text{rough}}$.

# 6 Conclusion

In the situation of a pointwise dipolar magnetization, we saw that the quality of the moment estimation is improved by the preliminary localization of the dipole. This first approximation step is currently being refined, using the relations derived in Section 4, and we study its influence on the quality of the moment estimation. We are also working on the analysis of the sensitivity of our estimation procedure with respect to the characteristics of the dipole and also to the data acquisition protocol. In particular, it could be helpful to have additional available measurements of the magnetic field by adding sensors at some other heights. Last, but not least, we plan to investigate situations with several pointwise dipoles, in order to better reflect the magnetic behavior of composite rocks.

**Acknowledgments:** This research is part of the ANR project MagLune.

# References

[1] M. Clerc, J. Leblond, J.-P. Marmorat, and T. Papadopoulos. Source localization using rational approximation on plane sections. *Inverse Problems*, 28, 2012.

[2] G. H. Golub and Ch. F. van Loan. *Matrix computations*. John Hopkins University Press, third edition, 1996.

[3] J. D. Jackson. *Classical electrodynamics*. John Wiley & Sons, third edition, 2001.
[4] J. C. Mosher and R. M. Leahy. Recursive MUSIC: a framework for EEG and MEG source localization. *IEEE Trans. Biomed. Eng.*, 45(11), 1998.

[5] M. Uehara, J. Gattacceca, Y. Quesnel, C. Lepaulard, E. A. Lima, M. Manfredi, and P. Rochelette. A spinner magnetometer for large Apollo lunar samples. *Review of Scientific Instruments*, *American Institute of Physics*, *88*(10), 2017.

[6] B. P. Weiss and S. M. Tikoo. The lunar dynamo. *Science*, 346(6214), 2014.