BTTB–RRCG method for downward continuation of potential field data

Journal of Applied Geophysics 126 (2016) 74–86

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

BTTB–RRCG method for downward continuation of potential field data Yile Zhang a,⁎, Yau Shu Wong a, Yuanfang Lin b a b

Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, AB T6G 2G1, Canada Department of Mathematics, East China Normal University, Shanghai 200062, China

a r t i c l e

i n f o

Article history: Received 15 May 2015 Received in revised form 6 January 2016 Accepted 15 January 2016 Available online 19 January 2016 Keywords: Downward continuation Block-Toeplitz Toeplitz-Block Conjugate gradient Reweighted regularization

a b s t r a c t This paper presents a conjugate gradient (CG) method for accurate and robust downward continuation of potential field data. Utilizing the Block-Toeplitz Toeplitz-Block (BTTB) structure, the storage requirement and the computational complexity can be significantly reduced. Unlike the wavenumber domain regularization methods based on fast Fourier transform, the BTTB-based conjugate gradient method induces little artifacts near the boundary. The application of a re-weighted regularization in a space domain significantly improves the stability of the CG scheme for noisy data. The synthetic data with different levels of added noise and real field data are used to validate the effectiveness of the proposed scheme, and the computed results are compared with those based on recently proposed wavenumber domain methods and the Taylor series method. The simulation results verify that the proposed scheme is superior to the existing methods considered in this study in terms of accuracy and robustness. The proposed scheme is a powerful computational tool capable of applications for large scale data with modest computational cost. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Downward continuation is frequently applied to enhance the potential field data. It provides geological information at low elevation by using the field data from high elevation. In recent years, as the aerogravity and magnetic surveys become more widely used in prospecting (Zhang et al., 2015), it is desirable to develop efficient and robust downward continuation methods to deal with large amounts of aeropotential field data. According to the physical law, the potential field data at higher elevation contains dim geophysical information, which makes the data less valuable. The potential field data can be enhanced by using a downward continuation technique, such that the potential field at lower elevation or even underground within the harmonic source-free region (Pašteka et al., 2012) can be effectively estimated. In the wavenumber domain (Fourier spectral domain), the continuation can be carried out by multiplying a continuation factor with the spectrum of the observation data. Unfortunately, the downward continuation factor grows rapidly as the continuation distance increases. The high frequency components including noise in the observation data will be amplified and thus resulting a severe polluted solution. Therefore, using downward continuation in a wavenumber domain is an inherently unstable process. Using appropriate filters or constrains, stable downward continuation can be constructed. Dean (1958)) proposed a method to constrain the high frequency components. The use of a Wiener filter is investigated in Clarke (1969) and Pawlowski (1995). Recently, Pašteka et al. (2012) propose a robust wavenumber ⁎ Corresponding author. E-mail address: [email protected] (Y. Zhang).

http://dx.doi.org/10.1016/j.jappgeo.2016.01.009 0926-9851/© 2016 Elsevier B.V. All rights reserved.

domain method where the filter is designed based on the characteristics of Tikhonov regularization, Zeng et al. (2013) use an adaptive iterative Tikhonov method to apply a Tikhonov filter in each iteration in a wavenumber domain. The advantage of a wavenumber domain method is that the downward continuation process can be accelerated by fast Fourier transform (FFT). It has been proved that an appropriate designed filter can guarantee the accuracy and stability of downward continuation even for noisy data (Pašteka et al., 2012; Zeng et al., 2013). Another type of downward continuation method is based on the Taylor expansion, where the potential field at one elevation can be expanded by the potential field and its vertical derivative terms at another elevation. The success of the Taylor series method depends on the accuracy and stability in computing the vertical derivative terms. Fedi and Florio (2002, 2001) propose the ISVD method, where the odd vertical derivatives can be computed in a stable way, and the even order vertical derivative can be efficiently computed by finite difference. Zhang et al. (2013) propose a truncated Taylor series iterative scheme to achieve robust and stable downward continuation. Ma et al. (2013) compute the downward continuation by adding an upward continuation and a second vertical derivative at the observation plane, and the scheme can also be converted to an iterative version. The Taylor series method is capable of providing very accurate solution when the data are relatively clean, moreover the iterative Taylor series method usually has a fast convergence rate. It should be noted that both the wavenumber domain methods and the Taylor series methods can be accelerated by FFT. However, the FFT itself inherently introduces an artifact, and the FFT-induced artifact can be seen in many existing methods, this is particularly obvious in some iterative wavenumber domain methods (Zeng et al., 2013). The FFT-induced error in the downward continuation process has already

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

