Q Hua

A Fichet

A Wilkie

A Compact Representation for Fluorescent Spectral Data

Keywords:

a) Rendering of fluorescent materials using a reradiation matrix (27 552 B for the 3 pure reradiation spectra). (b) Rendering of fluorescent materials using the proposed Gaussian Mixture Model representation (529 B for all 3 pure re-radiation spectra).

8 0 (c) ∆E * 00 between (a) and (b).

Figure 1: In this scene, we use a near monochromatic light for illumination. While non-fluorescent materials only reflect shades of the light colour, fluorescent surfaces also reflect other wavelengths leading to the colourful appearance of the three fluorescent spheres in this scene.

In daylight environment, fluorescent effects are generally responsible for brighter and more vivid colours: weakly visible radiation are reemitted to wavelength our eyes are more sensitive to. To handle such spectra, we commonly use a re-radiation matrix (a). These triangular matrices have an important memory overhead. We propose to represent this data more efficiently with a Gaussian Mixture Model (b). While being compact, this representation is also precise with limited difference compared with using the full re-radiation matrix (c).

Introduction

Fluorescence is a common phenomenon both in natural elements and in chemical compounds. It has an important impact on how we perceive a scene, both in terms of hue, saturation, and luminance: a part of non-visible or near visible light is transferred into a different part of the visible range, so that the effect is widely used to enhance material brightness. In addition to this, it can also have the effect of making colour appear more vivid than reflective colour can normally look. While the importance of this effect for appearance modelling in computer graphics has been known for years, it has not been used in practice much yet, as adding it to graphics workflows -both on the side of user interfaces, as well as within rendering software proper -is a demanding problem.

One of the issues standing in the way of widespread adoption is that representing fluorescent effects consumes significant memory, and causes non-negligible computational overhead. To fully describe fluorescent spectra, one normally uses re-radiation matrices [START_REF] Donaldson | Spectrophotometry of fluorescent pigments[END_REF]: these represent emission and reemission efficiency for a given set of wavelengths. This requires storing a densely sampled 2-dimensional dataset for each fluorophore used in a scene. While explicitly storing this data is manageable for few fluorescent elements, this does not scale for fine-grained material definitions such as fluorescent textures.

In recent years, there has been an increasing interest in the effect in the research community. In one recent example, it was used to enlarge the colour gamut during spectral uplifting [JWH *19], in a fashion analogous to the use of optical brighteners in real materials. However, for lack of better alternatives, the method presented in that paper had to store a full re-radiation matrix for all pixels in a texture: which in turn made it incapable of working with more than a single fairly low-resolution input image.

In this paper, we propose to efficiently represent re-radiation matrix by using a Gaussian Mixture Model (GMM) to represent the re-radiation mechanism. This allows not only a low memory footprint, but also enables efficient importance sampling of wavelength shifting events.

Previous work

Spectral rendering Many rendering pipelines still rely on RGB or tristimulus values to compute light transport. This allows a compact representation for reflectance and illuminants, and fast computation. However, the resulting colours are not accurate and cannot be used for predictive rendering applications. Also, this type of rendering cannot simulate wavelength dependent effects such as dispersion or fluorescence. With the increasing computational power and the efficient multiplexing approach proposed by [WND*14], spectral rendering is becoming increasingly popular.

Bi-spectral rendering While not widely used, the concept of fluorescence in computer graphics is not new. Glassner [START_REF] Glassner | A Model for Fluorescence and Phosphorescence[END_REF] introduced fluorescence and phosphorescence to the computer graphics community and proposed a method to support those effects in a Whitted ray tracer. Wilkie et al. [START_REF] Wilkie | Combined Rendering of Polarization and Fluorescence Effects[END_REF] proposed a combination of rendering fluorescent and polarisation effects. Despite using a path tracer, they relied on a defined number of spectral bands, the wavelength not being a part of the Monte-Carlo integration. Mojzík et al. [START_REF] Mojzík | Handling Fluorescence in a Uni-directional Spectral Path Tracer[END_REF] introduced a method for both handling fluorescent medium in a Monte Carlo path tracer and making the wavelength domain a part of the integrand, and adapting Hero Wavelength Spectral Sampling [WND*14] to support fluorescence.

Color gamut enlargement One of the challenges when using a spectral renderer is the need to have spectrally defined assets. Assets are commonly defined in a tristimulus space (RGB, XYZ, etc.) and cannot be used directly as the input of a spectral renderer. To allow usage of RGB assets in a spectral renderer, one needs to "uplift" the original RGB data to a spectral form. There is previous work dealing with non-fluorescent spectral uplifting [OYH; JH19; MSHD15]. Recently Jung et al. [JWH*19] proposed an uplifting pipeline that additionally also uses fluorescent spectra to enhance the colour gamut. An important drawback of this bi-spectral uplifting method is its memory requirement: when dealing with texture uplifting, each pixel potentially needed to store a full reradiation matrix. However, this was not an intrinsic drawback of their method, but an indicator that this part of fluorescent rendering technology needed improvement.

Bi-spectral material models Wilkie et al. [START_REF] Wilkie | A Reflectance Model for Diffuse Fluorescent Surfaces[END_REF] proposed one of the early bi-spectral models adapted to current rendering techniques. This model uses a layered BRDF with a diffuse fluorescent component. Later, Hullin et al. [HHA*10] proposed an efficient acquisition setup to capture Bi-spectral Bidirectional Reflectance and Re-radiation Distribution Functions (BRRDF). They guided their acquisition by a Principle Component Analysis to lower the acquisition time.

Gaussian representations in rendering

The general idea of using Gaussian distributions in rendering is very powerful, as they can approximate the distribution of samples in any shape, given a large amount of sample numbers. A number of papers have experimented with Gaussian representations in the area of physically-based rendering. Jakob et al.

Background

In this section, we introduce the core mathematical and physical context used in the remainder of this article. We present the fluorescent effect, its necessity in physically-based rendering, and the typical representation of fluorescent data in a fluorescent capable renderer. We also review the basic notation of Gaussian mixture model and the formal framework of expectation maximisation, which is one of the most effective methods to train GMM in the machine learning field.

Fluorescence

Fluorescence is the effect where a "fluorophore" absorbs photons and re-emits those at a lower energy state. In the context of rendering, this means a colour shift when light is interacting with such a molecule. In contrast to phosphorescence, fluorescent effects occur within a very short period of time, and can be considered as instantaneous in the context of computer graphics.

Fluorescence produces interesting visual effects: as the human vision system is only sensitive to a comparatively narrow range of wavelengths, fluorescent materials often are used to increase visibility and colourfulness. Energy from the ultra-violet region, which is at best barely visible for the human eyes, is transferred to a bandwidth with much higher sensitivity. Such fluorescent elements appear brighter than the environment despite the decrease of energy between the absorption and the reemission that still happens.

A fluorescent element is sometimes characterised by an absorption and a reemission spectrum (see Fig. 2a). The first defines which wavelenghts lead to a reemission event and the amplitude define The synthetic representation of the average absorption and reemission spectra is often used to characterise a fluorophore. It assumes the shape of the reemission is similar regardless of the radiating wavelenght but scaled in intensity. However, Kasha's rule does not hold for may fluorophore hence the use of reradiation matrices (b). Such matrix provides precise information about the reemitted spectra given a radiating wavelength but have a significant memory footprint.

the intensity of the reemission. The reemission spectrum defines the average reemission accross all radiating wavelengths. The shift between the two central frequencies of those spectra is the Stokes' shift. These spectra then only scale in intensity depending on the radiating wavelength according to Kasha's rule.

When following the Kasha's rule, this reemission spectrum is considered of the same shape regardless of the radiating wavelength but scaled in intensity. Kasha's rule does not always strictly hold, however. This leads to the need for re-radiation matrices, which characterise, for each radiating wavelength, a distribution of reemission wavelengths (see Fig. 2a). While more flexible and suited for compute graphics application, this representation has an important memory overhead.

In the context of path tracing, a renormalised reemission spectra can be seen as the PDF of wavelength shifting for a given incoming radiation λ i . In reverse, a renormalised absorption can be seen as the PDF of absorbing a wavelength λ i for a given outgoing radiation λo.

Gaussian Mixture Model

A Gaussian Mixture Model (GMM) is a parametric probabilistic model representing a dataset by a linear superposition of Gaussian functions. This linear superposition is represented by:

p(x) = K ∑ k=1 π k N (x | µ k , Σ k ) (1)
where:

p(x) is the density of mixture of Gaussians, N (x | µ k , Σ k ) is a single Gaussian density, π k is the mixing coefficient (weight) of each Gaussian, K
is the number of Gaussians. In this article, we use a 2-dimensional multivariate Gaussian distributions for each component:

N (x | µ, Σ) = 1 2π |Σ| e -1 2 (x-µ) T Σ -1 (x-µ) (2)
where: µ is the mean, Σ is the covariance matrix.

We model our fluorescent dataset as a mixture of Gaussian and determine their individual parameters to most accurately represent the input data.

Expectation Maximisation (EM) is a two-stage iterative algorithm to estimate the mixture model parameters. It tries to maximise likelihoods for parametric models. The E-step (expectation) performs estimation of the expected log-likelihood for completedata set. M-step (maximisation) maximses the expected complete log-likelihood. The E-step and M-step are run iteratively until the expectation converges. EM and its variants are widely used to fit parametric models. As detailed introduction of EM is beyond the scope of our work, we only review briefly the sketch of EM, for detailed study on EM please refer to [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF]. Vorba et al. [VKŠ*14] also show a detailed version of how to adapt EM to the rendering context.

Method overview

Our pipeline consists of two main stages (see Fig. 3):

A fitting phase were we use the Expectation Maximisation algorithm to represent the input re-radiation dataset as a Gaussian Mixture Model. We describe this part in Sec. 5. A rendering phase were we reconstruct the re-radiation on the fly and perform importance sampling based on the fitted model. We describe this part in Sec. 6.

The pipeline starts from a 2-dimensional fluorescent measurement dataset. To get a clearer representation for fitting, we first filter this dataset. Then we perform weighted Expectation Maximisiation to get a compacted representation of the distribution. We then load the fitted GMM representation in the renderer. For each fluorescent sampling event, we use conditional GMM sampling on the given wavelength. Depending on the path direction, we either sample an excitation wavelength or reemission wavelength. The sampled wavelength will then shift to reproduce the fluorescent effect.

