Satellite jitter inversion estimation and compensation method considering time delay integration effect

Abstract

The present application relates to a kind of satellite flutter inversion estimation and compensation method considering time delay integration effect, comprising the following steps: according to TDI effect, the calculation of platform flutter is converted into integral form, and flutter inversion estimation model is constructed;According to known imaging time interval, multispectral image matching is carried out, and imaging parallax is obtained, and according to the flutter inversion estimation model, the flutter information is obtained;According to the flutter information, the attitude flutter rotation matrix is constructed, so that the real attitude rotation matrix containing attitude flutter is constructed, and the attitude information is obtained;According to the attitude information, the corrected image is obtained by virtual re-imaging.Compared with the prior art, the present application has good flutter detection and compensation performance.

Description

Technical Field

[0001] This invention relates to the field of on-orbit satellite attitude flutter detection technology, and in particular to a satellite flutter inversion estimation and compensation method that takes into account the time delay integral effect. Background Technology

[0002] High-resolution optical remote sensing satellites have received high attention from countries around the world, driving the continuous development of high-resolution remote sensing technology. However, attitude jitter limits the geometric quality and imaging performance of high-resolution optical satellites. If jitter is not compensated for in a timely and effective manner, it will cause a decline in image quality, which in turn will affect image registration and change detection, image stitching, and elevation modeling.

[0003] Currently, many scholars research methods for detecting attitude jitter in orbiting satellites, often categorizing them into two types based on the data source: indirect detection methods based on non-attitude sensors and direct detection methods based on attitude sensors. Indirect methods rely on satellite imagery or mapping products, which are easy to implement; however, image quality and ground terrain information can affect the accuracy of subsequent matching and feature extraction, making it a passive jitter detection method. Direct detection methods rely on high-precision attitude sensors, such as star sensors, gyroscopes, or angular displacement trackers. However, due to hardware limitations, domestically produced attitude sensors suffer from low sampling frequency, poor attitude determination accuracy, and insufficient reliability. Furthermore, for satellites already in orbit, it is impossible to detect jitter using additional hardware sensors. Therefore, it is necessary to detect and compensate for satellite attitude jitter using existing attitude measurement equipment without increasing hardware costs, thereby improving the accuracy of attitude data.

[0004] Currently, time-delayed integrated charge-coupled devices (TDICCDs) have replaced conventional CCDs as the mainstream sensor in high-resolution optical sensor imaging systems. In the multi-stage integration imaging process of TDICCDs, the jitter offset caused by each stage of jitter changes over time; therefore, the jitter impact on the synthesized image is no longer the same as that on a single-stage integrated image. If the multi-stage integration process is ignored, the error in detecting and estimating jitter information will increase with the increase of the integration stage number and frequency. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the existing technology and provide a satellite flutter inversion estimation and compensation method that considers the time delay integral effect, which has good flutter detection and compensation performance.

[0006] The objective of this invention can be achieved through the following technical solutions:

[0007] A satellite flutter inversion estimation and compensation method considering the time delay integral effect includes the following steps:

[0008] S1: Based on the TDI effect, the calculation of platform flutter is transformed into an integral form, and a flutter inversion estimation model is constructed;

[0009] S2: Perform multispectral image matching based on the known imaging time interval to obtain imaging parallax, and obtain flutter information based on the flutter inversion estimation model;

[0010] S3: Construct an attitude flutter rotation matrix based on the flutter information, thereby constructing a true attitude rotation matrix containing attitude flutter and obtaining attitude information;

[0011] S4: Obtain the corrected image through virtual re-imaging based on the attitude information.

[0012] Furthermore, the calculation of platform flutter is transformed into an integral form expression as follows:

[0013]

[0014] In the formula, d TDI (t) represents the image shift of the integrated image at time t, N1 is the integration series of the image, T is the single integration time, and A,f, These are the flutter amplitude, frequency, and phase, t i It is the initial integration time for line i.

[0015] Furthermore, the expression for the flutter inversion estimation model is:

[0016]

[0017] In the formula, s(t) represents the imaging parallax, Δt is the imaging time interval, and N1 and N2 are the integral series of the reference image and the image to be matched.

