CN113902645B - A Method for Acquiring RPC Correction Parameters of Spaceborne SAR Images Based on Inverse RD Positioning Model - Google Patents

Abstract

A satellite-borne SAR image RPC correction parameter acquisition method based on an inverse RD positioning model belongs to the technical field of microwave remote sensing SAR image geometric correction, and solves the problem of how to design a complete satellite-borne SAR image RPC correction parameter generation method based on real terrain positioning for a satellite-borne SAR general geometric correction model; and performing least square accurate solution on the RPC model equation by adopting a spectrum correction iteration method to solve the problem of equation morbidity caused by excessive unknown parameters and strong correlation of the equation, so that the accuracy of the SAR image general geometric correction is ensured and improved.

Description

A Method for Acquiring RPC Correction Parameters of Spaceborne SAR Images Based on Inverse RD Positioning Model

technical field

The invention belongs to the technical field of geometric correction of microwave remote sensing SAR images, and relates to a method for acquiring RPC correction parameters of spaceborne SAR images based on an inverse RD positioning model.

Background technique

As more and more advanced system spaceborne SAR (Synthetic Aperture Radar) satellites are put into operation and put into use, they will play an increasingly important role in scientific research and economic construction. Due to the SAR oblique range imaging mechanism, it needs to be corrected to the ground coordinates to facilitate quantitative geometric analysis and image interpretation, so as to better serve applications such as terrain surveying, resource surveys, and geological disaster monitoring and early warning. Therefore, the geometric correction of high-precision spaceborne SAR is of great significance. At present, there are mainly two methods based on SAR rigorous imaging model and sensor-independent general geometric model. The image RPC (Rational Polynomial Coefficients, Rational Polynomial Coefficients) model is a sensor-independent general geometric correction model. By expressing pixel coordinates as ratio polynomials with ground space coordinates as independent variables, the orthographic projection of SAR images is realized. Since the RPC correction model avoids the detailed information of SAR imaging sensors and platforms, and in the case of reasonable selection of model input, the model fitting accuracy is close to the SAR rigorous imaging model, which has the advantages of simplicity, generality, confidentiality, and efficiency. Therefore, the sensor- and platform-independent RPC model is widely used in geometric correction of spaceborne SAR images.

Regarding the solution of RPC model parameters, since the general model contains seventy-eight rational polynomial coefficients (referred to as RPCs), a sufficient and reasonably distributed array of control points is required for iterative calculation. Regarding the generation of control point arrays, there are usually methods of generating virtual control points on multiple elevation surfaces based on the SAR rigorous imaging model and methods of directly using a large number of ground truth control points; It is often difficult to achieve in product generation; in addition, due to the large number of unknown RPCs and the strong correlation between them, the problem of equation ill-conditioning often occurs in the equation solution and needs to be considered more.

In the prior art, the Chinese invention patent application with the publication number CN103218780A and the publication date of July 24, 2013 "Orthorectification Method for Uncontrolled Spaceborne SAR Images Based on Inverse RD Positioning Model" discloses a method based on inverse RD positioning The spaceborne SAR image orthorectification method of the model, the present invention takes it as the technical premise, and provides a general RPC correction model parameter generation method in combination with the ground true height positioning, so as to serve the SAR general geometric correction with technical confidentiality characteristics.

Contents of the invention

The purpose of the present invention is to solve the deficiencies of the prior art, and to design a complete RPC correction parameter generation method for spaceborne SAR images based on real terrain positioning for the general geometric correction model of spaceborne SAR.

The present invention solves the above technical problems through the following technical solutions:

1, a kind of RPC correction parameter acquisition method based on the inverse RD positioning model of spaceborne SAR image, is characterized in that, comprises the following steps:

S1. Acquisition of SAR imaging parameters and satellite orbit data: acquisition of spaceborne SAR system imaging and image parameters, satellite orbit data and ground elevation data;

S2. Polynomial fitting of satellite orbit data: according to the satellite position and velocity data at multiple moments in the SAR imaging system and image parameters, fit the function curve of the satellite orbit with respect to the azimuth pixel;

S3. Positioning of the four corners of the SAR image: Construct the positioning equation of the slant range-Doppler-Earth ellipsoid model, and iteratively solve the ground coordinates of the four corners of the SAR image, which are the first azimuth moment and the corresponding time of the SAR image respectively. The W6S84 coordinates of the closest slant distance and the farthest slant distance corresponding to the last azimuth time respectively determine the range of the ground area of the SAR image;

S4. Ground three-dimensional grid establishment: according to the geographical coordinates of the four corners and the global ground DEM (Digital Elevation Model, digital elevation model) data, the elevation data information within the geographic range of the required image is obtained; it is evenly divided into M*N Grid, calculate the three-dimensional coordinate information of the highest, lowest and intermediate values of elevation in each grid, and obtain 3M*N ground control point data in total; then perform Cartesian coordinate conversion to obtain new control points under the geocentric coordinate system position vector group;

S5. Inverse RD (Range Doppler, slant range Doppler) image positioning of the ground control point array: construct the slant range error equation and the Doppler frequency error equation, and inverse each control point according to the satellite orbit pixel fitting curve The RD positioning equation is solved iteratively by Newton to obtain the image pixel position corresponding to each control point;

S6. Compliance prediction and regularization processing of the control point array: perform compliance prediction on the obtained image pixel position and its corresponding ground three-dimensional coordinate information, remove points that are not within the image range; perform data regularization processing , to normalize the ground coordinates and image coordinates to reduce rounding errors in the calculation process;

S7. RPC model equation construction and least squares solution: establish a ratio polynomial between ground coordinates and image coordinates, and deform it to obtain a nonlinear model about RPCs; for the ill-conditioned problem of the equation coefficient matrix, use the least squares based on the spectral correction iterative method Multiply the solution to obtain an unbiased estimate of the RPCs parameters;

S8. Model accuracy verification: use the solved RPCs model parameters to calculate the image coordinates corresponding to the checkpoints, and evaluate the image square accuracy of the RPC model parameters by making a difference with the calculation results of the rigorous imaging geometric model.

The innovations of the technical solution of the present invention are: 1) According to the SAR rigorous imaging model, the inverse RD positioning algorithm and real terrain information are used to uniformly sample the real space of the imaging area, and an array of image ground control points with reasonable distribution is generated for the input of the RPC equation model ; 2) The least squares exact solution of the RPC model equations is carried out by using the spectrum correction iterative method to solve the ill-conditioned problem of the equations caused by too many unknowns and strong correlations, so that the accuracy of the general geometric correction of SAR images can be guaranteed and improved.

As a further improvement of the technical solution of the present invention, the spaceborne SAR system imaging and image parameters described in step S1 include: radar wavelength, image size, imaging Doppler frequency information, imaging time T ₀ of the first line of the image, the first line of the image Pixel pulse delay t ₀ or image short distance R ₀ , pulse repetition frequency PRF, and range sampling frequency RSR; the satellite orbit data includes: position and velocity vector data at the sampling point.

As a further improvement of the technical solution of the present invention, the method for fitting the function curve of the satellite orbit about the azimuth pixel described in step S2 is as follows:

S21, the polynomial fitting model about the position is as follows:

Among them, m represents a pixel, x _s , y _s , and z _s represent the position of the satellite orbit corresponding to the pixel in the x direction, y direction, and z direction, and α _i , β _i , and γ _i are the satellite orbits to be fitted about The coefficients of the function curve polynomial of the pixels in the azimuth direction, i=1, 2, 3, 4, 5;

S22, the polynomial fitting model about speed is as follows:

Among them, vx _s , vy _s , vz _s represent the velocity of the satellite orbit corresponding to pixel m at positions x _s , y _s , z _s .

As a further improvement of the technical solution of the present invention, the slant range-Doppler-Earth ellipsoid model positioning equation described in step S3 is as follows:

其中，R_e、R_p分别为椭球体的半长轴、半短轴；h为目标点的高程，通常可设为成像区域的平均高程；

为目标点t的位置矢量坐标；1) Earth ellipsoid model equation:

Among them, R _e and R _p are the semi-major axis and semi-minor axis of the ellipsoid respectively; h is the elevation of the target point, which can usually be set as the average elevation of the imaging area;

is the position vector coordinate of the target point t;

其中，

是卫星位置矢量，

R为卫星到目标的斜距，即

2) Slope distance equation:

in,

is the satellite position vector,

R is the slant distance from the satellite to the target, that is

其中，λ为雷达波长，

卫星速度矢量，

为目标速度矢量，采用ECR(Earth CenteredRotation，地心坐标系)坐标系统定义时其值可简化为0。3) Doppler frequency equation:

where λ is the radar wavelength,

satellite velocity vector,

is the target velocity vector, and its value can be simplified to 0 when defined by the ECR (Earth Centered Rotation, earth-centered coordinate system) coordinate system.

As a further improvement of the technical solution of the present invention, the formula for Cartesian coordinate transformation described in step S4 is:

为地球第一偏心率。in,

is the first eccentricity of the earth.

As a further improvement of the technical solution of the present invention, the construction of the slant range error equation and the Doppler frequency error equation described in step S5, according to the satellite orbit pixel fitting curve, performs the inverse RD positioning equation Newton iterative solution for each control point, The method of obtaining the image pixel position corresponding to each control point is as follows:

(1) According to the satellite fitting curve, construct the slant range error equation and Doppler frequency error equation as follows:

Among them, f _s is the sampling frequency RSR, C is the speed of light, λ is the radar wavelength, and f _d (n) is the Doppler center frequency corresponding to the pixel. Combined into f _d (n) = f ₀ +f ₁ n+f ₂ n ² ;

(2) The Newton iterative method is used to solve the inverse slope range-Doppler equation about the pixel position (m, n) as follows:

Among them, (m ₀ , n ₀ ) is the preset initial value, then the kth group of solutions in the iterative process can be obtained as follows:

in,

可简单通过轨道拟合公式获得；Among them, the partial conduction coefficient

It can be obtained simply by orbital fitting formula;

(3) Through multiple iterative loops, when the difference accuracy of two adjacent solutions reaches the convergence value, the iterative solution is stopped and the target pixel position is obtained.

As a further improvement of the technical solution of the present invention, the data regularization process described in step S6, the method of normalizing the ground coordinates and image coordinates is as follows:

Standardize the image coordinates and ground coordinates of 3M*N control point arrays, let (P, L, H) be the ground coordinates after the ground three-dimensional coordinates are standardized, and (X, Y) be the images after the normalized row and column values of the image pixel positions Coordinates, the standardized formula is as follows:

和

为地面坐标的标准化参数，即偏移值和比例值；m_off，m_scale，n_offf和n_scale为像素坐标的标准化参数。in,

with

is the normalization parameter of ground coordinates, that is, offset value and scale value; m _off , m _scale , n _offff and n _scale are normalization parameters of pixel coordinates.

As a further improvement of the technical solution of the present invention, the calculation formulas of the offset value and the ratio value are as follows:

m _scale = max(|m(i, j) _max -m _off |, |m(i, j) _min -m _off |);

n _scale = max(|n(i, j) _max -n _off |, |n(i, j) _min -n _off |).

As a further improvement of the technical solution of the present invention, the ratio polynomial between the ground coordinates and the image coordinates described in the step S7 is established, and the deformation method for obtaining a nonlinear model about RPCs is as follows:

(1) Establish a ratio polynomial between ground coordinates and image coordinates, and associate the ground point coordinates (P, L, H) with their corresponding pixel point coordinates (X, Y) with a ratio polynomial:

Among them, Num _line (P, L, H) = a ₀ +a ₁ L+a ₂ P+a ₃ H+a ₄ LP+a ₅ LH+a ₆ PH+a ₇ L ² +a ₈ P ² +a ₉ H ² +a ₁₀ PLH+a ₁₁ L ³ +a ₁₂ LP ² +a ₁₃ LH ² +a ₁₄ L ² P+a ₁₅ P ³ +a ₁₆ PH ² +a ₁₇ L ² H+a ₁₈ P ² H +a ₁₉ H ³ ;

Den _line (P, L, H)＝b ₀ +b ₁ L+b ₂ P+b ₃ H+b ₄ LP+b ₅ LH+b ₆ PH+b ₇ L ² +b ₈ P ² +b ₉ H ² +b ₁₀ PLH+b ₁₁ L ³ +b ₁₂ LP ² +b ₁₃ LH ² +b ₁₄ L ² P+b ₁₅ P ³ +b ₁₆ PH ² +b ₁₇ L ² H+b ₁₈ P ² H+b ₁₉ H ³ ;

Num _pixel (P, L, H)＝c ₀ +c ₁ L+c ₂ P+c ₃ H+c ₄ LP+c ₅ LH+c ₆ PH+c ₇ L ² +c ₈ P ² +c ₉ H ² +c ₁₀ PLH+c ₁₁ L ³ +c ₁₂ LP ² +c ₁₃ LH ² +c ₁₄ L ² P+c ₁₅ P ³ +c ₁₆ PH ² +c ₁₇ L ² H+c ₁₈ P ² H+c ₁₉ H ³ ;

Den _pixel (P, L, H)＝d ₀ +d ₁ L+d ₂ P+d ₃ H+d ₄ LP+d ₅ LH+d ₆ PH+d ₇ L ² +d ₈ P ² +d ₉ H ² +d ₁₀ PLH+d ₁₁ L ³ +d ₁₂ LP ² +d ₁₃ LH ² +d ₁₄ L ² P+d ₁₅ P ³ +d ₁₆ PH ² +d ₁₇ L ² H+d ₁₈ P ² H+d ₁₉ H ³ ;

(a _i , b _i , _ci , d _i , i=0, 1, ..., 19) in the above formula are the coefficients RPCs of the rational function model, and the comparison value polynomial is transformed to obtain the rational function model for each RPCs Linear model, and express the error equation listed by the control points in matrix form;

F _Y =Num _line (P, L, H)-Y*Den _line (P, L, H) = 0;

F _X = Num _pixel (P, L, H)-X*Den _pixel (P, L, H) = 0;

From the above formula, we know that row and column RPCs are independent of each other and can be solved separately; taking row RPCs as an example, the first formula in the above formula can be written as:

F _Y ＝a ₀ +a ₁ L+a ₂ P+a ₃ H+a ₄ LP+a ₅ LH+a ₆ PH+a ₇ L ² +a ₈ P ² +a ₉ H ² +a ₁₀ PLH+a ₁₁ L ³ +a ₁₂ LP ² +a ₁₃ LH ² +a ₁₄ L ² P+a ₁₅ P ³ +a ₁₆ PH ² +a ₁₇ L ² H+a ₁₈ P ² H+a ₁₉ H ³ -Y*( b ₀ +b ₁ L+b ₂ P+b ₃ H+b ₄ LP+b ₅ LH+b ₆ PH+b ₇ L ² +b ₈ P ² +b ₉ H ² +b ₁₀ PLH+b ₁₁ L ³ +b ₁₂ LP ² +b ₁₃ LH ² +b ₁₄ L ² P+b ₁₅ P ³ +b ₁₆ PH ² +b ₁₇ L ² H+b ₁₈ P ² H+b ₁₉ H ³ );

(2) Simultaneously combine all observation data to establish the observation equation error equation group: F _Y = BV _Y + ε; where:

V _Y = [a ₀ a ₁ ... a ₁₉ b ₁ b ₂ ... b ₁₉ ] ^T ;

Among them, ε is a random error. According to the principle of least squares adjustment, the error equation can be transformed into a normal equation for solution: (B ^T B)V _Y ＝B ^T F _Y .

As a further improvement of the technical solution of the present invention, for the ill-conditioned problem of the equation coefficient matrix described in step S7, the unbiased estimation of the RPCs parameters is obtained by using the least squares solution based on the spectral correction iterative method;

(1) For the ill-conditioned problem of the coefficient matrix of the equation, the spectral correction iteration method is used to solve the normal equation, and the unbiased parameter estimation that meets the requirements can be obtained when the normal equation is well-conditioned or ill-conditioned. The normal equation is corrected by the spectrum correction iterative method as follows:

(B ^T B+E)V _Y =B ^T F _Y +V _Y ;

Among them, E is the identity matrix, and after k iterations, the parameters can be estimated as follows:

V _Y ^(k) = (B ^T B + E) ^-1 (B ^T F _Y + V _Y ^(k-1) );

Similarly, the kth parameter estimation in the column direction is obtained as follows:

V _X ^(k) = (A ^T A + E) - ¹ (A ^T F _X +V _X ^(k-1) );

in,

V _X = [c ₀ c ₁ ... c ₁₉ d ₁ d ₂ ... d ₁₉ ] ^T ;

Among them, F _X ＝c ₀ +c ₁ L+c ₂ P+c ₃ H+c ₄ LP+c ₅ LH+c ₆ PH+c ₇ L ² +c ₈ P ² +c ₉ H ² +c ₁₀ PLH +c ₁₁ L ³ +c ₁₂ LP ² +c ₁₃ LH ² +c ₁₄ L ² P+c ₁₅ P ³ +c ₁₆ PH ² +c ₁₇ L ² H+c ₁₈ P ² H+c ₁₉ H ³ -X *(d ₀ +d ₁ L+d ₂ P+d ₃ H+d ₄ LP+d ₅ LH+d ₆ PH+d ₇ L ² +d ₈ P ² +d ₉ H ² +d ₁₀ PLH+d ₁₁ L ³ +d ₁₂ LP ² +d ₁₃ LH ² +d ₁₄ L ² P+d ₁₅ P ³ +d ₁₆ PH ² +d ₁₇ L ² H+d ₁₈ P ² H+d ₁₉ H ³ );

(2) During the solution process, if the difference between the two result values before and after the iterative calculation is less than the specified limit, the calculation stops. Generally, the difference between the two estimates is less than 10 ^-5 to meet the accuracy requirements, and the model parameters RPCs are obtained to generate rational functions Calibration model.

The advantage of the present invention is that: the technical solution of the present invention uses the inverse RD positioning model and real DEM information to uniformly sample the real space of the imaging area, and generates an array of input control points with reasonable distribution for RPCs to accurately solve. The use of redundant information in the input array of the RPC model, on the other hand, can fit the real imaging geometric model of the imaging area to be corrected, thereby improving the fitting accuracy of the RPC model and ensuring the RPC geometric positioning accuracy of the SAR image.

Description of drawings

FIG. 1 is a flow chart of a method for acquiring RPC correction parameters of a spaceborne SAR image based on an inverse RD positioning model according to an embodiment of the present invention.

Fig. 2 is a schematic diagram of geocentric coordinates used for RD positioning model description of a method for acquiring RPC correction parameters of a spaceborne SAR image based on an inverse RD positioning model according to an embodiment of the present invention.

detailed description

In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the embodiments of the present invention. Obviously, the described embodiments are part of the present invention Examples, not all examples. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without creative efforts fall within the protection scope of the present invention.

The technical solution of the present invention will be further described below in conjunction with the accompanying drawings of the description and specific embodiments:

Embodiment one

As shown in Figure 1, the method for obtaining RPC correction parameters of spaceborne SAR images based on the inverse RD positioning model includes the following steps:

1. SAR imaging parameters and satellite orbit data acquisition

Obtain spaceborne SAR system imaging parameters, image parameters, satellite orbit data, and ground elevation data, etc.; spaceborne SAR satellite orbit data includes position and velocity vector data at sampling points; SAR system imaging and image parameters include radar wavelength and image size , imaging Doppler frequency information, imaging time T ₀ of the first line of the image, pulse delay t ₀ of the first pixel of the image or image short distance R ₀ , pulse repetition frequency PRF, range sampling frequency RSR, etc. Using the above parameters, the imaging time t=T ₀ +m/PRF+t ₀ /2+n/RSR corresponding to any point (m, n) on the image can be obtained.

2. Polynomial fitting of satellite orbit data

According to the satellite position and velocity data at multiple times in the SAR imaging system and image parameters, the function curve of the satellite orbit with respect to the azimuth pixel is fitted; for the solution of the SAR satellite position, usually the width of the entire SAR image in the azimuth direction is 10,000 lines The order of magnitude, each row corresponds to a calculation of satellite position and velocity. To this end, the satellite ephemeris can be expressed as a row parameter polynomial. Perform polynomial fitting on satellite orbit data with pixels as independent variables:

The polynomial fit model with respect to location is as follows:

Among them, m represents a pixel, x _s , y _s , and z _s represent the position of the satellite orbit corresponding to the pixel in the x direction, y direction, and z direction, and α _i , β _i , and γ _i are respectively the orbit curve polynomials to be fitted coefficient;

The polynomial fitting model for velocity is as follows:

Among them, vx _s , vy _s , vz _s represent the velocity of the satellite orbit corresponding to pixel m at positions x _s , y _s , z _s .

In actual processing, the fitting orders of the above polynomials can be properly chosen. In addition, in order to improve the calculation accuracy, the satellite position and velocity can also be fitted separately, and then the fitted polynomial coefficients can be averaged to obtain the polynomial coefficients after comprehensive processing.

3. Positioning of the four corners of the SAR image

Calculate the ground range corresponding to the four corner points of the image from the rigorous imaging geometric model of the satellite image. Using the satellite orbit positions corresponding to the pixels obtained in the previous step, the three positioning equations of the slant range-Doppler-Earth ellipsoid model are constructed, and the ground coordinates of the four corners of the SAR image are iteratively solved, that is, the first The WGS84 coordinates of the nearest slant distance and the farthest slant distance corresponding to the first azimuth time and the last azimuth time respectively determine the range of the ground area of the SAR image; at present, the system-level geometric correction method generally uses iterative solutions to solve the equations consisting of the slope range equation, Doppler Equations and nonlinear equations constructed by the earth model to realize spaceborne SAR positioning. This method first constructs the nonlinear equations of the slope-Doppler-Earth model based on the satellite orbit data, SAR system parameters and imaging processing results; then solves the pixel geographic location according to the average elevation, according to the SRTM DEM Extract the elevation value at the corresponding position, and then re-substitute the new elevation value into the above nonlinear equations to calculate the new pixel position. After several iterations, the accurate pixel position is obtained.

The satellite orbit data provided by commonly used spaceborne SAR data are defined in the ECR coordinate system, so the RD positioning model is also established on the ECR coordinate system for discussion here. As shown in Figure 2, the Z-axis of the ECR coordinate system coincides with the earth's rotation axis, the X-axis points from the center of the earth to the intersection of the Greenwich meridian and the equatorial plane, and the Y-axis passes through the center of the earth and is perpendicular to the ZXO plane, which is determined by the right-handed spiral criterion direction. If the spaceborne SAR is at position s, and the ground imaging target point is at position t, then in the ECR coordinate system O-XYZ, the RD positioning model consists of the following three equations:

其中，R_e、R_p分别为椭球体的半长轴、半短轴；h为目标点的高程，通常可设为成像区域的平均高程；

为目标点t的位置矢量坐标。1) Earth ellipsoid model equation:

Among them, R _e and R _p are the semi-major axis and semi-minor axis of the ellipsoid respectively; h is the elevation of the target point, which can usually be set as the average elevation of the imaging area;

is the position vector coordinate of the target point t.

其中，

是卫星位置矢量，

R为卫星到目标的斜距，即

2) Slope distance equation:

in,

is the satellite position vector,

R is the slant distance from the satellite to the target, that is

其中，λ为雷达波长，

卫星速度矢量，

为目标速度矢量，采用ECR坐标系统定义时其值可简化为0。3) Doppler frequency equation:

where λ is the radar wavelength,

satellite velocity vector,

is the target velocity vector, and its value can be simplified to 0 when defined by the ECR coordinate system.

4. Ground 3D grid establishment

According to the geographical coordinates of the four corner points, the required geographic range of the image is obtained, and then the elevation data information within the range is obtained based on the global ground DEM data; it is uniformly divided into grids such as M*N, and the highest elevation in each grid is calculated. Three-dimensional coordinate information at the lowest and intermediate values (lat _{i, j} , lon _{i, j} , hei _{i, j} , i=0, 1, 2,..., 3M; j=0, 1, 2,..., 3N) , that is, a total of 3M*N three-dimensional data of ground control points composed of three dimensions of longitude, latitude and height are obtained; at the same time, for the control point data set, its WGS84 coordinates are converted into Cartesian coordinates (x _t (i, j), y _t (i, j), z _t (i, j), i=0, 1, 2, ..., 3M; j = 0, 1, 2, ..., 3N), get the pixel in the New position vector control point data set under the central coordinate system;

为地球第一偏心率。in,

is the first eccentricity of the earth.

5. Inverse RD image positioning of ground control point array

Construct the inverse slant range-Doppler nonlinear equation, according to the satellite orbit pixel fitting curve, solve the inverse RD equation Newton iteratively for each control point, and obtain the image pixel position (m(i, j) corresponding to each control point , n(i,j), i=0,1,2,...,3M; j=0,1,2,...,3N);

According to the satellite fitting curve, the slant range error equation and the Doppler frequency error equation are constructed as follows:

Among them, f _s is the sampling frequency RSR, C is the speed of light, λ is the radar wavelength, and f _d (n) is the Doppler center frequency corresponding to the pixel. Combined into f _d (n) = f ₀ +f ₁ n+f ₂ n ² ;

The solution of the inverse slope range-Doppler equation about the pixel position (m, n) adopts the Newton iterative method as follows:

Among them, (m ₀ , n ₀ ) is the preset initial value, then the kth group of solutions in the iterative process can be obtained as follows:

in,

可简单通过轨道拟合公式获得。通过多次迭代循环，当相邻两个解的差值精度达到收敛值，即可停止迭代求解，获取目标像素位置。Among them, the partial conduction coefficient

It can be obtained simply by orbital fitting formula. Through multiple iterative cycles, when the difference accuracy of two adjacent solutions reaches the convergence value, the iterative solution can be stopped and the target pixel position can be obtained.

6. Compliance prediction and regular processing of control point arrays

For the obtained image pixel position (m(i, j), n(i, j)) and its corresponding ground three-dimensional coordinate information (x _t (i, j), y _t (i, j), z _t (i , j)) Pre-judgment of compliance, whether it is within the image range, remove points that are not within the image range, and perform data regularization processing, normalize the ground coordinates and image coordinates, in order to reduce the order of magnitude difference in the calculation process The rounding error introduced by too large;

Normalize the image coordinates and ground coordinates of the 3M*N array of control points. Let (P, L, H) be the ground coordinates after the standardization of the ground three-dimensional coordinates, and (X, Y) be the image coordinates after the normalization of the row and column values of the image pixel positions, and the normalization formula is as follows:

和

为地面坐标的标准化参数，也即偏移值和比例值；m_off，m_scale，n_off和n_scale为像素坐标的标准化参数。in,

with

is the normalization parameter of ground coordinates, that is, offset value and scale value; m _off , m _scale , n _off and n _scale are normalization parameters of pixel coordinates.

The calculation method of offset value and scale value is as follows:

m _scale = max(|m(i, j) _max -m _off |, |m(i, j) _min -m _off |);

n _scale = max(|n(i, j) _max -n _off |, |n(i, j) _min -n _off |).

7. RPC model equation construction and least square solution

Establish a ratio polynomial between ground coordinates and image coordinates, and associate the ground point coordinates (P, L, H) with their corresponding pixel coordinates (X, Y) with a ratio polynomial:

Among them, Num _line (P, L, H) = a ₀ +a ₁ L+a ₂ P+a ₃ H+a ₄ LP+a ₅ LH+a ₆ PH+a ₇ L ² +a ₈ P ² +a ₉ H ² +a ₁₀ PLH+a ₁₁ L ³ +a ₁₂ LP ² +a ₁₃ LH ² +a ₁₄ L ² P+a ₁₅ P ³ +a ₁₆ PH ² +a ₁₇ L ² H+a ₁₈ P ² H +a ₁₉ H ³ ;

Den _line (P, L, H)＝b ₀ +b ₁ L+b ₂ P+b ₃ H+b ₄ LP+b ₅ LH+b ₆ PH+b ₇ L ² +b ₈ P ² +b ₉ H ² +b ₁₀ PLH+b ₁₁ L ³ +b ₁₂ LP ² +b ₁₃ LH ² +b ₁₄ L ² P+b ₁₅ P ³ +b ₁₆ PH ² +b ₁₇ L ² H+b ₁₈ P ² H+b ₁₉ H ³ ;

Num _pixel (P, L, H)＝c ₀ +c ₁ L+c ₂ P+c ₃ H+c ₄ LP+c ₅ LH+c ₆ PH+c ₇ L ² +c ₈ P ² +c ₉ H ² +c ₁₀ PLH+c ₁₁ L ³ +c ₁₂ LP ² +c ₁₃ LH ² +c ₁₄ L ² P+c ₁₅ P ³ +c ₁₆ PH ² +c ₁₇ L ² H+c ₁₈ P ² H+c ₁₉ H ³ ;

Den _pixel (P, L, H)＝d ₀ +d ₁ L+d ₂ P+d ₃ H+d ₄ LP+d ₅ LH+d ₆ PH+d ₇ L ² +d ₈ P ² +d ₉ H ² +d ₁₀ PLH+d ₁₁ L ³ +d ₁₂ LP ² +d ₁₃ LH ² +d ₁₄ L ² P+d ₁₅ P ³ +d ₁₆ PH ² +d ₁₇ L ² H+d ₁₈ P ² H+d ₁₉ H ³ ;

(a _i , b _i , _ci , d _i , i=0, 1, ..., 19) in the above formula are the coefficients RPCs of the rational function model, and the comparison value polynomial is transformed to obtain the rational function model for each RPCs Linear model, and express the error equation listed by the control points in matrix form;

F _Y =Num _line (P, L, H)-Y*Den _line (P, L, H) = 0;

F _X = Num _pixel (P, L, H)-X*Den _pixel (P, L, H) = 0;

From the above formula, we know that row and column RPCs are independent of each other and can be solved separately; taking row RPCs as an example, the first formula in the above formula can be written as:

F _Y ＝a ₀ +a ₁ L+a ₂ P+a ₃ H+a ₄ LP+a ₅ LH+a ₆ PH+a ₇ L ² +a ₈ P ² +a ₉ H ² +a ₁₀ PLH+a ₁₁ L ³ +a ₁₂ LP ² +a ₁₃ LH ² +a ₁₄ L ² P+a ₁₅ P ³ +a ₁₆ PH ² +a ₁₇ L ² H+a ₁₈ P ² H+a ₁₉ H ³ -Y*( b ₀ +b ₁ L+b ₂ P+b ₃ H+b ₄ LP+b ₅ LH+b ₆ PH+b ₇ L ² +b ₈ P ² +b ₉ H ² +b ₁₀ PLH+b ₁₁ L ³ +b ₁₂ LP ² +b ₁₃ LH ² +b ₁₄ L ² P+b ₁₅ P ³ +b ₁₆ PH ² +b ₁₇ L ² H+b ₁₈ P ² H+b ₁₉ H ³ );

Simultaneously combine all the observation data to establish the observation equation error equation group: F _Y ＝ BV _Y + ε; where:

V _Y = [a ₀ a ₁ ... a ₁₉ b ₁ b ₂ ... b ₁₉ ] ^T ;

Among them, ε is a random error. According to the principle of least squares adjustment, the error equation can be transformed into a normal equation for solution: (B ^T B)V _Y ＝B ^T F _Y .

In the solution of RPC model, due to the singularity or near-singularity of the normal equation, the least squares adjustment often fails to converge. In order to overcome the ill-conditioned problem of the coefficient matrix, the spectral correction iterative method is used to solve the coefficient estimation equation, and the unbiased parameter estimation that meets the requirements can be obtained when the normal equation is well-conditioned or ill-conditioned. The normal equation is corrected by the spectrum correction iterative method as follows:

(B ^T B+E)V _Y =B ^T F _Y +V _Y ;

Among them, E is the identity matrix, and after k iterations, the parameters can be estimated as follows:

V _Y ^(k) = (B ^T B + E) ^-1 (B ^T F _Y + V _Y ^(k-1) );

Similarly, the kth parameter estimation in the column direction is obtained as follows:

V _X ^(k) = (A ^T A + E) ^-1 (A ^T F _X + V _X ^(k-1) );

in,

V _X = [c ₀ c ₁ ... c ₁₉ d ₁ d ₂ ... d ₁₉ ] ^T ;

Among them, F _X ＝c ₀ +c ₁ L+c ₂ P+c ₃ H+c ₄ LP+c ₅ LH+c ₆ PH+c ₇ L ² +c ₈ P ² +c ₉ H ² +c ₁₀ PLH +c ₁₁ L ³ +c ₁₂ LP ² +c ₁₃ LH ² +c ₁₄ L ² P+c ₁₅ P ³ +c ₁₆ PH ² +c ₁₇ L ² H+c ₁₈ P ² H+c ₁₉ H ³ -X *(d ₀ +d ₁ L+d ₂ P+d ₃ H+d ₄ LP+d ₅ LH+d ₆ PH+d ₇ L ² +d ₈ P ² +d ₉ H ² +d ₁₀ PLH+d ₁₁ L ³ +d ₁₂ LP ² +d ₁₃ LH ² +d ₁₄ L ² P+d ₁₅ P ³ +d ₁₆ PH ² +d ₁₇ L ² H+d ₁₈ P ² H+d ₁₉ H ³ );

During the solution process, if the difference between the two result values before and after the iterative calculation is less than the specified limit, the calculation stops. Generally, the difference between the two estimates is less than 10 ^-5 to meet the accuracy requirements, and the model parameters RPCs are obtained to generate a rational function correction model.

8. Model accuracy verification

The image coordinates corresponding to the checkpoints are calculated by using the solved RPC model parameters, and the image square accuracy of the RPC model parameters is evaluated by comparing the calculation results with the rigorous imaging geometric model, and the acquisition and accuracy analysis of the spaceborne SAR geometric correction RPCs parameters are completed.

The present invention provides an accurate acquisition method of spaceborne SAR image RPC correction parameters based on the inverse RD positioning model. The method provides a complete and robust RPCs parameter calculation process for spaceborne SAR geometric correction. The main optimization steps include: 1) according to The SAR rigorous imaging model uses the inverse RD positioning algorithm to generate an array of control points with reasonable distribution, which is used for RPC equation model input and accurate model parameter generation; Solve the ill-conditioned problem of the equation caused by too many unknowns and strong correlation, so that the precise geometric correction of SAR images can be guaranteed and improved. The method of the present invention borrows the DEM data of the real imaging area, realizes uniform sampling of the real three-dimensional space of the imaging area based on the inverse RD positioning model, and generates an array of ground control points with reasonable distribution, instead of conventional virtual mapping of multiple elevation profiles, and then Redundant information input is reduced, real terrain imaging projection information is increased, and the accuracy of calibration model fitting is improved. In addition, the method of the present invention is from the input of the control data group generated by the inverse RD transformation to the robust solution of the RPC equation and the output of the correction model coefficients.

The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: it can still be described in the foregoing embodiments Modifications are made to the recorded technical solutions, or equivalent replacements are made to some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims (10)

1. a method for obtaining the RPC correction parameters of spaceborne SAR images based on the inverse RD positioning model, is characterized in that, comprising the following steps: S1. Acquisition of SAR imaging parameters and satellite orbit data: acquisition of spaceborne SAR system imaging and image parameters, satellite orbit data and ground elevation data; S2. Polynomial fitting of satellite orbit data: according to the satellite position and velocity data at multiple moments in the SAR imaging system and image parameters, fit the function curve of the satellite orbit with respect to the azimuth pixel; S3. Positioning of the four corners of the SAR image: Construct the positioning equation of the slant range-Doppler-Earth ellipsoid model, and iteratively solve the ground coordinates of the four corners of the SAR image, which are the first azimuth moment and the corresponding time of the SAR image respectively. The WGS84 coordinates of the closest slant distance and the farthest slant distance corresponding to the last azimuth moment respectively, determine the range of the ground area of the SAR image; S4. Ground three-dimensional grid establishment: obtain the elevation data information within the geographic range of the required image according to the geographic coordinates of the four corner points and the global ground DEM data; divide it evenly into M*N grids, and calculate each grid The three-dimensional coordinate information of the highest, lowest, and intermediate values of the elevation in the field can obtain 3M*N ground control point data in total; and then carry out Cartesian coordinate conversion to obtain a new control point position vector group under the geocentric coordinate system; S5. Inverse RD image positioning of the ground control point array: construct the slant range error equation and the Doppler frequency error equation, according to the satellite orbit pixel fitting curve, solve the inverse RD positioning equation Newton iteratively for each control point, and obtain each The image pixel position corresponding to the control point; S6. Compliance prediction and regularization processing of the control point array: perform compliance prediction on the obtained image pixel position and its corresponding ground three-dimensional coordinate information, remove points that are not within the image range; perform data regularization processing , to normalize the ground coordinates and image coordinates to reduce rounding errors in the calculation process; S7. RPC model equation construction and least squares solution: establish a ratio polynomial between ground coordinates and image coordinates, and deform it to obtain a nonlinear model about RPCs; for the ill-conditioned problem of the equation coefficient matrix, use the least squares based on the spectral correction iterative method Multiply the solution to obtain an unbiased estimate of the RPCs parameters; S8. Model accuracy verification: use the solved RPCs model parameters to calculate the image coordinates corresponding to the checkpoints, and evaluate the image square accuracy of the RPC model parameters by making a difference with the calculation results of the rigorous imaging geometric model. 2. A method for obtaining RPC correction parameters of spaceborne SAR images based on the inverse RD positioning model according to claim 1, wherein the spaceborne SAR system imaging and image parameters described in step S1 include: radar wavelength, Image size, imaging Doppler frequency information, imaging time T ₀ of the first line of the image, pulse delay t ₀ of the first pixel of the image or image short distance R ₀ , pulse repetition frequency PRF, range sampling frequency RSR; the satellite The track data includes: the position and velocity vector data of the sampling point at the moment. 3. A kind of method for obtaining RPC correction parameters of spaceborne SAR images based on the inverse RD positioning model according to claim 2, characterized in that, the method of fitting the function curve of the satellite orbit about the azimuth pixel described in step S2 as follows: S21, the polynomial fitting model about the position is as follows:

Among them, m represents a pixel, x _s , y _s , and z _s represent the position of the satellite orbit corresponding to the pixel in the x direction, y direction, and z direction, and α _i , β _i , and γ _i are the satellite orbits to be fitted about The coefficients of the function curve polynomial of the pixels in the azimuth direction, i=1, 2, 3, 4, 5; S22, the polynomial fitting model about speed is as follows:

Among them, vx _s , vy _s , vz _s represent the velocity of the satellite orbit corresponding to pixel m at positions x _s , y _s , z _s . 4. A kind of method for obtaining RPC correction parameters of spaceborne SAR images based on the inverse RD positioning model according to claim 3, characterized in that, the construction of the slant range-Doppler-Earth ellipsoid model positioning described in step S3 The equation is as follows: 1) Earth ellipsoid model equation:

Among them, R _e and R _p are the semi-major axis and semi-minor axis of the ellipsoid respectively; h is the elevation of the target point, which can usually be set as the average elevation of the imaging area;

is the position vector coordinate of the target point t; 2) Slope distance equation:

in,

is the satellite position vector,

R is the slant distance from the satellite to the target, that is

3) Doppler frequency equation:

where λ is the radar wavelength,

satellite velocity vector,

is the target velocity vector, and its value can be simplified to 0 when defined by the ECR coordinate system. 5. a kind of space-borne SAR image RPC correction parameter acquisition method based on inverse RD positioning model according to claim 4, is characterized in that, the formula that carries out Cartesian coordinate transformation described in step S4 is:

in,

is the first eccentricity of the earth. 6. a kind of space-borne SAR image RPC correction parameter acquisition method based on the inverse RD positioning model according to claim 5, is characterized in that, described in the step S5, constructs the slant range error equation and the Doppler frequency error equation, According to the satellite orbit pixel fitting curve, the inverse RD positioning equation Newton iterative solution is performed for each control point, and the method of obtaining the image pixel position corresponding to each control point is as follows: (1) According to the satellite fitting curve, construct the slant range error equation and Doppler frequency error equation as follows:

Among them, f _s is the sampling frequency RSR, C is the speed of light, λ is the radar wavelength, and f _d (n) is the Doppler center frequency corresponding to the pixel. Combined into f _d (n) = f ₀ +f ₁ n+f ₂ n ² ; (2) The Newton iterative method is used to solve the inverse slope range-Doppler equation about the pixel position (m, n) as follows:

Among them, (m ₀ , n ₀ ) is the preset initial value, then the kth group of solutions in the iterative process can be obtained as follows:

in,

Among them, the partial conduction coefficient

It can be obtained simply by orbital fitting formula; (3) Through multiple iterative loops, when the difference accuracy of two adjacent solutions reaches the convergence value, the iterative solution is stopped and the target pixel position is obtained. 7. A kind of spaceborne SAR image RPC correction parameter acquisition method based on the inverse RD positioning model according to claim 6, is characterized in that, in step S6, carry out data regularization processing, return ground coordinates and image coordinates The unified method is as follows: Standardize the image coordinates and ground coordinates of 3M*N control point arrays, let (P, L, H) be the ground coordinates after the ground three-dimensional coordinates are standardized, and (X, Y) be the images after the normalized row and column values of the image pixel positions Coordinates, the standardized formula is as follows:

in,

with

is the normalization parameter of ground coordinates, that is, offset value and scale value; m _off , m _scale , n _off and n _scale are normalization parameters of pixel coordinates. 8. a kind of space-borne SAR image RPC correction parameter acquisition method based on inverse RD positioning model according to claim 7, is characterized in that, described offset value and proportional value calculation formula are as follows:

m _scale = max(|m(i, j) _max -m _off |, |m(i, j) _min -m _off |); n _scale = max(|n(i, j) _max -n _off |, |n(i, j) _min -n _off |). 9. A kind of space-borne SAR image RPC correction parameter acquisition method based on the inverse RD positioning model according to claim 8, characterized in that, the ratio polynomial between the establishment of ground coordinates and image coordinates described in step S7, and The method of deforming to obtain a nonlinear model of RPCs is as follows: (1) Establish a ratio polynomial between ground coordinates and image coordinates, and associate the ground point coordinates (P, L, H) with their corresponding pixel point coordinates (X, Y) with a ratio polynomial:

Among them, Num _line (P, L, H) = a ₀ +a ₁ L+a ₂ P+a ₃ H+a ₄ LP+a ₅ LH+a ₆ PH+a ₇ L ² +a ₈ P ² +a ₉ H ² +a ₁₀ PLH+a ₁₁ L ³ +a ₁₂ LP ² +a ₁₃ LH ² +a ₁₄ L ² P+a ₁₅ P ³ +a ₁₆ PH ² +a ₁₇ L ² H+a ₁₈ P ² H +a ₁₉ H ³ ; Den _line (P, L, H)＝b ₀ +b ₁ L+b ₂ P+b ₃ H+b ₄ LP+b ₅ LH+b ₆ PH+b ₇ L ² +b ₈ P ² +b ₉ H ² +b ₁₀ PLH+b ₁₁ L ³ +b ₁₂ LP ² +b ₁₃ LH ² +b ₁₄ L ² P+b ₁₅ P ³ +b ₁₆ PH ² +b ₁₇ L ² H+b ₁₈ P ² H+b ₁₉ H ³ ; Num _pixel (P, L, H)＝c ₀ +c ₁ L+c ₂ P+c ₃ H+c ₄ LP+c ₅ LH+c ₆ PH+c ₇ L ² +c ₈ P ² +c ₉ H ² +c ₁₀ PLH+c ₁₁ L ³ +c ₁₂ LP ² +c ₁₃ LH ² +c ₁₄ L ² P+c ₁₅ P ³ +c ₁₆ PH ² +c ₁₇ L ² H+c ₁₈ P ² H+c ₁₉ H ³ ; Den _pixel (P, L, H)＝d ₀ +d ₁ L+d ₂ P+d ₃ H+d ₄ LP+d ₅ LH+d ₆ PH+d ₇ L ² +d ₈ P ² +d ₉ H ² +d ₁₀ PLH+d ₁₁ L ³ +d ₁₂ LP ² +d ₁₃ LH ² +d ₁₄ L ² P+d ₁₅ P ³ +d ₁₆ PH ² +d ₁₇ L ² H+d ₁₈ P ² H+d ₁₉ H ³ ; (a _i , b _i , _ci , d _i , i=0, 1, ..., 19) in the above formula are the coefficients RPCs of the rational function model, and the comparison value polynomial is transformed to obtain the rational function model for each RPCs Linear model, and express the error equation listed by the control points in matrix form; F _Y =Num _line (P, L, H)-Y*Den _line (P, L, H) = 0; F _X = Num _pixel (P, L, H)-X*Den _pixel (P, L, H) = 0; From the above formula, we know that row and column RPCs are independent of each other and can be solved separately; taking row RPCs as an example, the first formula in the above formula can be written as: F _Y ＝a ₀ +a ₁ L+a ₂ P+a ₃ H+a ₄ LP+a ₅ LH+a ₆ PH+a ₇ L ² +a ₈ P ² +a ₉ H ² +a ₁₀ PLH+a ₁₁ L ³ +a ₁₂ LP ² +a ₁₃ LH ² +a ₁₄ L ² P+a ₁₅ P ³ +a ₁₆ PH ² +a ₁₇ L ² H+a ₁₈ P ² H+a ₁₉ H ³ -Y*( b ₀ +b ₁ L+b ₂ P+b ₃ H+b ₄ LP+b ₅ LH+b ₆ PH+b ₇ L ² +b ₈ P ² +b ₉ H ² +b ₁₀ PLH+b ₁₁ L ³ +b ₁₂ LP ² +b ₁₃ LH ² +b ₁₄ L ² P+b ₁₅ P ³ +b ₁₆ PH ² +b ₁₇ L ² H+b ₁₈ P ² H+b ₁₉ H ³ ); (2) Simultaneously combine all observation data to establish the observation equation error equation group: F _Y = BV _Y + ε; where:

V _Y = [a ₀ a ₁ ... a ₁₉ b ₁ b ₂ ... b ₁₉ ] ^T ;

Among them, ε is a random error. According to the principle of least squares adjustment, the error equation can be transformed into a normal equation for solution: (B ^T B)V _Y ＝B ^T F _Y . 10. A method for obtaining RPC correction parameters of spaceborne SAR images based on an inverse RD positioning model according to claim 9, characterized in that, for the ill-conditioned problem of the equation coefficient matrix described in step S7, an iterative method based on spectral correction is adopted The least squares solution of , to obtain the unbiased estimation of RPCs parameters; (1) For the ill-conditioned problem of the coefficient matrix of the equation, the normal equation is solved by using the spectral correction iterative method, and unbiased parameter estimates that meet the requirements can be obtained when the normal equation is well-conditioned or ill-conditioned, and the normal equation is corrected by the spectral correction iterative method as follows : (B ^T B+E)V _Y =B ^T F _Y +V _Y ; Among them, E is the identity matrix, and after k iterations, the parameters can be estimated as follows: V _Y ^(k) = (B ^T B + E) ^-1 (B ^T F _Y + V _Y ^(k-1) ); Similarly, the kth parameter estimation in the column direction is obtained as follows: V _X ^(k) = (A ^T A + E) ^-1 (A ^T F _X +V _X ^(k-1) ); in,

V _X = [c ₀ c ₁ ... c ₁₉ d ₁ d ₂ ... d ₁₉ ] ^T ;

Among them, F _X ＝c ₀ +c ₁ L+c ₂ P+c ₃ H+c ₄ LP+c ₅ LH+c ₆ PH+c ₇ L ² +c ₈ P ² +c ₉ H ² +c ₁₀ PLH +c ₁₁ L ³ +c ₁₂ LP ² +c ₁₃ LH ² +c ₁₄ L ² P+c ₁₅ P ³ +c ₁₆ PH ² +c ₁₇ L ² H+c ₁₈ P ² H+c ₁₉ H ³ -X *(d ₀ +d ₁ L+d ₂ P+d ₃ H+d ₄ LP+d ₅ LH+d ₆ PH+d ₇ L ² +d ₈ P ² +d ₉ H ² +d ₁₀ PLH+d ₁₁ L ³ +d ₁₂ LP ² +d ₁₃ LH ² +d ₁₄ L ² P+d ₁₅ P ³ +d ₁₆ PH ² +d ₁₇ L ² H+d ₁₈ P ² H+d ₁₉ H ³ ); (2) During the solution process, if the difference between the two result values before and after the iterative calculation is less than the specified limit, the calculation stops. Generally, the difference between the two estimates is less than 10 ^-5 to meet the accuracy requirements, and the model parameters RPCs are obtained to generate rational functions Calibration model.

Images