A corner reflector centi-pixel extraction method

Abstract

The application discloses a corner reflector one-thousandth pixel extraction method and relates to the technical field of satellite measurement. The method comprises the following steps: step S1, acquiring original data, wherein the original data comprises geodetic coordinates of a corner reflector point, a SAR image and imaging parameters of the SAR image; step S2, constructing a range-Doppler (R-D) model, and calculating initial image coordinates (x0, y0) of the corner reflector point according to the R-D model by using the acquired original data; and step S3, based on the initial image coordinates of the corner reflector point acquired in step S2, solving pixel coordinates (x3, y3) of the corner reflector point in the SAR image. The application provides a method for extracting a corner reflector one-thousandth pixel, which is used for giving any corner reflector one-thousandth pixel level coordinates, and providing support for accurate extraction of a phase center of the corner reflector, geometric calibration and interferometric measurement calibration.

Description

Technical Field

[0001] This invention belongs to the field of satellite measurement technology, and in particular relates to a method for extracting one-thousandth of a pixel from a corner reflector. Background Technology

[0002] Synthetic Aperture Radar (SAR) satellites, with their unique all-weather, all-day Earth imaging capabilities, play an irreplaceable role in fields such as natural resource management and emergency disaster mitigation. To ensure the high quality and accuracy of SAR satellite imaging, geometric calibration using ground corner reflectors (CRs) is required after launch to achieve extremely high levels of planar positioning accuracy. Accurate extraction of corner reflector coordinates from SAR imagery is fundamental to research in geometric calibration, interferometric calibration, and deformation monitoring.

[0003] Corner reflectors, with their unique geometry and strong echo characteristics, appear as bright spots with high signal-to-noise ratios in SAR imagery, making them easy to identify and locate. However, despite some progress in corner reflector design research, such as the selection of reflector surface materials, design processes, manufacturing tolerances, and angular accuracy, research on precise extraction techniques for corner reflectors remains insufficient, particularly lacking extraction techniques capable of achieving an accuracy of one-thousandth of a pixel.

[0004] Currently, corner reflector extraction primarily relies on the contrast between the strong echo signal they generate in an image and the weaker echo of the surrounding terrain. Common methods include identifying and matching corner reflectors in an image using window templates and calculating correlation coefficients to determine their precise location. While this method is computationally efficient, its accuracy is typically limited to the pixel level. Another technique identifies reflector points by calculating the coherence coefficient and amplitude differences within a window surrounding the corner reflector. Although this method offers higher accuracy, processing the massive amounts of data makes it prohibitively expensive.

[0005] The coarse identification method for corner reflectors based on the satellite-ground position mathematical model, combined with the strong echo characteristics of corner reflectors, can determine the approximate range of corner reflectors. However, the errors in the SAR system imaging process can lead to excessive coarse identification offsets. The increased influence of surrounding ground features limits the accuracy of corner reflector point extraction.

[0006] Furthermore, existing technologies generally lack effective quality assessment and verification methods, making it difficult to accurately evaluate the accuracy of the extracted results. This can lead to significant calculation errors in high-precision topographic mapping.

[0007] Therefore, in order to overcome the shortcomings of the above-mentioned technologies, this invention proposes a method for extracting one-thousandth of a pixel from a corner reflector. The method aims to provide the geodetic coordinates of any corner reflector and calculate the one-thousandth pixel-level coordinates of the corner reflector through technology, providing strong technical support for the accurate extraction, geometric calibration, and interferometric calibration of the phase center of the corner reflector. Summary of the Invention

[0008] To solve the above-mentioned technical problems, the present invention provides a method for extracting one-thousandth of a pixel from a corner reflector, which can achieve a pixel positioning accuracy of one-thousandth of a pixel for corner reflectors.

[0009] To solve the above-mentioned technical problems, the present invention is achieved through the following technical solution:

[0010] This invention relates to a method for extracting one-thousandth of a pixel from a corner reflector, comprising the following steps:

[0011] Step S1: Obtain raw data, which includes the geodetic coordinates of the corner reflector point, SAR image, and imaging parameters of the SAR image;

