Satellite-borne SAR image in-orbit radiation and geometric correction method based on parallel computing architecture

Abstract

The invention relates to the technical field of synthetic aperture radar signal processing, and discloses a satellite-borne SAR image in-orbit radiation and geometric correction method based on a parallel computing architecture, and the radiation correction method comprises the steps: obtaining ephemeris parameters of a satellite under a geocentric coordinate system at all azimuth moments; calculating a slant distance between a satellite beam center pointing scene center point and an antenna phase center at each azimuth moment and a slant distance between each sampling point and a satellite in an SAR echo distance direction at each azimuth moment; calculating satellite beam center pointing at each azimuth moment and a lower visual angle between the satellite and a sub-satellite point connecting line; calculating an antenna distance direction beam directional diagram and constructing a radiation correction function; and performing parallel correction on radiation distortion of the SAR image based on the radiation correction function. According to the method, the SAR image quality and the accurate positioning of the SAR image are effectively improved, the radiation fidelity of the SAR image is met, meanwhile, the efficiency of a data processing task is improved, and the task requirement for high real-time performance is considered.

Description

Technical Field

[0001] The present invention relates to the technical field of synthetic aperture radar signal processing, and in particular to an on-orbit radiation and geometric correction method for spaceborne SAR images based on a parallel computing architecture. Background Art

[0002] Synthetic Aperture Radar (SAR) is an active high-resolution microwave imaging system that can continuously observe the area of ​​interest over a large area. In recent years, SAR has developed towards multi-polarization and multi-imaging modes, playing an important role in ecological environmental surveys, military reconnaissance, and civilian fields. This requires that the imaging results of SAR images should not only have high resolution, but also have high radiometric fidelity and geometric fidelity to further improve the quality of SAR images to better meet the needs of a wide range of applications. The popularity of high-quality and high-resolution SAR images has also promoted the real-time demand for rapid extraction of remote sensing image information. However, the response time of traditional remote sensing systems based on satellite data transmitted to ground stations for processing is too long, making it difficult to meet the needs of various applications with high real-time requirements. Radiometric correction and geometric correction are two key steps in the on-orbit preprocessing of satellite-borne SAR images. Radiometric correction can correct the radiometric distortion of the entire SAR image scene and improve the quality of SAR images; geometric correction can achieve the correct association and effective recovery of geographic coordinates and original pixel information of the image, and complete the precise positioning of target information in SAR images. Therefore, how to quickly perform radiation correction and geometric correction on the entire SAR image scene without radiation calibration data and digital elevation model is an application demand that needs to be urgently solved to improve the real-time performance of on-orbit SAR image preprocessing, which has important application value to the development of the national economy and national defense. Summary of the invention

[0003] The technical problem solved by the present invention is to provide an on-orbit radiation and geometric correction method for spaceborne SAR images based on a parallel computing architecture, which effectively solves the technical problems of low SAR image quality and slow processing efficiency.

[0004] The above technical problems are solved by the following technical solutions:

[0005] A method for on-orbit radiation correction of spaceborne SAR images based on a parallel computing architecture, the method comprising:

[0006] Acquire the ephemeris parameters of the satellite in the geocentric coordinate system at each orientation time, wherein the ephemeris parameters include at least one of a position vector, a velocity vector, a three-axis attitude angle, a pulse repetition frequency, a wave position sampling start time, a signal pulse width, a signal wavelength, an antenna range aperture length, an antenna azimuth aperture length, an antenna beam pitch angle, an antenna beam azimuth angle, a range sampling point number, an azimuth sampling point number, and a signal sampling frequency;

[0007] Calculate the slant distance between the satellite beam center pointing to the scene center and the antenna phase center at each azimuth moment, as well as the slant distance between each sampling point and the satellite at each azimuth moment;

[0008] Based on the slant distance between the center of the satellite beam pointing to the scene center point and the antenna phase center at each azimuth moment, the slant distance between each sampling point and the satellite at each azimuth moment, and the ephemeris parameters, the downward viewing angle between the direction of the satellite beam center at each azimuth moment and the line connecting the satellite and the sub-satellite point is calculated;

[0009] Calculate the antenna range beam pattern based on the satellite beam center direction at each orientation and the downward viewing angle between the satellite and the sub-satellite point, and construct a radiation correction function;

[0010] The radiation distortion of the SAR image is corrected in parallel based on the radiation correction function.

[0011] Beneficial effects: By obtaining the ephemeris parameters of the satellite in the geocentric coordinate system at each azimuth time, the slant distance between the center of the satellite beam pointing to the scene center and the antenna phase center at each azimuth time and the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth time are calculated; further, the downward viewing angle between the satellite beam center pointing at each azimuth time and the line connecting the satellite and the sub-satellite point is calculated based on the slant distance and ephemeris parameters; then, the antenna distance beam pattern is calculated based on the downward viewing angle between the satellite beam center pointing at each azimuth time and the line connecting the satellite and the sub-satellite point, and a radiation correction function is constructed. Finally, the radiation distortion of the SAR image is corrected in parallel through the radiation correction function. The quality of SAR images is effectively improved, the radiation fidelity of SAR images is met, and the efficiency of data processing tasks is also improved, while taking into account the task requirements of high real-time performance.

[0012] In one of the embodiments, the parallel correction of the radiation distortion of the SAR image based on the radiation correction function includes:

[0013] Dividing the SAR image into a plurality of sub-SAR images;

[0014] Setting thread blocks corresponding to the plurality of sub-SAR images based on a CUDA programming model;

[0015] Based on the graphics processor, the plurality of thread blocks are operated in parallel to complete the parallel correction of the radiation distortion of the plurality of sub-SAR images.

[0016] In one embodiment, the three-axis attitude angle includes a ground roll angle, a ground pitch angle, and a ground yaw angle; the calculation based on the ephemeris parameters of the slant distance between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth moment and the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth moment includes:

