A method for satellite DSM matching precision mutual verification based on overlapping area

Abstract

The application provides a kind of satellite DSM matching precision mutual verification method based on overlapping area, belong to satellite DSM technical field, the application is by screening radiation quality qualified forward-looking, orthoview, rear view image and constructs overlapping area index table, utilizes group behavior adaptive streaming perception artificial intelligence model to remove overlapping area feature point gross error in real time with local topological conduction mechanism, combined with streaming block memory cache queue and dynamic perception pyramid down-sampling batch by batch extraction epipolar image disparity field, then disparity field is sent into smooth wavelet deep decomposition manifold artificial intelligence model and is carried out multi-scale decomposition with continuous boundary convergence control equation constraint splicing seam elevation gradient, finally all scene digital surface model splicing fusion and output elevation gradient continuous digital surface model result, solve the technical problem that elevation field is discontinuous at splicing seam when multi-scene satellite digital surface model is spliced in overlapping area, because gross error is removed and elevation gradient constraint is missing.

Description

Technical Field

[0001] This invention belongs to the field of satellite DSM technology, and specifically relates to a method for mutual verification of satellite DSM matching accuracy based on overlapping regions. Background Technology

[0002] In the field of satellite stereo mapping, multi-line array sensors such as ZY-3 extract parallax fields and reconstruct digital surface models from forward-looking, front-looking, and back-looking linear array images, becoming the mainstream technical approach for large-scale topographic mapping. In multi-image stitching, existing methods typically use regional network adjustment based on rational function models to constrain the geometric relationships between adjacent scenes, extract corresponding points in overlapping areas using traditional feature matching algorithms, and then eliminate inter-scene elevation deviations through least-squares adjustment. However, under complex conditions such as cloud cover, radiation anomalies, or severe topographic relief, the concentration of corresponding points in overlapping areas often introduces a large number of gross errors. Traditional global matrix inversion adjustment is extremely sensitive to these gross errors; a single gross error node can contaminate the entire batch of adjustment results through matrix propagation, leading to distortion of the elevation difference field. Simultaneously, existing stitching methods lack explicit constraint mechanisms on the elevation gradients on both sides of the stitching seam. Local elevation discontinuities at the stitching seam cannot be effectively perceived and corrected by the adjustment process, manifesting as periodic jagged edges and fault artifacts in high-resolution digital surface model results. In other words, existing technologies have technical problems such as discontinuity of the elevation field at the splicing seam when splicing multiple satellite digital surface models in overlapping areas due to the lack of gross error elimination and elevation gradient constraints. Summary of the Invention

[0003] In view of this, the present invention provides a method for mutual verification of satellite DSM matching accuracy based on overlapping regions, which can solve the technical problem in the prior art that the elevation field at the splicing seam is discontinuous due to gross error elimination and lack of elevation gradient constraints when splicing multiple satellite digital surface models in overlapping regions.

[0004] This invention is implemented as follows: This invention provides a method for mutual verification of satellite DSM matching accuracy based on overlapping regions, comprising the following steps:

[0005] Acquire multi-sensor calibrated images, select forward-looking, front-looking, and back-looking images with qualified radiation quality levels according to the radiation quality grading system, extract the vector boundary of the overlapping area range, and construct an overlapping area index table;

[0006] The overlapping area images are fed into the group behavior adaptive streaming perception model to extract the spatial coordinate features and grayscale features of each image in the overlapping area, remove gross errors in real time, and output the feature point set of the overlapping area after removing gross errors.

[0007] For the feature point set of overlapping areas after removing gross errors, image data is loaded into a streaming block memory cache queue. Combined with dynamic perception pyramid downsampling processing, the disparity fields of the two sets of epipolar images, namely the forward-looking and forward-looking back-looking images, are extracted in batches.

[0008] Two sets of disparity fields are fed into a stationary wavelet depth decomposition manifold model. High-frequency disparity texture and low-frequency elevation trend are extracted from multiple wavelet subbands. The features of each subband are fused by a cross-scale topological transmission matrix to output the digital surface model results after elevation gradient alignment in the splicing seam region.

[0009] Calculate the elevation difference field of the digital surface model of the overlapping area of ​​adjacent scenes, calculate the mutual verification accuracy index of the overlapping area based on the dynamic node transmission weight evaluation function, and trigger the corresponding adjustment compensation strategy when the mutual verification accuracy index falls into different intervals to update the rational function model parameters.

[0010] All scene digital surface models after updating rational function model parameters are stitched and fused together. The elevation gradient at the stitching seam is constrained by the continuous boundary convergence control equation to smoothly transition, and the final mutually verified digital surface model result is output.

[0011] DSM refers to Digital Surface Model.

[0012] The specific rules for determining the radiation quality grading system are as follows: data with a cloud coverage area of ​​no more than 5% and no radiation anomalies are rated as Level 1; data with a cloud coverage area of ​​more than 5% but no more than 20% and no radiation anomalies, or cloud cover of no more than 20% and slight radiation anomalies, are rated as Level 2; data with cloud cover of more than 20% but no more than 60% and radiation anomalies not exceeding slight levels, or cloud cover of no more than 60% and moderate radiation anomalies, are rated as Level 3; data with cloud cover of more than 60% or severe radiation anomalies are rated as Level 4; Level 1 data is preferred, and if Level 1 data is insufficient, it is supplemented by downgrading to the next level.

[0013] The overlapping area index table records the spatial overlap range of adjacent scene images and the corresponding scene number using vector polygons, which is used for subsequent batch scheduling of overlapping area image data.

[0014] The structure of the adaptive streaming perception model for group behavior is as follows: the front layer is a spatially orthogonal decoupled channel mapping layer, which maps the input multi-temporal spatial coordinates and grayscale features to mutually orthogonal feature subspaces; after the mapping layer is a dynamic sparse perception layer, which embeds a learnable mask matrix to perform propagation pruning on neuron nodes whose spatial uncertainty exceeds the spatial uncertainty threshold; the core inference layer adopts a multi-level self-circulating feedback architecture, which models the hidden layer neuron nodes as mimic particles, and the particles simulate group collision avoidance behavior through mutual repulsion force operators; the output layer outputs the gross error probability vector of each feature point, and points whose probability exceeds the gross error probability threshold are marked as gross errors and removed.

[0015] In the adaptive streaming perception model for group behavior, when the predicted elevation value of a certain track deviates from the average value of the current batch, the collision penalty term reduces the transmission weight of the corresponding node according to the quadratic function relationship of the deviation magnitude; the transmission weight of neurons between layers is adaptively allocated according to the memory usage ratio of the feature maps of each layer; the feature data between each layer flows in memory in a circular queue, and the queue depth is dynamically adjusted according to the amount of data in the current batch.

[0016] Specifically, the streaming block memory cache queue divides the image data into several sub-blocks according to spatial blocks, loads them into memory sequentially in a first-in-first-out queue, retains only the sub-block data required for the current processing in each batch, and releases it immediately after processing, keeping the peak memory usage within the memory usage ratio threshold.

[0017] Specifically, dynamic perception pyramid downsampling involves dynamically selecting the downsampling ratio based on the texture information entropy of the current sub-block during the image pyramid construction process. Sub-blocks with texture information entropy lower than the information entropy threshold adopt a larger downsampling ratio, while sub-blocks with texture information entropy higher than the information entropy threshold maintain a smaller downsampling ratio.

[0018] The structure of the stationary wavelet depth decomposition manifold model is as follows: the input layer receives two sets of epipolar image disparity fields, one for forward-looking and one for back-looking; the first level is a multi-layer stationary wavelet transform decomposition layer, which cancels the downsampling operation and keeps the feature maps at each scale consistent with the input resolution; each wavelet decomposition generates one low-frequency sub-band and three high-frequency sub-bands in three directions; each sub-band establishes an independent manifold feature propagation channel; the deepest layer introduces a continuous boundary convergence control equation; and the output layer reconstructs the digital surface model result by inverse transforming the features after aligning each sub-band.

