A probability wind speed prediction method based on physical guidance and multi-scale feature fusion

Abstract

This invention provides a probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion, belonging to the interdisciplinary field of meteorological forecasting and computer science. The technical solution includes the following steps: S1: Preprocessing the dataset; S2: Constructing physically guided interactive features; S3: Extracting spatial features of the data at multiple scales using the MSFF-ResNet module; S4: Mining the time-series dynamic patterns of the data using the FAMD-Liquid Time-Constant Networks module; S5: Concatenating features using a gated fusion mechanism, calculating weights, and inputting the results into an MLP for prediction; S6: Designing a Bayesian-optimized composite loss function; S7: Using the trained model to perform multi-task wind speed prediction on a test set. This invention improves the accuracy and physical interpretability of probabilistic wind speed prediction under complex terrain.

Description

Technical Field

[0001] This invention relates to the field of meteorological forecasting and computer science, and in particular to a probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion. Background Technology

[0002] With the global energy structure transformation and the development of meteorological forecasting technology, wind energy's strategic position as a clean and renewable energy source is becoming increasingly prominent. Short-term wind speed forecasting is playing a crucial role in wind farm power dispatch, grid stability, and energy market transactions. Current numerical weather prediction (NWP) relies on global data and atmospheric dynamic equations for solving problems. While it is the mainstream method for forecasting meteorological elements, it struggles to accurately capture local wind field structures, surface roughness effects, and thermal driving mechanisms under complex terrain, leading to significant systematic biases in wind speed forecasts. Traditional statistical downscaling methods and machine learning models attempt to improve forecast accuracy, but most models only provide deterministic point predictions and cannot quantify the uncertainty of the forecast, making it difficult to meet the practical needs of wind power dispatch for probabilistic and interval forecasts.

[0003] In recent years, deep learning methods have been widely applied in meteorology due to their superior nonlinear fitting and feature fusion capabilities. In particular, the hybrid architecture of convolutional neural networks (CNNs) and recurrent neural networks (RNNs) has begun to show potential in wind speed prediction. However, existing models often employ single-scale convolutional kernels for spatial feature extraction, resulting in insufficient representation of the spatial heterogeneity of wind speed under complex terrain. Furthermore, time series modeling often relies on fixed gating mechanisms, leading to slow responses to the strong non-stationarity of wind speed, and lacks adaptive strategies for forgetting and retaining historical information. Current research on physical information fusion largely remains at the level of simple feature concatenation, resulting in insufficient physical interpretability of the models. Moreover, it often focuses on point prediction tasks, neglecting the need for uncertainty quantification in wind speed prediction, thus limiting its application value in risk assessment and decision support.

[0004] How to solve the above problems is the key to this invention. Summary of the Invention

[0005] The purpose of this invention is to provide a probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion, which aims to improve the accuracy and physical interpretability of probabilistic wind speed prediction under complex terrain.

[0006] The core idea of ​​this invention is as follows: This invention proposes a probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion. By integrating meteorological parameters from the THORPEX Interactive Global Encyclopedia (TIGGE) project, topographic data from the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 3 (ASTER GDEM V3), and temporal parameters, a physically guided interactive feature is constructed. This method deeply explores the interaction relationship between meteorological and topographic parameters. A dual-branch architecture is used to extract spatiotemporal dynamic features, and a multi-scale feature fusion module (MSFF) is combined to extract spatial information at different scales. A feature-aware adaptive memory decay mechanism (FAMD) dynamically adjusts the decay rate of the hidden states, achieving adaptive modeling of the time-series dependence of wind speed. Furthermore, an adaptive composite loss function is designed, and Bayesian optimization is used to automatically search for the optimal weight combination, achieving multi-objective collaborative optimization.

[0007] This invention is achieved through the following measures: a probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion, comprising the following steps:

[0008] 1.1: The meteorological parameters, topographic data of the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 3 (ASTER GDEM V3), and time parameters of the THORPEX Interactive Global Encyclopedia (TIGGE) project were preprocessed. The dataset was constructed through steps such as missing value handling, outlier removal using the 3σ rule, spatiotemporal linear interpolation, and standardization. The dataset was then divided into training, validation, and test sets according to the year.

[0009] 1.2: Based on atmospheric boundary layer dynamics theory and topographic radiation modulation mechanism, physical information-guided interactive features are constructed using meteorological and topographic parameters;

[0010] 1.3: Integrating a multi-scale feature fusion module (MSFF) into the ResNet architecture enables collaborative modeling of information at different spatial scales;

[0011] 1.4: Introducing a Feature-Aware Adaptive Memory Decay Module (FAMD) into the Liquid Time-Constant Networks architecture, and using an adaptive memory gating mechanism to capture the nonlinear temporal evolution of wind speed;

[0012] 1.5: The dual-branch features are concatenated using a gated fusion mechanism (GatedFusion), and the fusion weights are adaptively calculated before being input into the MLP to achieve prediction;

[0013] 1.6: Design a composite loss function, automatically search for loss weights through Bayesian optimization, and combine it with an adaptive moment estimation (AdamW) optimizer with weight decay, cosine annealing learning rate scheduling, and gradient pruning strategy to optimize model training;

[0014] 1.7: Use the trained model (PIMFNet) based on physical information-guided feature construction and Bayesian optimized composite loss function to make predictions on the test set, and output the predicted point values, probability distribution and prediction interval.

[0015] Step 1.1 includes the following steps:

[0016] 2.1: Missing values ​​were processed for meteorological parameters, topographic data of the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 3 (ASTER GDEM V3), and time parameters of the THORPEX Interactive Global Encyclopedia (TIGGE) project. Linear interpolation was used in the time dimension, and linear interpolation was first performed in the spatial dimension, followed by nearest neighbor interpolation.

[0017] 2.2: The 3σ rule was used to remove outliers that deviated from the mean by more than 3 times the standard deviation, and all data were uniformly interpolated to a 0.25°×0.25° grid to match the four time points of 06:00, 12:00, 18:00 and 24:00 each day;

[0018] 2.3: The meteorological and topographic parameters are standardized, and finally the training set, validation set and test set are divided according to the year.

