Large-scale satellite stereoscopic image data processing method, medium and system

Abstract

The invention provides a large-scale satellite three-dimensional image data processing method, medium and system, and belongs to the technical field of satellite three-dimensional image data processing.The method comprises the steps that Gaussian filtering and wavelet noise reduction preprocessing is conducted on an image, and then cooperative computing of a GPU and a CPU is achieved based on a computing task distribution function. Through multi-scale feature extraction and matrix optimization, a stable feature matrix and a variable feature matrix are obtained. A quadtree spatial index is adopted to carry out data partitioning, and SIFT feature extraction and K-means clustering are combined to optimize feature representation. And establishing a stereoscopic image registration model, determining an optimal splicing sequence through a minimum spanning tree algorithm, realizing image splicing by adopting a progressive texture fusion method, and finally generating a high-precision three-dimensional earth surface model. According to the invention, the unification of processing efficiency and precision is realized, and the technical problem that the processing efficiency of large-scale stereoscopic image data is difficult to improve on the premise of ensuring the processing precision in the prior art is solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of satellite stereo image data processing. Specifically, it relates to a method, medium and system for processing large-scale satellite stereo image data. Background Art

[0002] Satellite stereo image data processing technology has wide applications in fields such as geographic information systems, environmental monitoring, and urban planning. Traditional satellite stereo image processing methods mainly rely on a single processing flow and adopt a serialized data processing method, including steps such as image preprocessing, feature extraction, image matching, and 3D reconstruction. In practical applications, common processing methods include image matching based on feature points, stereo matching based on regions, and 3D reconstruction based on semantic segmentation and other technologies. These methods perform well when processing small-scale data and can achieve basic image processing and 3D reconstruction functions.

[0003] However, with the rapid development of satellite remote sensing technology, the amount of stereo image data obtained has increased exponentially, and traditional processing methods are facing serious challenges. First, a single processing flow is difficult to fully utilize the hardware resources of modern computers, resulting in low processing efficiency; second, feature extraction and matching in the process of large-scale data processing are easily affected by noise interference, affecting processing accuracy; third, traditional image stitching methods have cumulative errors when processing a large number of image blocks, and it is difficult to ensure the geometric accuracy and visual continuity of the stitching results. Especially when processing high-resolution stereo image data, these problems are more prominent.

[0004] The existing technology mainly adopts simple data block processing and serial computing methods, lacking effective task scheduling mechanisms and feature optimization strategies. In large-scale data processing, this method cannot achieve the optimal allocation of computing resources and is also difficult to ensure the quality of processing results. In addition, existing image stitching technologies often ignore the stability analysis of features and dynamic parameter adjustment, resulting in problems such as obvious seams and geometric deformation during the stitching process. Therefore, how to improve the processing efficiency of large-scale stereo image data while ensuring processing accuracy has become an urgent technical problem to be solved.

[0005] In summary, there is a technical problem in the existing technology that it is difficult to improve the processing efficiency of large-scale stereo image data while ensuring processing accuracy. Summary of the Invention

[0006] In view of this, the present invention provides a method, medium and system for processing large-scale satellite stereo image data, which can solve the technical problem that it is difficult to improve the processing efficiency of large-scale stereo image data while ensuring processing accuracy in the existing technology.

[0007] The present invention is implemented as follows: A method for processing large-scale satellite stereo image data provided by the first aspect of the present invention includes the following steps: preprocessing satellite stereo image data to obtain preprocessed image data; distributing the preprocessed image data to a GPU processing stream and a CPU processing stream for parallel computing; extracting a multi-dimensional feature matrix of the preprocessed image data; performing feature decomposition using a matrix chain multiplication dynamic programming algorithm; establishing a spatial index to generate data blocks; extracting feature vectors and performing clustering optimization; establishing a registration model and an evaluation function; constructing a mosaic and determining a mosaic order; performing image mosaic and depth information reconstruction; wherein the collaborative computing between the GPU and the CPU is realized through a computing task allocation function, the mosaic order is optimized using a minimum spanning tree algorithm, and dynamic optimization is performed based on a matching quality evaluation function.

[0008] Among them, the step of preprocessing is specifically to perform Gaussian filtering and wavelet denoising on the satellite stereo image data. By setting the Gaussian kernel size to 5 pixels and the standard deviation to 1, selecting the db4 wavelet basis function for 3-layer wavelet decomposition, and using the soft threshold method for denoising, the preprocessed image data is obtained.

[0009] Among them, the step of allocating computing tasks is specifically to establish a computing task allocation function. The computing task allocation function includes a matrix operation identification value and a data stream identification value. When the data scale is greater than 1024×1024 pixels and matrix operations are involved, the matrix operation identification value is set to true and distributed to the GPU processing stream; when the data volume is greater than 100 megabytes and there are continuous data stream operations, the data stream identification value is set to true and distributed to the CPU processing stream.

[0010] Among them, the step of extracting a multi-dimensional feature matrix is specifically to perform 4-layer Haar wavelet decomposition on the preprocessed image data, and extract a geospatial feature matrix, a spectral feature matrix, and a texture feature matrix; use the matrix chain multiplication dynamic programming algorithm to determine the optimal computing order, and perform singular value decomposition to obtain a stable feature matrix and a variable feature matrix.

[0011] Among them, the step of extracting feature vectors and performing clustering optimization is specifically to use the scale-invariant feature transform algorithm to extract feature descriptors, construct an 8-layer scale space, with each group containing 3 scale layers; use the K-means clustering algorithm to cluster the feature descriptors, the number of cluster centers is 0.1 times the total number of feature points, and the number of iterations is 100. Each feature descriptor is replaced with the cluster center point closest in distance to form an optimized feature vector group.

[0012] Among them, the steps of establishing the registration model are specifically to perform feature matching based on the optimized feature vector group, set the distance ratio threshold to 0.8, use the random sample consensus algorithm to remove outliers, set the inlier threshold to 3 pixels, set the number of iterations to 1000 times, and solve the affine transformation parameters by the least squares method to establish the registration model; construct a matching quality evaluation function, and re-execute the construction of the registration model when the matching quality evaluation value is less than 0.75.

[0013] Among them, the steps of constructing the mosaic image are specifically to calculate the overlapping area ratio and the matching quality evaluation value between the image data blocks, and use weighted summation to determine the edge weight value, with the overlapping area ratio weight being 0.3 and the matching quality evaluation value weight being 0.7; use the Kruskal algorithm to solve the minimum spanning tree to determine the mosaic order.

