Polarimetric sar image simulation method based on complex terrain conditions of forest

Abstract

The application provides a polarization SAR image simulation method for forests under complex terrain conditions, comprising the following steps: (1) setting parameters of a polarization SAR imaging process, including radar imaging parameters, forest parameters and terrain parameters; (2) generating a corresponding terrain model according to the terrain parameters; (3) constructing a forest model according to the forest parameters; (4) calculating a polarization SAR image according to the forest parameter model, and finally simulating the polarization SAR image of the forest under the complex terrain conditions. The method can improve the simulation of the polarization SAR image of the forest under the complex terrain conditions, can be used for understanding the remote sensing mechanism, testing new applications and designing the inversion algorithm of the forest structure, and can further explore the influence of the terrain on the polarization SAR image imaging process.

Description

Technical Field

[0001] This invention relates to the field of SAR image simulation, and in particular to a method for simulating polarimetric SAR images of forests under complex terrain conditions. Background Technology

[0002] Polarimetric Synthetic Aperture Radar (POLSAR), also known as polarimetric SAR, is an active microwave remote sensing technology that has shown increasing advantages in forestry applications such as forest identification, extraction, and inversion of forest biological parameters. The amplitude and phase information of polarimetric SAR images are crucial for extracting vegetation biological parameters. Polarimetric SAR image simulation is a prerequisite and a necessary and critical processing step for conducting forest identification, extraction, and inversion of forest biological parameters using polarimetric SAR. In areas with complex terrain, changes in topography can easily cause variations in the amplitude values ​​and polarimetric azimuth angles of polarimetric SAR images. To fully understand the patterns of these variations under complex terrain conditions, it is necessary to develop simulators for polarimetric SAR images of forests under complex terrain conditions.

[0003] There are three typical SAR image simulation methods. The first is feature-based SAR image simulation, mainly used to simulate the geometric and radiometric features of SAR images. This method directly maps the target's scattered field information to the image domain to obtain a SAR grayscale image. While simple in principle and fast in simulation, this method can only simulate the backscattering intensity information of the object under study and lacks the ability to simulate phase information. The simulated image only has amplitude information and no phase information, resulting in low fidelity. The second method is SAR image-based simulation, which combines incoherent simulation with information extracted from actual SAR images to achieve a simulated image of a different imaging configuration. This method is fast and accurate, but it requires acquiring a real SAR image first, which is difficult to obtain. Furthermore, the viewing angle of the target in the simulated image cannot differ too much from the viewing angle in the original image, thus limiting the method's effectiveness.

[0004] The third method is SAR image simulation based on echo signals. This method recreates the operation of the SAR system, uses a coherent scattering model to simulate the process, and finally obtains the SAR image by simulating the original echo signal. This method focuses on simulating the imaging process, achieving high accuracy. Furthermore, the use of a coherent scattering model allows for the acquisition of phase information under various polarization modes, resulting in realistic simulated images. This method has become a hot topic in the development of polarimetric SAR image simulation.

[0005] In echo-based polarimetric SAR image simulation methods, high-fidelity forest polarimetric SAR images have been simulated by establishing a coherent forest model. However, during the analysis of polarimetric SAR images, it has been gradually discovered that terrain has a significant impact on the imaging process and subsequent applications. Existing simulators cannot simulate forest polarimetric SAR images under complex terrain conditions. Therefore, improving and optimizing forest polarimetric SAR image simulation methods to explore the influence of terrain slope on polarimetric SAR images has become an urgent problem to be solved. Summary of the Invention

[0006] This invention provides a polarimetric SAR image simulation method for forests under complex terrain conditions. The aim is to provide sufficiently realistic simulation test data under complex terrain conditions, overcoming the limitations of available data to achieve understanding of remote sensing mechanisms, testing new forest monitoring applications and designs, and forest structure inversion algorithms. The technical solution is as follows:

[0007] A method for simulating forests using polarimetric SAR images under complex terrain conditions includes the following steps:

[0008] S1: Set the parameters for the polarimetric SAR imaging process, including radar imaging parameters, forest parameters, and terrain parameters;

[0009] S2: Generate the corresponding terrain model based on terrain parameters;

[0010] S3: Construct a forest model based on forest parameters;

[0011] S4: Based on the constructed terrain model and forest model, calculate the polarimetric SAR image and finally output the polarimetric SAR image of the forest under complex terrain conditions.

[0012] Furthermore, in step S1, the parameters for setting the polarimetric SAR imaging process include:

[0013] Radar imaging parameters: platform height, radar incident angle, wavelength center frequency, azimuth resolution and range resolution;

