Modelling and simulation of multifractal star-shaped particles

TL;DR

This paper introduces a flexible stochastic modeling framework for 3D star-shaped particles with variable fractal dimensions, using anisotropic Gaussian fields and a novel simulation algorithm, demonstrated through multifractal Earth topography examples.

Contribution

It presents a new approach to model and simulate multifractal star-shaped particles with variable Hausdorff dimensions using anisotropic Gaussian fields and a spectral turning bands-based algorithm.

Findings

Successful simulation of multifractal Earth topography.

Demonstration of the model's ability to capture shape variability.

Introduction of a new simulation algorithm for complex fractal particles.

Abstract

The problem of constructing flexible stochastic models to describe the variability in shape of solid particles is challenging. Natural objects often exhibit mono- or multi-fractal features, i.e. irregular shapes and self-similar patterns. This paper presents a general framework for modelling three-dimensional star-shaped particles with a locally variable Hausdorff (or fractal) dimension. In our approach, the radial function of the particle is represented by an anisotropic Gaussian random field on the sphere. We additionally derive a simulation algorithm being parenthetical to the spectral turning bands method proposed in Euclidean spaces. We illustrate the use of our proposal through numerical examples, including a multifractal simulated version of the Earth topography.

Equations

$$C({\mathbf{x}}_1,{\mathbf{x}}_2)=K\{d({\mathbf{x}}_1,{\mathbf{x}}_2)\},{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{S}}^2,$$

$$\mathrm{var}\left\{\sum _{i=1}^k{a}_{i}Z({\mathbf{x}}_i)\right\}=\sum _{i=1}^k\sum _{j=1}^k{a}_i{a}_{j}C({\mathbf{x}}_i,{\mathbf{x}}_j)\ge 0.$$

$$K(d)=\sum _{n=0}^{\infty }{b}_n{P}_n(\cos d),0\le d\le \pi ,$$

$$b_n=\frac{2n+1}{2}{\int }_0^{\pi }{P}_n(\cos \xi )\sin (\xi )K(\xi )\text{d}\xi ,n\in {\mathbb{N}}_0.$$

$$K(0)-K(d)\sim c_0{d}^{\alpha },$$

$$\underset{t\to \infty }{\lim }\frac{f(rt)}{f(t)}=1.$$

$$K(0)-K(d)\sim 2^{-\alpha }{d}^{\alpha }f(1/d),$$

$$\sum _{k=n}^{\infty }{b}_k\sim f(n)n^{-\alpha },$$

$$b_n=({\tau }^2+n^2{)}^{-\nu -1/2},n\in {\mathbb{N}}_0,$$

$$\sum _{k=n}^{\infty }({\tau }^2+k^2{)}^{-\nu -1/2}\sim n^{-2\nu }(2\nu {)}^{-1},$$

$$Y=\bigcup _{\mathbf{x}\in {\mathbb{S}}^2}\{\mathbf{o}+r\mathbf{x}:0\le r\le Z(\mathbf{x})\}.$$

$$H^{\eta }(Y)=\underset{ϵ\to 0}{\lim }\inf \left\{\sum |B_i{|}^{\eta }:\{B_i:i=1,2,\dots \}\text{ is an }ϵ\text{-cover of }Y\right\}$$

$$C({\mathbf{x}}_1,{\mathbf{x}}_2)=\sum _{n=0}^{\infty }{\left\{b_n({\mathbf{x}}_1)b_n({\mathbf{x}}_2)\right\}}^{1/2}{P}_n({\mathbf{x}}_1^{\top }{\mathbf{x}}_2),{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{S}}^2,$$

$$C({\mathbf{x}}_1,{\mathbf{x}}_2)\approx \sum _{n=0}^{\infty }{b}_n({\mathbf{x}}_0)P_n({\mathbf{x}}_1^{\top }{\mathbf{x}}_2),{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{S}}^2.$$

$$b_n(\mathbf{x})=({\tau }^2+n^2{)}^{-\nu (\mathbf{x})-1/2},n\in {\mathbb{N}}_0,\mathbf{x}\in {\mathbb{S}}^2,$$

