Air quality index prediction method based on three-level progressive residual decomposition

Abstract

The invention discloses an air quality index prediction method based on three-level progressive residual decomposition, and relates to the technical field of air quality monitoring. Comprising the following steps: step 1, collecting air quality monitoring data of a target monitoring station, and arranging air quality index values at continuous monitoring moments according to a time sequence to construct an original air quality index sequence; 2, performing resonance sparse signal decomposition on the original air quality index sequence, and subtracting a first-stage trend component from the original air quality index sequence point by point to obtain a first-stage residual sequence; 3, performing pole symmetric mode decomposition on the first-stage residual error sequence, and adding the rest intrinsic mode components and the final residual error sequence point by point to serve as third-stage random components; and 4, respectively predicting the first-stage trend component, the second-stage periodic component and the third-stage random component. According to the method, accurate matching of the prediction model and the component characteristics is realized, and the precision and stability of air quality index prediction are effectively improved.

Description

Technical Field

[0001] This invention relates to the field of air quality monitoring technology, and in particular to an air quality index prediction method based on three-level progressive residual decomposition. Background Technology

[0002] In recent years, researchers have gradually recognized that air quality index (AQI) sequences contain multiple components at different time scales, including long-term seasonal trends, periodic fluctuations related to diurnal variations and traffic cycles, and random disturbances caused by sudden pollution events. These components differ significantly in frequency characteristics and variation patterns, making it difficult to adapt a single prediction model to the characteristics of all types of components simultaneously. Therefore, combinatorial prediction methods based on signal decomposition have attracted widespread attention. The core idea is to first decompose the original sequence into several subsequences with different characteristics, then establish a prediction model for each subsequence, and finally reconstruct the prediction results of each subsequence into the final predicted value.

[0003] Common signal decomposition methods include wavelet decomposition, empirical mode decomposition (EMD) and its improvements, and variational mode decomposition (VMD). Wavelet decomposition decomposes a signal into different frequency sub-bands based on predefined wavelet basis functions. The quality of the decomposition result depends on the choice of wavelet basis functions, and fixed frequency band division methods are difficult to adaptively match the actual frequency characteristics of the signal. Empirical mode decomposition (EMD) is a data-driven adaptive decomposition method that can decompose a signal into several intrinsic mode functions (EMFs), but it suffers from two prominent problems: mode aliasing and endpoint effects. Mode aliasing refers to the mixing of components from different time scales within the same EEM, leading to unclear physical meaning of the decomposition result. Endpoint effects refer to the propagation of decomposition errors from the endpoints to the middle of the sequence due to insufficient data support for envelope fitting at signal boundaries. Variational mode decomposition decomposes a signal into a predetermined number of modal components by constructing a variational optimization problem, which alleviates the mode aliasing problem to some extent. However, the number of decomposition modes needs to be specified in advance, and for complex signals containing both continuous oscillation components and transient impulse components, a single decomposition mechanism is insufficient to effectively separate the two types of components. Summary of the Invention

[0004] In view of this, the present invention provides an air quality index prediction method based on three-level progressive residual decomposition. By combining resonant sparse signal decomposition with pole-symmetric mode decomposition, a three-level progressive decomposition framework is constructed. Resonant sparse signal decomposition is used to accurately extract trend components based on the sparsity differences in the wavelet domain of different quality factors. The pole-symmetric extension strategy of pole-symmetric mode decomposition is used to effectively suppress the endpoint effect and separate periodic components from random components. Long short-term memory networks, gated recurrent unit networks, and one-dimensional convolutional neural networks are used for differentiated prediction based on the different characteristics of each component. This achieves accurate matching between the prediction model and the component characteristics, effectively improving the accuracy and stability of air quality index prediction.

[0005] The technical solution adopted in this invention is as follows:

[0006] The air quality index prediction method based on three-level progressive residual decomposition includes the following steps:

[0007] Step 1: Collect air quality monitoring data from the target monitoring station, preprocess the air quality monitoring data, calculate the air quality index value at each monitoring time, and arrange the air quality index values ​​of continuous monitoring times in chronological order to construct the original air quality index sequence.

[0008] Step 2: Perform resonance sparse signal decomposition on the original air quality index sequence. By constructing a high-quality factor tunable wavelet transform filter bank and a low-quality factor tunable wavelet transform filter bank, perform dual-channel wavelet decomposition on the original air quality index sequence. Use the sparse coefficient optimization solution based on split augmented Lagrange contraction to obtain the high resonance component. Use the high resonance component as the first-level trend component. Subtract the first-level trend component from the original air quality index sequence point by point to obtain the first-level residual sequence.

[0009] Step 3: Perform pole-symmetric mode decomposition on the first-level residual sequence. Separate multiple intrinsic mode components in turn through pole-symmetric extension and envelope mean extraction. Add the first two intrinsic mode components point by point to obtain the second-level periodic components. Add the remaining intrinsic mode components to the final residual sequence point by point to obtain the third-level random components.

[0010] Step 4: Predict the first-level trend component, the second-level periodic component, and the third-level random component respectively, and add the prediction results of each component point by point to obtain the air quality index prediction result.

[0011] Furthermore, in step one, the air quality monitoring data includes the concentrations of fine particulate matter, inhalable particulate matter, sulfur dioxide, nitrogen dioxide, carbon monoxide, and ozone; the preprocessing includes linear interpolation of missing values ​​and median filtering of outliers.

[0012] Furthermore, in step two, the construction process of the high-quality factor tunable wavelet transform filter bank is as follows: The high-quality factor is set to 4, and the high redundancy parameter is set to 3. Based on the high-quality factor, the passband cutoff frequency of the high-pass filter is calculated to be at position 0.8 of the normalized frequency axis, and the passband cutoff frequency of the low-pass filter is calculated to be at position 0.2 of the normalized frequency axis. A high-pass finite impulse response (FIR) filter with a second-order zero at the zero position of the frequency axis and a corresponding low-pass FIR filter are constructed, ensuring that the high-pass FIR filter and the low-pass FIR filter meet the complete reconstruction condition. The high-pass FIR filter and the low-pass FIR filter are cascaded according to the decomposition level determined by the high redundancy parameter to form a high-quality factor tunable wavelet transform filter bank.

[0013] Furthermore, the construction process of the low-quality factor tunable wavelet transform filter bank is as follows: set the low-quality factor value to 1, set the low redundancy parameter value to 3, calculate the passband cutoff frequency position of both the high-pass filter and the low-pass filter based on the low-quality factor value, and construct the high-pass finite impulse response filter and the low-pass finite impulse response filter in the same way as the high-quality factor tunable wavelet transform filter bank, and cascade them to form the low-quality factor tunable wavelet transform filter bank.

[0014] Furthermore, in step two, the execution process of dual-channel wavelet decomposition is as follows: the original air quality index sequence is passed sequentially through low-pass and high-pass filters of the filter bank. After each filtering stage, it is downsampled according to the factor determined by the redundancy parameter. The decomposition is carried out step by step until the preset decomposition level is reached, and the low-pass sub-band coefficient sequences and high-pass sub-band coefficient sequences of each stage are output. All sub-band coefficient sequences obtained by decomposition based on the high-quality factor tunable wavelet transform filter bank are combined into a high-resonance wavelet coefficient vector, and all sub-band coefficient sequences obtained by decomposition based on the low-quality factor tunable wavelet transform filter bank are combined into a low-resonance wavelet coefficient vector.

[0015] Furthermore, the optimization solution for sparse coefficients based on split augmented Lagrange contraction includes: initializing the high-resonance sparse coefficient vector, low-resonance sparse coefficient vector, first auxiliary variable vector, second auxiliary variable vector, first Lagrange multiplier vector, and second Lagrange multiplier vector as all zero vectors; setting the penalty parameter, high-resonance regularization parameter, low-resonance regularization parameter, convergence threshold, and maximum number of iterations; performing iterative optimization until convergence; and performing inverse wavelet reconstruction on the converged high-resonance sparse coefficient vector to obtain the high-resonance component.

[0016] Furthermore, each iteration of the iterative optimization performs the following operations in sequence: subtract the low-resonance component reconstructed based on the current low-resonance sparse coefficient vector from the original air quality index sequence and then perform forward wavelet decomposition to update the high-resonance sparse coefficient vector; subtract the high-resonance component reconstructed based on the current high-resonance sparse coefficient vector from the original air quality index sequence and then perform forward wavelet decomposition to update the low-resonance sparse coefficient vector; perform soft-threshold shrinkage updates on the first auxiliary variable vector and the second auxiliary variable vector respectively; update the first Lagrange multiplier vector and the second Lagrange multiplier vector; the soft-threshold shrinkage update is performed as follows: when the absolute value of an element is greater than the quotient of the regularization parameter and the penalty parameter, the element value is subtracted from the product of the quotient of the regularization parameter and the penalty parameter and the element sign to obtain the shrunken element value; when the absolute value of an element is less than or equal to the quotient of the regularization parameter and the penalty parameter, the shrunken element value is set to zero.

[0017] Furthermore, in step three, the pole-symmetric mode decomposition includes: performing extreme point detection on the first-level residual sequence, marking sampling points whose values ​​are simultaneously greater than the values ​​of the adjacent preceding and following sampling points as maxima, and marking sampling points whose values ​​are simultaneously less than the values ​​of the adjacent preceding and following sampling points as minima; forming a set of maxima in chronological order from all maxima and minima in chronological order from all minima; and performing pole-symmetric extension processing on the set of maxima and the set of minima, wherein the pole-symmetric extension processing is performed by mirroring the extreme points at the boundary with the first and last sampling points of the first-level residual sequence as the center of symmetry.

