Satellite stereo image flutter detection and elimination method

Abstract

The invention discloses a satellite stereo image flutter detection and elimination method comprising the following steps: S1, carrying out original image block adjustment and DSM extraction; s2, calculating a column direction flutter error sequence of a non-NAD image based on a stereo image; s3, carrying out flutter saliency detection on the non-NAD image; s4, calculating a row and column flutter error sequence of the non-NAD image based on the same-view-angle image; s5, calculating a row direction flutter error sequence of the non-NAD image based on the stereo image; and S6, eliminating flutter errors of the non-NAD image. The method has the advantages that the flutter value obtained by the same-view-angle image and the flutter value obtained by the three-dimensional image are comprehensively utilized, and image flutter error elimination is achieved.

Description

Technical Field

[0001] The invention relates to the field of photogrammetry and remote sensing technology, and in particular to a satellite stereo image chatter detection and elimination method. Background Art

[0002] Optical stereo mapping satellites use a satellite platform to carry more than two cameras with different installation angles or a single camera that swings with the platform to achieve imaging of the target area at more than two angles. The three-dimensional geometric model of the object can be restored through stereo images and orientation parameters determined by orbits, attitude parameters and camera geometric parameters. Among them, attitude accuracy is the key factor that determines the accuracy of stereo image geometric inversion. Attitude accuracy is usually evaluated by attitude stability, which is a statistical value of attitude changes over time and can be regarded as the superposition of sinusoidal curves of different frequencies, including low-frequency and high-frequency parts. The high-frequency part generally affects the geometric accuracy of stereo images.

[0003] The orientation parameters of satellite images are generally converted into a general rational function model through a rigorous geometric model. The rational function model is a fitting model that can only fit the low-frequency components in the attitude change. The high-frequency part is called jitter. If it cannot be eliminated, it will directly affect the quality of the geometric inversion of the stereo image. The present invention calls the attitude component that the rational function model cannot fit the attitude jitter residual. The attitude jitter residual can be converted into the image surface coordinate residual and decomposed into the row and column directions of the image. Among them, the coordinate residual in the column direction mainly corresponds to the upper and lower parallax of the quasi-epiconular image, which affects the quality of dense matching and even causes matching failure. The row direction residual directly corresponds to the error of elevation inversion. Therefore, accurately detecting and eliminating attitude jitter residuals is an important guarantee for the accuracy of geometric inversion. Summary of the invention

[0004] The object of the present invention is to provide a satellite stereoscopic image vibration detection and elimination method, so as to solve the above-mentioned problems existing in the prior art.

[0005] In order to achieve the above object, the technical solution adopted by the present invention is as follows: A satellite stereo image vibration detection and elimination method comprises the following steps: S1. Original image block adjustment and DSM extraction: Based on the reference base map, reference DEM, the original stereo image to be processed and the RPC file, match the image points with the same name in the NAD image and the reference base map and calculate the corresponding object plane coordinates; obtain the elevation according to the object plane coordinates of the image points with the same name and the reference DEM, and use the coordinates of the image points with the same name in each view to form control points together with the corresponding object coordinates; match the connection points between the NAD image and the images of other viewpoints, and perform regional block adjustment iteratively based on the connection points and control points to correct the RPC parameters of each scene image; obtain DSM based on the corrected RPC parameters and stereo images; S2. Calculation of column-wise dither error sequence based on stereo images for non-NAD images: The projection image of the non-NAD image is obtained, and the column coordinate difference between the non-NAD image and its projection image is calculated by matching the image points of the same name point by point; the column direction flutter image plane error of each row of the image is obtained based on the column coordinate difference, and the column direction flutter error sequence of the non-NAD image based on the stereo image is obtained based on the flutter image plane error of each row; S3. Chatter significance detection of non-NAD images: Performing a segmented Fourier transform on the column-wise jitter error sequence of the non-NAD image based on the stereo image to determine whether the non-NAD image has obvious jitter, thereby realizing jitter significance detection of the non-NAD image; S4. Calculation of row and column flutter error sequence of non-NAD images based on images with the same viewing angle: By matching the same-name image points between different viewing angles of non-jitter image pairs and between non-jitter image pairs and jitter image pairs, control points are obtained based on overlapping image pairs, reference base map and reference DEM, and regional network adjustment is performed based on non-jitter image pairs and jitter image pairs to extract the overlapping area DEM; the projection image is obtained based on the overlapping area DEM; the column coordinate difference and row coordinate difference between the non-NAD image and the projection image are calculated based on the same-name image points of the overlapping range, and the row and column direction jitter error sequence in the overlapping range is obtained; the row and column jitter error sequence of the non-NAD image based on the same-viewing angle image is obtained by resetting the jitter value; S5. Calculation of row-wise dither error sequence based on stereo images for non-NAD images: The non-NAD image is projected into the error DSM obtained by using the grid point elevation of the DSM and the reference DEM elevation to obtain the error DSM value; the elevation error sequence is obtained by using the error DSM, and the row-direction dither error sequence of the non-NAD image based on the stereo image is obtained by calculating the corresponding relationship between the elevation error and the pixel deviation; the row-column direction dither error sequence of the non-NAD image based on the stereo image is divided into the low-frequency dither component and the high-frequency dither component in the corresponding row and column directions respectively; S6, Non-NAD image jitter error elimination: Based on the anti-jitter calculation, a non-NAD image with the jitter eliminated, i.e., an anti-distorted image, is obtained; the anti-distorted coordinates of each image point of the anti-distorted image are calculated under the condition of considering or not considering the low-frequency jitter component and / or the high-frequency jitter component in step S4 and / or step S5; based on the anti-distorted coordinates, the interpolated grayscales of all image points of the anti-distorted image are obtained, and then a new anti-distorted image, i.e., an anti-distorted image with the jitter factor eliminated, is obtained.

[0006] Preferably, step S1 specifically includes the following contents: S11, using the reference base map, reference DEM, the original image to be processed and the RPC file as input data, matching the NAD image with the reference base map at evenly distributed image points of the same name (ri, ci)-(Ri, Ci); calculating the corresponding object plane coordinates (Li, Bi) from the row and column coordinates of the reference base map and the TFW file; Among them, the original image with the smallest observation angle in the original stereo image is called NAD image, and the images with other angles are non-NAD images, which are represented by BWD or FWD, BWD represents backsight, and FWD represents foresight; r and c represent the row and column coordinates of the original image, R and C represent the row and column coordinates of the reference base map, i=0, 1,…, N-1, representing the point number; Li and Bi represent longitude and latitude, respectively; S12, bilinearly interpolating the object plane coordinates of the same-name image points on the reference DEM to obtain the elevation Hi; matching the original image plane coordinates (ri, ci) of the same-name image points with the non-NAD image to obtain the same-name image points on the non-NAD image; the image plane coordinates of the same-name image points at each viewing angle corresponding to the object plane coordinates (Li, Bi, Hi) together constitute a control point, thereby obtaining N control points; S12, individually matching the evenly distributed connection points between the NAD image and the other view images, performing block adjustment on these connection points and the previously matched N control points, and modifying the RPC parameters of each view image by iterative solution; S13. The corrected RPC parameters and stereo images of each scene are used as basic data, and dense point clouds are obtained by dense matching and forward intersection of image points with the same name, and the point clouds are gridded to obtain DSM.