$$Z(\mathbf{x})=\sum _{j=1}^m{\lambda }_j(\mathbf{x})Z_j(\mathbf{x}),\mathbf{x}\in {\mathbb{S}}^2,$$

$$C({\mathbf{x}}_1,{\mathbf{x}}_2)=\sum _{j=1}^m{\lambda }_j({\mathbf{x}}_1){\lambda }_j({\mathbf{x}}_2)K_j\{d({\mathbf{x}}_1,{\mathbf{x}}_2)\},{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{S}}^2.$$

$$C({\mathbf{x}}_1,{\mathbf{x}}_2)=\sum _{n=0}^{\infty }\left\{\sum _{j=1}^m{\lambda }_j({\mathbf{x}}_1){\lambda }_j({\mathbf{x}}_2)p_{n,j}\right\}{P}_n({\mathbf{x}}_1^{\top }{\mathbf{x}}_2),{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{S}}^2.$$

$$Z(\mathbf{x})={\left\{\frac{b_{\kappa }(\mathbf{x})(2\kappa +1)}{a_{\kappa }}\right\}}^{1/2}{P}_{\kappa }({\mathbf{w}}^{\top }\mathbf{x}),\mathbf{x}\in {\mathbb{S}}^2,$$

$$\tilde{Z}(\mathbf{x})=\frac{1}{L^{1/2}}\sum _{\ell =1}^L{Z}_{\ell }(\mathbf{x}),\mathbf{x}\in {\mathbb{S}}^2,$$

$$\nu (\mathbf{x})=0.6+0.3\times [1+\exp \{-10({\phi }_x-\pi /2)\}{]}^{-1}.$$

$$\nu (\mathbf{x})=0.6+0.3\times \exp \left\{-20({\phi }_x-\pi /2{)}^2\right\}.$$

$$P_n({\mathbf{x}}_1^{\top }{\mathbf{x}}_2)=\sum _{m=-n}^n{Y}_{n,m}({\mathbf{x}}_1)\overline{Y_{n,m}({\mathbf{x}}_2)},{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{S}}^2,$$

$$\sum _{i=1}^k\sum _{j=1}^k{a}_i{a}_{j}C({\mathbf{x}}_i,{\mathbf{x}}_j)$$

$${\int }_{{\mathbb{S}}^2}{P}_n({\mathbf{w}}^{\top }{\mathbf{x}}_1)P_k({\mathbf{w}}^{\top }{\mathbf{x}}_2)U(\mathrm{d}\mathbf{w})=\frac{{\delta }_{n,k}}{2n+1}{P}_n({\mathbf{x}}_1^{\top }{\mathbf{x}}_2),{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{S}}^2,$$

$${\int }_{{\mathbb{S}}^2}{P}_n({\mathbf{w}}^{\top }\mathbf{x})U(\mathrm{d}\mathbf{w})=0,\mathbf{x}\in {\mathbb{S}}^2.$$

$$E\{Z(\mathbf{x})\}=\sum _{n=0}^{\infty }{a}_n{\left\{\frac{b_n(\mathbf{x})(2n+1)}{a_n}\right\}}^{1/2}{\int }_{{\mathbb{S}}^2}{P}_n({\mathbf{w}}^{\top }\mathbf{x})U(\text{d}\mathbf{\omega }).$$

$$\mathrm{cov}\{Z({\mathbf{x}}_1),Z({\mathbf{x}}_2)\}=\sum _{n=0}^{\infty }\{b_n({\mathbf{x}}_1)b_n({\mathbf{x}}_2){\}}^{1/2}(2n+1){\int }_{{\mathbb{S}}^2}{P}_n({\mathbf{w}}^{\top }{\mathbf{x}}_1)P_n({\mathbf{w}}^{\top }{\mathbf{x}}_2)U(\text{d}\mathbf{\omega }),{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{S}}^2.$$

Full text

Keywords: Covariance function; Earth topography; Hausdorff dimension; Legendre polynomials; Random Fields; Schoenberg sequence.

1 Introduction