[0018] Furthermore, pole-symmetric mode decomposition also includes: generating an upper envelope based on the extended set of maxima using cubic spline interpolation, and generating a lower envelope based on the extended set of minima using cubic spline interpolation; calculating the arithmetic mean of the upper and lower envelopes at each sampling time, and forming an envelope mean line by arranging all arithmetic means in chronological order; subtracting the envelope mean line from the current input sequence point by point to obtain intrinsic mode candidate components; when an intrinsic mode candidate component satisfies the intrinsic mode function condition, outputting the intrinsic mode candidate component as the intrinsic mode component. The modal function condition is that the absolute value of the difference between the number of extreme points and the number of zero crossings is less than or equal to 1. The first-level residual sequence is subtracted point by point from the current intrinsic modal component to obtain the intermediate residual sequence. The extreme point detection, extreme point symmetric extension processing, upper and lower envelope generation, envelope mean line calculation and intrinsic modal candidate component extraction are repeatedly performed on the intermediate residual sequence until the number of extreme points in the intermediate residual sequence is less than 3. The first two intrinsic modal components are added point by point to obtain the second-level periodic component, and the remaining intrinsic modal components are added point by point to the final residual sequence to obtain the third-level random component.

[0019] Furthermore, in step four, a long short-term memory network is used to predict the first-level trend component to obtain a trend prediction sequence, a gated recurrent unit network is used to predict the second-level periodic component to obtain a periodic prediction sequence, and a one-dimensional convolutional neural network is used to predict the third-level random component to obtain a random prediction sequence.

[0020] By adopting the above technical solution, the present invention has produced the following beneficial effects:

[0021] This invention employs resonant sparse signal decomposition as the first-level decomposition method. By constructing high-quality factor tunable wavelet transform filter banks and low-quality factor tunable wavelet transform filter banks, and utilizing the differences in sparsity performance of trend components and transient components in wavelet domains with different quality factors, it achieves accurate extraction of long-period trend components. Compared to traditional wavelet decomposition and empirical mode decomposition methods, resonant sparse signal decomposition can adaptively separate trend components with continuous oscillation characteristics from disturbance components with transient impact characteristics, avoiding the subjectivity of pre-specifying the number of decomposition modes, and the decomposition results have more explicit physical meaning.

[0022] This invention employs pole-symmetric mode decomposition to perform a secondary decomposition of the first-level residual sequence. By utilizing a pole-symmetric extension strategy, it effectively suppresses the endpoint effect problem commonly found in traditional empirical mode decomposition methods. Suppressing the endpoint effect prevents the decomposition error from spreading from the sequence boundaries towards the center, ensuring the stability and accuracy of the decomposition results throughout the entire sequence. Pole-symmetric mode decomposition can further separate the first-level residual sequence into components with periodic oscillation characteristics and irregular random fluctuation components, laying the foundation for subsequent hierarchical prediction.

[0023] This invention employs a progressive residual decomposition strategy, extending the traditional single-level decomposition to a three-level progressive decomposition. Each level of decomposition focuses on separating a specific type of component: the first level extracts trend components, the second level extracts periodic components, and the third level separates random components. This progressive decomposition strategy makes the characteristics of each component purer, avoids mutual interference between different types of components, and facilitates subsequent prediction models to fully learn the intrinsic patterns of each component.

[0024] This invention employs differentiated prediction models for the three components based on their different characteristics. The trend component uses a Long Short-Term Memory (LSTM) network to capture long-range dependencies; the periodic component uses a gated recurrent unit (GRU) network to model periodic changes; and the random component uses a one-dimensional convolutional neural network to extract local fluctuation patterns. This hierarchical prediction strategy fully leverages the advantages of various neural network models, achieving precise matching between the prediction model and the component characteristics, and outperforming traditional methods using a single model in overall prediction performance. Attached Figure Description

[0025] Figure 1This is a diagram illustrating the effect of the three-level progressive residual decomposition in an embodiment of the present invention.

[0026] Figure 2 This is a frequency response characteristic diagram of the tunable quality factor wavelet transform filter bank in the embodiment of the present invention.

[0027] Figure 3 This is a schematic diagram of envelope fitting during the pole-symmetric mode decomposition process in an embodiment of the present invention;

[0028] Figure 4 This is a characteristic curve of the soft threshold shrinkage function in an embodiment of the present invention;

[0029] Figure 5 This is a comparison chart of the stratified prediction and the air quality index reconstruction prediction results in this embodiment of the invention. Detailed Implementation

[0030] Any feature disclosed in this specification, unless otherwise stated, may be replaced by other equivalent or similar features. That is, unless otherwise stated, each feature is merely one example of a series of equivalent or similar features.

[0031] The air quality index prediction method based on three-level progressive residual decomposition includes the following steps:

[0032] Step 1: Collect air quality monitoring data from the target monitoring station, preprocess the air quality monitoring data, calculate the air quality index value at each monitoring time, and arrange the air quality index values ​​of continuous monitoring times in chronological order to construct the original air quality index sequence.

[0033] Step 2: Perform resonance sparse signal decomposition on the original air quality index sequence. By constructing a high-quality factor tunable wavelet transform filter bank and a low-quality factor tunable wavelet transform filter bank, perform dual-channel wavelet decomposition on the original air quality index sequence. Use the sparse coefficient optimization solution based on split augmented Lagrange contraction to obtain the high resonance component. Use the high resonance component as the first-level trend component. Subtract the first-level trend component from the original air quality index sequence point by point to obtain the first-level residual sequence.

[0034] Step 3: Perform pole-symmetric mode decomposition on the first-level residual sequence. Separate multiple intrinsic mode components in turn through pole-symmetric extension and envelope mean extraction. Add the first two intrinsic mode components point by point to obtain the second-level periodic components. Add the remaining intrinsic mode components to the final residual sequence point by point to obtain the third-level random components.

[0035] Step 4: Predict the first-level trend component, the second-level periodic component, and the third-level random component respectively, and add the prediction results of each component point by point to obtain the air quality index prediction result.

[0036] Accurate prediction of the Air Quality Index (AQI) relies on high-quality raw monitoring data. In this embodiment, air quality monitoring data from target monitoring stations is used as the data source for the prediction system. The air quality monitoring data is acquired through specialized sensor equipment deployed at the monitoring stations, which continuously samples and measures the concentration of various pollutants in the atmosphere at fixed time intervals.

[0037] Air quality monitoring data includes the concentrations of fine particulate matter (PM2.5), inhalable particulate matter (PM10), sulfur dioxide (SO2), nitrogen dioxide (NO2), carbon monoxide (CO2), and ozone. PM2.5 concentration refers to the mass concentration of particles with an aerodynamic equivalent diameter (AEM) of 2.5 micrometers or less, typically recorded in micrograms per cubic meter (µg / m³). Inhalable particulate matter concentration refers to the mass concentration of particles with an AEM of 10 micrometers or less. SO2, NO2, and ozone concentrations are all recorded in µg / m³, while CO2 concentration is recorded in milligrams per cubic meter (mg / m³). In a specific implementation scenario, monitoring stations collect concentration data for these six pollutants hourly, continuously for 720 hours to form a 30-day monitoring dataset. Due to potential communication interruptions, equipment malfunctions, or environmental interference during sensor operation, the collected air quality monitoring data often contains missing and outlier values. Missing values ​​indicate that concentration data for certain monitoring times were not recorded, while outlier values ​​indicate that the concentration values ​​deviate significantly from the normal range. If these data quality issues are not addressed, they will directly affect the accuracy of subsequent air quality index calculations and the reliability of prediction models.

[0038] Air quality monitoring data undergoes preprocessing to eliminate the impact of data quality issues. Preprocessing includes linear interpolation for missing values ​​and median filtering for outliers. The specific implementation of linear interpolation for missing values ​​is as follows: when a missing concentration data point is detected at a certain monitoring time, the nearest valid data point before and after the missing point is found. Assuming the time index of the missing point is... The index of the time of the previous valid data point is The corresponding concentration value The index of the time of the next valid data point is The corresponding concentration value The concentration value of the missing point is then calculated using linear interpolation, and the calculation method is as follows: ,in This represents the interpolated concentration value. This linear interpolation method is based on the characteristic that pollutant concentrations change relatively slowly over a short period of time, and can effectively recover reasonable estimates of missing data.

[0039] The specific execution method of median outlier filtering is as follows: For the concentration data at each monitoring time, a sliding window centered on that time is selected, with a window length of 5 sampling points. All concentration values ​​within the window are extracted and sorted by value. The value in the middle position after sorting is taken as the filtered output value. Median filtering can effectively suppress the influence of impulse outliers while maintaining the overall trend of concentration data. Compared with mean filtering, median filtering has a stronger ability to suppress isolated extreme outliers and will not cause significant shifts in the filtered output due to the presence of a single outlier. In an optional implementation, the sliding window length can be adjusted according to the data sampling frequency and the distribution characteristics of outliers, for example, set to 3 or 7 sampling points.

