CN116090203A - Satellite Flutter Inversion Estimation and Compensation Method Considering Time Delay Integration Effect - Google Patents

Abstract

The invention relates to a satellite flutter inversion estimation and compensation method considering the time-delay integration effect, comprising the following steps: according to the TDI effect, the calculation of the platform flutter is converted into an integral form, and a flutter inversion estimation model is constructed; Perform multispectral image matching with known imaging time intervals to obtain imaging parallax, and obtain flutter information according to the flutter inversion estimation model; construct an attitude flutter rotation matrix based on the flutter information, thereby constructing a real attitude rotation matrix including attitude flutter, Attitude information is obtained; a corrected image is obtained through virtual re-imaging according to the attitude information. The method was verified experimentally using the China Resources No. 3 (ZY‑3) satellite image, and the feasibility of the method was verified. Compared with the prior art, the invention has good flutter detection and compensation performance.

Description

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

Technical Field

The present invention relates to the technical field of on-orbit satellite attitude flutter detection, and in particular to a satellite flutter inversion estimation and compensation method considering time delay integration effect.

Background Art

High-resolution optical remote sensing satellites are highly valued by countries around the world, which has promoted the continuous development of high-resolution remote sensing technology. However, attitude flutter limits the geometric quality and imaging performance of high-resolution optical satellites. If the flutter is not compensated in a timely and effective manner, the image quality will be degraded, which will in turn affect image registration and change detection, image stitching and elevation modeling.

At present, many scholars are studying the in-orbit satellite attitude flutter detection method, which is often divided into two categories according to the detection data source: indirect detection method based on non-attitude sensors and direct detection method based on attitude sensors. The indirect method relies on satellite images or mapping products, etc. This method is easy to implement, but factors such as image quality and ground terrain information will affect the accuracy of subsequent matching and feature extraction. It is a passive flutter detection method. The direct detection method relies on high-precision attitude sensors, such as star sensors, gyroscopes or angular displacement trackers. However, due to the limitations of hardware technology, domestic attitude sensors have problems such as low sampling frequency, poor attitude accuracy, and insufficient reliability. At the same time, for satellites that are already in orbit, it is impossible to use additional hardware sensor equipment to detect their flutter. Therefore, it is necessary to detect and compensate for satellite attitude flutter under the conditions of existing attitude measurement equipment without increasing the cost of hardware equipment to improve the accuracy of attitude data.

At present, time-delayed integration charge-coupled devices (TDICCDs) have replaced conventional CCD devices and become the mainstream sensor in high-resolution optical sensor imaging systems. In the multi-level integration imaging process of TDICCD, the vibration offset caused by each level of vibration changes with time, so the vibration effect on the synthetic image is no longer the same as that on the single-level integration image. If the multi-level integration process is ignored, the error of the detected and estimated vibration information will increase with the increase of the number of integration levels and frequency.

Summary of the invention

The purpose of the present invention is to overcome the defects of the above-mentioned prior art and to provide a satellite flutter inversion estimation and compensation method taking into account the time delay integration effect, which has good flutter detection and compensation performance.

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

A satellite flutter inversion estimation and compensation method considering time delay integration effect comprises the following steps:

S1: According to the TDI effect, the calculation of platform flutter is converted into integral form and the flutter inversion estimation model is constructed;

S2: performing multispectral image matching according to a known imaging time interval, obtaining imaging parallax, and obtaining flutter information according to the flutter inversion estimation model;

S3: constructing a posture vibration rotation matrix according to the vibration information, thereby constructing a real posture rotation matrix including the posture vibration, and obtaining the posture information;

S4: Obtaining a corrected image through virtual re-imaging according to the posture information.

Furthermore, the calculation of platform flutter is converted into an integral form as follows:

分别是颤振幅值，频率和相位，t_i是i行初始积分时间。Where d _TDI (t) is the image shift of the integrated image at time t, N ₁ is the image integration series, T is the single integration time, A,f,

are the vibration amplitude, frequency and phase respectively, and _ti is the initial integration time of row i.

