Julien Bonnel
email: julien.bonnel@gipsa-lab.inpg.fr

Gregoire Le Touzé

Barbara Nicolas

Jérôme I Mars

Cedric Gervaise

Automatic and passive whale localization in shallow water using gunshots

This paper presents an automatic and passive localization algorithm for low frequency impulsive sources in shallow water. This algorithm is based on the normal mode theory which characterizes propagation in this configuration. It uses specific signal processing tools and time-frequency representations to automatically extract features of the propagation. Then, it uses the dispersive properties of the oceanic waveguide as an advantage to perform the localization. Only few hydrophones are needed and neither knowledge of the oceanic environment nor simulation of the propagation is required. The proposed method is successfully applied on North Atlantic Whale gunshots in the Bay of Fundy recorded with a network of three hydrophones.

I. INTRODUCTION

The study of marine mammals is a difficult task as most of the visual observations are closely restricted by weather, daytime and environment. However, passive acoustics which only consists in "listening" the acoustical environment could provide another tool to obtain spatial and temporal distribution of marine mammals [START_REF] Spiesberger | Passive localization of calling animals and sensing of their acoustic environment using acoustic tomography[END_REF]. This could be helpful to understand and protect these animals, especially when endangered species are concerned. Indeed with passive acoustics, animals are located thanks to their own calls, in opposition to classical SONAR methods where a signal has to be emitted. This allows a more autonomous system, requiring less energy, and save the sea fauna being disturbed by active acoustic signals. Moreover, it can be used when visual observations failed.

When a marine mammal emits a call, the sound travels from the animal position to one or several hydrophones. The received signal depends on the oceanic environment, and on the positions of both source and receiver. Using signal processing techniques, it is possible to extract features from this signal and use them to estimate the source localization [START_REF] Laurinolli | Localization of North Atlantic Right whale sounds in the bay of Fundy using a sunobuoy array[END_REF] [START_REF] Wiggins | Waveguide propagation allows range estimates for north pacific right whales in the Bering sea[END_REF]. Localization algorithms are usually based on an acoustic propagation model and require the knowledge of oceanic environment.

Right whales calls have been described in [START_REF] Vanderlaan | Taggart Characterization of North Atlantic right whale (Eubalaena glacialis) sounds in the Bay of Fundy[END_REF]. They are often low frequency calls, but have various waveforms: constant low-frequency, moan, up and down sweeping modulations, and gunshot. Gunshots are loud impulsive sounds from 10Hz to 20kHz lasting approximately 2ms. They are produced by lone males (or small groups), mainly in the bay of Fundy (Canada) and are probably implied in reproduction [START_REF] Parks | The gunshot sounds produced by male north atlantic right whales (Eubalaena glacialis) and its potential function in reproductive advertisement[END_REF]. As they are emitted near the surface, they could be used for an automatic alert system to avoid whales and ship collisions. Moreover, the bay of Fundy is a shallow water area with internal tides producing large and quick variation of sound speed profile [START_REF] Desharnais | A description of the workshop dataset[END_REF]. This implies that a robust localization algorithm has to be developped. This paper presents a passive localization method (in a 2D horizontal plane) for low-frequency transient signals (such as gunshots) in shallow water environment, using a sparse network of hydrophones. The first part of the paper will introduce the experimental data used for this study. The main ideas of the modal propagation model on which relies our scheme are presented in a second part. A third part will describe the algorithm itself, including the necessary signal processing and time-frequency tools. Finally, the method will be applied on real data and discussion will be done.

II. THE DATA SET

The dataset used in this paper comes from the 2003 Workshop on detection, localization and classification of marine mammals using passive acoustics. The acquisition system is composed by five Ocean Bottom moored Hydrophones. Their localization is given in Table I. They have a flat sensibility from 50Hz to 700Hz and the data were digitized using a 12-bits A/D converter with a sampling frequency of 1200Hz. The area around the OBH is shallow water, with bathymetry from 100 to 200 meters. Sound speed profiles were measured during this experiment: they were downward refractive or had a local minimum. As said before, they also presented quick temporal variations. The bottom structure is mainly composed of a first Lahave clay layer over a thick layer of Scotian drift [START_REF]Surficial geology -Eastern Gulf of Maine and Bay of Fundy, Map 4011-G Canadian HydrographicService[END_REF]. The Lahave clay layer is characteristic as its compression sound speed is lower than the sound speed in water. It implies a high dispersion for normal mode propagation. The dataset contains several right whale sounds recorded in the bay of Fundy between 2000 and 2002 [START_REF] Desharnais | A description of the workshop dataset[END_REF], including gunshots. The Fig. 1 presents an example of a recorded gunshot (S035-2 on hydrophone H) in the time and the time-frequency domails. The latter presents a multicomponent pattern, which is typical of a dispersive normal mode propagation: each component has its own time of arrival which depends on frequency. We give some details on this propagation in the following section. In our configuration (shallow water and low frequency signals), the most suitable propagation model is normal mode theory. In this case, in a range independent environment, for a frequency f , the transfer function H between a source at depth z s and a receiver at a depth z r separeted by a radial distance r is [START_REF] Jensen | Computational ocean acoustics[END_REF] :

H(f ) ≈ Q ∞ m=1 g m (z s )g m (z r ) e jkr(m,f )r k r (m, f )r (1)
where g m is the m th modal function, k r (m, f ) the radial wavenumber of mode m (which is supposed to be real as the evanescent modes are not taken in account), and Q = e jπ/4 ρ(zs) (with ρ(z s ) the water density at the source depth). Thus, the propagation is multicomponent. For each component of index m, phase speed v φ and group speed v g can be defined by :

v φ (m, f ) = 2πf k r (2) v g (m, f ) = 2π ∂f ∂k r (3)
The group velocity describes the propagation speed of energy. We can note that v φ and v g depend both on frequency f and mode index m. Consequently, each frequency of each mode will travel with its own speed, which is the definition of a dispersive propagation.

If a transient signal is emitted with a time frequency modulation t e (f ) (t e is the time of emission of the frequency f ), the time-frequency structure T F R of the receive signal after modal propagation is:

T F R(t, f ) = ∞ m=1 A(m, f, r, z s , z r )δ(t -t e (f ) - r v g (m, f ) ) (4
) where δ(t) is the dirac distribution describing the localization of the time-frequencies structures and A is the attenuation term describing their amplitude.

IV. THE LOCALIZATION ALGORITHM

The main idea is to take advantage of the dispersive behaviour of the oceanic waveguide to localize a transient emission. For a source s emitting a transient signal with an unknown time-frequency modulation t e (f ), the arrival time of the frequency f of the mode m measured after propagation at a hydrophone n is given by:

t r (m, n, f ) = t e (f ) + r(s, n) v g (m, f ) (5)

A. Estimation of the arrival times

The first step of the algorithm is to estimate all the t r . It is impossible in the time domain as the modes are overlapped. It is neither direct in the time-frequency domain because the modes are broadband and close from each others. Consequently, each mode has to be first filtered.

1) Warping operators:

As proposed in [START_REF] Le Touzé | Matched time-frequency representations and warping operator for modal filtering[END_REF], the pressure signal will be warped in order to have a better representation of the modal information. The warping is based on a model of the environment and is computed with an unitary equivalence approach [START_REF] Baraniuk | Unitary equivalence: a new twist on signal processing[END_REF]. Here, the used model is the isovelocity one. It is made of a homogeneous layer of fluid between perfectly reflecting boundaries. Of course, this modelisation is simplistic, and does not match to the real oceanic environment. However, it is useful as it does not require information of the environment and is enough efficient for our goal.

In the isovelocity case, the pressure signal is given by [START_REF] Jensen | Computational ocean acoustics[END_REF]:

p(t) = m g m (t)e j2πνc(m)ξ(t) (6)
with g m (t) describing the envelop of the m th mode, ν c (m) the cutoff frequency of the m th mode (depending only on m, on the constant velocity V of the water and of the depth D of the waveguide) and ξ(t) the general dispersivity function which is:

ξ(t) = t 2 - r 2 V 2 (7)
The warping operator is based on a deformation function ω(t) and is noted W ω . Its aim is to linearize the pressure signal. Consequently, W ω p(t) must be a sum of linear structures and follows the equations :

W ω p(t) = m |ω ′ (t)|C m e j2πνc(m)ξ[ω(t)] (8)
W ω p(t) = m |ω ′ (t)|C m e j2πνc(m)t (9)
The corresponding deformation function ω is [START_REF] Le Touzé | Matched time-frequency representations and warping operator for modal filtering[END_REF]:

ω(t) = ξ -1 (t) = t 2 + r 2 V 2 (10)
If this operator is applied on a pressure signal (from an isovelocity waveguide), the modes of the warped signal will be sinusoids. There is also an inverse operator W -1 ω . It is defined by

W -1 ω W ω p(t) = p(t)
and is linked to the deformation function w -1 (t).

2) Modal filtering: This operator could be applied on a real pressure signal. As the real waveguide is not an isovelocity one and as r is unknown, the warped modes are not perfect sinusoids. However, they are quite well separated in the timefrequency domain and can be easily filtered. Here, a simple threshold is sufficient to create masks on time-frequency domain and filtered them. Consequently, the filtering scheme is as follow :

