IDENTIFYING VARYING MAGNETIC ANOMALIES USING GEOMAGNETIC MONITORING

. We are concerned with the inverse problem of identifying magnetic anomalies with varying parameters beneath the Earth using geomagnetic monitoring. Observations of the change in Earth’s magnetic ﬁeld–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 ex- isting results in [9] by two of the authors to the more practical and challenging scenario with varying anomalies.


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 Σ c and Σ signifying the core and the Earth, respectively. It is assumed that both Σ c and Σ are bounded on simply connected C 2 domains in R 3 and Σ c ⊂⊂ Σ. Then Σ s := Σ \ Σ 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 x = (x j ) 3 j=1 ∈ R 3 denote the space variable. We characterize the electromagnetic (EM) medium by the electric permittivity (x), magnetic permeability µ(x) and electric conductivity σ(x) in R 3 . Let (x), µ(x) and σ(x) be all real-valued scalar L ∞ functions, such that (x) and µ(x) are positive, and σ(x) is nonnegative. Let (x) = 0 , µ(x) = µ 0 and σ(x) = 0 for x ∈ R 3 \ Σ, where the two positive constants 6412 YOUJUN DENG, HONGYU LIU AND WING-YAN TSUI 0 and µ 0 signify the EM parameters of the uniformly homogeneous free space R 3 \Σ. Summarizing the above discussion, the medium distribution is described by where and also in what follows, χ denotes the characteristic function. (1) characterizes the material configuration of the Earth's initial status.
Suppose there is a collection of magnetic anomalies that the sizes of magnetic anomalies become bigger or smaller after some time in the shell of the Earth. Let the anomalies be represented as D l ⊂ R 3 , l = 1, 2, . . . , l 0 , where D l , 1 ≤ l ≤ l 0 is a simply-connected Lipschitz domain such that the corresponding material parameters are described by l , µ l and σ l . It is assumed that l , µ l and σ l are all positive constants with µ l = µ 0 , 1 ≤ l ≤ l 0 . From a practical point of view, the size of the magnetic anomaly (D l ; l , µ l , σ l ), 1 ≤ l ≤ 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 D l = s l δΩ + z l , l = 1, 2, . . . , l 0 , where Ω is a bounded Lipschitz domain in R 3 with Ω ⊂⊂ Σ, and δ ∈ R + with δ 1. In (2), z l ∈ R 3 , l = 1, 2, . . . , l 0 signify the positions of the magnetic anomalies. Moreover, s l ∈ R + , l = 1, 2, . . . , l 0 signify the variation parameters of the magnetic anomalies such that s l = δ α l , l = 1, 2, . . . , l 0 , with α l > −1, l = 1, 2, . . . , l 0 . Note that α l > 0 means that the anomaly D l is decaying with a magnitude of 1 + α l , l = 1, 2, . . . , l 0 , and α l < 0 means that the anomaly D l is growing with a magnitude of 1 + α l , l = 1, 2, . . . , l 0 . In (2), Ω is referred to as the reference domain, which signifies the shape of the anomaly. In this paper, we are mainly concerned with recovering the positions, varying parameters as well as the material parameters of the anomalies by knowing their shapes. It is noted that in (2), we assume that the reference domain remains the same for all the anomalies D l , l = 1, 2, . . . , l 0 , and this mainly serves to ease the exposition of our subsequent study. We would like to emphasize the reference domains can also vary for different anomalies, and this shall become more evident in our subsequent analysis. In fact, one can see from (59) in what follows, the reference domain mainly contributes to the underlying polarization tensors (PTs), and as long as those PTs are a-priori known, all of our results hold for the more general case with varying but a-priori known reference domains. Nevertheless, we stick to the case with the anomalies given in (2) throughout the rest of the paper. Moreover, we assume that D l , l = 1, 2, . . . , l 0 are sparsely distributed and far away from each other and z l , l = 1, 2, . . . , l 0 are far away from ∂Σ such that x−z l s l δ, l = 1, 2, . . . , l 0 , for any x ∈ ∂Σ. With the presence of the magnetic anomalies (D l ; l , µ l , σ l ), l = 1, 2, . . . , l 0 in the shell of the Earth, the EM medium configuration in R 3 is then described by σ(x) =σ c (x)χ(Σ c ) + Let E(x) and H(x) respectively signify the Earth's electric and magnetic fields, which satisfy the following Maxwell system in the time-harmonic regime for (x, ω) ∈ R 3 × R + (cf. [5,9]), In (5), ω ∈ R + is the frequency,ρ(x) ∈ L 2 (Σ c ) stands for the charge density of the Earth core Σ c , and v ∈ L ∞ (Σ c ) 3 signifies a velocity filed in the Earth core.
Here and also in what follows,ρ and v have zero extensions to R 3 . v × µH is the so-called motional electromotive force generated by the rotation of the Earth Σ. The last equation in (5) is called the Silver-Müller radiation condition and it holds uniformly in all directionsx := x x ∈ S 2 for x ∈ R 3 \ {0}. The Silver-Müller radiation condition characterizes the outgoing electromagnetic waves (cf. [14]), and also guarantees the well-posedness of a physical solution to (5). There exists a unique pair of solutions (E, H) ∈ H loc (curl, R 3 ) 2 to (5) (cf. [15]). Here and also in what follows, for any bounded domain M ⊂ R 3 with a Lipschitz boundary ∂M, Define k 0 := ω √ 0 µ 0 to be the wavenumber with respect to a frequency ω ∈ R + . Then the far-field pattern is the asymptotic expansion of the corresponding scattered electric field or magnetic field as x → +∞ (cf. [7,17]) and which hold uniformly for all directionsx. Let (E, H) be the solution to the Maxwell system (5) associated to the medium configuration in (4) with the variation parameters α l = 0, l = 1, 2, . . . , l 0 , and (E s , H s ) be the solution to (5) associated to (4) with the variation parameters α l = 0, l = 1, 2, . . . , l 0 . That is, (E, H) is the geomagnetic field of the Earth before the variation of the magnetic anomalies, whereas (E s , H s ) is the geomagnetic field after the change of the anomalies. Let Γ be an open surface located away from Σ, 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 where Γ := x ∈ S 2 ;x = x x , x ∈ Γ}. 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 (E ∞ , H ∞ ) and the EM wave (E, H) (cf. [7] ), and similar for the electromagnetic fields (E ∞ s , H ∞ s ) and (E s , H s ), respectively. By the real analyticity of H ∞ and H ∞ s , the inverse problem (10) is equivalent to the following one, 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 E(x, t) and H(x, t), (x, t) ∈ R 3 × R + , denote the electric and magnetic fields of the Earth Σ, and they satisfy the following time-dependent Maxwell system for (x, t) ∈ R 3 × R + : where ρ(x, t) = ρ c (x, t)χ(Σ c ) ∈ H 1 (L 2 (Σ c ), R + ) stands for the charge density of the Earth core Σ c , and v(x) ∈ L ∞ (Σ c ) 3 describes the velocity filed of the fluid in the Earth core.
Let (E, H) be the geomagnetic field of the Earth before the variation of the magnetic anomalies, whereas (E s , H s ) be the geomagnetic field after the change of the anomalies. That is, (E, H) is the solution to (12) associated to the medium configuration in (4) with the variation parameters α l = 0, l = 1, 2, . . . , l 0 , and (E s , H s ) is the solution to (12) associated to (4) with the variation parameters α l = 0, l = 1, 2, . . . , l 0 In such a setup, the magnetic anomaly detection can be formulated as follows, By introducing the following temporal Fourier transform for (x, ω) ∈ R 3 × R + : and setting one can show by direct verifications that (12) is reduced to (5), and accordingly, the inverse problem (13) is converted into (10); we also refer to [9] 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 (13), or equivalently (11). First, we emphasize that in (11) or (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 Γ in (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 Γ. 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 (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 [16,19,6]. Hence, one can have the Fourier transforms of the geomagnetic fields H s and H separately. By a similar justification, the infinite time interval for the measurement in (13) can be replaced by a finite one [0, T 0 ] for some finite T 0 ∈ R + , since the outgoing geomagnetic variation eventually leaves the Earth after a finite time.
The magnetic anomaly detection was mathematically investigated [10,9] by two of the authors of the present article. In [9], a Maxwell system is used for the geomagnetic model, whereas in [10], a magnetohydrodynamic system is used for the geomagnetic model. Both [9] and [10] 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 [10,9] 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 (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 [10,9], 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 lowfrequency 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. 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 (10), or equivalently (11), we derive in this section several critical asymptotic expansions of the geomagnetic fields H and H s , respectively, associated with (4)- (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 ⊂ R 3 be a bounded Lipschitz domain. We introduce some function spaces on the boundary ∂B. Let us first define the space of tangential vector fields by In (15) and also in what follows, unless otherwise specified, ν signifies the exterior unit normal vector to the boundary of the concerned domain. We denote by ∇ ∂B · the standard surface divergence operator defined on L 2 T (∂B). Let H s (∂B) be the usual Sobolev space of order s ∈ R on ∂B. Define the normed spaces of tangential fields by endowed with the norms Φ TH(div,∂B) = Φ L 2 (∂B) + ∇ ∂B · Φ L 2 (∂B) , respectively. Let k ∈ R + to be a wavenumber. We recall that the outgoing fundamental solution Γ k to the partial differential operator (∆ + k 2 ) in R 3 is given by For a scalar density φ ∈ L 2 (∂B), let S k and K k B : H 1/2 (∂B) → H 1/2 (∂B) be the Neumann-Poincaré operator where p.v. stands for the Cauchy principle value. The single layer potential operator S k B satisfies the following trace formula where (K k B ) * is the adjoint operator of K k B with respect to L 2 inner product. In addition, for a density Φ ∈ TH(div, ∂B), we define the vectorial single layer potential by It is known that ∇ × A k B satisfies the jump formula (cf. [1,17]) where is understood in the sense of uniform convergence on ∂B and the boundary integral operator M k B : TH(div,∂B) → TH(div,∂B) is given by We also define L k B : TH (div,∂B) → TH (div,∂B) by where the formula ∇ × ∇× = −∆ + ∇∇· is used and shall also be frequently used in what follows. It is mentioned that I 2 ± M k B is invertible on TH(div,∂B) when k is sufficiently small (cf. e.g. [1,18]). In the following , if k = 0, we formally set Γ k introduced in (18) to be Γ 0 = − 1 4π x , and the other integral operators introduced above can also be formally defined when k = 0. Finally, we define k s := ω √ µ 0 s .

2.2.
Low-frequency asymptotic expansions of the magnetic fields. Through out the rest of the paper, we let (E 0 , H 0 ) be the solution to (1) and (5). That is, E 0 and 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 H 0 and H s . In the following, we define the wavenumbers ζ l , l = 1, 2, . . . , l 0 by ζ 2 l := ω 2 µ l γ l , γ l := l + iσ l /ω, where the sign of ζ l is chosen such that ζ l ≥ 0. For the sake of simplicity, we setD := Σ s \ l0 l =1 D l . We also let ν l be the exterior unit normal vector defined on ∂D l , l = 1, 2, . . . , l 0 .
We next deal with the low-frequency asymptotic analysis of the geomagnetic system (5). As remarked earlier, we shall basically follow the general strategy developed in [9], 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,8,9,10,11,12,13].
We first deal with H 0 , which can be represented by the following integral ansatz, There holds the following asymptotic expansion, ) Let H 0 be the solution to (1) and (5). Then for ω ∈ R + sufficiently small, one has Next we establish the asymptotic expansion of the magnetic field H associated with the system (4) and (5). By using the transmission conditions across ∂Σ s and ∂D l , l = 1, 2, . . . , l 0 , (E s , H s ) is the solution to the following transmission problem: where [ν × E s ] and [ν × H s ] denote the jumps of ν × E s and ν × H s along ∂Σ and ∂D l , namely, Using the potential theory, the solution (E s , H s ) can be represented by the following formulas (cf. Lemma 3.1 of [9] for related discussions as well), and where the fields (Ê 0 ,Ĥ 0 ) in (31) and (32) satisfy (4) and (5) in (R 3 \ Σ) D witĥ H 0 given in Lemma 3.3 in [9] (see Lemma 2.3 below) and they depend on the background fields (E 0 , H 0 ) and the interface condition on ∂Σ c . Note that (Ê 0 ,Ĥ 0 ) do not depend on the parameter α introduced in (3).
By transmission conditions across ∂Σ and ∂D l , l = 1, 2, . . . , l 0 , with the operators P k,k D,D , D k,k D,D and F k,k D,D defined by and the operators M k D,D and L k D,D defined by where k, k ∈ {k 0 , k s , ζ 1 , ζ 2 , . . . , ζ l0 }, D, D ∈ {Σ, D 1 , D 2 , . . . , D l0 } and µ (k) , γ (k) are respectively the parameters µ and γ, which are related to k. In other words, Based on the above integral representation, we present the asymptotic formula of the magnetic field H s with respect to the frequency. To that end, we first define Define l 0 -by-l 0 matrices type operators M D and N D on TH(div, ∂D 1 ) × · · · × TH(div, ∂D l0 ) by where l, m ∈ {1, 2, . . . , l 0 } and We are in the position of representing the asymptotic expansion result of the magnetic field H s with respect to low frequency. (4) and (5). Then for ω ∈ R + sufficiently small, there hold the following asymptotic expansions: where Θ, Ξ ∈ TH(div,∂Σ) satisfy and Ψ (0) l ∈ TH(div,∂D l ), l = 1, 2, . . . , l 0 satisfy where L γ D is defined by with the operator M D and N D defined in (38) and (39), respectively. Here the notation ⊗ stands for the Kronecker product and e l is the l 0 -dimensional Euclidean unit vector. Π l ∈ L 2 (∂D l ), l = 1, 2, . . . , l 0 are defined by where J µ D is defined on L 2 (∂D 1 ) × L 2 (∂D 2 ) × · · · × L 2 (∂D l0 ) given by with the operator K * D defined in (37) and C is defined by The parameter λ γ l , λ µ l and λ are defined by and Proof. The proof follows from that of Theorem 3.1 in [9] with some straightforward modifications.

2.3.
Asymptotic expansions of H s ith respect to the anomaly size. It is noted that if σ l , l = 1, 2, . . . , l 0 are not identically zero, then the leading-order term in (40) may still depend on ω. In the subsequent analysis, we make use of the leading-order term in the asymptotic low-frequency expansion in the representation formula (40) to further study the asymptotic expansion with respect to the size of the magnetic anomalies. We stress that the conductivities (and other electromagnetic parameters) of the magnetic anomalies are finite and are independent of the frequency and size of the anomalies. We first present the following lemma.

YOUJUN DENG, HONGYU LIU AND WING-YAN TSUI
We next introduce some notations for the further analysis. In what follows, we let H ) For any simply connected domain D l and a conservative gradient field A in R 3 , three holds the following relation: By [9] one can find that the leading-order term 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 H s with respect to the size of the magnetic anomalies.

YOUJUN DENG, HONGYU LIU AND WING-YAN TSUI
The proof is complete.
For notational convenience, in the sequel, we introduce the matrix P l by where M l and D l are defined in (62) and (61), respectively.
Starting from now on and throughout the rest of the paper, we always assume that P l , l = 1, 2, . . . l 0 are nonsingular. Lemmas 2.7 and 2.8 indicate that this condition is non-void. Let H (0) s be the leading term of H s . Then from (59), there holds, for x ∈ R 3 \ Σ, (85) On the other hand, we let H (0) be the leading term of H, with variation parameters α l = 0, l = 1, 2, . . . , l 0 . Then there holds, for x ∈ R 3 \ Σ, By combining (85) and (86) one thus has the following result: Theorem 2.9. Suppose D l , l = 1, 2, . . . , l 0 are defined in (2) with δ ∈ R + sufficiently small and s l = δ α l , l = 1, 2, . . . , l 0 with −1/4 < α l < 1/3. Let (E s , H s ) be the solution to (4) and (5). Then for x ∈ R 3 \ Σ, there holds the following asymptotic expansion, Remark 3. Two remarks are in order. First, the asymptotic formula (87) describes the leading term of the difference of the magnetic field before and after the growth (−1/4 < α l < 0) or shrinkage (0 < α l < 1/3) of the magnetic anomalies, l = 1, 2, . . . , l 0 . The terms P l H (0) 0 (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 < α l < 1/3 is to ensure that 3(1 + α l ) < 4 and 4(1 + α l ) > 3, which guarantees that the leading order term in (87) can be distinguished from the higher order term for all l = 1, 2, . . . , l 0 . However, there are still cases that two different variational anomalies can not be distinguished. 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 + α l1 ) = 3(1 + α l2 ), or 3α l2 − 4α l1 = 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 (13).
3. Further asymptotic analysis. In this section, we make further analyse on the leading-order term of the asymptotic expansion in (87). As remarked in Remark 3, we show that the auxiliary terms P l H (0) 0 (z l ) can be well approximated by the measurement data in (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 [9]) Let z ∈ R 3 be fixed. Let x ∈ ∂B R and suppose z < R. There holds the following asymptotic expansion whereẑ = z/ z andx = x/ x . Y m n is the spherical harmonics of order m and degree n.
In the sequel, we define (cf. Theorem 2.4.7 [17]) for n ∈ N ∪ {0} and m = −n, −n + 1, . . . , n − 1, n, and for n ∈ N and m = −n, −n + 1, . . . , n − 1, n. Note that N m n+1 (x), Q m n−1 (x) and T m n (x) are spherical harmonics of order n. We can further establish the following lemmas. [10]) Let z l ∈ R 3 be fixed. Let x ∈ ∂B R and suppose z < R. There holds the following asymptotic expansion

Lemma 3.2. (Equation (3.9) in
where A m n := (n + 1)(x∇ s Y m andẑ = z/ z andx = x/ x . Y m n is the spherical harmonics of order m and degree n.

YOUJUN DENG, HONGYU LIU AND WING-YAN TSUI
Using the above results, we can further estimate the leading order term in (87). By substituting (91) into (87), one has Lemma 3.3. Let P l be defined in (83). Then there holds where C 0 and N 2 (x) are 3-by-3 matrices given by and respectively. a m ,m n ,n is defined in (99).
Lemma 3.4. Let P l be defined in (83). Then there holds where D 0 and Q 0 (x) are 3-by-3 matrices given by and respectively. a m ,m n ,n is defined in (99). Proof. The proof of (115) follows from a similar argument to that in the proof of where a m ,m n ,n and b m ,m n ,n are defined in (99) and (100) respectively. By using (102) the orthogonality of the vectorial spherical harmonics, one has (cf. (3.17) in [10]) Note that if n = 2, then by (102), one has where d m ,m n ,n, is defined in (104). Together with (108), one has Similarly, for n = 0, one obtains that Using (112) and (104), one has Finally, by taking the inner product of (93) and Q m 0 (x) in L 2 (S) and by using (109) and (110), one has (124) By rearranging the terms in equation (114), we finally arrive at (94). The proof is complete.
The following expressions are the general forms of Lemmas 3.3 and 3.4 respectively. By (98) and taking the inner product of (93) and (89) in L 2 (S), one can obtain that where S stands for the unit sphere. Similarly, one can also obtain that (126) We remark that the coefficient of higher order O( x −2 ) in (94) is accurately given in (114), and the general form is given in (125). Similar argument can be applied to the spherical harmonics Q m n−1 . If the quantity z l / x in (114) 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 P l H (0) 0 (z l ) is easy to recover. In the following, we consider a particular case with Γ = ∂B R , with B R a sufficiently large central ball containing Σ.
Lemma 3.5. Suppose that Γ = ∂B R and let P l be defined in (83). Then there hold the following relationships, and l0 l=1 (δ 3α l − 1)P l H where ≈ signifies the difference between the left-hand and right-hand terms is O(δ).
Proof. The proof of (127) and (128) follows from a similar argument to that in the proof of Lemma 3.3 in [10]. By straightforward calculations, one can show that C 0 and D 0 are invertible. The proof is completed by directly using (94).
Proof. First, by using (129) and (130) and unique continuation, one has and Then from (85), one has l0 l=1 s (1) l Suppose R ∈ R + is sufficiently large such that Σ ⊂⊂ B R . By (135) one readily has On the other hand, it can be verified that are harmonic functions in R 3 \ B R , which decay at infinity. Using this together with (136), and the maximum principle of harmonic functions, one can obtain that and therefore By taking x sufficiently large and comparing the order of x −1 one has (131). The proof is complete.
4.1. Uniqueness in recovering a single anomaly. We first present the uniqueness result in recovering a single anomaly. Proof. By the unique continuation principle, one has from (129) that H s,1 = H s,2 in R 3 \ D By using (131), one has where Q(x) is defined by Straightforward calculations show that Q(x) is nonsingular for anyx ∈ S 2 . Sincê x ∈ S 2 is arbitrarily given, (147) implies that Then by direct calculations, one has 1 .
To recover the variation parameters, by the unique continuation principle again, one sees that from (130) 1 . We let n = 0 and use Lemma 4.1 in [9], there holdŝ which readily implies that 1 .

YOUJUN DENG, HONGYU LIU AND WING-YAN TSUI
Remark 4. We remark that for a single anomaly case, we do not need the assumption that −1/4 < α  l , l = 1, 2, . . . , 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 [9]. In the following, we only sketch it. Using the formula (136) and analysis similar to the proof of Lemma 4.1, one can show that l0 l=1 s (1) By straightforward calculations, one can further show that where z l = z (1) l + t z (2) l with t ∈ (0, 1). Note that F s (x) defined in (158) is also harmonic in R 3 \ l0 l=1 (z (1) l ∪ z (2) l ). By using the analytic continuation of harmonic functions, one thus has that F s (x) ≡ 0 ∈ R 3 . Define F s : and Then by comparing the types of poles of F s 1 and F s 2 , one immediately finds that F s 1 = 0 and F s 2 = 0 in R 3 . Since H By using (160) and F s 2 = 0, one has s