Furthermore, the expression of the flutter inversion estimation model is:

Where s(t) represents the imaging parallax, Δt is the imaging time interval, and N ₁ and N ₂ are the integral series of the reference image and the image to be matched.

Furthermore, in step S2, the process of performing multispectral image matching according to the known imaging time interval is specifically as follows:

S201: Use Wallis filtering to enhance the image, improve the image contrast, and translate the image according to the known imaging time interval to achieve initial registration;

S202: constructing a matching window in the image using the SVD-RANSC sub-pixel phase correlation method, and obtaining a sub-pixel offset between each matching window;

S203: Eliminate mismatched points in the matching window.

Furthermore, in step S203, removing the mismatched points in the matching window specifically includes:

Calculate the normalized mutual correlation coefficient between the matching windows and remove the matching points whose correlation coefficient is less than the preset correlation threshold;

The parallax of the matching point pairs in the same image row is counted, and the matching points with large offset are eliminated according to the three-fold mean error principle.

Furthermore, the value of the correlation threshold is within the range of 0.6-0.8.

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

According to the vibration information, the image vibration in the along-track and vertical-track directions is converted into attitude angle changes through sensor parameters, and the attitude vibration rotation matrix constructed according to the attitude angle change is used to construct a real attitude rotation matrix containing attitude vibration. The expression of the real attitude rotation matrix containing attitude vibration is:

R＝R _ss R _si

Where R is the true attitude rotation matrix including attitude flutter, _Rss is the attitude flutter rotation matrix, and _Rsi represents the attitude obtained by the traditional satellite sensor.

Furthermore, the calculation expression of the attitude angle change is:

分别为颤振引起的翻滚角和俯仰角变化，j_across(t)为t时刻沿轨方向的影像颤振，j_along(t)为t时刻垂轨方向的影像颤振。Where f is the focal length, p is the pixel size, β is the off-axis angle, Δω and

are the changes of roll angle and pitch angle caused by the vibration, j _across (t) is the image jitter along the track at time t, and j _along (t) is the image jitter in the vertical track direction at time t.

Furthermore, the attitude flutter rotation matrix R _ss is a positive definite matrix, which is approximated to:

表示某一时刻的三个姿态角，e是自然对数，I₃是单位阵，

表示

叉乘。In the formula,

represents the three attitude angles at a certain moment, e is the natural logarithm, I ₃ is the unit matrix,

express

Cross product.

Furthermore, the step S4 specifically includes:

S401: construct an imaging model using the smoothed posture data, and calculate the ground point coordinates (B, L, H) of the corrected image point (r', c') on the average elevation surface H by orthographic projection;

S402: construct an imaging model using the posture information including the chatter information, and reversely project the ground point coordinates (B, L, H) to obtain an image point (r, c) on the original image;

S403: The coordinate (r', c') on the corrected image and the point (r, c) on the original image correspond to the same point, and the grayscale information of the image point is obtained by bilinear interpolation;

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

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

Based on the imaging principle of TDICCD camera, the present invention constructs a satellite flutter accurate inversion estimation model considering the time delay integration effect, and carries out dense matching of multispectral images to obtain flutter information, updates the attitude according to the estimated flutter information, and finally combines the updated accurate attitude to generate a corrected image using virtual re-imaging to eliminate the influence of satellite flutter. Experiments using the image data and attitude data of the ZY-3 satellite show that the method has good flutter detection and compensation performance.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG1 is a schematic flow chart of a satellite flutter inversion estimation and compensation method considering time delay integration effect provided in an embodiment of the present invention;

FIG2 is a schematic diagram of an experimental image ((a) CCD1, (b) CCD2, (c) CCD3) provided in an embodiment of the present invention;

FIG3 is a schematic diagram of a posture comparison result of a CCD1 provided in an embodiment of the present invention;

FIG4 is a schematic diagram of a posture comparison result of a CCD2 provided in an embodiment of the present invention;

FIG5 is a schematic diagram of a posture comparison result of a CCD3 provided in an embodiment of the present invention;

