Alexander Mendler
email: alexander.mendler@ubc.ca

Michael Döhler

Carlos E Ventura

Sensor placement with optimal damage detectability for statistical damage detection

Keywords: Optimal sensor placement, structural health monitoring, ambient vibrations, Gaussian feature, statistical damage detection, Fisher information

Damage diagnosis based on global structural vibrations critically depends on the sensor layout, in particular when a small number of sensors is used for large structures under unknown excitation. This paper proposes a sensor placement strategy that yields an optimized sensor layout with maximum damage detectability in selected structural components. The optimization criterion is based on the Fisher information, which quantifies the information that the damage-sensitive feature carries on the design parameters of structural components, such as material constants or cross-sectional values. It is evaluated using a finite element model, and considers the statistical uncertainties of the damage-sensitive feature. The methodology is shown for the stochastic subspace-based damage detection method, but can be applied to any damage-sensitive feature whose distribution can be approximated as Gaussian. It is suitable to find the optimal layout for a fixed number of sensors and to choose an appropriate number of sensors. Since the Fisher information is defined component-wise, the sensor layout can be tuned to become more sensitive to damage in local structural components, such as damage hotspots, non-inspectable components, or components that are critical for the safety and serviceability of the structure. For proof of concept, the sensor layout on a laboratory beam is optimized based on numerical simulations, and it is showcased that the optimal sensor layout leads to the highest damage detectability for experimental data.

Introduction

Modern societies are heavily dependent upon structural and mechanical systems, such as bridges, power plants, off-shore platforms, or wind turbines. To guarantee the safety and serviceability of these structures, periodically scheduled on-site inspections are conducted by trained personnel, and structural health monitoring (SHM) systems are installed. This process of SHM involves the observation through permanently installed sensors, the extraction of damage-sensitive features for data reduction, and their subsequent statistical evaluation using smart algorithms [START_REF] Farrar | Structural health monitoring: A machine learning perspective[END_REF]. This way, structural damage can be detected online, triggering on-demand inspections after sudden and unexpected damages or extreme events, such as storms, floods, earthquakes, or impacts of any sort. In particular, output-only vibration-based methods based on the global system response to unknown excitation have been actively developed, because they can be applied under normal operating conditions and with relatively few vibration sensors. It goes without saying that the success of such methods heavily depends on a strategic sensor placement.

In recent literature reviews [START_REF] Yi | Methodology developments in sensor placement for health monitoring of civil infrastructures[END_REF][START_REF] Ostachowicz | Optimization of sensor placement for structural health monitoring: A review[END_REF][START_REF] Barthorpe | Emerging trends in optimal structural health monitoring system design: From sensor placement to system evaluation[END_REF], numerous optimization criteria can be found to optimize the sensor placement for output-only methods. Low-level performance criteria aim to precondition the signal and optimize the signal-to-noise ratio. The measured signal contains information on natural frequencies, damping ratios, and mode shapes. To reduce the noise components, the mode shape observability and amplitude are used as performance criteria, including metrics such as the mode shape summation plot [START_REF] De Clerck | Development of several new tools for pre-test[END_REF], the average driving point residue [START_REF] Chung | On-orbit sensor placement and system identification of space station with limited instrumentations[END_REF], the eigenvalue vector product [START_REF] Larson | A comparison of modal test planning techniques: Excitation and sensor placement using the NASA 8 bay truss[END_REF], or the kinetic energy per mode [START_REF] Kammer | Sensor placement for on-orbit modal identification and correlation of large space structures[END_REF]. The redundancy of information is also considered as a low-level criterion and could be quantified through the mutual information criterion [START_REF]Encyclopedia of structural health monitoring[END_REF]. More advanced criteria aim to optimize the quality of the modal system identification. On one hand, this second group of criteria aims to maximize the identifiability and distinguishability of modes quantified through the modal assurance criterion [START_REF] Allemang | A correlation coefficient for modal vector analysis[END_REF], the modal clarity criterion [START_REF] Soman | Multi-type, multi-sensor placement optimization for structural health monitoring of long span bridges[END_REF], and the singular value decomposition ratio (i.e. the condition number) [START_REF] Yi | Methodology developments in sensor placement for health monitoring of civil infrastructures[END_REF]. On the other hand, this group of methods attempts to reduce the uncertainty in the estimation through measures such as the information entropy [START_REF] Papadimitriou | Optimal sensor placement methodology for parametric identification of structural systems[END_REF][START_REF] Yuen | Optimal sensor placement methodology for identification with unmeasured excitation[END_REF][START_REF] Papadimitriou | The effect of prediction error correlation on optimal sensor placement in structural dynamics[END_REF], the covariance of the state prediction error during Kalman filtering [START_REF] Kumar | Optimal location of measurements for distributed parameter estimation[END_REF][START_REF] Tzoumas | Sensor placement for optimal kalman filtering: Fundamental limits, submodularity, and algorithms[END_REF], and the Fisher information on the modal contribution of sensor locations [START_REF] Kammer | Sensor placement for on-orbit modal identification and correlation of large space structures[END_REF]. In the latter case, optimization criteria are formed based on the matrix determinant [START_REF] Kammer | Sensor placement for on-orbit modal identification and correlation of large space structures[END_REF], the minimum singular value [START_REF] Reynier | Sensors location for updating problems[END_REF], the singular value decomposition ratio [START_REF] Kim | Optimal sensor placement for dynamic testing of large structures[END_REF], or the absolute values of the main diagonal terms. Criteria such as the kriging prediction error [START_REF] Papadimitriou | Optimal experimental design in stochastic structural dynamics[END_REF] and the representative leastsquare [START_REF] Li | On optimal sensor placement criterion for structural health monitoring with representative least squares method[END_REF] moreover allow for the minimization of uncertainties at unmeasured degrees of freedom. Recent literature reviews emphasize that there is a need to develop method-specific performance criteria that are tailored to the employed damage diagnosis method [START_REF] Yi | Methodology developments in sensor placement for health monitoring of civil infrastructures[END_REF][START_REF] Barthorpe | Emerging trends in optimal structural health monitoring system design: From sensor placement to system evaluation[END_REF]. Such high-level criteria that target the detectability of damage directly are rare [START_REF] Shi | Optimum sensor placement for structural damage detection[END_REF][START_REF] Lin | Experimental investigation on multi-objective multi-type sensor optimal placement for structural damage detection[END_REF], and typically require supervised training using numerical data from numerous damage states [START_REF] Guratzsch | Structural health monitoring sensor placement optimization under uncertainty[END_REF]. Moreover, none of the existing criteria appear to consider the detectability of changes in selected structural components, which is the focus of this paper. Considering that optimal damage detectability in all components with a limited number of sensors is unrealistic, the goal is to develop an optimization criterion that targets a prioritized set of damage scenarios. Such damage hotspots are typically known from the operational evaluation of the structure and should be considered during sensor placement.

Finding the optimal sensor layout is not straightforward, regardless of the chosen optimization criterion discussed above. In most cases, the number of possible sensor combinations is excessively large, and it is infeasible to evaluate all combinations in an exhaustive search (ES). To remedy this, smart optimization methods are applied. Most performance criteria for sensor placement are non-differentiable, meaning the relation between input and output cannot be described through analytical functions, and deterministic optimization methods (which evaluate the gradient or Hessian matrix) are rarely used [START_REF] Yi | Methodology developments in sensor placement for health monitoring of civil infrastructures[END_REF]. Sequential sensor placement (SSP) methods, where the optimal sensor configuration is found by systematically adding or removing a subset of channels in each iteration step [START_REF] Li | Extended minmac algorithm and comparison of sensor placement methods[END_REF], are efficient, at least for single-objective optimization [START_REF] Papadimitriou | Optimal sensor placement methodology for parametric identification of structural systems[END_REF]. However, the solution depends on the initial sensor choice and the number of sensors that are removed at the time, and there is no guarantee that the optimal solution can be found [START_REF] Friswell | Clustering of sensor locations using the effective independence method[END_REF]. Stochastic optimization methods are more appropriate, as they combine local search with coordinated global search techniques, and randomized local and global elements to overcome local minima. One of the most prominent groups are genetic algorithms (GAs) [START_REF] Holland | Adaptation in natural and artificial systems[END_REF]. Other methods include simulated annealing [START_REF] Van Laarhoven | Simulated annealing: Theory and applications[END_REF][START_REF] Corana | Minimizing multimodal functions of continuous variables with the simulated annealing algorithm[END_REF], particle swarm optimization [START_REF] Eberhart | A new optimizer using particle swarm theory[END_REF], the distributed wolf method [START_REF] Yi | Multiaxial sensor placement optimization in structural health monitoring using distributed wolf algorithm[END_REF], the ant colony optimization method [START_REF] Overton | Sensor optimisation using an ant colony metaphor[END_REF], the artificial bee colony method [START_REF] Scott | A bee swarm algorithm for optimising sensor distributions for impact detection on a composite panel[END_REF], the firefly algorithm [START_REF] Zhou | Sensor placement optimization in structural health monitoring using clusterin-cluster firefly algorithm[END_REF], the micro-habitat frog leaping method [START_REF] Feng | Acceleration sensor placement technique for vibration test in structural health monitoring using microhabitat frog-leaping algorithm[END_REF], and the monkey algorithm [START_REF] Yi | Health monitoring sensor placement optimization for canton tower using immune monkey algorithm[END_REF], with an overview in [START_REF] Ostachowicz | Optimization of sensor placement for structural health monitoring: A review[END_REF][START_REF] Barthorpe | Emerging trends in optimal structural health monitoring system design: From sensor placement to system evaluation[END_REF].

