Jointly reconstructing ground motion and resistivity for ERT-based slope stability monitoring

Summary

Electrical resistivity tomography (ERT) is increasingly being used to investigate unstable slopes and monitor the hydrogeological processes within. But movement of electrodes or incorrect placement of electrodes with respect to an assumed model can introduce significant resistivity artefacts into the reconstruction. In this work, we demonstrate a joint resistivity and electrode movement reconstruction algorithm within an iterative Gauss–Newton framework. We apply this to ERT monitoring data from an active slow-moving landslide in the UK. Results show fewer resistivity artefacts and suggest that electrode movement and resistivity can be reconstructed at the same time under certain conditions. A new 2.5-D formulation for the electrode position Jacobian is developed and is shown to give accurate numerical solutions when compared to the adjoint method on 3-D models. On large finite element meshes, the calculation time of the newly developed approach was also proven to be orders of magnitude faster than the 3-D adjoint method and addressed modelling errors in the 2-D perturbation and adjoint electrode position Jacobian.

Hydrogeophysics, Electrical resistivity tomography (ERT), Inverse theory, Tomography

1 INTRODUCTION

Electrical resistivity tomography (ERT) is a geophysical method where current is injected between pairs of electrodes attached to the surface of the earth while the voltage differences are measured between other electrode pairs. From this data, one may reconstruct a resistivity distribution for the near surface earth that best fits the available data. For ERT, relatively small electrode movements or boundary modelling errors can lead to reconstructions with resistivity artefacts so severe that the resulting image is difficult to interpret (Zhou & Dhalin 2003). Similar reconstruction artefacts have been observed for biomedical electrical impedance tomography (EIT) where the fundamental mathematical problem is identical to ERT (Adler et al. 1996, 2011; Boyle & Adler 2011).

Intermittent electrode movement is expected during long-term monitoring of an active landslide site. One could re-survey the electrode locations after each movement, but this would be time consuming and expensive, particularly for remote locations. Changes in landslide topography are frequently accompanied by changes in the resistivity distribution of the slope: changes in soil properties, such as water saturation, affect slope stability and resistivity.

A resistivity reconstruction that does not account for electrode movement, when movements have occurred, will tend to fit the data poorly. Attempts to force fit a resistivity distribution will, in most cases, lead to images with severe resistivity artefacts which, without further information, could be misinterpreted. Simultaneously adjusting electrode position and resistivity can allow for better data fit, reduced resistivity artefacts, and indications of electrode movement. The methods described here separate the confounding factors, electrode movement and resistivity changes, which have long hindered ERT landslide monitoring attempts.

Previous work on electrode movement has largely focused on 2-D electrode position changes in the plane of the electrodes. In 2-D, changes in boundary shape that are conformal do not result in artefacts (Boyle et al. 2012a). Modelling in 3-D is important because it reflects the finite size of the electrodes and the resulting current flow. For 3-D reconstructions, conformal deformations are limited to rotation, scaling and translation which would normally be identified through an appropriate site survey. 3-D resistivity reconstructions can be prohibitive to calculate for a given level of fidelity and numerical convergence. The so-called 2.5-D method (Dey & Morrison 1979) combines a 2-D parametrization of the resistivity and electrode positions with a 3-D model of current flow in the ground. The method can enable significant reductions in the computational requirements for meshing and current density calculations while producing accurate results with respect to equivalent 3-D models.

Electrode movement and resistivity have been reconstructed for biomedical problems where changes both electrode movement (±1 per cent of electrode spacing) and resistivity (±50 per cent) have been relatively mild, enabling Gauss–Newton difference solutions (Soleimani et al. 2006; Jehl et al. 2015), and for large deformations on cylindrical domains with surface-normal deformations (Dardé et al. 2013a,b; Hyvönen et al. 2014) in simulation and on tank data.

For geophysics, the resistivity contrasts are generally strong, varying by orders of magnitude. These strong contrasts combined with large electrode movements, much greater than 1 per cent of electrode spacing, necessitate an iterative solution because the combined effects of the electrode movement give highly nonlinear interactions. Similar approaches are also currently being developed by other researchers (Kim et al. 2014; Wagner et al. 2015; Loke et al. 2015, 2016, 2017). Results on field data have been somewhat mixed and data dependent.

One approach to account for electrode position modelling errors is to alternate between an electrode position update and resistivity update (Wilkinson et al. 2010, 2016). Such an orthogonal greedy optimization approach may be slow to converge, or even fail to converge, because descent over the optimization surface is restricted to movement parallel to the axes (Cormen et al. 1990). Ideally, both electrode movement and resistivity should be reconstructed simultaneously so that updates can descend in the globally optimal direction across the objective function.

