Industrial park air pollution source tracing method and device

Abstract

The invention discloses an industrial park air pollution source tracing method and device, and the method comprises the steps: carrying out the grid division of an industrial park, and carrying out the interpolation calculation of the pollutant concentration corresponding to each grid of the industrial park; constructing a first objective function, and fitting parameters based on the pollutant concentration corresponding to each grid of the industrial park to obtain a first fitting result; constructing a second objective function, and fitting the parameters based on the pollutant concentration corresponding to each grid of the industrial park and the first fitting result to obtain a second fitting result; and determining the position of the air pollution source in the industrial park based on the second fitting result. According to the method, the coordinates of the industrial park pollution source emission chimney can be accurately positioned, and meanwhile, the leakage source intensity, the leakage source height and the foreign pollutant concentration are calculated.

Description

Technical Field

[0001] The present invention relates to the technical field of pollution treatment, and in particular to a method and a device for tracing the source of air pollution in an industrial park. Background Art

[0002] Industrial waste gas emissions in industrial areas are the primary cause of air pollution, and the air quality problem in industrial areas with dense industrial production enterprises is particularly severe. Most industrial waste gas emissions are discharged through factory chimneys and can be regarded as single pollution source emissions. The industrial waste gas of factories in the park includes SO2, NOx, COD, ammonia nitrogen, VOCs, PM2.5, PM10, etc. The monitoring stations in the park can monitor the indicators of each pollutant separately. In addition to the industrial waste gas emissions of factories in the park, the source of pollutants may also be affected by external pollutants.

[0003] How to locate and trace the sources of air pollution in industrial parks is an urgent problem that needs to be solved. Summary of the invention

[0004] The present invention provides a method and device for tracing the air pollution sources in an industrial park, which can solve the technical problem of locating and tracing the air pollution sources in the industrial park.

[0005] In various method embodiments of the present invention, a method for tracing the source of air pollution in an industrial park includes:

[0006] Step S1: regularly obtain monitoring data of the industrial park, and establish an initial industrial park coordinate system Oxy for the industrial park, wherein the origin of the initial industrial park coordinate system is the geographical center of the industrial park, the horizontal coordinate is parallel to the latitude, and the vertical coordinate is parallel to the longitude; when new monitoring data is obtained, proceed to step S2; otherwise, proceed to step S1;

[0007] Step S2: taking the new monitoring data as the current monitoring data, based on the current monitoring data, taking the monitoring site with the largest pollutant concentration as the coordinate origin, taking the wind direction as the positive direction of the x-axis, and performing coordinate transformation on the initial industrial park coordinate system Oxy to obtain the current industrial park coordinate system Ox'y';

[0008] Step S3: Based on the current industrial park coordinate system Ox'y', the industrial park is divided into grids, and the pollutant concentration corresponding to each grid of the industrial park is interpolated and calculated;

[0009] Step S4: Construct a first objective function, based on the pollutant concentrations corresponding to each grid in the industrial park, and calculate the parameter θ i =(u i ,σ i ,ξ i ,w i) to perform fitting, and obtain a first fitting result;

[0010] Among them, Θ i is a parameter set, i=1, 2, …, n, n is the number of columns in the grid, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The variance of i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid;

[0011] Step S5: Construct a second objective function, based on the pollutant concentrations corresponding to each grid in the industrial park and the first fitting result, and adjust the parameters Г=(Q, H0, tx v ) is fitted to obtain the second fitting result, where Г is the parameter set, Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of pollution source;

[0012] Step S6: Based on the second fitting result, determine the location of the air pollution source in the industrial park.

[0013] Optionally, the step S2 includes:

[0014] The new monitoring data is used as the current monitoring data, and based on the current monitoring data, the coordinates (x p ,y p ) as the coordinate origin, and the wind direction as the positive direction of the x-axis;

[0015] The initial industrial park coordinate system Oxy is transformed to obtain the current industrial park coordinate system Ox'y'. The coordinate transformation formula is:

[0016]

[0017] Wherein, x and y are the coordinates of the initial industrial park coordinate system Oxy respectively;

[0018] In the current industrial park coordinate system Ox'y', the monitoring site coordinates are , β is the The wind direction angle of the first monitoring station is The included angle in the positive longitude of the monitoring stations.

[0019] Optionally, the step S3 includes:

[0020] calculate ;

[0021] The area [X1, X2] × [Y1, Y2] is divided into n × m uniform grids, and the coordinates of the grids are marked as (tx i ,ty j ), i=1,2,…,n, j=1,2,…,m, where:

[0022]

[0023] The coordinates calculated by inverse distance interpolation are (tx i ,ty j )'s corresponding pollutant concentration:

[0024] .

[0025] Optionally, the step S4 includes:

[0026] Step S41: Constructing the first fitting function ;

[0027] Where x=tx i , represents the straight line where the i-th column grid is located, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The corresponding variance, ξ i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid;

[0028] Step S42: Construct a first objective function, the first objective function is ,in,

[0029] ,

[0030] f1(θ i ),…,f m (Θ i ) are all intermediate vectors, c i1 ,…,c im They are calculated using the inverse distance interpolation method (tx i ,ty1),…,(tx i ,ty m ) point;

[0031] Step S43: iteratively solving the first objective function to obtain a first fitting result.

[0032] Optionally, the step S43 includes:

[0033] Step S431: Initialize θ i =(u i ,σ i ,ξ i ,w i )’s initial value Θ i (0) =(u i (0) ,σ i (0) ,ξ i (0) ,w i (0) ), and initialize J (0) , set the current iteration step k to 0, set the value range of the first threshold ε to 0<ε<1, the value range of the second threshold θ to θ>0, the value range of the third threshold γ1 and the fourth threshold γ2 to 0<γ1<1<γ2, the value range of the fifth threshold η1 and the sixth threshold η2 to 0<η1<η2<1, initialize θ (0) ,θ (0) is the initial value of parameter θ, J (0) is the initialized Jacobian matrix;