[0017] The number of sub-satellite point echo pulses and the SAR echo acquisition start time are calculated based on the three-axis attitude angle and the velocity vector and position vector of the satellite in the geocentric coordinate system. The calculation formula is shown in formula (1):

[0018]

[0019] Based on the echo acquisition start time and the number of echo pulses at the sub-satellite point, the slant distance between the satellite at the beam center and the beam pointing point and the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth time are calculated. The calculation formula is shown in formula (2):

[0020]

[0021] Where c is the speed of light; R j R is the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth moment; ref is the slant distance between the satellite beam center pointing to the scene center and the antenna phase center at each orientation moment; center is the start time of echo acquisition; M is the number of sub-satellite echo pulses; t0 is the start time of wave position sampling; PRF is the pulse repetition frequency; T p is the signal pulse width; f s is the signal sampling frequency; N r is the number of sampling points in the distance direction.

[0022] In one embodiment, the calculation of the downward viewing angle between the satellite beam center pointing to the center point of the scene and the antenna phase center at each azimuth moment based on the slant distance between the satellite beam center pointing to the center point of the scene and the antenna phase center at each azimuth moment, the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth moment, and the ephemeris parameters includes:

[0023] Calculate the distance between the satellite and the sub-satellite point based on the position vector of the satellite and the average radius of the earth;

[0024] The downward viewing angle between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth time is calculated based on the slant distance between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth time, the slant distance between each sampling point and the satellite at each azimuth time, the distance between the satellite and the sub-satellite point, and the average radius of the earth; the calculation formula is shown in formula (3):

[0025]

[0026] Among them, β j is the downward angle between the satellite beam center direction and the line connecting the satellite and the sub-satellite point at each azimuth; H is the distance between the satellite and the sub-satellite point; R e is the average radius of the Earth.

[0027] In one of the embodiments, the step of calculating the antenna range beam pattern based on the satellite beam center direction at each orientation time and the downward viewing angle between the line connecting the satellite and the sub-satellite point and constructing a radiation correction function includes:

[0028] Calculate the difference between the downward viewing angle between the satellite beam center pointing and the line connecting the satellite and the sub-satellite point at each orientation and the downward viewing angle corresponding to the scene center. The calculation formula is shown in formula (4):

[0029] Δβ=β j -β ref (4)

[0030] The antenna range beam pattern is calculated based on the difference Δβ, and the calculation formula is shown in formula (5):

[0031]

[0032] A radiation correction function is constructed based on the antenna range beam pattern, and the calculation formula is shown in formula (6):

[0033]

[0034] Among them, β j β is the downward angle between the satellite beam center direction and the line connecting the satellite and the sub-satellite point at each azimuth moment; ref is the lower viewing angle corresponding to the center of the scene; Δβ is the difference between the lower viewing angle corresponding to each sampling point of the SAR echo distance at each azimuth and the lower viewing angle corresponding to the center of the scene; W a (β j ,H) is the antenna distance beam pattern; L r is the antenna distance aperture length; λ is the signal wavelength; F(W a (β j,H)) is the radiation correction function; α is the adjustable scale coefficient after image normalization; n is a constant, n∈[0,4].

[0035] On the other hand, the present invention also provides an on-orbit geometric correction method for spaceborne SAR images based on a parallel computing architecture, wherein the method further comprises:

[0036] Calculate the coordinates of the ground scene point illuminated by the antenna beam and the Doppler center frequency of the satellite relative to the aiming point of the antenna beam based on the position vector, velocity vector and antenna parameters of the satellite in the geocentric coordinate system of the satellite;

[0037] Obtaining average elevation information of the target scene, and combining the Doppler center frequency in a range Doppler SAR positioning model to generate coordinates and longitude and latitude information in a geocentric coordinate system corresponding to each pixel point in the SAR image;

[0038] Establishing an index number mapping relationship between the coordinates in the geocentric coordinate system corresponding to each pixel point in the SAR image and the geographic longitude and latitude grid;

[0039] A kernel function is constructed based on the graphics processor to map the SAR image pixels in the imaging plane to the latitude and longitude grids in sequence to complete the geometric correction.

[0040] In one embodiment, the calculating the coordinates of the ground scene point illuminated by the antenna beam and the Doppler center frequency of the satellite relative to the aiming point of the antenna beam based on the position vector, velocity vector and antenna parameters of the satellite in the geocentric coordinate system of the satellite includes:

[0041] Calculating a rotation matrix from an inertial coordinate system to a satellite orbit coordinate system based on the ephemeris parameters of the satellite;

[0042] Based on the satellite's three-axis attitude angle and the satellite's rotation order, calculate the rotation matrix from the satellite's orbital coordinate system to the satellite's body coordinate system;

[0043] Calculate the rotation matrix from the satellite body coordinate system to the antenna coordinate system based on the antenna azimuth angle and antenna elevation angle in the antenna coordinate system;

[0044] Calculate the intersection of the aiming vector of the antenna beam center and the earth's surface in the geocentric coordinate system, obtain the coordinates of the ground scene point illuminated by the antenna beam, and calculate the velocity vector of the antenna beam in the geocentric coordinate system;

[0045] Based on the coordinates of the ground scene point illuminated by the antenna beam and the velocity vector of the antenna beam in the geocentric coordinate system, the Doppler center frequency of the satellite relative to the antenna beam aiming point is calculated.

[0046] In one embodiment, the obtaining of average elevation information of the target scene and combining the Doppler center frequency in a range Doppler SAR positioning model to generate coordinates and longitude and latitude information in a geocentric coordinate system corresponding to each pixel point in the SAR image include:

[0047] Based on the average elevation information of the target scene and the coordinates and longitude and latitude information of each pixel point in the SAR image in the geocentric coordinate system;

