A Downscaling Optimization Method for Mountain Wind Fields Based on High-Precision Simulation

Abstract

This invention discloses a downscaling optimization method for mountain wind fields based on high-precision simulation, belonging to the field of meteorological numerical simulation technology. The method includes acquiring topographic elevation data and meteorological monitoring data of the target mountain area; identifying abrupt slope changes and constructing a spiral sampling path; calculating aspect variation values ​​to determine the sampling point spacing and generate target sampling data; performing orthogonal decomposition on the sampling data to obtain frequency components; calculating the surface wave coefficient to construct boundary perturbation equations to obtain surface stress distribution; calculating airflow motion state based on surface stress distribution; solving the vorticity equation to obtain vorticity characteristics of the leeward region and calculating airflow motion correction parameters; and performing downscaling iterative calculations with the meteorological monitoring data to generate target area wind field data with a spatial resolution of hundreds of meters. This invention improves the simulation accuracy of local wind fields under complex terrain conditions.

Description

Technical Field

[0001] This invention relates to the field of meteorological numerical simulation technology, specifically to a method for downscaling optimization of mountain wind fields based on high-precision simulation. Background Technology

[0002] With the rapid development of the wind power industry, the demand for refined wind field simulation in mountainous wind farm site selection is increasing. Existing wind field simulation methods are mainly based on large-scale meteorological data, obtaining local wind field distribution through numerical calculations. Under complex terrain conditions, traditional methods struggle to accurately depict the impact of small-scale terrain changes on airflow, resulting in limited accuracy in wind field simulation.

[0003] Currently, commonly used methods for simulating wind fields in mountainous areas include statistical downscaling and dynamic downscaling. Statistical downscaling methods rely on historical observation data to establish statistical relationships, which makes them less applicable in mountainous areas where observation data is scarce. While dynamic downscaling methods consider the influence of topography, their computational grid resolution is relatively coarse, making it difficult to reflect the modulating effect of topographic features on airflow.

[0004] Existing technologies often employ fixed computational grids when dealing with airflow motion under complex terrain conditions, failing to adaptively adjust sampling density based on terrain features. This results in insufficient accuracy in airflow simulation near key terrain feature points. Furthermore, there is a lack of effective methods to characterize the leeward vortex effect commonly found in mountainous areas, making it difficult to accurately depict the motion characteristics of airflow in terrain-blocked regions. Summary of the Invention

[0005] The purpose of this invention is to provide a downscaling optimization method for mountain wind fields based on high-precision simulation. This method improves the simulation accuracy of local wind fields under complex terrain conditions through adaptive terrain sampling and airflow motion correction.

[0006] The method for downscaling and optimizing mountain wind fields based on high-precision simulation provided in this embodiment of the invention includes the following steps:

[0007] Obtain topographic elevation data and meteorological monitoring data for the target mountainous area;

[0008] Identify abrupt slope changes in terrain elevation data and construct a spiral sampling path using these abrupt slope changes as starting points.

[0009] Calculate the slope aspect variation value along the spiral sampling path, determine the sampling point spacing based on the slope aspect variation value, and generate target sampling data;

[0010] Orthogonal decomposition is performed on the target sampling data to obtain the frequency components of the surface undulation. The surface fluctuation coefficient is calculated by combining the frequency components with meteorological monitoring data.

[0011] The boundary perturbation equation is constructed using the surface wave coefficient, and the surface stress distribution is obtained by solving the boundary perturbation equation.

[0012] The airflow motion state is calculated based on the surface stress distribution, and the airflow phase transition position is extracted. A vorticity equation is constructed, and the vorticity equation is solved to obtain the vorticity intensity and range of action in the leeward region. Based on the vorticity intensity and range of action in the leeward region, the airflow motion correction parameters are calculated.

[0013] By performing downscaling iterative calculations on airflow motion correction parameters and meteorological monitoring data, target area wind field data with a spatial resolution of hundreds of meters is generated.

[0014] Furthermore, identifying abrupt slope changes in the terrain elevation data and using these abrupt changes as starting points to construct a spiral sampling path includes:

[0015] The terrain elevation data is divided into grid cells. The elevation difference between each grid point in the grid cell and its adjacent grid points in the horizontal and vertical directions is calculated. The slope value of the grid point is generated based on the elevation difference.

[0016] Grid points whose slope values ​​fall within a preset slope threshold range are statistically analyzed and marked as candidate points for slope abrupt changes.

[0017] Calculate the spacing between candidate points of slope abrupt change, divide the candidate points of slope abrupt change into multiple point sets based on the spacing, and extract the candidate point with the largest slope value from each point set as the slope abrupt change point;

[0018] Calculate the elevation change value of the slope abrupt change point along the orthogonal coordinate axis and diagonal direction, and determine the direction with the largest elevation change value as the sampling direction;

[0019] The starting angle of the spiral path is set based on the sampling direction. The spiral sampling path is constructed with the elevation change value as the center, using the slope change point as the center.

[0020] Furthermore, the slope aspect variation value is calculated along the spiral sampling path, and the sampling point spacing is determined based on the slope aspect variation value to generate target sampling data, including:

[0021] Obtain the location and elevation data of sampling points along the spiral sampling path;

[0022] Construct an observation unit consisting of three adjacent sampling points, calculate the elevation difference and horizontal distance between the central sampling point and the adjacent sampling points of the observation unit, and generate the slope aspect change value of the central sampling point;

[0023] A slope aspect change sequence is generated by sliding the observation unit along a spiral sampling path;

[0024] The difference in slope aspect change values ​​between adjacent observation units is calculated based on the slope aspect change sequence, and the central sampling point whose slope aspect change value difference exceeds the preset difference range is marked as a terrain sampling point;

[0025] Calculate the path distance between adjacent terrain sampling points, determine the sampling point spacing based on the path distance and the slope aspect change value of the terrain sampling points, regenerate the sampling point positions according to the sampling point spacing, and output the target sampling data.

[0026] Furthermore, the target sampling data is orthogonally decomposed to obtain the frequency components of the land surface undulations. The land surface fluctuation coefficient is calculated by combining the frequency components with meteorological monitoring data, including:

[0027] The target sampled data is decomposed into forward and reverse orthogonal components to obtain forward and reverse frequency components. The ratio of the forward and reverse frequency components is used to construct a frequency feature sequence.

[0028] Calculate the changing trend of the frequency feature sequence, divide the frequency feature sequence into multiple observation groups according to the inflection point of the changing trend, and calculate the degree of dispersion of the frequency features within each observation group.

[0029] Based on the degree of dispersion, the frequency clustering characteristics of the observation group are determined, and the curve of the frequency clustering characteristics changing with spatial location is used as the basis for partitioning to generate observation segments.

[0030] Meteorological monitoring data is acquired within the observation section, and a correspondence is established between the gradient changes and frequency clustering characteristics of the meteorological monitoring data to form a spatiotemporal feature sequence.

[0031] The fluctuation amplitude of frequency clustering features in the spatiotemporal feature sequence is calculated, and the weight ratio of frequency clustering features and meteorological gradient values ​​is dynamically adjusted according to the fluctuation amplitude to generate the surface fluctuation coefficient.

[0032] Furthermore, by constructing a boundary perturbation equation using the surface wave coefficient, and solving the boundary perturbation equation, the surface stress distribution is obtained, including:

[0033] An equilibrium state function is constructed based on the surface wave coefficient. The second derivative of the equilibrium state function is calculated. The second derivative is combined with the spatial gradient of the surface wave coefficient to form a boundary perturbation term.

[0034] Construct boundary perturbation equations that include dynamic and perturbation terms, where the dynamic term is determined by the time rate of change of the surface wave coefficient, and the perturbation term is determined by the spatial distribution of the boundary perturbation term;

[0035] The boundary perturbation equation is solved iteratively. In each iteration, the weight coefficients of the dynamic and perturbation terms are dynamically adjusted based on the residual value of the current solution until the residual value is less than the preset residual threshold to obtain a converged solution.

[0036] The surface stress tensor is calculated based on the converged solution, a stress propagation matrix is ​​constructed, the surface stress tensor is substituted into the stress propagation matrix, and the surface stress distribution is obtained through matrix recursion.

[0037] Furthermore, based on the surface stress distribution, the airflow motion state is calculated, and the airflow phase transition location is extracted to construct the vorticity equation, including:

[0038] Calculate the stress gradient matrix based on surface stress distribution data;

[0039] A velocity potential function is constructed based on the stress gradient matrix, and the spatial derivative of the velocity potential function is calculated to obtain the airflow velocity vector field.

[0040] The curl component is calculated for the airflow velocity vector field, and the divergence component is calculated for the stress gradient matrix. The curl component and the divergence component are combined to construct the vorticity eigenvector.

[0041] Based on the vortex eigenvectors, a vortex equation coefficient matrix is ​​constructed, and the vortex equation coefficient matrix is ​​then divided into blocks and diagonalized to obtain the vortex equation.

[0042] Furthermore, solving the vorticity equation yields the vorticity intensity and effective range in the leeward region. Based on the vorticity intensity and effective range in the leeward region, the airflow motion correction parameters are calculated, including:

[0043] Solve the vorticity equation to obtain the vorticity intensity distribution in the leeward region, and determine the boundary of the leeward region's effective range based on the gradient threshold of the vorticity intensity distribution.

[0044] Within the leeward region's effective range boundary, spatial discrete points are established, and the vorticity intensity correction coefficient at these spatial discrete points is calculated. The vorticity intensity correction coefficient is then combined with the vorticity intensity field to generate airflow motion correction parameters.

