Satellite Stereo Image Flutter Detection and Elimination Method

Abstract

The present invention discloses a method for detecting and eliminating satellite stereo image flutter, which includes the following steps: S1, original image block adjustment and DSM extraction; S2, calculation of column-direction flutter error sequence based on stereo images of non-NAD images; S3, flutter significance detection of non-NAD images; S4, calculation of row and column flutter error sequences based on same-viewpoint images of non-NAD images; S5, calculation of row-direction flutter error sequence based on stereo images of non-NAD images; S6, elimination of flutter error of non-NAD images. The advantages are: comprehensively using the flutter values obtained from same-viewpoint images and the flutter values obtained from stereo images to achieve the elimination of image flutter error.

Description

Technical Field

[0001] The present invention relates to the technical field of photogrammetry and remote sensing, and particularly to a method for detecting and eliminating satellite stereo image flutter. Background Art

[0002] An optical stereo mapping satellite forms images of a target area from two or more angles by carrying more than two cameras with different installation angles on a satellite platform or by swinging a single camera with the platform. The three-dimensional geometric model of the object space can be restored through the stereo images, orbit, attitude parameters, and orientation parameters determined by the camera geometric parameters. Among them, the attitude accuracy is the key factor determining the geometric inversion accuracy of the stereo images. The attitude accuracy is usually evaluated by the attitude stability, which is a statistical value of the attitude changing with time and can be regarded as the superposition of sine curves with different frequencies, including a low-frequency part and a high-frequency part. Generally, the high-frequency part affects the geometric accuracy of the 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 flutter. If it cannot be eliminated, it will directly affect the quality of the geometric inversion of the stereo images. In the present invention, the attitude component that cannot be fitted by the rational function model is called the attitude flutter residual. The attitude flutter residual can be converted into image plane coordinate residuals and decomposed into the row and column directions of the image. Among them, the coordinate residuals in the column direction mainly correspond to the vertical parallax of the quasi-epipolar images, affecting the dense matching quality and even causing matching failure. The row direction residuals directly correspond to the errors in elevation inversion. Therefore, accurately detecting and eliminating the attitude flutter residuals is an important guarantee for the geometric inversion accuracy. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for detecting and eliminating satellite stereo image flutter, so as to solve the foregoing problems existing in the prior art.

[0005] To achieve the above purpose, the technical solution adopted by the present invention is as follows:

[0006] A method for detecting and eliminating satellite stereo image flutter includes the following steps:

[0007] S1. Original image block adjustment and DSM extraction:

[0008] Based on the reference base map, reference DEM, the original stereo image to be processed and the RPC file, match the homologous image points of 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 homologous image points and the reference DEM, and construct control points by combining the image plane homologous image point coordinates of each view with the corresponding object coordinates; match the connection points between the NAD image and other view images, and perform block adjustment based on the connection points and control points to iteratively correct the RPC parameters of each scene image; obtain the DSM based on the corrected RPC parameters and the stereo image;

[0009] S2. Calculation of the column direction flutter error sequence of non-NAD images based on stereo images:

[0010] Obtain the projection image of the non-NAD image, and calculate the column coordinate difference between the non-NAD image and its projection image by point-by-point matching of the homologous image points; obtain the column direction flutter image plane error of each row of the image based on the column coordinate difference, and obtain the column direction flutter error sequence of the non-NAD image based on stereo images based on the flutter image plane error of each row;

[0011] S3. Flutter significance detection of non-NAD images:

[0012] Perform piecewise Fourier transform on the column direction flutter error sequence of the non-NAD image based on stereo images to determine whether there is obvious flutter in the non-NAD image, and realize the flutter significance detection of the non-NAD image;

[0013] S4. Calculation of the row and column flutter error sequence of non-NAD images based on images of the same view:

[0014] By matching the homologous image points between different view images of non-flutter image pairs and between non-flutter image pairs and flutter image pairs, obtain control points based on overlapping image pairs, reference base map and reference DEM, perform block adjustment based on non-flutter image pairs and flutter image pairs to extract the overlapping area DEM; obtain the projection image based on the overlapping area DEM; calculate the column coordinate difference and row coordinate difference between the non-NAD image and the projection image within the overlapping range based on the homologous image points to obtain the row and column direction flutter error sequence within the overlapping range; obtain the row and column flutter error sequence of the non-NAD image based on images of the same view by resetting the flutter value;

[0015] S5. Calculation of the row direction flutter error sequence of non-NAD images based on stereo images:

[0016] By projecting non-NAD images onto the error DSM obtained from the elevation of grid points using DSM and the elevation of the reference DEM, the error DSM values are obtained; the elevation error sequence is obtained using the error DSM, and the stereo-image-based row-direction flutter error sequence of the non-NAD image is obtained by calculating the correspondence between the elevation error and the pixel deviation amount; the stereo-image-based row and column direction flutter error sequences of the non-NAD image are respectively divided into low-frequency flutter components and high-frequency flutter components in the corresponding row and column directions.

[0017] S6. Elimination of non-NAD image flutter error:

[0018] Based on anti-flutter calculation, the non-NAD image with flutter eliminated, i.e., the anti-distortion image, is obtained; the anti-distortion coordinates of each image point of the anti-distortion image are calculated under the condition of considering or not considering the low-frequency flutter component and / or high-frequency flutter component in step S4 and / or step S5; based on the anti-distortion coordinates, the interpolated gray levels of all image points of the anti-distortion image are obtained, and then a new anti-distortion image, i.e., the anti-distortion image with flutter factors eliminated, is obtained.

[0019] Preferably, step S1 specifically includes the following content.

[0020] S11. Using the reference base map, reference DEM, the original image to be processed, and the RPC file as input data, the homologous image points (ri, ci)-(Ri, Ci) evenly distributed by matching the NAD image with the reference base map are obtained; the corresponding object plane coordinates (Li, Bi) are calculated from the row and column coordinates of the reference base map and the TFW file.

