Regional heating distribution prediction method based on multi-source data fusion

Abstract

The present application relates to the technical field of heat distribution, and more particularly to a regional heat distribution prediction method based on multi-source data fusion, comprising: constructing an original observation matrix with time dimension data as rows and feature data of each heat exchange station as columns, performing difference calculation on each element in the original observation matrix along the time dimension to obtain an incremental matrix, and extracting the absolute value median of each row of the incremental matrix plus a smoothing term as a scale factor of each row. The present application reconstructs a feature sequence from a low-rank matrix and a sparse matrix, generates a dimensionless space weight based on the topological distance of the heat supply network for weighted fusion, and obtains a fusion feature through limiting, inverse hyperbolic tangent processing and scale factor reverse amplification. The fusion feature, the current heat distribution amount and the future outdoor environment temperature are input into a pre-trained neural network, which can make the model comprehensively refer to historical changes, spatial correlation and meteorological conditions, thereby improving the accuracy and stability of heat distribution prediction.

Description

Technical Field

[0001] This invention belongs to the field of heating distribution technology, and in particular relates to a regional heating distribution prediction method based on multi-source data fusion. Background Technology

[0002] In the daily operation and scheduling management of district heating systems, accurately predicting the heating demand and allocation of each heat exchange station is crucial for achieving on-demand heating, optimizing energy allocation, reducing operational energy consumption, and ensuring user comfort. Over long-term operation, district heating systems continuously accumulate a large amount of data on the status characteristics of heat exchange stations, such as supply and return water temperatures, flow rates, pressures, heating loads, and historical allocations. This data not only reflects the operational trends of each heat exchange station over time but is also influenced by the heating network topology, hydraulic transmission paths, and heat coupling relationships between stations, resulting in spatial correlations and mutual constraints between different heat exchange stations. Furthermore, heating loads and allocations are driven by external climatic factors such as future outdoor temperatures and weather changes. Therefore, district heating operation data typically exhibits characteristics of high dimensionality, multiple sources, high noise levels, and spatiotemporal coupling.

[0003] In related technologies, for example, Chinese patent document with publication number CN120317455B discloses a method and system for predicting heat load for district heating. It generates a dual-modal base adjacency matrix through a gating mechanism, generates a hidden feature matrix and a dynamic propagation adjacency matrix based on diffusion graph convolution, fuses the two adjacency matrices to generate a gated adaptive adjacency matrix, extracts spatial features based on diffusion graph convolution, and extracts temporal features based on spatiotemporal embedding and multi-head attention mechanisms, thereby reducing prediction bias and improving the speed of heat load prediction while ensuring accuracy.

[0004] However, in actual heating operation scenarios, factors such as aging of on-site sensor equipment, fluctuations in acquisition accuracy, communication network delays, signal interference, and sudden changes in local pipeline network conditions can easily contaminate the raw observation data with environmental noise, anomalous changes, and redundant information. This type of data quality degradation weakens the expression of the true operating patterns and affects the subsequent model's identification of trends in heating demand changes. Existing conventional heating forecasting technologies often struggle to simultaneously achieve anomalous disturbance removal, effective fluctuation retention, spatiotemporal structure extraction, and multi-source feature fusion when processing such high-dimensional data subject to multiple disturbances, resulting in insufficient feature quality for the input prediction model. Especially when outdoor temperatures change rapidly, local station loads fluctuate abnormally, or pipeline network operating conditions are adjusted, prediction results are prone to deviation, leading to excessive or insufficient heating distribution, resulting in energy waste, reduced scheduling efficiency, or decreased user comfort. Therefore, how to reduce the impact of noise and outliers from raw observation data, extract key features that reflect temporal variation patterns and spatial topological relationships, and combine these with future outdoor environmental factors to predict the heating distribution of each heat exchange station is a technical problem that urgently needs to be solved in the field of regional heating scheduling and management. Summary of the Invention

[0005] To address the problems of high noise interference, difficulty in removing abnormal fluctuations, insufficient extraction of spatiotemporal correlation features, and low accuracy of heating distribution prediction in existing technologies, this invention proposes a regional heating distribution prediction method based on multi-source data fusion.

[0006] This invention proposes a regional heating distribution prediction method based on multi-source data fusion, comprising: S1: Constructing an original observation matrix with time-dimensional data as rows and feature data of each heat exchange station as columns; performing difference calculation on each element in the original observation matrix along the time dimension to obtain an increment matrix; extracting the median absolute value of each row of the increment matrix and adding a smoothing term as the scaling factor for each row; dividing each element in the increment matrix by the scaling factor of its row and processing it with the hyperbolic tangent function to obtain a compressed matrix; S2: Performing singular value decomposition on the compressed matrix to obtain a singular value sequence; taking the position where the ratio of adjacent singular values ​​in the singular value sequence first falls below a truncation threshold as the truncation rank; retaining the first truncation rank singular values; subtracting a set shrinkage threshold from the singular values ​​after the truncation rank to perform a shrinkage operation and setting negative values ​​to zero to reconstruct a low-rank matrix. S1: Subtract the low-rank matrix from the compressed matrix to obtain the residual matrix. Set the elements in the residual matrix whose absolute value column ranking is lower than the ranking threshold to zero to obtain the sparse matrix. S2: Add the corresponding elements of the low-rank matrix and the sparse matrix to obtain the reconstructed feature matrix. Extract each column of the reconstructed feature matrix and arrange them along the time dimension to obtain the feature sequence corresponding to each heat exchange station. For any heat exchange station, combine the spatial topological distance of the heating pipeline network between the heat exchange station and other heat exchange stations to perform weighted summation calculation on the feature sequence corresponding to each heat exchange station. Combine the scale factor to restore the dimension and amplitude to obtain the fusion feature of the heat exchange station. Input the fusion feature of each heat exchange station, the actual heat distribution of each heat exchange station at the current time and the future outdoor ambient temperature into the pre-trained neural network model to output the predicted value of the heat distribution of each heat exchange station.

