Satellite image geometric correction method based on low-frequency resampling

Abstract

The invention provides a satellite image geometric correction method based on low-frequency resampling, and the method comprises the steps: obtaining basic data of an existing L1-level synthetic aperture radar SAR image product, and calculating a universal horizontal axis Mercator UTM projection range corresponding to an L1-level SAR image according to the basic data; calculating a rotation angle of an observation area of the L1-level SAR image in a minimum enclosing rectangle according to the UTM projection range, equally dividing an observation beam range based on the basic data, and calculating a scaling corresponding to each equally divided interval; and carrying out low-frequency resampling on the pixel points of each equally-divided interval according to the corresponding scaling, and rotating the resampled image according to the rotation angle to obtain an L2-level SAR image product after geometric correction. According to the method, the time required by the geometric correction of the SAR image is greatly shortened, the calculation complexity is reduced, and the convenience and applicability of the geometric correction of the satellite image are improved.

Description

Technical Field

[0001] The present application relates to the technical field of satellite image processing, and in particular to a satellite image geometric correction method based on low-frequency resampling. Background Art

[0002] With the rapid development of aerospace technology, the launch cost of small satellites has been reduced year by year, and a large number of small satellites have been successfully equipped with artificial intelligence (AI) real-time processors, which can realize the analysis, imaging, radiation correction, geometric correction, target detection and ground object classification of synthetic aperture radar (SAR) satellite in-orbit data.

[0003] In the related art, the traditional SAR image geometric correction method usually needs to establish a rational polynomial coefficient model (RPC) based on anchored control points and digital elevation models (DEM) on the ground for correction. Although this type of correction method can correct the distance image distortion (for example, compression or stretching) of L1-level SAR images, the computational complexity of this type of correction method is high, the correction time is long, and the preset ground control points are required, which consumes high computing resources. Therefore, the above-mentioned SAR image geometric correction method is not suitable for deployment on the real-time AI processor of the on-orbit SAR satellite with limited computing resources. Summary of the invention

[0004] The present application aims to solve one of the technical problems in the related art at least to some extent.

[0005] To this end, the first purpose of this application is to propose a satellite image geometric correction method based on low-frequency resampling, which can avoid the calculation of RPC parameters, quickly implement SAR image geometric correction, greatly reduce the time required for SAR image geometric correction, reduce the calculation complexity, and improve the convenience and applicability of satellite image geometric correction.

[0006] The second object of this application is to propose a satellite image geometric correction system based on low-frequency resampling;

[0007] A third object of the present application is to provide a non-transitory computer-readable storage medium.

[0008] To achieve the above-mentioned purpose, a first embodiment of the present application is to propose a satellite image geometric correction method based on low-frequency resampling, the method comprising the following steps:

[0009] Obtaining basic data of an existing L1-level synthetic aperture radar SAR image product, and calculating a Universal Transverse Mercator (UTM) projection range corresponding to the L1-level SAR image according to the basic data;

[0010] Calculating the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle according to the UTM projection range, and dividing the observation beam range into equal parts based on the basic data, and calculating the scaling ratio corresponding to each equal-divided interval;

[0011] The pixel points between each of the equal partitions are low-frequency resampled according to the corresponding scaling ratio, and the resampled image is rotated according to the rotation angle to obtain a geometrically corrected L2-level SAR image product.

[0012] Optionally, in one embodiment of the present application, calculating the Universal Transverse Mercator UTM projection range corresponding to the L1-level SAR image based on the basic data includes: extracting the longitude and latitude coordinates of the four corner points of the observation area from the basic data, and calculating the longitude and latitude range of the minimum circumscribed rectangle of the observation area based on the longitude and latitude coordinates of the four corner points; mapping the longitude and latitude range of the minimum circumscribed rectangle to the UTM projection coordinates of the four boundary points of the minimum circumscribed rectangle; and calculating the actual ground lengths in the distance and azimuth directions respectively based on the UTM projection coordinates of the four boundary points.

[0013] Optionally, in one embodiment of the present application, the calculation of the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle includes: determining the angle between the satellite orbit corresponding to the L1-level SAR image and the prime meridian; when the angle is a positive value, calculating a first difference based on the UTM projection coordinates of the left boundary point and the UTM projection coordinates of the lower boundary point, and calculating the rotation angle based on the first difference and the actual ground length in the distance direction; when the angle is a negative value, calculating a second difference based on the UTM projection coordinates of the lower boundary point and the UTM projection coordinates of the right boundary point, and calculating the rotation angle based on the second difference and the actual ground length in the distance direction.

[0014] Optionally, in one embodiment of the present application, the observation beam range is divided into equal parts based on the basic data, and the scaling ratio corresponding to each equal-division interval is calculated, including: calculating the observation beam range according to the near-point observation incident angle and the far-point observation incident angle extracted from the basic data, and dividing the observation beam range into equal parts according to a preset equal-division value; calculating the ground distance set corresponding to each equal-division interval in combination with the SAR satellite earth observation model and the basic data; calculating the ratio set of the distance to the total pixel corresponding to each equal-division interval according to the ground distance set; calculating the ground distance of the central interval between each equal-division interval, and calculating the scaling ratio set corresponding to each equal-division interval according to the ground distance set and the ground distance of the central interval.

[0015] Optionally, in one embodiment of the present application, the SAR satellite earth observation model and the basic data are combined to calculate the set of ground distances corresponding to each equally divided interval, including: calculating the perigee electromagnetic wave transmission distance and the apogee electromagnetic wave transmission distance in the SAR satellite earth observation model based on the basic data and the angles between each equally divided interval; calculating the ground distance corresponding to the current equally divided interval according to the perigee electromagnetic wave transmission distance, the apogee electromagnetic wave transmission distance, the angle between the current equally divided interval, and the angle of the previous interval adjacent to the current equally divided interval.

