CN107966137A - A kind of satellite platform flutter detection method based on TDICCD splice regions image - Google Patents

Abstract

The invention discloses a kind of satellite platform flutter detection method of TDICCD splice regions image, comprise the following steps：1) acquisition of overlay chart picture, overlay chart picture is obtained using time delay integration charge coupling device splice region.2) two width overlay chart pictures are carried out accurate dense Stereo Matching processing, obtain homonymy matching point, it is poor to calculate opposite image space of the same target in two width overlay chart pictures by the calculating of opposite image space difference；3) estimation of satellite platform flutter, it is poor according to opposite image space, estimate satellite platform flutter.The present invention can effectively solve the problems, such as that flutter detection accuracy is low, be only capable of detecting some isolated Frequency points, for satellite platform provide a kind of high accuracy, in wider frequency section all standing flutter detection method, so as to improve the accuracy of satellite attitude acquisition ability and remote sensing image information interpretation.

Description

A satellite platform flutter detection method based on TDICCD stitching region images

technical field

The invention relates to the field of aerospace remote sensing imaging, in particular to a satellite platform flutter detection method based on a time-delayed and integration charge coupled device TDICCD (Time-delayed and Integration Charge Coupled Device TDICCD) stitching area image.

Background technique

Satellite attitude information is very important for the stability control of satellite attitude in the orbit operation stage and the extraction of ground image information in the later stage. As a common phenomenon in the process of satellite operation, flutter has a direct impact on the accuracy of satellite attitude data. With the support of high-performance satellite attitude sensors, remote sensing technology powers such as Europe and the United States can effectively detect flutter within a certain range, thereby obtaining high-precision satellite attitude data. However, the current domestic sensing technology is still immature, and Foreign countries have blocked sales in the high-performance attitude sensor market. The spaceborne attitude sensors currently used in my country have problems such as low measurement bandwidth and poor accuracy, and cannot effectively measure the flutter of satellite platforms, which restricts the development of high-resolution space cameras in my country. Imaging quality and the accuracy of later image information extraction. How to effectively detect the flutter of high-resolution remote sensing satellite platforms without relying on high-precision satellite attitude sensing equipment has become one of the key problems to be solved in the field of space remote sensing imaging in my country.

In recent years, flutter detection methods based on remote sensing images have attracted attention at home and abroad. As the primary condition for flutter detection, the acquisition of remote sensing images with overlapping regions has become the focus of domestic and foreign scholars. The University of Cambridge in the United Kingdom used an area array staring camera to obtain overlapping images. Chang Guang Institute and Wuhan University proposed to add an auxiliary area array imaging sensor on the focal plane to capture overlapping images. The University of Tokyo and Tongji University in Japan used multi-spectral cameras to image different spectral bands The sensor obtains overlapping scene images. Although there are many ways to acquire overlapping images, the flutter detection algorithm is roughly the same. Usually, spectrum analysis technology is used to solve the amplitude, frequency and phase of flutter, and then the trigonometric function is used to fit the flutter curve. This algorithm can concentrate energy It is difficult to detect the main flutter component at the isolated frequency point, but it is difficult to detect the secondary flutter component with lower energy and wider frequency band. In addition, the chatter caused by the unavoidable energy leakage phenomenon in The vibration calculation error cannot be ignored. Therefore, how to solve the problem of low accuracy of flutter detection in the existing technology and can only detect a few isolated frequency points, and propose a new flutter detection method with high precision and full coverage in a wider frequency band, has become the current high-resolution flutter detection method in my country. The urgent task of high-speed remote sensing satellites.

Contents of the invention

The purpose of the present invention is to address the above-mentioned defects in the prior art, and propose a satellite platform flutter detection method based on TDICCD splicing area images, in order to effectively solve the problem of low flutter detection accuracy and only able to detect several isolated frequency points , to provide a flutter detection method with high precision and full coverage in a wide frequency band for the satellite platform, thereby improving the satellite attitude detection capability and the accuracy of remote sensing image information extraction.

The present invention adopts following technical scheme for solving technical problems:

A kind of feature of the satellite platform flutter detection method based on the TDICCD splicing area image of the present invention is to carry out as follows:

Step 1. Acquisition of overlapping images:

Step 1.1, the satellite platform is equipped with a high-resolution space camera, and there are multiple TDICCDs arranged in a staggered arrangement on the focal plane of the high-resolution space camera, and two adjacent TDICCDs are selected arbitrarily, recorded as TDICCD ₁ and TDICCD ₂ , along the TDI direction, the number of rows between TDICCD ₁ and TDICCD ₂ is marked as L, and in the direction perpendicular to the TDI, there are several columns of overlapping pixels at the junction of TDICCD ₁ and TDICCD ₂ to form a splicing area, which is marked as OA ₁ and OA ₂ , the number of overlapping pixel columns in the stitching area OA ₁ and OA ₂ is denoted as C;

The starting moment of any shooting task of the high-resolution spatial camera is recorded as time 0, and the shooting duration is recorded as T _w , where T _w ∈ Q, T _w > 0, and Q represents a set of rational numbers;

Use the stitching areas OA ₁ and OA ₂ to take two time-sharing shots of all targets in the detection area during [0,T _w ], and obtain an image I ₁ =(I ₁ ¹ ,I ₁ ² ,...,I ₁ ^r ,…,I ₁ ^R ) ^T and I ₂ =(I ₂ ¹ ,I ₂ ² ,…,I ₂ ^r ,…,I ₂ ^R ) ^T , where is the total number of rows of the images I ₁ and I ₂ , T _r is the exposure time of each row of images, I ₁ ^r and I ₂ ^r are the rth row images of the images I ₁ and I ₂ respectively, T _r ∈ Q , T _r >0, 1≤r≤R;

