A multi-resolution four-dimensional meteorological data spatio-temporal fusion method

Abstract

The application relates to the field of data processing and discloses a multi-resolution four-dimensional meteorological data space-time fusion method, which comprises the following steps: acquiring multi-source meteorological data and uniformly mapping the multi-source meteorological data to a four-dimensional space-time grid framework to obtain a first data field and a second data field; identifying a spatial strong gradient area based on the second data field by using an adaptive multi-scale voting threshold value; if there is no strong gradient area, directly performing multi-resolution pyramid weighted fusion; if there is a strong gradient area, dividing the strong gradient area into a plurality of connected sub-areas which are continuous in gradient direction, calculating and dividing each sub-area into similar sub-regions, and respectively judging whether there is systematic spatial misplacement; for the sub-area with misplacement, constructing a global systematic misplacement field, and performing reverse deformation correction on the first data field; finally, constructing a multi-resolution pyramid representation of the registered data field and the second data field, performing weighted fusion by using an adaptive weight, and reconstructing to generate a four-dimensional meteorological fusion data body. The application can adaptively allocate computing resources according to data content, and accurately correct spatial misplacement in a strong gradient area.

Description

Technical Field

[0001] This invention relates to the field of data processing, specifically to a spatiotemporal fusion method for multi-resolution four-dimensional meteorological data. Background Technology

[0002] Meteorological data fusion integrates multi-source data from different observation platforms such as meteorological satellites, weather radars, and numerical weather prediction models into a unified, high-precision, and high spatiotemporal resolution data volume. It is a core technical component of numerical weather prediction, climate monitoring, and severe weather warning. However, due to inherent systematic errors, geometric positioning biases, differences in projection methods, and asynchronous data acquisition times among different sensors, multi-source meteorological data often exhibit varying degrees of spatial misalignment. This misalignment has a relatively small impact on the fusion results in regions with gentle gradients in meteorological elements (such as within homogeneous air masses), but it can lead to severe fusion artifacts in regions with strong gradients, such as frontal ghosting, blurred boundaries, and false double-front structures, directly reducing the accuracy of subsequent meteorological analysis and forecasting.

[0003] Existing meteorological data fusion methods are mainly divided into two categories: One type is image registration-based methods, which first estimate the full-field offset field using techniques such as cross-correlation, optical flow, or phase correlation, and then perform geometric correction on the data before fusion. This type has the following shortcomings: First, full-field registration is computationally intensive, and the registration results in uniform regions with poor texture are unreliable, yet still consume a lot of computational resources. Second, existing registration methods usually use a single similarity metric (such as cross-correlation) and a fixed registration strategy, which cannot adapt to the characteristics of huge local texture differences in meteorological fields, and are prone to producing erroneous offsets in low similarity regions.

[0004] Another type is to directly perform weighted fusion or variational fusion, without considering the spatial misalignment between data. Direct fusion completely ignores the misalignment problem, which will inevitably produce ghosting or blurring in regions of strong gradients. Summary of the Invention

[0005] The purpose of this invention is to provide a spatiotemporal fusion method for multi-resolution four-dimensional meteorological data, thereby solving at least one of the above-mentioned technical problems.

[0006] The objective of this invention can be achieved through the following technical solutions: A spatiotemporal fusion method for multi-resolution four-dimensional meteorological data includes the following steps: The multi-source meteorological data to be fused is acquired and uniformly mapped to a four-dimensional spatiotemporal grid framework to obtain the first data field and the second data field, where the four dimensions include a three-dimensional spatial dimension and a one-dimensional temporal dimension. Based on the second data field, by calculating the gradient magnitude of each grid point and using an adaptive multi-scale voting threshold, the region formed by grid points whose gradient magnitude exceeds the judgment condition is identified as a strong spatial gradient region; if there is no strong spatial gradient region, the first data field and the second data field are directly fused by multi-resolution pyramid weighted fusion to generate a four-dimensional meteorological fusion data volume. If a region of strong spatial gradient exists, then within the region of strong spatial gradient, it is determined whether there is a systematic spatial misalignment between the first data field and the second data field; If the systematic spatial misalignment exists, the systematic misalignment field of the first data field relative to the second data field is calculated, and the first data field is subjected to reverse deformation correction based on the multi-resolution pyramid using the systematic misalignment field to obtain the registered data field; if the systematic spatial misalignment does not exist, the first data field is directly used as the registered data field. A multi-resolution pyramid representation of the registered data field and the second data field is constructed. At each resolution level, an adaptive weighting function is used for weighted fusion, and then a four-dimensional meteorological fusion data volume is generated through pyramid reconstruction. The adaptive weighting function is as follows: at each resolution level, for each grid point, the weight of the first data field is equal to the normalized ratio of the sum of the reciprocal of the local variance of that point and the reciprocal of the local variance of the second data field.

[0007] As a further technical solution, the specific process for determining whether a systematic spatial misalignment exists includes: The spatial strong gradient region is divided into multiple connected sub-regions, and the following judgment process is performed on each connected sub-region: Calculate the local structural similarity map between the first data field and the second data field within the current connected sub-region, and further divide the sub-region into high similarity sub-regions, medium similarity sub-regions, and low similarity sub-regions based on the similarity value; For highly similar sub-regions, the offset field is calculated using the sliding window cross-correlation method. If the average amplitude of the offset field exceeds the preset spatial offset threshold, it is determined that there is a systematic misalignment in the connected sub-region. For similar sub-regions, the offset field is calculated using both the sliding window cross-correlation method and the phase difference method. If the difference between the offset fields obtained by the two methods is less than the preset difference tolerance, it is determined that there is a systematic misalignment in the connected sub-region, and the weighted average of the two is taken as the offset field of the sub-region. If the difference is greater than or equal to the tolerance, the sub-region is marked as an undetermined sub-region. For low similarity sub-regions, the sub-regions are directly marked as undetermined sub-regions; For all undetermined sub-regions, a unified multi-resolution cyclic iterative registration method is used for judgment: starting from the top of the pyramid, the initial offset field is set to zero field, and the offset field is refined layer by layer. If the converged offset field satisfies spatial continuity and the average amplitude exceeds the spatial offset threshold, it is determined that the corresponding connected sub-region has a systematic misalignment; otherwise, it is determined that it does not exist. If at least one connected sub-region is determined to have a systematic misalignment, then the spatial strong gradient region is ultimately determined to have a systematic spatial misalignment, and the offset field of all connected sub-regions with misalignment is output; otherwise, it is ultimately determined that there is no systematic spatial misalignment.

[0008] The tolerance for difference is determined as follows: Within the similar sub-regions, the local standard deviations of the first offset field obtained by the sliding window cross-correlation method and the second offset field obtained by the phase difference method are calculated respectively. The difference tolerance is set to a preset constant between 0.2 and 0.5 times the maximum value of the two local standard deviations.

