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

Abstract

The present invention provides a large-scale satellite stereo image data processing method, medium and system, which belongs to the field of satellite stereo image data processing technology. First, the image is pre-processed by Gaussian filtering and wavelet noise reduction, and then the collaborative calculation of the GPU and CPU is realized based on the calculation task allocation function. Through multi-scale feature extraction and matrix optimization, a stable feature matrix and a variable feature matrix are obtained. Quadtree spatial index is used for data segmentation, combined with SIFT feature extraction and K-means clustering to optimize feature representation. A stereo image registration model is established, the optimal splicing order is determined by the minimum spanning tree algorithm, and the progressive texture fusion method is used to realize image splicing, and finally a high-precision three-dimensional surface model is generated. The present invention achieves the unity of processing efficiency and accuracy, and solves the technical problem that the existing technology is difficult to improve the processing efficiency of large-scale stereo image data while ensuring processing accuracy.

Description

Technical Field

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

[0002] Satellite stereo imagery data processing technology has widespread applications in geographic information systems, environmental monitoring, urban planning, and other fields. Traditional satellite stereo imagery processing methods primarily rely on a single processing pipeline, employing a sequential data processing approach that includes image preprocessing, feature extraction, image matching, and 3D reconstruction. In practical applications, commonly used processing methods include feature point-based image matching, region-based stereo matching, and 3D reconstruction based on semantic segmentation. These methods perform well when processing small amounts of data and can achieve basic image processing and 3D reconstruction capabilities.

[0003] However, with the rapid development of satellite remote sensing technology, the amount of stereo image data acquired has grown exponentially, posing significant challenges to traditional processing methods. First, a single processing pipeline struggles to fully utilize the hardware resources of modern computers, resulting in low processing efficiency. Second, feature extraction and matching during large-scale data processing are susceptible to noise, impacting processing accuracy. Third, traditional image stitching methods suffer from cumulative errors when processing large numbers of image blocks, making it difficult to ensure geometric accuracy and visual continuity in the resulting stitching. These challenges are particularly pronounced when processing high-resolution stereo image data.

[0004] Existing technologies primarily rely on simple data block processing and serial computation, lacking effective task scheduling mechanisms and feature optimization strategies. In large-scale data processing, this approach fails to achieve optimal allocation of computing resources and cannot guarantee the quality of processing results. Furthermore, existing image stitching techniques often neglect feature stability analysis and dynamic parameter adjustment, leading to problems such as noticeable seams and geometric distortion during the stitching process. Therefore, improving the processing efficiency of large-scale stereo image data while ensuring processing accuracy has become a pressing technical challenge.

[0005] In summary, the prior art has a technical problem in that it is difficult to improve the processing efficiency of large-scale stereoscopic image data while ensuring processing accuracy. Summary of the Invention

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

[0007] The present invention is implemented as follows: In a first aspect, the present invention provides a large-scale satellite stereo image data processing method comprising the following steps: preprocessing satellite stereo image data to obtain preprocessed image data; allocating the preprocessed image data to a GPU processing flow and a CPU processing flow for parallel computing; extracting a multidimensional 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 splicing graph and determining a splicing order; performing image splicing and depth information reconstruction; wherein collaborative computing of the GPU and the CPU is achieved through a computing task allocation function, the splicing 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 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, the preprocessed image data is obtained.

[0009] Among them, the step of allocating computing tasks is specifically to establish a computing task allocation function, which 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 is a continuous data flow operation, the data flow identification value is set to true and allocated to the CPU processing flow.

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

[0011] Among them, the steps of extracting feature vectors and clustering optimization are specifically to extract feature descriptors using a scale-invariant feature transformation algorithm, construct an 8-layer scale space, and each group contains 3 scale layers; cluster the feature descriptors using a K-means clustering algorithm, with the number of cluster centers being 0.1 times the total number of feature points and the number of iterations being 100, and replacing each feature descriptor with the nearest cluster center point to form an optimized feature vector group.

[0012] Among them, the step of establishing the registration model is to perform feature matching based on the optimized feature vector group, set the distance ratio threshold to 0.8, use the random sampling consistency algorithm to eliminate outliers, set the inlier threshold to 3 pixels, and the number of iterations to 1000 times. The registration model is established by solving the affine transformation parameters using the least squares method; construct a matching quality evaluation function, and re-execute the registration model construction when the matching quality evaluation value is less than 0.75.

[0013] Among them, the steps of constructing the mosaic map are to calculate the overlapping area ratio and matching quality evaluation value between 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; the Kruskal algorithm is used to solve the minimum spanning tree to determine the mosaic order.

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

[0015] A second aspect of the present invention provides a computer-readable storage medium, wherein the computer-readable storage medium stores program instructions. When the program instructions are run in a computer, they are used to execute the above-mentioned large-scale satellite stereo image data processing method.

[0016] The third aspect of the present invention provides a large-scale satellite stereoscopic image data processing system, which includes the above-mentioned computer-readable storage medium. 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.

[0017] Compared to existing technologies, the present invention provides a method, medium, and system for processing large-scale satellite stereo imagery data. This method, by establishing a computational task allocation function, enables collaborative computing between GPUs and CPUs. This method utilizes multi-scale feature extraction and matrix optimization strategies, combined with a quadtree spatial index structure, to significantly improve data processing efficiency and accuracy.

[0018] In practical applications, the method of the present invention effectively suppresses image noise through a combination of Gaussian filtering and wavelet noise reduction preprocessing. It improves the reliability of feature matching through feature stability decomposition and dynamic parameter adjustment. It also employs a graph-structure-based optimal splicing sequence generation strategy to reduce the impact of cumulative error. Furthermore, the application of a progressive texture fusion method ensures the visual naturalness of the splicing results. The method of the present invention maintains stable processing performance and accuracy, especially when processing terabytes of stereoscopic image data.

