Multi-source data fusion method based on adaptive weighted spherical harmonic signals, atmospheric quality analysis method and remote sensing vegetation monitoring method

Abstract

The present application relates to a kind of multi-source data fusion method based on adaptive weighted spherical harmonic signal, atmospheric quality analysis method and remote sensing vegetation monitoring method, to realize the fusion of multi-source data with spherical distribution characteristics, comprising: obtaining the multi-source spatio-temporal observation data to be fused;Based on the preset time window length, the residual time series with the center of current time is calculated, the noise level of data source dimension is estimated, and the covariance matrix based on noise intensity and spatial position is constructed;Based on the covariance matrix, the comprehensive method equation is constructed, and the preliminary spherical harmonic coefficient fusion solution is solved;Based on the current spherical harmonic coefficient fusion solution, the variance component of data source dimension is calculated, the fusion method equation of current time is constructed, the spherical harmonic coefficient fusion solution is updated, and the method equation weight below the preset threshold is screened out and the step is repeated;The final spherical harmonic coefficient fusion solution is obtained, and multi-source data fusion is realized.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a multi-source data fusion method, an atmospheric quality analysis method, and a remote sensing vegetation monitoring method based on adaptive weighted spherical harmonic signals. Background Technology

[0002] Spherical harmonic functions (SHFCs) are a complete orthogonal function system on a sphere, commonly used to represent and analyze spatiotemporal data with spherical distribution characteristics. They are widely used in modeling and compression of the global atmosphere, ocean, and geophysical fields, as well as in remote sensing data processing and remote sensing image downscaling models. With the development of observation methods and model simulations, more and more applications require the simultaneous fusion of spatiotemporal data from different sources to extract more stable and reliable SHFC signals. However, most existing extraction methods are designed for single data sources, directly establishing observation equations and solving for coefficients by combining gridded or discretized data with SHFC basis functions, making it difficult to fully utilize the complementarity between multi-source data. Even in some fusion attempts, equal-weighted least squares or fixed-weighting methods are typically used, assuming that each data source has the same or a priori set noise level.

[0003] In existing solutions, Chen H et al.'s paper "Multi-Source Soil Moisture Data Fusion Based on Spherical Cap Harmonic Analysis and Helmert Variance Component Estimation in the Western US" uses a fusion method that does not consider the weight matrix. However, this assumption does not match reality: whether it is soil moisture retrieval or optical remote sensing products such as the normalized vegetation index, their errors are highly sensitive to underlying surface type, vegetation cover, soil texture and moisture state, topographic relief, solar-sensor observation geometry, cloud / thin cloud and aerosol conditions, etc., showing significant spatial non-stationarity and heteroscedasticity. Taking the normalized vegetation index as an example, the product usually provides quality indicators at the pixel level, clearly distinguishing different situations such as high quality, edge quality, and influence by clouds or strong aerosols; these indicators directly indicate that the error variance of adjacent pixels at the same time is not equal.

[0004] Furthermore, existing methods lack adaptive mechanisms based on observation residuals or time series statistics, and fail to effectively adjust weights according to the characteristics of the data itself, thus limiting the accuracy and robustness of spherical harmonic signal extraction in the context of multi-source data fusion.

[0005] Therefore, there is an urgent need for a multi-source data fusion method to improve the reliability and applicability of fusion estimation. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of the existing technology by providing a multi-source data fusion method, an atmospheric quality analysis method, and a remote sensing vegetation monitoring method based on adaptive weighted spherical harmonic signals. This aims to solve or partially solve the problem that it is difficult to characterize the spatiotemporal differences of data noise in the fusion of multi-source data with spherical distribution characteristics, resulting in insufficient accuracy and robustness of the fusion results.

[0007] The objective of this invention can be achieved through the following technical solutions:

[0008] One aspect of the present invention provides a multi-source data fusion method based on adaptive weighted spherical harmonic signals, used to achieve the fusion of multi-source data with spherical distribution characteristics, the method comprising the following steps:

[0009] Acquire multi-source spatiotemporal observation data to be fused;

[0010] Based on a preset time window length, the residual time series centered on the current moment is calculated, the noise level of the data source dimension is estimated, and a covariance matrix based on noise intensity and spatial location is constructed.

[0011] Based on the covariance matrix, a comprehensive normal equation is constructed, and a preliminary spherical harmonic coefficient fusion solution is obtained.

[0012] Based on the current spherical harmonic coefficient fusion solution, calculate the variance component of the data source dimension, construct the fusion normal equation at the current moment, update the spherical harmonic coefficient fusion solution, filter out the normal equation weights below the preset threshold and repeat this step until the preset conditions are met.