The objective of this paper is to develop a sensor placement strategy that optimizes the detectability of damage in selected structural components, and to combine it with a genetic algorithm to efficiently find a close-to-optimal sensor layout for large structures. A performance criterion is derived for statistical damage detection tests and damage-sensitive features that are (approximately) Gaussian distributed, such as the stochastic subspace-based residual [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF][START_REF] Döhler | Fault detection, isolation and quantification from Gaussian residuals with application to structural damage diagnosis[END_REF][START_REF] Mendler | A reliability-based approach to determine the minimum detectable damage for statistical damage detection[END_REF]. The statistical test considers the uncertainties in the feature and the sensitivity toward changes in structural components, calculated based on a finite element (FE) model. Using the Fisher information [START_REF] Mendler | A reliability-based approach to determine the minimum detectable damage for statistical damage detection[END_REF][START_REF] Mendler | Minimum diagnosable damage and optimal sensor placement for structural health monitoring[END_REF], the information from the sensitivity analysis and the uncertainty estimation are combined and a performance criterion is developed, which takes as input the desired detectable damage in percent of local structural design parameters. Particular features of the approach are that the sensor layout can be optimized based on the undamaged structure and that it can be tuned to maximize the damage detectability in selected structural components, e.g. damage hotspots, or other components that are critical for the safety and serviceability of the structure. The resulting sensor placement strategy is suitable to find an optimal sensor layout, and can help to choose an appropriate number of sensors. Keeping in mind that no data is available from real structures before sensors are installed, the developed sensor placement strategy is validated based on numerical simulations and verified with experimental data of a laboratory beam. The results showcase that the optimal sensor layout leads to the highest damage detectability.

The paper is organized as follows: Section 2 gives background on the statistical damage detection methods for which the sensor layout will be optimized. Section 3 proposes the sensor placement strategy. In Section 4, the approach is applied to a pin-supported beam, where the sensor layout is optimized based on simulated data, with an experimental proof of concept study in Section 5. Ultimately, the paper is concluded in Section 6.

Statistical Damage Detection

Damage-sensitive Residual

Suppose that a consistent estimate of a damage-sensitive feature vector x = [x 1 , . . . , xG ] T is available from the measurements of a monitored system, computed from a measurement record of sample size N . Further suppose that damage manifests itself in a deviation of the feature mean; then, a damage-sensitive residual vector ε can be defined by normalizing the feature by its mean vector m = E(x) in the reference state, with

ε = x -m. (1)
If the underlying measurements are assumed to be the realizations of a stationary stochastic process, the vector ζ = √ N ε is often asymptotically normal distributed [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]

ζ = √ N ε -→ N (0, Σ) (reference state) (2)
where the mean vector is zero, and Σ is the covariance of the Gaussian residual vector ζ, describing the dispersion about the mean due to stochastic loads and measurement noise. In other words, the distribution of the residual can be approximated by a normal distribution when computed on sufficiently large sample size N -a criterion that is met by many features from the literature, e.g.