[0007] By employing a complete methodological chain—time difference scaling factor standardization, hyperbolic tangent amplitude limiting, singular value decomposition low-rank reconstruction, residual sparsification, spatial topology weighted fusion, and scale reduction before inputting into a pre-trained neural network—systematic purification and structural reorganization are completed before the data enters the prediction model. This effectively distinguishes the overall heating operation trend from significant local fluctuations, eliminates sensor white noise and abrupt outliers, and explicitly models the hydraulic topological correlations between stations, thereby improving the quality of input features at the source.

[0008] Further, the compression matrix is ​​obtained by: subtracting the corresponding element of the previous row from the element of the current row in the original observation matrix to obtain the increment matrix; calculating the absolute value of each element in each row of the increment matrix, taking the median of the absolute values ​​in the row and adding a smoothing term as the scaling factor of the row; dividing the value of each element in the increment matrix by the scaling factor of the row to obtain the standardized ratio; and inputting the standardized ratio into the hyperbolic tangent function for mapping to obtain the compressed matrix after amplitude limiting.

[0009] This invention fully utilizes the statistical characteristic that the median is insensitive to outliers, so that large-scale abnormal increments caused by equipment failures or extreme weather are smoothly truncated to around ±1, while the load increments within the normal fluctuation range maintain approximately linear scaling. Compared with conventional extreme value or mean standardization, it is less likely to be skewed by extreme values, effectively solving the problems of large differences in numerical dimensions and strong interference from abrupt outliers among different heat exchange stations.

[0010] Furthermore, singular value decomposition (SVD) is performed on the compressed matrix to obtain a singular value sequence. This includes: performing SVD on the compressed matrix to obtain a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix; and extracting the non-zero elements from the singular value diagonal matrix and arranging them in descending order to form the singular value sequence. This clearly arranges the heating backbone pattern and random noise, which were originally mixed in high-dimensional features, along the energy level, making the main features representing the common heating control mode of the system and the weak features representing random environmental disturbances numerically distinguishable. This provides a regular and orderly data foundation for subsequently adaptively determining the amount of information to retain based on energy distribution.

[0011] Furthermore, the position where the ratio of adjacent singular values ​​in the singular value sequence first falls below a truncation threshold is taken as the truncation rank. This involves: calculating the ratio of two adjacent singular values ​​in the singular value sequence to form a ratio sequence; comparing each element in the ratio sequence with the set truncation threshold along the sequence direction, and determining the index position corresponding to the first occurrence of a ratio below the truncation threshold as the truncation rank. This method achieves adaptive determination of the truncation rank, dynamically capturing the relative attenuation inflection point between the backbone information and noise based on the actual distribution of the singular values ​​in the data, which is extremely large in the first few terms and then drops sharply. Compared with the scheme of artificially pre-setting a fixed retained dimension, it avoids the problem of residual noise due to truncation that is too shallow or loss of effective structure due to truncation that is too deep.

[0012] Furthermore, the low-rank matrix is ​​reconstructed by: dividing the singular value sequence into two parts by a truncated rank, keeping the first truncated rank singular values ​​unchanged; for singular values ​​after the truncated rank, subtracting a set shrinkage threshold, and setting the difference to zero if it is negative, to obtain an updated singular value sequence; filling the updated singular value sequence into a singular value matrix or diagonal matrix that matches the singular value decomposition form, and performing matrix multiplication reconstruction by combining the transposes of the original left and right singular vector matrices to obtain the low-rank matrix. By keeping the first truncated rank principal singular values ​​unchanged, applying a soft threshold to shrink the tail singular values ​​after the truncated rank, and setting negative values ​​to zero before reconstructing the low-rank matrix, this method ensures that the key features of the system under trunk temperature fluctuations and hydraulic transmission delays are not lost, and also progressively suppresses the tail components corresponding to random noise and weak disturbances.

[0013] Furthermore, the sparse matrix is ​​obtained by setting zero to elements in the residual matrix whose absolute value column ranking is lower than the ranking threshold. This process includes: extracting elements from the residual matrix column by column, calculating the absolute value of each element, and sorting them in ascending order of absolute value within their respective columns, using the ranking of each element as a ranking marker; calculating the ratio of the truncated rank to the total rank of the compressed matrix, using the product of this ratio and the column length as the ranking threshold; and checking the ranking markers of each element in the residual matrix, setting zero to elements whose ranking markers are lower than the ranking threshold, thus obtaining the sparse matrix. This method achieves adaptive sparsity linked to the degree of low-rank principal structure extraction, allowing random minor fluctuations to be eliminated while retaining real local anomalous loads with prominent amplitudes. Compared to using a fixed amplitude threshold to screen residuals, this method can automatically adjust the sparsity based on the proportion of principal structures, avoiding the accidental deletion of significant mutations or excessive residual noise.

[0014] Further, the fusion characteristics of the heat exchange station are obtained, including: calculating the spatial topological distance of the heating network between the heat exchange station and other heat exchange stations; converting the spatial topological distance of the heating network into a dimensionless topological distance; generating initial spatial weights based on the dimensionless topological distance; normalizing the initial spatial weights to obtain the dimensionless spatial weight coefficients between the heat exchange station and other heat exchange stations; using the dimensionless spatial weight coefficients to perform weighted summation calculations on the feature elements of the corresponding heat exchange station in each feature sequence; limiting the amplitude of the weighted summed feature sequence to the domain of the inverse hyperbolic tangent function to obtain the amplitude-limited fusion sequence; and after processing the amplitude-limited fusion sequence with the inverse hyperbolic tangent function, multiplying it row by row by the corresponding scaling factor to restore it to the fluctuation level of the corresponding feature to obtain the fusion characteristics corresponding to the heat exchange station. This method explicitly models the stronger thermal coupling and hydraulic correlation between nearby heat exchange stations and the target site, enabling the prediction to take into account both the temporal variation pattern and the spatial network interaction characteristics. It also accurately restores the relative fluctuations within the standardized space to the original fluctuation levels of each feature. Compared with the scheme that only performs time series modeling and ignores the spatial constraints between stations, it can more realistically depict the spatiotemporal coupling relationship of district heating.