Realistic modelling of three-dimensional star-shaped particles has recently attracted great interest in a wide variety of scientific disciplines such as astronomy, material science and medicine, among others (Stoyan and Stoyan,, 1994). The Hausdorff (or fractal) dimension, a concept coming from the vocabulary of irregular shapes and self-similar patterns, is of paramount importance in the analysis of solid particles because it serves as a powerful mathematical tool to quantify the degree of smoothness or roughness of the surface of the particle. Indeed, recent studies suggest that fractal geometry can be useful for describing the topography of celestial bodies in the solar system (Kucinskas et al.,, 1992), the mechanisms of tumour growth and angiogenesis (Sedivy and Mader,, 1997) and the morphological features of sand particles (Zhou et al.,, 2017).

In the last decades, appealing stochastic models for star-shaped objects with constant Hausdorff dimension (also called monofractal objects) have been proposed by various authors, including Kent et al., (2000), Hobolth, (2003), Ziegel, (2013) and Hansen et al., (2015). Nevertheless, it is not uncommon to observe natural phenomena with multifractal patterns. For instance, Gagnon et al., (2006) argue that the Earth topography is multifractal and propose a simulation method in planar domains. Dellino and Liotino, (2002) study the relevance of multifractal volcanic ash particles from the eruptions of Monte Pilato-Rocche Rosse. The strong evidence of these complex geometries certainly calls for the development of new methodologies.

In this paper we propose an approach for modelling three-dimensional multifractal star-shaped particles. We follow the flexible model presented by Hansen et al., (2015), where the radial function of the particle is represented by mean of a Gaussian random field on the sphere. Specifically, the covariance function of the random field determines the Hausdorff dimension of the surface of the particle. While Hansen et al., (2015) focused on the study of isotropic random fields, and the mentioned monofractal objects, we extend their approach to the anisotropic case by admitting locally adaptive spatial dependencies. As a result, our model allows for a Hausdorff dimension varying from place to place on the boundary of the particle.

We additionally derive an algorithm for high resolution simulations, which relies on a spectral representation of the covariance function. We illustrate the use of our proposal through numerical examples, including a simulated version of the Earth topography, where the discrepancy between the regularities of the continents and the seafloor, as discussed by Gagnon et al., (2006), has been taken into account.

The paper is organized as follows. In Section 2 we review preliminary results about isotropic random fields on the sphere and introduce the definition of the fractal index. We also describe the modelling of star-shaped random particles in terms of its radial function. Section 3 proposes a general framework for escaping from isotropy, allowing for random fields with locally adaptive regularity properties. The connections between our proposal and a kernel-based method are also discussed. In addition, the associated simulation algorithm is derived. Section 4 illustrates numerical examples, while Section 5 concludes the paper with a discussion.

2 Background

2.1 Isotropic Gaussian random fields on the sphere

Let $\mathbb{S}^{2}=\{\bm{x}\in\mathbb{R}^{3}:\bm{x}^{\top}\bm{x}=1\}$ be the two-dimensional unit sphere, where ⊤ denotes the transpose operator, and consider a real-valued random field, $\{Z(\bm{x}):\bm{x}\in\mathbb{S}^{2}\}$, with finite second order moments. We assume $Z(\bm{x})$ to be Gaussian, i.e. for all $k\in\mathbb{N}$ and $\bm{x}_{1},\ldots,\bm{x}_{k}\in\mathbb{S}^{2}$, the random vector $\{Z(\bm{x}_{1}),\ldots,Z(\bm{x}_{k})\}^{\top}$ follows a multivariate Gaussian distribution. Thus, $Z(\bm{x})$ is completely characterized by its mean function, $\mu(\bm{x})=E\{Z(\bm{x})\}$, $\bm{x}\in\mathbb{S}^{2}$, and by its covariance function $C(\bm{x}_{1},\bm{x}_{2})={\rm cov}\{Z(\bm{x}_{1}),Z(\bm{x}_{2})\}$, $\bm{x}_{1},\bm{x}_{2}\in\mathbb{S}^{2}$.