[0048] In the range Doppler SAR positioning model, based on the average elevation information and Doppler center frequency of the target scene, the coordinates of the four corner points in the SAR image corresponding to the geocentric coordinate system are calculated, and the coordinates are converted into longitude and latitude information; the longitudes Ω of the four corner points are {Ω left_up ,Ω left_down ,Ω right_up ,Ω right_down}, and the latitudes Φ are {Φ left_up , Φ left_down , Φ right_up , Φ right_down};

[0049] Compare the upper left corner point Ω left_up and the latitude value of the lower left corner Ω left_down , determine whether the scene pixels are consistent with the north latitude to south latitude direction;

[0050] In response to the scene pixel being consistent with the north latitude to south latitude direction, the longitude and latitude values ​​corresponding to all the pixel points of the SAR image are calculated based on the longitude and latitude of the four corner points and by constructing a kernel function.

[0051] On the other hand, the present invention further provides an electronic device, comprising: a memory, a processor, a communication interface and a communication bus, wherein the processor, the memory and the communication interface communicate with each other via the communication bus;

[0052] The memory is used to store at least one executable instruction, and the executable instruction enables the processor to execute the on-orbit radiation and geometric correction method of space-borne SAR images based on the parallel computing architecture as described above.

[0053] On the other hand, the present invention also provides a computer-readable storage medium, which stores computer instructions, and the computer instructions are used to enable a processor to implement the on-orbit radiation and geometric correction method of space-borne SAR images based on a parallel computing architecture as described above when executed. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] In order to more clearly illustrate the specific implementation methods of the present invention or the technical solutions in the prior art, the drawings required for use in the specific implementation methods or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some implementation methods of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0055] Figure 1 A schematic flow chart of an on-orbit radiation correction method for spaceborne SAR images based on a parallel computing architecture according to an embodiment of the present invention;

[0056] Figure 2 It is a flowchart of an on-orbit radiation correction method for spaceborne SAR images based on a parallel computing architecture according to another embodiment of the present invention;

[0057] Figure 3 A schematic flow chart of an on-orbit geometric correction method for spaceborne SAR images based on a parallel computing architecture according to an embodiment of the present invention;

[0058] Figure 4 A schematic flow chart of an on-orbit radiation and geometric correction method for spaceborne SAR images based on a parallel computing architecture according to an embodiment of the present invention;

[0059] Figure 5 It is a SAR image of a certain scene without radiation correction;

[0060] Figure 6 A SAR image that has been subjected to radiation correction by an on-orbit radiation correction method for spaceborne SAR images based on a parallel computing architecture according to an embodiment of the present invention;

[0061] Figure 7 is the original SAR image of a certain scene;

[0062] Figure 8 The SAR image is corrected by an on-orbit radiation and geometric correction method for spaceborne SAR images based on a parallel computing architecture according to an embodiment of the present invention;

[0063] Fig. 9 The figure is a schematic diagram of the structure of an electronic device according to an embodiment of the present invention. DETAILED DESCRIPTION

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

[0065] In the description of the present invention, it should be understood that 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 technical features indicated. Thus, the features defined as "first" and "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, unless otherwise specified, "plurality" means two or more. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0066] According to an embodiment of the present invention, Figure 1 , Figure 2 and Figure 4 As shown, a method for on-orbit radiation correction of spaceborne SAR images based on a parallel computing architecture is provided, comprising the following steps:

[0067] Step S100: Obtain the satellite's ephemeris parameters in the geocentric coordinate system at each position and time. The ephemeris parameters include the position vector P sat =[X sat , Y sat , Z sat ], velocity vector V sat =[Vx sat , Vx sat , Vz sat ], three-axis attitude angle, pulse repetition frequency PRF, wave position sampling start time t0, signal pulse width T p , signal wavelength λ, antenna distance aperture length L r 、Antenna azimuth aperture length L a , antenna beam elevation angle Antenna beam azimuth Number of sampling points in distance direction N r , the number of sampling points in azimuth direction N a Sum signal sampling frequency f s etc.; the three-axis attitude angle includes the ground roll angle θ Roll 、Pitch angle to ground θ Pitch and the yaw angle θ Yaw .

[0068] Step S200: Calculate the slant distance R between the satellite beam center pointing to the scene center and the antenna phase center at each orientation time. ref And the slant distance R between each sampling point and the satellite at each azimuth moment j ,j∈[1:N r ].

[0069] In step S200, specifically, the number of sub-satellite point echo pulses and the SAR echo acquisition start time are calculated based on the three-axis attitude angle and the velocity vector and position vector of the satellite in the geocentric coordinate system. The calculation formula is shown in formula (1):

[0070]

[0071] Based on the echo acquisition start time and the number of echo pulses at the sub-satellite point, the slant distance between the satellite at the beam center and the beam pointing point and the slant distance between each sampling point and the satellite in the SAR echo range at each azimuth time are calculated. The calculation formula is shown in formula (2):

[0072]

[0073] Where c is the speed of light; R is the slant distance between each sampling point and the satellite at each azimuth time; R ref is the slant distance between the satellite beam center pointing to the scene center and the antenna phase center at each orientation moment; center is the start time of echo acquisition; M is the number of sub-satellite echo pulses; t0 is the start time of wave position sampling; PRF is the pulse repetition frequency; T p is the signal pulse width; f s is the signal sampling frequency; N r is the number of sampling points in the distance direction.

[0074] Step S300, based on the slant distance between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth moment, the slant distance between each sampling point and the satellite at each azimuth moment, and the ephemeris parameters, calculate the downward angle β between the satellite beam center pointing to the satellite and the line connecting the satellite and the sub-satellite point at each azimuth moment j , j∈[1:N r ].

