CN106324571A - Fast Realization method for simulation 3D scene SAR radar echo based on forward method - Google Patents

Abstract

The invention discloses a fast Realization method for simulation 3D scene SAR radar echo based on forward method. The main idea of the fast Realization method comprises the following steps: acquiring digital elevation model (DEM) data within SAR radar beam range, wherein the data comprises K point targets; randomly distributing the K point targets on a 2D ground in the geographic coordinate system, and calculating DEM data within the SAR radar beam range after refinement processing; acquiring back scattering property distribution of the DEM data within the SAR radar beam range after refinement processing; utilizing an improved downward angle comparison method to determine shaded areas of Q small surface elements, calculating final actual back scattering coefficients of the Q small surface elements after shaded area determination, and then calculating SAR radar echo of the Q small surface elements. At this time, the whole 3D scene echo simulation is completed in a graphic processing unit (GPU), and finally an echo signal in the GPU is sent to a CPU processor to be output.

Description

Rapid implementation method of three-dimensional scene-simulated SAR radar echo based on forward method

Technical Field

The invention belongs to the technical field of radar signal processing, and particularly relates to a rapid implementation method of an SAR (synthetic aperture radar) echo imitating a three-dimensional scene based on a forward method, which is suitable for simulating a large-scene three-dimensional terrain SAR echo.

Background

Synthetic Aperture Radar (SAR) can observe a target or a scene all the day and all the weather, and has huge application space in disaster monitoring and resource exploration, especially in the military and civil fields. With the increasing development of Synthetic Aperture Radar (SAR) technology, SAR echo signals meeting specific parameters are required to be obtained when a novel SAR design scheme and an imaging algorithm are verified and evaluated, and the SAR echo signals cannot be directly obtained by radar recording echoes, so that SAR echo simulation work is required to be carried out; SAR echo simulation is divided into a reverse method and a forward method, wherein the reverse method is used for performing echo recording based on a two-dimensional SAR image, calculation of backscattering coefficients is simplified, elevation information of a scene is lost, finally obtained information amount is only a subset of a reference SAR image, and imaging effect is limited; the forward method is used for recording echoes based on three-dimensional data of an original terrain structure, can realize scene reconstruction, and has important effects in complex terrain detection and scene matching guidance. Compared with the reverse method, the forward method is complex, but has the advantage that the reverse method cannot be realized, so that the research on the forward echo simulation of the three-dimensional ground scene is very necessary.

However, when forward echo simulation is performed on a three-dimensional ground large-scene target, the first problem faced is how to effectively simulate the process of backscattering. There are many methods for calculating the electromagnetic scattering, such as physical optics method, geometric diffraction method, etc., all involve the concrete scattering mechanism, the solving process is relatively difficult; considering that the synthetic aperture radar often uses a ground scene as an observation object, the reference value of the backscattering coefficient can be estimated by adopting the measured data in combination with the factors such as the radar working frequency, the beam ray side deflection angle and the like, the calculated amount is huge, and in addition, the problems of how to judge shadow shielding, how to quickly realize echo simulation of a large scene and the like exist.

Disclosure of Invention

Aiming at the defects in the prior art, the invention aims to provide a rapid implementation method of the three-dimensional scene SAR (synthetic aperture radar) simulated radar echo based on the forward method, and the rapid implementation method of the three-dimensional scene SAR simulated radar echo based on the forward method can rapidly obtain the three-dimensional scene SAR radar echo with high signal-to-noise ratio.

In order to achieve the technical purpose, the invention is realized by adopting the following technical scheme.

A rapid implementation method of an SAR (synthetic aperture radar) echo imitating a three-dimensional scene based on a forward method comprises the following steps:

step 1, acquiring the aircraft flight speed V, SAR of an SAR (synthetic aperture radar) and the aircraft flight direction of a T, SAR radar synthetic aperture azimuth time and digital elevation model DEM (digital elevation model) data in the SAR beam range; the digital elevation model DEM data in the SAR radar beam range comprises K point targets; k is a natural number;

step 2, performing interpolation processing on the digital elevation model DEM data in the SAR radar beam range to obtain interpolated digital elevation model DEM data; the digital elevation model DEM data after interpolation processing comprises K point targets, then the K point targets are randomly distributed on a two-dimensional terrain in a geographic coordinate system, one point target is arbitrarily selected as a datum point target, the coordinate of the datum point target is marked as (x, y), z (x, y) is a digital elevation corresponding to the datum point target, an accumulated probability density function F (t) of the datum point target on a carrier of the SAR radar at the moment t is obtained through calculation, and then the roughness variance sigma of the three-dimensional terrain of the datum point target on the carrier of the SAR radar is obtained through calculation; t is a time variable, and K is a natural number;

step 3, according to the roughness variance sigma of the three-dimensional terrain of the reference point target on the carrier of the SAR, performing fine processing on the interpolated digital elevation model DEM data to obtain DEM data contained in the beam range of the SAR after fine processing;

step 4, dividing DEM data in the SAR radar beam range after the refining treatment into Q small surface elements, wherein the beam side deflection angle of the α th small surface element is theta_αThen, the beam side deflection angles are respectively calculated to be theta₁Backscatter coefficient sigma of the 1 st minor element₁The beam side angle is theta_QBackscattering coefficient sigma of the Q-th minor element_QAnd the backscattering coefficient sigma is measured₁To the backscattering coefficient sigma_QThe sum is used as the backscattering characteristic distribution of DEM data in the SAR radar beam range after the refining processing, wherein α∈ {1,2, …, Q }, and Q is a natural number;

step 5, after backscattering characteristic distribution of DEM data in the SAR beam range after refined processing is obtained, shadow region judgment is carried out on Q small surface elements by using an improved lower visual angle comparison method, and then the final actual backscattering coefficient of the Q small surface elements after shadow region judgment is obtained through calculation;

and 6, respectively calculating respective SAR (synthetic aperture radar) echoes of the Q small surface elements according to the final actual backscattering coefficients of the Q small surface elements after shadow area judgment, then accumulating the respective SAR echoes of the Q small surface elements by using a thread extrapolation technology to obtain the SAR echoes of the Q small surface elements, and completing echo simulation of the whole three-dimensional scene.

The invention has the beneficial effects that: the method utilizes the fractal Brownian motion model FBM to perform fractal interpolation on the original DEM data, ensures the simulation precision, and utilizes the small surface element model to calculate the backscattering coefficient; meanwhile, the three-dimensional terrain is considered, the shadow area is judged, and finally, a Graphic Processing Unit (GPU) architecture is introduced to carry out deep optimization, so that the simulation efficiency is greatly improved, and the real-time processing requirement of SAR radar echoes in a three-dimensional scene is met.

Drawings

The present invention will be described in further detail with reference to the accompanying drawings and specific embodiments.

FIG. 1 is a flow chart of a method for rapidly realizing an SAR (synthetic aperture radar) echo imitating a three-dimensional scene based on a forward method;

FIG. 2(a) is a schematic of raw DEM data;

