A method, medium, and system for real-scene 3D modeling based on stereo image pairs and satellite data.

Abstract

This invention provides a method, medium, and system for real-scene 3D modeling based on stereo image pair satellite data, belonging to the field of 3D modeling technology. The invention divides stereo image pair satellite images into geographic grids and constructs Gaussian pyramids to extract multi-scale features. It then inputs these features into a multi-scale matching fusion model based on the principles of two dam water storage problems to obtain corresponding point sets with linear time complexity. A two-layer game theory model is used to optimize resource scheduling of the image processing unit. Based on the corresponding point sets, a dense parallax field is calculated and a 3D point cloud is reconstructed. Finally, a 3D mesh model is generated through registration fusion and a moving cube algorithm. Texture mapping and illumination consistency processing are applied to the 3D mesh model to output a real-scene 3D model of the target area. This invention solves the technical problem of low feature matching efficiency for large-scale stereo image pair satellite images.

Description

Technical Field

[0001] This invention belongs to the field of 3D modeling technology, and more specifically, relates to a method, medium, and system for real-scene 3D modeling based on stereo image pairs and satellite data. Background Technology

[0002] Real-scene 3D modeling technology extracts and matches features from satellite imagery using stereo images, then reconstructs 3D point clouds through parallax calculation. This technology has been widely applied in urban planning, topographic mapping, disaster monitoring, and resource surveys. Traditional methods use scale-invariant feature transformation algorithms to extract key points, perform global search matching by calculating the Euclidean distance between feature descriptors, then use a random sampling consensus algorithm to remove mismatched points, and finally reconstruct the 3D scene based on triangulation principles. However, in practical applications, traditional methods face severe efficiency bottlenecks when processing large-scale satellite imagery. This is mainly due to the brute-force search strategy for feature matching, requiring a global comparison of every feature point in both left and right images. The computational complexity increases quadratically with the image size and the number of feature points. When processing satellite imagery covering hundreds of square kilometers, the feature matching stage often takes hours or even days, severely restricting the real-time performance and practicality of 3D modeling. Because existing technologies lack efficient feature matching search strategies, they cannot effectively reduce the computational complexity of the matching process. Furthermore, the lack of a reasonable resource scheduling mechanism in large-scale image processing tasks leads to low utilization of computational resources. In other words, existing technologies suffer from the technical problem of low efficiency in feature matching of large-scale stereo images to satellite imagery. Summary of the Invention

[0003] In view of this, the present invention provides a method, medium and system for real-scene 3D modeling based on stereo image pair satellite data, which can solve the technical problem of low efficiency in feature matching of large-scale stereo image pair satellite images in the prior art.

[0004] The present invention is implemented as follows: The first aspect of the present invention provides a method for real-scene 3D modeling based on stereo image pairs satellite data, including acquiring stereo image pairs satellite image data of a target area, performing radiometric and geometric corrections, establishing a geographic grid block index, dividing the area into several image processing units according to the geographic grid block index; constructing a Gaussian pyramid for each image processing unit to perform multi-scale decomposition and extract scale-invariant key points, calculating the directional gradient histogram descriptor to establish a multi-scale feature set; inputting the multi-scale feature set into a multi-scale matching and fusion model for feature matching to obtain the corresponding point set, wherein the multi-scale matching and fusion model utilizes two dams... The principle of water storage problem is to perform linear traversal of the histogram of directional gradients descriptor sequence through a two-pointer strategy. The pointer is guided to move according to the Gaussian difference response value, and the histogram of directional gradients descriptor pairs with feature similarity higher than the matching threshold are accumulated as effective matching capacity. A two-layer game model is used to schedule resources for image processing units. Based on the corresponding point set and camera parameters, a disparity map is calculated and scale-space consistency constraint processing is performed. A dense disparity field is obtained through multi-scale disparity propagation for 3D point cloud reconstruction. The 3D point cloud is registered and fused to generate a 3D mesh model. Texture mapping and illumination consistency processing are performed to output a realistic 3D model.

[0005] The geographic grid block index is divided according to latitude and longitude intervals, with the size of a single grid ranging from 0.5 degrees to 2 degrees, establishing a spatial mapping relationship between stereo image pairs, satellite image data, and geographic coordinates.

[0006] The Gaussian pyramid is constructed by performing multiple Gaussian filters and downsampling on the image processing unit, with 4 to 6 layers and a scale factor between adjacent layers ∈ [1.414, 1.732].

[0007] The scale-invariant key points are obtained by calculating local extrema in the scale space. Gaussian difference response values ​​are calculated for each candidate point, and candidate points with Gaussian difference response values ​​greater than an adaptive threshold are retained as scale-invariant key points.

[0008] The directional gradient histogram descriptor is constructed by the gradient direction distribution within the neighborhood of the statistically scale-invariant keypoint. The neighborhood of the scale-invariant keypoint is divided into 4×4 sub-regions, and gradient histograms in 8 directions are calculated in each sub-region, forming a 128-dimensional directional gradient histogram descriptor vector.

[0009] The multi-scale matching fusion model includes an input layer that receives directional gradient histogram descriptors from different scale levels, a feature encoding layer containing three residual convolutional blocks, a multi-scale fusion layer that uses a feature pyramid network structure to fuse features at different scales, and a matching head layer containing two fully connected layers that output matching confidence and keypoint coordinate offsets.

[0010] The multi-scale matching fusion model utilizes an imbalanced data processing mechanism based on adaptive sampling to dynamically adjust the sampling probability of each matching difficulty category according to the sample difficulty coefficient, and extracts samples according to the adjusted sampling probability in each training batch to construct the training batch.

[0011] The principle of the two dam water storage problems is implemented by treating the histogram of directional gradients (HDRs) descriptor sequence of the left image as the boundary of the first container and the histogram of directional gradients (HDRs) descriptor sequence of the right image as the boundary of the second container. Feature similarity is regarded as the water storage height between the first and second containers. The first and second pointers are initialized to point to the starting positions of the histogram of directional gradients (HDRs) descriptor sequences of the left and right images, respectively. If the Gaussian difference response value of the histogram of directional gradients corresponding to the first pointer is less than the Gaussian difference response value of the histogram of directional gradients corresponding to the second pointer, the first pointer is moved one position to the right; otherwise, the second pointer is moved one position to the right. During the pointer movement, the number of histogram of directional gradients (HDRs) descriptor pairs with feature similarity higher than the matching threshold is accumulated as the effective matching capacity.

[0012] The matching confidence is determined by calculating the ratio of the distance between the histogram descriptors of the directional gradients of corresponding point pairs in the corresponding point set to the distance between the histogram descriptors of the directional gradients of the suboptimal matching point. When the distance ratio is less than 0.8, the matching confidence is considered to be high.

[0013] In the two-layer game model, the upper-layer reconstruction accuracy optimization model aims at the point cloud density and surface smoothness of the 3D point cloud reconstruction, while the lower-layer computing resource allocation model aims at the processing time and memory usage of the image processing unit. The upper-layer reconstruction accuracy optimization model and the lower-layer computing resource allocation model are coupled through the processing priority parameter of the image processing unit.

[0014] The camera parameters include an intrinsic camera matrix and an extrinsic camera matrix. The intrinsic camera matrix contains the focal length and principal point coordinates, while the extrinsic camera matrix contains a rotation matrix and a translation vector. The disparity map is obtained by calculating the pixel coordinate difference between corresponding point pairs in the corresponding point set. The disparity value represents the positional offset of the same spatial point in the left and right images. The scale-space consistency constraint processing is performed by comparing the disparity value differences between adjacent scale levels. When the disparity value difference between adjacent scale levels exceeds the inter-scale difference threshold, the current disparity value is smoothed and filtered for correction.

[0015] The multi-scale disparity propagation is achieved by iteratively updating the dense disparity field. The initial disparity field is generated by interpolation of the corresponding point set. In each iteration, for each pixel in the dense disparity field, neighboring pixels with similar disparity values ​​are searched. The disparity values ​​of the neighboring pixels are weighted and averaged to update the disparity value of the pixel.

[0016] The registration and fusion is achieved by calculating the rigid transformation matrix between the 3D point clouds reconstructed by adjacent image processing units. The rigid transformation matrix includes a rotation matrix and a translation vector. The rigid transformation matrix is ​​solved by minimizing the spatial distance error between corresponding point pairs.

[0017] The texture mapping maps the color information of the stereo image pair satellite image data onto the surface of the three-dimensional mesh model. The corresponding texture color is obtained by calculating the projection coordinates of the vertices of the three-dimensional mesh model on the image plane. The illumination consistency processing adjusts the brightness and color of the stereo image pair satellite image data from different perspectives to ensure that the texture color maintains a smooth transition at the stitching point.

