A remote sensing image restoration method based on precise vibration detection and estimation

Abstract

The present invention relates to a remote sensing image restoration method based on precise detection and estimation of chatter. The method comprises: acquiring a multispectral remote sensing image and sensor parameters, using a dense matching algorithm to obtain parallax, directly separating chatter from the parallax in the vertical-to-track direction, first converting the parallax into a relative residual in the along-track direction, and then separating chatter from the relative residual; converting chatter information into a point spread function, estimating the point spread function for each row of the multispectral remote sensing image using the sensor's order and integration time; and restoring the image using an optimal windowed RL algorithm based on the point spread function. Compared with existing technologies, the present invention significantly improves the radiometric and geometric quality of the restored image, enhances local detail clarity, and corrects geometric offset.

Description

Technical Field

[0001] The invention relates to a remote sensing image restoration method based on precise vibration detection and estimation. Background Art

[0002] During flight, high-resolution remote sensing satellites are inevitably subject to the influence of external environmental factors (such as Earth's gravity and temperature) and internal factors (such as the operation of internal loads like attitude control systems), resulting in inevitable flutter. Flutter reduces camera attitude stability, causing temporal variations in imaging attitude and resulting in pointing angle errors, which in turn degrade the geometric and radiometric quality of remote sensing images. If flutter is not effectively eliminated, it can limit subsequent applications of remote sensing imagery, such as geometric positioning, surface change detection, image stitching, and digital elevation model (DEM) generation. With the rapid development of remote sensing satellite cameras, image spatial resolution is increasing, and the impact of attitude flutter on image quality is also increasing. Time-delayed integrated charge-coupled devices (TDICCDs) use multiple CCDs to scan the same object multiple times, utilizing time-delayed integration technology to generate images. This improves camera sensitivity and image signal-to-noise ratio, making it a mainstream imaging device. Therefore, eliminating the impact of attitude flutter on TDICCD remote sensing images is a highly significant task.

[0003] Satellite attitude flutter detection is fundamental to image restoration, and its detection accuracy determines the overall restoration effectiveness. There are two main types of flutter detection methods. One type directly acquires information about platform attitude changes based on high-performance attitude sensors (such as angular displacement sensors and angular accelerometers). The accuracy of this method depends on the frequency and accuracy of attitude measurements. This type of method has been successfully applied to the French Pleiades, ALOS, ZY-3, Yaogan-26, and Gaofen-9 satellites. The other type of method indirectly inverts satellite attitude flutter from remote sensing imagery, detecting flutter from the imagery and directly analyzing flutter information on the image plane. This eliminates the need for attitude-related transformations and is more conducive to flutter image restoration. This method has been applied to Terra, HiRISE, LRO, PLEIADES-HR, ZY-3, and GaoFen-1 02 satellites.

[0004] Although there have been some studies on the restoration of the vibration effects of remote sensing images, there are still some areas that need improvement: (1) The estimation of along-track vibration detection needs to consider the influence of terrain to improve the accuracy of along-track vibration detection estimation; (2) The remote sensing image restoration process needs to consider the influence of TDI. Summary of the Invention

[0005] The purpose of the present invention is to overcome the defects of the above-mentioned prior art and provide a remote sensing image restoration method based on precise vibration detection and estimation.

[0006] The purpose of the present invention can be achieved by the following technical solutions:

[0007] A remote sensing image restoration method based on precise chatter detection and estimation, comprising:

[0008] Flutter detection steps: Acquire multispectral remote sensing images and sensor parameters, use a dense matching algorithm to obtain parallax, and directly separate flutter from the parallax in the vertical direction. In the along-track direction, first convert the parallax into a relative residual, and then separate flutter from the relative residual.

[0009] Point spread function estimation step: converting the chatter information into a point spread function, and using the sensor's order and integration time to estimate the point spread function for each row of the multispectral remote sensing image;

[0010] Image restoration steps: restore the image using the optimal window RL algorithm based on the point spread function.

[0011] Furthermore, in the vibration detection step, the dense matching algorithm includes: obtaining integer pixel matching results of the multispectral remote sensing image through a normalized cross-correlation matching algorithm; and obtaining sub-pixel matching results using a PEF phase algorithm.

[0012] Furthermore, in the flutter detection step, the vertical track direction parallax is modeled, and the vertical track direction parallax is converted into vertical track flutter image shift using a first conversion formula;

[0013] The model of the vertical track parallax is:

[0014]

[0015] Among them, A u ,ω u , and b u are the amplitude, frequency, initial phase and trend term of the vertical track parallax respectively; the expression of the first conversion formula is:

[0016]

[0017]

[0018]

[0019] Among them, A f_Cross is the amplitude of vertical track flutter image shift, φ f_Cross is the initial phase of the vertical track flutter image shift, and Δt is the observation time difference between adjacent bands.

[0020] Furthermore, in the flutter detection step, the along-track parallax is modeled by relative residuals, and then the relative residuals are converted into along-track flutter image motion by a second conversion formula;

[0021] The model of the along-track parallax is:

[0022]

[0023]

[0024] Where s(t) represents the relative residual, g1(t) and g2(t) are the parallaxes along the track between the first and second bands and between the second and third bands, dp1 and dp2 represent the number of physical interval lines between the first and second bands and between the second and third bands, respectively. s ,ω s , and b s are the amplitude, frequency, initial phase and trend terms of the relative residual respectively;

[0025] The expression of the second conversion formula is:

[0026]

[0027]

[0028]

[0029] Among them, A f_Along is the amplitude of the along-track flutter image shift, φ f_Along is the initial phase of the along-track flutter image shift, Δt1 and Δt2 are the observation time differences between the first and second bands and between the second and third bands, respectively.

[0030] Furthermore, in the point spread function estimation step, the size of the point spread function is expressed as:

[0031]

[0032] Where r is the number of rows, c is the number of columns, a and b are the vibration amplitudes in the along-track and perpendicular-track directions respectively. Represents rounding up.

[0033] Furthermore, in the point spread function estimation step, the element values ​​in the point spread function represent the weights of the pixels at the corresponding positions. The weights of each position of the point spread function are updated by traversing each pixel value point, and the weights of four positions around the point spread function are determined using an inverse bilinear sampling method. During the traversal process, the weights are superimposed. After traversing each vibration value sampling point, the point spread function is normalized.

[0034] The coordinate calculation formula of the vibration value sampling point in the point spread function is as follows:

[0035]

[0036] Among them, x PSF and y PSF is the coordinate of the vibration value sampling point in the point spread function, Δx and Δy are the vibration values ​​in the along-track and perpendicular-track directions, respectively.

[0037] Furthermore, the optimal window RL algorithm is used to perform maximum likelihood estimation of the image based on the Poisson noise model. The iterative formula is as follows:

[0038]

[0039] Among them, g is the blurred image, f is the restored image, h is the point spread function, h T is the transpose of h.

[0040] Furthermore, before image restoration, a window function is applied to the image to be restored to suppress the ringing phenomenon at the image boundary.

[0041] Compared with the prior art, the present invention has the following beneficial effects:

[0042] The present invention proposes a time delay integration (TDI) remote sensing image restoration method based on precise vibration detection and estimation and an optimal window RL (Richardson-Lucy) algorithm. First, the vibration information is precisely detected based on multispectral images. Then, a point spread function estimation algorithm is used to convert the vibration information into a point spread function. Finally, the remote sensing image is restored using an optimal window RL algorithm based on the point spread function. The present invention considers the influence of terrain factors during along-track vibration detection to achieve precise along-track vibration estimation. Simultaneously, the optimal window RL algorithm is used in the image restoration step to effectively suppress ringing in edge images. This significantly improves the radiometric and geometric quality of the restored image, making local details clearer and correcting geometric offsets. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 It is a schematic diagram of the overall process of the present invention.

[0044] Figure 2a This is the vertical track flutter diagram of the experimental image.

[0045] Figure 2b This is the flutter map of the experimental image along the track.

[0046] Figure 3a It is the original multispectral image of the experiment.

[0047] Figure 3b This is the optimal window WNR restoration result image of the experimental image.

[0048] Figure 3c This is the experimental image and the restoration result image of this embodiment.

[0049] Figure 4a Restore the front disparity map for the experimental image.

[0050] Figure 4b It is the original and front disparity map of the experimental image. DETAILED DESCRIPTION

[0051] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0052] like Figure 1 As shown, this embodiment provides a remote sensing image restoration method based on precise chatter detection and estimation. It consists of three main steps: chatter detection, point spread function estimation, and image restoration. The first step uses a dense matching method to obtain parallax, and then directly separates chatter from the parallax in the vertical-track direction. Due to the influence of terrain in the along-track direction, the parallax is first converted to a relative residual, and then the chatter is separated from the relative residual. The second step converts the chatter information into a point spread function, and uses information such as the CCD stage and integration time to estimate the point spread function for each image line. The final step is to restore the image using an optimal window RL algorithm based on the point spread function. These three steps are described in detail below.