[0045] Furthermore, the airflow motion correction parameters are downscaled and iteratively calculated with meteorological monitoring data to generate wind field data for the target area with a spatial resolution of hundreds of meters, including:

[0046] The calculation region is divided into multiple layers based on the spatial resolution of the meteorological monitoring data. The calculation region is further subdivided into calculation units layer by layer. Interpolation is performed on the meteorological monitoring data to obtain the initial wind field data of the calculation unit.

[0047] Obtain the terrain height data of the calculation unit, calculate the slope change value of the terrain height data, and weight the terrain height data according to the slope change value to obtain the terrain enhancement factor;

[0048] The wind field correction is generated by multiplying the airflow motion correction parameter with the terrain enhancement factor, and the wind field correction is superimposed with the initial wind field data to obtain the corrected wind field data.

[0049] Based on the principle of divergence conservation, a conservation equation is constructed. The corrected wind field data is substituted into the mapped wind field data at the calculation node of the conservation equation, and the difference between the corrected wind field data and the mapped wind field data is calculated to obtain the wind field residual.

[0050] The wind field residuals are compared with the convergence threshold. For the calculation units that are greater than the convergence threshold, the steps of generating corrected wind field data and calculating wind field residuals are repeated until the wind field residuals of all calculation units are less than the convergence threshold.

[0051] Based on the converged wind field data, target area wind field data with a spatial resolution of hundreds of meters is generated.

[0052] One technical solution provided in this embodiment of the invention is an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the steps of the method described in any of the foregoing embodiments.

[0053] One technical solution provided in this embodiment of the invention is a computer-readable storage medium storing computer program instructions, which, when executed by a processor, implement the steps in the method described in any of the preceding claims.

[0054] This invention constructs a spiral sampling path using abrupt slope changes as the starting point, and adaptively adjusts the sampling point spacing based on slope aspect changes, thereby increasing the sampling density in terrain-featured areas and enhancing the accuracy of surface undulation characterization. It calculates the surface wave coefficient by obtaining frequency components through orthogonal decomposition, and accurately describes the impact of terrain on airflow motion by combining it with boundary perturbation equations. Based on surface stress distribution, it extracts the airflow phase transition location and solves the vorticity equation, achieving a quantitative characterization of vorticity characteristics in the leeward region and effectively depicting the airflow motion patterns under complex terrain conditions. Furthermore, it employs airflow motion correction parameters for downscaling and iterative calculations, improving the spatial resolution of local wind fields and significantly enhancing the simulation accuracy of mountain wind fields, providing reliable technical support for wind farm site selection and wind energy resource assessment. Attached Figure Description

[0055] Figure 1 A flowchart of a method for downscaling and optimizing mountain wind fields based on high-precision simulation provided in an embodiment of the present invention;

[0056] Figure 2 This is a schematic diagram of orthogonal decomposition frequency component analysis provided in an embodiment of the present invention;

[0057] Figure 3 A flowchart of wind field iterative correction based on terrain correction is provided for embodiments of the present invention. Detailed Implementation

[0058] like Figure 1 As shown, Figure 1 A flowchart of a method for downscaling and optimizing mountain wind fields based on high-precision simulation provided in an embodiment of the present invention is shown. The method includes the following steps:

[0059] Obtain topographic elevation data and meteorological monitoring data for the target mountainous area;

[0060] Identify abrupt slope changes in terrain elevation data and construct a spiral sampling path using these abrupt slope changes as starting points.

[0061] Calculate the slope aspect variation value along the spiral sampling path, determine the sampling point spacing based on the slope aspect variation value, and generate target sampling data;

[0062] Orthogonal decomposition is performed on the target sampling data to obtain the frequency components of the surface undulation. The surface fluctuation coefficient is calculated by combining the frequency components with meteorological monitoring data.

[0063] The boundary perturbation equation is constructed using the surface wave coefficient, and the surface stress distribution is obtained by solving the boundary perturbation equation.

[0064] The airflow motion state is calculated based on the surface stress distribution, and the airflow phase transition position is extracted. A vorticity equation is constructed, and the vorticity equation is solved to obtain the vorticity intensity and range of action in the leeward region. Based on the vorticity intensity and range of action in the leeward region, the airflow motion correction parameters are calculated.

[0065] By performing downscaling iterative calculations on airflow motion correction parameters and meteorological monitoring data, target area wind field data with a spatial resolution of hundreds of meters is generated.

[0066] Identifying abrupt slope changes in terrain elevation data and constructing a spiral sampling path using these abrupt changes as starting points includes:

[0067] The terrain elevation data is divided into grid cells. The elevation difference between each grid point in the grid cell and its adjacent grid points in the horizontal and vertical directions is calculated. The slope value of the grid point is generated based on the elevation difference.

[0068] Grid points whose slope values ​​fall within a preset slope threshold range are statistically analyzed and marked as candidate points for slope abrupt changes.

[0069] Calculate the spacing between candidate points of slope abrupt change, divide the candidate points of slope abrupt change into multiple point sets based on the spacing, and extract the candidate point with the largest slope value from each point set as the slope abrupt change point;

[0070] Calculate the elevation change value of the slope abrupt change point along the orthogonal coordinate axis and diagonal direction, and determine the direction with the largest elevation change value as the sampling direction;

[0071] The starting angle of the spiral path is set based on the sampling direction. The spiral sampling path is constructed with the elevation change value as the center, using the slope change point as the center.

[0072] First, acquire the terrain elevation data of the mountainous area and divide it into grid cells. These grid cells can be regular rectangular grids, typically 50 meters × 50 meters in size, constructing a grid that completely covers the target mountainous area. For each grid point, calculate the elevation difference between it and its adjacent grid points in the horizontal and vertical directions. For example, for the grid point with coordinates (i, j), calculate the elevation difference between it and the grid points with coordinates (i+1, j), (i-1, j), (i, j+1), and (i, j-1). The elevation difference can be obtained by subtracting the elevation values ​​of adjacent grid points from the current grid point's elevation value.

[0073] Based on the calculated elevation differences, a slope value is generated for each grid point. The slope value is calculated using the central difference method, dividing the elevation differences in the horizontal and vertical directions by the corresponding grid spacing to obtain the slope components in the two directions. The slope value of the grid point is obtained by calculating the square root of the sum of the squares of these two components. If the grid spacing is 50 meters, and the east-west elevation difference of a certain grid point is 30 meters and the north-south elevation difference is 40 meters, then the slope value of that point is approximately 1.0, indicating an average slope of 45 degrees.

[0074] Grid points whose slope values ​​fall within a preset slope threshold range are marked as candidate points for slope abrupt changes. The preset slope threshold range can be set from 0.5 to 2.0, corresponding to slope angles of approximately 26.5 to 63.4 degrees. When the terrain transitions from gentle to steep areas, slope values ​​change significantly; these points of change are often key locations influencing wind fields in mountainous regions. In the terrain elevation data, if a grid point has a slope value of 1.2, falling within the preset slope threshold range, it is marked as a candidate point for slope abrupt changes.

[0075] For the marked candidate points of slope abrupt change, calculate the distance between them. The distance is obtained by calculating the Euclidean distance between the two candidate points on the horizontal plane. Based on the distance, the candidate points of slope abrupt change are divided into multiple point sets. The division rule is: if the distance between two candidate points is less than a preset threshold (e.g., 100 meters), they are assigned to the same point set. From each point set, the candidate point with the largest slope value is extracted as the slope abrupt change point. For example, if a point set contains three candidate points with slope values ​​of 0.7, 1.3, and 0.9, the candidate point with a slope value of 1.3 is selected as the slope abrupt change point of that point set.

[0076] For each defined abrupt change in slope, calculate its elevation change along the orthogonal coordinate axes (east-west and north-south directions) and diagonally. The elevation change is obtained by calculating the absolute value of the elevation difference between the abrupt change point and a grid point at a specific distance in the corresponding direction. This specific distance can be set to 150 meters, and eight directions are calculated: east, west, south, north, northeast, northwest, southeast, and southwest. The direction with the largest elevation change is determined as the sampling direction. If a slope abrupt change point has the largest elevation change in the northwest direction (210 meters), then the northwest direction is determined as the sampling direction for that point.

[0077] The starting angle of the spiral path is set based on the determined sampling direction. The mapping relationship between the sampling direction and the starting angle is: East corresponds to 0 degrees, Northeast corresponds to 45 degrees, North corresponds to 90 degrees, and so on. Taking the abrupt change in slope as the center, the previously calculated elevation change value is used as the spiral sampling radius to construct a spiral sampling path. The parametric equation of the spiral path is in polar coordinate form, where the radial distance increases linearly with the angle. Starting from the initial angle, the sampling point is expanded every certain angle (e.g., 10 degrees) until one or more loops of sampling are completed.

[0078] Multiple sampling points are selected along the spiral path, with each sampling point spaced 10 degrees apart. The sampling radius is calculated using the formula r = r0 + kθ, where r represents the sampling radius, r0 is the initial radius (related to the elevation change value, such as 0.8 times the elevation change value), k is the growth coefficient (e.g., 0.5), and θ is the current angle. This generated spiral path effectively captures the influence of complex mountain terrain on the wind field. Terrain feature parameters, including elevation, slope, and aspect, are recorded at each sampling point.