[0009] As a further technical solution, the process for obtaining the difference value between the offset fields obtained by the two methods is as follows: The first step is to convert the first offset field obtained by the sliding window cross-correlation method and the second offset field obtained by the phase difference method into local direction fields and local phase fields respectively within the similar sub-regions, thereby obtaining the first direction field, the first phase field, the second direction field, and the second phase field. The second step is to calculate the cross-entropy between the first directional field and the second directional field, as well as the cross-entropy between the first phase field and the second phase field for each pixel; multiply the two cross-entropies and take the cube root to obtain the modal cross-entropy difference for each pixel. The third step is to use the local structural similarity map within the similar sub-regions as spatial weights to accumulate the modal cross-entropy differences to obtain the weighted accumulated differences. At the same time, the spatial variation coefficient of the weighted accumulated differences is calculated. If the spatial variation coefficient exceeds a preset threshold, the weighted accumulated differences are multiplied by the chaotic modulation factor. The value of the chaotic modulation factor is equal to the difference between the peak and valley values ​​of the local structural similarity map within the sub-region, divided by the mean. The fourth step is to perform a ratio calculation between the weighted cumulative difference after chaotic modulation and the joint information entropy of the two offset fields in the similar sub-region to obtain the difference value; wherein the joint information entropy is the product of the entropy of the first offset field and the entropy of the second offset field.

[0010] As a further technical solution, the process of obtaining the difference value between the offset fields obtained by the two methods also includes: Fifth, apply a random perturbation test to the difference values ​​obtained in step four: That is, within the similar sub-regions, several sub-windows are randomly selected, and a random offset is applied to the first offset field in each sub-window. Then, the cross-entropy change rate with the second offset field in that sub-window is recalculated. If the cross-entropy change rate of multiple sub-windows is positive and the average change rate exceeds the preset sensitivity threshold, then the difference value is determined to be reliable and output directly. Otherwise, multiply the difference value by the attenuation coefficient to obtain the final difference value.

[0011] As a further technical solution, when there are multiple connected sub-regions determined to have systematic misalignment, the method for calculating the systematic misalignment field is as follows: The offset fields of each connected sub-region are used as local constraints. For regions where no misalignment is detected, the offset is set to zero. Then, thin plate spline interpolation or radial basis function interpolation is used to generate a smooth and continuous global offset field covering the entire first data field, which serves as the systematic misalignment field.

[0012] As a further technical solution, the specific process of multi-resolution iterative registration is as follows: Starting from the top of the pyramid, the initial value of the offset field of the current layer is set to the previous transfer value (zero field at the top layer). The offset field of the current layer is used to pre-correct the first data field, and the cross-correlation peak value between the pre-corrected data field and the second data field in the undetermined sub-region is calculated. If the increase in the cross-correlation peak exceeds the decay tolerance coefficient of the increase in the previous iteration, the offset field is passed to the next fine layer. If the increase is less than the preset ratio of the initial cross-correlation peak value of the layer, then revert to the previous layer and apply spatial smoothing constraints before recalculating; Iterate until the change in offset field between two consecutive iterations is less than the moving average of the change. If the final cross-correlation peak is greater than the preset success threshold, then a systematic misalignment is determined to exist; otherwise, it is determined not to exist.

[0013] As a further technical solution, the process of dividing the spatial strong gradient region into multiple connected sub-regions includes: Calculate the gradient direction field of the second data field within the strong gradient region of the space; Using grid points whose gradient magnitude exceeds the adaptive threshold as seed points, the gradient direction consistency connectivity criterion is used for region growth: only when the gradient direction angle between adjacent grid points is less than the direction threshold and their gradient magnitudes are all higher than the adaptive threshold, they are assigned to the same sub-region. After growth is complete, unvisited strong gradient points are used as new seed points for repeated growth until all points are partitioned. Each connected region obtained by growth is treated as a connected sub-region. The gradient direction within the sub-region is continuous, and different sub-regions are separated by discontinuous directions or low gradient boundaries.

[0014] As a further technical solution, the following is also included before the region growth: For the second data field in the strong gradient region of space, the gradient magnitude field at three scales, namely the original resolution, downsampled by 2 times, and downsampled by 4 times, is calculated respectively. The adaptive threshold Ts at each scale is calculated independently. Then, the gradient magnitude fields at the three scales are unified back to the original resolution through bilinear interpolation. The number of votes v∈{0,1,2,3} for each grid point at each scale is taken as the number of votes for whether it exceeds the corresponding threshold Ts. The final gradient magnitude threshold criterion for each grid point is defined as follows: If v≥2, then the point is determined to satisfy the condition that the gradient magnitude is higher than the threshold. If v=1, then further check whether there are at least 2 points with v≥2 in the neighborhood of this point. If there are, it is determined that the condition is satisfied; otherwise, it is not satisfied. If v=0, then the condition is not met; The grid points that meet the above conditions are sorted according to their v values ​​from highest to lowest. Points with v=3 are designated as first-level seed points, points with v=2 are designated as second-level seed points, and points with v=1 that pass the neighborhood check are designated as third-level seed points. During regional growth, the process starts with the first-level seed point. Once all the first-level seed points have grown, the second- and third-level seed points are used as new seed points to continue growth.

[0015] As a further technical solution, the preset direction threshold is obtained as follows: During the region growing process, the local standard deviation of the gradient direction angle of the grid points already included is calculated in real time for the currently growing connected sub-regions. ; The current orientation threshold for this sub-region The expression is: ; Wherein, is the baseline threshold, with a value ranging from 5° to 15°. The minimum permissible orientation threshold ranges from 3° to 5°. The maximum allowable orientation threshold is set between 20° and 30°.

[0016] The beneficial effects of this invention are: (1) This invention first identifies strong gradient regions in space and decides whether to perform misalignment detection and correction accordingly. When there are no strong gradient regions, pyramid fusion is performed directly, realizing adaptive matching of computing resources and data content. Compared with the existing technology, which performs full-field registration or direct fusion indiscriminately in all regions, this invention can avoid invalid registration calculations in flat regions and concentrate computing power on strong gradient regions for accurate correction. It greatly improves computing efficiency while ensuring fusion quality, and is especially suitable for practical application scenarios where strong gradient features are sparsely distributed in large-scale meteorological fields. (2) This invention divides the strong gradient region in space into multiple connected sub-regions with continuous gradient directions, and further divides them into high, medium and low similarity sub-regions based on local structural similarity. For different layers, the sliding window cross-correlation method, the dual-method cross-validation method or the multi-resolution iterative registration method are used to achieve adaptive adaptation to the local texture heterogeneity in the meteorological field. Compared with the existing technology that uses a single registration method to run the entire field, this invention can quickly and efficiently judge the texture-rich sub-regions, and adopt a more robust global optimization strategy in the texture-poor or structurally complex sub-regions, avoiding the problem of local methods failing in difficult areas and improving the accuracy of systematic misalignment judgment. (3) This invention designs a difference quantification method based on the ratio of cross-entropy of direction field, cross-entropy of phase field, spatial variation coefficient, chaotic modulation factor and joint information entropy, which is used to measure the consistency between the offset fields obtained by the sliding window cross-correlation method and the phase difference method. Compared with the existing technology that only uses simple statistical quantities such as mean square error or correlation coefficient, this invention can distinguish between differences caused by random noise and systematic structural differences, and adaptively adjust the difference sensitivity according to the spatial distribution of local structural similarity. Finally, the difference threshold can be transferred across scenarios by normalizing the joint information entropy, providing an accurate and reliable quantitative basis for the credible judgment of similar sub-regions. Attached Figure Description