[0021] Among them, the original image with the smallest observation angle in the original stereo images is called the NAD image, and the images of the remaining perspectives are non-NAD images, denoted by BWD or FWD, where BWD represents the back view and FWD represents the front view; r and c respectively 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 respectively represent longitude and latitude.

[0022] S12. The elevation Hi is obtained by bilinear interpolation on the reference DEM according to the object plane coordinates of the homologous image points; the original image plane coordinates (ri, ci) of the homologous image points are matched with the non-NAD image to obtain the homologous image points on the non-NAD image; the object coordinates (Li, Bi, Hi) corresponding to the homologous image point coordinates in each perspective together form a control point, and thus N control points are obtained.

[0023] S12. The connection points evenly distributed between the NAD image and the images of the remaining perspectives are separately matched, and these connection points are used for regional network adjustment with the N control points matched previously, and the RPC parameters of each scene of the image are modified by iterative solution.

[0024] S13. Take the corrected RPC parameters and stereo images of each scene as basic data, obtain a dense point cloud through dense matching and forward intersection of homologous image points, and grid the point cloud to obtain a DSM.

[0025] Preferably, in step S12, each time an iteration is performed, homologous image points with errors exceeding a certain threshold are removed through the RANSAC algorithm according to the image-side error of the connection points. At the same time, control points with errors exceeding a certain threshold are removed through the RANSAC algorithm according to the deviation error between the forward intersection coordinates of the control points and the object-side coordinates of the control points. The threshold is gradually reduced with each iteration compared to 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.

[0026] Preferably, step S2 specifically includes the following content.

[0027] S21. Project each image point of the non-NAD perspective image onto the NAD image point by point through the corrected RPC via the DSM, and obtain the corresponding gray value through bilinear interpolation to obtain a projection image with the same length and width as this 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.

[0028] S22. If there is no attitude flutter residual in the image, the column-direction coordinate deviation between the non-NAD image and other projection images will present as a random quantity with a small order of magnitude; if there is an attitude flutter residual, the column-direction coordinate deviation between the non-NAD image and other projection images will present as a systematic quantity corresponding to the flutter amplitude; match the homologous image points of the non-NAD image and other projection images point by point, and calculate the column coordinate difference between the two; after sorting all the column coordinate differences of each row of the image, take the median value as the column-direction flutter image surface error of this row of the image; calculate the flutter image surface error of each row in turn to obtain the column-direction flutter error sequence of the non-NAD image based on the stereo image.

[0029] Preferably, step S3 is specifically as follows.

[0030] Perform a segmented Fourier transform on the column-direction flutter error sequence of the non-NAD image based on the stereo image. The segmented width WF is an integer power of 2; starting from 1, multiply by 2 continuously, and each time a number that is twice the original is obtained until the result is WF; if the number of rows of the non-NAD image is an integer FN times the segmented width, only take the number of segmented width × FN in the column-direction flutter error sequence of the non-NAD image for FN Fourier transforms. The result of each Fourier transform is represented by FD, and each element in FD is a complex number.

[0031] Set up a floating-point array FP with a length of WF, and initialize all elements 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 first element to the 64th element of FP as FA, and the sum of all elements between the i s to the i e as FJ, where i s is greater than 1, i e is less than 64, corresponding to the frequency band of the attitude flutter curve within a certain time period or other specific conditions; calculate the ratio ρ = FJ / FA. If ρ is greater than a certain threshold, it is considered that there is obvious flutter in the image.

[0032] Preferably, step S4 specifically includes the following content,

[0033] S41. Match the homologous image points between the images of different perspectives of the non-flutter image pairs, and at the same time match the homologous image points between the non-flutter image pairs and the flutter image pairs. All overlapping image pairs match the control points by using the method of step S1 through the reference base map and the reference DEM. Adjust the regional network of the non-flutter image pairs and the flutter image pairs together, and then extract the overlapping area DSM by using all non-flutter image pairs;

[0034] Among them, the NAD image and the non-NAD image obtained in the same orbit with an overlapping range exceeding 80% form an image pair; the image pair determined by the flutter significance detection in step S3 is determined as the flutter-significant image pair, called the flutter image pair; at the same time, there is more than one image pair with a certain overlap degree with the flutter image pair but without obvious flutter, called the non-flutter image pair;

[0035] S42. Assume that there is obvious flutter in the non-NAD image of a flutter image pair, and there is no obvious flutter in the M images of the same perspective overlapping with it; first determine the overlapping range between the non-NAD image and a certain overlapping image, and project each image point of the non-NAD image within the overlapping range onto the overlapping image point by point through the overlapping area DSM, and calculate the corresponding gray value by bilinear interpolation, so as to obtain a projection image with the same length and width as the overlapping range of the non-NAD image;

[0036] S43. Match the homologous image points of the non-NAD image and the overlapping image within the overlapping range point by point, calculate the column coordinate difference between the two according to the coordinates of the homologous image points, sort all the column coordinate differences of each row of the image and take the middle value as the column direction flutter error of this row of the image, and calculate the flutter error of each row in turn to obtain a column direction flutter error sequence within the overlapping range;

[0037] S44. Calculate the difference in row coordinates between the two based on the coordinates of the homologous image points. After sorting all the row coordinate differences of each row of the image, take the median value as the row - direction flutter error of that row of the image. Calculate the flutter error of each row in turn to obtain a sequence of row - direction flutter errors within the overlapping range;

[0038] S45. Use the same method to calculate the overlapping range of this non - NAD image and all other overlapping images, calculate the projected images within the overlapping range, and then calculate the sequences of flutter errors in both the column and row directions within the overlapping range;

[0039] S46. For the image rows with multiple flutter values, take the average as the flutter value of that row of the image. For the image rows without overlapping images, set their flutter values to 0. Thus, a sequence of row - direction flutter errors based on the same - perspective images and a sequence of column - direction flutter errors based on the same - perspective images that are consistent with the number of rows of this non - NAD image are obtained;

[0040] S47. Use mean filtering to divide the sequences of flutter errors in the row and column directions of the non - NAD image based on the same - perspective images into the corresponding low - frequency flutter components and high - frequency flutter components in the row and column directions respectively.

