Method for determining the spatial coordinates of a source that causes a dispersion phenomenon

Abstract

Method for determining coordinates in space of a source at the origin of a dispersion phenomenon in a domain, the dispersion phenomenon being approximated by a system of partial differential equations (20), the method comprising computer-implemented steps of obtaining measurements (60) each quantifying at a point in the space of the domain and at a given instant a state variable modeled by the system of partial differential equations at said point and at said instant and iterative minimization (70) of the differences between the measurements and the calculations of the system of equations to locate a position (80) of said source.The system of equations is expressed (50) by finite differences supported by radial basis functions, the system of equations being solved by least squares, using a least squares solver, the system of equations having undergone an adjoint quadratic solution by choosing radial basis functions from among polyharmonic spline functions, multi-quadric inverse functions, or Gaussian functions. Figure for the abbreviation: FIGURE 3.

Description

Title of the invention: Method for determining coordinates in space of a source causing a dispersion phenomenon Field of the invention

[0001] The invention relates to the field of modeling the dispersion of atmospheric pollution, in particular in an urban environment.

[0002] We are interested in the identification (location, shape estimation and amplitude estimation) of sources of atmospheric pollution in an urban environment, in particular the urban canopy, by having a set of localized measurements.

[0003] Thus, the problem addressed is the identification of pollution sources in urban environments by a network of low-cost sensors. The sources emit pollutants (gases, particles) which are dispersed by the wind. The complex geometry of urban buildings generates a wind field which is not homogeneous on the scale of a neighborhood.

[0004] Based on a choice of modeling and a set of available noisy measurements (sensor network or mobile sensors), the aim is to solve an inverse problem to estimate the position and characteristics of the sources. It is a problem of robust fusion of spatio-temporal data in the presence of spatial constraints and measurement uncertainties.

[0005] Known reliable methods are methods based on physical models where one searches for one or more unknown input terms applied to partial differential equations (PDEs) by minimizing the error between the measurements and the output of a mathematical model representing the physical system, the latter therefore being sensitive to unknown inputs which must be estimated.

[0006] Regarding the identification of pollution sources, several methods have been applied in the state of the art such as Kalman filtering and its variants (Defforge, 2019; Mons et al., 2017), the calculation of variations (Kumar et al., 2015) sometimes associated with the use of reduced models (Khodayi-mehr, 2019) or even particle filters (Septier et al., 2020).

[0007] These methods have the disadvantage of being expensive in terms of calculation, or in the context of using reduced models that are less expensive in terms of calculation, of requiring the creation of a learning database for each dispersion scenario.

[0008] Furthermore, atmospheric pollution is a three-dimensional phenomenon which needs to be modeled on fine and complex calculation spaces.

[0009] A specific tool: radial basis functions with finite differences in least squares

[0010] The source term estimation approach employed in the present invention uses LS-RBF-FD (Least-Squares Radial Basis Functions Finite Differences) functions to create an accurate numerical pollutant dispersion model. LS-RBF-FD type models have recently appeared in the literature (Tominec et al., 2021). Similar to RBF-FD which have been known for about 10 years (Flyer et al., 2012) they share their speed but have increased robustness in the calculation of solutions, particularly for the calculation of three-dimensional field variables. Air pollution modeling

[0011] LS-RBF-FDs are used in one embodiment on a pollutant transport equation called advection-diffusion and noted ADEDP. This equation models the transport of pollutants considered without reaction term and governs a scalar concentration field u: xeQ = Q°u3Q, t êThu(x, t). This equation is written,

[0012] [Eq. 1] of . tt du , tt du , t T du  + + + = ) + s, dx \ A dx / dy \ - dy ) dz \ 'dz / xeQ° 3.11 ■ n = 0 xgTn u(x, t) -0 xgTd u(x, t = 0) = u0 xg£Ï with te T the time, x GQC R3 the position, Ub Kj the components of the velocities and diffusivities of the mean wind and diffusion field,11 the pollutant concentration,s the pollution source term.

[0013] it represents the spatial domain where the advection-diffusion differential equation is applied while 3Q = FD UF^ represents the edges of the domain where boundary conditions are applied. We are interested in a situation where FD is the edge domain where Dirichlet conditions apply and FN zero flux conditions.

[0014] The meteorology, namely, the wind U and diffusion K fields, are obtained by a third-party model.

[0015] Advection / diffusion equations are difficult to model in an urban environment with spatial constraints and are computationally expensive. It is necessary to propose a solution to make the calculation accessible.

[0016] Writing the advection / diffusion equation in steady state, in the form of differential operators

[0017] The ADEDP is solved in a stationary regime. In addition, the physical model can be reformulated in a condensed manner by expressing the derivative operators applied on c as a generic differential operator L on Q° and B on the conditions limits in 9Q = Fp U such that,

[0018] Lc(x)=s(x), xeQ°; Bc(x)=g(x), xedn.

[0019] In ADEDP, s is the source term and g = 0.

[0020] RBFs are known to satisfactorily interpolate spaces (2D, 3D, 4D in particular).