In this work, joint electrode position and resistivity reconstruction methods were applied, in 2.5-D, to two data sets (Figs 1 and 2, line#1 and line #5) from an active landslide in the UK. An efficient and accurate method for calculating the electrode position Jacobian on a 2.5-D resistivity model was developed. This work is motivated by developments in Boyle & Adler (2010); Boyle (2010) for 2-D electrode movement, and Boyle et al. (2014); Boyle (2016); Boyle & Adler (2016) where this data set and preliminary results were presented. A 2-D cross-sectional Finite Element Method (FEM) model of the local slope topology was built, with electrodes modelled using the Complete Electrode Model (CEM) (Somersalo et al. 1992; Rücker & Günther 2011). We reconstruct electrode movement and resistivity in an iterative Gauss–Newton algorithm, showing that under certain conditions both can be simultaneously reconstructed.

Figure 1.

Slope failures at Hollin Hill; (a) rotational failures near the top of the slope above line#5, June 2015, (b) earth flow at the toe of a landslide where line#5 runs through mid-slope with electrodes throughout the slipped material, June 2015, (c) satellite image of the hillside (2016), showing four landslide ‘lobes’, five lines of 32 electrodes as of February 2014, and base station location.[(c) Satellite imagery ©2016 DigitalGlobe, Getmapping plc, Infoterra Ltd Bluesky.]

Figure 2.

Electrode locations; (a) electrode locations for five lines of thirty-two electrodes each, line#1 (blue) to line#5 (green) as of February 2014, (b) where electrodes were 10 mm × 170 mm spikes of stainless steel selected for its conductivity, low cost and corrosion resistance.

2 BACKGROUND

Resistivity imaging has been used in geophysical investigations of the behaviour and precursors of landslides and failure surfaces (Perrone et al. 2004; Lapenna et al. 2005; Lebourg et al. 2005; Jongmans & Garambois 2007; Naudetb et al. 2008; Sass et al. 2008). The technique is attractive because resistivity is strongly dependent on water saturation, fracturing, clay content and weathering which are all key factors in slope stability (Piegaria et al. 2009). Slopes may be monitored over time to observe changes in these key parameters using automated systems to collect and analyse data on a daily basis (Kuras et al. 2009; Lebourg et al. 2010; Supper et al. 2014). Difference images may show immediate changes in water saturation (Suzuki & Higashi 2001; Friedel et al. 2006; Jomard et al. 2007) but are limited in their ability to perform long-term monitoring due to background resistivity changes and electrode movements. As illustrated in Fig. 3, small amounts of electrode movement may introduce significant artefacts (Zhou & Dhalin 2003; Oldenborger et al. 2005; Wilkinson et al. 2008). These artefacts may be reduced by accounting for changes in electrode position (Uhlemann et al. 2017).

Figure 3.

Electrode movement artefacts; simulated reconstructions with a conductive and insulating target (rectangular and square outlines, respectively), each with two electrodes (electrode #2 and #12 of 32 electrodes numbered left-to-right at 5 m intervals) having electrode displacements of (a) 0 per cent, (b) 5 per cent and (c) 25 per cent of electrode spacing on a 2-D half-space reconstruction (40 dB Signal-to-Noise Ratio (SNR), λ = 0.01, Laplace regularization, Wenner measurement pattern). Single or well separated electrode location errors introduce characteristic ‘ringing’ artefacts that can overwhelm resistivity-based information. [Reproduced, in part, from Boyle et al. (2017), fig. 2., with author permission.]

An active landslide was identified in North Yorkshire, UK (54°06΄39.2″N, 0°57΄34.9″W) and has been monitored since 2008 (Chambers et al. 2011; Wilkinson et al. 2016). The landslide is a very slow to slow moving composite multiple earth slide-earth flow (Uhlemann et al. 2016b), where a central portion of the slope has moved downhill by up to 3.5 m a year in some instances (Fig. 1). The slope itself exposes four formations: the Dogger Formation (DF), Whitby Mudstone Formation (WMF), Staithes Sandstone Formation (SSF) and Redcar Mudstone Formation (RMF), from top to bottom. The interfaces between sedimentary layers lie horizontally, with a gentle 5° dip to the North, determined through comparison of material interfaces at surrounding exposed slopes in the region (Chambers et al. 2011). The WMF, as the name implies, is a mudstone clay-based rock that is highly weathered and prone to movement during peak water saturation periods at Hollin Hill, occurring annually in the winter through early spring wet-season. The underlying SSF and unweathered RMF are structurally more competent, and landsliding is postulated to occur within the WMF (Uhlemann et al. 2016b). The slope lies at an average angle of 14° over a change in elevation of 40 m. The landslide movement is determined by translational movements upslope of the WMF–SSF interface. These lead to a loss of support of the local slope of the back-scarp, causing rotational failure. The mass accumulation at the WMF–SSF interface, driven by these translational movements, and elevated pore pressures then cause flow movements which form lobes (Uhlemann et al. 2016a,b).

A grid of five rows of thirty-two permanently installed electrodes travelled along with these movements, shifting their positions relative to each other (Fig. 2). An automated impedance imaging survey was executed bi-daily and data were transmitted to a remote office for storage and analysis. In 2008–2009, a middle section of line#1 exhibited a translational failure with movements of up to 1.6 m. In 2013–2014, upper and middle segments of line#5 had rotational and translational failures of up to 3.5 m.

The dipole-dipole measurement protocol used for line#1 and line#5 are visualized in Figs 4 and 5, showing the sequence of quadrupolar measurements with stimulus dipoles in red and measurement dipoles in blue. A single difference measurement is captured as one row of the protocol in the figure. In the adjacent vertical strip chart, the horizontally aligned measurements at the initial time (green) is contrasted with the homogeneous resistivity estimate (red) of what those measurements would be, and the change in measurements from initial to final measurement (blue). The rightmost strip chart shows the estimated error for each measurement as estimated from differences between reciprocal pairs of data for the initial measurements (green) and final measurements (blue). The figure serves to illustrate three challenges with this data set. First, the measurements do not fit, or nearly fit, a homogeneous resistivity model. Second, the change between the initial and final measurements (resistivity and movement) is as much as that caused by the inhomogeneous resistivity distribution in the initial measurements (resistivity only) for some data. Third, the measurement error varies over time so that the initial and final measurements have different associated error estimates. The sequence of measurements for line#5 (Fig. 5) started at the top of the slope and ran to the base, but is equivalent to that shown for line#1 (Fig. 4). The final measurements for line#5 have significantly worse reciprocal errors, but are the best of many data sets examined for a post-movement data set on line#5.

Figure 4.

Dipole–dipole measurement protocol for line#1; March 2008 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements V_a (green) compared to homogeneous resistivity at 32.1 Ωm shown as V_a − V_h (red), and the change from initial to final measurements V_b − V_a (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements.

Figure 5.

Dipole–dipole measurement protocol for line#5; February 2013 and February 2014 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements V_a (green) compared to homogeneous resistivity at 26.1 Ωm shown as V_a − V_h (red), and the change from initial to final measurements V_b − V_a (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements.

3 PRIOR RECONSTRUCTIONS AT HOLLIN HILL

Electrode movement has been previously reconstructed for line#1 data (2008/2009) along with separate resistivity sections (Wilkinson et al. 2010). The algorithm used in that instance achieved estimates of within 0.2 m (4 per cent of electrode spacing) of the electrode’s true positions as measured by Real-Time Kinematic Global Positioning System (RTK GPS).¹ An initial resistivity-only reconstruction with known electrode positions gave a plausible distribution that was in good agreement with available geological evidence: borehole and auger data, evaluation of local geology, aerial LiDAR, differential GPS, lab correlation of representative samples to measured conductivities, piezometric pore pressure measurements, and in-situ rainfall and temperature records (Fig. 6; Chambers et al. 2011; Merritt et al. 2014). Resistivity reconstructions of data collected after movement exhibited artefacts. These artefacts were reduced when reconstructed movements were incorporated. The electrode movement was penalized in the upslope direction. The position Jacobian was estimated based on an analytic half-space model with a homogeneous resistivity assigned to each group of measurements based on electrode separation. The electrode movement was then reconstructed by minimizing

\begin{eqnarray} \arg \min \sqrt{\sum _i |e_i|^2} + \alpha \sum _j |m_j| + \beta \sum _j \theta (m_j)|m_j| \end{eqnarray}

(1)

for weighting factors α = 0.06 m⁻¹, β = 0.32 m⁻¹, Heaviside step function θ, movement m_j, and misfit e_i = r_b − r_a as the difference between predicted r_b and observed r_a apparent resistivity ratios. The apparent resistivity ratio was calculated as the ratio of the analytic half-space models (9) before and after movement

\begin{eqnarray} r = \frac{\rho ^{\prime }}{\rho } \left( \frac{ \frac{1}{AM^{\prime }} - \frac{1}{BM^{\prime }} - \frac{1}{AN^{\prime }} - \frac{1}{BN^{\prime }}}{\frac{1}{AM} - \frac{1}{BM} - \frac{1}{AN} - \frac{1}{BN}} \right) \end{eqnarray}

(2)

for homogeneous resistivity ρ, electrodes spaced AM, BM, AN, BN, and the updated locations and resistivity after movement indicated by primed ΄ variables. Dipole-dipole data for measurements n = 1 were discarded, as they were judged to be too dependent on transverse movements which were not reconstructed.

A similar approach was applied (Wilkinson et al. 2016), to reconstruct 2-D surface xy-movements for the whole electrode grid by allowing for transverse movements through an additional weighting term

\begin{eqnarray} \arg \min &&\sum _i e_i^2 + \alpha \sum _j ||m_j|| + \beta \sum _j \theta (m_j^{(y)})|m_j^{(y)}| \nonumber\\ &&\quad+ \, \gamma \sum _j \theta (m_j^{(x)})|m_j^{(x)}| \end{eqnarray}

(3)

facilitating balancing the sensitivities of transverse |$m_j^{(y)}$| and longitudinal |$m_j^{(x)}$| movements by adjusting the weighting ratio β/γ.

Figure 6.

[From Wilkinson et al. (2010), figure 2, x, z-coordinates corrected:] 2-D resistivity image inverted from the baseline data set (2008 March). The inferred boundaries between the Whitby (WMF), Staithes (SSF) and Redcar (RMF) formations are shown by dotted black lines. Stratigraphic logs of boreholes are shown in grey scale. The main scarp and slipped WMF material are indicated by the black arrows. [Reproduced from Wilkinson et al. (2010), fig. 2, for comparison with Fig. 11(a) in this work. The RES2DINV reconstruction region, elsewhere in this paper referred to as the ‘RES2DINV outline,’ is selected by the software based on a heuristic pseudo-section method (Loke 2017).]

4 METHOD

While the work of Wilkinson et al. (2010, 2016), described in the previous section, focused on a sequential inversion procedure, here we simultaneously reconstruct resistivity and electrode position in a Gauss–Newton iterative framework. A 2-D model of the slope profile was constructed with independent parameters for resistivity and electrode position (Fig. 7). The ℓ²-norm of the data discrepancy and regularization terms was minimized by balancing the sensitivity of the resistivity parameters ρ and electrode movement x against the regularization terms

\begin{eqnarray} (\hat{\mathbf {\rho }},\hat{{\mathbf {x}}}) &=& \arg \min _{\mathbf {\rho },{\mathbf {x}}} \ ||\mathcal {F}(\mathbf {\rho },{\mathbf {x}}) - {\mathbf {m}}||_{\mathbf {W}} \nonumber \\ &&\quad + \, ||\lambda \beta {\mathbf {R}}_\mathbf {\rho }\log _{10}(\mathbf {\rho }_\star /\mathbf {\rho })||_{2} \nonumber \\ && \quad + \, || \lambda \eta {\mathbf {R}}_{{\mathbf {x}}} ({\mathbf {x}}- {\mathbf {x}}_\star )||_{2}. \end{eqnarray}

(4)

The optimal solution |$(\hat{\mathbf {\rho }},\hat{{\mathbf {x}}})$| minimizes the data discrepancy between a forward model |$\mathcal {F}$| and measurements m combined with some regularization R_ρ and R_x acting to smooth an otherwise ill-posed and ill-conditioned inverse problem. The forward model was the parametrized 2.5-D model (Fig. 7) of resistivity ρ and longitudinal² electrode position x and producing an estimate of expected measurements given the measurement protocol. Measurement error was weighted W based on estimates of measurement reliability. The regularization terms penalized changes in resistivity and electrode position from prior estimates (ρ_⋆, x_⋆). The relative sensitivity of the two types of parameters, resistivity and electrode movement, were balanced by adjusting the ratio of the scalar β/η. The overall strength of the regularization was adjusted by scaling both terms by λ. The resistivity was solved under a positive conductivity σ constraint by converting to inverse log units log ₁₀σ = log ₁₀1/ρ.

Figure 7.

The 2.5-D forward model FEM parametrization; right-to-left expanding from the region surrounding a single electrode, to the scale of the electrode array, to the scale of the region surrounding the electrode array, the forward model is parametrized for electrode position x and resistivity ρ, resistive regions selected to demonstrate mesh structure.

The objective function (4) was solved using the well-known iterative Gauss–Newton approach (Nocedal & Wright 1999,  section 10.2). The Gauss–Newton approach starts from an initial estimate (ρ₀, x₀), estimates the local slope of the objective function as the Jacobian (J_σ, J_x) to determine a new search direction (δρ, δx), and then performs an approximate line search in that direction to estimate an optimal step length α. At each iteration, the parameters were updated

\begin{eqnarray} (\delta \mathbf {\rho }, \delta {\mathbf {x}}) = -({{\bf J}}^{{\sf T}} {\mathbf {W}}{{\bf J}}+ \lambda ^2 {\bf Q})^{-1}({{\bf J}}^{{\sf T}} {\mathbf {W}}{\mathbf {b}}\, + \lambda ^2 {\bf Q} {\bf c}) \end{eqnarray}

(5)

\begin{eqnarray} &&{(\mathbf {\rho }_{n+1}, {\mathbf {x}}_{n+1}) = (\mathbf {\rho }_n, {\mathbf {x}}_n) + \alpha (\delta \mathbf {\rho }, \delta {\mathbf {x}})} \\ &&{{\rm for }\, \,{\mathbf {b}}= \mathcal {F}(\mathbf {\rho }_n,{\mathbf {x}}_n) - {\mathbf {m}}}\nonumber \\ &&{{\bf c} = \left[ \begin{array}{c}\log _{10} 1/\mathbf {\rho }_n - \log _{10} 1/\mathbf {\rho }_\star \\ {\mathbf {x}}_n-{\mathbf {x}}_\star \end{array} \right] = \left[ \begin{array}{c}\log _{10} \mathbf {\rho }_\star /\mathbf {\rho }_n\nonumber\\ {\mathbf {x}}_n-{\mathbf {x}}_\star \end{array}\right]}\nonumber \\ &&{{{\bf J}}= \left[ \begin{array}{c}\frac{-1}{\mathbf {\rho }_n \ln (10)} \, {{\bf J}}_{\sigma _n} \\ {{\bf J}}_{\mathbf {x}}\end{array} \right]} \nonumber \\ &&{{\mathbf {R}}= \left[ \begin{array}{c{@}{\quad}c} \beta {\mathbf {R}}_\mathbf {\rho } & 0 \\ 0 & \eta {\mathbf {R}}_{\mathbf {x}}\end{array}\right] \qquad {\rm and \, \, \, {\mathbf {R}}^{{\sf T}} {\mathbf {R}}= {\bf Q}}}\nonumber \end{eqnarray}

(6)

based on the data discrepancy b and the distance from the prior estimate |$\bf c$| in combination with the regularization R, the slope of the objective function J and measurement weighting W. This formulation agrees with that of Boyle et al. (2017), but is modified to address a log resistivity parametrization. In contrast to a typical resistivity-only inversion, the reconstruction parameters, Jacobian and regularization R have been extended to incorporate the new electrode position parameters. The resistivity Jacobian has been calculated for conductivity |${{\bf J}}_{\sigma _n}$| and then scaled, which is exactly equivalent to calculating the Jacobian on the log of resistivity.

The movement Jacobian was found to be sensitive to resistivity changes in the small elements adjacent to the electrodes. To address this sensitivity, the log conductivity regularization combined a smoothing prior near electrodes with Tikhonov regularization away from the electrodes

\begin{eqnarray} {\mathbf {R}}_\rho = {\mathbf {I}}+ \nu {\mathbf {R}}_{\mathfrak {L}} \end{eqnarray}

(7)

\begin{eqnarray} \nu = \exp -\, \frac{|x_e - x_\ell |}{|\bar{x}_{\ell \ell }|} \end{eqnarray}

(8)

where the Laplace smoothing |${\mathbf {R}}_{\mathfrak {L}}$| was scaled ν by the distance between each FEM element centre x_e and the nearest electrode x_ℓ, scaled by the average distance between electrodes |$\bar{x}_{\ell \ell }$|⁠. This prior encourages small changes from the expected resistivity in regions with little information. In regions near the electrodes, changes will be pushed towards the spaces between electrodes rather than directly under the electrodes, as well as encouraging smooth transitions in resistivity near the electrodes. The regularization for movement was the Tikhonov prior R_x = I. In principle, there are correlated changes between resistivity and electrode movement (Kim et al. 2014) which may be partially accounted for by setting the off-diagonal blocks of the regularization matrix to non-zero values, but in practice these were not characterized and in the absence of a better guess were set to zero.

The forward model was constructed as a 2-D cross-section based on the original electrode locations and the mesh was then perturbed by PCHIP³ interpolation (Carlson & Fritsch 1989) for electrode displacements. Forward modelled measurements and the Jacobians were calculated using the 2.5-D method (Dey & Morrison 1979).

We have made use of the log conductivity to restrict the resistivity reconstruction to physically meaningful positive values. Experiments with the log movement constraint, to restrict electrodes to downslope movement, resulted in similar reconstructions to the ones presented here which used unrestricted electrode movement. For the log movement parametrization, the behaviour at each iteration was different to an unscaled movement due to the structure of the movements in this data set. Because each electrode that moved had a different magnitude of movement, the log movement reconstruction tended to solve for each electrode’s reconstructed displacement separately: one electrode per iteration. The apparent single electrode updates were actually an artefact of the log scaling, where smaller movements were reduced by orders of magnitude, so as to be inconsequential. Once the largest electrode placement error had been corrected, the next largest error would be addressed. We feel that this highlights the importance of careful selection of the reconstruction parametrization. It is possible that some Fourier decomposition or other basis of electrode movement with appropriate regularization might achieve greater reconstruction accuracy without artificially fixing any single electrode’s location or degrees of freedom for movement.

5 2.5-D POSITION JACOBIAN

The 2-D electrode position Jacobian suffers from significant errors (Fig. 8), when compared to data from a 3-D model, which makes it inappropriate for 3-D problems. A 3-D electrode position Jacobian becomes prohibitively expensive to calculate as the mesh density grows. A 2.5-D approach offers a compromise by restricting sensitivity parametrization and electrode positions to the plane, while achieving high fidelity to equivalent 3-D models, at a fixed multiple of the 2-D computational effort.

Figure 8.

Jacobian sensitivity |${\rm diag}\,\,({{\bf J}}_x^{{\sf T}} {{\bf J}}_x)^\frac{1}{2}$| on a 16-electrode, homogeneous (σ = 1) half-space model; 2-D rank-one electrode position Jacobian (Gómez-Laberge & Adler 2008) shows orders-of-magnitude error in sensitivity estimate, while 3-D analytic (9) and rank-one estimates (Gómez-Laberge & Adler 2008) are in close agreement with the 2.5-D (13) estimate.

Two alternate methods were evaluated before developing the 2.5-D position Jacobian: a 2.5-D perturbation method which was relatively slow, and an analytic model of movement which was restricted to a homogeneous resistivity and the Point Electrode Model (PEM). Motivated by the short comings of these two methods, we develop the 2.5-D position Jacobian which is efficient and accounts for resistivity variation in the model using the CEM.

The perturbation method used the underlying 2.5-D forward simulations of the full FEM resistivity model using the CEM. These movement perturbation calculations proved to be prohibitively slow because a mesh perturbation resulted in recalculations of the system matrices and a new inversion of that system matrix. The cost grows with number of electrodes and movement dimensions, so that for n = 32 electrodes, estimated in d = 1 dimensions, nd + 1 forward simulations were required at each iteration of the Gauss–Newton inverse solution. A line search typically required three to six forward simulations, meaning that the perturbation Jacobian required far more time to calculate than the rest of each iteration. When FEM mesh density was increased sufficiently to achieve good estimates of electrode position changes, the calculations took an unreasonable amount of time. Low mesh densities sped up the calculations but exhibited significant errors when compared to a 3-D perturbation solution at sufficient mesh density.

An alternate solution was implemented by adapting the half-space analytic PEM forward model. For a half-space, the potential difference V measured over a homogeneous resistivity ρ with current I driven on the stimulus electrodes, is given by

\begin{eqnarray} V = \frac{\rho I}{2 \pi } \left( \frac{1}{AM} - \frac{1}{BM} - \frac{1}{AN} + \frac{1}{BN} \right) \end{eqnarray}

(9)

where each distance AM, BM, AN, BN is between a stimulus electrode and a measurement electrode. The model may be applied for any arbitrary pair-wise electrode placement. A position Jacobian may be constructed by applying the electrode movement perturbation. A next logical step would be to take the derivative of (9) with respect to electrode position and build up a solution specific approximate block-wise model of resistivity away from electrodes (Wilkinson et al. 2010), though this was not implemented in this work. Electrode positions were captured from the FEM model. A homogeneous resistivity was assigned based on the average resistivity of the current FEM model. The electrode position Jacobian produced by the half-space analytic perturbation method was compared to the 2.5-D perturbation Jacobian under homogeneous conditions. Using the modified half-space analytic perturbations was much faster to calculate than the 2.5-D perturbation method and reasonably accurate: the effect of topography was found to be somewhat accounted for. The loss of accuracy due to changing electrode models (CEM to PEM) and using a homogeneous resistivity were not so disruptive as to change signs in the position Jacobian although magnitudes were inaccurate.

Motivated by the efficiency of the 2-D position Jacobian of Gómez-Laberge & Adler (2008), and the errors introduced by using a 2-D electrode position Jacobian, the 2.5-D position Jacobian was developed as the corollary of the 2.5-D conductivity Jacobian. In general, the 2.5-D forward solver is a well-known technique and is commonly used in geophysics ERT applications. An approximately half-space geometry, and a resistivity that is nearly uniform along one axis, fit well with a 2.5-D model, and occur naturally in many geological settings. By adapting the adjoint, or ‘standard method,’ of calculating the conductivity Jacobian by a rank-one update, to the 2.5-D technique, resistivity may be efficiently reconstructed. To reconstruct electrode movement, we desire a similar 2.5-D implementation for the electrode position Jacobian. In the following, we outline the 2.5-D conductivity Jacobian and present a new derivation for the 2.5-D electrode position Jacobian. We use the formulation of the 2-D conductivity Jacobian from Boyle et al. (2017) as a basis for these developments.

The 2-D and 2.5-D conductivity Jacobians J_σ were calculated

\begin{eqnarray} {{\bf J}}_{\sigma ,{\rm 2D}} = - {\mathcal {T}}{\mathbf {A}}^{-1} {\mathbf {C}}^{{\sf T}} {\mathbf {S}}\frac{\partial {\mathbf {D}}}{\partial \sigma _n} {\mathbf {C}}{\mathbf {X}} \end{eqnarray}

(10)

\begin{eqnarray} {{\bf J}}_{\sigma ,{\rm 2.5D}} = - \frac{2}{\pi } \int _k {\mathcal {T}}{\mathbf {A}}_k^{-1} {\mathbf {C}}^{{\sf T}} ({\mathbf {S}}+ k^2 {\mathbf {T}}) \frac{\partial {\mathbf {D}}}{\partial \sigma _n} {\mathbf {C}}\frac{{\mathbf {X}}_k}{2} \end{eqnarray}

(11)

for measurement selection |${\mathcal {T}}$|⁠, system matrix A, mesh connectivity matrix C, mesh shape functions S, a conductivity change ∂D/∂σ_n, and the nodal voltages X = A⁻¹B for stimulus B over an electrode modelled as a shunt in the y-direction when the 2-D FEM is meshed over the x–z plane. The system matrix |${\mathbf {A}}= {\mathbf {C}}^{{\sf T}} {\mathbf {S}}{\mathbf {D}}{\mathbf {C}}$| is assembled from a connectivity matrix C mapping global node numbers to element-local node numbers, the element shape functions S and the conductivity D per element.

For the 2.5-D position Jacobian, the system matrix A_k is specific to the spatial frequency k, as are the nodal voltages |${\mathbf {X}}_k = {\mathbf {A}}_k^{-1} {\mathbf {B}}$|⁠. A perturbation node is selected at row u, column v, affecting a linear interpolatory shape function E. The 2.5-D position Jacobian may be calculated as an extension of the 2-D Jacobian, in an analogous way to (10), (11), as

\begin{eqnarray} {{\bf J}}_{x,{\rm 2D}} = - {\mathcal {T}}{\mathbf {A}}^{-1} {\mathbf {C}}^{{\sf T}} \frac{\partial {\mathbf {S}}}{\partial x_n} {\mathbf {D}}{\mathbf {C}}{\mathbf {X}} \end{eqnarray}

(12)

\begin{eqnarray} {{\bf J}}_{x,{\rm 2.5D}} = - \frac{2}{\pi } \int _k {\mathcal {T}}{\mathbf {A}}_k^{-1} {\mathbf {C}}^{{\sf T}} \frac{\partial ({\mathbf {S}}+ k^2 {\mathbf {T}})}{\partial x_n} {\mathbf {D}}{\mathbf {C}}\frac{{\mathbf {X}}_k}{2} \end{eqnarray}

(13)

where the 2-D position Jacobian may be efficiently calculated using the rank-one update for the conductivity Jacobian (Gómez-Laberge & Adler 2008) with some new terms.

Again, based on the 2-D formulation from Boyle et al. (2017), the element shape functions for element (e) may be summarized as the element shape matrix

\begin{eqnarray} {\mathbf {S}}^{(e)} = \frac{1}{2|\det {\mathbf {E}}|}{\mathbf {E}}_{\setminus 1}^{{\sf T}} {\mathbf {E}}_{\setminus 1} \end{eqnarray}

(14)

for a shape matrix E and a row-reduced version E_∖1 where the top row of the matrix is removed. The shape matrix is distorted by having its nodes perturbed leading to the first-order estimate of element deformation

\begin{eqnarray} \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} &=& \frac{1}{2}\, \Bigg ( \frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} {\mathbf {E}}_{\setminus 1}^{{{\sf T}} } {\mathbf {E}}_{\setminus 1} \, \nonumber\\ &&\qquad + \frac{1}{|\det {\mathbf {E}}|} \bigg (\frac{\partial {\mathbf {E}}_{\setminus 1}^{{{\sf T}} }}{\partial x_n} {\mathbf {E}}_{\setminus 1} + {\mathbf {E}}_{\setminus 1}^{{{\sf T}} }\frac{\partial {\mathbf {E}}_{\setminus 1}}{\partial x_n} \bigg ) \Bigg ) \end{eqnarray}

