Hydroclimatic reconstruction method based on multi-layer nested coupling

Abstract

The application relates to the technical field of data processing, and discloses a hydrological climate element history sequence reconstruction method based on multi-layer nested coupling, which comprises the following steps: unifying the time resolution of multi-source substitute data and observation sequences; performing adaptive noise complete ensemble empirical mode decomposition, recombining to generate high-frequency layer and low-frequency layer control sequences; calculating cross wavelet power spectrum, extracting instantaneous phase difference sequences and wavelet square coherence sequences; constructing a nested coupling regression benchmark model, calculating dynamic high-low frequency regression coefficients and generating a two-dimensional system lookup table; using the two-dimensional system lookup table to obtain regression coefficients of a history period, combining with coherent logic gate control sequences generated by the wavelet square coherence sequences to process distorted data, and outputting reconstructed sequences. The application isolates the interference of signals of different frequency bands, adapts to the time-varying properties of climate system parameters, and improves the reconstruction quality of long-term trends and short-term fluctuations.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, specifically to a method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling. Background Technology

[0002] Reconstructing historical sequences of hydro-climate elements is fundamental to exploring past climate evolution patterns. Conventional reconstruction processes primarily utilize multi-source proxy data, establishing statistical transformation models with target climate elements within an instrument calibration period with available observational data. These models are then applied to historical periods without observational records for data extrapolation. However, the natural climate system is a non-stationary evolutionary process, where high-frequency fluctuations and low-frequency background signals are controlled by different driving mechanisms. Existing reconstruction techniques often establish mapping relationships between variables on a single time scale, failing to separate the inherent frequency characteristics of the input sequences. This leads to interference between signals in different frequency bands during modeling, reducing the accuracy of the reconstruction results in reconstructing long-term climate evolution trends and short-term oscillations.

[0003] Meanwhile, traditional reconstruction models often employ linear regression with constant coefficients for fitting. Due to the nonlinear response characteristics within the climate system, the sensitivity of proxy data to hydro-climate elements shifts with changes in the climate background phase, making it difficult for fixed model parameters to adapt to the time-varying nature of these physical processes. Furthermore, over long historical projection periods, proxy data can be affected by external environmental noise or reach physical boundary conditions for biological growth, leading to weakened or distorted climate indicative significance. Current technologies typically lack dynamic monitoring and control mechanisms for such distorted data. Distorted proxy signals are continuously introduced into the regression system, causing errors to propagate step-by-step in the reconstruction sequence, affecting the statistical uniformity of historical hydro-climate reconstruction data. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling. This method solves the problems of existing hydrological and climatic history reconstruction methods failing to separate high- and low-frequency signals and often using constant regression coefficients, which are unable to adapt to the non-stationary processes and time-varying properties of natural climate evolution, and lacking control mechanisms for distorted proxy data, which easily leads to error propagation.

[0005] To achieve the above objectives, the present invention provides a method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling, comprising the following steps: The first aspect of this invention provides a method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling, comprising the following steps: By unifying the temporal resolution of multi-source proxy data sets and hydrological and climatic element observation sequences, a standard input sequence matrix is ​​generated, dividing the calibration period of the device measurement and the historical extrapolation period. Adaptive noise-complete set empirical mode decomposition is performed on each data sequence in the standard input sequence matrix to separate and extract intrinsic mode functions and residual components, and then reassemble them to generate high-frequency layer control sequences and low-frequency layer control sequences. Calculate the cross-wavelet power spectrum of the high-frequency layer control sequence and the low-frequency layer control sequence, lock the dominant periodic frequency band, and extract the instantaneous phase difference sequence and wavelet square coherence sequence within the dominant periodic frequency band; Construct a nested coupled regression benchmark model that includes a benchmark intercept term, and calculate the dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function by combining the instantaneous phase difference sequence to generate a two-dimensional system lookup table; Historical high-frequency regression coefficients and historical low-frequency regression coefficients are obtained based on a two-dimensional system lookup table. A coherent logic gate control sequence is generated based on the wavelet squared coherence sequence. Combining the coherent logic gate control sequence, historical high-frequency regression coefficients, historical low-frequency regression coefficients, high-frequency layer control sequence, low-frequency layer control sequence, and benchmark intercept term, a reconstructed sequence of historical hydrological and climatic elements is output.

[0006] Furthermore, the process of generating a standard input sequence matrix by unifying the time resolution of the multi-source proxy data set and the hydrological and climatic element observation sequences, and dividing the instrument calibration period and the historical extrapolation period, includes: Extract the intersection interval between the hydrological and climatic element observation sequence and the multi-source proxy data set on the time axis, and define the intersection interval as the instrument calibration period; Extract historical time periods from the multi-source proxy dataset that are earlier than the instrument calibration period and are not covered by instrument observation data, and define these historical time periods as the historical projection period. Identify the original sampling intervals of each time series in the multi-source proxy dataset, and perform temporal forced alignment on the multi-source proxy dataset based on a set globally unified reference time step: When the original sampling interval is less than the baseline time step, mean smoothing downsampling is performed along the time axis according to the window span of the baseline time step. When the original sampling interval is greater than the reference time step or there are non-uniform sampling distribution characteristics, numerical interpolation reconstruction is performed along the time axis. After completing the forced alignment operation in the time domain, the time series of each group in the multi-source proxy data set are truncated to a unified start and end boundary in the time dimension and combined to construct a standard input sequence matrix.

[0007] Furthermore, the process of performing adaptive noise-complete set empirical mode decomposition on each data sequence in the standard input sequence matrix, separating and extracting intrinsic mode functions and residual components, and recombining them to generate high-frequency layer control sequences and low-frequency layer control sequences includes: Multiple sets of mutually independent zero-mean Gaussian white noise sequences are added to each data sequence in the standard input sequence matrix, where the standard deviation of the zero-mean Gaussian white noise sequence is the product of the set noise amplitude coefficient and the standard deviation of the corresponding data sequence. Adaptive noise-complete ensemble empirical mode decomposition is performed on the noisy sequence, and the ensemble mean is calculated at each time node. Multiple eigenmode functions with descending frequency characteristics and a residual component characterizing the long-term monotonic climate evolution trend are separated step by step. Calculate the average oscillation period of each intrinsic mode function on the time axis, and set a physical boundary period threshold to define the separation boundary between high-frequency fluctuation signals and low-frequency background signals; Iterate through all intrinsic mode functions, select intrinsic mode functions whose average oscillation period is less than the physical boundary period threshold, and sum and average the selected intrinsic mode functions across the dimensions of the cross-generational data set to obtain the high-frequency layer control sequence; The intrinsic mode functions with an average oscillation period greater than or equal to the physical boundary period threshold are selected. The selected intrinsic mode functions are added to the corresponding residual components, and then a summation and averaging calculation is performed on the dimension of the cross-generational data set to obtain the low-frequency layer control sequence.

[0008] Furthermore, the process of calculating the cross-wavelet power spectrum of the high-frequency layer control sequence and the low-frequency layer control sequence to lock the dominant periodic frequency band includes: Based on the complex Morlet wavelet function, continuous wavelet transform operations are performed on the high-frequency layer control sequence and the low-frequency layer control sequence respectively to obtain the high-frequency layer continuous wavelet transform spectrum and the low-frequency layer continuous wavelet transform spectrum. Calculate the complex product of the complex conjugate of the high-frequency layer continuous wavelet transform spectrum and the low-frequency layer continuous wavelet transform spectrum to generate a cross wavelet power spectrum that includes the scaling factor and translation factor. Within the data region that passes the red noise confidence test in the cross-wavelet power spectrum, find the energy extremum point with the largest absolute value, extract the scale factor corresponding to the energy extremum point as the dominant scale, and lock the dominant periodic frequency band.

[0009] Furthermore, the process of extracting the instantaneous phase difference sequence and the wavelet squared coherence sequence within the dominant periodic frequency band includes: Extract the one-dimensional complex data sequence corresponding to the cross wavelet power spectrum on the dominant scale, perform imaginary and real part separation operations on the one-dimensional complex data sequence respectively, obtain the one-dimensional imaginary part sequence and the one-dimensional real part sequence, calculate the arctangent function of the ratio of the value of the one-dimensional imaginary part sequence to the value of the one-dimensional real part sequence under the same translation factor, and obtain the instantaneous phase difference sequence. A two-dimensional smoothing operation, encompassing translational smoothing along the time axis and weighted smoothing along the scale axis, is performed on the absolute square of the cross wavelet power spectrum at the dominant scale. The smoothed result is then divided by the product of the absolute square of the smoothed high-frequency continuous wavelet transform spectrum at the dominant scale and the absolute square of the smoothed low-frequency continuous wavelet transform spectrum at the dominant scale to obtain the wavelet square coherence sequence.

[0010] Furthermore, the process of constructing a nested coupled regression benchmark model that includes a benchmark intercept term, calculating the dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function by combining the instantaneous phase difference sequence, and generating a two-dimensional system lookup table includes: Using the corresponding sequences of hydrological and climatic element observation sequences during the instrument calibration period as dependent variables, and the corresponding sequences of high-frequency layer control sequences and low-frequency layer control sequences during the instrument calibration period as independent variables, a nested coupled regression benchmark model containing prior high-frequency regression coefficients, prior low-frequency regression coefficients, and benchmark intercept terms is established. Based on the corresponding sequence of the instantaneous phase difference sequence during the instrument measurement and calibration period, piecewise sliding regression is performed on the nested coupled regression benchmark model to extract the maximum and minimum values ​​of the local high-frequency regression coefficients and the maximum and minimum values ​​of the local low-frequency regression coefficients, which constitute the high-frequency parameter modulation interval and the low-frequency parameter modulation interval, respectively. The dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function are expanded using truncated Fourier series with positive integer expansion orders from the third to the fifth order, and a Tikhonov cost function containing a data fitting term and a regularization penalty term is constructed. The data fitting term is the sum of squared residuals between the corresponding sequence of the hydrological and climatic element observation sequence during the instrument calibration period and the reconstructed sequence of the nested coupled regression benchmark model. The regularization penalty term is the square integral of the second derivative of the dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function as a function of the independent phase difference variable, multiplied by the regularization parameter. Using prior high-frequency and low-frequency regression coefficients as initial iterative estimates, and the high-frequency and low-frequency parameter modulation intervals as hard constraint boundaries for optimization, the minimum value of the Tikhonov cost function is solved to obtain the set of Fourier coefficients that meet the convergence conditions. The continuous dynamic high-frequency and low-frequency regression coefficient functions are output, which together constitute a bivariate parameter modulation surface and are used to generate a two-dimensional system lookup table.

