Yann Denieul

Daniel Alazard

Joel Bordeneuve

Clement Toussaint

Gilles Taquin

V Aerodynamic

Interactions of Aircraft Design and Control: Actuators Sizing and Optimization for an Unstable Blended Wing-Body

Keywords:

In this paper the problem of integrated design and control for a civil blended wingbody aircraft is addressed. Indeed this configuration faces remarkable challenges related to handling qualities: namely the aircraft configuration in this study features a strong longitudinal instability for some specific flight points. Moreover it may lack control efficiency despite large and redundant movables. Stabilizing such a configuration may then lead to high control surfaces rates, meaning significant energy penalty and installation mass for flight control actuators, as well as challenges for actuators space allocation. Those penalties should therefore be taken into account early in the conceptual design phase, instead of being checked afterwards. Our approach consists in simultaneously designing a stabilizing controller and actuators dynamic characteristics, namely their bandwidth, for a given aircraft configuration. Our method relies on latest developments on nonsmooth optimization techniques for robust control design. For any aircraft configuration, guaranteed stability performance, as well as optimized control surfaces allocation with respect to their relative efficiency and inertia, are obtained. From these results, a segregation among trim and maneuver control surfaces is obtained, with guarantee that the later ones are able to cope with aircraft instability. Then it is checked that remaining trim control surfaces are sufficient for equilibrating the aircraft in any flight condition. This approach allows for a fast prototyping of control surfaces and actuators for unstable configurations. Different control surfaces layouts are evaluated in order to show the flexibility of the method. Nomenclature α Angle of attack C D Drag coefficient C m Pitching moment Cm 0 Zero-lift pitching moment Cm q Pitching moment pitch rate gradient C L Lift coefficient C Lq Lift coefficient pitch rate gradient δm Control surface deflection

I. Introduction

Among other disruptive aircraft configurations, the Flying Wing (also known in litterature as Blended Wing-Body (BWB) [START_REF] Roman | Aerodynamic Design Challenges of the Blended-Wing-Body Subsonic Transport[END_REF] or more recently Hybrid Wing-Body 2 ) has been identified for years as a potential candidate for the future of civil aviation. [START_REF] Liebeck | Design of the Blended-Wing Body subsonic transport[END_REF] Most studies in this field focused either on the overall aircraft design problem, optimizing the aircraft with respect to "traditional" disciplines in the conceptual design phase, such as aerodynamics, propulsion, weights and performance 1 , 4 or on stability and control related challenges for a given configuration. [START_REF] Kozek | Modeling and Control for a Blended Wing Body Aircraft[END_REF] However this sequential approach, while valid for traditional aircraft configurations, appears as too restrictive for the BWB case. Indeed its unusual characteristics, such as a strong longitudinal instability, do not allow treating the handling qualities issues without considering active stabilisation in the early design phase 6 . [START_REF] Wildschek | Flight Dynamics and Control Related Challenges for Design of a Commercial Blended Wing Body Aircraft[END_REF] At Airbus previous research topics were focusing on plan form optimization for performance objectives on the one hand, and handling qualities resolution [START_REF] Saucez | Optimization of Engine Failure on a Flying Wing Configuration[END_REF] on the other hand. This meant checking regulation and in-house criteria for maneuvers, in flight and on ground equilibria. One of the conclusions of these studies was the need for active stabilisation. It was also concluded that the active stabilisation may induce high deflection rates on control surfaces; this is even more problematic as these surfaces may be significantly larger than usual movables. They may also lack pitch efficiency. As a result a major mass and energy penalty would result from flight control actuators required to stabilize the aircraft. What we propose to perform here is a synthesis of these couplings between aircraft design and control, as well as setting a process capable of sizing and optimizing the flight control actuators in a preliminary way for dynamic stability criteria. This paper is organized as follows: in section II the general problem of stabilizing an aircraft with redundant actuators is adressed, the aerodynamic model and the different control surfaces layouts are presented. Then section III presents the problem of integrated design and control, and the specific H 2 /H ∞ formulation used in this paper in order to simultaneously optimize the controller and the actuators bandwidth. The objective is to meet constraints on handling qualities requirements while minimizing the "cost" of actuators. Finally section IV discusses the results obtained for different control surfaces layouts in terms of achievable actuators bandwidth and remaining trim capability.

II. Problem Setup

II.A. Design Challenges associated to the Control of an Unstable Configuration

From a handling qualities point of view, an aircraft is said statically stable if its center of gravity (CG) [START_REF] Roskam | Airplane flight dynamics and automatic flight controls[END_REF] is located forwards its aerodynamic center a . The aerodynamic center of an aircraft is mainly driven by its lifting surfaces characteristics, ie wing location with respect to the fuselage, airfoils and twist. Whereas static stability, or at least neutral stability, used to be a design constraint at preliminary design stage, recent aircraft designs tend towards instability. From a conceptual design point of view, Relaxed Static Stability (RSS) may be seeked for different reasons:

• High levels of instability enable excellent maneuverability properties for fighter aircraft. For instance the X-36 tailless demonstrator [START_REF] Brinker | Flight Testing of Reconfigurable Control Law on the X-36 Tailless Aircraft[END_REF] or the Grumman X-29 forward-swept wing experimental aircraft (see figure 1(a)), are both designed with negative static stability margins (35%l for the X-29 11 ); as a result any perturbation from an initial longitudinal equilibrium makes the aircraft dynamically depart. If this motion is adequately controlled, superior longitudinal maneuverability is obtained.

• For civil commercial aircraft, RSS may be a desirable feature for cruise performance optimization. Indeed a conventional statically stable tail-aft configuration requires a downlift force on the aft horizontal stabilizer in order to trim the aircraft in cruise. This results in higher lift coefficient on the main wing, therefore in higher induced and wave drag . Reducing the static margin by shifting the CG backwards and diminishing the horizontal tailplane size is therefore beneficial from a lift-over-drag ratio (LoD) point of view b : the induced drag and wetted area are decreased (see for instance the A380 figure 1(b)).

Of course levels of instability for such configurations are not in the same order of magnitude as for super-maneuverable fighter aircraft.

• However for some configurations instability is no more a design requirement, but a consequence of an overall layout. On the BWB studied in this paper (see figure 1(c)), planform and airfoil design results of optimisation of the LoD with constraints on the zero-lift pitching moment Cm 0 to ensure a feasible take-off rotation despite the nose-down pitching moment of the high-mounted engines. [START_REF] Meheut | Aerodynamic Shape Optimizations of a Blended Wing Body Configuration for Several Wing Planforms[END_REF] The aerodynamic center then comes as an output of the optimization process. The specific lift distribution of the BWB, featuring a large lifting centerbody, implies generally a large lift produced at the forward part of the aircraft. This tends to move the aerodynamic center forward. On the contrary fuel located in the wing box tends to move the CG backwards, resulting in a negative static margin and a strong instability for some flight points (see on figure 2 the relative position of the CG and aerodynamic center for different masses and speeds). It should be noted however that instability is not inherent to the BWB configuration. In the work of Lyu and Martins for instance [START_REF] Lyu | Aerodynamic Shape Optimization of a Blended-Wing-Body Aircraft[END_REF] , [START_REF] Lyu | Aerodynamic Design Optimization Studies of a Blended-Wing-Body Aircraft[END_REF] an optimization is performed on the planform and profiles of a BWB with constraints on cruise trim and minimum static margin. The resulting geometry is optimal from a LoD point of view within the constraints -therefore the optimal LoD is slightly decreased with respect to optimal LoD without stability constraint-, and remains stable and equilibrated.

