A method and apparatus for natural gas load forecasting

Abstract

The present application relates to a kind of natural gas load prediction method and device, wherein, method includes: the time series data of natural gas load is collected, and the abnormal value in time series data is eliminated to obtain first time series data;Using improved integrated empirical mode decomposition method to decompose first time series data, obtain multiple intrinsic mode function and residual term;Each intrinsic mode component and residual term obtained by decomposition is combined with external variable respectively, and multiple enhanced data sequences are constructed;Information factor is calculated for each enhanced data sequence, and according to information factor, enhanced data sequence is divided into first subset and second subset;The data sequence in the first subset is input into Transformer encoder to extract the first feature containing global information, and the data sequence in the second subset is extracted by convolution processing second feature;After the first feature and the second feature are combined, input BiLSTM module, and output natural gas load prediction result.

Description

Technical Field

[0001] This invention relates to the technical field of natural gas load forecasting, and more particularly to a method and apparatus for natural gas load forecasting. Background Technology

[0002] Natural gas, as a relatively clean and efficient energy source, boasts numerous advantages such as high combustion efficiency and low pollution emissions, and is gradually becoming the dominant urban gas fuel. Accelerating the development of the natural gas industry not only helps optimize the energy structure, reduce dependence on traditional high-polluting energy sources, and lower carbon emissions, thus playing a vital role in environmental protection, but also meets the growing demand for high-quality energy from urban residents, promotes efficient industrial production, and fosters sustainable economic development. Therefore, accelerating the development of the natural gas industry has become a significant trend in the global energy sector.

[0003] In the construction of urban gas transmission and distribution systems, gas load forecasting plays a crucial role. Accurate gas load forecasting is essential for the scientific construction of urban gas transmission and distribution networks. By accurately predicting gas loads in different areas and at different times, it is possible to rationally plan the network layout, ensure the stability and efficiency of gas transmission, and avoid gas shortages or surpluses caused by unreasonable network design. For the design of gas storage and peak-shaving facilities, load forecasting provides key information. Based on the forecast results, appropriate gas storage capacity and peak-shaving capacity can be determined to effectively cope with peak and valley changes in gas consumption and ensure a stable supply of urban gas. In addition, the formulation of gas purchase plans by urban gas companies also relies heavily on gas load forecasting. Accurate forecasting helps companies rationally allocate gas purchases, optimize procurement costs, and improve the company's economic efficiency. Furthermore, in terms of intelligent operation and management of the pipeline network, load forecasting provides important data support for intelligent scheduling and equipment maintenance, contributing to improving the overall intelligence level and management efficiency of the pipeline network operation.

[0004] Currently, academic and industry research on natural gas load forecasting mainly employs three types of methods. The first type is statistical methods, primarily time series models. Time series models analyze and predict based on the time sequence of historical data, estimating future load by uncovering patterns in data changes over time. However, these methods often assume that the data is stationary and linear. In practical applications, natural gas load is affected by various complex factors such as seasons, climate, and economic activities, resulting in significant non-stationarity and nonlinear characteristics of the data. This limits the prediction accuracy of time series models when dealing with complex and changing load data.

[0005] The second category comprises artificial intelligence methods primarily based on fuzzy logic reasoning and neural network algorithms. Fuzzy logic reasoning can handle imprecise and fuzzy information, while neural network algorithms possess powerful nonlinear mapping capabilities, enabling them to automatically learn complex patterns in data. However, artificial intelligence methods typically require large amounts of high-quality data for training, and the quality and quantity of data directly impact model performance. In practice, obtaining sufficient and accurate natural gas load-related data is not easy. Furthermore, overfitting is prone to occur during model training, leading to poor generalization ability on new data and consequently affecting prediction accuracy.

[0006] The third category consists of other methods, primarily the grey model method. The grey model method is suitable for situations with limited data and incomplete or unclear information. It generates and processes raw data to uncover potential patterns between data points. However, the grey model requires a high degree of regularity in the data. When natural gas load is affected by sudden factors or abnormal situations, its prediction accuracy is significantly impacted, making it difficult to accurately reflect the true changes in load.

[0007] Therefore, although various methods exist for natural gas load forecasting, these studies all have certain limitations, and their forecasting accuracy needs further improvement. To better meet the needs of urban gas transmission and distribution system construction and operation management, there is an urgent need to develop a more accurate and reliable natural gas load forecasting method. Summary of the Invention

[0008] The technical problem to be solved by the present invention is to provide a method and apparatus for predicting natural gas load.