(15)

where x_n refers to a global node numbered n and affects all elements e connected to that node. The local shape functions of each element, for first-order interpolatory shape functions on a 2-D mesh, are

(16)

for a triangular element (blue) with three nodes p₁, p₂, p₃.

To calculate the partial derivatives of the first-order interpolatory shape functions, we make use of the matrix determinant lemma for an invertible square matrix H where

\begin{eqnarray} \det ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} ) = \det ({\bf H})(1+{\mathbf {v}}^{{\sf T}} {\bf H}^{-1}{\mathbf {u}}). \end{eqnarray}

(17)

The update uses the rank-one perturbation vectors u and v, selecting by row and column, to manipulate a single element of the matrix, a node of our mesh, by a small perturbation. A first-order approximation of the derivative of a determinant via a rank-one perturbation is then

\begin{eqnarray} \frac{\partial \det ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} ) }{\partial x_n} = \det ({\bf H})\, {\mathbf {v}}^{{\sf T}} {\bf H}^{-1}{\mathbf {u}}. \end{eqnarray}

(18)

To evaluate the change in our shape function’s determinant, we use the partial derivative of an absolute function ∂|H| = H∂H/|H| and the inverse determinant equivalence |$\det ({\bf H}^{-1})=\det ({\bf H})^{-1}$| so that