[0041] Preferably, step S5 specifically includes the following content:

[0042] S51. Project each grid point of the DSM to the reference DEM according to the coordinates to obtain the corresponding floating - point row and column coordinates. Use 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 this grid point. Calculate the elevation values of all grid points to obtain an error DSM;

[0043] S52. Project each image point of the non - NAD image to the DSM through the RPC after block adjustment. According to the coordinates of the projected points, use bilinear interpolation to obtain the error DSM value; after calculating the corresponding error DSM values for each non - NAD image, obtain an elevation error map with the same length and width as this non - NAD image; after sorting the elevation errors of each row of the non - NAD image, take the median value to obtain the average error of that row of the image; after calculating each row, obtain a sequence of elevation errors with the same length as the number of rows of this non - NAD image;

[0044] S53. Calculate the corresponding relationship between the elevation error and the pixel deviation amount, and obtain the sequence of flutter errors of the non - NAD image measured by the image plane distance based on the corresponding relationship, that is, the sequence of row - direction flutter errors of the non - NAD image based on the stereo images;

[0045] S54. Use mean filtering to divide the column - direction flutter error sequence of the non - NAD image based on the stereo image obtained in step S2 and the row - direction flutter error sequence of the non - NAD image based on the stereo image obtained in step S53 into the corresponding low - frequency flutter components and high - frequency flutter components in the row and column directions respectively.

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

[0047] Preferably, step S6 specifically includes the following contents.

[0048] S61. Perform anti - flutter calculation point - by - point to obtain the non - NAD image with flutter eliminated, that is, the anti - distorted image; For the image point in the r - th row and c - th column of the anti - distorted image, calculate the anti - distorted coordinates (r1, c1).

[0049] S62. According to the anti - distorted coordinates (r1, c1), perform bilinear interpolation or bicubic convolution interpolation on the anti - distorted image, and calculate the gray value and assign it to the point (r, c) of the anti - distorted image.

[0050] S63. After calculating the interpolated gray values of all image points of the anti - distorted image, obtain a new anti - distorted image, which is the anti - distorted image with flutter factors eliminated.

[0051] Preferably, in step S61,

[0052] Without considering the flutter value calculated in step S4, there are three ways to calculate the anti - distorted coordinates.

[0053] (1) Considering both low - frequency flutter and high - frequency flutter, then r1 = r - JFr r , c1 = r - JFc c ;

[0054] (2) Only considering low - frequency flutter, then r1 = r - JFLr r , c1 = r - JFLc c ;

[0055] (3) Only considering high-frequency flutter, then r1 = r - JFHr r and c1 = r - JFHc c ;

[0056] where JFr r and JFc c are respectively the row-direction flutter sequences of the r-th row and the column-direction flutter sequences of the c-th column of the stereo image based on non-NAD images; JFLr r and JFLc c are respectively the row-direction low-frequency flutter components of the r-th row and the column-direction low-frequency flutter components of the c-th column of the stereo image based on non-NAD images; JFHr r and JFHc c are respectively the row-direction high-frequency flutter components of the r-th row and the column-direction high-frequency flutter components of the c-th column of the stereo image based on non-NAD images;

[0057] Considering the flutter values calculated in step S4, there are three ways to calculate the anti-distortion coordinates,

[0058] (1) Considering both the low-frequency flutter and high-frequency flutter in step S4, then r1 = r - JFOr r and c1 = r - JFOc c ;

[0059] (2) Considering both the low-frequency flutter in step S5 and the high-frequency flutter in step S4, then:

[0060] r1 = r - (JFLr r + JFOHr r ), c1 = r - (JFLc c + JFOHc c );

[0061] (3) Only considering the high-frequency flutter in step S4, then r1 = r - JFOHr r and c1 = r - JFOHc c ;

[0062] where JFOr r and JFOc c are respectively the row-direction flutter sequences of the r-th row and the column-direction flutter sequences of the c-th column of the same-view image based on non-NAD images; JFOHr r and JFOHc c are respectively the high-frequency flutter components of the r-th row and the column-direction high-frequency flutter components of the c-th column of the same-view image based on non-NAD images.

[0063] The beneficial effects of the present invention are as follows: 1. By matching dense control points, connection points and performing block adjustment, the systematic error between the DSM obtained from stereo images and the reference DEM is eliminated. 2. Based on the initial DSM corresponding to the stereo images, the projection images are calculated, the deviation of each pixel between the original non-NAD image and its projection image is matched, and the median value of the deviation is calculated row by row, which can reduce the influence of matching errors and other factors on the statistical value of the flutter. 3. The image pairs with significant flutter are quickly determined by the spectral analysis of the attitude residuals in the column direction. 4. According to the elevation error between the DSM and the DEM, median filtering is performed row by row in the image space to reduce the interference of the DSM error factor on the statistical value of the flutter in the row direction. 5. By comprehensively using the flutter values obtained from the same-view images and the flutter values obtained from the stereo images, the flutter error of the images is eliminated. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 It is a flowchart of the method in an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0065] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are only used to explain the present invention, but not to limit the present invention.

[0066] As Figure 1 shown, in this embodiment, a method for detecting and eliminating flutter of satellite stereo images is provided, and the method includes the following six parts:

[0067] I. Block adjustment of the original images and DSM extraction

[0068] Based on the reference base map, reference DEM, the original stereo images to be processed and the RPC files, the homologous pixels of the NAD image and the reference base map are matched and the corresponding object plane coordinates are calculated; the elevation is obtained according to the object plane coordinates of the homologous pixels and the reference DEM, and the image plane homologous pixel coordinates and the corresponding object coordinates of each view are together used to form control points; the connection points between the NAD image and the images of other views are matched, and based on the connection points and the control points, block adjustment is performed to iteratively correct the RPC parameters of each scene image; the DSM is obtained based on the corrected RPC parameters and the stereo images.

[0069] Specifically, it includes the following content.