1) Using the isovelocity model, the recorded signal is warped: modes become nearly sinusoids (in the timefrequency domain, they are nearly horizontal lines). 2) A time-frequency representation of the warped signal is computed: Short Time Fourier Transform (STFT) is chosen as it is easily computable and allows filtering in the time-frequency domain using masks. 3) Masks are created to filter warped modes (in the timefrequency domain). 4) Each filtered (warped) mode expressed in the time domain with an inverse STFT. 5) Each warped mode are unwarped (in the time domain).

3) Final estimation of arrival times:

As each mode is filtered, it is easy to obtain the times of arrival t r (m, n, f ) (given by equation 5). For a given hydrophone n, the reallocated spectrogram of the mode m is computed [START_REF] Auger | Improving the Readability of Time-Frequency and Time-Scale Representations by the Reassignment Method[END_REF]. Indeed, a reallocated spectrogram reassigned position of time-frequency points in order to obtain thin structures. But even after reallocation, the time-frequency representation of the mode is not perfectly thin. Consequently, the instantaneous frequency of the mode is obtained by computing the skeleton of the shape of the mode. This is done for each mode on each hydrophone and gives all the t r (to have a better understanding of the whole procedure, figures will illustrate it on section V-A). Once the times of arrival are known, they could be used to localize the source.

B. A ratio only depending on source-receiver distance

By combining different arrival times, it is now possible to define a ratio which depends only on source-receivers distance [START_REF] Gervaise | Robust 2D Localization of low-frequency calls in shallow waters using modal propagation modelling[END_REF]. First, let us considere the arrival times of a frequency f of two modes m and m ′ on a single hydrophone n. By substrating them, the influence of the time of emmission t e (f ) disappears:

d(m, m ′ , n, f ) = t r (f, m, n) -t r (f, m ′ , n) (11)
d(m, m ′ , n, f ) = r sn ( 1 v g (m, f ) - 1 v g (m ′ , f ) ) (12)
As the group velocity is the same for all hydrophones, it is possible to combine the quantities d for two different hydrophones n and n ′ to get rid of them :

R(m, m ′ , n, n ′ , f ) = d(m, m ′ , n, f ) d(m, m ′ , n ′ , f ) = r sn r sn ′ (13)
The quantity R is computable only with the received measures. It does not depend on the waveguide properties nor on the emitted signal. Consequently, it is a suitable quantity to robustly estimate the localization of the source. Thus, in a horizontal plane, it restrains the position of the source to lie on a circle (or a line if R = 1).

C. A cost function to localize the source

Instead of considering geometrical properties, the localization of the source is estimated by minimising a cost function J, as in [START_REF] Gervaise | Robust 2D Localization of low-frequency calls in shallow waters using modal propagation modelling[END_REF]. With the extracted t r , all the possible 5-uplets (m, m ′ , n, n ′ , f ) are formed (with m = m ′ and n = n ′ ). Then, the source localization Ŝ(x s , y s ) is estimated by solving the following estimation problem, which reprensents the sum of the quadratic distance to the constraint:

Ŝ(x s , y s ) = arg min M (x,y) J[M (x, y)] (14)
with

J[M (x, y)] = N i=1 [d 2 M (x,y),n -R 2 (m i , m ′ i , n i , n ′ i , f i ) d 2 M (x,y),n ′ ] 2 (15
) where M (x, y) is the point where the cost function J is computed, d M (x,y),n is the distance between the hydrophone n and the point M (x, y), and N is the number of 5-uplets (m, m ′ , n, n ′ , f ) previously selected.

A. Presentation of the results

The localization scheme described on the previous section was applied to the recording 'S035-2' from the data set of the 2003 Workshop on detection, localization and classification of marine mammals using passive acoustics. It contains a gunshot of right whale recorded in the Bay of Fundy. For this work, only hydrophones H, E, and L were used. Thus, the gunshot was not properly recorded on hydrophone J, and the environment changes a lot around hydrophone C. As an example, the signal recorded on hydrophone H is given in Fig. 1. Its warped version is presented in Fig. 2. For all the hydrophones, the warping parameters have been arbitrary chosen with the values r = 8500 m and V = 1500 m/s. Even if it does not correspond to reality, the modes are nearly horizontal and well separated. On Fig. 2 it is noticeable that time and frequency scales have changed since it is in the warped domain. Modes are easily filtered. The mask created to filter the second one could be seen in Fig. 3. Once it is done, each mode is unwarped. A time-frequency representation (reallocated spectrogram) of the second mode is given on Fig. 4. Then, their instantaneous frequency is computed. It allows to obtain all the t r (f ) for this mode. On Fig. 5, the red lines (which are the instantaneous frequencies of each modes) are constitued by all the extracted t r (f ). Then, all the possible 5-uplets R(m, m ′ , n, n ′ , f ) are computed. As told before, only hydrophones H, E and L were considered. Different tests were done using only the first two modes, the first three modes, or the first four modes. For each 4-uplet (m, m ′ , n, n ′ ), 10 different frequencies were chosen linearly spaced in the biggest mutual frequency band of m and m ′ . Criterium J was computed with a step of 5 meters. Table II summarizes the localization's results depending on the number of modes used, while table III allows to compare the proposed method with other techniques on the same data [START_REF] Desharnais | Gillard Right whale localisation using a downhill simplex inversion scheme[END_REF] [6] [START_REF] Laurinolli | Hay Localization of right whale sounds in the workshop bay of fundy dataset by spectrogram cross-correlation and hyperbolic fixing[END_REF]. Fig. 6 presents the J criterium and the corresponding estimated position when only the first two modes were considered. when only the first two modes are considered so the results can be compared with the one obtained by Gervaise et al. in [START_REF] Gervaise | Robust 2D Localization of low-frequency calls in shallow waters using modal propagation modelling[END_REF] where only this two modes were used. We will discuss this point in section V-B. In the other cases, our method give the same kind of results (localization and shape of J).

B. Discussion

Our gunshot localization is compatible with the solutions obtained in the litterature, and similar results were obtained on S070-3 and 2013-1 recordings. Even if the method was designed for the localization of right whale in the Bay of Fundy, it could be used whenever one wants to localize low frequency transient sounds in shallow water, for example :

• North Pacific right whales and Humpback whales in the Bering sear [START_REF] Mc Donald | Calls recorded from North Pacific right whales (Eubalaena Japonica) in the eastern Bering sea[END_REF] • eastern North Pacific blue whales [START_REF] Oleson | Behavioral context of call production by eastern North Pacific blue whales[END_REF] • blue whale in the Saint Lawrence [START_REF] Berchock | St Laurence blue whale vocalization revisited: characterization of calls detected from 1998 to 2001[END_REF] 1) Analysis of the results: By looking at the shape of the criterium J in Fig. 6, it is noticeable that the precision of the estimated position of the gunshot is lower for direction x than for y. It is normal as the localisation of the gunshot is nearly outside the hydrophone network for the direction x. The results presented in table II show that x GS decreases as the number of modes used in the algorithm increases. It may be explained because the first modes are more energetic and easier to track, so the precision is higher. This is confirmed by Fig. 7, where the shape of J is sharper for only two modes. 2) Comparisons with other techniques: When the method is compared with the litterature, four points should be noted:

• Our method is an improvment of the one presented by Gervaise and al. in [START_REF] Gervaise | Robust 2D Localization of low-frequency calls in shallow waters using modal propagation modelling[END_REF] as it is fully automatic. Moreover, the frequency band of interest is now higher as it depends on the couple mode/receiver instead of being fixed. This allows to consider more modes. Consequently, it is difficult to know which method is the more accurate. We assume that our method is better when the environment is unknown or varying, but Desharnais' and Laurinolli's methods are probably more accurate when the environment is constant and well known. Fig. 8 allows a graphical comparison of our results with the three other methods discussed in this subsection.

• Our method is based on modal propagation, but does not require to run a normal mode propagation code [START_REF] Wiggins | Waveguide propagation allows range estimates for north pacific right whales in the Bering sea[END_REF]. It is a major advantage in such an environment which cannot be modelised with a simple Pekeris waveguide as there is a poorly compact first layer in the bottom and a time variable celerity profile.

VI. CONCLUSION

This paper presents a passive localization scheme for low frequency transient sounds in shallow water and its application to right whale gunshot in the Bay of Fundy, Canada. This method is fully automatic. Moreover, it does not require information of the environment properties, nor need to run a propagation code (whereas it is based on normal mode theory). As it uses relative time of arrival, it is not sensitive to clock's drift. All this properties make it suitable to obtain a real-time localization system, provided that a real time implementation of the time-frequency tool is possible. As perspective, the author woud like to apply the proposed method on a larger dataset to test the proposed approach, and also on simulations to determine its precision and its variability to the different parameters.

TimeFig. 1 .

 1 Fig. 1. Time representation and time-frequency representation of the recorded gunshot on OBH #H

Fig. 2 .

 2 Fig. 2. Time representation and time-frequency representation of the warped version of the gunshot recorded on OBH #H