FIG6 is a schematic diagram of the CCD1 parallax ((a) parallax in the along-track direction and (b) parallax in the perpendicular-track direction) of B1 and B2 before and after vibration compensation provided in an embodiment of the present invention;

FIG7 is a schematic diagram of the CCD2 parallax ((a) parallax in the along-track direction and (b) parallax in the perpendicular-track direction) of B1 and B2 before and after vibration compensation provided in an embodiment of the present invention;

FIG. 8 shows the parallax of the CCD3 of B1 and B2 before and after vibration compensation ((a) parallax in the along-track direction and (b) parallax in the perpendicular-track direction) provided in an embodiment of the present invention.

DETAILED DESCRIPTION

In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings here can be arranged and designed in various different configurations.

Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, it does not require further definition and explanation in the subsequent drawings.

In the description of the present invention, it should be noted that the terms "center", "up", "down", "left", "right", "vertical", "horizontal", "inside", "outside", etc. indicate directions or positional relationships based on the directions or positional relationships shown in the accompanying drawings, or are the directions or positional relationships in which the inventive product is usually placed when in use. They are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and therefore should not be understood as a limitation on the present invention.

It should be noted that the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include one or more of the features. In the description of this application, the meaning of "plurality" is two or more, unless otherwise clearly and specifically defined.

In addition, the terms "horizontal", "vertical" and the like do not mean that the components are required to be absolutely horizontal or suspended, but can be slightly tilted. For example, "horizontal" only 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.

Example 1

As shown in FIG1 , this embodiment provides a satellite flutter inversion estimation and compensation method considering the time delay integration effect, including the following steps:

S1: According to the TDI effect, the calculation of platform flutter is converted into integral form and the flutter inversion estimation model is constructed;

S2: performing multispectral image matching according to a known imaging time interval, obtaining imaging parallax, and obtaining flutter information according to the flutter inversion estimation model;

S3: constructing a posture vibration rotation matrix according to the vibration information, thereby constructing a real posture rotation matrix including the posture vibration, and obtaining the posture information;

S4: Obtaining a corrected image through virtual re-imaging according to the posture information.

Each step is described in detail below.

S1: According to the TDI effect, the calculation of platform flutter is converted into integral form and the flutter inversion estimation model is constructed;

In the multi-level integration process of TDICCD, the image offset caused by platform vibration changes with time. Assuming that the platform vibration is a simple sine wave, the vibration after TDI integration can be expressed as:

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

When the attitude of the satellite is affected by the vibration, there is a geometric displacement caused by the vibration in each image row, and the parallax value of each image row is no longer a constant, and there is a parallax difference value related to the geometric displacement value.

是颤振幅值频率相位，t_i是i行初始积分时间,N₁,N₂是参考影像和待匹配影像的积分级数,T是单一积分时间,d_TDI是积分影像像移，s(t)表示成像视差。Where Δt is the imaging time interval, A,f,

is the vibration amplitude, frequency and phase, _ti is the initial integration time of line i, _N1 , _N2 are the integration series of the reference image and the image to be matched, T is the single integration time, _dTDI is the integrated image shift, and s(t) represents the imaging parallax.

S2: performing multispectral image matching according to a known imaging time interval, obtaining imaging parallax, and obtaining flutter information according to the flutter inversion estimation model;

To obtain the chatter offset, a pixel-by-pixel dense matching method based on sub-pixel phase correlation with singular value decomposition and random sampling consistency (SVD-RANSC) is used.

S201: First, use Wallis filtering to enhance the image and improve the image contrast. According to the known imaging time interval, translate the image to achieve initial registration;

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.

S203: At the same time, possible mismatching points are eliminated, including the following two steps: calculating the normalized mutual correlation coefficient between matching windows, and removing matching points whose correlation coefficient is less than the correlation threshold (0.7 in the experiment); counting the disparity of matching point pairs in the same image row, and eliminating matching points with large offsets according to the three-fold mean error principle.

