A remote sensing image fusion method considering image distortion and a computer readable medium

Abstract

The application provides a remote sensing image fusion method considering image distortion and a computer readable medium. The application solves an intersection polygon of an image participating in fusion, constructs a minimum circumscribed rectangle of the intersection polygon according to a geographical range of the intersection polygon. After grid division of the minimum circumscribed rectangle, initial matching is performed, and affine transformation processing is performed on the corrected hyperspectral resolution image based on the matching result. The intersection polygon of the image participating in fusion and the corresponding minimum circumscribed rectangle are updated, grid division is performed to obtain a plurality of blocks, and fine matching processing is performed, and further polynomial correction is performed on the hyperspectral resolution remote sensing image subjected to the affine transformation processing based on the matching result. Each block of the image participating in fusion is subjected to fusion processing, so that a fused remote sensing image is obtained. The application can effectively solve the problem that the fusion effect is poor due to the internal distortion of the image, and can realize fusion between large-width remote sensing images.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, specifically relating to a remote sensing image fusion method that takes into account image distortion and a computer-readable medium. Background Technology

[0002] The fusion of high spatial resolution and high spectral resolution images is called remote sensing image spatial-spectral fusion. Because it can overcome the inherent contradiction between the spatial and spectral resolution of an imaging system, spatial-spectral fusion is a hot topic in the field of remote sensing image fusion. However, in spatial-spectral fusion, differences in the data sources of spatial and spectral images and registration errors between images can significantly affect fusion performance. Furthermore, the wide coverage and significant internal topographic relief of remote sensing images often result in varying local geometric deformations, i.e., internal distortion. In such cases, using the same registration model for the entire image makes it difficult to achieve high registration accuracy, further limiting the fusion effect. How to effectively improve registration accuracy is a pressing problem that spatial-spectral fusion needs to solve. Summary of the Invention

[0003] This invention addresses the shortcomings of existing technologies by providing a remote sensing image fusion method and computer-readable medium that takes into account image distortion. Considering that image distortion leads to low registration accuracy and thus affects the fusion effect, this invention combines coarse registration and block-based fine registration methods during fusion to propose a remote sensing image fusion method that takes into account image distortion. This method is suitable for fusing high spatial resolution and high spectral resolution images and can provide technical support for the further application of remote sensing technology.

[0004] To achieve the above objectives, the technical solution provided by this invention is a remote sensing image fusion method that takes into account image distortion, comprising the following steps:

[0005] Step 1: Acquire the original high spatial resolution remote sensing image and its corresponding RPC positioning model, and perform geometric correction to obtain the corrected high spatial resolution remote sensing image; acquire the high spectral resolution remote sensing image and its corresponding RPC positioning model, and perform geometric correction to obtain the corrected high spectral resolution remote sensing image.

[0006] Step 2: Calculate the intersection polygons of the corrected high spatial resolution image and the corrected hyperspectral resolution image, and obtain the geographic extent of the intersection polygons. Construct the minimum bounding rectangle of the intersection polygons based on the geographic extent of the intersection polygons. Perform grid division on the minimum bounding rectangle of the intersection polygons to obtain multiple blocks of minimum bounding rectangles;

[0007] Step 3: Map each block of the minimum bounding rectangle to the corrected high spatial resolution remote sensing image to obtain each corrected high spatial resolution remote sensing block image; map each block of the minimum bounding rectangle to the corrected high spectral resolution remote sensing image to obtain each corrected high spectral resolution remote sensing block image; and perform initial matching processing between each corrected high spatial resolution remote sensing block image and its corresponding corrected high spectral resolution remote sensing block image to obtain multiple pairs of corresponding points for each block.

[0008] Step 4: Construct an affine transformation model based on multiple sets of corresponding point pairs in multiple blocks, and perform affine transformation processing on the corrected hyperspectral resolution image according to the affine transformation model to obtain an affine transformed hyperspectral resolution image.

[0009] Step 5: Calculate the updated intersecting polygons of the corrected high spatial resolution image and the affine-transformed high spectral resolution image, and obtain the geographical extent of the updated intersecting polygons. Construct the minimum bounding rectangle of the updated intersecting polygons based on the geographical extent of the updated intersecting polygons. Perform grid division on the minimum bounding rectangle of the updated intersecting polygons to obtain multiple blocks of minimum bounding rectangles;

[0010] Step 6: Map each block of the minimum bounding rectangle to the corrected high spatial resolution remote sensing image to obtain each corrected high spatial resolution remote sensing block image; map each block of the minimum bounding rectangle to the affine transformed high spectral resolution remote sensing image to obtain each corrected high spectral resolution remote sensing block image; and perform fine matching processing between each corrected high spatial resolution remote sensing block image and the corresponding affine transformed high spectral resolution remote sensing block image to obtain multiple pairs of corresponding points for each block.

[0011] Step 7: Construct the polynomial transformation model of each block based on multiple pairs of corresponding points in each block. Combine the multiple pairs of corresponding points in each sub-block and use the least squares method to solve for the parameters of the polynomial transformation model of each sub-block, thus obtaining the polynomial transformation model of each block.

[0012] Step 8: Perform quadratic polynomial correction on the affine-transformed hyperspectral resolution remote sensing image based on the quadratic polynomial transformation model of each block to obtain the hyperspectral resolution remote sensing block image of each quadratic polynomial correction.

[0013] Step 9: Perform fusion processing on each corrected high spatial resolution remote sensing block image and the corresponding quadratic polynomial corrected high spectral resolution remote sensing block image to obtain a fused high spatial resolution and high spectral resolution remote sensing image.

[0014] Preferably, step 2 calculates the intersection polygon of the corrected high spatial resolution image and the corrected hyperspectral resolution image, and the specific process is as follows:

[0015] Construct corrected high spatial resolution image vector polygons based on the geographic coordinates of the four corner points of the corrected high spatial resolution image;

[0016] Construct corrected hyperspectral resolution image vector polygons based on the geographic coordinates of the four corner points of the corrected hyperspectral resolution image;

[0017] The intersection polygon between the corrected high spatial resolution image vector polygon and the corrected high spectral resolution image vector polygon is calculated using a spatial intersection calculation method.

[0018] Step 2 involves obtaining the geographical extent of the intersecting polygons and constructing the minimum bounding rectangle of the intersecting region based on the geographical extent of the intersecting polygons. The specific process is as follows:

[0019] Obtain the longitude and latitude of each point in the intersecting polygon in sequence;

[0020] To construct the longitude range of an intersecting polygon, we select the minimum and maximum longitude values ​​for each point within the polygon from among multiple longitudes.

[0021] [FMIN_LON, FMAX_LON]

[0022] Where FMIN_LON represents the minimum longitude of each point in the intersecting polygon, and FMAX_LON represents the maximum longitude of each point in the intersecting polygon;

[0023] To construct the dimensional range of an intersecting polygon, we select the minimum and maximum dimensional values ​​of each point within the polygon from among multiple points in the polygon, as follows:

[0024] [FMIN_LAT, FMAX_LAT]

[0025] Where FMIN_LAT represents the minimum latitude of each point in the intersecting polygon, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygon;

[0026] Use (FMAX_LON, FMIN_LAT), (FMAX_LON, FMAX_LAT), (FMAX_LON, FMAX_LAT), and (FMAX_LON, FMIN_LAT) as the four vertices of the minimum bounding rectangle of the intersecting polygon to construct the minimum bounding rectangle of the intersecting polygon.

[0027] Step 2 involves performing grid division on the minimum bounding rectangle of the intersecting polygons to obtain multiple blocks of minimum bounding rectangles, as detailed below:

[0028] Set the number of grid divisions in the longitude direction to nFBK_lon, and set the number of grid divisions in the latitude direction to nFBK_lat;

[0029] Calculate the longitude width and latitude height of each minimum bounding rectangle block as follows:

[0030]

[0031] Where ΔF_lon is the longitude width of each minimum bounding rectangle block, ΔF_lat is the latitude height of each minimum bounding rectangle block, nFBK_lon is the number of grid divisions in the longitude direction, nFBK_lat is the set number of grid divisions in the latitude direction, FMIN_LON represents the minimum longitude of each point in the intersecting polygon, and FMAX_LON represents the maximum longitude of each point in the intersecting polygon. Where FMIN_LAT represents the minimum latitude of each point in the intersecting polygon, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygon.