Step 1.2, extract the overlapping parts (I ₁ ¹ ,I ₁ ² ,...,I ₁ ^RL ) ^T and (I ₂ ^L+1 ,I ₂ ^L+2 ,...,I ₂ in the images I ₁ and I ₂ ^R ) ^T , respectively denoted as overlapping images O ₁ and O ₂ , the total number of columns of the overlapping images O ₁ and O ₂ is C, and the total number of rows is denoted as M=RL;

Step 2. Calculation of relative imaging position difference:

Step 2.1. Use the image registration method to perform line-by-line matching processing on the overlapping images O ₁ and O ₂ to obtain a set of registration points with the same name P={(p ₁ ¹ ,p ₂ ¹ ),(p ₁ ² ,p ₂ ² ),…,(p ₁ ^m ,p ₂ ^m ),…,(p ₁ ^M ,p ₂ ^M )}, where (p ₁ ^m ,p ₂ ^m ) is the mth of the overlapping images O ₁ and O ₂ A pair of registration points with the same name obtained after image registration, 1≤m≤M, and have: and are the row coordinates and column coordinates of point p ₁ ^m in the overlapping image O ₁ respectively, and are the row coordinates and column coordinates of the point p ₂ ^m in the overlapping image O ₂ respectively,

Step 2.2. Calculate the row coordinate difference and column coordinate difference of two points in each pair of registration points with the same name in the same-name registration point set P respectively, and obtain the relative imaging position difference set S ^|| ={s ^| along the TDI direction ^| (1), s ^|| (2), ..., s ^|| (m), ..., s ^|| (M)} and the set of relative imaging position differences perpendicular to the TDI direction S ^⊥ ＝ {s ^⊥ (1 ), s ^⊥ (2), ..., s ^⊥ (m), ..., s ^⊥ (M)}, where s ^|| (m) is the m-th line of the overlapping image O ₁ and O ₂ after image registration The resulting relative imaging position difference along the TDI direction, and s ^⊥ (m) is the relative imaging position difference perpendicular to the TDI direction obtained after registration of the m-th line of the overlapping images O ₁ and O ₂ , and

Step 2.3: Arrange the set of relative imaging position differences S ^|| in the TDI direction with every continuous L elements starting from the first element, so as to form a two-dimensional matrix S _∑ ^|| of N rows×L columns :

in

Step 2.4: Arrange the set of relative imaging position differences S ^⊥ in the direction perpendicular to the TDI with every continuous L elements starting from the first element, so as to form a two-dimensional matrix S _∑ ^⊥ of N rows×L columns:

in

Step 3. Estimation of satellite platform flutter:

Step 3.1, define the flutter along the TDI direction as:

Among them, g ^|| (v L+u) is the flutter along the TDI direction at time (v L+u)×T _r , g _|| (v L+u)∈Q, u∈Z, v∈ Z, 1≤u≤L, 0≤v≤N, Z represents a set of integers;

Define flutter in the direction perpendicular to TDI as:

Among them, g ^⊥ (v L+u) is the flutter perpendicular to the TDI direction at time (v L+u)×T _r , g _⊥ (v L+u)∈Q, u∈Z, v∈Z , 1≤u≤L, 0≤v≤N;

Define flutter as G _∑ ＝{g _ij } _(N+1)×L , where is the unit vector along the TDI direction, is a unit vector perpendicular to the TDI direction, 1≤i≤N+1, 1≤j≤L;

Define the elements in the first row [g ^|| (1), g ^|| (2), ..., g ^|| (L)] in the flutter G ^∑|| along the TDI direction as the flutter sub-blocks along the TDI direction G ₁ ^|| ;

Define the first row of elements [g ^⊥ (1), g ^⊥ (2), ..., g ^⊥ (L)] in the flutter G _∑ ^⊥ perpendicular to the TDI direction to be the flutter sub-block G ₁ perpendicular to the TDI direction ^⊥ ;

Define elements from the second row to the N+1th row in the flutter G _∑ ^|| along the TDI direction

is the flutter sub-block G ₂ ^|| along the TDI direction;

Define elements from the second row to the N+1th row in the flutter G _∑ ^⊥ perpendicular to the TDI direction

is the flutter sub-block G ₂ ^⊥ perpendicular to the TDI direction;

Step 3.2. Incorporate the on-orbit operation parameters of the satellite in the time period [T _r , L·T _r ] in the shooting task into the sub-satellite point-to-earth imaging model, and obtain the flutter sub-block G ₁ ^{| |} Low frequency components below 1Hz and the low-frequency component of the flutter sub-block G ₁ ^⊥ below 1Hz in the direction perpendicular to the TDI in Indicates the low-frequency component of flutter below 1 Hz along the TDI direction at time u×T _r , Indicates the low-frequency component of flutter below 1 Hz in the direction perpendicular to TDI at u×T _r ,

Step 3.3, use the flutter sub-block G ₁ ^|| and its low frequency component along the TDI direction The sum of the squares of the relative residuals between is used as the objective function H(G ₁ ^|| ) along the TDI direction, as shown in formula (1):

In formula (1), g ^|| (i) represents the flutter along the TDI direction at time i×T _r , Indicates the low-frequency component of flutter below 1 Hz along the TDI direction at j×T _r time;

Take the flutter sub-block G ₁ ^⊥ and its low frequency component in the direction perpendicular to TDI The sum of the squares of the relative residuals between is used as the objective function H(G ₁ ^⊥ ) in the direction perpendicular to the TDI, as shown in formula (2):

In formula (2), g ^⊥ (i) represents the flutter in the direction perpendicular to TDI at time i×T _r , Indicates the low-frequency component of flutter below 1 Hz in the direction perpendicular to TDI at j×T _r ;

