Identifying variations of magnetic anomalies using geomagnetic monitoring

TL;DR

This paper rigorously proves that geomagnetic monitoring can uniquely identify the locations, variation rates, and material properties of magnetic anomalies beneath the Earth, extending previous results to include anomalies with changing parameters.

Contribution

It introduces a rigorous mathematical framework for uniquely recovering both static and varying magnetic anomalies using geomagnetic data.

Findings

Unique recovery of anomaly locations and parameters

Extension of previous results to dynamic anomalies

Mathematical proof of identifiability in inverse problems

Abstract

We are concerned with the inverse problem of identifying magnetic anomalies with varing parameters beneath the Earth using geomagnetic monitoring. Observations of the change in Earth's magnetic field--the secular variation--provide information about the anomalies as well as their variations. In this paper, we rigorously establish the unique recovery results for this magnetic anomaly detection problem. We show that one can uniquely recover the locations, the variation parameters including the growth or decaying rates as well as their material parameters of the anomalies. This paper extends the existing results in [8] by two of the authors to varying anomalies.

Equations

$$\begin{array}{rl}\sigma (\mathbf{x})& ={\sigma }_c(\mathbf{x})\chi ({\Sigma }_c),\\ \mu (\mathbf{x})& =({\mu }_c(\mathbf{x})-{\mu }_0)\chi ({\Sigma }_c)+{\mu }_0,\\ ϵ(\mathbf{x})& =({ϵ}_c(\mathbf{x})-{ϵ}_0)+({ϵ}_s(\mathbf{x})-{ϵ}_0)\chi ({\Sigma }_s)+{ϵ}_0.\end{array}$$

$$D_l=s_l\delta \Omega +{\mathbf{z}}_l,l=1,2,\dots ,l_0,$$

$$s_l={\delta }^{{\alpha }_l},l=1,2,\dots ,l_0,$$

$$\begin{array}{rl}\sigma (\mathbf{x})=& {\sigma }_c(\mathbf{x})\chi ({\Sigma }_c)+\sum _{l=1}^{l_0}{\sigma }_l\chi (D_l),\\ \mu (\mathbf{x})=& \left({\mu }_c(\mathbf{x})-{\mu }_0\right)\chi ({\Sigma }_c)+\sum _{l=1}^{l_0}({\mu }_l-{\mu }_0)\chi (D_l)+{\mu }_0,\\ ϵ(\mathbf{x})=& \left({ϵ}_c(\mathbf{x})-{ϵ}_0\right)\chi ({\Sigma }_c)+\left({ϵ}_s(\mathbf{x})-{ϵ}_0\right)\chi \left({\Sigma }_s\setminus \overline{\bigcup _{l=1}^{l_0}{D}_l}\right)\\ & +\sum _{l=1}^{l_0}({ϵ}_l-{ϵ}_0)\chi (D_l)+{ϵ}_0.\end{array}$$

