Star group SAR geometric calibration method based on permanent scatterer

Abstract

The invention discloses a satellite group SAR geometric calibration method based on a permanent scatterer, and solves the problems that the processing precision of a traditional geometric calibration method in the prior art is limited by the number of calibrators, and the traditional geometric calibration method cannot be implemented in many scenes, the method comprises the following steps: constructing a satellite group SAR system, and obtaining an SAR image corresponding to each SAR satellite; according to the latitude and longitude coordinates of the vertexes of the SAR image, obtaining the row and column coordinates of the vertexes of the overlapping area in the SAR image coordinate system, and extracting a plurality of PS points in the overlapping area; calculating a PS row and column coordinate index of each PS point in each corresponding SAR image; calculating the slant distance measurement error correction amount of each satellite to correct the SAR image to obtain a corrected SAR image; observation scenes of different landforms can be adapted, and a geometric calibration result with stable and reliable performance is provided.

Description

Technical Field

[0001] The invention relates to the technical field of radar signal processing, and in particular to a constellation SAR geometric calibration method based on permanent scatterers. Background Art

[0002] The single SAR satellite geometric calibration technology estimates and corrects the slant range measurement error of a single SAR satellite based on the long-term stable ground target point in the observation scene, thereby improving the geometric positioning accuracy of a single SAR image.

[0003] The constellation SAR geometric calibration technology is a geometric calibration method proposed based on the geometric calibration of a single satellite to solve the problem of inconsistent slant range measurement accuracy between different satellites in the constellation SAR system. Its purpose is to solve the problems of plane position offset and misalignment between multiple SAR images in a large area. This technology can greatly improve the accuracy of seamless splicing of SAR images and lay the foundation for applications such as SAR image target detection and target interpretation.

[0004] Currently, the commonly used SAR geometric calibration methods are: calibration method based on active calibrator and calibration method based on passive calibrator. An active calibrator refers to a calibrator whose installation platform is equipped with a power supply device and a remote control device, and the rotation angle of the calibrator can be remotely controlled according to the beam illumination direction of the radar, so as to achieve the best calibration performance. A passive calibrator refers to a calibrator that is not equipped with a power supply device and can only adjust the installation angle manually.

[0005] For observation scenarios in flat terrain such as domestic urban areas, a sufficient number of active or passive calibrators can be manually deployed. However, it is difficult to ensure a sufficient number of calibrators in overseas areas and remote areas such as steep mountainous areas, deserts, and Gobi Deserts. The processing accuracy of traditional geometric calibration methods is limited by the number of calibrators and cannot be implemented in many scenarios. Summary of the invention

[0006] The present invention provides a constellation SAR geometric calibration method based on permanent scatterers, which solves the problem that the processing accuracy of the traditional geometric calibration method in the prior art is limited by the number of calibrators and cannot be implemented in many scenarios, and realizes the observation scene that can adapt to different terrains and provides stable and reliable geometric calibration results. At the same time, the geometric calibration technology is expanded from the traditional single satellite error correction to the system parameter correction of three or more SAR satellite constellations, which can meet new industry needs.

[0007] The present invention provides a constellation SAR geometric calibration method based on permanent scatterers, the method comprising:

[0008] Construct a constellation SAR system, according to each SAR satellite S in the constellation SAR system K,K≥3 to obtain the SAR image corresponding to each SAR satellite;

[0009] Determine the longitude and latitude coordinates corresponding to each vertex of each of the SAR images, determine the overlap area of ​​the multiple SAR images according to the longitude and latitude coordinates corresponding to each vertex of each of the SAR images, and determine the longitude and latitude range of the overlap area of ​​each vertex of the overlap area according to the overlap area;

[0010] Perform SAR reverse positioning on the longitude and latitude range of the overlapped area to obtain the SAR image row and column coordinates of each vertex of the overlapped area in the SAR image coordinate system, and establish a corresponding relationship between the longitude and latitude range of the overlapped area and the row and column coordinates (row, col) of the SAR image;

[0011] Determine the screening conditions, extract multiple PS points that meet all the screening conditions in the overlapped area, and respectively determine the PS row and column coordinate index of each PS point in the corresponding SAR images;

[0012] Constructing a constellation SAR geometric calibration equation, and inputting the PS row and column coordinate index corresponding to each PS point into the constellation SAR geometric calibration equation, and calculating the slant range measurement error correction of each satellite in the constellation SAR system;

[0013] The SAR image is corrected according to the slant range measurement error correction amount of each satellite and a correction formula to obtain a corrected SAR image.

[0014] In a possible implementation, the working mode of the constellation SAR system includes: an independent working mode of a single satellite or a collaborative working mode of multiple satellites.

[0015] In a possible implementation manner, determining the overlapping area of ​​the multiple SAR images according to the longitude and latitude coordinates corresponding to each vertex of each SAR image includes:

[0016] Determine the latitude and longitude coordinates corresponding to each vertex of any satellite in the constellation SAR system as the initial latitude and longitude range of the overlap area;

[0017] Traversing the remaining satellites in the constellation SAR system, finding the intersection of the longitude and latitude range corresponding to the satellite and the initial longitude and latitude range to obtain an intersection range;

[0018] The initial longitude and latitude range is updated to the intersection range until the traversal is completed to obtain the overlapping area of ​​multiple SAR images.

[0019] In a possible implementation, performing SAR reverse positioning on the longitude and latitude range of the overlapped area to obtain SAR image row and column coordinates of each vertex of the overlapped area in a SAR image coordinate system includes:

[0020] Determine the satellite position vector and the satellite velocity vector of the satellite corresponding to the SAR image, and determine the overlap vertex position vector of each vertex of the overlap area according to the latitude and longitude range of the overlap area;

[0021] Determine the first moment when the Doppler frequency is 0 using the satellite position vector, the satellite velocity vector and the coincident vertex position vector;

[0022] The SAR image row and column coordinates of each vertex of the overlap area in the SAR image coordinate system are calculated according to the first moment and the coordinate calculation formula.

[0023] In a possible implementation manner, the determining the first moment when the Doppler frequency is 0 by using the satellite position vector, the satellite velocity vector, and the coincident vertex position vector includes:

[0024] The first moment when the Doppler frequency is 0 is calculated using the satellite position vector, the satellite velocity vector and the coincident vertex position vector using a first calculation formula; wherein the first calculation formula is expressed as:

[0025]

[0026] Among them, v s represents the satellite velocity vector; λ represents the wavelength of the radar electromagnetic wave; p t represents the coincident vertex position vector of each vertex in the coincident area; p s represents the satellite position vector, f dc Indicates the Doppler frequency.

[0027] In a possible implementation, the coordinate calculation formula is expressed as:

[0028]

[0029]

[0030] Among them, zeroDopplerTime represents the first moment; ImageStartTime represents the radar startup time; PRT represents the radar pulse repetition period; row represents the SAR image row coordinate of the SAR image coordinate system; col represents the SAR image column coordinate of the SAR image coordinate system; p t represents the coincident vertex position vector of each vertex in the coincident area; p sRepresents the satellite position vector; NearRange represents the near-edge slant range of the observation scene; WidthSpace represents the range sampling interval.

[0031] In a possible implementation, the screening condition includes:

[0032] In SAR images obtained by different satellites observing the same scene from different angles at different times, the shape and geographical location of the target points are the same;

[0033] The PS point and the surrounding background target include a contour boundary line;

[0034] The shapes and brightness of PS points in SAR images obtained by satellites observing the same scene from different angles have the highest similarity.

[0035] In a possible implementation, constructing a constellation SAR geometric calibration equation, inputting the PS row and column coordinate index corresponding to each PS point into the constellation SAR geometric calibration equation, and calculating the slant range measurement error correction of each satellite in the constellation SAR system, includes:

[0036] Determine any one of the multiple SAR images as a reference SAR image, and loop through the multiple PS points;

[0037] Determine the first coordinate of each PS point in the rectangular coordinate system according to the row and column coordinate index of each PS point in the reference SAR image;

[0038] Determine the overlap vertex position vectors and the first coordinates of each vertex in the overlap area according to the latitude and longitude range of the overlap area, and calculate the row and column coordinate geometric observation values ​​of the mth PS point in each SAR image;

[0039] Obtaining the slant range measurement error according to the deviation between the PS row and column coordinate index of each satellite and the row and column coordinate geometric observation value and the constellation SAR geometric calibration equation;

[0040] The slant range measurement error is fitted and solved to obtain the slant range measurement error correction value of each satellite in the constellation SAR system.

[0041] In a possible implementation manner, determining the first coordinate of each PS point in the rectangular coordinate system according to the row and column coordinate index of each PS point in the reference SAR image includes:

[0042] Using the row and column coordinate index of each PS point in the reference SAR image, a first coordinate of each PS point in the rectangular coordinate system is obtained through a second calculation formula; wherein the second calculation formula is specifically expressed as:

[0043]

[0044] Among them, X S1 Indicates the X-axis coordinate of the reference satellite in the rectangular coordinate system; Y S1 Indicates the Y-axis coordinate of the reference satellite in the rectangular coordinate system; Z S1 Indicates the Z-axis coordinate of the reference satellite in the rectangular coordinate system; row m_1 Indicates the row index in the reference SAR image; col m_1 represents the column index in the reference SAR image; X m Indicates the X-axis value of the first coordinate; Y m Indicates the Y-axis value of the first coordinate; Z m represents the Z-axis value of the first coordinate; NR represents the near-edge slant distance of the observed scene; WS represents the range sampling interval; V x Indicates the X-axis velocity of the reference satellite; V y Indicates the Y-axis velocity of the reference satellite; V z represents the Z-axis speed of the reference satellite; λ represents the wavelength of the radar electromagnetic wave; R e represents the equatorial radius of the earth; h represents the surface elevation of the observation scene; R P represents the polar radius of the Earth, f dc Indicates the Doppler frequency.