[0075] In step S300, specifically including calculating the distance between the satellite and the sub-satellite point based on the position vector of the satellite and the average radius of the earth;

[0076] Based on the slant distance between the center of the satellite beam pointing to the scene and the phase center of the antenna at each azimuth moment, the slant distance between each sampling point and the satellite at each azimuth moment, the distance between the satellite and the sub-satellite point, and the average radius of the earth, the downward viewing angle between the direction of the satellite beam center and the line connecting the satellite and the sub-satellite point at each azimuth moment is calculated; the calculation formula is shown in formula (3):

[0077]

[0078] Where βi is the downward angle between the satellite beam center direction and the line connecting the satellite and the sub-satellite point at each azimuth; H is the distance between the satellite and the sub-satellite point; Re is the average radius of the Earth.

[0079] The calculation formula for calculating the distance H between the satellite and the sub-satellite point based on the satellite's position vector and the average radius of the earth is:

[0080] H = sqrt(X sat 2 +Y sat 2 +Z sat 2 )-R e .

[0081] Step S400: Calculate the antenna distance beam pattern W based on the satellite beam center direction at each orientation and the downward viewing angle between the satellite and the sub-satellite point. a (β, H) and construct the radiation correction function F(W a (β j , H)).

[0082] In step S400, specifically, the difference between the downward viewing angle between the satellite beam center direction and the line connecting the satellite and the sub-satellite point at each orientation time and the downward viewing angle corresponding to the scene center is calculated, and the calculation formula is shown in formula (4):

[0083] Δβ=β j -β ref (4)

[0084] The antenna range beam pattern is calculated based on the difference Δβ. The calculation formula is shown in formula (5):

[0085]

[0086] The radiation correction function is constructed based on the antenna range beam pattern. The calculation formula is shown in formula (6):

[0087]

[0088] Among them, β j β is the downward angle between the satellite beam center direction and the line connecting the satellite and the sub-satellite point at each azimuth moment; ref is the lower viewing angle corresponding to the center of the scene; Δβ is the difference between the lower viewing angle corresponding to each sampling point of the SAR echo distance at each azimuth and the lower viewing angle corresponding to the center of the scene; W a (β j , H) is the antenna distance beam pattern; L r is the antenna distance aperture length; λ is the signal wavelength; F(W a (β j, H)) is the radiation correction function; α is the adjustable scale factor after image normalization; n is a constant, n∈[0,4].

[0089] Step S500: Based on the radiation correction function F(W a (β j , H)) corrects the radiation distortion of SAR images.

[0090] Step S500 also includes the following steps:

[0091] Step S510, dividing the SAR image into a plurality of sub-SAR images;

[0092] Step S520: setting thread blocks corresponding to the multiple sub-SAR images based on the CUDA programming model;

[0093] Step S530: performing parallel operations on multiple thread blocks based on a graphics processor to complete parallel correction of radiation distortion of multiple sub-SAR images.

[0094] Specifically, first copy the SAR image from the host to the device; then a ×N r The SAR image is divided into n blocks to obtain n sub-SAR image blocks. Based on the parallel computing platform and programming model (Compute Unified Device Architecture, CUDA), the size of the two-dimensional grid corresponding to the n sub-SAR images is set to grid = {G a , G r} thread blocks, each 2D thread block size is block = {B a , B r}; Therefore, the number of pixels that each thread block can process is:

[0095]

[0096] Then, a graphics processing unit (GPU) is used to perform multi-threaded parallel operations on the divided multiple thread blocks, and each pixel unit is compared with the correction function F(W a (β j , H)) coefficients are multiplied, and multi-threaded parallel calculations are performed to realize the parallel correction of radiation distortion of SAR images at the block level, and complete the radiation correction of SAR images. Correcting the radiation distortion of SAR images in parallel through GPU can effectively improve the correction speed to meet the task requirements with high real-time performance.

[0097] In this embodiment, the satellite's ephemeris parameters in the geocentric coordinate system at each azimuth time are obtained to calculate the slant distance between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth time, and the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth time; the satellite beam center pointing at each azimuth time and the downward viewing angle between the satellite and the subsatellite point are further calculated according to the slant distance and ephemeris parameters; the antenna distance beam pattern is calculated according to the satellite beam center pointing at each azimuth time and the downward viewing angle between the satellite and the subsatellite point, and a radiation correction function is constructed, and finally the radiation distortion of the SAR image is corrected in parallel through the radiation correction function. The SAR image quality is effectively improved, the radiation fidelity of the SAR image is met, and the efficiency of the data processing task is also improved, and the task requirements with high real-time performance are taken into account.

[0098] like Figure 5 and Figure 6 As shown, Figure 5 is the imaging result of the SAR image without radiation correction; Figure 6 This is the imaging result of the SAR image corrected by the on-orbit radiation correction method for spaceborne SAR images based on the parallel computing architecture provided by the embodiment of the present invention. By comparison, it can be seen that the clarity of the SAR image after radiation correction is relatively high, the overall brightness is balanced, and the quality of the SAR image is greatly improved.

[0099] According to an embodiment of the present invention, on the other hand, Figure 3 , Figure 4 As shown, a method for on-orbit geometric correction of spaceborne SAR images based on a parallel computing architecture is provided, comprising the following steps:

[0100] Step S600, based on the satellite's position vector, velocity vector and antenna parameters in the satellite's geocentric coordinate system, calculate the coordinates of the ground scene point illuminated by the antenna beam and the Doppler center frequency of the satellite relative to the antenna beam aiming point.

[0101] In step S600, it specifically includes:

[0102] Calculate the rotation matrix from the inertial coordinate system to the satellite orbit coordinate system based on the satellite's ephemeris parameters;

[0103] Based on the satellite's three-axis attitude angle and the satellite's rotation order, calculate the rotation matrix from the satellite's orbital coordinate system to the satellite's body coordinate system;