[0015] Furthermore, the fusion features, the actual heat distribution of each heat exchange station at the current moment, and the future outdoor ambient temperature are input into the pre-trained neural network model. This includes: acquiring outdoor ambient temperature data for the future forecast period provided by the meteorological system; using the fusion features as historical state change variables, the actual heat distribution of each heat exchange station at the current moment as the current state variable, and the future outdoor ambient temperature as the external condition variable, and concatenating these data to form the model input features, which are then used to input the pre-trained neural network model.

[0016] Furthermore, the predicted heat distribution values ​​for each heat exchange station are output, including: inputting the model input features into a pre-trained neural network model with a time-series processing structure, and calculating the predicted heat distribution values ​​for each heat exchange station in the corresponding future time period through forward propagation within the model. By concatenating the fused features as historical state change variables, the actual heat distribution at the current moment as the current state variable, and the future outdoor ambient temperature as the external condition variable into the model input, this invention enables the prediction model to simultaneously obtain historical operating inertia, real-time scheduling benchmarks, and future meteorological driving information.

[0017] Furthermore, the pre-trained neural network model is a long short-term memory network model, whose structure includes an input layer, at least one LSTM hidden layer, a fully connected regression layer, and an output layer.

[0018] Beneficial effects: First, by constructing a feature observation matrix for heat exchange stations and performing time difference, scale factor extraction, and hyperbolic tangent limiting, the impact of abnormal fluctuations and extreme outliers on prediction results can be weakened. The truncated rank is determined based on the ratio of adjacent singular values, and a low-rank matrix is ​​reconstructed by combining singular value shrinkage. A sparse matrix is ​​obtained through residual position thresholding, allowing for relative differentiation between overall operating trends and significant local fluctuations, thus improving the quality of input features. Second, feature sequences are generated by reconstructing the low-rank matrix and sparse matrix. Dimensionless spatial weights are generated based on the topological distance of the heating network for weighted fusion. After limiting, inverse hyperbolic tangent processing, and scale factor inverse amplification, fused features are obtained. Inputting the fused features, current heat allocation, and future outdoor ambient temperature into a pre-trained neural network allows the model to comprehensively consider historical changes, spatial correlations, and meteorological conditions, thereby improving the accuracy and stability of heat allocation prediction. Attached Figure Description

[0019] Figure 1 This is a flowchart of the heat distribution prediction method in an embodiment of the present invention;

[0020] Figure 2 This is a comparison diagram of the hyperbolic tangent mapping effect in the embodiments of the present invention; Figure 3 This is a time-series comparison diagram of the matrix decomposition results after scale reduction in an embodiment of the present invention. Detailed Implementation

[0021] This invention proposes a regional heating distribution prediction method based on multi-source data fusion, such as... Figure 1 As shown, it includes: S1, construct the observation matrix, and obtain the compression matrix through difference and amplitude limiting.

[0022] Specifically, an original observation matrix is ​​constructed with time-dimension data as rows and characteristic data of each heat exchange station as columns. The increment matrix is ​​obtained by performing difference calculation on each element in the original observation matrix along the time dimension. The median absolute value of each row of the increment matrix is ​​extracted and a smoothing term is added as the scaling factor of each row. The elements in the increment matrix are divided by the scaling factor of their respective rows and processed by the hyperbolic tangent function to obtain the compressed matrix.

[0023] The historical operating feature datasets of each heat exchange station, including supply and return water temperature, flow rate, and pressure, are read. The features of different dimensions are standardized in a dimensionless manner, and a two-dimensional array is initialized. The consecutive timestamps are sorted in ascending order as the row index of the two-dimensional array, and the feature data corresponding to each heat exchange station is used as the column data of the two-dimensional array to fill the original observation matrix.

[0024] A first-order difference operation is performed along the vertical axis of the original observation matrix, which is the time dimension. This involves subtracting the row vector elements of the previous time step from the row vector elements of the next time step to obtain the increment matrix. For each row of data in the increment matrix, the absolute values ​​of all elements in that row are calculated, and the median of all absolute values ​​is obtained. A smoothing term to prevent division by zero is added to this median. This sum is used as the scaling factor for each row. Each element in the increment matrix is ​​then divided by the scaling factor of its row, and the hyperbolic tangent of this ratio is calculated. Outliers with maxima are smoothly mapped to a numerical range of -1 to 1, resulting in a compression matrix.

[0025] In an optional embodiment, an original observation matrix is ​​constructed with time-dimension data as rows and characteristic data of each heat exchange station as columns. The increment matrix is ​​obtained by differencing each element in the original observation matrix along the time dimension. The median absolute value of each row of the increment matrix is ​​extracted and a smoothing term is added as the scaling factor for each row. Each element in the increment matrix is ​​divided by the scaling factor of its row and processed using the hyperbolic tangent function to obtain a compressed matrix, including: The increment matrix is ​​obtained by subtracting the corresponding element of the previous row from the element of the current row in the original observation matrix. The absolute value of each element in the increment matrix is ​​calculated, and the median of the absolute values ​​in the row is taken and a smoothing term to prevent division by zero is added as the scaling factor of that row. The value of each element in the increment matrix is ​​divided by the scaling factor of the row to obtain the standardized ratio. The standardized ratio is then input into the hyperbolic tangent function for mapping to obtain the compressed matrix after amplitude limiting.