$$\{\begin{array}{ll}\nabla \times \mathbf{H}=-i\omega ϵ\mathbf{E}+\sigma (\mathbf{E}+\mu \mathbf{v}\times \mathbf{H})& \text{in }{\mathbb{R}}^3,\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \nabla \times \mathbf{E}=i\omega \mu \mathbf{H}& \text{in }{\mathbb{R}}^3,\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \nabla \cdot (\mu \mathbf{H})=0,\nabla \cdot (ϵ\mathbf{E})=\hat{\rho }& \text{in }{\mathbb{R}}^3,\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \underset{\Vert x\Vert \to \infty }{\lim }\Vert \mathbf{x}\Vert (\sqrt{{\mu }_0}\mathbf{H}\times \hat{\mathbf{x}}-\sqrt{{ϵ}_0}\mathbf{E})=0.& \end{array}$$

$$H(\mathrm{curl};\mathcal{M}):=\{U\in (L^2(\mathcal{M}){)}^3;\nabla \times U\in (L^2(\mathcal{M}){)}^3\}$$

$$H_{loc}(\mathrm{curl};X):=\{U{|}_{\mathcal{M}}\in H(\mathrm{curl};\mathcal{M});\mathcal{M}\text{ is any bounded subdomain of }X\}.$$

$$\mathbf{E}(\mathbf{x})=\frac{e^{i{k}_0\Vert \mathbf{x}\Vert }}{\Vert \mathbf{x}\Vert }{\mathbf{E}}^{\infty }(\hat{\mathbf{x}})+O\left(\frac{1}{\Vert \mathbf{x}{\Vert }^2}\right),$$

$$\mathbf{H}(\mathbf{x})=\frac{e^{i{k}_0\Vert \mathbf{x}\Vert }}{\Vert \mathbf{x}\Vert }{\mathbf{H}}^{\infty }(\hat{\mathbf{x}})+O(\frac{1}{\Vert \mathbf{x}{\Vert }^2}),$$

$$(\mathbf{H}_{s}^{\infty}(\hat{\mathbf{x}};\omega)-\mathbf{H}^{\infty}(\hat{x};\omega))\Big{|}_{(\hat{\mathbf{x}},\omega)\in\widehat{\Gamma}\times\mathbb{R}_{+}}\longrightarrow\bigcup_{l=1}^{l_{0}}(D_{l};\mu_{l},\sigma_{l},\alpha_{l},\mathbf{z}_{l}),$$

$$(\mathbf{H}_{s}(\mathbf{x},\omega)-\mathbf{H}(\mathbf{x},\omega))\Big{|}_{(\mathbf{x},\omega)\in\partial\Sigma\times\mathbb{R}_{+}}\rightarrow\bigcup_{l=1}^{l_{0}}(D_{l};\mu_{l},\sigma_{l},\alpha_{l},\mathbf{z}_{l}).$$

$$\{\begin{array}{l}\nabla \times \mathcal{H}(\mathbf{x},t)=ϵ(\mathbf{x}){\partial }_t\mathcal{E}(\mathbf{x},t)+\sigma (\mathbf{x})\left(\mathcal{E}(\mathbf{x},t)+\mu (\mathbf{x})\mathbf{v}\times \mathcal{H}(\mathbf{x},t)\right),\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \nabla \times \mathcal{E}(\mathbf{x},t)=-\mu (\mathbf{x}){\partial }_t\mathcal{H}(\mathbf{x},t),\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \nabla \cdot (\mu (\mathbf{x})\mathcal{H}(\mathbf{x},t))=0,\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \nabla \cdot (ϵ(\mathbf{x})\mathcal{E}(\mathbf{x},t))=\rho (\mathbf{x},t),\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \mathcal{E}(\mathbf{x},0)=\mathcal{H}(\mathbf{x},0)=0,\end{array}$$

$$(\mathcal{H}_{s}(\mathbf{x},t)-\mathcal{H}(\mathbf{x},t))\Big{|}_{(\mathbf{x},t)\in\Gamma\times\mathbb{R}_{+}}\rightarrow\bigcup_{l=1}^{l_{0}}(D_{l};\mu_{l},\sigma_{l},\alpha_{l},\mathbf{z}_{l}).$$

$$\mathbf{J}(\mathbf{x},\omega )={\mathcal{F}}_t(\mathcal{J}):=\frac{1}{2\pi }{\int }_0^{\infty }\mathcal{J}(\mathbf{x},t)e^{i\omega t}dt,\mathcal{J}=\mathcal{E}\text{ or }\mathcal{H},$$

$$\begin{array}{rl}\mathbf{E}={\mathcal{F}}_t(\mathcal{E}),\mathbf{H}={\mathcal{F}}_t(\mathcal{H}),& {\partial }_t\mathbf{E}=-i\omega {\mathcal{F}}_t(\mathcal{E}),{\partial }_t\mathbf{H}=-i\omega {\mathcal{F}}_t(\mathcal{H}),\\ {\mathbf{E}}_s={\mathcal{F}}_t({\mathcal{E}}_s),{\mathbf{H}}_s={\mathcal{F}}_t({\mathcal{H}}_s),{\partial }_t{\mathbf{E}}_s& =-i\omega {\mathcal{F}}_t({\mathcal{E}}_s),{\partial }_t{\mathbf{H}}_s=-i\omega {\mathcal{F}}_t({\mathcal{H}}_s),\hat{\rho }={\mathcal{F}}_t(\rho ),\end{array}$$

$$L_T^2(\partial B):=\{\Phi \in L^2(\partial B{)}^3,\nu \cdot \Phi =0\}.$$

$$\mathrm{TH}(\mathrm{div},\partial B)$$

$$\mathrm{TH}(\mathrm{curl},\partial B)$$

$$\Vert \Phi {\Vert }_{\text{TH(div,}\partial B)}$$

$${\Gamma }_k(\mathbf{x})=-\frac{e^{ik\Vert \mathbf{x}\Vert }}{4\pi \Vert \mathbf{x}\Vert },\mathbf{x}\in {\mathbb{R}}^3\text{and}\mathbf{x}\ne 0.$$

$$S_B^k[\varphi ](\mathbf{x}):={\int }_{\partial B}{\Gamma }_k(\mathbf{x}-\mathbf{y})\varphi (\mathbf{y})d{s}_{\mathbf{y}},\mathbf{x}\in {\mathbb{R}}^3\setminus \partial B,$$

$${\mathcal{K}}_B^k[\varphi ](\mathbf{x}):=p.v.{\int }_{\partial B}\frac{\partial {\Gamma }_k(\mathbf{x}-\mathbf{y})}{\partial {\nu }_{\mathbf{y}}}\varphi (\mathbf{y})d{s}_{\mathbf{y}},\mathbf{x}\in \partial B,$$

$$\frac{\partial}{\partial\nu}S^{k}_{B}[\phi]\Big{|}_{\pm}=\left(\pm\frac{1}{2}I+(\mathcal{K}^{k}_{B})^{*}\right)[\phi]\quad\text{on}\ \partial B,$$

$${\mathcal{A}}_B^k[\Phi ](\mathbf{x}):={\int }_{\partial B}{\Gamma }_k(\mathbf{x}-\mathbf{y})\Phi (\mathbf{y})d{s}_{\mathbf{y}},\mathbf{x}\in {\mathbb{R}}^3\setminus \partial B.$$

$$\nu\times\nabla\times\mathcal{A}_{B}^{k}[\Phi]\Big{|}_{\pm}=\mp\frac{\Phi}{2}+\mathcal{M}^{k}_{B}[\Phi]\quad\text{on }\partial B,$$

$$\forall\mathbf{x}\in\partial B,\quad\nu\times\nabla\times\mathcal{A}^{k}_{B}[\Phi]\Big{|}_{\pm}(\mathbf{x})=\lim_{t\rightarrow 0^{+}}\nu\times\nabla\times\mathcal{A}^{k}_{B}[\Phi](\mathbf{x}\pm t\nu),$$

$${\mathcal{M}}_B^k[\Phi ](\mathbf{x})=\text{p.v.}\nu \times \nabla \times {\int }_{\partial B}{\Gamma }_k(\mathbf{x}-\mathbf{y})\Phi (\mathbf{y})d{s}_{\mathbf{y}}.$$

$${\mathcal{L}}_B^k[\Phi ](\mathbf{x}):={\nu }_{\mathbf{x}}\times \nabla \times \nabla {\mathcal{A}}_B^k[\Phi ](\mathbf{x})={\nu }_{\mathbf{x}}\times (k^2{\mathcal{A}}_B^k[\Phi ](\Phi )+\nabla \nabla \cdot {\mathcal{A}}_B^k[\Phi ]),$$

$${\mathbf{H}}_0=\{\begin{array}{ll}-i{\omega }^{-1}(\nabla \times \nabla \times {\mathcal{A}}_{\Sigma }^{k_0}[{\Phi }_c]+\nabla \times \nabla \times {\mathcal{A}}_{\Sigma }^{k_0}[{\Phi }_s]+{\omega }^2{ϵ}_0\nabla \times {\mathcal{A}}_{\Sigma }^{k_0}[{\Psi }_s])& \text{in }{\mathbb{R}}^3\setminus \overline{\Sigma },\\ -i{\omega }^{-1}(\nabla \times \nabla \times {\mathcal{A}}_{\Sigma }^{k_s}[{\Phi }_c]+\nabla \times \nabla \times {\mathcal{A}}_{\Sigma }^{k_s}[{\Phi }_s]+{\omega }^2{ϵ}_0\nabla \times {\mathcal{A}}_{\Sigma }^{k_s}[{\Psi }_s])& \text{in }{\Sigma }_s,\end{array}$$

$${\mathbf{H}}_0=\nabla S_{{\Sigma }_c}^0(\frac{I}{2}+({\mathcal{K}}_{{\Sigma }_c}^0{)}^{\ast }{)}^{-1}[\nu \cdot {\mathbf{H}}_0{|}_{\partial {\Sigma }_c}^{+}]+\mathcal{O}(\omega )\text{in }{\mathbb{R}}^3\setminus {\overline{\Sigma }}_c.$$

Full text

Keywords: Geomagnetic monitoring, magnetic anomaly detection, varying anomalies, unique recovery

2010 Mathematics Subject Classification:  35Q60, 35J05, 31B10, 35R30, 78A40

1. Introduction

1.1. Mathematical setup

Focusing on the mathematics, but not the physics, we first present the mathematical setup of the inverse problem of identifying magnetic anomalies with varying parameters using geomagnetic monitoring. We define the Earth as a core-shell structure with $\Sigma_{c}$ and $\Sigma$ signifying the core and the Earth, respectively. It is assumed that both $\Sigma_{c}$ and $\Sigma$ are bounded on simply connected $C^{2}$ domains in $\mathbb{R}^{3}$ and $\Sigma_{c}\Subset\Sigma$. Then $\Sigma_{s}:=\Sigma\setminus\overline{\Sigma_{c}}$ signifies the shell of the Earth. We do not restrict the Earth and its core as concentric balls in our study. Throughout we let $\mathbf{x}=(x_{j})^{3}_{j=1}\in\mathbb{R}^{3}$ denote the space variable. We characterize the electromagnetic (EM) medium by the electric permittivity $\epsilon(\mathbf{x})$, magnetic permeability $\mu(\mathbf{x})$ and electric conductivity $\sigma(\mathbf{x})$ in $\mathbb{R}^{3}$. Let $\epsilon(\mathbf{x}),\ \mu(\mathbf{x})$ and $\sigma(\mathbf{x})$ be all real-valued scalar $L^{\infty}$ functions, such that $\epsilon(\mathbf{x})$ and $\mu(\mathbf{x})$ are positive, and $\sigma(\mathbf{x})$ is nonnegative. Let $\epsilon(\mathbf{x})=\epsilon_{0},\ \mu(\mathbf{x})=\mu_{0}$ and $\sigma(\mathbf{x})=0$ for $\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Sigma}$, where the two positive constants $\epsilon_{0}$ and $\mu_{0}$ signify the EM parameters of the uniformly homogeneous free space $\mathbb{R}^{3}\setminus\overline{\Sigma}$. Summarizing the above discussion, the medium distribution is described by

[TABLE]

where and also in what follows, $\chi$ denotes the characteristic function. (1.1) characterizes the material configuration of the Earth’s initial status.

Suppose there is a collection of magnetic anomalies which are growing or decaying in the shell of the Earth. Let the anomalies be represented as $D_{l}\subset\mathbb{R}^{3},l=1,2,\dots,l_{0},$ where $D_{l},1\leq l\leq l_{0}$ is a simply-connected Lipschitz domain such that the corresponding material parameters are described by $\epsilon_{l},\ \mu_{l}$ and $\sigma_{l}$. It is assumed that $\epsilon_{l},\ \mu_{l}$ and $\sigma_{l}$ are all positive constants with $\mu_{l}\neq\mu_{0},1\leq l\leq l_{0}$. From a practical point of view, the size of the magnetic anomaly $(D_{l};\epsilon_{l},\mu_{l},\sigma_{l}),1\leq l\leq l_{0}$ is much smaller than the size of the Earth. Hence, we can assume that the magnetic anomalies, with varying sizes, have the following forms

[TABLE]

where $\Omega$ is a bounded Lipschitz domain in $\mathbb{R}^{3}$ with $\Omega\Subset\Sigma$, and $\delta\in\mathbb{R}_{+}$ with $\delta\ll 1$. $\mathbf{z}_{l}\in\mathbb{R}^{3},l=1,2,\dots,l_{0}$ signify the positions of the magnetic anomalies. Moreover, $s_{l}\in\mathbb{R}_{+},l=1,2,\dots,l_{0}$ signify the variation parameters of the magnetic anomalies such that

[TABLE]

with $\alpha_{l}>-1,l=1,2,\dots,l_{0}$. Note that $\alpha_{l}>0$ means that the anomaly $D_{l}$ is decaying with a magnitude of $1+\alpha$, and $\alpha<0$ means that the anomaly $D_{l}$ is growing with a magnitude of $1+\alpha$, $l=1,2,\ldots,l_{0}$. Throughout the rest of the paper, we assume that the distribution of $D_{l},l=1,2,\dots,l_{0}$ is sparse, and $\mathbf{z}_{l},l=1,2,\dots,l_{0}$ are far away from $\partial\Sigma$ such that $\mathbf{x}-\mathbf{z}_{l}\gg s_{l}\delta,l=1,2,\dots,l_{0},$ for any $\mathbf{x}\in\partial\Sigma$. With the presence of the magnetic anomalies $(D_{l};\epsilon_{l},\mu_{l},\sigma_{l}),$ $l=1,2,\dots,l_{0}$ in the shell of the Earth, the EM medium configuration in $\mathbb{R}^{3}$ is then described by

[TABLE]

Let $\mathbf{E}(\mathbf{x})$ and $\mathbf{H}(\mathbf{x})$ respectively signify the Earth’s electric and magnetic fields, which satisfy the following Maxwell system in the time-harmonic regime for $(\mathbf{x},\omega)\in\mathbb{R}^{3}\times\mathbb{R}_{+}$ (cf. [5, 8]),

[TABLE]

In (1.5), $\omega\in\mathbb{R}_{+}$ is the frequency, $\hat{\rho}(\mathbf{x})\in L^{2}(\Sigma_{c})$ stands for the charge density of the Earth core $\Sigma_{c}$, and $\mathbf{v}\in L^{\infty}(\Sigma)$ is the fluid velocity of the Earth $\Sigma$. Here, $\hat{\rho}$ and $\mathbf{v}$ have the zero extension to $\mathbb{R}^{3}$. $\mathbf{v}\times\mu\mathcal{H}$ is the so-called motional electromotive force generated by the rotation of the Earth $\Sigma$. The last equation in (1.5) is called the Silver-Müller radiation condition and it holds uniformly in all directions $\hat{\mathbf{x}}:=\frac{\mathbf{x}}{\|\mathbf{x}\|}\in\mathbb{S}^{2}$ for $\mathbf{x}\in\mathbb{R}^{3}\setminus\{0\}$. The Silver-Müller radiation condition characterizes the outgoing electromagnetic waves (cf.[13]), and also guarantees the well-posedness of a physical solution to (1.5). There exists a unique pair of solutions $(\mathbf{E},\mathbf{H})\in H_{loc}(\mathrm{curl},\mathbb{R}^{3})^{2}$ to (1.5) (cf. [14]). Here and also in what follows, for any bounded domain $\mathcal{M}\subset\mathbb{R}^{3}$ with a Lipschitz boundary $\partial\mathcal{M}$,

[TABLE]

and

[TABLE]

Define $k_{0}:=\omega\sqrt{\epsilon_{0}\mu_{0}}$ to be the wavenumber with respect to a frequency $\omega\in\mathbb{R}_{+}$. Then the far-field pattern is the asymptotic expansion of the corresponding scattered electric field or magnetic field as $\|\mathbf{x}\|\rightarrow+\infty$ (cf. [6, 16])

[TABLE]

and

[TABLE]

which hold uniformly for all directions $\hat{\mathbf{x}}$.

Let $(\mathbf{E},\mathbf{H})$ be the solution to the Maxwell system (1.5) associated to the medium configuration in (1.4) with the variation parameters $\alpha_{l}=0,l=1,2,\dots,l_{0}$, and $(\mathbf{E}_{s},\mathbf{H}_{s})$ be the solution to (1.5) associated to (1.4) with the variation parameters $\alpha_{l}\neq 0,l=1,2,\dots,l_{0}$. That is, $(\mathbf{E},\mathbf{H})$ is the geomagnetic field of the Earth before the variation of the magnetic anomalies, whereas $(\mathbf{E}_{s},\mathbf{H}_{s})$ is the geomagnetic field after the change of the anomalies. Let $\Gamma$ be an open surface located away from $\Sigma$, and it signifies the location of the receivers for the geomagnetic monitoring.

With the above preparations, the inverse problem for the magnetic anomaly detection can be reformulated as

[TABLE]

where \widehat{\Gamma}:=\big{\{}\hat{\mathbf{x}}\in\mathbb{S}^{2};\hat{\mathbf{x}}=\frac{\mathbf{x}}{\|\mathbf{x}\|},\mathbf{x}\in\Gamma\}. That is, by monitoring the change of Earth’s geomagnetic field, one intends to recover the locations, material parameters as well as the size variations of the anomalies. By Rellich’s theorem, there is a one-to-one correspondence between the pair of the far-field patterns $(\mathbf{E}^{\infty},\mathbf{H}^{\infty})$ and the EM wave $(\mathbf{E},\mathbf{H})$ (cf. [6] ), and similar for the electromagnetic fields $(\mathbf{E}_{s}^{\infty},\mathbf{H}_{s}^{\infty})$ and $(\mathbf{E}_{s},\mathbf{H}_{s})$, respectively. By the real analyticity of $\mathbf{H}^{\infty}$ and $\mathbf{H}_{s}^{\infty}$, the inverse problem (1.10) is equivalent to the following one,

[TABLE]

1.2. Physical background and connection to existing results

Geomagnetic detections have been successfully applied in science and engineering including military submarine detection, minerals searching and geoexploration (cf. [20]). In fact, the operations of the Swarm satellites, a set of three probes sent by the European Space Agency (ESA) in the year 2013, are capable of measuring variations in the Earth’s magnetic field. These so-called secular variations can be detected by studying what happens beneath the Earth.

In geophysics, one widely accepted but relatively simpler model of Earth’s geomagnetic field is described by a time-dependent Maxwell model as follows. Let $\mathcal{E}(\mathbf{x},t)$ and $\mathcal{H}(\mathbf{x},t),(\mathbf{x},t)\in\mathbb{R}^{3}\times\mathbb{R}_{+}$, denote the electric and magnetic fields of the Earth $\Sigma$, and they satisfy the following time-dependent Maxwell system for $(\mathbf{x},t)\in\mathbb{R}^{3}\times\mathbb{R}_{+}$:

[TABLE]

where $\rho(\mathbf{x},t)=\rho_{c}(\mathbf{x},t)\chi(\Sigma_{c})\in H^{1}(L^{2}(\Sigma_{c}),\mathbb{R}_{+})$ stands for the charge density of the Earth core $\Sigma_{c}$, and $\mathbf{v}(\mathbf{x})\in L^{\infty}(\Sigma)$ is the fluid velocity of the Earth $\Sigma$.

Let $(\mathcal{E},\mathcal{H})$ be the geomagnetic field of the Earth before the variation of the magnetic anomalies, whereas $(\mathcal{E}_{s},\mathcal{H}_{s})$ be the geomagnetic field after the change of the anomalies. That is, $(\mathcal{E},\mathcal{H})$ is the solution to (1.12) associated to the medium configuration in (1.4) with the variation parameters $\alpha_{l}=0,l=1,2,\dots,l_{0}$, and $(\mathbf{E}_{s},\mathbf{H}_{s})$ is the solution to (1.12) associated to (1.4) with the variation parameters $\alpha_{l}\neq 0,l=1,2,\dots,l_{0}$ In such a setup, the magnetic anomaly detection can be formulated as follows,

[TABLE]

By introducing the following temporal Fourier transform for $(\mathbf{x},\omega)\in\mathbb{R}^{3}\times\mathbb{R}_{+}$:

[TABLE]

and setting

[TABLE]

one can show by direct verifications that (1.12) is reduced to (1.5), and accordingly, the inverse problem (1.13) is converted into (1.10); we also refer to [8] for more relevant discussion about this reduction.

In what follows, we would like to make a few remarks regarding the mathematical formulation of the inverse problem (1.13), or equivalently (1.11). First, we emphasize that in (1.11) or (1.13), we do not assume the medium configuration of the Earth’s core, the charge density and the fluid velocity to be a priori known. This make our study more practical, but more challenging on the other hand. Second, the measurement surface $\Gamma$ in (1.13) is in a scale much smaller than the Earth, and we are mainly concerned with the region where the corresponding geomagnetic variation can reach $\Gamma$. Hence, the rest part of the Earth’s medium configuration can be assumed to be the same to that of the aforementioned region. Therefore, it is unobjectionable to claim that the medium configuration in (1.4) is a realistic one. Third, the variation of the magnetic anomalies is assumed to take place in a very long time scale, i.e., the variation process occurs very slowly. In fact, one of the practical scenarios that motivates the current study is the discovery of a giant underground iron river in 2016 by observations of the change of Earth’s geomagnetic field; see [15, 17, 18]. Hence, one can have the Fourier transforms of the geomagnetic fields $\mathcal{H}_{s}$ and $\mathcal{H}$ separately. By a similar justification, the infinite time interval for the measurement in (1.13) can be replaced by a finite one $[0,T_{0}]$ for some finite $T_{0}\in\mathbb{R}_{+}$, since the outgoing geomagnetic variation eventually leaves the Earth after a finite time.

The magnetic anomaly detection was mathematically investigated [9, 8] by two of the authors of the present article. In [8], a Maxwell system is used for the the geomagnetic model, whereas in [9], a magnetohydrodynamic system is used for the geomagnetic model. Both [8] and [9] deal with the case that the geomagnetic variation is caused by the presence of certain magnetized anomalies, which e.g. corresponds to a specific practical scenario of submarine detection. This is in a sharp difference to the present study, where the geomagnetic variation is caused by the change of the magnetic anomalies. That is, the magnetic anomaly already exists, and its growth or shrinkage causes the geomagnetic variation. In order to establish the unique recovery results in such a new situation, we basically follow the general strategy developed in [9, 8] by combining low-frequency asymptotics and various linearization techniques. However, there are several new challenges that require some new technical developments. It is emphasized that in (1.11), we do not assume that the Earth’s geomagnetic field in the case with no anomaly is known a priori. Hence, one cannot naively reduce the current study to the one considered in [9, 8], where the knowledge of the free geomagnetic field is a critical ingredient. Finally, we would like to point out that our uniqueness argument is constructive, and it actually implies a construction procedure for the anomaly detection.

The rest of the paper is organized as follows. In Section 2, we derive the low-frequency asymptotic expansions of the geomagnetic fields as well as further asymptotic expansions of the static fields with respect to the variation parameters of the anomalies. We basically refine the asymptotic arguments in [8] by further analyzing the linearized fields. Section 3 is devoted to the recovery of the free geomagnetic field. In Section 4, we establish the main unique recovery results on identifying the positions, the variation parameters as well as the material parameters of the magnetic anomalies.

2. Asymptotic expansions of the geomagnetic fields

For our subsequent study of the inverse problem (1.10), or equivalently (1.11), we derive in this section several critical asymptotic expansions of the geomagnetic fields $\mathbf{H}$ and $\mathbf{H}_{s}$, respectively, associated with (1.4)–(1.5). These include the low-frequency asymptotic expansions of the geomagnetic fields in obtaining the corresponding static fields, and further asymptotic expansions of the static fields with respect to the anomaly sizes. We make essential use of tools from the layer potential theory, which we shall first collect in the following.

2.1. Layer potentials

We recall some preliminary knowledge on the layer potential techniques for the subsequent use. Let $B\subset\mathbb{R}^{3}$ be a bounded Lipschitz domain. We introduce some function spaces on the boundary $\partial B$. Let us first define the space of tangential vector fields by

[TABLE]

We denote by $\nabla_{\partial B}\cdot$ the standard surface divergence operator defined on $L^{2}_{T}(\partial B)$. Let $H^{s}(\partial B)$ be the usual Sobolev space of order $s\in\mathbb{R}$ on $\partial B$. Define the normed spaces of tangential fields by

[TABLE]

endowed with the norms

[TABLE]

respectively.

Let $k\in\mathbb{R}_{+}$ to be a wavenumber. We recall that the outgoing fundamental solution $\Gamma_{k}$ to the PDO $(\Delta+k^{2})$ in $\mathbb{R}^{3}$ is given by

[TABLE]

For a scalar density $\phi\in L^{2}(\partial B)$, let $S_{B}^{k}:H^{-1/2}(\partial B)\rightarrow H^{1}(\mathbb{R}^{3}\setminus\partial B)$ be the single layer potential associated with $B\subset\mathbb{R}^{3}$,

[TABLE]

and $\mathcal{K}_{B}^{k}:H^{1/2}(\partial B)\rightarrow H^{1/2}(\partial B)$ be the Neumann-Poincaré operator

[TABLE]

where p.v. stands for the Cauchy principle value. In (2.6) and also in what follows, unless otherwise specified, $\nu$ signifies the exterior unit normal vector to the boundary of the concerned domain. The single layer potential $S^{k}_{B}$ satisfies the following trace formula

[TABLE]

where $(\mathcal{K}^{k}_{B})^{*}$ is the adjoint operator of $\mathcal{K}^{k}_{B}$ with respect to $L^{2}$ inner product.

In addition, for a density $\Phi\in\mathrm{TH}(\mathrm{div},\partial B)$, we define the vectorial single layer potential by

[TABLE]

It is known that $\nabla\times\mathcal{A}^{k}_{B}$ satisfies the jump formula (see [1, 16])

[TABLE]

where

[TABLE]

is understood in the sense of uniform convergence on $\partial B$ and the boundary integral operator $\mathcal{M}^{k}_{B}:\text{TH(div,}\partial B)\rightarrow\text{TH(div,}\partial B)$ is given by

[TABLE]

We also define $\mathcal{L}^{k}_{B}$ : TH (div,$\partial B)\rightarrow$ TH (div,$\partial B$) by

[TABLE]

where the formula $\nabla\times\nabla\times=-\Delta+\nabla\nabla\cdot$ is used and shall also be frequently used in what follows.

It is mentioned that $\frac{I}{2}\pm\mathcal{M}^{k}_{B}$ is invertible on TH(div,$\partial B$) when $k$ is sufficiently small (see e.g. [1, 19]). In the following , if $k=0$, we formally set $k$ introduced in (2.4) to be $\Gamma_{0}=-\frac{1}{4\pi\|\mathbf{x}\|}$, and the other integral operators introduced above can also be formally defined when $k=0$. Finally, we define $k_{s}:=\omega\sqrt{\mu_{0}\epsilon_{s}}$.

2.2. Low-frequency asymptotic expansions of the magnetic fields

Through out the rest of the paper, we let $(\mathbf{E}_{0},\mathbf{H}_{0})$ be the solution to (1.1) and (1.5). That is, $\mathbf{E}_{0}$ and $\mathbf{H}_{0}$ are the free/background geomagnetic field of the Earth. In this subsection, we derive the low-frequency asymptotic expansions of the magnetic fields $\mathbf{H}_{0}$ and $\mathbf{H}_{s}$. In the following, we define the wavenumbers $\zeta_{l},l=1,2,\dots,l_{0}$ by $\zeta_{l}^{2}:=\omega^{2}\mu_{l}\gamma_{l},\gamma_{l}:=\epsilon_{l}+i\sigma_{l}/\omega$, where the sign of $\zeta_{l}$ is chosen such that $\Im\zeta_{l}\geq 0$. For the sake of simplicity, we set $\tilde{D}:=\Sigma_{s}\setminus\overline{\bigcup_{l{{}^{\prime}}=1}^{l_{0}}{D_{l{{}^{\prime}}}}}$. We also let $\nu_{l}$ be the exterior unit normal vector defined on $\partial D_{l},l=1,2,\dots,l_{0}$.

We next deal with the low-frequency asymptotic analysis of the geomagnetic system (1.5). As remarked earlier, we shall basically follow the general strategy developed in [8], but there are also several new technical developments in the current setup of out study. We would also like to mention in passing some related literature on the low-frequency asymptotic analysis of the Maxwell system [3, 4, 7, 8, 9, 10, 11, 12].

We first deal with $\mathbf{H}_{0}$, which can be represented by the following integral ansatz,

[TABLE]

where $(\Phi_{c},\Phi_{s},\Psi_{s})\in\text{TH(div,}\partial\Sigma_{c})\times\text{TH(div,}\partial\Sigma_{c})\times\text{TH(div,}\partial\Sigma_{c})$. There holds the following asymptotic expansion,

Lemma 2.1**.**

(Lemma 2.3 in [8]) Let $\mathbf{H}_{0}$ be the solution to (1.1) and (1.5). Then for $\omega\in\mathbb{R}_{+}$ sufficiently small, one has

[TABLE]

Next we establish the asymptotic expansion of the magnetic field $\mathbf{H}$ associated with the system (1.4) and (1.5). By using the transmission conditions across $\partial\Sigma$ and $\partial D_{l},l=1,2,\dots,l_{0}$ , $(\mathbf{E}_{s},\mathbf{H}_{s})$ is the solution to the following transmission problem:

[TABLE]

where $[\nu\times\mathbf{E}_{s}]$ and $[\nu\times\mathbf{H}_{s}]$ denote the jumps of $\nu\times\mathbf{E}_{s}$ and $\nu\times\mathbf{H}_{s}$ along $\partial\Sigma$ and $\partial{D_{l}}$, namely,

[TABLE]

Using the potential theory, the solution $(\mathbf{E}_{s},\mathbf{H}_{s})$ can be represented by the following formulas

[TABLE]

and

[TABLE]

where by transmission conditions across $\partial\Sigma$ and $\partial D_{l},l=1,2,\dots,l_{0}$ , $(\Phi_{0},\Psi_{0})\in\text{TH(div,}\partial\Sigma)\times\text{TH(div},\partial\Sigma)$ and $(\Phi_{l},\Psi_{l})\in\text{TH(div,}\partial D_{l})\times\text{TH(div},\partial D_{l}),l=1,2,\dots,l_{0}$ satisfy

[TABLE]

with the operators $\mathcal{P}^{k,k^{\prime}}_{D,D^{\prime}}$, $\mathcal{D}^{k,k^{\prime}}_{D,D^{\prime}}$ and $\mathcal{F}^{k,k^{\prime}}_{D,D^{\prime}}$ defined by

[TABLE]

and the operators $\mathcal{M}^{k}_{D,D^{\prime}}$ and $\mathcal{L}^{k}_{D,D^{\prime}}$ defined by

[TABLE]

where $k,k^{\prime}\in\{k_{0},k_{s},\zeta_{1},\zeta_{2},\dots,\zeta_{l_{0}}\},D,D{{}^{\prime}}\in\{\Sigma,D_{1},D_{2},\dots,D_{l_{0}}\}$ and $\mu^{(k)},\gamma^{(k)}$ are respectively the parameters $\mu$ and $\gamma$, which are related to $k$. In other words,

[TABLE]

The fields $(\hat{\mathbf{E}}_{0},\hat{\mathbf{H}}_{0})$ in (2.17) and (2.18) satisfy (1.4) and (1.5) in $(\mathbb{R}^{3}\setminus\overline{\Sigma})\bigcup\tilde{D}$ with $\hat{\mathbf{H}}_{0}$ given in Lemma 2.1 in [8] (see Lemma 2.2 below) and they depend on the background fields $({\mathbf{E}}_{0},{\mathbf{H}}_{0})$ and the boundary condition on $\partial\Sigma_{c}$.

Based on the above integral representation, we present the asymptotic formula of the magnetic field $\mathbf{H}_{s}$ with respect to the frequency. To that end, we first define an $l_{0}$-by-$l_{0}$ matrix type operator $\mathbb{K}^{*}_{D}$ on $L^{2}(\partial D_{1})\times L^{2}(\partial D_{2})\times\cdots\times L^{2}(\partial D_{l_{0}})$ by

[TABLE]

Define $l_{0}$-by-$l_{0}$ matrices type operators $\mathbb{M}_{D}$ and $\mathbb{N}_{D}$ on $\mathrm{TH(div},\partial D_{1})\times\cdots\times\mathrm{TH(div},\partial D_{l_{0}})$ by

[TABLE]

where $l,m\in\{1,2,\dots,l_{0}\}$ and

[TABLE]

Theorem 2.1**.**

Let $(\mathbf{E}_{s},\mathbf{H}_{s})$ be the solution to (1.4) and (1.5). Then for $\omega\in\mathbb{R}_{+}$ sufficiently small, there hold the following asymptotic expansions:

[TABLE]

where $\Theta,\Xi\in\mathrm{TH(div,}\partial\Sigma)$ satisfy

[TABLE]

and $\Psi^{(0)}_{l}\in\mathrm{TH(div,}\partial D_{l}),l=1,2,\dots,l_{0}$ satisfy

[TABLE]

where $\mathbb{L}^{\gamma}_{D}$ is defined by

[TABLE]

with the operator $\mathbb{M}_{D}$ and $\mathbb{N}_{D}$ defined in (2.24) and (2.25), respectively. Here the notation $\otimes$ stands for the Kronecker product and $\bm{e}_{l}$ is the $l_{0}$-dimensional Euclidean unit vector. $\Pi_{l}\in L^{2}(\partial D_{l}),l=1,2,\dots,l_{0}$ are defined by

[TABLE]

where $\mathbb{J}^{\mu}_{D}$ is defined on $L^{2}(\partial D_{1})\times L^{2}(\partial D_{2})\times\cdots\times L^{2}(\partial D_{l_{0}})$ given by

[TABLE]

with the operator $\mathbb{K}^{*}_{D}$ defined in (2.23) and $\mathbf{C}$ is defined by

[TABLE]

The parameter $\lambda_{\gamma_{l}}$, $\lambda_{\mu_{l}}$ and $\lambda_{\epsilon}$ are defined by

[TABLE]

and

[TABLE]

Proof.

The proof follows from that of Theorem 3.1 in [8] with some straightforward modifications. ∎

Lemma 2.2**.**

(Lemma 3.3 in [8]) $\hat{\mathbf{H}}_{0}$ introduced in (2.26) can be written as

[TABLE]

where $\varphi_{l},l=1,2,\dots,l_{0}$ satisfy

[TABLE]

with the operator $\mathbb{J}^{\epsilon}_{D}$ defined in (2.31) with $\lambda_{\mu_{l}}$ replaced by $\lambda_{\epsilon_{l}},l=1,2,\dots,l_{0}.$

2.3. Asymptotic expansions of $\mathbf{H}_{s}$ with respect to the anomaly size

It is noted that if $\sigma_{l},l=1,2,\dots,l_{0}$ are not identically zero, then the leading-order term in (2.26) may still depend on $\omega$. In the subsequent analysis, we make use of the leading-order term in the asymptotic low-frequency expansion in the representation formula (2.26) to further study the asymptotic expansion with respect to the size of the magnetic anomalies.

We first present the following lemma.

Lemma 2.3**.**

Suppose $D_{l},l=1,2,\dots,l_{0}$ are defined in (1.2) with $\delta\in\mathbb{R}_{+}$ sufficiently small and the corresponding variation parameters $s_{l}=\delta^{\alpha_{l}},l=1,2,\dots,l_{0}$ with $\alpha_{l}>-1$. Then we have

[TABLE]

where $\mathbb{K}^{*}_{\Omega},\mathbb{M}_{\Omega}$ are $l_{0}\times l_{0}$ matrix-valued operators defined by

[TABLE]

respectively. $\mathbb{K}^{*}_{D}$ and $\mathbb{M}_{D}$ are defined in (2.23) and (2.24), respectively.

Proof.

We only prove the second assertion in (2.37), and the first one can be proved in a similar manner. For any $\mathbf{x},\mathbf{y}\in\partial D_{l}$, we let $\tilde{\mathbf{x}}=(s_{l}\delta)^{-1}(\mathbf{x}-\mathbf{z}_{l}),\tilde{\mathbf{y}}=(s_{l}\delta)^{-1}(\mathbf{y}-\mathbf{z}_{l})\in\partial\Omega,l=1,2,\dots,l_{0}$. Define $\tilde{\Phi}(\tilde{\mathbf{y}})=\Phi(\mathbf{y})$. By using change of variables, one can show that there holds

[TABLE]

On the other hand, letting $\mathbf{x}\in\partial D_{m}$ and $\mathbf{y}\in\partial D_{l}$ and $\tilde{\mathbf{x}}=(s_{m}\delta)^{-1}(\mathbf{x}-\mathbf{z}_{m}),\tilde{\mathbf{y}}=(s_{l}\delta)^{-1}(\mathbf{y}-\mathbf{z}_{l})\in\partial\Omega$, where $l,m\in\{1,2,\dots,l_{0}\}$ and $l\neq m$, one can show that

[TABLE]

By Taylor expansion there holds

[TABLE]

Thus, for the case when $\alpha_{m}>\alpha_{l}$,

[TABLE]

By direct calculations, we can also verify the cases when $\alpha_{m}\leq\alpha_{l}$. Hence,

[TABLE]

By substituting (2.39) and (2.43) back into (2.24), one readily has (2.37).

The proof is complete. ∎

We introduce some notations for further analysis. In what follows, we let $\mathbf{H}^{(0)}_{0}$ be the leading-order term of $\mathbf{H}_{0}$ in (2.26), $\hat{\mathbf{H}}^{(0)}_{0}$ be the leading-order term of $\hat{\mathbf{H}}_{0}$ (and $\varphi_{l}^{(0)}$ is the leading-order term of $\varphi_{l},l=1,2,\dots,l_{0}$). Similarly, we let $\mathbf{H}^{(0)}$ and $\mathbf{H}^{(0)}_{s}$ be the leading-order term of $\mathbf{H}$ and $\mathbf{H}_{s}$, respectively. Recall the following important result: .

Lemma 2.4**.**

(Lemma 3.6 in [8]) For any simply connected domain $D_{l}$ and a conservative gradient field $\mathbf{A}$ in $\mathbb{R}^{3}$, three holds the following relation:

[TABLE]

By [8] one can find that the leading-order term $\mathbf{H}^{(0)}_{0}$ is a conservative gradient field. With the above preparations, we are in a position to derive the first main result, i.e., the asymptotic expansion of the geomagnetic field $\mathbf{H}_{s}$ with respect to the size of the magnetic anomalies.

Theorem 2.2**.**

Suppose $D_{l},l=1,2,\dots,l_{0}$ are defined in (1.2) with $\delta\in\mathbb{R}_{+}$ sufficiently small and $s_{l}=\delta^{\alpha_{l}},l=1,2,\dots,l_{0}$ with $\alpha_{l}>-1$. Let $(\mathbf{E}_{s},\mathbf{H}_{s})$ be the solution to (1.4) and (1.5). The for $\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Sigma}$, there holds the following asymptotic expansion,

[TABLE]

where $\mathbf{P}_{0}$ is defined by

[TABLE]

The polarization tensors $\mathbf{D}_{l}$ and $\bm{M_{l}}$ are $3\times 3$ matrices defined by

[TABLE]

and

[TABLE]

respectively, $l=1,2,\dots,l_{0}$. More specifically, setting $\mathbf{M}_{l}=((\mathbf{D}_{l})_{mn}),m,n=1,2,3,\tilde{\mathbf{y}}=(\tilde{\mathbf{y}}_{1},\tilde{\mathbf{y}}_{2},\tilde{\mathbf{y}}_{3})^{T}$ and $\nu_{l}=(\nu_{l}^{(1)},\nu_{l}^{(2)},\nu_{l}^{(3)})^{T}$, we have

[TABLE]

$\mathbf{P}_{0}$* and $\mathbf{M}_{l}$ have similar forms.*

Proof.

We note that either $\omega\gamma_{l}$ or $1/(\gamma_{l}-\epsilon_{l})$ is of order $\omega,l=1,2,\dots,l_{0}$, no matter $\sigma$ is zero or nonzero. One can immediately find that the second term in (2.30) is of order $\omega$. By (2.26), it then can be seen that the leading order term $\mathbf{H}^{(0)}_{s}$ has the following form.

[TABLE]

where $\Psi^{\prime}_{l}$ and $\Pi^{\prime}_{l}$ are defined by

[TABLE]

and

[TABLE]

respectively, $l=1,2,\dots,l_{0}$. It can be verified that

[TABLE]

Thus, $\nabla\times\mathcal{A}^{0}_{D_{l}}[\Psi^{\prime}_{l}]$ is a gradient field of harmonic function in $\mathbb{R}^{3}\setminus\overline{D_{l}}$. By using the relations

[TABLE]

and

[TABLE]

in $\mathbb{R}^{3}\setminus\overline{\Sigma}$, (see Lemma 3.2 in [8]) and (2.53), one can derive that

[TABLE]

As before, for $\mathbf{y}\in\partial D_{l}$, we let $\tilde{\mathbf{y}}=(s_{l}\delta)^{-1}(\mathbf{y}-\mathbf{z}_{l})\in\partial\Omega$, and define $\tilde{\Psi}^{\prime}_{l}:=\Psi^{\prime}_{l}(\mathbf{y}),\tilde{\Pi}^{\prime}_{l}:=\Pi^{\prime}_{l}(\mathbf{y}),l=1,2,\dots,l_{0}$ and $\tilde{\hat{\mathbf{H}}}^{(0)}_{0}(\tilde{\mathbf{y}}):=\hat{\mathbf{H}}^{(0)}_{0}(\mathbf{y}).$ Then by Lemma 2.3, one has

[TABLE]

and

[TABLE]

$l=1,2,\dots,l_{0}.$ Hence, by using Lemma 2.4, there holds

[TABLE]

On the other hand, by the Taylor expansion, there holds

[TABLE]

and so by using (2.35), (2.58) and (2.7), one has

[TABLE]

where $\tilde{\varphi}_{l}^{(0)}(\tilde{\mathbf{y}}):=\varphi_{l}^{(0)}(\mathbf{y})$ and $\varphi_{l}(\mathbf{y})=\varphi_{l}^{(0)}(\mathbf{y})+O(\omega)$. By (2.36) and (2.58) one has

[TABLE]

Thus, by combining (2.59), (2.60) and (2.59) one can derive that

[TABLE]

For $\mathbf{x}-\mathbf{z}_{l}\gg s_{l}\delta$ there holds

[TABLE]

Define $\tilde{Q}_{l}(\tilde{\mathbf{y}}):=Q_{l}(\mathbf{y})$, where $\tilde{Q}_{l}$ is given in (2.57). By using change of variables and (2.35), one has

[TABLE]

Finally, by using the definitions (2.46)-(2.48) and substituting (2.55)-(2.62) into (2.54), one obtains

[TABLE]

The first equatily of (2.63) is obtain by using the following facts:

[TABLE]

To prove (2.65), we first set

[TABLE]

By using the jump formula (2.7) and integration by parts, one can show that there holds

[TABLE]

We only prove the first assertion in (2.65), and the second one can be proved in a similar manner. The proof is complete. ∎

For notational convenience, in the sequel, we introduce the matrix $\mathbf{P}_{l}$ by

[TABLE]

where $\mathbf{M}_{l}$ and $\mathbf{D}_{l}$ are defined in (2.48) and (2.47), respectively.

Remark 2.1*.*

We remark that the matrix $\mathbf{P}_{l}$ is exactly same as the matrix obtained in [8]. Hence, it satisfies the following nonsingularity properties:

Lemma 2.5** (Lemma 3.7 in [8]).**

If $\sigma_{l}\neq 0,l=1,2,\dots,l_{0}$ and $\epsilon_{s}=\epsilon_{0}$, then $\mathbf{P}_{l}=\mu_{0}\mathbf{M}_{l}+O(\omega)$ is nonsingular.

Lemma 2.6** (Lemma 3.8 in [8]).**

Suppose $\Omega$ is a ball. Let $\mathbf{P}_{l}$ be defined in (2.68). If there holds

[TABLE]

then $\mathbf{P}_{l}$ is nonsingular.

Let $\mathbf{H}_{s}^{(0)}$ be the leading term of $\mathbf{H}_{s}$. Then from (2.72), there holds, for $\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Sigma},$