[0009] To achieve the above-mentioned objective, the present invention provides a natural gas load forecasting method, comprising the following steps: S1. Collect time-series data of natural gas load and remove load anomalies from the time-series data to obtain the first time-series data; S2. An improved integrated empirical mode decomposition method is used to decompose the first time series data to obtain multiple intrinsic mode functions and residual terms; S3. Combine the intrinsic modal components and residual terms obtained from the decomposition with the external variables to construct multiple enhanced data sequences; S4. Calculate the information factor for each augmented data sequence, and divide the augmented data sequence into a first subset and a second subset according to the relationship between the information factor and the preset threshold; S5. Input the data sequence in the first subset into the Transformer encoder to extract the first feature containing global information, and extract the second feature by convolution processing the data sequence in the second subset; S6. Combine the first feature and the second feature and input them into the BiLSTM module to output the natural gas load prediction result.

[0010] According to one aspect of the present invention, in step S1, the step of acquiring time-series data of natural gas load and removing load anomalies from the time-series data to obtain first time-series data includes using a neighborhood density-based method for anomaly detection and deleting the detected anomalies, which includes: Calculate the distance between target detection point P and its nth nearest neighbor. ; Using the selected nearest neighbor as the center, and the distance... A circular region is constructed with a radius of 1 / n. The number of neighboring points contained within the circular region is counted, and a weighting factor is calculated based on the number of neighboring points. ; Based on weighting factors Calculate the weighted average distance of the target detection point P, and calculate the weighted average distance of each nearest neighbor of the target detection point P; The neighborhood density of target detection point P is calculated based on the weighted average distance of target detection point P and the weighted average distance of its nearest neighbors. The anomaly factor of the target detection point P is calculated based on the domain density; When the abnormal factor is greater than a given first threshold, the target detection point P is marked as an abnormal point and deleted.

[0011] According to one aspect of the present invention, step S2, in which the first time-series data is decomposed using an improved integrated empirical mode decomposition method to obtain multiple intrinsic mode functions and residual terms, includes: S21. Randomly generated Gaussian white noise, uniformly distributed white noise, and Laplace noise are added to the first time series data respectively to obtain three new signals; S22. Perform EMD decomposition on the three new signals respectively to obtain three initial modal components; S23. Calculate the sequence complexity coefficient for each initial modal component; S24. The first intrinsic mode component is obtained by weighted averaging the generated initial mode components according to the sequence complexity coefficient. S25. Calculate the residual of the first time series data after removing the first intrinsic mode component. Repeat steps S31 to S35 with the residual as the new signal to be decomposed until the final residual is a monotonic function. Complete the decomposition of the first time series data and output all intrinsic mode components and the final residual as the residual term.

[0012] According to one aspect of the invention, in step S23, the sequence complexity coefficient for calculating each initial modal component is expressed as: ; in, The sequence complexity coefficient can also be represented as . , This represents the first intermediate parameter in the calculation process. This represents the second intermediate parameter in the calculation process. Indicates time The initial modal components can also be denoted as , Indicates time The initial modal components can also be denoted as , This represents the third intermediate parameter in the calculation process. This represents the fourth intermediate parameter in the calculation process. , Indicates two different times, and , .

[0013] According to one aspect of the present invention, step S4, which involves calculating an information factor for each enhanced data sequence and dividing the enhanced data sequence into a first subset and a second subset based on the relationship between the information factor and a preset threshold, includes: S41. Construct two operators for calculating the maximum difference of an enhanced data sequence at different time scales; wherein the two operators are the first operator and the second operator, respectively, and are expressed as: ; in, Indicates the first operator, Indicates the second operator, Represents an enhanced data sequence. , This indicates two different times; S42. For each time Statistical satisfaction time The quantity and recorded as and statistical satisfaction time The quantity and recorded as ,in, This indicates the first preset threshold. S43. Based on the obtained quantity and quantity The information factor for each enhanced data sequence is calculated; S44. Based on the information factor and a preset threshold, a size comparison is performed to split the enhanced data sequence into a first subset and a second subset.

[0014] According to one aspect of the invention, in step S4, in step S43, based on the obtained quantity and quantity In the step of calculating the information factor for each enhanced data sequence, the information factor is represented as: ; in, Information factors can also be represented as , This represents the fifth intermediate parameter in the calculation process. This represents the sixth intermediate parameter in the calculation process. This represents the maximum value of the time range of the enhanced data sequence.