modal parameters [START_REF] Mellinger | Variance estimation of modal parameters from output-only and input/output subspace-based system identification[END_REF][START_REF] Oliveira | Vibration-based damage detection in a wind turbine using 1 year of data[END_REF] or features that are formed based on the subspace properties of output covariance Hankel matrices [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF][START_REF] Greś | Subspace-based Mahalanobis damage detection robust to changes in excitation covariance[END_REF][START_REF] Viefhues | Statistical subspace-based damage detection with estimated reference[END_REF]. Some features from the literature [START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF][START_REF] Viefhues | Fault detection for linear parameter varying systems under changes in the process noise covariance[END_REF][START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: from automated OMA to damage detection[END_REF] are even robust to changes in environmental conditions. In Appendix A, an exemplary Gaussian residual vector is given, i.e. the stochastic subspace-based residual, which will be used to demonstrate the developments in this paper.

Damage Parametrization and Detection

For meaningful damage detection, the data-driven features from the previous section have to be linked to monitoring parameters θ that characterize the structural system, e.g. design parameters of local structural components. This way, structural damage can be defined as a deviation of the monitoring vector from its reference vector, θθ 0 , and the effect on the damage-sensitive feature can be described through the first-order Taylor expansion

E θ [ε] ≈ E θ 0 [ε] + J (θ -θ 0 ), (3)
where

J = ∂E θ [ε]
∂θ θ=θ 0 is the first-order sensitivity matrix (the Jacobian matrix), and E θ is the expected value of a random variable calculated under system parameter θ. Due to the formulation as a residual vector, the expected value in the reference state is zero, i.e. E θ 0 [ε] = 0. To characterize the residual vector in the damaged state, the following damage hypotheses about the system states are made [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF] H 0 : θ = θ 0 (reference state)

H 1 : θ = θ 0 + δ/ √ N (damaged state), (4)
where

δ = √ N (θ -θ 0 ) (5)
is the asymptotic change vector. Then, the distribution of the Gaussian residual vector in the reference and the damaged state for N → ∞ results in [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF]

ζ = √ N ε -→ N (0, Σ) (reference state) N (J δ, Σ) (damaged state), (6)
meaning damage manifests itself in a deviation in the mean of the Gaussian residual vector. Based on these properties, statistical hypothesis tests can be applied to decide if the residual corresponds to the reference or to the damaged state. The corresponding test statistic [START_REF] Döhler | Fault detection, isolation and quantification from Gaussian residuals with application to structural damage diagnosis[END_REF]

t = ζ T Σ -1 J J T Σ -1 J -1 J T Σ -1 ζ ∼ χ 2 (ν, 0) (reference state) χ 2 (ν, λ) (damaged state) (7)
is approximately χ 2 -distributed, where the distribution is characterized through the number of degrees of freedom ν = rank(J T Σ -1 J ) [START_REF] Kammer | Sensor placement for on-orbit modal identification and correlation of large space structures[END_REF] and the non-centrality

λ = δ T J T Σ -1 J δ. (9)
The non-centrality from Eq. ( 9) describes the mean test response, i.e. the shift in the mean value of the test distribution due to the change δ, see Fig. 1. The test statistic from Eq. ( 7) is a scalar value, and an alarm is issued if the test statistic exceeds a safety threshold value t > t crit . This threshold is typically calculated in the training state based on the acceptable probability of false alarms (PFA) [START_REF] Mendler | A reliability-based approach to determine the minimum detectable damage for statistical damage detection[END_REF].

Sensor Placement Strategy

The primary goal of the sensor placement is to optimize the detectability of local damages. To achieve this goal, it has to be defined what damage is, and its effect on the statistical detection test has to be quantified. A common approach is to discretize the structure, i.e. to divide it into finite elements (FEs) and to form structural components by assigning the same structural design parameter θ h to all FEs within one component. Then, the physical system can be characterized through a parameter vector θ = [θ 1 , . . . , θ H ] T ∈ R H , which typically contains material constants or cross-sectional values, and damage can be parametrized as a deviation in the parameter vector from the reference values θ 0 through the vector θθ 0 .

The effect of damage on the damage-sensitive residual vector can be captured with a sensitivity matrix, where each column describes the sensitivity (the first-order derivative) of the damagesensitive residual vector to one design parameter θ h with

J = ∂E θ [ε] ∂θ θ=θ 0 =    ∂ ∂θ 1 E θ [ε 1 ] ∂ ∂θ 2 E θ [ε 1 ] ∂ ∂θ H E θ [ε 1 ] . . . . . . • • • . . . ∂ ∂θ 1 E θ [ε G ] ∂ ∂θ 2 E θ [ε G ] ∂ ∂θ H E θ [ε G ]    θ=θ 0 . ( 10
)
In case of the subspace-based residual from Appendix A, the residual is first derived with respect to the modal parameters, and secondly, modal parameters are derived with respect to structural parameters θ in a FE model, so the derivatives can be combined using the chain rule [START_REF] Allahdadian | Towards robust statistical damage localization via model-based sensitivity clustering[END_REF]. Whether or not damage is detectable not only depends on the sensitivity J , but also the covariance matrix Σ, i.e. the uncertainty related to the estimation of the damage-sensitive residual. The Fisher information matrix

F = J T Σ -1 J =      F 11 F 21 . . . F H1 F 22 . . . F H2 . . . . . . (sym.) F HH      (11)
considers both quantities, and is a scaling factor for the mean test response from Eq. ( 9). The main diagonal values F hh quantify the information content of the residual on each individual design parameter, and thus, quantify the detectability of changes therein [START_REF] Mendler | A reliability-based approach to determine the minimum detectable damage for statistical damage detection[END_REF]. The off-diagonal values, on the other hand, describe how similar the effect of changes in different parameters is on the residual. Since the Fisher information changes with changing sensor layout, it can be used as a performance criterion for sensor placement optimization [START_REF] Basseville | Optimal sensor location for detecting changes in dynamical behavior[END_REF]. However, to ease the physical interpretation, the next section translates the Fisher information into a more intuitive measure.

Optimization Criterion

This section explains how the Fisher information can be translated into a meaningful performance criterion for damage detectability. For the following considerations, the effect of changes in a single parameter θ h is analyzed while assuming no change in the other parameters, simplyfying the mean test response from Eq. ( 9) to

λ = δ 2 h F hh . (12)
Formulating the asymptotic change vector from Eq. ( 5) for a single parameter change, with

δ h = √ N (θ h -θ 0 h
), and substituting it into Eq. ( 12), while replacing the sample size with the measurement duration and the sampling frequency N = T f s , yields an enhanced formula for the mean test response

λ = T • (θ h -θ 0 h ) 2 f s F hh . ( 13
)
This formula is more intuitive than Eq. ( 9) because it relates the mean test response to physical quantities, such as the measurement duration during testing T and the absolute parameter change θ h -θ 0 h squared. In particular, it reveals an important property of the statistical damage detection test: the mean test response increases with increased measurement duration when the parameter change is constant. Consequently, smaller damage can be detected with longer measurements.

Eq. ( 13) gives the user the opportunity to define the requirements toward the diagnosis result as follows: The mean test response λ is a measure for the detectability of damage. On one hand, it should exceed a minimum value λ > λ min , so the distributions are clearly separated, a damaged structure is correctly classified as damaged, and the probability of detection is close to 100%, see Fig. 1. On the other hand, there is no benefit in increasing the mean test response beyond the minimum value, as no further information is gained for damage detection. Hence, a target value λ = λ min will be set for the optimization. Furthermore, the parameter change θ h -θ 0 h that should be detectable can be defined by the user. When expressed in percent of a reference parameter ∆ h = (θ h -θ 0 h )/θ 0 h , it becomes an intuitive input parameter to prioritize different structural components during damage detection. With the mean test response and the parameter change in Eq. ( 13) being fixed through physical considerations, and F hh being determined for a given sensor placement, the required measurement duration follows as

T h = 1 (∆ h θ 0 h ) 2 λ min f s F hh . ( 14
)
The duration T h is the required measurement duration to detect a specified damage extent ∆ h with the desired test response λ min . An efficient sensor configuration leads to a high Fisher information and a low measurement duration, as less data is necessary to achieve the desired probability of detection. Consequently, it is a comparable measure for damage detectability in the different structural components h, and an appropriate criterion for performance-based sensor placement. More information on how to set the input parameters is given in the following list.

• Minimum non-centrality λ min : Fixing the mean test response to a minimum value λ = λ min is equivalent to requiring a minimum reliability regarding the SHM system. It can be fixed based on an acceptable probability of false alarms (PFA), a code-based reliability index β for assessing existing structures, and the problem-specific number of degrees of freedom ν = rank(F) of the χ 2 -distribution of the test statistic. The relationship between the minimum non-centrality and the number of degrees of freedom is explained in detail in the referenced literature [START_REF] Mendler | A reliability-based approach to determine the minimum detectable damage for statistical damage detection[END_REF] with the final results being summarized in Fig. 2 for varying reliability indices. Moreover, different reliability levels λ min,h can be defined for individual structural components, depending on the damage consequences on local or global level or the accessibility during on-site inspections. Higher reliability could, for example, be required for structural components that are inaccessible, or for components that could cause a total system collapse.

• Minimum detectable damage ∆ h : The relative parameter change is a tuning parameter for the damage diagnosis. It enables the user to select relevant design parameters θ h in the FE model and to define the minimum change in these parameters that has to be sensed by the damage detection method at a pre-defined reliability. This feature is pivotal for global damage detection based on a small number of modes and a sparse sensor layout: irrelevant parameter changes in the global structure can be neglected for the optimization, and the sensor placement can be tuned to become more sensitive to damage hotspots or changes in local components, which are critical for the safety or serviceability of the structure. The framework inherently considers the magnitude of the monitoring parameter and only a relative change is to be specified. Among all structural components h, with h = 1, . . . , H, the decisive measurement duration (i.e. the measurement duration for the entire structure) is chosen as the maximum measurement duration

T = max{T 1 , • • • , T H }, (15)
as only this duration guarantees the detectability of the specified damages in all considered components. It depends on the sensor configuration c, as well as the number of sensors r. Algorithm 1 summarizes the evaluation of the objective function T = T (c, r). It should be noted that simulated vibration data from the undamaged structure y ∈ R r×N 0 of sample size N 0 is used.

Algorithm 1 Objective function T Input: FE model with monitoring parameter vector θ 0 , sampling frequency f s , minimum detectable damages ∆ h , PFA, reliability index β 1: If vibration data is not available, generate the measurement matrix y ∈ R r×N 0 2: Estimate the covariance, sensitivity, and Fisher information matrix from Eq. (11) using y; 3: Evaluate the number of degrees of freedom from Eq. ( 8) 4: Set the minimum non-centrality λ min based on β and the PFA [START_REF] Mendler | A reliability-based approach to determine the minimum detectable damage for statistical damage detection[END_REF] 5: for h = 1 to H do

6:

Calculate T h using Eq. ( 14) 7: end for 8: Calculate T using Eq. [START_REF] Kumar | Optimal location of measurements for distributed parameter estimation[END_REF] This section showed that a more efficient sensor layout leads to a higher Fisher information F hh (i.e. higher sensitivity and lower uncertainty) and a shorter measurement duration. Consequently, the optimal sensor layout for a fixed number of sensors r can be found as the one with the minimal measurement duration

T opt (r) = min c {T (c, r)}. (16)
Ultimately, an upper bound for the measurement duration is introduced (i.e. the maximum allowable measurement duration) and the optimization problem is restrained to values below one.

min c f (c, r) f (c, r) = T (c, r) T thres s.t. f (c, r) ≤ 1 (17)
This is an efficient way to remove unsuitable sensor layouts, in which some of the considered damage scenarios are undetectable. Furthermore, it allows one to weigh the objective function if it is to be combined with others in a multi-objective optimization approach. The subsequent sections will show how to find an optimal sensor placement that minimizes T , as well as an appropriate number of sensors.

Optimization Method

Exhaustive Search

A straightforward way to find the minimum of the optimization function in Eq. ( 17) is the evaluation of the performance criterion for all possible sensor configurations c in an exhaustive search (ES), and selecting the optimal configuration as the one with the shortest duration. The exhaustive approach is only applicable if the number of possible sensor combinations N c is reasonably small, because the calculation of the objective function can be computationally expensive. For r sensors and r tot candidate locations, the number of combinations N c is

N c (r, r tot ) = r tot r = r tot ! (r tot -r)!r! . ( 18
)
After simulating a single measurement record with data from all candidate locations, the corresponding channels for each sensor configuration are extracted from the measurement matrix, and the objective function is calculated and stored together with the sensor configuration number c (see Algorithm 1). Ultimately, the optimal sensor configuration number c opt is determined according to Eq. ( 17) and stored together with the corresponding measurement duration T opt .

Genetic Algorithm

This section describes how to efficiently find a close-to-optimal solution by modifying a genetic algorithm (GA) [START_REF] Deb | Multi-objective optimization using evolutionary algorithms[END_REF][START_REF] Messac | Optimization in practice with MATLAB: For engineering students and professionals[END_REF]. GAs are inspired by the natural evolution process based on Darwin's survival of the fittest principle. In this context, the objective function is referred to as the fitness function, a single sensor layout is an individual that consists of multiple genes (i.e. channel numbers), and a population is an array of individuals. For sensor placement, the optimization problem is encoded through integer numbers, meaning one gene is an integer channel number and each sensor layout is a vector that includes r numbers. Evolution takes place within a population, i.e. a large set of individuals, where the parent generation contributes its genes to its children through cross-over, mutation, and elitism [START_REF] Deb | Multi-objective optimization using evolutionary algorithms[END_REF]. For the application in this paper, the standard GA is employed that is pre-implemented in the Global Optimization Toolbox in MATLAB but the cross-over and mutation functions are modified. They include a restraint that does not allow for multiple sensors to be placed at the same locations. The mutation rate linearly decreases, meaning initially the number of perturbed channel numbers is high, but with increasing number of iterations (i.e. generations), the intensity of the perturbation gradually reduces.

The optimization is based on a single vibration record at all candidate sensor locations and launches with a random initial population. For each generation, the corresponding channels are extracted from the measurement matrix, the fitness of each individual i = 1, 2, . . . , N i is calculated (see Algorithm 1) and scaled using Eq. [START_REF] Reynier | Sensors location for updating problems[END_REF]. Based on the elite count, the fittest individuals are passed on to the next generation without modification, and a cross-over fraction determines the number of cross-over and mutation children. The selection of parents for cross-over is made in a tournament selection [START_REF] Yang | Structural optimization by genetic algorithms with tournament selection[END_REF], meaning a user-defined number of individuals (e.g. two individuals) are randomly selected in the current population, and only the fitter one is allowed to reproduce. The algorithm has the population evolve from generation to generation with g = 1, 2, . . . , N g , and terminates if one of the stopping criteria is met (i.e., a maximum number of generations is exceeded, a maximum computation time has passed, or the optimal solution has not improved by more than a specified tolerance value over a certain period of time or number of generations). Upon termination, the algorithm outputs the ranked elite members of the current population as a set of optimal solutions.

Appropriate Number of Sensors

The proposed sensor placement strategy is suitable to find an optimal sensor layout as well as an appropriate number of sensors. This is because an increasing number of sensors increases the Fisher information (increased sensitivity or lower uncertainty), and shortens the measurement duration. In any case, the optimization is launched with the minimum possible number of sensors (for example r = 1). Then, the number of sensors is gradually increased (from one to two, three, . . . ) while plotting the optimal measurement duration T opt over the number of sensors r in a performance curve and storing the corresponding sensor configuration. One way to decide on an appropriate sensor number is to terminate the optimization once the optimal measurement duration T opt drops below the maximum allowable measurement duration T thres , with

T opt < T thres . ( 19
)
Another approach is to plot the optimal measurement duration over a wide range of sensor numbers and to evaluate the elbow point (i.e. the point with the highest curvature of the performance curve) beyond which the measurement duration does not significantly improve. Alternatively, a cost function could be defined that relates the information gain with instrumentation costs. It is important to note that the optimal sensor layout can be entirely different for a different number of channels, meaning it cannot be found through sequential sensor placement, where the optimal sensor configuration is achieved by systematically adding or removing a subset of channels. Ultimately, it should be considered that an increased number of sensors leads to a longer computation time for each iteration and possibly an increased number of sensor combinations.

Application to a Steel Beam

In the following, the sensor placement strategy is applied to a pin-supported beam from Fig. 3, where the selected damage-sensitive residual is the stochastic subspace-based residual vector, see Appendix A. All considerations in this section are made based on a numerical model and simulated vibration data, while validation studies based on experimental data are reported in Section 5.

Structure and Instrumentation

The structure under consideration is a pin-supported beam with a length of 4.11 m, see Fig. 3. It exhibits a total mass of 56.8 kg, a modulus of elasticity of E = 210,000 MPa, and a hollow structural steel cross-section, HSS 152x51x4.78 mm. The supports at both ends are metal bolts in the centreline of the profile, which only allow for a rotation about the weak axis. An input source with a total mass of 3.6 kg is placed 68.5 cm away from the support, injecting white noise vertically into the beam, and eight accelerometers, with a weight of 1.28 kg each, are placed every 0.46 cm. Moreover, a 3-D finite element model with 19 joints and 18 beam elements is created in ANSYS, and possible sensor locations are defined at the positions P1 -P8 from Fig. 3. The shaker and the sensors are modelled as point masses with no rotational mass. Since 18 finite elements are used, the parametrization is limited to the 18 E-Moduli, mass densities, and cross sectional values. In this study, every two consecutive elements are combined to one component by assigning the same

with m h = 6.31 kg. For damage diagnosis, the first four modes of vibration in the vertical direction are considered, see Fig. 4. The optimization algorithm requires a vibration record with data from all possible sensor locations, so a single transient analysis is run while storing the accelerations at P1 -P8. For excitation, a Gaussian white noise input signal is applied at a single joint (the shaker location) and in the vertical direction. The output signal is generated using a modal reduction approach and sampled at a frequency of f s = 330 Hz. The measurement duration in the reference state is fixed to 110 min and all other signal processing parameters are summarized in Table B.1.

Minimum Measurement Duration

In this section, the measurement duration is calculated for a fixed sensor layout, e.g. for a single vertical channel at sensor position P1. According to Eq. ( 14), the minimum measurement duration depends on the detectable damage, which is set to ∆ h = 5% (0.31 kg) for each beam component, and the minimum non-centrality (the minimum mean test response), which is fixed to λ min = 30. Next, the Fisher information from Eq. ( 11) is calculated based on vibration data from the undamaged structure and on the FE parameter sensitivities. Ultimately, the main diagonals of the Fisher information are extracted and the minimum measurement duration is calculated for each beam component, using Eq. ( 14), with the results being visualized in Fig. 5 (left side). The measurement duration in beam segment 1 and 9 are not identical because of the asymmetrical mass distribution (the shaker is placed on beam segment 3/9), which causes the mode shapes and the sensitivity matrix to be asymmetrical. The placement of the shaker also leads to more pronounced vibration amplitudes in one half of the beam and modified covariance properties for the respective beam elements. The decisive measurement duration is the maximum measurement duration (T = T 5 = 442 s), which is used as a criterion for sensor placement. For later considerations, the measurement duration is evaluated for two additional sensor layouts-one with two sensors at positions P1 and P8 (Fig. 5, centre), and another with three sensors at P1, P4, and P8 (Fig. 5,

Optimal Sensor Layout a) Exhaustive Search

This section demonstrates how the optimal sensor configuration can be found for a fixed number of sensors, starting with the optimal location for a single sensor, r = 1. The approach to evaluate the performance criterion for a fixed sensor layout is explained in Section 3.1. The sensor is moved from c = P1 to c = P2, then to c = P3, and all other possible sensor locations from Fig. 3 while evaluating the measurement duration T h (c, 1) for all structural components h. All input parameters remain unchanged except the Fisher information changes, leading to a different maximum measurement duration for each sensor layout. For comparison, six sensor layouts are plotted in Fig. 6 together with the corresponding measurement duration (bright grey bars). The maximum measurement duration T (c, 1) = max h {T h (c, 1)} for each layout (grey bar) is copied into a performance chart and ranked, see Fig. 7. The optimum sensor layout is Configuration 8, as it exhibits the shortest measurement duration T opt (dark grey bar). The sensor layouts with one sensor at position P3 or one sensor at P6 are not considered for damage diagnosis, because Mode 3 is unobservable at these locations, see Fig. 4. As a consequence, an unreasonably high measurement duration is obtained.

b) Genetic Algorithm