[0018] A second aspect of the present invention provides a computer-readable storage medium storing program instructions that, when executed in a computer, are used to perform the above-described method for real-scene 3D modeling based on stereo image pair satellite data.

[0019] A third aspect of the present invention provides a real-scene 3D modeling system based on stereo image pair satellite data, comprising the aforementioned computer-readable storage medium, wherein the system is a computer, the computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor for executing program instructions stored in the computer-readable storage medium.

[0020] This invention addresses the technical problem of low efficiency in feature matching for large-scale stereo image pairs by introducing a feature matching strategy based on the principle of a two-dam water storage problem and combining it with a two-level game theory model for resource scheduling optimization. Specifically, the principle of the two-dam water storage problem transforms feature matching into a geometric optimization problem of container water storage. A dual-pointer strategy performs linear traversal through an ordered sequence of feature descriptors, using Gaussian difference response values ​​to guide pointer movement. This avoids the redundant computation of comparing all feature point pairs one by one in traditional brute-force search methods, reducing the time complexity of the matching search from quadratic to linear, significantly reducing computation time while maintaining matching quality. Furthermore, the two-level game theory model, through collaborative decision-making between upper-level reconstruction accuracy optimization and lower-level computational resource allocation, determines processing priorities based on the contribution of each image processing unit to the overall reconstruction quality, achieving rational scheduling of computational resources and further improving the overall efficiency of large-scale image processing tasks. In summary, this invention solves the technical problem of low efficiency in feature matching for large-scale stereo image pairs mentioned in the background art. Attached Figure Description

[0021] Figure 1 This is a schematic diagram of the multi-scale feature point distribution in the embodiment.

[0022] Figure 2 This is a comparison chart of feature matching times in the embodiments.

[0023] Figure 3 This is a diagram illustrating the iterative convergence process of the two-layer game model in the embodiment.

[0024] Figure 4 This is a cross-sectional distribution diagram of the reconstructed three-dimensional point cloud in the embodiment. Detailed Implementation

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

[0026] The first aspect of this invention provides a method for real-scene 3D modeling based on stereo image pairs and satellite data, comprising the following steps:

[0027] S1. Acquire stereo image pairs satellite image data of the target area, perform radiometric and geometric corrections on the stereo image pairs satellite image data, establish a geographic grid block index, and divide the stereo image pairs satellite image data into several image processing units according to the geographic grid block index.

[0028] S2. Construct a Gaussian pyramid for each image processing unit to perform multi-scale decomposition, extract scale-invariant key points at each scale level, calculate the orientation gradient histogram descriptor of each scale-invariant key point, and establish a multi-scale feature set for the image processing unit.

[0029] S3. Input the multi-scale feature set into the multi-scale matching fusion model for feature matching to obtain the corresponding point set between stereo image pairs. Judge the matching quality based on the matching confidence of the corresponding point set. If the matching confidence is lower than the matching confidence threshold, adjust the scale level range of feature extraction and repeat step S2.

[0030] S4. A two-layer game model consisting of an upper-layer reconstruction accuracy optimization model and a lower-layer computing resource allocation model is used to schedule the reconstruction tasks of the image processing unit. The processing priority and computing resource allocation scheme of each image processing unit are determined according to the game equilibrium solution.

[0031] S5. Calculate the disparity map based on the corresponding point set and the camera parameters of the image processing unit, perform scale-space consistency constraint processing on the disparity map, obtain a dense disparity field through multi-scale disparity propagation, and reconstruct the 3D point cloud based on the dense disparity field and camera parameters.

[0032] S6. Register and fuse the 3D point clouds reconstructed by each image processing unit, generate a 3D mesh model using the moving cube algorithm, perform texture mapping and lighting consistency processing on the 3D mesh model, and output a real-world 3D model of the target area.

[0033] The geographic grid block index is used to establish a spatial mapping relationship between stereo image pair satellite image data and geographic coordinates. The input includes the geographic range and image resolution of the stereo image pair satellite image data, and the output is a correspondence table between grid number and geographic coordinate. The geographic grid block index is divided according to latitude and longitude intervals, and the size range of a single grid is 0.5 degrees to 2 degrees.

[0034] The Gaussian pyramid is constructed by performing multiple Gaussian filters and downsampling on the image processing unit. The number of layers in the Gaussian pyramid is 4 to 6, and the scale factor between adjacent layers is 1.414 to 1.732. The value of the scale factor is obtained by experimentally and statistically analyzing the feature extraction stability under different scale factors. The range of scale factors that makes the key point repetition detection rate greater than 85% while keeping the scale unchanged is selected as the value range of the scale factor. The Gaussian pyramid is used to obtain multi-scale image representation from coarse to fine.

[0035] The scale-invariant keypoints are obtained by calculating local extrema in the scale space. Gaussian difference response values ​​are calculated for each candidate point, and candidate points with Gaussian difference response values ​​greater than an adaptive threshold are retained as scale-invariant keypoints. The adaptive threshold is determined based on the image contrast at the current scale level.

[0036] The directional gradient histogram descriptor is constructed by the gradient direction distribution within the neighborhood of the statistically scale-invariant keypoint. The neighborhood of the scale-invariant keypoint is divided into 4×4 sub-regions, and gradient histograms in 8 directions are calculated in each sub-region, forming a 128-dimensional directional gradient histogram descriptor vector.

[0037] The structure of the multi-scale matching fusion model is as follows: the input layer receives directional gradient histogram descriptors from different scale levels; the feature encoding layer contains three residual convolutional blocks, each consisting of two convolutional layers and a skip connection, with a kernel size of 3×3 and channel numbers of 64, 128, and 256 respectively; the multi-scale fusion layer adopts a feature pyramid network structure, fusing features of different scales through top-down paths and lateral connections; the matching head layer contains two fully connected layers, the first fully connected layer having 512 neurons, and the second fully connected layer outputting matching confidence and keypoint coordinate offset; the fusion weight parameters in the multi-scale fusion layer of the multi-scale matching fusion model are determined based on the geographic grid block index, the Gaussian difference response value of scale-invariant keypoints, and the gradient change rate of the dense disparity field. The fusion weight parameters are normalized by the softmax function and multiplied with the feature map of the corresponding scale to achieve adaptive fusion.

[0038] The multi-scale matching fusion model utilizes an imbalanced data processing mechanism based on adaptive sampling. This mechanism is implemented as follows: First, the sample quantity distribution of different matching difficulty categories in the training dataset is statistically analyzed, classifying matching difficulty into three categories: low-texture region matching, repetitive texture region matching, and normal texture region matching. Then, the sample difficulty coefficient for each matching difficulty category is calculated, determined by the ratio of the average matching error rate of that category to the global average matching error rate. Next, the sampling probability of each matching difficulty category is dynamically adjusted based on the sample difficulty coefficient. For matching difficulty categories with a sample difficulty coefficient greater than 1.5, the sampling probability is increased to 2 to 3 times the original value; for matching difficulty categories with a sample difficulty coefficient less than 0.5, the sampling probability is decreased to 0.5 to 0.7 times the original value. The increase and decrease factors of the sampling probability are obtained by combining factors that minimize the recall variance of each matching difficulty category through grid search on the validation dataset. Finally, samples are extracted in each training batch according to the adjusted sampling probability to construct the training batch.

[0039] The multi-scale matching fusion model optimizes the feature matching process using the principle of two dam water storage problems. The principle of two dam water storage problems is implemented as follows: the histogram of directional gradients (HCG) descriptor sequence of the left image is considered the boundary of the first container, and the histogram of directional gradients (HCG) descriptor sequence of the right image is considered the boundary of the second container; feature similarity is considered as the water storage height between the first and second containers; the first pointer is initialized to point to the starting position of the histogram of directional gradients (HCG) descriptor sequence of the left image, and the second pointer is initialized to point to the starting position of the histogram of directional gradients (HCG) descriptor sequence of the right image; the histogram of directional gradients (HCG) descriptor corresponding to the first pointer position and the histogram of directional gradients (HCG) corresponding to the second pointer position are calculated. The feature similarity of the descriptors is used as the current water level. If the Gaussian difference response value of the directional gradient histogram descriptor corresponding to the first pointer is less than the Gaussian difference response value of the directional gradient histogram descriptor corresponding to the second pointer, the first pointer is moved one position to the right; otherwise, the second pointer is moved one position to the right. During the pointer movement, the number of directional gradient histogram descriptor pairs with feature similarity higher than the matching threshold is accumulated as the effective matching capacity. The above pointer movement process is repeated until the first pointer reaches the end of the directional gradient histogram descriptor sequence of the left image or the second pointer reaches the end of the directional gradient histogram descriptor sequence of the right image, thus obtaining the set of directional gradient histogram descriptor pairs corresponding to the maximum effective matching capacity.

