A method for automatically measuring the relative image shift in the overlap area of ​​CCD camera stitching focal planes

Abstract

The present invention relates to a method for automatically measuring the relative image shift of the overlap area of ​​the CCD camera splicing focal plane, comprising the steps of: cutting out the overlap areas on both sides of the column position where the image overlap area boundary is located to obtain a pair of sub-areas to be screened; extracting local feature points based on the Fast operator to screen out overlap sub-area pairs whose number of local feature points is greater than a feature point threshold; performing edge detection based on the Canny operator to screen out overlap sub-area pairs whose number of edge pixels is greater than an edge threshold; using the product of the number of local feature points and the number of edge pixels as a screening index, selecting the overlap sub-area pair corresponding to the maximum value of the screening index as the optimal overlap sub-area pair; performing phase correlation registration on the optimal overlap sub-area pair to obtain the relative image shift between the two images. The present invention can automatically, quickly and accurately measure the relative image shift of the overlap area, greatly improving the reliability of the measurement results and meeting the requirements of high-precision time-series monitoring of the relative image shift of the overlap area of ​​a large number of satellites.

Description

Technical Field

[0001] The invention belongs to the field of aerospace technology, and in particular relates to a method for automatically measuring the relative image displacement of an overlapping area of ​​a CCD camera splicing focal plane. Background Art

[0002] To address the limitations of single-chip detectors for high-resolution imaging, satellites equipped with CCD cameras typically employ a multi-detector stitching method to capture wide field-of-view images. Each sub-detector is overlapped and collinearly arranged in the focal plane, independently capturing images of the Earth. Due to varying focal plane positions and the influence of factors such as the payload's flight attitude and ground undulation, this non-collinear multi-detector imaging pattern results in ground points having varying offsets between detectors. Even when fully compensated for image motion caused by payload flight, this relative image motion can cause artifacts such as heading and lateral seams, image point overlap, and misalignment within the overlapped regions when stitching the image into a wide-width image, reducing resolution. Given the precise relative image motion of each overlapped region, stitching techniques can be used to create a single wide-width image from images captured from different fields of view. However, the relative image motion of the overlapped regions can change with subsequent satellite on-orbit commissioning, changes in satellite parameters, and aging of satellite firmware, necessitating high-precision, long-term monitoring.

[0003] With the rapid development of space satellites and the growing demand for remote sensing applications, aerospace technology has entered a phase of large-scale application in large satellite constellations, characterized by multi-satellite networking and multi-network integration. Conventional methods for measuring relative image motion in the overlap zone are unable to meet the large-scale, high-precision time-series monitoring requirements of satellite constellations. High-precision, automated measurement of relative image motion in the overlap zone of CCD camera stitching focal planes on optical remote sensing satellites is crucial for meeting the time-series monitoring requirements of large-scale overlap zones in networked constellations.

[0004] Conventional manual measurement of the relative image shift in the overlap zone requires visual measurement of each image one by one, which is time-consuming and labor-intensive, with low accuracy and significant subjective influence. In addition, there are many schemes for determining the image shift in the overlap zone through matching algorithms. Most of these schemes are obtained by matching the same-name points between the overlap pixels of two CCDs, requiring all pixels on each CCD to be read and matched. This method is time-consuming and has complex algorithms. In addition, when there are unsuitable ground objects in the overlap zone, the calculation results have relatively large errors, making it difficult to meet the requirements of high-precision time-series monitoring of overlap zone image shift for large numbers of satellites. Summary of the Invention

[0005] Aiming at the demand for time-series monitoring of the relative image motion of the overlap areas of large quantities of satellites equipped with CCD cameras when networking constellations, the present invention provides a method for automatically measuring the relative image motion of the overlap areas of the stitching focal plane of CCD cameras based on the idea of ​​first extracting the suitable measurement area and then aligning it. This method can automatically screen the areas most suitable for measuring the overlap areas in images of unknown shooting scenes and measure the relative image motion of the overlap areas, thereby obtaining high-precision sub-pixel-level time-series monitoring results of the relative image motion of the overlap areas of satellites.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A method for automatically measuring the relative image shift of the overlap area of ​​a CCD camera stitching focal plane comprises the following steps:

[0008] Step 1: According to the column positions of the adjacent overlapping area boundaries in the known remote sensing image, the preset widths are extended to the left and right to obtain the left overlapping area image and the right overlapping area image respectively. Then, the left overlapping area image and the right overlapping area image are cropped without overlap according to the set cropping step size to obtain the overlapping sub-area pairs to be screened;

[0009] Step 2: Calculate the number of local feature points for each overlapping sub-region pair based on the Fast operator, and use the local feature point threshold screening method to screen out overlapping sub-region pairs whose local feature points are greater than the feature point threshold;

[0010] Step 3: Perform edge detection on the overlapping sub-region pairs selected in step 2 based on the Canny operator to obtain the corresponding number of edge pixels, and use the edge operator threshold screening method to screen out overlapping sub-region pairs whose edge pixel number is greater than the edge threshold;

[0011] Step 4: Sort the overlapping sub-region pairs selected in step 3 using the product of the number of local feature points and the number of edge pixels as the screening index, and select the overlapping sub-region pair corresponding to the maximum value of the screening index as the optimal overlapping sub-region pair;

[0012] Step 5: Perform phase correlation registration on the optimal overlapping sub-region pair, and obtain the relative image shift between the two images after registration.

[0013] Technical effect of the present invention: The method for automatically measuring the relative image motion of the overlap area of ​​the spliced ​​focal plane of the CCD camera of the optical remote sensing satellite proposed in the present invention can realize automatic, rapid and high-precision measurement of the relative image motion of the overlap area, greatly improve the reliability of the measurement results, and meet the task requirements of high-precision time-series monitoring of the relative image motion of the overlap area of ​​a large number of satellites. Compared with the classic overlap area relative image motion measurement method, this method automatically determines whether the image is suitable for measurement, and selects a suitable measurement area on the image that is suitable for measurement, thereby obtaining high-precision measurement results, thereby realizing the measurement of the relative image motion of the overlap area of ​​a large number of satellites in the operation of the remote sensing constellation. The measurement method has high accuracy, stable measurement results, is easy to implement, and can be effectively applied in engineering practice. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 This is a flow chart of a method for automatically measuring the relative image shift in the overlap area of ​​a CCD camera stitching focal plane according to an embodiment of the present invention;

[0015] Figure 2 This is a flow chart of automatically measuring the relative image shift of the overlapped area of ​​a single image using the method of the present invention, taking an image formed by splicing three detectors of left, middle and right as an example;

[0016] Figure 3 The results are time series monitoring results of the relative image motion of the single-star overlap area obtained using the method of the present invention. DETAILED DESCRIPTION

[0017] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and preferred embodiments.

[0018] like Figure 1 FIG2 is a flow chart of a method for automatically measuring the relative image shift of the overlapped area of ​​the spliced ​​focal plane of a CCD camera of an optical remote sensing satellite proposed by the present invention. The method specifically includes the following steps:

[0019] Step 1: Crop the overlapping areas on both sides of the remote sensing image overlapping area boundary to obtain the sub-area pairs to be screened.

[0020] First, based on the column positions of the adjacent overlap regions in the known remote sensing image, the left and right overlap region images are cropped by extending the preset width horizontally to the left and right, respectively, to obtain the left overlap region image and the right overlap region image. Then, the left overlap region image and the right overlap region image are cropped vertically according to the set cropping step size without overlap to obtain the overlap sub-region pairs to be screened. Each overlap sub-region pair includes a sub-region cropped from the left overlap region image and a sub-region cropped from the corresponding right overlap region image.

[0021] Step 2: Extract corner points based on the Fast operator.

[0022] This step calculates the number of local feature points for each overlapping sub-area pair obtained by cropping in step 1 based on the Fast operator, and uses the local feature point threshold screening method to screen out overlapping sub-area pairs with local feature points greater than the threshold.