[0016] Optionally, in one embodiment of the present application, calculating the ratio set of distances to total pixels corresponding to each equally divided interval based on the ground distance set includes: calculating the sum of all ground distances in the ground distance set; dividing the ground distance corresponding to each of the equally divided intervals by the sum, to obtain the total distance pixels corresponding to each of the equally divided intervals; calculating the scaling ratio set corresponding to each of the equally divided intervals includes: dividing the ground distance corresponding to each of the equally divided intervals by the ground distance of the central interval, to obtain the scaling ratio corresponding to each of the equally divided intervals.

[0017] Optionally, in one embodiment of the present application, the low-frequency resampling of the pixel points between each of the equally divided areas according to the corresponding scaling ratio includes: calculating a set of pixel points between equally divided areas based on a set of ratios of the total pixels in the distance direction, the number of points in the azimuth direction between each equally divided area, and the number of pixel points in the distance direction between each equally divided area; and stretching the pixels between the near point interval and compressing the pixels between the far point interval in combination with the set of pixel points between the equally divided areas and the scaling ratio set.

[0018] To achieve the above purpose, the second embodiment of the present application further proposes a satellite image geometric correction system based on low-frequency resampling, comprising the following modules:

[0019] A calculation module is used to obtain basic data of existing L1-level synthetic aperture radar SAR image products, and calculate the universal transverse Mercator UTM projection range corresponding to the L1-level SAR image according to the basic data;

[0020] An equal division module, used for calculating the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle according to the UTM projection range, and equally dividing the observation beam range based on the basic data, and calculating the scaling ratio corresponding to each equally divided interval;

[0021] The correction module is used to perform low-frequency resampling of the pixel points between each of the equal partitions according to the corresponding scaling ratio, and rotate the resampled image according to the rotation angle to obtain a geometrically corrected L2-level SAR image product.

[0022] In order to implement the above-mentioned embodiments, the third aspect embodiment of the present application also proposes a non-temporary computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, the satellite image geometric correction method based on low-frequency resampling in the above-mentioned first aspect embodiment is implemented.

[0023] The technical solution provided by the embodiment of the present application brings at least the following beneficial effects: the present application first calculates the UTM projection range corresponding to the current image based on the basic information in the L1-level SAR image product; then calculates the rotation angle of the observation area in the minimum circumscribed rectangle according to the UTM projection range; then divides the interval of the effective observation beam range into N equal parts based on the basic information, and calculates the scaling ratio of the observation pixels corresponding to each equal-divided interval compared to the central observation point; finally, resamples the pixel points between each equal-divided interval according to the scaling ratio, rotates the resampled image according to the rotation angle, and obtains the corrected output result. Thus, the present application can avoid the steps of calculating RPC parameters, quickly implement SAR image geometric correction, and greatly reduce the time required for SAR image geometric correction on the basis of ensuring the accuracy of SAR image geometric correction, reduce the computational complexity of geometric correction and the consumed computing resources, and can be applied to the on-orbit data processing tasks of SAR satellite AI real-time processors, thereby improving the convenience and applicability of satellite image geometric correction.

[0024] Additional aspects and advantages of the present invention will be given in part in the following description and in part will be obvious from the following description, or will be learned through practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] The above and / or additional aspects and advantages of the present application will become apparent and easily understood from the following description of the embodiments in conjunction with the accompanying drawings, in which:

[0026] Figure 1A flowchart of a satellite image geometric correction method based on low-frequency resampling proposed in an embodiment of the present application;

[0027] Figure 2 A flowchart of a method for calculating a UTM projection range proposed in an embodiment of the present application;

[0028] Figure 3 A flowchart of a method for equally dividing a beam range and calculating a scaling ratio proposed in an embodiment of the present application;

[0029] Figure 4 A schematic diagram of a SAR satellite earth observation model proposed in an embodiment of the present application;

[0030] Figure 5 A schematic diagram of a geometric correction process of an L1-level SAR image proposed in an embodiment of the present application;

[0031] Figure 6 A comparative schematic diagram of geometric correction of L1-level data of a satellite proposed in an embodiment of the present application;

[0032] Figure 7 A comparative schematic diagram of another method of performing geometric correction on L1-level data of a satellite proposed in an embodiment of the present application;

[0033] Figure 8 A schematic diagram of the structure of a satellite image geometric correction system based on low-frequency resampling proposed in an embodiment of the present application. DETAILED DESCRIPTION

[0034] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, and should not be construed as limiting the present invention.

[0035] It should be noted that the present application proposes a satellite image geometric correction method based on low-frequency resampling for the on-orbit geometric correction task of SAR images. This task is an important part of the on-orbit data processing of SAR satellites. Generally, the on-orbit data processing flow of SAR satellites includes the following steps: 1) The ground transmitting station systematically plans the SAR satellite mission and uploads the mission instructions; 2) After the SAR satellite moves to the mission area, the SAR satellite payload is turned on and the target area is observed; 3) The SAR satellite antenna receives the echo data and stores the original echo data in the AI ​​real-time processor of the SAR satellite through the data transmission module; 4) The AI ​​real-time processor performs the data parsing task and obtains the L0-level data parsing result based on the original echo data; 5) The AI ​​real-time processor executes the on-orbit imaging instruction and generates the L1-level SAR image product (without geometric correction) based on the L0-level data; 6) The AI ​​real-time processor executes the on-orbit radiation correction and geometric correction instructions and generates the L2-level corrected SAR image based on the L1-level imaging product; 7) The AI ​​real-time processor performs downstream tasks such as target detection and object classification and generates the L3-level target-level product based on the L2-level corrected SAR image. This application mainly focuses on step 6 in the above-mentioned SAR satellite on-orbit data processing flow.