[0014] Among them, the steps of performing image mosaic and depth information reconstruction are specifically to use the minimum cut graph algorithm to calculate the optimal stitching line and set a transition band with a width of 32 pixels; perform local re-mosaic when the local area gray difference value is greater than 8; use the semi-global matching algorithm to calculate the disparity map, set the disparity search range from 0 to 64 pixels, and reconstruct the three-dimensional surface model by the principle of triangulation.

[0015] The second aspect of the present invention provides a computer-readable storage medium, in which program instructions are stored, and when the program instructions run on a computer, they are used to execute the above-mentioned method for processing large-scale satellite stereo image data.

[0016] The third aspect of the present invention provides a large-scale satellite stereo image data processing system, which includes the above-mentioned computer-readable storage medium. The system can be any one of a computer, a server, and a single-chip microcomputer. The computer-readable storage medium is set inside the system, and a microprocessor for executing the program instructions stored in the computer-readable storage medium is set inside the system.

[0017] Compared with the prior art, the present invention provides a method, medium, and system for processing large-scale satellite stereo image data. The method for processing large-scale satellite stereo image data proposed by the present invention realizes the collaborative computing of the GPU and the CPU by establishing a computing task allocation function. This method adopts a multi-scale feature extraction and matrix optimization strategy, combined with a quadtree spatial index structure, which significantly improves the efficiency and accuracy of data processing.

[0018] In practical applications, the method of the present invention effectively suppresses image noise through the combined preprocessing of Gaussian filtering and wavelet denoising; improves the reliability of feature matching through the stability decomposition of features and dynamic parameter adjustment; and adopts an optimal stitching sequence generation strategy based on a graph structure to reduce the influence of cumulative errors. At the same time, the application of the progressive texture fusion method ensures the visual naturalness of the stitching result. Especially when dealing with stereo image data on the TB scale, the method of the present invention can still maintain stable processing performance and accuracy.

[0019] The present invention successfully solves the problems of efficiency and accuracy in the processing of large-scale satellite stereo image data, which is mainly due to its innovative computing resource scheduling strategy and feature optimization processing flow. By allocating processing tasks with different characteristics to the most suitable computing units and adopting a multi-level feature extraction and optimization strategy, the unity of processing efficiency and accuracy is achieved, and the technical problem that it is difficult for the existing technology to improve the processing efficiency of large-scale stereo image data while ensuring processing accuracy is solved. Brief Description of the Drawings

[0020] Figure 1 It is a flowchart of the method provided by the present invention. Detailed Embodiments

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention.

[0022] As Figure 1 shown, it is a flowchart of a method for processing large-scale satellite stereo image data provided by the first aspect of the present invention. The method includes the following steps: S01. Receive stereo image data collected by multiple satellites, and perform Gaussian filtering and wavelet denoising processing on the stereo image data to obtain preprocessed image data; S02. Establish a computing task allocation function, where the computing task allocation function includes a matrix operation identification value and a data stream identification value. According to the computing task allocation function, divide the tasks of the preprocessed image data, and allocate the image block operation tasks with the matrix operation identification value being true to the GPU processing stream for parallel acceleration calculation, and allocate the continuous data processing tasks with the data stream identification value being true to the CPU processing stream for pipeline calculation; S03. Perform multi-scale wavelet decomposition on the preprocessed image data, extract three-dimensional spatial features, multi-spectral band features, and gray-level co-occurrence matrix features, and obtain a geospatial feature matrix, a spectral feature matrix, and a texture feature matrix; S04. Use the matrix chain multiplication dynamic programming algorithm to calculate the minimum number of multiplication operations for the geospatial feature matrix, spectral feature matrix, and texture feature matrix, and perform singular value decomposition according to the minimum number of multiplication operations to obtain a stable feature matrix and a variable feature matrix; S05. Use a quadtree structure to partition the space of the stable feature matrix, establish a spatial hierarchical index, and generate an image data block index table; S06. Perform block processing on the preprocessed image data according to the image data block index table to generate a number of image data blocks of regular size; S07. Use the scale-invariant feature transform algorithm to extract the feature descriptors of each image data block, and use the K-means clustering algorithm to cluster the feature descriptors to obtain the feature clustering center points; S08. Replace each feature descriptor with the feature clustering center point closest in distance to form an optimized feature vector group; S09. Calculate the image registration parameters based on the optimized feature vector group, including the rotation angle parameter, scaling ratio parameter, and translation vector parameter, and establish an affine transformation model as the stereo image registration model; S10. Calculate the correlation coefficient between the stereo image registration model and the variable feature matrix, and establish a registration parameter correction function; S11. Establish a matching quality evaluation function based on the root mean square error criterion, calculate the gradient direction consistency value and gray correlation coefficient between image blocks to obtain the matching quality evaluation value. When the matching quality evaluation value is less than the preset matching threshold, return to execute steps S09 to S11; S12. Construct an image data block mosaic graph, where the graph nodes represent the image data blocks, and the graph edge weight value is the weighted sum of the overlapping area between adjacent image data blocks and the matching quality evaluation value; S13. Apply the minimum spanning tree algorithm to the image data block mosaic graph to determine the optimal splicing order of the image data blocks and generate a splicing operation sequence; S14. Splice the image data blocks one by one according to the splicing operation sequence, and use the progressive texture fusion method to process the overlapping area to obtain a spliced image; S15. Calculate the gray difference value at the boundary of adjacent image blocks in the spliced image. When the gray difference value is greater than the preset splicing threshold, perform local re-splicing optimization on the boundary area; S16. Use a stereo matching algorithm to extract the disparity information of the spliced image, calculate the three-dimensional space coordinates, construct a depth information map, and use the triangular mesh reconstruction algorithm to generate a three-dimensional surface model.

[0023] The specific implementation manners of the above steps are described in detail below. The specific implementation manner of step S01 is as follows: For the received satellite stereo image data, first, Gaussian filtering is used to smooth the image. By setting the Gaussian kernel size to 5 pixels and the standard deviation to 1, the noise in the image is preliminarily suppressed. Gaussian filtering can effectively remove Gaussian noise while maintaining the edge features of the image. Then, the wavelet denoising algorithm is adopted. The db4 wavelet basis function is selected, and the image is decomposed into 3 layers by wavelet. The wavelet coefficients are processed by the soft threshold method, and the threshold is set to 3 times the standard deviation. The denoised image is obtained through wavelet reconstruction. The purpose of this step is to eliminate various noises generated during the satellite imaging process and improve the accuracy of subsequent processing. Through the combination of Gaussian filtering and wavelet denoising, the detailed information of the image can be better maintained while effectively suppressing noise interference. The result of this step is recorded as the preprocessed image data.