[0011] Furthermore, the process of converting the bivariate parameter modulation surface into a two-dimensional system lookup table and performing mesh discretization includes: Set the discrete resolution of the independent phase difference variables, divide the definition interval of the independent phase difference variables into equidistant grids based on the discrete resolution, and obtain the discrete phase difference sequence; Numerical calculations are performed by substituting the coordinates of each node in the discrete phase difference sequence into the dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function, respectively, and calculating the discrete values ​​of the high-frequency regression coefficient and the low-frequency regression coefficient at the corresponding node positions. Perform spline interpolation smoothing operation, use cubic spline interpolation algorithm to smooth the connection of the discretized node data matrix, encapsulate the discrete phase difference sequence after spline interpolation smoothing, the set of corresponding high-frequency regression coefficient discrete values ​​and the set of low-frequency regression coefficient discrete values, construct and store a two-dimensional system lookup table.

[0012] Furthermore, the process of obtaining historical high-frequency regression coefficients and historical low-frequency regression coefficients based on a two-dimensional system lookup table includes: The system iterates through each discrete time node in the historical projection period in chronological order, extracts the current value of the corresponding sequence in the historical projection period of the instantaneous phase difference sequence as the input key value, and calls the two-dimensional system lookup table for retrieval. Traverse downwards in the index column of the two-dimensional system lookup table, lock the two adjacent discrete phase difference nodes that contain the input key value in the numerical range, and extract the discrete values ​​of the high-frequency regression coefficient and the low-frequency regression coefficient corresponding to the two discrete nodes. The distance weighting ratio is calculated based on the absolute value of the difference between the current input key value and the values ​​of two discrete nodes. The extracted discrete values ​​are then weighted and averaged using a one-dimensional linear interpolation algorithm to output the historical high-frequency regression coefficient and historical low-frequency regression coefficient for the corresponding time node.

[0013] Furthermore, the process of generating a coherent logic gate control sequence based on the wavelet squared coherence sequence, and then combining the coherent logic gate control sequence to output the historical hydrological and climatic element reconstruction sequence includes: The wavelet squared coherence sequence is compared time-by-time with the corresponding sequence within the historical projection period to a set coherence threshold. When the corresponding sequence at the current time is greater than or equal to the coherence threshold, the coherence logic gate control sequence is set to 1; When the corresponding sequence at the current time is less than the coherence threshold, the coherence logic gate control sequence is set to 0. Generate a simulated residual sequence that follows a normal distribution and has a mean of zero. The variance of the simulated residual sequence is equal to the average of the sum of squares of the residuals of the nested coupled regression benchmark model obtained during the instrument calibration period. Perform numerical algebraic calculations for integrated nested coupled derivations: The product of the historical high-frequency regression coefficients and the corresponding sequences of the high-frequency layer control sequence in the historical projection period is added to the product of the historical low-frequency regression coefficients and the corresponding sequences of the low-frequency layer control sequence in the historical projection period. The sum of the products is then multiplied by the coherent logic gate control sequence, followed by the addition of the baseline intercept term. Finally, the product of the difference between the value 1 and the coherent logic gate control sequence and the simulated residual sequence is added to output the reconstructed sequence of hydrological and climatic elements in the historical period.

[0014] Furthermore, after outputting the reconstructed sequence of historical hydrological and climatic elements, the process also includes the splicing and output of the full-band sequence: An overlapping time window is set between the historical projection period and the instrument calibration period. The mean deviation between the reconstructed sequence of hydrological and climatic elements in the historical period and the corresponding sequence of observed hydrological and climatic elements in the instrument calibration period is calculated within the overlapping time window. Using the calculated mean deviation, translation compensation is performed on the reconstructed sequence of historical hydro-climate elements to obtain the historical reconstructed sequence after deviation correction. The bias-corrected historical reconstruction sequence and the corresponding sequence of hydrological and climatic element observation sequence during the instrument calibration period are merged in the forward order of time variables to generate a full-band hydrological and climatic reconstruction sequence. A weighted moving average algorithm is adopted, and the sliding step size is set to 3 to 7 time nodes. The low-pass filtering principle is applied to the data segments on both sides of the time boundary point of the full-band hydrological and climate reconstruction sequence. The processed sequence is then encapsulated into a standard data structure for output.

[0015] This invention provides a method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling. It has the following beneficial effects: 1. This invention employs adaptive noise-complete ensemble empirical mode decomposition (EMD) technology to separate and reconstruct multi-source proxy data from observation sequences into high-frequency and low-frequency control sequences. Unlike conventional methods that establish variable transformation relationships on a single time scale, this step adapts to the non-stationary evolution of natural climate systems, isolates high-frequency fluctuation signals from low-frequency background signals, avoids mutual interference between data from different frequency bands during the modeling process, and improves the reconstruction quality of the reconstructed sequences for long-term evolution trends and short-term fluctuations.

[0016] 2. This invention extracts the instantaneous phase difference sequence between the high-frequency and low-frequency layers and combines it with the Tikhonov cost function to solve for the Fourier coefficients, constructing a dynamic nested coupled regression benchmark model and a two-dimensional system lookup table. This mechanism replaces the traditional constant-coefficient linear regression model, adapting to the time-varying parameter properties caused by the nonlinear characteristics within the climate system, and outputs dynamic high- and low-frequency regression coefficients. Simultaneously, the introduction of the two-dimensional system lookup table reduces the complexity of real-time calculations and improves the execution efficiency of large-scale numerical calculations over historical projection periods.

[0017] 3. This invention establishes a control mechanism for distorted data by generating a coherent logic gate control sequence based on wavelet squared coherence sequences. When the climate indication significance of the substitute data is weakened in historical periods due to environmental noise or physical boundary limitations, the system sets the logic gate control sequence to 0 through coherence threshold comparison, stopping the introduction of substitute data for that period, and instead using a simulated residual sequence generated based on the variance of the instrument calibration period residuals for numerical filling. This step limits the propagation of error signals during the reconstruction process, ensuring the uniformity of the statistical distribution of the reconstructed hydrological and climate results for historical periods. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the system architecture of the present invention; Figure 2 This is a schematic diagram of the method flow of the present invention; Figure 3 This is a comparison chart of variance stability during the historical signal breakage period of the present invention. Figure 4 This is a comparison chart of the reconstruction accuracy index and system computational efficiency of the present invention. Detailed Implementation

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

[0020] Please see the appendix Figure 1 This invention provides a system for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling, including: a preprocessing module, a decoupling and layering module, a feature extraction module, a surface construction module, and a sequence reconstruction module.

[0021] The preprocessing module receives the raw data sequence from external input. The output of the preprocessing module is connected to the input of the decoupled layering module.

[0022] The decoupling layering module receives the preprocessed data and outputs high-frequency and low-frequency control sequences. Both the high-frequency and low-frequency outputs of the decoupling layering module are connected to the input of the feature extraction module.

[0023] The feature extraction module generates time-frequency feature parameters based on the high-frequency layer control sequence and the low-frequency layer control sequence. The first parameter output of the feature extraction module is connected to the surface construction module, and the second parameter output of the feature extraction module is connected to the sequence reconstruction module.

[0024] The surface construction module generates a multidimensional data mapping table based on the input time-frequency feature parameters. The data output of the surface construction module is connected to the sequence reconstruction module. The sequence reconstruction module outputs the time-series reconstruction results based on the input mapping table.

[0025] Please see the appendix Figure 2 This invention provides a method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling, comprising the following steps: S10: Obtain the time-overlapping multi-source proxy data set and the hydro-climate element observation sequence during the instrument calibration period. Use the preprocessing module to perform time resolution uniform processing on the input multi-source proxy data set and hydro-climate element observation sequence to generate a standard input sequence matrix. S20: The decoupling layering module is used to perform adaptive noise complete set empirical mode decomposition on each data sequence in the standard input sequence matrix, separate and extract the intrinsic mode functions and residual components, and reorganize the intrinsic mode functions according to the given average period limit to generate high-frequency layer control sequences and low-frequency layer control sequences. S30, the feature extraction module is used to calculate the cross wavelet power spectrum of the high-frequency layer control sequence and the low-frequency layer control sequence, and to extract the instantaneous phase difference sequence and wavelet square coherence sequence of the high-frequency layer control sequence and the low-frequency layer control sequence within the dominant periodic frequency band. S40, during the instrument testing and calibration period, a high-frequency sensitivity modulation model is constructed using the surface construction module. The Tikhonov regularization constraint mechanism is introduced to inversely calculate the high-frequency dynamic sensitivity coefficients at discrete time nodes. The amplitude of the low-frequency layer control sequence is used as the first coordinate axis variable, the instantaneous phase difference sequence is used as the second coordinate axis variable, and the high-frequency dynamic sensitivity coefficients are used as the third coordinate axis variable. Two-dimensional cubic spline interpolation is performed to generate a bivariate parameter modulation surface. The bivariate parameter modulation surface is converted into a lookup table of discrete meshes and stored in the system storage unit. S50, during the historical projection period, uses the sequence reconstruction module to extract the amplitude and instantaneous phase difference of the low-frequency layer control sequence at the current time node, reads the high-frequency sensitivity modulation parameters at the current time node from the lookup table, sets a consistency threshold, compares the wavelet squared coherence of the current time node with the consistency threshold, generates a weight function to control the high-frequency variance term, performs numerical algebra calculations based on the weight function, the high-frequency sensitivity modulation parameters, the amplitude of the low-frequency layer control sequence, and the current value of the high-frequency layer control sequence, and outputs a full-band fused hydrological and climatic element reconstruction sequence.

[0026] The specific implementation process of acquiring a multi-source proxy dataset with temporal overlap and hydro-climatic element observation sequences during the instrument calibration period, and using a preprocessing module to perform time-resolution unification processing on the input multi-source proxy dataset and hydro-climatic element observation sequences to generate a standard input sequence matrix includes the following: Because different types of proxy data obtained in nature often have different growth or deposition rates, their original records have inherent differences in sampling frequency. In order to meet the alignment requirements of the data dimensions of the subsequent decoupled hierarchical module, it is necessary to unify the multi-source heterogeneous time series data to the same time base, thereby constructing a standard input format for matrix operations.