FIG. 2(b) is a diagram showing the result of fractal interpolation;

FIG. 3(a) shows DEM data in SAR radar beam range at 2^gM×2^gA distribution map in an N-dimensional zero matrix;

FIG. 3(b) is a diagram illustrating interpolation results of odd-numbered coordinate points;

FIG. 3(c) is a schematic diagram of a three-dimensional result of DEM data fractal interpolation;

FIG. 4 is a schematic diagram of the spatial geometry of the small surface element divided in the present invention;

FIG. 5 is a schematic diagram of shadow occlusion using a downward view based;

FIG. 6(a) a schematic diagram of a three-dimensional simulation scene setup; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, and the unit is a sampling unit;

fig. 6(b) a diagram of the SAR imaging result corresponding to the trajectory 1; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, and the unit is a sampling unit;

fig. 6(c) a diagram of SAR imaging results corresponding to trace 2; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, and the unit is a sampling unit;

FIG. 7 is a schematic diagram of measured DEM data for a region simulated by measured data in the present invention; (ii) a The horizontal axis represents azimuth direction, the unit is a sampling unit, the vertical axis represents distance direction, the unit is a sampling unit, the height axis represents elevation direction, and the unit is a sampling unit;

fig. 8(a) a diagram of the SAR imaging result corresponding to the trajectory 1; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, and the unit is a sampling unit;

fig. 8(b) a schematic diagram of SAR imaging results corresponding to trace 2; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, and the unit is a sampling unit;

fig. 8(c) a diagram of SAR imaging results corresponding to trace 3; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, and the unit is a sampling unit;

FIG. 8(d) is a diagram showing SAR imaging results corresponding to trace 4; the horizontal axis represents the azimuth direction, the unit is a sampling unit, and the vertical axis represents the distance direction, the unit is a sampling unit.

Detailed Description

Referring to fig. 1, the invention is a flow chart of a method for rapidly realizing an SAR (synthetic aperture radar) echo imitating a three-dimensional scene based on a forward method; the method for quickly realizing the SAR (synthetic aperture radar) echo imitating the three-dimensional scene based on the forward method comprises the following steps of:

step 1, determining a graphic processing unit, and acquiring SAR (synthetic aperture radar) parameters from the graphic processing unit, wherein the SAR parameters comprise the aircraft flight speed V, SAR of an SAR and the aircraft flight direction of a synthetic aperture azimuth time T, SAR radar and Digital Elevation Model (DEM) data in the SAR beam range; digital Elevation Model (DEM) data within the SAR radar beam range contains K point targets; k is a natural number.

Specifically, according to the method for quickly realizing the SAR echo imitating the three-dimensional scene based on the forward method, SAR radar parameters are set in a CPU (central processing unit) processor, the SAR radar parameters comprise the aircraft flight speed V, SAR of the SAR radar and the aircraft flight direction of a synthetic aperture azimuth time T, SAR radar of the SAR radar and Digital Elevation Model (DEM) data contained in the SAR radar beam range, and accordingly, a physical memory space in a Graphic Processing Unit (GPU) is allocated.

The CPU processor is a central processing unit in a conventional computer, the CPU processor serves as a host end, and a Graphic Processing Unit (GPU) is a parallel data calculator; the CPU processor and the Graphic Processing Unit (GPU) form a cooperative processing working mode, the CPU processor is responsible for processing transactions with strong logicality and serial calculation, and the Graphic Processing Unit (GPU) is responsible for processing parallelized tasks and data calculation. The CPU processor and the Graphics Processing Unit (GPU) are connected via a new generation bus or interface standard (PCI-Express, PCI-E) bus. Wherein, the thread (thread) is the minimum unit of a Graphic Processing Unit (GPU) and is an execution unit which is relatively independently schedulable; the thread block in the GPU is composed of a plurality of threads, and all threads in the same thread block share one memory and are communicated or quickly synchronized with each other.

Then storing the aircraft flight speed V, SAR of the SAR radar and the aircraft flight direction of the SAR radar synthetic aperture azimuth time T, SAR radar and DEM data contained in the SAR radar beam range into a Graphic Processing Unit (GPU); referring to fig. 2(a), a diagram of raw DEM data; DEM data contained in the SAR radar beam range contains K point targets; using a shared memory in a Graphic Processing Unit (GPU) when the distance direction unit number of the SAR radar beam is less than 4096 points; and when the distance direction unit number of the SAR radar beam is larger than 4096, calling a global memory in a Graphic Processing Unit (GPU), wherein the global memory and the shared memory are high-speed memories of the GPU.

Step 2, performing interpolation processing on Digital Elevation Model (DEM) data in the SAR radar beam range to obtain interpolated Digital Elevation Model (DEM) data; the Digital Elevation Model (DEM) data after interpolation processing comprises K point targets, then the K point targets are randomly distributed on a two-dimensional terrain in a geographic coordinate system, one point target is arbitrarily selected as a datum point target, the coordinate of the datum point target is marked as (x, y), z (x, y) is a digital elevation corresponding to the datum point target, an accumulated probability density function F (t) of the datum point target on a carrier of the SAR radar at the moment t is obtained through calculation, and then roughness variance sigma of the three-dimensional terrain of the datum point target on the carrier of the SAR radar is obtained through calculation; t is a time variable, and K is a natural number.

Specifically, since Digital Elevation Model (DEM) data has statistical self-similarity (i.e., random fractal characteristics), the Digital Elevation Model (DEM) data satisfies the statistical characteristics of a Fractal Brownian Motion (FBM) model; carrying out interpolation processing on Digital Elevation Model (DEM) data contained in a SAR radar beam range by utilizing the statistical self-similarity of the Fractal Brownian Motion (FBM) model to obtain interpolated Digital Elevation Model (DEM) data; referring to fig. 2(b), a diagram of a fractal interpolation result is shown; the Digital Elevation Model (DEM) data after interpolation processing comprises K point targets and is a fine processing result of the Digital Elevation Model (DEM) data; then, the K point targets are randomly distributed on a two-dimensional terrain in a geographic coordinate system, one of the point targets is arbitrarily selected as a datum point target, the coordinate of the datum point target is recorded as (x, y), z (x, y) is a digital elevation corresponding to the datum point target, and an accumulated probability density function F (t) of the datum point target on an aircraft of the SAR radar at the moment t is obtained through calculation, wherein the expression of the function is as follows:

$F\left(t\right)=P\&lsqb;\frac{z(x+\mathrm{\&Delta;}x,y+\mathrm{\&Delta;}y)-z(x,y)}{\left|\right|(\mathrm{\&Delta;}x,\mathrm{\&Delta;}y)|{|}^H}<t\&rsqb;$