[0070] 1.1. Using the reference base map, reference DEM, the original image to be processed, and the RPC file as input data, match the uniformly distributed corresponding image points ((ri, ci)-(Ri, Ci), where r and c represent the row and column coordinates of the original image, and R and C represent the row and column coordinates of the reference base map, i = 0, 1, …, N - 1, representing the point sequence number) between the original image with the smallest viewing angle in the original stereo images (hereinafter referred to as the NAD image, and the images of the remaining viewing angles are represented by BWD or FWD, BWD represents the back view, FWD represents the front view. Taking the two-view stereo images of GF-7 as an example, the stereo images of GF-7 include the NAD image and the FWD image, while the three-view stereo images of satellites such as ZY-3 include the NAD image, the BWD image, and the FWD image) and the reference base map. The corresponding object plane coordinates (i.e., longitude Li, latitude Bi) can be calculated from R, C, and the TFW file of the reference base map.

[0071] 1.2. Bilinearly interpolate the elevation Hi on the reference DEM according to the object plane coordinates (longitude Li and latitude Bi) of the corresponding image points. The original image plane coordinates (ri, ci) of the corresponding image points are matched with the non-NAD images to obtain the corresponding image points (rFWDi, cFWDi) on the FWD image and the corresponding image points (rBWDi, cBWDi) on the BWD image (only the three-view stereo images have (rBWDi, cBWDi)). The image plane corresponding image point coordinates of each viewing angle and the corresponding object coordinates (Li, Bi, Hi) together form a control point, and thus N control points are obtained.

[0072] 1.3. Individually match the uniformly distributed connection points between the NAD image and the images of the remaining viewing angles, and perform block adjustment on these connection points and the previously matched N control points. Modify the RPC parameters of each scene image through iterative solution. Each time an iteration is performed, according to the image plane error of the connection points, the corresponding image points with errors exceeding a certain threshold are removed by the RANSAC algorithm. At the same time, according to the deviation error between the forward intersection coordinates of the control points and the object coordinates of the control points, the control points with errors exceeding a certain threshold are removed by the RANSAC algorithm. The threshold is gradually reduced with each iteration compared to 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.

[0073] 1.4. Using the corrected RPC of each scene and the stereo images as the basic data, obtain the dense point cloud through dense matching and forward intersection of corresponding image points, and grid the point cloud to obtain the DSM.

[0074] II. Calculation of the column-direction flutter error sequence based on stereo images for non-NAD images

[0075] Obtain the projection image of the non - NAD image, and calculate the column coordinate difference between the non - NAD image and its projection image by point - by - point matching of homologous image points; obtain the column - direction flutter image plane error of each row of the image based on the column coordinate difference, and obtain the column - direction flutter error sequence of the non - NAD image based on the stereo image based on the flutter image plane error of each row.

[0076] Specifically, it includes the following content:

[0077] 2.1. Project each image point of the non - NAD perspective image (such as FWD image or BWD image) onto the NAD image point - by - point through the corrected RPC via DSM, and obtain the corresponding gray value by bilinear interpolation to obtain a projection image with the same length and width as this 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).

[0078] 2.2. Taking the FWD image and its projection image as an example, if there is no attitude flutter residual in the image, the column - direction coordinate deviation between the FWD image and its projection image will present as a random quantity with a very small magnitude. If there is an attitude flutter residual, the column - direction coordinate deviation between the FWD image and its projection image will present as a systematic quantity corresponding to the flutter amplitude. Point - by - point match the homologous image points of the FWD image and its projection image, and calculate the column coordinate difference between the two. After sorting all the column coordinate differences of each row of the image, take the median value as the column - direction flutter image plane error of this row of the image. Calculate the flutter image plane error of each row in turn to obtain a column - direction flutter error sequence JFc i (i = 0, 1, …, HF - 1), that is, the column - direction flutter error sequence of the FWD image based on the stereo image, and HF is the number of rows of the FWD image.

[0079] 2.3. If there are non - NAD images of other perspectives, calculate the column - direction flutter error sequence of the remaining non - NAD perspective images based on the stereo image in the same way, such as the column - direction flutter error sequence JBc of the BWD image i (i = 0, 1, …, HB - 1), and HB is the number of rows of the BWD image.

[0080] III. Flutter significance detection of non - NAD images

[0081] Perform a segmented Fourier transform on the column - direction flutter error sequence of the non - NAD image based on the stereo image to determine whether there is obvious flutter in the non - NAD image, and realize the flutter significance detection of the non - NAD image.

[0082] Specifically, it includes the following content:

[0083] 3.1. Taking the FWD image as an example, perform segmented Fourier transform on JFc i (i = 0, 1, …, HF - 1), where the segmented width is a power of 2, denoted as WF, and assume WF is 1024. Starting from 1, continuously multiply by 2, each time obtaining a number that is twice the original, until the result is WF, and record the number of multiplications as w. Assume HF is an integer FN times WF, then only take WF×FN numbers in JFc i (i = 0, 1, …, HF - 1) to perform FN Fourier transforms. The result of each Fourier transform is denoted as FD i (i = 0, 1, …, WF - 1), and each element in FD is a complex number. Set up a floating-point array FP i (i = 0, 1, …, WF - 1) with a length of WF and all elements initialized to 0. After each Fourier transform, each element of FP i (i = 0, 1, …, WF - 1) is added with the modulus of the corresponding complex number in FD i (i = 0, 1, …, WF - 1), 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 first 64 elements of FP i (i = 0, 1, …, WF - 1) as FA, and the sum of all elements from the i s - th to the i e - th as FJ. Where i s is greater than 1, i e is less than 64, for example, when i s equals 17, i e equals 21, which are empirical values corresponding to the frequency band of the attitude flutter curve within a certain time period or other specific conditions. Calculate the ratio ρ = FJ / FA. If ρ is greater than a certain threshold, such as ρ > 0.3, then it is considered that there is obvious flutter in the image.

[0084] 3.2. Use the same method to perform flutter significance detection on other non - NAD images.

[0085] IV. Calculation of the row - column flutter error sequence of non - NAD images based on images of the same perspective

