Eduardo H Sanchez
email: eduardo.sanchez@irt-saintexupery.com

Mathieu Serrurier
email: mathieu.serrurier@irt-saintexupery.com

Mathias Ortner
email: mathias.ortner@airbus.com

Mutual Information Measure for Image Segmentation Using Few Labels

Keywords: Mutual information maximization, Weakly supervised learning, Similarity measure, Image Segmentation, Satellite datasets

Recently several models have been developed to reduce the annotation effort which is required to perform semantic segmentation. Instead of learning from pixel-level annotations, these models learn from cheaper annotations, e.g. image-level labels, scribbles or bounding boxes. However, most of these models cannot easily be adapted to new annotations e.g. new classes since it requires retraining the model. In this paper, we propose a similarity measure between pixels based on a mutual information objective to determine whether these pixels belong to the same class. The mutual information objective is learned in a fully unsupervised manner while the annotations (e.g. points or scribbles) are only used during test time. For a given image, the unlabeled pixels are classified by computing their nearest-neighbors in terms of mutual information from the set of labeled pixels. Experimental results are reported on the Potsdam dataset and Sentinel-2 data is used to provide a real world use case where a large amount of unlabeled satellite images is available but only a few pixels can be labeled. On the Potsdam dataset, our model achieves 70.22% mIoU and 87.17% accuracy outperforming the state-of-the-art weakly-supervised methods.

Introduction

Most of the successful models for semantic segmentation rely on a supervised learning approach [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]. Even though these models achieve remarkable results, the effort of collecting carefully annotated data to train these models make them impractical to use in many contexts. Generally, these models require a training dataset composed of images with pixel-level annotations, e.g. a class label is assigned to every pixel in the image. The task of labeling images is very time-consuming, e.g. the reported time to segment a single image from the PASCAL VOC 2012 dataset is around 240 seconds [START_REF] Bearman | What's the point: Semantic segmentation with point supervision[END_REF]. Consider the particular case of satellite data, many missions have been launched to observe the Earth producing massive amounts of satellite images which are absolutely impossible to be annotated by human operators. For example, each of the Sentinel-2 mission satellites [START_REF] Drusch | Sentinel-2: Esa's optical high-resolution mission for gmes operational services[END_REF] provides up to 1.6TB of images per day.

Different methods have been developed to reduce the need of carefully pixellevel annotations for large-scale data analysis. These methods propose a weakly supervised approach for semantic segmentation where the required annotations are less tedious to obtain than pixel-level annotations such as image-level annotations [START_REF] Oh | Exploiting saliency for object segmentation from image level labels[END_REF], points [START_REF] Bearman | What's the point: Semantic segmentation with point supervision[END_REF], scribbles [START_REF] Lin | Scribblesup: Scribble-supervised convolutional networks for semantic segmentation[END_REF] or bounding boxes [START_REF] Khoreva | Simple does it: Weakly supervised instance and semantic segmentation[END_REF]. These annotations are included during the training stage for learning semantic segmentation models. As a consequence, these models are not easily adaptable to new annotations (e.g. to refine the segmentation results or add new class labels) since retraining the models using these new annotations is required. For this purpose, few-shot learning techniques for semantic segmentation have been proposed [START_REF] Xian | Semantic projection network for zero-and few-label semantic segmentation[END_REF] but it still requires a significant number of labeled samples from seen classes to perform well on unseen classes. Additionally, these methods often produce suboptimal results without providing the user with an interactive way to make corrections without the need to retrain the model.

Recent work has focused on mutual information estimation and maximization to perform representation learning in an unsupervised manner [START_REF] Belghazi | Mutual information neural estimation[END_REF][START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF][START_REF] Oord | Representation learning with contrastive predictive coding[END_REF][START_REF] Ozair | Wasserstein dependency measure for representation learning[END_REF]. The main goal of these unsupervised approaches is to capture the most salient attributes of data to perform downstream tasks using the learned representations. Extensions of the previous models have been proposed using a self-supervised approach in order to capture the shared attributes from multiple views of a common context [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF][START_REF] Tian | Contrastive multiview coding[END_REF]. We think that designing a self-supervised task to learn suitable representations for semantic segmentation is an appealing idea. In particular, our work is inspired by these models [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF] to learn a similarity measure without supervision.

In this work, we take a step forward and propose a model that performs semantic segmentation by computing the similarity between pixels based on a mutual information approach without requiring annotations during training. Using an ideal similarity measure as distance metric, pixels belonging to the same class are close while simultaneously distant from pixels belonging to other classes. A very few pixel-level annotations are only used during test time. Our model computes the mutual information similarity between labeled pixels and unlabeled pixels and then performs a per-pixel nearest-neighbor search from the set of labeled pixels to classify the unlabeled pixels.

Our model provides several advantages. First, there is no need to retrain our model when new annotations are included since the similarity measure is learned using an unsupervised learning approach. Second, our model requires a small amount of annotated data which can be acquired in multiple formats e.g. points, scribbles, bounding boxes. Third, we propose a simple neural network architecture that achieves competitive semantic segmentation results while keeping a reasonable processing time.

The following contributions are made in this paper:

• We propose a model that combines a similarity measure based on mutual information between pixels using self-supervised techniques [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF] and a nearest-neighbor search to perform semantic segmentation. • We show that excellent results can be achieved by labeling less than 0.75% of the total number of pixels in an image. • We present quantitative results for image segmentation on the Potsdam dataset [START_REF] Ji | Invariant information clustering for unsupervised image classification and segmentation[END_REF] outperforming the state-of-the-art weakly-supervised methods and qualitative results on Sentinel-2 data [START_REF] Drusch | Sentinel-2: Esa's optical high-resolution mission for gmes operational services[END_REF] to show a real world use case.

• We analyze the impact of using multiple views via data augmentation techniques [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF] on the segmentation performance and we perform an ablation study to evaluate the contribution of each element of the model.

Related work

Image segmentation Exceptional results have been achieved by fully supervised models on semantic segmentation [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]. To reduce the annotation effort required by supervised learning settings, several methods have been proposed which use cheaper annotations e.g. points [START_REF] Bearman | What's the point: Semantic segmentation with point supervision[END_REF], scribbles [START_REF] Lin | Scribblesup: Scribble-supervised convolutional networks for semantic segmentation[END_REF], image annotations [START_REF] Oh | Exploiting saliency for object segmentation from image level labels[END_REF] or bounding boxes [START_REF] Khoreva | Simple does it: Weakly supervised instance and semantic segmentation[END_REF]. Labels provided by points or scribbles are then propagated to unlabeled pixels during training [START_REF] Bearman | What's the point: Semantic segmentation with point supervision[END_REF][START_REF] Lin | Scribblesup: Scribble-supervised convolutional networks for semantic segmentation[END_REF]. The main drawback is that these models are not easy to adapt to new annotations for refining the segmentation results or adding new class labels as it requires retraining the whole model. GrabCut [START_REF] Rother | grabcut" interactive foreground extraction using iterated graph cuts[END_REF] performs interactive image segmentation using a bounding box to separate foreground and background. On the other hand, Khoreva et al. [START_REF] Khoreva | Simple does it: Weakly supervised instance and semantic segmentation[END_REF] propose a semantic segmentation method requiring a costly recursive training where bounding boxes are refined iteratively. Recent work has been presented [START_REF] Xian | Semantic projection network for zero-and few-label semantic segmentation[END_REF] to segment classes containing few labels in the dataset. However, this method still requires many training examples from the known classes to perform well on the unknown classes.

Self-supervised learning In contrast to the prevalent paradigm based on generative or reconstructive models, recent work has been focused on mutual information maximization for representation learning. These models maximize the mutual information between an input and its representation. Mutual information is computed using different estimators based on the Kullback-Leibler [START_REF] Belghazi | Mutual information neural estimation[END_REF], Jensen-Shannon [START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF], Wasserstein [START_REF] Ozair | Wasserstein dependency measure for representation learning[END_REF] divergences or noise-contrastive estimation [START_REF] Oord | Representation learning with contrastive predictive coding[END_REF]. Interesting extensions of these mutual information based frameworks have been presented to capture the common attributes from paired images [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF][START_REF] Tian | Contrastive multiview coding[END_REF]. Learning representations that capture the most significant attributes of an image from multiple views is useful for semantic segmentation.

Deep metric learning Measuring the similarity between pixels is a useful tool for image segmentation since similar pixels under a given criterion belong to the same class while dissimilar pixels belong to different classes. Generally, raw pixels are mapped to a representation space by a deep neural network and then similarity between pixels is computed in the representation domain [START_REF] Chen | Blazingly fast video object segmentation with pixel-wise metric learning[END_REF][START_REF] Fathi | Semantic instance segmentation via deep metric learning[END_REF][START_REF] Sun | Neural diffusion distance for image segmentation[END_REF]. For instance, Sun et al. [START_REF] Sun | Neural diffusion distance for image segmentation[END_REF] propose a neural diffusion distance to perform segmentation. However, it requires labeled data during training to be consistent with a human criterion. For video segmentation, Chen et al. [START_REF] Chen | Blazingly fast video object segmentation with pixel-wise metric learning[END_REF] propose a metric based on the triplet loss [START_REF] Chechik | Large scale online learning of image similarity through ranking[END_REF] which is trained in a supervised manner.