[0019] In step 1.2, five physical interaction features are constructed based on the meteorological and terrain parameters in the dataset processed in step 1.1, including the following steps:

[0020] 3.1: Based on the theory of atmospheric boundary layer dynamics and the topographic radiation modulation mechanism, five physical information-guided interactive features are constructed using meteorological and topographic parameters;

[0021] 3.2: Constructing the wind speed-topographic relief interaction factor (WR_IF1), based on the drag effect of topographic roughness on wind speed, and using topographic relief to quantify the dynamic modulation effect of surface roughness on near-surface wind speed, the formula is as follows:

[0022]

[0023] in, This indicates the wind speed at a height of 10 meters. It represents the degree of terrain relief, reflecting the surface roughness characteristics. The standard deviation of the global terrain relief is used for normalization. Indicates the terrain drag coefficient;

[0024] 3.3: Based on the dynamic modulation effect of terrain slope on wind speed, a wind speed-slope coupling factor (WS_IF2) is constructed. The slope is used to quantify the influence of terrain inclination on wind acceleration or deceleration. The formula is as follows:

[0025]

[0026] in, Indicates the slope angle of the terrain. The tangent value represents the slope, reflecting the degree of terrain inclination. The effective projection coefficient of slope wind. This represents the slope modulation coefficient, whose value is based on the slope wind speed enhancement effect in topographic dynamics theory.

[0027] 3.4: Based on the theory of atmospheric temperature vertical lapse rate, an altitude-temperature lapse factor (ET_IF3) is constructed. The dry adiabatic lapse rate is used to correct the temperature variation with altitude at 2 meters, reflecting the thermodynamic influence of local thermal conditions on the wind field. The formula is as follows:

[0028]

[0029] in, This indicates the temperature at a height of 2 meters. Indicates the dry adiabatic lapse rate. This represents the elevation of the terrain; divide by 1000 to convert to kilometers. Represents an exponential function;

[0030] 3.5: Construct an altitude-pressure correction factor (EP_IF4). Based on the hydrostatic equation and the pressure-altitude formula, it uses atmospheric elevation to correct for the exponential decay of surface air pressure with altitude, reflecting the dynamic driving effect of the pressure gradient force on the wind field. The formula is as follows:

[0031]

[0032] in, Indicates ground air pressure. Indicates the elevation of the terrain. Indicates atmospheric elevation;

[0033] 3.6: Constructing the radiation-slope modulation factor (RS_IF5): Based on the topographic radiation modulation mechanism and Lambert's cosine law, the modulation effect of slope on solar shortwave radiation is utilized to quantify the driving effect of slope on local thermal circulation. The formula is as follows:

[0034]

[0035] in, Represents net shortwave radiation at the Earth's surface. Indicates the radiative modulation coefficient;

[0036] 3.7: Standardize all interaction features to ensure that the feature values ​​are within a reasonable range, avoid the impact of extreme values ​​on model training, and concatenate them into an interaction feature matrix.

[0037] In step 1.3, the spatiotemporal features of the data are extracted from the dataset processed in steps 1.1 and 1.2 by introducing a ResNet branch architecture with a multi-scale feature fusion module (MSFF). Specifically, this includes the following steps:

[0038] 4.1: Spatial features are extracted using the ResNet architecture, and a Multi-Scale Feature Fusion (MSFF) module is introduced. This module captures features at different spatial scales through four parallel convolutional branches, achieving collaborative modeling. The formula is as follows:

[0039]

[0040]

[0041]

[0042]

[0043] in, For the input feature map, Indicates batch size. Indicates the number of input channels. and These represent the height and width of the input space grid, respectively. Indicates the kernel size as Two-dimensional convolution operation, Indicates the first The output feature maps of each scale branch, with the number of output channels of all branches being uniformly set. Maintain spatial dimensions consistent with the input;

[0044] 4.2: Channel-dimensional concatenation is performed on the convolutional output feature maps of the four scales. Batch normalization and activation functions are then used to achieve feature fusion and improve training stability. The formula is as follows:

[0045]

[0046]

[0047] in, This represents a concatenation operation along the channel dimension, resulting in a feature map with dimensions of [dimensionality missing]. , This indicates batch normalization, used to eliminate feature distribution bias and accelerate model convergence. This represents a linear rectified activation function, and introduces nonlinear feature representation. This represents the output feature map after multi-scale fusion, with the same dimensions as the concatenated version.

[0048] 4.3: To enhance feature representation capabilities, multi-scale fused features are input into residual blocks to extract deep features. Channel and spatial weights are adaptively adjusted using a Convolutional Block Attention (CBAM) mechanism. Finally, the original features are fused through residual connections.

[0049] 4.4: A three-layer fully connected network is used to map the time parameters into a high-dimensional time embedding vector, enhancing the expressive power of time dependencies. The formula is as follows:

[0050]

[0051]

[0052]

[0053] in, Input a vector of time parameters, including year, month, day, hour, and seasonal features. and This represents the weights and biases of the first fully connected layer network. and This represents the weights and biases of the second-layer fully connected network. and This represents the weights and biases of the third-layer fully connected network. and This represents the output of the hidden layer, used to progressively increase the feature dimension. This is the final temporal embedding vector, with dimensions consistent with the number of channels in the spatial features.

[0054] In step 1.4, a Liquid Time-Constant Networks branch architecture with a Feature-Aware Adaptive Memory Decay Module (FAMD) is used to mine the time-series dynamic patterns of the data from the dataset processed in steps 1.1 and 1.2. Specifically, this includes the following steps:

[0055] 5.1: Modeling the continuous dynamic evolution of time series based on Liquid Time-Constant Networks, using ordinary differential equations (ODEs) to describe the evolution of hidden states over time, and solving numerically using the Runge-Kutta fourth-order (RK4) method, as shown in the following formula:

[0056]

[0057]

[0058]

[0059]

[0060]

[0061] in, Indicates time step The hidden state tensor, For time steps Input temporal features, To parameterize the ODE function, output the time derivative of the hidden state. Indicates the integration step size. , , , These are the four intermediate slope terms of the RK4 method. For time steps The updated hidden state;