[0032] The smallest bounding rectangle block in row j and column i is defined as:

[0033] FBK i，j

[0034] i∈[1,nFBK_lon],j∈[1,nFBK_lat]

[0035] Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_lat is the number of latitudinal grid divisions.

[0036] Calculate the minimum longitude, maximum longitude, minimum latitude, and maximum latitude of the minimum bounding rectangle of the minimum bounding rectangle block in row j and column i, as follows:

[0037]

[0038] in, Let represent the minimum longitude of the smallest bounding rectangle block in row j and column i. This represents the maximum longitude of the smallest bounding rectangular block in row j and column i. This represents the minimum latitude of the smallest bounding rectangle block in the j-th row and i-th column. This represents the maximum latitude of the smallest bounding rectangle block in the j-th row and i-th column;

[0039] FMIN_LON represents the minimum longitude of each point in the intersecting polygons, and FMAX_LON represents the maximum longitude of each point in the intersecting polygons; where FMIN_LAT represents the minimum latitude of each point in the intersecting polygons, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygons.

[0040] Will As four vertices, construct the geographic rectangle corresponding to the minimum bounding rectangle of the j-th row and i-th column;

[0041] Preferably, step 3, which maps the image to the corrected high spatial resolution remote sensing image, yields each corrected high spatial resolution remote sensing block image. The specific process is as follows:

[0042] Divide the four vertices into blocks based on the minimum bounding rectangle in row j and column i. The geographic range is composed of image content obtained from corrected high spatial resolution remote sensing imagery.

[0043] The corrected high spatial resolution remote sensing patch image in row j and column i is defined as:

[0044] HSPA_FBK i，j

[0045] i∈[1,nFBK_lon],j∈[1,nFBK_lat]

[0046] Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_lat is the number of latitude grid divisions.

[0047] Step 3, which maps the image to the corrected hyperspectral resolution remote sensing image, yields each hyperspectral resolution remote sensing block image. The specific process is as follows:

[0048] Divide the four vertices into blocks based on the minimum bounding rectangle in row j and column i. The geographic range is composed of image content obtained from corrected hyperspectral resolution remote sensing images.

[0049] The corrected hyperspectral resolution remote sensing patch image in row j and column i is defined as:

[0050] HSPE_FBK i，j

[0051] i∈[1,nFBK_lon],j∈[1,nFBK_lat]

[0052] Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_lat is the number of latitudinal grid divisions.

[0053] Step 3 involves initial matching to obtain multiple pairs of identically named points for each block. The specific process is as follows:

[0054] Corrected high spatial resolution remote sensing patch image HSPA_FBK in row j and column i i，j and the corresponding corrected hyperspectral resolution remote sensing patch image HSPE_FBK in row j and column i. i，j The SIFT+KNN matching algorithm was used to obtain the corrected high spatial resolution remote sensing patch image HSPA_FBK in the j-th row and i-th column. i，j Pixel coordinates and corrected hyperspectral resolution remote sensing patch image HSPE_FBK in row j and column i i，j Pixel coordinates.

[0055] Get the minimum longitude of the sub-block in row j and column i. minimum latitude The high spatial resolution remote sensing patch image HSPA_FBK in row j and column i is corrected. i，j Pixel coordinates are mapped to geographic coordinates of high spatial resolution remote sensing block image by summing the product of pixel coordinates and resolution with the minimum value.

[0056] Get the minimum longitude of the sub-block in row j and column i. minimum latitude The corrected hyperspectral resolution remote sensing patch image HSPE_FBK in row j and column i. i，j Pixel coordinates are used to map hyperspectral resolution remote sensing block images to geographic coordinates by summing the product of pixel coordinates and resolution with the minimum value.

[0057] Merge the geographic coordinates of the high spatial resolution remote sensing block image and the high spectral resolution remote sensing block image corresponding to the sub-block in row j and column i into multiple pairs of points with the same name in row j and column i.

[0058] Preferably, step 4 involves constructing an affine transformation model based on multiple sets of corresponding point pairs across multiple blocks. The specific process is as follows:

[0059] A global affine transformation model is constructed based on the geographic coordinates of the corresponding high spatial resolution remote sensing block image and the high spectral resolution remote sensing block image in all block-named point pairs, and the parameters of the affine transformation model are solved by the least squares method.

[0060] Step 4 involves performing an affine transformation on the corrected hyperspectral resolution image according to an affine transformation model. The specific process is as follows:

[0061] By combining the affine transformation model parameters (e0, e1, e2, f0, f1, f2), an indirect correction method is used to perform affine transformation correction on the corrected hyperspectral resolution image, resulting in an affine transformed hyperspectral resolution image.

[0062] Preferably, the intersection polygon updated by calculating the corrected high spatial resolution image and the affine-transformed high spectral resolution image in step 5 is as follows:

[0063] A corrected high spatial resolution image vector polygon is constructed based on the geographic coordinates of the four corner points of the corrected high spatial resolution image. An affine transformed high spectral resolution image vector polygon is then constructed based on the geographic coordinates of the four corner points of the affine transformed high spectral resolution image. The intersection polygon between the corrected high spatial resolution image vector polygon and the affine transformed high spectral resolution image vector polygon is calculated using a spatial intersection calculation method; this is the updated intersection polygon.

[0064] Step 5 involves obtaining the geographical extent of the updated intersecting polygons and constructing the minimum bounding rectangle of the intersecting region based on the updated geographical extent of the intersecting polygons. The specific process is as follows:

[0065] Sequentially obtain the longitude and latitude of each point in the updated intersecting polygon;

[0066] To construct the longitude range of the updated intersecting polygon, we filter out the minimum and maximum longitude values ​​of each point within the updated polygon from the longitudes of multiple points.

[0067] [SMIN_LON, SMAX_LON]

[0068] Where SMIN_LON represents the minimum longitude of each point in the updated intersecting polygon, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygon;

[0069] In the updated intersecting polygon, the minimum and maximum latitude values ​​of each point are selected from the latitudes of multiple points in the updated polygon. The latitude range of the updated intersecting polygon is then constructed, defined as follows:

[0070] [SMIN_LAT, SMAX_LAT]

[0071] Where SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygon, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygon;

[0072] Use (SMAX_LON, SMIN_LAT), (SMAX_LON, SMAX_LAT), (SMAX_LON, SMAX_LAT), (SMAX_LON, SMIN_LAT) as the four vertices of the minimum bounding rectangle of the updated intersecting polygon to construct the minimum bounding rectangle of the updated intersecting polygon.

[0073] Step 5 involves performing grid division on the minimum bounding rectangle of the updated intersecting polygons to obtain multiple blocks of minimum bounding rectangles, as detailed below:

[0074] Set the number of grid divisions in the longitude direction to nSBK_lon, and set the number of grid divisions in the latitude direction to nSBK_lat;

[0075] Calculate the longitude width and latitude height of each minimum bounding rectangle block as follows:

[0076]

[0077] Where ΔS_lon is the longitude width of each minimum bounding rectangle block, ΔS_lat is the latitude height of each minimum bounding rectangle block, nSBK_lon is the number of grid divisions in the longitude direction, nSBK_lat is the set number of grid divisions in the latitude direction, SMIN_LON represents the minimum longitude of each point in the updated intersecting polygon, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygon; where SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygon, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygon.

[0078] The smallest bounding rectangle block in row q and column p is defined as:

[0079] SBK p，q

[0080] p∈[1,nSBK_lon], q∈[1,nSBK_lat]

[0081] Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_lon is the number of directional grid divisions, and nSBK_lat is the number of latitudinal grid divisions.

[0082] Calculate the minimum longitude, maximum longitude, minimum latitude, and maximum latitude of the minimum bounding rectangle of the minimum bounding rectangle of the q-th row and p-th column, as follows:

[0083]

[0084] in, This represents the minimum longitude of the smallest bounding rectangle block in row q and column p. This represents the maximum longitude of the smallest bounding rectangle in row q and column p. This represents the minimum latitude of the smallest bounding rectangle block in row q and column p. This represents the maximum latitude of the smallest bounding rectangle in row q and column p.

[0085] SMIN_LON represents the minimum longitude of each point in the updated intersecting polygon, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygon; where SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygon, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygon; As four vertices, construct the geographic rectangular block corresponding to the minimum bounding rectangular block in row q and column p;