In this paper, we propose a model that performs image segmentation in a weakly-supervised setting. The segmentation procedure is split into two stages. First, the model learns a mapping function from the pixel domain to a representation domain which captures relevant attributes for image segmentation using a mutual information based framework that combine the approaches [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF]. After training the mapping function, we use a mutual information objective to measure the similarity between pixels. In contrast to the models [START_REF] Chen | Blazingly fast video object segmentation with pixel-wise metric learning[END_REF][START_REF] Fathi | Semantic instance segmentation via deep metric learning[END_REF][START_REF] Sun | Neural diffusion distance for image segmentation[END_REF], the similarity measure is learned in a completely unsupervised manner. Second, our model computes the mutual information similarity between labeled pixels provided by an operator and unlabeled pixels and then performs a per-pixel nearest-neighbor search from the set of labeled pixels to propagate the labels to unlabeled pixels. The labeled pixels are only used during test time instead of training time like the models [START_REF] Bearman | What's the point: Semantic segmentation with point supervision[END_REF][START_REF] Oh | Exploiting saliency for object segmentation from image level labels[END_REF][START_REF] Khoreva | Simple does it: Weakly supervised instance and semantic segmentation[END_REF][START_REF] Lin | Scribblesup: Scribble-supervised convolutional networks for semantic segmentation[END_REF].

Background

Mutual information

The mutual information between two random variables X ∈ X and Z ∈ Z is defined in Equation 1 where p(x, z) is the joint probability density function of X and Z while p(x) and p(z) are the marginal probability density functions of X and Z, respectively.

I(X, Z) = X Z p(x, z) log p(x, z) p(x)p(z) dxdz (1)
It is straightforward to see that I(X, Z) is defined as the Kullback-Leibler divergence between the joint probability distribution P XZ and the product of the marginal distributions P X P Z , i.e. I(X, Z)=D KL (P XZ ∥ P X P Z ). Generally, computing the mutual information between high dimensional variables is a difficult task since the distributions P XZ and P X P Z are unknown. Thus, some methods based on deep neural networks have recently been proposed [START_REF] Belghazi | Mutual information neural estimation[END_REF][START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF][START_REF] Oord | Representation learning with contrastive predictive coding[END_REF][START_REF] Ozair | Wasserstein dependency measure for representation learning[END_REF].

Representation learning

Equation 1 can be used as objective for unsupervised learning where X is a variable corresponding to a given input (image, speech, text, etc) and Z is the representation of X. The representation Z is extracted by an encoder function defined by a deep neural network of parameters ψ, E ψ : X → Z, i.e. Z = E ψ (X).

The Deep InfoMax framework [START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF] proposes a mutual information estimator Î(X,Z) based on the Jensen-Shannon divergence instead of the Kullback-Leibler divergence, i.e. I (JSD) (X, Z) = D JS (P XZ ∥ P X P Z ).

Intuitively, let X i and X j be two observations of X. Let Z i and Z j be the representations of X i and X j respectively extracted via E ψ . Therefore, (X i , Z i ) is an input-representation pair sampled from the joint probability density function p(x, z) while (X i , Z j ) is an input-representation pair sampled from the product of the marginal probability density functions p(x)p(z).

We define a discriminator function defined by a deep neural network of parameters ρ, D ρ : X × Z → [0, 1] which represents the probability of a sample (X, Z) coming from p(x, z) instead of p(x)p(z), i.e. the probability that Z is the representation of X. The discriminator D ρ and the encoder E ψ are trained to assign a high probability to samples from p(x, z) (close to 1) and a low probability to samples from p(x)p(z) (close to 0) as shown in Equation 2.

max E ψ ,Dρ Î(X, Z) = E p(x,z) [log D ρ (X, Z)] + E p(x)p(z) [log (1 -D ρ (X, Z))] (2)
By redefining the discriminator function [START_REF] Nowozin | f-gan: Training generative neural samplers using variational divergence minimization[END_REF]

D ρ (X, Z) = e -T θ (X,Z) 1+e -T θ (X,Z)
where T θ : X ×Z → R is called the statistics network, we obtain the mutual information objective proposed by the Deep InfoMax framework [START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF] in Equation 3.

max E ψ ,T θ Î(X, Z)=E p(x,z) -log 1+e -T θ (X,Z) -E p(x)p(z) log 1+e T θ (X,Z) (3)
Two mutual information objectives are proposed in the Deep InfoMax framework. Maximizing the mutual information between an input X and a representation Z is called global mutual information, i.e. L global θ,ψ (X, Z) = Î(X, Z). Additionally, maximizing the mutual information between patches of the image X represented by a feature map C ψ (X) of the encoder E ψ and a feature representation Z is called local mutual information i.e. L local ϕ,ψ (X, Z)

= i Î(C (i) ψ (X), Z).

Method

In this paper, we propose a model that combines the mutual information based methods [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF] to learn a suitable representation domain to measure the similarity between pixels. Our model is trained in a fully unsupervised manner by leveraging large amounts of unlabeled data. Sanchez et al. [START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF] extends the Deep InfoMax framework to separate the common information and the exclusive information for paired images. Bachman et al. [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF] use the Deep InfoMax framework to perform self-supervised representation learning by maximizing the mutual information between representations extracted from multiple views of a shared context, e.g. the context is provided by an image and the multiple views are generated via data augmentation techniques. Learning the common information between images [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF] provides a way to compute how similar these images are. In Section 4.1, we present the mutual information objective to learn the similarity measure and we explain how to use this similarity measure to perform image segmentation in Section 4.2.

Shared mutual information

To create a suitable representation domain for image segmentation, we propose to capture the common information between images of the same context (e.g. satellite images from the same forest) into a shared representation. By removing the particular information of each image, we create a representation that distills the class information which is useful for image segmentation. We propose to learn this shared representation by using the principle presented in [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF]. Let X and Y be two images of the same context and let S X and S Y be the respective shared representations extracted by an encoder E ψ . In order to enforce learning only the common information between images X and Y , the methods [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF][START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF] maximizes the mutual information between the image X and the representation S Y and similarly, between the image Y and the representation S X . In order to create pairs of images of the same context, we follow the approach of Bachman et al. [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF] and we use data augmentation techniques (rotation, flip, pixel shift, color jitter) to create a second image from a given image, i.e. X = f (Y ) where f is a data augmentation function. We use the objective function proposed by Sanchez et al. [START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF] since it is simpler to optimize. Equations 4 and 5 displays the global and local mutual information maximization objectives.

L global M I (X, Y ) = L global θ,ψ (X, S Y ) + L global θ,ψ (Y, S X ) (4)
L local M I (X, Y ) = L local ϕ,ψ (X, S Y )+L local ϕ,ψ (Y, S X ) (5)
Sanchez et al. [START_REF] Sanchez | Learning disentangled representations via mutual information estimation[END_REF] also includes a L 1 constraint to force the shared representations to be identical as shown in Equation 6. The final objective function is displayed in Equation 7, where α, β and γ are constant coefficient.

L 1 (X, Y ) = E p(sx,sy) [|S X -S Y |] (6) max ψ,θ,ϕ L shared = αL global M I (X, Y ) + βL local M I (X, Y ) -γL 1 (X, Y ) (7)

Mutual information as similarity measure

Similarly to Chen et al. [START_REF] Chen | Blazingly fast video object segmentation with pixel-wise metric learning[END_REF], we perform per-pixel retrieval to find the closest pixel from the reference pixel set using the learned representations. A k-nearestneighbors approach is used to determine the class of unlabeled pixels by propagating the information from labeled pixels. A common way of computing the distance between pixels is to measure the L 1 or L 2 distance between their corresponding representations [START_REF] Chen | Blazingly fast video object segmentation with pixel-wise metric learning[END_REF]. Alternatively, we propose to use the global and local mutual information objectives introduced in Section 3.2.

During training, the mutual information objective is computed using an image X and a different view of X generated via data augmentation techniques, i.e. Y = f (X). In contrast, during test time the mutual information objective is computed using two different images. Let X i and X j be two image patches centered at the pixels i and j respectively and let S Xi and S Xj be the shared representations provided by the encoder function. The similarity between pixels i and j is measured by computing L global θ,ψ (X, S X ) and L local ϕ,ψ (X, S X ). After training, our model is capable to predict whether a shared representation S Xi corresponds to the image X i . Since the shared representation S Xi contains the class information of X i , it provides a means to identify pixels belonging to the same class of X i . For example, consider that X i and X j are two different images (e.g. a satellite image from an urban area and another from an agricultural area), the mutual information objective L global θ,ψ (X, S X ) (or L local ϕ,ψ (X, S X )) achieves a high score since it is easy to distinguish both images. On the other hand, suppose that X i and X j are two similar images (e.g. satellite images from the same forest), the mutual information objective L global θ,ψ (X, S X ) (or L local ϕ,ψ (X, S X )) achieves a low score since it is hard to distinguish the images.

Implementation details

Our model is composed of three deep neural networks: the encoder E ψ , the global statistics network T θ and the local statistics network T ϕ . The architecture details are provided in Figure 1. Every network is trained from scratch by using randomly initialized weights. To optimize the objective L shared defined in Equation 7, we use the Adam optimizer with a learning rate of 0.0001, β 1 =0.9 and β 2 =0.999. We use a batch size of 512 images. Images pairs are created by applying data augmentation techniques (flip, rotation, pixel shift, color jitter). Our baseline model use the following coefficients to weight the terms of the objective function L shared : α=0.5, β=1.0 and γ=0.1. The size of the shared representation is z dim =10. The training algorithm was executed on a NVIDIA Tesla K80 GPU.

The training and image segmentation procedures are summarized in Algorithms 1 and 2. More details are provided in the additional material section.

Algorithm 1 Training algorithm.

1: Random initialization of model parameters ψ (0) , θ (0) , ϕ (0) . 2: for k = 1; k = k + 1; k < number of iterations do 3: Sample a batch of C image patches {X1, . . . , XC }. Image patches have a size s.

4:

Create a new view of Xi via a data augmentation technique Yi = f (Xi). 5:

Create a batch of C paired images X : {(X1, Y1), ..., (XC , YC )}.

6:

Create a batch of C unpaired images X by shuffling the Y dimension of X.

7:

Compute L (k) = L shared (X, X, ψ (k) , θ (k) , ϕ (k) ):

L (k) = α -X sp (-T θ (C ψ (Xi),E ψ (Yi))) -X sp (T θ (C ψ (Xi),E ψ (Yi))) -X sp (-T θ (C ψ (Yi),E ψ (Xi))) -X sp (T θ (C ψ (Yi),E ψ (Xi))) + β j -X sp -T (j) ϕ (C ψ (Xi),E ψ (Yi)) -X sp T (j) ϕ (C ψ (Xi),E ψ (Yi)) -X sp -T (j) ϕ (C ψ (Yi),E ψ (Xi)) -X sp T (j) ϕ (C ψ (Yi),E ψ (Xi)) -γ X (|E ψ (Xi) -E ψ (Yi)|)
where the softplus function is defined by sp(x) = (1+e x ) 8:

Update the parameters ψ (k+1) , θ (k+1) and ϕ (k+1) by gradient ascent of L (k) . 9: end for Algorithm 2 Image segmentation algorithm.

1: Select an image Xt from the dataset. Image Xt has a size t ≫ s. 2: Label a set of L pixels into N classes P = {(p1, c1), . . . , (pL, cL)} from Xt. 3: where pi defines the coordinates and ci is the class of the i-th labeled pixel. 4: for unlabeled pixel at qj ∈ Xt do 5:

for labeled pixel at pi ∈ P do 6:

Select the image patches Xj and Xi of size s centered at qj and pi.

7:

Extract the representations SX i = E ψ (Xi) and SX j = E ψ (Xj). 8:

Extract the feature maps CX i = C ψ (Xi) and CX j = C ψ (Xj). 9:

Create the image-representation sets X = {(CX i , SX i ), (CX j , SX j )} 10:

and X = {(CX i , SX j ), (CX j , SX i )} 11:

Compute the global/local mutual information between Xj and Xi:

12: Di = L global θ,ψ = -X sp (-T θ (CX k ,SX k )) -X sp (T θ (CX k ,SX k )) or 13: Di = L local ϕ,ψ = j -X sp -T (j) ϕ (CX k ,SX k ) -X sp T (j) ϕ (CX k ,SX k ) 14:
end for 15:

Assign the pixel qj the class ci * of the nearest pixel i * = arg mini{Di} L i=1 . 16: end for 5 Experiments

Datasets

Potsdam The Potsdam dataset [START_REF]International Society for Photogrammetry and Remote Sensing: Isprs 2d semantic labeling contest[END_REF] contains 8550 aerial images of the city of Potsdam. Each image has a size of t=200 × 200 pixels and is composed of four channels: red, green, blue and infrared (RGBI). The dataset is split into three parts: 3150 unlabeled images, 4545 training labeled images and 855 test labeled images. Images are labeled into 6 classes (road, car, vegetation, tree, building and clutter). Similarly to [START_REF] Ji | Invariant information clustering for unsupervised image classification and segmentation[END_REF], we also perform image segmentation using a 3label version by merging classes (road and car, vegetation and tree and building and clutter). Image patches of size s=13 × 13 pixels are randomly sampled from the unlabeled images to optimize the model objective (Equation 7) and the test labeled images to report the experimental results. We use this dataset to provide quantitative results and comparisons to other models.

Sentinel-2 We collected 100GB of Sentinel-2 time series [START_REF] Drusch | Sentinel-2: Esa's optical high-resolution mission for gmes operational services[END_REF] by selecting several regions of interest on the Earth's surface. Images are acquired at 13 spectral bands using different spatial resolutions. We use the RGBI bands which correspond to bands at 10m spatial resolution. Our dataset is composed of 4200 time series of 12 images acquired at different dates between 2016 and 2018. The size of each image is t=512 × 512 pixels. Image patches of size s=9 × 9 pixels are randomly sampled from these images. In addition to data augmentation techniques, the function f creates an image pair by selecting an image patch Y from the same location of X but on a different date. Since there are no labels available, we use this dataset to provide qualitative results in a real world use case where a huge amount of unlabeled data is available and a few annotated pixels are provided by a human operator. Data can be downloaded from the Sentinel Hub [7]. More dataset construction details are provided in the additional material section. 1. Segmentation results. Accuracy and mIoU for N classes, M points per class and z dim =10 using the global/local mutual information in the Potsdam dataset.

Mutual Metric N = 6 N = 3 information M = 1 M = 5 M = 10 M = 50 M = 1 M = 5 M = 10 M =

Image segmentation on Potsdam

Global and local mutual information We train our model as described in Section 4.3 using the unlabeled images of the Potsdam dataset. Image segmentation is performed on test images where M pixels per class are known. Typically, these annotated pixels are provided by a human operator. To simplify the evaluation, annotated pixels are simulated by randomly sampling M pixels per class from the ground truth. We use several values of M ∈ {1, 5, 10, 50} to evaluate the performance on image segmentation. An example of the impact of M on the segmentation results is shown in Figure 2. By using the learned mutual information based similarity measure, nearest neighbor search is applied to classify pixels into one of N ∈ {3, 6} classes. The performance is reported in terms of mean intersection over union (mIoU) and accuracy. To measure the pixel similarity, we use either the global mutual information objective or the local mutual information objective (see Algorithm 2). Results are reported in Table 1. As expected, the performance is improved as M increases. Our experiments suggest that using the local mutual information objective achieves a better performance than the global mutual information when a few pixels are annotated while the performance is similar when a larger amount of annotated pixels is provided (M = 50). Many segmentation examples are shown in Figure 3.

Model comparison

To provide a comparison, we perform image segmentation using different similarity measures to search the nearest neighbor of the unlabeled pixels. First, we compute the nearest neighbor using the L 1 distance between raw pixels. Secondly, we use the L 1 distance between the representations extracted from the VAE model [START_REF] Kingma | Auto-encoding variational bayes[END_REF], Deep InfoMax model [START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF] and our model. Similar images do not necessarily have to be close in the representation domain in terms of the L 1 distance. Therefore, a low performance is expected at image segmentation using the L 1 distance between representations. Finally, we use the mutual information objective of Deep InfoMax [START_REF] Hjelm | Learning deep representations by mutual information estimation and maximization[END_REF]. As Deep InfoMax representations keeps all the image information, i.e. more than just class information, we expect this representation to be less appropriate for image segmentation. Table 2 Ablation study We analyze two important factors in our model: the influence of data augmentation techniques to generate multiple views (pixel shift, color jitter, image flip and image rotation) and the importance of some model components, e.g. the statistics networks. Results are displayed in Table 3. Several conclusions can be drawn from our experiments. First, the model architecture can be simplified since the global statistics network can be removed (α = 0) without modifying the performance on image segmentation. The local statistics network plays the most important role during training as pointed out by Bachman et al. [START_REF] Bachman | Learning representations by maximizing mutual information across views[END_REF]. Second, removing the L 1 distance between shared representations (γ = 0) leads to a slightly reduction in the performance. Third, when the shared representations are not swapped in Equations 4 and 5 (no SSR) the performance drastically decreases since these representation contains more information than the class information required for image segmentation. Concerning the data augmentation techniques, we surprisingly notice that the performance remains the same by individually removing the color jitter, image rotation and image flip. We believe that the effect of the color jitter is ignored since it is an attribute which is not captured in the shared representation. Additionally, the impact of removing the image rotation or image flip is minimal due to the local information objective where the mutual information is maximized between the representation and image patches instead of the whole image. On the other hand, removing the pixel shift degrades the performance considerably. By removing the data augmentation techniques and not swapping the shared representations (no SSR + no DA) the performance is significantly degraded. We also study the impact of the representation space dimension without noticing significant differences between z dim =10 and z dim =32.

Image segmentation on Sentinel-2 time series

Since the Sentinel-2 mission does not provide pixel-level annotations for image segmentation, we perform only qualitative experiments. In contrast to the Potsdam case where the annotated pixels are randomly sampled from the available ground truth, now we ask a human operator to label M pixels per class for each image during test time. The reader must note that scribbles, points or bounding box can be used to annotate the pixels. As these pixels are annotated under a human criterion, these pixels carry more significant information than pixels randomly sampled from the ground truth and thus the quality of the segmentation results improves significantly using just a few well-selected pixels. As shown in Figure 5 as the number of pixels per class M increases, the segmentation results considerably improve. Nevertheless, the percentage of annotated pixels remains insignificant. For instance, 60 annotated pixels in a 512 × 512 pixel image represent less than 0.03% of the total number of pixels. Also the time required for image segmentation is reasonable, an image of 512×512 pixels with 60 annotated pixels takes around 33 seconds to be segmented.

Segmentation over the time Since we maximize the mutual information between images from the same time series, the learned representation ignores the temporal information. As a consequence, by annotating pixels from a single image our model is capable to segment the whole time series the image belongs to. In Figure 6, it can be seen that the segmentation results are coherent over the time. For instance, agricultural areas are belonging to the same class regardless of whether these areas are grown or harvested.

Segmentation over the space In the same manner we perform image segmentation over the time using a single image, our model is able to do it over the space. The annotated pixels provided by a human operator are generally used to perform image segmentation on the image these pixels are extracted from. We also use these annotated pixel to segment other images from the same area achieving satisfactory results as can be seen in Figure 7. In general, using annotated pixels from a single image we can perform image segmentation on images of the same area independently of the acquisition time.

Conclusion

