ERT2D: A 2D MATLAB Simulator for Electrical Resistance Tomography Alireza Aghasi
∗
Abstract ERT2D is a MATLAB simulator for electrical resistance tomography (ERT). The code provides the data vector and the model Jacobian matrix to be used in any gradient-based inversion scheme. The software is specifically designed for 2D geophysics examples, although with slight modifications the code would support other applications (e.g. medical). The simulator uses the finite difference technique to numerically solve the underlying partial differential equation (PDE) and employs an Adjoint technique to calculate the model Jacobian matrix.
1
Problem Setup
For an imaging domain Ω with surface boundary Γ, the forward model of interest motivated by several geophysical applications is ∇ ⋅ (σ(x)∇φ(x)) = s(x) α(x)
∂ φ(x) + β(x)φ(x) = 0 ∂n
in Ω,
(1)
on Γ.
(2)
In the context of electrical tomography as a basic modality, φ is the electrical potential, s the current distribution and σ is the electrical conductivity. All quantities are assumed to be functions of 2D spatial coordinates, x = [x y]T . Regarding the boundary condition (2), the notation ∂φ/∂n denotes the normal component of the gradient of φ on Γ, and α(x) and β(x) are functions defined on the surface Γ (i.e., x ∈ Γ) which are not simultaneously zero. A common geophysical problem associated with the model in (1)–(2) is shown in Figure 1. In this problem Γ consist of Γn , the interface between the earth and air ∗
School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA, USA Email contact:
[email protected]
1
Figure 1: The 2D modeling region for ERT. The imaging region of size `y by `x is surrounded by the sources and detectors. A Neumann boundary condition is considered for the top region (y = 0) and absorbing boundary condition is considered for the remaining boundaries. The dashed lines correspond to the dipole nodes used in every experiment. where we impose a zero current condition (β = 0, α = 1) and Γm , a surface with Robin boundary condition where the values α and β are chosen to model an infinite half-space (as if the imaging domain is extended to infinity from the left, right and bottom sides). A well established absorbing boundary condition (ABC) to model an infinite half space is the one proposed by Dey and Morrison [3]. For a current point source located at xs , i.e., s(x) = δ(x − xs ), the proposed ABC is in the form of [1] ∂ φ(x) + β(x)φ(x) = 0, ∂n where
(3)
(x − xs ) ⋅ n , x ∈ Γm (4) ∥x − xs ∥2 and n is the outward unit normal vector at x ∈ Γm . Clearly the ABC considered is dependent on the source position and therefore the boundary condition changes when we consider different source locations for different experiments. For simplicity here, we only consider real-valued conductivity (that is the ERT problems) although the approach can easily be adapted to handle complex valued conductivity (impedance tomography). β(x) =
2
1.1
Experiment Setup
The ERT test configuration for the problem is shown in Figure 1. The imaging domain of size `y × `x extends from −`x /2 to `x /2 in the x direction and from −`y to 0 in the y direction. The top side (y = 0) represents the air interface where a Neumann boundary condition is applied. At the remaining sides, absorbing boundary condition is applied to model an infinite half space. The sensor placement and configuration are also shown in the figure. In total ns sensors may be located on the periphery of the imaging domain. To make the set of electrical measurements denser, ne experiments may be carried out. In each experiment two sensors should act as the current dipole (connected via the lines in Figure 1) and the remaining sensors measure the corresponding electric potential. The dipoles are usually chosen in a cross medium configuration to make the data more sensitive to the objects within the central parts of the domain. It should be noted that ERT2D only supports electrical currents in the dipole form. With this setup the size of the total data vector becomes ne (ns − 2). Measurement data are generated by a finite difference solution of the Poisson’s equation, over the rectangular domain shown. For this purpose the imaging domain is discretized into m grids in the y direction and n grids in the x direction (Figure2). In fact, the proposed boundary condition in [3] models an infinite half space more accurately as the distance to the current source increases. For this purpose the program allows some grid extension in the x and y direction, here denoted as nx and ny . We would like to note that the grid extension in the y direction only takes place at the bottom side. As an implementation detail, it is discussed in [4] that to apply the aforementioned ABC, in a finite difference scheme a Dirichlet boundary condition may be considered and instead the conductivity of the boundary nodes may be scaled by β(x). ERT2D follows the same strategy in imposing the proposed ABC.
1.2
Adjoint Field Sensitivity Calculation
ERT2D may be employed by any gradient-based inversion scheme to invert for the conductivity values based on the observed data. ERT2D is equipped with the Adjoint technique to calculate the sensitivity of the sensor readings to the conductivity of each grid within the domain [2]. The details of the Adjoint technique are brought here as a general overview of the method. For s(x) = δ(x − xs ), i.e., a point source current at xs ∈ Ω, we denote by φs (x) the resulting potential over the domain Ω and consider the measured potential at xd ∈ Ω as φds = ∫ φs (x)δ(x − xd )dx. (5) Ω
3
Figure 2: The imaging domain is discretized into m grids in the y direction and n grids in the x direction. The program allows some grid extension in the x and y direction, here shown as nx and ny . To apply the proposed ABC, the conductivity values are scaled by β(x) at the boundary nodes. The scaling is shown by the varying values of the conductivity at the boundary nodes. Here the scaling corresponds to the current point source highlighted in the figure. As stated in [2], the variation of φds resulting from a perturbation δσ in the conductivity (i.e., the Fr´echet derivative of the measurements with respect to the conductivity) can be expressed as dφds [δσ] = ∫ δσ ∇φs ⋅ ∇φd dx, (6) dσ Ω where φd is the adjoint field that results from placing the current point source at xd . Using the superposition principle the result may be easily adapted to the dipole current source case. In a discrete representation, denoting ei as the electric potential at the i-th sensor resulted from a current point source at the j-th sensor, and σk the conductivity of the k-th pixel in the domain (k = 1, ⋯, nm), we have ∂ei = (∇φi ) ⋅ (∇φj ), ∂σk which would be central equation in calculating the Jacobian entries.
4
(7)
1.3
The MATLAB Code
ERT2D package consists of two main functions: 1. paramPackGenerator.m 2. ERT2D.m The package is also equipped with a demo file (Demo.m) for a sample simulation. The paramPackGenerator function generates a single struct as the output which contains matrices and variables that only need to be calculated once. This struct is used by the ERT2D to generate the data vector and the Jacobian matrix. The function is called in the following manner: ERTParams=paramPackGenerator(domainGridSize, sigmaBackground,... domainDimX, domainDimY, ExNo sides, ExNo bottom) The input variables are as follows: • domainGridSize: is the domain grid size as a 1 × 2 vector [m, n] where m and n are the vertical and horizontal number of grids, depicted in Figure 2. • sigmaBackground: is the background conductivity value, or basically the conductivity value near the boundaries. This is a positive scalar. The program assigns this conductivity to the extended grids. To apply the ABC this conductivity value is appropriately scaled at the extended boundaries to model an infinite half space. • domainDimX: is simply `x in meters. • domainDimX: is `y in meters. • ExNo sides: is the horizontal grid extension nex . • ExNo bottom: is the vertical grid extension ney . The output is ERTParams which contains several variables required by the ERT2D program. Many data files and matrices (e.g., the finite difference gradient matrices) that can be calculated only once are encapsulated inside ERTParams and enable ERT2D to be evaluated multiple times without a need to generate such matrices each time. This is of course a desirable feature in inversion programs that ERT2D needs to be called iteratively. The program also requires three setting files to be made available as .mat files:
5
• nodes.mat: is a matrix of size ns × 2 which contains the normalized x and y coordinates of each sensor in the domain. Note that the coordinates need to be normalized to the length of the domain in each direction. In the normalized domain the x values vary from −0.5 to 0.5 and the y values vary between −1 and 0. For instance if `x = 10 and `y = 6, a sensor positioned at the point [−5, −3] should be represented by the coordinates [−0.5, −0.5] in the data file nodes.mat. • dipole configuration.mat: contains the dipole architecture. This is a matrix of size ne × 2 that indicates which sensors are used as the dipole nodes at each experiment. For instance if a row of dipole configuration contains the numbers 5 and 7, it means the sensors on the fifth and seventh lines of the variable nodes are considered as dipoles in one of the experiments. • dipole voltages.mat: This is also a matrix of size ne × 2 that indicates the dipole voltages in each experiment. Normally, this is only a matrix of ±1 entries. The main program ERT2D is called in the following manner: [sens data Jacobian]=ERT2D(cond disb, evaluation, ERTParams) The input variables are as follows: • cond disb: is an m×n matrix which contains the values of electrical conductivity at each grid. The program is meant to calculate the sensor measurements and sensitivity matrix for such conductivity distribution. • evaluation: only takes values of 1 and 2. If 1, the output is only the sensor data and the Jacobian calculation is skipped (the Jacobian matrix is empty). If 2, both the sensor data as well as the Jacobian matrix are evaluated. • ERTParams: is the encapsulated variable generated by the paramPackGenerator function which allows passing all required inputs and settings via a single variable. The output variables are as follows: • sens data: is a vector of size ne (ns − 2) which contains all the sensor potential measurements for the given cond disb. The first ns −2 elements of sens data are the sensor readings in the first experiment, the second ns −2 elements correspond to the second experiment sensor readings and so on. • Jacobian: is a matrix of size ne (ns − 2) × nm where the j-th column corresponds to the sensitivity of the ne (ns − 2) sensor readings to the j-th pixel in the conductivity matrix. 6
Finally, Demo.m is a simple demo that models ERT in a given domain with some generic input files. The program illustrates the problem setting and shows the data for a test conductivity distribution. Also the Jacobian is evaluated and compared by once using infinitesimal perturbation of a certain pixel in the conductivity and once extracting the corresponding column of the Jacobian matrix generated by the program.
1.4
Parallel Processing
ERT2D code supports MATLAB parallel processing feature, which may significantly speed up the processing time. Enabling the parallel processing option is simply possible by applying the following command in the MATLAB command window: >> matlabpool open If MATLAB is able to successfully connect to multiple processing units, any calls to ERT2D afterwards would enjoy the parallel precessing feature.
Acknowledgement This program is based upon work supported by the National Science Foundation under grant no. EAR 0838313. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of theNational Science Foundation.
References [1] Alireza Aghasi, Parametric Shape Based Methods for Inverse Problems, PhD thesis, Tufts University, 2012. [2] A. Aghasi and E.L. Miller, Sensitivity calculations for Poisson’s equation via the adjoint field method, Geoscience and Remote Sensing Letters, IEEE, 9 (2012), pp. 237 –241. [3] A. Dey and HF Morrison, Resistivity modeling for arbitrarily shaped threedimensional structures, Geophysics, 44 (1979), pp. 753–780. [4] Jie Zhang, Randall L Mackie, and Theodore R Madden, 3-d resistivity forward modeling and inversion using conjugate gradients, Geophysics, 60 (1995), pp. 1313–1325.
7