[0021] Any function f : yg Q f^yj ef(Q^ c Rd (d being the dimension of the space) can be reconstructed by some basic functions (RBF) 4» centered on nodes fxt Nl = Xc Q- ^we denote f 'interpolating function of f then, 1 1J i= î

[0022] [Eq. 2] f(y) lly- x ill )> VyeQ 0

[0023] Where || ■ || denotes the Euclidean norm. It is known that | is a satisfactory approximation of f proving under the condition of regular spacing of the base centers The interpolation problem then admits a unique solution of coefficient vector a for different bases [(Fomberg and Flyer, 2015) Chapter 3, p. 39-47] in particular for the PHS bases (polyharmonic splines or polyharmonic splines) [(Barnett, 2015) p. 13 theorem 2] used here. The PHS have the expression,

[0024] [Eq. 3] ^y) - r2k-i? k > 13 with r - || y - Xj ||.

[0025] Inverse problem in the sense of estimating source terms of partial differential equations

[0026] Solving the inverse problem is considered as executing the workflow with model-based methods.

[0027] A vector of unknowns 0 is the variable to be estimated, and is initialized to the value 00. For example in the advection-diffusion equation in (1)0 can correspond to the entire source vector for example 0 = s but also if the source is unique 6 can correspond to a parameterization which is expressed with an amplitude and a position for example 0 = (as, xs).

[0028] The Um measurements sample the field variable u, for example the concentration field of the advection-diffusion equation. They are generally noisy, and the algorithm implementing the inverse method can be provided with a noise model if the noise of the sensors is characterized.

[0029] The model (dispersion model) of physical constraint m applies to the unknown inputs of model 0 so that ue = m(9). Thus the calculated field variable ue is of the same nature as um. We apply a measurement model h to match the measurements Uni, in particular their positions, with llï(0) to calculate a deviation (generally quadratic) of type r2 — [um-h( m ( 0) )]2 9U' Allows the expression of an error term.

[0030] Variational or Bayesian methods optimize this error term to return an estimate g, possibly with uncertainty ranges. The estimation generally requires several calls to the physical model m, which is computationally expensive.

[0031] The output of the workflow is estimated parameters and uncertainty if the method is Bayesian.

[0032] The inverse problem addressed is the reconstruction of the discretized source term s. To identify a source term, whether parameterized or not, we speak of Source term estimation. State of the art

[0033] Radial basis functions, their finite difference variants and their finite difference and least squares variants LS-RBF-FD are known. LS-RBF-FD are used here for air pollution modeling and for solving the inverse problem of pollution source estimation in combination with model reduction.

[0034] The article by (Khodayi-mehr, 2019) uses reduced models for source identification. Model reduction is achieved by principal component analysis (PCA). To build these models, the authors created several stationary dispersion scenarios by varying only the position of a pollution source. They used an adjoint variational method for their inverse problem, which makes multiple calls to the reduced model learned by PCA on two stationary two-dimensional case studies. But height transport is not negligible. It is then preferable to use three-dimensional transport models. In addition, PCA requires training data. If the physical domain changes, or if the weather or the shape of the pollutant sources changes, a database must be recalculated.

[0035] The article (Liu et al., 2019) focuses on air pollution modeling based on RBF-FD, but not resolved either three-dimensionally or in least squares. The authors use a two-dimensional RBF-FD unsteady physical pollutant transport model. They calculated the weather and reused the mesh as calculation points for their model. The calculation domain includes rectangular obstacles. A complex urban geometry would require treating the interior courtyards of buildings, which would have to be closed beforehand at the risk of creating model indeterminacies. In addition, height transport is not negligible in dispersion.

[0036] The article by (Tominec et al., 2021) uses the LS-RBF-FD technique. The LS- RBF-FD is a generalization of RBF-FD. It is based on two sets of calculation points X and Y at Nx and Ny points.

[0037] Y is informed by physics and contains the boundary condition constraints of geometric shapes. X is less dense than Y and covers the physical domain without

[0038]

[0039]

[0040]

[0041]