[0018] Furthermore, in step S2, the process of performing multispectral image matching based on the known imaging time interval is specifically as follows:

[0019] S201: Enhance the image using Wallis filtering to improve image contrast, and perform initial registration by translating the image according to the known imaging time interval;

[0020] S202: Construct matching windows in the image using the SVD-RANSC subpixel phase correlation method and obtain the subpixel offset between each matching window;

[0021] S203: Remove erroneous matches from the matching window.

[0022] Furthermore, in step S203, removing erroneous matches from the matching window specifically includes:

[0023] Calculate the normalized cross-correlation coefficient between matching windows and remove matching points with correlation coefficients less than a preset correlation threshold;

[0024] The disparity of matching point pairs in the same image row is statistically analyzed, and matching points with large offsets are removed according to the three-times standard error principle.

[0025] Furthermore, the relevant threshold value is within the range of 0.6-0.8.

[0026] Furthermore, in step S3, the process of constructing the attitude flutter rotation matrix is ​​specifically as follows:

[0027] Based on the flutter information, the image flutter along the rail and perpendicular to the rail is converted into attitude angle changes using sensor parameters. The attitude flutter rotation matrix is ​​then constructed based on these attitude angle changes, thereby creating a true attitude rotation matrix that includes attitude flutter. The expression for this true attitude rotation matrix is ​​as follows:

[0028] R = R ss R si

[0029] In the formula, R is the true attitude rotation matrix including attitude flutter, R ss R is the attitude flutter rotation matrix. si This indicates the attitude obtained by traditional satellite sensors.

[0030] Furthermore, the expression for calculating the attitude angle change is:

[0031]

[0032]

[0033] In the formula, f is the focal length, p is the pixel size, β is the off-axis angle, and Δω and These represent the changes in roll and pitch angles caused by flutter, j across (t) represents the image flutter along the track direction at time t, j along (t) represents the image flutter in the vertical direction at time t.

[0034] Furthermore, the attitude flutter rotation matrix R ss Given a positive definite matrix, a Taylor expansion approximates it as follows when the angle reaches the arcsecond level:

[0035]

[0036] In the formula, Let I represent the three attitude angles at a given moment, where e is the natural logarithm and I3 is the identity matrix. express Cross product.

[0037] Further, step S4 specifically includes:

[0038] S401: Construct an imaging model using smoothed attitude data, and calculate the ground coordinates (B,L,H) of the corrected image point (r',c') on the average elevation surface H through orthographic projection;

[0039] S402: Construct an imaging model using attitude information containing flutter information, and obtain the image point (r,c) on the original image by back-projecting the ground point coordinates (B,L,H);

[0040] S403: The coordinates (r', c') on the corrected image and the point (r, c) on the original image correspond to the same point. The grayscale information of the image point is obtained through bilinear interpolation.

[0041] S404: Repeat steps S401-S403 to obtain the grayscale value of each virtual pixel.

[0042] Compared with the prior art, the present invention has the following advantages:

[0043] This invention, based on the imaging principle of a TDICCD camera, constructs a precise inversion estimation model for satellite flutter that considers the time delay integral effect. It then performs dense matching on multispectral images to obtain flutter information, updates the attitude based on the estimated flutter information, and finally combines the updated precise attitude with virtual re-imaging to generate a corrected image, thus eliminating the influence of satellite flutter. Experiments using image and attitude data from the ZY-3 satellite demonstrate that this method has excellent flutter detection and compensation performance. Attached Figure Description

[0044] Figure 1 This is a flowchart illustrating a satellite flutter inversion estimation and compensation method considering time delay integration effect provided in an embodiment of the present invention.

[0045] Figure 2 This is a schematic diagram of an experimental image ((a) CCD1, (b) CCD2, (c) CCD3) provided in an embodiment of the present invention;

[0046] Figure 3 This is a schematic diagram of the attitude comparison result of a CCD1 provided in an embodiment of the present invention;

[0047] Figure 4 This is a schematic diagram of the attitude comparison result of a CCD2 provided in an embodiment of the present invention;