[0019] Specifically, the cross-scale topological transmission matrix involves establishing a learnable linear transformation operator between feature maps at different wavelet decomposition levels. Through matrix multiplication, the local disparity gradient information of the low-level high-frequency sub-band is transmitted to the high-level low-frequency sub-band, while the global elevation trend of the high-level low-frequency sub-band is transmitted in reverse to the low-level high-frequency sub-band, thus achieving bidirectional collaborative fusion of multi-scale features.

[0020] The continuous boundary convergence control equation is as follows: the mean square of the elevation gradient difference within the pixel band within the threshold range of the pixel band width on both sides of the splicing seam is used as the regularization loss term to ensure that the elevation output of the model at the splicing seam position satisfies the gradient continuity constraint.

[0021] The design of the dynamic node propagation weight evaluation function is as follows: the mutual verification accuracy index is calculated by weighting three indicators: the mean absolute value of the elevation difference of all pairs of corresponding points in the overlapping area, the standard deviation of the elevation difference, and the effective coverage rate of corresponding points. The weight coefficients of the three indicators are determined by parameter search on multiple sets of experimental data with the goal of minimizing the number of subsequent adjustment compensation iterations.

[0022] Specifically, the adjustment compensation strategy is as follows: when the mutual verification accuracy index is not lower than the second accuracy threshold, the matching quality is deemed qualified and adjustment compensation is not triggered; when the mutual verification accuracy index is not lower than the first accuracy threshold but lower than the second accuracy threshold, a mild compensation strategy is triggered to apply affine transformation correction to the rational function model parameters; when the mutual verification accuracy index is lower than the first accuracy threshold, a global regional network re-adjustment strategy is triggered to resolve the rational function model parameters of all scenes in the current batch.

[0023] The training dataset for the adaptive streaming perception model of group behavior is specifically as follows: multiple overlapping image pairs covering four typical landforms—plains, mountains, cities, and water areas—are selected from historical satellite stereo image archives. For each pair of overlapping images, the locations of coarse and non-coarse points are labeled manually and semi-automatically to form a binary label set. The training subset and validation subset are divided by stratified sampling according to landform type. Data augmentation operations such as radiative noise perturbation, cloud and snow mask simulation, and geometric micro-offset are applied to the training subset.

[0024] The training dataset for the stationary wavelet depth decomposition manifold model is specifically designed as follows: scene pairs containing obvious sensor stitching misalignment are selected from historical stereo images, and normal scene pairs without stitching misalignment are combined to form a training set with balanced positive and negative samples. LiDAR point clouds or high-precision airborne digital elevation models are used as the ground truth elevation references. The ground truth elevation gradient at the stitching seam of each sample is calculated. Data augmentation operations such as random radiation noise, parallax offset, and occlusion simulation are applied to the training set.

[0025] The search range for the spatial uncertainty threshold is 0.15–0.45; the search range for the gross error probability threshold is 0.50–0.80; the memory usage ratio threshold is set to control peak memory usage within 15%–35% of the total data volume; and the search range for the information entropy threshold is 1.2–3.5. The pixel bandwidth threshold ranges from 3 to 9 pixels; the search range for the first precision threshold is 0.40 to 0.60, and the search range for the second precision threshold is 0.70 to 0.85; in the transmission channels of the high-frequency subband and the low-frequency subband, the high-frequency channel uses depthwise separable convolution, and the low-frequency channel uses large receptive field convolution, with the number of basic channels ranging from 32 to 128 and the number of decomposition levels ranging from 3 to 5.

[0026] This invention introduces a group behavior adaptive streaming perception model to perform real-time gross error removal on feature points in overlapping regions. It replaces global matrix inversion with a local topological propagation mechanism, allowing the abnormal propagation of individual gross error nodes to be suppressed immediately by a collision penalty mechanism, thus preventing gross error contamination from spreading to the entire batch of adjustment results. Furthermore, this invention employs a stationary wavelet depth decomposition manifold model to perform multi-scale decomposition of the disparity field and uses the continuous boundary convergence control equation as a strong constraint regularization term. This directly penalizes the elevation gradient discontinuity on both sides of the stitching seam at the loss function level, enabling the model to automatically smooth elevation gradient jumps in the stitching seam region during the inference phase. In summary, this invention solves the technical problem mentioned in the background art of elevation field discontinuity at stitching seams in multi-scene satellite digital surface models due to the lack of gross error removal and elevation gradient constraints during overlapping region stitching. Attached Figure Description

[0027] Figure 1 This is a flowchart of the method of the present invention.

[0028] Figure 2 A comparison chart of the standard deviation of elevation difference before and after removing gross errors in overlapping areas of different landform types.

[0029] Figure 3 A schematic diagram of the subband conduction channel structure for a stationary wavelet depth decomposition manifold model.

[0030] Figure 4 The graph shows the variation of the elevation residual of the splicing seam in the validation set with the training rounds during the training of the stable wavelet depth decomposition manifold model.

[0031] Figure 5 This is a comparison chart of the mean absolute values ​​of the elevation differences in the overlapping areas of OL-05 before and after the global regional network re-adjustment.

[0032] Figure 6 This is a graph showing the trend of elevation gradient jump variables at the splicing seam before and after splicing and fusion.

[0033] Figure 7 The final output is an elevation field distribution map of a multi-scene stitched digital surface model. Detailed Implementation

[0034] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0035] like Figure 1 The diagram shows a flowchart of a method for mutual verification of satellite DSM matching accuracy based on overlapping regions provided by the present invention. This method includes the following steps:

[0036] S01. Acquire multi-view resource No. 3 sensor calibrated images, select forward-looking, front-looking and back-looking images with radiation quality level 1 according to the radiation quality classification system, extract the vector boundary of the overlapping area range, and construct the overlapping area index table.

[0037] S02. Input the overlapping area images into the group behavior adaptive streaming perception model, extract the spatial coordinate features and grayscale features of each image in the overlapping area, remove gross errors in real time, and output the feature point set of the overlapping area after removing gross errors.

[0038] S03. For the feature point set of overlapping areas after removing gross errors, load the ultra-large image data with a streaming block memory cache queue, and combine dynamic perception pyramid downsampling processing to extract the disparity field of the two sets of epipolar images, namely the forward-looking and forward-looking back-looking images, in batches.

[0039] S04. Input the two sets of disparity fields into the stationary wavelet depth decomposition manifold model, extract the high-frequency disparity texture and low-frequency elevation trend on multiple wavelet sub-bands respectively, fuse the features of each sub-band through the cross-scale topological transmission matrix, and output the digital surface model result after the elevation gradient of the splicing seam region is aligned.

[0040] S05. Calculate the elevation difference field of the digital surface model of the overlapping area of ​​adjacent scenes, calculate the mutual verification accuracy index of the overlapping area based on the dynamic node transmission weight evaluation function, and trigger the corresponding adjustment compensation strategy when the mutual verification accuracy index falls into different intervals to update the rational function model parameters.

[0041] S06. After updating the rational function model parameters, all scene digital surface models are spliced ​​and fused. The elevation gradient at the splicing seam is constrained by the continuous boundary convergence control equation to smoothly transition, and the final mutually verified digital surface model results are output.

[0042] The criteria for the radiation quality grading system are as follows: data with cloud coverage of no more than 5% and no radiation anomalies are rated as Level 1; data with cloud coverage of more than 5% but no more than 20% and no radiation anomalies, or cloud cover of no more than 20% and slight radiation anomalies, are rated as Level 2; data with cloud cover of more than 20% but no more than 60% and radiation anomalies not exceeding slight levels, or cloud cover of no more than 60% and moderate radiation anomalies, are rated as Level 3; data with cloud cover of more than 60% or severe radiation anomalies are rated as Level 4. This method prioritizes Level 1 data, and when Level 1 data is insufficient, it is supplemented by downgrading the data accordingly.

[0043] The overlapping area index table records the spatial overlap range of adjacent scene images and the corresponding scene number using vector polygons. This is used to schedule overlapping area image data in batches later, avoiding redundant loading of the entire scene data.