[0040] After preprocessing, the Air Quality Index (AQI) values ​​for each monitoring time are calculated according to the AQI technical specifications. The AQI is an evaluation index that integrates the concentrations of multiple pollutants into a single numerical value, providing a direct reflection of air quality. The calculation process first calculates the sub-index for each pollutant, then takes the maximum value among all sub-indices as the AQI value for that monitoring time. For calculating the sub-index for a single pollutant, the concentration range of that pollutant needs to be determined by referring to the concentration limit classification table. Taking fine particulate matter (PM2.5) as an example, when the 24-hour average concentration is between 0 and 35 micrograms per cubic meter, the corresponding sub-index range is 0 to 50; when the concentration is between 35 and 75 micrograms per cubic meter, the corresponding sub-index range is 50 to 100; and so on until the highest concentration range. After determining the concentration range, the specific sub-index value is calculated using linear interpolation. The calculation method is as follows: ,in Indicates pollutants The sub-index, Indicates pollutants The measured concentration value, and These represent the lower and upper limits of the concentration range, respectively. and These represent the lower and upper limits of the sub-index interval, respectively.

[0041] After calculating the sub-indices for each of the six pollutants, the maximum value is taken as the air quality index value at that monitoring time. ,in This indicates the air quality index value. to These represent sub-indices for fine particulate matter, inhalable particulate matter, sulfur dioxide, nitrogen dioxide, carbon monoxide, and ozone, respectively. This method of taking the maximum value ensures that the air quality index reflects the most prominent air pollution problems at present.

[0042] The air quality index values ​​from continuous monitoring times are arranged in chronological order to construct the original air quality index sequence. Let there be a total of... If there are multiple monitoring times, the original air quality index sequence can be represented as: ,in Indicates the first Air quality index values ​​at each monitoring time. This represents the total length of the sequence. In this embodiment, data is continuously collected for 720 hours at 1-hour sampling intervals. The value is 720. The original air quality index series serves as the input data for the subsequent three-stage progressive residual decomposition.

[0043] refer to Figure 1 The top of the image shows the original air quality index sequence. The horizontal axis represents time in hours, and the vertical axis represents the air quality index (AQI) value. The original sequence contains 480 sampling points, corresponding to hourly monitoring data over 20 consecutive days. The waveform reveals a clear composite characteristic in the original sequence, containing both long-term, gradual transitions and short-term, rapid oscillations, with the AQI value fluctuating between 20 and 180. The second row in the figure represents the first-order trend component. This component was extracted from the original sequence using the resonant sparse signal decomposition method. The trend component has a relatively smooth waveform, mainly reflecting the variation of the air quality index over time scales of several days to several weeks, with its values ​​concentrated between 60 and 140. The waveform of the trend component exhibits a periodic variation superimposed with a linear increase, similar to a sine wave, which is consistent with the physical mechanism by which air pollutant concentrations are affected by seasonal factors and long-term emission trends. The third row in the figure represents the second-order periodic component. This component is obtained by performing pole-symmetric mode decomposition on the first-level residual sequence and extracting and adding the first two eigenmode components. The periodic component oscillates around zero with an amplitude range of approximately ±40, exhibiting regular periodic fluctuations. The oscillation frequency of the waveform reveals two main components: a diurnal period and a semi-diurnal period, corresponding to the morning and evening peak variations in traffic flow and the periodic influence of human activity's diurnal rhythm on air pollutant concentrations. The bottom of the figure represents the third-level random component. This component is formed by adding the remaining intrinsic mode components generated by pole-symmetric mode decomposition to the final residual sequence. The random component also fluctuates around zero, but its changes lack a clear regularity, with an amplitude range of approximately ±20. This component mainly reflects short-term fluctuations in the air quality index caused by unpredictable factors such as sudden pollution events and random meteorological disturbances. Through three-level progressive residual decomposition, the original air quality index sequence was successfully separated into three components with different time-scale characteristics and physical meanings, satisfying... The decomposition relationship.

[0044] Air Quality Index (AQI) series, as a typical environmental monitoring time series, contains various components with different characteristics. From a signal analysis perspective, AQI series simultaneously exhibit long-term, slowly varying trend components and short-term transient impact components. Trend components reflect the evolution of air pollution over a longer time scale, typically related to seasonal variations and changes in regional total emissions, and are characterized by long duration and multiple oscillation periods. Transient impact components reflect short-term pollution fluctuations caused by sudden pollution events and extreme weather conditions, and are characterized by short duration and drastic changes. Resonant sparse signal decomposition is a signal processing technique that can effectively separate these two types of components. Its core idea is to utilize the sparsity differences of different components in the wavelet domain to achieve component separation. Due to its multiple oscillation periods, trend components exhibit good sparsity in the high-quality factor wavelet transform domain, meaning they can be represented with fewer non-zero coefficients; transient impact components, due to their short duration, exhibit good sparsity in the low-quality factor wavelet transform domain. By constructing two sets of wavelet transform filter banks with different quality factors and applying sparsity constraints for optimization, the two types of components can be effectively separated.

[0045] refer to Figure 2 In the graph, the horizontal axis represents the normalized frequency, ranging from 0 to 1, where a normalized frequency of 1 corresponds to half the sampling frequency. The vertical axis represents the amplitude response, also ranging from 0 to 1, indicating the degree to which the filter passes different frequency components. The upper left corner of the graph shows the frequency response curve of the high-pass filter in a high-quality factor tunable wavelet transform filter bank. (High-quality factor) Setting it to 4 corresponds to the high-pass filter's passband cutoff frequency being located at 0.8 on the normalized frequency axis. The curve shows that when the normalized frequency is below 0.7, the filter's amplitude response is close to zero, indicating effective suppression of low-frequency components; when the normalized frequency is between 0.7 and 0.8, the amplitude response smoothly transitions from zero to 1, forming a transition band; when the normalized frequency is above 0.8, the amplitude response remains at 1, indicating that high-frequency components pass completely. The curve has a second-order zero at the zero position on the frequency axis, ensuring complete suppression of the DC component. The upper right of the figure shows the frequency response curve of the low-pass filter in the high-quality factor tunable wavelet transform filter bank, with the passband cutoff frequency located at 0.2 on the normalized frequency axis. The low-pass filter and high-pass filter satisfy the complete reconstruction condition; their combined use can achieve lossless decomposition and reconstruction of the signal. The lower left of the figure shows the frequency response curve of the high-pass filter in the low-quality factor tunable wavelet transform filter bank. Low quality factor With a quality factor of 1, the passband cutoff frequency of the high-pass filter is located at 0.5 on the normalized frequency axis. Compared to the high-quality factor configuration, the low-quality factor filter has a wider transition band and relatively lower frequency selectivity. This characteristic makes it more suitable for matching transient impulse components with short durations and wide spectral distributions. The lower right corner of the figure shows the frequency response curve of the low-pass filter of the low-quality factor tunable wavelet transform filter bank. Its passband cutoff frequency is also located at 0.5 on the normalized frequency axis, forming a symmetrical frequency division with the corresponding high-pass filter. The different quality factor configurations of the two filter banks cause the trend components and transient components in the original signal to exhibit different sparsity characteristics in their respective wavelet transform domains, providing a basis for subsequent component separation based on sparsity differences.

[0046] Resonant sparse signal decomposition is performed on the original air quality index sequence to extract the first-order trend component and generate the first-order residual sequence. The first step of resonant sparse signal decomposition is to construct a high-quality factor tunable wavelet transform filter bank and a low-quality factor tunable wavelet transform filter bank.

[0047] Tunable quality factor wavelet transform (TMB) is a wavelet transform method that allows for flexible setting of the oscillation characteristics of wavelet functions. Unlike traditional wavelet transforms that use fixed wavelet basis functions, TMB controls the oscillation degree of wavelet basis functions by adjusting the quality factor parameter. The quality factor is defined as the ratio of the center frequency to the bandwidth of the wavelet function. A higher quality factor results in more oscillation periods, making it more suitable for representing components with continuous oscillation characteristics; a lower quality factor results in fewer oscillation periods, making it more suitable for representing transient pulse components. This tunability allows the same signal to exhibit different sparsity characteristics in the wavelet transform domain with different quality factors, providing a theoretical basis for component separation based on sparsity differences.

[0048] The construction process of the high-quality factor tunable wavelet transform filter bank is as follows. The high-quality factor is set to 4, which represents the ratio of the wavelet function's center frequency to its bandwidth. A value of 4 means the wavelet function has a relatively large oscillation period, suitable for matching the trend components in the air quality index sequence. The high redundancy parameter is set to 3. The redundancy parameter controls the overcompleteness of the wavelet transform. A value of 3 means the total number of coefficients after the transform is approximately three times the length of the original signal. Appropriate redundancy can improve the stability of signal reconstruction and the accuracy of decomposition.

[0049] The frequency response characteristics of the filters are calculated based on the quality factor. The passband cutoff frequency of the high-pass filter is set at 0.8 of the normalized frequency axis, and the passband cutoff frequency of the low-pass filter is set at 0.2 of the normalized frequency axis. The normalized frequency axis maps the actual frequency to a range of 0 to 1, where 1 corresponds to half the sampling frequency. Setting the high-pass filter cutoff frequency to 0.8 means that the filter mainly passes signal components with frequencies higher than 0.4 times the sampling frequency, and setting the low-pass filter cutoff frequency to 0.2 means that the filter mainly passes signal components with frequencies lower than 0.1 times the sampling frequency. This frequency division method forms a transition region between the high-pass and low-pass bands. The width of the transition region is directly related to the quality factor; a higher quality factor corresponds to a narrower transition band, making the frequency band division more refined.