[0086] By matching the homologous image points between different - perspective images of non - flutter image pairs and between non - flutter image pairs and flutter image pairs, obtaining control points based on overlapping image pairs, reference base maps, and reference DEMs, performing block adjustment based on non - flutter image pairs and flutter image pairs to extract the DEM of the overlapping area; obtaining the projected image based on the DEM of the overlapping area; calculating the column - coordinate difference and row - coordinate difference between the non - NAD image within the overlapping range and the projected image based on the homologous image points to obtain the row - column direction flutter error sequence within the overlapping range; obtaining the row - column flutter error sequence of non - NAD images based on images of the same perspective by resetting the flutter value.

[0087] Specifically, it includes the following content:

[0088] 4.1. Assume that a pair of images (an NAD image and a non-NAD image obtained in the same orbit with an overlap range exceeding 80% with it form a pair of images) passes through the flutter significance detection in step S3 and is determined to be a pair of images with significant flutter (abbreviated as flutter image pair). At the same time, there is more than one pair of images with a certain overlap degree with this pair of images (the overlap degree is calculated according to the overlap degree of the NAD images of each pair of images), but there is no significant flutter (hereinafter referred to as non-flutter image pair). Match the corresponding homologous image points between the different perspective images of these non-flutter image pairs, and at the same time match the corresponding homologous image points between the non-flutter image pairs and the flutter image pairs. All overlapping image pairs match control points by the method of step S1 using the reference base map and the reference DEM. Adjust the regional network of the non-flutter image pairs and the flutter image pairs together, and then extract the overlapping area DSM using all non-flutter image pairs.

[0089] 4.2. Assume that there is obvious flutter in the FWD image of a flutter image pair, and there is no obvious flutter in the M co-perspective images overlapping with it. First, determine the overlapping range of the FWD image and the i-th (i = 0, 1,..., M - 1) overlapping image (referred to as the FWD i image). Project each image point in the overlapping range of the FWD image onto the FWD i image point by point through the overlapping area DSM, and calculate the corresponding gray value by bilinear interpolation, so as to obtain a projected image with the same length and width as the overlapping range of the FWD image.

[0090] 4.3. Match the corresponding homologous image points of the FWD image and the FWD i projected image in the overlapping range point by point, and calculate the difference in column coordinates between the two according to the coordinates of the corresponding homologous image points. After sorting all the column coordinate differences of each row of the image, take the middle value as the flutter error in the column direction of this row of the image. Calculate the flutter error of each row in turn to obtain a flutter error sequence in the column direction within the overlapping range.

[0091] 4.4. Calculate the difference in row coordinates between the two according to the coordinates of the corresponding homologous image points. After sorting all the row coordinate differences of each row of the image, take the middle value as the flutter error in the row direction of this row of the image. Calculate the flutter error of each row in turn to obtain a flutter error sequence in the row direction within the overlapping range.

[0092] 4.5. Calculate the overlapping range of the FWD image and the 1st to the M - 1th overlapping co-perspective images in the same way, calculate the projected images within the overlapping range, and then calculate the flutter error sequences in the column and row directions within the overlapping range.

[0093] 4.6 After calculating the flutter error sequences in the column and row directions for all overlapping images, there may be multiple flutter values for some rows. Take the average as the flutter value for the images in that row. For some image rows without overlapping images, set their flutter values to 0. Thus, a row-direction flutter error sequence JFOr based on the same-view images with the same number of rows as the FWD image is obtained. i (i = 0, 1, …, HF - 1), and a column-direction flutter error sequence JFOc based on the same-view images. i (i = 0, 1, …, HF - 1), where HF is the number of rows of the FWD image.

[0094] 4.7 Calculate the flutter error sequences based on the same-view images for other non-NAD images in the same way, such as the flutter error sequences in the row and column directions of the BWD image. The row-direction flutter error sequence of the BWD image is JBOr i (i = 0, 1, …, HB - 1), and the column-direction flutter error sequence is JBOc i (i = 0, 1, …, HB - 1), where HB is the number of rows of the BWD image.

[0095] 4.8 Mean filtering divides the flutter into low-frequency and high-frequency components: Taking the FWD image as an example, perform mean filtering on the JFOc i (i = 0, 1, …, HF - 1) sequence in the column direction based on the same-view images. Set the filtering width to 200 to obtain the low-frequency flutter component JFOLc in the column direction i (i = 0, 1, …, HF - 1). Subtract each element of the JFOc i (i = 0, 1, …, HF - 1) sequence from the corresponding element of the JFOLc i (i = 0, 1, …, HF - 1) to obtain the high-frequency flutter component JFOHc in the column direction i (i = 0, 1, …, HF - 1).

[0096] 4.9 Obtain the low-frequency flutter component JFOLr in the row direction based on the same-view images and the high-frequency flutter component JFOHr i (i = 0, 1, …, HF - 1) and i (i = 0, 1, …, HF - 1) in the same way.

[0097] V. Calculation of the row-direction flutter error sequence based on stereo images for non-NAD images

[0098] By projecting non-NAD images onto the error DSM obtained from the elevation of grid points in the DSM and the elevation of the reference DEM, the error DSM values are obtained; using the error DSM to obtain an elevation error sequence, and by calculating the correspondence between the elevation error and the pixel deviation amount, the stereo-image-based row-direction flutter error sequence of the non-NAD image is obtained; the stereo-image-based row and column direction flutter error sequences of the non-NAD image are respectively divided into low-frequency flutter components and high-frequency flutter components in the corresponding row and column directions.

[0099] 5.1 Project each grid point of the DSM onto the reference DEM according to the coordinates to obtain the corresponding floating-point row and column coordinates, and bilinearly interpolate 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 this grid point. Calculate the elevation values of all grid points to obtain an error DSM.

[0100] 5.2 Project each pixel of the FWD image onto the DSM through the RPC after block adjustment correction, and bilinearly interpolate according to the projected point coordinates to obtain the error DSM value. 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 the FWD and take the median value to obtain the average error of this row of the image. After calculating each row, an elevation error sequence JFHr i (i = 0, 1, …, HF - 1), where HF is the number of rows of the FWD image.

