## 2 Theories and Methods
### 3D FWI Algorithm for GPR
#### 2.2.2 Composite Model Constraint Strategy  
In 3D FWI for GPR, the inversion results for the relative permittivity are usually better than those for the conductivity. The contour of the physical property interface is clearer and there are fewer artifacts. At this point, the conductivity inversion results show extremely high sensitivity during the inversion process and often present volume effects. The contours of the physical property interface are blurred, and there are many artifacts. For this reason, we introduce the composite model constraint FWI, loading the modified TV regularization for the permittivity and MSG regularization for the conductivity. The objective function of the model constraint term can be expressed as  
\[\Phi_{m}(\mathbf{m})=\lambda_{\varepsilon}TV(\mathbf{m}_{\varepsilon})+ \lambda_{\sigma}MSG\left(\mathbf{m}_{\sigma}-\mathbf{m}_{\sigma}^{0}\right), \tag{12}\]  
where \(\mathbf{m}_{\varepsilon}\) is the permittivity term, \(\mathbf{m}_{\sigma}\) is the conductivity term, and \(\mathbf{m}_{\sigma}^{0}\) is the prior model, which is the initial model set in this paper. \(\lambda_{\varepsilon}\) and \(\lambda_{\sigma}\), respectively, represent the trade-off factors of the model parameter constraint terms related to the permittivity and conductivity, and \(TV\) represents the modified total variation regularization operator based on the L1 norm [34]:  
\[TV(\mathbf{m}_{\varepsilon})=\|\
abla\mathbf{m}_{\varepsilon}\|_{1}=\sum_{ \varepsilon=1}^{NE}\sqrt{\left(\
abla_{x}m_{\varepsilon}\right)^{2}+\left( \
abla_{y}m_{\varepsilon}\right)^{2}+\left(\
abla_{z}m_{\varepsilon}\right)^ {2}+\beta_{\varepsilon}^{2}}. \tag{13}\]  
Among these terms, \(\beta_{\varepsilon}^{2}\) is usually set to a minimum value, ensuring that the function is differentiable when the partial derivative is zero; it is \(1\times 10^{-20}\) in this paper. The gradient calculation formula is expressed as  
\[\mathbf{g}_{m}(\mathbf{m}_{\varepsilon})=\
abla TV(\mathbf{m}_{\varepsilon} )=-\
abla\cdot\left(\frac{\
abla\mathbf{m}_{\varepsilon}}{\sqrt{\left( \
abla\mathbf{m}_{\varepsilon}\right)^{2}+\beta_{\varepsilon}^{2}}}\right), \tag{14}\]\(MSG\) represents the conductivity term constrained by MSG regularization. The objective function can be expressed as a pseudo-quadratic function of the model parameters in the inverse problem [34; 41]:  
\[MSG\left(\mathbf{m}_{\sigma}-\mathbf{m}_{\sigma}^{0}\right)=\left\|\
abla\frac{ \mathbf{m}_{\sigma}-\mathbf{m}_{\sigma}^{0}}{\sqrt{\left(\mathbf{m}_{\sigma}- \mathbf{m}_{\sigma}^{0}\right)^{2}+\beta_{\sigma}^{2}}}\right\|_{2}^{2}=\left\| \mathbf{W}_{m}\left(\mathbf{m}-\mathbf{m}^{0}\right)\right\|_{2}^{2}, \tag{15}\]  
where \(\beta_{\sigma}^{2}\) is the focusing parameter of conductivity. When using the MSG method, the selection dependence on the focusing parameter is relatively weak, but empirical optimization selection is still required. Meanwhile, \(\mathbf{W}_{m}\) can be expressed as the product of the first-order difference matrix \(\mathbf{W}_{\bigtriangledown}\) and the coefficient matrix \(\mathbf{W}\):  
\[\mathbf{W}_{m}=\mathbf{W}_{\bigtriangledown}\cdot\mathbf{W}, \tag{16}\]  
where the first-order difference matrix \(\mathbf{W}_{\bigtriangledown}\) in the one-dimensional parameter space can be expressed as follows. In the three-dimensional case, it is only necessary to calculate the gradients in the \(x\), \(y\), and \(z\) directions separately:  
\[\mathbf{W}_{\bigtriangledown}= \begin{bmatrix}0&0&0&\dots&0&0\\ -1&1&0&\dots&0&0\\ 0&-1&1&\dots&0&0\\ 0&0&-1&\ddots&\vdots&0\\ \vdots&\vdots&\vdots&\dots&1&\vdots\\ 0&0&0&\dots&-1&1\end{bmatrix} . \tag{17}\]  
The coefficient matrix \(\mathbf{W}\) is  
\[\mathbf{W}= \begin{bmatrix}\frac{1}{\sqrt{\left(m_{1}-m_{1}^{0}\right)^{2}+ \beta^{2}}}&\dots&0\\ \vdots&\ddots&\vdots\\ 0&\dots&\frac{1}{\sqrt{\left(m_{NE}-m_{NE}^{0}\right)^{2}+\beta^{2}}}\end{bmatrix} . \tag{18}\]  
First, we suppose that the model parameter of the current iteration is \(\mathbf{m}_{k}\).