Pollution source grid space tracing method

Abstract

The invention relates to the technical field of pollution source positioning, and particularly discloses a pollution source grid space tracing method. The method comprises the steps of obtaining three-dimensional wind speed vector data and pollutant concentration distribution data, calculating partial derivatives, constructing tensors and extracting a main strain direction, determining grid parameters based on a characteristic value ratio, aligning a grid long axis to adjust the size of a unit, generating a dynamic deformation grid layout through interpolation, calculating a vector included angle, and clustering and marking a wind direction sudden change region boundary. Boundary points are extracted for spline fitting, a grid main shaft is reconstructed to construct a transmission path grid, gradient is calculated based on concentration distribution, a transmission path is tracked, a concentration peak value area is identified, and a pollution source space positioning result is obtained. According to the method, the wind speed change and pollutant concentration information are fused, dynamic grid adjustment and wind direction mutation clustering are combined, the transmission path tracking precision and the spatial resolution are improved, grid division errors are overcome, stable recognition of pollutant migration rules is enhanced, and the accuracy and reliability of pollution source positioning are improved.

Description

Technical Field

[0001] This invention relates to the field of pollution source location technology, and in particular to a method for tracing pollution sources using a grid spatial approach. Background Technology

[0002] The field of pollution source location technology mainly studies methods for identifying and tracing pollution sources in environmental monitoring and pollution control. It obtains the concentration information of pollutants in environmental media through means such as air quality monitoring, water pollution detection, and soil sample analysis, and analyzes the distribution and diffusion of pollutants in combination with geospatial information to determine the spatial location and emission characteristics of pollution sources. This field covers pollutant monitoring data collection, spatial distribution modeling, diffusion path inference, and pollution source area location, and is widely used in pollution control of atmospheric, water, and soil environments.

[0003] Traditional methods for spatial source tracing of pollution sources involve collecting and comparing data from monitoring points within a defined monitoring area using a fixed grid. The location of the pollution source is determined by the spatial differences in pollutant concentration measurements within each grid cell. This is achieved through methods such as monitoring point deployment and sampling, grid-by-grid comparison of concentration values, and spatial overlay analysis based on geographic coordinates. In addition, there are pollution source inversion methods based on reverse trajectory analysis or adjoint models, which simulate the reverse diffusion path of pollutants using meteorological fields to estimate the source area location. There are also empirical source tracing methods based on statistical correlation or machine learning models, which utilize historical monitoring data and meteorological characteristics to establish pollution source prediction models. However, traditional methods generally suffer from high computational demands, strong dependence on meteorological input, or inability to dynamically respond to wind field changes. Especially under complex wind speed gradients or significant local disturbances, traditional inversion or empirical models struggle to accurately depict the actual transport path of pollutants, resulting in insufficient accuracy in pollution source location.

[0004] Fixed-grid spatial source tracing methods cannot effectively adjust grid division when wind fields undergo abrupt changes or gradient shifts, resulting in an inaccurate reflection of wind field changes. This limitation is particularly pronounced when wind speed direction changes significantly or terrain is complex, especially in complex environments such as urban areas with multiple wind corridors or mountainous regions. Fixed grids cannot adapt to airflow transitions and local disturbances, causing deviations in pollutant diffusion path calculations, affecting the accuracy of pollution source location, and increasing the uncertainty of environmental governance decisions. Existing technologies have limited research on grid adaptive mechanisms for the nonlinear spatial variation characteristics of wind fields, often employing static or semi-dynamic grid division methods that only update concentration distribution in the time dimension, lacking real-time responses to changes in wind speed tensor characteristics. Especially during pollutant source tracing, the inability to directionally adjust the grid structure in conjunction with the principal strain direction of the wind field leads to insufficient source tracing accuracy and poor path continuity. Furthermore, traditional methods often use empirical partitioning or fixed threshold division in areas of abrupt wind direction changes, ignoring the spatial continuity of the wind field structure and failing to accurately identify abrupt boundary changes. Summary of the Invention

[0005] To address the issue that fixed-grid spatial source tracing methods cannot effectively adjust grid division when wind fields undergo abrupt changes or gradient shifts, resulting in inaccurate reflection of wind field changes, this limitation is particularly pronounced when wind speed direction changes significantly or terrain is complex. Especially in complex environments such as urban areas with multiple wind corridors or mountainous regions, fixed grids cannot adapt to airflow transitions and local disturbances, causing deviations in pollutant diffusion path calculations, affecting the accuracy of pollution source location, and increasing the uncertainty of environmental governance decisions. This invention provides a method for spatial source tracing using a grid, comprising the following steps: To achieve the above objectives, the present invention adopts the following technical solution: a method for tracing pollution sources using a grid spatial approach, comprising the following steps: S1: Obtain three-dimensional wind speed vector data and pollutant concentration distribution data, calculate the partial derivative of wind speed components using the finite difference method, construct the wind speed change rate tensor, perform eigenvalue decomposition on the tensor to extract the principal strain direction vector, and calculate the ratio of principal and secondary eigenvalues ​​to obtain the wind field adaptive grid parameters. S2: Based on the wind field adaptive grid parameters, align the grid along the long axis, adjust the element size, calculate the angle difference of the principal strain direction vector between adjacent elements, and when the angle exceeds the threshold, use linear interpolation to perform principal axis direction interpolation to generate a dynamic deformation grid layout. S3: Call the principal strain direction vector in the dynamic deformation mesh layout, calculate the angle between adjacent vectors point by point, mark the wind direction change point when it exceeds the change threshold, and use the K-means clustering algorithm to spatially group the change points to form the boundary of the wind direction change area. S4: Extract the coordinates of the boundary points on both sides of the boundary of the wind direction change area, use the cubic spline interpolation function to fit the boundary curve of the buffer zone, reconstruct the grid principal axis direction in the dynamic deformation grid layout, and construct a pollutant transport path tracking grid.

[0006] As a further aspect of the present invention, the wind field adaptive grid parameters include principal axis direction, cell size ratio and directional consistency coefficient; the dynamic deformation grid layout includes grid orientation, local continuity and spatial resolution; the boundary of the wind direction change region includes boundary range, boundary shape and spatial partitioning; and the pollutant transport path tracking grid includes path extensibility, boundary constraint and directional continuity.

[0007] As a further aspect of the present invention, the specific steps of S1 are as follows: S101: Obtain three-dimensional wind speed vector data and pollutant concentration distribution data. Based on the wind speed component values ​​adjacent to the grid nodes, use the finite difference method to calculate the partial derivative of the wind speed component, construct the rate of change matrix at the grid node, and extract the partial derivative value of the wind speed component based on the matrix to obtain the partial derivative value of the wind speed component. S102: Call the partial derivative values ​​of the wind speed components, combine them to construct the wind speed change rate tensor, and perform eigenvalue decomposition on the grid nodes based on the tensor to extract the principal strain direction vectors. Analyze the strain direction information of each grid node through matrix operations to generate a principal strain direction vector group. S103: Based on the main strain direction vector group, calculate the ratio of the main eigenvalue to the secondary eigenvalue. For each grid node, perform a weighted comparison based on its eigenvalue ratio and pollutant concentration distribution data to generate wind field adaptive grid parameters.