[0034] Step S432: If , then the iteration ends and Θ i (k) As the first fitting result, where Θ i (k) =(u i (k) ,σ i (k) ,ξ i (k) ,w i (k) ),u i (k) , σ i (k) ,ξ i (k) , w i (k) That is the u obtained by fitting i , σ i ,ξ i , w i (estimated value or approximate value); otherwise, let k be k plus 1, and go to step S433; where J (k) is the Jacobian matrix corresponding to the kth iteration;

[0035] Step S433: Calculation , where I is the identity matrix, θ (k) is the kth iteration value of parameter θ, d(k) (θ (k) ) is the calculated value corresponding to the kth iteration, ,

[0036] calculate

[0037] in, , ρ is the intermediate variable,

[0038] ;

[0039] d (k) =d (k) (θ (k) );

[0040] Step S434: If , then let ;

[0041] like , then let θ (k+1) =θ (k) ;

[0042] like , then let θ (k+1) =γ1θ (k)

[0043] Step S435: Calculate Θ i (k+1) =θ i (k) +d (k) (θ (k) ), go to step S432.

[0044] Optionally, the step S5 includes:

[0045] Step S51: Calculation , tu is the estimated value of the vertical coordinate of the location of the air pollution source in the industrial park; calculation , is the estimated value of the concentration of foreign pollutants;

[0046] Step S52: construct a second fitting function, wherein the second fitting function is C(x, Q, H0, tx v ):

[0047]

[0048] Among them, the parameters to be fitted are Г=(Q,H0,tx v ), Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of the pollution source; μ is the wind speed, σ y , σz are the diffusion coefficients in the horizontal and vertical directions of the current industrial park coordinate system Ox'y', respectively, and the calculation formula is as follows:

[0049] , t y1 ,t y2 ,t z1 ,t z2 is the Bridges diffusion coefficient, the atmospheric stability type is determined according to meteorological conditions, and z is the diffusion parameter in the vertical upward direction perpendicular to the current industrial park coordinate system Ox'y';

[0050] Step S53: Construct a second objective function, the second objective function is ,in,

[0051]

[0052]

[0053] Among them, σ yi Indicates that tx i Substitute σ y The calculated result, σ zi Indicates that tx i Substitute σ z The calculation result is as follows:

[0054] ;

[0055] Step S54: iteratively solving the second objective function to obtain a second fitting result.

[0056] Optionally, the step S54 includes:

[0057] Step S541: Initialize Г=(Q,H0,tx v )’s initial value Г (0) =(Q (0) ,H0 (0) ,tx v (0) ), and initialize J1 (0) , set the current iteration step k1 to 0, set the value range of the seventh threshold ε1 to 0<ε1<1, the value range of the eighth threshold θ1 to θ1>0, the value range of the ninth threshold γ3 and the tenth threshold γ4 to 0<γ3<1<γ4, the value range of the eleventh threshold η3 and the twelfth threshold η4 to 0<η3<η4<1; initialize θ1 so that θ1 (0) >0,θ1 (0) is the initial value of parameter θ1;

[0058] Step S542: If , then the iteration ends and Г (k1) As the second fitting result, Г (k1) =(Q (k1) ,H0 (k1) ,tx v (k1) ), tx v (k1) As the second fitting result corresponding to tx v ; Otherwise, let k1 be k1 plus 1, and go to step S543; where J1 (k1) is the corresponding Jacobian matrix of the k1th iteration;

[0059] Step S543: Calculate , where I1 is the identity matrix, θ1 (k1) is the k1th iteration value of parameter θ1, d (k1) (θ (k1) ) is the calculated value corresponding to the k1th iteration;

[0060] ;

[0061]

[0062] calculate

[0063] in,

[0064]

[0065] d (k1) =d (k1) (θ1 (k1) );

[0066] Step S544: If , then let θ1 (k1+1) =γ4θ1 (k1) ;

[0067] like , then let θ1 (k1+1) =θ1 (k) ;

[0068] like , then let θ1 (k1+1) =γ3θ1 (k1) ;

[0069] Step S545: Calculate Г (k1+1) =G (k1) +d (k1) (θ1 (k1) ), enter step S542.

[0070] In the above-mentioned method embodiments of the present invention, a device for tracing the source of air pollution in an industrial park includes:

[0071] Data acquisition module: configured to periodically acquire monitoring data of the industrial park, establish an initial industrial park coordinate system Oxy for the industrial park, the origin of the initial industrial park coordinate system is the geographical center of the industrial park, the horizontal coordinate is parallel to the latitude, and the vertical coordinate is parallel to the longitude; when new monitoring data is acquired, the coordinate conversion module is triggered; otherwise, the data acquisition module is triggered;

[0072] Coordinate conversion module: configured to use the new monitoring data as current monitoring data, based on the current monitoring data, use the monitoring site with the largest pollutant concentration as the coordinate origin, use the wind direction as the positive direction of the x-axis, and perform coordinate transformation on the initial industrial park coordinate system Oxy to obtain the current industrial park coordinate system Ox'y';

[0073] A grid division module: configured to divide the industrial park into grids based on the current industrial park coordinate system Ox'y', and interpolate and calculate the pollutant concentration corresponding to each grid of the industrial park;

[0074] The first fitting module is configured to construct a first objective function, based on the pollutant concentrations corresponding to each grid of the industrial park, and to adjust the parameter θ i =(u i ,σ i ,ξ i ,w i ) to perform fitting, and obtain a first fitting result;

[0075] Among them, Θ i is a parameter set, i=1, 2, …, n, n is the number of columns in the grid, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The variance of i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid;

[0076] The second fitting module is configured to construct a second objective function, based on the pollutant concentrations corresponding to each grid of the industrial park and the first fitting result, to adjust the parameter Г=(Q,H0,tx v ) is fitted to obtain the second fitting result, where Г is the parameter set, Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of pollution source;

[0077] Pollution source determination module: configured to determine the location of the air pollution source in the industrial park based on the second fitting result.

