Postorino Hadrien
email: hadrien.postorino@ensam.eu

Rebillat Marc

Eric Monteiro

Mechbal Nazih

Towards an industrial deployment of PZT based SHM processes: A dedicated metamodel for Lamb wave propagation

Keywords: SHM, Lamb Waves Propagation, Metamodel, Inverse Problem

Numerical simulations of Structural Health Monitoring processes based on wave propagation can be very costly in terms of computation time, especially for complex aeronautic composite structures, and therefore strongly limits the deployment of industrial applications. Metamodels build a relatively simple relationship between inputs and outputs from a set of data and thus can overcome that difficulty. A metamodel based on radial basis functions interpolation is build in order to predict a Lamb Wave measurement on a damaged composite plate equipped by a network of 3 piezoelectric elements. The input parameters describe the position of the damage. This surrogate model is used to predict the measured signals for new damage configurations with a limited computational cost. Moreover, this metamodel is used in a reverse way to solve the inverse problem. A swarm particle optimisation algorithm seeks to find the position of a damage from a set of simulated signals.This approach allows us to identify correctly the damage localisation for an unknown configuration, providing therefore a new method for damage localization.

Introduction

In the past decades, composite materials have become omnipresent in the aeronautical, automotive, space and nautical industries [START_REF] Giurgiutiu | Structural health monitoring of aerospace composites[END_REF]. They outperform metallic materials in terms of mechanical and chemical resistance. They allow for the design of lightweight structures with unique mechanical and geometric characteristics. They offer the possibility of reducing the weight of aircraft and increasing their lifespan. However composite structures are also subject to damages caused by fatigue, defect during manufacturing or impact incidents. These damages can typically take the form of cracks, holes, delimitation or discontinuity. Non-destructive testing still plays a major role in the industry to ensure the integrity of a structure but needs the interruption of service of the plane for regular human visual inspections.

Thus, the creation of smart structures able to monitor and autonomously determine their structural health has become an important subject of research and innovation over the last decades [START_REF] Cawley | Structural health monitoring: Closing the gap between research and industrial deployment[END_REF]. These new control and monitoring methods are called Structural Health Monitoring (SHM). SHM has the main objective to determine the structural health in real time. It will then be possible to adapt the maintenance strategy to the actual state of the structure, allowing to optimize the costs and to guarantee more safety.

SHM systems should ideally meet the following requirements [START_REF] Ostachowicz | Optimization of sensor placement for structural health monitoring: a review[END_REF] : low cost, continuous evaluation, sensitive to small and different types of damage, insensitive to operating conditions, noise and environmental changes. In order to evaluate the maturation of SHM techniques, four levels are distinguished: (i) Damage Detection, (ii) Damage localisation, (iii) Damage quantification and (iv) Damage characterization.

The most advanced techniques in SHM reach level 2 with satisfaction. On the other hand, level 3 and 4 are still active research subjects. These techniques can be grouped into three categories : vibration monitoring, deformation monitoring and wave propagation monitoring. This paper will focus on Lamb wave propagation monitoring.

Lamb waves have the particularity to propagate in plates with low attenuation. When the wave propagates through a damage, it is modified by the damage. By studying and analyzing the wave propagation in the structure, it is then possible to determine the structural state of the structure. The generation and measurement of Lamb waves is done using a network of piezoelectric transducer distributed over the entire structure under study. These PZT transducers act alternately in actuator mode and in sensor mode. The data is then collected and processed to determine the state of the structure. For this purpose, the recorded signals of unknown state are compared with the previously recorded signals of a healthy state. Different values, called Damage Index, allow the comparison of the signals.

The SHM has been subject to a significant amount of works over the past several decades [START_REF] Cawley | Structural health monitoring: Closing the gap between research and industrial deployment[END_REF]. Applications already exist in civil engineering for the monitoring of buildings, civil structures, roads and railways. In the aeronautics industry, the industrial deployment of research advances has been slow to take place. Various explanations can be given, but one of the main issue is the difficulty to evaluate correctly the performance of SHM systems. The evaluation of SHM systems passed by the computation quantity of interest such as the Probability Of Detection or Probability Of Localisation [START_REF] Schubert Kabban | The probability of detection for structural health monitoring systems: Repeated measures data[END_REF]. This involves to collect a huge amount of data on structures with a lot of damage configurations in order to make statistical evaluation. A fully experimental approach is thus difficult and costly: Moriot [START_REF] Moriot | A modelbased approach for statistical assessment of detection and localization performance of guided wave-based imaging tech[END_REF] robustly eval-uated their SHM processes based on experimental data with a plate equipped by piezoelectric transducers.

An evaluation based on numerical data appears to be more feasible. However, to simulate SHM processes is not a trivial task because of the complexity of the aeronautic structures and the large diversity of damage. Numerical simulations can be costly in terms of computational time, but are necessary in order to evaluate the processes and determine their effectiveness.

The creation of a metamodel (or surrogate model) able to predict the effect of a damage on a structure would considerably reduce the computation time and allow to test SHM processes on a large number of damage configurations. Building a metamodel means looking for a relatively simple relationship between inputs and outputs. Many methods of metamodelisation exists for a very large variety of applications. In the context of Lamb waves SHM, Borate et al. [START_REF] Borate | Data-driven structural health monitoring approach using guided lamb wave responses[END_REF] recently used Proper Orthognal Decomposition and Artificial Neural Network in order to predict the position of a damage in a numerical plate and in a beam based on the measurement of the displacement and acceleration.

An inverse use of the metamodel is also possible in order to find back the parameters that lead to a result of a simulation. The use of metamodels to solve an inverse problem has been implemented by Zhang [START_REF] Zhang | Vibration-based delamination detection in curved composite plates[END_REF] to find the characteristics of a damage from the measurement of eigenfrequencies. This numerical and experimental study shows that the use of a metamodel combined with a genetic algorithm obtains better results than a neural network to identify the parameters of a damage. Joy [START_REF] Joy | Implementation of kriging surrogate models for delamination detection in composite structures[END_REF] has also shown that the use of a metamodel for inverse problem resolution is more efficient than a direct neural network.

In this paper, a numerical model of a simple aeronautic structure with one damage is proposed. The parameters describing the position (x dam , y dam ) of the dommage are used to generate a database with different damage configuration. Then a metamodel based on Radial Basis Function (RBF) interpolation is build in order to strongly reduce the computation time of the simulations for new damage positions. This surrogate model can be used to predict and generate a high quantity of data. It is then used to solve the inverse problem and find the position of a damage from SHM measurements. That lead to a new strategy for damage localisation.

Numerical Model

The finite element simulation is computed with the SDT software developed by SDTools.

The structure under study is a 300x300x2.4 mm plate in a composite material for the aeronautic industry. The plate has 8 plis [0 / 45 / 0 / 45 / 45 / 0 / 45 / 0]. The mechanical properties can be seen in the table 1. The plate is modeled by shell elements with an linear elastic orthotropic behavior.

The damage is modeled by a local cylindrical reduction of the all Young and Shear modulus across the plate thickness. The lower the ratio of the module it is, the more severe the damage. Here the diameter of the damage is fixed at 10 mm and the severity at 50%.

The damping introduced in the model is stiffness proportional damping with a coefficient β = 1.5e-8 for the mass term. Concerning the mechanical boundary conditions, free boundaries are chosen.

In order to monitor its structural health, the plate is equipped by 3 PZT. Each PZT can be used as actuator or sensor. The excitation signal is a burst with 5 cycles at 150 kHz for the central frequency (figure 1).

The plate is meshed by a regular grid. The excitation frequency is fixed at 150 kHz, which corresponds, according to the dispersion curves, to a wavelength λ ≈ 20mm for this material. The dispersion curve for lamb wave can be seen figure 2. The mesh size must be less than 10 % of the wavelength, that is to say less than 2 mm. A metamodel can be seen as a black box which for any input vector x associates a prediction Ŷ, such as Ŷ = M(x). The application M is computed during a training phase from a set of previously computed input and output. Let p be the number of inputs parameters and q the number of outputs. The set of n inputs and outputs can be assembled respectively in a matrix X and Y. The term X i j corresponds to the value taken by the j-th parameter during the i-th simulation. The two matrices X, Y contain the data used to train the metamodel. The data (X j , Y j ) associated to the i-th simulation are called a Snapshot or sampling point. The snapshots are computed beforehand from a numerical model, or measured from experimental data but can also be a mix of numerical and experimental data within the framework of a hybrid model [START_REF] Chinesta | Virtual, digital and hybrid twins: A new paradigm in data-based engineering and engineered data[END_REF].

Table 1: Mechanical properties of the composite material

Name E 1 E 2 ν 1 ν 2 ρ G 1 G 2 Value 64.
X =         X 1 1 • • • X n 1 . . . . . . . . . X 1 i • • • X n i . . . . . . . . . X 1 p • • • X n p         = x dam 1 • • • x dam i • • • x dam n y dam 1 • • • y dam i • • • y dam n and Y =         Y 1 1 • • • Y n 1 . . . • • • . . . Y 1 i • • • Y n i . . . • • • . . . Y 1 q • • • Y n q        
A surrogate model is build in order to predict a Lamb wave SHM control process (figure 3) based on the numerical model of the structure previously described. One damage is included on the structure. Two parameters (x dam , y dam ) describing the damage position are considered. The size and the severity of the damage are fixed. So the inputs are X = (x dam , y dam ) with x dam the horizontal position of the damage and y dam the vertical position of the damage. The Y output is the measured signal on a sensor for one actuator. Each way of actuator/sensor are considered independently : it is necessary to build one metamodel for each way of actuator/sensor.

Sampling

