Pollen concentration prediction method and system based on big data analysis

Abstract

The invention relates to the technical field of learning methods, in particular to a pollen concentration prediction method and system based on big data analysis, and the method comprises the following steps: inputting a diffusion model through an elevation mutation rate matrix to calculate path nodes, screening and optimizing a path node combination, fitting a trend line, extracting a main path set, and identifying a risk unit; and executing local prediction path updating, and outputting a corrected concentration prediction and early warning signal. According to the method, a path dynamic model is established by converting terrain elevation change into mutation rate parameters, path nodes are screened in combination with landform blocking frequency factors generated by elevation difference and a wind speed direction difference angle, a trend line is reconstructed by adopting a cubic spline interpolation method, and a main path is dynamically corrected based on the direction difference angle. A slope normal angle and a wind pressure linkage correction coefficient are fused to identify an abnormal area, the modeling precision of a model is improved by integrating a terrain obstruction frequency, a slope angle difference and a dynamic wind pressure parameter, and the nonlinear effect of terrain obstruction and meteorological parameters is synchronously reflected.

Description

Technical Field

[0001] The present invention relates to the technical field of learning methods, and particularly to a pollen concentration prediction method and system based on big data analysis. Background Art

[0002] The technical field of learning methods includes computer systems and methods based on specific computational models. The core of this field lies in training models through sample data to enable them to have the ability of automatic reasoning and judgment and apply them to specific tasks. This field covers model construction and training methods such as supervised learning, unsupervised learning, and reinforcement learning, relying on algorithm structures such as neural networks, support vector machines, and decision trees, and completing the model construction process through the mapping relationship between sample features and target outputs. Its application scope covers multiple directions such as natural language processing, image recognition, and predictive analysis, has a technical system that promotes model iteration and optimization through a large amount of data accumulation and feature extraction, and emphasizes the verification of model performance evaluation and generalization ability.

[0003] Among them, the pollen concentration prediction method refers to constructing a prediction model based on historical meteorological data and pollen concentration data to estimate the change trend of pollen concentration in a specific area in the future. This method mainly aims at the dynamic change of pollen concentration under the influence of different time and environmental factors. By collecting multi-dimensional meteorological parameters such as temperature, humidity, wind speed, and precipitation in historical data, combined with known pollen concentration values, a regression model in supervised learning is used to complete data training. After model training, by inputting real-time meteorological data into the trained regression model, the prediction output of pollen concentration in the future period is realized. In the model construction process, gradient boosting regression or support vector regression methods are generally used for training optimization, and the stability and prediction accuracy of the model are evaluated according to the cross-validation method.

[0004] Existing methods rely on the linear regression relationship between meteorological parameters and pollen concentration, ignoring the physical barrier effect of terrain elevation changes on the diffusion path. The model does not establish a quantitative correlation between terrain elevation difference and diffusion resistance, resulting in the deviation of the predicted path from the actual observation trajectory in the transition area between plains and mountains. Static regression processing does not consider the local wind direction disturbance caused by terrain undulation, and cannot dynamically correct the deviation when the actual diffusion path has an angle with the meteorological derivation path. The model training lacks the parameter integration of the slope angle and the main wind direction difference, and the pollen deposition phenomenon caused by terrain reflection in the steep slope area is not effectively identified. Risk determination is only based on a single dimension of the concentration change rate, without associating the dynamic coupling relationship between the wind pressure time series change and the terrain slope direction, and the spatial distribution characteristics of high-risk areas cannot be accurately captured under strong convective weather. Summary of the Invention

[0005] The object of the present invention is to solve the drawbacks existing in the prior art, and to propose a pollen concentration prediction method and system based on big data analysis.

[0006] To achieve the above object, the present invention adopts the following technical solutions: A pollen concentration prediction method based on big data analysis, comprising the following steps:

[0007] S1: Perform two-dimensional grid division through regional digital elevation model data, calculate the elevation difference between the center points of adjacent grid cells, generate the elevation mutation rate per unit distance based on the grid spacing, construct an elevation mutation rate matrix, and input the elevation mutation rate matrix into the pollen diffusion model;

[0008] S2: Call the pollen diffusion model to output a path node sequence, extract the elevation mutation rate per unit distance of adjacent node pairs, divide the sum by the total number of nodes to generate a topographic barrier frequency factor, and screen node pairs that meet the conditions that the elevation mutation rate is less than 50% of the original path segment and the wind direction difference angle is less than 30 degrees to generate an optimized path node combination;

[0009] S3: Input the optimized path node combination into the cubic spline interpolation method, fit to generate a diffusion path trend line, calculate the direction difference angle between the fitted path and the original path, and add the paths with a difference angle less than 15 degrees to the main path set;

[0010] S4: Calculate the slope aspect angle difference through the difference between the slope normal angle of the grid cell and the main wind direction angle, combine the wind pressure time series data to generate a wind pressure linkage correction coefficient, and mark the units in the main path set with a concentration change rate exceeding 40% and a positive increasing wind pressure linkage correction coefficient as risk units.

[0011] As a further solution of the present invention, the elevation mutation rate matrix specifically includes elevation difference, grid spacing, and elevation mutation rate per unit distance. The optimized path node combination includes the screened node pairs and the corrected concentration distribution parameters. The main path set specifically refers to the fitted path trend line, the direction difference angle determination result, and the spatio-temporal diffusion trajectory simulation parameters. The risk units include high concentration change rate units and positive increasing wind pressure linkage correction coefficient units.

[0012] As a further solution of the present invention, the acquisition steps of the pollen diffusion model are specifically as follows:

[0013] S101: Obtain regional digital elevation model data, set grid spacing parameters, divide the region into two-dimensional grid cells, extract the plane coordinates of the center points of each grid cell, and calculate the elevation values of the grid center points according to the spatial correspondence between the elevation point coordinates and the grid center points to generate a grid cell elevation data set;

[0014] S102: Invoke the elevation dataset of the grid cells, traverse the row and column indices of the grid cells, extract the adjacent grid cells in the four directions of up, down, left, and right, calculate the absolute value of the elevation difference between the current grid center point and the adjacent grid center points, divide the elevation difference in each direction by the grid spacing parameter to obtain the elevation mutation rate per unit distance;

[0015] S103: Invoke the elevation mutation rate per unit distance, arrange the data in the order of grid row and column indices, fill the mutation rate values of the missing boundary grids with zero, construct a determinant-structured matrix data, generate an elevation mutation rate matrix, and transfer it to the pollen diffusion model.

[0016] As a further solution of the present invention, the steps for obtaining the optimized path node combination are specifically as follows:

[0017] S201: Invoke the path node sequence output by the pollen diffusion model, extract the elevation difference and the straight-line distance of all adjacent node pairs, divide the absolute value of the elevation difference by the straight-line distance, calculate the elevation mutation rate per unit distance of each node pair, summarize all the mutation rate values, and generate an adjacent node pair elevation mutation rate set;

[0018] S202: Based on the adjacent node pair elevation mutation rate set, accumulate all the mutation rate values, divide the accumulated result by the total number of nodes in the node sequence, calculate the ratio of the total mutation rate to the number of nodes, and generate a topographic barrier frequency factor;

[0019] S203: Traverse the adjacent node pair elevation mutation rate set, compare the mutation rate of each node pair with 50% of the topographic barrier frequency factor, screen out the node pairs with mutation rates less than the threshold, synchronously extract the wind speed direction vectors of the screened node pairs, calculate the included angle between the vectors, eliminate the node pairs with an included angle exceeding 30 degrees, and reorganize the remaining node pairs in the original order to generate an optimized path node combination.

[0020] As a further solution of the present invention, the steps for obtaining the main path set are specifically as follows:

[0021] S301: Invoke the optimized path node combination, arrange the node coordinate sequence in the input order, use the cubic spline interpolation method to calculate the cubic polynomial coefficients between adjacent nodes, ensure that the function is continuous and second-order differentiable at the nodes, construct a piecewise cubic polynomial function, and generate a continuous and differentiable diffusion path trend line;

[0022] S302: Based on the diffusion path trend line, extract the original path tangent direction angle at the equally spaced sampling points of the path, calculate the difference between the fitted path tangent direction angles at the corresponding positions, introduce the path curvature difference, sampling point density, and path length scale factor, and use the formula:

[0023]

[0024] Calculate the curvature-distance weighted angular offset of multiple sampling points, accumulate it and then normalize it in combination with the path length factor to generate the direction difference angle;

[0025] Among them, θ new represents the direction difference angle, Δα i represents the difference between the tangent angles of the original path and the fitted path at the i-th sampling point, represents the Euclidean distance of the path segment between the i-th sampling point and the previous node, K i represents the absolute value of the curvature difference between the original path and the fitted path at the i-th sampling point, D i represents the unit length sampling density of the path segment where the i-th sampling point is located, R represents the ratio factor of the total length of the current path to the length of the reference path, and m is the total number of sampling points;