[0078] In the above method embodiments of the present invention, a computer-readable storage medium is provided, wherein a plurality of instructions are stored in the storage medium, and the plurality of instructions are used for a processor to load and execute the above method.

[0079] In the above-mentioned method embodiments of the present invention, an electronic device includes: a processor for executing multiple instructions; a memory for storing multiple instructions; wherein the multiple instructions are used to be stored by the memory, and loaded by the processor to execute the method as described above.

[0080] The present invention has the following advantages: the pollutant concentration at the grid point is calculated by using the inverse distance interpolation method, and the air pollution source of the industrial park is traced by fitting the unknown parameters through the twice LM (Levenberg-Marquardt) algorithm. The present invention can accurately locate the coordinates of the chimneys of the pollution source in the industrial park, and simultaneously calculate the leakage source strength, leakage source height, and external pollutant concentration.

[0081] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] The above and other purposes, features and advantages of the present invention will become more apparent by describing the embodiments of the present invention in more detail in conjunction with the accompanying drawings. The accompanying drawings are used to provide a further understanding of the embodiments of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings, the same reference numerals generally represent the same components or steps.

[0083] Figure 1 This is a flow chart of a method for tracing the source of air pollution in an industrial park according to the present invention;

[0084] Figure 2 It is another schematic diagram of the process of tracing the air pollution source of the industrial park of the present invention;

[0085] Figure 3 This is a schematic diagram of the method flow of the first fitting of the present invention;

[0086] Figure 4 This is a schematic flow chart of the second fitting method of the present invention;

[0087] Figure 5 This is a schematic diagram of the structure of the device for tracing the source of air pollution in an industrial park according to the present invention;

[0088] Figure 6This is a schematic diagram of the structure of the electronic equipment for tracing the air pollution source in the industrial park according to the present invention. DETAILED DESCRIPTION

[0089] Below, the example embodiments according to the present invention will be described in detail with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments of the present invention, and it should be understood that the present invention is not limited by the example embodiments described herein. It should be noted that unless otherwise specifically stated, the relative arrangement of the parts and steps, the numerical expressions and the numerical values ​​set forth in these embodiments do not limit the scope of the present invention.

[0090] Those skilled in the art can understand that the terms "first", "second", S1, S2, etc. in the embodiments of the present invention are only used to distinguish different steps, devices or modules, etc., and neither represent any specific technical meaning nor represent the necessary logical order between them. It should also be understood that in the embodiments of the present invention, "multiple" can refer to two or more, and "at least one" can refer to one, two or more. It should also be understood that for any component, data or structure mentioned in the embodiments of the present invention, in the absence of explicit limitation or contrary revelation given in the context, it can generally be understood as one or more. In addition, the term "and / or" in the present invention is only a description of the association relationship of the associated objects, indicating that there can be three relationships, for example, A and / or B can represent: A exists alone, A and B exist at the same time, and B exists alone. In addition, the character " / " in the present invention generally indicates that the associated objects before and after are in an "or" relationship. It should also be understood that the description of each embodiment of the present invention emphasizes the differences between the embodiments, and the same or similar parts can refer to each other. For the sake of brevity, they will not be repeated one by one. At the same time, it should be understood that for the convenience of description, the sizes of the various parts shown in the drawings are not drawn according to the actual proportional relationship. The following description of at least one exemplary embodiment is merely illustrative in nature and is in no way intended to limit the invention and its application or use. Technologies, methods, and devices known to those of ordinary skill in the relevant art may not be discussed in detail, but where appropriate, they should be considered part of the specification. It should be noted that similar reference numerals and letters represent similar items in the following figures, so once an item is defined in one figure, it does not need to be further discussed in subsequent figures.

[0091] Embodiments of the present invention can be applied to electronic devices such as terminal devices, computer systems, servers, etc., which can operate with many other general or special computing system environments or configurations. Examples of well-known terminal devices, computing systems, environments and / or configurations suitable for use with terminal devices, computer systems, servers, etc. include, but are not limited to: personal computer systems, server computer systems, thin clients, thick clients, handheld or laptop devices, microprocessor-based systems, set-top boxes, programmable consumer electronics, network personal computers, small computer systems, large computer systems, and distributed cloud computing technology environments including any of the above systems, etc. Terminal devices, computer systems, servers, etc. Electronic devices can be described in the general context of computer system executable instructions (such as program modules) executed by computer systems. Generally, program modules can include routines, programs, target programs, components, logic, data structures, etc., which perform specific tasks or implement specific abstract data types. Computer systems / servers can be implemented in distributed cloud computing environments, where tasks are performed by remote processing devices linked through a communication network. In a distributed cloud computing environment, program modules can be located on local or remote computing system storage media including storage devices.

[0092] Exemplary Methods

[0093] Figure 1 FIG. 1 is a flow chart of a method for tracing the source of air pollution in an industrial park provided by an exemplary embodiment of the present invention. Figure 1 As shown, the following steps are included:

[0094] Step S1: regularly obtain monitoring data of the industrial park, and establish an initial industrial park coordinate system Oxy for the industrial park, wherein the origin of the initial industrial park coordinate system is the geographical center of the industrial park, the horizontal coordinate is parallel to the latitude, and the vertical coordinate is parallel to the longitude; when new monitoring data is obtained, proceed to step S2; otherwise, proceed to step S1;

[0095] Step S2: taking the new monitoring data as the current monitoring data, based on the current monitoring data, taking the monitoring site with the largest pollutant concentration as the coordinate origin, taking the wind direction as the positive direction of the x-axis, and performing coordinate transformation on the initial industrial park coordinate system Oxy to obtain the current industrial park coordinate system Ox'y';

[0096] Step S3: Based on the current industrial park coordinate system Ox'y', the industrial park is divided into grids, and the pollutant concentration corresponding to each grid of the industrial park is interpolated and calculated;