The construction of a metamodel requires beforehand the computation of a database X, Y resulting from numerical simulations or experimental data. The choice of input variables X for the training is not trivial because it strongly influences the precision and the speed of the metamodel. This choice of parameters X is called Design Of Experiment, or Sampling within the framework of metamodels. Many DOE exist: random sampling, Latin hypercube Sampling [START_REF] Kleijnen | Kriging metamodeling in simulation: A review 192[END_REF], full grid sampling, scattered grid sampling, etc. In this work, the sampling strategy is based on a randomly perturbed regular grid [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF] and can be seen figure 4. 196 damage configurations are computed with 3 different actuators, so a total of 588 simulations.

Signal Compression

The choice of the data representation is significant to obtain an efficient metamodel but also to facilitate its training. The data representation must be both parsimonious and, ideally, also have a physical meaning. The question is in which format express the output matrix Y of the metamodel. The Singular Value Decomposition offers the possibility to compress data sparsingly and easily.

The SVD gives a base in which pro-

∑ n r i=1 λ i ∑ n i=1 λ i ≥ 98% (1)
where n s is the number of snapshots, n r the number of modes in the reduced basis and {λ i } n 1 are the eigenvalue [START_REF] Guo | Data-driven reduced order modeling for timedependent problems 345[END_REF]. For this criteria, the reduced basis counts 16 vectors (figure 5).

Radial Basis Function Interpolation

Different metamodeling methods have been tested and evaluated. The method that got the best results is the Radial Basis Function interpolation (RBF). RBF were deployed by Hardy [START_REF] Hardy | Multiquadric equations of topography and other irregular surfaces[END_REF] for the interpolation of multi-variate scattered data. Radial basis functions have the particularity to decrease monotonically with distance from a central point. Many radial basis functions exist. In that section, it is considered that the metamodel is build only for one output (q = 1). For multiple outputs metamodels, each input is treated independently. The RBF model computes a linear combination of radial basis function ϕ centered on each X i sampling point. Each center X i is associated with a weight ω i such that :

∀x ∈ R p , Ŷ(x) = n ∑ j=1 ω j ϕ( x -X j ) (2)
The coefficients ω i can be computed by applying the interpolation condition :

∀i ∈ [1 : n], Y i = n ∑ j=1 ω j ϕ( X i -X j ) (3)
resulting a linear system of equations:

Y = ϕω ( 4
)
where Y is the vector of values of the function at the sampling points, ω is the vector of weights of the radial basis function, and the matrix ϕ is also known as the Gram matrix:

ϕ =    ϕ( X 1 -X 1 ) • • • ϕ( X n -X 1 ) . . . . . . . . . ϕ( X 1 -X n )) • • • ϕ( X n -X n )    (5)
The weights ω are then computing by inverting the matrix ϕ or using a QR decomposition for large matrices. The main parameters for RBF metamodels are the choice of the kernel function and the value of the kernel width defined by σ. In the presented results, an Inverse MultiQuadric kernel function ϕ(r, σ)

= 1 √ 1+ r σ is used
with a width of 0.06.

Metamodel Evaluation

Of the 196 snapshots from the simulations, 75% are used for training and 25% are kept for validation. The error between the predicted signal Ŷ and the validation signal Ŷ is computed with a FFT based implementation of the maximum of the correlation. This error is also a Damage Index used for damage detection.

Error = 1 -max   IFFT FFT( Ŷ)FFT(Y) E Ŷ E Y   (6)
Here the FFT is the Fast Fourrier Transform, IFFT the Inverse Fast Fourier Transform, E the mean value, Ŷ the predicted signal and Y the reference signal. The results obtained by this surrogate model are presented in the figures 6-9. Figures 6 and7 compare the predicted and simulated signals with minimum and maximum error respectively at the sensor 2 with the actuator 1. There is a good match between the signals, indicating that the metamodel correctly predicts the effect of a damage. The figure 9 shows the evolution of the error with the distance to the nearest training point.

Inverse Problem resolution

Metamodels allow to predict very quickly the results of a long simulation. Thus, a surrogate models can be used to predict signals on a large number of damage configurations. An inverse use of the metamodel is possible in order to find back the parameters that lead to a result of a simulation with the help of an optimisation algorithm. The surrogate model is used to find the (x dam , y dam ) position of a damage from a set of simulated signals. The signals come from the validation database, they have not been used to train the metamodel. A swarm particle minimization agorithm [START_REF] Kennedy | Particle swarm optimization[END_REF] seeks to minimize the global error. The global error here is the sum of the cross correlation based error (equation 6) for each actuator/sensor configuration between the signals predicted by the metamodel and the signals to be found Error i,j

where Error i,j is the error between the predicted and simulated signal at the sensor j with the actuator i. The optimisation process can be seen figure 14. This method has been tested on 30 points that don't belong to the training set. The results are shown in figures 10-13. In the figure 10, the predicted position appears in red while the actual damage position is displayed in blue. The optimisation algorithm reached the convergence for all the test points. This method gives promising results but some points remain far from the actual position of the damage (points 1, 2, 8, and 14). The reasons for these discrepancies are being investigated.

Conclusion

The previously described surrogate model allows to quickly predict the effect of a damage on the structure. It could be used to predict a large amount of data in order to evaluate statistically the SHM algorithms or to generate a database for other machine learning approach in SHM. The resolution of inverse problem based on a surrogate model gives promising results for damage localisation and must be tested with noisy and experimental data. Other parameters -such as the size, the severity, the thickness -need to be integrated into the damage model in order to gain realism.

Moreother, the damage model used in this work is probably not the most realistic. It will be necessary to study other damage models and fit numerical results with measurements carried out on a real plate. The inverse problem method should be extend to a metamodel with more parameters such as size, severity or shape of the dommage. That could lead to an efficient strategy to monitor structural health.

Acknowledgement

This work was supported buy the French Government via the program Fonds Unique Interministériel (project Monarque).

Output signal actuator 1 sensor 2 Fig. 1 :Fig. 2 :

 212 Fig. 1: Example of input and output signal

Fig. 3 :Fig. 4 :

 34 Fig. 3: Input/Output of the metamodel for actuator i to sensor j

Fig. 5 :

 5 Fig. 5: Cumulative energy of SVD

Fig. 6 :Fig. 7 :

 67 Fig. 6: Comparison of predicted/simulated signal with minimum error for actuator 1 to sensor 2

Fig. 8 :Fig. 9 :(

 89 Fig. 8: Error histogram for actuator 1 to sensor 2

Fig. 10 :Fig. 11 :Fig. 13 :

 101113 Fig. 10: Damage position estimation

Fig. 14 :

 14 Fig. 14: Block Diagram for the inverse problem resolution with particle swarm optimisation


Hanane Ariouat
email: hanane.ariouat@ut-capitole.fr

Eric Andonoff
email: eric.andonoff@ut-capitole.fr

From Declarative Knowledge to Process-based Crisis Resolution: Application to Flood Management

Crisis resolution is often based on official government plans that provide guidelines. In real time, when a crisis occurs, one or several plans have to be chosen, merged, refined to meet the specific requirements of the crisis, and then launched. Plans are often in a textual format, which makes their interpretation ambiguous and error prone. Therefore, in real time, the coordination of stakeholders becomes difficult and time consuming. Given these drawbacks, the transformation of a plan into a process provides several advantages: i) an accurate and machinereadable specification of coordination of actions to be done in the field, ii) a better common understanding between stakeholders responsible for these actions and iii) a mean to analyze, simulate and evaluate the crisis response before launching it. The problem being addressed in this paper is "how to deduce a process for driving crisis resolution from business knowledge (plans, stakeholders and their capacities) and relevant facts observed in the impacted field". This paper presents first a meta-model for capturing business knowledge and crisis situation and then a deduction approach deriving a process in a BPMN-like format. Flood of the Loire in June 2016 serves as a support for approach experiment.

Introduction

In crisis management, several participating actors (stakeholders) have to act simultaneously and urgently to reduce the crisis and its impacts on the real world [START_REF] Tahir | A collaborative information system architecture for process-based crisis management[END_REF]. To achieve this common goal efficiently, these actors must collaborate, or at least act in a coordinated way in order to make their activities as efficient as possible. In France, in civilian crisis management such as flood, fire, chemical accident or terrorist attack, this coordination definition is under the responsibility of a control center, called crisis cell, headed by a Préfet or by the Ministre de l'Intérieur, depending on the crisis scale. The cell, which is composed of the representatives of different public organizations involved in crisis resolution and which may be geographically distributed, is in charge of applying governmental plans defined by the law and providing guidelines for action in response to the crisis [START_REF] Andonoff | Interaction protocols for human-driven crisis resolution processes[END_REF]. However, these plans provide general guidelines and have to be adapted to meet crisis requirements. In addition, the coordination between actors involved in crisis resolution is not specified and must be defined by the crisis cell.

This paper deals with the design of the coordination of actions of actors involved in the field to reduce the crisis. It recommends a process-based approach to address this issue. Indeed, there is an easy mapping between governmental plans and processes [START_REF] Ariouat | Do process-based systems support emergent, collaborative and flexible processes? Comparative analysis of current systems[END_REF]: actions and actors from governmental plans correspond to process activities and roles while coordination of actions is explicitly modelled in the process using coordination patterns such as sequence, alternative, or parallelism between activities [START_REF] Dumas | Fundamentals of Business Process Management[END_REF]. Moreover, a process-based approach provides an understandable, accurate and machine-readable specification of actions to be done in the field as well as a means for analyzing, simulating and evaluating the crisis response before launching it.

More precisely, the paper contributes to make the design of this coordination easier for crisis cells. It is based on the following cycle driving crisis resolution. This cycle, visualized in Fig. 1, is made of three steps.