[0045] In a possible implementation, the correction formula is expressed as:

[0046] R(row,col)=NR+WS·col+R cal (row,col);

[0047] Where NR represents the near-edge slant range of the observed scene; WS represents the range sampling interval; col represents the SAR image row coordinate in the SAR image coordinate system; row represents the SAR image column coordinate in the SAR image coordinate system; NR+WS·col represents the radar slant range measurement value expression before correction; R cal (row,col) represents the slant distance measurement error correction; R(row,col) represents the corrected slant distance.

[0048] One or more technical solutions provided in the present invention have at least the following technical effects or advantages:

[0049] Compared with the prior art, the present invention has the following advantages:

[0050] The present invention calculates the slant range measurement error correction amount relative to each satellite in the constellation SAR system by using the overlapped area of ​​the SAR images obtained from each satellite, corrects the SAR image according to the slant range measurement error correction amount, gets rid of the dependence on the calibrator, and can adapt to the observation data of various terrains and landforms. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 A flow chart of a constellation SAR geometric calibration method based on permanent scatterers provided in an embodiment of the present invention;

[0052] Figure 2 A 3-meter resolution SAR image of the Taihang Mountains taken by four satellites of the Nuwa constellation provided in an embodiment of the present invention;

[0053] Figure 3 A schematic diagram of extracting contour lines of overlapping areas of constellation SAR images provided by an embodiment of the present invention;

[0054] Figure 4 The permanent scatterer screening result of the constellation SAR image provided by the embodiment of the present invention;

[0055] Figure 5 A schematic diagram of the geometric calibration results of a constellation SAR based on permanent scatterers provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0056] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0057] The present invention provides a constellation SAR geometric calibration method based on permanent scatterers, such as Figure 1 As shown, the method includes the following steps S101 to S106.

[0058] S101, constructing a constellation SAR system, according to each SAR satellite S in the constellation SAR system K Obtain the SAR image corresponding to each SAR satellite;

[0059] Here, the working modes of the constellation SAR system include: an independent working mode of a single satellite or a collaborative working mode of multiple satellites.

[0060] To build a constellation SAR system, three or more SAR satellites can form a constellation. Each member satellite in the constellation can work independently or in coordination with other member satellites in the constellation. During the period when the radar is turned on, each member satellite can image and process the ground objects within the range of the radar beam. The collaborative working mode means that only one satellite in the constellation transmits a signal at the same time, and the other multiple satellites receive the echo signal of the satellite at the same time. This working mode is also called dual-base / multi-base imaging mode. Each satellite can collect a SAR image product. The constellation members are numbered from S1 to S K , K represents the total number of cluster members.

[0061] S102, determine the longitude and latitude coordinates corresponding to each vertex of each SAR image [S k lat n , S k lon n ], determine the overlap area of ​​multiple SAR images according to the latitude and longitude coordinates corresponding to each vertex of each SAR image, and determine the latitude and longitude range of the overlap area of ​​each vertex of the overlap area according to the overlap area; wherein the latitude and longitude range of the overlap area is expressed as: latitude range [LatMin, LatMax] and longitude range [LonMin, LonMax];

[0062] Specifically, in step S102, the longitude and latitude coordinates corresponding to each vertex of each SAR image are determined. k lat n , S k lon n ], determining the overlapping area of ​​multiple SAR images according to the longitude and latitude coordinates corresponding to each vertex of each SAR image, including the following steps S1021 to S1023.

[0063] S1021, determining the latitude and longitude coordinates corresponding to each vertex of any satellite in the constellation SAR system as the initial latitude and longitude range of the overlap area;

[0064] S1022, traversing the remaining satellites in the constellation SAR system, finding the intersection of the longitude and latitude range corresponding to the satellite and the initial longitude and latitude range, to obtain an intersection range;

[0065] S1023, updating the initial longitude and latitude range to the intersection range, until the traversal is completed, and obtaining the overlapping area of ​​the multiple SAR images.

[0066] For example, according to the latitude and longitude information of the four corners of the SAR image collected by each member satellite of the constellation [S k lat n , S k lon n], the latitude range [LatMin, LatMax] and longitude range [LonMin, LonMax] of the overlapping areas between the SAR images of each member satellite can be extracted.

[0067] Among them, S k Indicates the satellite number of the constellation member, and the value range of subscript k is [1, K]; lat n Represents the latitude coordinates of the four vertices of the rectangular frame of the SAR image. The value range of subscript n is [1, 4]; lon n Represents the longitude coordinates of the four vertices of the rectangular frame of the SAR image. The value range of the subscript n is [1, 4]. LatMin represents the minimum latitude of the overlapping area, LatMax represents the maximum latitude of the overlapping area, LonMin represents the minimum longitude of the overlapping area, and LonMax represents the maximum longitude of the overlapping area.

