Method and system for processing satellite images

Abstract

The present invention relates to a method and a system for processing satellite images. In particular, said method and system for processing images allow a super-resolution image to be obtained by means of fusing a plurality of satellite images acquired at different instants in time by a sensor on board an orbiting satellite.

Description

[0001] METHOD AND SYSTEM FOR PROCESSING SATELLITE IMAGES

[0002] DESCRIPTION

[0003] OBJECT OF THE INVENTION

[0004] The present invention relates to a method and a system for processing satellite images. In particular, the method and the image processing system enable obtaining a super-resolution image by fusing a plurality of satellite images acquired at different times by a sensor on board an orbiting satellite.

[0005] BACKGROUND OF THE INVENTION

[0006] A satellite image or satellite picture can be defined as the visual representation of the Earth's surface or any other celestial body captured by a sensor mounted on an artificial satellite in orbit.

[0007] Satellite images have a wide variety of applications, such as climate studies, natural disaster monitoring, natural resource exploration, urban planning, the study and exploration of celestial bodies, or military monitoring.

[0008] Typically, artificial satellites have sensors to collect information reflected by the surfaces of celestial bodies and, once they receive this information, they send it to a receiving equipment on Earth, where it is properly processed to obtain satellite images.

[0009] Satellite images can be captured at many different resolutions, such as spatial, spectral, temporal, radiometric, or geometric. Specifically, spatial resolution refers to the pixel size in the satellite image that represents a portion of the celestial body's surface.

[0010] High-resolution spatial imagery allows for detailed work at scale, but obtaining such high-resolution images requires significant computational resources. This necessitates specialized and advanced hardware and software on the receiving equipment. Currently available onboard satellites, the technical resources are insufficient to process the raw images acquired by the sensors and thus produce high-resolution satellite imagery. Achieving this goal would require more specialized and advanced hardware, but weight and space limitations on the satellite prevent the deployment of such resources.

[0011] Therefore, there is a need for a satellite image processing method that can be executed by a system with limited hardware, such as that found on board a satellite. This way, it would not be necessary to send the information acquired by the sensors on board the satellite to a ground receiver for processing; instead, the images would be sent already processed in real time to the ground receiver immediately after being acquired by the sensor.

[0012] DESCRIPTION OF THE INVENTION

[0013] The present invention proposes a solution to the above problems by means of a computer-implemented method for processing satellite images according to claim 1, a system for processing satellite images according to claim 11, a celestial body observation camera according to claim 14, and a satellite according to claim 15. Preferred embodiments of the invention are defined in the dependent claims.

[0014] A first inventive aspect provides a computer-implemented method for processing satellite images, the method comprising the steps of:

[0015] (a) receiving a plurality of satellite images acquired by a sensor on board a satellite and telemetry data, the telemetry data comprising a plurality of positions and attitudes of the satellite with respect to an inertial reference frame centered on a celestial body at a plurality of telemetry time instants;

[0016] (b) align a pair of consecutive satellite images from the plurality of satellite images by carrying out the following sub-steps:

[0017] (b.1) select a first satellite image and a second satellite image to align,

[0018] (b.2) Calculate an initial horizontal displacement and an initial vertical displacement of the second satellite image with respect to the first satellite image:

[0019] - estimating a first time instant in which the sensor acquired the first satellite image from the telemetry data;

[0020] - estimating the position and attitude of the satellite at the first estimated time instant from the telemetry data;

[0021] - estimating a second time instant in which the sensor acquired the second satellite image from the telemetry data;

[0022] - estimating the position and attitude of the satellite at the second time instant estimated from the telemetry data;

[0023] - calculating the initial horizontal displacement and the initial vertical displacement from the positions and attitudes of the satellite at the first time instant and at the second time instant;

[0024] (b.3) Starting from the initial horizontal displacement and the initial vertical displacement, calculate the total horizontal displacement and the total vertical displacement that must be applied to the second satellite image to minimize the distance between said second satellite image and the first satellite image, the distance measured according to a first pre-established standard, preferably the Euclidean standard;

[0025] (b.4) align the first satellite image and the second satellite image by applying the total horizontal displacement and the total vertical displacement to the second satellite image;

[0026] (c) repeat step (b) for all pairs of consecutive satellite images from the plurality of satellite images, wherein in each iteration of step (b):

[0027] - the first satellite image selected in sub-stage (b.1) corresponds to the second satellite image of the previous iteration;

[0028] - the total horizontal displacement obtained in sub-stage (b.3) must be added to the total horizontal displacement obtained in sub-stage (b.3) of the previous iteration; and

[0029] - the total vertical displacement obtained in sub-stage (b.3) must be added to the total vertical displacement obtained in sub-stage (b.3) of the previous iteration;

[0030] (d) select as the reference image the first satellite image selected in sub-stage (b.1) of the first iteration of stage (b), and obtain at least one common region between the reference image and each of the second satellite images shifted in sub-stage (b.4) of all iterations of stage (b), resulting in a set of common regions;

[0031] (e) obtain an oversampled fused image by merging the reference image and the set of common regions by carrying out the following sub-steps: (e.1) generate an oversampled image, wherein the resolution of the oversampled image is the resolution of the reference image multiplied by a predefined oversampling factor;

[0032] (e.2) for each pixel of the oversampled image, determine the distance from the center of that pixel to the center of at least one pixel in each region of the set of regions and to the center of at least one pixel in the reference image, wherein the distance is measured according to a second pre-established norm, preferably the Euclidean norm;

[0033] (e.3) for each pixel of the oversampled image, estimate a weight for at least one pixel of each region of the set of regions and for at least one pixel of the reference image based on the determined distances; and (e.4) obtain the fused oversampled image, wherein the value of each pixel of the fused oversampled image is the sum of the values ​​of at least one pixel of each region of the set of regions and of at least one pixel of the reference image weighted according to the weights estimated in the previous sub-step.

[0034] First, the method requires receiving two or more satellite images acquired by a sensor on board a satellite and the telemetry data acquired during the acquisition of the satellite images.

[0035] Throughout this document, it will be understood that telemetry data comprises a plurality of satellite positions and attitudes determined at a plurality of time moments.

[0036] In a realization, telemetry data, such as the attitude of the satellite, is determined by one or more of the following devices on board the satellite: star trackers, solar sensors, or a magnetometer.

[0037] Throughout this document, the attitude of a satellite will be understood as the orientation or pointing of said satellite with respect to an inertial reference system centered on a celestial body.

[0038] Throughout this document, the inertial reference system will be understood as that with respect to which an element is positioned and oriented in relation to a celestial body, the inertial reference system being centered on the celestial body, with its origin at the center of mass of said celestial body.

[0039] Throughout this document, a celestial body will be understood to be any natural object located in outer space. Examples of celestial bodies are planets, such as Earth; natural satellites, such as the Moon; or stars.

[0040] Next, the images are aligned in consecutive pairs; for example, by registering the first satellite image with the second, the second with the third, and so on until the last satellite image of the plurality of satellite images is reached.

[0041] Throughout this document, aligning two two-dimensional satellite images will be understood as finding the displacement of one relative to the other along the horizontal and vertical axes that minimizes the distance between the two images. Throughout this document, distances are measured using one or more pre-established standards, preferably the Euclidean standard.

[0042] The image alignment sub-stages according to the invention are implemented to be able to be executed with hardware limited in terms of processing capacity, such as that which can be found on board a satellite.

[0043] To achieve this, the alignment procedure relies on a prior calculation of horizontal and vertical displacements from the received telemetry data. In other words, the alignment of each pair of consecutive images does not start from scratch, but rather from an initial estimate that takes into account how much the sensor on the satellite has moved between the acquisition of the first image of the pair to be aligned and the acquisition of the second image of the pair to be aligned. This significantly reduces the execution time for image alignment, thus lowering the technical requirements of the system responsible for the alignment. Consequently, this stage of the method can be carried out using the existing hardware on a satellite.

[0044] In particular, image alignment (stage b) requires the execution of the following sub-stages.

[0045] The first step of the alignment, sub-stage (b.1), requires selecting a pair of images from the plurality of satellite images: a first satellite image and a second satellite image.

[0046] Next, in step (b.2), an initial horizontal displacement and an initial vertical displacement are calculated to be applied to the second satellite image with respect to the first satellite image. These initial displacements are the basis for the alignment procedure as a tai.

[0047] To calculate the initial horizontal and vertical displacements, the times at which the sensor acquired the satellite images to be aligned must first be estimated. Similarly, the satellite's position and attitude at these estimated times must be determined; that is, the position and attitude the satellite had when the sensor acquired the images to be aligned. This estimation is performed using telemetry data.

[0048] In a realization, to estimate these time points, as well as the satellite's position and attitude, it is necessary to interpolate the previously received telemetry data. In general, the image reception rate is much higher than the telemetry data reception rate. For example, the sensor might capture 50 images per second, while the satellite's position and attitude data, along with the acquisition time point, are received at an average frequency of 4 Hz, which can reach at most 10 Hz. Therefore, for the satellite images to be aligned, the satellite's position and attitude, as well as the acquisition time point, can be estimated by interpolating the received telemetry data.

[0049] Once the time instants and the attitudes and positions of the satellite are known, the method continues by calculating the initial horizontal and vertical displacements from the variations in the positions and attitudes of the satellite between the two estimated time instants.

[0050] The method continues in sub-step (b.3), where, starting from the initial horizontal and vertical displacements, the total horizontal and vertical displacements that must be applied to the second satellite image are calculated to minimize the distance between the second satellite image and the first satellite image. As mentioned previously, the distance is measured according to a pre-established standard, preferably the Euclidean standard. To finalize the alignment, in sub-step (b.4), the first and second satellite images are aligned by applying the total horizontal and vertical displacements obtained in sub-step (b.3) to the second satellite image.

[0051] Next, step (c) of the method is executed. In this step, the alignment defined in step (b) must be repeated for all pairs of consecutive images from the plurality of satellite images.

[0052] In each new iteration of stage (b), the satellite image selected in sub-stage (b.1) as the first satellite image must correspond to the second satellite image selected in sub-stage (b.1) of the previous iteration of stage (b). Thus, if in the first iteration of stage (b) the second satellite image was aligned with respect to the first satellite image, in the second iteration of stage (b) the third satellite image will be aligned with respect to the second satellite image, and so on until the last satellite image of the plurality of satellite images is aligned with respect to the penultimate satellite image of the plurality of satellite images.

[0053] Additionally, in each iteration of stage (b) other than the first iteration, the horizontal displacement obtained in sub-stage (b.3) must be added to the horizontal displacement obtained in sub-stage (b.3) of the previous iteration. Similarly, in each iteration of stage (b) other than the first iteration, the vertical displacement obtained in sub-stage (b.3) must be added to the vertical displacement obtained in sub-stage (b.3) of the previous iteration. In this way, all images are aligned with respect to the image selected as the first satellite image in the first iteration of stage (b) of the method.

[0054] For example, if aligning the second image of the plurality of satellite images with respect to the first image of that plurality results in a horizontal displacement X and a vertical displacement Y; and aligning the third image of the plurality of satellite images with respect to the second image results in a horizontal displacement X' and a vertical displacement Y', then the total horizontal displacement calculated in step (b.3) of the second iteration of step (b) will be X+X' and the total vertical displacement calculated in step (b.3) of the second iteration of step (b) will be Y+Y'. In this way, both the second and third images would be recorded with respect to the first image of the plurality of satellite images.

[0055] Once the alignment of all image pairs from the plurality of satellite images has been carried out, the method continues with the execution of stage (d). In this stage, the first satellite image selected in sub-stage (b.1) of the first iteration of stage (b) is selected as the reference image. That is, the satellite image with respect to which the rest of the satellite images have been aligned after performing the complete chain of image pair alignments.

[0056] Taking into account this reference image, stage (d) of the method requires obtaining a common region between the reference image and each of the second satellite images shifted in sub-stage (b.4) of all iterations of stage (b), resulting in a set of common regions.

[0057] To do this, one of the second shifted satellite images is compared with the reference image to find at least one common region between them. In this way, at least one common region is selected for the second shifted satellite image, which will contain the pixels of the common region that correspond to the pixels in the reference image.

[0058] This obtaining of reference regions is repeated for each second displaced satellite image in sub-stage (b.4) of each iteration of stage (b), thus obtaining a set of common regions (at least one for each second displaced satellite image).

[0059] The method concludes with obtaining a super-resolution fused image by merging the reference image and the set of common regions. This fusion requires several sub-steps:

[0060] First, in sub-step (e.1), an oversampled image is generated with a resolution equal to that of the reference image multiplied by a predefined oversampling factor. This oversampled image, for example, could be an image in which all pixels have a predetermined value, such as zero or one. Next (sub-step e.2), for each pixel in the oversampled image, the distance from its center to the center of one or more pixels in the images to be merged (common regions and reference image) is determined. Specifically, once the images to be merged are positioned, the distances from the center of a pixel in the oversampled image to the centers of one or more pixels in each common region and to the centers of one or more pixels in the reference image are calculated. The distance is calculated according to a pre-established standard, preferably the Euclidean standard.

[0061] In a realization, once the images to be merged are positioned, the distances are calculated from the center of a pixel of the oversampled image to the centers of one or more pixels of each common region that coincide and / or are adjacent or bordering to the pixel of the oversampled image and to the centers of the pixels of the reference image that coincide and / or are adjacent or bordering to the pixel of the oversampled image.

[0062] This action is repeated for all pixels of the oversampled image.

[0063] The method continues in sub-step (e.3) by estimating, for each pixel of the oversampled image, a weight for each of the pixels in the set of common regions and the reference image for which distances have been calculated with respect to the pixel of the oversampled image. This weight depends on the determined distance; for example, the smaller the calculated distance, the greater the estimated weight for that pixel in the set of common regions or the reference image.

[0064] Finally, the method obtains an oversampled fused image in step (e.4). This oversampled fused image is based on the oversampled image, in which the value of each pixel is replaced by the sum of the pixel values ​​from the set of common regions and the reference image for which the distance from that pixel in the oversampled image has been calculated, but weighted according to the weights estimated in step (e.3).

[0065] Thanks to this method, a super-resolution image is obtained from a set of lower-resolution satellite images using hardware with limited processing capacity, such as that found on board a satellite. Image fusion, which produces the super-resolution image, requires prior alignment of the original satellite images, a computationally expensive operation. However, with the method of the invention, the alignment is based on a prior estimation of the horizontal and vertical displacements to be applied to one image relative to the other, where this estimation is performed using telemetry data.Since image alignment does not start from scratch, but is merely an adjustment of the initial displacements, the computational calculations are greatly sped up, which translates into a lower requirement for hardware computational capacity and, therefore, makes it possible for the method to be executed on a system on board a satellite.

[0066] In one realization, in step (b.3), the displacement is determined by minimizing a linear cross-correlation between the first satellite image and the second satellite image.

[0067] In a particular embodiment, sub-step (b.2) the calculation of the initial horizontal displacement and the initial vertical displacement from the positions and attitudes of the satellite at the first time instant and at the second time instant is carried out by means of the following sub-steps:

[0068] (1) in the first iteration of stage (b), for each pair of consecutive telemetry time instants, t teil yt tei2 , from the telemetry data:

[0069] - Calculate the pointing vectors of the sensor's center pixel and at least two pixels from two sensor segments, which are perpendicular to each other and intersect at the sensor's center, with respect to the inertial reference system based on the satellite's position and attitude at the first telemetry time instant. and at the second telemetry time instant, t tei2 ;

