This manual describes the workflow, file formats, main parameters, and module roles for the RFDFWIM (Radar Frequency-Domain Full-Waveform Inversion using MATLAB) package.
-
Create models and acquisition
Runcreate_models_mkl.mto generate:- Permittivity model:
input/MYMOD_eps.dat - Conductivity model:
input/MYMOD_sig.dat - Source positions:
input/MYMOD_source.dat - Receiver positions:
input/MYMOD_receiver.dat
- Permittivity model:
-
Forward modelling (generate “observed” data)
RunGPRFM.m. It:- Reads the above input files
- Runs 2D FDFD forward modelling with CFS-PML
- Saves synthetic data as
obs/obs_data.mat(variableprecrnew) - Optionally produces wavefield and model figures in
output/
-
Full-waveform inversion
RunRFDFWI.m. It:- Loads
obs/obs_data.matas observed data - Builds a smooth initial model from the true model (Gaussian filter)
- Runs FWI (gradient + optional Hessian, L-BFGS, Tikhonov, Wolfe line search)
- Saves models, gradients, wavefields, and figures in
outputmyINV/
- Loads
This is an inverse crime setup: the same forward code generates the data and is used inside the inversion.
| Script | Role |
|---|---|
| create_models_mkl.m | Builds two-cross (2C) synthetic models and 4-sided acquisition geometry; writes MYMOD_*.dat into input/. |
| GPRFM.m | Forward modelling driver: reads input, runs FDFD for all shots/frequencies, saves obs/obs_data.mat and optional outputs in output/. |
| RFDFWI.m | Main FWI driver: reads observed data and input models, runs iterative FWI, writes all inversion outputs to outputmyINV/. |
| inp_GPRmodel.m | Defines model parameters used mainly for conductivity (grid, PML, file names, plotting). |
| inp_GPRmodel1.m | Defines model and inversion parameters (grid, frequencies, PML, CFS-PML, L-BFGS, Wolfe, Tikhonov, paths). |
| File | Role |
|---|---|
| forward_shotNEW.m | Core forward routine: for each frequency, builds system matrix with CFS-PML, solves for wavefield, samples at receivers. Returns wavefield te and data precrnew. |
| RHS_TE1.m | Builds right-hand side for TE mode for a given shot and frequency. |
| PML_9p_CFS_stag_para.m | CFS-PML parameters for 9-point stencil. |
| imp_A_TE_9p_cfs_PML_stag1_para.m | Assembles impedance matrix for TE, 9-point, CFS-PML, staggered grid (variant 1). |
| imp_A_TE_9p_cfs_PML_stag2_para.m | Same as above, staggered variant 2 (MKL). |
| read_acq_jan.m | Reads input/MYMOD_receiver.dat and input/MYMOD_source.dat; returns acquisition struct (e.g. acq.xrec, acq.yrec, acq.xshot, acq.yshot, acq.nxrec, acq.nyrec, acq.nxshot, acq.nyshot, acq.ntr, acq.nsrc). |
| shift_acq.m | Shifts acquisition coordinates when PML is used (extended grid). |
| readeps1.m | Loads permittivity model from input/MYMOD_eps.dat. |
| readsig1.m | Loads conductivity model from input/MYMOD_sig.dat. |
| extend_eps1JAN.m | Extends ε (and σ) model with PML padding. |
| extract_model_ext.m | Extracts physical domain from extended (PML) model. |
| dotp.m | Dot product utility (e.g. for inner products). |
| File | Role |
|---|---|
| grad_obj_hessMKLnew.m | Computes forward wavefield, adjoint wavefield, data residual, L2 objective, gradients (σ, ε), and pseudo-Hessian (σ, ε) for TE mode. |
| cal_radargMKLnew.m | Simulated data at receivers from forward wavefield. |
| cal_resjMKLnew.m | Residual and L2 contribution per frequency/shot; optional Wiener weighting. |
| RHS_adj_TEm.m | Right-hand side for adjoint TE solve. |
| RHS_TE_hess.m | Right-hand side for Hessian-related solve. |
| ass_grad_TEMKLnew.m | Assembles gradient contributions from wavefields. |
| copy_grad_frame_TE.m | Copies gradient into frame (e.g. for PML handling). |
| norm_matrix_TE.m | Norm of gradient matrix (for scaling). |
| scale_grad_TE.m | Scales material parameters (σ, ε) for inversion. |
| Tikhonov_cost_TE.m | Tikhonov regularisation term in cost function. |
| Tikhonov_grad_TE.m | Gradient of Tikhonov term. |
| steep_descent.m | Steepest-descent direction (negative gradient). |
| LBFGS_TEmNEW.m | L-BFGS quasi-Newton update; returns preconditioned descent directions for σ and ε. |
| pseudo_hess_shin_ass_TEmNEW.m | Pseudo-Hessian (Shin et al.) assembly for TE. |
| check_descent_TE.m | Checks if search direction is a descent direction; can trigger L-BFGS reset. |
| wolfe_TENEW.m | Wolfe line search; updates σ and ε and returns step length. |
| calc_mat_change_wolfe_multi_para_TE.m | Material parameter update for multi-parameter Wolfe. |
| plot_wavefld.m | Plots wavefield for figures. |
| dotp_matrix.m | Matrix inner product. |
-
MYMOD_eps.dat, MYMOD_sig.dat
ASCII matrices: one value per grid point (row/column order as increate_models_mkl.mandreadeps1/readsig1). Size must matchmodel1.nx,model1.ny(e.g. 180×180). -
MYMOD_receiver.dat, MYMOD_source.dat
Two columns: x (distance), y (depth) in metres. One row per receiver or source. Read byread_acq_jan.m.
- output/ — Forward modelling figures and exports (e.g. true model, wavefield).
- outputmyINV/ — Inversion: figures (e.g.
figs/), saved models (models/), gradients (grad/,grad_tikh/,gradm/), wavefields (wavefield/), Hessian-related (Hgrad/), and.matfiles (e.g.FinalL2.mat). - obs/ — Contains
obs_data.mat(variableprecrnew) produced by GPRFM.
- model1.nx, model1.ny — Number of grid points in x and y (e.g. 180).
- model1.dh — Grid spacing [m] (e.g. 0.05).
- model1.npml — PML thickness in grid points (e.g. 10).
- model1.free — Free surface flag (0 = no free surface in tests).
- model1.eps0 — Vacuum permittivity (8.85e-12).
- model1.sig0 — Reference conductivity for scaling (e.g. 5.6e-3).
- model1.mu0 — Magnetic constant.
In RFDFWI.m the frequency vector is set as:
f = [50e6 60e6 70e6 80e6 90e6 100e6 125e6 150e6 175e6 200e6]; % HzSimultaneous multi-frequency inversion (e.g. Lavoué et al., Layek & Sengupta).
- itermax — Max FWI iterations (e.g. 1500).
- model1.wolfe — 1 = use Wolfe line search; otherwise fixed step.
- model1.LAMBDA_1, model1.LAMBDA_2 — Tikhonov weights.
- model1.nlbfgs — L-BFGS history length.
- model1.C1, model1.C2 — Wolfe constants.
- model1.SRTRADIUS — Source taper radius [m].
- model1.sigr_low, model1.sigr_high, model1.epsr_low, model1.epsr_high — Bounds for σ and ε.
- model1.STAG — 1 = staggered variant 1; 2 = staggered variant 2.
- model1.ACQMY — 1 = 4-sided acquisition (as in create_models_mkl).
- Source/receiver files:
input/MYMOD_source.dat,input/MYMOD_receiver.dat(see read_acq_jan.m).
The main loop stops when any of the following holds:
- ratio = L2(iter) / L2(1) ≤ 5e-5
- iter = itermax
- diff ≤ PRO (relative change of L2, with PRO=0 in the default script)
Final models are written and figures exported when the loop exits.
- seismic.map — MATLAB colormap file in the root directory; used for model and wavefield plots (e.g.
load('seismic.map'); colormap(flipud(seismic))).
- Lavoué et al. (2014) — FWI of GPR data.
- Layek & Sengupta (2019, 2021, 2023) — Multi-parameter GPR FDFD FWI.
- Meles et al. (2011) — GPR FDTD adjoint, two-cross style setup.
- Shin et al. (2001) — Pseudo-Hessian.
- Nocedal & Wright (2006) — L-BFGS and Wolfe conditions.
For more detail on algorithms and parameters, see the header comments in the corresponding .m files and the README/INSTALLATION/COMMANDS documents.