Step 3.4, use the optimization algorithm to solve the minimum points of the objective functions H(G ₁ ^|| ) and H(G ₁ ^⊥ ) respectively and and respectively as the estimated values of the flutter sub-block G ₁ ^|| along the TDI direction and the flutter sub-block G ₁ ^⊥ perpendicular to the TDI direction, as shown in formulas (3) and (4):

Step 3.5. According to the push-broom imaging principle of the space camera and the imaging characteristics of the TDICCD stitching area, the flutter detection models along the TDI direction and perpendicular to the TDI direction are established as shown in formulas (5) and (6):

G ₂ ^|| ＝A·G ₁ ^|| +B·S _∑ ^|| (5)

G ₂ ^⊥ ＝A·G ₁ ^⊥ +B·S _∑ ^⊥ (6)

In formula (5) and formula (6), A and B both represent the coefficient matrix, and

Step 3.6. According to the flutter detection model along the TDI direction shown in formula (5), use formula (7) to calculate the estimated value of the flutter sub-block G ₂ ^|| along the TDI direction

According to the flutter detection model perpendicular to the TDI direction shown in formula (6), use formula (8) to calculate the estimated value of the flutter sub-block G ₂ ^⊥ perpendicular to the TDI direction

Step 3.7, according to the definition of flutter along the TDI direction, use formula (9) to synthesize the estimated value of the flutter G _∑ ^|| along the TDI direction

According to the definition of flutter in the direction perpendicular to TDI, the estimated value of flutter G _∑ ^⊥ in the direction perpendicular to TDI is synthesized by formula (10)

Compared with prior art, the beneficial effect of the present invention is:

1. In view of the low accuracy of flutter detection in the prior art and the problem that only a few isolated frequency points can be detected, the present invention proposes a satellite platform flutter detection method based on TDICCD splicing area images, which realizes the detection without relying on high-precision attitude transmission. In the case of sensing equipment and auxiliary imaging sensors, the flutter of satellite platforms can be detected with high precision and full coverage in a wide frequency band.

2. The flutter frequency components of the satellite platform are complex. The existing method uses the superposition of several trigonometric functions to fit the flutter, and can extract the main components of the flutter with concentrated energy, which are often several isolated frequency points. The secondary components occupying a large continuous frequency band are difficult to detect, resulting in a large gap between the flutter detection results and the real situation. By establishing a flutter detection model, the present invention clarifies the quantitative relationship between flutter and the relative imaging position difference, avoids the theoretical error introduced by the prior art for approximate fitting of flutter, and is a strong support for improving the accuracy of flutter detection .

3. Existing methods usually rely on spectrum analysis techniques to obtain the frequency, amplitude, and phase of trigonometric functions, and the unavoidable energy leakage phenomenon in spectrum analysis techniques will introduce significant frequency, amplitude, and phase estimation errors, eventually leading to chatter Probing errors cannot be ignored. Based on the flutter detection model, the present invention adopts the time-domain analysis technology to directly obtain the time-domain value of the flutter of the satellite platform at the time of shooting, avoiding the flutter detection error introduced by energy leakage in the prior art, and further improving the accuracy of flutter detection Spend.

4. The flutter at different positions on the satellite is different, even the flutter at the same position is constantly changing with time, and the flutter at the main imaging sensor on the focal plane of the space camera and at the moment of image capture is the most concerned by people . For non-staring imaging and non-multispectral cameras, the existing technology often detects flutter by adding an auxiliary area array imaging sensor near the main imaging sensor on the focal plane, and the detection result is at the auxiliary imaging sensor on the focal plane instead of the main There are differences between the two in the chatter at the imaging sensor, especially for space cameras with relatively small apertures, the chatter at different positions on the focal plane differs more significantly. In addition, the two types of imaging sensors work in different ways and cannot be exposed synchronously. The above two factors eventually lead to deviations between the flutter detection results and the required detection position and time. The present invention uses the image taken by the splicing area of the main imaging sensor for flutter detection, and the detection result is the flutter at the main imaging sensor and at the moment of image shooting, which overcomes the flutter detection result and the required detection position in the prior art The problem of deviation from time to time has laid a foundation for further improving the accuracy of flutter detection.

5. The frequency band occupied by satellite platform flutter is relatively wide. Existing methods can detect flutter at isolated frequency points where energy is concentrated, but full coverage detection in a wider frequency band has not yet been achieved. The present invention can increase the sampling frequency of flutter by increasing the intensity of image registration, thereby increasing the bandwidth of flutter detection, and the bandwidth of flutter detection can be increased up to 1/2 of the reciprocal of the exposure time of each row of images in the camera .

Description of drawings

Fig. 1 is a flow chart of the method of the present invention;

Fig. 2 is a schematic diagram of the method of taking points evenly at equal intervals in the present invention;

Fig. 3 is the sub-pixel registration flowchart of the combination of coarse registration and fine registration in the present invention;

Fig. 4 is a schematic diagram of imaging of two adjacent TDICCD splicing regions in the prior art;

Fig. 5 is a schematic diagram of the focal plane arrangement of a high-resolution TDICCD visible light camera in the prior art;

Fig. 6 is a visible light panchromatic remote sensing image captured by a high-resolution TDICCD visible light camera in orbit in the prior art;

Fig. 7a is a diagram of the flutter detection results of the satellite platform perpendicular to the TDI direction of the present invention;

Fig. 7b is a diagram of the flutter detection results of the satellite platform along the TDI direction of the present invention;

Figure 7c is a power spectral density diagram of the flutter detection results of the satellite platform perpendicular to the TDI direction of the present invention;

Fig. 7d is a power spectral density diagram of the flutter detection results of the satellite platform along the TDI direction according to the present invention.

Detailed ways