wherein, P [ ·]Representing a probability distribution function, z (x, y) representing the digital elevation corresponding to the reference point target, and (x, y) representing the coordinate of the reference point target on the two-dimensional terrain, wherein the cumulative probability density function F (t) of the reference point target on the carrier of the SAR radar at the time t obeys N (0, sigma)²) The sigma represents the roughness variance of the three-dimensional terrain, the variance sigma of a flat area is smaller, and the sigma value of the variance sigma of a steep area is larger; h represents the self-similarity parameter of the reference point target on the carrier of the SAR radar, 0<H<1, self-similarity parameters H of a reference point target on an aircraft of the SAR indicate the self-similarity degree of the terrain, and can be used for describing the roughness degree of the ground, the self-similarity parameters H of a smooth area are large, and the self-similarity parameters H of a rough area are small;representing a sample of the reference point object over the two-dimensional terrain, △ x representing the cross of the reference point object over the two-dimensional terrainCoordinate interval length, △ y represents the length of the vertical coordinate interval of the reference point target on the two-dimensional terrain, and D is the terrain fractal dimension of the reference point target on the two-dimensional terrain, and D is 3-H.

Then, calculating to obtain the fractal Brown model equivalent expression of the reference point target on the SAR radar carrier according to the cumulative probability density function F (t) of the reference point target on the SAR radar carrier at the time t by utilizing the statistical self-similarity of the Fractal Brown Motion (FBM) model:

E[|z(x+△x,y+△y)-z(x,y)|]·||(△x,△y)||^-H＝C

wherein E (·) represents a mathematical expectation, z (x, y) represents a digital elevation corresponding to the reference point target, Δ x represents an abscissa interval length of the reference point target on the two-dimensional terrain, and Δ y represents an ordinate interval length of the reference point target on the two-dimensional terrain; c represents the mean value of the random variable, H represents the self-similarity parameter of the target of the reference point on the carrier of the SAR radar, and 0< H < 1.

Taking logarithms of the reference point target on two sides of the fractal Brown model equivalent on the SAR radar carrier at the same time to obtain the logarithmic form of the fractal Brown model equivalent:

lnE [ | z (X +. DELTA.x, Y +. DELTA.y) -z (X, Y) | ] -Hln | (. DELTA.x, Y) | | lnC, and then a regression line equation with X ═ ln | (. DELTA.x, Y) | | | and Y ═ lnE [ | -z (X +. DELTA.x, Y +. DELTA.y) -z (X, Y) | ] as variables and with the self-similarity parameter H of the reference point target on the vehicle of the SAR radar as a slope is obtained, Δ X represents the abscissa interval length of the reference point target on the two-dimensional terrain, Δ Y represents the ordinate interval length of the reference point target on the two-dimensional terrain, and lnC represents the intercept on the Y axis of the regression line equation with the self-similarity parameter H of the reference point target on the vehicle of the SAR radar as a slope.

Then, the logarithm form of the equivalent expression of the fractal Brown model is subjected to

lnE[|z(x+△x,y+△y)-z(x,y)|]Linear regression model fitting is carried out on the-Hln | (△ x, △ y) | | | (lnC), and calculation is carried out to obtainThe roughness variance σ of the three-dimensional terrain of the reference point target on the vehicle of the SAR radar,

and 3, according to the roughness variance sigma of the three-dimensional terrain of the reference point target on the carrier of the SAR, performing fine processing on the interpolated Digital Elevation Model (DEM) data to obtain the DEM data contained in the beam range of the SAR after the fine processing.

Specifically, after determining a self-similarity parameter H of a reference point target on a carrier of the SAR radar and a roughness variance σ of a three-dimensional terrain of the reference point target on the carrier of the SAR radar, refining Digital Elevation Model (DEM) data after interpolation processing, wherein the refining processing is fractal interpolation processing, and the process is as follows:

3.1: create a 2^gM×2^gAn N-dimensional zero matrix, and then DEM data in the range of SAR radar beams is placed in the 2-dimensional zero matrix^gM×2^gThe row and column coordinates of the N-dimensional zero matrix are even numbers, and reach 2 containing K point targets^gM×2^gAn N-dimensional matrix.

Wherein g is an iteration index of the fractal interpolation, an initial value of g is 1, g ∈ {1,2, …, MAX }, MAX is a set maximum interpolation iteration, MAX is 10 in this embodiment, and as shown in fig. 3(a), fig. 3(a) shows that DEM data in a SAR radar beam range is 2^gM×2^gProfile in an N-dimensional zero matrix.

3.2: the 2 containing the K point targets is obtained by interpolation^gM×2^gAll row and column coordinates in the N matrix are height values at odd numbers, and then 2 is added^gM×2^gThe height values of all odd row-column coordinates in the N matrix and the 2 comprising the K point targets^gM×2^gN dimensional matrix as 2 after g iteration^gM×2^gAn N-dimensional interpolation matrix.

Specifically, the2 containing K point targets^gM×2^gAll row and column coordinates in the N matrix are height values at odd numbers, wherein 2 is to be^gM×2^gThe height values z (i, j) of odd row and column coordinates (i, j) in the N matrix are expressed as:

$z(i,j)=\frac{1}{4}\&lsqb;z(i-1,j-1)+z(i+1,j-1)+z(i-1,j+1)+z(i+1,j+1)\&rsqb;+\sqrt{1-2^{2H-2}}\left|\right|\mathrm{\&Delta;}x|{|}^H\mathrm{\&sigma;}\&CenterDot;G$

and then 2, calculating the target containing K points^gM×2^gData at positions where row and column coordinates (i, j) in the N-dimensional matrix are odd numbers are data at positions where the row and column coordinates (i, j) are odd numbers, and the data at the positions where the row and column coordinates (i, j) are odd numbers are data at positions 2 containing K point targets^gM×2^gThe row and column coordinates (i, j) in the N-dimensional matrix are height values at odd numbers, as shown in fig. 3(b), and fig. 3(b) is a diagram illustrating interpolation results of odd-number position coordinate points.

Wherein G represents a random variable subject to a standard normal distribution, and | △ x | represents a sample interval of the interpolated 2M × 2N-dimensional matrix, and (i, j) represents the 2M × N-dimensional matrix containing K point targets^gM×2^gThe coordinates of the rows and columns in the N-dimensional matrix are odd point target coordinates, and z (i-1, j-1) represents 2^gM×2^gThe height value at the row and column coordinates of (i-1, j-1) in the N-dimensional matrix, and z (i +1, j-1) represents 2^gM×2^gThe height value at the row and column coordinates of (i +1, j-1) in the N-dimensional matrix, and z (i-1, j +1) represents 2^gM×2^gThe height value at the row and column coordinates of (i-1, j +1) in the N-dimensional matrix, and z (i +1, j +1) represents 2^gM×2^gThe coordinates of the rows and columns in the N-dimensional matrix are height values at (i +1, j +1), i is a positive integer less than 2M and j is a positive integer less than 2N.

3.3: the 2 containing the K point targets is obtained by interpolation^gM×2^gOne of all row and column coordinates in the N-dimensional matrix is a height value at an odd number, and then 2 after 1 st iteration is carried out^gM×2^gN-dimensional interpolation matrix and said 2^gM×2^gOne of all row-column coordinates in the N-dimensional matrix is a height value at an odd number and is used as a final 2 after the g-th iteration^gM×2^gAn N-dimensional interpolation matrix is obtained, and then the final 2 of the g iteration is obtained_gM×2_gSample spacing for N-dimensional interpolation matrix