[0062] 5.2: Introducing a Feature-Aware Adaptive Memory Decay (FAMD) mechanism to dynamically adjust the weights of historical information. By calculating adaptive decay gating weights, the hidden states are selectively forgotten and retained, enhancing the model's responsiveness to sudden changes in wind speed. The formula is as follows:

[0063]

[0064]

[0065] in, This represents the adaptive decay gate weight. It is the sigmoid activation function. and Let represent the weight matrix and bias vector of the decay gate, respectively. This represents the concatenation operation between the hidden state and the input features along the feature dimension. This is an element-wise multiplication method that implements element-wise modulation of the hidden state by the gate weights. This represents the hidden state after adaptive decay;

[0066] 5.3: Construct an Ordinary Differential Equation (ODE) function to model the evolution of the hidden state. The decayed hidden state, current input, and time are concatenated and input into a multilayer perceptron. The time derivative of the hidden state is calculated using the following formula:

[0067]

[0068]

[0069] in, This is the concatenated input vector. , and These represent the weight matrices of a three-layer fully connected network. , and These represent the corresponding bias vectors. This is a dropout regularization operation to prevent the model from overfitting. The time derivative of the hidden state;

[0070] 5.4: After four time steps of sequence iteration, the final hidden state is reshaped into a spatial feature map format, and the temporal feature representation is output as follows:

[0071]

[0072] in, This is the hidden state at the last time step. This indicates a tensor reshaping operation. This represents the temporal characteristics of the final output, with the number of channels being... =240, matching the spatial feature dimension.

[0073] In step 1.5, the features of the dual-branch model (MSFF-ResNet and FAMD-Liquid Time-Constant Networks) are spliced ​​using a gated fusion mechanism (GatedFusion), and the fusion weights are adaptively calculated before being input into the MLP to achieve prediction. Specifically, this includes the following steps:

[0074] 6.1: Adaptively fuse spatial and temporal features through a gated fusion mechanism (GatedFusion);

[0075] 6.2: Input the fused features into the MultiTaskHead, and output the point prediction, interval prediction and probability prediction results respectively through a shared feature extraction layer and three independent task branches.

[0076] In step 1.6, during the model training phase, using the datasets processed in steps 1.1 and 1.2, an adaptive composite loss function is designed for the dual-branch models (MSFF-ResNet and FAMD-Liquid Time-ConstantNetworks) constructed in steps 1.3 and 1.4. This is combined with a gradient pruning strategy to optimize model training. Specifically, the steps include the following:

[0077] 7.1: Design a composite loss function that combines point prediction loss, interval prediction loss, and probabilistic prediction loss. The formula is as follows:

[0078]

[0079] in, Indicates point prediction loss. Let represent the probability score loss for consecutive sorting. Indicates quantile loss. , as well as These represent the loss weights obtained by optimizing the search using Bayesian methods.

[0080] 7.2: The model parameters are updated using the Adaptive Moment Estimator with Weight Decay (AdamW) optimizer, and cosine annealing learning rate scheduling is adopted, combined with gradient pruning strategy to optimize model training.

[0081] In step 1.7, the test set data is used to output the prediction results through a trained model (PIMFNet) that uses physical information-guided feature construction and Bayesian optimization of the composite loss function.

[0082] Meanwhile, the present invention proposes an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the computer program is executed, it implements the steps of the method described in the present invention.

[0083] Furthermore, the present invention proposes a computer-readable storage medium having a computer program stored thereon, the computer program being configured to implement the steps of the method described in the present invention when invoked by a processor.

[0084] Finally, the present invention provides a computer program product comprising a computer program / instructions that, when executed by a processor, implement the steps of the method described in the present invention.

[0085] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0086] 1. This invention addresses the technical bottlenecks of traditional numerical weather prediction (NWP) models, which struggle to quantify prediction uncertainties and single deep learning models are unable to simultaneously capture multi-scale terrain features and non-stationary temporal dynamics, by constructing a probabilistic wind speed prediction framework guided by physical information and fused with multi-scale features. This method not only significantly improves the accuracy of wind speed point prediction under complex terrain but also outputs probability distributions and prediction intervals, providing comprehensive decision support for wind power scheduling, risk assessment, and energy system optimization.

[0087] 2. This invention addresses the shortcomings of existing models in terms of insufficient physical interpretability and poor adaptability to complex terrain by constructing physically-guided interactive features. Traditional deep learning models often employ data-driven end-to-end learning methods, lacking explicit modeling of atmospheric physical processes, resulting in insufficient generalization ability under complex terrain and extreme weather conditions. This invention innovatively constructs five physically interactive features, embedding atmospheric dynamic constraints into the model input layer. This allows the model to follow physical laws while learning data-driven patterns, significantly improving the prediction accuracy and physical consistency of the model under complex terrain conditions.

[0088] 3. This invention introduces a multi-scale feature fusion module (MSFF) and a feature-aware adaptive memory decay module (FAMD). By using parallel multi-scale convolutional kernels to extract multi-level spatial features from local to regional levels, it solves the problem that single-scale convolution cannot simultaneously capture the influence of terrain at different spatial scales. Simultaneously, it dynamically adjusts the weights of historical information in the Liquid Time-Constant Networks, addressing the slow response of traditional recurrent neural networks to strong wind non-stationarity. Existing methods have significant limitations in modeling the spatial heterogeneity of wind speed under complex terrain, lacking adaptability in the forgetting and retention strategies of historical information. When wind speed changes abruptly, the model struggles to quickly adjust memory weights, leading to prediction lag. This invention innovatively designs four parallel convolutional branches, enabling the model to adaptively fuse spatial information at different scales, enhancing its representation ability for complex terrain. Adaptive memory decay gating is introduced into the ordinary differential equation (ODE) of the Liquid Time-Constant Networks to dynamically evaluate the importance of historical information, significantly improving the model's ability to model non-stationary time series.