Fig. 3 .Time

 3 Fig. 3. Mask created to extract the second mode of the warped version of the gunshot recorded on OBH #H

Fig. 4 .Fig. 5 .

 45 Fig. 4. Time representation and reallocated spectrogram of the second mode (after filtering and unwarping) of the gunshot recorded on the OBH #H

Fig. 6 .

 6 Fig. 6. Global view of the networks of OBH and results of the localization scheme using only two modes

• 4 Fig. 7 .

 47 Fig. 7. Criterium J for y = y S considering the two or the four first modes

Fig. 8 .

 8 Fig.8. Estimated position and J criterium using only the first two modes, and estimated positions with the methods presented in[START_REF] Desharnais | Gillard Right whale localisation using a downhill simplex inversion scheme[END_REF],[START_REF] Gervaise | Robust 2D Localization of low-frequency calls in shallow waters using modal propagation modelling[END_REF] and[START_REF] Laurinolli | Hay Localization of right whale sounds in the workshop bay of fundy dataset by spectrogram cross-correlation and hyperbolic fixing[END_REF]

TABLE III S035

 III -2, GUNSHOT LOCALIZATION WITH THE FIRST TWO MODES AND

	COMPARISON WITH OTHER METHODS
	Method	x GS (m)	y GS (m)
	Bonnel et al.	7886	-1008
	Gervaise et al. [6]	9225	-1248
	Desharnais et al. [3]	8884	-848
	Laurinolli et al. [9]	8950	-970


Ajinkya Kulkarni
email: ajinkya.kulkarni@loria.fr

Vincent Colotte
email: vincent.colotte@loria.fr

Denis Jouvet
email: denis.jouvet@loria.fr

Improving transfer of expressivity for end-to-end multispeaker text-to-speech synthesis

Keywords: TTS, metric learning, expressivity, transfer learning

The main goal of this work is to generate expressive speech in different speaker's voices for which no expressive speech data is available. The presented approach conditions Tacotron 2 speech synthesis with latent representations extracted from text, speaker identity, and reference expressive Mel spectrogram. We propose to use multiclass N-pair loss in the end-to-end multispeaker expressive Text-To-Speech (TTS) for improving the transfer of expressivity to the target speaker's voice. We have jointly trained the end-to-end (E2E) TTS with multiclass N-pair loss to discriminate between various emotions. This augmentation of the loss function during training paves the way to enhance the latent space representation of emotions.

We have experimented with two different neural network architectures for expressivity in the encoder, namely global style token (GST) and variational autoencoder (VAE). We transferred the expressivity using the mean of latent representation extracted from the expressivity encoder for each emotion. The obtained results show that adding multiclass N-pair loss based deep metric learning in the training process improves expressivity in the desired speaker's voice.

I. INTRODUCTION

The term expressivity in speech usually refers to the characteristics of speech, such as emotions, speaking style, relationship of speech with gestures, and facial expression. Throughout this paper, we consider only the emotional characteristics of expressivity in speech. The current end-to-end TTS systems heavily relies on a large amount of speech corpus used for training the system [START_REF] Wang | Tacotron: A fully end-to-end text-to-speech synthesis model[END_REF]. Therefore, to build expressive speech synthesis for a new speaker, one has to create a speech corpus with various emotions. It is inconvenient to record an expressive speech corpus every time somebody wants to build an expressive speech synthesis system for a new speaker's voice. Furthermore, creating an expressive speech corpus is laborious and expensive in terms of workload. Though, many approaches proposed to use audio-books, films, dialogues. to create expressive speech synthesis. However, labeling the expressions is not a trivial task due to a large number of possible variations in a single emotion [START_REF] Charfuelan | Expressive speech synthesis in MARY TTS using audiobook data and emotion[END_REF]. This creates a bottleneck in the development of expressive speech synthesis in a new speaker's voice.

There are numerous frameworks that have been proposed for the implementation of expressivity transfer either by inter-polation of latent representations of prosody or of prosody embedding [START_REF] Hsu | Hierarchical generative modeling for controllable speech synthesis[END_REF]- [START_REF] Lee | Robust and fine-grained prosody control of endto-end speech synthesis[END_REF]. For controlling expressivity, these approaches enhance the Tacotron based TTS system by the addition of advanced deep neural network architecture such as variational autoencoder (VAE) [START_REF] Zhang | Learning latent representations for style control and transfer in end-to-end speech synthesis[END_REF], Gaussian mixture VAE [START_REF] Hsu | Hierarchical generative modeling for controllable speech synthesis[END_REF], Global style token [START_REF] Wang | Style tokens: Unsupervised style modeling, control and transfer in end-to-end speech synthesis[END_REF], FLOW [START_REF] Aggarwal | Using vaes and normalizing flows for one-shot text-to-speech synthesis of expressive speech[END_REF].

The systems mentioned above have shown significant performance in controlling expressivity, but many approaches use audiobook emotive storytelling style. Besides, they don't work with emotions such as joy, sadness, happiness, fear, anger. However, few approaches addressed the usage of multiple emotions [START_REF] Tits | Visualization and interpretation of latent spaces for controlling expressive speech synthesis through audio analysis[END_REF]- [START_REF] Kulkarni | Transfer learning of the expressivity using FLOW metric learning in multispeaker text-to-speech synthesis[END_REF], [START_REF] Dahmani | Conditional variational autoencoder for text driven expressive audiovisual speech synthesis[END_REF] in TTS. For instance, Deep convolutional TTS (DCTTS) was trained with multiple emotions along with variational inference [START_REF] Tits | Visualization and interpretation of latent spaces for controlling expressive speech synthesis through audio analysis[END_REF]. Recently, metric learning framework was introduced in a parametric TTS system [START_REF] Kulkarni | Deep Variational Metric Learning For Transfer Of Expressivity In Multispeaker Text To Speech[END_REF], [START_REF] Kulkarni | Transfer learning of the expressivity using FLOW metric learning in multispeaker text-to-speech synthesis[END_REF] for transfer of expressivity. In [START_REF] Kulkarni | Deep Variational Metric Learning For Transfer Of Expressivity In Multispeaker Text To Speech[END_REF], the authors proposed to build a recurrent conditional variational autoencoder based acoustic model and is trained using multiclass N-pair loss as additional loss function [START_REF] Sohn | Improved deep metric learning with multi-class n-pair loss objective[END_REF]. This work was further extended by the addition of Inverse Autoregressive Flow (IAF) for implementing an encoder network of the acoustic model [START_REF] Kulkarni | Transfer learning of the expressivity using FLOW metric learning in multispeaker text-to-speech synthesis[END_REF].

The above mentioned approaches indicates that the addition of multiclass N-pair loss increases the perceived expressivity in the target speaker's voice. They have shown promising results for parametric systems but still, depending on the bottleneck step of the duration model for each emotion. In this paper, we present a novel metric learning framework [START_REF] Kaya | Deep metric learning: A survey[END_REF] jointly trained with a Tacotron 2 based end-to-end TTS system [START_REF] Shen | Natural tts synthesis by conditioning wavenet on mel spectrogram predictions[END_REF]. This results in enhancing the latent representation of expressivity for better transfer learning performance. The expressivity information learned by the encoder influences the alignment of synthesized speech generated through attention.

The paper is organized as follows, Section II introduces the proposed architecture and the loss function used for training the end-to-end TTS. Section III describes the speech corpora and the processing of data before training. Section IV presents the experimentation setup. This is followed by Section V, presents the results obtained. Finally, Section VI and Section VII concludes the paper and presents a discussion.

II. PROPOSED ARCHITECTURE

The Tacotron 2 consists of a recurrent sequence-to-sequence feature prediction network that maps character embeddings to mel-scale spectrograms. This recurrent sequence-to-sequence replaces the convolutional block highway gates as proposed in Tacotron system [START_REF] Wang | Tacotron: A fully end-to-end text-to-speech synthesis model[END_REF]. A more detailed explanation of Tacotron 2 architecture is presented in [START_REF] Shen | Natural tts synthesis by conditioning wavenet on mel spectrogram predictions[END_REF].

We extend the state-of-the-art Tacotron 2 model [START_REF] Shen | Natural tts synthesis by conditioning wavenet on mel spectrogram predictions[END_REF] based on sequence to sequence with attention module to implement an end-to-end multispeaker expressive TTS system. In order to work with multiple speakers and expressivity, we modify the Tacotron 2 approach by adding an expressivity encoder and a speaker encoder.

A. General Framework

The proposed architecture takes input as text, which is then converted into a sequence of phonemes. The text encoder processes sequence of phonemes after passing through convolutional layers proceeded by Bidirectional Long Short Term Memory (BLSTM) based recurrent neural network to get z t as a latent representation of text. The reference mel spectrogram is given to the expressivity encoder to extract the latent emotional information as expressive embedding z e . For enabling a multispeaker setting, we have provided the speaker identity to speaker encoder to create embedding, z s . The speaker encoder network maps the speaker index to nonlinear fixed dimensional speaker embedding.

Afterwards, z t , z e , and z s are concatenated and given as input to the location sensitive attention module, as illustrated in Fig. 1. This assists the end-to-end TTS to learn the alignment between the sequence of phonemes and the desired Mel spectrogram. The decoder network is composed of pre-net, BLSTM based recurrent network, and convolutional layer based post-net.