Using a GMM to represent the original dataset drastically reduces the memory footprint. Instead of relying on a tabulated triangular matrix, we only need to store few parameters representing parameters of the mixture, the re-radiation is then evaluated analytically. This reduction of memory footprint is even more prominent when the rendering technique rely on random sampling, as this also removes the need for additional tabulated CDFs used to sample the re-radiation matrix.

Fitting phase

In the fitting phase, we apply parametric density estimation to our fluorescence data. Later, we show that we can reconstruct the reradiation matrix from the GMM on the fly with high accuracy, instead of using the raw tabulated data.

Filtering

Very often, measured re-radiation data contains high frequency noise. In our pipeline, we remove all data under a certain threshold to pre-filter the data prior to the fitting phase. We found this simple approach appropriate for all the datasets we have analysed. For more challenging cases, e.g., lower signal-to-noise with less pronounced fluorescence or less precise measurement device, a more advanced method such as Fourier analysis may be used to analyse and remove the noise without impacting the other steps of the pipeline.

Weighted Expectation Maximisation

The re-radiation matrix is a 2-dimensional dataset where each entry corresponds to an attenuation. When each value is divided by the integral of the matrix, it represents a tabulated probability distribution function, each entry is a weight of this distribution.

In the standard EM-algorithm, each element corresponds to an observation. To emulate this, we could repeat each entry according to its weight in the re-radiation matrix. From our experiments, this strategy works well for fitting the data with different weights. However, this technique requires multiplying the weight by a large number to accurately represent the original data. This implies dealing with huge observation sets and consequently, a high memory footprint.

A more suitable approach is to use the weighted EMalgorithm [GAFH16; VKŠ*14]. All non-null elements of the matrix are used as observation outcomes and weighted according to their respective values.

We apply the weighted EM algorithm to fit Gaussian Mixture Model with our datasets using the pomegranate Python library [Sch]. Each data entry is associated with the normalized reradiation weights. The details of the derivation on weighted EM can be found in appendix.

Re-radiation reconstruction

We can reconstruct the re-radiation matrix from the Gaussian Mixture Model found previously. However, the Gaussian Mixture gives a probability density function, hence its integral is unity. The reradiation integral does not hold this property. So, we scale the Gaussian mixture to reconstruct the re-radiation later on:

Φ (λ i , λo) = S • p (λ i , λo) (3)
where: S is the scaling factor to apply, Φ (λ i , λo) is the re-radiation function represented by a matrix in its tabulated form, p (λ i , λo)

is the probability density function of the GMM defined in Eqn. 1, λ i , λo are the wavelengths at a given radiating and reemitting event.

There are different valid strategies to retrieve the scaling factor S. In this paper, we propose two solutions.

Error minimisation The first strategy consists in a minimisation process. S is set as a parameter to minimise the error between the reference dataset and the reconstructed mixture.

S = arg min S ∑ Λi ∑ Λo Φ measured (λ i , λo) -S • p (λ i , λo) (4)
Such minimisation process involves a choice of a norm function. But, this choice has little influence when dealing with fluorescence: the data have low dynamic range. Therefore, we use a L 2 norm in the result section i.e. the mean square error between the reference dataset and the reconstruction defined in Eqn. 3.

By definition, such strategy reduces the average reconstruction error for a given value but, we show that in Sec.7 that is provides less accurate rendering.

Integral ratio

The second strategy consist in defining S as the ratio between the integral of the original dataset and the GMM.

S = Λ Φ measured (λ i , λo)dλ i dλo Λ p(λ i , λo)dλ i dλo (5)
This method ensures the reconstructed dataset has the same fluorescent albedo as the original dataset. While this technique does not reduce error for each entry, it provides better rendering results: our eyes are more sensitive to brightness variation than to slight color shift as we show in Sec.7.

Rendering Phase

During the rendering phase, when a fluorescent event occurs, we first need to sample an in-shifting or an out-shifting wavelength depending on the ray direction. Then, we need to evaluate the attenuation given a set of λ i and λo wavelengths. The sampling phase relies on conditional probability: either λ i or λo is known, and we have to sample a λo or a λ i according to the ray direction, i.e., an in-shifting or out-shifting event.

Conditional GMM

To sample a random wavelength shift after a fluorescence event, we need to "sample a slice" on the multivariate Gaussian mixture. With a fixed input wavelength λ i , we construct the conditional parameters for the resulting Gaussian mixture and generate the random variable as the output wavelength λo.

In our 2-dimensional dataset, we have:

λ = λ i λo , µ k = µ ka µ k b , Σ k = σ kaa σ k ab σ k ba σ k bb ( 6
)
where λ is a sample point, µ k and Σ k are the GMM parameters for the k th Gaussian.

We define T, the inversion of k th covariance matrix as:

T k = Σ -1 k = τ kaa τ k ab τ k ba τ k bb (7)
The conditional probability that given λ i , we shift to λo in each single component k is:

p k (λ i → λo) = N λo | µ k b|a , τ -1 k bb (8) with : µ k b|a = µ k b -τ -1 k bb τ k ba λ i -µ ka (9)
Then, we get the probability for the entire conditional GMM from a given λ i to λo:

p(λ i → λo) = K ∑ k=1 π (λi) k p k (λ i → λo) (10)
with : π

(λi) k = π k N λo | µ ka , σ kaa ∑ k N λo | µ ka , σ kaa . (11)
Finally, we generate the random variable based on the new conditional GMM model, we first choose which Gaussian component we will sample by the conditional mixing coefficients π (λi) k , then we generate each sample following normal distribution using the chosen k th Gaussian:

µ (λi) k = µ k b|a , σ (λi) k = τ -1 k bb . (12)
We have only shown the model for sampling the light path, but it is straightforward to symmetrically apply the model for sampling on an eye path.

Importance Sampling

When a ray hits a fluorescent element, two events can occur:

• the ray is reflected with the same wavelength as the radiating ray, there is no fluorescence event, • the ray is re-radiated to a different wavelength, there is a fluorescent event happening.

We first determine if a fluorescent event occurs based on the ratio between the main diagonal and the total reflectance for the ray radiating wavelength. If a fluorescent event occurs, we importance sample the re-radiation. We then use the Gaussian mixture to importance sample this fluorescent event (see Alg. 1). The figure shows the expected distribution of the reemission wavelength λo when the material is radiated at λ i = 460 nm. When a fluorescent event occur, the outgoing wavelength λo is sampled according to this conditional PDF. Notice the measurement noise present in the original data lowers the highest value of the PDF after normalisation when compared to the analytic conditional PDF derived from the GMM. Figure 4 show an example of the conditional GMM for a specific radiating wavelength λ i . Given a radiating wavelength, the model can generate a distribution of reemitted wavelengths. We use this distribution to importance sample the reemitted wavelength λo.

The probability for the transition between the set of radiating wavelength λ i and re-emitting wavelength λo is:

p shift (λ i → λo) =    Φ(λi,λi) rt (λi) if λ i = λo 1 - Φ(λi,λi) rt (λi) • ∑ K k=1 π (λi) k N (λo | µ (λi) k , σ (λi) k rt (λ i ) = Λ Φ(λ i , λr)dλr (13)
where: Φ(λ i , λ i ) is the non-fluorescent reflectance, rt (λ i ) is the total reflected energy given λ i .

In our renderer, we use Hero Wavelength Spectral Sampling (HWSS) [WND *14]. The HWSS technique multiplexes multiple wavelengths in a Monte-Carlo (MC) sample. It differs from a simple multiplexing approach by using a "Hero" component which is used to take all decisions on the directional component. This technique drastically reduces colour noise in contexts where there is a wavelength dependency on the directional decisions such as participating medium. In the case of fluorescence, we allow each of the wavelength in the HWSS vector to be independently importance sampled. If done directly, we can find non-zero samples with a zero probability which is incorrect for a MC integrator. Mojzík et al. [START_REF] Mojzík | Handling Fluorescence in a Uni-directional Spectral Path Tracer[END_REF] introduced a balance term to avoid this artefact. We use the same term in our implementation:

p j (λ i j → λo j ) = N ∏ k=1,k = j p shift (λ i k → λo k ) (14)
where: p j (λ i j → λo j ) is the probability of j th HWSS sample, p shift (λ i → λo) is the probability of shifting from λ i to λo, N is the size of the HWSS vector.

Algorithm 1: Random wavelength shifting. Given a λ i , we either have no fluorescent event, i.e. λ i = λo or, we sample a λo after a fluorescent event.

input : λ i output: λo and pdf := p shift (λ i → λo)

ξ 1 ← U(0, 1); p diag ← Φ(λi,λi) rt (λi) ; if ξ 1 ≤ p diag then //
s ← 0; for k ← 1 to K do s ← s + π (λi) k ; if s ≤ ξ 2 then break; end end / * Conditional sampling with rejection method * / do λo ← N µ (λi) k , σ (λi) k ; while λo ≤ λ i ; / * Eqn. 13 * / pdf ← (1 -p diag ) • ∑ K k=1 π (λi) k N (λo | µ (λi) k , σ (λi) k
; end

Implementation details

To evaluate the re-radiation, we store for each Gaussian of the mixture, its mean, covariance matrix and weight. We also need to store the scaling factor S. We precompute the determinant and the inverse of the covariance matrix before rendering to speed-up evaluation of the mixture. We store the diagonal in its original tabulated form without any alteration.

To efficiently perform the conditional wavelength shifting we also precompute two additional tabulated values to retrieve π

(λi) k or π (λo)
k : the sum of rows and the sum of columns of the Gaussian Mixture. In our implementation, we use a 1 nm sampling rate.

While additional values are precomputed and represent an additional memory footprint, this is still far less than the memory required when using a tabulated CDF which are needed when using only tabulated measurements.

We provide as supplemental material a repository containing the source code for fitting the measured data with Gaussian Mixture Model. The repository also contains a patch for ART renderer [START_REF]The Advanced Rendering Toolkit[END_REF] to use the computed GMM as input to reproduce the rendering showcased in this paper. It is accessible here: https://gitlab.com/huaq_ /compact-representation-fluorescence.

Results and Discussions

Looking at the entirety of our pipeline, we have:

• A reflectance spectrum corresponding to diagonal, the original non-fluorescent part. • A set of parameters representing the Gaussian mixture, which are:

-a weighting factor (1 element), -a mean vector (2 elements), -a covariance matrix (2x2 elements).

To evaluate the accuracy of our reconstruction, we propose to compare the results of the reconstructed re-radiation with the original tabulated data. We use the perceptual comparison on both rendering and reconstruction, with ∆E * 00 on rendering using a Monte-Carlo path tracer and the Mean Square Error between the original and reconstructed re-radiation matrices.

In this section, we review the results and discuss the strengths and weaknesses of our model, we especially emphasising on the reconstruction accuracy.

Re-radiation reconstruction accuracy

We use Mean Square Error (MSE) to evaluate the accuracy of the reconstructed matrices on a per-wavelength basis (see Fig. 5). For all datasets, the MSE significantly drops at 3 Gaussians. Adding additional Gaussians then have little effect on improving the reconstruction accuracy. The MSE shows Pink Post-It dataset is being more challenging to reconstruct: the re-radiation spectrum behaves less like a Gaussian distribution than the two other datasets. Using a different set of basis functions, which can account for skewness would probably improve the accuracy for this specific dataset.

Overall, our method provides good accuracy. Even with a small number of Gaussians, we fit the original datasets well, with a maximum MSE at 1.73 • 10 -5 when using the integral ratio for computing the weighting factor. Gaussians in two different lighting conditions. We compare with two different strategies to scale the GMM for re-radiation reconstruction: minimizing the MSE and integral ratio. The original datasets are measured with the excitation wavelength ranging from 300 to 700 nm and reemission ranging from 380 to 700 nm sampled each at 10 nm. The original dataset are stored in 49 by 41 matrices. To store the pure re-radiation spectrum, each of the original dataset needs 1148 elements. The GMM representation requires 7 elements per Gaussian and a weighting factor to represent the pure-re-radiation. Even with a low number of Gaussian, we get an accurate reconstruction with a minimal memory footprint. Using the integral ratio for the weighting factor gives better results despite the higher MSE (see Fig. 5): our vision system better notice a variation of brightness than a slight wavelength difference. Using the integral ratio better conserves the overall light intensity coming from re-radiation while MSE minimisation better conserves the accuracy on a per-wavelength basis.

Rendering accuracy

To evaluate the rendering accuracy, we use CIE-∆E * 00 difference between the reference datasets (tabulated values) and the reconstructed re-radiations using our proposed method with a varying number of Gaussians on renderings using a path tracer. We use three different datasets: measured re-radiation for pink, yellow and green Post-It. These datasets come from the ART open source renderer [START_REF]The Advanced Rendering Toolkit[END_REF] and originate from Labsphere Inc. according to the metadata included in the measurement archive † .

The lighting conditions have a significant impact on the error. The fluorescent nature of the dataset makes an objective evaluation of the rendering accuracy difficult. We choose to use two drastically different lighting environments: a daylight scene using a skydome from Hosek and Wilkie [START_REF] Hosek | An Analytic Model for Full Spectral Sky-Dome Radiance[END_REF] (this model is appropriate to use here as it contains an ultra-violet component) and a near monochromatic light radiating at 450 nm.

Our technique performs well even with a limited number of Gaussians (see Fig. 6). ∆E * 00 are lower than 4 in a daylight environment. Even for challenging cases where most of the energy falls into the re-radiation, with near monochromatic 450 nm illumination, the difference remains low for green and yellow Post-It dataset. For the pink Post-It dataset, the difference is more visible in the near monochromatic lighting. This is expected as this dataset is more challenging to fit with Gaussians.

Weighting factor

We evaluated two strategies for computing the weighting factor utilised to scale the GMM to reconstruct the re-radiation matrices: minimisation process on the MSE between the original dataset and the scaled GMM and integral ratio (see Sec. 5.3).

As expected, with the minimisation process, we obtain a lower MSE than using the integral ratio. Nevertheless, a lower MSE does not imply better rendering accuracy. In fact, we observed the opposite. At first glance, this might seem surprising. However, the ∆E * 00 evaluates visible difference for the human vision system. Our eyes are more sensitive to variation of intensity than to a slight shift in the spectral domain. Using the integral ratio ensures having a similar luminous quantity falling into fluorescence. Note that we did not use any perceptual weighting for computing this ratio. This gives us the insight that the weighting factor could also be optimised to minimise the ∆E * 00 given a specific illuminant. The strategy for computing the scaling factor depends on the target application. The integral ratio is better suited for perceptual applications like rendering with tone mapping. If the accuracy of the spectral reconstruction of re-radiation is needed, then the minimisation of MSE is a better choice.

Memory footprint analysis

Our method is memory efficient and works with only a fraction of the memory required when using raw tabulated data. † The ART measured bispectral data are available from ART_Resources/SpectralData/FluorescentPostIt.ark The measured re-radiation dataset we have used are sampled at a 10nm rate, for radiating wavelength ranging from 300nm to 800nm and reflected wavelength ranging from 380 nm to 800 nm. They are stored in 41 × 49 matrices. In the most favourable case, we only have to store the upper triangular part of the re-radiation due to the nature of fluorescence shifting wavelengths to a lower energy state. If extreme high accuracy is desired, each matrix can be interpolated to a 2-dimensional matrix sampled at a 1 nm rate; this matrix uses 98000 Bytes.

Apart from doing away with having to use the raw re-radiation matrix data, in practice, effective sampling of wavelength shifting is desirable. This requires additional memory space: additional tabulated CDF are computed at the initialisation stage. This triples the memory footprint when exclusively using tabulated data.

With our method, we only need a small set of Gaussian parameters to both evaluate on the fly the re-radiation and to importance sample wavelength shifting. If rendering accuracy is desired, our method can be used solely to importance sample while sparing the storage of tabulated CDF.

We have shown in Fig. 6 our fitting provides a good reconstruction of the tabulated re-radiation matrix with only a few Gaussians. Assume we use 4 Gaussian for fitting, while the tabulated data are sampled at 1 nm rate, the memory difference of one fluorescent material between the tabulated method and our method is approximate 2985776 Bytes.

Fluorescent textures

We can use our method along with fluorescent texture. In [JWH*19] they use such kind of textures to uplift RGB textures. This technique uses fluorescence on some pixels to enlarge the colour gamut. The generated textures can be very large in mem-ory due to the potential large number of re-radiation matrices. For a 256 × 256 texture, where every pixel is fluorescent, it could require up to 3 × 256 × 256 re-radiation matrices (one used for the attenuation and the two others for precomputed CDFs).

We applied our method to compress spectral uplifted textures from [JWH *19]. Our method accurately renders the fluorescent texture with a ∆E * 00 ≤ 8 between the reference image and our method (see Fig. 7). The texture is 32 × 32 pixel with 255 fluorescent pixels after the uplifting. The reradiation matrices sizes are 500 by 500 elements. The memory consumption for the diagonal is 4.10 MB and an additional 254.49 MB for the re-radiation. In this rendering, we used 2 Gaussians per fluorescent pixel lowering the re-radiation memory footprint down to 30.60 KB. The total memory footprint of the fully tabulated bi-spectral 32 × 32 texture is 258.59 MB compared to the 4.13 MB needed by when using GMM for the re-radiation. This makes it 62.67 times more memory efficient for this specific dataset. This does not even account for the additional memory overhead of tabulated CDFs desirable for efficient importance sampling of the tabulated dataset. This ratio is even higher when increasing the resolution of the texture implying a larger number of fluorescent pixels.

Conclusion and Future Work

We have presented a compact and efficient method for compressing re-radiation matrices. This allows to drastically lower the memory footprint needed to store a fluorescent spectrum. The memory reduction is useful in the case of fluorescent texture is present, as the large fluorescent texture has a high memory requirement. Our proposed method also enables efficient importance sampling.

In future work, we want to develop a parametric re-radiation model suited for edition. The presented method, while being memory efficient and accurate, lacks the ability for an artist to intuitively edit fluorescence. We also want to evaluate additional metrics for computing the weighting factor used in re-radiation reconstruction depending on the targeted application: visual accuracy or absolute re-radiation difference.

Although we show that Gaussian Mixture can be used to compact fluorescence, due to the limitation of public fluorescent dataset in the graphics community, it is hard to say which parametric model is the best to suit the fluorescence. Add some skewness distribution can be an interesting work to the specific pink fluorescent dataset. The access to more fluorescent data and make them usable in computer graphics will be beneficial to the future work in fluorescent representation.

There is also possibility to try out other fitting methods for the mixture, such as using MAP-EM or non-linear fitting framework with different loss functions. We tried variational bayesian inference for the fitting of the mixture, however, the rendering results can not give us consistent observation. We need more research on that to have a decisive answer on which fitting method is best to the fluorescence data. For the particular problem in color shift in Fig. 6, we think different fitting method might help us with the color shift, but another more potential way to improve that will be to design a human perception metric instead of any distance based metrics.

  c 2021 The Author(s) Eurographics Proceedings c 2021 The Eurographics Association.

  [START_REF] Jakob | Progressive Expectation-Maximization for Hierarchical Volumetric Photon Mapping[END_REF] introduced how to represent radiance using Gaussian mixture and apply for volumetric media. Herholz et al. [HEV*] and Vorba et al. [VKŠ*14] showed how Gaussian mixture can help with importance sampling in the path guiding framework. Herholz et al. [HES*18] used the Gaussian Mixture to sample analytic BRDF models and measured BRDF data.

  A fluorescent element is often characterised by its absorption spectrum and its reemission spectrum. The offset between the two spectra is called the Stokes shift. Example of a re-radiation matrix. This matrix gives the reemission spectra for every radiating λ i wavelength. We removed the diagonal (λ i = λo) to emphasise on the fluorescent effect.

Figure 2 :

 2 Figure2: (a) The synthetic representation of the average absorption and reemission spectra is often used to characterise a fluorophore. It assumes the shape of the reemission is similar regardless of the radiating wavelenght but scaled in intensity. However, Kasha's rule does not hold for may fluorophore hence the use of reradiation matrices (b). Such matrix provides precise information about the reemitted spectra given a radiating wavelength but have a significant memory footprint.

Figure 3 :

 3 Figure3: Representation of our pipeline. We model the input reradiation with a GMM in the fitting phase and then use it in the renderer. In the rendering phase, we importance sample the new wavelength according to the GMM.

Figure 4 :

 4 Figure4: Example of a conditional PDF derived from the GMM model compared to the unfiltered yellow fluorescent dataset. The figure shows the expected distribution of the reemission wavelength λo when the material is radiated at λ i = 460 nm. When a fluorescent event occur, the outgoing wavelength λo is sampled according to this conditional PDF. Notice the measurement noise present in the original data lowers the highest value of the PDF after normalisation when compared to the analytic conditional PDF derived from the GMM.

Figure 5 :

 5 Figure5: Mean squared error between the measured datasets and the reconstruction with different number of Gaussians and different weighting factors. We use the same mixture for both techniques to compute the scaling factor. From 3 Gaussians, we observe a significant drop in the reconstruction error using the MSE metric. Adding additional Gaussians does not significantly reduce MSE.

Figure 6 :

 6 Figure6: Comparison of reference datasets (tabulated re-radiation) with the proposed GMM representation using a different number of Gaussians in two different lighting conditions. We compare with two different strategies to scale the GMM for re-radiation reconstruction: minimizing the MSE and integral ratio. The original datasets are measured with the excitation wavelength ranging from 300 to 700 nm and reemission ranging from 380 to 700 nm sampled each at 10 nm. The original dataset are stored in 49 by 41 matrices. To store the pure re-radiation spectrum, each of the original dataset needs 1148 elements. The GMM representation requires 7 elements per Gaussian and a weighting factor to represent the pure-re-radiation. Even with a low number of Gaussian, we get an accurate reconstruction with a minimal memory footprint. Using the integral ratio for the weighting factor gives better results despite the higher MSE (see Fig.5): our vision system better notice a variation of brightness than a slight wavelength difference. Using the integral ratio better conserves the overall light intensity coming from re-radiation while MSE minimisation better conserves the accuracy on a per-wavelength basis.

  (a) Bi-spectral tabulated texture.(b) The same bispectral using GMM for re-radiation.

  ∆E * 00 difference between (a) and (b)

Figure 7 :

 7 Figure 7: With our method, we can drastically reduce the memory footprint of bi-spectral textures. We use the bi-spectral uplifting technique from Jung et al. [JWH*19] to generate a bi-spectral texture from an RGB image. Notice that the bottom part of both images has visible fluorescence. The reference image and our method are rendered with 60k SPP in the same base renderer except for the fluorescent sampling process.

c 2021 The Author(s) Eurographics Proceedings c 2021 The Eurographics Association.

Acknowledgements

We acknowledge funding by the Czech Science Foundation grant 19-07626S, the EU H2020-MSCA-ITN-2020 grant number 956585 (PRIME), and the Charles University grant number SVV-260588.

Appendix A:

The deviation of weighted EM used in the paper.

Given an observed data set X = {x 1 , x 2 , ..., xn}, for each observation in X, there is a corresponding latent variable in the set Z = {z 1 , z 2 , ..., zn} which helps to model dependencies between variables. The data set X, Z is seen as a complete data set.

When using EM to fit a mixture of K Gaussians, the latent variable set Z is a set of N-dimensional vector. Each component zn is a K-dimensional binary random variable, where z nk denotes the k th element of zn. z nk ∈ {0, 1}, ∑ k z nk = 1. This latent variable vector represents which Gaussian component generated our sample xn with probability p (z nk = 1) = π k .

Assume we have a dataset X = {x 1 , ..., xn}, the corresponding weight vector for each sample in X is W = {w 1 , ..., wn}. The mixture of Gaussian is:

The responsibility function is thus:

The means, covariances and mixture coefficients used to maximise the expected complete-data log-likelihood are:


Jean-Claude Léon
email: jean-claude.leon@hmg.inpg.fr

Philippe Maurice Marin
email: philippe.marin2@hmg.inpg.fr

Gilles Foucault
email: gilles.foucault.1@ens.etsmtl.ca

J.-C Léon

Operators and criteria for integrating FEA in the design workflow: Toward a multi-resolution mechanical model

Keywords: shape simplification, mesh, polyhedral model, mechanical criterion, finite element accuracy

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Design models are used by all the actors of the design process and therefore contain numerous details. These models are often too refined for mechanical analyses and their direct use would generate too many finite elements. The adaptation of the model shape needs the removal of its details when their presence has either no or limited effects on its mechanical behavior while requiring an important local mesh density. Examples of these details include fillets, but also detailed entities such as holes, small blocks, etc.

Various software make it possible to automate this step partially. Several categories of approaches have been proposed to solve the problems involved by the preparation of F.E. models from CAD data. A first one addresses configurations were small features must be removed to get the geometric model more compatible with the size of the F.E. required [START_REF] Mobley | An object oriented approach to geometry defeaturing for Finite Element Meshing[END_REF][START_REF] Dabke | Using features to support finite element idealisation[END_REF][START_REF] Fran¸cois | 3D automatic remeshing applied to model modification[END_REF]. These approaches are strongly dependent on the modelling history of the part and work on the construction tree of the object and the removal of user-selected features. A second one starts with a polyhedral model of the part [START_REF] Schroeder | Decimation of triangle meshes[END_REF][START_REF] Cohen | Efficient model simplification with global error bounds[END_REF][START_REF] Fine | Automated generation of FEA models though idealization operators[END_REF]. In order to simplify the model, different adaptation functions work on the initial polyhedral model. They combine decimation process and removal of topological details. Another category of approaches is characterized by idealization treatments. Such operations are often required to transform a volume into an open surface to model a plate behaviour. Similar operations hold for transforming a volume feature into a line to model a beam behaviour of the structure.

The accuracy of F.E. computations is one of the main concerns of the users. The sources of the errors are multiple, errors of discretization, uncertainty about the boundary conditions and the behaviour law of the constitutive material, simplification of the shape, ... The quality of F.E. computation can be strongly influenced by the simplifications carried out on the shape. Appropriately choosing and monitoring these simplifications is therefore of primary importance. When the preparation of the model is manual, its quality depends on the engineer's know-how. For an automatic simplification process, the monitoring process uses geometrical criteria, curve, size [START_REF] Véron | Static polyhedron simplification using error measurements[END_REF]... In a priori step, geometrical criteria related to the mechanical properties of the problem can be added, variation of mass, volume, sections, centre of inertia [START_REF] Foucault | Mechanical criteria for the preparation of finite element models[END_REF]... A posteriori indicators can be used also and adaptive simplification process can be performed to define the most suited simplified model for each analysis case.

Section 2 presents the existing polyhedral model simplification algorithms we use to automatically prepare F.E. models. In section 3, different implemented criteria are listed and some examples illustrate their efficiency.

Simplification operators

Our approach uses an intermediate polyhedral model of the object. Using such a representation, we can integrate data from CAD models, pre-existing F.E. meshes or 3D scans. algorithm. The simplification process is based on an iterative vertex removal. Figure 1 shows the result of a simplification process.

Mechanical criteria

During a simplification process, a priori criteria are geometric ones but they can be related to mechanical property variations of mass, volume, section, centre of inertia. These indicators have been developed with our decimation process. We could provide the user either macro-scale information over the whole object or micro-scale information on the smallest possible entity of the geometrical model. Figure 2 shows two examples of such an indicator. A priori criteria are an interesting information for the users but they cannot quantify directly the real influence of geometric simplifications on a F.E. simulation. For example, the errors generated by a hole removal will depend on the dimension of the hole but also of its location over the component. This form feature can be in an area involving either low or high stresses. To take into account its position, a mechanical criterion needs information about the highly stressed areas and hence a sketch of the F.E.A. results. To gain an idea of these results requires to set up an a posteriori process. After a simulation on the simplified part, mechanical criteria are computed for evaluate the influence of each suppressed feature. Based on criteria, the engineer can validate the quality of his (resp. her) F.E. results. These criteria can also used in an adaptive process of simplification to refine the simplified part by adding some of the suppressed features according to their mechanical influence.

Such a process of simplification was used in [START_REF] Véron | Shape preserving polyhedral simplification with bounded error[END_REF]. The criterion used was the error of discretization on the simplified problem and the map of sizes of an optimal mesh for this problem. The program removes details if the size of the detail is lower than that given in this area by the map of sizes. Another criterion can be estimation of the influence of each simplification on the strain energy variation [START_REF] Marin | Influence of geometrical simplifications on the accuracy of finite element computation[END_REF]. In the framework of stationary linear problems, this variation is given by eq. ( 1) where σ 2 are the stresses of the simplified model, U 1 the displacements of the initial model and n the outward unit normal. Obviously, these displacements are unknown and can be estimated by local computations around the removed feature. On Fig. 3, we show an example of static problem that illustrates the efficiency of our criterion. The simplified model has been obtained by suppressing four features. In order to evaluate the efficiency of the criterion, we compute the solution and the strain energy on the initial part but also on the simplified one. We can see on table 1 that this criterion produces a good estimation of the influence of each simplification. During an adaptive process, the user computes the solution and the influence indicator on the simplified model. In this example, if the user wants to bound the accuracy to 1%, the simplified model must be refined by adding features 3 and 4. The design process of complex structures needs the resolution of multiple mechanical analyses. Each analysis requires the generation of an adequate simulation model. To define easily such models and integrate F.E.A. in the design workflow, it is necessary to set up efficient simplification operators and criteria for evaluating and validating the simplification process. In this paper, we have presented a set of treatments contributing to this integration. Developments are in progress to add other mechanical criteria, like a priori stiffness variation, a posteriori indicator for modal analysis, ... Work is also needed to automate the adaptive process using our a posteriori indicator.

All these treatments help transform CAD models from pure into multiresolution mechanical models efficient for mechanical simulation preparation.

Fig. 1 .

 1 Fig. 1. Scan of a crankcase (courtesy Tomoadour-TurboMéca), generation of a F.E. mesh (total time ˜2h), 427132 faces in the initial model and 11924 in the F.E. mesh.

Fig. 2 .

 2 Fig. 2. Two examples of a priori criteria, variation of volume on the left part and variation of sections on the right part.

f eature boundary σ 2 . 1 )Fig. 3 .

 213 Fig. 3. Initial CAD model on the left and associated simplified model for analysis

Table 1 .

 1 For each feature simplification, comparison between the real strain energy variation and our estimation of this variation.

	A = real variation B = Influence indicator efficiency A/B
	feature 1 -0.014 %	-0.009 %	1.56 %
	feature 2 0.27 %	0.26 %	0.96 %
	feature 3 -8.4 %	-8.3 %	1.01 %
	feature 4 16 %	18 %	0.89 %

Acknowledgement. Future developments will be conducted in the AIM@SHAPE NoE framework. This Network of Excellence involves fourteen research teams about Shapes and associated semantic and it is supported by the European Community.


Neha Aswal
email: nehaaswal96@gmail.com

Eshwar Kuncham

Subhamoy Sen
email: subhamoy@iitmandi.ac.in

Laurent Mevel
email: laurent.mevel@inria.fr

Robust Interacting Particle-Kalman Filter based structural damage estimation using dynamic strain measurements under non-stationary excitation -an experimental study

Keywords: Structural Health Monitoring, Interacting Filters, Particle Filter, Kalman Filter, Damage Detection

Sensor types and their positioning is a major factor in structural health monitoring (SHM) to ensure certainty in estimation. While acceleration has predominantly been employed for damage detection, they are known to be costly and not frame invariant (except for moderately accurate GPS based accelerometers). A thorough monitoring of a real life structure requires dense instrumentation which might become expensive with costly sensor types. Further, damages mostly occur at rare events, like seismic base excitation, for which typical accelerometers are not proper. This study employs strain as a cheaper alternative for damage sensitive measurement that is also frame invariant. An interacting filtering approach with particle and Kalman filters is employed that estimates structural health from measured dynamic strains. Further to account for extreme non-stationary events like seismic excitation, robustness against uncertain inputs is induced in the filtering environment following an output injection approach. The proposed algorithm is tested on a seven story-one bay frame model and a real experimental beam structure.

INTRODUCTION

The objective of this paper is to monitor structures, subjected to non-stationary input forces through the response measurements, recorded using a network of sensors. To ascertain safety in structures, damages due to strong forces or extreme service conditions should be detected immediately after their occurrence. Traditionally, structural health monitoring (SHM) research employs deterministic approaches for real time structural damage detection with no consideration about the uncertainties originating from unavoidable model inaccuracies, sensor noises and unknown external disturbances [START_REF] Charles R Farrar | Excitation methods for bridge structures[END_REF][START_REF] Cawley | The location of defects in structures from measurements of natural frequencies[END_REF][3]. This limits the utility of the methods for real field applications. Bayesian filtering technique has been proved to be an efficient alternative that can promptly detect anomaly in structures in the presence of such uncertainties [4][5][START_REF] Zghal | Modal parameter estimation using interacting kalman filter[END_REF][START_REF] Sen | Seismic-induced damage detection through parallel force and parameter estimation using an improved interacting particle-kalman filter[END_REF][START_REF] Sen | Online structural damage identification technique using constrained dual extended kalman filter[END_REF].

Efficiency of SHM approaches depends on the density and positioning of the measurement locations. Nevertheless, cost of sensors and the accessibility of the location where the sensors are supposed to be placed pose the major challenge in real life application. In reality, following the financial constraints, the measured data are often sparse compromising the precision in system health estimation. However, employing a suitably selected numerical model (typically Finite Element Models (FEM)), the auxiliary correlation information between these sparse and non-collocated measured data, obtained from the conceptual models, can additionally be incorporated [START_REF] Sen | Online structural damage identification technique using constrained dual extended kalman filter[END_REF]. While this model-based approach benefits by reducing the estimation uncertainty obvious for sparse dataset, it demands a precise numerical model. The involvement of a numerical model in estimation inherently introduces modelling uncertainty besides the uncertainty due to measurement noises. Bayesian filtering based techniques can simultaneously handle model and measurement uncertainties. Structural health is further parameterized in this attempt and subsequently estimated typically employing acceleration as measurement. Yet, the acceleration sensors are costly and their placement on the nodes (structural joints) are not always easy due to the accessibility reasons. Also, under the base excitation conditions, the accelerometer fails to provide a reference independent response as measurement. Due to these practical complexities involved with acceleration sensors, frame invariant and comparatively cheaper strain gauges can be considered as an alternative for real life SHM problems. However, this demands a numerical model to map the system states to corresponding strain measurement.

Also, for filtering-based SHM approaches typical assumption on stationarity and/or Gaussianity for ambient forcing does not hold true for extreme loading conditions, like, earthquakes or high wind/waves. Eventually, with a forced idealization about the input, the algorithms may perform poorly for forcing scenarios that are not Gaussian and stationary. Present study discusses one such possibility in which the structural damage is required to be estimated under non-stationary forcing. To induce robustness against input forcing, current study takes basis of [START_REF] Zhang | State estimation for stochastic time varying systems with disturbance rejection[END_REF]'s work in which an output injection approach is adopted in order to reject the input. The system estimation is performed using an interacting filtering strategy in which a Kalman filter estimates the system states while a Particle filter monitors the variation in health parameters.

Further, since this study adopts strain as measurement, a model is required to map system states (typically displacement and velocity) to corresponding strain. The strain-displacement relationship is therefore developed through an FEM approach taking basis on the shape functions. The proposal is validated on a seven story-one bay frame structure involving 21 members and 84 degreesof-freedom(dof ). Further, real experiments are performed on a cantilever beam to calibrate its material properties and detect the presence of damage in the beam.

STATE SPACE FORMULATION FOR THE SYS-

TEM DYNAMICS The governing differential equation for a typical linear time varying (LTV) mechanical system with (n×n) order mass, stiffness and damping matrices being M, K and C (n is number of dof ) can be represented as,

M(t)q(t) + C(t) q(t) + K(t)q(t) = f (t) -M(t)τ äg (t) (1)
q(t), q(t) and q(t) are the system's acceleration, velocity and displacement responses, respectively. f (t) is the external ambient force acting on the structure. äg (t) is the ground acceleration (as disturbance input) and τ is the location matrix that maps the effect of ground excitation to each of the nodes of the system.

With

F(t) = 0 n I n -M(t) -1 K(t) -M(t) -1 C(t) 2n×2n , B c (t) = 0 n M(t) -1 2n×m , E c = 0 τ 2n×l and x(t) = q T (t) qT (t)
T 2n×1 , the system dynamics, given in Equation (1), can be described in its continuous time as

ẋ(t) = F(t)x(t) + B c (t)u(t) + E c äg (t) + v(t) (2)
where x(t) denotes the unobserved system state. 0 n and I n are the n th order null and identity matrices respectively. u(t) m×1 , being the ambient forcing on the structure, can be modeled as a stationary white Gaussian noise (SWGN) of constant covariance Q m×m . v(t) is the process noise modeled as SWGN with covariance Q v 2n×2n . However, no SWGN model can be defined for äg l×1 . l is the number of channels for the disturbance input. The discrete time representation of the system is as follows,

x k = F k x k-1 + B k u k + E k äg k + v k (3) with F k , B k , E k , x k , u k , x k , v k
and äg k correspond to their continuous time entities described earlier obtained through zero-order-hold technique. The system's unobserved states are estimated using dynamic strain as measurement with the measurement equation defined as,

ε k = H k x k + w k (4)
where, H kp×2n is the time dependent measurement function that maps the system states x k to the corresponding strain measurements ε kp×1 . In reality, H k signifies the strain displacement relationship modeled using the finite element approach for all the members combined. This mapping takes basis on Equation (26) (Appendix) and combines for all measured members. w kp×1 denotes the measurement noise originating from the sensors and can be modeled as an SWGN of constant covariance R p×p .

3 SYSTEM ESTIMATION WITH ROBUST INTER-ACTING PARTICLE-KALMAN FILTER This article proposes an interacting strategy (also called as Rao-Blackwellisation [10-12, 8, 13, 7]) combining a Particle Filter (PF) with a Kalman Filter (KF), similar to IPKF algorithm given by [START_REF] Zghal | Modal parameter estimation using interacting kalman filter[END_REF], which further improvises the state filter with induced robustness through output injection following [START_REF] Zhang | State estimation for stochastic time varying systems with disturbance rejection[END_REF], discussed later. Overall, the approach employs a PF in which particles represent the posterior distribution of the health parameters as a numerical approximation, while KF deals with the uncertainty in the state space. PF being robust to nonlinear systems have been employed for parameter estimation while KF estimates the system states since a linear system is considered in this study [14]. Prior to demonstrate the proposed IPKF algorithm, the involved robustness aspect has been discussed in the following.

Robust Kalman filter

Typical KF requires either the exactly known inputs or an idealized SWGN model of it. For seismic excitation, none of these idealizations is valid which calls for induction of robustness within the KF algorithm [15][16][17][18]. A previously developed method [START_REF] Sen | Seismic-induced damage detection through parallel force and parameter estimation using an improved interacting particle-kalman filter[END_REF] achieves the robustness against unknown input forces by estimating the forces in real-time through an additional force filter. Yet, this additional computation makes the algorithm slow and less prompt from the damage detection perspective. The study employs [START_REF] Zhang | State estimation for stochastic time varying systems with disturbance rejection[END_REF]'s robust approach that rejects the impact of the unknown input äg k in the state evolution by suitably injecting a part of the measured output in the state transition model. Owing to the measurement equation (cf. Equation (4)), the following holds true for an arbitrary bounded matrix

G k ∈ R 2n×m , 0 = G k (y k -H k x k -w k ) (5)
Adding Equation (3) with Equation ( 5) and further setting L k = I 2n -G k H k , the state equation, Equation (3), can be modified as

x k = F k x k-1 + B k u k + E k äg k + v k + G k (y k -H k x k -w k ) = Fk x k-1 + Bk u k + Ẽk äg k + G k y k + ṽk (6) with Fk = L k F k , Bk = L k B k , Ẽk = L k E k and ṽk = v k -G k w k . If G k is chosen such that G k = E k (H k E k ) † ,
with † denoting Moore-Penrose Pseudo-inverse operation, Ẽk renders to a null matrix. Equation ( 6) can then be transformed to Equation [START_REF] Sen | Seismic-induced damage detection through parallel force and parameter estimation using an improved interacting particle-kalman filter[END_REF], decoupled from the unknown input äg k as

x k = Fk x k-1 + Bk u k + G k y k + ṽk (7)
Again ṽk can be idealized as an SWGN process. By this new formalism, Equation (3) has been transformed to Equation [START_REF] Sen | Seismic-induced damage detection through parallel force and parameter estimation using an improved interacting particle-kalman filter[END_REF] in order to reject the unknown input äg k with an appropriate injection of the known output y k . Equation (4) can still be used as the associated measurement equation. With this approach, the state filter can be employed without any explicit or statistical knowledge of the unknown input forces due to ground acceleration äg k .

PF-based parameter estimation

PF attempts system uncertainty propagation through a cloud of N p independent parameter particles

Ξ k = [ξ 1 , ξ 2 , • • • , ξ N ]
and their evolution [19,20]. At any arbitrary time step k, the evolution of an arbitrary particle ξ i k-1 is basically a random perturbation around its current position,

ξ j k = ξ j k-1 + N(δξ k , σ ξ k ) (8)
where a Gaussian blurring is performed on ξ j k-1 with a shift δξ k and a spread of σ ξ k 1 . The perturbed positions for the particles are then re-centered for consistency. This re-centring is actually a second step evolution in which particles are pushed towards the particle mean. For this, the same evolution strategy is adopted, as has been adopted in [START_REF] Sen | Seismic-induced damage detection through parallel force and parameter estimation using an improved interacting particle-kalman filter[END_REF]. The evolved particles are subsequently weighted based on their likelihood against the current measurement to find the posterior.

Nested KF-based state estimation

Inside the KF, state is estimated for each of the j th particle. The estimation propagates current estimation (x j k-1|k-1 ), conditioned on the current parameter estimate for each j th particle, to the predicted value, x j k|k-1 following Equation ( 9),

x j k|k-1 = F j k x j k-1|k-1 + G j k y k (9)
with corresponding state covariance, P j k|k-1 ,

P j k|k-1 = F j k P j k-1|k-1 F jT k + Bj k Q k BjT k + G j k R k G jT k ( 10
)
with Q k and R k being covariances of process and measurement noise. The corresponding measurement predictions y i,j k|k-1 can be obtained exploiting Equation (4). Associated innovation, j k|k-1 , can be obtained through estimating the de-

1 A + BN(µ, σ) means A + Bz where z follows N(µ, σ)
parture of prediction from the real measurement y k . Next, Kalman gain, K j k , is calculated as,

K j k = P j k|k-1 H T k H k P j k|k-1 H T k + R k (11)
The predicted states and state covariance are updated accordingly,

x j k|k = x j k|k-1 + K j k j k|k-1 (12)
P j k|k = (I -K j k H k )P j k|k-1 (13)

Particle approximation

Likelihood of each particle, i.e. L(ξ j k ), is further calculated based on the innovation mean j k|k-1 and co-variance, S j k , of each KF. The normalized weight for each j th particle is further updated as,

w (ξ j k ) = w (ξ j k-1 )L(ξ j k ) N j=1 w (ξ j k-1 )L(ξ j k ) (14) with L(ξ j k ) = (2π) n |S k | -1
e -0.5¯ jT S k -1 ¯ j . The particle approximations for the parameters and states are then estimated as:

x k|k = N j=1 w (ξ j k )x j k|k ; and θ k|k = N j=1 w (ξ j k )ξ j k (15)

NUMERICAL EXPERIMENT

The numerical investigation is performed with an FE model of a two dimensional seven-storey one-bay shear frame (cf. Figure 1) considered as the "test structure" for the experimentation (as a proxy for the real concrete structure). The relevant material and geometric details are presented in Table 1. The model is defined with six noded (three translational and three rotational) Euler Bernoulli beam elements. A 2% Rayleigh damping is further assumed for damping. The structure is excited at its base using north-south component of El-centro earthquake (Data source: http: //peer.berkeley.edu/research/motions/). Synthetic strain responses are further simulated from this model under ambient as well as base excitation in its healthy as well as damaged condition. The sampling frequency, time series length, time of arrival of seismic excitation, and time of damage occurrence has been maintained constantly at 50 Hz, 1024, 2 seconds and 3 seconds, respectively. The sensors are assumed to be patched onto the soffit of beam members only. Dynamic strain measurements collected from the members have been considered as response data which is further contaminated with stationary white Gaussian noise (SWGN) of signal-to-noise-ratio (snr ) 1% to replicate real life scenario.

It is obvious that with increasing size of particle pool, the estimation accuracy and precision is supposed to be enhanced. Yet, this comes at an enhanced computational cost as well, which ideally should be maintained at a manageable limit. The numerical experiment employed 3000 particle for PF. Figure 2 presents the adopted north-south component of the seismic ground motion time history taken from the El-Centro ground motion data. In parallel, the comparison between estimated strains against measured strains is also demonstrated. The strain estimation is done through reproducing the member strains from estimated states. It can be verified from Figure 2, that even after introduction of the non-stationary base excitation, the proposed algorithm perfectly estimated the strains as well as detected damage location and intensity with accuracy and precision. The real time damage detection is presented in Figure 3, where from the promptness and accuracy of detection can be verified.

EXPERIMENTAL STUDY -CANTILEVER BEAM

The proposed algorithm has been experimentally tested on an aluminium cantilever beam with material and geometric specifications shown in Table 2.

The beam is segmented into 4 elements with a focus to estimate health indices (hi s) for each of the segments. The first three elements of the beam starting from the fixed support are instrumented with strain gauges (Gauge ). An additional accelerometer is also attached to the tip of the beam to obtain frequency response function (FRF) for model calibration purpose (cf. Figure 4). Both the accelerometer and the strain gauges are sampled at a fixed sampling rate of 2000 Hz for very long time (30 sec). The cantilever beam is subjected to arbitrary as well as known input forces through an impact hammer. The known input force later facilitates in obtaining the FRF for the beam set up, required for model calibration.

Next, a primary finite element model for the undamaged experimental beam set up is created using Euler Bernoulli beam element with assumed material properties, which is later calibrated using the FRF obtained from modal domain decomposition of the measured accelerations with Fast Fourier Transform (FFT). This first phase of calibration was done manually targeting to bring the primary model sufficiently close to the reality. For this, initial material elasticity is adopted from tensile testings previously performed on the same material. This elasticity value is then uniformly altered for all the elements to arrive close to the experimentally obtained FRF.

Simultaneously, FFT of the strain gauge signals has been performed to cross-check the working of the assembled fullbridge strain gauges. Signals obtained from accelerometer and strain gauges led to the same experimental frequency (first two frequencies: ω 1 = 27.34 Hz and ω 2 = 187.40 Hz) of the undamaged cantilever beam which also ascertained the reliability of the full-bridge assembly. The manually calibrated model is then employed as the undamaged predictor model with the proposed r-IPKF algorithm along with the strain gauge measurements obtained from the undamaged structure to further calibrate the element elasticity in accordance to the experimental conditions. A set of four hi s has been defined within a range of 0 to 100% that gets multiplied with the assumed elasticity values in the primary model. These hi s are then estimated using the proposed algorithm as parameter states. Thus the health of the element can also be tracked using such hi s assuming 0 and 100 % denoting complete failure and healthy state of the element, respectively. 500 particles have been selected for PF with α = 0.98 while the initial distribution for each hi is considered as Gaussian, N (1, 0.01). The estimated hi s are presented in Figure 6 where minute calibrations can be observed to be implemented on element level elasticity. Once calibrated, the updated predictor model perfectly replicates the reality.

Next, damage is induced in the cantilever by removing a volume of 7 × 22 × 3mm 3 from the e2 element of the as shown in the schematic diagram (cf. Figure 7). The signals (cf. Figure 8) obtained from the strain gauges patched on the elements, e1 -e3, are then run through the r-IPKF. Figure 10, histograms of the input strain signals for the algorithm, demonstrates the non-stationary nature of the input forcing. The strain signals have a sampling rate of 100 Hz with a time series length of 2048. The result, thus obtained, shows a drop in pre-calibrated hi to hi = 83.5% for the second element, which corresponds roughly to a 16.5% damage in that element (cf. Figure 9) when compared to the initial rigidity value. The proposed method couples a particle filter for health parameter estimation with Kalman filter that estimates the system states from the measured noisy response data. The robustness against arbitrary input forcing in induced within the state filter following an output injection technique given by [START_REF] Zhang | State estimation for stochastic time varying systems with disturbance rejection[END_REF]. Strain has been adopted in this study for the measurement due to their cheaper installation cost and frame invariant nature. The proposed method has been validated using simulated strain response data from a numerical frame structure as well as on a real cantilever beam structure. While for the numerical experiment, nonstationary seismic base excitation has been employed as force, for the experiment, it has been ensured that the input forcing is not stationary or Gaussian. Both validation study establishes the promptness, robustness and efficacy of the proposed algorithm. 7 APPENDIX: STRAIN-DISPLACEMENT RELA-TION With strain as measurement, a measurement model is required to map the states to corresponding strains. The system dynamics defined in equation 1 has been solved using FEM approach and the nodal displacement values are obtained. With a strain-displacement relationship defined using the associated shape functions, strains at members can be obtained through interpolation. In the following, the pertinent strain-displacement relation for a frame structure is detailed. 3D Euler-Bernoulli beam has been considered as the elements for the FEM solution of the frame structure. Relevant schematic is presented in Figure 11. For each element, two transverse displacements, w y (x, t) and w z (x, t), are the field variables which can be interpolated at any arbitrary point x along the length of the beam (cf. Figure 11) using the nodal displacements and the associated shape functions. Further, deflection of a frame member at any direction can be idealized as summations of deflection due to its beam action and truss action, as: w y (x, t) = w b y (x, t) + w t x (x, t) w z (x, t) = w b z (x, t) + w t x (x, t)

The deflection transverse to the beam orientation due to beam, i.e., w b (x, t), action at any arbitrary distance x can be interpolated using the nodal displacements q l (t) defined in local coordinate system as:

w b (x, t) = w b y (x, t) w b z (x, t) = ψ b (x)q l (t) (17)
where, ψ b (x) is the interpolation matrix that transforms nodal displacements q l (t) to corresponding transverse displacements. q l (t) can further be retrieved from the nodal displacements defined in global degrees of freedom through coordinate transformation as:

q l (t) = Tq(t) ( 18
)
where T is the coordinate transformation matrix. ψ b (x) is constituted with the interpolation functions as: with η = 2x/L -1 and, nodal displacement q l (t) can be expanded as:

ψ b (x) = 0 ψ 1 (x) 0 0 0 ψ 2 (x) • • • 0 0 ψ 1 (x) 0 ψ 2 (x) 0 • • • • • • 0 ψ 3 (x) 0 0 0 ψ 4 (x) • • • 0 0 ψ 3 (x) 0 ψ 4 (x) 0 (19)
q l (t) = q e x q e y q e z q e θx q e θy q e θz • • •

• • • q f x q f y q f z q f θx q f θy q f θz T (21)
Similarly, slope and curvature at any arbitrary point x can be calculated as:

where,

w b (x, t) = ψ b (x) q l (t) ( 23
)
The curvature is further measured in the form of longitudinal strain (ε xx ) using a strain gauge patched at the surface of the beam. This gauge will register total strain due to both beam action and truss action. The strain due to bending in the beam is obtained from this curvature equation as:

ε xx,b (x, t) = -y 0 /κ b (x, t) (24)
where y 0 is the distance of extreme fiber from the neutral axis where the strain gauge is patched. The strain due to truss action can also be obtained as:

ε xx,t (x, t) = - 1 l (q e x -q f x ) ( 25
)
where l is the length of the element. The total strain can further be obtained as:

ε xx (x, t) = ε xx,b (x, t) + ε xx,t (x, t) (26)

Table 1 :Figure 1 :

 11 Figure 1: Schematic diagram of the experimental frame structure

Figure 2 :Figure 3 :

 23 Figure 2: Reproduced strain from the estimated states

Figure 4 :

 4 Figure 4: Experimental set up

Figure 6 :

 6 Figure 6: Model calibration with undamaged response

Figure 7 :

 7 Figure 7: Schematic model of the damaged cantilever beam

  Measured strain response at third sensor point

Figure 8 :Figure 9 :

 89 Figure 8: Measured strain response of the damaged cantilever beam

( a )

 a Histogram of measured strain response at first sensor point (b) Histogram of measured strain response at second sensor point (c) Histogram of measured strain response at third sensor point

Figure 10 :

 10 Figure 10: Histogram of measured strain response of the damaged cantilever beam

Figure 11 :

 11 Figure 11: Assumed Euler-Bernoulli Beam element with local and global degrees of freedom

  η) 2 (2 -η); ψ 4 (x) = L 8 (1 + η) 2 (η -1)(20)

) 3 2 w

 2 θ b (x) = w b (x, t) = ψ b (x)q l (t) κ b (x) = (1 + θ b (x, t) 2 b (x, t)


Neha Aswal
email: nehaaswal96@gmail.com

Subhamoy Sen
email: subhamoy@iitmandi.ac.in

Laurent Mevel
email: laurent.mevel@inria.fr

Damage Detection in Tensegrity using Interacting Particle-Ensemble Kalman Filter

Keywords: Tensegrity, Structural Health Monitoring, Interacting Particle Ensemble Kalman Filter, Damage Detection

Tensegrity structures form a special class of truss with dedicated cables and bars, that take tension and compression, respectively. To ensure equilibrium, the tensegrity members are required to be prestressed. Over prolonged usage, the cables may lose their prestress while bars may buckle, affecting the structural stiffness as well as its dynamic properties. The stiffness of tensegrities also vary with the load even in the absence of damage. This can potentially mask the effect of damage leading to a false impression of tensegrity health. This poses a major challenge in tensegrity health monitoring especially when the load is stochastic and unknown.

Present study develops a vibration based output-only time-domain approach for monitoring the health of any tensegrity in the presence of uncertainties due to ambient force and measurement noise. An Interacting Particle Ensemble Kalman Filter (IPEnKF) has been used that can efficiently monitor tensegrity health from contaminated response data. IPEnKF combines a bank of Ensemble Kalman Filters to estimate response states while running within a Particle Filter envelop that estimates a set of location based health parameters. Further to make damage detection cheaper, strain responses are used as measurements. The efficiency of the proposed methodology has been demonstrated using numerical experiments performed on a simplex tensegrity.

Introduction

Coined by Buckminster Fuller, the term tensegrity is a combination of the wordstension and integrity. As the name suggests, these structures derive their integrity from the prestress present in their members. From being introduced as a simple yet elegant piece of art by Snelson to being utilised in various complex and state of the art engineering structures like masts, stadium roofs, bridges, space and aerospace structures, tensegrity has come a long way in less than a century. Tensegrities are light-weight and help creating large column-free spaces. Many researchers have developed various design methodologies for complex tensegrities that are statically stable and can be constructed into reality. Tensegrities are self-stiffening structures. They can modify their stiffnesses through morphing their shape in order to carry any external load. Consequently, the vibrational response of the structure [START_REF] Ashwear | Natural frequencies describe the pre-stress in tensegrity structures[END_REF] gets affected by the external loading. It should therefore be noted that, stiffness alteration due to member prestress modification induced by force variability does not imply a damage in the structure.

Hence, to identify a possible structural anomaly in tensegrity, it is important to have an Structural Health Monitoring (SHM) approach that takes the inherent nature of tensegrity [START_REF] Sabouni-Zawadzka | Inherent properties of smart tensegrity structures[END_REF] into account, yet very less literature is available on this. Bhalla et. al. [START_REF] Bhalla | Damage assessment of tensegrity structures using piezo transducers[END_REF] successfully assessed the damage existing in tensegrity using dynamic strain measurements. Sychterz and Smith [START_REF] Sychterz | Using dynamic measurements to detect and locate ruptured cables on a tensegrity structure[END_REF] compared three methods to detect damage in tensegrity, namely, frequency analysis, error-domain model falsification (EDMF) using node position measurement and moving-window principal component analysis (MWPCA) using strain measurements. It was observed that the natural frequencies and the mode shapes cannot be used to detect and localize the damage in a tensegrity. The complexities involved with tensegrity mode shapes do not support damage identification. EDMF tends to become costly when tracking positions at submillimeter resolution and the results are sensitive to the amount of ambient uncertainty. MWPCA has an advantage as it uses inexpensive strain gauges and is efficient with noisy signal with high signal to noise ratio (SNR) but fails at low SNR levels.

Evidently, most of the works pertaining to tensegrity health monitoring is defined for deterministic scenarios. Nevertheless, for a real one, uncertainties arise during modelling error, ambient forcing and noise in the measurement which cannot be accounted for in a deterministic approach. In order to propose a practical SHM approach aimed to monitor real life tensegrities, it is necessary to incorporate the effects of uncertainties originating from unavoidable model inaccuracies, sensor noises and unknown external disturbances.

Bayesian filtering has proved its efficiency in detecting damage in structures in presence of uncertainties. Bayesian filtering approaches define the structural dynamics using a set of unobserved variables-termed as "states". The state dynamics and their impact on the measured responses are further defined using two separate probabilistic models, namely, process and measurement models/equations which individually account for uncertainties arising from ambient forcing and sensor noises. The inherent self-stiffening property of tensegrities [START_REF] Sabouni-Zawadzka | Inherent properties of smart tensegrity structures[END_REF] can be modelled by considering the geometric non-linearity in the dynamic model. This physical space dynamic model subjected to a stochastic input forcing is further transformed into a discrete time state space model to facilitate the estimation of the system states. A set of location based health parameters, θ k , to monitor the tensegrity health in terms of stiffness deterioration, is further added to the state definition as parameter states alongside the regular system states.

Evidently, the process equation is non-linear. In addition, the mapping of the mentioned health parameters to the corresponding measured responses is also non-linear. Handling these non-linearities simultaneously poses a major challenge in tensegrity SHM. Using an Interacting filtering technique can be a reliable option [START_REF] Ghanem | Health monitoring for strongly non-linear systems using the ensemble kalman filter[END_REF] for such system health estimation problems. This paper attempts to employ an Interacting Particle Ensemble Kalman Filtering (IPEnKF) strategy to estimate health of tensegrity in presence of process and measurement uncertainties. The major benefit with this approach is that the computational efficiency can easily be enhanced by exploiting the flexibility of parallelization of the complete computation, allowing to target a prompt detection of damage.

Geometric Non-linear Finite Element Model

Non-linearity in a structure may arise due to geometric (large deformations, etc.), material (plasticity, etc.) or boundary (contact) non-linearity. Geometric non-linearity is only considered in this study since strain-displacement relation, in a tensegrity, is prominently non-linear. Since tensegrity can typically be classified as a special case of truss, its finite element model (FEM) can be considered with a typical bar element defined with natural discretization [START_REF] Kebiche | Geometrical non-linear analysis of tensegrity systems[END_REF].

Fig. 1: LCS and GCS for bar element type

The transformation from global coordinate system (GCS) to local coordinate system (LCS) is performed using the transformation matrix, T as,

q l k = Tq g k (1) with, T = cosθ x cosθ y cosθ z 0 0 0 0 0 0 cosθ x cosθ y cosθ z , q l k = {q l 1 q l 2 } T k and q g k = {q g 1 q g 2 q g 3 q g 4 q g 5 q g 6 } T k .
Subscript k denotes the time step at which the stiffness is being defined. The deformation at any point in the element can further be defined using its nodal displacements and shape functions N 1 and N 2 as,

u * k = [N 1 N 2 ]q l k ( 2
)
where N 1 = (1 -r)/2 and N 2 = (1 + r)/2 are the shape functions defined in the natural coordinate system with r being the natural variable defined within the range -1 ≤ r ≤ 1. The strain field can further be obtained differentiating the displacement field including second order expansion as,

ε (e) k = ∂u * k ∂x + 1 2 ∂u * k ∂x 2 (3)
which leads to the following expansion including the linear and non-linear parts of the strain as, ε

(e) k = {B L q l k } linear + {B N L (q l k )q l k } non-linear (4)
where

B L = ∂N1 ∂x ∂N2 ∂x and B N L (q l k ) = 1 2 q l k T ∂N 1 /∂x ∂N 2 /∂x ∂N1 ∂x ∂N2 ∂x
are linear and non-linear strain-displacement matrices, respectively. Employing the principle of virtual work, the element stiffness matrix can be obtained by equating the work done by the internal forces to the virtual work done by the external forces.

V (e) ε (e) k σ (e) k dV = {δq l k } T F (e) k (5) with σ (e)
k being the member stress that can be obtained from the constitutive relation as σ

(e) k = E (e) ε (e)
k with E (e) being the constitutive matrix. With this, Eq. ( 5) can be expanded as,

V (e) δq g k T T T B T k E (e) B k Tδq g k dV = {δq g k } T F (e) k (6)
with

B k = B L + B N L (q l k ).
Imposing the non-trivial solution condition, {∂q l k } = 0 and further comparing with

F (e) k = K (e)
k δq g k , the element stiffness matrix in the natural coordinate system can be obtained as,

K (e) k = E (e) A (e) 1 -1 T T B T k B k T l (e) 2 dr (7)
assuming a uniform cross section A (e) over the entire length l (e) of the element e. The Gauss-Quadrature method has further been applied to find the numerical integration for the above integral (ng = 1). Assembling all the element stiffness matrices and applying the natural boundary conditions, one can get the global stiffness matrix K k . Similarly, the mass matrix M can be obtained by following the consistent mass matrix assumption.

Non-linear Response from Newmark-beta

Due to its superior accuracy, Newmark-beta is one of the popular methods to find the response of the non-linear systems. The method solves an incremental equilibrium equation (see Eq. ( 8)) to find non-linear structural response variables , i.e., ük , uk and u k , of the following governing differential equation,

M∆ü k + C k (u k )∆ uk + K k (u k )∆u k = ∆P k (8)
where, ük (n×1) , uk (n×1) and u k (n×1) are the acceleration, velocity and displacement responses, respectively and the ∆ operator signifies the corresponding increment over each time step. M (n×n) is the time invariant global mass and C k (n×n) and K k (n×n) are the time dependent damping and stiffness matrices of the structure, respectively. The time dependence of K k originates from its dependence on the displacement response u k . With Rayleigh damping assumption, C k being a function of K k is also time dependent. P k (n×1) is the external force acting on the structure. Clearly, the above dynamics is an implicit equation which requires an iterative solution at each step. However, in the context of tensegrity SHM, the same equation can be redefined using an explicit formulation with Newmark-beta explicit integration scheme without losing much on the precision yet achieving promptness in detection. Newmark-beta is a well established method and for the sake of brevity, it is not discussed here. Details of this method can be found in [START_REF] Chopra | Dynamics of structures, a primer[END_REF]. The method shows an unconditional stability for γ = 0.5 and β = 0.25 (average constant acceleration assumption), hence these values are utilised in this study.

State space formulation for system dynamics

Bayesian filtering employs a process equation to define the system dynamics in terms of dynamics of a set of unobserved variable, named as states. To estimate the time varying nature of the system due to damage, a set of location based health parameters, θ k , are additionally defined as parameter states. These health parameters monitor the residual elasticity after damage within a range of [0, 1] in which 0 and 1 mean complete and no damage, respectively. While the system states are estimated using a bank of EnKF (Ensemble Kalman Filter), the parameter states are estimated through an envelop PF (Particle Filter).

Within EnKF, the system dynamics is defined with displacement, velocity and acceleration responses in the state vector, i.e. X k = {u k uk ük } T . The previously detailed non-linear FEM of tensegrity is further employed to propagate these system states in time. Prior to this, the tensegrity stiffness needs to be re-calibrated considering the current nodal displacement that causes a change in the prestress. Integrated using the Guass-Quadrature method, Eq. ( 7) is discretized to find the discrete elemental stiffness matrix with B k-1 re-calibrated as B k using the current nodal displacements u k-1 . All the elemental stiffness matrices are further assembled to obtained the global stiffness matrix K k . Next, the system states are forwarded in time using the Newmark-beta integration scheme.

X k = f (X k-1 , K k , M, dt, P k ) + v k ( 9
)
where f denotes the intrinsic calculations involved in Newark-beta integration scheme. dt denotes the time step for which the simulation is performed. P k is the stochastic ambient force acting on the structure modelled as a white Gaussian noise (WGN) of distribution N(0, Q P k ). 1 The process noise v k accounts for the unavoidable inaccuracies in the prediction model f (•) and can be modelled as WGN of distribution

N(0, Q k ).
Next, the unobserved response states are mapped to the corresponding strain fields in several members using the non-linear strain measurement mapping equation. The dependence of the strain-displacement matrix B k on the displacement u l k can be demonstrated using the notation B(X k ) since u k is a subset of X k . The measurement equation can then be obtained from Eq. ( 4) by incorporating the partial measurement information and the sensor noises as,

y k = HB(X k )TX k + w k ( 10
)
where H is the selection matrix that defines which of the member strains are measured. w k accounts for measurement noises arising from the sensors modeled using WGN model of distribution N(0, R k ).

System estimation using IPEnKF algorithm

The system states are estimated using an IPEnKF algorithm which is a modification of the Interacting Particle-kalman Filter (IPKF) algorithm [START_REF] Sen | Seismic-induced damage detection through parallel force and parameter estimation using an improved interacting particle-kalman filter[END_REF]. With IPEnKF, the response states are estimated using a nested EnKF instead of KF to enable handling the non-linear systems, while an enveloping PF estimates the parameter θ k . All of the filtering steps are briefly described in the following. The details can be found in the article [START_REF] Sen | Seismic-induced damage detection through parallel force and parameter estimation using an improved interacting particle-kalman filter[END_REF].

Envelop parameter filter -PF

In order to avoid an explicit analytical integration over the entire parameter space, θ k , a particle approximation of this integration is approached with PF that propagates the parametric uncertainty thorough a cloud of N p independent particles

Ξ k = [ξ 1 k , ξ 2 k , • • • , ξ Np k ].
At each time step, the particle distribution gets updated depending on their likelihood against measured data and thus no assumption is needed on their distribution apriori.

Within the envelop PF, the particle representations of parameter estimates evolve assuming a Gaussian perturbation around the current estimate ξ j k-1 as,

ξ j k = ξ j k-1 + N(δξ k ; σ ξ k ) (11)
where a Gaussian blurring is performed on ξ j k-1 with a shift δξ k and a spread of σ ξ k . For each of the particles, the corresponding stiffness matrix, i.e., K j k-1|k-1 , is first calculated and supplied to each of the EnKF candidates.

Nested state filter -EnKF

In the nested EnKF, N e sets of state ensembles for each j th parameter particle are propagated through the system state propagation model (see Eq. ( 9)). However, since the tensegrities are functions of nodal displacements, the supplied stiffness matrix K j k-1|k-1 has to be further modified based on the current estimate of nodal displacement responses q g k-1|k-1 i,j which is a subset of the state estimate X i,j k-1|k-1 pertinent to j th particle and i th ensemble. Eq. ( 7) can be employed and can be expressed as,

K i,j k|k-1 = m(K j k-1|k-1 , X i,j k-1|k-1 ) (12)
where m(•) is the stiffness calibration function detailed in Eq. [START_REF] Chopra | Dynamics of structures, a primer[END_REF]. With the modified stiffness matrix K i,j k|k-1 and a realization for the ambient forcing, P i,j k , drawn from a known statistic (say N (0, Q f )), the prior state estimates X i,j k-1|k-1 are propagated to the next time step as X i,j k|k-1 ,

X i,j k|k-1 = f (X i,j k-1|k-1 , K i,j k|k-1 , M, dt, P i,j k ) + v i,j k (13)
The measurement function described in Eq. ( 10) is further employed to map each of the ensembles to the corresponding measurement prediction, y i,j k|k-1 from which innovation, i,j k|k-1 , for each ensemble, can be obtained from the corresponding real measurement y k . Next, the ensemble mean estimates for the propagated states (X j k|k-1 )

and predicted measurements (y j k|k-1 ) are obtained as means of corresponding entities for all ensembles. Cross-covariance between state and measurement prediction C j,XY k and the measurement prediction covariance C j,Y Y k can further be computed following [START_REF] Evensen | The ensemble kalman filter: Theoretical formulation and practical implementation[END_REF]. The innovation error covariance, S j k and EnKF gain, G j k , are then obtained as

S j k = C j,YY k + R k and G j k = C j,XY k (S j k ) -1
. With this gain, the state ensembles are updated as,

X i,j k|k = X i,j k|k-1 + G j k i,j k|k-1 (14)

Particle approximation

Likelihood of each particle, i.e. L(ξ j k ), is further calculated based on the innovation mean j k|k-1 and co-variance, S j k , of all the ensemble simulations running within it. The normalized weight for each j th particle is further updated as,

w (ξ j k ) = w (ξ j k-1 )L(ξ j k ) N j=1 w (ξ j k-1 )L(ξ j k ) (15) with L(ξ j k ) = (2π) n |S k | -1
e -0.5¯ jT S k -1 ¯ j . The particle approximations for the parameters and states are then estimated as:

x k|k = N j=1 w (ξ j k )x j k|k ; and θ k|k = N j=1 w (ξ j k )ξ j k (16)

Numerical Experiment

The proposed SHM approach has been tested on a simplex tensegrity. Simplex tensegrities (see Fig. 2) are a type of cylindrical tensegrity with 3 bars and 9 cables. For dynamic analysis, it is fixed at its bottom nodes (1-3) while ambient Gaussian force is applied in x-dir on the fourth node (see Fig. 2). 90% damage is induced in element 11 (bar) of the simplex. Nine strain gauges are attached on all the unrestricted members, i.e, elements 4-12, from where the measurements are collected. The initial statically stable coordinates (see Fig. 2) are obtained from [START_REF] Attig | Output-only modal identification of tensegrity structures[END_REF]. The stable form of the simplex is further excited with an ambient WGN forcing of variance 10 4 N throughout the simulation. Damage is initiated 0.5 s after the start of the simulation. The response from strain gauges is sampled at a fixed sampling frequency of 100 Hz for 5 seconds and further contaminated with a WGN of 1% SNR.

The initial distribution type for the parameter particles (health indices) is set to be Gaussian, with their mean set as 1 assuming an undamaged condition and a Fig. 2: Simplex Tensegrity Configuration standard deviation of 0.02. 2500 filter particles are selected for PF while 75 ensembles are chosen for EnKF. A better precision may however be expected with a higher number of particles or ensemble in exchange of a higher computational cost. Fig. 3 compares the effect of various SNR (1%,2% and 5%) on the proposed approach for simplex. It can be observed that the damage percentage estimated is fairly precise and prompt for all the SNR levels, but the promptness of detection is rather poor for 5% SNR. This, however, can be improved at higher computational cost. The method is observed to perform poorly for 10% SNR levels which demarcates the limitation of the approach.

Conclusion

A novel interacting filtering based damage detection technique has been proposed for tensegrity structures. EnKF is coupled with PF for simultaneous estimate of system states along with identification of a set of location based health parameters that monitor the stiffness deterioration. The health assessment is performed using ambient vibration response without knowing the input forces explicitly. The proposed method is found to be prompt and precise in detecting damages in a simplex tensegrity. The sensitivity of the proposed method under different SNR levels (1%.2%,5% and 10%) has been investigated and limiting SNR level (10%) is also identified.

Funding: This study was funded by Science & Engineering Research Board (SERB), New Delhi, India, through grant file no. ECR/2018/001464.

Fig. 3 :

 3 Fig. 3: Measurement noise sensitivity of proposed approach for simplex tensegrity

A + BN(µ; σ) means A + Bz where z follows N(µ; σ)