[0023] When calculating the number of local feature points based on the Fast operator for overlapping sub-areas, first extract the corner points as local feature points based on the Fast operator. For the central pixel p, draw a circle with a radius of 3 pixels with p as the center, and there are 16 pixels on the circumference. By considering the 16 pixels on the circular window near the pixel p, if there are n (n<16) consecutive pixels with a greater or lesser intensity than the central pixel p, such a center can be preliminarily determined as a corner point. Calculate the state S of each point separately p→x , we can now obtain a set of states for each position p in the image, where each state refers to the state of a pixel and its neighboring (circular) pixels at that position. For strong intensities, a threshold t is required.

[0024]

[0025] Among them, I p is the gray value of pixel p, I p→x is the grayscale mean of the 16 points around p, t is the threshold, and d, s, and b are constants.

[0026] After counting the number of extracted local feature points, the number F of local feature points of the overlapping sub-region pair is obtained.

[0027] After calculating the number of local feature points for each cropped overlapping sub-area pair, the local feature point threshold screening method is used to perform the next calculation for the overlapping sub-area pairs whose local feature points are greater than the feature point threshold. Otherwise, the overlapping sub-area pair is considered unsuitable for measuring the offset, and the offset is recorded as an invalid value. Finally, the overlapping sub-area pairs whose local feature points are greater than the feature point threshold are screened out.

[0028] Step 3: Perform edge detection based on the Canny operator.

[0029] This step uses the Canny operator to perform edge detection on each overlapping sub-region pair whose number of local feature points exceeds the feature point threshold, obtaining the corresponding number of edge pixels. Based on the obtained number of edge pixels corresponding to each overlapping sub-region pair, the edge operator threshold screening method is then used to select overlapping sub-region pairs whose number of edge pixels exceeds the edge threshold.

[0030] Specifically, edge detection is performed on each of the selected overlapping sub-region pairs, including the following steps:

[0031] Step 3-1: Perform Gaussian filtering on the overlapping sub-region pairs to obtain a filtered image.

[0032] Since the edges of images are easily affected by noise, it is usually necessary to filter the image to remove the noise. Gaussian filtering is performed on the image to smooth some non-edge areas with weak textures to obtain more accurate edges. Define the Gaussian convolution kernel The Gaussian function is discretely approximated, and the weighted average of the pixels around the pixel is calculated by selecting the appropriate convolution kernel size and intensity to obtain the final filtered image.

[0033] Step 3-2: Calculate the horizontal and vertical gradients of each pixel in the filtered image using the Canny operator, and finally calculate the magnitude and direction of the gradient of each pixel.

[0034] The gradient direction is perpendicular to the edge direction, and the Canny operator returns the horizontal gradient G of the pixel point (x, y). x and vertical gradient G y , the magnitude G and direction θ (expressed as angle values) of the gradient are:

[0035]

[0036] θ=arctan 2 (G y ,G x )

[0037] The gradient direction is perpendicular to the edge direction and is usually selected from eight different directions: horizontal (left, right), vertical (up, down), and diagonal (upper left, lower left, upper right, lower right). When calculating the gradient of each pixel, two features are obtained: the magnitude and direction of the gradient.

[0038] Step 3-3: Traverse the pixels in the filtered image, perform non-maximum suppression based on the gradient amplitude and direction of each pixel, and obtain the initial edge.

[0039] After obtaining the magnitude and direction of the gradient for each pixel, the pixels in the filtered image are traversed and non-maximum suppression is performed, effectively removing all non-edge points. In practice, the algorithm iterates through each pixel, determining whether the current pixel is the maximum value among the surrounding pixels with the same gradient direction. Based on this determination, the algorithm decides whether to suppress the pixel. For each pixel, if it is a local maximum in the positive or negative gradient direction, it is retained; if not, it is suppressed. After performing non-maximum suppression on all pixels, the initial edge is obtained.

[0040] Step 3-4: Use the double threshold method to judge and mark all pixels included in the initial edge, determine the edge pixels, and then obtain the number of edge pixels.

