CN110286374B

CN110286374B - Interference SAR image simulation method based on fractal Brownian motion

Info

Publication number: CN110286374B
Application number: CN201910225864.6A
Authority: CN
Inventors: 李涛; 唐新明; 陈乾福; 高小明
Original assignee: Ministry Of Natural Resources Land Satellite Remote Sensing Application Center
Current assignee: Ministry Of Natural Resources Land Satellite Remote Sensing Application Center
Priority date: 2019-03-25
Filing date: 2019-03-25
Publication date: 2021-05-28
Anticipated expiration: 2039-03-25
Also published as: CN110286374A

Abstract

The invention discloses an interference SAR image simulation method based on fractal Brownian motion, which is characterized by comprising the following steps of: step S1, simulating an external DEM of the SAR image based on the fractal Brownian motion (fBm); step S2, simulating a phase of the SAR image using the external DEM of the SAR image simulated in step S1; step S3, simulating the amplitude of the SAR image using the external DEM of the SAR image simulated in step S1; and step S4, simulating the SAR image based on the phase of the SAR image obtained through simulation in the step S2 and the amplitude of the SAR image obtained through simulation in the step S3. The parameters adopted in the simulation process of the invention have excellent engineering applicability, and can ensure that all error sources are not greatly related to the interference geometric state, so that researchers can obtain more accurate error simulation conclusions based on the invention.

Description

Interference SAR image simulation method based on fractal Brownian motion

Technical Field

The invention relates to the technical field of radar image processing, in particular to an interference SAR image simulation method based on fractal Brownian motion.

Background

Synthetic aperture radar interferometry (InSAR) technology has become one of the most effective global mapping tools. At present, the SRTM provides Digital Elevation Model (DEM) data in a global range of 56 degrees S-60 degrees N, the nominal precision of the data reaches 16m, the size of a grid is 30 meters, and the data plays an immeasurable role in global-scale scientific research and engineering application. In 2016, German space announced complete global DEM mapping with a 90% regional accuracy of 3.5 meters and a global 90% regional accuracy of 1.3 meters if the Antarctic region is not considered, the data has been available for fine military positioning and identification of military targets. The great potential of InSAR topographic mapping makes InSAR technology meet a new era in engineering application. However, the interference SAR satellite in China does not have business operation capability and is deficient in data. At present, only three SAR satellites available in China are high, but the three high are not specially designed for interference performance, so the interference performance is poor. The previous experiment results show that only less than 10 pairs of interference pairs can be formed in the high-resolution third image close to 300 scenes. Therefore, the method for obtaining the interference pair by using the analog simulation method and carrying out InSAR terrain mapping related technical research based on the interference pair is an important way for carrying out domestic SAR satellite index demonstration, data processing and measurement calibration research.

In the current stage of interference image simulation process, real DEM data, such as SRTM data, is mostly adopted to complete the generation of interference pair data. Although the simulation method can meet the basic requirements of interferometry to a certain extent, the highest spatial resolution of the SRTM DEM data is only 30 meters, the simulation requirement of high spatial resolution cannot be met, the higher-resolution DEM data cannot be acquired freely, and the simulation cost is extremely high.

Fractal Brownian motion (fBm) is a two-dimensional random motion that has self-similarity and long-range correlation. Similarly, all the geological information has better self-similarity within a certain distance in the geological research process. Particularly, for the terrain, although the difference of the terrain in a small scale is large, the similarity of the terrain is extremely high in a macroscopic scale, for example, the mountainous region generally extends for tens of kilometers or even hundreds of kilometers, and the eastern plain region of China extends for thousands of kilometers. This distribution of topography is very consistent with the characteristics of fBm.

fBm is mainly used for the simulation of backscattering coefficient in SAR images or the analysis of sea wave in SAR images, but is rarely used for the simulation process of interference images. fBm is a pure mathematical model that can characterize the statistical features of elevation in azimuth and distance directions during InSAR processing in a quantitative manner, and fBm is not used in the prior art during interferometric image simulation. The invention starts from a simulation level, and based on an fBm model, carries out terrain simulation, geometric parameter simulation, interference phase simulation, backscattering coefficient simulation and interference image simulation, provides an interference pair for interference processing, and completes full link simulation analysis of an image. The invention is intended to provide a simulation basis for the standing launch of domestic SAR satellites.

Disclosure of Invention

The invention provides an interference SAR image simulation method based on fractal Brownian motion, which comprises the following steps:

step S1, simulating an external DEM of the SAR image based on the fractal Brownian motion (fBm);

step S2, simulating a phase of the SAR image using the external DEM of the SAR image simulated in step S1;

step S3, simulating the amplitude of the SAR image using the external DEM of the SAR image simulated in step S1;

and step S4, simulating the SAR image based on the phase of the SAR image obtained through simulation in the step S2 and the amplitude of the SAR image obtained through simulation in the step S3.

Preferably, step S1 specifically includes the following sub-steps:

s1.1, selecting an interference image pair from the existing SAR images, designating one scene image as a main image and the other scene image as a slave image, and extracting imaging geometric parameters of the main image and the slave image;

s1.2, constructing a slant range-Doppler (R-D) model, and determining external DEM parameters corresponding to the main image based on the imaging geometric parameters of the main image and the constructed R-D model;

and S1.3, simulating the external DEM of the main image by using fractal Brownian motion (fBm).

Preferably, step S1.2 specifically comprises the following sub-steps:

and S1.2.1, constructing an inclined distance equation. As follows:

wherein R is₀Is the near-to-location slope distance, Δ R is the slope distance resolution, (X)_S,Y_S,Z_S) Is the satellite coordinate at the time of imaging, (X)_P,Y_P,Z_P) The three-dimensional coordinates of the ground points are obtained, and x is the distance direction pixel coordinates of the ground points in the main image;

step S1.2.2, construct the Doppler equation as follows:

wherein (VX)_S VY_S VZ_S) Satellite velocity at the time of imaging, (X)_S,Y_S,Z_S) As satellite coordinates of the imaging time, (X)_P,Y_P,Z_P) Is the three-dimensional coordinate of a ground point, f_dIs the Doppler value, λ is the radar wavelength, R₀Is a near-location slant range, Δ R is a slant range resolution, and x is a distance-direction pixel coordinate of a ground point in the main image;

step S1.2.3, construct an ellipsoid equation, as follows:

wherein (X)_P,Y_P,Z_P) Is the three-dimensional coordinate of the ground point, h is the elevation of the ground point, R_eAnd R_pRespectively the equator radius and the polar radius of the ellipsoid;

step S1.2.4, calculating the corner coordinates of the main image by using the slant range equation, the Doppler equation and the ellipsoid equation constructed in the step 1.2.1-1.2.3;

and S1.2.5, determining external DEM parameters of the main image according to the corner point coordinates of the main image.

Preferably, step S1.3 specifically comprises the following sub-steps:

step S1.3.1, calculating a cyclic block matrix capable of providing a two-dimensional matrix satisfying fBm characteristics;

s1.3.2, S1.3.2, calculating a random field having the same covariance as the cyclic blocking matrix;

and S1.3.3, correcting the random field acquired in S1.3.2 to obtain the final external DEM of the SAR image.

Preferably, step S1.3.2 specifically includes the following sub-steps:

step S1.3.2.1, calculating the eigenvalue of the cyclic block matrix;