[0050] A high-pass finite impulse response (FIR) filter with a second-order zero at the zero position on the frequency axis is constructed. The purpose of setting a zero at the zero position on the frequency axis in the high-pass filter is to ensure complete suppression of the DC component. A second-order zero means that the frequency response of the high-pass filter has not only zero amplitude at zero frequency, but also zero first derivative. This design makes the suppression of low-frequency components smoother and avoids drastic changes in frequency response in the extremely low-frequency region. The design of the FIR filter uses the window function method. First, the impulse response of the ideal high-pass filter is determined, and then truncated by multiplying it by a Hanning window function. The filter order is set to 64 to achieve a balance between frequency selectivity and computational complexity.

[0051] A corresponding low-pass finite impulse response (FIR) filter is constructed such that the high-pass FIR filter and the low-pass FIR filter satisfy the complete reconstruction condition. The complete reconstruction condition means that the original signal can be accurately recovered after the signal is decomposed by the filter bank and then reconstructed by the corresponding synthetic filter bank. Satisfying the complete reconstruction condition requires a specific constraint relationship between the frequency responses of the analytical filter and the synthetic filter. In this embodiment, the low-pass filter is designed based on the high-pass filter obtained through frequency inversion and sign alternation, ensuring that the two form an orthogonal mirror filter pair, thereby satisfying the complete reconstruction condition.

[0052] A high-pass finite impulse response (FIR) filter and a low-pass FIR filter are cascaded to form a high-quality factor tunable wavelet transform (HFT) filter bank, based on the decomposition levels determined by the high redundancy parameter. The decomposition levels correspond to the redundancy parameter and the quality factor; with a high redundancy parameter of 3 and a quality factor of 4, the decomposition levels are set to 10. The filter bank employs an iterative filter bank structure. The first stage divides the input signal into low-frequency and high-frequency subbands, and subsequent stages further decompose the low-frequency subband, forming a tree-like multi-resolution decomposition structure.

[0053] The construction process of a low-quality-factor tunable wavelet transform filter bank is similar to that of a high-quality-factor tunable wavelet transform filter bank, but different parameter configurations are used to match the characteristics of transient impulse components. The low-quality-factor is set to 1, indicating that the center frequency and bandwidth of the wavelet function are equal, and the wavelet function contains only a small number of oscillation periods, suitable for matching short-duration transient components. The low redundancy parameter is set to 3, maintaining the same level of redundancy as the high-quality-factor tunable wavelet transform filter bank.

[0054] The frequency response characteristics of the filter are calculated based on the low quality factor value. Since the low quality factor is 1, the passband cutoff frequency of both the high-pass and low-pass filters is set to 0.5 of the normalized frequency axis. This symmetrical frequency division ensures that the high-pass and low-pass bands each occupy half of the frequency axis, resulting in a wider transition band compared to the high quality factor configuration, making it suitable for handling transient components with a broad frequency response. A high-pass finite impulse response (FECT) filter and a low-pass finite impulse response (FECT) filter are constructed and cascaded in the same manner as the high quality factor tunable wavelet transform filter bank to form a low quality factor tunable wavelet transform filter bank.

[0055] In optional implementations, the high-quality factor value can be selected from the range of 3 to 6, the low-quality factor value from the range of 1 to 2, and the redundancy parameter value from the range of 2 to 4. The specific values ​​can be adjusted according to the actual characteristics of the air quality index sequence. When the periodicity of the trend component in the sequence is more obvious, the high-quality factor value can be appropriately increased; when the duration of the transient impact component is shorter, the low-quality factor value can be appropriately decreased.

[0056] After constructing the filter bank, the original air quality index sequence is decomposed using a dual-channel wavelet decomposition. The dual-channel wavelet decomposition process is as follows: the original air quality index sequence is passed sequentially through low-pass and high-pass filters of the filter bank. After each filtering stage, it is downsampled according to the factor determined by the redundancy parameter. This process is repeated step by step until the preset number of decomposition levels is reached, outputting the low-pass subband coefficient sequences and high-pass subband coefficient sequences of each stage.

[0057] Specifically, let the original air quality index sequence be... It is then input into the first stage of a high-quality factor-tunable wavelet transform filter bank, and passed through a low-pass filter. After filtering, downsampling is performed to obtain the first-stage low-pass subband coefficients, which are then passed through a high-pass filter. After filtering, downsampling is performed to obtain the first-level high-pass subband coefficients. The downsampling factor is determined by the redundancy parameter; when the redundancy parameter is 3, the downsampling factor is 2. The first-level low-pass subband coefficients are then input into the second-level filter for decomposition, and this process is iterated until all 10 levels of decomposition are completed. After decomposition, the high-pass subband coefficient sequences of each level and the last-level low-pass subband coefficient sequence are arranged sequentially to form a high-resonance wavelet coefficient vector, denoted as [vector]. .

[0058] The original air quality index sequence is decomposed using the same method through a low-quality-factor tunable wavelet transform filter bank to obtain a low-resonance wavelet coefficient vector, denoted as . Because the two filter banks have different quality factors, the trend component in the original air quality index sequence is reflected in the high-resonance wavelet coefficient vector. The coefficients exhibit a sparse distribution, meaning that most coefficients have small absolute values, while only a small number have large absolute values; the transient impulse component is present in the low-resonance wavelet coefficient vector. It exhibits a sparse distribution.

[0059] High-resonance components are obtained through sparse coefficient optimization based on split-augmented Lagrange contraction. The goal of sparse coefficient optimization is to find a solution that maximizes the sparsity of both the high-resonance and low-resonance wavelet coefficient vectors, while satisfying signal reconstruction constraints. Sparsity is achieved through the coefficient vectors... Norms are used for measurement. The norm is defined as the sum of the absolute values ​​of the elements in a vector. The formal formulation of the optimization problem is: finding a vector of highly resonant sparse coefficients. and low resonance sparse coefficient vector , making Minimize while satisfying constraints ,in express Norm operations This represents the high resonance regularization parameter. This represents the low resonance regularization parameter. This represents the inverse transform matrix corresponding to the high-quality factor tunable wavelet transform filter bank. This represents the inverse transform matrix corresponding to the low-quality-factor tunable wavelet transform filter bank.

[0060] The high-resonance regularization parameter and the low-resonance regularization parameter control the relative proportions of the two types of components in the decomposition results. When the high-resonance regularization parameter is relatively small, the optimization process tends to retain more high-resonance components; when the low-resonance regularization parameter is relatively small, the optimization process tends to retain more low-resonance components. In this embodiment, since the goal is to extract the trend component as the high-resonance component, the high-resonance regularization parameter is set to 0.1 and the low-resonance regularization parameter is set to 0.5, so that the optimization process applies a weaker sparsity constraint to the high-resonance component, thereby preserving the main structure of the trend component.

[0061] The split-augmented Lagrange contraction method transforms the original optimization problem into a series of easily solvable subproblems by introducing auxiliary variables, and gradually approaches the optimal solution through alternating iterations. The high-resonance sparse coefficient vector, low-resonance sparse coefficient vector, first auxiliary variable vector, second auxiliary variable vector, first Lagrange multiplier vector, and second Lagrange multiplier vector are initialized to all zeros. A penalty parameter of 0.5 is set; this parameter controls the strength of the auxiliary variable constraint, with a larger value indicating a stronger consistency constraint between the auxiliary and original variables. The convergence threshold is set to 0.000001, and the maximum number of iterations is set to 500.

[0062] Perform iterative optimization until convergence. Each iteration of the iterative optimization performs the following operations in sequence:

[0063] After subtracting the low-resonance component reconstructed based on the current low-resonance sparse coefficient vector from the original air quality index sequence, a forward wavelet decomposition is performed to update the high-resonance sparse coefficient vector. Specifically, the current low-resonance sparse coefficient vector is first updated by the inverse transform of a low-quality-factor tunable wavelet transform filter bank. Reconstructed into low resonance components Then, the difference between the original air quality index sequence and the low resonance component is calculated. The difference is then subjected to forward wavelet decomposition using a high-quality factor tunable wavelet transform filter bank to obtain the target coefficient vector. The updated high-resonance sparse coefficient vector is then calculated by combining the first auxiliary variable vector and the first Lagrange multiplier vector.

[0064] After subtracting the high-resonance component reconstructed based on the current high-resonance sparse coefficient vector from the original air quality index sequence, forward wavelet decomposition is performed to update the low-resonance sparse coefficient vector. The execution process is similar to that of updating the high-resonance sparse coefficient vector. First, the high-resonance component is reconstructed, and after calculating the residual, it is decomposed through a low-quality-factor tunable wavelet transform filter bank. The updated low-resonance sparse coefficient vector is then calculated by combining the second auxiliary variable vector and the second Lagrange multiplier vector.

[0065] Soft threshold shrinkage is performed on the first auxiliary variable vector and the second auxiliary variable vector respectively. Soft threshold shrinkage is a method to achieve... The key operation of norm sparsity constraints is to shrink coefficients with small absolute values ​​to zero, thereby promoting the sparsity of the coefficient vector. The soft-threshold shrinkage update is performed as follows: when the absolute value of an element is greater than the quotient of the regularization parameter and the penalty parameter, the element value is subtracted from the product of the quotient of the regularization parameter and the penalty parameter and the element's sign to obtain the shrunken element value; when the absolute value of an element is less than or equal to the quotient of the regularization parameter and the penalty parameter, the shrunken element value is set to zero.