[0040] The adaptive sampling-based imbalanced data processing mechanism dynamically adjusts the sampling probability of different matching difficulty categories, enabling the multi-scale matching fusion model to pay more attention to difficult-to-match regions during training. This avoids the multi-scale matching fusion model overfitting simple samples while ignoring difficult samples, thereby improving the matching accuracy and robustness of the multi-scale matching fusion model in low-texture and repetitive texture regions, enhancing the adaptability of the multi-scale matching fusion model to different scenes, reducing the matching failure rate caused by sample imbalance, and improving the overall integrity of 3D reconstruction. The principle of the two dam water storage problems described above uses a dual-pointer strategy to efficiently search the space of histogram of directional gradients (HGP) descriptors, avoiding redundant calculations in traditional brute-force matching methods. This reduces the time complexity of the matching search from quadratic to linear. Furthermore, by using the concept of water storage height, feature similarity is transformed into a geometrically intuitive capacity metric. This allows the feature matching process to naturally select HGP descriptor pairs with high and uniform feature similarity, improving the spatial distribution rationality of matching pairs, reducing the false matching rate, and providing more reliable basic data for subsequent disparity map calculations and 3D point cloud reconstruction. This significantly improves the processing efficiency and matching quality of large-scale stereo image pair satellite image data matching.

[0041] The steps for establishing the training dataset of the multi-scale matching fusion model specifically include: collecting stereo image pairs from satellite image data of different geographical environments, including urban built-up areas, mountainous terrain areas, plain farmland areas, and water bodies; manually annotating the collected stereo image pairs, including the location coordinates and matching confidence labels of the matching keypoint pairs; calculating the disparity ground truth based on the annotated matching keypoint pairs, and using the disparity ground truth as a supervision signal; performing data augmentation on the annotated data, including random rotation, random scaling, random brightness adjustment, and random noise addition; and dividing the augmented annotated data into training, validation, and test sets in an 8:1:1 ratio.

[0042] The specific steps for training the multi-scale matching fusion model include: using the Adam optimizer for parameter optimization, setting the initial learning rate to 0.001, and decreasing the learning rate to 0.5 times the original value every 10 training epochs; the loss function consists of matching confidence loss and coordinate offset loss, with the matching confidence loss using the cross-entropy loss function and the coordinate offset loss using the smoothed L1 loss function, and the weight ratio of matching confidence loss to coordinate offset loss being 1:2; the training batch size is set to 32, and the total number of training epochs is set to 80 to 100; during the training process, the performance of the multi-scale matching fusion model is evaluated on the validation set every 5 training epochs, and the matching accuracy and average coordinate error on the validation set are recorded; training is stopped when the matching accuracy on the validation set does not improve for 10 consecutive training epochs, and the parameters of the multi-scale matching fusion model with the best performance on the validation set are selected as the final multi-scale matching fusion model.

[0043] The matching confidence is determined by calculating the ratio of the distance between the histogram descriptors of the directional gradients of corresponding point pairs in the corresponding point set to the distance between the histogram descriptors of the directional gradients of the suboptimal matching point. The matching confidence threshold is obtained by statistically analyzing the matching accuracy curves under different distance ratios on the validation dataset. The maximum distance ratio that makes the matching accuracy greater than 95% is selected as the matching confidence threshold. When the distance ratio is less than 0.8, the matching confidence is considered high, and when the distance ratio is greater than or equal to 0.8, the matching confidence is considered low.

[0044] The upper-layer reconstruction accuracy optimization model aims at the point cloud density and surface smoothness of the 3D point cloud reconstruction, while the lower-layer computing resource allocation model aims at the processing time and memory usage of the image processing unit. The upper-layer reconstruction accuracy optimization model and the lower-layer computing resource allocation model are coupled through the processing priority parameter of the image processing unit.

[0045] The objective function of the upper-layer reconstruction accuracy optimization model is used to maximize the reconstruction quality of the 3D point cloud. The inputs include point cloud density, surface normal vector consistency, and the number of matching pairs, and the output is the reconstruction accuracy score. The objective function of the upper-layer reconstruction accuracy optimization model is expressed as follows: the product of the normalized value of point cloud density and the normalized value of surface normal vector consistency is used as the basic score term, and the square root of the normalized value of the number of matching pairs is used as the matching quality correction term. The basic score term and the matching quality correction term are added together and multiplied by the spatial weight factor of the image processing unit to obtain the reconstruction accuracy score. The constraints of the upper-layer reconstruction accuracy optimization model are that the normalized value of point cloud density must be greater than 0.3, the normalized value of surface normal vector consistency must be greater than 0.5, and the normalized value of the number of matching pairs must be greater than 0.4. The threshold values ​​of the normalized value of point cloud density, the normalized value of surface normal vector consistency, and the normalized value of the number of matching pairs are obtained by statistically analyzing the reconstruction quality distribution under different normalization values ​​on the historical reconstruction dataset. The minimum normalization value that makes the reconstruction quality score in the top 30% is selected as the threshold value of the constraints.

[0046] The objective function of the lower-level computing resource allocation model is to minimize the weighted sum of total processing time and memory usage. The inputs include the data volume, processing priority, and available computing resources of each image processing unit, and the output is a resource allocation scheme. The objective function of the lower-level computing resource allocation model is expressed as follows: the sum of the products of the processing time of each image processing unit and the reciprocal of its processing priority is used as the time term; the maximum ratio of the memory usage of each image processing unit to the available memory capacity is used as the memory term; and the weighted sum of the normalized values ​​of the time term and the memory term is used as the total cost. The weighting coefficients are 0.6 and 0.4. The constraints of the lower-level computing resource allocation model are that the total memory usage of all image processing units does not exceed 90% of the available memory capacity, and the processing time of a single image processing unit does not exceed 120% of the time limit.

[0047] The processing priority parameter serves as the output variable of the upper-level reconstruction accuracy optimization model and the input variable of the lower-level computing resource allocation model, forming a coupling relationship between the upper-level reconstruction accuracy optimization model and the lower-level computing resource allocation model. This coupling relationship achieves game equilibrium through iterative solution.

[0048] The camera parameters include a camera intrinsic matrix and a camera extrinsic matrix. The camera intrinsic matrix contains the focal length and principal point coordinates, while the camera extrinsic matrix contains the rotation matrix and translation vector. The camera parameters are obtained from the metadata of the stereo image pair satellite image data.

[0049] The disparity map is obtained by calculating the pixel coordinate difference between corresponding point pairs in the corresponding point set. The disparity value represents the positional offset of the same spatial point in the left and right images, and the range of the disparity value is from 0 to 20% of the image width.

[0050] The scale-space consistency constraint processing is performed by comparing the disparity value difference between adjacent scale levels. When the disparity value difference between adjacent scale levels is greater than the inter-scale difference threshold, the current disparity value is smoothed and corrected by filtering. The inter-scale difference threshold is determined based on the scale factor and image resolution.

[0051] The multi-scale disparity propagation is achieved by iteratively updating the dense disparity field. The initial disparity field is generated by interpolation of the corresponding point set. In each iteration, for each pixel in the dense disparity field, neighboring pixels with similar disparity values ​​are searched in the neighborhood of the pixel. The disparity values ​​of the neighboring pixels are weighted and averaged to update the disparity value of the pixel. The number of iterations is 5 to 10.

[0052] The 3D point cloud reconstruction calculates the 3D coordinates of spatial points based on the geometric relationship between the dense parallax field and camera parameters. The geometric relationship is established through the principle of similar triangles, and the depth value of the 3D coordinates is inversely proportional to the parallax value.

[0053] The registration and fusion is achieved by calculating the rigid transformation matrix between the 3D point clouds reconstructed by adjacent image processing units. The rigid transformation matrix includes a rotation matrix and a translation vector. The rigid transformation matrix is ​​solved by minimizing the spatial distance error between corresponding point pairs.

[0054] The moving cube algorithm establishes a regular voxel mesh in three-dimensional space, determines whether each voxel vertex is located inside or outside the three-dimensional surface, and generates triangular patches based on the inside and outside markings of the voxel vertices. The vertex positions of the triangular patches are determined by linear interpolation.

[0055] The texture mapping maps the color information of the stereo image pair satellite image data onto the surface of the three-dimensional mesh model. The corresponding texture color is obtained by calculating the projection coordinates of the vertices of the three-dimensional mesh model on the image plane. The illumination consistency processing adjusts the brightness and color of the stereo image pair satellite image data from different perspectives to ensure a smooth transition of texture color at the stitching point.

[0056] The specific implementation methods of the above steps are described in detail below.