[0048] Figure 5 This is a schematic diagram of the attitude comparison result of a CCD3 provided in an embodiment of the present invention;

[0049] Figure 6This is a schematic diagram of the CCD1 parallax of B1 and B2 before and after flutter compensation in an embodiment of the present invention ((a) parallax along the rail direction and (b) parallax perpendicular to the rail direction).

[0050] Figure 7 This is a schematic diagram of the CCD2 parallax of B1 and B2 before and after flutter compensation in an embodiment of the present invention ((a) parallax along the rail direction and (b) parallax perpendicular to the rail direction).

[0051] Figure 8 This invention provides a CCD3 parallax of B1 and B2 before and after flutter compensation ((a) parallax along the track direction and (b) parallax perpendicular to the track direction). Detailed Implementation

[0052] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0053] Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0054] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0055] In the description of this invention, it should be noted that the terms "center," "upper," "lower," "left," "right," "vertical," "horizontal," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship in which the product of this invention is usually placed during use. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0056] It should be noted that the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.

[0057] Furthermore, terms such as "horizontal" and "vertical" do not imply that components must be absolutely horizontal or suspended, but rather that they can be slightly tilted. For example, "horizontal" simply means that its direction is more horizontal than "vertical," and does not mean that the structure must be completely horizontal, but can be slightly tilted.

[0058] Example 1

[0059] like Figure 1 As shown, this embodiment provides a satellite flutter inversion estimation and compensation method considering the time delay integral effect, including the following steps:

[0060] S1: Based on the TDI effect, the calculation of platform flutter is transformed into an integral form, and a flutter inversion estimation model is constructed;

[0061] S2: Perform multispectral image matching based on the known imaging time interval to obtain imaging parallax, and obtain flutter information based on the flutter inversion estimation model;

[0062] S3: Construct an attitude flutter rotation matrix based on the flutter information, thereby constructing a true attitude rotation matrix containing attitude flutter and obtaining attitude information;

[0063] S4: Obtain the corrected image through virtual re-imaging based on the attitude information.

[0064] The following is a detailed description of each step.

[0065] S1: Based on the TDI effect, the calculation of platform flutter is transformed into an integral form, and a flutter inversion estimation model is constructed;

[0066] In the multi-stage integration process of TDICCD, the image shift caused by platform flutter changes with time. Assuming that the platform flutter is a simple sine wave, the flutter after TDI integration can be expressed as:

[0067]

[0068] For a linear array TDICCD pushbroom sensor with parallax observation, theoretically, when the satellite's attitude is stable and there is no flutter, the parallax value of each row of images in the parallax image pair should be a constant, determined by the physical distance of the sensor array, orbital altitude, and flight speed.

[0069] When the satellite's attitude is affected by flutter, there is a geometric displacement caused by flutter in each image row. The disparity value of each image row is no longer constant, and there is a disparity difference value related to the geometric displacement value.

[0070]

[0071] In the formula, Δt is the imaging time interval, A,f, It is the flutter amplitude, frequency, and phase, t i The initial integration time is denoted as i, N1 and N2 are the integration series of the reference image and the image to be matched, T is the single integration time, and d is the integration time. TDI It represents the integral image shift, and s(t) represents the imaging parallax.

[0072] S2: Perform multispectral image matching based on the known imaging time interval to obtain imaging parallax, and obtain flutter information based on the flutter inversion estimation model;

[0073] To obtain the flutter offset, a pixel-by-pixel dense matching method based on singular value decomposition and random sample consensus (SVD-RANSC) subpixel phase correlation was adopted.

[0074] S201: First, Wallis filtering is used to enhance the image and improve its contrast. Initial registration is achieved by translating the image based on the known imaging time interval.

[0075] S202: Then, the sub-pixel offset between each window is obtained using the SVD-RANSAC sub-pixel phase correlation method. SVD-RANSAC has higher reliability and stronger robustness than Hoge's method and the RANSC algorithm.

[0076] S203: Simultaneously, possible mismatched points are eliminated, including the following two steps: calculate the normalized cross-correlation coefficient between matching windows and remove matching points with a correlation coefficient less than the correlation threshold (0.7 in the experiment); and statistically analyze the disparity of matching point pairs in the same image row and eliminate matching points with large offsets according to the three-times standard error principle.