[0026] S303: According to the direction difference angle, set a preset angle threshold as the screening criterion, traverse all diffusion paths, extract the improved direction difference angle values of each path, compare the direction difference angle values with the preset angle threshold, screen out the paths that meet the conditions and store them in a set to generate the main path set.

[0027] As a further solution of the present invention, the acquisition step of the risk unit is specifically:

[0028] S401: Obtain the slope normal angle and the main wind direction angle of the grid unit, calculate the cosine value of the included angle between the slope normal vector and the main wind direction vector, use the result of the inverse cosine calculation as the angle difference, take the absolute value of the angle difference to eliminate the influence of direction symmetry, and compare the absolute value result with the preset wind direction deviation reference threshold to generate the slope aspect angle difference of the grid unit;

[0029] S402: Based on the slope aspect angle difference of the grid unit and the wind pressure time series data, use the formula:

[0030]

[0031] Calculate to obtain the dynamic correction factor, and combine it with the time window mean value of the wind pressure time series data to generate the wind pressure linkage correction coefficient;

[0032] Among them, C adj represents the wind pressure linkage correction coefficient, Φ represents the slope aspect angle difference of the grid unit, P t represents the wind pressure value at the t-th moment, represents the mean value of the wind pressure time series, n represents the total amount of time series data, H v represents the vertical height change rate of the unit, R s represents the ratio of the wind speeds of adjacent units, β represents the terrain curvature influence factor, S c represents the surface roughness coefficient, and τ represents the time decay factor;

[0033] S403: Call the concentration change rate data in the main path set, set 40% as the screening threshold, extract the set of units with values exceeding the threshold, synchronously call the air pressure linkage correction coefficient, use the first derivative to judge the continuous time series change direction, take the intersection of the units meeting the positive growth condition and the units with concentration exceeding the threshold, and generate a risk unit marking set.

[0034] As a further aspect of the present invention, the method further includes:

[0035] S5: Based on the position coordinates of the risk units, execute local area prediction path update using the Kalman filtering algorithm within the radius expansion area, input the updated path nodes into the pollen diffusion model, and output the corrected pollen concentration prediction value and short-term high concentration warning signal.

[0036] As a further aspect of the present invention, the steps for obtaining the corrected pollen concentration prediction value and short-term high concentration warning signal are specifically as follows:

[0037] S501: Based on the position coordinates of the risk units and the radius expansion area, call the position parameters and velocity components of the existing path nodes, combine the weight distribution of the wind speed gradient and terrain resistance parameters, substitute the position deviation value and velocity correction term into the state transition equation, perform recursive update of the error covariance matrix, and correct the positions of the path nodes by superimposing the Kalman gain coefficient to generate a dynamic path node set;

[0038] S502: Extract the coordinate time series and spatial distribution characteristics of the dynamic path node set, perform linear fitting on the node spacing and diffusion rate, verify the continuity of the diffusion direction according to the atmospheric stability coefficient, perform grid integration on the pollen diffusion path using the Euler method, and iteratively correct the concentration gradient by superimposing the sedimentation rate and resistance parameters to generate a regional predicted concentration distribution;

[0039] S503: According to the regional predicted concentration distribution, count the number of grid units with concentration values exceeding the terrain elevation correction threshold, calculate the product relationship between the grid coverage area and the pollen residence time, perform probability density matching on the fluctuation amplitude and the preset warning threshold, and generate the corrected pollen concentration prediction value and short-term high concentration warning signal.

[0040] A pollen concentration prediction system based on big data analysis, the pollen concentration prediction system based on big data analysis is used to execute the above-mentioned pollen concentration prediction method based on big data analysis, and the system includes:

[0041] A terrain grid processing module, which is used to generate a two-dimensional grid through a regional digital elevation model, calculate the elevation difference between the center points of adjacent units, generate a unit distance elevation mutation rate based on the grid spacing, construct an elevation mutation rate matrix and input it into the pollen diffusion model;

[0042] A path node generation module, which is used to call the pollen diffusion model to output a path node sequence, extract the elevation mutation rate matrix of adjacent node pairs, accumulate the elevation mutation rate per unit distance and divide it by the total number of nodes to generate a topographic barrier frequency factor, and screen node pairs with an elevation mutation rate less than 50% of the original path segment and a wind direction difference angle less than 30 degrees to generate an optimized path node combination;

[0043] A trend interpolation module, which is used to input the optimized path node combination into a cubic spline interpolation method to fit a path, calculate the direction difference angle between the fitted path and the original path, and screen paths with a difference angle less than 15 degrees to generate a main path set;

[0044] A risk calibration module, which is used to calculate the slope aspect angle difference based on the difference between the slope normal angle and the main wind direction angle, generate a wind pressure linkage correction coefficient in combination with the wind pressure time series data, and mark the unit with a concentration change rate exceeding 40% and a positive increasing correction coefficient in the main path set as the risk unit coordinate;

[0045] A prediction update module, which is used to expand the radius area based on the risk unit coordinates, update the local prediction path using the Kalman filter algorithm, input the updated path nodes into the pollen diffusion model, and output the corrected pollen concentration prediction value and early warning signal.

[0046] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0047] In the present invention, by converting the terrain elevation change into a quantifiable mutation rate parameter, a dynamic association model between terrain barriers and pollen diffusion paths is established. Based on the topographic barrier frequency factor generated from the elevation difference of grid cells, combined with the wind direction difference angle, path nodes are screened to eliminate the interference of terrain mutations and wind direction offsets on the diffusion path. The cubic spline interpolation method is used to reconstruct the path trend line, and the main path set is dynamically corrected through the direction difference angle threshold, enhancing the model's correction ability for terrain-induced path distortion. The linkage correction coefficient integrating the slope normal angle and wind pressure data is used to establish a coupling analysis mechanism for terrain slope aspect and wind pressure changes, accurately identifying areas with abnormal concentration fluctuations. The multi-dimensional fusion of topographic barrier frequency, slope aspect angle difference, and dynamic wind pressure parameters improves the modeling accuracy of the pollen diffusion path under complex terrains, enabling the prediction results to simultaneously reflect the non-linear effects of terrain barrier effects and real-time meteorological parameters. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 It is a schematic diagram of the working process of the present invention;

[0049] Figure 2 It is a flowchart of the acquisition steps of the pollen diffusion model of the present invention;

[0050] Figure 3Flow chart of the acquisition steps of the optimized path node combination of the present invention;

[0051] Figure 4 Flow chart of the acquisition steps of the main path set of the present invention;

[0052] Figure 5 Flow chart of the acquisition steps of the risk unit of the present invention;

[0053] Figure 6 Flow chart of the acquisition steps of the corrected pollen concentration prediction value and short-term high-concentration warning signal of the present invention. Detailed implementation manners

[0054] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention, and are not used to limit the present invention.

[0055] In the description of the present invention, it should be understood that the orientation or positional relationship indicated by the terms "length", "width", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the accompanying drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the present invention. In addition, in the description of the present invention, "a plurality of" means two or more, unless otherwise specifically defined.

[0056] Embodiment 1

[0057] Please refer to Figure 1 , the present invention provides a technical solution: a pollen concentration prediction method based on big data analysis, including the following steps:

[0058] S1: Perform two-dimensional grid division through regional digital elevation model data, calculate the elevation difference between the central points of adjacent grid cells, generate the elevation mutation rate per unit distance based on the grid spacing, construct an elevation mutation rate matrix, and input the elevation mutation rate matrix into the pollen diffusion model;

[0059] S2: Call the pollen diffusion model to output a path node sequence, extract the elevation mutation rate per unit distance of adjacent node pairs, accumulate and divide by the total number of nodes to generate a topographic barrier frequency factor, and screen node pairs that satisfy the elevation mutation rate less than 50% of the original path segment and the wind direction difference angle less than 30 degrees to generate an optimized path node combination;

[0060] S3: Input the optimized path node combination into the cubic spline interpolation method to fit and generate a diffusion path trend line, calculate the direction difference angle between the fitted path and the original path, and add the paths with a difference angle less than 15 degrees to the main path set;

[0061] S4: Calculate the slope aspect difference angle through the difference between the slope normal angle of the grid cell and the main wind direction angle, combine with the wind pressure time series data to generate a wind pressure linkage correction coefficient, and mark the units in the main path set with a concentration change rate exceeding 40% and a positive increasing wind pressure linkage correction coefficient as risk units;

[0062] S5: Based on the position coordinates of the risk units, use the Kalman filter algorithm to perform local area prediction path update within the radius expansion area, input the updated path nodes into the pollen diffusion model, and output the corrected pollen concentration prediction value and short-term high concentration warning signal.