[0017] The invention will now be further described with reference to the accompanying drawings.

[0018] Figure 1 This is a diagram illustrating the method steps of the present invention. Detailed Implementation

[0019] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] Please see Figure 1As shown, this invention is a spatiotemporal fusion method for multi-resolution four-dimensional meteorological data, comprising the following steps: The process involves acquiring multi-source meteorological data to be fused and mapping it uniformly to a four-dimensional spatiotemporal grid framework to obtain the first and second data fields. The four dimensions include a three-dimensional spatial dimension and a one-dimensional temporal dimension. Multi-source meteorological data (such as satellite, radar, and numerical weather prediction models) have different spatial resolutions, temporal resolutions, and coordinate reference systems. This step unifies the original data into a standard four-dimensional (three-dimensional spatial + temporal) grid framework, providing a comparable and computable benchmark for subsequent fusion. This eliminates the heterogeneity between data and enables data from different sensors to be compared and algebraically calculated point by point in the same spatiotemporal coordinate system, avoiding spurious errors introduced by coordinate inconsistencies.

[0021] Based on the second data field, by calculating the gradient magnitude of each grid point and applying an adaptive multi-scale voting threshold, the region formed by grid points whose gradient magnitude exceeds the judgment condition is identified as a strong spatial gradient region; the specific method of the adaptive multi-scale voting threshold is as follows: For the second data field in the strong gradient region of space, the gradient magnitude field at three scales, namely the original resolution, downsampled by 2 times, and downsampled by 4 times, is calculated respectively. The adaptive threshold Ts at each scale is calculated independently. At each scale, Ts = μ + g * σ, where μ and σ are the mean and standard deviation of the gradient magnitude field at that scale, respectively, and g is a preset coefficient with a value range of 1.0 to 2.0, and the default value is 1.5.

[0022] Region growing relies on seed points, but the determination of seed points (whether the gradient magnitude exceeds the adaptive threshold) presents a fundamental contradiction at a single scale: at small scales, the gradient magnitude is sensitive to noise and prone to false positives; at large scales, the gradient magnitude decreases after smoothing, making it easy to miss small strong gradient bands. In this scheme, a multi-scale voting mechanism is adopted to calculate the gradient magnitude and adaptive threshold Ts (Ts=μ+gσ, g defaults to 1.5) at three scales: the original resolution, downsampling by 2 times, and downsampling by 4 times. Three scales were chosen because the spatial scale of meteorological features typically ranges from a few kilometers to hundreds of kilometers, and 2x and 4x downsampling can cover medium and macro scales. The gradient magnitude fields at the three scales were then restored to the original resolution using bilinear interpolation. The number of times each grid point exceeded its respective scale threshold, ν (0–3), was counted. Only points with ν ≥ 2 were ultimately identified as strong gradient points, meaning they exhibited strong gradients at at least two scales, effectively suppressing single-scale noise. Points with ν = 1 needed to be checked to see if there were at least two points with ν ≥ 2 in their 3×3 neighborhood; only those with ν ≥ 2 were accepted. Weak gradient points were allowed as connection points, but they required strong gradient support from their surroundings. Points meeting the criteria were divided into three levels of seed points based on their ν values: ν = 3 was level one (most reliable), ν = 2 was level two, and ν = 1, passing the neighborhood check, was level three. During region growth, level one seed points were used first, and after all level one seed points were completed, level two and level three seed points were used as new seeds to continue growth. This priority order ensures that growth starts from the most reliable core region and expands outward, avoiding erroneous growth caused by noisy seed points. Compared with the single-scale threshold of existing technologies, this invention can simultaneously ensure high recall (no missed detections) and high precision (no false detections) in strong gradient regions, playing a robust role in identifying seed points.

[0023] Then, the gradient magnitude fields at the three scales are unified back to the original resolution through bilinear interpolation, and the number of votes v∈{0,1,2,3} for each grid point at each scale is taken to determine whether it exceeds the corresponding threshold Ts. The final gradient magnitude threshold criterion for each grid point is defined as follows: If v≥2, then the point is determined to satisfy the condition that the gradient magnitude is higher than the threshold. If v=1, then further check whether there are at least 2 points with v≥2 in the neighborhood (3×3 window) of this point. If there are, it is determined that the condition is satisfied; otherwise, it is not satisfied. If v=0, then the condition is not met; The grid points that meet the above conditions are sorted according to their v values ​​from highest to lowest. Points with v=3 are designated as first-level seed points, points with v=2 are designated as second-level seed points, and points with v=1 that pass the neighborhood check are designated as third-level seed points. During region growth, the process starts with the first-level seed points. Once all the first-level seed points have been grown, the second- and third-level seed points are used as new seed points to continue growing. This avoids noise points caused by excessively low threshold values ​​becoming invalid seeds, and also avoids missing key connection points in weak gradient bands due to excessively high threshold values.

[0024] Key structures in meteorological fields (such as fronts, jet streams, and convective clouds) typically correspond to regions with large gradient amplitudes. However, single-scale gradient thresholds are easily affected by noise or weak gradient bands, leading to missed or over-detection. This step independently calculates adaptive thresholds at three scales (original, 2x downsampled, and 4x downsampled), and then uses a voting mechanism to comprehensively determine whether each grid point belongs to a strong gradient region. The multi-scale voting can simultaneously capture gradient features at different scales, avoiding misjudging small-scale noise as strong gradients and preventing the omission of small strong gradient bands after large-scale smoothing. The three-level seed point (ν=3,2,1) sorting and growth strategy ensures the robustness of region growth and prevents region fragmentation due to inappropriate single-scale thresholds. This solves the problem that traditional fixed gradient thresholds or single-scale adaptive thresholds cannot accurately identify multi-scale and multi-morphological strong gradient regions in complex meteorological fields, especially the problem that key connection points within weak gradient bands are easily missed.

[0025] If there is no strong spatial gradient region, the first data field and the second data field are directly fused using multi-resolution pyramid weighted fusion to generate a four-dimensional meteorological fusion data volume. When the entire meteorological field is very flat (such as inside a uniform air mass) and there is no strong gradient region, spatial misalignment detection and correction become meaningless or may even introduce errors. This branch skips misalignment processing and enters pyramid fusion directly, which can avoid invalid calculations, improve the processing efficiency of the method in featureless regions, prevent erroneous correction due to false misalignment, and solve the problem that existing fusion methods still force registration when there are no significant textures or gradient regions, which may produce negative effects.