[0070] - obtain the projection of the sensor center pixel and at least two pixels of the sensor segments onto the surface of the celestial body by calculating the intersection of the pointing vectors of said pixels with a reference geometric figure of a representative geodetic model of the celestial body;

[0071] - define a projected reference system as the projection of at least two pixels of the sensor segments and the projection of the sensor center pixel onto the surface of the celestial body assuming that: either the surface of the celestial body is flat; or the projected reference system is an orthogonal reference system formed by the projections of the pixels of the sensor segments and the projection of the sensor center pixel; and the projection of the sensor center pixel is the center of the projected reference system;

[0072] - For the sensor's center pixel: calculate the effective height as the magnitude of its pointing vector, where the terminal point of the pointing vector is its intersection with the surface of the celestial body, and calculate the ground sampling distance, GSD, which represents the size of the projection of the sensor's center pixel onto the surface of the celestial body: sensor pixel size

[0073] GSD = effective height ■ - - - - — — - sensor focal length

[0074] - calculate a first linear velocity, vf of the projection of the pixel of the center of the sensor due to the rotation of the celestial body, where the first linear velocity v^ is referenced in the inertial reference system, and the first linear velocity v^ is calculated from the angular velocity, co, of the celestial body and the perpendicular distance from the projection of the pixel of the center of the sensor to the axis of rotation of the celestial body, R : where cJ R are vectors referenced in the inertial reference system;

[0075] - Calculate a second linear velocity, vj, of the projection of the pixel at the center of the sensor due to the movement of the satellite: where: o pos^ is the position vector of the projection of the pixel of the center of the sensor at the second telemetry time instant, t tet2, the vector being referenced in the inertial reference system; I pos is the position vector of the projection of the pixel of the center of the sensor at the first temporal instant of telemetry, the vector being referenced in the inertial reference frame; I vj is referenced in the inertial reference frame;

[0076] - Calculate the total linear speed, v^, of the projection of the pixel at the center of the sensor in pixels per second:

[0077] - reference the total linear speed obtained in the previous sub-stage in the projected reference system;

[0078] (2) In the first iteration of stage (b), calculate the horizontal linear velocities, v H (t), and vertical, v v(t), as a function of time by interpolating the horizontal and vertical components of the total linear velocities, v^, obtained for all pairs of consecutive telemetry time instants; and

[0079] (3) In any iteration of step (b), for the pair of consecutive satellite images to be aligned, estimate the initial horizontal displacement and the initial vertical displacement by integrating the horizontal and vertical linear velocities with respect to time, v H (t) yv K (t), respectively, between the first time instant and the second time instant.

[0080] In this realization, a way of calculating the initial horizontal displacement and the initial vertical displacement from telemetry data is detailed; in particular, from the positions and attitudes of the satellite at the first time instant and at the second time instant that were estimated in sub-stage (b.2) of the method.

[0081] In steps (1) and (2) of this implementation, the horizontal and vertical velocities are estimated as a function of time from the sensor's projection onto the surface of the celestial body, taking into account all telemetry data, that is, the plurality of telemetry time points and the corresponding positions and attitudes of the satellite at those telemetry time points. These steps (1) and (2) must be performed only once during the execution of the method, that is, in the first iteration of stage (b).

[0082] Once the horizontal and vertical velocities as a function of time are available, in step (3) of this implementation, the initial horizontal and vertical displacements that would have to be applied to a second satellite image to align it with a first satellite image are estimated. This step must be performed in every iteration of stage (b), so that different initial horizontal and vertical displacements are calculated for each pair of images to be aligned in each iteration of stage (b).

[0083] During step (1) of this realization, for each pair of consecutive telemetry time instants, t teil yt tei2 From the telemetry data, the pointing vectors of the sensor's center pixel and at least two pixels from two sensor segments that intersect at the sensor's center and are perpendicular are calculated. In particular, for each telemetry time instant, t teil yt tei2The satellite's positions and attitudes are identified according to the telemetry data, and knowing these positions and attitudes of the satellite, the pointing vectors are calculated.

[0084] Preferably, the two pixels in the segments belong to the sensor's edges; that is, to the set of sensor pixels that delimit Io above, below, left, and right. For example, one pixel from the top edge of the sensor and one pixel from the right edge of the sensor.

[0085] Throughout this document, the pointing vector of a sensor pixel will be understood as that vector originating at the center of the pixel, which is directed according to the optical direction of the pixel, which intersects the surface of the celestial body at a specific point and which are referenced in the inertial reference system.

[0086] The next step is to obtain the projection of the sensor's center pixel and the segment pixels onto the surface of the celestial body. To do this, the intersection of each pointing vector with a reference geometric figure from a geodetic model representative of the celestial body is calculated. For example, the WGS 84 reference ellipsoid if the celestial body is Earth, this ellipsoid belonging to a known geodetic model.

[0087] From the projections of the pixels of the segments and the projection of the pixel at the center of the sensor, a projected reference system is defined assuming that the surface of the celestial body is flat, that the projected reference system formed by the projections of the pixels of the segments and the projection of the central pixel is orthogonal, and that the projection of the pixel at the center of the sensor is the center of the projected reference system.

[0088] At this point, it's important to consider that, due to the curvature of some celestial bodies, such as Earth, and the satellite's pointing, the reference system formed by the projections of the segments intersecting at the sensor's center and the projection of the sensor's center pixel is distorted and not orthogonal. However, within the context of defining the projected reference system, these errors are negligible. Therefore, the assumptions are made that the celestial body's surface is flat and that the reference system is orthogonal, with its center being the projection of the sensor's center pixel.

[0089] Next, for the pixel at the center of the sensor, the effective height is calculated as the magnitude of its pointing vector, understanding that the terminal point of the pointing vector corresponds to the intersection of said pointing vector with the surface of the celestial body.

[0090] Throughout this document, the effective height of a pixel will be understood as the distance between the pixel of the sensor on board the satellite and its projection onto the surface of the celestial body.

[0091] Additionally, the ground sampling distance, GSD, is calculated, which represents the size of the projection of said pixel from the center of the sensor onto the surface of the celestial body using the following expression: sensor pixel size

[0092] GSD = effective height ■ - - - - — — - sensor focal length

[0093] The sensor pixel size and sensor focal length are known commercial parameters, which depend on the type of sensor analyzed.

[0094] The method continues by calculating a first linear velocity, vf, of the projection of the sensor's center pixel resulting from the celestial body's rotation. This velocity is referenced to the inertial reference frame and is calculated from the angular velocity, df, of the celestial body and the perpendicular distance from the projection of the sensor's center pixel to the celestial body's rotation axis, R. d) and R are vectors referenced in the inertial reference frame. For example, the angular velocity of the Earth is 7.27 x 10' 5 The perpendicular distance from the projection of the sensor's center pixel to the Earth's rotation axis is 6,378 kilometers at the Earth's equator. The next step in the method is to calculate a second linear velocity, vj, of the projection of the sensor's center pixel due to the movement of the satellite on which the sensor is mounted: pos and pos are the position vectors of the projections of the sensor center pixel at the second telemetry time instant, t tei2 , and in the first instant of telemetry,

[0095] Finally, the total linear speed, v^, of the projection of the sensor's center pixel is calculated, converted to pixels per second using the GSD:

[0096] This total linear speed This is referenced in the inertial reference frame, and this reference must be modified to the projected reference frame. For example, by projecting the horizontal and vertical components of the total linear velocity onto the projected reference frame.

[0097] This step (1) of the method must be performed for each pair of telemetry time points, thus obtaining a plurality of total linear velocities. In other words, the method requires a first calculation of the total linear velocity from telemetry time points one and two of the telemetry data, a second calculation of the total linear velocity from telemetry time points two and three of the telemetry data, and so on until a final calculation of the total linear velocity from the penultimate and final telemetry time points.

[0098] In step (2) of this realization, the total linear velocity as a function of time is calculated on the horizontal and vertical axes, that is, the horizontal linear velocities, v H (t), and vertical, v v(t). To do this, the horizontal and vertical components of the plurality of total linear velocities, v^, obtained in step (1) are interpolated for all pairs of consecutive telemetry time instants.

[0099] Finally, in step (3), for the pair of consecutive satellite images to be aligned acquired at a first and a second time instant, the method estimates the initial horizontal displacement and the initial vertical displacement by integrating the horizontal and vertical linear velocities with respect to time, v H (t) yv K (t), respectively, between the first time instant and the second time instant.

[0100] Advantageously, this implementation takes into account all telemetry data as well as the rotation of the celestial body to very accurately estimate the initial horizontal and vertical displacements that must be applied to the second satellite image to align it with the first satellite image. The more accurate this estimate, the faster the image alignment is performed and the fewer computational resources the alignment operation requires, which translates into the feasibility of executing the method using the hardware available on board a satellite.

[0101] In a particular embodiment, sub-step (b.3) is carried out by means of the following sub-steps:

[0102] (b.3. i) For the first satellite image, generate a first image pyramid by subsampling the first satellite image with N predefined subsampling factors, where N is greater than or equal to one, where the first image pyramid comprises N+1 levels, and where:

[0103] • The first level of the first image pyramid comprises the first satellite image, and

[0104] • the i-th level of the first image pyramid comprises an image subsampled from the first satellite image with the k-th subsampling factor, where i = 2,... ,N+1 and k = i-1 ;

[0105] (b.3.ii) For the second satellite image, generate a second image pyramid by subsampling the second satellite image with M predefined subsampling factors, where M is greater than or equal to one, where the second image pyramid comprises M+1 levels, and where:

[0106] • M is less than or equal to N,

[0107] • The first level of the second image pyramid comprises the second satellite image, and

[0108] • the j-th level of the second image pyramid comprises an image subsampled from the second satellite image with the m-th subsampling factor, where j = 2,... ,M+1 and m = j-1 ;

[0109] (b.3. iii) Starting from the initial horizontal displacement and the initial vertical displacement, determine the horizontal displacement and the vertical displacement that must be applied to the image of level M+1 of the second pyramid of images to minimize the distance between said image and the image of the first pyramid of images of the same level, the distance measured according to a third pre-established norm, preferably the Euclidean norm;

[0110] (b.3.iv) Starting from the horizontal displacement and the vertical displacement determined in the previous sub-stage, determine the horizontal displacement and the vertical displacement that must be applied to the image of the lower level of the second pyramid of images to minimize the distance between said image and the image of the first pyramid of images of the same level, the distance measured according to the third pre-established rule, and where:

[0111] • if that lower level is level 1 of the second image pyramid, set the horizontal offset as total horizontal offset and the vertical offset as total vertical offset; or

[0112] • If that lower level is not level 1 of the second image pyramid, repeat step (b.3.iv).

[0113] This implementation describes a particular way of aligning satellite images in an efficient manner in terms of execution time.

[0114] As previously mentioned, the image alignment sub-stages according to the invention are implemented to be executed with hardware limited in processing capacity, such as that found on board a satellite. To this end, this embodiment seeks to align images starting from versions with a lower resolution than the original, so that the required computation times are manageable by such hardware with limited processing capacity.

[0115] First (sub-stage b.3.i) a first image pyramid with N+1 levels is generated. The first level of the pyramid comprises the first satellite image, which is subsampled with one or more predefined subsampling factors to generate the images of the remaining levels of the first image pyramid.

[0116] The higher the level of the pyramid, the greater the subsampling factor used to generate the subsampled image contained in that level.

[0117] Throughout this document, it will be understood that subsampling an image involves reducing its spatial resolution by combining sets of pixels according to a subsampling factor. For example, subsampling an image with a subsampling factor of 2 means combining two adjacent horizontal pixels and two adjacent vertical pixels to generate a single pixel in the subsampled image.

[0118] The method continues in stage (b.3.ii) generating a second pyramid of M+1 level images from the second satellite image.

[0119] The first level of the second pyramid comprises the second satellite image itself, which is subsampled with one or more subsampling factors to generate the images of the remaining levels of the pyramid.

[0120] The higher the level of the pyramid, the greater the subsampling factor used to generate the subsampled image contained in that level.

[0121] The number of levels in the second image pyramid can be equal to or less than the number of levels in the first image pyramid.

[0122] In the next step of the method, sub-stage (b.3.iii), an alignment is performed between the lowest resolution image of the second image pyramid (i.e., the image at level M+1 of the second image pyramid) and the image of the first image pyramid of equal resolution (i.e., the image at level M+1 of the first image pyramid). This alignment is based on the initial horizontal and vertical displacements estimated from telemetry data.

[0123] The result of this alignment comprises the horizontal and vertical displacement that must be applied to the lower-resolution satellite image of the second pyramid so that the distance between it and the lower-resolution image of the first pyramid is minimized. As mentioned previously, the distance is measured according to a pre-established standard, preferably the Euclidean standard.

[0124] Next, the alignment of the images of the remaining pyramid levels is carried out sequentially in stage (b.3.iv). Specifically, an alignment is performed between the image of level M of the second pyramid and the image of level M of the first pyramid. This alignment is based on the horizontal and vertical displacements determined for the images of the upper level, thus greatly speeding up the process and requiring less computational time and resources.

[0125] Thus, the result of this alignment of sub-stage (b.3.iv) comprises the horizontal displacement and the vertical displacement that must be applied to the image of level M of the second pyramid so that the distance between this and the image of level M of the first pyramid is minimal.

[0126] Next, the method evaluates whether the last alignment carried out corresponds to the alignment of the images of the first level of the pyramids (that is, of the first and second image pyramids):

[0127] - If the last recording was carried out on the images of level 1 of the pyramids, the horizontal displacement determined in the last iteration of stage (b.3.iv) is established as the total horizontal displacement and the vertical displacement determined in the last iteration of stage (b.3.iv) as the total vertical displacement; and the method continues in stage (b.4); or

[0128] - Otherwise, step (b.3.iv) is repeated sequentially until the alignment of the level 1 pyramid images is carried out, at which point the method continues in step (b.3.iv).

[0129] It should be mentioned that, in each iteration of stage (b.3. iv), the alignment is performed on the images included in the pyramid level immediately below the level on which the last alignment was carried out, the alignment also starting from the result (horizontal and vertical displacement) determined for the last alignment; for example, if in the last iteration of stage (b.3.iv) the alignment was carried out on the images of level M of the pyramids, in the next iteration the alignment must be carried out on the images of level M-1 taking into account the result of the alignment of the images of level M, and so on.

[0130] In one implementation, for all iterations of step (b) other than the first, sub-step (b.3.i) of generating a first image pyramid from the first satellite image of the pair to be registered is omitted, since this first image pyramid corresponds to the second image pyramid generated in the previous iteration. As already mentioned, the image fusion that results in the super-resolution image requires prior alignment of the original satellite images, a computationally expensive operation. However, thanks to this implementation, the alignment starts from undersampled images, which reduces the computational cost of this action, and the result of this alignment serves as a starting point for subsequent alignments of the higher-resolution images in the generated pyramids, which greatly speeds up the computational calculations of these subsequent alignments.

[0131] In one realization, in steps (bSiii) and (b.3.iv), the horizontal and vertical displacements are determined by minimizing a linear cross-correlation between the first satellite image and the second satellite image of the corresponding pyramid levels.

[0132] In a particular embodiment:

[0133] - in the first iteration of stage (b), the number of levels in the first image pyramid is equal to the number of levels in the second image pyramid; and

[0134] - In the remaining iterations of stage (b), the number of levels of the remaining image pyramids is less than the number of levels of the image pyramids generated in the first iteration.