The decoder takes encoder outputs with attention vector as input to predict the Mel spectrogram frame by frame. The output from the previous frame is first passed through the prenet network. The pre-net network consists of fully connected layers with the ReLU activation function. This predicted mel spectrogram from the pre-net and recurrent network is further passed through the post-net network. The post-net network composes of 5 layers of a convolutional network. The postnet improves the overall reconstruction performance of Mel spectrogram by predicting residual to add to the predicted Mel spectrogram.

B. Expressivity encoder

We have experimented with two neural network architectures for the implementation of expressivity encoder namely, Global Style Token (GST) [START_REF] Wang | Style tokens: Unsupervised style modeling, control and transfer in end-to-end speech synthesis[END_REF] and Variational Autoencoder (VAE) [START_REF] Zhang | Learning latent representations for style control and transfer in end-to-end speech synthesis[END_REF]. The GST based expressivity encoder consists of a reference encoder, style attention, and style embedding. The reference encoder maps expressivity of variable length mel spectrogram into a fixed-length vector, which is passed to the style attention layer. This layer applies a multi-head attention module to extract the similarity between reference embedding and each token in style embedding as an output of expressivity encoder. In this work, style embedding z e represents the expressivity as a stylistic factor to learn from the reference embedding.

The second architecture for the expressivity encoder is VAE based which is composed of a reference encoder and two feedforward layers to generate mean and standard deviation of latent variable z e . The reference encoder in VAE generates a hidden output which is passed through the feed-forward layers to obtain latent variable z e . The z e is obtained by using a reparameterization trick applied with mean and standard deviation. VAE based framework suffers from Kublack Leibler (KL) annealing problem in which reconstruction loss is suppressed by KL loss term [START_REF] Bowman | Generating sentences from a continuous space[END_REF]. In KL annealing problem, after few initial epochs, KL divergence term suddenly goes close to zero.

To avoid this, additional weight (close to zero) is multiplied to KL loss and gradually increased over the training epoch.

C. Multiclass N-pair loss

The end-to-end TTS is jointly trained with the multiclass Npair loss function. This assists in predicting Mel spectrogram output in desired emotions. For the transfer of expressivity, we have used pre-computed means of latent variables of each emotion. Thus, during the inference phase, for a given mean of latent variables of emotion, the system transfer expressive attributes to the target speaker's voice. The latent space representation of unclustered emotions may lead to the poor transfer of expressivity. And for better performance of expressivity transfer, we need the tightly bounded representation of the latent variables of emotions. Therefore, we propose a novel deep metric learning framework implemented using multiclass N-pair loss to further enhance the expressivity representation.

Deep metric learning has gained popularity for solving discriminative tasks in computer vision and image processing domain [START_REF] Lin | Deep variational metric learning[END_REF]. The deep metric learning framework assists in the clustering of embeddings by reducing the distance between embeddings of positive classes and increasing the distance between each negative class. The multiclass N-pair loss has shown better performance than triplet loss or contrastive loss by considering embeddings of multiple negative classes [START_REF] Sohn | Improved deep metric learning with multi-class n-pair loss objective[END_REF]. In the training phase, the model needs to reduce the multiclass N-pair loss function as stated in Eq.1, in addition to the reconstruction loss and attention loss.

Loss N-pair = log(1 + N -1 i=1 exp(z e z - i -z e z + )) (1)
In our approach, the multiclass N-pair loss function is applied with respect to the N emotion classes. This loss criteria increases the intercluster distance from N -1 negative samples and decreases the intracluster distance between positive samples and training examples [START_REF] Kulkarni | Deep Variational Metric Learning For Transfer Of Expressivity In Multispeaker Text To Speech[END_REF]. The positive example refers to the latent variables from the same emotion class and negative samples corresponds to the examples of various emotion classes. We have provided the mean of latent variables of emotion for sampling the positive and the negative examples. For N classes, z + is a positive example and z - i is a negative example as stated in Eq.1.

III. DATA PREPARATION

We have used 4 French Female speech synthesis corpora for implementing an end-to-end multispeaker expressive TTS system. The speech corpora used are Lisa neutral speech corpus (approx. 3hrs, in house speech synthesis corpus), SIWIS, neutral speech corpus (approx. 5hrs) [START_REF] Yamagishi | The siwis French speech synthesis database[END_REF], Synpaflex speech corpus (approx. 7hrs) [START_REF] Sini | Synpaflex-corpus: An expressive French audiobooks corpus dedicated to expressive speech synthesis[END_REF], and Caroline expressive speech corpus [START_REF] Dahmani | Conditional variational autoencoder for text driven expressive audiovisual speech synthesis[END_REF]. Caroline's expressive speech corpus consists of 6 emotions namely joy, surprise, fear, anger, sadness, and disgust (approx. 1hr for each emotion). Besides expressive speech, Caroline speech corpus also has neutral speech recorded for approximately 3hrs. Each speech corpus is split into train, validation, and test sets in 80 : 10 : 10 ratio respectively. In Synpaflex corpus, expressive speech samples are available but due to an insufficient number of speech samples for each emotion as well as the unbalanced distribution of emotional speech samples, we have used only the neutral voice of Synpaflex in our work.

We have used a sampling rate of 16000 Hz and extracted Mel spectrograms as acoustic features to be predicted by the end-to-end TTS system. We have applied STFT with an FFT length of 1024, hop length of 256, a window size of 1024, and extracted Mel spectrograms using 80 Mel filters. As input features to the end-to-end TTS, we have used a sequence of phonemes extracted from the text. For French grapheme to phoneme conversion SOJA-TTS tool (developed internally in the team) is used as a front-end text processor.

IV. EXPERIMENTAL SETUP

For training the end-to-end TTS system, we have used the same model parameters as explained in [START_REF] Wang | Style tokens: Unsupervised style modeling, control and transfer in end-to-end speech synthesis[END_REF], [START_REF] Zhang | Learning latent representations for style control and transfer in end-to-end speech synthesis[END_REF], [START_REF] Shen | Natural tts synthesis by conditioning wavenet on mel spectrogram predictions[END_REF] for implementing the Tacotron 2 system and expressivity encoders based on GST, and VAE. We have used a 128 dimensional latent variable of expressivity for both GST and VAE. Before 150K training steps, we have applied the weight of 0.0001 in every 200 steps in order to reduce the KL annealing effect. This weight is increased by 0.00001 after every 500 steps.We have adopted a similar technique for fine-tuning with multiclass N-pair loss, for which until 150K training steps a weight of 0 is applied on the multiclass N-pair loss, and afterwards the weight is increased by 0.001 after every 200 steps.

We have incorporated Waveglow [START_REF] Prenger | Waveglow: A flow based generative network for speech synthesis[END_REF] based neural vocoder for synthesizing speech waveform and trained it on 4 French speech synthesis corpora mentioned in Section III. For evaluating the performance improvement obtained using the addition of multiclass N-pair loss, we have used an end-toend TTS model with GST and VAE as baseline models. We have compared the baseline models with an end-to-end TTS systems trained along with the multiclass N-pair loss for both expressivity encoders, GST, and VAE.

V. RESULTS

In this section, Table I and Table II includes proposed end-to-end TTS (E2E) models, along with evaluation scores obtained using the parametric multispeaker expressive TTS as stated in [START_REF] Kulkarni | Deep Variational Metric Learning For Transfer Of Expressivity In Multispeaker Text To Speech[END_REF], [START_REF] Kulkarni | Transfer learning of the expressivity using FLOW metric learning in multispeaker text-to-speech synthesis[END_REF]. We have presented the evaluation scores for parametric TTS as recurrent conditional variational autoencoder (RCVAE), RCVAE N-pair, and inverse autoregressive flow (IAF) N-pair.

A. Objective evaluation

We have conducted an objective evaluation using Mel Cepstrum Distortion (MCD), F0 Root Mean Squared Error (F0 RMSE), and Voiced-Unvoiced error (VUV) between reference speech samples and proposed E2E TTS systems. The objective evaluation results are presented in Table I GST N-pair model has shown superior performance compared to other models. We opt for subjective evaluation to measure the performance of transfer of expressivity, due to the unavailability of reference emotional speech samples for Lisa, Siwis, and Synpaflex speech corpora.

B. Subjective evaluation

At first, we have evaluated an end-to-end (E2E) multispeaker expressive TTS systems using Mean Opinion Score (MOS) [START_REF] Streijl | Mean opinion score (MOS) revisited: methods and applications, limitations and alternatives[END_REF] based listening test. In this work, we used the absolute category ranking scale. Each listener had to assign a score for synthesized speech utterance between scale 1 to 5 considering the intelligibility, naturalness, and quality of speech utterance. Suppose the speech quality is bad the listener will then assign the score 1 and if the speech quality is excellent then the listener will assign the score 5. Each listening test consists of 10 randomly selected speech files from the test set for each model. A total of 14 French listeners participated in this MOS test and results are displayed in Table II with an associated 95% confidence interval.