[0019] This invention successfully addresses the efficiency and accuracy challenges in processing large-scale satellite stereo imagery data, primarily due to its innovative computing resource scheduling strategy and feature optimization process. By assigning processing tasks with different characteristics to the most appropriate computing units and employing a multi-level feature extraction and optimization strategy, it achieves a balance between processing efficiency and accuracy, resolving the technical challenge of existing technologies in improving the processing efficiency of large-scale stereo imagery data while maintaining processing accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 A flow chart of the method provided by the present invention. DETAILED DESCRIPTION

[0021] In order to make the purpose, 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] like Figure 1 FIG. 1 is a flow chart of a method for processing large-scale satellite stereo image data provided by the first aspect of the present invention. The method comprises the following steps:

[0023] S01, receiving stereo image data collected by multiple satellites, performing Gaussian filtering and wavelet noise reduction processing on the stereo image data to obtain pre-processed image data;

[0024] S02. Establish a computing task allocation function, wherein the computing task allocation function includes a matrix operation identification value and a data flow identification value, and divide the pre-processed image data into tasks according to the computing task allocation function, allocating image block computing tasks with a true matrix operation identification value to a GPU processing flow for parallel accelerated computing, and allocating continuous data processing tasks with a true data flow identification value to a CPU processing flow for pipeline computing;

[0025] S03, performing multi-scale wavelet decomposition on the pre-processed image data, extracting three-dimensional spatial features, multi-spectral band features and gray-level co-occurrence matrix features, and obtaining a geographic spatial feature matrix, a spectral feature matrix and a texture feature matrix;

[0026] S04, using a matrix chain multiplication dynamic programming algorithm to calculate the minimum multiplication operation order of the geographic spatial feature matrix, the spectral feature matrix, and the texture feature matrix, and performing singular value decomposition according to the minimum multiplication operation order to obtain a stable feature matrix and a variable feature matrix;

[0027] S05, using a quadtree structure to perform spatial division on the stable feature matrix, establish a spatial hierarchical index, and generate an image data block index table;

[0028] S06, dividing the pre-processed image data into blocks according to the image data block index table to generate a number of image data blocks of regular sizes;

[0029] S07, extracting feature descriptors of each image data block using a scale-invariant feature transformation algorithm, clustering the feature descriptors using a K-means clustering algorithm to obtain feature cluster center points;

[0030] S08, replacing each of the feature descriptors with the nearest feature cluster center point to form an optimized feature vector group;

[0031] S09, calculating image registration parameters based on the optimized feature vector group, including rotation angle parameters, scaling parameters, and translation vector parameters, and establishing an affine transformation model as a stereo image registration model;

[0032] S10, calculating the correlation coefficient between the stereo image registration model and the change feature matrix, and establishing a registration parameter correction function;

[0033] S11. Establish a matching quality evaluation function based on the root mean square error criterion, calculate the gradient direction consistency value and the grayscale correlation coefficient between the image blocks, and obtain a matching quality evaluation value. When the matching quality evaluation value is less than a preset matching threshold, return to step S09 to step S11;

[0034] S12, constructing an image data block mosaic graph, wherein a graph node represents the image data block, and a graph edge weight value is a weighted sum of an overlapping area between adjacent image data blocks and the matching quality evaluation value;

[0035] S13, applying a minimum spanning tree algorithm to the image data block splicing graph to determine an optimal splicing order of the image data blocks and generate a splicing operation sequence;

[0036] S14, stitching the image data blocks one by one according to the stitching operation sequence, processing overlapping areas using a progressive texture fusion method, and obtaining a stitched image;

[0037] S15, calculating the grayscale difference value at the boundary of adjacent image blocks in the stitched image, and when the grayscale difference value is greater than a preset stitching threshold, performing local re-stitching optimization on the boundary area;

[0038] S16. Using a stereo matching algorithm to extract disparity information of the stitched image, calculating three-dimensional spatial coordinates, constructing a depth information map, and using a triangular mesh reconstruction algorithm to generate a three-dimensional surface model.

[0039] The specific implementation methods of the above steps are described in detail below. The specific implementation method of step S01 is: for the received satellite stereo image data, first use Gaussian filtering to smooth the image, and by setting the Gaussian kernel size to 5 pixels and the standard deviation to 1, preliminarily suppress the noise in the image. Gaussian filtering can effectively remove Gaussian noise while maintaining the edge features of the image. Then, a wavelet denoising algorithm is used, and the db4 wavelet basis function is selected to perform a three-layer wavelet decomposition on the image. The wavelet coefficients are processed using 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 types of noise generated during the satellite imaging process and improve the accuracy of subsequent processing. The combination of Gaussian filtering and wavelet denoising can better maintain the detailed information of the image while effectively suppressing noise interference. The result of this step is recorded as preprocessed image data.

[0040] The specific implementation of step S02 is as follows: First, a computational task allocation function is constructed, which includes two main judgment criteria: a matrix operation identification value and a data flow identification value. For the input preprocessed image data, its data scale and computational complexity are calculated. When the data scale is greater than 1024×1024 pixels and involves matrix operations, the matrix operation identification value is set to true; when the data processing process involves continuous data flow operations and the data volume is greater than 100 megabytes, the data flow identification value is set to true. Based on these two identification values, the computational tasks are allocated to different processing units: tasks with a true matrix operation identification value are assigned to the GPU, utilizing its parallel computing power to accelerate matrix operations; tasks with a true data flow identification value are assigned to the CPU, utilizing its efficient data pipeline processing capabilities. The purpose of this step is to achieve optimal allocation of computing resources and improve overall computational efficiency.