[0007] Preferably, in step S12, each iteration is performed based on the image-space error of the connection points and the RANSAC algorithm is used to eliminate the same-name image points whose errors exceed a certain threshold. At the same time, based on the deviation error between the front intersection coordinates of the control point and the object-space coordinates of the control point, the RANSAC algorithm is used to eliminate the control points whose errors exceed a certain threshold. Each iteration gradually reduces the threshold compared with the previous one until the threshold is less than the minimum threshold value. When the threshold is less than the minimum threshold value, the threshold is set to the minimum threshold value.

[0008] Preferably, step S2 specifically includes the following contents: S21, projecting each image point of the non-NAD perspective image onto the NAD image point by point through the calibrated RPC via the DSM, obtaining the corresponding grayscale value through bilinear interpolation, and obtaining a projection image with the same length and width as the non-NAD image; if it is a two-view stereo image, only one projection image is obtained; if it is a three-view stereo image, two projection images are obtained; S22. If there is no attitude jitter residual in the image, the column coordinate deviation between the non-NAD image and the other projected images will appear as a random quantity of an order of magnitude smaller; if there is an attitude jitter residual, the column coordinate deviation between the non-NAD image and the other projected images will appear as a systematic quantity corresponding to the jitter amplitude; match the image points with the same name of the non-NAD image and the other projected images point by point, and calculate the column coordinate difference between the two; sort all the column coordinate differences of each row of images and take the middle value as the column direction jitter image error of the row of images; calculate the jitter image error of each row in turn, and obtain the column direction jitter error sequence of the non-NAD image based on the stereo image.

[0009] Preferably, step S3 specifically comprises: Perform a piecewise Fourier transform on the column-direction dither error sequence based on the stereo image of the non-NAD image, where the segment width WF is an integer power of 2; take 1 as the base, and multiply by 2 continuously, and obtain a number that is twice the original number each time, until the result is WF; if the number of rows of the non-NAD image is an integer FN times the segment width, then only the number of segment width × FN in the column-direction dither error sequence of the non-NAD image is taken for FN times Fourier transform, and the result of each Fourier transform is represented by FD, and each element in FD is a complex number; Set a floating-point array FP with a length of WF and all elements initialized to 0. After each Fourier transform, add the modulus of the corresponding complex number in FD to each element of FP, that is, the square of the real part of the complex number plus the square of the imaginary part, and then take the square root of the sum. Calculate the sum of the 1st to 64th elements of FP as FA. The i-th s To the i e The sum of all elements between is FJ, where i s Greater than 1, i e Less than 64, corresponding to the frequency band of the attitude chatter curve within a certain period of time or under other specific conditions; calculate the ratio ρ=FJ / FA, if ρ is greater than a certain threshold, it is considered that the image has obvious chatter.

[0010] Preferably, step S4 specifically includes the following contents: S41, the same-name image points are matched between different viewing angle images of the non-jitter image pair, the same-name image points are matched between the non-jitter image pair and the jitter image pair simultaneously, all overlapping image pairs are matched control points by the method for step S1 with reference to base map and reference DEM, the non-jitter image pair and the jitter image pair are carried out block adjustment together, and then all non-jitter image pairs are utilized to extract overlapping area DSM; Among them, the NAD image and the non-NAD image acquired on the same track with an overlap range of more than 80% constitute an image pair; the image pair that has undergone the jitter significance detection in step S3 is determined to be an image pair with significant jitter, which is called a jitter image pair; and the image pair that has more than one image pair with a certain degree of overlap with the jitter image pair but does not have significant jitter is called a non-jitter image pair; S42, assuming that a non-NAD image of a tremor image pair has obvious tremor, and the M images with the same viewing angle that overlap with it do not have obvious tremor; first determine the overlapping range of the non-NAD image and a certain overlapping image, project each image point of the non-NAD image within the overlapping range to the overlapping image point by point through the overlapping area DSM, and calculate the corresponding grayscale value by bilinear interpolation, so as to obtain a projected image with the same length and width as the overlapping range of the non-NAD image; S43, matching the image points with the same name of the non-NAD image and the overlapping image in the overlapping range point by point, calculating the column coordinate difference between the two according to the coordinates of the image points with the same name, sorting all the column coordinate differences of each row of images and taking the middle value as the column direction chatter error of the row of images, calculating the chatter error of each row in turn, and obtaining a column direction chatter error sequence in the overlapping range; S44, calculating the row coordinate difference between the two according to the coordinates of the image points with the same name, sorting all the row coordinate differences of each row of the image and taking the middle value as the row direction jitter error of the row of the image, calculating the jitter error of each row in turn, and obtaining a row direction jitter error sequence within an overlapping range; S45, using the same method to calculate the overlapping range of the non-NAD image and all other overlapping images, calculate the projection image within the overlapping range, and then calculate the chatter error sequence in the column and row directions within the overlapping range; S46, for an image row with multiple jitter values, taking an average as the jitter value of the image row, and for an influence row without overlapping images, setting its jitter value to 0, thereby obtaining a row-wise jitter error sequence based on the same-viewing angle image and a column-wise jitter error sequence based on the same-viewing angle image that are consistent with the number of rows of the non-NAD image; S47, using mean filtering to separate the row and column direction chatter error sequence of the non-NAD image based on the same viewing angle image into low-frequency chatter components and high-frequency chatter components in the corresponding row and column directions.