[0024] The specific implementation manner of step S02 is as follows: First, a computing task allocation function is constructed. This function contains two main judgment criteria: the matrix operation identification value and the data stream identification value. For the input preprocessed image data, its data scale and computing complexity are calculated. When the data scale is greater than 1024×1024 pixels and matrix operations are involved, the matrix operation identification value is set to true; when there are continuous data stream operations during the data processing and the data volume is greater than 100 megabytes, the data stream identification value is set to true. According to these two identification values, the computing tasks are allocated to different processing units: the tasks with the matrix operation identification value being true are allocated to the GPU to accelerate the matrix operations using its parallel computing ability; the tasks with the data stream identification value being true are allocated to the CPU to utilize its efficient data pipeline processing ability. The purpose of this step is to achieve the optimal configuration of computing resources and improve the overall computing efficiency.

[0025] The specific implementation manner of step S03 is as follows: The discrete wavelet transform is used to perform multi-scale decomposition on the preprocessed image data. The Haar wavelet is selected as the basis function for 4-layer decomposition, and the low-frequency approximation component and high-frequency detail components are extracted respectively. In terms of spatial feature extraction, the horizontal and vertical gradients of the image are calculated using the gradient operator, and a geospatial feature matrix is constructed by combining the gradient magnitude and direction information. For spectral feature extraction, the principal component analysis method is used to reduce the dimension of the multi-spectral band data, and the principal components with a contribution rate greater than 95% are retained to form a spectral feature matrix. In terms of texture feature extraction, the gray-level co-occurrence matrix is calculated. The distance parameter is set to 1, and the direction angles are 0 degrees, 45 degrees, 90 degrees, and 135 degrees. Statistical features such as energy, entropy, contrast, and correlation are extracted to construct a texture feature matrix. The purpose of this step is to comprehensively extract the multi-dimensional feature information of the image and provide a reliable feature basis for subsequent processing.

[0026] The specific implementation of step S04 is as follows: First, the matrix chain multiplication dynamic programming algorithm is applied to determine the optimal matrix multiplication order. By calculating the computational costs of different multiplication orders, the order with the minimum total computational amount is selected. In specific implementation, a cost matrix and a split point matrix are constructed. The matrix dimension is the number of feature matrices. The matrices are filled through the dynamic programming method, and finally the optimal calculation order is obtained. Then, according to the optimal order, singular value decomposition is performed on each feature matrix. The singular value threshold is set to 0.1 times the maximum singular value. The singular values greater than the threshold and their corresponding singular vectors form the stable feature matrix, and the rest form the variable feature matrix. The purpose of this step is to optimize the matrix calculation efficiency and at the same time achieve the stable decomposition of features.

[0027] The specific implementation of step S05 is as follows: A quadtree structure is used to construct the spatial index for the stable feature matrix. First, the entire feature space is divided into four equal sub-regions. The feature density is calculated for each sub-region. When the feature density is greater than the preset threshold of 0.8, the sub-region is further divided into four equal parts until the preset minimum region size of 32×32 pixels is reached or the feature density is less than the threshold. During the construction process, a unique index identifier is assigned to each node, and its spatial range and the included feature information are recorded to form a hierarchical spatial index structure. Finally, an image data block index table is generated, which contains the spatial position, size, and corresponding feature information of each data block. The purpose of this step is to establish an efficient spatial retrieval mechanism to provide index support for subsequent block processing.

[0028] The specific implementation of step S06 is as follows: According to the spatial division information in the image data block index table, the preprocessed image data is block-processed. First, the basic block size is determined to be 256×256 pixels. Considering the overlapping requirements for subsequent processing, an overlapping region is set during block division, and the size of the overlapping region is 64 pixels. For each data block, its original image data, position information, and adjacency relationship are saved. During the block division process, for the regions at the image edges that are less than the basic block size, mirror filling is used to expand them to the standard size. The purpose of this step is to divide the large-scale image data into small blocks suitable for parallel processing, and at the same time, through the setting of the overlapping region, provide a reliable transition region for subsequent feature matching and image stitching.

[0029] The specific implementation of step S07 is as follows: Use the Scale-Invariant Feature Transform (SIFT) algorithm to extract the feature points of each image data block, including the position, scale, and orientation information of the feature points. In the feature point detection stage, construct a scale space with 8 scale levels, where each group contains 3 scale levels. Calculate the response of the Difference of Gaussian (DoG) operator, and through extreme point detection and edge response verification, filter out the stable feature points. For each feature point, calculate its main direction and generate a 128-dimensional feature descriptor. Then use the K-means clustering algorithm to cluster all the feature descriptors, set the number of cluster centers to 0.1 times the total number of feature points, use the Euclidean distance as the similarity metric, and set the number of iterations to 100 until the change in the position of the cluster centers is less than the preset threshold of 0.01 or the maximum number of iterations is reached. The purpose of this step is to extract the significant features of the image data block and achieve dimensionality reduction and normalization of the features through clustering.

[0030] The specific implementation of step S08 is as follows: For the feature descriptors of each image data block, calculate the Euclidean distance between them and all the cluster center points, and select the cluster center point with the smallest distance as the alternative representation of this feature descriptor. In this way, map the original 128-dimensional feature descriptors to the cluster center space to achieve the simplification and unification of the feature representation. During the replacement process, simultaneously record the distance values between each feature descriptor and the corresponding cluster center, and these distance values will be used for the subsequent reliability evaluation of feature matching. The finally obtained optimized feature vector group has a more compact representation form and better discriminability. The purpose of this step is to reduce the redundancy of the feature representation and improve the efficiency and stability of the subsequent matching process.