[0044] The adaptive streaming sensing model for group behavior is based on deep learning. The specific structure of the model is as follows: the front layer is a spatially orthogonal decoupling channel mapping layer, which maps the input multi-temporal spatial coordinates and grayscale features to mutually orthogonal feature subspaces, eliminating radiation coupling interference between different sensor channels; after the mapping layer is a dynamic sparse sensing layer, which embeds a learnable mask matrix. The mask weights are adaptively updated during training, addressing spatial uncertainties exceeding a threshold. Perform conduction pruning on the neuron nodes, threshold Iterative experiments were conducted on overlapping regions of multiple sets of real satellite imagery to verify the optimal combination of gross error removal accuracy and computational cost. An initial search range of 0.15–0.45 is recommended. The core inference layer employs a multi-level self-recurrent feedback architecture, modeling hidden layer neurons as mimicry particles. These particles simulate group collision avoidance behavior through mutual repulsion operators. When the predicted elevation value of a certain orbit deviates from the average of the current batch, the collision penalty term reduces the transmission weight of that node according to a quadratic function of the deviation magnitude. The transmission weights between neurons are adaptively allocated based on the memory usage ratio of each layer's feature maps. Layers with higher memory usage receive lower weight redundancy limits, thus compressing the overall peak memory usage. Feature data between layers circulates in memory using a circular queue, with the queue depth dynamically adjusted based on the current batch data volume to avoid memory overflow caused by large batches. The output layer outputs the gross error probability vector for each feature point; if the probability exceeds a threshold... Points are marked as gross errors and removed, threshold It is recommended to determine the search range of 0.50–0.80 by comprehensively balancing the false positive rate and the false negative rate in comparative experiments in different landform regions. The network should not rely on global matrix storage, but only complete the self-organized calibration of heterogeneous elevation data through local topology propagation, thereby avoiding deadlock bottlenecks during large-scale multi-track data parallel processing. The specific steps for establishing the training dataset for the adaptive streaming perception model of group behavior include: selecting multiple overlapping image pairs covering four typical landforms (plains, mountains, cities, and water bodies) from historical ZY-3 satellite stereo image archives; for each pair of overlapping image pairs, manually and semi-automatically labeling the locations of coarse and non-coarse points to form a binary label set; stratifying the label set into training and validation subsets by landform type with a ratio of 8:2; and applying data augmentation operations such as radiative noise perturbation, cloud and snow mask simulation, and geometric micro-offset to the training subset to expand sample diversity. The specific steps for training the adaptive streaming perception model for group behavior include: using the weighted sum of the cross-entropy loss function and the collision penalty term as the total loss; setting the weight of the collision penalty term to a small value in the early stages of training, and linearly increasing it to a preset upper limit with each training round; the upper limit value is recommended to be determined in the range of 0.05 to 0.30 through preliminary experiments; and using an adaptive moment estimation optimizer to update the parameters, with an initial learning rate recommended to be... ~ The selection criteria are between [a certain range], with the accuracy of gross error removal on the validation set as the convergence criterion; during training, a sparse regularization constraint is applied to the learnable mask matrix, and the regularization coefficient is suggested to be [a certain value]. ~ The relationship was determined experimentally; when the improvement in validation set accuracy was less than [a certain value] over five consecutive training epochs, the result was [a certain value]. Training ends when the time is right.

[0045] The adaptive streaming sensing model for group behavior models hidden layer neurons as a swarm of mutually repulsive mimicking particles. When an anomaly occurs in the elevation prediction of a single track, the collision penalty mechanism instantly reduces the propagation weight of that node, thus achieving local self-organized removal of gross errors in multi-track intersection areas without relying on global matrix solving. This local topological propagation mechanism fundamentally avoids the bottleneck of global matrix inversion in large-scale parallel adjustment, enabling the system to maintain stable convergence when processing high-concurrency multi-track data. Furthermore, the dynamic allocation strategy of GPU memory and system memory effectively suppresses the peak resource consumption in ultra-large image batch processing, comprehensively improving the real-time performance and reliability of gross error removal in overlapping areas.

[0046] The dynamic node propagation weight evaluation function is designed as follows: the average of the absolute values ​​of the elevation differences of all pairs of corresponding points within the overlapping area is used. Standard deviation of elevation difference Effective coverage of homonymous points The cross-validation accuracy index is calculated using a weighted average of three indicators. The weighting coefficients of the three indicators were determined through parameter search on multiple sets of experimental data with the objective of minimizing the number of subsequent adjustment and compensation iterations. The suggested ranges for the weighting coefficients are 0.3–0.5, 0.2–0.4, and 0.2–0.4, respectively. When the current overlapping area is deemed to have acceptable matching quality, adjustment compensation is not triggered; when When a mild compensation strategy is triggered, an affine transformation is applied to correct the parameters of the rational function model; when When this occurs, the global regional network re-adjustment strategy is triggered, and the rational function model parameters for all scenes in the current batch are re-solved. and The specific values ​​were determined through gradient experiments conducted in regions with different landform types, using the minimum residual of subsequent digital surface model stitching as the criterion for iterative analysis. The initial search range is 0.40 to 0.60. The initial search range is 0.70 to 0.85.

[0047] Among them, the streaming block memory cache queue refers to dividing the ultra-large image data into several sub-blocks according to spatial blocks, loading them into memory sequentially in a first-in-first-out queue manner, retaining only the sub-block data required for current processing in each batch, and releasing it immediately after processing is completed, thereby controlling the peak memory usage within 15% to 35% of the total data volume. It is recommended to determine the specific ratio through pre-experiments under different image widths and server memory configurations.

[0048] Dynamic perceptual pyramid downsampling refers to dynamically selecting the downsampling ratio based on the texture information entropy of the current sub-block during the image pyramid construction process. Sub-blocks with information entropy below the threshold adopt a larger downsampling ratio to reduce the amount of computation, while sub-blocks with information entropy above the threshold maintain a smaller downsampling ratio to preserve texture details. It is recommended that the information entropy threshold be determined through statistical analysis of images of various landform types, and the recommended search range is 1.2 to 3.5 bits / pixel.

[0049] The stationary wavelet depth decomposition manifold model is based on deep learning. The specific structure of the model is as follows: the input layer receives two sets of epipolar image disparity fields: forward-looking and forward-looking backward-looking. The first level is a multi-layer stationary wavelet transform decomposition layer, which eliminates the downsampling operation of traditional convolutional networks, maintaining the feature maps at each scale consistent with the input resolution, thus preserving translation invariance. Each wavelet decomposition level generates one low-frequency sub-band and three high-frequency sub-bands in three directions, with a total of 3-5 decomposition levels. The number of decomposition levels is recommended to be determined experimentally on test areas with different terrain undulations, using the minimum seam elevation residual as the criterion. Each sub-band establishes an independent manifold feature transmission channel. The high-frequency channel uses depth-separable convolution to extract microscopic disparity textures, while the low-frequency channel... The model employs large receptive field convolutions to capture macroscopic elevation trends; cross-scale topological transmission matrices achieve feature fusion between different sub-bands through learnable linear transformations, with matrix parameters updated end-to-end during training; the deepest layer introduces a continuous boundary convergence control equation, using the squared mean of the difference in elevation gradients on both sides of the seam as a regularization term in the loss function, forcing each sub-band to smoothly align its elevation gradient in the seam region, thereby eliminating periodic jaggedness and faults; the output layer reconstructs the final digital surface model result by inverse transforming the aligned features of each sub-band; the number of channels in each convolutional layer of the model increases exponentially with layer depth, and the specific number of channels is recommended to be determined through ablation experiments under computational resource constraints, with a suggested range of 32–128 channels. The steps for establishing the training dataset for the stationary wavelet depth decomposition manifold model specifically include: selecting scene pairs containing obvious sensor stitching misalignment from historical ZY-3 stereo imagery, and combining them with normal scene pairs without stitching misalignment to form a training set with balanced positive and negative samples; using LiDAR point clouds or high-precision airborne digital elevation models as the ground truth elevation reference, calculating the ground truth elevation gradient at the stitching seam of each sample; and applying data augmentation operations such as random radiometric noise, disparity offset, and occlusion simulation to the training set. The steps for training the stationary wavelet depth decomposition manifold model specifically include: the total loss is the weighted sum of the mean square error of elevation reconstruction and the continuous boundary convergence regularization term; the regularization weight is recommended to be determined through pre-experiments within the range of 0.01 to 0.10; a cosine annealing learning rate scheduling strategy is adopted, and the initial learning rate is recommended to be... ~ The selection is between; when the elevation residual of the splice seam in the validation set decreases by less than 10 consecutive training rounds. Training ends when the time is right.