In the next study, tri-axial sensors are used at the possible sensor locations P1 -P8 instead of uni-axial sensors, so r tot = 24 measurement channels are now available, see Fig. 3. Moreover, the number of channels is increased from one to r = 3 resulting in 2,024 different configurations according to Eq. ( 18). This is an appropriate case study to demonstrate the efficiency of the GA, as the number of sensor combinations is small enough for the ES method and large enough for the GA to become efficient.

First, the ES method is applied, and the layouts are ranked with respect to the measurement duration. The optimal performance is achieved if the vertical channels at the sensor locations P1, P4, and P8 are chosen, with the sensor layout plotted in Fig. 5 (right side). Secondly, the GA is launched while assigning the final rank (from the ES) to the elite members of each generation. For

Appropriate Number of Sensors

The proposed sensor placement scheme can also be used to determine an acceptable number of sensors. In the preceding section, the optimization scheme was explained to find the optimal sensor configuration for a fixed number of sensors, where the optimal solution contains two pieces of information, i.e. the optimal sensor position and the corresponding measurement duration. The same procedure is now repeated for a varying number of sensors r ∈ [1, . . . , 4]. For one to three sensors, the optimal sensor layout is plotted in Fig. 5. Plotting the optimal measurement duration T opt over the corresponding number of sensors is a particularly intuitive way to understand the effect that an added sensor has on the performance of the damage detection method, see Fig. 9. With an increasing number of sensors, the slope gradually decreases, and at some point, adding additional sensors hardly improves the damage detectability. The elbow point, where the curve begins to flatten, can be an indicator for an appropriate number of sensors. Alternatively, the first best sensor configuration could be chosen as one that exhibits a measurement duration below the maximum allowable measurement duration T thres , which was set to 240 s. In this case, two sensors seems to be an appropriate number of sensors.