[0079] This invention accurately identifies abrupt slope changes in mountainous terrain, capturing the most significant topographic features affecting wind fields. The spiral sampling path constructed based on these abrupt slope changes effectively covers wind field characteristics in areas of dramatic terrain variation, resulting in more accurate downscaling simulations of the wind field. Compared to traditional regular grid sampling methods, this method significantly improves computational efficiency, reduces invalid sampling points, and lowers computational resource consumption. By adaptively determining the sampling direction and spiral path parameters, the simulation results are enhanced to adapt to complex mountainous terrain.

[0080] Calculate the slope aspect variation value along the spiral sampling path, determine the sampling point spacing based on the slope aspect variation value, and generate target sampling data including:

[0081] Obtain the location and elevation data of sampling points along the spiral sampling path;

[0082] Construct an observation unit consisting of three adjacent sampling points, calculate the elevation difference and horizontal distance between the central sampling point and the adjacent sampling points of the observation unit, and generate the slope aspect change value of the central sampling point;

[0083] A slope aspect change sequence is generated by sliding the observation unit along a spiral sampling path;

[0084] The difference in slope aspect change values ​​between adjacent observation units is calculated based on the slope aspect change sequence, and the central sampling point whose slope aspect change value difference exceeds the preset difference range is marked as a terrain sampling point;

[0085] Calculate the path distance between adjacent terrain sampling points, determine the sampling point spacing based on the path distance and the slope aspect change value of the terrain sampling points, regenerate the sampling point positions according to the sampling point spacing, and output the target sampling data.

[0086] Obtain the location and elevation data of sampling points along the spiral sampling path. Extract the location coordinates and corresponding elevation values ​​of each sampling point along the constructed spiral sampling path. For each sampling point, record its two-dimensional coordinates on the horizontal plane and its corresponding terrain elevation value. Sampling points are typically distributed along the spiral path at regular angular intervals, such as one sampling point every 10 degrees. Elevation data can be extracted using a digital elevation model, with an accuracy typically of 5 meters. If a spiral path starts at a point of abrupt change in slope, with coordinates (500, 600) and elevation of 850 meters, a sampling point is set every 10 degrees along the spiral path, centered on this point. The first sampling point in the 0-degree direction is located at (650, 600) with an elevation of 820 meters.

[0087] An observation unit is constructed, consisting of three adjacent sampling points, to calculate the aspect change value of the central sampling point. Each observation unit contains three consecutive sampling points, labeled as front, middle, and rear. The elevation difference between the central sampling point and its two adjacent sampling points, as well as the horizontal distance between them, are calculated. The elevation difference is obtained by subtracting the elevation of the adjacent points from the elevation of the central point, and the horizontal distance is obtained by calculating the Euclidean distance between the two points on the horizontal plane. Based on the elevation difference and the horizontal distance, the aspect angles in the front-middle and middle-rear directions are calculated. The aspect angle represents the direction of terrain tilt, typically with true north as 0 degrees and increasing clockwise. The aspect change value is defined as the absolute value of the difference between the aspect angle in the front-middle direction and the aspect angle in the middle-rear direction. If the aspect angle in the front-middle direction is 30 degrees and the aspect angle in the middle-rear direction is 85 degrees, then the aspect change value of the central sampling point is 55 degrees.

[0088] A complete aspect change sequence is generated by sliding observation units along a spiral sampling path. Starting from the beginning of the spiral path, observation units are constructed sequentially according to the path order. Each time, one sampling point is moved forward, updating the three points in the observation unit and recalculating the aspect change value at the center point. This sliding window method obtains the aspect change values ​​for all sampling points except the first and last two, forming the aspect change sequence. For a spiral path containing 36 sampling points, 34 aspect change values ​​can be generated, constituting the aspect change sequence.

[0089] Based on the aspect change sequence, the difference in aspect change values ​​between adjacent observation units is calculated. For any two adjacent values ​​in the aspect change sequence, the absolute value of their difference is calculated as the aspect change value difference. Center sampling points where the aspect change value difference exceeds a preset range are marked as terrain sampling points. The preset range can be set to 15 degrees; that is, if the aspect change value difference between adjacent observation units is greater than 15 degrees, it is considered that there is a significant terrain change at that location, requiring focused sampling. In a set of aspect change sequences, if two adjacent aspect change values ​​are 25 degrees and 45 degrees respectively, with a difference of 20 degrees, exceeding the preset range, the corresponding center sampling point is marked as a terrain sampling point.

[0090] The path distance between adjacent terrain sampling points is calculated, i.e., the cumulative distance between two terrain sampling points is calculated along a spiral path. The path distance is obtained by summing the straight-line distances between adjacent sampling points. The sampling point spacing is determined based on the path distance and the slope aspect change value of the terrain sampling points. An adaptive strategy is adopted to determine the sampling point spacing: the larger the slope aspect change value, the more drastic the terrain change, and the smaller the sampling point spacing is set; conversely, the smaller the slope aspect change value, the larger the sampling point spacing is. The specific calculation method is: base spacing × inverse proportional factor of slope aspect change value. The base spacing can be set to 50 meters, and the inverse proportional factor can be defined as 100 ÷ slope aspect change value. If the slope aspect change value of a certain terrain sampling point is 50 degrees, then the corresponding sampling point spacing is 100 meters.

[0091] Based on the determined sampling point spacing, sampling point positions are regenerated between adjacent terrain sampling points. Starting from one terrain sampling point, new sampling points are evenly distributed along a spiral path according to the calculated sampling point spacing until the next terrain sampling point is reached. For each newly generated sampling point, its position coordinates and corresponding elevation value are recorded to form the final target sampling data. If the path distance between two adjacent terrain sampling points is 300 meters and the sampling point spacing is 100 meters, then two new sampling points are evenly distributed between them. Including the terrain sampling points at both ends, a total of four sampling points are used for wind field simulation of this path segment.

[0092] The final output target sampling data includes information such as the location coordinates, elevation values, and slope aspect values ​​of the sampling points. This information will be used in subsequent downscaling simulations of wind fields in mountainous areas. The output target sampling data is usually saved as a structured file for easy reading and processing by the wind field simulation model.

[0093] This invention optimizes the downscaling simulation of wind fields in mountainous areas through an adaptive sampling strategy, effectively improving simulation accuracy and computational efficiency. The aspect variation value, as an indicator of terrain complexity, accurately reflects the influence of mountainous terrain on the wind field, ensuring a high degree of match between the sampling point distribution and terrain features. The spiral sampling path design fully considers the directional characteristics of mountainous wind fields, capturing wind field changes in different directions. This method, by finely describing the complex terrain of mountainous areas, significantly improves the accuracy and reliability of wind field downscaling simulation, providing solid technical support for wind energy resource assessment and micrometeorological forecasting in mountainous areas.

[0094] Orthogonal decomposition is performed on the target sampling data to obtain the frequency components of the land surface undulation. The land surface fluctuation coefficient is calculated by combining the frequency components with meteorological monitoring data, including:

[0095] The target sampled data is decomposed into forward and reverse orthogonal components to obtain forward and reverse frequency components. The ratio of the forward and reverse frequency components is used to construct a frequency feature sequence.

[0096] Calculate the changing trend of the frequency feature sequence, divide the frequency feature sequence into multiple observation groups according to the inflection point of the changing trend, and calculate the degree of dispersion of the frequency features within each observation group.

[0097] Based on the degree of dispersion, the frequency clustering characteristics of the observation group are determined, and the curve of the frequency clustering characteristics changing with spatial location is used as the basis for partitioning to generate observation segments.

[0098] Meteorological monitoring data is acquired within the observation section, and a correspondence is established between the gradient changes and frequency clustering characteristics of the meteorological monitoring data to form a spatiotemporal feature sequence.

[0099] The fluctuation amplitude of frequency clustering features in the spatiotemporal feature sequence is calculated, and the weight ratio of frequency clustering features and meteorological gradient values ​​is dynamically adjusted according to the fluctuation amplitude to generate the surface fluctuation coefficient.

[0100] Orthogonal decomposition is performed on the target sampled data in both forward and reverse directions. Forward decomposition refers to the decomposition process from low frequency to high frequency, while reverse decomposition refers to the decomposition process from high frequency to low frequency. Orthogonal decomposition employs wavelet transform, using the elevation data of the target sampled points as input and selecting wavelet basis functions suitable for terrain feature analysis, such as the Daubechies wavelet basis. For elevation sequences arranged along a spiral path, they are treated as one-dimensional signals for decomposition. In the forward decomposition process, the low-frequency components, representing the macroscopic undulations of the terrain, are extracted first, followed by the progressive extraction of the high-frequency components, representing the microscopic details of the terrain. Five decomposition levels can be set in the forward decomposition to obtain frequency components at different scales. If a sampled point sequence contains 256 elevation values, the frequency components obtained after forward decomposition include five detail components and one approximate component. The detail components correspond to frequency features from high to low, and the approximate component corresponds to the lowest frequency feature.

[0101] The reverse decomposition process is similar to the forward decomposition, but the processing order is reversed. First, the high-frequency components are extracted, then the low-frequency components are extracted step by step. Five decomposition levels are set up to obtain the corresponding frequency components. The ratio of the forward frequency components to the reverse frequency components is used to construct a frequency feature sequence. The ratio is calculated by dividing the energy value of the forward frequency component at the corresponding level by the energy value of the reverse frequency component. The energy value is obtained by summing the squares of the frequency components. If the energy value of the third level in the forward decomposition is 0.025, and the energy value of the corresponding level in the reverse decomposition is 0.020, then the frequency feature value of that level is 1.25. The above calculation is repeated for all levels to form a complete frequency feature sequence.