[0066] Let the coefficient elements to be processed be The regularization parameter is The penalty parameter is Then the threshold for soft threshold shrinkage is .when When the element value after shrinking is ,in express The sign function, when When it is a positive number, it takes the value of 1; when When it is negative, the value is -1; when At that time, the value of the shrunk element is 0. For the first auxiliary variable vector, the threshold is the quotient of the high-resonance regularization parameter and the penalty parameter, i.e. For the second auxiliary variable vector, the threshold is the quotient of the low-resonance regularization parameter and the penalty parameter, i.e. The larger threshold of the second auxiliary variable vector means that more coefficients will be shrunk to zero, reflecting the design intention to impose stronger sparsity constraints on low-resonance components.

[0067] refer to Figure 4 The horizontal axis in the graph represents the input coefficient values. The value ranges from -3 to +3. The vertical axis represents the output coefficient value after soft thresholding. The value ranges from -2.5 to +2.5. The black dashed line in the figure is the reference line. This represents a linear relationship where input and output are equal without shrinkage. The blue solid line in the figure represents the threshold parameter. The corresponding soft-threshold shrinkage function curve, where the threshold is determined by the high-resonance regularization parameter. With penalty parameters The quotient is determined. The red solid line in the figure represents the threshold parameter. The corresponding soft-threshold contraction function curve, where the threshold is determined by the low-resonance regularization parameter. With penalty parameters The quotient is determined. The mathematical expression for the soft threshold shrinkage function is: when hour, ;when hour, Three characteristic regions of the soft-threshold contraction function can be observed from the curve shape. A dead zone is formed within the interval, shrinking all input coefficient values ​​within this interval to zero. This forces small-amplitude coefficients to zero, a key mechanism for promoting the sparsity of the coefficient vector. In this region, the output value equals the input value minus the threshold, and the curve shifts downwards relative to the reference line. Units. In In this region, the output value equals the input value plus a threshold, and the curve shifts upward relative to the reference line. Units. By comparing the two curves, it can be found that the threshold The corresponding dead zone width is the threshold. Five times that of the previous method, this means that more coefficients in the low-resonance sparse coefficient vector will be shrunk to zero, reflecting the design intent of imposing stronger sparsity constraints on the low-resonance component. During the iterative optimization process, the soft threshold shrinkage operation is applied alternately to the auxiliary variables corresponding to the high-resonance sparse coefficient vector and the low-resonance sparse coefficient vector, gradually causing the two coefficient vectors to tend towards a sparse distribution, ultimately achieving effective separation of the trend component and the transient component.

[0068] Update the first and second Lagrange multiplier vectors. The update of the Lagrange multiplier vectors is adjusted based on the difference between the current sparse coefficient vector and the auxiliary variable vector. The update of the first Lagrange multiplier vector is the current value plus the product of the penalty parameter and the difference between the high-resonance sparse coefficient vector and the first auxiliary variable vector. The update of the second Lagrange multiplier vector is the current value plus the product of the penalty parameter and the difference between the low-resonance sparse coefficient vector and the second auxiliary variable vector.

[0069] After each iteration, a convergence index is calculated. The convergence index is defined as the sum of the Euclidean norm of the difference between the current high-resonance sparse coefficient vector and the high-resonance sparse coefficient vector of the previous iteration, and the Euclidean norm of the difference between the current low-resonance sparse coefficient vector and the low-resonance sparse coefficient vector of the previous iteration. The Euclidean norm is defined as the square root of the sum of the squares of the elements in the vector. The iteration process terminates when the convergence index is less than the convergence threshold or the maximum number of iterations is reached. In typical operation of this embodiment, the algorithm usually converges after 150 to 200 iterations.

[0070] High-resonance components are obtained by performing inverse wavelet reconstruction on the converged high-resonance sparse coefficient vector. The inverse wavelet reconstruction process is as follows: the coefficients of each sub-band are upsampled according to the corresponding upsampling factor. Upsampling is achieved by inserting zero values ​​between adjacent coefficients. The upsampled coefficients are then filtered and accumulated sequentially through each level of comprehensive filters, reconstructing from the deepest level to the shallowest level, and finally outputting a reconstructed sequence of the same length as the original signal. The reconstructed high-resonance components contain the trend components with multi-oscillation periodic characteristics of the original air quality index sequence.

[0071] The high-resonance component is used as the first-level trend component. This first-level trend component reflects the gradual change pattern of the air quality index over a longer time scale; its waveform is relatively smooth and includes the main characteristics of periodic patterns such as daily and weekly variations. The first-level residual sequence is obtained by subtracting the first-level trend component point by point from the original air quality index sequence. Let the original air quality index sequence be... The first-level trend component is Then the first-order residual sequence is ,in The first-level residual sequence contains all components of the original sequence except for the trend component, and will serve as the input for the subsequent second-level decomposition.

[0072] After resonant sparse signal decomposition, the long-period trend components have been removed from the first-level residual sequence, but it still contains various components with different time-scale characteristics. From the perspective of the physical mechanism of air quality change, the first-level residual sequence mainly contains two types of components: one is quasi-periodic fluctuation components related to factors such as daily cycles and traffic peak cycles, which exhibit relatively regular oscillation patterns in the time series; the other is irregular fluctuation components caused by factors such as sudden pollution events and random meteorological disturbances, whose changes lack obvious regularity. In order to achieve accurate prediction of the air quality index, it is necessary to further separate these two types of components and adopt appropriate prediction methods for their respective characteristics.

[0073] Pole-symmetric mode decomposition (PSM) is an adaptive signal decomposition technique based on local features. It decomposes complex signals layer by layer into several intrinsic mode components with different time scales. Compared to traditional empirical mode decomposition (EMD) methods, PSM effectively suppresses boundary effects during the decomposition process by employing a pole-symmetric extension strategy at signal boundaries. Boundary effects refer to the deviation in envelope fitting near signal endpoints due to insufficient data support, which affects the accuracy of the entire decomposition result. Pole-symmetric extension utilizes the distribution pattern of the signal's own extrema to construct virtual extension points at the boundaries, enabling the envelope to achieve a reasonable fitting effect even in the boundary region.

[0074] refer to Figure 3In the graph, the horizontal axis represents time in hours, ranging from 0 to 10. The vertical axis represents signal amplitude, ranging from -35 to +35. The blue solid line in the graph represents the first-level residual sequence to be decomposed. This sequence is obtained by subtracting the first-level trend component from the original air quality index sequence, containing a mixture of periodic fluctuations and random disturbances. The waveform shows that the sequence exhibits complex multi-frequency oscillation characteristics, with multiple local maxima and minima. Red triangles indicate detected maxima, and green inverted triangles indicate detected minima. A maximum is a position where the signal value is simultaneously greater than both the preceding and following sample values, and a minimum is a position where the signal value is simultaneously less than both the preceding and following sample values. Within a 10-hour time span, the sequence contains approximately 20 maxima and 20 minima; the distribution density of these extreme points reflects the dominant oscillation frequency of the signal. The red dashed line in the diagram represents the upper envelope, generated by using cubic spline interpolation on the extended set of maxima. The upper envelope passes through all maxima, forming a smooth transition curve between them, and overall constitutes the envelope of the signal's local maxima. The green dashed line represents the lower envelope, generated by using the same cubic spline interpolation on the extended set of minima. The lower envelope passes through all minima, forming the envelope of the signal's local minima. The black solid line represents the envelope mean, composed of the arithmetic mean of the upper and lower envelopes at each sampling time, calculated as follows: ,in Indicates the upper envelope. The lower envelope is represented by the envelope mean. The envelope mean represents the local average level of the signal at various times, reflecting the relatively low-frequency trend components in the signal. Subtracting the envelope mean from the original sequence yields the candidate eigenmode components. When a candidate component satisfies the condition that the absolute value of the difference between the number of extrema and the number of zero crossings is less than or equal to 1, it is output as the eigenmode component.

[0075] Pole-symmetric mode decomposition (PSM) is performed on the first-level residual sequence, separating multiple intrinsic mode components sequentially through pole-symmetric extension and envelope mean extraction. The first step of PSM is to detect extreme points in the first-level residual sequence. Extreme points are the locations where the local trend of signal change reverses, including both maxima and minima. Maxima are the turning points where the signal changes from an upward trend to a downward trend, and minima are the turning points where the signal changes from a downward trend to an upward trend. The distribution characteristics of extreme points directly reflect the local oscillation characteristics of the signal and form the basis for subsequent envelope construction.

[0076] The specific execution method for extreme point detection is as follows: traverse every sampling point in the first-level residual sequence except for the first and last points, and for the first point in the sequence... The sampling point is compared with the adjacent sampling point. The value of the sampling point and the first sampling point The values ​​of the nth sampling point are compared. When the nth sampling point value is... The value of the sampling point is simultaneously greater than the value of the first sampling point. The value of the sampling point and the first sampling point When the value of the sampling point is obtained, the value of the first sampling point will be... The sampling point is marked as a local maximum. This condition, where the value is simultaneously greater than both its immediate and adjacent points, ensures that the point is indeed at a local maximum. When the... The value of the sampling point is simultaneously less than the value of the first sampling point. The value of the sampling point and the first sampling point When the value of the sampling point is obtained, the value of the first sampling point will be... Each sampling point is marked as a minimum point.