[0013] The final spherical harmonic coefficient fusion solution is obtained, realizing multi-source data fusion.

[0014] As a preferred technical solution, the process of calculating the residual time series centered on the current moment based on a preset time window length, estimating the noise level of the data source dimension, and constructing a covariance matrix based on noise intensity and spatial location includes the following steps:

[0015] Establish observation equations for the aforementioned multi-source spatiotemporal observation data, observation noise vector, and spherical harmonic coefficient vector;

[0016] For each epoch, calculate the corresponding spherical harmonic coefficients and the residuals as an estimate of the observation noise;

[0017] At multiple pre-selected epochs, the root mean square value of the residual is calculated for each observation location to obtain the noise level vector of each location in the spatiotemporal observation data corresponding to each data source.

[0018] Based on the noise level vector and the spatial distance between locations, a covariance matrix is ​​constructed.

[0019] As a preferred technical solution, the covariance matrix is:

[0020]

[0021]

[0022]

[0023]

[0024] in, For the first The covariance matrix of observation data from each data source. , These represent the location-dependent noise intensity and the spatial distance matrix, respectively. For kernel function, Indicates the first in the observation data With the Spatial distance between locations Represents a diagonal matrix. For the first Observation location in observation data from multiple data sources The noise level vector, The total number of observation locations, Number of epochs For the first The residual of each epoch, It is a set of epochs.

[0025] As a preferred technical solution, the combined normal equation and the preliminary spherical harmonic coefficient fusion solution obtained by solving it are as follows:

[0026]

[0027]

[0028]

[0029] in, Total number of data sources Indicates transpose. For the weight matrix, This is a preliminary solution for the fusion of spherical harmonic coefficients. , For the simplified terms defined, For the first The design vector of each data source is composed of spherical harmonic basis functions. is the spherical harmonic coefficient.

[0030] As a preferred technical solution, the process of updating the spherical harmonic coefficient fusion solution includes the following steps:

[0031] Based on the preliminary spherical harmonic coefficient fusion solution, the residual of the data source dimension is calculated;

[0032] Variance components of the data source dimension are calculated based on residuals;

[0033] Taking the reciprocal of the square difference component yields the weights of the normal equation;

[0034] The fusion normal equation for the current moment is constructed based on the weights of the normal equation, and the new fusion solution of the spherical harmonic coefficients is obtained by solving it.

[0035] As a preferred technical solution, the residuals and variances of the data source dimension are:

[0036]

[0037]

[0038]

[0039] in, , The first Residuals and variances of each data source For spatiotemporal observation data, The design matrix is ​​composed of spherical harmonic basis functions. For spherical harmonic coefficients, Indicates transpose. For the weight matrix, The number of observation locations in the observation data. For the number of data sources, For the defined simplification terms For the first The design vector of each data source is composed of spherical harmonic basis functions. Let be the trace of the matrix.

[0040] As the preferred technical solution, the fusion equations at the current moment and the new spherical harmonic coefficient fusion solution are:

[0041]

[0042]

[0043]

[0044] in, For the number of data sources, For the first Normal equation weights for each data source The design matrix is ​​composed of spherical harmonic basis functions. Indicates transpose. For the weight matrix, For spatiotemporal observation data, To obtain a new spherical harmonic coefficient fusion solution.

[0045] As a preferred technical solution, the preset condition is that the difference in unit weight variance between different data sources calculated based on the current spherical harmonic coefficient fusion solution is less than a preset value.

[0046] Another aspect of the present invention provides an atmospheric quality analysis method based on adaptive weighting of multi-source data and extraction of spherical harmonic signals, comprising the following steps:

[0047] Acquire multiple atmospheric reanalysis datasets;

[0048] Construct an atmospheric load model based on surface integral and vertical integral, and considering tidal rejection and reverse pressure correction;

[0049] Based on the atmospheric reanalysis dataset, the corresponding atmospheric load is calculated using the atmospheric load model. By removing the average within the observation period, the observations corresponding to each reanalysis dataset are obtained as multi-source spatiotemporal observation data to be fused.

[0050] The aforementioned multi-source data fusion method based on adaptive weighted spherical harmonic signals is used to fuse the multi-source spatiotemporal observation data to be fused, and the final spherical harmonic coefficient fusion solution is obtained as the best estimate of the non-tidal high-frequency atmospheric quality change in the spherical harmonic domain after fusion. Atmospheric quality analysis is realized based on the best estimate.

[0051] Another aspect of the present invention provides a remote sensing vegetation monitoring method based on multi-source data adaptive weighting and spherical harmonic signal extraction, comprising the following steps:

[0052] Acquire multiple remote sensing imaging datasets;

[0053] Preprocessing is performed on the remote sensing imaging dataset, including resolution unification and correction.

[0054] The preprocessed remote sensing imaging dataset is used as the multi-source spatiotemporal observation data to be fused. The aforementioned multi-source data fusion method based on adaptive weighted spherical harmonic signals is used to fuse the multi-source spatiotemporal observation data to be fused, and the final spherical harmonic coefficient fusion solution is obtained as the fused remote sensing image.

[0055] Vegetation monitoring is achieved based on fused remote sensing images.

[0056] Compared with the prior art, the present invention has at least one of the following beneficial effects:

[0057] (1) Realizing Multi-Source Data Fusion Based on Two-Layer Adaptive Weighting: This invention achieves a two-layer adaptive weighting mechanism by characterizing noise at both spatial and temporal levels. Specifically, at the spatial level, this invention estimates the noise level of the data source dimension, constructs a covariance matrix based on noise intensity and spatial location, and uses the root mean square value of the residuals to construct a location-related diagonal weight matrix, thereby suppressing the interference of high-noise regions on the results. At the temporal level, the overall noise level of the data source is estimated at a single epoch, and the weighting coefficients of the normal equation are dynamically adjusted accordingly. Compared with traditional equal-weight or fixed-weight matrix methods, this method can significantly improve the accuracy and robustness of spherical harmonic signal extraction in the multi-source data fusion process, and has high practical value.

[0058] (2) Strong adaptability in atmospheric quality analysis: This invention fuses observations obtained from multiple atmospheric reanalysis datasets. When analyzing the Earth's gravity field by satellite inversion, the fused solution can adaptively approach the dataset that performs best at different time periods.

[0059] (3) Strong applicability in remote sensing vegetation monitoring: This invention fuses multiple preprocessed remote sensing imaging datasets and prepares monitoring results. The fused solution obtained under complex terrain can clearly show the spatial distribution pattern. Attached Figure Description

[0060] Figure 1 This is a flowchart of the multi-source data fusion method based on adaptive weighted spherical harmonic signals in the embodiment;

[0061] Figure 2 The flowchart shows the atmospheric quality analysis method based on multi-source data adaptive weighting and spherical harmonic signal extraction in the embodiment.

[0062] Figure 3 This is a schematic diagram of the second-order normalized fusion weight time series from June to July 2020 in the example embodiment;

[0063] Figure 4 The correlation coefficients between some low-order terms of the ERA5 standalone solution, the CRA-40 standalone solution, and the fused solution obtained in this invention, and the reference atmospheric demixing product are shown in the examples.

[0064] Figure 5 This is a flowchart of the remote sensing vegetation monitoring method based on multi-source data adaptive weighting and spherical harmonic signal extraction in the embodiment.

[0065] Figure 6This is a schematic diagram illustrating the spatial distribution of the differences in MODIS, VIIRS, and fused solutions in the northern Andes Mountains in the example. Detailed Implementation

[0066] 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, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0067] Example 1

[0068] To address the problem that existing methods, which generally employ equal or fixed weight matrices, struggle to characterize the spatiotemporal differences in data noise, leading to insufficient accuracy and robustness in the fusion results, this embodiment provides a multi-source data fusion method based on adaptive weighted spherical harmonic signals for fusing multi-source data with spherical distribution characteristics. This method introduces dynamic weight matrix construction and a time-varying weighted fusion mechanism within the framework of spherical harmonic expansion, effectively utilizing the complementarity between multi-source data, thereby significantly improving the accuracy and stability of the estimation.

[0069] See Figure 1 The method includes the following steps:

[0070] Step S1: Obtain the multi-source spatiotemporal observation data to be fused.

[0071] Step S2: Based on the preset time window length, calculate the residual time series centered on the current moment, estimate the noise level of the data source dimension, and construct a covariance matrix based on noise intensity and spatial location.

[0072] At some point It has multi-source observation data. ,in , These represent the number of data sources and the number of observation locations in the observation data, respectively. The parameter to be estimated is the spherical harmonic coefficient vector. The two satisfy the observation equation:

[0073]

[0074] in The design matrix is ​​composed of spherical harmonic basis functions. for The observation noise vector at time step. This observation model provides the basic framework for inverting spherical harmonic coefficients from discrete grid data.

[0075] This step utilizes... For time series residuals of a specified length centered on a given time point, statistical estimates of the noise level are performed for each data source. Specifically, for data sources centered on a given time point... A specific moment (i.e., an epoch) within a time series of a specified length, centered at a given time point. First use the unit weight matrix Independently calculate spherical harmonic coefficients :

