CN114841066A - Flying dust particle drift trajectory model and modeling method thereof - Google Patents

Flying dust particle drift trajectory model and modeling method thereof Download PDF

Info

Publication number
CN114841066A
CN114841066A CN202210465886.1A CN202210465886A CN114841066A CN 114841066 A CN114841066 A CN 114841066A CN 202210465886 A CN202210465886 A CN 202210465886A CN 114841066 A CN114841066 A CN 114841066A
Authority
CN
China
Prior art keywords
field
data
model
physical
training
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210465886.1A
Other languages
Chinese (zh)
Inventor
杨良怀
黄玉霖
王静雷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University of Technology ZJUT
Original Assignee
Zhejiang University of Technology ZJUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University of Technology ZJUT filed Critical Zhejiang University of Technology ZJUT
Priority to CN202210465886.1A priority Critical patent/CN114841066A/en
Publication of CN114841066A publication Critical patent/CN114841066A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D21/00—Measuring or testing not otherwise provided for
    • G01D21/02—Measuring two or more variables by means not covered by a single other subclass
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Evolutionary Biology (AREA)
  • Computer Hardware Design (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Health & Medical Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Biomedical Technology (AREA)
  • Geometry (AREA)
  • Algebra (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Fluid Mechanics (AREA)
  • Medical Informatics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a flying dust particle drift track model and a modeling method thereof, wherein the model comprises the following steps: 1) a data preprocessing stage: preprocessing data; 2) a model pre-training stage: dividing grids according to the geographical position of the station equipment; 3) a model pre-training stage: model pre-training is carried out through simulation data and a traditional physical method; 4) a transfer learning fine adjustment stage: and predicting the real flying dust particle drift track by using the trained model. The invention introduces the concept of physical field, replaces the calculation of speed field by the network of codec structure, and optimizes the calculation of flow term, diffusion term and volume force term. Experiments show that the raise dust track prediction model provided by the invention has greatly improved calculation speed and calculation stability compared with the traditional method, and has better universality and interpretability compared with a neural network method.

Description