[0089] 4. This invention constructs a multi-task learning framework that simultaneously outputs point predictions, probability distributions, and prediction intervals, achieving comprehensive uncertainty quantification of wind speed. This solves the problem that traditional models only provide deterministic predictions and cannot meet risk assessment needs. Furthermore, it designs a Bayesian-optimized composite loss function to automatically balance multiple task objectives. Existing models mostly focus on point prediction tasks, outputting a single predicted value and using fixed-weighted loss functions, making it difficult to achieve optimal balance. This invention innovatively designs three parallel task branches and uses a Bayesian optimization algorithm to automatically search for loss weights, achieving adaptive balance of multiple task objectives. This provides technical support for efficient wind energy resource development, precise wind power scheduling, and energy system optimization. Attached Figure Description

[0090] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used together with the embodiments of the invention to explain the invention and do not constitute a limitation thereof.

[0091] Figure 1 This invention provides a system framework diagram for a probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion. Detailed Implementation

[0092] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. Of course, the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0093] Example 1: See Figure 1 As shown, this embodiment provides a probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion, including the following steps:

[0094] (1-1) The meteorological parameters, topographic data of the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 3 (ASTER GDEM V3) and time parameters of the THORPEX Interactive Global Encyclopedia (TIGGE) project were preprocessed. The dataset was constructed through missing value processing, outlier removal by 3σ rule, spatiotemporal linear interpolation and standardization, and the dataset was divided into training set, validation set and test set according to year.

[0095] (1-2) Based on the atmospheric boundary layer dynamics theory and topographic radiation modulation mechanism, physical information-guided interactive features are constructed using meteorological parameters and topographic parameters;

[0096] (1-3) Integrating a multi-scale feature fusion module (MSFF) into the ResNet architecture enables collaborative modeling of information at different spatial scales;

[0097] (1-4) In the Liquid Time-Constant Networks architecture, a Feature-Aware Adaptive Memory Decay Module (FAMD) is introduced to capture the nonlinear temporal evolution of wind speed using an adaptive memory gating mechanism;

[0098] (1-5) The dual-branch features are spliced ​​together through the gated fusion mechanism, and the fusion weights are adaptively calculated and then input into the MLP to achieve prediction;

[0099] (1-6) Design a composite loss function, automatically search for loss weights through Bayesian optimization, and combine it with an adaptive moment estimation (AdamW) optimizer with weight decay, cosine annealing learning rate scheduling and gradient pruning strategy to optimize model training.

[0100] (1-7) Use the trained model (PIMFNet) based on physical information-guided feature construction and Bayesian optimization of composite loss function to make predictions on the test set, and output the point prediction values, probability distribution and prediction interval.

[0101] The data preprocessing operation in step (1-1) includes the following steps:

[0102] (2-1) Missing values ​​were processed for meteorological parameters, topographic data of the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 3 (ASTER GDEM V3), and time parameters of the THORPEX Interactive Global Encyclopedia (TIGGE) project. Linear interpolation was used in the time dimension, and linear interpolation was first performed in the spatial dimension, followed by nearest neighbor interpolation.

[0103] (2-2) The 3σ rule was used to remove outliers that deviated from the mean by more than 3 times the standard deviation, and all data were uniformly interpolated to a 0.25°×0.25° grid to match the four time points of 06:00, 12:00, 18:00 and 24:00 every day;

[0104] (2-3) Standardize the meteorological and topographic parameters, and finally divide them into training set, validation set and test set according to year.

[0105] Step (1-2) constructs five physical interaction features based on the meteorological and terrain parameters in the dataset processed in step (1-1), specifically including the following steps:

[0106] (3-1) Based on the atmospheric boundary layer dynamics theory and topographic radiation modulation mechanism, five physical information-guided interactive features are constructed using meteorological parameters and topographic parameters;

[0107] (3-2) Construct the wind speed-topographic relief interaction factor (WR_IF1). Based on the drag effect of topographic roughness on wind speed, the dynamic modulation effect of surface roughness on near-surface wind speed is quantified using the relief degree. The formula is as follows:

[0108]

[0109] in, This indicates the wind speed at a height of 10 meters. It represents the degree of terrain relief, reflecting the surface roughness characteristics. The standard deviation of the global terrain relief is used for normalization. Indicates the terrain drag coefficient;

[0110] (3-3) Based on the dynamic modulation effect of terrain slope on wind speed, a wind speed-slope coupling factor (WS_IF2) is constructed. The slope is used to quantify the influence of terrain inclination on wind acceleration or deceleration. The formula is as follows:

[0111]

[0112] in, Indicates the slope angle of the terrain. The tangent value represents the slope, reflecting the degree of terrain inclination. The effective projection coefficient of slope wind. This represents the slope modulation coefficient, whose value is based on the slope wind speed enhancement effect in topographic dynamics theory.

[0113] (3-4) Based on the theory of atmospheric temperature vertical lapse rate, an altitude-temperature lapse factor (ET_IF3) is constructed. The dry adiabatic lapse rate is used to correct the temperature variation with altitude at 2 meters, reflecting the thermodynamic influence of local thermal conditions on the wind field. The formula is as follows:

[0114]

[0115] in, This indicates the temperature at a height of 2 meters. Indicates the dry adiabatic lapse rate. This represents the elevation of the terrain; divide by 1000 to convert to kilometers. Represents an exponential function;

[0116] (3-5) Construct an altitude-pressure correction factor (EP_IF4). Based on the static equation and the pressure-altitude formula, use atmospheric elevation to correct the exponential decay of surface air pressure with altitude, reflecting the dynamic driving effect of the pressure gradient force on the wind field. The formula is as follows:

[0117]

[0118] in, Indicates ground air pressure. Indicates the elevation of the terrain. Indicates atmospheric elevation;

[0119] (3-6) Construct the radiation-slope modulation factor (RS_IF5). Based on the topographic radiation modulation mechanism and Lambert's cosine law, utilize the modulation effect of slope on solar shortwave radiation to quantify the driving effect of slope on local thermal circulation. The formula is as follows:

[0120]

[0121] in, Represents net shortwave radiation at the Earth's surface. Indicates the radiative modulation coefficient;