[0135] In this realization, the number of levels in the first and second image pyramids generated in the first iteration of stage (b) is the same. For the remaining iterations of stage (b), the number of levels in the generated pyramids is lower than in the first iteration, in order to reduce the computational cost of the alignment.

[0136] In a particular embodiment, the method further comprises a step prior to step (b) of correcting the plurality of satellite images.

[0137] In this realization, prior to image capture, the images may undergo at least one correction process. This could include, for example, the removal of fixed and / or variable noise, compensation for the sensor's non-uniform light effects, optical correction, and / or geometric correction. Advantageously, the corrected images are of higher quality than the original or raw images, resulting in a super-resolution image of superior quality obtained through the complete execution of the method.

[0138] In a particular embodiment, the method comprises a stage prior to stage (b) of fixed noise elimination in the plurality of satellite images, the fixed noise elimination stage being carried out by means of the following sub-stages:

[0139] - determining a dark image by averaging a plurality of images acquired by the sensor in the absence of illumination; and

[0140] - subtracting the dark image from each of the satellite images.

[0141] In this realization, the correction stage of the method involves removing fixed noise from the raw satellite images.

[0142] Throughout this document, fixed noise will be understood as a pattern of interference or distortion that appears consistently in all satellite images captured by the same sensor. This is due to the sensor's dark current and various thermal phenomena, primarily caused by the electrical current flowing through the electronic components when light is not striking the sensor. This noise manifests as dots, lines, or patterns unrelated to the information the image is attempting to represent, which can degrade the image's visual quality and sharpness.

[0143] To correct for fixed noise, the method comprises the sub-stages of 1) determining a dark image and 2) subtracting that dark image from each satellite image.

[0144] To determine the dark image, multiple images are acquired in the absence of sensor illumination, and these multiple images are averaged. Image acquisition can be performed by the sensor on Earth prior to satellite launch, or directly by the sensor on board the satellite once it is in orbit.

[0145] In a particular embodiment, the method further comprises a step prior to step (b) of correcting the plurality of satellite images (1) to compensate for the non-uniform effect on the sensor light, the correction step being carried out by the following sub-steps: - acquiring a plurality of imageries of a uniformly illuminated object and subtracting, from each one, the dark image;

[0146] - determining a corrective image by averaging the plurality of images from the previous sub-stage;

[0147] - normalizing the corrective image; and

[0148] - dividing each of the fixed noise-free satellite images by the normalized corrective image.

[0149] The sensors on board satellites can have variations in their sensitivity in different areas, which can result in a final image with irregularities in light intensity.

[0150] In this realization, the correction stage of the method involves reducing or eliminating the non-uniform light effect that the sensor has on different parts of the image. This correction is known as "flat correction".

[0151] This correction involves, firstly, acquiring multiple images of a uniformly illuminated object, or flat images. The acquisition of these multiple flat images can be performed by the sensor on Earth, prior to the satellite's launch, or directly by the sensor on board the satellite once it is in orbit.

[0152] These flat images are then subtracted from the dark image according to the previous realization to eliminate the fixed noise that appears in them.

[0153] By averaging the noise-free flat images, a corrective image is obtained that must be normalized so that all values ​​are around unity.

[0154] Finally, the corrected image is obtained by dividing each of the fixed-noise satellite images by the corrected image. The fixed-noise satellite images are those from which the dark image has been subtracted.

[0155] Advantageously, flat correction can help minimize unwanted variations in light intensity, which in turn can contribute to reducing noise in the final image. In one implementation, the sensor comprises a plurality of filters covering a plurality of pixels. In such a case, normalizing noise-free flat images requires 1) calculating the average value of the pixels for each filter on the sensor, and 2) dividing the number of pixels in the fixed noise-free flat images corresponding to the filter area by that average value.

[0156] In a particular embodiment, the method comprises a stage prior to stage (b) of variable noise elimination in each image of the plurality of satellite images, the variable noise elimination stage being carried out by the following sub-stages: i) calculating the median of each row of the satellite image, and replacing the value of the pixels in the row that exceed a predefined threshold by the value of said median; ii) estimating the modulus and phase of the Fourier transform of the satellite image resulting from the previous sub-stage; iii) calculating the median of the modulus of the Fourier transform; iv) determining the pixels of the modulus of the Fourier transform that are outside a predetermined radius measured from the center of said modulus of the Fourier transform; v) replacing the values ​​of the pixels determined in the previous sub-stage by the median calculated in sub-stage iii);vi) obtain a smoothed satellite image by applying the inverse Fourier transform; vii) obtain a high-frequency image by subtracting the smoothed satellite image from the satellite image; viii) for each row of the satellite image: calculate the median of the corresponding row of the high-frequency image, and subtract this median from each pixel of the row of the satellite image.

[0157] In this realization, the correction stage of the method involves reducing or eliminating variable noise or read noise present in the raw satellite images. This correction can also be applied to images where fixed noise has been corrected or to images where flat correction has been applied.

[0158] The correction of variable noise, first of all, requires detecting pixels that exceed a predefined noise threshold. If such pixels are detected, their values ​​are replaced by the median of the row in the satellite image in which they are located.

[0159] Next, the modulus and phase of the Fourier transform of the resulting image are estimated, the median of the modulus is calculated, and the pixels at the corners of the modulus of the Fourier transform are replaced by this median.

[0160] The corners of the modulus should be understood as the pixels that are outside a predetermined radius measured from the center of the Fourier transform modulus. For example, the pixels that are outside a radius of 50 pixels.

[0161] The method continues by obtaining a smoothed satellite image (of low frequencies) by applying the inverse Fourier transform to the modified modulus in stage v) and the phase estimated in stage ii). This image is subtracted from the original satellite image, thus obtaining a high-frequency image.

[0162] Finally, the median of each row in the high-frequency image is calculated, and this median is subtracted from the corresponding row in the satellite image. For example, if the median for row 5 of the high-frequency image is X, each pixel in row 5 of the satellite image will have a value equal to its original value minus X.

[0163] In a particular embodiment, the method comprises a stage prior to stage (b) of optical distortion correction in the plurality of satellite images, the optical distortion correction stage being carried out by means of the following sub-stages:

[0164] - receiving a distortion profile of the sensor that comprises: o the position of the sensor pixel that acts as the optical distortion center, and or an optical distortion function that indicates the optical distortion value as a function of the distance to the distortion center;

[0165] - For each pixel of the sensor, calculate the distance from the pixel to the center of optical distortion, the distance measured according to a fourth pre-established standard, preferably the Euclidean standard;

[0166] - For each pixel of the sensor, calculate the optical distortion value using the optical distortion function and the distance calculated in the previous sub-stage; where the optical distortion value indicates the displacement that would have to be made in the pixel to obtain a fictitious image without optical distortion;

[0167] - generate a coordinate map that defines the relationship between the coordinates of the sensor pixels and the coordinates of the pixels of the fictitious image without optical distortion, where the coordinates of the pixels of the fictitious image without optical distortion are calculated by adding to the position of the sensor pixels their corresponding optical distortion values ​​obtained in the previous sub-stage and where the coordinates of the pixels of the satellite images coincide with the coordinates of the sensor pixels;

[0168] - obtain a plurality of corrected images by applying the coordinate map to the plurality of satellite images, where the coordinates of the pixels of the corrected images coincide with the coordinates of the pixels of the fictitious image without optical distortion, and the value of each pixel of a corrected image is obtained by interpolating the pixel values ​​of a satellite image whose coordinates correspond to the pixel coordinates of the corrected image according to the coordinate map.

[0169] Optical distortion is the variation in magnification across the image plane relative to its center. This phenomenon occurs in the sensors aboard the satellite, where the resulting satellite image is distorted or deviates from its original shape due to various irregularities, such as aberrations and blurring.

[0170] In particular, optical aberrations are deviations from the actual image that occur due to imperfections in the lenses or other optical elements of the sensor. These can include chromatic aberrations (differences in how lenses focus different wavelengths of light), spherical aberrations (caused by the inability of lenses to focus light evenly across all areas of the image), and other forms of optical aberrations that affect image quality.

[0171] Most sensors exhibit pincushion distortion, where image magnification increases with distance from the optical axis. The visible effect is that lines not passing through the center of the image curve inward, toward the center.

[0172] This optical distortion manifests itself after the fusion stage, since the plurality of common regions does not coincide completely, leading to the appearance of artifacts.

[0173] To correct optical distortion, a distortion profile must first be generated by calibrating the sensor. This distortion profile defines the position of the sensor pixel that acts as the center of optical distortion, which may or may not coincide with the physical center of the sensor. The distortion profile also includes an optical distortion function that relates the optical distortion value to the distance from the center of distortion.

[0174] Once the distortion profile is received, the distance of each pixel of the sensor to the optical distortion center is calculated, the distance measured according to a pre-established standard; preferably the Euclidean standard.

[0175] Based on the estimated distances and using the optical distortion function, an optical distortion value is calculated for each pixel of the sensor. This optical distortion value indicates the displacement that would have to be made to the pixel to obtain a fictitious image without optical distortion.

[0176] Throughout this document, a fictitious image without optical distortion will be understood as a matrix of a predetermined size, where each element of the matrix corresponds to a pixel of predetermined value; for example, zero or one.

[0177] Next, a coordinate map is generated that defines the relationship between the coordinates of the sensor pixels (which coincide with the coordinates of the satellite images acquired by that sensor and which exhibit optical distortion) and the coordinates of the pixels in the fictitious image without optical distortion. The coordinates of the pixels in the fictitious image without optical distortion are calculated by adding the corresponding optical distortion values ​​obtained in the previous sub-stage to the position of the sensor pixels.

[0178] Finally, the coordinate map is applied to the satellite images to obtain a plurality of corrected satellite images. The pixel coordinates of these corrected satellite images coincide with the pixel coordinates of the optically undistorted dummy image, and the sensor pixel coordinates coincide with the pixel coordinates of the satellite images. The coordinate map to be applied to the satellite images indicates which pixel or pixels of the distorted image (satellite image acquired by the sensor) should be used to calculate the value of each pixel of the corrected image (without optical distortion) by interpolation.

[0179] The coordinate map can also be applied to previously corrected images with no fixed and / or variable noise, and / or to images that have undergone flat correction. All of these images exhibit optical distortion.

[0180] In a particular embodiment, the method comprises a stage prior to stage (b) of geometric distortion correction in the plurality of satellite images, the geometric distortion correction stage being carried out by means of the following sub-stages:

[0181] - divide the plurality of satellite images into at least one set of satellite images, each of the sets of satellite images comprising a plurality of consecutive satellite images;

[0182] - Select a representative image for each set of satellite images; and for each representative image of each set of satellite images:

[0183] - estimate a representative time instant at which the sensor acquired the representative image from the telemetry data;

[0184] - estimate the position and attitude of the satellite at the estimated representative time instant from the telemetry data;

[0185] - calculate the pointing vectors of the sensor pixels with respect to the inertial reference system centered on a celestial body from the position and attitude of the satellite at the representative time instant;

[0186] - obtain the projection of each pixel of the sensor onto the surface of the celestial body by calculating the intersection of the pointing vectors of the sensor pixels with a reference geometric figure of a representative geodetic model of the celestial body;

[0187] - estimate the effective height of each pixel of the sensor by calculating the magnitude of its corresponding pointing vector, where the terminal point of each pointing vector is its intersection with the surface of the celestial body;

[0188] - Calculate the ground sampling distance, GSD, of each sensor pixel, which represents the size of the sensor pixel's projection onto the surface of the celestial body: sensor pixel size

[0189] GSD = effective height ■ - - - - — — - sensor focal length

[0190] - select the minimum GSD from among all the GSDs calculated for all the pixels of the sensor;

[0191] - For each pixel of the sensor, divide the distance between the projection of the pixel and the projection of the pixel from the center of the sensor by the minimum GSD, thus obtaining an equivalent number of pixels for that distance;

[0192] - generate a fictitious image without geometric distortion from the equivalent number of pixels calculated in the previous sub-stage for each pixel of the sensor;

[0193] - generate a coordinate map that defines the relationship between the coordinates of the sensor pixels and the coordinates of the pixels of the fictitious image without geometric distortion, where the coordinates of the satellite image pixels coincide with the coordinates of the sensor pixels;

[0194] - obtain a plurality of corrected images by applying the coordinate map to the set of satellite images, where the coordinates of the pixels of the corrected images coincide with the coordinates of the pixels of the fictitious image without geometric distortion, and the value of each pixel of a corrected image is obtained by an interpolation of the pixel values ​​of a satellite image whose coordinates correspond to the coordinates of the pixel of the corrected image according to the coordinate map, the satellite image being part of the set of satellite images to which the representative image belongs.

[0195] The geometric distortion present in satellite images of a celestial body occurs due to the sensor's position relative to the body's surface at the time of acquisition. These perspective variations during image acquisition pose a problem during image alignment and recording.

[0196] In images taken from a very low angle (that is, those in which parts of the surface being photographed are not directly below the sensor, but rather offset at a certain angle), the actual distance between two objects located in adjacent pixels varies considerably. Therefore, distortion correction geometry involves creating a new image in which the distance represented by each pixel is constant.

[0197] In this realization, the correction stage involves correcting the geometric distortion of the raw satellite images acquired by the sensor from a coordinate map. This correction can be applied to previously corrected images with no fixed and / or variable noise, to images that have undergone flat correction, and / or to images that have undergone optical distortion correction.

[0198] To expedite the distortion geometry correction process, the multiple satellite images are divided into one or more sets, for each of which a representative image is selected. A coordinate map is then generated for each of these sets, and this coordinate map is applied to all satellite images belonging to its corresponding set to correct the distortion geometry.

[0199] First, it is necessary to know, from the telemetry data, the exact time at which the satellite's onboard sensor acquired the representative satellite image of a set with distorted geometries, as well as the satellite's position and attitude relative to the inertial reference frame at that time. With this data, the off-nadir at which each representative image of each set was taken can be calculated.

[0200] In a realization, to estimate this representative time instant, as well as the satellite's position and attitude, telemetry data is interpolated. As mentioned previously, the image reception speed is much higher than the telemetry data reception speed. For example, cameras can capture 50 images per second, while the satellite's position and attitude data, along with the acquisition time instant, are received at an average frequency of 4 Hz, which can reach at most 10 Hz. Therefore, for each representative image in a set, the satellite's position and attitude, as well as the acquisition time instant, can be estimated by interpolating the received data.

[0201] Next, the pointing vectors of all the sensor pixels are calculated with respect to the inertial reference system, taking into account the position and attitude of the satellite at the representative time instant.

[0202] The next step is to obtain the projection of each sensor pixel onto the surface of the celestial body by calculating the intersection of each pointing vector with a reference geometric figure from a geodetic model representative of the celestial body. For example, the WGS 84 reference ellipsoid if the celestial body is Earth, this ellipsoid belonging to a known geodetic model.

[0203] The method continues by estimating the effective height of each pixel of the sensor by calculating the magnitude of its corresponding pointing vector, understanding that the terminal point of each vector is its intersection with the surface of the celestial body.

[0204] The pointing vectors of each pixel, and therefore the effective heights, are variable. The further the pixel's pointing vector is from nadir, the more variable it becomes.

[0205] Once the effective height of each pixel is known, its ground sampling distance (GSD) is calculated. This GSD represents the size of the sensor pixel's projection onto the surface of the celestial body.

[0206] GSD = effective height ■ - - - - — — - sensor focal length

[0207] The sensor pixel size and sensor focal length are known commercial parameters, which depend on the type of sensor analyzed.