Experimental Validation

Typically, the sensor placement optimization is performed before the structure is instrumented. Therefore, the optimization is performed based on simulated vibration data, as it is done in the previous section. What remains to be shown is that the optimized sensor layout (based on simulated data) actually leads to higher damage detectability using real data.

The experimental validation is performed on a laboratory steel beam that was already described in Section 4. In Section 4, the sensor layout was optimized for an optimal detectability across all beam components with T (c, r) = max{T 1 , • • • , T H }. However, this approach is not suitable for validation studies, because the masses are now only applied to component 5 or component 8, but the respective beam component might not lead to the decisive measurement duration. To remedy this, the performance criterion from Eq. ( 15) is modified to T = T 5 and T = T 8 . In other words, the optimization procedure is rerun first to optimize the detectability in the beam components 5 and 8, respectively, and the damage detection results for the ranked sensor layouts are compared.

For validation, the damage detection test from Eq. ( 7) is applied to measured vibration data before and after an extra mass of 0.31 kg is applied, which corresponds to a relative change of mass by 5%. The test is applied 100 times so the empirical distribution can be plotted in histograms and the mean test response can be determined, as visualized in Fig. 11. Whereas the minimum measurement duration (from simulated data) is used as a performance criterion for optimal sensor placement, the mean test response (based on real data) is used as a validation criterion for the actual performance.

Optimized Sensor Layout

The primary purpose of the optimization scheme is to maximize damage detectability. Maximum detectability should lead to a minimum measurement duration during the sensor placement optimization, as well as a maximum mean test response during the laboratory test. If two sensors are to be placed at positions P1 -P8, N c = 28 sensor combinations are possible. Mass on component 5: First, the damage scenario with an extra mass at midspan is considered. Fig. 12 visualizes the sensor layout for the best sensor layout (left), the 15th best configuration (centre), and the worst sensor layout (right), together with the empirical mean test response that was measured during the laboratory experiments. The respective mean test responses decrease from 45, to 31, and finally to 2. So, the worst configuration with sensors at P3 and P6 fails to detect the applied mass simply due to an inappropriate sensor placement. The result can be understood by looking at the mode shapes from Fig. 4. The sensor locations P3 and P6 are vibrational nodes of Mode 3, meaning no information from this mode of vibration is available for damage detection, resulting in a low Fisher information. In summary, a shorter measurement duration (based on simulated data) leads to a more pronounced mean test response for real data. This validates that the performance criterion is a measure for the detectability of damage and that the optimization can be performed based on simulated data.

Mass on component 8: Secondly, the scenario with a 5% extra mass near the supports is examined. Again, the best, the 15th best, and the worst sensor layouts are shown with corresponding mean test responses of 45, 31, and 20, so the sensor placement strategy successfully finds the optimal configuration. The reason why the sensor Configuration 15 in Fig. 13 is different than in Fig. 12 is because the measurement duration in individual components is used as a performance criteria instead of the maximum measurement duration over all components for validation (T = T 8 vs. T = T 5 ). It appears that the optimal sensor layout is the same for both damage scenarios (with sensors at P1 and P8), but an extra mass on Component 8 near the supports leads to a more pronounced response of the statistical test than with an extra mass on component 5 (λ = 45 vs. 27). This is consistent with the findings from Fig. 5 (middle), as the damage detectability (expressed through the required measurement duration) for component 8 is higher than for component 5 if the sensors are placed at positions P1 and P8.

Increased Number of Sensors

The results based on numerical simulations confirm the intuitive assumption that an increased number of sensors leads to higher damage detectability, so the measurement duration can also be used to decide on an appropriate number of sensors. To demonstrate this, the best configurations with one sensor, two sensors, and four sensors are illustrated in Fig. 14 together with the empirical distributions during the laboratory experiment. This figure clarifies that with an increasing number of sensors, the mean test response increases from over 30 to 41 and 46, which verifies that an increased number of sensors leads to higher damage detectability (Section 4.4 and Fig. 9).

Conclusion

This paper proposes a performance criterion for sensor placement optimization and explains how it can be combined with an adapted genetic algorithm (GA) to efficiently find a close-to-optimal solution for large structures. The primary purpose is to optimize the probability of detecting damages, where damage is defined as a deterioration of structural design parameters, such as material properties or cross-sectional values. Since damage detectability in all structural components with a limited number of sensors is unrealistic, the approach requires a prioritized set of damages scenarios whose detectability is optimized. Then, the approach is suitable to find an optimal sensor layout for a fixed number of sensors, as well as an appropriate number of sensors.

One aspect of the approach is that, although the damage diagnosis is based on global vibration measurements, it can be focused on the monitoring of local structural components that are critical for the safety or the serviceability of the structure. The proposed method requires the user to pre-select the structural parameters to be monitored and to define the detectable damage. This information tunes the sensor layout to become more sensitive to anticipated damage scenarios and yields a performance-based sensor layout. Another aspect is that the method considers the excitation and noise characteristics of the measured signals, captured through the covariance matrix of the damage-sensitive feature. Since real data is typically not available at all candidate positions, the covariance can be estimated based on simulated vibration data generated under realistic noise and excitation characteristics. Based on a laboratory steel beam, it is demonstrated that the findings based on simulated vibration records can be transferred to real structures, despite the noisy measurement environment and possibly different excitation levels. An interesting future study would be to take the leap from laboratory experiments to real-life structures; this will be the subject of future publications.

the matrix can be decomposed into its singular values and vectors The non-zero singular values in D 1 are ranked in descending order and cut off after the system order n = 2m, where m is the number of modes that are present in the signal. Correspondingly, the left singular vectors are split into the column space S 1 and the left null space S 0 . In the following, it is assumed that the matrix S 0 is obtained in the reference state. Due to the orthogonality of the two vector spaces, it holds that S T 0 H = 0 if and only if H also corresponds to the reference state, since damage alters the dynamic properties in the block Hankel matrix. Hence, a damage-sensitive residual can be formed by pre-multiplying the left null space from reference data to the Hankel matrix estimate from incoming data ε = vec S T 0 Ĥ .

H =         R 1 R 2 . . . R q R 2 R 3 . . .
(A.

2)