In particular, the 2 containing K point targets^gM×2^gOne of all row and column coordinates in the N-dimensional matrix is a height value at an odd number, wherein 2 is defined^gM×2^gOne of the row and column coordinates (i ', j') in the N matrixThe height value at odd number is recorded as z (i ', j'), and the expression is:

$\hat{z}(i^{\&prime;},j^{\&prime;})=\frac{1}{4}\&lsqb;z(i^{\&prime;},j^{\&prime;}-1)+z(i^{\&prime;}-1,j^{\&prime;})+z(i^{\&prime;}+1,j^{\&prime;})+z(i^{\&prime;},j^{\&prime;}+1)\&rsqb;+2^{-H/2}\sqrt{1-2^{2H-2}}\left|\right|\mathrm{\&Delta;}\hat{x}|{|}^H\mathrm{\&sigma;}\&times;G$

wherein G represents a random variable subject to a standard normal distribution;represents the final 2 after the g-th iteration^gM×2^gThe sample interval of the N-dimensional interpolation matrix, (i ', j') represents the 2-dimensional interpolation matrix containing K point targets^gM×2^gThe coordinates of rows and columns in the N-dimensional matrix have an odd number of target coordinates of a point, and z (i ', j' -1) represents 2^gM×2^gThe height value at the row-column coordinate (i ', j' -1) in the N-dimensional matrix, and z (i '-1, j') represents 2^gM×2^gThe height value at the row-column coordinate (i '-1, j') in the N-dimensional matrix, z (i '+1, j') represents 2^gM×2^gThe height value at the row and column coordinates (i '+1, j') in the N-dimensional matrix, z (i ', j' +1) represents 2^gM×2^gThe coordinates of the rows and columns in the N-dimensional matrix are height values at (i ', j' +1), i 'is a positive integer less than 2M, and j' is a positive integer less than 2N.

3.4 setting simulation accuracy | | | △ x_TIf final 2 after the g-th iteration^gM×2^gSample spacing for N-dimensional interpolation matrixGreater than or equal to the set simulation precision | | | △ x_TIf yes, add 1 to g and return to substep 3.2.

If final 2 after the g-th iteration^gM×2^gSample spacing for N-dimensional interpolation matrixLess than the set simulation precision | | | △ x_TIf l, then overlapThe generation is stopped and the final 2 after the g-th iteration at this time^gM×2^gK point targets and all height values contained in the N-dimensional interpolation matrix are used as DEM data contained in the refined SAR radar beam range, such as a DEM data fractal interpolation three-dimensional result schematic diagram shown in fig. 3 (c).

Step 4, dividing DEM data in the SAR radar beam range after the refining treatment into Q small surface elements, wherein the beam side deflection angle of the α th small surface element is theta_αThen, the beam side deflection angles are respectively calculated to be theta₁Backscatter coefficient sigma of the 1 st minor element₁The beam side angle is theta_QBackscattering coefficient sigma of the Q-th minor element_QAnd the backscattering coefficient sigma is measured₁To the backscattering coefficient sigma_QAnd the sum is used as the backscattering characteristic distribution of DEM data in the SAR radar beam range after the refining processing, wherein α∈ {1,2, …, Q } is natural numbers.

Specifically, the backscattering coefficient of the radio wave irradiation scene is obtained by combining experimental data and DEM data reflecting the three-dimensional structure of the real scene. In general, the three-dimensional ground can be regarded as being composed of a large number of plane facets, the electromagnetic scattering property of the ground scene is just the result of coherent superposition of backscatter fields of all the facets, and the concept of facet subdivision is the basis of an electromagnetic calculation method. The invention takes a quadrilateral surface element section as an example. Considering that the target irradiated by the SAR radar is a large-range three-dimensional ground scene, when the small surface elements are split, the size of the small surface elements is often larger than the signal wavelength but far smaller than the resolution unit of the SAR system, and the requirements are solved in the fine processing of DEM data. The spatial geometry of each small bin is determined by its center position and normal vector, and its backscatter coefficient can be approximated by the following method.

4.1, performing small surface element subdivision on DEM data in the SAR radar beam range after the refining treatment to obtain Q small surface elements, randomly distributing the Q small surface elements in a three-dimensional space, α∈ {1,2, …, Q }, wherein α represents the α th small surface element, the initial value of α is 1, and Q represents the refining treatmentThe number of small surface elements divided by DEM data in the wave beam range of the SAR radar is m × n scattering points, and figure 4 is a space geometric relationship schematic diagram of the small surface elements divided in the invention, wherein a coordinate systemIs the local coordinate system of the α th minor element,as a coordinateThe center of α th small surface element, the directions of X axis, Y axis and Z axis are the same as the geographic coordinate system, the point S is the radar loader position, and coordinates are setThe scattering point height of the α th minor plane is

4.2 calculating to obtain coordinatesThe plane equation of the plane element α is

$q_{i_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}},j_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}}=a\&times;i_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}{\mathrm{\&Delta;x}}_{\mathrm{\&alpha;}}+b\&times;j_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}{\mathrm{\&Delta;y}}_{\mathrm{\&alpha;}}+c$

Wherein a represents the slope of α th facet in the X-axis direction, b represents the slope of α th facet in the Y-axis direction, c represents the center height of α th facet, △ X_αDenotes the sampling interval of the α th bin in the X-axis direction, △ y_αIndicating the sampling interval of the α th small bin along the Y-axis direction,denotes the α th minor element in the coordinate system O_i”-the coordinates in XYZ and the coordinates in XYZ,representing coordinatesAt the center of the α th minor bin.

Then, a central optimal equation minF (a, b, c) of the small bin α is calculated by using a least square method, and the expression is as follows:

$\min F(a,b,c)=\underset{i^{\&prime;\&prime;}=0}{\overset{m-1}{\&Sigma;}}\underset{j^{\&prime;\&prime;}=0}{\overset{n-1}{\&Sigma;}}{(q_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}}-h_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}})}^2=\underset{i^{\&prime;\&prime;}=0}{\overset{m-1}{\&Sigma;}}\underset{j^{\&prime;\&prime;}=0}{\overset{n-1}{\&Sigma;}}{(a\&times;i^{\&prime;\&prime;}\mathrm{\&Delta;}x+b\&times;j^{\&prime;\&prime;}\mathrm{\&Delta;}y+c-h_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}})}^2$

wherein h is_i”,j”Denotes the true height of the scattering point of α th bin at coordinates (i ", j"), a denotes the slope of α th bin in the X-axis direction, b denotes the slope of α th bin in the Y-axis direction, c denotes the center height of α th bin, △ X denotes the sampling interval of α th bin in the X-axis direction, △ Y denotes the α th binThe sampling interval in the Y-axis direction, (i ", j") represents the α th bin in the coordinate system O_i”-coordinates in XYZ, O_i”Representing the center of the α th facet at coordinate (i ", j").