[0026] A raw observation matrix of size M×N is collected and constructed, where M represents the time step (e.g., collecting hourly data within one day, where M=24), and N represents the number of feature columns for each heat exchange station. When each heat exchange station selects one feature, N corresponds to the number of heat exchange stations, for example, N=50. When each heat exchange station contains multiple features, N corresponds to the number of feature columns resulting from the combination of the number of heat exchange stations and the number of feature types. The elements in the matrix correspond to the dimensionless standardized feature data. The raw observation matrix is ​​traversed along the row direction, and the elements of each heat exchange station corresponding to the t-1th row are subtracted from the t-th row to obtain an increment matrix of size (M-1)×N, consisting of the increment row vectors corresponding to each difference step, where 2≤t≤M. This increment matrix mainly represents the short-term fluctuations of the data. For each difference step in the increment matrix, extract N elements from the row vectors. Take their absolute values ​​and calculate the median using a sorting algorithm or a fast selection algorithm. Add a smoothing term to this median to prevent division by zero anomalies. This median serves as the scaling factor for that difference step, resulting in a scaling factor column vector of length M-1. The preferred range of values ​​for the smoothing term is... to In the example, the smoothing term can be set to .

[0027] After obtaining the scale factor corresponding to each row, a scalar division operation is performed on each element in the increment matrix. The increment value in row t and column j is divided by the scale factor corresponding to row t, resulting in the standardized ratio after scale normalization, where 1 ≤ j ≤ N. This reduces the negative impact of numerical scale differences and abrupt outliers among different heat exchange stations. The standardized ratios in the entire increment matrix are input into the hyperbolic tangent function for activation mapping. This function nonlinearly limits and continuously maps the range of input standardized ratios to a value range greater than -1 and less than 1. During this process, the load increment within the normal fluctuation range can maintain approximately linear scaling, while large-scale anomalous increments caused by equipment failure or extreme weather are truncated to around ±1, outputting a well-structured and smooth extreme value compressed matrix, such as... Figure 2 As shown.

[0028] S2, singular value decomposition reconstructs the low-rank matrix, and residual selection yields a sparse matrix.

[0029] Specifically, singular value decomposition is performed on the compressed matrix to obtain a sequence of singular values. The position where the ratio of adjacent singular values ​​in the sequence first falls below the truncation threshold is taken as the truncation rank. The first truncation rank of singular values ​​is retained. The singular values ​​after the truncation rank are subtracted from the set shrinkage threshold to perform a shrinkage operation, and negative values ​​are set to zero. A low-rank matrix is ​​reconstructed. The low-rank matrix is ​​subtracted from the compressed matrix to obtain the residual matrix. Elements in the residual matrix whose absolute value column order is lower than the order threshold are set to zero to obtain a sparse matrix.

[0030] The obtained compressed matrix is ​​subjected to singular value decomposition, which yields the left singular vector matrix, singular value matrix, right singular vector matrix, and a singular value sequence arranged in descending order. Starting from the second term of the singular value sequence, the ratio of the current term to the previous term is calculated iteratively to obtain a ratio sequence. The array index position where the value is first less than a set truncation threshold is found by traversing the ratio sequence. The integer value corresponding to the index position is extracted and determined as the truncation rank. If no ratio less than the truncation threshold is found, the total rank of the compressed matrix is ​​determined as the truncation rank. The first truncation rank of important singular values ​​in the singular value sequence is kept unchanged. For singular values ​​after the truncation rank that represent random noise or weak perturbation, a soft threshold shrinkage operation is performed by subtracting the set shrinkage threshold. All negative values ​​generated after the shrinkage operation are compared with zero and uniformly replaced with zero.

[0031] The processed singular values ​​are filled into a diagonal matrix that matches the singular value decomposition form to obtain the reconstructed singular value matrix. The continuous matrix product between the left singular vector matrix, the reconstructed singular value matrix, and the transpose of the right singular vector matrix is ​​calculated in sequence and remapped back to the original dimension to generate a low-rank matrix. The low-rank matrix is ​​subtracted from the initial compressed matrix to obtain the residual matrix representing the anomalous mutation components. For each column vector of the residual matrix, it is first transformed into a non-negative sequence and then sorted in ascending order to extract the position label of each element. The position threshold is obtained by multiplying the ratio of the truncated rank to the total rank of the compressed matrix by the column length. Then, the elements with position labels lower than the position threshold are assigned a value of zero, and only the significant residual elements with higher positions are retained. The sparse matrix is ​​then output.

[0032] In an optional embodiment, singular value decomposition is performed on the compressed matrix to obtain a sequence of singular values. The position where the ratio of adjacent singular values ​​in the sequence first falls below a truncation threshold is taken as the truncation rank, including: Singular value decomposition is performed on the compressed matrix to obtain a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix. Non-zero elements in the singular value diagonal matrix are extracted and arranged in descending order to form a singular value sequence. The ratio of the (i+1)th singular value to the ith singular value in the singular value sequence is calculated to form a ratio sequence. Elements in the ratio sequence are compared with a set truncation threshold one by one along the direction of the ratio sequence. The index position corresponding to the first occurrence of a ratio less than the truncation threshold is determined as the truncation rank.

[0033] When performing decomposition and noise reduction on the (M-1)×N dimensional compressed matrix obtained after hyperbolic tangent mapping, the compressed matrix is ​​calculated using a full-rank or truncated singular value decomposition algorithm. The output is an orthogonal left singular vector matrix of size (M-1)×(M-1), a diagonal matrix containing r non-zero singular values, where r is the true rank of the matrix and r≤min(M-1,N), and an orthogonal right singular vector matrix of size N×N.