[TABLE]

On the other hand, we let $\mathbf{H}^{(0)}$ be the leading term of $\mathbf{H}$, with variation parameters $\alpha_{l}=0,l=1,2,\dots,l_{0}$. Then there holds, for $\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Sigma},$

[TABLE]

By combining (2.70) and (2.71) one thus has the following result:

Theorem 2.3**.**

Suppose $D_{l},l=1,2,\dots,l_{0}$ are defined in (1.2) with $\delta\in\mathbb{R}_{+}$ sufficiently small and $s_{l}=\delta^{\alpha_{l}},l=1,2,\dots,l_{0}$ with $-1/4<\alpha_{l}<1/3$. Let $(\mathbf{E}_{s},\mathbf{H}_{s})$ be the solution to (1.4) and (1.5). Then for $\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Sigma}$, there holds the following asymptotic expansion,

[TABLE]

Remark 2.2*.*

Two remarks are in order. First, the asymptotic formula (2.72) describes the leading term of the difference of the magnetic field before and after the growth ($-1/4<\alpha_{l}<0$) or shrinkage ($0<\alpha_{l}<1/3$) of the magnetic anomalies, $l=1,2,\ldots,l_{0}$. The terms $\mathbf{P}_{l}\mathbf{H}_{0}^{(0)}(\mathbf{z}_{l})$ are not possible to calculate or measure directly, and we use the measurement of perturbation of magnetic fields to estimate the terms in the subsequent section. Second, the restriction $-1/4<\alpha_{l}<1/3$ is to ensure that $3(1+\alpha_{l})<4$ and $4(1+\alpha_{l})>3$, which guarantees that the leading order term in (2.72) can be distinguished from the higher order term for all $l=1,2,\ldots,l_{0}$. As an illustrative example, we suppose there are two anomalies, say $l_{1}$ and $l_{2}$, with one is growing ($l_{1}<0$) and the other is decaying ($l_{2}>0$). If there holds that $4(1+\alpha_{l_{1}})=3(1+\alpha_{l_{2}})$, or $3\alpha_{l_{2}}-4\alpha_{l_{1}}=1$, then the higher order term of the anomaly $l_{1}$ is the same as the leading order term of the anomaly $l_{2}$, which means that $l_{2}$ can not be distinguished, and thus the shrinkage of the anomaly $l_{2}$ can not be recovered by only using the measurement (1.13).