[0086] Preferably, step 6 involves mapping the image to the corrected high spatial resolution remote sensing image to obtain each corrected high spatial resolution remote sensing block image. The specific process is as follows:

[0087] Divide the four vertices into blocks based on the minimum bounding rectangle in row q and column p. The geographic range is composed of image content obtained from corrected high spatial resolution remote sensing imagery.

[0088] The corrected high spatial resolution remote sensing patch image in row q and column p is defined as:

[0089] HSPA_SBK p，q

[0090] p∈[1,nSBK_lon], q∈[1,nSBK_lat]

[0091] Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_lon is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction.

[0092] Step 6 describes mapping the hyperspectral resolution remote sensing image to an affine transformation to obtain each corrected high spatial resolution remote sensing block image. The specific process is as follows:

[0093] Divide the four vertices into blocks based on the minimum bounding rectangle in row q and column p. The geographic range is composed of image content obtained from corrected high spatial resolution remote sensing imagery.

[0094] The high spatial resolution remote sensing patch image of the affine transformation in row q and column p is defined as:

[0095] HSPE_SBK p，q

[0096] p∈[1,nSBK_lon], q∈[1,nSBK_lat]

[0097] Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_ion is the number of directional grid divisions, and nSBK_lat is the number of latitude grid divisions.

[0098] Step 6 describes obtaining multiple pairs of identical points for each block through fine-grained matching. The specific process is as follows:

[0099] Corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p p，q The validity determination is as follows:

[0100] If CNT_HSPA_SBK p，q ≥α, then

[0101] If CNT_HSPA_SBK p，q <α, then

[0102] Among them, CNT_HSPA_SBK p，q This represents the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The number of pixels with values ​​greater than 0, where α represents the threshold for determining the validity of the block. This represents the validity attribute of the corrected high spatial resolution remote sensing block image in row q and column p.

[0103] like Then the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The corresponding affine transformation of the hyperspectral resolution remote sensing patch image HSPE-SBK in row q and column p p，q No further matching will be performed.

[0104] like Then the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The corresponding affine transformation of the hyperspectral resolution remote sensing patch image HSPE_SBK in row q and column p p，q The matching process is as follows:

[0105] Corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p p，q The corresponding affine transformation of the hyperspectral resolution remote sensing patch image HSPE_SBK in row q and column p p，q The SIFT+RANSAC matching algorithm was used to obtain the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q Pixel coordinates and corrected hyperspectral resolution remote sensing patch image HSPE_SBK (row q, column p) p，q Pixel coordinates.

[0106] The pixel coordinates of the high spatial resolution remote sensing block image corresponding to the block in row q and column p and the pixel coordinates of the high spectral resolution remote sensing block image corresponding to the affine transformation are merged into multiple pairs of points with the same name in the block in row q and column p.

[0107] Preferably, step 7, which involves constructing a polynomial transformation model for each block based on multiple pairs of corresponding points, is as follows:

[0108] Obtain the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q property

[0109] like Then the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p is directly used. p，q The corresponding affine transformation of the hyperspectral resolution remote sensing patch image HSPE_SBK in row q and column p p，q Matched pairs of points with the same name are used to solve the polynomial model.

[0110] like Sequentially evaluate the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p-1. p-1， q attribute The corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p+1. p+1，q property The corrected high spatial resolution remote sensing patch image HSPA_SBK in row q-1 and column p. p，q-1 ,property The corrected high spatial resolution remote sensing patch image HSPA_SBK in row q+1 and column p.p，q+1 property If the value is 1, then the corresponding block will match multiple pairs of identical points as multiple pairs of identical points in the current sub-block of the qth row and pth column.

[0111] The polynomial transformation model for the block in row q and column p is constructed as follows:

[0112] Obtain the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q Matching the h-th pixel coordinates

[0113]

[0114] h∈[1, CT] p，q ]ph∈[1,nSBK_lon],q h ∈[1, nSBK_lat]

[0115] in, This represents the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The high spectral resolution remote sensing patch image HSPE_SBK with affine transformation of row q and column p p，q The number of matched points with the same name. h represents the sequence number, nSBK_lon is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction.

[0116] Obtain the hyperspectral resolution remote sensing patch image HSPE_SBK with affine transformation at row q and column p. p，q Matching the h-th pixel coordinates

[0117]

[0118] h∈[1, CT] p，q ]p h ∈[1, nSBK_lon], q h ∈[1, nSBK_lat]

[0119] in, This represents the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The high spectral resolution remote sensing patch image HSPE_SBK with affine transformation of row q and column p p，q The number of matched points with the same name. h represents the sequence number, nSBK_lon is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction.

[0120] The corrected high spatial resolution remote sensing patch image HSPA_SBK is sequentially calculated based on the q-th row and p-th column. p，qMatching the h-th pixel coordinates The hyperspectral resolution remote sensing patch image HSPE_SBK with affine transformation of row q and column p p，q Matching the h-th pixel coordinates Construct the h-th quadratic polynomial in row q and column p, thereby establishing the quadratic polynomial model of the entire block in row q and column p, and use least squares to solve for the corresponding quadratic polynomial model parameters.

[0121] Preferably, step 9 involves fusion processing to obtain a fused remote sensing image with high spatial resolution and high spectral resolution, as detailed below:

[0122] Set the weight of the blue band to B. w The weight of the green band is G. w The weight of the red band is R. w The near-infrared band weight is IR. w .

[0123] The pixel values ​​of each corrected high spatial resolution remote sensing block image and the corresponding quadratic polynomial corrected high spectral resolution remote sensing block image are acquired pixel by pixel. The blue-green-red band adjustment value ADJ1 and the near-infrared adjustment value ADJ2 corresponding to each pixel value are calculated as follows:

[0124]

[0125] HSPA DN The values ​​represent the corrected panchromatic pixel values ​​on the high spatial resolution remote sensing image patch. ADJ1 represents the adjustment value for the blue, green, and red bands for each pixel value; ADJ2 represents the near-infrared adjustment value for each pixel value; AVGW1 is the weighted average of the red, green, and blue bands for each pixel; and AVGW2 is the weighted average of the red, green, blue, and near-infrared bands for each pixel. The calculation method is as follows.

[0126]

[0127] Among them, B DN G represents the blue band pixel value in a hyperspectral resolution remote sensing patch image with quadratic polynomial correction. DN R represents the green band pixel value in a quadratic polynomial corrected hyperspectral resolution remote sensing patch image. DN IR represents the red band pixel value in a hyperspectral resolution remote sensing patch image with quadratic polynomial correction. DN This represents the near-infrared pixel value in a hyperspectral resolution remote sensing block image corrected by a quadratic polynomial.

[0128] The image value after fusion is calculated pixel-by-pixel using the blue-green-red band adjustment value ADJ1 corresponding to each pixel value and the near-infrared adjustment value ADJ2 corresponding to each pixel value, as follows:

[0129]

[0130] Among them, B FUS G represents the blue band pixel value of the fused image. FUS Indicates the green band pixel value of the fused image, R FUS Indicates the red band pixel value and IR of the fused image. FUS This represents the pixel value in the near-infrared band of the fused image.

[0131] The present invention also provides a computer-readable medium storing a computer program executed by an electronic device, which, when run on the electronic device, causes the electronic device to perform the steps of the remote sensing image fusion method that takes into account image distortion.

[0132] The advantages of this invention are that, based on an initial registration and block-based fine registration strategy, and by applying polynomial correction to each block individually, it effectively solves the problem of poor fusion results caused by image distortion, thus improving the fusion effect. Furthermore, the block-based approach requires less memory and enables the fusion of wide-swath remote sensing images. Attached Figure Description

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

[0134] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

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

[0136] The following is combined with Figure 1 The technical solution of this invention is a remote sensing image fusion method that takes into account image distortion, as detailed below:

[0137] Step 1: Acquire the original high spatial resolution remote sensing image and its corresponding RPC positioning model, and perform geometric correction to obtain the corrected high spatial resolution remote sensing image; acquire the high spectral resolution remote sensing image and its corresponding RPC positioning model, and perform geometric correction to obtain the corrected high spectral resolution remote sensing image.

[0138] Step 2: Calculate the intersection polygons of the corrected high spatial resolution image and the corrected hyperspectral resolution image, and obtain the geographic extent of the intersection polygons. Construct the minimum bounding rectangle of the intersection polygons based on the geographic extent of the intersection polygons. Perform grid division on the minimum bounding rectangle of the intersection polygons to obtain multiple blocks of minimum bounding rectangles;