step S1.3.2.2, generating a two-dimensional complex random field;

step S1.3.2.3, forming a characteristic value matrix by the characteristic values obtained in the step S1.3.2.1, multiplying the characteristic value matrix by the two-dimensional complex random field generated in the step S1.3.2.2, and calculating to obtain a random field with the same covariance as the cyclic block matrix;

step S1.3.2.4, performing Fourier transform on the S1.3.2.3 acquired random field, and extracting data of the positive frequency part of the random field after Fourier transform;

step S1.3.2.5, resolving an independent random field by using the positive frequency data obtained in step S1.3.2.4, and obtaining two independent random fields corresponding to the real part and the imaginary part of the positive frequency data respectively;

and S1.3.2.6, performing stretching transformation on the data of the two independent random fields to obtain a stretched and transformed random field.

Preferably, step S1.3.3 specifically includes the following sub-steps:

step S1.3.3.1, adding an incremental stationary property to the random field;

step S1.3.3.2, stretch transforming the data of the random field with the incremental stationary characteristic obtained in step S1.3.3.1;

step S1.3.3.3, adding the random field obtained in step S1.3.2.6 to the random field obtained in step S1.3.3.2;

step S1.3.3.4, stretch the random field obtained in step S1.3.3.3 to a specified elevation range.

Preferably, step S2 specifically includes the following sub-steps:

s2.1, simulating interference phases of the main image and the auxiliary image by using an external DEM of the SAR image;

s2.2, simulating the phase of the secondary image by using a random noise model;

step S2.3 of calculating a phase of the master image using the interference phase of the slave image obtained in step S2.1 and the phase of the slave image obtained in step S2.2; wherein:

the expression of the phase of the main image is

Wherein the content of the first and second substances,

for the phase of the slave image obtained in step S2.2,

the interference phase of each pixel in the pixel space of the main image,

is the noise phase.

Preferably, step S2.1 specifically comprises the following sub-steps:

step S2.1.1, converting the DEM obtained by simulation in the step S1 into SAR image space;

step S2.1.2, calculating the vertical baseline and the parallel baseline of each pixel point in the SAR image space by using the DEM converted in the step S2.1.1;

at step S2.1.3, the interference phase is calculated using the vertical and parallel baselines obtained at step S2.1.2.

Preferably, step S3 specifically includes the following sub-steps:

step S3.1, simulating the amplitude of the main image;

step S3.2, simulating the amplitude of the slave image;

wherein, step S3.1 specifically includes the following substeps:

step S3.1.1, calculating a local normal vector for each pixel in the primary image pixel space.

Local normal vector

Calculated using the following formula:

wherein the content of the first and second substances,

is a sight line direction vector corresponding to the pixel point (1, y),

is the sight line corresponding to the pixel point (0, y)The direction vector is a vector of the direction,

the sight line vector corresponding to the pixel point (1, y-1);

step S3.1.2, based on the local normal vector of the pixel point in the main image pixel space, determining whether the pixel point is a flat ground pixel; the method specifically comprises the following steps:

the angle between the local normal vector and the local coordinate vector can be represented as χ, which is expressed as:

wherein the content of the first and second substances,

is the unit vector of the local normal vector,

a unit vector which is a pixel point coordinate vector; preferably, when x<And judging the pixel point to be a flat pixel at 1 degree.

Step S3.1.3, calculating the nominal incident angle β of each pixel point as follows:

in the formula, R_HDistance of satellite to earth center, R₁Is the main image slant distance expressed as R₁＝R₀+ x.DELTA.R, wherein R₀Is the main image near-location slant distance, x is the x-direction coordinate of the pixel point, DeltaR is the slant distance resolution of the image, R_eIs the radius of curvature of the earth, and h is the elevation of the ground point.

Step S3.1.4, calculating the local incident angle θ of each pixel point based on the judgment of whether the pixel point is a flat pixel and the nominal incident angle of the pixel point, specifically as follows:

where β is the nominal angle of incidence obtained in step S3.1.3,

is a unit vector of the velocity of the satellite,

a unit vector that is a local normal vector;

step S3.1.5, calculating amplitude value A of each pixel point based on the judgment of whether the pixel point is flat pixel and the local incident angle of the pixel point_mThe method comprises the following steps:

wherein θ is the local incident angle of each pixel obtained in step S3.1.4, μ is the pixel size normalization parameter, and the expression is as follows:

wherein χ is an included angle between the local normal vector and the local coordinate vector, and θ is the local incident angle of each pixel point obtained in step S3.1.4.

Preferably, step S3 preferably further includes the following sub-steps:

step S3.1.6, calculating the amplitude value of the overlap area;

at step S3.1.7, shadow region amplitude values are calculated.

Preferably, step S4 includes the following sub-steps:

s4.1, simulating the main image;

for each pixel of the main image, the pixel value is expressed as a complex number, i.e.

F(x,y,m)＝a_m+b_mi

Wherein F (x, y, m) is the pixel value of each pixel of the main image, (x, y) is the coordinates of the pixel of the main image, m is the complex value of the pixel of the main image,

A_mthe amplitude of the main image simulated in step S3.1,

the phase of the simulated main image in the step 2.3;

step S4.2, simulating the slave image;

for each pixel of the slave image, its pixel value is expressed as a complex number, i.e.

F(x,y,s)＝a_s+b_si

Wherein F (x, y, s) is the pixel value of each pixel of the slave image, (x, y) is the coordinate of the slave image pixel, s is the complex value of the slave image pixel,

as is the slave image amplitude simulated in step S3.2,

the slave image phase simulated in step 2.2.

The parameters adopted in the simulation process are real TanDEM-X satellite interference parameters, the parameters have excellent engineering applicability, and all error sources can not be greatly related to the interference geometric state in the interference processing process, so that researchers can obtain a more accurate error simulation conclusion based on the method.

Drawings

Fig. 1 is a flowchart of an SAR interferometric image simulation method based on fractal brownian motion according to an embodiment of the present invention.

Fig. 2(a) and 2(b) are external digital elevation model fBm results of interferometric SAR images obtained based on fractal brownian motion simulation according to an embodiment of the present invention, where fig. 2(a) is the result obtained after real random field correction and fig. 2(b) is the result obtained after imaginary random field correction.

Fig. 3 is a schematic view of an observation geometry of InSAR in topographic mapping according to an embodiment of the present invention.

Fig. 4(a) is a rectangle obtained by connecting four corner points of a main image externally connected to the DEM of an interferometric SAR image obtained based on fractal brownian motion simulation according to an embodiment of the present invention; fig. 4(b) is a projected SAR master image coordinate space elevation map.

FIG. 5(a) is a simulated total phase according to an embodiment of the present invention; fig. 5(b) is a simulated terrain phase according to an embodiment of the present invention.

Fig. 6(a) is a simulated slave image phase according to an embodiment of the present invention, and fig. 6(b) is a simulated master image phase according to an embodiment of the present invention.

FIG. 7 is a diagram illustrating a method for determining flat pixels according to an embodiment of the invention.

Fig. 8(a) shows the master image amplitude simulation result, and fig. 8(b) shows the slave image amplitude simulation result.

Fig. 9(a) shows a result of superimposing the phase and amplitude of the master image, and fig. 9(b) shows a result of superimposing the phase and amplitude of the slave image.

Detailed Description

The invention is described in greater detail below, with reference to the above figures and using exemplary embodiments.