[0011] Preferably, step S5 specifically includes the following contents: S51, projecting each grid point of the DSM onto the reference DEM according to the coordinates to obtain the corresponding floating point row and column coordinates, performing bilinear interpolation to obtain the corresponding elevation of the reference DEM, subtracting the elevation of the DSM grid point from the elevation of the reference DEM to obtain the error value of the grid point, and calculating the elevation values ​​of all grid points to obtain an error DSM; S52, projecting each image point of the non-NAD image onto the DSM through the RPC corrected by the regional block adjustment, and obtaining the error DSM value according to the bilinear interpolation of the projection point coordinates; after calculating the corresponding error DSM value for each non-NAD image, obtaining an elevation error map with the same length and width as the non-NAD image; sorting the elevation errors of each row of non-NAD images and taking the median value to obtain the average error of the row of images; after calculating each row, obtaining an elevation error sequence with a length consistent with the number of rows of the non-NAD image; S53, calculating the correspondence between the elevation error and the pixel deviation, and obtaining a chatter error sequence of the non-NAD image measured by the image plane distance based on the correspondence, that is, a row-direction chatter error sequence of the non-NAD image based on the stereo image; S54, using mean filtering to separate the column-direction jitter error sequence based on the stereo image of the non-NAD image obtained in step S2 and the row-direction jitter error sequence based on the stereo image of the non-NAD image obtained in step S53 into low-frequency jitter components and high-frequency jitter components in corresponding row and column directions.

[0012] Preferably, step S53 is specifically as follows: assuming that the HEIGHT_OFF parameter of the NAD image correction RPC is Ho, and the NAD center row and column are rc and cc; calculating the corresponding non-NAD image row and column coordinates rf and cf according to rc, cc, Ho, the NAD image correction RPC, and the non-NAD image correction RPC; calculating the corresponding non-NAD image row and column coordinates rf1 and cf1 according to rc, cc, Ho+1, the NAD image correction RPC, and the non-NAD image correction RPC in the same way; calculating the image plane distance Df between the non-NAD image image plane coordinates (rf, cf) and the coordinates (rf1, cf1); multiplying each element of the elevation error sequence by Df to obtain the jitter error sequence of the non-NAD image measured by the image plane distance, that is, the row direction jitter error sequence of the non-NAD image based on the stereo image.

[0013] Preferably, step S6 specifically includes the following contents: S61, performing anti-chatter calculation point by point to obtain a non-NAD image with eliminated chatter, that is, an anti-distorted image; calculating an anti-distorted coordinate (r1, c1) for an image point in the rth row and cth column of the anti-distorted image; S62, performing bilinear interpolation or bicubic convolution interpolation on the dewarped image according to the dewarped coordinates (r1, c1), calculating a gray value and assigning it to the point (r, c) of the dewarped image; S63, after calculating the interpolation grayscales of all the pixels of the dewarped image, a new dewarped image is obtained, which is the dewarped image with the chatter factor eliminated.

[0014] Preferably, in step S61, Ignoring the flutter value calculated in step S4, there are three ways to calculate the anti-distortion coordinates: (1) Considering both low-frequency flutter and high-frequency flutter, r1=r-JFr r , c1=r-JFc c ; (2) Considering only low-frequency flutter, r1=r-JFLr r , c1=r-JFLc c ; (3) Consider only high-frequency flutter, then r1=r-JFHr r , c1=r-JFHc c ; Among them, JFr r and JFc c are the row-wise dither sequence of the rth row and the column-wise dither sequence of the cth column of the non-NAD image based on the stereo image; JFLr r and JFLc c are the row-wise low-frequency dither component of the rth row and the column-wise low-frequency dither component of the cth column of the non-NAD image based on the stereo image; JFHr r and JFHc c They are respectively a row-direction high-frequency dither component of the r-th row and a column-direction high-frequency dither component of the c-th column of the non-NAD image based on the stereo image; Considering the flutter value calculated in step S4, there are three ways to calculate the anti-distortion coordinates: (1) Considering both the low-frequency flutter and the high-frequency flutter in step S4, r1=r-JFOr r , c1=r-JFOc c ; (2) Considering the low-frequency vibration in step S5 and the high-frequency vibration in step S4 at the same time, we can obtain: r1=r-(JFLr r +JFOHr r ), c1=r-(JFLc c +JFOHc c ); (3) Considering only the high-frequency vibration of step S4, r1=r-JFOHr r , c1=r-JFOHc c ; Among them, JFOr r and JFOc c are the row-wise dither sequence of the rth row and the column-wise dither sequence of the cth column of the non-NAD image based on the same-viewing angle image; JFOHr r and JFOHc cThey are respectively the high-frequency dither component of the rth row and the column-wise high-frequency dither component of the cth column of the non-NAD image based on the same-viewing angle image.

[0015] The beneficial effects of the present invention are: 1. The present invention eliminates the systematic error between the DSM obtained by stereo images and the reference DEM by matching dense control points, connecting points and performing regional network adjustment. 2. The present invention calculates the projection image based on the initial DSM corresponding to the stereo image, matches the deviation of each image point between the original non-NAD image and its projection image, and calculates the median of the deviation by row, which can reduce the influence of matching errors and other factors on the statistics of the jitter value. 3. The present invention quickly determines the image pairs with significant jitter through the spectrum analysis of the attitude residual in the column direction. 4. The present invention performs median filtering by row in the image side according to the elevation error between the DSM and the DEM, thereby reducing the interference of the DSM error factor on the statistics of the jitter value in the row direction. 5. The present invention comprehensively utilizes the jitter value obtained from the same-viewing angle image and the jitter value obtained from the stereo image to eliminate the image jitter error. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 is a flow chart of a method in an embodiment of the present invention. DETAILED DESCRIPTION

[0017] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings. It should be understood that the specific implementation methods described herein are only used to explain the present invention and are not used to limit the present invention.

[0018] like Figure 1 As shown, in this embodiment, a satellite stereo image vibration detection and elimination method is provided, and the method includes the following six parts: 1. Original image block adjustment and DSM extraction Based on the reference base map, reference DEM, the original stereo image to be processed and the RPC file, match the image points with the same name in the NAD image and the reference base map and calculate the corresponding object plane coordinates; obtain the elevation according to the object plane coordinates of the image points with the same name and the reference DEM, and use the coordinates of the image points with the same name in each perspective and the corresponding object coordinates to form control points; match the connection points between the NAD image and the images of other perspectives, and perform regional block adjustment iteratively based on the connection points and control points to correct the RPC parameters of each scene image; obtain DSM based on the corrected RPC parameters and the stereo image.

