A Joint Numerical Simulation Method for Mine Seismic-Electric Coupling

1. Introduction

Coal mining in China is progressively extending to greater depths, where increasing geological complexity elevates the risk of water inrush accidents caused by concealed structures—such as collapse columns and faults—posing serious threats to safety and personnel. Precise characterization of the spatial distribution and water-richness of these hidden hazards is a prerequisite for scientific water hazard control and safe pressurized mining. Currently, borehole DC resistivity and seismic exploration have emerged as core techniques for detecting concealed hazards underground, owing to their complementary physical principles.

The mine DC resistivity method, based on contrasts in electrical conductivity, is highly sensitive to low-resistivity anomalies such as water-bearing bodies and fractured zones, enabling effective identification of water-rich structures in roofs and floors. However, constrained by underground observation systems and algorithmic limitations, it suffers from insufficient resolution, making it difficult to precisely delineate anomaly boundaries; this often leads to anomaly diffusion and false anomalies, particularly in scenarios with multiple coexisting targets. Mine seismic exploration, relying on seismic impedance contrasts, offers high resolution and signal-to-noise ratio, excelling at characterizing structural interfaces and rock mass integrity. Nevertheless, its response to variations in water content is weak, rendering it ineffective for distinguishing between dry and water-saturated structures. The inherent limitations of single geophysical methods result in significant non-uniqueness in data interpretation, failing to meet the demand for high-precision exploration in modern smart mines.

Joint inversion, which integrates the advantages of multiple methods to achieve data complementarity and mutual constraints, has become a key technology to overcome the bottlenecks of single-method detection. Since Vozoff and Jupp (1975) pioneered true joint inversion, related research has evolved from one-dimensional layered models to three-dimensional complex structures and from linear to nonlinear approaches [1]. Lines et al.(1988) employed the least-squares method for joint inversion and compared the merits of sequential and synchronous strategies [2]. In 1995, Hering et al(1995). proposed the first 1D joint inversion formulation for DC resistivity and seismic surface-wave data, laying the theoretical groundwork for subsequent seismic-electric inversion [3]. Kis et al.(2001) further realized the joint inversion of multiple wave types and DC apparent resistivity for 2D geological bodies [4]. Gallardo and Meju(2003) have conducted extensive systematic research in this field; their proposed cross-gradient constraint theory avoids presupposing empirical relationships between petrophysical parameters. Instead, it enhances the structural consistency of different property models to achieve joint inversion, significantly improving accuracy and reliability [5]. This method has since been widely applied in studies such as 3D joint inversion of potential field data.

In China, Wang Jialin and colleagues summarized the application conditions for joint inversion through extensive research [6]. Chen Xiao and Yu Peng introduced regularization and adaptive techniques into synchronous joint inversion, enhancing solution stability and computational efficiency [7]. Zhou Lifen employed cross-gradient theory to couple resistivity and seismic velocity data, achieving 2D synchronous joint inversion of magnetotelluric and seismic data [8]. Gao Ji and Zhang Haijiang developed a joint inversion scheme for seismic traveltime and DC resistivity data with cross-gradient constraints, effectively mitigating shadow zone issues in seismic body-wave tomography and DC resistivity imaging [9]. Peng Miao, Tan Handong et al. conducted three-dimensional joint inversion of magnetotelluric and seismic travel time data with cross-gradient constraints. The study demonstrated that joint inversion enhanced the recovery degree of deep velocity structure characteristics of the abnormal body [10].These domestic studies have confirmed the feasibility of applying cross-gradient theory to seismic-electric joint inversion. However, a critical assessment shows that existing studies primarily focus on general geophysical scenarios. Significant gaps remain regarding adaptation to the confined spaces and complex boundaries (e.g., roadway interference) of underground mines, and there is a lack of targeted research on typical hazardous structures like water-conducting collapse columns and faults. Moreover, the critical balance between inversion accuracy and computational efficiency under deep mining conditions remains unresolved. Concurrently, the advancement of intelligent mining has intensified electromagnetic interference, challenging the robustness of current methods against practical operational conditions.