[0031] The specific implementation of step S09 is as follows: Based on the optimized feature vector group, construct a stereo image registration model. First, use the nearest neighbor distance ratio method for feature matching, set the distance ratio threshold to 0.8, and filter out the reliable matching point pairs. Use the Random Sample Consensus (RANSAC) algorithm for outlier rejection, set the inlier threshold to 3 pixels, and the number of iterations to 1000 times to obtain an accurate set of matching points. Based on the set of matching points, calculate the affine transformation parameters, including the rotation angle, scaling ratio, and translation vector. Solve the transformation parameters by the least squares method to establish the geometric transformation relationship between the image pairs. The purpose of this step is to establish an accurate spatial correspondence relationship between the image blocks and provide a geometric transformation basis for the subsequent image stitching.

[0032] The specific implementation of step S10 is as follows: Calculate the relationship between the stereo image registration model and the change feature matrix. First, for each feature vector in the change feature matrix, use the registration model for transformation and calculate the difference in feature response before and after the transformation. The normalized cross-correlation method is used to calculate the correlation coefficient. On this basis, a registration parameter correction function is established. This function uses the weighted least squares method to determine the correction coefficient, and the weight value is determined according to the reliability of the feature response. When the correlation coefficient is lower than 0.7, the adaptive adjustment of the registration parameters is triggered. The purpose of this step is to establish a dynamic correction mechanism for the registration parameters and improve the adaptability and robustness of the registration model.

[0033] The specific implementation of step S11 is as follows: Construct a matching quality evaluation function based on the root mean square error criterion. This function comprehensively considers two aspects: the consistency of gradient directions and the gray-level correlation between image patches. The consistency of gradient directions is evaluated by calculating the difference in gradient directions at corresponding points, and the direction difference threshold is set to 30 degrees. The gray-level correlation is calculated using the normalized cross-correlation coefficient. Sub-blocks of size 32×32 pixels are selected in the overlapping area for correlation analysis. The weighted sum of the two evaluation indicators is used as the final matching quality evaluation value, with the weight ratio of the gradient direction consistency being 0.6 and the gray-level correlation being 0.4. When the matching quality evaluation value is less than the preset matching threshold of 0.75, it is necessary to return and re-execute the construction and parameter optimization of the registration model. The purpose of this step is to ensure the accuracy and reliability of image matching.

[0034] The specific implementation of step S12 is as follows: Construct a stitching graph representing the connection relationship between image data blocks. The nodes in the graph represent each image data block, and the edges represent the connection relationship between adjacent data blocks. For each pair of adjacent image data blocks, calculate the area ratio of their overlapping region and the matching quality evaluation value, and use the weighted summation method to determine the weight value of the edge, where the weight of the overlapping area ratio is 0.3 and the weight of the matching quality evaluation value is 0.7. When the area ratio of the overlapping region is less than 0.2 or the matching quality evaluation value is less than 0.6, no connection edge is established between the corresponding nodes. The purpose of this step is to establish the topological relationship between image data blocks and provide a graph structure support for the generation of the optimal stitching sequence.

[0035] The specific implementation of step S13 is as follows: Convert the image data block stitching graph into a minimum spanning tree problem and use the Kruskal algorithm to solve the minimum spanning tree. First, sort the weight values of all edges in the graph. Starting from the edge with the smallest weight, sequentially determine whether adding this edge will form a loop. If it does not form a loop, add this edge to the minimum spanning tree. During the process of adding edges, use the union-find data structure to maintain the connectivity between nodes to ensure that the finally generated is a tree structure. According to the direction of the edges of the minimum spanning tree, determine the stitching order of the image data blocks and generate a stitching operation sequence. The purpose of this step is to optimize the stitching order and minimize the cumulative error.

[0036] The specific implementation of step S14 is as follows: Perform the stitching of the image data blocks according to the stitching operation sequence, and use the progressive texture fusion method to process the overlapping areas. In the overlapping areas, first use the minimum graph cut algorithm to calculate the optimal stitching line. The cost function of the stitching line takes into account two factors: pixel value difference and gradient difference. Then, set a transition zone on both sides of the stitching line. The width of the transition zone is set to 32 pixels, and within the transition zone, a weighted fusion method is used for smooth transition, and the weight changes linearly as the distance of the pixel from the stitching line increases. For each stitching operation, update the global coordinate system information and the stitching status flag of the image. The purpose of this step is to achieve seamless stitching of the image and eliminate stitching traces.

[0037] The specific implementation of step S15 is as follows: Calculate the gray difference value at the boundary of adjacent image blocks in the stitched image. Use the sliding window method with a window size of 16×16 pixels to calculate the root mean square difference of the local area point by point along the boundary. When the gray difference value of a certain local area is greater than the preset stitching threshold of 8, mark this area as the area to be optimized. For the area to be optimized, recalculate the optimal stitching line of this area, expand the range of the transition zone to 64 pixels, and use a non-linear weight function for gray value fusion to achieve optimized stitching of the local area. The purpose of this step is to further improve the stitching effect and eliminate local discontinuities.

[0038] The specific implementation of step S16 is as follows: Use the stereo matching algorithm to extract the disparity information of the stitched image. First, perform epipolar rectification on the image, use the semi-global matching algorithm to calculate the disparity map, set the disparity search range from 0 to 64 pixels, the number of cost aggregation paths to 8, and the smoothness parameter to 0.8. According to the disparity information and camera parameters, calculate the three-dimensional space coordinates through the principle of triangulation to obtain the depth information map, and then use the triangular mesh reconstruction algorithm to generate the three-dimensional surface model.

[0039] The functions or calculation processes involved in the present invention are described in detail as follows: The Gaussian filter kernel function is expressed as follows: ; In the formula, is the standard deviation of the Gaussian kernel, and its value is 1; is the pixel coordinate position; is the base of the natural logarithm.

[0040] The wavelet denoising threshold function is expressed as follows: ; In the formula, is the wavelet coefficient; is the soft threshold, ; is the total number of image pixels; is the noise standard deviation; is the indicator function.

[0041] The calculation task allocation function is expressed as follows: ; In the formula, is the matrix operation complexity index, ranging from 0 to 1; is the data traffic index, ranging from 0 to 1; is the calculation dependency index, ranging from 0 to 1; is the weight coefficient, taking values 0.5, 0.3, and 0.2 respectively. When , the task is assigned to the GPU processing stream; otherwise, it is assigned to the CPU processing stream.

[0042] The geospatial feature matrix is expressed as follows: ; In the formula, is the row column pixel gradient magnitude.

[0043] The spectral feature matrix is expressed as follows: ; In the formula, is the original multispectral data; is the transformation matrix; is the residual matrix.

