A satellite image tremor error detection and correction method based on infinite filter fusion

Abstract

The invention relates to the field of satellite remote sensing data processing technology, and in particular to a satellite image tremor error detection and correction method based on infinite filtering fusion, comprising the following steps: S1, band separation of multispectral images: band separation of original multispectral images; S2, image enhancement processing: Wallis filtering image enhancement processing of multispectral images of each band; S3, dense matching and registration error acquisition of multispectral images: dense matching to obtain relative registration error, and dense ground control point matching to obtain absolute registration error; S4, H-infinity filtering fusion: selecting a combination to obtain a fusion result; S5, bicubic interpolation image compensation: using the fusion result to compensate the image by using a bicubic interpolation method. The invention makes full use of their respective advantages and avoids their disadvantages, thereby realizing more accurate calculation and estimation of the absolute jitter error of the satellite platform tremor.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite remote sensing data processing, and in particular to a satellite image tremor error detection and correction method based on infinite filter fusion. Background Art

[0002] Satellite platform tremor, or satellite jitter, has become one of the most important factors affecting the geometric quality of high-resolution optical satellite images. It is usually caused by dynamic structures on the satellite platform, temperature changes, orbit control operations, etc. Satellite jitter has been observed on many high-resolution satellites, including QuickBird, Beijing-1, ALOS, Pleiades, ZiYuan-02C (ZY1-02C) and ZiYuan-3 (ZY-3). During the early flight data collection of Ziyuan-3, it was found that the satellite jitter was 0.6-0.7Hz, resulting in internal periodic distortion of the three-line array camera image and the multispectral camera image.

[0003] In order to solve the problem caused by satellite platform tremor, the first step is to measure the platform tremor. How to accurately detect the tremor of the satellite platform becomes the key to solving the problem. The platform tremor of the satellite is usually a periodic motion or a motion composed of multiple periods and frequencies superimposed on the attitude pointing error. In order to reconstruct the jitter, the periodic parameters should be determined, including the frequency, amplitude and phase of the relative time-varying errors in two directions. The methods for detecting tremor of the satellite platform are mainly divided into two categories. One is based on attitude sensors to measure attitude fluctuations, and the other is based on the use of images, through parallax observation images (such as multispectral images, interlaced CCD (charge-coupled device) images and stereo image pairs), orthophotos or linear objects in images, to obtain the distortion caused by satellite jitter.

[0004] Since the first category is too dependent on high-performance hardware equipment, most researchers currently mainly use satellite images and their products to detect and estimate platform tremors. However, a single detection method is difficult to detect accurate results. The method using parallax images for detection is too dependent on model construction. Although it can detect accurate relative registration errors, the phase obtained through model construction is uncertain. The tremor detection method using high-precision reference images DOM can directly detect absolute registration errors and obtain relatively accurate phases. However, due to the limitations of fitting accuracy, this method has certain limitations in the accurate determination of amplitude. Studies have shown that the accuracy of satellite platform tremor detection not only depends on the accuracy of the amplitude, but the estimation of the phase is also crucial. Even if the amplitude measurement is accurate, if there is an error in the phase, the final detection result will be inaccurate. Summary of the invention

[0005] Based on the above purpose, the present invention provides a satellite image tremor error detection and correction method based on infinite filter fusion.

[0006] The satellite image tremor error detection and correction method based on infinite filter fusion includes the following steps:

[0007] S1, band separation of multispectral images: band separation of the original multispectral images to obtain multispectral images of each band;

[0008] S2, image enhancement processing: perform wallis filter image enhancement processing on the multispectral images of each band to obtain enhanced multispectral images;

[0009] S3, dense matching and registration error acquisition of multispectral images: Dense matching and registration error acquisition of multispectral images: Dense matching is performed between the bands of the enhanced multispectral images. Based on the normalized cross correlation coefficient method (NCC) matching and least squares matching, the nonlinear relative registration error caused by platform vibration between the images of different band combinations is obtained. At the same time, dense ground control point matching is performed between each band of the enhanced multispectral image and the high-precision reference image DOM (digital orthophoto) of the corresponding area. The linear error between each band image and the reference image DOM is removed through the affine transformation model to obtain the nonlinear absolute registration error caused by platform vibration;