[0014] Terrain parameters: azimuth slope, distance slope, surface roughness, and surface dryness; when inputting terrain parameters, it is possible to input forward and reverse slope parameters. When the input terrain slope parameter is greater than 0, the terrain plane tilts in the direction away from the sensor, forming a forward slope; when the input terrain slope parameter is less than 0, the terrain plane tilts in the direction towards the sensor, forming a reverse slope.

[0015] Forest parameters: tree species type, average tree height, forest area, and forest density.

[0016] Furthermore, when generating terrain models, it is possible to generate terrain models with positive slopes as well as terrain models with negative slopes.

[0017] The terrain surface constructed using discrete Fourier transform is as follows:

[0018]

[0019] Where -N≤m, n≤N, and N is the number of terms in the Fourier expansion. h mn Let L be the discrete coefficients of the random surface function, x be the image azimuth value, y be the image range value, i be the imaginary unit, and L be the discrete coefficients of the random surface function. x L represents the azimuth length of the image. y The image distance length;

[0020] When there is a directional slope s on the terrain surface x and distance slope s y This makes the surface height for any point in the SAR image plane equal to: xs x +ys y Then the normal vector of the terrain surface is:

[0021]

[0022] The normal vector of the terrain surface. and These represent the unit vectors along the x, y, and z axes in the ground-distance coordinate system; the normal vector of the terrain surface when the ground is tilted. unit vector in the z-axis direction The following relationship exists:

[0023]

[0024] θ t Let represent the angle between the terrain surface normal vector and the z-axis; therefore, the final height of the ground at the azimuth x-axis and the range y-axis is:

[0025] g(x,y)=xs x +ys y +h(x,y)cosθ t (4)

[0026]

[0027] Terrain in x g The direction has a slope of s x , in y g The direction has a slope of s y When s yWhen s = 0, it is a flat terrain; when s = 0, it is a flat terrain. y When >0, the terrain plane is y g The axis tilts in the positive direction, forming a positive slope; when s y <0, terrain plane in y direction g The axis tilts in the negative direction, forming a reverse slope; in x g Keeping the direction unchanged, by changing y g The slope in a certain direction creates complex terrain.

[0028] When the terrain slope changes, the length of the overlay and shadow range will change, thus affecting the geometry of the terrain model in step S2.

[0029] The formula for calculating the length of the shadow area of ​​the positive slope is:

[0030]

[0031] The formula for calculating the length of the overlapping area of ​​the positive slope is:

[0032]

[0033] The formula for calculating the length of the shaded area of ​​the reverse slope is:

[0034]

[0035] The formula for calculating the length of the overlapping area of ​​the reverse slope is:

[0036]

[0037] Where S is the length of the shadow area, L is the length of the overlapping area, h is the forest parameter input in step S1: average tree height, θ is the radar incident angle of the radar imaging parameter, and α is the range slope of the terrain parameter.

[0038] Furthermore, in step S3, based on the forest parameters set in step S1: tree species type, average tree height, forest area, and forest density, the parameters required for forest model construction are calculated: tree height, canopy radius, number of trees, and tree location.

[0039] Tree height calculation includes the following steps:

[0040] The height of each tree will be extracted based on the normal distribution of the input average tree height, with the standard deviation of the tree height being 5% of the average tree height;

[0041] After calculating the tree height, the crown radius is calculated based on the allometric growth equation determined by the tree species and tree height.

[0042] Then, based on the maximum number of trees calculated from the canopy radius and the input forest area, the final number of trees is determined by combining this with the input forest density.

[0043] Furthermore, in step S3, the tree position coordinates refer to the location of the bottom of each tree trunk, which is composed of the azimuth and distance directions in the terrain model in step S2. In step S3, based on the number of trees, the height of the trees, and the crown radius, the Monte Carlo technique is used in the terrain model in step S2 to generate the tree position coordinates randomly and avoid collisions, so as to achieve a realistic tree position distribution. Then, the main trunk, primary and secondary branches of each tree, as well as the crown full of tertiary branches and leaves, are constructed in the terrain model in step S2 to complete the establishment of the forest model.

[0044] Furthermore, in step S4, forest polarimetric SAR images under both forward and reverse slopes are generated. Based on the forest model established in step S3, polarimetric synthetic aperture radar image calculations are performed on it, which is divided into the following stages:

[0045] (1) Calculate the effective dielectric constant of the canopy and understory: The scattered field obtained by the Foldy-Lax equation is converted into the corresponding forward scattering coefficient, and then integrated to determine the dielectric constant of the canopy and understory; after determining the dielectric constant of the canopy and understory, each element is assigned the same average dielectric constant;