In this paper, we have proposed to use a mutual information based similarity measure to perform image segmentation. Our approach offers the advantage of learning the proposed similarity measure in an unsupervised manner leveraging large amounts of unlabeled data. Then, per-pixel nearest-neighbor search using the proposed similarity measure is carried out to assign classes to the unlabeled pixels from the labeled pixels provided by a human operator. In particular, we have studied the case of aerial/satellite data where massive amounts of unlabeled images are available while the annotations are scarce. In the Potsdam case, our experiments suggest that the local mutual information objective is useful to measure similarity between pixels. Our approach outperforms other approaches based on state-of-the-art methods demonstrating the usefulness of our learned representation domain. On the other hand, the ablation experiments show that the model can be further simplified as some data augmentation techniques are more relevant and the global mutual information objective can be removed. In the Sentinel-2 case, we have shown that image segmentation can be performed over the time and over the space using a very few amount of annotated pixels, e.g. labeled pixels are less than 0.002% of the total number of pixels of a time series and it can be achieved in a reasonable amount of time.

Fig. 1 .

 1 Fig.1. Network architecture. The encoder and statistics networks are implemented using convolutional and dense layers defined by the number of units, k: kernel size, f : feature maps, BN: batch normalization[START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] and activation function. The statistics networks and the encoder share weights: the input F of the statistics network is the output of the Conv 2 layer of the encoder, C ψ (X).

Fig. 2 .

 2 Fig. 2. Image segmentation examples. During test time, only M points per class randomly sampled from the ground truth are used to perform image segmentation. As M increases the accuracy and mIoU are improved. Best viewed in color and zoom-in.

Fig. 3 . 2 .

 32 Fig. 3. Image segmentation examples using M points per class randomly sampled from the ground truth. The image segmentation performance is improved as M increases.

Fig. 4 .

 4 Fig. 4. Image segmentation model comparison. Our model produce the closest predictions to the ground truth using N =6, M =50 and z dim =10 in the Potsdam dataset.

Fig. 5 .

 5 Fig. 5. Image segmentation examples in the Sentinel-2 dataset. A human operator identifies N classes in the satellite image and selects M pixels per label. a) Buenos Aires, Argentina; b) Valencia, Spain. Best viewed in color and zoom-in.

Fig. 6 .

 6 Fig. 6. Image segmentation over the time in the Sentinel-2 dataset. In the first column, the input image and the selected pixels are displayed for M =10. Our method is able to perform image segmentation with few labeled pixels on the entire time series the input image belongs to. The time series and the corresponding predictions are shown in the remaining columns for a) Toulouse, France; b) Valencia, Spain. Best viewed in color and zoom-in.

Fig. 7 .

 7 Fig. 7. Image segmentation over the space in the Sentinel-2 dataset. Selected pixels are not only useful for propagating the information from labeled pixels to unlabeled pixels in the same image but also in different images of the same area. a) Toulouse area, France; b) Tubarjal area, Saudi Arabia. Best viewed in color and zoom-in.


Dawid Machala
email: dawid.machala@isl.eu.

Floriane Collin
email: floriane.collin@univ-lorraine.fr

Marion Gilson
email: marion.gilson@univ-lorraine.fr

Marie Albisser
email: marie.albisser@isl.eu

Simona Dobre

Global Sensitivity Analysis for Modeling the Free-Flight Behavior of an Artillery Projectile

A vehicle's in-flight behavior can be represented by the Newton-Euler equations of motion: usually, such a model has a nonlinear and continuous description based on ordinary differential equations. The model structure can be altered using analytic transformations-such as transformation into a quasi-LPV structure-and, when the model is supposed to reflect an existing physical system, some of its parameters might be uncertain. This paper's objective is twofold: firstly, to investigate if a quasi-LPV model does accurately replicate the physical behavior of a large-caliber spin-stabilized projectile in free flight; and secondly, to verify if the available embedded sensors data could be used to determine non-identifiability of uncertain model parameters. In order to investigate these prospects, a permutation-based global sensitivity analysis is employed. The proposed study has highlighted the non-identifiable parameters of the model, either during the whole flight or some periods of it. These results can help to enhance the preparation of free flight experiments at the open-range test site of the French-German Research Institute of Saint Louis. Nomenclature M = Mach number V = airspeed, m s -1 α, β, α t = angle of attack, angle of sideslip, total angle of attack, deg C X , C N α , C y pα = axial force, normal force slope, and Magnus force slope coefficients C lp , C mα , C mq , C npα = roll damping, pitch moment slope, pitch damping and Magnus moment slope coefficients

I x , I t = longitudinal and lateral moments of inertia, kg m 2 u, v, w = linear velocities in x, y, and z body axis, m s -1 p, q, r = roll, pitch, and yaw rate, rad s -1

x E , y E , z E = center of gravity position, expressed with respect to the Earth reference frame, m φ, θ, ψ = Euler angles: roll, pitch, and yaw, deg A i , G i , H i = accelerometer, gyroscope, and magnetometer signal in i-th sensor axis

I. Introduction

A important trend in the development of 155 mm artillery ammunition aims towards higher precision and longer range, especially through the use of various guidance and navigation techniques to assure their capability of hitting even a small or moving target. Such a development relies on a well-investigated and -understood mathematical description of a projectile's behavior. A first-principles' model of a projectile is described by nonlinear Newton-Euler laws of motion, and depends on a set of aerodynamic coefficients that need to be quantified for each flying vehicle [START_REF] Zipfel | Modeling and Simulation of Aerospace Vehicle Dynamics[END_REF]. In other words, the knowledge of the aerodynamic characterisation allows to predict the trajectory and assess the dynamical behavior of a projectile, hence underlying the importance of their quantification.

The aerodynamic coefficients can be quantified using several complementary methods, such as Computational Fluid Dynamics simulations (e.g. in [START_REF] Hamel | CFD and Parametric Study on a 155 mm Artillery Shell Equipped with a Roll-Decoupled Course Correction Fuze[END_REF]), semi-empirical codes (e.g. PRODAS [START_REF] Bryson | Approach for Understanding Range Extension of Gliding Indirect Fire Munitions[END_REF] and Missile DATCOM [START_REF] Hamel | CFD and Parametric Study on a 155 mm Artillery Shell Equipped with a Roll-Decoupled Course Correction Fuze[END_REF]), or wind-tunnel experiments (e.g. in [START_REF] Strub | Modeling, Identification and Control of a Guided Projectile in a Wind Tunnel[END_REF]). Their results then serve as a first guess of the aerodynamic coefficients, later to be confirmed under a real-world experimentation-ballistic free-flight tests [START_REF] Albisser | Aerodynamic Coefficient Identification of a Space Vehicle from Multiple Free-Flight Tests[END_REF].

Free-flight tests pose several technical and scientific difficulties. It must be assured that the projectile will land in a dedicated safe zone, even though some of the initial state variables are unknown; or that the identified model structure will be useful for the guidance and navigation tasks, even though the physical model is nonlinear. Moreover, due to the fact that the experiments are performed in real-flight conditions, several sources of uncertainty are introduced to model parameters. For instance, initial state vector will differ from its nominal value due to stochastic effects induced by interior ballistics [START_REF] Carlucci | Ballistics: Theory and Design of Guns and Ammunition[END_REF]. Similarly, an initial guess of the aerodynamic coefficients, based on similar architectures in the literature, will likely be imprecise. These uncertainties influence the model dynamics, and thus propagate through the model into its output signals.

The influence that the uncertain model parameters have on the model output can be assessed by the means of sensitivity analysis that apportions the output uncertainty to its different sources [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF]. Such a knowledge offers the possibility of evaluating the influence of parameter uncertainties on the model behavior [START_REF] Saltelli | Sensitivity analysis in practice[END_REF][START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF][START_REF] Cukier | Nonlinear Sensitivity Analysis of Multiparameter Model Systems[END_REF]. When assessed, the 'least' influential parameters can be set to a nominal value, since they do not significantly affect the model dynamics, whereas the identification and modelling efforts can be concentrated on the most influential parameters. The methods of sensitivity analysis can be divided into local and global approaches. The local technique assess the uncertain parameters' influence around nominal operating points [START_REF] Saltelli | Sensitivity analysis in practice[END_REF], primarily using one-at-a-time sampling and differential analysis [START_REF] Hamby | A review of techniques for parameter sensitivity analysis of environmental models[END_REF].

However, they fail to characterize the behavior of the model far from the selected nominal condition, and-in the case of one-at-the-time sampling-explore only a limited area of the state-space [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF]. The global methods aim to explore the full space of uncertain model parameters. There are several strategies to assess the global sensitivity of the model.

Principally the variance-based analysis using so-called Sobol' indices is employed (also called Analysis of Variance, ANOVA) [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF][START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF], which allows to determine both the individual influence of each parameter and their interactions [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF][START_REF] De Rocquigny | Chapter 18. Sensitivity analysis methods[END_REF].

In aerospace applications, sensitivity analysis is usually conducted using local sensitivity analysis, as in [START_REF] Torres | Trajectory Accuracy Sensitivity to Modeling Factors[END_REF][START_REF] Wang | Sensitivity Analysis of Multidisciplinary Rotorcraft Simulations[END_REF][START_REF] York | Efficient Aircraft Multidisciplinary Design Optimization and Sensitivity Analysis via Signomial Programming[END_REF][START_REF] Dannenhoffer | Design Sensitivity Calculations Directly on CAD-based Geometry[END_REF][START_REF] Makhija | Spiral: A General Framework For Parameter Sensitivity Analysis[END_REF][START_REF] Albisser | Sensitivity Analysis for Global Parameter Identification. Application to Aerodynamic Coefficients[END_REF].