Let us introduce the geodesic distance on $\mathbb{S}^{2}$, which is the main ingredient to define the property of isotropy of a random field. For two locations, $\bm{x}_{1}$ and $\bm{x}_{2}$ in $\mathbb{S}^{2}$, their geodesic distance is defined as $d(\bm{x}_{1},\bm{x}_{2})=\arccos\{\bm{x}_{1}^{\top}\bm{x}_{2}\}\in[0,\pi]$. We shall equivalently use $d(\bm{x}_{1},\bm{x}_{2})$ or the shortcut $d$ to denote the geodesic distance, whenever no confusion arises. Following Marinucci and Peccati, (2011), the random field is called (weakly) isotropic if it has a constant mean, and if its covariance function can be written as

[TABLE]

for some continuous function $K:[0,\pi]\rightarrow\mathbb{R}$. Thus, the covariance function just depends on the geodesic distance or, equivalently, the inner product. It is common to call $K$ the isotropic part of the covariance function $C$ (see, e.g., Guella and Menegatto,, 2018). For Gaussian random fields, isotropy also implies that the probability distribution of $\{Z(\bm{x}_{1}),\ldots,Z(\bm{x}_{k})\}^{\top}$ is invariant under the group of rotations on $\mathbb{S}^{2}$ (see Marinucci and Peccati,, 2011).

Note that for all $k\in\mathbb{N}$, and for all systems of points $\bm{x}_{1},\ldots,\bm{x}_{k}\in\mathbb{S}^{2}$ and constants $a_{1},\ldots,a_{k}\in\mathbb{R}$, we have that

[TABLE]

Condition (2.2), called semi positive definiteness, is a necessary and sufficient condition for a covariance function. In his pioneering paper, Schoenberg, (1942) showed that $C$ as in (2.1) is semi positive definite if, and only if, its isotropic part $K$ has a series representation in the form

[TABLE]

where ${P}_{n}$ denotes the Legendre polynomial of degree $n$, and $\{b_{n}:n\in\mathbb{N}_{0}\}$ is a sequence of nonnegative coefficients, such that $\sum_{n=0}^{\infty}b_{n}<\infty$. Classical inversion formulas yield

[TABLE]

Following terminology introduced by Ziegel, (2014), we refer to this sequence as a Schoenberg sequence.

The covariance function is often specified to belong to a parametric family whose members are known to be semi positive definite. For a thorough review of positive definite functions on spheres and a vast list of parametric families, we refer the reader to Gneiting, (2013).

2.2 Fractal index

The asymptotic behaviour of the isotropic part $K$ near zero, quantified by the fractal index, regulates the degree of smoothness or roughness of the sample paths of the associated isotropic random field. Recent findings in this direction can be found in Hitczenko and Stein, (2012), Hansen et al., (2015), Lang and Schwab, (2015), Guinness and Fuentes, (2016) and Clarke et al., (2018).

Formally, a random field with covariance function as in (2.1) has fractal index $\alpha>0$ if there exists a constant $c_{0}>0$ such that

[TABLE]

as $d\downarrow 0$. The fractal index exists for most parametric families of covariance functions, in which case it is always true that $0<\alpha\leq 2$, where $\alpha=2$ and $\alpha\rightarrow 0$ correspond respectively to extreme smoothness and roughness of the sample paths.

Abelian and Tauberian theorems (Bingham,, 1978; Malyarenko,, 2004) relate the asymptotic behaviour of $K$ near zero to that of its Schoenberg sequence near infinity. In consequence, the fractal index can be characterized in terms of the decay of the Schoenberg sequence. Indeed, a function $f:(0,\infty)\rightarrow(0,\infty)$ is called slowly varying at infinity if, for all $r>0$,

[TABLE]

Then, Malyarenko, (2004) shows that, for $0<\alpha<2$,

[TABLE]

as $d\downarrow 0$, if, and only if,

[TABLE]

as $n\rightarrow\infty$. The implication from (2.3) to (2.2) is called Abelian, whereas the converse is called Tauberian. The following example plays a fundamental role throughout the manuscript.

Example 2.1**.**

The Legendre-Matérn covariance function proposed by Guinness and Fuentes, (2016) is characterized by the Schoenberg sequence

[TABLE]

where $\tau$ and $\nu$ are positive parameters. While $\tau$ regulates the practical range (the distance at which the covariance function reaches certain threshold) of the random field, $\nu$ is capable of controlling the smoothness of the sample paths. Indeed, since