[0019] Specifically include the following contents: 1.1. Using the reference base map, reference DEM, the original image to be processed and the RPC file as input data, the original image with the smallest observation angle in the original stereo image (hereinafter referred to as NAD image, the images of other angles are represented by BWD or FWD, BWD represents rear view, FWD represents front view, taking the two-view stereo image of Gaofen-7 as an example, Gaofen-7 stereo image includes NAD image and FWD image, while the three-view stereo image of satellites such as Ziyuan-3 includes NAD image, BWD image and FWD image) is matched with the reference base map to match the evenly distributed image points of the same name ((ri, ci)-(Ri, Ci), r and c represent the row and column coordinates of the original image, respectively, R and C represent the row and column coordinates of the reference base map, i=0, 1,…, N-1, represents the point number), and the corresponding object plane coordinates (i.e. longitude Li and latitude Bi) can be calculated from R and C and the TFW file of the reference base map.

[0020] 1.2. Obtain the elevation Hi by bilinear interpolation on the reference DEM based on the object plane coordinates (longitude Li and latitude Bi) of the same-name image points. The original image plane coordinates (ri, ci) of the same-name image points are matched with the non-NAD image to obtain the same-name image points (rFWDi, cFWDi) on the FWD image and the same-name image points (rBWDi, cBWDi) on the BWD image (only three-view stereo images have (rBWDi, cBWDi)). The image plane coordinates of the same-name image points of each view and the corresponding object plane coordinates (Li, Bi, Hi) together constitute a control point, thus obtaining N control points.

[0021] 1.3. Separately match the evenly distributed connection points between the NAD image and the other view images, perform regional block adjustment on these connection points and the previously matched N control points, and modify the RPC parameters of each scene image through iterative solution. Each iteration, according to the image error of the connection point, the RANSAC algorithm is used to remove the same-name image points whose errors exceed a certain threshold. At the same time, according to the deviation error between the intersection coordinates in front of the control point and the object coordinates of the control point, the RANSAC algorithm is used to remove the control points whose errors exceed a certain threshold. Each iteration gradually reduces the threshold compared with the previous one until the threshold is less than the minimum threshold value. When the threshold is less than the minimum threshold value, the threshold is set to the minimum threshold value.

[0022] 1.4. Take the corrected RPC and stereo images of each scene as basic data, obtain dense point cloud through dense matching and forward intersection of image points with the same name, and grid the point cloud to obtain DSM.

[0023] 2. Column-wise Flutter Error Sequence Calculation Based on Stereo Images for Non-NAD Images The projection image of the non-NAD image is obtained, and the column coordinate difference between the non-NAD image and its projection image is calculated by matching the image points of the same name point by point; the column direction flutter image plane error of each row of the image is obtained based on the column coordinate difference, and the column direction flutter error sequence of the non-NAD image based on the stereo image is obtained based on the flutter image plane error of each row.

[0024] Specifically include the following contents: 2.1. Project each pixel of the non-NAD perspective image (such as FWD image or BWD image) to the NAD image point by point through the calibrated RPC and DSM, and obtain the corresponding grayscale value through bilinear interpolation to obtain a projection image with the same length and width as the non-NAD image. If it is a two-view stereo image, only one projection image is obtained (such as the projection image corresponding to the FWD image). If it is a three-view stereo image, two projection images are obtained (the projection image corresponding to the FWD image and the projection image corresponding to the BWD image).

[0025] 2.2. Taking the FWD image and its projected image as an example, if there is no attitude jitter residual in the image, the column coordinate deviation between the FWD image and its projected image will appear as a random quantity with a very small order of magnitude. If there is an attitude jitter residual, the column coordinate deviation between the FWD image and its projected image will appear as a systematic quantity corresponding to the jitter amplitude. Match the same-name image points of the FWD image and its projected image point by point, and calculate the column coordinate difference between the two. Sort all the column coordinate differences of each row of the image and take the middle value as the column direction jitter image error of the row of the image. Calculate the jitter image error of each row in turn to obtain a column direction jitter error sequence JFc i (i=0,1,…,HF-1), i.e., the column-wise dither error sequence of the FWD image based on the stereo image, and HF is the number of FWD image rows.

[0026] 2.3. If there are non-NAD images from other perspectives, the same method is used to calculate the column-wise dither error sequence based on the stereo image for the remaining non-NAD perspective images, such as the column-wise dither error sequence JBc based on the stereo image for the BWD image. i (i=0,1,…,HB-1), HB is the number of BWD image rows.

[0027] 3. Chatter Significance Detection of Non-NAD Images The column-wise jitter error sequence of the non-NAD image based on the stereo image is subjected to piecewise Fourier transform to determine whether the non-NAD image has obvious jitter, thereby realizing jitter significance detection of the non-NAD image.

[0028] Specifically, it includes the following contents: 3.1. Taking FWD image as an example, JFc i(i=0,1,…,HF-1) performs a piecewise Fourier transform, with the piecewise width being an integer power of 2, denoted by WF, and assuming WF is 1024. Using 1 as the base, multiply by 2 continuously, each time getting a number that is twice the original number, until the result is WF, and record the number of consecutive multiplications as w. Assuming HF is an integer FN times WF, then only JFc is taken i The number of WF×FN in (i=0,1,…,HF-1) is subjected to FN Fourier transforms, and the result of each Fourier transform is expressed as FD i (i=0,1,…,WF-1) means that each element in FD is a complex number. Set a floating point array FP with a length of WF and all elements initialized to 0 i (i=0,1,…,WF-1), each time the Fourier transform is completed, FP i Add FD to each element of (i=0,1,…,WF-1) i (i=0,1,…,WF-1) corresponds to the modulus of the complex number, that is, the square of the real part plus the square of the imaginary part of the complex number, and then the square root of the sum. Calculate FP i (i=0,1,…,WF-1) The sum of the 1st to 64th elements is FA, and the i-th s To the i e The sum of all elements between is FJ. s Greater than 1, i e Less than 64, such as i s is equal to 17, i e Equal to 21, which is an empirical value, corresponding to the frequency band of the attitude chatter curve in a certain period of time or under other specific conditions. Calculate the ratio ρ=FJ / FA. If ρ is greater than a certain threshold, such as ρ>0.3, it is considered that the image has obvious chatter.

[0029] 3.2. The same method is used to detect the chatter significance of other non-NAD images.

[0030] 4. Calculation of row-column flutter error sequence based on same-viewing angle images for non-NAD images By matching the same-name image points between different viewing angles of non-jitter image pairs and between non-jitter image pairs and jitter image pairs, control points are obtained based on overlapping image pairs, reference base map and reference DEM, and regional block adjustment is performed based on the non-jitter image pairs and jitter image pairs to extract the overlapping area DEM; the projection image is obtained based on the overlapping area DEM; the column coordinate difference and row coordinate difference between the non-NAD image and the projection image are calculated based on the same-name image points in the overlapping range to obtain the row and column direction jitter error sequence within the overlapping range; the row and column jitter error sequence of the non-NAD image based on the same-viewing angle image is obtained by resetting the jitter value.