[0122] (3-7) Standardize all interaction features to ensure that the feature values ​​are within a reasonable range, avoid the impact of extreme values ​​on model training, and concatenate them into an interaction feature matrix.

[0123] In steps (1-3), a ResNet branch architecture with a multi-scale feature fusion module (MSFF) is introduced to extract the spatiotemporal features of the data from the dataset processed in steps (1-1) and (1-2). Specifically, this includes the following steps:

[0124] (4-1) Spatial features are extracted using the ResNet architecture, and a multi-scale feature fusion module (MSFF) is introduced. Features at different spatial scales are captured through four parallel convolutional branches to achieve collaborative modeling. The formula is as follows:

[0125]

[0126]

[0127]

[0128]

[0129] in, For the input feature map, Indicates batch size. Indicates the number of input channels. and These represent the height and width of the input space grid, respectively. Indicates the kernel size as Two-dimensional convolution operation, Indicates the first The output feature maps of each scale branch, with the number of output channels of all branches being uniformly set. Maintain spatial dimensions consistent with the input;

[0130] (4-2) The convolutional output feature maps of the four scales are concatenated along the channel dimension. Batch normalization and activation functions are used to achieve feature fusion and improve training stability. The formula is as follows:

[0131]

[0132]

[0133] in, This represents a concatenation operation along the channel dimension, resulting in a feature map with dimensions of [dimensionality missing]. , This indicates batch normalization, used to eliminate feature distribution bias and accelerate model convergence. This represents a linear rectified activation function, and introduces nonlinear feature representation. This represents the output feature map after multi-scale fusion, with the same dimensions as the concatenated version.

[0134] (4-3) To enhance feature representation, multi-scale fused features are input into residual blocks to extract deep features, and channel and spatial weights are adaptively adjusted through a convolutional block attention mechanism (CBAM). Finally, the original features are fused through residual connections.

[0135] (4-4) A three-layer fully connected network is used to map the time parameters into a high-dimensional time embedding vector, which enhances the ability to express time dependencies. The formula is as follows:

[0136]

[0137]

[0138]

[0139] in, Input a vector of time parameters, including year, month, day, hour, and seasonal features. and This represents the weights and biases of the first fully connected layer network. and This represents the weights and biases of the second-layer fully connected network. and This represents the weights and biases of the third-layer fully connected network. and This represents the output of the hidden layer, used to progressively increase the feature dimension. This is the final temporal embedding vector, with dimensions consistent with the number of channels in the spatial features.

[0140] In steps (1-4), a Liquid Time-Constant Networks branch architecture with a Feature-Aware Adaptive Memory Decay Module (FAMD) is used to mine the time-series dynamic patterns of the data from the dataset processed in steps (1-1) and (1-2). Specifically, this includes the following steps:

[0141] (5-1) The continuous dynamic evolution of time series is modeled based on Liquid Time-Constant Networks. Ordinary Differential Equations (ODEs) are used to describe the evolution of the hidden states over time. The Runge-Kutta fourth-order (RK4) method is used for numerical integration to solve the problem. The formula is as follows:

[0142]

[0143]

[0144]

[0145]

[0146]

[0147] in, Indicates time step The hidden state tensor, For time steps Input temporal features, To parameterize the ODE function, output the time derivative of the hidden state. Indicates the integration step size. , , , These are the four intermediate slope terms of the RK4 method. For time steps The updated hidden state;

[0148] (5-2) Introducing a Feature-Aware Adaptive Memory Decay (FAMD) mechanism to dynamically adjust the weights of historical information, and selectively forgetting and retaining hidden states by calculating adaptive decay gate weights, thereby enhancing the model's ability to respond to sudden changes in wind speed. The formula is as follows:

[0149]

[0150]

[0151] in, This represents the adaptive decay gate weight. It is the sigmoid activation function. and Let represent the weight matrix and bias vector of the decay gate, respectively. This represents the concatenation operation between the hidden state and the input features along the feature dimension. This is an element-wise multiplication method that implements element-wise modulation of the hidden state by the gate weights. This represents the hidden state after adaptive decay;

[0152] (5-3) Construct an ordinary differential equation (ODE) function to model the evolution of the hidden state. The decayed hidden state, the current input, and the time are embedded and concatenated before being input into the multilayer perceptron. The time derivative of the hidden state is calculated, as shown in the following formula:

[0153]

[0154]

[0155] in, This is the concatenated input vector. , and These represent the weight matrices of a three-layer fully connected network. , and These represent the corresponding bias vectors. This is a dropout regularization operation to prevent the model from overfitting. The time derivative of the hidden state;

[0156] (5-4) After four time steps of sequence iteration, the final hidden state is reshaped into a spatial feature map format, and the temporal feature representation is output as follows:

[0157]

[0158] in, This is the hidden state at the last time step. This indicates a tensor reshaping operation. This represents the temporal characteristics of the final output, with the number of channels being... =240, matching the spatial feature dimension.

[0159] In steps (1-5), the features of the dual-branch models (MSFF-ResNet and FAMD-Liquid Time-Constant Networks) are spliced ​​using a gated fusion mechanism (GatedFusion), and the fusion weights are adaptively calculated before being input into the MLP to achieve prediction. Specifically, the steps are as follows:

[0160] (6-1) Adaptively fuse spatial and temporal features through a gated fusion mechanism (GatedFusion);

[0161] (6-2) Input the fused features into the MultiTaskHead, and output the point prediction, interval prediction and probability prediction results respectively through the shared feature extraction layer and three independent task branches.

[0162] In steps (1-6), during the model training phase, using the datasets processed in steps (1-1) and (1-2), an adaptive composite loss function is designed for the dual-branch models (MSFF-ResNet and FAMD-Liquid Time-ConstantNetworks) constructed in steps (1-3) and (1-4), and the model training is optimized by combining a gradient pruning strategy. Specifically, the steps include the following:

[0163] (7-1) Design a composite loss function that combines point prediction loss, interval prediction loss, and probability prediction loss, as shown in the following formula:

[0164]

[0165] in, Indicates point prediction loss. Let represent the probability score loss for consecutive sorting. Indicates quantile loss. , as well as These represent the loss weights obtained by optimizing the search using Bayesian methods.