[0012] Step S2: Construct a distance-Doppler (RD) model and use the acquired raw data to calculate the initial pixel coordinates (x0, y0) of the corner reflector point based on the RD model;

[0013] Step S3: Based on the initial pixel coordinates (x0, y0) of the corner reflector point obtained in step S2, calculate the pixel coordinates (x3, y3) of the corner reflector point in the SAR image.

[0014] As a preferred embodiment of the present invention, step S2 specifically includes the following sub-steps:

[0015] Step S2.1: Transform the geodetic coordinates of the corner reflector point to the geocentric rectangular coordinate system;

[0016] Step S2.2: Obtain the SAR satellite state vector at any time by fitting the satellite orbit with the SAR satellite state vector;

[0017] Step S2.3: Calculate the azimuth imaging time t of the corner reflector point in the image. Specifically, this is calculated based on the start time, the time sampling interval, and the row number of the corner reflector point in the SAR image. The formula is: t = t0 + t i ·i;

[0018] Step S2.4: Construct the RD model;

[0019] Step S2.5: Solve the RD model using Newton's iteration method to find the initial pixel coordinates (x0, y0) of the corner reflector point.

[0020] As a preferred embodiment of the present invention, the conversion formula used for the conversion from geodetic coordinates to geocentric rectangular coordinates in step S2.1 is as follows:

[0021]

[0022] As a preferred embodiment of the present invention, the satellite orbit fitting in step S2.2 employs polynomial fitting; specifically, the polynomial is:

[0023]

[0024] As a preferred embodiment of the present invention, the RD model in step S2.4 is as follows:

[0025] Distance equation:

[0026] Doppler equations:

[0027] Earth ellipsoid equation:

[0028] As a preferred embodiment of the present invention, step S3 specifically includes the following sub-steps:

[0029] Step S3.1: Using the initial pixel coordinates (x0, y0) of the corner reflector point obtained in step S2 as the center, select a complex matrix with a window size of N1×N1;

[0030] Step S3.2: Calculate the amplitude value of the N1×N1 complex matrix to obtain its amplitude value matrix; specifically, calculate the amplitude value according to the following formula:

[0031]

[0032] Where Re(i,j) represents the real part of the cell value (i,j) and Im(i,j) represents the imaginary part of the cell value (i,j);

[0033] Step S3.3: When the amplitude value sig(i,j) mentioned in step S3.2 is the largest, its corresponding pixel coordinates are (x1,y1).

[0034] Step S3.4: Using the pixel coordinates (x1, y1) mentioned in step S3.3 as the center, select a complex matrix with a window size of 2N2×2N2;

[0035] Step S3.5: Upsample the complex matrix with a window size of 2N2×2N2 by M1 times (M1<1024) to obtain the upsampled and reconstructed complex matrix;

[0036] Step S3.6: Calculate the amplitude value of the upsampled and reconstructed complex matrix to obtain its amplitude value matrix sig1(i,j);

[0037] Step S3.7: When the amplitude value described in step S3.6 is the largest, the corresponding pixel coordinates are (x2, y2).

[0038] Step S3.8: Perform M2 upsampling on the amplitude value of the row (i.e., distance direction) where pixel coordinate x2 is located to obtain a one-dimensional complex matrix Dx2; perform M2 upsampling on the amplitude value of the column (i.e., azimuth direction) where pixel coordinate y2 is located to obtain a one-dimensional complex matrix Dy2.

[0039] Step S3.9, calculate the magnitude matrix of the one-dimensional complex matrix Dx2 as follows: The magnitude matrix of the one-dimensional complex matrix Dy2 is

[0040] Step S3.10, calculate the amplitude value matrix. and The decibel (dB) values ​​are respectively and Specifically, the decibel value is calculated using the following formula:

[0041]

[0042] Where max(*) represents the maximum amplitude value in the amplitude matrix;

[0043] Step S3.11, retrieve the distance-direction db value All the dB values ​​from the left 3dB to the right 3dB position constitute a new matrix.

[0044] Step S3.12, based on the distance to the db value Trigonometric function fitting is performed on the pixel coordinates and db values;