[0031] Specifically, it includes the following contents: 4.1. Assume that an image pair (NAD image and non-NAD image acquired on the same track with an overlap range of more than 80% constitute an image pair) is determined as an image pair with significant jitter (referred to as jitter image pair) after the jitter significance detection in step S3. At the same time, there is one or more image pairs with a certain degree of overlap with this image pair (the overlap degree is calculated according to the overlap degree of NAD images of each image pair), but there is no significant jitter (hereinafter referred to as non-jitter image pairs). Match the same-name image points between the different perspective images of these non-jitter image pairs, and match the same-name image points between non-jitter image pairs and jitter image pairs. All overlapping image pairs are matched with control points by the reference base map and reference DEM using the method of step S1. Perform regional block adjustment on the non-jitter image pairs and the jitter image pairs, and then use all the non-jitter image pairs to extract the overlapping area DSM.

[0032] 4.2. Assume that the FWD image of a tremor image pair has obvious tremor, and the M overlapping images of the same viewing angle do not have obvious tremor. First, determine the FWD image and the i-th (i=0,1,…,M-1) overlapping image (called FWD i The overlapping range of the FWD image) is projected point by point through the overlapping area DSM to the FWD i The corresponding grayscale value of the image is calculated by bilinear interpolation, thereby obtaining a projection image with the same length and width as the overlapping range of the FWD image.

[0033] 4.3. Point-by-point matching of overlapping range FWD images and FWD i For the image points with the same name in the projected image, the column coordinate difference between the two is calculated based on the coordinates of the image points with the same name. After sorting all the column coordinate differences of each row of the image, the middle value is taken as the column direction flutter error of the row of the image. The flutter error of each row is calculated in turn to obtain a column direction flutter error sequence within the overlapping range.

[0034] 4.4. Calculate the row coordinate difference between the two based on the coordinates of the same-named image points. Sort all the row coordinate differences of each row of the image and take the middle value as the row direction jitter error of the row of the image. Calculate the jitter error of each row in turn to obtain a row direction jitter error sequence within the overlapping range.

[0035] 4.5. The same method is used to calculate the overlapping range of the FWD image and the 1st to M-1th overlapping images of the same viewing angle, calculate the projection image within the overlapping range, and then calculate the vibration error sequence in the column and row directions within the overlapping range.

[0036] 4.6. After calculating the jitter error sequences in both the column and row directions corresponding to the overlapping images, some rows may have multiple jitter values, and the average is taken as the jitter value of the row image. Some image rows do not have overlapping images, and their jitter values ​​are set to 0. In this way, the row direction jitter error sequence JFOr based on the same viewing angle image that is consistent with the number of FWD image rows is obtained. i (i=0,1,…,HF-1), and the column-wise flutter error sequence JFOc based on the same-viewing angle image i (i=0,1,…,HF-1), HF is the number of FWD image rows.

[0037] 4.7. The same method is used to calculate the jitter error sequence of other non-NAD images based on the same viewing angle, such as the jitter error sequence of the BWD image in both row and column directions. The jitter error sequence in the BWD row direction is JBOr i (i=0,1,…,HB-1), the column-wise chatter error sequence is JBOc i (i=0,1,…,HB-1), HB is the number of BWD image rows.

[0038] 4.8. Mean filtering divides the chatter into low-frequency and high-frequency components: Taking FWD images as an example, the JFOc in the column direction of the same-viewing angle image is i The (i=0,1,…,HF-1) sequence is mean filtered with the filter width set to 200 to obtain the column-wise low-frequency flutter component JFOLc i (i=0,1,…,HF-1), change JFOc i Subtract JFOLc from each element of the sequence (i=0,1,…,HF-1) i (i=0,1,…,HF-1) corresponding elements, the column-wise high-frequency flutter component JFOHc is obtained i (i=0,1,…,HF-1).

[0039] 4.9. The same method is used to obtain the row-wise low-frequency jitter component JFOLr based on the same-viewing angle image. i (i=0,1,…,HF-1) and high frequency flutter component JFOHr i (i=0,1,…,HF-1).

[0040] 5. Calculation of row-wise flutter error sequence based on stereo images for non-NAD images The non-NAD image is projected into the error DSM obtained by using the grid point elevation of DSM and the reference DEM elevation to obtain the error DSM value; the elevation error sequence is obtained by using the error DSM, and the row direction dither error sequence of the non-NAD image based on the stereo image is obtained by calculating the corresponding relationship between the elevation error and the pixel deviation; the row and column direction dither error sequence of the non-NAD image based on the stereo image is divided into the low-frequency dither component and the high-frequency dither component in the corresponding row and column directions.

[0041] 5.1. Project each grid point of DSM to the reference DEM to obtain the corresponding floating point row and column coordinates. Perform bilinear interpolation to obtain the corresponding elevation of the reference DEM. Subtract the elevation of the DSM grid point from the elevation of the reference DEM to obtain the error value of the grid point. Calculate the elevation values ​​of all grid points to obtain an error DSM.

[0042] 5.2. Project each image point of the FWD image onto the DSM through the RPC corrected by the regional block adjustment, and obtain the error DSM value based on the bilinear interpolation of the projection point coordinates. After calculating the corresponding error DSM value for each FWD image, an elevation error map with the same length and width as the FWD image is obtained. Sort the elevation errors of each row of FWD and take the median to obtain the average error of the row of images. After each row is calculated, an elevation error sequence JFHr with a length consistent with the number of FWD rows is obtained. i (i=0,1,…,HF-1), HF is the number of FWD image rows.

[0043] 5.3. Calculate the correspondence between elevation error and pixel deviation: Assume that the HEIGHT_OFF parameter of NAD image correction RPC is Ho, and the NAD center row and column are rc and cc. Calculate the corresponding FWD image row and column coordinates rf and cf based on rc, cc and Ho and NAD image correction RPC and FWD image correction RPC. Similarly, calculate the corresponding FWD image row and column coordinates rf1 and cf1 based on rc, cc and Ho+1 and NAD image correction RPC and FWD image correction RPC. Calculate the image plane distance Df between the FWD image plane coordinates (rf, cf) and (rf1, cf1).

[0044] 5.4、JFHR i Each element of (i=0,1,…,HF-1) is multiplied by Df to obtain the flutter error sequence JFr measured by the image plane distance i (i=0,1,…,HF-1), which is the row-wise dither error sequence based on stereo images.

[0045] 5.5. The same method is used to calculate the dither error sequence based on stereo images of other non-NAD images, such as the dither error sequence based on stereo images of BWD images JBr i(i=0,1,…,HF-1).