[0057] The specific implementation of step S1 involves first acquiring stereo image pairs of the target area via a satellite remote sensing platform. This data includes multiple images of the same area taken from different perspectives. Next, radiometric correction is performed on the acquired raw image data to eliminate radiometric distortion caused by atmospheric scattering, uneven sensor response, and other factors, ensuring that the image's radiometric values ​​accurately reflect the reflectivity of ground features. Then, geometric correction is performed, establishing a mapping relationship between image coordinates and geographic coordinates using ground control points and satellite orbit parameters, eliminating geometric distortion caused by satellite attitude changes, Earth curvature, and terrain undulations. After correction, a geographic grid block index is established. This index divides the image data into regular grids according to latitude and longitude intervals, with the latitude and longitude dimensions of a single grid set between 0.5 and 2 degrees. This range balances computational efficiency and reconstruction accuracy. Finally, based on the established geographic grid block index, the complete stereo image pair data is divided into several independent image processing units, each corresponding to a geographic grid, facilitating subsequent parallel processing and resource scheduling.

[0058] The specific implementation of step S2 involves constructing a Gaussian pyramid structure for each image processing unit. This structure is generated by repeatedly performing Gaussian filtering and downsampling operations on the original image, resulting in 4 to 6 scale levels. The scale factor used between adjacent levels is set between 1.414 and 1.732. This factor range has been experimentally verified to maintain a repeat detection rate of over 85% for scale-invariant keypoints, ensuring the stability of feature extraction. At each level of the constructed Gaussian pyramid, scale-invariant keypoints are extracted by calculating local extrema in the scale space. Specifically, the Gaussian difference response value is calculated for each candidate pixel, reflecting the saliency of the pixel at the current scale. The response value is compared with a threshold adaptively determined based on the image contrast of the current level, and candidate points with response values ​​exceeding the threshold are retained as the final scale-invariant keypoints. For each extracted keypoint, the distribution characteristics of gradient directions are statistically analyzed within its neighborhood. The neighborhood is divided into 16 sub-regions of 4×4. The gradient intensity of 8 main directions is calculated in each sub-region and a histogram is constructed. Finally, a 128-dimensional directional gradient histogram descriptor vector is formed. This descriptor has rotation invariance and illumination robustness and can effectively describe the local texture features of the keypoint.

[0059] The specific implementation of step S3 involves inputting the multi-scale feature sets extracted by each image processing unit into a multi-scale matching fusion model for feature matching processing. This model employs a deep neural network architecture. The input layer receives histogram descriptors of oriented gradients at different scales. The feature encoding layer contains three residual convolutional blocks for extracting deep features. Each residual convolutional block consists of two 3×3 convolutional layers and skip connections, with channel numbers of 64, 128, and 256 respectively. The multi-scale fusion layer uses a feature pyramid network structure, achieving the fusion of features at different scales through top-down paths and lateral connections. The fusion weight parameters are adaptively determined based on the geographic grid block index, the Gaussian difference response value of the keypoints, and the gradient change rate of the dense disparity field. The weights are normalized using softmax and then multiplied by the corresponding feature map to complete the fusion. The matching head layer contains two fully connected layers. The first layer has 512 neurons, and the second layer outputs the matching confidence and keypoint coordinate offset. The model utilizes an imbalanced data processing mechanism based on adaptive sampling. First, it statistically analyzes the distribution of samples in the training set for three categories: low-texture region matching, repetitive texture region matching, and normal texture region matching. It then calculates the difficulty coefficient for each category, which is the ratio of the category's average matching error rate to the global average error rate. For difficult categories with a difficulty coefficient greater than 1.5, the sampling probability is increased by 2 to 3 times; for easy categories with a difficulty coefficient less than 0.5, the sampling probability is reduced to 0.5 to 0.7 times. The specific multiplier is determined by minimizing the variance of the recall rate for each category through grid search. The model also employs a two-pointer strategy to optimize the matching process, treating the descriptor sequences of the left and right images as the boundaries of two containers and the feature similarity as the water level. By moving the pointer, the number of descriptor pairs with feature similarity higher than the matching threshold is accumulated as the effective matching capacity. This method reduces the time complexity of the matching search from quadratic to linear, significantly improving processing efficiency. After obtaining the corresponding point set, the ratio of the descriptor distance of each pair of matching points to the distance of the second-best matching point is calculated as the matching confidence. A match is considered reliable when the ratio is less than 0.8. If the overall matching confidence is lower than the preset threshold, the scale level range of feature extraction is adjusted and step S2 is re-executed. The threshold is determined by selecting the maximum distance ratio that makes the matching accuracy greater than 95% through validation set statistics.

[0060] The specific implementation of step S4 involves constructing a two-layer game model consisting of an upper-layer reconstruction accuracy optimization model and a lower-layer computational resource allocation model to schedule resources for the reconstruction task of the image processing unit. The upper-layer model aims to maximize the quality of 3D point cloud reconstruction. Inputs include point cloud density, surface normal vector consistency, and the number of matching pairs. The product of the normalized point cloud density value and the normalized surface normal vector consistency value is used as the basic scoring term, and the square root of the normalized number of matching pairs is used as the matching quality correction term. The sum of these two values ​​is multiplied by the spatial weight factor of the image processing unit to obtain the reconstruction accuracy score. Constraints require that the normalized point cloud density value be greater than 0.3, the normalized surface normal vector consistency value be greater than 0.5, and the normalized number of matching pairs be greater than 0.4. These thresholds are determined by statistically selecting the minimum values ​​that place the reconstruction quality score within the top 30% based on historical data. The lower-level model aims to minimize total processing time and memory usage. Inputs include the data volume, processing priority, and available computing resources for each image processing unit. The time term is calculated by summing the products of each unit's processing time and the inverse of its priority, while the memory term is calculated by maximizing the ratio of memory usage to available memory for each unit. These two terms are normalized and then weighted and summed with coefficients of 0.6 and 0.4 to obtain the total cost. Constraints require that the total memory usage of all units does not exceed 90% of available memory, and the processing time of a single unit does not exceed 120% of the time limit. The processing priority parameter serves as both the output variable of the upper-level model and the input variable of the lower-level model, forming a coupling relationship between the two models. Game equilibrium is achieved through iterative solution, ultimately determining the processing priority and computing resource allocation scheme for each image processing unit.

[0061] The specific implementation of step S5 involves first obtaining camera parameters from the metadata of the stereo image pair data, including a camera intrinsic parameter matrix containing focal length and principal point coordinates, and a camera extrinsic parameter matrix containing rotation matrix and translation vector. Using these camera parameters and the corresponding point set obtained in step S3, a disparity map is calculated. The disparity value represents the difference in pixel coordinates of the same spatial point in the left and right images, typically ranging from 0 to 20% of the image width. The calculated disparity map is then subjected to scale-space consistency constraint processing. This processing is achieved by comparing the disparity value differences between adjacent scale levels. When the disparity value difference between adjacent levels exceeds the inter-scale difference threshold determined based on the scale factor and image resolution, a smoothing filter correction is performed on the current disparity value to eliminate inter-scale inconsistencies. Next, a dense disparity field is obtained through multi-scale disparity propagation. The initial disparity field is generated by interpolation from a sparse corresponding point set, and then it is updated iteratively 5 to 10 times. In each iteration, for each pixel in the disparity field, neighboring pixels with similar disparity values ​​are searched. The disparity values ​​of these neighboring pixels are then weighted and averaged to update the disparity value of the pixel, gradually spreading the sparse disparity to the entire image to form a dense disparity field. Finally, the 3D coordinates of the spatial points are calculated based on the geometric relationship between the dense disparity field and camera parameters. This geometric relationship is based on the principle of similar triangles, where the depth value of the 3D coordinates is inversely proportional to the disparity value. Through this relationship, each pixel on the 2D image is converted into point cloud coordinates in 3D space, completing the 3D point cloud reconstruction.

[0062] The specific implementation of step S6 involves first performing registration and fusion processing on the 3D point clouds reconstructed by each image processing unit. This processing is achieved by calculating the rigid transformation matrix between adjacent unit point clouds. The rigid transformation matrix includes a rotation matrix and a translation vector, which is obtained by minimizing the spatial distance error of corresponding point pairs, ensuring that adjacent point clouds are accurately aligned in a unified coordinate system. After registration, the point clouds of all units are fused into a complete 3D point cloud model. Next, a 3D mesh model is generated using the moving cube algorithm. This algorithm establishes a regular voxel mesh in 3D space, determines whether the eight vertices of each voxel are located inside or outside the 3D surface, and generates corresponding triangular faces by looking up a predefined triangular facet configuration table based on the inside and outside markings of the vertices. The vertex positions of the triangular faces are precisely determined by linear interpolation, ultimately forming a continuous mesh surface composed of triangular faces. Texture mapping is then performed on the generated 3D mesh model, mapping the color information of the stereo image data to the mesh surface. The texture color of the corresponding pixel is obtained by calculating the projection coordinates of the 3D mesh vertices on the image plane, giving the model a realistic appearance. Finally, illumination consistency processing is performed. This process adjusts the brightness and color of image data from different perspectives to ensure a smooth transition of texture colors at mesh joints, eliminating color differences and seams caused by variations in lighting conditions, and outputting the final realistic 3D model of the target area.