[0063] The elevation mutation rate matrix specifically includes elevation difference, grid spacing, and elevation mutation rate per unit distance. The optimized path node combination includes the selected node pairs and corrected concentration distribution parameters. The main path set specifically refers to the diffusion path trend line, direction difference angle determination result, and spatio-temporal diffusion trajectory simulation parameters. The risk units include high concentration change rate units and positive increasing wind pressure linkage correction coefficient units.

[0064] Please refer to Figure 2 , and the steps for obtaining the pollen diffusion model are specifically as follows:

[0065] S101: Obtain the regional digital elevation model data, set the grid spacing parameter, divide the region into two-dimensional grid cells, extract the plane coordinates of the center point of each grid cell, and calculate the elevation value of the grid cell center point according to the spatial correspondence between the elevation point coordinates and the grid center point to generate a grid cell elevation data set;

[0066] First, obtain the digital elevation model (DEM) data of the specified research area. The geographical scope of this area is set from 116°30'0” to 116°45'0” east longitude and from 40°00'0” to 40°15'0” north latitude. This DEM data is provided in raster file format, recording the altitude of each point on the ground, and its spatial resolution is determined to be 30 meters. Based on this data, set the spacing parameter of the two-dimensional grid to 100 meters. This setting aims to balance the requirements of calculation efficiency and terrain detail expression, and divide the entire research area into multiple 100-meter × 100-meter square grid cells.

[0067] Subsequently, systematically traverse all the divided grid cells. For each grid cell, accurately calculate the planar coordinates (X, Y) of its geometric center point, using the Universal Transverse Mercator (UTM) projection coordinate system. For the starting grid cell (index 0, 0) at the southwest corner of the study area, the calculated planar coordinates of its center point are (450050E, 4428050N). Next, based on the precise geographic coordinates (x, y, z) of each elevation point recorded in the DEM data, find the DEM elevation points that are spatially adjacent to the current grid center point (X, Y). Using the Inverse Distance Weighting (IDW) interpolation method, calculate the elevation value Z corresponding to the grid center point (X, Y) according to the coordinates and elevation values (z) of the adjacent DEM elevation points and their distances from the grid center point. Specifically, if the four closest DEM elevation points around a certain grid center point are

[0068] P1(450040, 4428060, 150m), P2(450065, 4428070, 152m), P3(450035, 4428030, 155m), P4(450055, 4428025, 153m), by calculating the distances from each point to the center point (450050, 4428050) and assigning corresponding weights, the interpolated elevation value of this center point is 152.8 meters.

[0069] Repeat the process of extracting the center point coordinates and calculating the elevation value for each grid cell in the study area. After processing all the grid cells, integrate the results into a structured data set. This data set contains the row and column indices (i, j) of each grid cell, the center point X coordinate, the center point Y coordinate, and the calculated center point elevation Z. This data set, namely the grid cell elevation data set, provides the basic data support for subsequent terrain analysis.

[0070] S102: Call the grid cell elevation data set, traverse the row and column indices of the grid cells, extract the four adjacent grid cells in the upper, lower, left, and right directions, calculate the absolute value of the elevation difference between the current grid center point and the adjacent grid center points, and divide the elevation difference in each direction by the grid spacing parameter to obtain the elevation mutation rate per unit distance;

[0071] Call the grid cell elevation data set. This data set stores each grid cell identified by the row and column indices (i, j) and its center point elevation value Z(i, j). The program starts the traversal process of all grid cells and accesses each cell through the row and column indices.

[0072] For the currently processed grid cell (i, j), the system identifies and extracts the grid cells in the four basic directions (up, down, left, and right) that are directly adjacent to it spatially: the upper cell (i - 1, j), the lower cell (i + 1, j), the left cell (i, j - 1), and the right cell (i, j + 1). Retrieve the elevation values of the centers of these adjacent cells from the dataset: Z(i - 1, j), Z(i + 1, j), Z(i, j - 1), Z(i, j + 1). Then, calculate the absolute value of the elevation difference between the elevation Z(i, j) of the center of the current grid cell and the elevations of the centers of these four adjacent grid cells. Calculate the absolute value of the elevation difference |ΔZ up | = |Z(i, j) - Z(i - 1, j)|. If the elevation Z(1, 1) of the current cell (1, 1) is 152.8 meters and the elevation Z(0, 1) of its upper cell (0, 1) is 160.2 meters, then |ΔZ up | = |152.8 - 160.2| = 7.4 meters. Similarly, calculate the lower |ΔZ down | = |Z(1, 1) - Z(2, 1)|, the left |ΔZ left | = |Z(1, 1) - Z(1, 0)|, and the right |ΔZ right | = |Z(1, 1) - Z(1, 2)| for the absolute value of the elevation difference.

[0073] Subsequently, divide the absolute value of the elevation difference calculated in each direction by the grid spacing parameter (100 meters) set in S101. This operation aims to obtain the elevation mutation rate per unit distance without units, reflecting the degree of vertical elevation change within a unit horizontal distance. Taking the upper direction as an example, its mutation rate Rate up = |ΔZ up | / 100 meters = 7.4 meters / 100 meters = 0.074. Perform the same calculation on the absolute values of the elevation differences for the lower, left, and right directions to obtain the corresponding mutation rates Rate down , Rate left , Rate right . Repeat this process for all non-boundary grid cells in the study area to calculate the elevation mutation rate per unit distance of each cell relative to its four neighbors. The value range distribution of the elevation mutation rate is used for subsequent analysis. Values below 0.01 indicate flat terrain, values between 0.01 and 0.1 indicate gentle to medium slopes, and values above 0.1 indicate steep terrain.

[0074] S103: Call the elevation mutation rate per unit distance, arrange the data in the order of grid row and column indices, fill the mutation rate values of the missing boundary grids with zero, construct a determinant-structured matrix data, generate an elevation mutation rate matrix, and transfer it to the pollen diffusion model.

[0075] Call the elevation mutation rate data per unit distance in four directions of each grid cell. Systematically organize and arrange these data strictly in the order of the row and column indices (i, j) of the grid.

[0076] Considering the grid cells located at the boundary of the study area, they will lack adjacent cells in one or two directions. Specifically, the cells in the first row have no upper adjacent cells, the last row has no lower cells, the first column has no left cells, and the last column has no right cells. For the missing neighbor directions of these boundary cells, the corresponding elevation mutation rate values are deterministically set to zero. For example, for the grid cell (0, 0), the mutation rates Rate up and Rate left are both assigned the value of 0. This processing method is applied to all boundary cells to ensure the integrity and consistency of the subsequent processed data structure.

[0077] Then, construct one or more structured matrix data. The dimensions of these matrices exactly match the number of grids divided in the study area. One implementation is to construct four independent matrices to store the mutation rate values of all cells in specific directions (up, down, left, right): M up , M down , M left , M right . In the M up matrix, the element M up [i, j] stores the up-direction mutation rate of the cell (i, j), and all element values in the first row of this matrix are 0. Similarly, generate the mutation rate matrices for the other three directions. These matrices together constitute the complete elevation mutation rate matrix data. Finally, generate a structured elevation mutation rate matrix and pass this matrix data as a key factor of terrain influence to the core calculation module of the pollen diffusion model.

[0078] Please refer to Figure 3 , the specific steps for obtaining the optimized path node combination are as follows:

[0079] S201: Call the sequence of path nodes output by the pollen diffusion model, extract the elevation difference and straight-line distance of all adjacent node pairs, divide the absolute value of the elevation difference by the straight-line distance, calculate the elevation mutation rate per unit distance of each node pair, summarize all the mutation rate values, and generate a set of elevation mutation rates of adjacent node pairs;

[0080] Call a specific sequence of path nodes simulated and output by the pollen diffusion model. This sequence consists of a series of nodes arranged in the order of diffusion time, and each node records its three-dimensional spatial coordinates (X, Y, Z). Select one of the paths P, and its node sequence is

[0081] P1(450100,4428100,150),P2(450180,4428160,155),P3(450250,4428200,153),P4(450330,4428230,162),P5(450400,4428280,160),P6(450450,4428350,165)(coordinate unit: meter).

[0082] The program iterates over all pairs of adjacent nodes in this sequence:

[0083] (P1, P2), (P2, P3), (P3, P4), (P4, P5), (P5, P6). For each pair of adjacent nodes, such as (P1, P2), first extract their elevation values z1 = 150 meters and z2 = 155 meters, and calculate the absolute value of the elevation difference between the two nodes |ΔZ1| = |155-150| = 5 meters. Next, extract the plane coordinates of the two nodes (x1, y1) = (450100, 4428100) and (x2, y2) = (450180, 4428160), and calculate the two-dimensional plane straight-line distance between them.