[0041] The specific implementation method of step S03 is: use discrete wavelet transform to perform multi-scale decomposition on the pre-processed image data, select Haar wavelet as the basis function, perform 4-layer decomposition, and extract low-frequency approximate components and high-frequency detail components respectively. In terms of spatial feature extraction, the gradient operator is used to calculate the horizontal and vertical gradients of the image, and the geographic space feature matrix is ​​constructed by combining the gradient amplitude and direction information. For spectral feature extraction, the principal component analysis method is used to reduce the dimension of the multi-spectral band data, retain the principal components with a contribution rate greater than 95%, and form a spectral feature matrix. In terms of texture feature extraction, the grayscale co-occurrence matrix is ​​calculated, the distance parameter is set to 1, the direction angle is set to 0 degrees, 45 degrees, 90 degrees and 135 degrees, and statistical features such as energy, entropy, contrast, and correlation are extracted to construct a texture feature matrix. This step aims to comprehensively extract the multi-dimensional feature information of the image and provide a reliable feature basis for subsequent processing.

[0042] The specific implementation method 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 cost of different multiplication orders, the order with the smallest total computational amount is selected. In the specific implementation, a cost matrix and a split point matrix are constructed. The matrix dimension is the number of feature matrices. The matrix is ​​filled by the dynamic programming method to finally obtain the optimal computational order. Then, according to the optimal order, a 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 constitute a stable feature matrix, and the rest constitute a variable feature matrix. The purpose of this step is to optimize the matrix calculation efficiency while achieving stable feature decomposition.

[0043] The specific implementation method of step S05 is: use a quadtree structure to construct a spatial index for the stable feature matrix. First, divide the entire feature space into four equal sub-regions, calculate the feature density for each sub-region, and when the feature density is greater than the preset threshold value of 0.8, continue to divide the sub-region 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 feature information it contains 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.

[0044] The specific implementation method of step S06 is: according to the spatial division information in the image data block index table, the pre-processed image data is divided into blocks. First, the basic block size is determined to be 256×256 pixels. Considering the overlapping requirements of subsequent processing, an overlapping area is set when dividing the blocks, and the overlapping area size is 64 pixels. For each data block, its original image data, position information and adjacency relationship are saved. In the block process, for areas where the image edge is less than the basic block size, a mirror filling method is used to expand to the standard size. The purpose of this step is to divide large-scale image data into small blocks suitable for parallel processing, and at the same time, the setting of the overlapping area provides a reliable transition area for subsequent feature matching and image stitching.

[0045] The specific implementation method of step S07 is: using the scale-invariant feature transformation algorithm to extract the feature points of each image data block, including the feature point position, scale and direction information. In the feature point detection stage, a scale space is constructed, the number of scale layers is set to 8, each group contains 3 scale layers, the Gaussian difference operator response is calculated, and stable feature points are screened out through extreme point detection and edge response test. For each feature point, its main direction is calculated to generate a 128-dimensional feature descriptor. Then, the K-means clustering algorithm is used to cluster all feature descriptors, the number of cluster centers is set to 0.1 times the total number of feature points, the Euclidean distance is used as the similarity metric, and the number of iterations is set to 100, until the change in the cluster center position 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.

[0046] The specific implementation method of step S08 is: for each feature descriptor of the image data block, calculate its Euclidean distance with all cluster center points, and select the cluster center point with the smallest distance as the alternative representation of the feature descriptor. In this way, the original 128-dimensional feature descriptor is mapped to the cluster center space to achieve simplification and unification of feature representation. During the replacement process, the distance value between each feature descriptor and the corresponding cluster center is recorded at the same time, and these distance values ​​will be used for subsequent feature matching reliability evaluation. The final optimized feature vector group has a more compact representation and better distinguishability. The purpose of this step is to reduce the redundancy of feature representation and improve the efficiency and stability of the subsequent matching process.

[0047] The specific implementation method of step S09 is: construct a stereo image registration model based on the optimized feature vector group, first use the nearest neighbor distance ratio method to perform feature matching, set the distance ratio threshold to 0.8, and screen out reliable matching point pairs. Use the random sampling consistency algorithm to eliminate external points, set the internal point threshold to 3 pixels, and iterate 1000 times to obtain an accurate set of matching points. Based on the matching point set, calculate the affine transformation parameters, including rotation angle, scaling ratio and translation vector. Solve the transformation parameters by the least squares method to establish a geometric transformation relationship between the image pairs. The purpose of this step is to establish an accurate spatial correspondence between image blocks and provide a geometric transformation basis for subsequent image stitching.

[0048] The specific implementation method of step S10 is: calculate the relationship between the stereo image registration model and the variable feature matrix, first use the registration model to transform each eigenvector in the variable feature matrix, and calculate the difference in feature response before and after the transformation. The normalized cross-correlation method is used to calculate the correlation coefficient, and on this basis, a registration parameter correction function is established. The 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 to improve the adaptability and robustness of the registration model.

[0049] The specific implementation method of step S11 is: constructing a matching quality evaluation function based on the root mean square error criterion, which comprehensively considers the two aspects of gradient direction consistency and grayscale correlation between image blocks. Gradient direction consistency is evaluated by calculating the difference in gradient direction at corresponding points, and the direction difference threshold is set to 30 degrees. Grayscale correlation is calculated using the normalized cross-correlation coefficient, and a sub-block of 32×32 pixels is 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, and the weight ratio is 0.6 for gradient direction consistency and 0.4 for grayscale correlation. When the matching quality evaluation value is less than the preset matching threshold of 0.75, it is necessary to return to re-execute the construction of the registration model and parameter optimization. The purpose of this step is to ensure the accuracy and reliability of image matching.

[0050] The specific implementation method of step S12 is: construct a splicing graph for representing the connection relationship between image data blocks, where 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, the area ratio of their overlapping areas and the matching quality evaluation value are calculated, and the weight value of the edge is determined by weighted summation, 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 area is less than 0.2 or the matching quality evaluation value is less than 0.6, no connecting edge is established between the corresponding nodes. The purpose of this step is to establish the topological relationship between the image data blocks and provide graph structure support for the generation of the optimal splicing sequence.