[0044] The singular value decomposition is expressed as follows: ; In the formula, is the input feature matrix; is the left and right singular vector matrix; is the singular value diagonal matrix; is the residual matrix.

[0045] The K-means clustering objective function is expressed as follows: ; In the formula, is the number of clusters; is the th cluster; is the th cluster center; is the feature vector.

[0046] The affine transformation matrix is expressed as follows: ; In the formula, is the transformation coefficient; is the translation component.

[0047] The registration parameter correction function is expressed as follows: ; In the formula, is the parameter correction amount; is the current matching quality; is the target matching quality; is the gain matrix; is the historical correction amount; is the historical weight coefficient, and its value range is from 0 to 1.

[0048] The matching quality evaluation function is expressed as follows: ; In the formula, is the image block to be matched; is the gradient operator; is the mean value; is the standard deviation; is the weight coefficient.

[0049] The edge weight calculation function is expressed as follows: ; In the formula, is the overlapping area; is the area of the image block; is the matching quality; is the depth continuity metric; is the weight coefficient.

[0050] The disparity calculation cost function is expressed as follows: ; In the formula, is the pixel matching cost; is the smoothness constraint term; is the temporal consistency constraint term; is the balance coefficient.

[0051] The derivation and establishment processes of multiple functions are described in detail below.

[0052] 1. Derivation and establishment process of the Gaussian filter kernel function: The Gaussian filter kernel function is derived based on the normal distribution in probability theory. The derivation steps are as follows: First, construct the basic form of the two-dimensional normal distribution: ; To ensure that the sum of the filter kernel weights is 1, determine the coefficient A: ; Finally, a standardized Gaussian filter kernel function is obtained: .

[0053] Parameter Through experimental verification, when a better balance can be achieved between denoising and preserving edge details.

[0054] 2. The derivation and establishment process of the wavelet denoising threshold function: Based on the Bayesian estimation theory, the derivation steps are as follows: Establish a noise model: ; where is the noisy wavelet coefficient, is the true signal, is the Gaussian noise; Minimize the mean square error: ; The soft threshold function is obtained by solving: .

[0055] Threshold is obtained by maximum likelihood estimation.

[0056] 3. The derivation and establishment process of the computing task allocation function: First, construct the matrix operation complexity index: ; In the formula, is the number of rows of the th matrix; is the number of columns of the th matrix; is the matrix operation type coefficient, taking 1 for addition and subtraction operations, 2 for multiplication and division operations, and 3 for inverse operations; is the reference operation amount, with a value of floating-point operations; is the number of matrices to be processed.

[0057] Construct the data traffic index: ; In the formula, is the size of the th data block; is the data access frequency coefficient, taking 1 for continuous access and 2 for random access; is the system memory capacity; is the number of data blocks.

[0058] Construct the data dependence degree index: ; In the formula, is an element of the dependence relationship matrix. If data block depends on data block then take 1, otherwise take 0; is the total number of data blocks.

[0059] The final task allocation function is expressed as: ; In the formula, is the weight coefficient, which is determined through the following steps: (1) Construct the judgment matrix: ; (2) Calculate the eigenvalue and eigenvector; (3) Normalize to obtain the weight coefficient: .

[0060] When , the task is assigned to the GPU processing stream; otherwise, it is assigned to the CPU processing stream.

[0061] The derivation and establishment process of the feature matrix: The spatial feature matrix is based on multi-scale gradient analysis: Calculate the horizontal and vertical gradients: ; Calculate the gradient amplitude: ; Form the feature matrix S.

[0062] The spectral feature matrix adopts the principal component analysis method: calculate the covariance matrix; eigenvalue decomposition; select the principal components to construct the transformation matrix W.

[0063] The derivation and establishment process of the registration parameter correction function: Designed based on the feedback control theory: (1) Establish the error model: ; (2) Introduce proportional-integral control: ; (3) Consider the influence of the historical correction amount: .

[0064] The parameter K is obtained through system identification, is determined through cross-validation.

[0065] 6. Derivation and establishment process of the matching quality evaluation function: Comprehensively consider gradient consistency and gray correlation: (1)Gradient consistency measurement: ; (2)Gray correlation measurement: ; (3)Obtain the evaluation function through weighted combination.

[0066] The weight coefficient is obtained through regression analysis of a large amount of experimental data.

[0067] The second aspect of the present invention provides a computer-readable storage medium, in which program instructions are stored. When the program instructions run on a computer, they are used to execute the above-mentioned method for processing large-scale satellite stereo image data.

[0068] The third aspect of the present invention provides a large-scale satellite stereo image data processing system, which includes the above-mentioned computer-readable storage medium. The system can be any one of a computer, a server, and a single-chip microcomputer. The computer-readable storage medium is set inside the system, and a microprocessor for executing the program instructions stored in the computer-readable storage medium is set inside the system.

[0069] Specifically, the principle of the present invention is: The technical principle of the present invention is based on the analysis of the characteristics of computing tasks and the optimization of multi-level features. First, by analyzing the computing characteristics of different processing tasks, a task allocation mechanism based on matrix operation identification values and data flow identification values is established. Matrix operation-intensive tasks are suitable for parallel processing on the GPU, while data flow-intensive tasks are more suitable for pipeline processing on the CPU. This allocation mechanism makes full use of the advantages of heterogeneous computing platforms.

[0070] In terms of feature processing, the present invention uses multi-scale wavelet decomposition to extract three-dimensional spatial features, multi-spectral band features, and texture features, which describe the characteristics of the image from different angles. The calculation order of features is optimized through the matrix chain multiplication dynamic programming algorithm, reducing the computational complexity. The stability decomposition of features ensures the retention of key features, and at the same time, the dynamic parameter adjustment is carried out by changing the feature matrix, improving the adaptability of the system.

[0071] In terms of image stitching, the present invention constructs a spatial index structure based on a quadtree to achieve efficient spatial retrieval. The optimal stitching sequence is determined through the minimum spanning tree algorithm, combined with the progressive texture fusion method, ensuring the geometric accuracy and visual quality of the stitching result. This multi-level optimization strategy ensures the efficiency and reliability of the entire processing flow.

[0072] A specific Embodiment 1 of the present invention is provided below. The specific implementation of each step in this Embodiment 1 is described in detail as follows.