[0166] (7-2) The model parameters are updated using the Adaptive Moment Estimator with Weight Decay (AdamW) optimizer, and the learning rate is scheduled using cosine annealing, combined with the gradient pruning strategy to optimize model training.

[0167] In steps (1-7), test set data is used to output prediction results through a trained model (PIMFNet) that uses physical information-guided feature construction and Bayesian optimization of the composite loss function.

[0168] On the same dataset, namely the eight meteorological parameters most correlated with wind speed data in ERA5 hourly data on single levels selected from the THORPEX Interactive Global Collection (TIGGE) project, the three topographic parameters most correlated with wind speed data in ERA5 hourly data on single levels selected from the ASTER GDEM V3 topographic data, two additional meteorological parameters (2m temperature and surface air pressure), five constructed physical interaction features (WR_IF1, WS_IF2, ET_IF3, EP_IF4, RS_IF5), and five time parameters, covering four time points each day from January 1, 2024 to December 31, 2024 (UTC 06:00, 12:00, 18:00, and 24:00). The point prediction performance of the proposed model based on physical information-guided feature construction and Bayesian optimization of the composite loss function, as well as existing models such as the Long Short-Term Memory (LSTM) network model, the European Centre for Medium-Range Weather Forecasts (ECMWF) model, the Transformer model, and the Gated Recurrent Unit (GRU) model, was evaluated. Seven performance indicators from the field of wind speed prediction research (i.e., root mean square error (RMSE), mean absolute error (MAE), correlation coefficient (R), relative root mean square error (rRMSE), relative mean absolute error (rMAE), percentage of samples with wind speed absolute error not exceeding 1 m / s (FA), and mean absolute percentage error (MAPE)) were used to automatically evaluate the quality of the model.

[0169] Table 1. Comparison of point prediction results between the method of this invention and other methods.

[0170]

[0171] Experiments show that the PIMFNet model proposed in this embodiment, which uses physically-guided feature construction and a Bayesian-optimized composite loss function, achieves more accurate deterministic wind speed prediction compared to benchmark methods. Specifically, the model in this embodiment effectively captures the physical evolution and nonlinear dynamic characteristics of wind speed sequences through a physically-guided feature construction approach combined with a Bayesian-optimized composite loss function. Furthermore, by constructing a dual-branch architecture of MSFF-ResNet and FAMD-LiquidTime-Constant Networks, the model exhibits excellent predictive capabilities under various geographical and meteorological conditions, surpassing these baseline methods. Specifically, for FA, the proposed method achieves at least a 9.27% ​​performance improvement; for RMSE, it reduces the error by at least 19.29%; for MAE, it reduces the error by at least 18.25%; for rRMSE, it reduces the error by at least 17.41%; for rMAE, it reduces the error by at least 15.11%; for MAPE, it reduces the error by at least 30.97%; and for R, it improves the performance by at least 2.15%. These results demonstrate the competitiveness of the proposed method.

[0172] Example 2: To compare the performance of the model of this invention with existing probabilistic wind speed prediction models for interval prediction, the same dataset as in Example 1 was selected for testing. Specifically, the dataset consisted of eight meteorological parameters with the highest correlation to wind speed data in ERA5 hourly data selected from the THORPEX Interactive Global Collection (TIGGE) project; three terrain parameters with the highest correlation to wind speed data in ERA5 hourly data selected from the ASTER GDEM V3 topographic data; two additional meteorological parameters (2m temperature and surface air pressure); and five constructed physical interaction features (WR_IF1, WS_IF2, ET_IF3, EP_IF4). The model consists of RS_IF5 and five time parameters, covering four time points each day from January 1, 2024 to December 31, 2024 (UTC 06:00, 12:00, 18:00, and 24:00). The proposed model (PIMFNet) based on physical information-guided feature construction and Bayesian optimization of the composite loss function, along with existing Long Short-Term Memory (LSTM) models, Deep Autoregressive (DeepAR) models, Transformer models, and Autoregressive Integrated Moving Average (ARIMA) models, are evaluated. Four performance metrics from the wind speed forecasting research field (i.e., coverage probability CP, mean width percentage MWP, mean coverage rate MC, and prediction interval normalized mean width PINAW) are used to automatically assess the model quality.

[0173] Table 2 Comparison of Interval Prediction Results between the Method of this Invention and Other Methods

[0174]

[0175] Experiments show that the model proposed in this embodiment, based on physical information-guided feature construction and Bayesian optimized composite loss function (PIMFNet), achieves more accurate wind speed interval prediction compared to the baseline method. Specifically, for CP, the method in this embodiment maintains a relatively high error level of 0.9763, ensuring the reliability of interval prediction coverage; for MWP, the method in this embodiment can reduce the error by at least 16.68%; for MC, the method in this embodiment can reduce the error by at least 18.19%; and for PINAW, the method in this embodiment can reduce the error by at least 16.73%. These results demonstrate the competitiveness of the proposed method in this embodiment.

[0176] Example 3: To compare the performance of the model of this invention with existing probabilistic wind speed prediction models in 2024 and various seasons for probabilistic prediction, the same dataset as in Example 1 was selected for testing. This dataset consisted of eight meteorological parameters most correlated with wind speed data in ERA5 hourly data selected from the THORPEX Interactive Global Collection (TIGGE) project; three topographic parameters most correlated with wind speed data in ERA5 hourly data selected from the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 3 (ASTER GDEM V3) topographic data; two additional meteorological parameters (2m temperature and surface air pressure); five constructed physical interaction features (WR_IF1, WS_IF2, ET_IF3, EP_IF4, RS_IF5); and five time parameters, covering four time points (UTC) each day from January 1, 2024 to December 31, 2024. (06:00, 12:00, 18:00, and 24:00). The proposed model based on physical information-guided feature construction and Bayesian optimization of the composite loss function (PIMFNet) and existing models such as Long Short-Term Memory (LSTM), Deep Autoregressive (DeepAR), Transformer, and Autoregressive Integrated Moving Average (ARIMA) are evaluated. The continuous ranking probability score (CRPS) index from the field of wind speed prediction research is used to automatically evaluate the quality of the models.