The inherent complementarity between mine DC resistivity and seismic exploration provides a natural foundation for addressing the challenge of accurately detecting concealed hazardous structures. Consequently, this paper focuses on developing a mine-adapted seismic-electric joint inversion method. We aim to optimize the structural coupling constraints and inversion algorithms. Through numerical simulations based on theoretical models of various anomalies and typical disaster-inducing structures, we seek to establish a joint inversion framework suited to complex mine conditions. This work is intended to provide a reliable technical pathway for the precise detection of concealed geological hazards, thereby offering robust support for smart mine construction and safety in production.

2. Individual Inversion

2.1. Mine DC Resistivity Method

2.1.1. Forward Modeling

The forward model is constructed using the finite element method (FEM) with an unstructured triangular mesh. The computational domain is discretized into triangular elements to accurately conform to complex geological geometries and roadway boundaries in the mine environment (Figure 1) [11].

Based on the minimum energy principle of electrostatics, the boundary value problem of the potential distribution governed by partial differential equations is transformed into a variational problem [12]. Discretization yields a large sparse linear algebraic system:

$$\nabla \cdot \left(\sigma \nabla \mathit{U}\right)=-2\mathit{I}\delta \left(A\right) \in \mathsf{\Omega },$$

(1)

$${\displaystyle \frac{\partial \mathit{U}}{\partial n}}=0 \in {\mathsf{\Gamma }}_s,$$

(2)

$${\displaystyle \frac{\partial \mathit{U}}{\partial n}}+{\displaystyle \frac{\cos \left(r,n\right)}{r}}\mathit{U}=0 \in {\mathsf{\Gamma }}_{\mathrm{\infty }}$$

(3)

where σ denotes the electrical conductivity of the subsurface medium, δ represents the Dirac delta function, n is the coordinate variable in the outward normal direction of the boundary, r is the position vector from the source point to a boundary point, Ω signifies the computational domain, and Γs, Γ∞ refer to the ground surface and subsurface boundaries, respectively.

Linear basis functions are adopted for interpolation within each triangular element. Integral results of all elements are assembled to derive the global matrix system:

$$\begin{array}{ll}F\left(U\right)& =\sum {\mathrm{F}}_e\left(U\right)=\sum {\displaystyle \frac{\sigma }{2}}{U}_e^{\mathrm{T}}\left(K_{1e}+K_{2e}\right){\mathit{U}}_e-{\mathit{U}}_{A}I\\ & ={\displaystyle \frac{1}{2}}{U}_e^{\mathrm{T}}KU-{\mathit{U}}^{\mathrm{T}}P\end{array}$$

(4)

where U denotes the global nodal potential vector. According to the variational extremum principle, setting the first variation of Equation (4) to zero generates the linear algebraic system:

$$\mathit{K}\mathit{U}=\mathit{P}$$

(5)

The symmetric successive over-relaxation (SSOR) preconditioned conjugate gradient (CG) method is adopted for efficient solution. An SSOR preconditioner is constructed to improve the condition number of the stiffness matrix K. Combined with sparse matrix storage techniques, an optimized SSOR-preconditioned iterative algorithm is developed referring to the incomplete Cholesky decomposition preconditioned biconjugate gradient strategy, and the nodal potential distribution U of the study area is ultimately obtained. This solution scheme guarantees computational accuracy while significantly enhancing the solving efficiency for large-scale 3D problems.

2.1.2. Inversion Algorithm

The inversion problem is formulated as a regularized least-squares optimization problem. The objective function is established based on the Occam inversion algorithm to balance data fitting accuracy and model smoothness [13]:

$$\mathsf{\Phi }=(\Delta \overrightarrow{d}-J\overrightarrow{m}{)}^T(\Delta \overrightarrow{d}-J\overrightarrow{m})+\lambda (R\overrightarrow{m}{)}^T\left(R\overrightarrow{m}\right)$$

(6)

where Δd is the residual vector between field observations and theoretical forward data, m denotes the model parameter vector, J represents the sensitivity (Jacobian) matrix, λ is the dynamically adjusted regularization parameter, and R signifies the model roughness matrix.

The model parameter update equation is expressed as:

$$\Delta \overrightarrow{m}=(J^{T}J+\lambda R^{T}R{)}^{-1}J\Delta \overrightarrow{d}$$

(7)

In DC resistivity inversion, the root mean square (RMS) error is a critical indicator for evaluating the consistency between model responses and field measurements. It quantifies the overall deviation between theoretical calculations and observed data, providing an objective criterion for inversion reliability. The inversion convergence is defined as the relative change rate of RMS error being less than 2%, and the RMS error is formulated as:

$$RMS=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{i=1}^N({\displaystyle \frac{d_i-f(\overrightarrow{m}{)}_i}{d_i}}{)}^2}\times 100$$

(8)

where N is the number of observation data points, ${\mathit{d}}_{\mathit{i}}$ is the $\mathit{i}$-th measured datum, and $\mathit{f}(\overrightarrow{\mathit{m}}{)}_{\mathit{i}}$ represents the theoretical response corresponding to the model parameters.

2.2. Mine Seismic Exploration

2.2.1. Forward Modeling

Unstructured mesh discretization is performed via Delaunay triangulation, employing squared slowness (the reciprocal of velocity squared) as the model parameter to adapt to the complex geological interfaces in the mine [14,15]. Based on Fermat’s principle of minimum traveltime, a curved-ray tracing method is utilized to simulate seismic wave propagation paths, iteratively correcting ray trajectories to align with true paths via the variational principle [16].

The governing forward equation system is expressed as:

$${\sum }_{j=1}^M{L}_{ij}{s}_j=t_i(i=\mathrm{1,2},\dots ,N)$$

(9)

In Eq. (9),${\mathit{L}}_{\mathit{i}\mathit{j}}$denotes the length of the$\mathit{i}$-th ray traversing the$\mathit{j}$-th grid cell,${\mathit{s}}_{\mathit{j}}$represents the squared slowness of the$\mathit{j}$-th cell,${\mathit{t}}_{\mathit{i}}$is the traveltime of the $\mathit{i}$-th ray,$\mathit{M}$is the total number of grid cells, and $\mathit{N}$is the total number of rays. Solving this system yields theoretical traveltime datasets.

2.2.2. Inversion Algorithm

The LSQR algorithm (least-squares QR decomposition) is adopted to solve large sparse linear systems. By exploiting matrix sparsity, the computational load is drastically reduced, memory consumption is minimized, and numerical instability inherent in conventional least-squares methods is avoided; the algorithm exhibits linear convergence characteristics. This method is based on the Lanczos bidiagonalization process, iteratively constructing orthogonal matrices to transform the original system into bidiagonal form, subsequently solved via QR decomposition techniques [17].

The inversion objective function is formulated as:

$$\min \Vert Ax=b{\Vert }_2$$

(10)

where$\mathbf{A}$is the coefficient matrix containing ray segment lengths, x denotes the squared slowness parameter vector, and b is the traveltime residual vector. The algorithm iteratively constructs orthogonal matrices and updates model parameters until the residuals decline to the observational error level and satisfy predefined accuracy criteria.

3. Joint Inversion

3.1. Joint Inversion Model

A synchronous joint multi-parameter inversion (SJIMM) framework is employed, constructing a target function that integrates a data-fitting term, a model regularization term, and a structural constraint term to achieve collaborative inversion of DC resistivity and seismic data [18]:

$${\phi }_{SJIMM}\left(\mathbf{m}\right)={\displaystyle {\sum }_{i=1}^A} {\phi }_{\left(d,i\right)}\left({\mathbf{m}}_{\left(j\right)}\right)+{\displaystyle {\sum }_{j=1}^B} {\mu }_{\left(j\right)}{\phi }_{\left(m,j\right)}\left({\mathbf{m}}_{\left(j\right)}\right)+{\displaystyle {\sum }_{l=1}^C} {\beta }_{\left(l\right)}{\psi }_{\left(l\right)}\left(\mathbf{m}\right)$$

(11)

where the parameters are defined as follows:

${\phi }_{\left(d,i\right)}\left(\mathit{m}\right)$: residual norm;

A: total number of independent observation datasets;

B: number of inverted petrophysical parameters;

C: number of coupling constraints between parameters;

$\mathbf{m}={\left[{\mathbf{m}}_{\left(1\right)},{\mathbf{m}}_{\left(j\right)},\dots ,{\mathbf{m}}_{\left(B\right)}\right]}^T$: model parameter vector;

${\mu }_{\left(j\right)}$: regularization parameter;

${\beta }_{\left(l\right)}$: weight coefficient of the$l$-th model constraint function;

${\psi }_{\left(l\right)}\left(\mathbf{m}\right)$:$l$-th model constraint function;

${\phi }_{\left(m,j\right)}\left({\mathbf{m}}_{\left(j\right)}\right)$: model regularization term.

When the inversion system involves only two datasets and two petrophysical parameters, the generalized joint inversion framework (Eq. (11)) simplifies to the specific form:

$$\begin{gathered} {\phi }_{SJIMM}\left(\mathbf{m}\right)={\phi }_{\left(d,1\right)}\left(\mathbf{m}\right)+{\phi }_{\left(d,2\right)}\left(\mathbf{m}\right) \\ +{\mu }_{\left(1\right)}\left({\phi }_{\left(m,1\right)}\left({\mathbf{m}}_{\left(1\right)}\right)+{\displaystyle \frac{{\mu }_{\left(2\right)}}{{\mu }_{\left(1\right)}}}{\phi }_{\left(m,2\right)}\left({\mathbf{m}}_{\left(2\right)}\right)+{\displaystyle {\sum }_{l=1}^c} {\displaystyle \frac{{\beta }_{\left(l\right)}}{{\mu }_{\left(1\right)}}}{\psi }_{\left(l\right)}\left(\mathbf{m}\right)\right) \end{gathered}$$

(12)

3.2. Structural Constraint Conditions

Cross-gradient constraints are employed to establish structural consistency between resistivity and velocity parameters, eliminating the need for presupposed empirical relationships between physical properties. This enhances inversion accuracy by strengthening the synchrony of dual property models [19]. The constraint function is expressed as:

$${\phi }_{CG}\left({\mathbf{m}}_{\left(1\right)},{\mathbf{m}}_{\left(2\right)}\right)={‖\nabla {\mathbf{m}}_{\left(1\right)}\times \nabla {\mathbf{m}}_{\left(2\right)}‖}_2^2$$

(13)

In numerical computation, the cross-gradient function for the $\mathit{i}$-th triangular mesh element is:

$${\phi }_{C{G}^i}={‖\begin{array}{ccc}\overrightarrow{x}& \overrightarrow{y}& \overrightarrow{z}\\ {\displaystyle \frac{\partial {m_1}^i}{\partial x}}& {\displaystyle \frac{\partial {m_1}^i}{\partial y}}& {\displaystyle \frac{\partial {m_1}^i}{\partial z}}\\ {\displaystyle \frac{\partial {m_2}^i}{\partial x}}& {\displaystyle \frac{\partial {m_2}^i}{\partial y}}& {\displaystyle \frac{\partial {m_2}^i}{\partial z}}\end{array}‖}^2={\left(t_{CG}^{x,i}\right)}^2+{\left(t_{CG}^{y,i}\right)}^2+{\left(t_{CG}^{z,i}\right)}^2$$