To solve the problems in the background art, the present invention provides an interferometric SAR image simulation method using fractional Brownian motion (fBm). The method adopts a fractal Brownian motion model to determine the DEM, and adopts interference geometric parameters to obtain the image interference pair, so that the method has strong practicability for image analog simulation.

The main purpose of the invention is to simulate DEM data by using fractional Brownian motion (fBm), and utilize interference geometric parameters to complete the simulation of interference images. The simulation result comprises two interfered SAR images.

As shown in fig. 1, first, the present invention firstly simulates the DEM, and during the simulation, the geometric parameters of the main image are used to determine the coverage area of the DEM and the grid size of the DEM, thereby determining the basic properties of the external DEM of the main image, and using fBm to simulate the external DEM. And then, accurately encoding the external DEM by using the geometric parameters of the main image again, and projecting the external DEM into an SLC image space to enable DEM data to correspond to SLC data one by one. Second, image phase is simulated using DEM. In the simulation process, it is considered that only phase noise is contained in the slave image, and therefore, the phase of the slave image is determined using white gaussian noise. And simulates the master image phase in conjunction with the slave image phase, noise, flat ground phase, and terrain phase. Third, the SLC image amplitude is simulated. During the simulation, the DEM data generated by fBm and the geometric parameters of the master image and the slave image are needed. Fourth, the simulated phase and amplitude are combined into a master-slave SLC image. And fifthly, performing conventional interference processing on the master-slave SLC image to acquire DEM data, performing cross validation on the DEM data and the simulated DEM, and giving a corresponding precision report.

The interference SAR image simulation method based on fractal Brownian motion specifically comprises the following steps.

And step S1, simulating an external Digital Elevation Model (DEM) of the interference SAR image based on the fractal Brownian motion (fBm).

Firstly, selecting a pair of interference data from the existing SAR images, and randomly appointing one scene image as a reference image, namely a main image; the other scene image is used as a slave image. And secondly, determining external DEM parameters of the main image, wherein the external DEM parameters comprise longitude and latitude, longitude and latitude resolution, row and column number, elevation range and the like of the upper left corner of the external DEM of the main image. Finally, an external DEM was simulated using fBm.

Specifically, step S1 includes the steps of:

step S1.1, selecting an interference image pair from the existing SAR images, designating one scene image as a main image and the other scene image as a slave image, and extracting imaging geometric parameters of the main image and the slave image. It should be noted that the master and slave images can be arbitrarily designated by those skilled in the art according to actual needs.

According to the preferred embodiment of the invention, an interference image pair of TanDEM-X in Rituchi county of Nauchi, Sichuan, is selected, the imaging time of the interference pair is 2013, 9 months and 12 days, an image corresponding to a TDX-1 sensor in the interference image pair is designated as a main image, and an image corresponding to a TSX-1 sensor in the interference image pair is designated as a slave image. The geometric parameters of the two images are extracted, and the imaging geometric parameters of the master image and the slave image are respectively shown in table 1 and table 2.

TABLE 1 imaging geometry of the principal image

TABLE 2 imaging geometry from images

S1.2, constructing a slant range-Doppler (R-D) model, and determining external DEM parameters corresponding to the main image based on the imaging geometric parameters of the main image and the constructed R-D model.

Specifically, the imaging geometric parameters of the main image given in step S1.1 (see table 1) are used, and coordinates of four corner points of the main image are determined by combining with an R-D model commonly used in the SAR imaging process. The R-D model is an important equation for SAR image positioning, and establishes pixel coordinates (X, y) and ground coordinates (X) of a certain point P by using the slant range and Doppler parameters provided in S1.1 and combining with ellipsoid parameters of ground points_P,Y_P,Z_P) The corresponding relation between them. The R-D model contains three submodels in total, namely a slope distance equation, a Doppler equation and an ellipsoid equation. Step S1.2 specifically includes the following steps:

and S1.2.1, constructing an inclined distance equation.

The slope equation expresses that the distance between the satellite coordinates and the ground coordinates at the imaging time is equal to the slope obtained in the header file. The skew equation is as follows:

wherein R is₀Is the near-to-location slope distance, Δ R is the slope distance resolution, (X)_S,Y_S,Z_S) Is the satellite coordinate at the time of imaging, (X)_P,Y_P,Z_P) The three-dimensional coordinates of the ground points are obtained, and x is the distance direction pixel coordinates of the ground points in the main image.

The resolution of each parameter is specifically as follows:

wherein the near point is inclined distance R₀And the slant resolution Δ R are known quantities, being the perigee slant and distance resolution of the main image in table 1, respectively.

Three-dimensional coordinates (X) of satellite_S,Y_S,Z_S) The specific calculation process is as follows:

the three-dimensional coordinates of the satellite are interpolated by adopting N state vectors of the satellite. Firstly, a polynomial relation between the satellite position and the satellite imaging time is constructed by using N state vectors, preferably, a 4 th-order polynomial is adopted, and taking a state vector of a satellite coordinate X axis of the imaging time as an example, the relation between the state vector and the satellite imaging time can be expressed as follows:

wherein (t)₁ t₂ … t_N) For the time corresponding to N state vectors, (a)₁ a₂ a₃ a₄ a₅) Is a polynomial coefficient, and the coefficient is,

for the X-axis component of the N state vectors at the time of satellite imaging, in S1.1, Table 1 gives the state vectors of the 19 main imagesAmount, so N ═ 19. Then in equation (2), the X-axis component of the state vector, the time for which 19 state vectors correspond, are known quantities, only (a)₁ a₂ a₃ a₄ a₅) Is an unknown quantity. Substituting the parameters of Table 1 into equation (2) to make

Then according to the least square principle, the method can obtain

Similarly, by substituting the state vector of the satellite coordinate Y axis at the time of imaging into equations (2) to (5), the fitting parameter in the Y direction can be obtained

And Z-direction fitting parameters

y is the azimuth pixel coordinate of the ground point in the main image, and the following relation exists between y and the imaging time t

t＝t₀+△t·(y-1) (8)

Wherein, t₀In the azimuth imaging start time, step S1.1, the azimuth start time given in table 1, Δ t is the azimuth imaging time interval, step S1.1, the azimuth one-line time given in table 1.

Thus, the slope equation is constructed as follows:

formula (7) is an example of formulas (1), (2), and (6), and the meanings of the variables in the formulas are the same as those in formulas (1), (2), and (6). To this end, the unknowns in the skew equation only contain two types, namely the ground coordinates (X)_P,Y_P,Z_P) And pixel coordinates (x, y).

Step S1.2.2, a Doppler equation is constructed.

The doppler equation describes the doppler effect due to the relative motion of the satellite and the ground. Satellite Velocity (VX) at imaging time_S VY_S VZ_S) And satellite side view direction vector ((X)_S-X_P) (Y_S-Y_P) (Z_S-Z_P) Cross product between) represents the corresponding doppler effect.

The specific expression of the constructed Doppler equation is as follows:

wherein (VX)_S VY_S VZ_S) Satellite velocity at the time of imaging, (X)_S,Y_S,Z_S) As satellite coordinates of the imaging time, (X)_P,Y_P,Z_P) Is the three-dimensional coordinate of a ground point, f_dIs the Doppler value, λ is the radar wavelength, R₀The method comprises the steps of obtaining a main image, obtaining a slope distance of a near point, obtaining a slope distance resolution, and obtaining distance direction pixel coordinates of a ground point in the main image.