[0027] Multiple sets of heterogeneous historical climate proxy records were collected through the data access interface of the preprocessing module to form a multi-source proxy data set. .in, It is a collection of multi-source substitute data. This is the time series of the first set of proxy data. This is the time series of the second set of substitute data. For the first Time series of group substitute data; The total number of substitute data sets, whose value is a positive integer greater than or equal to 2; Indicates the first Time series of group substitute data For index, its value range is Positive integers; For time variables, representing discrete time points. Multi-source proxy dataset. It covers data types with different response mechanisms. Multi-source proxy data sets. It includes biological and geological proxies. The specific physical manifestations of the biological proxies include tree-ring width sequences and tree-ring maximum latewood density sequences. The specific physical manifestations of the geological proxies include stalagmite oxygen isotope sequences, lacustrine sediment grain size sequences, and ice core accumulation rate sequences. Simultaneously, actual observational data from meteorological stations are collected to construct observational sequences for hydrological and climatic elements. ,in This is an observation sequence of hydrological and climatic elements. The specific data formats include instrument-measured precipitation sequences, instrument-measured average temperature sequences, and standardized precipitation evapotranspiration index sequences.

[0028] Analysis of hydrological and climatic element observation sequences Records the time span and simultaneously analyzes multi-source proxy data sets. Record the time span of each time series. Extract the observation sequences of hydrological and climatic elements. With multi-source proxy data sets The intersection interval on the time axis is defined as the instrument calibration period. ,in For instrument calibration, extract multi-source proxy datasets. Earlier than the instrument calibration period Furthermore, the historical period lacking instrumental observation data is designated as the historical projection period. ,in This is the period of historical deduction.

[0029] Identifying multi-source proxy datasets The original sampling intervals for each time series group are defined. For heterogeneous data with differing original sampling intervals, a globally unified reference time step is established. ,in This is the baseline time step. The time resolution is determined based on the target reconstruction requirements, taking a value of 1 year or 10 years, corresponding to an annual resolution scale or an interdecadal resolution scale, respectively. This is based on the baseline time step. Multi-source proxy datasets Perform a forced time-domain alignment operation. This applies when the original sampling interval of a sequence within a group is less than the baseline time step. At that time, according to the reference time step The window span is subjected to mean-smoothing downsampling along the time axis. When the original sampling interval of a sequence within a group is greater than the baseline time step... When non-uniform sampling distribution characteristics exist, numerical interpolation reconstruction is performed along the time axis. For resampling and basic numerical interpolation calculations of non-uniform time series, those skilled in the art can use cubic spline interpolation or piecewise linear interpolation operations according to the degree of data missing. The basic interpolation logic and mathematical solution principles are well-known technologies in this field and will not be elaborated here.

[0030] After completing the temporal alignment operation, the multi-source substitute data set will be... Each time series in the dataset is truncated to a uniform start and end boundary along the time dimension, and then combined to construct a standard input sequence matrix. The column vectors of the standard input sequence matrix correspond to independent, uniformly resolved proxy data, while the row vectors correspond to discrete reference time nodes. The preprocessing module combines the standard input sequence matrix with the hydrological and climatic element observation sequences that have undergone time boundary alignment. The parameters are output to the decoupled layered module as reference input parameters for subsequent frequency domain analysis.

[0031] The specific implementation process of performing adaptive noise complete set empirical mode decomposition on each data sequence in the standard input sequence matrix using a decoupling hierarchical module to separate and extract intrinsic mode functions and residual components includes the following: Climate proxy sequences obtained from nature possess nonlinear and non-stationary physical properties, containing multi-scale fluctuations caused by different climate-driving mechanisms. Directly fusing multi-source data in the time domain can lead to mutual interference between physical information at different scales. By performing adaptive noise-complete ensemble empirical mode decomposition, the time series is peeled away layer by layer in the frequency domain, which can solve the mode aliasing phenomenon in empirical mode decomposition and achieve the extraction of signals at a single physical scale.

[0032] The first one in the standard input sequence matrix Time series of group substitute data As the basic decomposition object, in For the first Time series of group substitute data For indexing, Use time as the variable. Set the number of times the set is executed. ,in The set number is defined, and its value range is set to 1. A positive integer. Set the noise amplitude coefficient. ,in This is the noise amplitude coefficient, and its value range is set to... .

[0033] Calculate the first Time series of group substitute data The standard deviation of . In the . Time series of group substitute data Add to Construct a set of mutually independent zero-mean Gaussian white noise sequences. A set of noise-added sequences, wherein the standard deviation of the added Gaussian white noise sequence is the noise amplitude coefficient. With the Time series of group substitute data The product of the standard deviations. For this Perform standard empirical mode decomposition on each noisy sequence to extract the first-order mode components obtained from the decomposition of each noisy sequence. The ensemble mean of the first-order modal components is calculated at each time node to generate the first-order modal component. Time series of group substitute data The corresponding first-order eigenmode function ,in For the first The first-order intrinsic mode function corresponding to the time series of the group of proxy data. Time series of group substitute data Subtract the first-order intrinsic mode function Obtain the first-order residual sequence.

[0034] In subsequent decomposition iterations, a new combination of input signals is constructed by adding first-order Gaussian white noise mode features generated by standard empirical mode decomposition to the next-order residual sequence, and local extrema and means are continued to be obtained. For the extreme point envelope fitting and local mean subtraction operations in multiple iterations, those skilled in the art can perform them according to conventional empirical mode decomposition rules; the basic envelope extraction and averaging logic are well-known techniques in the field and will not be elaborated upon here.

[0035] Repeat the above noise injection and mean calculation process to separate the noise level step by step. eigenmode functions ,in For the first The time series corresponding to the group of substitute data eigenmode functions of order; This is the index of the intrinsic mode function, and its range is... Positive integers; This represents the total number of intrinsic mode functions. Without prior specification, the iteration terminates when the total number of local extrema in the residual sequence during the iteration process is less than or equal to two. The output sequence at termination is the residual component. ,in For the first The residual components corresponding to the time series of the group of substitute data.

[0036] After completing the adaptive noise complete set empirical mode decomposition, the first Time series of group substitute data In mathematical structure, it is decomposed into a linear combination form as shown in the following formula: ; In the formula, For the first Time series of group substitute data; It is a time variable; This represents the total number of intrinsic mode functions; This is the index of the intrinsic mode function; For the first The time series corresponding to the group of substitute data eigenmode functions of order; For the first The residual components corresponding to the time series of the proxy data are used. Through this formula, the original non-stationary proxy series is transformed into a superposition of multiple eigenmode functions with descending frequency characteristics and a residual component representing the long-term evolution trend, providing a basic data source for subsequent high- and low-frequency boundary delineation.

[0037] The specific implementation process for reorganizing intrinsic mode functions based on given average period limits to generate high-frequency and low-frequency control sequences includes the following: In the Empirical Mode Decomposition (EMD) process, the resulting eigenmode functions exhibit a frequency-decreasing characteristic from high to low. To rigorously distinguish high-frequency fluctuation signals from low-frequency background signals at the physical level, a quantitative periodicity assessment of each eigenmode function is required. The average period is calculated for each order of eigenmode functions output from the EMD of the adaptive noise complete set. For the [specific example]... The time series corresponding to the group of substitute data eigenmode functions Calculate its average oscillation period on the time axis. .in, For the first The time series corresponding to the group of substitute data eigenmode functions of order 1 For indexing, The index of the intrinsic mode function. For time variables, This represents the average oscillation period. In practical implementation, this is achieved by traversing the... The time series corresponding to the group of substitute data eigenmode functions Extract all local maxima and local minima, and count the total number of extreme points. Divide the total time span of the sequence by the total number of extreme points and multiply by two to obtain the average oscillation period of the modal component. For the identification and statistical calculation of extreme points, those skilled in the art can perform the task according to the conventional rules for determining the derivative of discrete sequences. The basic extreme value statistical logic is a well-known technology in this field and will not be elaborated here.

[0038] Obtain the average oscillation period of all intrinsic mode functions. Then, set the physical boundary period threshold. ,in This is the physical boundary period threshold. Used to define the separation boundary between high-frequency fluctuation signals and low-frequency background signals. When the baseline time step of the multi-source proxy dataset is one year, the physical boundary period threshold is... The value is set to ten years, or based on the climate background state evolution cycle of the study area, it is set to a value greater than or equal to ten years and less than or equal to thirty years, so that interannual climate fluctuations below the decadal scale are classified as high-frequency signals, and background climate baselines at and above the decadal scale are classified as low-frequency signals.

[0039] Based on physical boundary period threshold Reorganize the intrinsic mode functions to generate high-frequency layer control sequences. ,in This is a high-frequency layer control sequence. Iterate through all the... Using time series data from a group of proxy data, the average oscillation period is selected. Less than the physical boundary period threshold The intrinsic mode functions (EMFs) are then summed and averaged across the dimensions of the cross-generational data set to obtain the high-frequency layer control sequence. The specific formula for calculating the combination is as follows: ; In the formula, It is a high-frequency layer control sequence; The total number of groups of substitute data; For indexing; It is a time variable; For the first The time series corresponding to the group of substitute data eigenmode functions of order; This is the index of the intrinsic mode function; The average oscillation period; The physical boundary periodic threshold. High-frequency layer control sequence. Characterizes the extreme value fluctuations of multi-source proxy data on an interannual scale.

[0040] Based on physical boundary period threshold The remaining intrinsic mode functions and residual components are recombined to generate low-frequency layer control sequences. ,in This is a low-frequency layer control sequence. Iterate through all the... Using time series data from a group of proxy data, the average oscillation period is selected. Greater than or equal to the physical boundary period threshold The intrinsic mode functions. The selected intrinsic mode functions are compared with the first... Residual components corresponding to the time series of the group of substitute data Add them together, where the residual components Physically, it characterizes the long-term monotonic climate evolution trend that was not extracted by empirical mode decomposition, forming the core component of the background climate baseline. Subsequently, a summation-averaging calculation is performed across the dimensions of the cross-generational data set to obtain the low-frequency layer control sequence. The specific formula for calculating the combination is as follows: ; In the formula, This is a low-frequency layer control sequence; The total number of groups of substitute data; For indexing; It is a time variable; For the first The time series corresponding to the group of substitute data eigenmode functions of order; This is the index of the intrinsic mode function; For the first The residual components corresponding to the time series of the group of substitute data; The average oscillation period; The physical boundary periodic threshold. Low-frequency layer control sequence. Characterize the long-term background evolution trend of multi-source proxy data on interdecadal and longer timescales.