[0077] All maxima are arranged in chronological order to form a set of maxima, and all minima are arranged in chronological order to form a set of minima. Let the length of the first-order residual sequence be... The set of maxima contains The set of local maxima and minima contains [a set of local maxima and minima]. There are several minimum points. In this embodiment, for a first-order residual sequence of length 720, the number of maximum points is typically about 80 to 120, and the number of minimum points is similar to the number of maximum points. The set of maximum points can be represented as... ,in Indicates the first The time index of a maximum point Indicates the first The values ​​of the maximum points. The set of minimum points is represented in a similar way.

[0078] After extreme point detection, extreme point symmetric extension is performed on the sets of maxima and minima. The purpose of extreme point symmetric extension is to construct virtual extended extreme points at the sequence boundaries, so that subsequent cubic spline interpolation can obtain stable fitting results in the boundary region. Without boundary extension, cubic spline interpolation will produce large fitting errors near the endpoints due to the lack of constraints, manifested as abnormal upward or downward tilting of the envelope at the endpoints. This error will gradually accumulate and spread towards the middle of the sequence as iterative screening proceeds.

[0079] The extreme point symmetric continuation process is performed by mirroring the extreme points at the boundary, using the first and last sampling points of the first-level residual sequence as the center of symmetry. For the continuation process on the left boundary of the sequence, the first sampling point of the first-level residual sequence is used as the center of symmetry. Let the time index of the first sampling point be... The value is Select the first three maxima from the set of maxima, and let their time indices be respectively... , , The values ​​are respectively , , By mirroring these three maxima points relative to the first sampling point, three left-side extended maxima points are obtained, with their time indices as follows: , , The numerical values ​​remain unchanged. The geometric meaning of mirror flipping is to symmetrically transform the original extreme points with the endpoints of the sequence as the axis, and the flipped points are symmetrically distributed with respect to the origin with respect to the endpoints.

[0080] The same approach is used to handle the left boundary extension of the minimum point set, as well as the right boundary extension of both the maximum and minimum point sets. For the right boundary extension, the last sampling point of the first-order residual sequence is used as the center of symmetry. The last three extreme points from each of the maximum and minimum point sets are selected and mirrored. Selecting three extreme points for extension is a proven and reasonable choice; too few extension points will lead to insufficient interpolation boundary constraints, while too many extension points will introduce excessive virtual information that interferes with the decomposition results.

[0081] The extended maxima, all maxima in the extended maxima set, and the extended maxima on the right are arranged chronologically to form the extended maxima set. Similarly, the extended minima, all minima in the extended minima set, and the extended minima on the right are arranged chronologically to form the extended minima set. The extended set of extrema adds three virtual extrema at each end of the sequence, providing sufficient boundary constraints for subsequent interpolation calculations.

[0082] The upper envelope is generated using cubic spline interpolation based on the extended set of maxima. Cubic spline interpolation is a numerical method for constructing smooth curves between adjacent data points. It possesses the desirable property of continuous second derivatives, generating visually smooth and natural fitted curves. Compared to linear interpolation, cubic spline interpolation avoids sharp angles at data points; compared to higher-order polynomial interpolation, cubic spline interpolation avoids endpoint oscillations caused by the Runge phenomenon.

[0083] The specific implementation method of cubic spline interpolation is as follows: a cubic polynomial curve segment is constructed between every two adjacent maxima in the extended set of maxima. Let the two adjacent maxima be... and In the interval Construct a cubic polynomial on ,in , , , The coefficients to be determined are: Let be the time variable. The function value of each curve segment at the extended maximum point must be equal to the ordinate value of that point, i.e. and This ensures that the interpolation curve passes through all the maximum points.

[0084] To ensure the smoothness of the entire interpolation curve, the first and second derivatives of adjacent curve segments at the connection points must be continuously equal. Continuity of the first derivative ensures that the curve has no sharp corners at the connection points, while continuity of the second derivative ensures a smooth transition in the curvature of the curve at the connection points. These continuity conditions form a system of linear equations concerning the coefficients of the curve segments. Since the number of equations is two fewer than the number of unknowns, boundary conditions are needed to ensure a unique solution to the system of equations. Natural boundary conditions are used, setting the second derivative values ​​at both ends to zero. These natural boundary conditions assume that the curve has no curvature at the endpoints, which is a computationally simple and physically meaningful boundary treatment method.

[0085] After solving the above system of linear equations to obtain the coefficients of all cubic polynomial curve segments, the curve segments are connected to form the upper envelope. The upper envelope passes through all the maxima and presents a smooth transition curve between the maxima, forming an overall envelope of the signal's local maxima. Using the same cubic spline interpolation method, the lower envelope is generated based on the extended set of minima. The lower envelope passes through all the minima, forming an envelope of the signal's local minima.

[0086] Calculate the arithmetic mean of the upper and lower envelopes at each sampling time, and then construct the envelope mean line by arranging all the arithmetic means in chronological order. Let the upper envelope be at the th sampling time... The value at each sampling time is The lower envelope is at the first The value at each sampling time is Then the envelope mean line is at the th The value at each sampling time is ,in This represents the envelope mean. The envelope mean line represents the local average level of the signal at various times, reflecting the low-frequency trend components of the signal.

[0087] The intrinsic mode candidate components are obtained by subtracting the envelope mean line point by point from the current input sequence. In the first screening iteration, the current input sequence is the first-level residual sequence. Let the current input sequence be... The envelope mean line is Then the candidate components of the intrinsic modes are Subtracting the envelope mean from the signal eliminates the local trend component of the signal, making the resulting sequence symmetrically distributed relative to the zero mean.

[0088] When an intrinsic mode candidate component satisfies the intrinsic mode function (IMF) condition, the IMF candidate component is output as the intrinsic mode component. The IMF condition is the criterion for determining whether a signal component has good oscillation characteristics. The IMF condition is that the absolute value of the difference between the number of extrema and the number of zero-crossings is less than or equal to 1. A zero-crossing is the position where the signal value crosses zero when it changes from positive to negative or from negative to positive. For an ideal oscillating signal, there must be a zero-crossing between adjacent maxima and minima; therefore, the number of extrema and the number of zero-crossings should be equal or at most differ by 1. When the absolute value of the difference between the number of extrema and the number of zero-crossings of an IMF candidate component is greater than 1, it indicates that the component has not yet completely separated the highest frequency oscillation component, and further screening iterations are needed.

[0089] When a candidate intrinsic mode component (IMC) does not satisfy the IMC condition, the IMC is used as a new input sequence, and the entire process of extreme point detection, extreme point symmetric extension, upper and lower envelope generation, envelope mean line calculation, and IMC extraction is repeated until the obtained IMC satisfies the IMC condition. This iterative screening process gradually removes residual low-frequency components from the IMC candidates, making them increasingly closer to the ideal IMC shape. In this embodiment, the extraction of each IMC typically requires 3 to 8 screening iterations.

[0090] After extracting the first intrinsic mode component (IMC), the first-level residual sequence is subtracted point by point from the current IMC to obtain the intermediate residual sequence. The highest frequency oscillation component is removed from the intermediate residual sequence, while lower frequency oscillation components and residual trends are retained. The complete process of extreme point detection, extreme point symmetric extension, upper and lower envelope generation, envelope mean calculation, and IMC candidate component extraction is repeated on the intermediate residual sequence to extract the second IMC. The updated intermediate residual sequence is then further decomposed, extracting the third and fourth IMC components sequentially, until the number of extreme points in the intermediate residual sequence is less than 3. When the number of extreme points is less than 3, it is no longer possible to construct an effective envelope for further decomposition. The intermediate residual sequence at this point is the final residual sequence, representing the monotonic trend components or constant components in the signal that cannot be further decomposed.

[0091] In this embodiment, pole-symmetric mode decomposition is performed on the first-order residual sequence of length 720, typically extracting 4 to 6 intrinsic mode components. The intrinsic mode components are arranged in the extraction order, with the first intrinsic mode component containing the highest frequency oscillation, and the oscillation frequencies of subsequent intrinsic mode components decreasing sequentially. From the physical meaning of air quality changes, the first intrinsic mode component typically corresponds to rapid fluctuations on an hourly scale, the second intrinsic mode component typically corresponds to periodic changes on a daily scale, and the third and subsequent intrinsic mode components correspond to fluctuations or irregular disturbances on longer time scales.

[0092] The first two eigenmode components are summed point by point to form the second-order periodic component. Let the first eigenmode component be... The second intrinsic mode component is Then the second-order periodic component is The rationale for selecting the first two intrinsic mode components to form the second-level periodic components is that the first two intrinsic mode components mainly contain oscillatory components related to regular factors such as daily cycles and traffic peak cycles. These components have relatively stable periodic characteristics and are suitable for prediction using recurrent neural networks that can capture sequence dependencies.

[0093] The remaining intrinsic mode components are added point-by-point to the final residual sequence to form the third-level random components. Let the total number of extracted... The intrinsic mode components, and the final residual sequence are: Then the third-level random component is the third to the third. The result is the point-by-point summation of each intrinsic modal component and the final residual sequence. The third-level random component contains irregular fluctuation components and decomposition residuals. The changes in these components lack obvious periodicity, but may have local pattern features, making them suitable for prediction using convolutional neural networks that can extract local features.

