Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 464793, 12 pages

http://dx.doi.org/10.1155/2015/464793

## Three-Dimensional Induced Polarization Parallel Inversion Using Nonlinear Conjugate Gradients Method

^{1}Key Laboratory of Geo-Detection, China University of Geosciences, Ministry of Education, Beijing 100083, China^{2}School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China^{3}Exploratory Drilling Corporation Well Logging Company, Daqing, Heilongjiang 163000, China

Received 19 March 2015; Accepted 27 April 2015

Academic Editor: Sergio Preidikman

Copyright © 2015 Huan Ma et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Four kinds of array of induced polarization (IP) methods (surface, borehole-surface, surface-borehole, and borehole-borehole) are widely used in resource exploration. However, due to the presence of large amounts of the sources, it will take much time to complete the inversion. In the paper, a new parallel algorithm is described which uses message passing interface (MPI) and graphics processing unit (GPU) to accelerate 3D inversion of these four methods. The forward finite differential equation is solved by ILU0 preconditioner and the conjugate gradient (CG) solver. The inverse problem is solved by nonlinear conjugate gradients (NLCG) iteration which is used to calculate one forward and two “pseudo-forward” modelings and update the direction, space, and model in turn. Because each source is independent in forward and “pseudo-forward” modelings, multiprocess modes are opened by calling MPI library. The iterative matrix solver within CULA is called in each process. Some tables and synthetic data examples illustrate that this parallel inversion algorithm is effective. Furthermore, we demonstrate that the joint inversion of surface and borehole data produces resistivity and chargeability results are superior to those obtained from inversions of individual surface data.

#### 1. Introduction

IP methods are important in geophysical electrical surveys. Surface exploration is used for detecting metallic and nonmetallic minerals, using various observational techniques. Borehole-surface, surface-borehole, and borehole-borehole exploration methods are also widely used in mining field contexts and give improved results for the prospecting for ores, oil, and gas.

Many previous resistivity and IP inversion methods have been developed, including 2D least-squares inversion [1], a 3D resistivity inversion using alpha center method [2], many linear and nonlinear IP inversion methods [3, 4], and a 3D IP conjugate inversion method [5].

Since 3D problems are very much larger than 1D and 2D problems, their solution is correspondingly more costly in computer time. Parallel MPI and GPU algorithms are generally used for numerical calculation. MPI is a library that is called by Fortran77, Fortran90, C, and C++ software and serves to communicate between processes. It has two basic patterns: the equal pattern and the principal-subordinate pattern [6].

GPUs are parallel processors with a large number of computation units and are superior to CPUs in both processing capability and memory bandwidth. They are cheaper and consume less power, and they are therefore used for large-scale, high-performance computation work [7, 8]. Programs are written in the familiar C, C++, Fortran, or an ever-expanding list of supported languages [9]. CULA is a set of GPU-accelerated linear algebra libraries utilizing the NVIDIA CUDA parallel computing platform to dramatically improve the speed of advanced mathematical computation. In other words, to use these functions an NVIDIA GPU with CUDA support is required. The library contains CULA Dense and CULA Sparse.

Solution is very slow for forward and “pseudo-forward” equations using the nonlinear conjugate gradients method. Therefore, a set of 3D IP forward and inversion parallel algorithms are needed for surface, borehole-surface, surface-borehole, and borehole-borehole techniques. This paper introduces a parallel algorithm developed by parallel programming to improve 3D NLCG IP inversion. First the theory of the NLCG algorithm for 3D IP inversion is presented. Second, implementation of the MPI and GPU parallel algorithm is briefly introduced. Finally, the efficiency of inversion codes is analyzed, with tables and synthetic data examples.

#### 2. Forward

The electrical potential for arbitrary conductivity in a half-space is given by the differential equation:where is the conductivity of rocks (S/m); denotes the potential (V) in absence of IP effect; is the electric current. The quantity on the right-hand side is in Ampere/meter^{3}, while currents are in Ampere; and is the Dirac delta function. The Neumann boundary condition is applied for the earth-air interface; and a mixed boundary condition is adopted at infinite boundaries. The equation can be written as [10]where is the angle between the radial distances from the source point to the outward normal spatial coordinate on the outer boundary. Then, the forward mapping operator is defined as

If ground is chargeable, then the potential , which is induced by the application of a constant current will be different from . According to Siegel’s formulation [11], with the conductivity replaced by , the effect of the chargeability of ground can be modeled by direct current resistivity forward mapping operator . Then, apparent chargeability can be defined by equationwhere and are the apparent resistivity computed by the potential and , respectively. So the apparent chargeability can be computed by two DC resistivity modelings with conductivity and .

##### 2.1. Finite Difference Method

Usually, there are two ways to discretize the differential equation: the finite difference method and the finite element method. The finite difference method was applied for 3D forward modeling in the present study. The matrix form of (1) can be written aswhere is the coefficient matrix; is the potential vector; and is a vector containing the locations of the positive and negative current sources. There are many numerical methods for solving linear equations; of these, the ILU0 preconditioner and the BICGSTAB iterative method [12] was chosen to solve (5) in serial algorithm. The ILU0 preconditioner is an incomplete LU factorization with zero fill-in [13]:where is the lower triangular and is the upper triangular.

##### 2.2. Forward Test

In order to test the accuracy of the forward algorithm, we compared our result of using the proposed BICGSTAB method with results from the preconditioned conjugate gradient (CGPC) algorithm [14]. The model in the present study consisted of a 10 Ω·m rectangular prism measuring 40 m × 40 m × 50 m embedded in a 100 Ω·m homogeneous half-space. The top of the model was located 30 m below ground surface. Four data acquisition methods were used for recording pole-pole data: surface-surface, surface-borehole, borehole-surface, and borehole-borehole (Figure 1). Their responses are shown in Figure 2, which indicates that the proposed BICGSTAB algorithm produced highly accurate results.