[0041] The decoupling layer module will generate high-frequency layer control sequences and low-frequency layer control sequence The output is sent to the feature extraction module. By recombining the components of different periodic attributes, the decoupled hierarchical module constructs a high- and low-frequency dual-layer control architecture in the data transmission flow.

[0042] The specific implementation process of calculating the cross-wavelet power spectrum of the high-frequency layer control sequence and the low-frequency layer control sequence using the feature extraction module and locking the dominant periodic frequency band includes the following: Since high-frequency fluctuations and low-frequency background baselines in the climate system are not isolated, they undergo dynamic response delays and energy coupling during their physical evolution. Traditional Fourier transforms can only expand signals in the global frequency domain and cannot capture the local changes in signals over time. To quantify the physical energy resonance intensity between high-frequency and low-frequency control sequences and their dynamic evolution over time in the two-dimensional time-frequency domain, it is necessary to construct a cross-wavelet power spectrum to provide a computational basis for subsequent extraction of phase difference and consistency indices.

[0043] Obtain the high-frequency layer control sequence output by the decoupling layer module With low-frequency layer control sequence .in, It is a high-frequency layer control sequence. It is a low-frequency layer control sequence. The time variable is used. The basis functions of the continuous wavelet transform are defined as complex Morlet wavelet functions. Complex Morlet wavelet functions contain complex structures and possess the mathematical property of simultaneously extracting amplitude and phase information from discrete time series. Based on the complex Morlet wavelet function, the high-frequency layer control sequence is processed... and low-frequency layer control sequence Perform continuous wavelet transform operations to obtain the high-frequency layer continuous wavelet transform spectrum. With low-frequency layer continuous wavelet transform spectrum .in, For high-frequency layer continuous wavelet transform spectrum, The low-frequency layer continuous wavelet transform spectrum, As a scale factor, Translation factor. Scale factor. The range of values ​​for is determined by the total time span of the input sequence and the minimum resolvable period defined by the Nyquist sampling theorem, and is physically proportional to the period of the wave signal. (Shift factor) Characterizing the translation position of the wavelet basis function on the time axis, it corresponds to a specific discrete-time variable in discrete sequence analysis. .

[0044] After obtaining the continuous wavelet transform spectrum of a single sequence, calculate the continuous wavelet transform spectrum of the high-frequency layer. With low-frequency layer continuous wavelet transform spectrum The complex product generates the cross wavelet power spectrum. Specifically, the high-frequency layer continuous wavelet transform spectrum The specific formula for multiplying by the complex conjugate of the low-frequency layer continuous wavelet transform spectrum is as follows: ; In the formula, The cross-wavelet power spectrum; The high-frequency layer is a continuous wavelet transform spectrum; It is the complex conjugate of the continuous wavelet transform spectrum of the low-frequency layer; Scale factor; The translation factor. Cross wavelet power spectrum. The absolute value of represents the common physical energy density of the high-frequency layer control sequence and the low-frequency layer control sequence under the corresponding scaling factor and translation factor. For edge effect elimination and red noise background spectrum testing in continuous wavelet transform, those skilled in the art can perform the tasks according to conventional wavelet influence cone and red noise confidence test rules. The basic edge correction and statistical test logic are well-known techniques in this field and will not be elaborated upon here.

[0045] Because the response bands of multi-source proxy data contain multiple weak noise extrema, in order to extract the dominant physical signals representing the main climate drivers, it is necessary to cross-wavelet power spectrum analysis. Mid-locking dominant scale ,in Dominant scale. Traverse the cross-wavelet power spectrum. All scale factors Translation factor The absolute value data matrix is ​​used. Within the valid data region that passes the red noise confidence test, the energy extremum point with the largest absolute value is found. The scale factor corresponding to this energy extremum point is extracted and used as the dominant scale. Dominant Scale During the instrument calibration and historical projection periods, it serves as a unified reference frequency band for subsequent one-dimensional data extraction of instantaneous phase difference sequences and wavelet squared coherence sequences. Through a dominant scale locking mechanism, random physical noise in non-dominant frequency bands is filtered out, ensuring that the time-frequency characteristic parameters received by the surface construction module reflect the energy coupling relationship between high and low frequency sequences.

[0046] The specific implementation process of extracting the instantaneous phase difference sequence and wavelet squared coherence sequence of the high-frequency layer control sequence and the low-frequency layer control sequence within the dominant periodic frequency band using the feature extraction module includes the following: There is a physical propagation time delay between high-frequency fluctuations and low-frequency background in the climate system. Extracting the instantaneous phase difference sequence can quantify the dynamic lead and lag relationship between the high-frequency layer control sequence and the low-frequency layer control sequence in the dominant frequency band.

[0047] Extracting the cross wavelet power spectrum In dominant scale The corresponding one-dimensional complex data sequence. For the cross wavelet power spectrum, As a scale factor, The translation factor is... The dominant scale is used. Since the cross-wavelet power spectrum is constructed based on the complex Morlet wavelet, its calculation result is in complex form. The imaginary and real parts are separated from this one-dimensional complex data sequence to obtain independent one-dimensional imaginary and real part sequences. The instantaneous phase difference sequence is obtained by calculating the arctangent function of the ratio of the imaginary to the real part values ​​under the same translation factor. ,in This is a sequence of instantaneous phase differences. The specific formula for its combination calculation is as follows: ; In the formula, It is a sequence of instantaneous phase differences; It is the arctangent function; Operators for extracting the imaginary part of complex numbers; Operators for extracting the real part of complex numbers; The cross-wavelet power spectrum at the dominant scale is the cross-wavelet power spectrum. Middle scale factor The complex sequence extracted at that time; As the dominant scale; This is the translation factor. Based on the mathematical properties of the arctangent function, the instantaneous phase difference sequence... The range of values ​​for is Radius. When the extracted value is positive, it indicates that the high-frequency layer control sequence leads the low-frequency layer control sequence in evolution phase; when the extracted value is negative, it indicates that the high-frequency layer control sequence lags the low-frequency layer control sequence in evolution phase.

[0048] In long-term records of multi-source proxy data, the response relationship between high- and low-frequency components exhibits a break along the time axis due to interference from local environmental noise. To quantify the physical consistency strength of high- and low-frequency signals in the local time domain and avoid introducing pure noise signals during historical extrapolation, a feature extraction module is used to calculate the wavelet squared coherence sequence. ,in This is a wavelet squared coherence sequence. Wavelet squared coherence is physically equivalent to the mapping of local correlation coefficients in the two-dimensional time-frequency domain, and is used to evaluate the degree of covariance between two sequences in the time-frequency space.

[0049] For dominant scale The absolute square of the cross-wavelet power spectrum is then subjected to two-dimensional smoothing, and divided by the product of the absolute square of the smoothed high-frequency continuous wavelet transform spectrum and the absolute square of the smoothed low-frequency continuous wavelet transform spectrum. The specific combined calculation formula is as follows: ; In the formula, It is a wavelet squared coherence sequence; For time-frequency smoothing operators; Cross wavelet power spectrum at the dominant scale; The high-frequency layer continuous wavelet transform spectrum at the dominant scale is the high-frequency layer continuous wavelet transform spectrum. Middle scale factor The complex sequence extracted at that time; The low-frequency layer continuous wavelet transform spectrum at the dominant scale is the low-frequency layer continuous wavelet transform spectrum. Middle scale factor The complex sequence extracted at that time; As the dominant scale; The translation factor; For absolute value operations. Wavelet squared coherence sequence. The numerical value is normalized to be between 0 and 1. The closer the value is to 1, the stronger the consistency of the physical evolution of high and low frequency signals at the corresponding time point.

[0050] If the time-frequency smoothing operator is not introduced in the wavelet coherence calculation The numerator and denominator will mathematically cancel each other out, resulting in a result that is always 1. Therefore, the time-frequency smoothing operator It covers translational smoothing mechanisms along the time axis and weighted smoothing mechanisms along the scale axis. For time-frequency smoothing operators... The specific implementation method can be carried out by those skilled in the art in accordance with the rule of combining time-domain window filters and scale-domain window filters. The basic two-dimensional moving weighted average filter logic and mathematical derivation are well-known technologies in this field and will not be elaborated here.

[0051] After the feature extraction module completes the parsing of the aforementioned basic parameter matrix, it extracts the instantaneous phase difference sequence containing physical delay information. Square coherence sequence The independent coherence input to the sequence reconstruction module provides temporal characteristic variables for the three-dimensional mesh construction of the bivariate parameter modulation surface and the injection of logic gate variance.

[0052] The specific implementation process of constructing a nested coupled regression benchmark model and extracting prior parameters during the instrument calibration period using the surface construction module includes the following: The response mechanisms of high-frequency extreme fluctuations and low-frequency background evolution in multi-source proxy data to hydrological and climatic elements differ. To quantify the intensity of this independent response, it is necessary to establish a baseline mapping relationship between high- and low-frequency components and hydrological and climatic elements, and extract prior sensitivity parameters, during the instrument calibration period when measured data is available, treating high- and low-frequency components as independent driving variables.

[0053] Instrument calibration period based on preprocessing stage Extract data segments within the corresponding time axis interval, where This is the instrument calibration period. The observation sequences of hydrological and climatic elements during this period will be acquired. High-frequency layer control sequence during instrument calibration period and the low-frequency layer control sequence during the instrument calibration period .in, For the calibration period of hydrological and climatic element observation sequences, This is the high-frequency layer control sequence during the calibration period. This is the low-frequency layer control sequence during the calibration period. It is a time variable.

[0054] The calibration period hydro-climatic element observation sequence As the dependent variable, the high-frequency layer control sequence during the calibration period and calibration period low-frequency layer control sequence As independent variables, a nested coupled regression benchmark model is established using the surface construction module. The nested coupled regression benchmark model is represented as a multivariate linear equation containing both high- and low-frequency input terms. Its specific combination calculation formula is as follows: ; In the formula, This is the observation sequence of hydrological and climatic elements during the calibration period; This is the high-frequency layer control sequence for the calibration period; This is the low-frequency layer control sequence during the calibration period; These are the prior high-frequency regression coefficients; These are the prior low-frequency regression coefficients; The baseline intercept term; The model residual sequence is used to characterize other physical environment noise that is not explained by the high-frequency layer control sequence and the low-frequency layer control sequence; It is a time variable.