[0104] Calculate the rotation matrix from the satellite body coordinate system to the antenna coordinate system based on the antenna azimuth angle and antenna elevation angle in the antenna coordinate system;

[0105] Calculate the intersection of the aiming vector of the antenna beam center and the earth's surface in the geocentric coordinate system, obtain the coordinates of the ground scene point illuminated by the antenna beam, and calculate the velocity vector of the antenna beam in the geocentric coordinate system;

[0106] Based on the coordinates of the ground scene point illuminated by the antenna beam and the velocity vector of the antenna beam in the geocentric coordinate system, the Doppler center frequency of the satellite relative to the antenna beam aiming point is calculated.

[0107] Specifically, according to the position vector P of the satellite in the geocentric coordinate system sat =[X sat , Y sat , Z sat ], velocity vector V sat =[Vx sat , Vx sat , Vz sat ], calculate the rotation matrix from the inertial coordinate system to the satellite orbit coordinate system Since the velocity vector V sat_ω The rotation of the earth needs to be taken into account, so V sat_ω =V sat +ω e ×P sat ; Among them, ω e is the angular velocity of the Earth’s rotation, that is, ω e =[0, 0, 7.3e-5].

[0108] According to the rolling angle θ Roll 、Pitch angle to ground θ pitch 、Yaw angle over the ground θ Yaw And the satellite's rotation sequence defines the calculation of the rotation matrix A from the satellite orbit coordinate system to the satellite body coordinate system ve ; Further according to the antenna beam elevation angle Antenna beam azimuth Calculate the rotation matrix A from the satellite body coordinate system to the antenna coordinate system ea Assume that the aiming line of the antenna beam center coincides with the Z axis of the antenna coordinate system, so the aiming vector P of the antenna beam center can be obtained e =[X e , Y e , Z e ]The expression in the central coordinate system:

[0109] P e =A ov ·A ve ·A ea [0, 0, 1] T

[0110] According to the slant distance R between the center of the satellite beam pointing to the center of the scene and the phase center of the antennaref and the aiming vector P at the center of the antenna beam e =[X e , Y e , Z e ] Calculate the coordinates P of the ground scene point illuminated by the antenna beam tar =[X tar , Y tar , Z tar ] and calculate the velocity vector V of the antenna beam in the geocentric coordinate system tar =[Vx tar , Vy tar , Vz tar ], the calculation formula is:

[0111]

[0112] Then, according to the coordinates of the ground scene point illuminated by the antenna beam and the velocity vector of the antenna beam in the geocentric coordinate system, the Doppler center frequency of the satellite relative to the antenna beam aiming point is calculated. The calculation formula is shown in formula (7):

[0113]

[0114] Among them, f d is the Doppler center frequency; P tar is the coordinate of the ground scene point illuminated by the antenna beam; P sat is the satellite position vector; V sat is the velocity vector of the satellite; V tar is the velocity vector of the antenna beam in the geocentric coordinate system; λ is the signal wavelength.

[0115] Furthermore, to calculate the distance R between the centroid of the antenna coordinate system and the intersection of the line of sight of the antenna beam center and the earth's surface, min and the number of echo pulses at the sub-satellite point M b , P e =[X e , Y e , Z e ] and P sat =[X sat , Y sat , Z sat ] into the ellipsoid equation to solve for R min and M b :

[0116]

[0117] Among them, a is the major axis of the earth; b is the minor axis of the earth.

[0118] Step S700: Obtain average elevation information of the target scene, and combine the Doppler center frequency in the range Doppler SAR positioning model to generate coordinates and longitude and latitude information in the geocentric coordinate system corresponding to each pixel point in the SAR image.

[0119] In step S700, specifically including: based on the average elevation information of the target scene and the coordinates and longitude and latitude information of each pixel point in the SAR image in the geocentric coordinate system;

[0120] In the range Doppler SAR positioning model, based on the average elevation information and Doppler center frequency of the target scene, the coordinates of the four corner points in the SAR image in the geocentric coordinate system are calculated and converted into longitude and latitude information; the longitudes Ω of the four corner points are {Ω left_up ,Ω left_down ,Ω right_up ,Ω right_down}, and the latitudes Φ are {Φ left_up , Φ left_down , Φ right_up , Φ right_down};

[0121] Compare the upper left corner point Ω left_up and the latitude value of the lower left corner Ω left_down , determine whether the scene pixels are consistent with the north latitude to south latitude direction;

[0122] In response to the scene pixels being consistent with the north latitude to south latitude direction, the longitude and latitude values ​​corresponding to all pixel points of the SAR image are calculated based on the longitude and latitude of the four corner points and by constructing a kernel function.

[0123] Specifically, in order to more clearly describe the relationship between the SAR image coordinates (i, j) and the latitude and longitude grid coordinates (m, n), the index of the SAR image row number i in azimuth and the index of the column number j in range are first converted into the slant range R and the azimuth time t a To express it, as shown below:

[0124]

[0125] Then, based on the earth model, slant range history and Doppler equation, the range Doppler positioning equation is established to calculate the position coordinates of each pixel point in the SAR image in the geocentric coordinate system. The expression is:

[0126]

[0127] Among them, h is the average elevation of the target scene, (X T ,Y T ,Z T ) is the geographical location coordinates corresponding to each pixel point of the SAR image to be solved.

[0128] Newton iteration method is further used to solve (X T ,Y T ,Z T ),set up Then the iterative formula for solving the range Doppler positioning equation can be expressed by the following generalized operator:

[0129]

[0130] in, represents the inverse matrix at the kth iteration solution, then:

[0131]