[0073] The specific implementation of step S01 is as follows: For the received satellite stereo image data, first perform smoothing processing on the image using Gaussian filtering. The Gaussian filtering kernel function is expressed as follows: ; In the formula, is the standard deviation of the Gaussian kernel, with a value of 1; is the pixel coordinate position; is the base of the natural logarithm. The implementation of Gaussian filtering adopts the separable kernel method. First, perform horizontal filtering, and then perform vertical filtering. The filtered image is expressed as: ; In the formula, is the filtered image; is the original image; represents the convolution operation. Then, adopt the wavelet denoising algorithm, select the db4 wavelet basis function, perform 3-layer wavelet decomposition on the image, and process the wavelet coefficients using the soft threshold method. The threshold function is expressed as follows: ; In the formula, is the wavelet coefficient; is the soft threshold, ; is the total number of image pixels; is the noise standard deviation; is the indicator function. The purpose of this step is to eliminate various noises generated during satellite imaging and improve the accuracy of subsequent processing.

[0074] The specific implementation of step S02 is as follows: Construct a computing task allocation function to achieve optimal allocation of computing resources by evaluating matrix operation complexity indicators, data traffic indicators, and data dependence degree indicators. The matrix operation complexity indicator is expressed as follows: ; In the formula, is the number of rows of the th matrix; is the number of columns of the th matrix; is the matrix operation type coefficient. Take 1 for addition and subtraction operations, 2 for multiplication and division operations, and 3 for inverse operations; is the reference operation amount, with a value of floating-point operations; is the number of matrices to be processed. The data traffic indicator is expressed as follows: ; In the formula, is the size of the th data block; is the data access frequency coefficient, taking 1 for consecutive access and 2 for random access; is the system memory capacity; is the number of data blocks. The data dependence degree index is expressed as follows: ; In the formula, is an element of the dependence relationship matrix. If data block depends on data block then it takes 1, otherwise it takes 0; is the total number of data blocks. The final task allocation function is expressed as: ; In the formula, is the weight coefficient, determined by the analytic hierarchy process. When , the task is assigned to the GPU processing stream, otherwise it is assigned to the CPU processing stream.

[0075] The specific implementation of step S03 is: perform multi-dimensional feature extraction on the preprocessed image data. First, use discrete wavelet transform for multi-scale decomposition. Select Haar wavelet as the basis function and perform 4-layer decomposition. The coefficient matrix after each layer of decomposition is expressed as follows: ; In the formula, is the number of decomposition layers; is the coefficient serial number. In terms of spatial feature extraction, use the Sobel operator to calculate the gradient information of the image: ; ; In the formula, are the gradients in the horizontal and vertical directions respectively; is the input image. The gradient magnitude and direction are calculated as follows: ; ; Construct a geospatial feature matrix: .

[0076] The specific implementation of step S04 is: use the matrix chain multiplication dynamic programming algorithm to determine the optimal calculation order of the feature matrix and perform singular value decomposition to achieve feature separation. First, construct a cost matrix for the matrix multiplication operation of the feature matrix. For the matrix sequence , where matrix The dimension of is ; In the formula, represents the minimum number of scalar multiplications required for calculating . After determining the optimal calculation order, perform singular value decomposition on each feature matrix: ; In the formula, is the input feature matrix; is the left and right singular vector matrix; is the singular value diagonal matrix; is the residual matrix. The extraction of the stable feature matrix is based on singular value threshold selection: ; In the formula, is the threshold coefficient, with a value of 0.1; is the maximum singular value.

[0077] The specific implementation of step S05 is: Based on the quadtree structure, construct a spatial index for the stable feature matrix. The feature density calculation formula during the spatial partitioning process is: ; In the formula, is the feature density of block ; is the weight of the th feature point in the block; is the area of the block; is the number of feature points in the block. The feature point weight calculation formula is: ; In the formula, is the distance from the feature point to the center of the block; is the scale parameter, with a value of 0.2 times the side length of the block. The block splitting condition is: ; In the formula, is the density threshold, with a value of 0.8; is the minimum block area, with a value of pixels.

[0078] The specific implementation of step S06 is: According to the spatial index, partition the preprocessed image data. The block size calculation formula is: ; In the formula, is the block side length; is the size of the GPU video memory; is the number of image channels. The size of the overlapping area between adjacent blocks is calculated as follows: ; In the formula, is the width of the overlapping area. For the area at the edge of the image that is less than a complete block, mirror filling is adopted: ; In the formula, are the width and height of the image respectively.

[0079] The specific implementation of step S07 is: using the Scale-Invariant Feature Transform (SIFT) algorithm to extract feature points. First, construct a scale space: ; In the formula, is a scale-variable Gaussian kernel; is the scale parameter. Calculate the Difference of Gaussian (DoG) images: ; In the formula, is the scale factor, and its value is . The local extreme detection conditions for feature points are: ; In the formula, is the detection threshold, and its value is 0.03; is the 26-neighborhood of the feature point. The main direction of the feature point is calculated based on the gradient histogram: ; .

[0080] The specific implementation of step S08 is: performing clustering optimization on the extracted feature descriptors. First, calculate the distance matrix between the feature descriptors: ; In the formula, is the -th component of the -th feature descriptor; is the descriptor dimension, and its value is 128. Use the K-means clustering algorithm for feature classification. The clustering objective function is: ; In the formula, is the number of clusters, and its value is 0.1 times the total number of feature points; is the -th cluster; is the The number of cluster centers. The update formula for the cluster centers is as follows: ; In the formula, is the number of samples in the th cluster. The clustering termination condition is: ; In the formula, is the convergence threshold, with a value of 0.01.

[0081] The specific implementation of step S09 is: Based on the optimized feature vector group, a stereo image registration model is constructed. First, feature matching is performed, and the nearest neighbor distance ratio method is used for preliminary screening: ; In the formula, are the nearest neighbor and the second nearest neighbor distances respectively; is the distance ratio threshold, with a value of 0.8. Solve the affine transformation matrix of the matching point pairs: ; Among them, the transformation parameters are solved by the least squares method: ; In the formula, is the corresponding point pair; is the number of matching point pairs.

[0082] The specific implementation of step S10 is: Establish a registration parameter correction function, and the correction amount of the registration parameters is calculated as follows: ; In the formula, is the current matching quality evaluation value; is the target matching quality value, with a value of 0.9; is the gain matrix, obtained through system identification: ; is the historical correction amount, and its update formula is: ; In the formula, is the smoothing coefficient, with a value of 0.7.