[0102] The trend of the frequency characteristic sequence is calculated, and inflection points are identified. The trend is obtained by calculating the first difference of the frequency characteristic sequence, where the first difference value is the difference between two adjacent frequency characteristic values. An inflection point is defined as the position where the first difference value changes from positive to negative or vice versa. Based on the inflection points, the frequency characteristic sequence is divided into multiple observation groups. If there are three inflection points in the frequency characteristic sequence, located at positions 45, 120, and 210, the entire sequence is divided into four observation groups, containing frequency characteristic values ​​at positions 1-45, 46-120, 121-210, and 211-256, respectively. The dispersion of the frequency characteristic values ​​within each observation group is calculated, and the dispersion is represented by the standard deviation. A larger standard deviation indicates more significant fluctuations in the frequency characteristics within that observation group, and a greater degree of topographic relief.

[0103] Based on the calculated dispersion, the frequency clustering characteristics of each observation group are determined. These characteristics describe the distribution of topographic relief at different frequencies. A dispersion-weighted average method is used to calculate the frequency clustering characteristic values. The characteristic values ​​of each frequency level are multiplied by their corresponding weighting coefficients and then summed; the weighting coefficients are proportional to the dispersion. If an observation group contains five frequency levels with characteristic values ​​of 1.25, 0.98, 1.42, 0.85, and 1.10, and a dispersion of 0.23, while another observation group has a dispersion of 0.15, then the first observation group has a relatively larger weighting coefficient, and its calculated frequency clustering characteristic value better reflects the actual frequency distribution. The frequency clustering characteristics are constructed as a curve representing the change in spatial location, serving as the basis for dividing the observation area. Inflection points on the curve indicate significant changes in topographic features and can be used as boundaries for sections. Based on the inflection points of the curve, the entire observation area is divided into several continuous observation sections.

[0104] Meteorological monitoring data is acquired within the defined observation area. This data includes parameters such as wind speed, wind direction, temperature, and humidity, obtained through meteorological monitoring stations deployed within the observation area. The spatial gradient of the meteorological monitoring data is calculated. Gradient change represents the rate of change of meteorological parameters with spatial location; for example, the wind speed gradient represents the change in wind speed per unit distance. A correspondence is established between the gradient changes of the meteorological monitoring data and frequency clustering characteristics, forming a spatiotemporal feature sequence. This correspondence is established through correlation analysis, calculating the correlation coefficient between the meteorological gradient and the frequency clustering characteristics; a higher correlation coefficient indicates a stronger correlation. The spatiotemporal feature sequence contains the frequency clustering characteristic values ​​and meteorological gradient values ​​for the corresponding locations, as well as their correlation index.

[0105] The fluctuation amplitude of frequency clustering features in a spatiotemporal feature sequence is calculated. The fluctuation amplitude is obtained by calculating the standard deviation of the frequency clustering feature over time. A larger standard deviation indicates a more significant change in the frequency clustering feature over time. The weight ratio of the frequency clustering feature and the meteorological gradient value is dynamically adjusted based on the fluctuation amplitude. In areas with large fluctuation amplitudes, the weight of the frequency clustering feature is relatively increased; in areas with small fluctuation amplitudes, the weight of the meteorological gradient value is relatively increased. The adjusted weight ratio is used to calculate the surface fluctuation coefficient, which is calculated by multiplying the frequency clustering feature by its weight and then adding the meteorological gradient value multiplied by its weight. If the frequency clustering feature value at a certain spatiotemporal point is 1.35, the fluctuation amplitude is 0.25, and the corresponding meteorological gradient value is 0.18, and the weight ratio calculated based on the fluctuation amplitude is 0.65:0.35, then the surface fluctuation coefficient is 1.35×0.65+0.18×0.35=0.9405.

[0106] This invention effectively captures the influence mechanism of complex mountainous terrain on wind fields through orthogonal decomposition and frequency feature analysis, significantly improving the accuracy and applicability of wind field simulation. By adaptively weighting the frequency clustering features, it accurately reflects the changing patterns of wind fields under different terrain conditions, effectively solving the problem of insufficient simulation accuracy of traditional methods in rapidly changing mountainous terrain areas. This method can adapt to various complex terrain conditions, especially demonstrating superior performance in areas with dramatic terrain gradient changes such as canyons and ridges, providing more reliable technical support for mountainous wind energy resource assessment, atmospheric pollution diffusion simulation, and meteorological disaster early warning.

[0107] like Figure 2 As shown, Figure 2 This diagram illustrates the orthogonal decomposition frequency component analysis of this embodiment, showcasing the extraction effects of different technical solutions on terrain frequency features. This technical solution uses the Daubechies wavelet basis for orthogonal decomposition, achieving frequency component ratios of 1.25, 0.98, 1.42, 0.85, and 1.10 at the five decomposition levels. Compared to existing technologies (FFT spectral analysis) and traditional methods (linear regression), it demonstrates a higher frequency feature capture capability. Particularly at the third level, the ratio of this technical solution reaches 1.42, which is 49.5% higher than the 0.95 of FFT spectral analysis and 73.2% higher than the 0.82 of linear regression. This indicates that this technical solution has a significant advantage in extracting mid-frequency components, more accurately characterizing the mid-frequency features of terrain undulations, which is crucial for the refined description of complex mountainous terrain. At the first and fifth levels, this technical solution also achieves ratios of 1.25 and 1.10, respectively, demonstrating its excellent performance in extracting high-frequency and low-frequency features, laying a solid foundation for the subsequent accurate calculation of the surface fluctuation coefficient.

[0108] Using the surface wave coefficient to construct a boundary perturbation equation, solving the boundary perturbation equation yields the surface stress distribution, including:

[0109] An equilibrium state function is constructed based on the surface wave coefficient. The second derivative of the equilibrium state function is calculated. The second derivative is combined with the spatial gradient of the surface wave coefficient to form a boundary perturbation term.

[0110] Construct boundary perturbation equations that include dynamic and perturbation terms, where the dynamic term is determined by the time rate of change of the surface wave coefficient, and the perturbation term is determined by the spatial distribution of the boundary perturbation term;

[0111] The boundary perturbation equation is solved iteratively. In each iteration, the weight coefficients of the dynamic and perturbation terms are dynamically adjusted based on the residual value of the current solution until the residual value is less than the preset residual threshold to obtain a converged solution.

[0112] The surface stress tensor is calculated based on the converged solution, a stress propagation matrix is ​​constructed, the surface stress tensor is substituted into the stress propagation matrix, and the surface stress distribution is obtained through matrix recursion.

[0113] First, an equilibrium state function is constructed based on the surface wave coefficient. This function describes the stable state characteristics of the land surface under wind influence and is a function of the surface wave coefficient. A polynomial fitting method is used to construct the equilibrium state function, employing the surface wave coefficient as the independent variable to construct a fifth-order polynomial: In this equation, F(C) represents the equilibrium state function, C represents the surface fluctuation coefficient, a0 represents the constant term coefficient, a1 represents the linear term coefficient, a2 represents the quadratic term coefficient, a3 represents the cubic term coefficient, a4 represents the quartic term coefficient, and a5 represents the quintic term coefficient. These coefficients are obtained by fitting historical observation data using the least squares method, and their values ​​reflect the actual response characteristics of the surface in a specific region under different wind field conditions. The combination of coefficients determines the overall shape of the equilibrium state function curve, thus accurately describing the stability characteristics of the surface under various wind field intensities.

[0114] If the surface fluctuation coefficient is 1.25, substituting it into the fitted polynomial, the calculated equilibrium state function value may be: F(1.25)=0.875F(1.25)=0.875F(1.25)=0.875.

[0115] After the equilibrium state function is constructed, its second derivative is calculated, which is the second partial derivative of the equilibrium state function with respect to the surface coordinates. The second derivative reflects the spatial acceleration variation of surface fluctuations and is an important indicator for describing the intensity of surface disturbances. The central difference method is used to calculate the second derivative to improve computational accuracy. ,in, Let represent the second partial derivative of the equilibrium state function with respect to coordinate x, and let represent the curvature of the function in that direction. Let F(x) represent the equilibrium state function value at the point where the x-coordinate moves forward one step, and let F(x) represent the equilibrium state function value at the current point. This represents the equilibrium state function value at a point where the x-coordinate moves back one step. This represents the spatial distance between adjacent calculation points. If the equilibrium state function value of a point is 0.875, and the equilibrium state function values ​​of its two adjacent points in the horizontal direction are 0.920 and 0.825 respectively, and the distance between the adjacent points is 10 meters, then the second derivative of this point in the horizontal direction is 0.000175m. −2 .

[0116] The second derivative of the equilibrium state function is combined with the spatial gradient of the surface wave coefficient to form the boundary perturbation term. The spatial gradient of the surface wave coefficient represents the rate of change of surface wave intensity in space and is calculated using the finite difference method. The combination method employs a weighted summation: Where B represents the boundary disturbance term, and represents the intensity of the disturbance at the surface boundary. The weighting coefficients represent the second derivative. The second derivative of the equilibrium state function is represented. The weighting coefficients represent the spatial gradient. The spatial gradient represents the rate of change of the surface undulation coefficient in space. The weighting coefficients are dynamically adjusted based on terrain complexity; in areas with high terrain complexity, the weight of the second derivative is relatively increased; in areas with low terrain complexity, the weight of the spatial gradient is relatively increased. If the second derivative at a certain point is 0.000175m... -2 The spatial gradient is 0.012m. -1 With a terrain complexity score of 0.7 and weighting coefficients of 0.6 and 0.4, the calculated boundary disturbance term is 0.000585m. −2 .