[0051] The specific implementation method of step S13 is as follows: the image data block splicing graph is converted into a minimum spanning tree problem, and the minimum spanning tree is solved using the Kruskal algorithm. First, the weight values ​​of all edges in the graph are sorted, starting with the edge with the smallest weight, and whether adding the edge will form a loop is determined in turn. If a loop will not be formed, the edge is added to the minimum spanning tree. In the process of adding edges, a union-find data structure is used to maintain the connectivity between nodes to ensure that a tree structure is ultimately generated. Based on the direction of the edges in the minimum spanning tree, the splicing order of the image data blocks is determined, and a splicing operation sequence is generated. The purpose of this step is to optimize the splicing order and minimize the cumulative error.

[0052] The specific implementation of step S14 is as follows: the image data blocks are stitched together according to the stitching operation sequence, and the overlapping areas are processed using a progressive texture fusion method. Within the overlapping area, the minimum graph cut algorithm is first used to calculate the optimal stitching line. The stitching cost function considers both pixel value differences and gradient differences. A transition zone is then set on both sides of the stitching line, with a width of 32 pixels. Within the transition zone, a weighted fusion method is used for smooth transitions, with the weights varying linearly with the distance from the pixel to the stitching line. For each stitching operation, the image's global coordinate system information and stitching status flag are updated. The purpose of this step is to achieve seamless image stitching and eliminate stitching artifacts.

[0053] Step S15 is specifically implemented by calculating the grayscale difference at the boundaries of adjacent image blocks in the stitched image. Using a sliding window method with a window size of 16×16 pixels, the root mean square difference of the local area is calculated point by point along the boundary. When the grayscale difference value of a local area exceeds the preset stitching threshold of 8, the area is marked as requiring optimization. For each area requiring optimization, the optimal stitching line is recalculated, and the transition zone is expanded to 64 pixels. A nonlinear weighting function is used to fuse the grayscale values ​​to achieve optimized stitching of the local area. This step aims to further improve the stitching effect and eliminate local discontinuities.

[0054] The specific implementation of step S16 is as follows: A stereo matching algorithm is used to extract disparity information from the stitched images. First, epipolar correction is performed on the images. A disparity map is calculated using a semi-global matching algorithm. The disparity search range is set to 0 to 64 pixels, the number of cost aggregation paths is set to 8, and the smoothness parameter is set to 0.8. Based on the disparity information and camera parameters, 3D spatial coordinates are calculated using triangulation principles to obtain a depth information map. A 3D surface model is then generated using a triangulated mesh reconstruction algorithm.

[0055] The functions or calculation processes involved in the present invention are described in detail as follows:

[0056] The Gaussian filter kernel function is expressed as follows:

[0057] ;

[0058] Where, is the Gaussian kernel standard deviation, which is 1; is the pixel coordinate position; is the base of natural logarithms.

[0059] The wavelet denoising threshold function is expressed as follows:

[0060] ;

[0061] Where, is the wavelet coefficient; is the soft threshold, ; is the total number of image pixels; is the noise standard deviation; is the indicator function.

[0062] The computing task allocation function is expressed as follows:

[0063] ;

[0064] Where, is the matrix operation complexity index, ranging from 0 to 1; is the data flow index, ranging from 0 to 1; To calculate the dependency index, it ranges from 0 to 1; is the weight coefficient, which takes values ​​of 0.5, 0.3, and 0.2 respectively. When , the task is assigned to the GPU processing flow; otherwise, it is assigned to the CPU processing flow.

[0065] The geospatial feature matrix is ​​represented as follows:

[0066] ;

[0067] Where, For the OK The gradient magnitude of the column pixels.

[0068] The spectral feature matrix is ​​expressed as follows:

[0069] ;

[0070] Where, is the original multispectral data; is the transformation matrix; is the residual matrix.

[0071] The singular value decomposition is expressed as follows:

[0072] ;

[0073] Where, is the input feature matrix; is the left and right singular vector matrix; is a diagonal matrix of singular values; is the residual matrix.

[0074] The K-means clustering objective function is expressed as follows:

[0075] ;

[0076] Where, is the number of clusters; For the clusters; For the cluster centers; is the feature vector.

[0077] The affine transformation matrix is ​​expressed as follows:

[0078] ;

[0079] Where, is the transformation coefficient; is the translation component.

[0080] The registration parameter correction function is expressed as follows:

[0081] ;

[0082] Where, is the parameter correction amount; is the current matching quality; Matching quality to the target; is the gain matrix; is the historical correction amount; is the historical weight coefficient, ranging from 0 to 1.

[0083] The matching quality evaluation function is expressed as follows:

[0084] ;

[0085] Where, is the image block to be matched; is the gradient operator; is the mean; is the standard deviation; is the weight coefficient.

[0086] The edge weight calculation function is expressed as follows:

[0087] ;

[0088] Where, is the overlapping area; is the image block area; For matching quality; is the depth continuity measure; is the weight coefficient.

[0089] The disparity calculation cost function is expressed as follows:

[0090] ;

[0091] Where, is the pixel matching cost; is the smoothness constraint; is the temporal consistency constraint; is the balance coefficient.

[0092] The following is a detailed description of the derivation and establishment process of multiple functions.

[0093] 1. The derivation and establishment process of Gaussian filter kernel function:

[0094] The Gaussian filter kernel function is derived based on the normal distribution in probability theory. The derivation steps are as follows:

[0095] First, construct the basic form of the two-dimensional normal distribution:

[0096] ;

[0097] To ensure that the sum of the filter kernel weights is 1, determine the coefficient A:

[0098] ;

[0099] Finally, the standardized Gaussian filter kernel function is obtained:

[0100] .

[0101] parameter Through experimental verification, when It can achieve a good balance between denoising and preserving edge details.