[0045] Step S3.13, find Find the coordinates of the peak point and the coordinates of the fitted value, and determine whether the absolute value of the difference between the coordinates x3 of the peak point and the coordinates x'3 of the fitted value reaches the order of 1 / 1000. Output the distance to the coordinates x3 of the peak point.

[0046] Step S3.14, azimuth db value Repeat steps S3.11-S3.13 to output the coordinates y3 corresponding to the azimuth peak point.

[0047] As a preferred embodiment of the present invention, step S3.5 specifically includes the following sub-steps:

[0048] Step S3.5.1: Perform a two-dimensional fast Fourier transform on the complex matrix obtained in S3.4 to convert it to the frequency domain;

[0049] Step S3.5.2: Upsample the complex matrix transferred to the frequency domain by M1 times to obtain the upsampled complex matrix.

[0050] Step S3.5.3: Reconstruct the upsampled complex matrix using a two-dimensional inverse Fourier transform. The reconstructed complex matrix D2[M1×2N2,M1×2N2] is obtained.

[0051] As a preferred embodiment of the present invention, step S3.52 specifically includes the following sub-steps:

[0052] Step S3.5.2.1: Calculate the amplitude value of the complex matrix D2[M1×2N2,M1×2N2] obtained in step S3.5.1. The sum of the amplitude values ​​in each column is then calculated to obtain the matrix tmp. row (i,1);

[0053] In step S3.5.2.2, locate step S3.5.2.1 tmp. row The row number i1 where the minimum value is located in (i,1);

[0054] Step S3.5.2.3: Construct a two-dimensional zero matrix O1[(M1-1)×2N2,2N2];

[0055] Step S3.5.2.4: The two-dimensional zero matrix obtained in step S3.5.2.3

[0056] O1[(M1-1)×2N2,2N2] is inserted into sig. D2 The (i,j)th row is obtained, and the matrix is ​​obtained.

[0057] Step S3.5.2.5, for the matrix in step S3.5.2.4 The sum of the amplitude values ​​in each row is used to obtain the matrix tmp. col (1,j);

[0058] Step S3.5.2.6: Locate the matrix from step S3.5.2.5. The column number j1 where the minimum value is located;

[0059] Step S3.5.2.7: Construct a two-dimensional zero matrix O2[M1×2N2,(M1-1)×2N2];

[0060] Step S3.5.2.8: The two-dimensional zero matrix obtained in step S3.5.2.7

[0061] O2[M1×2N2,(M1-1)×2N2] Insert Find the (j1+1)th column and obtain the matrix. The above eight sub-steps complete the first M1-fold upsampling process.

[0062] As a preferred embodiment of the present invention, step S3.8 specifically includes the following sub-steps:

[0063] Step S3.8.1, determine the upsampling factor M2; specifically, M2 is determined according to the following formula:

[0064] M2 = 1024 / M1

[0065] Step S3.8.2: Perform M2 upsampling on the amplitude value of the row (i.e., distance direction) where pixel coordinate x2 is located to obtain a one-dimensional complex matrix Dx2, and perform M2 upsampling on the amplitude value of the column (i.e., azimuth direction) where pixel coordinate y2 is located to obtain a one-dimensional complex matrix Dy2.

[0066] As a preferred embodiment of the present invention, step S3.8.2 specifically includes the following sub-steps:

[0067] Step S3.8.2.1: Construct a one-dimensional zero matrix and determine the size N3 of the one-dimensional zero matrix O; specifically, determine N3 according to equation (10).

[0068] N3=(2N2×M1)×M2

[0069] =1024*2N2

[0070] Step S3.8.2.2: Perform spectrum centering on the amplitude value of the row (i.e., distance direction) where pixel coordinate x2 is located, then perform fast Fourier transform on the amplitude value after spectrum centering, and finally perform spectrum centering again; this process is to ensure that the low-frequency part of the spectrum is located at the beginning of the array and the high-frequency part is located at the end of the array.

[0071] Step S3.8.2.3: Assign the result obtained in step S3.8.2.2 to the zero matrix O3. Within the range, we obtain the assigned matrix Dx2;