Only a few works deal with global sensitivity analysis (GSA). Among them, the prominent applications include the works of Huan et al. [START_REF] Huan | Global Sensitivity Analysis and Estimation of Model Error, Toward Uncertainty Quantification in Scramjet Computations[END_REF] who used a surrogate modeling approach to identify the most influential parameters in scramjet combustor simulations. Wei et al. [START_REF] Wei | Extended Monte Carlo Simulation for Parametric Global Sensitivity Analysis and Optimization[END_REF] employed a modified Monte-Carlo approach to calculate sensitivity indices in order to minimize the variance of model output for a wing-structure model. Jiang et al. [START_REF] Jiang | Multidisciplinary Statistical Sensitivity Analysis Considering Both Aleatory and Epistemic Uncertainties[END_REF] have used GSA in order to guide a design process of an aircraft by assessing which model parameters influence the most its final acquisition cost or maximum flight time. It is worth noting that the results of a local sensitivity analysis can be misleading or even significantly different than of the global study. For example, Hale et al. [START_REF] Hale | Sensitivity Analysis Methods for Systems with Epistemic Uncertainties[END_REF] have shown that local sensitivity analysis might wrongly indicate some of the model parameters as influential; as well as it significantly underestimates the influence of aerodynamic coefficients due to their nonlinearity. However, to authors' knowledge, only local sensitivity analysis has been applied in the existing studies of free-flight experiments in the domain of exterior ballistics.

In the context of the aforementioned experimental and practical difficulties associated with free-flight experiments, the authors investigate the possibilities of enhancing the traditional procedure for conducting the tests by the use of GSA. Potentially, this technique could be used as an effective simulation-based tool of visualizing physical knowledge in ballistic systems and detecting non-identifiable model parameters. Moreover, this technique could be seen as an simulation-based alternative to observability analysis and information theory, that are sometimes used in exterior ballistics to detect which model parameters should be identified, as shown, e.g. in [START_REF] Rogers | Effective use of magnetometer feedback for smart projectile applications[END_REF][START_REF] Rogers | A low-cost orientation estimator for smart projectiles using magnetometers and thermopiles[END_REF][START_REF] Sahu | Aeromechanics and Control of Projectile Roll Using Coupled Simulation Techniques[END_REF].

The following provides the problem statement: to attribute the sources of uncertainties on the simulated embedded sensor signals to the uncertain model parameters using GSA. For this purpose, a first-principles-based model has been transformed into a quasi-linear parameter varying (LPV) model of lower computational complexity. Such a model structure is selected due to its prevalence in the domain of control and navigation. GSA is then performed on the resulting quasi-LPV model. It is then checked if the implications of the sensitivity indices on the identifiability of uncertain parameters are coherent with the domain knowledge of exterior ballistics. The quasi-LPV model structure has been recently presented [START_REF] Machala | Quasi-LPV modelling of a projectile's behaviour in flight[END_REF] and its accuracy has been numerically confirmed [START_REF] Machala | The Analysis of Vehicle's In-Flight Behaviour Using Quasi-LPV and Nonlinear Models[END_REF].

The structure of the paper is as follows. In section II, the nonlinear and quasi-LPV model structures are introduced, and the sources of their uncertainty are explained. Section III provides an overview of the theory behind the GSA whereas its results when applied to the quasi-LPV model are summarized in section IV. Section V concludes the article by summarizing the main contributions and future prospects of this work.

II. Modeling

The nonlinear model, based on the 6-degrees-of-freedom Newton-Euler laws of motion [START_REF] Albisser | Identification of aerodynamic coefficients from free flight data[END_REF], and its representation in a quasi-LPV form [START_REF] Machala | Quasi-LPV modelling of a projectile's behaviour in flight[END_REF], are presented hereafter. Both models assume that the mass distribution is uniform and that the shape of projectile is axis-symmetrical.

A. Nonlinear model

The model of the in-flight dynamics is a nonlinear state-space continuous-time structure M:

M :            x(t) = f (x(t), C(x(t), p a )), x(0) = x 0 y(t) = g(x(t)), (1)
where the state vector x ∈ R 12 has 12 components corresponding to the translational velocities (u, v, w), rotational velocities (p, q, r), position of the center of gravity of the vehicle with respect to Earth (x E , y E , z E ), and the roll, pitch and yaw angles (φ, θ and ψ, respectively) [START_REF] Zipfel | Modeling and Simulation of Aerospace Vehicle Dynamics[END_REF]. The aerodynamic coefficients C are described by polynomial functions whose time invariant parameters are p a . These descriptions depend on the linear and angular velocities, as well as on the Mach number M and the total angle of attack α t , which are defined by:

M = V a , α t = arccos u V , ( 2
)
where a is the speed of sound, and V is the total velocity, i.e. V = √ u 2 + v 2 + w 2 [START_REF] Albisser | Identification of aerodynamic coefficients from free flight data[END_REF]. The output vector y represents the data obtained by on-board and on-ground measurement techniques available during free flight tests. Herein, embedded magnetometers, accelerometers and gyroscopes are considered.

The full form of the nonlinear model equations in the variable-roll reference frame is described by the following set of equations [START_REF] Albisser | Identification of aerodynamic coefficients from free flight data[END_REF]:

Translational dynamics             u v w             =             0 -r q r 0 -p -q p 0                         u v w             + g             -sθ sφcθ cφcθ             + 1/m             f x,aero f y,aero f z,aero             , (3)
Attitude dynamics

            p q r             =             0 p • r I t -I x I t p • q I x -I z I t             +             1 I x m x,aero 1
I y m y,aero 1
I z m z,aero             , ( 4
)
Translational kinematics             x E y E z E                         cθcψ sφsθcψ -sψcφ cφsθcψ + sφsψ cθsψ sφsθsψ + cψcφ cφsθsψ -sφcψ -sθ sφcθ cφcθ                         u v w             , (5)
Attitude kinematics

            φ θ ψ             =             1 sφtθ cφtθ 0 cφ -sφ 0 sφ 1 cθ cφ 1 cθ                         p q r             , (6)
where the forces and moments are given by:

            f x,aer o f y,aer o f z,aer o             = qS             C X 0 + C X a s 2 α t -C ypα d 2V w V p -C N α v V -C ypα d 2V v V p -C N α w V             , (7)
            m x,aer o m y,aer o m z,aer o             = qSd             C lp p C mα w V -C npα d 2V v V p + C mq d 2V q -C mα v V -C npα d 2V w V p + C mq d 2V r             . ( 8
)
The state equations depend on various physical properties, measured or determined for each projectile before the flight tests, such as the mass m, longitudinal and lateral moments of inertia I x , I t , reference surface S, and diameter d. The dynamic pressure is denoted by q = 1 2 ρV 2 , where ρ is the air density, also assumed to be known. Symbols s, c, t correspond to trigonometric sine, cosine, and tangent, respectively.

The aerodynamic coefficients can be divided into the force and moment coefficients. The forces are described via the axial force C X (decomposed as

C X = C X 0 + C X a s 2 α t ),

B. Quasi-LPV model

The proposed quasi-LPV model is based on analytic function substitution applied to the nonlinear model equations.

In order to simplify the computational complexity of the nonlinear model, and to aid the analytic substitution process, a coordinate-system transformation has been applied to the nonlinear model equations.

In aerospace modeling, the usual choice for a coordinate system location is to rigidly fix it on the vehicle's center of mass and to subject it to all the projectile's motions [START_REF] Zipfel | Modeling and Simulation of Aerospace Vehicle Dynamics[END_REF][START_REF] Wernert | Flight Dynamics Properties of 155 mm Spin-Stabilized Projectiles Analyzed in Different Body Frames[END_REF]. Alternatively, the model can be developed with respect to a frame that is fixed to the projectile's center of mass, but does not experience the rolling motion [START_REF] Wernert | Flight Dynamics Properties of 155 mm Spin-Stabilized Projectiles Analyzed in Different Body Frames[END_REF]. Such descriptions are called "variable-roll" and "non-rolling" reference frames, respectively. The latter formulation has an interesting property if the projectile is axis-symmetric (in y and z axes) and spin-stabilized: the nonlinear model equations become decoupled from its fastest changing variable, the roll angle φ [START_REF] Wernert | Flight Dynamics Properties of 155 mm Spin-Stabilized Projectiles Analyzed in Different Body Frames[END_REF]. As a result, the computational complexity diminishes.

The non-rolling frame has been applied for this purpose as early as in 1960s (e.g. in [START_REF] Barnett | Trajectory equations for a six-degree-of-freedom missile[END_REF]).

Transformation of dynamics state variables in variable-roll reference frame (3)-( 4) into their non-rolling equivalents can be described using the following transformation equations as:

            u v w             =             1 0 0 0 cφ -sφ 0 sφ cφ                         u v w             ,             p q r             =             1 0 0 0 cφ -sφ 0 sφ cφ                         p q r             , (9)
where the superscript • indicates a variable expressed with respect to the non-rolling frame. The kinematic state variables, i.e. Euler angles and positions, are equivalent in both frames.

A quasi-LPV model can be obtained in a non-rolling frame through a set of algebraic transformations, in accordance with the function substitution approach [START_REF] Pfifer | Quasi-LPV Model of a NDI-Controlled Missile Based on Function Substitution[END_REF]. The full derivation of the model has been presented in [START_REF] Machala | Quasi-LPV modelling of a projectile's behaviour in flight[END_REF], and results in the following set of equations:

x = A(x ) • x , (10)
where the non-rolling state vector x = u v w p q r θ ψ x E y E z E T ∈ R 11 , and the parameterdependent state matrix can be defined as follows:

A(x) =                                             0 -r q X 1 (u ,V,z E )-g(z E )sθ p 0 0 0 0 r X 3 (V, z E ) r tθ + X 2 (V, z E )p 0 0 0 0 0 -q -r tθ -X 2 (V, z E )p X 3 (V, z E ) g(z E )cθ p 0 0 0 0 0 0 0 X 4 (V, z E ) 0 0 0 0 0 X 7 (V, z E )p X 5 (V, z E ) -I x I s r -r r tθ p X 6 (V, z E ) 0 0 0 0 -X 5 (V, z E ) X 7 (V, z E )p I x I s q + q r tθ p 0 X 6 (V, z E ) 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 c -1 θ 0 0 cθcψ -sψ sθcψ 0 0 0 0 0 cθsψ cψ sθsψ 0 0 0 0 0 -sθ 0 cθ 0 0 0 0 0                                             , (11)
where:

X 1 (u, V, z E ) = -0.5m -1 ρ(z E )V 2 SC X (u, V) X 2 (V, z E ) = 0.25m -1 ρ(z E )SdC ypα (V) X 3 (V, z E ) = -0.5m -1 ρ(z E )V SC N α (V) X 4 (V, z E ) = 0.25ρ(z E )V Sd 2 C lp (V) X 5 (V, z E ) = 0.5ρ(z E )V SdC mα (V) X 6 (V, z E ) = 0.25ρ(z E )V Sd 2 C mq (V) X 7 (V, z E ) = 0.25ρ(z E )Sd 2 C npa (V).
The scheduling signals (i.e. the time-varying elements of the state matrix) of the quasi-LPV model consist of eight variables, which means that the model could be considered as a rather complex in control application-especially that the complexity of an LPV controller increases exponentially with the number of variables [START_REF] Kwiatkowski | Parameter Reduction for LPV Systems via Principal Components Analysis[END_REF]. However, such a structure can be deemed to be a good compromise between control and identification requirements. Firstly, should the complexity be reduced, fixing some of the scheduling signals will result in a control-oriented structure presented in [START_REF] Sève | Pitch/Yaw Channels Control Design for a 155mm Projectile with Rotating Canards, using a H-infinity Loop-Shaping Design Procedure[END_REF].

Secondly, should the model behavior be accurate for a wide variation of flight conditions, which is necessary for the identification step [START_REF] Albisser | Identification of aerodynamic coefficients from free flight data[END_REF], the model can be left 'as-is'. Out of the eight scheduling signals, five of them can be measured by available measurement techniques: the angular rates by gyroscopes, the altitude by a trajectory tracker, and the total velocity by a Doppler radar. The remaining scheduling signals, i.e. u , θ, ψ, can be simulated or estimated, as in [START_REF] Changey | Attitude Estimation of a projectile using Magnetometers and Accelerometers: Experimental Validation[END_REF].

For the purpose of conducting GSA, the case of a projectile equipped with a set of embedded sensors will be considered in flight simulations. The analysis will be, for the sake of simplicity, shown only for the embedded sensors.

Therefore, the following data are analyzed: the signals of axial A x and radial A y , A z accelerometer; radial G y , G z gyroscopes; and axial H x and radial H y , H z magnetometers. These sensors have been previously used and calibrated in ISL's works, and their equations with further details are presented in [START_REF] Albisser | Identification of aerodynamic coefficients from free flight data[END_REF][START_REF] Albisser | Aerodynamic Characterization of a New Concept of Long Range Projectiles from Free Flight Data[END_REF][START_REF] Machala | Free-Flight Projectile Behaviour: LPV modelling and global sensitivity analysis[END_REF][START_REF] Adam | Development of a gun-hardened IMU based on low-cost sensors: Hardware design and modeling[END_REF]. All of their signals are time-varying.

C. Sources of model uncertainty

When a given model structure is simulated, its trajectory depends on the initial state vector x 0 and on the set of model parameters, such as its aerodynamic coefficient parameters p a , mass, and inertia. In the following subsections, the set of nominal values of x 0 and p a will be discussed, together with the possible sources of uncertainty influencing them.

Uncertainties of initial state variables

State variables in the model are divided into these describing translational dynamics, attitude dynamics, translational kinematics, and attitude kinematics. Specifics of each influence the model uncertainty in a different way:

• Translational dynamics. In the case of free-flight experiments, as performed at the ISL's open-range test site (see [START_REF] Albisser | Identification of aerodynamic coefficients from free flight data[END_REF]), the flights are usually performed over a range of flight conditions (e.g. sub-, trans-, and super-sonic regimes) in order to sufficiently explore the projectile flight dynamics. Herein, a selection of four different initial velocities u 0 , ranging from 300 to 600 m s -1 has been chosen as the set of nominal values. Since the projectiles are launched from a rifled-bore powder gun, and each firing varies slightly depending on inaccuracies in powder amount as well as due to internal ballistic effects (such as differences in firing temperature, pressure distribution in-bore, etc [START_REF] Carlucci | Ballistics: Theory and Design of Guns and Ammunition[END_REF]), an uncertainty range of ±1% is assumed. The stochastic effects of internal ballistics on the remaining linear velocities are negligible, therefore v 0 , w 0 can be assumed to be zero, without any interval of uncertainty imposed upon them.

• Attitude dynamics. In the case of the initial angular rates, the effects of internal ballistics are non-negligible. The initial roll rate p 0 is a function of initial axial velocity u 0 and the twist rate of the launcher's bore [START_REF] Carlucci | Ballistics: Theory and Design of Guns and Ammunition[END_REF]. Herein, p 0 = 2.0283 • u 0 has been assumed, which renders p 0 a linear transformation of u 0 . Since u 0 is already considered as an uncertain parameter, no uncertainty is assumed on p 0 . For the initial yaw and pitch rates, the nominal value is q 0 = r 0 = 0 rad s -1 . However, their nominal values can be influenced by the stochastic uncertainties stemming from the internal ballistics [START_REF] Carlucci | Ballistics: Theory and Design of Guns and Ammunition[END_REF]. Due to a lack of more precise estimate, a significant interval of uncertainty of ±3 rad s -1 has been assumed on q 0 and r 0 .

• Attitude kinematics. The initial position of a projectile is equivalent to the position of the launcher. In this case study, it corresponds to x E 0 = y E 0 = 0 m with the initial altitude of z E 0 = -1.5 m, corresponding to the launcher height. Finally, concerning the initial attitude of the projectile, it is dependent on the target location, as the launcher aims to deliver the projectile to its terminal destination, which in the case of ISL's test site corresponds to the centre of a sandbay, located at [1000, 0, -6] m. The negative value of the altitude stems from the definition of the Earth reference frame: positive altitude axis aims towards the center of Earth.

• Translational kinematics. The initial attitude of the projectile has to be determined, as it has a direct influence of the trajectory and point of impact. In order to find the initial pitch θ 0 and yaw ψ 0 angles allowing to hit the target, a least squares estimation has been performed, using a set of fixed nominal values of the remaining state variables, i.e. only θ 0 , ψ 0 were estimated. The resulting values have been selected as the new nominal values for GSA.

An uncertainty interval of ±0.00049 rad has been assumed on θ 0 and ψ 0 that corresponds to the measurement precision of the launcher's elevation and azimuth at the ISL's test range. The initial roll angle φ 0 can be fixed to zero, as it does not influence the trajectory.

All the nominal values of state variables are summarized in Table 1, together with their imposed uncertainty intervals.

Uncertainties of model parameters

In the context of free-flight experiments performed at ISL, the in-flight behavior of a projectile is investigated for a flat-fire scenario. Under such circumstances, the axial force coefficient C X and normal force coefficient C N α have a high influence on projectile dynamics [START_REF] Carlucci | Ballistics: Theory and Design of Guns and Ammunition[END_REF]. The pitching moment is considered as well-known for the 155mm architecture.

Consequently, for this study only C X and C N α have been chosen as the objects of interest for GSA. The nominal values of all aerodynamic coefficients are known in a form of lookup tables, for which an initial guess is based on previous free-flight experiments, results from the literature, semi-empirical codes, CFD simulations, etc. Each of the methods offers a different way to quantify the projectile aerodynamics, hence it can also be considered as a source of uncertainty over the initial guess of the values of the aerodynamic coefficients. Herein, an additional 5% of uncertainty is considered.

Secondly, an uncertainty of the projectile's mass m is imposed. The choice of mass as a parameter of interest is on one hand realistic, since some mass uncertainty is to be expected [START_REF] Carlucci | Ballistics: Theory and Design of Guns and Ammunition[END_REF]; and on the other hand, mass uncertainty has only a negligible influence on the trajectory over a short-range firing. Therefore, it can be considered as an example of parameter which would be difficult to identify using embedded sensors.

All the uncertainty intervals of the uncertain model parameters are summarized in Table 2.

Table 2 Uncertain projectile parameters

Parameter Nominal value Uncertainty Unit

m 45 ±0.25% kg C X tabular ±5% - C N α tabular ±5% -
An underlying assumption permitting the use of variance-based techniques lays in the modality of the uncertainty distribution. Shall the distribution be, e.g. multi-modal, other techniques-such as entropy-based approach-would be preferred [START_REF] Auder | Safety, reliability and risk analysis : theory, methods and applications[END_REF]. An another assumption is the one of a uniform distribution of the uncertainty of model parameters. Even though such an assumption might in some cases change the results of sensitivity analysis (e.g. in [START_REF] Paleari | Sensitivity analysis of a sensitivity analysis: We are likely overlooking the impact of distributional assumptions[END_REF]), it is chosen here both for the sake of simplicity and due to the lack of parameter values around which the distribution could be easily centered, especially for the case of aerodynamic coefficients. Should such a value exist, normal or triangular distributions may become a better choice. In the present case, there is no a priori reason for any parameter value to be more probable than another one in the uncertainty interval.

III. Global Sensitivity Analysis theory