Then the unit normal vector of the plane where the α th facet is located is calculated to be

The incident field vector of the electric wave emitted from the radar phase center to the α th small element center isThe local incident angle of the radio wave on the α th small element is defined as the beam side deflection angle theta of the α th small element_α，SO_αIndicating the position of the radar carrier to the center length of the α th facet, → indicating the vector.

Let the actual projected area of the α th facet on the ground be S_αThen coordinate ofThe area of the α th minor element is Assuming that the experimental data value of the backscattering coefficient of the vegetation ground of the unit area vertically irradiated by the radio waves with the same property is sigma (0), the beam side deflection angle is calculated to be theta_αThe backscattering coefficient of the α th minor element is sigma_α，

${\mathrm{\&sigma;}}_{\mathrm{\&alpha;}}\&ap;\mathrm{\&sigma;}\left(0\right)\&times;{\hat{S}}_{\mathrm{\&alpha;}}\&times;{\cos }^2\left({\mathrm{\&theta;}}_{\mathrm{\&alpha;}}\right).$

4.3 add 1 to α and repeat substep 4.2 until a beam side-slip angle θ is obtained_QThe backscattering coefficient of the second Q-th small surface element Q is sigma_QAnd the beam side deflection angle obtained at this time is theta₁The backscattering coefficient of the 1 st minor element is sigma₁The beam side angle is theta_QThe backscattering coefficient of the second Q-th small surface element is sigma_QAs the backscatter characteristic distribution of DEM data within the SAR radar beam range after the refinement processing.

And 5, after backscattering characteristic distribution of DEM data in the SAR beam range after refined processing is obtained, judging the shadow region of the Q small surface elements by using an improved lower visual angle comparison method, and further calculating to obtain the final actual backscattering coefficient of the Q small surface elements after judgment of the shadow region.

Specifically, when the radio waves emitted by the radar transmitter are propagated along a straight line in space and blocked by a mountain peak, a tall building and other topographic structures in the propagation process, the electromagnetic waves cannot be received by the back surfaces of the topographic structures, the backscattering coefficients of the back surfaces cannot be obtained, and the shadow phenomenon appears on the SAR image. The judgment of the terrain structure shielding area is a quite complicated process, which can cause the rapid increase of the calculated amount in the echo simulation process, and particularly causes the judgment of a multiple shielding area along with the increase of the scene complexity, so the invention utilizes a lower viewing angle comparison method to judge the shadow area of the three-dimensional scene.

A shadow simulation schematic diagram of a three-dimensional terrain is shown in FIG. 5, and FIG. 5 is a schematic diagram of shadow occlusion using a downward viewing angle; and (3) a shadow simulation schematic diagram in a certain azimuth angle interval after SAR radar beam division.

5.1 initialization: and m 'belongs to {1,2, …, Q }, wherein m' represents the m 'th small element and also represents the index of the small element, the initial value of m' is 1, and Q represents the number of small elements of DEM data division in the SAR radar beam range after refinement.

5.2 randomly selecting the mth 'small surface element, and setting the ground distance of the mth' small surface element as P_m'The m' th small element is calculated to have a downward viewing angle of β_m'And further calculating the ground distance P of the small surface element with the ground distance smaller than the m' th small surface element in the Q-1 ground distances corresponding to the rest Q-1 small surface elements_m'Q 'minor panel, Q'<Then, the respective lower viewing angles of the Q ' small bins are calculated as Q-1, and Q ' lower viewing angles are obtained, and the maximum lower viewing angle of the Q ' lower viewing angles is recorded as β_maxIf β_m'≤β_maxIf yes, the mth small surface element is shielded; otherwise, the m' small surface element is not blocked; and then calculating to obtain a shielding judgment function of the m' small surface element:

wherein,is the azimuth angle of the m 'small element, and r is the ground distance of the m' small element. In step 4, calculating the backscattering coefficient of the mth small surface element to be sigma by a method for subdividing the small surface elements_m'Then the final backscattering coefficient of the mth small surface element after shadow region judgment is calculated to be

5.3; adding 1 to m', repeating the substep 5.2 until the final backscattering coefficient of the Q-th small surface element after the shadow region judgment is obtained to beAnd the final backscattering coefficient of the 1 st small surface element after the shadow region judgment is obtained at the momentThe final backscattering coefficient of the Q-th small surface element after shadow region judgment isAnd the final actual backscattering coefficient of the Q small surface elements after the shadow area judgment is obtained.

And 6, respectively calculating respective SAR (synthetic aperture radar) echoes of the Q small surface elements by using a conventional concentric circle algorithm according to the final actual backscattering coefficients of the Q small surface elements after shadow region judgment, then accumulating the respective SAR echoes of the Q small surface elements by using a thread extrapolation technology to obtain the SAR echoes of the Q small surface elements, completing echo simulation of the whole three-dimensional scene in a Graphic Processing Unit (GPU), and finally sending SAR echo signals in the GPU to a CPU (central processing unit) for outputting.

The effectiveness of the method of the invention is further verified by simulation experiments on the following measured data.

Simulation experiment of ideal three-dimensional model

In a simulation experiment for judging a shadow area, basic parameters of a simulation platform and an environmental parameter radar are shown in a table I; selecting a conical hill scene as shown in fig. 6(a) as a reference image, wherein the total number of pixels is 1024 × 1024; FIG. 6(a) a schematic diagram of a three-dimensional simulation scene setup; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, and the unit is a sampling unit; the scene parameters of the simulated conical hills are shown in table two.

Table-simulation platform and environmental parameters

Table two simulation conical hill scene parameter

By utilizing the rapid implementation method of the three-dimensional scene SAR echo simulation based on the forward method, the scene is scanned according to the track 1 and the track 2 in the simulation process. 1024 x 1024 point targets were scanned during the simulation. Then, respectively carrying out imaging processing on the generated echo data by using a conventional line-frequency scaling algorithm to obtain two imaging results shown in fig. 6(b) -6 (c), wherein the two images are both subjected to ground distance correction; wherein, fig. 6(b) shows the SAR imaging result obtained by scanning according to the track 1; FIG. 6(c) shows SAR imaging results obtained by trace 2 scanning; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, and the unit is a sampling unit;

from the simulation result, the SAR imaging result after the carrier of the SAR radar flies along the preset track 1 is shown in fig. 6(b), wherein a shadow area in the image is located right below the image, and the direction of the shadow is also right below, which corresponds to the irradiation of the track 1; the result of SAR imaging after the carrier of the SAR radar flies along the preset trajectory 2 is shown in fig. 6(c), where the shaded area is located at the lower right and towards the lower right corner, which corresponds to the irradiation of the trajectory 2.

(II) simulation experiment of actual DEM data