[0050] The stationary wavelet depth decomposition manifold model maintains translation invariance by canceling downsampling, enabling the microscopic parallax texture and macroscopic elevation trend of the stitching seam region to be modeled synchronously in the same resolution space. This avoids the problem of blurred stitching seam details caused by resolution loss in traditional downsampling convolutional networks. The continuous boundary convergence control equation, as a strong constraint regularization term, directly penalizes the elevation discontinuity on both sides of the stitching seam at the loss function level. This allows the model to automatically smooth the elevation gradient jump in the stitching seam region during the inference stage after training convergence. Thus, it achieves the restoration of the continuity of the elevation field at the stitching seam without relying on manual editing, significantly improving the overall consistency of the digital surface model of multi-sensor stitched images.

[0051] Among them, the cross-scale topological transmission matrix refers to the learnable linear transformation operator established between feature maps of different wavelet decomposition levels. Through matrix multiplication, the local disparity gradient information of the low-level high-frequency sub-band is transmitted to the high-level low-frequency sub-band, while the global elevation trend of the high-level low-frequency sub-band is transmitted in reverse to the low-level high-frequency sub-band, realizing the bidirectional collaborative fusion of multi-scale features.

[0052] Among them, the continuous boundary convergence control equation refers to using the square mean of the elevation gradient difference within a pixel band of a certain width on both sides of the splicing seam as the regularization loss term. It is recommended that the width be between 3 and 9 pixels and determined experimentally. The regularization term ensures that the elevation output of the model at the splicing seam position satisfies the gradient continuity constraint.

[0053] Among them, rational function model parameters refer to the rational polynomial coefficients that describe the mapping relationship between satellite image coordinates and ground object coordinates, and affine transformation correction refers to superimposing affine transformation compensation on the existing rational function model parameters to correct systematic geometric deviations.

[0054] As an optional implementation, the present invention also provides a system for implementing the above method on a computer. The system includes a processor and a readable storage medium. The readable storage medium stores program instructions. When the program instructions are run on the processor, they sequentially call a streaming block memory cache queue module, a group behavior adaptive streaming perception model inference module, and a stationary wavelet depth decomposition manifold model inference module to complete the entire processing flow from sensor-corrected image acquisition to the output of a mutually verified qualified digital surface model.

[0055] The specific implementation of step S01 is as follows: The purpose of this step is to screen images with acceptable radiation quality and construct an overlapping area index to establish a data foundation for subsequent processing. First, multiple scenes of Ziyuan-3 sensor-corrected images are acquired in batches from the satellite data archive. For each image, the cloud coverage area ratio and radiation anomaly level are calculated, and the images are graded according to the radiation quality grading system: images with a cloud coverage area of ​​no more than 5% and no radiation anomaly are graded as Level 1; images with a cloud coverage area greater than 5% but no more than 20% and no radiation anomaly, or cloud cover no more than 20% and slight radiation anomaly, are graded as Level 2, and so on up to Level 4. Level 1 images are retained first, and Level 2 and Level 3 data are added sequentially if the number is insufficient. For the screened images, spatial intersection operations are performed using the geographic coverage polygons in the image metadata to extract the overlapping area range between adjacent scenes and generate vector boundaries. The overlapping area index table records the spatial overlap range of each pair of adjacent scenes and the corresponding scene number using vector polygons, and is stored in descending order of overlap area. This is used for subsequent batch scheduling of overlapping area image data as needed, avoiding redundant loading of entire scene data.

[0056] The specific implementation of step S02 is as follows: The purpose of this step is to remove gross feature points from the overlapping area image in real time, providing reliable data for subsequent disparity extraction. The overlapping area image is input into the group behavior adaptive streaming perception model. The model's front layer is a spatial orthogonal decoupling channel mapping layer, which orthogonally maps multi-temporal spatial coordinates and grayscale features to eliminate radiation coupling interference between different sensor channels. The dynamic sparse perception layer embeds a learnable mask matrix to perform propagation pruning on neuron nodes whose spatial uncertainty exceeds the spatial uncertainty threshold (suggested search range 0.15~0.45). The core inference layer models the hidden layer neurons as mimic particles. The particles simulate group collision avoidance behavior through mutual repulsion force operators. When the elevation prediction value of a node deviates from the batch mean, the collision penalty term reduces the propagation weight of that node according to the quadratic function relationship of the deviation magnitude, thereby completing the self-organization suppression of gross nodes within the local topology range without the need for global matrix inversion. The feature data of each layer flows in memory in a circular queue, and the queue depth is dynamically adjusted to prevent memory overflow. The output layer outputs the gross error probability vector of each feature point. Points with probabilities exceeding the gross error probability threshold (suggested search range 0.50 to 0.80) are removed. Finally, the output layer is the feature point set of the overlapping region after removing gross errors.

[0057] The specific implementation of step S03 is as follows: The purpose of this step is to efficiently extract the disparity field of epipolar images under memory-constrained conditions. For the feature point set of overlapping regions after removing gross errors, a streaming block-based memory cache queue is used to divide the ultra-large image data into several sub-blocks according to spatial blocks. These sub-blocks are loaded into memory sequentially using a first-in-first-out queue. Each batch retains only the sub-block data required for current processing and releases it immediately after processing, keeping peak memory usage within 15% to 35% of the total data volume. During the pyramid construction stage, dynamic perceptual pyramid downsampling dynamically selects the downsampling ratio based on the texture entropy of the current sub-block: when the texture entropy is lower than the entropy threshold (recommended search range 1.2 to 3.5). Sub-blocks with a larger downsampling ratio are used to reduce computation, while sub-blocks with texture entropy higher than a threshold are kept at a smaller downsampling ratio to preserve texture details. Based on the above processing, epipolar stereo matching is performed on the forward-looking frontal image pairs and the frontal-looking backal image pairs in batches. A semi-global matching algorithm is used to calculate the cost aggregation disparity field, and two sets of epipolar image disparity fields are output as input for the next step of elevation reconstruction.

[0058] The specific implementation of step S04 is as follows: The purpose of this step is to reconstruct the elevation field from the disparity field and eliminate elevation gradient jumps at the seams. Two sets of disparity fields are input into a stationary wavelet depth decomposition manifold model. The first level of this model is a multi-layer stationary wavelet transform decomposition layer, eliminating the downsampling operation of traditional convolutional networks and maintaining the consistency of feature maps at each scale with the input resolution, thereby preserving translation invariance and avoiding blurring of seam details due to resolution loss. Each wavelet decomposition level generates one low-frequency sub-band and three high-frequency sub-bands in three directions, with a total of 3 to 5 decomposition levels. The number of decomposition levels is experimentally determined on test areas with different terrain undulations, using the minimum elevation residual at the seams as the criterion. The high-frequency sub-band uses depth-separable convolution to extract microscopic disparity textures, while the low-frequency sub-band uses large receptive field convolution to capture macroscopic elevation trends. A cross-scale topological transmission matrix enables bidirectional feature fusion between different sub-bands through a learnable linear transformation. It transmits local disparity gradient information from low-level high-frequency sub-bands to high-level low-frequency sub-bands, while simultaneously transmitting the global elevation trend from high-level low-frequency sub-bands back to low-level high-frequency sub-bands. The deepest layer introduces a continuous boundary convergence control equation, using the squared mean of the elevation gradient differences within 3–9 pixels on each side of the seam as the regularization loss term, forcing smooth alignment of elevation gradients in the seam region among the sub-bands. The output layer reconstructs the digital surface model result from the inverse transformation of the aligned features of each sub-band.