\begin{eqnarray} \frac{\partial \, |\det {\mathbf {E}}|^{-1} }{\partial x_n}= \frac{\partial \, |\det ({\mathbf {E}})^{-1}|}{\partial x_n} = \frac{\partial \, |\det ({\mathbf {E}}^{-1})|}{\partial x_n} \end{eqnarray}

(19)

\begin{eqnarray} \phantom{\frac{\partial \, |\det {\mathbf {E}}|^{-1} }{\partial x_n}}= \frac{\det ({\mathbf {E}}^{-1}) \, \frac{\partial (\det {\mathbf {E}}^{-1})}{\partial x_n} }{|\det ({\mathbf {E}}^{-1})|} \end{eqnarray}

(20)

and the partial derivative of the determinant

\begin{eqnarray} \frac{\partial \det ({\mathbf {E}}^{-1}) }{\partial x_n} = \det ({\mathbf {E}}^{-1})\, {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}} \end{eqnarray}

(21)

can be applied to (20) after reducing the determinants

\begin{eqnarray} \frac{\det ({\mathbf {E}}^{-1}) \det ({\mathbf {E}}^{-1})}{|\det ({\mathbf {E}}^{-1})|} = \frac{|\det {\mathbf {E}}|}{\det ({\mathbf {E}})^2} = \frac{1}{|\det {\mathbf {E}}|} \end{eqnarray}