[0036] However, the satellite image geometric correction method in the related art has high computational complexity, takes a long time to correct, and requires preset ground control points, which consumes high computing resources and is not suitable for deployment on on-orbit SAR satellite real-time AI processors with limited computing resources. To this end, this application proposes a satellite image geometric correction method based on low-frequency resampling. This method is based on a variety of basic information in L1-level SAR imaging products, and implements the strategy of maintaining azimuth, stretching the distance to the near point, and compressing the distance to the far point, to achieve rapid geometric correction of L1-level SAR images. This method is suitable for on-orbit data processing tasks of SAR satellite AI real-time processors, and greatly reduces the geometric correction time of SAR images without significantly reducing the geometric correction accuracy.

[0037] A satellite image geometric correction method and system based on low-frequency resampling proposed in an embodiment of the present invention will be described below with reference to the accompanying drawings.

[0038] Figure 1 A flowchart of a satellite image geometric correction method based on low-frequency resampling proposed in an embodiment of the present application, such as Figure 1 As shown, the method comprises the following steps:

[0039] Step S101, obtaining basic data of an existing L1-level synthetic aperture radar SAR image product, and calculating a Universal Transverse Mercator UTM projection range corresponding to the L1-level SAR image according to the basic data.

[0040] Specifically, the present application first obtains the uncorrected L1-level SAR image product generated in step 5 of the above-mentioned SAR satellite on-orbit data processing flow, and determines the basic information in the L1-level SAR imaging product. Then, based on the basic information in the L1-level SAR image product, the Universal Transverse Mercator (UTM) projection range corresponding to the current image is calculated.

[0041] Among them, the basic data of the acquired L1-level SAR image products contain various types of data, such as, but not limited to: the latitude and longitude range of the target area of ​​satellite image acquisition, the center point resolution in the range and azimuth directions, the number of sampling points in the range direction before geometric correction, the near-point observation incident angle, the far-point observation incident angle, the data width, the satellite orbit altitude, and the average radius of the earth, etc.

[0042] In order to more clearly illustrate the specific implementation process of calculating the UTM projection range of the present application, a calculation method proposed in an embodiment of the present application is exemplified below. Figure 2 A flow chart of a method for calculating a UTM projection range proposed in an embodiment of the present application, such as Figure 2 As shown, the method comprises the following steps:

[0043] Step S201, extracting the longitude and latitude coordinates of the four corner points of the observation area from the basic data, and calculating the longitude and latitude range of the minimum circumscribed rectangle of the observation area based on the longitude and latitude coordinates of the four corner points.

[0044] Specifically, the coordinate ranges of the longitude and latitude of the four corner points of the observation area corresponding to the current image are extracted from the basic data of the L1 SAR image, where the longitude and latitude coordinate values ​​of the four corner points: top left (TL), bottom left (BL), top right (TR), and bottom right (BR) are respectively recorded as (TL lon ,TL lat ),(BL lon ,BL lat ),(TR lon ,TR lat ),(BR lon ,BR lat ).

[0045] Then, the longitude and latitude coordinates of the four corner points are used as input data, and the longitude and latitude range of the minimum bounding rectangle of the current SAR image observation area is calculated based on the longitude and latitude coordinates of the four corner points (TL, BL, TR, BR). As a possible implementation method, the longitude and latitude range can be calculated by the following formula:

[0046] (llon ,l lat )=min(TL lon ,BL lon ,TR lon ,BR lon )

[0047] (R lon ,R lat )=max(TL lon ,BL lon ,TR lon ,BR lon )

[0048] (T lon ,T lat )=max(TL lat ,BL lat ,TR lat ,BR lat )

[0049] (B lon ,B lat )=min(TL lat ,BL lat ,TR lat ,BR lat )

[0050] Among them, (L lon ,L lat ), (R lon ,R lat ), (T lon ,T lat ), (B lon ,B lat ) represent the longitude and latitude of the left, right, upper and lower boundary points of the minimum circumscribed rectangle of the current observation area, min(·) represents the minimum value of the current longitude / longitude point set and its corresponding longitude / longitude value, and max(·) represents the maximum value of the current longitude / longitude point set and its corresponding longitude / latitude value.

[0051] Step S202: Map the latitude and longitude range of the minimum bounding rectangle to the UTM projection coordinates of the four boundary points of the minimum bounding rectangle.

[0052] Specifically, the latitude and longitude range of the minimum circumscribed rectangle of the observation area of ​​the current SAR image obtained in the previous step is mapped to UTM projection coordinates. As a possible implementation method, the UTM projection coordinates of the four boundary points can be calculated by the following formula:

[0053] (L utmx ,L utmy )=UTM(L lon ,L lat)

[0054] (R utmx ,R utmy )=UTM(R lon ,R lat )

[0055] (T utmx ,T utmy )=UTM(T lon ,T lat )

[0056] (B utmx ,B utmy )=UTM(B lon ,B kat )

[0057] Among them, UTM(·) means converting longitude and latitude coordinates into UTM projection coordinates, (L utmx ,L utmy ), (R utmx ,R utmy ), (T utmx ,T utmy ), (B utmx ,B utmy ) represent the UTM projection coordinate values ​​of the boundary points (left, right, top, and bottom) of the minimum circumscribed rectangle of the current image observation area.