In the satellite-borne SAR processing process, the Doppler effect is ensured to be 0 through on-satellite yaw guidance and ground imaging, which is called zero Doppler effect. Thus, in the present invention, the value on the right side of the equation is constantly 0. In the above equation, the unknowns contain only the satellite velocity at the time of imaging.

The satellite velocity solution at the time of imaging is as follows.

At step S1.2.2.1, a velocity vector for the satellite is calculated.

The satellite velocity vector is interpolated using N velocity vectors of the satellite. A polynomial relationship of satellite velocity and satellite imaging time is first constructed using N velocity vectors, preferably a 4 th order polynomial. The construction method is completely consistent with the process of calculating the three-dimensional coordinates of the satellite in the step S1.2.1, and details are not repeated here. Taking the main image as an example, the relationship between the satellite velocity at the imaging time and the imaging time can be obtained by the same calculation method as the satellite position at the imaging time in equations (2) to (7):

in the formula (VX)_S VY_S VZ_S) Y is the azimuth pixel coordinate of the ground point in the main image, and t is the imaging time.

Step S1.2.2.2, a Doppler equation is constructed.

In the formula (VX)_S VY_S VZ_S) Satellite velocity at the time of imaging, (X)_S,Y_S,Z_S) As satellite coordinates of the imaging time, (X)_P,Y_P,Z_P) Is the three-dimensional coordinates of the ground points.

To this end, the unknowns in the doppler equation only contain two types, namely the ground coordinates (X)_P,Y_P,Z_P) And azimuth and range pixel coordinates (x, y) of ground points in the main image.

Step S1.2.3, an ellipsoid equation is constructed.

The ellipsoid equation constructed is:

wherein the content of the first and second substances,(X_P,Y_P,Z_P) Is the three-dimensional coordinate of the ground point, h is the elevation of the ground point, R_eAnd R_pThe equatorial radius and the polar radius of the ellipsoid, respectively, in the WGS84 coordinate system the equatorial radius is 6378137.00m and the polar radius is 6356752.31 m. Thus, strictly speaking, if one wants to obtain accurate coordinates of ground points, one needs to ensure that the ground points are elevated by a known amount. In this embodiment, accurate elevation information is obtained using the simulation process in step S1.3. After accurate elevation information is obtained by simulation, the unknown parameters in the ellipsoid equation only have ground coordinates.

And S1.2.4, solving the corner coordinates of the main image by using the slant range equation, the Doppler equation and the ellipsoid equation constructed in the steps 1.2.1-1.2.3.

The pixel coordinates of the four corner points of the image correspond to four different ground coordinates. Three equations in steps S1.2.1, S1.2.2 and S1.2.3 are needed in the coordinate solving process. In the equation, the unknowns include three types, namely pixel coordinates, ground point coordinates, and ground point elevation. In this step, the pixel coordinates and the ground point elevation need to be determined to solve the ground point coordinates.

The data used for simulation is in the county of Ouchuan province, the terrain fluctuation is large, and the elevation is in the range of 3500-6000 meters. The exact elevation needs to be obtained in step S1.3 by simulation, where the summarized elevation values are used, preferably 3500 meters, with coordinates (0,0) in the upper left corner. The coordinates and the elevation are substituted into the equations (9), (12) and (13), so that the corresponding coordinate at the upper left corner is (-947867.6,5436136.9,3194792.1), and the coordinate is converted into the corresponding longitude and latitude which are (99.89089 degrees E,30.23596 degrees N).

Similarly, the elevation is set to be 3500 m, and longitude and latitude coordinates corresponding to the coordinates of the other three corner points are obtained by successive solution to be (100.05818 ° E,30.21120 ° N) (100.027203 ° E,30.05569 ° N) (99.86016 ° E,30.08055 ° N).

From step S1.2.4, the main image coverage ranges from 99.86016 ° E-100.05818 ° E,30.05569 ° N-30.23596 ° N. The external DEM range is slightly larger than the image range, and preferably, the following external DEM parameters are set:

TABLE 3 external DEM parameters

A standard fractal Brownian motion is defined as a random process B with an index of H (0< H <1) in a certain probability space, and the following conditions are met:

1. satisfying B (0,0) ═ 0 with a probability of 1, and B (x, y) is a continuous function of (x, y);

2. for any variable (x, y), (h, k) e R²With a two-dimensional increment B (x + h, y + k) -B (x, y) obeying the expectation of 0, and a variance of (h)²+k²)^HIs normally distributed with a probability satisfying

Wherein, (x, y) is independent variable of probability density distribution function, (H, k) is increment of independent variable, P is probability density function, B is Brownian motion process, H is Hurst parameter, exp is natural index, r is integral parameter,

representing the integral from minus infinity to z.

Vice versa, whenever the above conditions are met, it can be considered fBm. Methods of implementing fBm include the Hosking algorithm, Cholesky algorithm, Davies and Harte algorithm, and the like. Preferably, the present embodiment employs the Davies and Harte algorithm.

Step S1.3 specifically includes the following substeps:

at step S1.3.1, a circular blocking matrix is computed that is capable of providing a two-dimensional matrix that satisfies the 2 fBm conditions mentioned above.

The purpose of calculating the cyclic block matrix is to find a square matrix so that the square matrix provides a two-dimensional matrix that satisfies the condition of fBm. In general, the expression of the cyclic blocking matrix C is:

where γ (·) is the covariance function of the square matrix, M is the larger value between the number of rows S and the number of columns L of the DEM to be simulated, which in this embodiment is 8000 in table 3. Each cyclic block matrix C may be decomposed as follows:

C＝QΛQ^* (16)

wherein Q is an identity matrix, Q^*Is the complex conjugate of the transpose of the unit matrix Q, Λ is the diagonal matrix formed by the eigenvalues of the cyclic block matrix C, and any element in the unit matrix Q is defined as:

wherein the content of the first and second substances,

Q_j,kfor the element with coordinate (j, k) in the identity matrix Q, any element in the diagonal matrix Λ is defined as:

wherein the content of the first and second substances,

λ_kis the k-th eigenvalue, r, of the cyclic blocking matrix C_jThe number of the j +1 element of the first row of the cyclic block matrix C is exp is a natural index, M is a larger value between the row number S and the column number L of the DEM to be simulated, and in this embodiment, the row number and the column number of the DEM given in table 3 are 8000. Since the cyclic block matrix C is a positive definite symmetric matrix, the eigenvalue isPositive real numbers. If the eigenvalue of a certain matrix is

The eigenvector is the same as the cyclic block matrix C, so this matrix is still a positive definite real matrix. Such a matrix can be defined as

S＝QΛ^1/2Q^* (19)

This matrix eigenvalue will satisfy the second property of fBm, i.e., the two-dimensional delta is fit to a normal distribution, where Q is the identity matrix, Q is the complex conjugate of the transpose of the identity matrix Q, and Λ is the diagonal matrix formed by the eigenvalues of the cyclic block matrix C.

The core of this step is therefore to obtain a qualified cyclic block matrix. The matrix size is half the expected size. From S1.2.5, the size of the circular blocking matrix is 4000 × 4000. In order to construct the cyclic block matrix, row-column coordinate normalization is required to be firstly carried out, then parameters of the distance weight matrix are calculated, a single block matrix and a block cyclic matrix are constructed according to the parameters, and finally the characteristic value of the cyclic block matrix is calculated. Step S1.3.1 specifically includes the following substeps:

step S1.3.1.1, normalize the row and column coordinates of the monolithic matrix to be simulated.

The purpose of this step is to avoid the data overflow caused by too large value in the resolving process. The coordinate normalization formula is:

where (tx, ty) is the row-column coordinates of the monolithic matrix to be simulated after normalization, and R is a stretching factor, preferably, R may be 2. S is the image row number of the monolithic matrix to be simulated, L is the image row number, x is the column coordinate of the monolithic matrix to be simulated, and y is the row coordinate of the monolithic matrix to be simulated. In this embodiment, S is 4000 and L is 4000.

At step S1.3.1.2, parameters of the distance weight matrix are calculated.

The distance weight matrix parameters are input parameters of the distance weight matrix of the single block in step S1.3.2.3 to show the distance self-similarity of fBm, and the detailed distance weight matrix parameters are described below.

First, the Hurst parameter is determined fBm. The range of the Hurst parameter is 0< H <1, the smaller the H, the smaller the terrain difference, the larger the H, the larger the terrain difference, and when H is 0.5, the ordinary brownian motion is achieved. Preferably, H is 0.7.

Next, other parameters of fBm are determined,

when H is less than or equal to 0.75,

when H >0.75 of the total weight of the alloy,

wherein H is Hurst parameter, beta, alpha and c₁、c₂The intermediate parameter has no practical mathematical meaning, and is only convenient for subsequent formula description. In this embodiment, H is 0.7, corresponding β is 0, α is 0.7, and c is₁＝0.3,c₂0.7. These 4 parameters will be used in step S1.3.1.3 in succession.

At step S1.3.1.3, a monolithic matrix is computed using the parameters of step S1.3.1.2. The single block matrix is the basic construction of the cyclic block matrix and is the input parameter of step S1.3.1.4.

The purpose of computing the monolithic matrix is to provide a cyclic block matrix with self-similarity of distances. The specific process is described below.

Introducing an isotropic function r expressed as

Wherein (x, y) is the coordinate of each element in the block matrix to be simulated, (x)₀,y₀) For homogeneity in a monolithic matrix to be simulatedThe starting point coordinates of the function, i.e., (0, 0). The distance weight function of each element in the monolithic matrix to be simulated is expressed as

Wherein c is₁C2, α, and β are the distance weight parameters determined in step S1.3.1.2.

At step S1.3.1.4, a circular block matrix is computed.

And carrying out left-right mirroring and upper-lower mirroring on the distance weight function, wherein the expanded result is the cyclic block matrix. The expansion results are shown below to the right.

Where ρ is_slIs the element of the ith column and ith row of the distance weight matrix, S is the total column number of the distance weight matrix, and L is the total row number of the distance weight matrix. In this embodiment, S is 4000 and L is 4000. The final cyclic block matrix is a 8000 x 8000 matrix.

At step S1.3.2, a random field having the same covariance as the cyclic blocking matrix is calculated.

This step is to solve the identity matrix Q in equation (19) so that the identity matrix Q satisfies the definition of equation (17).

Step S1.3.2 specifically includes the following substeps:

at step S1.3.2.1, eigenvalues of the cyclic block matrix are calculated.

The eigenvalue of the cyclic blocking matrix is set to lambda, which satisfies the following formula according to the definition of the eigenvalue

Wherein E is an identity matrix, C is a cyclic block matrix, and lambda in the equation is solved to obtain a characteristic value of the matrix, and the implementation is implementedFor example, the matrix has 8000 eigenvalues, represented as (λ)₁ λ₂ …… λ₈₀₀₀) The eigenvalues provide preconditions for the creation of the matrix Q in equation (17) in step S1.3.2.

At step S1.3.2.2, a two-dimensional complex random field is generated.

Two-dimensional complex random field F_crThe expression for each element in (a) is:

F_cr(j,k)＝a+bi (27)

wherein, F_cr(j,k)Is the element of the kth row of the jth column in the two-dimensional complex random field, a is the real part of the element, b is the imaginary part of the element,

in the present embodiment, the random number matrix size is (8000 × 8000). The real and imaginary parts can thus each form a 8000 x 8000 matrix.

And S1.3.2.3, forming a characteristic value matrix by the characteristic values obtained in the step S1.3.2.1, multiplying the characteristic value matrix by the two-dimensional complex random field generated in the step S1.3.2.2, and calculating to obtain a random field with the same covariance as the cyclic block matrix.

Forming a characteristic value matrix by the characteristic values obtained in the step S1.3.2.1, and multiplying the characteristic value matrix by the two-dimensional complex random field to obtain a random field with the characteristics of a cyclic block matrix, namely the random field

Wherein, F_rFor random fields having the same covariance as the cyclic blocking matrix, F_crFor the two-dimensional complex random field generated in step S1.3.2.1, (λ)₁ λ₂ …… λ₈₀₀₀) Is the eigenvalue of the cyclic blocking matrix.

Step S1.3.2.4, for the random field F obtained S1.3.2.3_rFourier transform is carried out, and random field F after Fourier transform is extracted_rThe positive frequency part of (1).

FFT is carried out in the complex domain, and the expression is

F_(j,k)＝FFT(F_r)_(j,k),j＝1,2,...,8000,k＝1,2,...,8000 (29)

Wherein, F_(j,k)After Fourier transform, the element of the k-th row of the j-th column of the matrix, FFT is a fast Fourier transform function, F_rRandom fields with the same covariance as the circulant blocking matrix come from step S1.3.2.2. Splitting the matrix form of the above formula into element form, the transformation result of each element (j, k) can be expressed as:

wherein, F_(j,k)After Fourier transformation, the element of the k-th row of the j-th column of the matrix, F_rRandom fields with the same covariance as the circulant blocking matrix come from step S1.3.2.2. M is the matrix size of Fr, M and n are integral independent variables, exp [ alpha ], []In order to be an exponential function of the,

for the present embodiment, M is 8000, and for the above equation, the values before and after fourier transform are complex numbers, the lower left corner of the matrix after transform is a positive frequency part, and the upper right corner is a negative frequency part, and the spectral characteristics of the two parts are the same. We take only the data in the lower left corner positive frequency part for subsequent simulation analysis.

after the complex FFT, the real frequency space and the complex frequency space are orthogonal, and the random fields formed by the two are independent of each other, so that the real part and the imaginary part of the positive frequency part data in S1.3.2.4 constitute two independent random fields respectively.

And S1.3.2.6, performing stretch transformation on the data of the two independent random fields obtained in the step S1.3.2.5 to obtain a stretch transformed random field.

The step S1.3.2.5, where the acquired independent random field data is frequency domain data, has a very large value domain range, and is prone to causing data overflow in the calculation process. Therefore, data needs to be stretched to between [0,1 ]. The step S1.3.2.6 includes the following steps:

s1.3.2.6.1, subtracting the value at the origin of the independent random field to ensure that the origin of the independent random field is 0 to meet the first characteristic of fBm as described in step S1.3.

S1.3.2.6.2, data stretching is performed on the data of the two independent random fields respectively.

The stretching formula is:

wherein, F'_(j,k)Is the value of the k-th row of the j-th column after stretching, F_(e,k)Is the value of the k-th row of the jth column before stretching, F_minIs the minimum value of the modified independent random field matrix of S1.3.2.6.1, F_maxIs the maximum value of the matrix.

Step S1.3.3 specifically includes the following substeps:

at step S1.3.3.1, an incremental stationary property is added to the random field.

The purpose of random field correction is to make fBm incrementally smooth based on the noise fractal. The method for correcting the random field is to add a numerical value with increment smoothness property to each pixel value in the random field, and the numerical value can be expressed as:

where G is the delta-stationary matrix, G_(j,k)Is the element of the jth column and kth row of the delta-stationary matrix G, n₁And n₂Are respectively two random numbers, tx_jFor the sheet to be simulated in step S1.3.1.1X-coordinate, ty, of jth column after row-column coordinate normalization of block matrix_kY coordinate, c, of the k column normalized for the row-column coordinates of the monolithic matrix to be simulated in step S1.3.1.1₂The distance weight matrix parameters defined in step S1.3.1.2.

At step S1.3.3.2, the data obtained at step S1.3.3.1 for the random field with the incrementally stationary nature is stretch transformed. Specifically, the data is stretched to [0,1 ].

The stretching process is the same as step S1.3.2.6. And will not be described in detail herein.

At step S1.3.3.3, the random field obtained at step S1.3.2.5 is added to the random field obtained at step S1.3.3.2.

From step S1.2.5, the DEM elevation ranges are 3500-6000 m. Therefore, the fBm results were stretched to 3500-6000 m. And obtaining the external DEM of the main image.

Preferably, the fBm result generated by the imaginary part is selected for subsequent simulation.

The simulation results are shown in fig. 2(a) and 2 (b). Wherein fig. 2(a) is the circumscribed DEM result obtained after the real random field obtained in step S1.3.2.5 is corrected in step S1.3.3, and fig. 2(b) is the circumscribed DEM result obtained after the imaginary random field obtained in step S1.3.2.5 is corrected in step S1.3.3.

And step S2, simulating the phase of the SAR image by using the external DEM of the SAR image simulated in the step S1.

According to the theory of interference, the phase difference between the master image and the slave image constitutes the interference phase. Therefore, in order to complete the phase simulation of the master and slave images, it is necessary to complete the interference phase simulation of the master and slave images first, and the imaging geometric parameters of the master and slave images in step S1.1 and the external DEM of the simulation in step S1.3 are needed in the simulation process. The simulation from the image phase is then completed. And finally, completing the simulation of the phase of the main image.

The details of each step are as follows.

And S2.1, simulating interference phases of the main image and the auxiliary image by using the external DEM of the SAR image.

The interference geometry used in the calculation is shown in figure 3. Subsequent discussion will proceed based on this interference geometry. Wherein S is₁Antenna phase center position, S, corresponding to the main image₂Is the antenna phase center position corresponding to the slave image. B is the length of the base line, alpha is the inclination angle of the base line, theta is the side viewing angle, rho is the included angle between the base line vector and the LOS direction, B_⊥Is a vertical baseline, R and R plus delta R are respectively the distance from the phase center of the master-slave image antenna to a ground point P, H is the satellite height, R is the distance from the phase center of the master-slave image antenna to the ground point P_HDistance of satellite to earth center, R_hIs the distance from the ground point to the center of the earth, R_eIs the radius of curvature of the earth, and h is the elevation of the ground point. This figure is a basic schematic diagram of radar interference.

Step S2.1 specifically includes the following substeps:

and step S2.1.1, converting the DEM obtained by simulation in the step S1 into SAR image space.

And (3) projecting each coordinate point in the external DEM obtained in the step S2 to the SAR image pixel space by adopting the R-D models in the steps S1.2.1, S1.2.2 and S1.2.3, and keeping the pixel coordinates (x, y) of the points which satisfy that x is more than or equal to 1 and less than or equal to S and y is more than or equal to 1 and less than or equal to L, wherein S is the total column number of the DEM obtained by simulation, and L is the total row number. The elevation map after projection is shown in fig. 4, where fig. 4(a) is the external DEM obtained by simulation in step S1.3, the black frame line is the image range corresponding to the main image, and fig. 4(b) is the DEM converted into the SAR image space.

And S2.1.2, calculating the vertical baseline and the parallel baseline of each pixel point in the SAR image space by using the DEM converted in the step S2.1.1.

The interference phase mainly includes two types, i.e., a flat ground effect phase and a terrestrial phase. Wherein the flat ground effect phase is represented by a parallel baseline and the terrain phase is represented by a perpendicular baseline. This step therefore requires solving for both baseline parameters. The input value of this step is the DEM converted in step S2.1.1, and the output value is used to perform the phase solution of step S2.1.3.

Step S2.1.2 specifically includes the following substeps:

step S2.1.2.1, calculating the gaze direction vector of each pixel.

The vector of the line of sight for each pixel is expressed as follows:

wherein the content of the first and second substances,

is a sight line direction vector corresponding to the ground object at the pixel coordinate (x, y),

the phase center coordinate vector of the main image can be obtained by calculating the three-dimensional coordinate of the satellite in the step S1.2.1.

And for the coordinate vector of each pixel point, the coordinate vector is obtained by resolving through an R-D model, the resolving process is shown in steps S1.2.1-S1.2.4, and the elevation information is provided by the DEM converted in step S2.1.1.

Step S2.1.2.2, calculating the vertical line vector of each pixel point

Wherein the content of the first and second substances,

using satellite position vectors as satellite flight direction vectors

And the pixel point sight line direction vector obtained by calculation in the step S2.1.2.1

Obtained by resolving, i.e.

The solution process of the satellite position vector is seen in step s 1.2.1.

Step S2.1.2.3, calculating a baseline vector of each pixel point in the SAR image space.

Step S2.1.2.3 specifically includes the following substeps:

at step S2.1.2.3.1, a pixel vector for each point is calculated

Pixel coordinates in the slave image.

In the calculation, the slope equation in step s1.2.1 and the doppler equation in step S1.2.2 are used. I.e. two unknowns are solved using two equations. The least square or newton iteration method may be used in the calculation process, and is not described here.

Step S2.1.2.3.2, calculating a position vector corresponding to the pixel coordinates in the image

The solution process of the satellite position vector is seen in step s 1.2.1.

Step S2.1.2.3.3, calculating a baseline vector for each pixel point

Wherein the content of the first and second substances,

is the position vector of the satellite in the primary image,

is the position vector of the satellite from the image. Both satellite vectors can be obtained by resolving the three-dimensional coordinates of the satellite in the step S1.2.1.

Step S2.1.2.4, calculating the vertical baseline and the parallel baseline of each pixel point according to the baseline vector of each pixel point.

Vertical base line B of each pixel point_⊥And parallel base lines B_||The calculation can be made by:

wherein the content of the first and second substances,

for the line-of-sight vector of pixel point (x, y) calculated in step S2.1.2.1,

the vector of the pixel point perpendicular to the line of sight calculated for step S2.1.2.2.

Step S2.1.3 specifically includes the following substeps:

step S2.1.3.1, calculating the slant distance R of each pixel point of the secondary image₂。

Wherein rho is the included angle between the baseline vector and the radar visual line and is expressed as

cos^-1In the form of an inverse cosine function,

