Multifractal point processes and the spatial distribution of wildfires in French Mediterranean regions
MMultifractal point processes and the spatial distributionof wildfires in French Mediteranean regions
R. Ba¨ıle, J.F. Muzy ∗ , X. Silvani Laboratoire des Systemes Physiques de l’Environnement
CNRS UMR 6134 - University of Corsica
Campus Grimaldi, 20250 CORTE (France)
Abstract
We introduce a simple and wide class of multifractal point processes as Coxprocesses with a stochastic intensity corresponding to a multifractal measure.We then propose a maximum likelihood approach by means of a standardExpectation-Maximization procedure in order to estimate the distribution ofintensities at all scales and thus the scaling laws of its moments. The wild-fire distribution gathered in the Promethee French mediterranean wildfiredatabase is investigated within this framework. This allows us to computethe statistical moments associated with the spatial distribution of annuallikelihood of fire event occurence. We show that for each order q , these mo-ments display a well defined scaling behavior with a non-linear spectrum ofscaling exponents ζ q . From our study, it thus appears that the spatial distri-bution of the widlfire ignition annual risk can be described by a non-trivial,multifractal singularity spectrum and that this risk cannot be reduced toproviding a number of events per km . Our analysis is confirmed by a directspatial correlation estimation of the intensity logarithms whose the peculiarslowly decreasing shape corresponds to the hallmark of multifractal cascades.The multifractal features of the three regions that are studied appear to besimilar, constant over time and independent of the burnt areas associatedwith each fire event. Keywords:
Spatial point patterns Cox processes Multifractal measures ∗ Corresponding author
Email addresses: [email protected] (R. Ba¨ıle), [email protected] (J.F.Muzy), [email protected] (X. Silvani)
Preprint submitted to Elsevier June 4, 2020 a r X i v : . [ phy s i c s . s o c - ph ] J un atural fire Prom´eth´ee database Highlights • A multifractal spatial point process is defined as a spatial Cox processwith a multifractal random intensity measure. • The scaling and multifractal properties of this measure can be em-pirically estimated from samples by a maximum likelihood approachrelying on the EM method. • The annual wildfire ignition spatial distribution over three south FrenchMediterranean regions can be represented within this framework anddisplay strong multifractal scaling properties that appear to be similaraccross these regions. 2 . Introduction
All over the world, each year, wildfires are responsible of high hazard andrelated damages with strong impact on economic activity, biodiversity, de-crease in forest, soil degradation and greenhouse effects. Measuring and fore-casting wildland fire risk is therefore of prime importance for safety manage-ment services. This risk has mainly two components, namely the occurrenceprobability and the hazard that accounts for the severity of the consideredevent through e.g., the widlfire burned area and its expected damages. Inthis paper, we focus on the risk related to the ignition probability which hasbeen at the heart of a wide number of studies over the past decade (see e.g.Plucinski (2012); Prestemon et al. (2013); Mhawej et al. (2015) and referencestherein). Many of these studies tried to describe the temporal and spatialdependence in order to build tools and methods for evaluating or forecastingthe fire hazard and producing reliable operational information for help-to-decision making in prevention strategies (see, e.g., Genton et al. (2006); Xuand Schoenberg (2011); Ager et al. (2014); Rodrigues et al. (2014); Zhanget al. (2016)).Our purpose in this paper is to describe the spatial distribution of igni-tion risk at a coarser level and mainly to study its main statistical featuresthrough its scaling properties. As many natural hazards, forest fires involvenon-linear physical processes over a wide range of scales, yielding the sci-entific community to use tools and concepts from the physics of complexsystems to study them. Within this framework, universality, scaling andself-similarity have proven to be fruitful concepts to account for many quan-titative aspects of wildfire properties and to design pertinent phenomenolog-ical models. For instance, one of the main observation extensively studiedduring the last decade is the power-law behavior of the distribution of burnedareas (Malamud et al. (1998)) that has been considered within the theory ofself-organized criticality (SOC) (Bak et al. (1987)) or according to the mech-anism of Highly Optmimized Toterance (HOT) (Carlson and Doyle (2000)).According to the first scenario, natural fires are examples of self-organizedcritical systems (Ricotta et al. (1999, 2001); Turcotte and Malamud (2004))(i.e. a dynamical system “naturally” behaving near a critical state) while theHOT theory assumes that the system is somehow “optimized” by a naturalselection process in some state that makes it vulnerable to unusual condi-tions. More recently, several studies have revealed that scaling laws are alsoobserved on the temporal aspects of wildfire statistics, i.e. the statistics as-3ociated with inter event time distributions (Song et al. (2001); Telesca et al.(2005); Corral et al. (2008)).In the present study, we ambition to bring a contribution to this importantand promising field of research. In complement to the aforementioned studieson temporal distribution of wildfires, our aim is to address basic questionsrelated to their spatial distribution. For instance what is the dimension ofthe domain where fires are observed ? One can imagine that the ignitionlocations are distributed along few main roads or they can be uniformlyspread over a territory (e.g. a forestal zone): in the first case the fire hazardshould be measured in number of occurrences per kilometer while in the thirdcase it should be in number of events per kilometer squared. These simplesituations illustrate the fact that the estimation of the spatial properties offire distribution cannot be performed without the knowledge of the dimensionof the set supporting these events. On a more general ground, this shows thatstuduying the scaling properties of the spatial distribution of the ignition riskmeasure is an important issue. Unlike in the field of earthquakes (Molchanand Kronrod (2005)), this question remains, to a large extend, overlooked inthe literature devoted to wildland risks and has only been addressed by fewstudies. For example in Telesca et al. (2007); Tuia et al. (2008), the authorsstudy the fractal nature of fire distribution patterns in central Italy by themean of the correlation integral, box counting and sandbox methods andfound that ignitions are spread over a set a fractal dimension D c (cid:39) . N of fire ignitions in some area of size ε .The paper is organized as follows: in section 2, we define in which sensewe consider that a point process is multifractal and introduce the class ofmultifractal Cox processes. We then describe a maximum likelihood methodto estimate their multifractal scaling exponents. Our purpose is illustratedusing a simple example of a log-normal cascade spread over a Cantor set in1D. The application of this approach to wildfire data is provided in Sec. 3where we consider the spatial variations of the ignition annual rate in the”Prom´eth´ee” database for 3 regions in the French Mediterranean area. Aftera brief description of the database, we study the event clustering propertiesusing the Ripley L ( r ) function. We then study the spatial empirical distri-butions of annual number of ignitions and its multifractal scaling propertiesthat are compared for the 3 available regions. In the concluding section, wecomment our findings in terms of fire hazard and provide some prospects forfuture research.
2. Multifractal Cox processes and their estimation of their scalingproperties
Since a realization of a space-time point process over a bounded space-time interval, consists in a finite collection of points, it corresponds to set adimension D = 0. In that respect, when one refers to the fractal or multi-fractal nature of a point pattern there is a need to precisely define what oneexactly refers to. In Vere-Jones (1999), the author shows that the (multi-)fractality of a space-time point process can be either considered from thepoint of view of the scaling properties of the associated spatial intensity pro-cess or from its clustering properties. In the first case, one assumes thatthe process is ergodic in time and observed over a sufficiently long period sothat one can estimate the spatial fluctuations of the intensity measure (alsoreferred to as the expectation measure). In the second situation, the pro-cess is assumed to be homogeneous in space but events occur in a a strongly5orrelated way, meaning that the so-called Palm distribution is slowly de-creasing, i.e., behaves as a power-law. This latter approach is the one thatwas notaby taken in Ogata and Katsura (1991). On a general ground, onecan consider that the two previous features occur simultaneously and thatthe observed scaling properties is the intricate result of both intensity spatialinhomogeneities and event occurrence correlations.In this paper, we will mainly consider the first scenario proposed in Vere-Jones (1999) that basically consists in neglecting possible correlations be-tween event occurrence and focusing exclusively on the peculiar spatial fluc-tuations of the expectation measure that encodes all the non-trivial scal-ing features. More precisely, we will consider a point process dN ( (cid:126)x ) in R d ( d = 1 , , . . . , being the dimension of the embedding euclidian space) to bean homogeneous spatial Cox process that is a heterogeneous Poisson processwhich intensity distribution is itself an homogeneuous stochastic process (seeIllian et al. (2008) for the precise definition of a Cox process). This intensityis provided by a random measure Λ, i.e., we suppose that, for any spatialdomain B ⊂ R d : P rob { N ( B ) = n } = e − Λ( B ) Λ( B ) n n ! (1)with n ∈ N and Λ( B ) = (cid:90) B d Λ( (cid:126)x ) . (2)Let B ε ( (cid:126)x ) be a ball (or a square) of size ε centered at position (cid:126)x . If d Λ( (cid:126)x ) = ρ ( (cid:126)x ) d d x (where d d x stands for the Lebesgue measure in R d ), i.e., d Λ is de-scribed by a density function ρ ( (cid:126)x ), one has: ρ ( (cid:126)x ) = lim ε → Λ( B ε ( (cid:126)x )) ε d = lim ε → P rob { N ( B ε ( (cid:126)x )) = 1 } ε d (3)However, if the measure d Λ( (cid:126)x ) is singular as respect to the Lebesgue measure,the behavior Λ( B ε ( (cid:126)x )) ∼ ε d does not hold and a density ρ ( (cid:126)x ) cannot bedefined. In that case, one has generically:Λ( B ε ( (cid:126)x )) ∼ ε → + ε α ( (cid:126)x ) (4)where α ( (cid:126)x ) ∈ R + ∗ is the local singularity exponent of Λ at position (cid:126)x . This isprecisely the situation we want to account in this paper which occurs whenΛ is a multifractal measure, i.e., the local regularity is not the euclidiandimension ( α = d ) and strongly varies pointwise. The multifractal properties6f Λ can be described through the so-called multifractal formalism introducedin the context of turbulence Frisch and Parisi (1985) and widely used inmany areas ranging from chaotic dynamical systems to econo-physics. Thisformalism allows one to link the scaling properties of the measure with thestatistical distribution of its singularity exponents α ( (cid:126)x ). More precisely, onedefines, for q ≥
0, the partition function : Z ( q, ε ) = (cid:88) i,B ε ( (cid:126)x i ) ∈P ε Λ( B ε ( (cid:126)x i )) q (5)where P ε is a partition of S the observed support of Λ by boxes { B ε ( (cid:126)x i ) } i ofsize ε . From the small scale behavior of Z ( q, ε ), one defines the spectrum τ q of scaling exponents: Z ( q, ε ) ∼ ε → + ε τ q . (6)It is noteworthy that, from a practical point of view, ε → ε (cid:28) L where L is some well defined large scale in the problem called the integralscale . Eq. (5) is thus generally replaced by the following, more stringent,exact scaling relation that amounts to assuming a self-similarity property ofthe measure Λ: Z ( q, ε ) (cid:39) Z q (cid:16) εL (cid:17) τ q for 0 < ε ≤ L. (7)We can remark, that from definition (5), that τ = 0 (by the additivity of themeasure) and τ = − D c where D c is the fractal dimension (also called the“capacity”) of the set S that supports the measure Λ. In order to characterizethe distribution of local singularity exponents α , one introduces the so-calledsingularity spectrum f ( α ) defined as the fractal (Haussdorf) dimension ofthe iso-singulariy sets: f ( α ) = Dim H { (cid:126)x , α ( (cid:126)x ) = α } (8)Roughly speaking, this equation means that at scale ε , the number of boxeswhere Λ( B ε ( (cid:126)x i )) ∼ ε α is: N ( ε, α ) ∼ ε − f ( α ) . (9) In general Z ( q, ε ) can also be defined for q < dx ) that arevery sensitive to measurement errors and the size of the statistical sample. For q < f ( α ) and τ ( q ), as defined inresp. Eqs. (8) and (6), are Legendre transform each other: f ( α ) = min q ( qα − τ q ) τ q = min α ( qα − f ( α ))It results that q can be interpreted as a value of the derivative of f ( α )and conversely α is a value of the slope of τ q . Notice that in the particularcase when τ q is linear with a slope α , one recovers that fact that f ( α ) = α for the unique value α = τ (cid:48) q = α .As mentionned previously, they are many alternative formulations of theabove “box counting” method as for instance the correlation integral ap-proach or the sandbox method (Grassberger and Procaccia (1983); Feder(1998)). Instead of the partition functions (5), one can consider the follow-ing variant method that involves the the empirical moments of Λ( B ε ( (cid:126)x )): M ( q, ε ) = (cid:104) Λ( B ε ( (cid:126)x )) q (cid:105) = (cid:90) G ( z, ε ) z q dz (10)where we have denoted by (cid:104) . (cid:105) the empirial mean over the spatial support S of the measure Λ and G ( z, ε ) the spatial distribution of Λ( B ε ( (cid:126)x )) at (cid:126)x ∈ S .The scaling behavior of M ( q, ε ) defines the spectrum of exponents ζ q : M ( q, ε ) (cid:39) K q (cid:16) εL (cid:17) ζ q (11)which can be related to τ q as: ζ q = τ q + D c . (12)Indeed, since the total number of boxes at scale ε needed to cover S of isnothing but Z (0 , ε ), we have (cid:104) Λ( B ε ( (cid:126)x )) q (cid:105) (cid:39) Z ( q, ε ) Z (0 , ε )which leads (12), thanks to the scaling relationship (7) and to the equality τ = − D c . This entails notably: f ( α ) = D c + min q ( qα − ζ q ) (13)that relates the singularity spectrum to the spectrum ζ q of the scaling expo-nents of the empirical moments of Λ( B ε ).8 .2. A maximum-likelihood approach to estimate the multifractal propertiesof d Λ( (cid:126)x )In order to obtain the multifractal spectrum of a Cox point process, onethus need to estimate the empirical distribution of its expectation measurefluctuations at all scales, Λ( B ε ( (cid:126)x )). A standard way to do this is to use aparametric form of this distribution and, for each scale ε , estimate the param-eters by a maximum likelihood. However, if one observes only few realizationsof the inhomogeneous Poisson process associated with a given intensity mea-sure d Λ( (cid:126)x ), the latter is not directly observable since one gets, for each (cid:126)x andeach box size ε , only few samples of the random variables N ( B ε ( (cid:126)x )) drawnaccording to the Poisson law (1). If the parametrization of the intensitydistribution is a sum (or a “mixture”) of several simple distributions, theobservable distribution will be a mixture of compound Poisson random vari-ables. A standard way to estimate the parameters of this mixture in order tomaximize the (log-) likelihood is to use an Expectation-Maximization (EM)procedure (see e.g. Hastie et al. (2001)). To be more specific, let g ( z, Θ ) bea family of probability density functions over R + ∗ of parameters Θ and letus represent the distribution G ( z, ε ) of Λ( B ε ( (cid:126)x )) as the following mixture: G ( z, ε ) = J (cid:88) k =1 w k g ( z, Θ k ) . (14)where the dependence in the scale ε relies in all parameters J , w k ’s and Θ k ’sBy definition of N ( d(cid:126)x ) as an inhomogeneous Poisson process of intensity d Λ( (cid:126)x ), P ( n, ε ), the empirical distribution of observed events N ( B ε ( (cid:126)x )), canbe written as: P ( n, ε ) = P rob { N ( B ε ( (cid:126)x ) = n } = 1 n ! (cid:90) ∞ e − z z n G ( z, ε ) dz and therefore, by defining: h ( n, Θ ) = 1 n ! (cid:90) ∞ e − z z n g ( z, Θ ) dz . we have the finite mixture representation of P ( n, ε ): P ( n, ε ) = J (cid:88) k =1 w k h ( n, Θ k ) (15)9hich parameters { w k , Θ k } k =1 ,...,J can estimated using an EM method. Thehyperparameter J can be chosen using a BIC or AIC selection criterion.Thanks to Eq. (10), the partition function M ( q, ε ) is then estimated as: M ( q, ε ) = J (cid:88) k =1 w k (cid:90) z q g ( z, Θ k ) dz . (16)In practice, we will consider the distribution g ( z, Θ ) to be either a Gammadistribution or a log-Normal distribution that have both 2 parameters. Thecorresponding compound Poisson distributions h ( n, Θ ) are respectively Ne-gative Binomial distributions and Poisson-Log-Normal distributions. Noticethat in the latter case, no closed-form formula is available for h ( n, Θ ) and anumerical evaluation of the integral has to be performed. Let us mention thatour approach is precisely the one formerly developed in Hwa (1995); Jie andShaoshun (1997) in the context of particle physics where the authors proposedto filter out the “statistical fluctuations” (i.e. the Poisson random part) byestimating the “factorial moments of continuous order” that correspond tothe standard moments M ( q, ε ) associated with the distribution G ( z, ε ) of theintensity measure. For that purpose, they described the observed number ofevents distribution as a mixture of negative binomial distributions. R . In order to illustrate our approach, let us consider the simple example ofa Cox process in R ( d = 1) which intensity is a log-infinitely divisible ran-dom multiplicative cascade. Such cascade models represent the paradigm ofmultifractal processes with well defined exact scaling properties (Muzy et al.(2000); Bacry and Muzy (2003)). More precisely, we study a Cox processwith an intensity d Λ( x ) that is provided by the lacunary random cascademodel defined in Muzy and Ba¨ıle (2016). This model consists in building alog-infinitely divisible random cascade that is supported by a Cantor set ofarbitrary dimension D ≤
1. In our example, we consider a measure Λ( dx ) dis-tributed on a set of dimension D = 0 . λ = 0 .
05 and integral (i.e. maximum correla-tion) scale L = 512 (see Muzy and Ba¨ıle (2016) for the precise meaning ofthese parameters). 10
500 1000 1500 2000 2500 3000 3500 40000510 Λ ( d x ) ( a ) x N ( d x ) ( b ) Figure 1: (a) Sample of a log-Normal multifractal measure spread over a Cantor set ofdimension D = 0 .
8. The integral scale is L = 512 and the intermittency coefficient λ = 0 .
05. (b) Realization of an inhomogeneous Poisson process associated with themultifractal intensity displayed in (a).
Such a measure can be shown to be multifractal with a moment scalingexponents and singularity spectrum that are quadratic functions: ζ q = α q − λ q (17) f ( α ) = D − ( α − α ) λ with α = D + λ dx ) (where we choose numerically dx = 1) over 8 integralscales in represented in Fig. 1(a). A realization of the associated countingprocess N ( dx ) (each interval dx contains a random number drawn with aPoisson law of intensity Λ( dx )) is displayed in Fig. 1(b). The latter processcan be viewed as a “noisy” version of Λ( dx ) where Poisson statistical fluctua-tions are superimposed to the spatial log-normal cascade noise. If one wantsto characterize the genuine “risk”, i.e., the statistical properties of Λ( dx ),one has to get rid of these Poisson fluctuations. This is the purpose of the11 J . . . . . . . B i c S c o r e ( a r b i tr a r y un i t s ) Figure 2: Baysian information criterion (“BIC Score”) for the selection of the mixturenumber J for scales ε = 5 (blue), ε = 13 (green) and ε = 40 (orange). The units havebeen chosen and the range has been shifted for display purpose (in particular the minimumlikelihood has been arbitrary set to zero). One can see that J = 3 provides in each casethe best score. previously described maximum likelihood method relying on representations(14), (15) that aims at recovering G ( z, ε ), the distribution of Λ( B ε ( x )) from P ( n, ε ), the distribution of N ( B ε ( x )). We choose to represent P ( n, ε ) as amixture of J Negative Binomial distributions which amounts to representing G ( z, ε ) as a mixture of J Gamma distributions. We performed our numericalestimation using a sample N ( dx ) over a total length of 256 integral scales.Our analysis was performed over a range of scales ε such that 1 ≤ ε ≤ J can be determined using a BIC or AIC selec-tion criterion. As illustrated in Fig. 2, we found that, for all scales, J = 3achieves (or almost achieves) the maximization of the penalized likelihood.12
20 40 60 80 100 120 140 n . . . . P ( n , ε ) ( a ) z . . . . . G ( z , ε ) ( b ) Figure 3: (a) Negative Binomial mixture (solid line) and empirical ( • ) distribution of N ( B ε ( x )) for two scales ε = 13 and ε = 40. (b) Associated Gamma mixture distribution G ( z, ε ) as compared to the observed one computed directly from the log-Normal sampleof Λ( dx ). The performance of the parametric maximum likelihood method to fitthe empirical distributions of N ( B ε ( x )) at each scale is illustrated in Fig. 3for scales ε = 13 and ε = 40. One can see that in both cases, the observedempirical distributions are fitted fairly well by a Negative Binomial mixturewith J = 3. It results that G ( z, ε ), the original distributions of Λ( B ε ( x )),are also well fitted by the associated mixtures of Gamma distributions (Fig.3(b)). From these mixtures, M ( q, ε ) is then computed at each scale ε usingEq. (10) and from a linear fit a ln( M ( q, ε )) as a function of ln ε one getsan estimate of the spectrum ζ q . The estimation obtained from the sample oflength 256 L in our example is reported in Fig. 4(a). The shaded grey regionindicates an order of magnitude of the estimation error when one changes therange of scales used to perform the fit from smallest scales to largest ones.13 q ζ q ( a ) .
60 0 .
65 0 .
70 0 .
75 0 .
80 0 . α − . . . . . . f ( α ) ( b ) Figure 4: (a) Estimated ( • ) ζ q spectrum of moment scaling exponents as compared tothe analytical expression (17) (solid line). (b) Estimated f ( α ) singularity spectrum ( • )as obtained by the numerical Legendre transform of the ζ q estimation. The solid linerepresents the expected parabolic spectrum as given by Eq. (18). In both figures theshaded region stands for the observed variations in the estimated spectra when one changesthe range of scales to perform the linear fit in log-log representation of the moments. We can see that the numerical estimation procedure detailed below allowsone to recover very precisely the expected parabolic ζ q function. It results,using a numerical Legendre transform that one can estimate also quite wellthe shape of the singularity spectrum for this model (Fig. 4(b)).
3. The multifractal approach applied to the Prom´eth´ee database
The Prom´eth´ee database was created in 1973 in order to gather sev-eral informations relative to the wildfire in French Mediterranean regions.It reports many different features (year, region administrative number, firenumber, geographical coordinates, date, burnt surface, characteristics of thefirst fire fight action, nature of the damages, vegetal species, information rel-ative to the fire cause, ...). of the wildfire occurred in 15 French departments.The information flux sources come from various national services acting ineach region (fire fighters, forest managers, police, civil safety, army and airforce). Some efforts to homogenize these multiple data were performed inthe 80’s, especially through a new system of coordinates specically designedfor fire management: the DFCI coordinates (DFCI is the French acronym for14D´efense de la Forˆet Contre les Incendies”, i.e. it refers to all processes con-cerning forest defense facing to wildfires). This coordinate system consists ina set of nested grid layers with increasing resolution going for 100 × km to 2 × km at the finest resolution. The fire ignition locations are thereforeavailable with the rather poor spatial resolution of 4 km .In order to handle consistent data and to avoid biased results, we chooseto process the database restricted from the 1 st of January, 1992 to the 31 st ofDecember, 2018. We consider separately three main regions namely, “Cor-sica”, “Provence-Alpes-Cˆote d’Azur” (PACA) and “Languedoc-Roussilon”(LR). The main statistics of the database are summarized in table 3.1 wherewe see that the number of reported ignitions is larger in Corsica than in thetwo other regions whilst it is almost four times less wide.Region Surface ( km ) N tot N S> N S> N S> Corsica 8680 21073 8872 4666 1449PACA 31400 18213 10734 3472 1151LR 27376 12552 9155 5143 1396
Table 1: Main statistical features of the Prom´eth´ee database: we reported the overallsurface of each regions, the number of forest fire events observed from 1992 to 2018, andthe number of the events where the burnt area is greater than respectively 1000,10000 and50000 m . The fire ignition process over a given region can be considered as a spatio-temporal Point process and therefore the intensity (or expectation) measure d Λ( (cid:126)x, t ), i.e., the mean number of events between t and t + dt in a square of size d x located at (cid:126)x , depends on both (cid:126)x and t . Notably, as already discussed informer studies (see e.g. Zhang et al. (2014); Bajocco et al. (2017)), Λ( d x, dt )displays strong annual seasonality, the ignition rate being much larger in thedry summer season than in the winter season. Since our goal is to studythe statistical properties of the spatial fluctuations of the fire occurrencelikelihood, in order to avoid these seasonal effects and to consider a purelyspatial Point process, we focus in the paper on the annual ignition rate.Hereafter, d Λ( x, y ) will stand for the intensity associated with the annualnumber of events at some given location (cid:126)x = ( x, y ) over an infinitesimalsquare of size dxdy . Let us denote by S the support of d Λ, namely the setwhere d Λ( x, y ) is non vanishing. It can be empirically defined as the set ofcenters (cid:126)x k of the 4 km DFCI cells (the smallest available resolution) where15 .0 100.0 200.0 300.0 400.0 500.0 600.00.050.0100.0150.0200.0250.0300.0350.0400.0 0.00.030.130.320.611.0
Figure 5: Map (cid:98) Λ k of annual fire occurrence rate for the Southern French Mediteraneanregions as estimated from Prom´eth´ee database over the period 1992-2018. there has been at least one event between during the whole sample period1992-2018. Let S be the cardinal of S and let us denote by { Λ k } k =1 ,...,S the intensity associated with DFCI cell in S , i.e., Λ k = Λ( B ( (cid:126)x k )). Alongthe same line we will denote by N k the (random) number of events in cell k during a year: N k = N ( B ( (cid:126)x k )). In practice, a surrogate for Λ k can beobtained as the mean number of events per year observed during the whole27 years period: (cid:98) Λ k = 127 (cid:88) y =1992 N k ( y ) ≡ E ( N k ) (19)where we have denoted by N k ( y ) the realization of N k at year y and definedthe expectation E ( . ) as the average across all 27 years of the sample. InFigure 5, we reported such an annual fire occurrence mean density observedin each DFCI square in the French southern regions. Density levels areindicated in number of ignition per km and per year. We see that thedistribution appears very inhomogeneous in space, with obvious clustering16
10 20 30 40 50 60 r ( km ) L ( r ) Figure 6: Estimated Ripley inhomegeneous L ( r ) = (cid:112) K ( r ) function for Corsica,PACAand LR regions. The solid line stands for L ( r ) = r as expected for uncorrelated eventnumbers N k ’s. properties, high ignition levels being observed close to main road axes anddensely populated areas. In order to describe the spatial fluctuations of the annual number of fireoccurrence dN ( x, y ), we will use the approach described in Sec. 2, wherewe suppose that dN ( x, y ) is an Cox process, i.e., conditionally to randommultifractal spatial intensity d Λ( x, y ), dN ( x, y ) is an inhomogeneous Poissonprocess. This assumption notably implies that, for any given the intensityfunction, the observed number of ignitions during a given year over distinctareas are uncorrelated. This means that the observed spatial clustering ofignition locations is exclusively due to the spatial fluctuations of the intensityfield. In order to check for such a feature we follow the method proposedin Hering et al. (2009) where the authors define a inhomogeneous version of17he Ripley K -function that allows one to filter out the spatial dependenceof the intensity and test to remaining existing correlations. Accordingly, onedefines: K ( r ) = S − (cid:88) (cid:126)x k ∈S (cid:88) (cid:126)x j ∈ B r ( (cid:126)x k ) ∩S\ (cid:126)x k E [ N k N j ] w jk (cid:98) Λ k (cid:98) Λ j (20)where S , S , N k , (cid:98) Λ k and the expectation E ( . ) are defined in the previoussection and w kj is a Ripley edge correction factor designed to correct biasescaused by the edges of the domain (see Hering et al. (2009)) . Notice that inabsence of correlation, E [ N k N j ] = (cid:98) Λ k (cid:98) Λ j and therefore L ( r ) = (cid:112) K ( r ) = r .We have computed L ( r ) according to expression (20) for the 3 regions. As itcan be seen in Fig. 6, the plots of L ( r ) closely follow the straight line L ( r ) = r in the 3 cases which suggests that the random variables N k and N j are un-correlated. This result confirms the finding of Hering et al. (2009) and showsthat observations are compatible with a Cox process. Notice that it mightbe quite surprising that one does not observe any spatial anti-correlation be-tween fire occurrence events since one expects that after a wildfire anotherone cannot occur nearby within an already burnt area. However one has toremind that the minimal considered surface is our study is 4 km which isquite huge as respect to the typical wildfire burnt areas. Very large fires thatwill contribute to an effective anti-correlation are very few and statisticallyinsignificant. In this section we report the empirical results we obtained using the dataof the 3 regions when estimating the law of Λ( B ε ( (cid:126)x )) and then the multifractalscaling laws ζ q , f ( α ) following the method of Sec. 2.2. We chose to computethe probability distribution at points { (cid:126)x k } k =1 ,...,S that are the centers of theDFCI cells in the set S where there has been at least one event. Let usmention that, along the same line than in previous section, we have takeninto account the edge effects by applying a Ripley edge correction factor w ij Actually, we adapted the usual Ripley correction factor from circle geometry to squaregeometry. If d kj = max( | x j − x k | , | y j − y k | ) stands for the square distance, w kj representsthe fraction of DFCI squares m at distance d kj from (cid:126)x k that are in the studied region.
50 100 150 200 250 300 n . . . . . . . P ( n , ε ) n − − − − l n P ( n , ε ) n . . . . . P ( n , ε ) n − − − − − l n P ( n , ε ) ( a )( b ) ( c )( d )( a )( b ) ( c )( d )( a )( b ) ( c )( d ) Figure 7: Spatial probability distributions of N ( B ε ( (cid:126)x )) for ε = 6 , ,
22 km (and 38 km in(b,d)) in Corsica and LR regions. Panels (a) and (c) are in linear scales while (b) and (d)are in logarithmic scale in order to emphasize the tails of the distributions. The left plots(a,b) correspond to estimates obtained when accounting for all forest fire events Corsicaand the right plots (c,d) correspond to fire events of same type in LR region. Symbols ( • )represent the empirical data and the solid lines are the best fit obtained with a mixtureof J = 3 Negative Binomial distributions. when estimating, for each (cid:126)x k , N ( B ε ( (cid:126)x k )) from observations at smallest scale( ε = 2) N j : N ( B ε ( (cid:126)x k )) = (cid:88) (cid:126)x j ∈ B ε ( (cid:126)x k ) w − kj N j We estimate the empirical distributions of N ( B ε ( (cid:126)x k )) for ε in the range[2 , km of each region for all reported forest fire events and also when werestrict the sample to burnt area larger than 1000 , m . Foreach scale, we compute the best representation of the distribution in termsof a mixture Negative Binomial or Poisson-Log-normal distributions. Sinceboth approaches lead to similar results, we only report the results relative to19 ln ε l n M ( q , ε ) ln ε l n M ( q , ε ) ( a ) ( b ) Figure 8: Scaling properties of intensity moments M ( q, ε ). ln M ( q, ε ), as defined in Eq.(10), is plotted as a function of ln ε for q = 1 , . , , . . . , Negative Binomial distributions that are easier to handle from a numericalpoint of view. As illustrated in Sec. 2.3, the hyper-parameter J can bechosen, at each scale, by the mean of a BIC selection criterion. We foundthat in all cases J = 3 provides a near optimum choice for all scales.In Fig. 7, we have displayed the empirical distributions and their fit usinga Negative Binomial mixture, of annual number of fire occurrences in boxes ofdifferent sizes in the Corsica and Languedoc Roussillon regions (plots for thePACA region are similar). We can see that the shape of these distributionsstrongly depend on the considered scale and appears to be over-dispersed asrespect to a simple Poisson law. In both linear (Figs.7(a,c)) and logarithmicplots (Figs.7(b,d)), we can also see that the mixture model provides a verygood fit of the empirical data, around their maximum values as well as intheir tail behavior. The scaling behavior of the intensity measure moments20 q ζ q ( a ) . . . . . . α . . . . . . . . . f ( α ) ( b ) Figure 9: Multifractal spectra ζ q and f ( α ) for the three regions as obtained from theempirical scaling properties of M ( q, ε ). The grey shaded regions correspond to the observedfluctuations in spectra for Corsica data when one changes the fitting interval towards largeand small scales. Blue, orange and green symbols correpond respectively to Corsica, LRand PACA regions. can then be computed at all scales from these mixture estimations. FromEq. (10), one can thus compute the q order moment at all scales ε and itsscaling exponent ζ q . In Fig. 8, the estimated M ( q, ε ) are displayed in log-logrepresentation for q ∈ , . , , . . . , M ( q, ε ) is rather well modeled by power-lawover a range extending up to ε (cid:39) km . At larger scales, (cid:15) (cid:38) km , it seemsthat there is a regime change toward smaller exponent values. This featurecan be related to boundary effects that necessarily break the scaling laws atlarge scales but can also be explained by the existence of a large correlation(integral) scale L inherent to each multifractal process (see below).The values of ζ q and f ( α ) spectra we estimated for each of the 3 regionsof the database are reported in Fig. 9. One can see that all estimated spec-tra have a well pronounced strictly concave shape specific of a multifractal21 . . . . . . . ln( ε ) . . . . . ∂ l n M ( q , ε ) ∂ ε | q = Figure 10: Estimating the intermittency coefficient of the intensity of distribution. ∂ ln M ( q,ε ) ∂q (cid:12)(cid:12) q =0 is plotted as a function of the scale logarithm ln ε for Corsica (blue), LR(orange) and PACA (green) regions. The linear deacreasing behavior provides, accordingto Eq. (22) an estimation λ (cid:39) . measure. Moreover, these empirical functions are very close to each otherindicating a somehow universal character of the multifractal properties of d Λ( (cid:126)x ) in the French Mediterranean region. We notably find that for all 3distributions the forest fire ignitions occur on a set of dimension D c (cid:39) . α (cid:39)
2. This fractal dimension of the support of theignition intensity, close to D = 2, is greater than the vales D c (cid:39) . α ∈ [1 . , .
2] forall three regions. It is tempting to interpret these extreme values as being as-sociated respectively with line and simple area geometries while this remainsto be confirmed by a study focusing on local properties that would probablyrequest better resolved spatial data.In the log-normal case, as in the example consider in Sec. 3, the intensityof the multifractality is entirely characterized by the intermittency coefficient (Eqs. (17), (18)). On a general ground, this coefficient quantifying thenon-linear character of ζ q is usually defined as (see. e.g., Frisch (1995); Muzyand Ba¨ıle (2016)): λ = − ζ (cid:48)(cid:48) = − ∂ ζ q ∂q (cid:12)(cid:12)(cid:12)(cid:12) q =0 . (21)It results from Eq. (11) that ∂ ln M ( q, ε ) ∂q (cid:12)(cid:12)(cid:12)(cid:12) q =0 = − λ ln (cid:16) εL (cid:17) + V (22)with V = ∂ ln K q ∂q (cid:12)(cid:12) q =0 . This equation can be interpreted quite easily withina multiplicative cascade picture. Indeed, ∂ ln M ( q,ε ) ∂q (cid:12)(cid:12) q =0 is nothing but thevariance of ln Λ( B ε ( (cid:126)x )) as respect to its spatial fluctuations. When one goesfrom fine to large scales, this variance decreases as linear function of ln( ε ),meaning that each time one divides the resolution by, e.g., a factor 2, oneadds a random term ω of constant variance to ln Λ( B ε ( (cid:126)x )), i.e. the value ofΛ( B ε ( (cid:126)x )) is multiplied by a random factor W = e ω and thus it is a randommultiplicative cascade. The slope of this linear function, i.e., the quantity λ ln(2), thus corresponds to the variance of ω = ln W . Notice that withinthis picture, V is interpreted as the large scale variance of ln Λ( B ε ( (cid:126)x )) that,when ε ≥ L , no longer depends on ε . Eq. (22) provides a simple way todirectly estimating λ from empirical data. In Fig. 10 are reported thesecond order derivative of ln M ( q, ε ) around q = 0 (as approximated by afinite difference scheme with ∆ q = 0 .
2) for Corsica, LR and PACA regions.The three curves appear strikingly to be parallel with a slope λ (cid:39) .
4. Thismeans that the 3 spatial distributions of fire ignitions can be described bythe same cascade model. The precise value of the integral scale L is hard toestimate without the knowledge of V but one can evaluate the magnitudeorder of their ratio (provided V remains constant across regions).As emphasized in (Muzy et al. (2000); Bacry and Muzy (2003)), a simplemethod to estimate the integral scale and the intermittency coefficient is tostudy the spatial dependence covariance of ln Λ( B (cid:96) ( (cid:126)x )) that, according to themultifractal cascade picture, should behave when (cid:96) ≤ | (cid:126)x − (cid:126)x | ≤ L as: Cov { ln Λ( B (cid:96) ( (cid:126)x )) , ln Λ( B (cid:96) ( (cid:126)x )) } (cid:39) λ ln L | (cid:126)x − (cid:126)x | (23)23 . . . . . . . ln( ε ) . . . . . . . ˆ C l n N t o t ( ε ) Figure 11: Empirical covariance (Eq. (24)) of a the logarithmic intensity surrogate as afunction of the logarithm of the lag ε for Corsica (blue), LR (orange) and PACA (green)regions. A multiplicative bias correction calibrated using a synthetic cascade model asbeen applied to each curve. where B (cid:96) is a small box of size (cid:96) and the covariance has to be understood asunder the law of spatial fluctuations. This equation, which can be shown tobe a direct consequence of the existence of a mutliplicative cascade process(Muzy et al. (2000); Bacry and Muzy (2003)), means that when one plots thecovariance of the logarithm of the intensity as a function of the logarithm ofthe spatial distance, one gets a straight line of slope − λ and intercept ln L (see also Ba¨ıle and Muzy (2010)). Since the values of intensity field d Λ( (cid:126)x )are not observable (which precisely motivated the use of EM method forthe moment estimation) one cannot directly check the validity of Eq. (23).However, we can use the surrogate intensity (cid:98) Λ k introduced in Sec. 3.1 (Eq.(19)) by collecting all the ignitions events at a given spatial location over the24hole period of 27 years and estimate expression (23) as: (cid:98) C ln N tot ( ε ) = (cid:28) ln (cid:98) Λ( B ( (cid:126)x )) ln (cid:98) Λ( B ( (cid:126)x )) (cid:12)(cid:12)(cid:12)(cid:12) | (cid:126)x − (cid:126)x | = ε (cid:29) − (cid:68) ln (cid:98) Λ( B ( (cid:126)x )) (cid:69) (24)where again (cid:104) . (cid:105) has to be understood as a mean as respect to spatial posi-tions over the support of Λ. It is noteworthy that, when computing such anempirical covariance, one observes a bias as respect to the true intermittencycoefficient that depends on the amplitude of (cid:104) Λ k (cid:105) . Indeed, the greatest theintensity, the smallest the size of relative fluctuations of (cid:98) Λ k and therefore thesmallest this bias. The exact dependence of this bias as a function of (cid:104) Λ k (cid:105) can be hardly expressed analytically and has been calibrated using the toymodel described in Sec. 2.3. In Fig. 11, we have plotted the bias-correctedempirical covariance (24) as obtained from the three regions dat. It is strik-ing that in for all 3 regions, the logarithm of the local intensities appears tobe strongly spatially correlated over large distances with a correlation func-tion that decreases logarithmically, precisely as one expects for a randomcascade model (Eq. (23)). The values λ (cid:39) . L as the interceptof the empirical curves. We find respectively L Corsica (cid:39)
30 km, L LR (cid:39)
50 kmand L P ACA (cid:39)
90 km. Notice however that the uncertainty of these valuesis quite large due the number of possible biases (linked to the assumptionswe made, the quality of the data,....) and to the finite sample statisticalfluctuations. We can say that the order of magnitude of the integral scale isaround 50 km.To conclude our empirical study of multifractal properties of fire ignitionspatial disribution, we have checked that the spectra ζ q and f ( α ) are stableover the considered time period by performing estimations over a slidingwindows of 6 years from 1992 to 2012. Up to some fluctuations for thehighest singularities (i.e., the smallest values of α corresponding to large q values and thus to tail events) we observed that the results are consistent.We also checked that the observed scaling properties do not depend on theevent intensities, namely on the size of the considered wildfires. This isillustrated in Fig. 12 where are displayed the estimated multifractal spectrafor ignitions intensities in Corsica corresponding to all events, events with25 q ζ q ( a ) .
25 1 .
50 1 .
75 2 .
00 2 . α . . . . . . . f ( α ) ( b ) Figure 12: Estimation of (a) ζ q and (b) f ( α ) spectra for annual intensities of all fireevents in Corsica (blue symbols), events with a final burnt area greater than 1 ha (orangesymbols) and events with a burnt area greater than 5 ha (green symbols). final burnt area greater than 1 ha and events with a burnt area greater than5 ha. We see that all curves can be hardly distinguished and thus that thespatial distribution of ignition risk appears to be unrelated to the area of theburnt vegetation. This contrasts with the situation observed for earthquakedata where the scaling properties turn out to depend on the event magnitudes(see e.g., Molchan and Kronrod (2005)). These findings seem also to be incontradiction with the results reported by Telesca et al. (2007).-
4. Summary and prospects
In this paper we presented a new method to estimate the multifractalproperties of point patterns with clustering features when these latter resultfrom the spatial fluctuations of the expectation measure and not from pe-culiar correlations in the event occurrence likelihood. The paradigm of suchprocess is a spatial Cox processes with an intensity measure that is providedby a random cascade model. When only a few number of realizations areavailable so that the intensity measure remains unknown, we have shown26hat the moments of this measure at each scale can be still be estimatedthrough a maximum likelihood approach that consists on representing theobserved distribution number of events as a mixture of simpler distributions(like e.g. Negative Binomial distributions). The model calibration can thenbe performed using a classical Expectation Maximization procedure. Ourapproach has been validated on a simple toy model involving a log-normalmultifractal intensity lying on a random Cantor set.We have applied this framework to the annual wildfire ignition eventsof French Mediterranean regions gathered in the Prom´eth´ee database. Weshown that the ζ q exponent spectrum governing the power-law behavior ofthe order q moment of the intensity distribution behaves as a strictly concavefunction, the hallmark of multifractal processes. All the three studied regionsexhibit almost the same multifractal features: A dimension of the supportclose to D c (cid:39) . α (cid:39) .