[0139] Step 2 calculates the intersection polygon of the corrected high spatial resolution image and the corrected hyperspectral resolution image. The specific process is as follows:

[0140] Construct corrected high spatial resolution image vector polygons based on the geographic coordinates of the four corner points of the corrected high spatial resolution image;

[0141] Construct corrected hyperspectral resolution image vector polygons based on the geographic coordinates of the four corner points of the corrected hyperspectral resolution image;

[0142] The intersection polygon between the corrected high spatial resolution image vector polygon and the corrected high spectral resolution image vector polygon is calculated using a spatial intersection calculation method.

[0143] Step 2 involves obtaining the geographical extent of the intersecting polygons and constructing the minimum bounding rectangle of the intersecting region based on the geographical extent of the intersecting polygons. The specific process is as follows:

[0144] Obtain the longitude and latitude of each point in the intersecting polygon in sequence;

[0145] To construct the longitude range of an intersecting polygon, we select the minimum and maximum longitude values ​​for each point within the polygon from among multiple longitudes.

[0146] [FMIN_LON, FMAX_LON]

[0147] Where FMIN_LON represents the minimum longitude of each point in the intersecting polygon, and FMAX_LON represents the maximum longitude of each point in the intersecting polygon.

[0148] To construct the dimensional range of an intersecting polygon, we select the minimum and maximum dimensional values ​​of each point within the polygon from among multiple points in the polygon, as follows:

[0149] [FMIN_LAT, FMAX_LAT]

[0150] Where FMIN_LAT represents the minimum latitude of each point in the intersecting polygon, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygon;

[0151] Use (FMAX_LON, FMIN_LAT), (FMAX_LON, FMAX_LAT), (FMAX_LON, FMAX_LAT), (FMAX_LON, FMIN_LAT) as the four vertices of the minimum bounding rectangle of the intersecting polygon to construct the minimum bounding rectangle of the intersecting polygon.

[0152] Step 2 involves performing grid division on the minimum bounding rectangle of the intersecting polygons to obtain multiple blocks of minimum bounding rectangles, as detailed below:

[0153] Set the number of grid divisions in the longitude direction to nFBK_lon, and set the number of grid divisions in the latitude direction to nFBK_lat;

[0154] Calculate the longitude width and latitude height of each minimum bounding rectangle block as follows:

[0155]

[0156] Where ΔF_lon is the longitude width of each minimum bounding rectangle block, ΔF_lat is the latitude height of each minimum bounding rectangle block, nFBK_lon = 5 is the number of grid divisions in the longitude direction, nFBK_lat = 5 is the set number of grid divisions in the latitude direction, FMIN_LON represents the minimum longitude of each point in the intersecting polygon, FMAX_LON represents the maximum longitude of each point in the intersecting polygon, and FMIN_LAT represents the minimum latitude of each point in the intersecting polygon, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygon.

[0157] The smallest bounding rectangle block in row j and column i is defined as:

[0158] FBK i，j

[0159] i∈[1,nFBK_lon],j∈[1,nFBK_lat]

[0160] Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_lat is the number of latitudinal grid divisions.

[0161] Calculate the minimum longitude, maximum longitude, minimum latitude, and maximum latitude of the minimum bounding rectangle of the minimum bounding rectangle block in row j and column i, as follows:

[0162]

[0163] in, Let represent the minimum longitude of the smallest bounding rectangle block in row j and column i. This represents the maximum longitude of the smallest bounding rectangular block in row j and column i. Let represent the minimum dimension of the smallest bounding rectangle block in the j-th row and i-th column. This represents the maximum latitude of the smallest bounding rectangle block in the j-th row and i-th column;

[0164] FMIN_LON represents the minimum longitude of each point in the intersecting polygons, and FMAX_LON represents the maximum longitude of each point in the intersecting polygons; where FMIN_LAT represents the minimum latitude of each point in the intersecting polygons, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygons.

[0165] Will As four vertices, construct the geographic rectangle corresponding to the minimum bounding rectangle of the j-th row and i-th column;

[0166] Step 3: Map each block of the minimum bounding rectangle to the corrected high spatial resolution remote sensing image to obtain each corrected high spatial resolution remote sensing block image; map each block of the minimum bounding rectangle to the corrected high spectral resolution remote sensing image to obtain each corrected high spectral resolution remote sensing block image; and perform initial matching processing between each corrected high spatial resolution remote sensing block image and its corresponding corrected high spectral resolution remote sensing block image to obtain multiple pairs of corresponding points for each block.

[0167] Step 3, which maps the image to the corrected high spatial resolution remote sensing image, yields each corrected high spatial resolution remote sensing patch image. The specific process is as follows:

[0168] Divide the four vertices into blocks based on the minimum bounding rectangle in row j and column i. The geographic range is composed of image content obtained from corrected high spatial resolution remote sensing imagery.

[0169] The corrected high spatial resolution remote sensing patch image in row j and column i is defined as:

[0170] HSPA_FBK i，j

[0171] i∈[1,nFBK_lon],j∈[1,nFBK_lat]

[0172] Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_1at is the number of latitude grid divisions.

[0173] Step 3, which maps the image to the corrected hyperspectral resolution remote sensing image, yields each hyperspectral resolution remote sensing block image. The specific process is as follows:

[0174] Divide the four vertices into blocks based on the minimum bounding rectangle in row j and column i. The geographic range is composed of image content obtained from corrected hyperspectral resolution remote sensing images.

[0175] The corrected hyperspectral resolution remote sensing patch image in row j and column i is defined as:

[0176] HSPE_FBK i，j

[0177] i∈[1,nFBK_lon],j∈[1,nFBK_lat]

[0178] Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_lat is the number of latitudinal grid divisions.

[0179] Step 3 involves initial matching to obtain multiple pairs of identically named points for each block. The specific process is as follows:

[0180] Corrected high spatial resolution remote sensing patch image HSPA_FBK in row j and column i i，j and the corresponding corrected hyperspectral resolution remote sensing patch image HSPE_FBK in row j and column i. i，j The SIFT+KNN matching algorithm was used to obtain the corrected high spatial resolution remote sensing patch image HSPA_FBK in the j-th row and i-th column. i，j Pixel coordinates and corrected hyperspectral resolution remote sensing patch image HSPE_FBK in row j and column i i，j Pixel coordinates.

[0181] Get the minimum longitude of the sub-block in row j and column i. minimum latitude The high spatial resolution remote sensing patch image HSPA_FBK in row j and column i is corrected. i，jPixel coordinates are mapped to geographic coordinates of high spatial resolution remote sensing block image by summing the product of pixel coordinates and resolution with the minimum value.

[0182] Get the minimum longitude of the sub-block in row j and column i. minimum latitude The corrected hyperspectral resolution remote sensing patch image HSPE_FBK in row j and column i. i，j Pixel coordinates are used to map hyperspectral resolution remote sensing block images to geographic coordinates by summing the product of pixel coordinates and resolution with the minimum value.

[0183] Merge the geographic coordinates of the high spatial resolution remote sensing block image and the high spectral resolution remote sensing block image corresponding to the sub-block in row j and column i into multiple pairs of points with the same name in row j and column i.

[0184] Step 4: Construct an affine transformation model based on multiple sets of corresponding point pairs in multiple blocks, and perform affine transformation processing on the corrected hyperspectral resolution image according to the affine transformation model to obtain an affine transformed hyperspectral resolution image.

[0185] Step 4 describes the construction of an affine transformation model based on multiple sets of corresponding point pairs across multiple blocks. The specific process is as follows:

[0186] A global affine transformation model is constructed based on the geographic coordinates of the corresponding high spatial resolution remote sensing block image and the high spectral resolution remote sensing block image in all block-named point pairs, and the parameters of the affine transformation model are solved using the least squares method.

[0187] Step 4 involves performing an affine transformation on the corrected hyperspectral resolution image according to an affine transformation model. The specific process is as follows:

[0188] By combining the affine transformation model parameters (e0, e1, e2, f0, f1, f2), an indirect correction method is used to perform affine transformation correction on the corrected hyperspectral resolution image, resulting in an affine transformed hyperspectral resolution image.