[0010] S4, H-infinity filter fusion: select a combination, fuse the relative registration error obtained by dense matching with the absolute registration error obtained by dense ground control point matching through H-infinity filtering to obtain a fusion result;

[0011] S5, bicubic interpolation image compensation: the fusion result is used to compensate the image using the bicubic interpolation method.

[0012] Furthermore, the image enhancement processing in S2 performs wallis filter image enhancement processing on the multispectral images of each band respectively, enhances the pattern and features of the image, and reduces the linear grayscale difference between the two bands. The image enhancement processing calculation formula is:

[0013]

[0014] Among them, g(x,y) and f(x,y) are the original image and the filtered image respectively, r 1 and r 0 is the multiplicative coefficient and the additive coefficient, m f and f are the expected mean and expected standard deviation of the local window, respectively, g and gare the mean and standard deviation of the local window respectively, c is the expansion constant of the image standard deviation, and b is the brightness coefficient of the image.

[0015] Furthermore, the dense matching and registration error acquisition of the multispectral image in S3 includes:

[0016] S31, setting configuration parameters: setting the window size of the normalized correlation coefficient method NCC matching, setting the window of the least square method, and setting the correlation coefficient threshold;

[0017] S32, coarse matching: performing normalized correlation coefficient (NCC) matching between each band of the enhanced multispectral image to obtain the correlation coefficient;

[0018] S33, fine matching: Based on the result of coarse matching, the least square method is used to perform fine matching at the sub-pixel level, and equations are established for the points in the left and right matching windows. The affine transformation model is considered to optimize the affine deformation coefficient and the linear coefficient, and the fine relative registration error is obtained by calculating the fine coordinate difference of the same-name points;

[0019] S34, dense ground control point (GCP) matching: downsample the reference image DOM to the same resolution as the enhanced multispectral image bands, and perform dense matching with the reference image DOM for each band image, and select the correlation coefficient and the size of the matching window to ensure that sufficient control points (GCPs) are found in the reference image DOM and the corresponding same-name points are found in the band image;

[0020] S35, error calculation and correction: project the control points into the image space of the corresponding band image, calculate the error by comparing the coordinates of the control points in the band image with the coordinates of the corresponding points with the same name, refine the error using the affine transformation model, remove the translation, and obtain the difference value that changes with the line. Fit the difference value with a sine function to obtain the absolute registration error curve between each band image and the reference image DOM data.

[0021] Furthermore, the correlation coefficient calculation formula in S32 is:

[0022]

[0023] Among them, S x,y Indicates a local window of a certain size with x, y as the upper left corner coordinates in the matching band, and g x+r,y+c It represents a local window of a certain size with x+r, y+c as the upper left corner coordinate in the reference band. and They represent the grayscale mean of the corresponding band window, r and c represent the row and column displacements corresponding to the search area, and ρ(x,y) represents the correlation coefficient.

[0024] Furthermore, the equations established for the points in the left and right matching windows in S33 are expressed as:

[0025] g 1 (x 1 ,y 1 )+n 1 (x 1 ,y 1 )=h 0 +h 1 *g 2 (x 2 ,y 2 )+n 2 (x 2 ,y 2 );

[0026] Among them, x 1 ,y 1 ,x 2 ,y 2 are the matching window pixel coordinates of the matching band image and the reference band image, respectively, 1 (x 1 ,y 1 ) and g 2 (x 1 ,y 1 ) represent their grayscale, n(x 1 ,y 1 )、n(x 2 ,y 2 ) is the gray value noise, h 0 and h 1 represents the linear deformation coefficient;

[0027] The affine transformation model is expressed as:

[0028]

[0029] Among them, a 0 ,a 1 ,a 2 ,b 0 ,b 1 ,b 2 is the affine deformation coefficient.

[0030] Furthermore, the H-infinity filter fusion in S4 includes:

[0031] S41, select fusion object: determine the fusion object, select the combination with the shortest interval between bands for relative registration error, so that the result is more precise and the error is smaller, and select the data group with the highest control point quality and the best fitting effect for absolute registration error;