[0008] As a further aspect of the present invention, the specific steps of S2 are as follows: S201: Based on the wind field adaptive grid parameters, align the grid major axis direction, obtain the boundary vector of the grid cell, compare it with the main wind field direction, adjust the angle between the grid cell major axis and the main wind field direction, and correct the aspect ratio of the grid cell to obtain the cell orientation angle value. The unit orientation angle value is the angle between the long axis direction of the grid unit and the main direction of the wind field. S202: Call the unit orientation angle value, calculate the angle difference between the principal strain direction vectors of adjacent units, determine the angle difference according to the set angle threshold, record the units whose angle difference exceeds the threshold, and extract their angle difference value to obtain the angle difference between adjacent units. The angle threshold is the critical value for determining the difference in the main direction angle between adjacent units, and is set to 5°~15°. It is adaptively adjusted according to the scale of wind field changes and the requirements for simulation accuracy. S203: Based on the angle difference between adjacent units, linear interpolation is used to perform directional interpolation in the main axis direction, and the unit arrangement order is adjusted in combination with the coordinate position of the grid units to obtain a dynamic deformation grid layout.

[0009] As a further aspect of the present invention, the specific steps of S3 are as follows: S301: Call the principal strain direction vector in the dynamic deformation mesh layout, calculate the angle between adjacent position vectors point by point, and compare the calculated angle with the set mutation threshold. When the angle exceeds the threshold, it is marked as a mutation point, and an angle mutation mark value is generated. S302: Based on the angle mutation marker value, extract the coordinate information of the marked mutation points, use the K-means clustering algorithm, set the initial cluster center, perform iterative calculation based on the coordinate information of the mutation points in space, adjust the category division according to the deviation from the cluster center, and generate the mutation point clustering coefficient. S303: Based on the clustering coefficient of the mutation points, the coordinates of mutation points in the same category are subjected to boundary fitting calculation through a spatial fitting algorithm to divide the spatial boundary of the mutation points and form the boundary of the wind direction mutation area.

[0010] As a further aspect of the present invention, the angle mutation marker value refers to the marker data performed at the corresponding position when the included angle between adjacent position vectors exceeds a preset mutation threshold; The mutation point clustering coefficient refers to the data generated after classifying mutation points using a clustering algorithm; The mutation threshold is a preset value used to determine whether the angle between adjacent position vectors constitutes a mutation.

[0011] As a further aspect of the present invention, the specific steps of S4 are as follows: S401: Call the boundary of the wind direction change area to extract the coordinates of the boundary points on both sides, detect the extracted boundary point coordinates in sequence, calculate the spatial distance between adjacent points and sum them to obtain the average value, and generate the boundary point spacing value; S402: Call the boundary point spacing value, use the cubic spline interpolation function to fit the coordinates of the boundary points on both sides, calculate the tangential angle of the interpolation node and deduce the curvature, and generate the buffer zone curvature coefficient; S403: Based on the curvature coefficient of the buffer zone, reconstruct the main axis direction of the grid in the dynamic deformation grid layout, register with the grid cell boundary according to the geometric distribution characteristics of the fitted curve, calculate the offset angle value of the main axis direction relative to the curve direction within the cell, and establish a pollutant transport path tracking grid.

[0012] As a further aspect of the present invention, the boundary point spacing value refers to the value obtained by calculating the spatial distance between adjacent boundary points on the boundary of the wind direction change area and taking the average value. The curvature coefficient of the buffer zone refers to the curvature value calculated based on the change of tangential angle after fitting the boundary point coordinates through a cubic spline interpolation function.

[0013] As a further aspect of the present invention, the method further includes step S5: S5: Call the pollutant transport path tracking grid to allocate pollutant concentration values ​​to the unit center point, calculate the unit concentration gradient value along the principal strain direction vector, determine the pollutant transport path through gradient reverse tracking, identify the concentration peak point at the beginning of the path, and obtain the spatial positioning result of the pollution source. The spatial location results of the pollution sources include peak coordinates, spatial location, and source point distribution.

[0014] As a further aspect of the present invention, the specific steps of S5 are as follows: S501: Based on the pollutant transport path tracking grid, the pollutant concentration value is assigned to the center point of the unit. The concentration difference between multiple units is analyzed in combination with spatial coordinate parameters. The concentration distribution of the unit in the grid area is calculated based on the concentration difference between units, and the unit concentration distribution value is generated. S502: Call the unit concentration distribution value, and based on the spatial distance between the principal strain direction vector and the coordinates of adjacent units, use a differential method to compare the angle relationship between the concentration difference of adjacent units and the direction vector, calculate the concentration gradient of multiple units, and obtain the unit concentration gradient sequence. S503: Based on the unit concentration gradient sequence, call the concentration gradient values ​​of adjacent units to perform gradient reverse tracing, determine the transmission direction of the path nodes in spatial coordinates step by step, lock the starting end of the path and locate the coordinates of the peak point, and obtain the spatial positioning result of the pollution source.

[0015] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, by combining wind speed vector data and pollutant concentration distribution data to calculate the rate of change of wind speed and extracting the principal strain direction, dynamic adjustments can be made during grid construction, ensuring that the grid layout is optimized according to wind field changes, no longer relying on fixed grid division. Identification and clustering of abrupt wind direction changes help clarify the boundaries of unstable wind field regions. Furthermore, by fitting curves and reconstructing the grid principal axes, the tracking accuracy of pollutant transport paths is improved. This effectively avoids the errors caused by uniform division in traditional methods, improves spatial resolution, ensures stable tracking of pollutant migration patterns under complex wind field conditions, and enhances the accuracy and reliability of pollution source location. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 This is a schematic diagram of the steps of the present invention; Figure 2 This is a detailed schematic diagram of S1 of the present invention; Figure 3 This is a detailed schematic diagram of S2 of the present invention; Figure 4 This is a detailed schematic diagram of S3 of the present invention; Figure 5 This is a detailed schematic diagram of S4 of the present invention; Figure 6 This is a detailed schematic diagram of S5 of the present invention. Detailed Implementation

[0018] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0019] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0020] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the difference between them, they convey the same meaning. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the difference between them, they convey the same meaning.