[0083] The specific implementation of step S11 is: Construct a matching quality evaluation function, considering gradient consistency and gray correlation comprehensively: ; In the formula, They are the weights of gradient consistency and gray correlation, with values of 0.6 and 0.4 respectively. The evaluation area adopts a sliding window method: ; In the formula, is the window radius, with a value of 16 pixels.

[0084] The specific implementation of step S12 is: construct a mosaic of image data blocks, and the edge weight calculation function of the graph is: ; In the formula, is the overlapping area; is the area of the image block; is the matching quality; is the depth continuity metric: ; In the formula, is the number of pixels in the overlapping area; is the depth value of the corresponding point.

[0085] The specific implementation of step S13 is: based on the constructed mosaic of image data blocks, use the Kruskal algorithm to solve the minimum spanning tree. First, sort the edges in the graph in ascending order of weight, and the sorting criterion for the edges is: ; In the formula, is the edge in the graph; is the corresponding edge weight. Use the union-find data structure to maintain node connectivity, and the search function for the representative element of the set is: ; In the formula, is the parent node of node . The optimization function for the set merge operation is: ; In the formula, is the height of the tree where node is located.

[0086] The specific implementation of step S14 is: perform image stitching according to the stitching order determined by the minimum spanning tree, and use the progressive texture fusion method to process the overlapping area. First, use the minimum cut graph algorithm to calculate the optimal seam line, and the graph cut energy function is defined as: ; In the formula, is the data term, indicating the cost of pixel being marked as ; It is a smoothing term, representing the penalty for the difference in adjacent pixel labels: ; ; In the formula, is the smoothing coefficient, with a value of 10; is the truncation threshold, with a value of 30.

[0087] The specific implementation of step S15 is: locally optimize the spliced image and calculate the gray difference value of the boundary region: ; In the formula, is the local window size, with a value of 16. When the gray difference value is greater than the threshold, adopt adaptive weight fusion: ; The weight function is defined as: ; In the formula, is the distance from the pixel to the suture line; is the smoothing parameter, with a value of 32.

[0088] The specific implementation of step S16 is: adopt the semi-global matching algorithm to extract the disparity information and construct the cost function for disparity calculation: ; In the formula, is the pixel matching cost: ; is the smoothness constraint term: ; is the temporal consistency constraint term: ; In the formula, are the smoothness penalty coefficients, with values of 0.8 and 1.6 respectively; are the temporal constraint parameters, with values of 0.5 and 2 respectively. Finally, calculate the three-dimensional coordinates through triangulation: ; ; ; In the formula, is the camera focal length; is the baseline length; is the disparity value; is the principal point coordinate.

[0089] The parameters in all the above steps are determined through a large number of experiments. While ensuring the robustness of the algorithm, the balance between computational efficiency and accuracy is ensured. The entire processing flow fully considers the characteristics of satellite images and adopts a multi-level optimization strategy to achieve efficient and reliable stereo image processing.

[0090] To better understand and implement the present invention, Example 2 of a specific application scenario of the present invention is provided below: When a certain research institute conducts a project on urban 3D modeling, it needs to process high-resolution satellite stereo image data covering an area of about 100 square kilometers. This batch of data includes 5 main images and 5 auxiliary images, with a spatial resolution of 0.5 meters. The size of each image is 25000×20000 pixels, and the total data volume is about 60GB. The research team uses the method of the present invention to process this batch of data.

[0091] First, the research team preprocesses the received original image data. Gaussian filtering is performed using a Gaussian filter kernel with a size of 5×5 pixels. The Gaussian kernel function is: , where . Through Gaussian filtering, the random noise in the image is effectively suppressed, and the image signal-to-noise ratio is increased from the original 26dB to 31dB. Then, the db4 wavelet basis function is used to perform 3-layer wavelet decomposition on the image, and the soft threshold function: is used for coefficient processing, where . After wavelet reconstruction, the preprocessed image after noise reduction is obtained.

[0092] The research team configures 2 NVIDIA A100 GPUs and 1 dual-way Intel Xeon processor in the processing system. According to the computing task allocation function: , the feature extraction and matching tasks involving matrix operations are assigned to the GPU processing stream, and the data reading and image stitching tasks are assigned to the CPU processing stream. The specific task allocation is shown in Table 1.

[0093] Table 1: Computing Task Allocation Table

[0094] The research team performs multi-scale feature extraction on the preprocessed images. The geospatial feature matrix S is constructed by calculating the gradient magnitude, and the matrix dimension is 25000×20000. Principal component analysis is performed on the multi-spectral data to obtain the spectral feature matrix P, and the first 3 principal components with a contribution rate of 97.3% are retained. The energy, entropy, contrast, and correlation features of the gray-level co-occurrence matrix are calculated to construct the texture feature matrix.

[0095] The matrix chain multiplication dynamic programming algorithm is used to determine the optimal calculation order, and the singular value decomposition is performed on the feature matrix: Set the singular value threshold to 0.1 times the maximum singular value. The part greater than the threshold constitutes the stable feature matrix, and the rest constitutes the variable feature matrix. The stable feature matrix contains the main structural information of the image, with a dimension of 25000×1000.

[0096] The research team uses a quadtree structure to perform spatial partitioning on the stable feature matrix. Set the minimum area size to 32×32 pixels and the feature density threshold to 0.8. The generated spatial index tree has a depth of 8 layers and contains approximately 12000 leaf nodes. Based on the index structure, the image data is divided into 2048 basic data blocks, each with a size of 256×256 pixels, and a 64-pixel overlapping area is set between adjacent blocks.

[0097] Extract SIFT feature points for each image data block. On average, about 200 feature points are extracted per block. Use the K-means clustering algorithm to cluster the feature descriptors. The clustering objective function is: , and the number of clusters is set to 20. Replace the original feature descriptors with the nearest cluster centers to form an optimized feature vector group.

[0098] Construct an affine transformation model based on the optimized feature vector group. The transformation matrix is: . Screen the matching point pairs through the RANSAC algorithm, and set the inlier threshold to 3 pixels. Calculate the registration parameter correction amount: , where , and iterate and optimize until the matching quality evaluation value is greater than 0.75.

[0099] The research team constructs a mosaic of image data blocks and calculates the edge weights: , where , , . Apply the Kruskal algorithm to solve the minimum spanning tree and determine the mosaic order. A 32-pixel-wide progressive fusion transition zone is used in the overlapping area to smooth the mosaic boundary.