(22)

so that

\begin{eqnarray} \frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} &= \frac{ {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}}{|\det {\mathbf {E}}|} \end{eqnarray}

(23)

The partial derivative of the reduced shape matrix can be approximated using the Sherman–Morrison formula

\begin{eqnarray} ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} )^{-1} = {\bf H}^{-1} - \frac{{\bf H}^{-1} {\mathbf {u}}{\mathbf {v}}^{{\sf T}} {\bf H}^{-1}}{1+{\mathbf {v}}^{{\sf T}} {\bf H}^{-1} {\mathbf {u}}} \end{eqnarray}

(24)

to get a rank-one update

\begin{eqnarray} \frac{\partial {\mathbf {E}}_{\setminus 1} }{\partial x_n} = -({\mathbf {E}}{\mathbf {u}}{\mathbf {v}}^{{\sf T}} {\mathbf {E}})_{\setminus 1} \end{eqnarray}

(25)

for a small perturbation such that |${\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}\ll 1$|⁠.

To go from a 2-D solution to a 2.5-D solution, a correction term k²T appears in the system matrices so that the shape matrices S^(e) are extended to become S^(e) + k²T^(e) for spatial wave-number k.

\begin{eqnarray} {\mathbf {T}}^{(e)} = \frac{1}{2|\det {\mathbf {E}}|} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray}