[0026] If a region of strong spatial gradient exists, then within the region of strong spatial gradient, it is determined whether there is a systematic spatial misalignment between the first data field and the second data field; If the systematic spatial misalignment exists, the systematic misalignment field of the first data field relative to the second data field is calculated, and the first data field is subjected to reverse deformation correction based on the multi-resolution pyramid using the systematic misalignment field to obtain the registered data field; if the systematic spatial misalignment does not exist, the first data field is directly used as the registered data field. A multi-resolution pyramid representation of the registered data field and the second data field is constructed. At each resolution level, an adaptive weighting function is used for weighted fusion, and then a four-dimensional meteorological fusion data volume is generated through pyramid reconstruction. The adaptive weighting function is as follows: at each resolution level, for each grid point, the weight of the first data field is equal to the normalized ratio of the sum of the reciprocal of the local variance of that point and the reciprocal of the local variance of the second data field.

[0027] A hidden challenge in meteorological data fusion is that visual artifacts caused by spatial misalignment only appear in regions with strong gradients (such as fronts and dry lines). Within uniform air masses, even a shift of a few pixels is difficult for the human eye and subsequent numerical calculations to detect. Therefore, registration is not a necessity for the entire field, but rather a tool specific to regions with strong gradients. Existing technologies either force registration across the entire field (which is computationally intensive, and registration in flat areas is inherently unreliable due to the lack of texture leading to inaccurate shift estimation) or directly fuse data to ignore misalignment (ghosting occurs at fronts). In this scheme, strong gradient regions in space are first adaptively identified based on the second data field (reference field). If no strong gradient regions are found, it means that the entire field is a low-risk area, and the misalignment detection and correction are skipped directly, and multi-resolution pyramid fusion is entered. This saves computing power and avoids meaningless registration in textureless areas. If strong gradient regions exist, subsequent judgment is made. This framework of partitioning first and then making decisions essentially focuses computing resources on the spatial location with the highest information content. Compared with existing technologies that either cut registration outright or cut off registration outright, this invention can dynamically distribute data according to the data content, which improves the overall efficiency of the fusion system.

[0028] The specific process for determining whether a systematic spatial misalignment exists includes: The spatial strong gradient region is divided into multiple connected sub-regions, and the following judgment process is performed on each connected sub-region: Calculate the local structural similarity map between the first data field and the second data field within the current connected sub-region. The local structural similarity map uses the SSIM (Structural Similarity) algorithm to calculate the similarity of brightness, contrast, and structure for each pixel and its neighborhood. The product of these three factors is taken as the similarity value of the pixel, with a value range of [0,1]. Then, based on the similarity value, the sub-region is divided into high similarity sub-regions (similarity ≥ 0.7), medium similarity sub-regions (0.4 ≤ similarity < 0.7), and low similarity sub-regions (similarity < 0.4). For highly similar sub-regions, the offset field is calculated using the sliding window cross-correlation method. If the average amplitude of the offset field exceeds the preset spatial offset threshold, it is determined that there is a systematic misalignment in the connected sub-region. For similar sub-regions, the offset field is calculated using both the sliding window cross-correlation method and the phase difference method. If the difference between the offset fields obtained by the two methods is less than the preset difference tolerance, it is determined that there is a systematic misalignment in the connected sub-region, and the weighted average of the two is taken as the offset field of the sub-region. If the difference is greater than or equal to the tolerance, the sub-region is marked as an undetermined sub-region. For low similarity sub-regions, the sub-regions are directly marked as undetermined sub-regions; For all undetermined sub-regions, a unified multi-resolution cyclic iterative registration method is used for judgment: starting from the top of the pyramid, the initial offset field is set to zero field, and the offset field is refined layer by layer. If the converged offset field satisfies spatial continuity and the average amplitude exceeds the spatial offset threshold, it is determined that the corresponding connected sub-region has a systematic misalignment; otherwise, it is determined that it does not exist. The spatial continuity refers to the smoothness of local changes in the offset field. Specifically, it is determined as follows: the magnitude of the offset vector difference between any two adjacent grid points in the offset field is less than 1.5 times the standard deviation of the offset magnitude in that local area, and there are no isolated jump points in the entire field (a jump point is defined as the difference between the offset mean of its 8 neighboring areas and the local standard deviation is more than 3 times).

[0029] If at least one connected sub-region is determined to have a systematic misalignment, then the spatial strong gradient region is ultimately determined to have a systematic spatial misalignment, and the offset field of all connected sub-regions with misalignment is output; otherwise, it is ultimately determined that there is no systematic spatial misalignment.

[0030] The tolerance for difference is determined as follows: Within the similar sub-regions, the local standard deviations of the first offset field obtained by the sliding window cross-correlation method and the second offset field obtained by the phase difference method are calculated respectively. The difference tolerance is set to a preset constant between 0.2 and 0.5 times the maximum value of the two local standard deviations.

[0031] In regions with strong gradients, misalignment detection still faces challenges, namely, varying texture richness at different locations, making a single registration method inapplicable. For example, in the core region of a front (high similarity), the sliding window cross-correlation method is accurate and fast enough; however, at the edge or broken regions of a front (medium to low similarity), the cross-correlation method may be affected by texture changes within the window, leading to biases, while the phase difference method is sensitive to sub-pixel shifts but easily affected by phase wrapping interference. In this proposed solution, the strong gradient region is first divided into multiple connected sub-regions (each corresponding to a meteorological feature boundary) according to the gradient direction continuity. Then, the local structural similarity is calculated within each sub-region, further dividing it into high, medium, and low similarity sub-regions. High similarity sub-regions: Cross-correlation is used directly due to its high signal-to-noise ratio, allowing for rapid judgment. Medium similarity sub-regions: Both cross-correlation and phase difference methods are used simultaneously. If the results are close (difference less than the tolerance), a weighted average is taken as the offset field, as the two independent methods mutually verify each other, resulting in high reliability. If the difference is large, it indicates complex texture or occlusion in the region, making both cross-correlation and phase difference methods unreliable, and the region is marked as an undetermined sub-region. Low-similarity sub-regions are directly marked as undetermined sub-regions. All undetermined sub-regions are uniformly assigned to multi-resolution iterative registration. This scheme is a global optimization method that refines layer by layer from a coarse scale, avoiding local extrema. Finally, as long as one connected sub-region is determined to have a misalignment, the entire strong gradient region is considered to have a systematic misalignment, and the offset field of all misaligned sub-regions is output. The purpose of the above design is to handle simple regions with the lightest method, medium regions with cross-validation, and difficult regions with the most robust but expensive method, achieving Pareto optimality in accuracy and computational cost. Compared with existing technologies that use a single method to run across the entire field, this invention can adaptively allocate algorithm resources, playing a dual role in improving judgment accuracy and computational efficiency.

[0032] The process of obtaining the difference between the offset fields obtained by the two methods is as follows: The first step is to convert the first offset field obtained by the sliding window cross-correlation method and the second offset field obtained by the phase difference method into local direction fields and local phase fields respectively within the similar sub-regions, thereby obtaining the first direction field, the first phase field, the second direction field, and the second phase field. The local orientation field is calculated as follows: for each pixel, calculate the components (u, v) of the offset in the x and y directions, then the orientation angle θ = arctan(v / u) of that point, and the orientation field is the unit orientation vector (cosθ, sinθ) of each pixel; the local phase field is calculated as follows: treat the offset field as a complex field (real part is u, imaginary part is v), perform a Hilbert transform, and take the argument of the complex number as the local phase.