[0189] Step 5: Calculate the updated intersecting polygons of the corrected high spatial resolution image and the affine-transformed high spectral resolution image, and obtain the geographical extent of the updated intersecting polygons. Construct the minimum bounding rectangle of the updated intersecting polygons based on the geographical extent of the updated intersecting polygons. Perform grid division on the minimum bounding rectangle of the updated intersecting polygons to obtain multiple blocks of minimum bounding rectangles;

[0190] Step 5 involves calculating the intersecting polygon of the corrected high spatial resolution image and the affine-transformed high spectral resolution image after updating. The specific process is as follows:

[0191] A corrected high spatial resolution image vector polygon is constructed based on the geographic coordinates of the four corner points of the corrected high spatial resolution image. An affine transformed high spectral resolution image vector polygon is then constructed based on the geographic coordinates of the four corner points of the affine transformed high spectral resolution image. The intersection polygon between the corrected high spatial resolution image vector polygon and the affine transformed high spectral resolution image vector polygon is calculated using a spatial intersection calculation method; this is the updated intersection polygon.

[0192] Step 5 involves obtaining the geographical extent of the updated intersecting polygons and constructing the minimum bounding rectangle of the intersecting region based on the updated geographical extent of the intersecting polygons. The specific process is as follows:

[0193] Sequentially obtain the longitude and latitude of each point in the updated intersecting polygon;

[0194] To construct the longitude range of the updated intersecting polygon, we filter out the minimum and maximum longitude values ​​of each point within the updated polygon from the longitudes of multiple points.

[0195] [SMIN_LON, SMAX_LON]

[0196] Where SMIN_LON represents the minimum longitude of each point in the updated intersecting polygon, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygon;

[0197] In the updated intersecting polygon, the minimum and maximum latitude values ​​of each point are selected from the latitudes of multiple points in the updated polygon. The latitude range of the updated intersecting polygon is then constructed, defined as follows:

[0198] [SMIN_LAT, SMAX_LAT]

[0199] Where SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygon, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygon;

[0200] Use (SMAX_LON, SMIN_LAT), (SMAX_LON, SMAX_LAT), (SMAX_LON, SMAX_LAT), (SMAX_LON, SMIN_LAT) as the four vertices of the minimum bounding rectangle of the updated intersecting polygon to construct the minimum bounding rectangle of the updated intersecting polygon.

[0201] As described in step 5, the minimum bounding rectangle of the updated intersecting polygons is divided into multiple blocks of minimum bounding rectangles by grid partitioning, as follows:

[0202] Set the number of grid divisions in the longitude direction to nSBK_lon, and set the number of grid divisions in the latitude direction to nSBK_lat;

[0203] Calculate the longitude width and latitude height of each minimum bounding rectangle block as follows:

[0204]

[0205] Where ΔS_lon is the longitude width of each minimum bounding rectangle block, ΔS_lat is the latitude height of each minimum bounding rectangle block, nSBK_lon = 5 is the number of grid divisions in the longitude direction, nSBK_lat = 5 is the set number of grid divisions in the latitude direction, SMIN_LON represents the minimum longitude of each point in the updated intersecting polygon, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygon; where SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygon, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygon.

[0206] The smallest bounding rectangle block in row q and column p is defined as:

[0207] SBK p，q

[0208] p∈[1,nSBK_lon], q∈[1,nSBK_lat]

[0209] Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_lon is the number of directional grid divisions, and nSBK_lat is the number of latitudinal grid divisions.

[0210] Calculate the minimum longitude, maximum longitude, minimum latitude, and maximum latitude of the minimum bounding rectangle of the minimum bounding rectangle of the q-th row and p-th column, as follows:

[0211]

[0212] in, This represents the minimum longitude of the smallest bounding rectangle block in row q and column p. This represents the maximum longitude of the smallest bounding rectangle in row q and column p. This represents the minimum latitude of the smallest bounding rectangle block in row q and column p. This represents the maximum latitude of the smallest bounding rectangle in row q and column p.

[0213] SMIN_LON represents the minimum longitude of each point in the updated intersecting polygon, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygon; where SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygon, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygon; As four vertices, construct the geographic rectangular block corresponding to the minimum bounding rectangular block in row q and column p;

[0214] Step 6: Map each block of the minimum bounding rectangle to the corrected high spatial resolution remote sensing image to obtain each corrected high spatial resolution remote sensing block image; map each block of the minimum bounding rectangle to the affine transformed high spectral resolution remote sensing image to obtain each corrected high spectral resolution remote sensing block image; and perform fine matching processing between each corrected high spatial resolution remote sensing block image and the corresponding affine transformed high spectral resolution remote sensing block image to obtain multiple pairs of corresponding points for each block.

[0215] Step 6, which maps the image to the corrected high spatial resolution remote sensing image, yields each corrected high spatial resolution remote sensing patch image. The specific process is as follows:

[0216] Divide the four vertices into blocks based on the minimum bounding rectangle in row q and column p. The geographic range is composed of image content obtained from corrected high spatial resolution remote sensing imagery.

[0217] The corrected high spatial resolution remote sensing patch image in row q and column p is defined as:

[0218] HSPA_SBK p，q

[0219] p∈[1,nSBK_lon], q∈[1,nSBK_lat]

[0220] Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_ion is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction.

[0221] Step 6 describes mapping the hyperspectral resolution remote sensing image to an affine transformation to obtain each corrected high spatial resolution remote sensing block image. The specific process is as follows:

[0222] Divide the four vertices into blocks based on the minimum bounding rectangle in row q and column p. The geographic range is composed of image content obtained from corrected high spatial resolution remote sensing imagery.

[0223] The high spatial resolution remote sensing patch image of the affine transformation in row q and column p is defined as:

[0224] HSPE_SBK p，q

[0225] p∈[1,nSBK_lon], q∈[1,nSBK_lat]

[0226] Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_lon is the number of directional grid divisions, and nSBK_lat is the number of latitude grid divisions.

[0227] Step 6 describes obtaining multiple pairs of identical points for each block through fine-grained matching. The specific process is as follows:

[0228] Corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p p，q The validity determination is as follows:

[0229] If CNT_HSPA_SBK p，q ≥α, then

[0230] If CNT_HSPA_SBK p，q <α, then

[0231] Among them, CNT_HSPA_SBK p，q This represents the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The number of pixels with values ​​greater than 0, where α represents the threshold for determining the validity of the block (α = 50). This represents the validity attribute of the corrected high spatial resolution remote sensing block image in row q and column p.

[0232] like Then the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The corresponding affine transformation of the hyperspectral resolution remote sensing patch image HSPE_SBK in row q and column p p，q No further matching will be performed.

[0233] like Then the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The corresponding affine transformation of the hyperspectral resolution remote sensing patch image HSPE_SBK in row q and column p p，q The matching process is as follows:

[0234] Corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p p，q The corresponding affine transformation of the hyperspectral resolution remote sensing patch image HSPE_SBK in row q and column p p，q The SIFT+RANSAC matching algorithm was used to obtain the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q Pixel coordinates and corrected hyperspectral resolution remote sensing patch image HSPE_SBK (row q, column p) p，q Pixel coordinates.

[0235] The pixel coordinates of the high spatial resolution remote sensing block image corresponding to the block in row q and column p and the pixel coordinates of the high spectral resolution remote sensing block image corresponding to the affine transformation are merged into multiple pairs of points with the same name in the block in row q and column p.

[0236] Step 7: Construct the polynomial transformation model of each block based on multiple pairs of corresponding points in each block. Combine the multiple pairs of corresponding points in each sub-block and use the least squares method to solve for the parameters of the polynomial transformation model of each sub-block, thus obtaining the polynomial transformation model of each block.

[0237] Step 7 describes the construction of a polynomial transformation model for each block based on multiple pairs of corresponding points, as follows:

[0238] Obtain the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q property

[0239] like Then the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p is directly used. p，q The corresponding affine transformation of the hyperspectral resolution remote sensing patch image HSPE_SBK in row q and column p p，q Matched pairs of points with the same name are used to solve the polynomial model.

[0240] like Sequentially evaluate the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p-1. p-1，q property The corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p+1. p+1，q property The corrected high spatial resolution remote sensing patch image HSPA_SBK in row q-1 and column p. p，q-1 ,property The corrected high spatial resolution remote sensing patch image HSPA_SBK in row q+1 and column p. p，q+1 property If the value is 1, then the corresponding block will match multiple pairs of identical points as multiple pairs of identical points in the current sub-block of the qth row and pth column.