Sensitivity analysis answers the question of how can the uncertainty of the model output be apportioned to different sources of uncertainty in the model input [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF]. Such a knowledge offers the possibility of assessing the influence of parameter uncertainties on the model behavior, which can be applied for several purposes [START_REF] Saltelli | Sensitivity analysis in practice[END_REF][START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF][START_REF] Cukier | Nonlinear Sensitivity Analysis of Multiparameter Model Systems[END_REF], two of which are of interest for the present article:

• Test the understanding of the model, e.g. by testing the coherency between the results and the a priori knowledge in the field, in order to confirm the physical behavior of the simulated model.

• Provide insight into identifiability of the model's parameters, i.e. taking advantage of the knowledge that non-influential parameters are also non-identifiable [START_REF] Dobre | Limits of variance-based sensitivity analysis for non-identifiability testing in high dimensional dynamic models[END_REF]. Therefore, GSA might be used to answer the question if the uncertainty in sensor signals can be apportioned to any of the uncertain parameters -if not, they would be deemed non-identifiable.

In the context of GSA, the general model structure from Equation (1) has to be extended to allow for stochastic variables to be present in the model. Therefore, we will consider a modified structure M GS A , such that:

M GS A :            x(t, χ) = f (x(t, χ), C(x(t, χ), p a ( χ))) , x(0, χ) = x 0 y(t, χ) = g (x(t, χ)) , (12)
where

y(t) = [A x , A y , A z , G y , G z , H x , H y , H z ] T .
The model [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] depends on the initial state vector x(0, χ) and on the vector of parameters of aerodynamic coefficients p a ( χ). These two vectors contain elements that are uncertain and can be considered as random variables, dependent on an elementary random event χ. The vector of their uncertain elements is denoted as

Ω ∈ R s , Ω = [Ω 1 ( χ), Ω 2 ( χ), ..., Ω s ( χ)],
where Ω i ( χ) ∈ R. For the sake of simplicity, the notation ( χ)

will be dropped hereinafter.

A. Sensitivity indices

Sensitivity indices allow to quantify how much, on average, the output variance would be reduced if an uncertain parameter Ω i ∈ Ω was fixed to its constant nominal value. In other words, they allow to assess the effect that the uncertain parameter Ω i has on a j-th output of the model, noted as y j . Main sensitivity indices are defined for any given time instant, following [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF][START_REF] Anstett-Collin | Sensitivity analysis of complex models: Coping with dynamic and static inputs[END_REF], as:

S j i = V(E(y j |Ω i )) V(y j ) . ( 13
)
where V(y j ) is the total variance of y j and V(E(y j |Ω i )) is the conditional variance of y j given Ω i .

It is assumed that uncertain parameters Ω i are independent. The value of S i lies between 0 and 1. The closer to 1, the higher the contribution of Ω i to the total variance of the output. The sum of main sensitivity indices offers the following insight into the influence of uncertain parameters on the model, for a given output signal:

• i S i = 1 if the model is additive, indicating that the main sensitivity indices are sufficient to fully describe the model uncertainties.

• i S i < 1 otherwise, indicating the influence of interactions between model parameters [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF].

In the latter case, total sensitivity indices can be used to assess both the individual influence of a parameter Ω i , and its interactions with other parameters [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF]. The total sensitivity indices are defined as:

S j T i = E(V(y j |Ω ∼i )) V(y j ) . ( 14
)
where notation Ω ∼i indicates all the elements of Ω except for Ω i .

Similarly as in the case of the main sensitivity indices, the sum of all total indices for a given output signal:

• i S T i = 1 if the model is additive, • i S T i > 1 otherwise.
Overall, all the important rules for interpreting the sensitivity indices, with respect to a selected output signal, are summarized in Table 3.

Table 3 Rules for interpreting the obtained sensitivity indices.

Observation Interpretation Reference

S i ≈ 1

The parameter x i influences directly the output uncertainty.

[7] S i ≈ 0

The parameter x i does not influence directly the output uncertainty. [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF] s i S i ≈ 1

The model is additive, interactions are not present, i.e.

S T i = S i [43] S T i ≈ 1
The parameter x i influences (directly or via interactions) the output uncertainty.

[7]

S T i ≈ 0

The parameter x i does not influence the output uncertainty and is non-identifiable.

[7, 45]

B. Estimation of the sensitivity indices

Since the original Sobol' work from 1993, several methods have been developed for computing the variance-based sensitivity indices in more computationally efficient way, such as: the Fourier amplitude sensitivity test [START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF][START_REF] Tarantola | Variance-Based Sensitivity Indices of Computer Models With Dependent Inputs: the Fourier Amplitude Sensitivity Test[END_REF],

the surrogate models approach [START_REF] Ratto | Non-parametric estimation of conditional moments for sensitivity analysis[END_REF][START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF], and Monte Carlo simulation improved with various techniques to reduce its computational complexity [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF][START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF][START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF]. Therefore, the choice of sensitivity indices estimation method is not a trivial task.

In this paper, the approach of applying GSA, as opposed to the local analysis, is proposed. For this purpose, the Monte Carlo simulation-based method is suggested, as it allows a direct computation of sensitivity indices, and for its simplicity in implementation. The authors take advantage of the non-rolling frame transformation in order to reduce the simulation time of each model run.

In order to reduce the required number of simulations, a permutation-based approach proposed by Mara [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF] is employed. This approach consists of three steps:

1) A first model response vector

Y 1 = g(Ω 1
) is generated, based on the measurements function g of the first Sobol'

sample Ω 1 . The sample matrix has a size N × s, where N is the amount of Monte-Carlo samples, and s is the number of uncertain model parameters.

2) A second sample matrix Ω 2 is generated based on column-wise permutation of Ω 1 matrix: each i-th column (out of d columns in total) is permuted by a permutation vector r i . Hence, noting the permutation matrix as R = [r 1 , r 2 , ..., r s ], the second sample matrix is obtained as

Ω 2 = Ω 1 • R. The second model response vector is calculated as Y 2 = g(Ω 2 ).
3) Main sensitivity indices are calculated for each of j output signal with respect to each of s uncertain parameters, as:

Ŝj i = (N -1) -1 N k=1 (y r i k j -μj 1 )(y 2 k j -μj 2 ) σ j 1 σ j 2 , (15)
where y r i k j is k-th element of the j-th output signal obtained through reordering the y 1 k j accordingly with column-wise permutation r i . Symbols µ and σ denote the mean value and the standard deviation, respectively, and . indicates an estimated value.

For the purpose of upcoming GSA, eight parameters from Tables 1 and2 have been deemed as uncertain:

Ω = [u 0 , q 0 , r 0 , θ 0 , ψ 0 , m, C X , C N α ] T . Such a setting of uncertain parameters allows to evaluate their influence on the system dynamics, which in turn influences the sensors' measurements, as noted in Fig. 1. Therefore, GSA can be used to answer which parameters, and to what extent, do influence the measurement sensor equations. Moreover, the number of uncertain parameters is equal to s = 8, whereas each simulation of a single trajectory of the quasi-LPV model lasts less than three seconds [START_REF] Machala | The Analysis of Vehicle's In-Flight Behaviour Using Quasi-LPV and Nonlinear Models[END_REF]. Such a case further supports the selection of a variance-based methods as the basis of the global sensitivity assessment: for the cases when the number of parameters is lower than 20, and the CPU time per one run is lower than one minute, the choice of variance-based methods is suggested as a rule-of-a-thumb by De Rocquingy [START_REF] De Rocquigny | Chapter 18. Sensitivity analysis methods[END_REF], as well as by Saltelli et al. [START_REF] Saltelli | Sensitivity analysis in practice[END_REF].

IV. Simulation results

Simulations for the purpose of GSA using Mara's method have been performed based on the set of nominal conditions and their uncertainty intervals as summarized in Tables 1 and2. The analysis has been performed in a time-varying approach-as information on sensitivity could be hidden at some periods of time, while being readily visible in the others [START_REF] Pianosi | Understanding the time-varying importance of different uncertainty sources in hydrological modelling using global sensitivity analysis[END_REF]. The value of sensitivity indices is calculated every 2 ms, which is faster than the quasi-LPV The procedure of obtaining the results can be summarized as follows:

• Perform a sampling of the uncertain parameter space, using the uncertainty ranges as dictated by Tables 1 and2, while assuming that the samples are uniformly distributed. The sampling size, for which convergent results have been obtained, equals to two samples Ω 1 and Ω 2 , each of the same size, with the total sample size of N = 40000.

• Simulate the model trajectory (i.e. the values of the time-varying output signals) for each of the N samples over the distance of 1 km corresponding to the distance between the projectile launcher and the sandbay.

• Calculate the sensitivity indices S i , i = 1, 2, ..., s, based on Equation [START_REF] Wang | Sensitivity Analysis of Multidisciplinary Rotorcraft Simulations[END_REF] for every time instant.

A. Results interpretation for the full parameter vector

All the main sensitivity indices, calculated using Mara's method, have been plotted for each time instant (see Figure 2 and3).