[0084] Then, the calculated absolute value of the elevation difference |ΔZ1| is divided by the calculated straight-line distance Dist1 to obtain the unit distance elevation mutation rate Rate1 of this pair of adjacent nodes = |ΔZ1| / Dist1 = 5 meters / 100 meters = 0.05. Repeat this calculation process for all adjacent node pairs in the path sequence (P2, P3), (P3, P4), (P4, P5), (P5, P6) to obtain their respective mutation rates: These calculated mutation rate values are summed up to form the adjacent node pair elevation mutation rate set of the path {0.05, 0.0248, 0.1053, 0.0232, 0.0581}.

[0085] S202: Based on the elevation mutation rate set of adjacent node pairs, all mutation rate values are accumulated, the accumulated result is divided by the total number of nodes in the node sequence, the ratio of the total mutation rate to the number of nodes is calculated, and the topographic barrier frequency factor is generated;

[0086] The elevation mutation rate set of adjacent nodes based on path P

[0087] {0.05, 0.0248, 0.1053, 0.0232, 0.0581}. First, perform a cumulative summation operation on all mutation rate values in the set, and calculate SumRate = 0.05 + 0.0248 + 0.1053 + 0.0232 +

[0088] 0.0581 = 0.2614.

[0089] Next, obtain the total number of nodes N in the path node sequence. The path P contains a total of N = 6 nodes from P1 to P6. Then, divide the total sum of mutation rates SumRate calculated by the total number of nodes N in the node sequence, calculate this ratio, and define it as the topographic barrier frequency factor Factor = SumRate / N = 0.2614 / 6 ≈ 0.0436. This factor quantifies the average level of the overall terrain undulation of the path relative to the number of path nodes (representing a certain degree of path length or complexity).

[0090] S203: Traverse the set of elevation mutation rates of adjacent node pairs, compare the mutation rate of each node pair with 50% of the topographic barrier frequency factor, filter out the node pairs with mutation rates less than the threshold, synchronously extract the wind speed direction vectors of the filtered node pairs, calculate the included angle between the vectors, filter out the node pairs with an included angle exceeding 30 degrees, and reorganize the remaining node pairs in the original order to generate an optimized path node combination.

[0091] Traverse the set of elevation mutation rates of adjacent node pairs of path P

[0092] {0.05, 0.0248, 0.1053, 0.0232, 0.0581}, and call the topographic barrier frequency factor Factor ≈ 0.0436 calculated by S202. Calculate the screening threshold based on this factor. The threshold is set to 50% of the factor value, aiming to screen out the path segments with relatively gentle terrain. The threshold calculation process is: Threshold = 0.5 × Factor = 0.5 × 0.0436 = 0.0218. Compare the mutation rate Rate of each node pair in the set with this threshold Threshold. The screening criterion is Rate i ≤ Threshold. The comparison results are as follows: Rate1(P1, P2) = 0.05 > 0.0218 (filtered out) Rate2(P2, P3) = 0.0248 > 0.0218 (filtered out) i Rate3(P3, P4) = 0.1053 > 0.0218 (filtered out) Rate4(P4, P5) = 0.0232 > 0.0218 (filtered out) Rate5(P5, P6) = 0.0581 > 0.0218 (filtered out) In this example, the mutation rates of all node pairs are greater than 0.0218, and there are no node pairs that meet the conditions. To demonstrate the subsequent steps, adjust the threshold setting logic and use a fixed threshold. This threshold is obtained based on the analysis of a large number of paths and represents the limit of small terrain undulation. Set the fixed threshold Threshold

[0093] fixed ​= 0.03. Re-screening: Rate1(P1, P2) = 0.05 > 0.03 (rejected) Rate2(P2, P3) = 0.0248 <= 0.03 (keep P2 - P3) Rate3(P3, P4) = 0.1053 > 0.03 (rejected) Rate4(P4, P5) = 0.0232 <= 0.03 (keep P4 - P5) Rate5(P5, P6) = 0.0581 > 0.03 (rejected). The node pairs retained after screening are (P2, P3) and (P4, P5).

[0094] Next, process these screened node pairs and extract their corresponding wind speed direction vectors. The wind speed direction is obtained from the coupled meteorological model corresponding to each node position and time. Obtain the wind direction vectors V wind2 and V wind4 at nodes P2 and P4. At the same time, calculate the connection vectors of the screened node pairs themselves: Vec 23 (from P2 to P3) and Vec 45 (from P4 to P5). Set an angle threshold between vectors AngleThreshold = 30°. This threshold is based on the following considerations: maintaining a certain direction continuity of the path and avoiding overly abrupt turns. 30 degrees is an acceptable range of direction changes in typical mesoscale diffusion. Now check the potential continuous segments formed by the retained nodes in the original sequence. Since the retained are P2 - P3 and P4 - P5, which are not directly adjacent in the original sequence, there is no continuous segment for which the angle needs to be calculated. If the screening result is (P2, P3) and (P3, P4), then calculate the angle between Vec 23 and Vec 34 . If the angle is greater than 30°, then process according to the rules (for example, reject the latter node pair P3 - P4 and keep P2 - P3). In this example, since there is no continuous segment, no angle rejection is required.

[0095] Finally, re-organize the remaining node pairs in the order of the nodes in the original path sequence. Since only the non - continuous (P2, P3) and (P4, P5) are retained, they form two independent path segments. The final optimized path node combination generated is {Segment 1: [P2, P3], Segment 2: [P4, P5]}. If the screening result is P2 - P3 - P4, then the combination is [P2, P3, P4].

[0096] Please refer to Figure 4 for the specific steps to obtain the main path set:

[0097] S301: Call the optimized path node combination, arrange the node coordinate sequences in the input order, and use the cubic spline interpolation method to calculate the cubic polynomial coefficients between adjacent nodes to ensure that the function is continuous and second-order differentiable at the nodes, construct a piecewise cubic polynomial function, and generate a continuously differentiable diffusion path trend line;

[0098] Call the optimized path node combination. Take one segment as an example, the sequence is P2'(450180, 4428160, 155), P3'(450250, 4428200, 153). Arrange these node coordinates in the input order. Use the cubic spline interpolation method to construct a smoothly connected curve segment between adjacent nodes.

[0099] For the node pair (P2', P3'), find a parametric cubic polynomial curve segment S1(t) = (X(t), Y(t), Z(t)), where X(t) = a x1 t 3 + b x1 t 2 + c x1 t + d x1 , Y(t) and Z(t) are similar, and the parameter t varies from 0 to 1 (or use the cumulative chord length as the parameter). This curve segment must pass exactly through the starting point and the ending point: S1(0) = P2′ and S1(1) = P3′. If there are more nodes, such as P2', P3', P4', then at the internal node P3', it is required that the curve segment S1(t) connecting P2'-P3' and the curve segment S2(t) connecting P3'-P4' have continuous first-order derivatives (the tangent vectors are the same, S′1(t P3′ ) = S′2(t P3′ )) and continuous second-order derivatives (the curvature changes smoothly, S″1(t P3′ ) = S″2(t P3′ ))). These continuity conditions and the endpoint conditions together form a system of linear equations. Solve this system of equations (usually in combination with boundary conditions, such as setting the natural spline condition that the second-order derivative at the endpoints is zero) to obtain all the coefficients of each cubic polynomial (such as a x1 , b x1 , c x1 , d x1 ,..).

[0100] Connect all the calculated piecewise cubic polynomial functions S i (t) in sequence to form a smooth curve that passes through all the optimized nodes and is continuous in position, tangent, and curvature (except at the endpoints). This curve defines a continuously differentiable diffusion path trend line, providing a geometric basis for subsequent analysis.

[0101] S302: Based on the diffusion path trend line, extract the original path tangent direction angle at equally spaced sampling points on the path, calculate the difference in the tangent direction angles of the fitting paths at the corresponding positions, introduce the path curvature difference, sampling point density, and path length scale factor, and use the formula:

[0102]

[0103] Perform operations to obtain the curvature-distance weighted angle offset of multiple sampling points, accumulate them, and normalize them in combination with the path length factor to generate the direction difference angle;

[0104] where, θ new represents the direction difference angle, Δα i represents the difference in the tangent angles of the original path and the fitting path at the i-th sampling point, represents the Euclidean distance of the path segment between the i-th sampling point and the previous node, K i represents the absolute value of the curvature difference between the original path and the fitting path at the i-th sampling point, D i represents the unit length sampling density of the path segment where the i-th sampling point is located, R represents the ratio factor of the current path total length to the reference path length, and m is the total number of sampling points;