[0033] The second step is to calculate the cross-entropy between the first and second directional fields, and the cross-entropy between the first and second phase fields, for each pixel. The cross-entropy is calculated as follows: the directional angle is discretized into 36 intervals (each in 10° increments), and the probability distributions P and Q are obtained by statistically analyzing the proportion of pixels in each directional interval. The phase field is also discretized into 36 intervals before the cross-entropy is calculated.

[0034] Multiply the two cross-entropies and take the cube root to obtain the modal cross-entropy difference for each pixel. The third step is to use the local structural similarity map within the similar sub-regions as spatial weights to accumulate the modal cross-entropy differences in a weighted manner to obtain the weighted accumulated difference. At the same time, the spatial variation coefficient of the weighted accumulated difference is calculated. If the spatial variation coefficient exceeds a preset threshold (the preset threshold is 0.3), the weighted accumulated difference is multiplied by the chaotic modulation factor. The value of the chaotic modulation factor is equal to the difference between the peak and valley values ​​of the local structural similarity map within the sub-region, divided by the mean; the peak value refers to the maximum value of the local structural similarity values ​​of all grid points within the sub-region, the valley value refers to the minimum value, and the mean value refers to the arithmetic mean.

[0035] The fourth step is to perform a ratio calculation between the weighted cumulative difference after chaotic modulation and the joint information entropy of the two offset fields in the similar sub-region to obtain the difference value; wherein the joint information entropy is the product of the entropy of the first offset field and the entropy of the second offset field. The entropy of the offset field is estimated using an equal-width histogram. The offset amplitude range is divided into K equal intervals, where K is the cube root of the total number of grid points in the sub-region, rounded down. The percentage of pixels in each interval is then calculated. Then entropy .

[0036] Within similar subregions, how can we quantify the difference between two offset fields obtained by the cross-correlation method and the phase difference method? Traditional methods, such as mean squared error (MSE), can only reflect numerical differences but cannot distinguish between differences caused by random noise and systematic structural differences. For example, in regions with repetitive textures (such as periodic clouds), two offset fields may be numerically similar, but their orientation field distributions are completely different. This difference is structural and should be amplified. In uniform regions, however, numerical differences may be large, but the orientation field distribution is random. This difference is noise-driven and should be suppressed. In this scheme, the first step is to convert the offset fields into local orientation fields and local phase fields. The orientation field reflects the direction of the offset, and the phase field reflects the fine structure of the offset. The second step is to calculate the cross-entropy of the two orientation fields and the cross-entropy of the two phase fields. Cross-entropy measures the difference between two probability distributions. Here, the orientation angle is discretized into 36 intervals (each in 10° increments), and the distribution is obtained by statistically analyzing the percentage of pixels in each interval. The larger the cross-entropy, the more inconsistent the distributions of the two fields. The third step involves multiplying the two cross-entropies and taking the cube root to obtain the modal cross-entropy difference. Multiplying by the cube root balances the contributions of the direction and phase fields, preventing noise from one field from dominating the result. The local structural similarity map is used as a spatial weight for weighted accumulation, yielding a weighted accumulated difference. The reason for weighting is that in regions with high similarity, the two offset fields should be highly consistent, and the difference weight should be large; in regions with low similarity, even large differences do not necessarily indicate a real problem. Simultaneously, the spatial coefficient of variation (standard deviation / mean) of the weighted accumulated difference is calculated. If the coefficient of variation > 0.3, it indicates a highly uneven distribution of differences. In this case, a chaotic modulation factor is multiplied, which equals the peak-to-valley difference of the structural similarity map within the sub-region divided by the mean. Essentially, the unevenness of structural similarity is used to modulate the difference value: if the structural similarity changes drastically, it indicates a complex texture in the region, requiring amplification of the difference sensitivity; if the change is gradual, the original difference is maintained. The fourth step involves calculating the ratio between the modulated weighted accumulated difference and the joint information entropy (entropy product) of the two offset fields. The reason for using joint information entropy is that the larger the entropy (the more chaotic) of the offset field itself, the larger the absolute value of its difference will be. Dividing by the entropy product can normalize it, making the difference value unaffected by the complexity of the offset field itself, and facilitating the uniform setting of tolerance thresholds. The internal logic of the entire calculation chain is: starting from the difference in direction / phase structure, using local similarity weighting, using the coefficient of variation and chaotic modulation factor to adaptively adjust the sensitivity, and using joint information entropy for normalization, finally obtaining a difference value that is adaptive to texture changes, robust to noise, and comparable across scenes. Compared with the simple MSE or correlation coefficient of existing technologies, this invention can truly measure the structural consistency of the two registration results, playing a role in accurately judging whether similar sub-regions are reliable.

[0037] The process of obtaining the difference between the offset fields obtained by the two methods also includes: Fifth, apply a random perturbation test to the difference values ​​obtained in step four: That is, within the similar sub-regions, several sub-windows are randomly selected, and a random offset is applied to the first offset field in each sub-window. Then, the cross-entropy change rate with the second offset field in that sub-window is recalculated. The number of sub-windows is one percent of the total number of grid points in the sub-region and is at least 5. Each sub-window is 5×5 grid points in size. The magnitude of the random offset is uniformly distributed in the range of [-dmax, dmax], where dmax is 0.5 times the standard deviation of the first offset field in the sub-window.

[0038] If the cross-entropy change rate of multiple sub-windows is positive and the average change rate exceeds the preset sensitivity threshold, then the difference value is determined to be reliable and output directly. Otherwise, multiply the difference value by the attenuation coefficient to obtain the final difference value; The decay coefficient is equal to the power of the negative of the mean absolute value of the cross-entropy change rate of all sub-windows, with the natural constant e as the base, and is therefore always less than 1.

[0039] Because the difference values ​​can still be deceived by accidental consistency—that is, in regions of repetitive texture (such as large-scale ocean waves or uniform cloud streets), cross-correlation and phase difference methods may calculate numerically close offset fields due to periodic pattern matching—this closeness is fragile. Applying a small perturbation to the first offset field will break the match, and the cross-entropy will increase significantly. True systematic misalignment is different: if there is indeed an overall translation of the first data field relative to the second data field, then applying a small perturbation to the first offset field should result in a small change in cross-entropy (because the perturbation breaks the correct match). In this scheme, multiple sub-windows are randomly selected within the similar sub-regions (the number is 1% of the total number of grid points, but at least 5, and the window size is 5×5). A random offset (with an amplitude uniformly distributed within 0.5 times the standard deviation of the first offset field within the window) is applied to the first offset field within each sub-window, and the rate of change of cross-entropy between the first and second offset fields within that window is recalculated. If more than half of the sub-windows have a positive cross-entropy change rate (i.e., the difference increases after perturbation) and the average change rate exceeds the sensitivity threshold, it indicates that the original consistency is fragile and accidental, and the difference value is unreliable. In this case, the difference value is multiplied by a decay coefficient less than 1. Multiplying by a decay coefficient less than 1 further reduces the difference value. The smaller the difference between the offset fields obtained by the two methods, the easier it is to meet the condition of less than the preset difference tolerance, thus determining that there is a systematic misalignment in the similar sub-regions. This means that when random perturbation testing reveals that the original consistency is fragile and accidental (i.e., unreliable), this scheme actively adjusts the difference value to a smaller direction, prompting the sub-region to be judged as having a misalignment. Therefore, it is better to mark a region with complex texture and accidental matching as having a misalignment and send it to the subsequent multi-resolution cyclic iterative registration (this registration method is more robust to difficult regions) than to erroneously judge that there is no misalignment due to the false reliability of the difference value, and thus directly output an unreliable offset field. In other words, this attenuation mechanism implements a conservative strategy: in situations of uncertainty, it prioritizes triggering a more rigorous global registration process rather than blindly trusting the consistency of local bi-methods. Compared to existing technologies that directly use static difference values ​​without stability verification, this invention effectively identifies and handles the risk of accidental consistency in repetitive texture regions through the combination of random perturbation testing and attenuation coefficients, thereby improving the judgment of systematic misalignment.