[0041] Because the edges obtained through the above steps contain some virtual edges, a dual-threshold method is applied to determine the true edges. This method sets two thresholds, a high threshold and a low threshold, and determines the edge type based on the relationship between the gradient value of the current edge pixel and these two thresholds. Specifically, if the gradient value of the current edge pixel is greater than or equal to the high threshold, the current edge pixel is marked as a strong edge; if the gradient value of the current edge pixel is between the high and low thresholds, the current edge pixel is marked as a virtual edge and must be retained; if the gradient value of the current edge pixel is less than or equal to the low threshold, the current edge pixel is suppressed. After dual-threshold screening, the edge pixels are finally determined, and the number of edge pixels C is obtained.

[0042] After edge operator threshold screening, the next calculation is performed for overlapping sub-area pairs whose edge pixel number is greater than the edge threshold. Otherwise, the overlapping sub-area pair is considered unsuitable for offset measurement and the offset is recorded as an invalid value.

[0043] Step 4: Sort and select the optimal overlapping sub-region pairs.

[0044] After the above steps, the number of edge pixels C and the number of local feature points F of the overlapping sub-area pair are obtained. Next, for each overlapping sub-area pair whose number of local feature points F and number of edge pixels C both meet the corresponding thresholds, the product of the number of local feature points and the number of edge pixels is used as the screening index. The overlapping sub-area pairs are sorted in descending order of the screening index. After sorting, the sub-area pair with the largest index is selected as the optimal overlapping sub-area pair.

[0045] Set the edge threshold T c And feature point threshold T f For each overlapping sub-area, when the number of edge pixels C is lower than T c Or the number of feature point pixels F is less than T f When the overlapped sub-area is not suitable for measuring relative image shift, the overlapped sub-area is not measured. When there is no sub-area suitable for measurement in all overlapped sub-areas in a certain image, the overlapped sub-area offset measurement is not performed on the image. When the number of overlapped sub-areas suitable for measurement is not 0, the condition C>T is satisfied. c and F>T f When aligning the overlapping sub-areas, set T according to this step. CF = C*F is used as the basis for determining the degree of suitability of the overlapped sub-area for measurement, and T is obtained by sorting CF The largest overlapping sub-area pair is used as the optimal overlapping sub-area pair for offset measurement, thereby reducing the global error.

[0046] Step 5: Measure the relative image shift.

[0047] Phase correlation registration is performed on the optimal overlapping sub-area pair obtained in step 4. After registration, the relative image shift between the two images is obtained, and then the measured relative image shift is obtained.

[0048] The above-mentioned phase correlation registration method is used to calculate the relative image shift between the left and right optimal overlapping sub-areas in the optimal overlapping sub-area pair, and the specific steps include:

[0049] Step 5-1: First, apply the Hanning window function to the left and right optimal overlapping sub-regions in the screened optimal overlapping sub-region pair to remove the image boundary effect, and obtain images src1 and src2;

[0050] Step 5-2: Calculate the Fourier transform of the two images src1 and src2 respectively. The calculation formula is as follows:

[0051] G a =DFT(src1)

[0052] G b =DFT(src2}

[0053] Among them, DFT is two-dimensional Fourier transform processing, G a is the frequency domain transformation result of image src1, G b is the frequency domain transformation result of image src2;

[0054] Step 5-3: According to G a and G b Calculate the power spectrum R, the calculation formula is as follows:

[0055]

[0056] Step 5-4: Calculate the inverse Fourier transform of the power spectrum R. The calculation formula is as follows:

[0057] r = DFT -1 (R)

[0058] Where r is the phase-matched pulse function, DFT -1 It is a two-dimensional inverse Fourier transform process;

[0059] Step 5-5: Calculate the maximum value of r, i.e., the peak position, and calculate the sub-pixel precision position within a certain range of the window centered on the peak position. Finally, determine the offsets a and b. The calculation formula is as follows:

[0060]

[0061]

[0062] Where a is the relative displacement of the optimal overlapping sub-region pair in the X direction, b is the relative displacement of the optimal overlapping sub-region pair in the Y direction, f(i,j) is where i is the X coordinate in the frequency domain, j is the Y coordinate in the frequency domain, and S*S represents the set of pixels within the calculation window. The calculated relative displacements in the X and Y directions are used as the relative image motion.

