An atmospheric data fusion method based on wavelet transform

Abstract

The application provides a wavelet transform-based atmospheric data fusion method, and relates to the field of satellite data processing. The method takes FY-4B data as a unified space-time base plate, performs two-dimensional discrete wavelet transform on the spatially overlapped regions of FY-4B and FY-3 series data to obtain wavelet coefficients of each scale of each satellite data; performs time enhancement processing on the wavelet coefficients of each scale of the FY-4B data, and performs inverse discrete wavelet transform on the wavelet coefficients of each scale of each satellite data to reconstruct each satellite image signal; calculates the comprehensive scores of each satellite data to determine the initial fusion weight, and corrects the final fusion weight through an intelligent credibility switching mechanism; and performs weighted summation on each satellite image signal according to the final fusion weight to generate a surface reflectivity product. The application can fuse the complementary advantages of geostationary satellites and polar-orbiting satellites, and improve the spatio-temporal resolution and detail fidelity of multi-source satellite data fusion.

Description

Technical Field

[0001] This application relates to the field of satellite data processing technology, and in particular to an atmospheric data fusion method based on wavelet transform. Background Technology

[0002] Data fusion is a core technology for integrating multi-source observation and numerical model data to compensate for the spatiotemporal coverage deficiencies of single sensors. It is also a fundamental support for numerical weather prediction, severe weather monitoring, and climate research. However, for in-depth analysis of the environmental and climate effects of dust aerosols, most existing studies rely on single satellite observations or traditional interpolation methods for collaborative research. The inherent contradiction between the insufficient spatial resolution of my country's FY-4B geostationary satellite and the temporal coverage gaps of the FY-3D / FY-3F / FY-3E polar-orbiting satellites means that a single satellite cannot accurately capture the rapid transport and local outbreak processes of dust aerosols, further accumulating errors between simulated and observed data. In addition, traditional fusion methods often use smoothing and noise reduction when fusing multi-source satellite data, which can easily lead to the loss of fine features such as local edges and the spread of boundary outliers, making it difficult to balance the continuity of large-scale background with the fidelity of small-scale details. Summary of the Invention

[0003] This application provides an atmospheric data fusion method based on wavelet transform, aiming to conduct spatiotemporal complementary fusion of multi-source satellite data based on multi-scale analysis techniques using wavelet transform. It combines discrete wavelet multi-scale decomposition, multi-index adaptive weight allocation, and regional overlap fusion algorithms to overcome the inherent contradictions of smooth transition and detail loss in current multi-source satellite data fusion processes, such as insufficient spatial resolution of geostationary satellites and time coverage gaps of polar-orbiting satellites. This method combines the multi-scale time-frequency decomposition capability of discrete wavelet transform with the advantages of multi-index adaptive weight allocation. Furthermore, it introduces a component-level differential fusion mechanism and a sudden change mechanism to detect and dynamically correct the fusion weights of components at different scales. This will help provide a more reliable, refined, and continuous data solution for atmospheric remote sensing.

[0004] To achieve the above objectives, this application adopts the following technical solution: This application provides an atmospheric data fusion method based on wavelet transform, the method comprising: Acquire multi-source satellite data, which includes geostationary meteorological satellite data and polar-orbiting meteorological satellite data. The geostationary meteorological satellite data includes FY-4B data, and the polar-orbiting meteorological satellite data includes FY-3D data, FY-3F data, and FY-3E data. Using the FY-4B data as a unified spatiotemporal base, the spatial overlap region between the FY-3D, FY-3F, and FY-3E data and the FY-4B data is determined. Two-dimensional discrete wavelet transforms are then performed on the data within the spatial overlap region for each of the FY-4B, FY-3D, FY-3F, and FY-3E data, decomposing each satellite data step-by-step into low-frequency approximate components and multiple high-frequency detail components, thus obtaining the wavelet coefficients at each scale for the FY-4B data, the FY-3D data, the FY-3F data, and the FY-3E data. Time-enhanced processing is performed on the wavelet coefficients of each scale of the FY-4B data to obtain enhanced FY-4B wavelet coefficients at each scale. Inverse discrete wavelet transform is then performed on the enhanced FY-4B wavelet coefficients at each scale to reconstruct the enhanced FY-4B image signal. Inverse discrete wavelet transform is then performed on the wavelet coefficients of each scale of the FY-3D data, the FY-3F data, and the FY-3E data, respectively, to reconstruct the FY-3D image signal, FY-3F image signal, and FY-3E image signal registered with the FY-4B data. Based on a multi-dimensional dynamic weight calculation system, the comprehensive scores of the FY-4B data, the FY-3D data, the FY-3F data, and the FY-3E data are calculated respectively. Based on the comprehensive score of each satellite data, the initial fusion weight of each satellite data is calculated; The initial fusion weights are corrected according to the intelligent reliability switching mechanism to determine the final fusion weights of each satellite data. Based on the final fusion weights, the enhanced FY-4B image signal, the FY-3D image signal, the FY-3F image signal, and the FY-3E image signal are weighted and summed to generate a high spatiotemporal resolution surface reflectance product.