[0034] Extract all non-zero elements from the main diagonal of the diagonal matrix, ensuring they are arranged in descending order. This forms a sequence of singular values ​​containing important characteristic information. This sequence typically exhibits a few extremely large values ​​at the beginning, such as the first singular value being 85.3 and the second being 32.1, representing the core information of the system's heating pattern. The values ​​then drop sharply and tend to level off, such as the tenth singular value being 1.2 and the eleventh being 1.1, representing the characteristics of environmental random noise. Iterate through this descendingly ordered sequence of singular values ​​and calculate the relative attenuation ratio between the (i+1)th smaller singular value and the i-th adjacent larger singular value. This yields a ratio sequence containing r-1 elements. A preset truncation threshold is used, preferably between 0.05 and 0.25. For a typical urban heat exchange network, the truncation threshold is set to 0.15 in this example. Element comparison scanning is performed forward along the ratio sequence. If the ratio sequence first satisfies the condition when scanning the k-th element... If the relative decrease of two adjacent singular values ​​exceeds 85%, it indicates that the information content is low starting from the (k+1)th feature component, and it turns into Gaussian white noise and unstructured perturbation. Therefore, the index value k that triggers the judgment condition is specified as the optimal truncated rank on which the matrix low-rank reconstruction process is based.

[0035] In an optional embodiment, the first rank of singular values ​​is retained, and the singular values ​​after the rank truncation are subtracted from a set shrinkage threshold to perform a shrinkage operation, and negative values ​​are set to zero to reconstruct a low-rank matrix, including: The singular value sequence is divided into two parts by the truncated rank. The first truncated rank singular values ​​are kept unchanged. For the singular values ​​after the truncated rank, a set shrinkage threshold is subtracted. If the difference is negative, it is set to zero to obtain the updated singular value sequence. The updated singular value sequence is filled into a singular value matrix or diagonal matrix that matches the singular value decomposition form. The matrix is ​​then reconstructed by matrix multiplication using the transposes of the original left and right singular vector matrices to obtain the low-rank matrix.

[0036] Based on the determined truncation rank k as the dividing boundary, the original singular value sequence containing r elements is split into the main feature segment representing the common control mode of the heat exchange station community, i.e., index 1 to k, and the tail feature segment representing minor fluctuations and noise interference, i.e., index k+1 to r.

[0037] In practice, all larger singular values ​​from index 1 to k are kept unchanged to ensure that the system does not lose features when faced with trunk temperature fluctuations and hydraulic transmission delays. A soft-threshold shrinkage calculation is performed on the remaining smaller singular values ​​within the index range k+1 to r. Each element is subtracted from a set shrinkage threshold, which is calculated based on the variance of the data. The preferred value range is 0.5 to 2.0 times the average of the smaller singular values ​​before processing. For example, if the mean of the residual tail sequence is 0.8, the shrinkage threshold can be set to 1.2. After subtracting the shrinkage threshold, the difference is compared with 0 using a maximum comparison function. All negative differences caused by the subtraction operation are set to zero, retaining only the residual weak signals greater than the threshold. This forms an updated singular value sequence with sparse characteristics. After the shrinkage mapping is completed, the updated sparse singular value sequence is sequentially backfilled into the main diagonal positions corresponding to the initial all-zero matrix, constructing a purified and truncated diagonal matrix of new singular values.

[0038] Using standard matrix multiplication, the left singular vector matrix calculated and saved in the previous step and the right singular vector matrix after transpose are called. By performing continuous matrix multiplication reconstruction operations from left to right, this process removes redundant noise components and combines and superimposes the main directions of each order of features purified in the subspace to restore a low-rank matrix with a size of (M-1)×N. This low-rank matrix represents the regular heating distribution baseline drift state of all heat exchange station networks in the entire area as the outdoor ambient temperature changes, after removing high-frequency local extrema and sensor white noise.

[0039] In an optional embodiment, the residual matrix is ​​obtained by subtracting the low-rank matrix from the compressed matrix, and the elements in the residual matrix whose absolute value column order is lower than a certain threshold are set to zero to obtain a sparse matrix, including: Subtracting the low-rank matrix from the compressed matrix yields the corresponding residual matrix. Elements are extracted from the residual matrix column by column, and the absolute value of each element is calculated. The elements are then sorted in ascending order by absolute value within their respective columns. The position of each element is used as a position marker. The ratio of the truncated rank to the total rank of the compressed matrix is ​​calculated, and the product of this ratio and the column length is used as the position threshold. The position markers of each element in the residual matrix are checked one by one, and elements with position markers lower than the position threshold are set to zero, resulting in a sparse matrix.

[0040] By subtracting the reconstructed low-rank matrix of the same size from the compressed matrix of size (M-1)×N, a residual matrix containing only high-frequency interference, sudden sensor reading changes, and low-probability abnormal loads is obtained. N column vectors are extracted from the residual matrix, each column corresponding to M-1 time nodes. All elements within each column are then processed to their absolute values. A bubble sort or quicksort algorithm is applied within the column domain to sort the elements in ascending order of absolute value. The element with the smallest absolute value is assigned a rank of 1, and so on. The ascending rank index of all elements is assigned and used as the position marker for subsequent masking operations, so that random minor fluctuations correspond to lower position markers and sudden interferences correspond to higher position markers. Establish the ranking threshold calculation rule: Obtain the value of k / r according to the formula, where k is the determined truncation rank and r is the true rank of the matrix. This ratio represents the degree of extraction of the low-rank principal structure. Multiply this ratio by the total number of elements contained in each column, i.e., the time dimension length M-1, and perform a rounding operation to determine the ranking threshold for the column dimension truncation. For example, when k=35, r=40, and M-1=50, k / r=0.875, and the result after multiplication is 43.75, which can be rounded down to 43 as the ranking threshold. When the truncation rank is larger, the ranking threshold increases accordingly, and the number of small residual elements set to zero increases, thus making the residual matrix sparser.

[0041] Scan all elements in each column and perform a Boolean operation based on the element's corresponding rank marker and rank threshold. Elements with rank markers lower than the rank threshold are removed by rounding down to zero. Only significant residual elements with rank markers not lower than the rank threshold are retained in their original positions in the residual matrix. This step is used to filter out locally variable components with relatively prominent amplitudes from the residual matrix, resulting in a sparse matrix with a sparse distribution.

[0042] S3, the superimposed matrix fuses spatial features, and the model outputs the predicted value.