[0059] The specific implementation of step S05 is as follows: The purpose of this step is to quantify the matching accuracy of overlapping areas and trigger the corresponding adjustment compensation strategy. The elevation difference field of the digital surface model of the overlapping areas of adjacent scenes is calculated, and the mean absolute value of the elevation difference of all pairs of corresponding points is extracted. Standard deviation of elevation difference Effective coverage of homonymous points The mutual verification accuracy index is calculated using three weighted indicators based on the dynamic node propagation weight evaluation function. The suggested weighting coefficients for the three indicators are 0.3–0.5, 0.2–0.4, and 0.2–0.4, respectively, determined by searching multiple sets of experimental data with the objective of minimizing the number of subsequent adjustment and compensation iterations. When the accuracy is not lower than the second precision threshold (recommended search range 0.70–0.85), the matching quality is considered acceptable, and adjustment compensation is not triggered; when… When the accuracy is not lower than the first precision threshold (recommended search range 0.40–0.60) and lower than the second precision threshold, a mild compensation strategy is triggered to correct systematic geometric deviations by superimposing affine transformation compensation on the rational function model parameters; when When the accuracy falls below the first precision threshold, a global regional network re-adjustment strategy is triggered, resolving the rational function model parameters for all scenes in the current batch. The updated rational function model parameters are used for subsequent stitching and fusion.

[0060] The specific implementation of step S06 is as follows: The purpose of this step is to stitch together all scene digital surface models into a final result with a continuous elevation field. For all scene digital surface models after updating the rational function model parameters, stitching and fusion processing is performed batch by batch according to the scene number correspondence in the overlapping region index table. At the stitching seam, a continuous boundary convergence control equation is used to constrain the smooth transition of the elevation gradient. This equation uses the square mean of the elevation gradient difference within the pixel bands on both sides of the stitching seam as a penalty term, driving the elevation field on both sides of the stitching seam to converge to a continuous state at the gradient level, thereby eliminating periodic jagged edges and toposcopic artifacts. Finally, a mutually verified and qualified digital surface model result is output, and the continuity of the elevation gradient in the result is guaranteed globally by the continuous boundary convergence control equation.

[0061] It should be noted that the key technologies of this invention include: First, a local topological gross error removal mechanism for the adaptive streaming perception model of group behavior. By modeling neurons as mimicry particles and introducing a collision penalty term, the transmission weight of a single abnormal node is suppressed in real time, achieving self-organized removal of gross errors without the need for global matrix storage, fundamentally avoiding the problem of a single gross error node contaminating the entire batch of results in global matrix adjustment; Second, a non-subsampling multi-scale modeling mechanism for the stationary wavelet depth decomposition manifold model. By maintaining the feature maps at each scale consistent with the input resolution, microscopic parallax textures and macroscopic elevation trends are modeled collaboratively in the same resolution space. Combined with the explicit penalty of the seam gradient by the continuous boundary convergence control equation, the model can automatically smooth the elevation gradient jump in the seam region during the inference stage after training convergence; Third, a hierarchical adjustment compensation mechanism driven by a dynamic node transmission weight evaluation function. By quantifying the matching quality through three weighted indicators and mapping it to two strategies, namely mild affine correction or global re-adjustment, adaptive matching between adjustment cost and accuracy requirements is achieved. The three mechanisms work together to form a closed loop of gross error removal, elevation reconstruction and adjustment compensation, eliminating the root cause of elevation discontinuity at splice seams from three dimensions: data quality, model constraints and parameter updates.

[0062] It should be noted that in large-scale multi-track satellite stereo mapping scenarios, when multiple landform types such as plains, mountains, cities, and water bodies exist simultaneously in the same survey area, the texture information entropy of different landform regions varies greatly. A uniform downsampling ratio cannot simultaneously meet the computational efficiency requirements of low-texture areas and the detail preservation requirements of high-texture areas. This problem is particularly prominent in batch processing of ultra-wide-format images. The reason for the above technical problem is that traditional pyramid construction uses a fixed downsampling ratio. Low-texture areas waste computational resources due to excessive detail preservation, while high-texture areas lose the local features required for disparity matching due to oversampling. This causes accuracy jumps in the disparity field at landform boundaries, which in turn leads to local distortions in the elevation reconstruction results. The usual solution to the above technical problem is to adopt a uniform multi-resolution processing strategy for panoramic images or to manually delineate landform zones and set downsampling parameters for each zone. However, manually delineating zones is difficult to implement in automated processing workflows, and a uniform strategy cannot adapt to drastic changes in landforms within the same scene. Especially when mountains and water bodies coexist in the same sub-block, it is difficult to simultaneously meet the processing needs of both types of areas, regardless of whether a large or small downsampling ratio is used. This invention effectively solves this technical problem. Dynamic perceptual pyramid downsampling independently calculates the texture information entropy for each sub-block during pyramid construction and dynamically selects the downsampling ratio for that sub-block based on the comparison between the information entropy and a threshold. This allows low-texture sub-blocks to be downsampled at a large ratio to compress computation, while high-texture sub-blocks are downsampled at a small ratio to preserve details. This automatically adapts to the diverse terrain within the same scene during a single pyramid construction process. This mechanism, combined with a streaming block-based memory cache queue, ensures that the downsampling decision for each sub-block is completed during the loading phase, eliminating the need for pre-scanning of panoramic data. This achieves an adaptive balance between computational efficiency and disparity accuracy at the sub-block level, thereby eliminating the disparity accuracy jump problem in multi-terrain scenarios where a fixed downsampling strategy is used.

[0063] Specifically, the principle of this invention is as follows: The fundamental reason this invention can solve the aforementioned technical problems lies in the synergistic design of two core mechanisms. First, the localization of the gross error removal mechanism. Traditional global matrix adjustment incorporates all points of the same name into a unified linear equation system. Gross errors are distributed to other nodes through matrix coupling, leading to an overall shift in the adjustment results. This invention models hidden layer neurons as mutually repulsive pseudo-particles. When the elevation prediction value of a node deviates from the batch mean, the collision penalty term reduces the transmission weight of that node according to a quadratic function of the deviation magnitude, ensuring that the abnormal influence is absorbed locally rather than propagated globally. This mechanism does not rely on global matrix storage, fundamentally cutting off the cross-node contamination path of gross errors and ensuring the reliability of the adjustment input data. Second, the explicit constraint mechanism for the elevation gradient of the stitching seam. The stationary wavelet transform eliminates the downsampling operation of traditional convolutional networks, ensuring that the feature maps at each scale maintain the same input resolution. This allows microscopic parallax textures and macroscopic elevation trends to be modeled synchronously in the same resolution space, avoiding blurring of stitching seam details due to resolution loss. The continuous boundary convergence control equation uses the squared mean of the elevation gradient differences within the pixel bands on both sides of the seam as the regularization loss term. During training, it continuously penalizes gradient discontinuities, causing the model parameters to converge to a state that automatically smooths the seam. During inference, this constraint is implicitly embedded in the model weights, allowing gradient smoothing to be applied to any seam region without manual intervention. These two mechanisms work together: the former ensures the quality of the disparity field data entering the elevation reconstruction model, while the latter ensures the geometric continuity of the reconstruction results at the seam. This collaborative approach eliminates the root cause of elevation discontinuities at the seam from both the input data quality and output constraints perspectives.

[0064] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0065] The specific implementation of step S01 is as follows: After acquiring multi-scene calibrated images from the Resource 3 sensor, each image is graded according to the radiation quality grading system. Images with cloud coverage of no more than 5% and no radiation anomalies are rated as Level 1 and prioritized for inclusion in the processing flow; if Level 1 data is insufficient, it is downgraded and supplemented accordingly. Subsequently, the spatial overlap range of adjacent images is extracted, and the correspondence between each image number and the overlap boundary is recorded using vector polygons to construct an overlap area index table for subsequent batch scheduling of image data, avoiding redundant loading of entire scene data. This step mainly uses conventional geographic information processing operations and does not require additional calculation formulas.

[0066] The specific implementation of step S02 is as follows: The overlapping area image is fed into the group behavior adaptive streaming sensing model. This model is based on deep learning and has a four-layer structure. The front layer is a spatially orthogonal decoupling channel mapping layer, which maps the input multi-temporal spatial coordinates and grayscale features to mutually orthogonal feature subspaces to eliminate radiation coupling interference between different sensor channels. The mapping operation formula is expressed as follows:

[0067] ;