[0005] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, two-dimensional discrete wavelet transform is performed on the FY-4B, FY-3D, FY-3F, and FY-3E data within the spatially overlapping region, respectively, to decompose each satellite data into low-frequency approximate components and multiple high-frequency detail components, thereby obtaining wavelet coefficients at each scale for each satellite data, including: Use separate low-pass filters h and high-pass filter g Two-dimensional image signals were processed in both row and column directions. f ( x , yConvolution and downsampling are performed to decompose the two-dimensional image signal into a low-frequency approximate component. LL 1 and three high-frequency detail components, the three high-frequency detail components including a horizontal detail component. LH 1. Vertical detail components HL 1 and diagonal detail components HH 1. The calculation process is expressed as follows: in, LL 1( m , n ) represents the coefficient value of the low-frequency approximate component at position (m,n) after the first decomposition. LH 1( m , n Let be the coefficient value of the horizontal detail component at position (m,n) after the first decomposition. HL 1( m , n ) represents the coefficient value of the vertical detail component at position (m,n) after the first decomposition. HH 1( m , n ) represents the coefficient value of the diagonal detail component at position (m,n) after the first decomposition. i This is the row index variable for the filter coefficients. j For the column index variable of the filter coefficients, h ( i Let be the coefficient value of the low-pass filter at the i-th position. h ( j Let be the coefficient value of the low-pass filter at the j-th position. g ( i Let be the coefficient value of the high-pass filter at the i-th position. g ( j Let be the coefficient value of the high-pass filter at the j-th position. f (2 m - i ,2 n - j ) represents the pixel value of the original two-dimensional image at coordinates (2m-i, 2n-j); A recursive decomposition method is used to approximate the low-frequency components at each level. LL k The decomposition continues until the preset number of decomposition layers is reached, resulting in the final low-frequency approximation components and all high-frequency detail components at each layer. These final low-frequency approximation components and all high-frequency detail components at each layer are then used as wavelet coefficients for each scale of the corresponding satellite data.

[0006] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, the calculation formula for the inverse discrete wavelet transform is: in, LL 1(2 m - i ,2 n - j ) represents the coefficient value of the first-layer low-frequency approximation component at coordinates (2m-i, 2n-j). LH 1(2 m - i ,2 n - j ) represents the coefficient value of the first-level horizontal detail component at coordinates (2m-i, 2n-j). HL 1(2 m - i ,2 n - j ) represents the coefficient value of the first layer vertical detail component at coordinates (2m-i, 2n-j). HH 1(2 m - i ,2 n - j ) represents the coefficient value of the first-layer diagonal detail component at coordinates (2m-i, 2n-j). This is the satellite image signal reconstructed by inverse discrete wavelet transform.

[0007] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, the wavelet coefficients of FY-4B data at each scale are subjected to time enhancement processing to obtain enhanced FY-4B wavelet coefficients at each scale, including: Calculate the Gaussian weights at each time point within the time window. w t The calculation formula is: in, t 0 represents the current moment. T For all moments within the time window, exp is the natural exponential function. t It is a time variable; Based on the Gaussian weights, the wavelet coefficients at each scale are weighted and averaged using the following formula to obtain the standard-enhanced FY-4B wavelet coefficients at each scale. : in, C k,t For the first t At the [time]th moment k Wavelet coefficients at scale; For regions deemed untrustworthy by the aforementioned intelligent trust switching mechanism, median filtering is applied to the high-frequency detail coefficients to obtain strongly filtered and enhanced FY-4B wavelet coefficients at each scale. The formula for calculating the median filter is as follows: in, median This is the median function, which sorts the wavelet coefficients at all times within the time window and takes the median value.

[0008] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, based on a multi-dimensional dynamic weight calculation system, the comprehensive scores of the FY-4B data, the FY-3D data, the FY-3F data, and the FY-3E data are calculated respectively, including: Calculate hourly time stability index, daily time stability index, spatial information complexity index, reference matching index, and sudden change penalty coefficient; The comprehensive score of each satellite data is calculated based on the hourly time stability index, daily time stability index, spatial information complexity index, reference matching index, and sudden change penalty coefficient. The overall score of the FY-4B data is obtained by weighting the hourly time stability index, the daily time stability index, and the spatial information complexity index according to the first preset weight, the second preset weight, and the third preset weight, respectively, and the sum of the first preset weight, the second preset weight, and the third preset weight is 100%. The comprehensive score of any satellite data among the FY-3D, FY-3F, and FY-3E data is obtained by multiplying the daily time stability index, the spatial information complexity index, and the reference matching index by the fourth, fifth, and sixth preset weights, respectively, and then by the sudden change penalty coefficient of the corresponding satellite data. The sum of the fourth, fifth, and sixth preset weights is 100%.

[0009] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, the calculation of hourly time stability index, daily time stability index, spatial information complexity index, reference matching degree index, and sudden change penalty coefficient includes: Time stability index S t ( X The formula for calculating ) is: in, σ ( X ) is a data sequence X standard deviation μ( X ) is a data sequence X The mean value; the time stability index is applied to the data sequence within the hourly time window and the data sequence within the daily time window respectively, to obtain the hourly time stability index and the daily time stability index respectively; The formula for calculating the spatial information complexity index is as follows: in, V ( X ()( x , y Let be a uniform filtering function, representing the local variance at pixel (x,y) after applying the uniform filtering function to the image X. V For local variance images, Ss ( X ()( x , y ) is an index of the spatial information complexity of image X at pixel (x,y). x , y ) represents the pixel coordinates in the image, and min is the function to find the minimum value. i , j () represents the pixel coordinate offset within the local window. V ( X ()( i , j The value of the local variance image at offset coordinate (i,j) is given. This is a local window region, that is, a neighborhood range centered at (x,y); The formula for calculating the reference matching degree index is as follows: in, M ( X , Y ) represents the reference matching degree index, and max is the maximum value function. X The satellite data sequence to be evaluated. Y For reference data sequence; The formula for calculating the sudden change penalty coefficient is as follows: in, α To determine the severity of the punishment, β To detect sensitivity, d ( X , t 0) represents the current time. t The relative difference between a data value of 0 and the mean of the time window.

[0010] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, the formula for calculating the relative difference is: in, For the current moment t The data value of 0, The data value from the previous moment. This represents the data value at the next moment.

[0011] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, the initial fusion weight of each satellite data is calculated according to the comprehensive score of each satellite data, including: The initial fusion weights are calculated using the following formula. w i : in, s i For the first i The overall score of each satellite s 4b The overall score for the FY-4B satellite. s 3d The overall score for the FY-3D satellite. s 3f This is the overall score for the FY-3F satellite.

[0012] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, the initial fusion weights are corrected according to the intelligent reliability switching mechanism, including: The maximum relative difference between the FY-3 series data and the FY-4B data is calculated using the following formula. D max : in, FY 3 D , FY 4 B , FY 3 F , FY 3 E These are the numerical values ​​of their respective satellite data in the corresponding regions; when D max When the value is ≤0.25, the FY-3 series satellite data in the corresponding region is determined to be reliable. The initial fusion weight remains unchanged, and the initial fusion weight is used as the final fusion weight for that region. when D maxWhen the value is greater than 0.25, the FY-3 series satellite data in the corresponding area is determined to be unreliable. The fusion weight of the FY-4B data is forcibly increased, and the area is switched to a strong filtering enhancement mode. The increased fusion weight of the FY-4B data is a preset weight threshold. The fusion weights of the FY-3D, FY-3F, and FY-3E data in the area are redistributed according to the proportion of their initial fusion weights to obtain the final fusion weight of the area.

[0013] Preferably, in the above-mentioned atmospheric data fusion method based on wavelet transform, a high spatiotemporal resolution surface reflectance product is generated using the following formula. F : in, W 4b The final fusion weights for FY-4B data. FY 4 B enhanced For the enhanced FY-4B image signal, W 3d , W 3f , W 3e These are the final fusion weights for the FY-3D, FY-3F, and FY-3E data, respectively. FY 3 D aligned , FY 3 F aligned , FY 3 E aligned These are the FY-3D, FY-3F, and FY-3E image signals, respectively, registered with FY-4B data.

[0014] The atmospheric data fusion method based on wavelet transform provided in this application has at least the following beneficial effects: This application constructs a multi-source meteorological satellite data fusion scheme based on wavelet multi-scale analysis. Using FY-4B data as a unified spatiotemporal base, fusion calculations are performed only within the spatially overlapping regions of FY-3D, FY-3F, FY-3E, and FY-4B data. By introducing db4 wavelet three-level multi-scale decomposition and reconstruction technology into the multi-source satellite data fusion framework, two-dimensional discrete wavelet transform is performed on FY-4B and FY-3 series data. This separates the satellite observation signals into low-frequency approximate components and multiple high-frequency detail components, allowing large-scale background information and local texture details to be processed at their respective scales. This effectively avoids the spatial detail loss and boundary outlier diffusion problems caused by smoothing and noise reduction in traditional fusion methods.

[0015] Specifically, in terms of time enhancement, a weighted average of wavelet coefficients at various scales is applied by calculating Gaussian weights for each moment within the time window to achieve standard time enhancement. Simultaneously, for regions deemed unreliable based on credibility criteria, median filtering is used for strong filtering enhancement of high-frequency detail coefficients, effectively suppressing impulse noise and outliers and improving the robustness of the fusion results. Regarding weight allocation, hourly-level time stability indicators, daily-level time stability indicators, spatial information complexity indicators, reference matching degree indicators, and abrupt change penalty coefficients are calculated. Based on these multi-dimensional indicators, a comprehensive score for each satellite data is calculated, and the initial fusion weight for each satellite data is then calculated based on the comprehensive score, achieving adaptive optimization of the fusion weights. Furthermore, by calculating the maximum relative difference (Dmax) between the FY-3 series data and the FY-4B data, when Dmax is greater than 0.25, the FY-3 series data in that region is deemed unreliable. The fusion weight of the FY-4B data is forcibly increased, and the region is switched to a strong filtering enhancement mode. This effectively corrects fusion errors caused by time gaps and observation misalignments in polar-orbiting satellites, strictly ensuring the continuity of long-term data series.

[0016] This application combines the aforementioned regional overlap fusion strategy, db4 wavelet three-layer multi-scale decomposition and reconstruction, dual-mode time enhancement, multi-dimensional dynamic weight calculation, and intelligent credibility switching mechanism to balance the continuity of large-scale background with the fidelity of small-scale details, while also reducing the problems of poor spatiotemporal adaptability and data distortion in the fusion process. Attached Figure Description

[0017] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this application and, together with the description, serve to explain the principles of this application.

[0018] Figure 1 A technical roadmap for an atmospheric data fusion method based on wavelet transform provided for embodiments of this application; Figure 2 A flowchart illustrating an atmospheric data fusion method based on wavelet transform, provided as an embodiment of this application.

[0019] The accompanying drawings illustrate specific embodiments of this application, which will be described in more detail below. These drawings and descriptions are not intended to limit the scope of the concept in any way, but rather to illustrate the concepts of this application to those skilled in the art through reference to specific embodiments. Detailed Implementation

[0020] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. When the following description relates to the drawings, unless otherwise indicated, the same numbers in different drawings denote the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this application. Rather, they are merely examples of apparatuses and methods consistent with some aspects of this application as detailed in the appended claims.

[0021] The collection, storage, use, processing, transmission, provision, and disclosure of relevant data and information in the technical solution of this application all comply with the provisions of relevant laws and regulations and do not violate public order and good morals.

[0022] The technical solution of this application and how the technical solution of this application solves the above-mentioned technical problems are described in detail below with specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments. The embodiments of this application will now be described with reference to the accompanying drawings.

[0023] This application provides an atmospheric data fusion method based on wavelet transform. It aims to create a simpler data fusion method combining wavelet transform and data fusion, built using Python. Targeting visible light channel data from the FY-4B geostationary meteorological satellite and the FY-3D / FY-3F / FY3E polar-orbiting meteorological satellites, it achieves spatiotemporal complementary fusion of multi-source satellite data through wavelet transform multi-scale decomposition and reconstruction technology. The regional overlap fusion strategy, combined with intelligent weight allocation, time window filtering, and abrupt change detection mechanisms, generates high spatiotemporal resolution surface reflectance products at 15-minute intervals. This method can be widely applied in atmospheric environmental monitoring, high-frequency dynamic observation of severe weather such as sandstorms, cloud monitoring, aerosol inversion, and ecological environment change analysis, providing high-precision data support for weather forecasting and disaster prevention and mitigation.

[0024] like Figure 1The diagram illustrates the technical roadmap for an atmospheric data fusion method based on wavelet transform, provided in this application. This method addresses the complementary characteristics of the FY-4B geostationary meteorological satellite (high temporal resolution but insufficient spatial detail) and the FY-3D / FY-3F / FY3E polar-orbiting meteorological satellites (high spatial resolution but limited temporal coverage). It designs a hybrid framework using FY-4B as a unified spatiotemporal base and performing multi-scale fusion only in spatially overlapping regions. Specifically, it employs a db4 wavelet three-level decomposition multi-resolution analysis technique to decompose different satellite data into low-frequency approximate components and high-frequency detail components for separate processing. Simultaneously, it introduces multi-dimensional dynamic weight calculation, abrupt anomaly detection, and a reliability switching mechanism. While strictly ensuring the continuity of the time series, it maximizes the fusion of high spatial detail information from the FY3 series satellites, effectively solving the problems of boundary outlier diffusion and temporal discontinuity inherent in traditional fusion methods.

[0025] Specifically, please combine Figure 2 As shown, the atmospheric data fusion method based on wavelet transform includes the following steps S10-S70.

[0026] S10: Acquire multi-source satellite data, including geostationary meteorological satellite data and polar-orbiting meteorological satellite data. Geostationary meteorological satellite data includes FY-4B data, and polar-orbiting meteorological satellite data includes FY-3D data, FY-3F data, and FY-3E data.

[0027] In this embodiment, the multi-source satellite data includes geostationary meteorological satellite data and polar-orbiting meteorological satellite data. The geostationary meteorological satellite data is FY-4B data. The FY-4B satellite operates in geostationary orbit and can perform continuous high-frequency observations of a fixed area, offering high temporal resolution but relatively low spatial resolution. The polar-orbiting meteorological satellite data includes FY-3D, FY-3F, and FY-3E data. These satellites operate in sun-synchronous polar orbits, offering high spatial resolution but with gaps in temporal coverage. All of the above satellite data are visible light channel reflectance data that has undergone radiometric calibration and geometric correction, and can be stored in HDF or NetCDF format.

[0028] S20: Using FY-4B data as a unified spatiotemporal base, determine the spatial overlap region between FY-3D, FY-3F, and FY-3E data and FY-4B data; perform two-dimensional discrete wavelet transform on the data of FY-4B, FY-3D, FY-3F, and FY-3E data in the spatial overlap region, decompose each satellite data into low-frequency approximate components and multiple high-frequency detail components, and obtain the wavelet coefficients of each scale of FY-4B data, FY-3D data, FY-3F data, and FY-3E data.

[0029] In this embodiment, the spatially overlapping region refers to the area geographically covered by the FY-3D, FY-3F, and FY-3E polar-orbiting satellite data and the FY-4B geostationary satellite data. Two-dimensional discrete wavelet transform is a multi-scale signal analysis tool. It decomposes the image level by level through convolution and downsampling in the row and column directions using low-pass and high-pass filters, respectively, into low-frequency approximate components reflecting large-scale background trends and high-frequency detail components representing local details such as edges and textures. The low-frequency approximate components are thumbnails of the original image after low-pass filtering and downsampling, preserving the overall outline and main information of the image; the high-frequency detail components are the result of high-pass filtering of the original image, containing horizontal, vertical, and diagonal details. Wavelet coefficients at each scale refer to the sum of all low-frequency approximate coefficients and high-frequency detail coefficients obtained from the above decomposition at each decomposition level.

[0030] In some embodiments, for any two-dimensional image signal f ( x , y Its first-order two-dimensional discrete wavelet can be filtered by a separate low-pass filter. h and high-pass filter g Convolution and downsampling are performed in the row and column directions respectively, decomposing it into a low-frequency approximate component. LL 1 and three high-frequency detail components: Horizontal detail sharing LH 1. Vertical details sharing HL 1 and diagonal detail components HH 1. Its mathematical expression is: in, LL 1( m , n ) represents the coefficient value of the low-frequency approximate component at position (m,n) after the first decomposition. LH 1( m , n Let be the coefficient value of the horizontal detail component at position (m,n) after the first decomposition. HL 1( m , n ) represents the coefficient value of the vertical detail component at position (m,n) after the first decomposition. HH 1( m , n ) represents the coefficient value of the diagonal detail component at position (m,n) after the first decomposition. i This is the row index variable for the filter coefficients. j For the column index variable of the filter coefficients, h ( i Let be the coefficient value of the low-pass filter at the i-th position. h (j Let be the coefficient value of the low-pass filter at the j-th position. g ( i Let be the coefficient value of the high-pass filter at the i-th position. g ( j Let be the coefficient value of the high-pass filter at the j-th position. f (2 m - i ,2 n - j ) represents the pixel value of the original two-dimensional image at coordinates (2m-i, 2n-j).

[0031] A recursive decomposition method is used to approximate the low-frequency components at each level. LL k The decomposition continues until the preset number of decomposition layers L=3 is reached, ultimately yielding an LL3 low-frequency component and 3L high-frequency detail components. The final low-frequency approximation component and all high-frequency detail components at each level are used as wavelet coefficients for each scale of the corresponding satellite data.

[0032] S30: Perform time-enhanced processing on the wavelet coefficients of each scale of the FY-4B data to obtain enhanced FY-4B wavelet coefficients at each scale. Then, perform inverse discrete wavelet transform on the enhanced FY-4B wavelet coefficients at each scale to reconstruct the enhanced FY-4B image signal. Perform inverse discrete wavelet transform on the wavelet coefficients of each scale of the FY-3D data, the FY-3F data, and the FY-3E data to reconstruct the FY-3D image signal, FY-3F image signal, and FY-3E image signal registered with the FY-4B data.

[0033] In this embodiment, temporal enhancement processing leverages the good temporal continuity of geostationary satellite data, smoothing and optimizing the current moment's data using multiple frames of data within the time window before and after the current moment. Because the FY-4B satellite imaging time interval is short, the reflectivity changes in the same area within a short period are relatively gradual; multi-moment fusion can effectively suppress random noise. Gaussian weighted averaging assigns weights based on the temporal distance between each moment and the current moment; the closer the moment, the greater the weight, and the farther away, the smaller the weight. Median filtering uses the median of the wavelet coefficients at each moment within the time window as the filtering result, effectively suppressing impulse noise and outliers. Inverse discrete wavelet transform is the inverse process of wavelet decomposition; by upsampling and convolution operations on the wavelet coefficients at each scale, the frequency domain coefficients are restored to the spatial domain image signal.

[0034] In some embodiments, the corresponding inverse discrete wavelet transform (IDWT) reconstructs the original signal by upsampling and convolution of each scale coefficient, and its expression is: in, LL 1(2 m - i ,2 n - j ) represents the coefficient value of the first-layer low-frequency approximation component at coordinates (2m-i, 2n-j). LH 1(2 m - i ,2 n - j ) represents the coefficient value of the first-level horizontal detail component at coordinates (2m-i, 2n-j). HL 1(2 m - i ,2 n - j ) represents the coefficient value of the first layer vertical detail component at coordinates (2m-i, 2n-j). HH 1(2 m - i ,2 n - j ) represents the coefficient value of the first-layer diagonal detail component at coordinates (2m-i, 2n-j). This is the satellite image signal reconstructed by inverse discrete wavelet transform.

[0035] The db4 wavelet was chosen as the basis function because it has good orthogonality, tight support and vanishing moment characteristics, which can effectively separate noise while preserving image edge information, making it particularly suitable for multi-scale analysis of meteorological satellite reflectivity data; while the 3-layer decomposition layer is the optimal balance between spatial detail preservation and computational efficiency. Experiments show that this layer can effectively separate low-frequency information of atmospheric background and high-frequency detail information such as dust and cloud systems.

[0036] In some embodiments, to further improve the temporal stability and noise resistance of FY-4B data, a dual-mode wavelet time enhancement module is designed to perform multi-scale fusion enhancement on FY-4B data one hour before and after the current time. The standard enhancement mode uses a Gaussian weighted average strategy to fuse wavelet coefficients at each scale. First, the Gaussian weights for each time point within the time window are calculated: in t 0 represents the current moment. T For all moments within the time window, exp is the natural exponential function. t The time variable is used. Then, a weighted average is performed on the low-frequency approximation coefficients and high-frequency detail coefficients at each scale: in C k,t For the firstt At the [time]th moment k The wavelet coefficients at the scale are ultimately reconstructed using inverse wavelet transform to obtain the standard enhanced FY-4B data. For regions with low reliability, this embodiment provides a strong filtering mode that uses median filtering instead of mean filtering for high-frequency detail coefficients. Median filtering can effectively suppress impulse noise and outliers, prevent high-frequency noise in untrusted regions from being amplified, and significantly improve the robustness of the fusion results.

[0037] S40: Based on a multi-dimensional dynamic weight calculation system, calculate the comprehensive score of FY-4B data, FY-3D data, FY-3F data and FY-3E data respectively.

[0038] In this embodiment, the multi-dimensional dynamic weight calculation system is a framework for comprehensively evaluating satellite data quality from multiple perspectives, including five evaluation indicators: hourly time stability indicator, reflecting the degree of data fluctuation within an hourly time window (lower fluctuation indicates higher stability); daily time stability indicator, reflecting the degree of data fluctuation within a daily time window; spatial information complexity indicator, reflecting the richness of local texture details in the image (higher complexity indicates richer spatial detail information carried by the image); reference matching degree indicator, measuring the consistency between polar-orbiting satellite data and geostationary satellite data in the time series; and abrupt change penalty coefficient, used to suppress fusion errors caused by data mutations. The comprehensive score is a weighted score obtained based on the above indicators, used to characterize the reliability and contribution of each satellite data in the fusion process.

[0039] In some embodiments, this multi-dimensional dynamic weight calculation system comprehensively considers four aspects: data temporal stability, spatial information complexity, reference matching degree, and abrupt anomalies, assigning adaptive fusion weights to each satellite data set. The temporal stability index is measured by the reciprocal of the coefficient of variation, reflecting the degree of data fluctuation over time. in σ ( X ) is a data sequence X standard deviation μ ( X The mean is denoted as . In this embodiment, hourly and daily time stability were calculated respectively to capture the fluctuation characteristics of short and long time scales. Specifically, the time stability index was applied to the data sequences within hourly and daily time windows respectively to obtain the hourly and daily time stability indices.

[0040] Spatial information complexity is measured by the normalized variance within a 3x3 local window, reflecting the texture detail content of the image. Its calculation formula is as follows: in V ( X ()( x , y Let be a uniform filtering function, representing the local variance at pixel (x,y) after applying the uniform filtering function to the image X. V For local variance images, Ss ( X ()( x , y ) is an index of the spatial information complexity of image X at pixel (x,y). x , y ) represents the pixel coordinates in the image, and min is the function to find the minimum value. i , j () represents the pixel coordinate offset within the local window. V ( X ()( i , j The value of the local variance image at offset coordinate (i,j) is given. It is a local window region, that is, the neighborhood range centered at (x,y).

[0041] The reference matching index is used to calculate the average consistency between the FY-3 series data and the FY-4B reference data within a time window: in, M ( X , Y ) represents the reference matching degree index, and max is the maximum value function. X The satellite data sequence to be evaluated. Y This is a reference data sequence.

[0042] This reference matching index can effectively identify regions in the FY-3 series data that are inconsistent with the FY-4B time series. To suppress sudden outliers in the FY-3 data, this embodiment introduces a sudden change penalty coefficient, using the sigmoid function to perform a non-linear mapping of relative differences. The formula for calculating the sudden change penalty coefficient is as follows: in α To determine the severity of the penalty, this embodiment takes... α =2.0, β To determine the detection sensitivity, this embodiment takes... β =0.15, d ( X ,t 0) represents the current time. t The relative difference between the data value of 0 and the mean of the time window is considered, along with the gradient difference with adjacent time points. d ( X , t The formula for calculating 0) is as follows: in, For the current moment t The data value of 0, The data value from the previous moment. This represents the data value at the next moment.

[0043] The value range of the sudden change penalty coefficient is [0,1]. When the data undergoes a significant sudden change, the penalty coefficient drops rapidly, thereby reducing the fusion weight of abnormal data.

[0044] Based on the above four dimensions of indicators, the comprehensive scores of the three satellite data are calculated respectively. The comprehensive score of FY-4B is obtained by weighting 40% hourly stability, 30% daily stability, and 30% spatial information richness; the comprehensive scores of FY-3D, FY-3F, and FY-3E are obtained by weighting 30% daily stability, 30% spatial information richness, and 40% reference matching degree, and are multiplied by a sudden change penalty coefficient for correction.

[0045] S50: Calculate the initial fusion weight of each satellite data based on the comprehensive score of each satellite data.

[0046] In this embodiment, the initial fusion weight is a weight value obtained by normalizing the comprehensive score of each satellite data, and is used to measure the proportion of each satellite data in the final fusion result. The higher the comprehensive score of the satellite data, the greater its initial fusion weight, indicating that the data contributes more to the fusion.

[0047] In some embodiments, the initial fusion weights are calculated using the following formula. w i : in, s i For the first i The overall score of each satellite s 4b The overall score for the FY-4B satellite. s 3d The overall score for the FY-3D satellite. s 3f This is the overall score for the FY-3F satellite.

[0048] S60: Based on the intelligent reliability switching mechanism, the initial fusion weights are corrected to determine the final fusion weights of each satellite data.

[0049] In this embodiment, the intelligent reliability switching mechanism is a protection strategy that dynamically adjusts the fusion weights based on the matching degree between polar-orbiting satellite data and geostationary satellite data. When there is a significant deviation between polar-orbiting satellite data and geostationary satellite data in certain areas, it indicates that the polar-orbiting satellite data in that area may be unreliable due to factors such as cloud movement, calibration deviation, or transit boundaries. In this case, it is necessary to increase the weight ratio of geostationary satellite data to ensure the reliability of the fusion result. The final fusion weight is the weight value ultimately used by each satellite data in the actual fusion process, determined after reliability assessment and weight correction.

[0050] In some embodiments, to avoid distortion of the fusion results caused by excessive differences between FY-3 series satellite data and FY-4B data, an intelligent reliability switching mechanism is designed based on the current approach. First, the maximum relative difference between FY-3 series data and FY-4B data is calculated: in, FY 3 D , FY 4 B , FY 3 F , FY 3 E These represent the numerical values ​​of their respective satellite data in the corresponding regions.

[0051] when D max When the value is ≤0.25, the FY-3 series satellite data in the corresponding region is deemed reliable. The initial fusion weight remains unchanged, and the initial fusion weight is used as the final fusion weight for that region.

[0052] when D max When the value is greater than 0.25, the FY-3 series satellite data in that region is considered unreliable, and the corresponding FY-4B data weight is forcibly increased. The remaining weight is dynamically allocated to available FY-3 series satellite data. Simultaneously, in unreliable regions, the system automatically switches to FY-4B data in a strong filtering enhancement mode to further suppress the impact of noise and outliers. This mechanism effectively addresses the cloud movement mismatch caused by the difference between FY-3 data and FY-4B observation time, as well as potential calibration biases and overflight boundary issues within the FY-3 data.

[0053] S70: Based on the final fusion weights, the enhanced FY-4B image signal, FY-3D image signal, FY-3F image signal, and FY-3E image signal are weighted and summed to generate a high spatiotemporal resolution surface reflectance product.

[0054] In this embodiment, the surface reflectance product is the final output fusion result. It combines the high temporal resolution advantage of the FY-4B geostationary satellite and the high spatial resolution advantage of the FY-3 series polar-orbiting satellites, possessing a high-frequency update capability of once every 15 minutes and a fine spatial detail representation capability. The weighted summation can be achieved by multiplying the image signals of each satellite at the same pixel location by their respective final fusion weights and then summing them to obtain the fused reflectance value at that pixel location.

[0055] In some embodiments, the wavelet transform-based data fusion result is obtained by weighting it using the dynamic weights described above: in, W 4b The final fusion weights for FY-4B data. FY 4 B enhanced For the enhanced FY-4B image signal, W 3d , W 3f , W 3e These are the final fusion weights for the FY-3D, FY-3F, and FY-3E data, respectively. FY 3 D aligned , FY 3 F aligned , FY 3 E aligned These are the FY-3D, FY-3F, and FY-3E image signals, respectively, registered with FY-4B data. F It is a product for surface reflectivity.

[0056] For pixels missing from all satellite data, NaN values ​​are used for padding. This embodiment outputs not only the reflectivity data of the fused channels but also dynamic weight maps, confidence mask maps, and maximum difference maps for each satellite's corresponding channels to assess the quality of the fused data. The overall fusion process unfolds using a wavelet fusion paradigm of spatial matching, multi-scale decomposition, coefficient fusion, and reconstruction. Furthermore, the robustness of the model is better explained through a region overlap strategy and a dynamic switching mechanism.

[0057] In summary, this application uses FY-4B data as a unified spatiotemporal base and performs fusion calculations only within spatially overlapping regions. Through db4 wavelet three-level multi-scale decomposition and reconstruction, the observation signals from each satellite are separated into low-frequency approximate components and multiple high-frequency detail components for independent processing, effectively avoiding the problems of spatial detail loss and boundary outlier propagation in traditional fusion methods. Regarding temporal enhancement, standard time enhancement is achieved through Gaussian weighted averaging, and strong filtering enhancement is applied to unreliable regions using median filtering to suppress impulse noise and outliers. In terms of weight allocation, the comprehensive score and initial fusion weight of each satellite data are calculated by calculating hourly and daily time stability indices, spatial information complexity indices, reference matching indices, and abrupt change penalty coefficients, achieving multi-dimensional adaptive optimization of the fusion weights. Regarding credibility determination, the credibility of polar-orbiting satellite data is determined by calculating the maximum relative difference Dmax, and when unreliable data is found, the FY-4B weights are increased and a strong filtering mode is switched, effectively correcting fusion errors caused by time gaps and observation misalignments in polar-orbiting satellites.

[0058] Furthermore, this application achieves spatiotemporal complementary fusion of multi-source heterogeneous satellite data through decoupling processing of high and low frequency components, balancing large-scale background continuity with small-scale detail fidelity, and generating high spatiotemporal resolution surface reflectance products. In practical applications, parameters such as wavelet basis function type, decomposition level, time window length, local variance window size, abrupt change detection sensitivity, and confidence threshold can be flexibly adjusted according to specific needs to construct fusion scenarios with different levels of accuracy and efficiency.

[0059] In summary, the core of the technical solution of this application lies in the lightweight strategy of using FY-4B as a unified spatiotemporal base and performing fusion computation only in spatially overlapping regions, as well as the fusion framework that combines db4 wavelet three-layer multi-scale decomposition and reconstruction, multi-dimensional dynamic weight calculation, dual-mode time enhancement and credibility adaptive switching.

[0060] The above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of this application.

Claims

1. A wavelet transform-based atmospheric data fusion method, characterized in that, The method includes: Acquire multi-source satellite data, which includes geostationary meteorological satellite data and polar-orbiting meteorological satellite data. The geostationary meteorological satellite data includes FY-4B data, and the polar-orbiting meteorological satellite data includes FY-3D data, FY-3F data, and FY-3E data. Using the FY-4B data as a unified spatiotemporal base, the spatial overlap region between the FY-3D, FY-3F, and FY-3E data and the FY-4B data is determined. Two-dimensional discrete wavelet transforms are then performed on the data within the spatial overlap region for each of the FY-4B, FY-3D, FY-3F, and FY-3E data, decomposing each satellite data step-by-step into low-frequency approximate components and multiple high-frequency detail components, thus obtaining the wavelet coefficients at each scale for the FY-4B data, the FY-3D data, the FY-3F data, and the FY-3E data. Time-enhanced processing is performed on the wavelet coefficients of each scale of the FY-4B data to obtain enhanced FY-4B wavelet coefficients at each scale. Inverse discrete wavelet transform is then performed on the enhanced FY-4B wavelet coefficients at each scale to reconstruct the enhanced FY-4B image signal. Inverse discrete wavelet transform is then performed on the wavelet coefficients of each scale of the FY-3D data, the FY-3F data, and the FY-3E data, respectively, to reconstruct the FY-3D image signal, FY-3F image signal, and FY-3E image signal registered with the FY-4B data. Based on a multi-dimensional dynamic weight calculation system, the comprehensive scores of the FY-4B data, the FY-3D data, the FY-3F data, and the FY-3E data are calculated respectively. Based on the comprehensive score of each satellite data, the initial fusion weight of each satellite data is calculated; The initial fusion weights are corrected according to the intelligent reliability switching mechanism to determine the final fusion weights of each satellite data. Based on the final fusion weights, the enhanced FY-4B image signal, the FY-3D image signal, the FY-3F image signal, and the FY-3E image signal are weighted and summed to generate a high spatiotemporal resolution surface reflectance product.

2. The atmospheric data fusion method based on wavelet transform as described in claim 1, characterized in that, Two-dimensional discrete wavelet transforms were performed on the FY-4B, FY-3D, FY-3F, and FY-3E data within the spatially overlapping region. This decomposed each satellite data into low-frequency approximate components and multiple high-frequency detail components, yielding wavelet coefficients for each satellite data at various scales, including: Use separate low-pass filters h and high-pass filter g Two-dimensional image signals were processed in both row and column directions. f ( x , y Convolution and downsampling are performed to decompose the two-dimensional image signal into a low-frequency approximate component. LL 1 and three high-frequency detail components, the three high-frequency detail components including a horizontal detail component. LH 1. Vertical detail components HL 1 and diagonal detail components HH 1. The calculation process is expressed as follows: in, LL 1( m , n ) represents the coefficient value of the low-frequency approximate component at position (m,n) after the first decomposition. LH 1( m , n Let be the coefficient value of the horizontal detail component at position (m,n) after the first decomposition. HL 1( m , n ) represents the coefficient value of the vertical detail component at position (m,n) after the first decomposition. HH 1( m , n ) represents the coefficient value of the diagonal detail component at position (m,n) after the first decomposition. i This is the row index variable for the filter coefficients. j For the column index variable of the filter coefficients, h ( i Let be the coefficient value of the low-pass filter at the i-th position. h ( j Let be the coefficient value of the low-pass filter at the j-th position. g ( i Let be the coefficient value of the high-pass filter at the i-th position. g ( j Let be the coefficient value of the high-pass filter at the j-th position. f (2 m - i ,2 n - j ) represents the pixel value of the original two-dimensional image at coordinates (2m-i, 2n-j); A recursive decomposition method is used to approximate the low-frequency components at each level. LL k Continue the decomposition process until the preset number of decomposition layers is reached, to obtain the final low-frequency approximate components and all high-frequency detail components of each layer. Use the final low-frequency approximate components and all high-frequency detail components of each layer as wavelet coefficients of the corresponding satellite data at each scale.

3. The atmospheric data fusion method based on wavelet transform as described in claim 2, characterized in that, The formula for calculating the inverse discrete wavelet transform is as follows: in, LL 1(2 m - i ,2 n - j ) represents the coefficient value of the first-layer low-frequency approximation component at coordinates (2m-i, 2n-j). LH 1(2 m - i ,2 n - j ) represents the coefficient value of the first-level horizontal detail component at coordinates (2m-i, 2n-j). HL 1(2 m - i ,2 n - j ) represents the coefficient value of the first layer vertical detail component at coordinates (2m-i, 2n-j). HH 1(2 m - i ,2 n - j ) represents the coefficient value of the first-layer diagonal detail component at coordinates (2m-i, 2n-j). This is the satellite image signal reconstructed by inverse discrete wavelet transform.

4. The atmospheric data fusion method based on wavelet transform as described in claim 1, characterized in that, Time-enhanced processing was performed on the wavelet coefficients of FY-4B data at each scale to obtain enhanced FY-4B wavelet coefficients at each scale, including: Calculate the Gaussian weights at each time point within the time window. w t The calculation formula is: in, t 0 represents the current moment. T For all moments within the time window, exp is the natural exponential function. t It is a time variable; Based on the Gaussian weights, the wavelet coefficients at each scale are weighted and averaged using the following formula to obtain the standard-enhanced FY-4B wavelet coefficients at each scale. : in, C k,t For the first t At the [time]th moment k Wavelet coefficients at scale; For regions deemed untrustworthy by the aforementioned intelligent trust switching mechanism, median filtering is applied to the high-frequency detail coefficients to obtain strongly filtered and enhanced FY-4B wavelet coefficients at each scale. The formula for calculating the median filter is as follows: in, median This is the median function, which sorts the wavelet coefficients at all times within the time window and takes the median value.

5. The atmospheric data fusion method based on wavelet transform as described in claim 1, characterized in that, Based on a multi-dimensional dynamic weight calculation system, the comprehensive scores of the FY-4B data, FY-3D data, FY-3F data, and FY-3E data are calculated respectively, including: Calculate hourly time stability index, daily time stability index, spatial information complexity index, reference matching index, and sudden change penalty coefficient; The comprehensive score of each satellite data is calculated based on the hourly time stability index, daily time stability index, spatial information complexity index, reference matching index, and sudden change penalty coefficient. The overall score of the FY-4B data is obtained by weighting the hourly time stability index, the daily time stability index, and the spatial information complexity index according to the first preset weight, the second preset weight, and the third preset weight, respectively, and the sum of the first preset weight, the second preset weight, and the third preset weight is 100%. The comprehensive score of any satellite data among the FY-3D, FY-3F, and FY-3E data is obtained by multiplying the daily time stability index, the spatial information complexity index, and the reference matching index by the fourth, fifth, and sixth preset weights, respectively, and then by the sudden change penalty coefficient of the corresponding satellite data. The sum of the fourth, fifth, and sixth preset weights is 100%.

6. The atmospheric data fusion method based on wavelet transform as described in claim 5, characterized in that, The calculation of hourly time stability indices, daily time stability indices, spatial information complexity indices, reference matching indices, and sudden change penalty coefficients includes: Time stability index S t ( X The formula for calculating ) is: in, σ ( X () is a data sequence X standard deviation μ ( X () is a data sequence X The mean of the time stability index is obtained by applying the time stability index to the data sequence within the hourly time window and the data sequence within the daily time window, respectively. The formula for calculating the spatial information complexity index is as follows: in, V ( X ()( x , y Let be a uniform filtering function, representing the local variance at pixel (x,y) after applying the uniform filtering function to the image X. V For local variance images, Ss ( X ()( x , y ) is an index of the spatial information complexity of image X at pixel (x,y). x , y ) represents the pixel coordinates in the image, and min is the function to find the minimum value. i , j () represents the pixel coordinate offset within the local window. V ( X ()( i , j The value of the local variance image at offset coordinate (i,j) is given. This is a local window region, that is, a neighborhood range centered at (x,y); The formula for calculating the reference matching degree index is as follows: in, M ( X , Y ) is the reference matching index, and max is the maximum value function. X The satellite data sequence to be evaluated. Y For reference data sequence; The formula for calculating the sudden change penalty coefficient is as follows: in, α To determine the severity of the punishment, β To detect sensitivity, d ( X , t 0) represents the current time. t The relative difference between a data value of 0 and the mean of the time window.

7. The atmospheric data fusion method based on wavelet transform as described in claim 6, characterized in that, The formula for calculating the relative difference is: in, For the current moment t The data value of 0, The data value from the previous moment. This represents the data value at the next moment.

8. The atmospheric data fusion method based on wavelet transform as described in claim 1, characterized in that, Based on the comprehensive score of each satellite data, the initial fusion weight of each satellite data is calculated, including: The initial fusion weights are calculated using the following formula. w i : in, s i For the first i The overall score of each satellite s 4b The overall score for the FY-4B satellite. s 3d The overall score for the FY-3D satellite. s 3f This is the overall score for the FY-3F satellite.

9. The atmospheric data fusion method based on wavelet transform as described in claim 4, characterized in that, According to the intelligent trustworthiness switching mechanism, the initial fusion weights are corrected, including: The maximum relative difference between the FY-3 series data and the FY-4B data is calculated using the following formula. D max : in, FY 3 D , FY 4 B , FY 3 F , FY 3 E These are the numerical values ​​of their respective satellite data in the corresponding regions; when D max When the value is ≤0.25, the FY-3 series satellite data in the corresponding region is determined to be reliable. The initial fusion weight remains unchanged, and the initial fusion weight is used as the final fusion weight for that region. when D max When the value is greater than 0.25, the FY-3 series satellite data in the corresponding area is determined to be unreliable. The fusion weight of the FY-4B data is forcibly increased, and the area is switched to a strong filtering enhancement mode. The increased fusion weight of the FY-4B data is a preset weight threshold. The fusion weights of the FY-3D, FY-3F, and FY-3E data in the area are redistributed according to the proportion of their initial fusion weights to obtain the final fusion weight of the area.

10. The atmospheric data fusion method based on wavelet transform as described in claim 1, characterized in that, High spatiotemporal resolution surface reflectance products are generated using the following formula. F : in, W 4b The final fusion weights for FY-4B data. FY 4 B enhanced For the enhanced FY-4B image signal, W 3d , W 3f , W 3e These are the final fusion weights for the FY-3D, FY-3F, and FY-3E data, respectively. FY 3 D aligned , FY 3 F aligned , FY 3 E aligned These are the FY-3D, FY-3F, and FY-3E image signals, respectively, registered with FY-4B data.

Images