been studied by researchers in Trompat et al. (2003); Cooper (2004); and Pašteka et al. (2012). To resolve the difficulty, either extrapolation is needed to extend the original data (Pašteka et al., 2012), or a smaller window should be used to exclude the results near the boundary. For the Taylor series methods, besides the FFT-induced error, another problem is the robustness for the noisy data. Although the ISVD method (Fedi and Florio, 2002; Fedi and Florio, 2001) can be applied to compute the odd derivatives in a stable way, the even derivatives are still computed by the standard finite difference which is sensitive to the noise. Another iterative method such as that based on Taylor series has a similar problem (Zhang et al., 2013; Ma et al., 2013). Without a denoising procedure, it is hard to apply the Taylor series methods for the field data with relatively large noise. In summary, the FFT provides an efficient computation with the numerical complexity of order n logn, where n is the number of unknowns. To apply the FFT, a continuation process including the regularization is usually converted into the wavenumber domain. For this reason, regularized downward continuation in space domain has seldom been investigated. Zhang and Wong (2015) propose a numerical scheme for 3D gravity field inversion, where a special algebraic structure called Block-Toeplize Toeplize-Block (BTTB) matrix is utilized to make the scheme efficient. In this paper, we consider the conjugate gradient (CG) method utilizing the BTTB structure for downward continuation problem. The BTTB structure is derived from the downward continuation formulation in space domain, and it has the same numerical efficiency as the FFTbased methods. However, compared with the FFT-based methods, the proposed method induces very small artifact near the boundary, such that neither extrapolation nor tailoring process is required to reduce the boundary error. The characteristic of the BTTB structure allows the use of an iterative scheme without accumulating the error near the boundary. Combining the BTTB structure with the re-weighted regularized conjugate gradient method (BTTB–RRCG), a stable downward continuation method can be constructed. Here, all formulations are in time domain, such that various space domain regularization stabilizers can be applied. We compare the proposed computational scheme with other recently proposed schemes for downward continuation. The simulation results for synthetic and field data demonstrate that the proposed scheme is more accurate and robust for applications using clean and noisy data. In Section 2, the formulation of downward continuation is presented. We briefly introduce the Tikhonov regularized method (TR) (Pašteka et al., 2012), adaptive iterative Tikhonov method (AIT) (Zeng et al., 2013), and stable Taylor series methods (ITS) (Ma et al., 2013). Section 3 focuses on the proposed BTTB–RRCG scheme. In Section 4, synthetic field data are used to validate the proposed numerical scheme, and the result is compared with those obtained by the TR, AIT and ITS methods. A Gaussian noise from 0.1% to 5% of the maximum magnitude of the synthetic data is added to test the robustness. The error is analyzed by using RMS and the relative error in terms of L-2 norm and L-∞ norm. Particularly, the FFT-induced error near the boundary is investigated. In Section 5, we apply the proposed scheme to test cases with field data, and similar to the synthetic cases, the result is compared with other existing methods. 2. Mathematical background

h0  h ∞ ∞ ∫ ∫ h 2π ∞ ∞

T ðx0 ; y0 ; hÞ ðx 

lower elevation h, such that h0 N h. The downward continuation process is to seek T(x, y, h) at lower elevation from the potential field T(x, y, h0) at higher elevation. Denote the kernel as K, the integral Eq. (1) can be converted into the following convolution form ∞

∞

0

0

Tðx; y; h0 Þ ¼ ∫ ∞ ∫ ∞ Kðx  x0 ; y  y0 ; h0  hÞTðx0 ; y0 ; hÞdx dy ;

ð2Þ

which can be further simplified as Tðh0 Þ ¼ K  TðhÞ;

ð3Þ

where * denotes the convolution. According to the convolution theorem, F ðTðh0 ÞÞ ¼ F ðK  TðhÞÞ ¼ F ðKÞ  F ðTðhÞÞ;

ð4Þ

therefore, Tðh0 Þ ¼ F 1 ð F ðKÞ  F ðTðhÞÞÞ:

ð5Þ

Since ∞

∞

F ðKÞ ¼ ∫ ∞ ∫ ∞ Kðx; yÞe2πiðuxþvyÞ dxdy ¼ eðh0 hÞ

pffiffiffiffiffiffiffiffiffiffi u2 þv2

;

ð6Þ

denote T(h0) and T(h) by Th0 and Th, respectively, then Eq. (3) can be rewritten into the following matrix form T h0 ¼ F 1 Λ FT h ;

ð7Þ

where F and F-1 are the Fourier matrices corresponding to a 2D Fourier transform, and Λ is the continuation kernel K in a wavenumber domain pffiffiffiffiffiffiffiffiffiffi 2 2 given by Eq. (6). Consider h0 - h N 0, the kernel eðh0 hÞ u þv is stable, since the high frequency component can be compressed. This explains why the upward continuation is a stable process. According to Eq. (7), the most straightforward way to conduct a downward continuation is Th ¼ F1 Λ1 FTh0 ;

ð8Þ

pffiffiffiffiffiffiffiffiffiffi 2 2 where Λ-1 is given by eðh0 hÞ u þv . Obviously, since h0 - hN 0, the kernel pffiffiffiffiffiffiffiffiffiffi 2 2 given by eðh0 hÞ u þv will amplify all frequency components in Th0, such that the solution of Th will be polluted by the high frequency component or noise in Th0. Denote Th by T, Eq. (8) can be rewritten in a simplified form as ̂ ; T̂ ¼ Λ‐1 Th0

ð9Þ

where T̂ and Tĥ 0 are the potential field in wavenumber domain with heights h and h0. Therefore, a downward continuation is an inherently unstable process, and conventionally, there are mainly two approaches to resolve this issue: Tikhonov regularization in wavenumber domain and the Taylor series methods. 2.1. Wavenumber domain Tikhonov regularization method

The relationship between the potential field data at two observation planes is given by Blakely (1996): Tðx; y; h0 Þ ¼

75

x0 Þ2

2

þ ðy  y0 Þ2 þ ðh  h0 Þ

i3=2 ;

ð1Þ

where x and y are the horizontal coordinates, T(x, y, h0) is the observation field at higher elevation h0, and T(x, y, h) is the unknown field at

Let us denote the downward continuation formulation (Eq. (1)) into the following form: Th0 ¼ AT;

ð10Þ

where A is the upward continuation operator. As we have shown before, solving Eq. (10) is an ill-posed problem, which is equivalent to compute Eq. (8). Tikhonov and Arsenin (1977) propose an effective way to

76

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

resolve this problem. Instead of solving Eq. (10) directly, they convert the problem into the following minimization problem

(2013) introduced a method to remove the odd derivative from the Taylor formulation in the following form

n   o min jjWd AT  Th0 jj2 þ μjjWm T  Tref jj2 ;

Tðx; y; hÞ ≈ 2Tðx; y; h0 Þ  Tðx; y; h0 þ ΔhÞ þ

ð11Þ

where μ is the regularization parameter, Wd and Wm are the data weighting matrix and model weighting matrix, respectively, and Tref is the prior information. In a downward continuation, Tref is usually a zero vector. Let Wd and Wm be the identity matrix, then the solution of Eq. (11) can be given by  1 AT Th0 : TTik ¼ AT A þ μI

ð12Þ

Convert Eq. (12) into the wavenumber domain, we have (Zeng et al., 2013): ̂ ¼ TTik

Λ2 2