[0032] S42, determining the spatial state equation of fusion filtering: establishing a spatial state equation of fusion input values ​​representing absolute registration error and relative registration error;

[0033] S43, defining a fusion filtering system and a cost function: defining a fusion filtering system, and establishing a cost function for implementing fusion filtering;

[0034] S44, executing a fusion filtering algorithm: applying a fusion filtering algorithm according to a set cost function to ensure that the fusion filtering at each moment has a solution;

[0035] S45, result fitting: performing sinusoidal function fitting on all fused output results to obtain the fused absolute registration error value.

[0036] Furthermore, the spatial state equation of the fusion filter is expressed as:

[0037]

[0038] Among them, the two state quantities x in the state equation 1 and x 2 represents the jitter value of two different bands at the same location, k-1 and k represent the current imaging time and the next imaging time respectively, y k (1) and y k (2) represents the fused input value.

[0039] Furthermore, the fusion filtering system is expressed as:

[0040] X k =F k X k-1 +W k

[0041] y k =H k X k +V k

[0042] z k =L k X k ;

[0043] in is the state transfer matrix, is the measurement matrix, L k =

[10] , z k is the expected estimated variable, W k and V k is a Gaussian white noise term, [x 1,k ,x 2,k ] T =X k ,[w1,k ,w 2,k ] T =W k and [v 1,k ,v 2,k ] T =V k ;

[0044] The cost function for implementing fusion filtering is expressed as:

[0045]

[0046] Where P 0 , Q k and R k are the covariance of the initial state, system process noise and observation noise, respectively, S k represents a symmetric positive definite matrix.

[0047] Furthermore, the calculation formula of the fusion filtering algorithm is:

[0048]

[0049] in, is the fusion output, the parameter θ is the boundary of the specified cost function, and J 1 <1 / θ is given.

[0050] In order to ensure that the fusion filter has a solution at every moment, the following conditions must be met:

[0051]

[0052] Furthermore, the bicubic interpolation image compensation in S5 includes:

[0053] S51, calculating the tremor value: for each row of the tremor image, substituting the imaging time into the fitting function of the absolute registration error, and calculating the corresponding tremor value;

[0054] S52, determining the position of the interpolation point: for the pixel position to be interpolated, determining the data points around it, and selecting multiple neighboring points around it as references for interpolation;

[0055] S53, calculate weights: for multiple neighbor points in the interpolation area, calculate weights according to the distances between the neighbor points and the interpolation point. The weights are determined by the interpolation basis function, which determines the influence of each neighbor point on the interpolation point. The interpolation basis function is expressed as:

[0056]

[0057] S54, calculating an estimated value of the interpolation point: using the calculated weight and the values ​​of the neighboring points, calculating the estimated value of the interpolation point by a weighted average method;

[0058] S55, row level operation: for each row of pixels in the image, repeat steps S51 to S54 to generate a compensated image;

[0059] S56, effect evaluation: performing comparative analysis of the relative registration errors before and after compensation and the relative registration errors of the compensation results before and after fusion to evaluate the effect of fusion.

[0060] Beneficial effects of the present invention:

[0061] The present invention, through a satellite image tremor error detection and correction method based on infinite filter fusion, can effectively fuse multispectral tremor detection data from several sources, and make full use of their respective advantages and avoid their disadvantages, so as to achieve more accurate calculation and estimation of the absolute jitter error of the satellite platform tremor. The method is simple, efficient, easy to implement, with small estimation error and strong robustness. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0063] Figure 1 A schematic diagram of a detection and correction method flow chart of an embodiment of the present invention;

[0064] Figure 2 A schematic diagram of a relative registration error curve obtained by dense matching according to an embodiment of the present invention;

[0065] Figure 3 A schematic diagram of an absolute registration error curve obtained by dense ground control point matching according to an embodiment of the present invention;

[0066] Figure 4 is a schematic diagram of a fusion curve according to an embodiment of the present invention;

[0067] Figure 5 A schematic diagram of comparative analysis of relative registration errors before and after compensation according to an embodiment of the present invention;

