Two-dimensional rough surface modeling

In the study of this paper, in order to improve the accuracy of the model, the two-dimensional rough surface was generated with full consideration of the actual conditions of the soil ground and integrating the Monte Carlo method [Reference Caflisch16] and the Gaussian spectral function [Reference Bourlier, Bergine and Saillard17].

The surface height function of a two-dimensional rough surface is

(11)$$f( {x, \;y} ) = \displaystyle{1 \over {L_xL_y}}\sum\limits_{m = {-}\infty }^\infty {\sum\limits_{n = {-}\infty }^\infty {b_{mn}} \exp \left({\displaystyle{{\,j2\pi mx} \over {L_x}}} \right)} \exp \left({\displaystyle{{\,j2\pi ny} \over {L_y}}} \right)$$

In the formula, L _x, L _y represents the length of the rough surface in two dimensions, and b _m,n is calculated as

(12)$$b_{mn} = 2\pi \sqrt {L_xL_yW( {k_{xm}, \;k_{yn}} ) } \cdot \alpha $$

(13)$$\alpha {\rm} = \left\{{\matrix{ {\displaystyle{{N( {0, \;1} ) + jN( {0, \;1} ) } \over {\sqrt 2 }}} \hfill & {m\ne 0, \;{{N_x} / {2, \;n\ne 0, \;{{N_y} / 2}}}} \hfill \cr {N( {0, \;1} ) } \hfill & {m = 0, \;{{N_x} / {2\;{\rm or}\;n = 0, \;{{N_y} / 2}}}} \hfill \cr } } \right.$$

N(0, 1) is the random sampling points on the normal distribution, N _x and N _y are the number of sampling points on the rough surface in two dimensions, W(k _xm, k _yn) represents

(14)$$W( k_{xm}, \;k_{ym}) = \displaystyle{{l_xl_yh^2} \over {4\pi }}\exp \left({-\displaystyle{{k_x^2 l_x^2 } \over 4}-\displaystyle{{k_y^2 l_y^2 } \over 4}} \right)$$

where h is the root-mean-square height of the Gaussian rough surface, l _x and l _y are the correlation lengths of the rough surface in two dimensions,k _xm and k _yn are the discrete points of the spatial frequencies in two dimensions, respectively.

(15)$$k_{xm} = \displaystyle{{2\pi m} \over {L_x}}, \;\quad k_{yn} = \displaystyle{{2\pi n} \over {L_y}}$$

The schematic diagram of the two-dimensional rough surfaces with different roughness generated by the method in this paper is shown in Fig. 2.

Figure 2. Rough surface model. (a) h = 0.1 m,   l _x = l _y = 0.5 m; (b) h = 0.3 m,   l _x = l _y = 1.0 m.

Conical incident wave

A finite truncation of the infinitely extended rough surface is required in the numerical calculation. In order to eliminate the effects such as reflection and edge bypassing caused by the abrupt truncation of the rough surface edge, a conical incident wave is used instead of plane wave incidence.

The conical wave requires attention to the beam width, which determines the irradiation width of the incident wave. When it is too large, the size of the rough surface needs to be increased accordingly to ensure the irradiated area; when it is too small, the rough surface and the target cannot be fully irradiated and a reasonable compound scattering field cannot be obtained. Reasonable calculation results and faster efficiency can be obtained when g and rough surface size L satisfy the following relationship

(16)$$g \ge \displaystyle{6 \over {{( {\cos \theta_i} ) }^{1.5}}}, \;L = 4g$$

For TE incident waves, the incident electromagnetic field is

(17)$$\eqalign{{\boldsymbol E}_{inc}( {\boldsymbol r} ) & = \int_{-\infty }^{ + \infty } {\int_{-\infty }^{ + \infty } {\exp ( {{\rm j}k_xx + {\rm j}k_yy-{\rm j}k_zz} ) } } \cr & \quad \cdot E( {k_x, \;k_y} ) {\boldsymbol h}_i( {-k_z} ) dk_xdk_y}$$

(18)$$\eqalign{{\boldsymbol H}_{inc}( {\boldsymbol r} ) & = {-}\displaystyle{1 \over {\eta _0}}\int_{-\infty }^{ + \infty } {\int_{-\infty }^{ + \infty } {\exp ( {{\rm j}k_xx + {\rm j}k_yy-{\rm j}k_zz} ) } } \cr & \quad \cdot E( {k_x, \;k_y} ) {\boldsymbol v}_i( {-k_z} ) dk_xdk_y}$$

where E(k _x, k _y) is the incident spectrum, and the horizontal polarization direction ${\boldsymbol h}_i$ and vertical polarization direction ${\boldsymbol v}_i$ are denoted as

(19)$${\boldsymbol h}_i( {-k_z} ) = \displaystyle{1 \over {k_\rho }}( {{\boldsymbol x}k_y-{\boldsymbol y}k_x} ) $$

(20)$${\boldsymbol v}_i( {-k_z} ) = \displaystyle{{k_z} \over {k_0k_\rho }}( {{\boldsymbol x}k_x + {\boldsymbol y}k_y} ) + \displaystyle{{k_\rho } \over {k_0}}{\boldsymbol z}$$

η ₀ is the free space wave impedance, k _x, k _y, k _z are the spatial spectral domains in the x, y, z directions, respectively, $k_\rho = \sqrt {k_x^2 + k_y^2 }$.

For TM incident waves, the incident electromagnetic field is

(21)$$\eqalign{{\boldsymbol E}_{inc}( {\boldsymbol r} ) & = \int_{-\infty }^{ + \infty } {\int_{-\infty }^{ + \infty } {\exp ( {\,jk_xx + jk_yy-jk_zz} ) } } \cr & \quad \cdot E( {k_x, \;k_y} ) {\boldsymbol v}_i( {-k_z} ) dk_xdk_y}$$

(22)$$\eqalign{{\boldsymbol H}_{inc}( {\boldsymbol r} ) & = {-}\displaystyle{1 \over {\eta _0}}\int_{-\infty }^{ + \infty } {\int_{-\infty }^{ + \infty } {\exp ( {\,jk_xx + jk_yy-jk_zz} ) } } \cr & \quad \cdot E( {k_x, \;k_y} ) {\boldsymbol h}_i( {-k_z} ) dk_xdk_y}$$

The angles of the incident waves are set to $\theta _i = 60^\circ$,$\varphi _i = 0^\circ$, and the two-dimensional normalized conical incident wave amplitudes are shown in Fig. 3. From the figure, it can be seen that the magnetic field is strongest at the center of the rough surface, and the farther away from the center, the weaker the magnetic field is, which slowly goes to zero at the edge. Therefore, the errors caused by the truncation of the rough surface can be avoided in the numerical simulation.

Figure 3. Normalized conical wave schematic.

Modified small slope approximation

Based on the fact that the soil rough surface is anisotropic, the classical SSA method with appropriate modifications is used in this paper.

The core of the SSA method to solve the electromagnetic scattering problem is to solve for the scattering amplitude T. In the scattering field calculation, the solution equation for T is as follows

(23)$$T( {\boldsymbol q}, \;{\boldsymbol k}) = {-}\int_S {e^{{-}i{\boldsymbol q}\cdot {\boldsymbol x}}\hat{{\boldsymbol n}}} \cdot \nabla P_{\boldsymbol k}( {\boldsymbol x}) dS$$

where the point on the scattering cross section S is labeled by the two-dimensional vector x, $P_{\boldsymbol k}( {\boldsymbol r})$ satisfying the boundary conditions and the Helmholtz equation.

(24)$$( \nabla ^2 + k^2) P_{\boldsymbol k}( {\boldsymbol r}) = 0$$

Its complete solution can be obtained by integrating:

(25)$$P_{\boldsymbol k}( {\boldsymbol r}) = e^{i{\boldsymbol k}\cdot {\boldsymbol r}}-\displaystyle{1 \over {4\pi }}\int\limits_S {G_k( {\boldsymbol r}-{\boldsymbol x}) } \hat{{\boldsymbol n}}\cdot \nabla P_{\boldsymbol k}( {\boldsymbol x}) dS$$

(26)$$G_k( {\boldsymbol r}-{\boldsymbol x}) = e^{ik\vert {{\boldsymbol r}-{\boldsymbol x}} \vert }/\vert {{\boldsymbol r}-{\boldsymbol x}} \vert $$

where r is the surface element vector and $G_k( {\boldsymbol r}-{\boldsymbol x})$ represents the Green's function.

According to the above equation, the scattering amplitude value T of the soil rough surface can be found, which in turn leads to σ _suf:

(27)$$\sigma _{suf} = \left\langle {T{\rm \ast }T} \right\rangle {\rm /16}\pi ^ 2$$

Combine (23)(24) (25)(26)(27):

(28)$$\sigma _{suf} = \displaystyle{ 2 \over \pi }\left({\displaystyle{{k_zq_z} \over {Q_hQ_z}}} \right)^2\cdot I( \alpha ) $$

(29)$$\alpha = \displaystyle{{0.08\pi } \over {200[ \Gamma {( 7/4) }^2] }}\left[{{\left({\displaystyle{{u^2} \over g}} \right)}^{1/2}\displaystyle{{Q_z^2 } \over {Q_h^{3/2} }}} \right]$$

I(α) is the spectral integral of the rough surface and ${\boldsymbol k}-{\boldsymbol q} = ( Q_h, \;Q_z)$.Q _h and Q _z can be solved by the following equations.

(30)$$Q_h = k_0\cdot \sqrt {{\cos }^2\theta _{inc} + {\cos }^2\theta _{scat}-2\cos \theta _{inc}\cos \theta _{scat}\cos \phi _{bi}} $$

(31)$$Q_z = {-}k_0( {\sin \theta_{inc} + \sin \theta_{scat}} ) $$

It can be found that the soil rough surface scattering coefficient depends on the following integral.

(32)$$I( \alpha ) = \int_0^\infty {J_0( y) yB( y, \;\alpha ) {\rm d}y} $$

where

(33)$$J_0( y) = \sum\nolimits_{n = 0}^\infty {\displaystyle{{( {-}y^2/4) } \over {{( n!) }^2}}} ^n$$

(34)$$B( y, \;\alpha ) = \exp ( -\alpha y^{1.5}) $$

In the above equation, J ₀(y) is the first class zero-order Bessel function and B(y, α) is the calculated coefficient.

Combining the above equations, I(α) can be expanded into the following equation:

(35)$$I( \alpha ) = \displaystyle{2 \over 3}\alpha ^{{-}4/3}\sum\limits_{n = 0}^\infty {\displaystyle{{{( {-}n) }^n} \over {{( n!) }^2}}} \cdot \left({\displaystyle{1 \over {4\alpha^{4/3}}}} \right)^n\Gamma \left({\displaystyle{{4( n + 1) } \over 3}} \right)$$

The summation based on the traditional SSA requires manual truncation, and when α is large, the formula converges quickly and the error generated by manual truncation is small; when α is small, the formula converges slowly, which will lead to the error generated by truncation beyond the engineering range, which is uncontrollable in the actual calculation, so it is unrealistic to apply equation (35) for the solution. Therefore, modifications based on the traditional SSA method can be considered to speed up the calculation process and reduce the error.

Based on the above idea, we can find that Γ( − z) has a remainder −( − 1)ⁿ/n! at a non-negative integer n. Obviously, by applying the residue theorem, the series terms in I(α) can be turned into:

(36)$$\eqalign{& \sum\limits_{n = 0}^\infty {\displaystyle{{{( {-}n) }^n} \over {{( n!) }^2}}} \cdot \displaystyle{1 \over {4\alpha ^{4/3}}}\Gamma \left({\displaystyle{{4( n + 1) } \over 3}} \right)\cr = & -\displaystyle{1 \over {2\pi i}}\int_C {\displaystyle{{\Gamma ( {-}z) \Gamma \left[{\displaystyle{{4( n + 1) } \over 3}} \right]} \over {\Gamma ( z + 1) }}{\left({\displaystyle{1 \over {4\alpha^{4/3}}}} \right)}^zdz} } $$

Substituting (36) back into (35):

(37)$$\eqalign{I( \alpha ) & = {-}\displaystyle{2 \over \pi }\sum\limits_{n = 0}^N {\displaystyle{{{( -\alpha ) }^n} \over {( n!) }}} \cdot \left[{\Gamma \left({\displaystyle{4 \over 3}n + 1} \right)} \right]^2\cdot \sin \displaystyle{{3\pi n} \over 4}2^{3n/2} \cr & \pm R_N}$$

where

(38)$$R_N = \alpha ^{N + 0.5}A( N) B( N) $$

(39)$$A( N) = [ {2^{4.25 + 1.5N}e^{{-}1/2( N + 2.5) }{0.75}^{N + 1}} \cdot {x^{0.5( N + 3.5) }} ] /3\sqrt \pi $$

(40)$$\eqalign{& B( N) = \int_0^\infty {{( 1 + u^2) }^{{1 \over 4}( 3N + 3.5) }} \left[{{( 1 + \displaystyle{1 \over {4x}}) }^2 + u^2} \right]^{-( N + 1) /2}\cdot \cr & \exp \left[{-xu\left({2{\tan }^{{-}1}u-\displaystyle{4 \over 3}{\tan }^{{-}1}\displaystyle{u \over {1-1/x}}} \right)} \right]\cdot \displaystyle{{\cosh ( \pi xu) } \over {\cosh ( 4\pi xu/3) }}du} $$

In the above equation, B(N) decays exponentially, and the decay accelerates as u increases.

Thus, the modified SSA method can be used to calculate the scattering coefficient of soil rough surface to ensure the accuracy and calculation speed at the same time, which meets the engineering requirements.

Validation of algorithm

The main work of this subsection is to verify the correctness of the SSA algorithm used in this paper. The validation algorithm uses the MOM algorithm, which is convincingly accurate and precise in its calculation. In the validation, both SSA and MOM methods were used to calculate the backscattering coefficients for the same rough surface with the following settings of the rough surface parameters with reference to the arid region: dimension L _x × L _y = 10λ × 10λ, correlation length l _x × l _y = 0.8λ × 0.8λ, root mean square height h = 0.1λ, dielectric constant (13.61, − 0.03) and incident wave frequency f = 300MHz.

The calculation results of the two methods are plotted in Fig. 4, and the two curves match well, indicating that the algorithm is effective in calculating the rough surface of Gaussian media.

Figure 4. Electromagnetic Scattering Model.

Influence of different roughness

This subsection focuses on the effects of different soil roughness on the scattering characteristics in arid regions, and the backscattering coefficients were calculated for root-mean-square heights of 0.05, 0.10 and 0.15 m, respectively, and plotted in Fig. 5. The incident wave frequency is 300 MHz, the soil moisture content is 10.35%, and the HH and VV polarization results are as follows.

Figure 5. Validation of Algorithm.

Observing Fig. 6, in general, the scattering intensity decreases linearly with the increase of radar incidence angle, and the scattering intensity is maximum when the incident wave is perpendicular to the ground. From the calculation results of different roughness, as the soil roughness becomes larger, the backward scattering increases, and they are positively correlated. According to this property, the roughness of the region can be obtained by comparison. And the roughness of arid areas is often affected by desertification - the more severe the desertification, the more sand and gravel fill the gullies, the flatter the land surface, the lower the roughness, so to some extent it can reflect the degree of desertification in the area. In another way, the intensity of the backscattering coefficient in a certain area can be continuously detected to study whether the desertification in the area has been aggravated or weakened, and to provide a reference for the subsequent decision of regional management.

Figure 6. Curves with different root mean square. (a) HH polarization; (b) VV polarization.

Effect of soil moisture content

In arid areas where precipitation is scarce and evaporation is high, water scarcity is a major limiting factor for plant growth. However, arid areas have enough heat and they have the potential to become highly productive if they receive proper irrigation and fertilization. Therefore, this subsection calculates the scattering coefficients of soils with different moisture levels and discusses their scattering characteristics.

According to the actual measurement data in the southeastern part of Ejina Banner, Inner Mongolia, the soil moisture content varies greatly in different periods, with the highest value of 38.07% and the lowest value of 0.23%. In order to fully investigate the electromagnetic scattering characteristics in arid regions, the scattering results were calculated for two extreme water content conditions, as shown in Fig. 7.

Figure 7. Curves with different soil moisture content. (a) HH polarization; (b) VV polarization.

It was found that for the HH polarization results, the two curves could be clearly distinguished in the (0°, 85°) range, and their overlap was very small. While in the VV polarization mode, the scattering coefficients of soils with different moisture contents become closer and closer with increasing angles, and finally almost overlap. In general, the soil moisture content has a strong influence on the scattering coefficients, and the differences are greatest when the incident wave is perpendicular to a rough surface, which can be easily distinguished with the naked eye. It can be inferred from the simulated results of soil water content under two extreme conditions in the region that soil water content is an important factor affecting the electromagnetic scattering characteristics in arid regions. The variation of soil water content will lead to fluctuations of electromagnetic scattering intensity on the rough surface. The results of this paper can be combined with remote sensing techniques to invert soil water content based on radar scattering coefficients, which is important for soil water content studies such as irrigation control, ecological studies, vorticity covariance, and slope stabilization climate science.

The frequency effect

Different wavelengths of electromagnetic waves have different penetrating abilities and different sensitivities to rough surface parameters. Therefore, it is important to study the electromagnetic scattering characteristics of arid areas under different electromagnetic spectrum environments and find the electromagnetic “optimal solution” for environmental protection and ecological monitoring.

In this subsection of the study, the soil moisture content was 10.35% and the root mean square height of the rough surface was 0.05 m. The calculated results were plotted in Fig. 8.

Figure 8. Curves with different frequencies. (a) HH polarization; (b) VV polarization.

Through the curves we can find that, intuitively, the backscattering coefficient of the arid rough surface decreases continuously with the increase of the incident wave electromagnetic frequency. However, there are some differences between the results of HH polarization and VV polarization: ①For the VV polarization results, in the range of (0°, 90°), the three curves always do not intersect and are distinguished obviously; ②For the HH polarization results, in the angle of (0°, 48°), the three curves can be distinguished by the naked eye; in the range of (49°, 90°), the scattering coefficients of the three frequencies are extremely close to each other, and the magnitude cannot be determined by visual observation.

Therefore, when studying the effect of electromagnetic frequency on the scattering coefficient, better results can be obtained by choosing the incident wave with VV polarization.