[0117] A boundary perturbation equation is constructed, comprising a dynamic term and a disturbance term. The dynamic term is determined by the time rate of change of the surface wave coefficient, representing the time-varying characteristics of surface waves. The time rate of change is obtained by dividing the difference in the surface wave coefficient at consecutive time points by the time interval. ,in, The time rate of change of the surface wave coefficient represents the speed at which the wave coefficient changes over time. The surface fluctuation coefficient represents the time step forward. This represents the surface fluctuation coefficient at the current point in time. The time interval is 1 hour. If the surface fluctuation coefficient of a certain point at two adjacent time points is 1.25 and 1.32 respectively, then the time variation rate is 0.07.

[0118] The boundary perturbation equations are in the form of partial differential equations: ,in, The value represents the rate of change of surface stress over time, the stress state represents the speed at which it changes over time, and α represents the weighting coefficient of the dynamic term, indicating the proportion of the contribution of the rate of change over time to the stress change. β represents the time rate of change of the surface fluctuation coefficient, β represents the weighting coefficient of the disturbance term, β represents the contribution ratio of boundary disturbance to stress change, and B represents the boundary disturbance term.

[0119] The boundary perturbation equations are solved iteratively using the finite element method. The computational domain is discretized into a triangular mesh, and unknowns are defined at each mesh point. Based on the boundary perturbation equations and boundary conditions, a system of algebraic equations is established: [A]S = b, where [A] represents the coefficient matrix, S represents the unknown solution vector, and b represents the right-hand side term vector.

[0120] An improved conjugate gradient method is used to solve the system of equations. In each iteration, the weighting coefficients of the dynamic and perturbation terms are dynamically adjusted based on the residual value of the current solution. The residual value is defined as: r (k) =b-[A]S (k) Where b is the vector of the right-hand side of the system of algebraic equations, and r (k) Let S represent the residual vector of the k-th iteration, and let S represent the accuracy of the solution to the equation. (k) Let represent the solution vector for the k-th iteration.

[0121] The surface stress tensor is calculated based on the convergent solution obtained through iterative solving. The surface stress tensor describes the stress state experienced by the Earth's surface, including normal stress and shear stress components. To calculate the surface stress tensor, the convergent solution is first substituted into the stress-strain equation to obtain the strain components: ,in, Let f(S) represent the strain components, f(S) represent the degree of material deformation, f(S) represent the stress-strain transformation function, S represent the convergent solution obtained through iterative solving, and f(S) represent the stress state. Then, based on the constitutive relation of the material, the stress components are calculated. For linear elastic materials, stress is directly proportional to strain. ,in, The stress component represents the internal stress state of the material, and E represents the elastic modulus, indicating the stiffness characteristic of the material. This represents the strain component. If the strain component at a certain point is 0.002 and the elastic modulus of the surface material is 2000 MPa, then the corresponding stress component is 4 MPa.

[0122] After calculating each stress component, a stress tensor matrix is ​​constructed, with the diagonal elements of the matrix representing the normal stress components and the off-diagonal elements representing the shear stress components.

[0123] A stress propagation matrix is ​​constructed to describe the propagation law of stress on the Earth's surface. Based on stress wave propagation theory, the matrix considers the influence of terrain characteristics on stress propagation. Matrix elements represent the attenuation coefficient of stress propagating from the source point to the target point. The attenuation coefficient is related to the distance between the two points, the terrain relief, and material properties. The greater the distance, the more significant the attenuation; the greater the terrain relief, the faster the attenuation; and the greater the material damping, the faster the attenuation. If the distance from the source point to the target point is 100 meters, the terrain relief is 0.5, and the material damping is 0.2, the calculated attenuation coefficient may be 0.65. The surface stress tensor is substituted into the stress propagation matrix, and the stress value at each point is obtained through matrix multiplication. The above calculation is repeated for the entire computational domain to obtain the complete surface stress distribution.

[0124] The surface stress distribution is further optimized through matrix recursive operations. These operations consider the superposition effect of stress propagation, meaning that the stress state at a point is influenced not only by direct sources but also by stresses propagating from surrounding points. Starting from the boundary point, the recursive operation progresses inwards, calculating the stress state of one point at a time and using it as a new source point to influence surrounding points that have not yet been calculated. During the recursive process, stress values ​​are continuously updated until the stress values ​​at all points stabilize. The resulting surface stress distribution more accurately reflects the stress state of the surface under wind conditions.

[0125] This invention achieves a technological breakthrough in the refined simulation of wind fields under complex mountainous terrain conditions by constructing boundary perturbation equations and solving for surface stress distribution. The method fully considers the perturbation effect of surface fluctuations on the wind field and ensures high accuracy and stability of the calculation results through an iterative solution strategy that dynamically adjusts the weighting coefficients. The introduction of the stress propagation matrix effectively captures the spatial correlation of the wind field during its propagation on the surface, overcoming the limitations of traditional methods in stress propagation simulation.

[0126] The airflow motion state is calculated based on the surface stress distribution, and the airflow phase transition location is extracted. The vorticity equation is then constructed, including:

[0127] Calculate the stress gradient matrix based on surface stress distribution data;

[0128] A velocity potential function is constructed based on the stress gradient matrix, and the spatial derivative of the velocity potential function is calculated to obtain the airflow velocity vector field.

[0129] The curl component is calculated for the airflow velocity vector field, and the divergence component is calculated for the stress gradient matrix. The curl component and the divergence component are combined to construct the vorticity eigenvector.

[0130] Based on the vortex eigenvectors, a vortex equation coefficient matrix is ​​constructed, and the vortex equation coefficient matrix is ​​then divided into blocks and diagonalized to obtain the vortex equation.

[0131] The stress gradient matrix is ​​calculated based on surface stress distribution data. This matrix represents the rate of change of surface stress in space, reflecting the spatial non-uniformity of the surface stress state. The central difference method is used to calculate the stress gradient matrix, differentiating the surface stress distribution data in three spatial dimensions. In the horizontal direction, the stress values ​​of two adjacent grid points are selected, and the difference is calculated and divided by the distance between the two points to obtain the horizontal stress gradient. The same difference method is used to calculate the stress gradient in the vertical direction. For complex mountainous terrain, stress distribution typically exhibits significant spatial variation, resulting in larger stress gradient values. In ridge areas, the stress gradient is even more pronounced due to the abrupt changes in terrain. For a grid point in a mountainous area, the stress values ​​of two adjacent points in the horizontal direction are 4.2 MPa and 3.8 MPa, with a distance of 50 meters. Therefore, the stress gradient in this direction is 0.008 MPa / m. In the vertical direction, the stress values ​​of two adjacent points are 4.2 MPa and 3.5 MPa, with a height difference of 30 meters. Therefore, the stress gradient in the vertical direction is 0.023 MPa / m. The stress gradients in each direction are combined into a stress gradient matrix, and the matrix elements represent the rate of change of stress in different directions.

[0132] A velocity potential function is constructed based on the stress gradient matrix, describing the potential energy distribution of airflow. In constructing the velocity potential function, the stress gradient matrix is ​​converted into a velocity potential energy matrix, and the conversion process considers the influence mechanism of surface stress on airflow. Topographic parameters and airflow density parameters are introduced into the conversion formula to reflect the influence of topographic features and airflow characteristics on the velocity potential. For areas with larger stress gradients, the gradient of the velocity potential function is also larger, indicating that the airflow is driven by a stronger force. The topographic parameter is set to 0.65, and the airflow density parameter is set to 1.2 kg / m³. 3 For a point with a stress gradient of 0.008 MPa / m, the calculated velocity potential value is 0.00624 MPa. The spatial derivative of the velocity potential function is calculated to obtain the airflow velocity vector field. The spatial derivative represents the rate of change of the velocity potential function in space, pointing in the direction of the fastest decrease in the potential function, and its magnitude represents the numerical value of the rate of change. The spatial derivative is calculated using the finite difference method, differentiating in three spatial dimensions to obtain the three components of the velocity vector. For a point with a velocity potential value of 0.00624 MPa, its derivative in the horizontal direction is -0.00015 MPa / m, corresponding to a velocity component of 12.5 m / s; the derivative in the vertical direction is -0.00022 MPa / m, corresponding to a velocity component of 18.3 m / s. Combining the velocity components in the three directions constitutes the complete airflow velocity vector field.

[0133] The curl component is calculated for the airflow velocity vector field. Curl represents the rotational characteristics of airflow and is an important indicator for constructing vorticity characteristics. When calculating curl, the partial derivatives of each component of the velocity vector field are taken spatially, and a cross product is performed. The magnitude of the curl indicates the intensity of airflow rotation, and the direction indicates the orientation of the rotation axis. For complex mountainous terrain, the curl distribution often exhibits complex spatial variations due to airflow separation and convergence caused by topographic undulations. In ridge and canyon areas, the curl value is larger, indicating stronger airflow rotation. For points with velocity components of 12.5 m / s and 18.3 m / s, the calculated curl value is 0.35 / s. The divergence component is calculated for the stress gradient matrix. Divergence represents the divergence or convergence characteristics of airflow and reflects the local variation trend of airflow mass. When calculating divergence, the diagonal elements of the stress gradient matrix are summed to obtain the scalar divergence value. A positive divergence indicates airflow divergence, and a negative divergence indicates airflow convergence. For points with diagonal elements of the stress gradient matrix of 0.008, 0.005, and 0.023 MPa / m, the calculated divergence value is 0.036 MPa / m. The vorticity eigenvector is constructed by combining the curl and divergence components. This vorticity eigenvector contains information on the rotational and divergent characteristics of airflow motion and is a comprehensive index describing the airflow motion state. Each element of the vorticity eigenvector corresponds to one of the three curl components and one of the divergence components. For a point with a curl value of 0.35 / s and a divergence value of 0.036 MPa / m, the constructed vorticity eigenvector contains four elements, representing the curl components in three spatial directions and the divergence value, respectively.