[0102] 2. The derivation and establishment process of wavelet denoising threshold function:

[0103] Based on Bayesian estimation theory, the derivation steps are:

[0104] Build a noise model:

[0105] ;

[0106] in is the noisy wavelet coefficient, is the real signal, is Gaussian noise;

[0107] Minimize mean squared error:

[0108] ;

[0109] The soft threshold function is obtained by solving:

[0110] .

[0111] Threshold It is estimated using maximum likelihood.

[0112] 3. The derivation and establishment process of the calculation task allocation function:

[0113] First, construct the matrix operation complexity index:

[0114] ;

[0115] Where, For the The number of rows in the matrix; For the The number of columns in the matrix; The coefficient is the matrix operation type, the value is 1 for addition and subtraction, 2 for multiplication and division, and 3 for inverse operation; is the benchmark operation amount, and its value is floating-point operations; is the number of matrices to be processed.

[0116] Constructing data flow indicators:

[0117] ;

[0118] Where, For the The size of the data block; is the data access frequency coefficient, which is 1 for continuous access and 2 for random access; is the system memory capacity; is the number of data blocks.

[0119] Constructing data dependency indicators:

[0120] ;

[0121] Where, is the dependency matrix element, if the data block Dependent data blocks Then take 1, otherwise take 0; is the total number of data blocks.

[0122] The final task allocation function is expressed as:

[0123] ;

[0124] Where, is the weight coefficient, which is determined by the following steps:

[0125] (1) Constructing a judgment matrix:

[0126] ;

[0127] (2) Calculate eigenvalues ​​and eigenvectors;

[0128] (3) Normalization process to obtain the weight coefficient: .

[0129] when When , the task is assigned to the GPU processing flow; otherwise, it is assigned to the CPU processing flow.

[0130] The derivation and establishment process of the characteristic matrix:

[0131] The spatial feature matrix is ​​based on multi-scale gradient analysis:

[0132] Calculate the horizontal and vertical gradients:

[0133] ;

[0134] Calculate the gradient magnitude:

[0135] ;

[0136] Form the feature matrix S.

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

[0138] The derivation and establishment process of the registration parameter correction function:

[0139] Design based on feedback control theory:

[0140] (1) Establish error model:

[0141] ;

[0142] (2) Introducing proportional integral control:

[0143] ;

[0144] (3) Consider the impact of historical revisions:

[0145] .

[0146] Parameter K is obtained through system identification. Determined by cross validation.

[0147] 6. Derivation and establishment process of matching quality evaluation function:

[0148] Comprehensively consider gradient consistency and grayscale correlation:

[0149] (1) Gradient consistency measure:

[0150] ;

[0151] (2) Grayscale correlation measurement:

[0152] ;

[0153] (3) The evaluation function is obtained by weighted combination.

[0154] The weight coefficients are obtained through regression analysis of a large amount of experimental data.

[0155] A second aspect of the present invention provides a computer-readable storage medium, wherein the computer-readable storage medium stores program instructions. When the program instructions are run in a computer, they are used to execute the above-mentioned large-scale satellite stereo image data processing method.

[0156] The third aspect of the present invention provides a large-scale satellite stereoscopic image data processing system, which includes the above-mentioned computer-readable storage medium. 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.

[0157] Specifically, the present invention is based on the analysis of computing task characteristics and multi-level feature optimization. First, by analyzing the computing characteristics of different processing tasks, a task allocation mechanism based on matrix operation identifiers and data flow identifiers is established. Matrix operation-intensive tasks are better suited for parallel processing on GPUs, while data flow-intensive tasks are more suitable for pipeline processing on CPUs. This allocation mechanism fully leverages the advantages of heterogeneous computing platforms.

[0158] In terms of feature processing, this invention uses multiscale wavelet decomposition to extract three-dimensional spatial features, multispectral band features, and texture features, which describe image characteristics from different perspectives. A matrix chain multiplication dynamic programming algorithm optimizes the feature calculation sequence, reducing computational complexity. Stable feature decomposition ensures the retention of key features, while dynamic parameter adjustment through variable feature matrices enhances system adaptability.

[0159] For image stitching, this invention constructs a quadtree-based spatial index structure, enabling efficient spatial retrieval. The optimal stitching sequence is determined using a minimum spanning tree algorithm, combined with a progressive texture fusion method to ensure the geometric accuracy and visual quality of the stitching results. This multi-level optimization strategy ensures the efficiency and reliability of the entire processing process.

[0160] 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.

[0161] The specific implementation of step S01 is: for the received satellite stereo image data, firstly, a Gaussian filter is used to smooth the image. The Gaussian filter kernel function is expressed as follows:

[0162] ;

[0163] Where, is the Gaussian kernel standard deviation, which is 1; is the pixel coordinate position; is the base of the natural logarithm. Gaussian filtering is implemented using a separable kernel. First, horizontal filtering is performed, then vertical filtering is performed. The filtered image is represented as:

[0164] ;

[0165] Where, is the filtered image; is the original image; Represents the convolution operation. Then, the wavelet denoising algorithm is used, and the db4 wavelet basis function is selected to perform a three-layer wavelet decomposition on the image. The wavelet coefficients are processed using the soft threshold method. The threshold function is expressed as follows:

[0166] ;

[0167] Where, is the wavelet coefficient; is the soft threshold, ; is the total number of image pixels; is the noise standard deviation; The purpose of this step is to eliminate all kinds of noise generated during satellite imaging and improve the accuracy of subsequent processing.

[0168] The specific implementation of step S02 is: constructing a computing task allocation function to achieve optimal allocation of computing resources by evaluating the matrix operation complexity index, data flow index and data dependency index. The matrix operation complexity index is expressed as follows:

[0169] ;