[0241] The polynomial transformation model for the block in row q and column p is constructed as follows:

[0242] Obtain the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q Matching the h-th pixel coordinates

[0243]

[0244] h∈[1, CT] p，q ]ph∈[1,nSBK_lon], qh∈[1,nSBK_lat]

[0245] in, This represents the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The high spectral resolution remote sensing patch image HSPE_SBK with affine transformation of row q and column p p，q The number of matched points with the same name. h represents the sequence number, nSBK_ion is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction.

[0246] Obtain the hyperspectral resolution remote sensing patch image HSPE_SBK with affine transformation at row q and column p. p，q Matching the h-th pixel coordinates

[0247]

[0248] h∈[1, CT] p，q ]p h ∈[1, nSBK_lon], q h ∈[1, nSBK_lat]

[0249] in, This represents the corrected high spatial resolution remote sensing patch image HSPA_SBK in row q and column p. p，q The high spectral resolution remote sensing patch image HSPE_SBK with affine transformation of row q and column p p，q The number of matched points with the same name. h represents the sequence number, nSBK_lon is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction.

[0250] The corrected high spatial resolution remote sensing patch image HSPA_SBK is sequentially calculated based on the q-th row and p-th column. p，q Matching the h-th pixel coordinates The hyperspectral resolution remote sensing patch image HSPE_SBK with affine transformation of row q and column p p，q Matching the h-th pixel coordinates Construct the h-th quadratic polynomial in row q and column p, thereby establishing the quadratic polynomial model of the entire block in row q and column p, and use least squares to solve for the corresponding quadratic polynomial model parameters.

[0251] Step 8: Perform quadratic polynomial correction on the affine-transformed hyperspectral resolution remote sensing image based on the quadratic polynomial transformation model of each block to obtain the hyperspectral resolution remote sensing block image of each quadratic polynomial correction.

[0252] Step 9: Perform fusion processing on each corrected high spatial resolution remote sensing block image and the corresponding quadratic polynomial corrected high spectral resolution remote sensing block image to obtain a fused high spatial resolution and high spectral resolution remote sensing image.

[0253] Step 9 describes the fusion process to obtain a fused remote sensing image with high spatial resolution and high spectral resolution, as detailed below:

[0254] Set the weight of the blue band to B. w =0.75, green band weight is G w =0.85, red band weight is R w =0.6, the near-infrared band weight is IR w =0.3. Pixel values ​​of each corrected high spatial resolution remote sensing block image and its corresponding quadratic polynomial corrected high spectral resolution remote sensing block image are acquired pixel by pixel. The blue-green-red band adjustment value ADJ1 and the near-infrared adjustment value ADJ2 corresponding to each pixel value are calculated as follows:

[0255]

[0256] HSPA DN The values ​​represent the corrected panchromatic pixel values ​​on the high spatial resolution remote sensing image patch. ADJ1 represents the adjustment value for the blue, green, and red bands for each pixel value; ADJ2 represents the near-infrared adjustment value for each pixel value; AVGw1 is the weighted average of the red, green, and blue bands for each pixel; and AVGW2 is the weighted average of the red, green, blue, and near-infrared bands for each pixel. The calculation method is as follows.

[0257]

[0258] Among them, B DN G represents the blue band pixel value in a hyperspectral resolution remote sensing patch image with quadratic polynomial correction. DNR represents the green band pixel value in a quadratic polynomial corrected hyperspectral resolution remote sensing patch image. DN IR represents the red band pixel value in a hyperspectral resolution remote sensing patch image with quadratic polynomial correction. DN This represents the near-infrared pixel value in a hyperspectral resolution remote sensing block image corrected by a quadratic polynomial.

[0259] The image value after fusion is calculated pixel-by-pixel using the blue-green-red band adjustment value ADJ1 corresponding to each pixel value and the near-infrared adjustment value ADJ2 corresponding to each pixel value, as follows:

[0260]

[0261] Among them, B FUs G represents the blue band pixel value of the fused image. FUs Indicates the green band pixel value of the fused image, R FUs Indicates the red band pixel value and IR of the fused image. FUS This represents the pixel value in the near-infrared band of the fused image.

[0262] A specific embodiment of the present invention also provides a computer-readable medium.

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

[0264] The server workstation stores the computer program executed by the electronic device. When the computer program runs on the electronic device, it causes the electronic device to perform the steps of the remote sensing image fusion method that takes into account image distortion according to the embodiments of the present invention.

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

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

Claims

1. A remote sensing image fusion method that takes into account image distortion, characterized in that, Includes the following steps: Step 1: Acquire the original high spatial resolution remote sensing image and its corresponding RPC positioning model, and perform geometric correction to obtain the corrected high spatial resolution remote sensing image; acquire the high spectral resolution remote sensing image and its corresponding RPC positioning model, and perform geometric correction to obtain the corrected high spectral resolution remote sensing image. Step 2: Calculate the intersection polygons of the corrected high spatial resolution image and the corrected high spectral resolution image and obtain the geographic extent of the intersection polygons. Construct the minimum bounding rectangle of the intersection polygons based on the geographic extent of the intersection polygons. The minimum bounding rectangle of intersecting polygons is divided into grids to obtain multiple blocks of minimum bounding rectangles; Step 3: Map each block of the minimum bounding rectangle to the corrected high spatial resolution remote sensing image to obtain each corrected high spatial resolution remote sensing block image; map each block of the minimum bounding rectangle to the corrected high spectral resolution remote sensing image to obtain each corrected high spectral resolution remote sensing block image; and perform initial matching processing between each corrected high spatial resolution remote sensing block image and its corresponding corrected high spectral resolution remote sensing block image to obtain multiple pairs of corresponding points for each block. Step 4: Construct an affine transformation model based on multiple sets of corresponding point pairs in multiple blocks, and perform affine transformation processing on the corrected hyperspectral resolution image according to the affine transformation model to obtain an affine transformed hyperspectral resolution image. Step 5: Calculate the updated intersecting polygons of the corrected high spatial resolution image and the affine transformed high spectral resolution image, and obtain the geographical extent of the updated intersecting polygons. Construct the minimum bounding rectangle of the updated intersecting polygons based on the geographical extent of the updated intersecting polygons. The minimum bounding rectangle of the updated intersecting polygons is divided into multiple blocks of minimum bounding rectangles by performing grid division. Step 6: Map each block of the minimum bounding rectangle to the corrected high spatial resolution remote sensing image to obtain each corrected high spatial resolution remote sensing block image; map each block of the minimum bounding rectangle to the affine transformed high spectral resolution remote sensing image to obtain each corrected high spectral resolution remote sensing block image; and perform fine matching processing between each corrected high spatial resolution remote sensing block image and the corresponding affine transformed high spectral resolution remote sensing block image to obtain multiple pairs of corresponding points for each block. Step 7: Construct the polynomial transformation model of each block based on multiple pairs of corresponding points in each block. Combine the multiple pairs of corresponding points in each sub-block and use the least squares method to solve for the parameters of the polynomial transformation model of each sub-block, thus obtaining the polynomial transformation model of each block. Step 8: Perform quadratic polynomial correction on the affine-transformed hyperspectral resolution remote sensing image based on the quadratic polynomial transformation model of each block to obtain the hyperspectral resolution remote sensing block image of each quadratic polynomial correction. Step 9: Perform fusion processing on each corrected high spatial resolution remote sensing block image and the corresponding quadratic polynomial corrected high spectral resolution remote sensing block image to obtain a fused high spatial resolution and high spectral resolution remote sensing image. Step 7, which describes constructing the polynomial transformation model for each block based on multiple pairs of corresponding points, is as follows: Obtain the corrected high spatial resolution remote sensing patch image in row q and column p. property , like Then, the corrected high spatial resolution remote sensing block image in row q and column p is directly used. The hyperspectral resolution remote sensing patch image with the corresponding affine transformation of row q and column p. Matched pairs of points with the same name are used to solve the polynomial model; like The corrected high spatial resolution remote sensing patch image in row q and column p-1 is judged sequentially. property Corrected high spatial resolution remote sensing patch image in row q and column p+1 property Corrected high spatial resolution remote sensing patch image in row q-1 and column p ,property Corrected high spatial resolution remote sensing patch image in row q+1 and column p property If the value is 1, then the corresponding block will match multiple pairs of identical points as multiple pairs of identical points in the current sub-block of the qth row and pth column. The polynomial transformation model for the block in row q and column p is constructed as follows: Obtain the corrected high spatial resolution remote sensing patch image in row q and column p. Matching the h-th pixel coordinates in, This represents the corrected high spatial resolution remote sensing patch image in row q and column p. The hyperspectral resolution remote sensing patch image of row q and column p of affine transformation The number of matching points with the same name; Indicates the sequence number, nSBK_lon is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction; Obtain the hyperspectral resolution remote sensing patch image of the affine transformation of the q-th row and p-th column. Matching the h-th pixel coordinates in, This represents the corrected high spatial resolution remote sensing patch image in row q and column p. The hyperspectral resolution remote sensing patch image of row q and column p of affine transformation The number of matching points with the same name; Indicates the sequence number, nSBK_lon is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction; The high spatial resolution remote sensing patch image is corrected sequentially based on the q-th row and p-th column. Matching the h-th pixel coordinates The hyperspectral resolution remote sensing patch image with affine transformation of row q and column p Matching the h-th pixel coordinates Construct the h-th quadratic polynomial in row q and column p, thereby establishing the quadratic polynomial model of the entire block in row q and column p, and use least squares to solve for the corresponding quadratic polynomial model parameters.