[TABLE]

as $n\rightarrow\infty$, we conclude from (2.3) that for $0<\nu<1$ the associated random field has fractal index $\alpha=2\nu$. This result can be extended to the limit case $\nu=1$, using additional technical arguments in Bingham, (1978). It is also worth noting that for $\nu>1$ the random field is mean square differentiable, in which case the fractal index is $\alpha=2$, and we refer the reader to Guinness and Fuentes, (2016) for details.

2.3 Star-shaped particles

We represent a particle as a compact set $Y\subset\mathbb{R}^{3}$, being star-shaped with respect to an interior point $\bm{o}$, that is, for all $\bm{x}\in Y$, the line segment from $\bm{o}$ to $\bm{x}$ is contained in $Y$. The set $Y$ is completely characterized by its radial function, defined as $Z(\bm{x})=\max\{r\geq 0:\bm{o}+r\bm{x}\in Y\}$, for $\bm{x}\in\mathbb{S}^{2}$. Accordingly, $Y$ can be represented as

[TABLE]

Adopting the framework proposed by Hansen et al., (2015), we model the radial function $Z(\bm{x})$ as a Gaussian random field on $\mathbb{S}^{2}$. Since $Z(\bm{x})$ has potentially negative values, it might be necessary replace $Z(\bm{x})$ with $Z_{c}(\bm{x})=\max\{c,Z(\bm{x})\}$ for some $c>0$. The interior point $\bm{o}$ can be assumed to be the origin.

The regularity of the surface of $Y$ can be mathematically quantified by the Hausdorff dimension (Adler,, 1981). We now turn to a formal definition in terms of ball coverings (Hansen et al.,, 2015). For $\epsilon>0$, an $\epsilon$-cover of $Y$ is a countable collection $\{B_{i}:i=1,2,\ldots\}$ of balls $B_{i}\subset\mathbb{R}^{3}$ of diameter $|B_{i}|$ less than or equal to $\epsilon$ that covers $Y$. Let

[TABLE]

be the $\eta$-dimensional Hausdorff measure of $Y$. The Hausdorff dimension of $Y$ is the unique nonnegative number $\eta_{0}$ such that $H^{\eta}(Y)=\infty$ if $\eta<\eta_{0}$ and $H^{\eta}(Y)=0$ if $\eta>\eta_{0}$.

Hansen et al., (2015) show that in the special case when $Z(\bm{x})$ is isotropic, the Hausdorff dimension of the particle is $\eta_{0}=3-\alpha/2$ almost surely, where $\alpha$ is the fractal index of $Z(\bm{x})$ defined in (2.1). Hence, it is always true that $2\leq\eta_{0}<3$. For sets describing traditional smooth shapes, the Hausdorff dimension matches the conventional topological dimension $\eta_{0}=2$. In contrast, as the Hausdorff dimension exceeds the topological dimension, the set becomes progressively irregular.

We see that when the radial function is characterized by an isotropic random field, the Hausdorff dimension of the surface of the object is constant. Particles with spatially varying Hausdorff dimension shall be obtained by relaxing the hypothesis of isotropy.

3 Anisotropic Random Fields on the Sphere

3.1 Locally varying Schoenberg sequences

As discussed in the previous section, the construction of versatile models for star-shaped particles relies on the appropriate specification of anisotropic random field models. We propose an approach to escape from isotropy, based on spatially adaptive Schoenberg sequences. More precisely, we consider the covariance function

[TABLE]

where $\{b_{n}(\bm{x}):n\in\mathbb{N}_{0}\}$ is a sequence of nonnegative functions on $\mathbb{S}^{2}$, with $\sum_{n=0}^{\infty}b_{n}(\bm{x})<\infty$, for each $\bm{x}\in\mathbb{S}^{2}$. We call this sequence of functions an adaptive Schoenberg sequence.

It is straightforward to show that (3.1) yields a semi positive definite function (see Appendix A for a proper justification). Moreover, (3.1) is still positive definite if we replace $\left\{b_{n}(\bm{x}_{1})b_{n}(\bm{x}_{2})\right\}^{1/2}$ with any positive definite function on $\mathbb{S}^{2}\times\mathbb{S}^{2}$. However, we shall see that (3.1) is sufficiently general to achieve models with flexible fractal properties.