[0177] Table 3 Comparison of Probability Prediction Results between the Method of this Invention and Other Methods

[0178]

[0179] Experiments show that the model proposed in this embodiment, based on physical information-guided feature construction and Bayesian optimized composite loss function (PIMFNet), achieves more accurate wind speed probability prediction compared to the baseline method. Specifically, for the CRPS index throughout the year, the method of this embodiment can reduce the error by at least 5.37%; for the CRPS index in spring, the method can reduce the error by at least 6.13%; for the CRPS index in summer, the method can reduce the error by at least 5.82%; for the CRPS index in autumn, the method can reduce the error by at least 6.98%; and for the CRPS index in winter, the method can reduce the error by at least 8.40%. These results demonstrate the competitiveness of the proposed method.

[0180] Example 4: This example proposes an electronic system, including: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the method steps of the present invention.

[0181] Example 5: This example proposes a computer-readable storage medium storing a computer program thereon. When the computer program is executed by a processor, it implements the steps of the method described in this invention, which will not be repeated here.

[0182] Example 6: This example proposes a computer program product, including a computer program / instructions. When the computer program / instructions are executed by a processor, they implement the steps of the method described in this invention, which will not be repeated here.

[0183] It should be noted that the processing flow of embodiments 4-6 corresponds to the specific steps of the method provided in embodiment 1 of the present invention, and has the corresponding functional modules and beneficial effects of the method. Technical details not described in detail in this embodiment can be found in the method provided in embodiment 1 of the present invention.

[0184] The program code used to implement the methods of this application may be written in any combination of one or more programming languages. This program code may be provided to a processor or controller of a general-purpose computer, special-purpose computer, or other programmable data processing device, such that when executed by the processor or controller, the functions / operations specified in the flowcharts and / or block diagrams are implemented. The program code may be executed entirely on a machine, partially on a machine, as a standalone software package partially on a machine and partially on a remote machine, or entirely on a remote machine or server.

[0185] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion, characterized in that, Includes the following steps: S1: Preprocess the meteorological parameters, topographic data of ASTER GDEM V3 (Advanced Spaceborne Thermal Emission and Reflection Radiometer) from the THORPEX interactive global dataset TIGGE project, and time parameters. Construct the dataset through missing value handling, outlier removal using the 3σ rule, spatiotemporal linear interpolation, and standardization steps, and divide it into training, validation, and test sets according to year. S2: Based on atmospheric boundary layer dynamics theory and topographic radiation modulation mechanism, physical information-guided interactive features are constructed using meteorological and topographic parameters; S3: Integrating the multi-scale feature fusion module MSFF into the ResNet residual network architecture to achieve collaborative modeling of information at different spatial scales; S4: In the Liquid Time-Constant Networks architecture, a Feature Aware Adaptive Memory Decay (FAMD) module is introduced to capture the nonlinear temporal evolution of wind speed using an adaptive memory gating mechanism. S5: The dual-branch features are concatenated through the gated fusion mechanism GatedFusion, and the fusion weights are adaptively calculated before being input into the MLP to achieve prediction. S6: Design a composite loss function, automatically search for loss weights through Bayesian optimization, and combine it with the AdamW optimizer with weight decay, cosine annealing learning rate scheduling and gradient pruning strategy to optimize model training. S7: Use the trained PIMFNet model, which is based on physical information-guided feature construction and Bayesian optimization of the composite loss function, to make predictions on the test set, and output the predicted point values, probability distributions, and prediction intervals.

2. The probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion according to claim 1, characterized in that, S1 includes the following steps: S11: Missing values ​​are processed for meteorological parameters, topographic data of ASTER GDEM V3 (Advanced Spaceborne Thermal Emission and Reflection Radiometer) and time parameters of the THORPEX interactive global collection TIGGE project. Linear interpolation is used in the time dimension, and linear interpolation is first performed in the spatial dimension, followed by nearest neighbor interpolation. S12: Use the 3σ rule to remove outliers that deviate from the mean by more than 3 times the standard deviation, and interpolate all data to a 0.25°×0.25° grid to match the four time points of 06:00, 12:00, 18:00 and 24:00 each day; S13: Standardize the meteorological and topographic parameters, and finally divide them into training, validation, and test sets according to the year.

3. The probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion according to claim 1, characterized in that, S2 includes the following steps: S21: Based on atmospheric boundary layer dynamics theory and topographic radiation modulation mechanism, five physical information-guided interactive features are constructed using meteorological and topographic parameters; S22: Construct the wind speed-topographic relief interaction factor WR_IF1. Based on the drag effect of topographic roughness on wind speed, the dynamic modulation effect of surface roughness on near-surface wind speed is quantified using the terrain relief degree. The formula is as follows: ； in, This indicates the wind speed at a height of 10 meters. Relief indicates the degree of terrain relief, reflecting the surface roughness characteristics. The standard deviation of the global terrain relief is used for normalization. Indicates the terrain drag coefficient; S23: Based on the dynamic modulation effect of terrain slope on wind speed, a wind speed-slope coupling factor WS_IF2 is constructed. The slope is used to quantify the influence of terrain inclination on wind acceleration or deceleration. The formula is as follows: ； in, Indicates the slope angle of the terrain. The tangent value represents the slope, reflecting the degree of terrain inclination. The effective projection coefficient of slope wind. This represents the slope modulation coefficient, whose value is based on the slope wind speed enhancement effect in topographic dynamics theory. S24: Based on the theory of atmospheric temperature vertical lapse rate, an altitude-temperature lapse factor ET_IF3 is constructed. The dry adiabatic lapse rate is used to correct the temperature variation with altitude at 2 meters, reflecting the thermodynamic influence of local thermal conditions on the wind field. The formula is as follows: ； in, This indicates the temperature at a height of 2 meters. Indicates the dry adiabatic lapse rate. This represents the elevation of the terrain; divide by 1000 to convert to kilometers. Represents an exponential function; S25: Construct the altitude-pressure correction factor EP_IF4, based on the hydrostatic equation and the pressure-altitude formula. It uses atmospheric elevation to correct for the exponential decay of surface air pressure with altitude, reflecting the dynamic driving effect of the pressure gradient force on the wind field. The formula is as follows: ； in, Indicates ground air pressure. Indicates the elevation of the terrain. Indicates atmospheric elevation; S26: Construct the radiation-slope modulation factor RS_IF5. Based on the topographic radiation modulation mechanism and Lambert's cosine law, utilize the modulation effect of slope on solar shortwave radiation to quantify the driving effect of slope on local thermal circulation. The formula is as follows: ； in, Represents net shortwave radiation at the Earth's surface. Indicates the radiative modulation coefficient; S27: Standardize all interaction features to ensure that the feature values ​​are within a reasonable range, avoid the impact of extreme values ​​on model training, and concatenate them into an interaction feature matrix.