[0077] S3: Construct an attitude flutter rotation matrix based on the flutter information, thereby constructing a true attitude rotation matrix containing attitude flutter and obtaining attitude information;

[0078] Generally, for images captured at or near the nadir, vertical flutter is primarily affected by the roll angle, while along-track flutter is primarily affected by the pitch angle. Therefore, sensor parameters are used to convert along-track and vertical flutter into attitude angle changes.

[0079]

[0080] In the formula, f is the focal length, p is the pixel size, and β is the off-axis angle. These represent the roll and pitch angle changes caused by flutter, respectively. Based on the attitude and orbit information, an attitude flutter rotation matrix R is then constructed. ss , which represents the rotation matrix of the satellite body in the orbital system.

[0081] Therefore, by combining the flutter information estimated from the image, the true attitude rotation matrix R, which includes attitude flutter, can be obtained:

[0082] R = R ss R si (4)

[0083] Rotation matrix R si This indicates the attitude obtained by traditional satellite sensors, without recording high-frequency jitter information.

[0084] R ss It is a positive definite matrix, and through Taylor expansion, when the angle is small (on the order of arcseconds), it can usually be approximated as:

[0085]

[0086] Let I represent the three attitude angles at a given moment, where e is the natural logarithm and I3 is the identity matrix. express Cross product, therefore formula (4) can be expressed as matrix multiplication.

[0087]

[0088] S4: Obtain the corrected image through virtual re-imaging based on the attitude information.

[0089] After the pose is updated, we obtain the corrected image based on virtual re-imaging. The specific steps are as follows.

[0090] S401: First, an imaging model is constructed using smoothed attitude data. Then, the ground coordinates (B,L,H) of the corrected image point (r',c') on the average elevation surface H are calculated by orthographic projection.

[0091] S402: Then, an imaging model is constructed using the attitude containing flutter information, and the ground point coordinates (B,L,H) are back-projected to obtain the image point (r,c) on the original image;

[0092] S403: The coordinates (r', c') on the corrected image and the point (r, c) on the original image correspond to the same point, but the image point coordinates may not fall on the center of the pixel. Therefore, the grayscale information of the image point is obtained by bilinear interpolation.

[0093] S404: Repeat steps S401-S403 to obtain the grayscale value of each virtual pixel.

[0094] Experiments and Analysis

[0095] The image data used in this invention is ZY-3 multispectral image, which includes three CCDs, and each CCD image consists of four bands.

[0096] The updated pose data was obtained using the method described above. To verify the correctness of the algorithm, it was compared with the original pose of ZY-3, and the results are as follows. Figure 2-5 As shown, we can see that the results are consistent with the original pose data of different CCDs.

[0097] After generating the corrected image, flutter detection is performed again to check whether the flutter has been effectively compensated. Figure 6 , Figure 7 and Figure 8 The disparity before and after flutter compensation is shown. We can clearly see that the period disappears and the disparity concentrates near 0. After effectively correcting jitter distortion, the root mean square error in the vertical direction decreases from 0.3 pixels to 0.1 pixels, and the root mean square error along the rail direction decreases from 0.2 pixels to 0.1 pixels. These results demonstrate the effectiveness of our method.

[0098] in conclusion

[0099] This invention, based on the imaging principle of a TDICCD camera, constructs a precise inversion estimation model for satellite flutter that considers the time delay integral effect. It then performs dense matching on multispectral images to obtain flutter information, updates the attitude based on the estimated flutter information, and finally combines the updated precise attitude with virtual re-imaging to generate a corrected image, thus eliminating the influence of satellite flutter. Experiments using image and attitude data from the ZY-3 satellite demonstrate that this method has excellent flutter detection and compensation performance.

[0100] The preferred embodiments of the present invention have been described in detail above. It should be understood that those skilled in the art can make numerous modifications and variations based on the concept of the present invention without creative effort. Therefore, all technical solutions that can be obtained by those skilled in the art based on the concept of the present invention through logical analysis, reasoning, or limited experimentation on the basis of existing technology should be within the scope of protection defined by the claims.