[0055] For this nested coupled regression benchmark model, ordinary least squares (OLS) is used to solve for the unknown parameters in the equations. For parameter estimation using OLS, those skilled in the art can perform the task by minimizing the sum of squares of the model residual sequence; the underlying linear algebraic solution logic is well-known in the field and will not be elaborated upon here.

[0056] Solve the above equation to obtain the prior high-frequency regression coefficients. and prior low-frequency regression coefficients Prior high-frequency regression coefficients Characterizing the baseline sensitivity of hydro-climate elements to high-frequency fluctuations. Prior low-frequency regression coefficients. The baseline sensitivity of hydro-climate elements to low-frequency background is characterized.

[0057] Because the driving effect of high- and low-frequency components on climate elements in a climate system is not constant, but rather fluctuates non-stationarily with the degree of phase misalignment between the high- and low-frequency components, instantaneous phase difference sequences are extracted to determine the boundary range of parameter variation with phase. During the instrument calibration period Within the data segment, obtain the instantaneous phase difference sequence during the calibration period. .in, It is an instantaneous phase difference sequence. This is the instantaneous phase difference sequence during the calibration period. The shift factor is based on the instantaneous phase difference sequence during the calibration period. A piecewise moving regression was performed on the nested coupled regression benchmark model. The time span of the moving window was set, with a value ranging from 15 to 30 years (positive integers), specifically matched to the total length of the instrument calibration period data and the typical interdecadal cycle of local climate evolution. As the moving window shifted along the time axis, local high-frequency and low-frequency regression coefficients were calculated within each window. The maximum and minimum values ​​of the local high-frequency regression coefficients were extracted to form the high-frequency parameter modulation interval. The maximum and minimum values ​​of the local low-frequency regression coefficients were also extracted to form the low-frequency parameter modulation interval.

[0058] The surface construction module will use prior high-frequency regression coefficients Prior low-frequency regression coefficients Reference intercept item The corresponding high-frequency parameter modulation range and low-frequency parameter modulation range are stored as prior constraints to provide numerical boundaries for generating a bivariate parameter modulation surface covering the entire phase difference domain.

[0059] The specific implementation process of using the surface construction module to perform inverse solving within the bivariate parameter modulation interval, eliminating parameter singularities and generating a continuous modulation surface includes the following: After extracting prior parameters and local modulation intervals, relying solely on local data segments for piecewise parameter fitting can easily lead to parameter singularities when the phase difference sequence undergoes abrupt changes. This results in discontinuities in the parameter function within the phase space, causing mathematical overfitting. To obtain a globally continuous mapping relationship between high- and low-frequency regression coefficients and phase difference, and to avoid severe oscillations in the fitting curve caused by environmental noise interference, a Tikhonov regularization mechanism is introduced to perform inverse optimization.

[0060] Setting dynamic high-frequency regression coefficient function With dynamic low-frequency regression coefficient function ,in As an independent phase difference variable, its value range is set to... Radius. Since the phase difference possesses cyclic properties, a truncated Fourier series is used to expand these two dynamic coefficient functions, transforming the fitting process of the nonlinear curve into solving a set of discrete Fourier coefficients. In practice, the order of the truncated Fourier series expansion is set to a positive integer between third and fifth order to ensure that while capturing smooth low-frequency changes, high-frequency random disturbances are filtered out. The basic expansion and order reduction rules for the truncated Fourier series can be performed by those skilled in the art according to conventional signal harmonic analysis rules; these are well-known techniques in the field and will not be elaborated upon here.

[0061] Construct a Tikhonov cost function that includes a data fitting term and a regularization penalty term. The data fitting term is the sum of squared residuals between the observed data and the reconstructed sequence from the dual-input model during the instrument calibration period, and the regularization penalty term is the dynamic regression coefficient function varying with the independent phase difference variable. The square integral of the changing second derivative.

[0062] In the cost function In, regularization parameters It is used to balance data fitting accuracy and parametric surface smoothness. In practical implementation, the regularization parameter... The value of is determined based on the generalized cross-validation algorithm, and the regularization strength is locked by finding the extreme point that minimizes the variance of the prediction error. The integral penalty term in the cost function minimizes the curvature of the regression coefficient curve throughout the entire phase domain, thereby eliminating abrupt slope changes and singular extreme points caused by local data interference.

[0063] Solve the Tikhonov cost function using the interior-point method or the trust region reflection optimization algorithm. The minimum value of the coefficients is obtained to acquire the set of Fourier coefficients that satisfy the convergence condition. During the iterative solution process, the extracted prior high-frequency regression coefficients are... and prior low-frequency regression coefficients The initial iterative estimate is used as the constant term. The high-frequency parameter modulation interval and the low-frequency parameter modulation interval are used as hard constraint boundaries for optimization. When the dynamic regression coefficient value generated by the iteration exceeds the corresponding modulation interval, a boundary truncation operation is performed, resetting it to the critical value of the interval.

[0064] After convergence, the obtained Fourier coefficients are substituted into the original function form to output a continuous dynamic high-frequency regression coefficient function. With dynamic low-frequency regression coefficient function These two functions together form a bivariate parametric modulation surface in the three-dimensional space of the phase difference domain. This surface provides a numerical mapping relationship for the system to extract the corresponding high-frequency drive weights and low-frequency drive weights based on the input instantaneous phase difference values ​​during subsequent historical simulations.

[0065] The specific implementation process of using the surface construction module to discretize continuous parametric modulated surfaces and generate system lookup tables includes the following: The dynamic high-frequency regression coefficient function obtained by inverse solution using Tikhonov regularization With dynamic low-frequency regression coefficient function It is a continuous analytic expression. Where... It is a dynamic high-frequency regression coefficient function. This is a dynamic low-frequency regression coefficient function. These are independent phase difference variables. When performing long-term historical climate series extrapolation, analytically calculating continuous functions at each time point consumes significant computational resources. To improve the computational efficiency of parameter extraction, the continuous function model needs to be converted into a discrete data index structure. Discretizing the continuous function and creating a lookup table works by trading memory space for program execution time, replacing complex calculus calculations with table lookups during the extrapolation phase.

[0066] Set the discrete resolution of the phase difference domain and perform mesh discretization. Define independent phase difference variables. discrete resolution ,in Discrete resolution. The magnitude of this value determines the granularity of the parameter mapping, and its value is set to a real number between 0.01 and 0.05 radians, based on the system memory capacity and the accuracy requirements of historical extrapolation. This is based on the set discrete resolution. For independent phase difference variables Definition range Perform equidistant mesh generation to obtain the included Discrete phase difference sequence of nodes .in It is a discrete phase difference sequence. For node indexing, The total number of discrete nodes is calculated as the total length of the interval. Divide by discrete resolution Round down and add one, node index The value range is 1 to Positive integers between [a certain range].

[0067] Perform numerical computation and spline interpolation smoothing operations. Transform the discrete phase difference sequence... Substitute the coordinates of each node into the dynamic high-frequency regression coefficient function. With dynamic low-frequency regression coefficient function In the process, the discrete values ​​of the high-frequency regression coefficients at the corresponding node positions are obtained through calculation. Discrete values ​​of low-frequency regression coefficients ,in These are the discrete values ​​of the high-frequency regression coefficients. These are the discrete values ​​of the low-frequency regression coefficients. To ensure the continuity of the numerical derivatives between discrete nodes and avoid truncation errors within the lookup table step size interval, a cubic spline interpolation algorithm is used to smoothly connect the discretized node data matrix, providing a smooth transition between discrete nodes and preventing numerical jumps when no precise node is found. The construction of the basis functions and the setting of the boundary conditions in the cubic spline interpolation algorithm can be performed by those skilled in the art according to conventional piecewise polynomial fitting rules. The basic matrix chasing method solution logic is a well-known technique in this field and will not be elaborated here.

[0068] Construct and store a two-dimensional system lookup table. The discrete phase difference sequence, smoothed by spline interpolation... The corresponding high-frequency regression coefficient discrete values The set of low-frequency regression coefficient discrete values The set of data is encapsulated to generate a two-dimensional system lookup table. ,in This is a two-dimensional system lookup table. The mapping relationship at the data structure level satisfies the following formula: ; In the formula, This is the parameter output vector for the corresponding discrete node; It is a discrete phase difference sequence; These are the discrete values ​​of the high-frequency regression coefficients; These are the discrete values ​​of the low-frequency regression coefficients.

[0069] The generated two-dimensional system lookup table follows the steps described above. It is stored in the memory buffer area of ​​the sequence reconstruction module. During subsequent historical sequence extrapolation, when the sequence reconstruction module receives the instantaneous phase difference input variable at any historical moment, it no longer calls the original function for analytical calculation. Instead, it uses a rounding algorithm or a linear interpolation algorithm to look up the table in the two-dimensional system. Retrieve and output the corresponding discrete values ​​of high-frequency regression coefficients. Discrete values ​​of low-frequency regression coefficients By using a two-dimensional system lookup table, the integration and partial derivative operations of nonlinear functions are transformed into basic memory address mapping and data reading operations, reducing the underlying computing overhead of the system.

[0070] The specific implementation process of using the sequence reconstruction module to execute a dynamic parameter lookup mechanism and extract time-varying weight coefficients during the historical projection period includes the following: After constructing the two-dimensional system lookup table, the system enters the historical extrapolation stage where measured hydrological and climatic observation data is lacking. In this stage, the response intensity of high-frequency extreme fluctuations and low-frequency background evolution to hydrological elements is not a static constant, but rather dynamically shifts with the instantaneous phase difference between high and low-frequency components. Therefore, it is necessary to extract the corresponding dynamic regression coefficients based on the phase difference state at historical moments through a memory lookup table mechanism to construct a non-stationary extrapolation model. Transforming the complex climate-physical response relationship into computer memory read operations can reduce computational power consumption during long-term historical data iterations.

[0071] Historical projection period based on preprocessing stage Extract data segments within the corresponding time axis interval, where For the historical projection period, the data matrix output by the feature extraction module is extracted to obtain the high-frequency layer control sequences within the historical projection period. Low-frequency layer control sequences during the historical projection period and the instantaneous phase difference sequence during the historical projection period .in, This is a high-frequency layer control sequence from the historical period; This is a low-frequency layer control sequence from the historical period; This is a sequence of instantaneous phase differences over a historical period, derived from a wavelet shift factor, corresponding to specific time variables. ; It is a time variable.