This construction is convenient because when both $\bm{x}_{1}$ and $\bm{x}_{2}$ are near some fixed location $\bm{x}_{0}\in\mathbb{S}^{2}$, and $b_{n}(\cdot)$ is a sufficiently smooth function, for each $n\in\mathbb{N}_{0}$, we have

[TABLE]

As a result, covariance functions with arbitrary locally isotropic behaviours are possible, making this class of models attractive. In particular, the fractal index can vary from place to place on the spherical surface, in which case we say that the random field is multifractal.

Example 3.1**.**

A natural extension of the Legendre-Matérn model (2.4) is given by the adaptive Schoenberg sequence

[TABLE]

where $\nu(\cdot)$ is a positive function that controls the local regularity of the sample paths. In what follows, we call (3.2) the adaptive Legendre-Matérn model. Of course, $\tau$ may also be a locally adaptive parameter, however, we essentially focus on the fractal properties of the model.

3.2 Connections with a kernel-based method

We now show that the proposed model, based on adaptive Schoenberg sequences, offers great flexibility. Actually, it is capable of emulating covariance functions coming from kernel-based methods, widely used in the spatial analysis literature (Fuentes,, 2002; Nott and Dunsmuir,, 2002).

In a kernel-based approach, the random field $Z(\bm{x})$ on $\mathbb{S}^{2}$ is represented as a spatially weighted combination of isotropic random fields. This strategy permits to model dissimilar local dependency structures in different spatial zones. Let $D_{1},\ldots,D_{m}$ be a collection of subregions that cover $\mathbb{S}^{2}$, and let $\lambda_{j}(\bm{x})$ be a positive kernel function centered at the centroid of $D_{j}$, for all $j=1,\ldots,m$ (see Schreiner,, 1997). Consider the random field

[TABLE]

where $Z_{1}(\bm{x}),\ldots,Z_{m}(\bm{x})$ is a collection of independent isotropic random fields on $\mathbb{S}^{2}$. Suppose that the covariance function of $Z_{j}(\bm{x})$ has isotropic part $K_{j}$. Thus, $Z(\bm{x})$ has the following covariance function

[TABLE]

Since $K_{j}$ admits a representation of the form $K_{j}(d)=\sum_{n=0}^{\infty}p_{n,j}P_{n}(\cos d)$, for each $j=1,\ldots,m$, where $\{p_{n,j}:n\in\mathbb{N}_{0}\}$ is the associated $j$th Schoenberg sequence, the covariance function (3.1) can be written as

[TABLE]

We observe that (3.2) can be obtained as the sum of $m$ covariance functions of the form (3.1), where the $j$th adaptive Schoenberg sequence of type (3.1) can be taken as $b_{n,j}(\bm{x})=\lambda_{j}^{2}(\bm{x})p_{n,j}$.

3.3 Simulation algorithm

This section presents a method for simulating anisotropic Gaussian random fields on $\mathbb{S}^{2}$. The representation (3.1) allows for an immediate simulation procedure based on the adaptive Schoenberg sequence. The following proposition is crucial to develop the simulation algorithm.

Proposition 3.1**.**

Let $\kappa$ be a discrete random variable with ${\rm pr}(\kappa=n)=a_{n}$, $n\in\mathbb{N}_{0}$, where $\{a_{n}:n\in\mathbb{N}_{0}\}$ is a probability mass sequence with a support containing that of the adaptive Schoenberg sequence $\{b_{n}(\bm{x}):n\in\mathbb{N}_{0}\}$, for all $\bm{x}\in\mathbb{S}^{2}$, and let $\bm{w}$ be a random vector uniformly distributed on $\mathbb{S}^{2}$. Suppose that $\kappa$ and $\bm{w}$ are independent. Then, the random field

[TABLE]

has mean $\mu(\bm{x})=\{a_{0}b_{0}(\bm{x})\}^{1/2}$, and its covariance function is given by (3.1).