[0105] Based on the continuously differentiable diffusion path trend line (cubic spline curve), sample along this curve at equal arc length intervals, and set the sampling interval to 5 meters. This interval selection aims to capture the detailed changes in the path while controlling the computational amount. Obtain a series of sampling points s1, s2,..., s m . Assume that the calculated total path length is L total = 180.62 meters (including two segments P2 - P3 and P4 - P5, with lengths of 80.62 meters and 86.02 meters respectively. Assume only the P2 - P3 segment with a length of 80.62 meters is analyzed, then m = floor(80.62 / 5)+1 = 17 sampling points).

[0106] For the i-th sampling point s i , first determine its position on the original optimized path segment (such as P2 - P3), and obtain the tangent direction angle α orig,i of this original path segment, that is, the direction angle of the vector Vec 23 . Then, calculate the tangent direction angle α i of the cubic spline curve at the same point s fit,i . This angle is obtained by calculating the first derivative vector (X′(t), Y′(t), Z′(t)) of the spline function at the parameter value corresponding to s i , and then calculating the horizontal direction angle of this vector. Calculate the difference Δα i = α orig,i - α fit,i , and take its absolute value |Δαi |.

[0107] Meanwhile, calculate the curve distance of the i-th sampling point s i from the starting node (P2) of its affiliated path segment Calculate at the sampling point s i the curvature difference K between the original straight line segment (with curvature 0) and the fitted spline curve i = |0 - Curvature fit,i | = |Curvature fit,i |. The curvature of the fitted curve is calculated from the second derivative of the spline function. Determine the sampling density D per unit length of the path segment where the i-th sampling point is located i . Since it is equi-arc-length sampling, the number of sampling points per unit length (meter) is 1 / 5 = 0.2 points / meter, so D i = 0.2. Set a reference path length L base , which is based on the statistical analysis of the lengths of multiple simulated paths, taking the average value or representative length, and set L base = 200 meters. Calculate the total length L of the currently analyzed path segment (P2 - P3) total = 80.62 meters and the scale factor R of the reference length = L total / L base = 80.62 / 200 ≈ 0.403.

[0108] Substitute these parameters into the direction difference angle calculation formula:

[0109] Formula parameter assignment and calculation process (taking m = 3 sampling points as an example for simplified demonstration): To ensure the calculation is executable, all parameters involved in the operation need to use compatible expressions. The angle Δα i is in degrees (°), the distance is in meters (m), the curvature K i is in 1 / meter (m-1), the density D i is in points / meter (m-1), and the scale factor R is dimensionless. The input K of the exponential function i usually needs to be dimensionless, but here K i represents the curvature and has the unit of m-1. To be strictly dimensionless, it needs to be multiplied by a characteristic length, but in this formula structure, K i appears directly in the exponent, probably intending to express the influence degree of the curvature magnitude on the weight, and the value with the unit of m-1 is directly substituted numerically. The unit of the final result θ new will inherit from Δα i , that is, degrees (°).

[0110] Suppose the data of three sampling points calculated are as follows: Point

[0111] 1 (s1): |Δα1| = 1.5°, K1 = 0.002 m -1 , D1 = 0.2 m -1 Point

[0112] 2 (s2): |Δα2| = 0.8°, K2 = 0.001 m -1 , D2 = 0.2 m -1 Point

[0113] 3 (s3): |Δα3| = 2.1°, K3 = 0.003 m -1 , D3 = 0.2 m -1 Total length L of the path segment total = 80.62 m, R = 0.403.

[0114] Calculate the numerator (sum of weighted angular differences): Term 1: 1.5 × 5 × e 0.002 ≈ 7.5 × 1.002002 =

[0115] 7.515 Term 2: 0.8 × 10 × e 0.001 ≈ 8.0 × 1.001000 = 8.008 Term 3: 2.1 × 15 × e 0.003 ≈

[0116] 31.5 × 1.003004 = 31.595 Numerator ∑ = 7.515 + 8.008 + 31.595 = 47.118

[0117] Calculate the square root term in the denominator (square root of the sum of weighted squared distances): Term 1: 0.2 × 5 2 = 0.2 × 25 = 5 Term 2: 0.2 × 10 2 = 0.2 × 100 = 20 Term 3: 0.2 × 15 2 = 0.2 × 225 = 45 Sum inside the square root = 5 + 20 + 45 = 70 Square root term

[0118] Calculate the denominator: Denominator = square root term × R = 8.367 × 0.403 ≈ 3.372

[0119] Calculate θ new : θ new = Numerator / Denominator = 47.118 / 3.372 ≈ 13.97°

[0120] The benefit of the formula is that it takes into account the distance between sampling points on the path segment and the curvature difference K iFor the angular difference |Δα i |, it is weighted and combined with the sampling density D i and the total path length ratio R for normalization, so as to not only reflect the average angular deviation, but more comprehensively quantify the overall directional difference in geometric shape between the fitted path and the original (optimized) path segment, especially emphasizing the deviation contributions of long-distance segments and regions with high curvature changes.

[0121] This calculation finally obtains the directional difference angle θ new ≈13.97°. This value of 13.97° represents the overall directional deviation degree between the cubic spline fitting curve of this path segment (P2 - P3) and its original straight line segment, integrating local angular deviation, geometric features (distance, curvature) and path scale factors. Repeat this calculation for all optimized path segments (such as P4 - P5).

[0122] S303: Based on the directional difference angle, set a preset angle threshold as the screening criterion, traverse all diffusion paths, extract the improved directional difference angle values of each path, compare the directional difference angle values with the preset angle threshold, screen out the paths that meet the conditions and store them in a set to generate the main path set.

[0123] According to the directional difference angle θ new of each diffusion path (or path segment). Set a preset angle threshold Threshold θ as the criterion for screening the main path. The setting of this threshold refers to the comparison and analysis of historical simulation data and the measured pollen diffusion range, and the results show that the paths with θ new lower than 15° have a strong correlation with the main diffusion channels, so select Threshold θ = 15°.

[0124] The system traverses all path segments for which θ new has been calculated. Suppose the directional difference angles of two calculated path segments are θ new,P2P3 = 13.97° and θ new,P4P5 = 16.5° respectively. Compare the θ new value of each path segment with the preset threshold Threshold θ = 15°. The screening condition is θ new ≤ Threshold θ . For the path segment P2 - P3: 13.97° ≤ 15°, it meets the condition and is selected. For the path segment P4 - P5: 16.5° > 15°, it does not meet the condition and is excluded.

[0125] Store all the filtered path segments that meet the conditions (in this example, only the P2 - P3 segment) into a new set. This set constitutes the main path set. Generate the main path set {path segment P2 - P3}.

[0126] Please refer to Figure 5 , and the steps for obtaining the risk unit are specifically as follows:

[0127] S401: Obtain the slope normal angle and the main wind direction angle of the grid cell, calculate the cosine value of the included angle between the slope normal vector and the main wind direction vector, use the result of the inverse cosine calculation as the angle difference, take the absolute value of the angle difference to eliminate the influence of direction symmetry, and compare the absolute value result with the preset wind direction deviation reference threshold to generate the slope aspect angle difference of the grid cell;

[0128] Obtain the terrain information of a specific grid cell (i, j) in the study area and the prevailing wind direction at that time. By analyzing the DEM data generated in S101, calculate that the slope β = 10° and the slope aspect α = 135° (southeast slope) of this cell. Obtain the prevailing wind direction angle θ wind = 270° (due west) at the same moment from the meteorological data source.

[0129] Calculate the slope normal vector of this cell and the spatial included angle γ between the horizontal main wind direction vector The slope normal vector is determined by the slope and the slope aspect: The cosine value of the included angle Calculate the included angle γ = arccos(-0.1228). Since the input of the arccos function requires a numerical value, the angle unit needs to be unified. Here, direct numerical calculation is used. γ ≈ 97.06°.

[0130] This included angle γ represents the actual spatial included angle between the slope normal and the horizontal wind direction. Take its absolute value |γ| = 97.06°. Set a wind direction deviation reference threshold Threshold γ . This threshold is set at 90°, representing the state where the wind direction is orthogonal to the slope. Deviation from this reference indicates that the wind has a component parallel to the slope or against / along the slope. Compare the calculated absolute value of the angle difference |γ| with 90°. The deviation degree of 97.06° relative to 90° is |97.06 - 90| = 7.06°. This absolute value of the angle difference |γ| = 97.06° is defined as the slope aspect angle difference Φ of the grid cell for subsequent calculations. Generate the slope aspect angle difference Φ = 97.06° of this grid cell.

[0131] S402: Based on the slope aspect angle difference of the grid cell and the wind pressure time series data, use the formula:

[0132]

[0133] Obtain the dynamic correction factor through calculation, and combine it with the time window mean of the wind pressure time series data to generate the wind pressure linkage correction coefficient;