[0068] Figure 6 It is a schematic diagram of time domain comparative analysis of compensation results before and after fusion according to an embodiment of the present invention.

[0069] Figure 7 It is a schematic diagram of frequency domain comparison and analysis of compensation results before and after fusion according to an embodiment of the present invention;

[0070] Figure 8is a schematic diagram of relative registration error before compensation according to an embodiment of the present invention;

[0071] Fig. 9 It is a schematic diagram of relative registration error compensated before fusion according to an embodiment of the present invention;

[0072] Fig.10 It is a schematic diagram of relative registration error compensated after fusion according to an embodiment of the present invention. DETAILED DESCRIPTION

[0073] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with specific embodiments.

[0074] It should be noted that, unless otherwise defined, the technical terms or scientific terms used in the present invention should be understood by people with ordinary skills in the field to which the present invention belongs. The "first", "second" and similar words used in the present invention do not indicate any order, quantity or importance, but are only used to distinguish different components. "Include" or "comprise" and similar words mean that the elements or objects appearing before the word include the elements or objects listed after the word and their equivalents, without excluding other elements or objects. "Connect" or "connected" and similar words are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. "Up", "down", "left", "right" and the like are only used to indicate relative positional relationships. When the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0075] like Figure 1-Figure 10 As shown, the satellite image tremor error detection and correction method based on infinite filter fusion includes the following steps:

[0076] S1. Separate the multispectral image into spectral bands to obtain the visible light bands of blue, green, and red, which are marked as B1, B2, B3, and the infrared band B4 respectively.

[0077] S2. Perform wallis filter image enhancement processing on each band to obtain enhanced multispectral images. The detailed steps are as follows:

[0078] The B1, B2, and B3 bands are each subjected to a Wallis filter image enhancement process. The main purpose is to enhance the pattern and features of the image and reduce the linear grayscale difference between the two bands so as to have better reliability in subsequent image matching. Use the following formula:

[0079]

[0080] Among them, g(x,y) and f(x,y) are the original image and the filtered image respectively, r1 and r 0 is the multiplicative coefficient and the additive coefficient, m f and f are the expected mean and expected standard deviation of the local window, respectively, g and g are the mean and standard deviation of the local window, c is the expansion constant of the image standard deviation, and b is the brightness coefficient of the image. Input the original image g(x,y), m f 、s f , c and b, we can get the filtered image f(x,y). f 、s f According to the grayscale difference setting between the two bands, the value of c is generally 0.75-1, and the value of b can be set to 0.5-1.

[0081] S3. Select different band combinations for dense matching. There are three combinations: B2-B1, B3-B1 and B3-B2. The row offsets are 76, 140 and 64 respectively. Then, perform dense matching respectively. The detailed steps are as follows:

[0082] First, the window size of the normalized cross-correlation coefficient method NCC matching is set to 11, the window size of the least squares method is set to 9, and the correlation coefficient can be set to 0.85.

[0083] In order to improve the matching accuracy, a two-step sub-pixel matching method is introduced: one is a “coarse matching” step and the other is a subsequent “fine matching” step.

[0084] In the first step, the two bands are matched using the normalized correlation coefficient method (NCC), mainly using the following formula:

[0085]

[0086] Among them, S x,y Indicates a local window of a certain size with x, y as the upper left corner coordinates in the matching band, and g x+r,y+c It represents a local window of a certain size with x+r, y+c as the upper left corner coordinates in the reference band. and They represent the grayscale mean of the corresponding band window respectively, where r and c represent the row and column displacements corresponding to the search area. Usually, the search area selects eight neighborhoods, and ρ(x,y) represents the correlation coefficient, which ranges from [-1,1].

[0087] Next, the "fine matching" step is performed to obtain accurate disparity at the sub-pixel level based on the results obtained from the rough matching step, mainly using the least squares algorithm. Equations are established for the points in the left and right matching windows so that the points in the windows satisfy the equation as much as possible:

[0088] g 1 (x 1 ,y 1 )+n 1 (x 1 ,y 1 )=h 0 +h 1 *g 2 (x 2 ,y 2 )+n 2 (x 2 ,y 2 ) (2)