[0132] In each iteration (X T ,Y T ,Z T ), in order to speed up the convergence speed, the center point P of each row of pixels in the SAR image is selected tar =[X tar ,Y tar ,Z tar ] as the initial value. Calculate T after each iteration k+1 and T k The maximum difference in the x-axis, y-axis and z-axis directions, that is:

[0133]

[0134] When σ k+1 When it is less than 1e-7, stop the iteration and output (X T ,Y T ,Z T )’s solution.

[0135] The coordinates (X T ,Y T ,Z T ) into the following formula, and then find the longitude Φ corresponding to the pixel point (i, j). The calculation expression is:

[0136]

[0137] The calculation of the latitude Ω corresponding to the pixel point (i, j) requires k iterations until the latitude Ω value converges. The parameters required to calculate the latitude value include the earth radius R e , the extreme flattening of the reference ellipsoid and the radius of curvature of the reference ellipsoid N e , the formula for iterative calculation of latitude is as follows:

[0138]

[0139] Combined with the average elevation information of the target scene, the coordinates of the four corner points of the SAR image in the geocentric coordinate system are calculated based on the range Doppler SAR positioning model, and the coordinates are converted into longitude and latitude information, that is, the longitudes Ω of the four corner points are {Ω left_up ,Ω left_downn ,Ω right_up ,Ω right_down}, and the latitudes Φ are {Φ left_up , Φ left_down , Φ right_up , Φ right_down}.

[0140] Take the upper left corner point Ω of the SAR image left_up and the latitude value of the lower left corner Ω left_down Compare and determine whether the scene pixels are consistent with the north latitude to south latitude direction, and then the longitude and latitude values ​​corresponding to all pixels in the SAR image can be calculated by constructing a kernel function based on the longitude and latitude of these four corner points.

[0141] First, calculate the row slope {k row_Φ , k row_Ω} and the column slope {k col_Φ , k col_Ω} and the starting longitude and latitude values ​​{Ω start , Φ start}, if Ω left_up >Ω left_down ,but:

[0142]

[0143] If Ω left_up <Ω left_down ,but:

[0144]

[0145] The intervals of the longitude and latitude grids ΔΦ and ΔΩ are further set, and the number of distance points N and the number of azimuth points M of the longitude and longitude grids are calculated according to the largest longitude and latitude values ​​and the smallest longitude and latitude values ​​in the scene, so as to obtain the longitude and latitude values ​​corresponding to all pixels of the SAR image. The calculation formula is:

[0146]

[0147] Step S800: Establish an index number mapping relationship between the coordinates of each pixel point in the SAR image in the geocentric coordinate system and the geographic longitude and latitude grid.

[0148] Step S900: construct a kernel function based on a graphics processor to map the SAR image pixels in the imaging plane to the latitude and longitude grids in sequence to complete the geometric correction.

[0149] In steps S800 and S900, after the longitude and latitude grids are divided and the longitude and latitude values ​​of all pixels of the SAR image are obtained, a mapping relationship between the row and column numbers of the SAR image pixel point (i, j) and the longitude and latitude grid (m, n) is established. The image pixel point (i, j) is mapped one-to-one with the thread in the thread block, and the global memory offset of the thread block and thread index can be obtained by mapping the thread and block index to the matrix coordinates, and then the matrix coordinates are mapped to the storage unit of the global memory. Among them, the corresponding relationship between (i, j) and the index identifiers threadIdx.x, threadIdx.y, blockIdx.x, blockIdx.y, blockDim.x and blockDim.y is:

[0150]

[0151] Among them, threadIdx.x and threadIdx.y are the x-dimensional index and y-dimensional index of the thread in each thread block respectively; blockIdx.x and blockIdx.y are the x-dimensional index and y-dimensional index of the thread block in the thread grid respectively; blockDim.x and blockDim.y are the x-direction dimension size and y-direction dimension size of the thread block respectively.

[0152] The kernel function for constructing and calculating the latitude and longitude values ​​of the SAR image pixel point (i, j) and the mapping relationship between the row and column numbers of the latitude and longitude grid (m, n) is:

[0153]

[0154] Finally, the SAR image is copied from the host to the device, and the SAR image is divided into multiple sub-SAR images according to the GPU video memory size. The grid size and thread block size are defined, and the kernel function of the geometric correction mapping is executed in the GPU. All pixel points in the sub-SAR image are mapped to the longitude and latitude grid space, and the final SAR image result after radiation correction and geometric correction is output.

[0155] like Figure 7 , Figure 8 As shown, Figure 7 is the original SAR image of a scene, Figure 8 The SAR image is subjected to radiation correction and geometric correction by the on-orbit radiation and geometric correction method for spaceborne SAR images based on a parallel computing architecture provided by an embodiment of the present invention. Figure 7 and Figure 8 It can be seen that the SAR image after radiation correction and geometric correction has relatively high image quality and can ensure accurate positioning.

[0156] In summary, the on-orbit radiation and geometric correction method for satellite-borne SAR images based on a parallel computing architecture provided by an embodiment of the present invention not only improves the quality of SAR images and the precise positioning of SAR images, but also realizes the efficiency and real-time nature of data processing tasks, and helps the on-orbit SAR processor on board to quickly complete the tasks of radiation correction and geometric correction, and while satisfying the radiation fidelity of SAR image products and the precise positioning of SAR images, it can quickly respond to task requirements and take into account task requirements with high real-time requirements.

[0157] Fig. 9 The schematic diagram of the structure of an embodiment of an electronic device provided by an embodiment of the present invention is shown. The specific embodiment of the present invention does not limit the specific implementation of the electronic device.

[0158] like Fig. 9 As shown, the electronic device may include: a processor (processor) 502 , a communication interface (Communications Interface) 504 , a memory (memory) 506 , and a communication bus 508 .