[0063] It should be noted that the key technical ideas of this invention include an imbalanced data processing mechanism based on adaptive sampling, a two-pointer matching strategy for the water storage problem of two dams, and a resource scheduling mechanism based on a two-layer game model.

[0064] An imbalanced data processing mechanism based on adaptive sampling dynamically adjusts the sampling probabilities of different matching difficulty categories, allowing the model to pay more attention to difficult regions such as low-texture and repetitive textures during training, avoiding overfitting to easy samples while ignoring difficult ones. Compared to traditional balanced sampling methods, this mechanism significantly improves the matching accuracy and robustness of the model in complex scenes, reduces matching failures caused by uneven sample distribution, and enhances the integrity and reliability of the overall 3D reconstruction.

[0065] The dual-pointer matching strategy for two dam water storage problems transforms feature matching into a container water storage problem. By using two pointers to efficiently search the descriptor sequence space, it avoids the redundant computation of traditional brute-force matching. This strategy reduces the time complexity of the matching search from quadratic to linear. Simultaneously, by utilizing the concept of water storage height, it naturally selects matching pairs with high feature similarity and uniform distribution, improving the spatial distribution rationality of matching pairs, reducing the false matching rate, and providing more reliable basic data for subsequent disparity calculation and point cloud reconstruction.

[0066] The resource scheduling mechanism of the two-layer game theory model achieves optimal allocation of computing resources while ensuring reconstruction quality through a coupled game between upper-layer reconstruction accuracy optimization and lower-layer computational resource allocation. Compared with the traditional fixed-priority scheduling method, this mechanism can dynamically adjust the processing order and resource allocation according to the reconstruction difficulty and resource requirements of each image processing unit, balancing the contradiction between reconstruction accuracy and computational efficiency, and significantly improving the overall processing efficiency of large-scale regional 3D reconstruction.

[0067] The synergistic effect of three key technological approaches forms a complete technical system. The adaptive sampling mechanism improves the training quality of the feature matching model, enabling it to accurately identify feature points in various scenes. The dual-pointer matching strategy efficiently performs feature matching during the inference phase, quickly generating high-quality corresponding point sets. The two-layer game theory model optimizes resource scheduling at the global level, ensuring that each processing unit receives reasonable computing resources. These three elements work together to form a closed-loop optimization from data training and feature matching to resource scheduling, achieving significant improvements in processing efficiency, matching accuracy, and resource utilization compared to traditional methods. This makes it particularly suitable for high-precision real-world 3D modeling tasks covering large areas.

[0068] It should be noted that this invention also solves the following technical problem: how to balance the reconstruction accuracy requirements of different image processing units under limited computing resources during the 3D reconstruction of large-scale stereo image pairs from satellite images. Traditional methods typically employ uniform allocation strategies or simple priority sorting based on data volume, failing to consider the differences in contribution of different regions to the final reconstruction quality. This results in wasted computing resources in areas with minimal impact on overall accuracy, while complex terrain areas and densely built-up areas, which play a crucial role in reconstruction quality, suffer from insufficient resources, affecting accuracy. This invention constructs a two-layer game model. The upper layer establishes a reconstruction accuracy scoring function with point cloud density and surface smoothness as objectives, determining accuracy requirements based on the spatial weight factors and matching quality of each image processing unit. The lower layer establishes a resource allocation optimization function with processing time and memory usage as objectives. The two layers are coupled through processing priority parameters to form a Stackelberg game. By iteratively solving the game equilibrium, computing resources are allocated differentially according to accuracy requirements, ensuring that key areas receive sufficient resources to guarantee reconstruction quality. Simultaneously, all areas are processed under constraints, achieving optimal reconstruction quality under limited resources. Furthermore, this invention also addresses the technical problem of low feature matching accuracy in low-texture and repetitive texture regions. By employing an imbalanced data processing mechanism based on adaptive sampling, the sampling probability of different matching difficulty categories is dynamically adjusted according to the sample difficulty coefficient. This allows the multi-scale matching fusion model to pay more attention to low-texture regions such as water bodies and plains, and repetitive texture regions such as building clusters during training. This avoids the model overfitting simple samples while ignoring difficult samples, thereby improving the matching accuracy and robustness of the model in complex scenes, reducing matching failures caused by insufficient feature discrimination, and enhancing the integrity of 3D reconstruction.

[0069] A second aspect of the present invention provides a computer-readable storage medium storing program instructions that, when executed in a computer, are used to perform the above-described method for real-scene 3D modeling based on stereo image pair satellite data.

[0070] A third aspect of the present invention provides a real-scene 3D modeling system based on stereo image pair satellite data, comprising the aforementioned computer-readable storage medium. The system is any one of a computer, a server, or a microcontroller. The computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor that executes the program instructions stored in the computer-readable storage medium.

[0071] Specifically, the principle of this invention is as follows: The invention solves this technical problem by improving processing efficiency through algorithm-level complexity optimization and system-level resource scheduling. First, the principle of solving the two dam water storage problem breaks through the computational bottleneck of traditional matching methods. This method treats the feature descriptor sequences of the left and right images as the boundaries of two containers, and the feature similarity as the water storage height. It initializes two pointers to point to the starting positions of the two sequences, and determines the pointer movement direction based on the relationship between the Gaussian difference response values. During the movement, descriptor pairs with feature similarity higher than the matching threshold are accumulated as effective matching capacity. This strategy allows the algorithm to find the set of descriptor pairs corresponding to the maximum effective matching capacity with only one linear traversal. Compared to traditional methods that require comparing all possible feature point pairs, this method reduces the time complexity from quadratic to linear, directly translating into a significant reduction in processing time. Second, the geometric concept of water storage height not only simplifies similarity calculation but also naturally filters out spatially evenly distributed matching pairs by maximizing capacity, avoiding excessive concentration of matching points in local areas and providing more reliable basic data for subsequent disparity calculation. Furthermore, the multi-scale feature extraction and matching fusion mechanism obtains image representations at different scales by constructing Gaussian pyramids. The feature pyramid network integrates multi-scale information to enhance feature description capabilities, enabling the matching process to run in parallel at each scale level and fully utilizing multi-core computing resources. Finally, the two-layer game theory model models the reconstruction task of the image processing unit as a decision problem between upper and lower layers. The upper layer evaluates the reconstruction accuracy requirements based on point cloud density and surface smoothness, while the lower layer optimizes resource allocation based on processing time and memory usage. The two layers are coupled through processing priority parameters and iteratively solve for the Nash equilibrium, ensuring that computing resources are prioritized for key areas, avoiding resource waste, and achieving a systematic improvement in overall matching efficiency. These techniques work synergistically from both algorithm optimization and resource management perspectives, aligning with the technical logic of efficiency optimization based on the essence of computational complexity.

[0072] The following provides a specific embodiment 1 of the present invention. The specific implementation of step S1 in this embodiment 1 is the same as that described above, and will not be repeated in detail here. The specific implementation of other steps is described in detail below.

[0073] Step S2 involves constructing a Gaussian pyramid structure for each image processing unit. This structure is generated by repeatedly performing Gaussian filtering and downsampling operations on the original image, resulting in 4 to 6 scale levels. The scale factor between adjacent levels is set between 1.414 and 1.732. This factor range has been experimentally verified to maintain a repeat detection rate of over 85% for scale-invariant keypoints, ensuring the stability of feature extraction. Scale-invariant keypoints are extracted at each level of the constructed Gaussian pyramid by calculating local extrema in the scale space. Specifically, the Gaussian difference response value is calculated for each candidate pixel, and the formula for this response value is as follows:

[0074] ;

[0075] In the formula, This is the normalized Gaussian difference response value. and The spatial coordinates of the candidate pixels are in pixels. The current scale parameter is in pixels. The scale factor is set to 1.6 by default. This represents the Gaussian filter response value at a larger scale. This represents the Gaussian filter response value at the current scale. The maximum Gaussian response value at the current level is used for normalization. This formula detects local extrema of an image by calculating the difference in Gaussian filter responses at different scales. The normalization process eliminates the dimensional differences between scales, making the response values ​​at different scale levels comparable.

[0076] The formula for determining the adaptive threshold is expressed as follows:

[0077] ;

[0078] In the formula, For the first Layer-scale adaptive threshold The baseline threshold is empirically set at 0.03. The scale level index has a value range of 1 to 6. For the first Contrast of layered images, The reference contrast is typically taken as the median of the contrast across all levels. This formula dynamically adjusts the threshold based on the image contrast of the current level, increasing the threshold for high-contrast areas to reduce redundant feature points, and decreasing the threshold for low-contrast areas to ensure a sufficient number of feature points.