And carrying out simulation verification and correlation analysis by utilizing the actually measured DEM data of a certain area. The simulation still selects a simulation platform and environment parameters shown in table one, the basic information of the DEM data after interpolation selected in SAR echo simulation is shown in table three below, in order to compare the shadow effect, the basic information of the preset SAR radar carrier flight trajectory is shown in table four, H in the table represents the carrier flight height, beta represents the beam center line downward viewing angle, and theta represents the included angle between the flight trajectory and the X axis; the interpolated three-dimensional map of the actually measured DEM data is shown in FIG. 7, and FIG. 7 is a schematic diagram of the actually measured DEM data of a certain area simulated by the actually measured data in the invention; the horizontal axis represents the azimuth direction, the unit is a sampling unit, the vertical axis represents the distance direction, the unit is a sampling unit, the height axis represents the elevation direction, and the unit is a sampling unit.

Scene parameter of table three actually measured DEM data

Basic information of flight path of watch-four-carrier aircraft

In order to compare and distinguish the differences of the SAR radar echoes under different illumination angles, imaging results of the SAR echoes of the carrier of the SAR radar along the preset flight trajectories 1 to 4 are shown in fig. 8(a) to 8 (d); fig. 8(a) a diagram of the SAR imaging result corresponding to the trajectory 1; fig. 8(b) a schematic diagram of SAR imaging results corresponding to trace 2; fig. 8(c) a diagram of SAR imaging results corresponding to trace 3; the horizontal axis represents the azimuth direction, and fig. 8(d) is a diagram of the SAR imaging result corresponding to the track 4; horizontal axes in fig. 8(a) -8 (d) each represent an azimuth direction, and the unit is a sampling unit; the vertical axes each represent a distance direction, and the unit is a sampling unit.

Fig. 8(a) to 8(d) are imaging results of simulated SAR echoes of the aircraft in the flight trajectories 1 to 4, respectively, and all four images are corrected for the ground distance. Comparing fig. 8(a) with fig. 8(b) and fig. 8(c) with fig. 8(d), it can be found that although the SAR radar height and the downward viewing angle of the beam center are the same, the imaging result has detail difference due to different flying directions and different shadow effect; comparing fig. 8(a) with fig. 8(c) and fig. 8(b) with fig. 8(d), it can be found that the height of the carrier is not very large compared to the elevation value of the measured DEM data, the texture of the upper and lower images is basically the same, but the shadow effect is slightly different, and the shadow effect of the lower view angle of 60 ° is obvious, which fully illustrates the accuracy of the method of the present invention.

In conclusion, the simulation experiment verifies the correctness, the effectiveness and the reliability of the method.

It will be apparent to those skilled in the art that various changes and modifications may be made in the present invention without departing from the spirit and scope of the invention; thus, if such modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include such modifications and variations.

Claims (8)

1. A rapid implementation method of an SAR (synthetic aperture radar) echo imitating a three-dimensional scene based on a forward method is characterized by comprising the following steps:

step 1, acquiring the aircraft flight speed V, SAR of an SAR (synthetic aperture radar) and the aircraft flight direction of a T, SAR radar synthetic aperture azimuth time and digital elevation model DEM (digital elevation model) data in the SAR beam range; the digital elevation model DEM data in the SAR radar beam range comprises K point targets; k is a natural number;

step 2, performing interpolation processing on the digital elevation model DEM data in the SAR radar beam range to obtain interpolated digital elevation model DEM data; the digital elevation model DEM data after interpolation processing comprises K point targets, then the K point targets are randomly distributed on a two-dimensional terrain shape in a geographic coordinate system, one point target is arbitrarily selected as a datum point target, the coordinate of the datum point target is marked as (x, y), z (x, y) is a digital elevation corresponding to the datum point target, an accumulated probability density function F (t) of the datum point target on a carrier of the SAR radar at the moment t is obtained through calculation, and then the roughness variance sigma of the three-dimensional terrain of the datum point target on the carrier of the SAR radar is obtained through calculation; t is a time variable, and K is a natural number;

step 3, according to the roughness variance sigma of the three-dimensional terrain of the reference point target on the carrier of the SAR, performing fine processing on the interpolated digital elevation model DEM data to obtain DEM data contained in the beam range of the SAR after fine processing;

step 4, dividing DEM data in the SAR radar beam range after the refining treatment into Q small surface elements, wherein the beam side deflection angle of the α th small surface element is theta_αThen, the beam side deflection angles are respectively calculated to be theta₁Backscatter coefficient sigma of the 1 st minor element₁The beam side angle is theta_QBackscattering coefficient sigma of the Q-th minor element_QAnd the backscattering coefficient sigma is measured₁To the backscattering coefficient sigma_QThe sum is used as the backscattering characteristic distribution of DEM data in the SAR radar beam range after the refining processing, wherein α∈ {1,2, …, Q }, and Q is a natural number;

step 5, after backscattering characteristic distribution of DEM data in the SAR beam range after refined processing is obtained, shadow region judgment is carried out on Q small surface elements by using an improved lower visual angle comparison method, and then the final actual backscattering coefficient of the Q small surface elements after shadow region judgment is obtained through calculation;

and 6, respectively calculating respective SAR (synthetic aperture radar) echoes of the Q small surface elements according to the final actual backscattering coefficients of the Q small surface elements after shadow area judgment, then accumulating the respective SAR echoes of the Q small surface elements by using a thread extrapolation technology to obtain the SAR echoes of the Q small surface elements, and completing echo simulation of the whole three-dimensional scene.

2. The method as claimed in claim 1, wherein in step 2, the cumulative probability density function f (t) of the reference point target on the vehicle of the SAR radar at the time t is expressed as:

$F\left(t\right)=P\&lsqb;\frac{z(x+\mathrm{\&Delta;}x,y+\mathrm{\&Delta;}y)-z(x,y)}{\left|\right|(\mathrm{\&Delta;}x,\mathrm{\&Delta;}y)|{|}^H}<t\&rsqb;$

wherein, P [ ·]Representing a probability distribution function, z (x, y) representing a digital elevation corresponding to the fiducial target, (x, y) representing coordinates of the fiducial target on the two-dimensional terrain, △ x representing a length of an abscissa interval of the fiducial target on the two-dimensional terrain, △ y representing a length of an ordinate interval of the fiducial target on the two-dimensional terrainLength of coordinate interval; h denotes the self-similarity parameters of the reference point target on the aircraft of the SAR radar,representing a sample of the fiducial point target on the two-dimensional terrain.

3. The method as claimed in claim 1, wherein in step 2, the roughness variance σ of the three-dimensional terrain of the reference point target on the vehicle of the SAR radar is obtained by the following process:

calculating to obtain a fractal Brown model equivalent expression of the reference point target on the SAR radar carrier according to an accumulated probability density function F (t) of the reference point target on the SAR radar carrier at the time t:

E[|z(x+△x,y+△y)-z(x,y)|]·||(△x,△y)^-H＝C

wherein E (·) represents a mathematical expectation, z (x, y) represents a digital elevation corresponding to the reference point target, Δ x represents a length of an abscissa interval of the reference point target on the two-dimensional terrain, and Δ y represents a length of an ordinate interval of the reference point target on the two-dimensional terrain. C represents the mean value of the random variable, H represents the self-similarity parameter of the reference point target on the carrier of the SAR radar, and 0< H < 1;