Fig. 1. Crisis resolution cycle

The first step is the identification of relevant facts from data recorded in the field. The second step is the deduction of the corresponding Crisis Resolution Process (CRP) taking into account business knowledge of actors involved in crisis resolution. Finally the third step is the execution of relevant actions according to the CRP, which leads to the modification of the crisis situation.

The paper focuses on the deduction step, which is a fundamental step in the crisis resolution cycle. Deduction of CRPs has already been addressed in literature (e.g., [START_REF] Bénaben | Semantic issues in model-driven management of information system interoperability[END_REF][START_REF] Macé Ramète | Towards a predictive model for decision support in road crisis management[END_REF]), but the main drawback of existing contributions is that they only support deduction of services organized in sequence or in parallel. They are unable to deduce choices and thus they are unable to model alternative plans that crisis cells have to consider. In addition these solutions fail in explaining why services are selected, both in terms of facts observed (a service is selected to deal with a risk or damage), and in terms of services (a service is selected as it is necessary to another one).

The current paper addresses these two previous issues by providing two main contributions: a metamodel including knowledge to derive and justify the CRP deduction and an algorithm able to mine a CRP including alternatives plans. Also an experimentation of our approach on a real case study, namely the flood of the Loire in June 2016 is given.

In the reminder of this paper, we first place our CRP deduction approach with regard to related work. Sec. 3 presents our meta-model defining the concepts needed to represent both facts describing a given crisis and business knowledge related to crisis resolution. Sec. 4 is dedicated to the deduction of CRP. It presents the deduction approach and introduces the main algorithms implementing it. Sec. 5 details our experimentation through a specific instantiation of our meta-model. Finally the conclusion summarizes the paper and mentions some open issues.

Related work

We have found several contributions advocating a process-based approach in the field of crisis management and each of them has implemented algorithms for process deduction (e.g., [START_REF] Bénaben | Semantic issues in model-driven management of information system interoperability[END_REF][START_REF] Macé Ramète | Towards a predictive model for decision support in road crisis management[END_REF]). In these contributions, as in ours, the inferred process is expressed as a BPMN-like diagram. However, in each of these contributions, deduction differs from ours.

First the contribution described in [START_REF] Bénaben | Semantic issues in model-driven management of information system interoperability[END_REF] recommends a composition-based deduction from input and output of services that have to be deployed in the field. Both sequencing and parallelization of services can be deduced. More precisely a service having outputs corresponding to inputs of another one is represented by two BPMN tasks, corresponding to the considered services, linked together by a sequence flow. For instance, if all the outputs of service a correspond to all the inputs of service b, then services a and b are modelled as tasks connected by a sequence flow from a to b in the BPMN diagram. In addition, when a service has outputs corresponding to inputs of several services, these last services are modelled as parallel services. For instance, if outputs of service a (e.g., e1 and e2) correspond to all the inputs of service b and if the other outputs of service a (e.g., e3) correspond to all the inputs of service c, then services b and c are modelled as parallel tasks executed after service a.

Second, the contribution described in [6] recommends a knowledge-based deduction. More precisely this works recommends a specific OWL ontology for both crisis situation and crisis response modelling. It also introduces specific SWRL rules for both selection of services to be deployed in the field and their ordering as sequencing or parallelization. This second contribution is interesting because, unlike the first one ( [START_REF] Bénaben | Semantic issues in model-driven management of information system interoperability[END_REF]), which recommends a compositionbased deduction according to inputs and outputs of services as in web service composition, it takes into account knowledge of actors involved in the field introducing the notion of objective and linking objectives with services to be deployed to reach them. However, this knowledge is expressed within rules which are difficult to define and hard-coded in the ontology. In addition, these rules only deduce sequencing and parallelization of services: they are unable to deduce alternatives to reduce risk or deal with damage. However, deducing alternatives for driving crisis resolution is really relevant for crisis cells, which need to know the possible solutions before making decisions.

We also have found several contributions addressing process deduction in the field of process mining [START_REF] Van Der Aalst | Process Mining -Data Science in Action[END_REF]. The main one is the Alpha algorithm [START_REF] Van Der Aalst | Workflow Mining: Discovering process models from event logs[END_REF], which identifies process schemas from log file repositories. These log files record executed activities (services), actors performing them and corresponding timestamp. Alpha mines these log files to identify the coordination of executed activities by actors. Moreover, some works have addressed the deduction of activities to deal with a given situation. First, the contribution described in [START_REF] Schonenberg | Supporting flexible process through recommendations based on history[END_REF] recommends a deduction service that is implemented in the ProM tool and that aims at providing the user with the next possible activity to perform in a given situation. This service exploits process log repositories to match the running process with existing cases. In [START_REF] Ribeiro | A recommender system for process discovery[END_REF], the recommended system compares different process mining algorithms to identify the one which best fits a given situation. Comparison is based on measurements such as fitness and generalization, which allow the evaluation of the performance and the quality of the compared algorithms. Finally, [START_REF] Rangiha | A framework to capture and reuse process knowledge in business process design and execution using social tagging[END_REF] describes a recommender task system that uses social tagging to collect relevant information from discussions between process actors during process execution. Analysis of these tags allows the system for deducing new tasks when the same process must be executed again.

The approach recommended in this paper is in line with the one described in [START_REF] Macé Ramète | Towards a predictive model for decision support in road crisis management[END_REF]. As in [START_REF] Macé Ramète | Towards a predictive model for decision support in road crisis management[END_REF], we advocate a knowledge-based deduction considering both facts observed in the field and business knowledge of crisis actors.

More precisely business knowledge corresponds to already existing solution modelled as services offered by crisis actors and possible corresponding plans defining coordination between them. In addition, we also model relations of these services in terms of what do they require, what do they cause, why using one or another. These relations correspond to knowledge of crisis actors indicating how services must be performed in the field. Unlike [START_REF] Macé Ramète | Towards a predictive model for decision support in road crisis management[END_REF], we do not model this knowledge as a set of hardcoded rules in an ontology but rather as data stored in a database. Thus our approach is declarative, making knowledge management easier. Moreover, we fully exploit this knowledge as we deduce choices, in addition to sequencing and parallelization of services. On the other hand, we are able to explain why services participate in the CRP as they deal with risk or damage, or as they are dependent from other participating services.

Moreover, we also exploit the powerful conventional process mining algorithms, notably Alpha [START_REF] Van Der Aalst | Workflow Mining: Discovering process models from event logs[END_REF], which extract process schema from file log repositories. Such algorithms are also helpful in crisis management for deducing CRP. However our CRP deduction approach differs from the process mining one as the activities (services) and their coordination are no more extracted from log files but rather deduced by a matching process and by using pre-exiting relations between the services modeled in a crisis meta-model. Therefore, we start from scratch i.e. without execution cases (log files) and we guarantee to express all the possible scenarios. Finally, unlike [START_REF] Schonenberg | Supporting flexible process through recommendations based on history[END_REF][START_REF] Ribeiro | A recommender system for process discovery[END_REF][START_REF] Rangiha | A framework to capture and reuse process knowledge in business process design and execution using social tagging[END_REF], we focus on process deduction to deal with a situation as it is the essential requirement of crisis cells.

The crisis meta model

This section presents a meta-model for defining both crisis description (facts) and the required knowledge for crisis treatment (business knowledge), which are two fundamental dimensions that must be taken into account in crisis representation [START_REF] Bénaben | A metamodel for knowledge management in crisis management[END_REF]. This resulting meta-model is given in Fig. 2 as an UML class diagram. While facts are case dependent, knowledge are specific to a crisis domain (flood, forest fire, earthquake…). They enable a declarative description of the resolution process since its control structure, i.e. the coordination of the actions, is not explicitly described but derived/discovered from the facts and binary relations between actions. This aspect differentiates our proposition from other meta-models ( [START_REF] Bénaben | A meta-model and its ontology to guide crisis characterization and its collaborative management[END_REF]) where the process is explicitly described and requires a heavy work from the users. Also our metamodel includes both the concept of plan and services. A risk/damage could be associated to a set of individual services and/or to a resolution plan encapsulating a set of coherent services. This offers to the user flexibility to describe and organize in a modular way his knowledge. Fig. 2. Crisis meta-model

Facts representation

The facts are abstracted in the Risk/Damage class. This concept corresponds to an observed fact in the field, which can either be risk or damage. Damage is a negative situation affecting for instance population (e.g., flooded house with people inside), building (e.g., flooded school), road (e.g., cut-off road)…, while risk is the potential for damage. For each risk or damage, we store its type (risk or damage), and its nature, i.e. if it is already known or not. When it is known, it is linked to the knowledge base and more particularly to the corresponding Intrinsic Risk/Damage (relationship correspond), which gives access to the known solution for treating it. When it is unknown, the crisis cell has to specify at run time how to deal with this new risk or damage, indicating which services to be deployed in the field (relationship deal with). In addition, for each risk or damage, we store a specific property indicating if the risk or damage has priority or not. A priority risk or damage has to be considered in the deduction process when deducing the CRP, while a not priority risk or damage will not be taken into account by the current deduction; it will be taken into account later, when another deduction is made. Crisis cell members may change the value of this property according to the urgency of risk or damage.

Business knowledge representation

