## 2 Theories and Methods
### 3D FWI Algorithm for GPR  
#### 2.2.1 Misfit and Gradient  
The objective function of FWI can be expressed as the L2 norm of the difference between the simulated data and the observed data. Assuming that the numbers of excitation sources and receiving points are FWI is \(N_{s}\) and \(N_{r}\), respectively, and the total number of  
Figure 1: The flowchart of the 3D FWI algorithm for GPR with the composite model constraint strategy. frequency components involved in the inversion is \(N_{\omega}\), substituting (3) and (4) into (1) gives [11; 28]  
\[\Phi_{d}(\mathbf{m})=\frac{1}{2}\sum_{i=1}^{N_{\omega}}\sum_{j=1}^{N_{\omega}} \sum_{k=1}^{N_{\omega}}\eta_{k}\Big{[}\mathbf{P}_{i,j}\mathbf{A}_{k}(\mathbf{ m})^{-1}\mathbf{q}_{i,k}-\mathbf{d}_{i,j,k}^{obs}\Big{]}^{\dagger}\Big{[} \mathbf{P}_{i,j}\mathbf{A}_{k}(\mathbf{m})^{-1}\mathbf{q}_{i,k}-\mathbf{d}_{i,j,k}^{obs}\Big{]}, \tag{5}\]  
where \({}^{\dagger}\) represents the transpose (\({}^{\dagger}\)) conjugate (*) operator, which is used to ensure that the value of the objective function is the true (real number) norm. The subscripts \(i,j,\) and \(k\), respectively, indicate that the value is related to the excitation source, the receiving point, and the frequency, and \(\eta_{k}\) is the weighting factor for data of different frequencies [47]:  
\[\eta_{k}=\left[\left(\frac{1}{\omega_{k}^{2}}\right)/\left(\sum_{k=1}^{N_{ \omega}}\frac{1}{\omega_{k}^{2}}\right)\right]^{n}, \tag{6}\]  
where \(n\) is the frequency-weighted index factor. We conduct this study on both the permittivity and the conductivity parameters simultaneously. Generally, the range of the relative permittivity \(\varepsilon_{r}\) is between 1 and 81, while the conductivity \(\sigma\) is mostly distributed in the order of \(10^{-5}\) to \(10^{-3}\) for the common detection objects of GPR. Considering the difference in the order of magnitude between the two, the relative conductivity \(\sigma_{r}=\sigma/\sigma_{0}\) and the parameter scale adjustment factor \(\gamma\) are introduced to ensure that \(\varepsilon_{r}\) and \(\sigma_{r}\) are of the same order of magnitude [11; 28]. At this point, the model parameters can be expressed as  
\[ \begin{bmatrix}\varepsilon\\ \sigma\end{bmatrix} = \begin{bmatrix}\varepsilon_{0}\varepsilon_{r}\\ \sigma_{0}\sigma_{r}\end{bmatrix} = \begin{bmatrix}\varepsilon_{0}\mathbf{m}_{ \varepsilon}\\ \gamma\sigma_{0}\mathbf{m}_{\varphi}\end{bmatrix} , \tag{7}\]  
where \(\varepsilon_{0}=8.854187817\times 10^{-12}\) (F/m) is the vacuo permittivity; \(\sigma_{0}\) is the reference value of the conductivity, which is approximately equal to \(\varepsilon_{0}\omega_{0}\); \(\omega_{0}=2\pi f_{0}\) is the angular frequency; and \(f_{0}\) is the center frequency of the antenna. On this occasion, we adopt the normalization strategy to update the model's parameters, which can be written as  
\[\mathbf{m}= \begin{bmatrix}\mathbf{m}_{\varepsilon}\\ \mathbf{m}_{\sigma}\end{bmatrix} = \begin{bmatrix}\varepsilon_{r}\\ \sigma_{r}/\gamma\end{bmatrix} . \tag{8}\]  
The parameter scale adjustment factor \(\gamma\) determines the solving process of the inversion, considering whether the inversion is dominated by \(\varepsilon_{r}\) or \(\sigma_{r}\) to a certain extent. Due to the insensitivity of \(\sigma_{r}\) to inversion and its tendency to produce inaccurate inversion results, the \(\varepsilon_{r}\) inversion result can be used to guide the \(\sigma_{r}\) inversion result by selecting an appropriate parameter scale adjustment factor \(\gamma\).