[0042] need to be confined to the edges of the domain (X may have points that extend beyond the domain, or only inside the domain). Thus, in the manner of RBF-FD, it is possible to construct a differentiation matrix D which discretely models by RBF-FD the operators L and B defined in the interior of Q, Q and its edges d£l so that the solution of the system comparable to (2) is expressed by: 2 ü a = argmin || Dv - f || With (j £ rNx, and D a differentiation matrix Ny x Nx with, D U f= v k / k 0 if y, e 3 / 2 VK Where L and B are the differentiation matrices obtained by RBF-FD for the spatial operators L and B. And a point of Y.

[0043] The article by Tominec et al., 2021 uses uniform diffusion diffusive equations with mixed boundary conditions (Dirichlet, Neumann). No three-dimensional solutions are presented or discussed. Air pollution modeling is not mentioned. Inverse problems are not addressed.

[0044] The principle of RBF and in particular RBF-FD is presented among others in (Bamett, 2015) p. 2-13.

[0045] Creation of partial differential equation models supported by RBF

[0046] RBFs have spatial interpolation properties. If we consider a finite set YGQ which partitions Q and allows us to represent the variations of in equation (3) with ( i Ny the elements of Y, and we further consider ye Y= y. HjJpl two sets of computational nodes defined in the interior of the domain Q° and on its edge The first set being a set X of spanning nodes Q, with Nx = card(X) - these nodes are used later as RBF centers to perform local interpolation approximations in the manner of equation (2). The distribution of the nodes uses an energy minimization algorithm that uses the

[0047]

[0048]

[0049]

[0050]

[0051]

[0052]

[0053]

[0054]

[0055]

[0056] Halton pseudo-random sequence to cover the spatial domain satisfactorily. A set Y of evaluation nodes equipped with two subsets of nodes Yu-et corresponding respectively to the nodes inside the domains, i.e. in Y^ and those present on the edge, Y5£i. Y is generated in a similar way to X. We note = card(Y£ï) and Ny^ = card(Y3 J> with Ny = card(Y) = Nytf + Ny— Methods based on RBF-FD: The idea in this method is to approximate the constraints imposed on a variable u governed by a partial differential equation at a point in space y = yk G Y. Base centers / PDE points (X) neighboring the evaluation point are used to express these constraints. This set of neighboring points in X forms a set of points centered on yk G Y. This set is called Sy CX and called a stencil. Let xL G Syk, it is possible to locally interpolate a variable u GQ such that V Xyk G Syk c X, [Eq ' 7| uJ'xY) =E” =1 Wi(|>( Jx^-xJI ) \p.(x y k). With w> and \ the weighting coefficients of the approximation and the Pi of the monomials used to regularize the approximation to impose In , . ( We denote the local interpolating function of u. We note n- fv 1 the vector of u for all xyt e Sv • Local interpolation (7) UîW a\ i / i A can be vectorized such that, Where d* is the RBF collocation matrix of the Syk stencil such that g |x.|n — $y [Eq. 9] ¢( Il Xn-Xj || ) <l>( || Xn-Xn [Fig.l] An example of a stencil is seen in Figure 1. It is composed of X points, centered around yk and Y among other points. The vectors w and X are obtained by inverting the matrix a / v \, [Eq. 10] zœ\ . si / ujv U

[0057]

[0058]

[0059]

[0060]

[0061]

[0062]

[0063]

[0064]

[0065]

[0066]

[0067]

[0068]

[0069]

[0070]

[0071] In a similar way, one can locally interpolate a differential operator L on u at a point yk GY such that, [Eq ' 11] Read a (y) = There is - x ! He ) + ^i X i L Pi( y ) Equation 11 can be reformulated in the manner of equation 8 in [Eq - 12] Read a (yJ = (lo) T ( y k ) w + (LP^yj By introducing (10) into (12) we obtain, W 'UU((<1 (LPr(y k ))A-^^^^ \ / \ 0 / And by setting, [Eq. 14] / _ i T <lt\ \\t / t \ \ A-1, on obtient ( w •>' ) = ((L <t)(yJ) , (LPj(yJ)A’ after removing the polynomial terms that cancel, [Eq. 15] LUa^yJ = w^Ua^yJ Which is equivalent to Read a (y J = ( Ua ( ) ' ' • ' ' Ua ( X " k ) ) T ' X G Sy k C X To express the finite differences applied to the operator Es vectors of weights wL and \L are found by solving the system [Eq. 16] üU.L1 <P P W L Lpt oJLxL. ■(M>Xy=y k )' .(Lp)(y = y k ). It is thus possible to approximate each differential operator Lu^ yj V y GY as a weighted sum of n neighboring values ​​of u defined in X. It is possible to carry out this approximation on other differential operators such as that of boundary conditions B. Thus it is possible to create a discrete differentiation matrix D with dimensions Ny x Nx which for any yk GY expresses the differentiation weights which apply to it. The weights are stored row by row such that V ye Y, [Eq. 17] [0, ..., wp, ..., w^, ...0]Yy k eQ° [0, wp, ..., w§, ...0] if The 0s being added in the weight vector of D for the nodes of X not part of the stencil Sy . General definition of invention The invention is concerned with modeling, in particular inverse modeling of the urban air pollution and uses LS-RBF-FD modeling. LS-RBF-FD is used to express the direct air pollution model and a least squares solver to obtain the solution, especially for advection-diffusion type transport equations, where the RBF-FD method only does collocation.

[0072] The invention provides model reduction without a learning method, by choosing a small solution space X, which speeds up the calculation and is useful for inversion.

[0073] An advantage of this method is the reduction in the size of the unknown constituted by the source vector: the problem is thus better posed.

[0074] Generally speaking, RBFs allow the projecting of the sought variables onto small-dimensional spaces.

[0075] The invention also consists of an original and efficient adjoint quadratic resolution. It is implemented efficiently on a simple case study.

[0076] The invention relates to a novel modeling of atmospheric pollutant transport using the LS-RBF-FD method. This method is fast and has an inherent capacity to operate model reduction, which opens the way to tracing back to the source.

[0077] Thus, to improve existing source search techniques, a method is proposed for determining coordinates in space of a source at the origin of a dispersion phenomenon in a domain, the dispersion phenomenon being approximated by a system of partial differential equations, the method comprising steps, implemented by computer,

[0078] - obtaining measurements each quantifying at a point in the domain space and at at a given instant a state variable modeled by the system of partial differential equations at said point and at said instant and

[0079] - of iterative minimization of the differences between the measurements and the calculations of the system of equations to locate a position of said source.

[0080] The domain is typically a part of an urban canopy, i.e. a district of a city or a peri-urban area, whose atmosphere is modeled, up to approximately 50 m above the ground, for example. The diffusion source may be a stationary car emitting exhaust gases, for example.

[0081] The method is particular because the system of equations is expressed by finite differences supported by radial basis functions, the system of equations being solved in least squares, using a least squares solver, the system of equations having been the subject of an adjoint quadratic resolution by choosing radial basis functions from among polyharmonic spline functions, multi-quadric inverse functions, or Gaussian functions, with regularization polynomials of order 1 or higher.

[0082] We take advantage of the fact that LS-RBF-FDs have an oversampled discretization and allow a resolution of a least squares problem.

[0083] The invention is based on the distinction of two sets of calculation points X, Y, one of which has a cardinality greater than the other by a factor of at least 3.

[0084] The method used guarantees numerical stability on the processor that carries out the calculations.

[0085] Thanks to this method based on LS-RBF-FD, we can take a not very fine approximation space (X), while obtaining quality solutions. This choice makes it possible to reduce the size of the unknown (the values ​​of the source vector sought) and to accelerate the calculations.

[0086] It skillfully uses a phenomenon of stress relaxation, and has more evaluation points than unknowns. We are looking for a source term 5 which minimizes the deviation from measurements when we apply the physical model to the source term.

[0087] According to advantageous and optional characteristics

[0088] - the steps of obtaining localized measurements and minimizing the differences between measurements and calculations can be carried out in a spatial domain described with Neumann and Dirichlet boundary conditions.

[0089] - the domain can be characterized, besides by its limits, also by data of field within the domain.

[0090] - the dispersion phenomenon may be a dispersion of pollutants and that said transport is modeled by an advection-diffusion equation.

[0091] - the dispersion phenomenon can be a diffusion of heat from a point hot in an energy system.

[0092] - the dispersion phenomenon can be a sound diffusion in an environment complex.

[0093] - the least squares solver can be L-BFGS-B or LSQR.

[0094] - the polyharmonic spline function used can be PHS5, and beyond this specific function we favor polyharmonic spline functions of order at least 3, and preferably odd order. We can use other RBFs such as multiquadric inverse functions or Gaussian functions.

[0095] - the iterative minimization of the differences between the measurements and the calculations of the system of equations can be carried out to also determine an intensity (80) of said source. List of figures:

[0096] The invention will be better understood and other advantages will appear on reading the description which follows, given without limitation and thanks to the appended figures among which:

[0097] [Fig.l] is the representation of a stencil, or set of neighboring points in

[0098]

[0099]

[0100]

[0101]

[0102]

[0103]

[0104]

[0105]

[0106]

[0107]

[0108]

[0109]

[0110] [YES]

[0112]

[0113]

[0114] which approximate the constraints imposed at said point on a variable governed by a partial differential equation. [Fig.2] is the representation of the set of points in the domain, quite dense, for which the constraints of the boundary conditions and the applicable physical rules are given in the model. [Fig. 3] is a representation of an embodiment of a method according to the invention. [Fig.4] is the representation of results (left part) of the method, implemented for a particular case (represented right part). [Fig.5] is the representation of other results for the same particular case. Figures 6 and 7 are the representation of other results for the same particular case with prior use of simulated results. [Fig.8] is a representation of results obtained as part of an extension of the method. The particularity of the LS-RBF-FD and its advantages: LS-RBF-FD poses two distinct sets of calculations such that YX and card(Y) / card(X) > 3. The value of 3 or more is a threshold that has proven reasonable. The technique then consists of solving a least squares problem in the case of the stationary ADEDP presented, ^■'^u^argminllDu-fll 2 . «"•(«. o) el^xRMeR* le ü agglomerated source term defined in the coordinates of Y including the zero boundary conditions of Neumann and Dirichlet and the source term such that yJ«y k eY ff cÆ 0 if yk e Y3q The least squares formulation of equation 18 has several advantages: - It helps to relax the equality constraints on boundary conditions of type u = d(y) on FD and V u • n = t(y) on FN for the Neumann and Dirichlet conditions respectively. - The stability of one is encouraged by the fact that the problem in equation 18 is overdetermined, with more evaluation points in Y than unknowns in X. - The method is robust even on sparse computing grids. - The differentiation matrix D is very sparse (it contains a lot of 0s) and therefore the resolution is done in a short time. - It is possible to choose a very small solution space X, 5 times less dense than Y for example to project the model onto a space of particular small dimension-

[0115] fully suited to inverse problems. By the definition of D with D W= y ) if yg Y q = g 12 ok / •'k s (y k ) ^y k eY aQ , we can extract from D two sub-matrices, Dn« and D3q which stack respectively the LW for all yk G Yiy and gfy ] for all yk GY^. \ k /

[0116] Formulation (18) admits a unique solution when D is of full rank with [Eq. 20]Ua^D+f, with p+ _ the left pseudo-inverse of D. In practice the The problem in equation 18 is solved by a least-squares solver. The solvers L-BFGS-B (variant B of the Limited-memory BFGS algorithm - for Broyden-Fletcher-Goldfarb-Shanno) and LSQR (least-squares algorithm) were tested. The fastest one turned out to be LSQR. However, from equation 20 we can construct an interesting optimality condition for the inverse problem,

[0117] [Eq. 21] Eu,, -D^Yn-) = (1 with E =

[0118] [Fig.2] Figure 2 is a three-dimensional view of a Y domain generated with Ny = 15000. We see the interior knots in Q°, as well as Neuman knots and Dirichlet knots.

[0119] Thus, a direct Eulerian modeling of pollution has been presented which allows to calculate a dispersion of pollutants by giving as input a topology and heterogeneous meteorological data. The calculation is carried out by LS-RBF-FD. The solution to the problem of stationary equation 1 is calculated efficiently by a least squares solver. Inverse Modeling, an adjoint approach

[0120] The inverse problem consists of simultaneously estimating the source term § e RNi and the pollutant concentration ue RNx using a limited number of measurements. Ÿ RN' with Ns, the number of measurements located in Q°. If a source term s is known in X it is possible to estimate the source at any point zg Yq« using the interpolation formula of equation 16 evaluated with neutral differential operators (without application of derivative). Thus we can construct an interpolation operator Ws which makes it possible to approximate in X a source defined in YQ° such that, 8^) = ^8..

[0121] The measurement operator denoted C is similarly defined as an interpolation operator from the field space X to the measurement space such that

[0122] [Eq. 22] ÿ = Cû, with u, the model output, ÿ the model output brought back into the space of measurements located in Q°.

[0123] For the purpose of estimating the source term, the least squares problem is solved next - this is the regularized least squares error minimization between the output of the LS-RBF-FD model and the measurements -

[0124] [Eq. 23] • t / \ _ j!| En IP aA Ile l|2 mmj ( s, u ) - 2 || Lu - y mes || r , + ô IP II K , under Eû - = 0,

[0125] where R and B are measurement error and regularization covariance matrices respectively.

[0126] The optimal solution (s* if ) of this source term estimation problem is the solution of the linear system,

[0127] [Eq. 24] / KE \ / - ' / CTR4V \ I ] { u ] — J mes [F „0 \xJ -

[0128] where K = cVt Q = DLWsBTjD0o. Moreover s* = BVjDfVX* where À^st is the adjoint of û”. Considering that the matrix E is positive definite, a* = [E + QE*Kf QB^R 1 ^ et i' est donnée P" À* = - E 'C T R ï Cû - J

[0129] To convince ourselves of this, we introduce the Lagrangian of equation 23., read, a 4 = 4 II Cn - v || 2 , + 4 h || 2 ,+x T (In - d / W La ""æ de La s ran ® e \ / 2 H - my 1st r 1 2 HH g! \ Li J allows us to deduce the necessary conditions of optimality with:

[0130] = )+E2 = 0

[0131] [Eq.26] V = Q o = 0,

[0132] [Eq. 27] v £ = Q o Eü - DL^Bs = 0, AT 7

[0133] Where Equation 25 is the adjoint equation of the model. If E is positive definite, and thus invertible, \ can be explicitly calculated as _ E^C^ ÏCU*- V ï VV 7 mes / equation 26 we obtain — B'PJDqA •

[0134] And in equation 27 by introducing the expression of §*, we obtain Eû* _ Qp — (). Finally, ÿ-s _ [g + qe 'K^QE ^^R^ÿ by introducing the expression of E 4- QE! K is always invertible if E is positive definite because the symmetric matrix QE !K is at least positive semi-definite.

[0135] The invertibility of the matrix E is guaranteed if and only if the differentiation matrix D in equation 18 is of full rank. This condition can be satisfied by choosing a suitable cardinality of X and Y. The linear system of equation 24 of size 2NX X 2NX is solved without explicit matrix inversion to speed up calculations on large domains, notably by using a least squares LSQR solver.

[0136] [Fig.3] Thus the principle of the invention is represented in [Fig.3]. The data processor receives geometric information 10, a model of continuous partial differential equations 20, and input field data 30 - which are given for the interior of the domain.

[0137] In the main embodiment, this input field data is meteorology data, but in an embodiment where one would be interested in the diffusion of heat in a material, the density in the material, the electric charge, or even the salinity.

[0138] The geometric information 10, the continuous partial differential equation model 20, and the input field data 30 are raw input data of the model.

[0139] The geometric information 10 is processed so as to provide boundary conditions to the domain, during a step 40. Processed input data is thus obtained, by the creation of a calculation domain, including the boundary conditions.

[0140] During a step 50, a discretization of the partial differential equation model 20 is carried out, to discretize it by the radial basis functions with finite differences LSRBF-FD. This is carried out taking into account both the computational domain established during step 40, and the input field data 30. Discretized physical operators are obtained.

[0141] The method naturally includes obtaining measurements 60.

[0142] Then, on the basis of the measurements 60 and the discretized physical operators obtained during the step 40, the inverse model 70 is implemented, including an iterative resolution based on the explanation, in the processor and the associated memory, of the adjoint equation of the discretized operators by matrix inversion, the radial basis functions having been chosen for the purpose of making this inversion possible. The inverse model thus comprises a resolution of a least squares problem for the purposes of the iterative resolution: a value representative of an error between the estimates by the model and the measurements is minimized.

[0143] At the end of the inverse model 70, we obtain the coordinates - that is to say the positions - of the sources (or of the source, in the case of a single source) as well as optionally the estimated emission speeds. We also optionally obtain the intensities (often concentrations) at the different points of the domain. These are the results 80. Implementation example

[0144] [Fig.4] The three-dimensional physical domain consists of a domain comprising three buildings with a height of 10, 15 and 12 meters respectively. The domain dimensions are Lx — 205 mx Ly =120 mx Lz= 105 m. The three-dimensional domain is represented in the right part of Figure 4. Figure 2 which has been mentioned above is a representation of the distribution of the evaluation nodes Y.

[0145] We use the PHS5 radial basis functions, according to the equation

[0146] [Eq. 28] (|^y) - r5 with r= j|y-xj|.

[0147] Second-order regularization polynomials are used. The direct LS-RBF-FD model of equation 21 is used with 15000 evaluation nodes to ensure good resolution. The set X has much fewer nodes with Nx = 2060 which naturally leads to a model reduction that limits the size of the vectors of unknowns s and a. Neumann boundary conditions of zero flux are assigned to the ground and the edges of the buildings as well as to the ceiling of the simulation domain while Dirichlet conditions are imposed on the other edges of the domain. The stencil size is 90, the maximum order of regularization of the polynomials in equation 7 reaches 2. The matrices R and B are chosen as diagonals with R — rÏNs and B = blyjx where In. denotes an identity matrix of size N. X N. and where 1 and b are diagonal regularization coefficients. In the examples presented, r = 1, b = 0.15.

[0148] [Fig.4] [Fig.4] shows, on the left, the pollutant concentration field (shade) and the wind field (arrows) at 2 m height above the ground. On the right, the area studied is shown, comprising three buildings of different sizes. The source is at a height of about two meters.

[0149] Meteorological data are used as input information in the calculations. The mean wind fields are provided by the PMSS (Parallel-Micro-SWIFT-SPRAY) software, where its quasi-stationary Eulerian Micro-SWIFT-M wind field model is used, which solves the Navier-Stokes equations on the Reynolds average principle to calculate a mean wind field U and a diffusion K which reflects the turbulent fluctuations of winds that can be assimilated over the long term to a diffusive phenomenon. Meteorological data calculated by models other than Micro-SWIFT-M could in principle be suitable for the LS-RBF-FD model presented above.

[0150] The Ÿmcs measurements are simulated by the LS-RBF-FD “forward” model on a fine computational domain and with the addition of noise such that

[0151] [Eq. 29] y — yû + gobs with, êobs ~ N( 0, (J2) °ù ÿ° is the reference pollutant concentration returned by a fine LS-RBF-FD model with Ny = 50000, Nx = 15648 (see Figure 3, left) and ëobs, a Gaussian noise of variance <72, with a = 0.002. The true source term used to generate the reference concentration field is a Gaussian defined in Q°:

[0152] [Eq-30] S:xl_>asexpxeQ° with as, Xs, 0S = X [0.42L» 0.39Ly, 0.02LJ and [8 m, 8 m, 8 m]. The source term is shown in 2D section in [Fig.5].

[0153] Two inversion experiments were performed with additive noise and using 1500 and 200 measurements distributed throughout the three-dimensional computational domain.

[0154] [Fig.5] The estimated source terms, as well as the resulting concentration fields are presented in Figure 5. The source positions xs are well estimated in both cases with errors less than 5 m on the x ” Y plane for the 1500 measurements case in Part A of the figure while xs is displaced by 25 m along x and 10 m along Y in Part C of the figure. Some inversion artifacts are visible in Part C where the source reaches 0.5 units instead of 0 near the building located to the southeast. However the sources are estimated well enough to reconstruct pollution field patterns in Parts B and D, close to the reference field returned by the fine model visible in [Fig.2].

[0155] The resulting field from an estimation of the source term at 1500 measurements (part B of [Fig.5]) however reconstructs the field better in order of magnitude than with 200 measurements (part D of [Fig.5]) because more information is provided.

[0156] Several parameters can impact the quality of a source term estimate. Among them, the number of nodes in X and Y and the weight of the regularization in the ratio b / 1 impacts the pattern of the estimated source and its order of magnitude.

[0157] The computation time is reasonably short. Once the operator matrix D is computed, the estimation of the source and its resulting field are obtained in about 4 minutes on this example with an Intel Xeon E-2276M 6-core processor and 32 Gb of RAM running a still poorly optimized python program.

[0158] [Fig.5] In Figure 5, the source estimate with 1500 measurements is shown in quadrant A, a resulting concentration field in quadrant B, a source estimate with 200 measurements in quadrant C and a concentration field reconstructed from 200 measurements in quadrant D. The crosses correspond to the planar coordinates of the measurements. The results are displayed at a height of z = 2 m.

[0159] [Fig.6] The PMSS software can qualitatively calculate the dispersion of pollutants. Using the data returned by the PMSS particle model where the plume is generated with a point source (i.e. generated in a single mesh) and with many PMSS model particles, a quality reference field Ÿo can be created (on the left in [Fig.6]) and used to create a set of simulated measurements: here measurements all located on a domain plane at a height of 2.14 m (on the right in [Fig.6]). It is then possible to go back to the source with an even more reduced LS-RBF-FD estimation model (up to 626) points in X and obtain a correct estimate of the source and its location in about thirty seconds.

[0160] Thus, on the left in [Fig.6], we see the PMSS reference field used to create simulated measurements, the displayed slice being 2.14 m high. On the right we see the position of measurements simulated by PMSS: an entire slice of PMSS domain at 2.14 meters high, 860 measurements having been provided.

[0161] [Fig.7] Source estimation results with Nx = 626 and 1865 points are visible in Thus, and show that the source position is very precisely estimated, including with a very reduced LS-RBF-FD inversion model at 626 points in X (on the left in Figure 7). The use of a finer LS-RBF-FD approximation model in X (Nx = 1865) however avoids the presence of source artifacts (Thus, on the right in [Fig.7]).

[0162] Thus, on the left, in Figure 7, we see an adjoint source estimate with a very reduced LS-RBF-FD inversion model: Nx = 626 points. On the right, we have shown an adjoint source estimate with a reduced LS-RBF-FD inversion model at Nx = 626 points.

[0163] The adjoint method based on RBF-FD models therefore works with noisy measurements, coming from a very fine LS-RBF-FD model and produces estimates of source terms with 2060 unknowns having 1500 or 200 measurements dispersed in the domain.

[0164] In the case where the simulated measurements are provided by PMSS, a model different from the LS-RBF-FD model, professionally developed and validated in a wind tunnel (Moussafir et al., 2004) and where the measurements are located at approximately 2 meters altitude, the localization of the source is very fast and very precise with the inverse method presented above. Other embodiments

[0165] The inversion procedure can integrate Kalman filtering which then makes it possible to update the estimate by adding measurements without having to redo a least squares calculation on an entire batch of measurements. A Kalman ensemble filter can be used and makes it possible to avoid having to explicitly calculate inverses. This is particularly used when we are interested in the transient regime. Kalman ensemble filters are suitable for large-dimensional problems.

[0166] According to yet another variant, an iterative adjoint procedure is carried out: this makes it possible to iterate towards the solution without having to store heavy matrices in memory, and this is therefore suitable for large domains.

[0167] According to yet another, we set up an efficient implicit time advancement scheme in calculation and using a physical solver to evolve the concentration states in time:

[0168] [Fig.8] [Fig.8] presents in this context an example of unsteady LS-RBF-FD evolution.

[0169] The transient physical solver serves as the basis for transient identification perspectives. It is formulated as follows: the unsteady problem creates a problem of

[0170]

[0171]

[0172]

[0173]

[0174]

[0175]

[0176]

[0177]

[0178]

[0179]

[0180]

[0181]

[0182]

[0183]

[0184]

[0185] minimization similar to the stationary problem where the calculated solution, $a as a reminder, satisfies the minimization problem, [Eq. 31] min II Due-f II ù In the unstationary problem, we add to the minimization constraint the time evolution of the field variable u. [Eq. 32] min || dtü + Dû - f || To advance in time in a discrete calculation scheme we consider a time following t + 1 to express a discrete derivative by an implicit scheme, the unstationary problem being translated into, [Eq. 33] m|n |lp -f-Dû(t+ l)-f || û(t+D ■ 7 With, Do an RBF-FD operator with derivative order 0 (interpolator) of Yn» in X when it acts on nodes yk e Yn» and 0 otherwise. So that, [Eq. 34] / \ 1 w°fv,X) if V eYn. where w0 _ is a set of DJ y = ] V k' / J k -2 \J k / i * kl I \ / t 0 otherwise RBF-FD weights to perform a local interpolation of ü for a stencil centered at n - 1 neighbors in X. So we solve [Eq. 35] mjn || F + j) _ ^)] + + A f || , with 00+1) 01 ' [Eq. 36] [Dk =Lk, Dr " And so we solve [Eq. 37] mn || +DWt+ 1)-(-Wtj + f) || û(t+D 11 df / VM; u This expression is solvable with LSQR or LBFGS-B, conjugate gradients or other solvers. Inverse modeling based on LS-RBF-FD is applicable to any search for unknowns, source terms or field variables in partial differential equations. Thus, it can be implemented for identification (localization) of hot spots in batteries, source estimation in a NRBC-E context (e.g. biological / chemical contamination in a water network, a triggered leak, etc.) For the pollution problem, the method is useful for identifying pollution sources for risk-related services (firefighters, army, etc.). The invention can also be used to locate acoustic sources in complex environments. List of cited references

[0186] Barnett, G.A., 2015. A Robust RBF-FD Formulation based on Polyharmonic Splines and Polynomials. University of Colorado Boulder.

[0187] Defforge, C., 2019. Data assimilation for micrometeorological applications with the fluid dynamics model Code_Satume. Paris-Est, Paris.

[0188] Flyer, N., Lehto, E., Biaise, S., Wright, G.B., St-Cyr, A., 2012. A guide to RBF-generated finite différences for nonlinear transport: Shallow water simulations on a sphere. J. Comput. Phys. 231, 4078-4095. https: / / doi.Org / 10.1016 / j.jcp.2012.01.028

[0189] Fomberg, B., Flyer, N., 2015. A Primer on Radial Basis Functions with Applications to the Geosciences, Society for Industrial and Applied Mathematics. ed.

[0190] Khodayi-mehr, R., 2019. Model-Based Learning and Control of Advection-Diffusion Transport Using Mobile Robots 285.

[0191] Kumar, P., Feiz, A.-A., Singh, S.K., Ngae, P., Turbelin, G., 2015. Reconstruction of an atmospheric tracer source in an urban-like environment. J. Geophys. Res. Atmosphères 120, 12589-12604. https: / / doi.org / 10.1002 / 2015JD024110

[0192] Liu, J., Li, X., Hu, X., 2019. A RBF-based differential quadrature method for solving two-dimensional variable-order time fractional advection-diffusion équation. J. Comput. Phys. 384, 222-238. https: / / doi.Org / 10.1016 / j.jcp.2018.12.043

[0193] Mons, V., Margheri, L., Chassaing, J.-C., Sagaut, P., 2017. Data assimilation-based reconstruction of urban pollutant release characteristics. J. Wind Eng. Ind. Aerodyn. 169, 232-250. https: / / doi.Org / 10.1016 / j.jweia.2017.07.007

[0194] Moussafir, J., Oldrini, O., Tinarelli, G., Sontowski, J., 2004. A new operational approach to deal with dispersion around obstacles : the MSS (Micro Swiftt Spray) software suite. 5.

[0195] Septier, F., Armand, P., Duchenne, C., 2020. A bayesian inference procedure based on inverse dispersion modelling for source terni estimation in built-up environments. Atmos. Environ. 242, 117733. https: / / doi.Org / 10.1016 / j.atmosenv.2020.117733

[0196] Tominec, I., Larsson, E., Heryudono, A., 2021. A least squares radial basis function finite différence method with improved stability properties (No. arXiv:2003.03132). arXiv.< / l>

Claims

Claims

1. Method for determining coordinates in the space of a source at the origin of a dispersion phenomenon in a domain, the dispersion phenomenon being approximated by a system of partial differential equations (20), the method comprising steps, implemented by computer, - of obtaining measurements (60) each quantifying at a point in the space of the domain and at a given instant a state variable modeled by the system of partial differential equations at said point and at said instant and - of iterative minimization (70) of the differences between the measurements and the calculations of the system of equations to locate a position (80) of said source, the method being characterized in that the system of equations is expressed (50) by finite differences supported by radial basis functions, the system of equations being solved in least squares, using a least squares solver,the system of equations having been the subject of an adjoint quadratic resolution by choosing radial basis functions among polyharmonic spline functions, multi-quadric inverse functions, or Gaussian functions.,

2. Method for determining coordinates according to claim 1, characterized in that the steps of obtaining the localized measurements and minimizing the differences between the measurements and the calculations are carried out in a described spatial domain (40) with Neumann and Dirichlet boundary conditions.

3. Method for determining coordinates according to claim 1 or claim 2, characterized in that the domain is characterized, in addition to its boundaries, also by field data (30) within the domain.

4. Determination method according to one of claims 1 to 3, characterized in that the dispersion phenomenon is a dispersion of pollutants and that said transport is modeled by an advection-diffusion equation.

5. Determination method according to one of claims 1 to 3, characterized in that the dispersion phenomenon is a diffusion of heat from a hot spot in an energy system.

6. Determination method according to one of claims 1 to 3, characterized in that the dispersion phenomenon is a sound diffusion in a complex environment.

7. Determination method according to one of claims 1 to 6, characterized in that the least squares solver is L-BFGS-B or LSQR.

8. Determination method according to one of claims 1 to 7, characterized in that the radial basis functions are the polyharmonic spline functions PSH5.

9. Determination method according to one of claims 1 to 8, characterized in that the iterative minimization (70) of the differences between the measurements and the calculations of the system of equations is carried out to also determine an intensity (80) of said source.

Images