[0079] For each extracted keypoint, the distribution characteristics of gradient directions are statistically analyzed within its neighborhood. The neighborhood is divided into 16 sub-regions of 4×4. The gradient intensity of 8 main directions is calculated in each sub-region and a histogram is constructed. Finally, a 128-dimensional directional gradient histogram descriptor vector is formed. This descriptor has rotation invariance and illumination robustness and can effectively describe the local texture features of the keypoint.

[0080] The specific implementation of step S3 involves inputting the multi-scale feature sets extracted by each image processing unit into a multi-scale matching fusion model for feature matching processing. This model employs a deep neural network architecture. The input layer receives histogram descriptors of oriented gradients at different scales. The feature encoding layer contains three residual convolutional blocks for extracting deep features. Each residual convolutional block consists of two 3×3 convolutional layers and skip connections, with channel numbers of 64, 128, and 256 respectively. The multi-scale fusion layer uses a feature pyramid network structure, achieving the fusion of features at different scales through top-down paths and lateral connections. The formula for determining the fusion weight parameters is as follows:

[0081] ;

[0082] In the formula, For the first Normalized fusion weights at each scale level, For scale-level indexing, For the first The normalized coded value of the geographic grid block index corresponding to each scale level ranges from 0 to 1. For the first Normalized average Gaussian difference response values ​​of key points at each scale level For the first Normalized average gradient rate of change of dense parallax fields at each scale level , , The first The geographic index weight coefficients, response value weight coefficients, and gradient weight coefficients at each scale level are automatically learned through model training. This represents the total number of scale levels. It is a natural exponential function. The summation index variable is used. This formula uses the softmax function to transform the feature importance at multiple scale levels into normalized weights. The weight calculation comprehensively considers geographical location, feature saliency, and disparity gradient information, enabling the model to adaptively focus on features at different scales.

[0083] The matching head layer consists of two fully connected layers. The first layer has 512 neurons, and the second layer outputs the matching confidence score and keypoint coordinate offset. After obtaining the corresponding point set, the formula for calculating the matching confidence score for each pair of matching points is as follows:

[0084] ;

[0085] In the formula, To match confidence levels, The directional gradient histogram descriptor for the best-matching point pair is the Euclidean distance. The directional gradient histogram descriptor for suboptimal matching points is the Euclidean distance. The maximum possible distance describing a subspace is typically set to a value used for normalization. Because the descriptor is 128-dimensional. This formula evaluates the reliability of the match by comparing the distance ratio between the best match and the second-best match; the smaller the ratio, the more significant the advantage of the best match over the second-best match, and the more reliable the match. When A score less than 0.8 is considered to indicate a high confidence level in the match. A score greater than or equal to 0.8 is considered to indicate a low confidence level in the match.

[0086] The specific implementation of step S4 involves constructing a two-layer game model consisting of an upper-layer reconstruction accuracy optimization model and a lower-layer computing resource allocation model to schedule resources for the reconstruction task of the image processing unit. The objective function formula of the upper-layer reconstruction accuracy optimization model is expressed as follows:

[0087] ;

[0088] In the formula, For the first Reconstruction accuracy score for each image processing unit This is the index number of the image processing unit. For the first Normalized point cloud density values ​​of each image processing unit For the first The surface normal vector consistency normalization value of each image processing unit For the first Normalized value of the number of matching pairs per image processing unit For the first The spatial weighting factor for each image processing unit is determined based on geographic importance and ranges from 0.5 to 2.0. The formula for calculating the normalized point cloud density is as follows:

[0089] ;

[0090] In the formula, For the first The actual point cloud density of each image processing unit is expressed in points per square meter. The reference point cloud density is typically taken as 100 points per square meter. The formula for calculating the surface normal vector uniformity normalization value is as follows:

[0091] ;

[0092] In the formula, For the first The number of point clouds per image processing unit For the first The first image processing unit The angle between the normal vectors of a point cloud and its neighboring point clouds is expressed in radians. The value of the cosine of the included angle. The point index within the point cloud has a value range of 1 to 1. The formula for calculating the normalized value of the number of matching pairs is as follows:

[0093] ;

[0094] In the formula, For the first The actual number of matching pairs per image processing unit For the first The theoretical maximum number of matching pairs per image processing unit is estimated by the square root of the number of keypoints. The upper-level objective function is determined by the basic scoring term. Reflecting the density and surface quality of the point cloud, by matching quality correction terms. Considering the impact of the number of matching pairs on reconstruction quality, the square root form avoids the number of matching pairs excessively dominating the scoring, and the spatial weighting factor... This reflects the differences in importance between different regions.

[0095] The objective function formula for the lower-level computing resource allocation model is expressed as follows:

[0096] ;

[0097] In the formula, For the total cost of normalization, This represents the total number of image processing units. For the first The processing time of each image processing unit is measured in seconds. For the first The processing priority of each image processing unit ranges from 1 to 10. The maximum processing time for a single image processing unit is in seconds. For the first The memory footprint of each image processing unit is measured in megabytes. The available memory capacity is in megabytes (MB), and 0.6 and 0.4 are weighting coefficients for the time and memory terms, respectively. The time term in this formula reflects the principle of allocating resources to high-priority tasks first by weighting the processing time with the inverse of priority. The memory term reflects resource bottleneck constraints by using the maximum memory usage ratio. The combined cost is obtained by normalizing and then weighting the sum of the two terms. (Processing priority parameter) As the output variables of the upper-level model and the input variables of the lower-level model, a coupling relationship is formed between the two-level models, and the game equilibrium is achieved through iterative solution.

[0098] The specific implementation of step S5 involves first obtaining camera parameters from the metadata of the stereo image pair data, including the camera intrinsic parameter matrix containing focal length and principal point coordinates, and the camera extrinsic parameter matrix containing rotation matrix and translation vector. Using these camera parameters and the corresponding point set obtained in step S3, a disparity map is calculated. The disparity value calculation formula is expressed as follows:

[0099] ;

[0100] In the formula, Coordinates in the left image Normalized disparity value at that location, The x-coordinate of the matching point in the left image is in pixels. The x-coordinate of the corresponding matching point in the right image is in pixels. The ordinate of the matching point in the left image is in pixels. The image width is expressed in pixels for normalization. This formula calculates the parallax by measuring the difference in the horizontal coordinates of corresponding points in the left and right images, and the normalization process eliminates the influence of image size.

[0101] The formula for determining the inter-scale difference threshold in scale-space consistency constraint processing is as follows:

[0102] ;

[0103] In the formula, For the first Inter-scale difference threshold at the layer scale The baseline difference coefficient has an empirical value of 0.05. The scaling factor ranges from 1.414 to 1.732. The current image resolution is expressed in meters per pixel. The reference image resolution is typically set to 1 meter per pixel. This formula determines the allowable range of disparity values ​​between scales based on the scale factor and image resolution. When the disparity value difference between adjacent levels exceeds this threshold, a smoothing filter correction is performed.

[0104] The formula for updating disparity values ​​in multi-scale disparity propagation is expressed as follows:

[0105] ;

[0106] In the formula, For the first After the second iteration, the pixel The normalized disparity value, For the first After the second iteration, the neighboring pixels The normalized disparity value, For pixels The set of neighboring pixels, For neighboring pixels For the center pixel The weight, The iteration number index ranges from 1 to 10. The formula for calculating the neighborhood weight is as follows:

[0107] ;

[0108] In the formula, The disparity similarity control parameter is set to 0.1 by default. The disparity value update formula updates the disparity value of the current pixel by weighted averaging the disparity values ​​of neighboring pixels. The weights are determined based on disparity similarity, with neighboring pixels having greater weights, gradually spreading sparse disparity to the entire image to form a dense disparity field.

[0109] The formula for calculating the depth value in 3D point cloud reconstruction is as follows:

[0110] ;

[0111] In the formula, The depth value of a spatial point is in meters. The camera focal length is expressed in pixels. The baseline distance for the stereo image pair, in meters, is extracted from the translation vector of the camera extrinsic matrix. The absolute disparity value is expressed in pixels. The calculation yielded the result. This formula establishes an inverse relationship between disparity and depth based on the principle of similar triangles, and calculates the three-dimensional coordinates of a point in space using known camera parameters and disparity values.