This method appears as the spherical counterpart of the spectral turning bands simulation algorithm used in Euclidean spaces (see, e.g., Mantoglou and Wilson,, 1982 and Emery and Arroyo,, 2018). For a neater exposition, the proof of Proposition 3.1 is deferred to Appendix B.

Even though the random field ${Z}(\bm{x})$ in (3.1) has the predefined covariance function, it is clearly non-Gaussian distributed. By a central limit effect, an approximately Gaussian random field can be obtained by

[TABLE]

where ${Z}_{1}(\bm{x}),\ldots,{Z}_{L}(\bm{x})$ are $L$ independent copies from (3.1), and $L$ must be a large integer.

Basically, the method consists of an adequate weighted and rescaled combination of Legendre waves. The algorithm separates the choice of the adaptive Schoenberg sequence from the choice of the sequence $\{a_{n}:n\in\mathbb{N}_{0}\}$, which is equivalent to an important sampling technique. We observe that location and dispersion parameters can also be added in order to control the mean and variance of the sample paths. As discussed in Emery and Arroyo, (2018), the process time of this algorithm is proportional to the number $L$ of copies and to the number of target locations, and it turns to be considerably fast.

4 Numerical Examples

4.1 Simulated particles

We now illustrate simulated particles from our proposal. Let $\varphi_{x}\in[0,\pi]$ be the polar angle of $\bm{x}\in\mathbb{S}^{2}$ in the spherical coordinates system. We consider the adaptive Legendre-Matérn model (3.2), with $\tau^{2}=0.1$ and two different structures for $\nu(\bm{x})$:

(A)

Increasing Hausdorff dimension, from South to North, according to

[TABLE]

(B)

Hausdorff dimension with high values at the poles and lower values at the Equator, according to

[TABLE]

For each case, the Hausdorff dimension varies in terms of the polar angle, approximately in the range $(2.1,2.4)$. We set $b_{0}(\bm{x})\equiv 0$ and incorporate a constant mean. We have also rescaled the covariance function to obtain a constant variance. The resulting radial function has mean $E\{Z(\bm{x})\}=100$ and variance ${\rm var}\{Z(\bm{x})\}=10$. The random variable $\kappa$ in Proposition 3.1 is chosen to follow a negative binomial distribution with parameters $r=1$ (number of failures) and $p=0.1$ (success probability). For this distribution, the probability mass function decays slowly to zero, i.e. it favors the occurrence of large values of $\kappa$, which is a desirable property in the sampling process. In our experiment, we have considered a spatial grid, in terms of the polar and the azimuthal angles, with $2000\times 1000$ target points. Figure 4.1 shows the simulated particles, with $L=300$, for the cases (A) and (B), where it is seen that the simulations match the theoretical features.

4.2 Earth topography: continents versus seafloor roughness

The aim of this section is to provide a simulated version of the Earth topography. We consider an extension of the monofractal simulation given by Hansen et al., (2015). In our analysis, we take into account the different degrees of roughness present in the continents and the seafloor (Gagnon et al.,, 2006).

The mean function $\mu(\bm{x})$ is estimated from online data sources (see, e.g., data outputs from the Joint Institute for the Study of the Atmosphere and Ocean, Seattle, United States), by using spherical harmonics regression, which is the natural basis for the spherical geometry (Marinucci and Peccati,, 2011), i.e. we consider the regressors $P_{n}^{m}(\sin\varphi_{x})\cos(m\theta_{x})$ and $P_{n}^{m}(\sin\varphi_{x})\sin(m\theta_{x})$, where $\varphi_{x}$ and $\theta_{x}$ are the polar and azimuthal angles in the spherical coordinates system, and $P_{n}^{m}$ denotes the associated Legendre polynomial, for $n\in\mathbb{N}_{0}$ and $m=0,\ldots,n$. We set $n=60$ in the estimation setting.

While the mean function captures the large scale variation of the Earth topography, we must incorporate the small scale variation. According to Gagnon et al., (2006), we distinguish between the Hausdorff dimension of the seafloor ($\eta_{0,S}=2.68$) and the continents ($\eta_{0,C}=2.46$). We then consider a dichotomic Schoenberg sequence from the adaptive Legendre-Matérn model. We rescale the covariance, and the resulting realization is chosen to have a constant standard deviation of approximately $2640$ meters, obtained empirically from the online resources mentioned above. Figure 4.2 depicts the simulated topography of planet Earth, on a grid of longitudes and latitudes with resolution of $1^{\circ}\times 1^{\circ}$, where we have used the same setting as in Section 4.1.