[0208] For each projected pixel, the distance between the projection of that pixel and the projection of the sensor's central pixel is then calculated as the difference between the estimated coordinates for both pixel projections. These coordinates were estimated by the intersection of each pointing vector with the reference geometric figure of the geodetic model representing the celestial body.

[0209] The calculated distance is then divided by the minimum GSD calculated for all projected sensor pixels. This yields an equivalent number of pixels for each distance, where all pixels represent the same distance. From the estimated equivalent number of pixels for each distance, a simulated image without geometric distortion is generated. This image will contain more pixels than the actual sensor and the satellite images, especially as the nadir of the pointing decreases, because the geometric distortion will be greater.

[0210] Throughout the document, a fictitious image without geometric distortion will be understood as a matrix of a size in accordance with the calculated equivalent number of pixels, where each element of the matrix corresponds to a pixel of predetermined value; for example, zero or one.

[0211] Finally, from this simulated image without geometric distortion, a corrected image of the same size is obtained by interpolating the pixel values ​​of the satellite image. Specifically, the coordinate map indicates which pixel or pixels of the satellite image (the image acquired by the sensor that exhibits geometric distortion) should be used to calculate the value of each pixel in the interpolated corrected image. The coordinates of the pixels in the corrected image coincide with the coordinates of the pixels in the simulated image without geometric distortion, and the coordinates of the sensor pixels coincide with the coordinates of the pixels in the satellite images.

[0212] The coordinate map can also be applied to previously corrected images with no fixed and / or variable noise, and / or to images that have undergone flat correction, and / or to images corrected for optical distortion. All of these images exhibit geometric distortion.

[0213] In one embodiment, the method comprises a stage prior to stage (b) of correcting optical and geometric distortions in the plurality of satellite images, the stage of correcting optical and geometric distortions being carried out by means of the following sub-stages:

[0214] - receive a distortion profile from the sensor comprising: either the position of the sensor pixel that acts as the center of optical distortion, or an optical distortion function that indicates the value of optical distortion as a function of the distance to the center of distortion;

[0215] - For each pixel of the sensor, calculate the distance from the pixel to the center of optical distortion, the distance measured according to a fourth pre-established standard, preferably the Euclidean standard;

[0216] - For each pixel of the sensor, calculate the optical distortion value using the optical distortion function and the distance calculated in the previous sub-stage; where the optical distortion value indicates the displacement that would have to be made in the pixel to obtain a fictitious image without optical distortion;

[0217] - generate a first coordinate map that defines the relationship between the coordinates of the sensor pixels and the coordinates of the pixels of the fictitious image without optical distortion, where the coordinates of the pixels of the fictitious image without optical distortion are calculated by adding to the coordinates of the sensor pixels their corresponding optical distortion values ​​obtained in the previous sub-stage and where the coordinates of the pixels of the satellite images coincide with the coordinates of the sensor pixels;

[0218] - divide the plurality of satellite images into at least one set of satellite images, each of the sets of satellite images comprising a plurality of consecutive satellite images;

[0219] - Select a representative image for each set of satellite images; and for each representative image of each set of satellite images:

[0220] - estimate a representative time instant at which the sensor acquired the representative image from the telemetry data;

[0221] - estimate the position and attitude of the satellite at the estimated representative time instant from the telemetry data;

[0222] - calculate the pointing vectors of the sensor pixels with respect to the inertial reference system centered on the celestial body from the position and attitude of the satellite at the representative time instant;

[0223] - obtain the projection of each pixel of the sensor onto the surface of the celestial body by calculating the intersection of the pointing vectors of the sensor pixels with a reference geometric figure of a representative geodetic model of the celestial body;

[0224] - estimate the effective height of each pixel of the sensor by calculating the magnitude of its corresponding pointing vector, where the terminal point of each pointing vector is its intersection with the surface of the celestial body;

[0225] - Calculate the ground sampling distance, GSD, of each sensor pixel, which represents the size of the sensor pixel's projection onto the surface of the celestial body: sensor pixel size

[0226] GSD = effective height ■ - - - - — — - sensor focal length

[0227] - select the minimum GSD from among all the GSDs calculated for all the pixels of the sensor;

[0228] - For each pixel of the sensor, divide the distance between the projection of the pixel and the projection of the pixel from the center of the sensor by the minimum GSD, thus obtaining an equivalent number of pixels for that distance;

[0229] - generate a fictitious image without optical and geometric distortions from the equivalent number of pixels calculated in the previous sub-stage for each pixel of the sensor;

[0230] - generate a second coordinate map that defines the relationship between the pixel coordinates of a fictitious image without optical distortion and with geometric distortion, and the pixel coordinates of the fictitious image without optical and geometric distortion generated; where the pixel coordinates of the fictitious image without optical distortion and with geometric distortion coincide with the pixel coordinates of the sensor and with the pixel coordinates of the satellite images, and it is assumed that said image is free of optical distortion;

[0231] - generate a combined coordinate map, which defines the relation between the coordinates of the pixels of the fictitious image without optical and geometric distortions, and the coordinates of the sensor pixels, where:

[0232] (i) Evaluating the second coordinate map, the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion are identified that correspond to the coordinates of the pixels of the fictitious image without optical and geometric distortion;

[0233] (ii) Evaluating the first coordinate map, the coordinates of the pixels of the fictitious image without optical distortion are estimated by interpolation, which correspond to the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion identified in step (i),

[0234] (iii) Evaluating the first coordinate map, the coordinates of the sensor pixels that correspond to the coordinates of the optically undistorted fictitious image pixels estimated in step (ii) are determined by interpolation; and

[0235] (iv) the combined coordinate map is generated by establishing the relation between the sensor pixel coordinates determined in step (iii) and the pixel coordinates of the fictitious image without optical and geometric distortions; and

[0236] - obtain a plurality of corrected images by applying the combined coordinate map to the set of satellite images, where the coordinates of the pixels of the corrected images coincide with the coordinates of the pixels of the fictitious image without optical and geometric distortions, and the value of each pixel of a corrected image is obtained by combining the values ​​of the pixels of a satellite image whose coordinates correspond to the coordinates of the pixel of the corrected image according to the combined coordinate map, the satellite image being part of the set of satellite images to which the representative image belongs.

[0237] In this realization, the method stages that allow for the correction of optical and geometric distortions at the same time are described.

[0238] First, an initial coordinate map is generated following the steps described in the implementation that outlines the optical distortion correction in isolation. Specifically, this initial coordinate map indicates which pixel or pixels in the optically distorted image (the satellite image acquired by the sensor) should be used to calculate the value of each pixel in an image without optical distortion through interpolation.

[0239] In this case, the first coordinate map is not applied to the satellite images, but is used to subsequently generate a combined coordinate map.

[0240] Once the first coordinate map is obtained, the method continues with the realization steps in which the geometric distortion correction is defined to generate a second coordinate map.

[0241] Specifically, this second coordinate map indicates which pixel or pixels of a fictitious image with geometric distortion, but without optical distortion, we should take values ​​from to calculate the value of each pixel of a fictitious image without optical and geometric distortion by interpolation. For the generation of the second coordinate map, it is assumed that the images with geometric distortion have already been corrected for optical distortion. Throughout this document, a fictitious image without optical distortion but with geometric distortion will be understood as a matrix of a size corresponding to the calculated equivalent number of pixels, where each element of the matrix corresponds to a pixel with a predetermined value; for example, zero or one.

[0242] Throughout the document, a fictitious image without optical and geometric distortions will be understood as a matrix of a size in accordance with the calculated equivalent number of pixels, where each element of the matrix corresponds to a pixel of predetermined value; for example, zero or one.

[0243] The method continues by generating a combined coordinate map from the first coordinate map and the second coordinate map.

[0244] To generate this combined coordinate map, the following steps are performed.

[0245] (i) Evaluating the second coordinate map, the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion are identified that correspond to the coordinates of the pixels of the fictitious image without optical and geometric distortion.

[0246] (ii) Evaluating the first coordinate map, the coordinates of the pixels of the fictitious image without optical distortion (from the first coordinate map) that correspond to the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion (from the second coordinate map) identified in step (i) are estimated by interpolation.

[0247] (iii) Evaluating the first coordinate map, the coordinates of the sensor pixels that correspond to the coordinates of the optically undistorted fictitious image pixels estimated in step (ii) are determined by interpolation.

[0248] (iv) Finally, the combined coordinate map is generated by establishing the relation between the sensor pixel coordinates determined in step (iii) (from the first coordinate map) and the pixel coordinates of the fictitious image without optical and geometric distortions of the second coordinate map.

[0249] The distortion correction process concludes by obtaining a plurality of corrected images of the same size as the image without optical and geometric distortion. To achieve this, the combined coordinate map is applied to each satellite image, so that the value of each pixel in a corrected image (whose coordinates coincide with the coordinates of the pixels in the image without optical and geometric distortion) is obtained by combining the pixel values ​​of a satellite image (with coordinates coinciding with those of the sensor pixels) whose positions correspond to the pixel coordinates of the corrected image (or to the pixel coordinates of the fictitious image without optical and geometric distortion) according to the combined coordinate map.

[0250] Alternatively, the combined coordinate map can also be applied to previously corrected images without fixed and / or variable noise, and / or to images that have had a flat correction applied.

[0251] To expedite the correction of optical and geometric distortions, the plurality of satellite images is divided into one or more sets, for each of which a representative image is selected. Thus, the second coordinate map and the combined coordinate map are generated for each of these sets, with the combined coordinate map being applied to all satellite images belonging to the set for which it was generated.

[0252] Advantageously, creating a combined coordinate map only requires correcting the satellite images once, meaning that the interpolation of pixel values ​​from the satellite images occurs only once. This speeds up the calculations required for correction, since the optical and geometric distortion correction stage is faster than performing the optical and geometric distortion correction stages separately. Additionally, using a single combined coordinate map reduces system memory requirements.

[0253] In a particular embodiment:

[0254] - the estimation of the first time instant in which the sensor acquired the first satellite image is carried out by interpolating the plurality of telemetry time instants; and / or

[0255] - the estimation of the second time instant at which the sensor acquired the second satellite image is carried out by interpolating the plurality of telemetry time instants; and / or

[0256] - the estimation of the representative time instant at which the sensor acquired the representative satellite image is carried out by interpolating the plurality of telemetry time instants; and / or

[0257] - The estimation of the position and attitude of the satellite at the first time instant, and / or at the second time instant and / or at the estimated representative time instant is carried out by interpolating the plurality of positions and attitudes of the satellite from the telemetry data.

[0258] The speed at which images are received is generally much higher than the speed at which telemetry data is received. For example, the sensor can capture 50 images per second, while the satellite's position and attitude data, along with the acquisition time, are received at an average frequency of 4 Hz, which can reach at most 10 Hz.

[0259] Thus, for one and / or both images to be aligned, in this realization the position and attitude of the satellite are estimated, as well as the time of acquisition, by interpolating the telemetry data received.

[0260] Additionally or alternatively, for the representative image of each set of satellite images (according to the realizations in which the geometric distortion or the optical and geometric distortions are corrected) in this realization the position and attitude of the satellite is estimated, as well as the representative time instant of acquisition, by interpolating the telemetry data received.

[0261] A second inventive aspect provides a system for processing satellite images comprising a processor configured to carry out steps (e.2) and (e.3) of the method of the first inventive aspect and a programmable logic gate array, FPGA, configured to carry out the remaining steps of the method of the first inventive aspect.

[0262] Throughout this document, an FPGA (Field Programmable Gate Array) will be understood as a programmable device comprising logic blocks configurable in terms of their interconnection and functionality. The FPGA allows the real-time design and implementation of custom digital circuits, offering great flexibility to adapt to specific needs without having to change the physical chip. Its reprogrammability is particularly advantageous in environments requiring continuous hardware adaptability.

[0263] In one realization, the processor and the FPGA are configured to execute steps (a)-(e) of the method on board a satellite.

[0264] In one realization, the FPGA is in communication with the satellite sensor to receive the raw satellite images acquired by said satellite.

[0265] In one realization, the FPGA is additionally configured to deserialize the raw satellite images to translate them into pixels.

[0266] In one realization, the FPGA is additionally configured to compress satellite images.

[0267] In one realization, the system additionally comprises a storage element and the FPGA is configured to store the satellite images, the compressed satellite images, the reference image and / or the set of common regions in the storage element.

[0268] In one realization, the FPGA comprises at least one direct memory access module, DMA, configured to access the pair of images to be aligned and / or the reference image and / or the set of common regions previously stored in the system storage element.

[0269] In one realization, the FPGA comprises at least one write-only direct memory access (DMA) module configured to store, in the storage element, the offset calculated by the FPGA during the alignment of a pair of images.

[0270] In this implementation, the storage element is memory, a disk, or a database that communicates with the FPGA. In this implementation, the storage element is DDR memory (double data rate memory).

[0271] In a realization in which the displacement is determined by minimizing a linear cross-correlation between the first satellite image and the second satellite image to be aligned, the FPGA comprises a cross-correlation module configured to perform a linear cross-correlation between the images.

[0272] In a more particular realization, the cross-correlation module comprises a plurality of modules configured to perform the operations required by cross-correlation; for example, a module that applies a Vonann matrix, a module that calculates the Fourier transform in two dimensions, or a complex number multiplier module.

[0273] In a realization, the FPGA is configured to adjust the precision of the integer and fractional parts of the result of each operation it performs.

[0274] The complexity of the FPGA lies in its inability to perform floating-point operations. Therefore, for each operation, the FPGA, according to this implementation, adjusts the precision of the integer and fractional parts of the result to prevent overflow. Throughout this document, an overflow will be understood to occur when an attempt is made to create a numeric value that falls outside the range that can be represented with a predetermined number of digits, either greater than the maximum or less than the minimum representable value.

[0275] In one embodiment, the FPGA comprises a demultiplexer, a plurality of image and weight storage modules, and a fusion module, wherein: the demultiplexer is configured to receive, by means of at least one DMA memory access module, the reference image and the set of common regions together with the weights estimated in sub-step (e.3) of the method by the system processor; each image and weight storage module is configured to store and send the reference image or a common region from the set of common regions together with their weights to the fusion module; and the fusion module is configured to receive the weights and receive and fuse the reference image and the set of common regions based on those weights.

[0276] In a realization, each image and weight storage module is configured to store a plurality of rows from the reference image or common region, and to send part of the plurality of rows from the reference image or common region to the merging module while temporarily storing the remaining rows.

[0277] In this way, the FPGA advantageously stores and processes in parallel the lines of the reference image and the set of common regions.

[0278] In this implementation, each image and weight storage module is configured as a ping-pong buffer. Advantageously, this configuration allows storing part of the rows of the reference image or a common region while sending another part of the rows, thus avoiding idle time that unnecessarily lengthens the method's execution time.

[0279] In a realization, the image storage and weight modules are configured to send to the fusion module at the same time.

[0280] In one embodiment, the FPGA further comprises an output memory comprising a plurality of RAM memory blocks, e.g., BRAMs, configured to store a plurality of lines of the oversampled fused image and to send such lines to the processor; e.g., to store the oversampled fused image in a storage element or to send it to an additional system module or to a receiver located on Earth.

[0281] In one realization, the FPGA additionally comprises a special width converter configured to change the width of the data; for example, to go from 48 bits input to 64 bits output.

[0282] In one realization, the FPGA comprises at least one write-only direct memory access (DMA) module configured to send the oversampled fused image to a processor.

[0283] In a realization, the number of pixels in the reference image and the set of common regions that contribute to each pixel in the oversampled merged image is configurable.

[0284] In a realization, in the same clock cycle, a plurality of oversampled pixels of the fused image are calculated in parallel, for example, four pixels.