Business knowledge completes the facts with the following concepts: Intrinsic Risk/Damage, Plan, Service, Actor, Data, Choice, Condition and Type. Operational actions that can be executed in the field are modelled as services provided by crisis actors. For each service, we store its input and output data. Moreover, the rules relative to the use of these services are expressed as relations between services (relationship link). A relation between two services may be require, cause, or follow. Types require and cause define a strong relation among considered services, indicating that both services have to be executed one after the other: require indicates there is a precedence relation among them while cause indicates that there is a succession relation among them. At the opposite, the type follow defines a weak relation among considered services, indicating that a service will be performed after another, but not necessarily right after. In addition, we also have introduced another relation between services, namely the choice relation. The idea is to support alternative modelling, each alternative being a possible solution to deal with an issue. A condition defines when using this alternative. Finally services provided by actors may be used to deal with intrinsic risk, possibly as part of a plan, which corresponds to already specified set of actions to be undertaken to address an issue [START_REF] Tahir | A collaborative information system architecture for process-based crisis management[END_REF].

Deduction of crisis resolution processes

A Crisis Resolution Process (CRP) is a process driving crisis resolution. It includes the set of ordered actions (services) to be undertaken by crisis actors and their coordination.

Deduction principle

Basis of deduction are both facts observed in the field, which correspond to risks to be reduced or damages to be repaired, and business knowledge of actors, specified as services, and relations between these services. Deduction principle is given in Fig. 3 as a BPMN process diagram. This process includes three main steps.

Fig. 3. Deduction Principle process

The first step is the Service Matching step, which matches observed facts with business knowledge and more precisely, intrinsic risk and damage and corresponding services. The result of this step is a minimal set of services to be deployed. This set is then completed in the Service Expansion step. To do this, we exploit the relations between services to identify additional services to be deployed. The result of this expansion step is the set of services to be coordinated in the corresponding CRP. Finally, the Service Ordering step is responsible for ordering services w.r.t. their relation. It is visualized as a sub-process in Fig. 3. First, we build a matrix describing dependences existing between considered services from relations existing between them. As in conventional process mining algorithms [START_REF] Van Der Aalst | Process Mining -Data Science in Action[END_REF], we consider three types of dependences:

• causal dependence: a causal dependence between services a and b, denoted a b, indicates that service a has to be executed just before service b,

• parallel dependence: a parallel dependence between services a and b, denoted a || b, indicates that services a and b are executed in any order, • unrelated dependence: an unrelated dependence between services a and b, denoted a # b, indicates that services a and b are completely independent one from another, that is it does not exist any causal or parallel dependence between them. Then, from this dependence matrix, we build the corresponding Petri Net from which we derive the corresponding BPMN-like diagram. The Petri net serves as a support for CRP analysis, simulation and Algo. 1. algorithm ServiceExpansion validation, while the BPMN serves as a support for CRP execution. Note that the Petri net formalism has been chosen as it provides formal and executable specifications to analyze, simulate, check and validate the described process [START_REF] Van Der Aalst | The application of Petri Nets to Workflow Management[END_REF] while BPMN has been chosen as it is the language of the process engine that we use in this work.

Algorithms for CRP deduction

Different algorithms have been written to implement CRP deduction. We present below the two main ones.

The first algorithm is the algorithm ServiceExpansion, detailed below in Algo. 1. This algorithm implements service expansion using choice, cause and require relations as follows. The idea is to add services that are required to, consequence of, or alternative to each service obtained after matching. For that we use two sets of services, namely tobeExpanded, whose initial value is the set of services obtained after matching, and Expended, and we add to Expended both a service x from tobeExpanded and services connected to x by require, cause or choice relation.

The second algorithm is the algorithm ServiceCoordination, detailed in Algo. 2. This algorithm deduces the CRP to be deployed in the field using a process mining-based approach. More precisely this algorithm extends the Alpha algorithm [START_REF] Van Der Aalst | Workflow Mining: Discovering process models from event logs[END_REF], which is the key algorithm in process mining, and which needs to be revisited to address the deduction of CRP service coordination. Unlike Alpha [START_REF] Van Der Aalst | Workflow Mining: Discovering process models from event logs[END_REF] in which the services (activities) and their coordination are mined from log file repositories, we deduce services in line with facts observed in the field (cf. Algo 1) and define their coordination by using pre-exiting relations between deduced services (cf. Algo 2). Therefore, we start from scratch, i.e. without examples of execution (log files) and we guarantee to express all the possible scenarios. While in Alpha a dependency between two activities is deduced from their direct succession in traces recorded in log files, we generate it from preexisting relations between services. We also add artificial services to represent complex patterns such processes starting with parallel activities. Moreover, as in Alpha, the Petri net formalism supports process description in terms of places, transitions, corresponding to actions to be executed, and arcs, connecting places and transitions. However, as defended before, our algorithm also provides a BPMNlike representation of the CRP, which is more convenient to crisis cell members.

