D Maalouf

I Tamanaja

E Campos

A Chemori
email: ahmed.chemori@lirmm.fr

V Creuze

J Torres

R Lozano

From PD to Nonlinear Adaptive Depth-Control of a Tethered Autonomous Underwater Vehicle

Keywords: Underwater robotics, Depth control, Parameters uncertainty, PD controller, Adaptive control

This paper deals with an experimental comparison between a proportional derivative (PD) controller and an adaptive nonlinear state feedback one, both applied on a tethered autonomous underwater vehicle. The aim is to show the behavior of the closed loop system in the nominal case for each of these two controllers, and then to test their robustness towards some parameters changes. The PD on one hand has a good performance for systems with an unknown model. The adaptive control law on the other hand is known to adjust the unknown parameters of the plant in order to converge to the desired trajectory. This study shows experimental results performed using each of the above mentioned control laws.

INTRODUCTION

Underwater Vehicles have gained an increase interest in the last decades given the multitude of operations they can perform in various fields. Different control techniques have then been proposed and applied in order to deal with the various challenges arising from the nonlinearities and time varying behavior of the vehicle's dynamics. An overview of the main control techniques for underwater vehicles can be found in [START_REF] Fossen | Marine Control Systems:Guidance, Navigation and Control of Ships, Rigs and Underwater Vehicles[END_REF], [START_REF] Yildiz | State-of-the-Art System Solutions for Unmanned Underwater Vehicles[END_REF] and [START_REF] Yuh | Design and Control of Autonomous Underwater Robots: A Survey[END_REF]. Control design methods can be either linear or nonlinear. In the linear techniques PID controllers are still very common but they do not guarantee system's position stability given the highly nonlinear behavior of the underwater vehicle [START_REF] Yildiz | State-of-the-Art System Solutions for Unmanned Underwater Vehicles[END_REF]). Some more advanced linear techniques such as adding an acceleration feedback to the PID in order to make the vehicle less sensitive to external disturbances or using least square regulators to track time varying reference trajectories can be found in Fossen [2002]. These techniques can also be suitable for uncertain linear systems disturbed by additive white noise and the presence of some immeasurable states. However, limitations when using such techniques are still present and that is why more research has been oriented recently towards nonlinear control schemes. In this latter category, various methods were successfully applied to underwater vehicles such as sliding mode control: [START_REF] Fossen | Sliding control of MIMO nonlinear systems[END_REF] and [START_REF] Healey | Toward an improved understanding of thruster dynamics for underwater vehicles[END_REF] and nonlinear adaptive control: [START_REF] Antonelli | An Adaptive Law for Guidance and Control of Remotely Operated Vehicles[END_REF], [START_REF] Antonelli | A Novel Adaptive Control Law for Autonomous Underwater Vehicle[END_REF] and [START_REF] Fossen | Adaptive control of nonlinear underwater robotic systems[END_REF]. The former is a robust control scheme towards parameters uncertainties and external disturbances. Its main drawback is the chattering phenomenom which directly affects the thrusters leading to a high energetic consumption and possible damages. The latter method, which is model based, ensures the stability of the system through a suitable parameters estimation (i.e. not necessarily the "true" values of the parameters) [START_REF] Slotine | Applied Nonlinear Control[END_REF] ). Exper-imental results showing the efficiency of this method have been performed on the underwater vehicle ODIN in [START_REF] Antonelli | Adaptive control of an autonomous underwater vehicle: experimental results on ODIN[END_REF]. In [START_REF] Maalouf | A novel application of multivariable l1 adaptive control: From design to realtime implementation on an underwater vehicle[END_REF] the first implementation of an L 1 adaptive controller on an underwater vehicle was presented with experimental results. Other nonlinear methods based on intelligent control have also been proposed. They include for instance neural networks and fuzzy logic controllers. Neural networks is a powerful tool since it has a parallel structure, applicability to hardware implementation and multivariable nature. However, training time is too long especially in the presence of parameters changes which makes it hard to be implemented in real-time applications [START_REF] Shaw | Fuzzy control of industrial systems-theory and applications[END_REF]). Fuzzy logic controllers have been tested in simulations such as in [START_REF] Chang | Model-based fuzzy modeling and control for autonomous underwater vehicles in the horizontal plane[END_REF], but they require many trial and error cycles to achieve the desired performance [START_REF] Kim | A novel neuro-fuzzy controller for autonomous underwater vehicles[END_REF] ). To combine the advantages of the previous two mentioned controllers, a study supported by simulations has been provided in [START_REF] Kim | A novel neuro-fuzzy controller for autonomous underwater vehicles[END_REF], where a neuro-fuzzy controller is proposed but no experimental results have been carried out yet to validate the theory. We can find in literature some comparisons among control laws through simulations, [START_REF] Antonelli | On the Use of Adaptive/Integral Actions for Six-Degrees-of-Freedom Control of Autonomous Underwater Vehicles[END_REF] and [START_REF] Campa | Robust control of underwater vehicles: sliding mode control vs. mu synthesis[END_REF] but few studies perform these comparisons empirically. In [START_REF] Smallwood | The effect of model accuracy and thruster saturation on tracking performance of model based controllers for underwater robotic vehicles: experimental results[END_REF], adaptive and fixed-model based controllers have been compared with a PD controller in terms of tracking performance in the presence of thruster saturation and modeling errors. The presented experimental results show that model based controllers exhibit a bad performance in presence of modeling errors compared to the classical Proportional Derivative (PD) controller, while adaptive ones compensate properly these uncertainties. Thrusters saturation deteriorates the performance of all the controllers. In this paper we will experimentally compare a PD controller with an adaptive nonlinear state feedback one. We present experimental results for one degree of freedom, namely the depth.

lirmm-00778061, version 1 -18 Jan 2013

Our contribution lies in the presented experimental comparative study that shows the differences in the robot's behavior regarding each applied method in terms of settling time, and trajectory following tested in various conditions. On one hand, the PD controller was chosen because it is considered the basic one used among the non adaptive schemes and on the other hand, the nonlinear adaptive state feedback controller is the standard adaptive one. Comparing these two methods would give us an idea about the benefits of using an adaptive scheme. Two tests of robustness were carried out. The first one deals with the robustness towards parameters changes. It was technically performed by changing the buoyancy of the vehicle by adding floating balls on the top of the submarine as shown in Fig. 6. The second one deals with the robustness towards punctual external disturbances. The experiments have been conducted using the tethered autonomous vehicle Triton-PR entirely designed and built by the LAFMIA (Franco-Mexican Laboratory specialized in Informatics and Automatic Control). This study has been performed to highlight the importance of applying more advanced control techniques in underwater robotics applications and it can be easily generalized to more complex situations. The paper is organized as follows: in the second part of the paper we present the features of the vehicle we have used, and we remind the basic principles used in its modeling. In the second part, we present theoretical aspects of both controllers that we aim to compare. In the last part of the paper, we detail the real-time experiments that we have conducted, we compare the obtained results, and finally we conclude. The Triton-PR submarine (cf. Fig. 1) is an underactuated vehicle, whose propulsion system consists of five thrusters (the thrusters were built using DC motors with metal propellers). These actuators allow controlling the vehicle's orientation (yaw, pitch, roll) as well as its cartesian position (x and z). The yaw control is provided thanks to differential speed control of the thrusters 4 and 5 (cf. Fig. 1). Pitch control is obtained similarly using thrusters 1, 2 and 3, whereas the roll control employs thrusters 2 and 3. The translational motion along the z axis is regulated by decreasing or increasing the combined speed of thrusters 1, 2 and 3. Similarly, translational motion along the x axis is obtained thanks to thrusters 4 and 5. The different forces induced by the five thrusters are illustrated in Fig. 2. The prototype is equipped with various sensors for position and orientation measurements. An analog 6 DOF IMU (Inertial Measurement Unit) is used for roll and pitch, a compass module (HMC6352, with 0.5 degree heading resolution, and an I2C Interface) is used for the attitude, and a silicon pressure sensor (MPX5999D on-chip signal conditioned, 0 to 150 psi) for depth measurement. To process the sensors' data and insure communication between the control PC, the sensors and the power stage, a RabbitCore RCM3400 board is used. This latter includes a Rabbit 3000 microcontroller operating at 29.4 MHz.

TRITON-PR PROTOTYPE: DESCRIPTION AND DYNAMICS

Description of the experimental platform

Once data are processed by the PC, the computed values are sent to the rabbit to activate the motors' power (controlled by Pulse Width Modulation thanks to MD03 -24V 20A H Bridge Motor Drives with I2C Interface). Real time communication is provided through a RS232 link. Fig. 3 shows a schematic view summarizing the various components of the vehicle's hardware and their interactions.

Dynamic modeling of the system

Throughout the paper, the variables in bold represent matrices and the others represent scalars.

By considering the inertial generalized forces, the hydrodynamic effects, the gravity, and buoyancy contributions as well as the effects of the actuators (thrusters), the dynamic model of an underwater vehicle in matrix form, using the SNAME notation and the representation proposed in [START_REF] Fossen | Marine Control Systems:Guidance, Navigation and Control of Ships, Rigs and Underwater Vehicles[END_REF], is written as:

η η η = J(η) J(η) J(η)ν ν ν M ν M ν M ν +C(ν)ν C(ν)ν C(ν)ν +D(ν)ν D(ν)ν D(ν)ν +g(η) g(η) g(η) = τ τ τ (1)
where ν ν ν = [u, v, w, p, q, r] T ,η η η = [x, y, z, ϕ, ϑ , ψ] T are vectors of velocities (in the body-fixed frame) and position/Euler angles (in the earth-fixed frame) respectively. J(η) J(η) J(η) ∈ R 6×6 is the transformation transformation matrix mapping from the bodyfixed frame to the earth-fixed one. The model matrices M M M, C C C and D D D denote inertia (including added mass), Coriolis-centripetal (including added mass), and damping respectively, while g g g is a vector of gravitational/buoyancy forces. τ τ τ is the vector of forces and torques. The experiments have been performed in a small pool and hence, external disturbances such as currents were not taken into account in the model (1). In the case of our study, the vehicle used will be moving at low velocities, which makes the Coriolis terms negligible. Therefore the dynamics (1) can be rewritten as:

M ν M ν M ν +n(ν, η) n(ν, η) n(ν, η) = τ τ τ (2) with n(ν, η) n(ν, η) n(ν, η) = D(ν)ν D(ν)ν D(ν)ν +g(η) g(η) g(η)
Our vehicle is equipped with 5 thrusters controlling 5 degrees of freedom as described above. The control input dimension will therefore be limited to 5 and is given by:

τ = T Ku τ = T Ku τ = T Ku (3)
where u u u ∈ R 5×1 is the vector of control inputs, K K K ∈ R 5×5 is the force coefficient diagonal matrix, and T T T ∈ R 6×5 is the actuators configuration matrix. An identification test has been performed to determine the relationship between the input voltage to the thrusters and the generated force.

PROPOSED CONTROL SCHEMES

In this section the proposed two controllers are presented: the proportional derivative and the adaptive nonlinear state feedback. They are implemented on the prototype described in section 2 in order to perform trajectory following along the z axis. The full model of the system has been previously presented but in this study the controller will be applied to the dynamics of the vehicle in depth.

Proportional Derivative (PD) controller

This control strategy is based on separated saturation functions. The objective of this control law, based on the use of saturation functions, is to limit aggressive control signals sent to the thrusters that may damage them. Therefore, the control input in depth is given by:

τ z = -σ b z 2 (k z 2 ż) -σ b z 1 (k z 1 (z -z d )) (4)
where k z 1 , k z 2 are positive constants gains, representing respectively the proportional and the derivative gains. σ b z 1 and σ b z 2 are saturation functions. It has been proven in [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF] that this control input can ensure a global stabilization of the system and that there exists a time large enough such that (zz des ) → 0, ż → 0, as t → ∞. The gains of the PD controller were empirically tuned to minimize the ISTSE criterion (Integral of Time Multiplied by squared Error) while compromising between the stability and a fast convergence.

Adaptive Nonlinear State Feedback Controller

The adaptive state feedback controller is a state feedback controller with an adaptation part. It provides an online estimation of the unknown model parameters in order to ensure a proper trajectory following [START_REF] Fossen | Marine Control Systems:Guidance, Navigation and Control of Ships, Rigs and Underwater Vehicles[END_REF]. The control law is given by:

τ τ τ = Ma b Ma b Ma b + n(ν, η) n(ν, η) n(ν, η) (5)
where the hat symbol denotes the parameter estimates, a b a b a b is the body frame commanded acceleration and n(ν, η)

n(ν, η) n(ν, η) is the estimate of n(ν, η) n(ν, η) n(ν, η) in (2).
Since the dynamic model is linear in its parameters, the adaptive control law is then written as:

τ τ τ = Φ(a b , ν, η) θ Φ(a b , ν, η) θ Φ(a b , ν, η) θ (6)
where Φ Φ Φ is the regressor matrix and θ θ θ is the vector of the estimated parameters. To guarantee that error converges to zero, the commanded acceleration in the inertial frame a n a n a n is chosen as the following proportional derivative (PD) control:

a n a n a n = ηd ηd ηd -K d η K d η K d η -K p η K p η K p η (7) with η η η = η η η -η d η d η d , η η η the first derivative of η η η, η d η d η d
is the desired trajectory and ηd ηd ηd the desired acceleration Combining equations ( 1) and ( 7), the expression of a b a b a b , the acceleration in the body-fixed frame becomes:

a b a b a b = J -1 J -1 J -1 (a n a n a n -Jν Jν Jν) (8)
The parameter update law is given by: θ

θ θ = -ΓΦ T (a b , ν, η)J -1 y ΓΦ T (a b , ν, η)J -1 y ΓΦ T (a b , ν, η)J -1 y (9)
where Γ Γ Γ is the adaptation gain diagonal matrix, J J J the transformation matrix and y y y the combined error defined as: y y y = c 0 η η η + c 1 η η η (10) with c 0 and c 1 being constant gains chosen to ensure the stability of y y y. In [START_REF] Fossen | Marine Control Systems:Guidance, Navigation and Control of Ships, Rigs and Underwater Vehicles[END_REF] the complete proof of stability stating that the convergence of the position error to zero is guaranteed by applying Barbalat's lemma. Given the available sensors and actuators, we have chosen the case of trajectory tracking for motions in depth. Given equation (3), we get the vector of control input to be:

u u u = T -1 K -1 T -1 K -1 T -1 K -1 τ τ τ (11)
The (only) parameter to be estimated is (W -B), which represents the difference between the weight W and the buoyancy B. The control input of equation ( 11) will be limited to the degree of freedom around the depth and will be of dimension 1. It is important to emphasize that a good parameter convergence is more guaranteed to occur when the followed trajectory is rich enough to excite the parameters under study [START_REF] Slotine | Applied Nonlinear Control[END_REF]). However, parameters will converge to a set of values that allow trajectory following, but these values are not necessarily the "true" ones. Moreover, the parameter vector will be bounded but not necessarily convergent. The gains of this controller have been tuned similarly to the case of the PD controller in order to minimize the ISTSE criterion.

REAL-TIME EXPERIMENTAL RESULTS

In this section the obtained experimental results will be presented and discussed. They result from the application of the proposed controllers detailed in section 3, to the underwater vehicle testbed described in section 2. The different parameters of the proposed controllers are summarized in table 1. Three experimental scenarios are performed, namely: Scenario 1: Control of Triton-PR in nominal conditions, Scenario 2: Robustness towards modeling uncertainties, Scenario 3: External disturbances rejection.

Table 1. Parameters used in the experiments

Parameter Value

k z 1 2 k z 2 0.2 Γ 1 K p 0.8 K d 0.2 c o 0.2 c 1 0.1

Control of Triton-PR in nominal conditions

For this scenario, the vehicle is considered without any external disturbances. The obtained experimental results are depicted in Fig. 4 and 5 for both controllers, showing their tracking performance. as well as its corresponding reference trajectory (dashed line)

Figure 5 shows the rapid convergence of the system under the adaptive controller to the trajectory after a suitable gain adjustment while a clear overshoot and more relevant oscillations can be seen in the case of the PD controller. We can also see that the final immersion value obtained with the PD controller is affected by a permanent static error (3 cm).

Robustness towards parameters'change

The objective of this scenario is to test the robustness of both controllers towards parameters changes of the system. For that purpose one proposes to change the buoyancy by means of adding two small spherical balls of 5.1 cm diameter, fixed on the top of the vehicle as illustrated in Fig. 6. The added floatability to the system can then be estimated to be around 1.4 N.

Fig. 6. Two floating balls have been added to the vehicle in order to change its buoyancy.

Once the buoyancy has been changed, the two controllers applied in the first scenario were again tested on the Triton-PR. We can clearly note from Fig. 7 that the PD was not able to follow the trajectory given that the system changed due to the added buoyancy. The previous gains adjusted for the nominal case tested before were too small to allow the vehicle to reach the desired reference trajectory. By contrast with the PD, the adaptive control law was able to compensate this parameter change. Adaptation of the parameters lasted 120 seconds and then allowed a better trajectory following. adapted to the buoyancy change by converging to a different value. The delay noted in the convergence of the trajectory ( 120 seconds) in Fig. 7 coincides with the time needed for the parameter to converge to its steady state estimated value. This convergence time could be shortened by increasing the adaptive gain (however, more attention should be paid since this can lead to instability). Once the parameters have adapted, the vehicle's behavior will be the same as in the nominal case. The difference between the estimated values in the nominal case (-0.8N) and the modified case (-3.5N) does not correspond to the value (+1.4N) of the added buoyancy. This is natural as adaptive controllers do not necessarily need the true values of the estimated parameters to ensure good trajectory following [START_REF] Slotine | Applied Nonlinear Control[END_REF].

External disturbances rejection

The objective of this scenario is to test the performance of both controllers in presence of external disturbances. The idea is that after the vehicle reached its steady state, a punctual downward force is applied to it. Fig. 9 and 10 below show the recovery of the system after the application of this external disturbance in each of the above two studied cases (PD and Adaptive state feedback controllers). We notice that the PD controller is more sensitive to external disturbances and needs more time to recover than the adaptive controller. Furthermore, it was noticed that with the PD controller an important steady state error was observed. (solid line) obtained with the adaptive controller as well as the reference trajectory used (dashed line) to test its robustness towards a punctual disturbance at t = 170 seconds.

Remark:

The present paper is accompanied with a short video showing the obtained experimental results available at the following link: http://www.youtube.com/watch?v=oqEsXOrpPPk.

SUMMARY OF COMPARISON BETWEEN THE PROPOSED CONTROLLERS

From Fig. 4, 5, 7, 9, and 10, we can conclude that the PD controller exhibits important limitations in terms of robustness towards external disturbances and parameter change. Under nominal conditions a better trajectory tracking is observed when the adaptive nonlinear state feedback controller is used. In nominal conditions, the steady state error observed with the adaptive controller was around 4% compared to 7% for the PD one. When buoyancy was changed, no important difference was observed on the adaptive controller since the error was around 3% while 44% was noticed with the PD controller. Using a PID controller instead of the PD would cancel the steady state error. The gains are not to be adjusted again after the first tuning for both control laws. The performed experiments show clearly that with small model changes the adaptive controller is able to compensate for them and steer the controlled system to the desired trajectory. Some Matlab T M simulations have been performed previously and it was observed that larger changes (e.g +80% of buoyancy) would have induced saturation on the thrusters. Hence, the scenarios of the performed experiments were chosen to be in a reasonable range of model changes. Such a small increase in buoyancy can be noted for instance in real environments where salinity is likely to change.

CONCLUSION AND FUTURE WORK

An experimental comparison study between a PD controller and an adaptive nonlinear state feedback controller has been performed on an underwater vehicle in one degree of freedom.

Experimental results have shown that the adaptive controller is largely more robust than the PD controller, since it is able to insure a better trajectory following and has a better compensation for external disturbances. Most importantly, it can adapt online unknown or changing parameters and achieve convergence to the desired trajectory. When parameters are modified, the PD controller behaves poorly. The future work will include more parameters (e.g drag parameters) in the estimation vector of the adaptive controller and the implementation of a PID with a systematic way of tuning the gains for richer trajectories involving more degrees of freedom to induce a more complex nonlinear coupled dynamics.

Fig. 1 .

 1 Fig. 1. View of the Triton-PR submarine and its five thrusters.

Fig. 2 .

 2 Fig. 2. Position and orientation of the forces induced by the five thrusters.

Fig. 4 .

 4 Fig. 4. Experiment of a PD controller on Triton-PR: Time history of the measured z-axis (depth) position (solid line) as well as its corresponding reference trajectory (dashed line)

Fig. 5 .

 5 Fig. 5. Experiment of adaptive nonlinear state feedback controller on Triton-PR: Time history of the measured z-axis (depth) position (solid line) as well as its corresponding reference trajectory (dashed line).

Fig. 8 Fig. 7 .

 87 Fig.8shows the estimation of the parameter (W -B) by the adaptive controller for the nominal case and the case with additional buoyancy. One can note that this parameter varied and

  Fig. 8. Time history of the estimated parameter (W -B) resulting from the applied adaptive controller for both scenarios: nominal case and buoyancy change case.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Time history of the measured z-axis (depth) position (solid line) obtained with the PD controller as well as the reference trajectory (dashed line) used to test its robustness towards a punctual disturbance.

Table 2 .

 2 Table 2 below shows the root mean square tracking error at steady state for both scenarios: nominal one as well as the scenario with change in buoyancy. Root Mean Square Error at Steady State (in meters).

	Controller Nominal Case	Buoyancy Change
	PD	0.032	0.20
	Adaptive	0.017	0.011


F Orieux
email: orieux@iap.fr

J.-F Giovannelli

T Rodet
email: rodet@lss.supelec.fr

A Abergel
email: abergel@ias.u-psud.fr

Estimating hyperparameters and instrument parameters in regularized inversion Illustration for Herschel/SPIRE map making

Keywords: data analysis -methods, statistical -methods, numerical -techniques, image processing 1. Unsupervised myopic inversion

We describe regularized methods for image reconstruction and focus on the question of hyperparameter and instrument parameter estimation, i.e. unsupervised and myopic problems. We developed a Bayesian framework that is based on the posterior density for all unknown quantities, given the observations. This density is explored by a Markov chain Monte-Carlo sampling technique based on a Gibbs loop and including a Metropolis-Hastings step. The numerical evaluation relies on the SPIRE instrument of the Herschel observatory. Using simulated and real observations, we show that the hyperparameters and instrument parameters are correctly estimated, which opens up many perspectives for imaging in astrophysics.

The agreement of physical models and observations is a crucial question in astrophysics, however, observation instruments inevitably have defects and limitations (limited pass-band, nonzero response time, attenuation, error and uncertainty, etc.). Their inversion by numerical processing must, as far as possible, be based on an instrument model that includes a description of these defects and limitations. The difficulties of such inverse problems, and notably their often ill-posedness, were well identified several decades ago in various communities: signal and image processing and statistics, and also mathematical physics and astrophysics. It seems pertinent to take advantage of the knowledge amassed by these communities concerning both the analysis of the problems and their solutions.