4. The probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion according to claim 1, characterized in that, S3 includes the following steps: S31: Spatial features are extracted using the ResNet residual network architecture, and a multi-scale feature fusion module (MSFF) is introduced. Four parallel convolutional branches capture features at different spatial scales to achieve collaborative modeling, as shown in the following formula: ； ； ； ； in, For the input feature map, Indicates batch size. Indicates the number of input channels. and These represent the height and width of the input space grid, respectively. Indicates the kernel size as Two-dimensional convolution operation, Indicates the first The output feature maps of each scale branch, with the number of output channels of all branches being uniformly set. Maintain spatial dimensions consistent with the input; S32: Concatenate the convolutional output feature maps of the four scales along the channel dimension, and combine batch normalization and activation functions to achieve feature fusion and improve training stability. The formula is as follows: ； ； in, This represents a concatenation operation along the channel dimension, resulting in a feature map with dimensions of [dimensionality missing]. , This indicates batch normalization, used to eliminate feature distribution bias and accelerate model convergence. This represents a linear rectified activation function, and introduces nonlinear feature representation. This represents the output feature map after multi-scale fusion, with the same dimensions as the concatenated version. S33: To enhance feature representation capabilities, multi-scale fused features are input into residual blocks to extract deep features. Channel and spatial weights are adaptively adjusted using the Convolutional Block Attention (CBAM) mechanism, and finally, the original features are fused through residual connections. S34: A three-layer fully connected network maps time parameters into high-dimensional time embedding vectors, enhancing the expressive power of time dependencies. The formula is as follows: ； ； ； in, Input a vector of time parameters, including year, month, day, hour, and seasonal features. and This represents the weights and biases of the first fully connected layer network. and This represents the weights and biases of the second-layer fully connected network. and This represents the weights and biases of the third-layer fully connected network. and This represents the output of the hidden layer, used to progressively increase the feature dimension. This is the final temporal embedding vector, with dimensions consistent with the number of channels in the spatial features.

5. The probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion according to claim 1, characterized in that, S4 includes the following steps: S41: Modeling the continuous dynamic evolution of time series based on Liquid Time-Constant Networks, using ordinary differential equations (ODEs) to describe the evolution of hidden states over time, and solving numerically using the Runge-Kutta fourth-order RK4 method, as shown in the following formula: ； ； ； ； ； in, Indicates time step The hidden state tensor, For time steps Input temporal features, To parameterize the ODE function, output the time derivative of the hidden state. Indicates the integration step size. , , , These are the four intermediate slope terms of the RK4 method. For time steps The updated hidden state; S42: Introducing the Feature Aware Adaptive Memory Decay (FAMD) mechanism to dynamically adjust historical information weights. By calculating adaptive decay gate weights, hidden states are selectively forgotten and retained, enhancing the model's responsiveness to sudden changes in wind speed. The formula is as follows: ； ； in, This represents the adaptive decay gate weight. It is the sigmoid activation function. and Let represent the weight matrix and bias vector of the decay gate, respectively. This represents the concatenation operation between the hidden state and the input features along the feature dimension. This is an element-wise multiplication method that implements element-wise modulation of the hidden state by the gate weights. This represents the hidden state after adaptive decay; S43: Construct an ordinary differential equation (ODE) function to model the evolution of the hidden state. The decayed hidden state, current input, and time are concatenated and input into a multilayer perceptron. The time derivative of the hidden state is calculated using the following formula: ； ； in, This is the concatenated input vector. , and These represent the weight matrices of a three-layer fully connected network. , and These represent the corresponding bias vectors. This is a dropout regularization operation to prevent the model from overfitting. The time derivative of the hidden state; S44: After four time steps of sequence iteration, the final hidden state is reshaped into a spatial feature map format, and the temporal feature representation is output as follows: ； in, This is the hidden state at the last time step. This indicates a tensor reshaping operation. This represents the temporal characteristics of the final output, with the number of channels being... =240, matching the spatial feature dimension.

6. The probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion according to claim 1, characterized in that, S5 includes the following steps: S51: Adaptive fusion of spatial and temporal features through the gated fusion mechanism GatedFusion; S52: Input the fused features into the MultiTaskHead, and output point prediction, interval prediction and probability prediction results respectively through a shared feature extraction layer and three independent task branches.

7. The probabilistic wind speed prediction method based on physical guidance and multi-scale feature fusion according to claim 1, characterized in that, S6 includes the following steps: S61: Design a composite loss function that combines point prediction loss, interval prediction loss, and probabilistic prediction loss, as shown in the following formula: ； in, Indicates the point prediction loss. Denotes the probability score loss for consecutive sorting. Indicates quantile loss. , as well as These represent the loss weights obtained by optimizing the search using Bayesian methods. S62: The AdamW optimizer with weighted decay adaptive moment estimation updates the model parameters, and cosine annealing learning rate scheduling is used, combined with gradient pruning strategy to optimize model training.

8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the computer program is executed, it implements the steps of the method as described in any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, The computer program is configured to implement the steps of the method according to any one of claims 1 to 7 when invoked by a processor.

10. A computer program product comprising a computer program / instructions, characterized in that, When the computer program / instructions are executed by the processor, they implement the steps of the method according to any one of claims 1 to 7.

Images