[0063] The following example uses the image composed of the left, middle and right detectors as an example. Figure 2 The present invention is described in detail. The present invention provides a method for automatically measuring the relative image shift of the overlapped area of ​​the CCD camera splicing focal plane of an optical remote sensing satellite. The method calculates the relative image shift of the overlapped area of ​​the camera splicing focal plane in real time based on the satellite received image.

[0064] The flow chart for measuring the relative image shift in the overlapping area of ​​a single image is as follows: Figure 2 As shown, the following steps are included:

[0065] (1) Cropping overlapping area images of remote sensing images

[0066] For a single-scene, three-detector stitched remote sensing image, the left, center, and right overlap regions are determined by horizontally shifting the overlap region boundaries by a distance of 150 pixels to the left and right, respectively. Furthermore, N left-center overlap sub-region pairs and N center-right overlap sub-region pairs are cropped vertically with a step size of 1500 pixels to screen for optimal relative image shift measurement sub-regions.

[0067] (2) Calculate the number of local feature points

[0068] The number of Fast local feature points is calculated for each cropped overlapping sub-region pair. The local feature point threshold is used for screening. For overlapping sub-region pairs with a number of local feature points greater than the feature point threshold, the next calculation step is performed. Otherwise, the overlapping sub-region pair is considered unsuitable for offset measurement and the offset is recorded as an invalid value.

[0069] (3) Calculate the number of edge pixels

[0070] The number of Canny edge pixels is calculated for each overlapping sub-region whose number of local feature points meets the feature point threshold requirement. The edge operator threshold is used to screen overlapping sub-region pairs whose number of edge pixels exceeds the edge threshold. Otherwise, the overlapping sub-region pair is considered unsuitable for offset measurement and the offset is recorded as an invalid value.

[0071] (4) Sorting and selecting the optimal overlapping sub-area

[0072] For each overlapping sub-area pair whose number of local feature points and edge pixels meets the threshold requirements, the product of the number of local feature points and the number of edge pixels is used as the screening index. After sorting, the overlapping sub-area pair with the largest index is selected as the optimal overlapping sub-area.

[0073] (5) Measuring relative image shift

[0074] Phase correlation registration is performed on the obtained optimal measurement sub-area pair to obtain the relative displacement between the two images, and then the measured relative image displacement is obtained.