(14)

where$t_{CG}^{x,i}$, $t_{CG}^{y,i}$, and$t_{CG}^{z,i}$represent the cross-gradient components along the x,$\mathrm{y}$, and z directions, respectively, computed via the vector cross product of model parameter gradients. This design ensures synchronized variations of different petrophysical parameters at structural boundaries.

3.3. Inversion Solution Procedure

The Gauss–Newton optimization algorithm is utilized to minimize the joint inversion target function. The iterative update formula is:

$$\begin{gathered} \left({\mathbf{J}}^{\left(n\right)T}{{\mathbf{W}}_d}^T{\mathbf{W}}_d{\mathbf{J}}^{\left(n\right)}+\mathit{\mu }{{\mathbf{W}}_m}^T{\mathbf{W}}_m+\beta {\mathbf{H}}_{CG}^{\left(n\right)}\right)\mathsf{\Delta }{\mathbf{m}}^{\left(n\right)} \\ ={\mathbf{J}}^{\left(n\right)T}{{\mathbf{W}}_d}^T{\mathbf{W}}_d\left({\mathbf{d}}^{obs}-{\mathbf{f}}^{\left(n\right)}\left(\mathbf{m}\right)\right)-\mathit{\mu }{{\mathbf{W}}_m}^T{\mathbf{W}}_m\left({\mathbf{m}}^{\left(n\right)}-{\mathbf{m}}^{ref}\right)-{\mathbf{P}}_{CG}^{\left(n\right)} \end{gathered}$$

(15)

where ${\mathbf{W}}_d$ denotes the data weight matrix,${\mathbf{W}}_m$is the model weight matrix,${\mathbf{m}}^{ref}$represents the reference model,${\mathbf{H}}_{CG}$is the Hessian matrix of the cross-gradient term, and${\mathbf{P}}_{CG}$is the gradient vector of the cross-gradient term.

4. Numerical Simulations

4.1. Seismic-Electric Joint Inversion Numerical Simulations

The DC resistivity method is highly sensitive to subsurface water content and fluid distribution, enabling inference of rock connectivity and fluid properties from resistivity variations. Seismic exploration responds to rock elastic moduli and density, allowing characterization of geological structures and fragmentation degrees via wave velocity and attenuation analyses. Physically realistic models and survey geometries representative of underground mine conditions are constructed. Model domains, anomaly sizes, and spatial positions are defined to match practical mining scenarios. Models are discretized via unstructured meshes to accommodate complex mine geometries and roadway boundaries, with seismic sources and geophones deployed in the working face roadways for realistic data acquisition [20]. Specifically, mine seismic wave detection employs a multi-channel seismic transmission imaging technique, with geophone arrays deployed at 10 m intervals along the roadways on both sides of the working face. The observation system adopts an equal-interval source excitation and full-array reception mode, with seismic source points arranged at 50 m intervals. A single-point excitation and full-array reception strategy is implemented to acquire high-quality seismic wavefield data, ensuring sufficient spatial coverage and angular illumination of the target zone. Mine DC resistivity detection employs a dipole-dipole configuration for data acquisition, with current (source) and potential (measuring) electrodes deployed at 10 m intervals along the sidewalls of the roadways to obtain potential difference data. This configuration provides high sensitivity to lateral resistivity variations and is well suited for detecting water-bearing structures ahead of the working face. Based on the above observation systems, numerical simulation studies are carried out separately to analyze the individual inversion effects of the two geophysical methods, providing a basis for subsequent joint inversion and comprehensive interpretation.

Model 1: Multiple Isolated Anomalies

Model comprises three low-resistivity anomalies embedded in a high-resistivity background and three corresponding low-velocity anomalies in a high-velocity background.

Figure 2. DC resistivity and seismic exploration model of Model 1(Upper resistivity model, lower velocity model).