[0134] A coefficient matrix for the vorticity equation is constructed based on the vorticity eigenvectors. The vorticity equation describes the variation of airflow vorticity with time and space and is the core equation for simulating airflow motion. When constructing the coefficient matrix, the elements of the vorticity eigenvectors are organized according to the form of the vorticity equation. The elements of the coefficient matrix represent the coefficients of each term in the vorticity equation, reflecting the physical characteristics of airflow motion. For cases where the vorticity eigenvectors contain four elements, the constructed coefficient matrix is ​​a 4×4 square matrix. The diagonal elements of the matrix represent the self-evolution term of vorticity, and the off-diagonal elements represent the interaction terms between the vorticity components. The values ​​of the diagonal elements are proportional to the corresponding elements of the vorticity eigenvectors, with the proportionality coefficient determined according to the airflow characteristics. The values ​​of the off-diagonal elements are related to the interaction strength of the vorticity components; in complex terrain, the interaction is usually more significant. For points with vorticity eigenvector elements of 0.25 / s, 0.15 / s, 0.2 / s, and 0.036 MPa / m, the diagonal elements of the constructed coefficient matrix are 0.15, 0.09, 0.12, and 0.0216, respectively, while the off-diagonal elements are determined according to the interaction strength, such as 0.03 and 0.02.

[0135] The coefficient matrix of the vorticity equation is block-diagonalized to transform it into a more easily solvable form. The purpose of block-diagonalization is to convert the coefficient matrix into a diagonal block matrix, with each diagonal block corresponding to a similar type of vorticity characteristic, thus simplifying the solution process. Block-diagonalization employs a similarity transformation method, constructing an appropriate transformation matrix to transform the original coefficient matrix into a block-diagonal form. The construction of the transformation matrix is ​​based on the eigenvalue and eigenvector analysis of the coefficient matrix, dividing elements with similar eigenvalues ​​into blocks. For a 4×4 coefficient matrix, possible block divisions include two 2×2 diagonal blocks or four 1×1 diagonal blocks, depending on the distribution of the matrix's eigenvalues. For the above coefficient matrix, assuming the eigenvalues ​​are 0.16, 0.14, 0.1, and 0.02, where the first two eigenvalues ​​are similar and the last two are similar, the matrix can be divided into two 2×2 diagonal blocks. The block-diagonalized matrix form is simpler and easier to solve. Based on the block-diagonalized coefficient matrix, the complete vorticity equation is constructed. The vorticity equation includes time derivative terms, convection terms, and diffusion terms, which respectively describe the change of vorticity over time, the transport of vorticity by the airflow, and the diffusion process of vorticity. The coefficients of the equation are determined by the coefficient matrix after block diagonalization, and the solution of the equation reveals the distribution law of airflow vorticity with time and space.

[0136] This invention achieves a precise description of airflow motion under complex terrain by constructing a vorticity equation. It organically combines surface stress distribution with airflow dynamics characteristics, comprehensively capturing the motion patterns of airflow in mountainous terrain through the calculation of stress gradient matrix, velocity potential function, curl component, and divergence component. The construction of vorticity eigenvectors and the block diagonalization of the vorticity equation coefficient matrix effectively improve computational efficiency and simulation accuracy, showing significant advantages, especially in complex terrain areas such as ridges and canyons.

[0137] Solving the vorticity equation yields the vorticity intensity and range of influence in the leeward region. Based on this vorticity intensity and range, the airflow motion correction parameters are calculated, including:

[0138] Solve the vorticity equation to obtain the vorticity intensity distribution in the leeward region, and determine the boundary of the leeward region's effective range based on the gradient threshold of the vorticity intensity distribution.

[0139] Within the leeward region's effective range boundary, spatial discrete points are established, and the vorticity intensity correction coefficient at these spatial discrete points is calculated. The vorticity intensity correction coefficient is then combined with the vorticity intensity field to generate airflow motion correction parameters.

[0140] Solving the vorticity equation yields the vorticity intensity distribution in the leeward region. The vorticity equation, a partial differential equation describing the variation of airflow vorticity, involves the temporal evolution, spatial distribution, and interaction of vorticity with the airflow velocity field. The finite difference method is used to solve the vorticity equation, discretizing the continuous physical space into grid points and solving for the vorticity value at each grid point. The computational domain covers the entire mountainous terrain, and a non-uniform grid strategy is adopted, with higher grid density in areas of dramatic terrain changes and lower grid density in areas of gentle terrain. The leeward region refers to the side of the mountain facing away from the prevailing wind direction, where airflow separation occurs as it bypasses the mountain, forming a complex vortex structure. The vorticity intensity in the leeward region is typically higher, reflecting the intensity of airflow rotation. Appropriate initial and boundary conditions are set during the solution process. The initial conditions are based on vorticity distributions provided by meteorological observation data or larger-scale wind field simulations; the boundary conditions specify vorticity values ​​at the inflow boundary and use extrapolation conditions at the outflow boundary. For a mountain with a height of 800 meters and a width of 2000 meters, when the prevailing wind direction is westerly and the wind speed is 10 m / s, the vorticity intensity distribution in the leeward region is obtained by solving the vorticity equation. The calculation results show that within 500 meters of the eastern side of the mountain, the maximum vorticity intensity reaches 0.15 m / s, indicating that the airflow rotation in this area is relatively intense, forming a distinct vortex structure.

[0141] The leeward region's effective range boundary is determined based on the gradient threshold of vorticity intensity distribution. The vorticity intensity gradient represents the rate of change of vorticity in space; regions with larger gradients are typically the interfaces between different vortex structures. To determine the leeward region's effective range, the spatial gradient of vorticity intensity is first calculated, i.e., the derivatives of vorticity intensity in three spatial directions. The gradient is calculated using the central difference method. For the vorticity intensity value at a grid point, the difference between its adjacent points is calculated and divided by the distance to obtain the gradient component at that point. The gradient components in the three directions are combined to obtain the total gradient value. The vorticity intensity gradient represents the change in vorticity intensity per unit distance. A gradient threshold of 0.0002 is set; regions greater than this threshold are identified as the boundary of the leeward region's effective range. By connecting all points that satisfy the gradient threshold condition, a closed boundary curve is formed; the area enclosed by this curve is the leeward region's effective range. For the aforementioned mountain, the longitudinal depth of the leeward zone is approximately twice the mountain's height, or 1600 meters; the lateral width is approximately 1.5 times the mountain's width, or 3000 meters; and the vertical extension reaches 1.2 times the summit height, or 960 meters. The determination of this effective range takes into account aerodynamic characteristics and topographical influences, accurately reflecting the spatial distribution of the leeward vortex structure.

[0142] Spatially discrete points are established within the boundary of the leeward region's influence area to facilitate subsequent calculation of the vorticity intensity correction coefficient. The distribution density of these discrete points is related to the rate of change of vorticity intensity; the point density is higher in areas with rapid vorticity changes and lower in areas with gradual changes. An adaptive meshing technique is used to generate the discrete points, dynamically adjusting the point density distribution based on the vorticity intensity gradient value. Within the leeward region's influence area boundary, a basic mesh resolution of 50 meters is set, and the mesh resolution is locally refined based on the vorticity intensity gradient value; the larger the gradient value, the higher the mesh resolution. For areas with vorticity intensity gradient values ​​exceeding 0.0004, the mesh resolution is increased to 25 meters; for areas with gradient values ​​below 0.0001, the mesh resolution is reduced to 100 meters. This adaptive meshing strategy ensures a detailed description of key areas even with limited computational resources.

[0143] The correction coefficient for vorticity intensity at spatially discrete points is calculated. This correction coefficient reflects the degree of influence of vorticity on airflow motion. The calculation of the correction coefficient considers multiple factors, including vorticity intensity value, vorticity intensity gradient, the relative position of the discrete point to the mountain, and local airflow stability. Weighting factors are introduced into the calculation formula, with the weights of different factors determined according to their influence on airflow motion. The weight of the vorticity intensity value is 0.4, indicating that vorticity intensity itself is the main factor affecting airflow motion; the weight of the vorticity intensity gradient is 0.3, reflecting the influence of the non-uniformity of vorticity distribution on airflow motion; the weight of the relative position of the discrete point to the mountain is 0.2, considering the constraint effect of terrain on airflow motion; and the weight of local airflow stability is 0.1, representing the moderating effect of airflow stability on the influence of vorticity. For a discrete point with a vorticity intensity value of 0.12 / s, a vorticity intensity gradient of 0.0003, a distance of 500 meters from the mountain, and a local airflow stability index of 0.8, the calculated correction coefficient is 0.94. The correction factor typically ranges from 0.5 to 1.5; a larger value indicates a stronger influence of vorticity on airflow.