Figure 4.3 displays the Earth topography as a function of the longitude and latitude. It gives us a complete picture of the Earth map. The radial function (in kilometers) along the latitude $-30^{\circ}$ is also supplied. This example illustrates that continents are smoother than the seafloor, which is justified by tectonic forces and erosion-like processes (Gagnon et al.,, 2006).

5 Discussion

We have introduced a flexible framework for modelling and simulating three-dimensional multifractal star-shaped particles. The radial function of the object has been modelled by mean of an anisotropic Gaussian random field on the sphere, which is generated from a locally varying Schoenberg sequence. A simple adaptive version of the Legendre-Matérn covariance function has been proposed and employed to generate particles with place to place variable Hausdorff dimension. Our findings have been exemplified through numerical experiments, including an illustration of the Earth topography. Our approach may be used as a building block for more sophisticated models emulating the surface of celestial bodies.

A natural generalization of this work is to consider temporally dependent Schoenberg sequences, allowing for particles with a dynamically updated Hausdorff dimension. The findings of Berg and Porcu, (2017) might be useful here. Another interesting problem is the search for efficient methods for estimating locally variable Hausdorff dimensions. One challenging possibility is to adapt the tools reviewed in Gneiting et al., (2012). Advances in this direction also include the works of Anderes and Stein, (2011) and Ziegel, (2013). The extension to non-Gaussian particles can be tackled by using transformations of Gaussian random fields as in Xu and Genton, (2017). This topic is interesting because the probability distribution of the random field can influence some geometrical aspects of the particle. For instance, Hansen et al., (2015) show that Gamma-Lévy particles exhibit more pronounced spikes. From a computational viewpoint, one may be interested in proposing an improved version of the simulation algorithm by using parallel computing.

Finally, Gneiting, (2013) states, in his open Problem 15, that new methodologies involving anisotropic dependencies are also desirable in environmental and climatological phenomena (see also Hitczenko and Stein,, 2012 and Castruccio and Stein,, 2013). So we believe that our findings in Section 3 can be useful to develop new applications in various fields related to spatial analysis.

Appendices

A Positive definiteness of (3.1)

The set of spherical harmonic functions, $\{{Y}_{n,m}:n\in\mathbb{N}_{0},m=-n,\ldots,n\}$, form an orthogonal basis of the Hilbert space of square integrable functions on $\mathbb{S}^{2}$. Explicit expressions for the real and the imaginary parts of $Y_{m,m}$ has been used in Section 4. The addition theorem for spherical harmonic functions (Marinucci and Peccati,, 2011) establishes that

[TABLE]

where $\overline{c}$ denotes the complex conjugate of $c\in\mathbb{C}$. The semi positive definiteness of (3.1) is a direct consequence of the addition theorem. In fact, a straightforward calculation shows that

[TABLE]

where $|c|=(\overline{c}c)^{1/2}$ denotes the magnitude of $c$. The last expression is clearly nonnegative.

B Proof of Proposition 3.1

Before we state the proof of Proposition 3.1, we must introduce some properties of Legendre polynomials. Ziegel, (2014, Equation (2.4)) establishes that, for any $n,k\in\mathbb{N}_{0}$,

[TABLE]

where $U$ is the uniform probability measure on $\mathbb{S}^{2}$ and $\delta_{n,k}$ denotes the Kronecker delta. In particular, for any $n\in\mathbb{N}$, we have

[TABLE]

We now proceed with the proof of Proposition 3.1. Let ${Z}(\bm{x})$ be the random field defined in (3.1). The mean function is

[TABLE]

The property (B.2) implies that $E\{{Z}(\bm{x})\}=\{a_{0}b_{0}(\bm{x})\}^{1/2}$. On the other hand, the covariance function of ${Z}(\bm{x})$ is

[TABLE]

Using (B.1), we obtain the covariance function (3.1). The proof is completed.