Λ þμ

̂ Λ−1 T̂h0 : Λ−1 T̂h0 ¼ LTik

ð13Þ

Compare Eqs. (13) and (9), it can be seen that the R̂Tik is the Tikhonov regularization filter in a wavenumber domain, which is given by L̂Tik

pffiffiffiffiffiffiffiffiffiffi 2 2 e−2h u þv pffiffiffiffiffiffiffiffiffiffi ¼ : 2 2 e−2h u þv þ μ

2

∂ Tðx; y; h0 Þ 2 h0 ; ∂z2

ð18Þ

where T(x, y, h0 + Δh) is the upward continuation of the observation field. Since the upward continuation is a stable process, the downward continuation scheme ((18)) should be stable. Both the ISVD method and Taylor formulation ((18)) require computation for the second vertical derivative by a finite difference. From the Laplace equation, the second vertical derivative can be computed by the second order horizontal derivative in x and y direction. Although finite difference is an efficient way to compute the horizontal derivative (Fedi and Florio, 2002), it is sensitive to noisy data. In our simulation, assume the grid size is given by h, to reduce the error to the minimum, we apply a second order central scheme ((19)) and a forward scheme ((20)) to compute the second horizontal derivative in the interior and on the boundary. Hence the overall accuracy for the second derivative is of second order. 00

f ðxÞ ¼ 00

f ðxÞ ¼

f ðx þ hÞ  2 f ðxÞ þ f ðx  hÞ h

ð19Þ

2

f ðx þ 2hÞ  2 f ðx þ hÞ þ f ðxÞ

ð20Þ

2

h

ð14Þ

3. BTTB–RRCG iterative scheme Based on Eq. (14), Zeng et al. (2013) propose an iterative Tikhonov scheme in the following form: pffiffiffiffiffiffiffiffiffiffi −2h u2 þv2 e ̂ ̂ pffiffiffiffiffiffiffiffiffiffi þ ; Λ−1 Rn−1 T̂n ¼ Tn−1 2 2 e−2h u þv þ μ

ð15Þ

̂ Λ1 T̂h , Rn̂ ¼ Tĥ  Λ1 T̂n . with the initial value T0 given by T̂0 ¼ LTik 0 0 Actually, the Tikhonov regularization ((11)) may have different forms with different Tikhonov regularization stabilizers. Applying the Tikhonov formulation given in Tikhonov et al. (1969); Pašteka et al. (2012) propose another efficient Tikhonov regularization filter in a wavenumber domain in the form of ̂ ¼ LTik

1

1þ

μ ðu2

pffiffiffiffiffiffiffiffiffiffi : 2 2 þ v2 Þeh u þv

ð16Þ

Our numerical simulations show that as a one-step Tikhonov regularization filter, Eq. (16) is better than Eq. (14) in terms of robustness and accuracy, while the iterative version of Eq. (14) is slightly better than Eq. (16).

We now present a new approach for downward continuation. Instead of converting the downward continuation formulation into a wavenumber domain as in Eq. (6), we work with a space domain formulation and discretize Eq. (1) as T ðxðiÞ; yðiÞ; h0 Þ ¼

N X M X GðxðiÞ; yðiÞ; h; x0 ð jÞ; y0 ðkÞ; hÞT ðx0 ð jÞ; y0 ðkÞ; hÞΔxΔy;

ð21Þ

j¼1 k¼1

i ¼ 1; 2; ⋯; N  M;

where h0 and h are defined as before, N and M are the number of data grids in the x and y direction, Δ x and Δy are the grid interval in the x and y direction. Denote the data points on the lower plane indicated by Th(x V(i), y V(j), h), where i = 1 , ⋯ , N , j = 1 , ⋯ , M. Renumbering the data grids to Th(x V(l), y V(l), h), where l = 1 , ⋯ , N * M, which means that we rearrange the index of the data without changing the total number of data points. Thus, Eq. (21) can be rewritten as T ðxðiÞ; yðiÞ; h0 Þ ¼

NM X

GðxðiÞ; yðiÞ; z; x0 ðlÞ; y0 ðlÞ; h0 ÞT ðx0 ðlÞ; y0 ðlÞ; hÞΔxΔy;

l¼1

i ¼ 1; 2; ⋯; N  M; ð22Þ

2.2. Taylor series methods Differed from the FFT-based iterative method, the Taylor series method is used to express the potential field at the elevation h by the potential field at another elevation h0 as

where. Gði; l; hÞ ¼ h

2

∂Tðx; y; h0 Þ 1 ∂ Tðx; y; h0 Þ 2 Δh þ ⋯ Δh þ 2! ∂z ∂z2 m 1 ∂ Tðx; y; h0 Þ m Δh ; þ m! ∂zm

Tðx; y; hÞ ¼ Tðx; y; h0 Þ þ

ð17Þ

where Δh =h0 - h. The ISVD method proposed by Fedi and Florio (2002, 2001) uses vertical integrating the field in a wavenumber domain to compute the odd vertical derivative, which has a good the stability for noisy data. Ma et al.

h  h0 2

2

2

ðxðiÞ  xðlÞÞ þ ðyðiÞ  yðlÞÞ þ ðh  h0 Þ

i32 :

ð23Þ

By Eqs. (22) and (23), the original downward continuation formulation (Eq. (1)) can be approximated by the linear system T0 ¼ GT:

ð24Þ

Here, T0 is the discretized observation field, T is the unknown field data, and G is a (N × M) by (N × M) BTTB matrix generated

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

from the discretization ((22)). The BTTB matrix is given in the following form: 2

GMN

Gð0Þ 6 Gð1Þ 6 ¼6 6⋮ 4 GðM−2Þ GðM−1Þ

Gð−1Þ Gð0Þ Gð1Þ ⋯ GðM−2Þ

⋯ Gð−1Þ Gð0Þ ⋱ ⋯

Gð2−MÞ ⋯ ⋱ ⋱ Gð1Þ

3 Gð1−MÞ Gð2−MÞ 7 7 7; ⋮ 7 Gð−1Þ 5 Gð0Þ