[0058] Step S203, based on the UTM projection coordinates of the four boundary points, respectively calculate the actual ground lengths in the distance direction and the azimuth direction.

[0059] Specifically, based on the four boundary points (L utmx ,L utmy ), (R utmx ,R utmy ), (T utmx ,T utmy ), (B utmx ,B utmy )’s UTM projection coordinates, and calculate the actual ground lengths in the range direction (RD) and azimuth direction (AD) respectively.

[0060] As a possible implementation method, the calculation can be performed in different ways based on the difference in the angle between the satellite orbit corresponding to the image and the prime meridian. In the first example, when the angle between the satellite orbit corresponding to the current image and the prime meridian is positive, the actual ground length in the range direction is calculated using the following formula:

[0061]

[0062] The actual ground length in azimuth is calculated using the following formula:

[0063]

[0064] In the second example, when the angle between the satellite orbit corresponding to the current image and the prime meridian is negative, the actual ground length in the distance direction is calculated using the following formula:

[0065]

[0066] The actual ground length in azimuth is calculated using the following formula:

[0067]

[0068] Therefore, the calculation result of the UTM projection range output in this step includes the coordinates of the four boundary points of the UTM: utmx ,L utmy ), (R utmx ,R utmy ), (T utmx ,T utmy ), (B utmx ,B utmy ), as well as the actual distance length RD and the actual azimuth length AD.

[0069] Step S102, calculating the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle according to the UTM projection range, and dividing the observation beam range into equal parts based on basic data, and calculating the scaling ratio corresponding to each equal-division interval.

[0070] Specifically, the rotation angle of the minimum circumscribed rectangle of the observation area of ​​the current L1-level SAR image is calculated first. The input data of this process are the coordinates of the four UTM boundary points and the distance length RD calculated in step S101.

[0071] In one embodiment of the present application, the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle is calculated, including: determining the angle between the satellite orbit corresponding to the L1-level SAR image and the prime meridian; when the angle is positive, calculating a first difference based on the UTM projection coordinates of the left boundary point and the UTM projection coordinates of the lower boundary point, and calculating the rotation angle according to the first difference and the actual ground length in the distance direction; when the angle is negative, calculating a second difference based on the UTM projection coordinates of the lower boundary point and the UTM projection coordinates of the right boundary point, and calculating the rotation angle according to the second difference and the actual ground length in the distance direction.

[0072] Specifically, this embodiment calculates the rotation angle of the actual observation area in the minimum circumscribed rectangle based on four boundary points and the distance direction. In the first case, when the angle between the satellite orbit corresponding to the current image and the prime meridian is positive, the rotation angle is calculated by the following formula:

[0073] RA=sin -1 ((L utmy -B utmy ) / RD)

[0074] In the second case, when the angle between the satellite orbit corresponding to the current image and the prime meridian is negative, the rotation angle is calculated by the following formula:

[0075] RA=sin -1 ((B utmy -R utmy ) / RD)

[0076] Among them, RA represents the rotation angle of the actual observation area in the minimum circumscribed rectangle.

[0077] Furthermore, the effective observation beam range corresponding to the observation area of ​​the L1-level SAR image is divided into N equal intervals, and the scaling ratio of the observation pixel corresponding to each divided interval obtained after the equal division compared with the central observation point is calculated, where N is the observation beam angle equal division coefficient.

[0078] It can be understood that the present application divides the observation beam range into N equal parts, which is equivalent to dividing the angle between the near-point observation line and the far-point observation line in the SAR satellite earth observation model into equal parts. Relevant data can be extracted from the basic data and substituted into the SAR satellite earth observation model for calculation.

[0079] In order to more clearly illustrate the specific implementation process of calculating the UTM projection range of the present application, a method proposed in an embodiment of the present application is exemplified below. Figure 3 A flowchart of a method for equally dividing a beam range and calculating a scaling ratio proposed in an embodiment of the present application, such as Figure 3 As shown, the method comprises the following steps:

[0080] Step S301, calculate the observation beam range according to the near point observation incident angle and the far point observation incident angle extracted from the basic data, and divide the observation beam range into equal parts according to a preset equal division value.

[0081] The input data in the equal division and calculation method of this embodiment includes: the near point incident angle of the current image observation extracted from the basic data and the far point angle of incidence And the preset observation beam angle division coefficient N.

[0082] Specifically, based on the near point incident angle and the far point angle of incidence Calculate observation beam range Divide it into N equal parts to get the set of zoomed areas

[0083] As a possible implementation method, the observation beam range is first calculated by the following formula:

[0084]

[0085] The observation beam range is then divided equally using the following formula:

[0086]

[0087] Among them, k = 1, 2,…, N.

[0088] It should be noted that the elements in the zoom area set are obtained for each equally divided interval. The present application can flexibly adjust the geometric correction accuracy by adjusting the value of the observation beam angle equally divided coefficient N in combination with actual factors such as correction accuracy requirements and computing resources. Among them, the larger the value of N (not exceeding the number of sampling points in the actual distance), the higher the geometric correction accuracy, and the greater the calculation complexity. Conversely, the smaller the value of N (not less than 100), the lower the geometric correction accuracy, and the greater the calculation complexity.

[0089] Step S302: Calculate the ground distance set corresponding to each equally divided interval by combining the SAR satellite earth observation model and basic data.

[0090] Among them, the SAR satellite earth observation model proposed in this application is as follows Figure 4 As shown, the observed near point A corresponds to the near point incident angle, and the observed near point B corresponds to the far point incident angle.