A question then arises: for a given configuration, are we able to adequately control its unstable modes, and at which expense on the flight control system (FCS) cost and complexity? Indeed it was shown in previous a From a dynamic point of view, the appropriate point to be considered for stability is the maneuver point. If the CG is at the maneuver point position, then the short period mode is at limit of stability. This point lies backwards the aerodynamic center, due to Cmq and C Lq damping effects. On a conventional tail-aft configuration, Cmq effect can be significant due to tail lever arm. On the contrary on a BWB Cmq effect is very small, meaning that aerodynamic center and maneuver point are nearly equivalent.

b More precisely trim drag is minimal for the CG located at a so-called Center of Pressure -CoP -. For a CG located at the CoP no pitching moment is needed to balance the aircraft. The aircraft CoP lies backwards the aerodynamic center: therefore going towards this point tends to increase the aircraft instability. studies that despite constant advances in the field of Control, fundamental physical limitations remain when controlling an unstable plant. A seminal lecture by Stein [START_REF] Stein | The practical, physical (and sometimes dangerous) consequences of control must be respected, and the underlying principles must be clearly and well taught[END_REF] emphasises the fact that a controlled unstable mode remains only locally stable, and turns unstable again when reaching non-linearities such as actuators saturations and rate limits. Whatever the control method used, physical limitations such as actuators bandwidth and rate limits, sensors and processors sampling rates, mechanical structure modes prevent from controlling high frequency unstable modes. A consequence of this is illustrated in a paper by Rogers and Collins 11 about X-29 control synthesis. Two control methods, namely H 2 and H ∞ design are compared. From a theoretical point of view the H ∞ design exhibits superior performance, however actuators limitations are largely exceeded. When taking into account those limitations into the control design, both methods perform equally, but neither is able to meet the requirements. The aim of this study is to determine the required actuators dynamics that stabilize our configuration with appropriate performance, by taking advantage of the redundant actuators on the whole trailing edge, as can be seen on figure 3(a). Actuators sizing for BWB large control surfaces with possibly high deflection rates has been identified for years as a challenging task. [START_REF] Cameron | Control Allocation Challenges and Requirements for the Blended Wing Body[END_REF] Quoting a seminal paper on BWB configuration by Liebeck, 3 "If the BWB is designed with negative static margin (unstable), it will require active flight control with a high bandwidth, and the control system power required may be prohibitive". Indeed as stated by Garmendia et al., [START_REF] Daniel | Assessment of Electrically Actuated Redundant Control Surface Layouts for a Hybrid Wing Body Concept[END_REF] secondary power for FCS P F CS may be evaluated in a preliminary way by the equation 1:

P F CS = n controls i=1 HM max i . θmax i (1)
where HM max i and θmax i are the maximum hinge moment and maximum deflection rate of the i-th control surface respectively, and n controls is the number of control surfaces.

• On the one hand, from Roman et al. [START_REF] Roman | Aerodynamic Design Challenges of the Blended-Wing-Body Subsonic Transport[END_REF] we know that hinge moments are related to the scale of a control surface through a "square-cube law": control surface area increases as the square of the scale λ 2 , whereas hinge moment increases with the cube of this scale λ 3 . Large BWB control surfaces lead to high hinge moments requirements. However hinge moments computation is out of scope of this study.

• On the other hand, deflection rate is a direct consequence of a stability augmentation system (SAS) for an unstable configuration. It can be demonstrated that under the assumptions of a linear dynamic model with static state-feedback compensator K and actuator modeled by a first-order transfer function with bandwidth ω (for notations please refer to figure 7):

y act u act = ω ω + s (2)
The maximal deflection rate θmax for an initial perturbation X 0 on the aircraft states -e.g. a perturbation on the angle of attack -is obtained at the initial time and is equal to:

θmax = θ0 = ωK(X 0 -X eq ) (3)
This means that deflection rates are both a function of the actuators dynamics and the control law synthesis.

The aim of this study is then to perform a combined synthesis of control laws and actuators bandwidth sizing, in order to limit deflection rates induced by the SAS. Keeping in mind previous remarks on hinge moments, deflection rates should moreover be more limited for larger control surfaces.

II.B. The Blended Wing-Body Configuration

In this section the aerodynamic model as well as the three evaluated control surfaces layouts are described.

II.B.1. Aerodynamic Model

The flight dynamics model used in this study comprises only longitudinal equations written in the aerodynamic reference system R a (x a , y a , z a ):

m V = - 1 2 ρV 2 SC D -mg sin γ + F (4)
-mV γ = - 1 2 ρV 2 SC L + mg cos γ (5)
B q = 1 2 ρV 2 SlC m (6) θ = q (7) θ = α + γ (8)
where:

• concerning aircraft parameters, m and B denote the aircraft mass and inertia around the y a axis respectively, S and l are reference surface and length, corresponding to the wing aera and mean aerodynamic chord respectively.

• angles are classically defined as follows: γ, α and θ denote the flight path angle, angle of attack and aircraft pitch attitude respectively.

• concerning aerodynamic parameters, V denotes the aerodynamic speed, ρ denotes the volumetric mass, and C D , C L and C m denote the drag, lift and pitching moment coefficients respectively.

This model is fully non-linear, with embedded non-linearities in the aerodynamic coefficients C D , C L and C m with respect to variables such as Mach, angle of attack and control surfaces deflections. This model will be used for trim capacity evaluation of section IV.B. Aerodynamic coefficients for this flight dynamics model were obtained using different CFD methods. [START_REF] Meheut | Aerodynamic design of transonic flying wing configurations[END_REF] For control synthesis purpose, a linear state-space representation is obtained from Eq. 4 to 8 using finite differences method on different flight points. This state-space representation reads:

     δ V δ α δ q δ θ      =      x V x α x q x θ -z V -z α 1 -z q 0 0 m α m q 0 0 0 1 0           δV δα δq δθ      +      x δx x δmi 0 -z δmi 0 m δmi 0 0      ∆δx ∆δm i (9)
The states are δV , δα, δq, δθ, variations around an equilibrium of the airspeed, angle of attack, pitch rate, and attitude respectively. The different terms of all matrices are developed in the Appendix. Concerning the controls, ∆δx denotes the thrust command, and ∆δm i denotes the i-th control surface command, the different control surfaces layouts being developped in Section II.B.2.

II.B.2. Control Surfaces Layout

In this section three different control surfaces layouts are presented and their rationale is explained.

The problem of number and spacing of trailing edge control surfaces for BWB was already pointed out by Garmendia et al [START_REF] Daniel | Assessment of Electrically Actuated Redundant Control Surface Layouts for a Hybrid Wing Body Concept[END_REF] . [START_REF] Garmendia | Method for Evaluating Electrically Actuated Hybrid Wing Body Control Surface Layouts[END_REF] They provided an extensive view on all currently studied BWB configurations and controls layouts. Another study of interest was performed by Belschner: [START_REF] Belschner | Definition of systems architecture for flight, load and vibration control of blended wing body type aircraft[END_REF] an optimisation was run on failure cases considering electro-mechanical actuators (EMA) driving independant control surfaces. An optimal layout of 23 control surfaces resulted of the process. Failure case analysis is out of scope of our study, as is a sizing of a whole FCS. Rather we expect to rough out actuator dynamics, which would then serve as an input for a more precise FCS sizing.

The initial configuration evaluated here is visible on figure 3(a). Five control surfaces are spread along the half span, except a gap from elevon 1 to 2 -control surfaces are numbered from inboard to outboard -due to the presence of the engine pylon. The whole trailing edge of a BWB is usually devoted to control, as the lack of longitudinal lever arm should be compensated by large surfaces. On control surfaces 3, 4 and 5 the relative inboard and outboard chords were kept constant, so their absolute chords decrease due to the taper ratio of the outer wing. Their chords are moreover limited by a rear spar. The inboard chord of the second control is also limited by cargo and cabin considerations. (c) Iso-inertia layout. From a static point of view, previous studies6 demonstrated that total control surfaces efficiency of this configuration was sufficient to fulfill both longitudinal and lateral handling qualities requirements, provided all surfaces are used as multicontrol. Control allocation algorithms were implemented for that purpose. [START_REF] Denieul | Control Allocation for multicontrol surfaces applied to unconventional aircraft configurations[END_REF] So to keep comparison as fair as possible, total area devoted to control surfaces is kept as close as possible to the reference configuration when evaluating alternative layouts.

From this point, two layouts were compared to the initial one in order to test the fast prototyping method presented here. Elevons number is kept constant. The most inboard elevon remained constant due to hereabove mentioned constraints. Total span is still devoted to controls, only relative span -and chord when necessary -are changed. Two designs are presented here:

• an Iso-area design, visible on figure 3(b), where an equality constraint was set on control surfaces areas of 2, 3, 4 and 5. Relative chords of 3, 4 and 5 are set constant at 22%. Only the inboard chord of 2 has to be slightly decreased compared to the initial design in order to provide a feasible solution. The overall area is then a bit decreased.

• an Iso-inertia design, visible on figure 3(c), where total area is set equal to the initial total area in order to provide a feasible design. Then relative chords of 3, 4 and 5 are set as free parameters to fulfill the equality requirement on all control surfaces inertia with respect to their hinge.

All three layouts geometrical characteristics are summarized in table 1. The most significant difference between the derived layouts and the initial one is that the largest elevon, which is initially inboard -namely elevon 2 -is now outboard -namely elevon 5 -particularly for the "iso-inertia" design .

Span, m

Surface, m 2 Inertia, kg.m 2 Number Init. Iso-surf. Iso-inert. Init. Iso-surf. Iso-inert. Init. Iso-surf. Iso-inert.

II.B.3. Aerodynamics of the Derived Layouts Computation

Aerodynamic model of the derived layouts is described in this section. For sake of clarity the initial layout is called "reference" aircraft, and the two derived layouts -namely iso-area and iso-inertia layouts -are called "project" aircraft. The reference flight dynamics model described in section II.B.1 is used for all three configurations from section II.B.2 as planform and airfoils are kept constant for all configurations; therefore only control surfaces aerodynamic efficiencies need to be evaluated for the two derived layouts. For that purpose the Athena Vortex Lattice (AVL) software [START_REF] Drela | AVL-Aerodynamic Analysis, Trim Calculation, Dynamic Stability Analysis, Aircraft Configuration Development[END_REF] was used together with calibration factors coming from the supposedly known aerodynamic coefficients of the initial BWB design. More precisely the lift coefficient C L introduced in Eq.4 comprises a i-th control surface deflection dependency C Li which can be written as:

C Li = k N L δmi (α, δm i )C L δm i δm i ( 10
)
where k N L δmi (α, δm) accounts for loss of control surface efficiency as a function of the angle of attack α and control surface deflection δm i , and C L δm i is the lift gradient of the i-th control surface. k N L δmi (α, δm) is known for reference aircraft, and kept for project aircraft. Reference lift gradient C ref L δm i is supposed known from previous studies, and is compared with the lift gradient computed by AVL for the reference configuration C AV L L δm i (see figure 5(a)). From these data a calibration factor accounting for AVL lack of accuracy is computed:

∆C AV L Li = C ref L δm i C AV L L δm i (11)
Then an AVL computation is run on project geometry in order to compute the lift gradient of project aircraft C AV L proj L δm i . This gradient is finally calibrated using the previously computed calibration factor ∆C AV L Li of Eq.11 (see figure 5

Trim Capabilities Evaluation

• Compensator gains

• Actuators bandwidth

• Maneuverability / Trim segregation

C proj Li = k N L δmi ∆C AV L Li C AV L proj L δm i δm i (12)
A similar process is used for computing control surfaces pitching moment efficiencies of project configurations. More precisely the lift coefficient C m introduced in Eq.6 comprises a i-th control surface deflection dependency C mi which can be written as:

C mi = k N L δmi (α, δm i )C L δm i δm i (X CG -X Fi ) ( 13
)
where X CG denotes the x-wise CG position, and X Fi denotes the i-th control surface aerodynamic center x-wise position. However X Fi is not a direct output of AVL and has to be computed as follows:

X AV L Fi = X CG - C AV L m δm i C AV L L δm i (14)
This leads to computing the aerodynamic center calibration factor for the reference aircraft knowing the actual aerodynamic center X ref Fi for the i-th elevon:

∆X AV L Fi = X ref Fi X AV L Fi ( 15
)
Finally an AVL computation is run on project aircraft in order to compute the pitching moment gradient of project aircraft i-th elevon C AV L proj m δm i , and the resulting aerodynamic center is computed and calibrated as follows:

X AV L proj Fi = ∆X AV L Fi . X CG - C AV L proj m δm i C AV L proj L δm i (16)
This method, summarized on figure 4, combines the advantages of fast data generation through light CFD computation, and far better accuracy than AVL direct output through accurate knowledge on a reference configuration. Here it is applied only to control surfaces aerodynamic coefficients, but the process would be similar for computing any aerodynamic coefficient of a project aircraft with a baseline knowledge on a similar reference aircraft.

III. Integrated Design and Control

In this section the problem of integrated design and control is formulated.

III.A. Plant-Controller Optimization General Problem

The usual way of designing and stabilizing an aircraft is a sequential approach: first the aircraft is sized taking into account disciplines such as aerodynamics, weight, engine sizing, performance and open-loop handling qualities. Then when the overall aircraft design is mostly frozen control laws are designed in order to improve closed-loop handling qualities and passenger comfort, or to alleviate loads. The assumption for this sequential procedure is that the open-loop plant is adequately designed in order to guarantee good control performance. From an engineering point of view, this is generally true for classical tail-aft configurations, for this configuration was historically selected as being flyable without any control law. This is no more the case for BWB configuration. Also from a theoretical point of view, sequential plant-controller design was shown to lead to suboptimal configurations. [START_REF] Fathy | On the Coupling Between the Plant and Controller Optimization Problems[END_REF] The idea is then to simultaneously optimize the plant and controller to achieve a globally optimal design. This problem, known in litterature as Plant-Controller Optimization, Integrated Design and Control or Co-Design, was adressed in different domains such as spacecraft design for controlling flexible appendages, 23 chemistry 24 and underwater vehicles. [START_REF] Silvestre | Plant-Controller Optimization with Applications to Integrated Surface Sizing and Feedback Controller Design for Autonomous Underwater Vehicles[END_REF] In aircraft design several approaches were investigated. Already in the 70s' a study by the US Air Force showed the benefits of RSS for a control-configured vehicle in terms of operating weight empty (OWE) and gross weight (GW). More recently studies by Perez et al. [START_REF] Perez | Multidisciplinary optimization framework for control-configuration integration in aircraft conceptual design[END_REF] the problem of integrated design and control for spacecraft avionics was addressed using structured control laws and plant parameter by Alazard et al. [START_REF] Alazard | Avionics/Control co-design for large flexible space structures[END_REF] Later this formulation was extended to aircraft design [START_REF] Denieul | Integrated Design and Control of a Flying Wing Using Nonsmooth Optimization Techniques[END_REF] and gain-scheduling techniques. [START_REF] Lhachemi | A structured -based optimization approach for integrated plant and self-scheduled flight control system design[END_REF] This structured controller approach is used in the present paper. Section III.B formulates the control problem, and section III.C presents the co-design approach.

III.B. H 2 /H ∞ Formulation for Perturbation Rejection on the Acceleration Sensitivity Function

First it is important to remind that the aim of this study is not to provide directly implementable control laws, rather a quick prototyping of control laws and control surfaces layouts for different configurations. Therefore static state feedback u = -K.[α q θ] T of the angle-of-attack, pitch rate and pitch angle respectively, is used. This corresponds to some kind of augmented pitch damper: on most recent aircraft degraded control law feature only direct order of the pilot to control surfaces, in addition with a pitch rate feedback. However it was found that for such unstable configurations as the BWB pitch-rate feedback is not sufficient to properly stabilize the short-period oscillation. Therefore α and θ feedback were added, in order to add some feedback "stiffness". The u command vector is composed of ten elements: a different feedback order is sent to each of the ten control surfaces. This allows for a linear control allocation strategy of the control law optimizer: elevons are differently used according to their relative pitch efficiencies.

Aircraft Longitudinal Model The control problem formulation is based on a H ∞ weighting on the acceleration sensitivity function. This approach is extensively described in work by Fezans et al. [START_REF] Fezans | H∞ control design for generalized second order systems based on acceleration sensitivity function[END_REF] The rationale for using this approach in our problem is the assumption that from a sizing point of view, perturbation rejection is more stringent than command tracking. In other words, it is harder to reject unknown disturbances than to follow known commands. Therefore we are focusing on perturbation rejection -ie stabilizing the aircraft while facing external disturbance such as turbulence -instead of command tracking such as maneuver response to pilot input. This assumption would deserve to be checked, however this is out of scope of the study.

w s 2 +2ξωs+ω 2 s 2 z inf u z 2 α q θ qa q Weighting Function
We consider an unknown disturbance w that perturbates the aircraft pitch acceleration qa , leading to a resulting pitch acceleration q = qa + w (see figure 6). This is equivalent to considering a disturbing moment on pitch axis. The objective of the control law is then to limit the influence of perturbation w on closed-loop acceleration. In order to obtain the desired behaviour for closed-loop pitch acceleration a weighting function

W 1 = s 2 +2ξωs+ω 2 s 2
is specified on resulting pitch acceleration q. If this frequency template is fulfilled, the pitch motion will behave as a second-order system with specified damping ξ and pulsation ω. The H ∞ constraint can now be written as follows:

W 1 T w→ q ∞ ≤ γ ∞ ( 17
)
∞ being the H ∞ norm, and γ ∞ being a value slightly above 1. In practice we will assume γ ∞ = 1.5. However this problem is under-constrained, and without further specifications commands are allowed to raise to high values. To avoid this, a minimization objective is set on the H 2 norm of the transfer between the disturbance w and the command u. The mixed H 2 /H ∞ problem is then formulated as follows:

min K T w→u 2 subject to: T w→z inf ∞ ≤ γ ∞ ( 18
)
z inf being defined on figure 6. In other words: over the set of controllers satisfying closed-loop requirements for perturbation rejection, choose solution with minimal control energy.

III.C. Co-design on Actuators Bandwidth

In section III.B the initial control problem was introduced. In this section the co-design problem is formulated. The formulation presented here was more extensively developed in previous studies. [START_REF] Denieul | Integrated Design and Control of a Flying Wing Using Nonsmooth Optimization Techniques[END_REF] Aicraft Longitudinal Model

α q θ ω i ω i +s Actuators w z inf K opt u act y act = u s 2 +2ξωs+ω 2 s 2
Weighting Function q Actuators dynamics of each elevon is now modelled by a first order transfer of the form:

y act u act (s) = ω i ω i + s , i = 1...10 (19)
where y act and u act are the actuators output and input respectively, as defined on figure 7. ω i represents the i-th elevon bandwidth. This bandwidth is not a fixed parameter, but a variable which is optimized conjointly with control law synthesis. Ω = [ω 1 , . . . ω 10 ] T is the vector of design parameters. The aim of combined optimization is to minimise the actuators bandwidth required to properly stabilize the aircraft. Indeed as shown in section II, high actuators bandwidth lead to high deflection rates, as well as heavy and high power-consuming actuators. Moreover, smaller control surfaces should be allowed to move faster than bigger ones: therefore it was chosen to weight each bandwidth ω i by the associated control surface inertia J i with respect to its hinge in the objective formulation. The objective function associated to aircraft parameters is then: min Ω . . .

J i ω i . . . ∞ = min Ω max i |J i ω i | subject to: 0 ≤ Ω ≤ Ω max (20)
where Ω max is an upper bound for actuators bandwidth, eg. Ω max = 8Hz, and J = [J 1 , ...J 10 ] T is a vector of all control surfaces inertia. Finally the multiobjective optimization is written as a min-max of the two previously described objectives:

min K,Ω max i { 1 M axAmp T w→uact (s, K) 2 , |J i ω i |} subject to: T w→z inf ∞ ≤ γ ∞ , 0 ≤ Ω ≤ Ω max (21)
Please note that with the introduction of actuators dynamics the energy minimisation criterion acts no more on the transfer T w→u but on T w→uact . This was done on purpose, so that the control energy objective T w→uact (s, K) 2 mostly affects the gains K and not the bandwidth Ω. Said shortly, T w→uact (s, K) 2 is concerned with minimizing values of K, and min max |J i ω i | is concerned with minimizing values of Ω. A trade-off arises between K and Ω: same control performance may be achieved for higher gain values and smaller bandwidth and reciprocally.

In order to be compared and conjointly optimized, the two objective functions are normalized through the parameter M axAmp. This parameter has to be carefully chosen as it drives the trade-off between control and design requirements. Usually, in multi-objective control, design weights are seen only as degrees of freedom for the designer; no attention is paid on physical meaning of these weights. The situation here is different, for one of the objectives comprises physical variables. Therefore the weight M axAmp is calculated as follows:

M axAmp = k × T w→uact 0 2 JΩ 0 ∞ (22)
where:

• T w→uact 0 2 represents minimal energy control, computed with a linear quadratic (LQ) synthesis and no actuators dynamics -ie infinite bandwidth. As this represents an ideal situation, optimal value of the multiobjective problem will necessarily be above this value.

• JΩ 0 ∞ is the value of design objective for initial bandwidth, eg. Ω 0 = 3Hz. As we seek to mimimize this criterion, optimal value will be below this value.

• k represents the factor which we are willing to loose on control energy optimality, in order to gain on actuators bandwidth. Therefore we want k ≥ 1. Indeed in eq.( 21) objectives are normalized to 1, so when optimality is reached following relation is obtained:

T w→uact opt 2 JΩ opt ∞ = M axAmp (23)
So the k factor can be interpreted as:

k = T w→uact opt 2 T w→uact 0 2 . JΩ 0 ∞ JΩ opt ∞ (24)
For this study we found that k = 5 is providing good results for all control surfaces layouts.

From an implementation point of view, we chose to work with the systune and slTunable 37 routines from MATLAB TM for:

• It allows mixed H 2 /H ∞ synthesis and multiobjective optimization.

• It allows structured parameters.

• Bounds on the variables are easily applicable.

• Design and simulation model can be derived directly from a single block diagram file.

• Constraints specifications are not limited to specifications on frequencies, but may also handle pole placement constraints, which may be more suitable to handling qualities purpose. This is out of scope of the study but may be addressed in future work.

Finally this optimization was run on a single flight point, corresponding to the most sizing one from codesign procedure point of view. To find this sizing point, the process is run for the initial configuration on different flight points corresponding to different altitudes and Mach numbers. The sizing point is selected as being the flight point requiring maximal actuators bandwidth in order to stabilize the configuration with adequate performance. Interestingly, it was found that this sizing point, corresponding to a high-altitude low-speed flight point, does not correspond to the maximal absolute instability of the short period oscillation mode. Rather it combines a strong longitudinal instability with a lack of control surfaces efficiency due to low dynamic pressure.

IV. Results

In this section, results of the optimization problem defined in section II and III are presented, and trim capabilities of the different configurations are evaluated. For each layout an optimization process is performed, and performance of the three layouts are discussed. On table 2 numerical values of some meaningful norms introduced in section III.C are presented before and after optimisation process, for the three evaluated layouts. Results are consistent with what was stated in section III.C:

IV.A. Actuators Bandwidth Preliminary Sizing for Different Layouts

• T w→uact 2 increases between initial (LQ without actuators dynamics) and optimal design. This means that gains have to be increased as we move away from an ideal situation in order to stabilize the aircraft with required performance.

• JΩ ∞ is decreased from initial to optimal design.

• T w→z inf opt ∞ is below the γ ∞ value. This indicates that required performance specification is fulfilled for all designs. This may also be seen on figure 9, where frequency responses of the three closed-loop layouts are plotted, as well as the W 1 template for performance specification. Frequency responses are located below the weighting function W 1 , which indicates that performance is satisfied for all control laws.

Looking at optimized actuators bandwidth on figure 8, some significant differences among different layouts are observed.

• Initial layout features three fast outboard control surfaces, corresponding to the three smallest elevons.

Inboard elevons 1 and 2 have a slow dynamics.

• Iso-surface layout features one fast outboard elevon pair, and to a lower extend also elevon 4. All three inboard elevons feature slow dynamics. • Iso-inertia layout is the most balanced one, with only relatively low dynamics actuators. Small differences on numerical values are observed due to difference of pitch efficiencies of different elevonsoutboard control surfaces having more x-wise lever arm-but this is second order.

Most important is that control synthesis provides a control architecture with guaranteed performance, for reasonable gains increase -60% increase from an H 2 norm perspective-from ideal to optimized design. Resulting bandwidth are challenging but not unrealistic. Even though this study is too preliminary to draw some definitive conclusions, it seems that the initial layout is too unbalanced for being efficient, some elevons being very small and other very large.

IV.B. Evaluation of Trim Capabilities under Trim / Maneuverability Segregation Assumptions

The co-design process presented in section IV.A gives some insight about a frequential allocation among control surfaces: it stipulates that in order to provide sufficient perturbation rejection performance, some control surfaces need to move very fast while some other are allowed to have slower dynamics. The former ones may be called "maneuverability" control surfaces, while the latter ones may be called "trim" control surfaces. While the co-design process states that maneuverability control surfaces are sufficient to properly stabilize the aircraft, it is not set yet whether remaining trim control surfaces are sufficient for equilibrating the aircraft. The main interest of such a trim / maneuverability segregation would be to use different actuators technologies for these two functions. Indeed it was shown by Garmendia et al. [START_REF] Daniel | Uncertainty Quantification for the Actuation Power Requirements of a Hybrid Wing Body Configuration with Electrically Actuated Flight Control Surfaces[END_REF] allocating only trim to some elevons allowed for choosing irreversible actuators: hence no power penalty is paid in permanent flight once the control surface is deflected to its equilibrating position.

For the three layouts, trim control surfaces are selected as those with actuator bandwidth below 1Hz. For each configuration, elevons evaluated for trim are summarized in table 3.

• δm = δm min minimal control surfaces deflection allocated for trim.

Longitudinal free variables to solve flight mechanics equations are:

• F the thrust.

• α the angle of attack.

• X CG the x-wise CG location. On all evaluated layouts roll control is allocated to elevons 4 and 5. For the initial configuration and according to co-design results of figure 8(a), capacity of equilibrating the aircraft with only control surfaces 1 and 2 is evaluated. Results are presented on figure 10(a). It can be seen that trim on glide criterion is not challenging, for forward and backward limits respectively far exceed allowable CG positions. Same conclusion is valid for trim turn and trim stall forward criteria, which lie forward the W&CG diagram. However trim stall backward is not satisfied for some backward CG positions. Therefore a second computation is run using all elevons in for trimming the aircraft for this specific flight case. Indeed trim stall is a certification maneuver which needs to be performed with all trim capabilities. Results are presented on figure 10(b). Appart from all other criteria being relaxed, it is noteworthy that trim stall backward moved backward, as expected. However, for some specific combinations of mass and CG this criterion still lies inside the W&CG envelope. After investigations this means that α stall is not reached for those backward configurations; nevertheless required C Lmax is achieved thanks to C L contribution of positive deflections of the elevons for this flight case. This criterion is therefore considered as fulfilled.

IV.B.2. Initial Control Surfaces Layout

IV.B.3. Iso-Surface Control Surfaces Layout

Following what was presented in table 3, trim capacity of configuration iso-area is evaluated using elevons 1, 2 and 3 as trim devices. Results are presented on figure 11. One can see that all criteria are satisfied, except the trim stall backward criterion. As stated in section IV.B.2, this criterion allows for using all available control surfaces for trim. with such an allocation, the aircraft would be trimmable for backward CG at stall speed.

IV.B.4. Iso-Inertia Control Surfaces Layout

This last configuration features all control surfaces used both for maneuverability and trim, as presented on figure 8(c). Results of trim evaluation using all control surfaces are presented on figure 12. Once again all criteria but the trim stall backward are fulfilled. As stated previously the apparent lack of control efficiency for this maneuver comes from the formulation of the criterion, specifying a target α stall . However it was

V. Conclusion

A new method for fast prototyping of control surfaces actuators for unconventional unstable configurations is proposed. This method relies upon latest developments of nonsmooth optimization for structured controllers design. Stabilizing control law gains and actuators bandwidth are optimized in a single process. The result is guaranteed closed-loop performance as well as trim / maneuverability control surfaces segregation. In a last part of the study remaining capability of trim elevons is evaluated. The flexibility of the method is demonstrated through evaluation of three different control surfaces layouts for a blended wing-body. As a preliminary recommendation the authors suggest that control surfaces layouts with highly

  (a) Grumman X-29 experimental aircraft. (b) A380. (c) Airbus Blended Wing-Body project.

Figure 1 .

 1 Figure 1. Illustration on unstable aircraft configuration for different design purposes.

Figure 2 .

 2 Figure 2. Weight & CG diagram of the studied BWB, and aerodynamic center in low speed (blue) and high speed (red) in percentage of Mean Aerodynamic Chord (MAC).

  (a) Initial control surfaces layout. (b) Iso-area layout.

Figure 3 .

 3 Figure 3. Control surfaces layouts evaluated.

Figure 4 .

 4 Figure 4. General process of the study including computation of aerodynamic model for project aircraft.

  (a) Comparison of pitch gradient coefficients for reference aircraft: C ref L δm i (dotted) vs C AV L L δm i (plain). (b) Comparison of pitch gradient coefficients for project aircraft before and after calibration respectively: C AV L proj L δm i (plain) vs C proj L δm i (dotted) Projet aircraft here is "iso-area" layout.

Figure 5 .

 5 Figure 5. Comparison of pitch gradient coefficients C L δm i from AVL outputs with reference aircraft aerodynamic data (left) and with calibrated data (right). Indices denote elevon number.

Figure 6 .

 6 Figure 6. Standard problem of perturbation rejection of the acceleration sensitivity function

Figure 7 .

 7 Figure 7. Closed-loop problem of integrated design and control of state-feedback and actuators bandwidth. Open-loop aircraft dynamics is represented in green, weighting function in blue and tunable parameters in orange.

  (a) Initial control surfaces layout. (b) Iso-surface layout. (c) Iso-inertia layout.

Figure 8 .

 8 Figure 8. Control surfaces inertia (normalized), pitch efficiencies and actuators bandwidth preliminary sizing through integrated design and control procedure for different control surfaces layouts.

Figure 9 .

 9 Figure 9. Frequency responses of Tw→z inf for initial (blue), iso-surface (green) and iso-inertia (red) layouts. W1 weighting function is also represented.

  (a) Weight and CG diagram for initial layout, with elevons 1 and 2 used for trim. (b) Weight and CG diagram for initial layout, with all elevons used for trim.

Figure 10 .

 10 Figure 10. Weight and CG diagrams for initial layout.

Figure 11 .

 11 Figure 11. Weight and CG diagram for "iso-surface" layout, with elevons 1 2 and 3 used for trim.

Figure 12 .

 12 Figure 12. Weight and CG diagram for "iso-inertia" layout, with all elevons used for trim.

Table 1

 1

	1	2.7	2.7	2.7	11.9	11.9	11.9	760	760	760
	2	8.8	4.5	3.1	38.4	18.9	12	2430	1100	590
	3	7	4.6	1.7	17.3	18.9	7.9	350	1090	590
	4	7.3	7.6	3.1	11.1	18.9	12	60	390	590
	5	8.3	14.7	19.9	6.5	18.9	41.2	23	106	590

. Control surfaces characteristics for three different layouts: initial, iso-surface and iso-inertia. Control surfaces are numbered from most inboard (1) to most outboard

(5)

.

  (b)), giving:

		Reference Aircraft
		Geometry & Aerodynamics
	Project Aircraft	AVL	Project Aircraft
	Calibrated Aerodynamics	Aerodynamic Calculation	Geometry
	Computation of Aerodynamic
	Model for Project Aircraft
	Integrated Design and
	Control Law synthesis
	Specifications on control law:
	• Perturbation rejection
	• Minimal energy control
	• Minimal actuators bandwidth

  Such sophisticated dynamic controllers require reduction techniques in order to be practically implementable. With recent advances in the field of nonsmooth optimization techniques for structured H ∞ synthesis31 , 32

	incorporated a stability and control (S&C) module into a multidisciplinary
	optimization (MDO) process, showing improved aircraft design over a traditional process. Another approach
	of interest consists in integrating plant relevant parameters into a controller optimization. Plant and con-
	troller parameters are then optimized simultaneously in the controller optimization. Most formulations to
	our knowledge involve linear matrix inequalities (LMI) formulations of open and closed-loop handling quali-
	ties requirements. In the work of Niewoehner and Kaminer, a sequential optimization of control surfaces and
	controller is performed using LMI optimization for state-feedback controller 27 .

[START_REF] Niewoehner | Integrated aircraft-controller design using linear matrix inequalities[END_REF]

More recent work of Liao et al.

[START_REF] Liao | An LMI-based optimization approach for integrated plant/output-feedback controller design[END_REF]

involve H ∞ plant-controller optimization also based on LMI for a single aircraft parameter, namely the elevator size. This work is extended to mixed H 2 /H ∞ control with optimization of several aircraft parameters in a related study.

[START_REF] Liao | Mixed H 2 /H∞ Sub-Optimization Approach for Integrated Aircraft/Controller Design[END_REF]

However these approaches rely on classical H ∞ therory, and therefore lead to full-order controllers.

Table 2 .

 2 Different norms for initial and optimal parameters.

of 19 American Institute of Aeronautics and Astronautics

Acknowledgments

This work is part of a CIFRE PhD thesis in cooperation between ISAE-Supaero, ONERA and the Future Projects Office of Airbus Operations SAS.

B Aircraft Athena Vortex Lattice BWB Blended Wing-Body CG

Control surfaces layout Initial Iso-area Iso-inertia Elevons used for trim 1-2 1-2-3 1-2-3-4-5

Table 3. Control surfaces used for trim according to co-design results for the three layouts.

IV.B.1. Trim Criteria Description

CS-25 regulation states that it must be possible to trim the aircraft during every phase of the flight. Three relevant trim criteria are evaluated in this study: trim on glide, trim stall and trim turn.

1. Trim on glide: it must be possible to trim the aircraft on a 3 • glide path approach, for any combination of mass and CG. To evaluate this criterion eq. ( 4) to ( 8) are solved at the equilibrium for different masses with following constraints:

• V = V app the approach speed.

• γ = -3 • the flight path angle.

• δm = δm min , respectively δm = δm max , minimal, resp. maximal control surfaces deflection allocated for trim.

Free variables to solve flight mechanics equations are:

• F the thrust.

• α the angle of attack.

• X CG the x-wise CG location.

Please note that in our formulation the CG is not restricted to its allowable locations, as defined on W &CG diagram of figure 2 for instance. Rather a minimal -upwards -trim control surfaces deflection is imposed, and maximal forward X CG position that such a deflection is able to balance is computed. If this X CG position lies forward the allowable CG envelope, then trim elevons are sufficiently efficient. Similarly a minimal -downwards -trim deflection leads to a backward X CG position, which should lie backwards the allowable CG envelope.

2. Trim stall : The aircraft should remain trimmable and controllable near stall angle of attack, for any combination of mass and CG. Similarly to the trim on glide criterion, equilibrated flight mechanics equations are solved with following constraints:

• α = α stall the stall angle of attack.

• γ = 0 • .

• δm = δm min , respectively δm = δm max , minimal, resp. maximal control surfaces deflection allocated for trim.

Free variables to solve flight mechanics equations are:

• θ the pitch angle.

• V the air speed.

• X CG the x-wise CG location.

3. Trim turn: this criterion is no more a pure longitudinal one. The aim of this criterion is to check the aircraft manoeuvrability in the following case: a coordinate turn at 45 • , landing speed and the trim at minimum setting. Lateral flight dynamics equations are not described in this paper for sake of brevity, therefore only longitudinal constraints and free variables are listed below. However the actual equilibrium includes constraints on sideslip β, r and p yaw and roll rates respectively, and bank angle φ = 45 • . Longitudinal constraints are:

different elevons sizes should be avoided, unless a segregation between maneuverability and trim is performed, as it is investigated in this paper. Therefore from this study the so-called "iso-surface" layout seems promising compared to the initial one. Future work may include refined actuators modelling in order to assess relative weight and energy penalty for different configurations. Lateral control laws may also be included into the design process, for multi-control elevons may be sized by longitudinal-lateral maneuvers. A last study of interest would be to compare control surfaces layout with different number of elevons.

Appendix

Developing the coefficients in state-space matrices of section II.B.1 gives:


Maria Vakalopoulou

Christos Platias

Maria Papadomanolaki

Nikos Paragios

Konstantinos Karantzalos

SIMULTANEOUS REGISTRATION, SEGMENTATION AND CHANGE DETECTION FROM MULTISENSOR, MULTITEMPORAL SATELLITE IMAGE PAIRS

Keywords: Simultaneous registration, segmentation and change detection from multisensor, multitemporal satellite image pairs Multimodal, Deimos, Iris, Satellite, Video, Classification, Deep learning, CNN

HAL is

INTRODUCTION

Currently, there are more than forty satellites in operation offering remote sensing data solutions from multiple sensors with various spatial, spectral and temporal resolutions. By 2024, this number is expected to climb to 160, including earth observation satellites from both private enterprises and governments, excluding the numerous nanosatellites and cubesats. It is, therefore, critical to build the capacities for operational exploitation of these multisensor, multiresolution, multitemporal datasets towards adequate return on initial investments as well as efficient spatio-temporal environmental monitoring.

Along with data handling, serving and pre-processing issues, multimodal data fusion techniques hold a primary role in information extraction and exploitation. Despite the important research effort during the last decades ( [START_REF] Gmez-Chova | Multimodal Classification of Remote Sensing Images: A Review and Future Directions[END_REF] and the references therein) there are still important challenges to be addressed. Towards this direction, the Image Analysis and Data Fusion Committee of the IEEE Geoscience and Remote * Corresponding author's email: mariavak@central.ntua.gr Sensing Society organises annual data fusion contests [START_REF] Dalla Mura | Challenges and Opportunities of Multimodality and Data Fusion in Remote Sensing[END_REF] in order to highlight efficient approaches [START_REF] Longbotham | Multi-Modal Change Detection, Application to the Detection of Flooded Areas: Outcome of the 2009-2010 Data Fusion Contest[END_REF][START_REF] Berger | Multi-Modal and Multi-Temporal Data Fusion: Outcome of the 2012 GRSS Data Fusion Contest[END_REF][START_REF] Liao | Processing of Multiresolution Thermal Hyperspectral and Digital Color Data: Outcome of the 2014 IEEE GRSS Data Fusion Contest[END_REF] and way forward. In this paper, a novel methodology has been developed able to ingest information from multisensor datasets of different spatial, spectral and temporal resolutions as those of the 2016 contest.

In particular, we extended the formulation of [START_REF] Vakalopoulou | Simultaneous Registration and Change Detection in Multitemporal, Very High Resolution Remote Sensing Data[END_REF][START_REF] Vakalopoulou | Graph-based Registration, Change Detection and Classification in Very High Resolution Multitemporal Remote Sensing Data[END_REF] by adding another graph which is related to the segmentation problem. The proposed formulation jointly considers datadriven costs regarding the segmentation likelihoods (various classification approaches can be employed e.g., [START_REF] Lagrange | Benchmarking classification of earth-observation data: From learning explicit features to convolutional networks[END_REF]), registration metrics (e.g., similarity metrics) and change detection scores. These energies are efficiently coupled with local geometric constraints in the context of a higher order graph. Reduction methods are used to map this graph into a pairwise one which is then optimized using efficient linear programming. Promising experimental results indicating less than 2 pixels in mean displacement errors regarding the registration, above 77% in most cases regarding the segmentation completeness and correctness rates and around or above 70% in change detection demonstrate the high potentials of the developed method.

METHODOLOGY

Energy Formulation

Following the notations of [START_REF] Vakalopoulou | Simultaneous Registration and Change Detection in Multitemporal, Very High Resolution Remote Sensing Data[END_REF][START_REF] Vakalopoulou | Graph-based Registration, Change Detection and Classification in Very High Resolution Multitemporal Remote Sensing Data[END_REF], here, we add another graph G seg associated with the image segmentation problem. The first graph G reg involves nodes where the labels correspond to deformations vectors from the registration process (mapping between the source and the target images), the second G ch refers to nodes with binary labels expressing changes in the temporal domain, while the last G seg refers to the labels representing the segmentation structures being present in the image. The proposed energy function (1) consists of three terms and couples the three different graphs to one:

E reg,ch,seg (l reg p , l ch p , l seg p ) = E ch + E seg + E reg (1)
The labels for each node of the coupling graph will be l p = [l ch p , l seg p , l reg p ]; l ch p ∈ {0, 1} represent the labels for change detection, l seg p ∈ {0, 1} are the labels for binary segmentation and l reg p are the registration labels (with l reg p ∈ ∆ where ∆ = [d 1 , . . . , d n ] corresponds to all possible displacements). Concluding the label space can be summarized as

L ∈ {0, 1} × {0, 1} × ∆.
Both the registration and change detection terms are following the same formulations as in [START_REF] Vakalopoulou | Simultaneous Registration and Change Detection in Multitemporal, Very High Resolution Remote Sensing Data[END_REF][START_REF] Vakalopoulou | Graph-based Registration, Change Detection and Classification in Very High Resolution Multitemporal Remote Sensing Data[END_REF], while the goal of segmentation is to assign the correct segmentation label to each node of the target image. The segmentation graph, therefore, contains a potential term with the classification score for each class (2) and a pairwise term (3) which penalises different segmentation labels. Let us denote by v s,It (x) the feature vector for the target image at every point x and by Ψ l seg p (•) the classification score for each label l seg p . η is the projection function.

V seg (l seg p ) = Ω η(||x -p||)Ψ l seg p (v s,It (x))dx (2)
V pq,seg (l seg p , l seg q ) = ||l seg p -l seg q || (3)
It should be noted that for simplicity, the segmentation graph is defined only at the target image, and co-segmentation will be achieved through the coupling of the segmentation labels with the change detection ones.

Coupling the Energy Terms

The coupling between the three terms is performed by one potential term (4) which penalises different segmentation labels in the absence of change and reverse, between the source and the target images for all possible displacements.

V reg,ch,seg (l reg p , l ch p , l seg p ) = (1 -l ch p ) Ω η(||x -p||)Ψ l seg p (v s,Is (x + d l reg p ))dx + l ch p Ω η(||x -p||) (Ψ (1-l seg p ) (v s,Is (x + d l reg p )))dx (4)
Let us abuse the notation and consider a node with an index p ∈ G (we recall that the three graphs are identical) corresponding to the same node throughout the three graphs (G reg , G ch , G seg ). Then the global energy (5) can be expressed as follows:

E reg,ch,seg (l reg p , l ch p , l seg p ) = w 1 p∈G V reg,ch (l reg p , l ch p )+ w 2 p∈G V reg,ch,seg (l reg p , l ch p , l seg p ) + w 3 p∈Gseg V seg (l seg p )+ w 4 p∈Greg q∈N (p)
V pq,reg (l reg p , l reg q )+ w 5

p∈G ch q∈N (p)
V pq,ch (l ch p , l ch q )+ w 6 p∈Gseg q∈N (p)

V pq,seg (l seg p , l seg q )

(5) and each modality, a deep neural network classification approach was employed. Such deep learning frameworks have reported high classification accuracy rates for a number of cases [START_REF] Lagrange | Benchmarking classification of earth-observation data: From learning explicit features to convolutional networks[END_REF][START_REF] Mnih | Machine Learning for Aerial Image Labeling[END_REF][START_REF] Vakalopoulou | Building Detection in Very High Resolution Multispectral Data with Deep Learning Features[END_REF]. Training, testing and validation polygons were created under a semi-automated procedure and after an initial registration of all datasets to Google's image mosaic in order to relate every pixel to a Google's map color (i.e., terrain class). Spectral analysis on the derived numerous polygons and probabilities were employed (as in [START_REF] Mnih | Machine Learning for Aerial Image Labeling[END_REF]) in order to define 8 terrain classes i.e.,: Roads, Buildings, Building Shadows, Soil, Sea, Ship/vessels, Vegetation and Vegetation shadows. Clouds on the Deimos March image were manually annotated, while in the Iris video sequence/frame the class Vegetation shadows were merged with Vegetation.

EXPERIMENTAL RESULTS AND EVALUATION

The training for the two Deimos images (D1, D2) was performed on the large (12.760 by 11.000 pixels) overlapping region with 8 classes (e.g., Figure 1), while for the Iris (V) 7 classes were employed. Based on the derived polygons numerous patches of size 21x21 were created by centring each patch to annotated pixels [START_REF] Vakalopoulou | Building Detection in Very High Resolution Multispectral Data with Deep Learning Features[END_REF]. Approximately 200.000 randomly selected patches for each class were used for the satellite Deimos data, while 50.000 for the Iris video frames. The training procedure was accomplished with the use of a Convolutional Neural Network (CNN) which consisted of 10 layers. The network was made of two convolutional layers each followed by tangent and max pooling layers, and ending to two fully connected, a tangent and a linear layers. The model was trained with a learning rate of 1 for 40 epochs, while every 2 epochs the learning rate was reduced at half. Regarding the CNN architecture for the Iris video RGB frames, the overall implementation was the same, whereas the size of each patch was 3x21x21. Implementation details: Regarding the employed parameters, in a similar way as in [START_REF] Vakalopoulou | Simultaneous Registration and Change Detection in Multitemporal, Very High Resolution Remote Sensing Data[END_REF] we used 3 image and 4 grid levels. The parameter C was set to 1800. The Sum Absolute of Differences plus Gradient Inner Products (SADG) was used as a similarity metric. All the parameters and the weights had been set using grid search. For the optimisation procedure, the FastPD was employed 1 . Registration results: For the validation of the registration results (Figure 2), several Ground Control Points (GCPs) were manually collected in all resulting image pairs. In Table 1 the mean displacement errors for both axis and the distance in pixels are presented. It should be mentioned, that the reg-Table 1: Quantitative evaluation for the registration results. istration process didn't manage to address the largest relief displacements of the tallest buildings/skyscrapers of this part of Vancouver and these errors hindered both the segmentation and change detection results. All other building rooftops, roads, terrain classes were registered with a sub-pixel accuracy. A quantitative comparison is also provided in Table 1 (first four columns) based on [START_REF] Karantzalos | Efficient and Automated Multi-Modal Satellite Data Registration through MRFs and Linear Programming[END_REF].

Segmentation results: For the segmentation and change detection results, the quantitative evaluation of the framework was performed using the completeness, correctness and overall quality criteria at object level. The resulting true positives, false negatives and false positives were calculated on the validation dataset after the application of the developed framework. The framework was tested for three different classes i.e., Buildings, Ship/vessels and Vegetation. After the optimisation and based on the polygon of the class Sea derived automatically from Google map all segmented as Buildings objects in the Sea were ignored and respectively all Ship/vessels in the land, as well.

As one can observes in Figure 3, although the classification scores (e.g., Figure 1) constrain significantly the solution, the developed framework integrates both scores and similarities ameliorating the segmentation results in several image regions. The quantitative evaluation (Table 2) indicated that the detection completeness rates were above 78% (apart from the class Buildings in the D1-D2 pair) while the detection correctness rates were above 72% in all cases. The higher rates were for the class Vegetation indicating that the NIR Deimos- Change Detection results: Qualitatively, the same errors were observed on the change detection results (Figure 4) for both image pairs. Quantitative evaluation results (Table 3) indicated lower completeness and correctness rates than the segmentation task, as expected. This was mainly due to a number of false positives in the dense urban regions were the relief displacements were significant due to the tallest buildings and skyscrapers.

CONCLUSION

A novel generic graph-based framework was developed addressing simultaneous registration, segmentation and change detection in multisensor data of different spectral, spatial and temporal resolutions. The quite promising experimental results indicated: (i) for the registration task a mean displacement error of less than 2 pixels, (ii) for the segmentation task in most cases completeness and correctness rates above 77% and (iii) for the change detection around or above 70%. The main errors derived from the important relief displacements around the tallest buildings and from the quite similar reflectance spectra of the different man-made objects including ship/vessels.

  The evaluation of the developed framework was performed on the '2016 IEEE GRSS Data Fusion Contest' dataset which includes one Deimos-2 multispectral image acquired in March'15 [D1], one Deimos-2 multispectral image acquired in May'15 [D2] and one video sequence from the Iris RGB sensor acquired in July'15 [V]. The two Deimos images were radiometrically corrected and then pansharpened based on the standard High Pass Filter method, resulting into an overlapping image pair of approximately 12760 by 11000 pixels. The overlapping image pairs with frames from the Iris video sequence were approximately 4720 by 2680 pixels. In order to employ an additional image/map which could serve as a reference/target map while contributing on the automation of the subsequent training procedure, an image mosaic [G] as well as the corresponding map were downloaded from Google Map APIs. All raw unregistered data as well as several experimental results can be viewed here: http://users.ntua.gr/karank/Demos/DemoContest16.html Classification scores: In order to estimate the required in (2) and (4) classification scores Ψ l seg p (•) for each label l seg p

Fig. 1 :

 1 Fig. 1: Classification map with the dominating scores for the Deimos-2 May'15 (D2) image.

Fig. 2 :

 2 Fig. 2: Chessboard visualization with the resulted registered multisensor data.

Fig. 3 :

 3 Fig. 3: Segmentation results after the application of the developed framework on the D1-D2 pair. Results on D2.

Fig. 4 :

 4 Fig. 4: Change Detection from multitemporal, multi-sensor between: (a) a Deimos March'15 and a Deimos May'15 (D1-D2, left), (b) an Iris video sequence (first frame) and a Deimos May'15 (V-D2, right).

Table 3 :

 3 Quantitative evaluation for the change detection.2 band significantly contributed in class separation. Most segmentation errors were due to false alarms near the port, pier and ship wake on the sea, while Buildings and Roads were in certain cases confused.

Table 2 :

 2 Quantitative evaluation for the segmentation results

D2-D1: w 3 = 340, w 4 = 30, w 5 = 4.5, w 6 = 7.5, V-D2: w 3 = 340, w 4 = 40, w 5 = 10, w 6 = 14

ACKNOWLEDGMENTS

The authors would like to thank Deimos Imaging for acquiring and providing the data used in this study, and the IEEE GRSS Image Analysis and Data Fusion Technical Committee.

6. REFERENCES