[0285] In one particular implementation, the FPGA is additionally configured to perform, in parallel:

[0286] - at least one of the sub-stages of stage (b) of the method; and / or

[0287] - stage (c) of the method; and / or

[0288] - stage (d) of the method; and / or

[0289] - at least one of the sub-stages of stage (e) of the method.

[0290] Image alignment and image fusion are two computationally expensive processes. Therefore, to process satellite images in real time or near real time, beyond software acceleration techniques, hardware techniques that allow for high parallelization capabilities, such as FPGA technology, are required.

[0291] An FPGA allows different processes to be executed in parallel; in this way, in this realization, different stages and sub-stages of the method are executed in the FPGA in parallel to take advantage of its ability to accelerate computationally expensive operations, thus reducing the total processing time.

[0292] For example, the FPGA can parallelize the sub-stages of image fusion, in which multiplication and addition operations are performed on a large number of pixels, which represents a significant saving in the execution time of the fusion stage.

[0293] Another possible hardware solution to parallelize the execution of the method is the use of a GPU (Graphics Processing Unit), an electronic device that allows high-speed mathematical calculations or graphics processing.

[0294] However, the FPGA-based solution has lower power consumption than the GPU-based one, so the use of FPGAs is optimal when running the image processing method on board a satellite.

[0295] In a particular embodiment, the processor comprises at least one of the following threads: a first thread configured to carry out the fixed noise correction steps, and / or the variable noise correction steps, and / or the flat correction step of the method of the first inventive aspect; or

[0296] - a second thread configured to carry out the optical distortion correction stage of the method of the first inventive aspect or the optical and geometric distortion correction stage of the first inventive aspect; or

[0297] - a third thread configured to carry out the geometric distortion correction stage of the method of the first inventive aspect or the optical and geometric distortion correction stage of the first inventive aspect; or

[0298] - a fourth thread to command the FPGA to execute steps (a)-(e) of the method of the first inventive aspect; or

[0299] - a combination of any of the above; and where the processor threads are additionally configured to carry out the stages in parallel.

[0300] In this implementation, in addition to the FPGA, the system includes a multi-threaded processor. Each of these threads executes one or more stages of image correction, or commands the actions of the FPGA.

[0301] Specifically, the satellite images acquired by the sensors are processed by different processor threads to execute one or more image correction stages in parallel; for example, to remove fixed and / or variable noise, or to correct optical and / or geometric distortions. In addition, the processor has another thread to command the FPGA to execute stages (a)-(e) of the method in parallel with the execution of the correction stages.

[0302] In one embodiment, the processor utilizes both the second and third threads to perform the optical distortion correction and geometric distortion correction stages separately, i.e., each correction stage on a different thread. In an alternative embodiment, the processor utilizes either the second or third thread to perform the optical and geometric distortion correction stages simultaneously.

[0303] In a realization, the actions of the different processor threads are executed on board a satellite and in real time.

[0304] In one realization, the processor is additionally configured to manage radio data transmission, to send telemetry and / or to monitor temperatures.

[0305] In one implementation, the FPGA and the processor are part of a Multi-Processor System On-Chip (MPSoC).

[0306] Throughout this document, an MPSoC is defined as a multi-microprocessor chip that includes at least one microprocessor and one FPGA, which acts as the processing core. Optionally, the MPSoC includes a microcontroller.

[0307] In a realization, the MPSoCs is part of a System-in-Module (SOM) along with memory, and the communications and electronics necessary for its operation.

[0308] In one particular embodiment, the SOM is part of an on-board computer (OnBoard Computer, OBC).

[0309] In one realization, the MPSoC's microprocessor and / or FPGA are additionally configured to transmit the super-resolution images obtained by the FPGA to a receiver located on Earth.

[0310] In one particular embodiment, the system is implemented on a printed circuit board, PCB.

[0311] In one embodiment, the PCB comprises at least one communication interface. In one embodiment, the PCB is further configured to monitor electronic, physical, and software parameters, such as temperature, power supply, voltages, currents, and internal system clocks. In one embodiment, the PCB comprises a microcontroller configured to manage the power electronics and system redundancies.

[0312] A third inventive aspect provides a celestial body observation camera comprising:

[0313] - a computer comprising at least one system according to the second inventive aspect; and

[0314] - a sensor configured to acquire a plurality of satellite images. In one implementation, the camera's computer comprises a plurality of systems. Advantageously, the computer has redundancy, replicating the system two or more times for greater reliability.

[0315] In one implementation, the camera computer comprises a plurality of systems implemented on a PCB. Advantageously, the computer comprises a redundant PCB, replicated two or more times for greater reliability.

[0316] In one realization, the PCB board comprises at least two sockets to integrate two or more systems, which facilitates the maintainability and scalability of the camera.

[0317] In one implementation, the computer is an OBC comprising a PCB with two SOMs (System-on-Module). Each SOM includes a MPSoC and the memory, communications, and electronics necessary for its operation. Advantageously, the SOMs allow for modularity and versatility through standard connectors, facilitating interchangeability between boards.

[0318] A fourth inventive aspect provides a satellite comprising the celestial body observation camera of the third inventive aspect.

[0319] The payload of the satellite of the fourth inventive aspect is the celestial body observation camera, comprising the sensor configured to acquire satellite images, and a computer (OBC) comprising one or more systems according to the second inventive aspect.

[0320] The process of acquiring and processing satellite images begins with instructions from the ground regarding the geographic area to be photographed. Once the satellite reaches the specified location, the On-Board Control (OBC) is activated, orienting the satellite towards the area of ​​interest to begin image acquisition using a predetermined communication protocol, such as Cameralink or Aurora.

[0321] The images captured by the sensor are transmitted to the OBC, where the system of the second inventive aspect plays a crucial role. This system comprises an FPGA configured to receive and read the satellite images, deserialize them (translating them into pixels), and process them using the method of the first inventive aspect. Optionally, the processing carried out by the FPGA includes image compression and / or storage in a storage element (e.g., in a memory communicating with the FPGA).

[0322] Once the desired area has been captured, the OBC system shuts down. For image download, when the satellite flies over a reception area, the stored images are uploaded from the storage device and transmitted via radio.

[0323] The execution of the method of the first inventive aspect can be carried out either immediately after the acquisition of images, or immediately before their transmission.

[0324] Advantageously, the satellite of the invention facilitates image processing within the required timeframes and coordinates internal and external communications. This optimizes the acquisition, processing, and transmission of images from space, significantly improving the performance and quality of the captured images, with options both for maintaining images within the satellite's electronics until transmission and for direct transmission via compatible devices.

[0325] All features and / or method steps described in this specification (including the claims, description and drawings) may be combined in any combination, except for combinations of such features that are mutually exclusive.

[0326] DESCRIPTION OF THE DRAWINGS

[0327] These and other features and advantages of the invention will become clearer from the following detailed description of a preferred embodiment, given only as an illustrative and non-limiting example, with reference to the accompanying figures.

[0328] Figure 1 This figure shows a schematic of the computer-implemented method for processing satellite images according to a realization of the invention.

[0329] Figure 2 This figure illustrates an example of the realization of the invention in which the initial horizontal and vertical displacements between a pair of satellite images are calculated based on telemetry data.

[0330] Figure 3 This figure shows the alignment of a pair of satellite images based on image pyramids with different resolutions according to a realization of the invention.

[0331] Figures 4a-4b These figures show the generation of a pyramid of satellite images (Fig. 4a) as a function of different sub-sampling factors (Fig. 4b).

[0332] Figure 5 This figure shows the alignment of two pairs of satellite images based on image pyramids with different resolutions, the pyramids having a different number of levels, according to a realization of the invention.

[0333] Figures 6a-6c These figures show the effects of optical distortion and geometric distortion on satellite images.

[0334] Figure 7 This image shows a system for processing satellite images according to an embodiment of the invention.

[0335] Figure 8 This image shows a schematic of the architecture of an FPGA system to carry out the alignment of a pair of satellite images according to a realization of the invention.

[0336] Figure 9 This image shows a schematic of the architecture of an FPGA of the system to carry out the obtaining of the oversampled fused image according to a realization of the invention.

[0337] Figure 10 This image shows an observation camera comprising a plurality of systems for processing satellite images according to an embodiment of the invention.

[0338] Figure 11 This image shows a satellite comprising a plurality of systems for processing satellite images according to an embodiment of the invention.

[0339] DETAILED EXPLANATION OF THE INVENTION

[0340] Figure 1 shows a schematic of the computer-implemented method (100) for processing satellite images (1) according to an embodiment of the invention, the method (100) comprising the following steps:

[0341] Stage a

[0342] In step (a), a plurality of satellite images (1) acquired by a sensor (2) on board a satellite (3) are received. Additionally, telemetry data (4) are received, this telemetry data (4) comprising a plurality of positions and attitudes of the satellite (3) with respect to an inertial reference frame centered on a celestial body (14) at a plurality of telemetry time points. The sensor (2) on board the satellite (3) acquires the plurality of satellite images (1) during the course of the plurality of telemetry time points.

[0343] The number of satellite images (1) shown in this Figure 1 is three, but this number may be higher or lower.

[0344] Stage b

[0345] During stage (b), a first pair of consecutive satellite images from the plurality of satellite images (1) is aligned by carrying out the following sub-stages.

[0346] In sub-stage (b.1), a first satellite image (1.1), named as Im1 in Figure 1, and a second satellite image (1.2), named as Im2 in Figure 1, are selected. Both images are intended to be aligned.

[0347] In sub-stage (b.2), an initial horizontal displacement (8.1) and an initial vertical displacement (8.2) of the second satellite image (1.2) with respect to the first satellite image (1.1) are calculated. To do this: a first time instant is estimated, at which the sensor (2) acquired the first satellite image (1.1), from the telemetry data (4); the position and attitude of the satellite (3) at the first estimated time instant are estimated from the telemetry data (4); a second time instant is estimated, at which the sensor (2) acquired the second satellite image (1.2), from the telemetry data (4); the position and attitude of the satellite (3) at the second estimated time instant are estimated from the telemetry data (4); and the initial displacements (8.1, 8.2) are calculated from the positions and attitudes of the satellite (3) at the first time instant and at the second time instant.

[0348] In step (b.3) the initial displacements (8.1, 8.2) are refined; that is, a total horizontal displacement (9.1) and a total vertical displacement (9.2) are calculated, which must be applied to the second satellite image (1.2) to minimize the distance between this second satellite image (1.2) and the first satellite image (1.1). The distance must be measured according to a pre-established standard, preferably the Euclidean standard.

[0349] In step (b.4) the first satellite image (1.1) and the second satellite image (1.2) are aligned by applying the total displacements (9.1, 9.2) to the second satellite image (1.2).

[0350] Stage c

[0351] Once the first pair of images Im1 and Im2 has been aligned, method (100) continues in step (c) by repeating step (b) for all consecutive pairs of satellite images from the plurality of satellite images (1).

[0352] In each iteration of stage (b), the first satellite image (1.1) selected in sub-stage (b.1) corresponds to the second satellite image (1.2) of the previous iteration. This is shown in Figure 1, where image Im2 from the previous iteration becomes image Im1 in stage (c).

[0353] Furthermore, in each iteration of step (b), the total horizontal displacement (9.1) obtained in sub-step (b.3) must be added to the total horizontal displacement (9.1) obtained in sub-step (b.3) of the previous iteration, and the total vertical displacement (9.2) obtained in sub-step (b.3) must be added to the total vertical displacement (9.2) obtained in sub-step (b.3) of the previous iteration. This is necessary to ensure that all images (1) are aligned with respect to the same initial image, in this case, the image Im1 selected in the first iteration of step (b).

[0354] Stage d

[0355] Next, the method (100) continues in stage (d), selecting as the reference image (10) the first satellite image (1.1) selected in sub-stage (b.1) of the first iteration of stage (b).

[0356] Additionally, at least one common region is obtained between the reference image (10) and each of the second satellite images (1.2) that will have been shifted in sub-step (b.4) of all iterations of step (b). The result is a set of common regions (11). Method (100) ends with step (e), in which a fused image (13) with oversampling, that is, super-resolution, is obtained. To achieve this, the reference image (10) and the set of common regions (11) are fused.

[0357] This merger consists of several sub-stages:

[0358] First, in sub-step (e.1), an oversampled image (12) is generated, where the resolution of the oversampled image (12) is the resolution of the reference image (10) multiplied by a predefined oversampling factor. The pixel values ​​of this oversampled image (12), for example, are zeros or ones.

[0359] In sub-step (e.2), for each pixel (12.1) of the oversampled image (12), the distance from the center of that pixel (12.1) to the center of one or more pixels of each region of the set of regions (11) and to the center of one or more pixels of the reference image (10) is determined. The distance is measured according to a second pre-established norm, preferably the Euclidean norm.

[0360] In sub-step (e.3), for each pixel (12.1) of the oversampled image (12.1), a weight (p1, p2, and p3 in Figure 1) is estimated for the pixel(s) in each region of the set of regions (11) and for the pixel(s) in the reference image (10) based on the determined distances. In one example, the greater the estimated distance, the lower the weight assigned.

[0361] Finally, in sub-step (e.4) the oversampled fused image (13) is obtained. The value of each pixel (13.1) of the oversampled fused image (13) is calculated by summing the values ​​of the pixel(s) of each region of the set of regions (11) and of the pixel(s) of the reference image (10), weighting them according to the weights (p1, p2, p3) estimated in the previous sub-step.

[0362] Figure 2 illustrates an example of the realization of the invention in which the initial horizontal (8.1) and vertical (8.2) displacements between a pair of satellite images (1.1, 1.2) are calculated as a function of telemetry data (4); in particular, from the estimation of the positions and attitudes of the satellite (3) at the first time instant and at the second time instant.

[0363] This calculation of the initial displacements (8.1, 8.2) is carried out using the following steps:

[0364] Step 1

[0365] In the first iteration of stage (b), for each pair of consecutive telemetry time instants, t teil yt tei2 The following actions are performed.

[0366] The pointing vectors of the sensor center pixel (2) and of at least two pixels from two sensor segments (2), which intersect at the sensor center (2) and are perpendicular, are calculated with respect to the inertial reference system from the position (posi and pos2) and attitude (acti and act2) of the satellite (3) at the first telemetry time instant, and at the second telemetry time instant, t tei2 In other words, for each telemetry time instant, t teil yt tei2 The positions (posi and pos2) and attitudes (acti and act2) of the satellite (3) are identified according to the telemetry data (4) and, knowing these positions (posi and pos2) and attitudes (acti and act2) of the satellite (3), the pointing vectors are calculated.

[0367] As can be seen in Figure 2, the pointing vectors have been calculated for the center pixel of the sensor (2) and for two pixels of the segments that intersect at that center, these pixels being located at the upper and right ends of the sensor (2).

[0368] Next, the projection of the sensor center pixel (2) and the pixels of the sensor segments (2) onto the surface of the celestial body (14) is obtained by calculating the intersection of the pointing vectors of said pixels with a reference geometric figure of a representative geodetic model of the celestial body (14).

[0369] In Figure 2, the celestial body (14) is the Earth and the geometric figure is the WGS 84 reference ellipsoid. However, other celestial bodies (14) and other geodetic models are also included in the context of the invention.

[0370] Method (100) further defines a projected reference system (15) as the projection of the pixels of the sensor segments (2) and the projection of the sensor center pixel (2) onto the surface of the celestial body (14), assuming that the surface of the celestial body (14) is flat, that the reference system formed by the projections of the pixels of the sensor segments (2) and the projection of the sensor center pixel (2) is orthogonal, and that the projection of the sensor center pixel