[0159] The processor 502, the communication interface 504, and the memory 506 communicate with each other via a communication bus 508. The communication interface 504 is used to communicate with other devices such as a client or a network element of another server. The processor 502 is used to execute a program 510, which can specifically execute the relevant steps in the above-mentioned embodiment of the on-orbit radiation and geometric correction method for spaceborne SAR images based on a parallel computing architecture.

[0160] Specifically, the program 510 may include program code including computer-executable instructions.

[0161] The processor 502 may be a central processing unit (CPU), a graphics processing unit (GPU), or an application specific integrated circuit (ASIC), or one or more integrated circuits configured to implement the embodiments of the present invention. The one or more processors included in the electronic device may be processors of the same type, such as one or more CPUs; or may be processors of different types, such as one or more CPUs and one or more ASICs.

[0162] The memory 506 is used to store the program 510. The memory 506 may include a high-speed RAM memory, and may also include a non-volatile memory (non-volatile memory), such as at least one disk memory.

[0163] The program 510 can be specifically called by the processor 502 to enable the electronic device to execute the relevant steps in the above-mentioned embodiment of the on-orbit radiation and geometric correction method for spaceborne SAR images based on parallel computing architecture.

[0164] It can be understood by those skilled in the art that Fig. 9 The structure shown is only for illustration and does not limit the structure of the above-mentioned device. Fig. 9 More or fewer components as shown, or with Fig. 9 Different configurations are shown.

[0165] The embodiment of the present invention also provides a computer-readable storage medium. The method according to the embodiment of the present invention can be implemented in hardware, firmware, or can be implemented as a computer code that can be recorded in a storage medium, or can be implemented as a computer code that is originally stored in a remote storage medium or a non-temporary machine-readable storage medium and will be stored in a local storage medium through a network download, so that the method described herein can be stored in such software processing on a storage medium using a general-purpose computer, a dedicated processor, or programmable or dedicated hardware. Among them, the storage medium can be a magnetic disk, an optical disk, a read-only storage memory, a random access memory, a flash memory, a hard disk or a solid-state hard disk, etc.; further, the storage medium can also include a combination of the above types of memories. It can be understood that a computer, a processor, a microprocessor controller, or programmable hardware includes a storage component that can store or receive software or computer code. When the software or computer code is accessed and executed by a computer, a processor, or hardware, the method shown in the above embodiment is implemented.

[0166] In the specific contents of the above-mentioned specific implementation methods, the various technical features can be combined in any non-contradictory manner. In order to make the description concise, not all possible combinations of the above-mentioned technical features are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0167] The specific contents of the above specific embodiments only express several embodiments of the present invention, and the description is relatively specific and detailed, but it cannot be understood as limiting the scope of the patent of the present invention. It should be pointed out that for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention shall be based on the attached claims.

Claims

1. A method for on-orbit radiation correction of spaceborne SAR images based on a parallel computing architecture, characterized in that: The method comprises: Acquire the ephemeris parameters of the satellite in the geocentric coordinate system at each orientation time, wherein the ephemeris parameters include at least one of a position vector, a velocity vector, a three-axis attitude angle, a pulse repetition frequency, a wave position sampling start time, a signal pulse width, a signal wavelength, an antenna range aperture length, an antenna azimuth aperture length, an antenna beam pitch angle, an antenna beam azimuth angle, a range sampling point number, an azimuth sampling point number, and a signal sampling frequency; Calculate the slant distance between the satellite beam center pointing to the scene center and the antenna phase center at each azimuth moment, as well as the slant distance between each sampling point and the satellite at each azimuth moment; Based on the slant distance between the center of the satellite beam pointing to the scene center point and the antenna phase center at each azimuth moment, the slant distance between each sampling point and the satellite at each azimuth moment, and the ephemeris parameters, the downward viewing angle between the direction of the satellite beam center at each azimuth moment and the line connecting the satellite and the sub-satellite point is calculated; Calculate the antenna range beam pattern based on the satellite beam center direction at each orientation and the downward viewing angle between the satellite and the sub-satellite point, and construct a radiation correction function; The radiation distortion of the SAR image is corrected in parallel based on the radiation correction function.

2. The method according to claim 1, characterized in that The parallel correction of the radiation distortion of the SAR image based on the radiation correction function includes: Dividing the SAR image into a plurality of sub-SAR images; Setting thread blocks corresponding to the plurality of sub-SAR images based on a CUDA programming model; Based on the graphics processor, the plurality of thread blocks are operated in parallel to complete the parallel correction of the radiation distortion of the plurality of sub-SAR images.

3. The method according to claim 1, characterized in that: The three-axis attitude angle includes a rolling angle to the ground, a pitch angle to the ground, and a yaw angle to the ground; the calculation based on the ephemeris parameters of the slant distance between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth moment and the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth moment includes: The number of sub-satellite point echo pulses and the SAR echo acquisition start time are calculated based on the three-axis attitude angle and the velocity vector and position vector of the satellite in the geocentric coordinate system. The calculation formula is shown in formula (1): Based on the echo acquisition start time and the number of echo pulses at the sub-satellite point, the slant distance between the satellite at the beam center and the beam pointing point and the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth time are calculated. The calculation formula is shown in formula (2): Where c is the speed of light; R j R is the slant distance between each sampling point and the satellite in the SAR echo distance direction at each azimuth moment; ref is the slant distance between the satellite beam center pointing to the scene center and the antenna phase center at each orientation moment; center is the start time of echo acquisition; M is the number of sub-satellite echo pulses; t0 is the start time of wave position sampling; PRF is the pulse repetition frequency; T p is the signal pulse width; f s is the signal sampling frequency; N r is the number of sampling points in the distance direction.