[0091] Specifically, based on Figure 4 The SAR satellite earth observation model shown in the figure uses the relevant basic data to substitute into the calculation formula, and can calculate the actual ground distance corresponding to each interval according to the relative near-point incidence angle and the relative far-point incidence angle in each zoom interval, and then obtain the ground distance set {D1, D2, D3, ..., D N}.

[0092] In this embodiment, the SAR satellite earth observation model and basic data are combined to calculate the ground distance set corresponding to each equally divided interval, including: calculating the perigee electromagnetic wave transmission distance and the apogee electromagnetic wave transmission distance in the SAR satellite earth observation model based on the basic data and the angles between each equally divided interval; calculating the ground distance corresponding to the current equally divided interval according to the perigee electromagnetic wave transmission distance, the apogee electromagnetic wave transmission distance, the angle between the current equally divided interval, and the angle of the previous interval adjacent to the current equally divided interval.

[0093] Specifically, for the kth equally divided interval, the actual corresponding ground distance of the interval can be calculated by the following formula:

[0094]

[0095] Where k = 1, 2, ..., N. The various parameters in the above formula can be calculated by the following formula:

[0096]

[0097] Among them, SA means Figure 4 The mid-perigee electromagnetic wave transmission distance, SB represents Figure 4 The distance of electromagnetic wave transmission at mid-distance and far-distance, expressed by ∠SOA Figure 1 The midpoint angle, ∠SOB, represents Figure 1 The midpoint and apocenter angle, RE represents the average radius of the earth, SH represents Figure 4 Medium satellite orbit altitude. represents the angle between the current equal parts, represents the angle between the previous equal parts,

[0098] It is understandable that for the zoom area set For each equally divided interval in , the corresponding actual ground distance can be calculated by the above method, thus obtaining the ground distance set {D1, D2, D3, …, D N}.

[0099] Step S303, based on the ground distance set, calculate the ratio set of the distance to the total pixels corresponding to each equally divided interval.

[0100] Specifically, based on The corresponding ground distance set {D1,D2,D3,…,D N}, calculate the ratio set of distance to total pixels before geometric correction corresponding to each equal interval {R1, R2, R3, …, R N}.

[0101] In this embodiment, based on the ground distance set, the ratio set of the total distance pixels corresponding to each equally divided interval is calculated, including: calculating the sum of all ground distances in the ground distance set; dividing the ground distance corresponding to each equally divided interval by the sum to obtain the total distance pixels corresponding to each equally divided interval.

[0102] Specifically, the total distance pixels corresponding to the kth equally divided interval can be calculated by the following formula:

[0103]

[0104] Among them, k = 1, 2,…, N.

[0105] Step S304, calculating the ground distance of the center interval between each equally divided area, and calculating the zoom ratio set corresponding to each equally divided area according to the ground distance set and the ground distance of the center interval.

[0106] Specifically, based on the ground distance between the center intervals of each equal partition, the actual corresponding scaling ratio set {Z1, Z2, Z3, ..., Z N}.

[0107] In this embodiment, the set of scaling ratios corresponding to each equally divided interval is calculated, including: dividing the ground distance corresponding to each equally divided interval by the ground distance of the central interval to obtain the scaling ratio corresponding to each equally divided interval.

[0108] Specifically, the scaling ratio corresponding to the kth equally divided interval can be calculated by the following formula:

[0109]

[0110] Among them, D mid Indicates the ground distance corresponding to the center interval in the set of equally divided intervals. The center interval can be the interval at the most central position among multiple equally divided intervals.

[0111] Therefore, the output of this step includes: the ground distance set {D1, D2, D3, ..., D N}, the ratio set of distance to total pixels {R1,R2,R3,…,R N} and the scaling ratio set {Z1,Z2,Z3,…,Z N}.

[0112] Step S103, low-frequency resampling is performed on the pixel points between each equal partition according to the corresponding scaling ratio, and the resampled image is rotated according to the rotation angle to obtain a geometrically corrected L2-level SAR image product.

[0113] Specifically, the pixels between each equal partition are low-frequency resampled according to the corresponding scaling ratio calculated in step S102 to achieve fast geometric correction. The resampling may include operations such as stretching or compressing the pixels.

[0114] The input data of this step includes the L1-level SAR image product to be corrected and the ground distance set {D1, D2, D3, ..., D N} and the ratio set of distance to total pixels {R1,R2,R3,…,R N}.

[0115] In one embodiment of the present application, the pixel points between each equally divided area are low-frequency resampled according to the corresponding scaling ratio, including: calculating the pixel point set of the equally divided interval based on the ratio set of the total pixels in the distance direction, the number of points in the azimuth direction between each equally divided area, and the number of pixel points in the distance direction between each equally divided area; combining the pixel point set of the equally divided interval and the scaling ratio set, stretching the pixels in the near point interval and compressing the pixels in the far point interval.

[0116] Specifically, in this embodiment, when performing low-frequency resampling fast geometric correction, in the first step, based on the ratio set {R1, R2, R3, ..., R N}, select the matrix data between each equal partition for distance correction, and obtain the pixel point set of the equal partition interval before geometric correction as {RM1,RM2,RM3,…,RM N}.

[0117] in, P represents the number of azimuth points between the kth equal partitions (this value is fixed), Q k Indicates the distance between the equal partitions in pixels. is a parameter used to represent the dimension, Q k It can be calculated by the following formula:

[0118] Q k =R k ·RDN

[0119] Among them, RDN represents the total number of pixels in the range direction of the L1-level SAR image before geometric correction.