[0371] (2) corresponds to the center of the projected reference system (15)

[0372] Then, for the sensor's center pixel (2), the effective height (magnitude of its pointing vector, understanding that the vector's terminal point is its intersection with the surface of the celestial body) and the sampling distance to the ground, GSD, are calculated: sensor pixel size

[0373] GSD = effective height ■ - - - - — — - sensor focal length

[0374] Using the above data, a first linear velocity, vf, of the projection of the sensor's center pixel (2) due to the rotation of the celestial body (14) is calculated. The first linear velocity v^ is referenced in the inertial reference system and is calculated from the angular velocity, df, of the celestial body (14) and the perpendicular distance from the projection of the sensor's center pixel to the rotation axis of the celestial body. where cJ R are vectors referenced in the inertial reference system.

[0375] In this example, the angular velocity of the Earth is 7.27 x 10' 5 rad / sy R is 6.378 kilometers.

[0376] Secondly, a second linear speed, vj, of the projection of the pixel of the center of the sensor (2) due to the movement of the satellite (3) is calculated: where: o pos^ is the position vector of the projection of the pixel of the sensor center (2) at the second telemetry time instant, t tei2 , the vector being referenced in the inertial reference system; I pos is the position vector of the projection of the pixel of the center of the sensor (2) at the first time instant of telemetry, the vector being referenced in the inertial reference frame; I vj is referenced in the inertial reference frame. Finally, the total linear velocity, v^, of the projection of the sensor center pixel (2) is calculated in pixels per second and this total linear velocity is referenced in the projected reference frame (15):

[0377] To reference the total linear speed in the projected reference system (15), for example, the horizontal and vertical components of that total linear speed are projected onto the projected reference system (15).

[0378] This total linear velocity must be calculated for each pair of consecutive telemetry time points. That is, method (100) requires a first calculation of the total linear velocity taking into account telemetry time points one and two of the telemetry data (4), a second calculation of the total linear velocity taking into account telemetry time points two and three of the telemetry data (4), and so on until a final calculation of the total linear velocity taking into account the penultimate and final telemetry time points.

[0379] Step 2

[0380] In the first iteration of stage (b), the horizontal linear velocities, v, are calculated. H (t), and vertical, v v(t), as a function of time by interpolating the horizontal and vertical components of the total linear velocities, v^, obtained for all pairs of consecutive telemetry time instants.

[0381] Step 3

[0382] In any iteration of step (b), for the pair of consecutive satellite images to be aligned, the initial horizontal displacement (8.1) and the initial vertical displacement (8.2) are estimated by integrating the horizontal and vertical linear velocities with respect to time, v H (t) yv K (t), respectively, between the first time instant and the second time instant in which each satellite image of the pair of consecutive satellite images to be aligned was acquired.

[0383] Figure 3 shows the alignment of a pair of satellite images (1.1, 1.2) carried out in sub-step (b.3) of method (100), which is based on the alignment of image pyramids (16, 17) with different resolutions according to an embodiment of the invention. This particular embodiment of the satellite image alignment (1.1, 1.2) comprises the following sub-steps: Sub-step (b.3.i)

[0384] For the first satellite image (1.1), a first pyramid (16) of images is generated by subsampling the first satellite image (1.1) with N predefined subsampling factors (in this example N=3), where N is greater than or equal to one.

[0385] As can be seen in figure 3, the first image pyramid (16) comprises 4 levels (N+1 levels), where:

[0386] • the first level of the first image pyramid (16) comprises the first satellite image (1.1), and

[0387] • the i-th level of the first image pyramid (16) comprises an image subsampled from the first satellite image (1.1) with the k-th subsampling factor, where i = 2,... ,4 (that is, N+1) and k = i-1.

[0388] Sub-stage (b.3.ii)

[0389] For the second satellite image (1.2), a second image pyramid (17) is generated by subsampling the second satellite image (1.2) with M predefined subsampling factors (in this example M=3), where M is greater than or equal to one.

[0390] Figures 4a and 4b show the generation of this pyramid (17) of satellite images (Fig. 4a) as a function of the 3 different sub-sampling factors (Fig. 4b) fi, f2 and fs.

[0391] As can be seen in figure 4a, this second pyramid of images (17) comprises 4 levels (M+1 levels), where:

[0392] • M is equal to N (although it could be less),

[0393] • the first level of the second image pyramid (17) comprises the second satellite image (1.2), and

[0394] • the j-th level of the second image pyramid (17) comprises a subsampled image (1.2') from the second satellite image (1.2) with the m-th subsampling factor (fi, f2 or fs), where j = 2,... ,4 (i.e., M+1) and m = j-1.

[0395] Sub-stage (bSiii)

[0396] Starting from the initial horizontal displacement (8.1) and the initial vertical displacement (8.2), the horizontal displacement and the vertical displacement (Displacement HA / 1 in Figure 3) that must be applied to the image of level 4 (level M+1) of the second pyramid of images (17) are determined in order to minimize the distance between said image and the image of the first pyramid of images (16) of the same level 4, the distance measured according to a third pre-established norm, preferably the Euclidean norm

[0397] Sub-stage (b.3.iv)

[0398] Starting from the horizontal and vertical displacements determined in the previous sub-step (H / V Displacement 1), the horizontal and vertical displacements (H / V Displacement 2) that must be applied to the image of the lower level (level 3) of the second image pyramid (17) are determined to minimize the distance between this image and the image of the first image pyramid (16) of the same level 3, the distance measured according to the third pre-established rule. Since this lower level (level 3) is not level 1 of the second image pyramid (17), step (b.3.iv) is repeated.

[0399] Thus, starting from the horizontal and vertical displacements determined in the previous sub-step (H / V Displacement 2), the horizontal and vertical displacements (H / V Displacement 3) that must be applied to the image of the lower level (level 2) of the second image pyramid (17) are determined in order to minimize the distance between this image and the image of the first image pyramid (16) of the same level 2, the distance measured according to the third pre-established rule. Since this lower level (level 2) is not level 1 of the second image pyramid (17), step (b.3.iv) is repeated.

[0400] Thus, starting from the horizontal and vertical displacements determined in the previous sub-step (H / V Displacement 3), the horizontal and vertical displacements (H / V Displacement 4) that must be applied to the image of the lower level (level 1) of the second image pyramid (17) are determined to minimize the distance between this image and the image of the first image pyramid (16) of the same level 1, the distance measured according to the third pre-established rule. Since this lower level (level 1) is level 1 of the second image pyramid (17), the horizontal displacement of H / V Displacement 4 is established as the total horizontal displacement (9.1) and the vertical displacement of H / V Displacement 4 as the total vertical displacement (9.2). Figure 5 shows the alignments of a first pair of satellite images (1.1 , 1.2) and of a second pair of satellite images (1.1 , 1.2) based on image pyramids (16, 17) with different resolutions, the pyramids of the alignment of the first pair of images (1.1, 1.2) presenting the same number of levels and the pyramids of the alignment of the second pair of images (1.1, 1.2) presenting a different number of levels, according to an embodiment of the invention.

[0401] In particular, Figure 5 shows a first alignment of a first pair of consecutive images Im1 and Im2 during the first iteration of step (b) of method (100). In this first alignment, images 1 and 2 of the plurality of satellite images (1) are aligned, for each of which a pyramid (16, 17) with the same number of levels is generated; in this example, 4 levels (N = 3).

[0402] Additionally, Figure 5 shows a second alignment of a second pair of consecutive images Im1 and Im2 during the second iteration of step (b) of method (100). In this second alignment, images 2 and 3 of the plurality of satellite images (1) are aligned, for each of which a pyramid (16, 17) is generated with a different number of levels; in this example, 4 levels for the first pyramid (16) and 2 levels for the second pyramid (17).

[0403] It should be mentioned that the Im2 of the first alignment becomes the Im1 of the second alignment, and that in the second iteration of stage (b) of method (100), the first pyramid (16) would have already been generated in the previous iteration, so it would only be necessary to generate the second pyramid of the Im2 of the second iteration.

[0404] In one realization, method (100) additionally comprises a step prior to step (b) of correction of the plurality of satellite images (1).

[0405] In one embodiment, method (100) comprises a step prior to step (b) of fixed noise elimination in the plurality of satellite images (1), the fixed noise elimination step being carried out by means of the following sub-steps:

[0406] - determining a dark image by averaging a plurality of images acquired by the sensor (2) in the absence of Nomination; and

[0407] - subtracting the dark image from each of the satellite images.

[0408] In one embodiment, method (100) further comprises a step prior to step (b) of correction of the plurality of satellite images (1) to compensate for the non-uniform effect on the light of the sensor (2), the correction step being carried out by means of the following sub-steps:

[0409] - acquiring a plurality of images of a uniformly illuminated object and subtracting, from each one, the dark image;

[0410] - determining a corrective image by averaging the plurality of images from the previous sub-stage;

[0411] - normalizing the corrective image; and

[0412] - dividing each of the fixed noise-free satellite images by the normalized corrective image.

[0413] In one embodiment, the method (100) comprises a step prior to step (b) of variable noise elimination in each image of the plurality of satellite images (1), the variable noise elimination step being carried out by the following sub-steps: i) calculating the median of each row of the satellite image, and replacing the value of the pixels in the row that exceed a predefined threshold by the value of said median; ii) estimating the modulus and phase of the Fourier transform of the satellite image resulting from the previous sub-step; iii) calculating the median of the modulus of the Fourier transform; iv) determining the pixels of the modulus of the Fourier transform that are outside a predetermined radius measured from the center of said modulus of the Fourier transform; v) replacing the values ​​of the pixels determined in the previous sub-step by the median calculated in sub-step iii);vi) obtain a smoothed satellite image by applying the inverse Fourier transform; vii) obtain a high-frequency image by subtracting the smoothed satellite image from the satellite image; viii) for each row of the satellite image: calculate the median of the corresponding row of the high-frequency image, and subtract this median from each pixel of the row of the satellite image.

[0414] In one realization, method (100) comprises a stage prior to stage (b) of correction of the distortion geometries in the plurality of satellite images (1), the distortion geometries correction stage being carried out by means of the following sub-stages:

[0415] - divide the plurality of satellite images (1) into at least one set of satellite images, each of the sets of satellite images comprising a plurality of consecutive satellite images;

[0416] - Select a representative image for each set of satellite images; and for each representative image of each set of satellite images:

[0417] - estimate a representative time instant at which the sensor (2) acquired the representative image from the telemetry data (4);

[0418] - estimate the position and attitude of the satellite (3) at the estimated representative time instant from the telemetry data (4);

[0419] - calculate the pointing vectors of the sensor pixels (2) with respect to the inertial reference system centered on a celestial body (14) from the position and attitude of the satellite (3) at the representative time instant;

[0420] - obtain the projection of each pixel of the sensor (2) onto the surface of the celestial body (14) by calculating the intersection of the pointing vectors of the pixels of the sensor (2) with a reference geometric figure of a representative geodetic model of the celestial body (14);

[0421] - estimate the effective height of each pixel of the sensor (2) by calculating the magnitude of its corresponding pointing vector, where the terminal point of each pointing vector is its intersection with the surface of the celestial body (14);

[0422] - Calculate the ground sampling distance, GSD, of each sensor pixel (2), which represents the size of the projection of the sensor pixel (2) onto the surface of the celestial body (14): sensor pixel size GSD = effective height ■ - - - - — — - sensor focal length

[0423] - select the minimum GSD from among all the GSDs calculated for all pixels of the sensor (2);

[0424] - for each pixel of the sensor (2), divide the distance between the projection of the pixel and the projection of the pixel from the center of the sensor (2) by the minimum GSD, thus obtaining an equivalent number of pixels for that distance;

[0425] - generate a fictitious image without geometric distortion from the equivalent number of pixels calculated in the previous sub-stage for each pixel of the sensor (2);

[0426] - generate a coordinate map that defines the relationship between the coordinates of the sensor pixels (2) and the coordinates of the pixels of the fictitious image without geometric distortion, where the coordinates of the pixels of the satellite images coincide with the coordinates of the pixels of the sensor (2);

[0427] - obtain a plurality of corrected images by applying the coordinate map to the set of satellite images, where the coordinates of the pixels of the corrected images coincide with the coordinates of the pixels of the fictitious image without geometric distortion, and the value of each pixel of a corrected image is obtained by an interpolation of the pixel values ​​of a satellite image whose coordinates correspond to the coordinates of the pixel of the corrected image according to the coordinate map, the satellite image being part of the set of satellite images to which the representative image belongs.

[0428] In one embodiment, method (100) further comprises a step prior to step (b) of correcting optical and geometric distortions in the plurality of satellite images (1), the optical and geometric distortion correction step being carried out by means of the following sub-steps:

[0429] - receive a distortion profile from the sensor (2) comprising: either the position of the sensor pixel (2) acting as the optical distortion center, or an optical distortion function indicating the optical distortion value as a function of the distance to the distortion center;

[0430] - for each pixel of the sensor (2), calculate the distance of the pixel to the center of optical distortion, the distance measured according to a fourth pre-established standard, preferably the Euclidean standard;

[0431] - for each pixel of the sensor (2), calculate the optical distortion value using the optical distortion function and the distance calculated in the previous sub-stage; where the optical distortion value indicates the displacement that would have to be made in the pixel to obtain a fictitious image without optical distortion;

[0432] - generate a first coordinate map that defines the relationship between the coordinates of the sensor pixels (2) and the coordinates of the pixels of the fictitious image without optical distortion, where the coordinates of the pixels of the image without optical distortion are calculated by adding to the coordinates of the sensor pixels (2) their corresponding optical distortion values ​​obtained in the previous sub-stage and where the coordinates of the satellite image pixels (1) coincide with the coordinates of the sensor pixels (2);

[0433] - divide the plurality of satellite images (1) into at least one set of satellite images, each of the sets of satellite images comprising a plurality of consecutive satellite images;

[0434] - Select a representative image for each set of satellite images; and for each representative image of each set of satellite images:

[0435] - estimate a representative time instant at which the sensor (2) acquired the representative image from the telemetry data (4);

[0436] - estimate the position and attitude of the satellite (3) at the estimated representative time instant from the telemetry data (4);

[0437] - calculate the pointing vectors of the sensor pixels (2) with respect to the inertial reference system centered on the celestial body (14) from the position and attitude of the satellite (3) at the representative time instant;

[0438] - obtain the projection of each pixel of the sensor (2) onto the surface of the celestial body (14) by calculating the intersection of the pointing vectors of the pixels of the sensor (2) with a reference geometric figure of a representative geodetic model of the celestial body;

[0439] - estimate the effective height of each pixel of the sensor (2) by calculating the magnitude of its corresponding pointing vector, where the terminal point of each pointing vector is its intersection with the surface of the celestial body;

[0440] - Calculate the ground sampling distance, GSD, of each sensor pixel (2), which represents the size of the projection of the sensor pixel (2) onto the surface of the celestial body: sensor pixel size GSD = effective height ■ - - - - — — - sensor focal length

[0441] - select the minimum GSD from among all the GSDs calculated for all pixels of the sensor (2);

[0442] - for each pixel of the sensor (2), divide the distance between the projection of the pixel and the projection of the pixel of the center of the sensor (2) by the minimum GSD, thus obtaining an equivalent number of pixels for that distance;

[0443] - generate a fictitious image without optical and geometric distortions from the equivalent number of pixels calculated in the previous sub-stage for each pixel of the sensor (2);