The hat symbol (•) indicates that the block Hankel matrix is estimated based on a finite sample size N , using the sample covariance matrix Ri = 1 N -i N -i k=1 y k+i y T k . The vectorization operator vec(•) stacks the columns of the residual matrix to obtain a residual vector that fulfills the CLT from Eq. [START_REF] Chung | On-orbit sensor placement and system identification of space station with limited instrumentations[END_REF]

Figure 1 :

 1 Figure 1: Reliability considerations based on the χ 2 -distribution of the test statistic.

Figure 2 :

 2 Figure 2: Minimum non-centrality λmin based on the reliability index for assessing structural systems [39].

Figure 3 :

 3 Figure 3: Finite element model of the HSS beam including the possible sensor locations P1-P8

Figure 4 :

 4 Figure 4: Numerical mode shapes from ANSYS

Figure 5 :

 5 Figure 5: Measurement duration to detect a 5% damage in each mass parameter for a fixed sensor layout

Figure 8 :

 8 Figure 8: Convergence of the GA towards the optimal sensor layout, determined based on an ES

Figure 9 :

 9 Figure 9: Increased test response due to an increasing number of sensors

1 .

 1 All signal processing parameters are identical and listed in Table B.1. The instrumentation consists of an electro-dynamic shaker (Mini Smartshaker, Model K2004E01) and eight tri-axial seismic sensors (Tromino Tromographs) at the sensor locations P1 -P8, see Fig. 3, synchronized through radio communication. In this study, only vertical vibration measurements are considered. Since damage was defined as a change in mass, non-destructive tests using extra masses can be used to simulate different damage scenarios. Two different damage scenarios are considered: and extra mass on beam component 5, as well as an extra mass on beam component 8, see Fig.3.

Figure 10 :

 10 Figure 10: Photos of the laboratory experiment

2 Figure 12 : 20 Figure 13 :

 2122013 Figure 11: Validation scheme

46 Figure 14 :

 4614 Figure 14: Increased test response due to an increasing number of sensors (only the optimal configurations are shown)

  .

	Appendix B. Signal Processing Parameters
	Data		Segmenting	Processsing
	Measured quantity	acceleration Training segments 100 No. of sensors	varying
	Sampling frequency	330 Hz	Testing segments	100 Time lags	p + 1/q=20
	Reference data length	110 min	Samples/segment 9,900 System order	n=18
	Training/testing data	110 min	Duration/segment 30 s Block lengths	N b =20

Table B.1: Input parameter sheet

Acknowledgments

This work was supported by a strategic research grant from the Natural Sciences and Engineering Research Council of Canada, the German Academic Exchange Service, and a Mitacs Globalink Research Award.

example, if an elite member contains the optimum channel configuration from the ES, its rank is one, where the worst solution has a rank of 2,024. The population size is set to 50 and the number of generations to 20. The elite count is set to 10%, meaning the best 0.1 • 50 = 5 sensor configurations survive to the next generation without modification. The cross-over fraction is set to 25%, so a quarter of all children, other than elite children, are made up of cross-over children. The results are summarized in a convergence chart, see Fig. 8. The GA effectively eliminates all layouts that exhibit sensors at the degrees of freedom in transverse and longitudinal direction, as expected, since only vertical modes of vibration are considered. It finds the global minimum within five generations, so 5 • 50 = 250 calculations of the objective function. That means the computation is sped up by a factor of 2, 024/250 ≈ 8 in comparison to the ES. A close-to-optimal solution with a rank of 10 is already found within four generations and 200 calculation steps. In conclusion, the GA efficiently finds an appropriate sensor layout, and outperforms the ES even in this relatively small combinatorial problem.

Appendix A. Gaussian Residual Vectors

In this section, the stochastic subspace-based residual is recalled, as an exemplary residual vector that fulfills the CLT from Eq. ( 6). The residual vector is formed in the subspace of covariance functions (in the time domain), so it can be evaluated in quasi real-time. The output covariances R i = E[y k+i y T k ] evaluate the similarities of output signals y (displacements, velocitities, or accelerations) at different sensor locations and time lags i [START_REF] Brincker | Introduction to operational modal analysis[END_REF]. When arranged in a block Hankel matrix H, with a size that is determined by the time lag parameters p and q, the dynamic properties of a vibrating structure can be retrieved from its column space. Using singular value decomposition,


Alice Fontbonne
email: alice.fontbonne@protonmail.com

Hervé Sauer

François Goudail

Comparison of three methods for end-to-end optimization of hybrid optical/digital imaging systems with professional optical design software

Keywords: Complex lens, end-to-end optimization, hybrid system, optical system design

We compare different methods to optimize end-to-end a hybrid optical/digital system for best and most uniform performance over the field-of-view with the Synopsys ® CodeV ® lens design software. We have extended the native optimization capability of this software by implementing different methods that leverage the deconvolution during the optimization process of the hybrid optical system as a whole, including simultaneously their optical and digital image processing parts. We show that the joint optimization of the lens and the processing through a true restored-image quality criterion significantly enhances the final post-processed image quality, and allows to fine-tune the residual balancing between on-axis and peripheral fields with simple weighting coefficients.

INTRODUCTION

Many imaging systems consist today of a complex optical system followed by a digital processing algorithm. It seems therefore interesting to optimize the whole imaging chain -including postprocessing -instead of separating the optimization of the optical part (i.e, the lens) from that of the processing part, as is still regularly done. This approach is called "optical/digital codesign" or "end-to-end optimization". Two different approaches are currently investigated to implement this approach. The first one consists in developing a ray tracing software that computes both the image provided by the lens and its derivatives with respect to the lens parameters. This information is then sent to a neural network (NN) that processes the image. Both the lens and the NN parameters are optimized with a stochastic gradient algorithm [START_REF] Halé | End-to-end sensor and neural network design using differential ray tracing[END_REF][START_REF] Sun | End-to-end complex lens design with differentiable ray tracing[END_REF]. The objective of this method is to be totally automatic and to dispense with the experience of the optical designer. We chose an alternative approach, which consists in taking advantage of commercial lens design software (like Zemax ® OpticStudio ® , Synopsys ® CodeV ® , etc.), whose ray tracing and optimization algorithm qualities are well known, and which constitute a familiar environment for the lens designer. These software packages do not natively take into account the influence of post-processing algorithms on lens design, but this capability can be added to them in order to implement end-to-end optimization [START_REF] Stork | Theoretical foundations for joint digital-optical analysis of electro-optical imaging systems[END_REF][START_REF] Vettenburg | Holistic optical-digital hybrid-imaging design:wide-field reflective imaging[END_REF][START_REF] Wang | An integral design strategy combining optical system and image processing to obtain high resolution images[END_REF][START_REF] Burcklen | Joint digitaloptical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software[END_REF].

In this paper, we investigate and compare three methods for optimizing a hybrid system endto-end with CodeV: a conventional method (as a reference) that does not take digital processing into account at the optimization stage, a "mixed method" that implicitly takes it into account, and real "end-to-end" method that explicitly minimizes the mean square error (MSE) between the processed image and a "ground truth". Considering, as a use-case example, optimization of a Cooke triplet to obtain the best and most uniform image quality in the field-of-view (FoV), we show that the MSE-based method is the most efficient one, and that the MSE criterion may be tuned with simple weighting coefficients to change at will the balancing between the residual image quality discrepancy between on-axis and peripheral fields. In Sec. 2, we describe the details of these optimization methods. Then, in Sec. 3, we apply them to the optimization of a Cooke triplet with wide FoV and compare the obtained image performance on the basis of quantitative quality metrics and qualitative analysis of simulated images.

OPTIMIZATION METHODS

Let us consider the set Ψ = {ψ 1 , ψ 2 , ..., ψ K } of K field positions spanning the whole FoV. We denote hψ k (ν) the optical transfer function (OTF) of the imaging system at the position ψ k , ν the spatial frequency, S oo (ν) the statistical power spectral density (PSD) of the scene (with S oo (ν)dν = 1) and S nn (ν) the PSD of the noise (set so that the signal-to-noise ratio on the raw acquired image is 34 dB, with SN R = 10 log 10 S oo (ν)dν/ S nn (ν)dν ). For S oo (ν), as natural scenes are well represented with a power-law PSD [START_REF] Daniel | Origins of scaling in natural images[END_REF][START_REF] Van Der Schaaf | Modelling the power spectra of natural images: Statistics and information[END_REF], we use a generic ideal image model S oo (ν) ∝ ν -2.5 . The MSE between the final (restored) image and the ground truth at a position ψ k , averaged over image and noise realizations, can be expressed in the Fourier domain as [START_REF] Diaz | Increase in depth of field taking into account deconvolution by optimization of pupil mask[END_REF]:

M SE(ψ k ) = | hψ k (ν) wΨ (ν) -1| 2 S oo (ν) + | wΨ (ν)| 2 S nn (ν) dν (1)
By averaging the MSE over all the field positions, one obtains the average MSE defined as :

M SE mean (Ψ) = 1 K K k=1 M SE(ψ k ) (2)
The average Wiener filter is the linear deconvolution filter that minimizes M SE mean [START_REF] Diaz | Increase in depth of field taking into account deconvolution by optimization of pupil mask[END_REF]:

wΨ (ν) = 1 K K k=1 hψ k (ν) ⋆ 1 K K k=1 | hψ k (ν)| 2 + Snn(ν) Soo(ν) , (3)
In this article, we will consider hybrid systems whose optical part is a Cooke Triplet lens with large FoV (the general specifications of the aimed Cooke triplets are given in Tab. 1) and the processing algorithm is the average Wiener deconvolution filter defined Eq. 3. In addition, we will consider three methods for optimizing a hybrid system:

• The "conventional" method is natively implemented in CodeV. It consists in minimizing the quadratic sum of the spot-diagram diameters over a limited set of field positions laying in one meridional plane. This classical method assumes equal weights for the squared spotdiagram diameters at all the considered field positions, which generally favors the on-axis performance. The lens optimized with this method will be called "CCT" (for Conventional Cooke Triplet).

• The "MTF equalization" method, which implicitly takes into account that the image will be processed with the average Wiener filter. Since it is obvious that for this filter to be efficient, the MTFs at all field positions must be similar to each other and without nulling in the spatial frequency range of interest, we use a "surrogate" optimization criterion that favors equality between the MTF on-axis and the MTFs at a set of field positions, and enforces a lower bound on the value of the on-axis MTF to make deconvolution efficient. This method allows to impose a certain homogeneity of performance across the FoV [START_REF] Burcklen | Joint digitaloptical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software[END_REF].

The lens optimized with this method will be called "HCT" (for Homogeneous Cooke Triplet).

• The "MSE optimization" takes the image post-processing explicitly into account by minimizing the average MSE defined in Eq. 2. Moreover, we will see in the following that introducing weights for the fields ψ k makes it possible to achieve controlled balancing between on-axis image quality and off-axis image quality.

Focal length 50 mm Aperture F/4 FoV 40 deg (±20 deg)

Table 1: Optical specifications of the Cooke triplets.

The MTFs of the CCT at different positions in the field, for tangential (T) and radial (R) orientations, are displayed on Fig. 1(a). As it is common for conventional optical systems (with usual weighting of fields), the on-axis MTF is much higher than the other MTFs. In particular, the MTFs at maximal field positions present quasi-nullings at frequency 40 lp.mm -1 . Conversely, the MTFs of the HCT are quite close to each other (Fig. 1(b)). This homogeneity of the MTFs leads to reduce the average MSE of the hybrid system when post-processing with the average Wiener filter. On the other hand, the fact that the MTFs are significantly lower than the on-axis MTF of the CCT has an adverse effect on the MSE. Moreover, we can see in Fig. 1(b) that the MTF is the lowest for the on-axis field, which is not, in general, a preferred configuration. In the next section, we investigate how direct optimization of the MSE criterion, which we have implemented in the CodeV software (through a compiled extension written in the C language), allows to further increase the system performance [START_REF] Fontbonne | Conception conjointe combinaison optique / traitement : Une nouvelle approche de la conception optique de haut niveau[END_REF][START_REF] Fontbonne | Comparison of methods for end-toend co-optimization of optical systems and image processing with commercial lens design software[END_REF]. It will also be shown that some slight modifications of the average MSE expression to introduce weightings between fields, allow to fine-tune the balancing between image quality on-axis and at peripheral fields.

End-to-end optimization with weighted average MSE criterion

In order to give the optical designer some freedom to emphasize some parts of the field-of-view when it is relevant for a given application, we assign a weight α k to each field ψ k , and define the weighted average MSE as follows:

M SE p (Ψ) = 1 A K k=1 α k M SE(ψ k ) (4)
where A = K k=1 α k . This results in the following expression for the weighted average Wiener filter which will be used to optimally minimize the weighted average MSE:

wΨ (ν) = 1 A K k=1 α k hψ k (ν) ⋆ 1 A K k=1 α k | hψ k (ν)| 2 + Snn(ν) Soo(ν)
.

(

) 5
The set Ψ = {ψ 1 , ψ 2 , . . . , ψ 13 } is made of 13 fields well distributed on the sensor, where ψ 1 corresponds to the on-axis position, and is associated to the weight α 1 . Throughout this article, the value of α 1 = p will vary, and all the other weights α k (for k ̸ = 1) stay equal to 1. To investigate the efficiency of the MSE-based optimization criterion and the effect of tuning the on-axis weight α 1 , we analyze in this section systems optimized with different values of p. Namely, we consider the criteria M SE 1 , M SE 2 , M SE 3 and M SE 4 , constructed respectively with α 1 = 2, α 1 = 3 and α 1 = 4, i.e. multiplying the weight of the on-axis position by 2, 3, and 4, respectively. More precisely, the "MSE1 triplet" is obtained by minimizing the M SE 1 criterion with Code V, and using the HCT as the starting point of the optimization. Then, we obtain the "MSE2 triplet" by using the M SE 2 criterion and taking the MSE1 triplet as the starting point, then the "MSE3 triplet" by using the M SE 3 criterion and the MSE2 triplet as the starting point, and finally the "MSE4 triplet" by using the M SE 4 criterion and the MSE3 triplet as the starting point. In the following, we will evaluate the impact of this change of on-axis weight with various quality metrics.

MTFs and effective MTFs

The MTFs of these four MSE triplets are presented on the first column of Fig. 2. We can notice a common feature to the four graphs: the MTFs are close to each other and without any nulling, even though they are somewhat low. They share this characteristic with the HCT, that served as a starting point for the optimization of the "MSE1 triplet". Moreover, for both the HCT and the MSE1 triplet, it is the on-axis MTF that is lower than the others. Nevertheless, increasing the weight of the on-axis field position is effective in gradually rising the on-axis MTF (red curve). For instance, the on-axis MTF value at 20 lp.mm -1 is 0.21 for the MSE1 triplet, 0.27 for the MSE2 triplet, 0.28 for the MSE3 triplet and 0.31 for the MSE4 triplet. Interestingly, this improvement does not noticeably decrease the MTF values at peripheral field positions.

Let us now observe on the right column of Fig. 2, the effective MTFs of the MSE triplets. The effective MTF is the product of the MTF (at a given position in the FoV) with the frequency response of the deconvolution filter -which is, in our case, the average Wiener filter (Eq. 3, without weighting). In the ideal case, and in the absence of noise, the effective MTF would be equal to 1 at all spatial frequencies. A value greater than 1 at a given spatial frequency indicates over-contrast, while a value smaller than 1 indicates under-contrast relatively to the ideally restored image, which is the ideal image of the scene without any blur due to lens aberrations and diffraction. The effective MTF obtained with the average Wiener filter is not uniformly equal to 1, since it corresponds to the best compromise (in the M SE p sense) between good reconstruction of the useful part of the signal and mitigation of noise reinforcement. We can notice on Fig. 2(b) that in terms of effective MTF, the MSE1 triplet clearly performs poorer on-axis than in the peripheral positions of the FoV. With the MSE2 triplet, the onaxis curve rises strongly (Fig. 2(d)). Such an improvement is also visible between the MSE2 and the MSE3 triplets. For example, the effective MTF rise from 0.58 at 40 lp.mm -1 for the MSE2 triplet (Fig. 2(d)) up to 0.62 for the MSE3 triplet (Fig. 2(f)). Finally, the MSE4 criterion further improves the performance on the axis, keeping the effective MTF higher than 0.7 for spatial frequencies even larger than 40 lp.mm -1 . Therefore, rising the on-axis weight α 1 yields a constant improvement of the on-axis MTF and effective MTF without any significant degradation at peripheral positions in the FoV.