The main goal of this work is to transfer the emotion as expressive attributes to the target speaker's voice without altering the speaker's voice characteristics. As there is no possible way to extract quantitative results for evaluation of transfer of expressivity without reference to expressive speech samples, we opt for speaker MOS and expressive MOS as a qualitative measure for expressivity transfer.

Similarly, for speaker MOS as well, we instructed the listeners to assign the score between 1 (bad) and 5 (excellent) to the speech samples based on the speaker similarity between reference speaker speech and synthesized expressive speech. Likewise, for expressive MOS, listeners are directed to provide scores between 1 (bad) and 5 (excellent) depending on how synthesized speech utterance resembles the expressivity given in the reference speech utterance. A total of 14 French listeners performed both listening tests mentioned above, where each listener scored 18 speech utterances for each speaker-emotion pair and model. The results obtained through expressive MOS and speaker MOS are presented in Table II with associated 95% confidence intervals.

VI. DISCUSSION

From Table II, the obtained MOS scores for each E2E model are consistent with the objective evaluation score in Table I. The E2E GST N-pair model outperformed all models, thus usage of multi-head attention assists in speech synthesis. The E2E GST N-pair model performance shows that the addition of N-pair to expressivity encoder boosts the model performance which can also be seen with the performance of VAE based expressivity encoder.

The speaker MOS score of the E2E VAE N-pair model is higher than the other E2E models. The IAF N-pair based parametric TTS performed slightly better than the E2E VAE N-pair model for retaining speaker attributes. The IAF N-pair model uses x-vector based speaker embedding for creating speaker representation [START_REF] Kulkarni | Transfer learning of the expressivity using FLOW metric learning in multispeaker text-to-speech synthesis[END_REF]. This results in better speaker representation than when only speaker identity is provided for creating speaker embeddings.

The E2E VAE N-pair model obtains the highest expressive MOS score. This shows that the E2E VAE N-pair model can generalize a better emotion latent space using an expressivity encoder than the E2E GST N-pair model.

The proposed E2E models learn the duration information as alignment derived from attention vector as opposed to parametric TTS, where explicit duration model is required. From Table II, MOS scores for parametric TTS range between 2 and 3, while E2E models have greater than 3. The parametric TTS systems conduct expressivity transfer in acoustic space only, which lacks the interpolation in duration prediction. Thus, E2E TTS system not only influence the prosody of synthesized speech but also the alignment of synthesized speech for each emotion.

VII. CONCLUSION

We proposed to use multiclass N-pair loss on latent representation extracted using expressivity encoder to derive emotion as latent information. During the inference phase, we investigated the transfer of expressivity without the explicit need of reference Mel spectrogram. We used pre-computed means of latent variables of each emotion for expressivity transfer. The obtained results show that the performance of expressivity transfer is significantly improved "with" the addition of N-pair loss in comparison to "without" use of Npair loss. For transfer of expressivity VAE based expressivity encoder generalizes emotion representation better than GST. To our knowledge, the presented work is the first approach to use metric learning in an end-to-end multispeaker TTS system.

VIII. ACKNOWLEDGEMENT

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).

Fig. 1 .

 1 Fig. 1. End-to-end multispeaker expressive text to speech system

  . From TableI E2E

TABLE I OBJECTIVE

 I EVALUATION FOR AN END-TO-END TTS SYSTEM

	Model	MCD F0 RMSE VUV error
	E2E GST	5.12	24.10	10.41
	E2E VAE	5.29	24.24	10.72
	E2E GST N-pair	4.71	23.63	8.50
	E2E VAE N-pair	4.82	23.70	9.10

TABLE II SUBJECTIVE

 II EVALUATION OF AN END-TO-END TTS SYSTEM

	Model	MOS	Speaker	Expressive
			MOS	MOS
	RCVAE	2.62 ± 0.5	2.40 ± 0.3	1.53 ± 0.3
	RCVAE N-pair	2.97 ± 0.4	2.86 ± 0.3	1.93 ± 0.2
	IAF N-pair	3.02 ± 0.4 2.93 ± 0.4 2.03 ± 0.3
	E2E GST	3.51 ± 0.3	2.57 ± 0.2	3.05 ± 0.2
	E2E VAE	3.38 ± 0.4	2.71 ± 0.3	3.12 ± 0.2
	E2E GST N-pair	3.72 ± 0.4	2.65 ± 0.2	3.15 ± 0.4
	E2E VAE N-pair	3.47 ± 0.3	2.83 ± 0.3 3.33 ± 0.3


O Ezvan

X Zeng

R Ghanem

B Gencturk

Dynamic substructuring for multilevel spent nuclear fuel containers with high modal density

The vibration analysis of a multilevel spent nuclear fuel canister loaded with 68 fuel assemblies (FA) is considered. The multilevel aspect leads to large numbers of degrees of freedom (DOF) and vibration eigenmodes. The Craig-Bampton (CB) substructuring technique is well suited to perform the modal analysis in presence of few interface DOF between the components and the structural levels, as well as for tackling the high modal density. The Shift-Invert Lanczos solver is used for solving the large-scale CB generalized eigenvalue problem (GEP), whose large dimension results from the numerous component modes. A block factorization based on the Schur complement that exploits the sparsity of the CB matrices enables considerable speed-up. A CB model for the FA is also implemented as it exhibits two structural levels. Since many local vibration modes do not contribute much to the response, a strategy for ranking the importance of the component modes is developed, so that a significant portion of them are removed and the CB GEP lightened.

Introduction

This paper deals with the linear vibration analysis of a spent nuclear fuel (SNF) canister loaded with 68 identical fuel assemblies (FA). Each FA gathers 92 fuel rods that are wrapped together and maintained in parallel configuration by spacer grids. The loaded SNF canister being considered in vertical configuration, each FA is inserted in one of the basket cells of the canister, maintained by localized attachments at the top and simply supported by the bottom plate of the cylindrical canister. The detailed geometry of the FA requires a finer mesh in comparison to that of the canister and its basket. The structure is characterized by the presence of three distinct structural levels, namely (i) a main body made up of the canister and the basket as the upper level, (ii) the 68 FA as the intermediate level, and (iii) the 68 × 92 = 6256 fuel rods as the lower level. There is the repetition of a very large number of interconnected components, each of which having its own resonant vibrations. As a result, the multilevel structure exhibits a very large number of modes, among which many consist of local vibrations of small components. The vibration analysis relies on computing nearly half a million modes using a finite element (FE, [START_REF] Zienkiewicz | The finite element method[END_REF]) model with more than 130 million degrees of freedom (DOF). Fortunately, the localized attachments between the several components and scales allows an advantageous domain decomposition to be considered, with small interfaces. The vibration modes are the eigenvectors of the generalized eigenvalue problem (GEP) involving the mass and stiffness FE matrices [START_REF] Craig | Fundamentals of structural dynamics[END_REF][START_REF] Ohayon | Advanced computational vibroacoustics -reduced-order models and uncertainty quantification[END_REF]. For several reasons, the Craig-Bampton (CB) substructuring technique [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF][START_REF] De Klerk | General framework for dynamic substructuring: history, review and classification of techniques[END_REF] is considered as an approximate solver for this GEP. First, the substructural modes of the FA need not be computed 68 times. Second, the small interfaces between the components lead to a small number of interface coupling DOF. Third, the iterative solvers such as Shift-Invert Lanczos (SIL, [START_REF] Wilkinson | The algebraic eigenvalue problem[END_REF][START_REF] Parlett | The symmetric eigenvalue problem[END_REF][START_REF] Grimes | A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems[END_REF]) or the Subspace Iteration Method [START_REF] Bathe | Solution methods for eigenvalue problems in structural mechanics[END_REF][START_REF] Bathe | The subspace iteration method -revisited[END_REF][START_REF] Bathe | Finite element procedures[END_REF], although very efficient general-purpose solvers, turn out to be very time consuming in the present case.

The CB method consists in approximating the system modes as linear combinations of the component modes and interface static constraint modes through Galerkin projection of the FE GEP. The number of static constraint modes is given by the number of interface DOF, whereas the number of component modes is generally adjusted by considering a cutoff frequency. Even though the number of static constraint modes is small, the dimension of the CB GEP obtained through Galerkin projection is very large, as a result of the large number of component modes. Usually for the CB method, the CB GEP is solved directly for all the eigenpairs using for instance Householder tridiagonalization, with a negligible computational cost. Unfortunately, in the present case, the CB GEP is a large-scale GEP with sparse matrices, and for which a large portion of the eigenpairs are required. Although the SIL solver is better suited to solving a few eigenpairs of a sparse GEP, spectrum slicing with multiple SIL is a possible strategy for solving many eigenpairs. The strategy consists in dividing the frequency band into several slices and considering eigenvalue shifts given by the center of the slices [START_REF] Ericsson | The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems[END_REF]. Independent eigenvalue computations with the different shifts are then run in parallel, for which the eigenvectors are orthogonal by construction (upon convergence). For each shift, the matrix to be factored for use in the iterative process is different. However, the sparsity of the CB matrices allows an efficient factorization. With a small number of boundary (interface) DOF and a large number of substructure (modal) DOF, a block factorization based on the Schur complement (see [START_REF] Saad | Distributed Schur complement techniques for general sparse linear systems[END_REF][START_REF] Zhang | The Schur complement and its applications[END_REF]) of the interface coordinates enables considerable gains. Thus, the CB GEP is solved by multiple SIL with block factorization by Schur complement.