[0072] The sequence reconstruction module traverses the historical projection period in chronological order. Each discrete time node within At each point in time. Extracting the instantaneous phase difference sequence of the historical period The current value. Using this value as the input key, the system lookup table stored in the memory buffer is invoked. Perform a search, among which This is a lookup table for a two-dimensional system.

[0073] Due to the two-dimensional system lookup table The grid nodes are constructed using discrete phase difference sequences, with the step size determined by the discrete resolution. The decision is made because the input historical instantaneous phase difference values ​​are continuously calculated floating-point numbers, which in most cases lie between two adjacent discrete grid nodes, failing to hit the precise address index. To obtain smooth parameter values ​​and avoid step errors in parameter extraction, the system embeds a one-dimensional linear interpolation algorithm during the table lookup process. Specifically, the program searches the table in a two-dimensional system. Traverse downwards through the index column to locate the two adjacent discrete phase difference nodes within the numerical range that contain the current input key value. Extract the discrete values ​​of the high-frequency regression coefficients corresponding to these two discrete nodes. Discrete values ​​of low-frequency regression coefficients ,in These are the discrete values ​​of the high-frequency regression coefficients. These are the discrete values ​​of the low-frequency regression coefficients. Subsequently, the distance weighting ratio is calculated based on the absolute value of the difference between the current input key value and the values ​​of these two discrete nodes, and the extracted discrete values ​​are then weighted and averaged. For the allocation of distance weights and the multiplication-addition operations in the one-dimensional linear interpolation algorithm, those skilled in the art can perform the operation according to conventional numerical interpolation rules; the basic geometric ratio conversion is a well-known technique in the field and will not be elaborated upon here.

[0074] Through the above table lookup and interpolation calculations, the sequence reconstruction module outputs the corresponding time node. Historical high-frequency regression coefficients Low-frequency regression coefficients of historical periods .in, These are the high-frequency regression coefficients for the historical period. These are low-frequency regression coefficients for historical periods. As the time variable... As the process progresses, the table lookup operation is executed cyclically to generate data covering the entire historical projection period. Two sets of one-dimensional time series data were obtained. These two dynamic parameter sequences reflect the non-stationary changes in the sensitivity of hydrological and climatic elements to high-frequency fluctuations and low-frequency background as phase evolves during historical periods without observational data, providing time-varying weight coefficients that evolve over time for the nested coupled reconstruction model.

[0075] The sequence reconstruction module performs coherent logic gate interception and variance injection during the historical extrapolation period to complete the specific implementation process of nested coupled extrapolation, including the following: In long-term records of multi-source proxy data, the physical consistency between high-frequency and low-frequency signals exhibits localized breaks in time due to interference from local environmental noise. During these noise-dominated periods, the response relationship between high-frequency and low-frequency signals fails. At this time, the two sequences lack a physical covariance relationship. Substituting the control sequence into the model for calculation would result in the reconstruction results being dominated by environmental noise. To prevent the reconstruction curve from becoming a straight line due to a lack of driving force during the consistency failure period, thus causing overall sequence variance collapse, the system introduces a logic gate judgment mechanism based on energy resonance intensity to block calculations in the consistency failure interval and inject background variance conforming to statistical characteristics.

[0076] Extracting the wavelet squared coherence sequence output by the feature extraction module During the historical deduction period Within the data segment, obtain the historical wavelet squared coherence sequence. .in, It is a wavelet squared coherence sequence; This is the period of historical deduction; This is a historical wavelet squared coherence sequence, derived from a wavelet shift factor, corresponding to the time variable. ; This is a time variable. A coherence threshold is set to determine whether the signal has failed. ,in This is the coherence threshold. The numerical determination relies on Monte Carlo simulation testing of the red noise background spectrum, extracting the critical coherence value corresponding to the 95% confidence interval, which is a real number ranging from 0.6 to 0.85. The Monte Carlo simulation testing process for the red noise background spectrum can be performed by those skilled in the art according to conventional wavelet coherence significance testing rules, which is a well-known technique in the field and will not be elaborated upon here.

[0077] Construct a binarized coherent logic gate control sequence. Then, use the historical wavelet squared coherence sequence... With coherence threshold Perform time-sequential numerical comparisons. Construct a coherent logic gate control sequence. ,in This is a coherent logic gate control sequence. When the time variable... Below Greater than or equal to At that time, it is determined that the high and low frequency signals are in a state of physical resonance, and so on. The value is 1; when the time variable Below Less than At that time, it is determined that the signal response relationship at the current moment is broken, and so on. The value is 0.

[0078] The deduction of the integrated nested coupling of execution logic interception and variance injection. The sequence reconstruction module is based on dynamic parameters obtained from table lookup, the input driving sequence, and the coherent logic gate control sequence. The final reconstruction of the hydrological and climatic element sequences is performed. The system generates the simulation residual sequences. ,in This is a simulated residual sequence. It exhibits as random white noise following a normal distribution, with a mean of zero and a variance equal to the model residual sequence calculated by the surface construction module during the instrument calibration period. The average of the sum of squares. Using the above variables, the specific combined calculation formula for the nested coupling deduction is as follows: ; In the formula, This is a reconstructed sequence of historical hydrological and climatic elements; This is a coherent logic gate control sequence; These are high-frequency regression coefficients from historical periods; This is a high-frequency layer control sequence from the historical period; These are low-frequency regression coefficients for historical periods; This is a low-frequency layer control sequence from the historical period; The baseline intercept term; To simulate the residual sequence; It is a time variable.

[0079] Through the logical judgment and switching of the above formula, when the consistency of the climate driving signal is strong, that is... The value of is 1. The system uses dynamic high-frequency components and low-frequency components for multiply-add coupling reconstruction. Zero, simulated residual sequence The input weights are set to zero; when environmental noise causes a break in consistency, i.e. When the value is 0, the system blocks the coupling input of dynamic components and resets the model output value to the baseline intercept term. Superimposed simulated residual sequences This interception and injection mechanism can filter out purely noise-dominated signals while maintaining the variance stationarity of the reconstructed sequence on the time axis, outputting a historical hydrological and climatic element reconstructed sequence that conforms to statistical characteristics.

[0080] After completing the sequence calculation for the historical projection period using the sequence reconstruction module, the specific implementation process of splicing and outputting the full-band sequence includes the following: Obtain the reconstructed sequence of historical hydrological and climatic elements generated through comprehensive nested coupling inference. This sequence covers the historical projection period on the timeline. Extracting the hydro-climate element observation sequences retained during the calibration period from the preprocessing stage. This sequence covers the instrument calibration period on the timeline. .in, This is a reconstructed sequence of historical hydrological and climatic elements; This is the period of historical deduction; This is the observation sequence of hydrological and climatic elements during the calibration period; For instrument calibration period; It is a time variable.

[0081] Because the reconstructed model is based on statistical regression methods, the baseline level of its predictions may be shifted due to physical residuals not explained by the model. If historical hydrological and climatic elements are directly reconstructed into a series... With the hydro-climate element observation sequence during the calibration period When performing end-to-end splicing, numerical discontinuities can occur at the time boundaries. The system introduces a benchmark alignment mechanism based on overlapping time windows to eliminate such systematic biases by aligning the DC baseline.

[0082] During the historical deduction period Instrument calibration period An overlapping time window is set between the points, with a time span of 5 to 10 time steps including the points before and after the time boundary. A reconstructed sequence of historical hydrological and climatic elements is calculated. With the hydro-climate element observation sequence during the calibration period The mean deviation within this overlapping time window. Using the calculated mean deviation, the historical hydro-climate element sequences are reconstructed. Perform translation compensation to obtain the historical reconstruction sequence after deviation correction. .in, The historical reconstructed sequence after bias correction is numerically calculated by adding the aforementioned mean bias constant to the original sequence.

[0083] The historical reconstruction sequence after bias correction With the hydro-climate element observation sequence during the calibration period Along time variable The sequences are merged in a forward order to generate a full-band hydro-climate reconstruction sequence. The combined calculation formula for this splicing process is as follows: ; In the formula, This is a full-band hydro-climate reconstruction sequence; This is the historical reconstructed sequence after bias correction; This is the period of historical deduction; This is the observation sequence of hydrological and climatic elements during the calibration period; For instrument calibration period; The time variable is used. This full-band hydro-climate reconstruction sequence The time span covers the time domain from the beginning of the historical projection period to the end of the instrument measurement and calibration period.

[0084] To eliminate local high-frequency oscillations at the splicing boundary, a weighted moving average algorithm is used to perform smoothing filtering on the data segments on both sides of the time boundary. By setting a sliding step size, a low-pass filtering principle is applied to filter out abrupt noise in the boundary region, where the sliding step size is set to 3 to 7 time nodes. The weight allocation and calculation of the weighted moving average algorithm can be performed by those skilled in the art according to conventional signal smoothing filtering rules; its basic convolution operation logic is well-known in the field and will not be elaborated here.

[0085] Full-band hydro-climate reconstruction sequence after benchmark alignment and smoothing filtering It is encapsulated into a standard data structure. In the specific implementation, the system packages the sequence data and corresponding timestamp metadata, outputs it as a file in comma-separated value format, a common network data format, or a hierarchical data format, and stores it in a local database. The output is a full-band hydro-climate reconstruction sequence. It can be read by external climate assessment models or hydrological analysis systems to analyze the interdecadal evolution of the climate system.

[0086] This invention provides an electronic device for executing the aforementioned climate sequence reconstruction method. This electronic device, through the cooperation of physical components at the hardware level, provides computational resources and storage space for signal decomposition, surface construction, and nested coupled inference.

[0087] The electronic device includes at least one processor and a memory communicatively connected to the at least one processor. Data interaction and instruction transfer between the processor and the memory are achieved via a system bus. The system bus, serving as a data transmission channel between hardware components, includes a data bus, an address bus, and a control bus, used to transmit electrical signals and addressing information between various hardware components, ensuring that the processor can retrieve instructions and data to be processed from the memory.

[0088] As the computing core of electronic devices, processors can take the physical form of central processing units, graphics processing units, microprocessors, digital signal processors, or application-specific integrated circuits. Processors are used to parse computer instructions stored in memory and call upon computing resources to perform numerical computation tasks, including ensemble empirical mode decomposition, continuous wavelet transform, truncated Fourier series expansion, Tikhonov regularized extremum solving, and two-dimensional lookup table interpolation.