[0015] According to one aspect of the present invention, in step S5, where the data sequence in the first subset is input into the Transformer encoder to extract a first feature containing global information, the Transformer encoder includes an embedding layer, a multi-head attention layer, and a feedforward layer; the step then includes: The data sequences in the first subset are transformed into high-dimensional embedding representations through an embedding layer, and positional encoding based on sine and cosine functions is added to generate features. ; Features are processed through a multi-head attention layer. A linear transformation is performed to obtain the query matrix Q, the key matrix K, and the value matrix V. A scaling factor is introduced to calculate the self-attention representation, and a multi-head attention mechanism is used to calculate the multi-head attention representation. Features Features are obtained by performing residual connections with multi-head attention representations and layer normalization. ; The features are processed by an activation function in the feedforward layer. Activation yields features ; Features and characteristics The residuals are connected and layer normalization is performed to obtain the first feature.

[0016] According to one aspect of the present invention, in step S5, the step of extracting the second feature from the data sequence in the second subset through convolution processing involves performing a sliding window operation on the data sequence of the second subset using a convolution kernel. The vectors at each position in the data sequence are multiplied by the weight vectors at the corresponding positions of the convolution kernel, and the results are summed. After adding a bias term, the second feature is obtained. The second feature is then expressed as: ; in, The second characteristic can be abbreviated as , This represents the weight vector of the convolution kernel at position i. Indicates the bias term. This represents each data sequence of the second subset.

[0017] According to one aspect of the present invention, step S6, which involves combining the first feature and the second feature and inputting the result into the BiLSTM module to output the natural gas load prediction result, includes: Combined features are obtained based on the first and second features; The LSTM structure is used to process the combined features of the input in the forward direction to calculate the forward hidden state; The input combined features are arranged in reverse order of time series, and the reverse data sequence is processed using the same LSTM structure to obtain the reverse hidden state; The forward and reverse hidden states are merged, imported into the output layer, and an activation function is used to obtain the natural gas load prediction result.

[0018] To achieve the above-mentioned objective, the present invention provides a natural gas load forecasting device applied to the aforementioned natural gas load forecasting method, comprising: The data acquisition module is used to acquire natural gas load data and organize it into time-series data; The preprocessing module is used to detect and remove outliers from the time series data to generate the first time series data. The decomposition module is used to decompose the first time series data into multiple intrinsic mode functions and residual terms using an improved integrated empirical mode decomposition method. The sequence construction module is used to combine the intrinsic mode components and residual terms obtained from decomposition with external variables to construct multiple augmented data sequences. The classification module is used to calculate the information factor for each augmented data sequence and divide the augmented data sequence into a first subset and a second subset according to the relationship between the information factor and a preset threshold. The feature extraction module is used to input the data sequence in the first subset into the Transformer encoder to extract the first feature containing global information, and to extract the second feature from the data sequence in the second subset through convolution processing. The prediction module combines the first feature and the second feature and inputs them into the BiLSTM module to output the natural gas load prediction result.

[0019] According to one aspect of the present invention, an outlier detection method based on neighborhood density can adaptively identify and remove local outliers in the load data, avoiding interference from outliers on subsequent decomposition and modeling. Simultaneously, an improved integrated empirical mode decomposition method is employed, fusing Gaussian white noise, uniformly distributed white noise, and Laplace noise, and combining this with a weighted average based on sequence complexity coefficients. This effectively overcomes the mode aliasing problem in traditional EMD methods, making the intrinsic mode functions obtained from the decomposition more physically meaningful and representative.

[0020] According to one aspect of the present invention, this scheme combines intrinsic modal components and residual terms with external variables (such as environmental parameters, weekday type, etc.) to construct an enhanced data sequence, fully integrating the multi-source influencing factors of natural gas load. Based on this, an information factor is introduced to adaptively partition the enhanced sequence. Sequences with high complexity are assigned to a Transformer encoder to extract global information, while sequences with low complexity are extracted for local features through convolution operations. Furthermore, this classification mechanism effectively avoids the incompatibility problem of a unified model with heterogeneous data, significantly improving the targeting and efficiency of feature extraction.

[0021] According to one embodiment of the present invention, this scheme significantly improves the accuracy of load forecasting results by integrating the advantages of multiple models. Specifically, the Transformer encoder excels at capturing long-term dependencies and global interaction information, convolution operations efficiently extract locally stationary patterns, and the BiLSTM module further models the forward and reverse dynamic evolution of the time series. The fusion of these three structures enables the model to handle complex non-stationary load sequences while preserving detailed variation features. The final output load forecasting results effectively improve the load forecasting capability of this scheme in terms of accuracy, stability, and generalization ability.

[0022] According to one aspect of the present invention, this scheme forms an end-to-end prediction framework from data cleaning, signal decomposition, feature enhancement, adaptive partitioning to multi-mode feature extraction and fusion. Corresponding processing mechanisms are designed to address common problems in natural gas load data such as noise, missing data, anomalies, and nonlinear fluctuations, enabling the model to maintain stable prediction performance under different external variables (such as environmental parameters, different seasons, different climatic conditions, and different gas consumption scenarios), thus possessing good engineering practical value.