taking logarithms of the reference point target on two sides of the fractal Brown model equivalent on the SAR radar carrier at the same time to obtain the logarithmic form of the fractal Brown model equivalent:

ln E [ | z (X +. DELTA.x, Y +. DELTA.y) -z (X, Y) | ] -hn | (. DELTA.x, Y) | | ln C, and then a regression line equation with X ═ ln | (. DELTA.x, Y) | | and Y ═ ln E [ | z (X +. DELTA.x, Y +. DELTA.y) -z (X, Y) | ] as variables and with the self-similarity parameter H of the reference point target on the vehicle of the SAR radar as the slope is obtained, Δ X represents the horizontal coordinate interval length of the reference point target on the two-dimensional terrain, Δ Y represents the vertical coordinate interval length of the reference point target on the two-dimensional terrain, lnC represents the intercept of the regression line equation with the self-similarity parameter H of the reference point target on the vehicle of the SAR radar as the slope on the Y axis;

then, the logarithm form of the equivalent expression of the fractal Brown model is subjected to

ln E[|z(x+△x,y+△y)-z(x,y)|]Linear regression model fitting is carried out on H ln (△ x, △ y) | | | ln C, the roughness variance sigma of the three-dimensional terrain of the reference point target on the carrier of the SAR radar is obtained through calculation,

4. the method for rapidly realizing the forward method-based three-dimensional scene imitation SAR radar echo is characterized in that in step 3, the digital elevation model DEM data after interpolation processing is refined to obtain DEM data included in a wave beam range of the SAR radar after the refinement processing, and the process is as follows:

3.1: create a 2^gM×2^gAn N-dimensional zero matrix, and then DEM data in the range of SAR radar beams is placed in the 2-dimensional zero matrix^gM×2^gThe row and column coordinates of the N-dimensional zero matrix are even numbers, and reach 2 containing K point targets^gM×2^gAn N-dimensional matrix.

Wherein g is an iteration index of fractal interpolation, the initial value of g is 1, g belongs to {1,2, …, MAX }, and MAX is a set maximum interpolation iteration number;

3.2: the 2 containing the K point targets is obtained by interpolation^gM×2^gAll row and column coordinates in the N matrix are height values at odd numbers, and then 2 is added^gM×2^gThe height values of all odd row-column coordinates in the N matrix and the 2 comprising the K point targets^gM×2^gN dimensional matrix as 2 after g iteration^gM×2^gAn N-dimensional interpolation matrix;

3.3: the 2 containing the K point targets is obtained by interpolation^gM×2^gOne of all row and column coordinates in the N-dimensional matrix is a height value at an odd number, and then 2 after 1 st iteration is carried out^gM×2^gN-dimensional interpolation matrix and said 2^gM×2^gAll rows and columns in the N-dimensional matrixOne of the coordinates has a height value at an odd number as the final 2 after the g-th iteration^gM×2^gAn N-dimensional interpolation matrix is obtained, and then the final 2 of the g iteration is obtained^gM×2^gSample spacing for N-dimensional interpolation matrix

3.4 setting simulation accuracy | | | △ x_TIf final 2 after the g-th iteration^gM×2^gSample spacing for N-dimensional interpolation matrixGreater than or equal to the set simulation precision | | | △ x_TIf yes, adding 1 to g, and returning to substep 3.2;

if final 2 after the g-th iteration^gM×2^gSample spacing for N-dimensional interpolation matrixLess than the set simulation precision | | | △ x_TIf yes, stopping iteration and obtaining the final 2 after the g-th iteration^gM×2^gAnd taking the K point targets and all height values contained in the N-dimensional interpolation matrix as DEM data contained in the SAR radar beam range after the refining treatment.

5. The method for rapidly realizing the SAR (synthetic aperture radar) echo imitating the three-dimensional scene based on the forward method as claimed in claim 4, wherein the 2 target containing K points is 2^gM×2^gAll row and column coordinates in the N matrix are height values at odd numbers, wherein 2 is to be^gM×2^gThe height values z (i, j) of odd row and column coordinates (i, j) in the N matrix are expressed as:

$z(i,j)=\frac{1}{4}\&lsqb;z(i-1,j-1)+z(i+1,j-1)+z(i-1,j+1)+z(i+1,j+1)\&rsqb;+\sqrt{1-2^{2H-2}}\left|\right|\mathrm{\&Delta;}x|{|}^H\mathrm{\&sigma;}\&CenterDot;G$

wherein G represents a random variable subject to a standard normal distribution, and | △ x | represents a sample interval of the interpolated 2M × 2N-dimensional matrix, and (i, j) represents the 2M × N-dimensional matrix containing K point targets^gM×2^gThe coordinates of the rows and columns in the N-dimensional matrix are odd point target coordinates, and z (i-1, j-1) represents 2^gM×2^gThe height value at the row and column coordinates of (i-1, j-1) in the N-dimensional matrix, and z (i +1, j-1) represents 2^gM×2^gThe height value at the row and column coordinates of (i +1, j-1) in the N-dimensional matrix, and z (i-1, j +1) represents 2^gM×2^gThe height value at the row and column coordinates of (i-1, j +1) in the N-dimensional matrix, and z (i +1, j +1) represents 2^gM×2^gHeight value at row and column coordinates (i +1, j +1) in N-dimensional matrixI is a positive integer less than 2M and j is a positive integer less than 2N.

6. The method for rapidly realizing the SAR (synthetic aperture radar) echo imitating the three-dimensional scene based on the forward method as claimed in claim 4, wherein the 2 target containing K points is 2^gM×2^gOne of all row and column coordinates in the N-dimensional matrix is a height value at an odd number, wherein 2 is defined^gM×2^gThe height value of one odd-numbered row-column coordinate (i ', j') in the N matrix is recorded as z (i ', j'), and the expression is as follows:

$\hat{z}(i^{\&prime;},j^{\&prime;})=\frac{1}{4}\&lsqb;z(i^{\&prime;},j^{\&prime;}-1)+z(i^{\&prime;}-1,j^{\&prime;})+z(i^{\&prime;}+1,j^{\&prime;})+z(i^{\&prime;},j^{\&prime;}+1)\&rsqb;+2^{-H/2}\sqrt{1-2^{2H-2}}\left|\right|\mathrm{\&Delta;}\hat{x}|{|}^H\mathrm{\&sigma;}\&times;G$