ð25Þ

ðmÞ g 6 0ðmÞ 6 g1 6

ðmÞ g 1 ðmÞ g0 ðmÞ g1

GðmÞ ¼ 6 6⋮ 6 ðmÞ 4 g N2

⋯ ðmÞ g N2

ðmÞ

g N1

⋯ ðmÞ g 1 ðmÞ g0 ⋱ ⋯

ðmÞ g 2N

⋯ ⋱ ⋱ ðmÞ g1

ðmÞ g 1N ðmÞ g 2N

⋮ ðmÞ g 1 ðmÞ g0

ð26Þ

where gi is constant along its diagonals and the value is defined by Eq. (23) as Gðm; n; hÞ ¼ h

h  h0 2

m2 Δx2 þ n2 Δy2 þ ðh  h0 Þ

i32 ;

ð27Þ

where m = 0 , 1 , ⋯ , M - 1, n = 0 , 1 , ⋯ , N - 1. One important property of a BTTB matrix is that the first row and the first column contain all information of a given matrix. Consequently, for the matrix G given by Eq. (25), we only need to store the first row and first column, such that the storage requirement for the sensitivity matrix can be dramatically reduced from (N * M)2 to 2(N * M). It should be noted that the BTTB matrix is always a symmetric matrix regardless of the choice for Δx and Δy. Hence, the storage requirement for the sensitivity matrix is further reduced to N * M. Another important feature of the BTTB structure is that any BTTB matrix can be embedded into a Block-Circulant Circulant-Block (BCCB) matrix as (Chan and Jin, 2007):  Cn ¼

GMN 

  ; GMN

ð28Þ

such that the matrix–vector product for any BTTB matrix can be computed via 

GMN 

    GMN T  T ; ¼ GMN 0 †

ð29Þ

where × is the BTTB matrix determined by GMN, T is a vector as defined in Eq. (24), 0 is a zero vector with the same dimension as T, and † is the part to be dropped. The details of this operation can be found in Chan and Jin (2007). The embedding of the BTTB matrix in Eq. (28) is very critical. It is known (Davis, 1979) that the BCCB matrix C can be diagonalized by the Fourier matrix F and its conjugate transpose, i.e., C ¼ F  Λ̂ F;

ð30Þ

where F is the Fourier matrix. Recall that the F in Eq. (8) is the Fourier transform operator. Applying the Fourier transform operator to a given matrix, is equivalent to performing a premultiplication between the Fourier matrix and the given matrix (Chan and Jin, 2007). The Fourier matrix is given as: pffiffiffiffiffiffiffi 1 2πijk ð F n Þ j;k ¼ pffiffiffi e n ; i ¼ 1 n for 0 ≤j , k ≤ n- 1.

Eq. (30) is the eigenvalues of C, and it has a different dimension from Λ. The computation of Λ̂ in Eq. (30) is more efficient than computing the Λ in Eq. (6). For a n by n circulant matrix C, the Λ̂ can be computed via FFT by observing that the first column of Fn is p1ffiffin 1n , where ð1; 0; …; 0ÞT ∈ Rn , then by using Eq. (30), we have

3 7 7 7 7; m ¼ 0; 1; ⋯; M  1; 7 7 5

It should be noted that the Λ̂ in Eq. (30) is totally different from the Λ in Eq. (7). Recall in Eq. (7), the Λ is the downward continuation kernel in wavenumber domain, and it can be computed by Eq. (6). However, Λ̂ in

1n ¼ ð1; 1; …; 1ÞT ∈ Rn is the vector consisting of all ones. Let e1 ¼

in which each block G (m) is a Toeplitz matrix given by 2

77

ð31Þ

1 F n C n e1 ¼ pffiffiffi Λn 1n ; n

ð32Þ

which implies that Λ̂ can be computed by applying fast Fourier transform to C. The BTTB matrix has many other attractive properties including the construction for an efficient preconditioner, and recent work on preconditioners has been reported in Chan and Jin (2007). Now the downward continuation problem ((1)) is converted into solving the linear system ((24)). Similar to a wavenumber domain method, the Tikhonov regularization ((12)) can also be applied. In this study, we consider using the conjugate gradient type method to solve the regularization problem ((12)), since theoretically it has a rapid rate of convergence (Allen and Isaacson, 2011). Compared with the steepest decent method (Xu et al., 2007; Zeng et al., 2013), which can be regarded as the first order gradient method, the conjugate gradient method is a second order gradient method. Combining the BTTB property (Eqs. (30) and (29)) with the RRCG scheme in (Zhdanov, 2002), we develop a BTTB–RRCG scheme. The details of the standard RRCG scheme can be found in Zhdanov (2002): RRCG. For the minimization problem given in Eq. (11), denote the observation field Th0 by d, and the potential field at the target plane Th by m. Let m0 be an initial approximation, α0 be the initial regularization parameter. The re-weighted regularized conjugate gradient (RRCG) algorithm is given as follows. r0 ¼ Am0 −d; s0 ¼ Wm ðm0 −mref Þ ; Iα0 0 ¼ Iα0 ðm0 Þ ¼ AT W2d r0 þ α 0 Wm s0 ; for n ¼ 1; 2; 3⋯ rn ¼ Amn −d; sn ¼ Wm ðmn −mref Þ; Iαn n ¼ Iαn ðmn Þ ¼ AT W2d rn þ α n Wm sn ;    2 n−1  βαn n ¼ Iαn n 2 =Iαn−1 ; α

n−1 ~I αn ¼ Iαn n ~In−1 ; Iα0 0 ¼ Iα0 0 ; n

  h   i ~α n ¼ ~Iα n T Iαn = ~Iα n T AT W2 A þ αW 2 ~Iα n ; k n n n n d m n α

~α n ~I n ; γ ¼ ksnþ1 k2 =ksn k2 ; mnþ1 ¼ mn −k n n  α nþ1 ¼

αn α n =γ

if γ ≤1 i f γN1:

By combing the discretization process and the properties of the BTTB matrix with the RRCG algorithm, the BTTB–RRCG algorithm is given as follows: BTTB–RRCG. For the downward continuation problem ((10)), denote the observation field Th0 by d, and the potential field at the target plane Th by m.

78

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

Step 1. Discretize the upward continuation operator ((1)) by using Eq. (22). Step 2. Generate the matrix G in Eq. (24) in the form of Eq. (25), where each element in the forward continuation matrix is defined by Eq. (27), and a compact storage format is used to store G. Step 3. Choose appropriate regularization in time domain, and convert the problem ((24)) into the minimization problem ((11)). The regularization stabilizer is also in the form of the BTTB matrix. Step 4. Apply the RRCG algorithm mentioned above to the minimization problem. In each operation, the system matrix G is embedded into a BCCB matrix as in Eq. (29), such that the matrix–vector product can be conducted by using Eqs. (30) and (32). The details of the compact storage in step 2 can be found in Chan and Jin (2007), and the regularization stabilizer in terms of the BTTB structure has been reported in Zhang and Wong (2015). The initial regularization parameter α0 is chosen according to a trial and error method. It should be noted that when the field data is clean without noise, the BTTB–CG method should have a better performance than the BTTB–RRCG method, since the regularization itself will inevitably introduce a certain degree of distortion in the downward continuation solution. On the other hand, it has been shown that by applying a regularization, the stability of the gradient type methods can be greatly improved (Zhdanov, 2002; Zeng et al., 2013). In the next section, we will verify the proposed scheme by using synthetic and field data, and investigate the sensitivity of regularization parameter α in the proposed BTTB–RRCG method. 4. Simulation using synthetic data To validate the proposed BTTB–RRCG scheme, we now consider a test case with a synthetic magnetic field data. The computation is carried out by a laptop with i7-3630 CPU and 12 G RAM. The numerical results will be compared with those obtained by the Tikhonov regularized method (TR) (Pašteka et al., 2012), the adaptive iterative Tikhonov method (AIT) (Zeng et al., 2013) and the iterative Taylor series method (ITS). We consider two synthetic tests in this section. By tradition, the synthetic field is generated by a 3D synthetic susceptibility or density model. However, in the first synthetic test, we present another approach for the synthetic test. First, we design a pseudo-magnetic field on the ground level (h = 0 m), and we refer this to a pseudomagnetic field since it is not generated by a 3D susceptibility model. Then, we realize an upward continuation of the pseudo-magnetic field on the ground level to two different elevations h1 and h2 generating two synthetic fields Th1 and Th2, where h2 N h1. Finally, we conduct a downward continuation for Th2 with a continuation distance Δh = h2 - h1, such that an inferred potential field Th1 can be estimated. By comparing Th1 with Th1, and evaluating the statistics of the difference between Th1 and Th1, we study the characteristics of the downward continuation. Performing the synthetic downward continuation simulation in such way is reasonable, because we do not really need to compute the underground 3D susceptibility distribution. Moreover, it has been shown that with some assumptions, the 3D susceptibility distribution can be simplified into a 2D case (Zhang et al., 2015). Therefore, a pseudo-magnetic field distribution is sufficient. More importantly, we can now design the distribution of the magnetic field and include anomalies with various frequency components and anomalies near the boundary, such that the edge-effect can be investigated. However, in a tradition approach, the 3D susceptibility model is always in the interior of the model domain, and the generated field has a very small or zero value near the boundary. Consequently, there is almost no edge-effect in the computation. The synthetic field data is too simple to evaluate the performance of a downward continuation scheme for test cases with real field data, in which the fields are usually complicated near the boundary. The design of a pseudo-magnetic field on the ground

Fig. 1. Pseudo-magnetic field on the ground level (h=0 m).

level is also much more convenient than constructing a complicated 3D susceptibility distribution. The pseudo-magnetic field on the ground level is shown in Fig. 1, where Δx = Δy = 9.98 m. The synthetic field contains two semi-circle anomalies on the upper left and lower right boundaries, and there is no other anomaly near the boundary. In the interior of the synthetic field, there are two circle anomalies with sharp boundary variation and a swirl shape anomaly with two shape corners. The synthetic field is designed to include the field variation with different frequencies, and the two semi-circle anomalies on the boundary are used to investigate the edge-effect of the scheme. We now conducted an upward continuation to the synthetic field as shown in Fig. 1 with continuation distance Δh =50 m and Δh= 250 m, where the T50 and T250 are shown in Fig. 2(a) and (b), respectively. Then we carry out a downward continuation for the field at h = 250 m in Fig. 2 with a continuation distance Δh = 200 m, and the performance of the computational scheme can be evaluated by comparing the downward continuation field T50 and T50. Note that the field data is slightly perturbed by adding 0.005% Gaussian noise to the maximum of the magnitude of T250. For the regularization parameters, an optimal regularization parameter can be determined using the C - norm method (Pašteka et al., 2012). However, to compare with different regularization methods, we apply a straightforward procedure to evaluate the optimal regularization parameter μ in the TR, AIT and BTTB–CG methods by plotting the RMS of the downward continuation error with different μ, where the RMS is defined in Eq. (33). Since the exact downward continuation result is known, the value of μ is obviously optimal. The RMS vs μ curves for the TR, AIT and BTTB–RRCG methods are shown in Fig. 3. In this study, optimal parameters used in the simulation are μTR = 514, μAIT = 0.0157, and μRRCG =0.095. We have considered data with different noise levels in the range from 0% to 5%. However, the optimal value for μ is almost the same, since the optimal regularization parameter is not sensitive when the noise level is less than or equal to 5%, and this has been confirmed for the TR, AIT and BTTB–RRCG methods. Using the TR method as an example and introducing 0%, 2% and 4% noise to the data, we plot the RMS vs μ graph as shown in Fig. 4. It can be seen that the optimal regularization parameters for the data with different noise levels are between 500 and 600, which are very similar. The optimal iteration of AIT is 2 from a trial and error method. By plotting the residue given by e = || ATh - T0 ||2 vs iteration number N, we can estimate the iteration numbers in the BTTB–RRCG. For the ITS method, once the noise is added, the iteration will amplify the noise if no denoising filter is applied. Therefore, we only consider a non-

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