In the present embodiment, a kind of satellite platform flutter detection method based on TDICCD splicing area image, as shown in Figure 1, is carried out as follows:

Step 1. Acquisition of overlapping images:

Step 1.1. The satellite platform is equipped with a high-resolution space camera. On the focal plane of the high-resolution space camera, there are usually multiple staggered and spliced TDICCDs. Randomly select two adjacent TDICCDs, which are recorded as TDICCD ₁ and TDICCD ₂ , in the direction along the TDI, the number of rows between TDICCD ₁ and TDICCD ₂ is denoted as L, and in the direction perpendicular to the TDI, there are several columns of overlapping pixels at the junction of TDICCD ₁ and TDICCD ₂ to form a splicing area, which is denoted as OA ₁ and OA ₂ , the number of overlapping pixel columns in the mosaic area OA ₁ and OA ₂ is denoted as C;

The starting moment of any shooting task of the high-resolution spatial camera is recorded as time 0, and the shooting duration is recorded as T _w , where T _w ∈ Q, T _w >0; Q represents a set of rational numbers;

Use the splicing area OA ₁ and OA ₂ to take two time-sharing shots of all targets in the detection area during [0, T _w ], and obtain the image I ₁ = (I ₁ ¹ ,I ₁ ² ,...,I ₁ ^r ,…,I ₁ ^R ) ^T and I ₂ =(I ₂ ¹ ,I ₂ ² ,…,I ₂ ^r ,…,I ₂ ^R ) ^T , where is the total number of rows of images I ₁ and I ₂ , T _r is the exposure time of each row of images, I ₁ ^r and I ₂ ^r are the r-th row images of images I ₁ and I ₂ respectively, T _r ∈ Q, T _r >0,1≤r≤R;

Step 1.2, extract the overlapping parts (I ₁ ¹ ,I ₁ ² ,…,I ₁ ^RL ) ^T and (I ₂ ^L+1 ,I ₂ ^L+2 ,…,I ₂ ^R ) in images I ₁ and I ₂ ^T is recorded as overlapping images O ₁ and O ₂ respectively, the total number of columns of overlapping images O ₁ and O ₂ is C, the total number of rows is (RL), and M=RL;

Step 2. Calculation of relative imaging position difference:

Step 2.1. Perform line-by-line matching processing on the overlapped images O ₁ and O ₂ by using an image registration method. The specific process of line-by-line matching processing is shown in Figure 2. _In each line of images of the overlapping image O1, several registration points are selected in an evenly spaced manner, and a template of a certain size is selected as the reference template with each registration point as the center. Select the template at the corresponding position in the overlapping image _O2 as the template to be registered, and use the image registration method to match the two templates to obtain a pair of registration points with the same name. Perform the same processing on other registration points in the same row of images, and take the arithmetic mean value of the registration results of all the registration points in the same row as the registration result of the entire row of images.

The present invention adopts a sub-pixel registration method combining coarse registration and fine registration. The specific registration process is shown in Fig. The quasi-templates are u ₁ (m,n) and u ₂ (m,n) respectively, and the corresponding Fourier transforms are:

where ξ and η are spatial frequency variables. The cross power spectrum of template u ₁ (m,n) and u ₂ (m,n) is:

where U ₁ ^* (ξ,η) is the complex conjugate of U ₁ (ξ,η). Perform inverse Fourier transform on the cross-power spectrum C _ps (ξ,η) to obtain the two-dimensional impulse function δ(i,j), as shown in formula (15), the peak point (i _p ,j _p ) is the pixel-level coarse registration result.

IFFT2{C _ps (ξ,η)}＝IFFT2{e j2π ^(iξ+jη) }＝δ(i,j) (3)

The basic principle of the fine registration method based on cross-correlation two-dimensional surface fitting is: according to the rough registration results, select templates u _{1 (} m,n) and u _{2 (} mi _p , nj _p ), using the normalized cross-correlation coefficient as an evaluation index to measure the similarity between two templates, the expression is as follows:

Where C(u,v) is the normalized cross-correlation coefficient when the center points of the two templates are (u,v) apart. Use the least squares method to solve the expression of the quadratic surface of the cross-correlation coefficient. The peak of the surface is the maximum point of the normalized cross-correlation coefficient C(u,v), and the center points of the corresponding two templates are the same-name matching on time.

Use the above image registration method to perform line-by-line matching processing on overlapping images O ₁ and O ₂ , and obtain a set of registration points with the same name P={(p ₁ ¹ ,p ₂ ¹ ),(p ₁ ² ,p ₂ ² ),…, (p ₁ ^m ,p ₂ ^m ),,,…,(p ₁ ^M ,p ₂ ^M )}, where (p ₁ ^m ,p ₂ ^m ) is obtained after registration of the mth row of the overlapping images O ₁ and O ₂ A pair of registration points with the same name, 1≤m≤M, and have: and are the row coordinates and column coordinates of point p ₁ ^m in the overlapping image O ₁ respectively, and are the row coordinates and column coordinates of the point p ₂ ^m in the overlapping image O ₂ respectively,

Step 2.2, respectively calculate the row coordinate difference and column coordinate difference of two points in each pair of registration points with the same name in the same-name registration point set P, and obtain the relative imaging position difference set S ^|| = {s ^|| ( 1), s ^|| (2), ..., s ^|| (m), ..., s ^|| (M)} and the set of relative imaging position differences perpendicular to the TDI direction S ^⊥ ＝{s ^⊥ (1), s ^⊥ (2), ..., s ^⊥ (m), ..., s ^⊥ (M)}, where s ^|| (m) is obtained after registration of the m-th line of the overlapping images O ₁ and O ₂ The relative imaging position difference along the TDI direction, and s ^⊥ (m) is the relative imaging position difference perpendicular to the TDI direction obtained after registration of the m-th line of the overlapping images O ₁ and O ₂ , and