[0023] According to one aspect of the present invention, the high-precision load forecast results output by this approach can be directly used for load scheduling, gas storage and peak shaving planning, gas purchase plan formulation, and intelligent operation management of urban gas pipeline networks. This helps to reduce the operating costs of gas companies, improve energy utilization efficiency, and support the digital and intelligent transformation of natural gas systems. Attached Figure Description

[0024] Figure 1 This is a flowchart illustrating the steps of the natural gas load forecasting method of the present invention. Figure 2 This is a flowchart of the natural gas load forecasting method of the present invention. Detailed Implementation

[0025] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. The embodiments cannot be described in detail here, but the embodiments of the present invention are not limited to the following embodiments.

[0026] Combination Figure 1 and Figure 2 As shown, according to one embodiment of the present invention, a natural gas load forecasting method includes the following steps: S1. Collect time-series data of natural gas load and remove load anomalies from the time-series data to obtain the first time-series data; S2. An improved integrated empirical mode decomposition method is used to decompose the first time series data to obtain multiple intrinsic mode functions and residual terms; S3. Combine the intrinsic modal components and residual terms obtained from the decomposition with the external variables to construct multiple enhanced data sequences; S4. Calculate the information factor for each augmented data sequence, and divide the augmented data sequence into a first subset and a second subset according to the relationship between the information factor and the preset threshold; S5. Input the data sequence in the first subset into the Transformer encoder to extract the first feature containing global information, and process the data sequence in the second subset through convolution to extract the second feature; S6. Combine the first feature and the second feature and input them into the BiLSTM module to output the natural gas load prediction result.

[0027] According to one embodiment of the present invention, in step S1, the step of collecting time-series data of natural gas load uses natural gas load data and organizes it into time-series data in chronological order.

[0028] According to one embodiment of the present invention, in step S1, the step of acquiring time-series data of natural gas load and removing load anomalies from the time-series data to obtain the first time-series data includes using a neighborhood density-based method for anomaly detection and deleting the detected anomalies, which includes: Calculate the distance between target detection point P and its nth nearest neighbor. In this embodiment, the target detection point P is said to have n nearest neighbor points. Then the target detection point P and its nth nearest neighbor point Distance between Represented as: ; in, This represents the target detection point P and its nearest neighbors. The Euclidean distance between them .

[0029] Using the selected nearest neighbor as the center, and the distance... A circular region is constructed with a radius of 1 / n. The number of neighboring points contained within the circular region is counted, and a weighting factor is calculated based on the number of neighboring points. In this embodiment, the neighboring points contained within the statistical circular region are... The number of elements that satisfy the set Neighboring points The quantity, and this quantity is denoted as The weighting factor is calculated based on the number of neighboring points. Represented as: .

[0030] Based on weighting factors Calculate the weighted average distance of the target detection point P, and calculate the weighted average distance of each nearest neighbor of the target detection point P; in this embodiment, the weighted average distance of the target detection point P is expressed as: ; in, This represents the weighted average distance to the target detection point P; Similarly, the same calculation method can be used to calculate each nearest neighbor of the target detection point P. The weighted average distance, denoted as ; The neighborhood density of target detection point P is calculated based on the weighted average distance of target detection point P and the weighted average distance of its nearest neighbors; in this embodiment, the neighborhood density of target detection point P is expressed as: ; in, This represents the neighborhood density of the target detection point P at order n. Used to calculate all That is to say ; The anomaly factor of the target detection point P is calculated based on the neighborhood density; in this embodiment, the anomaly factor is expressed as: ; in, Indicates an abnormal factor.

[0031] When the anomaly factor is greater than a given first threshold, the target detection point P is marked as an anomaly and deleted.

[0032] In this embodiment, by selecting different target detection points P, after removing all outliers, a new time series data, namely the first time series data, can be obtained. ,in, , This represents the maximum value of the time range for the first time series data.

[0033] According to one embodiment of the present invention, step S2, in which the first time-series data is decomposed using an improved integrated empirical mode decomposition method to obtain multiple intrinsic mode functions and residual terms, includes: S21. Randomly generated Gaussian white noise, uniformly distributed white noise, and Laplace noise are added to the first time-series data respectively to obtain three new signals; in this embodiment, Gaussian white noise, uniformly distributed white noise, and Laplace noise are added to the first time-series data respectively, thereby generating three new time-series data, i.e., three new signals, which are represented as follows: .

[0034] S22. Perform EMD decomposition on the three new signals respectively to obtain three initial modal components; among them, the new signal based on EMD decomposition (i.e., integrated empirical mode decomposition method) can be represented as: ; in, Represents the initial modal components. This represents the residual components after EMD decomposition.

[0035] S23. Calculate the sequence complexity coefficient for each initial modal component; in this embodiment, the sequence complexity coefficient is expressed as: ; in, The sequence complexity coefficient can also be represented as . , This represents the first intermediate parameter in the calculation process. This represents the second intermediate parameter in the calculation process. Indicates time The initial modal components can also be denoted as , Indicates time The initial modal components can also be denoted as , This represents the third intermediate parameter in the calculation process. This represents the fourth intermediate parameter in the calculation process. , Indicates two different times, and , .

[0036] Furthermore, for each initial modal component Each can obtain a corresponding sequence complexity coefficient. .

[0037] S24. The generated initial modal components are weighted and averaged according to the sequence complexity coefficient to obtain the first intrinsic modal component; in this embodiment, the first intrinsic modal component is obtained by weighting and averaging the three generated initial modal components according to the sequence complexity coefficient, and is expressed as: ; in, This represents the first intrinsic mode component.

[0038] S25. Calculate the residual of the first time-series data after removing the first intrinsic mode component. Repeat steps S31 to S35 with the residual as the new signal to be decomposed until the final residual is a monotonic function, thus completing the decomposition of the first time-series data. Output all intrinsic mode components and the final residual as the residual term. In this embodiment, the residual of the first time-series data after removing the first intrinsic mode component is expressed as: ; in, This represents the residual of the first time series data after removing the first intrinsic mode component; By repeating steps S31 to S35 with the obtained residual as the new signal to be decomposed, the second intrinsic mode component can be obtained. and the second residual ; Repeat the above steps until the final residual is a monotonic function and cannot be further decomposed. At this point, the decomposition process ends, completing the decomposition of the first time series data. Output all intrinsic mode components and the final residual as the residual term; where the original first time series data... The final decomposition form is expressed as: ; in, This represents the intrinsic mode components obtained from the decomposition. This represents the final residual obtained from the decomposition. This indicates the number of intrinsic modal components.

[0039] According to one embodiment of the present invention, in step S3, the process of combining each intrinsic mode component and residual term obtained from the decomposition with external variables to construct multiple enhanced data sequences involves combining the intrinsic mode components obtained from the decomposition with external variables. and residuals Combined with external variables to construct enhanced data sequences ,in , , This represents the maximum value in the time range of the enhanced data sequence. This indicates the number of sequences in the enhanced data sequence; in this embodiment, external variables include environmental data (such as temperature), day of the week, date, etc.

[0040] Through the above settings, this scheme uses the improved EMD method (integrating three types of noise) to decompose the IMFs (intrinsic mode functions) and residuals, which are then forcibly "bound" to external variables such as temperature and day of the week, constructing an enhanced data sequence that integrates multiple factors. This makes the input data for prediction more reasonable, thereby transforming the simple problem of natural gas load prediction into a multi-factor prediction scheme based on comprehensive consideration of multiple factors, thus making the prediction results more accurate.

[0041] According to one embodiment of the present invention, step S4, which involves calculating an information factor for each enhanced data sequence and dividing the enhanced data sequence into a first subset and a second subset based on the relationship between the information factor and a preset threshold, includes: S41. Construct two operators for calculating the maximum difference of an enhanced data sequence at different time scales; wherein the two operators are the first operator and the second operator, respectively, and are expressed as: ; in, Indicates the first operator, Indicates the second operator, Represents an enhanced data sequence. , This indicates two different times; S42. For each time Statistical satisfaction time The quantity and recorded as and statistical satisfaction time The quantity and recorded as ,in, This indicates the first preset threshold. S43. Based on the obtained quantity and quantity The information factor for each enhanced data sequence is calculated; in this embodiment, the information factor is represented as: ; in, Information factors can also be represented as , This represents the fifth intermediate parameter in the calculation process. This represents the sixth intermediate parameter in the calculation process. This represents the maximum value within the time range of the enhanced data sequence. Therefore, for each enhanced data sequence... This will allow you to obtain a corresponding information factor. .

[0042] S44. Based on the information factor and a preset threshold, a size comparison is performed to split the enhanced data sequence into a first subset and a second subset. In this embodiment, based on the information factor... Is it greater than the preset threshold? This allows for the enhancement of data sequences. Divided into two subsets, namely the first subset: Second subset The first subset All information factors are greater than or equal to the preset threshold. The second subset All information factors are less than the preset threshold. , The number of sequences in the first subset. The number of sequences in the second subset.

