Ka-sar satellite standard product generation method and computer readable medium

Abstract

The application provides a Ka-SAR satellite standard product generation method and a computer readable medium. Firstly, an intensity map converted image is generated; the intensity map image is divided into sub-blocks, and local geometric correction is performed on the intensity map image; resampling is performed on the DOM according to the range of each sub-block; a heterogeneous matching algorithm is used to obtain matching points between the DOM resampling sub-blocks and the geometric correction sub-blocks, and the matching points are respectively mapped into control point geographic coordinates and control point image coordinates on the intensity map image; then, regional network adjustment is used to remove gross errors of the control points; the RFM model is updated using the control points after the gross error removal, so that an updated RFM model is obtained; the foregoing steps are iterated to obtain a final RFM model; the intensity map image is indirectly corrected using the final RFM model, so that a Ka-SAR satellite image standard product is obtained. The application has the advantages that the positioning accuracy of the Ka-SAR image is improved through matching control points with the DOM, and support is provided for fine application of the Ka-SAR image in water conservancy and other industries.

Description

Technical Field

[0001] This invention belongs to the field of microwave remote sensing technology, and particularly relates to a method for generating Ka-SAR satellite standard products and a computer-readable medium. Background Technology

[0002] With the rise of large-scale satellite applications, micro-SAR satellites have developed rapidly, providing a new way to obtain repeated observations in a short period of time. Ka-SAR satellites can be applied to marine environmental monitoring and flood detection, showing broad application prospects. However, current Ka-SAR satellites are mainly microsatellites, which are limited by their size. Compared with standard satellites, microsatellites have relatively poor geometric performance, limiting their application in photogrammetry. Traditional methods can improve the geometric positioning accuracy of satellite images by adding well-distributed ground control points, but this requires a lot of financial and human resources. To address this, a method for generating standard Ka-SAR satellite products is proposed. Combining the high geometric positioning accuracy of standard SAR satellites (such as Gaofen-3), high-precision control points are matched with the DOM generated by standard SAR satellites to improve the positioning accuracy of Ka-SAR images, thereby generating high-precision Ka-SAR satellite standard products, providing support for their refined applications in industries such as water conservancy. The method presented in this paper can provide a new approach for collaborative processing of multi-source SAR satellite image data. Summary of the Invention

[0003] To address the problem of poor initial positioning accuracy of multi-source SAR satellites, this invention proposes a method for generating standard Ka-SAR satellite products and a computer-readable medium.

[0004] The technical solution of this invention is a method for generating Ka-SAR satellite standard products, as detailed below:

[0005] Perform intensity map conversion to generate the corresponding intensity map image;

[0006] The intensity map image is divided into sub-blocks;

[0007] Calculate the output resolution for geometric correction based on the intensity map image information;

[0008] Local geometric correction is performed on the intensity map image to obtain multiple geometrically corrected sub-blocks of the intensity map image. Based on the range of each sub-block, the pixel values ​​of the corresponding region are resampled on the DOM to obtain multiple resampled sub-blocks of the DOM.

[0009] A heterogeneous matching algorithm is used to obtain the geographic coordinates of multiple pairs of points located in the DOM resampled sub-block and the geographic coordinates of the intensity map image geometrically corrected sub-block; the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block and the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image are then mapped sequentially.

[0010] Gross errors were eliminated by using regional network adjustment to obtain accurate 3D geographic coordinates of control points on the DOM and pixel coordinates of control points in intensity map imagery.

[0011] The RFM model is updated by using the precise geographic coordinates of control points on the DOM and the pixel coordinates of control points in the intensity map image.

[0012] The final RFM model is obtained through iterative updates. Using the final RFM model and auxiliary DEM data, the intensity map image is indirectly corrected to obtain the standard product of Ka-SAR satellite imagery.

[0013] The specific steps of this invention are as follows:

[0014] Step 1: Acquire the original Ka-SAR single-view complex image, auxiliary rational function model parameter file, C-band Gaofen-3 DOM, and digital elevation model;

[0015] Step 2: Perform intensity map transformation on the original Ka-SAR single-look complex image to generate the corresponding intensity map image;

[0016] Step 3: Obtain the pixel width and pixel height of the intensity map image, and divide the intensity map image into sub-blocks based on the pixel width and pixel height;

[0017] Step 4: Calculate the output resolution for geometric correction based on the intensity map image information;

[0018] Step 5: Based on each geographic sub-block, combined with DEM auxiliary data, perform local geometric correction on the intensity map image to obtain multiple geometrically corrected sub-blocks of the intensity map image. Based on the range of each sub-block, resample the pixel values ​​of the corresponding area on the DOM to obtain multiple resampled sub-blocks of the DOM.

[0019] Step 6: Using each DOM resampled sub-block as the reference image and the geometrically corrected sub-block of the intensity map image as the registration image, a heterogeneous matching algorithm is used to obtain multiple sets of corresponding point pairs located at the geographic coordinates of the DOM resampled sub-block and the geographic coordinates of the geometrically corrected sub-block of the intensity map image; the geographic coordinates of the DOM resampled sub-block are combined with the auxiliary DEM and mapped to the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM; the geographic coordinates of each geometrically corrected sub-block of the intensity map image are mapped to the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image.

[0020] Step 7: Use regional network adjustment to remove outliers by matching the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM with the pixel coordinates of the geographic coordinates of each geometrically corrected sub-block on the intensity map image, so as to obtain the accurate geographic 3D coordinates of the control points on the DOM and the pixel coordinates of the control points on the intensity map image.

[0021] Step 8: Update the RFM model using the accurate geographic coordinates of control points on the DOM and the pixel coordinates of control points in the intensity map image to obtain the updated RFM model;

[0022] Step 9: Replace the RFM model of the intensity map image with the accurate RFM model updated in Step 8, and repeat Steps 3-8 to obtain the final RFM model.

[0023] Step 10: Using the final RFM model and auxiliary DEM data, the intensity map image is indirectly corrected to obtain the standard product of Ka-SAR satellite imagery.

[0024] Preferably, step 2, which involves generating the corresponding intensity map image, is as follows:

[0025] The Ka-SAR single-look complex image is defined as

[0026] The intensity map image is defined as I intensity ;

[0027] The pixel values ​​corresponding to the intensity map in row p and column q are calculated as follows:

[0028]

[0029] in,

[0030] This represents the pixel value corresponding to the p-th row and q-th column of the first band in a Ka-SAR single-look complex image.

[0031] This represents the pixel value corresponding to the p-th row and q-th column of the second band in a Ka-SAR single-view complex image.

[0032] I intensity (p,q) represents the pixel value corresponding to the intensity map image in row p and column q;

[0033] As a preferred embodiment, step 3 involves dividing the intensity map image into sub-blocks, and the specific process is as follows:

[0034] Set the number of molecular blocks in the column direction to nGX, and set the number of molecular blocks in the row direction to nGY;

[0035] Calculate the column spacing and row spacing of each sub-block separately, as follows:

[0036]

[0037] Where CenX represents the spacing of each sub-block in the column direction, CenY represents the spacing of each sub-block in the row direction, Width represents the pixel width of the intensity map image, Height represents the pixel height of the intensity map image, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction.

[0038] The sub-block in row j and column i is defined as:

[0039] CEL i,j

[0040] i∈[1,nGX], j∈[1,nGY]

[0041] Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction;

[0042] Set the sub-block size to Kernel;

[0043] Calculate the top-left pixel coordinates of the sub-block in row j and column i, as follows:

[0044]

[0045]

[0046]

[0047] Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the top-left pixel coordinates of the sub-block in row j and column i. The column coordinate represents the top-left pixel coordinate of the sub-block in row j and column i. The row coordinate represents the top-left pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block.

[0048] The geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in row j and column i are defined as follows:

[0049]

[0050] in, This represents the longitude of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the latitude of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the elevation of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column;