[0101] 5.3 Calculate the correspondence between the elevation error and the pixel deviation amount: Assume that the HEIGHT_OFF parameter of the corrected RPC of the NAD image is Ho, and the center row and column of the NAD are rc and cc. Calculate the corresponding row and column coordinates rf and cf of the FWD image according to rc, cc, Ho and the corrected RPC of the NAD image and the corrected RPC of the FWD image. In the same way, calculate the corresponding row and column coordinates rf1 and cf1 of the FWD image according to rc, cc, Ho + 1 and the corrected RPC of the NAD image and the corrected RPC of the FWD image. Calculate the image plane distance Df between the FWD image plane coordinates (rf, cf) and (rf1, cf1).

[0102] 5.4 Multiply each element of JFHr i (i = 0, 1, …, HF - 1) by Df to obtain the flutter error sequence JFr i (i = 0, 1, …, HF - 1), that is, the stereo-image-based row-direction flutter error sequence.

[0103] 5.5 Calculate the stereo-image-based flutter error sequences of the remaining non-NAD images in the same way, such as the stereo-image-based flutter error sequence JBr of the BWD image i(i = 0, 1, …, HF - 1).

[0104] 5.6. Mean filtering divides flutter into low - frequency components and high - frequency components: Taking the FWD image as an example, perform mean filtering on its JFc i (i = 0, 1, …, HF - 1) sequence in the column direction based on the stereo image, set the filtering width to 200, and obtain the low - frequency flutter component JFLc i (i = 0, 1, …, HF - 1) in the column direction. Subtract each element of the JFc i (i = 0, 1, …, HF - 1) sequence from the corresponding element of JFLc i (i = 0, 1, …, HF - 1) to obtain the high - frequency flutter component JFHc i (i = 0, 1, …, HF - 1) in the column direction of the FWD image.

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

[0106] VI. Flutter error elimination for non - NAD images

[0107] Obtain the non - NAD image with flutter eliminated, i.e., the anti - distorted image, through anti - flutter calculation; calculate the anti - distorted coordinates of each image point of the anti - distorted image under the condition of considering or not considering the low - frequency flutter component and / or high - frequency flutter component in step S4 and / or step S5; obtain the interpolated gray value of all image points of the anti - distorted image based on the anti - distorted coordinates, and further obtain a new anti - distorted image, i.e., the anti - distorted image with flutter factors eliminated.

[0108] Specifically, it includes the following content.

[0109] 6.1. Perform anti - flutter calculation point by point to obtain the non - NAD image with flutter eliminated (hereinafter referred to as the anti - distorted image); taking the FWD image as an example, for the image point at the r - th row and c - th column of the anti - distorted image, calculate the anti - distorted coordinates (r1, c1). There are two cases for the anti - distorted coordinates:

[0110] A. Without considering the flutter value calculated in the fourth part, there are three calculation methods.

[0111] (1) Considering both low - frequency flutter and high - frequency flutter, then r1 = r - JFr r , c1 = r - JFc c .

[0112] (2) Only considering low - frequency flutter, then r1 = r - JFLr r , c1 = r - JFLcc .

[0113] (3) Only considering high-frequency flutter, then r1 = r - JFHr r , c1 = r - JFHc c .

[0114] B. Considering the flutter values calculated in the fourth part, there are three calculation methods,

[0115] (1) Considering both the low-frequency flutter and high-frequency flutter in step S4, then r1 = r - JFOr r , c1 = r - JFOc c .

[0116] (2) Considering both the low-frequency flutter in step S5 and the high-frequency flutter in step S4, then:

[0117] r1 = r - (JFLr r + JFOHr r ), c1 = r - (JFLc c + JFOHc c ).

[0118] (3) Only considering the high-frequency flutter in step S4, then r1 = r - JFOHr r , c1 = r - JFOHc c .

[0119] 6.2. Perform bilinear interpolation or bicubic convolution interpolation on the FWD image according to the anti-distortion coordinates (r1, c1), and calculate the gray value and assign it to the point (r, c) of the anti-distorted FWD image.

[0120] 6.3. After calculating the interpolated gray values of all the pixel points of the FWD anti-distorted image, a new FWD image is obtained, which is the anti-distorted image after eliminating the flutter factor.

[0121] 6.4. Calculate the anti-distorted images of the remaining non-NAD perspective images in the same way.

[0122] By adopting the above technical solutions disclosed in the present invention, the following beneficial effects are obtained:

[0123] The present invention provides a method for detecting and eliminating satellite stereo image flutter. By matching dense control points, tie points and performing block adjustment, the systematic error between the DSM obtained from the stereo image and the reference DEM is eliminated. 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 value of the deviation row by row, which can reduce the influence of factors such as matching error on the flutter value statistics. The present invention quickly determines the image pairs with significant flutter through the spectral analysis of the attitude residuals in the column direction. The present invention performs median filtering row by row in the image plane according to the elevation error between the DSM and the DEM, reducing the interference of the DSM error factor on the flutter value statistics in the row direction. The present invention comprehensively utilizes the flutter values obtained from the same-viewpoint images and the flutter values obtained from the stereo images to achieve the elimination of image flutter errors.