[0043] Specifically, the corresponding elements of the low-rank matrix and the sparse matrix are added together to obtain the reconstructed feature matrix. The columns of the reconstructed feature matrix are extracted and arranged along the time dimension to form a feature sequence. Combined with the spatial topological distance and scale factor of the heating network of the heat exchange station, the fused features are calculated. The fused features, the actual heat distribution of each heat exchange station at the current time and the future outdoor ambient temperature are input into the pre-trained neural network model to output the predicted value of the heat distribution of each heat exchange station.

[0044] The elements with the same row and column indices in the low-rank matrix and sparse matrix obtained in the previous steps are added together to generate a reconstructed feature matrix. This reconstructed feature matrix is ​​then traversed sequentially along the column index dimension. Each column, representing the historical time-series feature data of a specific heat exchange station, is extracted as a one-dimensional tensor and arranged in chronological order to form a continuous feature sequence. Based on the geographic information system node coordinates of the regional heating network system design drawings, the shortest path topological distance between the current predicted target heat exchange station and other spatially associated heat exchange stations along the actual pipeline route is calculated. This shortest path topological distance is converted into a dimensionless topological distance. The dimensionless topological distance is then transformed into initial spatial weights using the reciprocal of the distance with a smoothing term, an exponential decay function, or a Gaussian kernel function. Finally, the initial spatial weights are normalized to obtain dimensionless spatial weight coefficients.

[0045] Using a series of calculated spatial weight coefficients, a linear weighted summation operation is performed on the feature sequences corresponding to all the aforementioned associated heat exchange stations to integrate spatial interaction information. The one-dimensional result sequence obtained by weighted summation is first subjected to amplitude limiting processing to make it fall within the domain of the inverse hyperbolic tangent function, resulting in an amplitude-limited fusion sequence. After performing inverse hyperbolic tangent processing on the amplitude-limited fusion sequence, the data amplitude is inversely amplified by multiplying each element of the time step by the scale factor of the corresponding time step row saved above, generating a fusion feature representing the temporal variation law and the spatial pipeline network association interaction characteristics.

[0046] A long short-term memory network model was constructed. The network weight parameter file, which was pre-trained on a large historical heating scheduling dataset using the Adam optimization algorithm, was imported. The fusion features obtained from the above analysis were combined with the future outdoor ambient temperature values ​​obtained from the meteorological bureau data interface and the actual heating allocation values ​​of each heat exchange station at the current moment, which were transmitted back in real time from the field data acquisition system. The three heterogeneous data were spliced ​​and integrated into a unified input tensor along the feature dimension direction. This input tensor was passed to the pre-trained LSTM network model to trigger the forward propagation mechanism to calculate the deep spatiotemporal hidden state. The dimensionality reduction mapping output of the Linear layer at the end of the neural network model yielded the predicted heating allocation values ​​of each heat exchange station at the next scheduling cycle.

[0047] In an optional embodiment, the low-rank matrix and the sparse matrix are added element-wise to obtain a reconstructed feature matrix. The columns of the reconstructed feature matrix are extracted and arranged along the time dimension to form a feature sequence. Combined with the spatial topological distance and scale factor of the heat exchange station's heating network, the fused features are calculated, including: The low-rank matrix and the sparse matrix are added element-wise to generate a reconstructed feature matrix. Each column of the reconstructed feature matrix is ​​extracted and arranged in chronological order to form a feature sequence, yielding the feature sequence for each heat exchange station. For any heat exchange station, the spatial topological distance of the heating network between this station and all other heat exchange stations is calculated. This spatial topological distance is converted into a dimensionless topological distance. Initial spatial weights are generated based on the dimensionless topological distance and normalized to obtain the dimensionless spatial weight coefficients between this heat exchange station and all other heat exchange stations. The dimensionless spatial weight coefficients are used to perform a weighted summation of the feature elements of the corresponding heat exchange station in each feature sequence. The weighted summation feature sequence is then limited to the domain of the inverse hyperbolic tangent function to obtain a limited fusion sequence. After processing the limited fusion sequence with the inverse hyperbolic tangent function, it is multiplied by the corresponding scaling factor row-wise to restore the fluctuation level of the corresponding feature, thus obtaining the fusion feature corresponding to that heat exchange station.

[0048] The elements at corresponding positions in each row and column of the aforementioned low-rank matrix and sparse matrix of the same size are summed to reconstruct a feature matrix that takes into account both the main control mechanism of the system and the specific operating conditions of individual heat exchangers. One-dimensional column vectors of this reconstructed feature matrix are extracted station by station along the column direction. Each column contains M-1 time steps, which are aligned sequentially based on the order of occurrence to establish a feature sequence representing the longitudinal fluctuations of a single specific heat exchange station over time. For the target heat exchange station to be predicted, the geographic information and heating GIS topology system are retrieved to query the hydraulic extension topological distances between it and other surrounding heat exchange stations measured along the actual pipeline route. For example, the value ranges from 50m to 3000m.

[0049] Dividing the hydraulic extension topological distance by the maximum topological distance yields the dimensionless topological distance. The reciprocal of this dimensionless topological distance (after adding a preset smoothing term) is used to construct the initial spatial weights, which are then normalized to obtain dimensionless spatial weight coefficients. At the same time point, these dimensionless spatial weight coefficients are multiplied by the feature sequences of the corresponding neighboring stations, and a weighted summation aggregation operation is performed. This characterizes the network interaction features where closer heat exchange stations may have a stronger correlation with the target heat exchange station. After fusing the inter-station spatial hydraulic topological features, the time-series sequence array containing spatial fusion information is first limited to the domain of the inverse hyperbolic tangent function to obtain a limited fusion sequence. This limited fusion sequence is then passed into the inverse hyperbolic tangent function module, where an approximate inverse mapping restoration process expands and restores the new sequence data, confined to a compact limited range greater than -1 and less than 1, back to the original relative variation ratio range.