[0089] Among them, x 1 ,y 1 ,x 2 ,y 2 are the matching window pixel coordinates of the matching band image and the reference band image, respectively, 1 (x 1 ,y 1 ) and g 2 (x 1 ,y 1 ) represent their grayscale, n(x 1 ,y 1 )、n(x 2 ,y 2 ) is the gray value noise, h 0 and h 1 represents the linear deformation coefficient;

[0090] Considering the affine transformation model, then:

[0091]

[0092] Among them, a 0 a 1 a 2 b 0 b 1 b 2 is the affine deformation coefficient. After optimally solving the affine deformation coefficient and the linear coefficient under the condition of satisfying equation (2), the precise sub-pixel coordinates can be obtained. The precise relative registration error can be obtained by calculating the precise coordinate difference of two points of the same name. The final result is as follows: Figure 2 shown.

[0093] At the same time, the B1, B2, and B3 bands are matched with the high-precision reference image DOM with a spatial resolution of 1m for dense ground control points. The detailed steps are as follows:

[0094] First, the high-precision reference image DOM needs to be downsampled to the same resolution as the images of each band. Next, for each band image, it is densely matched with the reference image DOM, and the appropriate correlation coefficient and the size of the matching window are selected to ensure that enough points are found in the reference image DOM, and the corresponding points of the same name are found in the band image. These points found in the DOM are selected as control points (GCP). Subsequently, a large number of control points are projected into the image space of the corresponding band image using the rational function model (RFM). In this process, the corresponding rational polynomial coefficients (RPC) need to be selected for different bands. Through this process, the accurate coordinates of the control points in each band image can be obtained.

[0095] Next, the coordinates of the control points in the band image and their corresponding points of the same name are compared, the differences are calculated, and they are refined through affine modeling to eliminate linear errors. At the same time, the system translation is removed to obtain the difference values ​​that change with the line. Finally, these difference values ​​are fitted with a sine function to obtain the registration error (i.e., absolute registration error) curve between the images of each band and the DOM data. The results are as follows: Figure 3 shown.

[0096] S4. Select a suitable combination and fuse the relative registration error and the absolute registration error through H-infinity filtering to obtain the fusion result. The detailed steps are as follows:

[0097] First, determine the best fusion object, relative registration error and absolute registration error. The relative registration error selects the combination with the shortest interval between bands, which can make the result more precise and the error smaller. Here, the B2-B3 band can be selected, and the absolute registration error selects the group with the highest control point quality and the best fitting effect. Next, determine the spatial state equation of the fusion filter:

[0098]

[0099] Among them, the two state quantities x in the state equation 1 and x 2 represents the jitter value of two different bands at the same location, k-1 and k represent the current imaging time and the next imaging time respectively, y k (1) and y k (2) represents the fused input value.

[0100] By expressing [x 1,k ,x 2,k ] T =X k ,[w 1,k ,w 2,k ] T =W k and [v1,k ,v 2,k ] T =V k , the corresponding fusion filter system is:

[0101] X k =F k X k-1 +W k

[0102] y k =H k X k +V k

[0103] z k =L k X k ;

[0104] in is the state transfer matrix, is the measurement matrix, L k =

[10] , z k is the expected estimated variable, W k and V k is a Gaussian white noise term.

[0105] The cost function for implementing fusion filtering is:

[0106]

[0107] Where P 0 , Q k and R k are the covariance of the initial state, system process noise and observation noise, respectively, S k represents a symmetric positive definite matrix.

[0108] The fusion filtering algorithm is derived from the following formula:

[0109]

[0110] in, is the fusion output, the parameter θ is the boundary of the specified cost function, and J 1 <1 / θ is given. In order to ensure that the fusion filter has a solution at every moment, the following conditions need to be met:

[0111]

[0112] Finally, all the output results are merged Perform sine function fitting. The absolute registration error value after fusion can be obtained. Figure 4 shown.

[0113] S5, using the bicubic interpolation method to perform image tremor compensation on the fusion result obtained in S4, the detailed steps are as follows:

[0114] For the nth row of the tremor image, assume that its imaging time is t n , we can use t n Substituting into the fitting function of the absolute registration error, the tremor value Δd of the row is obtained. n Next, add the tremor value Δd to each pixel in the corresponding row n , for the position of the point that needs to be interpolated, that is, the position of a certain pixel, determine the data points around it and take the 16 neighboring points around it.

[0115] For the 16 neighboring points in the interpolation area, the distance between them and the interpolation point is calculated respectively, and the influence weight of the neighboring point on the interpolation point is determined according to the size of the distance. These weights can be expressed by the interpolation basis function:

[0116]

[0117] According to the interpolation weight and the value of the neighboring point, the estimated value of the interpolation point is calculated by weighted average or other methods. The above operation is performed on each row of pixels to obtain the compensated image. Figure 5-Figure 10 As shown, a comparative analysis of the relative registration errors before and after compensation and a comparative analysis of the relative registration errors of the compensation results before and after fusion were performed to evaluate the fusion effect.

[0118] Those skilled in the art should understand that the discussion of any of the above embodiments is merely illustrative and is not intended to imply that the scope of the present invention is limited to these examples. Under the concept of the present invention, the technical features in the above embodiments or different embodiments may be combined, the steps may be implemented in any order, and there are many other variations of the different aspects of the present invention as described above, which are not provided in detail for the sake of simplicity.

[0119] The present invention is intended to cover all such substitutions, modifications and variations that fall within the broad scope of the claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A satellite image tremor error detection and correction method based on infinite filter fusion, characterized in that: The following steps are involved: S1, band separation of multispectral images: band separation of the original multispectral images to obtain multispectral images of each band; S2, image enhancement processing: perform wallis filter image enhancement processing on the multispectral images of each band to obtain enhanced multispectral images; S3, dense matching and registration error acquisition of multispectral images: dense matching is performed between the enhanced multispectral image bands, and the nonlinear relative registration error caused by platform vibration between the images of different band combinations is obtained based on the normalized cross correlation coefficient method NCC matching and the least squares matching. At the same time, dense ground control point matching is performed between each band of the enhanced multispectral image and the high-precision reference image DOM of the corresponding area, and the linear error between each band image and the reference image DOM is removed by the affine transformation model to obtain the nonlinear absolute registration error caused by platform vibration. The dense matching and registration error acquisition of the multispectral image includes: S31, setting configuration parameters: setting the window size of the normalized correlation coefficient method NCC matching, setting the window of the least square method, and setting the correlation coefficient threshold; S32, coarse matching: performing normalized correlation coefficient (NCC) matching between each band of the enhanced multispectral image to obtain the correlation coefficient; S33, fine matching: Based on the result of coarse matching, the least square method is used to perform fine matching at the sub-pixel level, and equations are established for the points in the left and right matching windows. The affine transformation model is considered to optimize the affine deformation coefficient and the linear coefficient, and the fine relative registration error is obtained by calculating the fine coordinate difference of the same-name points; S34, dense ground control point matching: downsample the reference image DOM to the same resolution as the enhanced multispectral image in each band, perform dense matching with the reference image DOM for each band image, select the correlation coefficient and the size of the matching window, and find the corresponding points with the same name in the band image; S35, error calculation and correction: project the control point into the image space of the corresponding band image, calculate the error by comparing the coordinate difference of the control point in the band image with the coordinate difference of the corresponding point of the same name, refine the error by using the affine transformation model, remove the translation, obtain the difference value that changes with the line, and perform sinusoidal function fitting on the difference value to obtain the absolute registration error curve between each band image and the reference image data; S4, H-infinity filter fusion: Select a combination and fuse the relative registration error obtained by dense matching with the absolute registration error obtained by dense ground control point matching through H-infinity filtering to obtain the fusion result, including: S41, select fusion object: determine the fusion object, select the combination with the shortest interval between bands for relative registration error, and select the data group with the highest control point quality and the best fitting effect for absolute registration error; S42, determining the spatial state equation of fusion filtering: establishing a spatial state equation of fusion input values ​​representing absolute registration error and relative registration error; S43, defining a fusion filtering system and a cost function: defining a fusion filtering system, and establishing a cost function for implementing fusion filtering; S44, executing a fusion filtering algorithm: applying a fusion filtering algorithm according to a set cost function to ensure that the fusion filtering at each moment has a solution; S45, result fitting: performing sinusoidal function fitting on all fused output results to obtain the fused absolute registration error value; S5, bicubic interpolation image compensation: the fusion result is used to compensate the image using the bicubic interpolation method.