Claims

1. A satellite flutter inversion estimation and compensation method considering the time delay integral effect, characterized in that, Includes the following steps: S1: Based on the TDI effect, the calculation of platform flutter is transformed into an integral form, and a flutter inversion estimation model is constructed; The expression for the flutter inversion estimation model is as follows: In the formula, Indicates image parallax. It is the imaging time interval. It is the integral series of the reference image and the image to be matched. The integral image shift at time t, It is a single integration time. These are the flutter amplitude, frequency, and phase. yes Initial integration time; S2: Perform multispectral image matching based on the known imaging time interval to obtain imaging parallax, and obtain flutter information based on the flutter inversion estimation model; S3: Construct an attitude flutter rotation matrix based on the flutter information, thereby constructing a true attitude rotation matrix containing attitude flutter and obtaining attitude information; S4: Obtain the corrected image through virtual re-imaging based on the attitude information; Step S4 specifically includes: S401: Construct an imaging model using smoothed attitude data, and calculate the corrected image points through orthographic projection. On the average elevation surface Ground point coordinates ; S402: Construct an imaging model using attitude information including flutter information, ground point coordinates. Back projection obtains image points on the original image ; S403: Correct coordinates on the image and points on the original image For the same point, the grayscale information of the image point is obtained through bilinear interpolation; S404: Repeat steps S401-S403 to obtain the grayscale value of each virtual pixel.

2. The satellite flutter inversion estimation and compensation method considering time delay integral effect according to claim 1, characterized in that, The expression for converting platform flutter calculation into integral form is as follows: 。 3. The satellite flutter inversion estimation and compensation method considering time delay integral effect according to claim 1, characterized in that, In step S2, the process of performing multispectral image matching based on the known imaging time interval is as follows: S201: Enhance the image using Wallis filtering to improve image contrast, and perform initial registration by translating the image according to the known imaging time interval; S202: Construct matching windows in the image using the SVD-RANSC subpixel phase correlation method and obtain the subpixel offset between each matching window; S203: Remove erroneous matches from the matching window.

4. The satellite flutter inversion estimation and compensation method considering time delay integral effect according to claim 3, characterized in that, In step S203, removing erroneous matches from the matching window specifically includes: Calculate the normalized cross-correlation coefficient between matching windows and remove matching points with correlation coefficients less than a preset correlation threshold; The disparity of matching point pairs in the same image row is statistically analyzed, and matching points are eliminated according to the three-times standard error principle.

5. A satellite flutter inversion estimation and compensation method considering time delay integral effects according to claim 4, characterized in that, The relevant threshold value is within the range of 0.6-0.

8.

6. The satellite flutter inversion estimation and compensation method considering time delay integral effect according to claim 1, characterized in that, In step S3, the process of constructing the attitude flutter rotation matrix is ​​as follows: Based on the flutter information, the image flutter along the rail and perpendicular to the rail is converted into attitude angle changes using sensor parameters. The attitude flutter rotation matrix is ​​then constructed based on these attitude angle changes, thereby creating a true attitude rotation matrix that includes attitude flutter. The expression for this true attitude rotation matrix is ​​as follows: In the formula, This is the true attitude rotation matrix that includes attitude flutter. For attitude flutter rotation matrix, This indicates the attitude obtained by traditional satellite sensors.

7. A satellite flutter inversion estimation and compensation method considering time delay integration effect according to claim 6, characterized in that, The expression for calculating the attitude angle change is: In the formula, It's the focal length. It's the pixel size. It is the off-axis angle. and These represent the changes in roll and pitch angles caused by flutter. Let t be the image flutter along the track direction at time t. The image flutter in the vertical direction at time t.

8. A satellite flutter inversion estimation and compensation method considering time delay integration effect according to claim 6, characterized in that, The attitude flutter rotation matrix Given a positive definite matrix, a Taylor expansion approximates it as follows when the angle reaches the arcsecond level: In the formula, These represent the three attitude angles at a given moment. It is a unit array. express Cross product.

Images