[0046] 5.6. Mean filtering divides the chatter into low-frequency components and high-frequency components: Taking FWD images as an example, the JFc in the column direction of the stereo image is i The (i=0,1,…,HF-1) sequence is mean filtered, and the filter width is set to 200 to obtain the column-wise low-frequency vibration component JFLc i (i=0,1,…,HF-1), JFc i Subtract JFLc from each element of the sequence (i=0,1,…,HF-1) i (i=0,1,…,HF-1) corresponding elements, the high-frequency dither component JFHc in the column direction of the FWD image is obtained i (i=0,1,…,HF-1).

[0047] 5.7. The same method is used to obtain the low-frequency jitter component JFLr in the row direction of the FWD image based on the stereo image i (i=0,1,…,HF-1) and high frequency flutter component JFHr i (i=0,1,…,HF-1).

[0048] 6. Non-NAD Image Jitter Error Elimination Based on the anti-jitter calculation, a non-NAD image with the jitter eliminated, i.e., an anti-distorted image, is obtained; the anti-distorted coordinates of each image point of the anti-distorted image are calculated under the condition of considering or not considering the low-frequency jitter component and / or the high-frequency jitter component in step S4 and / or step S5; based on the anti-distorted coordinates, the interpolated grayscales of all image points of the anti-distorted image are obtained, and then a new anti-distorted image, i.e., an anti-distorted image with the jitter factor eliminated, is obtained.

[0049] Specifically, it includes the following contents: 6.1. Perform anti-chatter calculation point by point to obtain a non-NAD image with chatter eliminated (hereinafter referred to as an anti-distorted image); taking the FWD image as an example, calculate the anti-distorted coordinates (r1, c1) for the image point in the rth row and cth column of the anti-distorted image. There are two cases for the anti-distorted coordinates: A. Ignoring the flutter value calculated in the fourth part, there are three calculation methods: (1) Considering both low-frequency flutter and high-frequency flutter, r1=r-JFr r , c1=r-JFc c .

[0050] (2) Considering only low-frequency flutter, r1=r-JFLr r , c1=r-JFLc c .

[0051] (3) Consider only high-frequency flutter, then r1=r-JFHr r , c1=r-JFHc c .

[0052] B. Considering the flutter value calculated in the fourth part, there are three calculation methods: (1) Considering both the low-frequency flutter and the high-frequency flutter in step S4, r1=r-JFOr r , c1=r-JFOc c .

[0053] (2) Considering the low-frequency vibration in step S5 and the high-frequency vibration in step S4 at the same time, we can obtain: r1=r-(JFLr r +JFOHr r ), c1=r-(JFLc c +JFOHc c ).

[0054] (3) Considering only the high-frequency vibration of step S4, r1=r-JFOHr r , c1=r-JFOHc c .

[0055] 6.2. Perform bilinear interpolation or bicubic convolution interpolation on the FWD image according to the dewarped coordinates (r1, c1), calculate the grayscale value and assign it to the dewarped FWD image point (r, c).

[0056] 6.3. After calculating the interpolation grayscale of all pixels in the FWD anti-distortion image, a new FWD image is obtained, which is the anti-distortion image with the chatter factor eliminated.

[0057] 6.4. The same method is used to calculate the dewarped images of other non-NAD viewing angles.

[0058] By adopting the above technical solution disclosed in the present invention, the following beneficial effects are obtained: The present invention provides a satellite stereo image jitter detection and elimination method, which eliminates the systematic error between the DSM obtained by the stereo image and the reference DEM by matching dense control points, connecting points and performing regional network adjustment. The present invention calculates the projection image based on the initial DSM corresponding to the stereo image, matches the deviation of each image point between the original non-NAD image and its projection image, and calculates the median of the deviation by row, which can reduce the influence of factors such as matching error on the statistics of the jitter value. The present invention quickly determines the image pairs with significant jitter through the spectrum analysis of the attitude residual in the column direction. The present invention performs median filtering by row in the image plane according to the elevation error between the DSM and the DEM, and reduces the interference of the DSM error factor on the statistics of the jitter value in the row direction. The present invention comprehensively utilizes the jitter value obtained by the same-viewing angle image and the jitter value obtained by the stereo image to eliminate the image jitter error.

[0059] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be considered as the scope of protection of the present invention.

Claims

1. A satellite stereo image vibration detection and elimination method, characterized in that: The following steps are included: S1. Original image block adjustment and DSM extraction: Based on the reference base map, reference DEM, the original stereo image to be processed and the RPC file, match the image points with the same name in the NAD image and the reference base map and calculate the corresponding object plane coordinates; obtain the elevation according to the object plane coordinates of the image points with the same name and the reference DEM, and use the coordinates of the image points with the same name in each view to form control points together with the corresponding object coordinates; match the connection points between the NAD image and the images of other viewpoints, and perform regional block adjustment iteratively based on the connection points and control points to correct the RPC parameters of each scene image; obtain DSM based on the corrected RPC parameters and stereo images; S2. Calculation of column-wise dither error sequence based on stereo images for non-NAD images: The projection image of the non-NAD image is obtained, and the column coordinate difference between the non-NAD image and its projection image is calculated by matching the image points of the same name point by point; the column direction flutter image plane error of each row of the image is obtained based on the column coordinate difference, and the column direction flutter error sequence of the non-NAD image based on the stereo image is obtained based on the flutter image plane error of each row; S3. Chatter significance detection of non-NAD images: Performing a segmented Fourier transform on the column-wise jitter error sequence of the non-NAD image based on the stereo image to determine whether the non-NAD image has obvious jitter, thereby realizing jitter significance detection of the non-NAD image; S4. Calculation of row and column flutter error sequence of non-NAD images based on images with the same viewing angle: By matching the same-name image points between different viewing angles of non-jitter image pairs and between non-jitter image pairs and jitter image pairs, control points are obtained based on overlapping image pairs, reference base map and reference DEM, and regional block adjustment is performed based on non-jitter image pairs and jitter image pairs to extract the overlapping area DEM; The projection image is obtained based on the DEM of the overlapping area; the column coordinate difference and row coordinate difference between the non-NAD image and the projection image are calculated based on the same-name image points in the overlapping range to obtain the row and column direction dither error sequence within the overlapping range; the row and column dither error sequence of the non-NAD image based on the same-viewing angle image is obtained by resetting the dither value; S5. Calculation of row-wise dither error sequence based on stereo images for non-NAD images: The non-NAD image is projected into the error DSM obtained by using the grid point elevation of the DSM and the reference DEM elevation to obtain the error DSM value; the elevation error sequence is obtained by using the error DSM, and the row-direction dither error sequence of the non-NAD image based on the stereo image is obtained by calculating the corresponding relationship between the elevation error and the pixel deviation; the row-column direction dither error sequence of the non-NAD image based on the stereo image is divided into the low-frequency dither component and the high-frequency dither component in the corresponding row and column directions respectively; S6, Non-NAD Image Jitter Error Elimination: Based on the anti-jitter calculation, a non-NAD image with the jitter eliminated, i.e., an anti-distorted image, is obtained; the anti-distorted coordinates of each image point of the anti-distorted image are calculated under the condition of considering or not considering the low-frequency jitter component and / or the high-frequency jitter component in step S4 and / or step S5; based on the anti-distorted coordinates, the interpolated grayscales of all image points of the anti-distorted image are obtained, and then a new anti-distorted image, i.e., an anti-distorted image with the jitter factor eliminated, is obtained.