wherein G represents a random variable subject to a standard normal distribution;represents the final 2 after the g-th iteration^gM×2^gThe sample interval of the N-dimensional interpolation matrix, (i ', j') represents the 2-dimensional interpolation matrix containing K point targets^gM×2^gThe coordinates of rows and columns in the N-dimensional matrix have an odd number of target coordinates of a point, and z (i ', j' -1) represents 2^gM×2^gThe height value at the row-column coordinate (i ', j' -1) in the N-dimensional matrix, and z (i '-1, j') represents 2^gM×2^gThe height value at the row-column coordinate (i '-1, j') in the N-dimensional matrix, z (i '+1, j') represents 2^gM×2^gThe height value at the row and column coordinates (i '+1, j') in the N-dimensional matrix, z (i ', j' +1) represents 2^gM×2^gThe coordinates of the rows and columns in the N-dimensional matrix are height values at (i ', j' +1), i 'is a positive integer less than 2M, and j' is a positive integer less than 2N.

7. The method as claimed in claim 4, wherein in step 4, the distribution of backscattering characteristics of DEM data in the SAR radar beam range after the refining processing is obtained by the following steps:

4.1 making small surface element subdivision on DEM data in the SAR radar beam range after the refining treatment to obtain Q small surface elements, randomly distributing the Q small surface elements in a three-dimensional space, α∈ {1,2, …, Q }, wherein α represents the α th small surface element, the initial value of α is 1, Q represents the number of small surface elements divided by the DEM data in the SAR radar beam range after the refining treatment, each small surface element comprises m × n scattering points, and referring to FIG. 4, the spatial geometrical relationship of the small surface elements is shown, wherein a coordinate system is shown in the figureIs the local coordinate system of the α th minor element,as a coordinateThe center of α th small surface element, the directions of X axis, Y axis and Z axis are the same as the geographic coordinate system, the point S is the radar loader position, and coordinates are setThe scattering point height of the α th minor plane is

4.2 calculating to obtain coordinatesThe plane equation of the plane element α is

$q_{i_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}},j_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}}=a\&times;i_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}{\mathrm{\&Delta;x}}_{\mathrm{\&alpha;}}+b\&times;j_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}{\mathrm{\&Delta;y}}_{\mathrm{\&alpha;}}+c$

Wherein a represents the slope of α th facet in the X-axis direction, b represents the slope of α th facet in the Y-axis direction, c represents the center height of α th facet, △ X_αDenotes the sampling interval of the α th bin in the X-axis direction, △ y_αIndicating the sampling interval of the α th small bin along the Y-axis direction,denotes the α th minor element in the coordinate system O_i”-the coordinates in XYZ and the coordinates in XYZ,representing coordinatesAt the center of the α th minor panel;

then, a central optimal equation minF (a, b, c) of the small bin α is calculated by using a least square method, and the expression is as follows:

$\min F(a,b,c)=\underset{i^{\&prime;\&prime;}=0}{\overset{m-1}{\&Sigma;}}\underset{j^{\&prime;\&prime;}=0}{\overset{n-1}{\&Sigma;}}{(q_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}}-h_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}})}^2=\underset{i^{\&prime;\&prime;}=0}{\overset{m-1}{\&Sigma;}}\underset{j^{\&prime;\&prime;}=0}{\overset{n-1}{\&Sigma;}}{(a\&times;i^{\&prime;\&prime;}\mathrm{\&Delta;}x+b\&times;j^{\&prime;\&prime;}\mathrm{\&Delta;}y+c-h_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}})}^2$

wherein h is_i”,j”Denotes the position at coordinates (i ', j')α, a represents the slope of α th bin in the X-axis direction, b represents the slope of α th bin in the Y-axis direction, c represents the center height of α th bin, △ X represents the sampling interval of α th bin in the X-axis direction, △ Y represents the sampling interval of α th bin in the Y-axis direction, (i ", j") represents the sampling interval of α th bin in the coordinate system O_i”-coordinates in XYZ, O_i”Represents the center of the α th facet at coordinate (i ", j");

then the unit normal vector of the plane where the α th facet is located is calculated to be

The incident field vector of the electric wave emitted from the radar phase center to the α th small element center isThe local incident angle of the radio wave on the α th small element is defined as the beam side deflection angle theta of the α th small element_α，SO_αRepresents the length from the position of the radar carrier to the center of the α th facet, → represents a vector;

let the actual projected area of the α th facet on the ground be S_αThen coordinate ofThe area of the α th minor element is Assuming that the experimental data value of the backscattering coefficient of the vegetation ground of the unit area vertically irradiated by the radio waves with the same property is sigma (0), the beam side deflection angle is calculated to be theta_αThe backscattering coefficient of the α th minor element is sigma_α，

${\mathrm{\&sigma;}}_{\mathrm{\&alpha;}}\&ap;\mathrm{\&sigma;}\left(0\right)\&times;{\hat{S}}_{\mathrm{\&alpha;}}\&times;{\cos }^2\left({\mathrm{\&theta;}}_{\mathrm{\&alpha;}}\right).$

4.3 add 1 to α and repeat substep 4.2 until a beam side-slip angle θ is obtained_QThe backscattering coefficient of the second Q-th small surface element Q is sigma_QAnd the beam side deflection angle obtained at this time is theta₁The backscattering coefficient of the 1 st minor element is sigma₁The beam side angle is theta_QThe backscattering coefficient of the second Q-th small surface element is sigma_QAs the backscatter characteristic distribution of DEM data within the SAR radar beam range after the refinement processing.

8. The method for rapidly realizing the three-dimensional scene-imitated SAR radar echo based on the forward method according to claim 4, wherein in the step 5, the final actual backscattering coefficient of Q small surface elements after shadow judgment is obtained by the following process:

5.1 initialization: m 'belongs to {1,2, …, Q }, m' represents the mth 'small surface element and also represents the index of the small surface element, the initial value of m' is 1, and Q represents the number of small surface elements of DEM data division in the SAR radar beam range after refinement processing;

5.2 randomly selecting the m 'small surface element and setting the m' small surface elementThe ground distance of the small planar unit is P_m'The m' th small element is calculated to have a downward viewing angle of β_m'And further calculating the ground distance P of the small surface element with the ground distance smaller than the m' th small surface element in the Q-1 ground distances corresponding to the rest Q-1 small surface elements_mQ ' of ' Small Panel, Q '<Then, the respective lower viewing angles of the Q ' small bins are calculated as Q-1, and Q ' lower viewing angles are obtained, and the maximum lower viewing angle of the Q ' lower viewing angles is recorded as β_maxIf β_m'≤β_maxIf yes, the mth small surface element is shielded; otherwise, the m' small surface element is not blocked; and then calculating to obtain a shielding judgment function of the m' small surface element:

wherein,is the azimuth angle of the m 'small element, and r is the ground distance of the m' small element. In step 4, calculating the backscattering coefficient of the mth small surface element to be sigma by a method for subdividing the small surface elements_m'Then the final backscattering coefficient of the mth small surface element after shadow region judgment is calculated to be

5.3; adding 1 to m', repeating the substep 5.2 until the final backscattering coefficient of the Q-th small surface element after the shadow region judgment is obtained to beAnd the final backscattering coefficient of the 1 st small surface element after the shadow region judgment is obtained at the momentThe final backscattering coefficient of the Q-th small surface element after shadow region judgment isAnd the final actual backscattering coefficient of the Q small surface elements after the shadow area judgment is obtained.