S3: constructing an attitude flutter rotation matrix according to the flutter information, thereby constructing a true attitude rotation matrix including the attitude flutter, and obtaining attitude information;

Generally speaking, for imaging at or near the sub-satellite point, vertical track flutter is mainly affected by the roll angle, while along-track flutter is mainly affected by the pitch angle. Therefore, the image flutter in the along-track and vertical track directions is converted into attitude angle changes through sensor parameters.

分别为颤振引起的翻滚角和俯仰角变化，进而根据姿轨信息构建姿态颤振旋转矩阵R_ss，表示卫星本体在轨道系下的旋转矩阵。Where f is the focal length, p is the pixel size, and β is the off-axis angle.

They are the changes of roll angle and pitch angle caused by flutter respectively, and then the attitude flutter rotation matrix R _ss is constructed according to the attitude and orbit information, which represents the rotation matrix of the satellite body in the orbit system.

Therefore, combined with the vibration information estimated from the image, the true attitude rotation matrix R including the attitude vibration can be obtained:

R＝R _ss R _si (4)

The rotation matrix R _si represents the attitude obtained by the traditional satellite sensor and does not record high-frequency jitter information.

R _ss is a positive definite matrix. Through Taylor expansion, when the angle is small (arc second level), it can usually be approximated as:

表示某一时刻的三个姿态角，e是自然对数，I₃是单位阵，

表示

叉乘，因此将公式(4)表示为矩阵相乘的形式。

represents the three attitude angles at a certain moment, e is the natural logarithm, I ₃ is the unit matrix,

express

Cross product, so formula (4) is expressed as a matrix multiplication form.

S4: Obtain a corrected image through virtual re-imaging according to the posture information.

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

S401: First, an imaging model is constructed using the smoothed posture data, and the ground point coordinates (B, L, H) of the rectified image point (r', c') on the average elevation surface H are calculated by orthographic projection;

S402: Then, an imaging model is constructed using the posture including the chatter information, and the ground point coordinates (B, L, H) are reversely projected to obtain an image point (r, c) on the original image;

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 coordinates of the image point do not necessarily fall at the center of the pixel, so the grayscale information of the image point is obtained by bilinear interpolation;

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

Experiment and analysis

The image data used in the present invention is a ZY-3 multispectral image, which includes three CCDs, and each CCD image consists of four bands.

Based on the above method, the updated posture data was obtained. In order to verify the correctness of the algorithm, it was compared with the original posture of ZY-3. The results are shown in Figure 2-5. We can find that the results are consistent with the original posture data of different CCDs.

After generating the corrected image, the chatter detection is performed again to check whether the chatter is effectively compensated. Figures 6, 7, and 8 show the parallax before and after chatter compensation. We can clearly see that the period disappears and the parallax is concentrated around 0. After effectively correcting the jitter distortion, the root mean square error in the vertical direction drops from 0.3 pixels to 0.1 pixels, and the root mean square error in the along-track direction drops from 0.2 pixels to 0.1 pixels. The above results show the effectiveness of our method.

in conclusion

Based on the imaging principle of TDICCD camera, the present invention constructs a satellite flutter accurate inversion estimation model considering the time delay integration effect, and carries out dense matching of multispectral images to obtain flutter information, updates the attitude according to the estimated flutter information, and finally combines the updated accurate attitude to generate a corrected image using virtual re-imaging to eliminate the influence of satellite flutter. Experiments using the image data and attitude data of the ZY-3 satellite show that the method has good flutter detection and compensation performance.

The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims (10)

1. The satellite flutter inversion estimation and compensation method considering the time delay integral effect is characterized by comprising the following steps of:

s1: converting the calculation of the platform flutter into an integral form according to the TDI effect, and constructing a flutter inversion estimation model;

s2: carrying out multispectral image matching according to a known imaging time interval, obtaining imaging parallax, and obtaining flutter information according to the flutter inversion estimation model;

s3: constructing a posture flutter rotation matrix according to the flutter information, so as to construct a real posture rotation matrix containing posture flutter and obtain posture information;