[0043] The above settings reflect the complexity of the data through the constructed information factor. The higher the information factor, the more complex the data. Therefore, dividing the data into two parts and selecting appropriate feature processing methods for each part will make the prediction results more accurate.

[0044] According to one embodiment of the present invention, in step S5, where the data sequence in the first subset is input into the Transformer encoder to extract the first feature containing global information, the Transformer encoder includes an embedding layer, a multi-head attention layer, and a feedforward layer. The step then includes: The data sequences in the first subset are transformed into high-dimensional embedding representations through an embedding layer, and positional encoding based on sine and cosine functions is added to generate features. Specifically, the location coding used is as follows: ; in, Indicates the number of dimensions. Represents an even number of dimensions. Represents odd-numbered dimensions. Indicates the position of each data point in the data sequence. The total number of dimensions of the embedded representation. For the first Position in each dimension is The location encoding of the data.

[0045] Features are processed through a multi-head attention layer. A linear transformation is performed to obtain the query matrix Q, key matrix K, and value matrix V. A scaling factor is introduced to calculate the self-attention representation, and a multi-head attention mechanism is used to calculate the multi-head attention representation. In this embodiment, the query matrix Q, key matrix K, and value matrix V are represented as follows: ; in, , , It is a parameter matrix.

[0046] Furthermore, by introducing a scaling factor To avoid gradient vanishing due to excessively large dot products and to maintain gradient stability during training, the self-attention representation is calculated: ; in, This represents the normalization function used to generate the probability distribution.

[0047] Furthermore, based on the self-attention mechanism, a multi-head attention mechanism is calculated to obtain the multi-head attention representation, specifically: ; in, For the first A self-attention matrix, , For splicing operations, , , This is the weight matrix. The parameter matrix is ​​used to compress the attention matrix. This indicates a multi-head attention mechanism. For the number of attention, This represents the self-attention mechanism.

[0048] Features Features are obtained by performing residual connections with multi-head attention representations and layer normalization. ; The features are processed by an activation function in the feedforward layer. Activation yields features ; Features and characteristics The residuals are connected and layer normalization is performed to obtain the first feature.

[0049] According to one embodiment of the present invention, in step S5, the step of extracting the second feature from the data sequence in the second subset through convolution processing involves using a convolution kernel to perform a sliding window operation on the data sequence of the second subset. The vector at each position in the data sequence is multiplied by the weight vector at the corresponding position of the convolution kernel, and the sum is obtained. After adding a bias term, the second feature is obtained. The second feature is then expressed as: ; in, The second characteristic can be abbreviated as , This represents the weight vector of the convolution kernel at position i. Indicates the bias term. This represents each data sequence of the second subset.

[0050] By dividing the enhanced data sequence into different subsets for separate processing, this scheme achieves parallel extraction of features from different subsets. This not only effectively improves data processing efficiency but also ensures the extraction accuracy of different data features. Through Transformer, long-range dependencies and global information are captured, while convolution efficiently extracts local features, reducing computational overhead. This parallel processing meets the high-precision requirements for different types of data, thus contributing to improved prediction accuracy. In particular, this processing method fully utilizes the characteristics of data with varying complexities to achieve complementary feature extraction, enhancing the ability to characterize multi-scale and multi-frequency components of natural gas load.

[0051] According to one embodiment of the present invention, step S6, which involves combining the first feature and the second feature and inputting the result into the BiLSTM module to output the natural gas load prediction result, includes: Combined features are obtained based on the first and second features; An LSTM structure is used to process the combined features of the input in the forward direction to calculate the forward hidden state. In this embodiment, the hidden state is calculated based on the input data features and the hidden state of the previous time step. The forgetting state is calculated using the forgetting gate, input gate, memory cell, and output gate respectively. Input status Unit status and output status The calculation formula is as follows: ; in, This represents the obtained positive hidden state. For activation function, Current moment The combined features are obtained based on the combination of the first and second features. Represents the weight matrix. Represents the weight parameter. Represents the bias term. It is the hyperbolic tangent activation function.

[0052] The input combined features are arranged in reverse order of time series, and the reverse data sequence is processed using the same LSTM structure to obtain the reverse hidden state; The forward hidden state and the reverse hidden state are merged, imported into the output layer, and an activation function is used to obtain the natural gas load prediction result. In this embodiment, the state after merging the forward hidden state and the reverse hidden state is represented as follows: ; in, Indicates a positive hidden state. This indicates a reversed hidden state.

[0053] By combining the forward and reverse processing of bidirectional LSTM with the above settings, this scheme captures the bidirectional dependencies of time series, enabling it to make more comprehensive use of historical information and improve the accuracy of prediction results.