[0124] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. A method for detecting and eliminating satellite stereo image flutter, characterized in that: It includes the following steps: S1. Original image block adjustment and DSM extraction: Based on the reference base map, reference DEM, original stereo images to be processed and RPC files, match the homologous image points of 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 homologous image points and the reference DEM, and form control points together with the image plane coordinates of the homologous image points of each view; match the connection points between the NAD image and other view images, and perform block adjustment iteration based on the connection points and control points to correct the RPC parameters of each scene image; obtain the DSM based on the corrected RPC parameters and stereo images; S2. Calculation of column-direction flutter error sequence of non-NAD images based on stereo images: Obtain the projected image of the non-NAD image, and calculate the column coordinate difference between the non-NAD image and its projected image by point-by-point matching of the homologous image points; obtain the column-direction flutter image plane error of each row of images based on the column coordinate difference, and obtain the column-direction flutter error sequence of the non-NAD image based on stereo images based on the flutter image plane error of each row; S3. Flutter significance detection of non-NAD images: Perform piecewise Fourier transform on the column-direction flutter error sequence of the non-NAD image based on stereo images to determine whether there is obvious flutter in the non-NAD image, and realize the flutter significance detection of the non-NAD image; S4. Calculation of row-column flutter error sequence of non-NAD images based on same-view images: By matching the homologous image points between different view images of non-flutter image pairs and between non-flutter image pairs and flutter image pairs, obtain control points based on overlapping image pairs, reference base map and reference DEM, and perform block adjustment based on non-flutter image pairs and flutter image pairs to extract the overlapping area DEM; Obtain the projected image based on the overlapping area DEM; calculate the column coordinate difference and row coordinate difference between the non-NAD image within the overlapping range and the projected image based on the homologous image points, and obtain the row-column direction flutter error sequence within the overlapping range; obtain the row-column flutter error sequence of the non-NAD image based on same-view images by resetting the flutter value; S5. Calculation of row-direction flutter error sequence of non-NAD images based on stereo images: Project the non-NAD image into the error DSM obtained by using the grid point elevation of the DSM and the elevation of the reference DEM to obtain the error DSM value; use the error DSM to obtain the elevation error sequence, and obtain the row-direction flutter error sequence of the non-NAD image based on stereo images by calculating the corresponding relationship between the elevation error and the pixel deviation amount; divide the row-column direction flutter error sequence of the non-NAD image based on stereo images into low-frequency flutter components and high-frequency flutter components in the corresponding row and column directions respectively; S6. Elimination of flutter error of non-NAD images Obtain the anti-flutter non-NAD image based on anti-flutter calculation, that is, the anti-distortion image; calculate the anti-distortion coordinates of each image point of the anti-distortion image under the condition of considering or not considering the low-frequency flutter component and / or high-frequency flutter component in step S4 and / or step S5; obtain the interpolated gray value of all image points of the anti-distortion image based on the anti-distortion coordinates, and then obtain a new anti-distortion image, that is, the anti-distortion image eliminating the flutter factor.

2. The satellite stereo image flutter detection and elimination method according to claim 1, characterized in that: Step S1 specifically includes the following content: S11. Take the reference base map, reference DEM, the original image to be processed and the RPC file as input data, and match the evenly distributed homologous image points (ri, ci)-(Ri, Ci) between the NAD image and the reference base map; calculate 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 viewing angle in the original stereo image is called the NAD image, and the images of the remaining viewing angles are non-NAD images, denoted by BWD or FWD, where BWD represents the back view and FWD represents the front view; r and c respectively 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 respectively represent longitude and latitude. S12. Bilinearly interpolate the elevation Hi on the reference DEM according to the object plane coordinates of the homologous image points; match the original image plane coordinates (ri, ci) of the homologous image points with the non-NAD image to obtain the homologous image points on the non-NAD image; the object coordinates (Li, Bi, Hi) corresponding to the image-side homologous image points of each viewing angle together form a control point, and thus N control points are obtained. S12. Independently match the evenly distributed connection points between the NAD image and the images of the remaining viewing angles, perform 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. S13. Take the RPC parameters and stereo images of each corrected scene as basic data, obtain a dense point cloud through dense matching and forward intersection of homologous image points, and grid the point cloud to obtain the DSM.

3. The satellite stereo image flutter detection and elimination method according to claim 2, characterized in that: In step S12, each time an iteration is performed, homologous image points with errors exceeding the preset threshold are removed through the RANSAC algorithm according to the image-side error of the connection points. At the same time, control points with errors exceeding the preset threshold are removed through the RANSAC algorithm according to the deviation error between the forward intersection coordinates of the control points and the object coordinates of the control points. The threshold is gradually reduced each time 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 flutter detection and elimination method according to claim 3, characterized in that: Step S2 specifically includes the following content: S21. Project each image point of the non-NAD view image onto the NAD image point by point through the corrected RPC via the DSM, and bilinearly interpolate to obtain the corresponding gray value, obtaining a projection image with the same length and width as this 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 flutter residual in the image, the column-direction coordinate deviation between the non-NAD image and other projection images will present as a small random quantity of the same order of magnitude; if there is an attitude flutter residual, the column-direction coordinate deviation between the non-NAD image and other projection images will present as a systematic quantity corresponding to the flutter amplitude. Match the homologous image points of the non-NAD image and other projection images point by point, and calculate the column coordinate difference between the two. After sorting all the column coordinate differences of each row of the image, take the median value as the column-direction flutter image plane error of this row of the image. Calculate the flutter image plane error of each row in turn to obtain the column-direction flutter error sequence of the non-NAD image based on the stereo images.

5. The satellite stereo image flutter detection and elimination method according to claim 4, characterized in that: Step S3 specifically is as follows. Perform a segmented Fourier transform on the column-direction flutter error sequence of the non-NAD image based on the stereo images. The segmented width WF is an integer power of 2. Starting from 1, multiply by 2 continuously, and each time get a number that is twice the original until the result is WF. If the number of rows of the non-NAD image is an integer FN times the segmented width, then only take the number of elements equal to the segmented width × FN in the column-direction flutter error sequence of the non-NAD image for FN Fourier transforms. The result of each Fourier transform is represented by FD, and each element in FD is a complex number. Set up a floating-point array FP with a length of WF, where all elements are 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 first 64 elements of FP as FA, and the sum of all elements between the i s to the i e as FJ, where i s is greater than 1 and i e is less than 64, corresponding to the frequency band of the attitude flutter curve within the preset time period; Calculate the ratio ρ = FJ / FA. If ρ is greater than the preset threshold, it is considered that there is obvious flutter in this image.