[0068] In the formula, The feature matrix after orthogonal mapping has dimension . The orthogonal mapping weight matrix has dimensions of . It is learned end-to-end during the training process and satisfies orthogonal constraints. , The constraint is the identity matrix, and it is added to the training loss as a regularization term. The input feature matrix has a dimension of . It is composed of the spatial coordinates and grayscale features of each scene in the overlapping area. Dimension of orthogonal subspace Input channel number This represents the number of feature points.

[0069] The mapping layer is followed by a dynamic sparse sensing layer, which embeds a learnable mask matrix to handle spatial uncertainties exceeding a threshold. The neuron nodes are subjected to conduction pruning, and the pruning judgment formula is expressed as follows:

[0070] ;

[0071] In the formula, For the first Activation values ​​of each neuron node This is an estimate of the spatial uncertainty of the node, obtained by approximating the node activation variance during training. The conduction value after pruning For the first The learnable mask values ​​corresponding to each node are adaptively updated during training. To determine the spatial uncertainty pruning threshold, an initial search range of 0.15–0.45 is recommended. This threshold is determined by finding the optimal balance between validation set gross error removal accuracy and computational cost across multiple sets of real satellite imagery data in overlapping areas.

[0072] The core inference layer models the hidden layer neurons as mimicry particles, with these particles simulating group collision avoidance behavior through mutual repulsion operators. When the... Predicted elevation of the track Deviation from the current batch group average When a collision occurs, the collision penalty mechanism reduces the propagation weight of the node according to a quadratic function of the deviation magnitude, as expressed in the following formula:

[0073] ;

[0074] In the formula, After the collision penalty Node propagation weight The basic propagation weights for this node are obtained by initializing the training parameters. The collision penalty coefficient is recommended to be determined through preliminary experiments within the range of 0.05 to 0.30. For the first Elevation prediction values ​​of track nodes The average value of the current batch is calculated using the following formula: This represents the total number of nodes in the current batch. The inter-layer neuron propagation weights are adaptively allocated according to the memory usage ratio of each layer's feature map, as expressed by the following formula:

[0075] ;

[0076] In the formula, For the first Upper limit of propagation weights after layer adaptive allocation For the first Upper limit of layer weight redundancy For the first The ratio of memory usage of the layer feature map to the total memory usage satisfies The adjustment coefficient is typically set to 1.0. This is the layer number index. The total loss function is a weighted sum of the cross-entropy loss and the collision penalty term, expressed as follows:

[0077] ;

[0078] In the formula, Total loss Cross-entropy loss The collision penalty term weight is initially set to a small value during training, and then linearly increases to a preset upper limit with each training epoch. It is recommended that the upper limit be determined within the range of 0.05 to 0.30. This is the batch variance term for the current batch of node elevation predictions, used to quantify the degree of group deviation. The output layer outputs the gross error probability value for each feature point. The probability exceeds the threshold Points were marked as gross errors and discarded. The recommended search range for the gross error threshold is 0.50–0.80, determined through a comprehensive trade-off between false positive and false negative rates in comparative experiments across different landform types. Sparse regularization constraints are applied to the learnable mask matrix. Regularization coefficients are recommended to be in the range of 10 ... ~ The parameters were determined experimentally; an adaptive moment estimator optimizer was used for parameter updates, and the initial learning rate was suggested to be... ~ Select between; when the improvement in validation set accuracy is less than 5 consecutive training epochs. Training terminates at this point. The training dataset is selected from historical ZY-3 satellite stereo imagery, featuring multiple overlapping image pairs covering four typical landforms: plains, mountains, cities, and water areas. Gross and non-gross points are labeled using a combination of manual and semi-automatic methods. Stratified sampling is performed according to landform type, with a training subset to validation subset ratio of 8:2. Data augmentation operations such as radiative noise perturbation, cloud and snow mask simulation, and geometric micro-offset are applied.

[0079] The specific implementation of step S03 is as follows: For the feature point set of overlapping areas after removing gross errors, a streaming block-based memory cache queue is used to load the massive image data. The image is divided into several sub-blocks according to spatial blocks, and these sub-blocks are sequentially loaded into memory using a first-in-first-out queue. Each batch retains only the sub-blocks needed for current processing, and releases them immediately after processing. Peak memory usage is controlled within 15% to 35% of the total data volume. During the image pyramid construction process, dynamic perception pyramid downsampling is performed based on the texture information entropy of the current sub-block. Dynamically select the downsampling ratio. Below the threshold Sub-blocks with larger downsampling ratios are used to reduce computation, while sub-blocks with values ​​above a threshold are kept at smaller downsampling ratios to preserve texture details. The recommended search range for the texture information entropy threshold is 1.2–3.5 bits / pixel, determined through statistical analysis of images of various landform types. The disparity fields of the foreground and background epipolar images are extracted in batches, and two sets of disparity field data are output.

[0080] The specific implementation of step S04 is as follows: Two sets of disparity fields are fed into a stationary wavelet depth decomposition manifold model. After receiving the two sets of epipolar image disparity fields, the input layer performs multi-level stationary wavelet transform decomposition, canceling downsampling to maintain the resolution of feature maps at each scale consistent with the input. A total of 3-5 decomposition levels are set; the specific number of levels is recommended to be determined experimentally using the minimum seam elevation residual criterion. Each decomposition level generates a low-frequency sub-band. With three high-frequency sub-bands , , ,in Numbering is used for decomposition levels. High-frequency channels employ depthwise separable convolution to extract microscopic parallax textures, while low-frequency channels use large receptive field convolution to capture macroscopic elevation trends. Cross-scale topological transmission matrix. Feature fusion between different subbands is achieved using learnable linear transformations, as expressed in the following formula:

[0081] ;

[0082] In the formula, For the first Feature vector after level fusion For the first The multi-scale topological transmission matrix, whose dimensions are determined by the number of channels in each sub-band, is updated end-to-end during training. It transmits the local disparity gradient information of the low-level high-frequency sub-bands to the high-level low-frequency sub-bands through matrix multiplication, while simultaneously transmitting the global elevation trend of the high-level low-frequency sub-bands back to the low-level high-frequency sub-bands. For the first Low-frequency subband feature vector , , The first The high-frequency sub-band feature vectors are generated in three directions. The continuous boundary convergence control equation uses the squared mean of the elevation gradient difference within the pixel bands on both sides of the splicing seam as the regularization loss term, and the formula is as follows:

[0083] ;

[0084] In the formula, For splice seam regularization loss This refers to the total number of pixels included in the calculation within the pixel band on both sides of the seam. The recommended pixel band width is between 3 and 9 pixels, determined experimentally. and These are the left and right sides of the splice seam, respectively. The elevation gradient values ​​at each pixel are both of the same dimension, representing the rate of change of elevation with respect to pixel position. The total loss is a weighted sum of the mean square error of elevation reconstruction and the continuous boundary convergence regularization term, expressed in the following formula:

[0085] ;

[0086] In the formula, The total loss of the stationary wavelet depth decomposition manifold model Mean square error for elevation reconstruction For regularization weights, it is recommended to determine them through preliminary experiments within the range of 0.01 to 0.10. and Both are dimensionless normalized loss values, and their dimensions are consistent. A cosine annealing learning rate scheduling strategy is adopted, and the initial learning rate is recommended to be... ~ The selection is between; when the elevation residual of the splice seam in the validation set decreases by less than 10 consecutive training rounds. Training terminates at this point. The output layer reconstructs the final digital surface model result by inverse transforming the features after aligning each sub-band. The number of channels in each convolutional layer increases exponentially with the layer depth, and the recommended range for the basic number of channels is 32–128. The specific number of channels is determined through ablation experiments. The training dataset is selected from historical ZY-3 stereo imagery, containing balanced samples with obvious stitching seam misalignment and normal scene pairs. LiDAR point clouds or high-precision airborne digital elevation models are used as ground truth elevation references, and data augmentation operations such as random radiation noise, parallax shift, and occlusion simulation are applied.

[0087] The specific implementation of step S05 is as follows: After calculating the elevation difference field of the digital surface model of the overlapping area of ​​adjacent scenes, the mutual verification accuracy index is calculated based on the dynamic node transmission weight evaluation function. The formula is expressed as follows:

[0088] ;