[0051] Calculate the top-right pixel coordinates of the sub-block in row j and column i, as follows:

[0052]

[0053]

[0054]

[0055] Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. Represents the top-right pixel coordinates of the sub-block in row j and column i; The column coordinate represents the top-right pixel coordinate of the sub-block in row j and column i; The row coordinate represents the top-right pixel coordinate of the sub-block in row j and column i; Kernel represents the size of the sub-block.

[0056] Combining DEM-aided data and the RFM model, the 3D geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column are calculated using ray tracing, defined as:

[0057]

[0058] in, This represents the longitude of the geographic coordinates corresponding to the top-right pixel coordinates of the sub-block in row j and column i. This represents the latitude of the geographic coordinates corresponding to the top right pixel coordinates of the sub-block in the j-th row and i-th column; This represents the elevation of the geographic coordinates corresponding to the top-right pixel coordinates of the sub-block in the j-th row and i-th column;

[0059] Calculate the bottom-left pixel coordinates of the sub-block in row j and column i, as follows:

[0060]

[0061]

[0062]

[0063] Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the bottom-left pixel coordinates of the sub-block in row j and column i; The column coordinate represents the bottom left pixel coordinate of the sub-block in the j-th row and i-th column; The row coordinate represents the bottom left pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block.

[0064] The geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in row j and column i are defined as follows:

[0065]

[0066] in, This represents the longitude of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the latitude of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the elevation of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column;

[0067] Calculate the bottom right pixel coordinates of the sub-block in row j and column i, as follows:

[0068]

[0069]

[0070]

[0071] Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the bottom right pixel coordinates of the sub-block in row j and column i. The column coordinate represents the bottom right pixel coordinate of the sub-block in the j-th row and i-th column; The row coordinate represents the bottom right pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block.

[0072] The geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i are defined as follows:

[0073]

[0074] in, This represents the longitude of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i. This represents the latitude of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i. This represents the elevation of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in the j-th row and i-th column;

[0075] The minimum longitude among the longitudes of the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is used as the minimum longitude of the sub-block in the j-th row and i-th column, denoted as .

[0076] The minimum longitude among the geographic coordinates corresponding to the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is taken as the maximum longitude of the sub-block in the j-th row and i-th column, denoted as .

[0077] The minimum latitude of the geographic sub-block in the j-th row and i-th column is calculated from the latitudes of its top-left, top-right, bottom-left, and bottom-right pixel coordinates. This minimum latitude is denoted as the minimum latitude of the geographic sub-block in the j-th row and i-th column.

[0078] The minimum value among the latitudes of the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is taken as the maximum latitude of the sub-block in the j-th row and i-th column, denoted as .

[0079] Will Construct the sub-block of geography in the j-th row and i-th column in a counter-clockwise order.

[0080] As a preferred embodiment, step 4, which involves calculating the geometrically corrected output resolution based on the intensity map image information, is specifically as follows:

[0081] Read the intensity map image resolution Res using the GDAL library;

[0082] Pre-set the output resolution sampling factor threshold T for geometric correction;

[0083] The output resolution of the geometric correction is FRes = T × Res;

[0084] Preferably, step 5 involves performing local geometric correction on the intensity map image to obtain multiple geometrically corrected sub-blocks of the intensity map image, as detailed below:

[0085] According to the sub-block geographic sub-block CEL in row j and column i i,j Coordinates of the four corresponding vertices For this region, local geometric correction is performed using an indirect method based on the output resolution FRES of the geometrically corrected image, combined with the intensity map image RFM model and auxiliary DEM data, to obtain the geometrically corrected sub-block of the intensity map image, defined as Ortho_CEL. i,j , i∈[1,nGX], j∈[1,nGY];

[0086] Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction;

[0087] Step 5 involves resampling the pixel values ​​of the corresponding region on the DOM to obtain multiple resampled sub-blocks of the DOM, as detailed below:

[0088] According to the sub-block geographic sub-block CEL in row j and column i i,j The coordinates of the four corresponding vertices For this region, bilinear interpolation is used to resample it in the DOM according to the output resolution FRes, resulting in a DOM resampled sub-block, defined as DOM_CEL. i,j , i∈[1,nGX], j∈[1,nGY];

[0089] Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction;

[0090] Preferably, step 6 involves combining the geographic coordinates of the DOM resampled sub-blocks with the auxiliary DEM to map the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM.

[0091] The elevation of each DOM resampled sub-block is obtained from the DEM auxiliary data using bilinear interpolation based on the corresponding longitude and latitude positions of each DOM resampled sub-block.

[0092] Combine the geographic coordinates and elevation of each DOM resampled sub-block to construct the 3D geographic coordinates of the control points corresponding to each DOM resampled sub-block on the DOM.

[0093] Step 6 involves mapping the geographic coordinates of each geometrically corrected sub-block of the intensity map image to the corresponding pixel coordinates on the intensity map image:

[0094] The elevation of each geometrically corrected sub-block of the intensity map image is obtained from the DEM auxiliary data by bilinear interpolation based on the geographic coordinates of each geometrically corrected sub-block in the intensity map image at the corresponding longitude and latitude positions.

[0095] The geographic 3D coordinates of each geometrically corrected sub-block of the intensity map image are constructed by combining the geographic coordinates of each geometrically corrected sub-block of the intensity map image with the elevation of each geometrically corrected sub-block of the intensity map image.

[0096] By combining the geographic 3D coordinates of each geometrically corrected sub-block of the intensity map image, the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image are obtained by forward calculation based on the RFM model of the intensity map image.

[0097] Preferably, step 7, which involves using regional network adjustment to remove outliers by matching the geographic 3D coordinates of the control points corresponding to each DOM resampling sub-block on the DOM with the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image, specifically involves:

[0098] Step 7.1: RFM model construction;

[0099]

[0100] in, This represents the regularized pixel coordinates of the g-th control point;

[0101] Represents the regularized row coordinates of the g-th control point;

[0102] Represents the regularized column coordinates of the g-th control point;

[0103] Represents the regularized geographic coordinates of the g-th control point;

[0104] This represents the regularized geographic coordinates (longitude) of the g-th control point;

[0105] Represents the regularized geographic coordinates (latitude) of the g-th control point;

[0106] Represents the regularized geographic coordinate elevation of the g-th control point;

[0107] P1 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1;

[0108] P2 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1;

[0109] P3 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1;

[0110] P4 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1;

[0111] It can be expressed in functional form as Let it be F r g ;

[0112] It can be expressed in functional form as Let it be F c g ;