[0112] The specific implementation of step S6 involves first performing registration and fusion processing on the 3D point clouds reconstructed by each image processing unit. This processing is achieved by calculating the rigid transformation matrix between adjacent unit point clouds. The rigid transformation matrix includes a rotation matrix and a translation vector, which is obtained by minimizing the spatial distance error of corresponding point pairs, ensuring that adjacent point clouds are accurately aligned in a unified coordinate system. After registration, the point clouds of all units are fused into a complete 3D point cloud model. Next, a 3D mesh model is generated using the moving cube algorithm. This algorithm establishes a regular voxel mesh in 3D space, determines whether the eight vertices of each voxel are located inside or outside the 3D surface, and generates corresponding triangular faces by looking up a predefined triangular facet configuration table based on the inside and outside markings of the vertices. The vertex positions of the triangular faces are precisely determined by linear interpolation, ultimately forming a continuous mesh surface composed of triangular faces. Texture mapping is then performed on the generated 3D mesh model, mapping the color information of the stereo image data to the mesh surface. The texture color of the corresponding pixel is obtained by calculating the projection coordinates of the 3D mesh vertices on the image plane, giving the model a realistic appearance. Finally, illumination consistency processing is performed. This process adjusts the brightness and color of image data from different perspectives to ensure a smooth transition of texture colors at mesh joints, eliminating color differences and seams caused by variations in lighting conditions, and outputting the final realistic 3D model of the target area.

[0113] It needs to be explained that the Gaussian difference response value formula detects image feature points by calculating the difference of Gaussian filtering at different scales. Normalization ensures that the response values ​​at different scales have uniform dimensions and comparability, avoiding instability in feature point detection caused by scale differences. This improves the robustness and consistency of feature extraction, providing a high-quality keypoint set for subsequent feature matching. The adaptive threshold formula dynamically adjusts the feature point selection threshold based on the contrast of the current image level. This allows the algorithm to adaptively balance the quantity and quality of feature points under different contrast environments. For high-contrast regions, increasing the threshold reduces redundant feature points and lowers the computational burden; for low-contrast regions, decreasing the threshold ensures a sufficient number of feature points to maintain reconstruction integrity, significantly improving the algorithm's adaptability in complex scenes. The surface normal vector consistency normalization value formula quantifies surface smoothness by calculating the average of the cosine values ​​of the angle between the normal vectors of adjacent point clouds. Its functional form is:

[0114] ;

[0115] This formula transforms the angle into a similarity measure between 0 and 1 using a cosine function. The averaging operation integrates the consistency information of the normal vectors of all point pairs; a value closer to 1 indicates a smoother surface, while a smaller value indicates a rougher surface or more noise. The formula for normalizing the number of matching pairs achieves normalization by dividing the actual number of matching pairs by the theoretical maximum value, eliminating the influence of differences in the number of feature points in different regions and making the matching quality comparable. The objective function of the upper-layer reconstruction accuracy optimization model comprehensively evaluates the impact of point cloud density, surface quality, and the number of matching pairs on reconstruction accuracy through a combination of basic scoring terms and matching quality correction terms. The introduction of spatial weight factors reflects the differences in importance of different regions. This objective function provides a clear quality optimization objective for resource scheduling, guiding the system to prioritize regions with high reconstruction accuracy requirements, thereby improving the overall quality and integrity of the 3D model. The objective function of the lower-level computational resource allocation model balances the dual constraints of processing efficiency and resource consumption through a weighted combination of time and memory terms. The inverse priority weighting mechanism ensures that high-priority tasks receive resource allocation first, and the maximum memory consumption ratio reflects the impact of resource bottlenecks on system performance. This objective function achieves task scheduling optimization under limited computational resource constraints, significantly improving the efficiency of large-scale stereo image data processing. The disparity value calculation formula eliminates the influence of image size differences through normalization, ensuring that disparity values ​​of images of different sizes have uniform dimensions and comparability, simplifying the subsequent depth value calculation and 3D reconstruction data processing flow. The inter-scale difference threshold formula dynamically determines the tolerance range of disparity consistency constraints based on the scale factor and image resolution, enabling the algorithm to adaptively perform scale-space consistency checks according to the characteristics of different scale levels. This effectively eliminates the disparity discontinuity problem caused by scale changes, improving the smoothness and reliability of dense disparity fields.

[0116] To better understand and implement this invention, the following is a specific application scenario of the invention, Example 2: A mountainous urban border area contains complex topography, including densely built-up urban areas, hilly and mountainous regions, and some plains and farmland, with a total area of ​​approximately 285 square kilometers. Traditional aerial photogrammetry methods in this region are hampered by cloud cover and topographical undulations, making it difficult to obtain complete 3D models. The technical team decided to use stereo imagery to create realistic 3D models from satellite data to solve the challenge of 3D reconstruction of large-scale, complex terrain.

[0117] The technical team first acquired stereo satellite imagery data of the target area, with a resolution of 0.5m and a coverage area spanning 1.8 degrees longitude and 1.6 degrees latitude. Radiometric correction was performed on the acquired stereo satellite imagery data to eliminate radiometric distortion caused by atmospheric scattering and uneven sensor response. Then, geometric correction was performed to correct geometric deformation caused by satellite attitude deviation and Earth curvature. A geographic grid index was established for the corrected imagery data, dividing the entire area into 6 image processing units with a grid number and geographic coordinate range, spaced 1.2 degrees apart in latitude and longitude.

[0118] For each image processing unit, the technical team constructed a 5-layer Gaussian pyramid for multi-scale decomposition. Layer 0 represents the original image with a resolution of 0.5m; Layer 1, after Gaussian filtering and downsampling, has a resolution of 0.71m; Layer 2 has a resolution of 1.01m; Layer 3 has a resolution of 1.43m; and Layer 4 has a resolution of 2.03m. The scale factor between adjacent layers was set to 1.414, which resulted in an 87% repetition detection rate for scale-invariant keypoints. Scale-invariant keypoints were extracted at each scale level, and the Gaussian difference response value was calculated for each candidate point by calculating local extrema in the scale space. 18,542 candidate points were detected at Layer 0, and 12,368 scale-invariant keypoints were retained after adaptive thresholding. The adaptive threshold was determined based on the image contrast at the current scale level; the adaptive threshold for Layer 0 was 0.028, and the adaptive threshold for Layer 4 was reduced to 0.015.

[0119] For the extracted scale-invariant keypoints, histogram-oriented gradient (HAR) descriptors were calculated. The neighborhood of each keypoint was divided into 4×4 sub-regions, and histograms of gradients in eight directions were calculated within each sub-region, forming a 128-dimensional HAR descriptor vector. A total of 74,215 scale-invariant keypoints were extracted from the six image processing units, establishing a multi-scale feature set. (For example...) Figure 1 As shown, the distribution of feature points at different scale levels exhibits a clear hierarchical structure.

[0120] The technical team inputs a multi-scale feature set into a pre-trained multi-scale matching fusion model for feature matching. The model's training dataset contains stereo image pairs from 15 different geographical environments. After manual annotation and data augmentation, the training set contains 58,000 samples, the validation set contains 7,250 samples, and the test set contains 7,250 samples. The training process uses the Adam optimizer with an initial learning rate of 0.001, and the model converges after 92 epochs. The model achieves a matching accuracy of 96.8% on the validation set, with an average coordinate error of 0.73 pixels.

[0121] During feature matching, the model employs an imbalanced data processing mechanism based on adaptive sampling. The technical team statistically analyzed the sample size distribution for different matching difficulty categories in the training dataset, as shown in Table 1.

[0122] Table 1. Sample distribution of matching difficulty categories

[0123]

[0124] The sampling probability of each matching difficulty category is dynamically adjusted based on the sample difficulty coefficient. The sampling probability of low-texture region matching and repetitive texture region matching is increased to 2.5 times and 2.2 times the original, respectively, while the sampling probability of normal texture region matching is decreased to 0.6 times the original. This adjustment allows the model to pay more attention to difficult-to-match regions during training, and the recall variance of each matching difficulty category decreases from 0.0382 to 0.0089.

[0125] The model also employs the principle of two dam water storage problems to optimize the feature matching process. The histogram of oriented gradients (HOPs) descriptor sequences of the left and right images are considered as the boundaries of two containers, respectively, and a two-pointer strategy is used to search for matching pairs. The first and second pointers are initialized to point to the starting positions of the two sequences, respectively. The pointers are moved based on the comparison of Gaussian difference response values, accumulating the number of HOP descriptor pairs with a feature similarity higher than 0.75. For example... Figure 2 As shown, this strategy reduces the time complexity of matching searches from Reduce to The matching time for processing 74,215 feature points was reduced from 438 seconds in the traditional method to 52 seconds.

[0126] After feature matching was completed, the technical team obtained the corresponding point sets between stereo image pairs. The sizes of the corresponding point sets for the six image processing units were 8542, 9136, 7823, 8965, 9247, and 8418 corresponding point pairs, respectively. Matching quality was judged based on the matching confidence of the corresponding point sets. The matching confidence was determined by calculating the ratio of the distance between the histogram descriptors of the directional gradients of the corresponding point pairs to the distance of the second-best matching point. A matching confidence threshold of 0.8 was set; a distance ratio less than 0.8 was considered to indicate high matching confidence. Statistical results showed that the average matching confidence of the six image processing units were 0.68, 0.72, 0.65, 0.71, 0.74, and 0.69, respectively, all higher than the matching confidence threshold, indicating that the matching quality met the requirements.