[0053] 1. Flutter detection steps:

[0054] Flutter detection involves two steps: dense matching and disparity (or relative residual) separation. The dense matching method employed in this embodiment utilizes a coarse-to-fine matching strategy. First, a normalized cross-correlation matching algorithm is used to obtain integer-pixel matching results. Then, a Peak Evaluation Formula (PEF) phase correlation algorithm is used to obtain sub-pixel matching results.

[0055] Parallax modeling in vertical direction:

[0056]

[0057] A u ,ω u , and b u Represent the amplitude, frequency, initial phase and trend of parallax respectively. The vertical track parallax is converted into vertical track flutter image motion. The conversion formula is as follows:

[0058]

[0059]

[0060] in, A f_Cross is the amplitude of vertical track flutter image shift, φ f_Cross is the initial phase of the vertical track flutter image shift, and Δt is the observation time difference between adjacent bands.

[0061] Because terrain factors can produce a certain amount of parallax in the along-track direction, the along-track parallax obtained by dense matching cannot be directly converted into dither. Instead, the parallax between B1-B2 (from the first band to the second band) and B2-B3 (from the second band to the third band) needs to be converted into relative residuals.

[0062]

[0063] Where s(t) represents the relative residual, g1(t) and g2(t) are the parallaxes along the track between B1 and B2, and between B2 and B3, respectively. dp1 and dp2 represent the number of physical rows between B1 and B2, and between B2 and B3, respectively.