[0120] The second step is to divide the pixel points into equal intervals {RM1, RM2, RM3, …, RM N} and the scaling ratio set {Z1,Z2,Z3,…,Z N}, stretch the pixels in the near point interval and compress the pixels in the far point interval. The specific implementation process of stretching or compressing the relevant pixels is as follows: Figure 5 As shown. Thus, the geometrically corrected image can be obtained.

[0121] Furthermore, after obtaining the geometrically corrected image, the method further includes: a third step of vertically mirror-flipping the corrected image according to the data storage rule, and then rotating the flipped image by the rotation angle RA calculated in step S102, that is, rotating the image by the rotation angle of the actual observation area in the minimum circumscribed rectangle, and obtaining I sAR .

[0122] Finally, write the longitude and latitude ranges corresponding to the four corner points of the observation area boundary into I SAR The final L2 SAR image product IG with coordinate information is obtained from the data SAR .

[0123] In order to more clearly and intuitively describe the advantages of the satellite image geometric correction method based on low-frequency resampling of the present application compared with other correction methods, the following is an example of geometric correction of specific satellite image data in practical applications.

[0124] In one embodiment of the present application, taking the L1 level data of the GF3 satellite as an example, the two types of data of the satellite are geometrically corrected by three methods: geometric correction based on RPC in the related art, geometric correction based on RPC and DEM, and geometric correction based on the low-frequency resampling method proposed in the present application. The geometric correction results of the two types of data are shown as follows: Figure 6 and Figure 7 shown.

[0125] Specifically, the test is completed by inputting L1A-level data of 24354*24518 pixels (the resolution of the data is 5m, the width is 50km, and dual polarization), and the test results are shown in Table 1.

[0126] Table 1 Comparative analysis of geometric correction methods

[0127]

[0128] From Table 1 and Figure 6 and Figure 7 It can be seen that the correction method of the present application can greatly reduce the time required for geometric correction of SAR images without significantly reducing the geometric correction accuracy and meet the accuracy requirements of SAR image geometric correction. It can achieve rapid correction with low computational complexity and improve the efficiency of geometric correction of SAR images.

[0129] In summary, the satellite image geometric correction method based on low-frequency resampling in the embodiment of the present application first calculates the UTM projection range corresponding to the current image based on the basic information in the L1-level SAR image product; then calculates the rotation angle of the observation area in the minimum circumscribed rectangle according to the UTM projection range; then divides the interval of the effective observation beam range into N equal parts based on the basic information, and calculates the scaling ratio of the observation pixels corresponding to each equal-divided interval compared to the central observation point; finally, resamples the pixel points between each equal-divided interval according to the scaling ratio, rotates the resampled image according to the rotation angle, and obtains the corrected output result. Therefore, this method can avoid steps such as calculating RPC parameters, quickly implement SAR image geometric correction, and greatly reduce the time required for SAR image geometric correction on the basis of ensuring the accuracy of SAR image geometric correction, reduce the computational complexity of geometric correction and the consumed computing resources, and can be applied to the on-orbit data processing tasks of SAR satellite AI real-time processors, thereby improving the convenience and applicability of satellite image geometric correction.

[0130] In order to implement the above embodiment, the present application also proposes a satellite image geometric correction system based on low-frequency resampling. Figure 8 A schematic diagram of the structure of a satellite image geometric correction system based on low-frequency resampling proposed in an embodiment of the present application is shown in FIG. Figure 8 As shown, the system includes: a calculation module 100, an equal division module 200 and a correction module 300.

[0131] The calculation module 100 is used to obtain the basic data of the existing L1-level synthetic aperture radar SAR image product, and calculate the Universal Transverse Mercator UTM projection range corresponding to the L1-level SAR image according to the basic data.

[0132] The equal division module 200 is used to calculate the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle according to the UTM projection range, and to divide the observation beam range into equal parts based on basic data, and calculate the scaling ratio corresponding to each equal division interval.

[0133] The correction module 300 is used to perform low-frequency resampling of the pixel points between each equal partition according to the corresponding scaling ratio, and rotate the resampled image according to the rotation angle to obtain a geometrically corrected L2-level SAR image product.

[0134] Optionally, in one embodiment of the present application, the calculation module 100 is specifically used to: extract the longitude and latitude coordinates of the four corner points of the observation area from the basic data, and calculate the longitude and latitude range of the minimum enclosing rectangle of the observation area based on the longitude and latitude coordinates of the four corner points; map the longitude and latitude range of the minimum enclosing rectangle to the UTM projection coordinates of the four boundary points of the minimum enclosing rectangle; and calculate the actual ground length in the distance direction and the azimuth direction based on the UTM projection coordinates of the four boundary points.

[0135] Optionally, in one embodiment of the present application, the equal division module 200 is specifically used to: determine the angle between the satellite orbit corresponding to the L1-level SAR image and the prime meridian; when the angle is positive, calculate a first difference based on the UTM projection coordinates of the left boundary point and the UTM projection coordinates of the lower boundary point, and calculate the rotation angle based on the first difference and the actual ground length in the distance direction; when the angle is negative, calculate a second difference based on the UTM projection coordinates of the lower boundary point and the UTM projection coordinates of the right boundary point, and calculate the rotation angle based on the second difference and the actual ground length in the distance direction.

[0136] Optionally, in one embodiment of the present application, the equal division module 200 is specifically used to: calculate the observation beam range based on the near-point observation incident angle and the far-point observation incident angle extracted from the basic data, and divide the observation beam range into equal parts according to a preset equal division value; calculate the ground distance set corresponding to each equal division interval in combination with the SAR satellite earth observation model and the basic data; calculate the ratio set of the distance to the total pixel corresponding to each equal division interval based on the ground distance set; calculate the ground distance of the central interval between each equal division area, and calculate the scaling ratio set corresponding to each equal division interval based on the ground distance set and the ground distance of the central interval.