[0094] In an optional implementation, the number of intrinsic mode components constituting the second-level periodic component can be adjusted according to the actual characteristics of the air quality index sequence. When the frequency range of the periodic component in the sequence is wide, the first three or four intrinsic mode components can be included in the second-level periodic component; when the periodic component is relatively concentrated in the high-frequency region, only the first intrinsic mode component can be used as the second-level periodic component. The determination is based on the frequency characteristics of each intrinsic mode component and its correspondence with known periodic factors.

[0095] The first-level trend component, the second-level periodic component, and the third-level stochastic component are predicted separately. The first-level trend component is predicted using a Long Short-Term Memory (LSTM) network. LSM is a recurrent neural network structure specifically designed for processing sequential data. It solves the gradient vanishing problem in long-sequence modeling of traditional recurrent neural networks by introducing a gating mechanism. The changes in the first-level trend component are relatively gradual but span a long period, requiring the prediction model to capture dependencies over a long time. The long-range memory capability of LSM makes it an ideal choice for predicting the trend component.

[0096] The forget gate calculates a forgetting coefficient vector between 0 and 1 based on the current input and the hidden state from the previous time step. A coefficient close to 0 indicates that the memory at the corresponding position should be forgotten, while a coefficient close to 1 indicates that the memory at the corresponding position should be retained. The output of the forget gate is... The calculation yielded, where Indicates the forget gate at time The output vector, This represents the sigmoid activation function. The weight matrix represents the forget gate. This represents the hidden state vector from the previous time step. This represents the input vector at the current moment. This represents the bias vector of the forget gate. This indicates a vector concatenation operation.

[0097] The input gate determines which parts of the current input information should be written into memory. The calculation of the input gate consists of two parts: first, calculating the input gate coefficient vector... ,in This represents the input gate coefficient vector. This represents the weight matrix of the input gate. This represents the bias vector of the input gate; then the candidate memory vector is calculated. ,in Represents the candidate memory vector. This represents the hyperbolic tangent activation function. The weight matrix representing the candidate memory. This represents the bias vector of the candidate memory.

[0098] The update of the memory unit is performed by combining the outputs of the forget gate and the input gate, and the calculation method is as follows: ,in Indicates time The memory cell state vector, This represents the state vector of the memory cell at the previous time step. This represents an element-wise multiplication operation. The element-wise multiplication of the forget gate coefficient with the previous time step memory achieves selective forgetting, while the element-wise multiplication of the input gate coefficient with the candidate memory achieves selective writing.

[0099] The output gate determines which information from the memory cell should be output to the hidden state. The calculation of the output gate is as follows: ,in This represents the output gate coefficient vector. This represents the weight matrix of the output gate. This represents the bias vector of the output gate. The hidden state is updated as follows: ,in Indicates time The hidden state vector. The hidden state serves as the network's output representation at the current time step and is used for subsequent prediction calculations.

[0100] In this embodiment, the Long Short-Term Memory (LSTM) network used to predict the first-level trend component employs a two-layer stacked structure, with each layer containing 64 memory units. The network input is the numerical sequence of the first-level trend component over the past 48 hours, and the output is the trend prediction sequence for the next 24 hours. After the input sequence is processed by the two LSM layers, the hidden state at the last time step is mapped to a 24-dimensional output vector, i.e., the trend prediction sequence, through a fully connected layer. The network training uses mean squared error as the loss function, employs the Adam optimizer for parameter updates, sets the learning rate to 0.001, and sets the training epochs to 100.

[0101] The second-level periodic component is predicted using a gated recurrent unit (ROU) network. The RNU is a simplified variant of the Long Short-Term Memory (LSTM) network, merging the forget gate and input gate into an update gate, while introducing a reset gate to control the utilization of historical information. The RNU has fewer parameters than the LSM, trains faster, and demonstrates good performance in predicting medium-length sequences. The second-level periodic component exhibits relatively regular periodic oscillations, with dependencies in the sequence mainly concentrated within a few cycles. The computational efficiency advantage of the RNU makes it a suitable choice for predicting periodic components.

[0102] The gated loop unit network contains two gated units: an update gate and a reset gate. The update gate determines how much information from the previous hidden state should be passed to the current state, calculated as follows: ,in Indicates the update gate at time The output vector, This represents the weight matrix of the updated gate. This represents the bias vector of the update gate. The reset gate determines how much information from the previous hidden state should be used to compute the candidate hidden state; the calculation method is as follows: ,in Indicates the time to reset the door. The output vector, This represents the weight matrix of the reset gate. This represents the bias vector of the reset gate.

[0103] The calculation of the candidate hidden state is as follows: ,in Represents the candidate hidden state vector. The weight matrix representing the candidate hidden state. This represents the bias vector of the candidate hidden state. The element-wise multiplication of the reset gate coefficients with the previous hidden state allows the network to selectively ignore some historical information. The hidden state is updated as follows: The update gate coefficient controls the mixing ratio of the hidden state and the candidate hidden state in the previous time step.

[0104] In this embodiment, the gated recurrent unit network used to predict the second-level periodic components adopts a two-layer stacked structure, with each layer containing 48 hidden units. The network input is the numerical sequence of the second-level periodic components over the past 48 hours, and the output is the predicted periodic sequence for the next 24 hours. The network training configuration is similar to that of the Long Short-Term Memory network, employing the mean squared error loss function and the Adam optimizer.

[0105] The third-level random component is predicted using a one-dimensional convolutional neural network (CNN). The CNN extracts features by sliding convolutional kernels across the input sequence, effectively capturing local pattern features. While the third-level random component lacks significant long-range dependencies, it may exhibit short-term local fluctuation patterns. The CNN's ability to extract local features makes it suitable for processing sequence data with local structural characteristics.

[0106] The one-dimensional convolution operation is performed by sliding the convolution kernel along the time dimension of the input sequence and calculating the inner product of the kernel and the corresponding input region at each position. Let the input sequence be... The kernel length is The convolution kernel parameters are Then the first convolution output sequence The elements are ,in The first convolution output represents the... One element, This represents the bias parameter. Convolution operations are typically followed by an activation function to introduce a nonlinear transformation; commonly used activation functions include the modified linear unit (MRU).

[0107] In this embodiment, the one-dimensional convolutional neural network used to predict the third-level random component comprises three convolutional layers and two fully connected layers. The first convolutional layer uses 32 convolutional kernels with a kernel length of 5; the second convolutional layer uses 64 convolutional kernels with a kernel length of 3; and the third convolutional layer uses 64 convolutional kernels with a kernel length of 3. Each convolutional layer is followed by a modified linear unit activation function and a max-pooling layer with a pooling window length of 2. The features extracted by the convolutional layers are flattened and then input into the fully connected layers. The first fully connected layer contains 128 neurons, and the second fully connected layer outputs a 24-dimensional random prediction sequence. The network input is the numerical sequence of the third-level random component over the past 48 hours.

[0108] After forecasting the three components, the prediction results of each component are summed point by point to obtain the air quality index prediction result. Let the trend prediction sequence be... The periodic prediction sequence is The random prediction sequence is The predicted air quality index is: The point-by-point addition reconstruction method corresponds to the three-level progressive residual decomposition process. During decomposition, the original sequence is subtracted from each component step by step, and during reconstruction, the prediction results of each component are added point by point to restore the complete prediction sequence.

[0109] Figure 5 This chart compares the prediction results of hierarchical forecasting and air quality index reconstruction, demonstrating the complete process of using three different neural networks to predict the three-level components separately and then superimposing the prediction results of each component to obtain the final air quality index prediction value. The horizontal axis represents time in hours, ranging from -48 to +24, where 0 is the prediction starting point, the negative range corresponds to historical observation data, and the positive range corresponds to the prediction period. The first row shows the prediction results of the first-level trend component, and the vertical axis represents the trend component value, ranging from 80 to 110. The blue solid line represents the historical trend component data of the past 48 hours, which is extracted from the historical air quality index sequence by resonant sparse signal decomposition and used as the input to the Long Short-Term Memory (LSTM) network. The green solid line represents the actual trend component value for the next 24 hours, used to evaluate prediction accuracy. The red dashed line represents the trend prediction value output by the LTM network. The LTM network uses a two-layer stacked structure, with each layer containing 64 memory units, and the input is the trend component sequence of the past 48 hours. The output is a trend prediction sequence for the next 24 hours. The red semi-transparent area represents the prediction confidence interval, reflecting the range of uncertainty in the prediction results. A comparison of the curves shows that the Long Short-Term Memory (LSTM) network can capture the gradual changes in the trend component well, and the predicted curve generally follows the actual curve. The second row in the figure shows the prediction results for the second-level periodic components. The vertical axis represents the periodic component values, ranging from -40 to +40, with a zero-value horizontal line separating the positive and negative regions. The blue solid line represents the historical periodic component data for the past 48 hours, the green solid line represents the actual periodic component values ​​for the next 24 hours, and the red dashed line represents the periodic prediction values ​​output by the gated recurrent unit network.