[0144] The vorticity intensity correction coefficient is combined with the vorticity intensity field to generate airflow motion correction parameters. These parameters are comprehensive indices describing the airflow characteristics in the leeward region, used to correct large-scale wind field simulation results and achieve refined downscaling of the wind field. The combination method uses a weighted product: the product of the correction coefficient and the vorticity intensity, multiplied by a spatial decay function. The spatial decay function describes the decay of the correction parameters with increasing distance from the boundary, typically using an exponential decay form. The decay constant is determined based on the spatial scale of the leeward region; the larger the scale, the slower the decay. For the aforementioned leeward region, the decay constant is set to 0.002 / m. For a discrete point with a correction coefficient of 0.94, a vorticity intensity of 0.12 / s, and a distance of 800 meters from the boundary, the calculated airflow motion correction parameter is 0.094 / s. The distribution of the airflow motion correction parameters reflects the influence of the vortex structure in the leeward region on airflow motion, providing an important physical basis for subsequent wind field downscaling. By applying the correction parameters to the large-scale wind field simulation results, a refined simulation of the wind field under complex mountainous terrain was achieved, improving the simulation accuracy, especially the ability to characterize the complex vortex structure in the leeward region.

[0145] This invention achieves an accurate characterization of airflow characteristics in leeward regions under complex mountainous terrain by solving the vorticity equation and calculating airflow motion correction parameters. The introduction of airflow motion correction parameters enables wind field simulation results to accurately reflect the disturbance effect of terrain on airflow, especially providing a more precise description of the formation, development, and dissipation processes of vortex structures in leeward regions. This provides high-precision wind field data support for mountainous wind energy resource assessment, atmospheric pollutant diffusion prediction, and extreme weather event early warning.

[0146] like Figure 3 As shown, by performing downscaling iterative calculations on the airflow motion correction parameters and meteorological monitoring data, target area wind field data with a spatial resolution of hundreds of meters is generated, including:

[0147] The calculation region is divided into multiple layers based on the spatial resolution of the meteorological monitoring data. The calculation region is further subdivided into calculation units layer by layer. Interpolation is performed on the meteorological monitoring data to obtain the initial wind field data of the calculation unit.

[0148] Obtain the terrain height data of the calculation unit, calculate the slope change value of the terrain height data, and weight the terrain height data according to the slope change value to obtain the terrain enhancement factor;

[0149] The wind field correction is generated by multiplying the airflow motion correction parameter with the terrain enhancement factor, and the wind field correction is superimposed with the initial wind field data to obtain the corrected wind field data.

[0150] Based on the principle of divergence conservation, a conservation equation is constructed. The corrected wind field data is substituted into the mapped wind field data at the calculation node of the conservation equation, and the difference between the corrected wind field data and the mapped wind field data is calculated to obtain the wind field residual.

[0151] The wind field residuals are compared with the convergence threshold. For the calculation units that are greater than the convergence threshold, the steps of generating corrected wind field data and calculating wind field residuals are repeated until the wind field residuals of all calculation units are less than the convergence threshold.

[0152] Based on the converged wind field data, target area wind field data with a spatial resolution of hundreds of meters is generated.

[0153] The computational domain is divided into multiple layers based on the spatial resolution of meteorological monitoring data. Each computational domain is further subdivided into computational units, and interpolation is performed on the meteorological monitoring data to obtain the initial wind field data for each computational unit. Meteorological monitoring data typically has low spatial resolution, which cannot directly meet the requirements for refined simulation of wind fields in mountainous areas. The division of the multi-layered computational domain uses a nested grid method, gradually transitioning from a coarser outer grid to a finer inner grid. The outer grid covers the entire simulation area and its surroundings, with a resolution matching the meteorological monitoring data. The resolution of the inner grid gradually increases. For typical mountainous areas, the outer grid resolution is set to 5 kilometers, the middle layer to 1 kilometer, and the inner layer to 200 meters. The computational unit is the basic computational unit for wind field downscaling. Based on the target resolution requirements, the inner grid is further subdivided into computational units. Each computational unit is rectangular, with a side length on the order of hundreds of meters, the same as the target spatial resolution. Interpolation is performed on the meteorological monitoring data to obtain the initial wind field data for each computational unit. The interpolation method uses distance-weighted interpolation, where a weight coefficient is calculated based on the distance from the computational unit to the meteorological monitoring point, and the weight is proportional to the reciprocal of the distance. For a certain calculation unit, the wind speeds at three surrounding monitoring points are 5.2 m / s, 6.1 m / s, and 4.8 m / s, respectively, and the distances are 8.5 km, 6.2 km, and 9.1 km, respectively. The initial wind speed of the unit is calculated to be 5.5 m / s.

[0154] The process involves acquiring terrain elevation data for a computational unit, calculating the slope variation value of this data, and then weighting the elevation data based on the slope variation value to obtain a terrain enhancement factor. Terrain elevation data reflects the undulation of the land surface and is a key factor influencing wind field distribution in mountainous areas. The elevation data is acquired using a digital elevation model with a resolution of 30 meters to ensure the capture of subtle changes in the mountainous terrain. The slope variation value is calculated; slope represents the steepness of the terrain and is calculated as the ratio of the elevation difference between two adjacent points to the horizontal distance. The slope variation value indicates the rate of change of slope in space, reflecting the complexity of the terrain. For a specific computational unit, the surrounding terrain elevations are 520 meters, 580 meters, 540 meters, and 510 meters, respectively. The calculated average slope is 0.15, and the slope variation value is 0.002 / meter. The terrain elevation data is then weighted based on the slope variation value to obtain the terrain enhancement factor. The weighting method uses an exponential function. When the slope change is small, the terrain enhancement factor is close to 1; when the slope change is large, the terrain enhancement factor increases, enhancing the influence of terrain on the wind field. For a calculation unit with a slope change of 0.002 / m, the terrain enhancement factor is calculated to be 1.35.

[0155] The wind field correction is generated by multiplying the airflow correction parameter by the terrain enhancement factor. This correction is then superimposed on the initial wind field data to obtain the corrected wind field data. The airflow correction parameter reflects the influence of airflow rotation characteristics and vortex structure on the wind field, while the terrain enhancement factor represents the amplifying effect of terrain complexity on the wind field. Multiplying both comprehensively considers airflow dynamics and terrain influence, providing a more accurate description of wind field distribution in mountainous areas. For a calculation unit with an airflow correction parameter of 0.094 m / s and a terrain enhancement factor of 1.35, the wind speed correction component is calculated to be 1.2 m / s, and the wind direction correction component is calculated to be 15 degrees. The corrected wind field data is obtained by superimposing the wind field correction on the initial wind field data. Wind speed superposition uses vector addition, projecting the correction component onto the initial wind direction for vector superposition; wind direction superposition directly adds the correction angle to the initial wind direction. For a calculation unit with an initial wind speed of 5.5 m / s and a wind direction of 10 degrees east of north, the superimposed corrected wind speed is 6.4 m / s, and the corrected wind direction is 25 degrees east of north.

[0156] Based on the principle of divergence conservation, a conservation equation is constructed. The corrected wind field data is substituted into the mapped wind field data at the calculation nodes of the conservation equation, and the difference between the corrected and mapped wind field data is calculated to obtain the wind field residual. The principle of divergence conservation is a fundamental principle of fluid mechanics, stating that without considering mass source terms, the divergence of a fluid should be zero, i.e., the mass of airflow flowing into a region should be equal to the mass of airflow flowing out of that region. For each calculation unit, the airflow flux on its six faces is calculated; the sum of these fluxes should be zero. The corrected wind field data is substituted into the conservation equation to calculate the divergence value for each calculation unit. Based on the divergence value, the mapped wind field data satisfying divergence conservation is derived. For a calculation unit with a corrected wind speed of 6.4 m / s and a wind direction of 25 degrees east of north, the calculated divergence value is 0.0004 / s, and the mapped wind speed is 6.2 m / s with a wind direction of 22 degrees east of north. The difference between the corrected and mapped wind field data is calculated to obtain the wind field residual. For the above calculation unit, the wind speed residual is 0.2 m / s and the wind direction residual is 3 degrees.

[0157] The wind field residuals are compared with a convergence threshold. For computational units with residuals greater than the convergence threshold, the steps of generating corrected wind field data and calculating wind field residuals are repeated until the wind field residuals of all computational units are less than the convergence threshold. The convergence threshold is set at 0.1 m / s for wind speed residuals and 2 degrees for wind direction residuals. When the wind field residuals of a computational unit are less than these thresholds, the wind field calculation of that unit is considered to have converged. For computational units with wind field residuals greater than the convergence threshold, iterative calculation is required. During the iteration process, the mapped wind field data is used as the initial wind field data for the new round. The wind field correction and corrected wind field data are recalculated, and mapping and residual calculation are performed again. For computational units with wind speed residuals of 0.2 m / s and wind direction residuals of 3 degrees, iterative calculation is performed. After the second iteration, the wind speed residual is 0.08 m / s and the wind direction residual is 1.5 degrees, which meets the convergence condition, and the iteration terminates.

[0158] Based on the converged wind field data, target area wind field data with a spatial resolution of hundreds of meters is generated. The converged wind field data has already considered aerodynamic characteristics, topographic effects, and conservation constraints, and can accurately describe the wind field distribution under complex mountainous terrain. When generating the target area wind field data, the wind field data of each calculation unit is organized according to spatial location to form regular grid data. The data format includes three basic elements: spatial coordinates, wind speed, and wind direction angle, with a resolution of 100 meters, covering the entire target area. The data is represented in vector form, and the wind speed is decomposed into east-west and north-south components for easy subsequent analysis and application. The generated high-resolution wind field data can clearly show the local characteristics of the mountain wind field, such as the ridge acceleration effect, valley guiding effect, and leeward vortex, providing a high-precision data foundation for mountainous wind energy resource assessment and atmospheric pollution diffusion simulation.