as a baseline vectorUnit vector of (i.e.

For the baseline vector calculated in step S2.1.2.3,

is the length of the modulus of the baseline vector,

being unit vectors of the pixel's sight vector, i.e.

The sight line vector of the pixel point calculated in step S2.1.2.1, (x, y) are the pixel coordinates of the pixel point,

is the modulo length of the sight line vector. R₁Is the slant distance in the main image, expressed as R₁＝R₀+x△R，R₀The main image is the slope distance of the near point, x is the x-direction coordinate of the pixel point, Δ R is the slope distance resolution of the image, and both the slope distance of the near point and the slope distance resolution are described in table 1.

Step S2.1.3.2, calculating the total phase of each pixel

Wherein λ is radar wavelength, the embodiment uses X-band data with wavelength of 3cm, where R is₁Is the slant distance in the main image, expressed as R₁＝R₀+ x.DELTA.R, wherein R₀The main image is the slope distance of the near point, x is the x-direction coordinate of the pixel point, Δ R is the slope distance resolution of the image, and both the slope distance of the near point and the slope distance resolution are described in table 1. R₂For each image from the video image calculated in step S2.1.3.1Slope distance of prime point. In this embodiment, the calculated total phase is shown in fig. 5 (a).

Step S2.1.3.3, calculating the terrain phase of each pixel

The expression of the terrain phase is as follows:

wherein the content of the first and second substances,

for the total phase of each pixel obtained in step S2.1.3.2, λ is the radar wavelength, and this embodiment uses X-band data with a wavelength of 3cm and B_||The resulting parallel baselines for each pixel point are resolved for step S2.1.2.4. The calculated topographic phase is shown in fig. 5 (b).

Step S2.2, the phase of the slave image is simulated using a random noise model.

The phase of the slave image is entirely composed of random noise. According to the relationship between interference phase noise and coherence

Wherein the content of the first and second substances,

as phase noise

The variance of (a) and gamma are coherence, and gamma is generally led to be gamma in the topographic mapping process>0.9。N_LIn order to reduce image noise from multiple views of an image, the multiple views are generally set to 1 in order to maintain phase reality during simulation. In this case, the phase noise can be solved as

Phase noise of monoscopic images according to the law of error propagation

Is expressed as

Wherein the content of the first and second substances,

is the phase noise obtained by equation 42.

Therefore, random noise with a mean value of 0 and a standard deviation of 0.242 is generated, and the magnitude thereof is 8000 × 8000, as the slave image phase. Namely, it is

Wherein

In order to obtain a phase from the image,

to obey the phase of the normal distribution, N (0,0.242) represents a standard normal distribution with a mean of 0 and a variance of 0.242, and the result is shown in fig. 6 (a).

And a step S2.3 of calculating the phase of the master image using the interference phase of the slave image obtained in the step S2.1 and the phase of the slave image obtained in the step S2.2.

The phase of the main image is expressed as

Wherein the content of the first and second substances,

for the slave image phase obtained in step S2.2,

the interference phase for each pixel in the pixel space of the primary image, obtained from step S2.1.3.2,

the generation process for the noise phase is the same as step S2.2, and is not described again. The phase of the main image obtained by the simulation is shown in fig. 6 (b).

And step S3, simulating the amplitude of the SAR image by using the external DEM of the SAR image simulated in the step S1.

The SAR image is a side view image, and there are occlusions and shadows in the image. Thus, an image is composed of three main types of features, the first type being normal pixels composed of flat and non-flat pixels, the second type being pixels of the overlap area, and the third type being pixels of the shadow area. In the calculation process, the local normal vector, the image plane normal, the nominal incident angle, the local incident angle and the like of each pixel point in the pixel space of the master image and the slave image are used. The method comprises the following specific steps.

Step S3.1, the amplitude of the main image is simulated.

Local normal vector

Calculated using the following formula:

wherein the content of the first and second substances,

is a pixel point (1, y) corresponding line-of-sight direction vectors,

is the sight line vector corresponding to the pixel point (0, y),

is a sight line direction vector corresponding to the pixel point (1, y-1),

representing a cross product of the vector. The three line-of-sight vectors in the formula may all be obtained by step S2.1.2.1. This step provides unit vectors of local normal vectors to step S3.1.2, and further determines whether a pixel belongs to a flat, non-flat, overlapping, or shaded area.

Step S3.1.2, based on the local normal vector of the pixel in the primary image pixel space, determine whether the pixel is a flat pixel.

The angle between the local normal vector and the local coordinate vector can be represented as χ, which is expressed as

Wherein cos^-1In the form of an inverse cosine function,

a unit vector of a local normal vector, expressed as

The local normal vector for each pixel point calculated for step S3.1.1,

being the modulo length of the local normal vector,

the expression of the unit vector of the pixel point coordinate vector is the same as the expression, the coordinate vector of each speed limiting point is required to be obtained by calculation, the coordinate vector is obtained by calculation through an R-D model, the calculation process is shown in steps S1.2.1-S1.2.4, and the elevation information is provided by the DEM converted in step S2.1.1. When χ is smaller than a certain angle, it can be considered as a flat pixel (as shown in fig. 7). Preferably, let χ<And a flat pixel at 1 deg..

At step S3.1.3, the nominal angle of incidence β for each pixel is calculated.

In which the parameters have the same meanings as in step S2.1, cos^-1As an inverse cosine function, R_HDistance of satellite to earth center, R₁Is the main image slant distance expressed as R₁＝R₀+ x.DELTA.R, wherein R₀The main image is the slope distance of the near point, x is the x-direction coordinate of the pixel point, Δ R is the slope distance resolution of the image, and both the slope distance of the near point and the slope distance resolution are described in table 1. R_eIs the radius of curvature of the earth, and h is the elevation of the ground point.

Step S3.1.4, calculating the local angle of incidence θ of each pixel point based on the judgment of whether the pixel point is a flat pixel and the nominal angle of incidence of the pixel point.

Where β is the nominal angle of incidence, cos, obtained in step S3.1.3^-1In the form of an inverse cosine function,

which is the unit vector of satellite velocity, obtained by step S1.2.2.1,

the unit vector, which is the local normal vector, is obtained by step S3.1.3.

Step S3.1.5, calculating amplitude value A of each pixel point based on the judgment of whether the pixel point is flat pixel and the local incident angle of the pixel point_m。

The formula expression is as follows:

wherein χ is an included angle between the local normal vector and the local coordinate vector, and θ is the local incident angle of each pixel point obtained in the step S3.1.4

Thus, the present invention has completed the calculation of the amplitudes of most pixels in the pixel space of the main image. However, the SAR image is a side view imaging, and a phenomenon such as a specific overlap or shadow is generated in the imaging process. It is therefore necessary to calculate the amplitude values of the overlap region and the shadow region in the following steps.

At step S3.1.6, the overlap region amplitude value is calculated.

At step S3.1.6.1, the overlap point is determined.

Calculating the local slope angle psi of each point, which is expressed as

Wherein cos^-1In the form of an inverse cosine function,

a unit vector, which is a local normal vector, calculated by equation (48),

a unit vector of the image plane normal, with ψ greater than 90 °, is defined as the eclipse mask point. Unit vector of image plane normal

Is expressed as

Wherein, | N_imgL is the module length of the normal of the image plane, N_imgIs the normal of the image plane and has the expression of

The calculation method for the satellite velocity vector at the imaging time of the primary image is shown in step S1.2.2.1.

And the sight line vector corresponding to the pixel point (1, y-1) is obtained.

At step S3.1.6.2, the eclipse point amplitude value is calculated.

Finding the nearest valid pixel to the left of the eclipse mask point and its amplitude value (x)_l,y_l,A_ml) The nearest valid pixel on the right and its amplitude value (x)_r,y_r,A_mr) The amplitude A of this point (x, y) is obtained using linear interpolation_mI.e. by