[0021] In this embodiment of the invention, sometimes a subscript such as W1 may be written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0022] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0023] Please see Figure 1This invention provides a method for spatial tracing of pollution sources using a grid, comprising the following steps: S1: Obtain three-dimensional wind speed vector data and pollutant concentration distribution data, calculate the partial derivative of wind speed components using the finite difference method, construct the wind speed change rate tensor, perform eigenvalue decomposition on the tensor to extract the principal strain direction vector, and calculate the ratio of principal and secondary eigenvalues ​​to obtain the wind field adaptive grid parameters. S2: Align the long axis of the grid based on the wind field adaptive grid parameters, adjust the element size, calculate the angle difference of the principal strain direction vector of adjacent elements, and when the angle exceeds the threshold, use linear interpolation to perform principal axis direction interpolation to generate a dynamic deformation grid layout. S3: Call the principal strain direction vector in the dynamic deformation mesh layout, calculate the angle between adjacent vectors point by point, mark the wind direction change point when it exceeds the change threshold, and use the K-means clustering algorithm to spatially group the change points to form the boundary of the wind direction change area. S4: Extract the coordinates of the boundary points on both sides of the wind direction change zone boundary, use the cubic spline interpolation function to fit the boundary curve of the buffer zone, reconstruct the grid principal axis direction in the dynamic deformation grid layout, and construct a pollutant transport path tracking grid. S5: Call the pollutant transport path tracking grid to allocate pollutant concentration values ​​to the cell center point, calculate the cell concentration gradient value along the principal strain direction vector, determine the pollutant transport path through gradient back tracking, identify the concentration peak point at the beginning of the path, and obtain the spatial location result of the pollution source; The wind field adaptive grid parameters include principal axis direction, cell size ratio, and directional consistency coefficient; the dynamic deformation grid layout includes grid orientation, local continuity, and spatial resolution; the boundary of the wind direction change zone includes boundary range, boundary shape, and spatial partitioning; the pollutant transport path tracking grid includes path extension, boundary constraint, and directional continuity; and the spatial location results of pollution sources include peak coordinates, spatial location, and source point distribution.

[0024] Please see Figure 2 The specific steps of S1 are as follows: S101: Obtain three-dimensional wind speed vector data and pollutant concentration distribution data. Based on the wind speed component values ​​adjacent to the grid nodes, use the finite difference method to calculate the partial derivative of the wind speed component, construct the rate of change matrix at the grid node, and extract the partial derivative value of the wind speed component based on the matrix to obtain the partial derivative value of the wind speed component. Three-dimensional wind speed vector data and pollutant concentration distribution data are acquired through atmospheric monitoring equipment. Instantaneous wind speed measurements in three orthogonal directions are retrieved, and the wind speed components in the x, y, and z directions of each grid node are read. These components are then categorized according to the spatial grid number. For example, in a 10×10×5 three-dimensional cubic grid, the position corresponding to the number (i, j, k) has a wind speed component. , , As the initial input data, pollutant concentration data are then used at the same spatial coordinate points. Data is collected in the form of ppm measurements from sensors, which are directly correlated with the corresponding grid point numbers. For calculating the rate of change of wind speed components, wind speed values ​​from adjacent grid points are selected for each grid point, and spatial finite difference calculations are performed, for example... Quantity The rate of change of direction is taken as the ratio of the adjacent right node to the left node. Value difference, divided by the grid distance between the two points Similarly, along , Directions respectively , , Perform a similar operation on the components to form partial derivatives. , , Equal numerical values; in specific execution, for example, the node numbered (5, 5, 3) has its adjacent right side... m / s, adjacent left side m / s, grid spacing m, then The rate of change of direction s-1, similarly in If the direction m / s, m / s, grid spacing m, then the rate of change s -1 ,exist Direction such as m / s, m / s, spacing m, then the rate of change is s -1 The summation of the three directional rates of change for each of the above components constitutes the wind speed gradient matrix for that node. The structure is as follows: ; At node (5, 5, 3), assume that the distance is measured through adjacent nodes. The rate of change of component direction , , s -1 , The rate of change of component direction , , s -1 Then the matrix is: ; The above values ​​are the partial derivatives of the wind speed components. Components with larger drift values ​​indicate significant changes in wind speed in that direction. For actual monitoring scenarios, the numerical range can be divided using threshold grading methods, such as the absolute value of the rate of change. s -1 This is denoted as the lower interval. absolute value s -1 Let be the middle interval, and be the absolute value. s -1 Let this be the high interval; such as the partial derivative at this node. s -1 That is, falling into the high range indicates The directional airflow gradient changes dramatically and is closely related to the spatial transport of pollutants.

[0025] S102: Call the partial derivatives of the wind speed components, combine them to construct the wind speed change rate tensor, and perform eigenvalue decomposition on the grid nodes based on the tensor to extract the principal strain direction vectors. Analyze the strain direction information of each grid node through matrix operations to generate a principal strain direction vector group. When calling the partial derivatives of the wind speed components to construct the wind speed rate of change tensor, this matrix is ​​treated as a numerical representation of the velocity gradient field, and is filled one by one according to the node coordinates. Then, it is symmetricized to extract the strain rate tensor. That is, take and Average value: ; When substituting into the above formula, first find the transpose matrix. : ; By adding the matrices term by term, we get: ; Then perform the matrix division by 2 operation, and you get ; Next Eigenvalue decomposition yields three eigenvalues. , , And the corresponding eigenvectors. During eigenvalue decomposition, the determinant is calculated. Solve by algebraic operations For example, in this case, the numerical solution yielded: ; Pick The corresponding eigenvectors serve as the principal strain direction vectors. Solve the system of linear equations The normalized eigenvectors are obtained as follows: ; This vector represents the proportion of the principal strain direction in the x, y, and z components within a unit length of space, indicating the main stretching direction of the wind field at that node. It can be directly used for subsequent calculations of the spatial distribution relationship with pollutant concentrations. In practical applications, to avoid misclassifying small eigenvalues ​​close to zero as principal eigenvalues, an eigenvalue filtering threshold is set before solving the problem. ,For example s -1 (Referring to the historical statistical characteristic value scale of the wind field monitoring area of ​​this type, if the value is less than this, the characteristic significance is not significant.) In this example, all three characteristic values ​​are greater than this threshold, so they are all included in the primary and secondary judgment.

[0026] Repeat the above operations on multiple nodes to differentiate the nodes. The principal strain direction vectors are summarized and formed, and the data structure of this vector group can be defined as an array: ; in Corresponding to 3D grid coordinates, for example Nodes .

[0027] Table 1: Partial Examples of Principal Strain Direction Vector Groups (Unit: Dimensionless)

[0028] As shown in Table 2, the vector is composed of... The matrix eigenvalue decomposition yields vectors whose directions define the main strain orientation of the local wind field at each node, which are ultimately integrated into a main strain direction vector group.