[0050] Based on the pre-stored scale factor vector, the elements of the scale-limited restored sequence are multiplied by the corresponding scale factor according to the time row index, completing the scale inverse restoration within the standardized feature space. This restores the relative fluctuation ratio to the fluctuation magnitude of the corresponding feature. The result, which includes hydraulic topological attributes and historical change features, is output as the fused feature. The time-series comparison of the matrix decomposition results after scale restoration is shown below. Figure 3 As shown.

[0051] In an optional embodiment, the fusion features, the actual heat distribution of each heat exchange station at the current moment, and the future outdoor ambient temperature are input into a pre-trained neural network model, which outputs a predicted value of the heat distribution of each heat exchange station, including: The system acquires outdoor ambient temperature data for the future forecast period provided by the meteorological system; it uses the fusion features as historical state change variables, the actual heat distribution of each heat exchange station at the current moment as the current state variable, and the future outdoor ambient temperature as the external condition variable, and concatenates these data to form the model input features; it then inputs these model input features into a pre-trained neural network model with a time series processing structure, and through forward propagation calculation within the model, outputs the predicted heat distribution values ​​for each heat exchange station in the corresponding future period.

[0052] When preparing the basic data and external boundary conditions for prediction driven by deep learning networks, the continuous hourly future outdoor dry-bulb temperature data of the period to be predicted is obtained through the standard meteorological service API interface. This data serves as the driving force variable for determining the heat dissipation rate of the building system. For example, it is predicted that the ambient temperature will fluctuate between -5°C and -12°C in the next 6 hours.

[0053] In the input data fusion preparation layer, a splicing channel is opened, and the fusion features obtained from the previous processing are selected as the time-series state change feature group representing the hysteresis and historical load slip characteristics of the heat exchange system. The heat meter values ​​at the current real scale of the SCADA system control port of each heat exchange station are collected simultaneously as the current state variable representing the base state of the pipeline network operation starting point. Combined with the external condition variable of future outdoor ambient temperature, the three types of feature groups are spliced ​​along the feature channel dimension and merged to construct a dimensionally regular model comprehensive input tensor matrix. For example, the tensor dimension size is the number of samples × time step × number of multimodal features.

[0054] The pre-trained neural network model is preferably a long short-term memory network model, whose structure includes an input layer, at least one LSTM hidden layer, a fully connected regression layer, and an output layer. The model input is a tensor matrix formed by concatenating historical state change variables, current state variables, and external condition variables. The historical state change variables include the aforementioned fused features, the current state variables include the actual heat distribution of each heat exchange station at the current time, and the external condition variables include the future outdoor ambient temperature. The model output is the predicted value of the heat distribution of each heat exchange station within the target prediction window.

[0055] During the model training phase, historical heating operation data, historical outdoor ambient temperature data, and the actual heating allocation for each heat exchange station during the corresponding time period are acquired. The fused features obtained after the aforementioned feature processing, the current heating allocation, and the future outdoor ambient temperature are used to construct training samples according to time windows. The actual heating allocation within the target prediction window is used as the supervision label. The training samples are input into the Long Short-Term Memory network model. The heating load change trend and state dependency within the historical time window are extracted recursively through the input gate, forget gate, and output gate of the LSTM hidden layer. The prediction result is then obtained by mapping through a fully connected regression layer. A loss function is constructed based on the error between the prediction result and the supervision label. The backpropagation algorithm and the Adam optimization algorithm are used to iteratively update the model parameters until the loss function converges or reaches the preset training rounds. The trained model weight parameters are then saved.

[0056] During the model inference phase, a tensor matrix containing fused features, the current actual heat allocation, and the future outdoor ambient temperature is input into the long short-term memory network model that has been trained and whose weight parameters are locked. The input tensor is recursively calculated by the LSTM hidden layer to extract the temporal dependency features, and then the fully connected regression layer performs dimensionality reduction mapping. Finally, the output layer outputs the predicted heat allocation values ​​covering each heat exchange station under the jurisdiction of the current heating system and corresponding to the next target prediction window.

[0057] This experiment uses a dataset of 50 heat exchange stations in a northern city operating continuously for 120 days as experimental conditions. The data sampling frequency is once per hour, and it is divided into training, validation, and test sets in a 7:2:1 ratio. The experimental hardware environment is deployed on a deep learning computing server equipped with a high-end graphics processor. The evaluation metrics are mean absolute error and root mean square error. The experiment sets up five control schemes: a basic control group using only a pre-trained neural network, ablation group 1 which removes hyperbolic tangent limiting and scale inverse restoration operations, ablation group 2 which removes singular value decomposition truncation and sparse matrix reconstruction steps, ablation group 3 which abandons the spatial topology weight fusion mechanism of heat exchange stations, and a complete scheme group deploying all algorithm steps.

[0058] The mean absolute error (MAE) of the basic control group was 26.45 kW, and the root mean square error (RMSE) was 33.18 kW. Ablation group one, lacking amplitude constraints and scale reduction for extreme value fluctuations, had an MAE of 21.32 kW and an RMS of 27.65 kW. Ablation group two, without singular value decomposition truncation and sparse residual reconstruction, had an MAE of 18.54 kW and an RMS of 24.21 kW. Ablation group three, without introducing a spatial topology weighting fusion mechanism for heat exchange stations, had an MAE of 15.67 kW and an RMS of 20.13 kW. The complete scheme group achieved a lower prediction error under experimental conditions, with an MAE of 11.05 kW and an RMS of 14.36 kW.

Claims