2. The remote sensing image fusion method considering image distortion according to claim 1, characterized in that: Step 2 calculates the intersection polygon of the corrected high spatial resolution image and the corrected hyperspectral resolution image. The specific process is as follows: Construct corrected high spatial resolution image vector polygons based on the geographic coordinates of the four corner points of the corrected high spatial resolution image; Construct corrected hyperspectral resolution image vector polygons based on the geographic coordinates of the four corner points of the corrected hyperspectral resolution image; The intersection polygon between the corrected high spatial resolution image vector polygon and the corrected high spectral resolution image vector polygon is calculated using a spatial intersection calculation method. Step 2 involves obtaining the geographical extent of the intersecting polygons and constructing the minimum bounding rectangle of the intersecting region based on the geographical extent of the intersecting polygons. The specific process is as follows: Obtain the longitude and latitude of each point in the intersecting polygon in sequence; To construct the longitude range of an intersecting polygon, we select the minimum and maximum longitude values ​​for each point within the polygon from among multiple longitudes. [FMIN_LON, FMAX_LON] Where FMIN_LON represents the minimum longitude of each point in the intersecting polygon, and FMAX_LON represents the maximum longitude of each point in the intersecting polygon; To construct the dimensional range of an intersecting polygon, we select the minimum and maximum dimensional values ​​of each point within the polygon from among multiple points in the polygon, as follows: [FMIN_LAT, FMAX_LAT] Where FMIN_LAT represents the minimum latitude of each point in the intersecting polygon, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygon; Use (FMAX_LON, FMIN_LAT), (FMAX_LON, FMAX_LAT), (FMAX_LON, FMAX_LAT), (FMAX_LON, FMIN_LAT) as the four vertices of the minimum bounding rectangle of the intersecting polygon to construct the minimum bounding rectangle of the intersecting polygon; Step 2 involves performing grid division on the minimum bounding rectangle of the intersecting polygons to obtain multiple blocks of minimum bounding rectangles, as detailed below: Set the number of grid divisions in the longitude direction to nFBK_lon, and set the number of grid divisions in the latitude direction to nFBK_lat; Calculate the longitude width and latitude height of each minimum bounding rectangle block as follows: in, The longitude and width of each minimum bounding rectangle block. The latitude and height of each minimum bounding rectangle block are defined by nFBK_lon, the number of grid divisions in the longitude direction, and nFBK_lat, which sets the number of grid divisions in the latitude direction. FMIN_LON represents the minimum longitude of each point in the intersecting polygons, and FMAX_LON represents the maximum longitude of each point in the intersecting polygons. In addition, FMIN_LAT represents the minimum latitude of each point in the intersecting polygons, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygons. The smallest bounding rectangle block in row j and column i is defined as: ， Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_lat is the number of latitudinal grid divisions. Calculate the minimum longitude, maximum longitude, minimum latitude, and maximum latitude of the minimum bounding rectangle of the minimum bounding rectangle block in row j and column i, as follows: in, Let represent the minimum longitude of the smallest bounding rectangle block in row j and column i. This represents the maximum longitude of the smallest bounding rectangular block in row j and column i. This represents the minimum latitude of the smallest bounding rectangle block in the j-th row and i-th column. This represents the maximum latitude of the smallest bounding rectangle block in the j-th row and i-th column; FMIN_LON represents the minimum longitude of each point in the intersecting polygons, and FMAX_LON represents the maximum longitude of each point in the intersecting polygons; where FMIN_LAT represents the minimum latitude of each point in the intersecting polygons, and FMAX_LAT represents the maximum latitude of each point in the intersecting polygons. Will( , ), ( , ), ( , ), , ) are used as four vertices to construct the geographic rectangle corresponding to the minimum bounding rectangle of the j-th row and i-th column.

3. The remote sensing image fusion method considering image distortion according to claim 2, characterized in that: Step 3, which maps the image to the corrected high spatial resolution remote sensing image, yields each corrected high spatial resolution remote sensing patch image. The specific process is as follows: Divide the blocks into four vertices based on the minimum bounding rectangle in row j and column i. , ), ( , ), ( , ), , The geographic range is composed of corrected high spatial resolution remote sensing images, and the corresponding geographic range image content is obtained from the geographic range. The corrected high spatial resolution remote sensing patch image in row j and column i is defined as: ， Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_lat is the number of latitudinal grid divisions. Step 3, which maps the image to the corrected hyperspectral resolution remote sensing image, yields each hyperspectral resolution remote sensing block image. The specific process is as follows: Divide the blocks into four vertices based on the minimum bounding rectangle in row j and column i. , ), ( , ), ( , ), , The geographic range is composed of corrected hyperspectral resolution remote sensing images, and the corresponding geographic range image content is obtained from the geographic range. The corrected hyperspectral resolution remote sensing patch image in row j and column i is defined as: ， Where i represents the block number in the longitude direction, j represents the block number in the latitude direction, nFBK_lon is the number of directional grid divisions, and nFBK_lat is the number of latitudinal grid divisions. Step 3 involves initial matching to obtain multiple pairs of identically named points for each block. The specific process is as follows: Corrected high spatial resolution remote sensing patch image of row j and column i and the corresponding corrected hyperspectral resolution remote sensing patch image in row j and column i. The SIFT+KNN matching algorithm was used to obtain the corrected high spatial resolution remote sensing patch image in the j-th row and i-th column. Pixel coordinates and corrected hyperspectral resolution remote sensing patch image in row j and column i Pixel coordinates; Get the minimum longitude of the sub-block in row j and column i. minimum latitude The high spatial resolution remote sensing block image in row j and column i is corrected. Pixel coordinates are mapped to geographic coordinates of high spatial resolution remote sensing block image by summing the product of pixel coordinates and resolution with the minimum value. Get the minimum longitude of the sub-block in row j and column i. minimum latitude The corrected hyperspectral resolution remote sensing block image in row j and column i is... Pixel coordinates are used to map hyperspectral resolution remote sensing block images to geographic coordinates by summing the product of pixel coordinates and resolution with the minimum value. Merge the geographic coordinates of the high spatial resolution remote sensing block image and the high spectral resolution remote sensing block image corresponding to the sub-block in row j and column i into multiple pairs of points with the same name in row j and column i.

4. The remote sensing image fusion method considering image distortion according to claim 3, characterized in that: Step 4 describes the construction of an affine transformation model based on multiple sets of corresponding point pairs across multiple blocks. The specific process is as follows: A global affine transformation model is constructed based on the geographic coordinates of the corresponding high spatial resolution remote sensing block image and the high spectral resolution remote sensing block image in all block-named point pairs, and the parameters of the affine transformation model are solved by the least squares method. Step 4 involves performing an affine transformation on the corrected hyperspectral resolution image according to an affine transformation model. The specific process is as follows: Combined with affine transformation model parameters An indirect correction method is used to perform affine transformation correction on the corrected hyperspectral resolution image to obtain an affine transformed hyperspectral resolution image.