[0068] The specific implementation method is as follows:

[0069] The initial value of the overlap area latitude range [LatMin, LatMax] can be assigned to the latitude range of the S1 satellite image [S1Lat min , S1Lat max ], the initial value of the overlap area longitude range [LonMin, LonMax] can be assigned to the longitude range of the S1 satellite image [S1Lon min , S1Lon max ], then traverse the other K-1 member satellites, and update the values ​​of LatMin, LatMax, LonMin and LonMax by calculating the intersection of the longitude and latitude range sets. The update rule is:

[0070] If condition 1 is met: symbol represents the empty set, according to LatMin and S k Lat min The size relationship between them determines whether it needs to be updated. If LatMin k Lat min , let LatMin = S k Lat min , otherwise LatMin maintains the existing value and does not update; if LatMax>S k Lat max , let LatMax = S k Lat max Otherwise, LatMax maintains the existing value and does not update it.

[0071] If condition 1 is not true, then number S k ​The SAR images acquired by the satellite do not overlap with the SAR images of other satellites in the constellation. k The satellite does not participate in subsequent operations;

[0072] Secondly, judge condition 2: Is it true? If both conditions 1 and 2 are true, according to LonMin and S k Lon min The size relationship between them determines whether an update is needed. If LonMin k Lon min , let LonMin = S k Lon min Otherwise, LonMin maintains the current value and does not update; if LonMax>S k Lon max , let LonMax = S k Lon max Otherwise, LonMax maintains the current value and does not update.

[0073] If any of condition 1 and condition 2 is not true, then the judgment number S k The SAR images acquired by the satellite do not overlap with the SAR images of other satellites in the constellation. k The satellite no longer participates in subsequent operations.

[0074] S103, performing SAR reverse positioning on the latitude and longitude range of the overlapped area to obtain the SAR image row and column coordinates (row, col) of each vertex of the overlapped area in the SAR image coordinate system;

[0075] Specifically, in step S103, the longitude and latitude range of the overlapped area is subjected to SAR reverse positioning to obtain the SAR image row and column coordinates (row, col) of each vertex of the overlapped area in the SAR image coordinate system, including the following steps S1031 to S1033.

[0076] S1031, determining a satellite position vector and a satellite velocity vector of a satellite corresponding to the SAR image, and determining a coincidence vertex position vector of each vertex in the coincidence area according to the latitude and longitude range of the coincidence area;

[0077] S1032, using the satellite position vector, the satellite velocity vector and the coincidence vertex position vector to calculate the first moment when the Doppler frequency is 0 according to a first calculation formula;

[0078] Here, the first calculation formula is expressed as:

[0079]

[0080] Among them, v​s represents the satellite velocity vector, with the dimension of m / s; λ represents the wavelength of the radar electromagnetic wave, with the dimension of m; p t represents the coincident vertex position vector of each vertex in the coincident area; p s Represents the satellite position vector, with dimension m.

[0081] S1033, calculating the SAR image row and column coordinates (row, col) of each vertex of the overlapped area in the SAR image coordinate system according to the first moment and the coordinate calculation formula.

[0082] Here, the coordinate calculation formula is expressed as:

[0083]

[0084] Among them, zeroDopplerTime represents the first moment; ImageStartTime represents the radar startup time; PRT represents the radar pulse repetition period; row represents the SAR image row coordinate of the SAR image coordinate system; col represents the SAR image column coordinate of the SAR image coordinate system; p t represents the coincident vertex position vector of each vertex in the coincident area; p s Represents the satellite position vector; NearRange represents the near-edge slant range of the observation scene; WidthSpace represents the range sampling interval.

[0085] For example, in formula (1), f dc represents the Doppler frequency, the dimension is Hz; v s represents the satellite's velocity vector, with the dimension being m / s, p s The dimension is m; p t The overlapping vertex position vector of each vertex in the overlapping area has a dimension of m; λ represents the wavelength of the radar electromagnetic wave, and has a dimension of m.

[0086] According to formula (1), search for f dc The time when the Doppler time is equal to zero is recorded as zeroDopplerTime, abbreviated as ZDT. Substituting it into formula (2) we can get p t The row coordinate in the SAR image coordinate system is denoted as row. In formula (2), ImageStartTime represents the radar startup time, and its dimension is s. RPT represents the radar pulse repetition period, and its dimension is s. According to formula (3), p can be calculated. tIn the column coordinate of the SAR image coordinate system, it is denoted as col. In formula (3), NearRange represents the near-edge slant range (hereinafter abbreviated as NR) of the observed scene, and its dimension is m. WidthSpace represents the range sampling interval (hereinafter abbreviated as WS), and its dimension is m. The values ​​of the above two parameters can be found in the parameter description file attached to the SAR image.