2. The satellite image tremor error detection and correction method based on infinite filter fusion according to claim 1 is characterized in that: The image enhancement processing in S2 is to perform wallis filter image enhancement processing on the multispectral images of each band respectively, enhance the pattern and features of the image, and reduce the linear grayscale difference between the two bands. The image enhancement processing calculation formula is: ； in, and are the original image and the filtered image respectively. and are the multiplicative and additive coefficients, and are the expected mean and expected standard deviation of the local window, respectively. and are the mean and standard deviation of the local window, c is the expansion constant of the image standard deviation, and b is the brightness coefficient of the image.

3. The satellite image tremor error detection and correction method based on infinite filter fusion according to claim 1 is characterized in that: The correlation coefficient calculation formula in S32 is: ； in, Indicates matching bands with x , y As a local window with a certain size of upper left corner coordinates, It represents a local window of a certain size with x+r, y+c as the upper left corner coordinate in the reference band. and Respectively represent the grayscale mean of the corresponding band window, r and c represents the row and column displacements corresponding to the search area, Represents the correlation coefficient.

4. The satellite image tremor error detection and correction method based on infinite filter fusion according to claim 3 is characterized in that: The equations established for the points in the left and right matching windows in S33 are expressed as: ； in, , are the matching window pixel coordinates of the matching band image and the reference band image respectively, and Represent their grayscale respectively, , is the gray value noise, and represents the linear deformation coefficient; The affine transformation model is expressed as: ； in, is the affine deformation coefficient.

5. The satellite image tremor error detection and correction method based on infinite filter fusion according to claim 1 is characterized in that: The spatial state equation of the fusion filter is expressed as: ； ； Among them, the two state quantities in the state equation are x 1 and x 2 means the jitter value of two different bands at the same position. k -1 and k Respectively represent the current imaging time and the next imaging time, and Represents the input value of the fusion.

6. The satellite image tremor error detection and correction method based on infinite filter fusion according to claim 5 is characterized in that: The fusion filtering system is expressed as: ； in is the state transfer matrix, is the measurement matrix, , is the expected estimated variable, and is a Gaussian white noise term, , and ; The cost function for implementing fusion filtering is expressed as: ； in , and are the covariances of the initial state, system process noise and observation noise, respectively. represents a symmetric positive definite matrix.

7. The satellite image tremor error detection and correction method based on infinite filter fusion according to claim 6 is characterized in that: The calculation formula of the fusion filtering algorithm is: ； in, is the fusion output, parameter is the boundary of the specified cost function, given by <1 / given; In order to ensure that the fusion filter has a solution at every moment, the following conditions must be met: 。 8. The satellite image tremor error detection and correction method based on infinite filter fusion according to claim 7 is characterized in that: The bicubic interpolation image compensation in S5 includes: S51, calculating the tremor value: for each row of the tremor image, substituting the imaging time into the fitting function of the absolute registration error, and calculating the corresponding tremor value; S52, determining the position of the interpolation point: for the pixel position to be interpolated, determining the data points around it, and selecting multiple neighboring points around it as references for interpolation; S53, calculating weights: for multiple neighbor points in the interpolation area, weights are calculated according to the distances between the neighbor points and the interpolation point. The weights are determined by the interpolation basis function, which determines the influence of each neighbor point on the interpolation point. The interpolation basis function is expressed as: ； S54, calculating an estimated value of the interpolation point: using the calculated weight and the values ​​of the neighboring points, calculating the estimated value of the interpolation point by a weighted average method; S55, row level operation: for each row of pixels in the image, repeat steps S51 to S54 to generate a compensated image; S56, effect evaluation: performing comparative analysis of the relative registration errors before and after compensation and the relative registration errors of the compensation results before and after fusion to evaluate the effect of fusion.

Images