6. The method for detecting and eliminating satellite stereo image flutter according to claim 5, characterized in that: Step S4 specifically includes the following content. S41. Match the homologous image points between the different perspective images of the non-flutter image pairs, and at the same time match the homologous image points between the non-flutter image pairs and the flutter image pairs. For all overlapping image pairs, match the control points using the method in step S1 through the reference base map and reference DEM. Perform block adjustment on the non-flutter image pairs and the flutter image pairs together, and then extract the overlapping area DSM using all the non-flutter image pairs. Among them, the NAD image and the non-NAD image obtained from the same orbit with an overlapping range exceeding 80% form an image pair. The image pair determined by the flutter significance detection in step S3 is determined as the flutter-significant image pair, called the flutter image pair. At the same time, there is more than one image pair that has a preset overlap degree with the flutter image pair but does not have significant flutter, which is called the non-flutter image pair. S42. Assume that there is obvious flutter in the non-NAD image of a flutter image pair, and there is no obvious flutter in the M co-perspective images overlapping with it. First, determine the overlapping range between the non-NAD image and a certain overlapping image. Project each image point of the non-NAD image within the overlapping range onto the overlapping image point by point through the overlapping area DSM, and calculate the corresponding gray value by bilinear interpolation, so as to obtain a projection image with the same length and width as the overlapping range of the non-NAD image. S43. Match the homologous image points of the non-NAD image and the overlapping image within the overlapping range point by point. Calculate the column coordinate difference between the two according to the coordinates of the homologous image points. After sorting all the column coordinate differences of each row of the image, take the median value as the column-direction flutter error of this row of the image. Calculate the flutter error of each row in turn to obtain a column-direction flutter error sequence within the overlapping range. S44. Calculate the difference in row coordinates between the two based on the coordinates of homologous image points. After sorting all the row coordinate differences of each row of the image, take the median value as the row-direction flutter error of that row of the image. Calculate the flutter error of each row in turn to obtain a sequence of row-direction flutter errors within the overlapping range. S45. Calculate the overlapping range of this non-NAD image and all other overlapping images in the same way. Calculate the projected image within the overlapping range, and then calculate the sequences of flutter errors in the column and row directions within the overlapping range. S46. For the image rows with multiple flutter values, take the average as the flutter value of that row of the image. For the image rows without overlapping images, set their flutter values to 0. Thus, a sequence of row-direction flutter errors based on same-viewpoint images and a sequence of column-direction flutter errors based on same-viewpoint images that are consistent with the number of rows of this non-NAD image are obtained. S47. Use mean filtering to divide the sequences of flutter errors in the row and column directions of the non-NAD image based on same-viewpoint images into corresponding low-frequency and high-frequency flutter components in the row and column directions respectively.

7. The satellite stereo image flutter detection and elimination method according to claim 6, characterized in that: Step S5 specifically includes the following contents. S51. Project each grid point of the DSM to the reference DEM according to the coordinates 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 reference DEM from the elevation of the DSM grid point to obtain the error value of this grid point. Calculate the elevation values of all grid points to obtain an error DSM. S52. Project each image point of the non-NAD image to the DSM through the RPCs corrected by block adjustment. Bilinear interpolation based on the projected point coordinates to obtain the error DSM value. After calculating the corresponding error DSM values for each non-NAD image, an elevation error map with the same length and width as this non-NAD image is obtained. After sorting the elevation errors of each row of the non-NAD image and taking the median value, the average error of that row of the image is obtained. After calculating each row, an elevation error sequence with the same length as the number of rows of this non-NAD image is obtained. S53. Calculate the corresponding relationship between the elevation error and the pixel deviation amount, and obtain the sequence of flutter errors of the non-NAD image measured by the image plane distance based on this corresponding relationship, that is, the sequence of row-direction flutter errors of the non-NAD image based on stereo images. S54. Use mean filtering to divide the sequence of column-direction flutter errors of the non-NAD image based on stereo images obtained in step S2 and the sequence of row-direction flutter errors of the non-NAD image based on stereo images obtained in step S53 into corresponding low-frequency and high-frequency flutter components in the row and column directions respectively.

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

9. The satellite stereo image flutter detection and elimination method according to claim 8, characterized in that: Step S6 specifically includes the following contents: S61. Perform anti-flutter calculation point by point to obtain a non-NAD image with flutter eliminated, that is, an anti-distortion image; for the image point at the r-th row and c-th column of the anti-distortion image, calculate the anti-distortion coordinates (r1, c1). S62. Perform bilinear interpolation or bicubic convolution interpolation on the anti-distortion image according to the anti-distortion coordinates (r1, c1), and calculate the gray value and assign it to the point (r, c) of the anti-distortion image. S63. After calculating the interpolated gray levels of all image points of the anti-distortion image, obtain a new anti-distortion image, which is the anti-distortion image with flutter factors eliminated.

10. The method for satellite stereo image flutter detection and elimination according to claim 9, wherein: In step S61, Without considering 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, then r1 = r - JFr r , c1 = r - JFc c ; (2) Only considering low-frequency flutter, then r1 = r - JFLr r , c1 = r - JFLc c ; (3) Considering only high-frequency flutter, then r1 = r - JFHr r , c1 = r - JFHc c ; Among them, JFr r and JFc c are respectively the row-direction flutter sequence of the r-th row and the column-direction flutter sequence of the c-th column of the stereo image-based non-NAD image; JFLr r and JFLc c are respectively the row-direction low-frequency flutter component of the r-th row and the column-direction low-frequency flutter component of the c-th column of the stereo image-based non-NAD image; JFHr r and JFHc c are respectively the row-direction high-frequency flutter component of the r-th row and the column-direction high-frequency flutter component of the c-th column of the stereo image-based non-NAD 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, then r1 = r - JFOr r , c1 = r - JFOc c ; (2) Considering both the low-frequency flutter in step S5 and the high-frequency flutter in step S4, then: r1 = r - (JFLr r + JFOHr r ), c1 = r - (JFLc c + JFOHc c ); (3) Considering only the high-frequency flutter in step S4, then r1 = r - JFOHr r , c1 = r - JFOHc c ; Among them, JFOr r and JFOc c are respectively the row-direction flutter sequence of the r-th row and the column-direction flutter sequence of the c-th column of the same-view image based on the non-NAD image; JFOHr r and JFOHc c are respectively the high-frequency flutter component of the r-th row and the column-direction high-frequency flutter component of the c-th column of the same-view image based on the non-NAD image.

Images