[0113] Regularized geographic coordinates of the g-th control point Regularized pixel coordinates of the g-th control point The geographic coordinates of the g-th control point (X) g ,Y g Z g ) and the pixel coordinates of the g-th control point (R g C g The correspondence is as follows:

[0114]

[0115]

[0116] in,

[0117] This represents the regularized geographic coordinates (longitude) of the g-th control point;

[0118] Represents the regularized geographic coordinates (latitude) of the g-th control point;

[0119] Represents the regularized geographic coordinate elevation of the g-th control point;

[0120] X g This represents the longitude of the geographic coordinates of the g-th control point;

[0121] Y g This represents the latitude and geographic coordinates of the g-th control point.

[0122] Z g This represents the geographic coordinates and elevation of the g-th control point;

[0123] X o The translation parameters representing the longitude regularization of geographic coordinates are obtained from the RFM parameter file in step 1.

[0124] Y oThe translation parameter representing the latitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1;

[0125] Z o The translation parameters representing geographic coordinate elevation regularization are obtained from the RFM parameter file in step 1;

[0126] X s The scaling factor representing the longitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1;

[0127] Y s The scaling factor representing the latitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1.

[0128] Z s The scaling factor representing geographic coordinate elevation regularization is obtained from the RFM parameter file in step 1;

[0129] Represents the regularized row coordinates of the g-th control point;

[0130] Represents the regularized column coordinates of the g-th control point;

[0131] R g This represents the row coordinates of the g-th control point;

[0132] C g Represents the column coordinates of the g-th control point;

[0133] R o The translation parameters representing row coordinate regularization are obtained from the RFM parameter file in step 1;

[0134] C o The translation parameters representing column coordinate regularization are obtained from the RFM parameter file in step 1;

[0135] R s The scaling factor for row coordinate regularization is obtained from the RFM parameter file in step 1.

[0136] C s The scaling factor for column coordinate regularization is obtained from the RFM parameter file in step 1.

[0137] Step 7.2: Constructing the image-side compensation model based on the affine transformation model;

[0138] When using the image-based compensation scheme to compensate for systematic errors in RFM, the regularized geographic coordinates of the g-th control point described by RFM are... The pixel coordinates of the g-th control point (R) g C g Revised to:

[0139]

[0140] in,

[0141] R g This represents the row coordinates of the g-th control point;

[0142] C g Represents the column coordinates of the g-th control point;

[0143] This represents the regularized geographic coordinates (longitude) of the g-th control point;

[0144] Y n g Represents the regularized geographic coordinates (latitude) of the g-th control point;

[0145] Represents the regularized geographic coordinate elevation of the g-th control point;

[0146] P1 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1;

[0147] P2 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1;

[0148] P3 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1;

[0149] P4 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1;

[0150] R o The translation parameters representing row coordinate regularization are obtained from the RFM parameter file in step 1;

[0151] C o The translation parameters representing column coordinate regularization are obtained from the RFM parameter file in step 1;

[0152] R s The scaling factor for row coordinate regularization is obtained from the RFM parameter file in step 1.

[0153] C s The scaling factor for column coordinate regularization is obtained from the RFM parameter file in step 1.

[0154] (ΔR g ,ΔC g ) represents the pixel coordinates (R) of the g-th control point. g C gThe system error compensation value,

[0155] ΔR g The system error compensation value represents the row coordinate of the g-th control point;

[0156] ΔC g This represents the system error compensation value for the column coordinates of the g-th control point;

[0157] The calculation method is as follows;

[0158]

[0159] Where (e0,e1,e2,f0,f1,f2) are the affine transformation parameters of the intensity map image.

[0160] e0 is the first parameter of the row direction affine transformation, e1 is the second parameter of the row direction affine transformation, and e2 is the third parameter of the row direction affine transformation.

[0161] f0 is the first parameter of the column direction affine transformation, f1 is the second parameter of the column direction affine transformation, and f2 is the third parameter of the column direction affine transformation;

[0162] Step 7.3: Establish an error equation for each control point, and combine all error equations into a total error equation;

[0163] The error equation for the g-th control point is in the following form:

[0164]

[0165] Treating the affine transformation coefficients as unknowns, the error equation for the g-th control point is rewritten in matrix form:

[0166] V g =A g tl g

[0167] in,

[0168] This is the set of residual vectors of the row and column coordinate observations of the image of the g-th control point;

[0169] Let be the residual vector of the row coordinate observation values ​​of the image of the g-th control point;

[0170] The residual vector of the image column coordinate observation values ​​of the g-th control point;

[0171] This is the coefficient matrix corresponding to the g-th control point;

[0172] R gThis represents the row coordinates of the g-th control point;

[0173] C g Represents the column coordinates of the g-th control point;

[0174] t=[Δe0 Δe1 Δe2 Δf0 Δf1 Δf2] T This represents the incremental vector of error compensation parameters for the intensity map image coordinate system;

[0175] Δe0 is the increment of the first systematic error compensation parameter in the row direction;

[0176] Δe1 is the increment of the second systematic error compensation parameter in the row direction;

[0177] Δe2 is the increment of the third systematic error compensation parameter in the row direction;

[0178] Δf0 is the increment of the first systematic error compensation parameter in the column direction; Δf1 is the increment of the second systematic error compensation parameter in the column direction;

[0179] Δf2 is the increment of the third systematic error compensation parameter in the column direction;

[0180] in,

[0181] The difference between the row coordinates of the g-th control point and the geographic coordinates of the g-th control point, calculated by substituting the approximate values ​​of the affine transformation parameters into the RFM model.

[0182] The difference between the column coordinates of the g-th control point and the geographic coordinates of the g-th control point, calculated by substituting the approximate values ​​of the affine transformation parameters into the RFM model.

[0183] Equations were established for all control points in the manner described above and represented as matrices. The matrices were then combined to obtain the overall error equation.

[0184] Step 7.4: Normalize the total error equation and solve for the affine transformation parameters using the least squares method.

[0185] Step 7.5: Based on the solved affine transformation parameters, calculate the image-side coordinate error of each control point and calculate the mean error of the image-side coordinate error of the control points.

[0186] The image-side coordinate error of the g-th control point is:

[0187] The image-side coordinates of the g-th control point are calculated based on the solved affine transformation parameters.

[0188] The calculated row coordinates of the g-th control point are obtained based on the solved affine transformation parameters.

[0189] The column coordinates of the g-th control point are calculated based on the solved affine transformation parameters.

[0190] R g This represents the row coordinates of the g-th control point;

[0191] C g Represents the column coordinates of the g-th control point;

[0192] The formula for calculating the mean square error is as follows:

[0193]

[0194] N is the number of control points, and RMSE is the image square error.

[0195] The calculated row coordinates of the g-th control point are obtained based on the solved affine transformation parameters.

[0196] The column coordinates of the g-th control point are calculated based on the solved affine transformation parameters.

[0197] R g This represents the row coordinates of the g-th control point;

[0198] C g Represents the column coordinates of the g-th control point; Step 7.6: Based on the preset error threshold for the image coordinate error of the control points, control points whose image coordinate error is greater than the error threshold are removed; thus obtaining the accurate geographic 3D coordinates of the control points on the DOM and the pixel coordinates of the control points in the intensity map image.

[0199] The mean square error threshold of the image-side coordinate error of the preset control point is calculated as follows:

[0200] TH rmse =KTH×RMSE

[0201] Where KTH is the standard error threshold coefficient for the image-side coordinate error of the set control point; RMSE is the image-side standard error.

[0202] The present invention also provides a computer-readable medium storing a computer program executed by an electronic device, which, when run on the electronic device, performs the steps of the satellite standard product generation method.

[0203] The advantage of this invention lies in its full utilization of the high geometric positioning accuracy of the domestically produced C-band Gaofen-3 satellite. It uses Gaofen-3 global orthophoto maps with a planar accuracy better than 10m as control data, and matches the Luojia-2 imagery with it to obtain high-precision control points. Through high-precision geometric processing, the geometric positioning accuracy of the Luojia-2 Ka-SAR imagery is improved, effectively enhancing the planar geometric positioning accuracy from the original 650-950m range of the Luojia-2 strip pattern imagery to better than 5m. This provides support for the refined application of Luojia-2 imagery data in industries such as water conservancy and surveying. The method presented in this paper offers a new approach to the fusion of multi-source SAR data. Attached Figure Description

[0204] Figure 1 : Flowchart of the method according to an embodiment of the present invention.

[0205] Figure 2 The image above shows the result of overlaying the standard product of the Wujing Luojia No. 2 strip pattern image with the GF-3DOM base map in an embodiment of the present invention. Detailed Implementation

[0206] 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.

[0207] In specific implementation, the method proposed in the technical solution of this invention can be automatically executed by those skilled in the art using computer software technology. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of this invention and computer equipment including the computer program running the corresponding computer program, should also be within the protection scope of this invention.

[0208] The following describes a specific embodiment of the present invention, with reference to the accompanying drawings, which is a method for generating satellite standard products and a computer-readable medium, as follows:

[0209] Step 1: Acquire the original Ka-SAR single-view complex image, the auxiliary Rational Function Model (RFM) parameter file, the C-band Gaofen-3 DOM, and the Digital Elevation Model (DEM);

[0210] Step 2: Perform intensity map transformation on the original Ka-SAR single-look complex image to generate the corresponding intensity map image;

[0211] Step 2, which involves generating the corresponding intensity map image, is as follows:

[0212] The Ka-SAR single-look complex image is defined as

[0213] The intensity map image is defined as I intensity ;

[0214] The pixel values ​​corresponding to the intensity map in row p and column q are calculated as follows:

[0215]

[0216] in,

[0217] This represents the pixel value corresponding to the p-th row and q-th column of the first band in a Ka-SAR single-look complex image.

[0218] This represents the pixel value corresponding to the p-th row and q-th column of the second band in a Ka-SAR single-view complex image.

[0219] I intensity (p,q) represents the pixel value corresponding to the intensity map image in row p and column q;

[0220] Step 3: Obtain the pixel width and pixel height of the intensity map image, and divide the intensity map image into sub-blocks based on the pixel width and pixel height;

[0221] Step 3, which involves dividing the intensity map image into sub-blocks, is described in the following steps:

[0222] Set the number of molecular blocks in the column direction to nGX, and set the number of molecular blocks in the row direction to nGY;

[0223] Initial matching: nGX = 2, nGY = 15;

[0224] Fine matching: nGX=3, nGY=30;

[0225] Calculate the column spacing and row spacing of each sub-block separately, as follows:

[0226]

[0227] Where CenX represents the spacing of each sub-block in the column direction, CenY represents the spacing of each sub-block in the row direction, Width represents the pixel width of the intensity map image, Height represents the pixel height of the intensity map image, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction.

[0228] The sub-block in row j and column i is defined as:

[0229] CEL i,j

[0230] i∈[1,nGX], j∈[1,nGY]

[0231] Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction;

[0232] Set the sub-block size to Kernel;

[0233] Set to 2048 for initial matching and 1024 for fine matching;

[0234] Calculate the top-left pixel coordinates of the sub-block in row j and column i, as follows:

[0235]

[0236]

[0237]

[0238] Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the top-left pixel coordinates of the sub-block in row j and column i. The column coordinate represents the top-left pixel coordinate of the sub-block in row j and column i. The row coordinate represents the top-left pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block.

[0239] The geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in row j and column i are defined as follows:

[0240]

[0241] in, This represents the longitude of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the latitude of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the elevation of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column;

[0242] Calculate the top-right pixel coordinates of the sub-block in row j and column i, as follows:

[0243]

[0244]

[0245]

[0246] Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. Represents the top-right pixel coordinates of the sub-block in row j and column i; The column coordinate represents the top-right pixel coordinate of the sub-block in row j and column i; The row coordinate represents the top-right pixel coordinate of the sub-block in the j-th row and i-th column; Kernel represents the size of the sub-block.

[0247] Combining DEM-aided data and the RFM model, the 3D geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column are calculated using ray tracing, defined as:

[0248]

[0249] in, This represents the longitude of the geographic coordinates corresponding to the top right pixel coordinates of the sub-block in the j-th row and i-th column; This represents the latitude of the geographic coordinates corresponding to the top right pixel coordinates of the sub-block in row j and column i. This represents the elevation of the geographic coordinates corresponding to the top-right pixel coordinates of the sub-block in the j-th row and i-th column;

[0250] Calculate the bottom-left pixel coordinates of the sub-block in row j and column i, as follows:

[0251]

[0252]

[0253]

[0254] Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the bottom-left pixel coordinates of the sub-block in row j and column i; The column coordinate represents the bottom left pixel coordinate of the sub-block in the j-th row and i-th column; The row coordinate represents the bottom left pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block.

[0255] The geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in row j and column i are defined as follows:

[0256]

[0257] in, This represents the longitude of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the latitude of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the elevation of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column;

[0258] Calculate the bottom right pixel coordinates of the sub-block in row j and column i, as follows:

[0259]

[0260]

[0261]

[0262] Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the bottom right pixel coordinates of the sub-block in row j and column i. The column coordinate represents the bottom right pixel coordinate of the sub-block in the j-th row and i-th column; The row coordinate represents the bottom right pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block.

[0263] The geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i are defined as follows:

[0264]

[0265] in, This represents the longitude of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i. This represents the latitude of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i. This represents the elevation of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in the j-th row and i-th column;

[0266] The minimum longitude among the longitudes of the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is used as the minimum longitude of the sub-block in the j-th row and i-th column, denoted as .

[0267] The minimum longitude among the geographic coordinates corresponding to the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is taken as the maximum longitude of the sub-block in the j-th row and i-th column, denoted as .

[0268] The minimum latitude of the geographic sub-block in the j-th row and i-th column is calculated from the latitudes of its top-left, top-right, bottom-left, and bottom-right pixel coordinates. This minimum latitude is denoted as the minimum latitude of the geographic sub-block in the j-th row and i-th column.

[0269] The minimum value among the latitudes of the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is taken as the maximum latitude of the sub-block in the j-th row and i-th column, denoted as .

[0270] Will Construct the sub-block of geography in the j-th row and i-th column in a counter-clockwise order.

[0271] Step 4: Calculate the output resolution for geometric correction based on the intensity map image information;

[0272] Step 4, which calculates the output resolution for geometric correction based on the intensity map image information, is as follows:

[0273] Read the intensity map image resolution Res using the GDAL library;

[0274] Pre-set the output resolution sampling factor threshold T for geometric correction;

[0275] Set T=2 for initial matching; set T=3 for fine matching;

[0276] The output resolution of the geometric correction is FRes = T × Res;

[0277] Step 5: Based on each geographic sub-block, combined with DEM auxiliary data, perform local geometric correction on the intensity map image to obtain multiple geometrically corrected sub-blocks of the intensity map image. Based on the range of each sub-block, resample the pixel values ​​of the corresponding area on the DOM to obtain multiple resampled sub-blocks of the DOM.

[0278] Step 5 involves performing local geometric correction on the intensity map image to obtain multiple geometrically corrected sub-blocks of the intensity map image, as detailed below:

[0279] According to the sub-block geographic sub-block CEL in row j and column i i,j The coordinates of the four corresponding vertices For this region, local geometric correction is performed using an indirect method based on the output resolution FRES of the geometrically corrected image, combined with the intensity map image RFM model and auxiliary DEM data, to obtain the geometrically corrected sub-block of the intensity map image, defined as Ortho_CEL. i,j , i∈[1,nGX], j∈[1,nGY];

[0280] Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction;

[0281] Step 5 involves resampling the pixel values ​​of the corresponding region on the DOM to obtain multiple resampled sub-blocks of the DOM, as detailed below:

[0282] According to the sub-block geographic sub-block CEL in row j and column i i,j The coordinates of the four corresponding vertices For this region, bilinear interpolation is used to resample it in the DOM according to the output resolution FRes, resulting in a DOM resampled sub-block, defined as DOM_CEL. i,j , i∈[1,nGX], j∈[1,nGY];

[0283] Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction;

[0284] Step 6: Using each DOM resampled sub-block as the reference image and the geometrically corrected sub-block of the intensity map image as the registration image, a heterogeneous matching algorithm is used to obtain multiple sets of corresponding point pairs located at the geographic coordinates of the DOM resampled sub-block and the geographic coordinates of the geometrically corrected sub-block of the intensity map image; the geographic coordinates of the DOM resampled sub-block are combined with the auxiliary DEM and mapped to the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM; the geographic coordinates of each geometrically corrected sub-block of the intensity map image are mapped to the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image.

[0285] Step 6 describes mapping the geographic coordinates of the DOM resampled sub-blocks to the auxiliary DEM to obtain the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM.

[0286] The elevation of each DOM resampled sub-block is obtained from the DEM auxiliary data using bilinear interpolation based on the corresponding longitude and latitude positions of each DOM resampled sub-block.

[0287] Combine the geographic coordinates and elevation of each DOM resampled sub-block to construct the 3D geographic coordinates of the control points corresponding to each DOM resampled sub-block on the DOM.

[0288] Step 6 involves mapping the geographic coordinates of each geometrically corrected sub-block of the intensity map image to the corresponding pixel coordinates on the intensity map image:

[0289] The elevation of each geometrically corrected sub-block of the intensity map image is obtained from the DEM auxiliary data by bilinear interpolation based on the geographic coordinates of each geometrically corrected sub-block in the intensity map image at the corresponding longitude and latitude positions.

[0290] The geographic 3D coordinates of each geometrically corrected sub-block of the intensity map image are constructed by combining the geographic coordinates of each geometrically corrected sub-block of the intensity map image with the elevation of each geometrically corrected sub-block of the intensity map image.

[0291] By combining the geographic 3D coordinates of each geometrically corrected sub-block of the intensity map image, the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image are obtained by forward calculation based on the RFM model of the intensity map image.

[0292] Step 7: Use regional network adjustment to remove outliers by matching the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM with the pixel coordinates of the geographic coordinates of each geometrically corrected sub-block on the intensity map image, so as to obtain the accurate geographic 3D coordinates of the control points on the DOM and the pixel coordinates of the control points on the intensity map image.

[0293] Step 7, which involves using regional network adjustment to remove outliers by matching the geographic 3D coordinates of the control points corresponding to each DOM resampling sub-block on the DOM with the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image, specifically involves:

[0294] Step 7.1: RFM model construction;

[0295]

[0296] in, This represents the regularized pixel coordinates of the g-th control point;

[0297] Represents the regularized row coordinates of the g-th control point;

[0298] Represents the regularized column coordinates of the g-th control point;

[0299] Represents the regularized geographic coordinates of the g-th control point;

[0300] This represents the regularized geographic coordinates (longitude) of the g-th control point;

[0301] Y ng Represents the regularized geographic coordinates (latitude) of the g-th control point;

[0302] Represents the regularized geographic coordinate elevation of the g-th control point;

[0303] P1 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1;

[0304] P2 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1;

[0305] P3 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1;

[0306] P4 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1;

[0307] It can be expressed in functional form as Let it be F r g ;

[0308] It can be expressed in functional form as Let it be F c g ;

[0309] Regularized geographic coordinates of the g-th control point Regularized pixel coordinates of the g-th control point The geographic coordinates of the g-th control point (X) g ,Y g Z g ) and the pixel coordinates of the g-th control point (R g C g The correspondence is as follows:

[0310]

[0311]

[0312] in,

[0313] This represents the regularized geographic coordinates (longitude) of the g-th control point;

[0314] Represents the regularized geographic coordinates (latitude) of the g-th control point;

[0315] Represents the regularized geographic coordinate elevation of the g-th control point;

[0316] X g This represents the longitude of the geographic coordinates of the g-th control point;

[0317] Y g This represents the latitude and geographic coordinates of the g-th control point.

[0318] Z g This represents the geographic coordinates and elevation of the g-th control point;

[0319] X o The translation parameters representing the longitude regularization of geographic coordinates are obtained from the RFM parameter file in step 1.

[0320] Y o The translation parameter representing the latitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1;

[0321] Z o The translation parameters representing geographic coordinate elevation regularization are obtained from the RFM parameter file in step 1;

[0322] X s The scaling factor representing the longitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1;

[0323] Y s The scaling factor representing the latitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1.

[0324] Z s The scaling factor representing geographic coordinate elevation regularization is obtained from the RFM parameter file in step 1;

[0325] Represents the regularized row coordinates of the g-th control point;

[0326] Represents the regularized column coordinates of the g-th control point;

[0327] R g This represents the row coordinates of the g-th control point;

[0328] C g Represents the column coordinates of the g-th control point;

[0329] R o The translation parameters representing row coordinate regularization are obtained from the RFM parameter file in step 1;

[0330] C o The translation parameters representing column coordinate regularization are obtained from the RFM parameter file in step 1;

[0331] R sThe scaling factor for row coordinate regularization is obtained from the RFM parameter file in step 1.

[0332] C s The scaling factor for column coordinate regularization is obtained from the RFM parameter file in step 1.

[0333] Step 7.2: Constructing the image-side compensation model based on the affine transformation model;

[0334] When using the image-based compensation scheme to compensate for systematic errors in RFM, the regularized geographic coordinates of the g-th control point described by RFM are... The pixel coordinates of the g-th control point (R) g C g Revised to:

[0335]

[0336] in,

[0337] R g This represents the row coordinates of the g-th control point;

[0338] C g Represents the column coordinates of the g-th control point;

[0339] This represents the regularized geographic coordinates (longitude) of the g-th control point;

[0340] Y n g Represents the regularized geographic coordinates (latitude) of the g-th control point;

[0341] Represents the regularized geographic coordinate elevation of the g-th control point;

[0342] P1 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1;

[0343] P2 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1;

[0344] P3 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1;

[0345] P4 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1;

[0346] R o The translation parameters representing row coordinate regularization are obtained from the RFM parameter file in step 1;

[0347] C o The translation parameters representing column coordinate regularization are obtained from the RFM parameter file in step 1;

[0348] R s The scaling factor for row coordinate regularization is obtained from the RFM parameter file in step 1.

[0349] C s The scaling factor for column coordinate regularization is obtained from the RFM parameter file in step 1.

[0350] (ΔR g ,ΔC g ) represents the pixel coordinates (R) of the g-th control point. g C g The system error compensation value,

[0351] ΔR g The system error compensation value represents the row coordinate of the g-th control point;

[0352] ΔC g This represents the system error compensation value for the column coordinates of the g-th control point;

[0353] The calculation method is as follows;

[0354]

[0355] Where (e0,e1,e2,f0,f1,f2) are the affine transformation parameters of the intensity map image.

[0356] e0 is the first parameter of the row direction affine transformation, e1 is the second parameter of the row direction affine transformation, and e2 is the third parameter of the row direction affine transformation.

[0357] f0 is the first parameter of the column direction affine transformation, f1 is the second parameter of the column direction affine transformation, and f2 is the third parameter of the column direction affine transformation;

[0358] Step 7.3: Establish an error equation for each control point, and combine all error equations into a total error equation;

[0359] The error equation for the g-th control point is in the following form:

[0360]

[0361] Treating the affine transformation coefficients as unknowns, the error equation for the g-th control point is rewritten in matrix form:

[0362] V g =A g tl g

[0363] in,

[0364] This is the set of residual vectors of the row and column coordinate observations of the image of the g-th control point;

[0365] Let be the residual vector of the row coordinate observation values ​​of the image of the g-th control point;

[0366] The residual vector of the image column coordinate observation values ​​of the g-th control point;

[0367] This is the coefficient matrix corresponding to the g-th control point;

[0368] R g This represents the row coordinates of the g-th control point;

[0369] C g Represents the column coordinates of the g-th control point;

[0370] t=[Δe0 Δe1 Δe2 Δf0 Δf1 Δf2] T This represents the incremental vector of error compensation parameters for the intensity map image coordinate system.

[0371] Δe0 is the increment of the first systematic error compensation parameter in the row direction;

[0372] Δe1 is the increment of the second systematic error compensation parameter in the row direction;

[0373] Δe2 is the increment of the third systematic error compensation parameter in the row direction;

[0374] Δf0 is the increment of the first systematic error compensation parameter in the column direction; Δf1 is the increment of the second systematic error compensation parameter in the column direction;

[0375] Δf2 is the increment of the third systematic error compensation parameter in the column direction;

[0376] in,

[0377] The difference between the row coordinates of the g-th control point and the geographic coordinates of the g-th control point, calculated by substituting the approximate values ​​of the affine transformation parameters into the RFM model.

[0378] The difference between the column coordinates of the g-th control point and the geographic coordinates of the g-th control point, calculated by substituting the approximate values ​​of the affine transformation parameters into the RFM model.

[0379] Equations were established for all control points in the manner described above and represented as matrices. The matrices were then combined to obtain the overall error equation.

[0380] Step 7.4: Normalize the total error equation and solve for the affine transformation parameters using the least squares method.

[0381] Step 7.5: Based on the solved affine transformation parameters, calculate the image-side coordinate error of each control point and calculate the mean error of the image-side coordinate error of the control points.

[0382] The image-side coordinate error of the g-th control point is:

[0383] The image-side coordinates of the g-th control point are calculated based on the solved affine transformation parameters.

[0384] The calculated row coordinates of the g-th control point are obtained based on the solved affine transformation parameters.

[0385] The column coordinates of the g-th control point are calculated based on the solved affine transformation parameters.

[0386] R g This represents the row coordinates of the g-th control point;

[0387] C g Represents the column coordinates of the g-th control point;

[0388] The formula for calculating the mean square error is as follows:

[0389]

[0390] N is the number of control points, and RMSE is the image square error.

[0391] The calculated row coordinates of the g-th control point are obtained based on the solved affine transformation parameters.

[0392] The column coordinates of the g-th control point are calculated based on the solved affine transformation parameters.

[0393] R g This represents the row coordinates of the g-th control point;

[0394] C g Represents the column coordinates of the g-th control point; Step 7.6: Based on the preset error threshold for the image coordinate error of the control points, control points whose image coordinate error is greater than the error threshold are removed; thus obtaining the accurate geographic 3D coordinates of the control points on the DOM and the pixel coordinates of the control points in the intensity map image.

[0395] The mean square error threshold of the image-side coordinate error of the preset control point is calculated as follows:

[0396] TH rmse =KTH×RMSE

[0397] Where KTH is the standard error threshold coefficient for the image-side coordinate error of the set control point; RMSE is the image-side standard error.

[0398] Initial matching KTH=3; Fine matching KTH=2;

[0399] Step 8: Update the RFM model using the accurate geographic coordinates of control points on the DOM and the pixel coordinates of control points in the intensity map image to obtain the updated RFM model;

[0400] Step 9: Replace the RFM model of the intensity map image with the accurate RFM model updated in Step 8, and repeat Steps 3-8 to obtain the final RFM model.

[0401] Step 10: Using the final RFM model and auxiliary DEM data, the intensity map image is indirectly corrected to obtain the standard product of Ka-SAR satellite imagery.

[0402] Table 1. Detailed information on Luojia-2 strip pattern images.

[0403] Tianjin 2023.07.31 Ascend 4065×211250 39.3379°N, 116.5865°E Dongying 2023.07.24 Ascend 4237×187500 37.2733°N, 118.4765°E Wuhan 2023.07.18 Ascend 4661×237150 30.4455°N, 114.197°E Egypt 2023.12.05 Descending orbit 3440×111400 31.2345°N, 32.2427°E Brazil 2023.10.30 Descending orbit 4858×226400 -20.3799°N, -51.3707°E

[0404] Luojia-2 01 is a Ka-SAR microwave remote sensing satellite developed by Wuhan University. This satellite can acquire multi-angle, three-dimensional information; achieve Ka-band video SAR imaging; and can be applied to the monitoring and calculation of various water conservancy and hydropower elements such as water flow velocity, discharge, area, and deformation. The experiment selected five Luojia-2 strip mode images for verification, as shown in Table 1. Control points were matched between the Luojia-2 Ka-SAR images and the GF-3DOM. Subsequently, the positioning accuracy of the Luojia-2 strip mode images was verified and improved based on the control points, as shown in Table 2.

[0405] Table 2. Original localization and improved results of Luojia-2 strip imagery using this method.

[0406]

[0407] The zero-control checkpoint planar accuracy of the five imagery scenes ranges from 650 to 950 meters. After processing using the method described in this paper, the high-precision control points matched, using a 4-control configuration, achieve a checkpoint planar accuracy better than 5 meters, demonstrating a significant improvement in geometric positioning accuracy. This reflects the effectiveness of the proposed method and can provide support for high-precision processing of the Luojia-2 satellite. The overlay results of the standard product of the five Luojia-2 strip pattern imagery scenes with the GF-3 DOM basemap are shown below. Figure 2 As shown.

[0408] The computer-readable medium is a server workstation;

[0409] The server workstation stores computer programs executed by the electronic device. When the computer program runs on the electronic device, it causes the electronic device to execute the Ka-SAR satellite standard product generation method and steps described in this embodiment of the invention.

[0410] It should be understood that any parts not described in detail in this specification belong to the prior art.

[0411] It should be understood that the above description of the embodiments is quite detailed, but it should not be considered as a limitation on the scope of protection of this invention. Those skilled in the art can make substitutions or modifications under the guidance of this invention without departing from the scope of protection of the claims of this invention, and all such substitutions or modifications fall within the scope of protection of this invention. The scope of protection of this invention should be determined by the appended claims.

Claims

1. A method for generating standard Ka-SAR satellite products, characterized in that, include: Perform intensity map conversion to generate the corresponding intensity map image; The intensity map image is divided into sub-blocks; Calculate the output resolution for geometric correction based on the intensity map image information; Local geometric correction is performed on the intensity map image to obtain multiple geometrically corrected sub-blocks of the intensity map image. Based on the range of each sub-block, the pixel values ​​of the corresponding region are resampled on the DOM to obtain multiple resampled sub-blocks of the DOM. A heterogeneous matching algorithm is used to obtain the geographic coordinates of multiple pairs of points located in the DOM resampled sub-block and the geographic coordinates of the intensity map image geometrically corrected sub-block; the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block and the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image are then mapped sequentially. Gross errors were eliminated by using regional network adjustment to obtain accurate 3D geographic coordinates of control points on the DOM and pixel coordinates of control points in intensity map imagery. The RFM model is updated by using the precise geographic coordinates of control points on the DOM and the pixel coordinates of control points in the intensity map image. The final RFM model is obtained through iterative updates. Using the final RFM model and auxiliary DEM data, the intensity map image is indirectly corrected to obtain the standard product of Ka-SAR satellite imagery.

2. The method for generating satellite standard products according to claim 1, characterized in that, Includes the following steps: Step 1: Acquire the original Ka-SAR single-view complex image, auxiliary rational function model parameter file, C-band Gaofen-3 DOM, and digital elevation model; Step 2: Perform intensity map transformation on the original Ka-SAR single-look complex image to generate the corresponding intensity map image; Step 3: Obtain the pixel width and pixel height of the intensity map image, and divide the intensity map image into sub-blocks based on the pixel width and pixel height; Step 4: Calculate the output resolution for geometric correction based on the intensity map image information; Step 5: Based on each geographic sub-block, combined with DEM auxiliary data, perform local geometric correction on the intensity map image to obtain multiple geometrically corrected sub-blocks of the intensity map image. Based on the range of each sub-block, resample the pixel values ​​of the corresponding area on the DOM to obtain multiple resampled sub-blocks of the DOM. Step 6: Using each DOM resampled sub-block as the reference image and the geometrically corrected sub-block of the intensity map image as the registration image, a heterogeneous matching algorithm is used to obtain multiple sets of corresponding point pairs located at the geographic coordinates of the DOM resampled sub-block and the geographic coordinates of the geometrically corrected sub-block of the intensity map image; the geographic coordinates of the DOM resampled sub-block are combined with the auxiliary DEM and mapped to the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM; the geographic coordinates of each geometrically corrected sub-block of the intensity map image are mapped to the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image. Step 7: Use regional network adjustment to remove outliers by matching the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM with the pixel coordinates of the geographic coordinates of each geometrically corrected sub-block on the intensity map image, so as to obtain the accurate geographic 3D coordinates of the control points on the DOM and the pixel coordinates of the control points on the intensity map image. Step 8: Update the RFM model using the accurate geographic coordinates of control points on the DOM and the pixel coordinates of control points in the intensity map image to obtain the updated RFM model; Step 9: Replace the RFM model of the intensity map image with the updated and accurate RFM model from step 8, and repeat steps 3-8 to obtain the final RFM model. Step 10: Using the final RFM model and auxiliary DEM data, the intensity map image is indirectly corrected to obtain the standard product of Ka-SAR satellite imagery.

3. The method for generating satellite standard products according to claim 2, characterized in that, Step 2, which involves generating the corresponding intensity map image, is as follows: The Ka-SAR single-look complex image is defined as The intensity map image is defined as I intensity ; The pixel values ​​corresponding to the intensity map in row p and column q are calculated as follows: in, This represents the pixel value corresponding to the p-th row and q-th column of the first band in a Ka-SAR single-look complex image. This represents the pixel value corresponding to the p-th row and q-th column of the second band in a Ka-SAR single-view complex image. I intensity (p,q) represents the pixel value corresponding to the intensity map image in row p and column q.

4. The method for generating satellite standard products according to claim 3, characterized in that: Step 3, which involves dividing the intensity map image into sub-blocks, is described in the following steps: Set the number of molecular blocks in the column direction to nGX, and set the number of molecular blocks in the row direction to nGY; Calculate the column spacing and row spacing of each sub-block separately, as follows: Where CenX represents the spacing of each sub-block in the column direction, CenY represents the spacing of each sub-block in the row direction, Width represents the pixel width of the intensity map image, Height represents the pixel height of the intensity map image, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction. The sub-block in row j and column i is defined as: THE i,j i∈[1,nGX], j∈[1,nGY] Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction; Set the sub-block size to Kernel; Calculate the top-left pixel coordinates of the sub-block in row j and column i, as follows: Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the top-left pixel coordinates of the sub-block in row j and column i. The column coordinate represents the top-left pixel coordinate of the sub-block in row j and column i. The row coordinate represents the top-left pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block. The geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in row j and column i are defined as follows: in, This represents the longitude of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the latitude of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the elevation of the geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column; Calculate the top-right pixel coordinates of the sub-block in row j and column i, as follows: Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. Represents the top-right pixel coordinates of the sub-block in row j and column i; The column coordinate represents the top-right pixel coordinate of the sub-block in row j and column i; The row coordinate represents the top-right pixel coordinate of the sub-block in the j-th row and i-th column; Kernel represents the size of the sub-block. Combining DEM-aided data and the RFM model, the 3D geographic coordinates corresponding to the top-left pixel coordinates of the sub-block in the j-th row and i-th column are calculated using ray tracing, defined as: in, This represents the longitude of the geographic coordinates corresponding to the top right pixel coordinates of the sub-block in the j-th row and i-th column; This represents the latitude of the geographic coordinates corresponding to the top right pixel coordinates of the sub-block in row j and column i. This represents the elevation of the geographic coordinates corresponding to the top-right pixel coordinates of the sub-block in the j-th row and i-th column; Calculate the bottom-left pixel coordinates of the sub-block in row j and column i, as follows: Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the bottom-left pixel coordinates of the sub-block in row j and column i; The column coordinate represents the bottom left pixel coordinate of the sub-block in the j-th row and i-th column; The row coordinate represents the bottom left pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block. The geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in row j and column i are defined as follows: in, This represents the longitude of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the latitude of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column; This represents the elevation of the geographic coordinates corresponding to the bottom left pixel coordinates of the sub-block in the j-th row and i-th column; Calculate the bottom right pixel coordinates of the sub-block in row j and column i, as follows: Where CenX represents the column spacing of each sub-block, and CenY represents the row spacing of each sub-block. This represents the bottom right pixel coordinates of the sub-block in row j and column i. The column coordinate represents the bottom right pixel coordinate of the sub-block in row j and column i; The row coordinate represents the bottom right pixel coordinate of the sub-block in the j-th row and i-th column, and Kernel represents the size of the sub-block. The geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i are defined as follows: in, This represents the longitude of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i. This represents the latitude of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in row j and column i. This represents the elevation of the geographic coordinates corresponding to the bottom right pixel coordinates of the sub-block in the j-th row and i-th column; The minimum longitude among the longitudes of the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is used as the minimum longitude of the sub-block in the j-th row and i-th column, denoted as . The minimum longitude among the geographic coordinates corresponding to the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is taken as the maximum longitude of the sub-block in the j-th row and i-th column, denoted as . The minimum latitude of the geographic sub-block in the j-th row and i-th column is calculated from the latitudes of its top-left, top-right, bottom-left, and bottom-right pixel coordinates. This minimum latitude is denoted as the minimum latitude of the geographic sub-block in the j-th row and i-th column. The minimum value among the latitudes of the top-left, top-right, bottom-left, and bottom-right pixels of the sub-block in the j-th row and i-th column is taken as the maximum latitude of the sub-block in the j-th row and i-th column, denoted as . Will Construct the sub-block of geography in the j-th row and i-th column in a counter-clockwise order.

5. The method for generating satellite standard products according to claim 4, characterized in that: Step 4, which calculates the output resolution for geometric correction based on the intensity map image information, is as follows: Read the intensity map image resolution Res using the GDAL library; Pre-set the output resolution sampling factor threshold T for geometric correction; The output resolution of the geometric correction is FRes = T × Res.

6. The method for generating satellite standard products according to claim 5, characterized in that: Step 5 involves performing local geometric correction on the intensity map image to obtain multiple geometrically corrected sub-blocks of the intensity map image, as detailed below: According to the sub-block geographic sub-block CEL in row j and column i i,j The coordinates of the four corresponding vertices For this region, local geometric correction is performed using an indirect method based on the output resolution FRES of the geometrically corrected image, combined with the intensity map image RFM model and auxiliary DEM data, to obtain the geometrically corrected sub-block of the intensity map image, defined as Ortho_CEL. i,j , i∈[1,nGX], j∈[1,nGY]; Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction.

7. The method for generating satellite standard products according to claim 6, characterized in that: Step 5 involves resampling the pixel values ​​of the corresponding region on the DOM to obtain multiple resampled sub-blocks of the DOM, as detailed below: According to the sub-block geographic sub-block CEL in row j and column i i,j The coordinates of the four corresponding vertices For this region, bilinear interpolation is used to resample it in the DOM according to the output resolution FRes, resulting in a DOM resampled sub-block, defined as DOM_CEL. i,j , i∈[1,nGX], j∈[1,nGY]; Where i represents the sub-block number in the column direction, j represents the sub-block number in the row direction, nGX is the set number of sub-blocks in the column direction, and nGY is the set number of sub-blocks in the row direction.

8. The method for generating satellite standard products according to claim 7, characterized in that: Step 6 describes mapping the geographic coordinates of the DOM resampled sub-blocks to the auxiliary DEM to obtain the geographic 3D coordinates of the control points corresponding to each DOM resampled sub-block on the DOM. The elevation of each DOM resampled sub-block is obtained from the DEM auxiliary data using bilinear interpolation based on the corresponding longitude and latitude positions of each DOM resampled sub-block. Combine the geographic coordinates and elevation of each DOM resampled sub-block to construct the 3D geographic coordinates of the control points corresponding to each DOM resampled sub-block on the DOM. Step 6 involves mapping the geographic coordinates of each geometrically corrected sub-block of the intensity map image to the corresponding pixel coordinates on the intensity map image: The elevation of each geometrically corrected sub-block of the intensity map image is obtained from the DEM auxiliary data by bilinear interpolation based on the geographic coordinates of each geometrically corrected sub-block in the intensity map image at the corresponding longitude and latitude positions. The geographic 3D coordinates of each geometrically corrected sub-block of the intensity map image are constructed by combining the geographic coordinates of each geometrically corrected sub-block of the intensity map image with the elevation of each geometrically corrected sub-block of the intensity map image. By combining the geographic 3D coordinates of each geometrically corrected sub-block of the intensity map image, the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image are obtained by forward calculation based on the RFM model of the intensity map image.

9. The method for generating satellite standard products according to claim 8, characterized in that: Step 7, which involves using regional network adjustment to remove outliers by matching the geographic 3D coordinates of the control points corresponding to each DOM resampling sub-block on the DOM with the pixel coordinates corresponding to the geographic coordinates of each geometrically corrected sub-block on the intensity map image, specifically involves: Step 7.1: RFM model construction; in, This represents the regularized pixel coordinates of the g-th control point; Represents the regularized row coordinates of the g-th control point; Represents the regularized column coordinates of the g-th control point; Represents the regularized geographic coordinates of the g-th control point; This represents the regularized geographic coordinates (longitude) of the g-th control point; Represents the regularized geographic coordinates (latitude) of the g-th control point; Represents the regularized geographic coordinate elevation of the g-th control point; P1 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1; P2 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1; P3 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1; P4 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1; It can be expressed in functional form as Recorded as It can be expressed in functional form as Recorded as Regularized geographic coordinates of the g-th control point Regularized pixel coordinates of the g-th control point The geographic coordinates of the g-th control point (X) g ,Y g Z g ) and the pixel coordinates of the g-th control point (R g C g The correspondence is as follows: in, This represents the regularized geographic coordinates (longitude) of the g-th control point; Represents the regularized geographic coordinates (latitude) of the g-th control point; Represents the regularized geographic coordinate elevation of the g-th control point; X g This represents the longitude of the geographic coordinates of the g-th control point; Y g This represents the latitude and geographic coordinates of the g-th control point. Z g This represents the geographic coordinates and elevation of the g-th control point; X o The translation parameters representing the longitude regularization of geographic coordinates are obtained from the RFM parameter file in step 1. Y o The translation parameter representing the latitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1; Z o The translation parameters representing geographic coordinate elevation regularization are obtained from the RFM parameter file in step 1; X s The scaling factor representing the longitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1; Y s The scaling factor representing the latitude regularization of geographic coordinates is obtained from the RFM parameter file in step 1. Z s The scaling factor representing geographic coordinate elevation regularization is obtained from the RFM parameter file in step 1; Represents the regularized row coordinates of the g-th control point; Represents the regularized column coordinates of the g-th control point; R g This represents the row coordinates of the g-th control point; C g Represents the column coordinates of the g-th control point; R o The translation parameters representing row coordinate regularization are obtained from the RFM parameter file in step 1; C o The translation parameters representing column coordinate regularization are obtained from the RFM parameter file in step 1; R s The scaling factor for row coordinate regularization is obtained from the RFM parameter file in step 1. C s The scaling factor for column coordinate regularization is obtained from the RFM parameter file in step 1. Step 7.2: Constructing the image-side compensation model based on the affine transformation model; When using the image-based compensation scheme to compensate for systematic errors in RFM, the regularized geographic coordinates of the g-th control point described by RFM are... The pixel coordinates of the g-th control point (R) g C g Revised to: in, R g This represents the row coordinates of the g-th control point; C g Represents the column coordinates of the g-th control point; This represents the regularized geographic coordinates (longitude) of the g-th control point; Represents the regularized geographic coordinates (latitude) of the g-th control point; Represents the regularized geographic coordinate elevation of the g-th control point; P1 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1; P2 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized row coordinates, obtained from the RFM parameter file in step 1; P3 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1; P4 is the three-stage polynomial coefficient that connects the regularized geographic coordinates and the regularized column coordinates, obtained from the RFM parameter file in step 1; R o The translation parameters representing row coordinate regularization are obtained from the RFM parameter file in step 1; C o The translation parameters representing column coordinate regularization are obtained from the RFM parameter file in step 1; R s The scaling factor for row coordinate regularization is obtained from the RFM parameter file in step 1. C s The scaling factor for column coordinate regularization is obtained from the RFM parameter file in step 1. (ΔR g ,ΔC g ) represents the pixel coordinates (R) of the g-th control point. g C g The system error compensation value, ΔR g The system error compensation value represents the row coordinate of the g-th control point; ΔC g This represents the system error compensation value for the column coordinates of the g-th control point; The calculation method is as follows; Where (e0,e1,e2,f0,f1,f2) are the affine transformation parameters of the intensity map image; e0 is the first parameter of the row direction affine transformation, e1 is the second parameter of the row direction affine transformation, and e2 is the third parameter of the row direction affine transformation. f0 is the first parameter of the column direction affine transformation, f1 is the second parameter of the column direction affine transformation, and f2 is the third parameter of the column direction affine transformation; Step 7.3: Establish an error equation for each control point, and combine all error equations into a total error equation; The error equation for the g-th control point is in the following form: Treating the affine transformation coefficients as unknowns, the error equation for the g-th control point is rewritten in matrix form: V g =A g t-l g in, V g =[v Rg v Cg ] T This is the set of residual vectors of the row and column coordinate observations of the image of the g-th control point; v Rg Let be the residual vector of the row coordinate observation values ​​of the image of the g-th control point; v Cg The residual vector of the image column coordinate observation values ​​of the g-th control point; This is the coefficient matrix corresponding to the g-th control point; R g This represents the row coordinates of the g-th control point; C g Represents the column coordinates of the g-th control point; t=[Δe0 Δe1 Δe2 Δf0 Δf1 Δf2] T This represents the incremental vector of error compensation parameters for the intensity map image coordinate system; Δe0 is the increment of the first systematic error compensation parameter in the row direction; Δe1 is the increment of the second systematic error compensation parameter in the row direction; Δe2 is the increment of the third systematic error compensation parameter in the row direction; Δf0 is the increment of the first systematic error compensation parameter in the column direction; Δf1 is the increment of the second systematic error compensation parameter in the column direction; Δf2 is the increment of the third systematic error compensation parameter in the column direction; in, The difference between the row coordinates of the g-th control point and the geographic coordinates of the g-th control point, calculated by substituting the approximate values ​​of the affine transformation parameters into the RFM model; The difference between the column coordinates of the g-th control point and the geographic coordinates of the g-th control point, calculated by substituting the approximate values ​​of the affine transformation parameters into the RFM model; Equations were established for all control points in the manner described above and represented as matrices. The matrices were then combined to obtain the overall error equation. Step 7.4: Normalize the total error equation and solve for the affine transformation parameters using the least squares method. Step 7.5: Based on the solved affine transformation parameters, calculate the image-side coordinate error of each control point and calculate the mean error of the image-side coordinate error of the control points. The image-side coordinate error of the g-th control point is: The image-side coordinates of the g-th control point are calculated based on the solved affine transformation parameters. The calculated row coordinates of the g-th control point are obtained based on the solved affine transformation parameters. The column coordinates of the g-th control point are calculated based on the solved affine transformation parameters. R g This represents the row coordinates of the g-th control point; C g Represents the column coordinates of the g-th control point; The formula for calculating the mean square error is as follows: N is the number of control points, and RMSE is the image square error. The calculated row coordinates of the g-th control point are obtained based on the solved affine transformation parameters. The column coordinates of the g-th control point are calculated based on the solved affine transformation parameters. R g This represents the row coordinates of the g-th control point; C g Represents the column coordinates of the g-th control point; Step 7.6: Based on the preset error threshold for the image coordinate error of the control points, control points whose image coordinate error is greater than the error threshold are removed; thus obtaining the accurate geographic 3D coordinates of the control points on the DOM and the pixel coordinates of the control points in the intensity map image. The mean square error threshold of the image-side coordinate error of the preset control point is calculated as follows: TH rmse =KTH×RMSE Where KTH is the threshold coefficient for the image-side coordinate error of the set control point; RMSE is the image-side mean square error.

10. A computer-readable medium, characterized in that, It stores a computer program executed by an electronic device, which, when run on the electronic device, causes the electronic device to perform the steps of the method as described in any one of claims 1-9.

Images