[0029] S103: Based on the principal strain direction vector group, calculate the ratio of principal eigenvalues ​​to secondary eigenvalues. For each grid node, perform a weighted comparison with the pollutant concentration distribution data based on its eigenvalue ratio to generate wind field adaptive grid parameters. When performing the eigenvalue ratio calculation based on the principal strain direction vector group, for node (5, 5, 3), three sets of eigenvalues ​​are called. , , According to the formula for the ratio of principal eigenvalues ​​to secondary eigenvalues Divide the two values ​​directly to get This ratio indicates that the elongation rate of the node in the principal strain direction is 5 times that in the secondary strain direction. When comparing multiple nodes, it is necessary to consider the node's elongation rate in the principal strain direction. Values ​​form a sequence And set the ratio grading intervals, for example Classified into the low ratio range Included in the middle range. In this example, the high ratio interval is included. It belongs to the high ratio interval; then for each node, it is... The values ​​are compared with the pollutant concentration data using a weighted average. In this example, the pollutant concentration at node (5, 5, 3) is... The sampled value was 120 ppm, and the weighted comparison calculation method was as follows: ,in This is the concentration normalized value. For ratio weights, For concentration weighting, select , (Weights were set based on the importance ratios and concentrations in historical monitoring data). By regional concentration range ppm linear normalization, for example Substitute the values ​​into the weighted calculation formula: ; The weighted comparison value of this node is 3.04. For comparability analysis, [the following value is used]. Values ​​and preset wind field adaptability parameter thresholds Make comparisons, assume (This threshold is derived from the original data of wind field simulation and source tracing optimization, and is taken as the boundary value that can effectively distinguish between obvious adaptability and maladaptability.) Then this node... The node was determined to have significant wind field adaptability. Following this process, the ratio calculation and weighting were repeated for each grid node, and dense sampling was performed near nodes with extremely high ratios to supplement data and improve data accuracy. The node's... The values ​​form a set of overall wind field adaptability parameters. This parameter set corresponds to the three-dimensional distribution of the grid in space and can be directly used as input data for subsequent spatial pollution source tracing. In implementation, it is also necessary to ensure that the node concentration normalization range is calculated based on the global minimum and maximum values ​​of the same batch of monitoring datasets, avoiding inconsistent normalization results introduced by measurements in different time periods. For example, in another node (5, 6, 3), if... , but It falls within the medium-range, and if the concentration at this node... ppm Weighted calculation The value is below the threshold. If the wind field adaptability is not significant, then the nodes are compared one by one to complete the entire process of generating the wind field adaptability grid parameters from the principal strain direction vector group.

[0030] Please see Figure 3The specific steps of S2 are as follows: S201: Align the long axis of the grid with the wind field adaptive grid parameters, obtain the boundary vector of the grid cell, compare it with the main wind field direction, adjust the angle between the long axis of the grid cell and the main wind field direction, and correct the aspect ratio of the grid cell to obtain the cell orientation angle value. The element orientation angle value is the angle between the long axis direction of the grid element and the main direction of the wind field. Based on wind field adaptive grid parameters First, the three-dimensional mesh array cells are aligned along their major axis, that is, aligned in each mesh cell with... Based on the principal direction information, the principal strain direction vector is extracted. And as a reference direction for major axis optimization, such as nodes. , If the expected major axis direction of the node is consistent with this vector, then the boundary vector of the mesh cell can be obtained. , , , respectively representing the element along , , The three unit vectors of the boundary direction, in actual measurement , , The direction of the element's intrinsic major axis is defined as... (Initialization available edge) Unit vector of the axis ),Will and Calculate the included angle, and take the angle as: ; Substitute the example data: ; Vector magnitude: ; but ; Compare the calculated angle with the set principal direction angle tolerance. (Example set as) (Based on engineering experience and mesh deformation control requirements) comparison, if If so, it is determined that the major axis direction needs to be adjusted, in this example. Exceed Perform the adjustment operation; during the adjustment process, adjust the long axis of the unit. According to the rotation matrix Along the axis of rotation (Normalized as the direction of the rotation axis) Rotation angle After rotation Vector and The included angle represents the corrected alignment state; subsequently, the corrected aspect ratio of the unit is calculated, defining the original aspect ratio as... ,in and These are the actual lengths of the element along its major and minor axes, respectively. In this example, the measured lengths are... m、 m, then According to the rotation correction ratio: ; Updated aspect ratio ; Finally, after rotating Vector re-interaction By comparison, the final included angle is obtained: ; The corrected angle is This value is stored as the element orientation angle value of this node, making the three data points: corrected aspect ratio. Final angle and corresponding coordinates A set of orientation parameter data is generated, which is the result of the element orientation angle value of the node.

[0031] S202: Call the element orientation angle value, calculate the angle difference between the principal strain direction vectors of adjacent elements, judge the angle difference according to the set angle threshold, record the elements whose angle difference exceeds the threshold, and extract their angle difference value to obtain the angle difference between adjacent elements. The angle threshold is the critical value for judging the difference in the main direction angle between adjacent units, set to 5°~15°, and adaptively adjusted according to the scale of wind field changes and simulation accuracy requirements; Calling the unit orientation angle value sequence First, for each cell in the 3D mesh, locate its neighboring cells, for example, for nodes. Its neighboring nodes are The direction is acceptable and ,exist The direction is acceptable and ,exist The direction is acceptable and ,by Taking two adjacent nodes in a direction as an example, the angle difference calculation operation is performed. Nodes , Nodes Then the angle difference Similarly, in direction, node ,but ;exist direction, Nodes ,but ; Calculate the angle difference in each direction and the angle threshold. Comparison, Within the set range to The value is adaptively adjusted based on the scale of wind field changes and the required simulation accuracy. For example, when the grid scale in this area is small and more local variations need to be preserved, the value is set to... lower limit Select As the judgment value; when performing the angle difference judgment action, it is determined whether the threshold is exceeded for each direction: direction No record, direction Record the pair of units. direction Record the pair of elements; for the recorded adjacent elements, their specific angle differences also need to be extracted and stored, for example... Directional difference , Directional difference In actual computation, for each node This adjacent comparison and threshold determination process needs to be repeated, and the results need to be used to generate a difference set grouped by direction. and decision mark set ,in Indicates direction Exceeding the threshold Indicates direction The threshold has not been exceeded, for example, the node The result is the difference set (Unit: °) and Marker Set To reduce discontinuities in subsequent directional interpolation, a range check will be performed on the differences in the same direction during result processing. The interval is considered as the range of adaptation changes, and will be greater than The entries marked as range values ​​need to be processed separately in subsequent interpolation. For example, if the difference between adjacent directions of a node is... Then at this node The marker value is recorded as 2 separately for storage as a special distinction type; the above process is performed one by one on the grid cells to obtain the angle difference data between adjacent cells.