[0089] Memory, as a non-transitory computer-readable storage medium, is used to store software programs, instruction sets, and the code of various computing modules. The physical implementation of memory includes volatile and non-volatile memory, specifically encompassing random access memory (RAM), read-only memory (ROM), flash memory, solid-state drives (SSDs), or hard disk drives (HDDs). Internally, memory is divided into multiple memory buffer areas, used for persistently storing preprocessed multi-source proxy data, wavelet coherence matrices, and discretized two-dimensional system lookup tables. As well as the residual sequences and full-band hydro-climate reconstruction sequences generated at each extrapolation stage in the temporary cache. .in, For a two-dimensional system lookup table; This is a full-band hydro-climate reconstruction sequence; It is a time variable.

[0090] The electronic device also includes communication interfaces and input / output interfaces. The communication interface includes a wired Ethernet interface or a wireless radio frequency module, used to establish a network connection with an external meteorological database server, download and update multi-source substitute data and meteorological observation data, and transmit the calculated standard data structure file to a cloud storage center. The input / output interfaces are used to connect external devices, including a monitor, keyboard, and mouse, to receive initial configuration parameters input by the operator, such as independent phase difference variables. discrete resolution or coherence threshold The sequence evolution curves are displayed in a graphical interface. These are independent phase difference variables; For discrete resolution; This is the coherence threshold.

[0091] The memory stores computer instructions that can be executed by the processor. When the processor executes these instructions, the electronic device is configured to instantiate the aforementioned preprocessing module, feature extraction module, surface construction module, and sequence reconstruction module, sequentially performing data filtering, frequency band stripping, phase difference extraction, parametric modulation surface generation, and sequence deduction under logic gate control. During execution, the processor reads data from the memory, performs matrix operations and logical judgments corresponding to the aforementioned algorithms, and then writes the output results back to the designated address segment of the memory, realizing hardware-software interaction.

[0092] For the physical circuit design, instruction set calling mechanism, and driver programming of the underlying architecture of computer hardware, those skilled in the art can follow the conventional computer system integration rules. The basic microprocessor architecture and bus communication protocol are well-known technologies in this field and will not be elaborated here.

[0093] This invention provides a non-transitory computer-readable storage medium. The non-transitory computer-readable storage medium stores computer program instructions, which, when executed by a processor of an electronic device, implement the climate sequence reconstruction method detailed in the foregoing embodiments.

[0094] For the physical form of non-transitory computer-readable storage media, its underlying features include magnetic storage media, optical reading media, and solid-state semiconductor storage media. Specifically, the storage medium can be a portable universal serial bus flash drive, read-only optical disc, digital multifunction optical disc, hard disk drive, solid-state drive, erasable programmable read-only memory, or magnetic tape. These storage media store software program code in binary data within the physical medium by changing the local magnetic field polarity, surface optical reflectivity, or the charge state of semiconductor floating-gate transistors, preventing data loss due to power failure.

[0095] The computer program instructions stored in the storage medium are generated from source code written in a high-level programming language after compilation and linking by a compiler. When the actual simulation task starts, the underlying operating system reads the corresponding binary machine code from this non-transitory computer-readable storage medium via the input / output bus, loads it page by page into volatile memory, and queues and decodes it in the processor's instruction register. Based on the decoded opcode and operands, the processor calls the corresponding hardware arithmetic logic unit, sequentially performing ensemble empirical mode decomposition on the input multi-source proxy data, continuous wavelet transform on the extracted intrinsic mode functions, and utilizing the generated two-dimensional system lookup table. Perform parameter interpolation and integrated nested coupled calculations for the historical projection period, where... This is a lookup table for a two-dimensional system.

[0096] During instruction execution, the electronic device loads multi-source proxy data into memory, the processor performs mathematical operations according to the instruction logic, temporarily stores the generated intermediate variables, and finally generates and outputs a full-band hydro-climate reconstruction sequence. ,in This is a full-band hydro-climate reconstruction sequence. It is a time variable.

[0097] By encapsulating mathematical algorithms in a non-transitory computer-readable storage medium, the stripping of multi-band climate signals, surface fitting of dynamic regression coefficients, and sequence calculation processes based on coherent logic gate control are converted into standardized computer executable files. This encapsulation mechanism can be deployed and reused on different computer workstations or server arrays, providing a data processing environment for the evolution analysis of long-baseline hydroclimate.

[0098] For file system formatting standards, bad sector detection and management mechanisms, and underlying data read / write control protocols of computer-readable storage media, those skilled in the art can implement them in accordance with conventional storage device manufacturing and interface specifications. The basic data persistence principles are well-known technologies in the field and will not be elaborated here.

[0099] Specific application examples: To verify the actual effect of the hydrological and climatic element historical sequence reconstruction method based on multi-layer nested coupling provided by the present invention, this embodiment provides a set of real simulation and simulation experiments and effect comparison analysis.

[0100] Experimental Data Basis and Parameter Setting: This experiment selected a multi-source proxies dataset from a certain region spanning the past 1000 years (1001 to 2000 AD). (Includes tree-ring width and stalagmite oxygen isotope simulation data). Data from 1901 to 2000 (a total of 100 years) were selected as the instrument calibration period. Obtain the corresponding hydrological and climatic element observation sequences. (Precipitation); Data from 1001 to 1900 AD, a total of 900 years, was selected as the historical projection period. Reference time step The timeframe is set to one year. During the feature extraction stage, a decoupled hierarchical module is used to set the physical boundary periodic threshold. For a period of 10 years, high-frequency layer control sequences and low-frequency layer control sequences are separated. A continuous parameter function is constructed based on the Tikhonov regularization mechanism (with the penalty term expansion order set to third), and a discrete resolution is used. =0.05 radians generates a two-dimensional system lookup table Coherence threshold Set to 0.65 (corresponding to 95% red noise confidence level).

[0101] Comparison of scheme designs: To compare the results, the experiment introduced the traditional principal component multiple linear regression (TLR) model as the benchmark. The TLR model fits the proxy data to the observed sequence in the time domain without frequency domain decoupling, and the parameters are static constants over the entire time axis, lacking a logic gate variance injection mechanism.

[0102] Results Analysis and Technical Effects: See attached document Figure 3 During this period, due to environmental noise interference, the wavelet squared coherence between the multi-source substitute sequence and climate elements was locally lower than the set threshold. <0.65).

[0103] Comparison of variance stationarity: Depend on Figure 3 It can be seen that, due to the TLR model's use of global static regression and lack of consistency judgment, within the signal failure interval, the fluctuation amplitude of its reconstructed sequence (grey dashed line) decreases, becoming a curve fluctuating around the mean, exhibiting a phenomenon of reduced variance, and failing to reflect the interannual fluctuation characteristics of the natural climate system. This invention (black solid line) controls the sequence by triggering coherent logic gates. =0, stop inputting signals for that interval, and automatically inject an analog residual sequence. The reconstruction results maintained the baseline intercept term consistent with the calibration period. (like Figure 3 As shown on the vertical axis, this experiment set the precipitation baseline to 500 mm and retained the original high-frequency variance fluctuation characteristics, thus avoiding the problem of low variance in the historical reconstruction results.

[0104] See attached document Figure 4(The values ​​displayed on the vertical axis of "Calculation Time" in the figure have all been scaled by dividing by 10).