Step 2.3. Arrange the set of relative imaging position differences S ^|| along the TDI direction with every continuous L elements starting from the first element, so as to form a two-dimensional matrix S _∑ ^|| of N rows×L columns:

in

Step 2.4. Arrange the relative imaging position difference set S ^⊥ in the direction perpendicular to the TDI direction and arrange every continuous L elements in a row from the first element, thus forming a two-dimensional matrix S _∑ ^⊥ of N rows×L columns:

in

Step 3. Estimation of satellite platform flutter:

Step 3.1, define the flutter along the TDI direction as

Among them, g ^|| (v·L+u) is the flutter along the TDI direction at time (v·L+u)×T _r , g _|| (v·L+u)∈Q, u∈Z,v∈ Z, 1≤u≤L, 0≤v≤N, Z represents a set of integers;

Define flutter in the direction perpendicular to TDI as

Among them, g ^⊥ (v L+u) is the flutter perpendicular to the TDI direction at time (v L+u)×T _r , g _⊥ (v L+u)∈Q, u∈Z, v∈Z , 1≤u≤L, 0≤v≤N;

Define flutter as G _∑ ＝{g _ij } _(N+1)×L , where is the unit vector along the TDI direction, is a unit vector perpendicular to the TDI direction, 1≤i≤N+1, 1≤j≤L;

Define the elements in the first row [g ^|| (1), g ^|| (2), ..., g ^|| (L)] in the flutter G ^∑|| along the TDI direction as the flutter sub-blocks along the TDI direction G ₁ ^|| ;

Define the first row of elements [g ^⊥ (1), g ^⊥ (2), ..., g ^⊥ (L)] in the flutter G _∑ ^⊥ perpendicular to the TDI direction to be the flutter sub-block G ₁ perpendicular to the TDI direction ^⊥ ;

Define elements from the second row to the N+1th row in the flutter G _∑ ^|| along the TDI direction

is the flutter sub-block G ₂ ^|| along the TDI direction;

Define elements from the second row to the N+1th row in the flutter G _∑ ^⊥ perpendicular to the TDI direction

is the flutter sub-block G ₂ ^⊥ perpendicular to the TDI direction;

Step 3.2. Bring the on-orbit operation parameters of the satellite in the time period [T _r , L T _r ] into the document "WangJiaqi, Yu Ping, Yan Changxiang, Ren Jianyue, Hebin. Space optical remote sensor image motion velocity vector computational modeling, error budget and synthesis[J].Chinese Optics Letters,2005,(07):414-417." In the sub-satellite point-to-earth imaging model given, the flutter sub-block G ₁ ^|| Low frequency components below 1Hz and the low-frequency component of the flutter sub-block G ₁ ^⊥ below 1Hz in the direction perpendicular to the TDI in Indicates the low-frequency component of flutter below 1 Hz along the TDI direction at time u×T _r , Indicates the low-frequency component of flutter below 1 Hz in the direction perpendicular to TDI at u×T _r ,

Step 3.3. Flutter sub-block G ₁ ^|| and its low-frequency components along the TDI direction The sum of the squares of the relative residuals between is used as the objective function H(G ₁ ^|| ) along the TDI direction, as shown in formula (5):

In formula (5), g ^|| (i) represents the flutter along the TDI direction at time i×T _r , Indicates the low-frequency component of flutter below 1 Hz along the TDI direction at j×T _r time;

Flutter sub-block G ₁ ^⊥ and its low frequency component in the direction perpendicular to TDI The sum of the squares of the relative residuals between is used as the objective function H(G ₁ ^⊥ ) in the direction perpendicular to the TDI, as shown in formula (6):

In formula (6), g ^⊥ (i) represents the flutter in the direction perpendicular to TDI at time i×T _r , Indicates the low-frequency component of flutter below 1 Hz in the direction perpendicular to TDI at j×T _r ;

Step 3.4, use the conjugate gradient method in the optimization algorithm with the convergence rate and calculation amount in the middle to solve the minimum points of the objective functions H(G ₁ ^|| ) and H(G ₁ ^⊥ ) respectively and and respectively as the estimated values of the flutter sub-block G ₁ ^|| along the TDI direction and the flutter sub-block G ₁ ^⊥ perpendicular to the TDI direction, as shown in formulas (7) and (8):

Step 3.5. Deduce the flutter detection model according to the push-broom imaging principle of the space camera and the imaging characteristics of the TDICCD stitching area. The specific derivation process is as follows:

According to the principle of push-broom imaging, the same target will be imaged twice in time by the splicing of two adjacent TDICCDs, while the flutter of the satellite platform is constantly changing with time, and the splicing areas OA ₁ and OA ₂ will image the same target twice in time The flutter states of the corresponding satellite platforms may be different, which leads to differences in the imaging positions of the target in the stitching area OA ₁ and OA ₂ , as shown in Figure 4, taking the direction perpendicular to the TDI as an example, the relative imaging position difference is related to The vibration relationship can be expressed as:

s ^⊥ (t) = g ^⊥ (t+Δt)-g ^⊥ (t) (9)

In formula (9), g ^⊥ (t) is the flutter in the direction perpendicular to the TDI at time t, Δt is the time interval between the two imaging of the same target in the stitching area OA ₁ and OA ₂ , and Δt=L×T _r , L is the number of rows between two adjacent TDICCDs, T _r is the exposure time of each row of images, and s ^⊥ (t) is the relative imaging position difference perpendicular to the TDI direction at time t.

The formula (9) is discretized and sorted by taking the exposure time T _r of each row of images as the sampling period to obtain:

g ^⊥ (m+L) = g ^⊥ (m) + s ^⊥ (m) (10)