[0097] Step S4: Construct a first objective function, based on the pollutant concentrations corresponding to each grid in the industrial park, and calculate the parameter θi =(u i ,σ i ,ξ i ,w i ) to perform fitting, and obtain a first fitting result;

[0098] Among them, Θ i is a parameter set, i=1, 2, …, n, n is the number of columns in the grid, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The variance of i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid;

[0099] Step S5: Construct a second objective function, based on the pollutant concentrations corresponding to each grid in the industrial park and the first fitting result, and adjust the parameters Г=(Q, H0, tx v ) is fitted to obtain the second fitting result, where Г is the parameter set, Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of pollution source;

[0100] Step S6: Based on the second fitting result, determine the location of the air pollution source in the industrial park.

[0101] Since the emissions from industrial parks are divided into time periods and are continuous within the emission period, the method of the present invention is suitable for dealing with such situations in industrial parks.

[0102] The method of the present invention makes the following assumptions:

[0103] (1) The pollutant is a single point source and the leakage source strength is stable.

[0104] (2) The wind speed is constant.

[0105] (3) Diffusion is isotropic.

[0106] (4) Ignoring the terrain undulations, the ground can be represented by a plane.

[0107] (5) Pollutants will not penetrate into the ground.

[0108] The step S1, regularly obtaining monitoring data of the industrial park, includes:

[0109] Obtaining a topographic map of the industrial park; establishing an initial industrial park coordinate system Oxy for the industrial park;

[0110] The industrial park has S1 emission chimneys and L1 monitoring stations. The coordinates of the emission chimneys are (ChX s ,ChY s ), the coordinates of the monitoring stations and the pollutant concentrations are recorded as ,in, , ChX s ,ChY s are the horizontal and vertical coordinates of the sth emission chimney in the initial industrial park coordinate system Oxy, Respectively The horizontal and vertical coordinates of the monitoring stations and the concentration of pollutants monitored. The wind direction angle of the first monitoring station is The included angle in the positive longitude of the monitoring stations is β.

[0111] In the present invention, the current time data, wind direction angle data, wind speed data at a height of 10 meters, daytime radiation conditions, and nighttime cloud cover conditions are read. The daytime radiation conditions are divided into strong, medium, and weak, and the nighttime cloud cover conditions are divided into cloudy and cloudless.

[0112] The step S2 comprises:

[0113] The new monitoring data is used as the current monitoring data, and based on the current monitoring data, the coordinates (x p ,y p ) as the coordinate origin, and the wind direction as the positive direction of the x-axis;

[0114] The initial industrial park coordinate system Oxy is transformed to obtain the current industrial park coordinate system Ox'y'. The coordinate transformation formula is:

[0115]

[0116] Wherein, x and y are the coordinates of the initial industrial park coordinate system Oxy respectively;

[0117] In the current industrial park coordinate system Ox'y', the monitoring site coordinates and pollutant concentrations are recorded as .

[0118] The step S3 comprises:

[0119] calculate ;

[0120] The area [X1, X2] × [Y1, Y2] is divided into n × m uniform grids, and the coordinates of the grids are marked as (tx i ,ty j ), i=1,2,…,n, j=1,2,…,m, where:

[0121]

[0122] The coordinates calculated by inverse distance interpolation are (tx i ,ty j )'s corresponding pollutant concentration:

[0123]

[0124] Step S4: construct a first objective function, based on the pollutant concentrations corresponding to each grid in the industrial park, i =(u i ,σ i ,ξ i ,w i ) is fitted to obtain a first fitting result, including:

[0125] Step S41: Constructing the first fitting function ;

[0126] Where x=tx i , represents the straight line where the i-th column grid is located, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The corresponding variance, ξ i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid;

[0127] Step S42: Construct a first objective function, the first objective function is ,in,

[0128] ,

[0129] f1(θ i ),…,f m (Θ i ) are all intermediate vectors, c i1 ,…,c im They are calculated using the inverse distance interpolation method (tx i ,ty1),…,(tx i ,ty m ) point;

[0130] Step S43: iteratively solving the first objective function to obtain a first fitting result.

[0131] The step S43 comprises:

[0132] Step S431: Initialize θ i =(u i ,σ i ,ξ i ,w i )’s initial value Θ i (0) =(u i (0) ,σ i (0) ,ξ i (0) ,w i (0) ), and initialize J (0) , set the current iteration step k to 0, set the value range of the first threshold ε to 0<ε<1, the value range of the second threshold θ to θ>0, the value range of the third threshold γ1 and the fourth threshold γ2 to 0<γ1<1<γ2, the value range of the fifth threshold η1 and the sixth threshold η2 to 0<η1<η2<1, initialize θ (0) ,θ (0) is the initial value of parameter θ, J (0) is the initialized Jacobian matrix;

[0133] Step S432: If , then the iteration ends and Θ i (k) As the first fitting result, where Θ i (k) =(u i (k) ,σ i (k) ,ξ i (k) ,w i (k) ),u i (k) , σ i (k) ,ξ i (k) , w i (k) That is the u obtained by fitting i , σ i ,ξ i , w i ; Otherwise, let k be k plus 1, and go to step S433; where J (k) is the Jacobian matrix corresponding to the kth iteration;

[0134] Step S433: Calculation , where I is the identity matrix, θ (k)is the kth iteration value of parameter θ, d (k) (θ (k) ) is the calculated value corresponding to the kth iteration, ,

[0135] calculate

[0136] in, ,ρ is the intermediate variable,

[0137] ;

[0138] d (k) =d (k) (θ (k) );

[0139] Step S434: If , then let ;

[0140] like , then let θ (k+1) =θ (k) ;

[0141] like , then let θ (k+1) =γ1θ (k)

[0142] Step S435: Calculate Θ i (k+1) =θ i (k) +d (k) (θ (k) ), go to step S432.