[0054] According to one embodiment of the present invention, the present invention provides a natural gas load forecasting device applied to the aforementioned natural gas load forecasting method, comprising: The data acquisition module is used to acquire natural gas load data and organize it into time-series data; The preprocessing module is used to detect and remove outliers from the time series data to generate the first time series data. The decomposition module is used to decompose the first time series data into multiple intrinsic mode functions and residual terms using an improved integrated empirical mode decomposition method. The sequence construction module is used to combine the intrinsic mode components and residual terms obtained from decomposition with external variables to construct multiple augmented data sequences. The classification module is used to calculate the information factor for each augmented data sequence and divide the augmented data sequence into a first subset and a second subset according to the relationship between the information factor and a preset threshold. The feature extraction module is used to input the data sequence in the first subset into the Transformer encoder to extract the first feature containing global information, and to extract the second feature from the data sequence in the second subset through convolution processing. The prediction module combines the first feature and the second feature and inputs them into the BiLSTM module to output the natural gas load prediction result.

[0055] Specific limitations regarding the natural gas load forecasting device can be found in the limitations of the natural gas load forecasting method described above, and will not be repeated here. Each module in the aforementioned natural gas load forecasting device can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device in hardware form, or stored in the memory of a computer device in software form, so that the processor can call and execute the corresponding operations of each module.

[0056] In this embodiment, the memory may be, but is not limited to, Random Access Memory (RAM), Read Only Memory (ROM), Programmable Read-Only Memory (PROM), Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read-Only Memory (EEPROM), etc.

[0057] In this embodiment, the processor can be an integrated circuit chip with signal processing capabilities. The processor can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0058] The above description is merely an example of a specific solution of the present invention. For any devices and structures not described in detail herein, it should be understood that they are implemented using common devices and methods already available in the art.

[0059] The above description is merely one embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method of natural gas load forecasting, characterized by, Includes the following steps: S1. Collect time-series data of natural gas load and remove load anomalies from the time-series data to obtain the first time-series data; S2. An improved integrated empirical mode decomposition method is used to decompose the first time series data to obtain multiple intrinsic mode functions and residual terms; S3. Combine the intrinsic modal components and residual terms obtained from the decomposition with the external variables to construct multiple enhanced data sequences; S4. Calculate the information factor for each augmented data sequence, and divide the augmented data sequence into a first subset and a second subset according to the relationship between the information factor and the preset threshold; S5. Input the data sequence in the first subset into the Transformer encoder to extract the first feature containing global information, and extract the second feature by convolution processing the data sequence in the second subset; S6. Combine the first feature and the second feature and input them into the BiLSTM module to output the natural gas load prediction result.

2. The natural gas load forecasting method according to claim 1, characterized in that, In step S1, the step of collecting time-series data of natural gas load and removing load anomalies from the time-series data to obtain the first time-series data involves using a neighborhood density-based method for anomaly detection and deleting the detected anomalies. This includes: Calculate the distance between target detection point P and its nth nearest neighbor. ; Using the selected nearest neighbor as the center, and the distance... A circular region is constructed with a radius of 1 / n. The number of neighboring points contained within the circular region is counted, and a weighting factor is calculated based on the number of neighboring points. ; Based on weighting factors Calculate the weighted average distance of the target detection point P, and calculate the weighted average distance of each nearest neighbor of the target detection point P; The neighborhood density of target detection point P is calculated based on the weighted average distance of target detection point P and the weighted average distance of its nearest neighbors. The anomaly factor of the target detection point P is calculated based on the domain density; When the abnormal factor is greater than a given first threshold, the target detection point P is marked as an abnormal point and deleted.

3. The natural gas load forecasting method according to claim 2, characterized in that, Step S2, which involves using an improved integrated empirical mode decomposition method to decompose the first time-series data and obtain multiple intrinsic mode functions and residual terms, includes: S21. Randomly generated Gaussian white noise, uniformly distributed white noise, and Laplace noise are added to the first time series data respectively to obtain three new signals; S22. Perform EMD decomposition on the three new signals respectively to obtain three initial modal components; S23. Calculate the sequence complexity coefficient for each initial modal component; S24. The first intrinsic mode component is obtained by weighted averaging the generated initial mode components according to the sequence complexity coefficient. S25. Calculate the residual of the first time series data after removing the first intrinsic mode component. Repeat steps S31 to S35 with the residual as the new signal to be decomposed until the final residual is a monotonic function. Complete the decomposition of the first time series data and output all intrinsic mode components and the final residual as the residual term.