(26)

For the 2.5-D position Jacobian, this change adds a new partial derivative term ∂T^(e)/∂x_n

\begin{eqnarray} \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} \, \rightarrow \, \frac{\partial ({\mathbf {S}}^{(e)}+k^2{\mathbf {T}}^{(e)})}{\partial x_n} = \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} + k^2 \frac{\partial {\mathbf {T}}^{(e)}}{\partial x_n} \end{eqnarray}

(27)

where ∂S^(e)/∂x_n is already available from the 2-D calculations, and the additional term ∂T^(e)/∂x_n may be derived

\begin{eqnarray} \frac{\partial {\mathbf {T}}^{(e)} }{\partial x_n} = \frac{1}{2}\frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray}

(28)

but we already have |$\partial |\det {\mathbf {E}}|^{-1} / \partial x_n$| from deriving the partial derivatives of the S^(e) term (23) giving

\begin{eqnarray} \frac{\partial {\mathbf {T}}^{(e)} }{\partial x_n} = \frac{1}{2}\frac{ {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}}{|\det {\mathbf {E}}|} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray}

(29)

for a linear interpolatory shape function E perturbing a node at row u and column v.

This adjoint or rank-one perturbation method for the 2.5-D electrode position Jacobian may be calculated much more quickly than a direct perturbation method because the system matrices do not need to be recalculated and inverted to determine the change in measurements due to electrode movement.

The 2.5-D electrode position Jacobian (13) was found to be 25.9 times faster than the equivalent 3-D rank-one update method (Gómez-Laberge & Adler 2008), implementing (12) in 3-D, for mesh geometries used in this work (Intel Core i5-2500K 4-core processor at 3.30 GHz with 32 GB memory). 2-D electrode position Jacobian estimates differ significantly from 3-D solutions (Fig. 8), so have not been presented in Fig. 9. The computational cost of the 2-D rank-one update method for calculating the electrode position Jacobian (Gómez-Laberge & Adler 2008) was orders of magnitude faster than a naïve 2-D perturbation method. The 2-D rank-one update was 7.1 times faster than the 2.5-D method across most mesh sizes, which is accounted for by the numerical integration implied by (13). There are likely to be further gains from optimizing this implementation for multiple processing cores because key portions of the Jacobian calculation (14), (15), (26), (29) can be performed in parallel and the Jacobian typically consumes a significant portion of the total calculation time in each Gauss–Newton iteration (Boyle et al. 2012b). Mesh density was measured as the inverse of the maximum element height h (elements per metre) for both 2-D and 3-D meshes. The relative speed-up for a particular mesh density 1/h grows as a function of the number of FEM mesh nodes n and elements which must be calculated in the Jacobian where n = O(h⁻²) in 2-D and n = O(h⁻³) in 3-D; the larger the mesh the greater the benefit conferred by the 2.5-D Jacobian approach.

Figure 9.

The 2.5-D Jacobian calculations scale with FEM node count: the 2.5-D movement Jacobian speed advantage over a rank-one 3-D calculation (Gómez-Laberge & Adler 2008) grows as mesh density (FEM elements per metre) 1/h, where h is the maximum element height for the entire mesh. Error bars show max/min run times over 20 runs, the run time is closely related to the number of nodes in the FEM mesh where 2-D meshes have far fewer nodes for the same mesh density.

6 RESULTS

The column ℓ²-norm sensitivity (the diagonal of |${{\bf J}}^{{\sf T}} {{\bf J}}$|⁠) was plotted by replacing reconstructed resistivity with the log of estimated sensitivity. Sensitivity in these plots was expected to be greatest near the electrodes and diminish elsewhere. Simple Dirichlet boundary conditions away from the electrodes, used in these simulations, introduced errors, which were observable as variations in sensitivity at unexpected locations. We wished to model an approximately half-space forward model but unexpected increases in sensitivity near the sides and bottom were found to be caused by the boundary conditions which were deflecting current flow. Boundary condition errors can be corrected in a number of ways, the simplest of which is to increase the modelled domain until the error is small enough. One could, alternatively, implement appropriate ‘infinite elements’ at the boundary (Babus̆ka 1972). Another method is to estimate a ‘primary’ field for each stimulus using an analytic half-space model or very detailed one-time-use FEM mesh, and then calculate a ‘secondary’ field update on a smaller FEM with different resistivity as a correction (Günther et al. 2006). The primary-secondary type methods rely on small changes in resistivity far enough from the boundary to leave the ‘primary’ field largely unperturbed and it is not immediately obvious how this method may be affected by electrode movements or surface deformation without recalculating the primary field at each update. Neumann and mixed boundary conditions away from the electrodes were not considered. We have used an expanded model, in the interests of reliable results under deformed boundaries, at the expense of some lost computational efficiencies. Regions of high sensitivity were initially noted at depth where little sensitivity was expected. To determine how far the FEM model boundaries needed to be extended, an analytic PEM half-space model was compared to CEM homogeneous resistivity FEM simulated measurements. The model boundaries were extended approximately one electrode array length in each of the +x, −x and −z directions. The boundary extension reduced boundary condition related errors in simulated measurements to within measured noise levels and removed the artefacts from the sensitivity plots.