Figure 2. DC resistivity and seismic exploration model of Model 1(Upper resistivity model, lower velocity model).

The model is discretized via unstructured triangular meshes, and synthetic DC and seismic data are generated for inversion tests. Single-method inversions reveal that resistivity-only inversion produces distorted anomaly contours, pronounced smearing, and spurious anomalies, whereas seismic-only inversion approximates the anomaly shapes but exhibits significant vertical smearing.

Figure 3. Single inversion of model 1(Upper resistivity profile, lower velocity profile).

Figure 3. Single inversion of model 1(Upper resistivity profile, lower velocity profile).

Joint inversion rectifies anomaly contours, suppresses smearing, eliminates spurious anomalies, and reconstructs anomaly spatial features with high fidelity, markedly outperforming single-method inversions.

Figure 4. Joint inversion of Model 1.

Figure 4. Joint inversion of Model 1.

Model 2: Fractured Water-Conducting Channel

Model simulates a linear water-conducting fracture, featuring a single elongated low-resistivity anomaly and a matching low-velocity anomaly.

Figure 5. DC resistivity and seismic exploration model of Model 2.

Figure 5. DC resistivity and seismic exploration model of Model 2.

In the single resistivity inversion result, the low-resistivity anomaly appears elongated, its sharp terminations at both ends are not clearly identified, and spurious anomalies are present. In the single seismic inversion result, the sharp ends of the low-velocity anomaly are similarly not delineated accurately. (Figure 6)

Joint inversion refines anomaly extents, accurately captures sharp channel terminations, and restrains smearing, demonstrating high precision for linear hydrogeological conduit imaging.

Figure 7. Joint inversion of Model 2.

Figure 7. Joint inversion of Model 2.

4.2. Summary

The seismic-electric joint inversion is implemented following a standardized workflow: data input, model initialization, unstructured mesh discretization, optimization of the joint inversion objective function via the Gauss–Newton method, and model parameter update through linear search, until the convergence criteria are satisfied.

Typical concealed hazard-inducing geological structures are constructed as model anomalies, and a comparative analysis is conducted focusing on the performance of single-method inversions (mine DC resistivity method and mine seismic exploration) versus seismic-electric joint inversion. The comparison demonstrates that by integrating multiple geophysical datasets, joint inversion effectively reduces the solution non-uniqueness inherent in single-method inversions and improves the accuracy and reliability of inversion results.

5. Conclusions

This study addresses the limitations of low resolution and high non-uniqueness inherent in single geophysical methods for detecting concealed hazardous structures in underground mines. By combining theoretical derivation, algorithm development, and numerical validation, we establish a comprehensive seismic–electric joint inversion framework. The main conclusions are summarized as follows.

(1) Single-method inversion frameworks tailored to mining conditions have been developed. For the DC resistivity method, we couple unstructured-mesh finite-element forward modeling with Occam-style regularized inversion, which reliably reconstructs the positions and resistivity parameters of anomalous bodies. For the seismic method, we adopt curved-ray forward modeling combined with the LSQR algorithm, which accurately images velocity anomalies. These frameworks lay a solid foundation for the subsequent joint inversion.

(2) A cross-gradient-constrained seismic–electric joint inversion scheme has been formulated. The scheme unifies the data-fitting term, regularization term, and cross-gradient term within a single objective function. By enforcing structural consistency between the resistivity and velocity models without assuming any empirical petrophysical relationship, the proposed scheme effectively suppresses spurious anomalies and image smearing, thereby improving the reliability of the inverted solutions.

(3) The effectiveness of the joint inversion has been verified through numerical experiments on representative mining scenarios, including multiple isolated anomalies, fractured water-conducting channels, and water-bearing faults. Compared with the single-method inversions, the proposed approach substantially improves the morphological reconstruction of anomalies, eliminates spurious features, and reduces solution non-uniqueness, offering a practical tool for the precise detection of concealed hazard-inducing structures in coal mines.