[0032] S203: Based on the angle difference between adjacent elements, linear interpolation is used to perform directional interpolation in the principal axis direction, and the element arrangement order is adjusted in combination with the coordinate position of the grid elements to obtain a dynamic deformation grid layout. The set of angle differences between adjacent units and the set of judgment markers are called. First, for angle differences that are not less than a set threshold, the following criteria are used: And not greater than the maximum fit difference The node pairs perform directional interpolation, for example, the node exist Directional difference ,exist Directional difference All meet the interpolation conditions; in linear interpolation operations, consider two units that need to be interpolated. and The angles of the principal strain directions in a certain direction are respectively and The spatial distance between the two is calculated based on the difference in node coordinates and the grid step size, for example... Direction and Node distance is m, then in this direction Interpolation angle at each equal division point It can be calculated using the formula: ; in Number the interpolation points; of and of For example, if (Two interpolation points), then Interpolation angle: ; hour ; After the interpolation results are determined, the principal axis direction of each interpolation node is updated to the omnidirectional vector corresponding to the interpolation angle (obtained by rotating the initial direction of the major axis around that angle), and the aspect ratio of the corresponding node is updated. The value was adjusted to match the new orientation angle, and the calculation method remained the same. ,in The original rotation angle is replaced with the interpolated angle, and the calculation is recalculated to ensure that the geometry and orientation are consistent. This process requires consideration of the node's three-dimensional coordinate position. Adjust the unit arrangement order and renumber the nodes along the interpolation direction according to the rule of increasing or decreasing the interpolation result, for example, in Direction, from In the direction of interpolation, the node numbers should be arranged from smallest to largest; otherwise, the geometric topology in the interpolation region will be destroyed. (See this example.) The result of the three nodes in the direction is When, the sorting is interpolation nodes This permutation is recorded in the data structure; Direction, similarly of and of Intersect three evenly spaced points at 20m intervals, and calculate the angle difference. Increment at each 5m step The interpolation result is obtained. Fill in the corresponding intermediate nodes in sequence; Table 2: Example of interpolation calculation (unit: °)

[0033] As shown in Table 2, the interpolation calculation is strictly performed on node pairs that exceed the threshold but are within the preset upper limit. The angle result obtained by interpolation is bound to the spatial node position and used for mesh reordering. After completing the interpolation and arrangement adjustment in all directions, the updated node major axis direction and aspect ratio are used to form a new three-dimensional array structure, and the new direction vector, aspect ratio and coordinate position of each node are recorded in the mesh data file to obtain the complete dataset of the dynamic deformation mesh layout.

[0034] Please see Figure 4 The specific steps of S3 are as follows: S301: Call the principal strain direction vector in the dynamic deformation mesh layout, calculate the angle between adjacent vectors point by point, and compare the calculated angle with the set abrupt change threshold. When the angle exceeds the threshold, it is marked as an abrupt change point, and an angle abrupt change marker value is generated. To retrieve the principal strain direction vector of each node in the dynamic deformation mesh layout, first, for the nodes... Read its and adjacent nodes of When calculating the angle between the two, the vector dot product operation is performed first: ; in express Components, in this example ; Next, calculate the vector magnitude: ; ; Substitute the value into the formula for the included angle: ; Denominator: ; ratio: ; ; Similarly, calculate and The angle between the nodal direction vectors and The angle between the node direction vectors is calculated, and the above calculation steps are repeated to obtain the angle values ​​between adjacent elements in multiple directions. After completing the point-by-point calculation of adjacent node pairs, the results are compared with the mutation threshold. Compare, After selecting engineering parameters within the range The adjustment is adaptive, and in this example, it is based on the stability requirement of the wind field's main direction distribution. If the angle change is less than this value, it is considered to be in a gradual range; if it is greater than this value, it is considered to be a sudden change. Continue to use the aforementioned... and The calculation example shows that the included angle is... If it is determined to be a mutation point, then the record table between that node will be... The access direction is marked as 1 (mutation), otherwise if it is less than or equal to The node is then marked as 0 (gradual). After performing this judgment and marking on all node pairs, the mutation judgment results of the nodes are stored in the marking matrix according to the node number and direction component. ,in Value This matrix represents the mutation state relative to the node in that direction. The elements can take two values ​​{0, 1}, where 1 represents a mutation and 0 represents no mutation. In this example, for... Node, assuming in Direction angle (Mark 0) direction (Mark 1) direction (Mark 1), then the mutation marker value of this node and the corresponding three directions is To facilitate subsequent clustering processing, the positions of nodes marked as 1 need to be... With direction Export the coordinates of the mutation points together. Each record contains both the three-dimensional spatial location and the direction information of the mutation. Once the marking is completed, the set of mutation point coordinates corresponds to the set of angle mutation mark values.

[0035] S302: Based on the angle mutation marker value, extract the coordinate information of the marked mutation points, use the K-means clustering algorithm, set the initial cluster center, perform iterative calculation based on the coordinate information of the mutation points in space, adjust the category division according to the deviation from the cluster center, and generate the mutation point clustering coefficient; To access the set of coordinates for mutation points, first convert each data record in the set into actual spatial coordinates. This coordinate is obtained by converting the node number to the grid spacing, for example, node Grid spacing at The three directions are respectively m、 m、 m, coordinates of the origin ,but m, m, m, the spatial coordinates, are used as input variables for subsequent clustering calculations; when performing K-means clustering, the initial number of cluster centers must be set first. The principle for selecting this value is to ensure that the number of mutation points within each class is not less than the minimum number of cluster points. In this example, a total of 48 mutation points were measured, based on experience. ,but: ; Final decision As the initial number of cluster centers; then... The points in the middle are randomly selected in order. Use different locations as initial centers For example, points respectively , , , , Calculate the relationship between each mutation point and this... Euclidean distance between the centers: ; And assign the point to the nearest central category, mutation point. right The distance is: ; right distance m, the distance between the other centers is larger, therefore Classified The class to which the point belongs; after completing the initial class division of the points, calculate the new center coordinates for each class, the new center The calculation is the average of the coordinates of multiple types of points, that is: , , ; in To determine the number of points of this type, this example... The new center is obtained by averaging the 11 points of the class. The new center is used as the cluster center for the next iteration. The process of distance calculation, category adjustment, and center update is repeated until the center coordinates change between two adjacent iterations. ; Convergence threshold In this example, the resolution is set to 0.5m (this value is 2.5% of the minimum grid resolution of 20m). The clustering converges after 7 iterations, yielding the final clustering result. For each mutation point, its corresponding category number is recorded. and its original coordinates and direction Store them together and calculate the average distance between each point within a class and the center: ; As the clustering coefficient of the class, the clustering coefficient range in this example is [15.2m, 42.8m]. A smaller coefficient indicates that the distribution of points within the class is more concentrated. The final output "mutation point clustering coefficient" dataset includes the class number, center coordinates, clustering coefficient and set of mutation points within the class.