3. Further asymptotic analysis

In this section, we make further analyse on the leading-order term of the asymptotic expansion in (2.72). As remarked in Remark 2.2, we show that the auxiliary terms $\mathbf{P}_{l}\mathbf{H}_{0}^{(0)}(\mathbf{z}_{l})$ can be well approximated by the measurement data in (1.10). This shall be of critical importance if one intends to develop a practical reconstruction scheme by using our unique recovery results in the next section. We shall mainly make use of the orthogonality of spherical harmonics functions. Before that, we first present some axillary results.

Lemma 3.1**.**

(Lemma 3.9 in [8]) Let $\mathbf{z}\in\mathbb{R}^{3}$ be fixed. Let $\mathbf{x}\in\partial B_{R}$ and suppose $\|\mathbf{z}\|<R$. There holds the following asymptotic expansion

[TABLE]

where $\hat{\mathbf{z}}=\mathbf{z}/\|\mathbf{z}\|$ and $\hat{\mathbf{x}}=\mathbf{x}/\|\mathbf{x}\|$. $Y^{m}_{n}$ is the spherical harmonics of order $m$ and degree $n$.

Lemma 3.2**.**

(Equation (3.9) in [9]) Let $\mathbf{z}_{l}\in\mathbb{R}^{3}$ be fixed. Let $\mathbf{x}\in\partial B_{R}$ and suppose $\|\mathbf{z}\|<R$. There holds the following asymptotic expansion