Flying dust particle drift trajectory model and modeling method thereof
Technical Field
The invention relates to the technical field of flying dust trajectory prediction, in particular to a flying dust particle drift trajectory model and a modeling method thereof.
Background
As the rate of urban and social development continues to increase, more and more people are beginning to pay attention to environmental and residential issues. At this stage many developing countries are still facing the problem of air pollution and the latest GBD (global disease burden) analysis continues to identify air pollution as one of the most important risk factors leading to death and disability. Atmospheric particulates are a component of air pollution and are listed as the sixth most risky factor for early death.
Most of dust trajectory prediction models are based on traditional physical models or neural networks, most of traditional methods are based on experimental field research, numerical simulation research, meteorological field analysis research and data mining methods, and the traditional physical models are complex in calculation steps, huge in calculation amount and tend to be in an unstable state; neural network-based methods such as convolutional neural network and cyclic neural network prediction have their particular application areas and are also lacking in interpretability.
The invention introduces the concept of field, replaces the calculation of the velocity field in the traditional method by the network of the codec structure, and optimizes the calculation of flow term, diffusion term and volume force term. Experiments show that the dust trajectory prediction model provided by the invention has greatly improved accuracy and calculation speed compared with the traditional method, improves the calculation stability, and has better universality and interpretability compared with a neural network method.
Disclosure of Invention
The invention provides a flying dust particle drift trajectory model and a modeling method thereof, which mainly aim at the problems of heavy calculation task, lack of precision and difficulty in calculating a flying dust motion velocity field in the traditional method, improve the related traditional physical algorithm and combine a neural network method.
The technical scheme of the invention is as follows:
the dust emission data used in the present invention contains the following information: location information of equipment sites, equipment numbers, equipment longitudes, latitudes, equipment data states, temperatures, humidity, wind speeds, wind directions, PM2.5 concentrations, and pressures.
A flying dust particle drift track model and a modeling method thereof are characterized by comprising the following steps:
1) a data preprocessing stage: preprocessing data;
2) a model pre-training stage: dividing grids according to the geographical position of the station equipment;
3) a model pre-training stage: model pre-training is carried out through simulation data and a traditional physical method;
4) and (3) a transfer learning fine-tuning stage: and predicting the real drift track of the flying dust particles by using the trained model.
Further, 1) data preprocessing: in order to improve the accuracy of model prediction, the following data preprocessing work needs to be carried out:
1.1) data deduplication: removing repeated data uploaded in the acquisition process of the equipment through linked list query;
1.2) wind speed treatment: calculating to obtain a longitudinal component u of wind speed and a transverse component v of wind speed according to the wind speed and wind direction parameters in the raise dust data by taking geodetic longitude and latitude as coordinate axes;
1.2.1) calculating an included angle theta between the wind direction and the coordinate axis of the geodetic plane according to the wind direction;
1.2.2) calculating wind speed V calculating lateral component u and longitudinal component V:
the lateral component u ═ Vcos (θ); a longitudinal component v-Vsin (θ);
1.3) processing the longitude and latitude: in order to reduce the influence caused by the longitude and latitude errors of the station and improve the accuracy of a prediction result, the method adopts a Gaussian projection algorithm to convert (L, B) under a common geodetic coordinate into (x, y) under a Gaussian coordinate;
1.3.1) according to the longitude and latitude, calculating the longitude of the central meridian, the belt number and the longitude difference (the distance from a point to the central meridian), wherein the central meridian of 3 degrees is selected by the model, and the calculation process is as follows:
1.3.1.1) calculate the 3 ° band n 3 ：n 3 ＝int(L/3+0.5)；
1.3.1.2) calculating the 3 ° band n 3 Central meridian L of 0 ：L 0 ＝3n 3 ；
1.3.1.3) calculate the longitude difference l (the distance between the station and the central meridian): l ═ L-L 0 。
1.3.2) calculating the meridian arc length X according to the selection of the ellipsoid, and selecting a Clay-Laves-Fuji ellipsoid for calculation (the model can be applied to a 54 Beijing coordinate system):
X＝111134.861B-16036.480 sin 2B+16.828 sin 4B-0.022 sin 6B；
1.3.3) calculate the value of each parameter in the gaussian coordinate forward formula, wherein, B is the geodetic latitude, t is geodetic dimension tangent value, e is the ellipse eccentricity, N is the radius of curvature of the mortise-western circle, and a is the ellipsoid major radius parameter:
Figure BDA0003624071190000021
1.3.4) substituting the parameters obtained by the calculation into a Gaussian coordinate forward calculation formula to obtain coordinates (x, y) under a Gaussian coordinate system:
x＝X+N sin B cos Bl 2 +N sin B cos 3 B(5-t 2 -9e 2 cos 2 B)l 4
y＝N cos Bl+N cos 3 B(1-t 2 +e 2 cos 2 B)l 3 +N cos 5 B(5-18t 2 +t 4 )l 5
further, 2) grid processing:
in the dust raising movement, the mutual influence of adjacent stations on the space is very close, so that the dust raising law of a single station is not researched any more, and the distribution condition and the change law of the dust in the space are researched from the overall perspective through a speed field, a concentration field, a temperature field, a pressure field and a humidity field. To describe the physical field, further meshing according to the geographical location of the site device is required.
2.1) clustering analysis: in order to improve the accuracy of prediction, before grid division, DBSCAN clustering analysis is carried out according to the physical position of station equipment to obtain station equipment clustering clusters under Gaussian coordinates;
2.2) meshing: after a cluster with a better effect is obtained, dividing uniform square grids according to the geographical position of the station equipment in the cluster;
2.3) data interpolation: and performing linear interpolation according to the distance between the equipment stations and the dust data value, and finally storing the dust data in the grid in the form of a speed field, a concentration field, a temperature field, a pressure field and a humidity field.
Further, 3) model pre-training:
at present, two mainstream methods applied to predicting the drift track of dust particles are a forward calculation method based on a traditional physical method and a neural network method based on numerical simulation, but the two methods have certain disadvantages: the traditional physical method is too tedious in calculation and complex in steps; neural network methods lack interpretability and it is generally difficult to obtain a true velocity field as a training label. Therefore, the invention provides a flying dust particle drift track model and a modeling method thereof.
3.1) generating simulation data: simulation generation of T in divided grids 0 The dust simulation data at each moment comprises a speed field, a concentration field, a temperature field, a pressure field and a humidity field of the dust;
3.2) forward calculation based on the traditional physical method: on the basis of simulation data, updating a volume force term according to a traditional physical method (mainly according to an N-S Navier-Stokes equation), and carrying out forward calculation on the calculation of a flow term and a diffusion term to obtain a subsequent time node T 1 ，T 2 ，T 3 ，…，T n ，T n+1 The dust physical field augmentation dataset of (1):
Figure BDA0003624071190000041
3.2.1) N-S equation: the N-Snavier-Stokes equation (Navier-Stokes equations) is a fluid motion equation describing conservation of fluid momentum, and the specific formula is as follows:
Figure BDA0003624071190000042
wherein,
Figure BDA0003624071190000043
is the material derivative of the fluid with respect to the velocity field. The primary focus of the Euler view is the grid point in space, and hence the material derivative
Figure BDA0003624071190000044
Equal to the rate of change of fluid velocity over time at the grid point
Figure BDA0003624071190000045
And the rate of change u · u of the fluid under the action of the velocity field. Other parameter aspects: ρ is the fluid density; u is a velocity vector; p is pressure; f is the external force applied to the fluid per unit volume and μ is the viscosity of the fluid.
Figure BDA0003624071190000046
The above formula is the update process of the fluid state, W 0 Represents the initial state of the particle, and the state after the calculation of the volume force term (f) is W 1 Through convection term
Figure BDA0003624071190000047
Calculated resulting state W of 2 And finally by the diffusion term (μ +) 2 u) to obtain the state W 3 。
3.2.2) volume force term: the volume force term is calculated as f ═ ρ g (g is the acceleration of gravity and ρ is the fluid density).
3.2.3) convection term: in the fluid simulation, convection terms are automatically executed, and the essence of convection is that under the action of different speed fields, fluid micelles are mutually transferred between adjacent grid points and physical quantity values on the grid points are updated, and a given oneThe grid point is located at (x, y) and the physical quantity field value of the point is
Figure BDA0003624071190000048
The transverse component and the longitudinal component of the velocity field corresponding to the particles on the grid points are u and v respectively, and the physical quantity field value of the grid points after 1 time step delta t is solved
Figure BDA0003624071190000049
Comprises the following steps:
3.2.3.1) tracking 1 time step backwards to obtain the grid position (x-u delta t, y-v delta t) of the particle;
3.2.3.2) constructing a Catmull-Rom spline interpolation function according to the physical field grid to calculate and obtain the variation mean value of the particle physical quantity field
Figure BDA0003624071190000051
3.2.3.3) end up
Figure BDA0003624071190000052
The calculation formula of (2) is as follows:
Figure BDA0003624071190000053
3.2.4) diffusion term: the calculation of the diffusion term is another very important part of the fluid simulation. The calculation formula of the diffusion term of the fluid is a v ＝μ▽ 2 u, wherein a v Mu is the viscosity coefficient of the fluid, which is the acceleration of the particles under the action of the fluid's viscous force. Suppose a physical quantity at grid coordinates (x, y) at time n +1
Figure BDA0003624071190000054
It is known that the physical quantity at the nth moment needs to be solved
Figure BDA0003624071190000055
The specific calculation steps are as follows:
3.2.4.1) establishing a backward Euler equation
Figure BDA0003624071190000056
Wherein Δ t is the time step;
3.2.4.2) calculating the Laplace operator by adopting a center difference method:
Figure BDA0003624071190000057
wherein
Figure BDA0003624071190000058
Is composed of
Figure BDA0003624071190000059
Physical quantity around the grid, wherein delta x is grid width, and delta y is grid height, a square grid is used in the method, and inherent delta x is delta y;
3.2.4.3) order
Figure BDA00036240711900000510
The constant term is proposed to obtain:
Figure BDA00036240711900000511
3.2.4.4) converting the matrix into a matrix form, and solving the large-scale positive definite sparse matrix by adopting a conjugate gradient method.
Figure BDA00036240711900000512
(3.2.5) according to the volume force term, the convection term and the diffusion term, the subsequent time node T is obtained by updating and calculating 1 ，T 2 ，…，T n ，T n+1 Temporal raise dust physical field data
Figure BDA0003624071190000061
Figure BDA0003624071190000062
Thereby augmenting the data set
Figure BDA0003624071190000063
Pre-training is performed.
(3.3) pre-training: calculating T obtained in section (3.2) by a traditional physical method 1 ，T 2 ，…，T n ，T n+1 Dust physical information of time node
Figure BDA0003624071190000064
Inputting the data set into a neural network model U for pre-training:
(3.3.1) model selection: the invention introduces the concept of grid and field when researching the dust emission data, and the data set D is essentially a group of space distribution data expressing the time variation of various physical variable fields, so that the data set D can be analogized to pixel distribution data of an image. Therefore, compared with other networks, the U-net model with excellent performance in the image segmentation prediction field is more adaptive to the dust physical field data, so that the U-net model is selected as the pre-training model.
(3.3.2) training details: the input channel size is 32 × 32 × 6 × 2, where 32 × 32 is the grid size, and 6 × 2 represents 6 pieces of physical field data (velocity field lateral component, velocity field longitudinal component, concentration field, temperature field, pressure field, humidity field) and two adjacent time fields. The output channel is 32 × 32 × 3 × 1, where 32 × 32 is the grid size, and 3 × 1 represents 3 physical fields (velocity field lateral component, velocity field longitudinal component, and concentration field) and 1 temporal field generated by prediction.
Further, 4) migration learning and fine tuning:
and inputting the real data set into a model U which is subjected to pre-training, and adjusting the iteration number (epoch) training batch (batch-size) and the loss function (loss) to obtain a final flying dust particle drift prediction track.
(4.1) model migration: reserving all parameters in the pre-trained U-Net model, and taking the parameters as a training model of a real data set;
(4.2) data training: the input channel size is 32 × 32 × 4 × 2, where 32 × 32 is the grid size, and 4 × 2 represents 4 real physical field data (concentration field, temperature field, pressure field, humidity field) and two adjacent time fields. The output channel is 32 × 32 × 3 × 1, where 32 × 32 is the grid size, and 3 × 1 represents 3 physical fields (velocity field lateral component, velocity field longitudinal component, and concentration field) and 1 time field generated by prediction;
(4.3) model fine-tuning: the value of the adjustment parameter epoch is 240, the value of the training batch parameter batch size is 8, and L2 is set as a loss function, so that the final flying dust drift trajectory is obtained through training. Where the L2 loss function is the Least Squares Error (LSE) function, the purpose of using LSE is to minimize the true value y i And training result f (x) i ) The sum of squares of the differences D between L2 The function formula is as follows:
Figure BDA0003624071190000071
the invention has the following beneficial effects: compared with the traditional physical method, the model of the invention is superior to the traditional method in the aspects of accuracy, calculation time and stability, and has better universality and stronger interpretability compared with a numerical estimation method only using a neural network, so that the model can be more widely applied to most scenes.
Drawings
FIG. 1 is an overall architecture diagram of the present invention;
FIG. 2 is a schematic representation of Gaussian coordinates of the present invention;
FIG. 3 is a schematic diagram of a grid of the present invention; the white gradation of each small square in fig. 3 represents the size of the density field, and the higher the white gradation, the higher the dust density in the grid;
FIG. 4 is a schematic diagram of the semi-Lagrangian method of the present invention;
FIG. 5 is a calculated dust velocity field according to a conventional method of the present invention;
FIG. 6 is a dust concentration field calculated by the conventional method of the present invention; fig. 6 is composed of 32 × 32 small squares, the white scale of each small square represents the size of the concentration field, and the higher the white scale, the higher the dust concentration in the grid;
fig. 7 is a schematic diagram of a U-net network architecture used in the present invention.
Detailed Description
The invention is further described below with reference to the accompanying drawings.
Referring to fig. 1-7, a flying dust particle drift trajectory model and a modeling method thereof specifically include the following steps:
1) data preprocessing:
(1.1) removing the weight: checking the connection table, and screening out repeated dust raising data uploaded by equipment;
(1.2) wind speed processing:
(1.2.1) calculating an included angle theta between the wind direction and the coordinate axis of the earth plane according to the wind direction;
(1.2.2) calculating wind speed V calculating lateral component u and longitudinal component V:
the lateral component u ═ Vcos (θ); the longitudinal component v is Vsin (θ).
(1.3) Gaussian projection is carried out on the site coordinates: the objective of this step is to convert (L, B) in longitude and latitude coordinates into (x, y) in gaussian coordinates, and the specific steps can refer to fig. 2;
(1.3.1) calculating the longitude and the belt number of the central meridian and the longitude difference (the distance from a point to the central meridian) according to the longitude and the latitude, wherein the central meridian of 3 degrees is selected by the model, and the calculation process is as follows:
(1) calculating the 3 degree zone n 3 ：n 3 ＝int(L/3+0.5)；
(2) Calculating the 3 degree zone n 3 Central meridian L of 0 ：L 0 ＝3n 3 ；
(3) Calculate the longitude difference l (the distance of the station and the central meridian): l ═ L-L 0 ；
(1.3.2) calculating the meridian arc length X according to selection of an ellipsoid, and selecting a Clay-Laves-Fuji ellipsoid for calculation (the model can be applied to a 54 Beijing coordinate system):
X＝111134.861B-16036.480 sin 2B+16.828 sin 4B-0.022 sin 6B；
(1.3.3) calculating the value of each parameter in the Gaussian coordinate forward calculation formula, wherein B is the geodetic latitude, t is the tangent value of the geodetic dimension, e is the ellipse eccentricity, N is the curvature radius of the Mao-West circle, and a is the ellipsoid long radius parameter:
Figure BDA0003624071190000081
(1.3.4) substituting the parameters obtained by the calculation into a Gaussian coordinate forward calculation formula to obtain coordinates (x, y) under a Gaussian coordinate system:
x＝X+N sin B cos Bl 2 +N sin B cos 3 B(5-t 2 -9e 2 cos 2 B)l 4
y＝N cos Bl+N cos 3 B(1-t 2 +e 2 cos 2 B)l 3 +N cos 5 B(5-18t 2 +t 4 )l 5
2) grid processing:
(2.1) clustering analysis: obtaining a cluster C through a DBscan clustering algorithm, wherein the steps are as follows:
(1) marking all objects as univisified; randomly selecting an unvisited object P;
(2) if the epsilon field of P has at least MinPts objects, a new cluster C is created, and P is added to the cluster C;
(3) let N be the object set in epsilon field of P, if there are at least MinPts objects in epsilon field of P, add these objects to N; if P is not already a member of any cluster, P is added to cluster C.
(2.2) linear interpolation to obtain a uniform grid (the grid diagram can refer to the attached figure 3):
(2.2.1) firstly finding out 3 points around the interpolation point according to the Delaunay triangulation algorithm
Triangle list triangles (pseudo code as follows):
Figure BDA0003624071190000091
(2.2.2) calling a griddata function to perform data interpolation to obtain a dust raising speed field, a concentration field, a temperature field, a pressure field and a humidity field which are represented by a uniform grid.
3) Pre-training a model:
(3.1) simulation of initial data: 23 in 20 days and night of 2021 year 07 month 20: dust data at time 00 as T 0 Simulated initial data of time of day.
(3.2) forward calculation by a physical traditional method: the method mainly relies on the Navier-Stokes equations (Navier-Stokes equations) in the forward calculation of the traditional method, and is called N-S equations for short.
(3.2.1) volume force term: the updated volume force term f is ρ g (g is the acceleration of gravity and ρ is the density).
(3.2.2) convection term (convection process can refer to FIG. 4): given a grid point location as (x, y), the physical quantity field value of the point is
Figure BDA0003624071190000101
The transverse component and the longitudinal component of the velocity field corresponding to the particles on the grid points are u and v respectively, and the physical quantity field value of the grid points after 1 time step delta t is solved
Figure BDA0003624071190000102
Comprises the following steps:
(1) tracking 1 time step backwards to obtain the grid position (x-u delta t, y-v delta t) of the particle;
(2) constructing a Catmull-Rom spline interpolation function according to the physical field grid to calculate the variation mean value of the particle physical quantity field
Figure BDA0003624071190000103
(3) Finally, the product is processed
Figure BDA0003624071190000104
The calculation formula of (2) is as follows:
Figure BDA0003624071190000105
(3.2.3) diffusion term: the main purpose of this term is to solve for the viscosity term a of the fluid v ＝μ▽ 2 u, wherein a v Mu is the viscosity coefficient of the fluid, which is the acceleration of the particles under the action of the fluid's viscous force. Assume that at time n +1, the grid coordinates (x, y) are physicalMeasurement of
Figure BDA0003624071190000106
It is known that the physical quantity at the nth moment needs to be solved
Figure BDA0003624071190000107
The specific calculation steps are as follows:
(1) establishing a backward Euler equation
Figure BDA0003624071190000108
Wherein Δ t is the time step;
(2) calculating a Laplace operator by adopting a center difference method:
Figure BDA0003624071190000109
wherein
Figure BDA00036240711900001010
Is composed of
Figure BDA00036240711900001011
Physical quantity around the grid, wherein delta x is grid width, and delta y is grid height, a square grid is used in the method, and inherent delta x is delta y;
(3) order to
Figure BDA00036240711900001012
The constant term is proposed to obtain:
Figure BDA00036240711900001013
(4) and converting the matrix form into a matrix form, and solving the large-scale positive definite sparse matrix by adopting a conjugate gradient method.
Figure BDA0003624071190000111
(3.2.4) according to the volume force term, forFlow item, diffusion item, and subsequent time node T calculated by traditional method 1 ，T 2 ，T 3 ，…，T n ，T n+1 Temporal raise dust physical field data
Figure BDA0003624071190000112
Figure BDA0003624071190000113
Thereby augmenting the data set
Figure BDA0003624071190000114
Training is performed to obtain a pre-training model U, and the velocity field and the concentration field in the augmented data set can be referred to in the attached figures 5 and 6.
(3.3) pre-training process:
(3.3.1) network selection:
based on the characteristics of flying dust drift motion (in space, flying dust data of adjacent stations have strong correlation), the invention introduces the concepts of grids and fields when researching the flying dust data, and the data set D is essentially a group of spatial distribution data expressing the time variation of various physical variable fields, so that the data set D can be analogized to pixel distribution data of an image. Therefore, compared with other networks, the U-net model with excellent performance in the image segmentation prediction field is more adaptive to the dust physical field data in the overall structure, so that the U-net model is selected as the pre-training model (the specific structure diagram of the U-net model can refer to the attached figure 7).
(3.3.2) training details:
(1) the pre-training selects a U-Net model, the value of the parameter epoch is set to 240, the value of the parameter batch size is set to 8, the size of the data set is 2000 groups, and the proportion of the training set, the test set and the verification set is 12: 3: 1, the regression index selects the L2 loss function.
(2) Input channel size is 32 × 32 × 6 × 2: where 32 × 32 is the grid size and 6 × 2 represents 6 physical field data (velocity field lateral component, velocity field longitudinal component, concentration field, temperature field, pressure field, humidity field) and two adjacent time fields. The output channels are 32 × 32 × 3 × 1: where 32 × 32 is the grid size and 3 × 1 represents 3 physical fields (velocity field lateral component, velocity field longitudinal component, and concentration field) and 1 time field generated by prediction.
4) Model migration and fine tuning:
(4.1) model migration: reserving all parameters in the pre-trained U-Net model, and taking the parameters as a training model of a real data set;
(4.2) data training: the input channel size is 32 × 32 × 4 × 2, where 32 × 32 is the grid size, and 4 × 2 represents 4 real physical field data (concentration field, temperature field, pressure field, humidity field) and two adjacent time fields. The output channel is 32 × 32 × 3 × 1, where 32 × 32 is the grid size, and 3 × 1 represents 3 physical fields (velocity field transverse component, velocity field longitudinal component, and concentration field) and 1 temporal field generated by prediction;
(4.3) model fine-tuning: the value of the adjustment parameter epoch is 240, the value of the training batch parameter batch size is 8, and L2 is set as a loss function, so that the final flying dust drift trajectory is obtained through training. Where the L2 loss function is the Least Squares Error (LSE) function, the purpose of using LSE is to minimize the true value y i And training result f (x) i ) Sum of squares of difference between D L2 The function formula is as follows:
Figure BDA0003624071190000121
the invention provides a flying dust particle drift trajectory model and a modeling method thereof by combining a codec neural network and an optimized traditional method aiming at the problems that the traditional method is too large in calculated amount, poor in precision and stability and poor in interpretability of the neural network method.
The embodiments described in this specification are merely illustrative of implementations of the inventive concept and the scope of the present invention should not be considered limited to the specific forms set forth in the embodiments but rather by the equivalents thereof as may occur to those skilled in the art upon consideration of the present inventive concept.