Firstly, the initial translational and angular velocities are dominating the uncertainty of the output signals, as shown in Fig. 2 (a) for the example of the radial accelerometer A y . This result, although in-line with physical expectations, obscures the influence of other parameters on the output signals. That is especially visible for the case of the aerodynamic coefficients. For instance, as shown in Fig. 2 (b), it appears as if uncertainty around the nominal value of C X does not visibly influence the measurement signals other than the axial accelerometer A x . Such an effect would be worrying, as it could indicate non-identifiability of the aerodynamic coefficients by most of the sensors. However, this effect stems from the exponentially higher influence of initial velocities compared to other factors, rather than from the suspected lack of influence of the aerodynamic coefficients. Since u 0 , q 0 , r 0 almost totally dominate, only a little effect is visible from the other parameters as observed through sensitivity indices normalised by definition to a [0,1] range. Sensitivity indices are relative to each other, and hence neither can they be directly compared between outputs, nor will a clear influence of the uncertain parameter be visible if one of them logarithmically dominates. The uncertain translational and angular velocities influence the sensors uncertainty to such an extent, that the desired effects of Euler angles and aerodynamic coefficients uncertainty is invisible. The oscillations visible in in Fig. 2 (a) stem from the rotation of a spin-stabilised projectile, and change over time together with changing precession/nutation patterns typical for spin-stabilised projectiles [START_REF] Wernert | Flight Dynamics Properties of 155 mm Spin-Stabilized Projectiles Analyzed in Different Body Frames[END_REF]. Another explanation would be to conclude that the influence of uncertain parameters is expressed through their interactions, and thus the total sensitivity indices would have to be calculated. However, such an assumption can be easily disproved, since s i=1 S j i > 0.9 for almost all of the outputs j for the whole duration of the flight, which in accordance with Table 3 implies that the model is approximately additive (Fig. 3). Thus, the main indices S In order to be able to assess the influence of the uncertainty of the aerodynamic coefficients, as well as of the other less influential model parameters, it has been decided to repeat the analysis for a modified set of uncertain parameters, as explained below.

0 0.5 1 1.5 2 Time [s] 0 0.1 0.2 0

B. Results interpretation for the reduced parameter vector

In the sequel, a modified set of parameters is examined. This time, the initial velocities are considered only at their nominal values, hence Ω = [θ 0 , ψ 0 , m, C X , C N α ] T , (s = 5).

The results observed in such a setting, as shown in Fig. 4 to 6 are more promising, as it is now visible that they • The mass influence on sensor measurements are negligible. This effect is visible in Fig. 6 (a), since all sensitivities remain close to zero. Such a result also suggests that mass would be effectively impossible to identify using only the embedded sensors. (a) The uncertainty of θ 0 influences only the magnetometer measurements H x , H z , especially in the beginning. Any attempt to identify it using accelerometers would be futile. The observed effects are especially interesting in the context of identification of unknown parameters. Knowing that non-influential parameters are non-identifiable [START_REF] Dobre | Limits of variance-based sensitivity analysis for non-identifiability testing in high dimensional dynamic models[END_REF], one can easily confirm that some of the model parameters might be identifiable using only some of the sensors, and only in some periods of time. For example, any attempt to identify initial Euler angles θ 0 and ψ 0 using magnetometer measurements H x , H y , H z after a few seconds of flight would be futile. Moreover, the conformity of the results with the domain knowledge combined with the previously described numerical accuracy of the quasi-LPV model [START_REF] Machala | The Analysis of Vehicle's In-Flight Behaviour Using Quasi-LPV and Nonlinear Models[END_REF], suggests that the quasi-LPV model has been correctly developed, and effectively preserves the real behavior of a projectile. It can be also concluded that the model with reduced uncertain parameter vector becomes more additive, as can be seen in Fig. 6 Overall tendencies of the sensitivity indices have been summarised in Tables 4 and5.

It may be noted that convergence of sensitivity indices estimates is assured through the use of sufficiently large number of samples; that in turn is rendered possible through the use of computationally efficient quasi-LPV model developed in a non-rolling reference frame. Should the model complexity be higher and calculations prohibitively computationally expensive, an approach of calculating confidence intervals around the estimated could be utilised. An example of such a technique is presented in [START_REF] Gilquin | Making best use of permutations to compute sensitivity indices with replicated orthogonal arrays[END_REF]. Alternatively, different sampling technique could be applied, as it can lead to a better coverage of the parameter space and can accelerate the convergence of the sensitivity indices, without requiring a large number of samples [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF].

V. Conclusions

In the paper, a quasi-LPV reformulation of a nonlinear first-principles-based model of a projectile in free-flight is described, on top of which a global sensitivity analysis is applied (GSA). The suggested approach aims at using the proposed quasi-LPV model-which is much faster due to the applied non-rolling reference frame-rather than the nonlinear one, in order to shorten the duration of each simulation. For GSA, a permutation-based sampling is used, in order to reduce the total number of simulations.

Benefits of GSA are twofold. Firstly, it highlights the possible identifiability limitations. Since non-sensitive parameters are non-identifiable, the approach has allowed to determine which parameters of the model can be deemed non-identifiable, either during the whole flight (i.e. mass of the projectile), or during some periods of it (e.g. initial Euler angles after the first seconds of flight). This aspect of work can be considered as a practical approach to gain more insight into the parametric variations within the model. Secondly, in the context of the present work, GSA is considered average, oscillatory sinusoidal increase G y dominates, oscillatory dominates, oscillatory G z dominates, oscillatory dominates, oscillatory as a tool to analyse the behavior of the model. That, in turn, allows to confirm that the quasi-LPV model-which has previously been tested only for its numerical precision-correctly represents the expected physical behavior of a studied system.

The future work will concentrate on applying GSA as a standardized tool for performing pre-experimentation procedure of different projectile architectures in the context of free-flight experiments conducted at the French-German Research Institute of Saint-Louis. As such, it will be used in tandem with uncertainty analysis that aim to evaluate the maximum dispersion of the impact points of the projectiles.

ρ = air density, kg m - 3 f

 3 i,aer o , m i,aer o = atmospheric force and moment, acting in i-th axis x, x = variable in variable-roll frame, non-rolling frame X i = auxiliary nonlinear function x, y = state vector, output vector s = number of uncertain model parameters s•, c•, t• = sine, cosine, and tangent function of • S = reference surface, m 2 S j i = main effect Sobol sensitivity index for i-th uncertain parameter and j-th output S j T i = total effect Sobol sensitivity index for all terms involving the i-th uncertain parameter for j-th output Ω i = i-th uncertain model parameter Ω, Ω 1 , Ω 2 = vector of uncertain parameters, its first and second sample χ = an elementary random event d = reference diameter, m m = mass, kg

  the normal force slope C N α and Magnus force slope C ypα coefficients. Similarly, the moments are expressed through the roll damping C lp , pitch moment coefficient slope C mα , pitch damping C mq , and Magnus moment slope C npα coefficients [29].

Fig. 1

 1 Fig. 1 Diagram of the global sensitivity analysis.

  Sensitivity indices of each uncertain parameter with respect to A y . Influence of translational and angular velocities dominates over all the other parameters. Main sensitivity indices of the output sensors with respect to C X .The influence of model uncertainties on sensors other than the axial accelerometer A x is visibly low.

Fig. 2

 2 Fig. 2 Exemplary results: a proof of the domination of angular and translational velocities over all the other uncertain parameters.

Fig. 3

 3 Fig. 3 Sum of main sensitivity indices for each output signal. The results suggest the additive nature of the model.

  The uncertainty of φ 0 influences the magnetometer measurements H x , H y , especially in the beginning. Any attempt to identify it using accelerometers would be futile.

Fig. 4

 4 Fig. 4 Sensitivity indices with respect to the initial Euler angles.

  (b) in comparison with Fig. 3. Additive nature of the model depends heavily on the selection of uncertain parameters and different selection might have resulted in a non-additive structure-in such case, a separate calculation of total sensitivity indices would have been necessary, e.g. using standard Monte Carlo approach [54]. It is worth noting that the correctness of the obtained sensitivity indices evolution can be additionally tested by calculating the sensitivity indices for denser time-grids. For the studied model, results remain equivalent to those obtained for a lower density of time instances.

  C X influences especially the axial accelerometer A x , as well as the radial measurements A y , A z in the instances of time when they are not influenced by C N α . C N α influences especially the radial accelerometer A y , A z , with an increasing over time influence on the magnetometers H x , H y , H z .

Fig. 5

 5 Fig. 5 Sensitivity indices with respect to the aerodynamic coefficients.

  Main sensitivity indices of the output sensors with respect to m. Mass uncertainty influences the sensors only slightly. An attempt to identify it using in-flight measurements would be difficult. Sum of main sensitivity indices. Results are approximately equal to one, indicating additive nature of the model.

Fig. 6

 6 Fig. 6 Exemplary results: a proof of the additive nature of the model and non-identifiability of the mass.

Table 1 Uncertain initial state variables

 1

		Nominal value
	State	Scenario A	Scenario B Scenario C Scenario D Uncertainty	Unit
	u 0	300	400	500	600	±1%	m s -1
	q 0	0	0	0	0	±3	rad s -1
	r 0	0	0	0	0	±3	rad s -1
	θ 0	0.0602	0.0370	0.0251	0.0188	±0.00049	rad
	ψ 0	-0.0009	-0.0008	-0.0006	-0.0005	±0.00049	rad

Table 4 Sensitivity of uncertain state variables.

 4 A x A y A z H x H y H z G y G z -parameter non-influential and non-identifiable, ds -sensitivity decreasing over time.

	θ 0	ni	ni	ni	ds	ni	ds	ni	ni
	ψ 0 ni	ni	ni	ds	ds	ni	ni	ni
	m	ni	ni	ni	ni	ni	ni	ni	ni
	Legend: ni

Table 5 Sensitivity of uncertain model parameters.

 5

		C X	C N α
	A x	dominant	non-identifiable
	A y	sinusoidal decrease	sinusoidal increase
	A z dominates, oscillatory	sinusoidal increase
	H x	average, oscillatory	sinusoidal increase
	H y	average, oscillatory	sinusoidal increase
	H z

Funding Sources

This work is in part financed by the French Defense Procurement Agency (DGA, Direction Générale de l'Armement) of the French Ministry of the Armed Forces.