[0036] S303: Based on the clustering coefficient of the mutation points, the coordinates of the mutation points in the same category are used to perform boundary fitting calculations through a spatial fitting algorithm to divide the spatial boundaries of the mutation points and form the boundaries of the wind direction mutation area. Based on the set of mutation points and the corresponding clustering coefficients, the category numbers are first selected. The class is used as a calculation example, and the coordinates of the cluster center of this class are... m, number of points within the class The three-dimensional coordinates of a point are determined by clustering the coordinates of mutation points. For example, some points in this class are: , , Before performing spatial boundary fitting, import this type of point set into a three-dimensional coordinate array. And then Coordinate system normalization is performed, which involves subtracting the minimum value of each dimension from the coordinates of that dimension and dividing by the range of that dimension. This ensures consistent orientation and scale during the fitting process. In this example... The minimum value of the coordinates is 198m, and the maximum value is 224m. Then, for a certain point: ; right , Similarly, coordinates are normalized; then, a parametric method for spatial fitting is selected for boundary description, and the boundary function is set. Two-dimensional parameters Controlling 3D Surfaces The coordinates of the surface control points are determined using the least squares method to minimize the sum of the squared average distances between the surface and the set of abrupt change points; in actual calculations, the categories are considered. The radial vector from each point to the initial centroid of the surface is normalized and then multiplied by the cluster radius. m is used as the initial radius for fitting, and the radial component of the point is expressed using spherical coordinates. It means that, among them Radial length, It is the azimuth angle. The polar angle is the midpoint in this example. The displacement vector relative to the center is: m; Radial length: m; Azimuth ; Polar angle ; After converting all points to the spherical parameter space, press Divide the grid according to the rules and fit it. Surface function, averaging over each parameter grid point The values ​​are then returned to the Cartesian coordinate system to generate the boundary surface coordinate set. ,Right now: ; After the calculation is complete, the boundary of this class is represented by a set of points on a three-dimensional closed surface; this operation is performed sequentially on... Performed for each cluster category, in this example Each category corresponds to The calculated boundary closed surface datasets are denoted as follows: 15.2m, 20.5m, 33.1m, 27.8m, and 42.8m respectively. To allow for direct use in subsequent wind field visualization, the boundary point set needs to be divided into individual... The value of the cross section projected onto Plane, generate contour line description, to Output the cross-sectional curve equation as a variable, at height intervals. m generates layered boundaries, thus forming a spatial layered boundary set for the wind direction change zone. This data not only retains the three-dimensional surface coordinates, but also has a layered two-dimensional outer contour representation, which serves as the final output result of the wind direction change zone boundary.

[0037] Please see Figure 5 The specific steps of S4 are as follows: S401: Call the boundary of the wind direction change area to extract the coordinates of the boundary points on both sides, detect the extracted boundary point coordinates in sequence, calculate the spatial distance between adjacent points and sum them to obtain the average value, and generate the boundary point spacing value; To access the boundary of a wind direction change zone, first select the set of boundary points on both sides of a certain change zone. (Left boundary) and (Right boundary), assuming in At section m It contains 8 boundary points, whose coordinates are as follows: ...; It contains 8 points, and their coordinates are as follows: , ...; For points on a boundary line, it is necessary to ensure that it is a monotonically continuous boundary sequence. Therefore, we first base it on... or Sort by the trend of coordinate changes, so that adjacent points The physical adjacency relationship is consistent with the sequential adjacency, and then the spatial distance value of adjacent points is calculated one by one. The calculation formula is: ; For example in Inside, and The distance is: m, m, m; ; For the whole All on the boundary Performing the same operation on each pair of adjacent points yields a distance sequence. ;exist The distance between adjacent points is also calculated on the boundary, resulting in... Average the two sets of distance values ​​along the boundaries, and define the average spacing on the left: ; The value on the right Similarly, in this example, we can calculate... m、 m; to obtain the comprehensive boundary point spacing value ,use: ; Substitute the values This value serves as the input scale reference for subsequent calculations of the buffer zone curvature.

[0038] Table 3: Example of boundary point spacing calculation (unit: m)

[0039] As shown in Table 3, the spatial distance between adjacent points is obtained by summing and taking the square root of the squared differences of each coordinate. Then, the average value is used to calculate the comprehensive distance between boundary points. (Example) m is the result of the boundary point spacing value in this step.

[0040] S402: Call the boundary point spacing value, use the cubic spline interpolation function to fit the coordinates of the boundary points on both sides, calculate the tangential angle of the interpolation node and deduce the curvature, and generate the buffer zone curvature coefficient; To retrieve the boundary point spacing value, first, set the coordinates of the boundary points on both the left and right sides of the area where the wind direction changes abruptly. and In the process, a sequence of nodes for fitting is generated for each boundary line. Each of them All according to Indicating spatial location, the sequence order remains unchanged according to the order of adjacent points in S401; assuming it is on the left boundary. The CCP Each node, according to the node After sorting the coordinates in ascending order, record the spatial distance between each pair of adjacent nodes. Confirm the length of multiple segments and The proportional deviation does not exceed In this example, the deviation range obtained in the distance sequence is within If all parameters meet this standard, it means that spline interpolation fitting can be performed directly. During the cubic spline interpolation calculation, the arc length parameter of each node is first... From the cumulative distance formula Get, for example, the first node arc length m, 2nd node m, the 3rd node m; then targeting , , Fit cubic spline functions to each of the three dimensions. , , The interpolation coefficients within each spline segment are obtained by using the coordinates of the interval endpoints and the continuity conditions of the first and second derivatives, ensuring that the curve completely coincides with the boundary points at the arc length nodes; when calculating the tangential angle, the tangential vector of the interpolation node is first calculated: ; Then, the tangent vector is converted into azimuth and inclination angles using the arctangent and inverse cosine functions. In this example... The tangential vector at point m is calculated as follows: , module length azimuth ,inclination The calculation of curvature is based on the formula for space curves: ; in Let be the derivative of the tangent vector with respect to the arc length, representing the rate of change of the curve direction with respect to the arc length. In this example... The value at m is Its mold length curvature Similarly, every [time] along the entire boundary The curvature is calculated once at the arc length points of the spacing to obtain the curvature sequence. Then, the average value is taken to obtain the curvature coefficient of the buffer zone: ; For example, obtained at the left boundary. The right boundary is obtained The curvature coefficient of the overall buffer zone is then taken as the average of both sides. This result is the curvature coefficient of the buffer zone.

[0041] S403: Based on the curvature coefficient of the buffer zone, reconstruct the main axis direction of the grid in the dynamic deformation grid layout, register it with the grid cell boundary according to the geometric distribution characteristics of the fitted curve, calculate the offset angle value of the main axis direction relative to the curve direction within the cell, and establish a pollutant transport path tracking grid. Call the buffer band curvature coefficient First, the principal axis direction vector of each grid cell in the dynamic deformation grid layout is compared with the tangential direction of the fitted buffer zone center curve. Establish an index matching system based on spatial correspondence, which involves projecting the geometric center of each grid cell onto the nearest point on the buffer zone curve and recording the arc length parameter of the projected point. and corresponding tangential direction For a given unit, such as the center point The nearest point on the curve is m. m, the difference in Euclidean distance between the two m, the associated buffer band tangential vector in this Below is After obtaining the tangential direction, calculate the principal axis direction of the calculation unit. and The offset angle between them: ; For example, this unit The dot product of the two vectors is: ; The module lengths are respectively , Then the offset angle is: ; Since the value is slightly over 1, it is approximately equal to the numerical precision. This indicates that the main axis direction of the unit is highly aligned with the direction of the buffer band curve; if the offset angle of a unit is greater than the set alignment threshold... If so, the principal axis direction of the unit needs to be reconstructed to reduce the offset. In this example, let's assume... (Selected based on the grid's wind direction resolution requirements; values ​​less than this limit are considered aligned.) For elements with excessive offset, while maintaining the element's aspect ratio, the original principal axis... Around perpendicular to and Rotation axis To reduce the offset angle to within the threshold, while referencing Adjust the coefficient of change of rotation amplitude in space ,in This is the reference starting point value for the arc length of the projected point of this element, so that the attenuation of the rotation amplitude is more significant in regions with high curvature; after completing the orientation correction and reconstruction of all elements, the updated principal axis orientation vector set will be used. Replace the corresponding data in the original layout, and mark the path tracking attribute value for cells within the buffer zone in the grid cell properties. The remaining areas are marked The final generated global 3D mesh data contains four types of information: spatial coordinates, corrected principal axis direction, aspect ratio, and path tracing markers. This mesh is the pollutant transport path tracing mesh, which can be directly used as the data basis for subsequent flow field source tracing simulation.