[TABLE]

where

[TABLE]

and $\hat{\mathbf{z}}=\mathbf{z}/\|\mathbf{z}\|$ and $\hat{\mathbf{x}}=\mathbf{x}/\|\mathbf{x}\|$. $Y^{m}_{n}$ is the spherical harmonics of order $m$ and degree $n$.

Using the above results, we can further estimate the leading order term in (2.72). By substituting (3.2) into (2.72), one has

[TABLE]

In the sequel, we define

[TABLE]

for $n\in\mathbb{N}\cup{0}$ and $m=-n,-n+1,\dots,n-1,n,$ and

[TABLE]

for $n\in\mathbb{N}$ and $m=-n,-n+1,\dots,n-1,n.$ Note that $\mathbf{N}^{m}_{n+1}(\hat{\mathbf{x}})$, $\mathbf{Q}^{m}_{n-1}(\hat{\mathbf{x}})$ and $\mathbf{T}^{m}_{n}(\hat{\mathbf{x}})$ are spherical harmonics of order $n$. We can further establish the following lemmas.

Lemma 3.3**.**

Let $\mathbf{P}_{l}$ be defined in (2.68). Then there holds

[TABLE]

where $\mathbf{C}_{0}$ and $\mathbf{N}_{2}(\hat{\mathbf{x}})$ are 3-by-3 matrices given by