[0040] When there are multiple connected sub-regions that are determined to have systematic misalignment, the method for calculating the systematic misalignment field is as follows: The offset fields of each connected sub-region are used as local constraints. For regions where no misalignment is detected, the offset is set to zero. Then, thin plate spline interpolation or radial basis function interpolation is used to generate a smooth and continuous global offset field covering the entire first data field, which serves as the systematic misalignment field.

[0041] When multiple connected sub-regions are determined to have systematic misalignment, each sub-region has its own local offset field. These offset fields may differ in direction and magnitude, and may even be contradictory (e.g., one adjacent sub-region offsets to the right while the other offsets to the left). A simple approach in existing technologies is to directly stitch the offset fields of each sub-region together, setting the offset fields of undetected misalignment regions to zero, and then using linear interpolation at the boundaries. This results in discontinuous deformation fields, leading to tearing or folding at the sub-region boundaries during reverse deformation correction. In this scheme, the offset fields of each connected sub-region are used as local constraint points (i.e., at these points, the deformation value is known). For regions where misalignment is not detected (including all grid points outside of strong gradient regions), the offset is set to zero as a boundary constraint. Then, thin plate spline (TPS) interpolation or radial basis function (RBF) interpolation is used. TPS interpolation minimizes the bending energy function, generating a deformation field with a second-order continuous derivative over the entire domain. Physically, this is equivalent to the elastic deformation of an infinitely large thin plate at a given constraint point. This deformation field does not produce tearing or folding, and the transition between constraint points is natural and smooth. RBF interpolation, on the other hand, is achieved through a linear combination of radial basis functions, offering better local control.

[0042] The choice between these methods depends on the specific data characteristics: if the constraint points are sparse and uniformly distributed, TPS is better; if the constraint points are dense and local variations are drastic, RBF is more flexible. The global offset field generated by this method not only ensures high accuracy within sub-regions but also guarantees smooth transitions between sub-regions and between sub-regions and misaligned regions. Compared to the simple stitching or linear interpolation of existing technologies, this invention can generate a physically reasonable global deformation field, thus preventing the introduction of additional artifacts during reverse deformation correction.

[0043] The specific process of multi-resolution cyclic iterative registration is as follows: Starting from the top of the pyramid, the initial value of the offset field of the current layer is set to the previous transfer value (zero field at the top layer). The offset field of the current layer is used to pre-correct the first data field, and the cross-correlation peak value between the pre-corrected data field and the second data field in the undetermined sub-region is calculated. If the increase in the cross-correlation peak exceeds the decay tolerance coefficient of the increase in the previous iteration (initially set to 0.5), then the offset field is passed to the next fine layer; If the increase is less than the preset ratio of the initial cross-correlation peak value of the layer, such as 10%, then revert to the previous layer and apply spatial smoothing constraints before recalculating. Iterate until the change in offset field between two consecutive iterations is less than the moving average of the change. If the final cross-correlation peak value is greater than the preset success threshold (e.g., 0.7, i.e., the cross-correlation peak value), then a systematic misalignment is determined to exist; otherwise, it is determined not to exist.

[0044] In this scheme, the approximate offset is first estimated at a low resolution (coarse layer), and then refined layer by layer. The specific process is as follows: Starting from the top of the pyramid, the initial value of the offset field of the current layer is set to the value passed from the previous layer (the top layer is zero field). The first data field is pre-corrected using the current offset field, and then the cross-correlation peak value between the pre-corrected data field and the second data field in the undetermined sub-region is calculated. A higher cross-correlation peak value indicates better alignment. If the increase in the cross-correlation peak value after this iteration (compared to the previous iteration) exceeds the decay tolerance coefficient of the increase in the previous iteration (initially 0.5), it indicates effective convergence, and the propagation continues to the next fine layer. If the increase is less than 10% of the initial cross-correlation peak value of that layer, it indicates that convergence has stalled or diverged. At this point, it is backed up to the previous layer, and a spatial smoothing constraint (e.g., Gaussian filtering) is applied before recalculation. The smoothing constraint can suppress erroneous gradients caused by noise, forcing the offset field to adjust in a more continuous direction. Iteration continues until the change in the offset field between two adjacent iterations is less than the moving average of the changes (i.e., convergence). Finally, if the cross-correlation peak value is greater than 0.7, a reliable offset field is considered to have been found, and a misalignment is determined; otherwise, it is determined that no misalignment exists. This scheme employs an adaptive step size control through a decay tolerance coefficient and a backoff mechanism to prevent iterations from diverging under noise. Spatial smoothing constraints are applied during backoff, which is equivalent to introducing regularization when encountering difficulties. Compared with existing technologies that use fixed number of iterations or single-resolution gradient descent, this invention can automatically adapt to the difficulty of the data, achieving stable convergence in difficult regions.

[0045] The process of dividing the spatial strong gradient region into multiple connected sub-regions includes: Calculate the gradient direction field of the second data field within the strong gradient region of the space; Using grid points whose gradient magnitude exceeds the adaptive threshold as seed points, the gradient direction consistency connectivity criterion is used for region growth: only when the gradient direction angle between adjacent grid points is less than the direction threshold and their gradient magnitudes are all higher than the adaptive threshold, they are assigned to the same sub-region. The adaptive threshold is calculated as follows: calculate the mean and standard deviation of the gradient magnitude of all grid points in the strong gradient region of space, and set the adaptive threshold as the weighted sum of the mean and standard deviation. The weight coefficient is adaptively adjusted according to the skewness of the gradient magnitude: if the skewness is positive, the weight coefficient is 1.5; if the skewness is negative, the weight coefficient is 0.5.

[0046] After growth is complete, unvisited strong gradient points are used as new seed points for repeated growth until all points are partitioned. Each connected region obtained by growth is treated as a connected sub-region. The gradient direction within the sub-region is continuous, and different sub-regions are separated by discontinuous directions or low gradient boundaries.

[0047] A low gradient boundary is a strip-shaped region between two adjacent connected sub-regions, where the gradient magnitude of each grid point within the strip-shaped region is less than 0.5 times the adaptive threshold, and the width of the strip-shaped region is at least 2 grid points, and the length is at least 5 grid points. When such a low gradient strip-shaped region exists between two strong gradient sub-regions, the low gradient strip-shaped region serves as the low gradient boundary separating the two sub-regions.