Quantitative performance evaluation

The effective MTFs are not sufficient to fully describe the performance of a hybrid system because they do not show the noise reinforcement introduced by the deconvolution filter. A better way to understand the influence of the processing on the final image quality is to evaluate the local MSE at a specific position ψ k in the FoV. Let us first consider the MSE without deconvolution (which is obtained by substituting, in Eq. 1, wΨ (ν) with a function uniformly equal to 1). Figure 3(a) represents its evolution as a function of the position in the FoV for the different Cooke triplets considered in the article. We notice that the best MSE (without deconvolution) is obtained with the CCT, in particular on the axis, where its value is 8 times lower than that of the MSE1 triplet. This was expected since the optimization of this triplet does not take digital processing into account. The MSE values obtained with the MSE triplets and the HCT are thus larger and close to each other, which was also expected since they are optimized to work with a digital processing. Moreover, we can notice that the MSE variations are correlated with the variations of the MTF: if the MTF at a given position in the FoV has higher values than another (especially for low spatial frequencies), then the MSE a this position is lower. Looking at the details, we also observe that increasing the weight of the central field (i.e., passing progressively from M SE 1 to M SE 4 ) allows to slightly improve the average performance before deconvolution, and more particularly on the axis. We can notice that the improvement brought by deconvolution clearly depends on the lens on which it is used. It is quite low for CCT, substantial for HCT and very significant for the different MSE triplets. It is important to note that, as desired, the performance at the center of the FoV gets better and better when increasing the weight on it, i.e., when passing from M SE 1 to M SE 4 . On the other hand, the performance at intermediate positions in the FoV (from 5 • on) is slightly degraded, but the four MSE triplets have equal performance at the maximum field position (20 • ).

This slight degradation of the local MSE with the field angle results in a decrease in average performance (after deconvolution) of the MSE2, MSE3 and MSE4 triplets. The "Image quality" (IQ) defined as:

IQ = 10 × log 10 1 M SE mean (Ψ) , (6)
where M SE mean (Ψ) is defined in Eq. 2, can be used to quantitatively study performances over the FoV. It is expressed in dB and its value for the 4 considered MSE triplets is given in Tab. 3.2. The MSE1 triplet is, as expected, the best one (16 dB), as it has been optimized with the same criterion that we use to evaluate its performance. However, this performance is not so far from that of the MSE4 triplet (-0.3 dB), which, on the other hand, has a better performance on the axis (Fig. 2(g-h) and Fig. 3(b)). Moreover, the performance with deconvolution of all the MSE triplets is significantly better than that of the CCT and HCT, which is equal to 13.5 dB for both lenses.

Qualitative evaluation on simulated images

A complementary method to evaluate the benefit of weighting the optimization criterion to achieve better quality on-axis is to use image simulations. Figure 4 represents a subpart (close to the on-axis position) of the images obtained with the different MSE triplets. The images in the first column are obtained without deconvolution, and the second column is obtained with deconvolution (Fig. 4(i) being the ground truth, i.e. the ideal scene image). Without deconvolution, we can see the improvement of on-axis visual quality as the weight of on-axis position increases, since from Fig. 4(a) to Fig. 4(g), images are more and more contrasted. Nevertheless, to approach the sharpness of the ground truth (Fig. 4(i)), it is necessary to deconvolve (Fig. 4(second column)). After deconvolution, the visual quality is quite similar between the four MSE triplets, but we can still observe the expected improvement. In particular, the "stone lace work" in the dark background (bottom left part of the image) is slightly blurrier in Fig. 4(b), obtained with the MSE1 triplet, than in Fig. 4(f), obtained with the MSE3 triplet, and in Fig. 4(h), obtained with the MSE4 triplet. This fits well the theoretical conclusions drawn from Fig. 3(b).

Figure 5 represents the same data as in Fig. 4, for another subpart of the image located at a large field angle. We can notice once again the great improvement between the case without deconvolution (Fig. 5(first column)) and the case with deconvolution (Fig. 5(second column)). Moreover, we can observe that after deconvolution, there is very little visual difference between the four MSE triplets, which corroborates the result obtained Fig. 3 for positions in the field between 15 • and 20 • . To conclude, image simulations leads to similar conclusions than the theoretical predictions and show that optimization with weighted MSE criteria (M SE 2 , M SE 3 and M SE 4 ) can improve the on-axis performance without significantly decreasing the average performance over the whole FoV, and in particular the performance at the peripheral positions of the FoV.

CONCLUSION

We have compared three different methods to optimize hybrid systems based on a Cooke triplet and a digital deconvolution: A conventional one, as a reference, where the lens is optimized without taking into account the post-processing, which is introduced in a second step; a second one with optimization using MTF equalization criteria that implicitly take into account the image restoration by deconvolution; and a third one that truly optimizes simultaneously the lens and the deconvolution algorithm through a restored-image quality criterion (the MSE), which has been implemented in practice in the CodeV software. For the aim of getting the best and most uniform post-processed image quality all over the field-of-view of a Cooke triplet, we have shown that the second method performs better than the first, and that the third overperforms the second. Moreover, the MSE criterion can be simply tuned with weights in order to emphasize some parts of the FoV, allowing in particular to balance at will image quality in central and peripheral field positions. These results show that a true hybrid system co-optimization criterion can be efficiently implemented in a professional lens design software, while giving the optical designer the flexibility to adapt the optimization criterion to its applications and to use all the software classic tools for system analysis and enforcement of manufacturing constraints.

This work has several perspectives. For now, the MSE triplets remain quite close from the starting point of the optimization. It would be useful to improve the co-design tool using the weighted MSE criterion in order to better explore the optimization landscape and move farther away from starting point. It is also possible to envisage the implementation of other co-design criteria or other deconvolution algorithms in the lens design software CodeV.

Figure 1 :

 1 Figure 1: (a) MTF of the CCT (b) MTF of the HCT.

Figure 2 :

 2 Figure 2: (First column) MTFs (Second column) Effective MTFs, for (a-b) the MSE1 Cooke Triplet, (c-d) the MSE2 Cooke Triplet (e-f) the MSE3 Cooke Triplet and (g-h) the MSE4 Cooke Triplet.

Figure 3 :

 3 Figure 3: Local MSE as a function of the field position in the FoV (a) without deconvolution; (b) with deconvolution with the average Wiener filter (without weighting).

Figure 3 (

 3 Figure 3(b) displays the MSE obtained by the different triplets after deconvolution with the average Wiener filter. We can notice that the improvement brought by deconvolution clearly depends on the lens on which it is used. It is quite low for CCT, substantial for HCT and very significant for the different MSE triplets. It is important to note that, as desired, the performance at the center of the FoV gets better and better when increasing the weight on it, i.e., when passing from M SE 1 to M SE 4 . On the other hand, the performance at intermediate positions in the FoV (from 5 • on) is slightly degraded, but the four MSE triplets have equal performance at the maximum field position (20 • ).This slight degradation of the local MSE with the field angle results in a decrease in average performance (after deconvolution) of the MSE2, MSE3 and MSE4 triplets. The "Image quality" (IQ) defined as:

Figure 4 :

 4 Figure 4: Center part of an image simulated through a hybrid system, (First column) without deconvolution and (Second column) with deconvolution with an average Wiener filter. (a-b) MSE1 Cooke Triplet, (c-d) MSE2 Cooke Triplet (e-f) MSE3 Cooke Triplet and (g-h) MSE4 Cooke Triplet. (i) Ground truth. All simulations of deconvolution have been done with a 34 dB SNR.

Figure 5 :

 5 Figure 5: Edge part of an image simulated through a hybrid system, (First column) without deconvolution and (Second column) with deconvolution with an average Wiener filter. (a-b) MSE1 Cooke Triplet, (c-d) MSE2 Cooke Triplet (e-f) MSE3 Cooke Triplet and (g-h) MSE4 Cooke Triplet. (i) Ground truth. All simulations of deconvolution have been done with a 34 dB SNR.

Table 2 :

 2 IQ of Cooke triplets with or without deconvolution.

	Cooke Triplet Without deconvolution With deconvolution with an average Wiener filter
	MSE1	10.1 dB	16.0 dB
	MSE2	10.6 dB	15.9 dB
	MSE3	10.5 dB	15.9 dB
	MSE4	10.8 dB	15.7 dB