[0143] The present invention uses the LM (Levenberg-Marquardt) algorithm to find the parameter to be fitted as Θ i =(u i ,σ i ,ξ i ,w i ). In the present invention, the first objective function is transformed into , the reasoning basis is: to minimize the sum of squares of the errors between the fitted value and the true value. In the present invention, F(θ i ) and F'(Θ i ) is continuously differentiable. Let J(Θ i ) is the Jacobian matrix, and its structure is f1(Θ i ),f2(Θ i ),…,f m (Θ i ) respectively for Θ i =(ui ,σ i ,ξ i ,w i ) Find the matrix formed by partial derivatives:

[0144]

[0145] Get parameter u i ,σ i ,ξ i ,w i The approximate solution of i At this point, the function reaches its maximum value, namely (tx i ,u i ) is the straight line x=tx i The maximum pollutant concentration point is , the process is repeated n times to obtain all parameters u i ,σ i ,ξ i ,w i Approximate solution for (i=1,2,…,n).

[0146] The step S5 comprises:

[0147] Step S51: Calculation , tu is the estimated value of the vertical coordinate of the location of the air pollution source in the industrial park; calculation , is the estimated value of the concentration of foreign pollutants;

[0148] Step S52: construct a second fitting function, wherein the second fitting function is C(x, Q, H0, tx v ):

[0149]

[0150] Among them, the parameters to be fitted are Г=(Q,H0,tx v ), Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of the pollution source; μ is the wind speed, σ y , σ z are the diffusion coefficients in the horizontal and vertical directions of the current industrial park coordinate system Ox'y', respectively, and the calculation formula is as follows:

[0151] , t y1 ,t y2 ,t z1 ,t z2is the Bridges diffusion coefficient, the atmospheric stability type is determined according to meteorological conditions, and z is the diffusion parameter in the vertical upward direction perpendicular to the current industrial park coordinate system Ox'y';

[0152] Step S53: Construct a second objective function, the second objective function is ,in,

[0153]

[0154]

[0155] Among them, σ yi Indicates that tx i Substitute σ y The calculated result, σ zi Indicates that tx i Substitute σ z The calculation result is as follows:

[0156] ;

[0157] Step S54: iteratively solving the second objective function to obtain a second fitting result.

[0158] The step S54 comprises:

[0159] Step S541: Initialize Г=(Q,H0,tx v )’s initial value Г (0) =(Q (0) ,H0 (0) ,tx v (0) ), and initialize J1 (0) , set the current iteration step k1 to 0, set the value range of the seventh threshold ε1 to 0<ε1<1, the value range of the eighth threshold θ1 to θ1>0, the value range of the ninth threshold γ3 and the tenth threshold γ4 to 0<γ3<1<γ4, the value range of the eleventh threshold η3 and the twelfth threshold η4 to 0<η3<η4<1; initialize θ1 so that θ1 (0) >0,θ1 (0) is the initial value of parameter θ1;

[0160] Step S542: If , then the iteration ends and Г (k1) As the second fitting result, Г (k1) =(Q (k1) ,H0 (k1) ,tx v (k1) ), tx v (k1)As the second fitting result corresponding to tx v ; Otherwise, let k1 be k1 plus 1, and go to step S543; where J1 (k1) is the corresponding Jacobian matrix of the k1th iteration;

[0161] Step S543: Calculate , where I1 is the identity matrix, θ1 (k1) is the k1th iteration value of parameter θ1, d (k1) (θ (k1) ) is the calculated value corresponding to the k1th iteration;

[0162] ;

[0163]

[0164] calculate

[0165] in,

[0166]

[0167] d (k1) =d (k1) (θ1 (k1) );

[0168] Step S544: If , then let θ1 (k1+1) =γ4θ1 (k1) ;

[0169] like , then let θ1 (k1+1) =θ1 (k) ;

[0170] like , then let θ1 (k1+1) =γ3θ1 (k1) ;

[0171] Step S545: Calculate Г (k1+1) =G (k1) +d (k1) (θ1 (k1) ), enter step S542.

[0172] The present invention uses the LM (Levenberg-Marquardt) algorithm to find the parameters to be fitted as Γ=(Q, H0, tx v ). In the present invention, the second objective function is transformed into , the reasoning is: minimize the sum of squares of the errors between the fitted value and the true value. In the present invention, F(Г) and F'(Г) are continuously differentiable. J(Г) is the Jacobian matrix, and the structure is f1(Г), f2(Г), ..., f n (Г) respectively for Г=(Q,H0,tx v ) Find the matrix formed by partial derivatives:

[0173]

[0174] Get the parameters Q, H0, tx v The approximate solution of .

[0175] In the present invention, the atmospheric stability is determined according to Table 1, and the Bridges diffusion coefficient is determined according to Table 2.

[0176]

[0177] In the present invention, the LM (Levenberg-Marquardt) algorithm is used to find the parameters to be fitted.

[0178] The present invention reads the monitoring data of the industrial park and uses the Gaussian plume diffusion model principle:

[0179]

[0180] Through two times of parameter fitting (the first parameter fitting is repeated n times), the specific coordinates of the pollution source emission chimney, the leakage source strength, the leakage source height, and the concentration of foreign pollutants are calculated.

[0181] The step S6 comprises:

[0182] Get the tx corresponding to the second fitting result v ,tu;

[0183] The coordinates (tx v ,tu) as the fitting position of the pollution source in the current industrial park coordinate system Ox'y', let x'=tx v , y'=tu, calculate the inverse coordinate transformation according to the following formula:

[0184]

[0185] Assign x and y to x respectively. v ,u;

[0186] Get the coordinates (tx v ,tu) in the initial industrial park coordinate system Oxy, denoted as (x v ,u);

[0187] Calculate (x v,u) Distance from each emission chimney:

[0188]

[0189] The emission chimney with the shortest distance is the location of the air pollution source in the industrial park.

[0190] Exemplary Devices

[0191] Figure 5 FIG. 1 is a schematic diagram of a device for tracing the source of air pollution in an industrial park provided by an exemplary embodiment of the present invention. Figure 5 As shown, this embodiment includes:

[0192] Data acquisition module: configured to periodically acquire monitoring data of the industrial park, establish an initial industrial park coordinate system Oxy for the industrial park, the origin of the initial industrial park coordinate system is the geographical center of the industrial park, the horizontal coordinate is parallel to the latitude, and the vertical coordinate is parallel to the longitude; when new monitoring data is acquired, the coordinate conversion module is triggered; otherwise, the data acquisition module is triggered;

[0193] Coordinate conversion module: configured to use the new monitoring data as current monitoring data, based on the current monitoring data, use the monitoring site with the largest pollutant concentration as the coordinate origin, use the wind direction as the positive direction of the x-axis, and perform coordinate transformation on the initial industrial park coordinate system Oxy to obtain the current industrial park coordinate system Ox'y';

[0194] A grid division module: configured to divide the industrial park into grids based on the current industrial park coordinate system Ox'y', and fit the pollutant concentration corresponding to each grid of the industrial park;

[0195] The first fitting module is configured to construct a first objective function, based on the pollutant concentrations corresponding to each grid of the industrial park, and to adjust the parameter θ i =(u i ,σ i ,ξ i ,w i ) to perform fitting, and obtain a first fitting result;

[0196] Among them, Θ i is a parameter set, i=1, 2, …, n, n is the number of columns in the grid, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The variance of i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid;

[0197] The second fitting module is configured to construct a second objective function, based on the pollutant concentrations corresponding to each grid of the industrial park and the first fitting result, to adjust the parameter Г=(Q,H0,tx v ) is fitted to obtain the second fitting result, where Г is the parameter set, Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of pollution source;

[0198] Pollution source determination module: configured to determine the location of the air pollution source in the industrial park based on the second fitting result.

[0199] Exemplary Electronic Devices

[0200] Figure 6 This is the structure of an electronic device 60 provided by an exemplary embodiment of the present invention. The electronic device may be any one or both of the first device and the second device, or a stand-alone device independent of them, and the stand-alone device may communicate with the first device and the second device to receive the collected input signals from them. Figure 6 FIG. 1 is a block diagram of an electronic device according to an embodiment of the present disclosure. Figure 6 As shown, the electronic device includes one or more processors 61 and a memory 62 .

[0201] The processor 61 may be a central processing unit (CPU) or other forms of processing units having data processing capabilities and / or instruction execution capabilities, and may control other components in the electronic device to perform desired functions.

[0202] The memory 62 may include one or more computer program products, which may include various forms of computer-readable storage media, such as volatile memory and / or non-volatile memory. The volatile memory may include, for example, random access memory (RAM) and / or cache memory (cache), etc. The non-volatile memory may include, for example, read-only memory (ROM), hard disk, flash memory, etc. One or more computer program instructions may be stored on the computer-readable storage medium, and the processor 61 may run the program instructions to implement the methods of the software programs of the various embodiments of the present disclosure described above and / or other desired functions. In one example, the electronic device may also include: an input device 63 and an output device 64, which are interconnected through a bus system and / or other forms of connection mechanisms (not shown). In addition, the input device 63 may also include, for example, a keyboard, a mouse, etc. The output device 64 may output various information to the outside. The output device 64 may include, for example, a display, a speaker, a printer, a communication network and a remote output device connected thereto, etc.

[0203] Of course, to simplify, Figure 6 Only some of the components related to the present disclosure in the electronic device are shown, and components such as a bus, an input / output interface, etc. are omitted. In addition, according to specific application situations, the electronic device may further include any other appropriate components.

[0204] Exemplary computer program products and computer-readable storage media

[0205] In addition to the above-mentioned methods and devices, an embodiment of the present disclosure may also be a computer program product, which includes computer program instructions, which, when executed by a processor, enable the processor to execute the steps of the method according to various embodiments of the present disclosure described in the above-mentioned "Exemplary Method" section of this specification.

[0206] The computer program product may be written in any combination of one or more programming languages ​​to write program code for performing the operations of the disclosed embodiments, including object-oriented programming languages ​​such as Java, C++, etc., and conventional procedural programming languages ​​such as "C" or similar programming languages. The program code may be executed entirely on the user computing device, partially on the user device, as a separate software package, partially on the user computing device and partially on a remote computing device, or entirely on a remote computing device or server.

[0207] In addition, an embodiment of the present disclosure may also be a computer-readable storage medium having computer program instructions stored thereon, which, when executed by a processor, enable the processor to execute the steps of the method according to various embodiments of the present disclosure described in the above “Exemplary Method” section of this specification.

[0208] The computer readable storage medium can adopt any combination of one or more readable media. The readable medium can be a readable signal medium or a readable storage medium. The readable storage medium can include, for example, but is not limited to, a system, device or device of electricity, magnetism, light, electromagnetic, infrared, or semiconductor, or any combination of the above. More specific examples (non-exhaustive list) of readable storage media include: an electrical connection with one or more wires, a portable disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the above.

[0209] The basic principles of the present disclosure are described above in conjunction with specific embodiments. However, it should be noted that the advantages, strengths, effects, etc. mentioned in the present disclosure are only examples and not limitations, and it cannot be considered that these advantages, strengths, effects, etc. are required by each embodiment of the present disclosure. In addition, the specific details disclosed above are only for the purpose of illustration and ease of understanding, and are not limitations. The above details do not limit the present disclosure to the necessity of adopting the above specific details to be implemented.

[0210] Each embodiment in this specification is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the system embodiment, since it basically corresponds to the method embodiment, the description is relatively simple, and the relevant parts can be referred to the partial description of the method embodiment.

[0211] The block diagrams of the devices, apparatuses, equipment, and systems involved in this disclosure are only illustrative examples and are not intended to require or imply that they must be connected, arranged, and configured in the manner shown in the block diagrams. As will be appreciated by those skilled in the art, these devices, apparatuses, equipment, and systems can be connected, arranged, and configured in any manner. Words such as "including," "comprising," "having," and the like are open words, referring to "including but not limited to," and can be used interchangeably therewith. The words "or" and "and" used herein refer to the words "and / or," and can be used interchangeably therewith, unless the context clearly indicates otherwise. The word "such as" used herein refers to the phrase "such as but not limited to," and can be used interchangeably therewith.