[0064] The relative residual s(t) is modeled using a sine function plus a trend term.

[0065]

[0066] A s ,ω s , and b s are the amplitude, frequency, initial phase and trend terms of the relative residual respectively.

[0067] The conversion formula from the relative residual along the track to the flutter function is as follows:

[0068]

[0069]

[0070] in, A f_Along is the amplitude of the along-track flutter image shift, φ f_Along is the initial phase of the along-track flutter image shift, Δt1 and Δt2 are the observation time differences between B1-B2 and B2-B3, respectively.

[0071] 2. Point Spread Function Estimation Steps

[0072] The point spread function, also known as the blur kernel, is spatially varying, meaning that each row of the image has a different point spread function, depending on the effects of dithering during the integration time. Point spread function estimation involves two steps: determining the size and updating the weights.

[0073] 1) Determine the size

[0074] The chatter information within the integration time is extracted based on the row number, and the chatter sampling time point is determined by the integration series and other information. Then, the chatter sampling point value is analyzed to obtain its amplitude. Since chatter values ​​can be positive or negative, and the size of the point spread function is usually an odd number, the size of the point spread function is:

[0075]

[0076] Where r is the number of rows, c is the number of columns, a and b are the vibration amplitudes in the along-track and perpendicular-track directions respectively. After determining the size, each matrix element is initialized to 0.

[0077] 2) Update weights

[0078] The element value in the point spread function represents the weight of the pixel at the corresponding position. The weight of each position of the point spread function is updated by traversing each pixel value point. The coordinate calculation formula of the vibration value sampling point in the point spread function is as follows:

[0079]