[0105] Comparison of model evaluation metrics: During the instrument calibration period, the Pearson correlation coefficient between the method of this invention and the observed data ( The value was 0.82, which is better than the 0.61 of the TLR model; Cross-validation error reduction value ( The value was increased from 0.35 in the TLR model to 0.68. This was achieved through separation. and And construct the dynamic regression coefficient function and This invention separates the differentiated climate response mechanisms at high and low frequencies, thereby improving the numerical values ​​of the model's evaluation indicators.

[0106] Comparison of computation time: In a historical extrapolation spanning 900 years (900 iterations), analytically calculating continuous high / low-frequency regression coefficient functions would take approximately 1450 milliseconds (corresponding to...). Figure 4 The reference value is 145.0. This invention employs discretization processing, searching a table using a two-dimensional system. Combined with a one-dimensional linear interpolation algorithm, the total calculation time is approximately 120 milliseconds (corresponding to...). Figure 4 The numerical value of this invention (12.0) shortens the computation time. The lookup table mechanism reduces the consumption of computing resources in climate extrapolation applications.

Claims

1. A method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling, characterized in that, Includes the following steps: The multi-source proxy data set and the hydrological and climatic element observation sequence are unified with the time resolution to generate a standard input sequence matrix, which divides the calibration period of the device and the historical projection period. Adaptive noise complete set empirical mode decomposition is performed on each data sequence in the standard input sequence matrix to separate and extract intrinsic mode functions and residual components, and then reassemble them to generate high-frequency layer control sequences and low-frequency layer control sequences. Calculate the cross-wavelet power spectrum of the high-frequency layer control sequence and the low-frequency layer control sequence, lock the dominant periodic frequency band, and extract the instantaneous phase difference sequence and wavelet square coherence sequence within the dominant periodic frequency band; A nested coupled regression benchmark model containing a benchmark intercept term is constructed. The dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function are calculated by combining the instantaneous phase difference sequence, and a two-dimensional system lookup table is generated. Based on the two-dimensional system lookup table, the historical high-frequency regression coefficients and historical low-frequency regression coefficients are obtained. A coherent logic gate control sequence is generated according to the wavelet squared coherence sequence. Combining the coherent logic gate control sequence, the historical high-frequency regression coefficients, the historical low-frequency regression coefficients, the high-frequency layer control sequence, the low-frequency layer control sequence, and the benchmark intercept term, the historical hydrological and climatic element reconstruction sequence is output.

2. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 1, characterized in that, The process of generating a standard input sequence matrix by unifying the time resolution of multi-source proxy data sets and hydrological and climatic element observation sequences, and dividing the instrument calibration period and historical extrapolation period, includes: Extract the intersection interval between the hydrological and climatic element observation sequence and the multi-source proxy data set on the time axis, and define the intersection interval as the instrument calibration period. Extract the historical time period in the multi-source proxy data set that is earlier than the instrument calibration period and is not covered by instrument observation data, and define the historical time period as the historical projection period. Identify the original sampling interval of each time series in the multi-source proxy dataset, and perform a time-domain forced alignment operation on the multi-source proxy dataset according to a set globally unified reference time step: when the original sampling interval is less than the reference time step, perform mean smoothing downsampling along the time axis according to the window span of the reference time step; When the original sampling interval is greater than the reference time step or there are non-uniform sampling distribution characteristics, numerical interpolation reconstruction is performed along the time axis. After completing the forced alignment operation in the time domain, the time series of each group in the multi-source substitute data set are truncated to a unified start and end boundary in the time dimension and combined to construct the standard input sequence matrix.

3. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 1, characterized in that, The step of performing adaptive noise-complete set empirical mode decomposition on each data sequence in the standard input sequence matrix, separating and extracting intrinsic mode functions and residual components, and recombining them to generate high-frequency layer control sequences and low-frequency layer control sequences includes: Multiple sets of independent zero-mean Gaussian white noise sequences are added to each data sequence in the standard input sequence matrix. The standard deviation of the zero-mean Gaussian white noise sequence is the product of the set noise amplitude coefficient and the standard deviation of the corresponding data sequence. Adaptive noise complete set empirical mode decomposition is performed on the noise-added sequence and the set mean is calculated at the time node. Multiple eigenmode functions with descending frequency characteristics and a residual component characterizing the long-term monotonic climate evolution trend are separated step by step. Calculate the average oscillation period of each intrinsic mode function on the time axis, and set the physical boundary period threshold to define the separation boundary between high-frequency fluctuation signals and low-frequency background signals; Iterate through all intrinsic mode functions, select intrinsic mode functions whose average oscillation period is less than the physical boundary period threshold, and sum and average the selected intrinsic mode functions across the dimensions of the cross-generational data set to obtain the high-frequency layer control sequence; The intrinsic mode functions whose average oscillation period is greater than or equal to the physical boundary period threshold are selected. The selected intrinsic mode functions are added to the corresponding residual components, and then a summation and averaging calculation is performed on the dimension of the cross-generational data set to obtain the low-frequency layer control sequence.

4. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 1, characterized in that, The calculation of the cross-wavelet power spectrum of the high-frequency layer control sequence and the low-frequency layer control sequence to lock the dominant periodic frequency band includes: Based on the complex Morlet wavelet function, continuous wavelet transform operations are performed on the high-frequency layer control sequence and the low-frequency layer control sequence respectively to obtain the high-frequency layer continuous wavelet transform spectrum and the low-frequency layer continuous wavelet transform spectrum. Calculate the complex product of the complex conjugate of the high-frequency layer continuous wavelet transform spectrum and the low-frequency layer continuous wavelet transform spectrum to generate a cross wavelet power spectrum that includes a scaling factor and a translation factor. Within the valid data region that passes the red noise confidence test in the cross wavelet power spectrum, find the energy extremum point with the largest absolute value, extract the scale factor corresponding to the energy extremum point as the dominant scale, and lock the dominant periodic frequency band.

5. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 4, characterized in that, After locking the dominant periodic frequency band, the step of extracting the instantaneous phase difference sequence and the wavelet squared coherence sequence within the dominant periodic frequency band includes: Extract the one-dimensional complex data sequence corresponding to the cross wavelet power spectrum on the dominant scale, perform imaginary and real part separation operations on the one-dimensional complex data sequence to obtain the one-dimensional imaginary part sequence and the one-dimensional real part sequence, calculate the arctangent function of the ratio of the value of the one-dimensional imaginary part sequence to the value of the one-dimensional real part sequence under the same translation factor, and obtain the instantaneous phase difference sequence. A two-dimensional smoothing operation covering translational smoothing along the time axis and weighted smoothing along the scale axis is performed on the absolute square of the cross wavelet power spectrum at the dominant scale. The smoothed result is divided by the product of the absolute square of the smoothed high-frequency continuous wavelet transform spectrum at the dominant scale and the absolute square of the smoothed low-frequency continuous wavelet transform spectrum at the dominant scale to obtain the wavelet square coherence sequence.

6. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 1, characterized in that, The construction of a nested coupled regression benchmark model including a benchmark intercept term, combined with the calculation of dynamic high-frequency regression coefficient functions and dynamic low-frequency regression coefficient functions based on the instantaneous phase difference sequence, generates a two-dimensional system lookup table, including: Using the corresponding sequences of the hydrological and climatic element observation sequences during the instrument calibration period as dependent variables, and the corresponding sequences of the high-frequency layer control sequences and the corresponding sequences of the low-frequency layer control sequences during the instrument calibration period as independent variables, a nested coupled regression benchmark model is established, which includes prior high-frequency regression coefficients, prior low-frequency regression coefficients, and the benchmark intercept term. Based on the corresponding sequence of the instantaneous phase difference sequence during the instrument calibration period, piecewise sliding regression is performed on the nested coupled regression benchmark model to extract the maximum and minimum values ​​of the local high-frequency regression coefficients and the maximum and minimum values ​​of the local low-frequency regression coefficients, which respectively constitute the high-frequency parameter modulation interval and the low-frequency parameter modulation interval. The dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function are expanded using truncated Fourier series of positive integer order from third to fifth order. A Tikhonov cost function containing a data fitting term and a regularization penalty term is constructed. The data fitting term is the sum of squared residuals between the corresponding sequence of the hydrological and climate element observation sequence during the instrument calibration period and the reconstructed sequence of the nested coupled regression benchmark model. The regularization penalty term is the square integral of the second derivative of the dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function as a function of the independent phase difference variable, multiplied by the regularization parameter. Using the prior high-frequency regression coefficients and the prior low-frequency regression coefficients as initial iterative estimates, and the high-frequency parameter modulation interval and the low-frequency parameter modulation interval as hard constraint boundaries for optimization, the minimum value of the Tikhonov cost function is solved to obtain a set of Fourier coefficients that satisfy the convergence condition. The continuous dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function are output, which together constitute a bivariate parameter modulation surface and are used to generate the two-dimensional system lookup table.

7. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 6, characterized in that, After constructing the bivariate parameter modulation surface, the generation process of converting the bivariate parameter modulation surface into the two-dimensional system lookup table and performing mesh discretization includes: Set the discrete resolution of the independent phase difference variable, divide the definition interval of the independent phase difference variable into equidistant grids based on the discrete resolution, and obtain the discrete phase difference sequence; Numerical calculations are performed by substituting the coordinates of each node in the discrete phase difference sequence into the dynamic high-frequency regression coefficient function and the dynamic low-frequency regression coefficient function, respectively, and calculating the discrete values ​​of the high-frequency regression coefficient and the low-frequency regression coefficient at the corresponding node positions. Perform spline interpolation smoothing operation, use cubic spline interpolation algorithm to smooth the connection of the discretized node data matrix, encapsulate the discrete phase difference sequence after spline interpolation smoothing, the set of discrete values ​​of the corresponding high-frequency regression coefficients and the set of discrete values ​​of the low-frequency regression coefficients, and construct and store the two-dimensional system lookup table.

8. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 7, characterized in that, The process of obtaining historical high-frequency regression coefficients and historical low-frequency regression coefficients based on the two-dimensional system lookup table includes: The system iterates through each discrete time node in the historical projection period in chronological order, extracts the current value of the corresponding sequence of the instantaneous phase difference sequence in the historical projection period as the input key value, and calls the two-dimensional system lookup table for retrieval. Traverse downwards in the index column of the two-dimensional system lookup table to lock the two adjacent discrete phase difference nodes that contain the input key value in the numerical range, and extract the high-frequency regression coefficient discrete value and low-frequency regression coefficient discrete value corresponding to the two discrete nodes. The distance weighting ratio is calculated based on the absolute value of the difference between the current input key value and the values ​​of the two discrete nodes. The extracted discrete values ​​are then weighted and averaged using a one-dimensional linear interpolation algorithm to output the high-frequency regression coefficient and the low-frequency regression coefficient of the historical period for the corresponding time node.

9. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 1, characterized in that, The process involves generating a coherent logic gate control sequence based on the wavelet squared coherence sequence, and combining this sequence with the historical high-frequency regression coefficients, the historical low-frequency regression coefficients, the high-frequency layer control sequence, the low-frequency layer control sequence, and the baseline intercept term to output a reconstructed sequence of historical hydrological and climatic elements, including: The wavelet squared coherence sequence is compared with the corresponding sequence in the historical extrapolation period and the set coherence threshold in time-by-time: when the corresponding sequence at the current time is greater than or equal to the coherence threshold, the coherence logic gate control sequence is set to 1. When the corresponding sequence at the current moment is less than the coherence threshold, the coherence logic gate control sequence is set to 0. Generate a simulated residual sequence that follows a normal distribution and has a mean of zero. The variance of the simulated residual sequence is equal to the average of the sum of squares of the residuals of the nested coupled regression benchmark model calculated during the instrument calibration period. Perform numerical algebraic calculations for integrated nested coupled extrapolation: multiply the historical high-frequency regression coefficients and the corresponding sequences of the high-frequency layer control sequences within the historical extrapolation period by the product of the historical low-frequency regression coefficients and the corresponding sequences of the low-frequency layer control sequences within the historical extrapolation period, multiply the sum of the products by the coherent logic gate control sequence, add the benchmark intercept term, and finally add the product of the value 1, the difference between the coherent logic gate control sequence, and the simulated residual sequence to output the reconstructed sequence of hydrological and climatic elements for the historical period.

10. The method for reconstructing historical sequences of hydrological and climatic elements based on multi-layer nested coupling according to claim 9, characterized in that, After outputting the reconstructed sequence of historical hydrological and climatic elements, the process also includes the splicing and output of the full-band sequence: An overlapping time window is set between the historical projection period and the instrument calibration period, and the mean deviation between the reconstructed sequence of hydrological and climatic elements in the historical period and the corresponding sequence of the observed sequence of hydrological and climatic elements in the instrument calibration period is calculated in the overlapping time window. Using the calculated mean deviation, translation compensation is performed on the historical hydro-climate element reconstruction sequence to obtain the bias-corrected historical reconstruction sequence. The historical reconstructed sequence after deviation correction is merged with the corresponding sequence of the hydrological and climate element observation sequence within the instrument calibration period in the forward order of time variables to generate a full-band hydrological and climate reconstruction sequence. A weighted moving average algorithm is adopted, and the sliding step size is set to 3 to 7 time nodes. The low-pass filtering principle is applied to the data segments on both sides of the time boundary point of the full-band hydrological and climate reconstruction sequence. The processed sequence is then encapsulated into a standard data structure for output.

Images