[0075] By automatically measuring the relative image motion of the single-scene overlap area, the time-series monitoring of the relative image motion of the single-star overlap area can be achieved. Figure 3 The method proposed in this invention is used to obtain the time series monitoring results of the relative image shift of the overlap area of ​​the JLGF03D14 optical remote sensing image equipped with three CCD cameras with stitched focal planes, including the X direction ( Figure 3 (a)) and Y direction ( Figure 3 (b) Relative image shift and X direction of the center-right overlap area ( Figure 3 (c)) and Y direction ( Figure 3 The monitoring curve shows that, without adjusting the relative image shift in the external overlap area, the measurement results for different ground scenes captured by the satellite in different attitudes and working conditions can reach the sub-pixel level and are stable within a single pixel.

[0076] The method for automatically measuring the relative image motion of the overlap area of ​​the spliced ​​focal plane of the CCD camera of an optical remote sensing satellite proposed in the present invention can automatically, quickly, and accurately measure the relative image motion of the overlap area, greatly improving the reliability of the measurement results and meeting the requirements of high-precision time-series monitoring of the relative image motion of the overlap area of ​​a large number of satellites. Compared with the classic overlap area relative image motion measurement method, this method automatically determines whether the image is suitable for measurement and selects the appropriate measurement area on the image that is suitable for measurement, thereby obtaining high-precision measurement results, thereby realizing the measurement of the relative image motion of the overlap area of ​​a large number of satellites in the operation of the remote sensing constellation. This measurement method has high accuracy, stable measurement results, is easy to implement, and can be effectively applied in engineering practice.

[0077] The technical features of the above-mentioned embodiments can be combined arbitrarily. In order to make the description concise, not all possible combinations of the technical features in the above-mentioned embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0078] The above-described embodiments merely illustrate several implementations of the present invention, and while their descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the patent. It should be noted that a person skilled in the art would be able to make numerous variations and improvements without departing from the spirit of the present invention, all of which fall within the scope of protection of the present invention. Therefore, the scope of protection of the patent for this invention shall be determined by the appended claims.

Claims

1. A method for automatically measuring the relative image shift of the overlap area of ​​the CCD camera splicing focal plane, characterized in that: The following steps are involved: Step 1: According to the column positions of the adjacent overlapping area boundaries in the known remote sensing image, the preset widths are extended to the left and right to obtain the left overlapping area image and the right overlapping area image respectively. Then, the left overlapping area image and the right overlapping area image are cropped without overlap according to the set cropping step size to obtain the overlapping sub-area pairs to be screened; Step 2: Calculate the number of local feature points for each overlapping sub-region pair based on the Fast operator, and use the local feature point threshold screening method to screen out overlapping sub-region pairs whose local feature points are greater than the feature point threshold; Step 3: Perform edge detection on the overlapping sub-region pairs selected in step 2 based on the Canny operator to obtain the corresponding number of edge pixels, and use the edge operator threshold screening method to screen out overlapping sub-region pairs whose edge pixel number is greater than the edge threshold; Step 4: Sort the overlapping sub-region pairs selected in step 3 using the product of the number of local feature points and the number of edge pixels as the screening index, and select the overlapping sub-region pair corresponding to the maximum value of the screening index as the optimal overlapping sub-region pair; Step 5: Perform phase correlation registration on the optimal overlapping sub-region pair, and obtain the relative image shift between the two images after registration.

2. The method for automatically measuring the relative image shift of the overlap area of ​​the CCD camera splicing focal plane according to claim 1, characterized in that: In step 3, edge detection is performed on each overlapped sub-region pair selected, including the following steps: Step 3-1: Perform Gaussian filtering on the overlapping sub-region pairs to obtain a filtered image; Step 3-2: Calculate the horizontal and vertical gradients of each pixel in the filtered image using the Canny operator, and finally calculate the magnitude and direction of the gradient of each pixel; Step 3-3: Traverse the pixels in the filtered image, perform non-maximum suppression based on the gradient amplitude and direction of each pixel, and obtain the initial edge; Step 3-4: using a double threshold method to judge and mark all pixels included in the initial edge, determine the edge pixels, and then obtain the number of edge pixels.

3. The method for automatically measuring the relative image shift of the overlap area of ​​the CCD camera splicing focal plane according to claim 1, characterized in that: Step 5 includes the following steps: Step 5-1: Applying a Hanning window function to the optimal overlapped sub-region pair to remove image boundary effects, thereby obtaining images src1 and src2; Step 5-2: Calculate the Fourier transform of images src1 and src2 respectively. The calculation formula is as follows: G a =DFT(src1} G a =DFT(src2} Step 5-3: Calculate the power spectrum R based on the Fourier transform results. The calculation formula is as follows: Step 5-4: Calculate the inverse Fourier transform of the power spectrum R to obtain the phase-matched pulse function r, which is calculated as follows: r=DFT -1 (R) Step 5-5: Calculate the peak position of the phase-matched pulse function r, calculate the sub-pixel precision position in the window centered on the peak position, and ultimately determine the offsets a and b. The calculation formula is as follows: Where a is the relative displacement of the optimal overlapping sub-area pair in the X direction, b is the relative displacement of the optimal overlapping sub-area pair in the Y direction, f(i,j) is, i is the X-direction coordinate in the frequency domain, j is the Y-direction coordinate in the frequency domain, and S*S represents the set of pixels in the calculation window.

4. The method for automatically measuring the relative image shift of the overlap area of ​​the CCD camera splicing focal plane according to claim 1, characterized in that: The preset width is 150 pixels, and the cropping step is 1500 pixels.

Images