[0212] The method and apparatus of the present disclosure may be implemented in many ways. For example, the method and apparatus of the present disclosure may be implemented by software, hardware, firmware, or any combination of software, hardware, and firmware. The above order of steps for the method is for illustration only, and the steps of the method of the present disclosure are not limited to the order specifically described above, unless otherwise specifically stated. In addition, in some embodiments, the present disclosure may also be implemented as a program recorded in a recording medium, which includes machine-readable instructions for implementing the method according to the present disclosure. Therefore, the present disclosure also covers a recording medium storing a program for executing the method according to the present disclosure.

[0213] It should also be noted that in the apparatus, equipment and method of the present disclosure, each component or each step can be decomposed and / or recombined. These decompositions and / or recombinations should be regarded as equivalent schemes of the present disclosure. The above description of the disclosed aspects is provided to enable any technician in the field to make or use the present disclosure. Various modifications to these aspects are very obvious to those skilled in the art, and the general principles defined herein can be applied to other aspects without departing from the scope of the present disclosure. Therefore, the present disclosure is not intended to be limited to the aspects shown here, but to the widest scope consistent with the principles and novel features disclosed herein.

[0214] The above description has been given for the purpose of illustration and description. In addition, this description is not intended to limit the embodiments of the present disclosure to the forms disclosed herein. Although multiple example aspects and embodiments have been discussed above, those skilled in the art will recognize certain variations, modifications, changes, additions and sub-combinations thereof.

Claims

1. A method for tracing the source of air pollution in an industrial park, characterized in that: The method comprises: Step S1: regularly obtain monitoring data of the industrial park, and establish an initial industrial park coordinate system Oxy for the industrial park, wherein the origin of the initial industrial park coordinate system is the geographical center of the industrial park, the horizontal coordinate is parallel to the latitude, and the vertical coordinate is parallel to the longitude; when new monitoring data is obtained, proceed to step S2; otherwise, proceed to step S1; Step S2: taking the new monitoring data as the current monitoring data, based on the current monitoring data, taking the monitoring site with the largest pollutant concentration as the coordinate origin, taking the wind direction as the positive direction of the x-axis, and performing coordinate transformation on the initial industrial park coordinate system Oxy to obtain the current industrial park coordinate system Ox'y'; Step S3: Based on the current industrial park coordinate system Ox'y', the industrial park is divided into grids, and the pollutant concentration corresponding to each grid of the industrial park is interpolated and calculated; Step S4: Construct a first objective function, based on the pollutant concentrations corresponding to each grid in the industrial park, and calculate the parameter θ i =(u i ,σ i ,ξ i ,w i ) to perform fitting, and obtain a first fitting result; Among them, Θ i is a parameter set, i=1, 2, …, n, n is the number of columns in the grid, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The variance of i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid; Step S5: Construct a second objective function, based on the pollutant concentrations corresponding to each grid in the industrial park and the first fitting result, and adjust the parameters Г=(Q, H0, tx v ) is fitted to obtain the second fitting result, where Г is the parameter set, Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of pollution source; Step S6: Based on the second fitting result, determine the location of the air pollution source in the industrial park.

2. The method according to claim 1, characterized in that The step S2 comprises: The new monitoring data is used as the current monitoring data, and based on the current monitoring data, the coordinates (x p ,y p ) as the coordinate origin; The initial industrial park coordinate system Oxy is transformed to obtain the current industrial park coordinate system Ox'y'. The coordinate transformation formula is: ， Wherein, x and y are the coordinates of the initial industrial park coordinate system Oxy respectively; In the current industrial park coordinate system Ox'y', the monitoring site coordinates and pollutant concentrations are recorded as , β is the The wind direction angle of the first monitoring station is The included angle in the positive longitude of the monitoring stations.

3. The method according to claim 2, characterized in that The step S3 comprises: calculate ; The area [X1, X2] × [Y1, Y2] is divided into n × m uniform grids, and the coordinates of the grids are marked as (tx i ,ty j ), i=1,2,…,n, j=1,2,…,m, where: ， The coordinates calculated by inverse distance interpolation are (tx i ,ty j )'s corresponding pollutant concentration: 。 4. The method according to claim 3, characterized in that The step S4 comprises: Step S41: Constructing the first fitting function ; Where x=tx i , represents the straight line where the i-th column grid is located, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The corresponding variance, ξ i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid; Step S42: Construct a first objective function, the first objective function is ,in, ， ， f1(θ i ),…,f m (Θ i ) are all intermediate vectors, c i1 ,…,c im They are calculated using the inverse distance interpolation method (tx i ,ty1),…,(tx i ,ty m ) point; Step S43: iteratively solving the first objective function to obtain a first fitting result.

5. The method according to claim 4, characterized in that The step S43 comprises: Step S431: Initialize θ i =(u i ,σ i ,ξ i ,w i )’s initial value Θ i (0) =(u i (0) ,σ i (0) ,ξ i (0) ,w i (0) ), and initialize J (0) , set the current iteration step k to 0, set the value range of the first threshold ε to 0<ε<1, the value range of the second threshold θ to θ>0, the value range of the third threshold γ1 and the fourth threshold γ2 to 0<γ1<1<γ2, the value range of the fifth threshold η1 and the sixth threshold η2 to 0<η1<η2<1, initialize θ (0) ,θ (0) is the initial value of parameter θ, J (0) is the initialized Jacobian matrix; Step S432: If , then the iteration ends and Θ i (k) As the first fitting result, where Θ i (k) =(u i (k) ,σ i (k) ,ξ i (k) ,w i (k) ),u i (k) , σ i (k) ,ξ i (k) , w i (k) That is the u obtained by fitting i , σ i ,ξ i , w i ; Otherwise, let k be k plus 1, and go to step S433; where J (k) is the Jacobian matrix corresponding to the kth iteration; Step S433: Calculation , where I is the identity matrix, θ (k) is the kth iteration value of parameter θ, d (k) (θ (k) ) is the calculated value corresponding to the kth iteration, , , calculate , in, , ρ is the intermediate variable, ； d (k) =d (k) (i (k) ); Step S434: If , then let ; like , then let θ (k+1) =θ (k) ; like , then let θ (k+1) =γ1θ (k) Step S435: Calculate Θ i (k+1) =θ i (k) +d (k) (θ (k) ), go to step S432.