[0080] Among them, (x PSF ,y PSF ) is the coordinate of the vibration value sampling point in the point spread function, Δx and Δy are the vibration values ​​in the along-track and vertical-track directions respectively. Then the inverse bilinear sampling method is used to determine the point spread function (x PSF ,y PSF ) The weights of the four positions around . During the traversal process, the weights are superimposed, and after traversing each vibration value sampling point, the point spread function is normalized.

[0081] 3. Image Restoration Steps

[0082] The Richardson Lucy algorithm, also known as the RL algorithm, performs maximum likelihood estimation on images based on the Poisson noise model. The iterative formula is as follows:

[0083]

[0084] Among them, g is the blurred image, f is the restored image, h is the point spread function, h T is the transpose of h.

[0085] Before image restoration, this embodiment applies a windowing function to the image to be restored to suppress ringing at the image boundaries, i.e., g'(x,y) = g(x,y)·ω(x,y), where g'(x,y) is the windowed image. The window function is expressed as Equation (11). Its size matches the size of the image to be restored, with the center region being 1 and the edge regions being related to the size of the point spread function, with their values ​​gradually approaching 0 from the inside out.

[0086]

[0087] Verification experiment

[0088] I. Data Introduction

[0089] The experimental data used multispectral remote sensing images from the Ziyuan-3 satellite, with a ground resolution of 5.8 meters. The images contain four bands, with adjacent bands separated by a certain distance, forming a parallax structure. This example uses this parallax relationship to detect the vibration of the multispectral camera and perform image restoration based on the detected vibration information.

[0090] II. Flutter Detection and Analysis

[0091] The dense matching method is used to obtain the average disparity of each row of pixels between adjacent bands. In the vertical track direction, the spectrum analysis and adjustment method are used to fit the disparity to obtain the various parameter values ​​in the disparity curve expression. Finally, the disparity parameter value is converted into the vibration parameter value using formula (2) (3). In the along-track direction, the disparity between B1-B2 and B2-B3 in the along-track direction is converted into relative residual, and the vertical track direction disparity fitting method is used to fit the relative residual to obtain the amplitude, frequency, initial phase and trend term of the relative residual fitting curve. The parameter value of the relative residual is converted into the vibration parameter value in the along-track direction using formula (6) (7). The vibration parameter values ​​in the vertical track direction and along-track direction are shown in Table 1, and the vibration curve is shown in Figure 2a and Figure 2b shown.

[0092] Table 1 Experimental image chatter detection results

[0093]

[0094]

[0095] III. Image Restoration and Discussion

[0096] The detected vibration result is converted into a point spread function, and then the image is restored using the restoration method proposed in this embodiment. The restoration result is as follows: Figure 3a to Figure 3cAs shown in the figure. Regarding geometric quality, the improvement in geometric accuracy is verified by comparing the disparity values ​​before and after image restoration. The disparity value is the difference between the corresponding dither functions of the two images. If the disparity value after restoration is smaller than the disparity value before restoration, it indicates that the impact of dither on the image has been reduced to a certain extent. Regarding radiometric quality, typical radiometric quality evaluation indicators (average gradient, contrast, and autocorrelation coefficient) are used to evaluate the images before and after restoration. Table 2 shows a comparison of the radiometric quality of images in different bands before and after restoration. The results show that the restoration results of the method in this embodiment are superior to those using the optimal window Wiener filter method.

[0097] Table 2 Comparison of radiation quality of experimental images before and after restoration of different bands

[0098]

[0099] like Figure 4a The disparity map of B1 and B2 before restoration is shown, showing an obvious sine or cosine law; Figure 4b The disparity map of the restored images B1 and B2 is shown in Figure 1. The fluctuation range is significantly smaller than before restoration and there is no obvious period. This shows that the method of this embodiment can significantly improve the geometric quality of the image.

[0100] In summary, this embodiment proposes a TDI remote sensing image restoration method based on precise vibration detection and estimation and an optimal window RL algorithm. This method first detects vibration in multispectral remote sensing images, then converts the vibration information into a point spread function (PSF), and finally restores the remote sensing image using an optimal window RL algorithm. Experiments using images from the Ziyuan-3 satellite validated the effectiveness of this method. The experimental results demonstrate significant improvements in both radiometric and geometric image quality. The restored image quality surpasses that achieved by the optimal window WNR algorithm, with clearer local details and corrected geometric offsets.