[0170] Where, For the The number of rows in the matrix; For the The number of columns in the matrix; The coefficient is the matrix operation type, the value is 1 for addition and subtraction, 2 for multiplication and division, and 3 for inverse operation; is the benchmark operation amount, and its value is floating-point operations; is the number of matrices to be processed. The data flow index is expressed as follows:

[0171] ;

[0172] Where, For the The size of the data block; is the data access frequency coefficient, which is 1 for continuous access and 2 for random access; is the system memory capacity; is the number of data blocks. The data dependency index is expressed as follows:

[0173] ;

[0174] Where, is the dependency matrix element, if the data block Dependent data blocks Then take 1, otherwise take 0; is the total number of data blocks. The final task allocation function is expressed as:

[0175] ;

[0176] Where, is the weight coefficient, which is determined by the hierarchical analysis method. When , the task is assigned to the GPU processing flow, otherwise it is assigned to the CPU processing flow.

[0177] The specific implementation of step S03 is: multi-dimensional feature extraction is performed on the pre-processed image data. First, discrete wavelet transform is used to perform multi-scale decomposition. Haar wavelet is selected as the basis function and a 4-layer decomposition is performed. The coefficient matrix after each layer of decomposition is expressed as follows:

[0178] ;

[0179] Where, is the number of decomposition layers; is the coefficient number. In terms of spatial feature extraction, the Sobel operator is used to calculate the gradient information of the image:

[0180] ;

[0181] ;

[0182] Where, are the gradients in the horizontal and vertical directions respectively; is the input image. The gradient magnitude and direction are calculated as follows:

[0183] ;

[0184] ;

[0185] Construct a geospatial feature matrix:

[0186] .

[0187] The specific implementation of step S04 is: using matrix chain multiplication dynamic programming algorithm to determine the optimal calculation order of the feature matrix, and performing singular value decomposition to achieve feature separation. First, construct the feature matrix multiplication operation cost matrix, for the matrix sequence , where the matrix The dimension is , the solution formula for the sub-problem of the optimal calculation order is:

[0188] ;

[0189] Where, Represents calculation The minimum number of scalar multiplications required. After determining the optimal calculation order, perform singular value decomposition on each characteristic matrix:

[0190] ;

[0191] Where, is the input feature matrix; is the left and right singular vector matrix; is a diagonal matrix of singular values; is the residual matrix. The extraction of the stable feature matrix is ​​based on the singular value threshold selection:

[0192] ;

[0193] Where, is the threshold coefficient, which is 0.1; is the maximum singular value.

[0194] The specific implementation of step S05 is: constructing a spatial index for the stable feature matrix based on a quadtree structure, and the feature density calculation formula in the spatial partitioning process is:

[0195] ;

[0196] Where, For blocks The characteristic density of For the block The weight of feature points; is the block area; is the number of feature points in the block. The formula for calculating the feature point weight is:

[0197] ;

[0198] Where, is the distance from the feature point to the center of the block; is a scale parameter, which is set to 0.2 times the block side length. The block splitting condition is:

[0199] ;

[0200] Where, is the density threshold, which is 0.8; is the minimum block area, and its value is Pixels.

[0201] The specific implementation of step S06 is: dividing the pre-processed image data into blocks according to the spatial index, and the block size calculation formula is:

[0202] ;

[0203] Where, is the side length of the block; is the GPU memory size; is the number of image channels. The overlapping area size of adjacent blocks is calculated as follows:

[0204] ;

[0205] Where, is the width of the overlapped area. For areas where the edge of the image is less than a complete block, mirror filling is used:

[0206] ;

[0207] Where, are the width and height of the image respectively.

[0208] The specific implementation of step S07 is: using the scale-invariant feature transformation algorithm to extract feature points, first constructing the scale space:

[0209] ;

[0210] Where, is a Gaussian kernel with variable scale; Is the scale parameter. Gaussian difference image calculation:

[0211] ;

[0212] Where, is the scale factor, and its value is The local extreme value detection conditions of feature points are:

[0213] ;

[0214] Where, is the detection threshold, which 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:

[0215] ;

[0216] .

[0217] The specific implementation of step S08 is: cluster optimization is performed on the extracted feature descriptors, and the distance matrix between the feature descriptors is first calculated:

[0218] ;

[0219] Where, For the The first feature descriptor Quantity is the descriptor dimension, and its value is 128. K-means clustering algorithm is used for feature classification, and the clustering objective function is:

[0220] ;

[0221] Where, is the number of clusters, which is 0.1 times the total number of feature points; For the clusters; For the Cluster centers. The update formula of cluster centers is:

[0222] ;

[0223] Where, For the The number of samples in a cluster. The clustering termination condition is:

[0224] ;

[0225] Where, is the convergence threshold, and its value is 0.01.

[0226] The specific implementation of step S09 is: constructing a stereo image registration model based on the optimized feature vector group, first performing feature matching, and performing preliminary screening using the nearest neighbor distance ratio method:

[0227] ;

[0228] Where, are the nearest neighbor and next nearest neighbor distances respectively; is the distance ratio threshold, which is set to 0.8. The affine transformation matrix of the matching point pair is solved:

[0229] ;

[0230] The transformation parameters are solved by the least squares method:

[0231] ;

[0232] Where, are corresponding point pairs; is the number of matching point pairs.

[0233] The specific implementation of step S10 is: establishing a registration parameter correction function, and the correction amount of the registration parameter is calculated as follows:

[0234] ;

[0235] Where, is the current matching quality evaluation value; is the target matching quality value, which is 0.9; is the gain matrix, obtained through system identification:

[0236] ;

[0237] is the historical correction amount, and its update formula is:

[0238] ;

[0239] Where, is the smoothing coefficient, and its value is 0.7.

[0240] The specific implementation of step S11 is: constructing a matching quality evaluation function, comprehensively considering gradient consistency and grayscale correlation:

[0241] ;