[0042] Please see Figure 6 The specific steps of S5 are as follows: S501: Based on the pollutant transport path tracking grid, the pollutant concentration value is assigned to the center point of the cell. The concentration difference between multiple cells is analyzed in combination with spatial coordinate parameters. The concentration distribution of the cell within the grid area is calculated based on the concentration difference between cells, and the cell concentration distribution value is generated. Based on the pollutant transport path tracking grid, the dataset of monitored pollutant concentration values ​​is first... According to the spatial location of the measuring point With the geometric center of the grid cell The concentration is assigned based on the distance between the points. Specifically, for each unit center, the nearest measuring point is found, and its concentration value is assigned to that unit. If the center is equidistant from multiple measuring points, the arithmetic mean of the concentrations at those points is taken as the concentration value for that unit. For example, unit The center coordinates are m, and the coordinates of its three nearest monitoring points are as follows: 52.1 μg / m³ 54.3 μg / m³ and 53.4 μg / m³; The distances of the three particles to the center point of the unit are 1.56m, 1.25m, and 1.20m respectively. Therefore, the closest distance is 1.20m, corresponding to a concentration of 53.4 μg / m³. This value is taken as... After initial allocation of grid cell concentration values, the concentration differences across multiple cells are analyzed using the spatial coordinate parameters of the cells, specifically in any reference direction (here). , , In three directions, the absolute value of the concentration difference is calculated based on the distance difference between adjacent units: ; in For direction marking, for example in In direction, The concentration of m is 53.4 μg / m³, which Adjacent units in the direction If the concentration is 55.5 μg / m³, then μg / m³; in If the concentration of adjacent units is 51.0 μg / m³, then μg / m³; in If the concentration of adjacent units is 50.5 μg / m³, then μg / m³; The concentration difference between adjacent cell pairs in the entire grid is grouped and statistically analyzed according to direction to obtain a set of differences in three directions. Subsequently, by combining the concentration difference and coordinate position in each three-dimensional direction, the concentration distribution value of the element within the grid region is generated. That is, by normalizing the vector sum of the differences of a certain element in the three directions with its concentration value, the local weight coefficient of the element in the concentration distribution field is calculated. ; in This represents the maximum concentration across the entire unit; in this example, the maximum global concentration is 60.0 μg / m³. m-element, vector sum: ; but ; This indicates that the local concentration variation range accounts for approximately 7.18% of the maximum stable concentration; and The concentration distribution values ​​of this unit are recorded as follows: This calculation is performed sequentially across the entire domain and generates... Set of concentration distribution values ​​for each unit This set represents the concentration distribution data after completing this step.