[0137] Optionally, in one embodiment of the present application, the equal division module 200 is specifically used to: calculate the perigee electromagnetic wave transmission distance and the apogee electromagnetic wave transmission distance in the SAR satellite earth observation model based on basic data and the angles between each equal division; calculate the ground distance corresponding to the current equal division interval according to the perigee electromagnetic wave transmission distance, the apogee electromagnetic wave transmission distance, the angle between the current equal divisions, and the angle of the previous interval adjacent to the current equal division interval.

[0138] Optionally, in one embodiment of the present application, the equal division module 200 is specifically used to: calculate the sum of all ground distances in the ground distance set; divide the ground distance corresponding to each equal division interval by the sum to obtain the total distance pixels corresponding to each equal division interval; divide the ground distance corresponding to each equal division interval by the ground distance of the central interval to obtain the zoom ratio corresponding to each equal division interval.

[0139] Optionally, in one embodiment of the present application, the correction module 300 is specifically used to: calculate a set of pixel points in the equally divided interval based on a set of ratios of the total pixels in the distance direction, the number of points in the azimuth direction between each equally divided area, and the number of pixel points in the distance direction between each equally divided area; and stretch the pixels in the near point interval and compress the pixels in the far point interval in combination with the set of pixel points in the equally divided interval and the set of scaling ratios.

[0140] It should be noted that the above explanation of the embodiment of the satellite image geometric correction method based on low-frequency resampling is also applicable to the system of this embodiment, and will not be repeated here.

[0141] To sum up, the satellite image geometric correction system based on low-frequency resampling of the embodiment of the present application can avoid steps such as calculating RPC parameters and quickly implement SAR image geometric correction. On the basis of ensuring the accuracy of SAR image geometric correction, it greatly reduces the time required for SAR image geometric correction, reduces the computational complexity of geometric correction and the consumed computing resources, can be applied to the on-orbit data processing tasks of SAR satellite AI real-time processors, and improves the convenience and applicability of satellite image geometric correction.

[0142] In order to implement the above embodiments, the present application also proposes a non-temporary computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the satellite image geometric correction method based on low-frequency resampling as described in any of the above embodiments is implemented.

[0143] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present application. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, without contradiction.

[0144] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include at least one of the features. In the description of this application, the meaning of "plurality" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined.

[0145] Any process or method description in a flowchart or otherwise described herein may be understood to represent a module, fragment or portion of code comprising one or more executable instructions for implementing the steps of a custom logical function or process, and the scope of the preferred embodiments of the present application includes alternative implementations in which functions may not be performed in the order shown or discussed, including performing functions in a substantially simultaneous manner or in the reverse order depending on the functions involved, which should be understood by technicians in the technical field to which the embodiments of the present application belong.

[0146] The logic and / or steps represented in the flowchart or otherwise described herein, for example, can be considered as an ordered list of executable instructions for implementing logical functions, and can be embodied in any computer-readable medium for use by an instruction execution system, device or apparatus (such as a computer-based system, a system including a processor, or other system that can fetch instructions from an instruction execution system, device or apparatus and execute the instructions), or in combination with these instruction execution systems, devices or apparatuses. For the purpose of this specification, "computer-readable medium" can be any device that can contain, store, communicate, propagate or transmit a program for use by an instruction execution system, device or apparatus, or in combination with these instruction execution systems, devices or apparatuses. More specific examples of computer-readable media (a non-exhaustive list) include the following: an electrical connection with one or more wires (electronic device), a portable computer disk box (magnetic device), a random access memory (RAM), a read-only memory (ROM), an erasable and programmable read-only memory (EPROM or flash memory), a fiber optic device, and a portable compact disk read-only memory (CDROM). In addition, the computer-readable medium may even be paper or other suitable medium on which the program is printed, since the program may be obtained electronically, for example, by optically scanning the paper or other medium and then editing, interpreting or processing in other suitable ways if necessary, and then stored in a computer memory.

[0147] It should be understood that the various parts of the present application can be implemented by hardware, software, firmware or a combination thereof. In the above-mentioned embodiments, multiple steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented by hardware, as in another embodiment, it can be implemented by any one of the following technologies known in the art or their combination: a discrete logic circuit having a logic gate circuit for implementing a logic function for a data signal, a dedicated integrated circuit having a suitable combination of logic gate circuits, a programmable gate array (PGA), a field programmable gate array (FPGA), etc.

[0148] A person skilled in the art may understand that all or part of the steps in the method for implementing the above-mentioned embodiment may be completed by instructing related hardware through a program, and the program may be stored in a computer-readable storage medium, which, when executed, includes one or a combination of the steps of the method embodiment.

[0149] In addition, each functional unit in each embodiment of the present application may be integrated into a processing module, or each unit may exist physically separately, or two or more units may be integrated into one module. The above-mentioned integrated module may be implemented in the form of hardware or in the form of a software functional module. If the integrated module is implemented in the form of a software functional module and sold or used as an independent product, it may also be stored in a computer-readable storage medium.

[0150] The storage medium mentioned above may be a read-only memory, a magnetic disk or an optical disk, etc. Although the embodiments of the present application have been shown and described above, it can be understood that the above embodiments are exemplary and cannot be understood as limiting the present application. A person of ordinary skill in the art may change, modify, replace and modify the above embodiments within the scope of the present application.

Claims