2. The satellite stereo image chatter detection and elimination method according to claim 1, characterized in that: Step S1 specifically includes the following contents: S11, using the reference base map, reference DEM, the original image to be processed and the RPC file as input data, matching the NAD image with the reference base map at evenly distributed image points of the same name (ri, ci)-(Ri, Ci); calculating the corresponding object plane coordinates (Li, Bi) from the row and column coordinates of the reference base map and the TFW file; Among them, the original image with the smallest observation angle in the original stereo image is called NAD image, and the images with other angles are non-NAD images, which are represented by BWD or FWD, BWD represents backsight, and FWD represents foresight; r and c represent the row and column coordinates of the original image, R and C represent the row and column coordinates of the reference base map, i=0, 1,…, N-1, representing the point number; Li and Bi represent longitude and latitude, respectively; S12, bilinearly interpolating the object plane coordinates of the same-name image points on the reference DEM to obtain the elevation Hi; matching the original image plane coordinates (ri, ci) of the same-name image points with the non-NAD image to obtain the same-name image points on the non-NAD image; the image plane coordinates of the same-name image points at each viewing angle corresponding to the object plane coordinates (Li, Bi, Hi) together constitute a control point, thereby obtaining N control points; S12, individually matching the evenly distributed connection points between the NAD image and the other view images, performing block adjustment on these connection points and the previously matched N control points, and modifying the RPC parameters of each view image by iterative solution; S13. The corrected RPC parameters and stereo images of each scene are used as basic data, and dense point clouds are obtained by dense matching and forward intersection of image points with the same name, and the point clouds are gridded to obtain DSM.

3. The satellite stereo image chatter detection and elimination method according to claim 2, characterized in that: In step S12, each iteration is performed based on the image-space error of the connection points, and the RANSAC algorithm is used to eliminate the same-name image points whose errors exceed a certain threshold. At the same time, the RANSAC algorithm is used to eliminate the control points whose errors exceed a certain threshold according to the deviation error between the front intersection coordinates of the control points and the object-space coordinates of the control points. The threshold is gradually reduced compared with the previous iteration until the threshold is less than the minimum threshold value. When the threshold is less than the minimum threshold value, the threshold is set to the minimum threshold value.

4. The satellite stereo image vibration detection and elimination method according to claim 3, characterized in that: Step S2 specifically includes the following contents: S21, projecting each image point of the non-NAD perspective image onto the NAD image point by point through the calibrated RPC via the DSM, obtaining the corresponding grayscale value through bilinear interpolation, and obtaining a projection image with the same length and width as the non-NAD image; if it is a two-view stereo image, only one projection image is obtained; if it is a three-view stereo image, two projection images are obtained; S22. If there is no attitude jitter residual in the image, the column coordinate deviation between the non-NAD image and the other projected images will appear as a random quantity of an order of magnitude smaller; if there is an attitude jitter residual, the column coordinate deviation between the non-NAD image and the other projected images will appear as a systematic quantity corresponding to the jitter amplitude; match the image points with the same name of the non-NAD image and the other projected images point by point, and calculate the column coordinate difference between the two; sort all the column coordinate differences of each row of images and take the middle value as the column direction jitter image error of the row of images; calculate the jitter image error of each row in turn, and obtain the column direction jitter error sequence of the non-NAD image based on the stereo image.

5. The satellite stereo image vibration detection and elimination method according to claim 4, characterized in that: Step S3 specifically comprises: Perform a piecewise Fourier transform on the column-direction dither error sequence based on the stereo image of the non-NAD image, where the segment width WF is an integer power of 2; take 1 as the base, and multiply by 2 continuously, and obtain a number that is twice the original number each time, until the result is WF; if the number of rows of the non-NAD image is an integer FN times the segment width, then only the number of segment width × FN in the column-direction dither error sequence of the non-NAD image is taken for FN times Fourier transform, and the result of each Fourier transform is represented by FD, and each element in FD is a complex number; Set a floating-point array FP with a length of WF and all elements initialized to 0. After each Fourier transform, add the modulus of the corresponding complex number in FD to each element of FP, that is, the square of the real part of the complex number plus the square of the imaginary part, and then take the square root of the sum. Calculate the sum of the 1st to 64th elements of FP as FA. The i-th s To the i e The sum of all elements between is FJ, where i s Greater than 1, i e Less than 64, corresponding to the frequency band of the attitude flutter curve within a certain period of time or under other specific conditions; The ratio ρ=FJ / FA is calculated. If ρ is greater than a certain threshold, it is considered that the image has obvious chatter.

