## 3 Numerical Experiments  
To test the performance of different regularization strategies in dual-parameter FWI for GPR, we consider synthetic 3D models, as shown in Figure 2a,c, emulating a high-contrast subsurface scenario comprising four blocky anomalies. The dimensions of this model in the \(x\), \(y\), and \(z\) directions are 6.0 m, 6.0 m, and 6.0 m, respectively, using discrete cube elements with a size of 0.12 m \(\times\) 0.12 m \(\times\) 0.12 m. The outer domain around the calculation area is a 4-cell PML. The background model consists of three layers, which, from top to bottom, are the air layer and the first and second background layers; the interfaces between these layers are roughly located at 0 m and 2.4 m. The sizes of the block-shaped anomalous bodies of A, B, C, and D are 0.60 m \(\times\) 0.60 m \(\times\) 0.72 m, 0.60 m \(\times\) 0.60 m \(\times\) 0.72 m, 1.08 m \(\times\) 0.60 m \(\times\) 0.84 m, and 0.60 m \(\times\) 1.08 m \(\times\) 0.84 m, respectively. The electrical material properties of this model are summarized in Table 1.  
Figure 2: True models and initial models: (**a**) true relative permittivity model; (**b**) initial relative permittivity model; (**c**) true conductivity model; (**d**) initial conductivity model. A, B, C, and D in (**a**,**c**) denote the positions of four block-shaped anomalous bodies. Red circles in (**b**) and white \(\times\)-shaped markers in (**d**), respectively, represent the transmitters and receivers in the data acquisition setup for 3D GPR. As shown in Figure 2b,d, the initial models were the backgrounds of the true models. The observation system adopted a multi-offset mode with co-polarizations, and there were in total 10 \(\times\) 10 sources and 21 \(\times\) 21 receivers evenly distributed in the air (sources and receivers are dipoles in the \(x\) direction, and each receiver observes data from each source). In this case, the scaling factor \(\gamma\) in (8) was set to 1.8. As seen in the flowchart shown in Figure 1, we have to select several frequencies before performing 3D FWI. In our experiments, we assumed the dominant frequency to be 200 MHz and used 10 frequencies, namely 50, 60, 70, 80, 100, 120, 150, 180, 200, and 220 MHz, to perform the inversion. Figure 3 shows the behavior of the frequency-weighted index factor \(n\) in Equation (6) on different frequency components of GPR data. It indicates that the frequency weighting strategy can highlight the low-frequency components and suppress the high-frequency components in the FWI framework. Here, the frequency-weighted index factor \(n\) in (6) was set to 2. The frequency-dependent wavefields or functions were computed in parallel using MATLAB R2022a on a Dell Precision T7920 workstation with an Intel (R) Xeon (R) CPU E3-1280 v3 @ 2.70 GHz and with 512 GB memory. In addition, three iteration stopping criteria were set: the difference in the misfit between the first and the current iteration was less than 1 \(\times\) 10\({}^{-5}\) (preventing overfitting), the misfit did not decrease after the step line search was performed 20 times (avoiding local minima), or the iteration number was larger than 50 (controlling the computational cost). The inversion process was terminated when any of these three conditions was satisfied. Furthermore, to measure the reconstruction errors of the models, we adopt the relative root mean square (RRMS) function and the structural similarity index measure (SSIM) as  
\[RRMS_{m}=\frac{\sqrt{\sum\left(m^{inv}-m^{true}\right)^{2}}}{\sqrt{\sum\left(m^ {0}-m^{true}\right)^{2}}}, \tag{22}\]  
where \(m^{inv}\), \(m^{true}\), and \(m^{0}\) represent the values in the inverted, true, and initial models.