[TABLE]

and

[TABLE]

respectively. $\mathbf{a}^{m^{\prime},m}_{n^{\prime},n}$ is defined in (3.12).

Proof.

From (3.3), one has

[TABLE]

where $\xi\in\mathbb{R}^{3}$. Note that $\mathbf{A}^{m}_{n}(\hat{\mathbf{x}})\xi$ is a linear combination of the spherical harmonics $\{\mathbf{N}_{n+1}^{m}(\hat{\mathbf{x}})\}$ and $\{\mathbf{Q}_{n-1}^{m}(\hat{\mathbf{x}})\}$. We show that $\mathbf{N}_{2}^{m^{\prime}}(\hat{\mathbf{x}}),m^{\prime}=-1,0,1$ are orthogonal to $\mathbf{A}_{n}^{m}(\hat{\mathbf{x}})\xi$ for $\xi\in\mathbb{R}^{3}$ and $n\neq 0,2$, where $\mathbf{A}_{n}^{m}(\hat{\mathbf{x}})$ is defined in (3.3). By vector calculus identity and straightforward computations, one can obtain that

[TABLE]

where

[TABLE]

and

[TABLE]

By using the following elementary result (cf. [16])

[TABLE]

one finally obtains

[TABLE]

where

[TABLE]

and

[TABLE]

By using (3.10) the orthogonality of the vectorial spherical harmonics, one has

[TABLE]

Note that if $n=2$, then by $\eqref{Axi}$, one has

[TABLE]

where $\mathbf{c}_{n^{\prime},n}^{m^{\prime},m}$ and $\mathbf{d}_{n^{\prime},n}^{m^{\prime},m}$ are defined in (3.16) and (3.17), respectively. Thus,

[TABLE]

By using the orthogonality of $\mathbf{N}_{3}^{m}(\hat{\mathbf{x}})$ and $Y_{1}^{m^{\prime}},m,m^{\prime}\in\{-1,0,1\}$, one has

[TABLE]

which indicates $\mathbf{a}_{1,2}^{m^{\prime},m}=3\mathbf{b}_{1,2}^{m^{\prime},m}$, where $\mathbf{a}_{n^{\prime},n}^{m^{\prime},m}$ and $\mathbf{b}_{n^{\prime},n}^{m^{\prime},m}$ are defined in (3.12) and (3.13), respectively. Together with (3.20) and (3.16), one has