Claims (5)

1. A flying dust particle drift track model and a modeling method thereof are characterized by comprising the following steps:
1) a data preprocessing stage: preprocessing data;
2) a model pre-training stage: dividing grids according to the geographical position of the station equipment;
3) a model pre-training stage: model pre-training is carried out through simulation data and a traditional physical method;
4) and (3) a transfer learning fine-tuning stage: and predicting the real flying dust particle drift track by using the trained model.
2. The model of flying dust particle drift trajectory and its modeling method as claimed in claim 1, wherein said data preprocessing stage comprises the steps of:
1.1) data duplication stage: screening dust data repeatedly uploaded by site equipment through link table query;
1.2) wind speed treatment stage: according to the wind speed and the wind direction parameters in the dust raising data, calculating the transverse component and the longitudinal component of the wind speed, and specifically executing the following steps:
1.2.1) calculating an included angle theta between the wind direction and the coordinate axis of the geodetic plane according to the wind direction;
1.2.2) calculating the wind speed V:
u＝V cos(θ)，v＝V sin(θ)；
wherein u represents a transverse component and v represents a longitudinal component;
1.3) Gaussian projection phase of geographic coordinates: transforming longitude and latitude coordinates (L, B) of the station into Gaussian coordinates (x, y) through Gaussian projection calculation, wherein L is geodetic longitude, and B is geodetic latitude; the specific execution steps are as follows:
1.3.1) calculating the longitude, the belt number and the longitude difference of the central meridian according to the longitude and the latitude, wherein the central meridian with 3 degrees is selected as a model, and the method comprises the following steps:
1.3.1.1) calculating a 3 deg. band number n 3 ：n 3 ＝int(L/3+0.5)；
1.3.1.2) calculating the band number n of 3 ° 3 Central meridian L of 0 ：L 0 ＝3n 3 ；
1.3.1.3) calculate the longitude difference l (the distance between the station and the central meridian): l ═ L-L 0 ；
1.3.2) calculating the meridian arc length X according to the selection of an ellipsoid; the model selects a Clay-Laves-Fuji ellipsoid for calculation:
X＝111134.861B-16036.480sin2B+16.828sin4B-0.022sin6B；
1.3.3) calculate the value of each parameter in the gaussian coordinate forward formula, wherein, B is the geodetic latitude, t is geodetic dimension tangent value, e is the ellipse eccentricity, N is the radius of curvature of the mortise-western circle, and a is the ellipsoid major radius parameter:
Figure FDA0003624071180000011
1.3.4) solving the coordinates (x, y) under the Gaussian coordinate system by using a Gaussian coordinate forward formula:
x＝X+N sinB cosBl 2 +N sinB cos 3 B(5-t 2 -9e 2 cos 2 B)l 4
y＝N cosBl+N cos 3 B(1-t 2 +e 2 cos 2 B)l 3 +Ncos 5 B(5-18t 2 +t 4 )l 5 。
3. the flying dust particle drift trajectory model and the modeling method thereof as claimed in claim 1, wherein said mesh processing stage comprises the steps of:
2.1) clustering analysis stage: performing DBSCAN clustering analysis according to the physical positions of the station equipment before grid division to obtain station equipment clustering clusters under Gaussian coordinates;
2.2) a meshing stage: after the cluster with better effect is obtained in the step 2.1), dividing uniform square grids according to the geographical position of the station equipment in the cluster;
2.3) data interpolation stage: and performing linear interpolation according to the distance between the equipment stations and the dust data value, and finally storing the dust data in the grid in the form of a speed field, a concentration field, a temperature field, a pressure field and a humidity field.
4. The flying dust particle drift trajectory model and the modeling method thereof as claimed in claim 1, wherein the network training phase comprises the steps of:
3.1) simulation data phase: simulating and generating flying dust simulation data at the T0 moment in the divided grids, wherein the flying dust simulation data comprise a speed field, a concentration field, a temperature field, a pressure field and a humidity field of the flying dust;
3.2) forward computing stage of traditional physical method: on the basis of the simulation data, according to the traditional physical method, namely, according to the N-S Navie-Stokes equation, the calculation of the volume force item f, the convection item and the diffusion item is updated, forward prediction is carried out, and the dust physical field data set of the subsequent time nodes T1, T2, T3, …, Tn and Tn +1 is obtained
Figure FDA0003624071180000021
The specific execution steps are as follows:
3.2.1) volume force term f calculation: f is rho g, wherein g is the acceleration of gravity and rho is the density of the fluid;
3.2.2) convection term calculation: given a grid point location as (x, y), the physical quantity field value of the point is
Figure FDA0003624071180000022
The transverse component and the longitudinal component of the velocity field corresponding to the particles on the grid points are u and v respectively, and the physical quantity field value of the grid points after 1 time step delta t is solved
Figure FDA0003624071180000023
The related steps are as follows:
3.2.2.1) tracking 1 time step backwards to obtain the grid position (x-u delta t, y-v delta t) of the particle;
3.2.2.2) constructing a Catmull-Rom spline interpolation function according to the physical field grid to calculate and obtain the variation mean value of the particle physical quantity field
Figure FDA0003624071180000024
3.2.2.3) Final
Figure FDA0003624071180000025
The calculation formula of (2) is as follows:
Figure FDA0003624071180000026
3.2.3) diffusion term calculation: the main purpose of this step is to solve the viscosity term a of the fluid v ＝μ▽ 2 u, wherein a v The acceleration of the particles under the action of the fluid viscosity force is shown, and mu is the viscosity coefficient of the fluid; suppose a physical quantity at grid coordinates (x, y) at time n +1
Figure FDA0003624071180000031
It is known that the physical quantity at the nth time needs to be solved
Figure FDA0003624071180000032
The specific calculation steps are as follows:
3.2.3.1) establishing a backward Euler equation
Figure FDA0003624071180000033
Wherein Δ t is the time step;
3.2.3.2) calculating the Laplace operator by adopting a center difference method:
Figure FDA0003624071180000034
wherein
Figure FDA0003624071180000035
Is composed of
Figure FDA0003624071180000036
Physical quantity around the grid, wherein delta x is the grid width, and delta y is the grid height;
3.2.3.3) order
Figure FDA0003624071180000037
The constant term is proposed to obtain:
Figure FDA0003624071180000038
3.2.3.4) transforming the matrix into a matrix form, and solving the large-scale positive definite sparse matrix by adopting a conjugate gradient method:
Figure FDA0003624071180000039
3.3) generating data set stage: according to the volume force item, the convection item and the diffusion item, updating and calculating to obtain a subsequent time node T 1 ，T 2 ，…，T n ，T n+1 Temporal raise dust physical field data
Figure FDA00036240711800000310
And then from the data set
Figure FDA00036240711800000311
Pre-training is carried out;
3.4) pre-training phase: calculating T obtained in the step C3) by a traditional physical method 1 ，T 2 ，T 3 ，…，T n ，T n+1 Dust physical information of time node
Figure FDA00036240711800000312
Pre-training a neural network model U as a data set:
3.4.1) model selection: selecting U-Net as a training network;
3.4.2) training details: the size of an input channel is 32 multiplied by 6 multiplied by 2, wherein 32 multiplied by 32 is the size of a grid, 6 multiplied by 2 represents 2 adjacent time fields of 6 physical field data, and the 6 physical field data are respectively a velocity field transverse component, a velocity field longitudinal component, a concentration field, a temperature field, a pressure field and a humidity field; the output channel is 32 × 32 × 3 × 1, where 32 × 32 is the grid size, 3 × 1 represents 3 physical fields and 1 time field generated by prediction, and the 3 physical fields are a velocity field transverse component, a velocity field longitudinal component, and a concentration field, respectively.
5. The flying dust particle drift trajectory model and the modeling method thereof according to claim 4, wherein the migration training and fine tuning phase comprises the following steps:
4.1) model migration phase: reserving all parameters in the U-Net model which is subjected to pre-training in the step 3.4), and taking the parameters as a training model of a real data set;
4.2) data training phase: the size of an input channel is 32 multiplied by 4 multiplied by 2, wherein 32 multiplied by 32 is the size of a grid, 4 multiplied by 2 represents 4 real physical field data and 2 adjacent time fields, and the 4 real physical field data are a concentration field, a temperature field, a pressure field and a humidity field respectively; the output channel is 32 × 32 × 3 × 1, where 32 × 32 is the size of the grid, 3 × 1 represents 3 physical fields and 1 time field generated by prediction, and 3 physical fields are respectively a velocity field transverse component, a velocity field longitudinal component and a concentration field;
4.3) fine tuning stage of the model: the value of the adjustment parameter epoch is 240, the value of the training batch parameter batch is 8, and L2 is set as the loss function, and the prediction result is finally obtained. Wherein the loss function L2 employs a Least Squares Error (LSE) function, the purpose of using LSE is to minimize the true value y i And training result f (x) i ) The sum of squares of the differences D between L2 The function formula is
Figure FDA0003624071180000041
CN202210465886.1A 2022-04-29 2022-04-29 Flying dust particle drift trajectory model and modeling method thereof Pending CN114841066A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210465886.1A CN114841066A (en) 2022-04-29 2022-04-29 Flying dust particle drift trajectory model and modeling method thereof

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210465886.1A CN114841066A (en) 2022-04-29 2022-04-29 Flying dust particle drift trajectory model and modeling method thereof