4. The method according to claim 3, characterized in that The method of calculating the downward viewing angle between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth moment based on the slant distance between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth moment, the slant distance between each sampling point and the satellite in each azimuth moment, and the ephemeris parameters, comprises: Calculate the distance between the satellite and the sub-satellite point based on the position vector of the satellite and the average radius of the earth; The downward viewing angle between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth time is calculated based on the slant distance between the satellite beam center pointing to the scene center point and the antenna phase center at each azimuth time, the slant distance between each sampling point and the satellite at each azimuth time, the distance between the satellite and the sub-satellite point, and the average radius of the earth; the calculation formula is shown in formula (3): Among them, β j is the downward angle between the satellite beam center direction and the line connecting the satellite and the sub-satellite point at each azimuth; H is the distance between the satellite and the sub-satellite point; R e is the average radius of the Earth.

5. The method according to claim 1, characterized in that The step of calculating the antenna range beam pattern based on the satellite beam center direction at each orientation moment and the downward viewing angle between the satellite and the sub-satellite point and constructing a radiation correction function includes: Calculate the difference between the downward viewing angle between the satellite beam center pointing and the line connecting the satellite and the sub-satellite point at each orientation and the downward viewing angle corresponding to the scene center. The calculation formula is shown in formula (4): Δβ=β j -b ref (4) The antenna range beam pattern is calculated based on the difference Δβ, and the calculation formula is shown in formula (5): A radiation correction function is constructed based on the antenna range beam pattern, and the calculation formula is shown in formula (6): Among them, β j β is the downward angle between the satellite beam center direction and the line connecting the satellite and the sub-satellite point at each azimuth moment; ref is the lower viewing angle corresponding to the center of the scene; Δβ is the difference between the lower viewing angle corresponding to each sampling point of the SAR echo distance at each azimuth and the lower viewing angle corresponding to the center of the scene; W a (β j ,H) is the antenna distance beam pattern; L r is the antenna distance aperture length; λ is the signal wavelength; F(W a (β j ,H)) is the radiation correction function; α is the adjustable scale coefficient after image normalization; n is a constant, n∈[0,4].

6. A method for on-orbit geometric correction of spaceborne SAR images based on a parallel computing architecture, characterized in that: The method further comprises: Calculate the coordinates of the ground scene point illuminated by the antenna beam and the Doppler center frequency of the satellite relative to the aiming point of the antenna beam based on the satellite's position vector, velocity vector and antenna parameters in the satellite's geocentric coordinate system; Obtaining average elevation information of the target scene, and combining the Doppler center frequency in a range Doppler SAR positioning model to generate coordinates and longitude and latitude information in a geocentric coordinate system corresponding to each pixel point in the SAR image; Establishing an index number mapping relationship between the coordinates in the geocentric coordinate system corresponding to each pixel point in the SAR image and the geographic longitude and latitude grid; A kernel function is constructed based on the graphics processor to map the SAR image pixels in the imaging plane to the latitude and longitude grids in sequence to complete the geometric correction.

7. The method according to claim 6, characterized in that The calculating, based on the position vector, velocity vector and antenna parameters of the satellite in the geocentric coordinate system of the satellite, the coordinates of the ground scene point illuminated by the antenna beam and the Doppler center frequency of the satellite relative to the aiming point of the antenna beam comprises: Calculating a rotation matrix from an inertial coordinate system to a satellite orbit coordinate system based on the ephemeris parameters of the satellite; Based on the satellite's three-axis attitude angle and the satellite's rotation order, calculate the rotation matrix from the satellite's orbital coordinate system to the satellite's body coordinate system; Calculate the rotation matrix from the satellite body coordinate system to the antenna coordinate system based on the antenna azimuth angle and antenna elevation angle in the antenna coordinate system; Calculate the intersection of the aiming vector of the antenna beam center and the earth's surface in the geocentric coordinate system, obtain the coordinates of the ground scene point illuminated by the antenna beam, and calculate the velocity vector of the antenna beam in the geocentric coordinate system; Based on the coordinates of the ground scene point illuminated by the antenna beam and the velocity vector of the antenna beam in the geocentric coordinate system, the Doppler center frequency of the satellite relative to the antenna beam aiming point is calculated.

8. The method according to claim 6, characterized in that The method of acquiring the average elevation information of the target scene and combining the Doppler center frequency in the range Doppler SAR positioning model to generate the coordinates and longitude and latitude information in the geocentric coordinate system corresponding to each pixel point in the SAR image includes: Based on the average elevation information of the target scene and the coordinates and longitude and latitude information of each pixel point in the SAR image in the geocentric coordinate system; In the range Doppler SAR positioning model, based on the average elevation information and Doppler center frequency of the target scene, the coordinates of the four corner points in the SAR image corresponding to the geocentric coordinate system are calculated, and the coordinates are converted into longitude and latitude information; the longitudes Ω of the four corner points are {Ω left_up ,Ω left_down ,Ω right_up ,Ω right_down }, and the latitudes Φ are {Φ left_up , Φ left_down , Φ right_up , Φ right_down }; Compare the upper left corner point Ω left_up and the latitude value of the lower left corner Ω left_down , determine whether the scene pixels are consistent with the north latitude to south latitude direction; In response to the scene pixel being consistent with the north latitude to south latitude direction, the longitude and latitude values ​​corresponding to all the pixel points of the SAR image are calculated based on the longitude and latitude of the four corner points and by constructing a kernel function.

9. An electronic device, characterized in that: Also includes: A memory, a processor, a communication interface and a communication bus, wherein the processor, the memory and the communication interface communicate with each other via the communication bus; The memory is used to store at least one executable instruction, and the executable instruction enables the processor to execute the on-orbit radiation and geometric correction method of space-borne SAR images based on a parallel computing architecture according to any one of claims 1 to 8.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer instructions, and the computer instructions are used to enable a processor to implement the on-orbit radiation and geometric correction method for spaceborne SAR images based on a parallel computing architecture according to any one of claims 1 to 8 when executed.

Images