[TABLE]

Similarly, for $n=0$, one obtains that

[TABLE]

Thus,

[TABLE]

By using the orthogonality of $\mathbf{N}_{1}^{m}(\hat{\mathbf{x}})$ and $Y_{1}^{m^{\prime}},m,m^{\prime}\in\{-1,0,1\}$, one has

[TABLE]

which indicates $\mathbf{a}_{1,0}^{m^{\prime},0}=\mathbf{b}_{1,0}^{m^{\prime},0}$. Combining (3.24) and (3.16), one has

[TABLE]

Finally, by taking the inner product of (3.4) and $\mathbf{N}_{2}^{m^{\prime}}(\hat{\mathbf{x}})$ in $L^{2}(\mathbb{S})$ and by using (3.22) and (3.23), one has

[TABLE]

By rearranging the terms in equation (3.27), we finally arrive at (3.7), which completes the proof. ∎

Lemma 3.4**.**

Let $P_{l}$ be defined in (2.68). Then there holds

[TABLE]

where $\mathbf{D}_{0}$ and $\mathbf{Q}_{0}(\hat{\mathbf{x}})$ are 3-by-3 matrices given by

[TABLE]

and

[TABLE]

respectively. $\mathbf{a}^{m^{\prime},m}_{n^{\prime},n}$ is defined in (3.12).

Proof.

The proof of (3.28) follows from a similar argument to that in the proof of Lemma 3.3. We show that $\mathbf{Q}_{0}^{m^{\prime}}(\hat{\mathbf{x}}),m^{\prime}=-1,0,1$ are orthogonal to $\mathbf{A}_{n}^{m}(\hat{\mathbf{x}})\xi$ for $\xi\in\mathbb{R}^{3}$ and $n\neq 0,2$, where $\mathbf{A}_{n}^{m}(\hat{\mathbf{x}})$ is defined in (3.3). By vector calculus identity and straightforward computations, one can obtain that

[TABLE]

where $\mathbf{a}^{m^{\prime},m}_{n^{\prime},n}$ and $\mathbf{b}^{m^{\prime},m}_{n^{\prime},n}$ are defined in (3.12) and (3.13) respectively. By using (3.15) the orthogonality of the vectorial spherical harmonics, one has

[TABLE]

Note that if $n=2$, then by $\eqref{AxiCD}$, one has

[TABLE]

where $\mathbf{d}_{n^{\prime},n,}^{m^{\prime},m}$ is defined in (3.17). Together with (3.21), one has

[TABLE]

Similarly, for $n=0$, one obtains that

[TABLE]

Using (3.25) and (3.17), one has

[TABLE]

Finally, by taking the inner product of (3.4) and $\mathbf{Q}_{0}^{m^{\prime}}(\hat{\mathbf{x}})$ in $L^{2}(\mathbb{S})$ and by using (3.22) and (3.23), one has

[TABLE]

By rearranging the terms in equation (3.27), we finally arrive at (3.7). The proof is complete. ∎

The following expressions are the general forms of Lemmas 3.3 and 3.4 respectively. By (3.11) and taking the inner product of (3.4) and (3.5) in $L^{2}(\mathbb{S})$, one can obtain that

[TABLE]

where $\mathbb{S}$ stands for the unit sphere. Similarly, one can also obtain that

[TABLE]

We remark that the coefficient of higher order $\mathcal{O}(\|\mathbf{x}\|^{-2})$ in (3.7) is accurately given in (3.27), and the general form is given in (3.38). Similar argument can be applied to the spherical harmonics $\mathbf{Q}_{n-1}^{m}$. If the quantity $\|\mathbf{z}_{l}\|/\|\mathbf{x}\|$ in (3.27) is small enough, that is, the radius of the measurement surface is much bigger than the length of the position vector of the magnetized anomaly, then the unknown term $\mathbf{P}_{l}\mathbf{H}_{0}^{(0)}(\mathbf{z}_{l})$ is easy to recover. In the following, we consider a particular case with $\Gamma=\partial B_{R}$, with $B_{R}$ a sufficiently large central ball containing $\Sigma$.

Lemma 3.5**.**

Suppose that $\Gamma=\partial B_{R}$ and let $\mathbf{P}_{l}$ be defined in (2.68). Then there hold the following relationships,

[TABLE]

and

[TABLE]

Proof.

The proof of (3.40) and (3.41) follows from a similar argument to that in the proof of Lemma 3.3 in [9]. By straightforward calculations, one can show that $\mathbf{C}_{0}$ and $\mathbf{D}_{0}$ are invertible. The proof is completed by directly using (3.7). ∎

4. Unique recovery results for magnetic anomaly detections

In this section, we present the main unique recovery results in identifying the magnetic anomalies for the inverse problem (1.11), which are contained in Theorems 4.1 and 4.2. In what follows, we let $D_{l}^{(1)}$ and $D_{l}^{(2)},l=1,2,\dots,l_{0},$ be two sets of magnetic anomalies, which satisfy (1.2) with $\mathbf{z}_{l}$ replaced by $\mathbf{z}_{l}^{(1)}$ and $\mathbf{z}_{l}^{(2)}$, respectively, and $s_{l}=\delta^{\alpha}_{l}$ replaced by $s_{l}^{(1)}=\delta^{\alpha^{(1)}_{l}}$ and $s_{l}^{(2)}=\delta^{\alpha^{(2)}_{l}}$, respectively. The material parameters $\epsilon_{l},\sigma_{l},\gamma_{l}$ and $\mu_{l}$ are also replaced by $\epsilon_{l}^{(1)},\sigma_{l}^{(1)},\gamma_{l}^{(1)},\mu_{l}^{(1)}$ and $\epsilon_{l}^{(2)},\sigma_{l}^{(2)},\gamma_{l}^{(2)},\mu_{l}^{(2)}$, for $D_{l}^{(1)}$ and $D_{l}^{(2)},l=1,2,\dots,l_{0}$ respectively. Let $\mathbf{H}_{s,j},j=1,2$ be the solutions to (1.4) and (1.5) with $D_{l}$ replaced by $D_{l}^{(1)}$ and $D_{l}^{(2)}$, respectively. Let $\mathbf{H}_{j},j=1,2$ be the solutions to (1.1) and (1.5) with $D_{l}$ replaced by $D_{l}^{(1)}$ and $D_{l}^{(2)}$, respectively. Let $\mathbf{H}_{s,j}^{(0)}$ be the leading-order terms of $\mathbf{H}_{s,j}$ with respect to $\omega\ll 1$. Let $\mathbf{H}_{s}^{(0)}$ be the leading-order terms of $\mathbf{H}_{s}$ with respect to $\omega\ll 1$. Denote by $\mathbf{D}_{l}^{(1)},\mathbf{D}_{l}^{(2)},\mathbf{M}_{l}^{(1)},\mathbf{M}_{l}^{(2)},\mathbf{P}_{l}^{(1)}$ and $\mathbf{P}_{l}^{(2)}$ the polarization tensors for $D_{l}^{(1)}$ and $D_{l}^{(2)}$, respectively, $l=1,2,\dots,l_{0}$.

Before that, we first present an important axillary result.

Lemma 4.1**.**

Suppose $D_{l},l=1,2,\dots,l_{0}$ are defined in (1.2) with $\delta\in\mathbb{R}_{+}$ sufficiently small and $s_{l}=\delta^{\alpha_{l}},l=1,2,\dots,l_{0}$ with $-1/4<\alpha_{l}<1/3$. Suppose also that $3(\alpha_{l_{1}}+1)\neq 4(\alpha_{l_{2}}+1)$ for any $l_{1},l_{2}\in\{1,2,\ldots,l_{0}\}$. Let $(\mathbf{E}_{s},\mathbf{H}_{s})$ be the solution to (1.4) and (1.5). If there hold

[TABLE]

and

[TABLE]

then one has

[TABLE]

for any $n\in\mathbb{N}\cup\{0\}$, where

[TABLE]

Here, $\hat{\mathbf{z}}_{l}^{(j)}=\mathbf{z}_{l}^{(j)}/\|\mathbf{z}_{l}^{(j)}\|,l=1,2\dots,l_{0},j=1,2$ , $\hat{\mathbf{x}}=\mathbf{x}/\|\mathbf{x}\|$ and $Y_{n}^{m}$ is the spherical harmonics of order $m$ and degree $n$.

Proof.

First, by using (4.1) and (4.2) and unique continuation, one has

[TABLE]

and

[TABLE]

Then from (2.70), one has

[TABLE]

Suppose $R\in\mathbb{R}_{+}$ is sufficiently large such that $\Sigma\Subset B_{R}$. By (4.7) one readily has

[TABLE]

On the other hand, it can be verified that

[TABLE]

are harmonic functions in $\mathbb{R}^{3}\setminus\overline{B_{R}}$, which decay at infinity. Using this together with (4.8), and the maximum principle of harmonic functions, one can obtain that

[TABLE]

and therefore

[TABLE]

By substituting (3.1) and (4.11) one has

[TABLE]

By taking $\|\mathbf{x}\|$ sufficiently large and comparing the order of $\|\mathbf{x}\|^{-1}$ one has (4.3). The proof is complete. ∎

4.1. Uniqueness in recovering a single anomaly

We first present the uniqueness result in recovering a single anomaly.

Theorem 4.1**.**