The ill-posedness comes from a deficit of available information (and not only from a "simple numerical problem"), which becomes all the more marked as resolution requirements increase. The inversion methods must therefore take other information into account to compensate for the deficits in the observations: this is known as regularization. Each reconstruction method is thus specialised for a certain class of objects (point sources, diffuse emission, superposition of the two, etc.) according to the information accounted for. Consequently, in as much as it relies on various sources of information, each method is based on a trade-off, which usually requires the setting of hyperparameters, denoted by ξ in the following. The question of their automatic tuning, namely unsupervised inversion, has been extensively studied and numerous attempts investigate statistical approaches: approximated, pseudo or marginal likelihood, in a Bayesian or non-Bayesian sense, EM, SEM and SAEM algorithms, etc. The reader may consult papers such as (Zhou et al. 1997;de Figueiredo & Leitao 1997;[START_REF] Saquib | [END_REF]Descombes et al. 1999;[START_REF] Molina | [END_REF]Lanterman et al. 2000;Pascazio & Ferraiuolo 2003;[START_REF] Blanc | [END_REF]Chantas et al. 2007;[START_REF] Giovannelli | [END_REF]Babacan et al. 2010;Orieux et al. 2010a) and reference books such as (Winkler 2003, Part.VI), (Li 2001, Chap. 7) or (Idier 2008, Chap. 8). Alternative methods are based on the L-curve (Hansen 1992;[START_REF] Wiegelmann | [END_REF] or on generalised cross-validation (Golub et al. 1979;Fortier et al. 1993;Ocvirk et al. 2006).

The construction of maps of high resolution and accuracy relies on increasingly complex instruments. So, inversion methods require instrument models that faithfully reflect the physical reality to distinguish, in the observations, between what is caused by the instrument and what is due to the actual sky. Then, a second set of parameters comes into play: the instrument parameters, denoted by η in the following, such as lobe width, amplitude of secondary lobes, response time, or gain. Their values are of prime importance and their settings are generally based on dedicated observation and rely on models and/or calibrations that inevitably contain errors. For example, the lobe widths are usually determined from a specific observation in a spectral band of non-zero width; consequently the result depends on the source spectrum. Correction factors can be applied but, naturally, they also contain errors when the source spectrum is poorly known or unknown. In contrast, our aim is to achieve myopic inversion, i.e. to estimate the instrument Article published by EDP Sciences A83, page 1 of 13 parameters without dedicated observation. The question arises in various fields: optical imaging [START_REF] Pankajakshani | [END_REF], interferometry [START_REF] Thiébaut | Proc. SPIE: Astronomical Telescopes and Instrumentation[END_REF], satellite observation [START_REF] Jalobeanu | Proc. IEEE ICASSP[END_REF], magnetic resonance force microscopy (Dobigeon et al. 2009), fluorescence microscopy (Zhang et al. 2007), deconvolution (Orieux et al. 2010b), etc. A similar problem deals with non-parametric intrument response (blind inversion), for which the literature is also very abundant: (Mugnier et al. 2004;[START_REF] Thiébaut | [END_REF]Fusco et al. 1999;[START_REF] Conan | [END_REF] in astronomy and (Lam & Goodman 2000;[START_REF] Likas | [END_REF]Molina et al. 2006;[START_REF] Bishop | Proc. IEEE ICIP[END_REF][START_REF] Xu | [END_REF] in the signal-image literature represent examples. The present paper is devoted to parameter estimation for the instrument model developed in our previous paper (Orieux et al. 2012b), based on an accurate instrument model.

A threefold problem has to be solved: from a unique observation, estimate the hyperparameters the instrument parameters and the map. This is referred to as unsupervised and myopic inversion. From the methodological point of view, the proposed inversion method comes within a Bayesian approach (Idier 2008). In this family, we find the classic Wiener and Kalman methods that calculate the expectation or the maximizer of a posterior density. In an equivalent way, the Phillips-Twomey-Tikhonov methods calculate the minimizer of a least-squares criterion with quadratic penalization. These methods are based on a secondorder analysis (Gaussian models, quadratic criteria) and lead to linear processing. The work proposed here is in a similar methodological vein as far as estimating the map goes; however, the contribution concerns the estimation of the hyperparameters and instrument parameters. We resort to an entirely Bayesian approach (also called full-Bayes) that models the information for each variable (observations, unknown map as well as hyperparameters and instrument parameters) through a probability density. Based on an a posteriori distribution for all the unknown variables, the proposed method jointly estimates the instrument parameters, the hyperparameters, and the map of interest. Regarding experimental data processing, the present paper follows (Orieux et al. 2012b) on inversion for the SPIRE instrument onboard Herschel, which requires the hyperparameters to be fixed by hand and the instrument parameters to be known. The proposed method can automatically tune these parameters and may permit the systematic and automatic processing of large information streams coming from present and future space-based instruments (e.g. Herschel, Planck, JWST, etc.).

The paper is structured as follows. Section 2 introduces the notation and sets out the problem. Section 3 presents the inversion method: it introduces the prior densities and leads to the posterior density. Section 4 describes the computing method based on Markov chain Monte-Carlo (Gibbs) stochastic sampling algorithms. The work is essentially evaluated on simulated observations and on a first set of real observations in the context of the SPIRE instrument onboard Herschel. The results are presented in Sect. 5. Finally, some conclusions and perspectives are provided in Sect. 6.

Notation, instrument, and map models

To produce accurate and reliable maps, the inversion must exploit a description that represents the acquisition process as faithfully as possible. In this sense, the instrument model is based on a map of the sky noted X, which is naturally a function of continuous spatial variables (α, β) ∈ R 2 (and possibly a spectral variable λ ∈ R + );

and accurately describes the formation of a set of N discrete observations grouped together in a vector y ∈ R N .

A general description of the map of the sky as a function of continuous spatial variables can be written starting from a basic function ψ by combination and regular shifting:

X(α, β) = i j x i j ψ α -i δ α , β -j δ β . (1)
The function ψ must be chosen so that this decomposition can describe the maps of interest and is easy to handle. It may be, among other choices, a pixel indicator, a cardinal sine function, or a wavelet (although in the last case the function ψ and the coefficients also depend on a scaling parameter). Whatever the choice, the map of interest is finally represented by its coefficients x i j , the number of which is arbitrarily large and collected in x ∈ R M in what follows. In practice, we choose the Gaussian family as this greatly simplifies the (theoretical and numerical) calculations of the model outputs, including for complex models (Orieux et al. 2012b;Rodet et al. 2008).

The presented work is quite generic in the sense that it is not a priori attached to a specific instrument. It deals with a general linear instrument model that describes, at least to a fair approximation, the physics of the processes in play: optics, electrics, and thermodynamics. It also includes the passage from a continuous physical reality to a finite number of discrete observations. The instrument is then described by

y = A η x + n, (2)
i.e. a general linear model w.r.t. x (a special case of which is the convolutive model). This model shows the instrument parameters η ∈ R K that define the form of the instrument response. The component n = y -A η x represents the measuring and modelling errors additively. For the SPIRE instrument (Griffin et al. 2010) of the Herschel Space Observatory (Pilbratt et al. 2010) launched in May 2009, the paper of Orieux et al. (2012b) gives the details of the instrument model construction. The results of Sect. 5 are based on this instrument.

Probabilistic models and inversion

The proposed inversion is developed in the framework of Bayesian statistics. It relies on the posterior density p(x, ξ, η|y) for the unknown quantities x (image), ξ (hyperparameters), and η (instrument parameters) given y (observations). This density brings together the information about the unknowns in the sense that it attaches more or less confidence to each value of the triplet (x, ξ, η). A summary of this density in the form of a mean and a standard deviation will provide (1) a point estimate (the posterior mean) for the map of interest and the parameters, and (2) an indication of the associated uncertainty (the posterior standard deviation).

Remark 1. In statistical terms [START_REF] Robert | The Bayesian Choice, Statistiques et probabilités appliquées[END_REF], the posterior mean is an optimal estimator. More precisely, of all the possible estimators (whether Bayesian or not, empirical or not, a computation code, etc.), the posterior mean yields the minimum mean square error (MMSE)1 . Regarding first-order statistics, this estimator has, moreover, a zero mean bias.

The posterior density is deduced as the ratio of the joint density for all considered quantities p(x, ξ, η, y) and the marginal density for the observations p(y) by application of Bayes' rule p(x, ξ, η|y) = p(x, ξ, η, y) p(y)

• (3)
Seen as a function of the unknowns (x, ξ, η), this posterior density is proportional to the joint density:

p(x, ξ, η|y) ∝ p x, ξ, η, y .

This joint density is essential as all the other densities (marginal, conditional, prior, posterior, etc.) can be deduced from it. It can be factorised in various forms and, in preparatio for the developments to follow, we write

p(x, ξ, η, y) = p y|x, ξ, η p x, ξ, η = p(y|x, ξ, η)p x|ξ, η p ξ, η = p(y|x, ξ, η)p (x|ξ) p(ξ)p(η), (5)
including the fact that (1) the hyperparameters ξ and the instrument parameters η are a priori independent and (2) the object x and the instrument parameters η are also a priori independent.

The different probability densities will be defined in the following sections according to the information available on each set of variables and according to practical concerns about dealing with the probability densities and numerical computation time.

Modelling of errors and likelihood

The factor p(y|x, ξ, η) in Eq. ( 5) is the density for the observations y given the map x, the instrument parameters η, and the hyperparameters ξ, i.e. the likelihood of the unknowns attached to the observations. Given Eq. ( 2), the construction of this likelihood is based on the model for the error n. The analysis developed in this paper is essentially founded on its mean m n and its covariance matrix Σ n , and the proposed model is Gaussian:

n ∼ N(m n , Σ n ). ( 6
)
The choice of the Gaussian model is also justified via information property: based on the sole information of finite mean and covariance, the Gaussian density is the model that introduces the least information [START_REF] Kass | [END_REF]. This property is also mentioned as a maximum entropy property of the Gaussian density.

m n is a scalar that models a possible non-zero mean for the noise, such as an offset (in the proposed numerical evaluation of Sect. 5, one offset for each bolometer is introduced). Regarding the covariance matrix, to lighten the notation, we set Σ -1 n = γ n Π n : γ n is a scale factor (called precision, homogeneous to an inverse variance) and Π n contains the structure itself. For a stationary model Π -1 n has a Toeplitz structure, for an auto-regressive model Π n is a band matrix, for a independent model Π n is diagonal, for a white and stationary model Π n = I, the identity matrix. In the developments below, the structure of Π n is given while the scale factor γ n and the mean m n are unknown and included in the vector ξ. The results in Sect. 5 are presented for the case Π n = I; hence γ n is the inverse of the noise power.

Remark 2. The proposed developments account for characteristics of the error n that may differ from channel to channel, sensor to sensor, etc. This will be the case in Sect. 5: a mean and power of the noise will be assigned and estimated for each bolometer.

As the error n is Gaussian and additive (Eqs. ( 6) and (2)), the vector of observations y, given x, ξ, η, is also Gaussian y|x, ξ, η ∼ N m y| * , Σ y| * with mean

m y| * = A η x + m n (7)
and with the same covariance as n: Σ y| * = Σ n . So, the likelihood of the unknowns attached to the observations reads

p(y|x, ξ, η) = (2π) -N/2 γ N/2 n det [Π n ] 1/2 × exp - 1 2 γ n y -m y| * t Π n y -m y| * . (8)
It includes the information provided by the observations as the transform of a map x by the instrument, taking its parameters η and the noise parameters γ n and m n into consideration.

Prior density for the map and spatial regularity

The aim of this section is to introduce a prior density p(x|ξ) for the unknown map coefficients x based upon available information about the map X. The present work is mainly devoted to extended emissions. From a spatial standpoint, such maps are relatively regular, i.e. they involve positive correlation. From the spectral standpoint, the power is mainly located at relatively low frequencies. The Gaussian density includes these second-order properties in a simple way. This choice can also be justified based on a maximum entropy principle. Its main interest here is to result in a linear processing method. It is written in the form

p(x|γ x ) = (2π) -M/2 γ M/2 x det Π x 1/2 exp - 1 2 γ x x t Π x x , (9)
where γ x is a precision parameter (homogeneous to an inversevariance that controls the regularity strength) and Π x is a precision matrix (homogeneous to an inverse-covariance matrix that controls the regularity structure). When the precision γ x is low (strong prior variance), the regularity information is weakly taken into account. Conversely, when the precision γ x is high (weak prior variance), the penalization of non-regular maps is high, i.e. the regularity is strongly imposed.

The subsequent developments are devoted to the design of Π x to account for the desired regularity of the map. A simple regularity measure R c [X] of the map X is the energy of some of its derivatives. These derivatives can address the spatial variables (α, β) separately, can rely on cross derivatives and can intervene at various orders. This is the classical Philipps-Twomey-Thikonov penalization idea [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF]. It can also embed directional derivatives or any differential operator [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]). In the simplest and natural case, we choose

R c [X] = ∂X ∂α 2 + ∂X ∂β 2 ,
where u is the standard function norm2 . Given the decomposition (1), it is easy to establish the partial derivatives of X from the derivatives of ψ. In the direction α, by noting ψ ′ α = ∂ψ/∂α, we have

∂X ∂α 2 = i j i ′ j ′ x i j x i ′ j ′ R 2 ψ ′ α α -i ′ δ α , β -j ′ δ β × ψ ′ α α -i δ α , β -j δ β dα dβ, which brings out the autocorrelation Ψ α = ψ ′ α ⋆ ψ ′ α of the deriva- tive of ψ. We then have a quadratic form in x ∂X ∂α 2 = i j i ′ j ′ x i j x i ′ j ′ Ψ α (i ′ -i) δ α , ( j ′ -j) δ β = x t Ψ α x.
As the coefficients Ψ α (i ′i) δ α , ( j ′j) δ β depend only on the difference between indices, the matrix Ψ α has a Toeplitz structure and the computations amounts to a discrete convolution that can be efficiently implemented by the use of fast Fourier transform (FFT). Finally, by performing the same development in the β dimension, a global quadratic norm appears: R c [X] = x t (Ψ α + Ψ β )x and designs the precision matrix

Π x = Ψ α + Ψ β .
For more details and for a spectral interpretation, see Sect. 2.1 and Appendix A of (Orieux et al. 2012b).

Prior distribution for hyperparameters (hyperprior)

The hyperparameters are the unknown parameters of the densities for the error and for map Eqs. ( 8), ( 9) and they are collected in the vector ξ = [m n , γ n , γ x ]. It has been said that γ x and γ n are the precisions (scale parameters) and m n is a mean (position parameter) of Gaussian densities.

The choice of the prior distributions for these hyperparameters is driven by two requirements: (i) little information is available a priori on their values and their relations and (ii) the chosen distributions must lead to efficient algorithms (see Sect. 4.2). Following this line of thought, we choose a prior distribution determined by Jeffreys' principle3 : p(γ) = 1/γ for γ x and γ n and p(m n ) = 1. Moreover, regarding the triplet of hyperparameters ξ = [m n , γ x , γ n ], they are modelled as independent variables, since no information is available about their eventual relations. Finally

p(m n , γ x , γ n ) = 1/γ x γ n , (10)
has two advantages 1. First of all, the posterior conditional densities for γ x and γ n (resp. for m n ), as shown in Sect. 4.2, will be gamma densities (resp. Gaussian density), which will make the implementation easier. 2. This prior distribution is non-informative (which introduces a minimum of information on the value of the hyperparameters) in the sense that it is invariant by certain parameterization changes [START_REF] Robert | The Bayesian Choice, Statistiques et probabilités appliquées[END_REF][START_REF] Kass | [END_REF].

Prior density for the instrument parameters

The instrument parameter η operates in a complex nonlinear way in the description of the observations. In consequence, whatever the prior density, the conditional posterior density for η (see Sect. 4.3) will not have a standard form. The choice is thus purely oriented by the information on the instruments and the question that arises concerns the encoding of the available information in the form of a probability density. If we have no information except a minimum and a maximum value for a given parameter, the choice of a uniform density over the interval is a reasonable one. If we have a nominal value with an associated uncertainty and no other information, the most suitable choice is a Gaussian density. The rest of the development is valid whatever the choice, and we consider the Gaussian case in the following.

y η ρ µ m n Π n γ n x γ x Π x
In addition, having no information available about possible links among the various parameters, we take it that the parameters are, a priori, independent and thus

p(η) = K k=1 (2πρ k ) -1/2 exp - (η k -µ k ) 2 2ρ k • (11)
In practice and for the first results presented in Sect. 5, the means µ k and variances ρ k were taken from the SPIRE observer manual or were fixed ad-hoc at plausible values.

Posterior density: histograms, mean and standard deviation

The posterior density (3) for all the unknowns x, ξ, η, is deduced from the joint density ( 5) for all quantities concerned as follows:

p(x, ξ, η|y) ∝ p(y|x, ξ, η)p(x|ξ)p(ξ)p(η).

In this expression,

the density p(x|γ x ) for the unknown map x ∈ R M is Gaussian (Eq. ( 9)); -the distribution p(ξ) (for ξ = [m n , γ n , γ x ]) is a Jeffreys' distribution (Eq. ( 10)); -the instrument parameters η ∈ R K is modelled by a Gaussian density p(η) (Eq. ( 11)); -the density p(y|x, ξ, η) for the observations y ∈ R N given the rest of the variables (i.e. likelihood) is a Gaussian density (Eq. ( 8)) and it is a function of x through m y| * given by Eq. ( 7).

Finally, from Eqs. ( 8), ( 9), (10), and ( 11), the posterior density can be written

p(x, ξ, η|y) ∝ γ M/2-1 x γ N/2-1 n K k=1 (2πρ k ) -1 × exp - 1 2 (η k -µ k ) 2 ρ k exp - 1 2 γ x x t Π x x × exp - 1 2 γ n (y -m y| * ) t Π n (y -m y| * ) . ( 12
)
This density brings together all information about the unknowns, and the estimators and algorithms presented below are entirely based on it. However, it is too complex to be analyzed directly as a whole and the difficulty stems from (i) the dimension of x (of size M ∼ 10 5 in practice) and (ii) the joint presence of other parameters (hyperparameters ξ and instrument parameters η). Moreover, for the latter in particular, the dependence Table 1. Gibbs algorithm.

Initialize x (0) , ξ (0) , η (0) for q = 1, 2, . . . do (1) sample x (q) under p(x|ξ (q-1) , η (q-1) , y)

(2) sample ξ (q) under p(ξ|x (q) , η (q-1) , y)

(3) sample η (q) under p(η|x (q) , ξ (q) , y) end for is complicated and cannot be identified with a standard form. The proposed approach is to explore the posterior density by means of stochastic sampling [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF][START_REF] Gilks | Markov Chain Monte Carlo in practice[END_REF]. The idea is to produce a set of samples x (q) , ξ (q) , η (q) , for q = 1, 2, . . . , Q, drawn at random under the posterior density. It is then possible, for example, to deduce histograms that approximate marginal densities, together with means and standard deviations. This strategy is by no means new but interest in its practical use has revived in recent years as new forms of algorithms have been developed and computer power has increased.

Concerning the estimates themselves for the map, the hyperparameters and the instrument parameters, we choose the posterior mean (PM), as indicated at the beginning of Sect. 3 (see also Remark 1). We will also look at the dispersion around the mean through the posterior standard deviation (PSD) and the links among components through posterior correlations. Using the set of samples x (q) , ξ (q) , η (q) , for q = 1, 2, . . . , Q, the posterior mean µ P and the posterior covariance matrix Γ P are computed by

µ P ≈ 1 Q Q q=1
x(q) (13)

Γ P ≈ 1 Q Q q=1 x(q) -µ P x(q) -µ P t , (14)
where x denotes the column concatenation x = [x; ξ; η]. Practically, it is not possible to compute the entire covariance Γ P , but it is possible to compute its diagonal elements to characterize the marginal errors for each component (each pixel, hyperparameters, instrument parameters) and to compute a few nondiagonal elements to measure the correlations between components.

Exploration of posterior density and the computation algorithm

We have introduced an instrument model and various probability densities to define the posterior density that brings together the information on the map and the parameters (hyperparameters and instrument parameters). We have also defined the posterior mean (PM) as an estimate and the posterior standard deviation (PSD) as a measure of the uncertainty. We have then introduced the idea of computations via stochastic sampling. The developments in the present section concern the algorithm for computing these samples.

The production of samples of the posterior density for the set (x, ξ, η) is not possible directly because of the complexity of the density. We therefore use a Gibbs algorithm [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF][START_REF] Gilks | Markov Chain Monte Carlo in practice[END_REF], which breaks the problem down into three simpler subproblems: sampling x, ξ, and η separately. This is an iterative algorithm, described in Table 1: each variable x, ξ, and η is drawn under its conditional posterior density given the current value of the other two variables. For each of the three steps, this conditional posterior density can be deduced directly (up to a multiplying factor) from the posterior density ( 12): all we have to do is to keep only the factors depending on the variable of interest. This algorithm is a Markov chain Monte-Carlo (MCMC) algorithm [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF][START_REF] Gilks | Markov Chain Monte Carlo in practice[END_REF] and is known to give (after a certain time, called the burn-in time) samples under the posterior density.

The conditional density of the map coefficients x (Step (1), Table 1) is Gaussian (Sect. 4.1). For the precisions ξ (Step (2), Table 1), the conditional densities are gamma densities (Sect. 4.2). They will be sampled using standard existing numerical routines (e.g. in Matlab). In contrast, the conditional density of the instrument parameters η (Step (3), Table 1) has a much more complex nonstandard form, so that it cannot be directly sampled by existing routines. To overcome this difficulty, sampling was carried out by means of a Metropolis-Hastings step (Sect. 4.3).

Map sampling

The density for the map x conditionally on the other variables is deduced from ( 12) by extracting the factors depending on x:

p(x|y, γ x , γ n , η) ∝ exp 1 2 γ x x t Π x x + γ n (y -m y| * ) t Π n (y -m y| * ) . ( 15
)
Considering the expression for m y| * given by Eq. ( 7), the argument of this exponential is quadratic in x. We deduce that we have a Gaussian density and, by rearranging the argument, we can determine the covariance and the mean

Σ x| * = γ n A t η Π n A η + γ x Π x -1 (16)
m x| * = γ n Σ x| * A t η Π n y, (17)
Remark 3. For a fixed value of the hyperparameters and the instrument parameters, the map m x| * defined by Eqs. ( 16), ( 17) is the maximizer (and the mean) of the conditional posterior density (15). This is the regularized least-squares solution denoted x(µ) parameterized by µ = γ x /γ n . This corresponds to the solution defined in our previous paper (Orieux et al. 2012b). For a convolutive instrument model, it is the Wiener solution (also called Wiener-Hunt solution; Orieux et al. 2010b).

Step (1) of Table 1 consists of sampling this Gaussian but this operation is made very difficult by three elements: (i) the large size of the map; (ii) the correlation introduced by the instrument model and the prior density; and (iii) the absence of structure of the instrument model (invariance, sparse nature). This problem can be solved using several approaches. For a convolutive instrument model (2), it is possible to approximately diagonalise the correlation matrix by FFT, thus producing a sample for the cost of an FFT (Chellappa & Chatterjee 1985;[START_REF] Chellappa | Markov Random Fields: Theory and Application[END_REF]Geman & Yang 1995;[START_REF] Giovannelli | [END_REF]Orieux et al. 2010b). If where the inverse of the correlation matrix is sparse, a partially parallel Gibbs sampler may be particularly efficient (Winkler 2003, Chap. 8). In the present case, neither the correlation nor its inverse possess the required properties. A general solution relies on factorizing the correlation matrix (Cholesky decomposition, diagonalization, etc.) but the large size of the matrix (M × M with M ∼ 10 5 ) does not permit the required calculations to be performed here.

The proposed solution consists of constructing a criterion such that its minimizer is a sample under the desired posterior conditional density. To do this, we perturb the means of the noise component and of the map component by an additive component with covariance γ x Π x and γ n Π n . A perturbed regularized least squares criterion is then introduced

J(x) = γ n y -A η x t Π n y -A η x + γ x x -m x t Π x (x -m x ),
and it can be shown (see Orieux et al. 2012a) that its minimizer 18) is Gaussian and does indeed have the correlation and mean defined by ( 16) and ( 17). This very powerful result has already been used by (Féron 2006;Orieux et al. 2012b). In different forms, similar ideas have been introduced and used by (Rue 2001;[START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF]Lalanne et al. 2001;[START_REF] Tan | Proc. IEEE ICASSP[END_REF].

x = γ n A t η Π n A η + γ x Π x -1 × γ n A t η Π n y + γ x Π x m x (
Remark 4. For the non perturbed criterion ( y = y and m x = 0), we have the regularized least-squares solution Eqs. ( 16), ( 17), that was mentioned in Remark 3.

Remark 5. The approach described in Orieux et al. (2012a) involves the sampling of the prior density (9) that is not properly defined here: the matrix Π x does not penalise the mean of the map (it is of deficient rank). But, for the same reason, the solution (18) does not depend on the mean of the realization of the prior density. Therefore, the simulated sample can have an arbitrary mean value.

Hyperparameter sampling

To determine the posterior conditional density for γ x , we examine the posterior density (12), and only keep the factors where γ x appears, which gives

p(γ x |y, x, γ n , η) ∝ γ M/2-1 x exp - γ x 2 ||x|| 2 Π x ,
and we recognize a Gamma density (see Appendix A)

γ x ∼ G M/2 , 2/||x|| 2 Π x . (19)
Concerning γ n , we also refer to the posterior density (12) and find

p(γ n |y, x, γ x , η) ∝ γ N/2-1 n exp - 1 2 γ n ||y -m y| * || 2 Π n ,
which is also a Gamma density

γ n ∼ G N/2 , 2/||y -m y| * || 2 Π n . (20)
For both γ x and γ n the second parameter of the Gamma density introduces a quadratic norm (regularity of the map in Eq. ( 19), and goodness-of-fit in Eq. ( 21)), which can be easily computed.

Remark 6. An intuitive interpretation can be given to these results starting from the fact that the mean of the Gamma density is equal to the product of its parameters (see Appendix A), here N/ ym y| * 2 Π n for (21). In this sense, the conditional posterior mean is the inverse of the empirical variance of the residuals. Consequently, when the goodness-of-fit term is small, the mean of the density is large and so the sampled value of γ n is also high reporting a high precision, i.e. a weak variance (and vice versa). The same holds for the map regularity in relation with the mean of the density (19) given by M/ x 2 Π x . These observations support the coherence of the model and reinforce the prior choice for these hyperparameters as a convenient one. 1. The current sample at step q is η (q) and it is either replaced or not by the proposed sample η p .

(a) Draw a sample η p under a proposal density. (b) Compute the acceptation ratio ρ by Eq. ( 23). (c) Replace η (q) by η p (i.e. η (q+1) = η p ) with the probability min(1, ρ), otherwise keep η (q) (i.e. η (q+1) = η (q) ).

Regarding the mean of the noise, m n , it is a scalar whose posterior conditional density is also deduced from the posterior density (12) and from ( 7)

p(m n |y, x, γ n , γ x , η) ∝ exp - 1 2 γ n (m n -m r ) 2 ,
where m r is the empirical mean of the residuals y -A η x. We then have a Gaussian density

m n ∼ N (m r , γ n ) . (21)
In the numerical evaluation of Sect. 5, one such mean is estimated for each bolometer.

Remark 7. If we examine the relationships above, we see that p(γ x |y, x, γ n , η) = p(γ x |x), in other words, γ x and (y, γ n , η) are independent conditionally on x. Similarly, we note that p(γ n |y, x, γ x , η) = p(γ n |y, x, η), which means that γ n and γ x are independent conditionally on (y, x, η). In addition, m n is independant of γ x , given y, x, γ n , η.

Instrument parameter sampling

The last step (Step (3) of Table 1) is more complex. As for the other variables, the posterior conditional density can be deduced from the posterior density (12) by keeping the factors that bring in η. There are two of these: the likelihood and the prior density, and we thus have

p(η|y, x, ξ) ∝ K k=1 exp - 1 2 (η k -µ k ) 2 ρ k × exp - 1 2 γ n (y -A η x) t Π n (y -A η x) . (22)
However, this is not a usual density, notably because there is no simple mathematical form to represent the dependence of the observation w.r.t. η. Thus, Step (3) of Table 1 cannot be carried out directly with standard sampling routines and we resort to a Metropolis-Hastings step in a random-walk version [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF][START_REF] Gilks | Markov Chain Monte Carlo in practice[END_REF] described in Table 2. It can be briefly explained as follows. Because it is impossible to draw a sample directly under the conditional posterior density ( 22), a sample is drawn under another density (namely the proposal density), but is not systematically accepted. Acceptance or rejection is also random with a precisely defined probability (see Eq. ( 23)) to ensure that, at convergence, we have samples under the target density [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF][START_REF] Gilks | Markov Chain Monte Carlo in practice[END_REF]. The algorithm is divided into three sub-steps summarized in Table 2 and detailed here.

(a) Draw a proposal η p as a perturbation of the current value: η p = η (q) + ε, deduce the instrument matrix A η p , and the corresponding model output m p y = A η p x (q) + m n .

(b) Compute the acceptation ratio ρ = p(η p |y, x, ξ) p(η (q) |y, x, ξ) ,

based on the conditionnal posterior law ratio that compares the goodness-of-fit for the current parameter and the proposed one. (c) Accept or reject the proposal, at random, with probability min(1, ρ). To do so, draw u uniformly in [0, 1] and take q) otherwise.

η (q+1) = η p if u < min{1, ρ} η (
These three substeps are inserted instead of Step (3) of Table 1.

The algorithm can be explained as follows. Starting with a current value η (q) , the algorithm proposes a new value η p and compares the goodness-of-fit for the two values. When the proposed value improves the fit, ρ > 1 and η p is accepted. When the proposed value degrades the fit, η p can be accepted or rejected, with a probability that is higher or lower, depending on how weak the degradation is.

Remark 8. There are other more complex (and potentially more efficient) approaches for Metropolis-Hastings sampling. In particular, the proposal density can be adapted, e.g. directional random walk [START_REF] Vacar | Proc. IEEE ICASSP[END_REF]. They are not exploited here but are considered in the development perspectives.

Numerical results

The previous sections presented the approach for building the posterior density and for its exploration by stochastic sampling using a Gibbs algorithm including a Metropolis-Hastings step. The mean and standard deviation (SD) of the posterior density are numerically computed as empirical averages based on simulated samples, from relations ( 13) and ( 14). The developments below show the practicability of the proposed method (models, estimate and algorithm), and provide a first numerical evaluation.

Evaluation methodology

The evaluation is based on the SPIRE instrument (Griffin et al. 2010) of the Herschel Space Observatory (Pilbratt et al. 2010) launched in May 2009. It focuses on the PMW channel (centred around 350 µm) and the Large Map protocol in the nominal operating conditions: scan back and forth with constant speed (30 ′′ /s) over two almost perpendicular directions (88 • ). The scans are associated with a high sampling frequency (F s ≈ 30 Hz) providing spatially redundant observation and Fig. 2 shows the corresponding redundancy/pointing map. The spatial shift between basis functions (see Eq. ( 1)) is fixed at δ α = δ β = 2 ′′ , based on our earlier work [START_REF] Orieux | Proc. IEEE ICIP[END_REF](Orieux et al. , 2012b)), to obtain the best gain in resolution without important increase of the computational cost. The angular size of the reconstructed map is 20 ′ × 20 ′ , i.e. a map of 600 × 600 coefficients. The associated direct model, including the whole acquisition chain (scanning strategy, mirror, horns, wavelength filters, bolometers, and electronics) is detailed in our previous paper (Orieux et al. 2012b) and represented by Eq. ( 2) of the present paper.

The unsupervised method is assessed based on two synthetic maps of extended emission (the Galactic Cirrus (Fig. 3e) and a realization of the prior density Eq. ( 9)) as well as based on a real observation (reflection nebula NGC 7023, Fig. 7). The paper also proposes a first assessment of the unsupervised and myopic approach based on a synthetic map with broader spectral content (the Galactic Cirrus with point sources (Fig. 3e)).

In the simulated cases, a zero-mean white Gaussian noise is added to the model output. Moreover, in these cases, since the original map (the "sky truth" denoted by x * ) is known, the quality of the reconstruction (denoted by x) can be quantified through an error:

E = i, j |x * i j -x i j | 2 i, j |x * i j | 2 , ( 24
)
where only coefficients in the observed area are taken into account, allowing assessments and comparisons between methods.

Algorithm behaviour and general comments

As explained in Sect. 4, the algorithm provides a series of samples that form a Markov chain for hyperparameters, instrument parameter and map. The MCMC theory then ensures that it correctly explores the parameter space and produces a density of samples reflecting the posterior density. Practically, the algorithm has been executed for the unsupervised problem as well as for the unsupervised and myopic problem. It has been run several times (1) using identical initial conditions and (2) using different initial conditions. In both cases, the same qualitative and quantitative behaviour as presented here has been systematically observed. The computation time takes about one hour for the unsupervised (nonmyopic) case and about ten hours for the unsupervised and myopic case. The main computational cost is due to computating the instrument model output given by Eq. ( 2).

Figures 678present some typical elements of the algorithm operation: visualization of progression, convergence phase (burn-in period), stable phase, etc. The evolution of the chain is shown for hyperparameters (Figs. 6 and7) and for instrument parameter (Fig. 8). It is thus possible to grasp how the parameter space is explored.

Unsupervised approach

Assessment of map estimation

The qualitative and quantitative assessment of the reconstructed maps is presented here for the Galactic Cirrus; the first results are shown in Fig. 3.

-The unsupervised method (proposed method) is outlined in Fig. 3a. The hyperparameters are automatically set (without knowing of the sky truth). 3a), best map (Fig. 3b), naive map (Fig. 3d), and true map (Fig. 3e). Figure 3c shows a profile (marked by the white line in Fig. 3e) and Fig. 3f the spectrum (circular means of power spectra). Uncertainties are given in Fig. 4 and quantitative results are given in Table 3. Comments are given in Sect. 5.3.1.

-The best-supervised method is outlined in Fig. 3b. The hyperparameters are set by hand to minimize the error E (knowing the sky truth). It is referred to as the best map and was previously presented in our paper (Orieux et al. 2012b). -The naive map (coaddition) and the true map are shown in Figs. 3d ande. -In addition, Fig. 3c gives spatial profiles (vertical in the middle of the map) and Fig. 3f gives the spectral4 profiles.

As expected and shown in (Orieux et al. 2012b), the inversion based on an accurate instrument model considerably improves the quality of the map: see the proposed map and the best map compared to the naive map and to the true map. The proposed map is visually very similar to the true map. In particular, our method restores details of small spatial scales (with spectral extension from null to high frequency) that are invisible on the naive map but are present on the true map (see Fig. 3c). In addition, our method correctly restores the structures of large spatial scales (low frequencies) and also the mean level of the map (null frequency), i.e. the photometry.

To assess the pertinence of estimating γ n and γ x in terms of map quality, we compared the proposed map with the best map. They are visually very similar (see Figs. 3a andb). In the quantitative terms given by Table 3, the best map produces an Table 3. Comparison of reconstruction error E (see Eq. ( 24)) for the Galactic Cirrus (the error E only accounts for the observed area of the map).

Reconstruction error E Unsupervised ( γ x ) 0.016% Best supervised (γ best x ) 0.0129% Naive map 0.0435%

Notes. The proposed approach (which does not require knowing the true map) produces an error only very slightly higher than the best map (which does require knowing true map).

error E of 0.0129% and the proposed map produces an error E of 0.016%, which is only slightly higher. In other words, the proposed unsupervised method automatically (without knowing the sky truth) determines hyperparameters that produce a map almost as good as the best map (which requires knowing the sky truth). However, the proposed map shows a fine grainy texture that is visible neither in the true nor in the best map. This feature is also visible on the residual map of the coefficiens (Fig. 5). This is also observable in Fig. 3f: the spectrum of the proposed map passes above the spectrum of the true map in the spectral band 0.025-0.035 arsec -1 . This defect is related to a slight overevaluation of the observation contribution with respect to the prior contribution. It is referred to as under-regularization and yields an overamplification of the observation in this spectral band. This confirms the behaviour previously observed in deconvolution (Orieux et al. 2010b;[START_REF] Giovannelli | [END_REF] or noted for the maximum likelihood (Fortier et al. 1993). Nevertheless, it is remarkable that this defect is correctly notified by the PSD, as explained in the next paragraph. Indeed, the approach naturally provides a measure of reliability through the PSD shown in Fig. 4. Two zones can be seen in Fig. 4a, in accordance with Fig. 2: the central zone (where observations are available) and the peripheral zone (extrapolated from observations of the central zone and based on the prior regularity). The boundary between the two zones also exhibits the variation of the observation hit and scanning strategy notably well. In addition, the posterior standard deviation also illustrates the difference between the zones observed with our without cross scan. From a spatial standpoint, Fig. 4b shows an interval around the estimated map with plus/minus PSD. The main result is that the true map is clearly within the interval. In a similar way, from a spectral standpoint the results are given by Fig. 4c (in relation with Fig. 3f): the true spectrum is also within the interval. More specifically, incorrectly reconstructed in the spectral band (above 0.025 arcsec -1 ), the stronger PSD clearly shows that the estimated spectrum is certainly submitted to unsatisfactory errors or confidence.

Assessment of hyperparameter estimation

This section assesses the unsupervised capabilities through evaluating the hyperparameter estimation using the Galactic Cirrus and a realization of the prior for the map (which makes a true value γ * x available). Fig. 6 shows the chains and the histograms that approximate the marginal posterior densities p(γ n |y) and p(γ x |y). In both cases, the histogram is relatively narrow although the prior is a wide non-informative Jeffreys' distribution (see Eq. ( 10)). In other words, the observations are sufficiently informative to quantify noise and regularity level and the method is able to capture this information.

From a quantitative standpoint, results are given in Table 4. For the Galactic Cirrus and prior realization, the estimated values γ n are very similar to the true value γ * n (error is less than 1%). Moreover, the PSD are very low (0.40%). In the case of prior realization, the estimated value γ x is in the correct range but the error is larger (about 17%) and the PSD is 1.7%. This difference can be naturally explained by two elements: (i) the noise is added at the system output, so it is directly observed, whereas the map is at the system input, i.e. indirectly observed; (ii) the added noise is a realization of the prior density for the noise while the Galactic Cirrus is not a realization of the prior density for the map.

Real observation processing

This section proposes a first assessment for a real observation. It is based on the reflection nebula NGC 7023 acquired during the science demonstration phase of Herschel, which as been presented in (Abergel et al. 2010) and was processed in our previous paper (Orieux et al. 2012b). There and here, computations are made on the level-1 files processed using HIPE. In our previous paper (Orieux et al. 2012b), the offsets were removed in a pre-processing step and the regularization parameter was tuned by hand compromise between gain in resolution and overamplification of the observations. In contrast, here both are automatically tuned.

Figure 7 presents the evolution of the chains for the hyperparameters: Fig. 7a for the noise parameter γ n and Fig. 7b for the image parameter γ x . It is important to notice that the algorithm behaves in a very similar manner for the real observation and for the simulated observation (see Fig. 7 compared to Fig. 6). The figures also give an empirical indication of the algorithm operation: after a burn-in time (empirically less than about 500 iterations) the stationary state is attained and the chain remains in a steady state: the samples are drawn under the posterior density. Concerning the offsets, the chains begin in the steady state, thanks to a good initialization based on (Orieux et al. 2012b) results. All bolometer offsets behave in the same manner with two example illustrated Figs. 7d ande. The empircal mean of the offsets sample start to stabilize at approximately 500 samples.

Figure 7c shows the corresponding reconstructed map. Its quality is equivalent to the quality of the map restored by empirically tuning the hyperparameter presented in Orieux et al. (2012b), Fig. 8. In other words, the proposed unsupervised method automatically determines hyperparameters (noise power and offsets as well as sky power) that produce a map almost as good as the map produced by a hand-made hypermarameter tuning. In addition, the map remains far better than the naive map shown in Fig. 7f.

Myopic and unsupervised approach

The myopic and unsupervised question is a threefold problem that is much more ambitious: estimate the instrument parameter, the hyperparameters, and the map itself from a unique observation. In addition, the instrument parameter intervenes in a complex way in the description of the observations, and moreover, the problem is stated in a context that is doubly delicate: ill-posedness and high-resolution.

In Orieux et al. (2012b), the equivalent PSF has a Gaussian shape whose standard deviation is proportional to the wavelength: σ o (λ) = ηλ. It is then integrated w.r.t. the wavelength (to include the spectral extend) and w.r.t. the time parameter (to account for the bolometer response) to form the global instrument response. To test the method, we consider the instrument parameter η to be poorly known and introduce elements of the feasibility to estimate it.

The prior is the Gaussian density given by Eq. ( 11), with K = 1. Its mean is taken from the SPIRE observer manual µ = 2.96 × 10 4 [ ′′ /m] and its standard deviation is set to ρ = 10 4 , i.e. a relatively large uncertainty. It is about 33% of the mean and an equivalent prior interval is [0.96 × 10 4 , 4.96 × 10 4 ] in a two-standard-deviations sense. Two cases are investigated for the true value (used to simulated observations): η * 1 = 2.46 × 10 4 and η * 2 = 3.46 × 10 4 . The conditional posterior for η (Sect. 4.3) does not have a standard form and its sampling (step (3) of Table 1) relies on a Metropolis-Hastings sampler, itself based on a random-walk with a Gaussian excursion. The size of the excursion was chosen so that the acceptation rate is around 50%.

Two maps are used for the observation: the Galactic Cirrus and the Galactic Cirrus with point sources. In each case, the algorithm was run several times from identical and different initializations, and shows similar qualitative and quantitative behaviours as those in Fig. 8.

Nevertheless, as expected, the spectral content of the Galactic Cirrus is not sufficiently extended towards high frequencies to provide an excitation that is adequate for instrument identification. In contrast, the Galactic Cirrus with point sources is more extended and estimations are more accurate. Table 5 presents quantitative assessments. The main result is that the estimation error is about 6%. It is a remarkable result given the difficulty of the problem (triple problem, complex relations, ill-possedness, and high resolution) and given that the prior uncertainty is about 33%. In other words, the method is able to capture information about instrument parameter, jointly with noise level, regularity level, and map from a unique observation. However, the parameter η seems to be slightly underestimated which, we explain as follows. The input map (with point sources) presents a broad spectral extent whereas the prior favours spatially extended maps (dominated by relatively low frequencies), so the posterior advocates a narrower PSF to compensate for this spectral discrepancy.

Figures 9a andb show the related maps. They must be compared to the map restored with the true instrument parameter and the best hyperparameter presented in Fig. 7b of Orieux et al. (2012b) and in Fig. 9d here. They must also be compared to the true map and the naive map also given in Orieux et al. (2012b) and in Figs. 9c-e here. As previously, the proposed maps show a fine grainy texture but despite this defect, they remain similar to the true map. The quality of the proposed maps is similar to the quality of the map restored with the true instrument parameter and the best hyperparameter. In addition, several point sources of the true map are visible on the proposed maps but not on the naive map. In other words, the proposed method automatically determines instrument parameter and hyperparameters that produce a map almost as good as the best one and better than the naive map.

Conclusion

We described regularized methods for image reconstruction and focused on parameter estimation:

hyperparameters, which guide the trade-off between priorbased and observation-based information; -instrument parameter, which tunes the physical characteristics of the model of the acquisition system.

They were jointly estimated with the map of interest. We were therefore dealing with an unsupervised and myopic inverse problem.

The most delicate point is jointly handling the different types of variables and their interactions in direct terms but, above all, in inverse terms. From a methodology point of view, we worked in the framework of hierarchical full Bayes strategies that model the available information for each set of variables (map, hyperparameters, instrument parameter, and observations) under a probability density. We defined the posterior density, which gathers the information on the map of interest and the parameters, given the observations. We then defined the posterior mean as an estimate of the map and the posterior standard deviation as a measure of uncertainty, which gives an uncertainty map. This approach makes it possible to work in a global and consistent framework to solve the problem as a whole. It draws its inspiration from our earlier works on deconvolution (Orieux et al. 2010b) and adapts them to the case at hand. The posterior density was explored by stochastic sampling using a Gibbs algorithm. The sampling of the map was difficult: we are dealing with a large-sized multivariate normal density for which classical techniques do not apply. We overcame this difficulty by constructing a sample as the minimizer of a well-chosen perturbated criterion (Orieux et al. 2012a). Another problematic point is the instrument parameter sampling: we are dealing with a very complex, nonstandard density. This difficulty was overcome by means of a Metropolis-Hastings step. The estimate of the map as well as the parameters (posterior mean) and the uncertainties (posterior standard deviation) were calculated numerically as empirical averages based on the simulated samples.

We presented a first application of the developments (Bayesian estimation method and stochastic sampling algorithm) in a real context: the SPIRE instrument of the Herschel Space Observatory. The study was essentially performed on simulated observations and has also yielded some initial results on real observations. We concluded that the approach is applicable and enables joint estimation of the map, the hyperparameters, and the instrument parameter from a unique observation. We showed, among other results, that the quality of the proposed map is similar to that obtained when the instrument parameter is known and the hyperparameters are fixed by hand in a supervised way (using the sky truth). The method shows remarkable results given the difficulty of the problem. It seems to us that these initial results are particularly promising and worth developing. They may open up many new perspectives for imaging in astrophysics in a myopic and unsupervised framework.

Fig. 1 .

 1 Fig. 1. Graphical dependency representation (hierarchical structure). The round (square) nodes correspond to unknown (fixed) quantities. The directions of the arrows indicate the dependencies.

Fig. 2 .

 2 Fig. 2. Redundancy/pointing map associated with our experiment of five crossed scans.

Fig. 3 .

 3 Fig.3. Comparison of reconstructed map for the Galactic Cirrus: proposed map (Fig.3a), best map (Fig.3b), naive map (Fig.3d), and true map (Fig.3e). Figure3cshows a profile (marked by the white line in Fig.3e) and Fig.3fthe spectrum (circular means of power spectra). Uncertainties are given in Fig.4and quantitative results are given in Table3. Comments are given in Sect. 5.3.1.

Fig. 5 .

 5 Fig.4. PSD and quantification of uncertainties. Fig.4ashows the map of the PSD and Figs.4b and cshow the interval around the estimated map ±PSD and the true map. In the spatial domain, Fig.4bis a profile (marked by the white line on 3e) and in the spectral domain Fig.4cis the circular means of the power spectra.

Fig. 6 .

 6 Fig.6. Chains and histograms for γ n , Fig.6a, and γx , Fig.6b, for the Galactic Cirrus. The chains show the burn-in period (about 1000 iterations) and the steady state. The corresponding histograms are computed on steady state only.

  Fig.7. Results for real observation processing (reflection nebula NGC 7023). Chains for the noise parameter (γ n ) and for the image parameter (γ x ) in Figs.7a and b. The stationary state is attained after a burn-in time of about 500 iterations. Figure7cshows the corresponding map. Figures7d and eillustrate chains and marginal histogram for two bolometer offsets.

Fig. 9 .

 9 Fig. 9. Restoration of cirrus superimposed on point sources. The proposed maps must be compared to the maps restored with the true instrument parameter and the best hyperparameter and with the naive map.

Table 2 .

 2 Step of the Metropolis-Hastings sampler, which replaces Step (3) of Table

Table 4 .

 4 Hyperparameter estimation: true values, estimates, PSD and best values. Cirrus 10 6 1.009 × 10 6 4.07 × 10 3 -8.99 × 10 11 2.46 × 10 10 2.47 × 10 12 Prior 10 6 1.003 × 10 6 4.05 × 10 3 4 × 10 11 3.28 × 10 11 1.07 × 10 10 8.37 × 10 11

	γ * n	γ n	PSD	γ * x	γ x	PSD	γ best x

Table 5 .

 5 Quantitative evaluation of the estimation of the instrument parameter using the Galactic Cirrus with point sources.Notes. Prior mean and standard deviation are µ = 2.96 × 10 4 and ρ = 10 4 . Instrument parameter chain (myopic and unsupervised approach) for the of Galactic Cirrus with point sources. Left (right) part of the figure deals with case 1, i.e. η * 1 = 2.46 × 10 4 (case 2, i.e. η * 2 = 3.46 × 10 4 ). The horizontal axis gives the iteration index and the vertical range is the prior interval in a two-standard-deviations sense. The true value is shown by the straight line.

						Case	η *		η	η -η *	(η -η * )/η *	σ
						1			2.46 × 10 4 2.29 × 10 4 -1.65 × 10 3	6.7%	2.2 × 10 2
						2			3.46 × 10 4 3.27 × 10 4 -1.86 × 10 3	5.4%	2.9 × 10 2
	6	x 10	4		Chain		6	4 x 10			Chain
					True						True
	5				Mean		5				Mean
					Mean+-sigma					Mean+-sigma
	4						4
	3						3
	2						2
	0	500	1000	1500	2000	0	500	1000	1500	2000
		(a) Case 1, parameter chain		(b) Case 2, parameter chain
	Fig. 8.

The mean square error is the expected value of the squared norm of the difference between estimated value and true value. The expectation is under the distribution of the observation and the unknown. The MSE is the sum of the variance and the squared bias of the estimator (under the distribution of the observation and the unknown).

The function squared norm is defined by u 2 = u(α, β) 2 dα dβ.A83, page

of 13 A&A 549, A83 (2013)

It yields a non-informative prior distribution based on a key feature that it is invariant under reparameterization. It is deduced as the determinant of the Fisher information matrix.

This spectrum is computed from the FFT-2D of the map by averaging the coefficients in regularly spaced concentric rings. This gives a 1D spectrum containing the isotropic approximation of the spectral map properties.

Appendix A: Gamma probability density

The gamma pdf for γ > 0, with given parameters a > 0 and b > 0, is written

The following properties hold: mean is

= ab 2 and maximizer is b(a -1) if and only if a > 1.