[0087] S104, determine the screening conditions, extract multiple PS points that meet all the screening conditions in the overlapping area, and determine the PS row and column coordinate index (row m_k ,col m_k );

[0088] Specifically, in step S104, the screening conditions include:

[0089] (1) Position stability: In SAR images obtained by different satellites observing the same scene from different angles at different times, the shape and geographical location of the target point are the same;

[0090] (2) Clear contour: There is a contour boundary between the PS point and the surrounding background target, which facilitates extraction and recognition.

[0091] (3) Consistent scattering characteristics in all directions: The shapes and brightness of PS points in SAR images obtained by satellites observing the same scene from different angles have the highest similarity, which facilitates accurate matching of multiple SAR images.

[0092] S105, construct the constellation SAR geometric calibration equation, and index the PS row and column coordinates corresponding to each PS point (row m_k ,col m_k ) is input into the constellation SAR geometric calibration equation to calculate the slant range measurement error correction of each satellite in the constellation SAR system;

[0093] Specifically, in step S105, a constellation SAR geometric calibration equation is constructed based on the polynomial model, and the PS row and column coordinate index (row m_k ,col m_k ) is input into the constellation SAR geometric calibration equation to calculate the slant range measurement error correction amount of each satellite in the constellation SAR system, including the following steps S1051 to S1055.

[0094] S1051, determining any one of the multiple SAR images as a reference SAR image, and looping through multiple PS points;

[0095] S1052: Based on the row and column coordinate index of each PS point in the reference SAR image, m_1 ,colm_1 ) and the second calculation formula to obtain the first coordinate P of each PS point in the rectangular coordinate system m (X m , Y m , Z m );

[0096] S1053, determining the overlapped vertex position vectors and the first coordinate P of each vertex in the overlapped area according to the latitude and longitude range of the overlapped area m , calculate the row and column coordinate geometric observation values ​​of the mth PS point in each SAR image (row m_k_g ,col m_k_g );

[0097] S1054, according to the PS row and column coordinate index of each satellite (row m_k ,col m_k ) and the deviation between the row and column coordinate geometry observations and the constellation SAR geometry calibration equation to obtain the slant range measurement error;

[0098] S1055, fitting and solving the slant range measurement error to obtain the slant range measurement error correction value of each satellite in the constellation SAR system.

[0099] For example, PS points that meet the requirements are extracted from the overlapping areas of multiple SAR images, and the row and column coordinate indexes of the PS points in each SAR image are recorded. m_k ,col m_k );

[0100] Among them, row m_k Indicates the row coordinate index, col m_k Represents the column coordinate index, the dimension is pixel. The subscript m_k represents the pixel coordinate index of the mth PS point in the kth SAR image, m∈[1,M], k∈[1,K], M represents the total number of PS points, and K represents the total number of SAR images.

[0101] The specific implementation process is as follows:

[0102] Select any one of the K SAR images as the reference image (without loss of generality, the first SAR image is used as the reference image below), loop through the M PS points, and find the row and column coordinate index (row m_1 ,col m_1 ) and formulas (4)-(6) to calculate the coordinates (X m , Y m , Z m ), hereinafter referred to as P m .

[0103] The second calculation formula is specifically expressed as:

[0104]

[0105] Among them, X S1 Indicates the X-axis coordinate of the reference satellite in the rectangular coordinate system; Y S1 Indicates the Y-axis coordinate of the reference satellite in the rectangular coordinate system; Z S1 Indicates the Z-axis coordinate of the reference satellite in the rectangular coordinate system; row m_1 represents the row index in the reference SAR image; X m Indicates the X-axis value of the first coordinate; Y m Indicates the Y-axis value of the first coordinate; Z m represents the Z-axis value of the first coordinate; NR represents the near-edge slant distance of the observed scene; WS represents the range sampling interval; V x Indicates the X-axis velocity of the reference satellite; V y Indicates the Y-axis velocity of the reference satellite; V z represents the Z-axis speed of the reference satellite; λ represents the wavelength of the radar electromagnetic wave; R e represents the equatorial radius of the earth; h represents the surface elevation of the observation scene; R P Represents the Earth's polar radius.

[0106] In formula (4), the symbol X S1 , Y S1 , Z S1 represents the position coordinates of the reference satellite in the WGS84 space rectangular coordinate system; in formula (5), the symbol V x , V y , V z Indicates the speed information of the reference satellite. The value of the satellite position speed parameter can be found in the parameter file attached to the SAR image. dc Indicates the focused Doppler center frequency, which is usually 0 and has a dimension of Hz. e represents the equatorial radius of the earth, with a value of 6378137 and a dimension of m. h represents the surface elevation of the observation scene, which can be obtained from the reference DEM file and has a dimension of m. R p It represents the polar radius of the Earth, its value is 6357000, and its dimension is m.