5. The remote sensing image fusion method considering image distortion according to claim 4, characterized in that: Step 5 involves calculating the intersecting polygon of the corrected high spatial resolution image and the affine-transformed high spectral resolution image after updating. The specific process is as follows: Construct a corrected high spatial resolution image vector polygon based on the geographic coordinates of the four corner points of the corrected high spatial resolution image; construct an affine transformed high spectral resolution image vector polygon based on the geographic coordinates of the four corner points of the affine transformed high spectral resolution image; calculate the intersection polygon between the corrected high spatial resolution image vector polygon and the affine transformed high spectral resolution image vector polygon using a spatial intersection calculation method, which is the updated intersection polygon. Step 5 involves obtaining the geographical extent of the updated intersecting polygons and constructing the minimum bounding rectangle of the intersecting region based on the updated geographical extent of the intersecting polygons. The specific process is as follows: Sequentially obtain the longitude and latitude of each point in the updated intersecting polygon; To construct the longitude range of the updated intersecting polygon, we filter out the minimum and maximum longitude values ​​of each point within the updated polygon from the longitudes of multiple points. [SMIN_LON, SMAX_LON] Where SMIN_LON represents the minimum longitude of each point in the updated intersecting polygon, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygon; In the updated intersecting polygon, the minimum and maximum latitude values ​​of each point are selected from the latitudes of multiple points in the updated polygon. The latitude range of the updated intersecting polygon is then constructed, defined as follows: [SMIN_LAT, SMAX_LAT] Where SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygon, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygon; Use (SMAX_LON, SMIN_LAT), (SMAX_LON, SMAX_LAT), (SMAX_LON, SMAX_LAT), (SMAX_LON, SMIN_LAT) as the four vertices of the minimum bounding rectangle of the updated intersecting polygon to construct the minimum bounding rectangle of the updated intersecting polygon. Step 5 involves performing grid division on the minimum bounding rectangle of the updated intersecting polygons to obtain multiple blocks of minimum bounding rectangles, as detailed below: Set the number of grid divisions in the longitude direction to nSBK_lon, and set the number of grid divisions in the latitude direction to nSBK_lat; Calculate the longitude width and latitude height of each minimum bounding rectangle block as follows: in, The longitude and width of each minimum bounding rectangle block. The latitude and height of each minimum bounding rectangle block are defined by nSBK_lon, the number of grid divisions in the longitude direction, and nSBK_lat, which sets the number of grid divisions in the latitude direction. SMIN_LON represents the minimum longitude of each point in the updated intersecting polygons, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygons. SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygons, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygons. The smallest bounding rectangle block in row q and column p is defined as: ， Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_lon is the number of directional grid divisions, and nSBK_lat is the number of latitudinal grid divisions. Calculate the minimum longitude, maximum longitude, minimum latitude, and maximum latitude of the minimum bounding rectangle of the minimum bounding rectangle of the q-th row and p-th column, as follows: in, This represents the minimum longitude of the smallest bounding rectangle block in row q and column p. This represents the maximum longitude of the smallest bounding rectangle in row q and column p. This represents the minimum latitude of the smallest bounding rectangle block in row q and column p. This represents the maximum latitude of the smallest bounding rectangle in row q and column p. SMIN_LON represents the minimum longitude of each point in the updated intersecting polygon, and SMAX_LON represents the maximum longitude of each point in the updated intersecting polygon; where SMIN_LAT represents the minimum latitude of each point in the updated intersecting polygon, and SMAX_LAT represents the maximum latitude of each point in the updated intersecting polygon; ( , ), ( , ), ( , ), , ) are used as four vertices to construct the geographic rectangle corresponding to the smallest bounding rectangle of the q-th row and p-th column.

6. The remote sensing image fusion method considering image distortion according to claim 5, characterized in that: Step 6, which maps the image to the corrected high spatial resolution remote sensing image, yields each corrected high spatial resolution remote sensing patch image. The specific process is as follows: Divide the blocks into four vertices based on the smallest bounding rectangle in row q and column p. , ), ( , ), ( , ), , The geographic range is composed of corrected high spatial resolution remote sensing images, and the corresponding geographic range image content is obtained from the geographic range. The corrected high spatial resolution remote sensing patch image in row q and column p is defined as: ， Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_lon is the number of grid divisions in the longitude direction, and nSBK_lat is the number of grid divisions in the latitude direction; Step 6 describes mapping the hyperspectral resolution remote sensing image to an affine transformation to obtain each corrected high spatial resolution remote sensing block image. The specific process is as follows: Divide the blocks into four vertices based on the smallest bounding rectangle in row q and column p. , ), ( , ), ( , ), , The geographic range is composed of corrected high spatial resolution remote sensing images, and the corresponding geographic range image content is obtained from the geographic range. The high spatial resolution remote sensing patch image of the affine transformation in row q and column p is defined as: ， Where p represents the block number in the longitude direction, q represents the block number in the latitude direction, nSBK_lon is the number of directional grid divisions, and nSBK_lat is the number of latitude grid divisions.

7. The remote sensing image fusion method considering image distortion according to claim 6, characterized in that: Step 6 describes obtaining multiple pairs of identical points for each block through fine-grained matching. The specific process is as follows: Corrected high spatial resolution remote sensing patch image of row q and column p The validity determination is as follows: like ,but ; like ,but 0; in, This represents the corrected high spatial resolution remote sensing patch image in row q and column p. The number of pixels with values ​​greater than 0. This represents the threshold for determining the validity of the block. This represents the validity attribute of the corrected high spatial resolution remote sensing block image in row q and column p. like Then the corrected high spatial resolution remote sensing patch image in row q and column p is... The hyperspectral resolution remote sensing patch image with the corresponding affine transformation of row q and column p. No further matching will be performed; like Then the corrected high spatial resolution remote sensing patch image in row q and column p is... The hyperspectral resolution remote sensing patch image with the corresponding affine transformation of row q and column p. The matching process is as follows: Corrected high spatial resolution remote sensing patch image of row q and column p The hyperspectral resolution remote sensing patch image with the corresponding affine transformation of row q and column p. The SIFT+RANSAC matching algorithm was used to obtain the corrected high spatial resolution remote sensing patch image at row q and column p. Pixel coordinates and corrected hyperspectral resolution remote sensing patch image at row q and column p Pixel coordinates; The pixel coordinates of the high spatial resolution remote sensing block image corresponding to the block in row q and column p and the pixel coordinates of the high spectral resolution remote sensing block image corresponding to the affine transformation are merged into multiple pairs of points with the same name in the block in row q and column p.

8. The remote sensing image fusion method considering image distortion according to claim 7, characterized in that: Step 9 describes the fusion process to obtain a fused remote sensing image with high spatial resolution and high spectral resolution, as detailed below: Set the weight of the blue band to The weight of the green band is The red band weights are Near-infrared band weights are ; Pixel values ​​of each corrected high spatial resolution remote sensing block image and its corresponding quadratic polynomial corrected high spectral resolution remote sensing block image are acquired pixel by pixel, and the blue, green, and red band adjustment values ​​corresponding to each pixel value are calculated. Near-infrared adjustment value corresponding to each pixel value The calculation method is as follows: in, This represents the corrected panchromatic pixel value on the high spatial resolution remote sensing patch image. Adjustment values ​​for the blue, green, and red bands corresponding to each pixel value. Near-infrared adjustment value corresponding to each pixel value The weighted average of the red, green, and blue bands corresponding to each pixel. The weighted average of the red, green, blue, and near-infrared bands for each pixel is calculated as follows: in, This represents the blue band pixel value in a hyperspectral resolution remote sensing patch image with quadratic polynomial correction. This represents the green band pixel value in a hyperspectral resolution remote sensing patch image with quadratic polynomial correction. This represents the red band pixel value in a hyperspectral resolution remote sensing patch image with quadratic polynomial correction. This represents the near-infrared band pixel value in a hyperspectral resolution remote sensing block image with quadratic polynomial correction. The adjustment values ​​for the blue, green, and red bands corresponding to each pixel value Near-infrared adjustment value corresponding to each pixel value The fused image values ​​are calculated pixel-by-pixel, as follows: in, This represents the blue band pixel value of the fused image. Indicates the green band pixel value of the fused image. Indicates the red band pixel value of the fused image. This represents the pixel value in the near-infrared band of the fused image.

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

Images