[0159] This invention achieves refined simulation of wind fields in complex mountainous terrains by performing downscaling iterative calculations on airflow motion correction parameters and meteorological monitoring data. This method comprehensively considers airflow dynamics and terrain influences, overcoming the limitations of traditional wind field simulations in terms of spatial resolution and physical mechanism description through techniques such as multi-layer nested mesh generation, terrain enhancement factor calculation, wind field correction generation, and divergence conservation constraints.

[0160] One technical solution provided in this embodiment of the invention is an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the steps of the method described in any of the foregoing embodiments.

[0161] One technical solution provided in this embodiment of the invention is a computer-readable storage medium storing computer program instructions, which, when executed by a processor, implement the steps in the method described in any of the preceding claims.

[0162] The specific embodiments described above are preferred embodiments of the present invention and are not intended to limit the specific scope of the present invention. The scope of the present invention includes, but is not limited to, these specific embodiments. All equivalent changes made in accordance with the shape and structure of the present invention are within the protection scope of the present invention.

Claims

1. A method for downscaling and optimizing wind fields in mountainous areas based on high-precision simulation, characterized in that, Includes the following steps: Obtain topographic elevation data and meteorological monitoring data for the target mountainous area; Identify abrupt slope changes in terrain elevation data and construct a spiral sampling path using these abrupt slope changes as starting points. Calculate the slope aspect variation value along the spiral sampling path, determine the sampling point spacing based on the slope aspect variation value, and generate target sampling data; Orthogonal decomposition is performed on the target sampling data to obtain the frequency components of the surface undulation. The surface fluctuation coefficient is calculated by combining the frequency components with meteorological monitoring data. The boundary perturbation equation is constructed using the surface wave coefficient, and the surface stress distribution is obtained by solving the boundary perturbation equation. The airflow motion state is calculated based on the surface stress distribution, and the airflow phase transition position is extracted. A vorticity equation is constructed, and the vorticity equation is solved to obtain the vorticity intensity and range of action in the leeward region. Based on the vorticity intensity and range of action in the leeward region, the airflow motion correction parameters are calculated. By performing downscaling iterative calculations on airflow motion correction parameters and meteorological monitoring data, target area wind field data with a spatial resolution of hundreds of meters is generated. Identifying abrupt slope changes in terrain elevation data and constructing a spiral sampling path using these abrupt changes as starting points includes: The terrain elevation data is divided into grid cells. The elevation difference between each grid point in the grid cell and its adjacent grid points in the horizontal and vertical directions is calculated. The slope value of the grid point is generated based on the elevation difference. Grid points whose slope values ​​fall within a preset slope threshold range are statistically analyzed and marked as candidate points for slope abrupt changes. Calculate the spacing between candidate points of slope abrupt change, divide the candidate points of slope abrupt change into multiple point sets based on the spacing, and extract the candidate point with the largest slope value from each point set as the slope abrupt change point; Calculate the elevation change value along the orthogonal coordinate axis and diagonal direction at the slope abrupt change point, and determine the direction with the largest elevation change value as the sampling direction; The starting angle of the spiral path is set based on the sampling direction. The spiral sampling path is constructed with the elevation change value as the center, using the slope change point as the center.

2. The method according to claim 1, characterized in that, Calculate the slope aspect variation value along the spiral sampling path, determine the sampling point spacing based on the slope aspect variation value, and generate target sampling data including: Obtain the location and elevation data of sampling points along the spiral sampling path; Construct an observation unit consisting of three adjacent sampling points, calculate the elevation difference and horizontal distance between the central sampling point and the adjacent sampling points of the observation unit, and generate the slope aspect change value of the central sampling point; A slope aspect change sequence is generated by sliding the observation unit along a spiral sampling path; The difference in slope aspect change values ​​between adjacent observation units is calculated based on the slope aspect change sequence, and the central sampling point whose slope aspect change value difference exceeds the preset difference range is marked as a terrain sampling point; Calculate the path distance between adjacent terrain sampling points, determine the sampling point spacing based on the path distance and the slope aspect change value of the terrain sampling points, regenerate the sampling point positions according to the sampling point spacing, and output the target sampling data.

3. The method according to claim 1, characterized in that, Orthogonal decomposition is performed on the target sampling data to obtain the frequency components of the land surface undulation. The land surface fluctuation coefficient is calculated by combining the frequency components with meteorological monitoring data, including: The target sampled data is decomposed into forward and reverse orthogonal components to obtain forward and reverse frequency components. The ratio of the forward and reverse frequency components is used to construct a frequency feature sequence. Calculate the changing trend of the frequency feature sequence, divide the frequency feature sequence into multiple observation groups according to the inflection point of the changing trend, and calculate the degree of dispersion of the frequency features within each observation group. Based on the degree of dispersion, the frequency clustering characteristics of the observation group are determined, and the curve of the frequency clustering characteristics changing with spatial location is used as the basis for partitioning to generate observation segments. Meteorological monitoring data is acquired within the observation section, and a correspondence is established between the gradient changes and frequency clustering characteristics of the meteorological monitoring data to form a spatiotemporal feature sequence. The fluctuation amplitude of frequency clustering features in the spatiotemporal feature sequence is calculated, and the weight ratio of frequency clustering features and meteorological gradient values ​​is dynamically adjusted according to the fluctuation amplitude to generate the surface fluctuation coefficient.

4. The method according to claim 1, characterized in that, Using the surface wave coefficient to construct a boundary perturbation equation, solving the boundary perturbation equation yields the surface stress distribution, including: An equilibrium state function is constructed based on the surface wave coefficient. The second derivative of the equilibrium state function is calculated. The second derivative is combined with the spatial gradient of the surface wave coefficient to form a boundary perturbation term. Construct boundary perturbation equations that include dynamic and perturbation terms, where the dynamic term is determined by the time rate of change of the surface wave coefficient, and the perturbation term is determined by the spatial distribution of the boundary perturbation term; The boundary perturbation equation is solved iteratively. In each iteration, the weight coefficients of the dynamic and perturbation terms are dynamically adjusted based on the residual value of the current solution until the residual value is less than the preset residual threshold to obtain a converged solution. The surface stress tensor is calculated based on the converged solution, a stress propagation matrix is ​​constructed, the surface stress tensor is substituted into the stress propagation matrix, and the surface stress distribution is obtained through matrix recursion.

5. The method according to claim 1, characterized in that, The airflow motion state is calculated based on the surface stress distribution, and the airflow phase transition location is extracted. The vorticity equation is then constructed, including: Calculate the stress gradient matrix based on surface stress distribution data; A velocity potential function is constructed based on the stress gradient matrix, and the spatial derivative of the velocity potential function is calculated to obtain the airflow velocity vector field. The curl component is calculated for the airflow velocity vector field, and the divergence component is calculated for the stress gradient matrix. The curl component and the divergence component are combined to construct the vorticity eigenvector. Based on the vortex eigenvectors, a vortex equation coefficient matrix is ​​constructed, and the vortex equation coefficient matrix is ​​then divided into blocks and diagonalized to obtain the vortex equation.

6. The method according to claim 1, characterized in that, Solving the vorticity equation yields the vorticity intensity and range of influence in the leeward region. Based on this vorticity intensity and range, the airflow motion correction parameters are calculated, including: Solve the vorticity equation to obtain the vorticity intensity distribution in the leeward region, and determine the boundary of the leeward region's effective range based on the gradient threshold of the vorticity intensity distribution. Within the leeward region's effective range boundary, spatial discrete points are established, and the vorticity intensity correction coefficient at these spatial discrete points is calculated. The vorticity intensity correction coefficient is then combined with the vorticity intensity field to generate airflow motion correction parameters.

7. The method according to claim 1, characterized in that, By performing downscaling iterative calculations on airflow motion correction parameters and meteorological monitoring data, target area wind field data with a spatial resolution of hundreds of meters is generated, including: The calculation region is divided into multiple layers based on the spatial resolution of the meteorological monitoring data. The calculation region is further subdivided into calculation units layer by layer. Interpolation is performed on the meteorological monitoring data to obtain the initial wind field data of the calculation unit. Obtain the terrain height data of the calculation unit, calculate the slope change value of the terrain height data, and weight the terrain height data according to the slope change value to obtain the terrain enhancement factor; The wind field correction is generated by multiplying the airflow motion correction parameter with the terrain enhancement factor, and the wind field correction is superimposed with the initial wind field data to obtain the corrected wind field data. Based on the principle of divergence conservation, a conservation equation is constructed. The corrected wind field data is substituted into the mapped wind field data at the calculation node of the conservation equation, and the difference between the corrected wind field data and the mapped wind field data is calculated to obtain the wind field residual. The wind field residuals are compared with the convergence threshold. For the calculation units that are greater than the convergence threshold, the steps of generating corrected wind field data and calculating wind field residuals are repeated until the wind field residuals of all calculation units are less than the convergence threshold. Based on the converged wind field data, target area wind field data with a spatial resolution of hundreds of meters is generated.

8. An electronic device, characterized in that, include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the steps of the method as described in any one of claims 1 to 7.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer program instructions that, when executed by a processor, implement the steps of the method as described in any one of claims 1 to 7.

Images