[0134] Among them, C adj represents the wind pressure linkage correction coefficient, Φ represents the slope aspect angle difference of the grid cell, P t represents the wind pressure value at the t-th moment, represents the time series mean of the wind pressure, n represents the total amount of time series data, H v represents the vertical height change rate of the unit, R s represents the ratio of the wind speeds of adjacent units, β represents the terrain curvature influence factor, S c represents the surface roughness coefficient, τ represents the time decay factor;

[0135] Based on the slope aspect angle difference of the grid cell Φ = 97.06°, and combined with the wind pressure time series data recorded by the meteorological monitoring station near the cell.

[0136] Table 1 Wind pressure time series data table of the monitoring point near a certain unit in the past 24 hours

[0137] Time (t) Wind pressure PtPt (Pa) Time (t) Wind pressure PtPt (Pa) 1 101325 13 101380 2 101330 14 101375 3 101340 15 101370 4 101350 16 101360 5 101360 17 101355 6 101370 18 101350 7 101380 19 101345 8 101390 20 101340 9 101400 21 101335 10 101410 22 101330 11 101405 23 101325 12 101395 24 101320

[0138] As shown in Table 1, the wind pressure monitoring values (unit: Pascal, Pa) in the past 24 hours (n = 24) are provided. The calculation formula of the wind pressure linkage correction coefficient is used: Note: The original formula has been adjusted, and the wind pressure deviation sum has been normalized to make the numerator of the first term a dimensionless value to ensure the validity of adding it to the second term (dimensionless).

[0139] Obtaining and assigning values to the formula parameters: Φ′: A dimensionless factor derived from the slope aspect angle difference Φ, reflecting the influence degree of the relative relationship between the slope surface and the wind direction on the wind pressure. The conversion rule is set as Φ′ = |cos(Φ)|, representing the projection size of the slope normal on the wind direction. Φ = 97.06°, then Φ′ = |cos(97.06°)| = |-0.1228| = 0.1228. P t : The wind pressure value at the t-th moment, from Table 1. The mean of the wind pressure time series data. Calculated from the data in Table 1: n: The total amount of time series data, n = 24. The sum of the normalized wind pressure deviations. Calculate The absolute value is used in the formula, so it is |-0.0001026| = 0.0001026. H v: Unit vertical height change rate. It is calculated by analyzing the DEM elevation values of this unit and its neighboring units, representing the degree of terrain undulation within the unit. The calculated value of H v = 0.05 (unitless, m / m). R s : Ratio of wind speeds of adjacent units. Through the coupled wind field model, the ratio of the wind speed of the current unit to the wind speed of the neighboring unit in the upwind direction is obtained. The calculated value of R s = 0.9m / s / 0.8m / s = 1.125 (unitless). β: Topographic curvature influence factor. It is obtained by calculating the Laplace value (second derivative) of the DEM at this unit, reflecting the terrain convexity and concavity. A positive value indicates a convex landform. The calculated value of β = 0.1 (unitless). S c : Surface roughness coefficient. It is obtained from the standard aerodynamic roughness parameter table according to the land cover type of this unit (such as grassland). The obtained value of S c = 0.05 (unitless). τ: Time decay factor. It represents the decay rate of the influence of past information, and is set to 0.95, indicating that the weight of information decays by 5% after one hour. This value is based on typical estimates of the response time of atmospheric processes. The obtained value of τ = 0.95 (unitless).

[0140] Calculation process: Calculate the denominator of the first term: Calculate the numerator of the first term: (Here, use the absolute value of the original sum, that is It may be more reasonable, that is Adopt the latter for calculation.) Calculate the first term: 0.00103 / 1.2108 ≈ 0.00085 Calculate the second term: Calculate C adj : C adj = 0.00085 + 0.00526 = 0.00611

[0141] The advantage of the formula is that it integrates the interaction between the terrain slope and the wind direction (Φ′), the actual fluctuations of recent wind pressures (∑), the local terrain undulation (H v ), the relative wind speed change (R s ), the topographic curvature (β), and the surface roughness (S c ), and considers the time influence (τ), dynamically generating a comprehensive correction coefficient, thereby more precisely adjusting the calculations related to wind pressure, such as concentration correction or risk assessment.

[0142] The calculated wind pressure linkage correction coefficient C adj≈0.00611. This is a unitless correction factor, and its numerical value reflects the comprehensive adjustment degree of the wind pressure effect relative to the reference state under the current specific terrain, meteorological, and surface conditions.

[0143] S403: Call the concentration change rate data in the main path set, set 40% as the screening threshold, extract the set of cells with values exceeding the threshold, synchronously call the wind pressure linkage correction coefficient, use the first derivative to judge the continuous time series change direction, and take the intersection of the cells meeting the positive growth condition and the cells with concentration exceeding the threshold to generate a risk cell marking set.

[0144] Call the pollen concentration change rate data of the grid cells passed by the main path set. These data are provided by the prediction of the pollen diffusion model and represent the increase in concentration per unit time (hour) (unit: micrograms per cubic meter per hour, μg / m 3 / h). Set the screening threshold of the concentration change rate. This threshold is determined based on the retrospective analysis of historical high-concentration events. When the concentration change rate exceeds a specific percentage of the regional average rate, the risk increases significantly. Set the threshold as 40% of the maximum observed rate. If the maximum rate observed in the region is 10 μg / m 3 / h, then the threshold Threshold rate = 0.40 × 10 = 4 μg / m 3 / h.

[0145] Traverse all the grid cells passed by the main path segment P2 - P3. Extract the current concentration change rate value of each cell. Suppose the path passes through cell A and cell B, and their rates are 5 μg / m 3 / h and 3 μg / m 3 / h respectively. Compare these rates with the threshold Threshold rate . For cell A: 5 μg / m 3 / h > 4 μg / m 3 / h, which meets the condition and is selected. For cell B: 3 μg / m 3 / h < 4 μg / m 3 / h, which does not meet the condition. Form a set of cells with high concentration change rates {A}.

[0146] Meanwhile, obtain the time series data of the wind pressure linkage correction coefficient C adj (from the continuous calculation in S402). It is necessary to analyze the recent change trend of C adj . Calculate the change rate of C adj using the first-order backward difference: ΔC adj (t) = C adj (t) - C adj (t - 1). Judge whether ΔC adj (t) > 0 holds, that is, whether Cadj Whether it is in a positive growth trend. Let C of unit A at the current time t adj (t) = 0.00611, and C at the previous time adj (t - 1) = 0.00590, then ΔC adj (t) = 0.00611 - 0.00590 = 0.00021 > 0. Unit A meets the C adj positive growth condition. Suppose the path also passes through unit C, and its concentration change rate is 2 μg / m 3 / h (< threshold), but C adj is also growing (ΔC adj (t) = 0.0003 > 0). A set of positive growth units of the correction coefficient {A, C} is formed.

[0147] Finally, perform the intersection operation of the set of high concentration change rate units {A} and the set of positive growth units of the correction coefficient {A, C}. The set of units that meet both conditions {A} is obtained. These units are marked as risk units. A risk unit marking set {A} is generated.

[0148] Please refer to Figure 6 , and the specific steps for obtaining the corrected pollen concentration prediction value and the short-term high concentration warning signal are as follows:

[0149] S501: Based on the position coordinates and radius expansion area of the risk unit, call the position parameters and velocity components of the existing path nodes, combine the weight distribution of the wind speed gradient and terrain resistance parameters, substitute the position deviation value and velocity correction term into the state transition equation, perform the recursive update of the error covariance matrix, and superimpose the Kalman gain coefficient to correct the position of the path nodes to generate a dynamic path node set;

[0150] Based on the risk unit marking set {A}. Obtain the position coordinates (the center point coordinates) of risk unit A and the radius defining its influence range. This radius is determined according to the risk level and diffusion scale and is set to R = 50 meters. Call the position parameters of the nodes (i.e., P2, P3) of the path segment P2 - P3 in the main path set S303 and the velocity components calculated by the diffusion model.

[0151] Combine the wind speed gradient information (obtain the wind speed data of different height layers from the meteorological model and calculate the vertical gradient) and the terrain resistance parameters (the resistance coefficient calculated based on the surface coverage type and terrain complexity around unit A and itself). Assign weights to these two factors. The weights w gradient and w resistance are determined based on sensitivity analysis, reflecting their relative dominant role in particle movement in this area. Set w gradient = 0.6 and w resistance= 0.4. Calculate the position deviation value caused by the position information of risk unit A (such as the deviation vector between the predicted node position and the center of unit A), and the velocity correction term ΔV calculated from the wind speed gradient and terrain resistance (combined by weight).