s4: and obtaining a corrected image through virtual reimaging according to the attitude information.

2. The method for estimating and compensating satellite chatter as defined in claim 1, wherein said calculating the platform chatter is converted into an integral form by the following expression:

wherein d _TDI (t) is the integral image shift at time t, N ₁ The number of integration stages T for an image is a single integration time, a, f,

respectively the flutter amplitude, frequency and phase, t _i Is the i line initial integration time.

3. The method for satellite flutter inversion estimation and compensation taking into account time delay integration effect according to claim 2, wherein the expression of the flutter inversion estimation model is:

where s (t) represents imaging parallax, Δt is an imaging time interval, N ₁ ,N ₂ Is the number of integration stages of the reference image and the image to be matched.

4. The method for estimating and compensating satellite flutter inversion according to claim 1, wherein in step S2, the process of performing multispectral image matching according to the known imaging time interval is specifically as follows:

s201: enhancing the image by using Wallis filtering, improving the contrast of the image, and translating the image according to the known imaging time interval to realize primary registration;

s202: constructing matching windows in the image piece by utilizing an SVD-RANSC sub-pixel phase correlation method, and obtaining sub-pixel offset between each matching window;

s203: and eliminating mismatching points in the matching window.

5. The method for satellite flutter inversion estimation and compensation according to claim 4, wherein in the step S203, the step of eliminating the mismatching point in the matching window specifically comprises:

calculating normalized cross-correlation coefficients among the matching windows, and removing matching points with the correlation coefficients smaller than a preset correlation threshold;

and counting parallax errors of the matching point pairs of the same image line, and removing the matching point with larger offset according to the error principle of three times.

6. The method for estimating and compensating for satellite flutter as recited in claim 5 wherein said correlation threshold has a value within the range of 0.6-0.8.

7. The satellite flutter inversion estimation and compensation method considering the time delay integration effect according to claim 1, wherein in the step S3, the process of constructing the attitude flutter rotation matrix is specifically as follows:

according to the flutter information, converting the image flutter along the rail and the vertical rail into attitude angle change through sensor parameters, and constructing the attitude flutter rotation matrix according to the attitude angle change, thereby constructing a real attitude rotation matrix containing the attitude flutter, wherein the expression of the real attitude rotation matrix containing the attitude flutter is as follows:

R＝R _ss R _si

wherein R is a true gesture rotation matrix containing gesture vibration, R _ss For the posture flutter rotation matrix, R _si Representing pose obtained by a conventional satellite sensorA state.

8. The method for estimating and compensating for satellite flutter inversion with consideration of time delay integration effect according to claim 7, wherein said calculation expression of attitude angle change is:

where f is the focal length, p is the pixel size, and β is the off-axis angle, Δω, and

roll angle and pitch angle changes caused by chatter, j _across (t) is the image flutter along the track direction at the moment t, j _along And (t) is the image flutter in the vertical track direction at the moment t.

9. The method for satellite flutter inversion estimation and compensation taking into account time delay integration effect according to claim 7, wherein said attitude flutter rotation matrix R _ss For positive definite matrix, through taylor expansion, when the angle reaches the order of the angle second, the approximation is:

wherein θ represents three attitude angles at a certain time, e is a natural logarithm, I ₃ Is a unit array, [ theta ] x]Representing theta cross-multiplication.

10. The method for estimating and compensating satellite flutter inversion according to claim 1, wherein said step S4 comprises:

s401: constructing an imaging model by using the smooth attitude data, and calculating the ground point coordinates (B, L, H) of the corrected image points (r ', c') on the average elevation plane H by orthographic projection;

s402: constructing an imaging model by utilizing attitude information containing flutter information, and obtaining image points (r, c) on an original image by back projection of ground point coordinates (B, L, H);

s403: correcting coordinates (r ', c') on the image and points (r, c) on the original image to correspond to the same point, and acquiring gray information of the image point through bilinear interpolation;

s404: steps S401-S403 are repeated to obtain a gray value for each virtual image point.

Images