In formula (10), g ^⊥ (m) is the flutter perpendicular to the direction of TDI at time m×T _r , and g ^⊥ (m+L) is the flutter perpendicular to the direction of TDI at time (m+L)×T _r vibration, s ^⊥ (m) is the relative imaging position difference perpendicular to the TDI direction at m×T _r time, 1≤m≤M.

For the flutter g ^⊥ (m+L) perpendicular to the TDI direction, m=1,2,...,M starts from the first element and every consecutive L elements are divided into groups, expressed as:

In formula (11)

According to the definitions of the flutter sub-block G ₁ ^⊥ perpendicular to the TDI direction, the flutter sub-block G ₂ ^⊥ perpendicular to the TDI direction, and the relative imaging position difference S _∑ ^⊥ perpendicular to the TDI direction, formula (11) can be expressed as :

Let the coefficient matrix coefficient matrix Then formula (12) can be expressed as:

G ₂ ^⊥ ＝A·G ₁ ^⊥ +B·S _∑ ^⊥ (13)

Equation (13) is the flutter detection model perpendicular to the TDI direction. Similarly, the derivation of the flutter detection model along the TDI direction is shown in Equation (14):

G ₂ ^|| ＝A·G ₁ ^|| +B·S _∑ ^|| (14)

The meanings of G ₁ ^|| , G ₂ ^|| and S _∑ ^|| in formula (14) are the same as those defined in the description, and the values of coefficient matrices A and B are the same as above;

Step 3.6. According to the chatter detection model along the TDI direction shown in formula (14), use the estimated value of the chatter sub-block G ₁ ^|| along the TDI direction and the relative imaging position difference matrix S _∑ ^|| along the TDI direction to calculate the estimated value of the flutter sub-block G ₂ ^|| along the TDI direction As shown in formula (15):

According to the flutter detection model perpendicular to the TDI direction shown in equation (13), the estimated value of the flutter sub-block G ₁ ^⊥ in the direction perpendicular to the TDI is used and the relative imaging position difference matrix S _∑ ^⊥ in the direction perpendicular to the TDI to calculate the estimated value of the flutter sub-block G ₂ ^⊥ in the direction perpendicular to the TDI As shown in formula (16):

Step 3.7, according to the definition of chatter along the TDI direction, the estimated value of the chatter sub-block G ₁ ^|| along the TDI direction and the estimated value of the flutter sub-block G ₂ ^|| along the TDI direction Synthetic estimates of flutter G _∑ ^|| along the TDI direction As shown in formula (17):

According to the definition of flutter in the direction perpendicular to TDI, the estimated value of the flutter sub-block G ₁ ^⊥ in the direction perpendicular to TDI and the estimated value of the flutter sub-block G ₂ ^⊥ perpendicular to the TDI direction Estimated values of the synthesized flutter G _∑ ^⊥ in the direction perpendicular to the TDI As shown in formula (18):

The method proposed by the present invention will be described below by taking a certain high-resolution TDICCD camera in my country as an example. The arrangement of the focal plane of the camera is shown in Figure 5. Five TDICCDs are arranged in a staggered manner. In the direction along the TDI, the number of rows between any two adjacent TDICCDs is L=3028. In the direction perpendicular to the integration of the TDI, The number of pixel columns in the stitching area of any two adjacent TDICCDs is C=40, the duration of one shooting task is T _w =30s, and the exposure time of each row of images is T _r =73μs. In this experiment, the remote sensing images taken by the camera at 38.8436° to 39.4515° north latitude during the orbit ascent stage _were taken by TDICCD ₃ and TDICCD ₄ , and the method of the present invention was used to detect the satellite platform flutter during image shooting _. A series of images, as shown in Figure 6, the left side is part of the remote sensing image taken by TDICCD ₃ , the right side is part of the remote sensing image taken by TDICCD ₄ , and the dotted line is part of the overlapping images of the stitching area.

Step 1. Acquisition of overlapping images:

Image I ₁ =(I ₁ ¹ , _{I 1 2 ,} ..., _{I 1 411520 )} ^{T and} I ² = ⁽ I ₂ ¹ ,I _{2 2} ,..., I ₂ ⁴¹¹⁵²⁰ ) ^T , extract the overlapping parts (I ₁ ¹ ,I ₁ ² ,…,I ₁ ⁴⁰⁸⁴⁹² ) ^T and (I ₂ ³⁰²⁹ ,I ₂ ³⁰³⁰ ,…,I ₂ ⁴¹¹⁵²⁰ ) ^T in images I ₁ and I ₂ , respectively denoted as overlapping images O ₁ and O ₂ ;

Step 2. Calculation of relative imaging position difference:

The overlapped images O ₁ and O ₂ are matched line by line by image registration method. The line-by-line matching process is specifically as follows: In the overlapping image O1, each line _of the image has 40 pixels in total, and the 10th, 20th, and 30th pixels are respectively selected as registration points, and the 8 lines centered on each registration point *The 8-column neighborhood is the reference template, select the module at the corresponding position in the overlapping image _O2 as the template to be registered, and use the sub-pixel registration method combining coarse registration and fine registration to match the two templates , to obtain a pair of registration points with the same name, take the arithmetic mean of the registration results of the three registration points in the same image as the registration result of the entire row of images, and follow the above steps for the images from the 11th row to the _408482th row of the overlapping image O1 Perform line-by-line matching processing to obtain a set of registration points with the same name. Calculate the row coordinate difference and column coordinate difference of two points in each pair of registration points with the same name in the same-name registration point set respectively, and obtain the relative imaging position difference set along the TDI direction and the relative imaging position difference set perpendicular to the TDI direction;

Step 3. Estimation of satellite platform flutter:

Bring the on-orbit operation parameters of the satellite into the sub-satellite point-to-earth imaging model in the time period of [73μs, 3028×73μs] in the shooting task, and obtain the low-frequency components of the flutter sub-block G ₁ ^|| below 1 Hz along the TDI direction and the low-frequency component of the flutter sub-block G ₁ ^⊥ below 1Hz in the direction perpendicular to the TDI Dither sub-block G ₁ ^|| and its low-frequency component along the TDI direction The sum of the squares of the relative residuals between is used as the objective function H(G ₁ ^|| ) along the TDI direction to flutter the sub-block G ₁ ^⊥ and its low-frequency components perpendicular to the TDI direction The sum of the squares _of the relative residuals between is used as the objective function H(G ₁ ^⊥ ₎ perpendicular to the ^TDI direction, and the conjugate gradient method is used to ^solve the minimum value points and and respectively as the estimated values of the flutter sub-block G ₁ ^|| along the TDI direction and the flutter sub-block G ₁ ^⊥ perpendicular to the TDI direction, according to the flutter detection model along the TDI direction, using the flutter sub-block G ₁ along the TDI direction Estimated value of ^|| and the relative imaging position difference matrix S _∑ ^|| along the TDI direction to calculate the estimated value of the flutter sub-block G ₂ ^|| along the TDI direction According to the flutter detection model perpendicular to the TDI direction, using the estimated value of the flutter sub-block G ₁ ^⊥ perpendicular to the TDI direction and the relative imaging position difference matrix S _∑ ^⊥ in the direction perpendicular to the TDI to calculate the estimated value of the flutter sub-block G ₂ ^⊥ in the direction perpendicular to the TDI The estimated value of the flutter sub-block G ₁ ^|| along the TDI direction and the estimated value of the flutter sub-block G ₂ ^|| along the TDI direction Synthetic estimates of flutter G _∑ ^|| along the TDI direction The estimated value of the flutter sub-block G ₁ ^⊥ in the direction perpendicular to the TDI and the estimated value of the flutter sub-block G ₂ ^⊥ perpendicular to the TDI direction Estimated values of the synthesized flutter G _∑ ^⊥ in the direction perpendicular to the TDI

Flutter detection results are shown in Figure 7a, Figure 7b, Figure 7c, and Figure 7d. In the direction perpendicular to the TDI, there are vibration peak points at 0.24Hz, 4.54Hz, 9.08Hz, 13.53Hz, 18Hz and 22.48Hz; In the direction, in addition to the above peak point, there is a vibration peak at 0.13Hz. It shows that the method of the invention can effectively detect the flutter of the satellite platform.

For specific implementation, reference may be made to corresponding descriptions of the above methods.

Claims (1)

1. a satellite platform flutter detection method based on TDICCD stitching area image, its feature is carried out as follows: Step 1. Acquisition of overlapping images: Step 1.1, the satellite platform is equipped with a high-resolution space camera, and there are multiple TDICCDs arranged in a staggered arrangement on the focal plane of the high-resolution space camera, and two adjacent TDICCDs are selected arbitrarily, recorded as TDICCD ₁ and TDICCD ₂ , along the TDI direction, the number of rows between TDICCD ₁ and TDICCD ₂ is marked as L, and in the direction perpendicular to the TDI, there are several columns of overlapping pixels at the junction of TDICCD ₁ and TDICCD ₂ to form a splicing area, which is marked as OA ₁ and OA ₂ , the number of overlapping pixel columns in the stitching area OA ₁ and OA ₂ is denoted as C; The starting moment of any shooting task of the high-resolution spatial camera is recorded as time 0, and the shooting duration is recorded as T _w , where T _w ∈ Q, T _w > 0, and Q represents a set of rational numbers; Use the stitching areas OA ₁ and OA ₂ to take two time-sharing shots of all targets in the detection area during [0,T _w ], and obtain an image I ₁ =(I ₁ ¹ ,I ₁ ² ,...,I ₁ ^r ,…,I ₁ ^R ) ^T and I ₂ =(I ₂ ¹ ,I ₂ ² ,…,I ₂ ^r ,…,I ₂ ^R ) ^T , where is the total number of rows of the images I ₁ and I ₂ , T _r is the exposure time of each row of images, I ₁ ^r and I ₂ ^r are the rth row images of the images I ₁ and I ₂ respectively, T _r ∈ Q , T _r >0, 1≤r≤R; Step 1.2, extract the overlapping parts (I ₁ ¹ ,I ₁ ² ,...,I ₁ ^RL ) ^T and (I ₂ ^L+1 ,I ₂ ^L+2 ,...,I ₂ in the images I ₁ and I ₂ ^R ) ^T , respectively denoted as overlapping images O ₁ and O ₂ , the total number of columns of the overlapping images O ₁ and O ₂ is C, and the total number of rows is denoted as M=RL; Step 2. Calculation of relative imaging position difference: Step 2.1. Use the image registration method to perform line-by-line matching processing on the overlapping images O ₁ and O ₂ to obtain a set of registration points with the same name P={(p ₁ ¹ ,p ₂ ¹ ),(p ₁ ² ,p ₂ ² ),…,(p ₁ ^m ,p ₂ ^m ),…,(p ₁ ^M ,p ₂ ^M )}, where (p ₁ ^m ,p ₂ ^m ) is the mth of the overlapping images O ₁ and O ₂ A pair of registration points with the same name obtained after image registration, 1≤m≤M, and have: and are the row coordinates and column coordinates of point p ₁ ^m in the overlapping image O ₁ respectively, and are the row coordinates and column coordinates of the point p ₂ ^m in the overlapping image O ₂ respectively, Step 2.2. Calculate the row coordinate difference and column coordinate difference of two points in each pair of registration points with the same name in the same-name registration point set P respectively, and obtain the relative imaging position difference set S ^|| ={s ^| along the TDI direction ^| (1), s ^|| (2), ..., s ^|| (m), ..., s ^|| (M)} and the set of relative imaging position differences perpendicular to the TDI direction S ^⊥ ＝ {s ^⊥ (1 ), s ^⊥ (2), ..., s ^⊥ (m), ..., s ^⊥ (M)}, where s ^|| (m) is the m-th line of the overlapping image O ₁ and O ₂ after image registration The resulting relative imaging position difference along the TDI direction, and s ^⊥ (m) is the relative imaging position difference perpendicular to the TDI direction obtained after registration of the m-th line of the overlapping images O ₁ and O ₂ , and Step 2.3: Arrange the set of relative imaging position differences S ^|| in the TDI direction with every continuous L elements starting from the first element, so as to form a two-dimensional matrix S _∑ ^|| of N rows×L columns : in Step 2.4: Arrange the set of relative imaging position differences S ^⊥ in the direction perpendicular to the TDI with every continuous L elements starting from the first element, so as to form a two-dimensional matrix S _∑ ^⊥ of N rows×L columns: in Step 3. Estimation of satellite platform flutter: Step 3.1, define the flutter along the TDI direction as: Among them, g ^|| (v L+u) is the flutter along the TDI direction at time (v L+u)×T _r , g _|| (v L+u)∈Q, u∈Z, v∈ Z, 1≤u≤L, 0≤v≤N, Z represents a set of integers; Define flutter in the direction perpendicular to TDI as: Among them, g ^⊥ (v L+u) is the flutter perpendicular to the TDI direction at time (v L+u)×T _r , g _⊥ (v L+u)∈Q, u∈Z, v∈Z , 1≤u≤L, 0≤v≤N; Define flutter as G _∑ ＝{g _ij } _(N+1)×L , where is the unit vector along the TDI direction, is a unit vector perpendicular to the TDI direction, 1≤i≤N+1, 1≤j≤L; Define the first row of elements [g ^|| (1), g ^|| (2), ..., g ^|| (L)] in the flutter G _∑ ^|| along the TDI direction as the flutter sub-block along the TDI direction G ₁ ^|| ; Define the first row of elements [g ^⊥ (1), g ^⊥ (2), ..., g ^⊥ (L)] in the flutter G _∑ ^⊥ perpendicular to the TDI direction to be the flutter sub-block G ₁ perpendicular to the TDI direction ^⊥ ; Define elements from the second row to the N+1th row in the flutter G _∑ ^|| along the TDI direction is the flutter sub-block G ₂ ^|| along the TDI direction; Define elements from the second row to the N+1th row in the flutter G _∑ ^⊥ perpendicular to the TDI direction is the flutter sub-block G ₂ ^⊥ perpendicular to the TDI direction; Step 3.2. Incorporate the on-orbit operation parameters of the satellite in the time period [T _r , L·T _r ] in the shooting task into the sub-satellite point-to-earth imaging model, and obtain the flutter sub-block G ₁ ^{| |} Low frequency components below 1Hz and the low-frequency component of the flutter sub-block G ₁ ^⊥ below 1Hz in the direction perpendicular to the TDI in Indicates the low-frequency component of flutter below 1 Hz along the TDI direction at time u×T _r , Indicates the low-frequency component of flutter below 1 Hz in the direction perpendicular to TDI at u×T _r , Step 3.3, use the flutter sub-block G ₁ ^|| and its low frequency component along the TDI direction The sum of the squares of the relative residuals between is used as the objective function H(G ₁ ^|| ) along the TDI direction, as shown in formula (1): In formula (1), g ^|| (i) represents the flutter along the TDI direction at time i×T _r , Indicates the low-frequency component of flutter below 1 Hz along the TDI direction at j×T _r time; Take the flutter sub-block G ₁ ^⊥ and its low frequency component in the direction perpendicular to TDI The sum of the squares of the relative residuals between is used as the objective function H(G ₁ ^⊥ ) in the direction perpendicular to the TDI, as shown in formula (2): In formula (2), g ^⊥ (i) represents the flutter in the direction perpendicular to TDI at time i×T _r , Indicates the low-frequency component of flutter below 1 Hz in the direction perpendicular to TDI at j×T _r ; Step 3.4, use the optimization algorithm to solve the minimum points of the objective functions H(G ₁ ^|| ) and H(G ₁ ^⊥ ) respectively and and respectively as the estimated values of the flutter sub-block G ₁ ^|| along the TDI direction and the flutter sub-block G ₁ ^⊥ perpendicular to the TDI direction, as shown in formulas (3) and (4): Step 3.5. According to the push-broom imaging principle of the space camera and the imaging characteristics of the TDICCD stitching area, the flutter detection models along the TDI direction and perpendicular to the TDI direction are established as shown in formulas (5) and (6): G ₂ ^|| ＝A·G ₁ ^|| +B·S _∑ ^|| (5) G ₂ ^⊥ ＝A·G ₁ ^⊥ +B·S _∑ ^⊥ (6) In formula (5) and formula (6), A and B both represent the coefficient matrix, and Step 3.6. According to the flutter detection model along the TDI direction shown in formula (5), use formula (7) to calculate the estimated value of the flutter sub-block G ₂ ^|| along the TDI direction According to the flutter detection model perpendicular to the TDI direction shown in formula (6), use formula (8) to calculate the estimated value of the flutter sub-block G ₂ ^⊥ perpendicular to the TDI direction Step 3.7, according to the definition of flutter along the TDI direction, use formula (9) to synthesize the estimated value of the flutter G _∑ ^|| along the TDI direction According to the definition of flutter in the direction perpendicular to TDI, the estimated value of flutter G _∑ ^⊥ in the direction perpendicular to TDI is synthesized by formula (10)