1. A method for predicting regional heating distribution based on multi-source data fusion, characterized in that, include: S1: Construct an original observation matrix with time dimension data as rows and characteristic data of each heat exchange station as columns. Perform difference calculation on each element in the original observation matrix along the time dimension to obtain the increment matrix. Extract the median absolute value of each row of the increment matrix and add a smoothing term as the scaling factor of each row. Divide the elements in the increment matrix by the scaling factor of their respective rows and process them with the hyperbolic tangent function to obtain the compressed matrix. S2: Perform singular value decomposition on the compressed matrix to obtain a singular value sequence. The position where the ratio of adjacent singular values ​​in the sequence first falls below a truncation threshold is taken as the truncation rank. The first truncation rank of singular values ​​is retained. The singular values ​​after the truncation rank are subtracted from the set shrinkage threshold to perform a shrinkage operation, and negative values ​​are set to zero. This reconstructs a low-rank matrix. The compressed matrix is ​​subtracted from the low-rank matrix to obtain the residual matrix. Elements in the residual matrix whose absolute value column ranking is lower than the ranking threshold are set to zero to obtain a sparse matrix. This includes: extracting elements from the residual matrix column by column, calculating the absolute value of each element, and sorting them in ascending order of absolute value within their respective columns, using the ranking of each element as a ranking marker; calculating the ratio of the truncation rank to the total rank of the compressed matrix, using the product of the ratio and the column length as the ranking threshold; and checking the ranking marker of each element in the residual matrix, setting elements whose ranking marker is lower than the ranking threshold to zero to obtain the sparse matrix. S3: Add corresponding elements of the low-rank matrix and the sparse matrix to obtain the reconstructed feature matrix. Extract each column of the reconstructed feature matrix and arrange them along the time dimension to obtain the feature sequence corresponding to each heat exchange station. For any heat exchange station, combine the spatial topological distance of the heating network between the heat exchange station and other heat exchange stations to perform weighted summation calculation on the feature sequence corresponding to each heat exchange station, and combine the scale factor to restore the dimension and amplitude to obtain the fusion feature of the heat exchange station. Input the fusion feature of each heat exchange station, the actual heat distribution of each heat exchange station at the current moment and the future outdoor ambient temperature into the pre-trained neural network model to output the predicted value of the heat distribution of each heat exchange station.

2. The regional heating distribution prediction method based on multi-source data fusion according to claim 1, characterized in that, The compression matrix is ​​obtained, including: Subtract the corresponding element of the previous row from the element of the current row in the original observation matrix to obtain the increment matrix; calculate the absolute value of each element in each row of the increment matrix, take the median of the absolute values ​​in the row and add a smoothing term as the scaling factor of the row; divide the value of each element in the increment matrix by the scaling factor of the row to obtain the standardization ratio; input the standardization ratio into the hyperbolic tangent function for mapping to obtain the compressed matrix after amplitude limiting.

3. The regional heating distribution prediction method based on multi-source data fusion according to claim 1, characterized in that, Singular value decomposition of the compressed matrix yields a sequence of singular values, including: Singular value decomposition is performed on the compressed matrix to obtain a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix. The non-zero elements in the singular value diagonal matrix are extracted and arranged in descending order to form a singular value sequence.

4. The regional heating distribution prediction method based on multi-source data fusion according to claim 1, characterized in that, The position where the ratio of adjacent singular values ​​in a singular value sequence first falls below the truncation threshold is taken as the truncation rank, including: Calculate the ratio of two adjacent singular values ​​in a sequence of singular values ​​to form a ratio sequence; The elements in the ratio sequence are compared one by one along the sequence direction with the set truncation threshold. The index position corresponding to the first ratio that is less than the truncation threshold is determined as the truncation rank. If no ratio is less than the truncation threshold, the total rank of the compressed matrix is ​​determined as the truncation rank.

5. The regional heating distribution prediction method based on multi-source data fusion according to claim 1, characterized in that, The reconstructed low-rank matrix includes: The singular value sequence is divided into two parts by the truncated rank, keeping the first truncated rank singular values ​​unchanged. For singular values ​​after the truncated rank, a set shrinkage threshold is subtracted, and if the difference is negative, it is set to zero to obtain the updated singular value sequence. The updated singular value sequence is filled into a singular value matrix or diagonal matrix that matches the singular value decomposition form. The matrix is ​​then reconstructed by matrix multiplication using the transposes of the original left and right singular vector matrices to obtain the low-rank matrix.

6. The regional heating distribution prediction method based on multi-source data fusion according to claim 1, characterized in that, The fusion characteristics of the heat exchange station were obtained, including: The spatial topological distance of the heating network between the heat exchange station and other heat exchange stations is calculated. This spatial topological distance is then converted into a dimensionless topological distance. Initial spatial weights are generated based on these dimensionless topological distances and normalized to obtain the dimensionless spatial weight coefficients between the heat exchange station and other heat exchange stations. The dimensionless spatial weight coefficients are then used to perform a weighted summation of the feature elements of the corresponding heat exchange station in each feature sequence. The weighted summation feature sequence is then limited to the domain of the inverse hyperbolic tangent function to obtain a limited fusion sequence. After processing the limited fusion sequence with the inverse hyperbolic tangent function, it is multiplied by the corresponding scaling factor row by row to restore the fluctuation level of the corresponding feature, thus obtaining the fusion feature corresponding to the heat exchange station.

7. The regional heating distribution prediction method based on multi-source data fusion according to claim 1, characterized in that, The neural network model is pre-trained by inputting the fusion features, the actual heat distribution of each heat exchange station at the current moment, and the future outdoor ambient temperature, including: Obtain outdoor ambient temperature data for the future forecast period provided by the meteorological system; use the fused features as historical state change variables, the actual heat distribution of each heat exchange station at the current moment as the current state variable, and the future outdoor ambient temperature as the external condition variable, and stitch the data together to form the model input features, which are then used to input the pre-trained neural network model.

8. The regional heating distribution prediction method based on multi-source data fusion according to claim 1, characterized in that, Output the predicted heat distribution values ​​for each heat exchange station, including: The model input features are fed into a pre-trained neural network model with a time series processing structure. Through forward propagation calculation within the model, the predicted heat distribution values ​​for each heat exchange station in the future time period are output.

9. The regional heating distribution prediction method based on multi-source data fusion according to claim 8, characterized in that, The pre-trained neural network model is a long short-term memory network model, whose structure includes an input layer, at least one LSTM hidden layer, a fully connected regression layer, and an output layer.

Images