Forward and Inversion Combination Procedure

• The procedure is a computational method combining the scaled boundary finite element method (SBFEM) with full-waveform inversion to reconstruct spatial density from time-domain data in heterogeneous media.
• It employs semi-weak spatial discretization, adjoint-based gradient computation, and a two-level L-BFGS optimization to enhance accuracy and convergence.
• Innovations such as on-the-fly elemental condensation and parallel processing enable scalable performance and reduced computational cost.

The forward and inversion combination procedure refers to the integrated computational and algorithmic strategy for reconstructing spatial fields (such as density) from indirect time-domain measurements, specifically in the context of full-waveform inversion (FWI) employing the scaled boundary finite element method (SBFEM). This approach synthesizes a detailed semi-weak spatial-temporal discretization, adjoint-based gradient calculation, and advanced quasi-Newton optimization. It features key innovations in the condensed assembly of element-level matrices and hierarchical parameter updates, yielding high computational efficiency and parallelizability in large-scale heterogenous media (Daneshyar et al., 12 Jan 2025).

1. Governing Equations and Physical Setting

The objective is the estimation of a spatially varying density $\rho(x)$ in a domain $\Omega$ based on observations of a scalar wavefield $u(x,t)$ driven by a prescribed source $f(x,t)$. The heterogeneous acoustic wave equation is posed as: $$\mathcal{L}u = \nabla \cdot (\rho c^2 \nabla u) - \rho\,\partial_t^2 u = f(x,t)$$ for $x \in \Omega$, $t \in [0,T]$, with initial conditions $u(x,0) = 0$, $\partial_t u(x,0) = 0$, and free-surface boundary condition $\nabla u \cdot n = 0$ on $\partial\Omega$. Here, $c(x)$ denotes the spatially varying wave speed. This formulation ensures the full density-derivative terms are retained in the forward model, a necessity for high-fidelity inversion in media with non-constant density.

2. Scaled Boundary Semi-Weak Discretization

The SBFEM partitions the computational domain into polygonal elements, each parameterized by a radial coordinate $\xi \in [0,1]$ from a scaling center and a circumferential coordinate $\eta$. The spatial operators are mapped to this local coordinate system through chain-rule transformations of derivatives: $$\xi^2 E_0(\xi) \partial^2_\xi \mathbf{u} +\xi[E_0 + E_1^T - E_1 + \xi E_{0,\xi}] \partial_\xi\mathbf{u} - [E_2 - \xi E_{1,\xi}^T]\mathbf{u} - \xi^2 M_0(\xi) \partial_t^2 \mathbf{u} = 0.$$ The matrices $E_0(\xi), E_1(\xi), E_2(\xi), M_0(\xi)$ are boundary integrals over shape functions and Jacobian determinants, explicitly dependent on the local density and wave speed. Crucially, terms like $E_{0,\xi}$ and $E_{1,\xi}$ (radial derivatives of the coefficients) account for spatial variation of density, a capability absent in legacy SBFEM.

3. Adjoint-State Method and Fréchet Gradient Computation

Parameter estimation is cast as minimization of the data misfit functional: $$J[\rho] = \frac12 \sum_r \int_0^T (u(x_r, t) - d(x_r, t))^2 \, dt,$$ where $d(x_r, t)$ denotes receiver data. Introducing a Lagrange multiplier (adjoint wavefield) $\lambda(x,t)$ linked via the wave equation, one derives the adjoint problem: $$\nabla \cdot (\rho c^2 \nabla \lambda) - \rho \partial_t^2 \lambda = (u - d) \sum_r \delta(x - x_r)$$ with terminal conditions $\lambda(x,T) = 0$, $\partial_t \lambda(x,T) = 0$. The Fréchet derivative (sensitivity kernel) for $J$ with respect to $\rho$ is given by: $$\frac{\delta J}{\delta \rho}(x) = \int_0^T [ \partial_t \lambda \, \partial_t u - c^2 \nabla \lambda \cdot \nabla u ](x,t) \, dt.$$ Forward propagation is performed first, followed by time-reversed adjoint propagation. Zero-lag cross-correlations between $\lambda$ and $u$ at each spatial point accumulate the gradient with respect to $\rho$.

4. Radial Discretization and On-the-Fly Elemental Condensation

Temporal discretization of the system employs Newmark’s method; spatially, the radial coordinate is discretized using central finite differences on a grid $\xi_0=0 < \xi_1 < \cdots < \xi_n=1$: $$\partial_\xi \mathbf{u}_i \approx \frac{\mathbf{u}_{i+1} - \mathbf{u}_{i-1}}{2h}, \quad \partial^2_\xi \mathbf{u}_i \approx \frac{\mathbf{u}_{i+1} - 2\mathbf{u}_i + \mathbf{u}_{i-1}}{h^2}.$$ At each time step, the semi-discretized system yields a local tridiagonal linear system per element. An element-level Thomas algorithm (Gaussian elimination) serially eliminates interior unknowns, producing a $2 \times 2$ nodal system for the element’s boundary. The local dynamic stiffness matrices and equivalent nodal sources are assembled into the global FEM system. This “on-the-fly” assembly avoids introducing auxiliary variables in the global solve, reducing both memory and computational expense.

5. Two-Level L-BFGS Optimization and Parameter Hierarchy

The inversion cycle consists of: (i) assembling the full cost $J[\rho]$ via SBFEM forward simulation, (ii) computing the gradient with a single adjoint solve, and (iii) parameter update using a limited-memory BFGS scheme: $$\rho^{(k+1)} = \rho^{(k)} - \alpha_k H_k^{-1} \nabla J[\rho^{(k)}],$$ where $H_k^{-1}$ is the L-BFGS quasi-inverse Hessian constructed from the previous $m$ gradients and updates, and $\alpha_k$ is chosen by line search to ensure sufficient descent. The “two-level” optimization hierarchy refers to solving the PDEs and computing the gradient on a coarser mesh for speed, with the density parameters (and gradients) interpolated or mapped onto a finer grid for resolution enhancement.

6. Computational Decomposition and Parallel Efficiency

Each element’s dynamic stiffness is computed independently, reducing the local tridiagonal solves to $O(n)$ per element, and the global stiffness matrix assembly and solution to $O(N)$ per time step ($N=$ total nodes). Both steps are highly amenable to shared-memory parallelization: the element-local condensation, global matrix assembly, and time-stepping loops for both forward and adjoint propagation are all dispatchable over threads. On a 12-thread laptop (AMD Ryzen 7 PRO 6850U), the wall-clock for a full epoch (forward+adjoint+gradient+L-BFGS update) is $\approx2.4\,$s for a mesh with 2000 elements and 5 radial intervals per element, demonstrating robust scalability and strong performance relative to alternative formulations (Daneshyar et al., 12 Jan 2025).

7. Impact and Limitations

The fully coupled forward-inversion procedure in the SBFEM framework achieves accurate imaging in media with strong heterogeneity and spatially varying density, outperforming earlier methods (such as those neglecting density derivatives). The local condensation approach effectively decomposes the global computational load and enhances parallelizability. The method is suited for a range of scales and architectures, delivering high-fidelity solutions at substantially reduced computational cost. The main limitations pertain to the need for sufficiently fine radial discretization to resolve boundary effects and the adequacy of the SBFEM mesh for the spatial variability of $\rho(x)$. Use of the full, non-simplified differential operator is essential for inversion quality; neglecting density derivatives significantly degrades the solution (Daneshyar et al., 12 Jan 2025).