[0107] Let p in formula (1) t =P m , p s and v s Substituting the position vector and velocity vector of each member satellite except the reference image in turn, the row and column coordinate geometric observation value of the mth PS point in the SAR image of each member satellite can be obtained. m_k_g,col m_k_g ), the subscript m_k_g represents the geometric observation value of the mth PS point in the kth SAR image, m∈[1,M], k∈[2,K]. Compare the column coordinate index value col of the PS point extracted from the SAR image in S104 m_k The geometric observation value of the PS point calculated by S105 based on the satellite orbit position, velocity vector and other geometric parameters col m_k_g The deviation between them can be used to estimate the slant range measurement error of each member satellite in the constellation relative to the reference satellite. The specific expression is as follows:

[0108] R cal_m_k =(col m_k -col m_k_g )·WS (7)

[0109] Among them, R cal_m_k represents the estimated S based on a single PS point with sequence number m k The satellite's slant range measurement error.

[0110] The slant range measurement error is usually caused by the refraction of electromagnetic wave signals when they pass through the atmosphere, resulting in the bending of the propagation path. The atmosphere is a non-uniform medium, and the slant range measurement error caused by the atmospheric refraction effect is a gradual error. The slant range measurement error estimated only based on a single PS point in formula (7) is difficult to reflect the variation trend of the slant range measurement error of the entire SAR image. The present invention accurately describes the radar slant range error correction value R according to the constellation SAR geometric calibration equation. cal The expression of the change trend of the row and column coordinates of the SAR image is as follows:

[0111] R cal (row,col)=c0+c1·col+c2·row+c3·col 2 +c4·row 2 +c5 row col (8)

[0112] Where row and col represent the row coordinate index and column coordinate index of the SAR image respectively, c0~c5 represent the polynomial fitting coefficients, and the process of estimating the fitting coefficients of the slant range measurement error based on M PS points is as shown in formulas (9)-(11):

[0113]

[0114] Formula (9) can be rewritten into matrix form as follows:

[0115] A·c=R (10)

[0116] Among them, the matrix A contains the row and column coordinate information of M PS points in the kth SAR image, the vector c is composed of 6 fitting coefficients c0 to c5, and the vector R represents the observed value of the slant range measurement error of M PS points. The least squares method can be used to solve the vector c, as shown in formula (11):

[0117] c=(A T A) -1 A T R (11)

[0118] In formula (11), the superscript T represents a matrix transposition operation, and the superscript -1 represents a matrix inversion operation.

[0119] S106, correcting the SAR image according to the slant range measurement error correction amount of each satellite and the correction formula to obtain a corrected SAR image.

[0120] The correction formula is expressed as:

[0121] R(row,col)=NR+WS·col+R cal (row,col) (12)

[0122] Where NR represents the near-edge slant range of the observed scene; WS represents the range sampling interval; col represents the SAR image row coordinate in the SAR image coordinate system; row represents the SAR image column coordinate in the SAR image coordinate system; NR+WS·col represents the radar slant range measurement value expression before correction; R cal (row,col) represents the slant distance measurement error correction; R(row,col) represents the corrected slant distance.

[0123] By correcting the radar slant range of each SAR image, the consistency of the geometric positioning accuracy of the constellation SAR images can be improved.

[0124] Compared with the prior art, the present invention has the following advantages:

[0125] (1) Both the active calibration method and the passive calibration method need to satisfy the condition that the observation scene contains a sufficient number of calibrators. However, it is difficult to ensure a sufficient number of calibrators in remote areas such as overseas areas and steep mountainous areas, deserts, and Gobi Deserts. The present invention gets rid of the dependence on calibrators and constructs a geometric calibration method based on permanent scatterer targets in SAR images, which can adapt to observation data of various terrains and landforms.

[0126] (2) The present invention expands the geometric calibration technology from the traditional single-satellite error correction to the system parameter correction of three or more SAR satellite constellations, which can meet new industry needs.

[0127] (3) The present invention supports multiple working modes, including a single satellite independent working mode and a multi-satellite collaborative working mode.

[0128] In a simulation experiment provided by the present invention, Figure 2 (A 3-meter resolution SAR image of the Taihang Mountains taken by four satellites of the Nuwa constellation) shows how to collect multiple SAR images of the same observation scene through the constellation SAR system. The horizontal axis in the figure represents the range direction, and the vertical axis represents the azimuth direction.

[0129] Figure 3 (Schematic diagram of the overlap area extraction of the constellation SAR image) shows the longitude and latitude coordinates of the four corners of the S1 satellite SAR image, including the maximum latitude S1Lat max , maximum longitude S1Lon max , minimum latitude S1Lat min and minimum longitude S1Lon min , and the longitude and latitude coordinates of the four corners of the overlap area are obtained by calculating the intersection of the longitude and latitude sets of multiple SAR images, including the maximum latitude LatMax of the overlap area, the maximum longitude LonMax of the overlap area, the minimum latitude LatMin of the overlap area, and the minimum longitude LonMin of the overlap area. The red rectangle in the figure indicates that the edge contour line of the overlap area is obtained by connecting the longitude and latitude coordinates of the four corners of the overlap area.