Publications (1)

Publication Number Publication Date
CN114841066A true CN114841066A (en) 2022-08-02

Family

ID=82568616

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210465886.1A Pending CN114841066A (en) 2022-04-29 2022-04-29 Flying dust particle drift trajectory model and modeling method thereof

Country Status (1)

Country Link
CN (1) CN114841066A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115794447A (en) * 2023-02-07 2023-03-14 青岛哈尔滨工程大学创新发展中心 Grid data transmission method for multi-physical field coupling

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115794447A (en) * 2023-02-07 2023-03-14 青岛哈尔滨工程大学创新发展中心 Grid data transmission method for multi-physical field coupling
CN115794447B (en) * 2023-02-07 2023-04-28 青岛哈尔滨工程大学创新发展中心 Grid data transmission method for multi-physical field coupling

Similar Documents

Publication Publication Date Title
Hengl et al. Geomorphometry: concepts, software, applications
Technitis et al. From A to B, randomly: a point-to-point random trajectory generator for animal movement
CN112862090B (en) Air temperature forecasting method based on deep space-time neural network
CN107688906B (en) Multi-method fused transmission line meteorological element downscaling analysis system and method
Gullu Coordinate transformation by radial basis function neural network
CN105554873A (en) Wireless sensor network positioning algorithm based on PSO-GA-RBF-HOP
Konakoglu et al. 2D coordinate transformation using artificial neural networks
Ziggah et al. Performance evaluation of artificial neural networks for planimetric coordinate transformation—a case study, Ghana
CN114841066A (en) Flying dust particle drift trajectory model and modeling method thereof
CN112712169A (en) Model building method and application of full residual depth network based on graph convolution
CN111310344B (en) Method for considering coupling effect of wind field and fire field in forest fire spreading simulation
Wong et al. Graph neural network based surrogate model of physics simulations for geometry design
Zhao et al. Examining land-use change trends in yucheng district, Ya’an city, China, using ANN-CA modeling
Zhang et al. The CA model based on data assimilation
CN111914465B (en) Clustering and particle swarm optimization-based method for calibrating hydrologic parameters of non-data region
Skalova et al. Comparison of three regression models for determining water retention curves
CN116703008A (en) Traffic volume prediction method, equipment and medium for newly built highway
Peprah et al. Appraisal of methods for estimating orthometric heights–A case study in a mine
Girot “Cloudism”: Towards a new culture of making landscapes
Ma et al. A Method for Establishing Tropospheric Atmospheric Refractivity Profile Model Based on Multiquadric RBF and k-means Clustering
Wang et al. Real-time rapid prediction of variations of Earth’s rotational rate
CN108121869B (en) TIN _ DDM buffer surface rapid construction method based on rolling ball model
Guan et al. Rain fall predict and comparing research based on Arcgis and BP neural network
Mandelman et al. Evaluation of a finite element formulation for the shallow water equations with numerical smoothing in the Gulf of San Jorge
CN117556719B (en) Hydrodynamic model verification method based on natural tracer ions

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination