Joël Serra

Christophe Bouvet
email: christophe.bouvet@isae.fr

Bruno Castanié

Caroline Petiot

Experimental and numerical analysis of Carbon Fiber Reinforced Polymer notched coupons under tensile loading

Keywords:

a b s t r a c t

The behavior of composite laminates subject to notch based stress concentrations is difficult to apprehend, especially the mechanisms of damage progression leading to total failure. Numerical and experimental investigations were carried out on three different stacking sequences of notched, thin ply carbon/epoxy laminates. This paper presents a computational study of notched tensile tests (U-notch) using the Discrete Ply Modeling (DPM) method, which has already proved efficient on both in-plane and out-of-plane loading cases, such as pull through, low velocity impact and compression after impact. The specificities of this finite element model are its discrete nature (interface elements to model delamination and matrix cracks), the small number of parameters required, and its robustness. This work follows on from the study of open-hole tensile tests (same three layups) by the same authors [1] and analyzes the influence of layup and notch shape. Comparisons with experiments (using infrared technology) demonstrate that tensile strengths, and failure scenarios and patterns are predicted with acceptable accuracy.

Introduction

Composite materials are now largely used in most transport vehicles, especially in the aerospace and astronautics fields [START_REF] Gay | Composite materials: design and applications[END_REF]. To demonstrate the strength of a fuselage, aircraft manufacturers must prove that it can sustain extensive damage called a ''2-bay crack" (Fig. 1). Thus notched strength is one of the main design drivers for aeronautic composite structures. In order to evaluate the residual strength of composite structures with a hole or notch, two main types of method can be used; stress-based methods and energy-based methods.

Literature review

The stress-based methods use the stress field in the vicinity of the hole/notch and the energy-based methods use the evolution of the Strain Energy Release Rate (SERR) when the crack propagates. For example, a stress-based method, the well-known ''point stress" failure model [START_REF] Awerbuch | Notched strength of composite laminates: predictions and experiments-a review[END_REF][START_REF] Tan | Stress concentrations in laminated composites[END_REF] uses the stress at a distance d 0 from the hole/notch in order to apply a failure criterion and to evaluate the structure failure. The physical meaning of this distance d 0 is the damaged area at the hole/notch tip that smooths the stress field. A similar distance is used in the ''average stress" failure model [START_REF] Awerbuch | Notched strength of composite laminates: predictions and experiments-a review[END_REF][START_REF] Tan | Stress concentrations in laminated composites[END_REF]. For this model, the stress is averaged over a distance d 0 , called the ''average distance", in order to smooth the stress field calculated in the neighborhood of the hole/notch. This distance d 0 allows simple elastic Finite Element (FE) calculations to be used, and avoids the need for complex calculations taking the composite damage into account. These methods also make it possible to simulate the size effect [START_REF] Bazant | The scaling of structural strength[END_REF][START_REF] Bazant | Size effect and fracture characteristics of composite laminates[END_REF][START_REF] Wisnom | Scaling effects in notched composites[END_REF][START_REF] Xu | An experimental investigation into size effects in quasi-isotropic carbon/epoxy laminates with sharp and blunt notches[END_REF][START_REF] Camanho | Prediction of size effects in notched laminates using continuum damage mechanics[END_REF], even though d 0 can depend on several factors [START_REF] Kassapoglou | Modeling the effect of damage in composite structures: simplified approaches[END_REF], such as the stacking sequence, material, hole/notch size, etc. In particular, the size effect leads to the result that a small hole is less penalizing than a large hole (for homothetic size of structures). If the average distance is taken to be constant, d 0 will be proportionally greater for a small hole than for a large one.

Numerous experiments have been conducted to explore the physical meaning of the size effect [START_REF] Bazant | The scaling of structural strength[END_REF][START_REF] Bazant | Size effect and fracture characteristics of composite laminates[END_REF][START_REF] Wisnom | Scaling effects in notched composites[END_REF][START_REF] Xu | An experimental investigation into size effects in quasi-isotropic carbon/epoxy laminates with sharp and blunt notches[END_REF][START_REF] Camanho | Prediction of size effects in notched laminates using continuum damage mechanics[END_REF]. For example, Wisnom et al. [START_REF] Wisnom | Scaling effects in notched composites[END_REF] studied the size effect on isotropic specimens with thin plies. They showed that the strength of notched composite laminates decreased when the notch size increased and concluded that the hole size effect was triggered by the presence of non-critical damage at the ply level, such as fiber failure, ply splitting, delamination or matrix cracking in the vicinity of the hole, which smooths the stress concentration. The ratio of the size of this ''fracture process zone" to the specimen size explains the strength difference between small and large specimens.

Energy-based methods use the SERR (or the strain intensity factor) at the crack vicinity in order to evaluate whether the crack propagates and whether the propagation is stable [START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF]. The same concept of average distance can also be used to take account of the damage developing at the crack tip without performing complex calculations involving the composite damage [START_REF] Kannan | Failure assessment on central-sharp notched carbon/epoxy laminates[END_REF][START_REF] Khashaba | Fracture Behavior of Woven Composites Containing Various Cracks Geometry[END_REF][START_REF] Hochard | A ply scale non-local fibre rupture criterion for CFRP woven ply laminated structures[END_REF].

A lot of models in the literature are based on stress-based methods and energy-based methods with an average distance concept [START_REF] Awerbuch | Notched strength of composite laminates: predictions and experiments-a review[END_REF][START_REF] Tan | Stress concentrations in laminated composites[END_REF][START_REF] Kassapoglou | Modeling the effect of damage in composite structures: simplified approaches[END_REF][START_REF] Waddoups | Macroscopic fracture mechanics of advanced composite materials[END_REF][START_REF] Yeh | Analytical damage zone size predictions for notched laminated composites[END_REF][START_REF] Zhen | The D criterion theory in notched composite materials[END_REF][START_REF] Chen | Prediction of the strength of notched fiber-dominated composite laminates[END_REF][START_REF] Eriksson | Strength of tensile loaded graphite/epoxy laminates containing cracks, open and filled holes[END_REF][START_REF] Vieille | Prediction of the notched strength of woven-ply PolyPhenylene Sulfide thermoplastic composites at a constant high temperature by a physically-based model[END_REF]. However, it can be difficult to assess whether a stress-based method or an energy-based method is more suitable to evaluate the residual strength of a given holed/notched composite structure. In addition, the average distance (d 0 ) is often difficult to evaluate and can change with the structure size. Eventually, it is necessary to develop models able to evaluate the residual strength whatever the size and the shape of the hole/notch. Camanho et al. proposed an alternative approach to evaluate the residual strength of composite laminates with a hole or crack [START_REF] Camanho | Prediction of size effects in notched laminates using continuum damage mechanics[END_REF]. This analytical approach relies on a both stress-and energy-based method. This theory, initially introduced by Leguillon [START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavities -size effects[END_REF], predicts failure when both criteria are met. Compared to classic stress-based or energybased methods, the main interest of this approach is to evaluate the failure strength of the structure and the average distance at the same time, avoiding having to choose the average distance a priori. In the same way, Y. Mohammed et al. used cohesive laws [START_REF] Mohammed | Size effect analysis of open-hole glass fiber composite laminate using two-parameter cohesive laws[END_REF] to determine the strength of an open-hole composite laminate sample. This model, based on Dugdale's [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt's [START_REF] Barenblatt | The mathematical theory of equilibrium cracks brittle fracture[END_REF] cohesive law concept links the failure strength to the length of the process zone. These fast analytical and numerical methods are very useful and easy to use but often do not take the influence of the stacking sequence into account since the strain field is assumed to be constant within the laminate thickness (due to the plate model) [START_REF] Camanho | Numerical modelling of failure in advanced composite materials[END_REF]. However, experimental results show substantial influence of the stacking sequence on both failure scenarios and strengths [START_REF] Hallett | Experimental investigation of progressive damage and the effect of layup in notched tensile tests[END_REF][START_REF] Hessabi | Effects of stacking sequence on fracture mechanisms in quasi-isotropic carbon/epoxy laminates under tensile loading[END_REF]. For example, the final failure of composite laminates with thick plies is mainly driven by delamination propagation [START_REF] Wisnom | Scaling effects in notched composites[END_REF][START_REF] Wisnom | Modelling discrete failures in composites with interface elements[END_REF]. This stacking sequence effect was also observed in the present study, where the failure strength of plain specimens was observed to vary by about 30% between the different stacking sequences of laminates with the same number of plies in each direction (Table 2). Moreover, more complete numerical models are needed to consider the damage (delamination, matrix cracking, fiber failure, etc.) developing in the vicinity of the hole/notch and to predict the final failure strength. Many models have been proposed in the literature [START_REF] Serra | Scaling effect in notched composites: the discrete ply model approach[END_REF][START_REF] Awerbuch | Notched strength of composite laminates: predictions and experiments-a review[END_REF][START_REF] Camanho | Prediction of size effects in notched laminates using continuum damage mechanics[END_REF][START_REF] Kassapoglou | Modeling the effect of damage in composite structures: simplified approaches[END_REF][START_REF] Hochard | A ply scale non-local fibre rupture criterion for CFRP woven ply laminated structures[END_REF][START_REF] Wisnom | Modelling discrete failures in composites with interface elements[END_REF][START_REF] Adam | Discrete ply model of circular pull-through test of fasteners in laminates[END_REF][START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF][START_REF] Ladevèze | On a damage mesomodel for laminates: micro-meso relationships, possibilities and limits[END_REF][START_REF] Naik | Notched strength of fabric laminates I: prediction[END_REF][START_REF] Kim | Notched strength and fracture criterion in fabric composite plates containing a circular hole[END_REF][START_REF] Xiao | Damage and fracture of notched non-woven and woven composite laminates[END_REF][START_REF] Belmonte | Characterisation and modelling of the notched tensile fracture of woven quasi-isotropic GFRP laminates[END_REF][START_REF] Afaghi-Khatibi | An effective crack growth model for residual strength evaluation of composite laminates with circular holes[END_REF] to simulate damage and to predict final failure in composite structures. The 3 main types of damage (fiber failure, matrix cracking and delamination) must be taken into account in order to accurately simulate damage in composite laminates and, in particular, to simulate the process zone in the vicinity of the hole/notch. For example, Camanho et al. [START_REF] Camanho | Numerical modelling of failure in advanced composite materials[END_REF] developed a model, based on continuum damage mechanics, able to simulate the onset and propagation of the main damage types. Abisset et al. used the damage meso-model [START_REF] Ladevèze | On a damage mesomodel for laminates: micro-meso relationships, possibilities and limits[END_REF] to model diffuse damage and matrix cracking with a progressive damage law, whereas a fracture mechanics law was used for fiber failure [START_REF] Abisset | On the validation of a damage mesomodel for laminated composites by means of open-hole tensile tests on quasi-isotropic laminates[END_REF]. The failure modes of a notched specimen subjected to tensile loading were predicted with good accuracy, but a large number of parameters was needed, which could be difficult to evaluate. Cohesive interface elements are also often used to simulate delamination, and sometimes matrix cracking. For example, Pinho et al. developed a model able to predict fiber and matrix cracking using smeared crack models and delamination using interface elements [START_REF] Pinho | Modelling failure of laminated composites using physically-based failure models[END_REF][START_REF] Chen | Numerical analysis of size effects on openhole tensile composite laminates[END_REF]. ONERA is developing a model based on a multi-scale progressive failure approach, to describe the softening behavior of the fiber failure, and interface elements, to account for delamination [START_REF] Laurin | A multiscale hybrid approach for damage and final failure predictions of composite structures[END_REF]. Similarly, Ridha et al. built a progressive damage model based on the different damage mechanisms [START_REF] Ridha | Modelling complex progressive failure in notched composite laminates with varying sizes and stacking sequences[END_REF], in which continuum damage mechanics is used to represent the in-plane damage and cohesive elements for delamination. Cohesive interface elements are also sometimes used to represent matrix cracking. For example, Wisnom et al. used Weibull's approach to represent fiber failure, cohesive interface elements to model delamination, and other interface elements inside plies to model splitting at the edge of the hole [START_REF] Hallett | An experimental and numerical investigation into the damage mechanisms in notched composites[END_REF][START_REF] Sun | Interaction of inter-and intralaminar damage in scaled quasi-static indentation tests: Part 2-Numerical simulation[END_REF][START_REF] Wisnom | The role of delamination in failure of fibre-reinforced composites[END_REF]. This work raises the question of whether there is a need for cohesive interface elements to accurately simulate matrix cracking and its interaction with delamination [START_REF] Wisnom | The role of delamination in failure of fibre-reinforced composites[END_REF]. The Discrete Ply Model (DPM), used in this paper, was one of the first models to use cohesive interface elements to simulate matrix cracking and delamination [START_REF] Bouvet | Low velocity impact modeling in laminate composite panels with discrete interface elements[END_REF][START_REF] Ostré | Edge impact modeling on stiffened composite structures[END_REF]. Other models of the same type exist in the literature [START_REF] Bergan | Assessment of damage containment features of a full-scale PRSEUS fuselage panel[END_REF][START_REF] Wisnom | Scaling effects in notched composites[END_REF][START_REF] Lammerant | Modelling of the interaction between matrix cracks and delaminations during impact of composite plates[END_REF][START_REF] Aymerich | Prediction of impact-induced delamination in cross-ply composite laminates using cohesive interface elements[END_REF][START_REF] De Moura | Modelling the interaction between matrix cracking and delamination in carbon-epoxy laminates under low velocity impact[END_REF].

Purpose of the study

The objective of this paper is to show that it is possible to simulate two very different tests including stress concentrations using the same model. The first analysis deals with the open hole tensile test and the scaling effect. This work was presented in a previous paper [START_REF] Serra | Scaling effect in notched composites: the discrete ply model approach[END_REF] and showed that the DPM was able to simulate the open hole tensile test with 3 different hole diameters (1 mm, 3.175 mm and 6.35 mm). The failure of this test is clearly a stress-based failure, driven by the stress field in the vicinity of the hole and by the average distance d 0 . The scaling effect observed in this case can be seen as an interaction of the process zone with the hole size. The higher the ratio of the average distance d 0 to the hole diameter is, the more important is the smoothing of the stress around the notch and the higher is the failure strength. In order to simulate this test and to obtain the size effect numerically, it is necessary to simulate damage of different types, such as matrix cracking, delamination or fiber breakage, and the interaction between these different types of damage is of great importance if the process zone is to be simulated accurately.

The second analysis, which is the topic of the present paper, deals with the notch tensile test. It is the goal of the present work to study the development of damage during a notch tensile test and to compare it with that found in the open hole tensile test. The challenge is both experimental and numerical. Experimentally, it concerns following the crack propagation during the test. Infra-Red Thermography (IRT) was used for this purpose and showed very interesting results. Numerically, the challenge is to use the DPM to simulate a stress-based test (the open hole tensile test) and an energy-based test (the notch tensile test). In order to simulate these different tests, it is necessary to take the process zone, and in particular matrix cracking, delamination and their interaction, into account. It is of equal importance to take the crack propagation into account, and in particular the energy dissipated by the tensile fiber failure.

Analytical results

''U-notch" specimen failure is clearly energy-based failure that is driven, to a first approximation, by the evolution of the SERR when the crack propagates. The problem is not trivial and a process zone is created in front of the notch at the beginning of the test and also in front of the crack initiated from this notch propagation during the test. In the present case, the propagation is stable during the first part of the experiment. The stability of the crack propagation is classically driven by the evolution of the SERR with the crack length: if the curve is increasing, the propagation is unstable and, if the curve is decreasing, the propagation is stable. Nevertheless, this stability criterion should be considered with caution, because it is true only if the SERR curve versus the crack length is accurate. This curve can be difficult to evaluate because the process zone smooths the stress at the crack tip; it moves with the crack propagation and its size and morphology can change. The creation of the process zone explains the dissipation of energy during the crack propagation and the evolution of its size explains the R-curve: the farther the crack propagates, the larger is the size of the process zone, the higher is the energy needed to propagate the crack and the higher is the critical SERR.

In order to highlight the failure scenario of the 2 test types studied here, the open hole and the notched tensile test, the curves of the SERR (in the framework of linear fracture mechanics) are plotted versus the crack length in Fig. 1. These curves, and the physical meaning extracted from them, should be taken with caution because they were elaborated with a simple elastic FE calculation with plate theory and the process zone was not simulated.

It is very important to note that, in this figure, we chose to plot the SERR for the failure stress of each test. Different values of failure stress were obtained for the 4 experiments: 0.63 (dimensionless stress) for the 1 mm-diameter hole, 0.59 for the 3.175 mmdiameter hole, 0.53 for the 6.35 mm-diameter hole and 0.4 for the notched specimen (these stresses were evaluated as the average of the failure stress for the 3 stacking sequences given in Table 2). This choice made it possible to highlight whether the final failure was more stress-based or energy-based. It is clear in this figure that the SERR of the open hole tensile test was very low, and thus the final failure should be stress-based, while the SERR of the notched tensile test was near the fracture toughness of tensile fiber failure, and thus the final failure should be energy-based. Nevertheless, the fracture toughness of tensile fiber failure is difficult to evaluate and a few values are available in the literature. We considered about 80-100 N/mm for this type of material [START_REF] Laffan | Measurement of the in situ ply fracture toughness associated with mode I fibre tensile failure in FRP. Part II: size and lay-up effects[END_REF][START_REF] Pinho | Fracture toughness of the tensile and compressive fibre failure modes in laminated composites[END_REF].

A few additional remarks may help with the interpretation of the curves of the open hole specimens. The curve of the 1 mmdiameter hole tends towards high values when the crack length tends to 2 mm. This is, of course, because the width to diameter ratio equals 5, so the ligament size equals 2 mm. Similar results were obtained at 6.35 mm for the 3.175 mm-diameter hole and at 12.7 mm for the 6.35 mm-diameter hole. The SERR curves of the holed specimens are strongly increasing and explain why the final failure scenario of these tests is unstable and brutal. This brutal final failure was obtained experimentally, except for 2 samples of 6.35 mm-diameter, where IRT showed the propagation of a small, stable crack. These experimental results should be taken with caution because it is more difficult to observe crack propagation in small specimens than in large ones.

The curves of the 3 hole sizes (Fig. 2) also explain why small specimens are more resistant than larger ones. Furthermore it should be kept in mind that, to evaluate the failure using these curves, an average distance d 0 should be used and a value of d 0 of about 1 mm would be appropriate as a first approximation [START_REF] Awerbuch | Notched strength of composite laminates: predictions and experiments-a review[END_REF][START_REF] Tan | Stress concentrations in laminated composites[END_REF]. The difference between the SERR curves should be assessed bearing in mind that the failure stress was used for the FE calculations. Then, as the failure stress is higher for the small hole than for the medium sized one, the SERR curve of the small hole would have been lower if the calculation had been made with the failure stress of the 3.175 mm-diameter hole (the SERR varies as the square of the stress). A comparison of the SERR curves of the 3 holes with different diameters for a given stress could clearly explain the hierarchy of the failure stress, except that the failure is stress-based and not energy-based! Consequently, the curve of the SERR should not be used to evaluate the failure strength but the stress field should be used. (It gives similar conclusions.)

The open hole tensile test is, however, not the topic of this paper and this article focuses on analyzing notched specimens. The SERR curve provides us with very interesting information. This curve is strongly increasing for the first 0.4 mm of crack length and becomes flatter thereafter (Fig. 2). This strong increase is, of course, linked to the notch end radius, which was equal to 0.5 mm in this study. The SERR curve of a perfect crack was plotted for comparison and the 2 curves are very similar as the crack length is more than 0.5 mm. This result highlights why no significant influence of the notch end radius is found experimentally, if the radius is small enough. After this first strong increase, the SERR increases more gradually. In fact, during the experiment, the crack propaga-tion was stable up to a crack length of about 10 mm. This result is not surprising seeing the SERR curve, considering that this curve should be taken with caution because the FE calculation is elastic and the process zone is not taken into account (and should induce an R-curve effect).

Experimental work

Material and setup

The material investigated in the present study was Hexcel's T700-M21 carbon epoxy unidirectional laminate with a nominal ply thickness of 0.125 mm. Three symmetric stacking sequences of 13 plies were studied. They presented the same number of plies in each direction (0°, 90°and ±45°), and only the relative position of plies changed between layups:

-C3-1 [45/-45/X/X/X/90/0/90/X/X/X/-45/45] -C3-2 [X/X/X/X/0/90/0/90/0/X/X/X/X] -C3-3 [X/X/X/X/X/0/0/0/X/X/X/X/X] For confidentiality reasons, the draping sequences cannot be fully disclosed in this article. The sample studied was 180 mm wide and 300 mm long (not taking the grips into account) with a 30 mm central notch having an end notch radius of 0.5 mm (Fig. 3). The ratio of sample width to notch length was 6, so slightly higher than for holed specimens (W/D = 5), and the SERR curve was clearly higher because of this ratio and the notch shape (Fig. 2).

Tests were performed at ambient temperature (about 20 °C) using a 450 kN hydraulic machine at a displacement control speed of 0.02 mm/s. During the test, displacement was measured using an LVDT sensor and load was measured using the load machine's sensor. Two stereo-correlation cameras were also positioned on one side of the specimen and an IRT camera on the other side (Fig. 4). Images were analyzed with VIC 3D for the stereocorrelation and Altair Ò for the IRT measurements. Unfortunately, the coupons slipped in the machine grips during the tests, so the strains shown in this paper were evaluated using displacement from Digital Image Correlation (DIC) measures over a distance of 160 mm (Fig. 3).

Four specimens were tested to total failure for the C3-1 and C3-3 layups, and three for C3-2. Different configurations linked to the camera lenses (close-up/wide shot, summarized in Table 1), were used in order to obtain the maximum amount of information on the fracture scenario. Stress/strain curves (Fig. 6) were consequently derived from configurations in which the DIC cameras displayed a large part of the specimen (wide shot). For close-up configurations, only failure loads were determined and are schematized with arrows in Fig. 5. Experimental setup errors on DIC camera parameters for specimen n°9 and on the IRT camera parameters for specimen n°7 made it impossible to use these images.

Results

Similarly to the experimental results obtained in the previous study of open hole tension tests [START_REF] Serra | Scaling effect in notched composites: the discrete ply model approach[END_REF], the present experimental results (Table 2) showed little discrepancy for the three stacking sequences. In Table 2, ''Plain" refers to plain samples used to compare results with holed/notched samples, ''Small" refers to the 1 mm-diameter hole, ''Medium" to the 3.175 mm-diameter hole, ''Large" to the 6.35 mm-diameter hole, ''Notched" to the notched samples, r 1 to the mean remote failure stress and CV to the coefficient of variation. To respect privacy policy, the level of stress and strain was normalized by the maximum failure stress reached by one coupon of the C3-3 layup. On average, C3-2 and C3-3 layups exhibited a total failure load slightly larger than that of the C3-1 layup. Differences between stacking sequences were much reduced compared to plain and open-hole tension tests. In comparison with holed specimen tests, failure loads were much lower. No structural failure (B) -defined as the first occurrence of load or stiffness dropping by more than 5% -was detected before total failure, except in coupons n°1 and n°8.

2.2.1.

Oriented laminate C3-1 2.2.1.1. Failure scenarios and patterns. The first damage observed on the C3-1 layup from point A (Fig. 5) could be fiber/matrix debonding, matrix cracking, delamination or fiber breakage. The most energetic failure mode is fiber fracture. A local increase in temperature of about one degree (Fig. 6) seems to indicate such a failure mode [START_REF] Lisle | Damage of woven composite under translaminar cracking tests using infrared thermography[END_REF]. The others are less energetic: they increase the temperature locally by a few tenths of a degree [START_REF] Lisle | Damage of woven composite under translaminar cracking tests using infrared thermography[END_REF][START_REF] Lisle | Damage analysis and fracture toughness evaluation in a thin woven composite laminate under static tension using infrared thermography[END_REF]. Once the propagation had begun, emissions of heat were detected on the right side of the notch (close-up for coupons n°1 , 2 and 3, Fig. 6). They appeared to be relatively symmetrical on either side of the notch (wide shot for coupon n°4, Fig. 6). The damage then continued to progress up to a brutal total failure. The crack propagated in an explosive manner, oriented at a few degrees with respect to the horizontal. After rupture of the laminate, a symmetrical pull-out failure pattern was obtained (Fig. 7).

Crack propagation.

In order to allow the application of a new criterion for comparing the numerical results with the experimental results, two attempts were made to determine the evolu-tion of the crack length during the test. The first used DIC, the second, IRT.

Using the DIC technique, we determined the longitudinal strain (e yy ) field near the crack (Fig. 8) at different progressions of the tensile test (41, 43 and 45% in Fig. 9). Test progression (e average )was defined by the average strain of the specimen normalized by the maximum strain obtained (plain specimen of C3-3 layup). The longitudinal strain was compared with fiber failure strain (e T 0 , Table 3) and, assuming the strain to be constant throughout the thickness of the laminate, it was possible to estimate the extent of broken fibers and then the crack length (along x). The relative position indicated (Figs. 9 and 10) is the position/ length-of-ligament ratio; 0 corresponding to the notch edge and 1 to the free edge of the test coupon. Longitudinal strain (e yy ) is normalized by the tensile fiber failure strain (e T 0 Þ to become e normalized yy (Fig. 8). The purple area corresponds to strain (e normalized yy Þ below 1, so, to a first approximation, to the zone without fiber fracture (in the y-direction).

It is then possible to obtain, from Fig. 8, the evolution of the fiber failure extent as a function of the normalized average strain (e average ) enforced. This evolution is shown in Fig. 9. However, it is necessary to consider the evolution of the damage obtained with care since it was measured on the outside ply. Delamination between the various plies made it impossible to have a constant strain throughout the thickness of the laminate; a different strain of the 0°inside the thickness could be observed. Moreover, it is likely that the ply observed by the DIC cameras entered the inelas- tic domain and, from then on, a measured strain greater than the tensile fiber failure strain does not necessarily denote fracture.

Composite laminates were then monitored with the help of an IRT camera in order to determine notch extension using a different methodology. Fiber rupture is very energetic compared to other types of damage but it was possible to track the crack progression by listing all the successive positions of heat emissions (Fig. 6). Thus the results given by the two methods could be compared (Fig. 10). IRT emission monitoring showed that damage occurred early in the test. In Fig. 10, D corresponds to the stage of the tensile test when the DIC showed a longitudinal strain (near the notch) greater than the fiber failure strain and therefore gives an estimation of the crack propagation extent. Before D, notch extension must be small; hence it is hidden inside the laminate with delamination around, which explains the absence of detection by DIC. After D, damage (fiber rupture, matrix cracking and delamination) becomes considerable at the notch edge and generates a longitudinal strain greater than the tensile fiber failure strain on the outside ply.

The notch extension determination method using IRT seems to be more relevant than the DIC process since the IRT does not assume a constant strain throughout the thickness of the laminate. However, it must slightly overestimate the crack progression due to the propagation of heat in the laminate. When a fiber from a 0°ply far from the outside surface breaks, the heat emitted propagates towards the surface until it is detected by the IRT camera. During this period of time, it also propagates within the plane of the laminate.

Oriented laminates C3-2 and C3-3

Failure scenarios associated with layups C3-2 and C3-3 were very similar to the ones observed for layup C3-1. In Fig. 5, from A, damage progresses from both notch ends and leads to total failure of the laminate.

The failure patterns shown in Figs. 11 and 12a reveal 45°oriented crack propagations (''Z shaped") for the two stacking sequences. Closer examination suggests that 0°fiber rupture has actually propagated along a direction approximately orthogonal to the tensile direction (Fig. 12b).

Afterwards, notch length evolutions (Figs. 13 and14) were determined using the same IRT process as for C3-1 layup.

Crack propagation seems to initiate sooner for C3-1 layup than for C3-2 layup, and sooner for C3-2 layup than for C3-3 layup. Maximum stable crack propagation (just before brutal total failure) varies from 5% to 12% of the ligament (75 mm long) without significant differences between the three layups.

Conclusions drawn from the experiments

Regarding all tensile tests (''Plain", ''Small", ''Medium", ''Large" and ''Notched"), C3-3 layup appears to be the most resistant (it exhibits the highest failure stress, Fig. 15). Differences between layups are, however, reduced when failure stress decreases (edge effects become less critical).

Net failure stresses obtained for the ''Notched" specimens are much smaller than open-hole ones (Fig. 15). This discrepancy is mainly related to the difference in SERR between notches and holes (Fig. 2).

In order to assess notch sensitivity independently of the stacking sequence effect (20% failure stress difference for ''Plain" speci- mens), the failure stress of each configuration (''Small", ''Medium", ''Large" and ''Notched") divided by ''Plain" failure stress (for each layup) is summarized in Fig. 16. The layups C3-1 and C3-2 appear to be less sensitive to the presence of the cut than C3-3.

The justifications developed in [START_REF] Serra | Scaling effect in notched composites: the discrete ply model approach[END_REF] for the open-hole tension case seem to be valid here too. Edge effects, triggering delamination for ''Plain" specimens observed under tension, have the opposite effect around a notch (they smooth stress concentration). C3-3 layup is not much subject to free edge delamination (''Plain" specimens). Hence it does not smooth stress concentrations near the notch and has a lower ratio (Fig. 16) than C3-1 and C3-2 layups.

Numerical modelling

DPM principle

The Discrete Ply Model has already been described in several papers [START_REF] Adam | Discrete ply model of circular pull-through test of fasteners in laminates[END_REF][START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF][START_REF] Bouvet | Low velocity impact modeling in laminate composite panels with discrete interface elements[END_REF][START_REF] Ostré | Edge impact modeling on stiffened composite structures[END_REF] and the same characteristics were kept here as in [START_REF] Serra | Scaling effect in notched composites: the discrete ply model approach[END_REF] for scaled open-hole tension test modelling. The main features, such as mesh construction and behavior law, are briefly recalled below. The principle of the DPM model is to simulate the major failure modes observed in composites as follows (Fig. 17):

Delamination is taken into account using classical interface elements between two consecutive plies (or groups of plies) of different orientation, each ply being modeled with one volumic finite element in the thickness. Then, the damage in the delamination interface elements is classically driven by fracture mechanics. Matrix cracking is taken into account using interface elements normal to the transversal direction. The damage of the matrix cracking interface elements is driven using fracture criteria evaluated in the neighboring volumic finite elements. These interface elements induce complex meshing but make it possible to automatically simulate the interaction between intraand inter-laminar damage. This interaction is crucial to account for the complex damage morphology observed in composite structures.

Fibre fracture is considered using continuum damage mechanics with an original formulation linking the integration points of the volumic finite element to impose a constant SERR per unit area [START_REF] Ostré | Edge impact modeling on stiffened composite structures[END_REF]. This approach can be compared to methods using the characteristic element length which allows mesh-size-independent modelling [START_REF] Laffan | Measurement of the in situ ply fracture toughness associated with mode I fibre tensile failure in FRP. Part II: size and lay-up effects[END_REF][START_REF] Raimondo | A progressive failure model for mesh-size-independent FE analysis of composite laminates subject to lowvelocity impact damage[END_REF][START_REF] Bazant | Progressive Crack and band theory for fracture of concrete[END_REF]. First order volume elements (8 nodes per element) were selected in order to simulate bending phenomena correctly while using only one element per ply (thickness direction). Two nodes are sufficient to represent affine behavior within the element. Therefore the strain field within the laminate thickness is a piecewise defined function.

Meshing and calculation strategies

Applying a constant size mesh to the large dimension (180Â300 mm 2 ) coupons would have generated prohibitive calculation costs. To solve this issue, different types of elements were used:

The standard elements of the DPM elements, detailed in the previous paragraph (Fig. 17), for the zone near the notch (Fig. 18-zone 3) Volumic finite elements with reduced integration (and without interface elements) in zone 2 (Fig. 18-zone 2) Standard volumic elements with a homogenized laminate behavior far from the notch (Fig. 18-zone 1)

The last two types of elements, never damaged during the propagation of the crack emanating from the notch, made the calculation faster. Numerous iterations were necessary to estimate the three zone dimensions. The blue zone (3) was reduced as far as possible without disrupting notch propagation. Connection between zone 1 and zone 2 was realized through a kinematic coupling: nodes from plies at the border between zone 1 and 2 were constrained to follow the displacement of the corresponding zone 1 element face.

To further reduce calculation costs, planar and central symmetries were used. It was verified that modelling one fourth of the coupon gave results equivalent to when the whole specimen was meshed. To set planar symmetry corresponding to the laminate symmetry plane, since the three layups studied contained an odd number of plies, the middle ply was represented with a thickness of half a ply (0.0625 mm).

Since damage progression is substantial before total failure, a customized management of the distorted elements had to be implemented. First, the calculation was divided into two parts: one step including very little damage (corresponding to the part Fig. 16. Notch sensitivity for the three layups. Fig. 17. Modelling of composite damage with the different element types [START_REF] Bouvet | Low velocity impact modeling in laminate composite panels with discrete interface elements[END_REF]. of the stress/strain curve before A in Fig. 5) and the other from the initiation of damage propagation until total failure of the laminate. For the second part, the step time value was set to one third that of the first. The Abaqus ''Distortion control" function was exploited to restrain the distortion of damaged elements through the use of an artificial energy. It was verified a posteriori that this energy was less than 1% of the internal energy. Even with those modifications, near the notch, element distortion (the first elements to be entirely damaged) made the calculation impossible. It was then decided to suppress totally failed elements, i.e. elements for which the associated energy was totally dissipated (fiber failure and delamination, cf. Fig. 17) or when the stress criterion was met (matrix cracking, cf. Fig. 17). The specimens being subject to tensile stress, deletion of these elements did not have any major influence (once ''broken", elements do not withstand any load). It would not have been possible to draw this conclusion in a compression loading case (but distortion would not have been of the same magnitude).

Once all the necessary updates had been implemented, the calculation was run on 1.4Â10 6 elements and lasted approximately 50 h (20 CPUs) using the EOS high performance calculator at CAL-MIP. The DPM uses a small number of parameters (Table 3), all of which were provided by experimental tests [START_REF] Serra | Scaling effect in notched composites: the discrete ply model approach[END_REF][START_REF] Adam | Discrete ply model of circular pull-through test of fasteners in laminates[END_REF][START_REF] Bouvet | Low velocity impact modeling in laminate composite panels with discrete interface elements[END_REF].

Model validation and discussion

Stress/strain curves obtained from numerical simulations correlated well with the experimental ones (Fig. 5). C3-3 layup stiffness loss (A-B) was well expressed by numerical simulation but failure stress relative order was not predicted correctly because failure stresses were slightly overestimated for C3-1 and C3-3 layups (Fig. 19). However, as experimental values were very close with non-negligible coefficients of variation (Table 2), reproducing this relative order was not of major importance.

The global behavior of laminates, represented by stress/strain curves, seemed to be correctly reproduced. To validate the numerical model further, notch evolution was determined for the three layups with the help of IRT technology (cf. § 2.2) and compared with simulations. Only the coupons shot in a close-up configuration (Table 1) were put to use. The others (wide-shot) did not enable a sufficiently precise length determination to be obtained. Numerical failure patterns were also compared to experimental ones.

Oriented laminate C3-1

Damage evolution of notched laminate was very progressive from A (Fig. 5). A ''butterfly wings" type evolution of matrix cracking (Fig. 20) was obtained numerically (very similar to the one observed by [START_REF] Serra | Scaling effect in notched composites: the discrete ply model approach[END_REF][START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF] on open-hole tension configurations). The 0°/90°interface delamination and 0°fiber fracture seemed to propagate perpendicularly to the tensile direction until total failure (Fig. 21).

The numerical failure pattern of Fig. 22 is very similar to the experimental one. Notch propagation is limited by the dimensions of the damageable zone (zone 3, Fig. 18) but the free edge of the coupon during experimental tests.

The evolution of the C3-1 layup notch length determined by counting the numbers of completely damaged volume elements (central ply 0°fiber failure) in the numerical simulations was very consistent with that determined using IRT. Notch propagation initiation and its maximal amplitude before total failure (8% of the ligament), were predicted with good accuracy (Fig. 23).

C3-2 and C3-3 layups

Failure scenarios and patterns simulated with the help of the DPM for C3-2 and C3-3 layups were very similar to those obtained for layup C3-1. In consequence, the numerical model failed to reproduce the 45°external plies failure pattern (''Z shape", Figs. 11 and12). Notch length evolutions were determined both numerically and experimentally (IRT) for C3-2 (Fig. 24) and C3-3 (Fig. 25) layups. The numerical model again predicted notch evolution kinematics that were very consistent with experiment.

Conclusion

The purpose of the study described in this paper was to investigate the development of damage during a notched tensile test and to compare it with the open hole tensile test performed on the same three layups. Numerical results have been confronted with experimental data and the influence of both the notch shape and dimensions and the stacking sequence have been analyzed.

Notched tensile tests confirmed the tendency already observed on open-hole specimens [START_REF] Serra | Scaling effect in notched composites: the discrete ply model approach[END_REF]: even though the C3-3 layup was more sensitive to the presence of a cut than the other two (C3-1 and C3-2), it exhibited the highest failure stress. Strength differences (hole and U-notch) were related to the shape of the cut. The SERR at a Unotched edge is higher than in the vicinity of a cracked circular hole (Fig. 2). The failure of notched coupons is hence triggered earlier and propagates within the specimen (total failure) for a stress lower than that in open-hole coupons under tensile loading. Two methods were employed to record the behavior of notched coupons, DIC and IRT. Using IRT to determine the evolution of the notched length gave more accurate results than the process based on DIC.

Modifications of the numerical model were necessary to simulate notched coupons under tensile loading at a reasonable calculation cost. Simulation outputs were validated thanks to numerous experimental observations. The DPM thus proved its ability to simulate the behavior of notched (holes of different diameter and U-notch) coupons under tensile loading, then to simulate a stress-based test (the open hole tensile test) and an energybased test (the notch tensile test). This shows that the approach is able to account for the process zone (in particular matrix cracking, delamination and their interaction) and to account for the crack propagation, (in particular the energy dissipated by the tensile fiber failure). Finally, in order to validate this numerical model at a higher level of the tests pyramid (structural level) and for different loading cases, a new test-rig, designed at the Institute Clement Ader, was used to test plates of large dimensions (400Â400 mm 2 ) and the DPM was adapted to simulate these innovating structural tests.

Fig. 1 .

 1 Fig. 1. Location of the ''2-bay crack" in the fuselage [3].

Fig. 2 .

 2 Fig. 2. SERR versus crack length for holed and notched specimens (a) and zoom on the first 2 mm (b).

Fig. 3 .

 3 Fig. 3. Notched specimen geometry (dimensions in mm).

Fig. 4 .Fig. 5 .

 45 Fig. 4. Experimental setup for the notched quasi-static tension

Fig. 6 .

 6 Fig. 6. Damage progression coupon n°1 (wide shot) and n°4 observed by IRT camera.

Fig. 7 .

 7 Fig. 7. Fracture facies of C3-1 layup.

Fig. 8 .

 8 Fig. 8. Crack extension determined using DIC.

Fig. 9 .

 9 Fig. 9. Longitudinal strain measurement near the crack using DIC.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Comparison of the two methods used to determine crack length evolution.

Fig. 12 .

 12 Fig. 12. C3-3 layup external (a) and internal (b) failure pattern.

Fig. 13 .

 13 Fig. 13. Notch length evolution (C3-2 layup) based on IRT process.

Fig. 14 .

 14 Fig. 14. Notch length evolution (C3-3 layup) based on IRT process.

Fig. 15 .

 15 Fig.[START_REF] Hochard | A ply scale non-local fibre rupture criterion for CFRP woven ply laminated structures[END_REF]. Net failure stresses for all sizes and layups.

Fig. 18 .

 18 Fig. 18. Numerical model of a notch specimen under tensile stress.

Fig. 19 .

 19 Fig. 19. Failure stress -numerical and experimental results. Fig. 20. Numerical simulation of 90°ply matrix cracking.

Fig. 20 .

 20 Fig. 19. Failure stress -numerical and experimental results. Fig. 20. Numerical simulation of 90°ply matrix cracking.

Fig. 21 .Fig. 22 .

 2122 Fig. 21. Numerical simulations of 0°/90°interface delamination and central ply 0°fi ber rupture evolutions.

Fig. 23 .

 23 Fig. 23. Notch length evolution (C3-1 layup) -numerical and experimental results.

Fig. 24 .

 24 Fig. 24. Notch length evolution (C3-2 layup) -numerical and experimental results.

Fig.

  Fig. Notch length evolution (C3-3 layup) -numerical and experimental results.

Table 1

 1 Test matrix.

	Layups	C3-1	C3-2	C3-3
	Coupon	1, 2, 3, 4	5, 6, 7	8, 9, 10, 11
	number
	IRT Camera	Close-up (1, 2, 3)	Close-up (6)	Close-up (8, 9)
		Wide shot (4)	Wide shot (5)	Wide shot (10, 11)
	DIC cameras	Close-up (4)	Close-up (5)	Close-up (10, 11)
		Wide shot (1, 2, 3) Wide shot (6, 7) Wide shot (8)

Table 2

 2 Results of plain, open-hole and notched tensile tests.

	Specimens	C3-1		C3-2		C3-3
		r 1 ð%Þ	CV (%)	r 1 ð%Þ	CV (%)	r 1 ð%Þ	CV (%)
	''Plain"	81	2.5	79	3.7	99	1.5
	''Small"	61	4.4	58	1.6	71	9.4
	''Medium"	56	1.7	54	4.5	67	0.7
	''Large"	53	0.5	50	1.3	57	4.1
	''Notched"	38	2.5	41	5.0	40	1.7

Table 3

 3 Material properties of the T700/M21 used in the Discrete Ply Model.

	Elastic Properties
	E T 1	Tensile Young's modulus in fiber direction	130 GPa
	E C 1	direction Compressive Young's modulus in fiber	100 GPa
	E 2	Transverse Young's modulus	7.7 GPa
	m 12	Poisson ratio	0.3
	G 12	Shear modulus	5.0 GPa
	Matrix cracking
	Y T	Transverse tensile strength	60 MPa
	S L	In-plane shear strength	110 MPa

Acknowledgements

The research that led to the results presented above received funds from the French National Research Agency under the VER-TEX project MATETPRO program (ANR-12-RMNP-0001).

The academic authors gratefully acknowledge the support provided by Airbus Group Innovations. This work was granted access to the HPC resources of CALMIP under allocation 2012-P1026.


Ankit Chiplunkar

Michele Colombo

Emmanuel Rachelson

Joseph Morlier

Adding Flight Mechanics to Flight Loads Surrogate Model using Multi-Output Gaussian Processes

published or not. The documents may come   L'archive ouverte pluridisciplinaire

In this paper analytical methods to formally incorporate knowledge of physics-based equations between multiple outputs in a Gaussian Process (GP) model are presented. In Gaussian Processes a multi-output kernel is a covariance function over correlated outputs. Using a general framework for constructing auto-and cross-covariance functions that are consistent with the physical laws, physics-based relationships among several outputs can be imposed. Results of the proposed methodology for simulated data and measurement from flight tests are presented. The main contribution of this paper is the application and validation of our methodology on a dataset of flight tests, while imposing knowledge of flight mechanics into the model.

Nomenclature

I. Introduction

Precise values of flight loads are highly useful information at several stages of an aircraft's life cycle. During certification, aircraft manufacturers need experimental flight loads to verify and improve their the-oretical models. Monitoring aircraft usage through loads can provide a better estimate of the remaining time-to-failure of components against their prescribed fatigue life. Aircraft flight loads are often difficult and expensive to measure experimentally. This calls for a need to develop a theoretical framework to merge different types of measurements and produce a robust flight loads surrogate model.

A. Background

Identification of an accurate flight loads surrogate model over the complete flight domain is an expensive exercise, notably because one cannot perform extensive flight testing on all the flight configurations. Since the last decade, techniques such as Neural Networks 1,2 have been applied while recently Gaussian Process Regression 3 (also known as Kriging) has shown great promise in developing the surrogate models that infer loads at untested flight configurations.

It has been proved extensively that introduction of prior information 4 improves the prediction capabilities of a GP surrogate model, diminishes chances of over-fitting and helps in identifying outliers. In recent times, techniques such as Gradient Enhanced Kriging (GEK) [START_REF] Liu | Development of Gradient-Enhanced Kriging Approximations for Multidisciplinary Design Optimization[END_REF][START_REF] Solak | Derivative observations in Gaussian Process models of dynamic systems[END_REF] have been used successfully in various engineering systems to impose prior information of gradients or local linear models, thereby improving the performance of surrogate models.

Building surrogate models on multi-output data is also known as co-kriging 7 or multi-output GP. [START_REF] Boyle | Dependent Gaussian processes[END_REF] Modelling multi-outputs is particularly difficult because we need to construct auto-and cross-covariance functions between different outputs based on the relationships between them. Recently, it has been shown that relationships between multiple outputs can be imposed by encoding it in the covariance functions [START_REF] Constantinescu | Physics-Based Covariance Models for Gaussian Processes with Multiple Outputs[END_REF] to derive models consistent with relationship and data.

B. Contribution

In this paper we extend the framework of Gradient Enhanced Kriging to integral enhanced kriging, quadratic enhanced kriging or "any relationship between outputs" enhanced kriging. Using a general framework to encode relationships between outputs 9 , we will impose prior knowledge into the Gaussian Process Regression framework and apply it for predicting flight loads.

In this paper we assume one output to be independent while the remaining outputs are linked to the first output using physics-based relationships. For instance, lift coefficient (C z ) and spanwise aerodynamic force (k z ) have an integral relationship between them equation 21. As will be described in subsection II.B, the auto-and cross-covariance functions between these outputs are evaluated exactly if the physics-based relation between them is linear 10 or using approximate joint-covariance for non-linear relationship between the outputs. [START_REF] Constantinescu | Physics-Based Covariance Models for Gaussian Processes with Multiple Outputs[END_REF] The main contribution of this paper is the application and validation of our methodology on a dataset of flight tests, while imposing knowledge of flight mechanics into the model.

II. Theoretical Setup

Let us start by defining a P -dimensional input space and a D-dimensional output space. Such that {(x j i , y j i )} for j ∈ [1; n i ] is the training dataset for the i th output. Here, n i is the number of measurement points for the i th output, while x j i ∈ R P and y j i ∈ R. We next define x i = {x 1 i ; x 2 i ; . . . ; x ni i } and y i = {y 1 i ; y 2 i ; . . . ; y ni i } as the full matrices containing all the training points for the i th output such that x i ∈ R ni×P and y i ∈ R ni .

y j i = f (x j i ) + j i ( 1
)
The science of finding a mapping function f denoted by equation 1 between {(x j i , y j i )} is called surrogate modeling in engineering design or regression in machine learning. We would use the technique of kriging (engineering design) [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF] or Gaussian Process Regression (machine learning) [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] to learn this mapping function f . A Gaussian Process is a probability distribution over infinitely many points, such that any finite subset of the process is a multi-variate gaussian distribution. A GP can be fully parameterized by a mean and covariance function equation 2.

f (x) = GP (m(x, θ), k(x, x , θ)) (2)
A random draw from a multi-variate gaussian distribution yields a random vector around its mean with values correlated through the covariance matrix. Similarly, a random draw from a GP provides a random function around the mean function m(x, θ) and of the shape as defined by the covariance function k(x, x , θ) (calculated between points x and x ). Hence, GPs provide us a method to define a family of functions whose shape is defined by its covariance function 12 also called kernel.

A. Gaussian Process Regression (GPR)

For the case of a single-output GP we assume that D equals to one, so our dataset will be {(x 1 , y 1 )}. We first start by defining a white noise measurement error between the measured outputs y and mapping function f . = N (0, σ 2 noise I)

Here, σ noise is the magnitude of measurement error. While performing regression in a GPR framework we first define a family of functions also called a prior described in subsection II.A.1. The next step involves taking the data and eliminating all the functions in our prior which do not obey the observations, this step gives us the posterior mean and covariance described in subsection II.A.2. Finally, we can further improve our predictions by improving our initial family of functions, this is achieved by maximizing the marginal likelihood subsection II.A.3.

Prior

Any regression algorithm starts with an initial family of functions or hypothesis space. By initially defining a family of functions we indirectly encode initial assumptions or bias (continuity, differentiability etc.) into our regression algorithm. If an algorithm is not able to represent the mapping function f in its hypothesis space we will find the closest according to a certain metric.

f (x) = GP (0, k(x, x , θ)) (4)
Here, the representation of mapping function f equation 4 defines a family of functions also called prior. Without loss of generality, we can assume that the mean function is zero. [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] The covariance function or kernel k(x, x ) defines the shape of functions defined by the prior equation 4. In this paper we will compare two types of kernels Standard Exponential kernels and Matérn kernels.

Standard Exponential Kernel

A popular choice of covariance function is a Squared Exponential (SE) function equation 5, because it defines a family of highly smooth (infinitely differentiable) non-linear functions as shown in figure 1a.

k SE (x, x , θ) = θ 2 1 exp[- d 2 2θ 2 2 ] (5)
For the case of the SE kernel equation 5 the hyperparameters (θ = [θ 1 , θ 2 ]) are; amplitude θ 1 which defines average distance from mean and the length scale θ 2 which defines the wiggliness of functions. Here, d defines the absolute distance between points |x -x |.

Matérn Kernel

The Mateŕn kernel is the second most popular kernel after the squared exponential kernel. It provides the flexibility to define a family of functions with varying degree of differentiability. We will take the case of twice differentiable Mateŕn kernel throughout this article.

k M at2 (x, x , θ) = θ 2 1 1 + √ 5d θ 2 + 5d 2 3θ 2 2 exp - √ 5d θ 2 (6)
The hyperparameters (θ = [θ 1 , θ 2 ]) are; the amplitude θ 1 which defines average distance from mean and the length scale θ 2 which defines the wiggliness of functions. The Matérn kernel in equation 6 defines a family of functions differentiable only twice figure 1b. We can observe that draws from figure 1a are smoother when compared to draws from figure 1b for the same set of hyperparameters θ.

Posterior

Once we have defined an appropriate prior, we have a look at the data {(x 1 , y 1 )}. Using the Bayes rule we eliminate all the functions that are not passing through the observed data points. The predicted mean equation 7 and variance equation 8 can be analytically calculated [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] . Figure 2a shows the posterior distribution upon introduction of data to the prior distribution figure 1a. The mean function does not pass through all the observation points because of the assumed measured noise term σ noise = 0.05. Figure 2b shows the posterior distribution upon introduction of data to the prior distribution 1b. The mean function does not pass through all the observation points because of the initial noise term σ noise = 0.05. We can observe that the mean of figure 2a is smoother when compared to the mean of figure 2b. Moreover the variance estimates in the case of the Matérn kernel are higher than that of the SE kernel. Since twice differentiability is a less strict initial assumption and has more functions in its family, more functions are possible through the data.

m(y * ) = k x * x (K xx + σ 2 noise I) -1 y (7) Cov(y * ) = k x * x * -k x * x (K xx + σ 2 noise I) -1 k xx * (8)

Maximizing Marginal Likelihood

To obtain good prediction capabilities, one needs to start from a prior that explains well our dataset. This is achieved by optimizing the Marginal Likelihood (ML) equation 9. The probability that our dataset comes from a family of functions defined by the prior is called the ML. [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] Hence, when we maximize the ML we are actually finding the optimal θ or family of functions that best describe our dataset.

log(P(y | x, θ)) = - 1 2 y T [K xx + σ 2 noise I] -1 y -log K xx + σ 2 noise I - n 2 log(2π) ( 9
)
The ML is a trade-off between a data-fit term ( We observe that the mean of figure 3b passes through all the observed data points but is more complex (wiggly). The mean in figure 3a is a smooth function but does not fit the data properly. While figure 3c is a good compromise between the data-fit term and the model complexity term.

B. Multi-Output Gaussian Process Regression

Given a dataset for multiple outputs {(x i , y i )} for i ∈ [1; D] we define the joint output vector Y = [y 1 ; y 2 ; y 3 ; . . . ; y D ] such that all the output values are stacked one after the other. Similarly, we define the joint input matrix as X = [x 1 ; x 2 ; x 3 ; . . . ; x D ]. For the sake of simplicity, suppose we measure two outputs y 1 and y 2 with some error, while the true physical process is defined by latent variables f 1 and f 2 . Then the relation between the output function, measurement error and true physical process can be written as follows.

y 1 = f 1 + n1 y 2 = f 2 + n2 (10)
Let us take the case of an explicit relationship between the two latent variables equation 11.

f 1 = g (f 2 , x 1 ) (11)
Here g (.) ∈ C 2 is an operator defining the relation between f 1 and an independent latent variable f 2 . Operator g(.) can be a known physical equation or a computer code between the outputs. A GP prior in such a setting with 2 output variables is expressed in equation 12.

f 1 f 2 ∼ GP g(0) 0 , K 11 K 12 K 21 K 22 (12)
K 12 and K 21 are cross-covariances between the two inputs x 1 and x 2 . K 22 is the auto-covariance function of independent output, while K 11 is the auto-covariance of the dependent output variable. The full covariance matrix K XX is also called the joint-covariance. While, the joint error matrix will be denoted by Σ;

Σ = σ 2 n1 0 0 σ 2 n2 ( 13
)
Where, n1 and n2 are measurement error sampled from a white gaussian noise N (0, σ n1 ) and N (0, σ n2 ).

Multi-Output Prior

For a linear operator g(.) the joint-covariance matrix can be derived analytically, [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF] due to the affine property of Gaussian's.

f 1 f 2 ∼ GP g(0) 0 , g(g(K 22 , x 2 ), x 1 ) g(K 22 , x 1 ) g(K 22 , x 2 ) K 22 (14)
Using the known relation between outputs we have successfully correlated two GP priors from equation 11. This effectively means that when we randomly draw a function f 2 it will result in a correlated draw of f 1 such that the two draws satisfy the equation 11. We have effectively represented the covariance function K 11 in terms of the hyperparameters of covariance function K 22 using the known relation between outputs. However, a non-linear operation g(.) on a Gaussian distributed variable does not result in a Gaussian distributed variable. Hence, for the case of non-linear g (.) the above joint-covariance matrix as derived in Equation 14is not positive semi-definite, which prevents us from taking its inverse in equation 9. Therefore we use an approximate joint-covariance [START_REF] Constantinescu | Physics-Based Covariance Models for Gaussian Processes with Multiple Outputs[END_REF] for imposing non-linear relations equation 15.

f 1 f 2 ∼ GP g(0) 0 , LK 22 L T LK 22 K 22 L T K 22 (15)
Here, L = ∂g ∂y y 2 = ȳ2

is the Jacobian matrix of g (.) evaluated at the mean of independent output y 2 .

Equation 15 is a Taylor series expansion for approximating related kernels. Since a Taylor series expansion is constructed from derivatives of a function which are linear operations, the resulting approximated joint kernel is a Gaussian kernel. [START_REF] Constantinescu | Physics-Based Covariance Models for Gaussian Processes with Multiple Outputs[END_REF] Figure 4 shows the joint prior for two functions f 1 and f 2 given a relation f 1 = df 2 /dx between them. The pairs of similarly colored random draws follow the given equation. We can verify this fact by observing that f 1 tends to zero when f 2 tends to a maxima or minima. The draws from figure 4a are smoother when compared to draws from figure 4b.

Multi-Output Posterior

To obtain predictions on non-measured points (X * ), we now proceed in a similar manner as prescribed in section A. We optimize the ML to find the optimal θ and then find the predicted mean and covariance as shown in equation 16 and equation 17.

m(Y * ) = K X * X (K XX + Σ) -1 Y (16) Cov(Y * ) = K X * X * -K X * X (K XX + Σ) -1 K XX * ( 17
)
As mentioned earlier while calculating the posterior we apply the bayes rule. This effectively means that we eliminate all the functions that do not pass through measured points. Since we have encoded the relation between f 1 and f 2 in the joint-prior. When we eliminate not possible functions from the prior of f 1 , we consequently eliminate functions from the prior of f 2 and vice-versa. In the process we come up with a posterior which is consistent with the observations and initial relation. Figure 5 shows the joint posterior for two functions f 1 and f 2 given a relation f 1 = df 2 /dx and data between them. We have used same data as used in figure 2 but have reduced the number of measurement to 8 and derivative points to 8. Figure 5a shows the posterior for a SE kernel. At the positions where the value of derivative f 1 is given, we see a shrinking of variance in f 2 . Draws from the posterior and mean of f 2 both have value of derivative as set by f 1 . Figure 5b shows the posterior for SE kernel. At the positions where value of derivative f 1 is given, we see a shrinking of variance in f 2 . Draws from the posterior and mean of f 2 both have value of derivative as set by f 1 .

Here the two outputs need not be measured at the same input points, this will be demonstrated in subsection III.B. Another interesting observation is that while optimizing the marginal likelihood for figure 5, the noise terms for two outputs are different (σ n2 = 0.02; σ n1 = 0.1). This means that the two outputs have different values of fidelity, this looks very similar to co-kriging. In fact co-kriging can be seen as a special case of joint multi-output GPR. Such that the two outputs are linked by the relation f 1 = f 2 and the two error measurements are linked by the relation σ n1 >> σ n2 , the output y 1 is a low fidelity measurement when compared to output y 2 .

III. Results

In this section, we provide a numerical illustration to the theoretical derivations in the earlier sections. We start with a synthetic problem where we try to learn the model over quadratic relationships. We compare the cross-validation error values of independent GPR with that of multi-output joint GPR. Finally, we compare the performance of our methods on flight-loads estimation for a horizontal tail plane.

The basic toolbox used for this paper is GPML provided with "Gaussian Process for Machine Learning", 4 we generate covariance functions to handle relationships as described in equations 14 and 15 using the "Symbolic Math Toolbox" in MATLAB 2014b. All experiments were performed on an Intel quad-core processor with 4Gb RAM.

A. Quadratic relation on Synthetic Data

We take the case of a quadratic operator g(.) equation 18.

f 1 = f 2 2 (18)
To generate the data we randomly draw a single function of f 2 as described in equation 19 for 50 equally spaced inputs between [-1, 1]. The data for f 1 is the calculated using equation 18. The latent functions f are then corrupted according to equation 19 which gives us the outputs y s. We finally hide the observations for y 1 in the domain x = [-0.2, 0.2]. We thus have our dataset for testing the performance of quadratic relationship.

f 2 ∼ GP [0, K SE (0.2, 1)] σ n2 ∼ N [0, 0.1] σ n1 ∼ N [0, 1] (19)
K SE (0.2, 1) means squared exponential kernel with length scale 0.2 and variance 1. σ n2 and σ n1 are the white noises added to the latent functions f 1 and f 2 respectively. Since the quadratic relationship g(.) is non-linear in nature we use the equation 15 to calculate the auto-and cross-covariance functions as explained in equation 20.

K 12 = 2µ y2 (x 1 )K 22 K 11 = 4µ y2 (x 1 ) 2 K 22 (20)
Here, µ y2 (x 1 ) is the mean value of function y 2 calculated at the input points x 1 . For the case of quadratic relationship the jacobian L as described in equation 15 comes out to be 2µ y2 (x 1 ). For non-linear operators g(.) the joint-covariance prior depends on the mean value of f 2 . K 22 and K 11 are the auto-covariance functions calculated at the input points x 1 and x 2 as described in equation 15 and equation 12. K 12 is the cross-covariance function between the input points x 1 and x 2 . We can observe that the joint-covariance has been expressed as a function of covariance K 22 using equation 18.

Figure 6a shows the independent fit of two GPR. For the case of y 1 there is a huge difference between the real data and predicted mean at points where data is unavailable. Figure 6b shows the joint-GP regression. Joint-GP model gives better prediction even in the absence of data for y 1 because transfer of information is happening from observations of y 2 present at those locations. For the case of y 2 there is no significant improvement in the prediction of two methods. 1 shows comparison of Root Mean Square Error (RMSE). 10 sets of experiments were run for 85% of data as training set and 15% of data as test set, the training and test sets were chosen randomly. We learn the optimal set of hyperparameters for on training data all 10 sets of experiments. Finally RMSE values are evaluated with respect to the test set. We see a significant improvement in performance for the case of joint multi-output GPR.

B. Flight Mechanics on Flight Test Data

Procuring data in engineering design is a costly exercise, a high-fidelity CFD simulation runs for days and a flight test campaign costs millions. In this context it is important to be data efficient while creating models with data. We leverage the knowledge of available relations between several types of data and build a more robust, physically consistent and data efficient model.

In this section we conduct experiments, applying our approach on the flight loads data. We look at normalized data of a simple longitudinal maneuver. The maneuver is quasi-static which means that airplane is in equilibrium at all times and there are no dynamic effects observed by the aircraft. The two outputs in our case are coefficient of lift C z on the Horizontal Tail Plane (HTP) and spanwise lift k z . We know that the integral of spanwise pressure will be equal to the coefficient of lift.

C z (η) = ηroot η edge k Z (η)dη (21)
.

Here, η edge denotes the edge of horizontal tail span while η root represents root of tail. The above equation is linear in nature and hence we will use equation 14 to calculate the auto-and cross-covariance functions.

In practice, the coefficient of lift C z is calculated using strain gauges on the HTP and the spanwise force is measured using pressure tappings. Chordwise pressures are measured at multiple η locations on the HTP, which are integrated chordwise to give spanwise lift distribution. We are trying to merge data coming from two different calculation processes of flight mechanics and aerodynamics using a basic physics based relation and improve the accuracy of our model. Figure 7 shows a comparison between two joint multi-output GPR on K z and C z using initial assumptions regarding order of differentiability. Figure 7a shows the effect of SE kernel (infinite differentiability) on the joint model. The SE kernel is a strong assumption for type of lift distribution. To accomodate for the value of measured C z the error band increases and the boundary condition at HTP root is not satisfied. Moreover, we identify points in the figure 7a which fall out of the 95% confidence band. We can claim with 95% probability that if our lift distributions come from a family of highly differentiable functions these outliers do not satisfy the physics of the problem and can be termed as measurement errors.

Figure 7b shows the effect of twice differentiable Matérn kernel on the joint model. The twice differentiable assumption gives our family of functions enough flexibility to follow the spanwise lift data. We observe that error estimates for C z improves considerably while there is an increased uncertainity at the HTP root for K z . This loss in order of differentiability can be explained due to fuselage HTP interaction near HTP root. Since there is a discontinuity in the surface due to presence of fuselage, the pressure distribution near the HTP root cannot remain infinitely smooth. We can conclude that even though we are working with very basic assumptions of differentiability the choice of kernel is important in performing the predictions.

This procured model can be used as a robust surrogate model which is consistent with the physics based relations. The relation g(.) can be a known functional relationship or some transformation arising from a computer code. This surrogate model can be used in further use cases for example finding optimal placements of sensors so as to reduce the overall uncertainty in the system, or performing bayesian optimization or Efficient Global Optimization (EGO) using information coming from several outputs. Making the steps more data efficient and consistent with the physics of the system.

IV. Conclusion

This paper describes a general framework of joint multi-output GPR to find mapping function between a set of input and related outputs [START_REF] Constantinescu | Physics-Based Covariance Models for Gaussian Processes with Multiple Outputs[END_REF] . Using very basic assumptions such as order of differentiability and knowledge of relation between the outputs a joint family of functions can be constructed. Later functions inconsistent with the observed data are eliminated using bayes rule. This gives us a prediction which is consistent with the initial assumptions, initial relationship and observed data.

We initially provide empirical demonstrations of the framework on synthetic data related through a quadratic relationship. The approach shows improvements in prediction capabilities and handling of missing data. The second demonstration was performed on flight loads data and demonstrated that imparting prior flight mechanical relationships into GPR allows for a better surrogate flight loads model. Thus, we are now capable of creating a more robust, physically consistent and interpretable surrogate model for these loads, thereby greatly improving the identification process.

As we keep on adding more physical information into GPR, the number of points in the joint-covariance matrix increase. However GPR does not scale well to very large number of input points. Future work will investigate recent advances in scaling up GPR 14 in order to handle several relationships and possibly cover the whole flight test dataset.

  (a) Draws from a Standard Exponential (SE) Kernel (θ = [1, 0.2]). (b) Draws from a Mateŕn Kernel differentiable twice (θ = [1, 0.2])

Figure 1 :

 1 Figure 1: Draws from Gaussian Process Priors. Solid black line defines the mean function, shaded blue region defines 95% confidence interval (2σ) distance away from mean. The dashed lines are functions randomly drawn from the priors

  (a) Posterior from a Standard Exponential (SE) Kernel (θ = [1, 0.2]; σnoise = 0.05). (b) Posterior from a Mateŕn Kernel differentiable twice (θ = [1, 0.2]; σnoise = 0.05).

Figure 2 :

 2 Figure 2: Gaussian Process Posterior. Solid black line defines the mean function, shaded blue region defines 95% confidence interval (2σ) distance away from mean. The dashed lines are functions randomly drawn from the posterior

  1 2 y T K -1 xx y) and a model complexity term (log |K xx |). The optimization of ML(θ) provides the best compromise in terms of explaining the existing dataset {(x 1 , y 1 )} and the initial assumptions encoded in the prior. Hence to find the best set of hyperparameters we maximize the marginal likelihood equation 9 with respect to its hyperparameters. log(ML) = -35.3 [0.05, 0.5]; σnoise = 0.01) and data. log(ML) = -8.2 [0.05, 0.5]; σnoise = 0.01) and data. log(ML) = 6.04 (c) Posterior between SE prior with optimized hyperparameters (θ = [0.45, 0.17]; σnoise = 0.0014) and data.

Figure 3 :

 3 Figure 3: Posteriors for 3 different sets of hyperparameters

  (a) Prior from a joint SE Kernel such that (θ = [1, 0.2]) and f1 = df2/dx. (b) Prior from a joint Mateŕn Kernel differentiable twice such that (θ = [1, 0.2]) and f1 = df2/dx.

Figure 4 :

 4 Figure 4: Multi-Output Gaussian Process Prior. Solid black line defines the mean function, shaded regions define 95% confidence interval (2σ) distance away from mean. The dashed lines are functions randomly drawn from the prior. Pairs of similar colored dashed lines follow the given relation between f 1 and f 2 .

  (a) Posterior from a joint SE Kernel such that (θ = [0.45, 0.17]; σn2 = 0.02; σn1 = 0.1) and f1 = df2/dx. (b) Posterior from a joint Mateŕn Kernel differentiable twice such that (θ = [0.45, 0.17]; σn2 = 0.02; σn1 = 0.1 and f1 = df2/dx.

Figure 5 :

 5 Figure 5: Multi-Output Gaussian Process Posterior. Solid black line defines the mean function, shaded region defines 95% confidence interval (2σ) distance away from mean. The dashed lines are functions randomly drawn from the posterior

  (a) Independent GP Regression for the two outputs y1 and y2. For y1 the data is hidden from section x = [-0.2, 0.2]. We can observe the huge difference between the real data and the predicted mean values at zone with no data.(b) Joint-GP Regression for the two outputs y1 and y2 related through equation 18. For y1 the data is hidden from section x = [-0.2, 0.2]. We can observe the improved prediction between zone with no data because information is being shared between the two outputs.

Figure 6 :

 6 Figure 6: GP Regression on Quadratic relationship. The solid black line represents the predicted mean while the shaded area denotes 95 % confidence interval (2σ) uncertainity region. The dashed black line represents the real value of f 1 . For y 1 the data is hidden from section x = [-0.2, 0.2].

  (a) Joint multi-output GPR between Kz and Cz using a SE kernel and equation 21. Optimized hyperparameters (θ = [0.42, 0.08]; σnKz = 1.3 * 10 -4 ; σnCz = 0.0414) . (b) Joint multi-output GPR between Kz and Cz using a twice differentiable Matérn kernel and equation 21. Optimized hyperparameters (θ = [0.53, 0.08]; σnKz = 9.3 * 10 -7 ; σnCz = 0.018)

Figure 7 :

 7 Figure 7: Joint multi-output GPR on K z and C z relationship. Solid black line defines the mean function, shaded region defines 95% confidence interval (2σ) distance away from mean.

  Adding Flight Mechanics to Flight Loads Surrogate Model using Multi-Output Gaussian Processes Ankit Chiplunkar * and Michele Colombo † Université de Toulouse, Institut Clement Ader (ICA), CNRS, ISAE-SUPAERO, Toulouse, 31077, France

	Airbus Operations S.A.S., Toulouse, 31060, France
	Emmanuel Rachelson ‡
	ISAE-Supaero, Département d'Ingénierie des Systèmes Complexes (DISC), Toulouse, 31055, France
	Joseph Morlier §

Table 1 :

 1 mean RMSE errors for quadratic relationship RMSE y 1 RMSE y 2

	Independent Single-output GPR	1.74	0.14
	Joint Multi-output GPR	0.26	0.12
	Table


Laurent Demagistri

Danielle Mitja

Eric Delaître

Ehsan Yazdanparast

Hamidreza Shahbazkia

Michel Petit

Palm-trees detection with very high resolution images, comparison between Geoeye and Pléiades sensors

come

•

  Ce travail bénéficie du soutien du programme TOSCA du CNES, incluant des images Pléiades. Les auteurs remercient les agriculteurs de la communauté de Benfica pour leur aide, disponibilité et convivialité durant les phases de travail sur le terrain.ContexteLa gestion des forêts tropicales passe par la conservation de la biodiversité notamment grâce au maintien des arbres et des palmiers utiles dans les agrosystèmes. Cependant ces espèces sont susceptibles d'être surexploitées et de disparaître. Leur surveillance nécessite d'identifier les zones d'occurrence, de caractériser les impacts des actions anthropiques et d'en assurer le suivi, par des indicateurs idoines qui font de plus en plus appel à la télédétection. L'étude porte sur un projet d'installation de petits agriculteurs sur un front pionnier amazonien : le PA-Benfica d'une superficie de 10000 ha, localisé dans la commune d'Itupiranga, dans le sud-est de l'état du Pará. Cette étude est financée dans le cadre du Projet CNES-TOSCA : CIC-TOOB -Chaîne d'Intégration des Capteurs de Télédétection pour le suivi de l'Occupation du sol et l'extraction automatique d'objets ressources : le cas du palmier Babaçu en Amazonie (CONVENTION CNES /IRD N° 130318). Le babaçu est un palmier natif dans ces forêts tropicales. Il peut atteindre 30m de haut et le diamètre moyen de sa couronne est de 8m. Lors du défrichement de la forêt pour l'installation des cultures et des pâturages, de nombreux pieds de cette espèce sont conservés et résistent de nombreuses années. Dans les pâturages, les couronnes, leurs ombres et leurs troncs sont visibles sur les images satellites THR et il est possible de marquer chaque individu. Cependant même s'il représente 93,6 % des palmiers des pâturages, le babaçu est en mélange avec d'autres espèces ayant une couronne de dimension comparable. Un repérage de la couronne des palmiers est tout d'abord effectué par photo-interprétation par un expert, sur une image Geoeye (résolution : 0,5m, fusion XS + P) de 2009 suivant un protocole d'échantillonnage particulier découpant l'image en mailles régulières d'1 ha (sur les 17526 mailles, 2354 ont été échantillonnées). Seules les mailles 100% pâturage sont considérées dans cette étude, soit 1075 mailles. Elles constituent une base de palmiers photo-interprétés sur les images THR qui permettra de généraliser l'estimation des performances des algorithmes de détection d'objets. Un algorithme de reconnaissance de palmiers sur des images THR panchromatiques a été mis au point en milieu ouvert (pâturages). Des techniques de filtrage et de morphologie mathématique sont utilisées pour extraire des images les pixels correspondant à des ombres. Le modèle de classification des ombres utilise des données d'apprentissage : un jeu d'ombres de palmiers et un jeu d'autres ombres non-palmiers. 25 attributs de forme sont calculés pour construire ce modèle. 2 L'algorithme de reconnaissance du palmier dans les pâturages et l'ajustement du modèle Extraction des ombres • Lissage • Filtre adaptatif local • Filtrage morphologique (fermeture) • Extraction des composantes connexes • Seuillage sur la taille des objets Classification des ombres • Apprentissage à partir d'un jeu d'entrainement • Reconnaissance des palmiers à partir de est en cours de test sur des images Pléiades et Geoeye prises en juillet 2013 à 1 semaine d'intervalle. La différence de résolution intrinsèque des capteurs apparait nettement : contour des objets et « grain » de l'image Geoeye plus précis. Les derniers résultats obtenus après ajustement de l'algorithme de détection sur une image Geoeye 2009 donnent un taux de 75,5% de palmiers correctement détectés. (Taux de détection visé : 80%) Légende: objets détectés automatiquement vrai positif (1 | 1) faux négatif (1 | 0) faux positif (0 | 1) : 0,7  0,5 m Date : 2013-07-19 13:54 Angle d'Incidence : 15.75° Élévation solaire : 55.59° Geoeye 2013 Résolution : 0,4  0,5m Date : 2013-07-12 13:39 Angle d'Incidence : 10,96° Élévation solaire : 52,58° palmier Vu le nombre de mailles photo-interprétées, l'estimation des taux de détection sur la totalité de la zone d'étude implique d'automatiser la mise en correspondance des séries de points issus (1) de la photo-interprétation (repérage des couronnes) et (2) de l'algorithme de détection (sur les ombres). . Les points des deux séries sont appariés en fixant la direction () et le rayon (r) de recherche en fonction de la position moyenne des ombres dans l'image. Pour chaque point détecté (« ombre »), les candidats potentiels de la série de points « couronne » sont triés à l'aide d'une fonction de poids paramétrée par l'orientation du couple de points considérés et leur distance. Pour la détection des palmiers les images THR (résolution inférieure à 1m) sont bien adaptées, car elles permettent de visualiser les palmes. • Dans les pâturages la photo-interprétation passe par une visualisation de l'ensemble couronnes-ombre de couronne alors que l'algorithme de détection automatique est appliqué à l'ombre de la couronne • La validation automatique permet d'apparier les palmiers détectés aux palmiers issus de la photo-interprétation. • Le taux de détection visé est de 80%. Pour Geoeye, les résultats sont supérieurs à 75%. Le calcul pour Pléiades n'est pas finalisé. L'étude en cours prévoit : • Un nouvel algorithme pour la détection des palmiers en milieu fermé (forêts/jachères) • Une tentative de différenciation des espèces de palmiers. Conclusions et perspectives Palm-trees detection with very high resolution images, comparison between Geoeye and Pléiades sensors Laurent DEMAGISTRI, Danielle MITJA, Eric DELAITRE, Ehsan YAZDANPARAST, Hamidreza SHAHBAZKIA, Michel PETIT

Acknowledgments

This work benefits from the support of the CNES TOSCA program, including Pléiades images. The authors thanks the farmers of the Benfica-PA community for their help, availability and conviviality during on-site work.

Framework

Managing rainforests implies preservation of biodiversity particularly by maintaining, within agro-ecosystems, trees and palm-trees that are useful for the local people. These species are likely to be overexploited and to disappear. The monitoring of their extent, implies: identify areas of occurrence and characterize on these, the impacts of human activities. For this purpose, nowadays, Remote Sensing techniques play an important part to set up suitable indicators. This study focuses on a Brazilian project dedicated to the installation of small farmers in some places of Amazonian deforestation areas: the "PA-Benfica" (10 000 ha), in the Itupiranga municipality (southeast Pará). This study is funded by CNES-TOSCA Project "CIC-TOOB": « Chaîne d'Intégration des Capteurs de Télédétection pour le suivi de l'Occupation du sol et l'extraction automatique d'objets ressources : le cas du palmier Babaçu en Amazonie » (CNES/IRD agreement No. 130318).

Babassu is an endemic palm-tree in those rainforests. It can reach 30m high and its crown mean diameter is 8m. When clearing forests for new crops and pastures, lots of stems of this species are preserved and resist many years. Crowns of palm-trees, their trunks and shadows, in pastures areas, can be seen on very high resolution satellite images; it is possible to mark each single palmtree. However, even if it counts for 93.6% of palm-trees in pastures, the babassu is mixed with other species with similar crown size. Firstly, the crowns of palm-trees are detected and marked manually by an expert (photo-interpretation) on a GeoEye image (spatial resolution: 0.5m, merged XS+P). A particular sampling protocol is used: it is based on a regular grid with 1ha cells. Out of the 17526 cells, 2354 cells have been analyzed by the expert. For this study, only "100% pasture" cells are considered (1075). This data base of marked palm-trees by photo-interpretation of a VHR image will be used to verify and extrapolate the detection rate of the object recognition algorithms. Given the great number of analyzed cells for validation, estimating the good-detection rate for the whole studied area implies automating the matching between (1) manually marked palm-trees (the crowns) and (2) detected targets (the shadows).

Palm-trees and their detection by photo-interpretation

. The two sets of points are mapped together by setting a search direction () and a radius (r). Both parameters depends on the mean location of the shadows in the image.

For each detected point (shadow), the potential targets from the crown set are sorted using a weight function of the distance and the angle between the points.

Automatic validation

Center of gravity of the shadows ( 2

Conclusions and work in prospect