[0152] Integrate the position deviation value and the velocity correction term into the state transition equation of the Kalman filter as external inputs or correction quantities. The Kalman filter process includes: predicting the node state (position, velocity) and the error covariance matrix P at the next moment; calculating the Kalman gain K, which balances the uncertainty of the model prediction (from P) and the uncertainty of the observation (risk unit position information); using the risk unit position information as the "observation value" to update the state estimate (position and velocity) of the node and the error covariance matrix. Perform this recursive update step of the Kalman filter on the nodes (or all nodes) near risk unit A on the path segment P2 - P3. Through this process, the positions of the nodes are dynamically corrected to reflect the special influence of the risk unit area. Generate a dynamic path node set containing the corrected node positions and velocities.

[0153] S502: Extract the coordinate time series and spatial distribution characteristics of the dynamic path node set, linearly fit the node spacing and the diffusion rate, verify the continuity of the diffusion direction based on the atmospheric stability coefficient, perform grid integration on the pollen diffusion path using the Euler method, and iteratively correct the concentration gradient by superimposing the sedimentation rate and the resistance parameter to generate the regional predicted concentration distribution;

[0154] Extract the time series coordinate data P′ k (t) and the spatial distribution characteristics of the dynamic path node set. Analyze the variation relationship of the node spacing d(t) and the diffusion rate v(t) (calculated from the node displacement) over time, perform a linear fit v(t) = a·d(t) + b to obtain the coefficients a and b, which are used to characterize the diffusion characteristics.

[0155] Judge the current atmospheric stability level based on meteorological data (such as wind speed, solar radiation, cloud cover), and use the Pasquill - Gifford classification method to determine it as class B (unstable). Consult the atmospheric stability parameter table, and the quantization coefficient S corresponding to class B stab = 1.5. Use this coefficient S stab to verify the direction continuity of the dynamic path node sequence, and check whether the angular distribution of the adjacent node displacement vectors conforms to the typical diffusion cone angle range under class B stability.

[0156] The Eulerian-based gridded advection-diffusion integration model is adopted to calculate the spatio-temporal distribution of pollen concentration within the study area (the grid defined by S101). At each time step Δt, the concentration C(i,j,t+Δt) of the grid cell (i,j) is calculated. The calculation is based on C(i,j,t) and the fluxes into and out of the cell calculated through a dynamic set of path nodes (reflecting the modified advection paths and velocities). The formula structure is C(t+Δt) = C(t) + Δt·(Advection + Diffusion + Sources / Sinks). When calculating the source / sink term, the settling velocity V settling and the terrain / vegetation resistance parameter R drag of pollen are superimposed. The settling velocity V settling is obtained from literature or experimental data according to the physical properties (particle size, density) of the target pollen type (such as ragweed pollen), and is determined to be V settling = 0.02 m / s. The resistance parameter R drag is calculated based on the surface roughness of the grid cell (from S c ) and the terrain complexity. The settling effect is manifested as the removal of pollen from the bottom of the computational domain, and the resistance is achieved by modifying the effective wind speed in the flux calculation.

[0157] By iteratively performing this gridded integration calculation over multiple time steps and continuously correcting the concentration gradient according to the dynamically changing node set and environmental parameters, a refined spatial distribution map C(i,j) of pollen concentration in the study area at the target prediction time is finally generated.

[0158] S503: According to the regional predicted concentration distribution, count the number of grid cells whose concentration values exceed the terrain elevation correction threshold, calculate the product relationship between the grid coverage area and the pollen residence time, and perform probability density matching between the fluctuation amplitude and the preset warning threshold to generate a corrected pollen concentration prediction value and a short-term high-concentration warning signal.

[0159] Based on the regional predicted concentration distribution map C(i,j). Set a terrain elevation correction concentration threshold Threshold conc for preliminary screening of high-concentration areas. This threshold is set based on the health risk assessment criteria and the local environmental background value, and is determined to be 50 μg / m 3 . Count the number N conc of grid cells in the region where the concentration value C(i,j) exceeds Threshold high . Calculate the total area Area high covered by these high-concentration cells = N high ×(100m×100m).

[0160] Meanwhile, combining wind field and terrain data, estimate the average residence time T of pollen in these high-concentration areasresidence This estimate takes into account local circulation and sedimentation processes to obtain T residence = 30 minutes. Calculate the cumulative exposure risk index ExposureIndex = Area high ×T residence .

[0161] Analyze the time series of predicted concentrations and calculate the concentration fluctuation amplitude (such as standard deviation or peak concentration). Set a concentration threshold Threshold warning for triggering early warnings. This threshold represents the level that may cause significant health responses and is set to 80 μg / m 3 . Based on the uncertainty analysis of concentration prediction (or multi-model ensemble prediction), calculate the probability P(C > Threshold warning ). Obtain a probability value of P = 0.15. warning

[0162] Set the early warning trigger rule: When P(C > Threshold warning ) is greater than the probability threshold of 0.1 and ExposureIndex is greater than the risk level threshold of 10000 m 2 ·min, trigger an early warning. In this example, P = 0.15 > 0.1. If the calculated ExposureIndex = 12000 m 2 ·min > 10000 m 2 ·min, then both conditions are met.

[0163] Based on the above analysis results (area of high concentration region, residence time, exposure index, concentration volatility, probability of exceeding the threshold), generate the final corrected pollen concentration prediction value (which can be the maximum predicted concentration, average concentration or comprehensive risk index of key regions) and the corresponding short-term high concentration early warning signal (in this example, the early warning signal is "trigger").

[0164] A pollen concentration prediction system based on big data analysis is used to execute the above pollen concentration prediction method based on big data analysis. The system includes:

[0165] A terrain grid processing module for generating a two-dimensional grid through a regional digital elevation model, calculating the elevation difference between the center points of adjacent cells, generating the elevation mutation rate per unit distance based on the grid spacing, constructing an elevation mutation rate matrix and inputting it into the pollen diffusion model;

[0166] ​A path node generation module, which is used to call a pollen diffusion model to output a path node sequence, extract an elevation mutation rate matrix of adjacent node pairs, accumulate the elevation mutation rate per unit distance and divide it by the total number of nodes to generate a topographic barrier frequency factor, and screen node pairs with an elevation mutation rate less than 50% of the original path segment and a wind direction difference angle less than 30 degrees to generate an optimized path node combination;

[0167] A trend interpolation module, which is used to input the optimized path node combination into a cubic spline interpolation method to fit a path, calculate the direction difference angle between the fitted path and the original path, and screen paths with a difference angle less than 15 degrees to generate a main path set;

[0168] A risk calibration module, which is used to calculate the slope aspect angle difference based on the difference between the slope normal angle and the main wind direction angle, generate a wind pressure linkage correction coefficient in combination with the wind pressure time series data, and mark the unit coordinates with a concentration change rate exceeding 40% and a positive increasing correction coefficient in the main path set as risk unit coordinates;

[0169] A prediction update module, which is used to expand the radius area based on the risk unit coordinates, update the local prediction path using the Kalman filter algorithm, input the updated path nodes into the pollen diffusion model, and output the corrected pollen concentration prediction value and warning signal.

[0170] The above are only the preferred embodiments of the present invention, and do not limit the present invention in other forms. Any person skilled in the art may use the disclosed technical content to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention without departing from the technical solution content of the present invention still belong to the protection scope of the technical solution of the present invention.

Claims

1. A pollen concentration prediction method based on big data analysis, characterized in that, It includes the following steps: S1: Perform two-dimensional grid division through regional digital elevation model data, calculate the elevation difference between the central points of adjacent grid cells, generate the elevation mutation rate per unit distance based on the grid spacing, construct an elevation mutation rate matrix, and input the elevation mutation rate matrix into the pollen diffusion model; S2: Call the pollen diffusion model to output a path node sequence, extract the elevation mutation rate per unit distance of adjacent node pairs, accumulate and divide by the total number of nodes to generate a topographic barrier frequency factor, screen node pairs that meet the conditions that the elevation mutation rate is less than 50% of the original path segment and the wind direction difference angle is less than 30 degrees, and generate an optimized path node combination; S3: Input the optimized path node combination into the cubic spline interpolation method, fit to generate a diffusion path trend line, calculate the direction difference angle between the fitted path and the original path, and add paths with a difference angle less than 15 degrees to the main path set; S4: Calculate the slope aspect angle difference through the difference between the slope normal angle of the grid cell and the main wind direction angle, combine with the wind pressure time series data to generate a wind pressure linkage correction coefficient, and mark the units in the main path set with a concentration change rate exceeding 40% and a positive increasing wind pressure linkage correction coefficient as risk units.