The accuracy of the CB approximation of the system modes depends on the number of component modes kept. Several criteria can be used to quantify the error perpetrated by the CB approximation, such as convergence of the eigenvalues and modal error. However, presumably due to the high modal density of the multilevel structure, these criteria turned out to lack robustness. The objective being to calculate the vibration response, given by frequency response functions (FRF) between given input excitation DOF and output observation DOF, the CB error is quantified through an FRF error criterion. For one given FRF, it is defined as the root mean square error in decibel (RMS-dB) of the modulus of the response. It should be noted that the exact response is required for computing such an error measure, and that it is not available. Since the response can be written as modal superposition and since numerous modes have negligible contribution, it is proposed to substitute the true FRF with the FRF restrained to a subspace spanned by a small subset of dominant modes. Such FRF spanned by a subset of modes is referred to as "partial FRF". The few dominant modes are recomputed exactly by using the SIL solver with block factorization by Schur complement according to the domain decomposition (this differs from the Schur complement factorization of the CB matrices). Then the partial FRF can be compared (CB approximation against FE approximation). It has been shown that this procedure allows to characterize the accuracy of the CB approximation with a good confidence.

Once its accuracy is quantified and deemed sufficient, the CB model can be considered as a reference. Then, it is possible to decrease the cutoff frequency without altering the accuracy. This is explained by the fact that the cutoff frequency required for obtaining rank correspondence between the exact and CB eigenvalues in the validation step above can be very high due to the small eigenvalue separation (caused by the very high modal density), which needs to remain larger than the relative eigenvalue error (otherwise, mode swap is susceptible to occur, in which case the exact and CB partial FRF may not coincide). Since the structure is multilevel, another CB step (inner CB), nested in the first one, is well indicated. The FA exhibits numerous components (the fuel rods) attached to the main body by localized connections (through the spacer grids). The cost for computing the FA modes (component modes of the outer CB) is thus decreased.

The most computationally demanding operation is the large-scale CB GEP. Its large dimension is due to the numerous component modes, which correspond to local vibrations of small components. It has been observed that many of the system modes made up of such local vibrations do not contribute much to the response of the canister. In this research, the focus is on the prediction of the vibration response of the canister exterior. Therefore, the removal of component modes with local vibrations is investigated in order to decrease the dimension and cost of the CB GEP while keeping a satisfactory level of accuracy. For doing so, a substructural (component) modal importance measure is introduced, defined as the RMS-dB error re-sulting from the sole removal of the substructural mode for which the error is to be characterized. However, the calculation of the RMS-dB error (with respect to the CB model) requires calculating the FRF with the CB model. It is proposed to do so with coarse frequency sampling and observation grid, using direct numerical simulation (frequency-by-frequency) involving a dynamic stiffness matrix that can again be factored by Schur complement. Then, it is shown that the approximate response resulting from the removal of one given mode can be calculated through a rank-1 modification of the Schur complement. Using the Scherman-Morrison-Woodbury matrix identity, the response can be computed very efficiently. This enables a precise yet affordable characterization of the substructural modal importance. Then, the least important modes can be removed and the induced error examined with respect to the cost reduction.

The paper is organized as follows. In Section 2 some well known concepts for modal analysis in computational structural dynamics are briefly described, including modal projection, domain decomposition, Craig-Bampton method, and Shift-Invert Lanczos approach with spectrum slicing. In Section 3 the methodology proposed for performing an efficient multiscale modal analysis is presented. Section 4 presents developments toward substructural mode selection that allow the overall cost to be significantly reduced. Finally, Section 5 presents numerical results and Section 6 draws conclusions. The theoretical developments and numerical results can be found with more details in [START_REF] Ezvan | Multiscale modal analysis of fullyloaded spent nuclear fuel canisters[END_REF][START_REF] Ezvan | Dominant substructural vibration modes for fully-loaded spent nuclear fuel canisters[END_REF].

Modal analysis in computational structural dynamics 2.1 Modal projection

The structural dynamics is described in the frequency domain by the discretization U(ω) of the displacement field by using the FE method. For all circular frequency ω belonging to a given frequency band of analysis B = 2π×]0, f u ] (f u is the upper bound in Hz), the complex vector U(ω) is given by

-ω 2 [M] + iω[D] + [K] U(ω) = F(ω) , (1)
in which [M], [D], and [K] are respectively the mass, damping, and stiffness matrices and where vector F(ω) is the discretization of the external forces. A reduced representation of the response can be obtained through linear combination of the vibration eigenmodes ϕ such that

U(ω) = ne α=1 ϕ α q α (ω) , (2)
in which α denotes the mode number (the modes are sorting according to their eigenfrequency or eigenvalue, as described below) and n e denotes the number of modes kept, while quantity q α (ω) constitutes a so-called "generalized" coordinate. In matrix form, one can write

U(ω) = [Φ]q(ω) with [Φ] the modal matrix whose size is (N × n e )
with N the number of DOF. The vibration modes ϕ are the eigenvectors satisfying the GEP

[K]ϕ = λ[M]ϕ , (3)
with λ the associated eigenvalue and f = √ λ/2π the eigenfrequency in Hz.

Domain decomposition

The structural domain can be partitioned. The DOF that are connected to several subdomains or substructures are referred to as boundary DOF and denoted by letter B. The other DOF belong to the interior of the substructures are denoted by letter S. These substructure DOF are not coupled with any other DOF than those present in the same substructure. Therefore, the coupling blocks between substructures are zero and consequently, a DOF ordering according to the domain decomposition leads to block-diagonal submatrices

[K SS ], [M SS ], and [D SS ].
Based on this domain decomposition, the linear system [K]x = f can be parallelized over the substructure, as described hereinafter. Let x (resp. f) be partitioned based on the subvectors x B and x S (resp. f B and f S ) associated with the boundary DOF and the substructure DOF. Then, it can be shown that introducing the Schur complement

[S K BB ] = [K BB ] -[K BS ][K SS ] -1 [K SB ] leads to the solution x S = [K SS ] -1 (f S -[K SB ]x B ) , x B = [S K BB ] -1 f B -[K BS ][K SS ] -1 f S . (4)

Craig-Bampton method

The above equations can be obtained by projecting the linear equations onto a square matrix [B] such that

[B] = I S t SB 0 I B , (5)
in which

[t SB ] = -[K SS ] -1 [K SB ]
is the matrix of the static constraint modes that represent the static reactions of the interfaces between the substructures. The CB technique first introduces this change of coordinates and then, for dimension reduction, introduces a projection on the component modes

[Φ S ] gathered in a matrix [R] such that [R] = Φ S 0 0 I B , (6)
where matrix [Φ S ] is block diagonal with each block containing the modal matrix of its corresponding substructure. The CB matrices are obtained by Galerkin projection over the matrix

[V ] = [B][R] and yields the CB GEP [K]q = λ[M ]q (7)
whose eigenvectors q yield an approximation ϕ = [V ]q to the eigenvectors ϕ of the GEP of Eq. ( 3). The accuracy and the cost of the CB procedure are related to the modal truncations for the substructures. Denoting as n J the number of modes kept for substructure number J, as N the number of substructures, and as N B the number of boundary DOF, the dimension ν of the CB model is ν = N B + n S with n S = N J=1 n J .

Shift-Invert Lanczos approach

For solving a GEP of the type of Eq. ( 3) one possible strategy is to introduce an eigenvalue shift σ ≥ 0 and consider shifting as

([K] -σ[M]) ϕ = (λ -σ) [M]
ϕ and then inverting (useful for eigenvalue separation) as

([K] -σ[M]) -1 [M]ϕ = 1 λ -σ ϕ . (8)
The largest eigenvalues of the standard eigenvalue problem

[A]ϕ = µ ϕ with [A] = ([K] -σ[M]) -1 [M]
and µ -1 = λ -σ correspond to the eigenvalues λ that are closest to the shift σ. Due to the orthogonality of the eigenvectors ϕ, two sets of eigenvectors (modes) computed from two different shifts are mutually orthogonal. This is verified as long as the eigensolutions are converged. This way, a large number of modes can be computed by considering several shifts. Thus, the eigenvalue calculation can be parallelized and the complexity is kept linear with respect to the number of modes sought. For each eigenvalue analysis of the type of Eq. ( 8), the matrix

[H] = [K] -σ[M]
has to be factored for solving linear systems in the iterative Lanczos process. In this work, an LDL factorization (see [START_REF] Golub | Matrix computations[END_REF])

[P] T [H][P] = [L][D][L]
T is considered, with [P] a permutation matrix, [L] a lower-triangular matrix, and [D] a block-diagonal matrix with block size of 1 or 2. Dividing the frequency band (spectrum slicing) to perform the eigenvalue analysis of Eq. ( 3) requires a good knowledge of the modal density, to balance the number of modes per shift as well as to avoid computing the same modes twice. In this work, the matrix size makes the LDL factorization costly in terms of memory and CPU usage. Therefore, the domain decomposition of Section 2 is considered and the LDL factorization is performed blockwise based on the Schur complement (see Section 3.1).