[0046] (2) Constructing an attenuation grid: After determining the average dielectric constant, the attenuation function is sampled through a three-dimensional lattice. The scattering coefficient is specified to be the same as the polarization attenuation associated with the nearest grid point at the center of the scatterer. At the same time, attenuation calculation is performed according to the spatial location of the terrain model constructed in step S2 to provide a non-uniform attenuation grid.

[0047] (3) Calculation and output of polarimetric SAR images.

[0048] Furthermore, the calculation and output of polarimetric SAR images include the following steps:

[0049] S11: Calculate the backscattering coefficient and scattering matrix of the ground unaffected by vegetation: The terrain model generated in step S2 is divided into several small planes. Each plane has a different orientation and area depending on its location on the terrain model. The roughness of each plane is determined based on the surface roughness and Gaussian roughness model input in step S1. The backscattering coefficient of each small plane is then calculated using the Bragg scattering model. Finally, for each small plane, a speckle phase is assigned from a uniform random distribution on [0, 2π] to obtain its scattering matrix:

[0050]

[0051] S lS is the scattering matrix of the small planes divided by the terrain model. lhh S lhv S lvh and S lvv These are the complex scattering coefficients corresponding to the four polarization modes hh, hv, vh, and vv in the facets divided by the terrain model. and Let f be the amplitude term, f be a random number distributed exponentially, and A be the area of ​​the small plane. The backscattering coefficient for hh polarization is... For vv polarization, the backscattering coefficient, e iφ Let i be the phase term, i be the imaginary unit, and φ be a random angle generated according to a uniform distribution, where 0 ≤ φ ≤ 2π.

[0052] Then, the small plane is rotated around the radar to obtain the scattering matrix in the global coordinate system, with a rotation angle of β. n ,

[0053]

[0054] s x For azimuth slope, s y The range slope is θ, and the radar incident angle is θ.

[0055] Then, the ground scattering matrix under the global frame is calculated based on the scattering matrix of the small plane with respect to the rotation angle:

[0056]

[0057] S is the scattering matrix in the global coordinate system. hh S hv S vh and S vv These are the complex scattering coefficients for the four polarization modes hh, hv, vh, and vv in the global coordinate system, respectively, and S lhh S lhv S lvh and S lvv These are the complex scattering coefficients corresponding to the four polarization modes hh, hv, vh, and vv in the small planes divided by the terrain model, and β. n This is the rotation angle by which the small plane is transformed into the ground.

[0058] S12: Calculate the backscattering coefficients of the direct scattering of the short vegetation layer and the ground reflection: The forward scattering coefficients of the short vegetation layer and the ground reflection are calculated using the Foldy-Lax equation, and their scattering coefficients are scaled to produce the expected backscattering coefficients.

[0059] S13: Calculate the backscattering coefficients of the forest's direct scattering and ground reflection. Based on the location of each tree obtained in step S3, subdivide the unique distribution of the stem, primary and secondary branches created at each tree location into short cylindrical elements. Calculate the backscattering coefficients of these elements' direct scattering and ground reflection using the Rayleigh-Gans method. Then, calculate the backscattering coefficients of smaller canopy elements such as leaves and tertiary branches using the Foldy-Lax equation.

[0060] S14: Based on the calculated backscattering coefficients of the unvegetated ground, the direct scattering of the short vegetation layer, and the backscattering coefficients of the ground reflection, as well as the direct scattering and the backscattering coefficients of the forest, the effective scattering centers are calculated and the corresponding phases are assigned according to different scenarios. Their contributions are coherently added to the image, and finally, the simulated polarimetric SAR image is output.

[0061] The proposed method for simulating forests under complex terrain conditions using polarimetric SAR images can be used to study the impact of terrain on the imaging process of polarimetric SAR images, verify the feasibility and accuracy of the polarimetric interferometric SAR algorithm for retrieving forest height, and further be used to understand remote sensing mechanisms, test new forest monitoring applications, and design inversion algorithms for forest structure. Attached Figure Description

[0062] Figure 1 It generates an image diagram of the corresponding terrain based on the parameters;

[0063] Figure 2 This is a schematic diagram of the positive slope, where α represents the positive slope value;

[0064] Figure 3 This is a schematic diagram of the reverse slope, where α represents the reverse slope value;

[0065] Figure 4 It is the geometric composition of the terrain model;

[0066] Figure 5 This is a schematic diagram of the shadow length in the reverse slope diagram, where S is the shadow area length, α represents the reverse slope value, h is the tree height, and θ is the radar incident angle;