[0444] - generate a second coordinate map that defines the relationship between the coordinates of the pixels of a fictitious image without optical distortion and with geometric distortion, and the coordinates of the pixels of the fictitious image without optical and geometric distortion generated; where the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion coincide with the coordinates of the sensor pixels (2) and with the coordinates of the pixels of the satellite images (1), and it is assumed that said image is free of optical distortion;

[0445] - generate a combined coordinate map, which defines the relation between the coordinates of the pixels of the fictitious image without optical and geometric distortions, and the coordinates of the sensor pixels (2), where:

[0446] (i) Evaluating the second coordinate map, the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion are identified that correspond to the coordinates of the pixels of the fictitious image without optical and geometric distortion;

[0447] (ii) Evaluating the first coordinate map, the coordinates of the pixels of the fictitious image without optical distortion are estimated by interpolation, which correspond to the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion identified in step (i),

[0448] (iii) Evaluating the first coordinate map, the coordinates of the sensor pixels (2) that correspond to the coordinates of the optically undistorted fictitious image pixels estimated in step (ii) are determined by interpolation; and

[0449] (iv) the combined coordinate map is generated by establishing the relation between the sensor pixel coordinates (2) determined in step (iii) and the pixel coordinates of the fictitious image without optical and geometric distortions; and

[0450] - obtain a plurality of corrected images by applying the combined coordinate map to the set of satellite images, where the coordinates of the pixels of the corrected images coincide with the coordinates of the pixels of the fictitious image without optical distortions and geometries, and the value of each pixel of a corrected image is obtained by combining the values ​​of the pixels of a satellite image whose coordinates correspond to the coordinates of the pixel of the corrected image according to the combined coordinate map, the satellite image being part of the set of satellite images to which the representative image belongs.

[0451] Figures 6a-6c show the effects of optical distortion and geometric distortion on satellite images (1).

[0452] Figure 6a (left image) shows a pincushion distortion, in which the magnification of the satellite image (1) increases with distance from the optical axis. The visible effect is that lines not passing through the center of image (1) curve inwards towards the center. Therefore, correcting this type of distortion, according to previous realizations, involves correcting the magnification effect so that the lines do not curve inwards, as shown in the dashed line on the right side of Figure 6a.

[0453] The geometric distortion present in satellite images (1) of a celestial body (14) is due to the location of the sensor (2) relative to the surface of the celestial body (14) at the time of acquisition. Figure 6b shows how, in images taken very off-nadir (that is, those in which parts of the surface to be photographed are not located directly below the sensor, but are deviated at a certain angle, as in point (b) in Figure 6b), the actual distance between two objects located in neighboring pixels varies enormously.

[0454] Geometric distortion means that the images (1) have a different GSD depending on whether the sensor pixel (2) is located nadir (point (a) in figure 6b) or very off-nadir (point (b) in figure 6b).

[0455] Figure 6c (left) shows a simplified scheme of the GSD variations depending on the location of each pixel of the sensor (2).

[0456] The correction of geometric distortion, therefore, consists of creating a new image (right image in Figure 6c) in which the distance represented by each pixel is constant. The number of pixels in the corrected image will be greater than the number of pixels in the satellite image (1), and greater the further off-nadir the image was acquired.

[0457] In a realization of method (100):

[0458] - the estimation of the first time instant in which the sensor (2) acquired the first satellite image (1.1) is carried out by interpolating the plurality of telemetry time instants (4); and / or

[0459] - the estimation of the second time instant in which the sensor (2) acquired the second satellite image (1.2) is carried out by interpolating the plurality of telemetry time instants (4); and / or

[0460] - the estimation of the representative time instant at which the sensor (2) acquired the representative satellite image is carried out by interpolating the plurality of telemetry time instants (4); and / or

[0461] - The estimation of the position and attitude of the satellite (3) at the first time instant, and / or at the second time instant and / or at the estimated representative time instant is carried out by interpolating the plurality of positions and attitudes of the satellite (3) from the telemetry data (4).

[0462] Figure 7 shows a system (20) for processing satellite images (1) according to an embodiment of the invention. The system (20) for processing satellite images (1) comprises a processor (22) configured to carry out steps (e.2) and (e.3) of the method and a programmable logic gate array (21), FPGA, configured to carry out the remaining steps of the method.

[0463] In one realization, the FPGA (21) is additionally configured to perform, in parallel:

[0464] - at least one of the sub-stages of stage (b) of the method; and / or

[0465] - stage (c) of the method; and / or

[0466] - stage (d) of the method; and / or

[0467] - at least one of the sub-stages of stage (e) of the method.

[0468] Figure 8 shows a schematic of the FPGA architecture (21) of the system (20) for carrying out the alignment of a pair of satellite images (1.1, 1.2) according to an embodiment of the invention.

[0469] In this figure, complex numbers are shown followed or preceded by an asterisk. All other numbers are real numbers.

[0470] In this figure, modules with double borders can have overflow.

[0471] The FPGA (21) of the embodiment shown in this figure 8 comprises a direct memory access module, DMA, configured to access the image pair (1.1 , 1.2) to be aligned and a write-direct memory access module, DMA, configured to store, in a storage element, the offset calculated by the FPGA (21) during the alignment of an image pair (1.1 , 1.2).

[0472] In one implementation, the storage element is a memory, a disk, or a database that is in communication with the FPGA (21). In one implementation, the storage element is a DDR memory (double data rate memory).

[0473] In the implementation shown in Figure 8, the FPGA (21) determines the displacement by minimizing a linear cross-correlation between the first satellite image (1.1) and the second satellite image (1.2) to be aligned. Thus, the FPGA (21) comprises a cross-correlation module configured to perform a linear cross-correlation between the images (1.1, 1.2). In particular, the cross-correlation module comprises a plurality of sub-modules configured to perform the operations required by the cross-correlation; for example, a module that applies a Vonann matrix, a module that calculates the two-dimensional Fourier transform, or a complex number multiplier module.

[0474] Additionally, the FPGA (21) comprises modules to adjust the precision of the integer and fractional parts of the result of each operation it performs.

[0475] Figure 9 (left) shows a schematic of the FPGA architecture (21) of the system (20) for carrying out the oversampled fused image acquisition (13) according to an embodiment of the invention.

[0476] For this fusion operation, the FPGA (21) comprises two direct memory access (DMA) modules configured to access the reference image (10) and the set of common regions (11) previously stored in the system storage element (20).

[0477] The FPGA (21) further comprises a demultiplexer, a plurality of image and weight storage modules, and a fusion module, wherein: the demultiplexer is configured to receive, by means of at least one DMA memory access module, the reference image (10) and the set of common regions (11) together with the weights estimated in sub-step (e.3) of the method (100) by the processor (22); each image and weight storage module is configured to store and send the reference image (10) or a common region from the set of common regions (11) together with its weights to the fusion module; and the fusion module is configured to receive the weights and receive and fuse the reference image (10) and the set of common regions (11) based on said weights.

[0478] In one embodiment, each image and weight storage module is configured to store a plurality of rows from the reference image (10) or from a common region (11), and to send part of the plurality of rows from the reference image (10) or from the common region (11) to the merging module while temporarily storing the rest of the rows.

[0479] In one implementation, the image and weight storage modules are configured to send to the fusion module at the same time.

[0480] Figure 9 shows that the image storage and weight modules store four rows and send three rows to the fusion module along with the corresponding weights for the pixels of the reference image (10) and the set of common regions (11) to be fused. The number of pixels from the reference image (10) and the set of common regions (11) that contribute to the generation of a single pixel of the oversampled fused image is configurable. In the example in Figure 9, four pixels of the oversampled fused image are calculated in parallel during the same clock cycle.

[0481] The FPGA (21) of Figure 9 further comprises an output memory comprising a plurality of RAM memory blocks, e.g. BRAMs, configured to store two lines of the oversampled fused image and to send such lines to the processor; e.g., to store the oversampled fused image in a storage element or to send it to an additional system module or to a receiver located on Earth.

[0482] The FPGA (21) in Figure 9 further comprises a special width converter configured to change the data width of the communication interface used, in this case AXI Stream, to go from 48 bits input (4 pixels * 12 bits) to 64 bits output, a data size compatible with those supported by a DMA module to which said special width converter is connected.

[0483] The FPGA (21) in Figure 9 further comprises a write-only direct memory access (DMA) module configured to send the oversampled fused image to a processor.

[0484] Figure 9 (center and right) shows in detail the demultiplexer, the image and weight storage modules, the output memory structure, and the output memory buffer.

[0485] The system (20) of Figure 7 further comprises a processor (22) comprising at least one of the following threads:

[0486] - a first thread (22.1) configured to carry out the fixed noise correction stages, and / or the variable noise correction stages and / or the flat correction stage of method (100); or

[0487] - a second thread (22.2) configured to carry out the optical distortion correction stage of method (100) or the optical and geometric distortion correction stage of method (100); or

[0488] - a third thread (22.3) configured to carry out the geometric distortion correction stage of method (100) or the optical and geometric distortion correction stage of method (100); or

[0489] - a fourth thread (22.4) to command the FPGA to execute steps (a)-(e) of method (100); or

[0490] - a combination of any of the above; and wherein the processor threads (3) are further configured to carry out the stages in parallel. Figure 10 shows a camera (18) for observing celestial bodies (14) comprising:

[0491] - a computer (19) comprising two systems (20) as shown in Figure 7 (although this number may be greater or less); and

[0492] - a sensor (2) configured to acquire a plurality of satellite images (1).

[0493] Figure 11 shows a satellite (3) comprising the celestial body observation camera (19) (14) shown in Figure 10.

Claims

1. CLAIMS 1.- Method (100) implemented by computer for the processing of satellite imageries (1), the method comprising the steps of: (a) receiving a plurality of satellite imageries (1) acquired by a sensor (2) on board a satellite (3) and telemetry data (4), the telemetry data (4) comprising a plurality of positions and attitudes of the satellite (3) with respect to an inertial reference system centered on a celestial body (14) at a plurality of telemetry time instants; (b) aligning a pair of consecutive satellite images from the plurality of satellite images (1) by carrying out the following sub-steps: (b.1) select a first satellite image (1.1) and a second satellite image (1.2) to align, (b.2) Calculate an initial horizontal displacement (8.1) and an initial vertical displacement (8.2) of the second satellite image (1.2) with respect to the first satellite image (1.1): - estimating a first time instant in which the sensor (2) acquired the first satellite image (1.1) from the telemetry data (4); - estimating the position and attitude of the satellite (3) at the first estimated time instant from the telemetry data (4); - estimating a second time instant in which the sensor (2) acquired the second satellite image (1.2) from the telemetry data (4); - estimating the position and attitude of the satellite (3) at the second time instant estimated from the telemetry data (4); - calculating the initial horizontal displacement (8.1) and the initial vertical displacement (8.2) from the positions and attitudes of the satellite (3) at the first time instant and at the second time instant; (b.3) Starting from the initial horizontal displacement (8.1) and the initial vertical displacement (8.2), calculate the total horizontal displacement (9.1) and the total vertical displacement (9.2) that must be applied to the second satellite image (1.2) to minimize the distance between said second satellite image (1.2) and the first satellite image (1.1), the distance measured according to a first pre-established standard, preferably the Euclidean standard; (b.4) align the first satellite image (1.1) and the second satellite image (1.2) by applying the total horizontal displacement (9.1) and the total vertical displacement (9.2) to the second satellite image (1.2); (c) repeat step (b) for all pairs of consecutive satellite imageries of the plurality of satellite imageries (1), wherein in each iteration of step (b): - the first satellite image (1.1) selected in sub-stage (b.1) corresponds to the second satellite image (1.2) of the previous iteration; - to the total horizontal displacement (9.1) obtained in sub-stage (b.3) must be added the total horizontal displacement (9.1) obtained in sub-stage (b.3) of the previous iteration; and - to the total vertical displacement (9.2) obtained in sub-stage (b.3) must be added the total vertical displacement (9.2) obtained in sub-stage (b.3) of the previous iteration; (d) select as the reference image (10) the first satellite image (1.1) selected in sub-stage (b.1) of the first iteration of stage (b), and obtain at least one common region between the reference image (10) and each of the second satellite images (1.2) shifted in sub-stage (b.4) of all iterations of stage (b), resulting in a set of common regions (11); (e) obtain an oversampled fused image (13) by merging the reference image (10) and the set of common regions (11) by carrying out the following sub-steps: (e.1) generate an oversampled image (12), wherein the resolution of the oversampled image (12) is the resolution of the reference image (10) multiplied by a predefined oversampling factor; (e.2) for each pixel (12.1) of the oversampled image (12), determine the distance from the center of said pixel (12.1) to the center of at least one pixel in each region of the set of regions (11) and to the center of at least one pixel in the reference image (10), wherein the distance is measured according to a second pre-established norm, preferably the Euclidean norm; (e.3) for each pixel (12.1) of the oversampled image (12.1), estimate a weight for at least one pixel from each region of the set of regions (11) and for at least one pixel from the reference image (10) as a function of the determined distances; and (e.4) obtain the oversampled fused image (13), wherein the value of each pixel (13.1) of the oversampled fused image (13) is the sum of the values ​​of at least one pixel from each region of the set of regions (11) and of at least one pixel from the reference image (10) weighted according to the weights estimated in the previous sub-stage.

2. Method (100) according to the preceding claim, wherein in sub-step (b.2) the calculation of the initial horizontal displacement (8.1) and the initial vertical displacement (8.2) from the positions and attitudes of the satellite (3) at the first time instant and at the second time instant is carried out by means of the following sub-steps: (1) in the first iteration of stage (b), for each pair of consecutive telemetry time instants, t teil yt tei2 , from the telemetry data (4): - calculate the pointing vectors of the sensor center pixel (2) and of at least two pixels of two sensor segments (2), which are perpendicular to each other and intersect at the center of the sensor (2), with respect to the inertial reference system from the position and attitude of the satellite (3) at the first telemetry time instant, and at the second telemetry time instant, t tei2 ; - obtain the projection of the pixel of the center of the sensor (2) and of the at least two pixels of the segments of the sensor (2) onto the surface of the celestial body (14) by calculating the intersection of the pointing vectors of said pixels with a reference geometric figure of a geodetic model representative of the celestial body (14); - define a projected reference system (15) as the projection of at least two pixels of the sensor segments (2) and the projection of the pixel at the center of the sensor (2) onto the surface of the celestial body (14) assuming that: either the surface of the celestial body (14) is flat; or the projected reference system (15) is an orthogonal reference system formed by the projections of the pixels of the sensor segments (2) and the projection of the pixel at the center of the sensor (2); and the projection of the pixel at the center of the sensor (2) is the center of the projected reference system (15); - for the sensor center pixel (2): calculate the effective height as the magnitude of its pointing vector, where the terminal point of the pointing vector is its intersection with the surface of the celestial body (14), and calculate the ground sampling distance, GSD, which represents the size of the projection of the sensor center pixel (2) onto the surface of the celestial body (14): sensor pixel size GSD = effective height ■ - - - - — — - sensor focal length - calculate a first linear velocity, vf of the projection of the pixel of the center of the sensor (2) due to the rotation of the celestial body (14), where the first linear velocity v^ is referenced in the inertial reference system, and the first linear velocity It is calculated from the angular velocity, co, of the celestial body (14) and the perpendicular distance from the projection of the pixel of the sensor center to the rotation axis of the celestial body (14), R: where d) and R are vectors referenced in the inertial reference system; - calculate a second linear speed, vj, of the projection of the pixel of the center of the sensor (2) due to the movement of the satellite (3): where: o pos^ is the position vector of the projection of the pixel of the sensor center (2) at the second telemetry time instant, t tei2 , the vector being referenced in the inertial reference system; polf is the position vector of the projection of the pixel of the sensor center (2) at the first telemetry time instant, the vector being referenced in the inertial reference frame; I vj is referenced in the inertial reference frame; - Calculate the total linear speed, vf, of the projection of the pixel at the center of the sensor (2) in pixels per second: - reference the total linear velocity vf obtained in the previous sub-stage in the projected reference system (15); (2) In the first iteration of stage (b), calculate the horizontal linear velocities, v H (t), and vertical, v v (t), as a function of time by interpolating the horizontal and vertical components of the total linear velocities, vf, obtained for all pairs of consecutive telemetry time instants; and (3) In any iteration of step (b), for the pair of consecutive satellite images to be aligned, estimate the initial horizontal displacement (8.1) and the initial vertical displacement (8.2) by integrating the horizontal and vertical linear velocities with respect to time, v H (t) yv K(t), respectively, between the first time instant and the second time instant.