Multiscale modal analysis 3.1 Block factorization based on Schur complement

Using domain decomposition, the symmetric indefinite matrix [H] can be factored as

[P] T [H][P] = [L][D][L] T
with the following construction:

[L] = L S 0 P T B H BS H -1 SS P S L S L B , [D] = D S 0 0 D B , [P] = P S 0 0 P B , (9)
in which the matrices with index S satisfy the block LDL factorization

[P S ] T [H SS ][P S ] = [L S ][D S ][L S ] T
(this factorization can be parallelized over the substructures according to the block-diagonal structure of [H SS ]) and where the matrices with index B satisfy the LDL factorization

[P B ] T [S H BB ][P B ] = [L B ][D B ][L B ]
T of the Schur complement. Such block LDL factorization as described in (9) will be used in several instances of the proposed methodology. It should be noted that the inertia count I([H]) (number of negative eigenvalues of matrix [H]) associated with a shift σ can be obtained as

I([H]) = I([D]) = I([D B ]) + I([D S ]) = I([D B ]) + N J=1 I([D J ]) . ( 10
)

Large-scale Craig-Bampton generalized eigenvalue problem

For performing the eigenvalue analysis of Eq. ( 3) in an affordable manner, the CB GEP of Eq. ( 7) is considered instead. In order to obtain an accurate CB model, a high frequency truncation is considered that leads to a large value for the number n S of substructural (or component) modes. The small interfaces between components yield a small number N B of boundary DOFs. The CB matrices are written as

[K] = Λ S 0 0 K BB , [M ] = I S M T BS M BS M BB , (11)
in which [I S ] is the identity matrix of dimension n S and [Λ S ] is the diagonal matrix of the eigenvalues of the substructures, and where the submatrices [K BB ] and [M BB ] are rather dense and the submatrix [M BS ] is sparse with dense blocks (their expression can be deduced from the equality

[K] = [V ] T [K][V ] or [M ] = [V ] T [M][V ]
). Due to the large value of n S related to the high modal density arising from the local vibrations of the numerous components, the dimension ν of the CB matrices [K] and [M ] can be several million. In this context, the GEP of Eq. ( 7) cannot be solved by the usual direct solvers by transformation methods such as the Householder tridiagonalization. Instead, the only tractable possibility is to consider the solvers for large-scale sparse GEP such as SIL or SIM. This is possible because the CB matrices exhibit quite a high sparsity, that is amplified by the large ratio n S /N B (large number of substructural modes and small number of boundary DOF). In addition, the block structure allows an LDL factorization just like that of Eq. ( 9) and, given that [Λ S ] and [I S ] are diagonal, significant gains can be obtained by doing so. For following the SIL approach of Section 2.4, the LDL factorization

[P ] T [H][P ] = [L][D][L] T of [H] = [K] -σ[M ] is considered: [L] = I S 0 -σP T B M BS (Λ S -σI S ) -1 L B , [D] = Λ S -σI S 0 0 D B , [P ] = I S 0 0 P B , (12)
I([H]) = I([D]) = I([D B ]) + I([Λ S ] -σ[I S ]) = I([D B ]) + N J=1 I([Λ J ] -σ[I J ]) . ( 13
)

Accuracy verification

The modal density can be estimated using several shifts σ and computing, for each of them, the inertia count I([D B ]) through the LDL factorization of the Schur complement [S H BB ]. Since the number N B of boundary DOF is small, this is quite efficient. Then, several slices are introduced and parallel SIL is undertaken. For each GEP, the matrix factorization requires the LDL factorization of the Schur complement [S H BB ]. Then, the linear solves involve thin matrix products involving matrices with similar structure as [M BS ]. The CB GEP can thus be solved very efficiently for a very large number of eigenpairs. Very large truncations are first considered, to obtain a CB model that can be considered as error-free.

Concerning the error quantification, the usual error measures based on eigenvalue convergence or the modal error did not turn out to be adequate for the multilevel structure considered. With the modal error, too large values were obtained for some modes while considering very high truncations. Furthermore, because of the very high modal density, the relative eigenvalue error could remain as high as the relative difference between two distinct eigenvalues. In fact, the eigenvalues are susceptible to swap even though the CB model is already very accurate. Since the purpose of the computational model is to calculate the FRF, the error measure is defined through the FRF error. Let U ij (ω) denote the entry i, j of the true FRF matrix and let U ij (ω) denote one approximation of it, for instance obtained with the CB model. The associated error measure ij for this particular FRF is defined as

ij = 1 |B| B (u ij (ω) -u ij (ω)) 2 dω , (14)
= µ + k σ , µ = 1 N 2 I N I i=1 N I j=1 ij , k > 0 , σ 2 = 1 N 2 I -1 N I i=1 N I j=1 ( ij -µ ) 2 , ( 15
)
in which k is a parameter that indirectly controls the probability of obtaining an error greater than when randomly considering one FRF with error ij . The exact FRF U ij (ω) are not available and consequently, this error measure cannot directly be used to quantify the accuracy of the CB model. Since the response can be written as modal superposition and since numerous modes have negligible contribution, it is proposed to substitute the true FRF with the FRF restrained to a subspace spanned by a small subset of dominant modes. The dominant modes are defined as those for which the error induced by their sole removal is the greatest. Such FRF spanned by a subset of modes is referred to as "partial FRF". This way, it is possible to calculate the true "partial FRF" by computing the subset of dominant modes by using a solver as accurate as machine precision. To do so, the SIL approach of Section 2.4 is used along with domain decomposition for the linear solver as described in Section 2.2. The number n d of dominant modes should be such that n d n e and adjusted such that the eigenvalue computation of Eq. (3) be affordable. Then, the use of the error measure defined in Eqs. ( 14) and ( 15) allows for controlling the accuracy of the n d dominant modes calculated through the CB GEP of Eq. [START_REF] Parlett | The symmetric eigenvalue problem[END_REF]. It should be noted that these dominant modes are very much likely to cover the whole frequency band. These dominant modes do not have any reason to have better accuracy than all the other modes that have not been verified by such direct comparison. To make sure of it, one possibility is to examine the other error measures (modal error and eigenvalue convergence) and verify that the dominant modes do not behave differently than the others.

Loose truncation and multilevel Craig-Bampton

The FRF error described in Section 3.3 required that the modes recomputed with SIL correspond to the same modes computed by the CB approach. Otherwise the error due to the mismatch between the FRF would not be due to the inaccuracy of the CB model, and the accuracy verification would be biased and most likely exhibiting too large errors. The mode correspondence was thus necessary for the accuracy certification described in Section 3.3 and required large truncations (large value for n S ). Nevertheless, the CB model does not necessarily need that large a value for n S . The CB model for which the accuracy has been verified in Section 3.3 is used as a reference and then the FRF error (not that of the partial FRF but that of the true FRF) is used to study the convergence of the CB model with respect to the truncation of the substructural (or component) modes.

In fact, the accuracy verification of Section 3.3 can be seen as an offline step necessary for preparing an efficient procedure to perform the dynamic analysis for the same computational model but altered with few changes such as damage or randomness. Beside the loose truncation compared to that used in the accuracy verification step of Section 3.3, several other approximations are examined to decrease the cost of the computational model while keeping the same level of accuracy.

Since the multilevel structure presents three structural levels, it is also beneficial to consider a second CB approximation (referred to as "inner CB"), nested in the first one (referred to as "outer CB"), that is to say devoted to computing the component modes of the outer CB. The component modes of the outer CB are the FA modes and the modes of the so-called frame substructure constituted of the canister and its basket. The inner CB only considers the CB model of a FA (the 68 FA are nominally identical) as it is assumed that the frame substructure is not subject to change. Just like the FA are connected to the upper scale through localized attachments with the basket and the canister bottom plate, the fuel rods are connected to the FA skeleton through localized attachments with the spacer grids. The same methodology is used for performing this inner CB step.

Substructural mode selection

For one different instance of the computational model, that is to say e.g., either for one choice of the set of mechanical parameters describing the computational model or for one damage scenario, the FA modes and static constraint modes need to be recomputed (no change is applied to the canister and its basket). As a consequence, the CB model is changed and the CB GEP of Eq. ( 7) needs to be solved again. The main cost is that of the CB GEP and for decreasing it, the strategy proposed is to remove the FA modes that do not contribute much to the CB model. From one instance of the computational model to the other, the importance or dominance of the FA modes varies, in general. Therefore, it needs to be calculated "online", that is to say for each instance of the computational model. Following the same definition as in Section 3.3, the dominant FA modes are defined as those for which the error induced by their sole removal is the greatest. Nevertheless, in this case, such calculation must be less time consuming than the time that is gained by dealing with a CB GEP with a reduced dimension resulting from the removal of FA modes.

For quantifying the error due to removing one given FA mode, the FRF error is again used as in Section 3.3. The true FRF is that obtained with the CB model without FA mode removal and the approximate FRF is that with single FA mode removal. Again, the true FRF is not available. The strategy thus proposed is to calculate the true FRF using direct numerical simulation using the dynamic stiffness matrix