At step S3.1.7, the gray scale values for the shadow regions are calculated.

At step S3.1.7.1, shadow points are determined.

And defining the pixel points with the local incidence angle larger than 90 degrees as shadow points. The local angle of incidence is calculated by step S3.1.4.

Step S3.1.7.2, calculate the shadow spot amplitude value.

Find the nearest valid pixel to the left of the shadow point and its amplitude value (x)_l,y_l,A_ml) The nearest valid pixel on the right and its amplitude value (x)_r,y_r,A_mr) The amplitude A of this point (x, y) is obtained using linear interpolation_mI.e. by

And step S3.2, simulating the amplitude of the slave image.

From the image amplitude simulation process, the parameters provided in table 2 were used. The specific steps are similar to S3.1 and are not described here again. Fig. 8 shows simulation results using the master image and the slave image.

And S4.1, simulating the main image.

F(x,y,m)＝a_m+b_mi (58)

A_mthe amplitude of the main image simulated in step S3.1,

the phase of the main image simulated in step 2.3.

And S4.2, simulating the slave image.

F(x,y,s)＝a_s+b_si (59)

A_sfor the slave image amplitude simulated in step S3.2,

the slave image phase simulated in step 2.2.

Fig. 9 shows the result of phase and amplitude superposition, where fig. 9(a) is the master image result and fig. 9(b) is the slave image result.

In the prior art, SRTM data or known DEM data are mostly adopted for simulation, and error sources in the topographic surveying and mapping process cannot be discussed from a pure mathematic angle in the simulation process, so that the invention provides a simulation method based on fBm, on one hand, fBm can express the self-similarity of the data, which is consistent with the similarity of the terrain in a certain range; fBm, on the other hand, can express the randomness of the data, consistent with random fluctuations within the very small range of the terrain itself. Therefore, the mathematical characteristics of fBm can describe the terrain information well, which makes it possible to research the terrain characteristics from a theoretical level. In addition, in the traditional simulation process, parameters are generated by a simulation system instead of real parameters, the parameters adopted in the simulation process are real TanDEM-X satellite interference parameters, the parameters have excellent engineering applicability, and all error sources can be ensured not to be greatly related to the interference geometric state in the interference processing process, so that researchers can obtain a more accurate error simulation conclusion based on the method. The invention provides good reference for satellite mapping based on the interference SAR technology.

The above description is only a preferred embodiment of the present invention, and for those skilled in the art, the present invention should not be limited by the description of the present invention, which should be interpreted as a limitation.

Claims

1. An interference SAR image simulation method based on fractal Brownian motion is characterized by comprising the following steps:

step S4, simulating the SAR image based on the phase of the SAR image obtained by simulation in step S2 and the amplitude of the SAR image obtained by simulation in step S3; the step S1 specifically includes the following sub-steps:

s1.3, simulating an external DEM of the main image by using fractal Brownian motion (fBm);

wherein, step S1.2 specifically comprises the following substeps:

step S1.2.1, constructing an oblique distance equation as follows:

wherein R is₀Is the near-to-location slope distance, Δ R is the slope distance resolution, (X)_S,Y_S,Z_S) Is the time of imagingStar coordinate (X)_P,Y_P,Z_P) The three-dimensional coordinates of the ground points are obtained, and x is the distance direction pixel coordinates of the ground points in the main image;

step S1.2.2, construct the Doppler equation as follows:

wherein (VX)_S，VY_S，VZ_S) Satellite velocity at the time of imaging, (X)_S,Y_S,Z_S) As satellite coordinates of the imaging time, (X)_P,Y_P,Z_P) Is the three-dimensional coordinate of a ground point, f_dIs the Doppler value, λ is the radar wavelength, R₀Is a near-location slant range, Δ R is a slant range resolution, and x is a distance direction pixel coordinate of a ground point in the main image;

step S1.2.3, construct an ellipsoid equation, as follows:

2. The method according to claim 1, wherein step S1.3 comprises in particular the sub-steps of:

3. The method according to claim 2, wherein step S1.3.2 comprises the following sub-steps:

step S1.3.2.1, calculating the eigenvalue of the cyclic block matrix;

step S1.3.2.2, generating a two-dimensional complex random field;

step S1.3.2.4, performing Fourier transform on the S1.3.2.3 acquired random field, and extracting positive frequency data of the random field after Fourier transform;

4. The method according to claim 2, wherein step S1.3.3 comprises the following sub-steps:

step S1.3.3.1, adding an incremental stationary property to the random field;

5. The method according to claim 1, wherein step S2 comprises the following sub-steps:

the expression of the phase of the main image is

Wherein the content of the first and second substances,

for the phase of the slave image obtained in step S2.2,

the interference phase of each pixel in the pixel space of the main image,

is the noise phase.

6. The method according to claim 5, wherein step S2.1 comprises in particular the sub-steps of:

7. The method according to claim 1, wherein step S3 comprises the following sub-steps:

step S3.1, simulating the amplitude of the main image;

step S3.2, simulating the amplitude of the slave image;

wherein, step S3.1 specifically includes the following substeps:

step S3.1.1, calculating a local normal vector of each pixel point in the pixel space of the main image;

local normal vector

Calculated using the following formula:

wherein the content of the first and second substances,

is the vector of the sight line direction corresponding to the pixel point (1, y),

is the vector of the sight line direction corresponding to the pixel point (0, y),

the vector is a visual line vector corresponding to the pixel point (1, y-1);

wherein the content of the first and second substances,

is the unit vector of the local normal vector,

a unit vector which is a pixel point coordinate vector; preferably, when χ is less than 1 °, the pixel point is judged to be a flat pixel;

in the formula, R_HDistance of satellite to earth center, R₁Is the main image slant distance expressed as R₁＝R₀+ x Δ R, wherein R₀Is the main image near-location slant distance, x is the x-direction coordinate of the pixel point, Δ R is the slant distance resolution of the image, R_eThe curvature radius of the earth is shown, and h is the elevation of a ground point;

where β is the nominal angle of incidence obtained in step S3.1.3,

is a unit vector of the velocity of the satellite,

a unit vector that is a local normal vector;

step S3.1.5, based on whether the pixel point is flat ground or notJudging the pixel and the local incident angle of the pixel, and calculating the amplitude value A of each pixel_mThe method comprises the following steps:

8. The method according to claim 7, wherein step S3 preferably further comprises the sub-steps of:

step S3.1.6, calculating the amplitude value of the overlap area;

at step S3.1.7, shadow region amplitude values are calculated.

9. The method according to claim 1, wherein step S4 comprises the sub-steps of:

s4.1, simulating the main image;

F(x,y,m)＝a_m+b_mi

A_mthe amplitude of the main image simulated in step S3.1,

the phase of the simulated main image in the step 2.3;

step S4.2, simulating the slave image;

F(x,y,s)＝a_s+b_si

A_sfor the slave image amplitude simulated in step S3.2,

the slave image phase simulated in step 2.2.