[0101] The above describes in detail the preferred embodiments of the present invention. It should be understood that those skilled in the art can make numerous modifications and variations based on the concepts of the present invention without inventive effort. Therefore, any technical solutions that can be derived by those skilled in the art through logical analysis, reasoning, or limited experimentation based on the concepts of the present invention and the prior art should be within the scope of protection defined by the claims.

Claims

1. A remote sensing image restoration method based on precise vibration detection and estimation, characterized in that: include: Flutter detection steps: Acquire multispectral remote sensing images and sensor parameters, use a dense matching algorithm to obtain parallax, and directly separate flutter from the parallax in the vertical direction. In the along-track direction, first convert the parallax into a relative residual, and then separate flutter from the relative residual. Point spread function estimation step: converting the chatter information into a point spread function, and using the sensor's order and integration time to estimate the point spread function for each row of the multispectral remote sensing image; Image restoration steps: restore the image using the optimal window Richardson Lucy algorithm based on the point spread function; In the vibration detection step, the dense matching algorithm includes: obtaining integer pixel matching results of the multispectral remote sensing image by using a normalized cross-correlation matching algorithm; obtaining sub-pixel matching results by using a PEF phase algorithm; In the point spread function estimation step, the size of the point spread function is expressed as: Where r is the number of rows, c is the number of columns, a and b are the vibration amplitudes in the along-track and perpendicular-track directions respectively. represents rounding up; In the point spread function estimation step, the element values ​​in the point spread function represent the weights of the pixels at the corresponding positions. The weights of each position of the point spread function are updated by traversing each pixel value point, and the weights of four positions around the point spread function are determined using an inverse bilinear sampling method. During the traversal process, the weights are superimposed. After traversing each vibration value sampling point, the point spread function is normalized. The coordinate calculation formula of the vibration value sampling point in the point spread function is as follows: Among them, x PSF and y PSF is the coordinate of the vibration value sampling point in the point spread function, Δx and Δy are the vibration values ​​in the along-track and perpendicular-track directions respectively; The optimal window Richardson Lucy algorithm performs maximum likelihood estimation on the image based on the Poisson noise model. The iterative formula is as follows: Among them, g is the blurred image, f is the restored image, h is the point spread function, h T is the transpose of h; Before image restoration, a window function is added to the image to be restored to suppress the ringing phenomenon at the image boundary.

2. The remote sensing image restoration method based on precise vibration detection and estimation according to claim 1, characterized in that: In the flutter detection step, the vertical track direction parallax is modeled and the vertical track direction parallax is converted into vertical track flutter image motion by using a first conversion formula; The model of the vertical track parallax is: Among them, A u ,ω u , and b u are the amplitude, frequency, initial phase and trend term of the vertical track parallax respectively; The expression of the first conversion formula is: Among them, A f_Cross is the amplitude of vertical track flutter image shift, φ f_Cross is the initial phase of the vertical track flutter image shift, and Δt is the observation time difference between adjacent bands.

3. The remote sensing image restoration method based on precise vibration detection and estimation according to claim 1, characterized in that: In the flutter detection step, the along-track parallax is modeled using relative residuals, and then the relative residuals are converted into along-track flutter image motion using a second conversion formula; The model of the along-track parallax is: Where s(t) represents the relative residual, g1(t) and g2(t) are the parallaxes along the track between the first and second bands and between the second and third bands, dp1 and dp2 represent the number of physical interval lines between the first and second bands and between the second and third bands, respectively. s ,ω s , and b s are the amplitude, frequency, initial phase and trend terms of the relative residual respectively; The expression of the second conversion formula is: Among them, A f_Along is the amplitude of the along-track flutter image shift, φ f_Along is the initial phase of the along-track flutter image shift, Δt1 and Δt2 are the observation time differences between the first and second bands and between the second and third bands, respectively.

Images