[H(ω)] = -ω 2 [M ] + iω[D] + [K] , for which the damping model is hysteretic such that [D] = 2ξ(ω)
ω [K] with ξ(ω) a frequency-dependent damping ratio. This way, the dynamic stiffness matrix exhibits the same sparsity as the shifted matrix [H] of the CB GEP of Eq. ( 7). The approximate FRF then considers the same damping model and is obtained in considering the removal of one FA mode (thus decreasing the dimension of the CB model by one unit). Let [M βA ] be the matrix similar to [M BS ] but associated with the FA substructure and for which the rows are restricted to its N β boundary DOF. The matrix

[s H,A ββ (ω)] = [M βA ][H AA (ω)] -1 [M βA ] T with [H AA (ω)] = -ω 2 [I A ] + (1 + 2iξ(ω)) [Λ A ]
(dynamic stiffness matrix of the FA) includes the contribution of the FA modes to the Schur complement [S H BB (ω)] of the dynamic stiffness matrix [H(ω)] of the complete system. Let v β denote the column vector of matrix [M βA ] corresponding to the FA mode to be removed.

Then it can be shown that the contribution [ s H,A ββ (ω)] of the remaining FA modes to the Schur complement

[S H BB (ω)] is given by [ s H,A ββ (ω)] = [s H,A ββ (ω)] - 1 h(ω) v β v T β , (16)
in which h(ω) = -ω 2 +(1 + 2iξ(ω)) λ A with λ A the eigenvalue of the removed FA mode. It can be deduced that the modified Schur complement [ S H BB (ω)] resulting from the single FA mode removal can be written as

[ S H BB (ω)] = [S H BB (ω)] + w 4 h(ω) v B v T B , (17)
in which vector v B is the extension of vector v β with zeros in the entries corresponding to boundary DOF of other substructures. It can be seen that the single mode removal comes down to a rank-1 modification of the Schur complement. The Sherman-Morrison-Woodbury formula is then used to efficiently deduce the inverse of the modified Schur complement based on the inverse of the original Schur complement. This is quite useful, considered that there are as many modified Schur complements to consider as the number of FA modes (they are each removed one at a time). It can be shown that then inverse [ S H BB (ω)]

-1

of the modified Schur complement can be written as

[ S H BB (ω)] -1 = [S H BB (ω)] -1 - 1 a(ω) p B (ω)p T B (ω) , (18)
in which the frequency-dependent quantities a(ω

) = h(ω) ω 4 + v T B p B (ω) and p B (ω) = [S H BB (ω)] -1 v B then
allow the modified FRF to be obtained.

For each mode of each FA and for a given sampling of frequencies, the approximate FRF are computed and used along the true FRF to calculate the error measure defined in Eqs. ( 14) and [START_REF] Ezvan | Multiscale modal analysis of fullyloaded spent nuclear fuel canisters[END_REF]. The importance of a given FA mode is defined by the value of the error measure obtained through the sole removal of this substructural mode. Thus, the FA modes can be sorted according to their importance. Then, a given percentage of the most dominant FA mode is kept and the CB GEP solved. A convergence analysis is carried out to determine an adequate percentage of substructural mode filtering, with respect to both accuracy and efficiency.

Numerical results

The frame substructure (canister and its basket) has N F = 2,209,084 DOF, one FA has N A = 1,935,837 DOF, and there is a total of N B = 1140 boundary DOF. Therefore, the FE model is of dimension N = N B + N F + 68 × N A = 133,847,140 DOF. Preliminary calculations pertaining to the offline stage are carried out. Among them, the first n F = 211,057 modes of the frame substructure are calculated, as well as the first n A = 46,383 FA modes. These numbers correspond to a frequency truncation of 20 kHz, given that the modal analysis is to be performed up to 1.2 kHz. Such a high truncation is considered for obtaining a reference computational model. The static constraint modes for both these structures are quite cheap to compute, as there are only N B = 1140 of them for the frame substructures and N β = 27 for the FA. Then, the CB GEP of Eq. ( 7) is solved for all the eigenpairs below 1.2 kHz: there are n e = 458,910 of them. The dimension of the CB GEP is ν = N B + n F + 68 × n A = 3,366,241. Even though it exhibits higher numbers than for the frame substructure, the calculation is faster, thanks to the efficient block factorization.

The n e = 458,910 CB system modes are removed one at a time and the error measure defined in Eqs. ( 14) and ( 15) is calculated for each of them. This allows for the modes to be ranked according to their dominance. For accuracy certification as described in Section 3.3, the n d = 200 most dominant modes are considered.

For each of them, the approximate eigenvalue given by the CB model with 20 kHz truncation is available. For each of these eigenvalues, the corresponding eigenpair is recomputed up to machine precision using the SIL approach with domain decomposition as described in Sections 2.2 and 2.4, by using a frequency shift equal to the approximate eigenvalue. Such a precise shit is necessary because of the very high modal density. The matrix factorization necessary for the Lanczos iterative process also gives the inertia count as a byproduct, which allows the rank correspondence to be verified. Then the "partial FRF" restrained to these n d = 200 dominant modes are computed with both the exact modes and the CB modes, which then allows to calculate the FRF error defined in Eqs. ( 14) and [START_REF] Ezvan | Multiscale modal analysis of fullyloaded spent nuclear fuel canisters[END_REF]. Using k = 4 (the choice for this value is fixed for the rest of this paper), the decibel error is of 0.46 dB, which is very low. Therefore, the CB model with 20 kHz truncation can be considered as a reference.

Then, a convergence analysis with respect to the truncation frequency is carried out, which shows that a 4 kHz truncation is sufficient to reach a decibel error of 0.38 dB with respect to the reference computational model (CB model with 20 kHz truncation frequency). There are only n F = 16,623 frame modes and n A = 14,064 FA modes up to 4 kHz and consequently, the dimension of the CB model is decreased to ν = 974,115.

Then, the FA modes are approximated by introducing an inner CB model for the FA. Using a 6 kHz inner truncation allows for obtaining a decibel error of 0.41 dB, that is to say no noticeable error is introduced. Due to this CB approximation, the number of FA modes slightly decreases to n A = 14,047. Applying the methodology presented in Section 4, these modes are sorted according to their importance. For that, a coarse frequency sampling of 1000 frequency points instead of 3000 frequency points as used everywhere else in this paper is considered. Moreover, the number of FRF considered is decreased, such that the size of several of the frequency-dependent matrices involved in the substructural mode selection are reduced. Despite these approximations, the methodology delivers satisfying results, as a removal of 75% of the least important FA modes leads to a decibel error of 0.45 dB, which is still very low. Thanks to this filtering, the size of the CB model is decreased to ν = 256,579 with n A = 3512. As a remark, the sole removal of the one most dominant FA mode gives an error of 7.43 dB, which is quite significant.

Regarding the computational cost, for the FA mode selection, the elapsed time for the true FRF is ∼ 2 minutes on one compute node devoted to one FA, and the elapsted time for computing the 68 × 14,047 = 955,196 FRF with FA mode removal is ∼ 12 minutes (therefore, it is five orders of magnitude faster). The mode filtering allows the CB GEP to be solved in ∼ 21 minutes with dimension ν = 256,579 instead of being solved in ∼ 118 minutes with dimension ν = 974,115. It should also be noted that the computational cost per mode is four orders of magnitude less for the reference CB GEP with 20 kHz truncation than for the SIL approach with domain decomposition as described in Sections 2.2 and 2.4. The proposed methodology allows the vibration analysis to be performed in about ∼ 40 minutes whereas it is nearly intractable with a standard approach.

Conclusions

For the vibration analysis of a multilevel structure characterized by pseudo-periodicity arising from the repetition of numerous identical components, an efficient solver combining Craig-Bampton substructuring and Shift-Invert Lanczos solver with spectrum slicing allowed to tackle the high computational cost induced by the large number of DOF and the very large modal density characterized by nearly half a million modes. The proposed methodology took advantage of the small number of interface DOF between components, which allowed efficient matrix factorizations. In addition, capturing the relevant physics by removing 75% of the least significant substructural modes allowed the computational cost to be further reduced. The proposed methodology enables the practical simulation of the complex multilevel structure for studying various instances of damage or variation of the mechanical properties.

  in which [P B ], [L B ], and [D B ] satisfy the LDL factorization [P B ] T [S H BB ][P B ] = [L B ][D B ][L B ] T of the Schur complement [S H BB ]. It should be noted that the inertia count I([H]) (number of negative eigenvalues of matrix [H]) associated with a shift σ can be obtained as

  in which u ij (ω) = 20 log 10 |U ij (ω)| and u ij (ω) = 20 log 10 | U ij (ω)| are the moduli in decibel (dB) scale, and where |B| denotes the length of the frequency band of analysis. Error ij is the RMS-dB error. The global error measure is then given by

Acknowledgements

This research was supported by the United States Department of Energy through the Nuclear Energy University Program under the Contract No. DE-NE0008529. The findings presented herein are those of the authors and do not necessarily reflect the views of the sponsor.