[0110] The Gated Recurrent Unit (GRU) network employs a two-layer stacked structure, with each layer containing 48 hidden units. Compared to Long Short-Term Memory (LSTM) networks, it has fewer parameters and faster training speed. The periodic components exhibit regular oscillatory characteristics. By learning periodic patterns in historical sequences, the GRU network can effectively predict the oscillation waveforms of future periods. The third row in the figure shows the prediction results for the third-level random components. The vertical axis represents the random component values, ranging from -20 to +20. The blue solid line represents the historical random component data of the past 48 hours, the green solid line represents the actual random component values ​​for the next 24 hours, and the red dashed line represents the random prediction values ​​output by the one-dimensional convolutional neural network. The one-dimensional convolutional neural network contains three convolutional layers and two fully connected layers, extracting local features by sliding the convolutional kernel across the input sequence. While the changes in random components lack long-range dependencies, they exhibit short-term local fluctuation patterns. The local feature extraction capability of the one-dimensional convolutional neural network makes it suitable for processing such sequences. The fourth row in the figure shows the reconstructed prediction results for the air quality index. The vertical axis represents the air quality index values, ranging from 40 to 160. The blue solid line represents the historical air quality index series over the past 48 hours, the green solid line represents the actual air quality index value for the next 24 hours, and the red dashed line represents the reconstructed forecast value obtained by overlaying the three components. The reconstructed forecast value is calculated as follows: ,in This is a trend forecast value. This is a periodic forecast value. These are random predicted values. The purple double-headed arrows in the figure indicate the deviations between the predicted and actual values ​​at certain times, visually demonstrating the magnitude of the prediction error. From the overall comparison, it can be observed that the reconstructed prediction curve can track the changing trend of the actual air quality index well, verifying the effectiveness of the three-level progressive residual decomposition combined with the hierarchical prediction strategy.

[0111] This invention is not limited to the specific embodiments described above. The invention extends to any new feature or combination disclosed in this specification, as well as any new method or process step or combination disclosed herein.

Claims

1. An air quality index prediction method based on three-level progressive residual decomposition, characterized in that, Includes the following steps: Step 1: Collect air quality monitoring data from the target monitoring station, preprocess the air quality monitoring data, calculate the air quality index value at each monitoring time, and arrange the air quality index values ​​of continuous monitoring times in chronological order to construct the original air quality index sequence. Step 2: Perform resonance sparse signal decomposition on the original air quality index sequence. By constructing a high-quality factor tunable wavelet transform filter bank and a low-quality factor tunable wavelet transform filter bank, perform dual-channel wavelet decomposition on the original air quality index sequence. Use the sparse coefficient optimization solution based on split augmented Lagrange contraction to obtain the high resonance component. Use the high resonance component as the first-level trend component. Subtract the first-level trend component from the original air quality index sequence point by point to obtain the first-level residual sequence. Step 3: Perform pole-symmetric mode decomposition on the first-level residual sequence. Separate multiple intrinsic mode components in turn through pole-symmetric extension and envelope mean extraction. Add the first two intrinsic mode components point by point to obtain the second-level periodic components. Add the remaining intrinsic mode components to the final residual sequence point by point to obtain the third-level random components. Step 4: Predict the first-level trend component, the second-level periodic component, and the third-level random component respectively, and add the prediction results of each component point by point to obtain the air quality index prediction result.

2. The method according to claim 1, characterized in that, In step one, the air quality monitoring data includes the concentrations of fine particulate matter, inhalable particulate matter, sulfur dioxide, nitrogen dioxide, carbon monoxide, and ozone; the preprocessing includes linear interpolation of missing values ​​and median filtering of outliers.

3. The method according to claim 1, characterized in that, In step two, the construction process of the high-quality factor tunable wavelet transform filter bank is as follows: The high-quality factor is set to 4, and the high redundancy parameter is set to 3. Based on the high-quality factor, the passband cutoff frequency of the high-pass filter is calculated to be at position 0.8 of the normalized frequency axis, and the passband cutoff frequency of the low-pass filter is calculated to be at position 0.2 of the normalized frequency axis. A high-pass finite impulse response (FIR) filter with a second-order zero at the zero position of the frequency axis and a corresponding low-pass FIR filter are constructed, ensuring that the high-pass FIR filter and the low-pass FIR filter meet the complete reconstruction condition. The high-pass FIR filter and the low-pass FIR filter are cascaded according to the decomposition level determined by the high redundancy parameter to form a high-quality factor tunable wavelet transform filter bank.

4. The method according to claim 3, characterized in that, The construction process of the low-quality factor tunable wavelet transform filter bank is as follows: set the low-quality factor value to 1, set the low redundancy parameter value to 3, calculate the passband cutoff frequency position of both the high-pass and low-pass filters based on the low-quality factor value, and the normalized frequency axis position is 0.

5. Construct the high-pass finite impulse response filter and the low-pass finite impulse response filter in the same way as the high-quality factor tunable wavelet transform filter bank, and cascade them to form the low-quality factor tunable wavelet transform filter bank.

5. The method according to claim 1, characterized in that, In step two, the execution process of dual-channel wavelet decomposition is as follows: the original air quality index sequence is passed through the low-pass and high-pass filters of the filter bank in sequence. After each filtering stage, it is downsampled according to the factor determined by the redundancy parameter. The decomposition is carried out step by step until the preset decomposition level is reached, and the low-pass sub-band coefficient sequences and high-pass sub-band coefficient sequences of each stage are output. All sub-band coefficient sequences obtained by decomposition based on the high-quality factor tunable wavelet transform filter bank are combined into a high-resonance wavelet coefficient vector, and all sub-band coefficient sequences obtained by decomposition based on the low-quality factor tunable wavelet transform filter bank are combined into a low-resonance wavelet coefficient vector.

6. The method according to claim 5, characterized in that, The optimization solution for sparse coefficients based on split augmented Lagrange contraction includes: initializing the high-resonance sparse coefficient vector, low-resonance sparse coefficient vector, first auxiliary variable vector, second auxiliary variable vector, first Lagrange multiplier vector, and second Lagrange multiplier vector as all zero vectors; setting the penalty parameter, high-resonance regularization parameter, low-resonance regularization parameter, convergence threshold, and maximum number of iterations; performing iterative optimization until convergence; and performing inverse wavelet reconstruction on the converged high-resonance sparse coefficient vector to obtain the high-resonance component.

7. The method according to claim 6, characterized in that, Each iteration of the iterative optimization performs the following operations in sequence: Subtract the low-resonance component reconstructed based on the current low-resonance sparse coefficient vector from the original air quality index sequence, then perform forward wavelet decomposition to update the high-resonance sparse coefficient vector; Subtract the high-resonance component reconstructed based on the current high-resonance sparse coefficient vector from the original air quality index sequence, then perform forward wavelet decomposition to update the low-resonance sparse coefficient vector; Perform soft-threshold shrinkage updates on the first and second auxiliary variable vectors respectively; Update the first and second Lagrange multiplier vectors; The soft-threshold shrinkage update is performed as follows: When the absolute value of an element is greater than the quotient of the regularization parameter and the penalty parameter, the element value is subtracted from the product of the quotient of the regularization parameter and the penalty parameter and the element's sign to obtain the shrunken element value; when the absolute value of an element is less than or equal to the quotient of the regularization parameter and the penalty parameter, the shrunken element value is set to zero.

8. The method according to claim 1, characterized in that, In step three, the pole-symmetric mode decomposition includes: performing extreme point detection on the first-level residual sequence, marking sampling points whose values ​​are simultaneously greater than the values ​​of the adjacent preceding and following sampling points as maxima, and marking sampling points whose values ​​are simultaneously less than the values ​​of the adjacent preceding and following sampling points as minima; forming a set of maxima in chronological order from all maxima and minima in chronological order from all minima; and performing pole-symmetric extension processing on the sets of maxima and minima, wherein the pole-symmetric extension processing is performed by mirroring the extreme points at the boundary with the first and last sampling points of the first-level residual sequence as the symmetry center.

9. The method according to claim 8, characterized in that, The pole-symmetric mode decomposition also includes: generating an upper envelope based on the extended set of maxima using cubic spline interpolation, and generating a lower envelope based on the extended set of minima using cubic spline interpolation; calculating the arithmetic mean of the upper and lower envelopes at each sampling time, and forming an envelope mean line by arranging all arithmetic means in chronological order; subtracting the envelope mean line from the current input sequence point by point to obtain candidate intrinsic mode components; and outputting the candidate intrinsic mode components as intrinsic mode components when they satisfy the intrinsic mode function conditions. The condition is that the absolute value of the difference between the number of extreme points and the number of zero-crossing points is less than or equal to 1; the first-level residual sequence is subtracted point by point from the current intrinsic mode component to obtain the intermediate residual sequence. The extreme point detection, extreme point symmetric extension processing, upper and lower envelope generation, envelope mean line calculation and intrinsic mode candidate component extraction are repeatedly performed on the intermediate residual sequence until the number of extreme points in the intermediate residual sequence is less than 3; the first two intrinsic mode components are added point by point as the second-level periodic components, and the remaining intrinsic mode components are added point by point to the final residual sequence as the third-level random components.

10. The method according to claim 1, characterized in that, In step four, a long short-term memory network is used to predict the first-level trend component to obtain a trend prediction sequence, a gated recurrent unit network is used to predict the second-level periodic component to obtain a periodic prediction sequence, and a one-dimensional convolutional neural network is used to predict the third-level random component to obtain a random prediction sequence.

Images