[0089] In the formula, For mutual verification accuracy index The mean of the absolute values ​​of the elevation differences between all corresponding points within the overlapping area. Standard deviation of elevation difference The effective coverage rate of corresponding points ranges from 0 to 1. , , These are the weighting coefficients for the three indicators, with suggested ranges of 0.3–0.5, 0.2–0.4, and 0.2–0.4, respectively. These coefficients are determined through a parameter search aimed at minimizing the number of subsequent adjustment and compensation iterations. Because… and The dimension of is elevation units. It is a dimensionless ratio. and The reciprocal dimension with elevation units. For dimensionless coefficients, make The entire system is normalized to a dimensionless precision index. When When the current overlapping area is deemed to have acceptable matching quality, adjustment compensation is not triggered; when At this time, a mild compensation strategy is triggered, applying an affine transformation to correct the parameters of the rational function model, as expressed in the following formula:

[0090] ;

[0091] In the formula, The corrected rational polynomial coefficient vector The coefficient vector of the rational polynomial before correction The affine transformation compensation vector is obtained by least-squares fitting of the residuals of corresponding point pairs in the overlapping region. When When this occurs, the global regional network re-adjustment strategy is triggered, and the rational function model parameters for all scenes in the current batch are re-solved. and The threshold for the accuracy index classification. The initial search range is recommended to be 0.40 to 0.60. The initial search range is recommended to be 0.70 to 0.85. The specific value should be determined by conducting gradient experiments in different landform regions and using the minimum residual of subsequent digital surface model splicing as the criterion for iterative analysis.

[0092] The specific implementation of step S06 is as follows: All scene digital surface models after updating the rational function model parameters are stitched and fused together. The elevation gradient at the stitching seams is smoothly transitioned using the continuous boundary convergence control equation, and the final mutually verified digital surface model result is output. The relevant constraint equations have been detailed in step S04 and will not be repeated here.

[0093] To better understand and implement this invention, the following is a specific application scenario of this invention, Example 2:

[0094] This embodiment selects a multi-view satellite stereo image dataset of a transitional zone from a mountainous area to an urban area as the processing object. This dataset contains nine ZY-3 sensor-corrected images, covering an area with significant topographic relief and landform types including mountains, urban areas, and water bodies. There is a 10%–30% overlap between adjacent images, and the ground resolution of each image is 2.5m. The complete processing procedure is described below according to steps S01 to S06 of this invention.

[0095] First, step S01 is performed to assess the radiation quality level of each of the nine images according to the radiation quality grading system. The statistical results are shown in Table 1.

[0096] Table 1. Results of Radiometric Quality Assessment for Each Scene

[0097]

[0098] As shown in Table 1, five images (SC-01, SC-02, SC-04, SC-06, and SC-08) were rated as Level 1, SC-03 and SC-05 as Level 2, SC-07 as Level 3, and SC-09 was removed because its cloud cover exceeded 60% and it was rated as Level 4. Subsequent processing will be based on the remaining eight images. Spatial intersection operations were performed on the eight images to extract 11 pairs of overlapping areas between adjacent images. An overlapping area index table was constructed and sorted from largest to smallest for later use.

[0099] In step S02, 11 pairs of overlapping region images are sequentially fed into the adaptive streaming perception model for group behavior to remove gross errors. In this embodiment, the structure of the adaptive streaming perception model for group behavior is as follows: The input layer receives spatial coordinate features and grayscale features processed by the spatial orthogonal decoupling channel mapping layer, with a channel dimension of 8; the dynamic sparse perception layer embeds a learnable mask matrix, with a spatial uncertainty threshold set to 0.30. The mask weights of neurons exceeding the threshold approach zero after training convergence, and the corresponding node propagation is pruned during the inference phase; the core inference layer has 4 layers, with the number of neurons in each layer being 256, 128, 64, and 32 respectively. The neurons in each layer are modeled as mimic particles, and the collision penalty term weight is set to 0.05 in the early stage of training, increasing linearly to 0.20 with each training round; the output layer is a linear mapping from 32-dimensional to 2-dimensional, outputting a gross error probability vector, with a gross error probability threshold set to 0.65. The data flow is as follows: overlapping image fragments are mapped to an 8-dimensional orthogonal feature subspace through a spatial orthogonal decoupling channel mapping layer, then sent to a dynamic sparse perception layer for masking and weighting, and then transmitted through a mimicry particle collision penalty through a 4-layer core inference layer. Finally, the output layer outputs the gross error probability of each feature point, and points with a probability exceeding 0.65 are removed.

[0100] like Figure 2 As shown in the figure, the statistical results of the outlier removal rate of the adaptive streaming perception model for group behavior in this embodiment on 11 pairs of overlapping regions show that, due to more obstruction, the average outlier removal rate in mountainous areas was 18.3%, in urban areas it was 9.7%, and in water areas it was 24.1%. After model processing, the standard deviation of the elevation difference of the corresponding point sets in various landform regions decreased significantly after removal, and the trend of change is shown in Table 2:

[0101] Table 2. Comparison of Standard Deviation of Elevation Difference in Overlapping Areas Before and After Gross Error Removal

[0102]

[0103] As shown in Table 2, as gross errors are gradually eliminated, the standard deviation curve of elevation difference tends to flatten. This trend indicates that the local topological transmission mechanism effectively blocks the bias effect of abnormal nodes on the batch mean, making the retained set of corresponding points closer to the true elevation distribution in a statistical sense.

[0104] In step S03, image data is loaded using a streaming, block-based memory cache queue. The sub-block size is set to 1024×1024 pixels, and the queue depth is dynamically adjusted. Peak memory usage is controlled to within 22% of the total data volume. The information entropy threshold for dynamic sensing pyramid downsampling is set to 2.1. The average information entropy of the mountainous sub-blocks is 3.4. Above the threshold, a small downsampling ratio is maintained to preserve texture details; the average information entropy of the water area sub-blocks is 0.8. If the value is below the threshold, a large-scale downsampling is used to compress the computational load. For each batch of forward-looking and forward-looking back-looking image pairs, a semi-global matching algorithm is used to calculate the disparity, outputting two sets of disparity fields.

[0105] In step S04, the two sets of disparity fields are input into the stationary wavelet depth decomposition manifold model. In this embodiment, the model employs a 4-level stationary wavelet transform decomposition, with each level generating one low-frequency sub-band and three directional high-frequency sub-bands, resulting in a total of 16 sub-band propagation channels. For example... Figure 3 As shown, the high-frequency subband propagation channel employs depthwise separable convolution, with the number of channels increasing exponentially with layer depth, and a base number of 64 channels. The low-frequency subband propagation channel uses a large receptive field convolution with a 7×7 receptive field, responsible for capturing macroscopic elevation trends. The cross-scale topological propagation matrix uses a 64×64-dimensional learnable linear transformation operator to achieve bidirectional feature fusion between subbands at each level. Local disparity gradient information from low-level high-frequency subbands is propagated upwards to high-level low-frequency subbands, while the global elevation trend from high-level low-frequency subbands is propagated downwards back to low-level high-frequency subbands. The continuous boundary convergence control equation uses the squared mean of the elevation gradient differences within 5 pixels on each side of the seam as the regularization loss term, with a regularization weight of 0.05. The trend of the seam elevation residuals on the validation set with training epochs is shown below. Figure 4 As shown, when the decrease is less than 10 times in 10 consecutive training rounds Training terminates at m, and the final convergence result shows that the elevation residual of the splice seam in the validation set is 0.38 m.

[0106] Perform step S05 to calculate the elevation difference field of 11 pairs of overlapping region digital surface models, extract three indicators, and calculate the mutual verification accuracy index. The weighting coefficients for the three indicators were determined through parameter search as follows: the absolute value of the elevation difference and the mean value. The weight is 0.40, and the standard deviation of the elevation difference is... The weight is 0.35, and the effective coverage of corresponding points is... The weight is 0.25. The cross-validation accuracy index for each overlapping region pair. The calculation results are shown in Table 3:

[0107] Table 3. Statistics on mutual verification accuracy index and triggering strategy for each overlapping region.

[0108]