6. The satellite stereo image chatter detection and elimination method according to claim 5, characterized in that: Step S4 specifically includes the following contents: S41, the same-name image points are matched between different viewing angle images of the non-jitter image pair, the same-name image points are matched between the non-jitter image pair and the jitter image pair simultaneously, all overlapping image pairs are matched control points by the method for step S1 with reference to base map and reference DEM, the non-jitter image pair and the jitter image pair are carried out block adjustment together, and then all non-jitter image pairs are utilized to extract overlapping area DSM; Among them, the NAD image and the non-NAD image acquired on the same track with an overlap range of more than 80% constitute an image pair; the image pair that has undergone the jitter significance detection in step S3 is determined to be an image pair with significant jitter, which is called a jitter image pair; and the image pair that has more than one image pair with a certain degree of overlap with the jitter image pair but does not have significant jitter is called a non-jitter image pair; S42, assuming that a non-NAD image of a tremor image pair has obvious tremor, and the M images with the same viewing angle that overlap with it do not have obvious tremor; first determine the overlapping range of the non-NAD image and a certain overlapping image, project each image point of the non-NAD image within the overlapping range to the overlapping image point by point through the overlapping area DSM, and calculate the corresponding grayscale value by bilinear interpolation, so as to obtain a projected image with the same length and width as the overlapping range of the non-NAD image; S43, matching the image points with the same name of the non-NAD image and the overlapping image in the overlapping range point by point, calculating the column coordinate difference between the two according to the coordinates of the image points with the same name, sorting all the column coordinate differences of each row of images and taking the middle value as the column direction chatter error of the row of images, calculating the chatter error of each row in turn, and obtaining a column direction chatter error sequence in the overlapping range; S44, calculating the row coordinate difference between the two according to the coordinates of the image points with the same name, sorting all the row coordinate differences of each row of the image and taking the middle value as the row direction jitter error of the row of the image, calculating the jitter error of each row in turn, and obtaining a row direction jitter error sequence within an overlapping range; S45, using the same method to calculate the overlapping range of the non-NAD image and all other overlapping images, calculate the projection image within the overlapping range, and then calculate the chatter error sequence in the column and row directions within the overlapping range; S46, for an image row with multiple jitter values, taking an average as the jitter value of the image row, and for an influence row without overlapping images, setting its jitter value to 0, thereby obtaining a row-wise jitter error sequence based on the same-viewing angle image and a column-wise jitter error sequence based on the same-viewing angle image that are consistent with the number of rows of the non-NAD image; S47, using mean filtering to separate the row and column direction chatter error sequence of the non-NAD image based on the same viewing angle image into low-frequency chatter components and high-frequency chatter components in the corresponding row and column directions.

7. The satellite stereo image chatter detection and elimination method according to claim 6, characterized in that: Step S5 specifically includes the following contents: S51, projecting each grid point of the DSM onto the reference DEM according to the coordinates to obtain the corresponding floating point row and column coordinates, performing bilinear interpolation to obtain the corresponding elevation of the reference DEM, subtracting the elevation of the DSM grid point from the elevation of the reference DEM to obtain the error value of the grid point, and calculating the elevation values ​​of all grid points to obtain an error DSM; S52, projecting each image point of the non-NAD image onto the DSM through the RPC corrected by the regional block adjustment, and obtaining the error DSM value according to the bilinear interpolation of the projection point coordinates; after calculating the corresponding error DSM value for each non-NAD image, obtaining an elevation error map with the same length and width as the non-NAD image; sorting the elevation errors of each row of non-NAD images and taking the median value to obtain the average error of the row of images; after calculating each row, obtaining an elevation error sequence with a length consistent with the number of rows of the non-NAD image; S53, calculating the correspondence between the elevation error and the pixel deviation, and obtaining a chatter error sequence of the non-NAD image measured by the image plane distance based on the correspondence, that is, a row-direction chatter error sequence of the non-NAD image based on the stereo image; S54, using mean filtering to separate the column-direction jitter error sequence based on the stereo image of the non-NAD image obtained in step S2 and the row-direction jitter error sequence based on the stereo image of the non-NAD image obtained in step S53 into low-frequency jitter components and high-frequency jitter components in corresponding row and column directions.

8. The satellite stereo image chatter detection and elimination method according to claim 7, characterized in that: Specifically, step S53 is as follows: assuming that the HEIGHT_OFF parameter of the NAD image correction RPC is Ho, and the NAD center row and column are rc and cc; calculating the corresponding non-NAD image row and column coordinates rf and cf according to rc, cc, Ho, the NAD image correction RPC, and the non-NAD image correction RPC; calculating the corresponding non-NAD image row and column coordinates rf1 and cf1 according to rc, cc, Ho+1, the NAD image correction RPC, and the non-NAD image correction RPC in the same way; calculating the image plane coordinate (rf, cf) of the non-NAD image and the image plane distance Df between the coordinates (rf1, cf1); multiplying each element of the elevation error sequence by Df to obtain the jitter error sequence of the non-NAD image measured by the image plane distance, that is, the row direction jitter error sequence of the non-NAD image based on the stereo image.

9. The satellite stereo image chatter detection and elimination method according to claim 8, characterized in that: Step S6 specifically includes the following contents: S61, performing anti-chatter calculation point by point to obtain a non-NAD image with eliminated chatter, that is, an anti-distorted image; calculating an anti-distorted coordinate (r1, c1) for an image point in the rth row and cth column of the anti-distorted image; S62, performing bilinear interpolation or bicubic convolution interpolation on the dewarped image according to the dewarped coordinates (r1, c1), calculating a gray value and assigning it to the point (r, c) of the dewarped image; S63, after calculating the interpolation grayscales of all the pixels of the dewarped image, a new dewarped image is obtained, which is the dewarped image with the chatter factor eliminated.

10. The satellite stereo image chatter detection and elimination method according to claim 9, characterized in that: In step S61, Ignoring the flutter value calculated in step S4, there are three ways to calculate the anti-distortion coordinates: (1) Considering both low-frequency flutter and high-frequency flutter, r1=r-JFr r , c1=r-JFc c ; (2) Considering only low-frequency flutter, r1=r-JFLr r , c1=r-JFLc c ; (3) Consider only high-frequency flutter, then r1=r-JFHr r , c1=r-JFHc c ; Among them, JFr r and JFc c are the row-wise dither sequence of the rth row and the column-wise dither sequence of the cth column of the non-NAD image based on the stereo image; JFLr r and JFLc c are the row-wise low-frequency dither component of the rth row and the column-wise low-frequency dither component of the cth column of the non-NAD image based on the stereo image; JFHr r and JFHc c They are respectively a row-direction high-frequency dither component of the r-th row and a column-direction high-frequency dither component of the c-th column of the non-NAD image based on the stereo image; Considering the flutter value calculated in step S4, there are three ways to calculate the anti-distortion coordinates: (1) Considering both the low-frequency flutter and the high-frequency flutter in step S4, r1=r-JFOr r , c1=r-JFOc c ; (2) Considering the low-frequency vibration in step S5 and the high-frequency vibration in step S4 at the same time, we can obtain: r1=r-(JFLr r +JFOHr r ),c1=r-(JFLc c +JFOHc c ); (3) Considering only the high-frequency vibration of step S4, r1=r-JFOHr r , c1=r-JFOHc c ; Among them, JFOr r and JFOc c are the row-wise dither sequence of the rth row and the column-wise dither sequence of the cth column of the non-NAD image based on the same-viewing angle image; JFOHr r and JFOHc c They are respectively the high-frequency dither component of the rth row and the column-wise high-frequency dither component of the cth column of the non-NAD image based on the same-viewing angle image.

Images