[0100] Finally, use the semi-global stereo matching algorithm to calculate the disparity map. The cost function is: , where . Calculate the three-dimensional coordinates through triangulation to generate point cloud data with a density of about 4 points per square meter. Use the Poisson reconstruction algorithm to generate a three-dimensional mesh model, realizing the fine three-dimensional reconstruction of Haidian District.

[0101] Compared with traditional satellite stereo image processing methods, traditional methods mainly adopt serial processing methods and lack an effective task scheduling mechanism. During feature extraction and matching, raw feature descriptors are often directly used, resulting in large computational amounts and being easily affected by noise. When splicing images, a simple sequential splicing strategy is adopted, which is prone to cumulative errors. Processing image data within a range of 100 square kilometers usually takes more than 72 hours. The present invention shortens the processing time to 24 hours through the collaborative computing of GPU and CPU, combined with feature optimization and the generation strategy of the optimal splicing sequence, and ensures the processing accuracy. The planar accuracy of the 3D reconstruction result is better than 1 meter, and the elevation accuracy is better than 1.5 meters, meeting the requirements of urban fine 3D modeling. In terms of processing efficiency, the present invention is 3 times more efficient than traditional methods; in terms of processing accuracy, the planar and elevation accuracies are increased by 30% and 25% respectively. At the same time, the method of the present invention has better robustness and can adapt to the image processing requirements of different scenarios. Through the stability analysis of features and dynamic parameter adjustment, the problem of error accumulation in large-scale image processing is effectively solved.

[0102] It should be noted that the detailed explanations of the variables involved in the present invention are shown in Table 2 below.

[0103] Table 2 Variable Explanation Table

[0104] The above is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention can easily think of changes or substitutions, which should all be covered by the protection scope of the present invention.

Claims

1. A large-scale satellite stereo image data processing method, characterized in that: The following steps are involved: The satellite stereo image data is preprocessed to obtain preprocessed image data; the preprocessed image data is allocated to a GPU processing flow and a CPU processing flow for parallel computing; a multidimensional feature matrix of the preprocessed image data is extracted; a matrix chain multiplication dynamic programming algorithm is used for feature decomposition; a spatial index is established to generate data blocks; feature vectors are extracted and clustered for optimization; a registration model and an evaluation function are established; a splicing graph is constructed and a splicing order is determined; image splicing and depth information reconstruction are performed; wherein the collaborative computing of the GPU and the CPU is realized by computing a task allocation function, a minimum spanning tree algorithm is used to optimize the splicing order, and dynamic optimization is performed based on a matching quality evaluation function.

2. The large-scale satellite stereo image data processing method according to claim 1, characterized in that: The preprocessing step specifically performs Gaussian filtering and wavelet noise reduction on the satellite stereo image data. By setting the Gaussian kernel size to 5 pixels and the standard deviation to 1, the db4 wavelet basis function is selected for 3-layer wavelet decomposition, and the soft threshold method is used for noise reduction to obtain the preprocessed image data.

3. The large-scale satellite stereo image data processing method according to claim 1, characterized in that: The step of allocating computing tasks specifically includes establishing a computing task allocation function, wherein the computing task allocation function includes a matrix operation identification value and a data flow identification value. When the data scale is greater than 1024×1024 pixels and involves matrix operations, the matrix operation identification value is set to true and allocated to the GPU processing flow; when the data volume is greater than 100 megabytes and there are continuous data flow operations, the data flow identification value is set to true and allocated to the CPU processing flow.

4. The large-scale satellite stereo image data processing method according to claim 1, characterized in that: The step of extracting a multidimensional feature matrix specifically comprises performing a 4-layer Haar wavelet decomposition on the preprocessed image data to extract a geographic space feature matrix, a spectral feature matrix and a texture feature matrix; The matrix chain multiplication dynamic programming algorithm is used to determine the optimal calculation sequence, and singular value decomposition is performed to obtain a stable characteristic matrix and a variable characteristic matrix.

5. The large-scale satellite stereo image data processing method according to claim 1, characterized in that: The steps of extracting feature vectors and clustering optimization are as follows: specifically, a scale-invariant feature transformation algorithm is used to extract feature descriptors, an 8-layer scale space is constructed, and each group contains 3 scale layers; a K-means clustering algorithm is used to cluster the feature descriptors, the number of cluster centers is 0.1 times the total number of feature points, the number of iterations is 100, and each feature descriptor is replaced by the nearest cluster center point to form an optimized feature vector group.

6. The large-scale satellite stereo image data processing method according to claim 5, characterized in that: The steps of establishing the registration model are as follows: performing feature matching based on the optimized feature vector group, setting the distance ratio threshold to 0.8, using the random sampling consistency algorithm to eliminate outliers, setting the inlier threshold to 3 pixels, and the number of iterations to 1000 times, and solving the affine transformation parameters by the least squares method to establish the registration model; constructing a matching quality evaluation function, and re-executing the registration model construction when the matching quality evaluation value is less than 0.

75.

7. The large-scale satellite stereo image data processing method according to claim 1, characterized in that: The steps of constructing the mosaic map are to calculate the overlapping area ratio and matching quality assessment value between the image data blocks, and use weighted summation to determine the edge weight value, with the overlapping area ratio weight being 0.3 and the matching quality assessment value weight being 0.7; the Kruskal algorithm is used to solve the minimum spanning tree to determine the mosaic order.

8. The large-scale satellite stereo image data processing method according to claim 1, characterized in that: The steps of performing image stitching and depth information reconstruction are as follows: specifically, the minimum graph cut algorithm is used to calculate the optimal stitching line, and a transition zone with a width of 32 pixels is set; local re-stitching is performed when the grayscale difference value of the local area is greater than 8; the disparity map is calculated using the semi-global matching algorithm, and the disparity search range is set to 0 to 64 pixels, and the three-dimensional surface model is reconstructed through the triangulation principle.

9. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores program instructions, and when the program instructions are executed in a computer, they are used to execute the large-scale satellite stereoscopic image data processing method according to any one of claims 1 to 8.

10. A large-scale satellite stereoscopic image data processing system, characterized in that: The system comprises the computer-readable storage medium as claimed in claim 9, wherein the system is any one of a computer, a server, and a single-chip microcomputer, the computer-readable storage medium is arranged in the system, and the system is provided with a microprocessor for executing program instructions stored in the computer-readable storage medium.

Images