79

Fig. 2. The upward continuation of the pseudo-magnetic field in Fig. 1 to the elevation of (a) h=50 m and (b) h=250 m.

iterative version for ITS, and we rename the method as TS in the following. The downward continuation results for the field in Fig. 2(b) with 0.005% Gaussian noise by the TR, AIT, TS and BTTB–RRCG methods are illustrated in Fig. 5. The results shown in Fig. 5 confirm that all four methods produce a stable solution. The error distribution is computed by plotting the difference between the continuation results at h = 50 m in Fig. 5 and the accurate field at h = 50 m in Fig. 2(a). The error distribution by the four methods is shown in Fig. 6. To quantify the performance of these methods, we define the RMS, RE2 and RE∞ as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N X M u 1 X ðTcon −Treal Þ2 ; RMS ¼ t N  M i¼1 j¼1

ð33Þ

RE2 ¼

kTcon −Treal k2 ; kTreal k2

ð34Þ

RE∞ ¼

jjTcon  Treal jj∞ : jjTreal jj∞

ð35Þ

Table 1 reports the performance of TR, AIT, TS and BTTB–RRCG in terms of RMS, RE2, RE∞ and the computing time. From Fig. 6, it can be seen that the edge-effect is clearly evident on the upper and lower boundaries by the TR and AIT methods. In contrast, the edge effect is much smaller for the TS and the proposed BTTB–RRCG method. It is interesting to note that the RMS and RE 2 for the TS are relatively low, however, the RE∞ is quite high. From the error distribution illustrated in Fig. 6(c), the effect due to noise is noticeable. From Fig. 7(c) and Table 1, the BTTB–RRCG method

Fig. 3. RMS vs μ graph for (a) TR method; (b) AIT method; and (c) BTTB–RRCG method.

80

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

Fig. 4. RMS vs μ graph for TR method under the different levels of noise.

produces good results in terms of RMS, RE2 and RE∞. The BTTB–RRCG method requires more computing time than other methods, since it is an iterative scheme. However, the computing time per iteration is similar than the TR and AIT methods. Note that the numerical complexity of BTTB–RRCG is of order n log n, which is same as a FFT base method. It can be easily verified that the error cannot be reduced by increasing the field data resolution. By increasing the resolution, the

frequency in the field data will also increase leading to amplification of the noise. Moreover, the edge-effect is unavoidable by using the FFT method. For the Taylor series method without the use of a regularization, the TS is very sensitive to noise. Adding a 1% Gaussian noise to the field data, the RMS, RE2 and RE∞ of TS increase rapidly to 0.2470, 695% and 4200%, respectively. By incorporating a denoising preprocessing, the robustness of TS can be improved. However, the

Fig. 5. The downward continuation with Δh=200 m for the 0.005% noised magnetic field in Fig. 2(b) by using (a) TR; (b) AIT; (c) TS and (d) BTTB–RRCG methods.

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

81

Fig. 6. The downward continuation error distribution by (a) TR; (b) AIT; (c) TS and (d) BTTB–RRCG methods.

Table 1 Computational errors using TR, AIT, TS and BTTB–RRCG methods.

RMS RE2 RE∞ Computing time (s)

TR

AIT

TS

BTTB–RRCG

0.0234 6.25% 17.00% 0.065 s

0.0233 6.17% 16.03% 0.068 s

0.0202 4.63% 21.76% 0.092 s

0.0187 3.99% 5.74% 6.738 s

denoising itself usually induces new errors. Therefore, in the following, we only compare the proposed method with regularized TR and AIT methods.

To further investigate the robustness, we increase the Gaussian noise from 1% to 5% in the synthetic field data. The errors of downward continuation solutions by various methods are shown in Table 2 and Fig. 7. From Fig. 7 and Table 2, we conclude that the TR, AIT and BTTB– RRCG methods have a robust performance for noisy data, while the BTTB–RRCG method provides the most accurate solution in terms of RMS and RE2. Moreover, the computed solutions by BTTB–RRCG have the smallest RE∞ confirming that the edge-effect is small compared with the other two regularized methods. Different from the conventional downward continuation methods, the proposed BTTB–RRCG method is an iterative method in space domain, which is seldom investigated due to the computation workload.

Fig. 7. The downward continuation error distribution for 5% noised data by (a) TR; (b) AIT; and (c) BTTB–RRCG methods.

82

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

Table 2 Error of downward continuation by using TR, AIT and ITM methods with different levels of noise. % of noise

% % % % % %

RMS AIT

BTTB–RRCG

TR

AIT

BTTB–RRCG

TR

AIT

BTTB–RRCG

0.0234 0.0235 0.0235 0.0236 0.0239 0.0244

0.0233 0.0233 0.0233 0.0233 0.0235 0.0233

0.0187 0.0190 0.0188 0.0193 0.0198 0.0208

6.25% 6.29% 6.22% 6.34% 6.53% 6.80%

6.17% 6.17% 6.20% 6.17% 6.28% 6.20%

3.99% 4.12% 4.02% 4.26% 4.46% 4.95%

16.99% 17.04% 17.04% 17.48% 18.05% 18.95%

16.03% 16.00% 16.00% 15.87% 16.32% 14.91%

5.74% 6.38% 7.16% 7.86% 9.61% 10.71%

Table 3 Storage and computing time by conventional and BTTB–RRCG methods. Problem resolution

*128 *256 *512 *1024

RE∞

RE2

TR

Conventional CG method

BTTB–RRCG

Storage requirement

Computing time

Storage requirement

Computing time

1.47 GB 23.52 GB 376.32 GB 6000 GB

260 s 1.16 h 18.6 h 297 h

0.31 MB 1.22 MB 4.95 MB 30 MB