3. Method (100) according to any of the preceding claims, wherein sub-step (b.3) is carried out by means of the following sub-steps: (b.3.i) for the first satellite image (1.1), generate a first pyramid (16) of images by subsampling the first satellite image (1.1) with N predefined subsampling factors, where N is greater than or equal to one, wherein the first pyramid of images (16) comprises N+1 levels, and wherein: • the first level of the first image pyramid (16) comprises the first satellite image (1.1), and • the i-th level of the first image pyramid (16) comprises an image subsampled from the first satellite image (1.1) with the k-th subsampling factor, where i = 2,... ,N+1 and k = i-1 ; (b.3.ii) for the second satellite image (1.2), generate a second image pyramid (17) by subsampling the second satellite image (1.2) with M predefined subsampling factors, M being greater than or equal to one, wherein the second image pyramid (17) comprises M+1 levels, and wherein: • M is less than or equal to N, • the first level of the second image pyramid (17) comprises the second satellite image (1.2), and • the j-th level of the second image pyramid (17) comprises a subsampled image (1.2') from the second satellite image (1.2) with the m-th subsampling factor, where j = 2,... ,M+1 and m = j-1 ; (b.3.iii) Starting from the initial horizontal displacement (8.1) and the initial vertical displacement (8.2), determine the horizontal and vertical displacements that must be applied to the image of level M+1 of the second image pyramid (17) to minimize the distance between that image and the image of the first image pyramid (16) of the same level, the distance measured according to a third pre-established norm, preferably the Euclidean norm; (b.3.iv) Starting from the horizontal and vertical displacements determined in the previous sub-step, determine the horizontal and vertical displacements that must be applied to the image of the lower level of the second imagery pyramid (17) to minimize the distance between said image and the image of the first imagery pyramid (16) of the same level, the distance measured according to the third pre-established norm, and where: • if that lower level is level 1 of the second image pyramid (17), set the horizontal displacement as total horizontal displacement (9.1) and the vertical displacement as total vertical displacement (9.2); or • if that lower level is not level 1 of the second pyramid of images (17), repeat step (b.3.iv).

4. Method (100) according to any of the preceding claims, further comprising a step prior to step (b) of fixed noise elimination in the plurality of satellite images (1), the fixed noise elimination step being carried out by means of the following sub-steps: - determining a dark image by averaging a plurality of images acquired by the sensor (2) in the absence of Nomination; and - subtracting the dark image from each of the satellite images.

5. Method (100) according to the preceding claim, further comprising a step prior to step (b) of correcting the plurality of satellite images (1) to compensate for the non-uniform effect on the light of the sensor (2), the correction step being carried out by means of the following sub-steps: - acquiring a plurality of images of a uniformly illuminated object and subtracting, from each one, the dark image; - determining a corrective image by averaging the plurality of images from the previous sub-stage; - normalizing the corrective image; and - dividing each of the fixed noise-free satellite images by the normalized corrective image.

6. Method (100) according to any of the preceding claims, further comprising a step prior to step (b) of variable noise elimination in each image of the plurality of satellite images (1), the variable noise elimination step being carried out by the following sub-steps: i) calculating the median of each row of the satellite image, and replacing the value of the pixels in the row that exceed a predefined threshold by the value of said medians; ii) estimate the modulus and phase of the Fourier transform of the satellite image resulting from the previous sub-stage; iii) calculate the median of the modulus of the Fourier transform; iv) determine the pixels of the modulus of the Fourier transform that are outside a predetermined radius measured from the center of said modulus of the Fourier transform; v) replace the pixel values ​​determined in the previous sub-stage with the median calculated in sub-stage iii); vi) obtain a smoothed satellite image by applying the inverse Fourier transform; vii) obtain a high-frequency image by subtracting the smoothed satellite image from the satellite image; viii) for each row of the satellite image: calculate the median of the corresponding row of the high-frequency image, and subtract said median from each pixel of the row of the satellite image.

7. Method (100) according to any of the preceding claims, further comprising a step prior to step (b) of correcting the optical distortion in the plurality of satellite images (1), the optical distortion correction step being carried out by means of the following sub-steps: - receive a distortion profile from the sensor (2) comprising: either the position of the sensor pixel (2) acting as the optical distortion center, or an optical distortion function indicating the optical distortion value as a function of the distance to the distortion center; - for each pixel of the sensor (2), calculate the distance of the pixel to the center of optical distortion, the distance measured according to a fourth pre-established standard, preferably the Euclidean standard; - for each pixel of the sensor (2), calculate the optical distortion value using the optical distortion function and the distance calculated in the previous sub-stage; where the optical distortion value indicates the displacement that would have to be made in the pixel to obtain a fictitious image without optical distortion; - generate a coordinate map that defines the relation between the coordinates of the sensor pixels (2) and the coordinates of the pixels of the fictitious image without optical distortion, wherein the coordinates of the pixels of the fictitious image without optical distortion are calculated by adding to the coordinates of the sensor pixels (2) their corresponding optical distortion values ​​obtained in the previous sub-stage and wherein the coordinates of the pixels of the satellite images (1) coincide with the coordinates of the sensor pixels (2); - obtain a plurality of corrected images by applying the coordinate map to the plurality of satellite images (1), where the coordinates of the pixels of the corrected images coincide with the coordinates of the pixels of the fictitious image without optical distortion, and the value of each pixel of a corrected image is obtained by interpolating the pixel values ​​of a satellite image whose coordinates correspond to the pixel coordinates of the corrected image according to the coordinate map.

8. Method (100) according to the preceding claim, further comprising a step prior to step (b) of correcting the geometric distortion in the plurality of satellite images (1), the geometric distortion correction step being carried out by means of the following sub-steps: - divide the plurality of satellite images (1) into at least one set of satellite images, each of the sets of satellite images comprising a plurality of consecutive satellite images; - Select a representative image for each set of satellite images; and for each representative image of each set of satellite images: - estimate a representative time instant at which the sensor (2) acquired the representative image from the telemetry data (4); - estimate the position and attitude of the satellite (3) at the estimated representative time instant from the telemetry data (4); - calculate the pointing vectors of the sensor pixels (2) with respect to the inertial reference system centered on a celestial body (14) from the position and attitude of the satellite (3) at the representative time instant; - obtain the projection of each pixel of the sensor (2) onto the surface of the celestial body (14) by calculating the intersection of the pointing vectors of the pixels of the sensor (2) with a reference geometric figure of a representative geodetic model of the celestial body (14); - estimate the effective height of each pixel of the sensor (2) by calculating the magnitude of its corresponding pointing vector, where the terminal point of each pointing vector is its intersection with the surface of the celestial body (14); - calculate the ground sampling distance, GSD, of each sensor pixel (2), which represents the size of the projection of the sensor pixel (2) onto the surface of the celestial body (14): sensor pixel size GSD = effective height ■ - - - - — — - sensor focal length - select the minimum GSD from among all the GSDs calculated for all pixels of the sensor (2); - for each pixel of the sensor (2), divide the distance between the projection of the pixel and the projection of the pixel from the center of the sensor (2) by the minimum GSD, thus obtaining an equivalent number of pixels for that distance; - generate a fictitious image without geometric distortion from the equivalent number of pixels calculated in the previous sub-stage for each pixel of the sensor (2); - generate a coordinate map that defines the relation between the coordinates of the sensor pixels (2) and the coordinates of the pixels of the fictitious image without geometric distortion, where the coordinates of the pixels of the satellite images coincide with the coordinates of the pixels of the sensor (2); - obtain a plurality of corrected images by applying the coordinate map to the set of satellite images, where the coordinates of the pixels of the corrected images coincide with the coordinates of the pixels of the fictitious image without geometric distortion, and the value of each pixel of a corrected image is obtained by an interpolation of the pixel values ​​of a satellite image whose coordinates correspond to the coordinates of the pixel of the corrected image according to the coordinate map, the satellite image being part of the set of satellite images to which the representative image belongs.

9. Method (100) according to any of claims 1 to 6, which additionally It comprises a stage prior to stage (b) of correcting optical and geometric distortions in the plurality of satellite images (1), the stage of correcting optical and geometric distortions being carried out by means of the following sub-stages: - receive a distortion profile from the sensor (2) comprising: either the position of the sensor pixel (2) acting as the optical distortion center, or an optical distortion function indicating the optical distortion value as a function of the distance to the distortion center; - for each pixel of the sensor (2), calculate the distance of the pixel to the center of optical distortion, the distance measured according to a fourth pre-established standard, preferably the Euclidean standard; - for each pixel of the sensor (2), calculate the optical distortion value using the optical distortion function and the distance calculated in the previous sub-stage; where the optical distortion value indicates the displacement that would have to be made in the pixel to obtain a fictitious image without optical distortion; - generate a first coordinate map that defines the relation between the coordinates of the sensor pixels (2) and the coordinates of the pixels of the fictitious image without optical distortion, where the coordinates of the pixels of the image without optical distortion are calculated by adding to the coordinates of the sensor pixels (2) their corresponding optical distortion values ​​obtained in the previous sub-stage and where the coordinates of the satellite image pixels (1) coincide with the coordinates of the sensor pixels (2); - divide the plurality of satellite images (1) into at least one set of satellite images, each of the sets of satellite images comprising a plurality of consecutive satellite images; - Select a representative image for each set of satellite images; and for each representative image of each set of satellite images: - estimate a representative time instant at which the sensor (2) acquired the representative image from the telemetry data (4); - estimate the position and attitude of the satellite (3) at the estimated representative time instant from the telemetry data (4); - calculate the pointing vectors of the sensor pixels (2) with respect to the inertial reference system centered on the celestial body (14) from the position and attitude of the satellite (3) at the representative time instant; - obtain the projection of each pixel of the sensor (2) onto the surface of the celestial body (14) by calculating the intersection of the pointing vectors of the pixels of the sensor (2) with a reference geometric figure of a representative geodetic model of the celestial body; - estimate the effective height of each pixel of the sensor (2) by calculating the magnitude of its corresponding pointing vector, where the terminal point of each pointing vector is its intersection with the surface of the celestial body (14); - Calculate the ground sampling distance, GSD, of each pixel of the sensor (2), which represents the size of the projection of the pixel of the sensor (2) onto the surface of the celestial body: pixel size of the sensor GSD = effective height ■ - - - - — — - sensor focal length - select the minimum GSD from among all the GSDs calculated for all pixels of the sensor (2); - for each pixel of the sensor (2), divide the distance between the projection of the pixel and the projection of the pixel of the center of the sensor (2) by the minimum GSD, thus obtaining an equivalent number of pixels for that distance; - generate a fictitious image without optical and geometric distortions from the equivalent number of pixels calculated in the previous sub-stage for each pixel of the sensor (2); - generate a second coordinate map that defines the relation between the coordinates of the pixels of a fictitious image without optical distortion and with geometric distortion, and the coordinates of the pixels of the fictitious image without optical and geometric distortion generated; where the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion coincide with the coordinates of the sensor pixels (2) and with the coordinates of the pixels of the satellite images (1), and it is assumed that said image is free of optical distortion; - generate a combined coordinate map, which defines the relation between the coordinates of the pixels of the fictitious image without optical and geometric distortions, and the coordinates of the sensor pixels (2), where: (i) Evaluating the second coordinate map, the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion are identified that correspond to the coordinates of the pixels of the fictitious image without optical and geometric distortions; (ii) Evaluating the first coordinate map, the coordinates of the pixels of the fictitious image without optical distortion are estimated by interpolation, which correspond to the coordinates of the pixels of the fictitious image without optical distortion and with geometric distortion identified in step (i), (iii) Evaluating the first coordinate map, the coordinates of the sensor pixels (2) that correspond to the coordinates of the optically undistorted fictitious image pixels estimated in step (ii) are determined by interpolation; and (iv) the combined coordinate map is generated by establishing the relation between the sensor pixel coordinates (2) determined in step (iii) and the pixel coordinates of the fictitious image without optical and geometric distortions; and - obtain a plurality of corrected images by applying the combined coordinate map to the set of satellite images, where the coordinates of the pixels of the corrected images coincide with the coordinates of the pixels of the fictitious image without optical and geometric distortions, and the value of each pixel of a corrected image is obtained by combining the values ​​of the pixels of a satellite image whose coordinates correspond to the coordinates of the pixel of the corrected image according to the combined coordinate map, the satellite image being part of the set of satellite images to which the representative image belongs.

10. Method (100) according to any of the preceding claims, wherein - the estimation of the first time instant in which the sensor (2) acquired the first satellite image (1.1) is carried out by interpolating the plurality of telemetry time instants (4); and / or - the estimation of the second time instant in which the sensor (2) acquired the second satellite image (1.2) is carried out by interpolating the plurality of telemetry time instants (4); and / or - the estimation of the representative time instant at which the sensor (2) acquired the representative satellite image is carried out by interpolating the plurality of telemetry time instants (4); and / or - the estimation of the position and attitude of the satellite (3) at the first time instant, and / or at the second time instant and / or at the time instant Representative estimates are carried out by interpolating the plurality of satellite positions and attitudes (3) from the telemetry data (4). 11.- System (20) for processing satellite images comprising a processor (22) configured to carry out steps (e.2) and (e.3) of the method according to claims 1 to 3 or 10; and a programmable logic gate array (21), FPGA, configured to carry out the remaining steps of the method according to claims 1 to 3 or 10.

12. System (20) according to the preceding claim, wherein the FPGA (21) is further configured to perform, in parallel: - at least one of the sub-stages of stage (b) of the method; and / or - stage (c) of the method; and / or - stage (d) of the method; and / or - at least one of the sub-stages of stage (e) of the method.

13. System (20) according to any of claims 11 or 12, wherein the processor (22) comprises at least one of the following threads: - a first thread (22.1) configured to carry out the correction steps of the method according to any of claims 4 to 6; or - a second thread (22.2) configured to carry out the correction stage of the method according to claim 7 or 9; or - a third thread (22.3) configured to carry out the correction stage of the method according to claim 8 or 9; or - a fourth thread (22.4) to command the FPGA to execute the steps of the method according to any of claims 1 to 3; or - a combination of any of the above; and wherein the processor threads (3) are additionally configured to carry out the stages in parallel.

14. Camera (18) for observing celestial bodies comprising: - a computer (19) comprising at least one system (20) according to claims 11-13; and - a sensor (2) configured to acquire a plurality of satellite images (1). 15.- Satellite (3) comprising the celestial body observation camera (19) according to claim 14.

Images