[0048] The gradient direction consistency connectivity criterion further includes: mapping the gradient direction to the interval [0°, 180°), using the Euclidean distance of the direction cosine to replace the angle difference, and determining that the directions are consistent when the direction cosine Euclidean distance between adjacent grid points is less than the direction consistency threshold. The direction consistency threshold is the larger of 0.5 times the standard deviation of the gradient direction in the sub-region and a preset minimum threshold. During the region growing process, for boundary points that are determined to have the same direction but whose gradient magnitude is lower than the adaptive threshold, if there are at least two already grown high gradient neighboring points around them, then the boundary point is included in the current sub-region as a buffer. After growth is complete, morphological closing operations are performed on each connected sub-region to fill the internal voids, and isolated regions with an area smaller than the minimum connected area are removed. The minimum connected area is set as the product of one-thousandth of the total number of grids in the spatial strong gradient region and the area of ​​a single grid.

[0049] Dividing strong gradient regions into connected sub-regions presents several problems: First, two strong gradient bands may be very close together, with a weak gradient transition band in between. Traditional connected component labeling based on amplitude thresholds might merge them into the same region (because the amplitude of the weak gradient band is below the threshold and does not constitute a separator), even though they are actually two independent meteorological features (e.g., two fronts with different orientations). Second, the gradient direction within the same meteorological feature (e.g., a curved front) can change continuously. If the direction threshold is fixed, the curved section might be cut off. This scheme first calculates the gradient direction field, then uses points with gradient amplitudes exceeding an adaptive threshold as seed points. The region growth criterion requires that the directional angle be less than the directional threshold and that the amplitude be greater than the adaptive threshold. Note that "greater than" means that every point on the growth path must meet the amplitude condition, ensuring that growth does not cross low-gradient regions. However, for true low-gradient boundaries, even if there is only a small weak gradient band between two strong gradient regions, growth will naturally stop because the amplitude is below the threshold, thus achieving separation. Furthermore, a quantitative standard for low-gradient boundaries is defined: the gradient amplitude is less than 0.5 times the adaptive threshold, and the width is ≥2 grids and the length is ≥5 grids, making the separation criterion operable. In addition, the directional threshold is dynamically adjusted to adapt to curved boundaries. After growth, the gradient direction within each connected sub-region is continuous, corresponding to a meteorological feature boundary; different sub-regions are either separated by abrupt directional changes (e.g., the intersection of two fronts) or by low-gradient bands. Compared to existing technologies that rely solely on amplitude thresholds, this invention can correctly separate strong gradient bands with different orientations while ensuring the integrity of curved boundaries, thus accurately segmenting meteorological features.

[0050] The preset direction threshold is obtained as follows: During the region growing process, the local standard deviation of the gradient direction angle of the grid points already included is calculated in real time for the currently growing connected sub-regions. The calculation window for the local standard deviation is all grid points within the minimum bounding rectangle of the currently grown points in the connected sub-region; The current orientation threshold for this sub-region The expression is: ; in, The baseline threshold is set at 5° to 15°. The minimum permissible orientation threshold ranges from 3° to 5°. The maximum allowable orientation threshold is set between 20° and 30°.

[0051] If the orientation threshold in region growing is fixed, setting it too small will cause curved fronts to be cut into multiple small fragments; setting it too large will cause gradient bands in different directions to be incorrectly merged. This scheme allows the orientation threshold to adapt to the degree of change in the gradient direction of the points already grown in the current sub-region; specifically, it calculates the local standard deviation of the gradient orientation angle of the points already included in the grid in real time. (The calculation window is the smallest bounding rectangle of the currently grown points). This reflects the degree of curvature of the current sub-region: Small indicates that the direction is consistent and the area is close to a straight line; A large value indicates a significant change in direction, suggesting the region is curved. Direction threshold. The formula is .in It is the baseline threshold (5°~15°), when When it is very small, θ≈ (but was) Clamping (not lower than 3°–5°), maintaining a tight grip; when When it increases, Linear increase, but not exceeding (20°~30°), to prevent unlimited relaxation. The essence of this formula is: to use the local standard deviation as an estimate of directional change, adding it to the baseline threshold so that the directional threshold matches the actual curvature of the sub-region. For example, the curvature of a straight front... ≈0, ≈ =10°, the angle between adjacent grid points must be less than 10° for growth to be guaranteed; a curved arc-shaped front, It may reach 10°, then The maximum allowable angle between adjacent points is 20°, allowing for continuous growth of curved boundaries. Compared to existing technologies with fixed thresholds, this invention can adaptively handle both straight and curved features, enabling the region growing algorithm to maintain strict segmentation of straight features while accommodating complete connections of curved features.

[0052] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A spatiotemporal fusion method for multi-resolution four-dimensional meteorological data, characterized in that, Includes the following steps: The multi-source meteorological data to be fused is acquired and uniformly mapped to a four-dimensional spatiotemporal grid framework to obtain the first data field and the second data field, where the four dimensions include a three-dimensional spatial dimension and a one-dimensional temporal dimension. Based on the second data field, by calculating the gradient magnitude of each grid point and using an adaptive multi-scale voting threshold, the region formed by grid points whose gradient magnitude exceeds the judgment condition is identified as a strong spatial gradient region; if there is no strong spatial gradient region, the first data field and the second data field are directly fused by multi-resolution pyramid weighted fusion to generate a four-dimensional meteorological fusion data volume. If a region of strong spatial gradient exists, then within the region of strong spatial gradient, it is determined whether there is a systematic spatial misalignment between the first data field and the second data field; If the systematic spatial misalignment exists, the systematic misalignment field of the first data field relative to the second data field is calculated, and the first data field is subjected to reverse deformation correction based on the multi-resolution pyramid using the systematic misalignment field to obtain the registered data field; if the systematic spatial misalignment does not exist, the first data field is directly used as the registered data field. A multi-resolution pyramid representation of the registered data field and the second data field is constructed. An adaptive weighting function is used for weighted fusion at each resolution level, and then a four-dimensional meteorological fusion data volume is generated through pyramid reconstruction.