6. The method according to claim 4, characterized in that The step S5 comprises: Step S51: Calculation , tu is the estimated value of the vertical coordinate of the location of the air pollution source in the industrial park; calculation , is the estimated value of the concentration of foreign pollutants; Step S52: construct a second fitting function, wherein the second fitting function is C(x, Q, H0, tx v ): ， Among them, the parameters to be fitted are Г=(Q,H0,tx v ), Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of the pollution source; μ is the wind speed, σ y , σ z are the diffusion coefficients in the horizontal and vertical directions of the current industrial park coordinate system Ox'y', respectively, and the calculation formula is as follows: ， t y1 ,t y2 ,t z1 ,t z2 is the Bridges diffusion coefficient, the atmospheric stability type is determined according to meteorological conditions, and z is the diffusion parameter in the vertical upward direction perpendicular to the current industrial park coordinate system Ox'y'; Step S53: Construct a second objective function, the second objective function is ,in, ， ， Among them, σ yi Indicates that tx i Substitute σ y The calculated result, σ zi Indicates that tx i Substitute σ z The calculation result is as follows: ； Step S54: iteratively solving the second objective function to obtain a second fitting result.

7. The method according to claim 6, characterized in that The step S54 comprises: Step S541: Initialize Г=(Q,H0,tx v )’s initial value Г (0) =(Q (0) ,H0 (0) ,tx v (0) ), and initialize J1 (0) , set the current iteration step k1 to 0, set the value range of the seventh threshold ε1 to 0<ε1<1, the value range of the eighth threshold θ1 to θ1>0, the value range of the ninth threshold γ3 and the tenth threshold γ4 to 0<γ3<1<γ4, the value range of the eleventh threshold η3 and the twelfth threshold η4 to 0<η3<η4<1; initialize θ1 so that θ1 (0) >0,θ1 (0) is the initial value of parameter θ1; Step S542: If , then the iteration ends and Г (k1) As the second fitting result, Г (k1) =(Q (k1) ,H0 (k1) ,tx v (k1) ), tx v (k1) As the second fitting result corresponding to tx v ; Otherwise, let k1 be k1 plus 1, and go to step S543; where J1 (k1) is the corresponding Jacobian matrix of the k1th iteration; Step S543: Calculate , where I1 is the identity matrix, θ1 (k1) is the k1th iteration value of parameter θ1, d (k1) (θ (k1) ) is the calculated value corresponding to the k1th iteration; ； ， calculate , in, , ， d (k1) =d (k1) (θ1 (k1) ); Step S544: If , then let θ1 (k1+1) =γ4θ1 (k1) ; like , then let θ1 (k1+1) =θ1 (k) ; like , then let θ1 (k1+1) =γ3θ1 (k1) ; Step S545: Calculate Г (k1+1) =G (k1) +d (k1) (θ1 (k1) ), enter step S542.

8. A device for tracing the source of air pollution in an industrial park, characterized in that: The device comprises: Data acquisition module: configured to periodically acquire monitoring data of the industrial park, establish an initial industrial park coordinate system Oxy for the industrial park, the origin of the initial industrial park coordinate system is the geographical center of the industrial park, the horizontal coordinate is parallel to the latitude, and the vertical coordinate is parallel to the longitude; when new monitoring data is acquired, the coordinate conversion module is triggered; otherwise, the data acquisition module is triggered; Coordinate conversion module: configured to use the new monitoring data as current monitoring data, based on the current monitoring data, use the monitoring site with the largest pollutant concentration as the coordinate origin, use the wind direction as the positive direction of the x-axis, and perform coordinate transformation on the initial industrial park coordinate system Oxy to obtain the current industrial park coordinate system Ox'y'; A grid division module: configured to divide the industrial park into grids based on the current industrial park coordinate system Ox'y', and interpolate and calculate the pollutant concentration corresponding to each grid of the industrial park; The first fitting module is configured to construct a first objective function, based on the pollutant concentrations corresponding to each grid of the industrial park, and to adjust the parameter θ i =(u i ,σ i ,ξ i ,w i ) to perform fitting, and obtain a first fitting result; Among them, Θ i is a parameter set, i=1, 2, …, n, n is the number of columns in the grid, u i is the estimate of the ordinate corresponding to the maximum value of the pollutant concentration in the pollutant concentration corresponding to each grid in the i-th column, σ i for u i The variance of i is the parameter adjustment coefficient of the i-th column grid, w i is the estimated value of the external pollutant concentration in the i-th grid; The second fitting module is configured to construct a second objective function, based on the pollutant concentrations corresponding to each grid of the industrial park and the first fitting result, to adjust the parameter Г=(Q,H0,tx v ) is fitted to obtain the second fitting result, where Г is the parameter set, Q is the leakage source strength, H0 is the leakage source height, tx v is the horizontal coordinate of pollution source; Pollution source determination module: configured to determine the location of the air pollution source in the industrial park based on the second fitting result.

9. A computer-readable storage medium, characterized in that: The storage medium stores a plurality of instructions; the plurality of instructions are used for a processor to load and execute the method as described in any one of claims 1 to 7.

10. An electronic device, characterized in that: The electronic device comprises: a processor for executing a plurality of instructions; and a memory for storing a plurality of instructions; wherein the plurality of instructions are used to be stored in the memory and loaded and executed by the processor according to any one of the methods described in claims 1-7.

Images