[0127] Next, the technical team utilized a two-layer game theory model to schedule resources for the reconstruction task of the image processing unit. The upper-layer reconstruction accuracy optimization model targets the point cloud density and surface smoothness of the reconstructed 3D point cloud, while the lower-layer computational resource allocation model targets the processing time and memory usage of the image processing unit. For example... Figure 3 As shown, by processing the priority parameter, the upper-level model is coupled with the lower-level model, and the game equilibrium is reached after 8 iterations.

[0128] The objective function of the upper-layer reconstruction accuracy optimization model calculates the reconstruction accuracy score. The inputs include point cloud density, surface normal vector consistency, and the number of matching pairs. The initial evaluation parameters for the six image processing units are shown in Table 2.

[0129] Table 2 Initial Evaluation Parameters of Image Processing Unit

[0130]

[0131] Based on the objective function calculation, the reconstruction accuracy scores for each image processing unit are 0.73, 0.82, 0.61, 0.78, 0.87, and 0.68, respectively. The lower-level computing resource allocation model determines a resource allocation scheme based on these scores and available computing resources, with 24 available computing cores and a total memory capacity of 128GB. The processing priorities output by the model, from highest to lowest, are: Image Processing Unit 5, Image Processing Unit 2, Image Processing Unit 4, Image Processing Unit 1, Image Processing Unit 6, and Image Processing Unit 3, with allocated computing cores of 5, 5, 4, 4, 3, and 3, respectively.

[0132] The disparity map is calculated based on the camera parameters of the corresponding point set and image processing unit. The camera parameters are obtained from the metadata of the stereo image pair satellite image data. The camera intrinsic matrix contains a focal length of 12850 pixels and the principal point coordinates are the image center. The camera extrinsic matrix contains a rotation matrix and a translation vector. The disparity map is obtained by calculating the pixel coordinate difference between corresponding point pairs, with the disparity value ranging from 0 to 18% of the image width. Scale-space consistency constraints are applied to the disparity map, comparing the disparity values ​​at adjacent scale levels. When the difference exceeds the inter-scale difference threshold of 3.2 pixels, the current disparity value is smoothed and corrected using a filtering process.

[0133] A dense disparity field is obtained through multi-scale disparity propagation, with the initial disparity field generated by interpolation of the corresponding point set. In each iteration, for each pixel in the dense disparity field, neighboring pixels with similar disparity values ​​are searched. The neighborhood search radius is set to 7 pixels, and the similarity threshold is set to 1.5 pixels. The disparity value of the current pixel is updated by weighted averaging of the disparity values ​​of the neighboring pixels, with the weights determined based on spatial distance and disparity value differences. After 7 iterations, the dense disparity field converges, and the average disparity values ​​of the 6 image processing units are 142.6, 158.3, 135.8, 151.4, 164.7, and 139.2 pixels, respectively.

[0134] 3D point cloud reconstruction was performed based on a dense parallax field and camera parameters. Geometric relationships were established using the principle of similar triangles, where the depth value of the 3D coordinates is inversely proportional to the parallax value. The 3D point cloud sizes reconstructed by the six image processing units were 28.47 million, 31.52 million, 26.34 million, 29.86 million, 33.18 million, and 27.65 million points, respectively. Figure 4 As shown, the reconstructed 3D point cloud accurately reproduces the topographic features of the mountainous city boundary area.

[0135] The reconstructed 3D point clouds from each image processing unit are registered and fused. The rigid transformation matrix between adjacent image processing units is calculated, including the rotation matrix and translation vector. The rigid transformation matrix is ​​solved by minimizing the spatial distance error between corresponding point pairs, with the registration error controlled within 0.85m. After registration and fusion, a complete 3D point cloud model containing 187 million points is obtained.

[0136] A 3D mesh model was generated using the moving cube algorithm, establishing a regular voxel mesh in 3D space with a voxel side length of 1.2m. The algorithm determined whether each voxel vertex was inside or outside the 3D surface, generating triangular faces based on the vertex's inside / outside markings. The vertex positions of these triangular faces were determined using linear interpolation. The generated 3D mesh model contained 65.47 million triangular faces, with an average mesh side length of 1.8m.

[0137] Texture mapping and lighting consistency processing are performed on the 3D mesh model to map the color information of the stereo image pair satellite image data onto the surface of the 3D mesh model. The texture color is obtained by calculating the projection coordinates of the vertices of the 3D mesh model onto the image plane; these projection coordinates are calculated based on camera parameters and the 3D coordinates of the vertices. Lighting consistency processing adjusts the brightness and color of the stereo image pair satellite image data from different viewpoints to ensure a smooth transition of texture colors at the stitching points, controlling color differences within 15 RGB color difference values.

[0138] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for real-scene 3D modeling based on stereo image pairs and satellite data, characterized in that, The process includes: acquiring stereo image pairs of satellite image data of the target area, performing radiometric and geometric corrections, establishing a geographic grid block index, and dividing the image processing unit according to the geographic grid block index; constructing a Gaussian pyramid for each image processing unit, performing multi-scale decomposition, extracting scale-invariant keypoints, and calculating directional gradient histogram descriptors to establish a multi-scale feature set; inputting the multi-scale feature set into a multi-scale matching and fusion model for feature matching to obtain the corresponding point set; the multi-scale matching and fusion model uses the principle of two dams impounding water to perform linear traversal in the directional gradient histogram descriptor sequence through a two-pointer strategy, guiding the pointer movement based on the Gaussian difference response value, and accumulating directional gradient histogram descriptor pairs with feature similarity higher than the matching threshold as effective matching capacity; using a two-layer game model to schedule resources for the image processing unit; calculating a disparity map based on the corresponding point set and camera parameters and performing scale-space consistency constraint processing; obtaining a dense disparity field through multi-scale disparity propagation for 3D point cloud reconstruction; registering and fusing the 3D point cloud to generate a 3D mesh model; performing texture mapping and illumination consistency processing to output a realistic 3D model. The principle of the two dam water storage problem is implemented by treating the histogram of directional gradients (HDRs) descriptor sequence of the left image as the boundary of the first container and the histogram of directional gradients (HDRs) descriptor sequence of the right image as the boundary of the second container. Feature similarity is regarded as the water storage height between the first and second containers. The first and second pointers are initialized to point to the starting positions of the histogram of directional gradients (HDRs) descriptor sequences of the left and right images, respectively. If the Gaussian difference response value of the histogram of directional gradients (HDR) descriptor corresponding to the first pointer is less than the Gaussian difference response value of the histogram of directional gradients (HDR) descriptor corresponding to the second pointer, the first pointer is moved one position to the right; otherwise, the second pointer is moved one position to the right. During the pointer movement, the number of histogram of directional gradients (HDRs) descriptor pairs with feature similarity higher than the matching threshold is accumulated as the effective matching capacity.

2. The method according to claim 1, characterized in that, The geographic grid block index is divided according to latitude and longitude intervals, and the size of a single grid ranges from 0.5 degrees to 2 degrees.

3. The method according to claim 2, characterized in that, The Gaussian pyramid is constructed by performing multiple Gaussian filters and downsampling on the image processing unit.

4. The method according to claim 3, characterized in that, The scale-invariant keypoints are obtained by calculating local extrema in the scale space. Gaussian difference response values ​​are calculated for each candidate point, and candidate points with Gaussian difference response values ​​greater than an adaptive threshold are retained as scale-invariant keypoints.

5. The method according to claim 4, characterized in that, The directional gradient histogram descriptor is constructed by the gradient direction distribution within the neighborhood of key points with statistically invariant scale.

6. The method according to claim 5, characterized in that, The multi-scale matching fusion model includes an input layer that receives directional gradient histogram descriptors from different scale levels, a feature encoding layer containing three residual convolutional blocks, a multi-scale fusion layer that uses a feature pyramid network structure to fuse features at different scales, and a matching head layer containing two fully connected layers that output matching confidence and keypoint coordinate offsets.

7. The method according to claim 6, characterized in that, The multi-scale matching fusion model utilizes an imbalanced data processing mechanism based on adaptive sampling to dynamically adjust the sampling probability of each matching difficulty category according to the sample difficulty coefficient. In each training batch, samples are extracted according to the adjusted sampling probability to construct the training batch.

8. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores program instructions, which, when executed in a computer, are used to perform a real-scene 3D modeling method based on stereo image pair satellite data as described in any one of claims 1-7.

9. A real-scene 3D modeling system based on stereo image pairs and satellite data, characterized in that, The system comprises the computer-readable storage medium of claim 8, wherein the system is a computer, the computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor that executes program instructions stored in the computer-readable storage medium.

Images