TY - JOUR
AU1 - Liu,, Bin
AU2 - Pang,, Yonghao
AU3 - Mao,, Deqiang
AU4 - Wang,, Jing
AU5 - Liu,, Zhengyu
AU6 - Wang,, Ning
AU7 - Liu,, Shenhua
AU8 - Zhang,, Xinxin
AB - SUMMARY 4-D electrical resistivity tomography (ERT), an important geophysical method, is widely used to observe dynamic processes within static subsurface structures. However, because data acquisition and inversion consume large amounts of time, rapid changes that occur in the medium during a single acquisition cycle are difficult to detect in a timely manner via 4-D inversion. To address this issue, a scheme is proposed in this paper for restructuring continuously measured data sets and performing GPU-parallelized inversion. In this scheme, multiple reference time points are selected in an acquisition cycle, which allows all of the acquired data to be sequentially utilized in a 4-D inversion. In addition, the response of the 4-D inversion to changes in the medium has been enhanced by increasing the weight of new data being added dynamically to the inversion process. To improve the reliability of the inversion, our scheme uses actively varied time-regularization coefficients, which are adjusted according to the range of the changes in model resistivity; this range is predicted by taking the ratio between the independent inversion of the current data set and historical 4-D inversion model. Numerical simulations and experiments show that this new 4-D inversion method is able to locate and depict rapid changes in medium resistivity with a high level of accuracy. Electrical resistivity tomography (ERT), Downhole methods, Inverse theory, Numerical modelling, Time-series analysis 1 INTRODUCTION Time-lapse or 4-D electrical resistivity tomography (ERT) has been widely applied in research fields including environmental monitoring (Doetsch et al.2012; Loke et al.2014; Kuras et al.2016), embankments (Case 2012; Song et al.2016) landslides (Wilkinson et al. 2010, 2016; Uhlemann et al.2017), and hydro geophysics (Jayawickreme et al. 2010; Hilbich et al.2011; Hermans et al.2016). The technique has the advantage of providing images of subsurface resistivity distribution that change dynamically with time. Numerous scholars have conducted studies on time-lapse and 4-D resistivity inversion methods. For example, Daily et al. (1992, 2004) proposed a ratio inversion method using initial and proceeding data sets; this is an effective approach at highlighting changes in resistivity with respect to the background. This method was proven to be useful for monitoring soil-moisture dynamics. Loke (1999) and Tsourlos et al. (2003) used the model obtained from the inversion of the initial data set to constrain the inversion of the proceeding time-lapse data sets. This method achieved good results in monitoring the leachate generation and transport. LaBrecque & Yang (2001) added difference-constrained items to the proposed objective functions to minimize the differences between the two data sets and two inversion models; this caused the inversion to converge faster and the measured data to fit well. Hayley et al. (2011) proposed a simultaneous inversion method, which helped improve the model resolution and computational efficiency. In these methods, the time-lapse inversion usually assumes that the static subsurface model is time-invariant over a single acquisition cycle, indicating that changes in the subsurface media can be ignored during data acquisition over a single time step. However, significant changes, such as rapid fluid migration in highly permeable media, can occur in the subsurface medium during this process. In such cases, the changes in material properties during data acquisition cannot be ignored, as the resulting subsurface images become severely distorted (Kim et al.2009). Owing to the time-consuming inversion process and distorted images, it is difficult to capture the rapid changes in the underground medium in a timely and dynamic manner. In response to this issue, some scholars have studied the use of original data directly for monitoring imaging instead of resistivity inversion. Rucker et al. (2014) conducted the underground injection monitoring experiment by using the data of raw current and voltage (then converted them to transfer resistance and apparent resistivity). This method can be used to monitor the subsurface changes dynamically and qualitatively. The results showed that there was a good match in the direction and coverage of underground injection migration. However, in-depth inversion studies are needed to obtain physical parameters (such as resistivity) and more accurately describe the spatial location of the changes in the subsurface medium and the migration process over time. Recently, Loke et al. (2014) proposed a 4-D smoothness-constrained inversion method to handle time-dependent survey configurations using an L-curve method. Kim (2005) and Kim et al. (2009) proposed a 4-D inversion algorithm that simultaneously inverts entire geophysical monitoring data sets in which both the subsurface model and monitoring data are defined in a 4-D space–time domain. Regularization is used during the inversion in both the spatial and time domains to effectively reduce inversion artefacts and improve the stability of the inverse problem. Karaoulis et al. (2011, 2014) and Kim et al. (2013) studied the regularization parameters in 4-D resistivity inversion equations in more detail and proposed a scheme to dynamically adjust the time-regularization coefficient according to the range of the change in spatial resistivity within various monitoring stages; this technique provides a theoretical foundation for solving problems concerning rapid changes. Bellmunt et al. (2016) proposed an index, known as the in-panel/off-panel sensitivity, to evaluate the capability of each electrode array to resolve the plume migration direction and further optimized monitoring data sets suitable for rapid change processes. This paper focuses on the 4-D monitoring of rapid resistivity changes in underground media and proposes an improved 4-D resistivity inversion method. First, the new algorithm divides complete data sets from each acquisition cycle into multiple subsets of data by using multiple reference time points. These subsets of data are then restructured and dynamically substituted into the improved 4-D inversion equation (which applies increased weights to new data and actively adjusts the time-regularization coefficients) for the inversion calculation. Then, the GPU parallel algorithm is used to optimize the 4-D inversion equation processes. Finally, our method was verified via two sets of numerical simulation tests and a fast solute transport monitoring experiment. 2 METHODOLOGY 2.1 Novel scheme for the restructuring and inversion of resistivity monitoring data The use of 4-D resistivity inversion equations for resistivity monitoring data usually allows the simultaneous inversion of data obtained at different times in just one inversion step. However, because the acquisition and inversion of large 4-D data sets are highly time consuming, it is difficult to capture rapid changes occurring in the medium during an acquisition cycle via 4-D inversion, resulting in delayed interpretations and reduced image resolution. To solve this problem, the authors propose a scheme for the restructuring and inversion of continuously measured data sets. Multiple reference time points were selected from each acquisition cycle and used to divide the data set into multiple segments for active extraction. The extracted data sets were then restructured using the preceding data sets and substituted into the 4-D inversion equation. Thus, this method achieves rapid inversion and imaging during the monitoring processes through the execution of 4-D inversions at every reference time point, thus enhancing the responsiveness of the 4-D inversion to rapid resistivity changes in the underground medium. In most 4-D inversions, a single inversion is performed on a number of complete data sets obtained over several successive acquisition cycles; this produces a corresponding number of geoelectric models. In this work, each acquisition cycle is divided into multiple segments using multiple reference time points. Hence, the new data being incorporated into the inversion in each instance represent only a fraction of the complete acquisition cycle. The amount of data being used in these inversions is relatively small, resulting in poor-quality inversion models. To ameliorate this issue, our scheme restructures the data sets prior to inversion. In this process, subsets of data acquired between the reference time points of adjacent acquisition cycles in the same sequence are restructured to produce multiple data sets with consistent data structures (i.e. with the same ABMN electrode array configuration) and data volumes, allowing the model sensitivity structure to remain constant across a time-series (Robinson et al.2018). Fig. 1 illustrates our measurement data extraction and restructuring scheme, where |${t_k}$| represents a point in time and |$\Delta t$| represents a time interval. Each colour represents a complete data set obtained over a single acquisition cycle during continuous data acquisition. For example, a complete measurement data set |${{\boldsymbol{D}}^i}(i = 1,2,...)$| could consist of |${t_{k - 5}}\sim {t_{k - 1}}$| or |${t_{k - 1}}\sim {t_{k + 3}}$| in Fig. 1. In this example, four reference time points are selected in each acquisition cycle and used to split each measurement data set (⁠|${{\boldsymbol{D}}^i}$|⁠) into four equal subsets of data, which are denoted as. |${\boldsymbol{D}}{{\boldsymbol{s}}^k}(k = 1,2,...)$| In this way, any set of four adjacent subsets of data can be restructured to form a new data set with a complete electrode array. Here, a restructured data set is defined as |${\boldsymbol{D}}{{\boldsymbol{r}}^j}(j = 1,2,...)$|⁠. For example, the |${t_{k - 4}}\sim {t_k}$|⁠, |${t_{k - 3}}\sim {t_{k + 1}}$|⁠, and |${t_{k - 2}}\sim {t_{k + 2}}$| data sets in Fig. 1 correspond to three restructured data sets. Figure 1. Open in new tabDownload slide Illustration of the measurement data set splitting, extraction and restructuring scheme. Each group with identical colours represents a complete data set obtained from a single continuous acquisition cycle. In this example, we split each data set into four parts. As shown in the lower left corner, each rectangular block that carries an inscribed circle here represents a subset of data, |${{\boldsymbol D}}{s^k}( {k = 1,2,...} )$|⁠. Subsets of data with identical inscribed circles have the same measurement arrays. Restructured data sets, |${\boldsymbol{D}}{{\boldsymbol{r}}^j}(j = 1,2,...)$|⁠, are shown in the bottom-right section. Figure 1. Open in new tabDownload slide Illustration of the measurement data set splitting, extraction and restructuring scheme. Each group with identical colours represents a complete data set obtained from a single continuous acquisition cycle. In this example, we split each data set into four parts. As shown in the lower left corner, each rectangular block that carries an inscribed circle here represents a subset of data, |${{\boldsymbol D}}{s^k}( {k = 1,2,...} )$|⁠. Subsets of data with identical inscribed circles have the same measurement arrays. Restructured data sets, |${\boldsymbol{D}}{{\boldsymbol{r}}^j}(j = 1,2,...)$|⁠, are shown in the bottom-right section. During a 4-D resistivity inversion performed using our scheme, all of the restructured data sets obtained before some given point in time were substituted into the inversion equation, and a simultaneous inversion was performed. For example, at time |${t_k}$|⁠, the |$[ {{\boldsymbol{D}}{{\boldsymbol{r}}^1},{\boldsymbol{D}}{{\boldsymbol{r}}^2},...,{\boldsymbol{D}}{{\boldsymbol{r}}^j}} ]$| data set was substituted into the inversion equation, producing j inverted images. This was repeated at each reference time point. Hence, a 4-D inversion was performed after the acquisition of each |${{\boldsymbol D}}{s^k}$| subset of data. 2.2 Improved 4-D resistivity inversion algorithm 2.2.1 Improved 4-D resistivity inversion equation The objective function of the conventional 4-D resistivity inversion method may be expressed as (Zhang et al.2005; Kim et al.2009) $$\begin{eqnarray*} {\boldsymbol \Phi } = {{\boldsymbol \varphi }_d} + \lambda {{\boldsymbol \varphi }_m} + \mu {{\boldsymbol \varphi }_t}, \end{eqnarray*}$$ (1) where |${{\boldsymbol \varphi }_d}$| is the data misfit term, |${{\boldsymbol \varphi }_m}$| and |${{\boldsymbol \varphi }_t}$| are the spatial and temporal regularizations of the model parameters, respectively, and |$\lambda $| and |$\mu $| are the regularization coefficients that balance the contributions of the regularization terms to the objective function. In this work, |$\lambda $|= 0.04 and |$\mu $|= 1. The expression of the data misfit is $$\begin{eqnarray*} {{\boldsymbol \varphi }_d} = {({W_d}{\boldsymbol{e}})^T}{W_d}{{\boldsymbol e}}, \end{eqnarray*}$$ (2) where e is the misfit vector between the actual observed data and theoretical data obtained from forward modelling and is defined as |${{\boldsymbol e}} = {\boldsymbol{Dr}} - ({\bf G}({{\boldsymbol M}}) + J{{\boldsymbol \Delta}}{{\boldsymbol M}})$|⁠, where |${{\boldsymbol{Dr}}} = [ {{\boldsymbol{D}}{{\boldsymbol{r}}^1},{\boldsymbol{D}}{{\boldsymbol{r}}^2},...,{\boldsymbol{D}}{{\boldsymbol{r}}^{Nj}}} ]$| is the data vector and |${N_j}$| is the current number of restructured data sets. |${\bf G}$| represents the forward modelling response, and |${{\boldsymbol M}} = {[{{{\boldsymbol M}}^1},{{{\boldsymbol M}}^2}, \cdots ,{{{\boldsymbol M}}^{{N_j}}}]^T}$| represents the model parameters of the |${N_j}$| data sets, where |${{{\boldsymbol M}}^{_j}}(j = 1,2,...{N_j})$| has a one-to-one correspondence with |${\boldsymbol{D}}{{\boldsymbol{r}}^j}$|⁠. Furthermore, |${\boldsymbol{\Delta M}}$| represents the model perturbation vector for the iterations of the inversion, and J is the extended Jacobian matrix, where |$J = {\boldsymbol {diag}}[{J^1},{J^2}, \cdots ,{J^{{N_j}}}]$|⁠; |${J^j}(j = 1,2, \cdots ,{N_j})$| is the Jacobian matrix for model parameter |${{\boldsymbol{M}}^j}$|⁠. Here, data-weighting matrix |${W_d}$| is added to the data misfit term, and is expressed as $$\begin{eqnarray*} {W_d} = {\rm{diag}}[1, \cdots ,1,\overbrace {{\omega _d}, \cdots ,{\omega _d}}^{{{\boldsymbol Ns}}}], \end{eqnarray*}$$ (3) where |$Ns$| is the number of data points contained in each subset of data in the dynamic updating process, and the elements of |${W_d}$| represent the weights of the measured data in the 4-D inversion equation. In our dynamic measurement data extraction and restructuring scheme, the measurement data most recently inserted into the inversion (subset of data containing |$Ns$| data points) contain information about the changes that are occurring below the surface. The effects of new data on the 4-D inversion may be amplified or attenuated by modulating the data weights (i.e. by changing the value of |${\omega _d}$|⁠). The spatial regularization terms proposed by Kim et al. (2009) are $$\begin{eqnarray*} {{\boldsymbol \varphi }_m}{\rm{ = (}}{\partial ^2}{\boldsymbol{\Delta M}}{)^T}{\rm{(}}{\partial ^2}{\boldsymbol{\Delta M)}}. \end{eqnarray*}$$ (4) Karaoulis et al. (2011) improved temporal regularization |${{\boldsymbol \varphi }_t}$| by adding an active time-constrained diagonal matrix A, which is based on the predicted range of the resistivity changes in the model. The new time regularization is now represented by eq. (5), where |${C_T}$| is a matrix that constrains model smoothness in the time domain, as the use of fixed time-regularization coefficients may suppress the inversion response to rapid changes in the real medium resistivity: $$\begin{eqnarray*} {{\boldsymbol \varphi }_t} = \left\{ {{C_T}} \right.({{\boldsymbol M}}{\left. { + {{\boldsymbol {\Delta M}}})} \right\}^T}A{C_T}({{\boldsymbol M}} + {\boldsymbol{\Delta M}}). \end{eqnarray*}$$ (5) The minimization of objective functions (1)–(5) yields the following 4-D inversion equation (eq. 6), where C is the space regularization matrix, $$\begin{eqnarray*} &&{\rm{(}}{J^T}{W_d^T}{W_d}J + \lambda {C^T}C + \mu {C_T^T}{AC_T}{\rm{)}}{\boldsymbol{\Delta M}} \nonumber \\ &&\quad = {J^T}{W_d^T}{W_d}({\boldsymbol{Dr}} - {\bf G}{\rm{(}}{{\boldsymbol M}}{\rm{))}} - \mu {C_T^T}{AC_T}\left. {{\boldsymbol M}} \right\}. \end{eqnarray*}$$ (6) In this work, the preconditioned conjugate gradient method was used to solve the aforementioned inversion equation, which produces |${\boldsymbol{\Delta M}}$|⁠. |${\boldsymbol{M}}$| is then iteratively updated to obtain the final inversion result. In inversion processes, the result of the current inversion is usually used as the initial model in the next inversion. In other words, if the inversion at time |${t_k}$| is |${\boldsymbol{M}} = {[ {{{{{\skew{4}\hat{\boldsymbol M}}}}^1},{{{{ \skew{4}\hat{\boldsymbol M}}}}^2}, \cdots {{{{ \skew{4}\hat{\boldsymbol M}}}}^{Nj}}} ]^T}$|⁠, the initial model at time |${t_k}_{ + 1}$| is then |${[ {{{{{\skew{4}\hat{\boldsymbol M}}}}^1},{{{{\skew{4}\hat{\boldsymbol M}}}}^2}, \cdots {{{{\skew{4}\hat{\boldsymbol M}}}}^{Nj}},{{{{\skew{4}\hat{\boldsymbol M}}}}^{Nj}}} ]^T}$|⁠. Here, |${{{ \skew{4}\hat{\boldsymbol M}}}^{Nj}}$| is the initial value of both |${{{\boldsymbol M}}^{Nj}}$| and |${{{\boldsymbol M}}^{Nj + 1}}$|⁠. The use of known inversion results as reference models for future inversions helps to reduce the non-uniqueness of 4-D inversions. 2.2.2 Prediction of the range of resistivity changes in the model In the method by Karaoulis et al. (2011), one inversion iteration is performed on each independent data set, and the range of the resistivity changes is then estimated in each region of the model by comparing the inversion results (i.e. the resistivity distributions). This information is then used to adjust the values of the time-regularization coefficients. Based on this algorithm, we propose an improved method for adjusting the time-regularization coefficients, in which the time dependence of the model resistivity is determined by calculating the ratio of the independent inversion of the latest restructured data set to the historical 4-D inversion. This time dependence is then used to determine the active time-constrained diagonal matrix A for the next time point. Our reasoning is that the 4-D inversion of the previous time step is readily available at every time step after the application of our data extraction and restructuring scheme; therefore, the independent inversion of the previous time step is unnecessary. In a 4-D resistivity inversion, all of the data sets are treated as a single complete data set, and the inversion is performed by substituting all of the data sets into a single inversion equation; regularization constraints are then applied on the models at each time step using time-regularization coefficients. The changes in resistivity over time obtained via 4-D resistivity inversion are much more reliable than those derived from the independent inversion of each time step. The specific details of our scheme are as follows. Fig. 1 shows only one data set at |${t_4}$|⁠. An independent inversion was performed to obtain the result corresponding to |${t_4}$|⁠: |${{\boldsymbol{\skew{4}\hat{M}}}^4}$|⁠. For the subsequent |${t_{\rm{k}}}_{}$| time points |$(k = 5,6,...)$|⁠, new restructured data sets |${\boldsymbol{D}}{{\boldsymbol{r}}^j}$| were formed from the newly acquired subsets of data and the previous data sets. As the inversion results from the previous time steps ( |${{{\skew{4}\hat{\boldsymbol M}}}^4} - {{{\skew{4}\hat{\boldsymbol M}}}^{k - 1}}$| resistivity distributions corresponding to |${t_4}\sim {t_{k - 1}}$|⁠) were already available, the results of all time points from |${t_4}$| to |${t_{\rm{k}}}$| (⁠|${{\boldsymbol{\skew{4}\hat{M}}}^4} - {{\boldsymbol{\skew{4}\hat{M}}}^{\rm{k}}}$|⁠) can be obtained simply by performing an independent inversion on the newly obtained |${\boldsymbol{D}}{{\boldsymbol{r}}^j}$| to produce the result corresponding to |${t_{\rm{k}}}$|⁠: |${{\boldsymbol{\skew{4}\hat{M}}}^{\rm{k}}}$|⁠. The range of resistivity changes at each time position is then estimated using these inversions, which are subsequently used to determine the active time-constrained diagonal matrix A. All of the currently existing restructured data sets (⁠|${\boldsymbol{Dr}} = [ {{\boldsymbol{D}}{{\boldsymbol{r}}^1},{\boldsymbol{D}}{{\boldsymbol{r}}^2},...,{\boldsymbol{D}}{{\boldsymbol{r}}^j}} ]$|⁠) are then substituted into the 4-D inversion equation. The results of this inversion are subsequently used to determine the active time-constrained diagonal matrix A for the 4-D inversion in the next time step. In this study, we adopted the value scheme for A reported by Karaoulis et al. (2011). 2.3 Fast GPU parallelization-based algorithm for solving inversion equations All of the restructured data sets in our new 4-D resistivity inversion method must be simultaneously substituted into the inversion equation. However, this effectively multiplies the inversion computational load. Therefore, the computational efficiency of the 4-D inversion must be improved. GPU-parallelized inversions have already been explored for use in fields such as seismological exploration (Liu et al.2009), gravity exploration (Chen et al.2012), and full-waveform inversion (Liu et al.2015), and are highly effective in improving the efficiency of inversion calculations. In this work, CUDA-Fortran based GPU parallelization code was used to improve and optimize the inversion program by exploiting the characteristics of the finite element method and our 4-D resistivity inversion algorithm. Solving the linear equations in eq. (6) is the most time-consuming step in our 4-D inversion algorithm. During monitoring processes, the continuous acquisition of the measured data increases the dimensionality of the |$( {{J^{{\rm T}} }{W_d^T}} {W_d}J + \lambda {C^T}C + \mu {C_T^T}{AC_T} )$| coefficient matrix by several folds, commensurately increasing the computational workload. As the coefficient matrix is a symmetric positive-definite matrix, we adopted the Jacobi preconditioned conjugate gradient (JPCG) method, which is an iterative method for numerically solving large sparse linear systems. The JPCG method is well suited for GPU implementation as it is highly parallelizable and converges very quickly. Although the iterations of an inversion must be computed in series, large numbers of addition, multiplications, and inner product operations must be performed between the sparse matrices and vectors during each iteration. Hence, parallelization of these operations significantly improves the computational efficiency of the inversion equation. To illustrate our method in further detail, let us suppose that there is a symmetric positive-definite sparse linear system such that $$\begin{eqnarray*} {A_J}{\boldsymbol{x}} = {\boldsymbol{b}}, \end{eqnarray*}$$ (7) where |${A_J}$| is an n × n positive-definite sparse matrix, |${\boldsymbol{b}}$| is a constant vector and |${\boldsymbol{x}}$| is the vector being solved. As this equation is equivalent to eq. (6), $$\begin{eqnarray*} \left\{ \begin{array}{@{}l@{}} {A_J} = {J^T}{W_d^T}{W_d}J + \lambda {C^T}C + \mu {C_T^T}{AC_T}\\ {\boldsymbol{x}} = {\boldsymbol{\Delta M}}\\ {{\boldsymbol b}} = {J^T}{W_d^T}{W_d}{\rm{(}}{\boldsymbol{Dr}} - {\bf G}{\rm{(}}{{\boldsymbol M}}{\rm{))}} - \mu {C_T^T}{AC_T}{\boldsymbol{M}} \end{array} \right.. \end{eqnarray*}$$ (8) In the preconditioned conjugate gradient algorithm, an initial solution for the vector is first substituted into the iteration to transform the inversion equation into a series of linear combinations between vectors or matrices and matrix-vector products. All of the linear algebraic operations in this algorithm can be implemented in parallel by the CUDA Basic Linear Algebra Subprograms (CUBLAS) and the CUDA basic linear algebra subroutines for SPARSE matrices (CUSPARSE) provided by CUDA (Wolfe 2012; Ruetsch & Fatica 2014). In addition, we have adopted certain strategies to optimize the inversion computations in our program. First, matrix multiplications are converted into matrix-vector multiplications by redefining eq. (7) as |${A_J}x = {({W_d}{{\boldsymbol J}}{\rm{)}}^T}{\rm{((}}{W_d}J{\rm{)}}{\boldsymbol{x}}{\rm{)}} + \lambda {C^T}{\rm{(}}C{{\boldsymbol x}}{\rm{)}} + \mu {C_T^T}( {A{\rm{(}}{C_T}{{\boldsymbol x}}{\rm{)}}} )$| to minimize the generation of matrix multiplications. Some of the matrix multiplications have been avoided by switching the order of the calculations, which greatly reduces the computational workload. Second, large sparse matrices such as |${W_d}$|⁠, C, and J are stored as compressed 1-D arrays in the program. This saves storage space and reduces or eliminates the wastage of computational time on zero elements, reducing the computational time of the inversion algorithm. The hardware environment used in this study consisted of an NVIDIA GeForce GTX960 graphics card (or graphics processing unit, GPU) and an Intel® Core™ i5-6500 CPU running at a frequency of 3.2 GHz. This GPU supports Compute Capability 5.2 and is equipped with 1024 stream processors; it has a core frequency of 1051 MHz, a video memory of 2 GB, and a memory frequency of 7012 MHz. The memory bandwidth of the GPU is 232 GB s−1. Comparisons of time consumption before and after parallelization for dozens of numerical simulations using different amounts of data show more than 20-fold (on average) improvements in inversion computational efficiency after GPU parallelization. 3 NUMERICAL EXPERIMENTS 3.1 Impact of data-weighting coefficients A set of geoelectric models was devised for a rapidly changing low-resistivity body in order to examine the influence of data-weighting coefficients on the inversion results. |${G_1}{\rm{\sim }}{G_5}$| geoelectric models corresponding to nine time points (⁠|${t_0}\sim {t_9}$|⁠) are shown in Fig. 2. A 2-D finite element model was established using a quadrilateral mesh, and the inversion was carried out in an 8 m × 14 m (length (X) × depth (Z)) region. The edges of the grids were defined as 0.5 m × 0.5 m, and the background resistivity of the model was 800 Ω m. At the initial time point, the geoelectric model contained only a single low-resistivity body (1.0 m (length) × 2.0 m (height)) with a resistivity of 50 Ω m. After some time, a rapid change occurred in the model at |${t_4}$|⁠, and the low-resistivity body began to move quickly toward the positive X–negative Z direction (i.e. a rapid change occurred toward the bottom right direction of the model). Both the location and size of this mass changed over time. The shapes and sizes of the model at each reference time point are shown in Fig. 2, in which |${G_2}{\rm{\sim }}{G_5}$| correspond to reference time points |${t_5}\sim {t_8}$|⁠. The resistivity of the low-resistivity body was assumed to be constant during these rapid changes. Data were acquired during this simulation via 2-D cross-hole ERT in the space–time observation mode in which 28 electrodes (spaced 0.5 m apart) were installed on each survey line. We used a measurement array known as the Bipole_Bipole (AM-BN), where a complete data set contains 5625 independent data points. Figure 2. Open in new tabDownload slide A set of geoelectric models for a low-resistivity body experiencing rapid changes. The background resistivity is 800 Ω m and the resistivity of the abnormal body that changes to the lower right is set to 50 Ω m. The model remains stable until time |${t_4}$|⁠. Starting from time |${t_4}$|⁠, for every time interval, that is |$\Delta t = {t_{i + 1}} - {t_i}$|⁠, the abnormal body moves a rectangular block of the same resistivity value to the lower right. The rectangular block that carries an inscribed circle below represents a subset of data collected during each time interval. Figure 2. Open in new tabDownload slide A set of geoelectric models for a low-resistivity body experiencing rapid changes. The background resistivity is 800 Ω m and the resistivity of the abnormal body that changes to the lower right is set to 50 Ω m. The model remains stable until time |${t_4}$|⁠. Starting from time |${t_4}$|⁠, for every time interval, that is |$\Delta t = {t_{i + 1}} - {t_i}$|⁠, the abnormal body moves a rectangular block of the same resistivity value to the lower right. The rectangular block that carries an inscribed circle below represents a subset of data collected during each time interval. In this numerical simulation, each c≈omplete data set was partitioned into four subsets of data, which were defined as |${\boldsymbol{D}}{{\boldsymbol{s}}^k}$| (⁠|$k = 1,2,...,8$|⁠), as shown in Fig. 2. A restructured data set (⁠|${\boldsymbol{D}}{{\boldsymbol{r}}^j}( {j = 1,2,...,5} )$|⁠) was then formed by an arbitrary set of four consecutive subsets of data. First, the |${{\boldsymbol D}}{{{\boldsymbol r}}^1} = [{\boldsymbol{D}}{{\boldsymbol{s}}^1},{{\boldsymbol D}}{{{\boldsymbol s}}^2},{{\boldsymbol D}}{{{\boldsymbol s}}^3},{{\boldsymbol D}}{{{\boldsymbol s}}^4}]$| restructured data set was obtained at time |${t_4}$| via forward numerical modelling. This data set was then substituted into the 4-D inversion equation to produce the |${t_4}$| inversion model. After the |${{\boldsymbol D}}{{{\boldsymbol s}}^5}$| subset was acquired at |${t_5}$|⁠, the |${{\boldsymbol D}}{{{\boldsymbol r}}^2} = [{{\boldsymbol D}}{{{\boldsymbol s}}^2},{{\boldsymbol D}}{{{\boldsymbol s}}^3},{{\boldsymbol D}}{{{\boldsymbol s}}^4},{{\boldsymbol D}}{{{\boldsymbol s}}^5}]$| data set was obtained by restructuring the subsets of data. Again, the restructured data set was substituted into the 4-D inversion to obtain the inversions corresponding to |${t_4}$| and |${t_5}$|⁠. A 5 per cent Gaussian noise was added to the forward modelling data to produce a more realistic data set. To limit the number of variables in this simulation to one, a constant value of 0.1 was selected for the time-regularization coefficient, |${a_i}$|⁠. Six iterations were performed for each inversion. In this way, the five inversion images corresponding to |${t_4}\sim {t_8}$| were obtained by |${t_8}$|⁠. Four different data-weighting coefficients (⁠|${\omega _d}$|= 1.0, 1.6, 2.2 and 5.0) were selected for substitution into the 4-D inversion equation, producing the 7 × 5 (row × column) set of inverted images shown in Fig. 3. Each row in this figure represents the 4-D inversions obtained using a given data-weighting coefficient, while each column shows the results corresponding to a given time point with different data-weighting coefficients. Figure 3. Open in new tabDownload slide Inversion results of the geoelectric models in Fig. 2 with different data-weighting coefficients (⁠|${\omega _d}$|⁠). Each row represents the inverted images obtained with a given |${\omega _d}$|⁠, while each column shows the results obtained for a given time point using different |${\omega _d}$| values. The multiple black-edged rectangles represent the state of the actual geoelectric models at each time step. Figure 3. Open in new tabDownload slide Inversion results of the geoelectric models in Fig. 2 with different data-weighting coefficients (⁠|${\omega _d}$|⁠). Each row represents the inverted images obtained with a given |${\omega _d}$|⁠, while each column shows the results obtained for a given time point using different |${\omega _d}$| values. The multiple black-edged rectangles represent the state of the actual geoelectric models at each time step. We determined that the inverted images in the first column of Fig. 3, which correspond to the inversion results using different data-weighting coefficients, do not differ greatly from each other. This is reasonable as the distribution of resistivity in the underground medium is static before |${t_4}$|⁠. In these inverted images, the location and shape of the low-resistivity body are well defined and generally consistent with the original geoelectric model, |${G_1}$|⁠. As rapid changes occur at |${t_4}$| in the low-resistivity body, the |${\boldsymbol{D}}{{\boldsymbol{s}}^5}$| subset of data contains the most recent information about the responses in resistivity to these geological changes. This response information is enhanced by the use of different data-weighting coefficients at |${t_5}$|⁠. The inverted images at |${t_5}$| (the second column in Fig. 3) show that the use of different data-weighting coefficients begins to cause divergence in the imaging of rapid changes in the low-resistivity body within this model. As the data-weighting coefficients increase, the position and shape of the low-resistivity body in the inverted images (as the body migrates along the bottom-right direction) become increasingly similar to those in the |${G_2}$| geoelectric model. Although the inversions with |${\omega _d}$| = 1.0 and |${\omega _d}$| = 1.6 do show the low-resistivity body moving in the bottom-right direction, the spatial positions of the low-resistivity body at |${t_6}$| in these inversions are incorrect. Hence, smaller data-weighting coefficients result in inversions that reflect the trends in underground resistivity changes but provide incorrect depictions of the location and shape of the low-resistivity anomaly. These errors hinder the interpretation and imaging of the actual changes occurring below the surface. Next, we analysed the shapes of the low-resistivity anomalies in the inverted images. Columns 4 and 5 in Fig. 3 show that the |${\omega _d}$| = 2.2 inversion describe the position of the rapidly changing low-resistivity anomaly with a high level of accuracy. However, in the |${\omega _d}$| = 5.0 inversion, the inverted shapes of the low-resistivity anomaly have changed, increasing more than those in |${G_2}$|⁠. This indicates that the inversion is most effective when the data-weighting coefficients are restricted to a certain range. When excessively large weights are assigned to the dynamically updated data sets, the fluctuations in the overall data set increase dramatically, and the data misfit terms in the objective function can have very large weights, causing these terms to dominate the inversion. This will reduce the effectiveness of the time and space regularizations and ultimately degrade the resolution of the inverted images. In summary, an increase in the weights of the new subsets of data in the inversion objective function amplifies the change-related response signals in the data. This is beneficial for the early detection of anomalous changes in underground media and the timely capture of predictive information. The example described in this section reveals that appropriate data weighting accelerates the detection of rapid changes in underground resistivity distributions compared to inversions without data weighting ( |${\omega _d}$| = 1.0). 3.2 Rapid change models A new set of geoelectric models was prepared for a rapidly changing low-resistivity body to examine the ability of our 4-D inversion equation to image medium rapid resistivity changes. Fig. 4 illustrates the |${G_1}\sim {G_3}$| geoelectric models for this low-resistivity body corresponding to nine time points. The inverted region was 6 m × 14 m (horizontal width (X) × depth (Z)), and the initial geoelectric model (⁠|${G_1}$|⁠) contained only a small low-resistivity body. The model mode of observation and resistivity settings were the same as those used in the previous example. The model begins to change rapidly at |${t_4}$|⁠, and the low-resistivity body moves quickly in the negative z-direction (i.e. toward the lower part of the model). This process continues until |${t_6}$|⁠, at which point the position and size of the low-resistivity body both change; the geoelectric model corresponding to |${t_6}$| is |${G_2}$|⁠. A second rapid change occurs in the model at this instant (⁠|${t_6}$|⁠), causing the low-resistivity body to migrate rapidly toward the bottom right. The model then stabilizes after |${t_7}$| (which corresponds to geoelectric model |${G_3}$|⁠). The resistivity of the body remains constant (at 50 Ω m) throughout these changes. Figure 4. Open in new tabDownload slide A set of geoelectric models for a low-resistivity body experiencing rapid changes. The background resistivity is 800 Ω m and the resistivity of the abnormal body that changes downwards and to the lower right is set to 50 Ω m. The model remains stable until time |${t_4}$|⁠. From times |${t_4}$| to |${t_7}$|⁠, the abnormal body moves rapidly downwards and to the lower right, and the model does not change after time |${t_7}$|⁠. The rectangular block that carries an inscribed circle below represents a subset of data collected during each time interval. Figure 4. Open in new tabDownload slide A set of geoelectric models for a low-resistivity body experiencing rapid changes. The background resistivity is 800 Ω m and the resistivity of the abnormal body that changes downwards and to the lower right is set to 50 Ω m. The model remains stable until time |${t_4}$|⁠. From times |${t_4}$| to |${t_7}$|⁠, the abnormal body moves rapidly downwards and to the lower right, and the model does not change after time |${t_7}$|⁠. The rectangular block that carries an inscribed circle below represents a subset of data collected during each time interval. In this example, each complete data set was similarly segmented into four subsets of data, and a new subset of data |${\boldsymbol{D}}{{\boldsymbol{s}}^k}$| (k = 1, 2, L, 8 in this case) was acquired for each (⁠|${t_i} - {t_i}_{ - 1}$|⁠) time interval. First, a restructured data set |${{\boldsymbol D}}{{{\boldsymbol r}}^1} = [{{\boldsymbol D}}{{{\boldsymbol s}}^1},{{\boldsymbol D}}{{{\boldsymbol s}}^2},{{\boldsymbol D}}{{{\boldsymbol s}}^3},{{\boldsymbol D}}{{{\boldsymbol s}}^4}]$| was obtained at |${t_4}$| through forward numerical modelling. An independent inversion was performed on this data set to produce the inversion at |${t_4}$|⁠, which corresponds to the |${G_1}$| geoelectric model. |${{\boldsymbol D}}{{{\boldsymbol s}}^5}$| and |${{\boldsymbol D}}{{{\boldsymbol s}}^6}$| were obtained during |${t_4}\sim {t_6}$|⁠. The restructuring of these subsets of data produced |${{\boldsymbol D}}{{{\boldsymbol r}}^2} = [{{\boldsymbol D}}{{{\boldsymbol s}}^2},{{\boldsymbol D}}{{{\boldsymbol s}}^3},{{\boldsymbol D}}{{{\boldsymbol s}}^4},{{\boldsymbol D}}{{{\boldsymbol s}}^5}]$| and |${{\boldsymbol D}}{{{\boldsymbol r}}^3} = [{{\boldsymbol D}}{{{\boldsymbol s}}^3},{{\boldsymbol D}}{{{\boldsymbol s}}^4},{{\boldsymbol D}}{{{\boldsymbol s}}^5},{{\boldsymbol D}}{{{\boldsymbol s}}^6}]$|⁠, which were substituted into the 4-D inversion equation to yield the inversions corresponding to |${t_5}$| and |${t_6}$|⁠. Note that beginning at |${t_5}$|⁠, the previous 4-D inversion results and the independent inversion of the latest restructured data set were used to determine the weight of the time-regularization constraint, as described in Section 2.2.2 This regularization was then used in the current 4-D inversion. New 4-D inversions were performed using restructured data sets formed by the |${{\boldsymbol D}}{{{\boldsymbol s}}^7}\sim {{\boldsymbol D}}{{{\boldsymbol s}}^9}$| subsets of data, which were acquired from |${t_6}$| to |${t_9}$|⁠. Data-weighting coefficient |${\omega _d}$| = 1.9 was used throughout this example. The inverted results are shown in Fig. 5. Figure 5. Open in new tabDownload slide Inversion results of the geoelectric models in Fig. 4 given by three resistivity inversion methods: our new 4-D resistivity inversion algorithm, the conventional 4-D resistivity inversion and the independent inversion of each restructured data set. The 4-D inversions produced by our new algorithm are shown at each time step. The multiple black-edged rectangles represent the state of the actual geoelectric models at each time step. Figure 5. Open in new tabDownload slide Inversion results of the geoelectric models in Fig. 4 given by three resistivity inversion methods: our new 4-D resistivity inversion algorithm, the conventional 4-D resistivity inversion and the independent inversion of each restructured data set. The 4-D inversions produced by our new algorithm are shown at each time step. The multiple black-edged rectangles represent the state of the actual geoelectric models at each time step. We begin our analysis of Fig. 5 by examining the results of our new 4-D resistivity inversion algorithm. Row 6 shows the images obtained from the independent inversions. Here, the position and shape of the low-resistivity body are essentially identical to those in the original geoelectric model, |${G_1}$|⁠. As the low-resistivity body undergoes rapid changes at |${t_4}$|⁠, the |${{\bf D}}{{{\bf s}}^5}$| subset of data contains information about the changes between |${G_1}$| and |${G_2}$|⁠. The inverted image shown in the second column of row 2 indicates that our 4-D inversion algorithm is able to respond to these changes. The results in row 3 were obtained after the acquisition of |${{\boldsymbol D}}{{{\boldsymbol s}}^6}$|⁠. In this step, the weight of the time regularization in the 4-D inversion equation was defined by the ratios between the inversion results at three time points. Two of these inversions produced the results shown in row 2 (which were obtained using our new algorithm), while the third produced the independent inversion shown in row 3. The third column in row 3 shows that the inverted image becomes very similar to |${G_2}$| upon the addition of new information about the changes between |${G_1}$| and |${G_2}$| (i.e. the substitution of |${{\boldsymbol D}}{{{\boldsymbol s}}^6}$| into the 4-D inversion equation). Information pertaining to the changes between |${G_2}$| and |${G_3}$| was acquired at |${t_6}$| and |${t_7}$|⁠. The result in the fourth column of row 4, which closely resembles |${G_3}$|⁠, is surprising, as it was obtained using only one-quarter of a complete data set. In other words, by increasing the weights of new data, our new 4-D inversion algorithm was able to produce a result that closely approximates the model at the current point in time. The fifth row shows the final results of our new 4-D algorithm, illustrating the images acquired at each time step of the monitoring process; these images are useful for demonstrating the dynamic evolution of the low-resistivity body during its migration. In the next step, our new 4-D inversion algorithm, the independent inversion, and the conventional 4-D inversion algorithm were compared. In the inversion at |${t_5}$|⁠, neither the independent inversion nor the conventional 4-D inversion were able to capture the rapid changes occurring in the model. Since three-quarters of the data used at this time were obtained from forward modelling, the resulting image resembles |${G_1}$|⁠. By increasing the data weight of |${{\boldsymbol D}}{{{\boldsymbol s}}^5}$|⁠, our new algorithm was able to capture these rapid changes and provide a timely illustration of the current resistivity distribution. At |${t_7}$|⁠, the images produced by the independent inversion and the conventional 4-D inversion algorithm only vaguely resemble |${G_3}$|⁠. By comparison, our new algorithm produces a far more accurate description of the |${G_3}$| model in terms of the spatial position and shape of the low-resistivity body in the inverted image. 4 PHYSICAL MODEL TEST 4.1 Model design A physical model test was designed for fast solute transport monitoring and imaging. Data was collected using cross-hole ERT and processed using our new 4-D resistivity inversion algorithm. The aim of this experiment was to verify the viability and applicability of our 4-D resistivity inversion scheme. The overall design of the physical model is shown in Fig. 6, and photos taken during this experiment are shown in Figs 6(b) and (c). This model test was conducted at a farm in Zhangqiu, Shandong Province, China. A fillable area was prepared by excavating a 2.0 m × 1.5 m × 1.5 m (length × width × depth) area; one of the surfaces was lined with a wooden baffle, while the three remaining surfaces consisted of exposed soil directly connected to surrounding soil. This configuration approximates the infinite boundary requirement and helps to minimize boundary effects in the resistivity experiment. The fillable area was filled with a mixture of standardized sand particles (particle sizes ranging between 0.5 and 2.0 mm) and farm soils at a volumetric ratio of 1:2. We use sand to simulate the fast solution transport for its good permeability and use clay for the stability of the model when we finally cut the model. After the area was filled and tamped, the resistivity of the fill was found to range between 300 and 500 Ω m. During filling with the similitude materials, two probe holes (P1 and P2) and three water injection holes (P3, P4 and P5) were installed in the middle of the model. The water injection holes were spaced 0.2 m apart and the pipes were located in the section formed by the P1 and P2 detection holes. The water injection outlets were located at depths of 0.9 m (P3), 0.6 m (P4) and 0.3 m (P5) inside the model fill medium, as shown in Fig. 6(a). The two detection holes, P1 and P2, were spaced 1.0 m apart, and 14 copper electrodes were installed for cross-hole ERT in these holes with an interelectrode spacing of 0.1 m. One end of each electrode was connected to an external ERT instrument via cables drawn through the detection holes. Figure 6. Open in new tabDownload slide A fast solute transport monitoring experimental model. (a) The overall design of the experimental model in which P1 and P2 are the detection holes, each carrying 14 copper electrodes, and P3, P4 and P5 are used to simulate three water injection pipes with different outlet depths. (b) Demonstration of the testing process. The experimental area are approximately 2.0 m × 1.5 m × 1.5 m (length × width × depth) and one of the model surfaces is lined with a wooden baffle while the three remaining surfaces directly connect to surrounding soil. (c) The final experimental model was filled with a volume ratio of 1:2 of sand and soil mixture, and water was injected into P3, P4 and P5 for the monitoring experiments. Figure 6. Open in new tabDownload slide A fast solute transport monitoring experimental model. (a) The overall design of the experimental model in which P1 and P2 are the detection holes, each carrying 14 copper electrodes, and P3, P4 and P5 are used to simulate three water injection pipes with different outlet depths. (b) Demonstration of the testing process. The experimental area are approximately 2.0 m × 1.5 m × 1.5 m (length × width × depth) and one of the model surfaces is lined with a wooden baffle while the three remaining surfaces directly connect to surrounding soil. (c) The final experimental model was filled with a volume ratio of 1:2 of sand and soil mixture, and water was injected into P3, P4 and P5 for the monitoring experiments. We used a measurement array known as the Bipole_Bipole (AM-BN). A super high-density resistivity instrument of the type FlashRES-UNIVERSAL was used, and the on/off time is 0.5 and 0.3 s separately during the test. The acquisition of each complete data set contained 960 independent data and took approximately 12 min. The data directly obtained by the instrument acquisition was resistance, that is, the ratio of the potential difference |$\Delta {{\boldsymbol U}}$| between two potential electrodes of M and N and the current I, when the current electrode A was positively charged and the current electrode B was negatively charged. The water injection and monitoring experiment was divided into three stages, and the entirety of the experiment was performed over two days. Here, in the experimental procedures, water injection monitoring at P3 was an independent process while water injection monitoring at P4 and P5 was a temporally continuous process. The beginning of the injection at P5 corresponds to the appearance of a water body at a new position within the model space, as well as the initiation of solute dispersion processes. Constant-pressure water injection encourages relatively stable solute dispersion and migration since the spatial locations of the water injection points are at fixed points. The procedures for each experimental stage are described in Table 1. On the first day of the experiment, background resistivity data were collected prior to water injection. The result of an inversion using these data was used as the initial model in the next inversion. Brine (a mixture of salt and tap-water, with a concentration of approximately 50 g L−1 and resistivity of about 0.32 Ω m) was then injected continuously into P3 at a constant pressure (by maintaining the height of the water column), and simultaneous data acquisition was performed. On the second day of the experiment, background resistivity data were again acquired prior to the monitoring experiment, and an inversion of this data was used as the initial model in the next inversion. Brine was then injected at a constant pressure into P4 and P5 while the fill medium was monitored. Water injection was terminated after a final injection of 150 mL. Rhodamine B (1.0 g L−1) was used in the brine as a tracer during this experiment; the differences between the dyed region and the inverted images were then examined by removing sections from the experimental soil body after the end of the experiment. Table 1. Duration of and water injection volume used in each experimental stage. Experimental stage . Duration . Water injection volume (mL) . Description of experiment . Injection into P3 (Step 1) 14:40–14:52 0 No water injection 14:52–15:05 460 Constant-pressure water injection (The water column is kept at a constant height) 15:05–15:17 270 Constant-pressure water injection Injection into P4 (Step 2) 9:48–10:00 0 No water injection 10:00–10:13 410 Constant-pressure water injection 10:13–10:19 190 Constant-pressure water injection Injection into P5 (Step 3) 10:19–10:25 180 Constant-pressure water injection 10:25–10:38 230 Constant-pressure water injection 10:38–10:50 150 Non-constant-pressure water injection (Injection terminated after 150 ml) Experimental stage . Duration . Water injection volume (mL) . Description of experiment . Injection into P3 (Step 1) 14:40–14:52 0 No water injection 14:52–15:05 460 Constant-pressure water injection (The water column is kept at a constant height) 15:05–15:17 270 Constant-pressure water injection Injection into P4 (Step 2) 9:48–10:00 0 No water injection 10:00–10:13 410 Constant-pressure water injection 10:13–10:19 190 Constant-pressure water injection Injection into P5 (Step 3) 10:19–10:25 180 Constant-pressure water injection 10:25–10:38 230 Constant-pressure water injection 10:38–10:50 150 Non-constant-pressure water injection (Injection terminated after 150 ml) Open in new tab Table 1. Duration of and water injection volume used in each experimental stage. Experimental stage . Duration . Water injection volume (mL) . Description of experiment . Injection into P3 (Step 1) 14:40–14:52 0 No water injection 14:52–15:05 460 Constant-pressure water injection (The water column is kept at a constant height) 15:05–15:17 270 Constant-pressure water injection Injection into P4 (Step 2) 9:48–10:00 0 No water injection 10:00–10:13 410 Constant-pressure water injection 10:13–10:19 190 Constant-pressure water injection Injection into P5 (Step 3) 10:19–10:25 180 Constant-pressure water injection 10:25–10:38 230 Constant-pressure water injection 10:38–10:50 150 Non-constant-pressure water injection (Injection terminated after 150 ml) Experimental stage . Duration . Water injection volume (mL) . Description of experiment . Injection into P3 (Step 1) 14:40–14:52 0 No water injection 14:52–15:05 460 Constant-pressure water injection (The water column is kept at a constant height) 15:05–15:17 270 Constant-pressure water injection Injection into P4 (Step 2) 9:48–10:00 0 No water injection 10:00–10:13 410 Constant-pressure water injection 10:13–10:19 190 Constant-pressure water injection Injection into P5 (Step 3) 10:19–10:25 180 Constant-pressure water injection 10:25–10:38 230 Constant-pressure water injection 10:38–10:50 150 Non-constant-pressure water injection (Injection terminated after 150 ml) Open in new tab 4.2 Parameter configuration In this experiment, each complete set of observation data was divided into two parts for restructuring and inversion. As shown in Fig. 7, current A and potential electrodes M were placed in the first borehole and current B and potential electrodes N were placed in the second borehole. A and B were defined at an identical level during the acquisition process. For the first part, only odd current electrodes were used, whereas for the second part, we used only even current electrodes. Figure 7. Open in new tabDownload slide A complete data set was split into two parts in the physical model test. For the first part, only odd current electrodes were used, whereas the second part used only even current electrodes. Figure 7. Open in new tabDownload slide A complete data set was split into two parts in the physical model test. For the first part, only odd current electrodes were used, whereas the second part used only even current electrodes. Then, the restructured data configuration was recombined in the second part of the last complete data set with the newly acquired data (i.e. the first part of the next complete data set), as shown in Fig. 8. We observed that the restructured data configuration also includes two parts, that is, even and odd parts. Both the restructured and original complete data sets have identical configurations. As their data volume and structure are identical, only their data sequences are different. Therefore, both data sets have identical sensitivity. Figure 8. Open in new tabDownload slide Illustration of the data set restructuring scheme in the physical model test. Figure 8. Open in new tabDownload slide Illustration of the data set restructuring scheme in the physical model test. Throughout this experiment, data was transmitted in real-time to a server for inversion calculations. The inversion modelling was mainly performed in the 1 m (x) × 1.4 m (z) region between the two detection holes. The smallest grids in the finite element analysis were quadrilateral and 0.05 m × 0.05 m in size. Data-weighting coefficient |${\omega _d}$| = 2.2 was used throughout this model test. The GPU parallel inversion program is performed before the data collection in the next part is finished. The final results produced by our new 4-D inversion method are shown in Fig. 9. Figure 9. Open in new tabDownload slide Inverted images obtained for each experimental stage. Rows 1–3 correspond to the inversion results after water injection into P1–P3, respectively. The black solid ellipse with the dashed line marks the position of the water injection outlet. Figure 9. Open in new tabDownload slide Inverted images obtained for each experimental stage. Rows 1–3 correspond to the inversion results after water injection into P1–P3, respectively. The black solid ellipse with the dashed line marks the position of the water injection outlet. 4.3 Experimental results and analysis Fig. 9 delineates the continuously migrating and changing low-resistivity body. To analyse the fit of the data in the physical model test, we compared the inversion results of the forward modelling with the measured data, as shown in Fig. 10. The forward modelling results for a1–a5, b1–b5 and c1–c5, as shown in Fig. 10, correspond to the inversion results in Fig. 9. It can be seen that the data fit well according to the curve comparison in the figures. Figure 10. Open in new tabDownload slide Open in new tabDownload slide Comparison of the forward modelling results (model from Fig. 7 in the original script) and the corresponding measured data. Figure 10. Open in new tabDownload slide Open in new tabDownload slide Comparison of the forward modelling results (model from Fig. 7 in the original script) and the corresponding measured data. However, owing to factors, such as uneven soil and sand mixture and spatially varying degrees of compaction, the anomalies that did not change over time appeared in the background and subsequent inversion results. Therefore, ERT images are provided in terms of percentage changes in resistivity to eliminate anomalies and more clearly reflect the position, shape and size of the water body in the sandy medium throughout water migration and dispersion. The percentage change model is defined as (Bellmunt et al.2012) $$\begin{eqnarray*} \delta {{\boldsymbol{M}}^{{{\prime}}}} = \frac{{{{\skew{4}\hat{\boldsymbol M}}} - {{{\boldsymbol M}}^{\mathrm{ ref}}}}}{{{{{\boldsymbol M}}^{\mathrm{ ref}}}}} \times 100\,\mathrm{ per\, cent} , \end{eqnarray*}$$ (11) where |${{\skew{4}\hat{\boldsymbol M}}}$| is the inversion model parameter at the current point in time and |${{{\boldsymbol M}}^{\mathrm{ ref}}}$| is the model parameter of the reference model. For each water injection stage, the first inverted result in Fig. 9 was used as the reference model (a1, b1 and c1 in Fig. 9). The percentage changes between the reference model resistivities and subsequent inversion resistivities were then calculated to produce the results shown in Fig. 11. Figs 9 and 11 illustrate the changes in the resistivity distribution of the underground medium during the injection of water into the medium via P3, P4 and P5. Figure 11. Open in new tabDownload slide Percent change in inverted results for each experimental stage. These images were obtained by calculating the difference between the first image in each row of Fig. 9 and other images in the same row, and then calculating the per cent quotient of these differences. The black solid ellipse with the dashed line marks the position of the water injection outlets. Figure 11. Open in new tabDownload slide Percent change in inverted results for each experimental stage. These images were obtained by calculating the difference between the first image in each row of Fig. 9 and other images in the same row, and then calculating the per cent quotient of these differences. The black solid ellipse with the dashed line marks the position of the water injection outlets. As shown in Fig. 9, independent analyses were conducted on the inverted images produced via continuous monitoring during the migration and dispersion of the low-resistivity body in each experimental stage; this analysis shows that the dynamic changes in the low-resistivity body consist primarily of horizontal migration/dispersion and vertical (downwards) infiltration. The percentage resistivity change plots in Fig. 11 characterize the independent motions of the water body during each stage; these images represent the percentage change in the difference between each inverted image and the reference image (the first image in each stage). Independent analyses of the percentage change plots for each experimental stage show that low-resistivity body dispersion always began near the water injection points (i.e. near the water injection outlets in the filling medium). This demonstrates that our new 4-D resistivity inversion algorithm is able to accurately locate the position at which changes in medium resistivity first begin to occur. Fig. 11 shows that the changes in the shape and size of the low-resistivity body were gradual and smooth during each stage, which is consistent with the actual water body movement. It shows that our new algorithm was able to rapidly capture the low-resistivity anomaly changes in the underground. Fig. 12(a) shows the photograph of a vertical section cut along the probed profile. The red dashed lines in the photos delineate the regions dyed with Rhodamine B, while the black dashed lines delineate the areas around the dyed regions that contained relatively high levels of moisture. The dye dispersion range was limited, as the filling medium-sand mixture also adsorbed the dye. Regions with higher levels of moisture distinct from the surrounding sandy medium can be observed at the boundaries of the dyed regions; these regions are marked by black dashed lines in Fig. 12. It is our opinion that the dyed regions reflect (or slightly underestimate) the actual extent of water migration during the water injection monitoring experiment. The water body continued to move downwards and diffuse into the surroundings because of gravity and the heterogeneity of the soil when water injection was paused or stopped and throughout a short stabilization period. However, the dye began to dissociate from the water body during this time, preventing the dye from dispersing alongside the water body. Figure 12. Open in new tabDownload slide Photos of the actual sections cut into the model after the end of the experiment. (a) Vertical section 1, that is, the profile probed by the resistivity experiment. (b) Composite inverted image at the final time point in each experimental stage (a5_a1, b5_b1 and c5_c1 in Fig. 11). Figure 12. Open in new tabDownload slide Photos of the actual sections cut into the model after the end of the experiment. (a) Vertical section 1, that is, the profile probed by the resistivity experiment. (b) Composite inverted image at the final time point in each experimental stage (a5_a1, b5_b1 and c5_c1 in Fig. 11). Comparing the dyed regions marked by red dashed lines in Fig. 12(a) and the images inverted at the final time point in Fig. 12(b), it is clear that the extent of dispersion and outer contour shape of the water body (as indicated by the dyed regions around the outlet of each water injection pipe) are consistent with the shape shown in the inverted image. The solute dispersion from each water injection outlet did not overlap in reality; however, the inverted images show merging of zones into a single contiguous area. This is caused by the volume effects and regularization constraints intrinsic to the ERT method, which lead to overestimation of the size of the low-resistivity bodies in the inverted images and errors in the identification and location of body boundaries. 5 DISCUSSION Unlike most studies that focus on the monitoring of relatively stable underground fluid motion over long periods of time, this study explores the recognition and imaging of rapid changes in medium resistivity that occur over short periods of time: that is, rapid changes that occur in a medium within a single acquisition cycle. In these cases, the data set recorded within a single acquisition can be quite complex, as it simultaneously contains information about static geological structures and dynamic fluid motions. The corresponding spatial-temporal resistivity model is, consequently, highly complex. In the conventional 4-D resistivity inversion method, several continuous and complete data sets are inverted simultaneously, which makes it difficult to detect rapid resistivity changes occurring in the medium within the given acquisition cycle. Owing to aliasing (i.e. not containing enough data), the inversion results of resistivity monitoring usually encountered the problem of temporal smear (Rucker 2014). Currently, the best solution to this problem is to use rapid, multiple-channel resistivity systems. In this study, we used different data-weighting and adaptive temporal regularization techniques as important measures for solving the aforementioned problem. Specifically, multiple reference time points were selected to divide complete data sets from each acquisition cycle into multiple segments; these segments were then restructured and dynamically substituted into the improved 4-D inversion equation (which applied increased weighting factors for new data and actively adjusted the time-regularization coefficients) for inversion computations. By applying these strategies, the artefacts caused by the data of previous time could be suppressed to some extent, helping to model imaging at the current time. This facilitates the inversion of images over shorter time intervals and enhances the responsiveness of the 4-D inversions, allowing for the monitoring of rapid changes in underground medium resistivity. In this study, we performed 4-D stitching and achieved a series of good inversion models. In contrast, in the 3-D spatial tiling by Rucker et al. (2009), the results showed that it was difficult to stitch the individual models together even when there was a large overlap. We think the differences between our study and that by Rucker et al. (2009) may be as follows. (1) Achieving good inversion results is more difficult when using spatial tiling than when using stitching in time. For 4-D monitoring, the monitored spatial position is usually fixed, that is, the underground structure model does not change much in a certain time or shows a stable continuous change. Therefore, data stitching in time is dominant in 4-D stitching. Generally, the overall change in data during continuous monitoring is usually small, and thus it is easier to obtain a series of good inversion models. Model stitching is dominant for 2-D or 3-D stitching. In contrast, spatial tiling usually corresponds to the model of different spatial regions, and the change of the underground structure is relatively large. Moreover, the acquired data also change greatly. The difficultly in stitching individual models together could be attributed to a few reasons. First, as the model can only be set to a finite grid and cannot accurately simulate the subsurface resistivity structure, the error generated by modelling from different regions of individual inversion will be amplified. Second, the surrounding geology of the detection and inversion regions is usually not homogeneous, and thus anomalous bodies near different regions can have different effects on the results. Moreover, the site noise in different regions is not consistent and could have different impacts on the inversion results. Thus, it is more difficult to achieve good inversion results using spatial tiling than temporal stitching at a fixed region. (2) With the significant increase in computing power and memory, multiple models at different times can be inverted simultaneously, further simplifying 4-D model stitching. In the proposed algorithm, dynamic data weighting and temporal regularization were introduced, which are beneficial to improve the stable inversion of continuous changes between adjacent models. Note that the number of restructured data sets is proportional to the number of data set splits performed, that is, the number of inversion models for each reference time point can be increased by increasing the number of data set splits. However, this increases the number of inversion computations and the corresponding time required to complete the 4-D inversion. Actual monitoring and data processing/interpretation processes are limited by practical time efficiency requirements and the number of reference time points should be selected accordingly (i.e. the number of reference time points should be selected so as to fit the computational time available for data acquisition and inversion). We performed detailed statistics on the time consumption of inversion of the numerical experiments and physical model test before and after GPU parallelization separately, and obtained results, as shown in Table 2. Among them, N1–N5 represent the five stages |${t_4}\sim {t_{\rm{8}}}$| of the data-weighting numerical simulation experiments. The time in Table 2 is obtained by averaging the time of six iterations in the total inversion process; this is also the average time in terms of four different data-weighting coefficients. b1–b5 and c1–c5 correspond to the 10 stages in the physical model test, and the time is obtained by averaging the time of six iterations. The data comparison in the table clearly shows that the computational efficiency of the inversion solution is increased by nearly 20 times after GPU parallel acceleration. Note that the number of model grids in the table refers to the 4-D model grids, implying that as the monitoring process continues, the corresponding model is added at the new moment, and the number of model grids increases. Table 2. Statistics of time consumed by the inversion process before and after GPU parallelization. N1–N5 represent the five stages (⁠|${t_4}\sim {t_{\rm{8}}}$|⁠) of the data-weighting numerical simulation experiments; b1–b5 and c1–c5 correspond to the 10 stages in the physical model test. Trial . Number of model grids . Number of electrodes . Number of independent data . The average time of each iteration before GPU parallelization (s) . The average time of each iteration after GPU parallelization (s) . Acceleration factor . N1 448 56 5625 839.7 41.9 20.0 N2 896 56 7031 5399.5 273.8 19.7 N3 1344 56 8438 8146.0 397.3 20.5 N4 1792 56 9844 8704.5 438.7 19.8 N5 2240 56 11 250 9719.4 483 .9 20.1 b1 560 28 960 185.9 8.5 21.9 b2 1120 28 1440 404.4 18.1 22.3 b3 1680 28 1920 595.8 25.4 23.5 b4 2240 28 2400 587 34.3 17.1 b5 2800 28 2880 821.9 42.6 19.3 c1 3360 28 3360 941.6 51.2 18.4 c2 3920 28 3840 1100.1 60.8 18.1 c3 4480 28 4320 1212.5 68.5 17.7 c4 5040 28 4800 1806.7 77.1 23.4 c5 5600 28 5280 1925.5 85.6 22.5 Trial . Number of model grids . Number of electrodes . Number of independent data . The average time of each iteration before GPU parallelization (s) . The average time of each iteration after GPU parallelization (s) . Acceleration factor . N1 448 56 5625 839.7 41.9 20.0 N2 896 56 7031 5399.5 273.8 19.7 N3 1344 56 8438 8146.0 397.3 20.5 N4 1792 56 9844 8704.5 438.7 19.8 N5 2240 56 11 250 9719.4 483 .9 20.1 b1 560 28 960 185.9 8.5 21.9 b2 1120 28 1440 404.4 18.1 22.3 b3 1680 28 1920 595.8 25.4 23.5 b4 2240 28 2400 587 34.3 17.1 b5 2800 28 2880 821.9 42.6 19.3 c1 3360 28 3360 941.6 51.2 18.4 c2 3920 28 3840 1100.1 60.8 18.1 c3 4480 28 4320 1212.5 68.5 17.7 c4 5040 28 4800 1806.7 77.1 23.4 c5 5600 28 5280 1925.5 85.6 22.5 Open in new tab Table 2. Statistics of time consumed by the inversion process before and after GPU parallelization. N1–N5 represent the five stages (⁠|${t_4}\sim {t_{\rm{8}}}$|⁠) of the data-weighting numerical simulation experiments; b1–b5 and c1–c5 correspond to the 10 stages in the physical model test. Trial . Number of model grids . Number of electrodes . Number of independent data . The average time of each iteration before GPU parallelization (s) . The average time of each iteration after GPU parallelization (s) . Acceleration factor . N1 448 56 5625 839.7 41.9 20.0 N2 896 56 7031 5399.5 273.8 19.7 N3 1344 56 8438 8146.0 397.3 20.5 N4 1792 56 9844 8704.5 438.7 19.8 N5 2240 56 11 250 9719.4 483 .9 20.1 b1 560 28 960 185.9 8.5 21.9 b2 1120 28 1440 404.4 18.1 22.3 b3 1680 28 1920 595.8 25.4 23.5 b4 2240 28 2400 587 34.3 17.1 b5 2800 28 2880 821.9 42.6 19.3 c1 3360 28 3360 941.6 51.2 18.4 c2 3920 28 3840 1100.1 60.8 18.1 c3 4480 28 4320 1212.5 68.5 17.7 c4 5040 28 4800 1806.7 77.1 23.4 c5 5600 28 5280 1925.5 85.6 22.5 Trial . Number of model grids . Number of electrodes . Number of independent data . The average time of each iteration before GPU parallelization (s) . The average time of each iteration after GPU parallelization (s) . Acceleration factor . N1 448 56 5625 839.7 41.9 20.0 N2 896 56 7031 5399.5 273.8 19.7 N3 1344 56 8438 8146.0 397.3 20.5 N4 1792 56 9844 8704.5 438.7 19.8 N5 2240 56 11 250 9719.4 483 .9 20.1 b1 560 28 960 185.9 8.5 21.9 b2 1120 28 1440 404.4 18.1 22.3 b3 1680 28 1920 595.8 25.4 23.5 b4 2240 28 2400 587 34.3 17.1 b5 2800 28 2880 821.9 42.6 19.3 c1 3360 28 3360 941.6 51.2 18.4 c2 3920 28 3840 1100.1 60.8 18.1 c3 4480 28 4320 1212.5 68.5 17.7 c4 5040 28 4800 1806.7 77.1 23.4 c5 5600 28 5280 1925.5 85.6 22.5 Open in new tab In this study, we used a measurement array known as the Bipole_Bipole (AM-BN). However, this is not the optimal choice. Wilkinson et al. (2012) proposed an algorithm to optimize survey design, which can effectively reduce polarization effect interference on the data. If we use this method, we may reduce the effect of polarization on the 4-D inversion results without pre-processing the data. 6 CONCLUSIONS In this study, we propose an improved 4-D resistivity inversion method for monitoring rapid changes in medium resistivity based on the conventional 4-D inversion algorithm. Through the use of GPU-parallelized inversion and the restructuring of new measurement data, our algorithm is able to perform data extraction and inversion over short time intervals with a 20-fold improvement in computational efficiency and without compromising image resolution. The imaging of medium changes can be improved by increasing the weight of newly introduced measurement data in the 4-D inversion equation, as shown by numerical simulations. Numerical inversion experiments on geoelectric models of rapid changes in medium resistivity show that the active adjustment of time-regularization coefficients (using the predicted range of changes in model resistivity) can be used to highlight changes in medium resistivity and enhance the stability of the inversion process. The viability and effectiveness of our new algorithm was also validated via a field experiment in which solute dispersion was successfully monitored using our new method. In summary, our improved method for the 4-D inversion of resistivity monitoring data is capable of accurately and quickly imaging rapid changes in underground medium resistivity. NOTATIONS The following symbols are used in this paper: |${\boldsymbol \Phi }$| 4-D inversion objective function |${{\boldsymbol \varphi }_d}$| Data misfit term |${{\boldsymbol \varphi }_m}$| Space-regularization term |${{\boldsymbol \varphi }_t}$| Time-regularization term |$\lambda $| Weight of the space regularization term |$\mu $| Weight of the time-domain regularization weighting coefficient C Space-regularization matrix |${C_T}$| Time-regularization matrix |${W_d}$| Data-weighting matrix |${\omega _d}$| Weight of the newest subset of data in the inversion equation |${{\boldsymbol e}}$| Misfit vector between actual measured data and theoretical data |${{{\boldsymbol D}}^i}(i = 1,2,...)$| The ith collected complete data set |${\boldsymbol{D}}{{\boldsymbol{s}}^k}(k = 1,2,...)$| The kth segmented subset of data |${\boldsymbol{D}}{{\boldsymbol{r}}^j}(j = 1,2,...)$| The jth complete data set from the subset of data restructuring |${N_S}$| Number of elements in subset of data |${\boldsymbol{D}}{{\boldsymbol{s}}^k}$| |${\bf G}$| Forward modelling |${{{\boldsymbol M}}^j}(j = 1,2,...)$| A group of model parameters |${\boldsymbol{\Delta M}}$| Perturbation vector of the model parameters during the iteration of |${\boldsymbol{M}}$| |${{\skew{4}\hat{\boldsymbol M}}}$| Model parameters after inversion |${J^j}(j = 1,2,...)$| Jacobian matrix for each model J Extended matrix of |${J^j}$| in the time domain |${A_J}$| Coefficient matrix of the JPCG algorithm |${\boldsymbol{b}}$| Rightmost term of the JPCG algorithm |${\boldsymbol{x}}$| Vector to be solved in the JPCG algorithm |${\boldsymbol \Phi }$| 4-D inversion objective function |${{\boldsymbol \varphi }_d}$| Data misfit term |${{\boldsymbol \varphi }_m}$| Space-regularization term |${{\boldsymbol \varphi }_t}$| Time-regularization term |$\lambda $| Weight of the space regularization term |$\mu $| Weight of the time-domain regularization weighting coefficient C Space-regularization matrix |${C_T}$| Time-regularization matrix |${W_d}$| Data-weighting matrix |${\omega _d}$| Weight of the newest subset of data in the inversion equation |${{\boldsymbol e}}$| Misfit vector between actual measured data and theoretical data |${{{\boldsymbol D}}^i}(i = 1,2,...)$| The ith collected complete data set |${\boldsymbol{D}}{{\boldsymbol{s}}^k}(k = 1,2,...)$| The kth segmented subset of data |${\boldsymbol{D}}{{\boldsymbol{r}}^j}(j = 1,2,...)$| The jth complete data set from the subset of data restructuring |${N_S}$| Number of elements in subset of data |${\boldsymbol{D}}{{\boldsymbol{s}}^k}$| |${\bf G}$| Forward modelling |${{{\boldsymbol M}}^j}(j = 1,2,...)$| A group of model parameters |${\boldsymbol{\Delta M}}$| Perturbation vector of the model parameters during the iteration of |${\boldsymbol{M}}$| |${{\skew{4}\hat{\boldsymbol M}}}$| Model parameters after inversion |${J^j}(j = 1,2,...)$| Jacobian matrix for each model J Extended matrix of |${J^j}$| in the time domain |${A_J}$| Coefficient matrix of the JPCG algorithm |${\boldsymbol{b}}$| Rightmost term of the JPCG algorithm |${\boldsymbol{x}}$| Vector to be solved in the JPCG algorithm Open in new tab |${\boldsymbol \Phi }$| 4-D inversion objective function |${{\boldsymbol \varphi }_d}$| Data misfit term |${{\boldsymbol \varphi }_m}$| Space-regularization term |${{\boldsymbol \varphi }_t}$| Time-regularization term |$\lambda $| Weight of the space regularization term |$\mu $| Weight of the time-domain regularization weighting coefficient C Space-regularization matrix |${C_T}$| Time-regularization matrix |${W_d}$| Data-weighting matrix |${\omega _d}$| Weight of the newest subset of data in the inversion equation |${{\boldsymbol e}}$| Misfit vector between actual measured data and theoretical data |${{{\boldsymbol D}}^i}(i = 1,2,...)$| The ith collected complete data set |${\boldsymbol{D}}{{\boldsymbol{s}}^k}(k = 1,2,...)$| The kth segmented subset of data |${\boldsymbol{D}}{{\boldsymbol{r}}^j}(j = 1,2,...)$| The jth complete data set from the subset of data restructuring |${N_S}$| Number of elements in subset of data |${\boldsymbol{D}}{{\boldsymbol{s}}^k}$| |${\bf G}$| Forward modelling |${{{\boldsymbol M}}^j}(j = 1,2,...)$| A group of model parameters |${\boldsymbol{\Delta M}}$| Perturbation vector of the model parameters during the iteration of |${\boldsymbol{M}}$| |${{\skew{4}\hat{\boldsymbol M}}}$| Model parameters after inversion |${J^j}(j = 1,2,...)$| Jacobian matrix for each model J Extended matrix of |${J^j}$| in the time domain |${A_J}$| Coefficient matrix of the JPCG algorithm |${\boldsymbol{b}}$| Rightmost term of the JPCG algorithm |${\boldsymbol{x}}$| Vector to be solved in the JPCG algorithm |${\boldsymbol \Phi }$| 4-D inversion objective function |${{\boldsymbol \varphi }_d}$| Data misfit term |${{\boldsymbol \varphi }_m}$| Space-regularization term |${{\boldsymbol \varphi }_t}$| Time-regularization term |$\lambda $| Weight of the space regularization term |$\mu $| Weight of the time-domain regularization weighting coefficient C Space-regularization matrix |${C_T}$| Time-regularization matrix |${W_d}$| Data-weighting matrix |${\omega _d}$| Weight of the newest subset of data in the inversion equation |${{\boldsymbol e}}$| Misfit vector between actual measured data and theoretical data |${{{\boldsymbol D}}^i}(i = 1,2,...)$| The ith collected complete data set |${\boldsymbol{D}}{{\boldsymbol{s}}^k}(k = 1,2,...)$| The kth segmented subset of data |${\boldsymbol{D}}{{\boldsymbol{r}}^j}(j = 1,2,...)$| The jth complete data set from the subset of data restructuring |${N_S}$| Number of elements in subset of data |${\boldsymbol{D}}{{\boldsymbol{s}}^k}$| |${\bf G}$| Forward modelling |${{{\boldsymbol M}}^j}(j = 1,2,...)$| A group of model parameters |${\boldsymbol{\Delta M}}$| Perturbation vector of the model parameters during the iteration of |${\boldsymbol{M}}$| |${{\skew{4}\hat{\boldsymbol M}}}$| Model parameters after inversion |${J^j}(j = 1,2,...)$| Jacobian matrix for each model J Extended matrix of |${J^j}$| in the time domain |${A_J}$| Coefficient matrix of the JPCG algorithm |${\boldsymbol{b}}$| Rightmost term of the JPCG algorithm |${\boldsymbol{x}}$| Vector to be solved in the JPCG algorithm Open in new tab ACKNOWLEDGEMENTS This work was supported by the National Key Research and Development Plan (Nos. 2016YFC0401805, 2016YFC0401801 and 2016YFC0801604), the National Program on Key Basic Research Project of China (Nos. 2014CB046901 and 2015CB058101), the National Key Scientific Instrument and Equipment Development Project (No. 51327802), the National Natural Science Foundation of China (Nos. 51739007, 51479104 and 41502279), the Royal Academy of Engineering under the UK-China Industry Academia Partnership Programme (UK-CIAPP\314), the Shandong Province Key Research and Development Plan (No. 2016ZDJS02A01), the Henan Province Major Science and Technology Special Projects program (No. 161100211100), and the Shandong University Fundamental Research Fund (No. 2017JC002). The authors gratefully acknowledge the support listed above. We would like to thank reviewers for their constructive comments that have helped to improve this paper. REFERENCES Bellmunt F. , Marcuello A. , Ledo J. , Queralt P. , Falgàs E. , Benjumea B. , Velasco V. , Vázquez-Suñé E. , 2012 . Time-lapse Cross-hole electrical resistivity tomography monitoring effects of an urban tunnel , J. Appl. Geophys. , 87 ( 12 ), 60 – 70 . 10.1016/j.jappgeo.2012.09.003 Google Scholar Crossref Search ADS WorldCat Crossref Bellmunt F. , Marcuello A. , Ledo J. , Queralt P. , 2016 . Capability of cross-hole electrical configurations for monitoring rapid plume migration experiments , J. Appl. Geophys. , 124 , 73 – 82 . 10.1016/j.jappgeo.2015.11.010 Google Scholar Crossref Search ADS WorldCat Crossref Case J.S. , 2012 . Inspection of earthen embankment dams using time lapse electrical resistivity tomography, Master's thesis , University of Mississippi . Chen Z. , Meng X. , Guo L. , Liu G. , 2012 . GICUDA: a parallel program for 3D correlation imaging of large scale gravity and gravity gradiometry data on graphics processing units with CUDA , Comput. Geosci. , 46 ( 3 ), 119 – 128 . OpenURL Placeholder Text WorldCat Daily W. , Ramirez A. , Labrecque D. , Nitao J. , 1992 . Electrical resistivity tomography of vadose water movement , Water Resour. Res. , 28 ( 5 ), 1429 – 1442 . 10.1029/91WR03087 Google Scholar Crossref Search ADS WorldCat Crossref Daily W. , Ramirez A. , Binley A. , Labrecque D. , 2004 . Electrical resistance tomography , Leading Edge , 23 ( 5 ), 438 – 442 . Google Scholar Crossref Search ADS WorldCat Doetsch J. , Linde N. , Vogt T. , Binley A. , Green A.G. , 2012 . Imaging and quantifying salt-tracer transport in a riparian groundwater system by means of 3D ERT monitoring , Geophysics , 77 ( 5 ), B207 – B218 . 10.1190/geo2012-0046.1 Google Scholar Crossref Search ADS WorldCat Crossref Hayley K. , Pidlisecky A. , Bentley L.R. , 2011 . Simultaneous time-lapse electrical resistivity inversion , J. Appl. Geophys. , 75 ( 2 ), 401 – 411 . 10.1016/j.jappgeo.2011.06.035 Google Scholar Crossref Search ADS WorldCat Crossref Hermans T. , Oware E. , Caers J. , 2016 . Direct prediction of spatially and temporally varying physical properties from time‐lapse electrical resistance data , Water Resour. Res. , 52 ( 9 ), doi:10.1002/2016WR019126 . 10.1002/2016WR019126 OpenURL Placeholder Text WorldCat Crossref Hilbich C. , Fuss C. , Hauck C. , 2011 . Automated time-lapse ERT for improved process analysis and monitoring of frozen ground , Permafrost Periglacial Process. , 22 ( 4 ), 306 – 319 . Google Scholar Crossref Search ADS WorldCat Jayawickreme D.H. , Dam R.L.V. , Hyndman D.W. , 2010 . Hydrological consequences of land-cover change: quantifying the influence of plants on soil moisture with time-lapse electrical resistivity , Geophysics , 75 ( 75 ), WA43 – WA50 . 10.1190/1.3464760 Google Scholar Crossref Search ADS WorldCat Crossref Karaoulis M.C. , Kim J.H. , Tsourlos P.I. , 2011 . 4D active time constrained resistivity inversion , J. Appl. Geophys. , 73 ( 1 ), 25 – 34 . 10.1016/j.jappgeo.2010.11.002 Google Scholar Crossref Search ADS WorldCat Crossref Karaoulis M.C. , Tsourlos P. , Kim J.H. , Revil A. , 2014 . 4D time-lapse ERT inversion: introducing combined time and space constraints , Near Surf. Geophys. , 11 ( 1 ), 25 – 34 . 10.3997/1873-0604.2013004 Google Scholar Crossref Search ADS WorldCat Crossref Kim J.H. , 2005 . Four dimensional inversion of DC resistivity monitoring data , in Near Surface 2005–11th European Meeting of Environmental and Engineering Geophysics , doi:10.3997/2214-4609-pdb.13.A006 . OpenURL Placeholder Text WorldCat Kim J.H. , Yi M.J. , Park S.G. , Kim J.G. , 2009 . 4-D inversion of DC resistivity monitoring data acquired over a dynamically changing earth model , J. Appl. Geophys. , 68 ( 4 ), 522 – 532 . 10.1016/j.jappgeo.2009.03.002 Google Scholar Crossref Search ADS WorldCat Crossref Kim J.H. , Supper R. , Tsourlos P. , Yi M.J. , 2013 . Four-dimensional inversion of resistivity monitoring data through Lp norm minimizations , Geophys. J. Int. , 195 ( 3 ), 1640 – 1656 . 10.1093/gji/ggt324 Google Scholar Crossref Search ADS WorldCat Crossref Kuras O. et al. . , 2016 . Geoelectrical monitoring of simulated subsurface leakage to support high-hazard nuclear decommissioning at the Sellafield Site, UK , Sci. Total Environ. , 566–567 , 350 – 359 . 10.1016/j.scitotenv.2016.04.212 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Labrecque D.J. , Yang X. , 2001 . Difference inversion of ERT data: a fast inversion method for 3-D in situ monitoring , J. Environ. Eng. Geophys. , 6 ( 2 ), 316 – 321 . Google Scholar Crossref Search ADS WorldCat Liu G.F , Liu Q. , Li B. , Tong X.L , Liu H. , 2009 . GPU/CPU co-processing parallel computation for seismic data processing in oil and gas exploration , J. Prog. Geophys. , 24 ( 5 ), 1671 – 1678 . OpenURL Placeholder Text WorldCat Liu L. , Ding R. , Liu H. , Liu H. , 2015 . 3D hybrid-domain full waveform inversion on GPU , Comput. Geosci. , 83 ( C ), 27 – 36 . Google Scholar Crossref Search ADS WorldCat Loke M.H. , 1999 . Constrained time-lapse resistivity imaging inversion , in Proceedings of 14th EEGS Symp. Application of Geophysics to Engineering and Environmental Problems , Denver, Colorado, USA , EEM7 . OpenURL Placeholder Text WorldCat Loke M.H. , Dahlin T. , Rucker D.F. , 2014 . Smoothness-constrained time-lapse inversion of data from 3D resistivity surveys , Near Surf. Geophys. , 12 ( 1 ), 5 – 24 . 10.3997/1873-0604.2013025 Google Scholar Crossref Search ADS WorldCat Crossref Robinson J. , Buda A. , Collick A. , Shober A. , Ntarlagiannis D. , Slater L. , 2018 . Understanding hydrological connectivity of nutrient transport from agricultural fields to drainage ditches using electrical geophysics , in Proceedings of the 8th Int. Conf. Environmental and Engineering Geophysics , Hangzhou, China . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rucker D.F. , 2014 . Investigating motion blur and temporal aliasing from time-lapse electrical resistivity , J. Appl. Geophys. , 111 , 1 – 13 . 10.1016/j.jappgeo.2014.09.010 Google Scholar Crossref Search ADS WorldCat Crossref Rucker D.F. , Levitt M.T. , Greenwood W.J. , 2009 . Three-dimensional electrical resistivity model of a nuclear waste disposal site , J. Appl. Geophys. , 69 ( 3-4 ), 150 – 164 . 10.1016/j.jappgeo.2009.09.001 Google Scholar Crossref Search ADS WorldCat Crossref Rucker D.F. , Crook N. , Winterton J. , McNeill M. , Baldyga C.A. , Noonan G. , Fink J.B. , 2014 . Real-time electrical monitoring of reagent delivery during a subsurface amendment experiment . Near Surf. Geophys. , 12 ( 1 ), 151 – 163 . 10.3997/1873-0604.2013017 Google Scholar Crossref Search ADS WorldCat Crossref Ruetsch G. , Fatica M. , 2014 . CUDA Fortran for Scientists and Engineers: Best Practices for Efficient CUDA Fortran Programming , Morgan Kaufmann . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Song S.Y. , Kim B. , Nam M.J. , Lim S.K. , 2016 . An analysis of changes in resistivity of general reservoir dams based on 4D inversion of time-lapse resistivity data . Aseg Extended Abstracts , 2016 ( 1 ), 1 . 10.1071/ASEG2016ab214 Google Scholar Crossref Search ADS WorldCat Crossref Tsourlos P. , Ogilvy R. , Meldrum P. , Williams G. , 2003 . Time-lapse monitoring in single boreholes using electrical resistivity tomography , J. Environ. Eng. Geophys. , 8 ( 1 ), 1 – vii . Google Scholar Crossref Search ADS WorldCat Uhlemann S. et al. . , 2017 . Four-dimensional imaging of moisture dynamics during landslide reactivation , J. geophys. Res. , 122 ( 1 ), 398 – 418 . Google Scholar Crossref Search ADS WorldCat Wilkinson P.B. , Chambers J.E. , Meldrum P.I. , Gunn D.A. , Ogilvy R.D. , Kuras O. , 2010 . Predicting the movements of permanently installed electrodes on an active landslide using time‐lapse geoelectrical resistivity data only , Geophys. J. Int. , 183 ( 2 ), 543 – 556 . 10.1111/j.1365-246X.2010.04760.x Google Scholar Crossref Search ADS WorldCat Crossref Wilkinson P.B. , Loke M.H. , Meldrum P.I. , Chambers J.E. , Kuras O. , Gunn D.A. , Ogilvy R.D. , 2012 . Practical aspects of applied optimized survey design for electrical resistivity tomography . Geophys. J. Int. , 189 ( 1 ), 428 – 440 . 10.1111/j.1365-246X.2012.05372.x Google Scholar Crossref Search ADS WorldCat Crossref Wilkinson P.B. , Chambers J.E. , Uhlemann S. , Meldrum P.I. , Smith A. , Dixon N. , Loke M.H. , 2016 . Reconstruction of landslide movements by inversion of 4-D electrical resistivity tomography monitoring data , Geophys. Res. Lett. , 43 , 1166 – 1174 . 10.1002/2015GL067494 Google Scholar Crossref Search ADS WorldCat Crossref Wolfe M. , 2012 . Cuda Fortran Programming Guide and Reference , The Portland Group . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Zhang Y. , Ghodrati A. , Brooks D.H. , 2005 . An analytical comparison of three spatio-temporal regularization methods for dynamic linear inverse problems in a common statistical framework , Inverse Probl. , 21 ( 1 ), 357 – 382 . 10.1088/0266-5611/21/1/022 Google Scholar Crossref Search ADS WorldCat Crossref © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - A rapid four-dimensional resistivity data inversion method using temporal segmentation
JF - Geophysical Journal International
DO - 10.1093/gji/ggaa019
DA - 2020-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/a-rapid-four-dimensional-resistivity-data-inversion-method-using-SzTddaAIge
SP - 586
VL - 221
IS - 1
DP - DeepDyve
ER -