0.33 s 1.63 s 6.11 s 24.88 s

However, taking advantage of the BTTB structure makes the CG type methods as effective as wavenumber domain methods. Table 3 reports the computing time and storage requirement for the various data sizes using the conventional method and BTTB–RRCG method. The resolution of M* N implies that the field data is given by an M by N matrix. From Table 3, it is noted that as the size of the problem increases, both the computing cost and storage requirement for the conventional CG increase exponentially. In contrast, the computational complexity for the BTTB–RRCG increases linearly. The second synthetic test focuses on a density model as shown in Fig. 8, and the gravity field is generated by the synthetic density model. The density model consists of two dipping prisms underground, where the densities of the long prism and the short prism are 1.0 g/cm 3 and 0.8 g/cm 3. The depth from the ground to the top of the density anomaly is 100 m. The density anomalies are used to generate a gravity field at h = - 50 m and h = 200 m as illustrated in Fig. 9a and b respectively, where the grid interval is Δx = Δy = 20 m. Denote the gravity field at h= - 50 m and h =200 m by T-50 and T200, we add 2.5% Gaussian noise to T200 and then apply the TR, AIT and the proposed BTTB–RRCG to conduct the downward continuation to T200 with continuation distance

h = 250 m to compute the field T50 . The downward continuation is conducted to the underground, because within the harmonic sourcefree region, the downward continuation should be always feasible. The continuation results are shown in Fig. 10, and the error distribution is given in Fig. 11. By comparing T50 in Fig. 10 with T-50 in Fig. 9a, we also report the error in terms of RMS, RE∞ and RE2 in Table 4. From the downward continuation results in Table 4, the proposed BTTB–RRCG method is clearly more accurate than the TR and AIT methods in terms of RMS, RE∞ and RE2. More importantly, consider that for the density model in Fig. 10, all anomalies are positive which means generated gravity fields should be positive. However, in Fig. 10, both the TR and AIT methods induce negative values on the left side. However, the proposed BTTB–RRCG scheme perfectly preserves the positivity property, and has a much smaller boundary effect than other two methods. 5. Applications using real field data Now, we apply the proposed BTTB–RRCG scheme using real field data. The field data shown in Fig. 12a is the magnetic field distribution at the ground level, where Δx = Δy = 10 m. The upward continuation of the field data to h= 200 m is shown in Fig. 12(b), By adding 2.5% Gaussian noise to the potential field in Fig. 12(b), the TR, AIT and BTTB–RRCG methods are used for a downward computation with a continuation distance h = 200 m. The same regularization parameter for the synthetic data is used for the real field data applications. In real applications, the exact downward continuation results are always unknown, therefore we can only have an estimation to the value for the optimal regularization parameters. In our study, we use the optimal regularization parameters obtained in the synthetic case for the real field data. Recall that in the downward continuation formulation ((10)), the upward continuation operator A depends only on the

Fig. 8. Synthetic density model.

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

83

Fig. 9. Generated field by density model in Fig. 8.

Fig. 10. The downward continuation results by (a) TR; (b) AIT; and (c) BTTB–RRCG methods.

Fig. 11. The downward continuation error distribution for by (a) TR; (b) AIT; and (c) BTTB–RRCG methods.

Δx ,Δy, and h. Once these parameters are fixed, the spectral characteristic of the continuation operator is determined. In the first synthetic case, Δx = Δy = 9.98 m, and h = 200 m, while in the real field data, dx = Table 4 Error of downward continuation by using TR, AIT and ITM methods for synthetic gravity data.

RMS RE∞ RE2

TR

AIT

BTTB–RRCG

0.0648 24.17% 31.18%

0.0607 23.52% 31.58%

0.0430 10.67% 19.74%

dy = 10 m, and h = 200 m, which means that the spectral characteristics of the continuation operators are similar between the first synthetic case and real field case. The downward continuation results for the real field data are shown in Fig. 13, and their error distributions are illustrated in Fig. 14. From Figs. 14(a) and 16(b), we observe that the edge-effect is evident near the boundary for the TR and AIT methods. A simple procedure to improve the computed solution is to remove a layer near the boundary. Figs. 15 and 16 illustrate the solutions by tailoring 40 grids (i.e., 400 m) from the edge for the solutions shown in Figs. 13 and 14. It is important to note that the solution computed by the BTTB–RRCG

84

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

Fig. 12. (a) The real magnetic field data at h=0 m; (b) the upward continuation by Δ=200 m of the field data in Fig. 12a.

Fig. 13. The downward continuation by Δh=200 m for 2.5% noised real field data by (a) TR; (b) AIT; and (c) BTTB–RRCG methods.

Fig. 14. The downward continuation error distribution for 2.5% noised real field data by (a) TR; (b) AIT; and (c) BTTB–RRCG methods.

Fig. 15. The downward continuation by Δh=200 m after tailoring for 2.5% noised real field data by (a) TR; (b) AIT; and (c) BTTB–RRCG methods.

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

85

Fig. 16. The downward continuation error distribution after tailoring for 2.5% noised real field data by (a) TR; (b) AIT; and (c) BTTB–RRCG methods.

method is very stable and with little edge effect. The BTTB–RRCG scheme is an accurate scheme even without the use of any extrapolating and tailoring process. From the tailored computed solutions for downward continuation shown in Figs. 15 and 16, it is clear that the BTTB–RRCG method is robust and accurate for the real field data compared with the TR and AIT methods. By evaluating downward continuation results by the TR, AIT and BTTB–RRCG methods in terms of RMS, Table 5 reports the values of the L2 norm and L∞ norm. Fig. 17 displays the convergence rate of the BTTB–RRCG for applications using the field data application. It is noted that the residue becomes stable after 15 iterations, hence we can terminate the iteration much earlier than the prescribed 35 iterations. The robustness and accuracy of the proposed BTTB–RRCG method are verified by the simulation results based on test cases using synthetic and field data. The CG type method is a popular iterative technique for solving large scale linear equations, and it is particularly efficient for large sparse matrices. The main computational work per iteration is typically dependent on the matrix–vector product operation. In the present applications, the matrix A is a full matrix, and without taking advantage of the BTTB structure, it requires considerable computing resources. We compare the computational complexity using the conventional RRCG method and the BTTB–RRCG method. The results speak for itself. It is