2. The spatiotemporal fusion method for multi-resolution four-dimensional meteorological data according to claim 1, characterized in that, The specific process for determining whether a systematic spatial misalignment exists includes: The spatial strong gradient region is divided into multiple connected sub-regions, and the following judgment process is performed on each connected sub-region: Calculate the local structural similarity map between the first data field and the second data field within the current connected sub-region, and further divide the sub-region into high similarity sub-regions, medium similarity sub-regions, and low similarity sub-regions based on the similarity value; For highly similar sub-regions, the offset field is calculated using the sliding window cross-correlation method. If the average amplitude of the offset field exceeds the preset spatial offset threshold, it is determined that there is a systematic misalignment in the connected sub-region. For similar sub-regions, the offset field is calculated using both the sliding window cross-correlation method and the phase difference method. If the difference between the offset fields obtained by the two methods is less than the preset difference tolerance, it is determined that there is a systematic misalignment in the connected sub-region, and the weighted average of the two is taken as the offset field of the sub-region. If the difference is greater than or equal to the tolerance, the sub-region is marked as an undetermined sub-region. For low similarity sub-regions, the sub-regions are directly marked as undetermined sub-regions; For all undetermined sub-regions, a unified multi-resolution cyclic iterative registration method is used for judgment: starting from the top of the pyramid, the initial offset field is set to zero field, and the offset field is refined layer by layer. If the converged offset field satisfies spatial continuity and the average amplitude exceeds the spatial offset threshold, it is determined that the corresponding connected sub-region has a systematic misalignment; otherwise, it is determined that it does not exist. If at least one connected sub-region is determined to have a systematic misalignment, then the spatial strong gradient region is ultimately determined to have a systematic spatial misalignment, and the offset field of all connected sub-regions with misalignment is output; otherwise, it is ultimately determined that there is no systematic spatial misalignment.

3. The spatiotemporal fusion method for multi-resolution four-dimensional meteorological data according to claim 2, characterized in that, The process of obtaining the difference between the offset fields obtained by the two methods is as follows: The first step is to convert the first offset field obtained by the sliding window cross-correlation method and the second offset field obtained by the phase difference method into local direction fields and local phase fields respectively within the similar sub-regions, thereby obtaining the first direction field, the first phase field, the second direction field, and the second phase field. The second step is to calculate the cross-entropy between the first directional field and the second directional field, as well as the cross-entropy between the first phase field and the second phase field for each pixel; multiply the two cross-entropies and take the cube root to obtain the modal cross-entropy difference for each pixel. The third step is to use the local structural similarity map within the similar sub-regions as spatial weights to accumulate the modal cross-entropy differences to obtain the weighted accumulated differences. At the same time, the spatial variation coefficient of the weighted accumulated differences is calculated. If the spatial variation coefficient exceeds a preset threshold, the weighted accumulated differences are multiplied by the chaotic modulation factor. The value of the chaotic modulation factor is equal to the difference between the peak and valley values ​​of the local structural similarity map within the sub-region, divided by the mean. The fourth step is to perform a ratio calculation between the weighted cumulative difference after chaotic modulation and the joint information entropy of the two offset fields in the similar sub-region to obtain the difference value; wherein the joint information entropy is the product of the entropy of the first offset field and the entropy of the second offset field.

4. The spatiotemporal fusion method for multi-resolution four-dimensional meteorological data according to claim 3, characterized in that, The process of obtaining the difference between the offset fields obtained by the two methods also includes: Fifth, apply a random perturbation test to the difference values ​​obtained in step four: That is, within the similar sub-regions, several sub-windows are randomly selected, and a random offset is applied to the first offset field in each sub-window. Then, the cross-entropy change rate with the second offset field in that sub-window is recalculated. If the cross-entropy change rate of multiple sub-windows is positive and the average change rate exceeds the preset sensitivity threshold, then the difference value is determined to be reliable and output directly. Otherwise, multiply the difference value by the attenuation coefficient to obtain the final difference value.

5. The spatiotemporal fusion method for multi-resolution four-dimensional meteorological data according to claim 2 or 3, characterized in that, When there are multiple connected sub-regions that are determined to have systematic misalignment, the method for calculating the systematic misalignment field is as follows: The offset fields of each connected sub-region are used as local constraints. For regions where no misalignment is detected, the offset is set to zero. Then, thin plate spline interpolation or radial basis function interpolation is used to generate a smooth and continuous global offset field covering the entire first data field, which serves as the systematic misalignment field.

6. The spatiotemporal fusion method for multi-resolution four-dimensional meteorological data according to claim 2, characterized in that, The specific process of multi-resolution cyclic iterative registration is as follows: Starting from the top of the pyramid, the initial value of the offset field of the current layer is set to the previous value. The offset field of the current layer is used to pre-correct the first data field, and the cross-correlation peak value between the pre-corrected data field and the second data field in the undetermined sub-region is calculated. If the increase in the cross-correlation peak exceeds the decay tolerance coefficient of the increase in the previous iteration, the offset field is passed to the next fine layer. If the increase is less than the preset ratio of the initial cross-correlation peak value of the layer, then revert to the previous layer and apply spatial smoothing constraints before recalculating; Iterate until the change in offset field between two consecutive iterations is less than the moving average of the change. If the final cross-correlation peak is greater than the preset success threshold, then a systematic misalignment is determined to exist; otherwise, it is determined not to exist.

7. The spatiotemporal fusion method for multi-resolution four-dimensional meteorological data according to claim 2, characterized in that, The process of dividing the spatial strong gradient region into multiple connected sub-regions includes: Calculate the gradient direction field of the second data field within the strong gradient region of the space; Using grid points whose gradient magnitude exceeds the adaptive threshold as seed points, the gradient direction consistency connectivity criterion is used for region growth: only when the gradient direction angle between adjacent grid points is less than the direction threshold and their gradient magnitudes are all higher than the adaptive threshold, they are assigned to the same sub-region. After growth is complete, unvisited strong gradient points are used as new seed points for repeated growth until all points are partitioned. Each connected region obtained by growth is treated as a connected sub-region. The gradient direction within the sub-region is continuous, and different sub-regions are separated by discontinuous directions or low gradient boundaries.

8. The spatiotemporal fusion method for multi-resolution four-dimensional meteorological data according to claim 7, characterized in that, Prior to region growth, the following is also included: For the second data field in the strong gradient region of space, the gradient magnitude field at three scales, namely the original resolution, downsampled by 2 times, and downsampled by 4 times, is calculated respectively. The adaptive threshold Ts at each scale is calculated independently. Then, the gradient magnitude fields at the three scales are unified back to the original resolution through bilinear interpolation. The number of votes v∈{0,1,2,3} for each grid point at each scale is taken to determine whether it exceeds the corresponding threshold Ts. The final gradient magnitude threshold criterion for each grid point is defined as follows: If v≥2, then the point is determined to satisfy the condition that the gradient magnitude is higher than the threshold. If v=1, then further check whether there are at least 2 points with v≥2 in the neighborhood of this point. If there are, it is determined that the condition is satisfied; otherwise, it is not satisfied. If v=0, then the condition is not met; The grid points that meet the above conditions are sorted according to their v values ​​from highest to lowest. Points with v=3 are designated as first-level seed points, points with v=2 are designated as second-level seed points, and points with v=1 that pass the neighborhood check are designated as third-level seed points. During regional growth, the process starts with the first-level seed point. Once all the first-level seed points have grown, the second- and third-level seed points are used as new seed points to continue growth.

9. The spatiotemporal fusion method for multi-resolution four-dimensional meteorological data according to claim 8, characterized in that, The preset direction threshold is obtained as follows: During the region growing process, the local standard deviation of the gradient direction angle of the grid points already included is calculated in real time for the currently growing connected sub-regions. ; The current orientation threshold for this sub-region The expression is: ； in, The baseline threshold is set at 5° to 15°. The minimum permissible orientation threshold ranges from 3° to 5°. The maximum allowable orientation threshold is set between 20° and 30°.

Images