[0130] Figure 4 (Results of permanent scatterer screening in swarm SAR images) illustrates the screening of permanent scatterer targets that meet the requirements from the overlapping area of ​​multiple SAR images. The permanent scatterer targets in the figure are marked with a red "+". It can be seen that even in mountainous areas with steep terrain, the present invention can screen out a sufficient number of permanent scatterer targets, laying a foundation for improving the geometric calibration accuracy of swarm SAR.

[0131] Figure 5 (Geometric calibration results of constellation SAR based on permanent scatterers) illustrates the slant range correction of the constellation SAR system obtained by estimating the permanent scatterer targets in the overlapped area of ​​SAR images. In the figure, the horizontal axis represents the pixel coordinates in the range direction, the vertical axis represents the pixel coordinates in the azimuth direction, and the chromaticity value represents the slant range error correction of each pixel in the SAR image.

[0132] According to calculations, the constellation SAR geometric calibration method proposed in the present invention can improve the SAR image positioning accuracy from 8 to 10 meters before calibration to 3 to 4 meters. Moreover, the entire algorithm program only takes 2 minutes to process a 10,000×10,000 pixel satellite-borne SAR image, which is quite efficient.

[0133] The various embodiments in this specification are described in a progressive manner, and the same or similar parts between the various embodiments can be referred to each other, and each embodiment focuses on the differences from other embodiments. All or part of the present invention can be used in many general or special computer system environments or configurations. For example: personal computers, server computers, handheld devices or portable devices, tablet devices, mobile communication terminals, multi-processor systems, microprocessor-based systems, programmable electronic devices, network PCs, minicomputers, mainframe computers, distributed computing environments including any of the above systems or devices, etc.

[0134] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, a person skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some or all of the technical features thereof may be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the scope of the technical solutions of the present invention.

Claims

1. A constellation SAR geometric calibration method based on permanent scatterers, characterized in that: include: Construct a constellation SAR system, according to each SAR satellite S in the constellation SAR system K ,K≥3 to obtain the SAR image corresponding to each SAR satellite; Determine the longitude and latitude coordinates corresponding to each vertex of each of the SAR images, determine the overlap area of ​​the multiple SAR images according to the longitude and latitude coordinates corresponding to each vertex of each of the SAR images, and determine the longitude and latitude range of the overlap area of ​​each vertex of the overlap area according to the overlap area; Perform SAR reverse positioning on the latitude and longitude range of the overlapped area to obtain the SAR image row and column coordinates of each vertex of the overlapped area in the SAR image coordinate system; Determine the screening conditions, extract multiple PS points that meet all the screening conditions in the overlapped area, and respectively determine the PS row and column coordinate index of each PS point in the corresponding SAR images; Constructing a constellation SAR geometric calibration equation, and inputting the PS row and column coordinate index corresponding to each PS point into the constellation SAR geometric calibration equation, and calculating the slant range measurement error correction of each satellite in the constellation SAR system; The SAR image is corrected according to the slant range measurement error correction amount of each satellite and a correction formula to obtain a corrected SAR image.

2. The constellation SAR geometric calibration method based on permanent scatterers according to claim 1 is characterized in that: The working modes of the constellation SAR system include: an independent working mode of a single satellite or a collaborative working mode of multiple satellites.

3. The constellation SAR geometric calibration method based on permanent scatterers according to claim 1 is characterized in that: Determining the overlapping area of ​​the plurality of SAR images according to the latitude and longitude coordinates corresponding to each vertex of each SAR image includes: Determine the latitude and longitude coordinates corresponding to each vertex of any satellite in the constellation SAR system as the initial latitude and longitude range of the overlap area; Traversing the remaining satellites in the constellation SAR system, finding the intersection of the longitude and latitude range corresponding to the satellite and the initial longitude and latitude range to obtain an intersection range; The initial longitude and latitude range is updated to the intersection range until the traversal is completed to obtain the overlapping area of ​​multiple SAR images.

4. The constellation SAR geometric calibration method based on permanent scatterers according to claim 1 is characterized in that: The performing SAR reverse positioning on the latitude and longitude range of the overlapped area to obtain SAR image row and column coordinates of each vertex of the overlapped area in the SAR image coordinate system includes: Determine the satellite position vector and the satellite velocity vector of the satellite corresponding to the SAR image, and determine the overlap vertex position vector of each vertex of the overlap area according to the latitude and longitude range of the overlap area; Determine the first moment when the Doppler frequency is 0 using the satellite position vector, the satellite velocity vector and the coincident vertex position vector; The SAR image row and column coordinates of each vertex of the overlap area in the SAR image coordinate system are calculated according to the first moment and the coordinate calculation formula.