The resistivity sensitivity S was plotted relative to the maximum sensitivity,

\begin{eqnarray} {\bf S} = {\rm diag}({{\bf J}}_\sigma ^{{\sf T}} {\mathbf {W}}{{\bf J}}_\sigma )/{\bf V} \end{eqnarray}

(30)

using the conductivity Jacobian J_σ, measurement inverse covariance/weighting W, and element volumes |$\bf V$|⁠, for line#1 (March 2008) and line#5 (February 2013) using the surveyed locations (Fig. 10). The region near the electrodes has been presented with annotations matching Fig. 6, as well as an image of the complete model. The initial and final resistivities were independently reconstructed using the surveyed locations (Fig. 11) and the difference between initial and final was used to create the expected resistivity change (Figs 12e and f). The initial resistivity for line#1 closely matched those published in Wilkinson et al. (2010) (Fig. 6) and achieved a similar <1 per cent RMS measurement misfit relative to a homogeneous resistivity model. Qualitatively and quantitatively, the two reconstructions (Figs 6 and 11a) are very similar.

Figure 10.

Inverse model parametrization with colours showing relative sensitivity |${\bf S} = {\rm diag}({{\bf J}}_\sigma ^{{\sf T}} {\mathbf {W}}{{\bf J}}_\sigma )/{\bf V}$| as S/max (S) for homogeneous resistivity; (a) line#1 March 2008, and (b) line#5 February 2013; note the distinct difference in slope profile between (a) and (b).

Figure 11.

Resistivity-only reconstructions using true electrode positions as measured by RTK GPS; (a) line#1 March 2008 to (b) April 2009 (λ = 52.1, σ₀ = 31.3 Ωm), and (c) line#5 February 2013 to (d) February 2014 (λ = 42.7, σ₀ = 25.6 Ωm).

Figure 12.

Change in resistivity and electrode movement for joint movement reconstructions (λ_σ = 0.1, λ_x = 0.07), (a–d) line#1, March 2008 to April 2009, and (e–h) line#5, February 2013 to February 2014.

Electrode movements were initially reconstructed without allowing resistivity change using an independent implementation of the electrode movement iterative solver. The observed behaviour of the algorithm (Fig. 13) was to first minimize the error in the electrode spacing (corresponding to the largest measurement misfit), and then to approach a more ‘correct’ solution. Artefacts of this approach to the minimum over the optimization surface remain where large steps in the true electrode movement are reconstructed as a balanced step by adjacent electrodes that come close to the true displacement between electrodes near that movement. The joint inversion code was checked against this result by setting the movement-resistivity balance parameter β to strongly favour electrode movement. Reconstructions showed no resistivity change and movements that were in close agreement with the movement-only reconstruction code. Small variations still existed between the two results due to differences in the Gauss–Newton implementation and inexact line search. These variations were small with respect to the overall electrode movement solution. When electrode movements were reconstructed with resistivity changes (Fig. 12), some portion of the reconstructed electrode movement was lost in favour of reconstructed resistivity change.

Figure 13.

Electrode movement without allowing for resistivity changes, iterations for (a) line#1, March 2008 to April 2009, and (b) line#5 movement, February 2013 to February 2014.

Due to the large electrode movements, it was found helpful to perform a crude version of successive relaxation. The first three iterations of the Gauss–Newton reconstruction were performed with an electrode movement hyperparameter that was an order of magnitude larger than following iterations. Without this adaptation, the reconstructed electrode movements showed poor agreement with measured locations, presumably because the Gauss–Newton iterations were trapped in a local minimum which favoured constructing resistivity artefacts near the electrodes. Exploring the space of hyperparameters near the selected hyperparameter did not reveal one which achieved better electrode movement reconstruction.

7 DISCUSSION

Resistivity was reconstructed for measured initial and final electrode locations (Fig. 11) which serve as an ‘ideal’ reconstruction. Resistivity and electrode displacement were simultaneously reconstructed for a survey located on a slowly moving landslide. Results exhibit some measure of oscillatory artefacts in the reconstructed movement (Figs 12c and g).

The resistivity distribution for line#1 (March 2008–April 2009; Figs 11a and b) changed by a relatively small amount when using true electrode locations before and after movement. This would suggest that, beyond the ground motion at the surface, no structural changes in the near surface seem to have occurred. It seems plausible that the increased area of low resistivity WMF might be indicative of increased saturation of the soil which led to the translational slide of WMF material moving over SSF substrate at the surface.

The resistivity changes for line#5 (February 2013–February 2014; Figs 11c and d), show a significantly different distribution after ground movement, which is interesting given that the two electrode lines are within 40 m of each other. The line#5 measurements (Fig. 11d) occurred after a very wet summer and winter period where there was a lot of seepage at the base of the lobe causing the deeper reductions in resistivity between z = 40 and 60 m. At the surface of the lobes, resistivity increased due to cracking of the top layer. A cracked surface experienced accelerated evaporation over increased surface area, resulting in localized resistivity increase. To a lesser extent, areas affected by surface cracking also showed increased resistivity due to the change in topography: relatively little current would be conducted across the air gap in cracks, contributing to an average increase in bulk conductivity, but the cracks do increase the surface along which current flows resulting in an effective increase in resistivity. We take particular note of the change in SSF resistivity around x = 80 m which may have developed vertical connectivity between the overlying WMF and RMF below, allowing vertical drainage. The flow would be from the saturated WMF, along the WMF-SSF boundary to the surface, then to a region of vertical connectivity downward through the SSF (x = 80 m), and then into the RMF where it has pooled underground. This proposed flow path might also explain the increase in resistivity in the SSF at x = 60 m: if the vertical connectivity were in a roughly vertical plane, it would cut off the outer section of the SSF and that outer section would drain downwards into the RMF leading to an increase in resistivity. The deeper segment of the SSF (x > 80 m) would maintain its resistivity because the general connectivity and saturation have not changed by much. The resistivity change may also be induced by model error in electrode placement or topography: the 2.5-D model limits fidelity in some respects. The poor quality of the line#5 post-movement data, as measured by reciprocal error, may be the cause of these changes, though the locations of the reciprocal errors along the electrode array were distributed along the length of the array so that we expect no concentrated region of low sensitivity that could cause resistivity changes in the reconstruction such as those observed Fig. 11(d). Given this speculation, it would be interesting to investigate this potential vertical fault in the SSF. Perhaps it is an indication of a major ground movement, still to come, as the lower slope drainage has changed significantly. The change in drainage may also help to stabilize the lower slope by providing a drainage path at-depth which will allow surface material in the SSF to consolidate. This stabilizing effect has been observed on other nearby slopes in previous years. It has been postulated in Uhlemann et al. (2017) that reactivation of the slope at line#1 was stabilized due to slope movement which caused preferential flow paths to open, lowering pore pressures on the slip surface of the lobe, thereby stabilizing the lobe.