[0242] Where, are the weights of gradient consistency and grayscale correlation, with values ​​of 0.6 and 0.4 respectively. The evaluation area uses a sliding window method:

[0243] ;

[0244] Where, The window radius is 16 pixels.

[0245] The specific implementation of step S12 is: constructing an image data block mosaic graph, and the edge weight calculation function of the graph is:

[0246] ;

[0247] Where, is the overlapping area; is the image block area; For matching quality; is the depth continuity measure:

[0248] ;

[0249] Where, is the number of pixels in the overlapping area; is the depth value of the corresponding point.

[0250] The specific implementation of step S13 is: based on the constructed image data block mosaic graph, the minimum spanning tree is solved using the Kruskal algorithm. First, the edges in the graph are sorted in ascending order according to their weights. The edge sorting criteria are:

[0251] ;

[0252] Where, is an edge in the graph; is the corresponding edge weight. A union-find data structure is used to maintain node connectivity, and the lookup function for the set representative element is:

[0253] ;

[0254] Where, For nodes The parent node of . The optimization function of the set merge operation is:

[0255] ;

[0256] Where, For nodes The height of the tree.

[0257] 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 areas. First, use the minimum graph cut algorithm to calculate the optimal stitching line. The graph cut energy function is defined as:

[0258] ;

[0259] Where, is a data item representing a pixel Marked as the price; Is a smoothing term that represents the penalty for differences in adjacent pixel labels:

[0260] ;

[0261] ;

[0262] Where, is the smoothing coefficient, the value is 10; is the cutoff threshold, and its value is 30.

[0263] The specific implementation of step S15 is: perform local optimization on the spliced ​​image and calculate the grayscale difference value of the boundary area:

[0264] ;

[0265] Where, is the local window size, which is set to 16. When the grayscale difference value is greater than the threshold, adaptive weight fusion is used:

[0266] ;

[0267] The weight function is defined as:

[0268] ;

[0269] Where, is the distance from the pixel to the seam line; is the smoothing parameter, and its value is 32.

[0270] The specific implementation of step S16 is: using a semi-global matching algorithm to extract disparity information and constructing a cost function for disparity calculation:

[0271] ;

[0272] Where, is the pixel matching cost:

[0273] ;

[0274] is the smoothness constraint:

[0275] ;

[0276] is the timing consistency constraint:

[0277] ;

[0278] Where, is the smoothness penalty coefficient, with values ​​of 0.8 and 1.6 respectively; are timing constraint parameters, with values ​​of 0.5 and 2 respectively. Finally, the three-dimensional coordinates are calculated by triangulation:

[0279] ;

[0280] ;

[0281] ;

[0282] Where, is the camera focal length; is the baseline length; is the disparity value; The main point coordinates.

[0283] The parameters in all the above steps were determined through extensive experimental verification, ensuring a balance between computational efficiency and accuracy while maintaining algorithm robustness. The entire processing flow fully considers the characteristics of satellite imagery and employs a multi-level optimization strategy to achieve efficient and reliable stereo image processing.

[0284] In order to better understand and implement the present invention, Example 2 of a specific application scenario of the present invention is provided below:

[0285] A research institute, working on a 3D urban modeling project, needed to process high-resolution satellite stereo imagery covering approximately 100 square kilometers. This data set consisted of five primary and five auxiliary images, with a spatial resolution of 0.5 meters. Each image measured 25,000 × 20,000 pixels, totaling approximately 60GB. The research team used the method presented in this paper to process this data.

[0286] First, the research team preprocessed the received raw image data. A 5×5 pixel Gaussian filter kernel was used for smoothing. The Gaussian kernel function was: ,in Through Gaussian filtering, random noise in the image is effectively suppressed, and the image signal-to-noise ratio is improved from the original 26dB to 31dB. Then, the db4 wavelet basis function is used to perform a three-layer wavelet decomposition on the image, and a soft threshold function is used: Perform coefficient processing, where ,After wavelet reconstruction, the denoised preprocessed image is obtained.

[0287] The research team configured two NVIDIA A100 GPUs and one dual-core 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 flow, and the data reading and image stitching tasks are assigned to the CPU processing flow. The specific task allocation is shown in Table 1.

[0288] Table 1: Computational task allocation table

[0289]

[0290] The research team performed multi-scale feature extraction on the preprocessed imagery. They constructed a geospatial feature matrix S with dimensions of 25,000 × 20,000 by calculating gradient magnitudes. Principal component analysis was performed on the multispectral data to obtain a spectral feature matrix P, retaining the top three principal components with a contribution of 97.3%. The energy, entropy, contrast, and correlation characteristics of the gray-level co-occurrence matrix were calculated to construct a texture feature matrix.

[0291] The matrix chain multiplication dynamic programming algorithm is used to determine the optimal calculation order and perform singular value decomposition on the characteristic matrix: The singular value threshold is set 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 and has a dimension of 25000×1000.

[0292] The research team used a quadtree structure to spatially partition the stable feature matrix. The minimum region size was set to 32×32 pixels, and the feature density threshold was set to 0.8. The resulting spatial index tree was eight levels deep and contained approximately 12,000 leaf nodes. Based on this index structure, the image data was divided into 2048 basic data blocks, each 256×256 pixels in size, with a 64-pixel overlap between adjacent blocks.

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

[0294] The affine transformation model is constructed based on the optimized eigenvector group, and the transformation matrix is: . Use the RANSAC algorithm to select matching point pairs, and set the inlier threshold to 3 pixels. Calculate the correction amount of the registration parameters: ,in , iterative optimization until the matching quality evaluation value is greater than 0.75.

[0295] The research team constructed a mosaic of image data blocks and calculated edge weights: ,in , , The Kruskal algorithm was used to solve the minimum spanning tree and determine the stitching order. A 32-pixel-wide progressive blending transition zone was used in the overlapping area to smooth the stitching boundary.