Table 5 Error by TR, AIT and BTTB–RRCG methods for real field data.

RMS RE2 RE∞

clear that the utilizing the BTTB-structure dramatically reduces both the storage and computing cost, and a large scale downward continuation problem can be computed with a modest computing resource. 6. Conclusion In this paper, we propose an efficient approach by utilizing the BTTB structure with the re-weighted regularized conjugate gradient method. Not only does the BTTB structure greatly enhance the stability and robustness of the downward continuation computation, but the computation workload and storage requirement are also significantly reduced. We demonstrate that compared with the conventional wavenumber domain FFT-based methods, the BTTB–RRCG scheme induces negligible artifact near the boundary. The simulation results using the clean and highly perturbed synthetic field data as well as the real field data verify that the proposed method is more accurate and robust compared to the Taylor-series and wavenumber domain regularized methods. It is also easier to implement various regularization stabilizers directly in time domain instead of wavenumber domain. By plotting the residue at each iteration, we can quickly estimate the number of iterations required for the BTTB–RRCG method. The BTTB structure allows the design of an optimal preconditioner to further accelerate the convergence rate, in which the preconditioner is also in the form of a BTTB matrix. Acknowledgment

TR

AIT

BTTB–RRCG

0.0348 3.39% 6.02%

0.0357 3.55% 6.52%

0.0240 1.61% 3.37%

The research is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC 05664). References

Fig. 17. Convergence rate of BTTB–RRCG method in the real field data application, where e=||ATh -T0 ||2.

Allen, M.B., Isaacson, E.L., 2011. Numerical Analysis for Applied Science vol. 35. John Wiley & Sons. Blakely, R.J., 1996. Potential Theory in Gravity and Magnetic Applications. Cambridge University Press. Chan, R.H., Jin, X.Q., 2007. An Introduction to Iterative Toeplitz Solvers vol. 5. SIAM. Clarke, G.K., 1969. Optimum second-derivative and downward-continuation filters. Geophysics 34 (3), 424–437. http://dx.doi.org/10.1190/1.1440020. Cooper, G., 2004. The stable downward continuation of potential field data. Explor. Geophys. 35 (4), 260–265. http://dx.doi.org/10.1071/EG04260. Davis, P.J., 1979. Circulant Matrices. American Mathematical Soc. Dean, W.C., 1958. Frequency analysis for gravity and magnetic interpretation. Geophysics 23 (1), 97–127. http://dx.doi.org/10.1190/1.1438457. Fedi, M., Florio, G., 2001. Detection of potential fields source boundaries by enhanced horizontal derivative method. Geophys. Prospect. 49 (1), 40–58. http://dx.doi.org/10. 1046/j.1365-2478.2001.00235.x. Fedi, M., Florio, G., 2002. A stable downward continuation by using the ISVD method. Geophys. J. Int. 151 (1), 146–156. http://dx.doi.org/10.1046/j.1365-246X.2002. 01767.x. Ma, G., Liu, C., Huang, D., Li, L., 2013. A stable iterative downward continuation of potential field data. J. Appl. Geophys. 98, 205–211. http://dx.doi.org/10.1016/j.jappgeo. 2013.08.018.

86

Y. Zhang et al. / Journal of Applied Geophysics 126 (2016) 74–86

Pašteka, R., Karcol, R., Kušnirák, D., Mojzeš, A., 2012. Regcont: a Matlab based program for stable downward continuation of geophysical potential fields using Tikhonov regularization. Comput. Geosci. 49, 278–289. http://dx.doi.org/10.1016/j.cageo.2012.06. 010. Pawlowski, R.S., 1995. Preferential continuation for potential-field anomaly enhancement. Geophysics 60 (2), 390–398. http://dx.doi.org/10.1190/1.1443775. Tikhonov, A.N., Arsenin, V.Y., 1977. Solutions of Ill-posed Problems. Vh Winston. Tikhonov, A.N., Glasko, V.B., Litvinen, O.K., Melikhov, V.P., 1969. Analytic continuation of a potential in the direction of disturbing masses by the regularization method. Geophysics 34 (6), 1020. Trompat, H., Boschetti, F., Hornby, P., 2003. Improved downward continuation of potential field data. Explor. Geophys. 34 (4), 249–256. http://dx.doi.org/10.1071/EG03249. Xu, S.Z., Yang, J., Yang, C., Xiao, P., Chen, S., Guo, Z., 2007. The iteration method for downward continuation of a potential field from a horizontal plane. Geophys. Prospect. 55 (6), 883–889. http://dx.doi.org/10.1111/j.1365-2478.2007.00634.x.

Zeng, X., Li, X., Su, J., Liu, D., Zou, H., 2013. An adaptive iterative method for downward continuation of potential-field data from a horizontal plane. Geophysics 78 (4), J43–J52. http://dx.doi.org/10.1190/geo2012-0404.1. Zhang, Y., Wong, Y.S., 2015. BTTB-based numerical schemes for three-dimensional gravity field inversion. Geophys. J. Int. 203 (1), 243–256. http://dx.doi.org/10.1093/gji/ ggv301. Zhang, H., Ravat, D., Hu, X., 2013. An improved and stable downward continuation of potential field data: the truncated Taylor series iterative downward continuation method. Geophysics 78 (5), J75–J86. http://dx.doi.org/10.1190/geo2012-0463.1. Zhang, Y., Wong, Y.S., Deng, J., Lei, S., Lambert, J., 2015. Numerical inversion schemes for magnetization using aeromagnetic data. Int. J. Numer. Anal. Model. 12 (4), 684–703. Zhdanov, M., 2002. Geophysical Inverse Theory and Regularization Problems: Methods in Geochemistry and Geophysics. Elsevier, Amsterdam.