[0067] Figure 6 This is a schematic diagram of the reverse slope and the overlapping length, where L is the overlapping area length, α represents the reverse slope value, h is the tree height, and θ is the radar incident angle;

[0068] Figure 7The images are backscattering coefficient maps of simulated forest polarimetric SAR images with positive slope, where the radar imaging parameters and forest parameters are exactly the same, only the range slope parameter in the terrain parameters is different. (a), (b) and (c) are images with range slopes of 5°, 10° and 15°, respectively.

[0069] Figure 8 This is a backscattering coefficient map of a simulated forest polarimetric SAR image with reverse slope. The radar imaging parameters and forest parameters are exactly the same, except for the range slope parameter in the terrain parameters. (a), (b), and (c) are images with range slopes of -5°, -10°, and -15°, respectively.

[0070] Figure 9 This is a flowchart illustrating the present invention. Detailed Implementation

[0071] The technical solution of the present invention will be further described below with reference to the accompanying drawings, but it is not limited thereto. Any modifications or equivalent substitutions to the technical solution of the present invention that do not depart from the spirit and scope of the technical solution of the present invention should be covered within the protection scope of the present invention.

[0072] like Figure 1 As shown, the method for simulating polarimetric SAR images of forests under complex terrain conditions is designed for simulating polarimetric interferometric synthetic aperture radar images, such as... Figure 9 As shown, it includes the following steps:

[0073] S1: Set the parameters for the polarimetric synthetic aperture radar imaging process, including radar imaging parameters, forest parameters, and terrain parameters;

[0074] S2: Generate the corresponding terrain model based on terrain parameters;

[0075] S3: Construct a forest model based on forest parameters;

[0076] S4: Based on the constructed terrain model and forest model, calculate the polarimetric SAR image and finally simulate the polarimetric SAR image of the forest under complex terrain conditions.

[0077] In step S1, the parameters for the polarimetric SAR imaging process are set. The input parameters include the following three types:

[0078] Radar imaging parameters: Platform height: The height of the simulated radar sensor above the ground, in meters; Radar incident angle: The angle between the radar beam and a straight line on the vertical surface, in degrees; Wavelength center frequency: The frequency corresponding to the simulated radar band, in gigahertz; Azimuth resolution: The resolution along the flight path, in meters; Range resolution: The resolution perpendicular to the flight path, in meters.

[0079] Terrain parameters: Azimuth slope: Terrain slope along the flight path, in percentage; Distance slope: Terrain slope perpendicular to the flight path, in percentage; Surface roughness: Simulates the surface roughness of the ground, ranging from smooth to rough (0-10); and Ground dryness: Simulates the dryness and wetness of the ground, ranging from dry to wet (0-10).

[0080] Forest parameters: Tree species type: includes shrubs, coniferous forests and broad-leaved forests, represented by numbers 0, 1 and 2 respectively; Average tree height: the average height of trees in the simulated forest, in meters; Forest area: the area of ​​the simulated forest, in hectares; Forest density: the number of trees per unit area, in trees per hectare.

[0081] When inputting terrain parameters, not only positive slope parameters but also negative slope parameters can be entered. When the input terrain slope parameter is greater than 0, the terrain plane tilts away from the sensor, forming a positive slope; when the input terrain slope parameter is less than 0, the terrain plane tilts towards the sensor, forming a negative slope.

[0082] like Figure 1 As shown, the corresponding terrain is generated based on the parameters. In the synthetic aperture radar simulation, the imaging radar travels in the positive x-axis direction.

[0083] The terrain surface is constructed using discrete Fourier transform as follows:

[0084]

[0085] Where -N≤m, n≤N, and N is the number of terms in the Fourier expansion. h mn Let L be the discrete coefficients of the random surface function, x be the image azimuth value, y be the image range value, i be the imaginary unit, and L be the discrete coefficients of the random surface function. x L represents the azimuth length of the image. y The distance in the image is the length.

[0086]

[0087] l is the large-scale surface correlation length, and i and j are loop variables where 0 ≤ i, j ≤ N. Δx is the azimuth interval, Δy is the range interval, L x L represents the azimuth length of the image. y The distance in the image is the length.

[0088] When there is a directional slope s on the terrain surface x and distance slope s y This makes the surface height for any point in the SAR image plane equal to: xs x +ysy Then the normal vector of the terrain surface is:

[0089]

[0090] The normal vector of the terrain surface. and These represent the unit vectors along the x, y, and z axes in the ground coordinate system. The terrain surface normal vector is defined when the ground is tilted. unit vector in the z-axis direction The following relationship exists:

[0091]

[0092] θ t Let represent the angle between the terrain surface normal vector and the z-axis; therefore, the final height of the ground at the azimuth x-axis and the range y-axis is:

[0093] g(x,y)=xs x +ys y +h(x,y)cosθ t (5)

[0094]

[0095] Terrain in x g The direction has a slope of s x , in y g The direction has a slope of s y h(x,y) is the terrain surface function. This patent addresses x... g Keeping the direction unchanged, by changing y g The slope in a certain direction creates complex terrain.

[0096] Furthermore, in step S2, not only can things like... Figure 2 The terrain model with the positive slope shown can also generate models such as Figure 3 The terrain model shown is based on the reverse slope. Figure 4 As shown, changes in terrain slope simultaneously alter the length of overlay and shadow extents, thus affecting the geometry of the terrain model.

[0097] The formula for calculating the length of the shadow area of ​​the positive slope is:

[0098]

[0099] The formula for calculating the length of the overlapping area of ​​the positive slope is:

[0100]

[0101] Taking into full account the principles of radar imaging and the formation of shadows and overlays during image imaging, the formula for calculating the area and length of the shadow due to reverse slope is as follows:

[0102]

[0103] The formula for calculating the length of the overlapping area of ​​the reverse slope is:

[0104]

[0105] Among them, such as Figure 5 As shown, S is the length of the shaded area, α represents the reverse slope value, h is the tree height, and θ is the radar incident angle. Figure 6 As shown, L is the overlap area length, α represents the reverse slope value, h is the tree height, and θ is the radar incident angle.

[0106] Furthermore, in step S3, the input parameters for constructing the forest model include: tree species type, average tree height, forest area, and forest density, which are determined based on the forest parameters input in step S1; the forest parameters that need to be calculated include: tree height, canopy radius, number of trees, and tree location.

[0107] Based on the number of trees, their height, and canopy radius, Monte Carlo techniques are used in the terrain model of step S2 to randomly generate tree position coordinates without collisions, achieving a realistic tree location distribution. Specifically, this includes the following steps:

[0108] (1) Tree height calculation: The height of each tree will be extracted based on the normal distribution of the input average tree height, and the standard deviation of the tree height is 5% of the average tree height.

[0109] (2) After calculating the tree height, the crown radius is calculated based on the allometric growth equation determined by the tree species and tree height.

[0110] (3) The maximum number of trees is calculated based on the canopy radius and the input forest area, and the final number of trees is determined by combining the input forest density.

[0111] (4) Tree position generation: The position coordinates of the tree refer to the location of the bottom of each tree trunk, which is composed of the azimuth and distance values ​​in the terrain model in step S2; in step S3, based on the number of trees, the height of the trees and the radius of the tree crown, the Monte Carlo technique is used in the terrain model in step S2 to generate the tree position coordinates randomly and avoid collisions, so as to achieve a realistic tree position distribution.

[0112] Then, in the terrain model of step S2, the trunk of each tree is constructed. First-order branches are then constructed on the trunk, second-order branches are constructed on the first-order branches, and finally, a canopy filled with tertiary branches and leaves is constructed on the second-order branches, thus completing the forest model. The details are as follows:

[0113] After determining the location of the trees, a tree model is constructed. The trunk and branches at all levels are represented in the computer as curved conical cylinders, and the equation describing their shape is:

[0114]

[0115] b0 represents the starting point of the branch; l is the branch length; t is the branch position parameter value and 0≤t≤1, which is zero at the beginning of the branch and 1 at the end of the branch; unit vector. The initial direction of the branch is defined; in the absence of shape change, c(t) is the branch curvature when the branch is straight. It is the direction of curvature; d p It is the curvature factor, and 0 ≤ d p <1, used to control The direction of curvature. The canopy is filled with leaves and tertiary branches. The relevant parameters for constructing the canopy model include: size, dielectric constant, orientation, and volume fraction. These parameters are defined according to the tree species. The corresponding canopy type can be defined by inputting the tree species type in step S1.

[0116] The construction of the forest model is completed by determining the forest parameters in step S1, generating the ground model in step S2, calculating the relevant forest parameters in step S3, and constructing the vegetation model.

[0117] Furthermore, after completing the construction of the forest model, in step S4, the polarimetric synthetic aperture radar image is calculated and output:

[0118] (1) Calculate the effective dielectric constant of the tree canopy and understory;

[0119] The effective dielectric constant of the canopy and understory is calculated using the Foldy-Lax equation by averaging the forward scattering coefficients of its constituent elements. The Foldy-Lax equation is expressed as:

[0120] ψ=ψ inc +ψ s (12)

[0121] ψ represents the received scattered field, ψ inc ψ represents the scattered field of the incident wave. s It represents the sum of the scattered fields of all particles in the scattering object;

[0122]

[0123] j represents a single particle in the scattering object, and N represents the total number of particles in the scattering field. This represents the scattering field of particle j;

[0124]

[0125] T represents the excitation field of particle j. j G represents the transition operator for particle j, and G0 represents the Green's function;

[0126]

[0127] ψ inc T represents the scattered field of the incident wave, l represents a particle different from j, and T l The transition operator for particle l. This represents the excitation field of particle l.

[0128] (2) Construct a decay grid;

[0129] After determining the average dielectric constant, the attenuation function is sampled through a three-dimensional lattice. The scattering coefficient is specified to be the same as the polarization attenuation associated with the nearest grid point at the center of the scatterer. At the same time, attenuation calculation is performed according to the spatial location of the terrain model constructed in step S2 to provide a non-uniform attenuation grid.

[0130] (3) Calculation and output of polarimetric SAR images: Based on the calculated backscattering coefficient of the ground scattering without vegetation cover, the direct scattering amount of the short vegetation layer and the backscattering coefficient after secondary scattering by ground reflection, as well as the direct scattering amount of the forest and the backscattering coefficient after secondary scattering by ground reflection.

[0131] (a) Calculate the backscattering coefficient and scattering matrix of the ground that is not obscured by vegetation.

[0132] The terrain model generated in step S2 is divided into several small planes. Each plane has a different orientation and area depending on its location on the terrain model. The roughness of each plane is determined based on the surface roughness and Gaussian roughness model input in step S1. Then, the backscattering coefficient of each small plane is calculated using the Bragg scattering model. Finally, for each small plane, a speckle phase is assigned from a uniform random distribution on [0, 2π] to obtain its scattering matrix:

[0133]

[0134] S l S is the scattering matrix of the small planes divided by the terrain model. lhh S lhv S lvhand S lvv These are the complex scattering coefficients corresponding to the four polarization modes hh, hv, vh, and vv in the facets divided by the terrain model. and Let f be the amplitude term, f be a random number distributed exponentially, and A be the area of ​​the small plane. The backscattering coefficient for hh polarization is... For vv polarization, the backscattering coefficient, e iφ Let i be the phase term, i be the imaginary unit, and φ be a random angle generated according to a uniform distribution, where 0 ≤ φ ≤ 2π.

[0135] Then, the small plane is rotated around the radar to obtain the scattering matrix in the global coordinate system, with a rotation angle of β. n and

[0136]

[0137] s x For azimuth slope, s y θ represents the range slope, and θ is the radar incident angle.

[0138] Then, the ground scattering matrix under the global frame is calculated based on the scattering matrix of the small plane with respect to the rotation angle:

[0139]

[0140] S is the scattering matrix in the global coordinate system. hh S hv S vh and S vv These are the complex scattering coefficients for the four polarization modes hh, hv, vh, and vv in the global coordinate system, respectively, and S lhh S lhv S lvh and S lvv These are the complex scattering coefficients corresponding to the four polarization modes hh, hv, vh, and vv under the small planes divided by the terrain model, and β. n This is the rotation angle by which the small plane is transformed into the ground.

[0141] (b) Calculate the backscattering coefficients of the direct scattering from the short vegetation layer and the ground reflectance. The backscattering coefficients of the short vegetation layer and the ground reflectance were calculated using the Foldy-Lax equations, and their scattering coefficients were scaled to produce the expected backscattering coefficients.

[0142] (c) Calculate the backscattering coefficients of the forest direct scattering and ground reflection. Based on the location of each tree obtained in step S3, subdivide the unique distribution of the stem, primary and secondary branches created at each tree location into short cylindrical elements. Calculate the backscattering coefficients of these elements direct scattering and ground reflection using the Rayleigh-Gans method. Then calculate the backscattering coefficients of smaller canopy elements such as leaves and tertiary branches using the Foldy-Lax equation.