[0296] Finally, the semi-global stereo matching algorithm is used to calculate the disparity map, and the cost function is: ,in The 3D coordinates were calculated through triangulation to generate point cloud data with a density of approximately 4 points per square meter. A Poisson reconstruction algorithm was used to generate a 3D mesh model, achieving a detailed 3D reconstruction of the Haidian District.

[0297] Compared to traditional satellite stereo image processing methods, which primarily employ serial processing and lack an effective task scheduling mechanism, the original feature descriptors are often used directly during feature extraction and matching, resulting in high computational complexity and susceptibility to noise. A simple sequential stitching strategy is employed for image stitching, which can easily lead to cumulative errors. Processing image data covering an area of ​​100 square kilometers typically requires over 72 hours. This present invention reduces processing time to 24 hours while ensuring accuracy through the collaborative computing of the GPU and CPU, combined with feature optimization and an optimal stitching sequence generation strategy. The 3D reconstruction results achieve a planar accuracy better than 1 meter and an elevation accuracy better than 1.5 meters, meeting the requirements for detailed 3D modeling of cities. This method achieves a three-fold improvement in processing efficiency compared to traditional methods, while improving both planar and elevation accuracy by 30% and 25%, respectively. Furthermore, the present method exhibits enhanced robustness and adapts to the image processing requirements of diverse scenarios. Through feature stability analysis and dynamic parameter adjustment, it effectively addresses the error accumulation problem in large-scale image processing.

[0298] It should be noted that the variables involved in the present invention are explained in detail as shown in Table 2 below.

[0299] Table 2 Variable Explanation Table

[0300]

[0301] The above description is only a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with this technical field can easily think of changes or replacements within the technical scope disclosed by the present invention, which should be covered by the scope of protection 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 assigned to the GPU processing flow and the CPU processing flow for parallel computing; the multidimensional feature matrix of the preprocessed image data is extracted; the 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 mosaic graph is constructed and a mosaic order is determined; image mosaicking and depth information reconstruction are performed; wherein the collaborative computing of the GPU and the CPU is realized by the computing task allocation function, the mosaic order is optimized by the minimum spanning tree algorithm, and dynamic optimization is performed based on the matching quality evaluation function; wherein the step of establishing the registration model is specifically to perform feature matching based on the optimized feature vector, set the distance ratio threshold to 0.8, use the random sampling consistency algorithm to eliminate outliers, set the inlier threshold to 3 pixels, and the number of iterations to 1000 Second, the affine transformation parameters are solved by the least squares method to establish a registration model; a matching quality evaluation function is constructed, and the registration model construction is re-executed when the matching quality evaluation value is less than 0.75; the step of constructing the mosaic map is specifically to calculate the overlapping area ratio and matching quality evaluation value between the image data blocks, and the weighted sum is used to determine the edge weight value, the overlapping area ratio weight is 0.3, and the matching quality evaluation value weight is 0.7; the Kruskal algorithm is used to solve the minimum spanning tree to determine the mosaic order; the steps of performing image mosaicking and depth information reconstruction are specifically to calculate the optimal stitching line using the minimum graph cut algorithm and set a transition zone with a width of 32 pixels; local re-splicing 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, the disparity search range is set to 0 to 64 pixels, and the three-dimensional surface model is reconstructed by the triangulation principle; the computing task allocation function is expressed as follows: ； Where, is the matrix operation complexity index, is the data flow indicator, To calculate the dependency index, is the weight coefficient; when When , the task is assigned to the GPU processing flow; otherwise, it is assigned to the CPU processing flow; ； Where, For the The number of rows in the matrix; For the The number of columns in the matrix; The coefficient is the matrix operation type, the value is 1 for addition and subtraction, 2 for multiplication and division, and 3 for inverse operation; is the benchmark operation amount, and its value is floating-point operations; is the number of matrices to be processed; ； Where, For the The size of the data block; is the data access frequency coefficient, which is 1 for continuous access and 2 for random access; is the system memory capacity; is the number of data blocks; ； Where, is the dependency matrix element, if the data block Dependent data blocks Then take 1, otherwise take 0; is the total number of data blocks; The matching quality evaluation function is expressed as follows: ； Where, is the image block to be matched; is the gradient operator; is the mean; is the standard deviation; is the weight coefficient; is the pixel coordinate position; The steps of extracting feature vectors and clustering optimization are specifically as follows: extracting feature descriptors for each image data block using a scale-invariant feature transformation algorithm, clustering the feature descriptors using a K-means clustering algorithm, and obtaining feature cluster centers; wherein, extracting feature descriptors includes extracting feature points using a scale-invariant feature transformation algorithm, calculating the main direction of each feature point, and generating a 128-dimensional feature descriptor; the following is the feature point calculation process: first, constructing a scale space: ； Where, is a Gaussian kernel with variable scale; is the scale parameter, is the original image, is the pixel coordinate position, Gaussian difference image calculation: ； Where, is the scale factor, and the local extreme value detection condition of the feature point is: ； Where, is the detection threshold, For the 26 neighborhoods of the feature points, the main direction of the feature points is calculated based on the gradient histogram: ； ； Among them, the edge weight calculation function of the splicing graph is: ； Where, is the overlapping area; is the image block area; To match the quality, is the weight coefficient; is the depth continuity measure: ； Where, is the number of pixels in the overlapping area; is the depth value of the corresponding point.

2. The large-scale satellite stereo image data processing method according to claim 1, characterized in that: The preprocessing step specifically involves performing 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, selecting the db4 wavelet basis function for 3-layer wavelet decomposition, and using the soft threshold method for noise reduction, the preprocessed image data is obtained.

3. 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 pre-processed 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.

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

5. A large-scale satellite stereo image data processing system, characterized in that: The system comprises the computer-readable storage medium according to claim 4, 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