5. The method for geometric calibration of constellation SAR based on permanent scatterers according to claim 4, characterized in that: The method of determining the first moment when the Doppler frequency is 0 by using the satellite position vector, the satellite velocity vector and the coincident vertex position vector comprises: The first moment when the Doppler frequency is 0 is calculated using the satellite position vector, the satellite velocity vector and the coincident vertex position vector using a first calculation formula; wherein the first calculation formula is expressed as: Among them, v s represents the satellite velocity vector; λ represents the wavelength of the radar electromagnetic wave; p t represents the coincident vertex position vector of each vertex in the coincident area; p s represents the satellite position vector, f dc Indicates the Doppler frequency.

6. The constellation SAR geometric calibration method based on permanent scatterers according to claim 4 is characterized in that: The coordinate calculation formula is expressed as: Among them, zeroDopplerTime represents the first moment; ImageStartTime represents the radar startup time; PRT represents the radar pulse repetition period; row represents the SAR image row coordinate of the SAR image coordinate system; col represents the SAR image column coordinate of the SAR image coordinate system; p t represents the coincident vertex position vector of each vertex in the coincident area; p s Represents the satellite position vector; NearRange represents the near-edge slant range of the observation scene; WidthSpace represents the range sampling interval.

7. The constellation SAR geometric calibration method based on permanent scatterers according to claim 1 is characterized in that: The screening conditions include: In SAR images obtained by different satellites observing the same scene from different angles at different times, the shape and geographical location of the target points are the same; The PS point and the surrounding background target include a contour boundary line; The shapes and brightness of PS points in SAR images obtained by satellites observing the same scene from different angles have the highest similarity.

8. The constellation SAR geometric calibration method based on permanent scatterers according to claim 1 is characterized in that: The constellation SAR geometric calibration equation is constructed, and the PS row and column coordinate index corresponding to each PS point is input into the constellation SAR geometric calibration equation to calculate the slant range measurement error correction of each satellite in the constellation SAR system, including: Determine any one of the multiple SAR images as a reference SAR image, and loop through the multiple PS points; Determine the first coordinate of each PS point in the rectangular coordinate system according to the row and column coordinate index of each PS point in the reference SAR image; Determine the overlap vertex position vectors and the first coordinates of each vertex in the overlap area according to the latitude and longitude range of the overlap area, and calculate the row and column coordinate geometric observation values ​​of the mth PS point in each SAR image; Obtaining the slant range measurement error according to the deviation between the PS row and column coordinate index of each satellite and the row and column coordinate geometric observation value and the constellation SAR geometric calibration equation; The slant range measurement error is fitted and solved to obtain the slant range measurement error correction value of each satellite in the constellation SAR system.

9. The method for geometric calibration of constellation SAR based on permanent scatterers according to claim 8, characterized in that: Determining the first coordinate of each PS point in the rectangular coordinate system according to the row and column coordinate index of each PS point in the reference SAR image includes: Using the row and column coordinate index of each PS point in the reference SAR image, a first coordinate of each PS point in the rectangular coordinate system is obtained through a second calculation formula; wherein the second calculation formula is specifically expressed as: Among them, X S1 Indicates the X-axis coordinate of the reference satellite in the rectangular coordinate system; Y S1 Indicates the Y-axis coordinate of the reference satellite in the rectangular coordinate system; Z S1 Indicates the Z-axis coordinate of the reference satellite in the rectangular coordinate system; row m_1 Indicates the row index in the reference SAR image, col m_1 represents the column index in the reference SAR image; X m Indicates the X-axis value of the first coordinate; Y m Indicates the Y-axis value of the first coordinate; Z m represents the Z-axis value of the first coordinate; NR represents the near-edge slant distance of the observed scene; WS represents the range sampling interval; V x Indicates the X-axis velocity of the reference satellite; V y Indicates the Y-axis velocity of the reference satellite; V z represents the Z-axis speed of the reference satellite; λ represents the wavelength of the radar electromagnetic wave; R e represents the equatorial radius of the earth; h represents the surface elevation of the observation scene; R P represents the Earth's polar radius, f dc Indicates the Doppler frequency.

10. The method for geometric calibration of constellation SAR based on permanent scatterers according to claim 1, characterized in that: The correction formula is expressed as: R(row,col)=NR+WS·col+R cal (row,col); Where NR represents the near-edge slant range of the observed scene; WS represents the range sampling interval; col represents the SAR image row coordinate in the SAR image coordinate system; row represents the SAR image column coordinate in the SAR image coordinate system; NR+WS·col represents the radar slant range measurement value expression before correction; R cal (row,col) represents the slant distance measurement error correction; R(row,col) represents the corrected slant distance.

Images