[0072] In step S3.8.2.4, the Dx2 after being assigned a value in step S3.8.2.3 is first subjected to spectrum centering, then to inverse fast Fourier transform, and finally to spectrum centering again; this operation ensures that the format of the final result is consistent with the original signal.

[0073] Step S3.8.2.5: The above four sub-steps complete the process of upsampling the amplitude value of the row where pixel coordinate x2 is located (i.e., the distance direction) by M2; the amplitude value of the azimuth direction is resampled and steps S3.8.2.1-S3.8.2.5 are repeated.

[0074] The present invention has the following beneficial effects:

[0075] 1. Significantly improves the positioning accuracy of corner reflectors;

[0076] 2. Improve SAR satellite imaging quality: By accurately extracting the coordinates of corner reflectors, this technology can improve the imaging quality of SAR satellites, which is crucial for applications in natural resource management, emergency response and disaster reduction, and other fields.

[0077] 3. Enhanced positioning accuracy: This method, through innovative algorithms and data processing technology, can provide higher positioning accuracy and reduce errors caused by manual pricking.

[0078] 4. Supports high-precision applications: This technology is particularly suitable for applications requiring high-precision corner reflector positioning, such as satellite geometric calibration, interferometric calibration, and deformation monitoring;

[0079] 5. Reduced error impact: Compared with existing technologies, this invention has significant advantages in reducing the impact of corner reflector extraction errors, especially in applications such as high-precision topographic mapping;

[0080] 6. Improve data processing efficiency: Not only does it improve the accuracy of the calculation, but the method has also been optimized in terms of data processing efficiency, making it more efficient than existing technologies.

[0081] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0082] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0083] Figure 1 This is a flowchart of the method for extracting one-thousandth of a pixel from a mid-angle reflector;

[0084] Figure 2 This is a schematic diagram of the geometric relationship for synthetic aperture radar satellite positioning;

[0085] Figure 3 This is a schematic diagram of a corner reflector on a SAR image. Detailed Implementation

[0086] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0087] like Figure 1 The image shows a method for extracting one-thousandth of a pixel from a corner reflector, which includes the following steps:

[0088] Step S1: Obtain raw data, which includes the coordinates of the corner reflector point, SAR image, and imaging parameters of the SAR image;

[0089] Step S2: Construct a range-Doppler (RD) model and use the acquired raw data to calculate the initial image coordinates (x0, y0) of the corner reflector point based on the RD model;

[0090] like Figure 2 As shown, specifically, step S2 includes the following:

[0091] Step S2.1: Convert the geodetic coordinates of the corner reflector point to a geocentric rectangular coordinate system; specifically, the conversion formula from geodetic coordinates to a geocentric rectangular coordinate system is as follows:

[0092]

[0093] Among them, (X) P Y P Z P () represents the geocentric rectangular coordinates of the corner reflector point, N is the radius of curvature of the prime meridian, (B, L, H) are the latitude, longitude and geodetic height of the corner reflector point respectively, and e is the first eccentricity of the reference ellipsoid.

[0094] Step S2.2: By fitting the satellite orbit with the SAR satellite state vector, the SAR satellite state vector at any given time is obtained. The position and velocity vectors of the satellite sensor at any given time can be fitted using a polynomial; specifically, the polynomial is:

[0095]

[0096] Among them, (X) s Y s Z s (V) is the position vector to be fitted. xs Yys Z zs () represents the velocity vector to be fitted. The coefficient matrix of the position-fitting polynomial. is the coefficient matrix of the velocity fitting polynomial, and t is the azimuth imaging time of the upper corner reflector point in the image;

[0097] Step S2.3, calculate the azimuth imaging time of the upper corner reflector point in the image; specifically, calculate t according to equation (4):

[0098] t = t0 + t i ·i (4)

[0099] Where t0 is the azimuth start time of the SAR image, t i denoted as the azimuth time sampling interval of the SAR image, and i is the row number of the corner reflector point on the SAR image.

[0100] Step S2.4, construct the RD model, as shown below:

[0101] Distance equation:

[0102] Doppler equations:

[0103] Earth ellipsoid equation:

[0104] Among them, R gc =[X g Y g Z g ] T R represents the position vector of a point on the ground. sc =[X s Y s Z s ] T This represents the position vector of the satellite sensor, r0 is the initial slant range in the range direction of the SAR image, and m j Here, j is the range resolution of the SAR image, j is the column number of the corner reflector point on the SAR image, and f is the range resolution of the SAR image. d R is the Doppler frequency. e R represents the semi-major axis of the Earth's ellipsoid. p H represents the minor semi-axis of the Earth's ellipsoid. t This represents the elevation value of the ground point.

[0105] Step S2.5: Solve the RD model constructed in step S2.4 using Newton's iteration method to solve for the initial pixel coordinates of the corner reflector points;

[0106] Step S3: Based on the initial image coordinates of the corner reflector point obtained in step S2, the pixel coordinates (x3, y3) of the corner reflector point in the SAR image are calculated using the one-thousandth pixel extraction method.

[0107] like Figure 3 As shown, specifically, step S3 includes the following:

[0108] Step S3.1: Using the initial pixel coordinates (x0, y0) of the corner reflector point obtained in step S2 as the center, select a complex matrix with a window size of N1×N1;

[0109] Step S3.2: Calculate the amplitude value of the N1×N1 complex matrix to obtain its amplitude value matrix; specifically, calculate the amplitude value according to equation (8):

[0110]

[0111] Where Re(i,j) represents the real part of the cell value (i,j) and Im(i,j) represents the imaginary part of the cell value (i,j);

[0112] Step S3.3: When the amplitude value sig(i,j) mentioned in step S3.2 is the largest, its corresponding pixel coordinates are (x1,y1).

[0113] Step S3.4: Using the pixel coordinates (x1, y1) mentioned in step S3.3 as the center, select a complex matrix D1[2N2, 2N2] with a window size of 2N2×2N2;

[0114] Step S3.5: Upsample the complex matrix with a window size of 2N2×2N2 by M1 times (M1<1024) to obtain the upsampled and reconstructed complex matrix D2[M1×2N2,M1×2N2];

[0115] Step S3.5.1: Perform a two-dimensional fast Fourier transform on the complex matrix D2[M1×2N2,M1×2N2] obtained in S3.4 to transform it into the frequency domain;

[0116] Step S3.5.2: Upsample the complex matrix transferred to the frequency domain by M1 times to obtain the upsampled complex matrix.

[0117] Step S3.5.2.1: Calculate the amplitude value of the complex matrix D2[M1×2N2,M1×2N2] obtained in step S3.5.1. The sum of the amplitude values ​​in each column is then calculated to obtain the matrix tmp. row (i,1);

[0118] In step S3.5.2.2, locate step S3.5.2.1 tmp. row The row number i1 where the minimum value is located in (i,1);

[0119] Step S3.5.2.3: Construct a two-dimensional zero matrix O1[(M1-1)×2N2,2N2];

[0120] Step S3.5.2.4: The two-dimensional zero matrix obtained in step S3.5.2.3

[0121] O1[(M1-1)×2N2,2N2] Insert Find the (i1+1)th row in the matrix and obtain the matrix.

[0122] Step S3.5.2.5, for the matrix in step S3.5.2.4 The sum of the amplitude values ​​in each row is used to obtain the matrix tmp. col (1,j);

[0123] Step S3.5.2.6: Locate the matrix from step S3.5.2.5. The column number j1 where the minimum value is located;

[0124] Step S3.5.2.7: Construct a two-dimensional zero matrix O2[M1×2N2,(M1-1)×2N2];

[0125] Step S3.5.2.8: Insert the two-dimensional zero matrix O2[M1×2N2,(M1-1)×2N2] obtained in step S3.5.2.7. Find the (j1+1)th column and obtain the matrix. The above eight sub-steps complete the first M1-fold upsampling process;

[0126] Step S3.5.3, the matrix obtained in step S3.5.2 Using the two-dimensional inverse Fourier transform, we obtain the complex matrix D2[M1×2N2,M1×2N2];

[0127] Step S3.6: Calculate the amplitude value of the upsampled complex matrix D2[M1×2N2,M1×2N2] to obtain its amplitude value matrix sig1(i,j); specifically, calculate its amplitude value according to equation (8);

[0128] Step S3.7: When the amplitude value described in step S3.6 is the largest, the corresponding pixel coordinates are (x2, y2).

[0129] Step S3.8: Perform M2 upsampling on the amplitude value of the row (i.e., distance direction) where pixel coordinate x2 is located to obtain a one-dimensional complex matrix Dx2; perform M2 upsampling on the amplitude value of the column (i.e., azimuth direction) where pixel coordinate y2 is located to obtain a one-dimensional complex matrix Dy2.

[0130] Step S3.8.1, determine the upsampling factor M2; specifically, determine M2 according to equation (9):

[0131] M2 = 1024 / M1 (9)

[0132] Step S3.8.2: Perform M2 upsampling on the amplitude value of the row (i.e., distance direction) where pixel coordinate x2 is located to obtain a one-dimensional complex matrix Dx2; perform M2 upsampling on the amplitude value of the column (i.e., azimuth direction) where pixel coordinate y2 is located to obtain a one-dimensional complex matrix Dy2.

[0133] Step S3.8.2.1: Construct a one-dimensional zero matrix and determine the size N3 of the one-dimensional zero matrix O; specifically, determine N3 according to equation (10).

[0134]

[0135] Step S3.8.2.2: Perform spectrum centering on the amplitude value of the row (i.e., distance direction) where pixel coordinate x2 is located, then perform fast Fourier transform on the amplitude value after spectrum centering, and finally perform spectrum centering again; this process is to ensure that the low-frequency part of the spectrum is located at the beginning of the array and the high-frequency part is located at the end of the array.

[0136] Step S3.8.2.3: Assign the result obtained in step S3.8.2.2 to the zero matrix O3. Within the range, we obtain the assigned matrix Dx2;

[0137] In step S3.8.2.4, the Dx2 after being assigned a value in step S3.8.2.3 is first subjected to spectrum centering, then to inverse fast Fourier transform, and finally to spectrum centering again; this operation ensures that the format of the final result is consistent with the original signal.

[0138] Step S3.8.2.5: The above four sub-steps complete the process of upsampling the amplitude value of the row where pixel coordinate x2 is located (i.e., the distance direction) by M2; the amplitude value of the azimuth direction is resampled and steps S3.8.2.1-S3.8.2.5 are repeated.

[0139] Step S3.9, calculate the magnitude matrix of the one-dimensional complex matrix Dx2 as follows: The magnitude matrix of the one-dimensional complex matrix Dy2 is

[0140] Step S3.10, calculate the amplitude value matrix. and The decibel (dB) values ​​are respectively and Specifically, the decibel value is calculated using the following formula:

[0141]

[0142] Where max(*) represents the maximum amplitude value in the amplitude matrix;

[0143] Step S3.11, retrieve the distance-direction db value All the dB values ​​from the left 3dB to the right 3dB position constitute a new matrix.

[0144] Step S3.12, based on the distance to the db value Trigonometric function fitting is performed on the pixel coordinates and db values;

[0145] Step S3.13, find Find the coordinates of the peak point and the coordinates of the fitted value, and determine whether the absolute value of the difference between the coordinates x3 of the peak point and the coordinates x'3 of the fitted value reaches the order of 1 / 1000. Output the coordinates x3 of the actual value of the peak point.

[0146] Step S3.14, azimuth db value Repeat steps S3.11-S3.13 to output the coordinates y3 corresponding to the azimuth peak point.

[0147] In the description of this specification, references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0148] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for extracting one-thousandth of a pixel from a corner reflector, characterized in that, Includes the following steps: Step S1: Obtain raw data, which includes the geodetic coordinates of the corner reflector point, SAR image, and imaging parameters of the SAR image; Step S2: Construct a distance-Doppler model and use the acquired raw data to calculate the initial pixel coordinates (x0, y0) of the corner reflector point based on the distance-Doppler model; Step S3: Based on the initial pixel coordinates (x0, y0) of the corner reflector point obtained in step S2, calculate the pixel coordinates (x3, y3) of the corner reflector point in the SAR image: Step S3 specifically includes the following sub-steps: Step S3.1: Using the initial pixel coordinates (x0, y0) of the corner reflector point obtained in step S2 as the center, select a complex matrix with a window size of N1×N1; Step S3.2: Calculate the amplitude value of the N1×N1 complex matrix to obtain its amplitude value matrix; specifically, calculate the amplitude value according to the following formula: Where Re(i,j) represents the real part of the cell value (i,j) and Im(i,j) represents the imaginary part of the cell value (i,j); Step S3.3: When the amplitude value sig(i,j) mentioned in step S3.2 is the largest, its corresponding pixel coordinates are (x1,y1). Step S3.4: Using the pixel coordinates (x1, y1) mentioned in step S3.3 as the center, select a complex matrix with a window size of 2N2×2N2; Step S3.5: Upsample the complex matrix with a window size of 2N2×2N2 by M1 times (M1<1024) to obtain the upsampled and reconstructed complex matrix; Step S3.6: Calculate the amplitude value of the upsampled and reconstructed complex matrix to obtain its amplitude value matrix sig1(i,j); Step S3.7: When the amplitude value described in step S3.6 is the largest, the corresponding pixel coordinates are (x2, y2). Step S3.8: Perform M2 upsampling on the amplitude value of the row containing pixel coordinate x2 to obtain a one-dimensional complex matrix Dx2; perform M2 upsampling on the amplitude value of the column containing pixel coordinate y2 to obtain a one-dimensional complex matrix Dy2. Step S3.9, calculate the magnitude matrix of the one-dimensional complex matrix Dx2 as follows: The magnitude matrix of the one-dimensional complex matrix Dy2 is Step S3.10, calculate the amplitude value matrix. and The decibel values ​​are respectively and Specifically, the decibel value is calculated using the following formula: Where max(*) represents the maximum amplitude value in the amplitude matrix; Step S3.11, retrieve the distance-direction db value All the dB values ​​from the left 3dB to the right 3dB position constitute a new matrix. Step S3.12, based on the distance to the db value Trigonometric function fitting is performed on the pixel coordinates and db values; Step S3.13, find Find the coordinates of the peak point and the coordinates of the fitted value, and determine whether the absolute value of the difference between the coordinates x3 of the peak point and the coordinates x'3 of the fitted value reaches the order of 1 / 1000. Output the distance to the coordinates x3 of the peak point. Step S3.14, azimuth db value Repeat steps S3.11-S3.13 to output the coordinates y3 corresponding to the azimuth peak point.

2. The method for extracting one-thousandth of a pixel from a corner reflector according to claim 1, characterized in that, Step S2 specifically includes the following sub-steps: Step S2.1: Transform the geodetic coordinates of the corner reflector point to the geocentric rectangular coordinate system; Step S2.2: Fit the satellite orbit using the SAR satellite state vector to obtain the SAR satellite state vector at any time. Step S2.3: Calculate the azimuth imaging time t of the corner reflector point in the image. Specifically, this is calculated based on the start time, the time sampling interval, and the row number of the corner reflector point in the SAR image. The formula is: t = t0 + t i ·i; Step S2.4, construct the distance-Doppler model; Step S2.5: Solve the distance-Doppler model using Newton's iteration method to find the initial pixel coordinates (x0, y0) of the corner reflector point.

3. The method for extracting one-thousandth of a pixel from a corner reflector according to claim 2, characterized in that, The conversion formula used in step S2.1 to transform geodetic coordinates to geocentric rectangular coordinates is as follows:

4. The method for extracting one-thousandth of a pixel from a corner reflector according to claim 2, characterized in that, The satellite orbit fitting in step S2.2 employs polynomial fitting; specifically, the polynomial is:

5. The method for extracting one-thousandth of a pixel from a corner reflector according to claim 2, characterized in that, The distance-Doppler model in step S2.4 is as follows: Distance equation: Doppler equations: Earth's ellipsoid equation:

6. The method for extracting one-thousandth of a pixel from a corner reflector according to claim 1, characterized in that, Step S3.5 specifically includes the following sub-steps: Step S3.5.1: Perform a two-dimensional fast Fourier transform on the complex matrix obtained in S3.4 to convert it to the frequency domain; Step S3.5.2: Upsample the complex matrix transferred to the frequency domain by M1 times to obtain the upsampled complex matrix. Step S3.5.3: Reconstruct the upsampled complex matrix using a two-dimensional inverse Fourier transform. The reconstructed complex matrix D2[M1×2N2,M1×2N2] is obtained.

7. The method for extracting one-thousandth of a pixel from a corner reflector according to claim 6, characterized in that, Step S3.5.2 specifically includes the following sub-steps: Step S3.5.2.1: Calculate the amplitude value of the complex matrix D2[M1×2N2,M1×2N2] obtained in step S3.5.

1. The sum of the amplitude values ​​in each column is then calculated to obtain the matrix tmp. row (i,1); In step S3.5.2.2, locate step S3.5.2.1 tmp. row The row number i1 where the minimum value is located in (i,1); Step S3.5.2.3: Construct a two-dimensional zero matrix O1[(M1-1)×2N2,2N2]; Step S3.5.2.4: Insert the two-dimensional zero matrix O1[(M1-1)×2N2,2N2] obtained in step S3.5.2.

3. Find the (i1+1)th row in the matrix and obtain the matrix. Step S3.5.2.5, for the matrix in step S3.5.2.4 The sum of the amplitude values ​​in each row is used to obtain the matrix tmp. col (1,j); Step S3.5.2.6: Locate the matrix from step S3.5.2.

5. The column number j1 where the minimum value is located; Step S3.5.2.7: Construct a two-dimensional zero matrix O2[M1×2N2,(M1-1)×2N2]; Step S3.5.2.8: Insert the two-dimensional zero matrix O2[M1×2N2,(M1-1)×2N2] obtained in step S3.5.2.

7. Find the (j1+1)th column and obtain the matrix. The above eight sub-steps complete the first M1-fold upsampling process.

8. The method for extracting one-thousandth of a pixel from a corner reflector according to claim 1, characterized in that, Step S3.8 specifically includes the following sub-steps: Step S3.8.1, determine the upsampling factor M2; specifically, M2 is determined according to the following formula: M2 = 1024 / M1 Step S3.8.2: Perform M² upsampling on the amplitude value of the row containing pixel coordinate x2 to obtain a one-dimensional complex matrix Dx2; perform M² upsampling on the amplitude value of the column containing pixel coordinate y2 to obtain a one-dimensional complex matrix Dy.

2. .

9. The method for extracting one-thousandth of a pixel from a corner reflector according to claim 8, characterized in that, Step S3.8.2 specifically includes the following sub-steps: Step S3.8.2.1: Construct a one-dimensional zero matrix and determine the size N3 of the one-dimensional zero matrix O; specifically, determine N3 according to the following formula: N3=(2N2×M1)×M2=1024*2N2 Step S3.8.2.2: Perform spectrum centering on the amplitude value of the row containing pixel coordinate x2, then perform Fast Fourier Transform on the spectrum-centered amplitude value, and finally perform spectrum centering again; this process is to ensure that the low-frequency part of the spectrum is located at the beginning of the array and the high-frequency part is located at the end of the array. Step S3.8.2.3: Assign the result obtained in step S3.8.2.2 to the zero matrix O3. Within the range, we obtain the assigned matrix Dx2; In step S3.8.2.4, the Dx2 after being assigned a value in step S3.8.2.3 is first subjected to spectrum centering, then subjected to inverse fast Fourier transform, and finally subjected to spectrum centering again. This operation ensures that the format of the final result is consistent with the original signal; Step S3.8.2.5: The above four sub-steps complete the process of upsampling the amplitude value of the row where pixel coordinate x2 is located by M2; the amplitude value of the azimuth direction is resampled and steps S3.8.2.1-S3.8.2.5 are repeated.

Images