[0109] As shown in Table 3, the mutual verification accuracy index of OL-05 is 0.32, which is lower than the first accuracy threshold of 0.48, triggering a global regional network re-adjustment and requiring a re-solution of the rational function model parameters for all scenes in this batch; OL-03, OL-07, and OL-10... Values ​​falling between 0.48 and 0.72 trigger a mild compensation strategy, superimposed with affine transformation compensation; other regions... A value of 0.72 or higher indicates acceptable matching quality, and adjustment compensation will not be triggered. For example... Figure 5 As shown, after OL-05 triggers global re-adjustment, the average absolute value of its elevation difference decreases from 2.87 m to 0.61 m, indicating that global regional network re-adjustment can effectively converge residual errors when the systematic geometric deviation is large, while simple affine correction is difficult to cover all nonlinear components when the deviation exceeds a certain order of magnitude.

[0110] In step S06, the eight digital surface models are stitched and fused batch by batch according to the overlapping region index table. The continuous boundary convergence control equation continuously applies gradient smoothing constraints within a 5-pixel band on both sides of the stitching seam. The elevation gradient jump variable at the stitching seam decreases from an average of 1.74 m / pixel before stitching to 0.29 m / pixel after stitching. The elevation gradient curve tends to flatten out rather than steepen near the stitching seam, as shown below. Figure 6 As shown. This trend is explained from a technical perspective: the convergence control equation for continuous boundaries uses the mean squared gradient difference as a penalty term, imposing implicit gradient consistency constraints on the elevation output of pixel bands on both sides of the stitching seam during the inference phase. This smooths the jagged faults remaining in traditional stitching methods due to the lack of gradient constraints into a continuously transitioning elevation field. The final output digital surface model is shown below. Figure 7 As shown, the splice seam is not visible in the elevation field, and the elevation field is continuous as a whole, which meets the judgment criteria for mutual verification.

[0111] like Figures 2 to 7 The accompanying drawings comprehensively record the key processing results and data change trends of each step of the present invention, which helps to verify the applicability of the method under various terrain conditions.

[0112] It should be noted that the variables involved in this invention are explained in detail in Tables 4 and 5.

[0113] Table 4. Variable Explanation Table (Part 1)

[0114]

[0115] Table 5. Variable Explanation Table (Part Two)

[0116]

[0117] 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 mutual verification of satellite DSM matching accuracy based on overlapping regions, characterized in that, Includes the following steps: Acquire multi-sensor calibrated images, select forward-looking, front-looking, and back-looking images with qualified radiation quality levels according to the radiation quality grading system, extract the vector boundary of the overlapping area range, and construct an overlapping area index table; The overlapping area images are fed into the group behavior adaptive streaming perception model to extract the spatial coordinate features and grayscale features of each image in the overlapping area, remove gross errors in real time, and output the feature point set of the overlapping area after removing gross errors. For the feature point set of overlapping areas after removing gross errors, image data is loaded into a streaming block memory cache queue. Combined with dynamic perception pyramid downsampling processing, the disparity fields of the two sets of epipolar images, namely the forward-looking and forward-looking back-looking images, are extracted in batches. Two sets of disparity fields are fed into a stationary wavelet depth decomposition manifold model. High-frequency disparity texture and low-frequency elevation trend are extracted from multiple wavelet subbands. The features of each subband are fused by a cross-scale topological transmission matrix to output the digital surface model results after elevation gradient alignment in the splicing seam region. Calculate the elevation difference field of the digital surface model of the overlapping area of ​​adjacent scenes, calculate the mutual verification accuracy index of the overlapping area based on the dynamic node transmission weight evaluation function, and trigger the corresponding adjustment compensation strategy when the mutual verification accuracy index falls into different intervals to update the rational function model parameters. All scene digital surface models after updating rational function model parameters are stitched and fused together. The elevation gradient at the stitching seam is constrained by the continuous boundary convergence control equation to smoothly transition, and the final mutually verified digital surface model result is output.

2. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 1, characterized in that, The specific rules for determining the radiation quality grading system are as follows: data with a cloud coverage area of ​​no more than 5% and no radiation anomalies are rated as Level 1; data with a cloud coverage area of ​​more than 5% but no more than 20% and no radiation anomalies, or cloud cover of no more than 20% and slight radiation anomalies, are rated as Level 2; data with cloud cover of more than 20% but no more than 60% and radiation anomalies not exceeding slight levels, or cloud cover of no more than 60% and moderate radiation anomalies, are rated as Level 3; data with cloud cover of more than 60% or severe radiation anomalies are rated as Level 4; Level 1 data is preferred, and if Level 1 data is insufficient, it is supplemented by downgrading to the next level.

3. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 2, characterized in that, The overlapping area index table records the spatial overlap range of adjacent scene images and the corresponding scene number using vector polygons, which is used for subsequent batch scheduling of overlapping area image data.

4. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 3, characterized in that, The structure of the adaptive streaming perception model for group behavior is as follows: the front layer is a spatially orthogonal decoupling channel mapping layer, which maps the input multi-temporal spatial coordinates and grayscale features to mutually orthogonal feature subspaces; The mapping layer is followed by a dynamic sparse perception layer, which embeds a learnable mask matrix and performs propagation pruning on neurons whose spatial uncertainty exceeds the spatial uncertainty threshold. The core inference layer adopts a multi-level self-circulating feedback architecture, modeling the hidden layer neurons as mimic particles. The particles simulate group collision avoidance behavior through mutual repulsion force operators. The output layer outputs the gross error probability vector of each feature point. Points with probabilities exceeding the gross error probability threshold are marked as gross errors and removed.

5. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 4, characterized in that, In the adaptive streaming sensing model for group behavior, when the predicted elevation value of a certain track deviates from the average value of the current batch of groups, the collision penalty term reduces the transmission weight of the corresponding node according to the quadratic function relationship of the deviation magnitude. The transmission weights of neurons between layers are adaptively allocated according to the memory usage ratio of the feature maps of each layer; the feature data between layers flows in memory in a circular queue, and the queue depth is dynamically adjusted according to the amount of data in the current batch.

6. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 5, characterized in that, The streaming block memory cache queue is as follows: image data is divided into several sub-blocks according to spatial blocks, and loaded into memory sequentially in a first-in-first-out queue. Each batch only retains the sub-block data required for the current processing and releases it immediately after processing, keeping the peak memory usage within the memory usage ratio threshold.

7. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 6, characterized in that, Dynamic perceptual pyramid downsampling specifically involves dynamically selecting the downsampling ratio based on the texture information entropy of the current sub-block during the image pyramid construction process. Sub-blocks with texture information entropy lower than the information entropy threshold adopt a larger downsampling ratio, while sub-blocks with texture information entropy higher than the information entropy threshold maintain a smaller downsampling ratio.

8. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 7, characterized in that, The structure of the stationary wavelet depth decomposition manifold model is as follows: The input layer receives two sets of epipolar image disparity fields: forward-looking and forward-looking back-looking; the first level is a multi-layer stationary wavelet transform decomposition layer, which cancels the downsampling operation and maintains the feature maps at each scale consistent with the input resolution; each wavelet decomposition generates one low-frequency subband and three high-frequency subbands in three directions; each subband establishes an independent manifold feature propagation channel; the deepest layer introduces a continuous boundary convergence control equation; the output layer reconstructs the digital surface model result by inverse transforming the features after aligning each subband.

9. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 8, characterized in that, The cross-scale topological transmission matrix is ​​specifically designed to establish learnable linear transformation operators between feature maps at different wavelet decomposition levels. Through matrix multiplication, the local disparity gradient information of the low-level high-frequency subband is transmitted to the high-level low-frequency subband, while the global elevation trend of the high-level low-frequency subband is transmitted inversely to the low-level high-frequency subband, thereby achieving bidirectional collaborative fusion of multi-scale features.

10. The method for mutual verification of satellite DSM matching accuracy based on overlapping regions according to claim 9, characterized in that, The continuous boundary convergence control equation is as follows: the mean square of the elevation gradient difference within the pixel band within the threshold range of the pixel band width on both sides of the splicing seam is used as the regularization loss term to ensure that the elevation output of the model at the splicing seam position satisfies the gradient continuity constraint.

Images