2. The pollen concentration prediction method based on big data analysis according to claim 1, wherein The elevation mutation rate matrix specifically includes elevation difference, grid spacing, and elevation mutation rate per unit distance. The optimized path node combination includes the screened node pairs and the corrected concentration distribution parameters. The main path set specifically refers to the fitted path trend line, the direction difference angle determination result, and the spatio-temporal diffusion trajectory simulation parameters. The risk units include high concentration change rate units and positive increasing wind pressure linkage correction coefficient units.

3. The pollen concentration prediction method based on big data analysis according to claim 2, characterized in that, The acquisition steps of the pollen diffusion model are specifically as follows: S101: Obtain regional digital elevation model data, set grid spacing parameters, divide the region into two-dimensional grid cells, extract the plane coordinates of the central points of each grid cell, and calculate the elevation value of the grid center point according to the spatial correspondence between the elevation point coordinates and the grid center point to generate a grid cell elevation data set; S102: Call the grid cell elevation data set, traverse the row and column indexes of the grid cells, extract the adjacent grid cells in the upper, lower, left, and right four directions, calculate the absolute value of the elevation difference between the current grid center point and the adjacent grid center points, and divide the elevation difference in each direction by the grid spacing parameter to obtain the elevation mutation rate per unit distance; S103: Call the elevation mutation rate per unit distance, arrange the data in the order of grid row and column indexes, fill the mutation rate values of the missing boundary grids with zero, construct a determinant structured matrix data, generate an elevation mutation rate matrix, and transfer it to the pollen diffusion model.

4. The pollen concentration prediction method based on big data analysis according to claim 3, wherein The acquisition steps of the optimized path node combination are specifically as follows: S201: Call the path node sequence output by the pollen diffusion model, extract the elevation difference and straight-line distance of all adjacent node pairs, divide the absolute value of the elevation difference by the straight-line distance, calculate the elevation mutation rate per unit distance of each node pair, summarize all the mutation rate values, and generate an adjacent node pair elevation mutation rate set; S202: Based on the elevation mutation rate set of the adjacent node pairs, accumulate all the mutation rate values, divide the accumulated result by the total number of nodes in the node sequence, calculate the ratio of the total mutation rate to the number of nodes, and generate a topographic barrier frequency factor. S203: Traverse the elevation mutation rate set of the adjacent node pairs, compare the mutation rate of each node pair with 50% of the topographic barrier frequency factor, screen out the node pairs with mutation rates less than the threshold, synchronously extract the wind speed direction vectors of the screened node pairs, calculate the included angle between the vectors, eliminate the node pairs with included angles exceeding 30 degrees, and reorganize the remaining node pairs in the original order to generate an optimized path node combination.

5. The pollen concentration prediction method based on big data analysis according to claim 4, wherein The specific steps for obtaining the main path set are as follows: S301: Call the optimized path node combination, arrange the node coordinate sequence in the input order, use the cubic spline interpolation method to calculate the cubic polynomial coefficients between adjacent nodes, ensure that the function is continuous and second-order differentiable at the nodes, construct a piecewise cubic polynomial function, and generate a continuously differentiable diffusion path trend line. S302: Based on the diffusion path trend line, extract the original path tangent direction angle at the equally spaced sampling points on the path, calculate the difference between the fitted path tangent direction angles at the corresponding positions, introduce the path curvature difference, sampling point density, and path length proportionality factor, and use the formula: Operate to obtain the curvature-distance weighted angle offset of multiple sampling points, accumulate and normalize it in combination with the path length factor to generate a direction difference angle. Among them, θ new represents the direction difference angle, Δα i represents the tangent angle difference between the original path and the fitted path at the i-th sampling point, represents the Euclidean distance of the path segment between the i-th sampling point and the previous node, K i represents the absolute value of the curvature difference between the original path and the fitted path at the i-th sampling point, D i represents the unit length sampling density of the path segment where the i-th sampling point is located, R represents the ratio factor of the total length of the current path to the length of the reference path, and m is the total number of sampling points; S303: According to the direction difference angle, set a preset angle threshold as the screening criterion, traverse all the diffusion paths, extract the improved direction difference angle values of each path, compare the direction difference angle values with the preset angle threshold, screen out the paths that meet the conditions and store them in a set to generate a main path set.

6. The pollen concentration prediction method based on big data analysis according to claim 5, characterized in that The specific steps for obtaining the risk unit are as follows: S401: Obtain the slope normal angle and the main wind direction angle of the grid unit, calculate the cosine value of the included angle between the slope normal vector and the main wind direction vector, take the arccosine calculation result as the angle difference, take the absolute value of the angle difference to eliminate the influence of direction symmetry, and compare the absolute value result with the preset wind direction deviation reference threshold to generate the slope aspect angle difference of the grid unit. S402: Based on the slope aspect angle difference of the grid unit and the wind pressure time series data, use the formula: Operate to obtain a dynamic correction factor, and combine it with the time window mean of the wind pressure time series data to generate a wind pressure linkage correction coefficient. Among them, C adj represents the wind pressure linkage correction coefficient, Φ represents the slope direction angle difference of the grid unit, P t represents the wind pressure value at the t-th moment, represents the mean value of the wind pressure time series, n represents the total amount of time series data, H v represents the vertical height change rate of the unit, R s represents the ratio of the wind speed of adjacent units, β represents the terrain curvature influence factor, S c represents the surface roughness coefficient, τ represents the time decay factor; S403: Call the concentration change rate data in the main path set, set 40% as the screening threshold, extract the unit set with values exceeding the threshold, synchronously call the wind pressure linkage correction coefficient, use the first derivative to judge the continuous time series change direction, and take the intersection of the units that meet the positive growth condition and the units with concentrations exceeding the threshold to generate a risk unit mark set.

7. The pollen concentration prediction method based on big data analysis according to claim 6, characterized in that The method further includes: S5: Based on the position coordinates of the risk unit, use the Kalman filter algorithm to perform local area prediction path update within the radius expansion area, input the updated path nodes into the pollen diffusion model, and output the corrected pollen concentration prediction value and short-term high concentration warning signal.

8. The pollen concentration prediction method based on big data analysis according to claim 7, characterized in that The specific steps for obtaining the corrected pollen concentration prediction value and short-term high concentration warning signal are as follows: S501: Based on the position coordinates and radius expansion area of the risk unit, call the position parameters and velocity components of the existing path nodes, combine the weight distribution of the wind speed gradient and terrain resistance parameters, substitute the position deviation value and velocity correction term into the state transition equation, perform recursive update of the error covariance matrix, and superimpose the Kalman gain coefficient to correct the positions of the path nodes to generate a dynamic path node set; S502: Extract the coordinate time series and spatial distribution characteristics of the dynamic path node set, perform linear fitting on the node spacing and diffusion rate, verify the continuity of the diffusion direction based on the atmospheric stability coefficient, perform grid integration on the pollen diffusion path using the Euler method, and iteratively correct the concentration gradient by superimposing the sedimentation rate and resistance parameters to generate a regional predicted concentration distribution; S503: According to the regional predicted concentration distribution, count the number of grid cells whose concentration values exceed the terrain elevation correction threshold, calculate the product relationship between the grid coverage area and the pollen retention time, perform probability density matching on the fluctuation amplitude and the preset warning threshold to generate a corrected pollen concentration prediction value and a short-term high-concentration warning signal.

9. A pollen concentration prediction system based on big data analysis, characterized in that, The system is used to implement the pollen concentration prediction method based on big data analysis according to any one of claims 1-8. The system includes: A terrain grid processing module, which is used to generate a two-dimensional grid through a regional digital elevation model, calculate the elevation difference between the central points of adjacent units, generate a unit distance elevation mutation rate based on the grid spacing, construct an elevation mutation rate matrix and input it into the pollen diffusion model; A path node generation module, which is used to call the pollen diffusion model to output a path node sequence, extract the elevation mutation rate matrix of adjacent node pairs, accumulate the unit distance elevation mutation rate and divide it by the total number of nodes to generate a geomorphic barrier frequency factor, and screen the node pairs whose elevation mutation rate is less than 50% of the original path segment and the wind speed direction difference angle is less than 30 degrees to generate an optimized path node combination; A trend interpolation module, which is used to input the optimized path node combination into the cubic spline interpolation method to fit the path, calculate the direction difference angle between the fitted path and the original path, and screen the paths with a difference angle less than 15 degrees to generate a main path set; A risk calibration module, which is used to calculate the slope aspect angle difference based on the slope normal angle and the main wind direction angle difference, generate a wind pressure linkage correction coefficient in combination with the wind pressure time series data, and mark the units in the main path set whose concentration change rate exceeds 40% and the correction coefficient shows a positive increase as the risk unit coordinates; A prediction update module, which is used to expand the radius area based on the risk unit coordinates, update the local prediction path using the Kalman filter algorithm, input the updated path nodes into the pollen diffusion model, and output the corrected pollen concentration prediction value and warning signal.

Images