[0043] S502: Call the element concentration distribution value, and based on the spatial distance between the principal strain direction vector and the coordinates of adjacent elements, use the difference method to compare the angle relationship between the concentration difference of adjacent elements and the direction vector, calculate the concentration gradient of multiple elements, and obtain the element concentration gradient sequence. To access the element concentration distribution set, first select any element from the dynamically deformed mesh. To obtain its concentration Local weights Principal strain direction vector and three-dimensional coordinates Then for each spatially adjacent cell of that cell (exist (three directions), based on coordinate differences , , Calculate the spatial distance between the two: ; The directional concentration gradient components between adjacent cells are calculated using a difference formula: Its unit is μg / m³·m⁻¹; For example, reference unit m concentration μg / m³, adjacent Directional unit m concentration μg / m³, distance m, then: μg / m³·m⁻¹; Similarly, direction μg / m³, distance m, get μg / m³·m⁻¹; in direction μg / m³, distance m, get μg / m³·m⁻¹; Organize the three directional components into a gradient vector. Then, with the principal strain direction of the element Compare the included angles: ; in ; Length of the module: ; ; ; The included angle is correlated with a preset threshold. In comparison, this example takes When the included angle is less than or equal to this value, the direction of the concentration gradient can be considered to be consistent with the direction of the principal strain; perform the above calculation on the global element to form a gradient sequence. .

[0044] S503: Based on the unit concentration gradient sequence, call the concentration gradient values ​​of adjacent units to perform gradient reverse tracing, determine the transmission direction of path nodes in spatial coordinates step by step, lock the starting end of the path and locate the coordinates of the peak point, and obtain the spatial location result of the pollution source. Based on the concentration gradient sequence, a pollutant pathway tracing chain is constructed for each unit in the reverse direction of the gradient. Specifically, in each unit... In the process of extracting gradient vectors Unitize it into Then take the opposite direction. Starting from this direction vector, the next cell is searched step by step along the grid. , The selection criterion is that its geometric center is... The direction of the line connecting the center and The included angle does not exceed the set tracking direction threshold. (This example assumes) ), and its concentration value Higher than For example, reference unit m, Among adjacent cells, select the cell whose center position vector direction is closest to the opposite direction. m, calculate the included angle of direction: ; And concentration The concentration in μg / m³ is greater than 53.4 μg / m³, therefore it is identified as the next node in the tracking chain; this process is repeated, determining the transmission direction of each path node step by step, and continuously moving towards cells with higher concentrations until no adjacent cells satisfying the conditions can be found, at which point the peak point of the path is considered reached; the coordinates of the peak point are recorded. As a suspected source of pollution, this reverse tracing was repeated for multiple starting units across the entire region. Peak points were statistically analyzed and spatially clustered according to concentration. In this example, the multiple peak points were concentrated within the coordinate range. m( direction), m( direction), m( Within the direction, the geometric center of this region is determined as m, and finally the spatial location result of the pollution source is obtained.

[0045] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for spatial source tracing of pollution sources using a grid system, characterized in that, Includes the following steps: S1: Obtain three-dimensional wind speed vector data and pollutant concentration distribution data, use the finite difference method to calculate the partial derivative of wind speed components, construct the wind speed change rate tensor, perform eigenvalue decomposition on the tensor to extract the principal strain direction vector, and calculate the ratio of principal and secondary eigenvalues ​​to obtain the wind field adaptive grid parameters. S2: Based on the wind field adaptive grid parameters, align the grid along the long axis, adjust the element size, calculate the angle difference of the principal strain direction vector between adjacent elements, and when the angle exceeds the threshold, use linear interpolation to perform principal axis direction interpolation to generate a dynamic deformation grid layout. S3: Call the principal strain direction vector in the dynamic deformation mesh layout, calculate the angle between adjacent vectors point by point, mark the wind direction change point when it exceeds the change threshold, and use the K-means clustering algorithm to spatially group the change points to form the boundary of the wind direction change area. S4: Extract the coordinates of the boundary points on both sides of the boundary of the wind direction change area, use the cubic spline interpolation function to fit the boundary curve of the buffer zone, reconstruct the grid principal axis direction in the dynamic deformation grid layout, and construct a pollutant transport path tracking grid.

2. The method for spatial tracing of pollution sources according to claim 1, characterized in that, The wind field adaptive grid parameters include principal axis direction, cell size ratio, and directional consistency coefficient; the dynamic deformation grid layout includes grid orientation, local continuity, and spatial resolution; the boundary of the wind direction change region includes boundary range, boundary shape, and spatial partitioning; and the pollutant transport path tracking grid includes path extensibility, boundary constraints, and directional continuity.

3. The method for spatial tracing of pollution sources according to claim 1, characterized in that, The specific steps of S1 are as follows: S101: Obtain three-dimensional wind speed vector data and pollutant concentration distribution data. Based on the wind speed component values ​​adjacent to the grid nodes, use the finite difference method to calculate the partial derivative of the wind speed component, construct the rate of change matrix at the grid node, and extract the partial derivative value of the wind speed component based on the matrix to obtain the partial derivative value of the wind speed component. S102: Call the partial derivative values ​​of the wind speed components, combine them to construct the wind speed change rate tensor, and perform eigenvalue decomposition on the grid nodes based on the tensor to extract the principal strain direction vectors. Analyze the strain direction information of each grid node through matrix operations to generate a principal strain direction vector group. S103: Based on the main strain direction vector group, calculate the ratio of the main eigenvalue to the secondary eigenvalue. For each grid node, perform a weighted comparison based on its eigenvalue ratio and pollutant concentration distribution data to generate wind field adaptive grid parameters.

4. The method for spatial tracing of pollution sources according to claim 3, characterized in that, The specific steps of S2 are as follows: S201: Based on the wind field adaptive grid parameters, align the grid major axis direction, obtain the boundary vector of the grid cell, compare it with the main wind field direction, adjust the angle between the grid cell major axis and the main wind field direction, and correct the aspect ratio of the grid cell to obtain the cell orientation angle value. S202: Call the unit orientation angle value, calculate the angle difference between the principal strain direction vectors of adjacent units, determine the angle difference according to the set angle threshold, record the units whose angle difference exceeds the threshold, and extract their angle difference value to obtain the angle difference value between adjacent units. S203: Based on the angle difference between adjacent units, linear interpolation is used to perform directional interpolation in the main axis direction, and the unit arrangement order is adjusted in combination with the coordinate position of the grid units to obtain a dynamic deformation grid layout.

5. The method for spatial tracing of pollution sources according to claim 4, characterized in that, The specific steps for S3 are as follows: S301: Call the principal strain direction vector in the dynamic deformation mesh layout, calculate the angle between adjacent position vectors point by point, and compare the calculated angle with the set mutation threshold. When the angle exceeds the threshold, it is marked as a mutation point, and an angle mutation mark value is generated. S302: Based on the angle mutation marker value, extract the coordinate information of the marked mutation points, use the K-means clustering algorithm, set the initial cluster center, perform iterative calculation based on the coordinate information of the mutation points in space, adjust the category division according to the deviation from the cluster center, and generate the mutation point clustering coefficient. S303: Based on the clustering coefficient of the mutation points, the coordinates of mutation points in the same category are subjected to boundary fitting calculation through a spatial fitting algorithm to divide the spatial boundary of the mutation points and form the boundary of the wind direction mutation area.

6. The method for spatial tracing of pollution sources according to claim 5, characterized in that, The angle mutation marker value refers to the marker data at the corresponding position when the included angle between adjacent vectors exceeds a preset mutation threshold; The mutation point clustering coefficient refers to the data generated after classifying mutation points using a clustering algorithm; The mutation threshold is a preset value used to determine whether the angle between adjacent position vectors constitutes a mutation.

7. The method for spatial tracing of pollution sources according to claim 5, characterized in that, The specific steps of S4 are as follows: S401: Call the boundary of the wind direction change area to extract the coordinates of the boundary points on both sides, detect the extracted boundary point coordinates in sequence, calculate the spatial distance between adjacent points and sum them to obtain the average value, and generate the boundary point spacing value; S402: Call the boundary point spacing value, use the cubic spline interpolation function to fit the coordinates of the boundary points on both sides, calculate the tangential angle of the interpolation node and deduce the curvature, and generate the buffer zone curvature coefficient; S403: Based on the curvature coefficient of the buffer zone, reconstruct the main axis direction of the grid in the dynamic deformation grid layout, register with the grid cell boundary according to the geometric distribution characteristics of the fitted curve, calculate the offset angle value of the main axis direction relative to the curve direction within the cell, and establish a pollutant transport path tracking grid.

8. The method for spatial tracing of pollution sources according to claim 7, characterized in that, The boundary point spacing value refers to the value obtained by calculating the spatial distance between adjacent boundary points on the boundary of the wind direction change area and taking the average value; The curvature coefficient of the buffer zone refers to the curvature value calculated based on the change of tangential angle after fitting the boundary point coordinates through a cubic spline interpolation function.

9. The method for spatial tracing of pollution sources according to claim 1, characterized in that, The method also includes step S5: S5: Call the pollutant transport path tracking grid to allocate pollutant concentration values ​​to the unit center point, calculate the unit concentration gradient value along the principal strain direction vector, determine the pollutant transport path through gradient reverse tracking, identify the concentration peak point at the beginning of the path, and obtain the spatial positioning result of the pollution source. The spatial location results of the pollution sources include peak coordinates, spatial location, and source point distribution.

10. The method for spatial tracing of pollution sources according to claim 9, characterized in that, The specific steps of S5 are as follows: S501: Based on the pollutant transport path tracking grid, the pollutant concentration value is assigned to the center point of the unit. The concentration difference between multiple units is analyzed in combination with spatial coordinate parameters. The concentration distribution of the unit in the grid area is calculated based on the concentration difference between units, and the unit concentration distribution value is generated. S502: Call the unit concentration distribution value, and based on the spatial distance between the principal strain direction vector and the coordinates of adjacent units, use a differential method to compare the angle relationship between the concentration difference of adjacent units and the direction vector, calculate the concentration gradient of multiple units, and obtain the unit concentration gradient sequence. S503: Based on the unit concentration gradient sequence, call the concentration gradient values ​​of adjacent units to perform gradient reverse tracing, determine the transmission direction of the path nodes in spatial coordinates step by step, lock the starting end of the path and locate the coordinates of the peak point, and obtain the spatial positioning result of the pollution source.

Images