More precisely, our algorithm first builds the dependence matrix basis of the CRP deduction, but in a very different way from Alpha as we exploit business knowledge and not execution logs. To get into detail of this matrix building, for the set of services obtained after service expansion, we build causal dependences in the matrix from require, cause and follow relations (line #1 in algo 2). We also analyze relation between these services to eventually define new services which correspond to choices and specify unrelated dependences according to choice relations (line #2). Finally, parallel dependences between these services are deduced using the following rules (line # Then, the construction of the Petri net is fairly similar to Alpha. However, even in this similar part, we extend Alpha adding specific places and transitions Algo. 2. algorithm ServiceCoordination (adaptation of Alpha algorithm) to build processes possibly starting with parallelism or alternative. More precisely, as Alpha, we identify initial and final services, which are services to be executed respectively at the beginning and at the end of the CRP. Then, the novelty is to define two virtual transitions: Start and End. Start is connected to each initial service so that they could be performed after Start. Also, each final service is connected to the End transition, so that the End transition merges the results of the final services. Another important difference with Alpha is that we are able to deduce alternatives involving empty activities as we automatically add these activities when building the matrix. Thus we overcome some limitations of Alpha (e.g., [START_REF] De Medeiros | Process mining: extending the -algorithm to mine short loops[END_REF][START_REF] Wen | Mining process models with non-free-choice constructs[END_REF]). Finally, the part of the algorithm inspired by Alpha is (i) the determination of X, the minimum set of couples (Servicesa, Servicesb) for which, each sa in Servicesa has a causal dependence with each sb in Servicesb as well as sa and sb are unrelated (line #13), (ii) the determination of Y, which is a subset of X (line #14) and (iii) the aggregation of the final Petri net (line #18). All other algorithm lines are specific to our CRP deduction.

The resulting Petri net is then mapped into a BPMN-like diagram (line #19), which does not include BPMN pool and lanes to be more readable for crisis cell members. We do not detail this mapping as it is quite classic (e.g., plug-in PROM supports mapping to BPMN from Petri net [START_REF]ProM: Process mining toolkit[END_REF]), but we highlight its specificities in GéNéPi. Indeed, in GéNéPi, BPMN is not only a notation for CRP visualization but also the executable process language of Iterop, the process engine that supports CRP execution. Thus to obtain a fully executable specification, we have mapped flowing conditions, i.e. conditions attached to sequence flow flowing from open exclusive gateways to activities (i.e., services) in the BPMN-like diagram. More precisely, if use conditions of services are defined in the meta-model, then these use conditions are the flowing conditions. Otherwise, the algorithm automatically adds an out data to the activity preceding an open exclusive gateway, and defines for each sequence flow flowing from this open exclusive gateway a condition in which this out data is involved. Another interesting aspect in this mapping is the labelling of services with the facts they deal with. Thereby the algorithm labels each service with the facts justifying the selection of the service in the CRP, making it possible to determine whether or not all activities related to a fact are carried out or not. Thus it is possible to modify crisis situation deleting facts processed from the list of facts to be taken into account. Finally, we simplify the CRP in removing Start and End services, which were introduced for consistency reasons when building the Petri net, but which are no more useful in the BPMN. We also remove added services in the Petri net for syntactic reasons but useless in the BPMN.

Approach experiment

We have conducted an experiment considering the flood of the Loire in June 2016 as a case study. This experiment has been set up in collaboration with the crisis cell of Orléans, prefecture of department 45 in France, in the context of the GéNéPi project 1 , which aims at making civilian crisis (notably floods) management easier for crisis cells. Orléans being often deeply affected by Loire's floods, their mastering is of utmost importance. Members of the crisis cell were the Préfet, head of prefecture, the COD, which is the operational committee set up within the crisis cell and finally the representatives of the different actors acting in the field (e.g., DDT that are responsible for dykes supervision…). The experiment has focused on the simulation of several days of the last flood of the Loire in June 2016. This section introduces the considered case study and reports on the crisis meta-model instantiation and the CRP deduction.

Case study

Flood of the Loire in June 2016 lasted 12 days. In our experiment, in accordance with the crisis cell, we only focused on days 1, 3, 7 and 9, thus highlighting different interesting aspects of the response (i.e., the deduced CRP) recommended by our deduction algorithm. This paper only reports on a simplified version of day 7, is complex enough to illustrate the value added by our algorithm to deduce complex processes that include crisis cell decision making and hierarchical communication towards different ministries to which crisis cell is accountable.

On day 7, rainfall forecast is substantial. Loire level should rise significantly and major concerns of crisis cell are the following facts:

• risk of civilian casualties in nursing home Saint Pryvé Lake: the nursing home has to be evacuated, • risk of flooding of motorway A71: the motorway has to be partly cut off, • risk of dyke failure in Saint Pryvé Saint Mesmin: municipality of Saint Pryvé Saint Mesmin, next to Orléans, could be flooded and some districts of the municipality could be evacuated.

Day 7: facts and knowledge modelling

The three previous risks are modelled as Risk/Damage in the recommended meta-model. Their priority is 1, they are known and are linked to three Intrinsic Risk/Damage via the relationship correspond.

1 http://gind.mines-albi.fr/en/projet/genepi These intrinsic risks serve as a basis for defining the required knowledge to address dyke failures, flooding of nursing homes and flooding of motorways. We do not detail below the modelling of risks but we rather focus on the modelling of services (without detailing their in and out data) along with their relation, as they serve as a basis for matrix dependence building.

On Day 7, 17 services, offered by 6 actors, have to be selected to participate into the CRP. These services and the corresponding actors are given in Table 1. In addition, Table 2 shows existing relations between these services and Table 3 gives their use conditions. On day 7, we only have two use conditions respectively for services Decision-making for evacuation and Dyke supervision continuation. Both conditions involve two out data from service Dyke state evaluation, namely dykeState, whose value is low or solid, and dykeSupervision whose value is yes or no.

Day 7: deduced CRP

As explained before we deduce a CRP as a Petri net diagram and we map it into a BPMN one. Due to lack of space, we only provide the BPMN-like diagram in Fig. 4.

Both diagrams are built from the dependence matrix given in Table 4. This dependence matrix is built from knowledge stored in the meta-model (cf. Tables 1, 2 and3). More precisely, dependence matrix building exploits relations between services (require, cause, follow and choice). In this building, we also identify services that are useful for modelling choices. For instance, services 18 and 19 have been added automatically as we have use conditions for services 5 and 7 and no choice relation for each of them. As a consequence, the algorithm identifies choice relations between services 5 and 18 and between services 7 and 19 and produces the corresponding unrelated dependences. In addition, unrelated dependence is symmetric: both a#b and b#a are indicated in the matrix (e.g., 5#18 and 18#5). Moreover, causal dependences are indicated regardless the way of reading the matrix: in line or in column (e.g., 1 2 and 2 1).

From this dependence matrix we deduce the corresponding Petri net diagram. All services are explicitly represented as transitions, including the added services 18 and 19 even if they do not correspond to services offered by crisis actors. Note that the Petri net focuses on the behavioral dimension of crisis actors: only services to be executed and their synchronization is deduced. It is then possible for crisis cells to simulate, validate and analyze the behavior of the CRP. Then the corresponding BPMN diagram is mapped from the Petri net. In the mapping process we label outgoing sequence flows from exclusive gateways with conditions that have to be checked. For instance, as illustrated in Fig. 5, the outgoing sequence flow from alternative gateway to service Decision-making for evacuation of actor CrisisCell is labelled with the following condition defining when the sequence flow is executed: dykeState='low'. Moreover, in the BPMN diagram, some services are labelled with risks while others are not. Regarding labelled services, this means that the risk is effectively taken into account when all the corresponding services have been performed. For instance, the risk of dyke failure in Saint-Privé Saint-Mesmin will be taken into account, and thus will be removed from Risk/Damage in the meta-model, when the services Prepare for dyke supervision, Dykesupervision and Dyke state evaluation will be completed. Regarding non labelled services, theycorrespond to decision-making services (e.g., Decision-making for evacuation), which rather generate new risks, or to services implementing hierarchical communication towards different ministries to which crisis cells are accountable (e.g., Send decision report). Finally, as explained in Sec. 4, we also simplify the CRP removing services Start and End along with syntactic added services (e.g., services 18 and 19). The mapping result is given in Fig. 4.

Discussion and conclusion

This paper has addressed deduction of processes driving crisis resolution from relevant facts observed in the field and business knowledge of actors involved in crisis resolution. Deduction step is a key step in crisis management (and in GéNéPi) as it provides crisis cells with guidelines for crisis resolution. Moreover these guidelines are consistent with both facts observed and crisis actors' knowledge. The recommended solution includes (i) a meta-model supporting facts and knowledge modelling and (ii) a set of algorithms implementing crisis resolution process deduction. In addition, the paper has reported on an experiment conducted in collaboration with a crisis cell, considering the flood of the Loire in June 2016 as a case study. The experiment took place in the context of the GéNéPi project, funded by the French research national agency. The experiment enabled us to assess difficulties of crisis cell members in modelling knowledge as instance of our meta-model. Fact modelling was really easy while knowledge modelling was more touchy, mainly the modelling of services and relations between them. Regarding services, the issue was their identification. Indeed, the Préfet, head of crisis cell, rather needs a macroscopic vision of the actions executed in the field, while the representatives of the different actors acting in the field need a more microscopic view, in which their actions are detailed. In the experiment, after discussion among the crisis cell members, we have modelled services according to Préfet's needs. Regarding relations between services, require and cause relation were easily identified by crisis cells members themselves, but it took time to identify choice relations. In addition, the experiment enabled us to assess the match between the deduced CRPs and what crisis cell members were expecting. Feedbacks were very good, as for each crisis day, the deduced CRP has matched with what the crisis cell would do, taking into account every risk or damage and not forgetting even one thing. Above all we have heard from crisis cell members that the inclusion of choices and decision-making services matched the way in which a crisis cell works, which is strength of our solution.

On the other side, the recommended deduction algorithm is a step forward with respect to existing solutions. Regarding existing deduction algorithms, unlike [START_REF] Bénaben | Semantic issues in model-driven management of information system interoperability[END_REF][START_REF] Macé Ramète | Towards a predictive model for decision support in road crisis management[END_REF], we recommend a declarative approach for knowledge modelling, thus making knowledge management easier to crisis cell members. Moreover we fully exploit this knowledge: (i) we deduce choices along with corresponding conditions, in addition to sequencing and parallelization of services and (ii) we label services with corresponding risk or damage. Labelling services is very important as it explains service selection in the CRP (thus the Préfet can defend undertaken actions to hierarchical authorities to which he is accountable) and it makes possible the deletion of risk or damage in the meta model (as they have been addressed in the field). Regarding process mining algorithms, our recommended algorithm extends the Alpha algorithm [START_REF] Van Der Aalst | Workflow Mining: Discovering process models from event logs[END_REF] building processes possibly starting with parallelism or alternative. Moreover, our algorithm does not need as input any log file as it only exploits knowledge for crisis actors for CRP deduction from facts observed.

However, three main improvements are required in our work. The first one is related to the consistency of modelled knowledge, and more precisely the consistency of relations between services. We did not investigate this point and have planned to do it shortly.

The second one is related to the filtering of services before Petri net deduction. So far, this filtering boils down to a unique selection of each service, even if it is required two times in the CRP. We will also address this key point in the next future. Finally, we are also interested in mining the social dimension that exists between actors involved in crisis cell resolution to determine in what this dimension influences crisis coordination.

  3): If a b and a c and not (b # c) Then b || c If a || b and a c and not(b c) then b || c

Fig. 4 .

 4 Fig. 4. Deduced CRP as a BPMN-like diagram

Table 1 .

 1 Services required on day 7

Table 2 .

 2 Relations between services

Table 3 .

 3 Use conditions for services

Table 4 .

 4 Dependence matrix on day 7


Sergey Vilov

Bastien Arnal

Eliel Hojman

Yonina C Eldar

Ori Katz

Emmanuel Bossy
email: emmanuel.bossy@univ-grenoble-alpes.fr

Super-resolution photoacoustic and ultrasound imaging with sparse arrays

It has previously been demonstrated that model-based reconstruction methods relying on a priori knowledge of the imaging point spread function (PSF) coupled to sparsity priors on the object to image can provide super-resolution in photoacoustic (PA) or in ultrasound (US) imaging. Here, we experimentally show that such reconstruction also leads to super-resolution in both PA and US imaging with arrays having much less elements than used conventionally (sparse arrays). As a proof of concept, we obtained super-resolution PA and US cross-sectional images of microfluidic channels with only 8 elements of a 128-elements linear array using a reconstruction approach based on a linear propagation forward model and assuming sparsity of the imaged structure. Although the microchannels appear indistinguishable in the conventional delay-and-sum images obtained with all the 128 transducer elements, the applied sparsity-constrained model-based reconstruction provides super-resolution with down to only 8 elements. We also report simulation results showing that the minimal number of transducer elements required to obtain a correct reconstruction is fundamentally limited by the signal-to-noise ratio. The proposed method can be straigthforwardly applied to any transducer geometry, including 2D sparse arrays for 3D super-resolution PA and US imaging.

Introduction

Ultrasound [START_REF] Wells | Ultrasound imaging[END_REF] (US) and photoacoustic [START_REF] Beard | Biomedical photoacoustic imaging[END_REF] (PA) imaging are now widely applied biomedical imaging modalities. They both usually use multielement transducer arrays as ultrasonic detectors for acquiring acoustic signals. Developed for two-dimensional (2D) or cross-sectional imaging, linear transducer arrays are widely spread in research and clinical applications. For threedimensional (3D) single shot imaging, 2D array matrices should be used instead of linear arrays. However, the availability of 3D imaging equipment is first limited by the sophisticated fabrication process involved: probes for 3D imaging may have several thousands of elements to connect which makes the assembly of such probes technically difficult [START_REF] Fenster | Three-dimensional ultrasound imaging[END_REF] . Second, to control simultaneously as many elements as possible, sophisticated ultrasound electronics with a very large number of channels are needed. Decades ago, sparse arrays were proposed [START_REF] Turnbull | Beam steering with pulsed two-dimensional transducer arrays[END_REF][START_REF] Smith | High-speed ultrasound volumetric imaging system. i. transducer design and beam steering[END_REF] to reduce substantially the number of transducer elements required for 3D US imaging. The corresponding reduction scheme suggests using a selection of elements of a dense periodic array without changing the total transducer aperture. Such a selection can either be a random subset of elements [START_REF] Turnbull | Beam steering with pulsed two-dimensional transducer arrays[END_REF] or a defined pattern [START_REF] Smith | High-speed ultrasound volumetric imaging system. i. transducer design and beam steering[END_REF] . Experimental investigations of imaging performance of sparse arrays in US imaging showed that the sparse array can provide a diffraction-limited resolution similar [START_REF] Roux | 3d diverging waves with 2d sparse arrays: A feasibility study[END_REF][START_REF] Roux | Validation of optimal 2d sparse arrays in focused mode: Phantom experiments[END_REF] to that of the all-element array. Special reconstruction techniques have been recently proposed for sparse array US imaging. In particular, the convolutional beamforming algorithm was reported [START_REF] Cohen | Sparse convolutional beamforming for ultrasound imaging[END_REF] to provide better (although still diffraction-limited) resolution than standard delay-and-sum beamforming. It is worth noting that in radar technologies it has also been shown that advanced compressed sensing methods [START_REF] Eldar | Compressed sensing: theory and applications[END_REF] permit preserving the diffraction -limited resolution when using radars with a reduced number of elements [START_REF] Rossi | Spatial compressive sensing for mimo radar[END_REF][START_REF] Cohen | Summer: Sub-nyquist mimo radar[END_REF] .

Sparse arrays have also been widely applied in 3D photoacoustic tomography [START_REF] Ephrat | Three-dimensional photoacoustic imaging by sparse-array detection and iterative image reconstruction[END_REF][START_REF] Chaudhary | Comparison of reconstruction algorithms for sparse-array detection photoacoustic tomography[END_REF][START_REF] Roumeliotis | Singular value decomposition analysis of a photoacoustic imaging system and 3d imaging at 0.7 fps[END_REF][START_REF] Roumeliotis | 3d photoacoustic imaging using a staring-sparse array with 60 transducers[END_REF][START_REF] Wang | Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography[END_REF][START_REF] Dean-Ben | Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography[END_REF][START_REF] Han | Three-dimensional optoacoustic reconstruction using fast sparse representation[END_REF] . In most of these studies sparse arrays were used to obtain diffraction-limited images with model-based reconstruction [START_REF] Chaudhary | Comparison of reconstruction algorithms for sparse-array detection photoacoustic tomography[END_REF][START_REF] Roumeliotis | Singular value decomposition analysis of a photoacoustic imaging system and 3d imaging at 0.7 fps[END_REF][START_REF] Roumeliotis | 3d photoacoustic imaging using a staring-sparse array with 60 transducers[END_REF][START_REF] Wang | Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography[END_REF][START_REF] Dean-Ben | Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography[END_REF][START_REF] Han | Three-dimensional optoacoustic reconstruction using fast sparse representation[END_REF] . The model-based reconstruction approach considers a forward linear model expressed as S = A A AT , where S are the acquired signals, T is the object to reconstruct, and the propagation matrix A A A is a library containing the point spread function (PSF) at all points of the discretized imaging zone. The reconstruction consists in finding the object that minimizes a cost function defined as the sum of a fidelity term (taking into account the measurement data and the model) and a regularization term (taking into account measurement noise and prior knowledge on the object).

As the main novelty of this work, we investigate the possibility of using sparse ultrasound arrays in the context of superresolution imaging with the model-based approach. The classical resolution in both PA and US imaging of biological tissues at depth is limited by acoustic diffraction. Specifically, it is the absorption of ultrasound in tissues that limits the highest ultrasound frequency detectable at a given depth, and therefore the resolution via the diffraction phenomenon. In particular, the depth-to-resolution ratio turns out to be of the order of 100 for both US and PA imaging of biological tissues. As an example, for a 10 MHz US probe, that can be ordinarily used for US or PA imaging at a few centimeters depth in tissues, diffraction limits the resolution to 100-200 µm in typical imaging conditions [START_REF] Foster | Advances in ultrasound biomicroscopy[END_REF] .

Imaging beyond the diffraction limit has first been investigated in optical microscopy, resulting in the advent of Nobel-prize winning pioneering methods such as photo-activated localization microscopy (PALM) [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF] and stimulated emission depletion (STED) microscopy [START_REF] Hell | Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy[END_REF] . Following the aforementioned advances in optics, several super-resolution techniques have been proposed more recently for both US and PA imaging. Super-resolution methods include localization-based imaging, fluctuationbased reconstruction and model-based reconstruction. As for optical imaging, the localization approach in US [START_REF] Viessmann | Acoustic super-resolution with ultrasound and microbubbles[END_REF][START_REF] O'reilly | A super-resolution ultrasound method for brain vascular mapping[END_REF][START_REF] Desailly | Sono-activated ultrasound localization microscopy[END_REF][START_REF] Harput | 3-d super-resolution ultrasound imaging with a 2-d sparse array[END_REF] and PA [START_REF] Vilov | Overcoming the acoustic diffraction limit in photoacoustic imaging by the localization of flowing absorbers[END_REF][START_REF] Dean-Ben | Localization optoacoustic tomography[END_REF] imaging relies on detection of individual scatterers or absorbers. The idea is that the coordinates of a point-like source can be determined with a precision much better than the size of the imaging PSF provided that this PSF can be separated from those of the other sources in some parameter space. This separation condition imposes a low concentration of sources, and therefore the use of contrast agents for visualization of biological structures. Fluctuation-based PA [START_REF] Chaigne | Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination[END_REF][START_REF] Chaigne | Super-resolution photoacoustic imaging via flow-induced absorption fluctuations[END_REF] and US [START_REF] Bar-Zion | Fast vascular ultrasound imaging with enhanced spatial resolution and background rejection[END_REF] techniques, based on the principles of super-resolution optical fluctuation imaging (SOFI) exploit uncorrelated fluctuations from different sources. While fluctuation-based approaches eliminate the need for isolating individual sources, they remain limited in terms of spatial resolution improvement and temporal resolution. It was also shown that model-based reconstruction with sparsity constraints on the sample, an approach originating from the field of compressed sensing [START_REF] Eldar | Compressed sensing: theory and applications[END_REF] , could yield super-resolved US [START_REF] Tuysuzoglu | Sparsity driven ultrasound imaging[END_REF][START_REF] Zhao | Single image super-resolution of medical ultrasound images using a fast algorithm[END_REF] and PA 33 images. A major advantage of model-based reconstruction over the previously mentioned localization-based and fluctuation-based techniques is that by requiring (in principle) only a single-shot acquisition it permits a high temporal resolution. Some very recent works [START_REF] Hojman | Photoacoustic imaging beyond the acoustic diffraction-limit with dynamic speckle illumination and sparse joint support recovery[END_REF][START_REF] Murray | Super-resolution photoacoustic microscopy using blind structured illumination[END_REF][START_REF] Bar-Zion | Sparsity-based ultrasound superresolution hemodynamic imaging[END_REF] have proposed to mix fluctuation-based and model-based techniques to achieve super-resolution in US [START_REF] Bar-Zion | Sparsity-based ultrasound superresolution hemodynamic imaging[END_REF] and PA imaging [START_REF] Murray | Super-resolution photoacoustic microscopy using blind structured illumination[END_REF] . In PA imaging, it was also demonstrated 34 that super-resolution can be obtained via the joint support recovery through the model-based approach. In particular, the vector S is composed of data measured at several PA acquisitions, with each acquisition corresponding to a random speckle illumination of the same absorbing structure forming the joint support.

In US imaging, the use of 2D sparse arrays for super-resolution imaging has recently been demonstrated by use of a localisation approach and microbbubles [START_REF] Harput | 3-d super-resolution ultrasound imaging with a 2-d sparse array[END_REF] . Here, we investigate the possibility to perform both US and PA super-resolution imaging of sparse samples with a label-free model-based reconstruction approach using a sparse array. More specifically, we apply sparsity-constrained model-based reconstruction to perform 2D super-resolution imaging of sparse test samples with only 8 (out of 128) elements of a linear transducer array. We show that this reconstruction approach can be used for both PA and plane-wave US imaging with the same experimental setup. To build the matrix A A A that describes the forward model, we propose a novel method that involves only one PSF acquisition, as opposed to the measurement of the full set of PSFs in the imaging zone. We also report simulation results showing how the reconstruction quality is related to the number of transducer elements and the signal-to-noise ratio (SNR).

Results

To demonstrate that sparsity-constrained model-based reconstruction can provide super-resolution in sparse-array PA or US imaging, we carried out two proof-of-concept experiments (one for US and one for PA imaging). The goal of each experiment was to recover a super-resolved cross-sectional image of a sparse five-channel microfluidic sample by processing data received by different subsets of elements of a linear ultrasound probe. The imaging configuration for both experiments is shown in Fig. 1 and a picture of the setup is available in Supplementary Materials. Importantly, in both PA and US experiments, we used the same imaging equipment, the same acquisition geometry, samples of identical structure, and the same reconstruction method. Only the nature of the contrast and the way the ultrasound wave is generated were different in these experiments. In the PA experiment, the microfluidic channels were filled with absorbing liquid and illuminated by pulsed light. In the US experiment, the microfluidic channels were filled with air and a plane ultrasound wave was emitted by the transducer array. In both experiments, the measurement data consisted of the ultrasound signals detected by the array for a single shot excitation (light pulse or plane wave). Prior to imaging the five-channel samples, the PSF required to build the forward model was measured for each type of sample (absorber-filled or air-filled) in the corresponding imaging mode (PA or US).

The image reconstruction consisted in finding the object that minimizes a cost function. This cost function was based on the forward model (derived from the PSF calibration) and a L1-norm regularization term used to suppress the measurement noise and select a sparse object. This model-based reconstruction was then compared to the conventional delay-and-sum approach taken as a typical diffraction-limited reconstruction approach. Further details on the experimental and reconstruction methods can be found in the Methods section.

We first demonstrated for both PA and US imaging that sparsity-based reconstruction led to super-resolution images of sparse samples when using all the 128 available elements of the probe. The conventional delay-and-sum images are shown in Fig. 2 b and g. It can be noticed that the conventional reconstruction is affected by diffraction-limited resolution: two neighboring Then, as a major result of this work, we demonstrated that it could also be possible to obtain similar super-resolved images using a sparse array in the place of the initial 128-element dense array. To emulate sparse-array imaging, we applied the same reconstruction approach to data acquired by only a fraction of the probe's elements. To maintain the same conventional diffraction limit, the resulting probe aperture was kept constant by including the first and the last elements of the probe. The other elements were regularly distributed along the probe. The objects reconstructed with 16, 8 and 4 elements are compared to the full-probe (128 elements) reconstruction in Fig. 3. The images in Fig. 3 show that in both PA and US experiments the five channels could have been imaged using down to only 8 array elements. With less than N min = 8 elements, we were not able to reconstruct the imaged structure properly. Random distributions of the transducer elements were also considered, but no significant differences in the results were found as compared to the linearly distributed case (See Supplementary Figures 2 and3).

In order to shed some light on the parameters that condition the minimal number N min of transducer elements needed to yield a faithfull reconstruction, we carried out a series of numerical simulations. In these simulations, as a measure of the reconstruction quality, we studied the correlation C between the reconstructed object and a modelled ideal object, as a function of the number of transducer elements and the signal-to-noise ratio (SNR). In this work, we define the SNR as the ratio between the peak amplitude of the radio frequency (RF) signal and the standard deviation of the measurement noise computed over a signal-free region of the RF data. The simulation results are illustrated in Fig. 4. As could be intuitively expected, these results show that for any fixed number N of transducer elements involved in the reconstruction, the reconstruction image quality increases with the SNR. Moreover, the SNR that assures a given reconstruction quality approximately scales as the square root of the number of elements, apart from N = 2. This is illustrated for C = 0.8 in Fig. 4 with the dashed line following SNR∝ √ N. For N = 2, our simulations still predict the possibility to reconstruct the imaged object, but with a much stronger requirement on the SNR.

Discussion

Previous studies reported super-resolution images obtained via model-based reconstruction with sparsity priors on the imaged object, in the field of US imaging [START_REF] Tuysuzoglu | Sparsity driven ultrasound imaging[END_REF][START_REF] Zhao | Single image super-resolution of medical ultrasound images using a fast algorithm[END_REF] and more recently in PA imaging [START_REF] Egolf | Sparsity-based reconstruction for super-resolved limited-view photoacoustic computed tomography deep in a scattering medium[END_REF] . In this work, we first illustrated the generality of the sparsity-constrained model-based approach by applying this method to obtain experimental PA and US images of sparse test objects, with the same imaging equipment being used in the PA and US experiments. As a key result of the present work, we demonstrated that imaging sparse objects beyond the acoustic diffraction limit remains possible even when a very small number of transducer elements are used for the reconstruction, much below the number required by Nyquist spatial sampling. This suggests that in imaging with conventional number of transducer elements, a significant part of the measurement information might be redundant. By means of numerical simulations, we showed that if the SNR is very high then by reconstructing images with a perfectly known forward model it is possible to obtain the correct object reconstruction with only two transducer elements. According to the simulation results, more than two elements might be necessary to provide the correct reconstruction in the case of typical experimental values of the SNR: for a given SNR on each transducer element, using more elements may be considered equivalent to reducing the influence of the measurement noise. In practice, additional uncertainties of the forward model might further limit the reconstruction quality: in our experiments the minimal number of elements N = 8 required for the reconstruction with correlation C = 0.8 was above the number of elements N = 3 predicted by the simulations for the same SNR and the same correlation.

In addition, a certain difference in image reconstruction quality between PA images (left columns in Fig. 2 and3) and US images (right columns in Fig. 2 and3) can be noticed. In particular, the SNR in the US reconstruction images is lower than 4/10 in the PA images. In addition, the average interchannel distance in the US images (112 µm) slightly deviates from the true value (125 µm), whereas in the PA images the interchannel distance is restored correctly. This difference between the PA and US reconstruction images may result from several reasons. First, the SNR in the US experiment was lower than in the PA experiment. Second, the uncertainties of the forward model might have played a greater role in the US experiment.

As a novelty of our work, we derived the forward model by measuring only one PSF, rather than all PSFs in the field of view, which was done, for instance, in [START_REF] Egolf | Sparsity-based reconstruction for super-resolved limited-view photoacoustic computed tomography deep in a scattering medium[END_REF] . This approach assumes that PSFs differ from each other only by time delays that are calculated based on wave propagation in water. However, in our experiments this condition was not fully satisfied at least due to the presence of a PDMS layer between the sources and the receivers. As the US experiment involves a double crossing of this PDMS layer, it is likely that in the US experiment the forward model is less accurate than in the PA experiment, leading to a less faithful reconstruction image.

In conclusion, we demonstrated experimentally the possibility of performing both PA and US super-resolution imaging of sparse samples with a sparse array (with down to 8 transducer elements) by applying the model-based reconstruction approach. Although the results demonstrated here were obtained in 2D imaging with a 1D linear array, a major advantage of the proposed approach is that it can be applied to any transducer geometry, and therefore in 3D imaging. Provided that the relative positions of the transducer elements are known, the proposed method of constructing the forward model from a single PSF measurement remains valid. It was shown by reconstructing simulated data that the SNR yields a fundamental limit on the number of transducer elements needed to provide a faithful reconstruction. The relative influence of the SNR and the uncertainties of the forward model on the reconstruction quality as well as the ultimate resolution limit are to be further investigated in order to fully understand and probably predict the number of transducer elements required for super-resolution imaging. While we demonstrated that it is possible to obtain super-resolution images of sparse test samples with sparse array, the presented method should be further investigated with more realistic experimental data before firm conclusions can be drawn on its performance in the context of biomedical PA and US imaging.

Samples

The microfluidic samples were prepared with a standard soft-lithography manufacturing technology 37 using polydimethylsiloxane (PDMS). Each sample consisted of five hollow channels sandwiched between two layers of PDMS. The thickness of the upper layer (the layer placed closer to the US probe in the experiments) was around 180 µm. This thickness was chosen as 5/10 small as possible since the influence of the PDMS on wave propagation was neglected in the model. Each channel was 40 µm wide (x-direction) and 50 µm deep (z-direction). The centre-to-centre distance between neighbouring channels was L cc = 125 µm. This distance was deliberately chosen smaller than the estimated lateral full width of half maximum (FWHM) of the PSF (155 µm). To acquire the PSF and build the forward model, additional samples containing only one microfluidic channel were prepared, the only channel being 10 µm wide and 50 µm deep. For PA imaging, the channels were filled with an absorbing dye solution (Patent Blue V, absorption peak at 640nm 38 ) to provide photoacoustic contrast, while for US imaging the channels were filled with air to provide acoustic contrast.

Experimental protocol

The imaging configuration used in experiments is illustrated in Fig. 1. The five microfluidic channels were aligned to cross the xz imaging plane perpendicularly, and positioned at the distance z f = 15 mm from the linear ultrasound probe, the distance z f corresponding to the elevational focus of the probe. As the US probe we used a capacitive micromachined ultrasonic (CMUT) array (L22-8v, Verasonics, USA: N = 256 elements with 128 consecutive elements used in our experiments, pitch ≈ 100 µm, center frequency f c ≈ 15 MHz). The probe and sample were immersed in a water tank. In the PA experiment, signals were generated by illuminating the sample with a 5 ns laser pulse (λ = 532 nm, fluence ≈ 3 mJ/cm 2 ) generated by a frequency-doubled Nd:YAG laser (Spitlight DPSS 250, Innolas Laser GmbH, Krailling, Germany). In the US experiment, a plane ultrasound wave (15 MHz center frequency short pulse) was emitted by the probe, and ultrasound waves backscattered from the sample were were recorded. To control emission and reception of ultrasound waves, the probe was connected to multichannel acquisition electronics (High Frequency Vantage 256, Verasonics, USA).

The Vantage system had two 128-channel connectors for US transducers, each connector providing interface for synchronous control of 128 transducer elements. As the transducer used in this work had 256 elements wired to a single probe connector, only 128 of its elements could be accessed simultaneously with our system. Although the whole transducer aperture could in principle be accessed using the internal Vantage multiplexer, this would require two distinct acquisitions. To demonstrate the possibility of single-shot super-resolution as a strong point of the proposed reconstruction approach, we did not use the Vantage system's multiplexer and restricted our acquisitions to 128 transducer elements.

To reconstruct objects in both PA and US experiments, the PSF was measured by applying the acquisition protocol described above for samples containing a single microfluidic channel. The only channel was placed in the center of the field of view of the ultrasound probe. The conventional lateral resolution defined as the lateral full width at half maximum (FWHM) of the PSF was measured after standard delay-and-sum reconstruction and turned out to be around 155 µm in both PA and US experiments.

6/10

Forward model

Under the assumption of a linear dependence between the quantity to image and the measured radio frequency (RF) signals, the voltage signals measured by the acquisition electronics can be expressed for both PA and US imaging as

s(t i , r k ) = A[T 0 (r)], (1)
where s(t i , r k ) is the signal at time t k measured by the transducer element located at r k , T 0 (r) is the quantity to image, and A is a linear operator that takes into account both the ultrasound generation/propagation, and the transducer response. In PA imaging, T 0 (r) is proportional to the absorption coefficient µ a (r), under the assumption of homogeneous light illumination and a homogeneous Gruneisen paramater. In US imaging, T 0 (r) is related to the distribution of the backscattering coefficient, provided that the single scattering regime is valid. By discretizing the object to reconstruct on a grid, this forward problem may then be written in a matrix form as

S = A A AT 0 , (2)
where S m×1 is a vector with all the RF data (m = number of time samples × number of transducer elements), T n×1 0 is the discretized version of the quantity to reconstruct (n = number of grid points in the reconstruction zone), and A A A m×n is a matrix representing the linear operator. Each of the n columns of A A A represents the RF response for one of the points 1..n of the reconstruction zone, i.e. each column of A A A is the vector data corresponding to the signals from a single point source in the imaging zone. The matrix A A A therefore contains the responses to each point source and can be considered as a matrix of all the system point spread functions (PSF) in the reconstruction space.

The matrix A A A can be modelled theoretically or measured experimentally. In some previous works, A A A was measured for each point of the reconstruction zone [START_REF] Roumeliotis | Singular value decomposition analysis of a photoacoustic imaging system and 3d imaging at 0.7 fps[END_REF][START_REF] Roumeliotis | 3d photoacoustic imaging using a staring-sparse array with 60 transducers[END_REF][START_REF] Egolf | Sparsity-based reconstruction for super-resolved limited-view photoacoustic computed tomography deep in a scattering medium[END_REF][START_REF] Murray | Super-resolution photoacoustic microscopy using blind structured illumination[END_REF] . In our work, the matrix A A A is also obtained experimentally but we perform a measurement for only one point source. As a point source, we use a single isolated microfluidic channel. All the columns of A A A are then derived from this measurement data assuming that although the signals from two distinct point sources differ from each other by their arrival time, apart from this they are identical. In other words, it is assumed that the response of each transducer element is the same for all points in the reconstruction zone. In our experiments, we also neglect the presence of the upper PDMS layer and therefore only consider wave propagation in water, i.e. in a homogeneous isotropic medium with a constant speed of sound c = 1500 m/s.

Under the assumptions stated above, when a single PSF is acquired for a source placed at {x j , z j }, any column i of the matrix A A A can be derived from the acquired RF data by shifting the signals for each transducer element k with the following time delay :

∆t i, j,k = δ US z j -z i c + 1 c (x k -x i ) 2 + (z k -z i ) 2 -(x k -x j ) 2 + (z k -z j ) 2 , (3)
where {x k , z k } are the coordinates of element k of the transducer array, {x i , z i } are the coordinates of point i of the reconstruction grid. In US imaging, δ US = 1 to account for the travel time of the emitted plane wave (δ US = 0 in PA imaging).

Image reconstruction

To solve the inverse problem, i.e. to find an estimate T of the solution T 0 to the forward model described above, the sparsityconstrained minimization approach similar to those already reported in previous works [START_REF] Tuysuzoglu | Sparsity driven ultrasound imaging[END_REF][START_REF] Egolf | Sparsity-based reconstruction for super-resolved limited-view photoacoustic computed tomography deep in a scattering medium[END_REF][START_REF] Murray | Super-resolution photoacoustic microscopy using blind structured illumination[END_REF] was used. This approach consists in finding T by solving the following minimization problem:

T = argmin T {||S -A A AT || 2 2 + α 2 ||T || 1 }. (4)
In Eq. ( 4) the l 2 -term ||S -A A AT || 2 2 = ∑ all data points p (S -A A AT ) 2 p corresponds to the least square fitting of the probed quantity T to the acquired data S. The regularization l 1 -term α 2 ||T || 1 = α 2 ∑ p |T p | is used to minimize the influence of the noise in S while selecting a sparse solution to underdetermined system (2). Regarding typical scales involved in the reconstruction, for reconstruction with 128 transducer elements the length of the vector T n×1 was n = 793 (number of reconstruction points in the field of view), the length of the vector S m×1 was m = 4608 (the number of transducer elements multiplied by the number of time samples for each element), the rank of the matrix A A A was rank(A) = 346, this rank being computed using the standard MATLAB function with the default value of the tolerance parameter. The rank of the matrix A A A being inferior to the number of unknowns n, system (2) was undetermined. The minimization operation described by Eq.( 4) was performed using a fast iterative shrinkage-thresholding algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Palomar | Convex optimization in signal processing and communications[END_REF][START_REF] Vu | Fista implementation in matlab[END_REF] . The object to image was reconstructed over a cartesian grid with a L cc /10 = 12.5 µm step. As in previous works [START_REF] Tuysuzoglu | Sparsity driven ultrasound imaging[END_REF][START_REF] Egolf | Sparsity-based reconstruction for super-resolved limited-view photoacoustic computed tomography deep in a scattering medium[END_REF][START_REF] Murray | Super-resolution photoacoustic microscopy using blind structured illumination[END_REF] , the value of the regularization parameter α was determined heuristically such that the reconstruction image was qualitatively close to the real sample. A 2D Gaussian filter (kernel size σ = 12.5 µm) and interpolation were applied to improve image visualization. The final image definition corresponds to a grid step of 3.125µm. Each reconstruction image was linearly normalized at the final stage via division of each pixel value by the maximum of the intensity in the initial reconstruction image.

The gaussian filter was applied to the reconstructed images only to smooth them out, as usually done for other approaches that may lead to point-like reconstructed structures (such as localization based super-resolution). The size of the filter has an influence on the apparent size of individual structures (the microchannels in our case). L1-minimization leads to point-like structure, as illustrated in Figures 2c and2h, which show that the size of the channels is not reconstructed correctly. However, the filter does not influence the resolution (defined as the ability to distinguish individual structures), as long as the filter size is smaller than the resolution limit of the method (here close to the distance between neighbouring channels). Here, the filter kernel size σ = 12.5 µm is chosen to smooth out the reconstructed object without affecting the ability to distinguish channels, but it should be kept in mind that the apparent size of the channels on the filtered image depends on this choice.

Simulations

Numerical simulations were performed to produce test data corresponding to imaging five point sources. The stated assumptions on our forward model were strictly followed. Specifically, we modelled the signals received on each transducer element as a time-delayed version of the same signal (simulated for five point sources), based on the delay law described by Eq. 3. The detected signals had a central frequency and bandwidth corresponding to those used in experiments and sampled at the same frequency as in experiments. Simulation data with different SNR values was generated by varying the amplitude of the modelled PA signals. In all simulations, Gaussian noise with a zero mean and the rms of σ n = 30 was added to the detected signals. Such noise corresponded to the noise produced by the acquisition electronics in our experiments. The simulated distribution of sources was then reconstructed following the same methods as used to reconstruct images from the experimental data. As a metrics of the reconstruction quality, we computed a normalized spatial cross-correlation between each reconstructed object T and the ideal reconstruction T true .

The ideal reconstruction T true was modeled by asserting the value of 1 to the 5 cells of the reconstruction grid corresponding to the positions of the 5 simulated sources. Then, the filtering and interpolation used for the reconstruction images were also applied to T true . In Fig. 4, the correlation C for each SNR was estimated by averaging over 100 noise realizations. For each reconstructed image T the correlation C was computed using Eq.( 5):

C = ∑ n i=1 T (i) • T true (i) ∑ n i=1 T 2 (i) • ∑ n i=1 T 2 true (i)
.

(5)

Figure 1 .

 1 Figure 1. Experimental setup. Samples carrying 5 parallel microchannels are placed perpendicularly to the imaging plane. Each channel is 40 µm wide and 50 µm deep, the center-to-center interchannel distance being L cc = 125 µm. In the PA experiment, the sample is illuminated by a laser pulse and the resulting PA signals are detected by a linear US probe. In the US experiment, a plane wave is emitted in the sample's direction and the backscattered signals are collected. The violet arrows designate the propagation direction of the ultrasound plane wave.

Figure 2 .

 2 Figure 2. Experimental images reconstructed in the PA (b-e) and US (g-j) experiments alongside with the corresponding PSFs (a,f). All the 128 elements of the probe were used. (b,g) Conventional beamforming cross-sectional images of the sample. The microchannels are indistinguishable as the center-to-center distance L cc = 125 µm is smaller than the lateral FWHM of the PSF (155 µm). (c,h) Images obtained with model-based reconstruction. The reconstruction recovers five distinct regions corresponding to the microfluidic channels. (d,i) Reconstruction images (c) and (h) after smoothing out by a 2D spatial Gaussian filter (σ = 12.5 µm) and interpolating on a 3.125 µm grid. (e,j) Yellow line: normalized amplitude profile on the PSF envelope images (a) and (f); red line: normalized amplitude profile on the envelope images (b) and (g); blue line: normalized amplitude profile on the filtered reconstruction images (d) and (i).

Figure 3 .

 3 Figure 3. Images obtained with model-based reconstruction in PA and US experiments using N = 4, 8, 16 and 128 transducer elements regularly distributed along the probe aperture. A 2D spatial Gaussian filter (σ = 12.5 µm) and interpolation on a 3.125 µm grid were applied after the reconstruction.

Figure 4 .

 4 Figure 4. Simulation results. (a) -Correlation C as a function of the SNR and the number of transducer elements N, dotted lines: SNR PA = 150 in the PA experiment, SNR US = 83 in the US experiment, dashed line: fitting function SNR∝ √ N for C = 0.8, (b-d) -typical simulation images: (b) -N = 64, SNR = 16, C = 0.85, (c) -N = 16, SNR = 10, C = 0.57, (d) -N = 128, SNR = 0.8, C = 0.23.

3/10

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grants 681514-ERC-COGCOHERENCE, 646804-ERC-COG-BNYQ, 677909-ERC-STG-DYNA-MIC), and by the Israel Science Foundation (grants 0100101 and 1361/18). O.K. was supported by an Azrieli Faculty Fellowship.

Competing interests

The authors declare no competing interests.