[0076]

[0077] And calculate the residuals as an estimate of the observation noise:

[0078]

[0079] At each observation location in an epoch of sequence length K. Calculate the root mean square value of the residuals:

[0080]

[0081] Therefore, we obtain the first... Noise level vector of group observation data And based on this, construct the covariance matrix of the data source:

[0082]

[0083] in and They represent the first and second elements in the observation vector, respectively. With the One location, This indicates the spatial distance between these two locations (such as the distance between a great circle on a sphere). This is the kernel function. The weight matrix used for weighted least squares is:

[0084]

[0085] This construction simultaneously encodes the location-dependent noise intensity (by... (Given) and spatially related structures (by) (Given), thus transforming the stable spatial noise patterns in the time series into location-related weights, which can effectively suppress the influence of high-noise regions on the results and achieve adaptive weighting in the spatial domain.

[0086] It should be noted that the kernel function in the above formula... The choice can be made according to the needs: when the unit kernel function is selected, the diagonal weights are obtained; if spatial correlation needs to be considered, the Gaussian kernel or the exponential decay kernel can be selected. The Gaussian kernel is often used in the case where the noise gradually decreases with distance due to its smooth decay characteristics. For example, atmospheric pressure fields and terrestrial water storage fields usually have smooth and continuous spatial distributions, which are suitable for using the Gaussian kernel to characterize local correlation. The exponential decay kernel is more commonly used in the case of rapid decay or the existence of obvious correlation length thresholds. For example, the correlation of observations such as seismic wave residuals and radar echoes will disappear rapidly after exceeding a certain spatial scale. Such data are more suitable for modeling with the exponential decay kernel.

[0087] The purpose of this step is to extract stable spatial noise patterns from the time series and transform them into location-dependent weight matrices. This effectively suppresses the influence of high-noise regions on the results, achieving adaptive weighting in the spatial domain, rather than simple equal weighting.

[0088] It should be noted that the noise level of the data is not constant but evolves gradually over time. Especially in the early stages, due to relatively outdated observation methods and computational models, the overall noise level is usually significantly higher than that of data from recent years. To address this, this method supports setting the update cycle of the weight matrix as needed, such as re-estimating it every 3 months, 6 months, or 12 months, thereby achieving a mechanism for dynamically updating the weight matrix and better reflecting the true noise level of the data at different time periods.

[0089] Step S3: Construct a comprehensive normal equation based on the covariance matrix and solve for the preliminary spherical harmonic coefficient fusion solution.

[0090] At the present moment Assuming that the data sets are independent of each other, the comprehensive method equation can be obtained:

[0091]

[0092] At this moment, note:

[0093]

[0094] Then we have a preliminary solution for the fusion of spherical harmonic coefficients:

[0095]

[0096] Step S4: Based on the current spherical harmonic coefficient fusion solution, calculate the variance component of the data source dimension, construct the fusion normal equation at the current moment, update the spherical harmonic coefficient fusion solution, filter out the normal equation weights below the preset threshold, and repeat this step until the preset conditions are met.

[0097] The fusion solution obtained in step S3 may be due to the weighting of multiple data. This is unreasonable; therefore, this step utilizes the variance component estimation method to estimate the unit weighted variance of each group of data. Make an estimate and adjust the weights accordingly.

[0098] Based on the initial fusion solution, calculate the residuals for each group of data:

[0099]

[0100] And use this to estimate the variance components of each data point at the current moment:

[0101]

[0102] Let be the trace of the matrix. At this point, It is a scalar quantity used to characterize the overall noise intensity of the data source at a given moment. Its purpose is to dynamically reflect the changes in observation quality over time, avoiding the masking of short-term anomalies by long-term averaging.

[0103] The final spherical harmonic coefficients are calculated by weighted fusion of normal equations from different data sources. Specifically, the reciprocal of the noise variance of each data source is used as the weight of the normal equations.

[0104]

[0105] Finally, using the weights of the time-varying normal equation, the fusion normal equation for the current time step is constructed:

[0106]

[0107] And from this, a new fusion solution is obtained:

[0108]

[0109] At this point, it is still possible to further estimate the unit weight variance of each group of data. To determine whether the weights of different data are reasonable, if not, repeat the above steps until the unit weight variances of several groups of data are approximately equal.

[0110] Building upon this, this method further introduces a quality control mechanism: if the weight of a certain data source at the current location... Significantly lower than the preset threshold If a data source is deemed to have excessively high noise levels and poor performance at that particular moment, it will be automatically removed during the fusion process, thus preventing inferior data sources from interfering with the final result. This mechanism ensures that the fusion process retains the complementarity of multiple sources while further improving the robustness and reliability of the results.

[0111] Step S5 yields the final spherical harmonic coefficient fusion solution, achieving multi-source data fusion.

[0112] This method characterizes noise at both spatial and temporal levels, achieving a two-layer adaptive weighting mechanism. Spatially, it constructs a location-dependent diagonal weight matrix using the root mean square value of the residuals, thereby suppressing the interference of high-noise regions on the results. Temporally, it estimates the overall noise level of the data source at a single epoch and dynamically adjusts the weighting coefficients of the normal equation accordingly. Compared with traditional equal-weight or fixed-weight matrix methods, this method significantly improves the accuracy and robustness of spherical harmonic signal extraction during multi-source data fusion, demonstrating high practical value.

[0113] Example 2

[0114] Based on Example 1, this example provides an atmospheric quality analysis method based on adaptive weighting of multi-source data and extraction of spherical harmonic signals. See [link to example]. Figure 2 It includes the following steps:

[0115] Step S1: Obtain multiple atmospheric reanalysis datasets.

[0116] In this embodiment, two sets of reanalysis data, ERA5 and CRA-40, were selected as datasets from different sources, and atmospheric surface pressure was used as the basis for the datasets. Multi-layer specific humidity Multi-layer temperature Topographical elevation The four variables were unified to a three-hour time interval and a consistent spatial grid, and the surface pressure sequence was subjected to daily and semi-diurnal tidal removal and backpressure (IB) correction.

[0117] Step S2: Construct an atmospheric load model based on surface integral and vertical integral, and considering tidal rejection and reverse pressure correction.

[0118] Construct order-dependent atmospheric loading based on the definition of atmospheric vertical integral Its physical meaning is the non-tidal high-frequency atmospheric mass variation relative to the long-term average, among which , These are co-latitude (i.e., the remainder of latitude) and longitude, respectively. The following is the calculation process for atmospheric load. Atmospheric load is defined as:

[0119]

[0120] in, The average radius of the Earth Represents the distance from any altitude to the Earth's center; Atmospheric density is a property related to spatial location. Earth's center distance and time Related functions, The order of atmospheric load.

[0121] In the setting of geometric and gravity parameters, the Earth ellipsoid radius related to latitude is denoted as... The geoid undulates , the height of terrain The geometric height is denoted as At this point, the distance from any height to the Earth's center is written as:

[0122]

[0123] Both the ERA5 and CRA-40 datasets employ the hydrostatic approximation assumption, namely:

[0124]

[0125] In the formula, For pressure, The acceleration due to gravity at any height is calculated using the following formula:

[0126]

[0127] in, This is the gravitational acceleration value taken from the surface of the GRS80 ellipsoid. At this time, the expression for the atmospheric load is:

[0128]

[0129] In the formula, This refers to the atmospheric surface pressure.

[0130] At this point, the integral variable is pressure. The distance to the center of the earth, gravitational acceleration, and pressure are all functions of geometric height. Therefore, we need to discretize the formula to facilitate calculation.

[0131] Assuming all changes in Earth's atmospheric mass occur within a thin layer of zero thickness that is tightly bound to Earth's topography, the formula becomes:

[0132]

[0133] This formula is the formula for calculating surface integrals.

[0134] Furthermore, considering the vertical structure of the atmosphere, the atmospheric load is discretized using the vertical layer definitions of the atmospheric dataset. ERA5 provides variables in mixed Sigma layers, which are defined by the air pressure at the model boundaries (i.e., half-layers), and uses constants reflecting the vertical coordinates. and To achieve:

[0135]

[0136] In the formula Represents the number of floors. This represents the half-level, which is the interface between two levels. for Atmospheric pressure at the half-layer position.

[0137] For CRA-40, variables are provided at the isobaric surface, and the pressure is directly used in the calculation according to the standard isobaric surface. In meteorology, altitude is usually geopotential height. The potential height difference of each layer can be calculated using the following formula:

[0138] In the formula , where is a virtual temperature, The gas constant for dry air. , These are gas temperature and specific humidity, respectively. The standard gravitational acceleration recommended by the World Meteorological Organization is 9.80665. .

[0139] Therefore, the height of each floor can be calculated using the following formula:

[0140]

[0141] In the formula The maximum number of vertical layers defined by the atmospheric dataset. for The potential height at the half-floor location. The conversion between potential height and geometric height is as follows:

[0142]

[0143] This unifies the multi-level variables of ERA5 and CRA-40 into a geometrical altitude coordinate system. At this point, the expression for atmospheric load is discretized as follows:

[0144]

[0145] This formula is the formula for calculating the perpendicular integral.

[0146] Step S3: Based on the atmospheric reanalysis dataset, calculate the corresponding atmospheric load using the atmospheric load model. By removing the average within the observation period, obtain the observations corresponding to each reanalysis dataset as the multi-source spatiotemporal observation data to be fused.

[0147] A hybrid processing method is adopted:

[0148]

[0149] In the formula and These represent the surface integral and the perpendicular integral, respectively. The surface pressure before tidal rejection and backpressure correction. This indicates the surface pressure after these two treatments. This represents the load Love number. Subsequently... The long-term average from 2003 to 2014 was removed, resulting in Based on this, Transform into:

[0150]

[0151] In the formula, This is the average density of the Earth. This indicates the data source. Two sets can be obtained for ERA5 and CRA-40 respectively. .

[0152] Step S4: Using the multi-source data fusion method based on adaptive weighted spherical harmonic signals from Example 1, the multi-source spatiotemporal observation data to be fused are fused to obtain the final spherical harmonic coefficient fusion solution, which serves as the best estimate of the non-tidal high-frequency atmospheric quality change in the spherical harmonic domain after fusion. Atmospheric quality analysis is then performed based on the best estimate.

[0153] The observations obtained in step S3 Substitute the fusion method from Example 1 into the solution. At each time step... , establish the first The fusion equations of order are derived and variance components are estimated to determine the adaptive weights at corresponding time points. Iteration stops when the variance components of the two sets of data are approximately equal, yielding the spherical harmonic coefficient solution. .

[0154] The solution vector contains all spherical harmonic coefficients from 0 to the maximum order, but in the nth order processing of this embodiment, only the coefficient components corresponding to order n are meaningful output results.

[0155] The final result is the best estimate of the non-tidal high-frequency atmospheric mass change in the spherical harmonic domain after fusion. It can be directly used for demixing in satellite gravity missions, or the equivalent water height change can be obtained through transformation for the analysis of mass redistribution in events such as extreme rainfall.

[0156] See Figure 3 This is a second-order normalized fusion weight time series for the period from June to July 2020, where ERA5 serves as a constant reference line (with a value of 1), and CRA-40 varies with time at each 3-hour epoch. The vertical axis represents the relative value of the fusion weight time series with respect to the constant reference line, and the horizontal axis represents the time index (3-hour step size), with each increment of 1 indicating an increase of 3 hours.

[0157] See Figure 4 The correlation coefficient (vertical axis) between the ERA5-only solution, the CRA40-only solution, and the fusion solution obtained by this method (i.e., the horizontal axis in the figure) and the latest version of the atmospheric demixing product published by the German Geosciences Research Center is shown.

[0158] Table 1 shows the signal-to-noise ratio (SNR) of mass change in each watershed and the root mean square (RMS) of noise, used to characterize noise levels, obtained by applying the individual solutions of ERA5 and CRA-40, the fused solution obtained by this method, and the fused solution obtained by setting only a unit weight matrix without considering the dynamic weight matrix of the dataset, to the Earth's gravity field inversion from 2019 to 2021. It is evident that considering the dynamic weight matrix of the dataset, compared to the fused solution using the unit weight matrix, has a more significant suppression effect on the noise level of most watersheds.

[0159] Table 1. Root mean square noise values ​​of the time-varying gravity field and signal-to-noise ratio of the mass change signal in various watersheds.

[0160]

[0161] The results show that the fused solution can adaptively approach the dataset that performs best at different time periods.

[0162] Example 3

[0163] Based on Example 1, this example provides a remote sensing vegetation monitoring method based on multi-source data adaptive weighting and spherical harmonic signal extraction. See [link to example]. Figure 5 It includes the following steps:

[0164] Step S1: Obtain multiple remote sensing imaging datasets.

[0165] The data used in this embodiment is the Normalized Difference Vegetation Index (NDVI), which is calculated based on the reflectance of red and near-infrared bands. This index effectively reflects the surface vegetation cover and growth status, and is a typical data type for remote sensing vegetation monitoring. The selected data comes from two sources: firstly, the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard NASA's Terra satellite. This product has provided long-term continuous observation sequences since 2000, has a mature algorithm system, and is widely used in ecological and climate monitoring research; however, its spatial resolution is limited, and the sensor has experienced some aging effects after many years of operation. Secondly, the data comes from the Visible Infrared Imaging Radiometer (VIIRS) aboard the US Suomi National Polar Partnership (Suomi NPP) satellite. This product has provided data since 2012 and has advantages such as advanced sensor design, higher spatial resolution, and significant improvements in geometric distortion and observation angle effects; however, its observation time series is relatively short, making it difficult to support long-term trend analysis alone. To ensure the reliability of the fusion results, this embodiment removes low-quality image data affected by clouds, snow, strong aerosols, etc. during the use process through quality control. Thus, while taking into account the advantages of MODIS long time series and VIIRS spatial quality, the fusion of the two types of data is achieved, resulting in a more stable and accurate vegetation index time series product.

[0166] Step S2 involves preprocessing the remote sensing imaging dataset, which includes resolution unification and correction.

[0167] To improve the comparability of the two-source NDVI, both monthly products were uniformly sampled to a resolution of 0.25°×0.25°. A land mask was constructed based on the land-water markers of the quality layer, and aggregation statistics were performed only on effective land pixels. This eliminated the influence of the ocean while retaining physically reasonable negative signals at high latitudes. Subsequently, before fusion, the median difference between the two sets of NDVI data was calculated based on the common historical months (12 months) preceding the target month. Based on this, additive correction was applied to the MODIS with poorer spatial quality to suppress additive systematic bias in the robust estimation of outliers at the pixel scale. After the above quality control and systematic error correction, uncertainty assessment and weighted fusion were then carried out.

[0168] Step S3: The preprocessed remote sensing imaging dataset is used as the multi-source spatiotemporal observation data to be fused. The multi-source data fusion method based on adaptive weighted spherical harmonic signals in Example 1 is used to fuse the multi-source spatiotemporal observation data to be fused, and the final spherical harmonic coefficient fusion solution is obtained as the fused remote sensing image.

[0169] Using NDVI products provided by MODIS and VIIRS as multi-source observation data:

[0170]

[0171] Based on the multi-source data fusion method of Example 1, a normal equation is first constructed using a unit weight matrix to calculate the residual time series and construct a spatial noise pattern, creating a diagonal weight matrix based on observation data from 12 months of 2024. Subsequently, a fusion normal equation is constructed, and iterative estimation using variance components is performed to adjust the weights of different observation data. When the variance components of the two sets of data are:

[0172]

[0173] The iteration terminates when the value is less than the threshold (set to 0.05 in this embodiment).

[0174] Step S4: Vegetation monitoring is achieved based on the fused remote sensing images.

[0175] See Figure 6 The figure illustrates the spatial distribution differences between MODIS, VIIRS, and the fusion solution of this embodiment (i.e., the fusion solution and the unit weight matrix fusion solution in the figure) in the northern Andes Mountains. The gray levels in the figure represent the solution values ​​at the corresponding longitude (Lon) and latitude (Lat). It can be seen that in areas with dramatic topographic relief and complex slope and aspect, the NDVI retrieved by MODIS is significantly affected by factors such as illumination differences, topographic shadows, and anisotropy of the two-way reflectance distribution function, resulting in poor applicability. It is almost impossible to identify the distribution characteristics of low NDVI in the mountains in the figure. In contrast, the fusion solution, compared with the traditional unit weight matrix fusion method, better inherits the advantages of VIIRS in complex terrain and more clearly reveals the spatial distribution pattern of low NDVI in the Andes Mountains.

[0176] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in the present invention, and these modifications or substitutions should all be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A multi-source data fusion method based on adaptive weighted spherical harmonic signals, characterized in that, The application discloses a method for realizing fusion of multi-source data with spherical distribution characteristics. The method comprises the following steps: Obtaining multi-source spatio-temporal observation data to be fused; Based on a preset time window length, calculating a residual time series centered on a current time, estimating noise levels of data source dimensions, and constructing a covariance matrix based on noise intensity and spatial position; Based on the covariance matrix, constructing a comprehensive method equation, and solving a preliminary spherical harmonic coefficient fusion solution; Based on the current spherical harmonic coefficient fusion solution, calculating variance components of the data source dimensions, constructing a fusion method equation at the current time, updating the spherical harmonic coefficient fusion solution, screening out method equation weights below a preset threshold, and repeating the step until a preset condition is met; Obtaining a final spherical harmonic coefficient fusion solution, and realizing multi-source data fusion, wherein is the total number of data sources, denotes the transpose, is the weight matrix, is the covariance matrix of the observation data of the th data source, is the preliminary spherical harmonic coefficient fusion solution, , is the defined reduction term, is the design vector of the th data source composed of spherical harmonic basis functions, is the spherical harmonic coefficient, is the spatio-temporal observation data.

2. The multi-source data fusion method based on adaptive weighted spherical harmonic signals according to claim 1, characterized in that, The comprehensive method equation and the preliminary spherical harmonic coefficient fusion solution obtained by solving are as follows: Based on a preset time window length, calculating a residual time series centered on a current time, estimating noise levels of data source dimensions, and constructing a covariance matrix based on noise intensity and spatial position comprise the following steps: Establishing an observation equation about the multi-source spatio-temporal observation data, an observation noise vector and a spherical harmonic coefficient vector; For each epoch, calculating corresponding spherical harmonic coefficients, and residuals as observation noise estimates; On a plurality of preselected epochs, calculating root mean square values of residuals for each observation position, and obtaining noise level vectors of each position in spatio-temporal observation data corresponding to each data source; 3. The multi-source data fusion method based on adaptive weighted spherical harmonic signals according to claim 1, characterized in that, Based on the noise level vectors and spatial distances between positions, constructing a covariance matrix. wherein, , are the position-dependent noise strength, the spatial distance matrix, respectively, is the kernel function, denotes the spatial distance between the th and the th position in the observation data, denotes the diagonal matrix, is the noise level vector for the observation position in the observation data of the th data source, is the total number of observation positions, is the number of epochs, is the residual for the th epoch, is the set of epochs.

4. The multi-source data fusion method based on adaptive weighted spherical harmonic signals according to claim 1, characterized in that, The covariance matrix is as follows: The process of updating the spherical harmonic coefficient fusion solution comprises the following steps: Based on the preliminary spherical harmonic coefficient fusion solution, calculating residuals of the data source dimensions; Based on the residuals, calculating variance components of the data source dimensions; Taking inverses of the variance components, obtaining method equation weights; 5. The multi-source data fusion method based on adaptive weighted spherical harmonic signals according to claim 4, characterized in that, Based on the method equation weights, constructing a fusion method equation at the current time, and solving a new spherical harmonic coefficient fusion solution. wherein , are the residuals, variances, respectively, of the th data source, is a design matrix composed of spherical harmonic basis functions, is a spherical harmonic coefficient, denotes the transpose, is a weight matrix, is the number of observation positions in the observation data, is the number of data sources, is a defined reduction term, is a design vector composed of spherical harmonic basis functions of the th data source, is the trace of a matrix.

6. The multi-source data fusion method based on adaptive weighted spherical harmonic signals according to claim 4, characterized in that, The residuals and variance of the data source dimensions are as follows: wherein, is the number of data sources, is the normal equation weight of the th data source, is a design matrix composed of spherical harmonic basis functions, denotes the transpose, is a weight matrix, is the spatio-temporal observation data, is the new spherical harmonic coefficient fusion solution obtained.

7. The multi-source data fusion method based on adaptive weighted spherical harmonic signals according to claim 1, characterized in that, The fusion method equation at the current time and the new spherical harmonic coefficient fusion solution are as follows:

8. An atmospheric quality analysis method based on multi-source data adaptive weighting and spherical harmonic signal extraction, characterized in that, The preset condition is that a difference between unit weight variances of different data sources calculated based on the current spherical harmonic coefficient fusion solution is less than a preset value. The method comprises the following steps: Obtaining a plurality of atmospheric reanalysis data sets; Constructing an atmospheric loading model based on surface integration and vertical integration, and considering tide removal and inverse barometric correction; Based on the atmospheric reanalysis data sets, calculating corresponding atmospheric loads by using the atmospheric loading model, removing averages in an observation period, obtaining observation quantities corresponding to each reanalysis data set as multi-source spatio-temporal observation data to be fused; 9.A remote sensing vegetation monitoring method based on multi-source data adaptive weighting and spherical harmonic signal extraction, characterized in that, The multi-source spatio-temporal observation data to be fused are fused by using the multi-source data fusion method based on adaptive weighted spherical harmonic signals according to any one of claims 1-7, a final spherical harmonic coefficient fusion solution is obtained, which is a best estimation of fused non-tidal high-frequency atmospheric mass variation in a spherical harmonic domain, and atmospheric mass analysis is realized based on the best estimation. The method comprises the following steps: Obtaining a plurality of remote sensing imaging data sets; The remote sensing imaging data set is preprocessed, which includes resolution unification and correction processing; The preprocessed remote sensing imaging data set is taken as multi-source spatio-temporal observation data to be fused, and the multi-source data fusion method based on adaptive weighted spherical harmonic signals is used to fuse the multi-source spatio-temporal observation data to be fused, so that a final spherical harmonic coefficient fusion solution is obtained as a fused remote sensing image. Vegetation monitoring is realized based on the fused remote sensing image.

Images