1. A satellite image geometric correction method based on low-frequency resampling, characterized in that: The following steps are involved: Obtaining basic data of an existing L1-level synthetic aperture radar SAR image product, and calculating a Universal Transverse Mercator (UTM) projection range corresponding to the L1-level SAR image according to the basic data; Calculating the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle according to the UTM projection range, and dividing the observation beam range into equal parts based on the basic data, and calculating the scaling ratio corresponding to each equal-divided interval; The pixel points between each of the equal partitions are low-frequency resampled according to the corresponding scaling ratio, and the resampled image is rotated according to the rotation angle to obtain a geometrically corrected L2-level SAR image product.

2. The method according to claim 1, characterized in that The calculating the universal transverse Mercator UTM projection range corresponding to the L1-level SAR image according to the basic data includes: Extracting the longitude and latitude coordinates of the four corner points of the observation area from the basic data, and calculating the longitude and latitude range of the minimum circumscribed rectangle of the observation area based on the longitude and latitude coordinates of the four corner points; Mapping the latitude and longitude range of the minimum bounding rectangle into UTM projection coordinates of the four boundary points of the minimum bounding rectangle; Based on the UTM projection coordinates of the four boundary points, the actual ground lengths in the distance direction and the azimuth direction are calculated respectively.

3. The method according to claim 2, characterized in that The step of calculating the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle includes: Determine the angle between the satellite orbit corresponding to the L1-level SAR image and the prime meridian; When the angle is positive, a first difference is calculated based on the UTM projection coordinates of the left boundary point and the UTM projection coordinates of the lower boundary point, and the rotation angle is calculated according to the first difference and the actual ground length in the distance direction; When the angle is a negative value, a second difference is calculated based on the UTM projection coordinates of the lower boundary point and the UTM projection coordinates of the right boundary point, and the rotation angle is calculated according to the second difference and the actual ground length in the distance direction.

4. The method according to claim 1, characterized in that The step of equally dividing the observation beam range based on the basic data and calculating the scaling ratio corresponding to each equally divided interval includes: Calculate the observation beam range according to the near-point observation incident angle and the far-point observation incident angle extracted from the basic data, and divide the observation beam range into equal parts according to a preset equal division value; Combining the SAR satellite earth observation model and the basic data, calculating the ground distance set corresponding to each equally divided interval; According to the ground distance set, calculate the ratio set of distances corresponding to the equal-divided intervals to the total pixels; The ground distances of the central intervals between the equal divisions are calculated, and the set of zoom ratios corresponding to the equal divisions is calculated according to the ground distance set and the ground distances of the central intervals.

5. The method according to claim 4, characterized in that The combining of the SAR satellite earth observation model and the basic data to calculate the ground distance set corresponding to each equally divided interval includes: Based on the basic data and the angles between the equal partitions, calculating the perigee electromagnetic wave transmission distance and the apogee electromagnetic wave transmission distance in the SAR satellite earth observation model; The ground distance corresponding to the current equally divided interval is calculated based on the perigee electromagnetic wave transmission distance, the apogee electromagnetic wave transmission distance, the angle between the current equally divided intervals, and the angle of the previous interval adjacent to the current equally divided interval.

6. The method according to claim 4, characterized in that The calculating, based on the ground distance set, a ratio set of distances corresponding to the respective equally divided intervals to total pixels comprises: Calculating the sum of all ground distances in the ground distance set; Divide the ground distance corresponding to each of the equally divided intervals by the sum value to obtain the total distance pixels corresponding to each of the equally divided intervals; The calculating of the set of scaling ratios corresponding to each equally divided interval includes: The ground distance corresponding to each of the equally divided intervals is divided by the ground distance of the central interval to obtain the zoom ratio corresponding to each of the equally divided intervals.

7. The method according to claim 4, characterized in that The step of performing low-frequency resampling of the pixel points between each of the equal partitions according to the corresponding scaling ratio comprises: Calculate a pixel point set between the equal partitions based on the ratio set of the total pixels in the distance direction, the number of points in the azimuth direction between each equal partition, and the number of pixel points in the distance direction between each equal partition; In combination with the equally divided interval pixel point set and the scaling ratio set, the near point interval pixels are stretched and the far point interval pixels are compressed.

8. A satellite image geometric correction system based on low-frequency resampling, characterized in that: include: A calculation module is used to obtain basic data of existing L1-level synthetic aperture radar SAR image products, and calculate the universal transverse Mercator UTM projection range corresponding to the L1-level SAR image according to the basic data; An equal division module, used for calculating the rotation angle of the observation area of ​​the L1-level SAR image in the minimum circumscribed rectangle according to the UTM projection range, and equally dividing the observation beam range based on the basic data, and calculating the scaling ratio corresponding to each equally divided interval; The correction module is used to perform low-frequency resampling of the pixel points between each of the equal partitions according to the corresponding scaling ratio, and rotate the resampled image according to the rotation angle to obtain a geometrically corrected L2-level SAR image product.

9. The system according to claim 8, characterized in that The computing module is specifically used for: Extracting the longitude and latitude coordinates of the four corner points of the observation area from the basic data, and calculating the longitude and latitude range of the minimum circumscribed rectangle of the observation area based on the longitude and latitude coordinates of the four corner points; Mapping the latitude and longitude range of the minimum bounding rectangle into UTM projection coordinates of the four boundary points of the minimum bounding rectangle; The actual ground lengths in the distance direction and the azimuth direction are calculated respectively based on the UTM projection coordinates of the four boundary points.

10. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the satellite image geometric correction method based on low-frequency resampling as described in any one of claims 1 to 7 is implemented.

Images