For a half-space model with a homogeneous resistivity, electrode positions are not unique. A translation of the entire set of electrodes will give identical measurements. Likewise, a scaling of all electrode positions is equivalent to a change in the homogeneous resistivity. When conductivities are inhomogeneous the electrode locations are somewhat fixed by the locations of the inhomogeneities. Examples of electrode position non-uniqueness manifested itself in this data as large oscillations in the reconstructed electrode movement when no measures were taken to address the issue.

Fixing the location of three electrodes at the upslope and downslope ends of the electrode array (Figs 12c and g) nearly eliminated these oscillations. As seen in the line#5 data, this is not necessarily a correct assumption, as both the top and bottom of a landslide may move, leading to resistivity artefacts. We infer that fixing these electrode locations was sufficient to damp the reconstructed movement’s oscillations because it fixes the relationship between a stimulus current, a measured voltage, and electrode separation (distance). Smoothing-type regularization of resistivity, used in this work, then controls how strictly the selected scaling of electrode separation is enforced. For example, electrodes that have contracted together in a region might cause a conductive artefact to be reconstructed under those electrodes. Increasing the resistivity regularization may suppress these artefacts and cause movement to be reconstructed by contracting the local electrode spacing to account for smaller than expected voltage measurements in the region.

Both unscaled and log scaled electrode movement were tested and found to give results with similar absolute positional error. We have elected to present the unscaled electrode movement in our reconstructions as it is a less restrictive choice. One could imagine low angle slopes, wetlands or floodplains where the expected direction of movement would not be known a priori. Peat wetlands, muskeg, and most ground exposed to deep frost experience seasonal expansion and contraction, due to freeze-thaw cycles, water table changes, or water and gas accumulation and evaporation, which can result in uplift and ground shifting in directions other than downslope (Taber 1930; Hansell et al. 1983; Price 2003; Strack et al. 2006; Uhlemann et al. 2016c). A further reason to avoid dependence on the log movement scaling is the extension of this work to transverse movements where the restriction to movements only to one side of the array seems inappropriate. In the data sets examined here, there are transverse movements which were caused by material accumulating at the toe of the landslide and towards the edges of the earth flow. These transverse movements can be predicted for this particular data set based on the pre-existing topology: line#1 electrodes moved east, downhill into a gully, while the line#5 electrodes moved west, downhill into the same gully.

Reconstructions for line#1 generally matched the true electrode locations within 0.20 m for movements of up to 1.46 m with the exception of electrode #9 and the three electrodes #6–#8 at the step in electrode position (Fig. 12a). Compared to Wilkinson et al. (2010) (≤0.2 m position error), these results are marginally less accurate. It seems probable that our results for the line#1 data differ from those of Wilkinson et al. (2010) due to the simultaneous resistivity and electrode position reconstruction, presented here, which removes artificial ordering constraints that occur with the sequential method of Wilkinson et al. (2010). Another source of differences in our results with respect to Wilkinson et al. (2010) is the restriction to downslope movements by a log parametrization in Wilkinson et al. (2010). As mentioned previously, trials of this log parametrization method in our simultaneous resistivity and electrode position inversion did not lead to improved results. Our results for the line#5 data appear to closely correlate with Wilkinson et al. (2016), where reconstructed electrode locations were generally within 0.2 m excepting some electrodes with errors up to 1.0 m, though results are not presented in as much detail in that case. In contrast to Wilkinson et al. (2016), our reconstructions do not address movements transverse to the electrode line.

Movement reconstructions for line #5 do not appear to be particularly accurate, perhaps due to the more significant resistivity changes inferred in the reconstruction and more widespread translational failure of the slope which shifted electrodes over most of the resistivity structure (Fig. 12g). These might be addressed by identifying the covariance between movement and resistivity change within the joint reconstruction algorithm regularization. It is also possible that some of the error in reconstructed electrode position may be due to the FEM discretization. An approach such as the Fréchet method for electrode movement may help to alleviate such issues, though in general it produces the same solutions as a 3-D rank-one update method (Dardé et al. 2012; Boyle et al. 2017). Adjusting the relationship between resistivity and movement regularization β caused greater electrode displacement error as resistivity regularization was reduced. These movement magnitudes represented movement of up to 32 per cent of the average 4.73 m electrode spacing, exceeding the joint resistivity-movement methods of Soleimani et al. (2006) which was limited to movements of approximately 1 per cent of electrode spacing.

8 CONCLUSIONS

This work demonstrated the practical application of a joint electrode movement-resistivity reconstruction using an iterative Gauss–Newton regularized solver. The electrode position Jacobian was calculated on the current resistivity at each iteration. Reconstructions show reasonable agreement with RTK GPS measured electrode locations, available resistivity estimates and geological structure.

The initial reconstructed resistivity model, used as a starting point for the electrode movement and resistivity reconstruction, was in close agreement with prior work (Figs 6 and 11a) (Wilkinson et al. 2010). Reconstructed changes in resistivity (Fig. 12) showed considerable variation, particularly around the region at the toe of the landslide. These changes in resistivity could be indicative of water saturation changes due to water seepage at the toe of the landslide or other geological causes. Another possibility is that the resistivity changes represent artefacts due to transverse and normal components of electrode movement which were not accounted for in these reconstructions.

We note that, in general, even when electrode displacements were not entirely accurate compared to true electrode positions, the error in the estimated and true distances between electrode positions was quite accurate after one or two Gauss–Newton iterations. Errors in electrode spacing were distributed fairly evenly across the electrode array, after which the displacements shifted towards their true positions in most cases. This suggests that a parametrization for electrode movement that encompasses electrode spacing may lead to improved outcomes. The results suggest that electrode grids are effective not only for resistivity monitoring, but also as a means of ground motion detection which may provide a cost-effective approach for landslide monitoring.

We are encouraged by these results and expect that with new protocols which measure between electrode lines, higher quality electrode position reconstructions will be possible which incorporate normal, lateral, and transverse electrode movements, as observed in the data sets presented here.

Acknowledgements

Identification of the Hollin Hill site, geology, field work and electrical measurements were carried out by the British Geological Survey. Thanks go out to the researchers at the British Geological Survey, Geophysical Tomography Team for generously sharing their time and data. They were supported by the Natural Environment Research Council (NERC) and their contributions to this work are published with the permission of the Executive Director of the British Geological Survey. In fond memory of Steve Gibson, and with thanks to Josie Gibson.

Code for the algorithms presented in this work extend EIDORS, an open source Matlab toolkit for impedance imaging inverse problems (Adler & Lionheart 2006; Adler et al. 2015, 2017), and use NetGen for meshing (Schöberl 1997).

This work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC).

Footnotes

1

Leica System 1200 RTK GPS in ‘kinematic mode’ real-time correction achieves as much as 10 mm (Root Mean Squared (RMS)) horizontal and 20 mm vertical accuracies (Merritt 2014).

2

Longitudinal movement being movement inline with the electrodes and along the surface.

3

Matlab interp1(X,Y,Xq,‘pchip’) 1-D interpolant based on Hermite derivatives.