Suppose $l_{0}=1$ and $\sigma_{1}^{(1)},\sigma_{1}^{(2)}\neq 0$. If there hold (4.1) and (4.2) and $\alpha_{1}^{(j)}>-1$, $j=1,2$, then $\mathbf{z}_{1}^{(1)}=\mathbf{z}_{1}^{(2)},\ \alpha_{1}^{(1)}=\alpha_{1}^{(2)}$ and $\mu_{1}^{(1)}=\mu_{1}^{(2)}$.

Proof.

By the unique continuation principle, one has from (4.1) that $\mathbf{H}_{s,1}=\mathbf{H}_{s,2}$ in $\mathbb{R}^{3}\setminus D_{1}^{(1)}\cup D_{1}^{(2)}$. First, we let $n=0$ and $l_{0}=1$ and note that $Y_{0}^{0}=\frac{1}{2\sqrt{\pi}}$. By using (4.4) and (3.5), one has

[TABLE]

Hence by using (4.3), there holds

[TABLE]

which readily implies that

[TABLE]

or in the form with respect to $\alpha_{l}^{(1)}$ and $\alpha_{l}^{(2)}$

[TABLE]

In what follows, we define $\mathbf{C}_{s,1}^{(j)}:=s_{1}^{(j)}\mathbf{P}_{1}^{(j)}\mathbf{H}_{0}^{(0)}(\mathbf{z}_{1}^{(j)}),j=1,2$. Then $\mathbf{C}_{s,1}^{(1)}=\mathbf{C}_{s,1}^{(2)}$. We claim that

[TABLE]

By using Lemma 2.3 in [8], one has $\mathbf{H}_{0}^{(0)}=\nabla u_{0}$, where $u_{0}$ is a non-constant harmonic function. Hence, by the maximum principle of harmonic functions, $\nabla u_{0}(\mathbf{z}_{1}^{(j)})\neq 0,j=1,2$. This together with the assumption that $\mathbf{P}_{1}^{(j)},j=1,2,$ are nonsingular matrices and $\delta\in\mathbb{R}_{+}$, one has $\mathbf{C}_{s,1}^{(j)}\neq 0,j=1,2.$

Now by setting $n=1$, we define

[TABLE]

Let $(\hat{\mathbf{x}},\bm{e}_{\phi},\bm{e}_{\theta})$ be the triplet of the local orthogonal unit vectors, where $\phi$ and $\theta$ depend on $\hat{\mathbf{x}}$. Note that

[TABLE]

By using (4.3), one has

[TABLE]

where $\bm{Q}(\hat{\mathbf{x}})$ is defined by

[TABLE]

Straightforward calculations show that $\bm{Q}(\hat{\mathbf{x}})$ is nonsingular for any $\hat{\mathbf{x}}\in\mathbb{S}^{2}$. Since $\hat{\mathbf{x}}\in\mathbb{S}^{2}$ is arbitrarily given, (4.19) implies that

[TABLE]

Then by direct calculations, one has

[TABLE]

To recover the variation parameters, by the unique continuation principle again, one sees that from (4.2) $\mathbf{H}_{1}=\mathbf{H}_{2}\in\mathbb{R}^{3}\setminus D_{1}^{(1)}\cup D_{1}^{(2)}$. We let $n=0$ and use Lemma 4.1 in [8], there holds

[TABLE]

which readily implies that

[TABLE]

By substituting (4.24) in (4.15) or (4.16) one has

[TABLE]

We thus have $D_{1}^{(1)}=D_{1}^{(2)}$. Hence, in the sequel, we let $D_{1}:=D_{1}^{(1)}=D_{1}^{(2)}$. Clearly, $\mathbf{H}_{s,1}=\mathbf{H}_{s,2}$ in $\mathbb{R}\setminus\overline{D_{1}}$. Since $\sigma_{1}^{(1)},\sigma_{1}^{(2)}\neq 0$, from (2.26), one can see that

[TABLE]

By using the jump formula one further has that

[TABLE]

Using the fact that $\frac{I}{2}+(\mathcal{K}^{0}_{D_{1}})^{*}$ is invertible on $L^{2}(\partial D_{1})$ and some elementary calculations, one sees from (4.27) that

[TABLE]

Note that $\nu\cdot\hat{\mathbf{H}}_{0}^{(0)}\in L^{2}_{0}(\partial D_{1})$, where $L^{2}_{0}(\partial D_{1})$ is a subset of $L^{2}(\partial D_{1})$ with zero average on $\partial D_{1}$. Since $\nu_{1}\cdot\hat{\mathbf{H}}_{0}^{(0)}\neq 0$, thus from (2.35) that $\nu_{1}\cdot\hat{\mathbf{H}}_{0}^{(0)}\neq 0$ on $\partial D_{1}$. Using this and the fact that $-\frac{I}{2}+(\mathcal{K}^{0}_{D_{1}})^{*}$ is invertible on $L^{2}_{0}(\partial D_{1})$, we finally from (4.28) that $\mu_{1}^{(1)}=\mu_{1}^{(2)}$. The proof is complete. ∎

Remark 4.1*.*

We remark that for a single anomaly case, we do not need the assumption that $-1/4<\alpha_{1}^{(j)}<1/3$ but only need $\alpha_{1}^{(j)}>-1$, $j=1,2$. In fact, one can easily see from (2.70) that the expansions of different orders with respect to $\delta$ can be distinguished exactly for one single anomaly.

4.2. Uniqueness in recovering multiple anomalies

Theorem 4.2**.**

If there hold (4.1) and (4.2) and there satisfy $-1/4<\alpha_{l}^{(j)}<1/3$ and $3(\alpha_{l_{1}}^{(j)}+1)\neq 4(\alpha_{l_{2}}^{(j)}+1)$, $j=1,2$, $l=1,2,\dots,l_{0}$ for any $l_{1},l_{2}\in\{1,2,\ldots,l_{0}\}$, then $\mathbf{z}_{l}^{(1)}=\mathbf{z}_{l}^{(2)}$, $\alpha_{l}^{(1)}=\alpha_{l}^{(2)}$ and $\mu_{l}^{(1)}=\mu_{l}^{(2)}$, $l=1,2,\dots,l_{0}$.

Proof.

With our earlier preparatiosn, the proof follows from a similar argument to that of Theorem 7.8 in [2] and a similar idea to that of Theorem 4.2 in [8]. In the following, we only sketch it. Using the formula (4.8) and analysis similar to the proof of Lemma 4.1, one can show that

[TABLE]

By straightforward calculations, one can further show that

[TABLE]

where $\mathbf{z}_{l}^{\prime}=\mathbf{z}_{l}^{(1)}+t^{\prime}\mathbf{z}_{l}^{(2)}$ with $t^{\prime}\in(0,1)$. Note that $F^{s}(\mathbf{x})$ defined in (4.30) is also harmonic in $\mathbb{R}^{3}\setminus\bigcup_{l=1}^{l_{0}}(\mathbf{z}_{l}^{(1)}\cup\mathbf{z}_{l}^{(2)})$. By using the analytic continuation of harmonic functions, one thus has that $F^{s}(\mathbf{x})\equiv 0\in\mathbb{R}^{3}$. Define $F^{s}:=F_{1}^{s}+F_{2}^{s}$, where

[TABLE]

and

[TABLE]

Then by comparing the types of poles of $F_{1}^{s}$ and $F_{2}^{s}$, one immediately finds that $F_{1}^{s}=0$ and $F_{2}^{s}=0$ in $\mathbb{R}^{3}$. Since $\mathbf{H}_{0}^{(0)}(\mathbf{x})$ does not vanish for $\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Sigma}$, then one has $s_{l}^{(1)}\mathbf{P}_{l}^{(1)}\mathbf{H}_{0}^{(0)}(\mathbf{z}_{l}^{(1)})\neq 0$. Hence one sees that

[TABLE]

By using (4.32) and $F^{s}_{2}=0$, one has

[TABLE]

Again, by similar argument to above and using the proof of Theorem 4.2 in [8],

[TABLE]

Hence we have

[TABLE]

Thus by unique continuation one further has

[TABLE]

To prove that $\mu_{l}^{(1)}=\mu_{l}^{(2)}$, one only needs to apply the following transmission conditions:

[TABLE]

for $j=1,2$ and use (4.37). The proof then follows similarly to [9]. ∎

5. Concluding Remark

In this paper, we consider the identification of magnetic anomalies with varing parameters beneath the Earth by using geomagnetic monitoring. The work is an extension of the existing results in [8] on identification of magnetized anomalies (size and material parameters). In this paper, the magnetic anomalies are supposed to grow or shrink when time elapses and our main aim is to identify the anomalies (size, material parameters, varying scale) by observation of the geomagnetic variation–the so-called secular variation. In order to avoid the use of the background geomagnetic field, which is not assumed to be known from a practical viewpoint, new technique strategy is developed in this paper. Our results provide a rigorous mathematical theory to some existing applications in the geophysics community.

Acknowledgment

The work of Y. Deng was supported by NSF grant of China No. 11601528, NSF grant of Hunan No. 2017JJ3432 and No. 2018JJ3622, PSCF and RFEB of Hunan No. 18YBQ077 and No. 18B337, Innovation-Driven Project of Central South University, No. 2018CX041. The work of H. Liu was supported by the FRG and startup grants from Hong Kong Baptist University, Hong Kong RGC General Research Funds, 12302017 and 12301218.