[0143] (d) Based on the calculated backscattering coefficients of ground scattering without vegetation obstruction, direct scattering from short vegetation layers, and ground reflection, as well as the backscattering coefficients of direct scattering from forests and ground reflection, the backscattering coefficients corresponding to hh polarization are used. Calculate the amplitude term for the scene corresponding to the polarization mode. Then, based on different scenarios, the angle φ of the effective scattering center allocation is calculated to obtain the corresponding phase term e. iφ As shown in formula (13). Then, according to formula (18), the scattering matrix S and complex scattering coefficient S of the corresponding scene are obtained. hh S hv S vh and S vv The contributions are then accumulated into an image accumulator. The image accumulator is a function of image position and complex scattering coefficients. By accumulating the complex scattering coefficients of different scenes at the same position, the complex scattering coefficients of all positions can be obtained to represent the complex scattering coefficients and scattering matrix of the entire image.

[0144] (4) Output simulated polarimetric SAR images, such as Figure 7 and Figure 8 As shown.

[0145] The proposed method for simulating forests under complex terrain conditions using polarimetric SAR images can be used to study the impact of terrain on the imaging process of polarimetric SAR images, verify the feasibility and accuracy of the polarimetric interferometric SAR algorithm for retrieving forest height, and further be used to understand remote sensing mechanisms, test new forest monitoring applications, and design inversion algorithms for forest structure.

Claims

1. A method for simulating forests under complex terrain conditions using polarimetric SAR images, comprising the following steps: S1: Set the parameters for the polarimetric SAR imaging process, including radar imaging parameters, forest parameters, and terrain parameters; S2: Generate the corresponding terrain model based on terrain parameters; S3: Construct a forest model based on forest parameters; S4: Based on the constructed terrain model and forest model, calculate the polarimetric SAR image and finally output the polarimetric SAR image of the forest under complex terrain conditions. In step S1, the parameters for setting the polarimetric SAR imaging process include: Radar imaging parameters: platform height, radar incident angle, wavelength center frequency, azimuth resolution and range resolution; Terrain parameters: azimuth slope, distance slope, surface roughness, and surface dryness; when inputting terrain parameters, it is possible to input forward and reverse slope parameters. When the input terrain slope parameter is greater than 0, the terrain plane tilts in the direction away from the sensor, forming a forward slope; when the input terrain slope parameter is less than 0, the terrain plane tilts in the direction towards the sensor, forming a reverse slope. Forest parameters: tree species type, average tree height, forest area, and forest density; When generating terrain models, it can generate terrain models with positive slopes as well as terrain models with negative slopes. The terrain surface constructed using discrete Fourier transform is as follows: ； in , Let be the number of terms in the Fourier expansion. , , The discrete coefficients of the random surface function, This represents the azimuth value of the image. This represents the image distance dimension. The imaginary unit, The image azimuth length. The image distance length; When there is a directional slope on the terrain surface and distance slope This makes the surface height for any point in the SAR image plane equal to: Then the normal vector of the terrain surface is: ； The normal vector of the terrain surface. , and These represent the unit vectors along the x, y, and z axes in the SAR ground distance coordinate system; the terrain surface normal vector when the ground is tilted. unit vector in the z-axis direction The following relationship exists: ； Let represent the angle between the terrain surface normal vector and the z-axis; therefore, the final height of the ground at the azimuth x-axis and the range y-axis is: ； ； The terrain is Slope in direction ,exist Slope in direction ;when At that time, it was flat terrain; when Time terrain plane The axis tilts in the positive direction, forming a positive slope; when Time terrain plane The axis tilts in the negative direction, forming a reverse slope; in Keeping the direction unchanged, by changing The slope in a certain direction creates complex terrain.

2. The method for simulating forests under complex terrain conditions using polarimetric SAR images according to claim 1, characterized in that: When the terrain slope changes, the length of the overlay and shadow range will change, thus affecting the geometry of the terrain model in step S2. The formula for calculating the length of the shadow area of ​​the positive slope is: ； The formula for calculating the length of the overlapping area of ​​the positive slope is: ； The formula for calculating the length of the shaded area of ​​the reverse slope is: ； The formula for calculating the length of the overlapping area of ​​the reverse slope is: ； in The length of the shaded area. The overlapping area is the length. The forest parameters input in step S1 are: average tree height, The radar incident angle is a radar imaging parameter. The distance slope is a terrain parameter.

3. The method for simulating forests under complex terrain conditions using polarimetric SAR images according to claim 1, characterized in that: In step S3, based on the forest parameters set in step S1: tree species type, average tree height, forest area, and forest density, the parameters required for forest model construction are calculated: tree height, canopy radius, number of trees, and tree location.

4. The method for simulating forests under complex terrain conditions according to claim 3, characterized in that: Tree height calculation includes the following steps: The height of each tree will be extracted based on a normal distribution of the input average tree height, with a standard deviation of 5% of the average tree height. After calculating the tree height, the crown radius is calculated based on the allometric growth equation determined by the tree species and tree height. Then, based on the maximum number of trees calculated from the canopy radius and the input forest area, the final number of trees is determined by combining this with the input forest density.