4. The natural gas load forecasting method according to claim 3, characterized in that, In step S23, the sequence complexity coefficient for calculating each initial modal component is expressed as follows: in, The sequence complexity coefficient can also be represented as . , This represents the first intermediate parameter in the calculation process. This represents the second intermediate parameter in the calculation process. Indicates time The initial modal components can also be denoted as , Indicates time The initial modal components can also be denoted as , This represents the third intermediate parameter in the calculation process. This represents the fourth intermediate parameter in the calculation process. , Indicates two different times, and , .

5. The natural gas load forecasting method according to claim 3 or 4, characterized in that, Step S4, which involves calculating the information factor for each enhanced data sequence and dividing the enhanced data sequence into a first subset and a second subset based on the relationship between the information factor and a preset threshold, includes: S41. Construct two operators for calculating the maximum difference of an enhanced data sequence at different time scales; wherein the two operators are the first operator and the second operator, respectively, and are expressed as: in, Indicates the first operator, Indicates the second operator, Represents an enhanced data sequence. , This indicates two different times; S42. For each time Statistical satisfaction Time The quantity and recorded as and statistical satisfaction Time The quantity and recorded as ,in, This indicates the first preset threshold. S43. Based on the obtained quantity and quantity The information factor for each enhanced data sequence was calculated. S44. Based on the information factor and a preset threshold, a size comparison is performed to split the enhanced data sequence into a first subset and a second subset.

6. The natural gas load forecasting method according to claim 5, characterized in that, In step S4, in step S43, based on the obtained quantity and quantity In the step of calculating the information factor for each enhanced data sequence, the information factor is represented as: in, Information factors can also be represented as , This represents the fifth intermediate parameter in the calculation process. This represents the sixth intermediate parameter in the calculation process. This represents the maximum value of the time range of the enhanced data sequence.

7. The natural gas load forecasting method according to claim 6, characterized in that, In step S5, where the data sequence in the first subset is input into the Transformer encoder to extract the first feature containing global information, the Transformer encoder includes an embedding layer, a multi-head attention layer, and a feedforward layer. The steps include: The data sequences in the first subset are transformed into high-dimensional embedding representations through an embedding layer, and positional encoding based on sine and cosine functions is added to generate features. ; Obtain the query matrix Q, key matrix K, and value matrix V through linear transformation of the features Calculate the self-attention representation by introducing a scaling factor, and perform multi-head attention mechanism calculation to obtain the multi-head attention representation; Features Features are obtained by performing residual connections with multi-head attention representations and layer normalization. ; The features are processed by an activation function in the feedforward layer. Activation yields features ; Features and characteristics The residuals are connected and layer normalization is performed to obtain the first feature.

8. The natural gas load forecasting method according to claim 7, characterized in that, In step S5, the step of extracting the second feature from the data sequence in the second subset through convolution processing involves using a sliding window operation on the data sequence of the second subset using a convolution kernel. The vector at each position in the data sequence is multiplied by the weight vector at the corresponding position of the convolution kernel, and the sum is calculated. After adding a bias term, the second feature is obtained. The second feature is then represented as: in, The second characteristic can be abbreviated as , This represents the weight vector of the convolution kernel at position i. Indicates the bias term. This represents each data sequence of the second subset.

9. The natural gas load forecasting method according to claim 8, characterized in that, Step S6, which involves combining the first feature and the second feature and inputting the result into the BiLSTM module to output the natural gas load prediction, includes: Combined features are obtained based on the first and second features; The LSTM structure is used to process the combined features of the input in the forward direction to calculate the forward hidden state; The input combined features are arranged in reverse order of time series, and the reverse data sequence is processed using the same LSTM structure to obtain the reverse hidden state; The forward and reverse hidden states are merged, imported into the output layer, and an activation function is used to obtain the natural gas load prediction result.

10. A natural gas load forecasting device applied to the natural gas load forecasting method according to any one of claims 1 to 9, characterized in that, include: The data acquisition module is used to acquire natural gas load data and organize it into time-series data; The preprocessing module is used to detect and remove outliers from the time series data to generate the first time series data. The decomposition module is used to decompose the first time series data into multiple intrinsic mode functions and residual terms using an improved integrated empirical mode decomposition method. The sequence construction module is used to combine the intrinsic mode components and residual terms obtained from decomposition with external variables to construct multiple augmented data sequences. The classification module is used to calculate the information factor for each augmented data sequence and divide the augmented data sequence into a first subset and a second subset according to the relationship between the information factor and a preset threshold. The feature extraction module is used to input the data sequence in the first subset into the Transformer encoder to extract the first feature containing global information, and to extract the second feature from the data sequence in the second subset through convolution processing. The prediction module combines the first feature and the second feature and inputs them into the BiLSTM module to output the natural gas load prediction result.

Images