5. The method for simulating forests under complex terrain conditions according to claim 3, characterized in that: In step S3, the tree position coordinates refer to the location of the bottom of each tree trunk, which is composed of the azimuth and distance directions in the terrain model in step S2. In step S3, based on the number of trees, the height of the trees, and the crown radius, the Monte Carlo technique is used in the terrain model in step S2 to generate the tree position coordinates randomly and avoid collisions, so as to achieve a realistic tree position distribution. Then, the main trunk, primary and secondary branches of each tree, as well as the crown full of tertiary branches and leaves, are constructed in the terrain model in step S2 to complete the establishment of the forest model.

6. The method for simulating forests under complex terrain conditions using polarimetric SAR images according to claim 5, characterized in that: In step S4, polarimetric SAR images of the forest under both forward and reverse slopes are generated. Based on the forest model established in step S3, polarimetric synthetic aperture radar image calculations are performed on it, which is divided into the following stages: (1) Calculate the effective dielectric constant of the canopy and understory: The scattered field obtained by the Foldy-Lax equation is converted into the corresponding forward scattering coefficient, and then integrated to determine the dielectric constant of the canopy and understory; After determining the dielectric constant of the canopy and understory, each element is assigned the same average dielectric constant; (2) Constructing an attenuation grid: After determining the average dielectric constant, the attenuation function is sampled through a three-dimensional lattice. The scattering coefficient is specified to be the same as the polarization attenuation associated with the nearest grid point at the center of the scatterer. At the same time, attenuation calculation is performed according to the spatial location of the terrain model constructed in step S2 to provide a non-uniform attenuation grid. (3) Calculation and output of polarimetric SAR images.

7. The method for simulating forests under complex terrain conditions according to claim 6, characterized in that: The calculation and output of polarimetric SAR images include the following steps: S11: Calculate the backscattering coefficient and scattering matrix of the ground unaffected by vegetation: The terrain model generated in step S2 is divided into several small planes. Each plane has a different orientation and area depending on its location on the terrain model. The roughness of each plane is determined based on the surface roughness and Gaussian roughness model input in step S1. The backscattering coefficient of each small plane is then calculated using the Bragg scattering model. Finally, for each small plane, a speckle phase is assigned from a uniform random distribution on [0, 2π] to obtain its scattering matrix: ； The scattering matrix of the small planes divided by the terrain model. , , and These are small planes divided by the terrain model. , , and Complex scattering coefficients corresponding to the four polarization modes and For the amplitude term, For random numbers based on an exponential distribution, Let the area of ​​the smaller plane be . for The backscattering coefficient of polarization, for The backscattering coefficient of polarization, For phase terms, The imaginary unit, For random angles generated based on a uniform distribution and ; Then, the small plane is rotated around the radar to obtain the scattering matrix in the global coordinate system, with the rotation angle being... , ； This refers to the azimuth slope. For distance slope, The radar incident angle; Then, the ground scattering matrix under the global frame is calculated based on the scattering matrix of the small plane with respect to the rotation angle: ； The scattering matrix in the global coordinate system. , , and In the global coordinate system , , and Complex scattering coefficients for four polarization modes , , and These are small planes divided by the terrain model. , , and Complex scattering coefficients corresponding to the four polarization modes The rotation angle for the small plane to transform into the ground; S12: Calculate the backscattering coefficients of the direct scattering of the short vegetation layer and the ground reflection: The forward scattering coefficients of the short vegetation layer and the ground reflection are calculated using the Foldy-Lax equation, and their scattering coefficients are scaled to produce the expected backscattering coefficients. S13: Calculate the backscattering coefficients of the forest's direct scattering and ground reflection. Based on the location of each tree obtained in step S3, subdivide the unique distribution of the stem, primary and secondary branches created for each tree location into short cylindrical elements. Calculate the backscattering coefficients of these elements' direct scattering and ground reflection using the Rayleigh-Gans method. Then, calculate the backscattering coefficients of the leaves and smaller canopy elements of the tertiary branches using the Foldy-Lax equation. S14: Based on the calculated backscattering coefficients of the unvegetated ground, the direct scattering of the short vegetation layer, and the backscattering coefficients of the ground reflection, as well as the direct scattering and the backscattering coefficients of the forest, the effective scattering centers are calculated and the corresponding phases are assigned according to different scenarios. Their contributions are coherently added to the image, and finally, the simulated polarimetric SAR image is output.

Images